From 4b0f52c0dd9df564d58140e1d717a2da218f0139 Mon Sep 17 00:00:00 2001 From: Gilles Sadowski Date: Sun, 9 Jun 2019 09:56:52 +0200 Subject: [PATCH] Use "BrentSolver" implementation from "Commons Numbers". --- pom.xml | 6 + .../math4/analysis/solvers/BrentSolver.java | 155 ++---------------- .../analysis/solvers/BrentSolverTest.java | 2 +- 3 files changed, 22 insertions(+), 141 deletions(-) diff --git a/pom.xml b/pom.xml index fcee9462c..e89f2cae1 100644 --- a/pom.xml +++ b/pom.xml @@ -413,6 +413,12 @@ 1.0-SNAPSHOT + + org.apache.commons + commons-numbers-rootfinder + 1.0-SNAPSHOT + + org.apache.commons commons-rng-client-api diff --git a/src/main/java/org/apache/commons/math4/analysis/solvers/BrentSolver.java b/src/main/java/org/apache/commons/math4/analysis/solvers/BrentSolver.java index 9e1e55193..1e0b3b844 100644 --- a/src/main/java/org/apache/commons/math4/analysis/solvers/BrentSolver.java +++ b/src/main/java/org/apache/commons/math4/analysis/solvers/BrentSolver.java @@ -94,150 +94,25 @@ public class BrentSolver extends AbstractUnivariateSolver { throws NoBracketingException, TooManyEvaluationsException, NumberIsTooLargeException { - double min = getMin(); - double max = getMax(); + final double min = getMin(); + final double max = getMax(); final double initial = getStartValue(); - final double functionValueAccuracy = getFunctionValueAccuracy(); - verifySequence(min, initial, max); + final org.apache.commons.numbers.rootfinder.BrentSolver rf = + new org.apache.commons.numbers.rootfinder.BrentSolver(getRelativeAccuracy(), + getAbsoluteAccuracy(), + getFunctionValueAccuracy()); - // Return the initial guess if it is good enough. - double yInitial = computeObjectiveValue(initial); - if (FastMath.abs(yInitial) <= functionValueAccuracy) { - return initial; + double root = Double.NaN; + try { + root = rf.findRoot(arg -> computeObjectiveValue(arg), + min, initial, max); + } catch (IllegalArgumentException e) { + // Redundant calls in order to throw the expected exceptions. + verifySequence(min, initial, max); + verifyBracketing(min, max); } - // 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); - } - - /** - * Search for a zero inside the provided interval. - * This implementation is based on the algorithm described at page 58 of - * the book - *
- * Algorithms for Minimization Without Derivatives, - * Richard P. Brent, - * Dover 0-486-41998-3 - *
- * - * @param lo Lower bound of the search interval. - * @param hi Higher bound of the search interval. - * @param fLo Function value at the lower bound of the search interval. - * @param fHi Function value at the higher bound of the search interval. - * @return the value where the function is zero. - */ - private double brent(double lo, double hi, - double fLo, double fHi) { - double a = lo; - double fa = fLo; - double b = hi; - double fb = fHi; - double c = a; - double fc = fa; - double d = b - a; - double e = d; - - final double t = getAbsoluteAccuracy(); - final double eps = getRelativeAccuracy(); - - while (true) { - if (FastMath.abs(fc) < FastMath.abs(fb)) { - a = b; - b = c; - c = a; - fa = fb; - fb = fc; - fc = fa; - } - - final double tol = 2 * eps * FastMath.abs(b) + t; - final double m = 0.5 * (c - b); - - if (FastMath.abs(m) <= tol || - Precision.equals(fb, 0)) { - return b; - } - if (FastMath.abs(e) < tol || - FastMath.abs(fa) <= FastMath.abs(fb)) { - // Force bisection. - d = m; - e = d; - } else { - double s = fb / fa; - double p; - double q; - // The equality test (a == c) is intentional, - // it is part of the original Brent's method and - // it should NOT be replaced by proximity test. - if (a == c) { - // Linear interpolation. - p = 2 * m * s; - q = 1 - s; - } else { - // Inverse quadratic interpolation. - q = fa / fc; - final double r = fb / fc; - p = s * (2 * m * q * (q - r) - (b - a) * (r - 1)); - q = (q - 1) * (r - 1) * (s - 1); - } - if (p > 0) { - q = -q; - } else { - p = -p; - } - s = e; - e = d; - if (p >= 1.5 * m * q - FastMath.abs(tol * q) || - p >= FastMath.abs(0.5 * s * q)) { - // Inverse quadratic interpolation gives a value - // in the wrong direction, or progress is slow. - // Fall back to bisection. - d = m; - e = d; - } else { - d = p / q; - } - } - a = b; - fa = fb; - - if (FastMath.abs(d) > tol) { - b += d; - } else if (m > 0) { - b += tol; - } else { - b -= tol; - } - fb = computeObjectiveValue(b); - if ((fb > 0 && fc > 0) || - (fb <= 0 && fc <= 0)) { - c = a; - fc = fa; - d = b - a; - e = d; - } - } + return root; } } diff --git a/src/test/java/org/apache/commons/math4/analysis/solvers/BrentSolverTest.java b/src/test/java/org/apache/commons/math4/analysis/solvers/BrentSolverTest.java index 3da021ccb..cc17d0c73 100644 --- a/src/test/java/org/apache/commons/math4/analysis/solvers/BrentSolverTest.java +++ b/src/test/java/org/apache/commons/math4/analysis/solvers/BrentSolverTest.java @@ -265,6 +265,6 @@ public final class BrentSolverTest { BrentSolver solver = new BrentSolver(); final double result = solver.solve(99, f, 1, 1e30, 1 + 1e-10); - Assert.assertEquals(804.93558250, result, 1e-8); + Assert.assertEquals(804.93558250, result, solver.getAbsoluteAccuracy()); } }