package de.intevation.flys.artifacts.model.fixings;

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

import de.intevation.flys.artifacts.math.Outlier;

import de.intevation.flys.artifacts.math.Outlier.IndexedValue;
import de.intevation.flys.artifacts.math.Outlier.Outliers;

import org.apache.commons.math.MathException;

import org.apache.commons.math.optimization.fitting.CurveFitter;

import org.apache.commons.math.optimization.general.LevenbergMarquardtOptimizer;

import gnu.trove.TDoubleArrayList;

import org.apache.log4j.Logger;

import java.util.ArrayList;
import java.util.List;

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

    /** Use instance of this factory to find meta infos for outliers. */
    public interface QWFactory {

        QW create(double q, double w);

    } // interface QWFactory

    protected Function      function;
    protected QWFactory     qwFactory;
    protected double        chiSqr;
    protected double []     parameters;
    protected ArrayList<QW> removed;


    public Fitting() {
        removed = new ArrayList<QW>();
    }

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

    public Fitting(Function function, QWFactory qwFactory) {
        this(function);
        this.qwFactory = qwFactory;
    }

    public Function getFunction() {
        return function;
    }

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

    public double getChiSquare() {
        return chiSqr;
    }

    public void reset() {
        chiSqr     = 0.0;
        parameters = null;
        removed.clear();
    }

    public boolean hasOutliers() {
        return !removed.isEmpty();
    }

    public List<QW> getOutliers() {
        return removed;
    }

    public QW [] outliersToArray() {
        return removed.toArray(new QW[removed.size()]);
    }

    public double [] getParameters() {
        return parameters;
    }

    public boolean fit(double [] qs, double [] ws) {

        TDoubleArrayList xs = new TDoubleArrayList(qs.length);
        TDoubleArrayList ys = new TDoubleArrayList(ws.length);

        for (int i = 0; i < qs.length; ++i) {
            if (!Double.isNaN(qs[i]) && !Double.isNaN(ws[i])) {
                xs.add(qs[i]);
                ys.add(ws[i]);
            }
            else {
                log.warn("remove invalid value " + qs[i] + " " + ws[i]);
            }
        }

        if (xs.size() < 2) {
            return false;
        }

        LevenbergMarquardtOptimizer lmo = new LevenbergMarquardtOptimizer();

        double [] parameters;

        List<IndexedValue> inputs = new ArrayList<IndexedValue>(xs.size());

        for (;;) {
            CurveFitter cf = new CurveFitter(lmo);

            for (int i = 0, N = xs.size(); i < N; ++i) {
                cf.addObservedPoint(xs.getQuick(i), ys.getQuick(i));
            }

            try {
                parameters = cf.fit(function, function.getInitialGuess());
            }
            catch (MathException me) {
                log.warn(me);
                return false;
            }

            if (qwFactory == null) {
                break;
            }

            inputs.clear();

            // This is the paraterized function for a given km.
            de.intevation.flys.artifacts.math.Function instance =
                function.instantiate(parameters);

            for (int i = 0, N = xs.size(); i < N; ++i) {
                double y = instance.value(xs.getQuick(i));
                if (Double.isNaN(y)) {
                    continue;
                }
                inputs.add(new IndexedValue(i, ys.getQuick(i) - y));
            }

            Outliers outliers = Outlier.findOutliers(inputs);

            if (!outliers.hasOutliers()) {
                break;
            }

            List<IndexedValue> rem = outliers.getRemoved();

            for (int i = rem.size()-1; i >= 0; --i) {
                int idx = rem.get(i).getIndex();
                removed.add(
                    qwFactory.create(
                        xs.getQuick(idx), ys.getQuick(idx)));
                xs.remove(idx);
                ys.remove(idx);
            }
        }

        chiSqr = lmo.getChiSquare();

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