Mercurial > dive4elements > river

comparison flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java @ 686:3dc61e00385e facets-slt

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merged with trunk and introduced hashing of computed values. flys-artifacts/branches/facets-slt@2126 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Ingo Weinzierl <ingo.weinzierl@intevation.de>
date Wed, 15 Jun 2011 15:28:54 +0000
parents 913b52064449
children eab5e5089d77
comparison
equal deleted inserted replaced
667:434146596838 686:3dc61e00385e
15 import org.apache.commons.math.analysis.interpolation.SplineInterpolator; 15 import org.apache.commons.math.analysis.interpolation.SplineInterpolator;
16 16
17 import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction; 17 import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
18 18
19 import org.apache.commons.math.ArgumentOutsideDomainException; 19 import org.apache.commons.math.ArgumentOutsideDomainException;
20
21 import org.apache.commons.math.exception.MathIllegalArgumentException;
20 22
21 public class WstValueTable 23 public class WstValueTable
22 implements Serializable 24 implements Serializable
23 { 25 {
24 private static Logger log = Logger.getLogger(WstValueTable.class); 26 private static Logger log = Logger.getLogger(WstValueTable.class);
105 public void interpolateW( 107 public void interpolateW(
106 Row other, 108 Row other,
107 double km, 109 double km,
108 double [] iqs, 110 double [] iqs,
109 double [] ows, 111 double [] ows,
110 WstValueTable table 112 WstValueTable table,
113 Calculation errors
111 ) { 114 ) {
112 double kmWeight = Linear.factor(km, this.km, other.km); 115 double kmWeight = Linear.factor(km, this.km, other.km);
113 116
114 QPosition qPosition = new QPosition(); 117 QPosition qPosition = new QPosition();
115 118
116 for (int i = 0; i < iqs.length; ++i) { 119 for (int i = 0; i < iqs.length; ++i) {
117 if (table.getQPosition(km, iqs[i], qPosition) != null) { 120 if (table.getQPosition(km, iqs[i], qPosition) != null) {
118 double wt = getW(qPosition); 121 double wt = getW(qPosition);
119 double wo = other.getW(qPosition); 122 double wo = other.getW(qPosition);
120 ows[i] = Linear.weight(kmWeight, wt, wo); 123 if (Double.isNaN(wt) || Double.isNaN(wo)) {
124 if (errors != null) {
125 // TODO: I18N
126 errors.addProblem(
127 km, "cannot find w for q = " + iqs[i]);
128 }
129 ows[i] = Double.NaN;
130 }
131 else {
132 ows[i] = Linear.weight(kmWeight, wt, wo);
133 }
121 } 134 }
122 else { 135 else {
136 if (errors != null) {
137 // TODO: I18N
138 errors.addProblem(km, "cannot find q = " + iqs[i]);
139 }
123 ows[i] = Double.NaN; 140 ows[i] = Double.NaN;
124 } 141 }
125 } 142 }
126 } 143 }
127 144
128 public double [][] interpolateWQ( 145 public double [][] interpolateWQ(
129 Row other, 146 Row other,
130 double km, 147 double km,
131 int steps, 148 int steps,
132 WstValueTable table 149 WstValueTable table,
150 Calculation errors
133 ) { 151 ) {
134 int W = Math.min(ws.length, other.ws.length); 152 int W = Math.min(ws.length, other.ws.length);
135 153
136 if (W < 1) { 154 if (W < 1) {
155 if (errors != null) {
156 // TODO: I18N
157 errors.addProblem("no ws found");
158 }
137 return new double[2][0]; 159 return new double[2][0];
138 } 160 }
139 161
140 double factor = Linear.factor(km, this.km, other.km); 162 double factor = Linear.factor(km, this.km, other.km);
141 163
150 double wqs = Linear.weight( 172 double wqs = Linear.weight(
151 factor, 173 factor,
152 table.getQIndex(i, km), 174 table.getQIndex(i, km),
153 table.getQIndex(i, other.km)); 175 table.getQIndex(i, other.km));
154 176
177 if (Double.isNaN(wws) || Double.isNaN(wqs)) {
178 if (errors != null) {
179 // TODO: I18N
180 errors.addProblem(km, "cannot find w or q");
181 }
182 }
183 else {
184 if (wqs < minQ) minQ = wqs;
185 if (wqs > maxQ) maxQ = wqs;
186 }
187
155 splineW[i] = wws; 188 splineW[i] = wws;
156 splineQ[i] = wqs; 189 splineQ[i] = wqs;
157 if (wqs < minQ) minQ = wqs;
158 if (wqs > maxQ) maxQ = wqs;
159 } 190 }
160 191
161 double stepWidth = (maxQ - minQ)/steps; 192 double stepWidth = (maxQ - minQ)/steps;
162 193
163 SplineInterpolator interpolator = new SplineInterpolator(); 194 SplineInterpolator interpolator = new SplineInterpolator();
164 PolynomialSplineFunction spline = 195 PolynomialSplineFunction spline;
165 interpolator.interpolate(splineQ, splineW); 196
197 try {
198 spline = interpolator.interpolate(splineQ, splineW);
199 }
200 catch (MathIllegalArgumentException miae) {
201 if (errors != null) {
202 // TODO: I18N
203 errors.addProblem(km, "spline creation failed");
204 }
205 log.error("spline creation failed");
206 return new double[2][0];
207 }
166 208
167 double [] outWs = new double[steps]; 209 double [] outWs = new double[steps];
168 double [] outQs = new double[steps]; 210 double [] outQs = new double[steps];
211
212 Arrays.fill(outWs, Double.NaN);
213 Arrays.fill(outQs, Double.NaN);
169 214
170 try { 215 try {
171 double q = minQ; 216 double q = minQ;
172 for (int i = 0; i < outWs.length; ++i, q += stepWidth) { 217 for (int i = 0; i < outWs.length; ++i, q += stepWidth) {
173 outWs[i] = spline.value(outQs[i] = q); 218 outWs[i] = spline.value(outQs[i] = q);
174 } 219 }
175 } 220 }
176 catch (ArgumentOutsideDomainException aode) { 221 catch (ArgumentOutsideDomainException aode) {
177 log.error("Spline interpolation failed.", aode); 222 if (errors != null) {
223 // TODO: I18N
224 errors.addProblem(km, "spline interpolation failed");
225 }
226 log.error("spline interpolation failed", aode);
178 } 227 }
179 228
180 return new double [][] { outWs, outQs }; 229 return new double [][] { outWs, outQs };
181 } 230 }
182 231
183 public double [][] interpolateWQ(int steps, WstValueTable table) { 232 public double [][] interpolateWQ(
184 233 int steps,
234 WstValueTable table,
235 Calculation errors
236 ) {
185 int W = ws.length; 237 int W = ws.length;
186 238
187 if (W < 1) { 239 if (W < 1) {
240 if (errors != null) {
241 // TODO: I18N
242 errors.addProblem(km, "no ws found");
243 }
188 return new double[2][0]; 244 return new double[2][0];
189 } 245 }
190 246
191 double [] splineQ = new double[W]; 247 double [] splineQ = new double[W];
192 248
193 double minQ = Double.MAX_VALUE; 249 double minQ = Double.MAX_VALUE;
194 double maxQ = -Double.MAX_VALUE; 250 double maxQ = -Double.MAX_VALUE;
195 251
196 QPosition qPosition = new QPosition();
197
198 for (int i = 0; i < W; ++i) { 252 for (int i = 0; i < W; ++i) {
199 splineQ[i] = table.getQIndex(i, km); 253 double sq = table.getQIndex(i, km);
200 if (splineQ[i] < minQ) minQ = splineQ[i]; 254 if (Double.isNaN(sq)) {
201 if (splineQ[i] > maxQ) maxQ = splineQ[i]; 255 if (errors != null) {
256 // TODO: I18N
257 errors.addProblem(
258 km, "no q found in " + (i+1) + " column");
259 }
260 }
261 else {
262 if (sq < minQ) minQ = sq;
263 if (sq > maxQ) maxQ = sq;
264 }
265 splineQ[i] = sq;
202 } 266 }
203 267
204 double stepWidth = (maxQ - minQ)/steps; 268 double stepWidth = (maxQ - minQ)/steps;
205 269
206 SplineInterpolator interpolator = new SplineInterpolator(); 270 SplineInterpolator interpolator = new SplineInterpolator();
207 271
208 PolynomialSplineFunction spline = 272 PolynomialSplineFunction spline;
209 interpolator.interpolate(splineQ, ws); 273
274 try {
275 spline = interpolator.interpolate(splineQ, ws);
276 }
277 catch (MathIllegalArgumentException miae) {
278 if (errors != null) {
279 // TODO: I18N
280 errors.addProblem(km, "spline creation failed");
281 }
282 log.error("spline creation failed");
283 return new double[2][0];
284 }
210 285
211 double [] outWs = new double[steps]; 286 double [] outWs = new double[steps];
212 double [] outQs = new double[steps]; 287 double [] outQs = new double[steps];
288
289 Arrays.fill(outWs, Double.NaN);
290 Arrays.fill(outQs, Double.NaN);
213 291
214 try { 292 try {
215 double q = minQ; 293 double q = minQ;
216 for (int i = 0; i < outWs.length; ++i, q += stepWidth) { 294 for (int i = 0; i < outWs.length; ++i, q += stepWidth) {
217 outWs[i] = spline.value(outQs[i] = q); 295 outWs[i] = spline.value(outQs[i] = q);
218 } 296 }
219 } 297 }
220 catch (ArgumentOutsideDomainException aode) { 298 catch (ArgumentOutsideDomainException aode) {
221 log.error("Spline interpolation failed.", aode); 299 if (errors != null) {
300 // TODO: I18N
301 errors.addProblem(km, "spline interpolation failed");
302 }
303 log.error("spline interpolation failed.", aode);
222 } 304 }
223 305
224 return new double [][] { outWs, outQs }; 306 return new double [][] { outWs, outQs };
225 } 307 }
226 308
280 public void sortRows() { 362 public void sortRows() {
281 Collections.sort(rows); 363 Collections.sort(rows);
282 } 364 }
283 365
284 public double [] interpolateW(double km, double [] qs, double [] ws) { 366 public double [] interpolateW(double km, double [] qs, double [] ws) {
285 367 return interpolateW(km, qs, ws, null);
368 }
369
370 public double [] interpolateW(
371 double km,
372 double [] qs,
373 double [] ws,
374 Calculation errors
375 ) {
286 int rowIndex = Collections.binarySearch(rows, new Row(km)); 376 int rowIndex = Collections.binarySearch(rows, new Row(km));
287 377
288 QPosition qPosition = new QPosition(); 378 QPosition qPosition = new QPosition();
289 379
290 if (rowIndex >= 0) { // direct row match 380 if (rowIndex >= 0) { // direct row match
291 Row row = rows.get(rowIndex); 381 Row row = rows.get(rowIndex);
292 for (int i = 0; i < qs.length; ++i) { 382 for (int i = 0; i < qs.length; ++i) {
293 ws[i] = getQPosition(km, qs[i], qPosition) != null 383 if (getQPosition(km, qs[i], qPosition) == null) {
294 ? row.getW(qPosition) 384 if (errors != null) {
295 : Double.NaN; 385 // TODO: I18N
386 errors.addProblem(km, "cannot find q = " + qs[i]);
387 }
388 ws[i] = Double.NaN;
389 }
390 else {
391 if (Double.isNaN(ws[i] = row.getW(qPosition))
392 && errors != null) {
393 // TODO: I18N
394 errors.addProblem(
395 km, "cannot find w for q = " + qs[i]);
396 }
397 }
296 } 398 }
297 } 399 }
298 else { // needs bilinear interpolation 400 else { // needs bilinear interpolation
299 rowIndex = -rowIndex -1; 401 rowIndex = -rowIndex -1;
300 402
301 if (rowIndex < 1 || rowIndex >= rows.size()) { 403 if (rowIndex < 1 || rowIndex >= rows.size()) {
302 // do not extrapolate 404 // do not extrapolate
303 Arrays.fill(ws, Double.NaN); 405 Arrays.fill(ws, Double.NaN);
406 if (errors != null) {
407 // TODO: I18N
408 errors.addProblem(km, "km not found");
409 }
304 } 410 }
305 else { 411 else {
306 Row r1 = rows.get(rowIndex-1); 412 Row r1 = rows.get(rowIndex-1);
307 Row r2 = rows.get(rowIndex); 413 Row r2 = rows.get(rowIndex);
308 r1.interpolateW(r2, km, qs, ws, this); 414 r1.interpolateW(r2, km, qs, ws, this, errors);
309 } 415 }
310 } 416 }
311 417
312 return ws; 418 return ws;
313 } 419 }
314 420
315 public double [][] interpolateWQ(double km) { 421 public double [][] interpolateWQ(double km) {
316 return interpolateWQ(km, DEFAULT_Q_STEPS); 422 return interpolateWQ(km, null);
317 } 423 }
318 424
319 public double [][] interpolateWQ(double km, int steps) { 425 public double [][] interpolateWQ(double km, Calculation errors) {
426 return interpolateWQ(km, DEFAULT_Q_STEPS, errors);
427 }
428
429 public double [][] interpolateWQ(double km, int steps, Calculation errors) {
320 430
321 int rowIndex = Collections.binarySearch(rows, new Row(km)); 431 int rowIndex = Collections.binarySearch(rows, new Row(km));
322 432
323 if (rowIndex >= 0) { // direct row match 433 if (rowIndex >= 0) { // direct row match
324 Row row = rows.get(rowIndex); 434 Row row = rows.get(rowIndex);
325 return row.interpolateWQ(steps, this); 435 return row.interpolateWQ(steps, this, errors);
326 } 436 }
327 437
328 rowIndex = -rowIndex -1; 438 rowIndex = -rowIndex -1;
329 439
330 if (rowIndex < 1 || rowIndex >= rows.size()) { 440 if (rowIndex < 1 || rowIndex >= rows.size()) {
331 // do not extrapolate 441 // do not extrapolate
442 if (errors != null) {
443 // TODO: I18N
444 errors.addProblem(km, "km not found");
445 }
332 return new double[2][0]; 446 return new double[2][0];
333 } 447 }
334 448
335 Row r1 = rows.get(rowIndex-1); 449 Row r1 = rows.get(rowIndex-1);
336 Row r2 = rows.get(rowIndex); 450 Row r2 = rows.get(rowIndex);
337 451
338 return r1.interpolateWQ(r2, km, steps, this); 452 return r1.interpolateWQ(r2, km, steps, this, errors);
339 } 453 }
340 454
341 public boolean interpolate( 455 public boolean interpolate(
342 double km, 456 double km,
343 double [] out, 457 double [] out,
378 492
379 return !Double.isNaN(out[0]); 493 return !Double.isNaN(out[0]);
380 } 494 }
381 495
382 public QPosition interpolate( 496 public QPosition interpolate(
383 double q, 497 double q,
384 double referenceKm, 498 double referenceKm,
385 double [] kms, 499 double [] kms,
386 double [] ws, 500 double [] ws,
387 double [] qs 501 double [] qs,
502 Calculation errors
388 ) { 503 ) {
389 return interpolate(q, referenceKm, kms, ws, qs, 0, kms.length); 504 return interpolate(
505 q, referenceKm, kms, ws, qs, 0, kms.length, errors);
390 } 506 }
391 507
392 public QPosition interpolate( 508 public QPosition interpolate(
393 double q, 509 double q,
394 double referenceKm, 510 double referenceKm,
395 double [] kms, 511 double [] kms,
396 double [] ws, 512 double [] ws,
397 double [] qs, 513 double [] qs,
398 int startIndex, 514 int startIndex,
399 int length 515 int length,
516 Calculation errors
400 ) { 517 ) {
401 QPosition qPosition = getQPosition(referenceKm, q); 518 QPosition qPosition = getQPosition(referenceKm, q);
402 519
403 if (qPosition == null) { 520 if (qPosition == null) {
404 // we cannot locate q at km 521 // we cannot locate q at km
522 Arrays.fill(ws, Double.NaN);
523 Arrays.fill(qs, Double.NaN);
524 if (errors != null) {
525 errors.addProblem(referenceKm, "cannot find q");
526 }
405 return null; 527 return null;
406 } 528 }
407 529
408 Row kmKey = new Row(); 530 Row kmKey = new Row();
409 531
410 int R1 = rows.size()-1; 532 int R1 = rows.size()-1;
411 533
412 for (int i = startIndex, end = startIndex+length; i < end; ++i) { 534 for (int i = startIndex, end = startIndex+length; i < end; ++i) {
535
536 if (Double.isNaN(qs[i] = getQ(qPosition, kms[i]))) {
537 if (errors != null) {
538 errors.addProblem(kms[i], "cannot find q");
539 }
540 ws[i] = Double.NaN;
541 continue;
542 }
543
413 kmKey.km = kms[i]; 544 kmKey.km = kms[i];
414
415 int rowIndex = Collections.binarySearch(rows, kmKey); 545 int rowIndex = Collections.binarySearch(rows, kmKey);
416
417 qs[i] = getQ(qPosition, kms[i]);
418 546
419 if (rowIndex >= 0) { 547 if (rowIndex >= 0) {
420 // direct row match 548 // direct row match
421 ws[i] = rows.get(rowIndex).getW(qPosition); 549 if (Double.isNaN(ws[i] = rows.get(rowIndex).getW(qPosition))
550 && errors != null) {
551 errors.addProblem(kms[i], "cannot find w");
552 }
422 continue; 553 continue;
423 } 554 }
424 555
425 rowIndex = -rowIndex -1; 556 rowIndex = -rowIndex -1;
426 557
427 if (rowIndex < 1 || rowIndex >= R1) { 558 if (rowIndex < 1 || rowIndex >= R1) {
428 // do not extrapolate 559 // do not extrapolate
560 if (errors != null) {
561 errors.addProblem(kms[i], "cannot find km");
562 }
429 ws[i] = Double.NaN; 563 ws[i] = Double.NaN;
430 continue; 564 continue;
431 } 565 }
432 Row r1 = rows.get(rowIndex-1); 566 Row r1 = rows.get(rowIndex-1);
433 Row r2 = rows.get(rowIndex); 567 Row r2 = rows.get(rowIndex);
434 568
435 ws[i] = r1.getW(r2, kms[i], qPosition); 569 if (Double.isNaN(ws[i] = r1.getW(r2, kms[i], qPosition))
570 && errors != null) {
571 errors.addProblem(kms[i], "cannot find w");
572 }
436 } 573 }
437 574
438 return qPosition; 575 return qPosition;
439 } 576 }
440 577
http://dive4elements.wald.intevation.org