Mercurial > dive4elements > river

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

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
comparison
equal deleted inserted replaced
9645:eb1a29fe823f 9646:0380717105ba
1 /** Copyright (C) 2017 by Bundesanstalt für Gewässerkunde
2 * Software engineering by
3 * Björnsen Beratende Ingenieure GmbH
4 * Dr. Schumacher Ingenieurbüro für Wasser und Umwelt
5 *
6 * This file is Free Software under the GNU AGPL (>=v3)
7 * and comes with ABSOLUTELY NO WARRANTY! Check out the
8 * documentation coming with Dive4Elements River for details.
9 */
10 package org.dive4elements.river.artifacts.model.fixings.fitting;
11
12 import org.apache.commons.math3.optim.MaxEval;
13 import org.apache.commons.math3.optim.MaxIter;
14 import org.apache.commons.math3.optim.nonlinear.scalar.GoalType;
15 import org.apache.commons.math3.optim.univariate.BrentOptimizer;
16 import org.apache.commons.math3.optim.univariate.SearchInterval;
17 import org.apache.commons.math3.optim.univariate.UnivariateObjectiveFunction;
18 import org.apache.commons.math3.optim.univariate.UnivariatePointValuePair;
19
20 /**
21 * Alternative solver for W = a * ln(b + m*Q)
22 * Basic approach is to optimize a, by locally solving b and m by linear regression within the error function.
23 *
24 * @author Gernot Belger
25 */
26 public final class LinearLogLinearizedFitting {
27
28 public static final class Result {
29
30 private final double a;
31
32 private final double b;
33
34 private final double m;
35
36 private final double error;
37
38 private final int bestAindex;
39
40 private final double chiSquare;
41
42 public Result(final double a, final double b, final double m, final double error, final int bestAindex, final double chiSquare) {
43 this.a = a;
44 this.b = b;
45 this.m = m;
46 this.error = error;
47 this.bestAindex = bestAindex;
48 this.chiSquare = chiSquare;
49 }
50
51 public double getA() {
52 return this.a;
53 }
54
55 public double getB() {
56 return this.b;
57 }
58
59 public double getM() {
60 return this.m;
61 }
62
63 public double getError() {
64 return this.error;
65 }
66
67 public int getBestAindex() {
68 return this.bestAindex;
69 }
70
71 public double getChiSquare() {
72 return this.chiSquare;
73 }
74
75 public Result withBestAIndex(final int i) {
76 return new Result(this.a, this.b, this.m, this.error, i, this.chiSquare);
77 }
78 }
79
80 private final double[] obsDischarges;
81 private final double[] obsWaterlevels;
82 private final TargetFunction targetFunction;
83 private final LinearRegressionSqrtErrorFunction errorFunction;
84
85 public LinearLogLinearizedFitting(final double[] obsDischarges, final double[] obsWaterlevels) {
86 this.obsDischarges = obsDischarges;
87 this.obsWaterlevels = obsWaterlevels;
88
89 this.targetFunction = new TargetFunction(this.obsDischarges);
90 this.errorFunction = new LinearRegressionSqrtErrorFunction(this.obsDischarges, this.obsWaterlevels, this.targetFunction);
91 }
92
93 public Result optimize() {
94
95 return estimateALinear();
96 }
97
98 private Result estimateALinear() {
99
100 double bestA = Double.NaN;
101 double leastError = Double.POSITIVE_INFINITY;
102 int bestAindex = -1;
103
104 // FIXME: a sollte nicht so groß werden
105 for (int i = 0; i < 20; i++) {
106
107 final double aLow = Math.pow(10, i - 1);
108 final double aStart = Math.pow(10, i);
109 final double aHigh = Math.pow(10, i + 1);
110
111 final double[] result = linearOptimize(aLow, aHigh, aStart, this.errorFunction);
112
113 final double a = result[0];
114 final double error = result[1];
115
116 if (error < leastError) {
117 leastError = error;
118 bestA = a;
119 bestAindex = i;
120 }
121 }
122
123 final double[] parameters = this.errorFunction.calc_parameters_trans(bestA);
124 final double b = parameters[0];
125 final double m = parameters[1];
126
127 // FIXME: noch post-optimieren?
128
129 /* calculate chi square */
130 final ChiSquare chiSquareFunc = new ChiSquare(this.obsWaterlevels);
131 final double[] waterlevels = this.targetFunction.calc_stages(bestAindex, b, m);
132 final double chiSquare = chiSquareFunc.evaluate(waterlevels);
133
134 return new Result(bestA, b, m, leastError, bestAindex, chiSquare);
135 }
136
137 private static double[] linearOptimize(final double aLow, final double aHigh, final double aStart, final LinearRegressionSqrtErrorFunction errorFunction) {
138
139 final BrentOptimizer optimizer = new BrentOptimizer(1e-10, 1e-12);
140
141 final SearchInterval searchInterval = new SearchInterval(aLow, aHigh, aStart);
142 final UnivariateObjectiveFunction function = new UnivariateObjectiveFunction(errorFunction);
143
144 final MaxEval maxEval = new MaxEval(Integer.MAX_VALUE);
145 final MaxIter maxIter = new MaxIter(Integer.MAX_VALUE);
146
147 final UnivariatePointValuePair result = optimizer.optimize(searchInterval, function, GoalType.MINIMIZE, maxEval, maxIter);
148
149 final double aEstimation = result.getPoint();
150 final double error = result.getValue();
151
152 return new double[] { aEstimation, error };
153 }
154 }
http://dive4elements.wald.intevation.org