changeset 6152:0587819960c3

Waterlevel differences & bed height differences: Add new model LinearInterpolated intented to unify the two very similiar calculations. The focus of the current implementation is correctness and not speed! The fact that the data sets more mostly sorted by station is not exploited. Doing so would improve performance significantly.
author Sascha L. Teichmann <teichmann@intevation.de>
date Sun, 02 Jun 2013 17:52:53 +0200
parents 63d1e2a9b311
children 12af732c9d0f
files artifacts/src/main/java/org/dive4elements/river/artifacts/model/LinearInterpolated.java
diffstat 1 files changed, 208 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/LinearInterpolated.java	Sun Jun 02 17:52:53 2013 +0200
@@ -0,0 +1,208 @@
+/* Copyright (C) 2011, 2012, 2013 by Bundesanstalt für Gewässerkunde
+ * Software engineering by Intevation GmbH
+ *
+ * This file is Free Software under the GNU AGPL (>=v3)
+ * and comes with ABSOLUTELY NO WARRANTY! Check out the
+ * documentation coming with Dive4Elements River for details.
+ */
+
+package org.dive4elements.river.artifacts.model;
+
+import gnu.trove.TDoubleArrayList;
+
+import java.io.Serializable;
+
+import java.util.ArrayList;
+import java.util.List;
+import java.util.Set;
+import java.util.TreeSet;
+
+import org.apache.commons.math.stat.descriptive.moment.StandardDeviation;
+
+import org.dive4elements.river.artifacts.math.Linear;
+
+import org.dive4elements.river.utils.EpsilonComparator;
+
+public class LinearInterpolated
+implements   Serializable
+{
+    public static final double EPSILON = 1e-5;
+    public static final double MULTIPLE_STD_DEV = 4d;
+
+    public static final EpsilonComparator CMP = new EpsilonComparator(EPSILON);
+
+    private TDoubleArrayList xs;
+    private TDoubleArrayList ys;
+
+    private List<Range> gaps;
+
+    public interface Operator {
+        double evaluate(double y1, double y2);
+    } // interface Operator
+
+    public static final Operator SUB = new Operator() {
+        @Override
+        public double evaluate(double y1, double y2) {
+            return y1 - y2;
+        }
+    };
+
+    public static final Operator MAX = new Operator() {
+        @Override
+        public double evaluate(double y1, double y2) {
+            return Math.max(y1, y2);
+        }
+    };
+
+    public LinearInterpolated() {
+        xs = new TDoubleArrayList();
+        ys = new TDoubleArrayList();
+    }
+
+    public LinearInterpolated(int capacity) {
+        xs = new TDoubleArrayList(capacity);
+        ys = new TDoubleArrayList(capacity);
+    }
+
+    public void add(double x, double y) {
+        xs.add(x);
+        ys.add(y);
+    }
+
+    public int size() {
+        return xs.size();
+    }
+
+    public void pointsInRange(double from, double to, Set<Double> points) {
+        if (from > to) {
+            double t = from;
+            from = to;
+            to = t;
+        }
+        for (int i = 0, S = size(); i < S; ++i) {
+            double x = xs.getQuick(i);
+            if (x >= from && x <= to) {
+                points.add(x);
+            }
+        }
+    }
+
+    public boolean inGap(double x) {
+        if (gaps != null) {
+            for (Range gap: gaps) {
+                if (gap.inside(x)) {
+                    return true;
+                }
+            }
+        }
+        return false;
+    }
+
+    public void detectGaps(double threshold) {
+        List<Range> gabs = new ArrayList<Range>();
+        for (int i = 1, S = size(); i < S; ++i) {
+            double x0 = xs.getQuick(i-1);
+            double x1 = xs.getQuick(i);
+            double minX, maxX;
+            if (x0 < x1) { minX = x0; maxX = x1; }
+            else         { minX = x1; maxX = x0; }
+            double diff = maxX - minX - 2d*EPSILON;
+            if (diff > threshold) {
+                gaps.add(new Range(minX+EPSILON, maxX-EPSILON));
+            }
+        }
+        this.gaps = gaps.isEmpty() ? null : gabs;
+    }
+
+    public void resetGaps() {
+        gaps = null;
+    }
+
+    public double guessGapThreshold() {
+        return guessGapThreshold(MULTIPLE_STD_DEV);
+    }
+
+    public double guessGapThreshold(double scale) {
+        int S = size();
+        if (S < 5) { // Too less points.
+            return Double.MAX_VALUE;
+        }
+
+        StandardDeviation s = new StandardDeviation();
+
+        for (int i = 1; i < S; ++i) {
+            double diff = Math.abs(xs.getQuick(i-1) - xs.getQuick(i));
+            s.increment(diff);
+        }
+
+        return scale*s.getResult();
+    }
+
+    public double value(double x) {
+        for (int i = 0, S = size(); i < S; ++i) {
+            double x1 = xs.getQuick(i);
+            if (Math.abs(x1 - x) < EPSILON) {
+                return ys.getQuick(i);
+            }
+            if (i > 0) {
+                double x0 = xs.getQuick(i-1);
+                double minX, maxX;
+                if (x0 < x1) { minX = x0; maxX = x1; }
+                else         { minX = x1; maxX = x0; }
+                if (x > minX && x < maxX) {
+                    return Linear.linear(
+                        x,
+                        x0, x1,
+                        ys.getQuick(i-1), ys.getQuick(i));
+                }
+            }
+        }
+        return Double.NaN;
+    }
+
+    public LinearInterpolated sub(
+        LinearInterpolated other, 
+        double             from,
+        double             to
+    ) {
+        return apply(SUB, other, from, to);
+    }
+
+    public LinearInterpolated max(
+        LinearInterpolated other, 
+        double             from,
+        double             to
+    ) {
+        return apply(MAX, other, from, to);
+    }
+
+    public LinearInterpolated apply(
+        Operator           operator,
+        LinearInterpolated other,
+        double             from,
+        double             to
+    ) {
+        Set<Double> points = new TreeSet<Double>(CMP);
+        points.add(from);
+        points.add(to);
+
+        this .pointsInRange(from, to, points);
+        other.pointsInRange(from, to, points);
+
+        LinearInterpolated result = new LinearInterpolated();
+
+        for (double x: points) {
+            if (!inGap(x) && !other.inGap(x)) {
+                double y1 = this .value(x);
+                double y2 = other.value(x);
+                double y  = operator.evaluate(y1, y2);
+                if (!Double.isNaN(y)) {
+                    result.add(x, y);
+                }
+            }
+        }
+
+        return result;
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :

http://dive4elements.wald.intevation.org