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 final 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 final class QPosition {

        protected int    index;
        protected double weight;

        public QPosition() {
        }

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

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

    } // class Position

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

        public Row() {
        }

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

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

        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 void interpolateW(
            Row           other,
            double        km,
            double []     iqs,
            double []     ows,
            WstValueTable table
        ) {
            double kmWeight = factor(km, this.km, other.km);

            QPosition qPosition = new QPosition();

            for (int i = 0; i < iqs.length; ++i) {
                if (table.getQPosition(km, iqs[i], qPosition) != null) {
                    double wt =       getW(qPosition);
                    double wo = other.getW(qPosition);
                    ows[i] = weight(kmWeight, wt, wo);
                }
                else {
                    ows[i] = Double.NaN;
                }
            }
        }

        public double [][] interpolateWQ(
            Row           other,
            double        km,
            int           steps,
            WstValueTable table
        ) {
            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,
                    table.getQIndex(i, km),
                    table.getQIndex(i, other.km));

                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, WstValueTable table) {

            int W = ws.length;

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

            double [] splineQ = new double[W];

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

            QPosition qPosition = new QPosition();

            for (int i = 0; i < W; ++i) {
                splineQ[i] = table.getQIndex(i, km);
                if (splineQ[i] < minQ) minQ = splineQ[i];
                if (splineQ[i] > maxQ) maxQ = splineQ[i];
            }

            double stepWidth = (maxQ - minQ)/steps;

            SplineInterpolator interpolator = new SplineInterpolator();

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

            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 getW(QPosition qPosition) {
            int    index  = qPosition.index;
            double weight = qPosition.weight;

            return weight == 1.0
                ? ws[index]
                : weight(weight, ws[index-1], ws[index]);
        }

        public double getW(
            Row       other,
            double    km,
            QPosition qPosition
        ) {
            double kmWeight = factor(km, this.km, other.km);

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

            double tw, ow;

            if (weight == 1.0) {
                tw = ws[index];
                ow = other.ws[index];
            }
            else {
                tw = weight(weight, ws[index-1], ws[index]); 
                ow = weight(weight, other.ws[index-1], other.ws[index]); 
            }

            return 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 WstValueTable(Column [] columns, List<Row> rows) {
        this.columns = columns;
        this.rows    = rows;
    }

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

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

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

        QPosition qPosition = new QPosition();

        if (rowIndex >= 0) { // direct row match
            Row row = rows.get(rowIndex);
            for (int i = 0; i < qs.length; ++i) {
                ws[i] = getQPosition(km, qs[i], qPosition) != null
                    ? row.getW(qPosition)
                    : Double.NaN;
            }
        }
        else { // needs bilinear interpolation
            rowIndex = -rowIndex -1;

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

        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, this);
        }

        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, this);
    }

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

        Row kmKey = new Row();

        QPosition nPosition = new QPosition();

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

            qs[i] = remap.eval(kms[i], getQ(qPosition, kms[i]));

            if (getQPosition(kms[i], qs[i], nPosition) == null) {
                ws[i] = Double.NaN;
                continue;
            }

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

            if (rowIndex >= 0) {
                // direct row match
                ws[i] = rows.get(rowIndex).getW(nPosition);
                continue;
            }

            rowIndex = -rowIndex -1;

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

    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 = getQPosition(kms[referenceIndex], q);

        if (qPosition == null) {
            // we cannot locate q at km
            return null;
        }

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

            rowIndex = Collections.binarySearch(rows, kmKey);

            qs[i] = getQ(qPosition, kms[i]);

            if (rowIndex >= 0) {
                // direct row match
                ws[i] = rows.get(rowIndex).getW(qPosition);
                continue;
            }

            rowIndex = -rowIndex -1;

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

            ws[i] = r1.getW(r2, kms[i], qPosition);
        }

        return qPosition;
    }

    public QPosition getQPosition(double km, double q) {
        return getQPosition(km, q, new QPosition());
    }

    public QPosition getQPosition(double km, double q, QPosition qPosition) {

        if (columns.length == 0) {
            return null;
        }

        double qLast = columns[0].getQRangeTree().findQ(km);

        if (Math.abs(qLast - q) < 0.00001) {
            return qPosition.set(0, 1d);
        }

        for (int i = 1; i < columns.length; ++i) {
            double qCurrent = columns[i].getQRangeTree().findQ(km);
            if (Math.abs(qCurrent - q) < 0.00001) {
                return qPosition.set(i, 1d);
            }

            double qMin, qMax;
            if (qLast < qCurrent) { qMin = qLast; qMax = qCurrent; }
            else                  { qMin = qCurrent; qMax = qLast; }

            if (q > qMin && q < qMax) {
                double weight = factor(q, qLast, qCurrent);
                return qPosition.set(i, weight);
            }
            qLast = qCurrent;
        }

        return null;
    }

    public double getQIndex(int index, double km) {
        return columns[index].getQRangeTree().findQ(km);
    }

    public double getQ(QPosition qPosition, double km) {
        int    index  = qPosition.index;
        double weight = qPosition.weight;

        if (weight == 1d) {
            return columns[index].getQRangeTree().findQ(km);
        }
        double q1 = columns[index-1].getQRangeTree().findQ(km);
        double q2 = columns[index  ].getQRangeTree().findQ(km);
        return weight(weight, q1, q2);
    }

    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;
        return a + factor*(b-a);
    }
}
// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :