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

@Test
public void testSumDeltaMinusDeltaSum() {

    final int ulps = 3;
    for (int i = 0; i < SUM_DELTA_MINUS_DELTA_SUM_REF.length; i++) {
        final double[] ref = SUM_DELTA_MINUS_DELTA_SUM_REF[i];
        final double a = ref[0];
        final double b = ref[1];
        final double expected = ref[2];
        final double actual = sumDeltaMinusDeltaSum(a, b);
        final double tol = ulps * FastMath.ulp(expected);
        final StringBuilder builder = new StringBuilder();
        builder.append(a).append(", ").append(b);
        Assert.assertEquals(builder.toString(), expected, actual, tol);
    }
}
 

@Test
public void testLogGamma() {
    final int ulps = 3;
    for (int i = 0; i < LOG_GAMMA_REF.length; i++) {
        final double[] data = LOG_GAMMA_REF[i];
        final double x = data[0];
        final double expected = data[1];
        final double actual = Gamma.logGamma(x);
        final double tol;
        if (expected == 0.0) {
            tol = 1E-15;
        } else {
            tol = ulps * FastMath.ulp(expected);
        }
        Assert.assertEquals(Double.toString(x), expected, actual, tol);
    }
}
 

/**
 * Returns true iff there is a column that is a scalar multiple of column
 * in searchMatrix (modulo tolerance)
 */
private boolean isIncludedColumn(double[] column, RealMatrix searchMatrix,
        double tolerance) {
    boolean found = false;
    int i = 0;
    while (!found && i < searchMatrix.getColumnDimension()) {
        double multiplier = 1.0;
        boolean matching = true;
        int j = 0;
        while (matching && j < searchMatrix.getRowDimension()) {
            double colEntry = searchMatrix.getEntry(j, i);
            // Use the first entry where both are non-zero as scalar
            if (FastMath.abs(multiplier - 1.0) <= FastMath.ulp(1.0) && FastMath.abs(colEntry) > 1E-14
                    && FastMath.abs(column[j]) > 1e-14) {
                multiplier = colEntry / column[j];
            }
            if (FastMath.abs(column[j] * multiplier - colEntry) > tolerance) {
                matching = false;
            }
            j++;
        }
        found = matching;
        i++;
    }
    return found;
}
 

/** Check the integration span.
 * @param equations set of differential equations
 * @param t target time for the integration
 * @exception NumberIsTooSmallException if integration span is too small
 * @exception DimensionMismatchException if adaptive step size integrators
 * tolerance arrays dimensions are not compatible with equations settings
 */
protected void sanityChecks(final ExpandableStatefulODE equations, final double t)
    throws NumberIsTooSmallException, DimensionMismatchException {

    final double threshold = 1000 * FastMath.ulp(FastMath.max(FastMath.abs(equations.getTime()),
                                                              FastMath.abs(t)));
    final double dt = FastMath.abs(equations.getTime() - t);
    if (dt <= threshold) {
        throw new NumberIsTooSmallException(LocalizedFormats.TOO_SMALL_INTEGRATION_INTERVAL,
                                            dt, threshold, false);
    }

}
 

@Test
public void testLogBeta() {
    final int ulps = 3;
    for (int i = 0; i < LOG_BETA_REF.length; i++) {
        final double[] ref = LOG_BETA_REF[i];
        final double a = ref[0];
        final double b = ref[1];
        final double expected = ref[2];
        final double actual = Beta.logBeta(a, b);
        final double tol = ulps * FastMath.ulp(expected);
        final StringBuilder builder = new StringBuilder();
        builder.append(a).append(", ").append(b);
        Assert.assertEquals(builder.toString(), expected, actual, tol);
    }
}
 

@Test
public void testLogGammaSum() {
    final int ulps = 2;
    for (int i = 0; i < LOG_GAMMA_SUM_REF.length; i++) {
        final double[] ref = LOG_GAMMA_SUM_REF[i];
        final double a = ref[0];
        final double b = ref[1];
        final double expected = ref[2];
        final double actual = logGammaSum(a, b);
        final double tol = ulps * FastMath.ulp(expected);
        final StringBuilder builder = new StringBuilder();
        builder.append(a).append(", ").append(b);
        Assert.assertEquals(builder.toString(), expected, actual, tol);
    }
}
 

@Test
public void testLogGamma1p() {

    final int ulps = 3;
    for (int i = 0; i < LOG_GAMMA1P_REF.length; i++) {
        final double[] ref = LOG_GAMMA1P_REF[i];
        final double x = ref[0];
        final double expected = ref[1];
        final double actual = Gamma.logGamma1p(x);
        final double tol = ulps * FastMath.ulp(expected);
        Assert.assertEquals(Double.toString(x), expected, actual, tol);
    }
}
 

/** Check the integration span.
 * @param equations set of differential equations
 * @param t target time for the integration
 * @exception NumberIsTooSmallException if integration span is too small
 * @exception DimensionMismatchException if adaptive step size integrators
 * tolerance arrays dimensions are not compatible with equations settings
 */
protected void sanityChecks(final ExpandableStatefulODE equations, final double t)
    throws NumberIsTooSmallException, DimensionMismatchException {

    final double threshold = 1000 * FastMath.ulp(FastMath.max(FastMath.abs(equations.getTime()),
                                                              FastMath.abs(t)));
    final double dt = FastMath.abs(equations.getTime() - t);
    if (dt <= threshold) {
        throw new NumberIsTooSmallException(LocalizedFormats.TOO_SMALL_INTEGRATION_INTERVAL,
                                            dt, threshold, false);
    }

}
 

@Test
public void testLogBeta() {
    final int ulps = 3;
    for (int i = 0; i < LOG_BETA_REF.length; i++) {
        final double[] ref = LOG_BETA_REF[i];
        final double a = ref[0];
        final double b = ref[1];
        final double expected = ref[2];
        final double actual = Beta.logBeta(a, b);
        final double tol = ulps * FastMath.ulp(expected);
        final StringBuilder builder = new StringBuilder();
        builder.append(a).append(", ").append(b);
        Assert.assertEquals(builder.toString(), expected, actual, tol);
    }
}
 

@Test
public void testLogGammaMinusLogGammaSum() {
    final int ulps = 4;
    for (int i = 0; i < LOG_GAMMA_MINUS_LOG_GAMMA_SUM_REF.length; i++) {
        final double[] ref = LOG_GAMMA_MINUS_LOG_GAMMA_SUM_REF[i];
        final double a = ref[0];
        final double b = ref[1];
        final double expected = ref[2];
        final double actual = logGammaMinusLogGammaSum(a, b);
        final double tol = ulps * FastMath.ulp(expected);
        final StringBuilder builder = new StringBuilder();
        builder.append(a).append(", ").append(b);
        Assert.assertEquals(builder.toString(), expected, actual, tol);
    }
}
 

@Test
public void testLogGamma1p() {

    final int ulps = 3;
    for (int i = 0; i < LOG_GAMMA1P_REF.length; i++) {
        final double[] ref = LOG_GAMMA1P_REF[i];
        final double x = ref[0];
        final double expected = ref[1];
        final double actual = Gamma.logGamma1p(x);
        final double tol = ulps * FastMath.ulp(expected);
        Assert.assertEquals(Double.toString(x), expected, actual, tol);
    }
}
 

@Test
public void testLogGammaSum() {
    final int ulps = 2;
    for (int i = 0; i < LOG_GAMMA_SUM_REF.length; i++) {
        final double[] ref = LOG_GAMMA_SUM_REF[i];
        final double a = ref[0];
        final double b = ref[1];
        final double expected = ref[2];
        final double actual = logGammaSum(a, b);
        final double tol = ulps * FastMath.ulp(expected);
        final StringBuilder builder = new StringBuilder();
        builder.append(a).append(", ").append(b);
        Assert.assertEquals(builder.toString(), expected, actual, tol);
    }
}
 

@Test
public void testLogGammaMinusLogGammaSum() {
    final int ulps = 4;
    for (int i = 0; i < LOG_GAMMA_MINUS_LOG_GAMMA_SUM_REF.length; i++) {
        final double[] ref = LOG_GAMMA_MINUS_LOG_GAMMA_SUM_REF[i];
        final double a = ref[0];
        final double b = ref[1];
        final double expected = ref[2];
        final double actual = logGammaMinusLogGammaSum(a, b);
        final double tol = ulps * FastMath.ulp(expected);
        final StringBuilder builder = new StringBuilder();
        builder.append(a).append(", ").append(b);
        Assert.assertEquals(builder.toString(), expected, actual, tol);
    }
}
 

/**
 * Build a differentiator with number of points and step size when independent variable is bounded.
 * <p>
 * When the independent variable is bounded (tLower &lt; t &lt; tUpper), the sampling
 * points used for differentiation will be adapted to ensure the constraint holds
 * even near the boundaries. This means the sample will not be centered anymore in
 * these cases. At an extreme case, computing derivatives exactly at the lower bound
 * will lead the sample to be entirely on the right side of the derivation point.
 * </p>
 * <p>
 * Note that the boundaries are considered to be excluded for function evaluation.
 * </p>
 * <p>
 * Beware that wrong settings for the finite differences differentiator
 * can lead to highly unstable and inaccurate results, especially for
 * high derivation orders. Using very small step sizes is often a
 * <em>bad</em> idea.
 * </p>
 * @param nbPoints number of points to use
 * @param stepSize step size (gap between each point)
 * @param tLower lower bound for independent variable (may be {@code Double.NEGATIVE_INFINITY}
 * if there are no lower bounds)
 * @param tUpper upper bound for independent variable (may be {@code Double.POSITIVE_INFINITY}
 * if there are no upper bounds)
 * @exception NotPositiveException if {@code stepsize <= 0} (note that
 * {@link NotPositiveException} extends {@link NumberIsTooSmallException})
 * @exception NumberIsTooSmallException {@code nbPoint <= 1}
 * @exception NumberIsTooLargeException {@code stepSize * (nbPoints - 1) >= tUpper - tLower}
 */
public FiniteDifferencesDifferentiator(final int nbPoints, final double stepSize,
                                       final double tLower, final double tUpper)
        throws NotPositiveException, NumberIsTooSmallException, NumberIsTooLargeException {

    if (nbPoints <= 1) {
        throw new NumberIsTooSmallException(stepSize, 1, false);
    }
    this.nbPoints = nbPoints;

    if (stepSize <= 0) {
        throw new NotPositiveException(stepSize);
    }
    this.stepSize = stepSize;

    halfSampleSpan = 0.5 * stepSize * (nbPoints - 1);
    if (2 * halfSampleSpan >= tUpper - tLower) {
        throw new NumberIsTooLargeException(2 * halfSampleSpan, tUpper - tLower, false);
    }
    final double safety = FastMath.ulp(halfSampleSpan);
    this.tMin = tLower + halfSampleSpan + safety;
    this.tMax = tUpper - halfSampleSpan - safety;

}
 

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

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

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

/**
 * Build a differentiator with number of points and step size when independent variable is bounded.
 * <p>
 * When the independent variable is bounded (tLower &lt; t &lt; tUpper), the sampling
 * points used for differentiation will be adapted to ensure the constraint holds
 * even near the boundaries. This means the sample will not be centered anymore in
 * these cases. At an extreme case, computing derivatives exactly at the lower bound
 * will lead the sample to be entirely on the right side of the derivation point.
 * </p>
 * <p>
 * Note that the boundaries are considered to be excluded for function evaluation.
 * </p>
 * <p>
 * Beware that wrong settings for the finite differences differentiator
 * can lead to highly unstable and inaccurate results, especially for
 * high derivation orders. Using very small step sizes is often a
 * <em>bad</em> idea.
 * </p>
 * @param nbPoints number of points to use
 * @param stepSize step size (gap between each point)
 * @param tLower lower bound for independent variable (may be {@code Double.NEGATIVE_INFINITY}
 * if there are no lower bounds)
 * @param tUpper upper bound for independent variable (may be {@code Double.POSITIVE_INFINITY}
 * if there are no upper bounds)
 * @exception NotPositiveException if {@code stepsize <= 0} (note that
 * {@link NotPositiveException} extends {@link NumberIsTooSmallException})
 * @exception NumberIsTooSmallException {@code nbPoint <= 1}
 * @exception NumberIsTooLargeException {@code stepSize * (nbPoints - 1) >= tUpper - tLower}
 */
public FiniteDifferencesDifferentiator(final int nbPoints, final double stepSize,
                                       final double tLower, final double tUpper)
        throws NotPositiveException, NumberIsTooSmallException, NumberIsTooLargeException {

    if (nbPoints <= 1) {
        throw new NumberIsTooSmallException(stepSize, 1, false);
    }
    this.nbPoints = nbPoints;

    if (stepSize <= 0) {
        throw new NotPositiveException(stepSize);
    }
    this.stepSize = stepSize;

    halfSampleSpan = 0.5 * stepSize * (nbPoints - 1);
    if (2 * halfSampleSpan >= tUpper - tLower) {
        throw new NumberIsTooLargeException(2 * halfSampleSpan, tUpper - tLower, false);
    }
    final double safety = FastMath.ulp(halfSampleSpan);
    this.tMin = tLower + halfSampleSpan + safety;
    this.tMax = tUpper - halfSampleSpan - safety;

}
 

