changeset 2418:899ca89f497e

Another partial fix for flys/issue499: Do the W to Q conversions needed for 'W am Pegel' correctly. flys-artifacts/trunk@4052 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Tue, 14 Feb 2012 16:48:13 +0000
parents e5fa3cbbe3ae
children 98a350bb91a9
files flys-artifacts/ChangeLog flys-artifacts/src/main/java/de/intevation/flys/artifacts/WINFOArtifact.java flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/Calculation4.java flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/Calculation6.java flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/DischargeTables.java
diffstat 5 files changed, 123 insertions(+), 101 deletions(-) [+]
line wrap: on
line diff
--- a/flys-artifacts/ChangeLog	Tue Feb 14 09:02:51 2012 +0000
+++ b/flys-artifacts/ChangeLog	Tue Feb 14 16:48:13 2012 +0000
@@ -1,3 +1,25 @@
+2012-02-14  Sascha L. Teichmann	<sascha.teichmann@intevation.de>
+
+	Another partial fix for flys/issue499: Do the W to Q conversions
+	needed for "W am Pegel" correctly.
+
+	* src/main/java/de/intevation/flys/artifacts/model/DischargeTables.java:
+	  Repaired getQsForW(): The mapping from W to Q is not unique! There
+	  could be more then one Q having the the same W.
+	  Ws are not strictly monoton/sorted so doing a binary search on this
+	  is just a fail. We now scan them linearly.
+
+	  XXX: The whole class is mess. The scaling stuff is a stupid
+	  and there is no caching.
+
+	* src/main/java/de/intevation/flys/artifacts/WINFOArtifact.java:
+	  Fetch the master discharge table for converting Ws to Qs. Handle
+	  the case that there are more Qs for a given W.
+	  
+	* src/main/java/de/intevation/flys/artifacts/model/Calculation6.java,
+	  src/main/java/de/intevation/flys/artifacts/model/Calculation4.java:
+	  Adjusted to new semantic.
+
 2012-02-13	Felix Wolfsteller	<felix.wolfsteller@intevation.de>
 
 	Partial Fix flys/issue500: text-orientation for texts.
--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/WINFOArtifact.java	Tue Feb 14 09:02:51 2012 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/WINFOArtifact.java	Tue Feb 14 16:48:13 2012 +0000
@@ -48,6 +48,7 @@
 import de.intevation.flys.model.FastCrossSectionLine;
 import de.intevation.flys.model.Gauge;
 import de.intevation.flys.model.River;
+import de.intevation.flys.model.DischargeTable;
 
 import de.intevation.flys.utils.DoubleUtil;
 import de.intevation.flys.utils.FLYSUtils;
@@ -925,36 +926,44 @@
             logger.debug("convert w->q with gauge '" + g.getName() + "'");
         }
 
-        DischargeTables dt = new DischargeTables(r.getName(), g.getName());
-        Map<String, double [][]>  tmp = dt.getValues();
+        DischargeTable dt = g.fetchMasterDischargeTable();
 
-        double[][] values = tmp.get(g.getName());
-        double[]   qs     = new double[ws.length];
+        if (dt == null) {
+            logger.warn("No master discharge table found for gauge '"
+                + g.getName() + "'");
+            return null;
+        }
+
+        double [][] values = DischargeTables.loadDischargeTableValues(dt, 1);
 
         TDoubleArrayList wsOut = new TDoubleArrayList(ws.length);
         TDoubleArrayList qsOut = new TDoubleArrayList(ws.length);
 
