Mercurial > dive4elements > river
view flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java @ 1839:c3b44da22a3a
Minor refactoring in preparation of generation of interpolators (WstValueTable) for static WQKms.
flys-artifacts/trunk@3175 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author | Felix Wolfsteller <felix.wolfsteller@intevation.de> |
---|---|
date | Mon, 07 Nov 2011 12:37:41 +0000 |
parents | 853cd2120d69 |
children | ddd425858169 |
line wrap: on
line source
package de.intevation.flys.artifacts.model; import java.io.Serializable; import de.intevation.flys.artifacts.math.Linear; import de.intevation.flys.artifacts.math.Function; import java.util.Arrays; import java.util.ArrayList; import java.util.List; import java.util.Collections; import org.apache.log4j.Logger; import org.apache.commons.math.analysis.interpolation.SplineInterpolator; import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction; import org.apache.commons.math.ArgumentOutsideDomainException; import org.apache.commons.math.exception.MathIllegalArgumentException; /** * W, Q and km data from database 'wsts' spiced with interpolation algorithms. */ public class WstValueTable implements Serializable { private static Logger log = Logger.getLogger(WstValueTable.class); public static final int DEFAULT_Q_STEPS = 500; /** * A Column in the table, typically representing one measurement session. */ public static final class Column implements Serializable { protected String name; protected QRangeTree qRangeTree; public Column() { } public Column(String name) { this.name = name; } public String getName() { return name; } public void setName(String name) { this.name = name; } public QRangeTree getQRangeTree() { return qRangeTree; } public void setQRangeTree(QRangeTree qRangeTree) { this.qRangeTree = qRangeTree; } } // class Column /** * A (weighted) position used for interpolation. */ public static final class QPosition { protected int index; protected double weight; public QPosition() { } 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 /** * A row, typically a position where measurements were taken. */ public static final class Row implements Serializable, Comparable<Row> { double km; double [] ws; public Row() { } public Row(double km) { this.km = km; } public Row(double km, double [] ws) { this(km); this.ws = ws; } /** * Compare according to place of measurement (km). */ public int compareTo(Row other) { double d = km - other.km; if (d < -0.0001) return -1; if (d > 0.0001) return +1; return 0; } public void interpolateW( Row other, double km, double [] iqs, double [] ows, WstValueTable table, Calculation errors ) { double kmWeight = Linear.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); if (Double.isNaN(wt) || Double.isNaN(wo)) { if (errors != null) { // TODO: I18N errors.addProblem( km, "cannot find w for q = " + iqs[i]); } ows[i] = Double.NaN; } else { ows[i] = Linear.weight(kmWeight, wt, wo); } } else { if (errors != null) { // TODO: I18N errors.addProblem(km, "cannot find q = " + iqs[i]); } ows[i] = Double.NaN; } } } public double [][] interpolateWQ( Row other, double km, int steps, WstValueTable table, Calculation errors ) { int W = Math.min(ws.length, other.ws.length); if (W < 1) { if (errors != null) { // TODO: I18N errors.addProblem("no ws found"); } return new double[2][0]; } double factor = Linear.factor(km, this.km, other.km); double [] splineQ = new double[W]; double [] splineW = new double[W]; double minQ = Double.MAX_VALUE; double maxQ = -Double.MAX_VALUE; for (int i = 0; i < W; ++i) { double wws = Linear.weight(factor, ws[i], other.ws[i]); double wqs = Linear.weight( factor, table.getQIndex(i, km), table.getQIndex(i, other.km)); if (Double.isNaN(wws) || Double.isNaN(wqs)) { if (errors != null) { // TODO: I18N errors.addProblem(km, "cannot find w or q"); } } else { if (wqs < minQ) minQ = wqs; if (wqs > maxQ) maxQ = wqs; } splineW[i] = wws; splineQ[i] = wqs; } double stepWidth = (maxQ - minQ)/steps; SplineInterpolator interpolator = new SplineInterpolator(); PolynomialSplineFunction spline; try { spline = interpolator.interpolate(splineQ, splineW); } catch (MathIllegalArgumentException miae) { if (errors != null) { // TODO: I18N errors.addProblem(km, "spline creation failed"); } log.error("spline creation failed"); return new double[2][0]; } double [] outWs = new double[steps]; double [] outQs = new double[steps]; Arrays.fill(outWs, Double.NaN); Arrays.fill(outQs, Double.NaN); try { double q = minQ; for (int i = 0; i < outWs.length; ++i, q += stepWidth) { outWs[i] = spline.value(outQs[i] = q); } } catch (ArgumentOutsideDomainException aode) { if (errors != null) { // TODO: I18N errors.addProblem(km, "spline interpolation failed"); } log.error("spline interpolation failed", aode); } return new double [][] { outWs, outQs }; } public double [][] interpolateWQ( int steps, WstValueTable table, Calculation errors ) { int W = ws.length; if (W < 1) { if (errors != null) { // TODO: I18N errors.addProblem(km, "no ws found"); } return new double[2][0]; } double [] splineQ = new double[W]; double minQ = Double.MAX_VALUE; double maxQ = -Double.MAX_VALUE; for (int i = 0; i < W; ++i) { double sq = table.getQIndex(i, km); if (Double.isNaN(sq)) { if (errors != null) { // TODO: I18N errors.addProblem( km, "no q found in " + (i+1) + " column"); } } else { if (sq < minQ) minQ = sq; if (sq > maxQ) maxQ = sq; } splineQ[i] = sq; } double stepWidth = (maxQ - minQ)/steps; SplineInterpolator interpolator = new SplineInterpolator(); PolynomialSplineFunction spline; try { spline = interpolator.interpolate(splineQ, ws); } catch (MathIllegalArgumentException miae) { if (errors != null) { // TODO: I18N errors.addProblem(km, "spline creation failed"); } log.error("spline creation failed"); return new double[2][0]; } double [] outWs = new double[steps]; double [] outQs = new double[steps]; Arrays.fill(outWs, Double.NaN); Arrays.fill(outQs, Double.NaN); try { double q = minQ; for (int i = 0; i < outWs.length; ++i, q += stepWidth) { outWs[i] = spline.value(outQs[i] = q); } } catch (ArgumentOutsideDomainException aode) { if (errors != null) { // TODO: I18N errors.addProblem(km, "spline interpolation failed"); } log.error("spline interpolation failed.", aode); } return new double [][] { outWs, outQs }; } public double getW(QPosition qPosition) { int index = qPosition.index; double weight = qPosition.weight; return weight == 1.0 ? ws[index] : Linear.weight(weight, ws[index-1], ws[index]); } public double getW( Row other, double km, QPosition qPosition ) { double kmWeight = Linear.factor(km, this.km, other.km); int index = qPosition.index; double weight = qPosition.weight; double tw, ow; if (weight == 1.0) { tw = ws[index]; ow = other.ws[index]; } else { tw = Linear.weight(weight, ws[index-1], ws[index]); ow = Linear.weight(weight, other.ws[index-1], other.ws[index]); } return Linear.weight(kmWeight, tw, ow); } } // class Row /** Rows in table. */ protected List<Row> rows; /** Columns in table. */ protected Column [] columns; public WstValueTable() { rows = new ArrayList<Row>(); } public WstValueTable(Column [] columns) { this(); this.columns = columns; } public WstValueTable(Column [] columns, List<Row> rows) { this.columns = columns; this.rows = rows; } /** Sort rows (by km). */ public void sortRows() { Collections.sort(rows); } public double [] interpolateW(double km, double [] qs, double [] ws) { return interpolateW(km, qs, ws, null); } /** * @param ws (output parameter), gets returned. * @return output parameter ws. */ public double [] interpolateW( double km, double [] qs, double [] ws, Calculation errors ) { 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) { if (getQPosition(km, qs[i], qPosition) == null) { if (errors != null) { // TODO: I18N errors.addProblem(km, "cannot find q = " + qs[i]); } ws[i] = Double.NaN; } else { if (Double.isNaN(ws[i] = row.getW(qPosition)) && errors != null) { // TODO: I18N errors.addProblem( km, "cannot find w for q = " + qs[i]); } } } } else { // needs bilinear interpolation rowIndex = -rowIndex -1; if (rowIndex < 1 || rowIndex >= rows.size()) { // do not extrapolate Arrays.fill(ws, Double.NaN); if (errors != null) { // TODO: I18N errors.addProblem(km, "km not found"); } } else { Row r1 = rows.get(rowIndex-1); Row r2 = rows.get(rowIndex); r1.interpolateW(r2, km, qs, ws, this, errors); } } return ws; } /** * Interpolate W and Q values at a given km. */ public double [][] interpolateWQ(double km) { return interpolateWQ(km, null); } /** * Interpolate W and Q values at a given km. */ public double [][] interpolateWQ(double km, Calculation errors) { return interpolateWQ(km, DEFAULT_Q_STEPS, errors); } public double [][] interpolateWQ(double km, int steps, Calculation errors) { int rowIndex = Collections.binarySearch(rows, new Row(km)); if (rowIndex >= 0) { // direct row match Row row = rows.get(rowIndex); return row.interpolateWQ(steps, this, errors); } rowIndex = -rowIndex -1; if (rowIndex < 1 || rowIndex >= rows.size()) { // do not extrapolate if (errors != null) { // TODO: I18N errors.addProblem(km, "km not found"); } return new double[2][0]; } Row r1 = rows.get(rowIndex-1); Row r2 = rows.get(rowIndex); return r1.interpolateWQ(r2, km, steps, this, errors); } public boolean interpolate( double km, double [] out, QPosition qPosition, Function qFunction ) { int R1 = rows.size()-1; out[1] = qFunction.value(getQ(qPosition, km)); if (Double.isNaN(out[1])) { return false; } QPosition nPosition = new QPosition(); if (getQPosition(km, out[1], nPosition) == null) { return false; } int rowIndex = Collections.binarySearch(rows, new Row(km)); if (rowIndex >= 0) { // direct row match out[0] = rows.get(rowIndex).getW(nPosition); return !Double.isNaN(out[0]); } rowIndex = -rowIndex -1; if (rowIndex < 1 || rowIndex > R1) { // do not extrapolate return false; } Row r1 = rows.get(rowIndex-1); Row r2 = rows.get(rowIndex); out[0] = r1.getW(r2, km, nPosition); return !Double.isNaN(out[0]); } /** * Look up interpolation of a Q at given positions. * * @param q the non-interpolated Q value. * @param referenceKm the reference km (e.g. gauge position). * @param kms positions for which to interpolate. * @param ws (output) resulting interpolated ws. * @param qs (output) resulting interpolated qs. * @param errors calculation object to store errors. */ public QPosition interpolate( double q, double referenceKm, double [] kms, double [] ws, double [] qs, Calculation errors ) { return interpolate( q, referenceKm, kms, ws, qs, 0, kms.length, errors); } public QPosition interpolate( double q, double referenceKm, double [] kms, double [] ws, double [] qs, int startIndex, int length, Calculation errors ) { QPosition qPosition = getQPosition(referenceKm, q); if (qPosition == null) { // we cannot locate q at km Arrays.fill(ws, Double.NaN); Arrays.fill(qs, Double.NaN); if (errors != null) { errors.addProblem(referenceKm, "cannot find q"); } return null; } Row kmKey = new Row(); int R1 = rows.size()-1; for (int i = startIndex, end = startIndex+length; i < end; ++i) { if (Double.isNaN(qs[i] = getQ(qPosition, kms[i]))) { if (errors != null) { errors.addProblem(kms[i], "cannot find q"); } ws[i] = Double.NaN; continue; } kmKey.km = kms[i]; int rowIndex = Collections.binarySearch(rows, kmKey); if (rowIndex >= 0) { // direct row match if (Double.isNaN(ws[i] = rows.get(rowIndex).getW(qPosition)) && errors != null) { errors.addProblem(kms[i], "cannot find w"); } continue; } rowIndex = -rowIndex -1; if (rowIndex < 1 || rowIndex > R1) { // do not extrapolate if (errors != null) { errors.addProblem(kms[i], "cannot find km"); } ws[i] = Double.NaN; continue; } Row r1 = rows.get(rowIndex-1); Row r2 = rows.get(rowIndex); if (Double.isNaN(ws[i] = r1.getW(r2, kms[i], qPosition)) && errors != null) { errors.addProblem(kms[i], "cannot find w"); } } 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 = Linear.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 Linear.weight(weight, q1, q2); } } // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :