package de.intevation.flys.artifacts.model;

import java.io.Serializable;

import de.intevation.flys.model.River;
import de.intevation.flys.model.Wst;
import de.intevation.flys.model.WstColumn;

import de.intevation.flys.backend.SessionHolder;

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

import org.hibernate.Session;
import org.hibernate.Query;
import org.hibernate.SQLQuery;

import org.hibernate.type.StandardBasicTypes;

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

    // TODO: put this into a property file
    public static final String SQL_POS_WQ = 
        "SELECT position, w, q, column_pos" +
        "    FROM wst_value_table"          +
        "    WHERE wst_id = :wst_id";

    public static final double DEFAULT_STEP_WIDTH = 0.01;

    public static class Column
    implements          Serializable
    {
        protected String name;

        public Column() {
        }

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

        public String getName() {
            return name;
        }

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

    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 double getWForKM(Row other, int index, double km) {
            double w1  = ws[index];
            double w2  = other.ws[index];
            double km1 = this.km;
            double km2 = other.km;
            return linear(km, km1, km2, w1, w2);
        }

        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.0) return -1;
            if (d > 0.0) return +1;
            return 0;
        }

        public double [][] cloneWQs() {
            return new double [][] {
                (double [])ws.clone(),
                (double [])qs.clone() };
        }

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

            int W1 =       ascendingWs();
            int W2 = other.ascendingWs();

            int W = Math.min(W1, W2);

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

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

            double minW = weight(factor, ws[0],   other.ws[0]);
            double maxW = weight(factor, ws[W-1], other.ws[W-1]);

            if (minW > maxW) {
                double t = minW;
                minW = maxW;
                maxW = t;
            }
            double [] x = new double[W];
            double [] y = new double[W];

            for (int i = 0; i < W; ++i) {
                x[i] = weight(factor, ws[i], other.ws[i]);
                y[i] = weight(factor, qs[i], other.qs[i]);
            }

            SplineInterpolator interpolator = new SplineInterpolator();
            PolynomialSplineFunction spline = interpolator.interpolate(x, y);

            double [] outWs =
                new double[(int)Math.ceil((maxW - minW)/stepWidth)];
            double [] outQs =
                new double[outWs.length];

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

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

        }

        public double [][] interpolateWQ(double stepWidth) {
            int W = ascendingWs(); // ignore back jumps

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

            double [] x = new double[W];
            double [] y = new double[W];

            for (int i = 0; i < W; ++i) {
                x[i] = ws[i];
                y[i] = qs[i];
            }

            SplineInterpolator interpolator = new SplineInterpolator();

            PolynomialSplineFunction spline = interpolator.interpolate(x, y);

            double minW = ws[0];
            double maxW = ws[W-1];

            double [] outWs =
                new double[(int)Math.ceil((maxW - minW)/stepWidth)];
            double [] outQs =
                new double[outWs.length];

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

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

        public int ascendingWs() {
            if (ws.length < 2) {
                return ws.length;
            }

            int idx = 1;

            for (; idx < ws.length; ++idx) {
                if (Double.isNaN(ws[idx]) || ws[idx] < ws[idx-1]) {
                    return idx;
                }
            }

            return idx;
        }

        public double [][] weightWQs(Row other, double km) {
            int W = Math.min(ws.length, other.ws.length);
            double factor = factor(km, this.km, other.km);

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

            for (int i = 0; i < W; ++i) {
                outWs[i] = weight(factor, ws[i], other.ws[i]);
                outQs[i] = weight(factor, qs[i], other.qs[i]);
            }

            return new double [][] { outWs, outQs };
        }
    }
    // class Row

    protected List<Row> rows;

    protected Column [] columns;

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

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


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

        double [] ws = new double[qs.length];

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

    public double [][] interpolateWQ(
        double  km,
        double  stepWidth, 
        boolean checkAscending
    ) {
        int rowIndex = Collections.binarySearch(rows, new Row(km));

        if (rowIndex >= 0) { // direct row match
            Row row = rows.get(rowIndex);
            return checkAscending
                ? row.interpolateWQ(stepWidth)
                : row.cloneWQs();
        }

        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 checkAscending
            ? r1.interpolateWQ(r2, km, stepWidth)
            : r1.weightWQs(r2, km);
    }

    public static WstValueTable getTable(River river) {
        return getTable(river, 0);
    }

    public static WstValueTable getTable(River river, int kind) {

        Session session = SessionHolder.HOLDER.get();

        Query query = session.createQuery(
            "from Wst where river=:river and kind=:kind");
        query.setParameter("river", river);
        query.setInteger("kind",    kind);

        List<Wst> wsts = query.list();

        if (wsts.isEmpty()) {
            return null;
        }

        Wst wst = wsts.get(0);

        // TODO: Do this sorting at database level
        List<WstColumn> wstColumns = new ArrayList(wst.getColumns());
        Collections.sort(wstColumns, new Comparator<WstColumn>() {
            public int compare(WstColumn a, WstColumn b) {
                int pa = a.getPosition();
                int pb = b.getPosition();
                if (pa < pb) return -1;
                if (pa > pb) return +1;
                return 0;
            }
        });

        Column [] columns = new Column[wstColumns.size()];
        for (int i = 0; i < columns.length; ++i) {
            columns[i] = new Column(wstColumns.get(i).getName());
        }

        // using native SQL here to avoid myriad of small objects.
        SQLQuery sqlQuery = session.createSQLQuery(SQL_POS_WQ)
            .addScalar("position",   StandardBasicTypes.DOUBLE)
            .addScalar("w",          StandardBasicTypes.DOUBLE)
            .addScalar("q",          StandardBasicTypes.DOUBLE)
            .addScalar("column_pos", StandardBasicTypes.INTEGER);

        sqlQuery.setInteger("wst_id", wst.getId());

        WstValueTable valueTable = new WstValueTable(columns);

        int lastColumnNo = -1;
        Row row          = null;

        Double lastQ     = -Double.MAX_VALUE;
        boolean qSorted  = true;

        for (Object  r: sqlQuery.list()) {
            Object[] result = (Object[]) r;

            double km    = (Double) result[0];
            Double w     = (Double) result[1];
            Double q     = (Double) result[2];
            int columnNo = (Integer)result[3];

            if (columnNo > lastColumnNo) { // new row
                if (row != null) {
                    row.qSorted = qSorted;
                    valueTable.rows.add(row);
                }
                row = new Row(
                    km,
                    new double[columnNo+1],
                    new double[columnNo+1]);
                lastQ = -Double.MAX_VALUE;
                qSorted = true;
            }

            row.ws[columnNo] = w != null ? w : Double.NaN;
            row.qs[columnNo] = q != null ? q : Double.NaN;

            if (qSorted && (q == null || lastQ > q)) {
                qSorted = false;
            }
            lastQ = q;

            lastColumnNo = columnNo;
        }

        if (row != null) {
            valueTable.rows.add(row);
        }

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