Mercurial > dive4elements > river
comparison 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 |
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 } |