g@9646: /** Copyright (C) 2017 by Bundesanstalt für Gewässerkunde
g@9646:  * Software engineering by
g@9646:  *  Björnsen Beratende Ingenieure GmbH
g@9646:  *  Dr. Schumacher Ingenieurbüro für Wasser und Umwelt
g@9646:  *
g@9646:  * This file is Free Software under the GNU AGPL (>=v3)
g@9646:  * and comes with ABSOLUTELY NO WARRANTY! Check out the
g@9646:  * documentation coming with Dive4Elements River for details.
g@9646:  */
g@9646: package org.dive4elements.river.artifacts.model.fixings.fitting;
g@9646: 
g@9646: import org.apache.commons.math3.optim.MaxEval;
g@9646: import org.apache.commons.math3.optim.MaxIter;
g@9646: import org.apache.commons.math3.optim.nonlinear.scalar.GoalType;
g@9646: import org.apache.commons.math3.optim.univariate.BrentOptimizer;
g@9646: import org.apache.commons.math3.optim.univariate.SearchInterval;
g@9646: import org.apache.commons.math3.optim.univariate.UnivariateObjectiveFunction;
g@9646: import org.apache.commons.math3.optim.univariate.UnivariatePointValuePair;
g@9646: 
g@9646: /**
g@9646:  * Alternative solver for W = a * ln(b + m*Q)
g@9646:  * Basic approach is to optimize a, by locally solving b and m by linear regression within the error function.
g@9646:  *
g@9646:  * @author Gernot Belger
g@9646:  */
g@9646: public final class LinearLogLinearizedFitting {
g@9646: 
g@9646:     public static final class Result {
g@9646: 
g@9646:         private final double a;
g@9646: 
g@9646:         private final double b;
g@9646: 
g@9646:         private final double m;
g@9646: 
g@9646:         private final double error;
g@9646: 
g@9646:         private final int bestAindex;
g@9646: 
g@9646:         private final double chiSquare;
g@9646: 
g@9646:         public Result(final double a, final double b, final double m, final double error, final int bestAindex, final double chiSquare) {
g@9646:             this.a = a;
g@9646:             this.b = b;
g@9646:             this.m = m;
g@9646:             this.error = error;
g@9646:             this.bestAindex = bestAindex;
g@9646:             this.chiSquare = chiSquare;
g@9646:         }
g@9646: 
g@9646:         public double getA() {
g@9646:             return this.a;
g@9646:         }
g@9646: 
g@9646:         public double getB() {
g@9646:             return this.b;
g@9646:         }
g@9646: 
g@9646:         public double getM() {
g@9646:             return this.m;
g@9646:         }
g@9646: 
g@9646:         public double getError() {
g@9646:             return this.error;
g@9646:         }
g@9646: 
g@9646:         public int getBestAindex() {
g@9646:             return this.bestAindex;
g@9646:         }
g@9646: 
g@9646:         public double getChiSquare() {
g@9646:             return this.chiSquare;
g@9646:         }
g@9646: 
g@9646:         public Result withBestAIndex(final int i) {
g@9646:             return new Result(this.a, this.b, this.m, this.error, i, this.chiSquare);
g@9646:         }
g@9646:     }
g@9646: 
g@9646:     private final double[] obsDischarges;
g@9646:     private final double[] obsWaterlevels;
g@9646:     private final TargetFunction targetFunction;
g@9646:     private final LinearRegressionSqrtErrorFunction errorFunction;
g@9646: 
g@9646:     public LinearLogLinearizedFitting(final double[] obsDischarges, final double[] obsWaterlevels) {
g@9646:         this.obsDischarges = obsDischarges;
g@9646:         this.obsWaterlevels = obsWaterlevels;
g@9646: 
g@9646:         this.targetFunction = new TargetFunction(this.obsDischarges);
g@9646:         this.errorFunction = new LinearRegressionSqrtErrorFunction(this.obsDischarges, this.obsWaterlevels, this.targetFunction);
g@9646:     }
g@9646: 
g@9646:     public Result optimize() {
g@9646: 
g@9646:         return estimateALinear();
g@9646:     }
g@9646: 
g@9646:     private Result estimateALinear() {
g@9646: 
g@9646:         double bestA = Double.NaN;
g@9646:         double leastError = Double.POSITIVE_INFINITY;
g@9646:         int bestAindex = -1;
g@9646: 
g@9646:         // FIXME: a sollte nicht so groß werden
g@9646:         for (int i = 0; i < 20; i++) {
g@9646: 
g@9646:             final double aLow = Math.pow(10, i - 1);
g@9646:             final double aStart = Math.pow(10, i);
g@9646:             final double aHigh = Math.pow(10, i + 1);
g@9646: 
g@9646:             final double[] result = linearOptimize(aLow, aHigh, aStart, this.errorFunction);
g@9646: 
g@9646:             final double a = result[0];
g@9646:             final double error = result[1];
g@9646: 
g@9646:             if (error < leastError) {
g@9646:                 leastError = error;
g@9646:                 bestA = a;
g@9646:                 bestAindex = i;
g@9646:             }
g@9646:         }
g@9646: 
g@9646:         final double[] parameters = this.errorFunction.calc_parameters_trans(bestA);
g@9646:         final double b = parameters[0];
g@9646:         final double m = parameters[1];
g@9646: 
g@9646:         // FIXME: noch post-optimieren?
g@9646: 
g@9646:         /* calculate chi square */
g@9646:         final ChiSquare chiSquareFunc = new ChiSquare(this.obsWaterlevels);
g@9646:         final double[] waterlevels = this.targetFunction.calc_stages(bestAindex, b, m);
g@9646:         final double chiSquare = chiSquareFunc.evaluate(waterlevels);
g@9646: 
g@9646:         return new Result(bestA, b, m, leastError, bestAindex, chiSquare);
g@9646:     }
g@9646: 
g@9646:     private static double[] linearOptimize(final double aLow, final double aHigh, final double aStart, final LinearRegressionSqrtErrorFunction errorFunction) {
g@9646: 
g@9646:         final BrentOptimizer optimizer = new BrentOptimizer(1e-10, 1e-12);
g@9646: 
g@9646:         final SearchInterval searchInterval = new SearchInterval(aLow, aHigh, aStart);
g@9646:         final UnivariateObjectiveFunction function = new UnivariateObjectiveFunction(errorFunction);
g@9646: 
g@9646:         final MaxEval maxEval = new MaxEval(Integer.MAX_VALUE);
g@9646:         final MaxIter maxIter = new MaxIter(Integer.MAX_VALUE);
g@9646: 
g@9646:         final UnivariatePointValuePair result = optimizer.optimize(searchInterval, function, GoalType.MINIMIZE, maxEval, maxIter);
g@9646: 
g@9646:         final double aEstimation = result.getPoint();
g@9646:         final double error = result.getValue();
g@9646: 
g@9646:         return new double[] { aEstimation, error };
g@9646:     }
g@9646: }