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: 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, sascha@633: WstValueTable table 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); sascha@655: ows[i] = Linear.weight(kmWeight, wt, wo); sascha@633: } sascha@633: else { 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, sascha@633: WstValueTable table sascha@633: ) { sascha@395: int W = Math.min(ws.length, other.ws.length); sascha@339: sascha@339: if (W < 1) { 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: sascha@395: splineW[i] = wws; sascha@395: splineQ[i] = wqs; sascha@395: if (wqs < minQ) minQ = wqs; sascha@395: if (wqs > maxQ) maxQ = wqs; sascha@339: } sascha@339: sascha@395: double stepWidth = (maxQ - minQ)/steps; sascha@395: sascha@339: SplineInterpolator interpolator = new SplineInterpolator(); sascha@395: PolynomialSplineFunction spline = sascha@395: interpolator.interpolate(splineQ, splineW); sascha@339: sascha@395: double [] outWs = new double[steps]; sascha@395: double [] outQs = new double[steps]; sascha@339: 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) { sascha@339: log.error("Spline interpolation failed.", aode); sascha@339: } sascha@339: sascha@339: return new double [][] { outWs, outQs }; sascha@339: } sascha@339: sascha@633: public double [][] interpolateWQ(int steps, WstValueTable table) { sascha@395: sascha@395: int W = ws.length; sascha@339: sascha@339: if (W < 1) { sascha@339: return new double[2][0]; sascha@339: } sascha@339: sascha@395: double [] splineQ = new double[W]; sascha@395: sascha@395: double minQ = Double.MAX_VALUE; sascha@395: double maxQ = -Double.MAX_VALUE; sascha@339: sascha@633: QPosition qPosition = new QPosition(); sascha@633: sascha@339: for (int i = 0; i < W; ++i) { sascha@633: splineQ[i] = table.getQIndex(i, km); sascha@633: if (splineQ[i] < minQ) minQ = splineQ[i]; sascha@633: if (splineQ[i] > maxQ) maxQ = splineQ[i]; sascha@339: } sascha@339: sascha@395: double stepWidth = (maxQ - minQ)/steps; sascha@395: sascha@339: SplineInterpolator interpolator = new SplineInterpolator(); sascha@339: sascha@395: PolynomialSplineFunction spline = sascha@633: interpolator.interpolate(splineQ, ws); sascha@339: sascha@395: double [] outWs = new double[steps]; sascha@395: double [] outQs = new double[steps]; sascha@339: 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) { sascha@339: 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@655: tw = Linear.weight(weight, ws[index-1], ws[index]); sascha@655: 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) { sascha@335: 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) { sascha@633: ws[i] = getQPosition(km, qs[i], qPosition) != null sascha@633: ? row.getW(qPosition) sascha@633: : Double.NaN; 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); sascha@335: } sascha@335: else { sascha@633: Row r1 = rows.get(rowIndex-1); sascha@633: Row r2 = rows.get(rowIndex); sascha@633: r1.interpolateW(r2, km, qs, ws, this); sascha@335: } sascha@335: } sascha@335: sascha@335: return ws; sascha@335: } sascha@335: sascha@339: public double [][] interpolateWQ(double km) { sascha@395: return interpolateWQ(km, DEFAULT_Q_STEPS); sascha@339: } sascha@339: sascha@395: public double [][] interpolateWQ(double km, int steps) { 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); sascha@633: return row.interpolateWQ(steps, this); sascha@339: } sascha@339: sascha@339: rowIndex = -rowIndex -1; sascha@339: sascha@339: if (rowIndex < 1 || rowIndex >= rows.size()) { sascha@339: // do not extrapolate 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: sascha@633: return r1.interpolateWQ(r2, km, steps, this); 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@655: 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: sascha@427: public QPosition interpolate( sascha@426: double q, sascha@645: double referenceKm, sascha@426: double [] kms, sascha@426: double [] ws, sascha@426: double [] qs sascha@426: ) { sascha@649: return interpolate(q, referenceKm, kms, ws, qs, 0, kms.length); sascha@649: } sascha@649: sascha@649: public QPosition interpolate( sascha@649: double q, sascha@649: double referenceKm, sascha@649: double [] kms, sascha@649: double [] ws, sascha@649: double [] qs, sascha@649: int startIndex, sascha@649: int length sascha@649: ) { sascha@645: QPosition qPosition = getQPosition(referenceKm, q); sascha@426: sascha@426: if (qPosition == null) { sascha@633: // we cannot locate q at km 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@426: kmKey.km = kms[i]; sascha@458: sascha@645: int rowIndex = Collections.binarySearch(rows, kmKey); sascha@426: sascha@633: qs[i] = getQ(qPosition, kms[i]); sascha@633: sascha@426: if (rowIndex >= 0) { sascha@426: // direct row match sascha@633: ws[i] = rows.get(rowIndex).getW(qPosition); sascha@426: continue; sascha@426: } sascha@426: sascha@426: rowIndex = -rowIndex -1; sascha@426: sascha@426: if (rowIndex < 1 || rowIndex >= R1) { sascha@426: // do not extrapolate 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: sascha@633: ws[i] = r1.getW(r2, kms[i], qPosition); 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 :