Java Code Examples for org.apache.commons.math3.util.FastMath#ceil()

/**
 * Calculates the exact value of {@code P(D_n < d)} using method described
 * in [1] and {@link org.apache.commons.math3.fraction.BigFraction} (see
 * above).
 *
 * @param d statistic
 * @return the two-sided probability of {@code P(D_n < d)}
 * @throws MathArithmeticException if algorithm fails to convert {@code h}
 * to a {@link org.apache.commons.math3.fraction.BigFraction} in expressing
 * {@code d} as {@code (k - h) / m} for integer {@code k, m} and
 * {@code 0 <= h < 1}.
 */
private double exactK(double d) throws MathArithmeticException {

    final int k = (int) FastMath.ceil(n * d);

    final FieldMatrix<BigFraction> H = this.createH(d);
    final FieldMatrix<BigFraction> Hpower = H.power(n);

    BigFraction pFrac = Hpower.getEntry(k - 1, k - 1);

    for (int i = 1; i <= n; ++i) {
        pFrac = pFrac.multiply(i).divide(n);
    }

    /*
     * BigFraction.doubleValue converts numerator to double and the
     * denominator to double and divides afterwards. That gives NaN quite
     * easy. This does not (scale is the number of digits):
     */
    return pFrac.bigDecimalValue(20, BigDecimal.ROUND_HALF_UP).doubleValue();
}
 

/**
 * Calculates the exact value of {@code P(D_n < d)} using method described
 * in [1] and {@link org.apache.commons.math3.fraction.BigFraction} (see
 * above).
 *
 * @param d statistic
 * @return the two-sided probability of {@code P(D_n < d)}
 * @throws MathArithmeticException if algorithm fails to convert {@code h}
 * to a {@link org.apache.commons.math3.fraction.BigFraction} in expressing
 * {@code d} as {@code (k - h) / m} for integer {@code k, m} and
 * {@code 0 <= h < 1}.
 */
private double exactK(double d) throws MathArithmeticException {

    final int k = (int) FastMath.ceil(n * d);

    final FieldMatrix<BigFraction> H = this.createH(d);
    final FieldMatrix<BigFraction> Hpower = H.power(n);

    BigFraction pFrac = Hpower.getEntry(k - 1, k - 1);

    for (int i = 1; i <= n; ++i) {
        pFrac = pFrac.multiply(i).divide(n);
    }

    /*
     * BigFraction.doubleValue converts numerator to double and the
     * denominator to double and divides afterwards. That gives NaN quite
     * easy. This does not (scale is the number of digits):
     */
    return pFrac.bigDecimalValue(20, BigDecimal.ROUND_HALF_UP).doubleValue();
}
 

public double cdf(double mean, double std, double x) {
    double zscore = (x - mean) / std;

    if (zscore > MAXZSCORE) {
        return 1.0;
    }
    if (zscore < MINZSCORE) {
        return 0.0;
    }

    // Interpolate the CDF
    double exactEntry = zscore / GRANULARITY + LUT_OFFSET;
    int lowerEntry = (int) FastMath.floor(exactEntry);
    int higherEntry = (int) FastMath.ceil(exactEntry);

    if (lowerEntry == higherEntry) {
        return CDF_LUT[lowerEntry];
    } else {
        return CDF_LUT[lowerEntry] +
                ((exactEntry - lowerEntry) * (CDF_LUT[higherEntry] - CDF_LUT[lowerEntry]));
    }
}
 

/**
 * Calculates {@code P(D_n < d)} using method described in [1] and doubles
 * (see above).
 *
 * @param d statistic
 * @return the two-sided probability of {@code P(D_n < d)}
 * @throws MathArithmeticException if algorithm fails to convert {@code h}
 * to a {@link org.apache.commons.math3.fraction.BigFraction} in expressing
 * {@code d} as {@code (k - h) / m} for integer {@code k, m} and
 * {@code 0 <= h < 1}.
 */
private double roundedK(double d) throws MathArithmeticException {

    final int k = (int) FastMath.ceil(n * d);
    final FieldMatrix<BigFraction> HBigFraction = this.createH(d);
    final int m = HBigFraction.getRowDimension();

    /*
     * Here the rounding part comes into play: use
     * RealMatrix instead of FieldMatrix<BigFraction>
     */
    final RealMatrix H = new Array2DRowRealMatrix(m, m);

    for (int i = 0; i < m; ++i) {
        for (int j = 0; j < m; ++j) {
            H.setEntry(i, j, HBigFraction.getEntry(i, j).doubleValue());
        }
    }

    final RealMatrix Hpower = H.power(n);

    double pFrac = Hpower.getEntry(k - 1, k - 1);

    for (int i = 1; i <= n; ++i) {
        pFrac *= (double) i / (double) n;
    }

    return pFrac;
}
 

