view artifacts/src/main/java/org/dive4elements/river/utils/DoubleUtil.java @ 8587:07c9ac22f611

(issue1755) Generalise BedQuality result handling The bedquality calculation now produces a result for each time period which has BedQualityResultValues for each specific result type. Formally this was split up in density, porosity and diameter classes with some bedload diameter classes mixed in for extra fun. The intent of this commit is to allow more shared code and generic access patterns to the BedQuality results.
author Andre Heinecke <andre.heinecke@intevation.de>
date Wed, 18 Mar 2015 18:42:08 +0100
parents 073ea4bcea58
children b486812f4f14
line wrap: on
line source
/* Copyright (C) 2011, 2012, 2013 by Bundesanstalt für Gewässerkunde
 * Software engineering by Intevation GmbH
 *
 * This file is Free Software under the GNU AGPL (>=v3)
 * and comes with ABSOLUTELY NO WARRANTY! Check out the
 * documentation coming with Dive4Elements River for details.
 */

package org.dive4elements.river.utils;

import org.dive4elements.river.artifacts.math.Linear;

import org.apache.commons.math.analysis.interpolation.LinearInterpolator;
import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;

import org.apache.commons.math.exception.DimensionMismatchException;
import org.apache.commons.math.exception.NumberIsTooSmallException;
import org.apache.commons.math.exception.NonMonotonousSequenceException;

import gnu.trove.TDoubleArrayList;

import java.util.Arrays;
import java.util.Comparator;

import org.apache.log4j.Logger;

/** Utils to deal with Double precision values. */
public class DoubleUtil
{
    /** Private log. */
    private static Logger log = Logger.getLogger(DoubleUtil.class);

    public static final double DEFAULT_STEP_PRECISION = 1e6;

    public static final Comparator<double []> DOUBLE_PAIR_CMP =
        new Comparator<double []>() {
            @Override
            public int compare(double [] a, double [] b) {
                double diff = a[0] - b[0];
                if (diff < 0d) return -1;
                if (diff > 0d) return +1;
                return 0;
            }
        };

    /** EPSILON for comparison of double precision values. */
    public static final double EPSILON = 1e-4;

    private DoubleUtil() {
    }

    public static final double [] explode(double from, double to, double step) {
        return explode(from, to, step, DEFAULT_STEP_PRECISION);
    }

    public static final double round(double x, double precision) {
        return Math.round(x * precision)/precision;
    }

    public static final double round(double x) {
        return Math.round(x * DEFAULT_STEP_PRECISION)/DEFAULT_STEP_PRECISION;
    }

    /**
     * Returns array with values from parameter from to to with given step width.
     * from and to are included.
     */
    public static final double [] explode(
        double from,
        double to,
        double step,
        double precision
    ) {
        double lower = from;

        double diff = to - from;
        double tmp  = diff / step;
        int    num = (int)Math.abs(Math.ceil(tmp)) + 1;

        if (num < 1) {
            return new double[0];
        }

        double [] values = new double[num];

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

        double max = Math.max(from, to);

        for (int idx = 0; idx < num; idx++) {
            if (lower - max > EPSILON) {
                return Arrays.copyOfRange(values, 0, idx);
            }

            values[idx] = round(lower, precision);
            lower      += step;
        }

        return values;
    }

    public static final double interpolateSorted(
        double [] xs,
        double [] ys,
        double x
    ) {
        int lo = 0, hi = xs.length-1;

        int mid = -1;

        while (lo <= hi) {
            mid = (lo + hi) >> 1;
            double mx = xs[mid];
                 if (x < mx) hi = mid - 1;
            else if (x > mx) lo = mid + 1;
            else return ys[mid];
        }
        if (mid < lo) {
            return lo >= xs.length
                ? Double.NaN
                : Linear.linear(x, xs[mid], xs[mid+1], ys[mid], ys[mid+1]);
        }
        return hi < 0
            ? Double.NaN
            : Linear.linear(x, xs[mid-1], xs[mid], ys[mid-1], ys[mid]);
    }

    public static final boolean isIncreasing(double [] array) {
        int inc = 0;
        int dec = 0;
        int sweet = (array.length-1)/2;
        for (int i = 1; i < array.length; ++i) {
            if (array[i] > array[i-1]) {
                if (++inc > sweet) {
                    return true;
                }
            }
            else if (++dec > sweet) {
                return false;
            }
        }
        return inc > sweet;
    }

    public static final double [] swap(double [] array) {
        int lo = 0;
        int hi = array.length-1;
        while (hi > lo) {
            double t  = array[lo];
            array[lo] = array[hi];
            array[hi] = t;
            ++lo;
            --hi;
        }

        return array;
    }

    public static final double [] swapClone(double [] in) {
        double [] out = new double[in.length];

        for (int j = out.length-1, i = 0; j >= 0;) {
            out[j--] = in[i++];
        }

        return out;
    }

    public static final double [] sumDiffs(double [] in) {
        double [] out = new double[in.length];

        for (int i = 1; i < out.length; ++i) {
            out[i] = out[i-1] + Math.abs(in[i-1] - in[i]);
        }

        return out;
    }

    public static final double sum(double [] values) {
        double sum = 0.0;
        for (double value: values) {
            sum += value;
        }
        return sum;
    }

    public static final double [] fill(int N, double value) {
        double [] result = new double[N];
        Arrays.fill(result, value);
        return result;
    }


    /** Use with parseSegments. */
    public interface SegmentCallback {
        void newSegment(double from, double to, double [] values);
    }


    /** Call callback for every string split by colon.
     * Expected format FROM:TO:VALUE1,VALUE2,VALUE3*/
    public static final void parseSegments(
        String          input,
        SegmentCallback callback
    ) {
        TDoubleArrayList vs = new TDoubleArrayList();

        for (String segmentStr: input.split(":")) {
            String [] parts = segmentStr.split(";");
            if (parts.length < 3) {
                log.warn("invalid segment: '" + segmentStr + "'");
                continue;
            }
            try {
                double from = Double.parseDouble(parts[0].trim());
                double to   = Double.parseDouble(parts[1].trim());

                vs.resetQuick();

                for (String valueStr: parts[3].split(",")) {
                    vs.add(round(Double.parseDouble(valueStr.trim())));
                }

                callback.newSegment(from, to, vs.toNativeArray());
            }
            catch (NumberFormatException nfe) {
                log.warn("invalid segment: '" + segmentStr + "'");
            }
        }
    }

    public static final boolean isValid(double [][] data) {
        for (double [] ds: data) {
            if (!isValid(ds)) {
                return false;
            }
        }
        return true;
    }

    public static final boolean isValid(double [] data) {
        for (double d: data) {
            if (Double.isNaN(d)) {
                return false;
            }
        }
        return true;
    }

    public static final boolean isNaN(double [] values) {
        for (double value: values) {
            if (!Double.isNaN(value)) {
                return false;
            }
        }
        return true;
    }

    /** In an array of doubles, search and return the maximum value. */
    public static final double maxInArray(double[] values) {
        double max = - Double.MAX_VALUE;
        for (double d: values) {
            if (d > max) max = d;
        }
        return max;
    }



    /** Sort a and b with a as key. b is ordered accordingly */
    public static final void sortByFirst(double [] a, double [] b) {
        // XXX: Not efficient but bulletproof.
        double [][] pairs = new double[a.length][2];
        for (int i = 0; i < a.length; ++i) {
            double [] p = pairs[i];
            p[0] = a[i];
            p[1] = b[i];

        }
        Arrays.sort(pairs, DOUBLE_PAIR_CMP);
        for (int i = 0; i < a.length; ++i) {
            double [] p = pairs[i];
            a[i] = p[0];
            b[i] = p[1];
        }
    }

    public static void removeNaNs(TDoubleArrayList [] arrays) {

        int dest = 0;

        int A = arrays.length;
        int N = arrays[0].size();

        OUTER: for (int i = 0; i < N; ++i) {
            for (int j = 0; j < A; ++j) {
                TDoubleArrayList a = arrays[j];
                double v = a.getQuick(i);
                if (Double.isNaN(v)) {
                    continue OUTER;
                }
                a.setQuick(dest, v);
            }
            ++dest;
        }

        if (dest < N) {
            for (int i = 0; i < A; ++i) {
                arrays[i].remove(dest, N-dest);
            }
        }
    }

    /** Convieniance function for results to get an interpolator.
     * This is basically a static wrapper to for LinearInterpolator.interpolate
     * with error handling. Returns null on error.*/
    public static PolynomialSplineFunction getLinearInterpolator(double[] x, double[] y) {
        LinearInterpolator lpol = new LinearInterpolator();
        try {
            return lpol.interpolate(x, y);
        } catch (DimensionMismatchException e) {
            log.error("KMs and Result values have different sizes. Failed to interpolate: " +
                    e.getMessage());
        } catch (NonMonotonousSequenceException e) {
            log.error("KMs are not monotonous. Failed to interpolate: " + e.getMessage());
        } catch (NumberIsTooSmallException e) {
            log.error("Result is to small. Failed to interpolate: " + e.getMessage());
        }
        return null;
    }

    public static PolynomialSplineFunction getLinearInterpolator(TDoubleArrayList x, TDoubleArrayList y) {
        return getLinearInterpolator(x.toNativeArray(), y.toNativeArray());
    }
}
// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :

http://dive4elements.wald.intevation.org