--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/extreme/ExtremeCalculation.java	Thu Apr 25 15:23:37 2013 +0200
@@ -0,0 +1,418 @@
+package org.dive4elements.river.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 org.dive4elements.river.artifacts.access.ExtremeAccess;
+
+import org.dive4elements.river.artifacts.math.Linear;
+//import org.dive4elements.river.artifacts.math.Utils;
+
+import org.dive4elements.river.artifacts.math.fitting.Function;
+import org.dive4elements.river.artifacts.math.fitting.FunctionFactory;
+
+import org.dive4elements.river.artifacts.model.Calculation;
+import org.dive4elements.river.artifacts.model.CalculationResult;
+import org.dive4elements.river.artifacts.model.RangeWithValues;
+import org.dive4elements.river.artifacts.model.RiverFactory;
+import org.dive4elements.river.artifacts.model.WQKms;
+import org.dive4elements.river.artifacts.model.WstValueTable;
+import org.dive4elements.river.artifacts.model.WstValueTableFactory;
+
+import org.dive4elements.river.model.River;
+
+import org.dive4elements.river.utils.DoubleUtil;
+import org.dive4elements.river.utils.KMIndex;
+
+import gnu.trove.TDoubleArrayList;
+
+import java.util.List;
+
+import java.awt.geom.Line2D;
+
+/** Calculate extrapolated W. */
+public class ExtremeCalculation
+extends      Calculation
+{
+    private static final Log log =
+        LogFactory.getLog(ExtremeCalculation.class);
+
+    protected String                river;
+    protected String                function;
+    protected double                from;
+    protected double                to;
+    protected double                step;
+    protected double                percent;
+    protected List<RangeWithValues> ranges;
+
+    public ExtremeCalculation() {
+    }
+
+    public ExtremeCalculation(ExtremeAccess access) {
+        String                river    = access.getRiver();
+        String                function = access.getFunction();
+        Double                from     = access.getFrom();
+        Double                to       = access.getTo();
+        Double                step     = access.getStep();
+        Double                percent  = access.getPercent();
+        List<RangeWithValues> ranges   = access.getRanges();
+
+        if (river == null) {
+            // TODO: i18n
+            addProblem("extreme.no.river");
+        }
+
+        if (function == null) {
+            // TODO: i18n
+            addProblem("extreme.no.function");
+        }
+
+        if (from == null) {
+            // TODO: i18n
+            addProblem("extreme.no.from");
+        }
+
+        if (to == null) {
+            // TODO: i18n
+            addProblem("extreme.no.to");
+        }
+
+        if (step == null) {
+            // TODO: i18n
+            addProblem("extreme.no.step");
+        }
+
+        if (percent == null) {
+            // TODO: i18n
+            addProblem("extreme.no.percent");
+        }
+
+        if (ranges == null) {
+            // TODO: i18n
+            addProblem("extreme.no.ranges");
+        }
+
+        if (!hasProblems()) {
+            this.river    = river;
+            this.function = function;
+            this.from     = Math.min(from, to);
+            this.to       = Math.max(from, to);
+            this.step     = Math.max(0.001d, Math.abs(step)/1000d);
+            this.percent  = Math.max(0d, Math.min(100d, percent));
+            this.ranges   = ranges;
+        }
+    }
+
+
+    /** Calculate an extreme curve (extrapolate). */
+    public CalculationResult calculate() {
+
+        WstValueTable wst = null;
+
+        River river = RiverFactory.getRiver(this.river);
+        if (river == null) {
+            // TODO: i18n
+            addProblem("extreme.no.such.river", this.river);
+        }
+        else {
+            wst = WstValueTableFactory.getTable(river);
+            if (wst == null) {
+                // TODO: i18n
+                addProblem("extreme.no.wst.table");
+            }
+        }
+
+        Function function =
+            FunctionFactory.getInstance().getFunction(this.function);
+        if (function == null) {
+            // TODO: i18n
+            addProblem("extreme.no.such.function", this.function);
+        }
+
+        return hasProblems()
+            ? new CalculationResult(this)
+            : innerCalculate(wst, function);
+    }
+
+
+    /** Name of wqkms like W(5000,6000) */
+    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(
+        WstValueTable wst,
+        Function      function
+    ) {
+        RangeWithValues range = null;
+
+        double [] chiSqr = { 0d };
+
+        KMIndex<Curve> curves = new KMIndex<Curve>();
+        WQKms [] wqkms = allocWQKms();
+
+        boolean debug = log.isDebugEnabled();
+
+        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);
+            }
+
+            boolean foundRange = false;
+
+            if (range == null || !range.inside(km)) {
+                for (RangeWithValues r: ranges) {
+                    if (r.inside(km)) {
+                        range = r;
+                        foundRange = true;
+                        break;
+                    }
+                }
+                // TODO: i18n
+                if (!foundRange) {
+                    addProblem(km, "extreme.no.range.inner");
+                    continue;
+                }
+            }
+
+            double [][] wqs = wst.interpolateTabulated(km);
+            if (wqs == null) {
+                // TODO: i18n
+                addProblem(km, "extreme.no.raw.data");
+                continue;
+            }
+
+            // XXX: This should not be necessary for model data.
+            if (!DoubleUtil.isValid(wqs)) {
+                // TODO: i18n
+                addProblem(km, "extreme.invalid.data");
+                continue;
+            }
+
+            double [][] fitWQs = extractPointsToFit(wqs);
+            if (fitWQs == null) {
+                // TODO: i18n
+                addProblem(km, "extreme.too.less.points");
+                continue;
+            }
+
+            double [] coeffs = doFitting(function, fitWQs, chiSqr);
+            if (coeffs == null) {
+                // TODO: i18n
+                addProblem(km, "extreme.fitting.failed");
+                continue;
+            }
+
+            Curve curve = new Curve(
+                wqs[1], wqs[0],
+                function.getName(),
+                coeffs,
+                chiSqr[0]);
+
+            curves.add(km, 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(km, "extreme.evaluate.failed", values[i]);
+                }
+                else {
+                    wqkms[i].add(w, q, km);
+                }
+            }
+        }
+
+        ExtremeResult result = new ExtremeResult(curves, wqkms);
+        return new CalculationResult(result, this);
+    }
+
+    protected 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 = 0; i < ws.length; ++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) {
+            chiSqr[0] = lmo.getChiSquare();
+        }
+        return coeffs;
+    }
+
+    protected double [][] extractPointsToFit(double [][] wqs) {
+
+        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];
+
+        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 lastDir = -2;
+
+        for (int i = N-3; i >= 0; --i) {
+            double q = qs[i];
+            double w = ws[i];
+
+            if ((ascending && w > w1) || (!ascending && w < w1)) {
+                break;
+            }
+
+            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;
+            }
+            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) {
+            return;
+        }
+
+        double minQ = qs.getQuick(0);
+        double maxQ = qs.getQuick(N-1);
+        double factor = Math.min(Math.max(0d, percent/100d), 1d);
+        double cutQ = Linear.weight(factor, minQ, maxQ);
+        int cutIndex = 0;
+        for (; cutIndex < N; ++cutIndex) {
+            double q = qs.getQuick(cutIndex);
+            if (minQ < maxQ) {
+                if (q > cutQ) {
+                    break;
+                }
+            }
+            else {
+                if (q < cutQ) {
+                    break;
+                }
+            }
+        }
+        ws.remove(0, cutIndex);
+        qs.remove(0, cutIndex);
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf-8 :