Mercurial > dive4elements > river

comparison flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/ExtremeCalculation.java @ 4385:10e277c2fe0f flys-2.9.4

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Some fixes to make the calculation work.
author Sascha L. Teichmann <teichmann@intevation.de>
date Sun, 04 Nov 2012 23:23:07 +0100
parents d095b4267772
children 73200f5462fa
comparison
equal deleted inserted replaced
4384:5fb7efba8144 4385:10e277c2fe0f
10 import org.apache.commons.math.optimization.general.LevenbergMarquardtOptimizer; 10 import org.apache.commons.math.optimization.general.LevenbergMarquardtOptimizer;
11 11
12 import de.intevation.flys.artifacts.access.ExtremeAccess; 12 import de.intevation.flys.artifacts.access.ExtremeAccess;
13 13
14 import de.intevation.flys.artifacts.math.Linear; 14 import de.intevation.flys.artifacts.math.Linear;
15 import de.intevation.flys.artifacts.math.Utils; 15 //import de.intevation.flys.artifacts.math.Utils;
16 16
17 import de.intevation.flys.artifacts.math.fitting.Function; 17 import de.intevation.flys.artifacts.math.fitting.Function;
18 import de.intevation.flys.artifacts.math.fitting.FunctionFactory; 18 import de.intevation.flys.artifacts.math.fitting.FunctionFactory;
19 19
20 import de.intevation.flys.artifacts.model.Calculation; 20 import de.intevation.flys.artifacts.model.Calculation;
31 import de.intevation.flys.utils.KMIndex; 31 import de.intevation.flys.utils.KMIndex;
32 32
33 import gnu.trove.TDoubleArrayList; 33 import gnu.trove.TDoubleArrayList;
34 34
35 import java.util.List; 35 import java.util.List;
36
37 import java.awt.geom.Line2D;
36 38
37 /** Calculate extrapolated W. */ 39 /** Calculate extrapolated W. */
38 public class ExtremeCalculation 40 public class ExtremeCalculation
39 extends Calculation 41 extends Calculation
40 { 42 {
183 double [] chiSqr = { 0d }; 185 double [] chiSqr = { 0d };
184 186
185 KMIndex<Curve> curves = new KMIndex<Curve>(); 187 KMIndex<Curve> curves = new KMIndex<Curve>();
186 WQKms [] wqkms = allocWQKms(); 188 WQKms [] wqkms = allocWQKms();
187 189
188 for (double km = from; km <= to; km += step) { 190 boolean debug = log.isDebugEnabled();
189 double currentKm = DoubleUtil.round(km); 191
190 192 from = DoubleUtil.round(from);
191 if (range == null || !range.inside(currentKm)) { 193 to = DoubleUtil.round(to);
194
195 for (double km = from; km <= to; km = DoubleUtil.round(km+step)) {
196
197 if (debug) {
198 log.debug("km: " + km);
199 }
200
201 if (range == null || !range.inside(km)) {
192 for (RangeWithValues r: ranges) { 202 for (RangeWithValues r: ranges) {
193 if (r.inside(currentKm)) { 203 if (r.inside(km)) {
194 range = r; 204 range = r;
195 break; 205 break;
196 } 206 }
197 } 207 }
198 // TODO: i18n 208 // TODO: i18n
199 addProblem(currentKm, "extreme.no.range.inner"); 209 addProblem(km, "extreme.no.range.inner");
200 continue; 210 continue;
201 } 211 }
202 212
203 double [][] wqs = wst.interpolateTabulated(currentKm); 213 double [][] wqs = wst.interpolateTabulated(km);
204 if (wqs == null) { 214 if (wqs == null) {
205 // TODO: i18n 215 // TODO: i18n
206 addProblem(currentKm, "extreme.no.raw.data"); 216 addProblem(km, "extreme.no.raw.data");
207 continue; 217 continue;
208 } 218 }
209 219
210 // XXX: This should not be necessary for model data. 220 // XXX: This should not be necessary for model data.
211 if (!DoubleUtil.isValid(wqs)) { 221 if (!DoubleUtil.isValid(wqs)) {
212 // TODO: i18n 222 // TODO: i18n
213 addProblem(currentKm, "extreme.invalid.data"); 223 addProblem(km, "extreme.invalid.data");
214 continue; 224 continue;
215 } 225 }
216 226
217 double [][] fitWQs = extractPointsToFit(wqs); 227 double [][] fitWQs = extractPointsToFit(wqs);
218 if (fitWQs == null) { 228 if (fitWQs == null) {
219 // TODO: i18n 229 // TODO: i18n
220 addProblem(currentKm, "extreme.too.less.points"); 230 addProblem(km, "extreme.too.less.points");
221 continue; 231 continue;
222 } 232 }
223 233
224 double [] coeffs = doFitting(function, wqs, chiSqr); 234 double [] coeffs = doFitting(function, fitWQs, chiSqr);
225 if (coeffs == null) { 235 if (coeffs == null) {
226 // TODO: i18n 236 // TODO: i18n
227 addProblem(currentKm, "extreme.fitting.failed"); 237 addProblem(km, "extreme.fitting.failed");
238 continue;
228 } 239 }
229 240
230 Curve curve = new Curve( 241 Curve curve = new Curve(
231 wqs[1], wqs[0], 242 wqs[1], wqs[0],
232 function.getName(), 243 function.getName(),
233 coeffs, 244 coeffs,
234 chiSqr[0]); 245 chiSqr[0]);
235 246
236 curves.add(currentKm, curve); 247 curves.add(km, curve);
237 248
238 double [] values = range.getValues(); 249 double [] values = range.getValues();
239 250
240 int V = Math.min(values.length, wqkms.length); 251 int V = Math.min(values.length, wqkms.length);
241 for (int i = 0; i < V; ++i) { 252 for (int i = 0; i < V; ++i) {
242 double q = values[i]; 253 double q = values[i];
243 double w = curve.value(q); 254 double w = curve.value(q);
244 if (Double.isNaN(w)) { 255 if (Double.isNaN(w)) {
245 // TODO: i18n 256 // TODO: i18n
246 addProblem(currentKm, "extreme.evaluate.failed", values[i]); 257 addProblem(km, "extreme.evaluate.failed", values[i]);
247 } 258 }
248 else { 259 else {
249 wqkms[i].add(w, q, currentKm); 260 wqkms[i].add(w, q, km);
250 } 261 }
251 } 262 }
252 } 263 }
253 264
254 ExtremeResult result = new ExtremeResult(curves, wqkms); 265 ExtremeResult result = new ExtremeResult(curves, wqkms);
273 lmo.setOrthoTolerance(tolerance); 284 lmo.setOrthoTolerance(tolerance);
274 lmo.setParRelativeTolerance(tolerance); 285 lmo.setParRelativeTolerance(tolerance);
275 286
276 CurveFitter cf = new CurveFitter(lmo); 287 CurveFitter cf = new CurveFitter(lmo);
277 288
278 for (int i = ws.length-1; i >= 0; --i) { 289 for (int i = 0; i < ws.length; ++i) {
279 cf.addObservedPoint(qs[i], ws[i]); 290 cf.addObservedPoint(qs[i], ws[i]);
280 } 291 }
281 292
282 try { 293 try {
283 coeffs = cf.fit(function, function.getInitialGuess()); 294 coeffs = cf.fit(function, function.getInitialGuess());
287 if (log.isDebugEnabled()) { 298 if (log.isDebugEnabled()) {
288 log.debug("tolerance " + tolerance + " + failed."); 299 log.debug("tolerance " + tolerance + " + failed.");
289 } 300 }
290 } 301 }
291 } 302 }
292 if (coeffs == null) { 303 if (coeffs != null) {
293 return null; 304 chiSqr[0] = lmo.getChiSquare();
294 } 305 }
295 chiSqr[0] = lmo.getChiSquare();
296 return coeffs; 306 return coeffs;
297 } 307 }
298 308
299 protected double [][] extractPointsToFit(double [][] wqs) { 309 protected double [][] extractPointsToFit(double [][] wqs) {
300 TDoubleArrayList ows = new TDoubleArrayList();
301 TDoubleArrayList oqs = new TDoubleArrayList();
302 310
303 double [] ws = wqs[0]; 311 double [] ws = wqs[0];
304 double [] qs = wqs[1]; 312 double [] qs = wqs[1];
305 313
306 int N = Math.min(ws.length, qs.length); 314 int N = Math.min(ws.length, qs.length);
313 double q2 = qs[N-1]; 321 double q2 = qs[N-1];
314 double w2 = ws[N-1]; 322 double w2 = ws[N-1];
315 double q1 = qs[N-2]; 323 double q1 = qs[N-2];
316 double w1 = ws[N-2]; 324 double w1 = ws[N-2];
317 325
326 boolean ascending = w2 > w1;
327
328 TDoubleArrayList ows = new TDoubleArrayList();
329 TDoubleArrayList oqs = new TDoubleArrayList();
330
318 oqs.add(q2); oqs.add(q1); 331 oqs.add(q2); oqs.add(q1);
319 ows.add(w2); ows.add(w1); 332 ows.add(w2); ows.add(w1);
320 333
321 int orientation = Utils.epsilonEquals(w1, w1) 334 int lastDir = -2;
322 ? 0
323 : w1 > w2
324 ? -1
325 : +1;
326 335
327 for (int i = N-3; i >= 0; --i) { 336 for (int i = N-3; i >= 0; --i) {
328 double cq = qs[i]; 337 double q = qs[i];
329 double cw = qs[i]; 338 double w = ws[i];
330 int newOrientation = Utils.relativeCCW( 339
331 q2, w2, 340 if ((ascending && w > w1) || (!ascending && w < w1)) {
332 q1, w1,
333 cq, cw);
334
335 if (newOrientation != 0 && newOrientation != orientation) {
336 break; 341 break;
337 } 342 }
338 oqs.add(cq); 343
339 ows.add(cw); 344 int dir = Line2D.relativeCCW(q2, w2, q1, w1, q, w);
340 345 //int dir = Utils.relativeCCW(q2, w2, q1, w1, q, w);
341 if (newOrientation != 0) { 346 if (lastDir == -2) {
342 // rotate last out 347 lastDir = dir;
343 q2 = q1; 348 }
344 w2 = w1; 349 else if (lastDir != dir) {
345 } 350 break;
346 q1 = cq; 351 }
347 w1 = cw; 352
353 oqs.add(q);
354 ows.add(w);
355 w2 = w1;
356 q2 = q1;
357 w1 = w;
358 q1 = q;
348 } 359 }
349 360
350 oqs.reverse(); 361 oqs.reverse();
351 ows.reverse(); 362 ows.reverse();
363
364 boolean debug = log.isDebugEnabled();
365 if (debug) {
366 log.debug("from table: " + N);
367 log.debug("after trim: " + oqs.size());
368 }
369
352 cutPercent(ows, oqs); 370 cutPercent(ows, oqs);
371
372 if (debug) {
373 log.debug("after percent cut: " + oqs.size());
374 }
353 375
354 return new double [][] { 376 return new double [][] {
355 ows.toNativeArray(), 377 ows.toNativeArray(),
356 oqs.toNativeArray() 378 oqs.toNativeArray()
357 }; 379 };
358 } 380 }
381
359 382
360 protected void cutPercent(TDoubleArrayList ws, TDoubleArrayList qs) { 383 protected void cutPercent(TDoubleArrayList ws, TDoubleArrayList qs) {
361 int N = qs.size(); 384 int N = qs.size();
362 if (percent <= 0d || N == 0) { 385 if (percent <= 0d || N == 0) {
363 return; 386 return;
http://dive4elements.wald.intevation.org