sascha@3188: package de.intevation.flys.artifacts.model.sq;
sascha@3188: 
sascha@3188: import de.intevation.flys.artifacts.math.fitting.Function;
sascha@3188: 
sascha@3188: import java.util.ArrayList;
sascha@3188: import java.util.Iterator;
sascha@3188: import java.util.List;
sascha@3188: 
sascha@3188: import org.apache.commons.math.MathException;
sascha@3188: 
sascha@3188: import org.apache.commons.math.optimization.fitting.CurveFitter;
sascha@3188: 
sascha@3188: import org.apache.commons.math.optimization.general.LevenbergMarquardtOptimizer;
sascha@3188: 
sascha@3188: import org.apache.log4j.Logger;
sascha@3188: 
sascha@3188: public class Fitting
sascha@3188: {
sascha@3188:     private static Logger log = Logger.getLogger(Fitting.class);
sascha@3188: 
sascha@3188:     protected Function function;
sascha@3188: 
sascha@3188:     protected double [] parameters;
sascha@3188: 
sascha@3188:     protected double stdDevFactor;
sascha@3188: 
sascha@3188:     protected double standardDeviation;
sascha@3188: 
sascha@3188:     protected double chiSqr;
sascha@3188: 
sascha@3190:     protected SQ [] remaining;
sascha@3188: 
sascha@3190:     protected List<SQ []> outliers;
sascha@3188: 
sascha@3188:     public Fitting() {
sascha@3188:     }
sascha@3188: 
sascha@3188:     public Fitting(Function function, double stdDevFactor) {
sascha@3188:         this.function     = function;
sascha@3188:         this.stdDevFactor = stdDevFactor;
sascha@3188:     }
sascha@3188: 
sascha@3188:     public Function getFunction() {
sascha@3188:         return function;
sascha@3188:     }
sascha@3188: 
sascha@3188:     public void setFunction(Function function) {
sascha@3188:         this.function = function;
sascha@3188:     }
sascha@3188: 
sascha@3188:     public double [] getParameters() {
sascha@3188:         return parameters;
sascha@3188:     }
sascha@3188: 
sascha@3188:     public void setParameters(double [] parameters) {
sascha@3188:         this.parameters = parameters;
sascha@3188:     }
sascha@3188: 
sascha@3188:     public double getStdDevFactor() {
sascha@3188:         return stdDevFactor;
sascha@3188:     }
sascha@3188: 
sascha@3188:     public void setStdDevFactor(double stdDevFactor) {
sascha@3188:         this.stdDevFactor = stdDevFactor;
sascha@3188:     }
sascha@3188: 
sascha@3188:     public double getStandardDeviation() {
sascha@3188:         return standardDeviation;
sascha@3188:     }
sascha@3188: 
sascha@3188:     public void setStandardDeviation(double standardDeviation) {
sascha@3188:         this.standardDeviation = standardDeviation;
sascha@3188:     }
sascha@3188: 
sascha@3188:     public double getChiSqr() {
sascha@3188:         return chiSqr;
sascha@3188:     }
sascha@3188: 
sascha@3188:     public void setChiSqr(double chiSqr) {
sascha@3188:         this.chiSqr = chiSqr;
sascha@3188:     }
sascha@3188: 
sascha@3190:     public SQ [] getRemaining() {
sascha@3188:         return remaining;
sascha@3188:     }
sascha@3188: 
sascha@3190:     public void setRemaining(SQ [] remaining) {
sascha@3188:         this.remaining = remaining;
sascha@3188:     }
sascha@3188: 
sascha@3190:     public List<SQ []> getOutliers() {
sascha@3188:         return outliers;
sascha@3188:     }
sascha@3188: 
sascha@3190:     public void setOutliers(List<SQ []> outliers) {
sascha@3188:         this.outliers = outliers;
sascha@3188:     }
sascha@3188: 
sascha@3188:     public void reset() {
sascha@3188:         outliers          = null;
sascha@3188:         remaining         = null;
sascha@3188:         parameters        = null;
sascha@3188:         standardDeviation = 0d;
sascha@3188:         standardDeviation = 0d;
sascha@3188:         chiSqr            = 0d;
sascha@3188:     }
sascha@3188: 
sascha@3188:     protected static final List<SQ> onlyValid(List<SQ> sqs) {
sascha@3188: 
sascha@3188:         List<SQ> good = new ArrayList<SQ>(sqs.size());
sascha@3188: 
sascha@3188:         for (SQ sq: sqs) {
sascha@3188:             if (sq.isValid()) {
sascha@3188:                 good.add(sq);
sascha@3188:             }
sascha@3188:         }
sascha@3188: 
sascha@3188:         return good;
sascha@3188:     }
sascha@3188: 
sascha@3188:     public boolean fit(List<SQ> sqs) {
sascha@3188: 
sascha@3188:         sqs = onlyValid(sqs);
sascha@3188: 
sascha@3188:         if (sqs.size() < 2) {
sascha@3188:             log.warn("Too less points for fitting.");
sascha@3188:             return false;
sascha@3188:         }
sascha@3188: 
sascha@3188:         final LevenbergMarquardtOptimizer lmo =
sascha@3188:             new LevenbergMarquardtOptimizer();
sascha@3188: 
sascha@3188:         CurveFitter cf = new CurveFitter(lmo);
sascha@3188: 
sascha@3188:         for (SQ sq: sqs) {
sascha@3188:             cf.addObservedPoint(sq.getQ(), sq.getS());
sascha@3188:         }
sascha@3188: 
sascha@3188:         try {
sascha@3188:             parameters = cf.fit(function, function.getInitialGuess());
sascha@3188:         }
sascha@3188:         catch (MathException me) {
sascha@3188:             log.warn(me);
sascha@3188:             return false;
sascha@3188:         }
sascha@3188: 
sascha@3188:         chiSqr = lmo.getChiSquare();
sascha@3188: 
sascha@3188:         final de.intevation.flys.artifacts.math.Function [] instance = {
sascha@3188:             function.instantiate(parameters)
sascha@3188:         };
sascha@3188: 
sascha@3188:         try {
sascha@3188:             remaining = Outlier.detectOutliers(
sascha@3188:                 new Outlier.Callback() {
sascha@3188: 
sascha@3188:                     List<List<SQ>> outliers =
sascha@3188:                         new ArrayList<List<SQ>>();
sascha@3188: 
sascha@3188:                     int currentIteration;
sascha@3188: 
sascha@3188:                     @Override
sascha@3188:                     public double eval(SQ sq) {
sascha@3188:                         double s = instance[0].value(sq.q);
sascha@3188:                         return s - sq.s;
sascha@3188:                     }
sascha@3188: 
sascha@3188:                     @Override
sascha@3188:                     public void iteration(int i) {
sascha@3188:                         currentIteration = i;
sascha@3188:                     }
sascha@3188: 
sascha@3188:                     @Override
sascha@3188:                     public void outlier(SQ sq) {
sascha@3188:                         if (currentIteration > outliers.size()) {
sascha@3188:                             outliers.add(new ArrayList<SQ>(2));
sascha@3188:                         }
sascha@3188:                         outliers.get(currentIteration-1).add(sq);
sascha@3188:                     }
sascha@3188: 
sascha@3188:                     @Override
sascha@3188:                     public void standardDeviation(double stdDev) {
sascha@3188:                         setStandardDeviation(stdDev);
sascha@3188:                     }
sascha@3188: 
sascha@3188:                     @Override
sascha@3188:                     public void reinitialize(Iterator<SQ> good)
sascha@3188:                     throws MathException
sascha@3188:                     {
sascha@3188:                         CurveFitter cf = new CurveFitter(lmo);
sascha@3188:                         while (good.hasNext()) {
sascha@3188:                             SQ sq = good.next();
sascha@3188:                             cf.addObservedPoint(sq.getQ(), sq.getS());
sascha@3188:                         }
sascha@3188: 
sascha@3188:                         parameters = cf.fit(
sascha@3188:                             function, function.getInitialGuess());
sascha@3188: 
sascha@3188:                         instance[0] = function.instantiate(parameters);
sascha@3188: 
sascha@3188:                         chiSqr = lmo.getChiSquare();
sascha@3188:                     }
sascha@3188: 
sascha@3188:                     @Override
sascha@3188:                     public void finished() {
sascha@3190:                         List<SQ []> result =
sascha@3190:                             new ArrayList<SQ []>(outliers.size());
sascha@3190: 
sascha@3190:                         for (List<SQ> ols: outliers) {
sascha@3190:                             result.add(ols.toArray(new SQ[ols.size()]));
sascha@3190:                         }
sascha@3190: 
sascha@3190:                         setOutliers(result);
sascha@3188:                     }
sascha@3188:                 },
sascha@3188:                 sqs,
sascha@3188:                 stdDevFactor);
sascha@3188:         }
sascha@3188:         catch (MathException me) {
sascha@3188:             log.warn(me);
sascha@3188:             return false;
sascha@3188:         }
sascha@3188: 
sascha@3188:         return true;
sascha@3188:     }
sascha@3188: }
sascha@3188: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :