View Javadoc
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    *      http://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  package org.apache.commons.math3.dfp;
18  
19  
20  import org.apache.commons.math3.analysis.solvers.AllowedSolution;
21  import org.apache.commons.math3.exception.MathInternalError;
22  import org.apache.commons.math3.exception.NoBracketingException;
23  import org.apache.commons.math3.exception.NullArgumentException;
24  import org.apache.commons.math3.exception.NumberIsTooSmallException;
25  import org.apache.commons.math3.util.Incrementor;
26  import org.apache.commons.math3.util.MathUtils;
27  
28  /**
29   * This class implements a modification of the <a
30   * href="http://mathworld.wolfram.com/BrentsMethod.html"> Brent algorithm</a>.
31   * <p>
32   * The changes with respect to the original Brent algorithm are:
33   * <ul>
34   *   <li>the returned value is chosen in the current interval according
35   *   to user specified {@link AllowedSolution},</li>
36   *   <li>the maximal order for the invert polynomial root search is
37   *   user-specified instead of being invert quadratic only</li>
38   * </ul>
39   * </p>
40   * The given interval must bracket the root.
41   *
42   * @version $Id: BracketingNthOrderBrentSolverDFP.java 1416643 2012-12-03 19:37:14Z tn $
43   */
44  public class BracketingNthOrderBrentSolverDFP {
45  
46     /** Maximal aging triggering an attempt to balance the bracketing interval. */
47      private static final int MAXIMAL_AGING = 2;
48  
49      /** Maximal order. */
50      private final int maximalOrder;
51  
52      /** Function value accuracy. */
53      private final Dfp functionValueAccuracy;
54  
55      /** Absolute accuracy. */
56      private final Dfp absoluteAccuracy;
57  
58      /** Relative accuracy. */
59      private final Dfp relativeAccuracy;
60  
61      /** Evaluations counter. */
62      private final Incrementor evaluations = new Incrementor();
63  
64      /**
65       * Construct a solver.
66       *
67       * @param relativeAccuracy Relative accuracy.
68       * @param absoluteAccuracy Absolute accuracy.
69       * @param functionValueAccuracy Function value accuracy.
70       * @param maximalOrder maximal order.
71       * @exception NumberIsTooSmallException if maximal order is lower than 2
72       */
73      public BracketingNthOrderBrentSolverDFP(final Dfp relativeAccuracy,
74                                              final Dfp absoluteAccuracy,
75                                              final Dfp functionValueAccuracy,
76                                              final int maximalOrder)
77          throws NumberIsTooSmallException {
78          if (maximalOrder < 2) {
79              throw new NumberIsTooSmallException(maximalOrder, 2, true);
80          }
81          this.maximalOrder = maximalOrder;
82          this.absoluteAccuracy = absoluteAccuracy;
83          this.relativeAccuracy = relativeAccuracy;
84          this.functionValueAccuracy = functionValueAccuracy;
85      }
86  
87      /** Get the maximal order.
88       * @return maximal order
89       */
90      public int getMaximalOrder() {
91          return maximalOrder;
92      }
93  
94      /**
95       * Get the maximal number of function evaluations.
96       *
97       * @return the maximal number of function evaluations.
98       */
99      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 }