Mercurial > dive4elements > river

comparison flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/Calculation4.java @ 3812:f788d2d901d6

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