Mercurial > dive4elements > river
comparison artifacts/src/main/java/org/dive4elements/river/artifacts/model/LinearInterpolated.java @ 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 | |
children | 12af732c9d0f |
comparison
equal
deleted
inserted
replaced
6151:63d1e2a9b311 | 6152:0587819960c3 |
---|---|
1 /* Copyright (C) 2011, 2012, 2013 by Bundesanstalt für Gewässerkunde | |
2 * Software engineering by Intevation GmbH | |
3 * | |
4 * This file is Free Software under the GNU AGPL (>=v3) | |
5 * and comes with ABSOLUTELY NO WARRANTY! Check out the | |
6 * documentation coming with Dive4Elements River for details. | |
7 */ | |
8 | |
9 package org.dive4elements.river.artifacts.model; | |
10 | |
11 import gnu.trove.TDoubleArrayList; | |
12 | |
13 import java.io.Serializable; | |
14 | |
15 import java.util.ArrayList; | |
16 import java.util.List; | |
17 import java.util.Set; | |
18 import java.util.TreeSet; | |
19 | |
20 import org.apache.commons.math.stat.descriptive.moment.StandardDeviation; | |
21 | |
22 import org.dive4elements.river.artifacts.math.Linear; | |
23 | |
24 import org.dive4elements.river.utils.EpsilonComparator; | |
25 | |
26 public class LinearInterpolated | |
27 implements Serializable | |
28 { | |
29 public static final double EPSILON = 1e-5; | |
30 public static final double MULTIPLE_STD_DEV = 4d; | |
31 | |
32 public static final EpsilonComparator CMP = new EpsilonComparator(EPSILON); | |
33 | |
34 private TDoubleArrayList xs; | |
35 private TDoubleArrayList ys; | |
36 | |
37 private List<Range> gaps; | |
38 | |
39 public interface Operator { | |
40 double evaluate(double y1, double y2); | |
41 } // interface Operator | |
42 | |
43 public static final Operator SUB = new Operator() { | |
44 @Override | |
45 public double evaluate(double y1, double y2) { | |
46 return y1 - y2; | |
47 } | |
48 }; | |
49 | |
50 public static final Operator MAX = new Operator() { | |
51 @Override | |
52 public double evaluate(double y1, double y2) { | |
53 return Math.max(y1, y2); | |
54 } | |
55 }; | |
56 | |
57 public LinearInterpolated() { | |
58 xs = new TDoubleArrayList(); | |
59 ys = new TDoubleArrayList(); | |
60 } | |
61 | |
62 public LinearInterpolated(int capacity) { | |
63 xs = new TDoubleArrayList(capacity); | |
64 ys = new TDoubleArrayList(capacity); | |
65 } | |
66 | |
67 public void add(double x, double y) { | |
68 xs.add(x); | |
69 ys.add(y); | |
70 } | |
71 | |
72 public int size() { | |
73 return xs.size(); | |
74 } | |
75 | |
76 public void pointsInRange(double from, double to, Set<Double> points) { | |
77 if (from > to) { | |
78 double t = from; | |
79 from = to; | |
80 to = t; | |
81 } | |
82 for (int i = 0, S = size(); i < S; ++i) { | |
83 double x = xs.getQuick(i); | |
84 if (x >= from && x <= to) { | |
85 points.add(x); | |
86 } | |
87 } | |
88 } | |
89 | |
90 public boolean inGap(double x) { | |
91 if (gaps != null) { | |
92 for (Range gap: gaps) { | |
93 if (gap.inside(x)) { | |
94 return true; | |
95 } | |
96 } | |
97 } | |
98 return false; | |
99 } | |
100 | |
101 public void detectGaps(double threshold) { | |
102 List<Range> gabs = new ArrayList<Range>(); | |
103 for (int i = 1, S = size(); i < S; ++i) { | |
104 double x0 = xs.getQuick(i-1); | |
105 double x1 = xs.getQuick(i); | |
106 double minX, maxX; | |
107 if (x0 < x1) { minX = x0; maxX = x1; } | |
108 else { minX = x1; maxX = x0; } | |
109 double diff = maxX - minX - 2d*EPSILON; | |
110 if (diff > threshold) { | |
111 gaps.add(new Range(minX+EPSILON, maxX-EPSILON)); | |
112 } | |
113 } | |
114 this.gaps = gaps.isEmpty() ? null : gabs; | |
115 } | |
116 | |
117 public void resetGaps() { | |
118 gaps = null; | |
119 } | |
120 | |
121 public double guessGapThreshold() { | |
122 return guessGapThreshold(MULTIPLE_STD_DEV); | |
123 } | |
124 | |
125 public double guessGapThreshold(double scale) { | |
126 int S = size(); | |
127 if (S < 5) { // Too less points. | |
128 return Double.MAX_VALUE; | |
129 } | |
130 | |
131 StandardDeviation s = new StandardDeviation(); | |
132 | |
133 for (int i = 1; i < S; ++i) { | |
134 double diff = Math.abs(xs.getQuick(i-1) - xs.getQuick(i)); | |
135 s.increment(diff); | |
136 } | |
137 | |
138 return scale*s.getResult(); | |
139 } | |
140 | |
141 public double value(double x) { | |
142 for (int i = 0, S = size(); i < S; ++i) { | |
143 double x1 = xs.getQuick(i); | |
144 if (Math.abs(x1 - x) < EPSILON) { | |
145 return ys.getQuick(i); | |
146 } | |
147 if (i > 0) { | |
148 double x0 = xs.getQuick(i-1); | |
149 double minX, maxX; | |
150 if (x0 < x1) { minX = x0; maxX = x1; } | |
151 else { minX = x1; maxX = x0; } | |
152 if (x > minX && x < maxX) { | |
153 return Linear.linear( | |
154 x, | |
155 x0, x1, | |
156 ys.getQuick(i-1), ys.getQuick(i)); | |
157 } | |
158 } | |
159 } | |
160 return Double.NaN; | |
161 } | |
162 | |
163 public LinearInterpolated sub( | |
164 LinearInterpolated other, | |
165 double from, | |
166 double to | |
167 ) { | |
168 return apply(SUB, other, from, to); | |
169 } | |
170 | |
171 public LinearInterpolated max( | |
172 LinearInterpolated other, | |
173 double from, | |
174 double to | |
175 ) { | |
176 return apply(MAX, other, from, to); | |
177 } | |
178 | |
179 public LinearInterpolated apply( | |
180 Operator operator, | |
181 LinearInterpolated other, | |
182 double from, | |
183 double to | |
184 ) { | |
185 Set<Double> points = new TreeSet<Double>(CMP); | |
186 points.add(from); | |
187 points.add(to); | |
188 | |
189 this .pointsInRange(from, to, points); | |
190 other.pointsInRange(from, to, points); | |
191 | |
192 LinearInterpolated result = new LinearInterpolated(); | |
193 | |
194 for (double x: points) { | |
195 if (!inGap(x) && !other.inGap(x)) { | |
196 double y1 = this .value(x); | |
197 double y2 = other.value(x); | |
198 double y = operator.evaluate(y1, y2); | |
199 if (!Double.isNaN(y)) { | |
200 result.add(x, y); | |
201 } | |
202 } | |
203 } | |
204 | |
205 return result; | |
206 } | |
207 } | |
208 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 : |