sascha@326: package de.intevation.flys.artifacts.model; sascha@326: sascha@326: import java.io.Serializable; sascha@326: sascha@326: import de.intevation.flys.model.River; sascha@326: import de.intevation.flys.model.Wst; sascha@326: import de.intevation.flys.model.WstColumn; sascha@326: sascha@326: import de.intevation.flys.backend.SessionHolder; sascha@326: sascha@335: import java.util.Arrays; sascha@326: import java.util.ArrayList; sascha@326: import java.util.Comparator; sascha@326: import java.util.List; sascha@326: import java.util.Collections; sascha@326: import java.util.Iterator; 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: import org.hibernate.Session; sascha@326: import org.hibernate.Query; sascha@326: import org.hibernate.SQLQuery; sascha@326: sascha@326: import org.hibernate.type.StandardBasicTypes; sascha@326: sascha@326: public class WstValueTable sascha@326: implements Serializable sascha@326: { sascha@339: private static Logger log = Logger.getLogger(WstValueTable.class); sascha@326: sascha@326: // TODO: put this into a property file sascha@326: public static final String SQL_POS_WQ = sascha@327: "SELECT position, w, q, column_pos" + sascha@327: " FROM wst_value_table" + sascha@327: " WHERE wst_id = :wst_id"; sascha@326: sascha@339: public static final double DEFAULT_STEP_WIDTH = 0.01; sascha@339: sascha@326: public static class Column sascha@326: implements Serializable sascha@326: { sascha@326: protected String name; sascha@326: 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@326: } sascha@326: // class Column sascha@326: sascha@326: public static class Row sascha@335: implements Serializable, Comparable sascha@326: { sascha@326: double km; sascha@335: double [] ws; sascha@335: double [] qs; sascha@335: boolean qSorted; 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@335: public Row(double km, double [] ws, double [] qs) { sascha@335: this(km); sascha@335: this.ws = ws; sascha@335: this.qs = qs; sascha@335: } sascha@335: sascha@335: public static final double linear( sascha@335: double x, sascha@335: double x1, double x2, sascha@335: double y1, double y2 sascha@335: ) { sascha@335: // y1 = m*x1 + b sascha@335: // y2 = m*x2 + b sascha@335: // y2 - y1 = m*(x2 - x1) sascha@335: // m = (y2 - y1)/(x2 - x1) # x2 != x1 sascha@335: // b = y1 - m*x1 sascha@335: sascha@335: if (x1 == x2) { sascha@335: return 0.5*(y1 + y2); sascha@335: } sascha@335: double m = (y2 - y1)/(x2 - x1); sascha@335: double b = y1 - m*x1; sascha@335: return x*m + b; sascha@335: } sascha@335: sascha@335: public static final double factor(double x, double p1, double p2) { sascha@335: // 0 = m*p1 + b <=> b = -m*p1 sascha@335: // 1 = m*p2 + b sascha@335: // 1 = m*(p2 - p1) sascha@335: // m = 1/(p2 - p1) # p1 != p2 sascha@335: // f(x) = x/(p2-p1) - p1/(p2-p1) <=> (x-p1)/(p2-p1) sascha@335: sascha@335: return p1 == p2 ? 0.0 : (x-p1)/(p2-p1); sascha@335: } sascha@335: sascha@335: public static final double weight(double factor, double a, double b) { sascha@335: return (1.0-factor)*a + factor*b; sascha@335: } sascha@335: sascha@335: public double getWForKM(Row other, int index, double km) { sascha@335: double w1 = ws[index]; sascha@335: double w2 = other.ws[index]; sascha@335: double km1 = this.km; sascha@335: double km2 = other.km; sascha@335: return linear(km, km1, km2, w1, w2); sascha@335: } sascha@335: sascha@335: public void interpolateW( sascha@335: Row other, sascha@335: double km, sascha@335: double [] input, sascha@335: double [] output sascha@335: ) { sascha@335: double factor = factor(km, this.km, other.km); sascha@335: sascha@335: for (int i = 0; i < input.length; ++i) { sascha@335: double q = input[i]; sascha@335: int idx1 = getQIndex(q); sascha@335: int idx2 = other.getQIndex(q); sascha@335: sascha@335: double w1 = idx1 >= 0 sascha@335: ? ws[idx1] sascha@335: : interpolateW(-idx1-1, q); sascha@335: sascha@335: double w2 = idx2 >= 0 sascha@335: ? other.ws[idx2] sascha@335: : other.interpolateW(-idx2-1, q); sascha@335: sascha@335: output[i] = weight(factor, w1, w2); sascha@335: } sascha@335: } sascha@335: sascha@335: public double interpolateW(int idx, double q) { sascha@335: return idx < 1 || idx >= qs.length sascha@335: ? Double.NaN // do not extrapolate sascha@335: : linear(q, qs[idx-1], qs[idx], ws[idx-1], ws[idx]); sascha@335: } sascha@335: sascha@335: public double interpolateW(double q) { sascha@335: if (Double.isNaN(q)) return Double.NaN; sascha@335: int index = getQIndex(q); ingo@360: return index >= 0 ? ws[index] : interpolateW(-index -1, q); sascha@335: } sascha@335: sascha@335: public int getQIndex(double q) { sascha@335: return qSorted ? binaryQIndex(q) : linearQIndex(q); sascha@335: } sascha@335: sascha@335: public int binaryQIndex(double q) { sascha@335: return Arrays.binarySearch(qs, q); sascha@335: } sascha@335: sascha@335: public int linearQIndex(double q) { sascha@335: switch (qs.length) { sascha@335: case 0: break; sascha@335: case 1: sascha@335: if (qs[0] == q) return 0; sascha@335: if (qs[0] < q) return -(1+1); sascha@335: return -(0+1); sascha@335: default: sascha@335: for (int i = 1; i < qs.length; ++i) { sascha@335: double qa = qs[i-1]; sascha@335: double qb = qs[i]; sascha@335: if (qa == q) return i-1; sascha@335: if (qb == q) return i; sascha@335: if (qa > qb) { double t = qa; qa = qb; qb = t; } sascha@335: if (q > qa && q < qb) return -(i+1); sascha@335: } sascha@335: return -qs.length; sascha@335: } sascha@335: sascha@335: return -1; sascha@335: } sascha@335: sascha@335: public int compareTo(Row other) { sascha@335: double d = km - other.km; sascha@335: if (d < 0.0) return -1; sascha@335: if (d > 0.0) return +1; sascha@335: return 0; sascha@326: } sascha@339: sascha@339: public double [][] cloneWQs() { sascha@339: return new double [][] { sascha@339: (double [])ws.clone(), sascha@344: (double [])qs.clone() }; sascha@339: } sascha@339: sascha@339: public double [][] interpolateWQ(Row other, double km, double stepWidth) { sascha@339: sascha@339: int W1 = ascendingWs(); sascha@339: int W2 = other.ascendingWs(); sascha@339: sascha@339: int W = Math.min(W1, W2); sascha@339: sascha@339: if (W < 1) { sascha@339: return new double[2][0]; sascha@339: } sascha@339: sascha@339: double factor = factor(km, this.km, other.km); sascha@339: sascha@339: double minW = weight(factor, ws[0], other.ws[0]); sascha@339: double maxW = weight(factor, ws[W-1], other.ws[W-1]); sascha@339: sascha@339: if (minW > maxW) { sascha@339: double t = minW; sascha@339: minW = maxW; sascha@339: maxW = t; sascha@339: } sascha@339: double [] x = new double[W]; sascha@339: double [] y = new double[W]; sascha@339: sascha@339: for (int i = 0; i < W; ++i) { sascha@339: x[i] = weight(factor, ws[i], other.ws[i]); sascha@339: y[i] = weight(factor, qs[i], other.qs[i]); sascha@339: } sascha@339: sascha@339: SplineInterpolator interpolator = new SplineInterpolator(); sascha@339: PolynomialSplineFunction spline = interpolator.interpolate(x, y); sascha@339: sascha@339: double [] outWs = sascha@339: new double[(int)Math.ceil((maxW - minW)/stepWidth)]; sascha@339: double [] outQs = sascha@339: new double[outWs.length]; sascha@339: sascha@339: try { sascha@339: double w = minW; sascha@339: for (int i = 0; i < outWs.length; ++i, w += stepWidth) { sascha@339: outQs[i] = spline.value(outWs[i] = w); 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@339: sascha@339: public double [][] interpolateWQ(double stepWidth) { sascha@339: int W = ascendingWs(); // ignore back jumps sascha@339: sascha@339: if (W < 1) { sascha@339: return new double[2][0]; sascha@339: } sascha@339: sascha@339: double [] x = new double[W]; sascha@339: double [] y = new double[W]; sascha@339: sascha@339: for (int i = 0; i < W; ++i) { sascha@339: x[i] = ws[i]; sascha@339: y[i] = qs[i]; sascha@339: } sascha@339: sascha@339: SplineInterpolator interpolator = new SplineInterpolator(); sascha@339: sascha@339: PolynomialSplineFunction spline = interpolator.interpolate(x, y); sascha@339: sascha@339: double minW = ws[0]; sascha@339: double maxW = ws[W-1]; sascha@339: sascha@339: double [] outWs = sascha@339: new double[(int)Math.ceil((maxW - minW)/stepWidth)]; sascha@339: double [] outQs = sascha@339: new double[outWs.length]; sascha@339: sascha@339: try { sascha@339: double w = minW; sascha@339: for (int i = 0; i < outWs.length; ++i, w += stepWidth) { sascha@339: outQs[i] = spline.value(outWs[i] = w); 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@339: public int ascendingWs() { sascha@339: if (ws.length < 2) { sascha@339: return ws.length; sascha@339: } sascha@339: sascha@339: int idx = 1; sascha@339: sascha@339: for (; idx < ws.length; ++idx) { sascha@339: if (Double.isNaN(ws[idx]) || ws[idx] < ws[idx-1]) { sascha@339: return idx; sascha@339: } sascha@339: } sascha@339: sascha@339: return idx; sascha@339: } sascha@339: sascha@339: public double [][] weightWQs(Row other, double km) { sascha@339: int W = Math.min(ws.length, other.ws.length); sascha@339: double factor = factor(km, this.km, other.km); sascha@339: sascha@339: double [] outWs = new double[W]; sascha@339: double [] outQs = new double[W]; sascha@339: sascha@339: for (int i = 0; i < W; ++i) { sascha@339: outWs[i] = weight(factor, ws[i], other.ws[i]); sascha@339: outQs[i] = weight(factor, qs[i], other.qs[i]); sascha@339: } sascha@339: sascha@339: return new double [][] { outWs, outQs }; sascha@339: } sascha@326: } sascha@326: // 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@335: sascha@335: public double [] interpolateW(double km, double [] qs) { sascha@335: sascha@335: double [] ws = new double[qs.length]; sascha@335: sascha@335: int rowIndex = Collections.binarySearch(rows, new Row(km)); sascha@335: 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@335: ws[i] = row.interpolateW(qs[i]); 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@335: rows.get(rowIndex-1).interpolateW( sascha@335: rows.get(rowIndex), sascha@335: km, qs, ws); sascha@335: } sascha@335: } sascha@335: sascha@335: return ws; sascha@335: } sascha@335: sascha@339: sascha@339: public double [][] interpolateWQ(double km) { sascha@339: return interpolateWQ(km, DEFAULT_STEP_WIDTH, true); sascha@339: } sascha@339: sascha@339: public double [][] interpolateWQ( sascha@339: double km, sascha@339: double stepWidth, sascha@339: boolean checkAscending sascha@339: ) { 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@339: return checkAscending sascha@339: ? row.interpolateWQ(stepWidth) sascha@339: : row.cloneWQs(); 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@339: return checkAscending sascha@339: ? r1.interpolateWQ(r2, km, stepWidth) sascha@339: : r1.weightWQs(r2, km); sascha@339: } sascha@339: sascha@326: public static WstValueTable getTable(River river) { sascha@326: return getTable(river, 0); sascha@326: } sascha@326: sascha@326: public static WstValueTable getTable(River river, int kind) { sascha@326: sascha@326: Session session = SessionHolder.HOLDER.get(); sascha@326: sascha@326: Query query = session.createQuery( sascha@326: "from Wst where river=:river and kind=:kind"); sascha@326: query.setParameter("river", river); sascha@326: query.setInteger("kind", kind); sascha@326: sascha@326: List wsts = query.list(); sascha@326: sascha@326: if (wsts.isEmpty()) { sascha@326: return null; sascha@326: } sascha@326: sascha@326: Wst wst = wsts.get(0); sascha@326: sascha@326: // TODO: Do this sorting at database level sascha@326: List wstColumns = new ArrayList(wst.getColumns()); sascha@326: Collections.sort(wstColumns, new Comparator() { sascha@326: public int compare(WstColumn a, WstColumn b) { sascha@326: int pa = a.getPosition(); sascha@326: int pb = b.getPosition(); sascha@326: if (pa < pb) return -1; sascha@326: if (pa > pb) return +1; sascha@326: return 0; sascha@326: } sascha@326: }); sascha@326: sascha@326: Column [] columns = new Column[wstColumns.size()]; sascha@326: for (int i = 0; i < columns.length; ++i) { sascha@326: columns[i] = new Column(wstColumns.get(i).getName()); sascha@326: } sascha@326: sascha@326: // using native SQL here to avoid myriad of small objects. sascha@326: SQLQuery sqlQuery = session.createSQLQuery(SQL_POS_WQ) sascha@326: .addScalar("position", StandardBasicTypes.DOUBLE) sascha@326: .addScalar("w", StandardBasicTypes.DOUBLE) sascha@326: .addScalar("q", StandardBasicTypes.DOUBLE) sascha@326: .addScalar("column_pos", StandardBasicTypes.INTEGER); sascha@326: sascha@326: sqlQuery.setInteger("wst_id", wst.getId()); sascha@326: sascha@326: WstValueTable valueTable = new WstValueTable(columns); sascha@326: sascha@326: int lastColumnNo = -1; sascha@326: Row row = null; sascha@326: sascha@335: Double lastQ = -Double.MAX_VALUE; sascha@335: boolean qSorted = true; sascha@335: ingo@360: for (Object r: sqlQuery.list()) { ingo@360: Object[] result = (Object[]) r; ingo@360: sascha@326: double km = (Double) result[0]; sascha@326: Double w = (Double) result[1]; sascha@326: Double q = (Double) result[2]; sascha@326: int columnNo = (Integer)result[3]; sascha@326: sascha@326: if (columnNo > lastColumnNo) { // new row sascha@326: if (row != null) { sascha@335: row.qSorted = qSorted; sascha@326: valueTable.rows.add(row); sascha@326: } sascha@335: row = new Row( sascha@335: km, sascha@335: new double[columnNo+1], sascha@335: new double[columnNo+1]); sascha@335: lastQ = -Double.MAX_VALUE; sascha@335: qSorted = true; sascha@326: } sascha@326: sascha@335: row.ws[columnNo] = w != null ? w : Double.NaN; sascha@335: row.qs[columnNo] = q != null ? q : Double.NaN; sascha@335: sascha@335: if (qSorted && (q == null || lastQ > q)) { sascha@335: qSorted = false; sascha@335: } sascha@335: lastQ = q; sascha@326: sascha@326: lastColumnNo = columnNo; sascha@326: } sascha@326: sascha@326: if (row != null) { sascha@326: valueTable.rows.add(row); sascha@326: } sascha@326: sascha@326: return valueTable; sascha@326: } sascha@326: } sascha@326: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :