view flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/Calculation4.java @ 4798:39885bdfc6fc

Added calculation of the "Umhuellende" to calculation of "W fuer ungleichwertige Abfluesse". This is done by figuring out the WST columns that imfold the data and then do simple "gleichwertige" calculations from the start of the interval. This is too much because only the Qs are needed for the "Umhuellende".
author Sascha L. Teichmann <teichmann@intevation.de>
date Sun, 13 Jan 2013 16:18:28 +0100
parents 22790758b132
children e3532c7850b9
line wrap: on
line source
package de.intevation.flys.artifacts.model;

import de.intevation.flys.artifacts.math.BackJumpCorrector;
import de.intevation.flys.artifacts.math.Function;
import de.intevation.flys.artifacts.math.Identity;
import de.intevation.flys.artifacts.math.Linear;

import de.intevation.flys.artifacts.model.WstValueTable.QPosition;

import de.intevation.flys.model.River;

import de.intevation.flys.utils.DoubleUtil;

import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

import org.apache.log4j.Logger;

public class Calculation4
extends      Calculation
{
    private static Logger logger = Logger.getLogger(Calculation4.class);

    public static final double MINIMAL_STEP_WIDTH = 1e-5;

    protected List<Segment> segments;

    protected boolean       isQ;

    public Calculation4() {
    }

    public Calculation4(List<Segment> segments, River river, boolean isQ) {

        this.segments = segments;
        this.isQ      = isQ;

        Segment.setReferencePointConvertQ(segments, river, isQ, this);
    }

    public CalculationResult calculate(
        WstValueTable table,
        double from, double to, double step
    ) {
        boolean debug = logger.isDebugEnabled();

        if (debug) {
            logger.debug(
                "calculate from " + from + " to " + to + " step " + step);
            logger.debug("# segments: " + segments.size());
            for (Segment segment: segments) {
                logger.debug("  " + segment);
            }
        }

        if (segments.isEmpty()) {
            logger.debug("no segments found");
            addProblem("no.segments.found");
            return new CalculationResult(new WQKms[0], this);
        }

        int numResults = segments.get(0).values.length;

        if (numResults < 1) {
            logger.debug("no values given");
            addProblem("no.values.given");
            return new CalculationResult(new WQKms[0], this);
        }


        WQKms [] results = new WQKms[numResults];
        for (int i = 0; i < results.length; ++i) {
            results[i] = new WQKms();
        }

        if (Math.abs(step) < MINIMAL_STEP_WIDTH) {
            step = MINIMAL_STEP_WIDTH;
        }

        if (from > to) {
            step = -step;
        }

        QPosition [] qPositions = new QPosition[numResults];

        Function [] functions = new Function[numResults];

        double [] out = new double[2];

        Segment sentinel = new Segment(Double.MAX_VALUE);
        Segment s1 = sentinel, s2 = sentinel;

        for (double pos = from;
             from < to ? pos <= to : pos >= to;
             pos = DoubleUtil.round(pos + step)
        ) {
            if (pos < s1.referencePoint || pos > s2.referencePoint) {
                if (debug) {
                    logger.debug("need to find new interval for " + pos);
                }
                // find new interval
                if (pos <= segments.get(0).referencePoint) {
                    // before first segment -> "gleichwertig"
                    if (debug) {
                        logger.debug("before first segment -> gleichwertig");
                    }
                    Segment   first  = segments.get(0);
                    double [] values = first.values;
                    double    refPos = first.referencePoint;
                    for (int i = 0; i < qPositions.length; ++i) {
                        qPositions[i] = table.getQPosition(
                            refPos, values[i]);
                    }
                    sentinel.setReferencePoint(-Double.MAX_VALUE);
                    s1 = sentinel;
                    s2 = segments.get(0);
                    Arrays.fill(functions, Identity.IDENTITY);
                }
                else if (pos >= segments.get(segments.size()-1).referencePoint) {
                    // after last segment -> "gleichwertig"
                    if (debug) {
                        logger.debug("after last segment -> gleichwertig");
                    }
                    Segment   last   = segments.get(segments.size()-1);
                    double [] values = last.values;
                    double    refPos = last.referencePoint;
                    for (int i = 0; i < qPositions.length; ++i) {
                        qPositions[i] = table.getQPosition(
                            refPos, values[i]);
                    }
                    sentinel.setReferencePoint(Double.MAX_VALUE);
                    s1 = last;
                    s2 = sentinel;
                    Arrays.fill(functions, Identity.IDENTITY);
                }
                else { // "ungleichwertig"
                    // find matching interval
                    if (debug) {
                        logger.debug("in segments -> ungleichwertig");
                    }
                    s1 = s2 = null;
                    for (int i = 1, N = segments.size(); i < N; ++i) {
                        Segment si1 = segments.get(i-1);
                        Segment si  = segments.get(i);
                        if (debug) {
                            logger.debug("check " + pos + " in " +
                                si1.referencePoint + " - " + si.referencePoint);
                        }
                        if (pos >= si1.referencePoint
                        &&  pos <= si. referencePoint) {
                            s1 = si1;
                            s2 = si;
                            break;
                        }
                    }

                    if (s1 == null) {
                        throw new IllegalStateException("no interval found");
                    }

                    Segment anchor, free;

                    if (from > to) { anchor = s1; free = s2; }
                    else           { anchor = s2; free = s1; }

                    // build transforms based on "gleichwertiger" phase
                    for (int i = 0; i < qPositions.length; ++i) {
                        QPosition qi = table.getQPosition(
                            anchor.referencePoint,
                            anchor.values[i]);

                        if ((qPositions[i] = qi) == null) {
                            addProblem(pos, "cannot.find.q", anchor.values[i]);
                            functions[i] = Identity.IDENTITY;
                        }
                        else {
                            double qA = table.getQ(qi, anchor.referencePoint);
                            double qF = table.getQ(qi, free  .referencePoint);

                            functions[i] = Double.isNaN(qA) || Double.isNaN(qF)
                                ? Identity.IDENTITY
                                : new Linear(
                                    qA, qF,
                                    anchor.values[i], free.values[i]);

                            if (debug) {
                                logger.debug(
                                    anchor.referencePoint + ": " +
                                    qA + " -> " + functions[i].value(qA) +
                                    " / " + free.referencePoint + ": " +
                                    qF + " -> " + functions[i].value(qF));
                            }
                        }
                    } // build transforms
                } // "ungleichwertiges" interval
            } // find matching interval

            for (int i = 0; i < qPositions.length; ++i) {
                QPosition qPosition = qPositions[i];

                if (qPosition == null) {
                    continue;
                }

                if (table.interpolate(pos, out, qPosition, functions[i])) {
                    results[i].add(out[0], out[1], pos);
                }
                else {
                    addProblem(pos, "cannot.interpolate.w.q");
                }
            }
        }

        // Backjump correction
        for (int i = 0; i < results.length; ++i) {
            BackJumpCorrector bjc = new BackJumpCorrector();

            double [] ws  = results[i].getWs();
            double [] kms = results[i].getKms();

            if (bjc.doCorrection(kms, ws, this)) {
                results[i] = new WQCKms(results[i], bjc.getCorrected());
            }
        }

        // name the curves
        for (int i = 0; i < results.length; ++i) {
            results[i].setName(createName(i));
        }

        // Generate the "Umhuellende".
        results = generateInfolding(table, results, from, to, step);

        return new CalculationResult(results, this);
    }

    protected WQKms [] generateInfolding(
        WstValueTable wst,
        WQKms []      results,
        double        from,
        double        to,
        double        step
    ) {
        WstValueTable.Column [] columns = wst.getColumns();

        InfoldingColumns ic = new InfoldingColumns(columns);
        ic.markInfoldingColumns(results);

        List<WQKms> infoldings = new ArrayList<WQKms>();

        boolean [] infoldingColumns = ic.getInfoldingColumns();
        for (int i = 0; i < infoldingColumns.length; ++i) {
            if (infoldingColumns[i]) {
                continue;
            }
            double q = columns[i].getQRangeTree().findQ(from);
            if (Double.isNaN(q)) {
                addProblem(from, "cannot.find.q");
                continue;
            }
            double [] kms = DoubleUtil.explode(from, to, step);
            double [] oqs = new double[kms.length];
            double [] ows = new double[kms.length];
            boolean success =
                wst.interpolate(q, from, kms, ows, oqs, this) != null;

            if (success) {
                // TODO: generate a better name. I18N.
                String name = "Umh\u00fcllende " + columns[i].getName();
                WQKms wqkms = new WQKms(kms, oqs, ows, name);
                infoldings.add(wqkms);
            }
        }

        int N = infoldings.size();
        if (N > 0) {
            WQKms [] newResults = new WQKms[results.length + N];
            System.arraycopy(results, 0, newResults, 0, results.length);
            for (int i = 0; i < N; ++i) {
                newResults[i+results.length] = infoldings.get(i);
            }
            results = newResults;
        }

        return results;
    }

    protected String createName(int index) {
        // TODO: i18n
        StringBuilder sb = new StringBuilder(isQ ? "Q" : "W");
        sb.append(" benutzerdefiniert (");
        for (int i = 0, N = segments.size(); i < N; ++i) {
            if (i > 0) {
                sb.append("; ");
            }
            Segment segment = segments.get(i);
            sb.append((segment.backup != null
                ? segment.backup
                : segment.values)[index]);
        }
        sb.append(')');
        return sb.toString();
    }
}
// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :

http://dive4elements.wald.intevation.org