sascha@625: package de.intevation.flys.artifacts.model; sascha@625: sascha@625: import java.io.Serializable; sascha@625: sascha@625: import java.util.List; sascha@2611: import java.util.ArrayList; sascha@625: sascha@633: import org.apache.log4j.Logger; sascha@625: sascha@625: public class QRangeTree sascha@625: implements Serializable sascha@625: { sascha@633: private static Logger log = Logger.getLogger(QRangeTree.class); sascha@633: sascha@2611: public static final double EPSILON = 1e-4; sascha@2611: sascha@625: public static class Node sascha@625: implements Serializable sascha@625: { sascha@625: Node left; sascha@625: Node right; sascha@625: Node prev; sascha@625: Node next; sascha@625: sascha@625: double a; sascha@625: double b; sascha@625: double q; sascha@625: sascha@625: public Node() { sascha@625: } sascha@625: sascha@625: public Node(double a, double b, double q) { sascha@625: this.a = a; sascha@625: this.b = b; sascha@625: this.q = q; sascha@625: } sascha@625: sascha@625: protected final double interpolatePrev(double pos) { sascha@625: /* sascha@625: f(prev.b) = prev.q sascha@625: f(a) = q sascha@625: sascha@625: prev.q = m*prev.b + n sascha@625: q = m*a + n <=> n = q - m*a sascha@625: sascha@625: q - prev.q = m*(a - prev.b) sascha@625: sascha@625: m = (q - prev.q)/(a - prev.b) # a != prev.b sascha@625: */ sascha@625: sascha@625: if (a == prev.b) { sascha@625: return 0.5*(q + prev.q); sascha@625: } sascha@625: double m = (q - prev.q)/(a - prev.b); sascha@625: double n = q - m*a; sascha@625: return m*pos + n; sascha@625: } sascha@625: sascha@625: protected final double interpolateNext(double pos) { sascha@625: /* sascha@625: f(next.a) = next.q sascha@625: f(b) = q sascha@625: sascha@625: next.q = m*next.a + n sascha@625: q = m*b + n <=> n = q - m*b sascha@625: sascha@625: q - next.q = m*(b - next.a) sascha@625: m = (q - next.q)/(b - next.a) # b != next.a sascha@625: */ sascha@625: sascha@625: if (b == next.a) { sascha@625: return 0.5*(q + next.q); sascha@625: } sascha@625: double m = (q - next.q)/(b - next.a); sascha@625: double n = q - m*b; sascha@625: return m*pos + n; sascha@625: } sascha@625: sascha@625: public double findQ(double pos) { sascha@633: sascha@625: Node current = this; sascha@625: for (;;) { sascha@625: if (pos < current.a) { sascha@625: if (current.left != null) { sascha@625: current = current.left; sascha@625: continue; sascha@625: } sascha@633: return current.prev != null sascha@625: ? current.interpolatePrev(pos) sascha@625: : Double.NaN; sascha@625: } sascha@625: if (pos > current.b) { sascha@625: if (current.right != null) { sascha@625: current = current.right; sascha@625: continue; sascha@625: } sascha@633: return current.next != null sascha@625: ? current.interpolateNext(pos) sascha@625: : Double.NaN; sascha@625: } sascha@625: return current.q; sascha@625: } sascha@625: } teichmann@4797: teichmann@4797: public Node findNode(double pos) { teichmann@4797: Node current = this; teichmann@4797: while (current != null) { teichmann@4797: if (pos < current.a) { teichmann@4797: current = current.left; teichmann@4797: } teichmann@4797: else if (pos > current.b) { teichmann@4797: current = current.right; teichmann@4797: } teichmann@4797: return current; teichmann@4797: } teichmann@4797: return null; teichmann@4797: } teichmann@4797: teichmann@4797: public boolean contains(double c) { teichmann@4797: return c >= a && c <= b; teichmann@4797: } sascha@625: } // class Node sascha@625: teichmann@4797: /** Class to cache the last found tree leaf in a search for Q. teichmann@4797: * Its likely that a neighbored pos search teichmann@4797: * results in using the same leaf node. So teichmann@4797: * caching this leaf will minimize expensive teichmann@4797: * tree traversals. teichmann@4797: * Modeled as inner class because the QRangeTree teichmann@4797: * itself is a shared data structure. teichmann@4797: * Using this class omits interpolation between teichmann@4797: * leaves. teichmann@4797: */ teichmann@4797: public final class QuickQFinder { teichmann@4797: teichmann@4797: private Node last; teichmann@4797: teichmann@4797: public QuickQFinder() { teichmann@4797: } teichmann@4797: teichmann@4797: public double findQ(double pos) { teichmann@4797: if (last != null && last.contains(pos)) { teichmann@4797: return last.q; teichmann@4797: } teichmann@4797: last = QRangeTree.this.findNode(pos); teichmann@4797: return last != null ? last.q : Double.NaN; teichmann@4797: } teichmann@4821: teichmann@4821: public double [] findQs(double [] kms, Calculation report) { teichmann@4821: return findQs(kms, new double[kms.length], report); teichmann@4821: } teichmann@4821: teichmann@4821: public double [] findQs( teichmann@5282: double [] kms, teichmann@5282: double [] qs, teichmann@4821: Calculation report teichmann@4821: ) { teichmann@4821: for (int i = 0; i < kms.length; ++i) { teichmann@4821: if (Double.isNaN(qs[i] = findQ(kms[i]))) { teichmann@4821: report.addProblem(kms[i], "cannot.find.q"); teichmann@4821: } teichmann@4821: } teichmann@4821: return qs; teichmann@4821: } teichmann@4797: } // class QuickQFinder teichmann@4797: sascha@625: protected Node root; sascha@625: sascha@625: public QRangeTree() { sascha@625: } sascha@625: sascha@2609: public static final class AccessQAB { sascha@2609: private int startIndex; sascha@2609: sascha@2609: public AccessQAB(int startIndex) { sascha@2609: this.startIndex = startIndex; sascha@2609: } sascha@2609: sascha@2609: public Double getQ(Object [] row) { sascha@2609: return (Double)row[startIndex]; sascha@2609: } sascha@2609: sascha@2609: public Double getA(Object [] row) { sascha@2609: return (Double)row[startIndex+1]; sascha@2609: } sascha@2609: sascha@2609: public Double getB(Object [] row) { sascha@2609: return (Double)row[startIndex+2]; sascha@2609: } sascha@2609: } sascha@2609: sascha@2609: public static final AccessQAB WITH_COLUMN = new AccessQAB(1); sascha@2609: public static final AccessQAB WITHOUT_COLUMN = new AccessQAB(0); sascha@2609: sascha@632: /** wstQRanges need to be sorted by range.a */ sascha@632: public QRangeTree(List qRanges, int start, int stop) { sascha@2609: this(qRanges, WITH_COLUMN, start, stop); sascha@2609: } sascha@2609: sascha@2609: public QRangeTree( sascha@2609: List qRanges, sascha@3076: AccessQAB accessQAB, sascha@2609: int start, sascha@2609: int stop sascha@2609: ) { sascha@632: if (stop <= start) { sascha@632: return; sascha@632: } sascha@632: sascha@2611: int N = stop-start; sascha@2611: sascha@2611: List nodes = new ArrayList(N); sascha@2611: sascha@2611: Node last = null; sascha@2611: sascha@2611: for (int i = 0; i < N; ++i) { felix@1890: Object [] qRange = qRanges.get(start + i); sascha@2609: Double q = accessQAB.getQ(qRange); sascha@2609: Double a = accessQAB.getA(qRange); sascha@2609: Double b = accessQAB.getB(qRange); sascha@632: sascha@2611: double av = a != null ? a.doubleValue() : -Double.MAX_VALUE; sascha@2611: double bv = b != null ? b.doubleValue() : Double.MAX_VALUE; sascha@2611: double qv = q.doubleValue(); sascha@2611: sascha@2611: // If nodes are directly neighbored and Qs are the same sascha@2611: // join them. sascha@2611: if (last != null sascha@2611: && Math.abs(last.b - av) < EPSILON sascha@2611: && Math.abs(last.q - qv) < EPSILON) { sascha@2611: last.b = bv; sascha@2611: } sascha@2611: else { sascha@2611: nodes.add(last = new Node(av, bv, qv)); sascha@2611: } sascha@632: } sascha@632: sascha@2612: if (log.isDebugEnabled()) { sascha@2612: log.debug("Before/after nodes join: " + sascha@2612: N + "/" + nodes.size()); sascha@2612: } sascha@2612: sascha@632: root = wireTree(nodes); sascha@632: } sascha@632: sascha@2611: protected static Node wireTree(List nodes) { sascha@2611: int N = nodes.size(); sascha@2611: for (int i = 0; i < N; ++i) { sascha@2611: Node node = nodes.get(i); sascha@2611: if (i > 0 ) node.prev = nodes.get(i-1); sascha@2611: if (i < N-1) node.next = nodes.get(i+1); sascha@625: } sascha@625: sascha@2611: return buildTree(nodes, 0, N-1); sascha@625: } sascha@625: sascha@2611: protected static Node buildTree(List nodes, int lo, int hi) { sascha@632: sascha@632: if (lo > hi) { sascha@632: return null; sascha@632: } sascha@632: sascha@625: int mid = (lo + hi) >> 1; sascha@2611: Node parent = nodes.get(mid); sascha@625: sascha@632: parent.left = buildTree(nodes, lo, mid-1); sascha@632: parent.right = buildTree(nodes, mid+1, hi); sascha@625: sascha@625: return parent; sascha@625: } sascha@625: aheinecke@5423: public double averageQ() { aheinecke@5423: double sum = 0d; aheinecke@5423: int n = 0; aheinecke@5423: for (Node node = head(); node != null; node = node.next) { aheinecke@5423: sum += node.q; aheinecke@5423: ++n; aheinecke@5423: } aheinecke@5423: return sum/n; aheinecke@5423: } aheinecke@5423: aheinecke@5423: public double maxQ() { aheinecke@5423: double max = -Double.MAX_VALUE; aheinecke@5423: for (Node node = head(); node != null; node = node.next) { aheinecke@5423: if (node.q > max) { aheinecke@5423: max = node.q; aheinecke@5423: } aheinecke@5423: } aheinecke@5423: return max; aheinecke@5423: } aheinecke@5423: sascha@625: public double findQ(double pos) { sascha@625: return root != null ? root.findQ(pos) : Double.NaN; sascha@625: } sascha@632: teichmann@4797: public Node findNode(double pos) { teichmann@4797: return root != null ? root.findNode(pos) : null; teichmann@4797: } teichmann@4797: sascha@3743: protected Node head() { sascha@3743: Node head = root; sascha@3743: while (head.left != null) { sascha@3743: head = head.left; sascha@3743: } sascha@3743: return head; sascha@3743: } sascha@3743: teichmann@4796: public boolean intersectsQRange(double qMin, double qMax) { teichmann@4796: if (qMin > qMax) { teichmann@4796: double t = qMin; teichmann@4796: qMin = qMax; teichmann@4796: qMax = t; teichmann@4796: } teichmann@4796: for (Node curr = head(); curr != null; curr = curr.next) { teichmann@4796: if (curr.q >= qMin || curr.q <= qMax) { teichmann@4796: return true; teichmann@4796: } teichmann@4796: } teichmann@4796: return false; teichmann@4796: } teichmann@4796: sascha@3743: public List findSegments(double a, double b) { sascha@3743: if (a > b) { double t = a; a = b; b = t; } sascha@3743: return findSegments(new Range(a, b)); sascha@3743: } sascha@3743: sascha@3743: public List findSegments(Range range) { sascha@3743: List segments = new ArrayList(); sascha@3743: sascha@3743: // Linear scan should be good enough here. sascha@3743: for (Node curr = head(); curr != null; curr = curr.next) { sascha@3743: if (!range.disjoint(curr.a, curr.b)) { sascha@3743: Range r = new Range(curr.a, curr.b); sascha@3743: if (r.clip(range)) { sascha@3743: segments.add(r); sascha@3743: } sascha@3743: } sascha@3743: } sascha@3743: sascha@3743: return segments; sascha@3743: } sascha@3743: ingo@2792: @Override ingo@2792: public String toString() { ingo@2792: StringBuilder sb = new StringBuilder(); ingo@2792: inorder(root, sb); ingo@2792: return sb.toString(); ingo@2792: } ingo@2792: ingo@2792: protected static void inorder(Node node, StringBuilder sb) { ingo@2792: if (node != null) { ingo@2792: inorder(node.left, sb); ingo@2792: sb.append('[') ingo@2792: .append(node.a) ingo@2792: .append(", ") ingo@2792: .append(node.b) ingo@2792: .append(": ") ingo@2792: .append(node.q) ingo@2792: .append(']'); ingo@2792: inorder(node.right, sb); ingo@2792: } ingo@2792: } ingo@2792: sascha@632: private static final String name(Object o) { sascha@632: return String.valueOf(System.identityHashCode(o) & 0xffffffffL); sascha@632: } sascha@632: sascha@632: public String toGraph() { sascha@632: StringBuilder sb = new StringBuilder(); sascha@632: sb.append("subgraph c"); sascha@632: sb.append(name(this)); sascha@632: sb.append(" {\n"); sascha@632: if (root != null) { sascha@1020: java.util.Deque stack = new java.util.ArrayDeque(); sascha@632: stack.push(root); sascha@1020: while (!stack.isEmpty()) { sascha@632: Node current = stack.pop(); sascha@632: String name = "n" + name(current); sascha@632: sb.append(name); sascha@632: sb.append(" [label=\""); sascha@632: sb.append(current.a).append(", ").append(current.b); sascha@632: sb.append(": ").append(current.q).append("\"]\n"); sascha@632: if (current.left != null) { sascha@632: String leftName = name(current.left); sascha@632: sb.append(name).append(" -- n").append(leftName).append("\n"); sascha@632: stack.push(current.left); sascha@632: } sascha@632: if (current.right != null) { sascha@632: String rightName = name(current.right); sascha@632: sb.append(name).append(" -- n").append(rightName).append("\n"); sascha@632: stack.push(current.right); sascha@632: } sascha@632: } sascha@632: } sascha@632: sb.append("}\n"); sascha@632: return sb.toString(); sascha@632: } sascha@625: } sascha@625: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :