changeset 655:913b52064449

Refactored version of "Berechnung 4" flys-artifacts/trunk@2053 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Sun, 05 Jun 2011 18:24:46 +0000
parents bbc966c81809
children baea7981477a
files flys-artifacts/ChangeLog flys-artifacts/src/main/java/de/intevation/flys/artifacts/FLYSArtifact.java flys-artifacts/src/main/java/de/intevation/flys/artifacts/WINFOArtifact.java flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/Function.java flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/Identity.java flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/Linear.java flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/LinearRemap.java flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/Calculation4.java flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/Segment.java flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WQCKms.java flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java flys-artifacts/src/main/java/de/intevation/flys/exports/ChartInfoGenerator.java flys-artifacts/src/main/java/de/intevation/flys/utils/DoubleUtil.java
diffstat 13 files changed, 621 insertions(+), 409 deletions(-) [+]
line wrap: on
line diff
--- a/flys-artifacts/ChangeLog	Fri Jun 03 10:21:13 2011 +0000
+++ b/flys-artifacts/ChangeLog	Sun Jun 05 18:24:46 2011 +0000
@@ -1,3 +1,51 @@
+2011-06-05  Sascha L. Teichmann <sascha.teichmann@intevation.de>
+
+	Refactored version of "Berechnung 4"
+
+	* src/main/java/de/intevation/flys/artifacts/model/Segment.java:
+	  Added instance fields for a reference point (= location of gauge)
+	  and backup of values (needed for naming).
+
+	* src/main/java/de/intevation/flys/artifacts/model/WQCKms.java:
+	  Added a constructor to be created from a WQKms. This is helpful
+	  if a WQKms is replaced by a back jump correction.
+
+	* src/main/java/de/intevation/flys/artifacts/model/Calculation4.java:
+	  New. Outfactored version of "W bei ungleichmaessigen Abflusslaengsschnitt".
+	  Much cleaner now and it should have a better handling of the corner
+	  cases.
+
+	* src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java:
+	  Removed the linear interpolation stuff. It is now in Linear. Removed
+	  the LinearRemap interpolation method because it is not needed any
+	  longer. Added a method to interpolate a given km with a given
+	  function.
+
+	* src/main/java/de/intevation/flys/artifacts/WINFOArtifact.java:
+	  Removed the old calc 4 and used the new one.
+
+	* src/main/java/de/intevation/flys/artifacts/math/LinearRemap.java:
+	  Deleted. Not needed any longer.
+
+	* src/main/java/de/intevation/flys/artifacts/math/Function.java:
+	  New. Interface for a uni-variate real function.
+
+	* src/main/java/de/intevation/flys/artifacts/math/Identity.java:
+	  New. Implements Function with f(x) = x
+
+	* src/main/java/de/intevation/flys/artifacts/math/Linear.java:
+	  New. Implements Function with f(x) = m*x + b
+
+	* src/main/java/de/intevation/flys/artifacts/FLYSArtifact.java:
+	  Factored some stuff out to DoubleUtil. Removed some dead code.
+	  Does some rounding correct.
+
+	* src/main/java/de/intevation/flys/utils/DoubleUtil.java: New.
+	  Centralized utils surrounding common double operations.
+
+	* src/main/java/de/intevation/flys/exports/ChartInfoGenerator.java:
+	  Removed superfluous imports.
+
 2011-06-03  Ingo Weinzierl <ingo@intevation.de>
 
 	  flys/issue90(Diagramm: Trennung der Diagrammfläche und Achsen aufheben)
--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/FLYSArtifact.java	Fri Jun 03 10:21:13 2011 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/FLYSArtifact.java	Sun Jun 05 18:24:46 2011 +0000
@@ -33,6 +33,8 @@
 import de.intevation.artifactdatabase.state.StateEngine;
 import de.intevation.artifactdatabase.transition.TransitionEngine;
 
+import de.intevation.flys.utils.DoubleUtil;
+
 import de.intevation.flys.model.Gauge;
 import de.intevation.flys.model.Range;
 import de.intevation.flys.model.River;
@@ -79,8 +81,6 @@
     /** The default step width between the start end end kilometer.*/
     public static final double DEFAULT_KM_STEPS = 0.1;
 
-    public static final double DEFAULT_PRECISION = 1e6;
-
 
     /** The identifier of the current state. */
     protected String currentStateId;
@@ -611,6 +611,36 @@
         }
     }
 
+    public double [] getFromToStep() {
+        if (!isRange()) {
+            return null;
+        }
+        double [] fromTo = getDistance();
+
+        if (fromTo == null) {
+            return null;
+        }
+
+        StateData dStep = getData("ld_step");
+        if (dStep == null) {
+            return null;
+        }
+
+        double [] result = new double[3];
+        result[0] = fromTo[0];
+        result[1] = fromTo[1];
+
+        try {
+            String step = (String)dStep.getValue();
+            result[2] = DoubleUtil.round(Double.parseDouble(step) / 1000d);
+        }
+        catch (NumberFormatException nfe) {
+            return null;
+        }
+
+        return result;
+    }
+
 
     /**
      * Returns the gauge based on the current distance and river.
@@ -967,36 +997,9 @@
         double to, 
         double step
     ) {
-        return getExplodedValues(from, to, step, DEFAULT_PRECISION);
+        return DoubleUtil.explode(from, to, step);
     }
 
-    public static double[] getExplodedValues(
-        double from, 
-        double to, 
-        double step,
-        double precision
-    ) {
-        double lower = from;
-
-        double diff = to - from;
-        double tmp  = diff / step;
-        int    num = (int)Math.abs(Math.ceil(tmp)) + 1;
-
-        double [] values = new double[num];
-
-        if (from > to) {
-            step = -step;
-        }
-
-        for (int idx = 0; idx < num; idx++) {
-            values[idx] = Math.round(lower * precision)/precision;
-            lower      += step;
-        }
-
-        return values;
-    }
-
-
     /**
      * Method to dump the artifacts state/data.
      */
--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/WINFOArtifact.java	Fri Jun 03 10:21:13 2011 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/WINFOArtifact.java	Sun Jun 05 18:24:46 2011 +0000
@@ -28,16 +28,15 @@
 
 import de.intevation.flys.artifacts.states.DefaultState;
 import de.intevation.flys.artifacts.context.FLYSContext;
-import de.intevation.flys.artifacts.math.BackJumpCorrector;
+
 import de.intevation.flys.artifacts.model.MainValuesFactory;
-import de.intevation.flys.artifacts.model.WQCKms;
 import de.intevation.flys.artifacts.model.WQDay;
 import de.intevation.flys.artifacts.model.WQKms;
 import de.intevation.flys.artifacts.model.WstValueTable;
 import de.intevation.flys.artifacts.model.WstValueTable.QPosition;
 import de.intevation.flys.artifacts.model.WstValueTableFactory;
-
-import de.intevation.flys.artifacts.math.LinearRemap;
+import de.intevation.flys.artifacts.model.Calculation4;
+import de.intevation.flys.artifacts.model.Segment;
 
 import gnu.trove.TDoubleArrayList;
 
@@ -583,204 +582,31 @@
             return new WQKms[0];
         }
 
-        boolean kmUp = river.getKmUp();
-
-        WstValueTable wst = WstValueTableFactory.getTable(river);
-        if (wst == null) {
+        WstValueTable table = WstValueTableFactory.getTable(river);
+        if (table == null) {
             logger.error("No wst found for selected river.");
             return new WQKms[0];
         }
 
-        double [][] segments = getRanges();
+        List<Segment> segments = getSegments();
 
-        if (logger.isDebugEnabled()) {
-            logger.debug("segments ----------------- enter");
-            for (double [] segment: segments) {
-                logger.debug("  " + joinDoubles(segment));
-            }
-            logger.debug("segments ----------------- leave");
-        }
-
-        if (segments.length < 1) {
-            logger.warn("no segments given");
+        if (segments == null) {
+            logger.error("Cannot create segments.");
             return new WQKms[0];
         }
 
-        if (segments.length == 1) {
-            // fall back to normal "Wasserstand/Wasserspiegellage" calculation
-            double [] qs = toQs(segments[0]);
-            if (qs == null) {
-                logger.warn("no qs given");
-                return new WQKms[0];
-            }
-            if (qs.length == 1) {
-                double [] kms = getKms(segments[0]);
-                HashSet<Integer> failed = new HashSet<Integer>();
-                WQKms [] results = computeWaterlevelData
-                    (kms, qs, wst, kmUp, failed);
-                setWaterlevelNames(
-                    results, qs, isQ() ? "Q" : "W", failed);
-                return results;
-            }
-        }
-
-        // more than one segment
-
-        double [] boundKms = extractBoundsKm(river, segments);
-
-
-        if (logger.isDebugEnabled()) {
-            logger.debug("bound kms: " + joinDoubles(boundKms));
-        }
-
-        double [][] iqs = null;
+        double [] range = getFromToStep();
 
-        for (int i = 0; i < segments.length; ++i) {
-            double [] iqsi = toQs(segments[i]);
-            if (iqsi == null) {
-                logger.warn("iqsi == null");
-                return new WQKms[0];
-            }
-
-            if (iqs == null) {
-                iqs = new double[iqsi.length][boundKms.length];
-            }
-            else if (iqs.length != iqsi.length) {
-                logger.warn("iqsi.logger != iqs.length: "
-                    + iqsi.length + " " + iqsi.length);
-                return new WQKms[0];
-            }
-
-            if (logger.isDebugEnabled()) {
-                logger.debug("segments qs[ " + i + "]: " + joinDoubles(iqsi));
-            }
-
-            for (int j = 0; j < iqs.length; ++j) {
-                iqs[j][i] = iqsi[j];
-            }
-        }
-
-        if (logger.isDebugEnabled()) {
-            for (int i = 0; i < iqs.length; ++i) {
-                logger.debug("iqs[" + i + "]: " + joinDoubles(iqs[i]));
-            }
+        if (range == null) {
+            logger.error("Cannot figure out range.");
+            return new WQKms[0];
         }
 
-        double [] boundWs = new double[boundKms.length];
-        double [] boundQs = new double[boundKms.length];
-
-        double [] okms = getKms(new double [] {
-            boundKms[0], boundKms[boundKms.length-1] });
-
-        ArrayList<WQKms> results = new ArrayList<WQKms>();
-
-        int referenceIndex = kmUp ?  0 : boundKms.length-1;
-
-        HashSet<Integer> failed = new HashSet<Integer>();
-
-        for (int i = 0; i < iqs.length; ++i) {
-            double [] iqsi = iqs[i];
-
-            QPosition qPosition = wst.interpolate(
-                iqsi[0],
-                boundKms[referenceIndex],
-                boundKms, boundWs, boundQs);
-
-            if (qPosition == null) {
-                logger.warn("interpolation failed for " + iqsi[i]);
-                failed.add(i);
-                continue;
-            }
-
-            LinearRemap remap = new LinearRemap();
-
-            for (int j = 1; j < boundKms.length; ++j) {
-                remap.add(
-                    boundKms[j-1], boundKms[j],
-                    boundQs[j-1],  iqsi[j-1],
-                    boundQs[j],    iqsi[j]);
-            }
-
-            double [] oqs = new double[okms.length];
-            double [] ows = new double[okms.length];
-
-            wst.interpolate(okms, ows, oqs, qPosition, remap);
-
-            BackJumpCorrector bjc = new BackJumpCorrector();
-            if (bjc.doCorrection(okms, ows)) {
-                logger.debug("Discharge longitudinal section has backjumps.");
-                results.add(new WQCKms(okms, oqs, ows, bjc.getCorrected()));
-            }
-            else {
-                results.add(new WQKms(okms, oqs, ows));
-            }
-        }
-
-        WQKms [] wqkms = results.toArray(new WQKms[results.size()]);
-
-        setDischargeLongitudinalSectionNames(
-            wqkms, iqs, isQ() ? "Q" : "W", failed);
+        Calculation4 calc4 = new Calculation4(segments, river, isQ());
 
-        return wqkms;
-    }
-
-    protected static String joinDoubles(double [] x) {
-        if (x == null) {
-            return "";
-        }
-        StringBuilder sb = new StringBuilder();
-        for (int i = 0; i < x.length; ++i) {
-            if (i > 0) sb.append(", ");
-            sb.append(x[i]);
-        }
-        return sb.toString();
-    }
-
-    protected double [] toQs(double [] range) {
-        double [] qs = getQs(range);
-        if (qs == null) {
-            logger.debug("Determine Q values based on a set of W values.");
-            double [] ws = getWs(range);
-            qs = getQsForWs(ws);
-        }
-        return qs;
-    }
-
+        WQKms [] results = calc4.calculate(table, range[0], range[1], range[2]);
 
-    /**
-     * Sets the name for discharge longitudinal section curves where each WQKms
-     * in <i>r</i> represents a column.
-     */
-    public static void setDischargeLongitudinalSectionNames(
-        WQKms  []    wqkms,
-        double [][]  iqs,
-        String       wq,
-        Set<Integer> failed
-    ) {
-        logger.debug("WINFOArtifact.setDischargeLongitudinalSectionNames");
-
-        // TODO: I18N
-
-        int pos = 0;
-
-        for (int j = 0; j < iqs.length; ++j) {
-            if (failed.contains(j)) {
-                continue;
-            }
-            StringBuilder sb = new StringBuilder(wq)
-                .append(" benutzerdefiniert (");
-
-            double [] iqsi = iqs[j];
-            for (int i = 0; i < iqsi.length; i++) {
-                if (i > 0) {
-                    sb.append("; ");
-                }
-                sb.append(iqsi[i]);
-            }
-            sb.append(")");
-
-            wqkms[pos++].setName(sb.toString());
-        }
+        return results;
     }
 }
 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/Function.java	Sun Jun 05 18:24:46 2011 +0000
@@ -0,0 +1,6 @@
+package de.intevation.flys.artifacts.math;
+
+public interface Function {
+    double value(double x);
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/Identity.java	Sun Jun 05 18:24:46 2011 +0000
@@ -0,0 +1,15 @@
+package de.intevation.flys.artifacts.math;
+
+public final class Identity 
+implements         Function
+{
+    public static final Identity IDENTITY = new Identity();
+
+    public Identity() {
+    }
+
+    public double value(double x) {
+        return x;
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/Linear.java	Sun Jun 05 18:24:46 2011 +0000
@@ -0,0 +1,68 @@
+package de.intevation.flys.artifacts.math;
+
+public final class Linear 
+implements         Function
+{
+    private double m;
+    private double b;
+
+    public Linear(
+        double x1, double x2,
+        double y1, double y2
+    ) {
+        // y1 = m*x1 + b
+        // y2 = m*x2 + b
+        // y2 - y1 = m*(x2 - x1)
+        // m = (y2 - y1)/(x2 - x1) # x2 != x1
+        // b = y1 - m*x1
+
+        if (x1 == x2) {
+            m = 0;
+            b = 0.5*(y1 + y2);
+        }
+        else {
+            m = (y2 - y1)/(x2 - x1);
+            b = y1 - m*x1;
+        }
+    }
+
+    public static final double linear(
+        double x,
+        double x1, double x2,
+        double y1, double y2
+    ) {
+        // y1 = m*x1 + b
+        // y2 = m*x2 + b
+        // y2 - y1 = m*(x2 - x1)
+        // m = (y2 - y1)/(x2 - x1) # x2 != x1
+        // b = y1 - m*x1
+
+        if (x1 == x2) {
+            return 0.5*(y1 + y2);
+        }
+        double m = (y2 - y1)/(x2 - x1);
+        double b = y1 - m*x1;
+        return x*m + b;
+    }
+
+    @Override
+    public double value(double x) {
+        return m*x + b;
+    }
+
+    public static final double factor(double x, double p1, double p2) {
+        // 0 = m*p1 + b <=> b = -m*p1
+        // 1 = m*p2 + b
+        // 1 = m*(p2 - p1)
+        // m = 1/(p2 - p1) # p1 != p2
+        // f(x) = x/(p2-p1) - p1/(p2-p1) <=> (x-p1)/(p2-p1)
+
+        return p1 == p2 ? 0.0 : (x-p1)/(p2-p1);
+    }
+
+    public static final double weight(double factor, double a, double b) {
+        //return (1.0-factor)*a + factor*b;
+        return a + factor*(b-a);
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/LinearRemap.java	Fri Jun 03 10:21:13 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,92 +0,0 @@
-package de.intevation.flys.artifacts.math;
-
-import org.apache.log4j.Logger;
-
-public class LinearRemap
-{
-    private static Logger logger = Logger.getLogger(LinearRemap.class);
-
-    public static class Segment {
-
-        protected Segment next;
-
-        protected double from;
-        protected double to;
-
-        protected double m;
-        protected double b;
-
-        public Segment() {
-        }
-
-        public Segment(
-            double from, double to, 
-            double m,    double b,
-            Segment next
-        ) {
-            this.from = from;
-            this.to   = to;
-            this.m    = m;
-            this.b    = b;
-        }
-
-        public double eval(double x) {
-            return m*x + b;
-        }
-    } // class Segment
-
-    protected Segment head;
-
-    public LinearRemap() {
-    }
-
-    public void add(
-        double from, double to,
-        double x1,   double y1,
-        double x2,   double y2
-    ) {
-        // y1 = m*x1 + b <=> b = y1 - m*x1
-        // y2 = m*x2 + b
-        // y2 - y1 = m*(x2 - x1)
-        // m = (y2 - y1)/(x2 - x1)
-
-        double m, b;
-
-        if (x2 == x1) {
-            m = 0.0;
-            b = 0.5*(y2 + y1);
-        }
-        else {
-            m = (y2 - y1)/(x2 - x1);
-            b = y1 - m*x1;
-        }
-
-        if (from > to) { double t = from; from = to; to = t; }
-
-        head = new Segment(from, to, m, b, head);
-
-        if (logger.isDebugEnabled()) {
-            logger.debug("LinearRemap.add --------- enter");
-            logger.debug("  range: [" + from + ", " + to + "]");
-            logger.debug("  " + x1 + " -> " + y1 + " (" + head.eval(x1) + ")");
-            logger.debug("  " + x2 + " -> " + y2 + " (" + head.eval(x2) + ")");
-            logger.debug("LinearRemap.add --------- leave");
-        }
-    }
-
-    public double eval(double pos, double x) {
-        Segment current = head;
-
-        while (current != null) {
-
-            if (pos >= current.from && pos <= current.to) {
-                return current.eval(x);
-            }
-
-            current = current.next;
-        }
-
-        return Double.NaN;
-    }
-}
-// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/Calculation4.java	Sun Jun 05 18:24:46 2011 +0000
@@ -0,0 +1,287 @@
+package de.intevation.flys.artifacts.model;
+
+import java.util.List;
+import java.util.Arrays;
+import java.util.Comparator;
+import java.util.Collections;
+
+import de.intevation.flys.utils.DoubleUtil;
+
+import de.intevation.flys.model.River;
+import de.intevation.flys.model.Gauge;
+
+import de.intevation.flys.artifacts.model.WstValueTable.QPosition;
+
+import de.intevation.flys.artifacts.math.Function;
+import de.intevation.flys.artifacts.math.Linear;
+import de.intevation.flys.artifacts.math.Identity;
+import de.intevation.flys.artifacts.math.BackJumpCorrector;
+
+import org.apache.log4j.Logger;
+
+public class Calculation4
+{
+    private static Logger logger = Logger.getLogger(Calculation4.class);
+
+    public static final double MINIMAL_STEP_WIDTH = 1e-5;
+
+    public static final Comparator<Segment> REF_CMP =
+        new Comparator<Segment>() {
+            public int compare(Segment a, Segment b) {
+                double d = a.referencePoint - b.referencePoint;
+                if (d < 0d) return -1;
+                return d > 0d ? +1 : 0;
+            }
+        };
+
+    protected List<Segment> segments;
+
+    protected boolean       isQ;
+
+    public Calculation4() {
+    }
+
+    public Calculation4(List<Segment> segments, River river, boolean isQ) {
+
+        this.segments = segments;
+        this.isQ      = isQ;
+
+        int numResults    = 1;
+
+        // assign reference points
+        for (Segment segment: segments) {
+            Gauge gauge = river.determineGauge(
+                segment.getFrom(), segment.getTo());
+
+            segment.setReferencePoint(gauge != null
+                ? gauge.getStation().doubleValue()
+                : 0.5*(segment.getFrom() + segment.getTo()));
+
+            double [] values = segment.values;
+
+            if (numResults == -1) {
+                numResults = values.length;
+            }
+            else if (numResults != values.length) {
+                throw new IllegalArgumentException("wrong length of values");
+            }
+
+            // convert to Q if needed
+            if (!isQ && gauge != null) {
+                double [][] table = new DischargeTables(
+                    river.getName(), gauge.getName()).getFirstTable();
+
+                // need the original values for naming
+                segment.backup();
+
+                for (int i = 0; i < values.length; ++i) {
+                    values[i] = DischargeTables.getQForW(table, values[i]);
+                }
+            }
+        } // for all segments
+
+        Collections.sort(segments, REF_CMP);
+    }
+
+    public WQKms [] calculate(
+        WstValueTable table,
+        double from, double to, double step
+    ) {
+        boolean debug = logger.isDebugEnabled();
+
+        if (debug) {
+            logger.debug(
+                "calculate from " + from + " to " + to + " step " + step);
+            logger.debug("# segments: " + segments.size());
+        }
+
+        if (segments.isEmpty()) {
+            logger.debug("no segments found");
+            return new WQKms[0];
+        }
+
+        int numResults = segments.get(0).values.length;
+
+        if (numResults < 1) {
+            logger.debug("no values given");
+            return new WQKms[0];
+        }
+
+
+        WQKms [] results = new WQKms[numResults];
+                    if (debug) {
+                        logger.debug("after last segment -> gleichwertig");
+                    }
+        for (int i = 0; i < results.length; ++i) {
+            results[i] = new WQKms();
+        }
+
+        if (Math.abs(step) < MINIMAL_STEP_WIDTH) {
+            step = MINIMAL_STEP_WIDTH;
+        }
+
+        if (from > to) {
+            step = -step;
+        }
+
+        QPosition [] qPositions = new QPosition[numResults];
+
+        Function [] functions = new Function[numResults];
+
+        double [] out = new double[2];
+
+        Segment sentinel = new Segment(Double.MAX_VALUE);
+        Segment s1 = sentinel, s2 = sentinel;
+
+        for (double pos = from;
+             from < to ? pos <= to : pos >= to;
+             pos = DoubleUtil.round(pos + step)
+        ) {
+            if (pos < s1.referencePoint || pos > s2.referencePoint) {
+                if (debug) {
+                    logger.debug("need to find new interval for " + pos);
+                }
+                // find new interval
+                if (pos <= segments.get(0).referencePoint) {
+                    // before first segment -> "gleichwertig"
+                    if (debug) {
+                        logger.debug("before first segment -> gleichwertig");
+                    }
+                    Segment   first  = segments.get(0);
+                    double [] values = first.values;
+                    double    refPos = first.referencePoint;
+                    for (int i = 0; i < qPositions.length; ++i) {
+                        qPositions[i] = table.getQPosition(
+                            refPos, values[i]);
+                    }
+                    sentinel.setReferencePoint(-Double.MAX_VALUE);
+                    s1 = sentinel;
+                    s2 = segments.get(0);
+                    Arrays.fill(functions, Identity.IDENTITY);
+                }
+                else if (pos >= segments.get(segments.size()-1).referencePoint) {
+                    // after last segment -> "gleichwertig"
+                    if (debug) {
+                        logger.debug("after last segment -> gleichwertig");
+                    }
+                    Segment   last   = segments.get(segments.size()-1);
+                    double [] values = last.values;
+                    double    refPos = last.referencePoint;
+                    for (int i = 0; i < qPositions.length; ++i) {
+                        qPositions[i] = table.getQPosition(
+                            refPos, values[i]);
+                    }
+                    sentinel.setReferencePoint(Double.MAX_VALUE);
+                    s1 = last;
+                    s2 = sentinel;
+                    Arrays.fill(functions, Identity.IDENTITY);
+                }
+                else { // "ungleichwertig"
+                    // find matching interval
+                    if (debug) {
+                        logger.debug("in segments -> ungleichwertig");
+                    }
+                    s1 = s2 = null;
+                    for (int i = 1, N = segments.size(); i < N; ++i) {
+                        Segment si1 = segments.get(i-1);
+                        Segment si  = segments.get(i);
+                        if (debug) {
+                            logger.debug("check " + pos + " in " + 
+                                si1.referencePoint + " - " + si.referencePoint);
+                        }
+                        if (pos >= si1.referencePoint 
+                        &&  pos <= si. referencePoint) {
+                            s1 = si1;
+                            s2 = si;
+                            break;
+                        }
+                    }
+
+                    if (s1 == null) {
+                        throw new IllegalStateException("no interval found");
+                    }
+
+                    Segment anchor, free;
+
+                    if (from > to) { anchor = s1; free = s2; }
+                    else           { anchor = s2; free = s1; }
+
+                    // build transforms based on "gleichwertiger" phase
+                    for (int i = 0; i < qPositions.length; ++i) {
+                        QPosition qi = table.getQPosition(
+                            anchor.referencePoint,
+                            anchor.values[i]);
+
+                        if ((qPositions[i] = qi) == null) {
+                            // TODO: error report
+                            functions[i] = Identity.IDENTITY;
+                        }
+                        else {
+                            double qA = table.getQ(qi, anchor.referencePoint);
+                            double qF = table.getQ(qi, free  .referencePoint);
+
+                            functions[i] = Double.isNaN(qA) || Double.isNaN(qF)
+                                ? Identity.IDENTITY
+                                : new Linear(
+                                    qA, qF,
+                                    anchor.values[i], free.values[i]);
+                        }
+                    } // build transforms
+                } // "ungleichwertiges" interval
+            } // find matching interval
+
+            for (int i = 0; i < qPositions.length; ++i) {
+                QPosition qPosition = qPositions[i];
+
+                if (qPosition == null) {
+                    continue;
+                }
+
+                if (table.interpolate(pos, out, qPosition, functions[i])) {
+                    results[i].add(out[0], out[1], pos);
+                }
+                else {
+                    // TODO: error report
+                }
+            }
+        }
+
+        // Backjump correction
+        for (int i = 0; i < results.length; ++i) {
+            BackJumpCorrector bjc = new BackJumpCorrector();
+
+            double [] ws  = results[i].getWs();
+            double [] kms = results[i].getKms();
+
+            if (bjc.doCorrection(kms, ws)) {
+                results[i] = new WQCKms(results[i], bjc.getCorrected());
+            }
+            // TODO: error report
+        }
+
+        // name the curves
+        for (int i = 0; i < results.length; ++i) {
+            results[i].setName(createName(i));
+        }
+
+        return results;
+    }
+
+    protected String createName(int index) {
+        // TODO: i18n
+        StringBuilder sb = new StringBuilder(isQ ? "Q" : "W");
+        sb.append(" benutzerdefiniert (");
+        for (int i = 0, N = segments.size(); i < N; ++i) {
+            if (i > 0) {
+                sb.append("; ");
+            }
+            Segment segment = segments.get(i);
+            sb.append((segment.backup != null
+                ? segment.backup
+                : segment.values)[index]);
+        }
+        sb.append(')');
+        return sb.toString();
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/Segment.java	Fri Jun 03 10:21:13 2011 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/Segment.java	Sun Jun 05 18:24:46 2011 +0000
@@ -9,6 +9,8 @@
 
 import gnu.trove.TDoubleArrayList;
 
+import de.intevation.flys.utils.DoubleUtil;
+
 public class Segment
 implements   Serializable
 {
@@ -17,10 +19,16 @@
     protected double    from;
     protected double    to;
     protected double [] values;
+    protected double [] backup;
+    protected double    referencePoint;
 
     public Segment() {
     }
 
+    public Segment(double referencePoint) {
+        this.referencePoint = referencePoint;
+    }
+
     public Segment(double from, double to, double [] values) {
         this.from   = from;
         this.to     = to;
@@ -43,6 +51,42 @@
         return sb.toString();
     }
 
+    public void setFrom(double from) {
+        this.from = from;
+    }
+
+    public void backup() {
+        backup = (double [])values.clone();
+    }
+
+    public double getFrom() {
+        return from;
+    }
+
+    public void setTo(double to) {
+        this.to = to;
+    }
+
+    public double getTo() {
+        return to;
+    }
+
+    public void setValues(double [] values) {
+        this.values = values;
+    }
+
+    public double [] getValues() {
+        return values;
+    }
+
+    public void setReferencePoint(double referencePoint) {
+        this.referencePoint = referencePoint;
+    }
+
+    public double getReferencePoint() {
+        return referencePoint;
+    }
+
     public static List<Segment> parseSegments(String input) {
 
         ArrayList<Segment> segments = new ArrayList<Segment>();
@@ -62,7 +106,8 @@
                 vs.clear();
 
                 for (String valueStr: parts[2].split(",")) {
-                    vs.add(Double.parseDouble(valueStr.trim()));
+                    vs.add(DoubleUtil.round(
+                        Double.parseDouble(valueStr.trim())));
                 }
 
                 double [] values = vs.toNativeArray();
--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WQCKms.java	Fri Jun 03 10:21:13 2011 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WQCKms.java	Sun Jun 05 18:24:46 2011 +0000
@@ -13,6 +13,16 @@
 
     protected TDoubleArrayList cw;
 
+    public WQCKms() {
+    }
+
+    public WQCKms(WQKms other, double [] cws) {
+        this.w   = other.w;
+        this.q   = other.q;
+        this.kms = other.kms;
+        this.cw  = new TDoubleArrayList(cws);
+    }
+
 
     public WQCKms(double[] kms, double[] qs, double[] ws, double[] cws) {
         super(kms, qs, ws);
--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java	Fri Jun 03 10:21:13 2011 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java	Sun Jun 05 18:24:46 2011 +0000
@@ -2,8 +2,8 @@
 
 import java.io.Serializable;
 
-
-import de.intevation.flys.artifacts.math.LinearRemap;
+import de.intevation.flys.artifacts.math.Linear;
+import de.intevation.flys.artifacts.math.Function;
 
 import java.util.Arrays;
 import java.util.ArrayList;
@@ -109,7 +109,7 @@
             double []     ows,
             WstValueTable table
         ) {
-            double kmWeight = factor(km, this.km, other.km);
+            double kmWeight = Linear.factor(km, this.km, other.km);
 
             QPosition qPosition = new QPosition();
 
@@ -117,7 +117,7 @@
                 if (table.getQPosition(km, iqs[i], qPosition) != null) {
                     double wt =       getW(qPosition);
                     double wo = other.getW(qPosition);
-                    ows[i] = weight(kmWeight, wt, wo);
+                    ows[i] = Linear.weight(kmWeight, wt, wo);
                 }
                 else {
                     ows[i] = Double.NaN;
@@ -137,7 +137,7 @@
                 return new double[2][0];
             }
 
-            double factor = factor(km, this.km, other.km);
+            double factor = Linear.factor(km, this.km, other.km);
 
             double [] splineQ = new double[W];
             double [] splineW = new double[W];
@@ -146,8 +146,8 @@
             double maxQ = -Double.MAX_VALUE;
 
             for (int i = 0; i < W; ++i) {
-                double wws = weight(factor, ws[i], other.ws[i]);
-                double wqs = weight(
+                double wws = Linear.weight(factor, ws[i], other.ws[i]);
+                double wqs = Linear.weight(
                     factor,
                     table.getQIndex(i, km),
                     table.getQIndex(i, other.km));
@@ -231,7 +231,7 @@
 
             return weight == 1.0
                 ? ws[index]
-                : weight(weight, ws[index-1], ws[index]);
+                : Linear.weight(weight, ws[index-1], ws[index]);
         }
 
         public double getW(
@@ -239,7 +239,7 @@
             double    km,
             QPosition qPosition
         ) {
-            double kmWeight = factor(km, this.km, other.km);
+            double kmWeight = Linear.factor(km, this.km, other.km);
 
             int    index  = qPosition.index;
             double weight = qPosition.weight;
@@ -251,11 +251,11 @@
                 ow = other.ws[index];
             }
             else {
-                tw = weight(weight, ws[index-1], ws[index]); 
-                ow = weight(weight, other.ws[index-1], other.ws[index]); 
+                tw = Linear.weight(weight, ws[index-1], ws[index]); 
+                ow = Linear.weight(weight, other.ws[index-1], other.ws[index]); 
             }
 
-            return weight(kmWeight, tw, ow);
+            return Linear.weight(kmWeight, tw, ow);
         }
     } // class Row
 
@@ -338,60 +338,45 @@
         return r1.interpolateWQ(r2, km, steps, this);
     }
 
-    public void interpolate(
-        double []   kms,
-        double []   ws,
-        double []   qs,
-        QPosition   qPosition,
-        LinearRemap remap
-    ) {
-        interpolate(kms, ws, qs, qPosition, remap, 0, kms.length);
-    }
-
-    public void interpolate(
-        double []   kms,
-        double []   ws,
-        double []   qs,
-        QPosition   qPosition,
-        LinearRemap remap,
-        int         startIndex,
-        int         length
+    public boolean interpolate(
+        double    km,
+        double [] out,
+        QPosition qPosition,
+        Function  qFunction
     ) {
         int R1 = rows.size()-1;
 
-        Row kmKey = new Row();
+        out[1] = qFunction.value(getQ(qPosition, km));
+
+        if (Double.isNaN(out[1])) {
+            return false;
+        }
 
         QPosition nPosition = new QPosition();
-
-        for (int i = startIndex, end = startIndex+length; i < end; ++i) {
-            kmKey.km = kms[i];
-
-            qs[i] = remap.eval(kms[i], getQ(qPosition, kms[i]));
-
-            if (getQPosition(kms[i], qs[i], nPosition) == null) {
-                ws[i] = Double.NaN;
-                continue;
-            }
-
-            int rowIndex = Collections.binarySearch(rows, kmKey);
+        if (getQPosition(km, out[1], nPosition) == null) {
+            return false;
+        }
 
-            if (rowIndex >= 0) {
-                // direct row match
-                ws[i] = rows.get(rowIndex).getW(nPosition);
-                continue;
-            }
-
-            rowIndex = -rowIndex -1;
+        int rowIndex = Collections.binarySearch(rows, new Row(km));
 
-            if (rowIndex < 1 || rowIndex >= R1) {
-                // do not extrapolate
-                ws[i] = Double.NaN;
-                continue;
-            }
-            Row r1 = rows.get(rowIndex-1);
-            Row r2 = rows.get(rowIndex);
-            ws[i] = r1.getW(r2, kms[i], nPosition);
+        if (rowIndex >= 0) {
+            // direct row match
+            out[0] = rows.get(rowIndex).getW(nPosition);
+            return !Double.isNaN(out[0]);
         }
+
+        rowIndex = -rowIndex -1;
+
+        if (rowIndex < 1 || rowIndex >= R1) {
+            // do not extrapolate
+            return false;
+        }
+
+        Row r1 = rows.get(rowIndex-1);
+        Row r2 = rows.get(rowIndex);
+        out[0] = r1.getW(r2, km, nPosition);
+
+        return !Double.isNaN(out[0]);
     }
 
     public QPosition interpolate(
@@ -480,7 +465,7 @@
             else                  { qMin = qCurrent; qMax = qLast; }
 
             if (q > qMin && q < qMax) {
-                double weight = factor(q, qLast, qCurrent);
+                double weight = Linear.factor(q, qLast, qCurrent);
                 return qPosition.set(i, weight);
             }
             qLast = qCurrent;
@@ -502,41 +487,8 @@
         }
         double q1 = columns[index-1].getQRangeTree().findQ(km);
         double q2 = columns[index  ].getQRangeTree().findQ(km);
-        return weight(weight, q1, q2);
+        return Linear.weight(weight, q1, q2);
     }
 
-    public static final double linear(
-        double x,
-        double x1, double x2,
-        double y1, double y2
-    ) {
-        // y1 = m*x1 + b
-        // y2 = m*x2 + b
-        // y2 - y1 = m*(x2 - x1)
-        // m = (y2 - y1)/(x2 - x1) # x2 != x1
-        // b = y1 - m*x1
-
-        if (x1 == x2) {
-            return 0.5*(y1 + y2);
-        }
-        double m = (y2 - y1)/(x2 - x1);
-        double b = y1 - m*x1;
-        return x*m + b;
-    }
-
-    public static final double factor(double x, double p1, double p2) {
-        // 0 = m*p1 + b <=> b = -m*p1
-        // 1 = m*p2 + b
-        // 1 = m*(p2 - p1)
-        // m = 1/(p2 - p1) # p1 != p2
-        // f(x) = x/(p2-p1) - p1/(p2-p1) <=> (x-p1)/(p2-p1)
-
-        return p1 == p2 ? 0.0 : (x-p1)/(p2-p1);
-    }
-
-    public static final double weight(double factor, double a, double b) {
-        //return (1.0-factor)*a + factor*b;
-        return a + factor*(b-a);
-    }
 }
 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
--- a/flys-artifacts/src/main/java/de/intevation/flys/exports/ChartInfoGenerator.java	Fri Jun 03 10:21:13 2011 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/exports/ChartInfoGenerator.java	Sun Jun 05 18:24:46 2011 +0000
@@ -1,6 +1,5 @@
 package de.intevation.flys.exports;
 
-import java.awt.Color;
 import java.awt.Transparency;
 import java.io.IOException;
 import java.io.OutputStream;
@@ -9,11 +8,8 @@
 
 import org.apache.log4j.Logger;
 
-import org.jfree.chart.ChartFactory;
 import org.jfree.chart.ChartRenderingInfo;
 import org.jfree.chart.JFreeChart;
-import org.jfree.chart.plot.PlotOrientation;
-import org.jfree.chart.plot.XYPlot;
 
 import de.intevation.artifacts.Artifact;
 import de.intevation.artifacts.CallContext;
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/utils/DoubleUtil.java	Sun Jun 05 18:24:46 2011 +0000
@@ -0,0 +1,48 @@
+package de.intevation.flys.utils;
+
+public class DoubleUtil
+{
+    public static final double DEFAULT_STEP_PRECISION = 1e6;
+
+    private DoubleUtil() {
+    }
+
+    public static final double [] explode(double from, double to, double step) {
+        return explode(from, to, step, DEFAULT_STEP_PRECISION);
+    }
+
+    public static final double round(double x, double precision) {
+        return Math.round(x * precision)/precision;
+    }
+
+    public static final double round(double x) {
+        return Math.round(x * DEFAULT_STEP_PRECISION)/DEFAULT_STEP_PRECISION;
+    }
+
+    public static final double [] explode(
+        double from,
+        double to,
+        double step,
+        double precision
+    ) {
+        double lower = from;
+
+        double diff = to - from;
+        double tmp  = diff / step;
+        int    num = (int)Math.abs(Math.ceil(tmp)) + 1;
+
+        double [] values = new double[num];
+
+        if (from > to) {
+            step = -step;
+        }
+
+        for (int idx = 0; idx < num; idx++) {
+            values[idx] = round(lower, precision);
+            lower      += step;
+        }
+
+        return values;
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :

http://dive4elements.wald.intevation.org