Mercurial > dive4elements > river

comparison flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java @ 655:913b52064449

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Refactored version of "Berechnung 4" flys-artifacts/trunk@2053 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Sun, 05 Jun 2011 18:24:46 +0000
parents 44175d4720f8
children c501f27c1f71 3dc61e00385e
comparison
equal deleted inserted replaced
654:bbc966c81809 655:913b52064449
1 package de.intevation.flys.artifacts.model; 1 package de.intevation.flys.artifacts.model;
2 2
3 import java.io.Serializable; 3 import java.io.Serializable;
4 4
5 5 import de.intevation.flys.artifacts.math.Linear;
6 import de.intevation.flys.artifacts.math.LinearRemap; 6 import de.intevation.flys.artifacts.math.Function;
7 7
8 import java.util.Arrays; 8 import java.util.Arrays;
9 import java.util.ArrayList; 9 import java.util.ArrayList;
10 import java.util.List; 10 import java.util.List;
11 import java.util.Collections; 11 import java.util.Collections;
107 double km, 107 double km,
108 double [] iqs, 108 double [] iqs,
109 double [] ows, 109 double [] ows,
110 WstValueTable table 110 WstValueTable table
111 ) { 111 ) {
112 double kmWeight = factor(km, this.km, other.km); 112 double kmWeight = Linear.factor(km, this.km, other.km);
113 113
114 QPosition qPosition = new QPosition(); 114 QPosition qPosition = new QPosition();
115 115
116 for (int i = 0; i < iqs.length; ++i) { 116 for (int i = 0; i < iqs.length; ++i) {
117 if (table.getQPosition(km, iqs[i], qPosition) != null) { 117 if (table.getQPosition(km, iqs[i], qPosition) != null) {
118 double wt = getW(qPosition); 118 double wt = getW(qPosition);
119 double wo = other.getW(qPosition); 119 double wo = other.getW(qPosition);
120 ows[i] = weight(kmWeight, wt, wo); 120 ows[i] = Linear.weight(kmWeight, wt, wo);
121 } 121 }
122 else { 122 else {
123 ows[i] = Double.NaN; 123 ows[i] = Double.NaN;
124 } 124 }
125 } 125 }
135 135
136 if (W < 1) { 136 if (W < 1) {
137 return new double[2][0]; 137 return new double[2][0];
138 } 138 }
139 139
140 double factor = factor(km, this.km, other.km); 140 double factor = Linear.factor(km, this.km, other.km);
141 141
142 double [] splineQ = new double[W]; 142 double [] splineQ = new double[W];
143 double [] splineW = new double[W]; 143 double [] splineW = new double[W];
144 144
145 double minQ = Double.MAX_VALUE; 145 double minQ = Double.MAX_VALUE;
146 double maxQ = -Double.MAX_VALUE; 146 double maxQ = -Double.MAX_VALUE;
147 147
148 for (int i = 0; i < W; ++i) { 148 for (int i = 0; i < W; ++i) {
149 double wws = weight(factor, ws[i], other.ws[i]); 149 double wws = Linear.weight(factor, ws[i], other.ws[i]);
150 double wqs = weight( 150 double wqs = Linear.weight(
151 factor, 151 factor,
152 table.getQIndex(i, km), 152 table.getQIndex(i, km),
153 table.getQIndex(i, other.km)); 153 table.getQIndex(i, other.km));
154 154
155 splineW[i] = wws; 155 splineW[i] = wws;
229 int index = qPosition.index; 229 int index = qPosition.index;
230 double weight = qPosition.weight; 230 double weight = qPosition.weight;
231 231
232 return weight == 1.0 232 return weight == 1.0
233 ? ws[index] 233 ? ws[index]
234 : weight(weight, ws[index-1], ws[index]); 234 : Linear.weight(weight, ws[index-1], ws[index]);
235 } 235 }
236 236
237 public double getW( 237 public double getW(
238 Row other, 238 Row other,
239 double km, 239 double km,
240 QPosition qPosition 240 QPosition qPosition
241 ) { 241 ) {
242 double kmWeight = factor(km, this.km, other.km); 242 double kmWeight = Linear.factor(km, this.km, other.km);
243 243
244 int index = qPosition.index; 244 int index = qPosition.index;
245 double weight = qPosition.weight; 245 double weight = qPosition.weight;
246 246
247 double tw, ow; 247 double tw, ow;
249 if (weight == 1.0) { 249 if (weight == 1.0) {
250 tw = ws[index]; 250 tw = ws[index];
251 ow = other.ws[index]; 251 ow = other.ws[index];
252 } 252 }
253 else { 253 else {
254 tw = weight(weight, ws[index-1], ws[index]); 254 tw = Linear.weight(weight, ws[index-1], ws[index]);
255 ow = weight(weight, other.ws[index-1], other.ws[index]); 255 ow = Linear.weight(weight, other.ws[index-1], other.ws[index]);
256 } 256 }
257 257
258 return weight(kmWeight, tw, ow); 258 return Linear.weight(kmWeight, tw, ow);
259 } 259 }
260 } // class Row 260 } // class Row
261 261
262 protected List<Row> rows; 262 protected List<Row> rows;
263 263
336 Row r2 = rows.get(rowIndex); 336 Row r2 = rows.get(rowIndex);
337 337
338 return r1.interpolateWQ(r2, km, steps, this); 338 return r1.interpolateWQ(r2, km, steps, this);
339 } 339 }
340 340
341 public void interpolate( 341 public boolean interpolate(
342 double [] kms, 342 double km,
343 double [] ws, 343 double [] out,
344 double [] qs, 344 QPosition qPosition,
345 QPosition qPosition, 345 Function qFunction
346 LinearRemap remap
347 ) {
348 interpolate(kms, ws, qs, qPosition, remap, 0, kms.length);
349 }
350
351 public void interpolate(
352 double [] kms,
353 double [] ws,
354 double [] qs,
355 QPosition qPosition,
356 LinearRemap remap,
357 int startIndex,
358 int length
359 ) { 346 ) {
360 int R1 = rows.size()-1; 347 int R1 = rows.size()-1;
361 348
362 Row kmKey = new Row(); 349 out[1] = qFunction.value(getQ(qPosition, km));
350
351 if (Double.isNaN(out[1])) {
352 return false;
353 }
363 354
364 QPosition nPosition = new QPosition(); 355 QPosition nPosition = new QPosition();
365 356 if (getQPosition(km, out[1], nPosition) == null) {
366 for (int i = startIndex, end = startIndex+length; i < end; ++i) { 357 return false;
367 kmKey.km = kms[i]; 358 }
368 359
369 qs[i] = remap.eval(kms[i], getQ(qPosition, kms[i])); 360 int rowIndex = Collections.binarySearch(rows, new Row(km));
370 361
371 if (getQPosition(kms[i], qs[i], nPosition) == null) { 362 if (rowIndex >= 0) {
372 ws[i] = Double.NaN; 363 // direct row match
373 continue; 364 out[0] = rows.get(rowIndex).getW(nPosition);
374 } 365 return !Double.isNaN(out[0]);
375 366 }
376 int rowIndex = Collections.binarySearch(rows, kmKey); 367
377 368 rowIndex = -rowIndex -1;
378 if (rowIndex >= 0) { 369
379 // direct row match 370 if (rowIndex < 1 || rowIndex >= R1) {
380 ws[i] = rows.get(rowIndex).getW(nPosition); 371 // do not extrapolate
381 continue; 372 return false;
382 } 373 }
383 374
384 rowIndex = -rowIndex -1; 375 Row r1 = rows.get(rowIndex-1);
385 376 Row r2 = rows.get(rowIndex);
386 if (rowIndex < 1 || rowIndex >= R1) { 377 out[0] = r1.getW(r2, km, nPosition);
387 // do not extrapolate 378
388 ws[i] = Double.NaN; 379 return !Double.isNaN(out[0]);
389 continue;
390 }
391 Row r1 = rows.get(rowIndex-1);
392 Row r2 = rows.get(rowIndex);
393 ws[i] = r1.getW(r2, kms[i], nPosition);
394 }
395 } 380 }
396 381
397 public QPosition interpolate( 382 public QPosition interpolate(
398 double q, 383 double q,
399 double referenceKm, 384 double referenceKm,
478 double qMin, qMax; 463 double qMin, qMax;
479 if (qLast < qCurrent) { qMin = qLast; qMax = qCurrent; } 464 if (qLast < qCurrent) { qMin = qLast; qMax = qCurrent; }
480 else { qMin = qCurrent; qMax = qLast; } 465 else { qMin = qCurrent; qMax = qLast; }
481 466
482 if (q > qMin && q < qMax) { 467 if (q > qMin && q < qMax) {
483 double weight = factor(q, qLast, qCurrent); 468 double weight = Linear.factor(q, qLast, qCurrent);
484 return qPosition.set(i, weight); 469 return qPosition.set(i, weight);
485 } 470 }
486 qLast = qCurrent; 471 qLast = qCurrent;
487 } 472 }
488 473
500 if (weight == 1d) { 485 if (weight == 1d) {
501 return columns[index].getQRangeTree().findQ(km); 486 return columns[index].getQRangeTree().findQ(km);
502 } 487 }
503 double q1 = columns[index-1].getQRangeTree().findQ(km); 488 double q1 = columns[index-1].getQRangeTree().findQ(km);
504 double q2 = columns[index ].getQRangeTree().findQ(km); 489 double q2 = columns[index ].getQRangeTree().findQ(km);
505 return weight(weight, q1, q2); 490 return Linear.weight(weight, q1, q2);
506 } 491 }
507 492
508 public static final double linear(
509 double x,
510 double x1, double x2,
511 double y1, double y2
512 ) {
513 // y1 = m*x1 + b
514 // y2 = m*x2 + b
515 // y2 - y1 = m*(x2 - x1)
516 // m = (y2 - y1)/(x2 - x1) # x2 != x1
517 // b = y1 - m*x1
518
519 if (x1 == x2) {
520 return 0.5*(y1 + y2);
521 }
522 double m = (y2 - y1)/(x2 - x1);
523 double b = y1 - m*x1;
524 return x*m + b;
525 }
526
527 public static final double factor(double x, double p1, double p2) {
528 // 0 = m*p1 + b <=> b = -m*p1
529 // 1 = m*p2 + b
530 // 1 = m*(p2 - p1)
531 // m = 1/(p2 - p1) # p1 != p2
532 // f(x) = x/(p2-p1) - p1/(p2-p1) <=> (x-p1)/(p2-p1)
533
534 return p1 == p2 ? 0.0 : (x-p1)/(p2-p1);
535 }
536
537 public static final double weight(double factor, double a, double b) {
538 //return (1.0-factor)*a + factor*b;
539 return a + factor*(b-a);
540 }
541 } 493 }
542 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 : 494 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
http://dive4elements.wald.intevation.org