view flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/DischargeTables.java @ 2230:59af81364eb1

Improved the 'historical discharge' calculation: implemented findValueForW(). flys-artifacts/trunk@3872 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Ingo Weinzierl <ingo.weinzierl@intevation.de>
date Wed, 01 Feb 2012 15:41:11 +0000
parents 6b0ae0f2cae6
children 899ca89f497e
line wrap: on
line source
package de.intevation.flys.artifacts.model;

import java.util.List;
import java.util.Map;
import java.util.HashMap;
import java.util.Arrays;
import java.util.Comparator;

import java.io.Serializable;

import org.hibernate.Session;
import org.hibernate.Query;

import org.apache.log4j.Logger;

import de.intevation.flys.backend.SessionHolder;
import de.intevation.flys.model.Gauge;
import de.intevation.flys.model.DischargeTable;
import de.intevation.flys.model.DischargeTableValue;

public class DischargeTables
implements   Serializable
{
    private static Logger log = Logger.getLogger(DischargeTables.class);

    public static final double DEFAULT_SCALE = 100.0;

    public static final int MASTER = 0;

    protected List<String> gaugeNames;

    protected String riverName;

    protected double scale;

    protected int    kind;

    protected Map<String, double [][]> values;

    public DischargeTables() {
    }

    public DischargeTables(String riverName, String gaugeName) {
        this(riverName, gaugeName, MASTER);
    }

    public DischargeTables(String riverName, String gaugeName, int kind) {
        this(riverName, new String [] { gaugeName }, kind);
    }

    public DischargeTables(String riverName, String [] gaugeNames) {
        this(riverName, gaugeNames, MASTER);
    }

    public DischargeTables(String riverName, String [] gaugeNames, int kind) {
        this(riverName, Arrays.asList(gaugeNames), kind);
    }

    public DischargeTables(
        String       riverName,
        List<String> gaugeNames,
        int          kind
    ) {
        scale           = Double.NaN;
        this.kind       = kind;
        this.riverName  = riverName;
        this.gaugeNames = gaugeNames;
    }

    public double [][] getFirstTable() {
        return getFirstTable(DEFAULT_SCALE);
    }

    public double [][] getFirstTable(double scale) {
        Map<String, double [][]> values = getValues(scale);
        for (double [][] table: values.values()) {
            return table;
        }
        return null;
    }

    public Map<String, double [][]> getValues() {
        return getValues(DEFAULT_SCALE);
    }

    public Map<String, double [][]> getValues(double scale) {
        if (values == null || scale != this.scale) {
            values = loadValues(scale);
            this.scale = scale;
        }
        return values;
    }

    /**
     * Returns mapping of gauge name to values.
     */
    protected Map<String, double [][]> loadValues(double scale) {
        Map<String, double [][]> values = new HashMap<String, double [][]>();

        Session session = SessionHolder.HOLDER.get();

        Query gaugeQuery = session.createQuery(
            "from Gauge where name=:gauge and river.name=:river");
        gaugeQuery.setParameter("river", riverName);

        for (String gaugeName: gaugeNames) {
            gaugeQuery.setParameter("gauge", gaugeName);
            List<Gauge> gauges = gaugeQuery.list();
            if (gauges.isEmpty()) {
                log.warn(
                    "no gauge '"+gaugeName+"' at river '"+riverName+"'");
                continue;
            }
            Gauge gauge = gauges.get(0);

            List<DischargeTable> tables = gauge.getDischargeTables();

            if (tables.isEmpty()) {
                log.warn(
                    "no discharge table for gauge '" + gaugeName + "'");
                continue;
            }

            // TODO: Filter by time interval
            DischargeTable table = tables.get(0);

            boolean qSorted = true;

            final double[][] vs =
                loadDischargeTableValues(table, scale, qSorted);

            if (!qSorted) {
                log.debug("need to sort by q values");
                // TODO: Do this db level.
                // XXX: This is so ugly :-(
                Integer [] indices = new Integer[vs[0].length];
                for (int i = 0; i < indices.length; ++i) {
                    indices[i] = i;
                }

                Arrays.sort(indices, new Comparator<Integer>() {
                    public int compare(Integer a, Integer b) {
                        double va = vs[1][a];
                        double vb = vs[1][b];
                        double d = va - vb;
                        if (d < 0.0) return -1;
                        if (d > 0.0) return +1;
                        return 0;
                    }
                });

                double [][] vs2 = new double[2][vs[0].length];
                for (int i = 0; i < indices.length; ++i) {
                    vs2[0][i] = vs[0][indices[i]];
                    vs2[1][i] = vs[1][indices[i]];
                }
                values.put(gaugeName, vs2);
            }
            else {
                values.put(gaugeName, vs);
            }

        }

        return values;
    }


    /**
     * @param table The discharge table
     * @param scale The scale factor to adjust W and Q values.
     * @param qSorted Boolean flag that is used to evaluate if the values should
     * be sorted.
     *
     * @return the values of a discharge table.
     */
    public static double[][] loadDischargeTableValues(
        DischargeTable table,
        double         scale,
        boolean        qSorted
    ) {
        List<DischargeTableValue> dtvs = table.getDischargeTableValues();

        final double [][] vs = new double[2][dtvs.size()];

        double lastQ = -Double.MAX_VALUE;

        int idx = 0;
        for (DischargeTableValue dtv: dtvs) {
            double q = dtv.getQ().doubleValue();
            vs[0][idx] = q * scale;
            vs[1][idx] = dtv.getW().doubleValue() * scale;
            ++idx;

            if (qSorted && lastQ > q) {
                qSorted = false;
            }
            lastQ = q;
        }

        return vs;
    }


    public static double getQForW(double [][] values, double w) {

        boolean debug = log.isDebugEnabled();

        if (debug) {
            log.debug("calculating getQForW(" + w + ")");
        }

        int index = Arrays.binarySearch(values[1], w);
        if (index >= 0) {
            return values[0][index];
        }

        index = -index - 1; // insert position

        if (index < 1 || index >= values[0].length) {
            // do not extraploate
            if (debug) {
                log.debug("we do not extrapolate: NaN");
            }
            return Double.NaN;
        }

        double w1 = values[1][index-1];
        double w2 = values[1][index  ];
        double q1 = values[0][index-1];
        double q2 = values[0][index  ];

        // q1 = m*w1 + b
        // q2 = m*w2 + b
        // q2 - q1 = m*(w2 - w1)
        // m = (q2 - q1)/(w2 - w1) # w2 != w1
        // b = q1 - m*w1

        double q;
        if (w1 == w2) {
            q = 0.5*(q1 + q2);
            if (debug) {
                log.debug("same w1 and w2: " + w1);
            }
        }
        else {
            double m = (q2 - q1)/(w2 - w1);
            double b = q1 - m*w1;
            q = w*m + b;
        }
        if (debug) {
            log.debug("Q(" + w + ") = " + q);
        }
        return q;
    }
}
// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :

http://dive4elements.wald.intevation.org