[Dive4elements-commits] [PATCH] First complete but untested version of the 'Auslagerung extremer Wasserspiegellagen' calculation
Wald Commits
scm-commit at wald.intevation.org
Thu Oct 25 17:25:47 CEST 2012
# HG changeset patch
# User Sascha L. Teichmann <teichmann at intevation.de>
# Date 1351178737 -7200
# Node ID 5cc9453456a7ee7bdfffab3d7aa004e858ee6db7
# Parent 2c6e571f366ae3cd488549cd67b1aacab3f79269
First complete but untested version of the 'Auslagerung extremer Wasserspiegellagen' calculation.
diff -r 2c6e571f366a -r 5cc9453456a7 flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/Utils.java
--- /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);
+ }
+}
diff -r 2c6e571f366a -r 5cc9453456a7 flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/Curve.java
--- 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) {
diff -r 2c6e571f366a -r 5cc9453456a7 flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/ExtremeCalculation.java
--- 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 :
diff -r 2c6e571f366a -r 5cc9453456a7 flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/extreme/ExtremeResult.java
--- 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 :
More information about the Dive4elements-commits
mailing list