1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * https://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22 package org.hipparchus.analysis.solvers;
23
24
25 import org.hipparchus.analysis.UnivariateFunction;
26 import org.hipparchus.exception.LocalizedCoreFormats;
27 import org.hipparchus.exception.MathIllegalArgumentException;
28 import org.hipparchus.exception.MathIllegalStateException;
29 import org.hipparchus.util.FastMath;
30 import org.hipparchus.util.Precision;
31
32 /**
33 * This class implements a modification of the <a
34 * href="http://mathworld.wolfram.com/BrentsMethod.html"> Brent algorithm</a>.
35 * <p>
36 * The changes with respect to the original Brent algorithm are:
37 * <ul>
38 * <li>the returned value is chosen in the current interval according
39 * to user specified {@link AllowedSolution},</li>
40 * <li>the maximal order for the invert polynomial root search is
41 * user-specified instead of being invert quadratic only</li>
42 * </ul><p>
43 * The given interval must bracket the root.</p>
44 *
45 */
46 public class BracketingNthOrderBrentSolver
47 extends AbstractUnivariateSolver
48 implements BracketedUnivariateSolver<UnivariateFunction> {
49
50 /** Default absolute accuracy. */
51 private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
52
53 /** Default maximal order. */
54 private static final int DEFAULT_MAXIMAL_ORDER = 5;
55
56 /** Maximal aging triggering an attempt to balance the bracketing interval. */
57 private static final int MAXIMAL_AGING = 2;
58
59 /** Reduction factor for attempts to balance the bracketing interval. */
60 private static final double REDUCTION_FACTOR = 1.0 / 16.0;
61
62 /** Maximal order. */
63 private final int maximalOrder;
64
65 /** The kinds of solutions that the algorithm may accept. */
66 private AllowedSolution allowed;
67
68 /**
69 * Construct a solver with default accuracy and maximal order (1e-6 and 5 respectively)
70 */
71 public BracketingNthOrderBrentSolver() {
72 this(DEFAULT_ABSOLUTE_ACCURACY, DEFAULT_MAXIMAL_ORDER);
73 }
74
75 /**
76 * Construct a solver.
77 *
78 * @param absoluteAccuracy Absolute accuracy.
79 * @param maximalOrder maximal order.
80 * @exception MathIllegalArgumentException if maximal order is lower than 2
81 */
82 public BracketingNthOrderBrentSolver(final double absoluteAccuracy,
83 final int maximalOrder)
84 throws MathIllegalArgumentException {
85 super(absoluteAccuracy);
86 if (maximalOrder < 2) {
87 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
88 maximalOrder, 2);
89 }
90 this.maximalOrder = maximalOrder;
91 this.allowed = AllowedSolution.ANY_SIDE;
92 }
93
94 /**
95 * Construct a solver.
96 *
97 * @param relativeAccuracy Relative accuracy.
98 * @param absoluteAccuracy Absolute accuracy.
99 * @param maximalOrder maximal order.
100 * @exception MathIllegalArgumentException if maximal order is lower than 2
101 */
102 public BracketingNthOrderBrentSolver(final double relativeAccuracy,
103 final double absoluteAccuracy,
104 final int maximalOrder)
105 throws MathIllegalArgumentException {
106 super(relativeAccuracy, absoluteAccuracy);
107 if (maximalOrder < 2) {
108 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
109 maximalOrder, 2);
110 }
111 this.maximalOrder = maximalOrder;
112 this.allowed = AllowedSolution.ANY_SIDE;
113 }
114
115 /**
116 * Construct a solver.
117 *
118 * @param relativeAccuracy Relative accuracy.
119 * @param absoluteAccuracy Absolute accuracy.
120 * @param functionValueAccuracy Function value accuracy.
121 * @param maximalOrder maximal order.
122 * @exception MathIllegalArgumentException if maximal order is lower than 2
123 */
124 public BracketingNthOrderBrentSolver(final double relativeAccuracy,
125 final double absoluteAccuracy,
126 final double functionValueAccuracy,
127 final int maximalOrder)
128 throws MathIllegalArgumentException {
129 super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
130 if (maximalOrder < 2) {
131 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
132 maximalOrder, 2);
133 }
134 this.maximalOrder = maximalOrder;
135 this.allowed = AllowedSolution.ANY_SIDE;
136 }
137
138 /** Get the maximal order.
139 * @return maximal order
140 */
141 public int getMaximalOrder() {
142 return maximalOrder;
143 }
144
145 /**
146 * {@inheritDoc}
147 */
148 @Override
149 protected double doSolve() {
150 return doSolveInterval().getSide(allowed);
151 }
152
153 /**
154 * Find a root and return the containing interval.
155 *
156 * @return an interval containing the root such that both end points meet the
157 * convergence criteria.
158 */
159 protected Interval doSolveInterval() {
160 // prepare arrays with the first points
161 final double[] x = new double[maximalOrder + 1];
162 final double[] y = new double[maximalOrder + 1];
163 x[0] = getMin();
164 x[1] = getStartValue();
165 x[2] = getMax();
166 verifyInterval(x[0], x[2]);
167 if (x[1] < x[0] || x[2] < x[1]) {
168 throw new MathIllegalArgumentException(
169 LocalizedCoreFormats.START_POINT_NOT_IN_INTERVAL, x[1], x[0], x[2]);
170 }
171
172 // evaluate initial guess
173 y[1] = computeObjectiveValue(x[1]);
174 if (y[1] == 0.0) {
175 // return the initial guess if it is a perfect root.
176 return new Interval(x[1], y[1], x[1], y[1]);
177 }
178
179 // evaluate first endpoint
180 y[0] = computeObjectiveValue(x[0]);
181 if (y[0] == 0.0) {
182 // return the first endpoint if it is a perfect root.
183 return new Interval(x[0], y[0], x[0], y[0]);
184 }
185
186 int nbPoints;
187 int signChangeIndex;
188 if (y[0] * y[1] < 0) {
189
190 // reduce interval if it brackets the root
191 nbPoints = 2;
192 signChangeIndex = 1;
193
194 } else {
195
196 // evaluate second endpoint
197 y[2] = computeObjectiveValue(x[2]);
198 if (y[2] == 0.0) {
199 // return the second endpoint if it is a perfect root.
200 return new Interval(x[2], y[2], x[2], y[2]);
201 }
202
203 if (y[1] * y[2] < 0) {
204 // use all computed point as a start sampling array for solving
205 nbPoints = 3;
206 signChangeIndex = 2;
207 } else {
208 throw new MathIllegalArgumentException(LocalizedCoreFormats.NOT_BRACKETING_INTERVAL,
209 x[0], x[2], y[0], y[2]);
210 }
211
212 }
213
214 // prepare a work array for inverse polynomial interpolation
215 final double[] tmpX = new double[x.length];
216
217 // current tightest bracketing of the root
218 double xA = x[signChangeIndex - 1];
219 double yA = y[signChangeIndex - 1];
220 double absYA = FastMath.abs(yA);
221 int agingA = 0;
222 double xB = x[signChangeIndex];
223 double yB = y[signChangeIndex];
224 double absYB = FastMath.abs(yB);
225 int agingB = 0;
226
227 // search loop
228 while (true) {
229
230 // check convergence of bracketing interval
231 final double xTol = getAbsoluteAccuracy() +
232 getRelativeAccuracy() * FastMath.max(FastMath.abs(xA), FastMath.abs(xB));
233 if (xB - xA <= xTol ||
234 FastMath.max(absYA, absYB) < getFunctionValueAccuracy() ||
235 Precision.equals(xA, xB, 1)) {
236 return new Interval(xA, yA, xB, yB);
237 }
238
239 // target for the next evaluation point
240 double targetY;
241 if (agingA >= MAXIMAL_AGING) {
242 // we keep updating the high bracket, try to compensate this
243 final int p = agingA - MAXIMAL_AGING;
244 final double weightA = (1 << p) - 1;
245 final double weightB = p + 1;
246 targetY = (weightA * yA - weightB * REDUCTION_FACTOR * yB) / (weightA + weightB);
247 } else if (agingB >= MAXIMAL_AGING) {
248 // we keep updating the low bracket, try to compensate this
249 final int p = agingB - MAXIMAL_AGING;
250 final double weightA = p + 1;
251 final double weightB = (1 << p) - 1;
252 targetY = (weightB * yB - weightA * REDUCTION_FACTOR * yA) / (weightA + weightB);
253 } else {
254 // bracketing is balanced, try to find the root itself
255 targetY = 0;
256 }
257
258 // make a few attempts to guess a root,
259 double nextX;
260 int start = 0;
261 int end = nbPoints;
262 do {
263
264 // guess a value for current target, using inverse polynomial interpolation
265 System.arraycopy(x, start, tmpX, start, end - start);
266 nextX = guessX(targetY, tmpX, y, start, end);
267
268 if (!((nextX > xA) && (nextX < xB))) {
269 // the guessed root is not strictly inside of the tightest bracketing interval
270
271 // the guessed root is either not strictly inside the interval or it
272 // is a NaN (which occurs when some sampling points share the same y)
273 // we try again with a lower interpolation order
274 if (signChangeIndex - start >= end - signChangeIndex) {
275 // we have more points before the sign change, drop the lowest point
276 ++start;
277 } else {
278 // we have more points after sign change, drop the highest point
279 --end;
280 }
281
282 // we need to do one more attempt
283 nextX = Double.NaN;
284
285 }
286
287 } while (Double.isNaN(nextX) && (end - start > 1));
288
289 if (Double.isNaN(nextX)) {
290 // fall back to bisection
291 nextX = xA + 0.5 * (xB - xA);
292 start = signChangeIndex - 1;
293 end = signChangeIndex;
294 }
295
296 // evaluate the function at the guessed root
297 final double nextY = computeObjectiveValue(nextX);
298 if (nextY == 0.0 || FastMath.abs(nextY) < getFunctionValueAccuracy() && allowed == AllowedSolution.ANY_SIDE) {
299 // we have either:
300 // - an exact root, so we don't we don't need to bother about the allowed solutions setting
301 // - or an approximate root and we know allowed solutions setting if to retrieve the value closest to zero
302 return new Interval(nextX, nextY, nextX, nextY);
303 }
304
305 if ((nbPoints > 2) && (end - start != nbPoints)) {
306
307 // we have been forced to ignore some points to keep bracketing,
308 // they are probably too far from the root, drop them from now on
309 nbPoints = end - start;
310 System.arraycopy(x, start, x, 0, nbPoints);
311 System.arraycopy(y, start, y, 0, nbPoints);
312 signChangeIndex -= start;
313
314 } else if (nbPoints == x.length) {
315
316 // we have to drop one point in order to insert the new one
317 nbPoints--;
318
319 // keep the tightest bracketing interval as centered as possible
320 if (signChangeIndex >= (x.length + 1) / 2) {
321 // we drop the lowest point, we have to shift the arrays and the index
322 System.arraycopy(x, 1, x, 0, nbPoints);
323 System.arraycopy(y, 1, y, 0, nbPoints);
324 --signChangeIndex;
325 }
326
327 }
328
329 // insert the last computed point
330 //(by construction, we know it lies inside the tightest bracketing interval)
331 System.arraycopy(x, signChangeIndex, x, signChangeIndex + 1, nbPoints - signChangeIndex);
332 x[signChangeIndex] = nextX;
333 System.arraycopy(y, signChangeIndex, y, signChangeIndex + 1, nbPoints - signChangeIndex);
334 y[signChangeIndex] = nextY;
335 ++nbPoints;
336
337 // update the bracketing interval
338 if (nextY * yA <= 0) {
339 // the sign change occurs before the inserted point
340 xB = nextX;
341 yB = nextY;
342 absYB = FastMath.abs(yB);
343 ++agingA;
344 agingB = 0;
345 } else {
346 // the sign change occurs after the inserted point
347 xA = nextX;
348 yA = nextY;
349 absYA = FastMath.abs(yA);
350 agingA = 0;
351 ++agingB;
352
353 // update the sign change index
354 signChangeIndex++;
355
356 }
357
358 }
359
360 }
361
362 /** Guess an x value by n<sup>th</sup> order inverse polynomial interpolation.
363 * <p>
364 * The x value is guessed by evaluating polynomial Q(y) at y = targetY, where Q
365 * is built such that for all considered points (x<sub>i</sub>, y<sub>i</sub>),
366 * Q(y<sub>i</sub>) = x<sub>i</sub>.
367 * </p>
368 * @param targetY target value for y
369 * @param x reference points abscissas for interpolation,
370 * note that this array <em>is</em> modified during computation
371 * @param y reference points ordinates for interpolation
372 * @param start start index of the points to consider (inclusive)
373 * @param end end index of the points to consider (exclusive)
374 * @return guessed root (will be a NaN if two points share the same y)
375 */
376 private double guessX(final double targetY, final double[] x, final double[] y,
377 final int start, final int end) {
378
379 // compute Q Newton coefficients by divided differences
380 for (int i = start; i < end - 1; ++i) {
381 final int delta = i + 1 - start;
382 for (int j = end - 1; j > i; --j) {
383 x[j] = (x[j] - x[j-1]) / (y[j] - y[j - delta]);
384 }
385 }
386
387 // evaluate Q(targetY)
388 double x0 = 0;
389 for (int j = end - 1; j >= start; --j) {
390 x0 = x[j] + x0 * (targetY - y[j]);
391 }
392
393 return x0;
394
395 }
396
397 /** {@inheritDoc} */
398 @Override
399 public double solve(int maxEval, UnivariateFunction f, double min,
400 double max, AllowedSolution allowedSolution)
401 throws MathIllegalArgumentException, MathIllegalStateException {
402 this.allowed = allowedSolution;
403 return super.solve(maxEval, f, min, max);
404 }
405
406 /** {@inheritDoc} */
407 @Override
408 public double solve(int maxEval, UnivariateFunction f, double min,
409 double max, double startValue,
410 AllowedSolution allowedSolution)
411 throws MathIllegalArgumentException, MathIllegalStateException {
412 this.allowed = allowedSolution;
413 return super.solve(maxEval, f, min, max, startValue);
414
415 }
416
417 /** {@inheritDoc} */
418 @Override
419 public Interval solveInterval(final int maxEval,
420 final UnivariateFunction f,
421 final double min,
422 final double max,
423 final double startValue)
424 throws MathIllegalArgumentException, MathIllegalStateException {
425 setup(maxEval, f, min, max, startValue);
426 this.allowed = null;
427 return doSolveInterval();
428 }
429
430 }