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@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@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: 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@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: 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: 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: } ingo@686: log.error("spline creation failed"); 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: 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: } ingo@686: log.error("spline creation failed"); 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: sascha@326: protected List rows; sascha@326: 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: sascha@458: public void sortRows() { sascha@458: Collections.sort(rows); sascha@458: } sascha@458: 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: sascha@339: public double [][] interpolateWQ(double km) { ingo@686: return interpolateWQ(km, null); sascha@339: } sascha@339: ingo@686: public double [][] interpolateWQ(double km, Calculation errors) { ingo@686: return interpolateWQ(km, DEFAULT_Q_STEPS, errors); ingo@686: } ingo@686: 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: } sascha@633: 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: 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@633: sascha@326: } sascha@326: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :