teichmann@5863: /* Copyright (C) 2011, 2012, 2013 by Bundesanstalt für Gewässerkunde teichmann@5863: * Software engineering by Intevation GmbH teichmann@5863: * teichmann@5994: * This file is Free Software under the GNU AGPL (>=v3) teichmann@5863: * and comes with ABSOLUTELY NO WARRANTY! Check out the teichmann@5994: * documentation coming with Dive4Elements River for details. teichmann@5863: */ teichmann@5863: teichmann@5831: package org.dive4elements.river.artifacts.model.fixings; sascha@3011: sascha@3107: import java.util.ArrayList; gernotbelger@9415: import java.util.Date; sascha@3107: import java.util.List; sascha@3011: sascha@3011: import org.apache.commons.math.MathException; sascha@3011: import org.apache.commons.math.optimization.fitting.CurveFitter; sascha@3011: import org.apache.commons.math.optimization.general.LevenbergMarquardtOptimizer; sascha@3107: import org.apache.commons.math.stat.descriptive.moment.StandardDeviation; rrenkert@4794: import org.apache.log4j.Logger; teichmann@5831: import org.dive4elements.river.artifacts.math.GrubbsOutlier; teichmann@5831: import org.dive4elements.river.artifacts.math.fitting.Function; sascha@3011: gernotbelger@9348: public class Fitting { sascha@3011: private static Logger log = Logger.getLogger(Fitting.class); sascha@3011: d@9627: private static double[] previousParameters = null; // new double[] { Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE }; Dann schlägt die d@9627: // Kurvenanpassung fehl! Gernot hat recht! d@9627: sascha@3011: /** Use instance of this factory to find meta infos for outliers. */ sascha@3096: public interface QWDFactory { gernotbelger@9360: QWD create(double q, double w, double deltaW, boolean isOutlier); sascha@3011: } sascha@3011: gernotbelger@9415: private static final class FittingData { gernotbelger@9415: public final FixingColumnWithData event; gernotbelger@9415: public final double w; gernotbelger@9415: public final double q; gernotbelger@9415: public final boolean isInterpolated; gernotbelger@9415: gernotbelger@9415: public FittingData(final FixingColumnWithData event, final double w, final double q, final boolean isInterpolated) { gernotbelger@9415: this.event = event; gernotbelger@9415: this.w = w; gernotbelger@9415: this.q = q; gernotbelger@9415: this.isInterpolated = isInterpolated; gernotbelger@9415: } gernotbelger@9415: } gernotbelger@9415: gernotbelger@9360: private final double chiSqr; sascha@3022: gernotbelger@9360: private final double[] parameters; sascha@3011: gernotbelger@9360: private final double standardDeviation; sascha@3011: gernotbelger@9415: private final double maxQ; gernotbelger@9360: gernotbelger@9415: public Fitting(final double[] parameters, final double standardDeviation, final double chiSqr, final double maxQ) { gernotbelger@9360: this.parameters = parameters; gernotbelger@9360: this.standardDeviation = standardDeviation; gernotbelger@9360: this.chiSqr = chiSqr; gernotbelger@9415: this.maxQ = maxQ; sascha@3022: } sascha@3022: sascha@3011: public double getChiSquare() { gernotbelger@9415: return this.chiSqr; sascha@3022: } sascha@3022: sascha@3065: public double getMaxQ() { gernotbelger@9415: return this.maxQ; sascha@3065: } sascha@3065: gernotbelger@9348: public double[] getParameters() { gernotbelger@9415: return this.parameters; sascha@3011: } sascha@3011: sascha@3107: public double getStandardDeviation() { gernotbelger@9415: return this.standardDeviation; sascha@3107: } sascha@3107: gernotbelger@9415: public static Fitting fit(final FixResultColumns resultColumns, final double km, final Function function, final boolean checkOutliers, gernotbelger@9415: final List eventColumns) { sascha@3011: gernotbelger@9415: final int numEvents = eventColumns.size(); gernotbelger@9415: final List data = new ArrayList<>(numEvents); sascha@3011: gernotbelger@9415: final double[] wTemp = new double[1]; gernotbelger@9415: gernotbelger@9415: for (final FixingColumnWithData event : eventColumns) { gernotbelger@9415: gernotbelger@9415: final double q = event.getQ(km); gernotbelger@9415: if (!Double.isNaN(q)) { gernotbelger@9415: final boolean isInterpolated = !event.getW(km, wTemp, 0); gernotbelger@9415: final double w = wTemp[0]; gernotbelger@9415: gernotbelger@9415: if (!Double.isNaN(w)) gernotbelger@9415: data.add(new FittingData(event, w, q, isInterpolated)); sascha@3011: } sascha@3011: } sascha@3011: gernotbelger@9415: if (data.size() < 2) { gernotbelger@9415: log.warn("Not enough data for fitting."); gernotbelger@9360: return null; sascha@3011: } sascha@3011: gernotbelger@9415: final List outliers = new ArrayList<>(data.size()); sascha@3011: teichmann@5831: org.dive4elements.river.artifacts.math.Function instance = null; sascha@3202: LevenbergMarquardtOptimizer lmo = null; gernotbelger@9360: double[] parameters = null; sascha@3202: sascha@3011: for (;;) { gernotbelger@9415: sascha@3202: parameters = null; gernotbelger@9415: aheinecke@7300: for (double tolerance = 1e-10; tolerance < 1e-1; tolerance *= 10d) { sascha@3011: sascha@3202: lmo = new LevenbergMarquardtOptimizer(); sascha@3202: lmo.setCostRelativeTolerance(tolerance); sascha@3202: lmo.setOrthoTolerance(tolerance); sascha@3202: lmo.setParRelativeTolerance(tolerance); sascha@3202: gernotbelger@9415: try { gernotbelger@9415: final CurveFitter cf = new CurveFitter(lmo); sascha@3202: gernotbelger@9415: for (final FittingData fittingData : data) gernotbelger@9415: cf.addObservedPoint(fittingData.q, fittingData.w); sascha@3202: d@9627: parameters = cf.fit(function, // previousParameters != null ? previousParameters : d@9627: function.getInitialGuess()); d@9627: sascha@3202: break; sascha@3202: } gernotbelger@9415: catch (final MathException me) { d@9627: sascha@3202: if (log.isDebugEnabled()) { gernotbelger@9360: log.debug("tolerance " + tolerance + " + failed.", me); sascha@3202: } sascha@3202: } sascha@3011: } gernotbelger@9360: sascha@3202: if (parameters == null) { aheinecke@7300: /* gernotbelger@9348: * log.debug("Parameters is null"); gernotbelger@9348: * for (int i = 0, N = xs.size(); i < N; ++i) { gernotbelger@9348: * log.debug("DATA: " + xs.getQuick(i) + " " + ys.getQuick(i)); gernotbelger@9348: * } gernotbelger@9348: */ gernotbelger@9360: return null; sascha@3011: } d@9627: previousParameters = parameters; gernotbelger@9415: // This is the parameterized function for a given km. raimund@3127: instance = function.instantiate(parameters); raimund@3127: gernotbelger@9360: if (!checkOutliers) sascha@3011: break; sascha@3011: gernotbelger@9415: /* find the outlier */ gernotbelger@9415: final List inputs = new ArrayList<>(data.size()); gernotbelger@9415: for (final FittingData fittingData : data) { sascha@3011: gernotbelger@9415: double y = instance.value(fittingData.q); gernotbelger@9415: if (Double.isNaN(y)) sascha@3565: y = Double.MAX_VALUE; gernotbelger@9415: gernotbelger@9415: inputs.add(Double.valueOf(fittingData.w - y)); sascha@3011: } sascha@3011: gernotbelger@9360: final Integer outlier = GrubbsOutlier.findOutlier(inputs); gernotbelger@9360: if (outlier == null) sascha@3011: break; sascha@3011: gernotbelger@9360: final int idx = outlier.intValue(); ingo@3066: gernotbelger@9415: // outliers.add(qwdFactory.create(xs.getQuick(idx), ys.getQuick(idx), Double.NaN, true)); gernotbelger@9415: final FittingData removed = data.remove(idx); gernotbelger@9415: outliers.add(removed); sascha@3022: } sascha@3022: gernotbelger@9415: /* now build result data */ gernotbelger@9415: final List qwds = new ArrayList<>(data.size()); gernotbelger@9360: gernotbelger@9415: /* calculate dW of outliers against the resulting function and add them to results */ gernotbelger@9415: for (final FittingData outlier : outliers) { gernotbelger@9415: final QWD qwd = createQWD(outlier, instance, true); gernotbelger@9360: qwds.add(qwd); gernotbelger@9415: resultColumns.addQWD(outlier.event, km, qwd); gernotbelger@9415: } gernotbelger@9415: gernotbelger@9415: /* gernotbelger@9415: * calculate dW of used values against the resulting function and add them to results , also calculate standard * gernotbelger@9415: * deviation gernotbelger@9415: */ gernotbelger@9415: final StandardDeviation stdDev = new StandardDeviation(); gernotbelger@9415: double maxQ = -Double.MAX_VALUE; gernotbelger@9415: gernotbelger@9415: for (final FittingData fittingData : data) { gernotbelger@9415: gernotbelger@9415: final QWD qwd = createQWD(fittingData, instance, false); gernotbelger@9415: qwds.add(qwd); gernotbelger@9415: resultColumns.addQWD(fittingData.event, km, qwd); gernotbelger@9415: gernotbelger@9415: stdDev.increment(qwd.getDeltaW()); gernotbelger@9415: gernotbelger@9415: final double q = qwd.getQ(); gernotbelger@9415: if (q > maxQ) gernotbelger@9415: maxQ = q; gernotbelger@9360: } gernotbelger@9360: gernotbelger@9360: final double standardDeviation = stdDev.getResult(); gernotbelger@9360: gernotbelger@9360: final double chiSqr = lmo.getChiSquare(); gernotbelger@9360: gernotbelger@9415: return new Fitting(parameters, standardDeviation, chiSqr, maxQ); gernotbelger@9415: } gernotbelger@9415: gernotbelger@9415: private static QWD createQWD(final FittingData data, final org.dive4elements.river.artifacts.math.Function function, final boolean isOutlier) { gernotbelger@9415: gernotbelger@9415: final FixingColumnWithData event = data.event; gernotbelger@9415: final Date date = event.getDate(); gernotbelger@9415: final boolean isInterpolated = data.isInterpolated; gernotbelger@9415: gernotbelger@9415: final double w = data.w; gernotbelger@9415: final double q = data.q; gernotbelger@9415: final double dw = (w - function.value(q)) * 100.0; gernotbelger@9415: gernotbelger@9415: return new QWD(q, w, date, isInterpolated, dw, isOutlier); sascha@3011: } gernotbelger@9348: }