/**
 * {@inheritDoc}This method in particular for R_1 uses ceil(pos-0.5)
 */
@Override
protected double estimate(final double[] values,
                          final int[] pivotsHeap, final double pos,
                          final int length, final KthSelector kthSelector) {
    return super.estimate(values, pivotsHeap, FastMath.ceil(pos - 0.5), length, kthSelector);
}
 

/**
 * Calculates {@code P(D_n < d)} using method described in [1] and doubles
 * (see above).
 *
 * @param d statistic
 * @return the two-sided probability of {@code P(D_n < d)}
 * @throws MathArithmeticException if algorithm fails to convert {@code h}
 * to a {@link org.apache.commons.math3.fraction.BigFraction} in expressing
 * {@code d} as {@code (k - h) / m} for integer {@code k, m} and
 * {@code 0 <= h < 1}.
 */
private double roundedK(double d) throws MathArithmeticException {

    final int k = (int) FastMath.ceil(n * d);
    final FieldMatrix<BigFraction> HBigFraction = this.createH(d);
    final int m = HBigFraction.getRowDimension();

    /*
     * Here the rounding part comes into play: use
     * RealMatrix instead of FieldMatrix<BigFraction>
     */
    final RealMatrix H = new Array2DRowRealMatrix(m, m);

    for (int i = 0; i < m; ++i) {
        for (int j = 0; j < m; ++j) {
            H.setEntry(i, j, HBigFraction.getEntry(i, j).doubleValue());
        }
    }

    final RealMatrix Hpower = H.power(n);

    double pFrac = Hpower.getEntry(k - 1, k - 1);

    for (int i = 1; i <= n; ++i) {
        pFrac *= (double) i / (double) n;
    }

    return pFrac;
}
 

/**
 * {@inheritDoc}This method in particular for R_1 uses ceil(pos-0.5)
 */
@Override
protected double estimate(final double[] values,
                          final int[] pivotsHeap, final double pos,
                          final int length, final KthSelector kthSelector) {
    return super.estimate(values, pivotsHeap, FastMath.ceil(pos - 0.5), length, kthSelector);
}
 

/** Get the smallest whole number larger than instance.
 * @return ceil(this)
 */
public DerivativeStructure ceil() {
    return new DerivativeStructure(compiler.getFreeParameters(),
                                   compiler.getOrder(),
                                   FastMath.ceil(data[0]));
}
 

/**
 * Helper method to create a multivariate normal mixture model which can be
 * used to initialize {@link #fit(MixtureMultivariateRealDistribution)}.
 *
 * This method uses the data supplied to the constructor to try to determine
 * a good mixture model at which to start the fit, but it is not guaranteed
 * to supply a model which will find the optimal solution or even converge.
 *
 * @param data Data to estimate distribution
 * @param numComponents Number of components for estimated mixture
 * @return Multivariate normal mixture model estimated from the data
 * @throws NumberIsTooLargeException if {@code numComponents\ is greater
 * than the number of data rows.
 * @throws NumberIsTooSmallException if {@code numComponents < 2}.
 * @throws NotStrictlyPositiveException if data has less than 2 rows
 * @throws DimensionMismatchException if rows of data have different numbers
 *             of columns
 * @see #fit
 */
public static MixtureMultivariateNormalDistribution estimate(final double[][] data,
                                                             final int numComponents)
    throws NotStrictlyPositiveException,
           DimensionMismatchException {
    if (data.length < 2) {
        throw new NotStrictlyPositiveException(data.length);
    }
    if (numComponents < 2) {
        throw new NumberIsTooSmallException(numComponents, 2, true);
    }
    if (numComponents > data.length) {
        throw new NumberIsTooLargeException(numComponents, data.length, true);
    }

    final int numRows = data.length;
    final int numCols = data[0].length;

    // sort the data
    final DataRow[] sortedData = new DataRow[numRows];
    for (int i = 0; i < numRows; i++) {
        sortedData[i] = new DataRow(data[i]);
    }
    Arrays.sort(sortedData);

    final int totalBins = numComponents;

    // uniform weight for each bin
    final double weight = 1d / totalBins;

    // components of mixture model to be created
    final List<Pair<Double, MultivariateNormalDistribution>> components =
            new ArrayList<Pair<Double, MultivariateNormalDistribution>>();

    // create a component based on data in each bin
    for (int binNumber = 1; binNumber <= totalBins; binNumber++) {
        // minimum index from sorted data for this bin
        final int minIndex
            = (int) FastMath.max(0,
                                 FastMath.floor((binNumber - 1) * numRows / totalBins));

        // maximum index from sorted data for this bin
        final int maxIndex
            = (int) FastMath.ceil(binNumber * numRows / numComponents) - 1;

        // number of data records that will be in this bin
        final int numBinRows = maxIndex - minIndex + 1;

        // data for this bin
        final double[][] binData = new double[numBinRows][numCols];

        // mean of each column for the data in the this bin
        final double[] columnMeans = new double[numCols];

        // populate bin and create component
        for (int i = minIndex, iBin = 0; i <= maxIndex; i++, iBin++) {
            for (int j = 0; j < numCols; j++) {
                final double val = sortedData[i].getRow()[j];
                columnMeans[j] += val;
                binData[iBin][j] = val;
            }
        }

        MathArrays.scaleInPlace(1d / numBinRows, columnMeans);

        // covariance matrix for this bin
        final double[][] covMat
            = new Covariance(binData).getCovarianceMatrix().getData();
        final MultivariateNormalDistribution mvn
            = new MultivariateNormalDistribution(columnMeans, covMat);

        components.add(new Pair<Double, MultivariateNormalDistribution>(weight, mvn));
    }

    return new MixtureMultivariateNormalDistribution(components);
}
 

