Mercurial > dive4elements > river

comparison flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java @ 3812:f788d2d901d6

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