teichmann@5863: /* Copyright (C) 2011, 2012, 2013 by Bundesanstalt für Gewässerkunde
teichmann@5863:  * Software engineering by Intevation GmbH
teichmann@5863:  *
teichmann@5994:  * This file is Free Software under the GNU AGPL (>=v3)
teichmann@5863:  * and comes with ABSOLUTELY NO WARRANTY! Check out the
teichmann@5994:  * documentation coming with Dive4Elements River for details.
teichmann@5863:  */
teichmann@5863: 
teichmann@5831: package org.dive4elements.river.artifacts.math;
sascha@343: 
sascha@343: import java.util.ArrayList;
sascha@343: import java.util.List;
sascha@343: 
sascha@343: import java.io.Serializable;
sascha@343: 
sascha@343: import org.apache.commons.math.analysis.interpolation.SplineInterpolator;
sascha@343: 
sascha@343: import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
sascha@343: 
sascha@343: import org.apache.commons.math.ArgumentOutsideDomainException;
sascha@343: 
sascha@649: import org.apache.commons.math.exception.MathIllegalArgumentException;
sascha@649: 
sascha@343: import org.apache.log4j.Logger;
sascha@343: 
teichmann@5831: import org.dive4elements.river.artifacts.model.Calculation;
sascha@674: 
teichmann@5831: import org.dive4elements.river.utils.DoubleUtil;
sascha@1672: 
sascha@343: public class BackJumpCorrector
sascha@343: implements   Serializable
sascha@343: {
sascha@343:     private static Logger log = Logger.getLogger(BackJumpCorrector.class);
sascha@343: 
sascha@343:     protected ArrayList<Double> backjumps;
sascha@343: 
sascha@343:     protected double [] corrected;
sascha@343: 
sascha@343:     public BackJumpCorrector() {
sascha@343:         backjumps = new ArrayList<Double>();
sascha@343:     }
sascha@343: 
sascha@343:     public boolean hasBackJumps() {
sascha@343:         return !backjumps.isEmpty();
sascha@343:     }
sascha@343: 
sascha@343:     public List<Double> getBackJumps() {
sascha@343:         return backjumps;
sascha@343:     }
sascha@343: 
sascha@343:     public double [] getCorrected() {
sascha@343:         return corrected;
sascha@343:     }
sascha@343: 
sascha@674:     public boolean doCorrection(
sascha@674:         double []   km,
sascha@674:         double []   ws,
sascha@674:         Calculation errors
sascha@674:     ) {
sascha@1672:         boolean wsUp = DoubleUtil.isIncreasing(ws);
sascha@642: 
sascha@642:         if (wsUp) {
sascha@1672:             km = DoubleUtil.swapClone(km);
sascha@1672:             ws = DoubleUtil.swapClone(ws);
sascha@406:         }
sascha@406: 
sascha@1672:         boolean kmUp = DoubleUtil.isIncreasing(km);
sascha@642: 
sascha@642:         if (!kmUp) {
sascha@1672:             km = DoubleUtil.sumDiffs(km);
sascha@642:         }
sascha@642: 
sascha@642:         if (log.isDebugEnabled()) {
sascha@642:             log.debug("BackJumpCorrector.doCorrection ------- enter");
sascha@1672:             log.debug("  km increasing: " + DoubleUtil.isIncreasing(km));
sascha@1672:             log.debug("  ws increasing: " + DoubleUtil.isIncreasing(ws));
sascha@642:             log.debug("BackJumpCorrector.doCorrection ------- leave");
sascha@406:         }
sascha@406: 
sascha@674:         boolean hasBackJumps = doCorrectionClean(km, ws, errors);
sascha@406: 
sascha@642:         if (hasBackJumps && wsUp) {
sascha@406:             // mirror back
sascha@1672:             DoubleUtil.swap(corrected);
sascha@406:         }
sascha@406: 
sascha@406:         return hasBackJumps;
sascha@406:     }
sascha@406: 
sascha@674:     protected boolean doCorrectionClean(
sascha@674:         double []   km,
sascha@674:         double []   ws,
sascha@674:         Calculation errors
sascha@674:     ) {
sascha@343:         int N = km.length;
sascha@343: 
sascha@343:         if (N != ws.length) {
sascha@343:             throw new IllegalArgumentException("km.length != ws.length");
sascha@343:         }
sascha@343: 
sascha@343:         if (N < 2) {
sascha@343:             return false;
sascha@343:         }
sascha@343: 
sascha@343:         SplineInterpolator interpolator = null;
sascha@343: 
sascha@343:         for (int i = 1; i < N; ++i) {
sascha@343:             if (ws[i] <= ws[i-1]) {
sascha@343:                 // no back jump
sascha@343:                 continue;
sascha@343:             }
sascha@343:             backjumps.add(km[i]);
sascha@343: 
sascha@343:             if (corrected == null) {
sascha@343:                 // lazy cloning
sascha@343:                 ws = corrected = (double [])ws.clone();
sascha@343:             }
sascha@343: 
sascha@1674:             double above = aboveWaterKM(km, ws, i);
sascha@343: 
sascha@1674:             if (Double.isNaN(above)) { // run over start km
sascha@1674:                 // fill all previous
sascha@343:                 for (int j = 0; j < i; ++j) {
sascha@343:                     ws[j] = ws[i];
sascha@343:                 }
sascha@343:                 continue;
sascha@343:             }
sascha@343: 
sascha@1674:             double distance = Math.abs(km[i] - above);
sascha@1674: 
sascha@343:             double quarterDistance = 0.25*distance;
sascha@343: 
sascha@1674:             double start = above - quarterDistance;
sascha@343: 
sascha@1674:             double startHeight = DoubleUtil.interpolateSorted(km, ws, start);
sascha@1674: 
sascha@1674:             if (Double.isNaN(startHeight)) {
sascha@1674:                 // run over start km
sascha@1674:                 startHeight = ws[0];
sascha@343:             }
sascha@343: 
sascha@1674:             double between = above + quarterDistance;
sascha@343: 
sascha@1676:             double aboveHeight = ws[i] + 0.25*(startHeight - ws[i]);
sascha@343: 
sascha@1676:             double [] x = { start,  above,  between };
sascha@1676:             double [] y = { startHeight, aboveHeight, ws[i] };
sascha@343: 
sascha@642:             if (log.isDebugEnabled()) {
sascha@642:                 for (int j = 0; j < x.length; ++j) {
sascha@642:                     log.debug("   " + x[j] + " -> " + y[j]);
sascha@642:                 }
sascha@642:             }
sascha@642: 
sascha@1674:             if (interpolator == null) {
sascha@1674:                 interpolator = new SplineInterpolator();
sascha@1674:             }
sascha@1674: 
sascha@649:             PolynomialSplineFunction spline;
sascha@742: 
sascha@649:             try {
sascha@649:                 spline = interpolator.interpolate(x, y);
sascha@649:             }
sascha@649:             catch (MathIllegalArgumentException miae) {
sascha@2166:                 errors.addProblem("spline.creation.failed");
sascha@649:                 log.error(miae);
sascha@649:                 continue;
sascha@649:             }
sascha@343: 
sascha@343:             try {
sascha@343:                 if (log.isDebugEnabled()) {
sascha@343:                     log.debug("spline points:");
sascha@343:                     for (int j = 0; j < x.length; ++j) {
sascha@343:                         log.debug(x[j] + " " + y[j] + " " + spline.value(x[j]));
sascha@343:                     }
sascha@343:                 }
sascha@343: 
sascha@1676:                 int j = i-1;
sascha@1676: 
sascha@1676:                 for (; j >= 0 && km[j] >= between; --j) {
sascha@1676:                     ws[j] = ws[i];
sascha@1676:                 }
sascha@1676: 
sascha@1676:                 for (; j >= 0 && ws[j] < startHeight; --j) {
sascha@1674:                     ws[j] = spline.value(km[j]);
sascha@343:                 }
sascha@343:             }
sascha@343:             catch (ArgumentOutsideDomainException aode) {
sascha@2166:                 errors.addProblem("spline.interpolation.failed");
sascha@343:                 log.error("spline interpolation failed", aode);
sascha@343:             }
sascha@343:         } // for all km
sascha@343: 
sascha@343:         return !backjumps.isEmpty();
sascha@343:     }
sascha@1674: 
sascha@1674: 
sascha@1674:     protected static double aboveWaterKM(
sascha@3076:         double [] km,
sascha@3076:         double [] ws,
sascha@1674:         int       wIndex
sascha@1674:     ) {
sascha@1674:         double w = ws[wIndex];
sascha@1674: 
sascha@1674:         while (--wIndex >= 0) {
sascha@1674:             // still under water
sascha@1674:             if (ws[wIndex] < w) continue;
sascha@1674: 
sascha@1674:             if (ws[wIndex] > w) {
sascha@1674:                 // f(ws[wIndex])   = km[wIndex]
sascha@1674:                 // f(ws[wIndex+1]) = km[wIndex+1]
sascha@1674:                 return Linear.linear(
sascha@1674:                     w,
sascha@1674:                     ws[wIndex], ws[wIndex+1],
sascha@1674:                     km[wIndex], km[wIndex+1]);
sascha@1674:             }
sascha@1674:             else {
sascha@1674:                 return km[wIndex];
sascha@1674:             }
sascha@1674:         }
sascha@1674: 
sascha@1674:         return Double.NaN;
sascha@1674:     }
sascha@343: }
sascha@343: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :