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: sascha@326: import java.io.Serializable; sascha@326: teichmann@5831: import org.dive4elements.river.artifacts.math.Linear; teichmann@5831: import org.dive4elements.river.artifacts.math.Function; teichmann@6434: import org.dive4elements.river.utils.DoubleUtil; 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: sascha@1937: import gnu.trove.TDoubleArrayList; sascha@1937: tom@8662: import static org.dive4elements.river.backend.utils.EpsilonComparator.CMP; tom@8662: 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: 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: teichmann@4797: 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: 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( sascha@1939: PolynomialSplineFunction spline, sascha@3076: double [] splineQs, sascha@1939: double [] splineWs sascha@1939: ) { sascha@1939: this.spline = spline; sascha@1939: this.splineQs = splineQs; sascha@1939: this.splineWs = splineWs; sascha@1939: } sascha@1939: sascha@1939: public double [][] sample( sascha@3076: int numSamples, sascha@3076: double km, sascha@1939: Calculation errors sascha@1939: ) { sascha@1939: double minQ = getQMin(); sascha@1939: double maxQ = getQMax(); sascha@1939: sascha@1939: double [] outWs = new double[numSamples]; sascha@1939: double [] outQs = new double[numSamples]; sascha@1939: sascha@1939: Arrays.fill(outWs, Double.NaN); sascha@1939: Arrays.fill(outQs, Double.NaN); sascha@1939: sascha@1939: 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) { sascha@1939: outWs[i] = spline.value(outQs[i] = q); sascha@1939: } sascha@1939: } sascha@1939: catch (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() { sascha@1939: return Math.min(splineQs[0], splineQs[splineQs.length-1]); sascha@1939: } sascha@1939: sascha@1939: public double getQMax() { sascha@1939: return Math.max(splineQs[0], splineQs[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: sascha@1939: double [] indices = new double[splineQs.length]; sascha@1939: for (int i = 0; i < indices.length; ++i) { sascha@1939: indices[i] = i; sascha@1939: } sascha@1939: sascha@1939: try { sascha@1939: SplineInterpolator interpolator = new SplineInterpolator(); sascha@1939: return interpolator.interpolate(indices, splineQs); sascha@1939: } sascha@1939: catch (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 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: /** tom@8704: * Sort Qs and Ws for this Row over Q. tom@8704: */ tom@8704: private double[][] getSortedWQ( tom@8704: WstValueTable table, tom@8704: Calculation errors tom@8704: ) { tom@8704: int W = ws.length; tom@8704: tom@8704: if (W < 1) { tom@8704: if (errors != null) { tom@8704: errors.addProblem(km, "no.ws.found"); tom@8704: } tom@8704: return new double[][] {{Double.NaN}, {Double.NaN}}; tom@8704: } tom@8704: tom@8704: double [] sortedWs = ws; tom@8704: double [] sortedQs = new double[W]; tom@8704: tom@8704: for (int i=0; i < W; ++i) { tom@8704: double q = table.getQIndex(i, km); tom@8704: if (Double.isNaN(q) && errors != null) { tom@8704: errors.addProblem( tom@8704: 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( tom@8704: Row other, tom@8704: double km, tom@8704: WstValueTable table, tom@8704: Calculation errors tom@8704: ) { tom@8704: int W = Math.min(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: tom@8704: double factor = Linear.factor(km, this.km, other.km); tom@8704: tom@8704: double [] sortedQs = new double[W]; tom@8704: double [] sortedWs = new double[W]; tom@8704: tom@8704: for (int i = 0; i < W; ++i) { tom@8704: double wws = Linear.weight(factor, ws[i], other.ws[i]); tom@8704: 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: */ tom@8704: private double[] findQsForW(double w, double[][] sortedWQ) { tom@8704: int W = sortedWQ[0].length; tom@8704: tom@8704: double[] sortedQs = sortedWQ[1]; tom@8704: double[] sortedWs = sortedWQ[0]; tom@8704: tom@8704: TDoubleArrayList qs = new TDoubleArrayList(); tom@8704: tom@8704: if (W > 0 && Math.abs(sortedWs[0]-w) < W_EPSILON) { tom@8704: 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) { tom@8704: double w2 = sortedWs[i]; tom@8704: if (Double.isNaN(w2)) { tom@8704: continue; tom@8704: } tom@8704: if (Math.abs(w2-w) < W_EPSILON) { tom@8704: double q = sortedQs[i]; tom@8704: if (!Double.isNaN(q)) { tom@8704: qs.add(q); tom@8704: } tom@8704: continue; tom@8704: } tom@8704: 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: tom@8704: double q1 = sortedQs[i-1]; tom@8704: double q2 = sortedQs[i]; tom@8704: if (Double.isNaN(q1) || Double.isNaN(q2)) { tom@8704: continue; tom@8704: } tom@8704: tom@8704: 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: */ sascha@335: public int compareTo(Row other) { tom@8662: return CMP.compare(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( 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: errors.addProblem( sascha@2166: 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( sascha@1938: WstValueTable table, sascha@1938: Calculation errors sascha@1938: ) { tom@8704: double[][] sortedWQ = getSortedWQ(table, errors); teichmann@7000: sascha@1938: try { sascha@1938: SplineInterpolator interpolator = new SplineInterpolator(); sascha@1938: PolynomialSplineFunction spline = tom@8704: interpolator.interpolate(sortedWQ[1], sortedWQ[0]); sascha@1938: tom@8704: return new SplineFunction(spline, sortedWQ[1], ws); sascha@1938: } sascha@1938: catch (MathIllegalArgumentException miae) { sascha@1938: if (errors != null) { sascha@2166: errors.addProblem(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( sascha@633: Row other, sascha@633: double km, ingo@686: WstValueTable table, ingo@686: Calculation errors sascha@633: ) { tom@8704: double[][] sortedWQ = getSortedWQ(other, km, table, errors); teichmann@6434: sascha@339: SplineInterpolator interpolator = new SplineInterpolator(); ingo@686: ingo@686: try { sascha@1938: PolynomialSplineFunction spline = tom@8704: interpolator.interpolate(sortedWQ[1], sortedWQ[0]); sascha@1938: tom@8704: return new SplineFunction(spline, sortedWQ[1], sortedWQ[0]); ingo@686: } ingo@686: catch (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( sascha@1938: Row other, sascha@1938: double km, sascha@1938: int steps, sascha@1938: WstValueTable table, sascha@1938: Calculation errors sascha@1938: ) { sascha@1938: SplineFunction sf = createSpline(other, km, table, errors); sascha@339: sascha@1938: return sf != null sascha@1938: ? sf.sample(steps, km, errors) sascha@1938: : new double[2][0]; 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@1938: SplineFunction sf = createSpline(table, errors); sascha@339: sascha@1938: return sf != null sascha@1938: ? sf.sample(steps, km, errors) sascha@1938: : new double[2][0]; 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@1937: tom@8704: public double [] findQsForW( tom@8704: double w, tom@8704: WstValueTable table, tom@8704: Calculation errors tom@8704: ) { tom@8704: log.debug("Find Qs for given W at tabulated km " + km); tom@8704: return findQsForW(w, getSortedWQ(table, errors)); sascha@1937: } sascha@1937: sascha@1937: public double [] findQsForW( sascha@1937: Row other, sascha@1937: double w, sascha@3076: double km, tom@8704: WstValueTable table, tom@8704: Calculation errors sascha@1937: ) { 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: ingo@2423: public double [] getMinMaxW(double [] result) { ingo@2423: double minW = Double.MAX_VALUE; ingo@2423: double maxW = -Double.MAX_VALUE; ingo@2423: for (int i = 0; i < ws.length; ++i) { ingo@2423: double w = 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: ingo@2423: public double [] getMinMaxW(Row other, double km, double [] result) { ingo@2423: double [] m1 = this .getMinMaxW(new double [2]); ingo@2423: double [] m2 = other.getMinMaxW(new double [2]); ingo@2423: 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 rows; sascha@326: felix@1838: /** Columns in table. */ sascha@326: protected Column [] columns; sascha@326: sascha@326: public WstValueTable() { sascha@326: rows = new ArrayList(); sascha@326: } sascha@326: sascha@326: public WstValueTable(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: */ sascha@633: public WstValueTable(Column [] columns, List rows) { sascha@633: this.columns = columns; sascha@633: this.rows = rows; sascha@633: } sascha@633: teichmann@4797: public Column [] getColumns() { teichmann@4797: return 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: */ 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) { 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)) ingo@686: && errors != null) { ingo@686: errors.addProblem( sascha@2166: 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) { sascha@2166: 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: sascha@2559: public double [] getMinMaxQ(double km) { sascha@2559: return getMinMaxQ(km, new double [2]); sascha@2559: } sascha@2559: sascha@2559: public double [] getMinMaxQ(double km, double [] result) { sascha@2559: double minQ = Double.MAX_VALUE; sascha@2559: double maxQ = -Double.MAX_VALUE; sascha@2559: sascha@2559: for (int i = 0; i < columns.length; ++i) { sascha@2559: double q = 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) { sascha@2560: 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) { sascha@2560: 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 sascha@2560: ? new double [] { minQ, maxQ } sascha@2560: : null; sascha@2560: } sascha@2559: ingo@2423: public double [] getMinMaxW(double km) { ingo@2423: return getMinMaxW(km, new double [2]); ingo@2423: ingo@2423: } ingo@2423: public double [] getMinMaxW(double km, double [] result) { ingo@2423: int rowIndex = Collections.binarySearch(rows, new Row(km)); ingo@2423: ingo@2423: if (rowIndex >= 0) { ingo@2423: return rows.get(rowIndex).getMinMaxW(result); ingo@2423: } ingo@2423: ingo@2423: rowIndex = -rowIndex -1; ingo@2423: ingo@2423: if (rowIndex < 1 || rowIndex >= rows.size()) { ingo@2423: // do not extrapolate ingo@2423: return null; ingo@2423: } ingo@2423: ingo@2423: Row r1 = rows.get(rowIndex-1); ingo@2423: Row r2 = 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) { ingo@2423: 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) { ingo@2423: 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 ingo@2423: ? new double [] { minW, maxW } ingo@2423: : null; ingo@2423: } ingo@2423: 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@3265: * felix@3265: * @param errors where to store errors. felix@3265: * felix@3265: * @return double double array, first index Ws, second Qs. felix@1838: */ ingo@686: public double [][] interpolateWQ(double km, 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( tom@8856: double km, tom@8856: int steps, tom@8856: Calculation errors tom@8856: ) { 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) { sascha@2166: 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: 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( 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) { sascha@2166: errors.addProblem(referenceKm, "cannot.find.q", 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) { 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]; 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) { 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: } 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) { 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: */ 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@1919: log.debug("WstValueTable.interpolateWQColumnwise"); 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@2248: felix@2248: // TODO Beyond definition, we could stop more clever. sascha@1937: if (rowIndex >= rows.size()) { felix@1919: rowIndex = rows.size() -1; sascha@1937: } 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); 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--) { felix@1888: if (!Double.isNaN(rows.get(before).ws[col])) { felix@1888: rowBefore = rows.get(before); felix@1888: break; felix@1888: } felix@1888: } felix@2248: if (rowBefore != null) { tom@8856: for (int after = rowIndex, R = rows.size(); tom@8856: after < R; tom@8856: after++ tom@8856: ) { felix@2248: if (!Double.isNaN(rows.get(after).ws[col])) { felix@2248: rowAfter = 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: tom@8704: public double [] findQsForW(double km, double w, Calculation errors) { sascha@1937: sascha@1937: int rowIndex = Collections.binarySearch(rows, new Row(km)); sascha@1937: sascha@1937: if (rowIndex >= 0) { tom@8704: return rows.get(rowIndex).findQsForW(w, this, errors); sascha@1937: } sascha@1937: sascha@1937: rowIndex = -rowIndex - 1; sascha@1937: sascha@1937: if (rowIndex < 1 || rowIndex >= rows.size()) { sascha@1939: // Do not extrapolate. sascha@1937: return new double[0]; sascha@1937: } sascha@1937: sascha@1939: // Needs bilinear interpolation. sascha@1937: Row r1 = rows.get(rowIndex-1); sascha@1937: Row r2 = rows.get(rowIndex); sascha@1937: tom@8704: return r1.findQsForW(r2, w, km, this, errors); sascha@1937: } sascha@1937: sascha@1939: protected SplineFunction createSpline(double km, Calculation errors) { sascha@1939: sascha@1939: int rowIndex = Collections.binarySearch(rows, new Row(km)); sascha@1939: sascha@1939: if (rowIndex >= 0) { sascha@1939: SplineFunction sf = 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: sascha@1939: if (rowIndex < 1 || rowIndex >= 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. sascha@1939: Row r1 = rows.get(rowIndex-1); sascha@1939: Row r2 = rows.get(rowIndex); sascha@1939: sascha@1939: 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( sascha@3076: double km1, sascha@1939: double km2, sascha@2186: Calculation errors sascha@2186: ) { 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: ingo@2330: ErrorHandler(Calculation errors) { ingo@2330: this.errors = errors; ingo@2330: } ingo@2330: ingo@2330: void error(double km, String key, Object ... args) { ingo@2330: if (errors != null && !hasErrors) { ingo@2330: hasErrors = true; ingo@2330: 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( sascha@3076: double km1, sascha@2186: double km2, sascha@1939: int numSamples, sascha@1939: Calculation errors sascha@1939: ) { sascha@1939: SplineFunction sf1 = createSpline(km1, errors); sascha@1939: if (sf1 == null) { sascha@1939: return new double[2][0]; sascha@1939: } sascha@1939: felix@2248: SplineFunction sf2 = createSpline(km2, errors); sascha@1939: if (sf2 == null) { sascha@1939: return new double[2][0]; sascha@1939: } sascha@1939: sascha@1939: 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: sascha@2404: 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: sascha@2316: int N = Math.min(sf1.splineQs.length, sf2.splineQs.length); sascha@1939: double stepWidth = N/(double)numSamples; sascha@1939: sascha@1939: PolynomialSplineFunction qW1 = sf1.spline; sascha@1939: PolynomialSplineFunction qW2 = sf2.spline; sascha@1939: sascha@2316: TDoubleArrayList ws1 = new TDoubleArrayList(numSamples); sascha@2316: TDoubleArrayList ws2 = new TDoubleArrayList(numSamples); sascha@2316: TDoubleArrayList qs1 = new TDoubleArrayList(numSamples); sascha@2316: TDoubleArrayList qs2 = new TDoubleArrayList(numSamples); sascha@1939: ingo@2330: 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: } sascha@1939: catch (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: } ingo@2330: catch (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: } ingo@2330: catch (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: } ingo@2330: catch (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: 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@3736: sascha@3736: public double [][] interpolateTabulated(double km) { sascha@3736: return interpolateTabulated(km, new double[2][columns.length]); sascha@3736: } sascha@3736: sascha@3736: public double [][] interpolateTabulated(double km, double [][] result) { sascha@3736: sascha@3736: int rowIndex = Collections.binarySearch(rows, new Row(km)); sascha@3736: sascha@3736: if (rowIndex >= 0) { sascha@3736: // Direct hit -> copy ws. sascha@3736: Row row = rows.get(rowIndex); sascha@3736: System.arraycopy( sascha@3736: row.ws, 0, result[0], 0, sascha@3736: Math.min(row.ws.length, result[0].length)); sascha@3736: } sascha@3736: else { sascha@3736: rowIndex = -rowIndex -1; sascha@3736: if (rowIndex < 1 || rowIndex >= rows.size()) { sascha@3736: // Out of bounds. sascha@3736: return null; sascha@3736: } sascha@3736: // Interpolate ws. sascha@3736: Row r1 = rows.get(rowIndex-1); sascha@3736: Row r2 = rows.get(rowIndex); sascha@3736: double factor = Linear.factor(km, r1.km, r2.km); sascha@3736: Linear.weight(factor, r1.ws, r2.ws, result[0]); sascha@3736: } sascha@3736: sascha@3736: double [] qs = result[1]; sascha@3736: for (int i = Math.min(qs.length, columns.length)-1; i >= 0; --i) { sascha@3736: qs[i] = 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() { felix@6923: for (Column column: 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: felix@6923: if (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?) */ sascha@3743: public List findSegments(double km1, double km2) { sascha@3743: return columns.length != 0 sascha@3743: ? columns[columns.length-1].getQRangeTree().findSegments(km1, km2) sascha@3743: : Collections.emptyList(); sascha@3743: } sascha@326: } sascha@326: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :