/** 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 };
    }
}