/**
 * {@inheritDoc}
 *
 * The default implementation returns
 * <ul>
 * <li>{@link #getSupportLowerBound()} for {@code p = 0},</li>
 * <li>{@link #getSupportUpperBound()} for {@code p = 1}, and</li>
 * <li>{@link #solveInverseCumulativeProbability(double, int, int)} for
 *     {@code 0 < p < 1}.</li>
 * </ul>
 */
public int inverseCumulativeProbability(final double p) throws OutOfRangeException {
    if (p < 0.0 || p > 1.0) {
        throw new OutOfRangeException(p, 0, 1);
    }

    int lower = getSupportLowerBound();
    if (p == 0.0) {
        return lower;
    }
    if (lower == Integer.MIN_VALUE) {
        if (checkedCumulativeProbability(lower) >= p) {
            return lower;
        }
    } else {
        lower -= 1; // this ensures cumulativeProbability(lower) < p, which
                    // is important for the solving step
    }

    int upper = getSupportUpperBound();
    if (p == 1.0) {
        return upper;
    }

    // use the one-sided Chebyshev inequality to narrow the bracket
    // cf. AbstractRealDistribution.inverseCumulativeProbability(double)
    final double mu = getNumericalMean();
    final double sigma = FastMath.sqrt(getNumericalVariance());
    final boolean chebyshevApplies = !(Double.isInfinite(mu) || Double.isNaN(mu) ||
            Double.isInfinite(sigma) || Double.isNaN(sigma) || sigma == 0.0);
    if (chebyshevApplies) {
        double k = FastMath.sqrt((1.0 - p) / p);
        double tmp = mu - k * sigma;
        if (tmp > lower) {
            lower = ((int) FastMath.ceil(tmp)) - 1;
        }
        k = 1.0 / k;
        tmp = mu + k * sigma;
        if (tmp < upper) {
            upper = ((int) FastMath.ceil(tmp)) - 1;
        }
    }

    return solveInverseCumulativeProbability(p, lower, upper);
}
 

public static double[] vectCeilWrite(double[] a, int[] aix, int ai, int alen, int len) {
	double[] c = allocVector(len, true);
	for( int j = ai; j < ai+alen; j++ )
		c[aix[j]] = FastMath.ceil(a[j]);
	return c;
}
 

/** {@inheritDoc} */
public double value(double x) {
    return FastMath.ceil(x);
}
 

public static void vectCeilAdd(double[] a, double[] c, int[] aix, int ai, int ci, int alen, int len) {
	for( int j = ai; j < ai+alen; j++ )
		c[ci + aix[j]] += FastMath.ceil(a[j]);
}
 

public static void vectCeilAdd(double[] a, double[] c, int ai, int ci, int len) {
	for( int j = ai; j < ai+len; j++, ci++)
		c[ci] +=  FastMath.ceil(a[j]);
}
 

