/* 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 static org.dive4elements.river.backend.utils.EpsilonComparator.CMP;

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

import org.apache.commons.math.ArgumentOutsideDomainException;
import org.apache.commons.math.analysis.interpolation.SplineInterpolator;
import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
import org.apache.commons.math.exception.MathIllegalArgumentException;
import org.apache.log4j.Logger;
import org.dive4elements.river.artifacts.math.Function;
import org.dive4elements.river.artifacts.math.Linear;
import org.dive4elements.river.utils.DoubleUtil;

import gnu.trove.TDoubleArrayList;

/**
 * W, Q and km data from database 'wsts' spiced with interpolation algorithms.
 */
public class WstValueTable
implements   Serializable
{
    private static Logger log = Logger.getLogger(WstValueTable.class);

    public static final int DEFAULT_Q_STEPS = 500;

    public static final int RELATE_WS_SAMPLES = 200;

    public static final double W_EPSILON = 0.000001;

    /**
     * A Column in the table, typically representing one measurement session.
     */
    public static final class Column
    implements                Serializable
    {
        protected String name;

        protected QRangeTree qRangeTree;

        public Column() {
        }

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

        public String getName() {
            return this.name;
        }

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

        public QRangeTree getQRangeTree() {
            return this.qRangeTree;
        }

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

    /**
     * A (weighted) position used for interpolation.
     */
    public static final class QPosition {

        protected int    index;
        protected double weight;

        public QPosition() {
        }

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

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

    } // class Position

    public static final class SplineFunction {

        public PolynomialSplineFunction spline;
        public double []                splineQs;
        public double []                splineWs;

        public SplineFunction(
                final PolynomialSplineFunction spline,
                final double []                splineQs,
                final double []                splineWs
                ) {
            this.spline   = spline;
            this.splineQs = splineQs;
            this.splineWs = splineWs;
        }

        public double [][] sample(
                final int         numSamples,
                final double      km,
                final Calculation errors
                ) {
            final double minQ = getQMin();
            final double maxQ = getQMax();

            final double [] outWs = new double[numSamples];
            final double [] outQs = new double[numSamples];

            Arrays.fill(outWs, Double.NaN);
            Arrays.fill(outQs, Double.NaN);

            final double stepWidth = (maxQ - minQ)/numSamples;

            try {
                double q = minQ;
                for (int i = 0; i < outWs.length; ++i, q += stepWidth) {
                    outWs[i] = this.spline.value(outQs[i] = q);
                }
            }
            catch (final ArgumentOutsideDomainException aode) {
                if (errors != null) {
                    errors.addProblem(km, "spline.interpolation.failed");
                }
                log.error("spline interpolation failed.", aode);
            }

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

        public double getQMin() {
            return Math.min(this.splineQs[0], this.splineQs[this.splineQs.length-1]);
        }

        public double getQMax() {
            return Math.max(this.splineQs[0], this.splineQs[this.splineQs.length-1]);
        }

        /** Constructs a continues index between the columns to Qs. */
        public PolynomialSplineFunction createIndexQRelation() {

            final double [] indices = new double[this.splineQs.length];
            for (int i = 0; i < indices.length; ++i) {
                indices[i] = i;
            }

            try {
                final SplineInterpolator interpolator = new SplineInterpolator();
                return interpolator.interpolate(indices, this.splineQs);
            }
            catch (final MathIllegalArgumentException miae) {
                // Ignore me!
            }
            return null;
        }
    } // class SplineFunction

    /**
     * A row, typically a position where measurements were taken.
     */
    public static final class Row
    implements                Serializable, Comparable<Row>
    {
        double    km;
        double [] ws;

        public Row() {
        }

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

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

        /**
         * Sort Qs and Ws for this Row over Q.
         */
        private double[][] getSortedWQ(
                final WstValueTable table,
                final Calculation   errors
                ) {
            final int W = this.ws.length;

            if (W < 1) {
                if (errors != null) {
                    errors.addProblem(this.km, "no.ws.found");
                }
                return new double[][] {{Double.NaN}, {Double.NaN}};
            }

            final double [] sortedWs = this.ws;
            final double [] sortedQs = new double[W];

            for (int i=0; i < W; ++i) {
                final double q = table.getQIndex(i, this.km);
                if (Double.isNaN(q) && errors != null) {
                    errors.addProblem(
                            this.km, "no.q.found.in.column", i+1);
                }
                sortedQs[i] = q;
            }

            DoubleUtil.sortByFirst(sortedQs, sortedWs);

            return new double[][] { sortedWs, sortedQs };
        }

        /**
         * Return an array of Qs and Ws for the given km between
         * this Row and another, sorted over Q.
         */
        private double[][] getSortedWQ(
                final Row           other,
                final double        km,
                final WstValueTable table,
                final Calculation   errors
                ) {
            final int W = Math.min(this.ws.length, other.ws.length);

            if (W < 1) {
                if (errors != null) {
                    errors.addProblem("no.ws.found");
                }
                return new double[][] {{Double.NaN}, {Double.NaN}};
            }

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

            final double [] sortedQs = new double[W];
            final double [] sortedWs = new double[W];

            for (int i = 0; i < W; ++i) {
                final double wws = Linear.weight(factor, this.ws[i], other.ws[i]);
                final double wqs = table.getQIndex(i, km);

                if (Double.isNaN(wws) || Double.isNaN(wqs)) {
                    if (errors != null) {
                        errors.addProblem(km, "cannot.find.w.or.q");
                    }
                }

                sortedWs[i] = wws;
                sortedQs[i] = wqs;
            }

            DoubleUtil.sortByFirst(sortedQs, sortedWs);

            return new double[][] { sortedWs, sortedQs };
        }

        /**
         * Find Qs matching w in an array of Qs and Ws sorted over Q.
         */
        private double[] findQsForW(final double w, final double[][] sortedWQ) {
            final int W = sortedWQ[0].length;

            final double[] sortedQs = sortedWQ[1];
            final double[] sortedWs = sortedWQ[0];

            final TDoubleArrayList qs = new TDoubleArrayList();

            if (W > 0 && Math.abs(sortedWs[0]-w) < W_EPSILON) {
                final double q = sortedQs[0];
                if (!Double.isNaN(q)) {
                    qs.add(q);
                }
            }

            for (int i = 1; i < W; ++i) {
                final double w2 = sortedWs[i];
                if (Double.isNaN(w2)) {
                    continue;
                }
                if (Math.abs(w2-w) < W_EPSILON) {
                    final double q = sortedQs[i];
                    if (!Double.isNaN(q)) {
                        qs.add(q);
                    }
                    continue;
                }
                final double w1 = sortedWs[i-1];
                if (Double.isNaN(w1)) {
                    continue;
                }

                if (w < Math.min(w1, w2) || w > Math.max(w1, w2)) {
                    continue;
                }

                final double q1 = sortedQs[i-1];
                final double q2 = sortedQs[i];
                if (Double.isNaN(q1) || Double.isNaN(q2)) {
                    continue;
                }

                final double q = Linear.linear(w, w1, w2, q1, q2);
                qs.add(q);
            }

            return qs.toNativeArray();
        }

        /**
         * Compare according to place of measurement (km).
         */
        @Override
        public int compareTo(final Row other) {
            return CMP.compare(this.km, other.km);
        }

        /**
         * Interpolate Ws, given Qs and a km.
         *
         * @param iqs Given ("input") Qs.
         * @param ows Resulting ("output") Ws.
         * @param table Table of which to use data for interpolation.
         */
        public void interpolateW(
                final Row           other,
                final double        km,
                final double []     iqs,
                final double []     ows,
                final WstValueTable table,
                final Calculation   errors
                ) {
            final double kmWeight = Linear.factor(km, this.km, other.km);

            final QPosition qPosition = new QPosition();

            for (int i = 0; i < iqs.length; ++i) {
                if (table.getQPosition(km, iqs[i], qPosition) != null) {
                    final double wt =       getW(qPosition);
                    final double wo = other.getW(qPosition);
                    if (Double.isNaN(wt) || Double.isNaN(wo)) {
                        if (errors != null) {
                            errors.addProblem(
                                    km, "cannot.find.w.for.q", iqs[i]);
                        }
                        ows[i] = Double.NaN;
                    }
                    else {
                        ows[i] = Linear.weight(kmWeight, wt, wo);
                    }
                }
                else {
                    if (errors != null) {
                        errors.addProblem(km, "cannot.find.q", iqs[i]);
                    }
                    ows[i] = Double.NaN;
                }
            }
        }


        public SplineFunction createSpline(
                final WstValueTable table,
                final Calculation   errors
                ) {
            final double[][] sortedWQ = getSortedWQ(table, errors);

            try {
                final SplineInterpolator interpolator = new SplineInterpolator();
                final PolynomialSplineFunction spline =
                        interpolator.interpolate(sortedWQ[1], sortedWQ[0]);

                return new SplineFunction(spline, sortedWQ[1], this.ws);
            }
            catch (final MathIllegalArgumentException miae) {
                if (errors != null) {
                    errors.addProblem(this.km, "spline.creation.failed");
                }
                log.error("spline creation failed", miae);
            }
            return null;
        }

        public SplineFunction createSpline(
                final Row           other,
                final double        km,
                final WstValueTable table,
                final Calculation   errors
                ) {
            final double[][] sortedWQ = getSortedWQ(other, km, table, errors);

            final SplineInterpolator interpolator = new SplineInterpolator();

            try {
                final PolynomialSplineFunction spline =
                        interpolator.interpolate(sortedWQ[1], sortedWQ[0]);

                return new SplineFunction(spline, sortedWQ[1], sortedWQ[0]);
            }
            catch (final MathIllegalArgumentException miae) {
                if (errors != null) {
                    errors.addProblem(km, "spline.creation.failed");
                }
                log.error("spline creation failed", miae);
            }

            return null;
        }

        public double [][] interpolateWQ(
                final Row           other,
                final double        km,
                final int           steps,
                final WstValueTable table,
                final Calculation   errors
                ) {
            final SplineFunction sf = createSpline(other, km, table, errors);

            return sf != null
                    ? sf.sample(steps, km, errors)
                            : new double[2][0];
        }


        public double [][] interpolateWQ(
                final int           steps,
                final WstValueTable table,
                final Calculation   errors
                ) {
            final SplineFunction sf = createSpline(table, errors);

            return sf != null
                    ? sf.sample(steps, this.km, errors)
                            : new double[2][0];
        }


        public double getW(final QPosition qPosition) {
            final int    index  = qPosition.index;
            final double weight = qPosition.weight;

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

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

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

            double tw, ow;

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

            return Linear.weight(kmWeight, tw, ow);
        }

        public double [] findQsForW(
                final double        w,
                final WstValueTable table,
                final Calculation   errors
                ) {
            log.debug("Find Qs for given W at tabulated km " + this.km);
            return findQsForW(w, getSortedWQ(table, errors));
        }

        public double [] findQsForW(
                final Row           other,
                final double        w,
                final double        km,
                final WstValueTable table,
                final Calculation   errors
                ) {
            log.debug("Find Qs for given W at non-tabulated km " + km);
            return findQsForW(w, getSortedWQ(other, km, table, errors));
        }

        public double [] getMinMaxW(final double [] result) {
            double minW =  Double.MAX_VALUE;
            double maxW = -Double.MAX_VALUE;
            for (int i = 0; i < this.ws.length; ++i) {
                final double w = this.ws[i];
                if (w < minW) minW = w;
                if (w > maxW) maxW = w;
            }
            result[0] = minW;
            result[1] = maxW;
            return result;
        }

        public double [] getMinMaxW(final Row other, final double km, final double [] result) {
            final double [] m1 = this .getMinMaxW(new double [2]);
            final double [] m2 = other.getMinMaxW(new double [2]);
            final double factor = Linear.factor(km, this.km, other.km);
            result[0] = Linear.weight(factor, m1[0], m2[0]);
            result[1] = Linear.weight(factor, m1[1], m2[1]);
            return result;
        }
    } // class Row

    /** Rows in table. */
    protected List<Row> rows;

    /** Columns in table. */
    protected Column [] columns;

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

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

    /**
     * @param columns The WST-columns.
     * @param rows A list of Rows that must be sorted by km.
     */
    public WstValueTable(final Column [] columns, final List<Row> rows) {
        this.columns = columns;
        this.rows    = rows;
    }

    public Column [] getColumns() {
        return this.columns;
    }

    /**
     * @param km Given kilometer.
     * @param qs Given Q values.
     * @param ws output parameter.
     */
    public double [] interpolateW(final double km, final double [] qs, final double [] ws) {
        return interpolateW(km, qs, ws, null);
    }


    /**
     * @param ws (output parameter), gets returned.
     * @return output parameter ws.
     */
    public double [] interpolateW(
            final double      km,
            final double []   qs,
            final double []   ws,
            final Calculation errors
            ) {
        int rowIndex = Collections.binarySearch(this.rows, new Row(km));

        final QPosition qPosition = new QPosition();

        if (rowIndex >= 0) { // direct row match
            final Row row = this.rows.get(rowIndex);
            for (int i = 0; i < qs.length; ++i) {
                if (getQPosition(km, qs[i], qPosition) == null) {
                    if (errors != null) {
                        errors.addProblem(km, "cannot.find.q", qs[i]);
                    }
                    ws[i] = Double.NaN;
                }
                else {
                    if (Double.isNaN(ws[i] = row.getW(qPosition))
                            && errors != null) {
                        errors.addProblem(
                                km, "cannot.find.w.for.q", qs[i]);
                    }
                }
            }
        }
        else { // needs bilinear interpolation
            rowIndex = -rowIndex -1;

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

        return ws;
    }

    /**
     * Interpolates a W for a km using a Q column position
     */
    public double interpolateW(final double km, final QPosition qPosition, final Calculation errors) {

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

        if (rowIndex >= 0) { // direct row match
            final Row row = this.rows.get(rowIndex);
            return row.getW(qPosition);
        }
        else { // needs bilinear interpolation
            rowIndex = -rowIndex - 1;
            if ((rowIndex <= 0) || (rowIndex >= this.rows.size()))
                return Double.NaN;
            else {
                final Row r1 = this.rows.get(rowIndex - 1);
                final Row r2 = this.rows.get(rowIndex);
                final double w1 = r1.getW(qPosition);
                final double w2 = r2.getW(qPosition);
                if (Double.isNaN(w1) || Double.isNaN(w2))
                    return Double.NaN;
                else {
                    final double kmWeight = Linear.factor(km, r1.km, r2.km);
                    return Linear.weight(kmWeight, w1, w2);
                }
            }
        }
    }

    public double [] getMinMaxQ(final double km) {
        return getMinMaxQ(km, new double [2]);
    }

    public double [] getMinMaxQ(final double km, final double [] result) {
        double minQ =  Double.MAX_VALUE;
        double maxQ = -Double.MAX_VALUE;

        for (int i = 0; i < this.columns.length; ++i) {
            final double q = this.columns[i].getQRangeTree().findQ(km);
            if (!Double.isNaN(q)) {
                if (q < minQ) minQ = q;
                if (q > maxQ) maxQ = q;
            }
        }

        if (minQ < Double.MAX_VALUE) {
            result[0] = minQ;
            result[1] = maxQ;
            return result;
        }

        return null;
    }

    public double [] getMinMaxQ(double from, double to, double step) {
        final double [] result = new double[2];

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

        if (from > to) {
            final double tmp = from;
            from = to;
            to = tmp;
        }

        step = Math.max(Math.abs(step), 0.0001);

        double d = from;
        for (; d <= to; d += step) {
            if (getMinMaxQ(d, result) != null) {
                if (result[0] < minQ) minQ = result[0];
                if (result[1] > maxQ) maxQ = result[1];
            }
        }

        if (d != to) {
            if (getMinMaxQ(to, result) != null) {
                if (result[0] < minQ) minQ = result[0];
                if (result[1] > maxQ) maxQ = result[1];
            }
        }

        return minQ < Double.MAX_VALUE
                ? new double [] { minQ, maxQ }
        : null;
    }

    public double [] getMinMaxW(final double km) {
        return getMinMaxW(km, new double [2]);

    }
    public double [] getMinMaxW(final double km, final double [] result) {
        int rowIndex = Collections.binarySearch(this.rows, new Row(km));

        if (rowIndex >= 0) {
            return this.rows.get(rowIndex).getMinMaxW(result);
        }

        rowIndex = -rowIndex -1;

        if (rowIndex < 1 || rowIndex >= this.rows.size()) {
            // do not extrapolate
            return null;
        }

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

        return r1.getMinMaxW(r2, km, result);
    }

    public double [] getMinMaxW(double from, double to, double step) {
        final double [] result = new double[2];

        double minW =  Double.MAX_VALUE;
        double maxW = -Double.MAX_VALUE;

        if (from > to) {
            final double tmp = from;
            from = to;
            to = tmp;
        }

        step = Math.max(Math.abs(step), 0.0001);

        double d = from;
        for (; d <= to; d += step) {
            if (getMinMaxW(d, result) != null) {
                if (result[0] < minW) minW = result[0];
                if (result[1] > maxW) maxW = result[1];
            }
        }

        if (d != to) {
            if (getMinMaxW(to, result) != null) {
                if (result[0] < minW) minW = result[0];
                if (result[1] > maxW) maxW = result[1];
            }
        }

        return minW < Double.MAX_VALUE
                ? new double [] { minW, maxW }
        : null;
    }

    /**
     * Interpolate W and Q values at a given km.
     */
    public double [][] interpolateWQ(final double km) {
        return interpolateWQ(km, null);
    }

    /**
     * Interpolate W and Q values at a given km.
     *
     * @param errors where to store errors.
     *
     * @return double double array, first index Ws, second Qs.
     */
    public double [][] interpolateWQ(final double km, final Calculation errors) {
        return interpolateWQ(km, DEFAULT_Q_STEPS, errors);
    }


    /**
     * Interpolate W and Q values at a given km.
     */
    public double [][] interpolateWQ(
            final double km,
            final int steps,
            final Calculation errors
            ) {
        int rowIndex = Collections.binarySearch(this.rows, new Row(km));

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

        rowIndex = -rowIndex -1;

        if (rowIndex < 1 || rowIndex >= this.rows.size()) {
            // do not extrapolate
            if (errors != null) {
                errors.addProblem(km, "km.not.found");
            }
            return new double[2][0];
        }

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

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

    public boolean interpolate(
            final double    km,
            final double [] out,
            final QPosition qPosition,
            final Function  qFunction
            ) {
        final int R1 = this.rows.size()-1;

        out[1] = qFunction.value(getQ(qPosition, km));

        if (Double.isNaN(out[1])) {
            return false;
        }

        final QPosition nPosition = new QPosition();
        if (getQPosition(km, out[1], nPosition) == null) {
            return false;
        }

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

        if (rowIndex >= 0) {
            // direct row match
            out[0] = this.rows.get(rowIndex).getW(nPosition);
            return !Double.isNaN(out[0]);
        }

        rowIndex = -rowIndex -1;

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

        final Row r1 = this.rows.get(rowIndex-1);
        final Row r2 = this.rows.get(rowIndex);
        out[0] = r1.getW(r2, km, nPosition);

        return !Double.isNaN(out[0]);
    }


    /**
     * Look up interpolation of a Q at given positions.
     *
     * @param q           the non-interpolated Q value.
     * @param referenceKm the reference km (e.g. gauge position).
     * @param kms         positions for which to interpolate.
     * @param ws          (output) resulting interpolated ws.
     * @param qs          (output) resulting interpolated qs.
     * @param errors      calculation object to store errors.
     */
    public QPosition interpolate(
            final double      q,
            final double      referenceKm,
            final double []   kms,
            final double []   ws,
            final double []   qs,
            final Calculation errors
            ) {
        return interpolate(
                q, referenceKm, kms, ws, qs, 0, kms.length, errors);
    }

    /**
     * Interpolate Q at given positions.
     * @param kms positions for which to calculate qs and ws
     * @param ws [out] calculated ws for kms
     * @param qs [out] looked up qs for kms.
     */
    public QPosition interpolate(
            final double      q,
            final double      referenceKm,
            final double []   kms,
            final double []   ws,
            final double []   qs,
            final int         startIndex,
            final int         length,
            final Calculation errors
            ) {
        final QPosition qPosition = getQPosition(referenceKm, q);

        if (qPosition == null) {
            // we cannot locate q at km
            Arrays.fill(ws, Double.NaN);
            Arrays.fill(qs, Double.NaN);
            if (errors != null) {
                errors.addProblem(referenceKm, "cannot.find.q", q);
            }
            return null;
        }

        final Row kmKey = new Row();

        final int R1 = this.rows.size()-1;

        for (int i = startIndex, end = startIndex+length; i < end; ++i) {

            if (Double.isNaN(qs[i] = getQ(qPosition, kms[i]))) {
                if (errors != null) {
                    errors.addProblem(kms[i], "cannot.find.q", q);
                }
                ws[i] = Double.NaN;
                continue;
            }

            kmKey.km = kms[i];
            int rowIndex = Collections.binarySearch(this.rows, kmKey);

            if (rowIndex >= 0) {
                // direct row match
                if (Double.isNaN(ws[i] = this.rows.get(rowIndex).getW(qPosition))
                        && errors != null) {
                    errors.addProblem(kms[i], "cannot.find.w.for.q", q);
                }
                continue;
            }

            rowIndex = -rowIndex -1;

            if (rowIndex < 1 || rowIndex > R1) {
                // do not extrapolate
                if (errors != null) {
                    errors.addProblem(kms[i], "km.not.found");
                }
                ws[i] = Double.NaN;
                continue;
            }
            final Row r1 = this.rows.get(rowIndex-1);
            final Row r2 = this.rows.get(rowIndex);

            if (Double.isNaN(ws[i] = r1.getW(r2, kms[i], qPosition))
                    && errors != null) {
                errors.addProblem(kms[i], "cannot.find.w.for.q", q);
            }
        }

        return qPosition;
    }

    /**
     * Linearly interpolate w at a km at a column of two rows.
     *
     * @param km   position for which to interpolate.
     * @param row1 first row.
     * @param row2 second row.
     * @param col  column-index at which to look.
     *
     * @return Linearly interpolated w, NaN if one of the given rows was null.
     */
    public static double linearW(final double km, final Row row1, final Row row2, final int col) {
        if (row1 == null || row2 == null) {
            return Double.NaN;
        }

        return Linear.linear(km,
                row1.km, row2.km,
                row1.ws[col], row2.ws[col]);
    }

    /**
     * Do interpolation/lookup of W and Q within columns (i.e. ignoring values
     * of other columns).
     * @param km position (km) at which to interpolate/lookup.
     * @return [[q0, q1, .. qx] , [w0, w1, .. wx]] (can contain NaNs)
     */
    public double [][] interpolateWQColumnwise(final double km) {
        log.debug("WstValueTable.interpolateWQColumnwise");
        final double [] qs = new double[this.columns.length];
        final double [] ws = new double[this.columns.length];

        // Find out row from where we will start searching.
        int rowIndex = Collections.binarySearch(this.rows, new Row(km));

        if (rowIndex < 0) {
            rowIndex = -rowIndex -1;
        }

        // TODO Beyond definition, we could stop more clever.
        if (rowIndex >= this.rows.size()) {
            rowIndex = this.rows.size() -1;
        }

        final Row startRow = this.rows.get(rowIndex);

        for (int col = 0; col < this.columns.length; col++) {
            qs[col] = this.columns[col].getQRangeTree().findQ(km);
            if (startRow.km == km && !Double.isNaN(startRow.ws[col])) {
                // Great. W is defined at km.
                ws[col] = startRow.ws[col];
                continue;
            }

            // Search neighbouring rows that define w at this col.
            Row rowBefore = null;
            Row rowAfter  = null;
            for (int before = rowIndex -1; before >= 0; before--) {
                if (!Double.isNaN(this.rows.get(before).ws[col])) {
                    rowBefore = this.rows.get(before);
                    break;
                }
            }
            if (rowBefore != null) {
                for (int after = rowIndex, R = this.rows.size();
                        after < R;
                        after++
                        ) {
                    if (!Double.isNaN(this.rows.get(after).ws[col])) {
                        rowAfter = this.rows.get(after);
                        break;
                    }
                }
            }

            ws[col] = linearW(km, rowBefore, rowAfter, col);
        }

        return new double [][] {qs, ws};
    }

    public double [] findQsForW(final double km, final double w, final Calculation errors) {

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

        if (rowIndex >= 0) {
            return this.rows.get(rowIndex).findQsForW(w, this, errors);
        }

        rowIndex = -rowIndex - 1;

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

        // Needs bilinear interpolation.
        final Row r1 = this.rows.get(rowIndex-1);
        final Row r2 = this.rows.get(rowIndex);

        return r1.findQsForW(r2, w, km, this, errors);
    }

    protected SplineFunction createSpline(final double km, final Calculation errors) {

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

        if (rowIndex >= 0) {
            final SplineFunction sf = this.rows.get(rowIndex).createSpline(this, errors);
            if (sf == null && errors != null) {
                errors.addProblem(km, "cannot.create.wq.relation");
            }
            return sf;
        }

        rowIndex = -rowIndex - 1;

        if (rowIndex < 1 || rowIndex >= this.rows.size()) {
            // Do not extrapolate.
            if (errors != null) {
                errors.addProblem(km, "km.not.found");
            }
            return null;
        }

        // Needs bilinear interpolation.
        final Row r1 = this.rows.get(rowIndex-1);
        final Row r2 = this.rows.get(rowIndex);

        final SplineFunction sf = r1.createSpline(r2, km, this, errors);
        if (sf == null && errors != null) {
            errors.addProblem(km, "cannot.create.wq.relation");
        }

        return sf;
    }

    /** 'Bezugslinienverfahren' */
    public double [][] relateWs(
            final double      km1,
            final double      km2,
            final Calculation errors
            ) {
        return relateWs(km1, km2, RELATE_WS_SAMPLES, errors);
    }

    private static class ErrorHandler {

        boolean     hasErrors;
        Calculation errors;

        ErrorHandler(final Calculation errors) {
            this.errors = errors;
        }

        void error(final double km, final String key, final Object ... args) {
            if (this.errors != null && !this.hasErrors) {
                this.hasErrors = true;
                this.errors.addProblem(km, key, args);
            }
        }
    } // class ErrorHandler


    /* TODO: Add optimized methods of relateWs to relate one
     *       start km to many end kms. The index generation/spline stuff for
     *       the start km is always the same.
     */
    public double [][] relateWs(
            final double      km1,
            final double      km2,
            final int         numSamples,
            final Calculation errors
            ) {
        final SplineFunction sf1 = createSpline(km1, errors);
        if (sf1 == null) {
            return new double[2][0];
        }

        final SplineFunction sf2 = createSpline(km2, errors);
        if (sf2 == null) {
            return new double[2][0];
        }

        final PolynomialSplineFunction iQ1 = sf1.createIndexQRelation();
        if (iQ1 == null) {
            if (errors != null) {
                errors.addProblem(km1, "cannot.create.index.q.relation");
            }
            return new double[2][0];
        }

        final PolynomialSplineFunction iQ2 = sf2.createIndexQRelation();
        if (iQ2 == null) {
            if (errors != null) {
                errors.addProblem(km2, "cannot.create.index.q.relation");
            }
            return new double[2][0];
        }

        final int N = Math.min(sf1.splineQs.length, sf2.splineQs.length);
        final double stepWidth = N/(double)numSamples;

        final PolynomialSplineFunction qW1 = sf1.spline;
        final PolynomialSplineFunction qW2 = sf2.spline;

        final TDoubleArrayList ws1 = new TDoubleArrayList(numSamples);
        final TDoubleArrayList ws2 = new TDoubleArrayList(numSamples);
        final TDoubleArrayList qs1 = new TDoubleArrayList(numSamples);
        final TDoubleArrayList qs2 = new TDoubleArrayList(numSamples);

        final ErrorHandler err = new ErrorHandler(errors);

        int i = 0;
        for (double p = 0d; p <= N-1; p += stepWidth, ++i) {

            double q1;
            try {
                q1 = iQ1.value(p);
            }
            catch (final ArgumentOutsideDomainException aode) {
                err.error(km1, "w.w.qkm1.failed", p);
                continue;
            }

            double w1;
            try {
                w1 = qW1.value(q1);
            }
            catch (final ArgumentOutsideDomainException aode) {
                err.error(km1, "w.w.wkm1.failed", q1, p);
                continue;
            }

            double q2;
            try {
                q2 = iQ2.value(p);
            }
            catch (final ArgumentOutsideDomainException aode) {
                err.error(km2, "w.w.qkm2.failed", p);
                continue;
            }

            double w2;
            try {
                w2 = qW2.value(q2);
            }
            catch (final ArgumentOutsideDomainException aode) {
                err.error(km2, "w.w.wkm2.failed", q2, p);
                continue;
            }

            ws1.add(w1);
            ws2.add(w2);
            qs1.add(q1);
            qs2.add(q2);
        }

        return new double [][] {
            ws1.toNativeArray(),
            qs1.toNativeArray(),
            ws2.toNativeArray(),
            qs2.toNativeArray() };
    }

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

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

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

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

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

        for (int i = 1; i < this.columns.length; ++i) {
            final double qCurrent = this.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) {
                final double weight = Linear.factor(q, qLast, qCurrent);
                return qPosition.set(i, weight);
            }
            qLast = qCurrent;
        }

        return null;
    }

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

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

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

    public double [][] interpolateTabulated(final double km) {
        return interpolateTabulated(km, new double[2][this.columns.length]);
    }

    public double [][] interpolateTabulated(final double km, final double [][] result) {

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

        if (rowIndex >= 0) {
            // Direct hit -> copy ws.
            final Row row = this.rows.get(rowIndex);
            System.arraycopy(
                    row.ws, 0, result[0], 0,
                    Math.min(row.ws.length, result[0].length));
        }
        else {
            rowIndex = -rowIndex -1;
            if (rowIndex < 1 || rowIndex >= this.rows.size()) {
                // Out of bounds.
                return null;
            }
            // Interpolate ws.
            final Row r1 = this.rows.get(rowIndex-1);
            final Row r2 = this.rows.get(rowIndex);
            final double factor = Linear.factor(km, r1.km, r2.km);
            Linear.weight(factor, r1.ws, r2.ws, result[0]);
        }

        final double [] qs = result[1];
        for (int i = Math.min(qs.length, this.columns.length)-1; i >= 0; --i) {
            qs[i] = this.columns[i].getQRangeTree().findQ(km);
        }
        return result;
    }


    /** True if no QRange is given or Q equals zero. */
    public boolean hasEmptyQ() {
        for (final Column column: this.columns) {
            if (column.getQRangeTree() == null) {
                return true;
            }
            else {
                if (Math.abs(column.getQRangeTree().maxQ()) <= 0.01d) {
                    return true;
                }
            }
        }

        if (this.columns.length == 0) {
            log.warn("No columns in WstValueTable.");
        }

        return false;
    }


    /** Find ranges that are between km1 and km2 (inclusive?) */
    public List<Range> findSegments(final double km1, final double km2) {
        return this.columns.length != 0
                ? this.columns[this.columns.length-1].getQRangeTree().findSegments(km1, km2)
                        : Collections.<Range>emptyList();
    }
}
// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :