teichmann@5863: /* Copyright (C) 2011, 2012, 2013 by Bundesanstalt für Gewässerkunde
teichmann@5863:  * Software engineering by Intevation GmbH
teichmann@5863:  *
teichmann@5994:  * This file is Free Software under the GNU AGPL (>=v3)
teichmann@5863:  * and comes with ABSOLUTELY NO WARRANTY! Check out the
teichmann@5994:  * documentation coming with Dive4Elements River for details.
teichmann@5863:  */
teichmann@5863: 
teichmann@5831: package org.dive4elements.river.artifacts.model;
sascha@326: 
mschaefer@9286: import static org.dive4elements.river.backend.utils.EpsilonComparator.CMP;
sascha@326: 
mschaefer@9286: import java.io.Serializable;
mschaefer@9286: import java.util.ArrayList;
mschaefer@9286: import java.util.Arrays;
mschaefer@9286: import java.util.Collections;
mschaefer@9286: import java.util.List;
mschaefer@9286: 
mschaefer@9286: import org.apache.commons.math.ArgumentOutsideDomainException;
mschaefer@9286: import org.apache.commons.math.analysis.interpolation.SplineInterpolator;
mschaefer@9286: import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
mschaefer@9286: import org.apache.commons.math.exception.MathIllegalArgumentException;
mschaefer@9286: import org.apache.log4j.Logger;
mschaefer@9286: import org.dive4elements.river.artifacts.math.Function;
teichmann@5831: import org.dive4elements.river.artifacts.math.Linear;
teichmann@6434: import org.dive4elements.river.utils.DoubleUtil;
sascha@427: 
sascha@1937: import gnu.trove.TDoubleArrayList;
sascha@1937: 
felix@1838: /**
felix@1838:  * W, Q and km data from database 'wsts' spiced with interpolation algorithms.
felix@1838:  */
sascha@326: public class WstValueTable
sascha@326: implements   Serializable
sascha@326: {
sascha@339:     private static Logger log = Logger.getLogger(WstValueTable.class);
sascha@326: 
sascha@395:     public static final int DEFAULT_Q_STEPS = 500;
sascha@339: 
sascha@2253:     public static final int RELATE_WS_SAMPLES = 200;
sascha@2186: 
tom@8664:     public static final double W_EPSILON = 0.000001;
tom@8664: 
felix@1721:     /**
felix@1721:      * A Column in the table, typically representing one measurement session.
felix@1721:      */
sascha@633:     public static final class Column
sascha@633:     implements                Serializable
sascha@326:     {
sascha@326:         protected String name;
sascha@326: 
sascha@632:         protected QRangeTree qRangeTree;
sascha@632: 
sascha@326:         public Column() {
sascha@326:         }
sascha@326: 
mschaefer@9286:         public Column(final String name) {
sascha@326:             this.name = name;
sascha@326:         }
sascha@326: 
sascha@326:         public String getName() {
mschaefer@9286:             return this.name;
sascha@326:         }
sascha@326: 
mschaefer@9286:         public void setName(final String name) {
sascha@326:             this.name = name;
sascha@326:         }
sascha@632: 
teichmann@4797:         public QRangeTree getQRangeTree() {
mschaefer@9286:             return this.qRangeTree;
sascha@632:         }
sascha@632: 
mschaefer@9286:         public void setQRangeTree(final QRangeTree qRangeTree) {
sascha@632:             this.qRangeTree = qRangeTree;
sascha@632:         }
sascha@633:     } // class Column
sascha@326: 
felix@1721:     /**
felix@1721:      * A (weighted) position used for interpolation.
felix@1721:      */
sascha@633:     public static final class QPosition {
sascha@426: 
sascha@426:         protected int    index;
sascha@426:         protected double weight;
sascha@426: 
sascha@426:         public QPosition() {
sascha@426:         }
sascha@426: 
mschaefer@9286:         public QPosition(final int index, final double weight) {
sascha@426:             this.index  = index;
sascha@426:             this.weight = weight;
sascha@426:         }
sascha@426: 
mschaefer@9286:         public QPosition set(final int index, final double weight) {
sascha@633:             this.index  = index;
sascha@633:             this.weight = weight;
sascha@633:             return this;
sascha@633:         }
sascha@633: 
sascha@426:     } // class Position
sascha@426: 
sascha@1939:     public static final class SplineFunction {
sascha@1939: 
sascha@1939:         public PolynomialSplineFunction spline;
sascha@1939:         public double []                splineQs;
sascha@1939:         public double []                splineWs;
sascha@1939: 
sascha@1939:         public SplineFunction(
mschaefer@9286:                 final PolynomialSplineFunction spline,
mschaefer@9286:                 final double []                splineQs,
mschaefer@9286:                 final double []                splineWs
mschaefer@9286:                 ) {
sascha@1939:             this.spline   = spline;
sascha@1939:             this.splineQs = splineQs;
sascha@1939:             this.splineWs = splineWs;
sascha@1939:         }
sascha@1939: 
sascha@1939:         public double [][] sample(
mschaefer@9286:                 final int         numSamples,
mschaefer@9286:                 final double      km,
mschaefer@9286:                 final Calculation errors
mschaefer@9286:                 ) {
mschaefer@9286:             final double minQ = getQMin();
mschaefer@9286:             final double maxQ = getQMax();
sascha@1939: 
mschaefer@9286:             final double [] outWs = new double[numSamples];
mschaefer@9286:             final double [] outQs = new double[numSamples];
sascha@1939: 
sascha@1939:             Arrays.fill(outWs, Double.NaN);
sascha@1939:             Arrays.fill(outQs, Double.NaN);
sascha@1939: 
mschaefer@9286:             final double stepWidth = (maxQ - minQ)/numSamples;
sascha@1939: 
sascha@1939:             try {
sascha@1939:                 double q = minQ;
sascha@1939:                 for (int i = 0; i < outWs.length; ++i, q += stepWidth) {
mschaefer@9286:                     outWs[i] = this.spline.value(outQs[i] = q);
sascha@1939:                 }
sascha@1939:             }
mschaefer@9286:             catch (final ArgumentOutsideDomainException aode) {
sascha@1939:                 if (errors != null) {
sascha@2166:                     errors.addProblem(km, "spline.interpolation.failed");
sascha@1939:                 }
sascha@1939:                 log.error("spline interpolation failed.", aode);
sascha@1939:             }
sascha@1939: 
sascha@1939:             return new double [][] { outWs, outQs };
sascha@1939:         }
sascha@1939: 
sascha@1939:         public double getQMin() {
mschaefer@9286:             return Math.min(this.splineQs[0], this.splineQs[this.splineQs.length-1]);
sascha@1939:         }
sascha@1939: 
sascha@1939:         public double getQMax() {
mschaefer@9286:             return Math.max(this.splineQs[0], this.splineQs[this.splineQs.length-1]);
sascha@1939:         }
sascha@1939: 
sascha@1939:         /** Constructs a continues index between the columns to Qs. */
sascha@1939:         public PolynomialSplineFunction createIndexQRelation() {
sascha@1939: 
mschaefer@9286:             final double [] indices = new double[this.splineQs.length];
sascha@1939:             for (int i = 0; i < indices.length; ++i) {
sascha@1939:                 indices[i] = i;
sascha@1939:             }
sascha@1939: 
sascha@1939:             try {
mschaefer@9286:                 final SplineInterpolator interpolator = new SplineInterpolator();
mschaefer@9286:                 return interpolator.interpolate(indices, this.splineQs);
sascha@1939:             }
mschaefer@9286:             catch (final MathIllegalArgumentException miae) {
sascha@1939:                 // Ignore me!
sascha@1939:             }
sascha@1939:             return null;
sascha@1939:         }
sascha@1939:     } // class SplineFunction
sascha@1939: 
felix@1721:     /**
felix@1838:      * A row, typically a position where measurements were taken.
felix@1721:      */
sascha@633:     public static final class Row
sascha@633:     implements                Serializable, Comparable<Row>
sascha@326:     {
sascha@326:         double    km;
sascha@335:         double [] ws;
sascha@326: 
sascha@326:         public Row() {
sascha@326:         }
sascha@326: 
mschaefer@9286:         public Row(final double km) {
sascha@326:             this.km = km;
sascha@335:         }
sascha@335: 
mschaefer@9286:         public Row(final double km, final double [] ws) {
sascha@335:             this(km);
sascha@335:             this.ws = ws;
sascha@335:         }
sascha@335: 
felix@1838:         /**
tom@8704:          * Sort Qs and Ws for this Row over Q.
tom@8704:          */
tom@8704:         private double[][] getSortedWQ(
mschaefer@9286:                 final WstValueTable table,
mschaefer@9286:                 final Calculation   errors
mschaefer@9286:                 ) {
mschaefer@9286:             final int W = this.ws.length;
tom@8704: 
tom@8704:             if (W < 1) {
tom@8704:                 if (errors != null) {
mschaefer@9286:                     errors.addProblem(this.km, "no.ws.found");
tom@8704:                 }
tom@8704:                 return new double[][] {{Double.NaN}, {Double.NaN}};
tom@8704:             }
tom@8704: 
mschaefer@9286:             final double [] sortedWs = this.ws;
mschaefer@9286:             final double [] sortedQs = new double[W];
tom@8704: 
tom@8704:             for (int i=0; i < W; ++i) {
mschaefer@9286:                 final double q = table.getQIndex(i, this.km);
tom@8704:                 if (Double.isNaN(q) && errors != null) {
tom@8704:                     errors.addProblem(
mschaefer@9286:                             this.km, "no.q.found.in.column", i+1);
tom@8704:                 }
tom@8704:                 sortedQs[i] = q;
tom@8704:             }
tom@8704: 
tom@8704:             DoubleUtil.sortByFirst(sortedQs, sortedWs);
tom@8704: 
tom@8704:             return new double[][] { sortedWs, sortedQs };
tom@8704:         }
tom@8704: 
tom@8704:         /**
tom@8704:          * Return an array of Qs and Ws for the given km between
tom@8704:          * this Row and another, sorted over Q.
tom@8704:          */
tom@8704:         private double[][] getSortedWQ(
mschaefer@9286:                 final Row           other,
mschaefer@9286:                 final double        km,
mschaefer@9286:                 final WstValueTable table,
mschaefer@9286:                 final Calculation   errors
mschaefer@9286:                 ) {
mschaefer@9286:             final int W = Math.min(this.ws.length, other.ws.length);
tom@8704: 
tom@8704:             if (W < 1) {
tom@8704:                 if (errors != null) {
tom@8704:                     errors.addProblem("no.ws.found");
tom@8704:                 }
tom@8704:                 return new double[][] {{Double.NaN}, {Double.NaN}};
tom@8704:             }
tom@8704: 
mschaefer@9286:             final double factor = Linear.factor(km, this.km, other.km);
tom@8704: 
mschaefer@9286:             final double [] sortedQs = new double[W];
mschaefer@9286:             final double [] sortedWs = new double[W];
tom@8704: 
tom@8704:             for (int i = 0; i < W; ++i) {
mschaefer@9286:                 final double wws = Linear.weight(factor, this.ws[i], other.ws[i]);
mschaefer@9286:                 final double wqs = table.getQIndex(i, km);
tom@8704: 
tom@8704:                 if (Double.isNaN(wws) || Double.isNaN(wqs)) {
tom@8704:                     if (errors != null) {
tom@8704:                         errors.addProblem(km, "cannot.find.w.or.q");
tom@8704:                     }
tom@8704:                 }
tom@8704: 
tom@8704:                 sortedWs[i] = wws;
tom@8704:                 sortedQs[i] = wqs;
tom@8704:             }
tom@8704: 
tom@8704:             DoubleUtil.sortByFirst(sortedQs, sortedWs);
tom@8704: 
tom@8704:             return new double[][] { sortedWs, sortedQs };
tom@8704:         }
tom@8704: 
tom@8704:         /**
tom@8704:          * Find Qs matching w in an array of Qs and Ws sorted over Q.
tom@8704:          */
mschaefer@9286:         private double[] findQsForW(final double w, final double[][] sortedWQ) {
mschaefer@9286:             final int W = sortedWQ[0].length;
tom@8704: 
mschaefer@9286:             final double[] sortedQs = sortedWQ[1];
mschaefer@9286:             final double[] sortedWs = sortedWQ[0];
tom@8704: 
mschaefer@9286:             final TDoubleArrayList qs = new TDoubleArrayList();
tom@8704: 
tom@8704:             if (W > 0 && Math.abs(sortedWs[0]-w) < W_EPSILON) {
mschaefer@9286:                 final double q = sortedQs[0];
tom@8704:                 if (!Double.isNaN(q)) {
tom@8704:                     qs.add(q);
tom@8704:                 }
tom@8704:             }
tom@8704: 
tom@8704:             for (int i = 1; i < W; ++i) {
mschaefer@9286:                 final double w2 = sortedWs[i];
tom@8704:                 if (Double.isNaN(w2)) {
tom@8704:                     continue;
tom@8704:                 }
tom@8704:                 if (Math.abs(w2-w) < W_EPSILON) {
mschaefer@9286:                     final double q = sortedQs[i];
tom@8704:                     if (!Double.isNaN(q)) {
tom@8704:                         qs.add(q);
tom@8704:                     }
tom@8704:                     continue;
tom@8704:                 }
mschaefer@9286:                 final double w1 = sortedWs[i-1];
tom@8704:                 if (Double.isNaN(w1)) {
tom@8704:                     continue;
tom@8704:                 }
tom@8704: 
tom@8704:                 if (w < Math.min(w1, w2) || w > Math.max(w1, w2)) {
tom@8704:                     continue;
tom@8704:                 }
tom@8704: 
mschaefer@9286:                 final double q1 = sortedQs[i-1];
mschaefer@9286:                 final double q2 = sortedQs[i];
tom@8704:                 if (Double.isNaN(q1) || Double.isNaN(q2)) {
tom@8704:                     continue;
tom@8704:                 }
tom@8704: 
mschaefer@9286:                 final double q = Linear.linear(w, w1, w2, q1, q2);
tom@8704:                 qs.add(q);
tom@8704:             }
tom@8704: 
tom@8704:             return qs.toNativeArray();
tom@8704:         }
tom@8704: 
tom@8704:         /**
felix@1838:          * Compare according to place of measurement (km).
felix@1838:          */
mschaefer@9286:         @Override
mschaefer@9286:         public int compareTo(final Row other) {
mschaefer@9286:             return CMP.compare(this.km, other.km);
sascha@326:         }
sascha@339: 
felix@1887:         /**
felix@1887:          * Interpolate Ws, given Qs and a km.
felix@1887:          *
felix@1887:          * @param iqs Given ("input") Qs.
felix@1887:          * @param ows Resulting ("output") Ws.
felix@1887:          * @param table Table of which to use data for interpolation.
felix@1887:          */
sascha@633:         public void interpolateW(
mschaefer@9286:                 final Row           other,
mschaefer@9286:                 final double        km,
mschaefer@9286:                 final double []     iqs,
mschaefer@9286:                 final double []     ows,
mschaefer@9286:                 final WstValueTable table,
mschaefer@9286:                 final Calculation   errors
mschaefer@9286:                 ) {
mschaefer@9286:             final double kmWeight = Linear.factor(km, this.km, other.km);
sascha@339: 
mschaefer@9286:             final QPosition qPosition = new QPosition();
sascha@633: 
sascha@633:             for (int i = 0; i < iqs.length; ++i) {
sascha@633:                 if (table.getQPosition(km, iqs[i], qPosition) != null) {
mschaefer@9286:                     final double wt =       getW(qPosition);
mschaefer@9286:                     final double wo = other.getW(qPosition);
ingo@686:                     if (Double.isNaN(wt) || Double.isNaN(wo)) {
ingo@686:                         if (errors != null) {
ingo@686:                             errors.addProblem(
mschaefer@9286:                                     km, "cannot.find.w.for.q", iqs[i]);
ingo@686:                         }
ingo@686:                         ows[i] = Double.NaN;
ingo@686:                     }
ingo@686:                     else {
ingo@686:                         ows[i] = Linear.weight(kmWeight, wt, wo);
ingo@686:                     }
sascha@633:                 }
sascha@633:                 else {
ingo@686:                     if (errors != null) {
sascha@2166:                         errors.addProblem(km, "cannot.find.q", iqs[i]);
ingo@686:                     }
sascha@633:                     ows[i] = Double.NaN;
sascha@633:                 }
sascha@633:             }
sascha@633:         }
sascha@633: 
sascha@1938: 
sascha@1938:         public SplineFunction createSpline(
mschaefer@9286:                 final WstValueTable table,
mschaefer@9286:                 final Calculation   errors
mschaefer@9286:                 ) {
mschaefer@9286:             final double[][] sortedWQ = getSortedWQ(table, errors);
teichmann@7000: 
sascha@1938:             try {
mschaefer@9286:                 final SplineInterpolator interpolator = new SplineInterpolator();
mschaefer@9286:                 final PolynomialSplineFunction spline =
mschaefer@9286:                         interpolator.interpolate(sortedWQ[1], sortedWQ[0]);
sascha@1938: 
mschaefer@9286:                 return new SplineFunction(spline, sortedWQ[1], this.ws);
sascha@1938:             }
mschaefer@9286:             catch (final MathIllegalArgumentException miae) {
sascha@1938:                 if (errors != null) {
mschaefer@9286:                     errors.addProblem(this.km, "spline.creation.failed");
sascha@1938:                 }
sascha@1938:                 log.error("spline creation failed", miae);
sascha@1938:             }
sascha@1938:             return null;
sascha@1938:         }
sascha@1938: 
sascha@1938:         public SplineFunction createSpline(
mschaefer@9286:                 final Row           other,
mschaefer@9286:                 final double        km,
mschaefer@9286:                 final WstValueTable table,
mschaefer@9286:                 final Calculation   errors
mschaefer@9286:                 ) {
mschaefer@9286:             final double[][] sortedWQ = getSortedWQ(other, km, table, errors);
teichmann@6434: 
mschaefer@9286:             final SplineInterpolator interpolator = new SplineInterpolator();
ingo@686: 
ingo@686:             try {
mschaefer@9286:                 final PolynomialSplineFunction spline =
mschaefer@9286:                         interpolator.interpolate(sortedWQ[1], sortedWQ[0]);
sascha@1938: 
tom@8704:                 return new SplineFunction(spline, sortedWQ[1], sortedWQ[0]);
ingo@686:             }
mschaefer@9286:             catch (final MathIllegalArgumentException miae) {
ingo@686:                 if (errors != null) {
sascha@2166:                     errors.addProblem(km, "spline.creation.failed");
ingo@686:                 }
felix@1860:                 log.error("spline creation failed", miae);
ingo@686:             }
sascha@339: 
sascha@1938:             return null;
sascha@1938:         }
ingo@686: 
sascha@1938:         public double [][] interpolateWQ(
mschaefer@9286:                 final Row           other,
mschaefer@9286:                 final double        km,
mschaefer@9286:                 final int           steps,
mschaefer@9286:                 final WstValueTable table,
mschaefer@9286:                 final Calculation   errors
mschaefer@9286:                 ) {
mschaefer@9286:             final SplineFunction sf = createSpline(other, km, table, errors);
sascha@339: 
sascha@1938:             return sf != null
mschaefer@9286:                     ? sf.sample(steps, km, errors)
mschaefer@9286:                             : new double[2][0];
sascha@339:         }
sascha@339: 
felix@1721: 
ingo@686:         public double [][] interpolateWQ(
mschaefer@9286:                 final int           steps,
mschaefer@9286:                 final WstValueTable table,
mschaefer@9286:                 final Calculation   errors
mschaefer@9286:                 ) {
mschaefer@9286:             final SplineFunction sf = createSpline(table, errors);
sascha@339: 
sascha@1938:             return sf != null
mschaefer@9286:                     ? sf.sample(steps, this.km, errors)
mschaefer@9286:                             : new double[2][0];
sascha@339:         }
sascha@339: 
sascha@426: 
mschaefer@9286:         public double getW(final QPosition qPosition) {
mschaefer@9286:             final int    index  = qPosition.index;
mschaefer@9286:             final double weight = qPosition.weight;
sascha@426: 
sascha@633:             return weight == 1.0
mschaefer@9286:                     ? this.ws[index]
mschaefer@9286:                             : Linear.weight(weight, this.ws[index-1], this.ws[index]);
sascha@426:         }
sascha@426: 
sascha@633:         public double getW(
mschaefer@9286:                 final Row       other,
mschaefer@9286:                 final double    km,
mschaefer@9286:                 final QPosition qPosition
mschaefer@9286:                 ) {
mschaefer@9286:             final double kmWeight = Linear.factor(km, this.km, other.km);
sascha@426: 
mschaefer@9286:             final int    index  = qPosition.index;
mschaefer@9286:             final double weight = qPosition.weight;
sascha@426: 
sascha@633:             double tw, ow;
sascha@426: 
sascha@633:             if (weight == 1.0) {
mschaefer@9286:                 tw = this.ws[index];
sascha@633:                 ow = other.ws[index];
sascha@426:             }
sascha@426:             else {
mschaefer@9286:                 tw = Linear.weight(weight, this.ws[index-1], this.ws[index]);
sascha@742:                 ow = Linear.weight(weight, other.ws[index-1], other.ws[index]);
sascha@426:             }
sascha@426: 
sascha@655:             return Linear.weight(kmWeight, tw, ow);
sascha@427:         }
sascha@1937: 
tom@8704:         public double [] findQsForW(
mschaefer@9286:                 final double        w,
mschaefer@9286:                 final WstValueTable table,
mschaefer@9286:                 final Calculation   errors
mschaefer@9286:                 ) {
mschaefer@9286:             log.debug("Find Qs for given W at tabulated km " + this.km);
tom@8704:             return findQsForW(w, getSortedWQ(table, errors));
sascha@1937:         }
sascha@1937: 
sascha@1937:         public double [] findQsForW(
mschaefer@9286:                 final Row           other,
mschaefer@9286:                 final double        w,
mschaefer@9286:                 final double        km,
mschaefer@9286:                 final WstValueTable table,
mschaefer@9286:                 final Calculation   errors
mschaefer@9286:                 ) {
tom@8704:             log.debug("Find Qs for given W at non-tabulated km " + km);
tom@8704:             return findQsForW(w, getSortedWQ(other, km, table, errors));
sascha@1937:         }
ingo@2423: 
mschaefer@9286:         public double [] getMinMaxW(final double [] result) {
ingo@2423:             double minW =  Double.MAX_VALUE;
ingo@2423:             double maxW = -Double.MAX_VALUE;
mschaefer@9286:             for (int i = 0; i < this.ws.length; ++i) {
mschaefer@9286:                 final double w = this.ws[i];
ingo@2423:                 if (w < minW) minW = w;
ingo@2423:                 if (w > maxW) maxW = w;
ingo@2423:             }
ingo@2423:             result[0] = minW;
ingo@2423:             result[1] = maxW;
ingo@2423:             return result;
ingo@2423:         }
ingo@2423: 
mschaefer@9286:         public double [] getMinMaxW(final Row other, final double km, final double [] result) {
mschaefer@9286:             final double [] m1 = this .getMinMaxW(new double [2]);
mschaefer@9286:             final double [] m2 = other.getMinMaxW(new double [2]);
mschaefer@9286:             final double factor = Linear.factor(km, this.km, other.km);
ingo@2423:             result[0] = Linear.weight(factor, m1[0], m2[0]);
ingo@2423:             result[1] = Linear.weight(factor, m1[1], m2[1]);
ingo@2423:             return result;
ingo@2423:         }
sascha@633:     } // class Row
sascha@326: 
felix@1838:     /** Rows in table. */
sascha@326:     protected List<Row> rows;
sascha@326: 
felix@1838:     /** Columns in table. */
sascha@326:     protected Column [] columns;
sascha@326: 
sascha@326:     public WstValueTable() {
mschaefer@9286:         this.rows = new ArrayList<>();
sascha@326:     }
sascha@326: 
mschaefer@9286:     public WstValueTable(final Column [] columns) {
sascha@326:         this();
sascha@326:         this.columns = columns;
sascha@326:     }
sascha@326: 
tom@8663:     /**
tom@8663:      * @param columns The WST-columns.
tom@8663:      * @param rows A list of Rows that must be sorted by km.
tom@8663:      */
mschaefer@9286:     public WstValueTable(final Column [] columns, final List<Row> rows) {
sascha@633:         this.columns = columns;
sascha@633:         this.rows    = rows;
sascha@633:     }
sascha@633: 
teichmann@4797:     public Column [] getColumns() {
mschaefer@9286:         return this.columns;
teichmann@4797:     }
teichmann@4797: 
felix@1887:     /**
felix@1887:      * @param km Given kilometer.
felix@1887:      * @param qs Given Q values.
felix@1887:      * @param ws output parameter.
felix@1887:      */
mschaefer@9286:     public double [] interpolateW(final double km, final double [] qs, final double [] ws) {
ingo@686:         return interpolateW(km, qs, ws, null);
ingo@686:     }
sascha@335: 
felix@1098: 
felix@1098:     /**
felix@1098:      * @param ws (output parameter), gets returned.
felix@1098:      * @return output parameter ws.
felix@1098:      */
ingo@686:     public double [] interpolateW(
mschaefer@9286:             final double      km,
mschaefer@9286:             final double []   qs,
mschaefer@9286:             final double []   ws,
mschaefer@9286:             final Calculation errors
mschaefer@9286:             ) {
mschaefer@9286:         int rowIndex = Collections.binarySearch(this.rows, new Row(km));
sascha@335: 
mschaefer@9286:         final QPosition qPosition = new QPosition();
sascha@633: 
sascha@335:         if (rowIndex >= 0) { // direct row match
mschaefer@9286:             final Row row = this.rows.get(rowIndex);
sascha@335:             for (int i = 0; i < qs.length; ++i) {
ingo@686:                 if (getQPosition(km, qs[i], qPosition) == null) {
ingo@686:                     if (errors != null) {
sascha@2166:                         errors.addProblem(km, "cannot.find.q", qs[i]);
ingo@686:                     }
ingo@686:                     ws[i] = Double.NaN;
ingo@686:                 }
ingo@686:                 else {
ingo@686:                     if (Double.isNaN(ws[i] = row.getW(qPosition))
mschaefer@9286:                             && errors != null) {
ingo@686:                         errors.addProblem(
mschaefer@9286:                                 km, "cannot.find.w.for.q", qs[i]);
ingo@686:                     }
ingo@686:                 }
sascha@335:             }
sascha@335:         }
sascha@335:         else { // needs bilinear interpolation
sascha@335:             rowIndex = -rowIndex -1;
sascha@335: 
mschaefer@9286:             if (rowIndex < 1 || rowIndex >= this.rows.size()) {
sascha@335:                 // do not extrapolate
sascha@335:                 Arrays.fill(ws, Double.NaN);
ingo@686:                 if (errors != null) {
sascha@2166:                     errors.addProblem(km, "km.not.found");
ingo@686:                 }
sascha@335:             }
sascha@335:             else {
mschaefer@9286:                 final Row r1 = this.rows.get(rowIndex-1);
mschaefer@9286:                 final Row r2 = this.rows.get(rowIndex);
ingo@686:                 r1.interpolateW(r2, km, qs, ws, this, errors);
sascha@335:             }
sascha@335:         }
sascha@335: 
sascha@335:         return ws;
sascha@335:     }
sascha@335: 
mschaefer@9286:     /**
mschaefer@9286:      * Interpolates a W for a km using a Q column position
mschaefer@9286:      */
mschaefer@9286:     public double interpolateW(final double km, final QPosition qPosition, final Calculation errors) {
mschaefer@9286: 
mschaefer@9286:         int rowIndex = Collections.binarySearch(this.rows, new Row(km));
mschaefer@9286: 
mschaefer@9286:         if (rowIndex >= 0) { // direct row match
mschaefer@9286:             final Row row = this.rows.get(rowIndex);
mschaefer@9286:             return row.getW(qPosition);
mschaefer@9286:         }
mschaefer@9286:         else { // needs bilinear interpolation
mschaefer@9286:             rowIndex = -rowIndex - 1;
mschaefer@9286:             if ((rowIndex <= 0) || (rowIndex >= this.rows.size()))
mschaefer@9286:                 return Double.NaN;
mschaefer@9286:             else {
mschaefer@9286:                 final Row r1 = this.rows.get(rowIndex - 1);
mschaefer@9286:                 final Row r2 = this.rows.get(rowIndex);
mschaefer@9286:                 final double w1 = r1.getW(qPosition);
mschaefer@9286:                 final double w2 = r2.getW(qPosition);
mschaefer@9286:                 if (Double.isNaN(w1) || Double.isNaN(w2))
mschaefer@9286:                     return Double.NaN;
mschaefer@9286:                 else {
mschaefer@9286:                     final double kmWeight = Linear.factor(km, r1.km, r2.km);
mschaefer@9286:                     return Linear.weight(kmWeight, w1, w2);
mschaefer@9286:                 }
mschaefer@9286:             }
mschaefer@9286:         }
mschaefer@9286:     }
mschaefer@9286: 
mschaefer@9286:     public double [] getMinMaxQ(final double km) {
sascha@2559:         return getMinMaxQ(km, new double [2]);
sascha@2559:     }
sascha@2559: 
mschaefer@9286:     public double [] getMinMaxQ(final double km, final double [] result) {
sascha@2559:         double minQ =  Double.MAX_VALUE;
sascha@2559:         double maxQ = -Double.MAX_VALUE;
sascha@2559: 
mschaefer@9286:         for (int i = 0; i < this.columns.length; ++i) {
mschaefer@9286:             final double q = this.columns[i].getQRangeTree().findQ(km);
sascha@2559:             if (!Double.isNaN(q)) {
sascha@2559:                 if (q < minQ) minQ = q;
sascha@2559:                 if (q > maxQ) maxQ = q;
sascha@2559:             }
sascha@2559:         }
sascha@2559: 
sascha@2559:         if (minQ < Double.MAX_VALUE) {
sascha@2559:             result[0] = minQ;
sascha@2559:             result[1] = maxQ;
sascha@2559:             return result;
sascha@2559:         }
sascha@2559: 
sascha@2559:         return null;
sascha@2559:     }
sascha@2559: 
sascha@2560:     public double [] getMinMaxQ(double from, double to, double step) {
mschaefer@9286:         final double [] result = new double[2];
sascha@2560: 
sascha@2560:         double minQ =  Double.MAX_VALUE;
sascha@2560:         double maxQ = -Double.MAX_VALUE;
sascha@2560: 
sascha@2560:         if (from > to) {
mschaefer@9286:             final double tmp = from;
sascha@2560:             from = to;
sascha@2560:             to = tmp;
sascha@2560:         }
sascha@2560: 
sascha@2560:         step = Math.max(Math.abs(step), 0.0001);
sascha@2560: 
sascha@2560:         double d = from;
sascha@2560:         for (; d <= to; d += step) {
sascha@2560:             if (getMinMaxQ(d, result) != null) {
sascha@2560:                 if (result[0] < minQ) minQ = result[0];
sascha@2560:                 if (result[1] > maxQ) maxQ = result[1];
sascha@2560:             }
sascha@2560:         }
sascha@2560: 
sascha@2560:         if (d != to) {
sascha@2560:             if (getMinMaxQ(to, result) != null) {
sascha@2560:                 if (result[0] < minQ) minQ = result[0];
sascha@2560:                 if (result[1] > maxQ) maxQ = result[1];
sascha@2560:             }
sascha@2560:         }
sascha@2560: 
sascha@2560:         return minQ < Double.MAX_VALUE
mschaefer@9286:                 ? new double [] { minQ, maxQ }
mschaefer@9286:         : null;
sascha@2560:     }
sascha@2559: 
mschaefer@9286:     public double [] getMinMaxW(final double km) {
ingo@2423:         return getMinMaxW(km, new double [2]);
ingo@2423: 
ingo@2423:     }
mschaefer@9286:     public double [] getMinMaxW(final double km, final double [] result) {
mschaefer@9286:         int rowIndex = Collections.binarySearch(this.rows, new Row(km));
ingo@2423: 
ingo@2423:         if (rowIndex >= 0) {
mschaefer@9286:             return this.rows.get(rowIndex).getMinMaxW(result);
ingo@2423:         }
ingo@2423: 
ingo@2423:         rowIndex = -rowIndex -1;
ingo@2423: 
mschaefer@9286:         if (rowIndex < 1 || rowIndex >= this.rows.size()) {
ingo@2423:             // do not extrapolate
ingo@2423:             return null;
ingo@2423:         }
ingo@2423: 
mschaefer@9286:         final Row r1 = this.rows.get(rowIndex-1);
mschaefer@9286:         final Row r2 = this.rows.get(rowIndex);
ingo@2423: 
ingo@2423:         return r1.getMinMaxW(r2, km, result);
ingo@2423:     }
ingo@2423: 
ingo@2423:     public double [] getMinMaxW(double from, double to, double step) {
mschaefer@9286:         final double [] result = new double[2];
ingo@2423: 
ingo@2423:         double minW =  Double.MAX_VALUE;
ingo@2423:         double maxW = -Double.MAX_VALUE;
ingo@2423: 
ingo@2423:         if (from > to) {
mschaefer@9286:             final double tmp = from;
ingo@2423:             from = to;
ingo@2423:             to = tmp;
ingo@2423:         }
ingo@2423: 
ingo@2423:         step = Math.max(Math.abs(step), 0.0001);
ingo@2423: 
ingo@2423:         double d = from;
ingo@2423:         for (; d <= to; d += step) {
ingo@2423:             if (getMinMaxW(d, result) != null) {
ingo@2423:                 if (result[0] < minW) minW = result[0];
ingo@2423:                 if (result[1] > maxW) maxW = result[1];
ingo@2423:             }
ingo@2423:         }
ingo@2423: 
ingo@2423:         if (d != to) {
ingo@2423:             if (getMinMaxW(to, result) != null) {
ingo@2423:                 if (result[0] < minW) minW = result[0];
ingo@2423:                 if (result[1] > maxW) maxW = result[1];
ingo@2423:             }
ingo@2423:         }
ingo@2423: 
ingo@2423:         return minW < Double.MAX_VALUE
mschaefer@9286:                 ? new double [] { minW, maxW }
mschaefer@9286:         : null;
ingo@2423:     }
ingo@2423: 
felix@1838:     /**
felix@1838:      * Interpolate W and Q values at a given km.
felix@1838:      */
mschaefer@9286:     public double [][] interpolateWQ(final double km) {
ingo@686:         return interpolateWQ(km, null);
sascha@339:     }
sascha@339: 
felix@1838:     /**
felix@1838:      * Interpolate W and Q values at a given km.
felix@3265:      *
felix@3265:      * @param errors where to store errors.
felix@3265:      *
felix@3265:      * @return double double array, first index Ws, second Qs.
felix@1838:      */
mschaefer@9286:     public double [][] interpolateWQ(final double km, final Calculation errors) {
ingo@686:         return interpolateWQ(km, DEFAULT_Q_STEPS, errors);
ingo@686:     }
ingo@686: 
ingo@686: 
felix@1860:     /**
felix@1860:      * Interpolate W and Q values at a given km.
felix@1860:      */
tom@8856:     public double [][] interpolateWQ(
mschaefer@9286:             final double km,
mschaefer@9286:             final int steps,
mschaefer@9286:             final Calculation errors
mschaefer@9286:             ) {
mschaefer@9286:         int rowIndex = Collections.binarySearch(this.rows, new Row(km));
sascha@339: 
sascha@339:         if (rowIndex >= 0) { // direct row match
mschaefer@9286:             final Row row = this.rows.get(rowIndex);
ingo@686:             return row.interpolateWQ(steps, this, errors);
sascha@339:         }
sascha@339: 
sascha@339:         rowIndex = -rowIndex -1;
sascha@339: 
mschaefer@9286:         if (rowIndex < 1 || rowIndex >= this.rows.size()) {
sascha@339:             // do not extrapolate
ingo@686:             if (errors != null) {
sascha@2166:                 errors.addProblem(km, "km.not.found");
ingo@686:             }
sascha@339:             return new double[2][0];
sascha@339:         }
sascha@339: 
mschaefer@9286:         final Row r1 = this.rows.get(rowIndex-1);
mschaefer@9286:         final Row r2 = this.rows.get(rowIndex);
sascha@339: 
ingo@686:         return r1.interpolateWQ(r2, km, steps, this, errors);
sascha@339:     }
sascha@339: 
sascha@655:     public boolean interpolate(
mschaefer@9286:             final double    km,
mschaefer@9286:             final double [] out,
mschaefer@9286:             final QPosition qPosition,
mschaefer@9286:             final Function  qFunction
mschaefer@9286:             ) {
mschaefer@9286:         final int R1 = this.rows.size()-1;
sascha@427: 
sascha@655:         out[1] = qFunction.value(getQ(qPosition, km));
sascha@655: 
sascha@655:         if (Double.isNaN(out[1])) {
sascha@655:             return false;
sascha@655:         }
sascha@633: 
mschaefer@9286:         final QPosition nPosition = new QPosition();
sascha@655:         if (getQPosition(km, out[1], nPosition) == null) {
sascha@655:             return false;
sascha@655:         }
sascha@427: 
mschaefer@9286:         int rowIndex = Collections.binarySearch(this.rows, new Row(km));
sascha@427: 
sascha@655:         if (rowIndex >= 0) {
sascha@655:             // direct row match
mschaefer@9286:             out[0] = this.rows.get(rowIndex).getW(nPosition);
sascha@655:             return !Double.isNaN(out[0]);
sascha@427:         }
sascha@655: 
sascha@655:         rowIndex = -rowIndex -1;
sascha@655: 
sascha@1681:         if (rowIndex < 1 || rowIndex > R1) {
sascha@655:             // do not extrapolate
sascha@655:             return false;
sascha@655:         }
sascha@655: 
mschaefer@9286:         final Row r1 = this.rows.get(rowIndex-1);
mschaefer@9286:         final Row r2 = this.rows.get(rowIndex);
sascha@655:         out[0] = r1.getW(r2, km, nPosition);
sascha@655: 
sascha@655:         return !Double.isNaN(out[0]);
sascha@427:     }
sascha@427: 
felix@1098: 
felix@1098:     /**
felix@1098:      * Look up interpolation of a Q at given positions.
felix@1098:      *
felix@1098:      * @param q           the non-interpolated Q value.
felix@1098:      * @param referenceKm the reference km (e.g. gauge position).
felix@1098:      * @param kms         positions for which to interpolate.
felix@1098:      * @param ws          (output) resulting interpolated ws.
felix@1098:      * @param qs          (output) resulting interpolated qs.
felix@1098:      * @param errors      calculation object to store errors.
felix@1098:      */
sascha@427:     public QPosition interpolate(
mschaefer@9286:             final double      q,
mschaefer@9286:             final double      referenceKm,
mschaefer@9286:             final double []   kms,
mschaefer@9286:             final double []   ws,
mschaefer@9286:             final double []   qs,
mschaefer@9286:             final Calculation errors
mschaefer@9286:             ) {
ingo@686:         return interpolate(
mschaefer@9286:                 q, referenceKm, kms, ws, qs, 0, kms.length, errors);
sascha@649:     }
sascha@649: 
felix@6711:     /**
felix@6711:      * Interpolate Q at given positions.
felix@6711:      * @param kms positions for which to calculate qs and ws
felix@6711:      * @param ws [out] calculated ws for kms
felix@6711:      * @param qs [out] looked up qs for kms.
felix@6711:      */
sascha@649:     public QPosition interpolate(
mschaefer@9286:             final double      q,
mschaefer@9286:             final double      referenceKm,
mschaefer@9286:             final double []   kms,
mschaefer@9286:             final double []   ws,
mschaefer@9286:             final double []   qs,
mschaefer@9286:             final int         startIndex,
mschaefer@9286:             final int         length,
mschaefer@9286:             final Calculation errors
mschaefer@9286:             ) {
mschaefer@9286:         final QPosition qPosition = getQPosition(referenceKm, q);
sascha@426: 
sascha@426:         if (qPosition == null) {
sascha@633:             // we cannot locate q at km
ingo@686:             Arrays.fill(ws, Double.NaN);
ingo@686:             Arrays.fill(qs, Double.NaN);
ingo@686:             if (errors != null) {
sascha@2166:                 errors.addProblem(referenceKm, "cannot.find.q", q);
ingo@686:             }
sascha@427:             return null;
sascha@426:         }
sascha@426: 
mschaefer@9286:         final Row kmKey = new Row();
sascha@645: 
mschaefer@9286:         final int R1 = this.rows.size()-1;
sascha@645: 
sascha@649:         for (int i = startIndex, end = startIndex+length; i < end; ++i) {
sascha@458: 
ingo@686:             if (Double.isNaN(qs[i] = getQ(qPosition, kms[i]))) {
ingo@686:                 if (errors != null) {
sascha@2166:                     errors.addProblem(kms[i], "cannot.find.q", q);
ingo@686:                 }
ingo@686:                 ws[i] = Double.NaN;
ingo@686:                 continue;
ingo@686:             }
ingo@686: 
ingo@686:             kmKey.km = kms[i];
mschaefer@9286:             int rowIndex = Collections.binarySearch(this.rows, kmKey);
sascha@426: 
sascha@426:             if (rowIndex >= 0) {
sascha@426:                 // direct row match
mschaefer@9286:                 if (Double.isNaN(ws[i] = this.rows.get(rowIndex).getW(qPosition))
mschaefer@9286:                         && errors != null) {
sascha@2166:                     errors.addProblem(kms[i], "cannot.find.w.for.q", q);
ingo@686:                 }
sascha@426:                 continue;
sascha@426:             }
sascha@426: 
sascha@426:             rowIndex = -rowIndex -1;
sascha@426: 
sascha@1681:             if (rowIndex < 1 || rowIndex > R1) {
sascha@426:                 // do not extrapolate
ingo@686:                 if (errors != null) {
sascha@2166:                     errors.addProblem(kms[i], "km.not.found");
ingo@686:                 }
felix@1888:                 ws[i] = Double.NaN;
sascha@426:                 continue;
sascha@426:             }
mschaefer@9286:             final Row r1 = this.rows.get(rowIndex-1);
mschaefer@9286:             final Row r2 = this.rows.get(rowIndex);
sascha@633: 
ingo@686:             if (Double.isNaN(ws[i] = r1.getW(r2, kms[i], qPosition))
mschaefer@9286:                     && errors != null) {
sascha@2166:                 errors.addProblem(kms[i], "cannot.find.w.for.q", q);
ingo@686:             }
sascha@426:         }
sascha@426: 
sascha@427:         return qPosition;
sascha@426:     }
sascha@633: 
felix@1888:     /**
felix@1888:      * Linearly interpolate w at a km at a column of two rows.
felix@1888:      *
felix@1888:      * @param km   position for which to interpolate.
felix@1888:      * @param row1 first row.
felix@1888:      * @param row2 second row.
felix@1888:      * @param col  column-index at which to look.
felix@1888:      *
felix@1888:      * @return Linearly interpolated w, NaN if one of the given rows was null.
felix@1888:      */
mschaefer@9286:     public static double linearW(final double km, final Row row1, final Row row2, final int col) {
felix@1888:         if (row1 == null || row2 == null) {
felix@1888:             return Double.NaN;
felix@1888:         }
felix@1888: 
felix@1888:         return Linear.linear(km,
mschaefer@9286:                 row1.km, row2.km,
mschaefer@9286:                 row1.ws[col], row2.ws[col]);
felix@1888:     }
felix@1888: 
felix@1888:     /**
felix@1888:      * Do interpolation/lookup of W and Q within columns (i.e. ignoring values
felix@1888:      * of other columns).
felix@1888:      * @param km position (km) at which to interpolate/lookup.
felix@1888:      * @return [[q0, q1, .. qx] , [w0, w1, .. wx]] (can contain NaNs)
felix@1888:      */
mschaefer@9286:     public double [][] interpolateWQColumnwise(final double km) {
felix@1919:         log.debug("WstValueTable.interpolateWQColumnwise");
mschaefer@9286:         final double [] qs = new double[this.columns.length];
mschaefer@9286:         final double [] ws = new double[this.columns.length];
felix@1888: 
felix@1888:         // Find out row from where we will start searching.
mschaefer@9286:         int rowIndex = Collections.binarySearch(this.rows, new Row(km));
felix@1888: 
felix@1888:         if (rowIndex < 0) {
felix@1888:             rowIndex = -rowIndex -1;
felix@1888:         }
felix@2248: 
felix@2248:         // TODO Beyond definition, we could stop more clever.
mschaefer@9286:         if (rowIndex >= this.rows.size()) {
mschaefer@9286:             rowIndex = this.rows.size() -1;
sascha@1937:         }
felix@1888: 
mschaefer@9286:         final Row startRow = this.rows.get(rowIndex);
felix@1888: 
mschaefer@9286:         for (int col = 0; col < this.columns.length; col++) {
mschaefer@9286:             qs[col] = this.columns[col].getQRangeTree().findQ(km);
teichmann@7015:             if (startRow.km == km && !Double.isNaN(startRow.ws[col])) {
felix@1888:                 // Great. W is defined at km.
felix@1888:                 ws[col] = startRow.ws[col];
felix@1888:                 continue;
felix@1888:             }
felix@1888: 
felix@1888:             // Search neighbouring rows that define w at this col.
felix@1888:             Row rowBefore = null;
felix@1888:             Row rowAfter  = null;
felix@1899:             for (int before = rowIndex -1; before >= 0; before--) {
mschaefer@9286:                 if (!Double.isNaN(this.rows.get(before).ws[col])) {
mschaefer@9286:                     rowBefore = this.rows.get(before);
felix@1888:                     break;
felix@1888:                 }
felix@1888:             }
felix@2248:             if (rowBefore != null) {
mschaefer@9286:                 for (int after = rowIndex, R = this.rows.size();
mschaefer@9286:                         after < R;
mschaefer@9286:                         after++
mschaefer@9286:                         ) {
mschaefer@9286:                     if (!Double.isNaN(this.rows.get(after).ws[col])) {
mschaefer@9286:                         rowAfter = this.rows.get(after);
felix@2248:                         break;
felix@2248:                     }
felix@1888:                 }
felix@1888:             }
felix@1888: 
felix@1888:             ws[col] = linearW(km, rowBefore, rowAfter, col);
felix@1888:         }
felix@1888: 
felix@1888:         return new double [][] {qs, ws};
felix@1888:     }
felix@1888: 
mschaefer@9286:     public double [] findQsForW(final double km, final double w, final Calculation errors) {
sascha@1937: 
mschaefer@9286:         int rowIndex = Collections.binarySearch(this.rows, new Row(km));
sascha@1937: 
sascha@1937:         if (rowIndex >= 0) {
mschaefer@9286:             return this.rows.get(rowIndex).findQsForW(w, this, errors);
sascha@1937:         }
sascha@1937: 
sascha@1937:         rowIndex = -rowIndex - 1;
sascha@1937: 
mschaefer@9286:         if (rowIndex < 1 || rowIndex >= this.rows.size()) {
sascha@1939:             // Do not extrapolate.
sascha@1937:             return new double[0];
sascha@1937:         }
sascha@1937: 
sascha@1939:         // Needs bilinear interpolation.
mschaefer@9286:         final Row r1 = this.rows.get(rowIndex-1);
mschaefer@9286:         final Row r2 = this.rows.get(rowIndex);
sascha@1937: 
tom@8704:         return r1.findQsForW(r2, w, km, this, errors);
sascha@1937:     }
sascha@1937: 
mschaefer@9286:     protected SplineFunction createSpline(final double km, final Calculation errors) {
sascha@1939: 
mschaefer@9286:         int rowIndex = Collections.binarySearch(this.rows, new Row(km));
sascha@1939: 
sascha@1939:         if (rowIndex >= 0) {
mschaefer@9286:             final SplineFunction sf = this.rows.get(rowIndex).createSpline(this, errors);
sascha@1939:             if (sf == null && errors != null) {
sascha@2166:                 errors.addProblem(km, "cannot.create.wq.relation");
sascha@1939:             }
sascha@1939:             return sf;
sascha@1939:         }
sascha@1939: 
sascha@1939:         rowIndex = -rowIndex - 1;
sascha@1939: 
mschaefer@9286:         if (rowIndex < 1 || rowIndex >= this.rows.size()) {
sascha@1939:             // Do not extrapolate.
sascha@1939:             if (errors != null) {
sascha@2166:                 errors.addProblem(km, "km.not.found");
sascha@1939:             }
sascha@1939:             return null;
sascha@1939:         }
sascha@1939: 
sascha@1939:         // Needs bilinear interpolation.
mschaefer@9286:         final Row r1 = this.rows.get(rowIndex-1);
mschaefer@9286:         final Row r2 = this.rows.get(rowIndex);
sascha@1939: 
mschaefer@9286:         final SplineFunction sf = r1.createSpline(r2, km, this, errors);
sascha@1939:         if (sf == null && errors != null) {
sascha@2166:             errors.addProblem(km, "cannot.create.wq.relation");
sascha@1939:         }
sascha@1939: 
sascha@1939:         return sf;
sascha@1939:     }
sascha@1939: 
sascha@1939:     /** 'Bezugslinienverfahren' */
sascha@1939:     public double [][] relateWs(
mschaefer@9286:             final double      km1,
mschaefer@9286:             final double      km2,
mschaefer@9286:             final Calculation errors
mschaefer@9286:             ) {
sascha@2186:         return relateWs(km1, km2, RELATE_WS_SAMPLES, errors);
sascha@2186:     }
sascha@2186: 
ingo@2330:     private static class ErrorHandler {
ingo@2330: 
ingo@2330:         boolean     hasErrors;
ingo@2330:         Calculation errors;
ingo@2330: 
mschaefer@9286:         ErrorHandler(final Calculation errors) {
ingo@2330:             this.errors = errors;
ingo@2330:         }
ingo@2330: 
mschaefer@9286:         void error(final double km, final String key, final Object ... args) {
mschaefer@9286:             if (this.errors != null && !this.hasErrors) {
mschaefer@9286:                 this.hasErrors = true;
mschaefer@9286:                 this.errors.addProblem(km, key, args);
ingo@2330:             }
ingo@2330:         }
ingo@2330:     } // class ErrorHandler
ingo@2330: 
ingo@2330: 
sascha@2186:     /* TODO: Add optimized methods of relateWs to relate one
sascha@3076:      *       start km to many end kms. The index generation/spline stuff for
sascha@2186:      *       the start km is always the same.
sascha@2186:      */
sascha@2186:     public double [][] relateWs(
mschaefer@9286:             final double      km1,
mschaefer@9286:             final double      km2,
mschaefer@9286:             final int         numSamples,
mschaefer@9286:             final Calculation errors
mschaefer@9286:             ) {
mschaefer@9286:         final SplineFunction sf1 = createSpline(km1, errors);
sascha@1939:         if (sf1 == null) {
sascha@1939:             return new double[2][0];
sascha@1939:         }
sascha@1939: 
mschaefer@9286:         final SplineFunction sf2 = createSpline(km2, errors);
sascha@1939:         if (sf2 == null) {
sascha@1939:             return new double[2][0];
sascha@1939:         }
sascha@1939: 
mschaefer@9286:         final PolynomialSplineFunction iQ1 = sf1.createIndexQRelation();
sascha@1939:         if (iQ1 == null) {
sascha@1939:             if (errors != null) {
sascha@2166:                 errors.addProblem(km1, "cannot.create.index.q.relation");
sascha@1939:             }
sascha@1939:             return new double[2][0];
sascha@1939:         }
sascha@1939: 
mschaefer@9286:         final PolynomialSplineFunction iQ2 = sf2.createIndexQRelation();
sascha@1939:         if (iQ2 == null) {
sascha@1939:             if (errors != null) {
sascha@2166:                 errors.addProblem(km2, "cannot.create.index.q.relation");
sascha@1939:             }
sascha@1939:             return new double[2][0];
sascha@1939:         }
sascha@1939: 
mschaefer@9286:         final int N = Math.min(sf1.splineQs.length, sf2.splineQs.length);
mschaefer@9286:         final double stepWidth = N/(double)numSamples;
sascha@1939: 
mschaefer@9286:         final PolynomialSplineFunction qW1 = sf1.spline;
mschaefer@9286:         final PolynomialSplineFunction qW2 = sf2.spline;
sascha@1939: 
mschaefer@9286:         final TDoubleArrayList ws1 = new TDoubleArrayList(numSamples);
mschaefer@9286:         final TDoubleArrayList ws2 = new TDoubleArrayList(numSamples);
mschaefer@9286:         final TDoubleArrayList qs1 = new TDoubleArrayList(numSamples);
mschaefer@9286:         final TDoubleArrayList qs2 = new TDoubleArrayList(numSamples);
mschaefer@9286: 
mschaefer@9286:         final ErrorHandler err = new ErrorHandler(errors);
sascha@1939: 
sascha@2316:         int i = 0;
sascha@2316:         for (double p = 0d; p <= N-1; p += stepWidth, ++i) {
ingo@2330: 
ingo@2330:             double q1;
sascha@1939:             try {
ingo@2330:                 q1 = iQ1.value(p);
sascha@1939:             }
mschaefer@9286:             catch (final ArgumentOutsideDomainException aode) {
ingo@2330:                 err.error(km1, "w.w.qkm1.failed", p);
ingo@2330:                 continue;
sascha@1939:             }
ingo@2330: 
ingo@2330:             double w1;
ingo@2330:             try {
ingo@2330:                 w1 = qW1.value(q1);
ingo@2330:             }
mschaefer@9286:             catch (final ArgumentOutsideDomainException aode) {
sascha@2403:                 err.error(km1, "w.w.wkm1.failed", q1, p);
ingo@2330:                 continue;
ingo@2330:             }
ingo@2330: 
ingo@2330:             double q2;
ingo@2330:             try {
ingo@2330:                 q2 = iQ2.value(p);
ingo@2330:             }
mschaefer@9286:             catch (final ArgumentOutsideDomainException aode) {
ingo@2330:                 err.error(km2, "w.w.qkm2.failed", p);
ingo@2330:                 continue;
ingo@2330:             }
ingo@2330: 
ingo@2330:             double w2;
ingo@2330:             try {
ingo@2330:                 w2 = qW2.value(q2);
ingo@2330:             }
mschaefer@9286:             catch (final ArgumentOutsideDomainException aode) {
sascha@2403:                 err.error(km2, "w.w.wkm2.failed", q2, p);
ingo@2330:                 continue;
ingo@2330:             }
ingo@2330: 
ingo@2330:             ws1.add(w1);
ingo@2330:             ws2.add(w2);
ingo@2330:             qs1.add(q1);
ingo@2330:             qs2.add(q2);
sascha@1939:         }
sascha@1939: 
sascha@2316:         return new double [][] {
sascha@2316:             ws1.toNativeArray(),
sascha@3076:             qs1.toNativeArray(),
sascha@2316:             ws2.toNativeArray(),
sascha@2316:             qs2.toNativeArray() };
sascha@1939:     }
sascha@1939: 
mschaefer@9286:     public QPosition getQPosition(final double km, final double q) {
sascha@633:         return getQPosition(km, q, new QPosition());
sascha@633:     }
sascha@633: 
mschaefer@9286:     public QPosition getQPosition(final double km, final double q, final QPosition qPosition) {
sascha@633: 
mschaefer@9286:         if (this.columns.length == 0) {
sascha@633:             return null;
sascha@633:         }
sascha@633: 
mschaefer@9286:         double qLast = this.columns[0].getQRangeTree().findQ(km);
sascha@633: 
sascha@633:         if (Math.abs(qLast - q) < 0.00001) {
sascha@633:             return qPosition.set(0, 1d);
sascha@633:         }
sascha@633: 
mschaefer@9286:         for (int i = 1; i < this.columns.length; ++i) {
mschaefer@9286:             final double qCurrent = this.columns[i].getQRangeTree().findQ(km);
sascha@633:             if (Math.abs(qCurrent - q) < 0.00001) {
sascha@633:                 return qPosition.set(i, 1d);
sascha@633:             }
sascha@633: 
sascha@633:             double qMin, qMax;
sascha@633:             if (qLast < qCurrent) { qMin = qLast; qMax = qCurrent; }
sascha@633:             else                  { qMin = qCurrent; qMax = qLast; }
sascha@633: 
sascha@633:             if (q > qMin && q < qMax) {
mschaefer@9286:                 final double weight = Linear.factor(q, qLast, qCurrent);
sascha@633:                 return qPosition.set(i, weight);
sascha@633:             }
sascha@633:             qLast = qCurrent;
sascha@633:         }
sascha@633: 
sascha@633:         return null;
sascha@633:     }
sascha@633: 
mschaefer@9286:     public double getQIndex(final int index, final double km) {
mschaefer@9286:         return this.columns[index].getQRangeTree().findQ(km);
sascha@633:     }
sascha@633: 
mschaefer@9286:     public double getQ(final QPosition qPosition, final double km) {
mschaefer@9286:         final int    index  = qPosition.index;
mschaefer@9286:         final double weight = qPosition.weight;
sascha@633: 
sascha@633:         if (weight == 1d) {
mschaefer@9286:             return this.columns[index].getQRangeTree().findQ(km);
sascha@633:         }
mschaefer@9286:         final double q1 = this.columns[index-1].getQRangeTree().findQ(km);
mschaefer@9286:         final double q2 = this.columns[index  ].getQRangeTree().findQ(km);
sascha@655:         return Linear.weight(weight, q1, q2);
sascha@633:     }
sascha@3736: 
mschaefer@9286:     public double [][] interpolateTabulated(final double km) {
mschaefer@9286:         return interpolateTabulated(km, new double[2][this.columns.length]);
sascha@3736:     }
sascha@3736: 
mschaefer@9286:     public double [][] interpolateTabulated(final double km, final double [][] result) {
sascha@3736: 
mschaefer@9286:         int rowIndex = Collections.binarySearch(this.rows, new Row(km));
sascha@3736: 
sascha@3736:         if (rowIndex >= 0) {
sascha@3736:             // Direct hit -> copy ws.
mschaefer@9286:             final Row row = this.rows.get(rowIndex);
sascha@3736:             System.arraycopy(
mschaefer@9286:                     row.ws, 0, result[0], 0,
mschaefer@9286:                     Math.min(row.ws.length, result[0].length));
sascha@3736:         }
sascha@3736:         else {
sascha@3736:             rowIndex = -rowIndex -1;
mschaefer@9286:             if (rowIndex < 1 || rowIndex >= this.rows.size()) {
sascha@3736:                 // Out of bounds.
sascha@3736:                 return null;
sascha@3736:             }
sascha@3736:             // Interpolate ws.
mschaefer@9286:             final Row r1 = this.rows.get(rowIndex-1);
mschaefer@9286:             final Row r2 = this.rows.get(rowIndex);
mschaefer@9286:             final double factor = Linear.factor(km, r1.km, r2.km);
sascha@3736:             Linear.weight(factor, r1.ws, r2.ws, result[0]);
sascha@3736:         }
sascha@3736: 
mschaefer@9286:         final double [] qs = result[1];
mschaefer@9286:         for (int i = Math.min(qs.length, this.columns.length)-1; i >= 0; --i) {
mschaefer@9286:             qs[i] = this.columns[i].getQRangeTree().findQ(km);
sascha@3736:         }
sascha@3736:         return result;
sascha@3736:     }
sascha@3743: 
felix@4119: 
felix@6923:     /** True if no QRange is given or Q equals zero. */
felix@6923:     public boolean hasEmptyQ() {
mschaefer@9286:         for (final Column column: this.columns) {
felix@6923:             if (column.getQRangeTree() == null) {
felix@6923:                 return true;
felix@6923:             }
felix@6923:             else {
felix@6923:                 if (Math.abs(column.getQRangeTree().maxQ()) <= 0.01d) {
felix@6923:                     return true;
felix@6923:                 }
felix@6923:             }
felix@6923:         }
felix@6923: 
mschaefer@9286:         if (this.columns.length == 0) {
felix@6923:             log.warn("No columns in WstValueTable.");
felix@6923:         }
felix@6923: 
felix@6923:         return false;
felix@6923:     }
felix@6923: 
felix@6923: 
felix@4119:     /** Find ranges that are between km1 and km2 (inclusive?) */
mschaefer@9286:     public List<Range> findSegments(final double km1, final double km2) {
mschaefer@9286:         return this.columns.length != 0
mschaefer@9286:                 ? this.columns[this.columns.length-1].getQRangeTree().findSegments(km1, km2)
mschaefer@9286:                         : Collections.<Range>emptyList();
sascha@3743:     }
sascha@326: }
sascha@326: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :