/* Copyright (C) 2011, 2012, 2013 by Bundesanstalt für Gewässerkunde
 * Software engineering by Intevation GmbH
 *
 * This file is Free Software under the GNU AGPL (>=v3)
 * and comes with ABSOLUTELY NO WARRANTY! Check out the
 * documentation coming with Dive4Elements River for details.
 */

package org.dive4elements.river.artifacts.model;

import gnu.trove.TDoubleArrayList;

import java.io.Serializable;

import java.util.ArrayList;
import java.util.List;
import java.util.Set;
import java.util.TreeSet;

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

import org.dive4elements.river.artifacts.math.Linear;

import org.dive4elements.river.utils.EpsilonComparator;

public class LinearInterpolated
implements   Serializable
{
    public static final double EPSILON = 1e-5;
    public static final double MULTIPLE_STD_DEV = 4d;

    public static final EpsilonComparator CMP = new EpsilonComparator(EPSILON);

    private TDoubleArrayList xs;
    private TDoubleArrayList ys;

    private List<Range> gaps;

    public interface Operator {
        double evaluate(double y1, double y2);
    } // interface Operator

    public static final Operator SUB = new Operator() {
        @Override
        public double evaluate(double y1, double y2) {
            return y1 - y2;
        }
    };

    public static final Operator MAX = new Operator() {
        @Override
        public double evaluate(double y1, double y2) {
            return Math.max(y1, y2);
        }
    };

    public LinearInterpolated() {
        xs = new TDoubleArrayList();
        ys = new TDoubleArrayList();
    }

    public LinearInterpolated(int capacity) {
        xs = new TDoubleArrayList(capacity);
        ys = new TDoubleArrayList(capacity);
    }

    public void add(double x, double y) {
        xs.add(x);
        ys.add(y);
    }

    public int size() {
        return xs.size();
    }

    public void pointsInRange(double from, double to, Set<Double> points) {
        if (from > to) {
            double t = from;
            from = to;
            to = t;
        }
        for (int i = 0, S = size(); i < S; ++i) {
            double x = xs.getQuick(i);
            if (x >= from && x <= to) {
                points.add(x);
            }
        }
    }

    public boolean inGap(double x) {
        if (gaps != null) {
            for (Range gap: gaps) {
                if (gap.inside(x)) {
                    return true;
                }
            }
        }
        return false;
    }

    public void detectGaps(double threshold) {
        List<Range> gabs = new ArrayList<Range>();
        for (int i = 1, S = size(); i < S; ++i) {
            double x0 = xs.getQuick(i-1);
            double x1 = xs.getQuick(i);
            double minX, maxX;
            if (x0 < x1) { minX = x0; maxX = x1; }
            else         { minX = x1; maxX = x0; }
            double diff = maxX - minX - 2d*EPSILON;
            if (diff > threshold) {
                gaps.add(new Range(minX+EPSILON, maxX-EPSILON));
            }
        }
        this.gaps = gaps.isEmpty() ? null : gabs;
    }

    public void resetGaps() {
        gaps = null;
    }

    public double guessGapThreshold() {
        return guessGapThreshold(MULTIPLE_STD_DEV);
    }

    public double guessGapThreshold(double scale) {
        int S = size();
        if (S < 5) { // Too less points.
            return Double.MAX_VALUE;
        }

        StandardDeviation s = new StandardDeviation();

        for (int i = 1; i < S; ++i) {
            double diff = Math.abs(xs.getQuick(i-1) - xs.getQuick(i));
            s.increment(diff);
        }

        return scale*s.getResult();
    }

    public double value(double x) {
        for (int i = 0, S = size(); i < S; ++i) {
            double x1 = xs.getQuick(i);
            if (Math.abs(x1 - x) < EPSILON) {
                return ys.getQuick(i);
            }
            if (i > 0) {
                double x0 = xs.getQuick(i-1);
                double minX, maxX;
                if (x0 < x1) { minX = x0; maxX = x1; }
                else         { minX = x1; maxX = x0; }
                if (x > minX && x < maxX) {
                    return Linear.linear(
                        x,
                        x0, x1,
                        ys.getQuick(i-1), ys.getQuick(i));
                }
            }
        }
        return Double.NaN;
    }

    public LinearInterpolated sub(
        LinearInterpolated other, 
        double             from,
        double             to
    ) {
        return apply(SUB, other, from, to);
    }

    public LinearInterpolated max(
        LinearInterpolated other, 
        double             from,
        double             to
    ) {
        return apply(MAX, other, from, to);
    }

    public boolean intersect(LinearInterpolated other) {
        if (xs.isEmpty() || other.xs.isEmpty()) {
            return false;
        }

        double tMax = xs.max();
        double oMin = other.xs.min();
        if (tMax < oMin) {
            return false;
        }

        double tMin = xs.min();
        double oMax = other.xs.max();
        return !(tMin > oMax);
    }

    public LinearInterpolated apply(
        Operator           operator,
        LinearInterpolated other,
        double             from,
        double             to
    ) {
        LinearInterpolated result = new LinearInterpolated();
        if (!intersect(other)) {
            return result;
        }

        Set<Double> points = new TreeSet<Double>(CMP);
        points.add(from);
        points.add(to);

        this .pointsInRange(from, to, points);
        other.pointsInRange(from, to, points);


        for (double x: points) {
            if (!inGap(x) && !other.inGap(x)) {
                double y1 = this .value(x);
                double y2 = other.value(x);
                double y  = operator.evaluate(y1, y2);
                if (!Double.isNaN(y)) {
                    result.add(x, y);
                }
            }
        }

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