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.analysis.MultivariateFunction;
g@9646: import org.apache.commons.math3.analysis.UnivariateFunction;
g@9646: import org.apache.commons.math3.stat.descriptive.DescriptiveStatistics;
g@9646: import org.apache.commons.math3.stat.regression.SimpleRegression;
g@9646: 
g@9646: /**
g@9646:  * @author Gernot Belger
g@9646:  *
g@9646:  */
g@9646: final class LinearRegressionSqrtErrorFunction implements MultivariateFunction, UnivariateFunction {
g@9646: 
g@9646:     private final double[] obsDischarges;
g@9646: 
g@9646:     private final double[] obsWaterlevels;
g@9646: 
g@9646:     private final double waterlevelMean;
g@9646: 
g@9646:     private final double waterlevelDeviation;
g@9646: 
g@9646:     private final double[] standardizedWaterlevel;
g@9646: 
g@9646:     private final SQRTFunction sqrtFunction;
g@9646: 
g@9646:     private final TargetFunction targetFunction;
g@9646: 
g@9646:     public LinearRegressionSqrtErrorFunction(final double[] obsDischarges, final double[] obsWaterlevels, final TargetFunction targetFunction) {
g@9646: 
g@9646:         this.obsDischarges = obsDischarges;
g@9646:         this.obsWaterlevels = obsWaterlevels;
g@9646: 
g@9646:         this.sqrtFunction = new SQRTFunction(obsWaterlevels);
g@9646:         this.targetFunction = targetFunction;
g@9646: 
g@9646:         final DescriptiveStatistics stats = new DescriptiveStatistics(obsWaterlevels);
g@9646: 
g@9646:         // Compute mean and standard deviation
g@9646:         this.waterlevelMean = stats.getMean();
g@9646:         this.waterlevelDeviation = stats.getStandardDeviation();
g@9646:         // -mittelwert durch standardabweichung /z-transformation)
g@9646: 
g@9646:         // initialize the standardizedSample, which has the same length as the sample
g@9646:         this.standardizedWaterlevel = new double[obsWaterlevels.length];
g@9646: 
g@9646:         for (int i = 0; i < obsWaterlevels.length; i++)
g@9646:             this.standardizedWaterlevel[i] = (obsWaterlevels[i] - this.waterlevelMean) / this.waterlevelDeviation;
g@9646:     }
g@9646: 
g@9646:     @Override
g@9646:     public double value(final double[] point) {
g@9646: 
g@9646:         final double a = point[0];
g@9646: 
g@9646:         return value(a);
g@9646:     }
g@9646: 
g@9646:     @Override
g@9646:     public double value(final double a) {
g@9646: 
g@9646:         final double error = calc_sqrt_trans(a);
g@9646:         if (Double.isNaN(error))
g@9646:             return Double.POSITIVE_INFINITY;
g@9646: 
g@9646:         return error;
g@9646:     }
g@9646: 
g@9646:     private double calc_sqrt_trans(final double a) {
g@9646: 
g@9646:         final double[] waterlevels = calc_stage_trans(a);
g@9646: 
g@9646:         return this.sqrtFunction.calc_sqrt(waterlevels);
g@9646:     }
g@9646: 
g@9646:     private double[] calc_stage_trans(final double a) {
g@9646: 
g@9646:         final double[] params = calc_parameters_trans(a);
g@9646: 
g@9646:         final double b = params[0];
g@9646:         final double m = params[1];
g@9646: 
g@9646:         return this.targetFunction.calc_stages(a, b, m);
g@9646:     }
g@9646: 
g@9646:     public double[] calc_parameters_trans(final double a) {
g@9646: 
g@9646:         // stages_obs_trans = trans_stages(stages_obs, parameter_a);
g@9646:         final double[] transWaterlevels = trans_stages(this.obsWaterlevels, a);
g@9646:         final double[] discharges = this.obsDischarges;
g@9646: 
g@9646:         // FIXME: z-nromalizing?
g@9646: 
g@9646:         // printArray(discharges);
g@9646:         // printArray(transWaterlevels);
g@9646: 
g@9646:         // final DescriptiveStatistics stats = new DescriptiveStatistics(transWaterlevels);
g@9646:         //
g@9646:         // // Compute mean and standard deviation
g@9646:         // final double mean = stats.getMean();
g@9646:         // final double deviation = stats.getStandardDeviation();
g@9646:         //
g@9646:         // final double[] normWaterlevel = new double[transWaterlevels.length];
g@9646:         // for (int i = 0; i < transWaterlevels.length; i++)
g@9646:         // normWaterlevel[i] = (transWaterlevels[i] - mean) / deviation;
g@9646: 
g@9646:         final SimpleRegression regression = new SimpleRegression(true);
g@9646:         for (int i = 0; i < discharges.length; i++)
g@9646:             regression.addData(discharges[i], transWaterlevels[i]);
g@9646: 
g@9646:         final double normSlope = regression.getSlope();
g@9646:         final double normIntercept = regression.getIntercept();
g@9646: 
g@9646:         // unnormalize
g@9646:         // final double slope = normSlope * deviation;
g@9646:         // final double intercept = normIntercept + mean;
g@9646:         final double slope = normSlope;
g@9646:         final double intercept = normIntercept;
g@9646: 
g@9646:         return new double[] { intercept, slope };
g@9646:     }
g@9646: 
g@9646:     // private void printArray(final double[] discharges) {
g@9646:     // for (final double d : discharges) {
g@9646:     // System.out.print(d);
g@9646:     // System.out.print(", ");
g@9646:     // }
g@9646:     // System.out.println(" ");
g@9646:     // }
g@9646: 
g@9646:     private double[] trans_stages(final double[] waterlevels, final double a) {
g@9646: 
g@9646:         final double[] transWaterlevels = new double[waterlevels.length];
g@9646:         for (int i = 0; i < transWaterlevels.length; i++)
g@9646:             transWaterlevels[i] = Math.exp(waterlevels[i] / a);
g@9646: 
g@9646:         return transWaterlevels;
g@9646:     }
g@9646: }