Mercurial > dive4elements > river
diff artifacts/src/main/java/org/dive4elements/river/artifacts/model/extreme/ExtremeCalculation.java @ 5838:5aa05a7a34b7
Rename modules to more fitting names.
author | Sascha L. Teichmann <teichmann@intevation.de> |
---|---|
date | Thu, 25 Apr 2013 15:23:37 +0200 |
parents | flys-artifacts/src/main/java/org/dive4elements/river/artifacts/model/extreme/ExtremeCalculation.java@bd047b71ab37 |
children | 4897a58c8746 |
line wrap: on
line diff
--- /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 :