Mercurial > dive4elements > river

comparison artifacts/src/main/java/org/dive4elements/river/artifacts/model/WstValueTable.java @ 9286:bcbae86ce7b3

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