Mercurial > dive4elements > river

comparison artifacts/src/main/java/org/dive4elements/river/artifacts/model/WstValueTable.java @ 8704:93a31cfb18c0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
(issue1787) Sort WQ data locally over Q bevor lookup of Qs for given W.
author Tom Gottfried <tom@intevation.de>
date Thu, 23 Apr 2015 19:13:33 +0200
parents b9e5fa7f7a28
children 5e38e2924c07
comparison
equal deleted inserted replaced
8703:e4f9e2316e92 8704:93a31cfb18c0
197 } 197 }
198 198
199 public Row(double km, double [] ws) { 199 public Row(double km, double [] ws) {
200 this(km); 200 this(km);
201 this.ws = ws; 201 this.ws = ws;
202 }
203
204 /**
205 * Sort Qs and Ws for this Row over Q.
206 */
207 private double[][] getSortedWQ(
208 WstValueTable table,
209 Calculation errors
210 ) {
211 int W = ws.length;
212
213 if (W < 1) {
214 if (errors != null) {
215 errors.addProblem(km, "no.ws.found");
216 }
217 return new double[][] {{Double.NaN}, {Double.NaN}};
218 }
219
220 double [] sortedWs = ws;
221 double [] sortedQs = new double[W];
222
223 for (int i=0; i < W; ++i) {
224 double q = table.getQIndex(i, km);
225 if (Double.isNaN(q) && errors != null) {
226 errors.addProblem(
227 km, "no.q.found.in.column", i+1);
228 }
229 sortedQs[i] = q;
230 }
231
232 DoubleUtil.sortByFirst(sortedQs, sortedWs);
233
234 return new double[][] { sortedWs, sortedQs };
235 }
236
237 /**
238 * Return an array of Qs and Ws for the given km between
239 * this Row and another, sorted over Q.
240 */
241 private double[][] getSortedWQ(
242 Row other,
243 double km,
244 WstValueTable table,
245 Calculation errors
246 ) {
247 int W = Math.min(ws.length, other.ws.length);
248
249 if (W < 1) {
250 if (errors != null) {
251 errors.addProblem("no.ws.found");
252 }
253 return new double[][] {{Double.NaN}, {Double.NaN}};
254 }
255
256 double factor = Linear.factor(km, this.km, other.km);
257
258 double [] sortedQs = new double[W];
259 double [] sortedWs = new double[W];
260
261 for (int i = 0; i < W; ++i) {
262 double wws = Linear.weight(factor, ws[i], other.ws[i]);
263 double wqs = table.getQIndex(i, km);
264
265 if (Double.isNaN(wws) || Double.isNaN(wqs)) {
266 if (errors != null) {
267 errors.addProblem(km, "cannot.find.w.or.q");
268 }
269 }
270
271 sortedWs[i] = wws;
272 sortedQs[i] = wqs;
273 }
274
275 DoubleUtil.sortByFirst(sortedQs, sortedWs);
276
277 return new double[][] { sortedWs, sortedQs };
278 }
279
280 /**
281 * Find Qs matching w in an array of Qs and Ws sorted over Q.
282 */
283 private double[] findQsForW(double w, double[][] sortedWQ) {
284 int W = sortedWQ[0].length;
285
286 double[] sortedQs = sortedWQ[1];
287 double[] sortedWs = sortedWQ[0];
288
289 TDoubleArrayList qs = new TDoubleArrayList();
290
291 if (W > 0 && Math.abs(sortedWs[0]-w) < W_EPSILON) {
292 double q = sortedQs[0];
293 if (!Double.isNaN(q)) {
294 qs.add(q);
295 }
296 }
297
298 for (int i = 1; i < W; ++i) {
299 double w2 = sortedWs[i];
300 if (Double.isNaN(w2)) {
301 continue;
302 }
303 if (Math.abs(w2-w) < W_EPSILON) {
304 double q = sortedQs[i];
305 if (!Double.isNaN(q)) {
306 qs.add(q);
307 }
308 continue;
309 }
310 double w1 = sortedWs[i-1];
311 if (Double.isNaN(w1)) {
312 continue;
313 }
314
315 if (w < Math.min(w1, w2) || w > Math.max(w1, w2)) {
316 continue;
317 }
318
319 double q1 = sortedQs[i-1];
320 double q2 = sortedQs[i];
321 if (Double.isNaN(q1) || Double.isNaN(q2)) {
322 continue;
323 }
324
325 double q = Linear.linear(w, w1, w2, q1, q2);
326 qs.add(q);
327 }
328
329 return qs.toNativeArray();
202 } 330 }
203 331
204 /** 332 /**
205 * Compare according to place of measurement (km). 333 * Compare according to place of measurement (km).
206 */ 334 */
254 382
255 public SplineFunction createSpline( 383 public SplineFunction createSpline(
256 WstValueTable table, 384 WstValueTable table,
257 Calculation errors 385 Calculation errors
258 ) { 386 ) {
259 int W = ws.length; 387 double[][] sortedWQ = getSortedWQ(table, errors);
260
261 if (W < 1) {
262 if (errors != null) {
263 errors.addProblem(km, "no.ws.found");
264 }
265 return null;
266 }
267
268 double [] splineQs = new double[W];
269 double [] splineWs = new double[W];
270
271 for (int i = 0; i < W; ++i) {
272 double sq = table.getQIndex(i, km);
273 if (Double.isNaN(sq) && errors != null) {
274 errors.addProblem(
275 km, "no.q.found.in.column", (i+1));
276 }
277 splineQs[i] = sq;
278 splineWs[i] = ws[i];
279 }
280
281 DoubleUtil.sortByFirst(splineQs, splineWs);
282 388
283 try { 389 try {
284 SplineInterpolator interpolator = new SplineInterpolator(); 390 SplineInterpolator interpolator = new SplineInterpolator();
285 PolynomialSplineFunction spline = 391 PolynomialSplineFunction spline =
286 interpolator.interpolate(splineQs, splineWs); 392 interpolator.interpolate(sortedWQ[1], sortedWQ[0]);
287 393
288 return new SplineFunction(spline, splineQs, ws); 394 return new SplineFunction(spline, sortedWQ[1], ws);
289 } 395 }
290 catch (MathIllegalArgumentException miae) { 396 catch (MathIllegalArgumentException miae) {
291 if (errors != null) { 397 if (errors != null) {
292 errors.addProblem(km, "spline.creation.failed"); 398 errors.addProblem(km, "spline.creation.failed");
293 } 399 }
300 Row other, 406 Row other,
301 double km, 407 double km,
302 WstValueTable table, 408 WstValueTable table,
303 Calculation errors 409 Calculation errors
304 ) { 410 ) {
305 int W = Math.min(ws.length, other.ws.length); 411 double[][] sortedWQ = getSortedWQ(other, km, table, errors);
306
307 if (W < 1) {
308 if (errors != null) {
309 errors.addProblem("no.ws.found");
310 }
311 return null;
312 }
313
314 double factor = Linear.factor(km, this.km, other.km);
315
316 double [] splineQs = new double[W];
317 double [] splineWs = new double[W];
318
319 for (int i = 0; i < W; ++i) {
320 double wws = Linear.weight(factor, ws[i], other.ws[i]);
321 double wqs = Linear.weight(
322 factor,
323 table.getQIndex(i, km),
324 table.getQIndex(i, other.km));
325
326 if (Double.isNaN(wws) || Double.isNaN(wqs)) {
327 if (errors != null) {
328 errors.addProblem(km, "cannot.find.w.or.q");
329 }
330 }
331
332 splineWs[i] = wws;
333 splineQs[i] = wqs;
334 }
335
336 DoubleUtil.sortByFirst(splineQs, splineWs);
337 412
338 SplineInterpolator interpolator = new SplineInterpolator(); 413 SplineInterpolator interpolator = new SplineInterpolator();
339 414
340 try { 415 try {
341 PolynomialSplineFunction spline = 416 PolynomialSplineFunction spline =
342 interpolator.interpolate(splineQs, splineWs); 417 interpolator.interpolate(sortedWQ[1], sortedWQ[0]);
343 418
344 return new SplineFunction(spline, splineQs, splineWs); 419 return new SplineFunction(spline, sortedWQ[1], sortedWQ[0]);
345 } 420 }
346 catch (MathIllegalArgumentException miae) { 421 catch (MathIllegalArgumentException miae) {
347 if (errors != null) { 422 if (errors != null) {
348 errors.addProblem(km, "spline.creation.failed"); 423 errors.addProblem(km, "spline.creation.failed");
349 } 424 }
412 } 487 }
413 488
414 return Linear.weight(kmWeight, tw, ow); 489 return Linear.weight(kmWeight, tw, ow);
415 } 490 }
416 491
417 public double [] findQsForW(double w, WstValueTable table) { 492 public double [] findQsForW(
418 493 double w,
419 TDoubleArrayList qs = new TDoubleArrayList(); 494 WstValueTable table,
420 495 Calculation errors
421 if (ws.length > 0 && Math.abs(ws[0]-w) < W_EPSILON) { 496 ) {
422 double q = table.getQIndex(0, km); 497 log.debug("Find Qs for given W at tabulated km " + km);
423 if (!Double.isNaN(q)) { 498 return findQsForW(w, getSortedWQ(table, errors));
424 qs.add(q);
425 }
426 }
427
428 for (int i = 1; i < ws.length; ++i) {
429 double w2 = ws[i];
430 if (Double.isNaN(w2)) {
431 continue;
432 }
433 if (Math.abs(w2-w) < W_EPSILON) {
434 double q = table.getQIndex(i, km);
435 if (!Double.isNaN(q)) {
436 qs.add(q);
437 }
438 continue;
439 }
440 double w1 = ws[i-1];
441 if (Double.isNaN(w1)) {
442 continue;
443 }
444
445 if (w < Math.min(w1, w2) || w > Math.max(w1, w2)) {
446 continue;
447 }
448
449 double q1 = table.getQIndex(i-1, km);
450 double q2 = table.getQIndex(i, km);
451 if (Double.isNaN(q1) || Double.isNaN(q2)) {
452 continue;
453 }
454
455 double q = Linear.linear(w, w1, w2, q1, q2);
456 qs.add(q);
457 }
458
459 return qs.toNativeArray();
460 } 499 }
461 500
462 public double [] findQsForW( 501 public double [] findQsForW(
463 Row other, 502 Row other,
464 double w, 503 double w,
465 double km, 504 double km,
466 WstValueTable table 505 WstValueTable table,
506 Calculation errors
467 ) { 507 ) {
468 TDoubleArrayList qs = new TDoubleArrayList(); 508 log.debug("Find Qs for given W at non-tabulated km " + km);
469 509 return findQsForW(w, getSortedWQ(other, km, table, errors));
470 double factor = Linear.factor(km, this.km, other.km);
471
472 if (ws.length > 0) {
473 double wt = Linear.weight(factor, ws[0], other.ws[0]);
474 if (Math.abs(wt-w) < W_EPSILON) {
475 double q = table.getQIndex(0, km);
476 if (!Double.isNaN(q)) {
477 qs.add(q);
478 }
479 }
480 }
481
482 for (int i = 1; i < ws.length; ++i) {
483 double w2 = Linear.weight(factor, ws[i], other.ws[i]);
484 if (Double.isNaN(w2)) {
485 continue;
486 }
487 if (Math.abs(w2-w) < W_EPSILON) {
488 double q = table.getQIndex(i, km);
489 if (!Double.isNaN(q)) {
490 qs.add(q);
491 }
492 continue;
493 }
494
495 double w1 = Linear.weight(factor, ws[i-1], other.ws[i-1]);
496 if (Double.isNaN(w1)) {
497 continue;
498 }
499
500 if (w < Math.min(w1, w2) || w > Math.max(w1, w2)) {
501 continue;
502 }
503
504 double q1 = table.getQIndex(i-1, km);
505 double q2 = table.getQIndex(i, km);
506 if (Double.isNaN(q1) || Double.isNaN(q2)) {
507 continue;
508 }
509
510 double q = Linear.linear(w, w1, w2, q1, q2);
511 qs.add(q);
512 }
513
514 return qs.toNativeArray();
515 } 510 }
516 511
517 public double [] getMinMaxW(double [] result) { 512 public double [] getMinMaxW(double [] result) {
518 double minW = Double.MAX_VALUE; 513 double minW = Double.MAX_VALUE;
519 double maxW = -Double.MAX_VALUE; 514 double maxW = -Double.MAX_VALUE;
1006 } 1001 }
1007 1002
1008 return new double [][] {qs, ws}; 1003 return new double [][] {qs, ws};
1009 } 1004 }
1010 1005
1011 public double [] findQsForW(double km, double w) { 1006 public double [] findQsForW(double km, double w, Calculation errors) {
1012 1007
1013 int rowIndex = Collections.binarySearch(rows, new Row(km)); 1008 int rowIndex = Collections.binarySearch(rows, new Row(km));
1014 1009
1015 if (rowIndex >= 0) { 1010 if (rowIndex >= 0) {
1016 return rows.get(rowIndex).findQsForW(w, this); 1011 return rows.get(rowIndex).findQsForW(w, this, errors);
1017 } 1012 }
1018 1013
1019 rowIndex = -rowIndex - 1; 1014 rowIndex = -rowIndex - 1;
1020 1015
1021 if (rowIndex < 1 || rowIndex >= rows.size()) { 1016 if (rowIndex < 1 || rowIndex >= rows.size()) {
1025 1020
1026 // Needs bilinear interpolation. 1021 // Needs bilinear interpolation.
1027 Row r1 = rows.get(rowIndex-1); 1022 Row r1 = rows.get(rowIndex-1);
1028 Row r2 = rows.get(rowIndex); 1023 Row r2 = rows.get(rowIndex);
1029 1024
1030 return r1.findQsForW(r2, w, km, this); 1025 return r1.findQsForW(r2, w, km, this, errors);
1031 } 1026 }
1032 1027
1033 protected SplineFunction createSpline(double km, Calculation errors) { 1028 protected SplineFunction createSpline(double km, Calculation errors) {
1034 1029
1035 int rowIndex = Collections.binarySearch(rows, new Row(km)); 1030 int rowIndex = Collections.binarySearch(rows, new Row(km));
http://dive4elements.wald.intevation.org