view flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/DischargeTables.java @ 333:67b3f54188aa

Discharge tables: Sorting of q values was done wrong. flys-artifacts/trunk@1731 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Tue, 19 Apr 2011 14:48:28 +0000
parents 0b2358bc716d
children b7c8df643dc4
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;
    }

    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 and kind=:kind");
        gaugeQuery.setParameter("river", riverName);
        gaugeQuery.setInteger("kind", kind);

        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);

            List<DischargeTableValue> dtvs =
                table.getDischargeTableValues();

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

            boolean qSorted = true;

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

            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[0][indices[i]];
                }
                values.put(gaugeName, vs2);
            }
            else {
                values.put(gaugeName, vs);
            }

        }

        return values;
    }

    public static boolean getWForQ(
        double [][] values, 
        double      q,  
        double []   result
    ) {
        int index = Arrays.binarySearch(values[1], q);
        if (index >= 0) {
            result[0] = values[1][index];
            return true;
        }

        index = -index - 1; // insert position

        if (index == 0 || index+1 >= values.length) {
            // do not extraploate
            return false;
        }

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

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

        if (q1 == q2) {
            result[0] = 0.5*(w1 + w2);
        }
        else {
            double m = (w2 - w1)/(q2 - q1);
            double b = w1 - m*q1;
            result[0] = q*m + b;
        }

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

http://dive4elements.wald.intevation.org