Mercurial > dive4elements > river
view 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 source
/** 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 }; } }