package de.intevation.gnv.math;

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

import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.MathException;

import org.apache.commons.math.analysis.UnivariateRealFunction;

import org.apache.commons.math.analysis.interpolation.SplineInterpolator;
import org.apache.commons.math.analysis.interpolation.UnivariateRealInterpolator;

import org.apache.log4j.Logger;

/**
 * @author <a href="mailto:ingo.weinzierl@intevation.de">Ingo Weinzierl</a>
 * @author <a href="mailto:sascha.teichmann@intevation.de">Sascha L. Teichmann</a>
 */
public class XYColumn
extends      Point2d
implements   UnivariateRealFunction
{
    private static Logger log = Logger.getLogger(XYColumn.class);

    protected List<HeightValue> values;

    protected transient UnivariateRealFunction curve;

    public XYColumn() {
        values = new ArrayList<HeightValue>();
    }

    public XYColumn(double x, double y, int i, int j) {
        super(x, y, i, j);
        values = new ArrayList<HeightValue>();
    }

    public void add(HeightValue value) {
        values.add(value);
    }

    public List<HeightValue> getValues() {
        return values;
    }

    public double value(double depth) {
        try {
            if (curve != null) {
                HeightValue h = values.get(0);
                // extrapolate beyond boundaries by repeating
                if (depth > h.z) return h.v;
                h = values.get(values.size()-1);
                if (depth < h.z) return h.v;
                return curve.value(depth);
            }
        }
        catch (FunctionEvaluationException fee) {
            log.error("evaluation failed", fee);
        }

        return Double.NaN;
    }

    public boolean prepare(XYDepth xyDepth) {
        if (curve == null) {
            int N = values.size();
            if (N == 0) {
                log.error("no points for interpolation");
                return false;
            }

            if (N == 1) {
                // only one value -> constant function
                curve = new ConstantFunction(values.get(0).v);
            }
            else { // more than on value
                double depth = xyDepth.depth(this);
                Collections.sort(values, HeightValue.INV_Z_COMPARATOR);

                // if there is no value at 0 repeat first value
                HeightValue first = values.get(0);
                if (first.z < 0d) {
                    values.add(0, new HeightValue(0d, first.v, first.k-1));
                    ++N;
                }

                // if there is no value at depth repeat last value
                HeightValue last = values.get(N-1);
                if (last.z > depth) {
                    values.add(new HeightValue(depth, last.v, last.k+1));
                    ++N;
                }
                if (N < 3) { // interpolate linear
                    first = values.get(0);
                    last  = values.get(N-1);
                    curve = new LinearFunction.Univariate(
                        first.z, first.v,
                        last.z,  last.v);
                }
                else { // higher degree interpolation
                    double [] z = new double[N];
                    double [] v = new double[N];
                    for (int i = 0; i < N; ++i) {
                        HeightValue h = values.get(N-1-i);
                        z[i] = h.z;
                        v[i] = h.v;
                    }
                    try {
                        curve = getInterpolator().interpolate(z, v);
                    }
                    catch (MathException me) {
                        log.error("interpolation failed", me);
                        return false;
                    }
                }
            }
        }
        return true;
    }

    protected UnivariateRealInterpolator getInterpolator() {
        return new SplineInterpolator();
    }
}
// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :