# HG changeset patch
# User Sascha L. Teichmann <sascha.teichmann@intevation.de>
# Date 1330352190 0
# Node ID 0dd58ab7e1183dfce138593d6ac39bd191f93989
# Parent  53e8bffbe06cafd730b936a258f78790c4ebe3c1
Added functions to be used for fitting in the "Fixierungsanalyse" and "Extremwertermittlung".

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

diff -r 53e8bffbe06c -r 0dd58ab7e118 flys-artifacts/ChangeLog
--- a/flys-artifacts/ChangeLog	Mon Feb 27 10:50:14 2012 +0000
+++ b/flys-artifacts/ChangeLog	Mon Feb 27 14:16:30 2012 +0000
@@ -1,3 +1,58 @@
+2012-02-27  Sascha L. Teichmann	<sascha.teichmann@intevation.de>
+
+	Added functions to be used for fitting in 
+	the "Fixierungsanalyse" and "Extremwertermittlung".
+
+	* src/main/java/de/intevation/flys/artifacts/math/fitting/Function.java: New.
+	  Abstract base class for functions to fit. Provides the name of the function,
+	  a short description, the names of the parameters and an initial parameter guess
+	  for the fit processe. Sub classes have to overwite the function evaluation and
+	  the partial derivative of the function in respect to the parameters.
+
+	  TODO: Add a meachnism for the inverse function (needed for AT export) and
+	  the first derivative (needed for the respective diagram).
+
+	* src/main/java/de/intevation/flys/artifacts/math/fitting/FunctionFactory.java:
+	  New. Factory to fetch a function by its name.
+
+	 * src/main/java/de/intevation/flys/artifacts/math/fitting/Exp.java: New.
+	   exp: W(Q) = m * a^Q + b
+
+	 * src/main/java/de/intevation/flys/artifacts/math/fitting/Quad.java: New.
+	   quad: W(Q) = n*Q^2 + m*Q + b
+
+	 * src/main/java/de/intevation/flys/artifacts/math/fitting/Linear.java: New.
+	   linear: W(Q) = m*Q + b
+
+	 * src/main/java/de/intevation/flys/artifacts/math/fitting/LogLinear.java: New.
+	   log-linear: W(Q) = a*ln(m*Q + b)
+
+	 * src/main/java/de/intevation/flys/artifacts/math/fitting/Log.java: New
+	   log: W(Q) = m*ln(Q + b)
+
+	 * src/main/java/de/intevation/flys/artifacts/math/fitting/Pow.java: New.
+	   pow: W(Q) = a*Q^c + d
+
+	   !!! This power function is new in the pool of functions to be fit. !!!
+	   See my mail "Manuelle Punkte in der Fixierungsanalyse" 2011-10-27 for details.
+	   The function exp-new found in the old FLYS function pool is omitted 
+	   because it is worthless and was maybe never used.
+
+	* src/main/java/de/intevation/flys/artifacts/math/fitting/App.java: New.
+	  Small test driver to check if the fitting is working. The points to
+	  fit are read from stdin the function to fit is determined by the
+	  system property 'function'. Example usage:
+
+	    $ mvn -e \
+	    -Dfunction=linear \
+	    -Dexec.mainClass=de.intevation.flys.artifacts.math.fitting.App exec:java <<EOF
+	    357.390696917 7546.72096163
+	    61.4291036312 1334.54835721
+	    799.962128234 16836.7698076
+	    126.52761023 2703.69789985
+	    900.448553398 18955.0578748
+	    EOF
+
 2012-02-27	Felix Wolfsteller	<felix.wolfsteller@intevation.de>
 
 	* src/main/resources/messages_de.properties:
