view flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/Calculation4.java @ 3304:02d5731b43a2

SQ relation: Build a valid SQResult. flys-artifacts/trunk@4989 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Fri, 13 Jul 2012 11:12:35 +0000
parents a0d9a99a5d17
children cb11919cccf9
line wrap: on
line source
package de.intevation.flys.artifacts.model;

import java.util.List;
import java.util.Arrays;
import java.util.Comparator;
import java.util.Collections;

import de.intevation.flys.utils.DoubleUtil;

import de.intevation.flys.model.River;
import de.intevation.flys.model.Gauge;
import de.intevation.flys.model.DischargeTable;

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

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

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;

    public static final Comparator<Segment> REF_CMP =
        new Comparator<Segment>() {
            public int compare(Segment a, Segment b) {
                double d = a.referencePoint - b.referencePoint;
                if (d < 0d) return -1;
                return d > 0d ? +1 : 0;
            }
        };

    protected List<Segment> segments;

    protected boolean       isQ;

    public Calculation4() {
    }

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

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

        int numResults = -1;

        // assign reference points
        for (Segment segment: segments) {
            Gauge gauge = river.maxOverlap(segment.getFrom(), segment.getTo());

            if (gauge == null) {
                logger.warn("no gauge found. Defaults to mid point.");
                segment.setReferencePoint(
                    0.5*(segment.getFrom()+segment.getTo()));
            }
            else {
                double ref = gauge.getStation().doubleValue();
                logger.debug(
                    "reference gauge: " + gauge.getName() +
                    " (km " + ref + ")");
                segment.setReferencePoint(ref);
            }

            double [] values = segment.values;

            if (numResults == -1) {
                numResults = values.length;
            }
            else if (numResults != values.length) {
                throw new IllegalArgumentException("wrong length of values");
            }

            // convert to Q if needed
            if (!isQ && gauge != null) {

                DischargeTable dt = gauge.fetchMasterDischargeTable();

                double [][] table =
                    DischargeTables.loadDischargeTableValues(dt, 1);

                // need the original values for naming
                segment.backup();

                for (int i = 0; i < values.length; ++i) {
                    double w = values[i] / 100.0;
                    double [] qs = DischargeTables.getQsForW(table, w);
                    if (qs.length == 0) {
                        logger.warn("No Qs found for W = " + values[i]);
                        addProblem("cannot.find.w.for.q", values[i]);
                        values[i] = Double.NaN;
                    }
                    else {
                        values[i] = qs[0];
                        if (qs.length > 1) {
                            logger.warn(
                                "More than one Q found for W = " + values[i]);
                        }
                    }
                }
            }
        } // for all segments

        Collections.sort(segments, REF_CMP);
    }

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

        return new CalculationResult(results, this);
    }

    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