view flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/fixings/FixCalculation.java @ 2792:fe987587ebc9

Merged revisions 4539-4540,4543,4545-4546 via svnmerge from file:///home/clients/bsh/bsh-generischer-viewer/Material/SVN/flys-artifacts/trunk ........ r4539 | teichmann | 2012-05-27 20:02:13 +0200 (So, 27 Mai 2012) | 1 line FixA: Added forgotten csv/report facets/generators to conf. ........ r4540 | teichmann | 2012-05-27 20:11:31 +0200 (So, 27 Mai 2012) | 1 line FixA: Fixed class cast bug in report facet. ........ r4543 | teichmann | 2012-05-28 20:35:01 +0200 (Mo, 28 Mai 2012) | 1 line FixA: Added facet to return delta w/t as CSV ........ r4545 | teichmann | 2012-05-28 22:59:27 +0200 (Mo, 28 Mai 2012) | 1 line FixA: Made Delta W/t calculation finally work ........ r4546 | teichmann | 2012-05-28 23:34:24 +0200 (Mo, 28 Mai 2012) | 1 line FixA: corrected fitting (Q->W instead W->Q). ........ flys-artifacts/tags/2.7@4547 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Ingo Weinzierl <ingo.weinzierl@intevation.de>
date Tue, 29 May 2012 04:58:29 +0000
parents 063b784b60b4
children 3cda41b5eb23
line wrap: on
line source
package de.intevation.flys.artifacts.model.fixings;

import de.intevation.flys.artifacts.FixationArtifactAccess;

import de.intevation.flys.artifacts.math.fitting.Function;
import de.intevation.flys.artifacts.math.fitting.FunctionFactory;

import de.intevation.flys.artifacts.model.Calculation;
import de.intevation.flys.artifacts.model.CalculationResult;
import de.intevation.flys.artifacts.model.FixingsColumn;
import de.intevation.flys.artifacts.model.FixingsColumnFactory;
import de.intevation.flys.artifacts.model.FixingsOverview;
import de.intevation.flys.artifacts.model.FixingsOverviewFactory;
import de.intevation.flys.artifacts.model.Parameters;

import de.intevation.flys.artifacts.model.FixingsOverview.Range;
import de.intevation.flys.artifacts.model.FixingsOverview.Fixing;
import de.intevation.flys.artifacts.model.FixingsOverview.IdFilter;
import de.intevation.flys.artifacts.model.FixingsOverview.AndFilter;
import de.intevation.flys.artifacts.model.FixingsOverview.DateRangeFilter;
import de.intevation.flys.artifacts.model.FixingsOverview.SectorRangeFilter;

import de.intevation.flys.utils.DoubleUtil;

import java.util.ArrayList;
import java.util.List;
import java.util.Date;

import org.apache.commons.math.MathException;

import org.apache.commons.math.optimization.fitting.CurveFitter;

import org.apache.commons.math.optimization.general.LevenbergMarquardtOptimizer;

import org.apache.log4j.Logger;