/**
 * Build a differentiator with number of points and step size when independent variable is bounded.
 * <p>
 * When the independent variable is bounded (tLower &lt; t &lt; tUpper), the sampling
 * points used for differentiation will be adapted to ensure the constraint holds
 * even near the boundaries. This means the sample will not be centered anymore in
 * these cases. At an extreme case, computing derivatives exactly at the lower bound
 * will lead the sample to be entirely on the right side of the derivation point.
 * </p>
 * <p>
 * Note that the boundaries are considered to be excluded for function evaluation.
 * </p>
 * <p>
 * Beware that wrong settings for the finite differences differentiator
 * can lead to highly unstable and inaccurate results, especially for
 * high derivation orders. Using very small step sizes is often a
 * <em>bad</em> idea.
 * </p>
 * @param nbPoints number of points to use
 * @param stepSize step size (gap between each point)
 * @param tLower lower bound for independent variable (may be {@code Double.NEGATIVE_INFINITY}
 * if there are no lower bounds)
 * @param tUpper upper bound for independent variable (may be {@code Double.POSITIVE_INFINITY}
 * if there are no upper bounds)
 * @exception NotPositiveException if {@code stepsize <= 0} (note that
 * {@link NotPositiveException} extends {@link NumberIsTooSmallException})
 * @exception NumberIsTooSmallException {@code nbPoint <= 1}
 * @exception NumberIsTooLargeException {@code stepSize * (nbPoints - 1) >= tUpper - tLower}
 */
public FiniteDifferencesDifferentiator(final int nbPoints, final double stepSize,
                                       final double tLower, final double tUpper)
        throws NotPositiveException, NumberIsTooSmallException, NumberIsTooLargeException {

    if (nbPoints <= 1) {
        throw new NumberIsTooSmallException(stepSize, 1, false);
    }
    this.nbPoints = nbPoints;

    if (stepSize <= 0) {
        throw new NotPositiveException(stepSize);
    }
    this.stepSize = stepSize;

    halfSampleSpan = 0.5 * stepSize * (nbPoints - 1);
    if (2 * halfSampleSpan >= tUpper - tLower) {
        throw new NumberIsTooLargeException(2 * halfSampleSpan, tUpper - tLower, false);
    }
    final double safety = FastMath.ulp(halfSampleSpan);
    this.tMin = tLower + halfSampleSpan + safety;
    this.tMax = tUpper - halfSampleSpan - safety;

}
 

public static SummaryStatistics assessAccuracy(final Method method,
                                               final DataInputStream in,
                                               final DataOutputStream out)
    throws IOException, IllegalAccessException, IllegalArgumentException,
    InvocationTargetException {

    if (method.getReturnType() != Double.TYPE) {
        throw new IllegalArgumentException("method must return a double");
    }

    final Class<?>[] types = method.getParameterTypes();
    for (int i = 0; i < types.length; i++) {
        if (!types[i].isPrimitive()) {
            final StringBuilder builder = new StringBuilder();
            builder.append("argument #").append(i + 1)
                .append(" of method ").append(method.getName())
                .append("must be of primitive of type");
            throw new IllegalArgumentException(builder.toString());
        }
    }

    final SummaryStatistics stat = new SummaryStatistics();
    final Object[] parameters = new Object[types.length];
    while (true) {
        try {
            for (int i = 0; i < parameters.length; i++) {
                parameters[i] = readAndWritePrimitiveValue(in, out,
                                                           types[i]);
            }
            final double expected = in.readDouble();
            if (FastMath.abs(expected) > 1E-16) {
                final Object value = method.invoke(null, parameters);
                final double actual = ((Double) value).doubleValue();
                final double err = FastMath.abs(actual - expected);
                final double ulps = err / FastMath.ulp(expected);
                out.writeDouble(expected);
                out.writeDouble(actual);
                out.writeDouble(ulps);
                stat.addValue(ulps);
            }
        } catch (EOFException e) {
            break;
        }
    }
    return stat;
}