teichmann@6152: /* Copyright (C) 2011, 2012, 2013 by Bundesanstalt für Gewässerkunde
teichmann@6152:  * Software engineering by Intevation GmbH
teichmann@6152:  *
teichmann@6152:  * This file is Free Software under the GNU AGPL (>=v3)
teichmann@6152:  * and comes with ABSOLUTELY NO WARRANTY! Check out the
teichmann@6152:  * documentation coming with Dive4Elements River for details.
teichmann@6152:  */
teichmann@6152: 
teichmann@6152: package org.dive4elements.river.artifacts.model;
teichmann@6152: 
teichmann@6152: import gnu.trove.TDoubleArrayList;
teichmann@6152: 
teichmann@6152: import java.io.Serializable;
teichmann@6152: 
teichmann@6152: import java.util.ArrayList;
teichmann@6152: import java.util.List;
teichmann@6152: import java.util.Set;
teichmann@6152: import java.util.TreeSet;
teichmann@6152: 
teichmann@6152: import org.apache.commons.math.stat.descriptive.moment.StandardDeviation;
teichmann@6152: 
teichmann@6152: import org.dive4elements.river.artifacts.math.Linear;
teichmann@6152: 
teichmann@8186: import org.dive4elements.river.backend.utils.EpsilonComparator;
teichmann@6152: 
teichmann@6152: public class LinearInterpolated
teichmann@6152: implements   Serializable
teichmann@6152: {
teichmann@6152:     public static final double EPSILON = 1e-5;
teichmann@6152:     public static final double MULTIPLE_STD_DEV = 4d;
teichmann@6152: 
teichmann@6152:     public static final EpsilonComparator CMP = new EpsilonComparator(EPSILON);
teichmann@6152: 
teichmann@6152:     private TDoubleArrayList xs;
teichmann@6152:     private TDoubleArrayList ys;
teichmann@6152: 
teichmann@6152:     private List<Range> gaps;
teichmann@6152: 
teichmann@6152:     public interface Operator {
teichmann@6152:         double evaluate(double y1, double y2);
teichmann@6152:     } // interface Operator
teichmann@6152: 
teichmann@6152:     public static final Operator SUB = new Operator() {
teichmann@6152:         @Override
teichmann@6152:         public double evaluate(double y1, double y2) {
teichmann@6152:             return y1 - y2;
teichmann@6152:         }
teichmann@6152:     };
teichmann@6152: 
teichmann@6152:     public static final Operator MAX = new Operator() {
teichmann@6152:         @Override
teichmann@6152:         public double evaluate(double y1, double y2) {
teichmann@6152:             return Math.max(y1, y2);
teichmann@6152:         }
teichmann@6152:     };
teichmann@6152: 
teichmann@6152:     public LinearInterpolated() {
teichmann@6152:         xs = new TDoubleArrayList();
teichmann@6152:         ys = new TDoubleArrayList();
teichmann@6152:     }
teichmann@6152: 
teichmann@6152:     public LinearInterpolated(int capacity) {
teichmann@6152:         xs = new TDoubleArrayList(capacity);
teichmann@6152:         ys = new TDoubleArrayList(capacity);
teichmann@6152:     }
teichmann@6152: 
teichmann@6152:     public void add(double x, double y) {
teichmann@6152:         xs.add(x);
teichmann@6152:         ys.add(y);
teichmann@6152:     }
teichmann@6152: 
teichmann@6152:     public int size() {
teichmann@6152:         return xs.size();
teichmann@6152:     }
teichmann@6152: 
teichmann@6152:     public void pointsInRange(double from, double to, Set<Double> points) {
teichmann@6152:         if (from > to) {
teichmann@6152:             double t = from;
teichmann@6152:             from = to;
teichmann@6152:             to = t;
teichmann@6152:         }
teichmann@6152:         for (int i = 0, S = size(); i < S; ++i) {
teichmann@6152:             double x = xs.getQuick(i);
teichmann@6152:             if (x >= from && x <= to) {
teichmann@6152:                 points.add(x);
teichmann@6152:             }
teichmann@6152:         }
teichmann@6152:     }
teichmann@6152: 
teichmann@6152:     public boolean inGap(double x) {
teichmann@6152:         if (gaps != null) {
teichmann@6152:             for (Range gap: gaps) {
teichmann@6152:                 if (gap.inside(x)) {
teichmann@6152:                     return true;
teichmann@6152:                 }
teichmann@6152:             }
teichmann@6152:         }
teichmann@6152:         return false;
teichmann@6152:     }
teichmann@6152: 
teichmann@6152:     public void detectGaps(double threshold) {
teichmann@6152:         List<Range> gabs = new ArrayList<Range>();
teichmann@6152:         for (int i = 1, S = size(); i < S; ++i) {
teichmann@6152:             double x0 = xs.getQuick(i-1);
teichmann@6152:             double x1 = xs.getQuick(i);
teichmann@6152:             double minX, maxX;
teichmann@6152:             if (x0 < x1) { minX = x0; maxX = x1; }
teichmann@6152:             else         { minX = x1; maxX = x0; }
teichmann@6152:             double diff = maxX - minX - 2d*EPSILON;
teichmann@6152:             if (diff > threshold) {
teichmann@6152:                 gaps.add(new Range(minX+EPSILON, maxX-EPSILON));
teichmann@6152:             }
teichmann@6152:         }
teichmann@6152:         this.gaps = gaps.isEmpty() ? null : gabs;
teichmann@6152:     }
teichmann@6152: 
teichmann@6152:     public void resetGaps() {
teichmann@6152:         gaps = null;
teichmann@6152:     }
teichmann@6152: 
teichmann@6152:     public double guessGapThreshold() {
teichmann@6152:         return guessGapThreshold(MULTIPLE_STD_DEV);
teichmann@6152:     }
teichmann@6152: 
teichmann@6152:     public double guessGapThreshold(double scale) {
teichmann@6152:         int S = size();
teichmann@6152:         if (S < 5) { // Too less points.
teichmann@6152:             return Double.MAX_VALUE;
teichmann@6152:         }
teichmann@6152: 
teichmann@6152:         StandardDeviation s = new StandardDeviation();
teichmann@6152: 
teichmann@6152:         for (int i = 1; i < S; ++i) {
teichmann@6152:             double diff = Math.abs(xs.getQuick(i-1) - xs.getQuick(i));
teichmann@6152:             s.increment(diff);
teichmann@6152:         }
teichmann@6152: 
teichmann@6152:         return scale*s.getResult();
teichmann@6152:     }
teichmann@6152: 
teichmann@6152:     public double value(double x) {
teichmann@6152:         for (int i = 0, S = size(); i < S; ++i) {
teichmann@6152:             double x1 = xs.getQuick(i);
teichmann@6152:             if (Math.abs(x1 - x) < EPSILON) {
teichmann@6152:                 return ys.getQuick(i);
teichmann@6152:             }
teichmann@6152:             if (i > 0) {
teichmann@6152:                 double x0 = xs.getQuick(i-1);
teichmann@6152:                 double minX, maxX;
teichmann@6152:                 if (x0 < x1) { minX = x0; maxX = x1; }
teichmann@6152:                 else         { minX = x1; maxX = x0; }
teichmann@6152:                 if (x > minX && x < maxX) {
teichmann@6152:                     return Linear.linear(
teichmann@6152:                         x,
teichmann@6152:                         x0, x1,
teichmann@6152:                         ys.getQuick(i-1), ys.getQuick(i));
teichmann@6152:                 }
teichmann@6152:             }
teichmann@6152:         }
teichmann@6152:         return Double.NaN;
teichmann@6152:     }
teichmann@6152: 
teichmann@6152:     public LinearInterpolated sub(
teichmann@7248:         LinearInterpolated other,
teichmann@6152:         double             from,
teichmann@6152:         double             to
teichmann@6152:     ) {
teichmann@6152:         return apply(SUB, other, from, to);
teichmann@6152:     }
teichmann@6152: 
teichmann@6152:     public LinearInterpolated max(
teichmann@7248:         LinearInterpolated other,
teichmann@6152:         double             from,
teichmann@6152:         double             to
teichmann@6152:     ) {
teichmann@6152:         return apply(MAX, other, from, to);
teichmann@6152:     }
teichmann@6152: 
teichmann@6153:     public boolean intersect(LinearInterpolated other) {
teichmann@6153:         if (xs.isEmpty() || other.xs.isEmpty()) {
teichmann@6153:             return false;
teichmann@6153:         }
teichmann@6153: 
teichmann@6153:         double tMax = xs.max();
teichmann@6153:         double oMin = other.xs.min();
teichmann@6153:         if (tMax < oMin) {
teichmann@6153:             return false;
teichmann@6153:         }
teichmann@6153: 
teichmann@6153:         double tMin = xs.min();
teichmann@6153:         double oMax = other.xs.max();
teichmann@6153:         return !(tMin > oMax);
teichmann@6153:     }
teichmann@6153: 
teichmann@6152:     public LinearInterpolated apply(
teichmann@6152:         Operator           operator,
teichmann@6152:         LinearInterpolated other,
teichmann@6152:         double             from,
teichmann@6152:         double             to
teichmann@6152:     ) {
teichmann@6153:         LinearInterpolated result = new LinearInterpolated();
teichmann@6153:         if (!intersect(other)) {
teichmann@6153:             return result;
teichmann@6153:         }
teichmann@6153: 
teichmann@6152:         Set<Double> points = new TreeSet<Double>(CMP);
teichmann@6152:         points.add(from);
teichmann@6152:         points.add(to);
teichmann@6152: 
teichmann@6152:         this .pointsInRange(from, to, points);
teichmann@6152:         other.pointsInRange(from, to, points);
teichmann@6152: 
teichmann@6152: 
teichmann@6152:         for (double x: points) {
teichmann@6152:             if (!inGap(x) && !other.inGap(x)) {
teichmann@6152:                 double y1 = this .value(x);
teichmann@6152:                 double y2 = other.value(x);
teichmann@6152:                 double y  = operator.evaluate(y1, y2);
teichmann@6152:                 if (!Double.isNaN(y)) {
teichmann@6152:                     result.add(x, y);
teichmann@6152:                 }
teichmann@6152:             }
teichmann@6152:         }
teichmann@6152: 
teichmann@6152:         return result;
teichmann@6152:     }
teichmann@6152: }
teichmann@6152: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :