001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017 package org.apache.commons.math3.analysis.solvers;
018
019
020 import org.apache.commons.math3.analysis.UnivariateFunction;
021 import org.apache.commons.math3.exception.MathInternalError;
022 import org.apache.commons.math3.exception.NoBracketingException;
023 import org.apache.commons.math3.exception.NumberIsTooSmallException;
024 import org.apache.commons.math3.exception.NumberIsTooLargeException;
025 import org.apache.commons.math3.exception.TooManyEvaluationsException;
026 import org.apache.commons.math3.util.FastMath;
027 import org.apache.commons.math3.util.Precision;
028
029 /**
030 * This class implements a modification of the <a
031 * href="http://mathworld.wolfram.com/BrentsMethod.html"> Brent algorithm</a>.
032 * <p>
033 * The changes with respect to the original Brent algorithm are:
034 * <ul>
035 * <li>the returned value is chosen in the current interval according
036 * to user specified {@link AllowedSolution},</li>
037 * <li>the maximal order for the invert polynomial root search is
038 * user-specified instead of being invert quadratic only</li>
039 * </ul>
040 * </p>
041 * The given interval must bracket the root.
042 *
043 * @version $Id: BracketingNthOrderBrentSolver.java 1379560 2012-08-31 19:40:30Z erans $
044 */
045 public class BracketingNthOrderBrentSolver
046 extends AbstractUnivariateSolver
047 implements BracketedUnivariateSolver<UnivariateFunction> {
048
049 /** Default absolute accuracy. */
050 private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
051
052 /** Default maximal order. */
053 private static final int DEFAULT_MAXIMAL_ORDER = 5;
054
055 /** Maximal aging triggering an attempt to balance the bracketing interval. */
056 private static final int MAXIMAL_AGING = 2;
057
058 /** Reduction factor for attempts to balance the bracketing interval. */
059 private static final double REDUCTION_FACTOR = 1.0 / 16.0;
060
061 /** Maximal order. */
062 private final int maximalOrder;
063
064 /** The kinds of solutions that the algorithm may accept. */
065 private AllowedSolution allowed;
066
067 /**
068 * Construct a solver with default accuracy and maximal order (1e-6 and 5 respectively)
069 */
070 public BracketingNthOrderBrentSolver() {
071 this(DEFAULT_ABSOLUTE_ACCURACY, DEFAULT_MAXIMAL_ORDER);
072 }
073
074 /**
075 * Construct a solver.
076 *
077 * @param absoluteAccuracy Absolute accuracy.
078 * @param maximalOrder maximal order.
079 * @exception NumberIsTooSmallException if maximal order is lower than 2
080 */
081 public BracketingNthOrderBrentSolver(final double absoluteAccuracy,
082 final int maximalOrder)
083 throws NumberIsTooSmallException {
084 super(absoluteAccuracy);
085 if (maximalOrder < 2) {
086 throw new NumberIsTooSmallException(maximalOrder, 2, true);
087 }
088 this.maximalOrder = maximalOrder;
089 this.allowed = AllowedSolution.ANY_SIDE;
090 }
091
092 /**
093 * Construct a solver.
094 *
095 * @param relativeAccuracy Relative accuracy.
096 * @param absoluteAccuracy Absolute accuracy.
097 * @param maximalOrder maximal order.
098 * @exception NumberIsTooSmallException if maximal order is lower than 2
099 */
100 public BracketingNthOrderBrentSolver(final double relativeAccuracy,
101 final double absoluteAccuracy,
102 final int maximalOrder)
103 throws NumberIsTooSmallException {
104 super(relativeAccuracy, absoluteAccuracy);
105 if (maximalOrder < 2) {
106 throw new NumberIsTooSmallException(maximalOrder, 2, true);
107 }
108 this.maximalOrder = maximalOrder;
109 this.allowed = AllowedSolution.ANY_SIDE;
110 }
111
112 /**
113 * Construct a solver.
114 *
115 * @param relativeAccuracy Relative accuracy.
116 * @param absoluteAccuracy Absolute accuracy.
117 * @param functionValueAccuracy Function value accuracy.
118 * @param maximalOrder maximal order.
119 * @exception NumberIsTooSmallException if maximal order is lower than 2
120 */
121 public BracketingNthOrderBrentSolver(final double relativeAccuracy,
122 final double absoluteAccuracy,
123 final double functionValueAccuracy,
124 final int maximalOrder)
125 throws NumberIsTooSmallException {
126 super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
127 if (maximalOrder < 2) {
128 throw new NumberIsTooSmallException(maximalOrder, 2, true);
129 }
130 this.maximalOrder = maximalOrder;
131 this.allowed = AllowedSolution.ANY_SIDE;
132 }
133
134 /** Get the maximal order.
135 * @return maximal order
136 */
137 public int getMaximalOrder() {
138 return maximalOrder;
139 }
140
141 /**
142 * {@inheritDoc}
143 */
144 @Override
145 protected double doSolve()
146 throws TooManyEvaluationsException,
147 NumberIsTooLargeException,
148 NoBracketingException {
149 // prepare arrays with the first points
150 final double[] x = new double[maximalOrder + 1];
151 final double[] y = new double[maximalOrder + 1];
152 x[0] = getMin();
153 x[1] = getStartValue();
154 x[2] = getMax();
155 verifySequence(x[0], x[1], x[2]);
156
157 // evaluate initial guess
158 y[1] = computeObjectiveValue(x[1]);
159 if (Precision.equals(y[1], 0.0, 1)) {
160 // return the initial guess if it is a perfect root.
161 return x[1];
162 }
163
164 // evaluate first endpoint
165 y[0] = computeObjectiveValue(x[0]);
166 if (Precision.equals(y[0], 0.0, 1)) {
167 // return the first endpoint if it is a perfect root.
168 return x[0];
169 }
170
171 int nbPoints;
172 int signChangeIndex;
173 if (y[0] * y[1] < 0) {
174
175 // reduce interval if it brackets the root
176 nbPoints = 2;
177 signChangeIndex = 1;
178
179 } else {
180
181 // evaluate second endpoint
182 y[2] = computeObjectiveValue(x[2]);
183 if (Precision.equals(y[2], 0.0, 1)) {
184 // return the second endpoint if it is a perfect root.
185 return x[2];
186 }
187
188 if (y[1] * y[2] < 0) {
189 // use all computed point as a start sampling array for solving
190 nbPoints = 3;
191 signChangeIndex = 2;
192 } else {
193 throw new NoBracketingException(x[0], x[2], y[0], y[2]);
194 }
195
196 }
197
198 // prepare a work array for inverse polynomial interpolation
199 final double[] tmpX = new double[x.length];
200
201 // current tightest bracketing of the root
202 double xA = x[signChangeIndex - 1];
203 double yA = y[signChangeIndex - 1];
204 double absYA = FastMath.abs(yA);
205 int agingA = 0;
206 double xB = x[signChangeIndex];
207 double yB = y[signChangeIndex];
208 double absYB = FastMath.abs(yB);
209 int agingB = 0;
210
211 // search loop
212 while (true) {
213
214 // check convergence of bracketing interval
215 final double xTol = getAbsoluteAccuracy() +
216 getRelativeAccuracy() * FastMath.max(FastMath.abs(xA), FastMath.abs(xB));
217 if (((xB - xA) <= xTol) || (FastMath.max(absYA, absYB) < getFunctionValueAccuracy())) {
218 switch (allowed) {
219 case ANY_SIDE :
220 return absYA < absYB ? xA : xB;
221 case LEFT_SIDE :
222 return xA;
223 case RIGHT_SIDE :
224 return xB;
225 case BELOW_SIDE :
226 return (yA <= 0) ? xA : xB;
227 case ABOVE_SIDE :
228 return (yA < 0) ? xB : xA;
229 default :
230 // this should never happen
231 throw new MathInternalError();
232 }
233 }
234
235 // target for the next evaluation point
236 double targetY;
237 if (agingA >= MAXIMAL_AGING) {
238 // we keep updating the high bracket, try to compensate this
239 final int p = agingA - MAXIMAL_AGING;
240 final double weightA = (1 << p) - 1;
241 final double weightB = p + 1;
242 targetY = (weightA * yA - weightB * REDUCTION_FACTOR * yB) / (weightA + weightB);
243 } else if (agingB >= MAXIMAL_AGING) {
244 // we keep updating the low bracket, try to compensate this
245 final int p = agingB - MAXIMAL_AGING;
246 final double weightA = p + 1;
247 final double weightB = (1 << p) - 1;
248 targetY = (weightB * yB - weightA * REDUCTION_FACTOR * yA) / (weightA + weightB);
249 } else {
250 // bracketing is balanced, try to find the root itself
251 targetY = 0;
252 }
253
254 // make a few attempts to guess a root,
255 double nextX;
256 int start = 0;
257 int end = nbPoints;
258 do {
259
260 // guess a value for current target, using inverse polynomial interpolation
261 System.arraycopy(x, start, tmpX, start, end - start);
262 nextX = guessX(targetY, tmpX, y, start, end);
263
264 if (!((nextX > xA) && (nextX < xB))) {
265 // the guessed root is not strictly inside of the tightest bracketing interval
266
267 // the guessed root is either not strictly inside the interval or it
268 // is a NaN (which occurs when some sampling points share the same y)
269 // we try again with a lower interpolation order
270 if (signChangeIndex - start >= end - signChangeIndex) {
271 // we have more points before the sign change, drop the lowest point
272 ++start;
273 } else {
274 // we have more points after sign change, drop the highest point
275 --end;
276 }
277
278 // we need to do one more attempt
279 nextX = Double.NaN;
280
281 }
282
283 } while (Double.isNaN(nextX) && (end - start > 1));
284
285 if (Double.isNaN(nextX)) {
286 // fall back to bisection
287 nextX = xA + 0.5 * (xB - xA);
288 start = signChangeIndex - 1;
289 end = signChangeIndex;
290 }
291
292 // evaluate the function at the guessed root
293 final double nextY = computeObjectiveValue(nextX);
294 if (Precision.equals(nextY, 0.0, 1)) {
295 // we have found an exact root, since it is not an approximation
296 // we don't need to bother about the allowed solutions setting
297 return nextX;
298 }
299
300 if ((nbPoints > 2) && (end - start != nbPoints)) {
301
302 // we have been forced to ignore some points to keep bracketing,
303 // they are probably too far from the root, drop them from now on
304 nbPoints = end - start;
305 System.arraycopy(x, start, x, 0, nbPoints);
306 System.arraycopy(y, start, y, 0, nbPoints);
307 signChangeIndex -= start;
308
309 } else if (nbPoints == x.length) {
310
311 // we have to drop one point in order to insert the new one
312 nbPoints--;
313
314 // keep the tightest bracketing interval as centered as possible
315 if (signChangeIndex >= (x.length + 1) / 2) {
316 // we drop the lowest point, we have to shift the arrays and the index
317 System.arraycopy(x, 1, x, 0, nbPoints);
318 System.arraycopy(y, 1, y, 0, nbPoints);
319 --signChangeIndex;
320 }
321
322 }
323
324 // insert the last computed point
325 //(by construction, we know it lies inside the tightest bracketing interval)
326 System.arraycopy(x, signChangeIndex, x, signChangeIndex + 1, nbPoints - signChangeIndex);
327 x[signChangeIndex] = nextX;
328 System.arraycopy(y, signChangeIndex, y, signChangeIndex + 1, nbPoints - signChangeIndex);
329 y[signChangeIndex] = nextY;
330 ++nbPoints;
331
332 // update the bracketing interval
333 if (nextY * yA <= 0) {
334 // the sign change occurs before the inserted point
335 xB = nextX;
336 yB = nextY;
337 absYB = FastMath.abs(yB);
338 ++agingA;
339 agingB = 0;
340 } else {
341 // the sign change occurs after the inserted point
342 xA = nextX;
343 yA = nextY;
344 absYA = FastMath.abs(yA);
345 agingA = 0;
346 ++agingB;
347
348 // update the sign change index
349 signChangeIndex++;
350
351 }
352
353 }
354
355 }
356
357 /** Guess an x value by n<sup>th</sup> order inverse polynomial interpolation.
358 * <p>
359 * The x value is guessed by evaluating polynomial Q(y) at y = targetY, where Q
360 * is built such that for all considered points (x<sub>i</sub>, y<sub>i</sub>),
361 * Q(y<sub>i</sub>) = x<sub>i</sub>.
362 * </p>
363 * @param targetY target value for y
364 * @param x reference points abscissas for interpolation,
365 * note that this array <em>is</em> modified during computation
366 * @param y reference points ordinates for interpolation
367 * @param start start index of the points to consider (inclusive)
368 * @param end end index of the points to consider (exclusive)
369 * @return guessed root (will be a NaN if two points share the same y)
370 */
371 private double guessX(final double targetY, final double[] x, final double[] y,
372 final int start, final int end) {
373
374 // compute Q Newton coefficients by divided differences
375 for (int i = start; i < end - 1; ++i) {
376 final int delta = i + 1 - start;
377 for (int j = end - 1; j > i; --j) {
378 x[j] = (x[j] - x[j-1]) / (y[j] - y[j - delta]);
379 }
380 }
381
382 // evaluate Q(targetY)
383 double x0 = 0;
384 for (int j = end - 1; j >= start; --j) {
385 x0 = x[j] + x0 * (targetY - y[j]);
386 }
387
388 return x0;
389
390 }
391
392 /** {@inheritDoc} */
393 public double solve(int maxEval, UnivariateFunction f, double min,
394 double max, AllowedSolution allowedSolution)
395 throws TooManyEvaluationsException,
396 NumberIsTooLargeException,
397 NoBracketingException {
398 this.allowed = allowedSolution;
399 return super.solve(maxEval, f, min, max);
400 }
401
402 /** {@inheritDoc} */
403 public double solve(int maxEval, UnivariateFunction f, double min,
404 double max, double startValue,
405 AllowedSolution allowedSolution)
406 throws TooManyEvaluationsException,
407 NumberIsTooLargeException,
408 NoBracketingException {
409 this.allowed = allowedSolution;
410 return super.solve(maxEval, f, min, max, startValue);
411 }
412
413 }