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.fixings;
sascha@3411: 
teichmann@5831: import org.dive4elements.river.artifacts.access.FixAnalysisAccess;
sascha@3411: 
teichmann@5831: import org.dive4elements.river.artifacts.math.fitting.Function;
sascha@3411: 
teichmann@5831: import org.dive4elements.river.artifacts.model.CalculationResult;
teichmann@5831: import org.dive4elements.river.artifacts.model.DateRange;
sascha@3411: 
teichmann@5831: import org.dive4elements.river.artifacts.model.FixingsOverview.AndFilter;
teichmann@5831: import org.dive4elements.river.artifacts.model.FixingsOverview.DateRangeFilter;
teichmann@5831: 
teichmann@5831: import org.dive4elements.river.artifacts.model.FixingsOverview.Fixing.Filter;
teichmann@5831: 
teichmann@5831: import org.dive4elements.river.artifacts.model.FixingsOverview.Fixing;
teichmann@5831: import org.dive4elements.river.artifacts.model.FixingsOverview.IdsFilter;
teichmann@5831: import org.dive4elements.river.artifacts.model.FixingsOverview.KmFilter;
teichmann@5831: import org.dive4elements.river.artifacts.model.FixingsOverview.SectorFilter;
teichmann@5831: 
teichmann@5831: import org.dive4elements.river.artifacts.model.FixingsOverview;
teichmann@5831: import org.dive4elements.river.artifacts.model.Parameters;
teichmann@5831: import org.dive4elements.river.artifacts.model.Range;
teichmann@5831: 
teichmann@5831: import org.dive4elements.river.utils.DateAverager;
teichmann@5831: import org.dive4elements.river.utils.KMIndex;
sascha@3411: 
sascha@3609: import gnu.trove.TIntIntHashMap;
sascha@3609: 
sascha@3411: import java.util.ArrayList;
sascha@3411: import java.util.Date;
sascha@3411: import java.util.List;
sascha@3411: 
sascha@3411: import org.apache.commons.math.stat.descriptive.moment.StandardDeviation;
sascha@3411: 
sascha@3411: import org.apache.log4j.Logger;
sascha@3411: 
sascha@3411: public class FixAnalysisCalculation
sascha@3419: extends      FixCalculation
sascha@3411: {
sascha@3411:     private static Logger log = Logger.getLogger(FixAnalysisCalculation.class);
sascha@3411: 
sascha@3411:     protected DateRange    referencePeriod;
sascha@3411:     protected DateRange [] analysisPeriods;
sascha@3411: 
sascha@3411:     public FixAnalysisCalculation() {
sascha@3411:     }
sascha@3411: 
sascha@3411:     public FixAnalysisCalculation(FixAnalysisAccess access) {
sascha@3419:         super(access);
sascha@3411: 
sascha@3411:         DateRange    referencePeriod = access.getReferencePeriod();
sascha@3411:         DateRange [] analysisPeriods = access.getAnalysisPeriods();
sascha@3411: 
sascha@3411:         if (referencePeriod == null) {
sascha@3411:             addProblem("fix.missing.reference.period");
sascha@3411:         }
sascha@3411: 
sascha@3411:         if (analysisPeriods == null || analysisPeriods.length < 1) {
sascha@3411:             addProblem("fix.missing.analysis.periods");
sascha@3411:         }
sascha@3411: 
sascha@3411:         if (!hasProblems()) {
sascha@3411:             this.referencePeriod = referencePeriod;
sascha@3411:             this.analysisPeriods = analysisPeriods;
sascha@3411:         }
sascha@3411:     }
sascha@3411: 
sascha@3437:     @Override
sascha@3437:     public CalculationResult innerCalculate(
sascha@3437:         FixingsOverview overview,
sascha@3437:         Function        func
sascha@3437:     ) {
sascha@3603:         ColumnCache cc = new ColumnCache();
sascha@3603: 
sascha@3603:         FitResult fitResult = doFitting(overview, cc, func);
sascha@3411: 
sascha@3435:         if (fitResult == null) {
sascha@3411:             return new CalculationResult(this);
sascha@3411:         }
sascha@3411: 
sascha@3435:         KMIndex<AnalysisPeriod []> analysisPeriods =
sascha@3435:             calculateAnalysisPeriods(
sascha@3435:                 func,
sascha@3435:                 fitResult.getParameters(),
sascha@3603:                 overview,
sascha@3603:                 cc);
sascha@3411: 
sascha@3411:         analysisPeriods.sort();
sascha@3411: 
sascha@3435:         FixAnalysisResult far = new FixAnalysisResult(
sascha@3435:             fitResult.getParameters(),
sascha@3435:             fitResult.getReferenced(),
sascha@3435:             fitResult.getOutliers(),
sascha@3411:             analysisPeriods);
sascha@3411: 
teichmann@6875:         // Workaraound to deal with same dates in data set
teichmann@6875:         far.makeReferenceEventsDatesUnique();
teichmann@6877:         far.remapReferenceIndicesToRank();
teichmann@6875: 
teichmann@6877:         far.makeAnalysisEventsUnique();
teichmann@6877:         // TODO: remapping
teichmann@6875: 
sascha@3435:         return new CalculationResult(far, this);
sascha@3411:     }
sascha@3411: 
sascha@3434:     @Override
sascha@3434:     protected Filter createFilter() {
sascha@3434:         Filter ids = super.createFilter();
sascha@3411:         DateRangeFilter rdf = new DateRangeFilter(
sascha@3411:             referencePeriod.getFrom(),
sascha@3411:             referencePeriod.getTo());
sascha@3434:         return new AndFilter().add(rdf).add(ids);
sascha@3411:     }
sascha@3411: 
sascha@3411:     protected KMIndex<AnalysisPeriod []> calculateAnalysisPeriods(
sascha@3411:         Function        function,
sascha@3411:         Parameters      parameters,
sascha@3603:         FixingsOverview overview,
sascha@3603:         ColumnCache cc
sascha@3411:     ) {
sascha@3411:         Range range = new Range(from, to);
sascha@3411: 
sascha@3411:         int kmIndex   = parameters.columnIndex("km");
sascha@3411:         int maxQIndex = parameters.columnIndex("max_q");
sascha@3411: 
sascha@3411:         double [] wq = new double[2];
sascha@3411: 
sascha@3411:         int [] parameterIndices =
sascha@3411:             parameters.columnIndices(function.getParameterNames());
sascha@3411: 
sascha@3411:         double [] parameterValues = new double[parameterIndices.length];
sascha@3411: 
sascha@3411:         DateAverager dateAverager = new DateAverager();
sascha@3411: 
sascha@3411:         KMIndex<AnalysisPeriod []> results =
sascha@3411:             new KMIndex<AnalysisPeriod []>(parameters.size());
sascha@3411: 
sascha@3411:         IdsFilter idsFilter = new IdsFilter(events);
sascha@3411: 
sascha@3609:         TIntIntHashMap [] col2indices =
sascha@3609:             new TIntIntHashMap[analysisPeriods.length];
sascha@3605: 
sascha@3607:         for (int i = 0; i < analysisPeriods.length; ++i) {
sascha@3609:             col2indices[i] = new TIntIntHashMap();
sascha@3607:         }
sascha@3605: 
sascha@3411:         for (int row = 0, R = parameters.size(); row < R; ++row) {
sascha@3411:             double km = parameters.get(row, kmIndex);
sascha@3411:             parameters.get(row, parameterIndices, parameterValues);
sascha@3411: 
sascha@3411:             // This is the paraterized function for a given km.
teichmann@5831:             org.dive4elements.river.artifacts.math.Function instance =
sascha@3411:                 function.instantiate(parameterValues);
sascha@3411: 
sascha@3411:             KmFilter kmFilter = new KmFilter(km);
sascha@3411: 
sascha@3411:             ArrayList<AnalysisPeriod> periodResults =
sascha@3411:                 new ArrayList<AnalysisPeriod>(analysisPeriods.length);
sascha@3411: 
sascha@3607:             for (int ap = 0; ap < analysisPeriods.length; ++ap) {
sascha@3607:                 DateRange analysisPeriod = analysisPeriods[ap];
sascha@3609:                 TIntIntHashMap col2index = col2indices[ap];
sascha@3607: 
sascha@3411:                 DateRangeFilter drf = new DateRangeFilter(
sascha@3411:                     analysisPeriod.getFrom(),
sascha@3411:                     analysisPeriod.getTo());
sascha@3411: 
sascha@3411:                 QWD []    qSectorAverages = new QWD[4];
sascha@3411:                 double [] qSectorStdDevs  = new double[4];
sascha@3411: 
sascha@3411:                 ArrayList<QWD> allQWDs = new ArrayList<QWD>();
sascha@3411: 
sascha@3411:                 // for all Q sectors.
felix@6560:                 for (int qSector = qSectorStart; qSector <= qSectorEnd; ++qSector) {
sascha@3411: 
sascha@3411:                     Filter filter = new AndFilter()
sascha@3411:                         .add(kmFilter)
sascha@3411:                         .add(new SectorFilter(qSector))
sascha@3411:                         .add(drf)
sascha@3411:                         .add(idsFilter);
sascha@3411: 
sascha@3411:                     List<Fixing.Column> metas = overview.filter(range, filter);
sascha@3411: 
sascha@3411:                     if (metas.isEmpty()) {
sascha@3411:                         // No fixings for km and analysis period
sascha@3411:                         continue;
sascha@3411:                     }
sascha@3411: 
sascha@3411:                     double sumQ = 0.0;
sascha@3411:                     double sumW = 0.0;
sascha@3411: 
sascha@3411:                     StandardDeviation stdDev = new StandardDeviation();
sascha@3411: 
sascha@3411:                     List<QWD> qwds = new ArrayList<QWD>(metas.size());
sascha@3411: 
sascha@3411:                     dateAverager.clear();
sascha@3411: 
sascha@3411:                     for (Fixing.Column meta: metas) {
sascha@3411:                         if (meta.findQSector(km) != qSector) {
sascha@3411:                             // Ignore not matching sectors.
sascha@3411:                             continue;
sascha@3411:                         }
sascha@3411: 
sascha@3411:                         Column column = cc.getColumn(meta);
sascha@3411:                         if (column == null || !column.getQW(km, wq)) {
sascha@3411:                             continue;
sascha@3411:                         }
sascha@3411: 
sascha@3411:                         double fw = instance.value(wq[1]);
sascha@3411:                         if (Double.isNaN(fw)) {
sascha@3411:                             continue;
sascha@3411:                         }
sascha@3411: 
sascha@3411:                         double dw = (wq[0] - fw)*100.0;
sascha@3411: 
sascha@3411:                         stdDev.increment(dw);
sascha@3411: 
sascha@3411:                         Date date = column.getDate();
sascha@3411:                         String description = column.getDescription();
sascha@3411: 
sascha@3411:                         QWD qwd = new QWD(
sascha@3604:                             wq[1], wq[0],
sascha@3604:                             description,
sascha@3604:                             date, true,
sascha@3609:                             dw, getIndex(col2index, column.getIndex()));
sascha@3411: 
sascha@3411:                         qwds.add(qwd);
sascha@3411: 
sascha@3411:                         sumW += wq[0];
sascha@3411:                         sumQ += wq[1];
sascha@3411: 
sascha@3411:                         dateAverager.add(date);
sascha@3411:                     }
sascha@3411: 
sascha@3411:                     // Calulate average per Q sector.
sascha@3411:                     int N = qwds.size();
sascha@3411:                     if (N > 0) {
sascha@3411:                         allQWDs.addAll(qwds);
sascha@3411:                         double avgW = sumW / N;
sascha@3411:                         double avgQ = sumQ / N;
sascha@3411: 
sascha@3411:                         double avgFw = instance.value(avgQ);
sascha@3411:                         if (!Double.isNaN(avgFw)) {
sascha@3411:                             double avgDw = (avgW - avgFw)*100.0;
sascha@3411:                             Date avgDate = dateAverager.getAverage();
sascha@3411: 
sascha@3411:                             String avgDescription = "avg.deltawt." + qSector;
sascha@3411: 
sascha@3411:                             QWD avgQWD = new QWD(
sascha@3604:                                 avgQ, avgW, avgDescription, avgDate, true, avgDw, 0);
sascha@3411: 
sascha@3411:                             qSectorAverages[qSector] = avgQWD;
sascha@3411:                         }
sascha@3411:                         qSectorStdDevs[qSector] = stdDev.getResult();
sascha@3411:                     }
sascha@3411:                     else {
sascha@3411:                         qSectorStdDevs[qSector] = Double.NaN;
sascha@3411:                     }
sascha@3411:                 } // for all Q sectors
sascha@3411: 
sascha@3411:                 QWD [] aqwds = allQWDs.toArray(new QWD[allQWDs.size()]);
sascha@3411: 
sascha@3411:                 AnalysisPeriod periodResult = new AnalysisPeriod(
sascha@3411:                     analysisPeriod,
sascha@3411:                     aqwds,
sascha@3411:                     qSectorAverages,
sascha@3411:                     qSectorStdDevs);
sascha@3411:                 periodResults.add(periodResult);
sascha@3411:             }
sascha@3411: 
sascha@3411:             double maxQ = -Double.MAX_VALUE;
sascha@3411:             for (AnalysisPeriod ap: periodResults) {
sascha@3411:                 double q = ap.getMaxQ();
sascha@3411:                 if (q > maxQ) {
sascha@3411:                     maxQ = q;
sascha@3411:                 }
sascha@3411:             }
sascha@3411: 
sascha@3411:             double oldMaxQ = parameters.get(row, maxQIndex);
sascha@3411:             if (oldMaxQ < maxQ) {
sascha@3411:                 parameters.set(row, maxQIndex, maxQ);
sascha@3411:             }
sascha@3411: 
sascha@3411:             results.add(km, periodResults.toArray(
sascha@3411:                 new AnalysisPeriod[periodResults.size()]));
sascha@3411:         }
sascha@3411: 
sascha@3411:         return results;
sascha@3411:     }
sascha@3609: 
felix@6574:     /** Returns the mapped value of colIdx or the size of the hashmap. */
sascha@3609:     private static final int getIndex(TIntIntHashMap map, int colIdx) {
sascha@3609:         if (map.containsKey(colIdx)) {
sascha@3609:             return map.get(colIdx);
sascha@3609:         }
sascha@3609:         int index = map.size();
sascha@3609:         map.put(colIdx, index);
sascha@3609:         return index;
sascha@3609:     }
sascha@3411: }
sascha@3411: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :