Mercurial > dive4elements > river
view flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/BackJumpCorrector.java @ 4922:c30771dbf1bd
CollectionAttribute: Comment.
author | Felix Wolfsteller <felix.wolfsteller@intevation.de> |
---|---|
date | Thu, 31 Jan 2013 09:28:18 +0100 |
parents | 5642a83420f2 |
children |
line wrap: on
line source
package de.intevation.flys.artifacts.math; import java.util.ArrayList; import java.util.List; import java.io.Serializable; import org.apache.commons.math.analysis.interpolation.SplineInterpolator; import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction; import org.apache.commons.math.ArgumentOutsideDomainException; import org.apache.commons.math.exception.MathIllegalArgumentException; import org.apache.log4j.Logger; import de.intevation.flys.artifacts.model.Calculation; import de.intevation.flys.utils.DoubleUtil; public class BackJumpCorrector implements Serializable { private static Logger log = Logger.getLogger(BackJumpCorrector.class); protected ArrayList<Double> backjumps; protected double [] corrected; public BackJumpCorrector() { backjumps = new ArrayList<Double>(); } public boolean hasBackJumps() { return !backjumps.isEmpty(); } public List<Double> getBackJumps() { return backjumps; } public double [] getCorrected() { return corrected; } public boolean doCorrection( double [] km, double [] ws, Calculation errors ) { boolean wsUp = DoubleUtil.isIncreasing(ws); if (wsUp) { km = DoubleUtil.swapClone(km); ws = DoubleUtil.swapClone(ws); } boolean kmUp = DoubleUtil.isIncreasing(km); if (!kmUp) { km = DoubleUtil.sumDiffs(km); } if (log.isDebugEnabled()) { log.debug("BackJumpCorrector.doCorrection ------- enter"); log.debug(" km increasing: " + DoubleUtil.isIncreasing(km)); log.debug(" ws increasing: " + DoubleUtil.isIncreasing(ws)); log.debug("BackJumpCorrector.doCorrection ------- leave"); } boolean hasBackJumps = doCorrectionClean(km, ws, errors); if (hasBackJumps && wsUp) { // mirror back DoubleUtil.swap(corrected); } return hasBackJumps; } protected boolean doCorrectionClean( double [] km, double [] ws, Calculation errors ) { int N = km.length; if (N != ws.length) { throw new IllegalArgumentException("km.length != ws.length"); } if (N < 2) { return false; } SplineInterpolator interpolator = null; for (int i = 1; i < N; ++i) { if (ws[i] <= ws[i-1]) { // no back jump continue; } backjumps.add(km[i]); if (corrected == null) { // lazy cloning ws = corrected = (double [])ws.clone(); } double above = aboveWaterKM(km, ws, i); if (Double.isNaN(above)) { // run over start km // fill all previous for (int j = 0; j < i; ++j) { ws[j] = ws[i]; } continue; } double distance = Math.abs(km[i] - above); double quarterDistance = 0.25*distance; double start = above - quarterDistance; double startHeight = DoubleUtil.interpolateSorted(km, ws, start); if (Double.isNaN(startHeight)) { // run over start km startHeight = ws[0]; } double between = above + quarterDistance; double aboveHeight = ws[i] + 0.25*(startHeight - ws[i]); double [] x = { start, above, between }; double [] y = { startHeight, aboveHeight, ws[i] }; if (log.isDebugEnabled()) { for (int j = 0; j < x.length; ++j) { log.debug(" " + x[j] + " -> " + y[j]); } } if (interpolator == null) { interpolator = new SplineInterpolator(); } PolynomialSplineFunction spline; try { spline = interpolator.interpolate(x, y); } catch (MathIllegalArgumentException miae) { errors.addProblem("spline.creation.failed"); log.error(miae); continue; } try { if (log.isDebugEnabled()) { log.debug("spline points:"); for (int j = 0; j < x.length; ++j) { log.debug(x[j] + " " + y[j] + " " + spline.value(x[j])); } } int j = i-1; for (; j >= 0 && km[j] >= between; --j) { ws[j] = ws[i]; } for (; j >= 0 && ws[j] < startHeight; --j) { ws[j] = spline.value(km[j]); } } catch (ArgumentOutsideDomainException aode) { errors.addProblem("spline.interpolation.failed"); log.error("spline interpolation failed", aode); } } // for all km return !backjumps.isEmpty(); } protected static double aboveWaterKM( double [] km, double [] ws, int wIndex ) { double w = ws[wIndex]; while (--wIndex >= 0) { // still under water if (ws[wIndex] < w) continue; if (ws[wIndex] > w) { // f(ws[wIndex]) = km[wIndex] // f(ws[wIndex+1]) = km[wIndex+1] return Linear.linear( w, ws[wIndex], ws[wIndex+1], km[wIndex], km[wIndex+1]); } else { return km[wIndex]; } } return Double.NaN; } } // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :