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<Row>
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);
sascha@335:             return index >= 0 ? ws[index] : interpolateW(-index -1);
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<Row> rows;
sascha@326: 
sascha@326:     protected Column [] columns;
sascha@326: 
sascha@326:     public WstValueTable() {
sascha@326:         rows = new ArrayList<Row>();
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<Wst> 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<WstColumn> wstColumns = new ArrayList(wst.getColumns());
sascha@326:         Collections.sort(wstColumns, new Comparator<WstColumn>() {
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: 
sascha@326:         for (Iterator<Object []> iter = sqlQuery.iterate(); iter.hasNext();) {
sascha@326:             Object [] result = iter.next();
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 :