Mercurial > dive4elements > river

comparison flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java @ 2424:092e519ff461

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merged flys-artifacts/2.6.1
author Thomas Arendsen Hein <thomas@intevation.de>
date Fri, 28 Sep 2012 12:14:26 +0200
parents ac528b883b47
children 8490faba00e7
comparison
equal deleted inserted replaced
2392:8112ec686a9a 2424:092e519ff461
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 [] getMinMaxW(double km) {
611 return getMinMaxW(km, new double [2]);
612
613 }
614 public double [] getMinMaxW(double km, double [] result) {
615 int rowIndex = Collections.binarySearch(rows, new Row(km));
616
617 if (rowIndex >= 0) {
618 return rows.get(rowIndex).getMinMaxW(result);
619 }
620
621 rowIndex = -rowIndex -1;
622
623 if (rowIndex < 1 || rowIndex >= rows.size()) {
624 // do not extrapolate
625 return null;
626 }
627
628 Row r1 = rows.get(rowIndex-1);
629 Row r2 = rows.get(rowIndex);
630
631 return r1.getMinMaxW(r2, km, result);
632 }
633
634 public double [] getMinMaxW(double from, double to, double step) {
635 double [] result = new double[2];
636
637 double minW = Double.MAX_VALUE;
638 double maxW = -Double.MAX_VALUE;
639
640 if (from > to) {
641 double tmp = from;
642 from = to;
643 to = tmp;
644 }
645
646 step = Math.max(Math.abs(step), 0.0001);
647
648 double d = from;
649 for (; d <= to; d += step) {
650 if (getMinMaxW(d, result) != null) {
651 if (result[0] < minW) minW = result[0];
652 if (result[1] > maxW) maxW = result[1];
653 }
654 }
655
656 if (d != to) {
657 if (getMinMaxW(to, result) != null) {
658 if (result[0] < minW) minW = result[0];
659 if (result[1] > maxW) maxW = result[1];
660 }
661 }
662
663 return minW < Double.MAX_VALUE
664 ? new double [] { minW, maxW }
665 : null;
666 }
667
668 /**
669 * Interpolate W and Q values at a given km.
670 */
671 public double [][] interpolateWQ(double km) {
672 return interpolateWQ(km, null);
673 }
674
675 /**
676 * Interpolate W and Q values at a given km.
677 */
678 public double [][] interpolateWQ(double km, Calculation errors) {
679 return interpolateWQ(km, DEFAULT_Q_STEPS, errors);
680 }
681
682 /**
683 * Interpolate W and Q values at a given km.
684 */
685 public double [][] interpolateWQ(double km, int steps, Calculation errors) {
686
687 int rowIndex = Collections.binarySearch(rows, new Row(km));
688
689 if (rowIndex >= 0) { // direct row match
690 Row row = rows.get(rowIndex);
691 return row.interpolateWQ(steps, this, errors);
692 }
693
694 rowIndex = -rowIndex -1;
695
696 if (rowIndex < 1 || rowIndex >= rows.size()) {
697 // do not extrapolate
698 if (errors != null) {
699 errors.addProblem(km, "km.not.found");
700 }
701 return new double[2][0];
702 }
703
704 Row r1 = rows.get(rowIndex-1);
705 Row r2 = rows.get(rowIndex);
706
707 return r1.interpolateWQ(r2, km, steps, this, errors);
708 }
709
710 public boolean interpolate(
711 double km,
712 double [] out,
713 QPosition qPosition,
714 Function qFunction
715 ) {
716 int R1 = rows.size()-1;
717
718 out[1] = qFunction.value(getQ(qPosition, km));
719
720 if (Double.isNaN(out[1])) {
721 return false;
722 }
723
724 QPosition nPosition = new QPosition();
725 if (getQPosition(km, out[1], nPosition) == null) {
726 return false;
727 }
728
729 int rowIndex = Collections.binarySearch(rows, new Row(km));
730
731 if (rowIndex >= 0) {
732 // direct row match
733 out[0] = rows.get(rowIndex).getW(nPosition);
734 return !Double.isNaN(out[0]);
735 }
736
737 rowIndex = -rowIndex -1;
738
739 if (rowIndex < 1 || rowIndex > R1) {
740 // do not extrapolate
741 return false;
742 }
743
744 Row r1 = rows.get(rowIndex-1);
745 Row r2 = rows.get(rowIndex);
746 out[0] = r1.getW(r2, km, nPosition);
747
748 return !Double.isNaN(out[0]);
749 }
750
751
752 /**
753 * Look up interpolation of a Q at given positions.
754 *
755 * @param q the non-interpolated Q value.
756 * @param referenceKm the reference km (e.g. gauge position).
757 * @param kms positions for which to interpolate.
758 * @param ws (output) resulting interpolated ws.
759 * @param qs (output) resulting interpolated qs.
760 * @param errors calculation object to store errors.
761 */
762 public QPosition interpolate(
763 double q,
764 double referenceKm,
765 double [] kms,
766 double [] ws,
767 double [] qs,
768 Calculation errors
769 ) {
770 return interpolate(
771 q, referenceKm, kms, ws, qs, 0, kms.length, errors);
772 }
773
774 public QPosition interpolate(
775 double q,
776 double referenceKm,
777 double [] kms,
778 double [] ws,
779 double [] qs,
780 int startIndex,
781 int length,
782 Calculation errors
783 ) {
784 QPosition qPosition = getQPosition(referenceKm, q);
785
786 if (qPosition == null) {
787 // we cannot locate q at km
788 Arrays.fill(ws, Double.NaN);
789 Arrays.fill(qs, Double.NaN);
790 if (errors != null) {
791 errors.addProblem(referenceKm, "cannot.find.q", q);
792 }
793 return null;
794 }
795
796 Row kmKey = new Row();
797
798 int R1 = rows.size()-1;
799
800 for (int i = startIndex, end = startIndex+length; i < end; ++i) {
801
802 if (Double.isNaN(qs[i] = getQ(qPosition, kms[i]))) {
803 if (errors != null) {
804 errors.addProblem(kms[i], "cannot.find.q", q);
805 }
806 ws[i] = Double.NaN;
807 continue;
808 }
809
810 kmKey.km = kms[i];
811 int rowIndex = Collections.binarySearch(rows, kmKey);
812
813 if (rowIndex >= 0) {
814 // direct row match
815 if (Double.isNaN(ws[i] = rows.get(rowIndex).getW(qPosition))
816 && errors != null) {
817 errors.addProblem(kms[i], "cannot.find.w.for.q", q);
818 }
819 continue;
820 }
821
822 rowIndex = -rowIndex -1;
823
824 if (rowIndex < 1 || rowIndex > R1) {
825 // do not extrapolate
826 if (errors != null) {
827 errors.addProblem(kms[i], "km.not.found");
828 }
829 ws[i] = Double.NaN;
830 continue;
831 }
832 Row r1 = rows.get(rowIndex-1);
833 Row r2 = rows.get(rowIndex);
834
835 if (Double.isNaN(ws[i] = r1.getW(r2, kms[i], qPosition))
836 && errors != null) {
837 errors.addProblem(kms[i], "cannot.find.w.for.q", q);
838 }
839 }
840
841 return qPosition;
842 }
843
844 /**
845 * Linearly interpolate w at a km at a column of two rows.
846 *
847 * @param km position for which to interpolate.
848 * @param row1 first row.
849 * @param row2 second row.
850 * @param col column-index at which to look.
851 *
852 * @return Linearly interpolated w, NaN if one of the given rows was null.
853 */
854 public static double linearW(double km, Row row1, Row row2, int col) {
855 if (row1 == null || row2 == null) {
856 return Double.NaN;
857 }
858
859 return Linear.linear(km,
860 row1.km, row2.km,
861 row1.ws[col], row2.ws[col]);
862 }
863
864 /**
865 * Do interpolation/lookup of W and Q within columns (i.e. ignoring values
866 * of other columns).
867 * @param km position (km) at which to interpolate/lookup.
868 * @return [[q0, q1, .. qx] , [w0, w1, .. wx]] (can contain NaNs)
869 */
870 public double [][] interpolateWQColumnwise(double km) {
871 log.debug("WstValueTable.interpolateWQColumnwise");
872 double [] qs = new double[columns.length];
873 double [] ws = new double[columns.length];
874
875 // Find out row from where we will start searching.
876 int rowIndex = Collections.binarySearch(rows, new Row(km));
877
878 if (rowIndex < 0) {
879 rowIndex = -rowIndex -1;
880 }
881
882 // TODO Beyond definition, we could stop more clever.
883 if (rowIndex >= rows.size()) {
884 rowIndex = rows.size() -1;
885 }
886
887 Row startRow = rows.get(rowIndex);
888
889 for (int col = 0; col < columns.length; col++) {
890 qs[col] = columns[col].getQRangeTree().findQ(km);
891 if (startRow.km == km && startRow.ws[col] != Double.NaN) {
892 // Great. W is defined at km.
893 ws[col] = startRow.ws[col];
894 continue;
895 }
896
897 // Search neighbouring rows that define w at this col.
898 Row rowBefore = null;
899 Row rowAfter = null;
900 for (int before = rowIndex -1; before >= 0; before--) {
901 if (!Double.isNaN(rows.get(before).ws[col])) {
902 rowBefore = rows.get(before);
903 break;
904 }
905 }
906 if (rowBefore != null) {
907 for (int after = rowIndex, R = rows.size(); after < R; after++) {
908 if (!Double.isNaN(rows.get(after).ws[col])) {
909 rowAfter = rows.get(after);
910 break;
911 }
912 }
913 }
914
915 ws[col] = linearW(km, rowBefore, rowAfter, col);
916 }
917
918 return new double [][] {qs, ws};
919 }
920
921 public double [] findQsForW(double km, double w) {
922
923 int rowIndex = Collections.binarySearch(rows, new Row(km));
924
925 if (rowIndex >= 0) {
926 return rows.get(rowIndex).findQsForW(w, this);
927 }
928
929 rowIndex = -rowIndex - 1;
930
931 if (rowIndex < 1 || rowIndex >= rows.size()) {
932 // Do not extrapolate.
933 return new double[0];
934 }
935
936 // Needs bilinear interpolation.
937 Row r1 = rows.get(rowIndex-1);
938 Row r2 = rows.get(rowIndex);
939
940 return r1.findQsForW(r2, w, km, this);
941 }
942
943 protected SplineFunction createSpline(double km, Calculation errors) {
944
945 int rowIndex = Collections.binarySearch(rows, new Row(km));
946
947 if (rowIndex >= 0) {
948 SplineFunction sf = rows.get(rowIndex).createSpline(this, errors);
949 if (sf == null && errors != null) {
950 errors.addProblem(km, "cannot.create.wq.relation");
951 }
952 return sf;
953 }
954
955 rowIndex = -rowIndex - 1;
956
957 if (rowIndex < 1 || rowIndex >= rows.size()) {
958 // Do not extrapolate.
959 if (errors != null) {
960 errors.addProblem(km, "km.not.found");
961 }
962 return null;
963 }
964
965 // Needs bilinear interpolation.
966 Row r1 = rows.get(rowIndex-1);
967 Row r2 = rows.get(rowIndex);
968
969 SplineFunction sf = r1.createSpline(r2, km, this, errors);
970 if (sf == null && errors != null) {
971 errors.addProblem(km, "cannot.create.wq.relation");
972 }
973
974 return sf;
975 }
976
977 /** 'Bezugslinienverfahren' */
978 public double [][] relateWs(
979 double km1,
980 double km2,
981 Calculation errors
982 ) {
983 return relateWs(km1, km2, RELATE_WS_SAMPLES, errors);
984 }
985
986 private static class ErrorHandler {
987
988 boolean hasErrors;
989 Calculation errors;
990
991 ErrorHandler(Calculation errors) {
992 this.errors = errors;
993 }
994
995 void error(double km, String key, Object ... args) {
996 if (errors != null && !hasErrors) {
997 hasErrors = true;
998 errors.addProblem(km, key, args);
999 }
1000 }
1001 } // class ErrorHandler
1002
1003
1004 /* TODO: Add optimized methods of relateWs to relate one
1005 * start km to many end kms. The index generation/spline stuff for
1006 * the start km is always the same.
1007 */
1008 public double [][] relateWs(
1009 double km1,
1010 double km2,
1011 int numSamples,
1012 Calculation errors
1013 ) {
1014 SplineFunction sf1 = createSpline(km1, errors);
1015 if (sf1 == null) {
1016 return new double[2][0];
1017 }
1018
1019 SplineFunction sf2 = createSpline(km2, errors);
1020 if (sf2 == null) {
1021 return new double[2][0];
1022 }
1023
1024 PolynomialSplineFunction iQ1 = sf1.createIndexQRelation();
1025 if (iQ1 == null) {
1026 if (errors != null) {
1027 errors.addProblem(km1, "cannot.create.index.q.relation");
1028 }
1029 return new double[2][0];
1030 }
1031
1032 PolynomialSplineFunction iQ2 = sf2.createIndexQRelation();
1033 if (iQ2 == null) {
1034 if (errors != null) {
1035 errors.addProblem(km2, "cannot.create.index.q.relation");
1036 }
1037 return new double[2][0];
1038 }
1039
1040 int N = Math.min(sf1.splineQs.length, sf2.splineQs.length);
1041 double stepWidth = N/(double)numSamples;
1042
1043 PolynomialSplineFunction qW1 = sf1.spline;
1044 PolynomialSplineFunction qW2 = sf2.spline;
1045
1046 TDoubleArrayList ws1 = new TDoubleArrayList(numSamples);
1047 TDoubleArrayList ws2 = new TDoubleArrayList(numSamples);
1048 TDoubleArrayList qs1 = new TDoubleArrayList(numSamples);
1049 TDoubleArrayList qs2 = new TDoubleArrayList(numSamples);
1050
1051 ErrorHandler err = new ErrorHandler(errors);
1052
1053 int i = 0;
1054 for (double p = 0d; p <= N-1; p += stepWidth, ++i) {
1055
1056 double q1;
1057 try {
1058 q1 = iQ1.value(p);
1059 }
1060 catch (ArgumentOutsideDomainException aode) {
1061 err.error(km1, "w.w.qkm1.failed", p);
1062 continue;
1063 }
1064
1065 double w1;
1066 try {
1067 w1 = qW1.value(q1);
1068 }
1069 catch (ArgumentOutsideDomainException aode) {
1070 err.error(km1, "w.w.wkm1.failed", q1, p);
1071 continue;
1072 }
1073
1074 double q2;
1075 try {
1076 q2 = iQ2.value(p);
1077 }
1078 catch (ArgumentOutsideDomainException aode) {
1079 err.error(km2, "w.w.qkm2.failed", p);
1080 continue;
1081 }
1082
1083 double w2;
1084 try {
1085 w2 = qW2.value(q2);
1086 }
1087 catch (ArgumentOutsideDomainException aode) {
1088 err.error(km2, "w.w.wkm2.failed", q2, p);
1089 continue;
1090 }
1091
1092 ws1.add(w1);
1093 ws2.add(w2);
1094 qs1.add(q1);
1095 qs2.add(q2);
1096 }
1097
1098 return new double [][] {
1099 ws1.toNativeArray(),
1100 qs1.toNativeArray(),
1101 ws2.toNativeArray(),
1102 qs2.toNativeArray() };
1103 }
1104
1105 public QPosition getQPosition(double km, double q) {
1106 return getQPosition(km, q, new QPosition());
1107 }
1108
1109 public QPosition getQPosition(double km, double q, QPosition qPosition) {
1110
1111 if (columns.length == 0) {
1112 return null;
1113 }
1114
1115 double qLast = columns[0].getQRangeTree().findQ(km);
1116
1117 if (Math.abs(qLast - q) < 0.00001) {
1118 return qPosition.set(0, 1d);
1119 }
1120
1121 for (int i = 1; i < columns.length; ++i) {
1122 double qCurrent = columns[i].getQRangeTree().findQ(km);
1123 if (Math.abs(qCurrent - q) < 0.00001) {
1124 return qPosition.set(i, 1d);
1125 }
1126
1127 double qMin, qMax;
1128 if (qLast < qCurrent) { qMin = qLast; qMax = qCurrent; }
1129 else { qMin = qCurrent; qMax = qLast; }
1130
1131 if (q > qMin && q < qMax) {
1132 double weight = Linear.factor(q, qLast, qCurrent);
1133 return qPosition.set(i, weight);
1134 }
1135 qLast = qCurrent;
1136 }
1137
1138 return null;
1139 }
1140
1141 public double getQIndex(int index, double km) {
1142 return columns[index].getQRangeTree().findQ(km);
1143 }
1144
1145 public double getQ(QPosition qPosition, double km) {
1146 int index = qPosition.index;
1147 double weight = qPosition.weight;
1148
1149 if (weight == 1d) {
1150 return columns[index].getQRangeTree().findQ(km);
1151 }
1152 double q1 = columns[index-1].getQRangeTree().findQ(km);
1153 double q2 = columns[index ].getQRangeTree().findQ(km);
1154 return Linear.weight(weight, q1, q2);
1155 }
1156 }
1157 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
http://dive4elements.wald.intevation.org