view flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/sq/Fitting.java @ 4187:21f4e4b79121

Refactor GaugeDischargeCurveFacet to be able to set a facet name For adding another output of the GaugeDischargeCurveArtifact it is necessary to provide to facet instances with different names. Therefore the GaugeDischargeCurveFacet is extended to set the facet name in the constructor.
author Björn Ricks <bjoern.ricks@intevation.de>
date Fri, 19 Oct 2012 13:25:49 +0200
parents b8b1280606c2
children aaf810d4ec82
line wrap: on
line source
package de.intevation.flys.artifacts.model.sq;

import de.intevation.flys.artifacts.math.fitting.Function;

import java.util.ArrayList;
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;

public class Fitting
implements   Outlier.Callback
{
    private static Logger log = Logger.getLogger(Fitting.class);

    public interface Callback {

        void afterIteration(
            double [] parameters,
            SQ []     measurements,
            SQ []     outliers,
            double    standardDeviation,
            double    chiSqr);
    } // interfacte

    protected Function function;

    protected double [] coeffs;

    protected de.intevation.flys.artifacts.math.Function instance;

    protected double stdDevFactor;
    protected double chiSqr;

    protected Callback callback;

    public Fitting() {
    }

    public Fitting(Function function, double stdDevFactor) {
        this();
        this.function     = function;
        this.stdDevFactor = stdDevFactor;
    }

    public Function getFunction() {
        return function;
    }

    public void setFunction(Function function) {
        this.function = function;
    }

    public double getStdDevFactor() {
        return stdDevFactor;
    }

    public void setStdDevFactor(double stdDevFactor) {
        this.stdDevFactor = stdDevFactor;
    }

    @Override
    public void initialize(List<SQ> sqs) throws MathException {

        LevenbergMarquardtOptimizer lmo =
            new LevenbergMarquardtOptimizer();

        CurveFitter cf = new CurveFitter(lmo);
        for (SQ sq: sqs) {
            cf.addObservedPoint(sq.getQ(), sq.getS());
        }

        coeffs = cf.fit(
            function, function.getInitialGuess());

        instance = function.instantiate(coeffs);

        chiSqr = lmo.getChiSquare();
    }

    @Override
    public double eval(SQ sq) {
        double s = instance.value(sq.q);
        return sq.s - s;
    }

    @Override
    public void iterationFinished(
        double   standardDeviation,
        SQ       outlier,
        List<SQ> remainings
    ) {
        if (log.isDebugEnabled()) {
            log.debug("iterationFinished ----");
            log.debug(" num remainings: " + remainings.size());
            log.debug(" has outlier: " + outlier != null);
            log.debug(" standardDeviation: " + standardDeviation);
            log.debug(" Chi^2: " + chiSqr);
            log.debug("---- iterationFinished");
        }
        callback.afterIteration(
            coeffs,
            remainings.toArray(new SQ[remainings.size()]),
            outlier != null ? new SQ [] { outlier} : new SQ [] {},
            standardDeviation,
            chiSqr);
    }

    protected static final List<SQ> onlyValid(List<SQ> sqs) {

        List<SQ> good = new ArrayList<SQ>(sqs.size());

        for (SQ sq: sqs) {
            if (sq.isValid()) {
                good.add(sq);
            }
        }

        return good;
    }

    public boolean fit(List<SQ> sqs, Callback callback) {

        sqs = onlyValid(sqs);

        if (sqs.size() < 2) {
            log.warn("Too less points for fitting.");
            return false;
        }

        this.callback = callback;

        try {
            Outlier.detectOutliers(this, sqs, stdDevFactor);
        }
        catch (MathException me) {
            log.warn(me);
            return false;
        }

        return true;
    }
}
// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :

http://dive4elements.wald.intevation.org