# HG changeset patch
# User Sascha L. Teichmann <teichmann@intevation.de>
# Date 1351178737 -7200
# Node ID 5cc9453456a7ee7bdfffab3d7aa004e858ee6db7
# Parent  2c6e571f366ae3cd488549cd67b1aacab3f79269
First complete but untested version of the 'Auslagerung extremer Wasserspiegellagen' calculation.

diff -r 2c6e571f366a -r 5cc9453456a7 flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/Utils.java
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/Utils.java	Thu Oct 25 17:25:37 2012 +0200
@@ -0,0 +1,61 @@
+package de.intevation.flys.artifacts.math;
+
+
+public final class Utils {
+
+    public static final double EPSILON = 1e-3;
+
+    private Utils() {
+    }
+
+    public static final boolean epsilonEquals(double a, double b) {
+        return epsilonEquals(a, b, EPSILON);
+    }
+
+    public static final boolean epsilonEquals(double a, double b, double eps) {
+        return Math.abs(a - b) < eps;
+    }
+
+    public static int relativeCCW(
+        double x1, double y1,
+        double x2, double y2,
+        double px, double py
+    ) {
+        if ((epsilonEquals(x1, x2) && epsilonEquals(y1, y2))
+        || ((epsilonEquals(x1, px) && epsilonEquals(y1, py)))) {
+            return 0; // Coincident points.
+        }
+        // Translate to the origin.
+        x2 -= x1;
+        y2 -= y1;
+        px -= x1;
+        py -= y1;
+        double slope2 = y2 / x2;
+        double slopep = py / px;
+        if (epsilonEquals(slope2, slopep) 
+        || (epsilonEquals(x2, 0.0) && epsilonEquals(px, 0.0))) {
+            return y2 > EPSILON // Colinear.
+                ? (py < -EPSILON ? -1 : py > y2 ? 1 : 0)
+                : (py > -EPSILON ? -1 : py < y2 ? 1 : 0);
+        }
+        if (x2 >= EPSILON && slope2 >= EPSILON) {
+            return px >= EPSILON // Quadrant 1.
+                ? (slope2 > slopep ? 1 : -1)
+                : (slope2 < slopep ? 1 : -1);
+        }
+
+        if (y2 > EPSILON) {
+            return px < -EPSILON // Quadrant 2.
+                ? (slope2 > slopep ? 1 : -1)
+                : (slope2 < slopep ? 1 : -1);
+        }
+        if (slope2 >= EPSILON) {
+            return px >= EPSILON // Quadrant 3.
+                ? (slope2 < slopep ? 1 : -1)
+                : (slope2 > slopep ? 1 : -1);
+        }
+        return px < -EPSILON // Quadrant 4.
+            ? (slope2 < slopep ? 1 : -1)
+            : (slope2 > slopep ? 1 : -1);
+    }
+}
diff -r 2c6e571f366a -r 5cc9453456a7 flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/Curve.java
--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/Curve.java	Thu Oct 25 15:17:01 2012 +0200
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/Curve.java	Thu Oct 25 17:25:37 2012 +0200
@@ -25,6 +25,7 @@
     protected double [] ws;
     protected String    function;
     protected double [] coeffs;
+    protected double    chiSqr;
 
     // The spline is pretty heavy weight so cache it with a soft ref only.
     protected transient SoftReference<Function> spline;
@@ -37,7 +38,8 @@
         double [] qs,
         double [] ws,
         String    function,
-        double [] coeffs
+        double [] coeffs,
+        double    chiSqr
     ) {
         this.qs       = qs;
         this.ws       = ws;
@@ -81,6 +83,24 @@
         return extrapolation;
     }
 
