package de.intevation.flys.artifacts.math;

import org.apache.commons.math.MathException;

import org.apache.commons.math.stat.descriptive.moment.Mean;
import org.apache.commons.math.stat.descriptive.moment.StandardDeviation;

import org.apache.commons.math.distribution.TDistributionImpl;

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

import org.apache.log4j.Logger;

public class Outlier
{
    public static final double DEFAULT_ALPHA = 0.05;

    private static Logger log = Logger.getLogger(Outlier.class);

    public static class IndexedValue
    implements          Comparable<IndexedValue>
    {
        protected int    index;
        protected double value;

        public IndexedValue() {
        }

        public IndexedValue(int index, double value) {
            this.index = index;
            this.value = value;
        }

        public int getIndex() {
            return index;
        }

        public void setIndex(int index) {
            this.index = index;
        }

        public double getValue() {
            return value;
        }

        public void setValue(double value) {
            this.value = value;
        }

        @Override
        public int compareTo(IndexedValue other) {
            int diff = index - other.index;
            if (index < 0) return -1;
            return index > 0 ? +1 : 0;
        }
    } // class IndexedValue

    public static class Outliers {

        protected List<IndexedValue> retained;
        protected List<IndexedValue> removed;

        public Outliers() {
        }

        public Outliers(
            List<IndexedValue> retained,
            List<IndexedValue> removed
        ) {
            this.retained = retained;
            this.removed  = removed;
        }

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

        public List<IndexedValue> getRetained() {
            return retained;
        }

        public void setRetained(List<IndexedValue> retained) {
            this.retained = retained;
        }

        public List<IndexedValue> getRemoved() {
            return removed;
        }

        public void setRemoved(List<IndexedValue> removed) {
            this.removed = removed;
        }
    } // class Outliers

    public Outlier() {
    }

    public static Outliers findOutliers(List<IndexedValue> inputValues) {
        return findOutliers(inputValues, DEFAULT_ALPHA);
    }

    public static Outliers findOutliers(
        List<IndexedValue> inputValues,
        double             alpha
    ) {
        ArrayList<IndexedValue> outliers = new ArrayList<IndexedValue>();

        ArrayList<IndexedValue> values =
            new ArrayList<IndexedValue>(inputValues);

        for (;;) {
            int N = values.size();

            if (N < 4) {
                break;
            }

            Mean mean = new Mean();
            StandardDeviation std = new StandardDeviation();

            for (IndexedValue value: values) {
                mean.increment(value.getValue());
                std .increment(value.getValue());
            }

            double m = mean.getResult();
            double s = std.getResult();

            double maxZ = -Double.MAX_VALUE;
            int iv = -1;
            for (int i = N-1; i >= 0; --i) {
                IndexedValue v = values.get(i);
                double z = Math.abs(m - v.getValue())/s;
                if (z > maxZ) {
                    maxZ = z;
                    iv = i;
                }
            }

            double t = Math.sqrt((N*(N-2)*maxZ*maxZ)
                /((N-1)*(N-1) - N*maxZ*maxZ));

            TDistributionImpl tdist = new TDistributionImpl(N-2);

            try {
                double p = tdist.cumulativeProbability(t);

                if (p < alpha) {
                    outliers.add(values.get(iv));
                    values.remove(iv);
                }
                else {
                    break;
                }
            }
            catch (MathException me) {
                log.error(me);
            }
        }

        Collections.sort(outliers);

        return new Outliers(values, outliers);
    }
}
// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :