sascha@326: package de.intevation.flys.artifacts.model;
sascha@326: 
sascha@326: import java.io.Serializable;
sascha@326: 
sascha@655: import de.intevation.flys.artifacts.math.Linear;
sascha@655: import de.intevation.flys.artifacts.math.Function;
sascha@427: 
sascha@335: import java.util.Arrays;
sascha@326: import java.util.ArrayList;
sascha@326: import java.util.List;
sascha@326: import java.util.Collections;
sascha@326: 
sascha@339: import org.apache.log4j.Logger;
sascha@339: 
sascha@339: import org.apache.commons.math.analysis.interpolation.SplineInterpolator;
sascha@339: 
sascha@339: import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
sascha@339: 
sascha@339: import org.apache.commons.math.ArgumentOutsideDomainException;
sascha@339: 
ingo@686: import org.apache.commons.math.exception.MathIllegalArgumentException;
ingo@686: 
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: 
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: 
sascha@326:         public Column(String name) {
sascha@326:             this.name = name;
sascha@326:         }
sascha@326: 
sascha@326:         public String getName() {
sascha@326:             return name;
sascha@326:         }
sascha@326: 
sascha@326:         public void setName(String name) {
sascha@326:             this.name = name;
sascha@326:         }
sascha@632: 
sascha@632:         public QRangeTree getQRangeTree() {
sascha@632:             return qRangeTree;
sascha@632:         }
sascha@632: 
sascha@632:         public void setQRangeTree(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: 
sascha@633:         public QPosition(int index, double weight) {
sascha@426:             this.index  = index;
sascha@426:             this.weight = weight;
sascha@426:         }
sascha@426: 
sascha@633:         public QPosition set(int index, 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: 
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: 
sascha@335:         public Row(double km) {
sascha@326:             this.km = km;
sascha@335:         }
sascha@335: 
sascha@633:         public Row(double km, double [] ws) {
sascha@335:             this(km);
sascha@335:             this.ws = ws;
sascha@335:         }
sascha@335: 
felix@1838:         /**
felix@1838:          * Compare according to place of measurement (km).
felix@1838:          */
sascha@335:         public int compareTo(Row other) {
sascha@335:             double d = km - other.km;
sascha@458:             if (d < -0.0001) return -1;
sascha@458:             if (d >  0.0001) return +1;
sascha@335:             return 0;
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(
sascha@633:             Row           other,
sascha@633:             double        km,
sascha@633:             double []     iqs,
sascha@633:             double []     ows,
ingo@686:             WstValueTable table,
ingo@686:             Calculation   errors
sascha@633:         ) {
sascha@655:             double kmWeight = Linear.factor(km, this.km, other.km);
sascha@339: 
sascha@633:             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) {
sascha@633:                     double wt =       getW(qPosition);
sascha@633:                     double wo = other.getW(qPosition);
ingo@686:                     if (Double.isNaN(wt) || Double.isNaN(wo)) {
ingo@686:                         if (errors != null) {
ingo@686:                             // TODO: I18N
ingo@686:                             errors.addProblem(
ingo@686:                                 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) {
ingo@686:                         // TODO: I18N
ingo@686:                         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@633:         public double [][] interpolateWQ(
sascha@633:             Row           other,
sascha@633:             double        km,
sascha@633:             int           steps,
ingo@686:             WstValueTable table,
ingo@686:             Calculation   errors
sascha@633:         ) {
sascha@395:             int W = Math.min(ws.length, other.ws.length);
sascha@339: 
sascha@339:             if (W < 1) {
ingo@686:                 if (errors != null) {
ingo@686:                     // TODO: I18N
ingo@686:                     errors.addProblem("no ws found");
ingo@686:                 }
sascha@339:                 return new double[2][0];
sascha@339:             }
sascha@339: 
sascha@655:             double factor = Linear.factor(km, this.km, other.km);
sascha@339: 
sascha@395:             double [] splineQ = new double[W];
sascha@395:             double [] splineW = new double[W];
sascha@339: 
sascha@395:             double minQ =  Double.MAX_VALUE;
sascha@395:             double maxQ = -Double.MAX_VALUE;
sascha@339: 
sascha@339:             for (int i = 0; i < W; ++i) {
sascha@655:                 double wws = Linear.weight(factor, ws[i], other.ws[i]);
sascha@655:                 double wqs = Linear.weight(
sascha@633:                     factor,
sascha@633:                     table.getQIndex(i, km),
sascha@633:                     table.getQIndex(i, other.km));
sascha@633: 
ingo@686:                 if (Double.isNaN(wws) || Double.isNaN(wqs)) {
ingo@686:                     if (errors != null) {
ingo@686:                         // TODO: I18N
ingo@686:                         errors.addProblem(km, "cannot find w or q");
ingo@686:                     }
ingo@686:                 }
ingo@686:                 else {
ingo@686:                     if (wqs < minQ) minQ = wqs;
ingo@686:                     if (wqs > maxQ) maxQ = wqs;
ingo@686:                 }
ingo@686: 
sascha@395:                 splineW[i] = wws;
sascha@395:                 splineQ[i] = wqs;
sascha@339:             }
sascha@339: 
sascha@395:             double stepWidth = (maxQ - minQ)/steps;
sascha@395: 
sascha@339:             SplineInterpolator interpolator = new SplineInterpolator();
ingo@686:             PolynomialSplineFunction spline;
ingo@686: 
ingo@686:             try {
ingo@686:                 spline = interpolator.interpolate(splineQ, splineW);
ingo@686:             }
ingo@686:             catch (MathIllegalArgumentException miae) {
ingo@686:                 if (errors != null) {
ingo@686:                     // TODO: I18N
ingo@686:                     errors.addProblem(km, "spline creation failed");
ingo@686:                 }
felix@1860:                 log.error("spline creation failed", miae);
ingo@686:                 return new double[2][0];
ingo@686:             }
sascha@339: 
sascha@395:             double [] outWs = new double[steps];
sascha@395:             double [] outQs = new double[steps];
sascha@339: 
ingo@686:             Arrays.fill(outWs, Double.NaN);
ingo@686:             Arrays.fill(outQs, Double.NaN);
ingo@686: 
sascha@339:             try {
sascha@395:                 double q = minQ;
sascha@395:                 for (int i = 0; i < outWs.length; ++i, q += stepWidth) {
sascha@395:                     outWs[i] = spline.value(outQs[i] = q);
sascha@339:                 }
sascha@339:             }
sascha@339:             catch (ArgumentOutsideDomainException aode) {
ingo@686:                 if (errors != null) {
ingo@686:                     // TODO: I18N
ingo@686:                     errors.addProblem(km, "spline interpolation failed");
ingo@686:                 }
ingo@686:                 log.error("spline interpolation failed", aode);
sascha@339:             }
sascha@339: 
sascha@339:             return new double [][] { outWs, outQs };
sascha@339:         }
sascha@339: 
felix@1721: 
ingo@686:         public double [][] interpolateWQ(
sascha@742:             int           steps,
ingo@686:             WstValueTable table,
ingo@686:             Calculation   errors
ingo@686:         ) {
sascha@395:             int W = ws.length;
sascha@339: 
sascha@339:             if (W < 1) {
ingo@686:                 if (errors != null) {
ingo@686:                     // TODO: I18N
ingo@686:                     errors.addProblem(km, "no ws found");
ingo@686:                 }
sascha@339:                 return new double[2][0];
sascha@339:             }
sascha@339: 
sascha@395:             double [] splineQ = new double[W];
sascha@395: 
sascha@742:             double minQ =  Double.MAX_VALUE;
sascha@742:             double maxQ = -Double.MAX_VALUE;
sascha@339: 
sascha@339:             for (int i = 0; i < W; ++i) {
ingo@686:                 double sq = table.getQIndex(i, km);
ingo@686:                 if (Double.isNaN(sq)) {
ingo@686:                     if (errors != null) {
ingo@686:                         // TODO: I18N
ingo@686:                         errors.addProblem(
ingo@686:                             km, "no q found in " + (i+1) + " column");
ingo@686:                     }
ingo@686:                 }
ingo@686:                 else {
ingo@686:                     if (sq < minQ) minQ = sq;
ingo@686:                     if (sq > maxQ) maxQ = sq;
ingo@686:                 }
ingo@686:                 splineQ[i] = sq;
sascha@339:             }
sascha@339: 
sascha@395:             double stepWidth = (maxQ - minQ)/steps;
sascha@395: 
sascha@339:             SplineInterpolator interpolator = new SplineInterpolator();
sascha@339: 
ingo@686:             PolynomialSplineFunction spline;
sascha@742: 
ingo@686:             try {
ingo@686:                 spline = interpolator.interpolate(splineQ, ws);
ingo@686:             }
ingo@686:             catch (MathIllegalArgumentException miae) {
ingo@686:                 if (errors != null) {
ingo@686:                     // TODO: I18N
ingo@686:                     errors.addProblem(km, "spline creation failed");
ingo@686:                 }
felix@1860:                 log.error("spline creation failed", miae);
ingo@686:                 return new double[2][0];
ingo@686:             }
sascha@339: 
sascha@395:             double [] outWs = new double[steps];
sascha@395:             double [] outQs = new double[steps];
sascha@339: 
ingo@686:             Arrays.fill(outWs, Double.NaN);
ingo@686:             Arrays.fill(outQs, Double.NaN);
ingo@686: 
sascha@339:             try {
sascha@395:                 double q = minQ;
sascha@395:                 for (int i = 0; i < outWs.length; ++i, q += stepWidth) {
sascha@395:                     outWs[i] = spline.value(outQs[i] = q);
sascha@339:                 }
sascha@339:             }
sascha@339:             catch (ArgumentOutsideDomainException aode) {
ingo@686:                 if (errors != null) {
ingo@686:                     // TODO: I18N
ingo@686:                     errors.addProblem(km, "spline interpolation failed");
ingo@686:                 }
ingo@686:                 log.error("spline interpolation failed.", aode);
sascha@339:             }
sascha@339: 
sascha@339:             return new double [][] { outWs, outQs };
sascha@339:         }
sascha@339: 
sascha@426: 
sascha@633:         public double getW(QPosition qPosition) {
sascha@633:             int    index  = qPosition.index;
sascha@633:             double weight = qPosition.weight;
sascha@426: 
sascha@633:             return weight == 1.0
sascha@633:                 ? ws[index]
sascha@655:                 : Linear.weight(weight, ws[index-1], ws[index]);
sascha@426:         }
sascha@426: 
sascha@633:         public double getW(
sascha@426:             Row       other,
sascha@426:             double    km,
sascha@633:             QPosition qPosition
sascha@426:         ) {
sascha@655:             double kmWeight = Linear.factor(km, this.km, other.km);
sascha@426: 
sascha@633:             int    index  = qPosition.index;
sascha@633:             double weight = qPosition.weight;
sascha@426: 
sascha@633:             double tw, ow;
sascha@426: 
sascha@633:             if (weight == 1.0) {
sascha@633:                 tw = ws[index];
sascha@633:                 ow = other.ws[index];
sascha@426:             }
sascha@426:             else {
sascha@742:                 tw = Linear.weight(weight, ws[index-1], 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@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() {
sascha@326:         rows = new ArrayList<Row>();
sascha@326:     }
sascha@326: 
sascha@326:     public WstValueTable(Column [] columns) {
sascha@326:         this();
sascha@326:         this.columns = columns;
sascha@326:     }
sascha@326: 
sascha@633:     public WstValueTable(Column [] columns, List<Row> rows) {
sascha@633:         this.columns = columns;
sascha@633:         this.rows    = rows;
sascha@633:     }
sascha@633: 
felix@1838:     /** Sort rows (by km). */
sascha@458:     public void sortRows() {
sascha@458:         Collections.sort(rows);
sascha@458:     }
sascha@458: 
felix@1887:     /**
felix@1887:      * @param km Given kilometer.
felix@1887:      * @param qs Given Q values.
felix@1887:      * @param ws output parameter.
felix@1887:      */
sascha@380:     public double [] interpolateW(double km, double [] qs, 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(
ingo@686:         double      km,
sascha@742:         double []   qs,
ingo@686:         double []   ws,
ingo@686:         Calculation errors
ingo@686:     ) {
sascha@335:         int rowIndex = Collections.binarySearch(rows, new Row(km));
sascha@335: 
sascha@633:         QPosition qPosition = new QPosition();
sascha@633: 
sascha@335:         if (rowIndex >= 0) { // direct row match
sascha@335:             Row row = 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) {
ingo@686:                         // TODO: I18N
ingo@686:                         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))
ingo@686:                     && errors != null) {
ingo@686:                         // TODO: I18N
ingo@686:                         errors.addProblem(
ingo@686:                             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: 
sascha@335:             if (rowIndex < 1 || rowIndex >= rows.size()) {
sascha@335:                 // do not extrapolate
sascha@335:                 Arrays.fill(ws, Double.NaN);
ingo@686:                 if (errors != null) {
ingo@686:                     // TODO: I18N
ingo@686:                     errors.addProblem(km, "km not found");
ingo@686:                 }
sascha@335:             }
sascha@335:             else {
sascha@633:                 Row r1 = rows.get(rowIndex-1);
sascha@633:                 Row r2 = 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: 
felix@1838:     /**
felix@1838:      * Interpolate W and Q values at a given km.
felix@1838:      */
sascha@339:     public double [][] interpolateWQ(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@1838:      */
ingo@686:     public double [][] interpolateWQ(double km, Calculation errors) {
ingo@686:         return interpolateWQ(km, DEFAULT_Q_STEPS, errors);
ingo@686:     }
ingo@686: 
felix@1860:     /**
felix@1860:      * Interpolate W and Q values at a given km.
felix@1860:      */
ingo@686:     public double [][] interpolateWQ(double km, int steps, Calculation errors) {
sascha@395: 
sascha@339:         int rowIndex = Collections.binarySearch(rows, new Row(km));
sascha@339: 
sascha@339:         if (rowIndex >= 0) { // direct row match
sascha@339:             Row row = rows.get(rowIndex);
ingo@686:             return row.interpolateWQ(steps, this, errors);
sascha@339:         }
sascha@339: 
sascha@339:         rowIndex = -rowIndex -1;
sascha@339: 
sascha@339:         if (rowIndex < 1 || rowIndex >= rows.size()) {
sascha@339:             // do not extrapolate
ingo@686:             if (errors != null) {
ingo@686:                 // TODO: I18N
ingo@686:                 errors.addProblem(km, "km not found");
ingo@686:             }
sascha@339:             return new double[2][0];
sascha@339:         }
sascha@339: 
sascha@339:         Row r1 = rows.get(rowIndex-1);
sascha@339:         Row r2 = rows.get(rowIndex);
sascha@339: 
ingo@686:         return r1.interpolateWQ(r2, km, steps, this, errors);
sascha@339:     }
sascha@339: 
sascha@655:     public boolean interpolate(
sascha@655:         double    km,
sascha@655:         double [] out,
sascha@655:         QPosition qPosition,
sascha@655:         Function  qFunction
sascha@649:     ) {
sascha@427:         int R1 = 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: 
sascha@633:         QPosition nPosition = new QPosition();
sascha@655:         if (getQPosition(km, out[1], nPosition) == null) {
sascha@655:             return false;
sascha@655:         }
sascha@427: 
sascha@655:         int rowIndex = Collections.binarySearch(rows, new Row(km));
sascha@427: 
sascha@655:         if (rowIndex >= 0) {
sascha@655:             // direct row match
sascha@655:             out[0] = 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: 
sascha@655:         Row r1 = rows.get(rowIndex-1);
sascha@655:         Row r2 = 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(
ingo@686:         double      q,
ingo@686:         double      referenceKm,
ingo@686:         double []   kms,
ingo@686:         double []   ws,
ingo@686:         double []   qs,
ingo@686:         Calculation errors
sascha@426:     ) {
ingo@686:         return interpolate(
ingo@686:             q, referenceKm, kms, ws, qs, 0, kms.length, errors);
sascha@649:     }
sascha@649: 
sascha@649:     public QPosition interpolate(
ingo@686:         double      q,
ingo@686:         double      referenceKm,
ingo@686:         double []   kms,
ingo@686:         double []   ws,
ingo@686:         double []   qs,
ingo@686:         int         startIndex,
ingo@686:         int         length,
ingo@686:         Calculation errors
sascha@649:     ) {
sascha@645:         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) {
ingo@686:                 errors.addProblem(referenceKm, "cannot find q");
ingo@686:             }
sascha@427:             return null;
sascha@426:         }
sascha@426: 
sascha@645:         Row kmKey = new Row();
sascha@645: 
sascha@645:         int R1 = 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) {
ingo@686:                     errors.addProblem(kms[i], "cannot find q");
ingo@686:                 }
ingo@686:                 ws[i] = Double.NaN;
ingo@686:                 continue;
ingo@686:             }
ingo@686: 
ingo@686:             kmKey.km = kms[i];
sascha@645:             int rowIndex = Collections.binarySearch(rows, kmKey);
sascha@426: 
sascha@426:             if (rowIndex >= 0) {
sascha@426:                 // direct row match
ingo@686:                 if (Double.isNaN(ws[i] = rows.get(rowIndex).getW(qPosition))
ingo@686:                 && errors != null) {
ingo@686:                     errors.addProblem(kms[i], "cannot find w");
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) {
ingo@686:                     errors.addProblem(kms[i], "cannot find km");
ingo@686:                 }
felix@1888:                 ws[i] = Double.NaN;
sascha@426:                 continue;
sascha@426:             }
sascha@426:             Row r1 = rows.get(rowIndex-1);
sascha@426:             Row r2 = rows.get(rowIndex);
sascha@633: 
ingo@686:             if (Double.isNaN(ws[i] = r1.getW(r2, kms[i], qPosition))
ingo@686:             && errors != null) {
ingo@686:                 errors.addProblem(kms[i], "cannot find w");
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:      */
felix@1888:     public static double linearW(double km, Row row1, Row row2, int col) {
felix@1888:         if (row1 == null || row2 == null) {
felix@1888:             return Double.NaN;
felix@1888:         }
felix@1888: 
felix@1888:         return Linear.linear(km,
felix@1888:             row1.km, row2.km,
felix@1888:             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:      */
felix@1888:     public double [][] interpolateWQColumnwise(double km) {
felix@1888:         double [] qs = new double[columns.length];
felix@1888:         double [] ws = new double[columns.length];
felix@1888: 
felix@1888:         // Find out row from where we will start searching.
felix@1888:         int rowIndex = Collections.binarySearch(rows, new Row(km));
felix@1888: 
felix@1888:         if (rowIndex < 0) {
felix@1888:             rowIndex = -rowIndex -1;
felix@1888:         }
felix@1888: 
felix@1888:         Row startRow = rows.get(rowIndex);
felix@1888: 
felix@1888:         for (int col = 0; col < columns.length; col++) {
felix@1888:             qs[col] = columns[col].getQRangeTree().findQ(km);
felix@1888:             if (startRow.km == km && startRow.ws[col] != Double.NaN) {
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--) {
felix@1888:                 if (!Double.isNaN(rows.get(before).ws[col])) {
felix@1888:                     rowBefore = rows.get(before);
felix@1888:                     break;
felix@1888:                 }
felix@1888:             }
felix@1888:             for (int after = rowIndex; after < rows.size(); after++) {
felix@1888:                 if (!Double.isNaN(rows.get(after).ws[col])) {
felix@1888:                     rowAfter = rows.get(after);
felix@1888:                     break;
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: 
sascha@633:     public QPosition getQPosition(double km, double q) {
sascha@633:         return getQPosition(km, q, new QPosition());
sascha@633:     }
sascha@633: 
sascha@633:     public QPosition getQPosition(double km, double q, QPosition qPosition) {
sascha@633: 
sascha@633:         if (columns.length == 0) {
sascha@633:             return null;
sascha@633:         }
sascha@633: 
sascha@633:         double qLast = 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: 
sascha@633:         for (int i = 1; i < columns.length; ++i) {
sascha@633:             double qCurrent = 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) {
sascha@655:                 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: 
sascha@633:     public double getQIndex(int index, double km) {
sascha@633:         return columns[index].getQRangeTree().findQ(km);
sascha@633:     }
sascha@633: 
sascha@633:     public double getQ(QPosition qPosition, double km) {
sascha@633:         int    index  = qPosition.index;
sascha@633:         double weight = qPosition.weight;
sascha@633: 
sascha@633:         if (weight == 1d) {
sascha@633:             return columns[index].getQRangeTree().findQ(km);
sascha@633:         }
sascha@633:         double q1 = columns[index-1].getQRangeTree().findQ(km);
sascha@633:         double q2 = columns[index  ].getQRangeTree().findQ(km);
sascha@655:         return Linear.weight(weight, q1, q2);
sascha@633:     }
sascha@326: }
sascha@326: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :