Mercurial > dive4elements > river

comparison flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/ExtremeCalculation.java @ 4259:5cc9453456a7

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First complete but untested version of the 'Auslagerung extremer Wasserspiegellagen' calculation.
author Sascha L. Teichmann <teichmann@intevation.de>
date Thu, 25 Oct 2012 17:25:37 +0200
parents c63f0b4ac1b4
children 10c1a413a1e0
comparison
equal deleted inserted replaced
4258:2c6e571f366a 4259:5cc9453456a7
1 package de.intevation.flys.artifacts.model.extreme; 1 package de.intevation.flys.artifacts.model.extreme;
2 2
3 import org.apache.commons.logging.Log;
4 import org.apache.commons.logging.LogFactory;
5
6 import org.apache.commons.math.MathException;
7
8 import org.apache.commons.math.optimization.fitting.CurveFitter;
9
10 import org.apache.commons.math.optimization.general.LevenbergMarquardtOptimizer;
11
3 import de.intevation.flys.artifacts.access.ExtremeAccess; 12 import de.intevation.flys.artifacts.access.ExtremeAccess;
13
14 import de.intevation.flys.artifacts.math.Utils;
4 15
5 import de.intevation.flys.artifacts.math.fitting.Function; 16 import de.intevation.flys.artifacts.math.fitting.Function;
6 import de.intevation.flys.artifacts.math.fitting.FunctionFactory; 17 import de.intevation.flys.artifacts.math.fitting.FunctionFactory;
7 18
8 import de.intevation.flys.artifacts.model.Calculation; 19 import de.intevation.flys.artifacts.model.Calculation;
9 import de.intevation.flys.artifacts.model.CalculationResult; 20 import de.intevation.flys.artifacts.model.CalculationResult;
10 import de.intevation.flys.artifacts.model.RangeWithValues; 21 import de.intevation.flys.artifacts.model.RangeWithValues;
11 import de.intevation.flys.artifacts.model.RiverFactory; 22 import de.intevation.flys.artifacts.model.RiverFactory;
23 import de.intevation.flys.artifacts.model.WQKms;
12 import de.intevation.flys.artifacts.model.WstValueTable; 24 import de.intevation.flys.artifacts.model.WstValueTable;
13 import de.intevation.flys.artifacts.model.WstValueTableFactory; 25 import de.intevation.flys.artifacts.model.WstValueTableFactory;
14 26
15 import de.intevation.flys.model.River; 27 import de.intevation.flys.model.River;
16 28
17 import de.intevation.flys.utils.DoubleUtil; 29 import de.intevation.flys.utils.DoubleUtil;
30 import de.intevation.flys.utils.KMIndex;
31
32 import gnu.trove.TDoubleArrayList;
18 33
19 import java.util.List; 34 import java.util.List;
20 35
21 /** Calculate extrapolated W. */ 36 /** Calculate extrapolated W. */
22 public class ExtremeCalculation 37 public class ExtremeCalculation
23 extends Calculation 38 extends Calculation
24 { 39 {
40 private static final Log log =
41 LogFactory.getLog(ExtremeCalculation.class);
42
25 protected String river; 43 protected String river;
26 protected String function; 44 protected String function;
27 protected double from; 45 protected double from;
28 protected double to; 46 protected double to;
29 protected double step; 47 protected double step;
117 return hasProblems() 135 return hasProblems()
118 ? new CalculationResult(this) 136 ? new CalculationResult(this)
119 : innerCalculate(wst, function); 137 : innerCalculate(wst, function);
120 } 138 }
121 139
140 protected String wqkmsName(int i) {
141 StringBuilder sb = new StringBuilder("W(");
142 boolean already = false;
143 for (RangeWithValues r: ranges) {
144 double [] values = r.getValues();
145 if (i < values.length) {
146 if (already) {
147 sb.append(", ");
148 }
149 else {
150 already = true;
151 }
152 // TODO: i18n
153 sb.append(values[i]);
154 }
155 }
156 return sb.append(')').toString();
157 }
158
159 protected WQKms [] allocWQKms() {
160 int max = 0;
161 for (RangeWithValues r: ranges) {
162 double [] values = r.getValues();
163 if (values.length > max) {
164 max = values.length;
165 }
166 }
167 WQKms [] wqkms = new WQKms[max];
168 for (int i = 0; i < max; ++i) {
169 wqkms[i] = new WQKms(wqkmsName(i));
170 }
171 return wqkms;
172 }
173
122 174
123 /** Calculate an extreme curve (extrapolate). */ 175 /** Calculate an extreme curve (extrapolate). */
124 protected CalculationResult innerCalculate( 176 protected CalculationResult innerCalculate(
125 WstValueTable wst, 177 WstValueTable wst,
126 Function function 178 Function function
127 ) { 179 ) {
128 RangeWithValues range = null; 180 RangeWithValues range = null;
129 181
130 KMs: for (double km = from; km <= to; km += step) { 182 double [] chiSqr = { 0d };
183
184 KMIndex<Curve> curves = new KMIndex<Curve>();
185 WQKms [] wqkms = allocWQKms();
186
187 for (double km = from; km <= to; km += step) {
131 double currentKm = DoubleUtil.round(km); 188 double currentKm = DoubleUtil.round(km);
132 189
133 if (range == null || !range.inside(currentKm)) { 190 if (range == null || !range.inside(currentKm)) {
134 for (RangeWithValues r: ranges) { 191 for (RangeWithValues r: ranges) {
135 if (r.inside(currentKm)) { 192 if (r.inside(currentKm)) {
137 break; 194 break;
138 } 195 }
139 } 196 }
140 // TODO: i18n 197 // TODO: i18n
141 addProblem(currentKm, "extreme.no.range"); 198 addProblem(currentKm, "extreme.no.range");
142 continue KMs; 199 continue;
143 } 200 }
144 201
145 double [][] wqs = wst.interpolateTabulated(currentKm); 202 double [][] wqs = wst.interpolateTabulated(currentKm);
146 if (wqs == null) { 203 if (wqs == null) {
147 // TODO: i18n 204 // TODO: i18n
153 if (!DoubleUtil.isValid(wqs)) { 210 if (!DoubleUtil.isValid(wqs)) {
154 // TODO: i18n 211 // TODO: i18n
155 addProblem(currentKm, "extreme.invalid.data"); 212 addProblem(currentKm, "extreme.invalid.data");
156 continue; 213 continue;
157 } 214 }
158 // TODO: Implement extraction of points for curve fitting. 215
159 // TODO: Implement curve fitting. 216 double [][] fitWQs = extractPointsToFit(wqs);
160 // TODO: Implement generating Curve object per km. 217 if (fitWQs == null) {
161 } 218 // TODO: i18n
162 219 addProblem(currentKm, "extreme.too.less.points");
163 ExtremeResult result = new ExtremeResult(); 220 continue;
221 }
222
223 double [] coeffs = doFitting(function, wqs, chiSqr);
224 if (coeffs == null) {
225 // TODO: i18n
226 addProblem(currentKm, "extreme.fitting.failed");
227 }
228
229 Curve curve = new Curve(
230 wqs[1], wqs[0],
231 function.getName(),
232 coeffs,
233 chiSqr[0]);
234
235 curves.add(currentKm, curve);
236
237 double [] values = range.getValues();
238
239 int V = Math.min(values.length, wqkms.length);
240 for (int i = 0; i < V; ++i) {
241 double q = values[i];
242 double w = curve.value(q);
243 if (Double.isNaN(w)) {
244 // TODO: i18n
245 addProblem(currentKm, "extreme.evaluate.failed", values[i]);
246 }
247 else {
248 wqkms[i].add(w, q);
249 }
250 }
251 }
252
253 ExtremeResult result = new ExtremeResult(curves, wqkms);
164 return new CalculationResult(result, this); 254 return new CalculationResult(result, this);
255 }
256
257 protected static double [] doFitting(
258 Function function,
259 double [][] wqs,
260 double [] chiSqr
261 ) {
262 LevenbergMarquardtOptimizer lmo = null;
263
264 double [] coeffs = null;
265
266 double [] ws = wqs[0];
267 double [] qs = wqs[1];
268
269 for (double tolerance = 1e-10; tolerance < 1e-3; tolerance *= 10d) {
270 lmo = new LevenbergMarquardtOptimizer();
271 lmo.setCostRelativeTolerance(tolerance);
272 lmo.setOrthoTolerance(tolerance);
273 lmo.setParRelativeTolerance(tolerance);
274
275 CurveFitter cf = new CurveFitter(lmo);
276
277 for (int i = ws.length-1; i >= 0; --i) {
278 cf.addObservedPoint(qs[i], ws[i]);
279 }
280
281 try {
282 coeffs = cf.fit(function, function.getInitialGuess());
283 break;
284 }
285 catch (MathException me) {
286 if (log.isDebugEnabled()) {
287 log.debug("tolerance " + tolerance + " + failed.");
288 }
289 }
290 }
291 if (coeffs == null) {
292 return null;
293 }
294 chiSqr[0] = lmo.getChiSquare();
295 return coeffs;
296 }
297
298 protected static double [][] extractPointsToFit(double [][] wqs) {
299 TDoubleArrayList ows = new TDoubleArrayList();
300 TDoubleArrayList oqs = new TDoubleArrayList();
301
302 double [] ws = wqs[0];
303 double [] qs = wqs[1];
304
305 int N = Math.min(ws.length, qs.length);
306
307 if (N < 2) {
308 log.warn("Too less points for fitting");
309 return null;
310 }
311
312 double q2 = qs[N-1];
313 double w2 = ws[N-1];
314 double q1 = qs[N-2];
315 double w1 = ws[N-2];
316
317 oqs.add(q2); oqs.add(q1);
318 ows.add(w2); ows.add(w1);
319
320 int orientation = Utils.epsilonEquals(w1, w1)
321 ? 0
322 : w1 > w2
323 ? -1
324 : +1;
325
326 for (int i = N-3; i >= 0; --i) {
327 double cq = qs[i];
328 double cw = qs[i];
329 int newOrientation = Utils.relativeCCW(
330 q2, w2,
331 q1, w1,
332 cq, cw);
333
334 if (newOrientation != 0 && newOrientation != orientation) {
335 break;
336 }
337 oqs.add(cq);
338 ows.add(cw);
339
340 if (newOrientation != 0) {
341 // rotate last out
342 q2 = q1;
343 w2 = w1;
344 }
345 q1 = cq;
346 w1 = cw;
347 }
348
349 // XXX: Not really needed for fitting.
350 // oqs.reverse();
351 // ows.reverse();
352
353 return new double [][] {
354 ows.toNativeArray(),
355 oqs.toNativeArray()
356 };
165 } 357 }
166 } 358 }
167 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf-8 : 359 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf-8 :
http://dive4elements.wald.intevation.org