Mercurial > dive4elements > river

changeset 8581:073ea4bcea58

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
(issue1755) Interpolate BedQuality Results This adds an interpolation function to each various bedQuality result class. Imho this is ok as the interpolation function can be seen as part of the result. The interpolation function is initalized on first use and can be accessed through get.*Interpol functions.
author Andre Heinecke <andre.heinecke@intevation.de>
date Mon, 16 Mar 2015 15:36:38 +0100
parents d9f038b8e2ce
children 17a3030bbda2
files artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/BedDiameterResult.java artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/BedParametersResult.java artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/BedloadDiameterResult.java artifacts/src/main/java/org/dive4elements/river/exports/minfo/BedQualityExporter.java artifacts/src/main/java/org/dive4elements/river/utils/DoubleUtil.java
diffstat 5 files changed, 225 insertions(+), 44 deletions(-) [+]
line wrap: on
line diff
--- a/artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/BedDiameterResult.java	Mon Mar 16 14:37:30 2015 +0100
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/BedDiameterResult.java	Mon Mar 16 15:36:38 2015 +0100
@@ -8,8 +8,13 @@
 
 package org.dive4elements.river.artifacts.model.minfo;
 
+import org.dive4elements.river.utils.DoubleUtil;
+
 import gnu.trove.TDoubleArrayList;
 
+import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
+
+import org.apache.commons.math.ArgumentOutsideDomainException;
 
 public class BedDiameterResult
 extends BedQualityDiameterResult
@@ -17,6 +22,11 @@
     protected TDoubleArrayList diameterCap;
     protected TDoubleArrayList diameterSub;
 
+    protected PolynomialSplineFunction interpolSub;
+    protected PolynomialSplineFunction interpolCap;
+    protected boolean nonInterpolSub;
+    protected boolean nonInterpolCap;
+
     public BedDiameterResult (
         String type,
         TDoubleArrayList diameterCap,
@@ -26,6 +36,10 @@
         super(type, km);
         this.diameterCap = diameterCap;
         this.diameterSub = diameterSub;
+        interpolSub = null;
+        nonInterpolSub = false;
+        interpolCap = null;
+        nonInterpolCap = false;
     }
 
     public double getDiameterCap(int ndx) {
@@ -49,6 +63,25 @@
         return Double.NaN;
     }
 
+    public double getDiameterCapInterpol(double km) {
+        if (nonInterpolCap) {
+            return Double.NaN;
+        }
+        if (interpolCap == null) {
+            interpolCap = DoubleUtil.getLinearInterpolator(kms, diameterCap);
+            if (interpolCap == null) {
+                nonInterpolCap = true;
+                return Double.NaN;
+            }
+        }
+        try {
+            return interpolCap.value(km);
+        } catch (ArgumentOutsideDomainException e) {
+            /* This is expected for many results. */
+            return Double.NaN;
+        }
+    }
+
     public double getDiameterSub(double km) {
         if (kms.indexOf(km) >= 0) {
             return diameterSub.get(kms.indexOf(km));
@@ -56,6 +89,25 @@
         return Double.NaN;
     }
 
+    public double getDiameterSubInterpol(double km) {
+        if (nonInterpolSub) {
+            return Double.NaN;
+        }
+        if (interpolSub == null) {
+            interpolSub = DoubleUtil.getLinearInterpolator(kms, diameterSub);
+            if (interpolSub == null) {
+                nonInterpolSub = true;
+                return Double.NaN;
+            }
+        }
+        try {
+            return interpolSub.value(km);
+        } catch (ArgumentOutsideDomainException e) {
+            /* This is expected for many results. */
+            return Double.NaN;
+        }
+    }
+
     public double[][] getDiameterCapData() {
         return new double[][] {
             kms.toNativeArray(),
--- a/artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/BedParametersResult.java	Mon Mar 16 14:37:30 2015 +0100
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/BedParametersResult.java	Mon Mar 16 15:36:38 2015 +0100
@@ -9,6 +9,9 @@
 package org.dive4elements.river.artifacts.model.minfo;
 
 import gnu.trove.TDoubleArrayList;
+import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
+import org.dive4elements.river.utils.DoubleUtil;
+import org.apache.commons.math.ArgumentOutsideDomainException;
 
 import java.io.Serializable;
 
@@ -22,6 +25,16 @@
     protected TDoubleArrayList loadDensitySub;
     protected TDoubleArrayList kms;
 
+    protected PolynomialSplineFunction interpolPoroSub;
+    protected PolynomialSplineFunction interpolPoroCap;
+    protected PolynomialSplineFunction interpolDensSub;
+    protected PolynomialSplineFunction interpolDensCap;
+
+    protected boolean nonInterpolPoroSub;
+    protected boolean nonInterpolPoroCap;
+    protected boolean nonInterpolDensSub;
+    protected boolean nonInterpolDensCap;
+
     public BedParametersResult() {
 
     }
@@ -38,6 +51,16 @@
         this.porositySub = porositySub;
         this.loadDensityCap = densityCap;
         this.loadDensitySub = densitySub;
+
+        PolynomialSplineFunction interpolPoroSub = null;
+        PolynomialSplineFunction interpolPoroCap = null;
+        PolynomialSplineFunction interpolDensSub = null;
+        PolynomialSplineFunction interpolDensCap = null;
+
+        nonInterpolPoroSub = false;
+        nonInterpolPoroCap = false;
+        nonInterpolDensSub = false;
+        nonInterpolDensCap = false;
     }
 
     public double getPorosityCap(int ndx) {
@@ -111,4 +134,81 @@
             loadDensitySub.toNativeArray()
         };
     }
+
+    public double getPorositySubInterpol(double km) {
+        if (nonInterpolPoroSub) {
+            return Double.NaN;
+        }
+        if (interpolPoroSub == null) {
+            interpolPoroSub = DoubleUtil.getLinearInterpolator(kms, porositySub);
+            if (interpolPoroSub == null) {
+                nonInterpolPoroSub = true;
+                return Double.NaN;
+            }
+        }
+        try {
+            return interpolPoroSub.value(km);
+        } catch (ArgumentOutsideDomainException e) {
+            /* This is expected for many results. */
+            return Double.NaN;
+        }
+    }
+
+    public double getPorosityCapInterpol(double km) {
+        if (nonInterpolPoroCap) {
+            return Double.NaN;
+        }
+        if (interpolPoroCap == null) {
+            interpolPoroCap = DoubleUtil.getLinearInterpolator(kms, porosityCap);
+            if (interpolPoroCap == null) {
+                nonInterpolPoroCap = true;
+                return Double.NaN;
+            }
+        }
+        try {
+            return interpolPoroCap.value(km);
+        } catch (ArgumentOutsideDomainException e) {
+            /* This is expected for many results. */
+            return Double.NaN;
+        }
+    }
+
+    public double getDensitySubInterpol(double km) {
+        if (nonInterpolDensSub) {
+            return Double.NaN;
+        }
+        if (interpolDensSub == null) {
+            interpolDensSub = DoubleUtil.getLinearInterpolator(kms, loadDensitySub);
+            if (interpolDensSub == null) {
+                nonInterpolDensSub = true;
+                return Double.NaN;
+            }
+        }
+        try {
+            return interpolDensSub.value(km);
+        } catch (ArgumentOutsideDomainException e) {
+            /* This is expected for many results. */
+            return Double.NaN;
+        }
+    }
+
+    public double getDensityCapInterpol(double km) {
+        if (nonInterpolDensCap) {
+            return Double.NaN;
+        }
+        if (interpolDensCap == null) {
+            interpolDensCap = DoubleUtil.getLinearInterpolator(kms, loadDensityCap);
+            if (interpolDensCap == null) {
+                nonInterpolDensCap = true;
+                return Double.NaN;
+            }
+        }
+        try {
+            return interpolDensCap.value(km);
+        } catch (ArgumentOutsideDomainException e) {
+            /* This is expected for many results. */
+            return Double.NaN;
+        }
+    }
+
 }
--- a/artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/BedloadDiameterResult.java	Mon Mar 16 14:37:30 2015 +0100
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/BedloadDiameterResult.java	Mon Mar 16 15:36:38 2015 +0100
@@ -9,13 +9,22 @@
 package org.dive4elements.river.artifacts.model.minfo;
 
 import org.dive4elements.river.artifacts.model.DateRange;
+import org.dive4elements.river.utils.DoubleUtil;
+
 import gnu.trove.TDoubleArrayList;
 
+import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
+
+import org.apache.commons.math.ArgumentOutsideDomainException;
 
 public class BedloadDiameterResult
 extends BedQualityDiameterResult
 {
     protected TDoubleArrayList diameter;
+    protected PolynomialSplineFunction interpol;
+
+    /** Set to true if this result can't be interpolated.*/
+    protected boolean nonInterpolResult;
 
     public BedloadDiameterResult(
         String type,
@@ -25,6 +34,8 @@
     ) {
         super (type, km);
         this.diameter = diameter;
+        interpol = null;
+        nonInterpolResult = false;
     }
 
     public double getDiameter(int ndx) {
@@ -41,6 +52,25 @@
         return Double.NaN;
     }
 
+    public double getDiameterInterpol(double km) {
+        if (nonInterpolResult) {
+            return Double.NaN;
+        }
+        if (interpol == null) {
+            interpol = DoubleUtil.getLinearInterpolator(kms, diameter);
+            if (interpol == null) {
+                nonInterpolResult = true;
+                return Double.NaN;
+            }
+        }
+        try {
+            return interpol.value(km);
+        } catch (ArgumentOutsideDomainException e) {
+            /* This is expected for many results. */
+            return Double.NaN;
+        }
+    }
+
     public double[][] getDiameterData() {
         return new double[][] {
             kms.toNativeArray(),
--- a/artifacts/src/main/java/org/dive4elements/river/exports/minfo/BedQualityExporter.java	Mon Mar 16 14:37:30 2015 +0100
+++ b/artifacts/src/main/java/org/dive4elements/river/exports/minfo/BedQualityExporter.java	Mon Mar 16 15:36:38 2015 +0100
@@ -1,4 +1,4 @@
-/* Copyright (C) 2011, 2012, 2013 by Bundesanstalt für Gewässerkunde
+/* Copyright (C) 2011, 2012, 2013, 2015 by Bundesanstalt für Gewässerkunde
  * Software engineering by Intevation GmbH
  *
  * This file is Free Software under the GNU AGPL (>=v3)
@@ -29,6 +29,9 @@
 import org.dive4elements.river.exports.AbstractExporter;
 import org.dive4elements.river.utils.Formatter;
 
+import org.dive4elements.river.artifacts.access.RangeAccess;
+import org.dive4elements.river.artifacts.D4EArtifact;
+
 
 public class BedQualityExporter
 extends AbstractExporter
@@ -59,43 +62,16 @@
         results = new BedQualityResult[0];
     }
 
-    /** Populate param kms with kms from param beds and loads, return it. */
-    private TDoubleArrayList populateKmList(BedDiameterResult[] beds,
-        BedloadDiameterResult[] loads,
-        TDoubleArrayList kms) {
-        for (int j = 0; j < beds.length; j++) {
-            TDoubleArrayList bkms = beds[j].getKms();
-            for (int k = 0, K = bkms.size(); k < K; k++) {
-                double km = bkms.get(k);
-                if (!kms.contains(km)) { // XXX: O(N^2)
-                    kms.add(km);
-                }
-            }
-        }
-        for (int j = 0; j < loads.length; j++) {
-            TDoubleArrayList lkms = loads[j].getKms();
-            for (int k = 0, L = lkms.size(); k < L; k++) {
-                double km = lkms.get(k);
-                if (!kms.contains(km)) { // XXX: O(N^2)
-                    kms.add(km);
-                }
-            }
-        }
-        return kms;
-    }
-
     /** Create double[] containing the data for rows in csv. */
     private List<double[]> createDataRows() {
-        // Calculate how many columns and rows we need.
-        TDoubleArrayList kms = new TDoubleArrayList();
+
+        double[] kms = new RangeAccess((D4EArtifact) master).getKmSteps();
 
         int cols = 1;
         for (BedQualityResult result : results) {
             BedDiameterResult[]     beds  = result.getBedResults();
             BedloadDiameterResult[] loads = result.getBedloadResults();
 
-            kms = populateKmList(beds, loads, kms);
-
             cols += beds.length * 2;
             if (beds.length > 0) {
                 cols += 4;
@@ -103,22 +79,21 @@
             cols += loads.length;
         }
 
-        kms.sort();
-
-        List<double[]> rows = new ArrayList<double[]>(kms.size());
-        for (int i = 0, K = kms.size(); i < K; i++) {
+        List<double[]> rows = new ArrayList<double[]>(kms.length);
+        for (double km: kms) {
             double[] row = new double[cols];
-            double km = kms.get(i);
             row[0] = km;
             for (int j = 0; j < results.length; j++) {
+
                 BedloadDiameterResult[] loads = results[j].getBedloadResults();
 
                 for(int k = 0; k < loads.length; k++) {
                     // k + 1: shift km column.
                     // j* loads.length: shift periods.
                     row[(k + 1) + (j * loads.length)] =
-                        loads[k].getDiameter(km);
+                        loads[k].getDiameterInterpol(km);
                 }
+
                 BedDiameterResult[] beds = results[j].getBedResults();
                 if (beds.length == 0) {
                     continue;
@@ -128,17 +103,18 @@
                     // j * beds.length * 2: shift periods.
                     // loads.length * results.length: shift bed load columns.
                     int ndx = (k * 2 + 1) + (j * beds.length * 2) + (loads.length * results.length);
-                    row[ndx] = beds[k].getDiameterCap(km);
-                    row[ndx + 1] = beds[k].getDiameterSub(km);
+                    row[ndx] = beds[k].getDiameterCapInterpol(km);
+                    row[ndx + 1] = beds[k].getDiameterSubInterpol(km);
                 }
+
                 BedParametersResult[] params = results[j].getParameters();
                 for(int k = 0; k < params.length; k++) {
                     // loads.length + (beds.lenght * 2) * (j + 1): shift bed and bedload columns.
                     int ndx = 1 + (loads.length + (beds.length * 2) * (j + 1));
-                    row[ndx] = params[k].getLoadDensityCap(km);
-                    row[ndx + 1] = params[k].getLoadDensitySub(km);
-                    row[ndx + 2] = params[k].getPorosityCap(km);
-                    row[ndx + 3] = params[k].getPorositySub(km);
+                    row[ndx] = params[k].getDensityCapInterpol(km);
+                    row[ndx + 1] = params[k].getDensitySubInterpol(km);
+                    row[ndx + 2] = params[k].getPorosityCapInterpol(km);
+                    row[ndx + 3] = params[k].getPorositySubInterpol(km);
                 }
             }
             rows.add(row);
@@ -149,7 +125,6 @@
 
     @Override
     protected void writeCSVData(CSVWriter writer) throws IOException {
-        // TODO Auto-generated method stub
         writeCSVHeader(writer);
 
         NumberFormat nf = Formatter.getFormatter(context, 1, 3);
@@ -177,7 +152,6 @@
 
     @Override
     protected void addData(Object data) {
-        // TODO Auto-generated method stub
         log.debug("addData()");
         if (!(data instanceof CalculationResult)) {
             log.warn("Invalid data type.");
--- a/artifacts/src/main/java/org/dive4elements/river/utils/DoubleUtil.java	Mon Mar 16 14:37:30 2015 +0100
+++ b/artifacts/src/main/java/org/dive4elements/river/utils/DoubleUtil.java	Mon Mar 16 15:36:38 2015 +0100
@@ -10,6 +10,13 @@
 
 import org.dive4elements.river.artifacts.math.Linear;
 
+import org.apache.commons.math.analysis.interpolation.LinearInterpolator;
+import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
+
+import org.apache.commons.math.exception.DimensionMismatchException;
+import org.apache.commons.math.exception.NumberIsTooSmallException;
+import org.apache.commons.math.exception.NonMonotonousSequenceException;
+
 import gnu.trove.TDoubleArrayList;
 
 import java.util.Arrays;
@@ -305,5 +312,23 @@
             }
         }
     }
+
+    /** Convieniance function for results to get an interpolator.
+     * This is basically a static wrapper to for LinearInterpolator.interpolate
+     * with error handling. Returns null on error.*/
+    public static PolynomialSplineFunction getLinearInterpolator(TDoubleArrayList x, TDoubleArrayList y) {
+        LinearInterpolator lpol = new LinearInterpolator();
+        try {
+            return lpol.interpolate(x.toNativeArray(), y.toNativeArray());
+        } catch (DimensionMismatchException e) {
+            log.error("KMs and Result values have different sizes. Failed to interpolate: " +
+                    e.getMessage());
+        } catch (NonMonotonousSequenceException e) {
+            log.error("KMs are not monotonous. Failed to interpolate: " + e.getMessage());
+        } catch (NumberIsTooSmallException e) {
+            log.error("Result is to small. Failed to interpolate: " + e.getMessage());
+        }
+        return null;
+    }
 }
 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
http://dive4elements.wald.intevation.org