Mercurial > dive4elements > river
diff flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java @ 633:d08f77e7f7e8
WST value table: Qs are now stored in ranges for each column.
flys-artifacts/trunk@2006 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Wed, 25 May 2011 15:31:25 +0000 |
parents | 07640ab913fd |
children | 433f67a076aa |
line wrap: on
line diff
--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java Tue May 24 14:46:45 2011 +0000 +++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java Wed May 25 15:31:25 2011 +0000 @@ -25,8 +25,8 @@ public static final int DEFAULT_Q_STEPS = 500; - public static class Column - implements Serializable + public static final class Column + implements Serializable { protected String name; @@ -54,32 +54,34 @@ public void setQRangeTree(QRangeTree qRangeTree) { this.qRangeTree = qRangeTree; } - } - // class Column + } // class Column - public static class QPosition { + public static final class QPosition { - protected double q; protected int index; protected double weight; public QPosition() { } - public QPosition(double q, int index, double weight) { + public QPosition(int index, double weight) { this.index = index; this.weight = weight; } + public QPosition set(int index, double weight) { + this.index = index; + this.weight = weight; + return this; + } + } // class Position - public static class Row - implements Serializable, Comparable<Row> + public static final class Row + implements Serializable, Comparable<Row> { double km; double [] ws; - double [] qs; - boolean qSorted; public Row() { } @@ -88,110 +90,9 @@ this.km = km; } - public Row(double km, double [] ws, double [] qs) { + public Row(double km, double [] ws) { this(km); this.ws = ws; - this.qs = qs; - } - - public static final double linear( - double x, - double x1, double x2, - double y1, double y2 - ) { - // y1 = m*x1 + b - // y2 = m*x2 + b - // y2 - y1 = m*(x2 - x1) - // m = (y2 - y1)/(x2 - x1) # x2 != x1 - // b = y1 - m*x1 - - if (x1 == x2) { - return 0.5*(y1 + y2); - } - double m = (y2 - y1)/(x2 - x1); - double b = y1 - m*x1; - return x*m + b; - } - - public static final double factor(double x, double p1, double p2) { - // 0 = m*p1 + b <=> b = -m*p1 - // 1 = m*p2 + b - // 1 = m*(p2 - p1) - // m = 1/(p2 - p1) # p1 != p2 - // f(x) = x/(p2-p1) - p1/(p2-p1) <=> (x-p1)/(p2-p1) - - return p1 == p2 ? 0.0 : (x-p1)/(p2-p1); - } - - public static final double weight(double factor, double a, double b) { - return (1.0-factor)*a + factor*b; - } - - public void interpolateW( - Row other, - double km, - double [] input, - double [] output - ) { - double factor = factor(km, this.km, other.km); - - for (int i = 0; i < input.length; ++i) { - double q = input[i]; - int idx1 = getQIndex(q); - int idx2 = other.getQIndex(q); - - double w1 = idx1 >= 0 - ? ws[idx1] - : interpolateW(-idx1-1, q); - - double w2 = idx2 >= 0 - ? other.ws[idx2] - : other.interpolateW(-idx2-1, q); - - output[i] = weight(factor, w1, w2); - } - } - - public double interpolateW(int idx, double q) { - return idx < 1 || idx >= qs.length - ? Double.NaN // do not extrapolate - : linear(q, qs[idx-1], qs[idx], ws[idx-1], ws[idx]); - } - - public double interpolateW(double q) { - if (Double.isNaN(q)) return Double.NaN; - int index = getQIndex(q); - return index >= 0 ? ws[index] : interpolateW(-index -1, q); - } - - public int getQIndex(double q) { - return qSorted ? binaryQIndex(q) : linearQIndex(q); - } - - public int binaryQIndex(double q) { - return Arrays.binarySearch(qs, q); - } - - public int linearQIndex(double q) { - switch (qs.length) { - case 0: break; - case 1: - if (qs[0] == q) return 0; - if (qs[0] < q) return -(1+1); - return -(0+1); - default: - for (int i = 1; i < qs.length; ++i) { - double qa = qs[i-1]; - double qb = qs[i]; - if (qa == q) return i-1; - if (qb == q) return i; - if (qa > qb) { double t = qa; qa = qb; qb = t; } - if (q > qa && q < qb) return -(i+1); - } - return -qs.length; - } - - return -1; } public int compareTo(Row other) { @@ -201,8 +102,35 @@ return 0; } - public double [][] interpolateWQ(Row other, double km, int steps) { + public void interpolateW( + Row other, + double km, + double [] iqs, + double [] ows, + WstValueTable table + ) { + double kmWeight = factor(km, this.km, other.km); + QPosition qPosition = new QPosition(); + + for (int i = 0; i < iqs.length; ++i) { + if (table.getQPosition(km, iqs[i], qPosition) != null) { + double wt = getW(qPosition); + double wo = other.getW(qPosition); + ows[i] = weight(kmWeight, wt, wo); + } + else { + ows[i] = Double.NaN; + } + } + } + + public double [][] interpolateWQ( + Row other, + double km, + int steps, + WstValueTable table + ) { int W = Math.min(ws.length, other.ws.length); if (W < 1) { @@ -219,7 +147,11 @@ for (int i = 0; i < W; ++i) { double wws = weight(factor, ws[i], other.ws[i]); - double wqs = weight(factor, qs[i], other.qs[i]); + double wqs = weight( + factor, + table.getQIndex(i, km), + table.getQIndex(i, other.km)); + splineW[i] = wws; splineQ[i] = wqs; if (wqs < minQ) minQ = wqs; @@ -248,7 +180,7 @@ return new double [][] { outWs, outQs }; } - public double [][] interpolateWQ(int steps) { + public double [][] interpolateWQ(int steps, WstValueTable table) { int W = ws.length; @@ -256,17 +188,17 @@ return new double[2][0]; } - double [] splineW = new double[W]; double [] splineQ = new double[W]; double minQ = Double.MAX_VALUE; double maxQ = -Double.MAX_VALUE; + QPosition qPosition = new QPosition(); + for (int i = 0; i < W; ++i) { - splineW[i] = ws[i]; - splineQ[i] = qs[i]; - if (qs[i] < minQ) minQ = qs[i]; - if (qs[i] > maxQ) maxQ = qs[i]; + splineQ[i] = table.getQIndex(i, km); + if (splineQ[i] < minQ) minQ = splineQ[i]; + if (splineQ[i] > maxQ) maxQ = splineQ[i]; } double stepWidth = (maxQ - minQ)/steps; @@ -274,7 +206,7 @@ SplineInterpolator interpolator = new SplineInterpolator(); PolynomialSplineFunction spline = - interpolator.interpolate(splineQ, splineW); + interpolator.interpolate(splineQ, ws); double [] outWs = new double[steps]; double [] outQs = new double[steps]; @@ -292,153 +224,40 @@ return new double [][] { outWs, outQs }; } - public QPosition getQPosition(double q) { - int qIndex = getQIndex(q); - - if (qIndex >= 0) { - // direct column match - return new QPosition(q, qIndex, 1.0); - } - qIndex = -qIndex -1; + public double getW(QPosition qPosition) { + int index = qPosition.index; + double weight = qPosition.weight; - if (qIndex < 0 || qIndex >= qs.length-1) { - // do not extrapolate - return null; - } - - return new QPosition( - q, qIndex, factor(q, qs[qIndex-1], qs[qIndex])); + return weight == 1.0 + ? ws[index] + : weight(weight, ws[index-1], ws[index]); } - public QPosition getQPosition(Row other, double km, double q) { - - QPosition qpt = getQPosition(q); - QPosition qpo = other.getQPosition(q); - - if (qpt == null || qpo == null) { - return null; - } - - double kmWeight = factor(km, this.km, other.km); - - // XXX: Index interpolation is a bit sticky here! - - int index = (int)Math.round( - weight(kmWeight, qpt.index, qpo.index)); - - double weight = weight(kmWeight, qpt.weight, qpo.weight); - - return new QPosition(q, index, weight); - } - - public void storeWQ( - QPosition qPosition, - double [] ows, - double [] oqs, - int index - ) { - int qIdx = qPosition.index; - double qWeight = qPosition.weight; - - if (qWeight == 1.0) { - oqs[index] = qs[qIdx]; - ows[index] = ws[qIdx]; - } - else { - oqs[index] = weight(qWeight, qs[qIdx-1], qs[qIdx]); - ows[index] = weight(qWeight, ws[qIdx-1], ws[qIdx]); - } - } - - public void storeWQ( + public double getW( Row other, double km, - QPosition qPosition, - double [] ows, - double [] oqs, - int index + QPosition qPosition ) { double kmWeight = factor(km, this.km, other.km); - double tq, tw; - double oq, ow; + int index = qPosition.index; + double weight = qPosition.weight; - int qIdx = qPosition.index; - double qWeight = qPosition.weight; + double tw, ow; - if (qWeight == 1.0) { - tq = qs[qIdx]; - tw = ws[qIdx]; - oq = other.qs[qIdx]; - ow = other.ws[qIdx]; + if (weight == 1.0) { + tw = ws[index]; + ow = other.ws[index]; } else { - tq = weight(qWeight, qs[qIdx-1], qs[qIdx]); - tw = weight(qWeight, ws[qIdx-1], ws[qIdx]); - oq = weight(qWeight, other.qs[qIdx-1], other.qs[qIdx]); - ow = weight(qWeight, other.ws[qIdx-1], other.ws[qIdx]); + tw = weight(weight, ws[index-1], ws[index]); + ow = weight(weight, other.ws[index-1], other.ws[index]); } - oqs[index] = weight(kmWeight, tq, oq); - ows[index] = weight(kmWeight, tw, ow); - } - - public void storeWQ( - QPosition qPosition, - LinearRemap remap, - double [] ows, - double [] oqs, - int index - ) { - int qIdx = qPosition.index; - double qWeight = qPosition.weight; - - oqs[index] = remap.eval(km, qWeight == 1.0 - ? qs[qIdx] - : weight(qWeight, qs[qIdx-1], qs[qIdx])); - - ows[index] = interpolateW(oqs[index]); + return weight(kmWeight, tw, ow); } - - public void storeWQ( - Row other, - double km, - QPosition qPosition, - LinearRemap remap, - double [] ows, - double [] oqs, - int index - ) { - double kmWeight = factor(km, this.km, other.km); - - double tq, tw; - double oq, ow; - - int qIdx = qPosition.index; - double qWeight = qPosition.weight; - - if (qWeight == 1.0) { - tq = remap.eval(this.km, qs[qIdx]); - oq = remap.eval(other.km, other.qs[qIdx]); - } - else { - tq = remap.eval( - this.km, - weight(qWeight, qs[qIdx-1], qs[qIdx])); - oq = remap.eval( - other.km, - weight(qWeight, other.qs[qIdx-1], other.qs[qIdx])); - } - - tw = interpolateW(tq); - ow = other.interpolateW(oq); - - oqs[index] = weight(kmWeight, tq, oq); - ows[index] = weight(kmWeight, tw, ow); - } - } - // class Row + } // class Row protected List<Row> rows; @@ -453,23 +272,27 @@ this.columns = columns; } + public WstValueTable(Column [] columns, List<Row> rows) { + this.columns = columns; + this.rows = rows; + } + public void sortRows() { Collections.sort(rows); } - - public double [] interpolateW(double km, double [] qs) { - return interpolateW(km, qs, new double[qs.length]); - } - public double [] interpolateW(double km, double [] qs, double [] ws) { int rowIndex = Collections.binarySearch(rows, new Row(km)); + QPosition qPosition = new QPosition(); + if (rowIndex >= 0) { // direct row match Row row = rows.get(rowIndex); for (int i = 0; i < qs.length; ++i) { - ws[i] = row.interpolateW(qs[i]); + ws[i] = getQPosition(km, qs[i], qPosition) != null + ? row.getW(qPosition) + : Double.NaN; } } else { // needs bilinear interpolation @@ -480,16 +303,15 @@ Arrays.fill(ws, Double.NaN); } else { - rows.get(rowIndex-1).interpolateW( - rows.get(rowIndex), - km, qs, ws); + Row r1 = rows.get(rowIndex-1); + Row r2 = rows.get(rowIndex); + r1.interpolateW(r2, km, qs, ws, this); } } return ws; } - public double [][] interpolateWQ(double km) { return interpolateWQ(km, DEFAULT_Q_STEPS); } @@ -500,7 +322,7 @@ if (rowIndex >= 0) { // direct row match Row row = rows.get(rowIndex); - return row.interpolateWQ(steps); + return row.interpolateWQ(steps, this); } rowIndex = -rowIndex -1; @@ -513,7 +335,7 @@ Row r1 = rows.get(rowIndex-1); Row r2 = rows.get(rowIndex); - return r1.interpolateWQ(r2, km, steps); + return r1.interpolateWQ(r2, km, steps, this); } public void interpolate( @@ -526,13 +348,24 @@ int R1 = rows.size()-1; Row kmKey = new Row(); + + QPosition nPosition = new QPosition(); + for (int i = 0; i < kms.length; ++i) { kmKey.km = kms[i]; + + qs[i] = remap.eval(kms[i], getQ(qPosition, kms[i])); + + if (getQPosition(kms[i], qs[i], nPosition) == null) { + ws[i] = Double.NaN; + continue; + } + int rowIndex = Collections.binarySearch(rows, kmKey); if (rowIndex >= 0) { // direct row match - rows.get(rowIndex).storeWQ(qPosition, remap, ws, qs, i); + ws[i] = rows.get(rowIndex).getW(nPosition); continue; } @@ -540,12 +373,12 @@ if (rowIndex < 1 || rowIndex >= R1) { // do not extrapolate - ws[i] = qs[i] = Double.NaN; + ws[i] = Double.NaN; continue; } Row r1 = rows.get(rowIndex-1); Row r2 = rows.get(rowIndex); - r1.storeWQ(r2, kms[i], qPosition, remap, ws, qs, i); + ws[i] = r1.getW(r2, kms[i], nPosition); } } @@ -562,28 +395,10 @@ int R1 = rows.size()-1; - QPosition qPosition = null; - - if (rowIndex >= 0) { // direct row match - qPosition = rows.get(rowIndex).getQPosition(q); - } - else { - rowIndex = -rowIndex -1; - - if (rowIndex < 1 || rowIndex >= R1) { - // do not extrapolate - return null; - } - - // interpolation case - Row r1 = rows.get(rowIndex-1); - Row r2 = rows.get(rowIndex); - - qPosition = r1.getQPosition(r2, kms[referenceIndex], q); - } + QPosition qPosition = getQPosition(kms[referenceIndex], q); if (qPosition == null) { - // we cannot locate q in reference row + // we cannot locate q at km return null; } @@ -592,9 +407,11 @@ rowIndex = Collections.binarySearch(rows, kmKey); + qs[i] = getQ(qPosition, kms[i]); + if (rowIndex >= 0) { // direct row match - rows.get(rowIndex).storeWQ(qPosition, ws, qs, i); + ws[i] = rows.get(rowIndex).getW(qPosition); continue; } @@ -602,15 +419,102 @@ if (rowIndex < 1 || rowIndex >= R1) { // do not extrapolate - ws[i] = qs[i] = Double.NaN; + ws[i] = Double.NaN; continue; } Row r1 = rows.get(rowIndex-1); Row r2 = rows.get(rowIndex); - r1.storeWQ(r2, kms[i], qPosition, ws, qs, i); + + ws[i] = r1.getW(r2, kms[i], qPosition); } return qPosition; } + + public QPosition getQPosition(double km, double q) { + return getQPosition(km, q, new QPosition()); + } + + public QPosition getQPosition(double km, double q, QPosition qPosition) { + + if (columns.length == 0) { + return null; + } + + double qLast = columns[0].getQRangeTree().findQ(km); + + if (Math.abs(qLast - q) < 0.00001) { + return qPosition.set(0, 1d); + } + + for (int i = 1; i < columns.length; ++i) { + double qCurrent = columns[i].getQRangeTree().findQ(km); + if (Math.abs(qCurrent - q) < 0.00001) { + return qPosition.set(i, 1d); + } + + double qMin, qMax; + if (qLast < qCurrent) { qMin = qLast; qMax = qCurrent; } + else { qMin = qCurrent; qMax = qLast; } + + if (q > qMin && q < qMax) { + double weight = factor(q, qLast, qCurrent); + return qPosition.set(i, weight); + } + qLast = qCurrent; + } + + return null; + } + + public double getQIndex(int index, double km) { + return columns[index].getQRangeTree().findQ(km); + } + + public double getQ(QPosition qPosition, double km) { + int index = qPosition.index; + double weight = qPosition.weight; + + if (weight == 1d) { + return columns[index].getQRangeTree().findQ(km); + } + double q1 = columns[index-1].getQRangeTree().findQ(km); + double q2 = columns[index ].getQRangeTree().findQ(km); + return weight(weight, q1, q2); + } + + public static final double linear( + double x, + double x1, double x2, + double y1, double y2 + ) { + // y1 = m*x1 + b + // y2 = m*x2 + b + // y2 - y1 = m*(x2 - x1) + // m = (y2 - y1)/(x2 - x1) # x2 != x1 + // b = y1 - m*x1 + + if (x1 == x2) { + return 0.5*(y1 + y2); + } + double m = (y2 - y1)/(x2 - x1); + double b = y1 - m*x1; + return x*m + b; + } + + public static final double factor(double x, double p1, double p2) { + // 0 = m*p1 + b <=> b = -m*p1 + // 1 = m*p2 + b + // 1 = m*(p2 - p1) + // m = 1/(p2 - p1) # p1 != p2 + // f(x) = x/(p2-p1) - p1/(p2-p1) <=> (x-p1)/(p2-p1) + + return p1 == p2 ? 0.0 : (x-p1)/(p2-p1); + } + + public static final double weight(double factor, double a, double b) { + //return (1.0-factor)*a + factor*b; + return a + factor*(b-a); + } } // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :