sascha@3775: package de.intevation.flys.artifacts.model.extreme;
sascha@3775: 
teichmann@4259: import org.apache.commons.logging.Log;
teichmann@4259: import org.apache.commons.logging.LogFactory;
teichmann@4259: 
teichmann@4259: import org.apache.commons.math.MathException;
teichmann@4259: 
teichmann@4259: import org.apache.commons.math.optimization.fitting.CurveFitter;
teichmann@4259: 
teichmann@4259: import org.apache.commons.math.optimization.general.LevenbergMarquardtOptimizer;
teichmann@4259: 
sascha@3775: import de.intevation.flys.artifacts.access.ExtremeAccess;
sascha@3775: 
teichmann@4355: import de.intevation.flys.artifacts.math.Linear;
teichmann@4385: //import de.intevation.flys.artifacts.math.Utils;
teichmann@4259: 
ingo@3785: import de.intevation.flys.artifacts.math.fitting.Function;
ingo@3785: import de.intevation.flys.artifacts.math.fitting.FunctionFactory;
ingo@3785: 
sascha@3775: import de.intevation.flys.artifacts.model.Calculation;
sascha@3775: import de.intevation.flys.artifacts.model.CalculationResult;
ingo@3785: import de.intevation.flys.artifacts.model.RangeWithValues;
ingo@3785: import de.intevation.flys.artifacts.model.RiverFactory;
teichmann@4259: import de.intevation.flys.artifacts.model.WQKms;
ingo@3785: import de.intevation.flys.artifacts.model.WstValueTable;
ingo@3785: import de.intevation.flys.artifacts.model.WstValueTableFactory;
ingo@3785: 
ingo@3785: import de.intevation.flys.model.River;
ingo@3785: 
ingo@3785: import de.intevation.flys.utils.DoubleUtil;
teichmann@4259: import de.intevation.flys.utils.KMIndex;
teichmann@4259: 
teichmann@4259: import gnu.trove.TDoubleArrayList;
ingo@3785: 
ingo@3785: import java.util.List;
sascha@3775: 
teichmann@4385: import java.awt.geom.Line2D;
teichmann@4385: 
felix@4040: /** Calculate extrapolated W. */
sascha@3775: public class ExtremeCalculation
sascha@3775: extends      Calculation
sascha@3775: {
teichmann@4259:     private static final Log log =
teichmann@4259:         LogFactory.getLog(ExtremeCalculation.class);
teichmann@4259: 
ingo@3785:     protected String                river;
ingo@3785:     protected String                function;
ingo@3785:     protected double                from;
ingo@3785:     protected double                to;
ingo@3785:     protected double                step;
ingo@3785:     protected double                percent;
ingo@3785:     protected List<RangeWithValues> ranges;
ingo@3785: 
sascha@3775:     public ExtremeCalculation() {
sascha@3775:     }
sascha@3775: 
sascha@3775:     public ExtremeCalculation(ExtremeAccess access) {
ingo@3785:         String                river    = access.getRiver();
ingo@3785:         String                function = access.getFunction();
ingo@3785:         Double                from     = access.getFrom();
ingo@3785:         Double                to       = access.getTo();
ingo@3785:         Double                step     = access.getStep();
ingo@3785:         Double                percent  = access.getPercent();
ingo@3785:         List<RangeWithValues> ranges   = access.getRanges();
ingo@3785: 
ingo@3785:         if (river == null) {
ingo@3785:             // TODO: i18n
ingo@3785:             addProblem("extreme.no.river");
ingo@3785:         }
ingo@3785: 
ingo@3785:         if (function == null) {
ingo@3785:             // TODO: i18n
ingo@3785:             addProblem("extreme.no.function");
ingo@3785:         }
ingo@3785: 
ingo@3785:         if (from == null) {
ingo@3785:             // TODO: i18n
ingo@3785:             addProblem("extreme.no.from");
ingo@3785:         }
ingo@3785: 
ingo@3785:         if (to == null) {
ingo@3785:             // TODO: i18n
ingo@3785:             addProblem("extreme.no.to");
ingo@3785:         }
ingo@3785: 
ingo@3785:         if (step == null) {
ingo@3785:             // TODO: i18n
ingo@3785:             addProblem("extreme.no.step");
ingo@3785:         }
ingo@3785: 
ingo@3785:         if (percent == null) {
ingo@3785:             // TODO: i18n
ingo@3785:             addProblem("extreme.no.percent");
ingo@3785:         }
ingo@3785: 
ingo@3785:         if (ranges == null) {
ingo@3785:             // TODO: i18n
ingo@3785:             addProblem("extreme.no.ranges");
ingo@3785:         }
ingo@3785: 
ingo@3785:         if (!hasProblems()) {
ingo@3785:             this.river    = river;
ingo@3785:             this.function = function;
ingo@3785:             this.from     = Math.min(from, to);
ingo@3785:             this.to       = Math.max(from, to);
ingo@3785:             this.step     = Math.max(0.001d, Math.abs(step)/1000d);
ingo@3785:             this.percent  = Math.max(0d, Math.min(100d, percent));
ingo@3785:             this.ranges   = ranges;
ingo@3785:         }
sascha@3775:     }
sascha@3775: 
felix@4059: 
felix@4059:     /** Calculate an extreme curve (extrapolate). */
sascha@3775:     public CalculationResult calculate() {
ingo@3785: 
ingo@3785:         WstValueTable wst = null;
ingo@3785: 
ingo@3785:         River river = RiverFactory.getRiver(this.river);
ingo@3785:         if (river == null) {
ingo@3785:             // TODO: i18n
ingo@3785:             addProblem("extreme.no.such.river", this.river);
sascha@3775:         }
ingo@3785:         else {
ingo@3785:             wst = WstValueTableFactory.getTable(river);
ingo@3785:             if (wst == null) {
ingo@3785:                 // TODO: i18n
ingo@3785:                 addProblem("extreme.no.wst.table");
ingo@3785:             }
ingo@3785:         }
ingo@3785: 
ingo@3785:         Function function =
ingo@3785:             FunctionFactory.getInstance().getFunction(this.function);
ingo@3785:         if (function == null) {
ingo@3785:             // TODO: i18n
ingo@3785:             addProblem("extreme.no.such.function", this.function);
ingo@3785:         }
ingo@3785: 
ingo@3785:         return hasProblems()
ingo@3785:             ? new CalculationResult(this)
ingo@3785:             : innerCalculate(wst, function);
ingo@3785:     }
ingo@3785: 
teichmann@4259:     protected String wqkmsName(int i) {
teichmann@4259:         StringBuilder sb = new StringBuilder("W(");
teichmann@4259:         boolean already = false;
teichmann@4259:         for (RangeWithValues r: ranges) {
teichmann@4259:             double [] values = r.getValues();
teichmann@4259:             if (i < values.length) {
teichmann@4259:                 if (already) {
teichmann@4259:                     sb.append(", ");
teichmann@4259:                 }
teichmann@4259:                 else {
teichmann@4259:                     already = true;
teichmann@4259:                 }
teichmann@4259:                 // TODO: i18n
teichmann@4259:                 sb.append(values[i]);
teichmann@4259:             }
teichmann@4259:         }
teichmann@4259:         return sb.append(')').toString();
teichmann@4259:     }
teichmann@4259: 
teichmann@4259:     protected WQKms [] allocWQKms() {
teichmann@4259:         int max = 0;
teichmann@4259:         for (RangeWithValues r: ranges) {
teichmann@4259:             double [] values = r.getValues();
teichmann@4259:             if (values.length > max) {
teichmann@4259:                 max = values.length;
teichmann@4259:             }
teichmann@4259:         }
teichmann@4259:         WQKms [] wqkms = new WQKms[max];
teichmann@4259:         for (int i = 0; i < max; ++i) {
teichmann@4259:             wqkms[i] = new WQKms(wqkmsName(i));
teichmann@4259:         }
teichmann@4259:         return wqkms;
teichmann@4259:     }
teichmann@4259: 
felix@4059: 
felix@4059:     /** Calculate an extreme curve (extrapolate). */
ingo@3785:     protected CalculationResult innerCalculate(
ingo@3785:         WstValueTable wst,
ingo@3785:         Function      function
ingo@3785:     ) {
ingo@3785:         RangeWithValues range = null;
ingo@3785: 
teichmann@4259:         double [] chiSqr = { 0d };
teichmann@4259: 
teichmann@4259:         KMIndex<Curve> curves = new KMIndex<Curve>();
teichmann@4259:         WQKms [] wqkms = allocWQKms();
teichmann@4259: 
teichmann@4385:         boolean debug = log.isDebugEnabled();
ingo@3785: 
teichmann@4385:         from = DoubleUtil.round(from);
teichmann@4385:         to = DoubleUtil.round(to);
teichmann@4385: 
teichmann@4385:         for (double km = from; km <= to; km = DoubleUtil.round(km+step)) {
teichmann@4385: 
teichmann@4385:             if (debug) {
teichmann@4385:                 log.debug("km: " + km);
teichmann@4385:             }
teichmann@4385: 
teichmann@4385:             if (range == null || !range.inside(km)) {
ingo@3785:                 for (RangeWithValues r: ranges) {
teichmann@4385:                     if (r.inside(km)) {
ingo@3785:                         range = r;
ingo@3785:                         break;
ingo@3785:                     }
ingo@3785:                 }
ingo@3785:                 // TODO: i18n
teichmann@4385:                 addProblem(km, "extreme.no.range.inner");
teichmann@4259:                 continue;
ingo@3785:             }
ingo@3785: 
teichmann@4385:             double [][] wqs = wst.interpolateTabulated(km);
ingo@3785:             if (wqs == null) {
ingo@3785:                 // TODO: i18n
teichmann@4385:                 addProblem(km, "extreme.no.raw.data");
ingo@3785:                 continue;
ingo@3785:             }
ingo@3785: 
ingo@3785:             // XXX: This should not be necessary for model data.
ingo@3785:             if (!DoubleUtil.isValid(wqs)) {
ingo@3785:                 // TODO: i18n
teichmann@4385:                 addProblem(km, "extreme.invalid.data");
ingo@3785:                 continue;
ingo@3785:             }
teichmann@4259: 
teichmann@4259:             double [][] fitWQs = extractPointsToFit(wqs);
teichmann@4259:             if (fitWQs == null) {
teichmann@4259:                 // TODO: i18n
teichmann@4385:                 addProblem(km, "extreme.too.less.points");
teichmann@4259:                 continue;
teichmann@4259:             }
teichmann@4259: 
teichmann@4385:             double [] coeffs = doFitting(function, fitWQs, chiSqr);
teichmann@4259:             if (coeffs == null) {
teichmann@4259:                 // TODO: i18n
teichmann@4385:                 addProblem(km, "extreme.fitting.failed");
teichmann@4385:                 continue;
teichmann@4259:             }
teichmann@4259: 
teichmann@4259:             Curve curve = new Curve(
teichmann@4259:                 wqs[1], wqs[0],
teichmann@4259:                 function.getName(),
teichmann@4259:                 coeffs,
teichmann@4259:                 chiSqr[0]);
teichmann@4259: 
teichmann@4385:             curves.add(km, curve);
teichmann@4259: 
teichmann@4259:             double [] values = range.getValues();
teichmann@4259: 
teichmann@4259:             int V = Math.min(values.length, wqkms.length);
teichmann@4259:             for (int i = 0; i < V; ++i) {
teichmann@4259:                 double q = values[i];
teichmann@4259:                 double w = curve.value(q);
teichmann@4259:                 if (Double.isNaN(w)) {
teichmann@4259:                     // TODO: i18n
teichmann@4385:                     addProblem(km, "extreme.evaluate.failed", values[i]);
teichmann@4259:                 }
teichmann@4259:                 else {
teichmann@4385:                     wqkms[i].add(w, q, km);
teichmann@4259:                 }
teichmann@4259:             }
ingo@3785:         }
ingo@3785: 
teichmann@4259:         ExtremeResult result = new ExtremeResult(curves, wqkms);
sascha@3775:         return new CalculationResult(result, this);
sascha@3775:     }
teichmann@4259: 
teichmann@4355:     protected double [] doFitting(
teichmann@4259:         Function    function,
teichmann@4259:         double [][] wqs,
teichmann@4259:         double []   chiSqr
teichmann@4259:     ) {
teichmann@4259:         LevenbergMarquardtOptimizer lmo = null;
teichmann@4259: 
teichmann@4259:         double [] coeffs = null;
teichmann@4259: 
teichmann@4259:         double [] ws = wqs[0];
teichmann@4259:         double [] qs = wqs[1];
teichmann@4259: 
teichmann@4259:         for (double tolerance = 1e-10; tolerance < 1e-3; tolerance *= 10d) {
teichmann@4259:             lmo = new LevenbergMarquardtOptimizer();
teichmann@4259:             lmo.setCostRelativeTolerance(tolerance);
teichmann@4259:             lmo.setOrthoTolerance(tolerance);
teichmann@4259:             lmo.setParRelativeTolerance(tolerance);
teichmann@4259: 
teichmann@4259:             CurveFitter cf = new CurveFitter(lmo);
teichmann@4259: 
teichmann@4385:             for (int i = 0; i < ws.length; ++i) {
teichmann@4259:                 cf.addObservedPoint(qs[i], ws[i]);
teichmann@4259:             }
teichmann@4259: 
teichmann@4259:             try {
teichmann@4259:                 coeffs = cf.fit(function, function.getInitialGuess());
teichmann@4259:                 break;
teichmann@4259:             }
teichmann@4259:             catch (MathException me) {
teichmann@4259:                 if (log.isDebugEnabled()) {
teichmann@4259:                     log.debug("tolerance " + tolerance + " + failed.");
teichmann@4259:                 }
teichmann@4259:             }
teichmann@4259:         }
teichmann@4385:         if (coeffs != null) {
teichmann@4385:             chiSqr[0] = lmo.getChiSquare();
teichmann@4259:         }
teichmann@4259:         return coeffs;
teichmann@4259:     }
teichmann@4259: 
teichmann@4355:     protected double [][] extractPointsToFit(double [][] wqs) {
teichmann@4259: 
teichmann@4259:         double [] ws = wqs[0];
teichmann@4259:         double [] qs = wqs[1];
teichmann@4259: 
teichmann@4259:         int N = Math.min(ws.length, qs.length);
teichmann@4259: 
teichmann@4259:         if (N < 2) {
teichmann@4259:             log.warn("Too less points for fitting");
teichmann@4259:             return null;
teichmann@4259:         }
teichmann@4259: 
teichmann@4259:         double q2 = qs[N-1];
teichmann@4259:         double w2 = ws[N-1];
teichmann@4259:         double q1 = qs[N-2];
teichmann@4259:         double w1 = ws[N-2];
teichmann@4259: 
teichmann@4385:         boolean ascending = w2 > w1;
teichmann@4385: 
teichmann@4385:         TDoubleArrayList ows = new TDoubleArrayList();
teichmann@4385:         TDoubleArrayList oqs = new TDoubleArrayList();
teichmann@4385: 
teichmann@4259:         oqs.add(q2); oqs.add(q1);
teichmann@4259:         ows.add(w2); ows.add(w1);
teichmann@4259: 
teichmann@4385:         int lastDir = -2;
teichmann@4259: 
teichmann@4259:         for (int i = N-3; i >= 0; --i) {
teichmann@4385:             double q = qs[i];
teichmann@4385:             double w = ws[i];
teichmann@4259: 
teichmann@4385:             if ((ascending && w > w1) || (!ascending && w < w1)) {
teichmann@4259:                 break;
teichmann@4259:             }
teichmann@4259: 
teichmann@4385:             int dir = Line2D.relativeCCW(q2, w2, q1, w1, q, w);
teichmann@4385:             //int dir = Utils.relativeCCW(q2, w2, q1, w1, q, w);
teichmann@4385:             if (lastDir == -2) {
teichmann@4385:                 lastDir = dir;
teichmann@4259:             }
teichmann@4385:             else if (lastDir != dir) {
teichmann@4385:                 break;
teichmann@4385:             }
teichmann@4385: 
teichmann@4385:             oqs.add(q);
teichmann@4385:             ows.add(w);
teichmann@4385:             w2 = w1;
teichmann@4385:             q2 = q1;
teichmann@4385:             w1 = w;
teichmann@4385:             q1 = q;
teichmann@4259:         }
teichmann@4259: 
teichmann@4355:         oqs.reverse();
teichmann@4355:         ows.reverse();
teichmann@4385: 
teichmann@4385:         boolean debug = log.isDebugEnabled();
teichmann@4385:         if (debug) {
teichmann@4385:             log.debug("from table: " + N);
teichmann@4385:             log.debug("after trim: " + oqs.size());
teichmann@4385:         }
teichmann@4385: 
teichmann@4355:         cutPercent(ows, oqs);
teichmann@4259: 
teichmann@4385:         if (debug) {
teichmann@4385:             log.debug("after percent cut: " + oqs.size());
teichmann@4385:         }
teichmann@4385: 
teichmann@4259:         return new double [][] {
teichmann@4259:             ows.toNativeArray(),
teichmann@4259:             oqs.toNativeArray()
teichmann@4259:         };
teichmann@4259:     }
teichmann@4355: 
teichmann@4385: 
teichmann@4355:     protected void cutPercent(TDoubleArrayList ws, TDoubleArrayList qs) {
teichmann@4355:         int N = qs.size();
teichmann@4355:         if (percent <= 0d || N == 0) {
teichmann@4355:             return;
teichmann@4355:         }
teichmann@4355: 
teichmann@4355:         double minQ = qs.getQuick(0);
teichmann@4355:         double maxQ = qs.getQuick(N-1);
teichmann@4355:         double factor = Math.min(Math.max(0d, percent/100d), 1d);
teichmann@4355:         double cutQ = Linear.weight(factor, minQ, maxQ);
teichmann@4355:         int cutIndex = 0;
teichmann@4355:         for (; cutIndex < N; ++cutIndex) {
teichmann@4355:             double q = qs.getQuick(cutIndex);
teichmann@4355:             if (minQ < maxQ) {
teichmann@4355:                 if (q > cutQ) {
teichmann@4355:                     break;
teichmann@4355:                 }
teichmann@4355:             }
teichmann@4355:             else {
teichmann@4355:                 if (q < cutQ) {
teichmann@4355:                     break;
teichmann@4355:                 }
teichmann@4355:             }
teichmann@4355:         }
teichmann@4355:         ws.remove(0, cutIndex);
teichmann@4355:         qs.remove(0, cutIndex);
teichmann@4355:     }
sascha@3775: }
sascha@3775: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf-8 :