Java Code Examples for org.apache.commons.math3.analysis.UnivariateFunction#value()

private void checkOrthogonality(final PolynomialFunction p1,
                                final PolynomialFunction p2,
                                final UnivariateFunction weight,
                                final double a, final double b,
                                final double nonZeroThreshold,
                                final double zeroThreshold) {
    UnivariateFunction f = new UnivariateFunction() {
        public double value(double x) {
            return weight.value(x) * p1.value(x) * p2.value(x);
        }
    };
    double dotProduct =
            new IterativeLegendreGaussIntegrator(5, 1.0e-9, 1.0e-8, 2, 15).integrate(1000000, f, a, b);
    if (p1.degree() == p2.degree()) {
        // integral should be non-zero
        Assert.assertTrue("I(" + p1.degree() + ", " + p2.degree() + ") = "+ dotProduct,
                          FastMath.abs(dotProduct) > nonZeroThreshold);
    } else {
        // integral should be zero
        Assert.assertEquals("I(" + p1.degree() + ", " + p2.degree() + ") = "+ dotProduct,
                            0.0, FastMath.abs(dotProduct), zeroThreshold);
    }
}
 

@Test
public void testShortcut() {
    final Sinc s = new Sinc();
    final UnivariateFunction f = new UnivariateFunction() {
        public double value(double x) {
            Dfp dfpX = new DfpField(25).newDfp(x);
            return DfpMath.sin(dfpX).divide(dfpX).toDouble();
        }
    };

    for (double x = 1e-30; x < 1e10; x *= 2) {
        final double fX = f.value(x);
        final double sX = s.value(x);
        Assert.assertEquals("x=" + x, fX, sX, 2.0e-16);
    }
}
 

/**
 * Test of parameters for the interpolator.
 */
@Test
public void testParameters() {
    UnivariateInterpolator interpolator = new NevilleInterpolator();

    try {
        // bad abscissas array
        double x[] = { 1.0, 2.0, 2.0, 4.0 };
        double y[] = { 0.0, 4.0, 4.0, 2.5 };
        UnivariateFunction p = interpolator.interpolate(x, y);
        p.value(0.0);
        Assert.fail("Expecting NonMonotonicSequenceException - bad abscissas array");
    } catch (NonMonotonicSequenceException ex) {
        // expected
    }
}
 

/**
 * Check that the endpoints specify an interval and the end points
 * bracket a root.
 *
 * @param function Function.
 * @param lower Lower endpoint.
 * @param upper Upper endpoint.
 * @throws NoBracketingException if the function has the same sign at the
 * endpoints.
 * @throws NullArgumentException if {@code function} is {@code null}.
 */
public static void verifyBracketing(UnivariateFunction function,
                                    final double lower,
                                    final double upper)
    throws NullArgumentException,
           NoBracketingException {
    if (function == null) {
        throw new NullArgumentException(LocalizedFormats.FUNCTION);
    }
    verifyInterval(lower, upper);
    if (!isBracketing(function, lower, upper)) {
        throw new NoBracketingException(lower, upper,
                                        function.value(lower),
                                        function.value(upper));
    }
}
 

/**
 * @param f Function.
 * @param x Argument.
 * @return {@code f(x)}
 * @throws TooManyEvaluationsException if the maximal number of evaluations is
 * exceeded.
 */
private double eval(UnivariateFunction f, double x) {
    try {
        evaluations.incrementCount();
    } catch (MaxCountExceededException e) {
        throw new TooManyEvaluationsException(e.getMax());
    }
    return f.value(x);
}
 

@Test
public void testInverse() {
    final UnivariateFunction inv = new Inverse();
    final UnivariateFunction log = new Log();

    final double lo = 12.34;
    final double hi = 456.78;

    final GaussIntegrator integrator = factory.legendre(60, lo, hi);
    final double s = integrator.integrate(inv);
    final double expected = log.value(hi) - log.value(lo);
    // System.out.println("s=" + s + " e=" + expected);
    Assert.assertEquals(expected, s, 1e-14);
}
 

@Test
public void testMoreThanOnePeriodCoverage() {
    final int n = 30;
    final double[] xval = new double[n];
    final double[] yval = new double[n];
    final double period = 12.3;
    final double offset = 45.67;

    double delta = period / 2;
    for (int i = 0; i < n; i++) {
        delta += 10 * period / n * rng.nextDouble();
        xval[i] = offset + delta;
        yval[i] = FastMath.sin(xval[i]);
    }

    final UnivariateInterpolator interP
        = new UnivariatePeriodicInterpolator(new LinearInterpolator(),
                                                 period, 1);
    final UnivariateFunction fP = interP.interpolate(xval, yval);

    // Test interpolation outside the sample data interval.
    for (int i = 0; i < n; i++) {
        final double xIn = offset + rng.nextDouble() * period;
        final double xOut = xIn + rng.nextInt(123456789) * period;
        final double yIn = fP.value(xIn);
        final double yOut = fP.value(xOut);

        Assert.assertEquals(yIn, yOut, 1e-6);
    }
}
 

/**
 * @param f Function.
 * @param x Argument.
 * @return {@code f(x)}
 * @throws TooManyEvaluationsException if the maximal number of evaluations is
 * exceeded.
 */
private double eval(UnivariateFunction f, double x) {
    try {
        evaluations.incrementCount();
    } catch (MaxCountExceededException e) {
        throw new TooManyEvaluationsException(e.getMax());
    }
    return f.value(x);
}
 

/**
 * @param f Function.
 * @param x Argument.
 * @return {@code f(x)}
 * @throws TooManyEvaluationsException if the maximal number of evaluations is
 * exceeded.
 */
private double eval(UnivariateFunction f, double x) {
    try {
        evaluations.incrementCount();
    } catch (MaxCountExceededException e) {
        throw new TooManyEvaluationsException(e.getMax());
    }
    return f.value(x);
}
 

@Test
public void testComparison() {
    final Sqrt s = new Sqrt();
    final UnivariateFunction f = new UnivariateFunction() {
            public double value(double x) {
                return Math.sqrt(x);
            }
        };

    for (double x = 1e-30; x < 1e10; x *= 2) {
        final double fX = f.value(x);
        final double sX = s.value(x);
        Assert.assertEquals("x=" + x, fX, sX, 0);
    }
}
 

/**
 * Returns an estimate of the integral of {@code f(x) * w(x)},
 * where {@code w} is a weight function that depends on the actual
 * flavor of the Gauss integration scheme.
 * The algorithm uses the points and associated weights, as passed
 * to the {@link #GaussIntegrator(double[],double[]) constructor}.
 *
 * @param f Function to integrate.
 * @return the integral of the weighted function.
 */
public double integrate(UnivariateFunction f) {
    double s = 0;
    double c = 0;
    for (int i = 0; i < points.length; i++) {
        final double x = points[i];
        final double w = weights[i];
        final double y = w * f.value(x) - c;
        final double t = s + y;
        c = (t - s) - y;
        s = t;
    }
    return s;
}
 

@Test
public void testShortcut() {
    final Sinc s = new Sinc();
    final UnivariateFunction f = new UnivariateFunction() {
            public double value(double x) {
                return FastMath.sin(x) / x;
            }
        };

    for (double x = 1e-30; x < 1e10; x *= 2) {
        final double fX = f.value(x);
        final double sX = s.value(x);
        Assert.assertEquals("x=" + x, fX, sX, 0);
    }
}
 

@Test
public void testMoreThanOnePeriodCoverage() {
    final int n = 30;
    final double[] xval = new double[n];
    final double[] yval = new double[n];
    final double period = 12.3;
    final double offset = 45.67;

    double delta = period / 2;
    for (int i = 0; i < n; i++) {
        delta += 10 * period / n * rng.nextDouble();
        xval[i] = offset + delta;
        yval[i] = FastMath.sin(xval[i]);
    }

    final UnivariateInterpolator interP
        = new UnivariatePeriodicInterpolator(new LinearInterpolator(),
                                                 period, 1);
    final UnivariateFunction fP = interP.interpolate(xval, yval);

    // Test interpolation outside the sample data interval.
    for (int i = 0; i < n; i++) {
        final double xIn = offset + rng.nextDouble() * period;
        final double xOut = xIn + rng.nextInt(123456789) * period;
        final double yIn = fP.value(xIn);
        final double yOut = fP.value(xOut);

        Assert.assertEquals(yIn, yOut, 1e-6);
    }
}
 

@Test
public void testMoreThanOnePeriodCoverage() {
    final int n = 30;
    final double[] xval = new double[n];
    final double[] yval = new double[n];
    final double period = 12.3;
    final double offset = 45.67;

    double delta = period / 2;
    for (int i = 0; i < n; i++) {
        delta += 10 * period / n * rng.nextDouble();
        xval[i] = offset + delta;
        yval[i] = FastMath.sin(xval[i]);
    }

    final UnivariateInterpolator interP
        = new UnivariatePeriodicInterpolator(new LinearInterpolator(),
                                                 period, 1);
    final UnivariateFunction fP = interP.interpolate(xval, yval);

    // Test interpolation outside the sample data interval.
    for (int i = 0; i < n; i++) {
        final double xIn = offset + rng.nextDouble() * period;
        final double xOut = xIn + rng.nextInt(123456789) * period;
        final double yIn = fP.value(xIn);
        final double yOut = fP.value(xOut);

        Assert.assertEquals(yIn, yOut, 1e-6);
    }
}
 

/**
 * @param f Function.
 * @param x Argument.
 * @return {@code f(x)}
 * @throws TooManyEvaluationsException if the maximal number of evaluations is
 * exceeded.
 */
private double eval(UnivariateFunction f, double x) {
    try {
        evaluations.incrementCount();
    } catch (MaxCountExceededException e) {
        throw new TooManyEvaluationsException(e.getMax());
    }
    return f.value(x);
}
 

@Test
public void testShortcut() {
    final Sinc s = new Sinc();
    final UnivariateFunction f = new UnivariateFunction() {
            public double value(double x) {
                return FastMath.sin(x) / x;
            }
        };

    for (double x = 1e-30; x < 1e10; x *= 2) {
        final double fX = f.value(x);
        final double sX = s.value(x);
        Assert.assertEquals("x=" + x, fX, sX, 0);
    }
}
 

@Test
public void testSine() {
    final int n = 30;
    final double[] xval = new double[n];
    final double[] yval = new double[n];
    final double period = 12.3;
    final double offset = 45.67;

    double delta = 0;
    for (int i = 0; i < n; i++) {
        delta += rng.nextDouble() * period / n;
        xval[i] = offset + delta;
        yval[i] = FastMath.sin(xval[i]);
    }

    final UnivariateInterpolator inter = new LinearInterpolator();
    final UnivariateFunction f = inter.interpolate(xval, yval);

    final UnivariateInterpolator interP
        = new UnivariatePeriodicInterpolator(new LinearInterpolator(),
                                                 period, 1);
    final UnivariateFunction fP = interP.interpolate(xval, yval);

    // Comparing with original interpolation algorithm.
    final double xMin = xval[0];
    final double xMax = xval[n - 1];
    for (int i = 0; i < n; i++) {
        final double x = xMin + (xMax - xMin) * rng.nextDouble();
        final double y = f.value(x);
        final double yP = fP.value(x);

        Assert.assertEquals("x=" + x, y, yP, Math.ulp(1d));
    }

    // Test interpolation outside the primary interval.
    for (int i = 0; i < n; i++) {
        final double xIn = offset + rng.nextDouble() * period;
        final double xOut = xIn + rng.nextInt(123456789) * period;
        final double yIn = fP.value(xIn);
        final double yOut = fP.value(xOut);

        Assert.assertEquals(yIn, yOut, 1e-7);
    }
}
 

/** Force a root found by a non-bracketing solver to lie on a specified side,
 * as if the solver was a bracketing one.
 * @param maxEval maximal number of new evaluations of the function
 * (evaluations already done for finding the root should have already been subtracted
 * from this number)
 * @param f function to solve
 * @param bracketing bracketing solver to use for shifting the root
 * @param baseRoot original root found by a previous non-bracketing solver
 * @param min minimal bound of the search interval
 * @param max maximal bound of the search interval
 * @param allowedSolution the kind of solutions that the root-finding algorithm may
 * accept as solutions.
 * @return a root approximation, on the specified side of the exact root
 * @throws NoBracketingException if the function has the same sign at the
 * endpoints.
 */
public static double forceSide(final int maxEval, final UnivariateFunction f,
                               final BracketedUnivariateSolver<UnivariateFunction> bracketing,
                               final double baseRoot, final double min, final double max,
                               final AllowedSolution allowedSolution)
    throws NoBracketingException {

    if (allowedSolution == AllowedSolution.ANY_SIDE) {
        // no further bracketing required
        return baseRoot;
    }

    // find a very small interval bracketing the root
    final double step = FastMath.max(bracketing.getAbsoluteAccuracy(),
                                     FastMath.abs(baseRoot * bracketing.getRelativeAccuracy()));
    double xLo        = FastMath.max(min, baseRoot - step);
    double fLo        = f.value(xLo);
    double xHi        = FastMath.min(max, baseRoot + step);
    double fHi        = f.value(xHi);
    int remainingEval = maxEval - 2;
    while (remainingEval > 0) {

        if ((fLo >= 0 && fHi <= 0) || (fLo <= 0 && fHi >= 0)) {
            // compute the root on the selected side
            return bracketing.solve(remainingEval, f, xLo, xHi, baseRoot, allowedSolution);
        }

        // try increasing the interval
        boolean changeLo = false;
        boolean changeHi = false;
        if (fLo < fHi) {
            // increasing function
            if (fLo >= 0) {
                changeLo = true;
            } else {
                changeHi = true;
            }
        } else if (fLo > fHi) {
            // decreasing function
            if (fLo <= 0) {
                changeLo = true;
            } else {
                changeHi = true;
            }
        } else {
            // unknown variation
            changeLo = true;
            changeHi = true;
        }

        // update the lower bound
        if (changeLo) {
            xLo = FastMath.max(min, xLo - step);
            fLo  = f.value(xLo);
            remainingEval--;
        }

        // update the higher bound
        if (changeHi) {
            xHi = FastMath.min(max, xHi + step);
            fHi  = f.value(xHi);
            remainingEval--;
        }

    }

    throw new NoBracketingException(LocalizedFormats.FAILED_BRACKETING,
                                    xLo, xHi, fLo, fHi,
                                    maxEval - remainingEval, maxEval, baseRoot,
                                    min, max);

}
 

/**
 * Check whether the interval bounds bracket a root. That is, if the
 * values at the endpoints are not equal to zero, then the function takes
 * opposite signs at the endpoints.
 *
 * @param function Function.
 * @param lower Lower endpoint.
 * @param upper Upper endpoint.
 * @return {@code true} if the function values have opposite signs at the
 * given points.
 * @throws NullArgumentException if {@code function} is {@code null}.
 */
public static boolean isBracketing(UnivariateFunction function,
                                   final double lower,
                                   final double upper)
    throws NullArgumentException {
    if (function == null) {
        throw new NullArgumentException(LocalizedFormats.FUNCTION);
    }
    final double fLo = function.value(lower);
    final double fHi = function.value(upper);
    return (fLo >= 0 && fHi <= 0) || (fLo <= 0 && fHi >= 0);
}
 

/**
 * Check whether the interval bounds bracket a root. That is, if the
 * values at the endpoints are not equal to zero, then the function takes
 * opposite signs at the endpoints.
 *
 * @param function Function.
 * @param lower Lower endpoint.
 * @param upper Upper endpoint.
 * @return {@code true} if the function values have opposite signs at the
 * given points.
 * @throws NullArgumentException if {@code function} is {@code null}.
 */
public static boolean isBracketing(UnivariateFunction function,
                                   final double lower,
                                   final double upper)
    throws NullArgumentException {
    if (function == null) {
        throw new NullArgumentException(LocalizedFormats.FUNCTION);
    }
    final double fLo = function.value(lower);
    final double fHi = function.value(upper);
    return (fLo >= 0 && fHi <= 0) || (fLo <= 0 && fHi >= 0);
}