Use "BrentSolver" implementation from "Commons Numbers".
This commit is contained in:
Gilles Sadowski
parent
49469f756d
commit
4b0f52c0dd
6
pom.xml
6
pom.xml
|
|
@ -413,6 +413,12 @@
|
|||
<version>1.0-SNAPSHOT</version>
|
||||
</dependency>
|
||||
|
||||
<dependency>
|
||||
<groupId>org.apache.commons</groupId>
|
||||
<artifactId>commons-numbers-rootfinder</artifactId>
|
||||
<version>1.0-SNAPSHOT</version>
|
||||
</dependency>
|
||||
|
||||
<dependency>
|
||||
<groupId>org.apache.commons</groupId>
|
||||
<artifactId>commons-rng-client-api</artifactId>
|
||||
|
|
|
|||
|
|
@ -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
|
||||
* <blockquote>
|
||||
* <b>Algorithms for Minimization Without Derivatives</b>,
|
||||
* <it>Richard P. Brent</it>,
|
||||
* Dover 0-486-41998-3
|
||||
* </blockquote>
|
||||
*
|
||||
* @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;
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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());
|
||||
}
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in New Issue