Mercurial > dive4elements > river

comparison artifacts/src/main/java/org/dive4elements/river/artifacts/math/BackJumpCorrector.java @ 5838:5aa05a7a34b7

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Rename modules to more fitting names.
author Sascha L. Teichmann <teichmann@intevation.de>
date Thu, 25 Apr 2013 15:23:37 +0200
parents flys-artifacts/src/main/java/org/dive4elements/river/artifacts/math/BackJumpCorrector.java@bd047b71ab37
children 4897a58c8746
comparison
equal deleted inserted replaced
5837:d9901a08d0a6 5838:5aa05a7a34b7
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 :
http://dive4elements.wald.intevation.org