mschaefer@9295: /** Copyright (C) 2017 by Bundesanstalt für Gewässerkunde
mschaefer@9295:  * Software engineering by
mschaefer@9295:  *  Björnsen Beratende Ingenieure GmbH
mschaefer@9295:  *  Dr. Schumacher Ingenieurbüro für Wasser und Umwelt
mschaefer@9295:  *
mschaefer@9295:  * This file is Free Software under the GNU AGPL (>=v3)
mschaefer@9295:  * and comes with ABSOLUTELY NO WARRANTY! Check out the
mschaefer@9295:  * documentation coming with Dive4Elements River for details.
mschaefer@9295:  */
mschaefer@9295: package org.dive4elements.river.artifacts.uinfo.salix;
mschaefer@9295: 
mschaefer@9375: import java.math.BigDecimal;
mschaefer@9295: import java.util.ArrayList;
mschaefer@9397: import java.util.Collection;
mschaefer@9295: import java.util.HashMap;
mschaefer@9295: import java.util.List;
mschaefer@9295: import java.util.Map;
mschaefer@9316: import java.util.Map.Entry;
mschaefer@9309: import java.util.NavigableMap;
mschaefer@9295: 
mschaefer@9295: import org.dive4elements.river.artifacts.WINFOArtifact;
mschaefer@9295: import org.dive4elements.river.artifacts.access.ComputationRangeAccess;
mschaefer@9397: import org.dive4elements.river.artifacts.common.AbstractResultType;
mschaefer@9295: import org.dive4elements.river.artifacts.common.GeneralResultType;
mschaefer@9295: import org.dive4elements.river.artifacts.common.ResultRow;
mschaefer@9295: import org.dive4elements.river.artifacts.model.Calculation;
mschaefer@9295: import org.dive4elements.river.artifacts.model.WstValueTable;
mschaefer@9295: import org.dive4elements.river.artifacts.model.WstValueTable.QPosition;
mschaefer@9295: import org.dive4elements.river.artifacts.model.WstValueTableFactory;
mschaefer@9295: import org.dive4elements.river.artifacts.sinfo.common.RiverInfoProvider;
mschaefer@9295: import org.dive4elements.river.artifacts.sinfo.tkhstate.WinfoArtifactWrapper;
mschaefer@9295: import org.dive4elements.river.artifacts.uinfo.UINFOArtifact;
gernotbelger@9328: import org.dive4elements.river.artifacts.uinfo.common.UInfoResultType;
mschaefer@9309: import org.dive4elements.river.artifacts.uinfo.salix.SalixLineAccess.ScenarioType;
mschaefer@9295: import org.dive4elements.river.model.Gauge;
mschaefer@9295: import org.dive4elements.river.model.MainValue;
mschaefer@9295: import org.dive4elements.river.model.MainValueType.MainValueTypeKey;
mschaefer@9361: import org.dive4elements.river.utils.Formatter;
mschaefer@9295: 
mschaefer@9295: /**
mschaefer@9295:  * Calculation of the result rows of the u-info salix line calc mode
mschaefer@9295:  *
mschaefer@9295:  * @author Matthias Schäfer
mschaefer@9295:  */
gernotbelger@9321: final class SalixLineCalculator {
mschaefer@9295: 
mschaefer@9375:     private static final BigDecimal SALIX_DISTANCE = new BigDecimal("2.31");
mschaefer@9295:     private final List<ResultRow> rows = new ArrayList<>();
mschaefer@9295: 
mschaefer@9295:     private final RiverInfoProvider riverInfoProvider;
mschaefer@9295: 
mschaefer@9295:     private final Map<Gauge, QPosition> gaugeMwPos;
mschaefer@9295:     private final Map<Gauge, QPosition> gaugeMnwPos;
mschaefer@9295:     private final Map<Gauge, QPosition> gaugeMhwPos;
mschaefer@9430:     private final Map<Gauge, QPosition> gaugeHw5Pos;
mschaefer@9368:     private QPosition refGaugeMwPos;
mschaefer@9368:     private QPosition refGaugeMnwPos;
mschaefer@9368:     private QPosition refGaugeMhwPos;
mschaefer@9430:     private QPosition refGaugeHw5Pos;
mschaefer@9397:     private Gauge firstGauge;
mschaefer@9295: 
mschaefer@9295:     private Calculation problems;
mschaefer@9295: 
mschaefer@9295:     private WstValueTable wst;
mschaefer@9295: 
gernotbelger@9321:     public SalixLineCalculator(final RiverInfoProvider riverInfoProvider) {
mschaefer@9295:         this.riverInfoProvider = riverInfoProvider;
mschaefer@9295:         this.gaugeMwPos = new HashMap<>();
mschaefer@9295:         this.gaugeMnwPos = new HashMap<>();
mschaefer@9295:         this.gaugeMhwPos = new HashMap<>();
mschaefer@9430:         this.gaugeHw5Pos = new HashMap<>();
mschaefer@9295:     }
mschaefer@9295: 
mschaefer@9295:     /**
mschaefer@9295:      * Calculate the salix line result rows
mschaefer@9295:      */
mschaefer@9309:     public void execute(final Calculation problems, final UINFOArtifact uinfo, final NavigableMap<Double, List<Double>> rangeScenarios,
mschaefer@9361:             final ScenarioType scenarioType, final String[] scenarioLabels, final String rangeString, final String additionalString,
mschaefer@9361:             final SalixLineCalculationResults results) {
mschaefer@9295: 
mschaefer@9295:         this.problems = problems;
mschaefer@9295:         this.wst = WstValueTableFactory.getTable(this.riverInfoProvider.getRiver());
mschaefer@9295: 
mschaefer@9295:         fetchGaugeMainValuePositions();
mschaefer@9295: 
mschaefer@9295:         final WINFOArtifact winfo = new WinfoArtifactWrapper(uinfo);
mschaefer@9295:         winfo.addStringData("ld_mode", "distance");
mschaefer@9295:         winfo.addStringData("ld_step", "100");
mschaefer@9295:         for (final double station : new ComputationRangeAccess(winfo).getKms()) {
gernotbelger@9321:             this.rows.add(createRow(station, rangeScenarios));
mschaefer@9295:         }
mschaefer@9309:         if (scenarioType == ScenarioType.REGIONAL)
mschaefer@9361:             results.addResult(new SalixLineCalculationRegionalResult("Salix-regional", scenarioLabels, rangeString, additionalString, this.rows), problems);
mschaefer@9309:         else if (scenarioType == ScenarioType.SUPRAREGIONAL)
mschaefer@9361:             results.addResult(new SalixLineCalculationSupraRegionalResult("Salix-supra", scenarioLabels, rangeString, additionalString, this.rows), problems);
mschaefer@9309:         else if (scenarioType == ScenarioType.HISTORICAL)
mschaefer@9361:             results.addResult(new SalixLineCalculationHistoricalResult("Salix-hist", scenarioLabels, rangeString, additionalString, this.rows), problems);
mschaefer@9309:         else
mschaefer@9361:             results.addResult(new SalixLineCalculationResult("Salix-simple", this.rows), problems);
mschaefer@9295:     }
mschaefer@9295: 
mschaefer@9295:     /**
mschaefer@9361:      * Fetch MQ, MNQ and MHQ of all gauges and determine the wst QPosition for each one
mschaefer@9295:      */
mschaefer@9295:     private void fetchGaugeMainValuePositions() {
mschaefer@9295:         this.gaugeMwPos.clear();
mschaefer@9295:         this.gaugeMnwPos.clear();
mschaefer@9295:         this.gaugeMhwPos.clear();
mschaefer@9430:         this.gaugeHw5Pos.clear();
mschaefer@9430: 
mschaefer@9397:         this.firstGauge = null;
mschaefer@9295:         for (final Gauge gauge : this.riverInfoProvider.getGauges()) {
mschaefer@9295:             this.gaugeMwPos.put(gauge, null);
mschaefer@9295:             this.gaugeMnwPos.put(gauge, null);
mschaefer@9295:             this.gaugeMhwPos.put(gauge, null);
mschaefer@9430:             this.gaugeHw5Pos.put(gauge, null);
mschaefer@9295:             final double gaugeKm = gauge.getStation().doubleValue();
mschaefer@9361:             for (final MainValue mv : MainValue.getValuesOfGaugeAndType(gauge, MainValueTypeKey.Q)) {
mschaefer@9361:                 if (mv.getMainValue().getName().equalsIgnoreCase("mq"))
mschaefer@9361:                     this.gaugeMwPos.put(gauge, this.wst.getQPosition(gaugeKm, mv.getValue().doubleValue()));
mschaefer@9361:                 else if (mv.getMainValue().getName().equalsIgnoreCase("mnq"))
mschaefer@9361:                     this.gaugeMnwPos.put(gauge, this.wst.getQPosition(gaugeKm, mv.getValue().doubleValue()));
mschaefer@9361:                 else if (mv.getMainValue().getName().equalsIgnoreCase("mhq"))
mschaefer@9361:                     this.gaugeMhwPos.put(gauge, this.wst.getQPosition(gaugeKm, mv.getValue().doubleValue()));
mschaefer@9430:                 else if (mv.getMainValue().getName().equalsIgnoreCase("hq5"))
mschaefer@9430:                     this.gaugeHw5Pos.put(gauge, this.wst.getQPosition(gaugeKm, mv.getValue().doubleValue()));
mschaefer@9295:             }
mschaefer@9397:             if (this.firstGauge == null) {
mschaefer@9368:                 this.refGaugeMwPos = this.gaugeMwPos.get(gauge);
mschaefer@9368:                 this.refGaugeMnwPos = this.gaugeMnwPos.get(gauge);
mschaefer@9368:                 this.refGaugeMhwPos = this.gaugeMhwPos.get(gauge);
mschaefer@9430:                 this.refGaugeHw5Pos = this.gaugeHw5Pos.get(gauge);
mschaefer@9397:                 this.firstGauge = gauge;
mschaefer@9368:             }
mschaefer@9368:         }
mschaefer@9368:         if (this.refGaugeMwPos == null)
mschaefer@9368:             this.problems.addProblem("uinfo_salix_calc.warning.missing_mq");
mschaefer@9368:         else {
mschaefer@9368:             if (this.refGaugeMhwPos == null)
mschaefer@9368:                 this.problems.addProblem("uinfo_salix_calc.warning.missing_mhq");
mschaefer@9368:             if (this.refGaugeMnwPos == null)
mschaefer@9368:                 this.problems.addProblem("uinfo_salix_calc.warning.missing_mnq");
mschaefer@9295:         }
mschaefer@9295:     }
mschaefer@9295: 
mschaefer@9295:     /**
mschaefer@9295:      * Create a result row for a station and its gauge, and add w-q-values as selected
mschaefer@9295:      */
gernotbelger@9321:     private ResultRow createRow(final double station, final NavigableMap<Double, List<Double>> rangeScenarios) {
mschaefer@9295: 
mschaefer@9295:         final ResultRow row = ResultRow.create();
mschaefer@9295:         row.putValue(GeneralResultType.station, station);
mschaefer@9394:         // Find station's gauge (obsolete version which calculates gauge-wise)
mschaefer@9368:         // final Gauge gauge = this.riverInfoProvider.getGauge(station, true);
mschaefer@9316:         // Interpolate mnw, mw, and mhw
mschaefer@9368:         // final double mnw = interpolateW(station, this.gaugeMnwPos.get(gauge));
mschaefer@9368:         // final double mw = interpolateW(station, this.gaugeMwPos.get(gauge));
mschaefer@9368:         // final double mhw = interpolateW(station, this.gaugeMhwPos.get(gauge));
mschaefer@9368:         final double mnw = interpolateW(station, this.refGaugeMnwPos);
mschaefer@9368:         final double mw = interpolateW(station, this.refGaugeMwPos);
mschaefer@9368:         final double mhw = interpolateW(station, this.refGaugeMhwPos);
mschaefer@9430:         final double hw5 = interpolateW(station, this.refGaugeHw5Pos);
gernotbelger@9429:         row.putValue(UInfoResultType.waterlevelMNW, mnw);
gernotbelger@9429:         row.putValue(UInfoResultType.waterlevelMW, mw);
gernotbelger@9429:         row.putValue(UInfoResultType.waterlevelMHW, mhw);
mschaefer@9430:         row.putValue(UInfoResultType.waterlevelMH5, hw5);
gernotbelger@9429: 
mschaefer@9316:         // Calc salix-line and mw-mnw
mschaefer@9295:         row.putValue(UInfoResultType.salixline, calcSalix(mhw, mw));
mschaefer@9361:         row.putValue(UInfoResultType.salix_mw_mnw, calcMwmnw(mw, mnw));
mschaefer@9430:         row.putValue(UInfoResultType.salixw, mhw - SALIX_DISTANCE.doubleValue());
mschaefer@9316:         // Calc scenario values (always all scenario types set, Result variant extracts the fields needed)
mschaefer@9316:         final List<SalixScenario> scenarios = new ArrayList<>();
mschaefer@9394:         final List<Double> deltaws = getDeltaWs(station, rangeScenarios);
mschaefer@9394:         for (final Double deltaw : deltaws) {
mschaefer@9394:             if (deltaw != null) {
mschaefer@9394:                 final double salix = calcSalix(mhw, mw + deltaw);
mschaefer@9430:                 scenarios.add(new SalixScenario((int) (deltaw * 100), salix, mhw + deltaw - SALIX_DISTANCE.doubleValue()));
mschaefer@9316:             }
mschaefer@9361:             else {
mschaefer@9361:                 scenarios.add(null);
mschaefer@9361:             }
mschaefer@9316:         }
mschaefer@9361:         row.putValue(UInfoResultType.customMultiRowColSalixScenarios, scenarios);
mschaefer@9397:         row.putValue(GeneralResultType.gaugeLabel, this.riverInfoProvider.findGauge(station));
mschaefer@9295:         return row;
mschaefer@9295:     }
mschaefer@9295: 
mschaefer@9295:     /**
mschaefer@9295:      * Interpolates the W for a station with a fixed (virtual) wst column position
mschaefer@9295:      */
mschaefer@9295:     private double interpolateW(final double station, final QPosition qPosition) {
mschaefer@9295:         if (qPosition != null)
mschaefer@9295:             return this.wst.interpolateW(station, qPosition, this.problems);
gernotbelger@9321:         return Double.NaN;
mschaefer@9295:     }
mschaefer@9295: 
mschaefer@9295:     /**
mschaefer@9295:      * Calculates the salix value
mschaefer@9295:      */
mschaefer@9295:     private double calcSalix(final double mhw, final double mw) {
mschaefer@9382:         if (Double.isNaN(mw) || Double.isInfinite(mw) || Double.isNaN(mhw) || Double.isInfinite(mhw))
mschaefer@9382:             return mw; // preserving NaN or Infinity
mschaefer@9375:         return Formatter.roundW(mhw).subtract(SALIX_DISTANCE).subtract(Formatter.roundW(mw)).doubleValue();
mschaefer@9295:     }
mschaefer@9295: 
mschaefer@9295:     /**
mschaefer@9295:      * Calculates the inverse MW-MNW difference
mschaefer@9295:      */
mschaefer@9295:     private double calcMwmnw(final double mw, final double mnw) {
mschaefer@9382:         if (Double.isNaN(mw) || Double.isInfinite(mw) || Double.isNaN(mnw) || Double.isInfinite(mnw))
mschaefer@9382:             return mnw - mw; // preserving NaN or Inifinity
mschaefer@9375:         return Formatter.roundW(mnw).subtract(Formatter.roundW(mw)).doubleValue();
mschaefer@9295:     }
mschaefer@9316: 
mschaefer@9316:     /**
mschaefer@9394:      * Gets the station-specific list of delta-ws of the active scenario, at least with one null item in any case
mschaefer@9316:      */
mschaefer@9394:     private List<Double> getDeltaWs(final double station, final NavigableMap<Double, List<Double>> rangeScenarios) {
gernotbelger@9321:         final Entry<Double, List<Double>> stationScenarios = rangeScenarios.floorEntry(station);
mschaefer@9394:         if (stationScenarios != null) {
mschaefer@9394:             return stationScenarios.getValue();
mschaefer@9394:         }
mschaefer@9394:         final List<Double> noScen = new ArrayList<>();
mschaefer@9394:         noScen.add(null);
mschaefer@9394:         return noScen;
mschaefer@9316:     }
mschaefer@9397: 
mschaefer@9397:     /**
mschaefer@9397:      * Find and return a height (iota, w main value) of a station in a previously calculated result
mschaefer@9397:      */
mschaefer@9397:     public double fetchStationHeight(final Calculation problems, final double station, final AbstractResultType resultType,
mschaefer@9397:             final SalixLineCalculationResult result) {
mschaefer@9397: 
mschaefer@9397:         // Search the station in the previously calculated result rows
mschaefer@9397:         final ResultRow stationRow = searchStation(station, result.getRows());
mschaefer@9397:         if (stationRow != null)
mschaefer@9397:             return stationRow.getDoubleValue(resultType);
mschaefer@9397:         return Double.NaN;
mschaefer@9397:     }
mschaefer@9397: 
mschaefer@9397:     /**
mschaefer@9397:      * Searches the row of a station in a result rows collection
mschaefer@9397:      */
mschaefer@9397:     private ResultRow searchStation(final double station, final Collection<ResultRow> rows) {
mschaefer@9397:         for (final ResultRow row : rows) {
mschaefer@9397:             if (row.getDoubleValue(GeneralResultType.station) > station + 0.0001)
mschaefer@9397:                 return row;
mschaefer@9397:         }
mschaefer@9397:         return null;
mschaefer@9397:     }
gernotbelger@9321: }