Mercurial > dive4elements > river
comparison flys-artifacts/src/main/java/org/dive4elements/river/artifacts/math/BackJumpCorrector.java @ 5831:bd047b71ab37
Repaired internal references
author | Sascha L. Teichmann <teichmann@intevation.de> |
---|---|
date | Thu, 25 Apr 2013 12:06:39 +0200 |
parents | flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/BackJumpCorrector.java@5642a83420f2 |
children |
comparison
equal
deleted
inserted
replaced
5830:160f53ee0870 | 5831:bd047b71ab37 |
---|---|
1 package org.dive4elements.river.artifacts.math; | |
2 | |
3 import java.util.ArrayList; | |
4 import java.util.List; | |
5 | |
6 import java.io.Serializable; | |
7 | |
8 import org.apache.commons.math.analysis.interpolation.SplineInterpolator; | |
9 | |
10 import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction; | |
11 | |
12 import org.apache.commons.math.ArgumentOutsideDomainException; | |
13 | |
14 import org.apache.commons.math.exception.MathIllegalArgumentException; | |
15 | |
16 import org.apache.log4j.Logger; | |
17 | |
18 import org.dive4elements.river.artifacts.model.Calculation; | |
19 | |
20 import org.dive4elements.river.utils.DoubleUtil; | |
21 | |
22 public class BackJumpCorrector | |
23 implements Serializable | |
24 { | |
25 private static Logger log = Logger.getLogger(BackJumpCorrector.class); | |
26 | |
27 protected ArrayList<Double> backjumps; | |
28 | |
29 protected double [] corrected; | |
30 | |
31 public BackJumpCorrector() { | |
32 backjumps = new ArrayList<Double>(); | |
33 } | |
34 | |
35 public boolean hasBackJumps() { | |
36 return !backjumps.isEmpty(); | |
37 } | |
38 | |
39 public List<Double> getBackJumps() { | |
40 return backjumps; | |
41 } | |
42 | |
43 public double [] getCorrected() { | |
44 return corrected; | |
45 } | |
46 | |
47 public boolean doCorrection( | |
48 double [] km, | |
49 double [] ws, | |
50 Calculation errors | |
51 ) { | |
52 boolean wsUp = DoubleUtil.isIncreasing(ws); | |
53 | |
54 if (wsUp) { | |
55 km = DoubleUtil.swapClone(km); | |
56 ws = DoubleUtil.swapClone(ws); | |
57 } | |
58 | |
59 boolean kmUp = DoubleUtil.isIncreasing(km); | |
60 | |
61 if (!kmUp) { | |
62 km = DoubleUtil.sumDiffs(km); | |
63 } | |
64 | |
65 if (log.isDebugEnabled()) { | |
66 log.debug("BackJumpCorrector.doCorrection ------- enter"); | |
67 log.debug(" km increasing: " + DoubleUtil.isIncreasing(km)); | |
68 log.debug(" ws increasing: " + DoubleUtil.isIncreasing(ws)); | |
69 log.debug("BackJumpCorrector.doCorrection ------- leave"); | |
70 } | |
71 | |
72 boolean hasBackJumps = doCorrectionClean(km, ws, errors); | |
73 | |
74 if (hasBackJumps && wsUp) { | |
75 // mirror back | |
76 DoubleUtil.swap(corrected); | |
77 } | |
78 | |
79 return hasBackJumps; | |
80 } | |
81 | |
82 protected boolean doCorrectionClean( | |
83 double [] km, | |
84 double [] ws, | |
85 Calculation errors | |
86 ) { | |
87 int N = km.length; | |
88 | |
89 if (N != ws.length) { | |
90 throw new IllegalArgumentException("km.length != ws.length"); | |
91 } | |
92 | |
93 if (N < 2) { | |
94 return false; | |
95 } | |
96 | |
97 SplineInterpolator interpolator = null; | |
98 | |
99 for (int i = 1; i < N; ++i) { | |
100 if (ws[i] <= ws[i-1]) { | |
101 // no back jump | |
102 continue; | |
103 } | |
104 backjumps.add(km[i]); | |
105 | |
106 if (corrected == null) { | |
107 // lazy cloning | |
108 ws = corrected = (double [])ws.clone(); | |
109 } | |
110 | |
111 double above = aboveWaterKM(km, ws, i); | |
112 | |
113 if (Double.isNaN(above)) { // run over start km | |
114 // fill all previous | |
115 for (int j = 0; j < i; ++j) { | |
116 ws[j] = ws[i]; | |
117 } | |
118 continue; | |
119 } | |
120 | |
121 double distance = Math.abs(km[i] - above); | |
122 | |
123 double quarterDistance = 0.25*distance; | |
124 | |
125 double start = above - quarterDistance; | |
126 | |
127 double startHeight = DoubleUtil.interpolateSorted(km, ws, start); | |
128 | |
129 if (Double.isNaN(startHeight)) { | |
130 // run over start km | |
131 startHeight = ws[0]; | |
132 } | |
133 | |
134 double between = above + quarterDistance; | |
135 | |
136 double aboveHeight = ws[i] + 0.25*(startHeight - ws[i]); | |
137 | |
138 double [] x = { start, above, between }; | |
139 double [] y = { startHeight, aboveHeight, ws[i] }; | |
140 | |
141 if (log.isDebugEnabled()) { | |
142 for (int j = 0; j < x.length; ++j) { | |
143 log.debug(" " + x[j] + " -> " + y[j]); | |
144 } | |
145 } | |
146 | |
147 if (interpolator == null) { | |
148 interpolator = new SplineInterpolator(); | |
149 } | |
150 | |
151 PolynomialSplineFunction spline; | |
152 | |
153 try { | |
154 spline = interpolator.interpolate(x, y); | |
155 } | |
156 catch (MathIllegalArgumentException miae) { | |
157 errors.addProblem("spline.creation.failed"); | |
158 log.error(miae); | |
159 continue; | |
160 } | |
161 | |
162 try { | |
163 if (log.isDebugEnabled()) { | |
164 log.debug("spline points:"); | |
165 for (int j = 0; j < x.length; ++j) { | |
166 log.debug(x[j] + " " + y[j] + " " + spline.value(x[j])); | |
167 } | |
168 } | |
169 | |
170 int j = i-1; | |
171 | |
172 for (; j >= 0 && km[j] >= between; --j) { | |
173 ws[j] = ws[i]; | |
174 } | |
175 | |
176 for (; j >= 0 && ws[j] < startHeight; --j) { | |
177 ws[j] = spline.value(km[j]); | |
178 } | |
179 } | |
180 catch (ArgumentOutsideDomainException aode) { | |
181 errors.addProblem("spline.interpolation.failed"); | |
182 log.error("spline interpolation failed", aode); | |
183 } | |
184 } // for all km | |
185 | |
186 return !backjumps.isEmpty(); | |
187 } | |
188 | |
189 | |
190 protected static double aboveWaterKM( | |
191 double [] km, | |
192 double [] ws, | |
193 int wIndex | |
194 ) { | |
195 double w = ws[wIndex]; | |
196 | |
197 while (--wIndex >= 0) { | |
198 // still under water | |
199 if (ws[wIndex] < w) continue; | |
200 | |
201 if (ws[wIndex] > w) { | |
202 // f(ws[wIndex]) = km[wIndex] | |
203 // f(ws[wIndex+1]) = km[wIndex+1] | |
204 return Linear.linear( | |
205 w, | |
206 ws[wIndex], ws[wIndex+1], | |
207 km[wIndex], km[wIndex+1]); | |
208 } | |
209 else { | |
210 return km[wIndex]; | |
211 } | |
212 } | |
213 | |
214 return Double.NaN; | |
215 } | |
216 } | |
217 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 : |