Mercurial > dive4elements > river
view flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/ExtremeCalculation.java @ 4655:cd44d28d0fbc
Move the access to artifact data to the Access object
Use BedHeightAccess class to receive the data from the artifact. This abstracts
the data access from the actual artifact.
author | Björn Ricks <bjoern.ricks@intevation.de> |
---|---|
date | Tue, 11 Dec 2012 09:44:04 +0100 |
parents | 5b8919ef601d |
children |
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; 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 :