Mercurial > dive4elements > river

comparison flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java @ 3818:dc18457b1cef

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merged flys-artifacts/pre2.7-2012-03-16
author Thomas Arendsen Hein <thomas@intevation.de>
date Fri, 28 Sep 2012 12:14:59 +0200
parents 0d8c97374dc9
children 5642a83420f2
comparison
equal deleted inserted replaced
2456:60ab1054069d 3818:dc18457b1cef
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 public double [][] interpolateWQ(double km, Calculation errors) {
738 return interpolateWQ(km, DEFAULT_Q_STEPS, errors);
739 }
740
741 /**
742 * Interpolate W and Q values at a given km.
743 */
744 public double [][] interpolateWQ(double km, int steps, Calculation errors) {
745
746 int rowIndex = Collections.binarySearch(rows, new Row(km));
747
748 if (rowIndex >= 0) { // direct row match
749 Row row = rows.get(rowIndex);
750 return row.interpolateWQ(steps, this, errors);
751 }
752
753 rowIndex = -rowIndex -1;
754
755 if (rowIndex < 1 || rowIndex >= rows.size()) {
756 // do not extrapolate
757 if (errors != null) {
758 errors.addProblem(km, "km.not.found");
759 }
760 return new double[2][0];
761 }
762
763 Row r1 = rows.get(rowIndex-1);
764 Row r2 = rows.get(rowIndex);
765
766 return r1.interpolateWQ(r2, km, steps, this, errors);
767 }
768
769 public boolean interpolate(
770 double km,
771 double [] out,
772 QPosition qPosition,
773 Function qFunction
774 ) {
775 int R1 = rows.size()-1;
776
777 out[1] = qFunction.value(getQ(qPosition, km));
778
779 if (Double.isNaN(out[1])) {
780 return false;
781 }
782
783 QPosition nPosition = new QPosition();
784 if (getQPosition(km, out[1], nPosition) == null) {
785 return false;
786 }
787
788 int rowIndex = Collections.binarySearch(rows, new Row(km));
789
790 if (rowIndex >= 0) {
791 // direct row match
792 out[0] = rows.get(rowIndex).getW(nPosition);
793 return !Double.isNaN(out[0]);
794 }
795
796 rowIndex = -rowIndex -1;
797
798 if (rowIndex < 1 || rowIndex > R1) {
799 // do not extrapolate
800 return false;
801 }
802
803 Row r1 = rows.get(rowIndex-1);
804 Row r2 = rows.get(rowIndex);
805 out[0] = r1.getW(r2, km, nPosition);
806
807 return !Double.isNaN(out[0]);
808 }
809
810
811 /**
812 * Look up interpolation of a Q at given positions.
813 *
814 * @param q the non-interpolated Q value.
815 * @param referenceKm the reference km (e.g. gauge position).
816 * @param kms positions for which to interpolate.
817 * @param ws (output) resulting interpolated ws.
818 * @param qs (output) resulting interpolated qs.
819 * @param errors calculation object to store errors.
820 */
821 public QPosition interpolate(
822 double q,
823 double referenceKm,
824 double [] kms,
825 double [] ws,
826 double [] qs,
827 Calculation errors
828 ) {
829 return interpolate(
830 q, referenceKm, kms, ws, qs, 0, kms.length, errors);
831 }
832
833 public QPosition interpolate(
834 double q,
835 double referenceKm,
836 double [] kms,
837 double [] ws,
838 double [] qs,
839 int startIndex,
840 int length,
841 Calculation errors
842 ) {
843 QPosition qPosition = getQPosition(referenceKm, q);
844
845 if (qPosition == null) {
846 // we cannot locate q at km
847 Arrays.fill(ws, Double.NaN);
848 Arrays.fill(qs, Double.NaN);
849 if (errors != null) {
850 errors.addProblem(referenceKm, "cannot.find.q", q);
851 }
852 return null;
853 }
854
855 Row kmKey = new Row();
856
857 int R1 = rows.size()-1;
858
859 for (int i = startIndex, end = startIndex+length; i < end; ++i) {
860
861 if (Double.isNaN(qs[i] = getQ(qPosition, kms[i]))) {
862 if (errors != null) {
863 errors.addProblem(kms[i], "cannot.find.q", q);
864 }
865 ws[i] = Double.NaN;
866 continue;
867 }
868
869 kmKey.km = kms[i];
870 int rowIndex = Collections.binarySearch(rows, kmKey);
871
872 if (rowIndex >= 0) {
873 // direct row match
874 if (Double.isNaN(ws[i] = rows.get(rowIndex).getW(qPosition))
875 && errors != null) {
876 errors.addProblem(kms[i], "cannot.find.w.for.q", q);
877 }
878 continue;
879 }
880
881 rowIndex = -rowIndex -1;
882
883 if (rowIndex < 1 || rowIndex > R1) {
884 // do not extrapolate
885 if (errors != null) {
886 errors.addProblem(kms[i], "km.not.found");
887 }
888 ws[i] = Double.NaN;
889 continue;
890 }
891 Row r1 = rows.get(rowIndex-1);
892 Row r2 = rows.get(rowIndex);
893
894 if (Double.isNaN(ws[i] = r1.getW(r2, kms[i], qPosition))
895 && errors != null) {
896 errors.addProblem(kms[i], "cannot.find.w.for.q", q);
897 }
898 }
899
900 return qPosition;
901 }
902
903 /**
904 * Linearly interpolate w at a km at a column of two rows.
905 *
906 * @param km position for which to interpolate.
907 * @param row1 first row.
908 * @param row2 second row.
909 * @param col column-index at which to look.
910 *
911 * @return Linearly interpolated w, NaN if one of the given rows was null.
912 */
913 public static double linearW(double km, Row row1, Row row2, int col) {
914 if (row1 == null || row2 == null) {
915 return Double.NaN;
916 }
917
918 return Linear.linear(km,
919 row1.km, row2.km,
920 row1.ws[col], row2.ws[col]);
921 }
922
923 /**
924 * Do interpolation/lookup of W and Q within columns (i.e. ignoring values
925 * of other columns).
926 * @param km position (km) at which to interpolate/lookup.
927 * @return [[q0, q1, .. qx] , [w0, w1, .. wx]] (can contain NaNs)
928 */
929 public double [][] interpolateWQColumnwise(double km) {
930 log.debug("WstValueTable.interpolateWQColumnwise");
931 double [] qs = new double[columns.length];
932 double [] ws = new double[columns.length];
933
934 // Find out row from where we will start searching.
935 int rowIndex = Collections.binarySearch(rows, new Row(km));
936
937 if (rowIndex < 0) {
938 rowIndex = -rowIndex -1;
939 }
940
941 // TODO Beyond definition, we could stop more clever.
942 if (rowIndex >= rows.size()) {
943 rowIndex = rows.size() -1;
944 }
945
946 Row startRow = rows.get(rowIndex);
947
948 for (int col = 0; col < columns.length; col++) {
949 qs[col] = columns[col].getQRangeTree().findQ(km);
950 if (startRow.km == km && startRow.ws[col] != Double.NaN) {
951 // Great. W is defined at km.
952 ws[col] = startRow.ws[col];
953 continue;
954 }
955
956 // Search neighbouring rows that define w at this col.
957 Row rowBefore = null;
958 Row rowAfter = null;
959 for (int before = rowIndex -1; before >= 0; before--) {
960 if (!Double.isNaN(rows.get(before).ws[col])) {
961 rowBefore = rows.get(before);
962 break;
963 }
964 }
965 if (rowBefore != null) {
966 for (int after = rowIndex, R = rows.size(); after < R; after++) {
967 if (!Double.isNaN(rows.get(after).ws[col])) {
968 rowAfter = rows.get(after);
969 break;
970 }
971 }
972 }
973
974 ws[col] = linearW(km, rowBefore, rowAfter, col);
975 }
976
977 return new double [][] {qs, ws};
978 }
979
980 public double [] findQsForW(double km, double w) {
981
982 int rowIndex = Collections.binarySearch(rows, new Row(km));
983
984 if (rowIndex >= 0) {
985 return rows.get(rowIndex).findQsForW(w, this);
986 }
987
988 rowIndex = -rowIndex - 1;
989
990 if (rowIndex < 1 || rowIndex >= rows.size()) {
991 // Do not extrapolate.
992 return new double[0];
993 }
994
995 // Needs bilinear interpolation.
996 Row r1 = rows.get(rowIndex-1);
997 Row r2 = rows.get(rowIndex);
998
999 return r1.findQsForW(r2, w, km, this);
1000 }
1001
1002 protected SplineFunction createSpline(double km, Calculation errors) {
1003
1004 int rowIndex = Collections.binarySearch(rows, new Row(km));
1005
1006 if (rowIndex >= 0) {
1007 SplineFunction sf = rows.get(rowIndex).createSpline(this, errors);
1008 if (sf == null && errors != null) {
1009 errors.addProblem(km, "cannot.create.wq.relation");
1010 }
1011 return sf;
1012 }
1013
1014 rowIndex = -rowIndex - 1;
1015
1016 if (rowIndex < 1 || rowIndex >= rows.size()) {
1017 // Do not extrapolate.
1018 if (errors != null) {
1019 errors.addProblem(km, "km.not.found");
1020 }
1021 return null;
1022 }
1023
1024 // Needs bilinear interpolation.
1025 Row r1 = rows.get(rowIndex-1);
1026 Row r2 = rows.get(rowIndex);
1027
1028 SplineFunction sf = r1.createSpline(r2, km, this, errors);
1029 if (sf == null && errors != null) {
1030 errors.addProblem(km, "cannot.create.wq.relation");
1031 }
1032
1033 return sf;
1034 }
1035
1036 /** 'Bezugslinienverfahren' */
1037 public double [][] relateWs(
1038 double km1,
1039 double km2,
1040 Calculation errors
1041 ) {
1042 return relateWs(km1, km2, RELATE_WS_SAMPLES, errors);
1043 }
1044
1045 private static class ErrorHandler {
1046
1047 boolean hasErrors;
1048 Calculation errors;
1049
1050 ErrorHandler(Calculation errors) {
1051 this.errors = errors;
1052 }
1053
1054 void error(double km, String key, Object ... args) {
1055 if (errors != null && !hasErrors) {
1056 hasErrors = true;
1057 errors.addProblem(km, key, args);
1058 }
1059 }
1060 } // class ErrorHandler
1061
1062
1063 /* TODO: Add optimized methods of relateWs to relate one
1064 * start km to many end kms. The index generation/spline stuff for
1065 * the start km is always the same.
1066 */
1067 public double [][] relateWs(
1068 double km1,
1069 double km2,
1070 int numSamples,
1071 Calculation errors
1072 ) {
1073 SplineFunction sf1 = createSpline(km1, errors);
1074 if (sf1 == null) {
1075 return new double[2][0];
1076 }
1077
1078 SplineFunction sf2 = createSpline(km2, errors);
1079 if (sf2 == null) {
1080 return new double[2][0];
1081 }
1082
1083 PolynomialSplineFunction iQ1 = sf1.createIndexQRelation();
1084 if (iQ1 == null) {
1085 if (errors != null) {
1086 errors.addProblem(km1, "cannot.create.index.q.relation");
1087 }
1088 return new double[2][0];
1089 }
1090
1091 PolynomialSplineFunction iQ2 = sf2.createIndexQRelation();
1092 if (iQ2 == null) {
1093 if (errors != null) {
1094 errors.addProblem(km2, "cannot.create.index.q.relation");
1095 }
1096 return new double[2][0];
1097 }
1098
1099 int N = Math.min(sf1.splineQs.length, sf2.splineQs.length);
1100 double stepWidth = N/(double)numSamples;
1101
1102 PolynomialSplineFunction qW1 = sf1.spline;
1103 PolynomialSplineFunction qW2 = sf2.spline;
1104
1105 TDoubleArrayList ws1 = new TDoubleArrayList(numSamples);
1106 TDoubleArrayList ws2 = new TDoubleArrayList(numSamples);
1107 TDoubleArrayList qs1 = new TDoubleArrayList(numSamples);
1108 TDoubleArrayList qs2 = new TDoubleArrayList(numSamples);
1109
1110 ErrorHandler err = new ErrorHandler(errors);
1111
1112 int i = 0;
1113 for (double p = 0d; p <= N-1; p += stepWidth, ++i) {
1114
1115 double q1;
1116 try {
1117 q1 = iQ1.value(p);
1118 }
1119 catch (ArgumentOutsideDomainException aode) {
1120 err.error(km1, "w.w.qkm1.failed", p);
1121 continue;
1122 }
1123
1124 double w1;
1125 try {
1126 w1 = qW1.value(q1);
1127 }
1128 catch (ArgumentOutsideDomainException aode) {
1129 err.error(km1, "w.w.wkm1.failed", q1, p);
1130 continue;
1131 }
1132
1133 double q2;
1134 try {
1135 q2 = iQ2.value(p);
1136 }
1137 catch (ArgumentOutsideDomainException aode) {
1138 err.error(km2, "w.w.qkm2.failed", p);
1139 continue;
1140 }
1141
1142 double w2;
1143 try {
1144 w2 = qW2.value(q2);
1145 }
1146 catch (ArgumentOutsideDomainException aode) {
1147 err.error(km2, "w.w.wkm2.failed", q2, p);
1148 continue;
1149 }
1150
1151 ws1.add(w1);
1152 ws2.add(w2);
1153 qs1.add(q1);
1154 qs2.add(q2);
1155 }
1156
1157 return new double [][] {
1158 ws1.toNativeArray(),
1159 qs1.toNativeArray(),
1160 ws2.toNativeArray(),
1161 qs2.toNativeArray() };
1162 }
1163
1164 public QPosition getQPosition(double km, double q) {
1165 return getQPosition(km, q, new QPosition());
1166 }
1167
1168 public QPosition getQPosition(double km, double q, QPosition qPosition) {
1169
1170 if (columns.length == 0) {
1171 return null;
1172 }
1173
1174 double qLast = columns[0].getQRangeTree().findQ(km);
1175
1176 if (Math.abs(qLast - q) < 0.00001) {
1177 return qPosition.set(0, 1d);
1178 }
1179
1180 for (int i = 1; i < columns.length; ++i) {
1181 double qCurrent = columns[i].getQRangeTree().findQ(km);
1182 if (Math.abs(qCurrent - q) < 0.00001) {
1183 return qPosition.set(i, 1d);
1184 }
1185
1186 double qMin, qMax;
1187 if (qLast < qCurrent) { qMin = qLast; qMax = qCurrent; }
1188 else { qMin = qCurrent; qMax = qLast; }
1189
1190 if (q > qMin && q < qMax) {
1191 double weight = Linear.factor(q, qLast, qCurrent);
1192 return qPosition.set(i, weight);
1193 }
1194 qLast = qCurrent;
1195 }
1196
1197 return null;
1198 }
1199
1200 public double getQIndex(int index, double km) {
1201 return columns[index].getQRangeTree().findQ(km);
1202 }
1203
1204 public double getQ(QPosition qPosition, double km) {
1205 int index = qPosition.index;
1206 double weight = qPosition.weight;
1207
1208 if (weight == 1d) {
1209 return columns[index].getQRangeTree().findQ(km);
1210 }
1211 double q1 = columns[index-1].getQRangeTree().findQ(km);
1212 double q2 = columns[index ].getQRangeTree().findQ(km);
1213 return Linear.weight(weight, q1, q2);
1214 }
1215 }
1216 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
http://dive4elements.wald.intevation.org