--- 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 :