Mercurial > dive4elements > river
changeset 2569:0dd58ab7e118
Added functions to be used for fitting in the "Fixierungsanalyse" and "Extremwertermittlung".
flys-artifacts/trunk@4095 c6561f87-3c4e-4783-a992-168aeb5c3f6f
line wrap: on
line diff
--- 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:
--- /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 :
--- /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 :
--- /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 :
--- /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 :
--- /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 :
--- /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 :
--- /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 :
--- /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 :
--- /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 :
--- 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 :