Mercurial > dive4elements > river
changeset 4259:5cc9453456a7
First complete but untested version of the 'Auslagerung extremer Wasserspiegellagen' calculation.
author | Sascha L. Teichmann <teichmann@intevation.de> |
---|---|
date | Thu, 25 Oct 2012 17:25:37 +0200 |
parents | 2c6e571f366a |
children | e4a415773b0a 270f3ac8a82e |
files | flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/Utils.java flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/Curve.java flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/ExtremeCalculation.java flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/ExtremeResult.java |
diffstat | 4 files changed, 348 insertions(+), 7 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/Utils.java Thu Oct 25 17:25:37 2012 +0200 @@ -0,0 +1,61 @@ +package de.intevation.flys.artifacts.math; + + +public final class Utils { + + public static final double EPSILON = 1e-3; + + private Utils() { + } + + public static final boolean epsilonEquals(double a, double b) { + return epsilonEquals(a, b, EPSILON); + } + + public static final boolean epsilonEquals(double a, double b, double eps) { + return Math.abs(a - b) < eps; + } + + public static int relativeCCW( + double x1, double y1, + double x2, double y2, + double px, double py + ) { + if ((epsilonEquals(x1, x2) && epsilonEquals(y1, y2)) + || ((epsilonEquals(x1, px) && epsilonEquals(y1, py)))) { + return 0; // Coincident points. + } + // Translate to the origin. + x2 -= x1; + y2 -= y1; + px -= x1; + py -= y1; + double slope2 = y2 / x2; + double slopep = py / px; + if (epsilonEquals(slope2, slopep) + || (epsilonEquals(x2, 0.0) && epsilonEquals(px, 0.0))) { + return y2 > EPSILON // Colinear. + ? (py < -EPSILON ? -1 : py > y2 ? 1 : 0) + : (py > -EPSILON ? -1 : py < y2 ? 1 : 0); + } + if (x2 >= EPSILON && slope2 >= EPSILON) { + return px >= EPSILON // Quadrant 1. + ? (slope2 > slopep ? 1 : -1) + : (slope2 < slopep ? 1 : -1); + } + + if (y2 > EPSILON) { + return px < -EPSILON // Quadrant 2. + ? (slope2 > slopep ? 1 : -1) + : (slope2 < slopep ? 1 : -1); + } + if (slope2 >= EPSILON) { + return px >= EPSILON // Quadrant 3. + ? (slope2 < slopep ? 1 : -1) + : (slope2 > slopep ? 1 : -1); + } + return px < -EPSILON // Quadrant 4. + ? (slope2 < slopep ? 1 : -1) + : (slope2 > slopep ? 1 : -1); + } +}
--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/Curve.java Thu Oct 25 15:17:01 2012 +0200 +++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/Curve.java Thu Oct 25 17:25:37 2012 +0200 @@ -25,6 +25,7 @@ protected double [] ws; protected String function; protected double [] coeffs; + protected double chiSqr; // The spline is pretty heavy weight so cache it with a soft ref only. protected transient SoftReference<Function> spline; @@ -37,7 +38,8 @@ double [] qs, double [] ws, String function, - double [] coeffs + double [] coeffs, + double chiSqr ) { this.qs = qs; this.ws = ws; @@ -81,6 +83,24 @@ return extrapolation; } + /** + * Gets the chiSqr for this instance. + * + * @return The chiSqr. + */ + public double getChiSqr() { + return this.chiSqr; + } + + /** + * Sets the chiSqr for this instance. + * + * @param chiSqr The chiSqr. + */ + public void setChiSqr(double chiSqr) { + this.chiSqr = chiSqr; + } + protected synchronized Function getSpline() { Function sp; if (spline != null) {
--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/ExtremeCalculation.java Thu Oct 25 15:17:01 2012 +0200 +++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/ExtremeCalculation.java Thu Oct 25 17:25:37 2012 +0200 @@ -1,7 +1,18 @@ 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; @@ -9,12 +20,16 @@ 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; @@ -22,6 +37,9 @@ public class ExtremeCalculation extends Calculation { + private static final Log log = + LogFactory.getLog(ExtremeCalculation.class); + protected String river; protected String function; protected double from; @@ -119,6 +137,40 @@ : 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( @@ -127,7 +179,12 @@ ) { RangeWithValues range = null; - KMs: for (double km = from; km <= to; km += step) { + 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)) { @@ -139,7 +196,7 @@ } // TODO: i18n addProblem(currentKm, "extreme.no.range"); - continue KMs; + continue; } double [][] wqs = wst.interpolateTabulated(currentKm); @@ -155,13 +212,148 @@ addProblem(currentKm, "extreme.invalid.data"); continue; } - // TODO: Implement extraction of points for curve fitting. - // TODO: Implement curve fitting. - // TODO: Implement generating Curve object per km. + + 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(); + 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 :
--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/ExtremeResult.java Thu Oct 25 15:17:01 2012 +0200 +++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/ExtremeResult.java Thu Oct 25 17:25:37 2012 +0200 @@ -2,10 +2,78 @@ import java.io.Serializable; +import de.intevation.flys.artifacts.model.WQKms; + +import de.intevation.flys.utils.KMIndex; + public class ExtremeResult implements Serializable { + protected KMIndex<Curve> curves; + protected WQKms [] wqkms; + public ExtremeResult() { } + + public ExtremeResult(KMIndex<Curve> curves, WQKms [] wqkms) { + this.curves = curves; + this.wqkms = wqkms; + } + + /** + * Gets the curves for this instance. + * + * @return The curves. + */ + public KMIndex<Curve> getCurves() { + return this.curves; + } + + /** + * Sets the curves for this instance. + * + * @param curves The curves. + */ + public void setCurves(KMIndex<Curve> curves) { + this.curves = curves; + } + + /** + * Gets the wqkms for this instance. + * + * @return The wqkms. + */ + public WQKms[] getWQKms() { + return this.wqkms; + } + + /** + * Gets the wqkms for this instance. + * + * @param index The index to get. + * @return The wqkms. + */ + public WQKms getWQKms(int index) { + return this.wqkms[index]; + } + + /** + * Sets the wqkms for this instance. + * + * @param wqkms The wqkms. + */ + public void setWQKms(WQKms[] wqkms) { + this.wqkms = wqkms; + } + + /** + * Sets the wqkms for this instance. + * + * @param index The index to set. + * @param wqkms The wqkms. + */ + public void setWQKms(int index, WQKms wqkms) { + this.wqkms[index] = wqkms; + } } // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :