# HG changeset patch
# User Sascha L. Teichmann <sascha.teichmann@intevation.de>
# Date 1303466247 0
# Node ID f165c7d5d6db0fcc1d5bded0365d6655b5b80f26
# Parent  f72c63713099d56959201526d59be595ddfdde00
Implementation of the "Ruecksprungkorrektur" to be done in "W fuer angepassten Abflusslaengschnitt".



flys-artifacts/trunk@1743 c6561f87-3c4e-4783-a992-168aeb5c3f6f

diff -r f72c63713099 -r f165c7d5d6db flys-artifacts/ChangeLog
--- a/flys-artifacts/ChangeLog	Thu Apr 21 08:49:58 2011 +0000
+++ b/flys-artifacts/ChangeLog	Fri Apr 22 09:57:27 2011 +0000
@@ -1,3 +1,34 @@
+2011-04-22	Sascha L. Teichmann	<sascha.teichmann@intevation.de>
+
+	* src/main/java/de/intevation/flys/artifacts/math/BackJumpCorrector.java:
+
+	  Implementation of the "Ruecksprungkorrektur" to be done in
+	  "W fuer angepassten Abflusslaengschnitt".
+
+	  All tests show the expected results. In some corner cases the
+	  algorithm described in the "Anwenderhandbuch" chapter 3.3.4.3 "Korrektur"
+	  has some definition shortcomings:
+
+	  a - What should happend when you cannot find point 2 because
+	      you cannot step back one quarter from point 3 because there
+		  is no data there any more (river too short in this direction)?
+		  The implemented algorithm raises point 3' only to an
+		  according factor. E.g. If you can step back the whole quarter
+		  distance the elevation is the full quarter. If you can
+		  step back only the half of the quarter the elevation is
+		  only an eighth.
+
+	  b - If the water heights between point 2 and 3 are constant then
+	      the algorithm will produce a spline interpolation that
+		  lowers those values. Is this intended?
+
+	  For real data the back jumps are expected to be more in the middle
+	  of the distance ranges so the corner cases are maybe not so
+	  important.
+
+	* src/main/java/de/intevation/flys/artifacts/states/RiverSelect.java:
+	  Removed superfluous import.
+
 2011-04-21  Ingo Weinzierl <ingo@intevation.de>
 
 	* src/main/java/de/intevation/flys/themes/ThemeFactory.java: Removed debug
diff -r f72c63713099 -r f165c7d5d6db flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/BackJumpCorrector.java
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/BackJumpCorrector.java	Fri Apr 22 09:57:27 2011 +0000
@@ -0,0 +1,188 @@
+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.log4j.Logger;
+
+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) {
+        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();
+            }
+
+            PolynomialSplineFunction spline = interpolator.interpolate(x, y);
+
+            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) {
+                log.error("spline interpolation failed", aode);
+            }
+        } // for all km
+
+        return !backjumps.isEmpty();
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
diff -r f72c63713099 -r f165c7d5d6db flys-artifacts/src/main/java/de/intevation/flys/artifacts/states/RiverSelect.java
--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/states/RiverSelect.java	Thu Apr 21 08:49:58 2011 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/states/RiverSelect.java	Fri Apr 22 09:57:27 2011 +0000
@@ -18,7 +18,6 @@
 
 import de.intevation.flys.model.River;
 
-import de.intevation.flys.artifacts.FLYSArtifact;
 import de.intevation.flys.artifacts.model.RiverFactory;
 import de.intevation.flys.artifacts.resources.Resources;