package de.intevation.flys.artifacts.model;

import java.io.Serializable;


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

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

import org.apache.log4j.Logger;

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

import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;

import org.apache.commons.math.ArgumentOutsideDomainException;

public class WstValueTable
implements   Serializable
{
    private static Logger log = Logger.getLogger(WstValueTable.class);

    public static final int DEFAULT_Q_STEPS = 500;

    public static class Column
    implements          Serializable
    {
        protected String name;

        protected QRangeTree qRangeTree;

        public Column() {
        }

        public Column(String name) {
            this.name = name;
        }

        public String getName() {
            return name;
        }

        public void setName(String name) {
            this.name = name;
        }

        public QRangeTree getQRangeTree() {
            return qRangeTree;
        }

        public void setQRangeTree(QRangeTree qRangeTree) {
            this.qRangeTree = qRangeTree;
        }
    }
    // class Column

    public static class QPosition {

        protected double q;
        protected int    index;
        protected double weight;

        public QPosition() {
        }

        public QPosition(double q, int index, double weight) {
            this.index  = index;
            this.weight = weight;
        }

    } // class Position

    public static class Row
    implements          Serializable, Comparable<Row>
    {
        double    km;
        double [] ws;
        double [] qs;
        boolean   qSorted;

        public Row() {
        }

        public Row(double km) {
            this.km = km;
        }

        public Row(double km, double [] ws, double [] qs) {
            this(km);
            this.ws = ws;
            this.qs = qs;
        }

        public static final double linear(
            double x,
            double x1, double x2,
            double y1, double y2
        ) {
            // y1 = m*x1 + b
            // y2 = m*x2 + b
            // y2 - y1 = m*(x2 - x1)
            // m = (y2 - y1)/(x2 - x1) # x2 != x1
            // b = y1 - m*x1

            if (x1 == x2) {
                return 0.5*(y1 + y2);
            }
            double m = (y2 - y1)/(x2 - x1);
            double b = y1 - m*x1;
            return x*m + b;
        }

        public static final double factor(double x, double p1, double p2) {
            // 0 = m*p1 + b <=> b = -m*p1
            // 1 = m*p2 + b
            // 1 = m*(p2 - p1)
            // m = 1/(p2 - p1) # p1 != p2
            // f(x) = x/(p2-p1) - p1/(p2-p1) <=> (x-p1)/(p2-p1)

            return p1 == p2 ? 0.0 : (x-p1)/(p2-p1);
        }

        public static final double weight(double factor, double a, double b) {
            return (1.0-factor)*a + factor*b;
        }

        public void interpolateW(
            Row       other,
            double    km,
            double [] input,
            double [] output
        ) {
            double factor = factor(km, this.km, other.km);

            for (int i = 0; i < input.length; ++i) {
                double q = input[i];
                int idx1 =       getQIndex(q);
                int idx2 = other.getQIndex(q);

                double w1 = idx1 >= 0
                    ? ws[idx1]
                    : interpolateW(-idx1-1, q);

                double w2 = idx2 >= 0
                    ? other.ws[idx2]
                    : other.interpolateW(-idx2-1, q);

                output[i] = weight(factor, w1, w2);
            }
        }

        public double interpolateW(int idx, double q) {
            return idx < 1 || idx >= qs.length
                ? Double.NaN // do not extrapolate
                : linear(q, qs[idx-1], qs[idx], ws[idx-1], ws[idx]);
        }

        public double interpolateW(double q) {
            if (Double.isNaN(q)) return Double.NaN;
            int index = getQIndex(q);
            return index >= 0 ? ws[index] : interpolateW(-index -1, q);
        }

        public int getQIndex(double q) {
            return qSorted ? binaryQIndex(q) : linearQIndex(q);
        }

        public int binaryQIndex(double q) {
            return Arrays.binarySearch(qs, q);
        }

        public int linearQIndex(double q) {
            switch (qs.length) {
                case 0: break;
                case 1:
                    if (qs[0] == q) return 0;
                    if (qs[0] <  q) return -(1+1);
                    return -(0+1);
                default:
                    for (int i = 1; i < qs.length; ++i) {
                        double qa = qs[i-1];
                        double qb = qs[i];
                        if (qa == q) return i-1;
                        if (qb == q) return i;
                        if (qa > qb) { double t = qa; qa = qb; qb = t; }
                        if (q > qa && q < qb) return -(i+1);
                    }
                    return -qs.length;
            }

            return -1;
        }

        public int compareTo(Row other) {
            double d = km - other.km;
            if (d < -0.0001) return -1;
            if (d >  0.0001) return +1;
            return 0;
        }

        public double [][] interpolateWQ(Row other, double km, int steps) {

            int W = Math.min(ws.length, other.ws.length);

            if (W < 1) {
                return new double[2][0];
            }

            double factor = factor(km, this.km, other.km);

            double [] splineQ = new double[W];
            double [] splineW = new double[W];

            double minQ =  Double.MAX_VALUE;
            double maxQ = -Double.MAX_VALUE;

            for (int i = 0; i < W; ++i) {
                double wws = weight(factor, ws[i], other.ws[i]);
                double wqs = weight(factor, qs[i], other.qs[i]);
                splineW[i] = wws;
                splineQ[i] = wqs;
                if (wqs < minQ) minQ = wqs;
                if (wqs > maxQ) maxQ = wqs;
            }

            double stepWidth = (maxQ - minQ)/steps;

            SplineInterpolator interpolator = new SplineInterpolator();
            PolynomialSplineFunction spline =
                interpolator.interpolate(splineQ, splineW);

            double [] outWs = new double[steps];
            double [] outQs = new double[steps];

            try {
                double q = minQ;
                for (int i = 0; i < outWs.length; ++i, q += stepWidth) {
                    outWs[i] = spline.value(outQs[i] = q);
                }
            }
            catch (ArgumentOutsideDomainException aode) {
                log.error("Spline interpolation failed.", aode);
            }

            return new double [][] { outWs, outQs };
        }

        public double [][] interpolateWQ(int steps) {

            int W = ws.length;

            if (W < 1) {
                return new double[2][0];
            }

            double [] splineW = new double[W];
            double [] splineQ = new double[W];

            double minQ =  Double.MAX_VALUE; 
            double maxQ = -Double.MAX_VALUE; 

            for (int i = 0; i < W; ++i) {
                splineW[i] = ws[i];
                splineQ[i] = qs[i];
                if (qs[i] < minQ) minQ = qs[i];
                if (qs[i] > maxQ) maxQ = qs[i];
            }

            double stepWidth = (maxQ - minQ)/steps;

            SplineInterpolator interpolator = new SplineInterpolator();

            PolynomialSplineFunction spline =
                interpolator.interpolate(splineQ, splineW);

            double [] outWs = new double[steps];
            double [] outQs = new double[steps];

            try {
                double q = minQ;
                for (int i = 0; i < outWs.length; ++i, q += stepWidth) {
                    outWs[i] = spline.value(outQs[i] = q);
                }
            }
            catch (ArgumentOutsideDomainException aode) {
                log.error("Spline interpolation failed.", aode);
            }

            return new double [][] { outWs, outQs };
        }

        public QPosition getQPosition(double q) {
            int qIndex = getQIndex(q);

            if (qIndex >= 0) {
                // direct column match
                return new QPosition(q, qIndex, 1.0);
            }

            qIndex = -qIndex -1;

            if (qIndex < 0 || qIndex >= qs.length-1) {
                // do not extrapolate
                return null;
            }

            return new QPosition(
                q, qIndex, factor(q, qs[qIndex-1], qs[qIndex]));
        }

        public QPosition getQPosition(Row other, double km, double q) {

            QPosition qpt = getQPosition(q);
            QPosition qpo = other.getQPosition(q);

            if (qpt == null || qpo == null) {
                return null;
            }

            double kmWeight = factor(km, this.km, other.km);

            // XXX: Index interpolation is a bit sticky here!

            int index = (int)Math.round(
                weight(kmWeight, qpt.index, qpo.index));

            double weight = weight(kmWeight, qpt.weight, qpo.weight);

            return new QPosition(q, index, weight);
        }

        public void storeWQ(
            QPosition qPosition,
            double [] ows,
            double [] oqs,
            int       index
        ) {
            int    qIdx    = qPosition.index;
            double qWeight = qPosition.weight;

            if (qWeight == 1.0) {
                oqs[index] = qs[qIdx];
                ows[index] = ws[qIdx];
            }
            else {
                oqs[index] = weight(qWeight, qs[qIdx-1], qs[qIdx]); 
                ows[index] = weight(qWeight, ws[qIdx-1], ws[qIdx]); 
            }
        }

        public void storeWQ(
            Row       other,
            double    km,
            QPosition qPosition,
            double [] ows,
            double [] oqs,
            int       index
        ) {
            double kmWeight = factor(km, this.km, other.km);

            double tq, tw;
            double oq, ow;

            int    qIdx    = qPosition.index;
            double qWeight = qPosition.weight;

            if (qWeight == 1.0) {
                tq = qs[qIdx];
                tw = ws[qIdx];
                oq = other.qs[qIdx];
                ow = other.ws[qIdx];
            }
            else {
                tq = weight(qWeight, qs[qIdx-1], qs[qIdx]); 
                tw = weight(qWeight, ws[qIdx-1], ws[qIdx]); 
                oq = weight(qWeight, other.qs[qIdx-1], other.qs[qIdx]); 
                ow = weight(qWeight, other.ws[qIdx-1], other.ws[qIdx]); 
            }

            oqs[index] = weight(kmWeight, tq, oq);
            ows[index] = weight(kmWeight, tw, ow);
        }

        public void storeWQ(
            QPosition   qPosition,
            LinearRemap remap,
            double []   ows,
            double []   oqs,
            int         index
        ) {
            int    qIdx    = qPosition.index;
            double qWeight = qPosition.weight;

            oqs[index] = remap.eval(km, qWeight == 1.0
                ? qs[qIdx]
                : weight(qWeight, qs[qIdx-1], qs[qIdx])); 

            ows[index] = interpolateW(oqs[index]);
        }

        public void storeWQ(
            Row         other,
            double      km,
            QPosition   qPosition,
            LinearRemap remap,
            double []   ows,
            double []   oqs,
            int         index
        ) {
            double kmWeight = factor(km, this.km, other.km);

            double tq, tw;
            double oq, ow;

            int    qIdx    = qPosition.index;
            double qWeight = qPosition.weight;

            if (qWeight == 1.0) {
                tq = remap.eval(this.km,  qs[qIdx]);
                oq = remap.eval(other.km, other.qs[qIdx]);
            }
            else {
                tq = remap.eval(
                    this.km,
                    weight(qWeight, qs[qIdx-1], qs[qIdx]));
                oq = remap.eval(
                    other.km,
                    weight(qWeight, other.qs[qIdx-1], other.qs[qIdx])); 
            }

            tw = interpolateW(tq);
            ow = other.interpolateW(oq);

            oqs[index] = weight(kmWeight, tq, oq);
            ows[index] = weight(kmWeight, tw, ow);
        }
    }
    // class Row

    protected List<Row> rows;

    protected Column [] columns;

    public WstValueTable() {
        rows = new ArrayList<Row>();
    }

    public WstValueTable(Column [] columns) {
        this();
        this.columns = columns;
    }

    public void sortRows() {
        Collections.sort(rows);
    }


    public double [] interpolateW(double km, double [] qs) {
        return interpolateW(km, qs, new double[qs.length]);
    }

    public double [] interpolateW(double km, double [] qs, double [] ws) {

        int rowIndex = Collections.binarySearch(rows, new Row(km));

        if (rowIndex >= 0) { // direct row match
            Row row = rows.get(rowIndex);
            for (int i = 0; i < qs.length; ++i) {
                ws[i] = row.interpolateW(qs[i]);
            }
        }
        else { // needs bilinear interpolation
            rowIndex = -rowIndex -1;

            if (rowIndex < 1 || rowIndex >= rows.size()) {
                // do not extrapolate
                Arrays.fill(ws, Double.NaN);
            }
            else {
                rows.get(rowIndex-1).interpolateW(
                    rows.get(rowIndex),
                    km, qs, ws);
            }
        }

        return ws;
    }


    public double [][] interpolateWQ(double km) {
        return interpolateWQ(km, DEFAULT_Q_STEPS);
    }

    public double [][] interpolateWQ(double km, int steps) {

        int rowIndex = Collections.binarySearch(rows, new Row(km));

        if (rowIndex >= 0) { // direct row match
            Row row = rows.get(rowIndex);
            return row.interpolateWQ(steps);
        }

        rowIndex = -rowIndex -1;

        if (rowIndex < 1 || rowIndex >= rows.size()) {
            // do not extrapolate
            return new double[2][0];
        }

        Row r1 = rows.get(rowIndex-1);
        Row r2 = rows.get(rowIndex);

        return r1.interpolateWQ(r2, km, steps);
    }

    public void interpolate(
        double []   kms,
        double []   ws,
        double []   qs,
        QPosition   qPosition,
        LinearRemap remap
    ) {
        int R1 = rows.size()-1;

        Row kmKey = new Row();
        for (int i = 0; i < kms.length; ++i) {
            kmKey.km = kms[i];
            int rowIndex = Collections.binarySearch(rows, kmKey);

            if (rowIndex >= 0) {
                // direct row match
                rows.get(rowIndex).storeWQ(qPosition, remap, ws, qs, i);
                continue;
            }

            rowIndex = -rowIndex -1;

            if (rowIndex < 1 || rowIndex >= R1) {
                // do not extrapolate
                ws[i] = qs[i] = Double.NaN;
                continue;
            }
            Row r1 = rows.get(rowIndex-1);
            Row r2 = rows.get(rowIndex);
            r1.storeWQ(r2, kms[i], qPosition, remap, ws, qs, i);
        }
    }

    public QPosition interpolate(
        double    q,
        int       referenceIndex,
        double [] kms,
        double [] ws,
        double [] qs
    ) {
        Row kmKey = new Row(kms[referenceIndex]);

        int rowIndex = Collections.binarySearch(rows, kmKey);

        int R1 = rows.size()-1;

        QPosition qPosition = null;

        if (rowIndex >= 0) { // direct row match
            qPosition = rows.get(rowIndex).getQPosition(q);
        }
        else {
            rowIndex = -rowIndex -1;

            if (rowIndex < 1 || rowIndex >= R1) {
                // do not extrapolate
                return null;
            }

            // interpolation case
            Row r1 = rows.get(rowIndex-1);
            Row r2 = rows.get(rowIndex);

            qPosition = r1.getQPosition(r2, kms[referenceIndex], q);
        }

        if (qPosition == null) {
            // we cannot locate q in reference row
            return null;
        }

        for (int i = 0; i < kms.length; ++i) {
            kmKey.km = kms[i];

            rowIndex = Collections.binarySearch(rows, kmKey);

            if (rowIndex >= 0) {
                // direct row match
                rows.get(rowIndex).storeWQ(qPosition, ws, qs, i);
                continue;
            }

            rowIndex = -rowIndex -1;

            if (rowIndex < 1 || rowIndex >= R1) {
                // do not extrapolate
                ws[i] = qs[i] = Double.NaN;
                continue;
            }
            Row r1 = rows.get(rowIndex-1);
            Row r2 = rows.get(rowIndex);
            r1.storeWQ(r2, kms[i], qPosition, ws, qs, i);
        }

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