Mercurial > dive4elements > river
view flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/ExtremeCalculation.java @ 4357:11cf29eb002c
FacetTypes: New type: exreme_wq_curve_base (for the Stuetz-Q-Poits).
author | Felix Wolfsteller <felix.wolfsteller@intevation.de> |
---|---|
date | Fri, 02 Nov 2012 11:58:11 +0100 |
parents | d095b4267772 |
children | 10e277c2fe0f |
line wrap: on
line source
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.Linear; import de.intevation.flys.artifacts.math.Utils; import de.intevation.flys.artifacts.math.fitting.Function; import de.intevation.flys.artifacts.math.fitting.FunctionFactory; import de.intevation.flys.artifacts.model.Calculation; 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; /** 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); } 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(); for (double km = from; km <= to; km += step) { double currentKm = DoubleUtil.round(km); if (range == null || !range.inside(currentKm)) { for (RangeWithValues r: ranges) { if (r.inside(currentKm)) { range = r; break; } } // TODO: i18n addProblem(currentKm, "extreme.no.range.inner"); continue; } double [][] wqs = wst.interpolateTabulated(currentKm); if (wqs == null) { // TODO: i18n addProblem(currentKm, "extreme.no.raw.data"); continue; } // XXX: This should not be necessary for model data. if (!DoubleUtil.isValid(wqs)) { // TODO: i18n addProblem(currentKm, "extreme.invalid.data"); continue; } 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, currentKm); } } } 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 = 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 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; } oqs.reverse(); ows.reverse(); cutPercent(ows, oqs); 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 :