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: 
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;
gernotbelger@9099: import org.dive4elements.artifacts.common.utils.StringUtils;
gernotbelger@9099: import org.dive4elements.river.artifacts.access.FixAccess;
gernotbelger@9099: import org.dive4elements.river.artifacts.math.fitting.Function;
gernotbelger@9099: import org.dive4elements.river.artifacts.math.fitting.FunctionFactory;
gernotbelger@9099: import org.dive4elements.river.artifacts.model.Calculation;
gernotbelger@9099: import org.dive4elements.river.artifacts.model.CalculationResult;
gernotbelger@9099: import org.dive4elements.river.artifacts.model.FixingsColumn;
gernotbelger@9099: import org.dive4elements.river.artifacts.model.FixingsColumnFactory;
gernotbelger@9099: import org.dive4elements.river.artifacts.model.FixingsOverview;
gernotbelger@9099: import org.dive4elements.river.artifacts.model.FixingsOverview.Fixing;
gernotbelger@9099: import org.dive4elements.river.artifacts.model.FixingsOverview.Fixing.Filter;
gernotbelger@9099: import org.dive4elements.river.artifacts.model.FixingsOverview.IdsFilter;
gernotbelger@9099: import org.dive4elements.river.artifacts.model.FixingsOverviewFactory;
gernotbelger@9099: import org.dive4elements.river.artifacts.model.Parameters;
gernotbelger@9099: import org.dive4elements.river.utils.DoubleUtil;
gernotbelger@9099: import org.dive4elements.river.utils.KMIndex;
sascha@2729: 
felix@5150: /** Calculation base class for fix. */
gernotbelger@9099: public abstract class FixCalculation extends Calculation {
sascha@2786:     private static Logger log = Logger.getLogger(FixCalculation.class);
sascha@2729: 
sascha@3435:     public static final double EPSILON = 1e-4;
sascha@3435: 
gernotbelger@9099:     public static final String[] STANDARD_COLUMNS = { "km", "chi_sqr", "max_q", "std-dev" };
sascha@3435: 
sascha@3435:     protected static class FitResult {
sascha@3435: 
gernotbelger@9363:         private final Parameters parameters;
gernotbelger@9363:         private final KMIndex<QWD[]> fixings;
sascha@3435: 
gernotbelger@9363:         public FitResult(final Parameters parameters, final KMIndex<QWD[]> fixings) {
sascha@3435:             this.parameters = parameters;
gernotbelger@9363:             this.fixings = fixings;
sascha@3435:         }
sascha@3435: 
sascha@3435:         public Parameters getParameters() {
gernotbelger@9099:             return this.parameters;
sascha@3435:         }
sascha@3435: 
gernotbelger@9363:         public KMIndex<QWD[]> getFixings() {
gernotbelger@9363:             return this.fixings;
sascha@3435:         }
gernotbelger@9363:     }
sascha@3435: 
gernotbelger@9099:     /**
gernotbelger@9099:      * Helper class to bundle the meta information of a column
gernotbelger@9099:      * 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;
gernotbelger@9099:         protected int index;
sascha@3434: 
sascha@3434:         public Column() {
sascha@3434:         }
sascha@3434: 
gernotbelger@9099:         public Column(final Fixing.Column meta, final FixingsColumn data, final int index) {
gernotbelger@9099:             this.meta = meta;
gernotbelger@9099:             this.data = data;
sascha@3603:             this.index = index;
sascha@3434:         }
sascha@3434: 
sascha@3434:         public Date getDate() {
gernotbelger@9099:             return this.meta.getStartTime();
sascha@3434:         }
sascha@3434: 
sascha@3434:         public String getDescription() {
gernotbelger@9099:             return this.meta.getDescription();
sascha@3434:         }
sascha@3434: 
sascha@3603:         public int getIndex() {
gernotbelger@9099:             return this.index;
sascha@3603:         }
sascha@3603: 
teichmann@6875:         public int getId() {
gernotbelger@9099:             return this.meta.getId();
teichmann@6875:         }
teichmann@6875: 
gernotbelger@9099:         public boolean getQW(final double km, final double[] qs, final double[] ws, final int index) {
gernotbelger@9099:             qs[index] = this.data.getQ(km);
gernotbelger@9099:             return this.data.getW(km, ws, index);
sascha@3434:         }
sascha@3434: 
gernotbelger@9099:         public boolean getQW(final double km, final double[] wq) {
gernotbelger@9099:             this.data.getW(km, wq, 0);
gernotbelger@9099:             if (Double.isNaN(wq[0]))
gernotbelger@9099:                 return false;
gernotbelger@9099:             wq[1] = this.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() {
gernotbelger@9099:             this.columns = new HashMap<>();
sascha@3434:         }
sascha@3434: 
gernotbelger@9099:         public Column getColumn(final Fixing.Column meta) {
gernotbelger@9099:             final Integer key = meta.getId();
gernotbelger@9099:             Column column = this.columns.get(key);
sascha@3434:             if (column == null) {
gernotbelger@9099:                 final FixingsColumn data = FixingsColumnFactory.getInstance().getColumnData(meta);
sascha@3434:                 if (data != null) {
gernotbelger@9099:                     column = new Column(meta, data, this.columns.size());
gernotbelger@9099:                     this.columns.put(key, column);
sascha@3434:                 }
sascha@3434:             }
sascha@3434:             return column;
sascha@3434:         }
sascha@3434:     } // class ColumnCache
sascha@3434: 
gernotbelger@9099:     protected String river;
gernotbelger@9099:     protected double from;
gernotbelger@9099:     protected double to;
gernotbelger@9099:     protected double step;
sascha@3419:     protected boolean preprocessing;
gernotbelger@9099:     protected String function;
gernotbelger@9099:     protected int[] events;
gernotbelger@9099:     protected int qSectorStart;
gernotbelger@9099:     protected int qSectorEnd;
sascha@2729: 
sascha@2729:     public FixCalculation() {
sascha@2729:     }
sascha@2729: 
gernotbelger@9099:     public FixCalculation(final FixAccess access) {
gernotbelger@9099:         final String river = access.getRiverName();
gernotbelger@9099:         final Double from = access.getLowerKm();
gernotbelger@9099:         final Double to = access.getUpperKm();
gernotbelger@9099:         final Double step = access.getStep();
gernotbelger@9099:         final String function = access.getFunction();
gernotbelger@9099:         final int[] events = access.getEvents();
gernotbelger@9099:         final Integer qSectorStart = access.getQSectorStart();
gernotbelger@9099:         final Integer qSectorEnd = access.getQSectorEnd();
gernotbelger@9099:         final 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()) {
gernotbelger@9099:             this.river = river;
gernotbelger@9099:             this.from = from;
gernotbelger@9099:             this.to = to;
gernotbelger@9099:             this.step = step;
gernotbelger@9099:             this.function = function;
gernotbelger@9099:             this.events = events;
gernotbelger@9099:             this.qSectorStart = qSectorStart;
gernotbelger@9099:             this.qSectorEnd = qSectorEnd;
sascha@3419:             this.preprocessing = preprocessing;
sascha@2729:         }
sascha@2729:     }
sascha@2729: 
gernotbelger@9099:     protected static String toString(final String[] parameterNames, final double[] values) {
gernotbelger@9099:         final StringBuilder sb = new StringBuilder();
sascha@3434:         for (int i = 0; i < parameterNames.length; ++i) {
gernotbelger@9099:             if (i > 0)
gernotbelger@9099:                 sb.append(", ");
sascha@3434:             sb.append(parameterNames[i]).append(": ").append(values[i]);
sascha@3434:         }
sascha@3434:         return sb.toString();
sascha@3434:     }
sascha@2729: 
gernotbelger@9099:     /**
gernotbelger@9099:      * Create filter to accept only the chosen events.
gernotbelger@9099:      * This factored out out to be overwritten.
teichmann@5741:      */
teichmann@5741:     protected Filter createFilter() {
gernotbelger@9099:         return new IdsFilter(this.events);
sascha@3434:     }
sascha@3434: 
gernotbelger@9099:     protected List<Column> getEventColumns(final FixingsOverview overview, final ColumnCache cc) {
gernotbelger@9099:         final Filter filter = createFilter();
sascha@3434: 
gernotbelger@9099:         final List<Fixing.Column> metas = overview.filter(null, filter);
sascha@3434: 
gernotbelger@9099:         final List<Column> columns = new ArrayList<>(metas.size());
sascha@3434: 
gernotbelger@9099:         for (final Fixing.Column meta : metas) {
gernotbelger@9099: 
gernotbelger@9099:             final Column data = cc.getColumn(meta);
sascha@3434:             if (data == null) {
sascha@3434:                 addProblem("fix.cannot.load.data");
gernotbelger@9099:             } 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.
gernotbelger@9099:     protected FitResult doFitting(final FixingsOverview overview, final ColumnCache cc, final Function func) {
gernotbelger@9099:         final 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: 
gernotbelger@9099:         final double[] qs = new double[eventColumns.size()];
gernotbelger@9099:         final double[] ws = new double[qs.length];
gernotbelger@9099:         final boolean[] interpolated = new boolean[ws.length];
sascha@3435: 
gernotbelger@9099:         final Fitting.QWDFactory qwdFactory = new Fitting.QWDFactory() {
sascha@3435:             @Override
gernotbelger@9363:             public QWD create(final double q, final double w, double deltaW, boolean isOutlier) {
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) {
gernotbelger@9099:                     if (Math.abs(qs[i] - q) < EPSILON && Math.abs(ws[i] - w) < EPSILON) {
gernotbelger@9099:                         final Column column = eventColumns.get(i);
gernotbelger@9363:                         return new QWD(qs[i], ws[i], column.getDescription(), column.getDate(), interpolated[i], deltaW, column.getId(), isOutlier); // Use database id here
sascha@3435:                     }
sascha@3435:                 }
sascha@3435:                 log.warn("cannot find column for (" + q + ", " + w + ")");
gernotbelger@9356:                 return new QWD(q, w, isOutlier);
sascha@3435:             }
sascha@3435:         };
sascha@3435: 
gernotbelger@9099:         final String[] parameterNames = func.getParameterNames();
sascha@3435: 
gernotbelger@9099:         final Parameters results = new Parameters(StringUtils.join(STANDARD_COLUMNS, parameterNames));
sascha@3435: 
sascha@3435:         boolean invalid = false;
sascha@3435: 
gernotbelger@9099:         final double[] kms = DoubleUtil.explode(this.from, this.to, this.step / 1000.0);
sascha@3435: 
sascha@3435:         if (debug) {
sascha@3435:             log.debug("number of kms: " + kms.length);
sascha@3435:         }
sascha@3435: 
gernotbelger@9363:         final KMIndex<QWD[]> fixings = new KMIndex<>();
sascha@3435: 
gernotbelger@9099:         final int kmIndex = results.columnIndex("km");
gernotbelger@9099:         final int chiSqrIndex = results.columnIndex("chi_sqr");
gernotbelger@9099:         final int maxQIndex = results.columnIndex("max_q");
gernotbelger@9099:         final int stdDevIndex = results.columnIndex("std-dev");
gernotbelger@9099:         final int[] parameterIndices = results.columnIndices(parameterNames);
sascha@3435: 
sascha@3435:         int numFailed = 0;
sascha@3435: 
gernotbelger@9363:         for (final double km : kms) {
sascha@3435: 
sascha@3435:             // Fill Qs and Ws from event columns.
gernotbelger@9363:             for (int j = 0; j < ws.length; ++j) 
sascha@3435:                 interpolated[j] = !eventColumns.get(j).getQW(km, qs, ws, j);
sascha@3435: 
gernotbelger@9363:             final Fitting fitting = Fitting.fit(func, qwdFactory, this.preprocessing, qs, ws);
gernotbelger@9363:             if (fitting == null) {
aheinecke@7300:                 log.debug("Fitting for km: " + km + " failed");
sascha@3435:                 ++numFailed;
sascha@3435:                 addProblem(km, "fix.fitting.failed");
sascha@3435:                 continue;
sascha@3435:             }
gernotbelger@9363:             
gernotbelger@9363:             final QWD[] fixingsArray = fitting.getFixingsArray();
gernotbelger@9363:             fixings.add(km, fixingsArray);
sascha@3435: 
gernotbelger@9099:             final int row = results.newRow();
gernotbelger@9099:             final 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());
gernotbelger@9099:             results.set(row, maxQIndex, fitting.getMaxQ());
sascha@3435:             invalid |= results.set(row, parameterIndices, values);
sascha@3435: 
sascha@3435:             if (debug) {
gernotbelger@9099:                 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: 
gernotbelger@9363:         fixings.sort();
sascha@3435: 
gernotbelger@9363:         return new FitResult(results, fixings);
sascha@3435:     }
sascha@3437: 
sascha@3437:     public CalculationResult calculate() {
gernotbelger@9099:         final FixingsOverview overview = FixingsOverviewFactory.getOverview(this.river);
sascha@2729: 
sascha@2729:         if (overview == null) {
sascha@2729:             addProblem("fix.no.overview.available");
sascha@2729:         }
sascha@2729: 
gernotbelger@9099:         final Function func = FunctionFactory.getInstance().getFunction(this.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:         }
gernotbelger@9099:         final CalculationResult result = innerCalculate(overview, func);
sascha@2729: 
teichmann@7525:         if (result != null) {
teichmann@7525:             // Workaraound to deal with same dates in data set
gernotbelger@9099:             final Object o = result.getData();
teichmann@7525:             if (o instanceof FixResult) {
gernotbelger@9099:                 final FixResult fr = (FixResult) o;
gernotbelger@9363:                 fr.makeEventsDatesUnique();
gernotbelger@9363:                 fr.remapEventIndicesToRank();
teichmann@7525:             }
teichmann@7525:         }
teichmann@7525: 
teichmann@7525:         return result;
sascha@2744:     }
sascha@2744: 
gernotbelger@9099:     protected abstract CalculationResult innerCalculate(FixingsOverview overview, Function function);
gernotbelger@9363: }