sascha@343: package de.intevation.flys.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@343: import org.apache.log4j.Logger;
sascha@343: 
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@642:     public boolean doCorrection(double [] km, double [] ws) {
sascha@642: 
sascha@642:         boolean wsUp = isIncreasing(ws);
sascha@642: 
sascha@642:         if (wsUp) {
sascha@642:             km = swapClone(km);
sascha@642:             ws = swapClone(ws);
sascha@406:         }
sascha@406: 
sascha@642:         boolean kmUp = isIncreasing(km);
sascha@642: 
sascha@642:         if (!kmUp) {
sascha@642:             km = sumDiffs(km);
sascha@642:         }
sascha@642: 
sascha@642:         if (log.isDebugEnabled()) {
sascha@642:             log.debug("BackJumpCorrector.doCorrection ------- enter");
sascha@642:             log.debug("  km increasing: " + isIncreasing(km));
sascha@642:             log.debug("  ws increasing: " + isIncreasing(ws));
sascha@642:             log.debug("BackJumpCorrector.doCorrection ------- leave");
sascha@406:         }
sascha@406: 
sascha@406:         boolean hasBackJumps = doCorrectionClean(km, ws);
sascha@406: 
sascha@642:         if (hasBackJumps && wsUp) {
sascha@406:             // mirror back
sascha@642:             swap(corrected);
sascha@406:         }
sascha@406: 
sascha@406:         return hasBackJumps;
sascha@406:     }
sascha@406: 
sascha@406:     protected boolean doCorrectionClean(double [] km, double [] ws) {
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@343:             int    back = i-2;
sascha@343:             double distance = Math.abs(km[i] - km[i-1]);
sascha@343:             double rest = 0.0;
sascha@343:             double ikm  = 0.0;
sascha@343: 
sascha@343:             while (back > -1) {
sascha@343:                 if (ws[back] < ws[i]) { // under water
sascha@343:                     // continue scanning backwards
sascha@343:                     distance += Math.abs(km[back+1] - km[back]);
sascha@343:                     --back;
sascha@343:                     continue;
sascha@343:                 }
sascha@343:                 if (ws[back] > ws[i]) { // over water
sascha@343:                     // need to calculate intersection
sascha@343:                     log.debug("intersection case");
sascha@343:                     double m = (km[back+1]-km[back])/(ws[back+1]-ws[back]);
sascha@343:                     double b = km[back]-ws[back]*m;
sascha@343:                     ikm = m*ws[i] + b;
sascha@343:                     distance += Math.abs(ikm - km[back+1]);
sascha@343:                     rest = Math.abs(km[back] - ikm);
sascha@343:                 }
sascha@343:                 else {
sascha@343:                     // equals height at ws[back]
sascha@343:                     log.debug("exact height match");
sascha@343:                     distance += Math.abs(km[back+1] - km[back]);
sascha@343:                     ikm = km[back];
sascha@343:                 }
sascha@343:                 break;
sascha@343:             }
sascha@343: 
sascha@343:             if (back < 0) {
sascha@642:                 //log.debug("run over left border");
sascha@343:                 // fill with same as ws[i]
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@343:             double quarterDistance = 0.25*distance;
sascha@343: 
sascha@343:             // Now further seek back for the max height
sascha@343: 
sascha@343:             double restDistance = Math.max(0.0, quarterDistance - rest);
sascha@343:             --back;
sascha@343: 
sascha@343:             double mkmw = ws[i];
sascha@343:             double mkm  = km[0];
sascha@343: 
sascha@343:             while (back > -1) {
sascha@343:                 double d = Math.abs(km[back+1] - km[back]);
sascha@343:                 restDistance -= d;
sascha@343:                 if (restDistance > 0.0) {
sascha@343:                     --back;
sascha@343:                     continue;
sascha@343:                 }
sascha@343:                 if (restDistance < 0.0) {
sascha@343:                     // need to calculate intersection
sascha@343:                     if (km[back] == km[back+1]) { // should not happen
sascha@343:                         mkm = km[back];
sascha@343:                         mkmw = 0.5*(ws[back] + ws[back+1]);
sascha@343:                     }
sascha@343:                     else {
sascha@343:                         double m = (ws[back+1]-ws[back])/(km[back+1]-km[back]);
sascha@343:                         double b = ws[back] - km[back]*m;
sascha@343:                         mkm = km[back] + restDistance;
sascha@343:                         mkmw = m*mkm + b;
sascha@343:                     }
sascha@343:                 }
sascha@343:                 else {
sascha@343:                     // exact match
sascha@343:                     mkm  = km[back];
sascha@343:                     mkmw = ws[back];
sascha@343:                 }
sascha@343:                 break;
sascha@343:             }
sascha@343: 
sascha@343:             double factor = back >= 0 && Math.abs(restDistance) < 1e-4
sascha@343:                 ? 1.0
sascha@343:                 : 1.0 - Math.min(1, Math.max(0, restDistance/quarterDistance));
sascha@343: 
sascha@343:             double ikmw = factor*0.25*(mkmw-ws[i]) + ws[i];
sascha@343: 
sascha@343:             double end = ikm + quarterDistance*factor;
sascha@343: 
sascha@343:             double [] x = { mkm,  ikm,  end   };
sascha@343:             double [] y = { mkmw, ikmw, ws[i] };
sascha@343: 
sascha@343:             if (interpolator == null) {
sascha@343:                 interpolator = new SplineInterpolator();
sascha@343:             }
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@343:             PolynomialSplineFunction spline = interpolator.interpolate(x, y);
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@343:                 for (back = Math.max(back, 0);
sascha@343:                      back < i && km[back] < end; 
sascha@343:                      ++back
sascha@343:                 ) {
sascha@343:                     // to 3/4 point fill with spline values
sascha@343:                     ws[back] = spline.value(km[back]);
sascha@343:                 }
sascha@343:                 while (back < i) {
sascha@343:                     // fill the rest with ws[i]
sascha@343:                     ws[back++] = ws[i];
sascha@343:                 }
sascha@343:             }
sascha@343:             catch (ArgumentOutsideDomainException aode) {
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@642: 
sascha@642:     public static final boolean isIncreasing(double [] array) {
sascha@642:         int inc = 0;
sascha@642:         int dec = 0;
sascha@642:         int sweet = (array.length-1)/2;
sascha@642:         for (int i = 1; i < array.length; ++i) {
sascha@642:             if (array[i] > array[i-1]) {
sascha@642:                 if (++inc > sweet) {
sascha@642:                     return true;
sascha@642:                 }
sascha@642:             }
sascha@642:             else if (++dec > sweet) {
sascha@642:                 return false;
sascha@642:             }
sascha@642:         }
sascha@642:         return inc > sweet;
sascha@642:     }
sascha@642: 
sascha@642:     public static final double [] swap(double [] array) {
sascha@642:         int lo = 0;
sascha@642:         int hi = array.length-1;
sascha@642:         while (hi > lo) {
sascha@642:             double t  = array[lo];
sascha@642:             array[lo] = array[hi];
sascha@642:             array[hi] = t;
sascha@642:             ++lo;
sascha@642:             --hi;
sascha@642:         }
sascha@642: 
sascha@642:         return array;
sascha@642:     }
sascha@642: 
sascha@642:     public static final double [] swapClone(double [] array) {
sascha@642:         double [] out = new double[array.length];
sascha@642:         int lo = 0;
sascha@642: 
sascha@642:         int hi = array.length-1;
sascha@642:         while (hi > lo) {
sascha@642:             out[lo] = array[hi];
sascha@642:             out[hi] = array[lo];
sascha@642:             ++lo;
sascha@642:             --hi;
sascha@642:         }
sascha@642:         return out;
sascha@642:     }
sascha@642: 
sascha@642:     public static final double [] sumDiffs(double [] in) {
sascha@642:         double [] out = new double[in.length];
sascha@642: 
sascha@642:         if (out.length > 0) {
sascha@642:             out[0] = 0d;
sascha@642:         }
sascha@642: 
sascha@642:         for (int i = 1; i < out.length; ++i) {
sascha@642:             out[i] = out[i-1] + Math.abs(in[i-1] - in[i]);
sascha@642:         }
sascha@642: 
sascha@642:         return out;
sascha@642:     }
sascha@343: }
sascha@343: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :