[Dive4elements-commits] [PATCH] Some fixes to make the calculation work

[Wald Commits]

# HG changeset patch
# User Sascha L. Teichmann <teichmann at intevation.de>
# Date 1352067787 -3600
# Node ID 10e277c2fe0f800f170e849c285574af6fe64ceb
# Parent  5fb7efba81441ed5633939546da58330fe671f4f
Some fixes to make the calculation work.

diff -r 5fb7efba8144 -r 10e277c2fe0f flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/ExtremeCalculation.java
--- 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) {