+    /**
+     * Gets the chiSqr for this instance.
+     *
+     * @return The chiSqr.
+     */
+    public double getChiSqr() {
+        return this.chiSqr;
+    }
+
+    /**
+     * Sets the chiSqr for this instance.
+     *
+     * @param chiSqr The chiSqr.
+     */
+    public void setChiSqr(double chiSqr) {
+        this.chiSqr = chiSqr;
+    }
+
     protected synchronized Function getSpline() {
         Function sp;
         if (spline != null) {
diff -r 2c6e571f366a -r 5cc9453456a7 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	Thu Oct 25 15:17:01 2012 +0200
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/ExtremeCalculation.java	Thu Oct 25 17:25:37 2012 +0200
@@ -1,7 +1,18 @@
 package de.intevation.flys.artifacts.model.extreme;
 
+import org.apache.commons.logging.Log;
+import org.apache.commons.logging.LogFactory;
+
+import org.apache.commons.math.MathException;
+
+import org.apache.commons.math.optimization.fitting.CurveFitter;
+
+import org.apache.commons.math.optimization.general.LevenbergMarquardtOptimizer;
+
 import de.intevation.flys.artifacts.access.ExtremeAccess;
 
+import de.intevation.flys.artifacts.math.Utils;
+
 import de.intevation.flys.artifacts.math.fitting.Function;
 import de.intevation.flys.artifacts.math.fitting.FunctionFactory;
 
@@ -9,12 +20,16 @@
 import de.intevation.flys.artifacts.model.CalculationResult;
 import de.intevation.flys.artifacts.model.RangeWithValues;
 import de.intevation.flys.artifacts.model.RiverFactory;
+import de.intevation.flys.artifacts.model.WQKms;
 import de.intevation.flys.artifacts.model.WstValueTable;
 import de.intevation.flys.artifacts.model.WstValueTableFactory;
 
 import de.intevation.flys.model.River;
 
 import de.intevation.flys.utils.DoubleUtil;
+import de.intevation.flys.utils.KMIndex;
+
+import gnu.trove.TDoubleArrayList;
 
 import java.util.List;
 
@@ -22,6 +37,9 @@
 public class ExtremeCalculation
 extends      Calculation
 {
+    private static final Log log =
+        LogFactory.getLog(ExtremeCalculation.class);
+
     protected String                river;
     protected String                function;
     protected double                from;
@@ -119,6 +137,40 @@
             : innerCalculate(wst, function);
     }
 
+    protected String wqkmsName(int i) {
+        StringBuilder sb = new StringBuilder("W(");
+        boolean already = false;
+        for (RangeWithValues r: ranges) {
+            double [] values = r.getValues();
+            if (i < values.length) {
+                if (already) {
+                    sb.append(", ");
+                }
+                else {
+                    already = true;
+                }
+                // TODO: i18n
+                sb.append(values[i]);
+            }
+        }
+        return sb.append(')').toString();
+    }
+
+    protected WQKms [] allocWQKms() {
+        int max = 0;
+        for (RangeWithValues r: ranges) {
+            double [] values = r.getValues();
+            if (values.length > max) {
+                max = values.length;
+            }
+        }
+        WQKms [] wqkms = new WQKms[max];
+        for (int i = 0; i < max; ++i) {
+            wqkms[i] = new WQKms(wqkmsName(i));
+        }
+        return wqkms;
+    }
+
 
     /** Calculate an extreme curve (extrapolate). */
     protected CalculationResult innerCalculate(
@@ -127,7 +179,12 @@
     ) {
         RangeWithValues range = null;
 
-        KMs: for (double km = from; km <= to; km += step) {
+        double [] chiSqr = { 0d };
+
+        KMIndex<Curve> curves = new KMIndex<Curve>();
+        WQKms [] wqkms = allocWQKms();
+
+        for (double km = from; km <= to; km += step) {
             double currentKm = DoubleUtil.round(km);
 
             if (range == null || !range.inside(currentKm)) {
@@ -139,7 +196,7 @@
                 }
                 // TODO: i18n
                 addProblem(currentKm, "extreme.no.range");
-                continue KMs;
+                continue;
             }
 
             double [][] wqs = wst.interpolateTabulated(currentKm);
@@ -155,13 +212,148 @@
                 addProblem(currentKm, "extreme.invalid.data");
                 continue;
             }
-            // TODO: Implement extraction of points for curve fitting.
-            // TODO: Implement curve fitting.
-            // TODO: Implement generating Curve object per km.
+
+            double [][] fitWQs = extractPointsToFit(wqs);
+            if (fitWQs == null) {
+                // TODO: i18n
+                addProblem(currentKm, "extreme.too.less.points");
+                continue;
+            }
+
+            double [] coeffs = doFitting(function, wqs, chiSqr);
+            if (coeffs == null) {
+                // TODO: i18n
+                addProblem(currentKm, "extreme.fitting.failed");
+            }
+
+            Curve curve = new Curve(
+                wqs[1], wqs[0],
+                function.getName(),
+                coeffs,
+                chiSqr[0]);
+
+            curves.add(currentKm, curve);
+
+            double [] values = range.getValues();
+
+            int V = Math.min(values.length, wqkms.length);
+            for (int i = 0; i < V; ++i) {
+                double q = values[i];
+                double w = curve.value(q);
+                if (Double.isNaN(w)) {
+                    // TODO: i18n
+                    addProblem(currentKm, "extreme.evaluate.failed", values[i]);
+                }
+                else {
+                    wqkms[i].add(w, q);
+                }
+            }
         }
 
-        ExtremeResult result = new ExtremeResult();
+        ExtremeResult result = new ExtremeResult(curves, wqkms);
         return new CalculationResult(result, this);
     }
+
+    protected static double [] doFitting(
+        Function    function,
+        double [][] wqs,
+        double []   chiSqr
+    ) {
+        LevenbergMarquardtOptimizer lmo = null;
+
+        double [] coeffs = null;
+
+        double [] ws = wqs[0];
+        double [] qs = wqs[1];
+
+        for (double tolerance = 1e-10; tolerance < 1e-3; tolerance *= 10d) {
+            lmo = new LevenbergMarquardtOptimizer();
+            lmo.setCostRelativeTolerance(tolerance);
+            lmo.setOrthoTolerance(tolerance);
+            lmo.setParRelativeTolerance(tolerance);
+
+            CurveFitter cf = new CurveFitter(lmo);
+
+            for (int i = ws.length-1; i >= 0; --i) {
+                cf.addObservedPoint(qs[i], ws[i]);
+            }
+
+            try {
+                coeffs = cf.fit(function, function.getInitialGuess());
+                break;
+            }
+            catch (MathException me) {
+                if (log.isDebugEnabled()) {
+                    log.debug("tolerance " + tolerance + " + failed.");
+                }
+            }
+        }
+        if (coeffs == null) {
+            return null;
+        }
+        chiSqr[0] = lmo.getChiSquare();
+        return coeffs;
+    }
+
+    protected static double [][] extractPointsToFit(double [][] wqs) {
+        TDoubleArrayList ows = new TDoubleArrayList();
+        TDoubleArrayList oqs = new TDoubleArrayList();
+
+        double [] ws = wqs[0];
+        double [] qs = wqs[1];
+
+        int N = Math.min(ws.length, qs.length);
+
+        if (N < 2) {
+            log.warn("Too less points for fitting");
+            return null;
+        }
+
+        double q2 = qs[N-1];
+        double w2 = ws[N-1];
+        double q1 = qs[N-2];
+        double w1 = ws[N-2];
+
+        oqs.add(q2); oqs.add(q1);
+        ows.add(w2); ows.add(w1);
+
+        int orientation = Utils.epsilonEquals(w1, w1)
+            ? 0
+            : w1 > w2
+                ? -1
+                : +1;
+
+        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);
+
+            if (newOrientation != 0 && newOrientation != orientation) {
+                break;
+            }
+            oqs.add(cq);
+            ows.add(cw);
+
+            if (newOrientation != 0) {
+                // rotate last out
+                q2 = q1;
+                w2 = w1;
+            }
+            q1 = cq;
+            w1 = cw;
+        }
+
+        // XXX: Not really needed for fitting.
+        // oqs.reverse();
+        // ows.reverse();
+
+        return new double [][] {
+            ows.toNativeArray(),
+            oqs.toNativeArray()
+        };
+    }
 }
 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf-8 :
