view flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/fixings/FixCalculation.java @ 3165:6b13d7e7b403

FixA: Added more debug output. flys-artifacts/trunk@4778 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Mon, 25 Jun 2012 14:36:46 +0000
parents dad513d5ce37
children 1b9f791937c3
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.AndFilter;
import de.intevation.flys.artifacts.model.FixingsOverview.DateRangeFilter;

import de.intevation.flys.artifacts.model.FixingsOverview.Fixing.Filter;

import de.intevation.flys.artifacts.model.FixingsOverview.Fixing;
import de.intevation.flys.artifacts.model.FixingsOverview.IdsFilter;
import de.intevation.flys.artifacts.model.FixingsOverview.KmFilter;
import de.intevation.flys.artifacts.model.FixingsOverview.SectorFilter;
import de.intevation.flys.artifacts.model.FixingsOverview.SectorRangeFilter;

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.Range;

import de.intevation.flys.utils.DateAverager;
import de.intevation.flys.utils.DoubleUtil;
import de.intevation.flys.utils.KMIndex;

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

import org.apache.commons.math.stat.descriptive.moment.StandardDeviation;

import org.apache.log4j.Logger;

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

    public static final double EPSILON = 1e-4;

    protected String       river;
    protected double       from;
    protected double       to;
    protected double       step;
    protected boolean      preprocessing;
    protected String       function;
    protected int []       events;
    protected DateRange    referencePeriod;
    protected DateRange [] 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();
        DateRange    referencePeriod = access.getReferencePeriod();
        DateRange [] analysisPeriods = access.getAnalysisPeriods();
        Integer      qSectorStart    = access.getQSectorStart();
        Integer      qSectorEnd      = access.getQSectorEnd();
        Boolean      preprocessing   = access.getPreprocessing();

        if (river == null) {
            addProblem("fix.missing.river");
        }

        if (from == null) {
            addProblem("fix.missing.from");
        }

        if (to == null) {
            addProblem("fix.missing.to");
        }

        if (step == null) {
            addProblem("fix.missing.step");
        }

        if (function == null) {
            addProblem("fix.missing.function");
        }

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

        if (referencePeriod == null) {
            addProblem("fix.missing.reference.period");
        }

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

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

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

        if (preprocessing == null) {
            addProblem("fix.missing.preprocessing");
        }

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

    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) {
            addProblem("fix.invalid.function.name");
        }

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

        final List<Column> eventColumns = getEventColumns(overview);

        if (eventColumns.size() < 2) {
            addProblem("fix.too.less.data.columns");
            return new CalculationResult(this);
        }

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

        final double  [] qs = new double[eventColumns.size()];
        final double  [] ws = new double[qs.length];
        final boolean [] interpolated = new boolean[ws.length];

        Fitting.QWDFactory qwdFactory = new Fitting.QWDFactory() {
            @Override
            public QWD create(double q, double w) {
                // Check all the event columns for close match
                // and take the description and the date from meta.
                for (int i = 0; i < qs.length; ++i) {
                    if (Math.abs(qs[i]-q) < EPSILON
                    &&  Math.abs(ws[i]-w) < EPSILON) {
                        Column column = eventColumns.get(i);
                        return new QWD(
                            qs[i], ws[i],
                            column.getDescription(),
                            column.getDate(),
                            interpolated[i],
                            0d);
                    }
                }
                log.warn("cannot find column for (" + q + ", " + w + ")");
                return new QWD(q, w);
            }
        };

        Fitting fitting = new Fitting(func, qwdFactory, preprocessing);

        String [] parameterNames = func.getParameterNames();

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

        boolean invalid = false;

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

        KMIndex<QW []>  outliers   = new KMIndex<QW []>();
        KMIndex<QWD []> referenced = new KMIndex<QWD []>(kms.length);

        int kmIndex             = results.columnIndex("km");
        int chiSqrIndex         = results.columnIndex("chi_sqr");
        int maxQIndex           = results.columnIndex("max_q");
        int stdDevIndex         = results.columnIndex("std-dev");
        int [] parameterIndices = results.columnIndices(parameterNames);

        int numFailed = 0;

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

            // Fill Qs and Ws from event columns.
            for (int j = 0; j < ws.length; ++j) {
                interpolated[j] = !eventColumns.get(j).getQW(km, qs, ws, j);
            }

            fitting.reset();

            if (!fitting.fit(qs, ws)) {
                ++numFailed;
                addProblem(km, "fix.fitting.failed");
                continue;
            }

            referenced.add(km, fitting.referencedToArray());

            if (fitting.hasOutliers()) {
                outliers.add(km, fitting.outliersToArray());
            }

            int row = results.newRow();
            double [] values = fitting.getParameters();

            results.set(row, kmIndex, km);
            results.set(row, chiSqrIndex, fitting.getChiSquare());
            results.set(row, stdDevIndex, fitting.getStandardDeviation());
            results.set(row, maxQIndex,   fitting.getMaxQ());
            invalid |= results.set(row, parameterIndices, values);

            if (debug) {
                log.debug("km: "+km+" " + toString(parameterNames, values));
            }
        }

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

        if (invalid) {
            addProblem("fix.invalid.values");
            results.removeNaNs();
        }

        KMIndex<AnalysisPeriod []> analysisPeriods =
            calculateAnalysisPeriods(func, results, overview);

        outliers.sort();
        referenced.sort();
        analysisPeriods.sort();

        FixResult fr = new FixResult(
            results,
            referenced, outliers,
            analysisPeriods);

        return new CalculationResult(fr, this);
    }

    protected String toString(String [] parameterNames, double [] values) {
        StringBuilder sb = new StringBuilder();
        for (int i = 0; i < parameterNames.length; ++i) {
            if (i > 0) sb.append(", ");
            sb.append(parameterNames[i]).append(": ").append(values[i]);
        }
        return sb.toString();
    }

    protected List<Column> getEventColumns(FixingsOverview overview) {

        FixingsColumnFactory fcf = FixingsColumnFactory.getInstance();

        IdsFilter ids = new IdsFilter(events);
        DateRangeFilter rdf = new DateRangeFilter(
            referencePeriod.getFrom(),
            referencePeriod.getTo());
        Filter filter = new AndFilter().add(rdf).add(ids);

        List<Fixing.Column> metas = overview.filter(null, filter);

        List<Column> columns = new ArrayList<Column>(metas.size());

        for (Fixing.Column meta: metas) {

            FixingsColumn data = fcf.getColumnData(meta);
            if (data == null) {
                addProblem("fix.cannot.load.data");
            }
            else {
                columns.add(new Column(meta, data));
            }
        }

        return columns;
    }


    protected KMIndex<AnalysisPeriod []> calculateAnalysisPeriods(
        Function        function,
        Parameters      parameters,
        FixingsOverview overview
    ) {
        Range range = new Range(from, to);

        int kmIndex   = parameters.columnIndex("km");
        int maxQIndex = parameters.columnIndex("max_q");

        ColumnCache cc = new ColumnCache();

        double [] wq = new double[2];

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

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

        DateAverager dateAverager = new DateAverager();

        KMIndex<AnalysisPeriod []> results =
            new KMIndex<AnalysisPeriod []>(parameters.size());

        IdsFilter idsFilter = new IdsFilter(events);

        for (int row = 0, R = parameters.size(); row < R; ++row) {
            double km = parameters.get(row, kmIndex);
            parameters.get(row, parameterIndices, parameterValues);

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

            KmFilter kmFilter = new KmFilter(km);

            ArrayList<AnalysisPeriod> periodResults =
                new ArrayList<AnalysisPeriod>(analysisPeriods.length);

            for (DateRange analysisPeriod: analysisPeriods) {
                DateRangeFilter drf = new DateRangeFilter(
                    analysisPeriod.getFrom(),
                    analysisPeriod.getTo());

                QWD []    qSectorAverages = new QWD[4];
                double [] qSectorStdDevs  = new double[4];

                ArrayList<QWD> allQWDs = new ArrayList<QWD>();

                // for all Q sectors.
                for (int qSector = qSectorStart; qSector < qSectorEnd; ++qSector) {

                    Filter filter = new AndFilter()
                        .add(kmFilter)
                        .add(new SectorFilter(qSector))
                        .add(drf)
                        .add(idsFilter);

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

                    if (metas.isEmpty()) {
                        // No fixings for km and analysis period
                        continue;
                    }

                    double sumQ = 0.0;
                    double sumW = 0.0;

                    StandardDeviation stdDev = new StandardDeviation();

                    List<QWD> qwds = new ArrayList<QWD>(metas.size());

                    dateAverager.clear();

                    for (Fixing.Column meta: metas) {
                        if (meta.findQSector(km) != qSector) {
                            // Ignore not matching sectors.
                            continue;
                        }

                        Column column = cc.getColumn(meta);
                        if (column == null || !column.getQW(km, wq)) {
                            continue;
                        }

                        double fw = instance.value(wq[1]);
                        if (Double.isNaN(fw)) {
                            continue;
                        }

                        double dw = (wq[0] - fw)*100.0;

                        stdDev.increment(dw);

                        Date date = column.getDate();
                        String description = column.getDescription();

                        QWD qwd = new QWD(
                            wq[1], wq[0], description, date, true, dw);

                        qwds.add(qwd);

                        sumW += wq[0];
                        sumQ += wq[1];

                        dateAverager.add(date);
                    }

                    // Calulate average per Q sector.
                    int N = qwds.size();
                    if (N > 0) {
                        allQWDs.addAll(qwds);
                        double avgW = sumW / N;
                        double avgQ = sumQ / N;

                        double avgFw = instance.value(avgQ);
                        if (!Double.isNaN(avgFw)) {
                            double avgDw = (avgW - avgFw)*100.0;
                            Date avgDate = dateAverager.getAverage();

                            String avgDescription = "avg.deltawt." + qSector;

                            QWD avgQWD = new QWD(
                                avgQ, avgW, avgDescription, avgDate, true, avgDw);

                            qSectorAverages[qSector] = avgQWD;
                        }
                        qSectorStdDevs[qSector] = stdDev.getResult();
                    }
                    else {
                        qSectorStdDevs[qSector] = Double.NaN;
                    }
                } // for all Q sectors

                QWD [] aqwds = allQWDs.toArray(new QWD[allQWDs.size()]);

                AnalysisPeriod periodResult = new AnalysisPeriod(
                    analysisPeriod,
                    aqwds,
                    qSectorAverages,
                    qSectorStdDevs);
                periodResults.add(periodResult);
            }

            double maxQ = -Double.MAX_VALUE;
            for (AnalysisPeriod ap: periodResults) {
                double q = ap.getMaxQ();
                if (q > maxQ) {
                    maxQ = q;
                }
            }

            double oldMaxQ = parameters.get(row, maxQIndex);
            if (oldMaxQ < maxQ) {
                parameters.set(row, maxQIndex, maxQ);
            }

            results.add(km, periodResults.toArray(
                new AnalysisPeriod[periodResults.size()]));
        }

        return results;
    }

    /** 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;
        }

        public Date getDate() {
            return meta.getStartTime();
        }

        public String getDescription() {
            return meta.getDescription();
        }

        public boolean getQW(
            double    km,
            double [] qs,
            double [] ws,
            int       index
        ) {
            qs[index] = data.getQ(km);
            return data.getW(km, ws, index);
        }

        public boolean getQW(double km, double [] wq) {
            data.getW(km, wq, 0);
            if (Double.isNaN(wq[0])) return false;
            wq[1] = data.getQ(km);
            return !Double.isNaN(wq[1]);
        }
    } // class Column


    /**
     * Helper class to find the data belonging to meta info more quickly.
     */
    public static class ColumnCache {

        protected Map<Integer, Column> columns;

        public ColumnCache() {
            columns = new HashMap<Integer, Column>();
        }

        public Column getColumn(Fixing.Column meta) {
            Integer key = meta.getId();
            Column column = columns.get(key);
            if (column == null) {
                FixingsColumn data = FixingsColumnFactory
                    .getInstance()
                    .getColumnData(meta);
                if (data != null) {
                    column = new Column(meta, data);
                    columns.put(key, column);
                }
            }
            return column;
        }
    } // class ColumnCache

    /** 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.
            DateRange period = analysisPeriods[i];

            AndFilter filter = new AndFilter();

            DateRangeFilter dateRangeFilter =
                new DateRangeFilter(period.getFrom(), period.getTo());

            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 + 4];
        result[0] = "km";
        result[1] = "chi_sqr";
        result[2] = "max_q";
        result[3] = "std-dev";
        System.arraycopy(parameters, 0, result, 4, parameters.length);
        return result;
    }
}
// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :

http://dive4elements.wald.intevation.org