teichmann@5863: /* Copyright (C) 2011, 2012, 2013 by Bundesanstalt für Gewässerkunde
teichmann@5863:  * Software engineering by Intevation GmbH
teichmann@5863:  *
teichmann@5994:  * This file is Free Software under the GNU AGPL (>=v3)
teichmann@5863:  * and comes with ABSOLUTELY NO WARRANTY! Check out the
teichmann@5994:  * documentation coming with Dive4Elements River for details.
teichmann@5863:  */
teichmann@5863: 
teichmann@5831: package org.dive4elements.river.artifacts.model.minfo;
raimund@3760: 
raimund@3769: import gnu.trove.TDoubleArrayList;
raimund@3769: 
raimund@3760: import java.util.LinkedList;
raimund@3760: import java.util.List;
raimund@3769: import java.util.Map;
raimund@3760: 
raimund@3760: import org.apache.log4j.Logger;
raimund@3760: 
teichmann@5831: import org.dive4elements.river.artifacts.access.BedQualityAccess;
teichmann@5831: import org.dive4elements.river.artifacts.model.Calculation;
teichmann@5831: import org.dive4elements.river.artifacts.model.CalculationResult;
teichmann@5831: import org.dive4elements.river.artifacts.model.DateRange;
teichmann@5831: import org.dive4elements.river.backend.SedDBSessionHolder;
raimund@3760: 
raimund@3760: 
raimund@3760: public class BedQualityCalculation extends Calculation {
raimund@3760: 
raimund@3760:     private static final Logger logger = Logger
raimund@3760:         .getLogger(BedQualityCalculation.class);
raimund@3760: 
raimund@3760:     protected String river;
raimund@3760:     protected double from;
raimund@3760:     protected double to;
raimund@3769:     protected List<String> bedDiameter;
raimund@3769:     protected List<String> bedloadDiameter;
raimund@3760:     protected List<DateRange> ranges;
raimund@3760: 
raimund@3760:     public BedQualityCalculation() {
raimund@3760:     }
raimund@3760: 
raimund@3760:     public CalculationResult calculate(BedQualityAccess access) {
raimund@3760:         logger.info("BedQualityCalculation.calculate");
raimund@3760: 
felix@7261:         String river = access.getRiverName();
raimund@3760:         Double from = access.getFrom();
raimund@3760:         Double to = access.getTo();
raimund@3769:         List<String> bedDiameter = access.getBedDiameter();
raimund@3769:         List<String> bedloadDiameter = access.getBedloadDiameter();
raimund@3760:         List<DateRange> ranges = access.getDateRanges();
raimund@3760: 
raimund@3760:         if (river == null) {
raimund@3760:             // TODO: i18n
raimund@3760:             addProblem("minfo.missing.river");
raimund@3760:         }
raimund@3760: 
raimund@3760:         if (from == null) {
raimund@3760:             // TODO: i18n
raimund@3760:             addProblem("minfo.missing.from");
raimund@3760:         }
raimund@3760: 
raimund@3760:         if (to == null) {
raimund@3760:             // TODO: i18n
raimund@3760:             addProblem("minfo.missing.to");
raimund@3760:         }
raimund@3760: 
raimund@3760:         if (ranges == null) {
raimund@3760:             // TODO: i18n
raimund@3760:             addProblem("minfo.missing.periods");
raimund@3760:         }
raimund@3760: 
raimund@3760:         if (!hasProblems()) {
raimund@3760:             this.river = river;
raimund@3760:             this.from = from;
raimund@3760:             this.to = to;
raimund@3760:             this.ranges = ranges;
raimund@3769:             this.bedDiameter = bedDiameter;
raimund@3769:             this.bedloadDiameter = bedloadDiameter;
raimund@3760: 
raimund@3760:             SedDBSessionHolder.acquire();
raimund@3760:             try {
raimund@3760:                 return internalCalculate();
raimund@3760:             }
raimund@3760:             finally {
raimund@3760:                 SedDBSessionHolder.release();
raimund@3760:             }
raimund@3760:         }
raimund@3760: 
raimund@3760:         return new CalculationResult();
raimund@3760:     }
raimund@3760: 
raimund@3760:     protected CalculationResult internalCalculate() {
raimund@3760: 
raimund@3760:         List<BedQualityResult> results = new LinkedList<BedQualityResult>();
raimund@3760:         // Calculate for all time periods.
raimund@3760:         for (DateRange dr : ranges) {
raimund@3769:             QualityMeasurements loadMeasurements =
ingo@3785:                 QualityMeasurementFactory.getBedloadMeasurements(
raimund@3769:                     river,
raimund@3769:                     from,
raimund@3769:                     to,
raimund@3769:                     dr.getFrom(),
raimund@3769:                     dr.getTo());
raimund@3769:             QualityMeasurements bedMeasurements =
raimund@3769:                 QualityMeasurementFactory.getBedMeasurements(
raimund@3769:                     river,
raimund@3769:                     from,
raimund@3769:                     to,
raimund@3769:                     dr.getFrom(),
raimund@3769:                     dr.getTo());
ingo@3784:             BedQualityResult result = new BedQualityResult();
ingo@3784:             result.setDateRange(dr);
raimund@3780:             if (bedDiameter != null) {
ingo@3784:                 result.add(calculateBedParameter(bedMeasurements, dr));
raimund@3780:                 for (String bd : bedDiameter) {
ingo@3784:                     BedDiameterResult bedResult =
raimund@3780:                         calculateBed(bedMeasurements, bd, dr);
sascha@3772: 
raimund@3780:                     // Avoid adding empty result sets.
raimund@3780:                     if (!bedResult.isEmpty()) {
ingo@3784:                         result.add(bedResult);
raimund@3780:                     }
raimund@3769:                 }
raimund@3769:             }
raimund@3780:             if (bedloadDiameter != null) {
raimund@3780:                 for (String bld : bedloadDiameter) {
ingo@3784:                     BedloadDiameterResult loadResult =
raimund@3780:                         calculateBedload(loadMeasurements, bld, dr);
ingo@3784:                     result.add(loadResult);
raimund@3780:                 }
raimund@3769:             }
ingo@3784:             results.add(result);
raimund@3760:         }
raimund@3760: 
raimund@3760:         return new CalculationResult(
raimund@3760:             results.toArray(new BedQualityResult[results.size()]), this);
raimund@3760:     }
raimund@3760: 
ingo@3784:     private BedParametersResult calculateBedParameter(
raimund@3769:         QualityMeasurements qm,
ingo@3784:         DateRange dr
raimund@3769:     ) {
raimund@3769:         List<Double> kms = qm.getKms();
ingo@3784:         QualityMeasurements capFiltered = filterCapMeasurements(qm);
ingo@3784:         QualityMeasurements subFiltered = filterSubMeasurements(qm);
raimund@3769:         TDoubleArrayList location = new TDoubleArrayList();
raimund@3769:         TDoubleArrayList porosityCap = new TDoubleArrayList();
raimund@3769:         TDoubleArrayList porositySub = new TDoubleArrayList();
raimund@3769:         TDoubleArrayList densityCap = new TDoubleArrayList();
raimund@3769:         TDoubleArrayList densitySub = new TDoubleArrayList();
ingo@3785: 
ingo@3784:         for(double km : kms) {
ingo@3784:             double[] pCap = calculatePorosity(capFiltered, km);
ingo@3784:             double[] pSub = calculatePorosity(subFiltered, km);
raimund@3769:             double[] dCap = calculateDensity(capFiltered, pCap);
raimund@3769:             double[] dSub = calculateDensity(subFiltered, pSub);
raimund@3769: 
raimund@3769:             double pCapRes = 0d;
raimund@3769:             double pSubRes = 0d;
raimund@3769:             double dCapRes = 0d;
raimund@3769:             double dSubRes = 0d;
raimund@3769:             for (int i = 0; i < pCap.length; i++) {
raimund@3769:                 pCapRes += pCap[i];
raimund@3769:                 dCapRes += dCap[i];
raimund@3769:             }
raimund@3769:             for (int i = 0; i < pSub.length; i++) {
raimund@3769:                 pSubRes += pSub[i];
raimund@3769:                 dSubRes += dSub[i];
raimund@3769:             }
raimund@3769:             location.add(km);
raimund@3769:             porosityCap.add((pCapRes / pCap.length) * 100 );
raimund@3769:             porositySub.add((pSubRes / pSub.length) * 100);
raimund@3769:             densityCap.add((dCapRes / dCap.length) / 1000);
raimund@3769:             densitySub.add((dSubRes / dSub.length) / 1000);
ingo@3785: 
raimund@3769:         }
ingo@3784: 
ingo@3784:         return new BedParametersResult(
raimund@3769:             location,
raimund@3769:             porosityCap,
raimund@3769:             porositySub,
raimund@3769:             densityCap,
raimund@3769:             densitySub);
raimund@3760:     }
raimund@3760: 
ingo@3784:     protected BedDiameterResult calculateBed(
ingo@3784:         QualityMeasurements qm,
ingo@3784:         String diameter,
ingo@3784:         DateRange range
ingo@3784:     ) {
ingo@3784:         List<Double> kms = qm.getKms();
ingo@3784:         TDoubleArrayList location = new TDoubleArrayList();
ingo@3784:         TDoubleArrayList avDiameterCap = new TDoubleArrayList();
ingo@3784:         TDoubleArrayList avDiameterSub = new TDoubleArrayList();
ingo@3784:         for (double km : kms) {
ingo@3784:             //Filter cap and sub measurements.
ingo@3784:             QualityMeasurements capFiltered = filterCapMeasurements(qm);
ingo@3784:             QualityMeasurements subFiltered = filterSubMeasurements(qm);
ingo@3784: 
ingo@3784:             List<QualityMeasurement> cm = capFiltered.getMeasurements(km);
ingo@3784:             List<QualityMeasurement> sm = subFiltered.getMeasurements(km);
ingo@3784: 
ingo@3784:             double avCap = calculateAverage(cm, diameter);
ingo@3784:             double avSub = calculateAverage(sm, diameter);
ingo@3784:             location.add(km);
rrenkert@4642:             avDiameterCap.add(avCap * 1000);// bring to mm.
rrenkert@4642:             avDiameterSub.add(avSub * 1000);
ingo@3784:         }
ingo@3784:         return new BedDiameterResult(
ingo@3784:             diameter,
ingo@3784:             avDiameterCap,
ingo@3784:             avDiameterSub,
ingo@3784:             location);
ingo@3784:     }
ingo@3784: 
raimund@3769:     private double[] calculateDensity(
raimund@3769:         QualityMeasurements capFiltered,
raimund@3769:         double[] porosity
raimund@3769:     ) {
raimund@3769:         double[] density = new double[porosity.length];
raimund@3769:         for (int i = 0; i < porosity.length; i++) {
raimund@3769:             density[i] = (1 - porosity[i]) * 2650;
raimund@3769:         }
raimund@3769:         return density;
raimund@3769:     }
raimund@3769: 
raimund@3769:     private double[] calculatePorosity(
raimund@3769:         QualityMeasurements capFiltered,
ingo@3784:         double km
raimund@3769:     ) {
raimund@3769:         List<QualityMeasurement> list = capFiltered.getMeasurements(km);
raimund@3769:         double[] results = new double[list.size()];
raimund@3769:         int i = 0;
raimund@3769:         for (QualityMeasurement qm : list) {
raimund@3769:             double deviation = calculateDeviation(qm);
raimund@3769:             double p = calculateP(qm);
tom@7670:             double porosity = 0.362 - 0.088 * deviation + 0.219 * p;
raimund@3769:             results[i] = porosity;
raimund@3769:             i++;
raimund@3769:         }
raimund@3769: 
raimund@3769:         return results;
raimund@3769:     }
raimund@3769: 
ingo@3784:     protected BedloadDiameterResult calculateBedload(
raimund@3769:         QualityMeasurements qm,
raimund@3769:         String diameter,
raimund@3769:         DateRange range
raimund@3769:     ) {
raimund@3769:         List<Double> kms = qm.getKms();
raimund@3769:         TDoubleArrayList location = new TDoubleArrayList();
raimund@3769:         TDoubleArrayList avDiameter = new TDoubleArrayList();
raimund@3769:         for (double km : kms) {
raimund@3769:             List<QualityMeasurement> measurements = qm.getMeasurements(km);
raimund@3769:             double mid = calculateAverage(measurements, diameter);
raimund@3769:             location.add(km);
rrenkert@6257:             avDiameter.add(mid * 1000);
raimund@3769:         }
ingo@3784:         return new BedloadDiameterResult(
raimund@3769:             diameter,
raimund@3769:             avDiameter,
raimund@3769:             location,
raimund@3769:             range);
raimund@3769:     }
raimund@3769: 
raimund@3769:     protected double calculateAverage(
raimund@3769:         List<QualityMeasurement> list,
raimund@3769:         String diameter
raimund@3769:     ) {
raimund@3769:         double av = 0;
raimund@3769:         for (QualityMeasurement qm : list) {
raimund@3769:             av += qm.getDiameter(diameter);
raimund@3769:         }
raimund@3769:         return av/list.size();
raimund@3769:     }
raimund@3769: 
raimund@3769:     protected QualityMeasurements filterCapMeasurements(
raimund@3769:         QualityMeasurements qms
raimund@3769:     ) {
raimund@3769:         List<QualityMeasurement> result = new LinkedList<QualityMeasurement>();
raimund@3769:         for (QualityMeasurement qm : qms.getMeasurements()) {
raimund@3769:             if (qm.getDepth1() == 0d && qm.getDepth2() <= 0.3) {
raimund@3769:                 result.add(qm);
raimund@3769:             }
raimund@3769:         }
raimund@3769:         return new QualityMeasurements(result);
raimund@3769:     }
raimund@3769: 
raimund@3769:     protected QualityMeasurements filterSubMeasurements(
raimund@3769:         QualityMeasurements qms
raimund@3769:     ) {
raimund@3769:         List<QualityMeasurement> result = new LinkedList<QualityMeasurement>();
raimund@3769:         for (QualityMeasurement qm : qms.getMeasurements()) {
raimund@3769:             if (qm.getDepth1() > 0d && qm.getDepth2() <= 0.5) {
raimund@3769:                 result.add(qm);
raimund@3769:             }
raimund@3769:         }
raimund@3769:         return new QualityMeasurements(result);
raimund@3769:     }
raimund@3769: 
raimund@3769:     public double calculateDeviation(QualityMeasurement qm) {
raimund@3769:         Map<String, Double> dm = qm.getAllDiameter();
teichmann@6151:         int size = dm.size();
teichmann@6151: 
teichmann@6151:         double phiM    = 0;
teichmann@6151:         double [] phis = new double[size];
teichmann@6151:         double [] ps   = new double[size];
teichmann@6151:         double scale   = -1d/Math.log(2d);
teichmann@6151: 
raimund@3769:         int i = 0;
teichmann@6151:         for (Map.Entry<String, Double> entry: dm.entrySet()) {
teichmann@6151:             double phi = scale*Math.log(entry.getValue());
teichmann@6151:             double p   = calculateWeight(qm, entry.getKey());
teichmann@6151:             phiM += phi * p;
teichmann@6151:             ps[i]   = p;
raimund@3769:             phis[i] = phi;
raimund@3769:             i++;
raimund@3769:         }
raimund@3769: 
raimund@3769:         double sig = 0d;
teichmann@6151:         for (i = 0; i < size; i++) {
raimund@3769:             sig += ps[i] * Math.exp(phis[i] - phiM);
raimund@3769:         }
raimund@3769:         double deviation = Math.sqrt(sig);
raimund@3769:         return deviation;
raimund@3769:     }
raimund@3769: 
raimund@3769:     protected double calculateP(QualityMeasurement qm) {
raimund@3769:         return calculateWeight(qm, "dmin");
raimund@3769:     }
raimund@3769: 
raimund@3769:     public double calculateWeight(QualityMeasurement qm, String diameter) {
raimund@3769:         Map<String, Double> dm = qm.getAllDiameter();
raimund@3769:         double value = qm.getDiameter(diameter);
raimund@3769: 
raimund@3769:         double sum = 0d;
raimund@3769:         for (Double d : dm.values()) {
teichmann@6150:             sum += d.doubleValue();
raimund@3769:         }
raimund@3769:         double weight = sum/100*value;
raimund@3769:         return weight;
raimund@3760:     }
raimund@3760: }
felix@4806: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :