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: 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: 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@1907: 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@1939: 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@1939: int numSamples, sascha@1939: 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@1939: // TODO: I18N sascha@1939: 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: /** 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@1938: sascha@1938: public SplineFunction createSpline( sascha@1938: WstValueTable table, sascha@1938: Calculation errors sascha@1938: ) { sascha@1938: int W = ws.length; sascha@1938: sascha@1938: if (W < 1) { sascha@1938: if (errors != null) { sascha@1938: // TODO: I18N sascha@1938: errors.addProblem(km, "no ws found"); sascha@1938: } sascha@1938: return null; sascha@1938: } sascha@1938: sascha@1938: double [] splineQs = new double[W]; sascha@1938: sascha@1938: for (int i = 0; i < W; ++i) { sascha@1938: double sq = table.getQIndex(i, km); sascha@1938: if (Double.isNaN(sq) && errors != null) { sascha@1938: // TODO: I18N sascha@1938: errors.addProblem( sascha@1938: km, "no q found in " + (i+1) + " column"); sascha@1938: } sascha@1938: splineQs[i] = sq; sascha@1938: } sascha@1938: sascha@1938: try { sascha@1938: SplineInterpolator interpolator = new SplineInterpolator(); sascha@1938: PolynomialSplineFunction spline = sascha@1938: interpolator.interpolate(splineQs, ws); sascha@1938: sascha@1938: return new SplineFunction(spline, splineQs, ws); sascha@1938: } sascha@1938: catch (MathIllegalArgumentException miae) { sascha@1938: if (errors != null) { sascha@1938: // TODO: I18N sascha@1938: 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: ) { 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@1938: return null; sascha@339: } sascha@339: sascha@655: double factor = Linear.factor(km, this.km, other.km); sascha@339: sascha@1938: double [] splineQs = new double[W]; sascha@1938: double [] splineWs = new double[W]; 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: sascha@1938: splineWs[i] = wws; sascha@1938: splineQs[i] = wqs; sascha@339: } sascha@339: sascha@339: SplineInterpolator interpolator = new SplineInterpolator(); ingo@686: ingo@686: try { sascha@1938: PolynomialSplineFunction spline = sascha@1938: interpolator.interpolate(splineQs, splineWs); sascha@1938: sascha@1938: return new SplineFunction(spline, splineQs, splineWs); 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: } 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: sascha@1937: public double [] findQsForW(double w, WstValueTable table) { sascha@1937: sascha@1937: TDoubleArrayList qs = new TDoubleArrayList(); sascha@1937: sascha@1937: if (ws.length > 0 && Math.abs(ws[0]-w) < 0.000001) { sascha@1937: double q = table.getQIndex(0, km); sascha@1937: if (!Double.isNaN(q)) { sascha@1937: qs.add(q); sascha@1937: } sascha@1937: } sascha@1937: sascha@1937: for (int i = 1; i < ws.length; ++i) { sascha@1937: double w2 = ws[i]; sascha@1937: if (Double.isNaN(w2)) { sascha@1937: continue; sascha@1937: } sascha@1937: if (Math.abs(w2-w) < 0.000001) { sascha@1937: double q = table.getQIndex(i, km); sascha@1937: if (!Double.isNaN(q)) { sascha@1937: qs.add(q); sascha@1937: } sascha@1937: continue; sascha@1937: } sascha@1937: double w1 = ws[i-1]; sascha@1937: if (Double.isNaN(w1)) { sascha@1937: continue; sascha@1937: } sascha@1937: sascha@1937: if (w < Math.min(w1, w2) || w > Math.max(w1, w2)) { sascha@1937: continue; sascha@1937: } sascha@1937: sascha@1937: double q1 = table.getQIndex(i-1, km); sascha@1937: double q2 = table.getQIndex(i, km); sascha@1937: if (Double.isNaN(q1) || Double.isNaN(q2)) { sascha@1937: continue; sascha@1937: } sascha@1937: sascha@1937: double q = Linear.linear(w, w1, w2, q1, q2); sascha@1937: qs.add(q); sascha@1937: } sascha@1937: sascha@1937: return qs.toNativeArray(); sascha@1937: } sascha@1937: sascha@1937: public double [] findQsForW( sascha@1937: Row other, sascha@1937: double w, sascha@1937: double km, sascha@1937: WstValueTable table sascha@1937: ) { sascha@1937: TDoubleArrayList qs = new TDoubleArrayList(); sascha@1937: sascha@1937: double factor = Linear.factor(km, this.km, other.km); sascha@1937: sascha@1937: if (ws.length > 0) { sascha@1937: double wt = Linear.weight(factor, ws[0], other.ws[0]); sascha@1937: if (!Double.isNaN(wt)) { sascha@1937: double q = table.getQIndex(0, km); sascha@1937: if (!Double.isNaN(q)) { sascha@1937: qs.add(q); sascha@1937: } sascha@1937: } sascha@1937: } sascha@1937: sascha@1937: for (int i = 1; i < ws.length; ++i) { sascha@1937: double w2 = Linear.weight(factor, ws[i], other.ws[i]); sascha@1937: if (Double.isNaN(w2)) { sascha@1937: continue; sascha@1937: } sascha@1937: if (Math.abs(w2-w) < 0.000001) { sascha@1937: double q = table.getQIndex(i, km); sascha@1937: if (!Double.isNaN(q)) { sascha@1937: qs.add(q); sascha@1937: } sascha@1937: continue; sascha@1937: } sascha@1937: double w1 = Linear.weight(factor, ws[i-1], other.ws[i-1]); sascha@1937: if (Double.isNaN(w1)) { sascha@1937: continue; sascha@1937: } sascha@1937: sascha@1937: if (w < Math.min(w1, w2) || w > Math.max(w1, w2)) { sascha@1937: continue; sascha@1937: } sascha@1937: sascha@1937: double q1 = table.getQIndex(i-1, km); sascha@1937: double q2 = table.getQIndex(i, km); sascha@1937: if (Double.isNaN(q1) || Double.isNaN(q2)) { sascha@1937: continue; sascha@1937: } sascha@1937: sascha@1937: double q = Linear.linear(w, w1, w2, q1, q2); sascha@1937: qs.add(q); sascha@1937: } sascha@1937: sascha@1937: return qs.toNativeArray(); sascha@1937: } 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: sascha@633: public WstValueTable(Column [] columns, List 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@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: } 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); 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: } sascha@1907: for (int after = rowIndex, R = rows.size(); after < R; 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@1937: public double [] findQsForW(double km, double w) { sascha@1937: sascha@1937: int rowIndex = Collections.binarySearch(rows, new Row(km)); sascha@1937: sascha@1937: if (rowIndex >= 0) { sascha@1937: return rows.get(rowIndex).findQsForW(w, this); 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: sascha@1937: return r1.findQsForW(r2, w, km, this); 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@1939: // TODO: I18N sascha@1939: errors.addProblem(km, "cannot create w/q 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@1939: // TODO: I18N sascha@1939: 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@1939: // TODO: I18N sascha@1939: errors.addProblem(km, "cannot create w/q relation"); sascha@1939: } sascha@1939: sascha@1939: return sf; sascha@1939: } sascha@1939: sascha@1939: /** 'Bezugslinienverfahren' */ sascha@1939: public double [][] relateWs( sascha@1939: double km1, sascha@1939: 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: sascha@1939: SplineFunction sf2 = createSpline(km1, 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@1939: // TODO: I18N sascha@1939: errors.addProblem(km1, "cannot create index/q relation"); sascha@1939: } sascha@1939: return new double[2][0]; sascha@1939: } sascha@1939: sascha@1939: PolynomialSplineFunction iQ2 = sf1.createIndexQRelation(); sascha@1939: if (iQ2 == null) { sascha@1939: if (errors != null) { sascha@1939: // TODO: I18N sascha@1939: errors.addProblem(km2, "cannot create index/q relation"); sascha@1939: } sascha@1939: return new double[2][0]; sascha@1939: } sascha@1939: sascha@1939: int N = sf1.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@1939: double [] ws1 = new double[numSamples]; sascha@1939: double [] ws2 = new double[numSamples]; sascha@1939: sascha@1939: Arrays.fill(ws1, Double.NaN); sascha@1939: Arrays.fill(ws2, Double.NaN); sascha@1939: sascha@1939: boolean hadErrors = false; sascha@1939: sascha@1939: double p = 0d; sascha@1939: for (int i = 0; i < numSamples; ++i, p += stepWidth) { sascha@1939: try { sascha@1939: double q1 = iQ1.value(p); sascha@1939: double w1 = qW1.value(q1); sascha@1939: double q2 = iQ2.value(p); sascha@1939: double w2 = qW2.value(q2); sascha@1939: ws1[i] = w1; sascha@1939: ws2[i] = w2; sascha@1939: } sascha@1939: catch (ArgumentOutsideDomainException aode) { sascha@1939: if (!hadErrors) { sascha@1939: // XXX: I'm not sure if this really can happen sascha@1939: // and if we should report this more than once. sascha@1939: hadErrors = true; sascha@1939: if (errors != null) { sascha@1939: // TODO: I18N sascha@1939: log.debug("W~W failed", aode); sascha@1939: errors.addProblem("W~W failed"); sascha@1939: } sascha@1939: } sascha@1939: } sascha@1939: } sascha@1939: sascha@1939: return new double [][] { ws1, ws2 }; 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@326: } sascha@326: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :