Mercurial > dive4elements > river

comparison flys-artifacts/src/main/java/org/dive4elements/river/artifacts/model/WstValueTable.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/WstValueTable.java@43e69af28b3c
children
comparison
equal deleted inserted replaced
5830:160f53ee0870 5831:bd047b71ab37
1 package org.dive4elements.river.artifacts.model;
2
3 import java.io.Serializable;
4
5 import org.dive4elements.river.artifacts.math.Linear;
6 import org.dive4elements.river.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 public Column [] getColumns() {
544 return columns;
545 }
546
547 /** Sort rows (by km). */
548 public void sortRows() {
549 Collections.sort(rows);
550 }
551
552 /**
553 * @param km Given kilometer.
554 * @param qs Given Q values.
555 * @param ws output parameter.
556 */
557 public double [] interpolateW(double km, double [] qs, double [] ws) {
558 return interpolateW(km, qs, ws, null);
559 }
560
561
562 /**
563 * @param ws (output parameter), gets returned.
564 * @return output parameter ws.
565 */
566 public double [] interpolateW(
567 double km,
568 double [] qs,
569 double [] ws,
570 Calculation errors
571 ) {
572 int rowIndex = Collections.binarySearch(rows, new Row(km));
573
574 QPosition qPosition = new QPosition();
575
576 if (rowIndex >= 0) { // direct row match
577 Row row = rows.get(rowIndex);
578 for (int i = 0; i < qs.length; ++i) {
579 if (getQPosition(km, qs[i], qPosition) == null) {
580 if (errors != null) {
581 errors.addProblem(km, "cannot.find.q", qs[i]);
582 }
583 ws[i] = Double.NaN;
584 }
585 else {
586 if (Double.isNaN(ws[i] = row.getW(qPosition))
587 && errors != null) {
588 errors.addProblem(
589 km, "cannot.find.w.for.q", qs[i]);
590 }
591 }
592 }
593 }
594 else { // needs bilinear interpolation
595 rowIndex = -rowIndex -1;
596
597 if (rowIndex < 1 || rowIndex >= rows.size()) {
598 // do not extrapolate
599 Arrays.fill(ws, Double.NaN);
600 if (errors != null) {
601 errors.addProblem(km, "km.not.found");
602 }
603 }
604 else {
605 Row r1 = rows.get(rowIndex-1);
606 Row r2 = rows.get(rowIndex);
607 r1.interpolateW(r2, km, qs, ws, this, errors);
608 }
609 }
610
611 return ws;
612 }
613
614 public double [] getMinMaxQ(double km) {
615 return getMinMaxQ(km, new double [2]);
616 }
617
618 public double [] getMinMaxQ(double km, double [] result) {
619 double minQ = Double.MAX_VALUE;
620 double maxQ = -Double.MAX_VALUE;
621
622 for (int i = 0; i < columns.length; ++i) {
623 double q = columns[i].getQRangeTree().findQ(km);
624 if (!Double.isNaN(q)) {
625 if (q < minQ) minQ = q;
626 if (q > maxQ) maxQ = q;
627 }
628 }
629
630 if (minQ < Double.MAX_VALUE) {
631 result[0] = minQ;
632 result[1] = maxQ;
633 return result;
634 }
635
636 return null;
637 }
638
639 public double [] getMinMaxQ(double from, double to, double step) {
640 double [] result = new double[2];
641
642 double minQ = Double.MAX_VALUE;
643 double maxQ = -Double.MAX_VALUE;
644
645 if (from > to) {
646 double tmp = from;
647 from = to;
648 to = tmp;
649 }
650
651 step = Math.max(Math.abs(step), 0.0001);
652
653 double d = from;
654 for (; d <= to; d += step) {
655 if (getMinMaxQ(d, result) != null) {
656 if (result[0] < minQ) minQ = result[0];
657 if (result[1] > maxQ) maxQ = result[1];
658 }
659 }
660
661 if (d != to) {
662 if (getMinMaxQ(to, result) != null) {
663 if (result[0] < minQ) minQ = result[0];
664 if (result[1] > maxQ) maxQ = result[1];
665 }
666 }
667
668 return minQ < Double.MAX_VALUE
669 ? new double [] { minQ, maxQ }
670 : null;
671 }
672
673 public double [] getMinMaxW(double km) {
674 return getMinMaxW(km, new double [2]);
675
676 }
677 public double [] getMinMaxW(double km, double [] result) {
678 int rowIndex = Collections.binarySearch(rows, new Row(km));
679
680 if (rowIndex >= 0) {
681 return rows.get(rowIndex).getMinMaxW(result);
682 }
683
684 rowIndex = -rowIndex -1;
685
686 if (rowIndex < 1 || rowIndex >= rows.size()) {
687 // do not extrapolate
688 return null;
689 }
690
691 Row r1 = rows.get(rowIndex-1);
692 Row r2 = rows.get(rowIndex);
693
694 return r1.getMinMaxW(r2, km, result);
695 }
696
697 public double [] getMinMaxW(double from, double to, double step) {
698 double [] result = new double[2];
699
700 double minW = Double.MAX_VALUE;
701 double maxW = -Double.MAX_VALUE;
702
703 if (from > to) {
704 double tmp = from;
705 from = to;
706 to = tmp;
707 }
708
709 step = Math.max(Math.abs(step), 0.0001);
710
711 double d = from;
712 for (; d <= to; d += step) {
713 if (getMinMaxW(d, result) != null) {
714 if (result[0] < minW) minW = result[0];
715 if (result[1] > maxW) maxW = result[1];
716 }
717 }
718
719 if (d != to) {
720 if (getMinMaxW(to, result) != null) {
721 if (result[0] < minW) minW = result[0];
722 if (result[1] > maxW) maxW = result[1];
723 }
724 }
725
726 return minW < Double.MAX_VALUE
727 ? new double [] { minW, maxW }
728 : null;
729 }
730
731 /**
732 * Interpolate W and Q values at a given km.
733 */
734 public double [][] interpolateWQ(double km) {
735 return interpolateWQ(km, null);
736 }
737
738 /**
739 * Interpolate W and Q values at a given km.
740 *
741 * @param errors where to store errors.
742 *
743 * @return double double array, first index Ws, second Qs.
744 */
745 public double [][] interpolateWQ(double km, Calculation errors) {
746 return interpolateWQ(km, DEFAULT_Q_STEPS, errors);
747 }
748
749
750 /**
751 * Interpolate W and Q values at a given km.
752 */
753 public double [][] interpolateWQ(double km, int steps, Calculation errors) {
754
755 int rowIndex = Collections.binarySearch(rows, new Row(km));
756
757 if (rowIndex >= 0) { // direct row match
758 Row row = rows.get(rowIndex);
759 return row.interpolateWQ(steps, this, errors);
760 }
761
762 rowIndex = -rowIndex -1;
763
764 if (rowIndex < 1 || rowIndex >= rows.size()) {
765 // do not extrapolate
766 if (errors != null) {
767 errors.addProblem(km, "km.not.found");
768 }
769 return new double[2][0];
770 }
771
772 Row r1 = rows.get(rowIndex-1);
773 Row r2 = rows.get(rowIndex);
774
775 return r1.interpolateWQ(r2, km, steps, this, errors);
776 }
777
778 public boolean interpolate(
779 double km,
780 double [] out,
781 QPosition qPosition,
782 Function qFunction
783 ) {
784 int R1 = rows.size()-1;
785
786 out[1] = qFunction.value(getQ(qPosition, km));
787
788 if (Double.isNaN(out[1])) {
789 return false;
790 }
791
792 QPosition nPosition = new QPosition();
793 if (getQPosition(km, out[1], nPosition) == null) {
794 return false;
795 }
796
797 int rowIndex = Collections.binarySearch(rows, new Row(km));
798
799 if (rowIndex >= 0) {
800 // direct row match
801 out[0] = rows.get(rowIndex).getW(nPosition);
802 return !Double.isNaN(out[0]);
803 }
804
805 rowIndex = -rowIndex -1;
806
807 if (rowIndex < 1 || rowIndex > R1) {
808 // do not extrapolate
809 return false;
810 }
811
812 Row r1 = rows.get(rowIndex-1);
813 Row r2 = rows.get(rowIndex);
814 out[0] = r1.getW(r2, km, nPosition);
815
816 return !Double.isNaN(out[0]);
817 }
818
819
820 /**
821 * Look up interpolation of a Q at given positions.
822 *
823 * @param q the non-interpolated Q value.
824 * @param referenceKm the reference km (e.g. gauge position).
825 * @param kms positions for which to interpolate.
826 * @param ws (output) resulting interpolated ws.
827 * @param qs (output) resulting interpolated qs.
828 * @param errors calculation object to store errors.
829 */
830 public QPosition interpolate(
831 double q,
832 double referenceKm,
833 double [] kms,
834 double [] ws,
835 double [] qs,
836 Calculation errors
837 ) {
838 return interpolate(
839 q, referenceKm, kms, ws, qs, 0, kms.length, errors);
840 }
841
842 public QPosition interpolate(
843 double q,
844 double referenceKm,
845 double [] kms,
846 double [] ws,
847 double [] qs,
848 int startIndex,
849 int length,
850 Calculation errors
851 ) {
852 QPosition qPosition = getQPosition(referenceKm, q);
853
854 if (qPosition == null) {
855 // we cannot locate q at km
856 Arrays.fill(ws, Double.NaN);
857 Arrays.fill(qs, Double.NaN);
858 if (errors != null) {
859 errors.addProblem(referenceKm, "cannot.find.q", q);
860 }
861 return null;
862 }
863
864 Row kmKey = new Row();
865
866 int R1 = rows.size()-1;
867
868 for (int i = startIndex, end = startIndex+length; i < end; ++i) {
869
870 if (Double.isNaN(qs[i] = getQ(qPosition, kms[i]))) {
871 if (errors != null) {
872 errors.addProblem(kms[i], "cannot.find.q", q);
873 }
874 ws[i] = Double.NaN;
875 continue;
876 }
877
878 kmKey.km = kms[i];
879 int rowIndex = Collections.binarySearch(rows, kmKey);
880
881 if (rowIndex >= 0) {
882 // direct row match
883 if (Double.isNaN(ws[i] = rows.get(rowIndex).getW(qPosition))
884 && errors != null) {
885 errors.addProblem(kms[i], "cannot.find.w.for.q", q);
886 }
887 continue;
888 }
889
890 rowIndex = -rowIndex -1;
891
892 if (rowIndex < 1 || rowIndex > R1) {
893 // do not extrapolate
894 if (errors != null) {
895 errors.addProblem(kms[i], "km.not.found");
896 }
897 ws[i] = Double.NaN;
898 continue;
899 }
900 Row r1 = rows.get(rowIndex-1);
901 Row r2 = rows.get(rowIndex);
902
903 if (Double.isNaN(ws[i] = r1.getW(r2, kms[i], qPosition))
904 && errors != null) {
905 errors.addProblem(kms[i], "cannot.find.w.for.q", q);
906 }
907 }
908
909 return qPosition;
910 }
911
912 /**
913 * Linearly interpolate w at a km at a column of two rows.
914 *
915 * @param km position for which to interpolate.
916 * @param row1 first row.
917 * @param row2 second row.
918 * @param col column-index at which to look.
919 *
920 * @return Linearly interpolated w, NaN if one of the given rows was null.
921 */
922 public static double linearW(double km, Row row1, Row row2, int col) {
923 if (row1 == null || row2 == null) {
924 return Double.NaN;
925 }
926
927 return Linear.linear(km,
928 row1.km, row2.km,
929 row1.ws[col], row2.ws[col]);
930 }
931
932 /**
933 * Do interpolation/lookup of W and Q within columns (i.e. ignoring values
934 * of other columns).
935 * @param km position (km) at which to interpolate/lookup.
936 * @return [[q0, q1, .. qx] , [w0, w1, .. wx]] (can contain NaNs)
937 */
938 public double [][] interpolateWQColumnwise(double km) {
939 log.debug("WstValueTable.interpolateWQColumnwise");
940 double [] qs = new double[columns.length];
941 double [] ws = new double[columns.length];
942
943 // Find out row from where we will start searching.
944 int rowIndex = Collections.binarySearch(rows, new Row(km));
945
946 if (rowIndex < 0) {
947 rowIndex = -rowIndex -1;
948 }
949
950 // TODO Beyond definition, we could stop more clever.
951 if (rowIndex >= rows.size()) {
952 rowIndex = rows.size() -1;
953 }
954
955 Row startRow = rows.get(rowIndex);
956
957 for (int col = 0; col < columns.length; col++) {
958 qs[col] = columns[col].getQRangeTree().findQ(km);
959 if (startRow.km == km && startRow.ws[col] != Double.NaN) {
960 // Great. W is defined at km.
961 ws[col] = startRow.ws[col];
962 continue;
963 }
964
965 // Search neighbouring rows that define w at this col.
966 Row rowBefore = null;
967 Row rowAfter = null;
968 for (int before = rowIndex -1; before >= 0; before--) {
969 if (!Double.isNaN(rows.get(before).ws[col])) {
970 rowBefore = rows.get(before);
971 break;
972 }
973 }
974 if (rowBefore != null) {
975 for (int after = rowIndex, R = rows.size(); after < R; after++) {
976 if (!Double.isNaN(rows.get(after).ws[col])) {
977 rowAfter = rows.get(after);
978 break;
979 }
980 }
981 }
982
983 ws[col] = linearW(km, rowBefore, rowAfter, col);
984 }
985
986 return new double [][] {qs, ws};
987 }
988
989 public double [] findQsForW(double km, double w) {
990
991 int rowIndex = Collections.binarySearch(rows, new Row(km));
992
993 if (rowIndex >= 0) {
994 return rows.get(rowIndex).findQsForW(w, this);
995 }
996
997 rowIndex = -rowIndex - 1;
998
999 if (rowIndex < 1 || rowIndex >= rows.size()) {
1000 // Do not extrapolate.
1001 return new double[0];
1002 }
1003
1004 // Needs bilinear interpolation.
1005 Row r1 = rows.get(rowIndex-1);
1006 Row r2 = rows.get(rowIndex);
1007
1008 return r1.findQsForW(r2, w, km, this);
1009 }
1010
1011 protected SplineFunction createSpline(double km, Calculation errors) {
1012
1013 int rowIndex = Collections.binarySearch(rows, new Row(km));
1014
1015 if (rowIndex >= 0) {
1016 SplineFunction sf = rows.get(rowIndex).createSpline(this, errors);
1017 if (sf == null && errors != null) {
1018 errors.addProblem(km, "cannot.create.wq.relation");
1019 }
1020 return sf;
1021 }
1022
1023 rowIndex = -rowIndex - 1;
1024
1025 if (rowIndex < 1 || rowIndex >= rows.size()) {
1026 // Do not extrapolate.
1027 if (errors != null) {
1028 errors.addProblem(km, "km.not.found");
1029 }
1030 return null;
1031 }
1032
1033 // Needs bilinear interpolation.
1034 Row r1 = rows.get(rowIndex-1);
1035 Row r2 = rows.get(rowIndex);
1036
1037 SplineFunction sf = r1.createSpline(r2, km, this, errors);
1038 if (sf == null && errors != null) {
1039 errors.addProblem(km, "cannot.create.wq.relation");
1040 }
1041
1042 return sf;
1043 }
1044
1045 /** 'Bezugslinienverfahren' */
1046 public double [][] relateWs(
1047 double km1,
1048 double km2,
1049 Calculation errors
1050 ) {
1051 return relateWs(km1, km2, RELATE_WS_SAMPLES, errors);
1052 }
1053
1054 private static class ErrorHandler {
1055
1056 boolean hasErrors;
1057 Calculation errors;
1058
1059 ErrorHandler(Calculation errors) {
1060 this.errors = errors;
1061 }
1062
1063 void error(double km, String key, Object ... args) {
1064 if (errors != null && !hasErrors) {
1065 hasErrors = true;
1066 errors.addProblem(km, key, args);
1067 }
1068 }
1069 } // class ErrorHandler
1070
1071
1072 /* TODO: Add optimized methods of relateWs to relate one
1073 * start km to many end kms. The index generation/spline stuff for
1074 * the start km is always the same.
1075 */
1076 public double [][] relateWs(
1077 double km1,
1078 double km2,
1079 int numSamples,
1080 Calculation errors
1081 ) {
1082 SplineFunction sf1 = createSpline(km1, errors);
1083 if (sf1 == null) {
1084 return new double[2][0];
1085 }
1086
1087 SplineFunction sf2 = createSpline(km2, errors);
1088 if (sf2 == null) {
1089 return new double[2][0];
1090 }
1091
1092 PolynomialSplineFunction iQ1 = sf1.createIndexQRelation();
1093 if (iQ1 == null) {
1094 if (errors != null) {
1095 errors.addProblem(km1, "cannot.create.index.q.relation");
1096 }
1097 return new double[2][0];
1098 }
1099
1100 PolynomialSplineFunction iQ2 = sf2.createIndexQRelation();
1101 if (iQ2 == null) {
1102 if (errors != null) {
1103 errors.addProblem(km2, "cannot.create.index.q.relation");
1104 }
1105 return new double[2][0];
1106 }
1107
1108 int N = Math.min(sf1.splineQs.length, sf2.splineQs.length);
1109 double stepWidth = N/(double)numSamples;
1110
1111 PolynomialSplineFunction qW1 = sf1.spline;
1112 PolynomialSplineFunction qW2 = sf2.spline;
1113
1114 TDoubleArrayList ws1 = new TDoubleArrayList(numSamples);
1115 TDoubleArrayList ws2 = new TDoubleArrayList(numSamples);
1116 TDoubleArrayList qs1 = new TDoubleArrayList(numSamples);
1117 TDoubleArrayList qs2 = new TDoubleArrayList(numSamples);
1118
1119 ErrorHandler err = new ErrorHandler(errors);
1120
1121 int i = 0;
1122 for (double p = 0d; p <= N-1; p += stepWidth, ++i) {
1123
1124 double q1;
1125 try {
1126 q1 = iQ1.value(p);
1127 }
1128 catch (ArgumentOutsideDomainException aode) {
1129 err.error(km1, "w.w.qkm1.failed", p);
1130 continue;
1131 }
1132
1133 double w1;
1134 try {
1135 w1 = qW1.value(q1);
1136 }
1137 catch (ArgumentOutsideDomainException aode) {
1138 err.error(km1, "w.w.wkm1.failed", q1, p);
1139 continue;
1140 }
1141
1142 double q2;
1143 try {
1144 q2 = iQ2.value(p);
1145 }
1146 catch (ArgumentOutsideDomainException aode) {
1147 err.error(km2, "w.w.qkm2.failed", p);
1148 continue;
1149 }
1150
1151 double w2;
1152 try {
1153 w2 = qW2.value(q2);
1154 }
1155 catch (ArgumentOutsideDomainException aode) {
1156 err.error(km2, "w.w.wkm2.failed", q2, p);
1157 continue;
1158 }
1159
1160 ws1.add(w1);
1161 ws2.add(w2);
1162 qs1.add(q1);
1163 qs2.add(q2);
1164 }
1165
1166 return new double [][] {
1167 ws1.toNativeArray(),
1168 qs1.toNativeArray(),
1169 ws2.toNativeArray(),
1170 qs2.toNativeArray() };
1171 }
1172
1173 public QPosition getQPosition(double km, double q) {
1174 return getQPosition(km, q, new QPosition());
1175 }
1176
1177 public QPosition getQPosition(double km, double q, QPosition qPosition) {
1178
1179 if (columns.length == 0) {
1180 return null;
1181 }
1182
1183 double qLast = columns[0].getQRangeTree().findQ(km);
1184
1185 if (Math.abs(qLast - q) < 0.00001) {
1186 return qPosition.set(0, 1d);
1187 }
1188
1189 for (int i = 1; i < columns.length; ++i) {
1190 double qCurrent = columns[i].getQRangeTree().findQ(km);
1191 if (Math.abs(qCurrent - q) < 0.00001) {
1192 return qPosition.set(i, 1d);
1193 }
1194
1195 double qMin, qMax;
1196 if (qLast < qCurrent) { qMin = qLast; qMax = qCurrent; }
1197 else { qMin = qCurrent; qMax = qLast; }
1198
1199 if (q > qMin && q < qMax) {
1200 double weight = Linear.factor(q, qLast, qCurrent);
1201 return qPosition.set(i, weight);
1202 }
1203 qLast = qCurrent;
1204 }
1205
1206 return null;
1207 }
1208
1209 public double getQIndex(int index, double km) {
1210 return columns[index].getQRangeTree().findQ(km);
1211 }
1212
1213 public double getQ(QPosition qPosition, double km) {
1214 int index = qPosition.index;
1215 double weight = qPosition.weight;
1216
1217 if (weight == 1d) {
1218 return columns[index].getQRangeTree().findQ(km);
1219 }
1220 double q1 = columns[index-1].getQRangeTree().findQ(km);
1221 double q2 = columns[index ].getQRangeTree().findQ(km);
1222 return Linear.weight(weight, q1, q2);
1223 }
1224
1225 public double [][] interpolateTabulated(double km) {
1226 return interpolateTabulated(km, new double[2][columns.length]);
1227 }
1228
1229 public double [][] interpolateTabulated(double km, double [][] result) {
1230
1231 int rowIndex = Collections.binarySearch(rows, new Row(km));
1232
1233 if (rowIndex >= 0) {
1234 // Direct hit -> copy ws.
1235 Row row = rows.get(rowIndex);
1236 System.arraycopy(
1237 row.ws, 0, result[0], 0,
1238 Math.min(row.ws.length, result[0].length));
1239 }
1240 else {
1241 rowIndex = -rowIndex -1;
1242 if (rowIndex < 1 || rowIndex >= rows.size()) {
1243 // Out of bounds.
1244 return null;
1245 }
1246 // Interpolate ws.
1247 Row r1 = rows.get(rowIndex-1);
1248 Row r2 = rows.get(rowIndex);
1249 double factor = Linear.factor(km, r1.km, r2.km);
1250 Linear.weight(factor, r1.ws, r2.ws, result[0]);
1251 }
1252
1253 double [] qs = result[1];
1254 for (int i = Math.min(qs.length, columns.length)-1; i >= 0; --i) {
1255 qs[i] = columns[i].getQRangeTree().findQ(km);
1256 }
1257 return result;
1258 }
1259
1260
1261 /** Find ranges that are between km1 and km2 (inclusive?) */
1262 public List<Range> findSegments(double km1, double km2) {
1263 return columns.length != 0
1264 ? columns[columns.length-1].getQRangeTree().findSegments(km1, km2)
1265 : Collections.<Range>emptyList();
1266 }
1267 }
1268 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
http://dive4elements.wald.intevation.org