diff -r 2c6e571f366a -r 5cc9453456a7 flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/ExtremeResult.java
--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/ExtremeResult.java	Thu Oct 25 15:17:01 2012 +0200
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/ExtremeResult.java	Thu Oct 25 17:25:37 2012 +0200
@@ -2,10 +2,78 @@
 
 import java.io.Serializable;
 
+import de.intevation.flys.artifacts.model.WQKms;
+
+import de.intevation.flys.utils.KMIndex;
+
 public class ExtremeResult
 implements   Serializable
 {
+    protected KMIndex<Curve> curves;
+    protected WQKms [] wqkms;
+
     public ExtremeResult() {
     }
+
+    public ExtremeResult(KMIndex<Curve> curves, WQKms [] wqkms) {
+        this.curves = curves;
+        this.wqkms = wqkms;
+    }
+
+    /**
+     * Gets the curves for this instance.
+     *
+     * @return The curves.
+     */
+    public KMIndex<Curve> getCurves() {
+        return this.curves;
+    }
+
+    /**
+     * Sets the curves for this instance.
+     *
+     * @param curves The curves.
+     */
+    public void setCurves(KMIndex<Curve> curves) {
+        this.curves = curves;
+    }
+
+    /**
+     * Gets the wqkms for this instance.
+     *
+     * @return The wqkms.
+     */
+    public WQKms[] getWQKms() {
+        return this.wqkms;
+    }
+
+    /**
+     * Gets the wqkms for this instance.
+     *
+     * @param index The index to get.
+     * @return The wqkms.
+     */
+    public WQKms getWQKms(int index) {
+        return this.wqkms[index];
+    }
+
+    /**
+     * Sets the wqkms for this instance.
+     *
+     * @param wqkms The wqkms.
+     */
+    public void setWQKms(WQKms[] wqkms) {
+        this.wqkms = wqkms;
+    }
+
+    /**
+     * Sets the wqkms for this instance.
+     *
+     * @param index The index to set.
+     * @param wqkms The wqkms.
+     */
+    public void setWQKms(int index, WQKms wqkms) {
+        this.wqkms[index] = wqkms;
+    }
 }
 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :