org.apache.commons.math.exception.NoBracketingException Java Examples

/**
 * This method attempts to find two values a and b satisfying <ul>
 * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
 * <li> <code> f(a) * f(b) <= 0 </code> </li>
 * </ul>
 * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
 * and <code>b</code> bracket a root of f.
 * <p>
 * The algorithm starts by setting
 * <code>a := initial -1; b := initial +1,</code> examines the value of the
 * function at <code>a</code> and <code>b</code> and keeps moving
 * the endpoints out by one unit each time through a loop that terminates
 * when one of the following happens: <ul>
 * <li> <code> f(a) * f(b) <= 0 </code> --  success!</li>
 * <li> <code> a = lower </code> and <code> b = upper</code>
 * -- ConvergenceException </li>
 * <li> <code> maximumIterations</code> iterations elapse
 * -- ConvergenceException </li></ul></p>
 *
 * @param function Function.
 * @param initial Initial midpoint of interval being expanded to
 * bracket a root.
 * @param lowerBound Lower bound (a is never lower than this value).
 * @param upperBound Upper bound (b never is greater than this
 * value).
 * @param maximumIterations Maximum number of iterations to perform
 * @return a two element array holding a and b.
 * @throws NoBracketingException if the algorithm fails to find a and b
 * satisfying the desired conditions.
 * @throws IllegalArgumentException if function is null, maximumIterations
 * is not positive, or initial is not between lowerBound and upperBound.
 */
public static double[] bracket(UnivariateRealFunction function,
                               double initial,
                               double lowerBound, double upperBound,
                               int maximumIterations)  {
    if (function == null) {
        throw new NullArgumentException(LocalizedFormats.FUNCTION);
    }
    if (maximumIterations <= 0)  {
        throw new NotStrictlyPositiveException(LocalizedFormats.INVALID_MAX_ITERATIONS, maximumIterations);
    }
    verifySequence(lowerBound, initial, upperBound);

    double a = initial;
    double b = initial;
    double fa;
    double fb;
    int numIterations = 0;

    do {
        a = FastMath.max(a - 1.0, lowerBound);
        b = FastMath.min(b + 1.0, upperBound);
        fa = function.value(a);

        fb = function.value(b);
        ++numIterations;
    } while ((fa * fb > 0.0) && (numIterations < maximumIterations) &&
            ((a > lowerBound) || (b < upperBound)));

    if (fa * fb > 0.0) {
        throw new NoBracketingException(LocalizedFormats.FAILED_BRACKETING,
                                        a, b, fa, fb,
                                        numIterations, maximumIterations, initial,
                                        lowerBound, upperBound);
    }

    return new double[] {a, b};
}
 

/**
 * {@inheritDoc}
 */
@Override
public double doSolve() {
    double min = getMin();
    double max = getMax();
    double initial = getStartValue();
    final double functionValueAccuracy = getFunctionValueAccuracy();

    verifySequence(min, initial, max);

    // Return the initial guess if it is good enough.
    double yInitial = computeObjectiveValue(initial);
    if (FastMath.abs(yInitial) <= functionValueAccuracy) {
        return initial;
    }

    // Return the first endpoint if it is good enough.
    double yMin = computeObjectiveValue(min);
    if (FastMath.abs(yMin) <= functionValueAccuracy) {
        return min;
    }

    // Reduce interval if min and initial bracket the root.
    if (yInitial * yMin < 0) {
        return laguerre(min, initial, yMin, yInitial);
    }

    // Return the second endpoint if it is good enough.
    double yMax = computeObjectiveValue(max);
    if (FastMath.abs(yMax) <= functionValueAccuracy) {
        return max;
    }

    // Reduce interval if initial and max bracket the root.
    if (yInitial * yMax < 0) {
        return laguerre(initial, max, yInitial, yMax);
    }

    throw new NoBracketingException(min, max, yMin, yMax);
}
 

/**
 * 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 function has the same sign at the
 * endpoints.
 */
public static void verifyBracketing(UnivariateRealFunction function,
                                    final double lower,
                                    final double upper) {
    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));
    }
}
 

/**
 * This method attempts to find two values a and b satisfying <ul>
 * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
 * <li> <code> f(a) * f(b) <= 0 </code> </li>
 * </ul>
 * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
 * and <code>b</code> bracket a root of f.
 * <p>
 * The algorithm starts by setting
 * <code>a := initial -1; b := initial +1,</code> examines the value of the
 * function at <code>a</code> and <code>b</code> and keeps moving
 * the endpoints out by one unit each time through a loop that terminates
 * when one of the following happens: <ul>
 * <li> <code> f(a) * f(b) <= 0 </code> --  success!</li>
 * <li> <code> a = lower </code> and <code> b = upper</code>
 * -- ConvergenceException </li>
 * <li> <code> maximumIterations</code> iterations elapse
 * -- ConvergenceException </li></ul></p>
 *
 * @param function Function.
 * @param initial Initial midpoint of interval being expanded to
 * bracket a root.
 * @param lowerBound Lower bound (a is never lower than this value).
 * @param upperBound Upper bound (b never is greater than this
 * value).
 * @param maximumIterations Maximum number of iterations to perform
 * @return a two element array holding a and b.
 * @throws NoBracketingException if the algorithm fails to find a and b
 * satisfying the desired conditions.
 * @throws IllegalArgumentException if function is null, maximumIterations
 * is not positive, or initial is not between lowerBound and upperBound.
 */
public static double[] bracket(UnivariateRealFunction function,
                               double initial,
                               double lowerBound, double upperBound,
                               int maximumIterations)  {
    if (function == null) {
        throw new NullArgumentException(LocalizedFormats.FUNCTION);
    }
    if (maximumIterations <= 0)  {
        throw new NotStrictlyPositiveException(LocalizedFormats.INVALID_MAX_ITERATIONS, maximumIterations);
    }
    verifySequence(lowerBound, initial, upperBound);

    double a = initial;
    double b = initial;
    double fa;
    double fb;
    int numIterations = 0;

    do {
        a = FastMath.max(a - 1.0, lowerBound);
        b = FastMath.min(b + 1.0, upperBound);
        fa = function.value(a);

        fb = function.value(b);
        ++numIterations;
    } while ((fa * fb > 0.0) && (numIterations < maximumIterations) &&
            ((a > lowerBound) || (b < upperBound)));

    if (fa * fb > 0.0) {
        throw new NoBracketingException(LocalizedFormats.FAILED_BRACKETING,
                                        a, b, fa, fb,
                                        numIterations, maximumIterations, initial,
                                        lowerBound, upperBound);
    }

    return new double[] {a, b};
}
 

/**
 * {@inheritDoc}
 */
@Override
protected double doSolve() {
    double min = getMin();
    double max = getMax();
    final double initial = getStartValue();
    final double functionValueAccuracy = getFunctionValueAccuracy();

    verifySequence(min, initial, max);

    // Return the initial guess if it is good enough.
    double yInitial = computeObjectiveValue(initial);
    if (FastMath.abs(yInitial) <= functionValueAccuracy) {
        return initial;
    }

    // Return the first endpoint if it is good enough.
    double yMin = computeObjectiveValue(min);
    if (FastMath.abs(yMin) <= functionValueAccuracy) {
        return min;
    }

    // Reduce interval if min and initial bracket the root.
    if (yInitial * yMin < 0) {
        return brent(min, initial, yMin, yInitial);
    }

    // Return the second endpoint if it is good enough.
    double yMax = computeObjectiveValue(max);
    if (FastMath.abs(yMax) <= functionValueAccuracy) {
        return max;
    }

    // Reduce interval if initial and max bracket the root.
    if (yInitial * yMax < 0) {
        return brent(initial, max, yInitial, yMax);
    }

    throw new NoBracketingException(min, max, yMin, yMax);
}
 

/**
 * {@inheritDoc}
 */
@Override
public double doSolve() {
    double min = getMin();
    double max = getMax();
    double initial = getStartValue();
    final double functionValueAccuracy = getFunctionValueAccuracy();

    verifySequence(min, initial, max);

    // Return the initial guess if it is good enough.
    double yInitial = computeObjectiveValue(initial);
    if (FastMath.abs(yInitial) <= functionValueAccuracy) {
        return initial;
    }

    // Return the first endpoint if it is good enough.
    double yMin = computeObjectiveValue(min);
    if (FastMath.abs(yMin) <= functionValueAccuracy) {
        return min;
    }

    // Reduce interval if min and initial bracket the root.
    if (yInitial * yMin < 0) {
        return laguerre(min, initial, yMin, yInitial);
    }

    // Return the second endpoint if it is good enough.
    double yMax = computeObjectiveValue(max);
    if (FastMath.abs(yMax) <= functionValueAccuracy) {
        return max;
    }

    // Reduce interval if initial and max bracket the root.
    if (yInitial * yMax < 0) {
        return laguerre(initial, max, yInitial, yMax);
    }

    throw new NoBracketingException(min, max, yMin, yMax);
}
 

/**
 * 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 function has the same sign at the
 * endpoints.
 */
public static void verifyBracketing(UnivariateRealFunction function,
                                    final double lower,
                                    final double upper) {
    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));
    }
}
 

/**
 * This method attempts to find two values a and b satisfying <ul>
 * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
 * <li> <code> f(a) * f(b) <= 0 </code> </li>
 * </ul>
 * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
 * and <code>b</code> bracket a root of f.
 * <p>
 * The algorithm starts by setting
 * <code>a := initial -1; b := initial +1,</code> examines the value of the
 * function at <code>a</code> and <code>b</code> and keeps moving
 * the endpoints out by one unit each time through a loop that terminates
 * when one of the following happens: <ul>
 * <li> <code> f(a) * f(b) <= 0 </code> --  success!</li>
 * <li> <code> a = lower </code> and <code> b = upper</code>
 * -- ConvergenceException </li>
 * <li> <code> maximumIterations</code> iterations elapse
 * -- ConvergenceException </li></ul></p>
 *
 * @param function Function.
 * @param initial Initial midpoint of interval being expanded to
 * bracket a root.
 * @param lowerBound Lower bound (a is never lower than this value).
 * @param upperBound Upper bound (b never is greater than this
 * value).
 * @param maximumIterations Maximum number of iterations to perform
 * @return a two element array holding a and b.
 * @throws NoBracketingException if the algorithm fails to find a and b
 * satisfying the desired conditions.
 * @throws IllegalArgumentException if function is null, maximumIterations
 * is not positive, or initial is not between lowerBound and upperBound.
 */
public static double[] bracket(UnivariateRealFunction function,
                               double initial,
                               double lowerBound, double upperBound,
                               int maximumIterations)  {
    if (function == null) {
        throw new NullArgumentException(LocalizedFormats.FUNCTION);
    }
    if (maximumIterations <= 0)  {
        throw new NotStrictlyPositiveException(LocalizedFormats.INVALID_MAX_ITERATIONS, maximumIterations);
    }
    verifySequence(lowerBound, initial, upperBound);

    double a = initial;
    double b = initial;
    double fa;
    double fb;
    int numIterations = 0;

    do {
        a = FastMath.max(a - 1.0, lowerBound);
        b = FastMath.min(b + 1.0, upperBound);
        fa = function.value(a);

        fb = function.value(b);
        ++numIterations;
    } while ((fa * fb > 0.0) && (numIterations < maximumIterations) &&
            ((a > lowerBound) || (b < upperBound)));

    if (fa * fb > 0.0) {
        throw new NoBracketingException(LocalizedFormats.FAILED_BRACKETING,
                                        a, b, fa, fb,
                                        numIterations, maximumIterations, initial,
                                        lowerBound, upperBound);
    }

    return new double[] {a, b};
}
 

/**
 * {@inheritDoc}
 */
@Override
protected double doSolve() {
    double min = getMin();
    double max = getMax();
    final double initial = getStartValue();
    final double functionValueAccuracy = getFunctionValueAccuracy();

    verifySequence(min, initial, max);

    // Return the initial guess if it is good enough.
    double yInitial = computeObjectiveValue(initial);
    if (FastMath.abs(yInitial) <= functionValueAccuracy) {
        return initial;
    }

    // Return the first endpoint if it is good enough.
    double yMin = computeObjectiveValue(min);
    if (FastMath.abs(yMin) <= functionValueAccuracy) {
        return min;
    }

    // Reduce interval if min and initial bracket the root.
    if (yInitial * yMin < 0) {
        return brent(min, initial, yMin, yInitial);
    }

    // Return the second endpoint if it is good enough.
    double yMax = computeObjectiveValue(max);
    if (FastMath.abs(yMax) <= functionValueAccuracy) {
        return max;
    }

    // Reduce interval if initial and max bracket the root.
    if (yInitial * yMax < 0) {
        return brent(initial, max, yInitial, yMax);
    }

    throw new NoBracketingException(min, max, yMin, yMax);
}
 

/**
 * {@inheritDoc}
 */
@Override
public double doSolve() {
    double min = getMin();
    double max = getMax();
    double initial = getStartValue();
    final double functionValueAccuracy = getFunctionValueAccuracy();

    verifySequence(min, initial, max);

    // Return the initial guess if it is good enough.
    double yInitial = computeObjectiveValue(initial);
    if (FastMath.abs(yInitial) <= functionValueAccuracy) {
        return initial;
    }

    // Return the first endpoint if it is good enough.
    double yMin = computeObjectiveValue(min);
    if (FastMath.abs(yMin) <= functionValueAccuracy) {
        return min;
    }

    // Reduce interval if min and initial bracket the root.
    if (yInitial * yMin < 0) {
        return laguerre(min, initial, yMin, yInitial);
    }

    // Return the second endpoint if it is good enough.
    double yMax = computeObjectiveValue(max);
    if (FastMath.abs(yMax) <= functionValueAccuracy) {
        return max;
    }

    // Reduce interval if initial and max bracket the root.
    if (yInitial * yMax < 0) {
        return laguerre(initial, max, yInitial, yMax);
    }

    throw new NoBracketingException(min, max, yMin, yMax);
}
 

/**
 * Check that the endpoints specify an interval and the function takes
 * opposite signs at the endpoints.
 *
 * @param function Function.
 * @param lower Lower endpoint.
 * @param upper Upper endpoint.
 * @throws NoBracketingException if function has the same sign at the
 * endpoints.
 */
public static void verifyBracketing(UnivariateRealFunction function,
                                    final double lower,
                                    final double upper) {
    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));
    }
}
 

/**
 * {@inheritDoc}
 */
@Override
protected double doSolve() {
    double min = getMin();
    double max = getMax();
    final double initial = getStartValue();
    final double functionValueAccuracy = getFunctionValueAccuracy();

    verifySequence(min, initial, max);

    // Return the initial guess if it is good enough.
    double yInitial = computeObjectiveValue(initial);
    if (FastMath.abs(yInitial) <= functionValueAccuracy) {
        return initial;
    }

    // Return the first endpoint if it is good enough.
    double yMin = computeObjectiveValue(min);
    if (FastMath.abs(yMin) <= functionValueAccuracy) {
        return min;
    }

    // Reduce interval if min and initial bracket the root.
    if (yInitial * yMin < 0) {
        return brent(min, initial, yMin, yInitial);
    }

    // Return the second endpoint if it is good enough.
    double yMax = computeObjectiveValue(max);
    if (FastMath.abs(yMax) <= functionValueAccuracy) {
        return max;
    }

    // Reduce interval if initial and max bracket the root.
    if (yInitial * yMax < 0) {
        return brent(initial, max, yInitial, yMax);
    }

    throw new NoBracketingException(min, max, yMin, yMax);
}
 

/** 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
 */
public static double forceSide(final int maxEval, final UnivariateRealFunction f,
                               final BracketedUnivariateRealSolver<UnivariateRealFunction> bracketing,
                               final double baseRoot, final double min, final double max,
                               final AllowedSolution allowedSolution) {

    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);

}
 

/**
 * {@inheritDoc}
 */
@Override
protected double doSolve() {
    final double min = getMin();
    final double max = getMax();

    verifyInterval(min, max);

    final double relativeAccuracy = getRelativeAccuracy();
    final double absoluteAccuracy = getAbsoluteAccuracy();
    final double functionValueAccuracy = getFunctionValueAccuracy();

    // x2 is the last root approximation
    // x is the new approximation and new x2 for next round
    // x0 < x1 < x2 does not hold here

    double x0 = min;
    double y0 = computeObjectiveValue(x0);
    if (FastMath.abs(y0) < functionValueAccuracy) {
        return x0;
    }
    double x1 = max;
    double y1 = computeObjectiveValue(x1);
    if (FastMath.abs(y1) < functionValueAccuracy) {
        return x1;
    }

    if(y0 * y1 > 0) {
        throw new NoBracketingException(x0, x1, y0, y1);
    }

    double x2 = 0.5 * (x0 + x1);
    double y2 = computeObjectiveValue(x2);

    double oldx = Double.POSITIVE_INFINITY;
    while (true) {
        // quadratic interpolation through x0, x1, x2
        final double q = (x2 - x1) / (x1 - x0);
        final double a = q * (y2 - (1 + q) * y1 + q * y0);
        final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
        final double c = (1 + q) * y2;
        final double delta = b * b - 4 * a * c;
        double x;
        final double denominator;
        if (delta >= 0.0) {
            // choose a denominator larger in magnitude
            double dplus = b + FastMath.sqrt(delta);
            double dminus = b - FastMath.sqrt(delta);
            denominator = FastMath.abs(dplus) > FastMath.abs(dminus) ? dplus : dminus;
        } else {
            // take the modulus of (B +/- FastMath.sqrt(delta))
            denominator = FastMath.sqrt(b * b - delta);
        }
        if (denominator != 0) {
            x = x2 - 2.0 * c * (x2 - x1) / denominator;
            // perturb x if it exactly coincides with x1 or x2
            // the equality tests here are intentional
            while (x == x1 || x == x2) {
                x += absoluteAccuracy;
            }
        } else {
            // extremely rare case, get a random number to skip it
            x = min + FastMath.random() * (max - min);
            oldx = Double.POSITIVE_INFINITY;
        }
        final double y = computeObjectiveValue(x);

        // check for convergence
        final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
        if (FastMath.abs(x - oldx) <= tolerance ||
            FastMath.abs(y) <= functionValueAccuracy) {
            return x;
        }

        // prepare the next iteration
        x0 = x1;
        y0 = y1;
        x1 = x2;
        y1 = y2;
        x2 = x;
        y2 = y;
        oldx = x;
    }
}
 

/**
 * {@inheritDoc}
 */
@Override
protected double doSolve() {
    final double min = getMin();
    final double max = getMax();

    verifyInterval(min, max);

    final double relativeAccuracy = getRelativeAccuracy();
    final double absoluteAccuracy = getAbsoluteAccuracy();
    final double functionValueAccuracy = getFunctionValueAccuracy();

    // x2 is the last root approximation
    // x is the new approximation and new x2 for next round
    // x0 < x1 < x2 does not hold here

    double x0 = min;
    double y0 = computeObjectiveValue(x0);
    if (FastMath.abs(y0) < functionValueAccuracy) {
        return x0;
    }
    double x1 = max;
    double y1 = computeObjectiveValue(x1);
    if (FastMath.abs(y1) < functionValueAccuracy) {
        return x1;
    }

    if(y0 * y1 > 0) {
        throw new NoBracketingException(x0, x1, y0, y1);
    }

    double x2 = 0.5 * (x0 + x1);
    double y2 = computeObjectiveValue(x2);

    double oldx = Double.POSITIVE_INFINITY;
    while (true) {
        // quadratic interpolation through x0, x1, x2
        final double q = (x2 - x1) / (x1 - x0);
        final double a = q * (y2 - (1 + q) * y1 + q * y0);
        final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
        final double c = (1 + q) * y2;
        final double delta = b * b - 4 * a * c;
        double x;
        final double denominator;
        if (delta >= 0.0) {
            // choose a denominator larger in magnitude
            double dplus = b + FastMath.sqrt(delta);
            double dminus = b - FastMath.sqrt(delta);
            denominator = FastMath.abs(dplus) > FastMath.abs(dminus) ? dplus : dminus;
        } else {
            // take the modulus of (B +/- FastMath.sqrt(delta))
            denominator = FastMath.sqrt(b * b - delta);
        }
        if (denominator != 0) {
            x = x2 - 2.0 * c * (x2 - x1) / denominator;
            // perturb x if it exactly coincides with x1 or x2
            // the equality tests here are intentional
            while (x == x1 || x == x2) {
                x += absoluteAccuracy;
            }
        } else {
            // extremely rare case, get a random number to skip it
            x = min + FastMath.random() * (max - min);
            oldx = Double.POSITIVE_INFINITY;
        }
        final double y = computeObjectiveValue(x);

        // check for convergence
        final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
        if (FastMath.abs(x - oldx) <= tolerance ||
            FastMath.abs(y) <= functionValueAccuracy) {
            return x;
        }

        // prepare the next iteration
        x0 = x1;
        y0 = y1;
        x1 = x2;
        y1 = y2;
        x2 = x;
        y2 = y;
        oldx = x;
    }
}
 

/**
 * {@inheritDoc}
 */
@Override
protected double doSolve() {
    double min = getMin();
    double max = getMax();
    verifyInterval(min, max);

    final double functionValueAccuracy = getFunctionValueAccuracy();

    // Index 0 is the old approximation for the root.
    // Index 1 is the last calculated approximation  for the root.
    // Index 2 is a bracket for the root with respect to x0.
    // OldDelta is the length of the bracketing interval of the last
    // iteration.
    double x0 = min;
    double x1 = max;

    double y0 = computeObjectiveValue(x0);
    // return the first endpoint if it is good enough
    if (FastMath.abs(y0) <= functionValueAccuracy) {
        return x0;
    }

    // return the second endpoint if it is good enough
    double y1 = computeObjectiveValue(x1);
    if (FastMath.abs(y1) <= functionValueAccuracy) {
        return x1;
    }

    // Verify bracketing
    if (y0 * y1 >= 0) {
        throw new NoBracketingException(min, max, y0, y1);
    }

    final double absoluteAccuracy = getAbsoluteAccuracy();
    final double relativeAccuracy = getRelativeAccuracy();

    double x2 = x0;
    double y2 = y0;
    double oldDelta = x2 - x1;
    while (true) {
        if (FastMath.abs(y2) < FastMath.abs(y1)) {
            x0 = x1;
            x1 = x2;
            x2 = x0;
            y0 = y1;
            y1 = y2;
            y2 = y0;
        }
        if (FastMath.abs(y1) <= functionValueAccuracy) {
            return x1;
        }
        if (FastMath.abs(oldDelta) < FastMath.max(relativeAccuracy * FastMath.abs(x1),
                                                  absoluteAccuracy)) {
            return x1;
        }
        double delta;
        if (FastMath.abs(y1) > FastMath.abs(y0)) {
            // Function value increased in last iteration. Force bisection.
            delta = 0.5 * oldDelta;
        } else {
            delta = (x0 - x1) / (1 - y0 / y1);
            if (delta / oldDelta > 1) {
                // New approximation falls outside bracket.
                // Fall back to bisection.
                delta = 0.5 * oldDelta;
            }
        }
        x0 = x1;
        y0 = y1;
        x1 = x1 + delta;
        y1 = computeObjectiveValue(x1);
        if ((y1 > 0) == (y2 > 0)) {
            // New bracket is (x0,x1).
            x2 = x0;
            y2 = y0;
        }
        oldDelta = x2 - x1;
    }
}
 

/** 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
 */
public static double forceSide(final int maxEval, final UnivariateRealFunction f,
                               final BracketedUnivariateRealSolver<UnivariateRealFunction> bracketing,
                               final double baseRoot, final double min, final double max,
                               final AllowedSolution allowedSolution) {

    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);

}
 

/**
 * {@inheritDoc}
 */
@Override
protected double doSolve() {
    final double min = getMin();
    final double max = getMax();

    verifyInterval(min, max);

    final double relativeAccuracy = getRelativeAccuracy();
    final double absoluteAccuracy = getAbsoluteAccuracy();
    final double functionValueAccuracy = getFunctionValueAccuracy();

    // x2 is the last root approximation
    // x is the new approximation and new x2 for next round
    // x0 < x1 < x2 does not hold here

    double x0 = min;
    double y0 = computeObjectiveValue(x0);
    if (FastMath.abs(y0) < functionValueAccuracy) {
        return x0;
    }
    double x1 = max;
    double y1 = computeObjectiveValue(x1);
    if (FastMath.abs(y1) < functionValueAccuracy) {
        return x1;
    }

    if(y0 * y1 > 0) {
        throw new NoBracketingException(x0, x1, y0, y1);
    }

    double x2 = 0.5 * (x0 + x1);
    double y2 = computeObjectiveValue(x2);

    double oldx = Double.POSITIVE_INFINITY;
    while (true) {
        // quadratic interpolation through x0, x1, x2
        final double q = (x2 - x1) / (x1 - x0);
        final double a = q * (y2 - (1 + q) * y1 + q * y0);
        final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
        final double c = (1 + q) * y2;
        final double delta = b * b - 4 * a * c;
        double x;
        final double denominator;
        if (delta >= 0.0) {
            // choose a denominator larger in magnitude
            double dplus = b + FastMath.sqrt(delta);
            double dminus = b - FastMath.sqrt(delta);
            denominator = FastMath.abs(dplus) > FastMath.abs(dminus) ? dplus : dminus;
        } else {
            // take the modulus of (B +/- FastMath.sqrt(delta))
            denominator = FastMath.sqrt(b * b - delta);
        }
        if (denominator != 0) {
            x = x2 - 2.0 * c * (x2 - x1) / denominator;
            // perturb x if it exactly coincides with x1 or x2
            // the equality tests here are intentional
            while (x == x1 || x == x2) {
                x += absoluteAccuracy;
            }
        } else {
            // extremely rare case, get a random number to skip it
            x = min + FastMath.random() * (max - min);
            oldx = Double.POSITIVE_INFINITY;
        }
        final double y = computeObjectiveValue(x);

        // check for convergence
        final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
        if (FastMath.abs(x - oldx) <= tolerance ||
            FastMath.abs(y) <= functionValueAccuracy) {
            return x;
        }

        // prepare the next iteration
        x0 = x1;
        y0 = y1;
        x1 = x2;
        y1 = y2;
        x2 = x;
        y2 = y;
        oldx = x;
    }
}
 

/**
 * Method for implementing actual optimization algorithms in derived
 * classes.
 *
 * @return the root.
 * @throws TooManyEvaluationsException if the maximal number of evaluations
 * is exceeded.
 * @throws NoBracketingException if the initial search interval does not bracket
 * a root and the solver requires it.
 */
protected abstract double doSolve()
    throws TooManyEvaluationsException, NoBracketingException;
 

/**
 * Method for implementing actual optimization algorithms in derived
 * classes.
 *
 * @return the root.
 * @throws TooManyEvaluationsException if the maximal number of evaluations
 * is exceeded.
 * @throws NoBracketingException if the initial search interval does not bracket
 * a root and the solver requires it.
 */
protected abstract double doSolve()
    throws TooManyEvaluationsException, NoBracketingException;