Mercurial > dive4elements > river
view flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/ExtremeCalculation.java @ 4282:8b4988815974
Added marker for Ws and Qs in Historical Discharge WQ charts.
Therefore, the XYChartGenerator got two new methods addDomainMarker(Marker, boolean) and addValueMarker(Marker, boolean).
The boolean parameters determine, if the marker should be visible or not. This is analogous to addAxisSeries(XYSeries, int, boolean).
author | Ingo Weinzierl <ingo.weinzierl@intevation.de> |
---|---|
date | Mon, 29 Oct 2012 05:59:27 +0100 |
parents | 5cc9453456a7 |
children | 10c1a413a1e0 |
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.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"); 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); } } } ExtremeResult result = new ExtremeResult(curves, wqkms); return new CalculationResult(result, this); } protected static 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 static 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; } // XXX: Not really needed for fitting. // oqs.reverse(); // ows.reverse(); return new double [][] { ows.toNativeArray(), oqs.toNativeArray() }; } } // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf-8 :