Mercurial > dive4elements > river
changeset 9646:0380717105ba
Implemented alternative fitting strategy for Log-Linear function.
line wrap: on
line diff
--- a/artifacts/pom.xml Mon Dec 02 17:21:43 2019 +0100 +++ b/artifacts/pom.xml Mon Dec 02 17:56:15 2019 +0100 @@ -133,7 +133,7 @@ <groupId>org.apache.commons</groupId> <artifactId>commons-math</artifactId> <version>2.2</version> - </dependency> + </dependency> <dependency> <groupId>com.h2database</groupId> <artifactId>h2</artifactId> @@ -195,6 +195,16 @@ <artifactId>groovy-all</artifactId> <version>1.6.0</version> </dependency> + <dependency> + <groupId>org.apache.commons</groupId> + <artifactId>commons-math3</artifactId> + <version>3.6.1</version> + </dependency> + <dependency> + <groupId>org.hipparchus</groupId> + <artifactId>hipparchus-optim</artifactId> + <version>1.6</version> + </dependency> </dependencies> <repositories> <repository>
--- a/artifacts/src/main/java/org/dive4elements/river/artifacts/math/fitting/LogLinearAlternative.java Mon Dec 02 17:21:43 2019 +0100 +++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/math/fitting/LogLinearAlternative.java Mon Dec 02 17:56:15 2019 +0100 @@ -15,7 +15,7 @@ */ public class LogLinearAlternative extends AbstractLogLinear { - public static final String NAME = "log-linear-alternative"; + public static final String NAME = "log-linear-linearisiert"; public static final Function INSTANCE = new LogLinearAlternative(); public LogLinearAlternative() {
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/AbstractFittingOperation.java Mon Dec 02 17:56:15 2019 +0100 @@ -0,0 +1,86 @@ +/** 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; + +import java.util.List; + +import org.apache.commons.math3.stat.descriptive.moment.StandardDeviation; +import org.dive4elements.river.artifacts.math.fitting.Function; + +/** + * @author Gernot Belger + */ +abstract class AbstractFittingOperation { + + private final Function function; + + public AbstractFittingOperation(final Function function) { + this.function = function; + } + + protected final Function getFunction() { + return this.function; + } + + protected final Fitting createFitting(final double[] parameters, final double chiSquare, final List<FittingData> data) { + + final org.dive4elements.river.artifacts.math.Function instance = this.function.instantiate(parameters); + final double[] waterlevels = calculateWaterlevels(instance, data); + final double maxQ = calculateMaxQ(data); + final double stdDev = calculateStandardDeviation(waterlevels, data); + + return new Fitting(parameters, stdDev, chiSquare, maxQ); + } + + private double calculateMaxQ(final List<FittingData> data) { + + double maxQ = -Double.MAX_VALUE; + + for (int i = 0; i < data.size(); i++) { + + final FittingData fittingData = data.get(i); + + if (fittingData.q > maxQ) + maxQ = fittingData.q; + } + + return maxQ; + } + + private double calculateStandardDeviation(final double[] waterlevels, final List<FittingData> data) { + + final StandardDeviation std = new StandardDeviation(); + + for (int i = 0; i < data.size(); i++) { + + final FittingData fittingData = data.get(i); + + final double diff = (waterlevels[i] - fittingData.q) * 100.0; // * 100 for whatever reason + + std.increment(diff); + } + + return std.getResult(); + } + + private double[] calculateWaterlevels(final org.dive4elements.river.artifacts.math.Function instance, final List<FittingData> data) { + + final double[] waterlevels = new double[data.size()]; + + for (int i = 0; i < data.size(); i++) { + + final FittingData fittingData = data.get(i); + + waterlevels[i] = instance.value(fittingData.q); + } + + return waterlevels; + } +} \ No newline at end of file
--- a/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/Fitting.java Mon Dec 02 17:21:43 2019 +0100 +++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/Fitting.java Mon Dec 02 17:56:15 2019 +0100 @@ -12,13 +12,10 @@ import java.util.Date; import java.util.List; -import org.apache.commons.math.MathException; -import org.apache.commons.math.optimization.fitting.CurveFitter; -import org.apache.commons.math.optimization.general.LevenbergMarquardtOptimizer; -import org.apache.commons.math.stat.descriptive.moment.StandardDeviation; import org.apache.log4j.Logger; import org.dive4elements.river.artifacts.math.GrubbsOutlier; import org.dive4elements.river.artifacts.math.fitting.Function; +import org.dive4elements.river.artifacts.math.fitting.LogLinearAlternative; public class Fitting { private static Logger log = Logger.getLogger(Fitting.class); @@ -28,20 +25,6 @@ QWD create(double q, double w, double deltaW, boolean isOutlier); } - private static final class FittingData { - public final FixingColumnWithData event; - public final double w; - public final double q; - public final boolean isInterpolated; - - public FittingData(final FixingColumnWithData event, final double w, final double q, final boolean isInterpolated) { - this.event = event; - this.w = w; - this.q = q; - this.isInterpolated = isInterpolated; - } - } - private final double chiSqr; private final double[] parameters; @@ -101,37 +84,14 @@ final List<FittingData> outliers = new ArrayList<>(data.size()); org.dive4elements.river.artifacts.math.Function instance = null; - LevenbergMarquardtOptimizer lmo = null; - double[] parameters = null; + Fitting fitting = null; + final IFittingOperation operation = createOperation(function); for (;;) { - parameters = null; - - for (double tolerance = 1e-10; tolerance < 1e-1; tolerance *= 10d) { - - lmo = new LevenbergMarquardtOptimizer(); - lmo.setCostRelativeTolerance(tolerance); - lmo.setOrthoTolerance(tolerance); - lmo.setParRelativeTolerance(tolerance); - - try { - final CurveFitter cf = new CurveFitter(lmo); + fitting = operation.execute(data); - for (final FittingData fittingData : data) - cf.addObservedPoint(fittingData.q, fittingData.w); - - parameters = cf.fit(function, function.getInitialGuess()); - break; - } - catch (final MathException me) { - if (log.isDebugEnabled()) { - log.debug("tolerance " + tolerance + " + failed.", me); - } - } - } - - if (parameters == null) { + if (fitting == null) { /* * log.debug("Parameters is null"); * for (int i = 0, N = xs.size(); i < N; ++i) { @@ -142,6 +102,7 @@ } // This is the parameterized function for a given km. + final double[] parameters = fitting.getParameters(); instance = function.instantiate(parameters); if (!checkOutliers) @@ -180,30 +141,24 @@ } /* - * calculate dW of used values against the resulting function and add them to results , also calculate standard * - * deviation + * calculate dW of used values against the resulting function and add them to results */ - final StandardDeviation stdDev = new StandardDeviation(); - double maxQ = -Double.MAX_VALUE; - for (final FittingData fittingData : data) { final QWD qwd = createQWD(fittingData, instance, false); qwds.add(qwd); resultColumns.addQWD(fittingData.event, km, qwd); - - stdDev.increment(qwd.getDeltaW()); - - final double q = qwd.getQ(); - if (q > maxQ) - maxQ = q; } - final double standardDeviation = stdDev.getResult(); + return fitting; + } - final double chiSqr = lmo.getChiSquare(); + private static IFittingOperation createOperation(final Function function) { - return new Fitting(parameters, standardDeviation, chiSqr, maxQ); + if (function instanceof LogLinearAlternative) + return new LogLinearFittingOperation((LogLinearAlternative) function); + + return new LevenbergMarquardtFittingOperation(function); } private static QWD createQWD(final FittingData data, final org.dive4elements.river.artifacts.math.Function function, final boolean isOutlier) {
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/FittingData.java Mon Dec 02 17:56:15 2019 +0100 @@ -0,0 +1,24 @@ +/** 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; + +final class FittingData { + public final FixingColumnWithData event; + public final double w; + public final double q; + public final boolean isInterpolated; + + public FittingData(final FixingColumnWithData event, final double w, final double q, final boolean isInterpolated) { + this.event = event; + this.w = w; + this.q = q; + this.isInterpolated = isInterpolated; + } +} \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/IFittingOperation.java Mon Dec 02 17:56:15 2019 +0100 @@ -0,0 +1,24 @@ +/** 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; + +import java.util.List; + +/** + * Abstraction for different fitting approaches depending on target function. + * + * @author Gernot Belger + * + */ +public interface IFittingOperation { + + Fitting execute(List<FittingData> data); + +} \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/LevenbergMarquardtFittingOperation.java Mon Dec 02 17:56:15 2019 +0100 @@ -0,0 +1,62 @@ +/** 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; + +import java.util.List; + +import org.apache.commons.math.MathException; +import org.apache.commons.math.optimization.fitting.CurveFitter; +import org.apache.commons.math.optimization.general.LevenbergMarquardtOptimizer; +import org.apache.log4j.Logger; +import org.dive4elements.river.artifacts.math.fitting.Function; + +/** + * @author Gernot Belger + */ +final class LevenbergMarquardtFittingOperation extends AbstractFittingOperation implements IFittingOperation { + + private static Logger log = Logger.getLogger(LevenbergMarquardtFittingOperation.class); + + public LevenbergMarquardtFittingOperation(final Function function) { + super(function); + } + + @Override + public Fitting execute(final List<FittingData> data) { + + final Function function = getFunction(); + + for (double tolerance = 1e-10; tolerance < 1e-1; tolerance *= 10d) { + + final LevenbergMarquardtOptimizer lmo = new LevenbergMarquardtOptimizer(); + lmo.setCostRelativeTolerance(tolerance); + lmo.setOrthoTolerance(tolerance); + lmo.setParRelativeTolerance(tolerance); + + try { + final CurveFitter cf = new CurveFitter(lmo); + + for (final FittingData fittingData : data) + cf.addObservedPoint(fittingData.q, fittingData.w); + + final double[] parameters = cf.fit(function, function.getInitialGuess()); + + return createFitting(parameters, lmo.getChiSquare(), data); + } + catch (final MathException me) { + if (log.isDebugEnabled()) { + log.debug("tolerance " + tolerance + " + failed.", me); + } + } + } + + return null; + } +} \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/LogLinearFittingOperation.java Mon Dec 02 17:56:15 2019 +0100 @@ -0,0 +1,52 @@ +/** 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; + +import java.util.List; + +import org.dive4elements.river.artifacts.math.fitting.LogLinearAlternative; +import org.dive4elements.river.artifacts.model.fixings.fitting.LinearLogLinearizedFitting; +import org.dive4elements.river.artifacts.model.fixings.fitting.LinearLogLinearizedFitting.Result; + +/** + * @author Gernot Belger + */ +final class LogLinearFittingOperation extends AbstractFittingOperation implements IFittingOperation { + + public LogLinearFittingOperation(final LogLinearAlternative function) { + super(function); + } + + @Override + public Fitting execute(final List<FittingData> data) { + + final double[] obsDischarges = new double[data.size()]; + final double[] obsWaterlevels = new double[data.size()]; + + for (int i = 0; i < data.size(); i++) { + final FittingData fittingData = data.get(i); + + obsDischarges[i] = fittingData.q; + obsWaterlevels[i] = fittingData.w; + } + + final LinearLogLinearizedFitting fitting = new LinearLogLinearizedFitting(obsDischarges, obsWaterlevels); + + final Result result = fitting.optimize(); + final double a = result.getA(); + final double b = result.getB(); + final double m = result.getM(); + final double chiSquare = result.getChiSquare(); + + final double[] parameters = new double[] { a, m, b }; + + return createFitting(parameters, chiSquare, data); + } +} \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/ArraySubstraction.java Mon Dec 02 17:56:15 2019 +0100 @@ -0,0 +1,30 @@ +/** 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; + +/** + * @author Gernot Belger + */ +final class ArraySubstraction { + + private final double[] minuend; + + public ArraySubstraction(final double[] minuend) { + this.minuend = minuend; + } + + public double[] evaluate(final double[] subtrahend) { + + final double[] difference = new double[this.minuend.length]; + for (int i = 0; i < difference.length; i++) + difference[i] = this.minuend[i] - subtrahend[i]; + return difference; + } +} \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/ChiSquare.java Mon Dec 02 17:56:15 2019 +0100 @@ -0,0 +1,52 @@ +/** 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 java.util.Arrays; + +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.stat.descriptive.AbstractUnivariateStatistic; +import org.apache.commons.math3.stat.descriptive.UnivariateStatistic; + +/** + * @author Gernot Belger + */ +final class ChiSquare extends AbstractUnivariateStatistic { + + private final double[] observation; + + public ChiSquare(final double[] observation) { + this.observation = observation; + } + + @Override + public UnivariateStatistic copy() { + /* stateless, hence we can return this */ + return this; + } + + @Override + public double evaluate(final double[] values, final int begin, final int length) throws MathIllegalArgumentException { + test(values, begin, length, true); + + final double[] residualsWeights = new double[length]; + Arrays.fill(residualsWeights, 0, length - 1, 1.0); + + double cost = 0; + + for (int i = begin; i < begin + length; i++) { + final double residual = values[i + begin] - this.observation[i + begin]; + // final double wresiduals = residual * FastMath.sqrt(residualsWeights[i]); + cost += residualsWeights[i] * residual * residual; + } + + return cost; + } +} \ No newline at end of file
--- /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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/LinearRegressionSqrtErrorFunction.java Mon Dec 02 17:56:15 2019 +0100 @@ -0,0 +1,147 @@ +/** 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.analysis.MultivariateFunction; +import org.apache.commons.math3.analysis.UnivariateFunction; +import org.apache.commons.math3.stat.descriptive.DescriptiveStatistics; +import org.apache.commons.math3.stat.regression.SimpleRegression; + +/** + * @author Gernot Belger + * + */ +final class LinearRegressionSqrtErrorFunction implements MultivariateFunction, UnivariateFunction { + + private final double[] obsDischarges; + + private final double[] obsWaterlevels; + + private final double waterlevelMean; + + private final double waterlevelDeviation; + + private final double[] standardizedWaterlevel; + + private final SQRTFunction sqrtFunction; + + private final TargetFunction targetFunction; + + public LinearRegressionSqrtErrorFunction(final double[] obsDischarges, final double[] obsWaterlevels, final TargetFunction targetFunction) { + + this.obsDischarges = obsDischarges; + this.obsWaterlevels = obsWaterlevels; + + this.sqrtFunction = new SQRTFunction(obsWaterlevels); + this.targetFunction = targetFunction; + + final DescriptiveStatistics stats = new DescriptiveStatistics(obsWaterlevels); + + // Compute mean and standard deviation + this.waterlevelMean = stats.getMean(); + this.waterlevelDeviation = stats.getStandardDeviation(); + // -mittelwert durch standardabweichung /z-transformation) + + // initialize the standardizedSample, which has the same length as the sample + this.standardizedWaterlevel = new double[obsWaterlevels.length]; + + for (int i = 0; i < obsWaterlevels.length; i++) + this.standardizedWaterlevel[i] = (obsWaterlevels[i] - this.waterlevelMean) / this.waterlevelDeviation; + } + + @Override + public double value(final double[] point) { + + final double a = point[0]; + + return value(a); + } + + @Override + public double value(final double a) { + + final double error = calc_sqrt_trans(a); + if (Double.isNaN(error)) + return Double.POSITIVE_INFINITY; + + return error; + } + + private double calc_sqrt_trans(final double a) { + + final double[] waterlevels = calc_stage_trans(a); + + return this.sqrtFunction.calc_sqrt(waterlevels); + } + + private double[] calc_stage_trans(final double a) { + + final double[] params = calc_parameters_trans(a); + + final double b = params[0]; + final double m = params[1]; + + return this.targetFunction.calc_stages(a, b, m); + } + + public double[] calc_parameters_trans(final double a) { + + // stages_obs_trans = trans_stages(stages_obs, parameter_a); + final double[] transWaterlevels = trans_stages(this.obsWaterlevels, a); + final double[] discharges = this.obsDischarges; + + // FIXME: z-nromalizing? + + // printArray(discharges); + // printArray(transWaterlevels); + + // final DescriptiveStatistics stats = new DescriptiveStatistics(transWaterlevels); + // + // // Compute mean and standard deviation + // final double mean = stats.getMean(); + // final double deviation = stats.getStandardDeviation(); + // + // final double[] normWaterlevel = new double[transWaterlevels.length]; + // for (int i = 0; i < transWaterlevels.length; i++) + // normWaterlevel[i] = (transWaterlevels[i] - mean) / deviation; + + final SimpleRegression regression = new SimpleRegression(true); + for (int i = 0; i < discharges.length; i++) + regression.addData(discharges[i], transWaterlevels[i]); + + final double normSlope = regression.getSlope(); + final double normIntercept = regression.getIntercept(); + + // unnormalize + // final double slope = normSlope * deviation; + // final double intercept = normIntercept + mean; + final double slope = normSlope; + final double intercept = normIntercept; + + return new double[] { intercept, slope }; + } + + // private void printArray(final double[] discharges) { + // for (final double d : discharges) { + // System.out.print(d); + // System.out.print(", "); + // } + // System.out.println(" "); + // } + + private double[] trans_stages(final double[] waterlevels, final double a) { + + final double[] transWaterlevels = new double[waterlevels.length]; + for (int i = 0; i < transWaterlevels.length; i++) + transWaterlevels[i] = Math.exp(waterlevels[i] / a); + + return transWaterlevels; + } +} \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/LinearizedFittingTest.java Mon Dec 02 17:56:15 2019 +0100 @@ -0,0 +1,228 @@ +/** 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 java.io.BufferedInputStream; +import java.io.IOException; +import java.io.InputStreamReader; +import java.math.BigDecimal; +import java.util.ArrayList; +import java.util.List; +import java.util.Map.Entry; +import java.util.NavigableMap; +import java.util.TreeMap; + +import org.apache.commons.lang.StringUtils; +import org.apache.commons.math3.analysis.MultivariateFunction; +import org.apache.commons.math3.optim.InitialGuess; +import org.apache.commons.math3.optim.MaxEval; +import org.apache.commons.math3.optim.MaxIter; +import org.apache.commons.math3.optim.PointValuePair; +import org.apache.commons.math3.optim.SimpleBounds; +import org.apache.commons.math3.optim.nonlinear.scalar.GoalType; +import org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunction; +import org.apache.commons.math3.optim.nonlinear.scalar.noderiv.BOBYQAOptimizer; +import org.apache.commons.math3.util.Pair; +import org.dive4elements.river.artifacts.model.fixings.fitting.LinearLogLinearizedFitting.Result; + +import au.com.bytecode.opencsv.CSVReader; + +/** + * @author Gernot Belger + */ +public class LinearizedFittingTest { + + public static void main(final String[] args) throws IOException { + + // read test data + final NavigableMap<BigDecimal, List<Pair<BigDecimal, BigDecimal>>> testData = readTestData(); + + for (final Entry<BigDecimal, List<Pair<BigDecimal, BigDecimal>>> entry : testData.entrySet()) { + + final BigDecimal station = entry.getKey(); + final List<Pair<BigDecimal, BigDecimal>> testSample = entry.getValue(); + optimizeA(station, testSample); + } + } + + private static void optimizeA(final BigDecimal station, final List<Pair<BigDecimal, BigDecimal>> testSample) { + + /* extrakt observations */ + final double[] obsDischarges = new double[testSample.size()]; + final double[] obsWaterlevels = new double[testSample.size()]; + for (int i = 0; i < obsWaterlevels.length; i++) { + obsDischarges[i] = testSample.get(i).getKey().doubleValue(); + obsWaterlevels[i] = testSample.get(i).getValue().doubleValue(); + } + + final TargetFunction targetFunction = new TargetFunction(obsDischarges); + final SqrtErrorFunction sqrtErrorFunction = new SqrtErrorFunction(obsWaterlevels, targetFunction); + + final Result directEstimation = estimateADirect(sqrtErrorFunction); + final Result linearEstimation = new LinearLogLinearizedFitting(obsDischarges, obsWaterlevels).optimize(); + + final Result directPostEstimation = postOptimize(directEstimation, sqrtErrorFunction); + final Result linearPostEstimation = postOptimize(linearEstimation, sqrtErrorFunction); + + printResult(station, testSample, directEstimation); + printResult(station, testSample, linearEstimation); + printResult(station, testSample, directPostEstimation); + printResult(station, testSample, linearPostEstimation); + } + + private static Result postOptimize(final Result estimation, final SqrtErrorFunction errorFunction) { + + final double a = estimation.getA(); + final double b = estimation.getB(); + final double m = estimation.getM(); + // final double error = estimation.getError(); + final double aIndex = estimation.getBestAindex(); + + final double aLow = Math.pow(10, aIndex - 1); + final double aHigh = Math.pow(10, aIndex + 1); + + return directOptimize(aLow, aHigh, a, b, m, errorFunction); + } + + private static void printResult(final BigDecimal station, final List<Pair<BigDecimal, BigDecimal>> testSample, final Result optimize) { + + final double a = optimize.getA(); + final double b = optimize.getB(); + final double m = optimize.getM(); + final double error = optimize.getError(); + + System.out.format("%s %.10f %.10f %.10f %.10f", station, a, b, m, error); + + for (final Pair<BigDecimal, BigDecimal> entry : testSample) { + + final BigDecimal waterlevel = entry.getSecond(); + System.out.print(' '); + System.out.print(waterlevel); + } + System.out.println(); + } + + private static Result estimateADirect(final SqrtErrorFunction sqrtErrorFunction) { + + Result best = null; + double leastError = Double.POSITIVE_INFINITY; + + // iteration über a von 10^^0 bis 10^^20 + for (int i = 0; i < 20; i++) { + + final double aStart = Math.pow(10, i); + final double aLow = Math.pow(10, i - 1); + final double aHigh = Math.pow(10, i + 1); + + final Result result = directOptimize(aLow, aHigh, aStart, 1, 1, sqrtErrorFunction); + final double error = result.getError(); + + if (error < leastError) { + leastError = error; + best = result.withBestAIndex(i); + } + } + + return best; + } + + private static Result directOptimize(final double aLow, final double aHigh, final double a, final double b, final double m, + final SqrtErrorFunction sqrtErrorFunction) { + + // n = 3 + // [n+2, (n+1)(n+2)/2] + // --> [5, 10] + final int interpolationPoints = 10; + final BOBYQAOptimizer optimizer = new BOBYQAOptimizer(interpolationPoints); + + /* optimization data */ + final MultivariateFunction function = new MultivariateFunction() { + + @Override + public double value(final double[] point) { + return sqrtErrorFunction.value(point[0], point[1], point[2]); + } + }; + + final MaxEval maxEval = new MaxEval(Integer.MAX_VALUE); + final MaxIter maxIter = new MaxIter(Integer.MAX_VALUE); + + final SimpleBounds bounds = new SimpleBounds(new double[] { aLow, -1e3, 0 }, new double[] { aHigh, 1e3, 1e3 }); + final double[] startValues = new double[] { a, b, m }; + final PointValuePair result = optimizer.optimize(GoalType.MINIMIZE, new ObjectiveFunction(function), new InitialGuess(startValues), bounds, maxEval, + maxIter); + + final Double error = result.getValue(); + final double[] point = result.getPoint(); + + final double aEstimation = point[0]; + final double bEstimation = point[1]; + final double mEstimation = point[2]; + + return new Result(aEstimation, bEstimation, mEstimation, error, -1, Double.NaN); + } + + private static NavigableMap<BigDecimal, List<Pair<BigDecimal, BigDecimal>>> readTestData() throws IOException { + + final NavigableMap<BigDecimal, List<Pair<BigDecimal, BigDecimal>>> data = new TreeMap<>(); + + try (final CSVReader reader = new CSVReader(new InputStreamReader(new BufferedInputStream(LinearizedFittingTest.class.getResourceAsStream("testdata.txt"))), + '\t')) { + + final String[] header = reader.readNext(); + if (header == null) + throw new IllegalStateException(); + + if (header.length < 2) + throw new IllegalStateException(); + + while (true) { + final String[] line = reader.readNext(); + if (line == null) + break; + + if (line.length != header.length) + throw new IllegalStateException(); + + final BigDecimal discharge = parseDecimal(line[0]); + if (discharge == null) + continue; + + for (int column = 1; column < line.length; column++) { + + final BigDecimal station = parseDecimal(header[column]); + if (station == null) + continue; + + if (!data.containsKey(station)) + data.put(station, new ArrayList<Pair<BigDecimal, BigDecimal>>()); + final List<Pair<BigDecimal, BigDecimal>> points = data.get(station); + + final BigDecimal waterlevel = parseDecimal(line[column]); + if (waterlevel != null) + points.add(Pair.create(discharge, waterlevel)); + } + } + } + + return data; + } + + private static BigDecimal parseDecimal(final String token) { + + if (StringUtils.isBlank(token)) + return null; + + if ("nan".equalsIgnoreCase(token)) + return null; + + return new BigDecimal(token); + } +} \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/SQRTFunction.java Mon Dec 02 17:56:15 2019 +0100 @@ -0,0 +1,30 @@ +/** 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.stat.StatUtils; + +/** + * @author Gernot Belger + */ +final class SQRTFunction { + private final ArraySubstraction substraction; + + public SQRTFunction(final double[] observations) { + this.substraction = new ArraySubstraction(observations); + } + + public double calc_sqrt(final double[] values) { + + final double[] difference = this.substraction.evaluate(values); + + return StatUtils.sumSq(difference); + } +} \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/SqrtErrorFunction.java Mon Dec 02 17:56:15 2019 +0100 @@ -0,0 +1,36 @@ +/** 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; + +/** + * @author Gernot Belger + */ +final class SqrtErrorFunction { + + private final SQRTFunction sqrtFunction; + + private final TargetFunction targetFunction; + + public SqrtErrorFunction(final double[] obsWaterlevels, final TargetFunction targetFunction) { + + this.targetFunction = targetFunction; + this.sqrtFunction = new SQRTFunction(obsWaterlevels); + } + + public double value(final double a, final double b, final double m) { + return calc_sqrt_trans(a, b, m); + } + + private double calc_sqrt_trans(final double a, final double b, final double m) { + + final double[] waterlevels = this.targetFunction.calc_stages(a, b, m); + return this.sqrtFunction.calc_sqrt(waterlevels); + } +} \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/TargetFunction.java Mon Dec 02 17:56:15 2019 +0100 @@ -0,0 +1,31 @@ +/** 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; + +/** + * @author Gernot Belger + */ +final class TargetFunction { + + private final double[] discharge; + + public TargetFunction(final double[] discharge) { + this.discharge = discharge; + } + + public double[] calc_stages(final double a, final double b, final double m) { + + final double[] waterlevels = new double[this.discharge.length]; + for (int i = 0; i < this.discharge.length; i++) + waterlevels[i] = a * Math.log(b + m * this.discharge[i]); + + return waterlevels; + } +} \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/fitting/testdata.txt Mon Dec 02 17:56:15 2019 +0100 @@ -0,0 +1,7 @@ +Q 541 541.1 541.2 541.3 541.4 541.5 541.6 541.7 541.8 541.9 542 542.1 542.2 542.3 542.4 542.5 542.6 542.7 542.8 542.9 543 +750.0 71.48 71.425 71.37 71.305 71.24 71.18 71.12 71.07 71.02 70.94 70.86 70.725 70.59 70.53 70.39 70.35 70.31 70.29 70.27 70.235 70.2 +750.0 71.46 71.405 71.35 71.295 nan nan nan 71.035 70.95 70.86 70.77 70.675 70.58 70.49 70.47 70.385 70.3 70.24 70.18 70.145 70.11 +1555.0 72.31 72.24 72.17 72.09 72.01 71.93 71.85 71.74 71.63 71.545 71.46 71.395 71.33 71.28 71.22 71.18 71.14 71.12 71.1 71.035 70.97 +1583.0 72.24 72.178 72.116 72.054 71.992 71.93 71.852 71.774 71.696 71.618 71.54 71.474 71.408 71.342 71.276 71.21 71.176 71.142 71.108 71.07 71.04 +1618.0 72.38 72.305 72.23 72.145 72.06 71.98 71.90 71.805 71.71 71.61 71.50 71.42 71.34 71.295 71.25 71.23 71.21 71.185 71.16 71.11 71.06 +5682.0 75.99 75.97 75.96 75.94 75.93 75.91 75.88 75.86 75.83 75.81 75.78 75.76 75.73 75.71 75.68 75.66 75.62 75.59 75.55 75.52 75.48
--- a/gwt-client/src/main/java/org/dive4elements/river/client/client/ui/fixation/FixFunctionSelect.java Mon Dec 02 17:21:43 2019 +0100 +++ b/gwt-client/src/main/java/org/dive4elements/river/client/client/ui/fixation/FixFunctionSelect.java Mon Dec 02 17:56:15 2019 +0100 @@ -33,7 +33,7 @@ funcDesc.put("log", "W(Q) = m*ln(Q + b)"); funcDesc.put("linear", "W(Q) = m * Q + b"); funcDesc.put("log-linear", "W(Q) = a*ln(m*Q+b)"); - funcDesc.put("log-linear-alternative", "W(Q) = a*ln(m*Q+b) (linearisiert)"); + funcDesc.put("log-linear-linearisiert", "W(Q) = a*ln(m*Q+b) (linearisiert)"); funcDesc.put("exp", "W(Q) = m * a^Q + b"); funcDesc.put("quad", "W(Q) = n*Q^2+m*Q+b"); funcDesc.put("pow", "W(Q) = a * Q^c + d");