public class FixCalculation
extends      Calculation
{
    private static Logger log = Logger.getLogger(FixCalculation.class);

    protected String    river;
    protected double    from;
    protected double    to;
    protected double    step;
    protected boolean   preprocessing;
    protected String    function;
    protected int []    events;
    protected long [][] analysisPeriods;
    protected int       qSectorStart;
    protected int       qSectorEnd;

    public FixCalculation() {
    }

    public FixCalculation(FixationArtifactAccess access) {

        String    river           = access.getRiver();
        Double    from            = access.getFrom();
        Double    to              = access.getTo();
        Double    step            = access.getStep();
        String    function        = access.getFunction();
        int []    events          = access.getEvents();
        long [][] analysisPeriods = access.getAnalysisPeriods();
        Integer   qSectorStart    = access.getQSectorStart();
        Integer   qSectorEnd      = access.getQSectorEnd();

        if (river == null) {
            // TODO: i18n
            addProblem("fix.missing.river");
        }

        if (from == null) {
            // TODO: i18n
            addProblem("fix.missing.from");
        }

        if (to == null) {
            // TODO: i18n
            addProblem("fix.missing.to");
        }

        if (step == null) {
            // TODO: i18n
            addProblem("fix.missing.step");
        }

        if (function == null) {
            // TODO: i18n
            addProblem("fix.missing.function");
        }

        if (events == null || events.length < 1) {
            // TODO: i18n
            addProblem("fix.missing.events");
        }

        if (analysisPeriods == null || analysisPeriods.length < 1) {
            // TODO: i18n
            addProblem("fix.missing.analysis.periods");
        }

        if (qSectorStart == null) {
            // TODO: i18n
            addProblem("fix.missing.qstart.sector");
        }

        if (qSectorEnd == null) {
            // TODO: i18n
            addProblem("fix.missing.qend.sector");
        }

        if (!hasProblems()) {
            this.river           = river;
            this.from            = from;
            this.to              = to;
            this.step            = step;
            this.function        = function;
            this.events          = events;
            this.analysisPeriods = analysisPeriods;
            this.qSectorStart    = qSectorStart;
            this.qSectorEnd      = qSectorEnd;
        }
    }

    public CalculationResult calculate() {

        boolean debug = log.isDebugEnabled();

        FixingsOverview overview =
            FixingsOverviewFactory.getOverview(river);

        if (overview == null) {
            addProblem("fix.no.overview.available");
        }

        Function func = FunctionFactory.getInstance()
            .getFunction(function);

        if (func == null) {
            // TODO: i18n
            addProblem("fix.invalid.function.name");
        }

        if (hasProblems()) {
            return new CalculationResult(this);
        }

        FixingsColumnFactory fcf = FixingsColumnFactory.getInstance();

        List<FixingsColumn> dataColumns =
            new ArrayList<FixingsColumn>(events.length);

        for (int eventId: events) {
            IdFilter idFilter = new IdFilter(eventId);

            List<Fixing.Column> columns = overview.filter(null, idFilter);
            if (columns.isEmpty()) {
                // TODO: i18n
                addProblem("fix.missing.column", eventId);
                continue;
            }
            FixingsColumn dataColumn = fcf.getColumnData(columns.get(0));
            if (dataColumn == null) {
                // TODO: i18n
                addProblem("fix.cannot.load.data", eventId);
                continue;
            }
            dataColumns.add(dataColumn);
        }

        if (dataColumns.size() < 2) {
            // TODO: i18n
            addProblem("fix.too.less.data.columns");
            return new CalculationResult(this);
        }

        double [] kms = DoubleUtil.explode(from, to, step / 1000.0);

        double [] ws = new double[dataColumns.size()];
        double [] qs = new double[ws.length];

        String [] parameterNames = func.getParameterNames();

        Parameters results =
            new Parameters(createColumnNames(parameterNames));

        boolean invalid = false;

        if (debug) {
            log.debug("number of kms: " + kms.length);
        }

        int numFailed = 0;

        for (int i = 0; i < kms.length; ++i) {
            double km = kms[i];

            for (int j = 0; j < ws.length; ++j) {
                FixingsColumn column = dataColumns.get(j);
                qs[j] = column.getQ(km);
                boolean interpolated = column.getW(km, ws, j);
                // TODO: mark as interpolated.
            }

            // TODO: Do preprocessing here!
            double [] parameters = fit(func, km, ws, qs);
            if (parameters == null) { // Problems are reported already.
                ++numFailed;
                continue;
            }

            int row = results.newRow();

            results.set(row, "km", km);
            for (int j = 0; j < parameters.length; ++j) {
                if (Double.isNaN(parameters[j])) {
                    invalid = true;
                }
                else {
                    results.set(row, parameterNames[j], parameters[j]);
                }
            }
            // TODO: Calculate statistics, too!
        }

        if (debug) {
            log.debug("success: " + (kms.length - numFailed));
            log.debug("failed: " + numFailed);
        }

        if (invalid) {
            // TODO: i18n
            addProblem("fix.invalid.values");
            results.removeNaNs();
        }

        // Calculate Delta W/t
        DeltaWTsKM deltaWTsKM = calculateDeltaWTs(
            func,
            overview,
            results);

        FixResult fr = new FixResult(results, deltaWTsKM);

        return new CalculationResult(fr, this);
    }

    public DeltaWTsKM calculateDeltaWTs(
        Function        function,
        FixingsOverview overview,
        Parameters      results
    ) {
        boolean debug = log.isDebugEnabled();

        DeltaWTsKM deltaWTsKM = new DeltaWTsKM(results.size());

        Column [][] analysisColumns = getAnalysisColumns(overview);

        int [] parameterIndices = 
            results.columnIndices(function.getParameterNames());

        double [] parameterValues =
            new double[parameterIndices.length];

        double [] ow = new double[1];

        int kmIdx = results.columnIndex("km");

        for (int i = 0, N = results.size(); i < N; ++i) {
            double km = results.get(i, kmIdx);
            results.get(i, parameterIndices, parameterValues);

            DeltaWTsKM.KM dwtkm = new DeltaWTsKM.KM(km);
            deltaWTsKM.add(dwtkm);

            // This is the paraterized function for a given km.
            de.intevation.flys.artifacts.math.Function instance =
                function.instantiate(parameterValues);

            // Evaluate all columns for all analysis periods.
            for (int j = 0; j < analysisColumns.length; ++j) {
                Column [] periodColumns = analysisColumns[j];

                int failedQ = 0;
                int failedW = 0;
                int failedC = 0;

                for (int k = 0; k < periodColumns.length; ++k) {
                    Column pc = periodColumns[k];

                    // Q from real data.
                    double q = pc.data.getQ(km);
                    if (Double.isNaN(q)) {
                        ++failedQ;
                        continue;
                    }

                    // Calculate W from function.
                    double nw = instance.value(q);
                    if (Double.isNaN(nw)) {
                        ++failedC;
                        continue;
                    }

                    // W from real data.
                    pc.data.getW(km, ow);

                    if (Double.isNaN(ow[0])) {
                        ++failedW;
                        continue;
                    }

                    double deltaW = (ow[0] - nw)*100.0; // in cm

                    DeltaWT deltaWT = new DeltaWT(
                        deltaW,
                        pc.meta.getStartTime(),
                        pc.meta.getDescription());

                    dwtkm.add(deltaWT);
                }
                if (debug) {
                    log.debug("failed W: " + failedW);
                    log.debug("failed Q: " + failedQ);
                    log.debug("failed C: " + failedC);
                    log.debug("input size: " + periodColumns.length);
                    log.debug("outpt size: " + dwtkm.size());
                }
            }
        }

        return deltaWTsKM;
    }

    /** Helper class to bundle the meta information of a column 
     *  and the real data.
     */
    protected static class Column {

        protected Fixing.Column meta;
        protected FixingsColumn data;

        public Column() {
        }

        public Column(Fixing.Column meta, FixingsColumn data) {
            this.meta = meta;
            this.data = data;
        }
    } // class Column

    /** Fetch meta and data columns for analysis periods. */
    protected Column [][] getAnalysisColumns(FixingsOverview overview) {

        boolean debug = log.isDebugEnabled();
        if (debug) {
            log.debug("number analysis periods: " + analysisPeriods.length);
        }

        Column columns [][] = new Column[analysisPeriods.length][];

        Range range = new Range(from, to);
        SectorRangeFilter sectorRangeFilter =
            new SectorRangeFilter(qSectorStart, qSectorEnd);

        FixingsColumnFactory fcf = FixingsColumnFactory.getInstance();

        for (int i = 0; i < columns.length; ++i) {

            // Construct filter for period.
            long [] period = analysisPeriods[i];

            AndFilter filter = new AndFilter();

            DateRangeFilter dateRangeFilter =
                new DateRangeFilter(
                    new Date(period[0]),
                    new Date(period[1]));

            filter.add(dateRangeFilter);
            filter.add(sectorRangeFilter);

            List<Fixing.Column> metaCols = overview.filter(range, filter);

            if (debug) {
                log.debug("number of filtered columns: " + metaCols.size());
            }

            ArrayList<Column> cols = new ArrayList<Column>(metaCols.size());

            // Only use columns which have data.
            for (Fixing.Column meta: metaCols) {
                FixingsColumn data = fcf.getColumnData(meta);
                if (data != null) {
                    cols.add(new Column(meta, data));
                }
            }

            if (debug) {
                log.debug("failed loading: " + (metaCols.size()-cols.size()));
            }
            columns[i] = cols.toArray(new Column[cols.size()]);
        }

        return columns;
    }

    protected static String [] createColumnNames(String [] parameters) {
        String [] result = new String[parameters.length + 1];
        result[0] = "km";
        // TODO: Add statistic columns, too.
        System.arraycopy(parameters, 0, result, 1, parameters.length);
        return result;
    }

    protected double [] fit(
        Function  function,
        double    km,
        double [] ws, 
        double [] qs
    ) {
        LevenbergMarquardtOptimizer lmo = new LevenbergMarquardtOptimizer();
        CurveFitter cf = new CurveFitter(lmo);

        for (int i = 0; i < ws.length; ++i) {
            if (!Double.isNaN(ws[i]) && !Double.isNaN(qs[i])) {
                cf.addObservedPoint(qs[i], ws[i]);
            }
        }

        try {
            return cf.fit(function, function.getInitialGuess());
        }
        catch (MathException me) {
            log.warn(me, me);
            addProblem(km, "fix.fitting.failed");
        }

        return null;
    }
}
// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :

http://dive4elements.wald.intevation.org