+        boolean generatedWs = false;
+
         for (int i = 0; i < ws.length; i++) {
             if (Double.isNaN(ws[i])) {
                 logger.warn("W is NaN: ignored");
                 continue;
             }
-            double w = ws[i] / 100.0;
-            double q = dt.getQForW(values, w);
-            if (Double.isNaN(q)) {
-                logger.warn("No Q found for W = " + ws[i]);
-                continue;
+            double w = ws[i] / 100d;
+            double [] qs = DischargeTables.getQsForW(values, w);
+
+            if (qs.length == 0) {
+                logger.warn("No Qs found for W = " + ws[i]);
             }
-            wsOut.add(w);
-            qsOut.add(q);
-            if (debug) {
-                logger.debug("w: " + w + " -> q: " + q);
+            else {
+                for (double q: qs) {
+                    wsOut.add(ws[i]);
+                    qsOut.add(q);
+                }
             }
+            generatedWs |= qs.length != 1;
         }
 
         return new double [][] {
             qsOut.toNativeArray(),
-            wsOut.toNativeArray()
+            generatedWs ? wsOut.toNativeArray() : null
         };
     }
 
--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/Calculation4.java	Tue Feb 14 09:02:51 2012 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/Calculation4.java	Tue Feb 14 16:48:13 2012 +0000
@@ -84,7 +84,19 @@
                 segment.backup();
 
                 for (int i = 0; i < values.length; ++i) {
-                    values[i] = DischargeTables.getQForW(table, values[i]);
+                    double w = values[i] * 100;
+                    double [] qs = DischargeTables.getQsForW(table, w);
+                    if (qs.length == 0) {
+                        logger.warn("No Qs found for W = " + values[i]);
+                        values[i] = Double.NaN;
+                    }
+                    else {
+                        values[i] = qs[0];
+                        if (qs.length > 1) {
+                            logger.warn(
+                                "More than one Q found for W = " + values[i]);
+                        }
+                    }
                 }
             }
         } // for all segments
