Mercurial > dive4elements > river
comparison flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java @ 633:d08f77e7f7e8
WST value table: Qs are now stored in ranges for each column.
flys-artifacts/trunk@2006 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Wed, 25 May 2011 15:31:25 +0000 |
parents | 07640ab913fd |
children | 433f67a076aa |
comparison
equal
deleted
inserted
replaced
632:07640ab913fd | 633:d08f77e7f7e8 |
---|---|
23 { | 23 { |
24 private static Logger log = Logger.getLogger(WstValueTable.class); | 24 private static Logger log = Logger.getLogger(WstValueTable.class); |
25 | 25 |
26 public static final int DEFAULT_Q_STEPS = 500; | 26 public static final int DEFAULT_Q_STEPS = 500; |
27 | 27 |
28 public static class Column | 28 public static final class Column |
29 implements Serializable | 29 implements Serializable |
30 { | 30 { |
31 protected String name; | 31 protected String name; |
32 | 32 |
33 protected QRangeTree qRangeTree; | 33 protected QRangeTree qRangeTree; |
34 | 34 |
52 } | 52 } |
53 | 53 |
54 public void setQRangeTree(QRangeTree qRangeTree) { | 54 public void setQRangeTree(QRangeTree qRangeTree) { |
55 this.qRangeTree = qRangeTree; | 55 this.qRangeTree = qRangeTree; |
56 } | 56 } |
57 } | 57 } // class Column |
58 // class Column | 58 |
59 | 59 public static final class QPosition { |
60 public static class QPosition { | 60 |
61 | |
62 protected double q; | |
63 protected int index; | 61 protected int index; |
64 protected double weight; | 62 protected double weight; |
65 | 63 |
66 public QPosition() { | 64 public QPosition() { |
67 } | 65 } |
68 | 66 |
69 public QPosition(double q, int index, double weight) { | 67 public QPosition(int index, double weight) { |
70 this.index = index; | 68 this.index = index; |
71 this.weight = weight; | 69 this.weight = weight; |
72 } | 70 } |
73 | 71 |
72 public QPosition set(int index, double weight) { | |
73 this.index = index; | |
74 this.weight = weight; | |
75 return this; | |
76 } | |
77 | |
74 } // class Position | 78 } // class Position |
75 | 79 |
76 public static class Row | 80 public static final class Row |
77 implements Serializable, Comparable<Row> | 81 implements Serializable, Comparable<Row> |
78 { | 82 { |
79 double km; | 83 double km; |
80 double [] ws; | 84 double [] ws; |
81 double [] qs; | |
82 boolean qSorted; | |
83 | 85 |
84 public Row() { | 86 public Row() { |
85 } | 87 } |
86 | 88 |
87 public Row(double km) { | 89 public Row(double km) { |
88 this.km = km; | 90 this.km = km; |
89 } | 91 } |
90 | 92 |
91 public Row(double km, double [] ws, double [] qs) { | 93 public Row(double km, double [] ws) { |
92 this(km); | 94 this(km); |
93 this.ws = ws; | 95 this.ws = ws; |
94 this.qs = qs; | |
95 } | |
96 | |
97 public static final double linear( | |
98 double x, | |
99 double x1, double x2, | |
100 double y1, double y2 | |
101 ) { | |
102 // y1 = m*x1 + b | |
103 // y2 = m*x2 + b | |
104 // y2 - y1 = m*(x2 - x1) | |
105 // m = (y2 - y1)/(x2 - x1) # x2 != x1 | |
106 // b = y1 - m*x1 | |
107 | |
108 if (x1 == x2) { | |
109 return 0.5*(y1 + y2); | |
110 } | |
111 double m = (y2 - y1)/(x2 - x1); | |
112 double b = y1 - m*x1; | |
113 return x*m + b; | |
114 } | |
115 | |
116 public static final double factor(double x, double p1, double p2) { | |
117 // 0 = m*p1 + b <=> b = -m*p1 | |
118 // 1 = m*p2 + b | |
119 // 1 = m*(p2 - p1) | |
120 // m = 1/(p2 - p1) # p1 != p2 | |
121 // f(x) = x/(p2-p1) - p1/(p2-p1) <=> (x-p1)/(p2-p1) | |
122 | |
123 return p1 == p2 ? 0.0 : (x-p1)/(p2-p1); | |
124 } | |
125 | |
126 public static final double weight(double factor, double a, double b) { | |
127 return (1.0-factor)*a + factor*b; | |
128 } | |
129 | |
130 public void interpolateW( | |
131 Row other, | |
132 double km, | |
133 double [] input, | |
134 double [] output | |
135 ) { | |
136 double factor = factor(km, this.km, other.km); | |
137 | |
138 for (int i = 0; i < input.length; ++i) { | |
139 double q = input[i]; | |
140 int idx1 = getQIndex(q); | |
141 int idx2 = other.getQIndex(q); | |
142 | |
143 double w1 = idx1 >= 0 | |
144 ? ws[idx1] | |
145 : interpolateW(-idx1-1, q); | |
146 | |
147 double w2 = idx2 >= 0 | |
148 ? other.ws[idx2] | |
149 : other.interpolateW(-idx2-1, q); | |
150 | |
151 output[i] = weight(factor, w1, w2); | |
152 } | |
153 } | |
154 | |
155 public double interpolateW(int idx, double q) { | |
156 return idx < 1 || idx >= qs.length | |
157 ? Double.NaN // do not extrapolate | |
158 : linear(q, qs[idx-1], qs[idx], ws[idx-1], ws[idx]); | |
159 } | |
160 | |
161 public double interpolateW(double q) { | |
162 if (Double.isNaN(q)) return Double.NaN; | |
163 int index = getQIndex(q); | |
164 return index >= 0 ? ws[index] : interpolateW(-index -1, q); | |
165 } | |
166 | |
167 public int getQIndex(double q) { | |
168 return qSorted ? binaryQIndex(q) : linearQIndex(q); | |
169 } | |
170 | |
171 public int binaryQIndex(double q) { | |
172 return Arrays.binarySearch(qs, q); | |
173 } | |
174 | |
175 public int linearQIndex(double q) { | |
176 switch (qs.length) { | |
177 case 0: break; | |
178 case 1: | |
179 if (qs[0] == q) return 0; | |
180 if (qs[0] < q) return -(1+1); | |
181 return -(0+1); | |
182 default: | |
183 for (int i = 1; i < qs.length; ++i) { | |
184 double qa = qs[i-1]; | |
185 double qb = qs[i]; | |
186 if (qa == q) return i-1; | |
187 if (qb == q) return i; | |
188 if (qa > qb) { double t = qa; qa = qb; qb = t; } | |
189 if (q > qa && q < qb) return -(i+1); | |
190 } | |
191 return -qs.length; | |
192 } | |
193 | |
194 return -1; | |
195 } | 96 } |
196 | 97 |
197 public int compareTo(Row other) { | 98 public int compareTo(Row other) { |
198 double d = km - other.km; | 99 double d = km - other.km; |
199 if (d < -0.0001) return -1; | 100 if (d < -0.0001) return -1; |
200 if (d > 0.0001) return +1; | 101 if (d > 0.0001) return +1; |
201 return 0; | 102 return 0; |
202 } | 103 } |
203 | 104 |
204 public double [][] interpolateWQ(Row other, double km, int steps) { | 105 public void interpolateW( |
205 | 106 Row other, |
107 double km, | |
108 double [] iqs, | |
109 double [] ows, | |
110 WstValueTable table | |
111 ) { | |
112 double kmWeight = factor(km, this.km, other.km); | |
113 | |
114 QPosition qPosition = new QPosition(); | |
115 | |
116 for (int i = 0; i < iqs.length; ++i) { | |
117 if (table.getQPosition(km, iqs[i], qPosition) != null) { | |
118 double wt = getW(qPosition); | |
119 double wo = other.getW(qPosition); | |
120 ows[i] = weight(kmWeight, wt, wo); | |
121 } | |
122 else { | |
123 ows[i] = Double.NaN; | |
124 } | |
125 } | |
126 } | |
127 | |
128 public double [][] interpolateWQ( | |
129 Row other, | |
130 double km, | |
131 int steps, | |
132 WstValueTable table | |
133 ) { | |
206 int W = Math.min(ws.length, other.ws.length); | 134 int W = Math.min(ws.length, other.ws.length); |
207 | 135 |
208 if (W < 1) { | 136 if (W < 1) { |
209 return new double[2][0]; | 137 return new double[2][0]; |
210 } | 138 } |
217 double minQ = Double.MAX_VALUE; | 145 double minQ = Double.MAX_VALUE; |
218 double maxQ = -Double.MAX_VALUE; | 146 double maxQ = -Double.MAX_VALUE; |
219 | 147 |
220 for (int i = 0; i < W; ++i) { | 148 for (int i = 0; i < W; ++i) { |
221 double wws = weight(factor, ws[i], other.ws[i]); | 149 double wws = weight(factor, ws[i], other.ws[i]); |
222 double wqs = weight(factor, qs[i], other.qs[i]); | 150 double wqs = weight( |
151 factor, | |
152 table.getQIndex(i, km), | |
153 table.getQIndex(i, other.km)); | |
154 | |
223 splineW[i] = wws; | 155 splineW[i] = wws; |
224 splineQ[i] = wqs; | 156 splineQ[i] = wqs; |
225 if (wqs < minQ) minQ = wqs; | 157 if (wqs < minQ) minQ = wqs; |
226 if (wqs > maxQ) maxQ = wqs; | 158 if (wqs > maxQ) maxQ = wqs; |
227 } | 159 } |
246 } | 178 } |
247 | 179 |
248 return new double [][] { outWs, outQs }; | 180 return new double [][] { outWs, outQs }; |
249 } | 181 } |
250 | 182 |
251 public double [][] interpolateWQ(int steps) { | 183 public double [][] interpolateWQ(int steps, WstValueTable table) { |
252 | 184 |
253 int W = ws.length; | 185 int W = ws.length; |
254 | 186 |
255 if (W < 1) { | 187 if (W < 1) { |
256 return new double[2][0]; | 188 return new double[2][0]; |
257 } | 189 } |
258 | 190 |
259 double [] splineW = new double[W]; | |
260 double [] splineQ = new double[W]; | 191 double [] splineQ = new double[W]; |
261 | 192 |
262 double minQ = Double.MAX_VALUE; | 193 double minQ = Double.MAX_VALUE; |
263 double maxQ = -Double.MAX_VALUE; | 194 double maxQ = -Double.MAX_VALUE; |
264 | 195 |
196 QPosition qPosition = new QPosition(); | |
197 | |
265 for (int i = 0; i < W; ++i) { | 198 for (int i = 0; i < W; ++i) { |
266 splineW[i] = ws[i]; | 199 splineQ[i] = table.getQIndex(i, km); |
267 splineQ[i] = qs[i]; | 200 if (splineQ[i] < minQ) minQ = splineQ[i]; |
268 if (qs[i] < minQ) minQ = qs[i]; | 201 if (splineQ[i] > maxQ) maxQ = splineQ[i]; |
269 if (qs[i] > maxQ) maxQ = qs[i]; | |
270 } | 202 } |
271 | 203 |
272 double stepWidth = (maxQ - minQ)/steps; | 204 double stepWidth = (maxQ - minQ)/steps; |
273 | 205 |
274 SplineInterpolator interpolator = new SplineInterpolator(); | 206 SplineInterpolator interpolator = new SplineInterpolator(); |
275 | 207 |
276 PolynomialSplineFunction spline = | 208 PolynomialSplineFunction spline = |
277 interpolator.interpolate(splineQ, splineW); | 209 interpolator.interpolate(splineQ, ws); |
278 | 210 |
279 double [] outWs = new double[steps]; | 211 double [] outWs = new double[steps]; |
280 double [] outQs = new double[steps]; | 212 double [] outQs = new double[steps]; |
281 | 213 |
282 try { | 214 try { |
290 } | 222 } |
291 | 223 |
292 return new double [][] { outWs, outQs }; | 224 return new double [][] { outWs, outQs }; |
293 } | 225 } |
294 | 226 |
295 public QPosition getQPosition(double q) { | 227 |
296 int qIndex = getQIndex(q); | 228 public double getW(QPosition qPosition) { |
297 | 229 int index = qPosition.index; |
298 if (qIndex >= 0) { | 230 double weight = qPosition.weight; |
299 // direct column match | 231 |
300 return new QPosition(q, qIndex, 1.0); | 232 return weight == 1.0 |
301 } | 233 ? ws[index] |
302 | 234 : weight(weight, ws[index-1], ws[index]); |
303 qIndex = -qIndex -1; | 235 } |
304 | 236 |
305 if (qIndex < 0 || qIndex >= qs.length-1) { | 237 public double getW( |
306 // do not extrapolate | |
307 return null; | |
308 } | |
309 | |
310 return new QPosition( | |
311 q, qIndex, factor(q, qs[qIndex-1], qs[qIndex])); | |
312 } | |
313 | |
314 public QPosition getQPosition(Row other, double km, double q) { | |
315 | |
316 QPosition qpt = getQPosition(q); | |
317 QPosition qpo = other.getQPosition(q); | |
318 | |
319 if (qpt == null || qpo == null) { | |
320 return null; | |
321 } | |
322 | |
323 double kmWeight = factor(km, this.km, other.km); | |
324 | |
325 // XXX: Index interpolation is a bit sticky here! | |
326 | |
327 int index = (int)Math.round( | |
328 weight(kmWeight, qpt.index, qpo.index)); | |
329 | |
330 double weight = weight(kmWeight, qpt.weight, qpo.weight); | |
331 | |
332 return new QPosition(q, index, weight); | |
333 } | |
334 | |
335 public void storeWQ( | |
336 QPosition qPosition, | |
337 double [] ows, | |
338 double [] oqs, | |
339 int index | |
340 ) { | |
341 int qIdx = qPosition.index; | |
342 double qWeight = qPosition.weight; | |
343 | |
344 if (qWeight == 1.0) { | |
345 oqs[index] = qs[qIdx]; | |
346 ows[index] = ws[qIdx]; | |
347 } | |
348 else { | |
349 oqs[index] = weight(qWeight, qs[qIdx-1], qs[qIdx]); | |
350 ows[index] = weight(qWeight, ws[qIdx-1], ws[qIdx]); | |
351 } | |
352 } | |
353 | |
354 public void storeWQ( | |
355 Row other, | 238 Row other, |
356 double km, | 239 double km, |
357 QPosition qPosition, | 240 QPosition qPosition |
358 double [] ows, | |
359 double [] oqs, | |
360 int index | |
361 ) { | 241 ) { |
362 double kmWeight = factor(km, this.km, other.km); | 242 double kmWeight = factor(km, this.km, other.km); |
363 | 243 |
364 double tq, tw; | 244 int index = qPosition.index; |
365 double oq, ow; | 245 double weight = qPosition.weight; |
366 | 246 |
367 int qIdx = qPosition.index; | 247 double tw, ow; |
368 double qWeight = qPosition.weight; | 248 |
369 | 249 if (weight == 1.0) { |
370 if (qWeight == 1.0) { | 250 tw = ws[index]; |
371 tq = qs[qIdx]; | 251 ow = other.ws[index]; |
372 tw = ws[qIdx]; | |
373 oq = other.qs[qIdx]; | |
374 ow = other.ws[qIdx]; | |
375 } | 252 } |
376 else { | 253 else { |
377 tq = weight(qWeight, qs[qIdx-1], qs[qIdx]); | 254 tw = weight(weight, ws[index-1], ws[index]); |
378 tw = weight(qWeight, ws[qIdx-1], ws[qIdx]); | 255 ow = weight(weight, other.ws[index-1], other.ws[index]); |
379 oq = weight(qWeight, other.qs[qIdx-1], other.qs[qIdx]); | 256 } |
380 ow = weight(qWeight, other.ws[qIdx-1], other.ws[qIdx]); | 257 |
381 } | 258 return weight(kmWeight, tw, ow); |
382 | 259 } |
383 oqs[index] = weight(kmWeight, tq, oq); | 260 } // class Row |
384 ows[index] = weight(kmWeight, tw, ow); | |
385 } | |
386 | |
387 public void storeWQ( | |
388 QPosition qPosition, | |
389 LinearRemap remap, | |
390 double [] ows, | |
391 double [] oqs, | |
392 int index | |
393 ) { | |
394 int qIdx = qPosition.index; | |
395 double qWeight = qPosition.weight; | |
396 | |
397 oqs[index] = remap.eval(km, qWeight == 1.0 | |
398 ? qs[qIdx] | |
399 : weight(qWeight, qs[qIdx-1], qs[qIdx])); | |
400 | |
401 ows[index] = interpolateW(oqs[index]); | |
402 } | |
403 | |
404 public void storeWQ( | |
405 Row other, | |
406 double km, | |
407 QPosition qPosition, | |
408 LinearRemap remap, | |
409 double [] ows, | |
410 double [] oqs, | |
411 int index | |
412 ) { | |
413 double kmWeight = factor(km, this.km, other.km); | |
414 | |
415 double tq, tw; | |
416 double oq, ow; | |
417 | |
418 int qIdx = qPosition.index; | |
419 double qWeight = qPosition.weight; | |
420 | |
421 if (qWeight == 1.0) { | |
422 tq = remap.eval(this.km, qs[qIdx]); | |
423 oq = remap.eval(other.km, other.qs[qIdx]); | |
424 } | |
425 else { | |
426 tq = remap.eval( | |
427 this.km, | |
428 weight(qWeight, qs[qIdx-1], qs[qIdx])); | |
429 oq = remap.eval( | |
430 other.km, | |
431 weight(qWeight, other.qs[qIdx-1], other.qs[qIdx])); | |
432 } | |
433 | |
434 tw = interpolateW(tq); | |
435 ow = other.interpolateW(oq); | |
436 | |
437 oqs[index] = weight(kmWeight, tq, oq); | |
438 ows[index] = weight(kmWeight, tw, ow); | |
439 } | |
440 } | |
441 // class Row | |
442 | 261 |
443 protected List<Row> rows; | 262 protected List<Row> rows; |
444 | 263 |
445 protected Column [] columns; | 264 protected Column [] columns; |
446 | 265 |
451 public WstValueTable(Column [] columns) { | 270 public WstValueTable(Column [] columns) { |
452 this(); | 271 this(); |
453 this.columns = columns; | 272 this.columns = columns; |
454 } | 273 } |
455 | 274 |
275 public WstValueTable(Column [] columns, List<Row> rows) { | |
276 this.columns = columns; | |
277 this.rows = rows; | |
278 } | |
279 | |
456 public void sortRows() { | 280 public void sortRows() { |
457 Collections.sort(rows); | 281 Collections.sort(rows); |
458 } | 282 } |
459 | 283 |
460 | |
461 public double [] interpolateW(double km, double [] qs) { | |
462 return interpolateW(km, qs, new double[qs.length]); | |
463 } | |
464 | |
465 public double [] interpolateW(double km, double [] qs, double [] ws) { | 284 public double [] interpolateW(double km, double [] qs, double [] ws) { |
466 | 285 |
467 int rowIndex = Collections.binarySearch(rows, new Row(km)); | 286 int rowIndex = Collections.binarySearch(rows, new Row(km)); |
287 | |
288 QPosition qPosition = new QPosition(); | |
468 | 289 |
469 if (rowIndex >= 0) { // direct row match | 290 if (rowIndex >= 0) { // direct row match |
470 Row row = rows.get(rowIndex); | 291 Row row = rows.get(rowIndex); |
471 for (int i = 0; i < qs.length; ++i) { | 292 for (int i = 0; i < qs.length; ++i) { |
472 ws[i] = row.interpolateW(qs[i]); | 293 ws[i] = getQPosition(km, qs[i], qPosition) != null |
294 ? row.getW(qPosition) | |
295 : Double.NaN; | |
473 } | 296 } |
474 } | 297 } |
475 else { // needs bilinear interpolation | 298 else { // needs bilinear interpolation |
476 rowIndex = -rowIndex -1; | 299 rowIndex = -rowIndex -1; |
477 | 300 |
478 if (rowIndex < 1 || rowIndex >= rows.size()) { | 301 if (rowIndex < 1 || rowIndex >= rows.size()) { |
479 // do not extrapolate | 302 // do not extrapolate |
480 Arrays.fill(ws, Double.NaN); | 303 Arrays.fill(ws, Double.NaN); |
481 } | 304 } |
482 else { | 305 else { |
483 rows.get(rowIndex-1).interpolateW( | 306 Row r1 = rows.get(rowIndex-1); |
484 rows.get(rowIndex), | 307 Row r2 = rows.get(rowIndex); |
485 km, qs, ws); | 308 r1.interpolateW(r2, km, qs, ws, this); |
486 } | 309 } |
487 } | 310 } |
488 | 311 |
489 return ws; | 312 return ws; |
490 } | 313 } |
491 | |
492 | 314 |
493 public double [][] interpolateWQ(double km) { | 315 public double [][] interpolateWQ(double km) { |
494 return interpolateWQ(km, DEFAULT_Q_STEPS); | 316 return interpolateWQ(km, DEFAULT_Q_STEPS); |
495 } | 317 } |
496 | 318 |
498 | 320 |
499 int rowIndex = Collections.binarySearch(rows, new Row(km)); | 321 int rowIndex = Collections.binarySearch(rows, new Row(km)); |
500 | 322 |
501 if (rowIndex >= 0) { // direct row match | 323 if (rowIndex >= 0) { // direct row match |
502 Row row = rows.get(rowIndex); | 324 Row row = rows.get(rowIndex); |
503 return row.interpolateWQ(steps); | 325 return row.interpolateWQ(steps, this); |
504 } | 326 } |
505 | 327 |
506 rowIndex = -rowIndex -1; | 328 rowIndex = -rowIndex -1; |
507 | 329 |
508 if (rowIndex < 1 || rowIndex >= rows.size()) { | 330 if (rowIndex < 1 || rowIndex >= rows.size()) { |
511 } | 333 } |
512 | 334 |
513 Row r1 = rows.get(rowIndex-1); | 335 Row r1 = rows.get(rowIndex-1); |
514 Row r2 = rows.get(rowIndex); | 336 Row r2 = rows.get(rowIndex); |
515 | 337 |
516 return r1.interpolateWQ(r2, km, steps); | 338 return r1.interpolateWQ(r2, km, steps, this); |
517 } | 339 } |
518 | 340 |
519 public void interpolate( | 341 public void interpolate( |
520 double [] kms, | 342 double [] kms, |
521 double [] ws, | 343 double [] ws, |
524 LinearRemap remap | 346 LinearRemap remap |
525 ) { | 347 ) { |
526 int R1 = rows.size()-1; | 348 int R1 = rows.size()-1; |
527 | 349 |
528 Row kmKey = new Row(); | 350 Row kmKey = new Row(); |
351 | |
352 QPosition nPosition = new QPosition(); | |
353 | |
529 for (int i = 0; i < kms.length; ++i) { | 354 for (int i = 0; i < kms.length; ++i) { |
530 kmKey.km = kms[i]; | 355 kmKey.km = kms[i]; |
356 | |
357 qs[i] = remap.eval(kms[i], getQ(qPosition, kms[i])); | |
358 | |
359 if (getQPosition(kms[i], qs[i], nPosition) == null) { | |
360 ws[i] = Double.NaN; | |
361 continue; | |
362 } | |
363 | |
531 int rowIndex = Collections.binarySearch(rows, kmKey); | 364 int rowIndex = Collections.binarySearch(rows, kmKey); |
532 | 365 |
533 if (rowIndex >= 0) { | 366 if (rowIndex >= 0) { |
534 // direct row match | 367 // direct row match |
535 rows.get(rowIndex).storeWQ(qPosition, remap, ws, qs, i); | 368 ws[i] = rows.get(rowIndex).getW(nPosition); |
536 continue; | 369 continue; |
537 } | 370 } |
538 | 371 |
539 rowIndex = -rowIndex -1; | 372 rowIndex = -rowIndex -1; |
540 | 373 |
541 if (rowIndex < 1 || rowIndex >= R1) { | 374 if (rowIndex < 1 || rowIndex >= R1) { |
542 // do not extrapolate | 375 // do not extrapolate |
543 ws[i] = qs[i] = Double.NaN; | 376 ws[i] = Double.NaN; |
544 continue; | 377 continue; |
545 } | 378 } |
546 Row r1 = rows.get(rowIndex-1); | 379 Row r1 = rows.get(rowIndex-1); |
547 Row r2 = rows.get(rowIndex); | 380 Row r2 = rows.get(rowIndex); |
548 r1.storeWQ(r2, kms[i], qPosition, remap, ws, qs, i); | 381 ws[i] = r1.getW(r2, kms[i], nPosition); |
549 } | 382 } |
550 } | 383 } |
551 | 384 |
552 public QPosition interpolate( | 385 public QPosition interpolate( |
553 double q, | 386 double q, |
560 | 393 |
561 int rowIndex = Collections.binarySearch(rows, kmKey); | 394 int rowIndex = Collections.binarySearch(rows, kmKey); |
562 | 395 |
563 int R1 = rows.size()-1; | 396 int R1 = rows.size()-1; |
564 | 397 |
565 QPosition qPosition = null; | 398 QPosition qPosition = getQPosition(kms[referenceIndex], q); |
566 | 399 |
567 if (rowIndex >= 0) { // direct row match | 400 if (qPosition == null) { |
568 qPosition = rows.get(rowIndex).getQPosition(q); | 401 // we cannot locate q at km |
569 } | 402 return null; |
570 else { | 403 } |
404 | |
405 for (int i = 0; i < kms.length; ++i) { | |
406 kmKey.km = kms[i]; | |
407 | |
408 rowIndex = Collections.binarySearch(rows, kmKey); | |
409 | |
410 qs[i] = getQ(qPosition, kms[i]); | |
411 | |
412 if (rowIndex >= 0) { | |
413 // direct row match | |
414 ws[i] = rows.get(rowIndex).getW(qPosition); | |
415 continue; | |
416 } | |
417 | |
571 rowIndex = -rowIndex -1; | 418 rowIndex = -rowIndex -1; |
572 | 419 |
573 if (rowIndex < 1 || rowIndex >= R1) { | 420 if (rowIndex < 1 || rowIndex >= R1) { |
574 // do not extrapolate | 421 // do not extrapolate |
575 return null; | 422 ws[i] = Double.NaN; |
576 } | 423 continue; |
577 | 424 } |
578 // interpolation case | |
579 Row r1 = rows.get(rowIndex-1); | 425 Row r1 = rows.get(rowIndex-1); |
580 Row r2 = rows.get(rowIndex); | 426 Row r2 = rows.get(rowIndex); |
581 | 427 |
582 qPosition = r1.getQPosition(r2, kms[referenceIndex], q); | 428 ws[i] = r1.getW(r2, kms[i], qPosition); |
583 } | 429 } |
584 | 430 |
585 if (qPosition == null) { | 431 return qPosition; |
586 // we cannot locate q in reference row | 432 } |
433 | |
434 public QPosition getQPosition(double km, double q) { | |
435 return getQPosition(km, q, new QPosition()); | |
436 } | |
437 | |
438 public QPosition getQPosition(double km, double q, QPosition qPosition) { | |
439 | |
440 if (columns.length == 0) { | |
587 return null; | 441 return null; |
588 } | 442 } |
589 | 443 |
590 for (int i = 0; i < kms.length; ++i) { | 444 double qLast = columns[0].getQRangeTree().findQ(km); |
591 kmKey.km = kms[i]; | 445 |
592 | 446 if (Math.abs(qLast - q) < 0.00001) { |
593 rowIndex = Collections.binarySearch(rows, kmKey); | 447 return qPosition.set(0, 1d); |
594 | 448 } |
595 if (rowIndex >= 0) { | 449 |
596 // direct row match | 450 for (int i = 1; i < columns.length; ++i) { |
597 rows.get(rowIndex).storeWQ(qPosition, ws, qs, i); | 451 double qCurrent = columns[i].getQRangeTree().findQ(km); |
598 continue; | 452 if (Math.abs(qCurrent - q) < 0.00001) { |
599 } | 453 return qPosition.set(i, 1d); |
600 | 454 } |
601 rowIndex = -rowIndex -1; | 455 |
602 | 456 double qMin, qMax; |
603 if (rowIndex < 1 || rowIndex >= R1) { | 457 if (qLast < qCurrent) { qMin = qLast; qMax = qCurrent; } |
604 // do not extrapolate | 458 else { qMin = qCurrent; qMax = qLast; } |
605 ws[i] = qs[i] = Double.NaN; | 459 |
606 continue; | 460 if (q > qMin && q < qMax) { |
607 } | 461 double weight = factor(q, qLast, qCurrent); |
608 Row r1 = rows.get(rowIndex-1); | 462 return qPosition.set(i, weight); |
609 Row r2 = rows.get(rowIndex); | 463 } |
610 r1.storeWQ(r2, kms[i], qPosition, ws, qs, i); | 464 qLast = qCurrent; |
611 } | 465 } |
612 | 466 |
613 return qPosition; | 467 return null; |
468 } | |
469 | |
470 public double getQIndex(int index, double km) { | |
471 return columns[index].getQRangeTree().findQ(km); | |
472 } | |
473 | |
474 public double getQ(QPosition qPosition, double km) { | |
475 int index = qPosition.index; | |
476 double weight = qPosition.weight; | |
477 | |
478 if (weight == 1d) { | |
479 return columns[index].getQRangeTree().findQ(km); | |
480 } | |
481 double q1 = columns[index-1].getQRangeTree().findQ(km); | |
482 double q2 = columns[index ].getQRangeTree().findQ(km); | |
483 return weight(weight, q1, q2); | |
484 } | |
485 | |
486 public static final double linear( | |
487 double x, | |
488 double x1, double x2, | |
489 double y1, double y2 | |
490 ) { | |
491 // y1 = m*x1 + b | |
492 // y2 = m*x2 + b | |
493 // y2 - y1 = m*(x2 - x1) | |
494 // m = (y2 - y1)/(x2 - x1) # x2 != x1 | |
495 // b = y1 - m*x1 | |
496 | |
497 if (x1 == x2) { | |
498 return 0.5*(y1 + y2); | |
499 } | |
500 double m = (y2 - y1)/(x2 - x1); | |
501 double b = y1 - m*x1; | |
502 return x*m + b; | |
503 } | |
504 | |
505 public static final double factor(double x, double p1, double p2) { | |
506 // 0 = m*p1 + b <=> b = -m*p1 | |
507 // 1 = m*p2 + b | |
508 // 1 = m*(p2 - p1) | |
509 // m = 1/(p2 - p1) # p1 != p2 | |
510 // f(x) = x/(p2-p1) - p1/(p2-p1) <=> (x-p1)/(p2-p1) | |
511 | |
512 return p1 == p2 ? 0.0 : (x-p1)/(p2-p1); | |
513 } | |
514 | |
515 public static final double weight(double factor, double a, double b) { | |
516 //return (1.0-factor)*a + factor*b; | |
517 return a + factor*(b-a); | |
614 } | 518 } |
615 } | 519 } |
616 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 : | 520 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 : |