Mercurial > dive4elements > river

comparison flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java @ 3468:f37e7e8907cb

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merged flys-artifacts/2.8.1
author Thomas Arendsen Hein <thomas@intevation.de>
date Fri, 28 Sep 2012 12:14:39 +0200
parents 4c4ec9e9650a
children da1969b05292
comparison
equal deleted inserted replaced
3387:5ffad8bde8ad 3468:f37e7e8907cb
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 import gnu.trove.TDoubleArrayList;
24
25 /**
26 * W, Q and km data from database 'wsts' spiced with interpolation algorithms.
27 */
28 public class WstValueTable
29 implements Serializable
30 {
31 private static Logger log = Logger.getLogger(WstValueTable.class);
32
33 public static final int DEFAULT_Q_STEPS = 500;
34
35 public static final int RELATE_WS_SAMPLES = 200;
36
37 /**
38 * A Column in the table, typically representing one measurement session.
39 */
40 public static final class Column
41 implements Serializable
42 {
43 protected String name;
44
45 protected QRangeTree qRangeTree;
46
47 public Column() {
48 }
49
50 public Column(String name) {
51 this.name = name;
52 }
53
54 public String getName() {
55 return name;
56 }
57
58 public void setName(String name) {
59 this.name = name;
60 }
61
62 public QRangeTree getQRangeTree() {
63 return qRangeTree;
64 }
65
66 public void setQRangeTree(QRangeTree qRangeTree) {
67 this.qRangeTree = qRangeTree;
68 }
69 } // class Column
70
71 /**
72 * A (weighted) position used for interpolation.
73 */
74 public static final class QPosition {
75
76 protected int index;
77 protected double weight;
78
79 public QPosition() {
80 }
81
82 public QPosition(int index, double weight) {
83 this.index = index;
84 this.weight = weight;
85 }
86
87 public QPosition set(int index, double weight) {
88 this.index = index;
89 this.weight = weight;
90 return this;
91 }
92
93 } // class Position
94
95 public static final class SplineFunction {
96
97 public PolynomialSplineFunction spline;
98 public double [] splineQs;
99 public double [] splineWs;
100
101 public SplineFunction(
102 PolynomialSplineFunction spline,
103 double [] splineQs,
104 double [] splineWs
105 ) {
106 this.spline = spline;
107 this.splineQs = splineQs;
108 this.splineWs = splineWs;
109 }
110
111 public double [][] sample(
112 int numSamples,
113 double km,
114 Calculation errors
115 ) {
116 double minQ = getQMin();
117 double maxQ = getQMax();
118
119 double [] outWs = new double[numSamples];
120 double [] outQs = new double[numSamples];
121
122 Arrays.fill(outWs, Double.NaN);
123 Arrays.fill(outQs, Double.NaN);
124
125 double stepWidth = (maxQ - minQ)/numSamples;
126
127 try {
128 double q = minQ;
129 for (int i = 0; i < outWs.length; ++i, q += stepWidth) {
130 outWs[i] = spline.value(outQs[i] = q);
131 }
132 }
133 catch (ArgumentOutsideDomainException aode) {
134 if (errors != null) {
135 errors.addProblem(km, "spline.interpolation.failed");
136 }
137 log.error("spline interpolation failed.", aode);
138 }
139
140 return new double [][] { outWs, outQs };
141 }
142
143 public double getQMin() {
144 return Math.min(splineQs[0], splineQs[splineQs.length-1]);
145 }
146
147 public double getQMax() {
148 return Math.max(splineQs[0], splineQs[splineQs.length-1]);
149 }
150
151 /** Constructs a continues index between the columns to Qs. */
152 public PolynomialSplineFunction createIndexQRelation() {
153
154 double [] indices = new double[splineQs.length];
155 for (int i = 0; i < indices.length; ++i) {
156 indices[i] = i;
157 }
158
159 try {
160 SplineInterpolator interpolator = new SplineInterpolator();
161 return interpolator.interpolate(indices, splineQs);
162 }
163 catch (MathIllegalArgumentException miae) {
164 // Ignore me!
165 }
166 return null;
167 }
168 } // class SplineFunction
169
170 /**
171 * A row, typically a position where measurements were taken.
172 */
173 public static final class Row
174 implements Serializable, Comparable<Row>
175 {
176 double km;
177 double [] ws;
178
179 public Row() {
180 }
181
182 public Row(double km) {
183 this.km = km;
184 }
185
186 public Row(double km, double [] ws) {
187 this(km);
188 this.ws = ws;
189 }
190
191 /**
192 * Compare according to place of measurement (km).
193 */
194 public int compareTo(Row other) {
195 double d = km - other.km;
196 if (d < -0.0001) return -1;
197 if (d > 0.0001) return +1;
198 return 0;
199 }
200
201 /**
202 * Interpolate Ws, given Qs and a km.
203 *
204 * @param iqs Given ("input") Qs.
205 * @param ows Resulting ("output") Ws.
206 * @param table Table of which to use data for interpolation.
207 */
208 public void interpolateW(
209 Row other,
210 double km,
211 double [] iqs,
212 double [] ows,
213 WstValueTable table,
214 Calculation errors
215 ) {
216 double kmWeight = Linear.factor(km, this.km, other.km);
217
218 QPosition qPosition = new QPosition();
219
220 for (int i = 0; i < iqs.length; ++i) {
221 if (table.getQPosition(km, iqs[i], qPosition) != null) {
222 double wt = getW(qPosition);
223 double wo = other.getW(qPosition);
224 if (Double.isNaN(wt) || Double.isNaN(wo)) {
225 if (errors != null) {
226 errors.addProblem(
227 km, "cannot.find.w.for.q", iqs[i]);
228 }
229 ows[i] = Double.NaN;
230 }
231 else {
232 ows[i] = Linear.weight(kmWeight, wt, wo);
233 }
234 }
235 else {
236 if (errors != null) {
237 errors.addProblem(km, "cannot.find.q", iqs[i]);
238 }
239 ows[i] = Double.NaN;
240 }
241 }
242 }
243
244
245 public SplineFunction createSpline(
246 WstValueTable table,
247 Calculation errors
248 ) {
249 int W = ws.length;
250
251 if (W < 1) {
252 if (errors != null) {
253 errors.addProblem(km, "no.ws.found");
254 }
255 return null;
256 }
257
258 double [] splineQs = new double[W];
259
260 for (int i = 0; i < W; ++i) {
261 double sq = table.getQIndex(i, km);
262 if (Double.isNaN(sq) && errors != null) {
263 errors.addProblem(
264 km, "no.q.found.in.column", (i+1));
265 }
266 splineQs[i] = sq;
267 }
268
269 try {
270 SplineInterpolator interpolator = new SplineInterpolator();
271 PolynomialSplineFunction spline =
272 interpolator.interpolate(splineQs, ws);
273
274 return new SplineFunction(spline, splineQs, ws);
275 }
276 catch (MathIllegalArgumentException miae) {
277 if (errors != null) {
278 errors.addProblem(km, "spline.creation.failed");
279 }
280 log.error("spline creation failed", miae);
281 }
282 return null;
283 }
284
285 public SplineFunction createSpline(
286 Row other,
287 double km,
288 WstValueTable table,
289 Calculation errors
290 ) {
291 int W = Math.min(ws.length, other.ws.length);
292
293 if (W < 1) {
294 if (errors != null) {
295 errors.addProblem("no.ws.found");
296 }
297 return null;
298 }
299
300 double factor = Linear.factor(km, this.km, other.km);
301
302 double [] splineQs = new double[W];
303 double [] splineWs = new double[W];
304
305 for (int i = 0; i < W; ++i) {
306 double wws = Linear.weight(factor, ws[i], other.ws[i]);
307 double wqs = Linear.weight(
308 factor,
309 table.getQIndex(i, km),
310 table.getQIndex(i, other.km));
311
312 if (Double.isNaN(wws) || Double.isNaN(wqs)) {
313 if (errors != null) {
314 errors.addProblem(km, "cannot.find.w.or.q");
315 }
316 }
317
318 splineWs[i] = wws;
319 splineQs[i] = wqs;
320 }
321
322 SplineInterpolator interpolator = new SplineInterpolator();
323
324 try {
325 PolynomialSplineFunction spline =
326 interpolator.interpolate(splineQs, splineWs);
327
328 return new SplineFunction(spline, splineQs, splineWs);
329 }
330 catch (MathIllegalArgumentException miae) {
331 if (errors != null) {
332 errors.addProblem(km, "spline.creation.failed");
333 }
334 log.error("spline creation failed", miae);
335 }
336
337 return null;
338 }
339
340 public double [][] interpolateWQ(
341 Row other,
342 double km,
343 int steps,
344 WstValueTable table,
345 Calculation errors
346 ) {
347 SplineFunction sf = createSpline(other, km, table, errors);
348
349 return sf != null
350 ? sf.sample(steps, km, errors)
351 : new double[2][0];
352 }
353
354
355 public double [][] interpolateWQ(
356 int steps,
357 WstValueTable table,
358 Calculation errors
359 ) {
360 SplineFunction sf = createSpline(table, errors);
361
362 return sf != null
363 ? sf.sample(steps, km, errors)
364 : new double[2][0];
365 }
366
367
368 public double getW(QPosition qPosition) {
369 int index = qPosition.index;
370 double weight = qPosition.weight;
371
372 return weight == 1.0
373 ? ws[index]
374 : Linear.weight(weight, ws[index-1], ws[index]);
375 }
376
377 public double getW(
378 Row other,
379 double km,
380 QPosition qPosition
381 ) {
382 double kmWeight = Linear.factor(km, this.km, other.km);
383
384 int index = qPosition.index;
385 double weight = qPosition.weight;
386
387 double tw, ow;
388
389 if (weight == 1.0) {
390 tw = ws[index];
391 ow = other.ws[index];
392 }
393 else {
394 tw = Linear.weight(weight, ws[index-1], ws[index]);
395 ow = Linear.weight(weight, other.ws[index-1], other.ws[index]);
396 }
397
398 return Linear.weight(kmWeight, tw, ow);
399 }
400
401 public double [] findQsForW(double w, WstValueTable table) {
402
403 TDoubleArrayList qs = new TDoubleArrayList();
404
405 if (ws.length > 0 && Math.abs(ws[0]-w) < 0.000001) {
406 double q = table.getQIndex(0, km);
407 if (!Double.isNaN(q)) {
408 qs.add(q);
409 }
410 }
411
412 for (int i = 1; i < ws.length; ++i) {
413 double w2 = ws[i];
414 if (Double.isNaN(w2)) {
415 continue;
416 }
417 if (Math.abs(w2-w) < 0.000001) {
418 double q = table.getQIndex(i, km);
419 if (!Double.isNaN(q)) {
420 qs.add(q);
421 }
422 continue;
423 }
424 double w1 = ws[i-1];
425 if (Double.isNaN(w1)) {
426 continue;
427 }
428
429 if (w < Math.min(w1, w2) || w > Math.max(w1, w2)) {
430 continue;
431 }
432
433 double q1 = table.getQIndex(i-1, km);
434 double q2 = table.getQIndex(i, km);
435 if (Double.isNaN(q1) || Double.isNaN(q2)) {
436 continue;
437 }
438
439 double q = Linear.linear(w, w1, w2, q1, q2);
440 qs.add(q);
441 }
442
443 return qs.toNativeArray();
444 }
445
446 public double [] findQsForW(
447 Row other,
448 double w,
449 double km,
450 WstValueTable table
451 ) {
452 TDoubleArrayList qs = new TDoubleArrayList();
453
454 double factor = Linear.factor(km, this.km, other.km);
455
456 if (ws.length > 0) {
457 double wt = Linear.weight(factor, ws[0], other.ws[0]);
458 if (!Double.isNaN(wt)) {
459 double q = table.getQIndex(0, km);
460 if (!Double.isNaN(q)) {
461 qs.add(q);
462 }
463 }
464 }
465
466 for (int i = 1; i < ws.length; ++i) {
467 double w2 = Linear.weight(factor, ws[i], other.ws[i]);
468 if (Double.isNaN(w2)) {
469 continue;
470 }
471 if (Math.abs(w2-w) < 0.000001) {
472 double q = table.getQIndex(i, km);
473 if (!Double.isNaN(q)) {
474 qs.add(q);
475 }
476 continue;
477 }
478 double w1 = Linear.weight(factor, ws[i-1], other.ws[i-1]);
479 if (Double.isNaN(w1)) {
480 continue;
481 }
482
483 if (w < Math.min(w1, w2) || w > Math.max(w1, w2)) {
484 continue;
485 }
486
487 double q1 = table.getQIndex(i-1, km);
488 double q2 = table.getQIndex(i, km);
489 if (Double.isNaN(q1) || Double.isNaN(q2)) {
490 continue;
491 }
492
493 double q = Linear.linear(w, w1, w2, q1, q2);
494 qs.add(q);
495 }
496
497 return qs.toNativeArray();
498 }
499
500 public double [] getMinMaxW(double [] result) {
501 double minW = Double.MAX_VALUE;
502 double maxW = -Double.MAX_VALUE;
503 for (int i = 0; i < ws.length; ++i) {
504 double w = ws[i];
505 if (w < minW) minW = w;
506 if (w > maxW) maxW = w;
507 }
508 result[0] = minW;
509 result[1] = maxW;
510 return result;
511 }
512
513 public double [] getMinMaxW(Row other, double km, double [] result) {
514 double [] m1 = this .getMinMaxW(new double [2]);
515 double [] m2 = other.getMinMaxW(new double [2]);
516 double factor = Linear.factor(km, this.km, other.km);
517 result[0] = Linear.weight(factor, m1[0], m2[0]);
518 result[1] = Linear.weight(factor, m1[1], m2[1]);
519 return result;
520 }
521 } // class Row
522
523 /** Rows in table. */
524 protected List<Row> rows;
525
526 /** Columns in table. */
527 protected Column [] columns;
528
529 public WstValueTable() {
530 rows = new ArrayList<Row>();
531 }
532
533 public WstValueTable(Column [] columns) {
534 this();
535 this.columns = columns;
536 }
537
538 public WstValueTable(Column [] columns, List<Row> rows) {
539 this.columns = columns;
540 this.rows = rows;
541 }
542
543 /** Sort rows (by km). */
544 public void sortRows() {
545 Collections.sort(rows);
546 }
547
548 /**
549 * @param km Given kilometer.
550 * @param qs Given Q values.
551 * @param ws output parameter.
552 */
553 public double [] interpolateW(double km, double [] qs, double [] ws) {
554 return interpolateW(km, qs, ws, null);
555 }
556
557
558 /**
559 * @param ws (output parameter), gets returned.
560 * @return output parameter ws.
561 */
562 public double [] interpolateW(
563 double km,
564 double [] qs,
565 double [] ws,
566 Calculation errors
567 ) {
568 int rowIndex = Collections.binarySearch(rows, new Row(km));
569
570 QPosition qPosition = new QPosition();
571
572 if (rowIndex >= 0) { // direct row match
573 Row row = rows.get(rowIndex);
574 for (int i = 0; i < qs.length; ++i) {
575 if (getQPosition(km, qs[i], qPosition) == null) {
576 if (errors != null) {
577 errors.addProblem(km, "cannot.find.q", qs[i]);
578 }
579 ws[i] = Double.NaN;
580 }
581 else {
582 if (Double.isNaN(ws[i] = row.getW(qPosition))
583 && errors != null) {
584 errors.addProblem(
585 km, "cannot.find.w.for.q", qs[i]);
586 }
587 }
588 }
589 }
590 else { // needs bilinear interpolation
591 rowIndex = -rowIndex -1;
592
593 if (rowIndex < 1 || rowIndex >= rows.size()) {
594 // do not extrapolate
595 Arrays.fill(ws, Double.NaN);
596 if (errors != null) {
597 errors.addProblem(km, "km.not.found");
598 }
599 }
600 else {
601 Row r1 = rows.get(rowIndex-1);
602 Row r2 = rows.get(rowIndex);
603 r1.interpolateW(r2, km, qs, ws, this, errors);
604 }
605 }
606
607 return ws;
608 }
609
610 public double [] getMinMaxQ(double km) {
611 return getMinMaxQ(km, new double [2]);
612 }
613
614 public double [] getMinMaxQ(double km, double [] result) {
615 double minQ = Double.MAX_VALUE;
616 double maxQ = -Double.MAX_VALUE;
617
618 for (int i = 0; i < columns.length; ++i) {
619 double q = columns[i].getQRangeTree().findQ(km);
620 if (!Double.isNaN(q)) {
621 if (q < minQ) minQ = q;
622 if (q > maxQ) maxQ = q;
623 }
624 }
625
626 if (minQ < Double.MAX_VALUE) {
627 result[0] = minQ;
628 result[1] = maxQ;
629 return result;
630 }
631
632 return null;
633 }
634
635 public double [] getMinMaxQ(double from, double to, double step) {
636 double [] result = new double[2];
637
638 double minQ = Double.MAX_VALUE;
639 double maxQ = -Double.MAX_VALUE;
640
641 if (from > to) {
642 double tmp = from;
643 from = to;
644 to = tmp;
645 }
646
647 step = Math.max(Math.abs(step), 0.0001);
648
649 double d = from;
650 for (; d <= to; d += step) {
651 if (getMinMaxQ(d, result) != null) {
652 if (result[0] < minQ) minQ = result[0];
653 if (result[1] > maxQ) maxQ = result[1];
654 }
655 }
656
657 if (d != to) {
658 if (getMinMaxQ(to, result) != null) {
659 if (result[0] < minQ) minQ = result[0];
660 if (result[1] > maxQ) maxQ = result[1];
661 }
662 }
663
664 return minQ < Double.MAX_VALUE
665 ? new double [] { minQ, maxQ }
666 : null;
667 }
668
669 public double [] getMinMaxW(double km) {
670 return getMinMaxW(km, new double [2]);
671
672 }
673 public double [] getMinMaxW(double km, double [] result) {
674 int rowIndex = Collections.binarySearch(rows, new Row(km));
675
676 if (rowIndex >= 0) {
677 return rows.get(rowIndex).getMinMaxW(result);
678 }
679
680 rowIndex = -rowIndex -1;
681
682 if (rowIndex < 1 || rowIndex >= rows.size()) {
683 // do not extrapolate
684 return null;
685 }
686
687 Row r1 = rows.get(rowIndex-1);
688 Row r2 = rows.get(rowIndex);
689
690 return r1.getMinMaxW(r2, km, result);
691 }
692
693 public double [] getMinMaxW(double from, double to, double step) {
694 double [] result = new double[2];
695
696 double minW = Double.MAX_VALUE;
697 double maxW = -Double.MAX_VALUE;
698
699 if (from > to) {
700 double tmp = from;
701 from = to;
702 to = tmp;
703 }
704
705 step = Math.max(Math.abs(step), 0.0001);
706
707 double d = from;
708 for (; d <= to; d += step) {
709 if (getMinMaxW(d, result) != null) {
710 if (result[0] < minW) minW = result[0];
711 if (result[1] > maxW) maxW = result[1];
712 }
713 }
714
715 if (d != to) {
716 if (getMinMaxW(to, result) != null) {
717 if (result[0] < minW) minW = result[0];
718 if (result[1] > maxW) maxW = result[1];
719 }
720 }
721
722 return minW < Double.MAX_VALUE
723 ? new double [] { minW, maxW }
724 : null;
725 }
726
727 /**
728 * Interpolate W and Q values at a given km.
729 */
730 public double [][] interpolateWQ(double km) {
731 return interpolateWQ(km, null);
732 }
733
734 /**
735 * Interpolate W and Q values at a given km.
736 *
737 * @param errors where to store errors.
738 *
739 * @return double double array, first index Ws, second Qs.
740 */
741 public double [][] interpolateWQ(double km, Calculation errors) {
742 return interpolateWQ(km, DEFAULT_Q_STEPS, errors);
743 }
744
745
746 /**
747 * Interpolate W and Q values at a given km.
748 */
749 public double [][] interpolateWQ(double km, int steps, Calculation errors) {
750
751 int rowIndex = Collections.binarySearch(rows, new Row(km));
752
753 if (rowIndex >= 0) { // direct row match
754 Row row = rows.get(rowIndex);
755 return row.interpolateWQ(steps, this, errors);
756 }
757
758 rowIndex = -rowIndex -1;
759
760 if (rowIndex < 1 || rowIndex >= rows.size()) {
761 // do not extrapolate
762 if (errors != null) {
763 errors.addProblem(km, "km.not.found");
764 }
765 return new double[2][0];
766 }
767
768 Row r1 = rows.get(rowIndex-1);
769 Row r2 = rows.get(rowIndex);
770
771 return r1.interpolateWQ(r2, km, steps, this, errors);
772 }
773
774 public boolean interpolate(
775 double km,
776 double [] out,
777 QPosition qPosition,
778 Function qFunction
779 ) {
780 int R1 = rows.size()-1;
781
782 out[1] = qFunction.value(getQ(qPosition, km));
783
784 if (Double.isNaN(out[1])) {
785 return false;
786 }
787
788 QPosition nPosition = new QPosition();
789 if (getQPosition(km, out[1], nPosition) == null) {
790 return false;
791 }
792
793 int rowIndex = Collections.binarySearch(rows, new Row(km));
794
795 if (rowIndex >= 0) {
796 // direct row match
797 out[0] = rows.get(rowIndex).getW(nPosition);
798 return !Double.isNaN(out[0]);
799 }
800
801 rowIndex = -rowIndex -1;
802
803 if (rowIndex < 1 || rowIndex > R1) {
804 // do not extrapolate
805 return false;
806 }
807
808 Row r1 = rows.get(rowIndex-1);
809 Row r2 = rows.get(rowIndex);
810 out[0] = r1.getW(r2, km, nPosition);
811
812 return !Double.isNaN(out[0]);
813 }
814
815
816 /**
817 * Look up interpolation of a Q at given positions.
818 *
819 * @param q the non-interpolated Q value.
820 * @param referenceKm the reference km (e.g. gauge position).
821 * @param kms positions for which to interpolate.
822 * @param ws (output) resulting interpolated ws.
823 * @param qs (output) resulting interpolated qs.
824 * @param errors calculation object to store errors.
825 */
826 public QPosition interpolate(
827 double q,
828 double referenceKm,
829 double [] kms,
830 double [] ws,
831 double [] qs,
832 Calculation errors
833 ) {
834 return interpolate(
835 q, referenceKm, kms, ws, qs, 0, kms.length, errors);
836 }
837
838 public QPosition interpolate(
839 double q,
840 double referenceKm,
841 double [] kms,
842 double [] ws,
843 double [] qs,
844 int startIndex,
845 int length,
846 Calculation errors
847 ) {
848 QPosition qPosition = getQPosition(referenceKm, q);
849
850 if (qPosition == null) {
851 // we cannot locate q at km
852 Arrays.fill(ws, Double.NaN);
853 Arrays.fill(qs, Double.NaN);
854 if (errors != null) {
855 errors.addProblem(referenceKm, "cannot.find.q", q);
856 }
857 return null;
858 }
859
860 Row kmKey = new Row();
861
862 int R1 = rows.size()-1;
863
864 for (int i = startIndex, end = startIndex+length; i < end; ++i) {
865
866 if (Double.isNaN(qs[i] = getQ(qPosition, kms[i]))) {
867 if (errors != null) {
868 errors.addProblem(kms[i], "cannot.find.q", q);
869 }
870 ws[i] = Double.NaN;
871 continue;
872 }
873
874 kmKey.km = kms[i];
875 int rowIndex = Collections.binarySearch(rows, kmKey);
876
877 if (rowIndex >= 0) {
878 // direct row match
879 if (Double.isNaN(ws[i] = rows.get(rowIndex).getW(qPosition))
880 && errors != null) {
881 errors.addProblem(kms[i], "cannot.find.w.for.q", q);
882 }
883 continue;
884 }
885
886 rowIndex = -rowIndex -1;
887
888 if (rowIndex < 1 || rowIndex > R1) {
889 // do not extrapolate
890 if (errors != null) {
891 errors.addProblem(kms[i], "km.not.found");
892 }
893 ws[i] = Double.NaN;
894 continue;
895 }
896 Row r1 = rows.get(rowIndex-1);
897 Row r2 = rows.get(rowIndex);
898
899 if (Double.isNaN(ws[i] = r1.getW(r2, kms[i], qPosition))
900 && errors != null) {
901 errors.addProblem(kms[i], "cannot.find.w.for.q", q);
902 }
903 }
904
905 return qPosition;
906 }
907
908 /**
909 * Linearly interpolate w at a km at a column of two rows.
910 *
911 * @param km position for which to interpolate.
912 * @param row1 first row.
913 * @param row2 second row.
914 * @param col column-index at which to look.
915 *
916 * @return Linearly interpolated w, NaN if one of the given rows was null.
917 */
918 public static double linearW(double km, Row row1, Row row2, int col) {
919 if (row1 == null || row2 == null) {
920 return Double.NaN;
921 }
922
923 return Linear.linear(km,
924 row1.km, row2.km,
925 row1.ws[col], row2.ws[col]);
926 }
927
928 /**
929 * Do interpolation/lookup of W and Q within columns (i.e. ignoring values
930 * of other columns).
931 * @param km position (km) at which to interpolate/lookup.
932 * @return [[q0, q1, .. qx] , [w0, w1, .. wx]] (can contain NaNs)
933 */
934 public double [][] interpolateWQColumnwise(double km) {
935 log.debug("WstValueTable.interpolateWQColumnwise");
936 double [] qs = new double[columns.length];
937 double [] ws = new double[columns.length];
938
939 // Find out row from where we will start searching.
940 int rowIndex = Collections.binarySearch(rows, new Row(km));
941
942 if (rowIndex < 0) {
943 rowIndex = -rowIndex -1;
944 }
945
946 // TODO Beyond definition, we could stop more clever.
947 if (rowIndex >= rows.size()) {
948 rowIndex = rows.size() -1;
949 }
950
951 Row startRow = rows.get(rowIndex);
952
953 for (int col = 0; col < columns.length; col++) {
954 qs[col] = columns[col].getQRangeTree().findQ(km);
955 if (startRow.km == km && startRow.ws[col] != Double.NaN) {
956 // Great. W is defined at km.
957 ws[col] = startRow.ws[col];
958 continue;
959 }
960
961 // Search neighbouring rows that define w at this col.
962 Row rowBefore = null;
963 Row rowAfter = null;
964 for (int before = rowIndex -1; before >= 0; before--) {
965 if (!Double.isNaN(rows.get(before).ws[col])) {
966 rowBefore = rows.get(before);
967 break;
968 }
969 }
970 if (rowBefore != null) {
971 for (int after = rowIndex, R = rows.size(); after < R; after++) {
972 if (!Double.isNaN(rows.get(after).ws[col])) {
973 rowAfter = rows.get(after);
974 break;
975 }
976 }
977 }
978
979 ws[col] = linearW(km, rowBefore, rowAfter, col);
980 }
981
982 return new double [][] {qs, ws};
983 }
984
985 public double [] findQsForW(double km, double w) {
986
987 int rowIndex = Collections.binarySearch(rows, new Row(km));
988
989 if (rowIndex >= 0) {
990 return rows.get(rowIndex).findQsForW(w, this);
991 }
992
993 rowIndex = -rowIndex - 1;
994
995 if (rowIndex < 1 || rowIndex >= rows.size()) {
996 // Do not extrapolate.
997 return new double[0];
998 }
999
1000 // Needs bilinear interpolation.
1001 Row r1 = rows.get(rowIndex-1);
1002 Row r2 = rows.get(rowIndex);
1003
1004 return r1.findQsForW(r2, w, km, this);
1005 }
1006
1007 protected SplineFunction createSpline(double km, Calculation errors) {
1008
1009 int rowIndex = Collections.binarySearch(rows, new Row(km));
1010
1011 if (rowIndex >= 0) {
1012 SplineFunction sf = rows.get(rowIndex).createSpline(this, errors);
1013 if (sf == null && errors != null) {
1014 errors.addProblem(km, "cannot.create.wq.relation");
1015 }
1016 return sf;
1017 }
1018
1019 rowIndex = -rowIndex - 1;
1020
1021 if (rowIndex < 1 || rowIndex >= rows.size()) {
1022 // Do not extrapolate.
1023 if (errors != null) {
1024 errors.addProblem(km, "km.not.found");
1025 }
1026 return null;
1027 }
1028
1029 // Needs bilinear interpolation.
1030 Row r1 = rows.get(rowIndex-1);
1031 Row r2 = rows.get(rowIndex);
1032
1033 SplineFunction sf = r1.createSpline(r2, km, this, errors);
1034 if (sf == null && errors != null) {
1035 errors.addProblem(km, "cannot.create.wq.relation");
1036 }
1037
1038 return sf;
1039 }
1040
1041 /** 'Bezugslinienverfahren' */
1042 public double [][] relateWs(
1043 double km1,
1044 double km2,
1045 Calculation errors
1046 ) {
1047 return relateWs(km1, km2, RELATE_WS_SAMPLES, errors);
1048 }
1049
1050 private static class ErrorHandler {
1051
1052 boolean hasErrors;
1053 Calculation errors;
1054
1055 ErrorHandler(Calculation errors) {
1056 this.errors = errors;
1057 }
1058
1059 void error(double km, String key, Object ... args) {
1060 if (errors != null && !hasErrors) {
1061 hasErrors = true;
1062 errors.addProblem(km, key, args);
1063 }
1064 }
1065 } // class ErrorHandler
1066
1067
1068 /* TODO: Add optimized methods of relateWs to relate one
1069 * start km to many end kms. The index generation/spline stuff for
1070 * the start km is always the same.
1071 */
1072 public double [][] relateWs(
1073 double km1,
1074 double km2,
1075 int numSamples,
1076 Calculation errors
1077 ) {
1078 SplineFunction sf1 = createSpline(km1, errors);
1079 if (sf1 == null) {
1080 return new double[2][0];
1081 }
1082
1083 SplineFunction sf2 = createSpline(km2, errors);
1084 if (sf2 == null) {
1085 return new double[2][0];
1086 }
1087
1088 PolynomialSplineFunction iQ1 = sf1.createIndexQRelation();
1089 if (iQ1 == null) {
1090 if (errors != null) {
1091 errors.addProblem(km1, "cannot.create.index.q.relation");
1092 }
1093 return new double[2][0];
1094 }
1095
1096 PolynomialSplineFunction iQ2 = sf2.createIndexQRelation();
1097 if (iQ2 == null) {
1098 if (errors != null) {
1099 errors.addProblem(km2, "cannot.create.index.q.relation");
1100 }
1101 return new double[2][0];
1102 }
1103
1104 int N = Math.min(sf1.splineQs.length, sf2.splineQs.length);
1105 double stepWidth = N/(double)numSamples;
1106
1107 PolynomialSplineFunction qW1 = sf1.spline;
1108 PolynomialSplineFunction qW2 = sf2.spline;
1109
1110 TDoubleArrayList ws1 = new TDoubleArrayList(numSamples);
1111 TDoubleArrayList ws2 = new TDoubleArrayList(numSamples);
1112 TDoubleArrayList qs1 = new TDoubleArrayList(numSamples);
1113 TDoubleArrayList qs2 = new TDoubleArrayList(numSamples);
1114
1115 ErrorHandler err = new ErrorHandler(errors);
1116
1117 int i = 0;
1118 for (double p = 0d; p <= N-1; p += stepWidth, ++i) {
1119
1120 double q1;
1121 try {
1122 q1 = iQ1.value(p);
1123 }
1124 catch (ArgumentOutsideDomainException aode) {
1125 err.error(km1, "w.w.qkm1.failed", p);
1126 continue;
1127 }
1128
1129 double w1;
1130 try {
1131 w1 = qW1.value(q1);
1132 }
1133 catch (ArgumentOutsideDomainException aode) {
1134 err.error(km1, "w.w.wkm1.failed", q1, p);
1135 continue;
1136 }
1137
1138 double q2;
1139 try {
1140 q2 = iQ2.value(p);
1141 }
1142 catch (ArgumentOutsideDomainException aode) {
1143 err.error(km2, "w.w.qkm2.failed", p);
1144 continue;
1145 }
1146
1147 double w2;
1148 try {
1149 w2 = qW2.value(q2);
1150 }
1151 catch (ArgumentOutsideDomainException aode) {
1152 err.error(km2, "w.w.wkm2.failed", q2, p);
1153 continue;
1154 }
1155
1156 ws1.add(w1);
1157 ws2.add(w2);
1158 qs1.add(q1);
1159 qs2.add(q2);
1160 }
1161
1162 return new double [][] {
1163 ws1.toNativeArray(),
1164 qs1.toNativeArray(),
1165 ws2.toNativeArray(),
1166 qs2.toNativeArray() };
1167 }
1168
1169 public QPosition getQPosition(double km, double q) {
1170 return getQPosition(km, q, new QPosition());
1171 }
1172
1173 public QPosition getQPosition(double km, double q, QPosition qPosition) {
1174
1175 if (columns.length == 0) {
1176 return null;
1177 }
1178
1179 double qLast = columns[0].getQRangeTree().findQ(km);
1180
1181 if (Math.abs(qLast - q) < 0.00001) {
1182 return qPosition.set(0, 1d);
1183 }
1184
1185 for (int i = 1; i < columns.length; ++i) {
1186 double qCurrent = columns[i].getQRangeTree().findQ(km);
1187 if (Math.abs(qCurrent - q) < 0.00001) {
1188 return qPosition.set(i, 1d);
1189 }
1190
1191 double qMin, qMax;
1192 if (qLast < qCurrent) { qMin = qLast; qMax = qCurrent; }
1193 else { qMin = qCurrent; qMax = qLast; }
1194
1195 if (q > qMin && q < qMax) {
1196 double weight = Linear.factor(q, qLast, qCurrent);
1197 return qPosition.set(i, weight);
1198 }
1199 qLast = qCurrent;
1200 }
1201
1202 return null;
1203 }
1204
1205 public double getQIndex(int index, double km) {
1206 return columns[index].getQRangeTree().findQ(km);
1207 }
1208
1209 public double getQ(QPosition qPosition, double km) {
1210 int index = qPosition.index;
1211 double weight = qPosition.weight;
1212
1213 if (weight == 1d) {
1214 return columns[index].getQRangeTree().findQ(km);
1215 }
1216 double q1 = columns[index-1].getQRangeTree().findQ(km);
1217 double q2 = columns[index ].getQRangeTree().findQ(km);
1218 return Linear.weight(weight, q1, q2);
1219 }
1220 }
1221 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
http://dive4elements.wald.intevation.org