Mercurial > dive4elements > river
view flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/BackJumpCorrector.java @ 999:b0218f21c664
Wire datacage with backend of artifact database.
flys-artifacts/trunk@2438 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Mon, 01 Aug 2011 14:58:15 +0000 |
parents | c09c9e05ecfa |
children | 0b6dac664bbb |
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; 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 = isIncreasing(ws); if (wsUp) { km = swapClone(km); ws = swapClone(ws); } boolean kmUp = isIncreasing(km); if (!kmUp) { km = sumDiffs(km); } if (log.isDebugEnabled()) { log.debug("BackJumpCorrector.doCorrection ------- enter"); log.debug(" km increasing: " + isIncreasing(km)); log.debug(" ws increasing: " + isIncreasing(ws)); log.debug("BackJumpCorrector.doCorrection ------- leave"); } boolean hasBackJumps = doCorrectionClean(km, ws, errors); if (hasBackJumps && wsUp) { // mirror back 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(); } int back = i-2; double distance = Math.abs(km[i] - km[i-1]); double rest = 0.0; double ikm = 0.0; while (back > -1) { if (ws[back] < ws[i]) { // under water // continue scanning backwards distance += Math.abs(km[back+1] - km[back]); --back; continue; } if (ws[back] > ws[i]) { // over water // need to calculate intersection log.debug("intersection case"); double m = (km[back+1]-km[back])/(ws[back+1]-ws[back]); double b = km[back]-ws[back]*m; ikm = m*ws[i] + b; distance += Math.abs(ikm - km[back+1]); rest = Math.abs(km[back] - ikm); } else { // equals height at ws[back] log.debug("exact height match"); distance += Math.abs(km[back+1] - km[back]); ikm = km[back]; } break; } if (back <= 0) { //log.debug("run over left border"); // fill with same as ws[i] for (int j = 0; j < i; ++j) { ws[j] = ws[i]; } continue; } double quarterDistance = 0.25*distance; // Now further seek back for the max height double restDistance = Math.max(0.0, quarterDistance - rest); --back; double mkmw = ws[i]; double mkm = km[0]; while (back > -1) { double d = Math.abs(km[back+1] - km[back]); restDistance -= d; if (restDistance > 0.0) { --back; continue; } if (restDistance < 0.0) { // need to calculate intersection if (km[back] == km[back+1]) { // should not happen mkm = km[back]; mkmw = 0.5*(ws[back] + ws[back+1]); } else { double m = (ws[back+1]-ws[back])/(km[back+1]-km[back]); double b = ws[back] - km[back]*m; mkm = km[back] + restDistance; mkmw = m*mkm + b; } } else { // exact match mkm = km[back]; mkmw = ws[back]; } break; } double factor = back >= 0 && Math.abs(restDistance) < 1e-4 ? 1.0 : 1.0 - Math.min(1, Math.max(0, restDistance/quarterDistance)); double ikmw = factor*0.25*(mkmw-ws[i]) + ws[i]; double end = ikm + quarterDistance*factor; double [] x = { mkm, ikm, end }; double [] y = { mkmw, ikmw, ws[i] }; if (interpolator == null) { interpolator = new SplineInterpolator(); } if (log.isDebugEnabled()) { for (int j = 0; j < x.length; ++j) { log.debug(" " + x[j] + " -> " + y[j]); } } PolynomialSplineFunction spline; try { spline = interpolator.interpolate(x, y); } catch (MathIllegalArgumentException miae) { // TODO: I18N errors.addProblem("creating spline interpolation 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])); } } for (back = Math.max(back, 0); back < i && km[back] < end; ++back ) { // to 3/4 point fill with spline values ws[back] = spline.value(km[back]); } while (back < i) { // fill the rest with ws[i] ws[back++] = ws[i]; } } catch (ArgumentOutsideDomainException aode) { // TODO: I18N errors.addProblem("spline interpolation failed."); log.error("spline interpolation failed", aode); } } // for all km return !backjumps.isEmpty(); } public static final boolean isIncreasing(double [] array) { int inc = 0; int dec = 0; int sweet = (array.length-1)/2; for (int i = 1; i < array.length; ++i) { if (array[i] > array[i-1]) { if (++inc > sweet) { return true; } } else if (++dec > sweet) { return false; } } return inc > sweet; } public static final double [] swap(double [] array) { int lo = 0; int hi = array.length-1; while (hi > lo) { double t = array[lo]; array[lo] = array[hi]; array[hi] = t; ++lo; --hi; } return array; } public static final double [] swapClone(double [] in) { double [] out = new double[in.length]; for (int j = out.length-1, i = 0; j >= 0;) { out[j--] = in[i++]; } return out; } public static final double [] sumDiffs(double [] in) { double [] out = new double[in.length]; for (int i = 1; i < out.length; ++i) { out[i] = out[i-1] + Math.abs(in[i-1] - in[i]); } return out; } } // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :