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@2729: 
teichmann@5831: import org.dive4elements.artifacts.common.utils.StringUtils;
sascha@3434: 
teichmann@5831: import org.dive4elements.river.artifacts.access.FixAccess;
sascha@3434: 
teichmann@5831: import org.dive4elements.river.artifacts.math.fitting.Function;
teichmann@5831: import org.dive4elements.river.artifacts.math.fitting.FunctionFactory;
sascha@2729: 
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.FixingsColumn;
teichmann@5831: import org.dive4elements.river.artifacts.model.FixingsColumnFactory;
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: 
teichmann@5831: import org.dive4elements.river.artifacts.model.FixingsOverview;
teichmann@5831: import org.dive4elements.river.artifacts.model.FixingsOverviewFactory;
teichmann@5831: import org.dive4elements.river.artifacts.model.Parameters;
teichmann@5831: 
teichmann@5831: import org.dive4elements.river.utils.DoubleUtil;
teichmann@5831: import org.dive4elements.river.utils.KMIndex;
sascha@2729: 
sascha@2729: import java.util.ArrayList;
sascha@2744: import java.util.Date;
sascha@3434: import java.util.HashMap;
sascha@3434: import java.util.List;
sascha@3434: import java.util.Map;
sascha@2729: 
sascha@2729: import org.apache.log4j.Logger;
sascha@2729: 
felix@5150: /** Calculation base class for fix. */
sascha@3437: public abstract class FixCalculation
sascha@3437: extends               Calculation
sascha@2729: {
sascha@2786:     private static Logger log = Logger.getLogger(FixCalculation.class);
sascha@2729: 
sascha@3435:     public static final double EPSILON = 1e-4;
sascha@3435: 
sascha@3435:     public static final String [] STANDARD_COLUMNS = {
sascha@3435:         "km", "chi_sqr", "max_q", "std-dev"
sascha@3435:     };
sascha@3435: 
sascha@3435:     protected static class FitResult {
sascha@3435: 
sascha@3435:         protected Parameters      parameters;
sascha@3435:         protected KMIndex<QWD []> referenced;
sascha@3729:         protected KMIndex<QWI []> outliers;
sascha@3435: 
sascha@3435:         public FitResult() {
sascha@3435:         }
sascha@3435: 
sascha@3435:         public FitResult(
sascha@3435:             Parameters      parameters,
sascha@3435:             KMIndex<QWD []> referenced,
sascha@3729:             KMIndex<QWI []> outliers
sascha@3435:         ) {
sascha@3435:             this.parameters = parameters;
sascha@3435:             this.referenced = referenced;
sascha@3435:             this.outliers   = outliers;
sascha@3435:         }
sascha@3435: 
sascha@3435:         public Parameters getParameters() {
sascha@3435:             return parameters;
sascha@3435:         }
sascha@3435: 
sascha@3435:         public KMIndex<QWD []> getReferenced() {
sascha@3435:             return referenced;
sascha@3435:         }
sascha@3435: 
sascha@3729:         public KMIndex<QWI []> getOutliers() {
sascha@3435:             return outliers;
sascha@3435:         }
sascha@3435:     } // class FitResult
sascha@3435: 
sascha@3434:     /** Helper class to bundle the meta information of a column
sascha@3434:      *  and the real data.
sascha@3434:      */
sascha@3434:     protected static class Column {
sascha@3434: 
sascha@3434:         protected Fixing.Column meta;
sascha@3434:         protected FixingsColumn data;
sascha@3603:         protected int           index;
sascha@3434: 
sascha@3434:         public Column() {
sascha@3434:         }
sascha@3434: 
sascha@3603:         public Column(Fixing.Column meta, FixingsColumn data, int index) {
sascha@3603:             this.meta  = meta;
sascha@3603:             this.data  = data;
sascha@3603:             this.index = index;
sascha@3434:         }
sascha@3434: 
sascha@3434:         public Date getDate() {
sascha@3434:             return meta.getStartTime();
sascha@3434:         }
sascha@3434: 
sascha@3434:         public String getDescription() {
sascha@3434:             return meta.getDescription();
sascha@3434:         }
sascha@3434: 
sascha@3603:         public int getIndex() {
sascha@3603:             return index;
sascha@3603:         }
sascha@3603: 
teichmann@6875:         public int getId() {
teichmann@6875:             return meta.getId();
teichmann@6875:         }
teichmann@6875: 
sascha@3434:         public boolean getQW(
sascha@3434:             double    km,
sascha@3434:             double [] qs,
sascha@3434:             double [] ws,
sascha@3434:             int       index
sascha@3434:         ) {
sascha@3434:             qs[index] = data.getQ(km);
sascha@3434:             return data.getW(km, ws, index);
sascha@3434:         }
sascha@3434: 
sascha@3434:         public boolean getQW(double km, double [] wq) {
sascha@3434:             data.getW(km, wq, 0);
sascha@3434:             if (Double.isNaN(wq[0])) return false;
sascha@3434:             wq[1] = data.getQ(km);
sascha@3434:             return !Double.isNaN(wq[1]);
sascha@3434:         }
sascha@3434:     } // class Column
sascha@3434: 
sascha@3434:     /**
sascha@3434:      * Helper class to find the data belonging to meta info more quickly.
sascha@3434:      */
sascha@3434:     protected static class ColumnCache {
sascha@3434: 
sascha@3434:         protected Map<Integer, Column> columns;
sascha@3434: 
sascha@3434:         public ColumnCache() {
sascha@3434:             columns = new HashMap<Integer, Column>();
sascha@3434:         }
sascha@3434: 
sascha@3434:         public Column getColumn(Fixing.Column meta) {
sascha@3434:             Integer key = meta.getId();
sascha@3434:             Column column = columns.get(key);
sascha@3434:             if (column == null) {
sascha@3434:                 FixingsColumn data = FixingsColumnFactory
sascha@3434:                     .getInstance()
sascha@3434:                     .getColumnData(meta);
sascha@3434:                 if (data != null) {
sascha@3603:                     column = new Column(meta, data, columns.size());
sascha@3434:                     columns.put(key, column);
sascha@3434:                 }
sascha@3434:             }
sascha@3434:             return column;
sascha@3434:         }
sascha@3434:     } // class ColumnCache
sascha@3434: 
sascha@3434: 
sascha@3419:     protected String  river;
sascha@3419:     protected double  from;
sascha@3419:     protected double  to;
sascha@3419:     protected double  step;
sascha@3419:     protected boolean preprocessing;
sascha@3419:     protected String  function;
sascha@3419:     protected int []  events;
sascha@3419:     protected int     qSectorStart;
sascha@3419:     protected int     qSectorEnd;
sascha@2729: 
sascha@2729:     public FixCalculation() {
sascha@2729:     }
sascha@2729: 
sascha@3419:     public FixCalculation(FixAccess access) {
felix@7261:         String  river         = access.getRiverName();
sascha@3451:         Double  from          = access.getFrom();
sascha@3451:         Double  to            = access.getTo();
sascha@3451:         Double  step          = access.getStep();
sascha@3451:         String  function      = access.getFunction();
sascha@3451:         int []  events        = access.getEvents();
sascha@3451:         Integer qSectorStart  = access.getQSectorStart();
sascha@3451:         Integer qSectorEnd    = access.getQSectorEnd();
sascha@3451:         Boolean preprocessing = access.getPreprocessing();
sascha@2729: 
sascha@2729:         if (river == null) {
sascha@2729:             addProblem("fix.missing.river");
sascha@2729:         }
sascha@2729: 
sascha@2729:         if (from == null) {
sascha@2729:             addProblem("fix.missing.from");
sascha@2729:         }
sascha@2729: 
sascha@2729:         if (to == null) {
sascha@2729:             addProblem("fix.missing.to");
sascha@2729:         }
sascha@2729: 
sascha@2729:         if (step == null) {
sascha@2729:             addProblem("fix.missing.step");
sascha@2729:         }
sascha@2729: 
sascha@2729:         if (function == null) {
sascha@2729:             addProblem("fix.missing.function");
sascha@2729:         }
sascha@2729: 
sascha@2729:         if (events == null || events.length < 1) {
sascha@2729:             addProblem("fix.missing.events");
sascha@2729:         }
sascha@2729: 
sascha@2744:         if (qSectorStart == null) {
sascha@2744:             addProblem("fix.missing.qstart.sector");
sascha@2744:         }
sascha@2744: 
sascha@2744:         if (qSectorEnd == null) {
sascha@2744:             addProblem("fix.missing.qend.sector");
sascha@2744:         }
sascha@2744: 
sascha@3419:         if (preprocessing == null) {
sascha@3419:             addProblem("fix.missing.preprocessing");
sascha@3419:         }
sascha@3419: 
sascha@2729:         if (!hasProblems()) {
sascha@3419:             this.river         = river;
sascha@3419:             this.from          = from;
sascha@3419:             this.to            = to;
sascha@3419:             this.step          = step;
sascha@3419:             this.function      = function;
sascha@3419:             this.events        = events;
sascha@3419:             this.qSectorStart  = qSectorStart;
sascha@3419:             this.qSectorEnd    = qSectorEnd;
sascha@3419:             this.preprocessing = preprocessing;
sascha@2729:         }
sascha@2729:     }
sascha@2729: 
sascha@3434:     protected static String toString(
sascha@3434:         String [] parameterNames,
sascha@3434:         double [] values
sascha@3434:     ) {
sascha@3434:         StringBuilder sb = new StringBuilder();
sascha@3434:         for (int i = 0; i < parameterNames.length; ++i) {
sascha@3434:             if (i > 0) sb.append(", ");
sascha@3434:             sb.append(parameterNames[i]).append(": ").append(values[i]);
sascha@3434:         }
sascha@3434:         return sb.toString();
sascha@3434:     }
sascha@2729: 
felix@5733: 
teichmann@5994:     /** Create filter to accept only the chosen events.
teichmann@5994:      *  This factored out out to be overwritten.
teichmann@5741:      */
teichmann@5741:     protected Filter createFilter() {
sascha@3434:         return new IdsFilter(events);
sascha@3434:     }
sascha@3434: 
sascha@3603:     protected List<Column> getEventColumns(
sascha@3603:         FixingsOverview overview,
sascha@3603:         ColumnCache     cc
sascha@3603:     ) {
sascha@3434:         FixingsColumnFactory fcf = FixingsColumnFactory.getInstance();
sascha@3434: 
teichmann@5741:         Filter filter = createFilter();
sascha@3434: 
sascha@3434:         List<Fixing.Column> metas = overview.filter(null, filter);
sascha@3434: 
sascha@3434:         List<Column> columns = new ArrayList<Column>(metas.size());
sascha@3434: 
sascha@3434:         for (Fixing.Column meta: metas) {
sascha@3434: 
sascha@3603:             Column data = cc.getColumn(meta);
sascha@3434:             if (data == null) {
sascha@3434:                 addProblem("fix.cannot.load.data");
sascha@3434:             }
sascha@3434:             else {
sascha@3603:                 columns.add(data);
sascha@3434:             }
sascha@3434:         }
sascha@3434: 
sascha@3434:         return columns;
sascha@3419:     }
sascha@3435: 
felix@5733:     // Fit a function to the given points from fixation.
sascha@3603:     protected FitResult doFitting(
sascha@3603:         FixingsOverview overview,
sascha@3603:         ColumnCache     cc,
sascha@3603:         Function        func
sascha@3603:     ) {
ingo@2792:         boolean debug = log.isDebugEnabled();
ingo@2792: 
sascha@3603:         final List<Column> eventColumns = getEventColumns(overview, cc);
sascha@3435: 
sascha@3435:         if (eventColumns.size() < 2) {
sascha@3435:             addProblem("fix.too.less.data.columns");
sascha@3435:             return null;
sascha@3435:         }
sascha@3435: 
sascha@3435:         final double  [] qs = new double[eventColumns.size()];
sascha@3435:         final double  [] ws = new double[qs.length];
sascha@3435:         final boolean [] interpolated = new boolean[ws.length];
sascha@3435: 
sascha@3435:         Fitting.QWDFactory qwdFactory = new Fitting.QWDFactory() {
sascha@3435:             @Override
sascha@3435:             public QWD create(double q, double w) {
sascha@3435:                 // Check all the event columns for close match
sascha@3435:                 // and take the description and the date from meta.
sascha@3435:                 for (int i = 0; i < qs.length; ++i) {
sascha@3435:                     if (Math.abs(qs[i]-q) < EPSILON
sascha@3435:                     &&  Math.abs(ws[i]-w) < EPSILON) {
sascha@3435:                         Column column = eventColumns.get(i);
sascha@3435:                         return new QWD(
sascha@3435:                             qs[i], ws[i],
sascha@3435:                             column.getDescription(),
sascha@3435:                             column.getDate(),
sascha@3435:                             interpolated[i],
sascha@3604:                             0d,
teichmann@6875:                             column.getId()); // Use database id here
sascha@3435:                     }
sascha@3435:                 }
sascha@3435:                 log.warn("cannot find column for (" + q + ", " + w + ")");
sascha@3435:                 return new QWD(q, w);
sascha@3435:             }
sascha@3435:         };
sascha@3435: 
sascha@3435:         Fitting fitting = new Fitting(func, qwdFactory, preprocessing);
sascha@3435: 
sascha@3435:         String [] parameterNames = func.getParameterNames();
sascha@3435: 
sascha@3435:         Parameters results =
sascha@3435:             new Parameters(
sascha@3435:                 StringUtils.join(STANDARD_COLUMNS, parameterNames));
sascha@3435: 
sascha@3435:         boolean invalid = false;
sascha@3435: 
sascha@3435:         double [] kms = DoubleUtil.explode(from, to, step / 1000.0);
sascha@3435: 
sascha@3435:         if (debug) {
sascha@3435:             log.debug("number of kms: " + kms.length);
sascha@3435:         }
sascha@3435: 
sascha@3729:         KMIndex<QWI []> outliers   = new KMIndex<QWI []>();
sascha@3435:         KMIndex<QWD []> referenced = new KMIndex<QWD []>(kms.length);
sascha@3435: 
sascha@3435:         int kmIndex             = results.columnIndex("km");
sascha@3435:         int chiSqrIndex         = results.columnIndex("chi_sqr");
sascha@3435:         int maxQIndex           = results.columnIndex("max_q");
sascha@3435:         int stdDevIndex         = results.columnIndex("std-dev");
sascha@3435:         int [] parameterIndices = results.columnIndices(parameterNames);
sascha@3435: 
sascha@3435:         int numFailed = 0;
sascha@3435: 
sascha@3435:         for (int i = 0; i < kms.length; ++i) {
sascha@3435:             double km = kms[i];
sascha@3435: 
sascha@3435:             // Fill Qs and Ws from event columns.
sascha@3435:             for (int j = 0; j < ws.length; ++j) {
sascha@3435:                 interpolated[j] = !eventColumns.get(j).getQW(km, qs, ws, j);
sascha@3435:             }
sascha@3435: 
sascha@3435:             fitting.reset();
sascha@3435: 
sascha@3435:             if (!fitting.fit(qs, ws)) {
aheinecke@7300:                 log.debug("Fitting for km: " + km + " failed");
sascha@3435:                 ++numFailed;
sascha@3435:                 addProblem(km, "fix.fitting.failed");
sascha@3435:                 continue;
sascha@3435:             }
sascha@3435: 
teichmann@6875:             QWD [] refs = fitting.referencedToArray();
teichmann@6875: 
teichmann@6875:             referenced.add(km, refs);
sascha@3435: 
sascha@3435:             if (fitting.hasOutliers()) {
sascha@3435:                 outliers.add(km, fitting.outliersToArray());
sascha@3435:             }
sascha@3435: 
sascha@3435:             int row = results.newRow();
sascha@3435:             double [] values = fitting.getParameters();
sascha@3435: 
sascha@3435:             results.set(row, kmIndex, km);
sascha@3435:             results.set(row, chiSqrIndex, fitting.getChiSquare());
sascha@3435:             results.set(row, stdDevIndex, fitting.getStandardDeviation());
sascha@3435:             results.set(row, maxQIndex,   fitting.getMaxQ());
sascha@3435:             invalid |= results.set(row, parameterIndices, values);
sascha@3435: 
sascha@3435:             if (debug) {
sascha@3435:                 log.debug("km: "+km+" " + toString(parameterNames, values));
sascha@3435:             }
sascha@3435:         }
sascha@3435: 
sascha@3435:         if (debug) {
sascha@3435:             log.debug("success: " + (kms.length - numFailed));
sascha@3435:             log.debug("failed: " + numFailed);
sascha@3435:         }
sascha@3435: 
sascha@3435:         if (invalid) {
sascha@3435:             addProblem("fix.invalid.values");
sascha@3435:             results.removeNaNs();
sascha@3435:         }
sascha@3435: 
sascha@3435:         outliers.sort();
sascha@3435:         referenced.sort();
sascha@3435: 
sascha@3435:         return new FitResult(
sascha@3435:             results,
sascha@3435:             referenced,
sascha@3435:             outliers);
sascha@3435:     }
sascha@3437: 
sascha@3437:     public CalculationResult calculate() {
sascha@2729:         FixingsOverview overview =
sascha@2729:             FixingsOverviewFactory.getOverview(river);
sascha@2729: 
sascha@2729:         if (overview == null) {
sascha@2729:             addProblem("fix.no.overview.available");
sascha@2729:         }
sascha@2729: 
sascha@2729:         Function func = FunctionFactory.getInstance()
sascha@2729:             .getFunction(function);
sascha@2729: 
sascha@2729:         if (func == null) {
sascha@2729:             addProblem("fix.invalid.function.name");
sascha@2729:         }
sascha@2729: 
sascha@2729:         if (hasProblems()) {
sascha@2729:             return new CalculationResult(this);
sascha@2729:         }
teichmann@7525:         CalculationResult result = innerCalculate(overview, func);
sascha@2729: 
teichmann@7525:         if (result != null) {
teichmann@7525:             // Workaraound to deal with same dates in data set
teichmann@7525:             Object o = result.getData();
teichmann@7525:             if (o instanceof FixResult) {
teichmann@7525:                 FixResult fr = (FixResult)o;
teichmann@7525:                 fr.makeReferenceEventsDatesUnique();
teichmann@7525:                 fr.remapReferenceIndicesToRank();
teichmann@7525:             }
teichmann@7525:         }
teichmann@7525: 
teichmann@7525:         return result;
sascha@2744:     }
sascha@2744: 
sascha@3437:     protected abstract CalculationResult innerCalculate(
sascha@2744:         FixingsOverview overview,
sascha@3437:         Function        function
sascha@3437:     );
sascha@2729: }
sascha@2729: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :