diff artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/LinearLogLinearizedFitting.java @ 9646:0380717105ba

Implemented alternative fitting strategy for Log-Linear function.
author Gernot Belger <g.belger@bjoernsen.de>
date Mon, 02 Dec 2019 17:56:15 +0100
parents
children
line wrap: on
line diff
--- /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

http://dive4elements.wald.intevation.org