mschaefer@9176: /** Copyright (C) 2017 by Bundesanstalt für Gewässerkunde
mschaefer@9176:  * Software engineering by
mschaefer@9176:  *  Björnsen Beratende Ingenieure GmbH
mschaefer@9176:  *  Dr. Schumacher Ingenieurbüro für Wasser und Umwelt
mschaefer@9176:  *
mschaefer@9176:  * This file is Free Software under the GNU AGPL (>=v3)
mschaefer@9176:  * and comes with ABSOLUTELY NO WARRANTY! Check out the
mschaefer@9176:  * documentation coming with Dive4Elements River for details.
mschaefer@9176:  */
mschaefer@9176: package org.dive4elements.river.artifacts.sinfo.common;
mschaefer@9176: 
mschaefer@9176: import java.util.ArrayList;
mschaefer@9176: import java.util.List;
mschaefer@9176: import java.util.Map.Entry;
mschaefer@9176: import java.util.NavigableMap;
mschaefer@9176: import java.util.TreeMap;
mschaefer@9176: 
mschaefer@9176: import org.apache.commons.math.FunctionEvaluationException;
mschaefer@9176: import org.apache.commons.math.analysis.interpolation.LinearInterpolator;
mschaefer@9176: import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
mschaefer@9176: import org.dive4elements.river.artifacts.math.Linear;
mschaefer@9176: import org.dive4elements.river.artifacts.model.Calculation;
mschaefer@9176: import org.dive4elements.river.backend.SessionHolder;
mschaefer@9176: import org.dive4elements.river.model.River;
mschaefer@9176: import org.hibernate.SQLQuery;
mschaefer@9176: import org.hibernate.Session;
mschaefer@9176: import org.hibernate.type.StandardBasicTypes;
mschaefer@9176: 
mschaefer@9176: import gnu.trove.TDoubleArrayList;
mschaefer@9176: 
mschaefer@9176: /**
mschaefer@9176:  * Loading and search/interpolation of a W/Q base table of a river
mschaefer@9176:  *
mschaefer@9176:  * @author Matthias Schäfer
mschaefer@9176:  *
mschaefer@9176:  */
mschaefer@9176: public final class WQBaseTableFinder {
mschaefer@9176: 
mschaefer@9176:     /***** FIELDS *****/
mschaefer@9176: 
mschaefer@9176:     // private static Logger log = Logger.getLogger(WQTableFinder.class);
mschaefer@9176: 
mschaefer@9176:     private final River river;
mschaefer@9176: 
mschaefer@9176:     private Calculation problems;
mschaefer@9176: 
mschaefer@9176:     private final List<String> columnNames;
mschaefer@9176: 
mschaefer@9176:     private final NavigableMap<Double, PolynomialSplineFunction> kmWs;
mschaefer@9176: 
mschaefer@9176:     private final NavigableMap<Double, PolynomialSplineFunction> kmQs;
mschaefer@9176: 
mschaefer@9176:     private static final String SQLCOLUMNS = "SELECT wc.position AS colindex, wc.name AS qzone"
mschaefer@9176:             + " FROM wsts w"
mschaefer@9176:             + "  INNER JOIN wst_columns wc ON w.id=wc.wst_id"
mschaefer@9176:             + " WHERE w.river_id = :river_id"
mschaefer@9176:             + "  AND w.kind = 0"
mschaefer@9176:             + " ORDER BY wc.position ASC";
mschaefer@9176: 
mschaefer@9176:     private static final String SQLMAIN = "SELECT wcv.position AS station, wc.position AS colindex, wcv.w, wqr.q"
mschaefer@9176:             + " FROM wsts w"
mschaefer@9176:             + "  INNER JOIN wst_columns wc ON w.id=wc.wst_id"
mschaefer@9176:             + "  INNER JOIN wst_column_values wcv ON wc.id=wcv.wst_column_id"
mschaefer@9176:             + "  INNER JOIN wst_column_q_ranges wcqr ON wc.id=wcqr.wst_column_id"
mschaefer@9176:             + "  INNER JOIN wst_q_ranges wqr ON wcqr.wst_q_range_id=wqr.id"
mschaefer@9176:             + "  INNER JOIN ranges r ON wqr.range_id=r.id AND wcv.position BETWEEN r.a AND r.b"
mschaefer@9176:             + " WHERE w.river_id = :river_id"
mschaefer@9176:             + "  AND w.kind = 0"
mschaefer@9176:             + "  AND r.river_id = :river_id"
mschaefer@9176:             + "  AND wcv.position BETWEEN :kmfrom - 1 AND :kmto + 1" // some tolerance for start and end of list
mschaefer@9176:             + " ORDER BY wcv.position ASC, "
mschaefer@9176:             + "  wc.position ASC";
mschaefer@9176: 
mschaefer@9176: 
mschaefer@9176:     /***** CONSTRUCTORS *****/
mschaefer@9176: 
mschaefer@9176:     private WQBaseTableFinder(final River river, final Calculation problems, final List<Object[]> colnames, final List<Object[]> rows) {
mschaefer@9176:         this.river = river;
mschaefer@9176:         this.problems = problems;
mschaefer@9176:         this.columnNames = new ArrayList<>();
mschaefer@9176:         for (final Object[] colname : colnames)
mschaefer@9176:             this.columnNames.add(colname[1].toString());
mschaefer@9176:         this.kmWs = new TreeMap<>();
mschaefer@9176:         this.kmQs = new TreeMap<>();
mschaefer@9176:         final TDoubleArrayList ws = new TDoubleArrayList();
mschaefer@9176:         final TDoubleArrayList qs = new TDoubleArrayList();
mschaefer@9176:         double km = Double.NaN;
mschaefer@9176:         for (int i = 0; i <= rows.size() - 1; i++) {
mschaefer@9176:             if (ws.isEmpty() || ((double) rows.get(i)[0] <= km + 0.0001)) {
mschaefer@9176:                 if (ws.isEmpty())
mschaefer@9176:                     km = (double) rows.get(i)[0];
mschaefer@9176:                 ws.add((double) rows.get(i)[2]);
mschaefer@9176:                 qs.add((double) rows.get(i)[3]);
mschaefer@9176:             }
mschaefer@9176:             if ((i == rows.size() - 1) || ((double) rows.get(i)[0] > km + 0.0001)) {
mschaefer@9176:                 try {
mschaefer@9176:                     this.kmWs.put(km, new LinearInterpolator().interpolate(ws.toNativeArray(), qs.toNativeArray()));
mschaefer@9176:                     this.kmQs.put(km, new LinearInterpolator().interpolate(qs.toNativeArray(), ws.toNativeArray()));
mschaefer@9176:                 }
mschaefer@9176:                 catch (final Exception e) {
mschaefer@9176:                     if (this.problems != null) {
mschaefer@9176:                         this.problems.addProblem(km, "wq_base_data.missing");
mschaefer@9176:                         // Report only once
mschaefer@9176:                         this.problems = null;
mschaefer@9176:                     }
mschaefer@9176:                 }
mschaefer@9176:                 ws.clear();
mschaefer@9176:                 qs.clear();
mschaefer@9176:                 km = (double) rows.get(i)[0];
mschaefer@9176:                 ws.add((double) rows.get(i)[2]);
mschaefer@9176:                 qs.add((double) rows.get(i)[3]);
mschaefer@9176:             }
mschaefer@9176:         }
mschaefer@9176:     }
mschaefer@9176: 
mschaefer@9176: 
mschaefer@9176:     /***** METHODS *****/
mschaefer@9176: 
mschaefer@9176:     /**
mschaefer@9176:      * Loads the the W/Q tables of a km range of a river
mschaefer@9176:      *
mschaefer@9176:      * @return The W/Q table finder of the river, or null
mschaefer@9176:      */
mschaefer@9176:     public static WQBaseTableFinder loadValues(final River river, final double fromKm, final double toKm, final Calculation problems) {
mschaefer@9176:         final Session session = SessionHolder.HOLDER.get();
mschaefer@9176:         final SQLQuery colQuery = session.createSQLQuery(SQLCOLUMNS)
mschaefer@9176:                 .addScalar("colindex", StandardBasicTypes.INTEGER)
mschaefer@9176:                 .addScalar("qzone", StandardBasicTypes.STRING);
mschaefer@9176:         colQuery.setParameter("river_id", river.getId());
mschaefer@9176:         final List<Object[]> colnames = colQuery.list();
mschaefer@9176:         if ((colnames == null) || colnames.isEmpty()) {
mschaefer@9176:             problems.addProblem("wq_table.missing");
mschaefer@9176:             return null;
mschaefer@9176:         }
mschaefer@9176:         final SQLQuery wqQuery = session.createSQLQuery(SQLMAIN)
mschaefer@9176:                 .addScalar("station", StandardBasicTypes.DOUBLE)
mschaefer@9176:                 .addScalar("colindex", StandardBasicTypes.INTEGER)
mschaefer@9176:                 .addScalar("w", StandardBasicTypes.DOUBLE)
mschaefer@9176:                 .addScalar("q", StandardBasicTypes.DOUBLE);
mschaefer@9176:         wqQuery.setParameter("river_id", river.getId());
mschaefer@9176:         wqQuery.setParameter("kmfrom", fromKm);
mschaefer@9176:         wqQuery.setParameter("kmto", toKm);
mschaefer@9176:         final List<Object[]> rows = wqQuery.list();
mschaefer@9176:         if ((rows != null) && !rows.isEmpty())
mschaefer@9176:             return new WQBaseTableFinder(river, problems, colnames, rows);
mschaefer@9176:         else {
mschaefer@9176:             problems.addProblem("wq_table.missing");
mschaefer@9176:             return null;
mschaefer@9176:         }
mschaefer@9176:     }
mschaefer@9176: 
mschaefer@9176:     /**
mschaefer@9176:      * If this provider may return valid data at all.
mschaefer@9176:      */
mschaefer@9176:     public boolean isValid() {
mschaefer@9176:         return (this.kmWs != null);
mschaefer@9176:     }
mschaefer@9176: 
mschaefer@9176:     /**
mschaefer@9176:      * Discharge for a W
mschaefer@9176:      *
mschaefer@9176:      * @param station
mschaefer@9176:      *            station to find or interpolate
mschaefer@9176:      * @param w
mschaefer@9176:      *            W in m+NN or m+NHN
mschaefer@9176:      * @return Q, or NegInf for w less than all, or PosInf for w greater then all, or NaN in case of exception
mschaefer@9176:      */
mschaefer@9176:     public double getDischarge(final double station, final double w) {
mschaefer@9176:         if (this.kmWs == null)
mschaefer@9176:             return Double.NaN;
mschaefer@9176:         final Entry<Double, PolynomialSplineFunction> lowerEntry = this.kmWs.floorEntry(Double.valueOf(station));
mschaefer@9176:         final Entry<Double, PolynomialSplineFunction> upperEntry = this.kmWs.ceilingEntry(Double.valueOf(station));
mschaefer@9176:         if ((lowerEntry == null) || (upperEntry == null))
mschaefer@9176:             return Double.NaN;
mschaefer@9176:         try {
mschaefer@9176:             final double lowerQ = lowerEntry.getValue().value(w);
mschaefer@9176:             final double upperQ = upperEntry.getValue().value(w);
mschaefer@9176:             return Linear.linear(station, lowerEntry.getKey().doubleValue(), upperEntry.getKey().doubleValue(), lowerQ, upperQ);
mschaefer@9176:         }
mschaefer@9176:         catch (@SuppressWarnings("unused") final FunctionEvaluationException e) {
mschaefer@9176:             // ignore exception because this can/will happen regularly
mschaefer@9176:             return Double.NaN;
mschaefer@9176:         }
mschaefer@9176:     }
mschaefer@9176: }