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

import org.apache.log4j.Logger;

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

    public static class 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;
        }
    } // class IndexedValue

    public Outlier() {
    }

    public static List<IndexedValue> 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);
                }
            }
            catch (MathException me) {
                log.error(me);
            }
        }


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