@Override
public double execute (double in) {
	switch(bFunc) {
		case SIN:    return FASTMATH ? FastMath.sin(in) : Math.sin(in);
		case COS:    return FASTMATH ? FastMath.cos(in) : Math.cos(in);
		case TAN:    return FASTMATH ? FastMath.tan(in) : Math.tan(in);
		case ASIN:   return FASTMATH ? FastMath.asin(in) : Math.asin(in);
		case ACOS:   return FASTMATH ? FastMath.acos(in) : Math.acos(in);
		case ATAN:   return Math.atan(in); //faster in Math
		// FastMath.*h is faster 98% of time than Math.*h in initial micro-benchmarks
		case SINH:   return FASTMATH ? FastMath.sinh(in) : Math.sinh(in);
		case COSH:   return FASTMATH ? FastMath.cosh(in) : Math.cosh(in);
		case TANH:   return FASTMATH ? FastMath.tanh(in) : Math.tanh(in);
		case CEIL:   return FASTMATH ? FastMath.ceil(in) : Math.ceil(in);
		case FLOOR:  return FASTMATH ? FastMath.floor(in) : Math.floor(in);
		case LOG:    return Math.log(in); //faster in Math
		case LOG_NZ: return (in==0) ? 0 : Math.log(in); //faster in Math
		case ABS:    return Math.abs(in); //no need for FastMath
		case SIGN:   return FASTMATH ? FastMath.signum(in) : Math.signum(in);
		case SQRT:   return Math.sqrt(in); //faster in Math
		case EXP:    return FASTMATH ? FastMath.exp(in) : Math.exp(in);
		case ROUND: return Math.round(in); //no need for FastMath
		
		case PLOGP:
			if (in == 0.0)
				return 0.0;
			else if (in < 0)
				return Double.NaN;
			else //faster in Math
				return in * Math.log(in);
		
		case SPROP:
			//sample proportion: P*(1-P)
			return in * (1 - in); 

		case SIGMOID:
			//sigmoid: 1/(1+exp(-x))
			return FASTMATH ? 1 / (1 + FastMath.exp(-in))  : 1 / (1 + Math.exp(-in));
		
		case ISNA: return Double.isNaN(in) ? 1 : 0;
		case ISNAN: return Double.isNaN(in) ? 1 : 0;
		case ISINF: return Double.isInfinite(in) ? 1 : 0;
		
		default:
			throw new DMLRuntimeException("Builtin.execute(): Unknown operation: " + bFunc);
	}
}
 

/**
 * Procedure rounds the given value to the given number of significant digits see
 * http://stackoverflow.com/questions/202302
 *
 * Note: The maximum double value in Java is on the order of 10^308, while the minimum value is on the order of
 * 10^-324. Therefore, you can run into trouble when applying the function roundToSignificantFigures to something
 * that's within a few powers of ten of Double.MIN_VALUE.
 *
 * Consequently, the variable magnitude may become Infinity, and it's all garbage from then on out. Fortunately,
 * this is not an insurmountable problem: it is only the factor magnitude that's overflowing. What really matters is
 * the product num * magnitude, and that does not overflow. One way of resolving this is by breaking up the
 * multiplication by the factor magintude into two steps.
 */
public static double roundToNumberOfSignificantDigits(double num, int n)
{
    if (num == 0)
    {
        return 0;
    }

    try
    {
        return new BigDecimal(num).round(new MathContext(n, RoundingMode.HALF_EVEN)).doubleValue();
    }
    catch (ArithmeticException ex)
    {
        // nothing to do
    }

    final double maxPowerOfTen = FastMath.floor(FastMath.log10(Double.MAX_VALUE));
    final double d = FastMath.ceil(FastMath.log10(num < 0 ? -num : num));
    final int power = n - (int) d;

    double firstMagnitudeFactor = 1.0;
    double secondMagnitudeFactor = 1.0;
    if (power > maxPowerOfTen)
    {
        firstMagnitudeFactor = FastMath.pow(10.0, maxPowerOfTen);
        secondMagnitudeFactor = FastMath.pow(10.0, (double) power - maxPowerOfTen);
    }
    else
    {
        firstMagnitudeFactor = FastMath.pow(10.0, (double) power);
    }

    double toBeRounded = num * firstMagnitudeFactor;
    toBeRounded *= secondMagnitudeFactor;

    final long shifted = FastMath.round(toBeRounded);
    double rounded = ((double) shifted) / firstMagnitudeFactor;
    rounded /= secondMagnitudeFactor;
    return rounded;
}
 

public static double[] vectCeilWrite(double[] a, int[] aix, int ai, int alen, int len) {
	double[] c = allocVector(len, true);
	for( int j = ai; j < ai+alen; j++ )
		c[aix[j]] = FastMath.ceil(a[j]);
	return c;
}
 

public static double[] vectCeilWrite(double[] a, int ai, int len) {
	double[] c = allocVector(len, false);
	for( int j = 0; j < len; j++, ai++)
		c[j] = FastMath.ceil(a[ai]);
	return c;
}
 

public static void vectCeilAdd(double[] a, double[] c, int[] aix, int ai, int ci, int alen, int len) {
	for( int j = ai; j < ai+alen; j++ )
		c[ci + aix[j]] += FastMath.ceil(a[j]);
}
 

/** {@inheritDoc} */
public double value(double x) {
    return FastMath.ceil(x);
}