Mercurial > dive4elements > river

comparison flys-artifacts/src/main/java/org/dive4elements/river/artifacts/model/extreme/ExtremeCalculation.java @ 5831:bd047b71ab37

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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/model/extreme/ExtremeCalculation.java@5b8919ef601d
children
comparison
equal deleted inserted replaced
5830:160f53ee0870 5831:bd047b71ab37
1 package org.dive4elements.river.artifacts.model.extreme;
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
12 import org.dive4elements.river.artifacts.access.ExtremeAccess;
13
14 import org.dive4elements.river.artifacts.math.Linear;
15 //import org.dive4elements.river.artifacts.math.Utils;
16
17 import org.dive4elements.river.artifacts.math.fitting.Function;
18 import org.dive4elements.river.artifacts.math.fitting.FunctionFactory;
19
20 import org.dive4elements.river.artifacts.model.Calculation;
21 import org.dive4elements.river.artifacts.model.CalculationResult;
22 import org.dive4elements.river.artifacts.model.RangeWithValues;
23 import org.dive4elements.river.artifacts.model.RiverFactory;
24 import org.dive4elements.river.artifacts.model.WQKms;
25 import org.dive4elements.river.artifacts.model.WstValueTable;
26 import org.dive4elements.river.artifacts.model.WstValueTableFactory;
27
28 import org.dive4elements.river.model.River;
29
30 import org.dive4elements.river.utils.DoubleUtil;
31 import org.dive4elements.river.utils.KMIndex;
32
33 import gnu.trove.TDoubleArrayList;
34
35 import java.util.List;
36
37 import java.awt.geom.Line2D;
38
39 /** Calculate extrapolated W. */
40 public class ExtremeCalculation
41 extends Calculation
42 {
43 private static final Log log =
44 LogFactory.getLog(ExtremeCalculation.class);
45
46 protected String river;
47 protected String function;
48 protected double from;
49 protected double to;
50 protected double step;
51 protected double percent;
52 protected List<RangeWithValues> ranges;
53
54 public ExtremeCalculation() {
55 }
56
57 public ExtremeCalculation(ExtremeAccess access) {
58 String river = access.getRiver();
59 String function = access.getFunction();
60 Double from = access.getFrom();
61 Double to = access.getTo();
62 Double step = access.getStep();
63 Double percent = access.getPercent();
64 List<RangeWithValues> ranges = access.getRanges();
65
66 if (river == null) {
67 // TODO: i18n
68 addProblem("extreme.no.river");
69 }
70
71 if (function == null) {
72 // TODO: i18n
73 addProblem("extreme.no.function");
74 }
75
76 if (from == null) {
77 // TODO: i18n
78 addProblem("extreme.no.from");
79 }
80
81 if (to == null) {
82 // TODO: i18n
83 addProblem("extreme.no.to");
84 }
85
86 if (step == null) {
87 // TODO: i18n
88 addProblem("extreme.no.step");
89 }
90
91 if (percent == null) {
92 // TODO: i18n
93 addProblem("extreme.no.percent");
94 }
95
96 if (ranges == null) {
97 // TODO: i18n
98 addProblem("extreme.no.ranges");
99 }
100
101 if (!hasProblems()) {
102 this.river = river;
103 this.function = function;
104 this.from = Math.min(from, to);
105 this.to = Math.max(from, to);
106 this.step = Math.max(0.001d, Math.abs(step)/1000d);
107 this.percent = Math.max(0d, Math.min(100d, percent));
108 this.ranges = ranges;
109 }
110 }
111
112
113 /** Calculate an extreme curve (extrapolate). */
114 public CalculationResult calculate() {
115
116 WstValueTable wst = null;
117
118 River river = RiverFactory.getRiver(this.river);
119 if (river == null) {
120 // TODO: i18n
121 addProblem("extreme.no.such.river", this.river);
122 }
123 else {
124 wst = WstValueTableFactory.getTable(river);
125 if (wst == null) {
126 // TODO: i18n
127 addProblem("extreme.no.wst.table");
128 }
129 }
130
131 Function function =
132 FunctionFactory.getInstance().getFunction(this.function);
133 if (function == null) {
134 // TODO: i18n
135 addProblem("extreme.no.such.function", this.function);
136 }
137
138 return hasProblems()
139 ? new CalculationResult(this)
140 : innerCalculate(wst, function);
141 }
142
143
144 /** Name of wqkms like W(5000,6000) */
145 protected String wqkmsName(int i) {
146 StringBuilder sb = new StringBuilder("W(");
147 boolean already = false;
148 for (RangeWithValues r: ranges) {
149 double [] values = r.getValues();
150 if (i < values.length) {
151 if (already) {
152 sb.append(", ");
153 }
154 else {
155 already = true;
156 }
157 // TODO: i18n
158 sb.append(values[i]);
159 }
160 }
161 return sb.append(')').toString();
162 }
163
164 protected WQKms [] allocWQKms() {
165 int max = 0;
166 for (RangeWithValues r: ranges) {
167 double [] values = r.getValues();
168 if (values.length > max) {
169 max = values.length;
170 }
171 }
172 WQKms [] wqkms = new WQKms[max];
173 for (int i = 0; i < max; ++i) {
174 wqkms[i] = new WQKms(wqkmsName(i));
175 }
176 return wqkms;
177 }
178
179
180 /** Calculate an extreme curve (extrapolate). */
181 protected CalculationResult innerCalculate(
182 WstValueTable wst,
183 Function function
184 ) {
185 RangeWithValues range = null;
186
187 double [] chiSqr = { 0d };
188
189 KMIndex<Curve> curves = new KMIndex<Curve>();
190 WQKms [] wqkms = allocWQKms();
191
192 boolean debug = log.isDebugEnabled();
193
194 from = DoubleUtil.round(from);
195 to = DoubleUtil.round(to);
196
197 for (double km = from; km <= to; km = DoubleUtil.round(km+step)) {
198
199 if (debug) {
200 log.debug("km: " + km);
201 }
202
203 boolean foundRange = false;
204
205 if (range == null || !range.inside(km)) {
206 for (RangeWithValues r: ranges) {
207 if (r.inside(km)) {
208 range = r;
209 foundRange = true;
210 break;
211 }
212 }
213 // TODO: i18n
214 if (!foundRange) {
215 addProblem(km, "extreme.no.range.inner");
216 continue;
217 }
218 }
219
220 double [][] wqs = wst.interpolateTabulated(km);
221 if (wqs == null) {
222 // TODO: i18n
223 addProblem(km, "extreme.no.raw.data");
224 continue;
225 }
226
227 // XXX: This should not be necessary for model data.
228 if (!DoubleUtil.isValid(wqs)) {
229 // TODO: i18n
230 addProblem(km, "extreme.invalid.data");
231 continue;
232 }
233
234 double [][] fitWQs = extractPointsToFit(wqs);
235 if (fitWQs == null) {
236 // TODO: i18n
237 addProblem(km, "extreme.too.less.points");
238 continue;
239 }
240
241 double [] coeffs = doFitting(function, fitWQs, chiSqr);
242 if (coeffs == null) {
243 // TODO: i18n
244 addProblem(km, "extreme.fitting.failed");
245 continue;
246 }
247
248 Curve curve = new Curve(
249 wqs[1], wqs[0],
250 function.getName(),
251 coeffs,
252 chiSqr[0]);
253
254 curves.add(km, curve);
255
256 double [] values = range.getValues();
257
258 int V = Math.min(values.length, wqkms.length);
259 for (int i = 0; i < V; ++i) {
260 double q = values[i];
261 double w = curve.value(q);
262 if (Double.isNaN(w)) {
263 // TODO: i18n
264 addProblem(km, "extreme.evaluate.failed", values[i]);
265 }
266 else {
267 wqkms[i].add(w, q, km);
268 }
269 }
270 }
271
272 ExtremeResult result = new ExtremeResult(curves, wqkms);
273 return new CalculationResult(result, this);
274 }
275
276 protected double [] doFitting(
277 Function function,
278 double [][] wqs,
279 double [] chiSqr
280 ) {
281 LevenbergMarquardtOptimizer lmo = null;
282
283 double [] coeffs = null;
284
285 double [] ws = wqs[0];
286 double [] qs = wqs[1];
287
288 for (double tolerance = 1e-10; tolerance < 1e-3; tolerance *= 10d) {
289 lmo = new LevenbergMarquardtOptimizer();
290 lmo.setCostRelativeTolerance(tolerance);
291 lmo.setOrthoTolerance(tolerance);
292 lmo.setParRelativeTolerance(tolerance);
293
294 CurveFitter cf = new CurveFitter(lmo);
295
296 for (int i = 0; i < ws.length; ++i) {
297 cf.addObservedPoint(qs[i], ws[i]);
298 }
299
300 try {
301 coeffs = cf.fit(function, function.getInitialGuess());
302 break;
303 }
304 catch (MathException me) {
305 if (log.isDebugEnabled()) {
306 log.debug("tolerance " + tolerance + " + failed.");
307 }
308 }
309 }
310 if (coeffs != null) {
311 chiSqr[0] = lmo.getChiSquare();
312 }
313 return coeffs;
314 }
315
316 protected double [][] extractPointsToFit(double [][] wqs) {
317
318 double [] ws = wqs[0];
319 double [] qs = wqs[1];
320
321 int N = Math.min(ws.length, qs.length);
322
323 if (N < 2) {
324 log.warn("Too less points for fitting");
325 return null;
326 }
327
328 double q2 = qs[N-1];
329 double w2 = ws[N-1];
330 double q1 = qs[N-2];
331 double w1 = ws[N-2];
332
333 boolean ascending = w2 > w1;
334
335 TDoubleArrayList ows = new TDoubleArrayList();
336 TDoubleArrayList oqs = new TDoubleArrayList();
337
338 oqs.add(q2); oqs.add(q1);
339 ows.add(w2); ows.add(w1);
340
341 int lastDir = -2;
342
343 for (int i = N-3; i >= 0; --i) {
344 double q = qs[i];
345 double w = ws[i];
346
347 if ((ascending && w > w1) || (!ascending && w < w1)) {
348 break;
349 }
350
351 int dir = Line2D.relativeCCW(q2, w2, q1, w1, q, w);
352 //int dir = Utils.relativeCCW(q2, w2, q1, w1, q, w);
353 if (lastDir == -2) {
354 lastDir = dir;
355 }
356 else if (lastDir != dir) {
357 break;
358 }
359
360 oqs.add(q);
361 ows.add(w);
362 w2 = w1;
363 q2 = q1;
364 w1 = w;
365 q1 = q;
366 }
367
368 oqs.reverse();
369 ows.reverse();
370
371 boolean debug = log.isDebugEnabled();
372 if (debug) {
373 log.debug("from table: " + N);
374 log.debug("after trim: " + oqs.size());
375 }
376
377 cutPercent(ows, oqs);
378
379 if (debug) {
380 log.debug("after percent cut: " + oqs.size());
381 }
382
383 return new double [][] {
384 ows.toNativeArray(),
385 oqs.toNativeArray()
386 };
387 }
388
389
390 protected void cutPercent(TDoubleArrayList ws, TDoubleArrayList qs) {
391 int N = qs.size();
392 if (percent <= 0d || N == 0) {
393 return;
394 }
395
396 double minQ = qs.getQuick(0);
397 double maxQ = qs.getQuick(N-1);
398 double factor = Math.min(Math.max(0d, percent/100d), 1d);
399 double cutQ = Linear.weight(factor, minQ, maxQ);
400 int cutIndex = 0;
401 for (; cutIndex < N; ++cutIndex) {
402 double q = qs.getQuick(cutIndex);
403 if (minQ < maxQ) {
404 if (q > cutQ) {
405 break;
406 }
407 }
408 else {
409 if (q < cutQ) {
410 break;
411 }
412 }
413 }
414 ws.remove(0, cutIndex);
415 qs.remove(0, cutIndex);
416 }
417 }
418 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf-8 :
http://dive4elements.wald.intevation.org