--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/Calculation6.java	Tue Feb 14 09:02:51 2012 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/Calculation6.java	Tue Feb 14 16:48:13 2012 +0000
@@ -282,13 +282,14 @@
 
 
     protected double findValueForW(DischargeTable dt, double w) {
-        double[][] vs = DischargeTables.loadDischargeTableValues(dt,SCALE,true);
-        return DischargeTables.getQForW(vs, w);
+        double[][] vs = DischargeTables.loadDischargeTableValues(dt, SCALE);
+        double [] qs = DischargeTables.getQsForW(vs, w);
+        return qs.length == 0 ? Double.NaN : qs[0];
     }
 
 
     protected double findValueForQ(DischargeTable dt, double q) {
-        double[][] vs = DischargeTables.loadDischargeTableValues(dt,SCALE,true);
+        double[][] vs = DischargeTables.loadDischargeTableValues(dt, SCALE);
         logger.warn("TODO: IMPLEMENT ME!");
 
         return 10;
--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/DischargeTables.java	Tue Feb 14 09:02:51 2012 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/DischargeTables.java	Tue Feb 14 16:48:13 2012 +0000
@@ -4,7 +4,6 @@
 import java.util.Map;
 import java.util.HashMap;
 import java.util.Arrays;
-import java.util.Comparator;
 
 import java.io.Serializable;
 
@@ -18,6 +17,8 @@
 import de.intevation.flys.model.DischargeTable;
 import de.intevation.flys.model.DischargeTableValue;
 
+import gnu.trove.TDoubleArrayList;
+
 public class DischargeTables
 implements   Serializable
 {
@@ -124,42 +125,9 @@
             // TODO: Filter by time interval
             DischargeTable table = tables.get(0);
 
-            boolean qSorted = true;
-
-            final double[][] vs =
-                loadDischargeTableValues(table, scale, qSorted);
-
-            if (!qSorted) {
-                log.debug("need to sort by q values");
-                // TODO: Do this db level.
-                // XXX: This is so ugly :-(
-                Integer [] indices = new Integer[vs[0].length];
-                for (int i = 0; i < indices.length; ++i) {
-                    indices[i] = i;
-                }
+            double [][] vs = loadDischargeTableValues(table, scale);
 
-                Arrays.sort(indices, new Comparator<Integer>() {
-                    public int compare(Integer a, Integer b) {
-                        double va = vs[1][a];
-                        double vb = vs[1][b];
-                        double d = va - vb;
-                        if (d < 0.0) return -1;
-                        if (d > 0.0) return +1;
-                        return 0;
-                    }
-                });
-
-                double [][] vs2 = new double[2][vs[0].length];
-                for (int i = 0; i < indices.length; ++i) {
-                    vs2[0][i] = vs[0][indices[i]];
-                    vs2[1][i] = vs[1][indices[i]];
-                }
-                values.put(gaugeName, vs2);
-            }
-            else {
-                values.put(gaugeName, vs);
-            }
-
+            values.put(gaugeName, vs);
         }
 
         return values;
@@ -176,82 +144,92 @@
      */
     public static double[][] loadDischargeTableValues(
         DischargeTable table,
-        double         scale,
-        boolean        qSorted
+        double         scale
     ) {
         List<DischargeTableValue> dtvs = table.getDischargeTableValues();
 
         final double [][] vs = new double[2][dtvs.size()];
 
-        double lastQ = -Double.MAX_VALUE;
-
         int idx = 0;
         for (DischargeTableValue dtv: dtvs) {
             double q = dtv.getQ().doubleValue();
             vs[0][idx] = q * scale;
             vs[1][idx] = dtv.getW().doubleValue() * scale;
             ++idx;
-
-            if (qSorted && lastQ > q) {
-                qSorted = false;
-            }
-            lastQ = q;
         }
 
         return vs;
     }
 
+    private static final double EPSILON = 1e-5;
 
-    public static double getQForW(double [][] values, double w) {
+    private static final boolean epsEquals(double a, double b) {
+        return Math.abs(a - b) < EPSILON;
+    }
+
+    private static final boolean between(double a, double b, double x) {
+        if (a > b) { double t = a; a = b; b = t; }
+        return x > a && x < b;
+    }
+
+    public static double [] getQsForW(double [][] values, double w) {
 
         boolean debug = log.isDebugEnabled();
 
         if (debug) {
-            log.debug("calculating getQForW(" + w + ")");
-        }
-
-        int index = Arrays.binarySearch(values[1], w);
-        if (index >= 0) {
-            return values[0][index];
-        }
-
-        index = -index - 1; // insert position
-
-        if (index < 1 || index >= values[0].length) {
-            // do not extraploate
-            if (debug) {
-                log.debug("we do not extrapolate: NaN");
-            }
-            return Double.NaN;
+            log.debug("getQsForW: W = " + w);
         }
 
-        double w1 = values[1][index-1];
-        double w2 = values[1][index  ];
-        double q1 = values[0][index-1];
-        double q2 = values[0][index  ];
+        double [] qs = values[0];
+        double [] ws = values[1];
 
-        // q1 = m*w1 + b
-        // q2 = m*w2 + b
-        // q2 - q1 = m*(w2 - w1)
-        // m = (q2 - q1)/(w2 - w1) # w2 != w1
-        // b = q1 - m*w1
+        int N = Math.min(qs.length, ws.length);
 
-        double q;
-        if (w1 == w2) {
-            q = 0.5*(q1 + q2);
+        if (N == 0) {
             if (debug) {
-                log.debug("same w1 and w2: " + w1);
+                log.debug("Q(" + w + ") = []");
+            }
+            return new double [0];
+        }
+
+        TDoubleArrayList outQs = new TDoubleArrayList();
+
+        if (epsEquals(ws[0], w)) {
+            outQs.add(qs[0]);
+        }
+
+        for (int i = 1; i < N; ++i) {
+            if (epsEquals(ws[i], w)) {
+                outQs.add(qs[i]);
+            }
+            else if (between(ws[i-1], ws[i], w)) {
+                double w1 = ws[i-1];
+                double w2 = ws[i];
+                double q1 = qs[i-1];
+                double q2 = qs[i];
+
+                // q1 = m*w1 + b
+                // q2 = m*w2 + b
+                // q2 - q1 = m*(w2 - w1)
+                // m = (q2 - q1)/(w2 - w1) # w2 != w1
+                // b = q1 - m*w1
+                // w1 != w2
+
+                double m = (q2 - q1)/(w2 - w1);
+                double b = q1 - m*w1;
+                double q = w*m + b;
+
+                outQs.add(q);
             }
         }
-        else {
-            double m = (q2 - q1)/(w2 - w1);
-            double b = q1 - m*w1;
-            q = w*m + b;
+
+        double [] result = outQs.toNativeArray();
+
+        if (debug) {
+            log.debug("Q(" + w + ") = " + Arrays.toString(result));
         }
-        if (debug) {
-            log.debug("Q(" + w + ") = " + q);
-        }
-        return q;
+
+        return result;
     }
 }
 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :

http://dive4elements.wald.intevation.org