Mercurial > dive4elements > river
changeset 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 | 5fb7efba8144 |
children | 0f93da769082 |
files | flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/ExtremeCalculation.java |
diffstat | 1 files changed, 64 insertions(+), 41 deletions(-) [+] |
line wrap: on
line diff
--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/ExtremeCalculation.java Fri Nov 02 17:48:53 2012 +0100 +++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/ExtremeCalculation.java Sun Nov 04 23:23:07 2012 +0100 @@ -12,7 +12,7 @@ import de.intevation.flys.artifacts.access.ExtremeAccess; import de.intevation.flys.artifacts.math.Linear; -import de.intevation.flys.artifacts.math.Utils; +//import de.intevation.flys.artifacts.math.Utils; import de.intevation.flys.artifacts.math.fitting.Function; import de.intevation.flys.artifacts.math.fitting.FunctionFactory; @@ -34,6 +34,8 @@ import java.util.List; +import java.awt.geom.Line2D; + /** Calculate extrapolated W. */ public class ExtremeCalculation extends Calculation @@ -185,46 +187,55 @@ KMIndex<Curve> curves = new KMIndex<Curve>(); WQKms [] wqkms = allocWQKms(); - for (double km = from; km <= to; km += step) { - double currentKm = DoubleUtil.round(km); + boolean debug = log.isDebugEnabled(); - if (range == null || !range.inside(currentKm)) { + from = DoubleUtil.round(from); + to = DoubleUtil.round(to); + + for (double km = from; km <= to; km = DoubleUtil.round(km+step)) { + + if (debug) { + log.debug("km: " + km); + } + + if (range == null || !range.inside(km)) { for (RangeWithValues r: ranges) { - if (r.inside(currentKm)) { + if (r.inside(km)) { range = r; break; } } // TODO: i18n - addProblem(currentKm, "extreme.no.range.inner"); + addProblem(km, "extreme.no.range.inner"); continue; } - double [][] wqs = wst.interpolateTabulated(currentKm); + double [][] wqs = wst.interpolateTabulated(km); if (wqs == null) { // TODO: i18n - addProblem(currentKm, "extreme.no.raw.data"); + addProblem(km, "extreme.no.raw.data"); continue; } // XXX: This should not be necessary for model data. if (!DoubleUtil.isValid(wqs)) { // TODO: i18n - addProblem(currentKm, "extreme.invalid.data"); + addProblem(km, "extreme.invalid.data"); continue; } double [][] fitWQs = extractPointsToFit(wqs); if (fitWQs == null) { // TODO: i18n - addProblem(currentKm, "extreme.too.less.points"); + addProblem(km, "extreme.too.less.points"); continue; } - double [] coeffs = doFitting(function, wqs, chiSqr); + double [] coeffs = doFitting(function, fitWQs, chiSqr); if (coeffs == null) { // TODO: i18n - addProblem(currentKm, "extreme.fitting.failed"); + addProblem(km, "extreme.fitting.failed"); + continue; } Curve curve = new Curve( @@ -233,7 +244,7 @@ coeffs, chiSqr[0]); - curves.add(currentKm, curve); + curves.add(km, curve); double [] values = range.getValues(); @@ -243,10 +254,10 @@ double w = curve.value(q); if (Double.isNaN(w)) { // TODO: i18n - addProblem(currentKm, "extreme.evaluate.failed", values[i]); + addProblem(km, "extreme.evaluate.failed", values[i]); } else { - wqkms[i].add(w, q, currentKm); + wqkms[i].add(w, q, km); } } } @@ -275,7 +286,7 @@ CurveFitter cf = new CurveFitter(lmo); - for (int i = ws.length-1; i >= 0; --i) { + for (int i = 0; i < ws.length; ++i) { cf.addObservedPoint(qs[i], ws[i]); } @@ -289,16 +300,13 @@ } } } - if (coeffs == null) { - return null; + if (coeffs != null) { + chiSqr[0] = lmo.getChiSquare(); } - chiSqr[0] = lmo.getChiSquare(); return coeffs; } protected double [][] extractPointsToFit(double [][] wqs) { - TDoubleArrayList ows = new TDoubleArrayList(); - TDoubleArrayList oqs = new TDoubleArrayList(); double [] ws = wqs[0]; double [] qs = wqs[1]; @@ -315,48 +323,63 @@ double q1 = qs[N-2]; double w1 = ws[N-2]; + boolean ascending = w2 > w1; + + TDoubleArrayList ows = new TDoubleArrayList(); + TDoubleArrayList oqs = new TDoubleArrayList(); + oqs.add(q2); oqs.add(q1); ows.add(w2); ows.add(w1); - int orientation = Utils.epsilonEquals(w1, w1) - ? 0 - : w1 > w2 - ? -1 - : +1; + int lastDir = -2; for (int i = N-3; i >= 0; --i) { - double cq = qs[i]; - double cw = qs[i]; - int newOrientation = Utils.relativeCCW( - q2, w2, - q1, w1, - cq, cw); + double q = qs[i]; + double w = ws[i]; - if (newOrientation != 0 && newOrientation != orientation) { + if ((ascending && w > w1) || (!ascending && w < w1)) { break; } - oqs.add(cq); - ows.add(cw); - if (newOrientation != 0) { - // rotate last out - q2 = q1; - w2 = w1; + int dir = Line2D.relativeCCW(q2, w2, q1, w1, q, w); + //int dir = Utils.relativeCCW(q2, w2, q1, w1, q, w); + if (lastDir == -2) { + lastDir = dir; } - q1 = cq; - w1 = cw; + else if (lastDir != dir) { + break; + } + + oqs.add(q); + ows.add(w); + w2 = w1; + q2 = q1; + w1 = w; + q1 = q; } oqs.reverse(); ows.reverse(); + + boolean debug = log.isDebugEnabled(); + if (debug) { + log.debug("from table: " + N); + log.debug("after trim: " + oqs.size()); + } + cutPercent(ows, oqs); + if (debug) { + log.debug("after percent cut: " + oqs.size()); + } + return new double [][] { ows.toNativeArray(), oqs.toNativeArray() }; } + protected void cutPercent(TDoubleArrayList ws, TDoubleArrayList qs) { int N = qs.size(); if (percent <= 0d || N == 0) {