Mercurial > dive4elements > river
comparison flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/ExtremeCalculation.java @ 4385:10e277c2fe0f flys-2.9.4
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; |