Mercurial > dive4elements > river

comparison flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java @ 1190:f514894ec2fd

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merged flys-artifacts/2.5
author Thomas Arendsen Hein <thomas@intevation.de>
date Fri, 28 Sep 2012 12:14:17 +0200
parents 1ea7eb72aaa6
children e5f7f25a511c
comparison
equal deleted inserted replaced
917:b48c36076e17 1190:f514894ec2fd
1 package de.intevation.flys.artifacts.model;
2
3 import java.io.Serializable;
4
5 import de.intevation.flys.artifacts.math.Linear;
6 import de.intevation.flys.artifacts.math.Function;
7
8 import java.util.Arrays;
9 import java.util.ArrayList;
10 import java.util.List;
11 import java.util.Collections;
12
13 import org.apache.log4j.Logger;
14
15 import org.apache.commons.math.analysis.interpolation.SplineInterpolator;
16
17 import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
18
19 import org.apache.commons.math.ArgumentOutsideDomainException;
20
21 import org.apache.commons.math.exception.MathIllegalArgumentException;
22
23 public class WstValueTable
24 implements Serializable
25 {
26 private static Logger log = Logger.getLogger(WstValueTable.class);
27
28 public static final int DEFAULT_Q_STEPS = 500;
29
30 public static final class Column
31 implements Serializable
32 {
33 protected String name;
34
35 protected QRangeTree qRangeTree;
36
37 public Column() {
38 }
39
40 public Column(String name) {
41 this.name = name;
42 }
43
44 public String getName() {
45 return name;
46 }
47
48 public void setName(String name) {
49 this.name = name;
50 }
51
52 public QRangeTree getQRangeTree() {
53 return qRangeTree;
54 }
55
56 public void setQRangeTree(QRangeTree qRangeTree) {
57 this.qRangeTree = qRangeTree;
58 }
59 } // class Column
60
61 public static final class QPosition {
62
63 protected int index;
64 protected double weight;
65
66 public QPosition() {
67 }
68
69 public QPosition(int index, double weight) {
70 this.index = index;
71 this.weight = weight;
72 }
73
74 public QPosition set(int index, double weight) {
75 this.index = index;
76 this.weight = weight;
77 return this;
78 }
79
80 } // class Position
81
82 public static final class Row
83 implements Serializable, Comparable<Row>
84 {
85 double km;
86 double [] ws;
87
88 public Row() {
89 }
90
91 public Row(double km) {
92 this.km = km;
93 }
94
95 public Row(double km, double [] ws) {
96 this(km);
97 this.ws = ws;
98 }
99
100 public int compareTo(Row other) {
101 double d = km - other.km;
102 if (d < -0.0001) return -1;
103 if (d > 0.0001) return +1;
104 return 0;
105 }
106
107 public void interpolateW(
108 Row other,
109 double km,
110 double [] iqs,
111 double [] ows,
112 WstValueTable table,
113 Calculation errors
114 ) {
115 double kmWeight = Linear.factor(km, this.km, other.km);
116
117 QPosition qPosition = new QPosition();
118
119 for (int i = 0; i < iqs.length; ++i) {
120 if (table.getQPosition(km, iqs[i], qPosition) != null) {
121 double wt = getW(qPosition);
122 double wo = other.getW(qPosition);
123 if (Double.isNaN(wt) || Double.isNaN(wo)) {
124 if (errors != null) {
125 // TODO: I18N
126 errors.addProblem(
127 km, "cannot find w for q = " + iqs[i]);
128 }
129 ows[i] = Double.NaN;
130 }
131 else {
132 ows[i] = Linear.weight(kmWeight, wt, wo);
133 }
134 }
135 else {
136 if (errors != null) {
137 // TODO: I18N
138 errors.addProblem(km, "cannot find q = " + iqs[i]);
139 }
140 ows[i] = Double.NaN;
141 }
142 }
143 }
144
145 public double [][] interpolateWQ(
146 Row other,
147 double km,
148 int steps,
149 WstValueTable table,
150 Calculation errors
151 ) {
152 int W = Math.min(ws.length, other.ws.length);
153
154 if (W < 1) {
155 if (errors != null) {
156 // TODO: I18N
157 errors.addProblem("no ws found");
158 }
159 return new double[2][0];
160 }
161
162 double factor = Linear.factor(km, this.km, other.km);
163
164 double [] splineQ = new double[W];
165 double [] splineW = new double[W];
166
167 double minQ = Double.MAX_VALUE;
168 double maxQ = -Double.MAX_VALUE;
169
170 for (int i = 0; i < W; ++i) {
171 double wws = Linear.weight(factor, ws[i], other.ws[i]);
172 double wqs = Linear.weight(
173 factor,
174 table.getQIndex(i, km),
175 table.getQIndex(i, other.km));
176
177 if (Double.isNaN(wws) || Double.isNaN(wqs)) {
178 if (errors != null) {
179 // TODO: I18N
180 errors.addProblem(km, "cannot find w or q");
181 }
182 }
183 else {
184 if (wqs < minQ) minQ = wqs;
185 if (wqs > maxQ) maxQ = wqs;
186 }
187
188 splineW[i] = wws;
189 splineQ[i] = wqs;
190 }
191
192 double stepWidth = (maxQ - minQ)/steps;
193
194 SplineInterpolator interpolator = new SplineInterpolator();
195 PolynomialSplineFunction spline;
196
197 try {
198 spline = interpolator.interpolate(splineQ, splineW);
199 }
200 catch (MathIllegalArgumentException miae) {
201 if (errors != null) {
202 // TODO: I18N
203 errors.addProblem(km, "spline creation failed");
204 }
205 log.error("spline creation failed");
206 return new double[2][0];
207 }
208
209 double [] outWs = new double[steps];
210 double [] outQs = new double[steps];
211
212 Arrays.fill(outWs, Double.NaN);
213 Arrays.fill(outQs, Double.NaN);
214
215 try {
216 double q = minQ;
217 for (int i = 0; i < outWs.length; ++i, q += stepWidth) {
218 outWs[i] = spline.value(outQs[i] = q);
219 }
220 }
221 catch (ArgumentOutsideDomainException aode) {
222 if (errors != null) {
223 // TODO: I18N
224 errors.addProblem(km, "spline interpolation failed");
225 }
226 log.error("spline interpolation failed", aode);
227 }
228
229 return new double [][] { outWs, outQs };
230 }
231
232 public double [][] interpolateWQ(
233 int steps,
234 WstValueTable table,
235 Calculation errors
236 ) {
237 int W = ws.length;
238
239 if (W < 1) {
240 if (errors != null) {
241 // TODO: I18N
242 errors.addProblem(km, "no ws found");
243 }
244 return new double[2][0];
245 }
246
247 double [] splineQ = new double[W];
248
249 double minQ = Double.MAX_VALUE;
250 double maxQ = -Double.MAX_VALUE;
251
252 for (int i = 0; i < W; ++i) {
253 double sq = table.getQIndex(i, km);
254 if (Double.isNaN(sq)) {
255 if (errors != null) {
256 // TODO: I18N
257 errors.addProblem(
258 km, "no q found in " + (i+1) + " column");
259 }
260 }
261 else {
262 if (sq < minQ) minQ = sq;
263 if (sq > maxQ) maxQ = sq;
264 }
265 splineQ[i] = sq;
266 }
267
268 double stepWidth = (maxQ - minQ)/steps;
269
270 SplineInterpolator interpolator = new SplineInterpolator();
271
272 PolynomialSplineFunction spline;
273
274 try {
275 spline = interpolator.interpolate(splineQ, ws);
276 }
277 catch (MathIllegalArgumentException miae) {
278 if (errors != null) {
279 // TODO: I18N
280 errors.addProblem(km, "spline creation failed");
281 }
282 log.error("spline creation failed");
283 return new double[2][0];
284 }
285
286 double [] outWs = new double[steps];
287 double [] outQs = new double[steps];
288
289 Arrays.fill(outWs, Double.NaN);
290 Arrays.fill(outQs, Double.NaN);
291
292 try {
293 double q = minQ;
294 for (int i = 0; i < outWs.length; ++i, q += stepWidth) {
295 outWs[i] = spline.value(outQs[i] = q);
296 }
297 }
298 catch (ArgumentOutsideDomainException aode) {
299 if (errors != null) {
300 // TODO: I18N
301 errors.addProblem(km, "spline interpolation failed");
302 }
303 log.error("spline interpolation failed.", aode);
304 }
305
306 return new double [][] { outWs, outQs };
307 }
308
309
310 public double getW(QPosition qPosition) {
311 int index = qPosition.index;
312 double weight = qPosition.weight;
313
314 return weight == 1.0
315 ? ws[index]
316 : Linear.weight(weight, ws[index-1], ws[index]);
317 }
318
319 public double getW(
320 Row other,
321 double km,
322 QPosition qPosition
323 ) {
324 double kmWeight = Linear.factor(km, this.km, other.km);
325
326 int index = qPosition.index;
327 double weight = qPosition.weight;
328
329 double tw, ow;
330
331 if (weight == 1.0) {
332 tw = ws[index];
333 ow = other.ws[index];
334 }
335 else {
336 tw = Linear.weight(weight, ws[index-1], ws[index]);
337 ow = Linear.weight(weight, other.ws[index-1], other.ws[index]);
338 }
339
340 return Linear.weight(kmWeight, tw, ow);
341 }
342 } // class Row
343
344 protected List<Row> rows;
345
346 protected Column [] columns;
347
348 public WstValueTable() {
349 rows = new ArrayList<Row>();
350 }
351
352 public WstValueTable(Column [] columns) {
353 this();
354 this.columns = columns;
355 }
356
357 public WstValueTable(Column [] columns, List<Row> rows) {
358 this.columns = columns;
359 this.rows = rows;
360 }
361
362 public void sortRows() {
363 Collections.sort(rows);
364 }
365
366 public double [] interpolateW(double km, double [] qs, double [] ws) {
367 return interpolateW(km, qs, ws, null);
368 }
369
370
371 /**
372 * @param ws (output parameter), gets returned.
373 * @return output parameter ws.
374 */
375 public double [] interpolateW(
376 double km,
377 double [] qs,
378 double [] ws,
379 Calculation errors
380 ) {
381 int rowIndex = Collections.binarySearch(rows, new Row(km));
382
383 QPosition qPosition = new QPosition();
384
385 if (rowIndex >= 0) { // direct row match
386 Row row = rows.get(rowIndex);
387 for (int i = 0; i < qs.length; ++i) {
388 if (getQPosition(km, qs[i], qPosition) == null) {
389 if (errors != null) {
390 // TODO: I18N
391 errors.addProblem(km, "cannot find q = " + qs[i]);
392 }
393 ws[i] = Double.NaN;
394 }
395 else {
396 if (Double.isNaN(ws[i] = row.getW(qPosition))
397 && errors != null) {
398 // TODO: I18N
399 errors.addProblem(
400 km, "cannot find w for q = " + qs[i]);
401 }
402 }
403 }
404 }
405 else { // needs bilinear interpolation
406 rowIndex = -rowIndex -1;
407
408 if (rowIndex < 1 || rowIndex >= rows.size()) {
409 // do not extrapolate
410 Arrays.fill(ws, Double.NaN);
411 if (errors != null) {
412 // TODO: I18N
413 errors.addProblem(km, "km not found");
414 }
415 }
416 else {
417 Row r1 = rows.get(rowIndex-1);
418 Row r2 = rows.get(rowIndex);
419 r1.interpolateW(r2, km, qs, ws, this, errors);
420 }
421 }
422
423 return ws;
424 }
425
426 public double [][] interpolateWQ(double km) {
427 return interpolateWQ(km, null);
428 }
429
430 public double [][] interpolateWQ(double km, Calculation errors) {
431 return interpolateWQ(km, DEFAULT_Q_STEPS, errors);
432 }
433
434 public double [][] interpolateWQ(double km, int steps, Calculation errors) {
435
436 int rowIndex = Collections.binarySearch(rows, new Row(km));
437
438 if (rowIndex >= 0) { // direct row match
439 Row row = rows.get(rowIndex);
440 return row.interpolateWQ(steps, this, errors);
441 }
442
443 rowIndex = -rowIndex -1;
444
445 if (rowIndex < 1 || rowIndex >= rows.size()) {
446 // do not extrapolate
447 if (errors != null) {
448 // TODO: I18N
449 errors.addProblem(km, "km not found");
450 }
451 return new double[2][0];
452 }
453
454 Row r1 = rows.get(rowIndex-1);
455 Row r2 = rows.get(rowIndex);
456
457 return r1.interpolateWQ(r2, km, steps, this, errors);
458 }
459
460 public boolean interpolate(
461 double km,
462 double [] out,
463 QPosition qPosition,
464 Function qFunction
465 ) {
466 int R1 = rows.size()-1;
467
468 out[1] = qFunction.value(getQ(qPosition, km));
469
470 if (Double.isNaN(out[1])) {
471 return false;
472 }
473
474 QPosition nPosition = new QPosition();
475 if (getQPosition(km, out[1], nPosition) == null) {
476 return false;
477 }
478
479 int rowIndex = Collections.binarySearch(rows, new Row(km));
480
481 if (rowIndex >= 0) {
482 // direct row match
483 out[0] = rows.get(rowIndex).getW(nPosition);
484 return !Double.isNaN(out[0]);
485 }
486
487 rowIndex = -rowIndex -1;
488
489 if (rowIndex < 1 || rowIndex >= R1) {
490 // do not extrapolate
491 return false;
492 }
493
494 Row r1 = rows.get(rowIndex-1);
495 Row r2 = rows.get(rowIndex);
496 out[0] = r1.getW(r2, km, nPosition);
497
498 return !Double.isNaN(out[0]);
499 }
500
501
502 /**
503 * Look up interpolation of a Q at given positions.
504 *
505 * @param q the non-interpolated Q value.
506 * @param referenceKm the reference km (e.g. gauge position).
507 * @param kms positions for which to interpolate.
508 * @param ws (output) resulting interpolated ws.
509 * @param qs (output) resulting interpolated qs.
510 * @param errors calculation object to store errors.
511 */
512 public QPosition interpolate(
513 double q,
514 double referenceKm,
515 double [] kms,
516 double [] ws,
517 double [] qs,
518 Calculation errors
519 ) {
520 return interpolate(
521 q, referenceKm, kms, ws, qs, 0, kms.length, errors);
522 }
523
524 public QPosition interpolate(
525 double q,
526 double referenceKm,
527 double [] kms,
528 double [] ws,
529 double [] qs,
530 int startIndex,
531 int length,
532 Calculation errors
533 ) {
534 QPosition qPosition = getQPosition(referenceKm, q);
535
536 if (qPosition == null) {
537 // we cannot locate q at km
538 Arrays.fill(ws, Double.NaN);
539 Arrays.fill(qs, Double.NaN);
540 if (errors != null) {
541 errors.addProblem(referenceKm, "cannot find q");
542 }
543 return null;
544 }
545
546 Row kmKey = new Row();
547
548 int R1 = rows.size()-1;
549
550 for (int i = startIndex, end = startIndex+length; i < end; ++i) {
551
552 if (Double.isNaN(qs[i] = getQ(qPosition, kms[i]))) {
553 if (errors != null) {
554 errors.addProblem(kms[i], "cannot find q");
555 }
556 ws[i] = Double.NaN;
557 continue;
558 }
559
560 kmKey.km = kms[i];
561 int rowIndex = Collections.binarySearch(rows, kmKey);
562
563 if (rowIndex >= 0) {
564 // direct row match
565 if (Double.isNaN(ws[i] = rows.get(rowIndex).getW(qPosition))
566 && errors != null) {
567 errors.addProblem(kms[i], "cannot find w");
568 }
569 continue;
570 }
571
572 rowIndex = -rowIndex -1;
573
574 if (rowIndex < 1 || rowIndex >= R1) {
575 // do not extrapolate
576 if (errors != null) {
577 errors.addProblem(kms[i], "cannot find km");
578 }
579 ws[i] = Double.NaN;
580 continue;
581 }
582 Row r1 = rows.get(rowIndex-1);
583 Row r2 = rows.get(rowIndex);
584
585 if (Double.isNaN(ws[i] = r1.getW(r2, kms[i], qPosition))
586 && errors != null) {
587 errors.addProblem(kms[i], "cannot find w");
588 }
589 }
590
591 return qPosition;
592 }
593
594 public QPosition getQPosition(double km, double q) {
595 return getQPosition(km, q, new QPosition());
596 }
597
598 public QPosition getQPosition(double km, double q, QPosition qPosition) {
599
600 if (columns.length == 0) {
601 return null;
602 }
603
604 double qLast = columns[0].getQRangeTree().findQ(km);
605
606 if (Math.abs(qLast - q) < 0.00001) {
607 return qPosition.set(0, 1d);
608 }
609
610 for (int i = 1; i < columns.length; ++i) {
611 double qCurrent = columns[i].getQRangeTree().findQ(km);
612 if (Math.abs(qCurrent - q) < 0.00001) {
613 return qPosition.set(i, 1d);
614 }
615
616 double qMin, qMax;
617 if (qLast < qCurrent) { qMin = qLast; qMax = qCurrent; }
618 else { qMin = qCurrent; qMax = qLast; }
619
620 if (q > qMin && q < qMax) {
621 double weight = Linear.factor(q, qLast, qCurrent);
622 return qPosition.set(i, weight);
623 }
624 qLast = qCurrent;
625 }
626
627 return null;
628 }
629
630 public double getQIndex(int index, double km) {
631 return columns[index].getQRangeTree().findQ(km);
632 }
633
634 public double getQ(QPosition qPosition, double km) {
635 int index = qPosition.index;
636 double weight = qPosition.weight;
637
638 if (weight == 1d) {
639 return columns[index].getQRangeTree().findQ(km);
640 }
641 double q1 = columns[index-1].getQRangeTree().findQ(km);
642 double q2 = columns[index ].getQRangeTree().findQ(km);
643 return Linear.weight(weight, q1, q2);
644 }
645
646 }
647 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
http://dive4elements.wald.intevation.org