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.dfp;
018
019
020 import org.apache.commons.math3.analysis.solvers.AllowedSolution;
021 import org.apache.commons.math3.exception.MathInternalError;
022 import org.apache.commons.math3.exception.NoBracketingException;
023 import org.apache.commons.math3.exception.NullArgumentException;
024 import org.apache.commons.math3.exception.NumberIsTooSmallException;
025 import org.apache.commons.math3.util.Incrementor;
026 import org.apache.commons.math3.util.MathUtils;
027
028 /**
029 * This class implements a modification of the <a
030 * href="http://mathworld.wolfram.com/BrentsMethod.html"> Brent algorithm</a>.
031 * <p>
032 * The changes with respect to the original Brent algorithm are:
033 * <ul>
034 * <li>the returned value is chosen in the current interval according
035 * to user specified {@link AllowedSolution},</li>
036 * <li>the maximal order for the invert polynomial root search is
037 * user-specified instead of being invert quadratic only</li>
038 * </ul>
039 * </p>
040 * The given interval must bracket the root.
041 *
042 * @version $Id: BracketingNthOrderBrentSolverDFP.java 1416643 2012-12-03 19:37:14Z tn $
043 */
044 public class BracketingNthOrderBrentSolverDFP {
045
046 /** Maximal aging triggering an attempt to balance the bracketing interval. */
047 private static final int MAXIMAL_AGING = 2;
048
049 /** Maximal order. */
050 private final int maximalOrder;
051
052 /** Function value accuracy. */
053 private final Dfp functionValueAccuracy;
054
055 /** Absolute accuracy. */
056 private final Dfp absoluteAccuracy;
057
058 /** Relative accuracy. */
059 private final Dfp relativeAccuracy;
060
061 /** Evaluations counter. */
062 private final Incrementor evaluations = new Incrementor();
063
064 /**
065 * Construct a solver.
066 *
067 * @param relativeAccuracy Relative accuracy.
068 * @param absoluteAccuracy Absolute accuracy.
069 * @param functionValueAccuracy Function value accuracy.
070 * @param maximalOrder maximal order.
071 * @exception NumberIsTooSmallException if maximal order is lower than 2
072 */
073 public BracketingNthOrderBrentSolverDFP(final Dfp relativeAccuracy,
074 final Dfp absoluteAccuracy,
075 final Dfp functionValueAccuracy,
076 final int maximalOrder)
077 throws NumberIsTooSmallException {
078 if (maximalOrder < 2) {
079 throw new NumberIsTooSmallException(maximalOrder, 2, true);
080 }
081 this.maximalOrder = maximalOrder;
082 this.absoluteAccuracy = absoluteAccuracy;
083 this.relativeAccuracy = relativeAccuracy;
084 this.functionValueAccuracy = functionValueAccuracy;
085 }
086
087 /** Get the maximal order.
088 * @return maximal order
089 */
090 public int getMaximalOrder() {
091 return maximalOrder;
092 }
093
094 /**
095 * Get the maximal number of function evaluations.
096 *
097 * @return the maximal number of function evaluations.
098 */
099 public int getMaxEvaluations() {
100 return evaluations.getMaximalCount();
101 }
102
103 /**
104 * Get the number of evaluations of the objective function.
105 * The number of evaluations corresponds to the last call to the
106 * {@code optimize} method. It is 0 if the method has not been
107 * called yet.
108 *
109 * @return the number of evaluations of the objective function.
110 */
111 public int getEvaluations() {
112 return evaluations.getCount();
113 }
114
115 /**
116 * Get the absolute accuracy.
117 * @return absolute accuracy
118 */
119 public Dfp getAbsoluteAccuracy() {
120 return absoluteAccuracy;
121 }
122
123 /**
124 * Get the relative accuracy.
125 * @return relative accuracy
126 */
127 public Dfp getRelativeAccuracy() {
128 return relativeAccuracy;
129 }
130
131 /**
132 * Get the function accuracy.
133 * @return function accuracy
134 */
135 public Dfp getFunctionValueAccuracy() {
136 return functionValueAccuracy;
137 }
138
139 /**
140 * Solve for a zero in the given interval.
141 * A solver may require that the interval brackets a single zero root.
142 * Solvers that do require bracketing should be able to handle the case
143 * where one of the endpoints is itself a root.
144 *
145 * @param maxEval Maximum number of evaluations.
146 * @param f Function to solve.
147 * @param min Lower bound for the interval.
148 * @param max Upper bound for the interval.
149 * @param allowedSolution The kind of solutions that the root-finding algorithm may
150 * accept as solutions.
151 * @return a value where the function is zero.
152 * @exception NullArgumentException if f is null.
153 * @exception NoBracketingException if root cannot be bracketed
154 */
155 public Dfp solve(final int maxEval, final UnivariateDfpFunction f,
156 final Dfp min, final Dfp max, final AllowedSolution allowedSolution)
157 throws NullArgumentException, NoBracketingException {
158 return solve(maxEval, f, min, max, min.add(max).divide(2), allowedSolution);
159 }
160
161 /**
162 * Solve for a zero in the given interval, start at {@code startValue}.
163 * A solver may require that the interval brackets a single zero root.
164 * Solvers that do require bracketing should be able to handle the case
165 * where one of the endpoints is itself a root.
166 *
167 * @param maxEval Maximum number of evaluations.
168 * @param f Function to solve.
169 * @param min Lower bound for the interval.
170 * @param max Upper bound for the interval.
171 * @param startValue Start value to use.
172 * @param allowedSolution The kind of solutions that the root-finding algorithm may
173 * accept as solutions.
174 * @return a value where the function is zero.
175 * @exception NullArgumentException if f is null.
176 * @exception NoBracketingException if root cannot be bracketed
177 */
178 public Dfp solve(final int maxEval, final UnivariateDfpFunction f,
179 final Dfp min, final Dfp max, final Dfp startValue,
180 final AllowedSolution allowedSolution)
181 throws NullArgumentException, NoBracketingException {
182
183 // Checks.
184 MathUtils.checkNotNull(f);
185
186 // Reset.
187 evaluations.setMaximalCount(maxEval);
188 evaluations.resetCount();
189 Dfp zero = startValue.getZero();
190 Dfp nan = zero.newInstance((byte) 1, Dfp.QNAN);
191
192 // prepare arrays with the first points
193 final Dfp[] x = new Dfp[maximalOrder + 1];
194 final Dfp[] y = new Dfp[maximalOrder + 1];
195 x[0] = min;
196 x[1] = startValue;
197 x[2] = max;
198
199 // evaluate initial guess
200 evaluations.incrementCount();
201 y[1] = f.value(x[1]);
202 if (y[1].isZero()) {
203 // return the initial guess if it is a perfect root.
204 return x[1];
205 }
206
207 // evaluate first endpoint
208 evaluations.incrementCount();
209 y[0] = f.value(x[0]);
210 if (y[0].isZero()) {
211 // return the first endpoint if it is a perfect root.
212 return x[0];
213 }
214
215 int nbPoints;
216 int signChangeIndex;
217 if (y[0].multiply(y[1]).negativeOrNull()) {
218
219 // reduce interval if it brackets the root
220 nbPoints = 2;
221 signChangeIndex = 1;
222
223 } else {
224
225 // evaluate second endpoint
226 evaluations.incrementCount();
227 y[2] = f.value(x[2]);
228 if (y[2].isZero()) {
229 // return the second endpoint if it is a perfect root.
230 return x[2];
231 }
232
233 if (y[1].multiply(y[2]).negativeOrNull()) {
234 // use all computed point as a start sampling array for solving
235 nbPoints = 3;
236 signChangeIndex = 2;
237 } else {
238 throw new NoBracketingException(x[0].toDouble(), x[2].toDouble(),
239 y[0].toDouble(), y[2].toDouble());
240 }
241
242 }
243
244 // prepare a work array for inverse polynomial interpolation
245 final Dfp[] tmpX = new Dfp[x.length];
246
247 // current tightest bracketing of the root
248 Dfp xA = x[signChangeIndex - 1];
249 Dfp yA = y[signChangeIndex - 1];
250 Dfp absXA = xA.abs();
251 Dfp absYA = yA.abs();
252 int agingA = 0;
253 Dfp xB = x[signChangeIndex];
254 Dfp yB = y[signChangeIndex];
255 Dfp absXB = xB.abs();
256 Dfp absYB = yB.abs();
257 int agingB = 0;
258
259 // search loop
260 while (true) {
261
262 // check convergence of bracketing interval
263 Dfp maxX = absXA.lessThan(absXB) ? absXB : absXA;
264 Dfp maxY = absYA.lessThan(absYB) ? absYB : absYA;
265 final Dfp xTol = absoluteAccuracy.add(relativeAccuracy.multiply(maxX));
266 if (xB.subtract(xA).subtract(xTol).negativeOrNull() ||
267 maxY.lessThan(functionValueAccuracy)) {
268 switch (allowedSolution) {
269 case ANY_SIDE :
270 return absYA.lessThan(absYB) ? xA : xB;
271 case LEFT_SIDE :
272 return xA;
273 case RIGHT_SIDE :
274 return xB;
275 case BELOW_SIDE :
276 return yA.lessThan(zero) ? xA : xB;
277 case ABOVE_SIDE :
278 return yA.lessThan(zero) ? xB : xA;
279 default :
280 // this should never happen
281 throw new MathInternalError(null);
282 }
283 }
284
285 // target for the next evaluation point
286 Dfp targetY;
287 if (agingA >= MAXIMAL_AGING) {
288 // we keep updating the high bracket, try to compensate this
289 targetY = yB.divide(16).negate();
290 } else if (agingB >= MAXIMAL_AGING) {
291 // we keep updating the low bracket, try to compensate this
292 targetY = yA.divide(16).negate();
293 } else {
294 // bracketing is balanced, try to find the root itself
295 targetY = zero;
296 }
297
298 // make a few attempts to guess a root,
299 Dfp nextX;
300 int start = 0;
301 int end = nbPoints;
302 do {
303
304 // guess a value for current target, using inverse polynomial interpolation
305 System.arraycopy(x, start, tmpX, start, end - start);
306 nextX = guessX(targetY, tmpX, y, start, end);
307
308 if (!(nextX.greaterThan(xA) && nextX.lessThan(xB))) {
309 // the guessed root is not strictly inside of the tightest bracketing interval
310
311 // the guessed root is either not strictly inside the interval or it
312 // is a NaN (which occurs when some sampling points share the same y)
313 // we try again with a lower interpolation order
314 if (signChangeIndex - start >= end - signChangeIndex) {
315 // we have more points before the sign change, drop the lowest point
316 ++start;
317 } else {
318 // we have more points after sign change, drop the highest point
319 --end;
320 }
321
322 // we need to do one more attempt
323 nextX = nan;
324
325 }
326
327 } while (nextX.isNaN() && (end - start > 1));
328
329 if (nextX.isNaN()) {
330 // fall back to bisection
331 nextX = xA.add(xB.subtract(xA).divide(2));
332 start = signChangeIndex - 1;
333 end = signChangeIndex;
334 }
335
336 // evaluate the function at the guessed root
337 evaluations.incrementCount();
338 final Dfp nextY = f.value(nextX);
339 if (nextY.isZero()) {
340 // we have found an exact root, since it is not an approximation
341 // we don't need to bother about the allowed solutions setting
342 return nextX;
343 }
344
345 if ((nbPoints > 2) && (end - start != nbPoints)) {
346
347 // we have been forced to ignore some points to keep bracketing,
348 // they are probably too far from the root, drop them from now on
349 nbPoints = end - start;
350 System.arraycopy(x, start, x, 0, nbPoints);
351 System.arraycopy(y, start, y, 0, nbPoints);
352 signChangeIndex -= start;
353
354 } else if (nbPoints == x.length) {
355
356 // we have to drop one point in order to insert the new one
357 nbPoints--;
358
359 // keep the tightest bracketing interval as centered as possible
360 if (signChangeIndex >= (x.length + 1) / 2) {
361 // we drop the lowest point, we have to shift the arrays and the index
362 System.arraycopy(x, 1, x, 0, nbPoints);
363 System.arraycopy(y, 1, y, 0, nbPoints);
364 --signChangeIndex;
365 }
366
367 }
368
369 // insert the last computed point
370 //(by construction, we know it lies inside the tightest bracketing interval)
371 System.arraycopy(x, signChangeIndex, x, signChangeIndex + 1, nbPoints - signChangeIndex);
372 x[signChangeIndex] = nextX;
373 System.arraycopy(y, signChangeIndex, y, signChangeIndex + 1, nbPoints - signChangeIndex);
374 y[signChangeIndex] = nextY;
375 ++nbPoints;
376
377 // update the bracketing interval
378 if (nextY.multiply(yA).negativeOrNull()) {
379 // the sign change occurs before the inserted point
380 xB = nextX;
381 yB = nextY;
382 absYB = yB.abs();
383 ++agingA;
384 agingB = 0;
385 } else {
386 // the sign change occurs after the inserted point
387 xA = nextX;
388 yA = nextY;
389 absYA = yA.abs();
390 agingA = 0;
391 ++agingB;
392
393 // update the sign change index
394 signChangeIndex++;
395
396 }
397
398 }
399
400 }
401
402 /** Guess an x value by n<sup>th</sup> order inverse polynomial interpolation.
403 * <p>
404 * The x value is guessed by evaluating polynomial Q(y) at y = targetY, where Q
405 * is built such that for all considered points (x<sub>i</sub>, y<sub>i</sub>),
406 * Q(y<sub>i</sub>) = x<sub>i</sub>.
407 * </p>
408 * @param targetY target value for y
409 * @param x reference points abscissas for interpolation,
410 * note that this array <em>is</em> modified during computation
411 * @param y reference points ordinates for interpolation
412 * @param start start index of the points to consider (inclusive)
413 * @param end end index of the points to consider (exclusive)
414 * @return guessed root (will be a NaN if two points share the same y)
415 */
416 private Dfp guessX(final Dfp targetY, final Dfp[] x, final Dfp[] y,
417 final int start, final int end) {
418
419 // compute Q Newton coefficients by divided differences
420 for (int i = start; i < end - 1; ++i) {
421 final int delta = i + 1 - start;
422 for (int j = end - 1; j > i; --j) {
423 x[j] = x[j].subtract(x[j-1]).divide(y[j].subtract(y[j - delta]));
424 }
425 }
426
427 // evaluate Q(targetY)
428 Dfp x0 = targetY.getZero();
429 for (int j = end - 1; j >= start; --j) {
430 x0 = x[j].add(x0.multiply(targetY.subtract(y[j])));
431 }
432
433 return x0;
434
435 }
436
437 }