comparison flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java @ 335:64cfbd631f29

Implemented the "Wasserstand/Wasserspiegellagen" calculation. flys-artifacts/trunk@1733 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Tue, 19 Apr 2011 17:43:39 +0000
parents e09634fbf6bc
children 4509ba8fae68
comparison
equal deleted inserted replaced
334:b7c8df643dc4 335:64cfbd631f29
6 import de.intevation.flys.model.Wst; 6 import de.intevation.flys.model.Wst;
7 import de.intevation.flys.model.WstColumn; 7 import de.intevation.flys.model.WstColumn;
8 8
9 import de.intevation.flys.backend.SessionHolder; 9 import de.intevation.flys.backend.SessionHolder;
10 10
11 import java.util.Arrays;
11 import java.util.ArrayList; 12 import java.util.ArrayList;
12 import java.util.Comparator; 13 import java.util.Comparator;
13 import java.util.List; 14 import java.util.List;
14 import java.util.Collections; 15 import java.util.Collections;
15 import java.util.Iterator; 16 import java.util.Iterator;
51 } 52 }
52 } 53 }
53 // class Column 54 // class Column
54 55
55 public static class Row 56 public static class Row
56 implements Serializable 57 implements Serializable, Comparable<Row>
57 { 58 {
58 double km; 59 double km;
59 double [] wq; 60 double [] ws;
61 double [] qs;
62 boolean qSorted;
60 63
61 public Row() { 64 public Row() {
62 } 65 }
63 66
64 public Row(double km, double [] wq) { 67 public Row(double km) {
65 this.km = km; 68 this.km = km;
66 this.wq = wq; 69 }
70
71 public Row(double km, double [] ws, double [] qs) {
72 this(km);
73 this.ws = ws;
74 this.qs = qs;
75 }
76
77 public static final double linear(
78 double x,
79 double x1, double x2,
80 double y1, double y2
81 ) {
82 // y1 = m*x1 + b
83 // y2 = m*x2 + b
84 // y2 - y1 = m*(x2 - x1)
85 // m = (y2 - y1)/(x2 - x1) # x2 != x1
86 // b = y1 - m*x1
87
88 if (x1 == x2) {
89 return 0.5*(y1 + y2);
90 }
91 double m = (y2 - y1)/(x2 - x1);
92 double b = y1 - m*x1;
93 return x*m + b;
94 }
95
96 public static final double factor(double x, double p1, double p2) {
97 // 0 = m*p1 + b <=> b = -m*p1
98 // 1 = m*p2 + b
99 // 1 = m*(p2 - p1)
100 // m = 1/(p2 - p1) # p1 != p2
101 // f(x) = x/(p2-p1) - p1/(p2-p1) <=> (x-p1)/(p2-p1)
102
103 return p1 == p2 ? 0.0 : (x-p1)/(p2-p1);
104 }
105
106 public static final double weight(double factor, double a, double b) {
107 return (1.0-factor)*a + factor*b;
108 }
109
110 public double getWForKM(Row other, int index, double km) {
111 double w1 = ws[index];
112 double w2 = other.ws[index];
113 double km1 = this.km;
114 double km2 = other.km;
115 return linear(km, km1, km2, w1, w2);
116 }
117
118 public void interpolateW(
119 Row other,
120 double km,
121 double [] input,
122 double [] output
123 ) {
124 double factor = factor(km, this.km, other.km);
125
126 for (int i = 0; i < input.length; ++i) {
127 double q = input[i];
128 int idx1 = getQIndex(q);
129 int idx2 = other.getQIndex(q);
130
131 double w1 = idx1 >= 0
132 ? ws[idx1]
133 : interpolateW(-idx1-1, q);
134
135 double w2 = idx2 >= 0
136 ? other.ws[idx2]
137 : other.interpolateW(-idx2-1, q);
138
139 output[i] = weight(factor, w1, w2);
140 }
141 }
142
143 public double interpolateW(int idx, double q) {
144 return idx < 1 || idx >= qs.length
145 ? Double.NaN // do not extrapolate
146 : linear(q, qs[idx-1], qs[idx], ws[idx-1], ws[idx]);
147 }
148
149 public double interpolateW(double q) {
150 if (Double.isNaN(q)) return Double.NaN;
151 int index = getQIndex(q);
152 return index >= 0 ? ws[index] : interpolateW(-index -1);
153 }
154
155 public int getQIndex(double q) {
156 return qSorted ? binaryQIndex(q) : linearQIndex(q);
157 }
158
159 public int binaryQIndex(double q) {
160 return Arrays.binarySearch(qs, q);
161 }
162
163 public int linearQIndex(double q) {
164 switch (qs.length) {
165 case 0: break;
166 case 1:
167 if (qs[0] == q) return 0;
168 if (qs[0] < q) return -(1+1);
169 return -(0+1);
170 default:
171 for (int i = 1; i < qs.length; ++i) {
172 double qa = qs[i-1];
173 double qb = qs[i];
174 if (qa == q) return i-1;
175 if (qb == q) return i;
176 if (qa > qb) { double t = qa; qa = qb; qb = t; }
177 if (q > qa && q < qb) return -(i+1);
178 }
179 return -qs.length;
180 }
181
182 return -1;
183 }
184
185 public int compareTo(Row other) {
186 double d = km - other.km;
187 if (d < 0.0) return -1;
188 if (d > 0.0) return +1;
189 return 0;
67 } 190 }
68 } 191 }
69 // class Row 192 // class Row
70 193
71 protected List<Row> rows; 194 protected List<Row> rows;
77 } 200 }
78 201
79 public WstValueTable(Column [] columns) { 202 public WstValueTable(Column [] columns) {
80 this(); 203 this();
81 this.columns = columns; 204 this.columns = columns;
205 }
206
207
208 public double [] interpolateW(double km, double [] qs) {
209
210 double [] ws = new double[qs.length];
211
212 int rowIndex = Collections.binarySearch(rows, new Row(km));
213
214 if (rowIndex >= 0) { // direct row match
215 Row row = rows.get(rowIndex);
216 for (int i = 0; i < qs.length; ++i) {
217 ws[i] = row.interpolateW(qs[i]);
218 }
219 }
220 else { // needs bilinear interpolation
221 rowIndex = -rowIndex -1;
222
223 if (rowIndex < 1 || rowIndex >= rows.size()) {
224 // do not extrapolate
225 Arrays.fill(ws, Double.NaN);
226 }
227 else {
228 rows.get(rowIndex-1).interpolateW(
229 rows.get(rowIndex),
230 km, qs, ws);
231 }
232 }
233
234 return ws;
82 } 235 }
83 236
84 public static WstValueTable getTable(River river) { 237 public static WstValueTable getTable(River river) {
85 return getTable(river, 0); 238 return getTable(river, 0);
86 } 239 }
131 WstValueTable valueTable = new WstValueTable(columns); 284 WstValueTable valueTable = new WstValueTable(columns);
132 285
133 int lastColumnNo = -1; 286 int lastColumnNo = -1;
134 Row row = null; 287 Row row = null;
135 288
289 Double lastQ = -Double.MAX_VALUE;
290 boolean qSorted = true;
291
136 for (Iterator<Object []> iter = sqlQuery.iterate(); iter.hasNext();) { 292 for (Iterator<Object []> iter = sqlQuery.iterate(); iter.hasNext();) {
137 Object [] result = iter.next(); 293 Object [] result = iter.next();
138 double km = (Double) result[0]; 294 double km = (Double) result[0];
139 Double w = (Double) result[1]; 295 Double w = (Double) result[1];
140 Double q = (Double) result[2]; 296 Double q = (Double) result[2];
141 int columnNo = (Integer)result[3]; 297 int columnNo = (Integer)result[3];
142 298
143 if (columnNo > lastColumnNo) { // new row 299 if (columnNo > lastColumnNo) { // new row
144 if (row != null) { 300 if (row != null) {
301 row.qSorted = qSorted;
145 valueTable.rows.add(row); 302 valueTable.rows.add(row);
146 } 303 }
147 row = new Row(km, new double[(columnNo+1) << 1]); 304 row = new Row(
148 } 305 km,
149 int idx = columnNo << 1; 306 new double[columnNo+1],
150 307 new double[columnNo+1]);
151 row.wq[idx ] = w != null ? w : Double.NaN; 308 lastQ = -Double.MAX_VALUE;
152 row.wq[idx+1] = q != null ? q : Double.NaN; 309 qSorted = true;
310 }
311
312 row.ws[columnNo] = w != null ? w : Double.NaN;
313 row.qs[columnNo] = q != null ? q : Double.NaN;
314
315 if (qSorted && (q == null || lastQ > q)) {
316 qSorted = false;
317 }
318 lastQ = q;
153 319
154 lastColumnNo = columnNo; 320 lastColumnNo = columnNo;
155 } 321 }
156 322
157 if (row != null) { 323 if (row != null) {

http://dive4elements.wald.intevation.org