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