teichmann@5863: /* Copyright (C) 2011, 2012, 2013 by Bundesanstalt für Gewässerkunde
teichmann@5863:  * Software engineering by Intevation GmbH
teichmann@5863:  *
teichmann@5994:  * This file is Free Software under the GNU AGPL (>=v3)
teichmann@5863:  * and comes with ABSOLUTELY NO WARRANTY! Check out the
teichmann@5994:  * documentation coming with Dive4Elements River for details.
teichmann@5863:  */
teichmann@5863: 
teichmann@5831: package org.dive4elements.river.artifacts.model;
teichmann@4821: 
teichmann@4798: import java.util.ArrayList;
sascha@3441: import java.util.Arrays;
sascha@3441: import java.util.List;
sascha@655: 
sascha@655: import org.apache.log4j.Logger;
gernotbelger@9479: import org.dive4elements.artifacts.CallMeta;
gernotbelger@9479: import org.dive4elements.river.artifacts.access.Calculation4Access;
gernotbelger@9479: import org.dive4elements.river.artifacts.math.BackJumpCorrector;
gernotbelger@9479: import org.dive4elements.river.artifacts.math.Function;
gernotbelger@9479: import org.dive4elements.river.artifacts.math.Identity;
gernotbelger@9479: import org.dive4elements.river.artifacts.math.Linear;
gernotbelger@9479: import org.dive4elements.river.artifacts.model.WstValueTable.QPosition;
gernotbelger@9479: import org.dive4elements.river.artifacts.resources.Resources;
gernotbelger@9479: import org.dive4elements.river.model.River;
gernotbelger@9479: import org.dive4elements.river.utils.DoubleUtil;
sascha@655: 
gernotbelger@9479: public class Calculation4 extends Calculation {
teichmann@8202:     private static Logger log = Logger.getLogger(Calculation4.class);
sascha@655: 
sascha@655:     public static final double MINIMAL_STEP_WIDTH = 1e-5;
sascha@655: 
sascha@655:     protected List<Segment> segments;
sascha@655: 
teichmann@4812:     protected boolean isQ;
gernotbelger@9479:     protected double from;
gernotbelger@9479:     protected double to;
gernotbelger@9479:     protected double step;
gernotbelger@9479:     protected String river;
sascha@655: 
sascha@655:     public Calculation4() {
sascha@655:     }
sascha@655: 
gernotbelger@9479:     public Calculation4(final Calculation4Access access) {
teichmann@8202:         log.debug("Calculation4Access.cnst");
gernotbelger@9479:         final String river = access.getRiverName();
gernotbelger@9479:         final List<Segment> segments = access.getSegments();
gernotbelger@9479:         final double[] range = access.getFromToStep();
gernotbelger@9479:         final boolean isQ = access.isQ();
sascha@655: 
teichmann@4812:         if (river == null) {
teichmann@4812:             addProblem("no.river.selected");
teichmann@4812:         }
sascha@655: 
teichmann@4812:         if (range == null) {
teichmann@4812:             addProblem("no.range.found");
teichmann@4812:         }
teichmann@4812: 
teichmann@4812:         if (segments == null || segments.isEmpty()) {
teichmann@4812:             addProblem("cannot.create.segments");
teichmann@4812:         }
teichmann@4812: 
teichmann@4812:         if (!hasProblems()) {
gernotbelger@9479:             this.river = river;
teichmann@4812:             this.segments = segments;
gernotbelger@9479:             this.from = range[0];
gernotbelger@9479:             this.to = range[1];
gernotbelger@9479:             this.step = range[2];
gernotbelger@9479:             this.isQ = isQ;
teichmann@4812:         }
sascha@655:     }
sascha@655: 
gernotbelger@9479:     public CalculationResult calculate(final CallMeta meta) {
teichmann@4812:         if (hasProblems()) {
teichmann@4812:             return new CalculationResult(new WQKms[0], this);
teichmann@4812:         }
teichmann@4812: 
teichmann@4812:         WstValueTable table = null;
gernotbelger@9479:         final River r = RiverFactory.getRiver(this.river);
teichmann@4812:         if (r == null) {
teichmann@4812:             addProblem("no.river.found");
gernotbelger@9479:         } else {
teichmann@4812:             table = WstValueTableFactory.getTable(r);
teichmann@4812:             if (table == null) {
teichmann@4812:                 addProblem("no.wst.for.river");
gernotbelger@9479:             } else {
gernotbelger@9479:                 Segment.setReferencePointConvertQ(this.segments, r, this.isQ, this);
teichmann@4812:             }
teichmann@4812:         }
teichmann@4812: 
gernotbelger@9479:         return hasProblems() ? new CalculationResult(new WQKms[0], this) : innerCalculate(table, meta);
teichmann@4812:     }
teichmann@4812: 
gernotbelger@9479:     protected CalculationResult innerCalculate(final WstValueTable table, final CallMeta meta) {
gernotbelger@9479:         final boolean debug = log.isDebugEnabled();
sascha@655: 
sascha@655:         if (debug) {
gernotbelger@9479:             log.debug("calculate from " + this.from + " to " + this.to + " step " + this.step);
gernotbelger@9479:             log.debug("# segments: " + this.segments.size());
gernotbelger@9479:             for (final Segment segment : this.segments) {
teichmann@8202:                 log.debug("  " + segment);
ingo@686:             }
sascha@655:         }
sascha@655: 
gernotbelger@9479:         final int numResults = this.segments.get(0).values.length;
sascha@655: 
sascha@655:         if (numResults < 1) {
teichmann@8202:             log.debug("no values given");
sascha@2166:             addProblem("no.values.given");
sascha@709:             return new CalculationResult(new WQKms[0], this);
sascha@655:         }
sascha@655: 
gernotbelger@9479:         final WQKms[] results = new WQKms[numResults];
sascha@655:         for (int i = 0; i < results.length; ++i) {
sascha@655:             results[i] = new WQKms();
sascha@655:         }
sascha@655: 
gernotbelger@9479:         if (Math.abs(this.step) < MINIMAL_STEP_WIDTH) {
gernotbelger@9479:             this.step = MINIMAL_STEP_WIDTH;
sascha@655:         }
sascha@655: 
gernotbelger@9479:         if (this.from > this.to) {
gernotbelger@9479:             this.step = -this.step;
gernotbelger@9479:         }
sascha@655: 
gernotbelger@9479:         final QPosition[] qPositions = new QPosition[numResults];
sascha@655: 
gernotbelger@9479:         final Function[] functions = new Function[numResults];
sascha@655: 
gernotbelger@9479:         final double[] out = new double[2];
gernotbelger@9479: 
gernotbelger@9479:         final Segment sentinel = new Segment(Double.MAX_VALUE);
sascha@655:         Segment s1 = sentinel, s2 = sentinel;
sascha@655: 
gernotbelger@9479:         for (double pos = this.from; this.from < this.to ? pos <= this.to : pos >= this.to; pos = DoubleUtil.round(pos + this.step)) {
sascha@655:             if (pos < s1.referencePoint || pos > s2.referencePoint) {
sascha@655:                 if (debug) {
teichmann@8202:                     log.debug("need to find new interval for " + pos);
sascha@655:                 }
sascha@655:                 // find new interval
gernotbelger@9479:                 if (pos <= this.segments.get(0).referencePoint) {
sascha@655:                     // before first segment -> "gleichwertig"
sascha@655:                     if (debug) {
teichmann@8202:                         log.debug("before first segment -> gleichwertig");
sascha@655:                     }
gernotbelger@9479:                     final Segment first = this.segments.get(0);
gernotbelger@9479:                     final double[] values = first.values;
gernotbelger@9479:                     final double refPos = first.referencePoint;
sascha@655:                     for (int i = 0; i < qPositions.length; ++i) {
gernotbelger@9479:                         qPositions[i] = table.getQPosition(refPos, values[i]);
sascha@655:                     }
sascha@655:                     sentinel.setReferencePoint(-Double.MAX_VALUE);
sascha@655:                     s1 = sentinel;
gernotbelger@9479:                     s2 = this.segments.get(0);
sascha@655:                     Arrays.fill(functions, Identity.IDENTITY);
gernotbelger@9479:                 } else if (pos >= this.segments.get(this.segments.size() - 1).referencePoint) {
sascha@655:                     // after last segment -> "gleichwertig"
sascha@655:                     if (debug) {
teichmann@8202:                         log.debug("after last segment -> gleichwertig");
sascha@655:                     }
gernotbelger@9479:                     final Segment last = this.segments.get(this.segments.size() - 1);
gernotbelger@9479:                     final double[] values = last.values;
gernotbelger@9479:                     final double refPos = last.referencePoint;
sascha@655:                     for (int i = 0; i < qPositions.length; ++i) {
gernotbelger@9479:                         qPositions[i] = table.getQPosition(refPos, values[i]);
sascha@655:                     }
sascha@655:                     sentinel.setReferencePoint(Double.MAX_VALUE);
sascha@655:                     s1 = last;
sascha@655:                     s2 = sentinel;
sascha@655:                     Arrays.fill(functions, Identity.IDENTITY);
gernotbelger@9479:                 } else { // "ungleichwertig"
gernotbelger@9479:                          // find matching interval
sascha@655:                     if (debug) {
teichmann@8202:                         log.debug("in segments -> ungleichwertig");
sascha@655:                     }
sascha@655:                     s1 = s2 = null;
gernotbelger@9479:                     for (int i = 1, N = this.segments.size(); i < N; ++i) {
gernotbelger@9479:                         final Segment si1 = this.segments.get(i - 1);
gernotbelger@9479:                         final Segment si = this.segments.get(i);
sascha@655:                         if (debug) {
gernotbelger@9479:                             log.debug("check " + pos + " in " + si1.referencePoint + " - " + si.referencePoint);
sascha@655:                         }
gernotbelger@9479:                         if (pos >= si1.referencePoint && pos <= si.referencePoint) {
sascha@655:                             s1 = si1;
sascha@655:                             s2 = si;
sascha@655:                             break;
sascha@655:                         }
sascha@655:                     }
sascha@655: 
sascha@655:                     if (s1 == null) {
sascha@655:                         throw new IllegalStateException("no interval found");
sascha@655:                     }
sascha@655: 
sascha@655:                     Segment anchor, free;
sascha@655: 
gernotbelger@9479:                     if (this.from > this.to) {
gernotbelger@9479:                         anchor = s1;
gernotbelger@9479:                         free = s2;
gernotbelger@9479:                     } else {
gernotbelger@9479:                         anchor = s2;
gernotbelger@9479:                         free = s1;
gernotbelger@9479:                     }
sascha@655: 
sascha@655:                     // build transforms based on "gleichwertiger" phase
sascha@655:                     for (int i = 0; i < qPositions.length; ++i) {
gernotbelger@9479:                         final QPosition qi = table.getQPosition(anchor.referencePoint, anchor.values[i]);
sascha@655: 
sascha@655:                         if ((qPositions[i] = qi) == null) {
sascha@2166:                             addProblem(pos, "cannot.find.q", anchor.values[i]);
sascha@655:                             functions[i] = Identity.IDENTITY;
gernotbelger@9479:                         } else {
gernotbelger@9479:                             final double qA = table.getQ(qi, anchor.referencePoint);
gernotbelger@9479:                             final double qF = table.getQ(qi, free.referencePoint);
sascha@655: 
gernotbelger@9479:                             functions[i] = Double.isNaN(qA) || Double.isNaN(qF) ? Identity.IDENTITY : new Linear(qA, qF, anchor.values[i], free.values[i]);
ingo@686: 
ingo@686:                             if (debug) {
gernotbelger@9479:                                 log.debug(anchor.referencePoint + ": " + qA + " -> " + functions[i].value(qA) + " / " + free.referencePoint + ": " + qF + " -> "
gernotbelger@9479:                                         + functions[i].value(qF));
ingo@686:                             }
sascha@655:                         }
sascha@655:                     } // build transforms
sascha@655:                 } // "ungleichwertiges" interval
sascha@655:             } // find matching interval
sascha@655: 
sascha@655:             for (int i = 0; i < qPositions.length; ++i) {
gernotbelger@9479:                 final QPosition qPosition = qPositions[i];
sascha@655: 
sascha@655:                 if (qPosition == null) {
sascha@655:                     continue;
sascha@655:                 }
sascha@655: 
sascha@655:                 if (table.interpolate(pos, out, qPosition, functions[i])) {
sascha@655:                     results[i].add(out[0], out[1], pos);
gernotbelger@9479:                 } else {
sascha@2166:                     addProblem(pos, "cannot.interpolate.w.q");
sascha@655:                 }
sascha@655:             }
sascha@655:         }
gernotbelger@9479:         final String custom = Resources.format(meta, "common.custom");
sascha@655: 
sascha@655:         // Backjump correction
sascha@655:         for (int i = 0; i < results.length; ++i) {
gernotbelger@9479:             final BackJumpCorrector bjc = new BackJumpCorrector();
sascha@655: 
gernotbelger@9479:             final double[] ws = results[i].getWs();
gernotbelger@9479:             final double[] kms = results[i].getKms();
sascha@655: 
ingo@686:             if (bjc.doCorrection(kms, ws, this)) {
sascha@655:                 results[i] = new WQCKms(results[i], bjc.getCorrected());
sascha@655:             }
sascha@655:         }
sascha@655: 
felix@4967:         // Name the curves.
sascha@655:         for (int i = 0; i < results.length; ++i) {
gernotbelger@9479:             results[i].setName(createName(i, custom));
sascha@655:         }
sascha@655: 
teichmann@4798:         // Generate the "Umhuellende".
gernotbelger@9479:         final ConstantWQKms[] infoldings = generateInfolding(table, results, this.from, this.to, this.step);
teichmann@4821: 
teichmann@4821:         // TODO: Use qkms in a new result type.
gernotbelger@9479:         final WQKms[] newResults = new WQKms[results.length + infoldings.length];
gernotbelger@9479:         System.arraycopy(results, 0, newResults, 0, results.length);
gernotbelger@9479:         System.arraycopy(infoldings, 0, newResults, results.length, infoldings.length);
teichmann@4798: 
teichmann@4835:         return new CalculationResult(newResults, this);
sascha@655:     }
sascha@655: 
gernotbelger@9479:     protected ConstantWQKms[] generateInfolding(final WstValueTable wst, final WQKms[] results, final double from, final double to, final double step) {
gernotbelger@9479:         final WstValueTable.Column[] columns = wst.getColumns();
teichmann@4798: 
gernotbelger@9479:         final InfoldingColumns ic = new InfoldingColumns(columns);
teichmann@4798:         ic.markInfoldingColumns(results);
teichmann@4798: 
gernotbelger@9479:         final List<ConstantWQKms> infoldings = new ArrayList<>();
teichmann@4798: 
gernotbelger@9479:         final boolean[] infoldingColumns = ic.getInfoldingColumns();
teichmann@4821: 
gernotbelger@9479:         double[] kms = null;
gernotbelger@9479:         double[] ws = null;
teichmann@4821: 
teichmann@4798:         for (int i = 0; i < infoldingColumns.length; ++i) {
teichmann@4799:             if (!infoldingColumns[i]) {
teichmann@4798:                 continue;
teichmann@4798:             }
teichmann@4821: 
teichmann@4835:             if (kms == null) {
teichmann@4835:                 kms = DoubleUtil.explode(from, to, step);
gernotbelger@9479:                 ws = new double[kms.length];
teichmann@4835:             }
teichmann@4821: 
gernotbelger@9479:             final QRangeTree.QuickQFinder qf = columns[i].getQRangeTree().new QuickQFinder();
teichmann@4821: 
gernotbelger@9479:             final int numProblemsBefore = numProblems();
gernotbelger@9479:             final double[] qs = qf.findQs(kms, this);
teichmann@4821: 
gernotbelger@9479:             final String name = columns[i].getName();
gernotbelger@9479:             final ConstantWQKms infolding = new ConstantWQKms(kms, qs, ws, name);
teichmann@4821: 
teichmann@4821:             if (numProblems() > numProblemsBefore) {
teichmann@4835:                 infolding.removeNaNs();
teichmann@4798:             }
teichmann@4798: 
teichmann@4835:             infoldings.add(infolding);
teichmann@4798:         }
teichmann@4798: 
teichmann@7254:         for (int i = 0, I = infoldings.size(); i < I; i++) {
gernotbelger@9479:             final ConstantWQKms infolding = infoldings.get(i);
gernotbelger@9479:             final String name = infolding.getName();
rrenkert@5138:             // TODO: i18n
rrenkert@5138:             if (i == 0) {
teichmann@7254:                 infolding.setName("untere Umh\u00fcllende " + name);
gernotbelger@9479:             } else if (i == I - 1) {
teichmann@7254:                 infolding.setName("obere Umh\u00fcllende " + name);
gernotbelger@9479:             } else {
teichmann@7254:                 infolding.setName("geschnitten " + name);
rrenkert@5138:             }
rrenkert@5138:         }
teichmann@7254: 
teichmann@4835:         return infoldings.toArray(new ConstantWQKms[infoldings.size()]);
teichmann@4798:     }
teichmann@4798: 
felix@4967:     // TODO: issue1109/2, merge with FixRealizingCalculation
gernotbelger@9479:     protected String createName(final int index, final String custom) {
gernotbelger@9479: 
gernotbelger@9479:         final StringBuilder sb = new StringBuilder(this.isQ ? "Q" : "W");
gernotbelger@9479:         sb.append(" ").append(custom).append(" (");
gernotbelger@9479:         for (int i = 0, N = this.segments.size(); i < N; ++i) {
sascha@655:             if (i > 0) {
sascha@655:                 sb.append("; ");
sascha@655:             }
gernotbelger@9479:             final Segment segment = this.segments.get(i);
gernotbelger@9479:             sb.append((segment.backup != null ? segment.backup : segment.values)[index]);
sascha@655:         }
sascha@655:         sb.append(')');
sascha@655:         return sb.toString();
sascha@655:     }
sascha@655: }