001: /*
002: * PoissonClusterer.java
003: *
004: * Created on April 19, 2006, 2:41 PM
005: */
006:
007: package org.geotools.graph.util.delaunay;
008:
009: import java.util.Collection;
010: import java.util.Iterator;
011: import java.util.Vector;
012: import org.geotools.feature.Feature;
013: import org.geotools.filter.Expression;
014: import org.geotools.graph.structure.Graph;
015: import org.geotools.graph.structure.basic.BasicGraph;
016:
017: /**
018: *
019: * @author jfc173
020: */
021: public class PoissonClusterer {
022:
023: private static double threshold = 1.0E-10;
024:
025: /** Creates a new instance of PoissonClusterer */
026: public PoissonClusterer() {
027: }
028:
029: public static Graph findClusters(Graph incoming, Expression base,
030: Expression target, double meanRate, int distance) {
031: Collection nodes = incoming.getNodes();
032: Iterator nodeIt = nodes.iterator();
033: Vector clusterNodes = new Vector();
034: Vector clusterEdges = new Vector();
035: System.out.println("x, y, actual, expected, probability");
036: while (nodeIt.hasNext()) {
037: DelaunayNode next = (DelaunayNode) nodeIt.next();
038: Feature nextFeature = next.getFeature();
039:
040: Object baseObj = base.getValue(nextFeature);
041: if (!(baseObj instanceof Number)) {
042: throw new RuntimeException("Expression " + base
043: + " must evaluate to a number on feature "
044: + nextFeature);
045: }
046: Object targetObj = target.getValue(nextFeature);
047: if (!(targetObj instanceof Number)) {
048: throw new RuntimeException("Expression " + target
049: + " must evaluate to a number on feature "
050: + nextFeature);
051: }
052: double totalBase = ((Number) baseObj).doubleValue();
053: double totalTarget = ((Number) targetObj).doubleValue();
054:
055: Collection newEdges = new Vector();
056: Vector newNodes = new Vector();
057: newNodes.add(next);
058:
059: if (distance == 1) {
060:
061: newEdges = next.getEdges();
062: // System.out.println("this node has " + newEdges.size() + " edges");
063: Iterator edgeIt = newEdges.iterator();
064: Vector removals = new Vector();
065: while (edgeIt.hasNext()) {
066: DelaunayEdge nextEdge = (DelaunayEdge) edgeIt
067: .next();
068: if (nextEdge.getEuclideanDistance() > 30) {
069: removals.add(nextEdge);
070: } else {
071: DelaunayNode neighbor = (DelaunayNode) nextEdge
072: .getOtherNode(next);
073: if (neighbor == null) {
074: throw new RuntimeException(
075: "We have a problem. "
076: + next
077: + " and "
078: + neighbor
079: + " should be neighbors via "
080: + nextEdge
081: + ", but aren't.");
082: }
083: Feature neighborFeature = neighbor.getFeature();
084: newNodes.add(neighbor);
085:
086: Object neighborsBaseObj = base
087: .getValue(nextFeature);
088: if (!(baseObj instanceof Number)) {
089: throw new RuntimeException(
090: "Expression "
091: + base
092: + " must evaluate to a number on feature "
093: + neighborFeature);
094: }
095: Object neighborsTargetObj = target
096: .getValue(nextFeature);
097: if (!(targetObj instanceof Number)) {
098: throw new RuntimeException(
099: "Expression "
100: + target
101: + " must evaluate to a number on feature "
102: + neighborFeature);
103: }
104: totalBase = totalBase
105: + ((Number) baseObj).doubleValue();
106: totalTarget = totalTarget
107: + ((Number) targetObj).doubleValue();
108: }
109: }
110: newEdges.removeAll(removals);
111: } else {
112: for (int i = 0; i <= distance; i++) {
113: Iterator nodeIt2 = newNodes.iterator();
114: Vector nodesToAdd = new Vector();
115: Vector edgesToAdd = new Vector();
116: while (nodeIt2.hasNext()) {
117: DelaunayNode next2 = (DelaunayNode) nodeIt2
118: .next();
119: // System.out.println("expanding from " + next2);
120: Collection edges = next2.getEdges();
121: // System.out.println("its edges are " + edges);
122: newEdges.addAll(edges);
123: Iterator another = edges.iterator();
124: while (another.hasNext()) {
125: DelaunayEdge nextEdge = (DelaunayEdge) another
126: .next();
127: DelaunayNode farNode = (DelaunayNode) nextEdge
128: .getOtherNode(next2);
129: // System.out.println("checking " + farNode + " in " + newNodes);
130: if (!(newNodes.contains(farNode))) {
131: nodesToAdd.add(farNode);
132: edgesToAdd.add(nextEdge);
133: // System.out.println("adding node " + farNode + " and edge " + nextEdge);
134: }
135: }
136: }
137: newNodes.addAll(nodesToAdd);
138: newEdges.addAll(edgesToAdd);
139: }
140:
141: // System.out.println("I've got " + newNodes + " and ");
142: // System.out.println(newEdges);
143:
144: totalBase = totalTarget = 0;
145: Iterator newNodeIt = newNodes.iterator();
146: while (newNodeIt.hasNext()) {
147: DelaunayNode nextNode = (DelaunayNode) newNodeIt
148: .next();
149: Feature nextFeature2 = nextNode.getFeature();
150: Object neighborsBaseObj = base
151: .getValue(nextFeature2);
152: if (!(baseObj instanceof Number)) {
153: throw new RuntimeException(
154: "Expression "
155: + base
156: + " must evaluate to a number on feature "
157: + nextFeature2);
158: }
159: Object neighborsTargetObj = target
160: .getValue(nextFeature2);
161: if (!(targetObj instanceof Number)) {
162: throw new RuntimeException(
163: "Expression "
164: + target
165: + " must evaluate to a number on feature "
166: + nextFeature2);
167: }
168: totalBase = totalBase
169: + ((Number) baseObj).doubleValue();
170: totalTarget = totalTarget
171: + ((Number) targetObj).doubleValue();
172: }
173: }
174: double expectedTarget = meanRate * totalBase;
175:
176: double top = ((Math.pow(Math.E, (0 - expectedTarget)) * (Math
177: .pow(expectedTarget, totalTarget))));
178: double bottom = fact((int) Math.round(totalTarget));
179: double poissonProb = top / bottom;
180: // System.out.println("testing " + newNodes);
181: // System.out.println("testing " + newEdges);
182: System.out.println(next.getCoordinate().x + ", "
183: + next.getCoordinate().y + ", " + totalTarget
184: + ", " + expectedTarget + ", " + poissonProb);
185:
186: if (poissonProb < threshold) {
187: clusterNodes.addAll(newNodes);
188: clusterEdges.addAll(newEdges);
189: }
190: } //end while (nodeIt.hasNext())
191:
192: return new BasicGraph(clusterNodes, clusterEdges);
193:
194: }
195:
196: private static double iterFact(int i, int f) {
197: if ((i == 0) || (i == 1)) {
198: return f;
199: } else {
200: return iterFact(i - 1, i * f);
201: }
202: }
203:
204: public static double fact(int i) {
205: return iterFact(i, 1);
206: }
207:
208: }
|