Mercurial > dive4elements > river

comparison flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/Calculation4.java @ 3318:dbe2f85bf160

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merged flys-artifacts/2.8
author Thomas Arendsen Hein <thomas@intevation.de>
date Fri, 28 Sep 2012 12:14:35 +0200
parents a0d9a99a5d17
children cb11919cccf9
comparison
equal deleted inserted replaced
2987:98c7a46ec5ae 3318:dbe2f85bf160
1 package de.intevation.flys.artifacts.model;
2
3 import java.util.List;
4 import java.util.Arrays;
5 import java.util.Comparator;
6 import java.util.Collections;
7
8 import de.intevation.flys.utils.DoubleUtil;
9
10 import de.intevation.flys.model.River;
11 import de.intevation.flys.model.Gauge;
12 import de.intevation.flys.model.DischargeTable;
13
14 import de.intevation.flys.artifacts.model.WstValueTable.QPosition;
15
16 import de.intevation.flys.artifacts.math.Function;
17 import de.intevation.flys.artifacts.math.Linear;
18 import de.intevation.flys.artifacts.math.Identity;
19 import de.intevation.flys.artifacts.math.BackJumpCorrector;
20
21 import org.apache.log4j.Logger;
22
23 public class Calculation4
24 extends Calculation
25 {
26 private static Logger logger = Logger.getLogger(Calculation4.class);
27
28 public static final double MINIMAL_STEP_WIDTH = 1e-5;
29
30 public static final Comparator<Segment> REF_CMP =
31 new Comparator<Segment>() {
32 public int compare(Segment a, Segment b) {
33 double d = a.referencePoint - b.referencePoint;
34 if (d < 0d) return -1;
35 return d > 0d ? +1 : 0;
36 }
37 };
38
39 protected List<Segment> segments;
40
41 protected boolean isQ;
42
43 public Calculation4() {
44 }
45
46 public Calculation4(List<Segment> segments, River river, boolean isQ) {
47
48 this.segments = segments;
49 this.isQ = isQ;
50
51 int numResults = -1;
52
53 // assign reference points
54 for (Segment segment: segments) {
55 Gauge gauge = river.maxOverlap(segment.getFrom(), segment.getTo());
56
57 if (gauge == null) {
58 logger.warn("no gauge found. Defaults to mid point.");
59 segment.setReferencePoint(
60 0.5*(segment.getFrom()+segment.getTo()));
61 }
62 else {
63 double ref = gauge.getStation().doubleValue();
64 logger.debug(
65 "reference gauge: " + gauge.getName() +
66 " (km " + ref + ")");
67 segment.setReferencePoint(ref);
68 }
69
70 double [] values = segment.values;
71
72 if (numResults == -1) {
73 numResults = values.length;
74 }
75 else if (numResults != values.length) {
76 throw new IllegalArgumentException("wrong length of values");
77 }
78
79 // convert to Q if needed
80 if (!isQ && gauge != null) {
81
82 DischargeTable dt = gauge.fetchMasterDischargeTable();
83
84 double [][] table =
85 DischargeTables.loadDischargeTableValues(dt, 1);
86
87 // need the original values for naming
88 segment.backup();
89
90 for (int i = 0; i < values.length; ++i) {
91 double w = values[i] / 100.0;
92 double [] qs = DischargeTables.getQsForW(table, w);
93 if (qs.length == 0) {
94 logger.warn("No Qs found for W = " + values[i]);
95 addProblem("cannot.find.w.for.q", values[i]);
96 values[i] = Double.NaN;
97 }
98 else {
99 values[i] = qs[0];
100 if (qs.length > 1) {
101 logger.warn(
102 "More than one Q found for W = " + values[i]);
103 }
104 }
105 }
106 }
107 } // for all segments
108
109 Collections.sort(segments, REF_CMP);
110 }
111
112 public CalculationResult calculate(
113 WstValueTable table,
114 double from, double to, double step
115 ) {
116 boolean debug = logger.isDebugEnabled();
117
118 if (debug) {
119 logger.debug(
120 "calculate from " + from + " to " + to + " step " + step);
121 logger.debug("# segments: " + segments.size());
122 for (Segment segment: segments) {
123 logger.debug(" " + segment);
124 }
125 }
126
127 if (segments.isEmpty()) {
128 logger.debug("no segments found");
129 addProblem("no.segments.found");
130 return new CalculationResult(new WQKms[0], this);
131 }
132
133 int numResults = segments.get(0).values.length;
134
135 if (numResults < 1) {
136 logger.debug("no values given");
137 addProblem("no.values.given");
138 return new CalculationResult(new WQKms[0], this);
139 }
140
141
142 WQKms [] results = new WQKms[numResults];
143 for (int i = 0; i < results.length; ++i) {
144 results[i] = new WQKms();
145 }
146
147 if (Math.abs(step) < MINIMAL_STEP_WIDTH) {
148 step = MINIMAL_STEP_WIDTH;
149 }
150
151 if (from > to) {
152 step = -step;
153 }
154
155 QPosition [] qPositions = new QPosition[numResults];
156
157 Function [] functions = new Function[numResults];
158
159 double [] out = new double[2];
160
161 Segment sentinel = new Segment(Double.MAX_VALUE);
162 Segment s1 = sentinel, s2 = sentinel;
163
164 for (double pos = from;
165 from < to ? pos <= to : pos >= to;
166 pos = DoubleUtil.round(pos + step)
167 ) {
168 if (pos < s1.referencePoint || pos > s2.referencePoint) {
169 if (debug) {
170 logger.debug("need to find new interval for " + pos);
171 }
172 // find new interval
173 if (pos <= segments.get(0).referencePoint) {
174 // before first segment -> "gleichwertig"
175 if (debug) {
176 logger.debug("before first segment -> gleichwertig");
177 }
178 Segment first = segments.get(0);
179 double [] values = first.values;
180 double refPos = first.referencePoint;
181 for (int i = 0; i < qPositions.length; ++i) {
182 qPositions[i] = table.getQPosition(
183 refPos, values[i]);
184 }
185 sentinel.setReferencePoint(-Double.MAX_VALUE);
186 s1 = sentinel;
187 s2 = segments.get(0);
188 Arrays.fill(functions, Identity.IDENTITY);
189 }
190 else if (pos >= segments.get(segments.size()-1).referencePoint) {
191 // after last segment -> "gleichwertig"
192 if (debug) {
193 logger.debug("after last segment -> gleichwertig");
194 }
195 Segment last = segments.get(segments.size()-1);
196 double [] values = last.values;
197 double refPos = last.referencePoint;
198 for (int i = 0; i < qPositions.length; ++i) {
199 qPositions[i] = table.getQPosition(
200 refPos, values[i]);
201 }
202 sentinel.setReferencePoint(Double.MAX_VALUE);
203 s1 = last;
204 s2 = sentinel;
205 Arrays.fill(functions, Identity.IDENTITY);
206 }
207 else { // "ungleichwertig"
208 // find matching interval
209 if (debug) {
210 logger.debug("in segments -> ungleichwertig");
211 }
212 s1 = s2 = null;
213 for (int i = 1, N = segments.size(); i < N; ++i) {
214 Segment si1 = segments.get(i-1);
215 Segment si = segments.get(i);
216 if (debug) {
217 logger.debug("check " + pos + " in " +
218 si1.referencePoint + " - " + si.referencePoint);
219 }
220 if (pos >= si1.referencePoint
221 && pos <= si. referencePoint) {
222 s1 = si1;
223 s2 = si;
224 break;
225 }
226 }
227
228 if (s1 == null) {
229 throw new IllegalStateException("no interval found");
230 }
231
232 Segment anchor, free;
233
234 if (from > to) { anchor = s1; free = s2; }
235 else { anchor = s2; free = s1; }
236
237 // build transforms based on "gleichwertiger" phase
238 for (int i = 0; i < qPositions.length; ++i) {
239 QPosition qi = table.getQPosition(
240 anchor.referencePoint,
241 anchor.values[i]);
242
243 if ((qPositions[i] = qi) == null) {
244 addProblem(pos, "cannot.find.q", anchor.values[i]);
245 functions[i] = Identity.IDENTITY;
246 }
247 else {
248 double qA = table.getQ(qi, anchor.referencePoint);
249 double qF = table.getQ(qi, free .referencePoint);
250
251 functions[i] = Double.isNaN(qA) || Double.isNaN(qF)
252 ? Identity.IDENTITY
253 : new Linear(
254 qA, qF,
255 anchor.values[i], free.values[i]);
256
257 if (debug) {
258 logger.debug(
259 anchor.referencePoint + ": " +
260 qA + " -> " + functions[i].value(qA) +
261 " / " + free.referencePoint + ": " +
262 qF + " -> " + functions[i].value(qF));
263 }
264 }
265 } // build transforms
266 } // "ungleichwertiges" interval
267 } // find matching interval
268
269 for (int i = 0; i < qPositions.length; ++i) {
270 QPosition qPosition = qPositions[i];
271
272 if (qPosition == null) {
273 continue;
274 }
275
276 if (table.interpolate(pos, out, qPosition, functions[i])) {
277 results[i].add(out[0], out[1], pos);
278 }
279 else {
280 addProblem(pos, "cannot.interpolate.w.q");
281 }
282 }
283 }
284
285 // Backjump correction
286 for (int i = 0; i < results.length; ++i) {
287 BackJumpCorrector bjc = new BackJumpCorrector();
288
289 double [] ws = results[i].getWs();
290 double [] kms = results[i].getKms();
291
292 if (bjc.doCorrection(kms, ws, this)) {
293 results[i] = new WQCKms(results[i], bjc.getCorrected());
294 }
295 }
296
297 // name the curves
298 for (int i = 0; i < results.length; ++i) {
299 results[i].setName(createName(i));
300 }
301
302 return new CalculationResult(results, this);
303 }
304
305 protected String createName(int index) {
306 // TODO: I18N
307 StringBuilder sb = new StringBuilder(isQ ? "Q" : "W");
308 sb.append(" benutzerdefiniert (");
309 for (int i = 0, N = segments.size(); i < N; ++i) {
310 if (i > 0) {
311 sb.append("; ");
312 }
313 Segment segment = segments.get(i);
314 sb.append((segment.backup != null
315 ? segment.backup
316 : segment.values)[index]);
317 }
318 sb.append(')');
319 return sb.toString();
320 }
321 }
322 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
http://dive4elements.wald.intevation.org