[PATCH] Implemented alternative fitting strategy for Log-Linear function
Wald Commits
scm-commit at wald.intevation.org
Mon Dec 2 17:56:30 CET 2019
# HG changeset patch
# User Gernot Belger <g.belger at bjoernsen.de>
# Date 1575305775 -3600
# Mon Dec 02 17:56:15 2019 +0100
# Node ID 0380717105bad5b29f62228402858948db6723b0
# Parent eb1a29fe823f77d05541a50d0b4ffa43ddfb0034
Implemented alternative fitting strategy for Log-Linear function.
diff -r eb1a29fe823f -r 0380717105ba artifacts/pom.xml
--- a/artifacts/pom.xml Mon Dec 02 17:21:43 2019 +0100
+++ b/artifacts/pom.xml Mon Dec 02 17:56:15 2019 +0100
@@ -133,7 +133,7 @@
<groupId>org.apache.commons</groupId>
<artifactId>commons-math</artifactId>
<version>2.2</version>
- </dependency>
+ </dependency>
<dependency>
<groupId>com.h2database</groupId>
<artifactId>h2</artifactId>
@@ -195,6 +195,16 @@
<artifactId>groovy-all</artifactId>
<version>1.6.0</version>
</dependency>
+ <dependency>
+ <groupId>org.apache.commons</groupId>
+ <artifactId>commons-math3</artifactId>
+ <version>3.6.1</version>
+ </dependency>
+ <dependency>
+ <groupId>org.hipparchus</groupId>
+ <artifactId>hipparchus-optim</artifactId>
+ <version>1.6</version>
+ </dependency>
</dependencies>
<repositories>
<repository>
diff -r eb1a29fe823f -r 0380717105ba artifacts/src/main/java/org/dive4elements/river/artifacts/math/fitting/LogLinearAlternative.java
--- a/artifacts/src/main/java/org/dive4elements/river/artifacts/math/fitting/LogLinearAlternative.java Mon Dec 02 17:21:43 2019 +0100
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/math/fitting/LogLinearAlternative.java Mon Dec 02 17:56:15 2019 +0100
@@ -15,7 +15,7 @@
*/
public class LogLinearAlternative extends AbstractLogLinear {
- public static final String NAME = "log-linear-alternative";
+ public static final String NAME = "log-linear-linearisiert";
public static final Function INSTANCE = new LogLinearAlternative();
public LogLinearAlternative() {
diff -r eb1a29fe823f -r 0380717105ba artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/AbstractFittingOperation.java
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/AbstractFittingOperation.java Mon Dec 02 17:56:15 2019 +0100
@@ -0,0 +1,86 @@
+/** Copyright (C) 2017 by Bundesanstalt für Gewässerkunde
+ * Software engineering by
+ * Björnsen Beratende Ingenieure GmbH
+ * Dr. Schumacher Ingenieurbüro für Wasser und Umwelt
+ *
+ * This file is Free Software under the GNU AGPL (>=v3)
+ * and comes with ABSOLUTELY NO WARRANTY! Check out the
+ * documentation coming with Dive4Elements River for details.
+ */
+package org.dive4elements.river.artifacts.model.fixings;
+
+import java.util.List;
+
+import org.apache.commons.math3.stat.descriptive.moment.StandardDeviation;
+import org.dive4elements.river.artifacts.math.fitting.Function;
+
+/**
+ * @author Gernot Belger
+ */
+abstract class AbstractFittingOperation {
+
+ private final Function function;
+
+ public AbstractFittingOperation(final Function function) {
+ this.function = function;
+ }
+
+ protected final Function getFunction() {
+ return this.function;
+ }
+
+ protected final Fitting createFitting(final double[] parameters, final double chiSquare, final List<FittingData> data) {
+
+ final org.dive4elements.river.artifacts.math.Function instance = this.function.instantiate(parameters);
+ final double[] waterlevels = calculateWaterlevels(instance, data);
+ final double maxQ = calculateMaxQ(data);
+ final double stdDev = calculateStandardDeviation(waterlevels, data);
+
+ return new Fitting(parameters, stdDev, chiSquare, maxQ);
+ }
+
+ private double calculateMaxQ(final List<FittingData> data) {
+
+ double maxQ = -Double.MAX_VALUE;
+
+ for (int i = 0; i < data.size(); i++) {
+
+ final FittingData fittingData = data.get(i);
+
+ if (fittingData.q > maxQ)
+ maxQ = fittingData.q;
+ }
+
+ return maxQ;
+ }
+
+ private double calculateStandardDeviation(final double[] waterlevels, final List<FittingData> data) {
+
+ final StandardDeviation std = new StandardDeviation();
+
+ for (int i = 0; i < data.size(); i++) {
+
+ final FittingData fittingData = data.get(i);
+
+ final double diff = (waterlevels[i] - fittingData.q) * 100.0; // * 100 for whatever reason
+
+ std.increment(diff);
+ }
+
+ return std.getResult();
+ }
+
+ private double[] calculateWaterlevels(final org.dive4elements.river.artifacts.math.Function instance, final List<FittingData> data) {
+
+ final double[] waterlevels = new double[data.size()];
+
+ for (int i = 0; i < data.size(); i++) {
+
+ final FittingData fittingData = data.get(i);
+
+ waterlevels[i] = instance.value(fittingData.q);
+ }
+
+ return waterlevels;
+ }
+}
\ No newline at end of file
diff -r eb1a29fe823f -r 0380717105ba artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/Fitting.java
--- a/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/Fitting.java Mon Dec 02 17:21:43 2019 +0100
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/Fitting.java Mon Dec 02 17:56:15 2019 +0100
@@ -12,13 +12,10 @@
import java.util.Date;
import java.util.List;
-import org.apache.commons.math.MathException;
-import org.apache.commons.math.optimization.fitting.CurveFitter;
-import org.apache.commons.math.optimization.general.LevenbergMarquardtOptimizer;
-import org.apache.commons.math.stat.descriptive.moment.StandardDeviation;
import org.apache.log4j.Logger;
import org.dive4elements.river.artifacts.math.GrubbsOutlier;
import org.dive4elements.river.artifacts.math.fitting.Function;
+import org.dive4elements.river.artifacts.math.fitting.LogLinearAlternative;
public class Fitting {
private static Logger log = Logger.getLogger(Fitting.class);
@@ -28,20 +25,6 @@
QWD create(double q, double w, double deltaW, boolean isOutlier);
}
- private static final class FittingData {
- public final FixingColumnWithData event;
- public final double w;
- public final double q;
- public final boolean isInterpolated;
-
- public FittingData(final FixingColumnWithData event, final double w, final double q, final boolean isInterpolated) {
- this.event = event;
- this.w = w;
- this.q = q;
- this.isInterpolated = isInterpolated;
- }
- }
-
private final double chiSqr;
private final double[] parameters;
@@ -101,37 +84,14 @@
final List<FittingData> outliers = new ArrayList<>(data.size());
org.dive4elements.river.artifacts.math.Function instance = null;
- LevenbergMarquardtOptimizer lmo = null;
- double[] parameters = null;
+ Fitting fitting = null;
+ final IFittingOperation operation = createOperation(function);
for (;;) {
- parameters = null;
-
- for (double tolerance = 1e-10; tolerance < 1e-1; tolerance *= 10d) {
-
- lmo = new LevenbergMarquardtOptimizer();
- lmo.setCostRelativeTolerance(tolerance);
- lmo.setOrthoTolerance(tolerance);
- lmo.setParRelativeTolerance(tolerance);
-
- try {
- final CurveFitter cf = new CurveFitter(lmo);
+ fitting = operation.execute(data);
- for (final FittingData fittingData : data)
- cf.addObservedPoint(fittingData.q, fittingData.w);
-
- parameters = cf.fit(function, function.getInitialGuess());
- break;
- }
- catch (final MathException me) {
- if (log.isDebugEnabled()) {
- log.debug("tolerance " + tolerance + " + failed.", me);
- }
- }
- }
-
- if (parameters == null) {
+ if (fitting == null) {
/*
* log.debug("Parameters is null");
* for (int i = 0, N = xs.size(); i < N; ++i) {
@@ -142,6 +102,7 @@
}
// This is the parameterized function for a given km.
+ final double[] parameters = fitting.getParameters();
instance = function.instantiate(parameters);
if (!checkOutliers)
@@ -180,30 +141,24 @@
}
/*
- * calculate dW of used values against the resulting function and add them to results , also calculate standard *
- * deviation
+ * calculate dW of used values against the resulting function and add them to results
*/
- final StandardDeviation stdDev = new StandardDeviation();
- double maxQ = -Double.MAX_VALUE;
-
for (final FittingData fittingData : data) {
final QWD qwd = createQWD(fittingData, instance, false);
qwds.add(qwd);
resultColumns.addQWD(fittingData.event, km, qwd);
-
- stdDev.increment(qwd.getDeltaW());
-
- final double q = qwd.getQ();
- if (q > maxQ)
- maxQ = q;
}
- final double standardDeviation = stdDev.getResult();
+ return fitting;
+ }
- final double chiSqr = lmo.getChiSquare();
+ private static IFittingOperation createOperation(final Function function) {
- return new Fitting(parameters, standardDeviation, chiSqr, maxQ);
+ if (function instanceof LogLinearAlternative)
+ return new LogLinearFittingOperation((LogLinearAlternative) function);
+
+ return new LevenbergMarquardtFittingOperation(function);
}
private static QWD createQWD(final FittingData data, final org.dive4elements.river.artifacts.math.Function function, final boolean isOutlier) {
diff -r eb1a29fe823f -r 0380717105ba artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/FittingData.java
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/FittingData.java Mon Dec 02 17:56:15 2019 +0100
@@ -0,0 +1,24 @@
+/** Copyright (C) 2017 by Bundesanstalt für Gewässerkunde
+ * Software engineering by
+ * Björnsen Beratende Ingenieure GmbH
+ * Dr. Schumacher Ingenieurbüro für Wasser und Umwelt
+ *
+ * This file is Free Software under the GNU AGPL (>=v3)
+ * and comes with ABSOLUTELY NO WARRANTY! Check out the
+ * documentation coming with Dive4Elements River for details.
+ */
+package org.dive4elements.river.artifacts.model.fixings;
+
+final class FittingData {
+ public final FixingColumnWithData event;
+ public final double w;
+ public final double q;
+ public final boolean isInterpolated;
+
+ public FittingData(final FixingColumnWithData event, final double w, final double q, final boolean isInterpolated) {
+ this.event = event;
+ this.w = w;
+ this.q = q;
+ this.isInterpolated = isInterpolated;
+ }
+}
\ No newline at end of file
diff -r eb1a29fe823f -r 0380717105ba artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/IFittingOperation.java
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/IFittingOperation.java Mon Dec 02 17:56:15 2019 +0100
@@ -0,0 +1,24 @@
+/** Copyright (C) 2017 by Bundesanstalt für Gewässerkunde
+ * Software engineering by
+ * Björnsen Beratende Ingenieure GmbH
+ * Dr. Schumacher Ingenieurbüro für Wasser und Umwelt
+ *
+ * This file is Free Software under the GNU AGPL (>=v3)
+ * and comes with ABSOLUTELY NO WARRANTY! Check out the
+ * documentation coming with Dive4Elements River for details.
+ */
+package org.dive4elements.river.artifacts.model.fixings;
+
+import java.util.List;
+
+/**
+ * Abstraction for different fitting approaches depending on target function.
+ *
+ * @author Gernot Belger
+ *
+ */
+public interface IFittingOperation {
+
+ Fitting execute(List<FittingData> data);
+
+}
\ No newline at end of file
diff -r eb1a29fe823f -r 0380717105ba artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/LevenbergMarquardtFittingOperation.java
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/LevenbergMarquardtFittingOperation.java Mon Dec 02 17:56:15 2019 +0100
@@ -0,0 +1,62 @@
+/** Copyright (C) 2017 by Bundesanstalt für Gewässerkunde
+ * Software engineering by
+ * Björnsen Beratende Ingenieure GmbH
+ * Dr. Schumacher Ingenieurbüro für Wasser und Umwelt
+ *
+ * This file is Free Software under the GNU AGPL (>=v3)
+ * and comes with ABSOLUTELY NO WARRANTY! Check out the
+ * documentation coming with Dive4Elements River for details.
+ */
+package org.dive4elements.river.artifacts.model.fixings;
+
+import java.util.List;
+
+import org.apache.commons.math.MathException;
+import org.apache.commons.math.optimization.fitting.CurveFitter;
+import org.apache.commons.math.optimization.general.LevenbergMarquardtOptimizer;
+import org.apache.log4j.Logger;
+import org.dive4elements.river.artifacts.math.fitting.Function;
+
+/**
+ * @author Gernot Belger
+ */
+final class LevenbergMarquardtFittingOperation extends AbstractFittingOperation implements IFittingOperation {
+
+ private static Logger log = Logger.getLogger(LevenbergMarquardtFittingOperation.class);
+
+ public LevenbergMarquardtFittingOperation(final Function function) {
+ super(function);
+ }
+
+ @Override
+ public Fitting execute(final List<FittingData> data) {
+
+ final Function function = getFunction();
+
+ for (double tolerance = 1e-10; tolerance < 1e-1; tolerance *= 10d) {
+
+ final LevenbergMarquardtOptimizer lmo = new LevenbergMarquardtOptimizer();
+ lmo.setCostRelativeTolerance(tolerance);
+ lmo.setOrthoTolerance(tolerance);
+ lmo.setParRelativeTolerance(tolerance);
+
+ try {
+ final CurveFitter cf = new CurveFitter(lmo);
+
+ for (final FittingData fittingData : data)
+ cf.addObservedPoint(fittingData.q, fittingData.w);
+
+ final double[] parameters = cf.fit(function, function.getInitialGuess());
+
+ return createFitting(parameters, lmo.getChiSquare(), data);
+ }
+ catch (final MathException me) {
+ if (log.isDebugEnabled()) {
+ log.debug("tolerance " + tolerance + " + failed.", me);
+ }
+ }
+ }
+
+ return null;
+ }
+}
\ No newline at end of file
diff -r eb1a29fe823f -r 0380717105ba artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/LogLinearFittingOperation.java
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/LogLinearFittingOperation.java Mon Dec 02 17:56:15 2019 +0100
@@ -0,0 +1,52 @@
+/** Copyright (C) 2017 by Bundesanstalt für Gewässerkunde
+ * Software engineering by
+ * Björnsen Beratende Ingenieure GmbH
+ * Dr. Schumacher Ingenieurbüro für Wasser und Umwelt
+ *
+ * This file is Free Software under the GNU AGPL (>=v3)
+ * and comes with ABSOLUTELY NO WARRANTY! Check out the
+ * documentation coming with Dive4Elements River for details.
+ */
+package org.dive4elements.river.artifacts.model.fixings;
+
+import java.util.List;
+
+import org.dive4elements.river.artifacts.math.fitting.LogLinearAlternative;
+import org.dive4elements.river.artifacts.model.fixings.fitting.LinearLogLinearizedFitting;
+import org.dive4elements.river.artifacts.model.fixings.fitting.LinearLogLinearizedFitting.Result;
+
+/**
+ * @author Gernot Belger
+ */
+final class LogLinearFittingOperation extends AbstractFittingOperation implements IFittingOperation {
+
+ public LogLinearFittingOperation(final LogLinearAlternative function) {
+ super(function);
+ }
+
+ @Override
+ public Fitting execute(final List<FittingData> data) {
+
+ final double[] obsDischarges = new double[data.size()];
+ final double[] obsWaterlevels = new double[data.size()];
+
+ for (int i = 0; i < data.size(); i++) {
+ final FittingData fittingData = data.get(i);
+
+ obsDischarges[i] = fittingData.q;
+ obsWaterlevels[i] = fittingData.w;
+ }
+
+ final LinearLogLinearizedFitting fitting = new LinearLogLinearizedFitting(obsDischarges, obsWaterlevels);
+
+ final Result result = fitting.optimize();
+ final double a = result.getA();
+ final double b = result.getB();
+ final double m = result.getM();
+ final double chiSquare = result.getChiSquare();
+
+ final double[] parameters = new double[] { a, m, b };
+
+ return createFitting(parameters, chiSquare, data);
+ }
+}
\ No newline at end of file
diff -r eb1a29fe823f -r 0380717105ba artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/ArraySubstraction.java
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/ArraySubstraction.java Mon Dec 02 17:56:15 2019 +0100
@@ -0,0 +1,30 @@
+/** Copyright (C) 2017 by Bundesanstalt für Gewässerkunde
+ * Software engineering by
+ * Björnsen Beratende Ingenieure GmbH
+ * Dr. Schumacher Ingenieurbüro für Wasser und Umwelt
+ *
+ * This file is Free Software under the GNU AGPL (>=v3)
+ * and comes with ABSOLUTELY NO WARRANTY! Check out the
+ * documentation coming with Dive4Elements River for details.
+ */
+package org.dive4elements.river.artifacts.model.fixings.fitting;
+
+/**
+ * @author Gernot Belger
+ */
+final class ArraySubstraction {
+
+ private final double[] minuend;
+
+ public ArraySubstraction(final double[] minuend) {
+ this.minuend = minuend;
+ }
+
+ public double[] evaluate(final double[] subtrahend) {
+
+ final double[] difference = new double[this.minuend.length];
+ for (int i = 0; i < difference.length; i++)
+ difference[i] = this.minuend[i] - subtrahend[i];
+ return difference;
+ }
+}
\ No newline at end of file
diff -r eb1a29fe823f -r 0380717105ba artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/ChiSquare.java
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/ChiSquare.java Mon Dec 02 17:56:15 2019 +0100
@@ -0,0 +1,52 @@
+/** Copyright (C) 2017 by Bundesanstalt für Gewässerkunde
+ * Software engineering by
+ * Björnsen Beratende Ingenieure GmbH
+ * Dr. Schumacher Ingenieurbüro für Wasser und Umwelt
+ *
+ * This file is Free Software under the GNU AGPL (>=v3)
+ * and comes with ABSOLUTELY NO WARRANTY! Check out the
+ * documentation coming with Dive4Elements River for details.
+ */
+package org.dive4elements.river.artifacts.model.fixings.fitting;
+
+import java.util.Arrays;
+
+import org.apache.commons.math3.exception.MathIllegalArgumentException;
+import org.apache.commons.math3.stat.descriptive.AbstractUnivariateStatistic;
+import org.apache.commons.math3.stat.descriptive.UnivariateStatistic;
+
+/**
+ * @author Gernot Belger
+ */
+final class ChiSquare extends AbstractUnivariateStatistic {
+
+ private final double[] observation;
+
+ public ChiSquare(final double[] observation) {
+ this.observation = observation;
+ }
+
+ @Override
+ public UnivariateStatistic copy() {
+ /* stateless, hence we can return this */
+ return this;
+ }
+
+ @Override
+ public double evaluate(final double[] values, final int begin, final int length) throws MathIllegalArgumentException {
+ test(values, begin, length, true);
+
+ final double[] residualsWeights = new double[length];
+ Arrays.fill(residualsWeights, 0, length - 1, 1.0);
+
+ double cost = 0;
+
+ for (int i = begin; i < begin + length; i++) {
+ final double residual = values[i + begin] - this.observation[i + begin];
+ // final double wresiduals = residual * FastMath.sqrt(residualsWeights[i]);
+ cost += residualsWeights[i] * residual * residual;
+ }
+
+ return cost;
+ }
+}
\ No newline at end of file
diff -r eb1a29fe823f -r 0380717105ba artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/LinearLogLinearizedFitting.java
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/LinearLogLinearizedFitting.java Mon Dec 02 17:56:15 2019 +0100
@@ -0,0 +1,154 @@
+/** Copyright (C) 2017 by Bundesanstalt für Gewässerkunde
+ * Software engineering by
+ * Björnsen Beratende Ingenieure GmbH
+ * Dr. Schumacher Ingenieurbüro für Wasser und Umwelt
+ *
+ * This file is Free Software under the GNU AGPL (>=v3)
+ * and comes with ABSOLUTELY NO WARRANTY! Check out the
+ * documentation coming with Dive4Elements River for details.
+ */
+package org.dive4elements.river.artifacts.model.fixings.fitting;
+
+import org.apache.commons.math3.optim.MaxEval;
+import org.apache.commons.math3.optim.MaxIter;
+import org.apache.commons.math3.optim.nonlinear.scalar.GoalType;
+import org.apache.commons.math3.optim.univariate.BrentOptimizer;
+import org.apache.commons.math3.optim.univariate.SearchInterval;
+import org.apache.commons.math3.optim.univariate.UnivariateObjectiveFunction;
+import org.apache.commons.math3.optim.univariate.UnivariatePointValuePair;
+
+/**
+ * Alternative solver for W = a * ln(b + m*Q)
+ * Basic approach is to optimize a, by locally solving b and m by linear regression within the error function.
+ *
+ * @author Gernot Belger
+ */
+public final class LinearLogLinearizedFitting {
+
+ public static final class Result {
+
+ private final double a;
+
+ private final double b;
+
+ private final double m;
+
+ private final double error;
+
+ private final int bestAindex;
+
+ private final double chiSquare;
+
+ public Result(final double a, final double b, final double m, final double error, final int bestAindex, final double chiSquare) {
+ this.a = a;
+ this.b = b;
+ this.m = m;
+ this.error = error;
+ this.bestAindex = bestAindex;
+ this.chiSquare = chiSquare;
+ }
+
+ public double getA() {
+ return this.a;
+ }
+
+ public double getB() {
+ return this.b;
+ }
+
+ public double getM() {
+ return this.m;
+ }
+
+ public double getError() {
+ return this.error;
+ }
+
+ public int getBestAindex() {
+ return this.bestAindex;
+ }
+
+ public double getChiSquare() {
+ return this.chiSquare;
+ }
+
+ public Result withBestAIndex(final int i) {
+ return new Result(this.a, this.b, this.m, this.error, i, this.chiSquare);
+ }
+ }
+
+ private final double[] obsDischarges;
+ private final double[] obsWaterlevels;
+ private final TargetFunction targetFunction;
+ private final LinearRegressionSqrtErrorFunction errorFunction;
+
+ public LinearLogLinearizedFitting(final double[] obsDischarges, final double[] obsWaterlevels) {
+ this.obsDischarges = obsDischarges;
+ this.obsWaterlevels = obsWaterlevels;
+
+ this.targetFunction = new TargetFunction(this.obsDischarges);
+ this.errorFunction = new LinearRegressionSqrtErrorFunction(this.obsDischarges, this.obsWaterlevels, this.targetFunction);
+ }
+
+ public Result optimize() {
+
+ return estimateALinear();
+ }
+
+ private Result estimateALinear() {
+
+ double bestA = Double.NaN;
+ double leastError = Double.POSITIVE_INFINITY;
+ int bestAindex = -1;
+
+ // FIXME: a sollte nicht so groß werden
+ for (int i = 0; i < 20; i++) {
+
+ final double aLow = Math.pow(10, i - 1);
+ final double aStart = Math.pow(10, i);
+ final double aHigh = Math.pow(10, i + 1);
+
+ final double[] result = linearOptimize(aLow, aHigh, aStart, this.errorFunction);
+
+ final double a = result[0];
+ final double error = result[1];
+
+ if (error < leastError) {
+ leastError = error;
+ bestA = a;
+ bestAindex = i;
+ }
+ }
+
+ final double[] parameters = this.errorFunction.calc_parameters_trans(bestA);
+ final double b = parameters[0];
+ final double m = parameters[1];
+
+ // FIXME: noch post-optimieren?
+
+ /* calculate chi square */
+ final ChiSquare chiSquareFunc = new ChiSquare(this.obsWaterlevels);
+ final double[] waterlevels = this.targetFunction.calc_stages(bestAindex, b, m);
+ final double chiSquare = chiSquareFunc.evaluate(waterlevels);
+
+ return new Result(bestA, b, m, leastError, bestAindex, chiSquare);
+ }
+
+ private static double[] linearOptimize(final double aLow, final double aHigh, final double aStart, final LinearRegressionSqrtErrorFunction errorFunction) {
+
+ final BrentOptimizer optimizer = new BrentOptimizer(1e-10, 1e-12);
+
+ final SearchInterval searchInterval = new SearchInterval(aLow, aHigh, aStart);
+ final UnivariateObjectiveFunction function = new UnivariateObjectiveFunction(errorFunction);
+
+ final MaxEval maxEval = new MaxEval(Integer.MAX_VALUE);
+ final MaxIter maxIter = new MaxIter(Integer.MAX_VALUE);
+
+ final UnivariatePointValuePair result = optimizer.optimize(searchInterval, function, GoalType.MINIMIZE, maxEval, maxIter);
+
+ final double aEstimation = result.getPoint();
+ final double error = result.getValue();
+
+ return new double[] { aEstimation, error };
+ }
+}
\ No newline at end of file
diff -r eb1a29fe823f -r 0380717105ba artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/LinearRegressionSqrtErrorFunction.java
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/LinearRegressionSqrtErrorFunction.java Mon Dec 02 17:56:15 2019 +0100
@@ -0,0 +1,147 @@
+/** Copyright (C) 2017 by Bundesanstalt für Gewässerkunde
+ * Software engineering by
+ * Björnsen Beratende Ingenieure GmbH
+ * Dr. Schumacher Ingenieurbüro für Wasser und Umwelt
+ *
+ * This file is Free Software under the GNU AGPL (>=v3)
+ * and comes with ABSOLUTELY NO WARRANTY! Check out the
+ * documentation coming with Dive4Elements River for details.
+ */
+package org.dive4elements.river.artifacts.model.fixings.fitting;
+
+import org.apache.commons.math3.analysis.MultivariateFunction;
+import org.apache.commons.math3.analysis.UnivariateFunction;
+import org.apache.commons.math3.stat.descriptive.DescriptiveStatistics;
+import org.apache.commons.math3.stat.regression.SimpleRegression;
+
+/**
+ * @author Gernot Belger
+ *
+ */
+final class LinearRegressionSqrtErrorFunction implements MultivariateFunction, UnivariateFunction {
+
+ private final double[] obsDischarges;
+
+ private final double[] obsWaterlevels;
+
+ private final double waterlevelMean;
+
+ private final double waterlevelDeviation;
+
+ private final double[] standardizedWaterlevel;
+
+ private final SQRTFunction sqrtFunction;
+
+ private final TargetFunction targetFunction;
+
+ public LinearRegressionSqrtErrorFunction(final double[] obsDischarges, final double[] obsWaterlevels, final TargetFunction targetFunction) {
+
+ this.obsDischarges = obsDischarges;
+ this.obsWaterlevels = obsWaterlevels;
+
+ this.sqrtFunction = new SQRTFunction(obsWaterlevels);
+ this.targetFunction = targetFunction;
+
+ final DescriptiveStatistics stats = new DescriptiveStatistics(obsWaterlevels);
+
+ // Compute mean and standard deviation
+ this.waterlevelMean = stats.getMean();
+ this.waterlevelDeviation = stats.getStandardDeviation();
+ // -mittelwert durch standardabweichung /z-transformation)
+
+ // initialize the standardizedSample, which has the same length as the sample
+ this.standardizedWaterlevel = new double[obsWaterlevels.length];
+
+ for (int i = 0; i < obsWaterlevels.length; i++)
+ this.standardizedWaterlevel[i] = (obsWaterlevels[i] - this.waterlevelMean) / this.waterlevelDeviation;
+ }
+
+ @Override
+ public double value(final double[] point) {
+
+ final double a = point[0];
+
+ return value(a);
+ }
+
+ @Override
+ public double value(final double a) {
+
+ final double error = calc_sqrt_trans(a);
+ if (Double.isNaN(error))
+ return Double.POSITIVE_INFINITY;
+
+ return error;
+ }
+
+ private double calc_sqrt_trans(final double a) {
+
+ final double[] waterlevels = calc_stage_trans(a);
+
+ return this.sqrtFunction.calc_sqrt(waterlevels);
+ }
+
+ private double[] calc_stage_trans(final double a) {
+
+ final double[] params = calc_parameters_trans(a);
+
+ final double b = params[0];
+ final double m = params[1];
+
+ return this.targetFunction.calc_stages(a, b, m);
+ }
+
+ public double[] calc_parameters_trans(final double a) {
+
+ // stages_obs_trans = trans_stages(stages_obs, parameter_a);
+ final double[] transWaterlevels = trans_stages(this.obsWaterlevels, a);
+ final double[] discharges = this.obsDischarges;
+
+ // FIXME: z-nromalizing?
+
+ // printArray(discharges);
+ // printArray(transWaterlevels);
+
+ // final DescriptiveStatistics stats = new DescriptiveStatistics(transWaterlevels);
+ //
+ // // Compute mean and standard deviation
+ // final double mean = stats.getMean();
+ // final double deviation = stats.getStandardDeviation();
+ //
+ // final double[] normWaterlevel = new double[transWaterlevels.length];
+ // for (int i = 0; i < transWaterlevels.length; i++)
+ // normWaterlevel[i] = (transWaterlevels[i] - mean) / deviation;
+
+ final SimpleRegression regression = new SimpleRegression(true);
+ for (int i = 0; i < discharges.length; i++)
+ regression.addData(discharges[i], transWaterlevels[i]);
+
+ final double normSlope = regression.getSlope();
+ final double normIntercept = regression.getIntercept();
+
+ // unnormalize
+ // final double slope = normSlope * deviation;
+ // final double intercept = normIntercept + mean;
+ final double slope = normSlope;
+ final double intercept = normIntercept;
+
+ return new double[] { intercept, slope };
+ }
+
+ // private void printArray(final double[] discharges) {
+ // for (final double d : discharges) {
+ // System.out.print(d);
+ // System.out.print(", ");
+ // }
+ // System.out.println(" ");
+ // }
+
+ private double[] trans_stages(final double[] waterlevels, final double a) {
+
+ final double[] transWaterlevels = new double[waterlevels.length];
+ for (int i = 0; i < transWaterlevels.length; i++)
+ transWaterlevels[i] = Math.exp(waterlevels[i] / a);
+
+ return transWaterlevels;
+ }
+}
\ No newline at end of file
diff -r eb1a29fe823f -r 0380717105ba artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/LinearizedFittingTest.java
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/LinearizedFittingTest.java Mon Dec 02 17:56:15 2019 +0100
@@ -0,0 +1,228 @@
+/** Copyright (C) 2017 by Bundesanstalt für Gewässerkunde
+ * Software engineering by
+ * Björnsen Beratende Ingenieure GmbH
+ * Dr. Schumacher Ingenieurbüro für Wasser und Umwelt
+ *
+ * This file is Free Software under the GNU AGPL (>=v3)
+ * and comes with ABSOLUTELY NO WARRANTY! Check out the
+ * documentation coming with Dive4Elements River for details.
+ */
+package org.dive4elements.river.artifacts.model.fixings.fitting;
+
+import java.io.BufferedInputStream;
+import java.io.IOException;
+import java.io.InputStreamReader;
+import java.math.BigDecimal;
+import java.util.ArrayList;
+import java.util.List;
+import java.util.Map.Entry;
+import java.util.NavigableMap;
+import java.util.TreeMap;
+
+import org.apache.commons.lang.StringUtils;
+import org.apache.commons.math3.analysis.MultivariateFunction;
+import org.apache.commons.math3.optim.InitialGuess;
+import org.apache.commons.math3.optim.MaxEval;
+import org.apache.commons.math3.optim.MaxIter;
+import org.apache.commons.math3.optim.PointValuePair;
+import org.apache.commons.math3.optim.SimpleBounds;
+import org.apache.commons.math3.optim.nonlinear.scalar.GoalType;
+import org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunction;
+import org.apache.commons.math3.optim.nonlinear.scalar.noderiv.BOBYQAOptimizer;
+import org.apache.commons.math3.util.Pair;
+import org.dive4elements.river.artifacts.model.fixings.fitting.LinearLogLinearizedFitting.Result;
+
+import au.com.bytecode.opencsv.CSVReader;
+
+/**
+ * @author Gernot Belger
+ */
+public class LinearizedFittingTest {
+
+ public static void main(final String[] args) throws IOException {
+
+ // read test data
+ final NavigableMap<BigDecimal, List<Pair<BigDecimal, BigDecimal>>> testData = readTestData();
+
+ for (final Entry<BigDecimal, List<Pair<BigDecimal, BigDecimal>>> entry : testData.entrySet()) {
+
+ final BigDecimal station = entry.getKey();
+ final List<Pair<BigDecimal, BigDecimal>> testSample = entry.getValue();
+ optimizeA(station, testSample);
+ }
+ }
+
+ private static void optimizeA(final BigDecimal station, final List<Pair<BigDecimal, BigDecimal>> testSample) {
+
+ /* extrakt observations */
+ final double[] obsDischarges = new double[testSample.size()];
+ final double[] obsWaterlevels = new double[testSample.size()];
+ for (int i = 0; i < obsWaterlevels.length; i++) {
+ obsDischarges[i] = testSample.get(i).getKey().doubleValue();
+ obsWaterlevels[i] = testSample.get(i).getValue().doubleValue();
+ }
+
+ final TargetFunction targetFunction = new TargetFunction(obsDischarges);
+ final SqrtErrorFunction sqrtErrorFunction = new SqrtErrorFunction(obsWaterlevels, targetFunction);
+
+ final Result directEstimation = estimateADirect(sqrtErrorFunction);
+ final Result linearEstimation = new LinearLogLinearizedFitting(obsDischarges, obsWaterlevels).optimize();
+
+ final Result directPostEstimation = postOptimize(directEstimation, sqrtErrorFunction);
+ final Result linearPostEstimation = postOptimize(linearEstimation, sqrtErrorFunction);
+
+ printResult(station, testSample, directEstimation);
+ printResult(station, testSample, linearEstimation);
+ printResult(station, testSample, directPostEstimation);
+ printResult(station, testSample, linearPostEstimation);
+ }
+
+ private static Result postOptimize(final Result estimation, final SqrtErrorFunction errorFunction) {
+
+ final double a = estimation.getA();
+ final double b = estimation.getB();
+ final double m = estimation.getM();
+ // final double error = estimation.getError();
+ final double aIndex = estimation.getBestAindex();
+
+ final double aLow = Math.pow(10, aIndex - 1);
+ final double aHigh = Math.pow(10, aIndex + 1);
+
+ return directOptimize(aLow, aHigh, a, b, m, errorFunction);
+ }
+
+ private static void printResult(final BigDecimal station, final List<Pair<BigDecimal, BigDecimal>> testSample, final Result optimize) {
+
+ final double a = optimize.getA();
+ final double b = optimize.getB();
+ final double m = optimize.getM();
+ final double error = optimize.getError();
+
+ System.out.format("%s %.10f %.10f %.10f %.10f", station, a, b, m, error);
+
+ for (final Pair<BigDecimal, BigDecimal> entry : testSample) {
+
+ final BigDecimal waterlevel = entry.getSecond();
+ System.out.print(' ');
+ System.out.print(waterlevel);
+ }
+ System.out.println();
+ }
+
+ private static Result estimateADirect(final SqrtErrorFunction sqrtErrorFunction) {
+
+ Result best = null;
+ double leastError = Double.POSITIVE_INFINITY;
+
+ // iteration über a von 10^^0 bis 10^^20
+ for (int i = 0; i < 20; i++) {
+
+ final double aStart = Math.pow(10, i);
+ final double aLow = Math.pow(10, i - 1);
+ final double aHigh = Math.pow(10, i + 1);
+
+ final Result result = directOptimize(aLow, aHigh, aStart, 1, 1, sqrtErrorFunction);
+ final double error = result.getError();
+
+ if (error < leastError) {
+ leastError = error;
+ best = result.withBestAIndex(i);
+ }
+ }
+
+ return best;
+ }
+
+ private static Result directOptimize(final double aLow, final double aHigh, final double a, final double b, final double m,
+ final SqrtErrorFunction sqrtErrorFunction) {
+
+ // n = 3
+ // [n+2, (n+1)(n+2)/2]
+ // --> [5, 10]
+ final int interpolationPoints = 10;
+ final BOBYQAOptimizer optimizer = new BOBYQAOptimizer(interpolationPoints);
+
+ /* optimization data */
+ final MultivariateFunction function = new MultivariateFunction() {
+
+ @Override
+ public double value(final double[] point) {
+ return sqrtErrorFunction.value(point[0], point[1], point[2]);
+ }
+ };
+
+ final MaxEval maxEval = new MaxEval(Integer.MAX_VALUE);
+ final MaxIter maxIter = new MaxIter(Integer.MAX_VALUE);
+
+ final SimpleBounds bounds = new SimpleBounds(new double[] { aLow, -1e3, 0 }, new double[] { aHigh, 1e3, 1e3 });
+ final double[] startValues = new double[] { a, b, m };
+ final PointValuePair result = optimizer.optimize(GoalType.MINIMIZE, new ObjectiveFunction(function), new InitialGuess(startValues), bounds, maxEval,
+ maxIter);
+
+ final Double error = result.getValue();
+ final double[] point = result.getPoint();
+
+ final double aEstimation = point[0];
+ final double bEstimation = point[1];
+ final double mEstimation = point[2];
+
+ return new Result(aEstimation, bEstimation, mEstimation, error, -1, Double.NaN);
+ }
+
+ private static NavigableMap<BigDecimal, List<Pair<BigDecimal, BigDecimal>>> readTestData() throws IOException {
+
+ final NavigableMap<BigDecimal, List<Pair<BigDecimal, BigDecimal>>> data = new TreeMap<>();
+
+ try (final CSVReader reader = new CSVReader(new InputStreamReader(new BufferedInputStream(LinearizedFittingTest.class.getResourceAsStream("testdata.txt"))),
+ '\t')) {
+
+ final String[] header = reader.readNext();
+ if (header == null)
+ throw new IllegalStateException();
+
+ if (header.length < 2)
+ throw new IllegalStateException();
+
+ while (true) {
+ final String[] line = reader.readNext();
+ if (line == null)
+ break;
+
+ if (line.length != header.length)
+ throw new IllegalStateException();
+
+ final BigDecimal discharge = parseDecimal(line[0]);
+ if (discharge == null)
+ continue;
+
+ for (int column = 1; column < line.length; column++) {
+
+ final BigDecimal station = parseDecimal(header[column]);
+ if (station == null)
+ continue;
+
+ if (!data.containsKey(station))
+ data.put(station, new ArrayList<Pair<BigDecimal, BigDecimal>>());
+ final List<Pair<BigDecimal, BigDecimal>> points = data.get(station);
+
+ final BigDecimal waterlevel = parseDecimal(line[column]);
+ if (waterlevel != null)
+ points.add(Pair.create(discharge, waterlevel));
+ }
+ }
+ }
+
+ return data;
+ }
+
+ private static BigDecimal parseDecimal(final String token) {
+
+ if (StringUtils.isBlank(token))
+ return null;
+
+ if ("nan".equalsIgnoreCase(token))
+ return null;
+
+ return new BigDecimal(token);
+ }
+}
\ No newline at end of file
diff -r eb1a29fe823f -r 0380717105ba artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/SQRTFunction.java
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/SQRTFunction.java Mon Dec 02 17:56:15 2019 +0100
@@ -0,0 +1,30 @@
+/** Copyright (C) 2017 by Bundesanstalt für Gewässerkunde
+ * Software engineering by
+ * Björnsen Beratende Ingenieure GmbH
+ * Dr. Schumacher Ingenieurbüro für Wasser und Umwelt
+ *
+ * This file is Free Software under the GNU AGPL (>=v3)
+ * and comes with ABSOLUTELY NO WARRANTY! Check out the
+ * documentation coming with Dive4Elements River for details.
+ */
+package org.dive4elements.river.artifacts.model.fixings.fitting;
+
+import org.apache.commons.math3.stat.StatUtils;
+
+/**
+ * @author Gernot Belger
+ */
+final class SQRTFunction {
+ private final ArraySubstraction substraction;
+
+ public SQRTFunction(final double[] observations) {
+ this.substraction = new ArraySubstraction(observations);
+ }
+
+ public double calc_sqrt(final double[] values) {
+
+ final double[] difference = this.substraction.evaluate(values);
+
+ return StatUtils.sumSq(difference);
+ }
+}
\ No newline at end of file
diff -r eb1a29fe823f -r 0380717105ba artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/SqrtErrorFunction.java
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/SqrtErrorFunction.java Mon Dec 02 17:56:15 2019 +0100
@@ -0,0 +1,36 @@
+/** Copyright (C) 2017 by Bundesanstalt für Gewässerkunde
+ * Software engineering by
+ * Björnsen Beratende Ingenieure GmbH
+ * Dr. Schumacher Ingenieurbüro für Wasser und Umwelt
+ *
+ * This file is Free Software under the GNU AGPL (>=v3)
+ * and comes with ABSOLUTELY NO WARRANTY! Check out the
+ * documentation coming with Dive4Elements River for details.
+ */
+package org.dive4elements.river.artifacts.model.fixings.fitting;
+
+/**
+ * @author Gernot Belger
+ */
+final class SqrtErrorFunction {
+
+ private final SQRTFunction sqrtFunction;
+
+ private final TargetFunction targetFunction;
+
+ public SqrtErrorFunction(final double[] obsWaterlevels, final TargetFunction targetFunction) {
+
+ this.targetFunction = targetFunction;
+ this.sqrtFunction = new SQRTFunction(obsWaterlevels);
+ }
+
+ public double value(final double a, final double b, final double m) {
+ return calc_sqrt_trans(a, b, m);
+ }
+
+ private double calc_sqrt_trans(final double a, final double b, final double m) {
+
+ final double[] waterlevels = this.targetFunction.calc_stages(a, b, m);
+ return this.sqrtFunction.calc_sqrt(waterlevels);
+ }
+}
\ No newline at end of file
diff -r eb1a29fe823f -r 0380717105ba artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/TargetFunction.java
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/TargetFunction.java Mon Dec 02 17:56:15 2019 +0100
@@ -0,0 +1,31 @@
+/** Copyright (C) 2017 by Bundesanstalt für Gewässerkunde
+ * Software engineering by
+ * Björnsen Beratende Ingenieure GmbH
+ * Dr. Schumacher Ingenieurbüro für Wasser und Umwelt
+ *
+ * This file is Free Software under the GNU AGPL (>=v3)
+ * and comes with ABSOLUTELY NO WARRANTY! Check out the
+ * documentation coming with Dive4Elements River for details.
+ */
+package org.dive4elements.river.artifacts.model.fixings.fitting;
+
+/**
+ * @author Gernot Belger
+ */
+final class TargetFunction {
+
+ private final double[] discharge;
+
+ public TargetFunction(final double[] discharge) {
+ this.discharge = discharge;
+ }
+
+ public double[] calc_stages(final double a, final double b, final double m) {
+
+ final double[] waterlevels = new double[this.discharge.length];
+ for (int i = 0; i < this.discharge.length; i++)
+ waterlevels[i] = a * Math.log(b + m * this.discharge[i]);
+
+ return waterlevels;
+ }
+}
\ No newline at end of file
diff -r eb1a29fe823f -r 0380717105ba artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/testdata.txt
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/testdata.txt Mon Dec 02 17:56:15 2019 +0100
@@ -0,0 +1,7 @@
+Q 541 541.1 541.2 541.3 541.4 541.5 541.6 541.7 541.8 541.9 542 542.1 542.2 542.3 542.4 542.5 542.6 542.7 542.8 542.9 543
+750.0 71.48 71.425 71.37 71.305 71.24 71.18 71.12 71.07 71.02 70.94 70.86 70.725 70.59 70.53 70.39 70.35 70.31 70.29 70.27 70.235 70.2
+750.0 71.46 71.405 71.35 71.295 nan nan nan 71.035 70.95 70.86 70.77 70.675 70.58 70.49 70.47 70.385 70.3 70.24 70.18 70.145 70.11
+1555.0 72.31 72.24 72.17 72.09 72.01 71.93 71.85 71.74 71.63 71.545 71.46 71.395 71.33 71.28 71.22 71.18 71.14 71.12 71.1 71.035 70.97
+1583.0 72.24 72.178 72.116 72.054 71.992 71.93 71.852 71.774 71.696 71.618 71.54 71.474 71.408 71.342 71.276 71.21 71.176 71.142 71.108 71.07 71.04
+1618.0 72.38 72.305 72.23 72.145 72.06 71.98 71.90 71.805 71.71 71.61 71.50 71.42 71.34 71.295 71.25 71.23 71.21 71.185 71.16 71.11 71.06
+5682.0 75.99 75.97 75.96 75.94 75.93 75.91 75.88 75.86 75.83 75.81 75.78 75.76 75.73 75.71 75.68 75.66 75.62 75.59 75.55 75.52 75.48
diff -r eb1a29fe823f -r 0380717105ba gwt-client/src/main/java/org/dive4elements/river/client/client/ui/fixation/FixFunctionSelect.java
--- a/gwt-client/src/main/java/org/dive4elements/river/client/client/ui/fixation/FixFunctionSelect.java Mon Dec 02 17:21:43 2019 +0100
+++ b/gwt-client/src/main/java/org/dive4elements/river/client/client/ui/fixation/FixFunctionSelect.java Mon Dec 02 17:56:15 2019 +0100
@@ -33,7 +33,7 @@
funcDesc.put("log", "W(Q) = m*ln(Q + b)");
funcDesc.put("linear", "W(Q) = m * Q + b");
funcDesc.put("log-linear", "W(Q) = a*ln(m*Q+b)");
- funcDesc.put("log-linear-alternative", "W(Q) = a*ln(m*Q+b) (linearisiert)");
+ funcDesc.put("log-linear-linearisiert", "W(Q) = a*ln(m*Q+b) (linearisiert)");
funcDesc.put("exp", "W(Q) = m * a^Q + b");
funcDesc.put("quad", "W(Q) = n*Q^2+m*Q+b");
funcDesc.put("pow", "W(Q) = a * Q^c + d");
More information about the Dive4Elements-commits
mailing list