diff -r 53e8bffbe06c -r 0dd58ab7e118 flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/fitting/App.java
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/fitting/App.java	Mon Feb 27 14:16:30 2012 +0000
@@ -0,0 +1,127 @@
+package de.intevation.flys.artifacts.math.fitting;
+
+import java.util.List;
+import java.util.Map;
+import java.util.ArrayList;
+import java.util.TreeMap;
+import java.util.Comparator;
+
+import java.io.IOException;
+import java.io.BufferedReader;
+import java.io.Reader;
+import java.io.InputStreamReader;
+
+import org.apache.commons.math.optimization.fitting.CurveFitter;
+
+import org.apache.commons.math.optimization.general.LevenbergMarquardtOptimizer;
+
+import org.apache.commons.math.MathException;
+
+public class App
+{
+    public static final double EPS = 1e-5;
+
+    public static final String FUNCTION_NAME =
+        System.getProperty("function", "linear");
+
+    public static final Comparator<Double> EPS_CMP =
+        new Comparator<Double>()  {
+            @Override
+            public int compare(Double a, Double b) {
+                double diff = a - b;
+                if (diff < -EPS) return -1;
+                if (diff >  EPS) return +1;
+                return 0;
+            }
+        };
+
+    public static final List<Double []>readPoints(Reader reader) 
+    throws IOException
+    {
+        Map<Double, Double> map = new TreeMap<Double, Double>(EPS_CMP);
+
+        BufferedReader input = new BufferedReader(reader);
+
+        String line;
+        while ((line = input.readLine()) != null) {
+            if ((line = line.trim()).length() == 0 || line.startsWith("#")) {
+                continue;
+            }
+
+            String [] parts = line.split("\\s+");
+
+            if (parts.length < 2) {
+                continue;
+            }
+
+            try {
+                Double x = Double.valueOf(parts[0]);
+                Double y = Double.valueOf(parts[1]);
+
+                Double old = map.put(x, y);
+
+                if (old != null) {
+                    System.err.println("duplicate x: " + x);
+                }
+            }
+            catch (NumberFormatException nfe) {
+                nfe.printStackTrace();
+            }
+        }
+
+        List<Double []> list = new ArrayList<Double []>(map.size());
+
+        for (Map.Entry<Double, Double> entry: map.entrySet()) {
+            list.add(new Double [] { entry.getKey(), entry.getValue() });
+        }
+
+        return list;
+    }
+
+    public static void main(String [] args) {
+
+        Function function = FunctionFactory
+            .getInstance()
+            .getFunction(FUNCTION_NAME);
+
+        if (function == null) {
+            System.err.println("Cannot find function '" + FUNCTION_NAME + "'.");
+            System.exit(1);
+        }
+
+        List<Double []> points = null;
+
+        try {
+            points = readPoints(new InputStreamReader(System.in));
+        }
+        catch (IOException ioe) {
+            ioe.printStackTrace();
+            System.exit(1);
+        }
+
+        LevenbergMarquardtOptimizer lmo = new LevenbergMarquardtOptimizer();
+
+        CurveFitter cf = new CurveFitter(lmo);
+
+        for (Double [] point: points) {
+            cf.addObservedPoint(point[0], point[1]);
+        }
+
+        double [] parameters = null;
+
+        try {
+            parameters = cf.fit(function, function.getInitialGuess());
+        }
+        catch (MathException me) {
+            me.printStackTrace();
+            System.exit(1);
+        }
+
+        String [] parameterNames = function.getParameterNames();
+
+        for (int i = 0; i < parameterNames.length; ++i) {
+            System.err.println(parameterNames[i] + ": " + parameters[i]);
+        }
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
diff -r 53e8bffbe06c -r 0dd58ab7e118 flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/fitting/Exp.java
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/fitting/Exp.java	Mon Feb 27 14:16:30 2012 +0000
@@ -0,0 +1,27 @@
+package de.intevation.flys.artifacts.math.fitting;
+
+public class Exp
+extends      Function
+{
+    public Exp() {
+        super(
+            "exp",
+            "W(Q) = m * a^Q + b",
+            new String [] { "m", "a", "b" });
+    }
+
+    @Override
+    public double value(double x, double [] parameters) {
+        return parameters[0]*Math.pow(parameters[1], x) + parameters[2];
+    }
+
+    @Override
+    public double [] gradient(double x, double [] parameters) {
+        return new double [] {
+            Math.pow(parameters[1], x),
+            Math.pow(parameters[1], x-1d)*x,
+            1d
+        };
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
diff -r 53e8bffbe06c -r 0dd58ab7e118 flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/fitting/Function.java
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/fitting/Function.java	Mon Feb 27 14:16:30 2012 +0000
@@ -0,0 +1,56 @@
+package de.intevation.flys.artifacts.math.fitting;
+
+import org.apache.commons.math.optimization.fitting.ParametricRealFunction;
+
+import de.intevation.flys.utils.DoubleUtil;
+
+public abstract class Function
+implements            ParametricRealFunction
+{
+    protected String    name;
+    protected String    description;
+    protected String [] parameterNames;
+    protected double [] initialGuess;
+
+    public Function() {
+    }
+
+    public Function(
+        String    name,
+        String    description,
+        String [] parameterNames
+    ) {
+        this(name,
+            description,
+            parameterNames,
+            DoubleUtil.fill(parameterNames.length, 1d));
+    }
+
+    public Function(
+        String    name,
+        String    description,
+        String [] parameterNames,
+        double [] initialGuess
+    ) {
+        this.name           = name;
+        this.parameterNames = parameterNames;
+        this.initialGuess   = initialGuess;
+    }
+
+    public String getName() {
+        return name;
+    }
+
+    public String getDescription() {
+        return description;
+    }
+
+    public String [] getParameterNames() {
+        return parameterNames;
+    }
+
+    public double [] getInitialGuess() {
+        return initialGuess;
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
diff -r 53e8bffbe06c -r 0dd58ab7e118 flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/fitting/FunctionFactory.java
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/fitting/FunctionFactory.java	Mon Feb 27 14:16:30 2012 +0000
@@ -0,0 +1,43 @@
+package de.intevation.flys.artifacts.math.fitting;
+
+import java.util.Collection;
+import java.util.LinkedHashMap;
+import java.util.Map;
+
+public final class FunctionFactory
+{
+    private static FunctionFactory instance;
+
+    private Map<String, Function> functions;
+
+    private FunctionFactory() {
+        functions = new LinkedHashMap<String, Function>();
+
+        registerFunction(new Log());
+        registerFunction(new Linear());
+        registerFunction(new LogLinear());
+        registerFunction(new Exp());
+        registerFunction(new Quad());
+        registerFunction(new Pow());
+    }
+
+    public static synchronized FunctionFactory getInstance() {
+        if (instance == null) {
+            instance = new FunctionFactory();
+        }
+        return instance;
+    }
+
+    public Function getFunction(String name) {
+        return functions.get(name);
+    }
+
+    public void registerFunction(Function function) {
+        functions.put(function.getName(), function);
+    }
+
+    public Collection<Function> getFunctions() {
+        return functions.values();
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
diff -r 53e8bffbe06c -r 0dd58ab7e118 flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/fitting/Linear.java
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/fitting/Linear.java	Mon Feb 27 14:16:30 2012 +0000
@@ -0,0 +1,20 @@
+package de.intevation.flys.artifacts.math.fitting;
+
+public class Linear
+extends      Function
+{
+    public Linear() {
+        super("linear", "W(Q) = m*Q + b", new String [] { "m", "b" });
+    }
+
+    @Override
+    public double value(double x, double [] parameters) {
+        return x*parameters[0] + parameters[1];
+    }
+
+    @Override
+    public double [] gradient(double x, double [] parameters) {
+        return new double [] { x, 1d };
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
diff -r 53e8bffbe06c -r 0dd58ab7e118 flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/fitting/Log.java
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/fitting/Log.java	Mon Feb 27 14:16:30 2012 +0000
@@ -0,0 +1,23 @@
+package de.intevation.flys.artifacts.math.fitting;
+
+public class Log
+extends      Function
+{
+    public Log() {
+        super("log",  "W(Q) = m*ln(Q + b)", new String [] { "m", "b" });
+    }
+
+    @Override
+    public double value(double x, double [] parameters) {
+        return parameters[0]*Math.log(x + parameters[1]);
+    }
+
+    @Override
+    public double [] gradient(double x, double [] parameters) {
+        return new double [] {
+            Math.log(x + parameters[1]),
+            parameters[0]/(x + parameters[0])
+        };
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
diff -r 53e8bffbe06c -r 0dd58ab7e118 flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/fitting/LogLinear.java
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/fitting/LogLinear.java	Mon Feb 27 14:16:30 2012 +0000
@@ -0,0 +1,28 @@
+package de.intevation.flys.artifacts.math.fitting;
+
+public class LogLinear
+extends      Function
+{
+    public LogLinear() {
+        super(
+            "log-linear",
+            "W(Q) = a*ln(m*Q + b)",
+            new String [] { "a", "m", "b" });
+    }
+
+    @Override
+    public double value(double x, double [] parameters) {
+        return parameters[0]*Math.log(parameters[1]*x + parameters[2]);
+    }
+
+    @Override
+    public double [] gradient(double x, double [] parameters) {
+        double l = parameters[1]*x + parameters[2];
+        return new double [] {
+            Math.log(l),
+            parameters[0]*x/l,
+            parameters[0]/l
+        };
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
diff -r 53e8bffbe06c -r 0dd58ab7e118 flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/fitting/Pow.java
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/fitting/Pow.java	Mon Feb 27 14:16:30 2012 +0000
@@ -0,0 +1,30 @@
+package de.intevation.flys.artifacts.math.fitting;
+
+public class Pow
+extends      Function
+{
+    public Pow() {
+        super(
+            "pow",
+            "W(Q) = a*Q^c + d",
+            new String [] { "a", "c", "d" });
+    }
+
+    @Override
+    public double value(double x, double [] parameters) {
+        return parameters[0]*Math.pow(x, parameters[1]) + parameters[2];
+    }
+
+    @Override
+    public double [] gradient(double x, double [] parameters) {
+        double a   = parameters[0];
+        double c   = parameters[1];
+        double x_c = Math.pow(x, c);
+        return new double [] {
+            x_c,
+            a*x_c*Math.log(x),
+            1d
+        };
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
diff -r 53e8bffbe06c -r 0dd58ab7e118 flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/fitting/Quad.java
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/fitting/Quad.java	Mon Feb 27 14:16:30 2012 +0000
@@ -0,0 +1,24 @@
+package de.intevation.flys.artifacts.math.fitting;
+
+public class Quad
+extends      Function
+{
+    public Quad() {
+        super(
+            "quad",
+            "W(Q) = n*Q^2 + m*Q + b",
+            new String [] { "n", "m", "b" });
+    }
+
+    @Override
+    public double value(double x, double [] parameters) {
+        // n*Q^2 + m*Q + b <=> Q*(n*Q + m) + b
+        return x*(parameters[0]*x + parameters[1]) + parameters[2];
+    }
+
+    @Override
+    public double [] gradient(double x, double [] parameters) {
+        return new double [] { x*x, x, 1d };
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
diff -r 53e8bffbe06c -r 0dd58ab7e118 flys-artifacts/src/main/java/de/intevation/flys/utils/DoubleUtil.java
--- a/flys-artifacts/src/main/java/de/intevation/flys/utils/DoubleUtil.java	Mon Feb 27 10:50:14 2012 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/utils/DoubleUtil.java	Mon Feb 27 14:16:30 2012 +0000
@@ -132,5 +132,11 @@
 
         return out;
     }
+
+    public static final double [] fill(int N, double value) {
+        double [] result = new double[N];
+        Arrays.fill(result, value);
+        return result;
+    }
 }
 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :