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   */
43  public class BracketingNthOrderBrentSolverDFP {
44  
45     /** Maximal aging triggering an attempt to balance the bracketing interval. */
46      private static final int MAXIMAL_AGING = 2;
47  
48      /** Maximal order. */
49      private final int maximalOrder;
50  
51      /** Function value accuracy. */
52      private final Dfp functionValueAccuracy;
53  
54      /** Absolute accuracy. */
55      private final Dfp absoluteAccuracy;
56  
57      /** Relative accuracy. */
58      private final Dfp relativeAccuracy;
59  
60      /** Evaluations counter. */
61      private final Incrementor evaluations = new Incrementor();
62  
63      /**
64       * Construct a solver.
65       *
66       * @param relativeAccuracy Relative accuracy.
67       * @param absoluteAccuracy Absolute accuracy.
68       * @param functionValueAccuracy Function value accuracy.
69       * @param maximalOrder maximal order.
70       * @exception NumberIsTooSmallException if maximal order is lower than 2
71       */
72      public BracketingNthOrderBrentSolverDFP(final Dfp relativeAccuracy,
73                                              final Dfp absoluteAccuracy,
74                                              final Dfp functionValueAccuracy,
75                                              final int maximalOrder)
76          throws NumberIsTooSmallException {
77          if (maximalOrder < 2) {
78              throw new NumberIsTooSmallException(maximalOrder, 2, true);
79          }
80          this.maximalOrder = maximalOrder;
81          this.absoluteAccuracy = absoluteAccuracy;
82          this.relativeAccuracy = relativeAccuracy;
83          this.functionValueAccuracy = functionValueAccuracy;
84      }
85  
86      /** Get the maximal order.
87       * @return maximal order
88       */
89      public int getMaximalOrder() {
90          return maximalOrder;
91      }
92  
93      /**
94       * Get the maximal number of function evaluations.
95       *
96       * @return the maximal number of function evaluations.
97       */
98      public int getMaxEvaluations() {
99          return evaluations.getMaximalCount();
100     }
101 
102     /**
103      * Get the number of evaluations of the objective function.
104      * The number of evaluations corresponds to the last call to the
105      * {@code optimize} method. It is 0 if the method has not been
106      * called yet.
107      *
108      * @return the number of evaluations of the objective function.
109      */
110     public int getEvaluations() {
111         return evaluations.getCount();
112     }
113 
114     /**
115      * Get the absolute accuracy.
116      * @return absolute accuracy
117      */
118     public Dfp getAbsoluteAccuracy() {
119         return absoluteAccuracy;
120     }
121 
122     /**
123      * Get the relative accuracy.
124      * @return relative accuracy
125      */
126     public Dfp getRelativeAccuracy() {
127         return relativeAccuracy;
128     }
129 
130     /**
131      * Get the function accuracy.
132      * @return function accuracy
133      */
134     public Dfp getFunctionValueAccuracy() {
135         return functionValueAccuracy;
136     }
137 
138     /**
139      * Solve for a zero in the given interval.
140      * A solver may require that the interval brackets a single zero root.
141      * Solvers that do require bracketing should be able to handle the case
142      * where one of the endpoints is itself a root.
143      *
144      * @param maxEval Maximum number of evaluations.
145      * @param f Function to solve.
146      * @param min Lower bound for the interval.
147      * @param max Upper bound for the interval.
148      * @param allowedSolution The kind of solutions that the root-finding algorithm may
149      * accept as solutions.
150      * @return a value where the function is zero.
151      * @exception NullArgumentException if f is null.
152      * @exception NoBracketingException if root cannot be bracketed
153      */
154     public Dfp solve(final int maxEval, final UnivariateDfpFunction f,
155                      final Dfp min, final Dfp max, final AllowedSolution allowedSolution)
156         throws NullArgumentException, NoBracketingException {
157         return solve(maxEval, f, min, max, min.add(max).divide(2), allowedSolution);
158     }
159 
160     /**
161      * Solve for a zero in the given interval, start at {@code startValue}.
162      * A solver may require that the interval brackets a single zero root.
163      * Solvers that do require bracketing should be able to handle the case
164      * where one of the endpoints is itself a root.
165      *
166      * @param maxEval Maximum number of evaluations.
167      * @param f Function to solve.
168      * @param min Lower bound for the interval.
169      * @param max Upper bound for the interval.
170      * @param startValue Start value to use.
171      * @param allowedSolution The kind of solutions that the root-finding algorithm may
172      * accept as solutions.
173      * @return a value where the function is zero.
174      * @exception NullArgumentException if f is null.
175      * @exception NoBracketingException if root cannot be bracketed
176      */
177     public Dfp solve(final int maxEval, final UnivariateDfpFunction f,
178                      final Dfp min, final Dfp max, final Dfp startValue,
179                      final AllowedSolution allowedSolution)
180         throws NullArgumentException, NoBracketingException {
181 
182         // Checks.
183         MathUtils.checkNotNull(f);
184 
185         // Reset.
186         evaluations.setMaximalCount(maxEval);
187         evaluations.resetCount();
188         Dfp zero = startValue.getZero();
189         Dfp nan  = zero.newInstance((byte) 1, Dfp.QNAN);
190 
191         // prepare arrays with the first points
192         final Dfp[] x = new Dfp[maximalOrder + 1];
193         final Dfp[] y = new Dfp[maximalOrder + 1];
194         x[0] = min;
195         x[1] = startValue;
196         x[2] = max;
197 
198         // evaluate initial guess
199         evaluations.incrementCount();
200         y[1] = f.value(x[1]);
201         if (y[1].isZero()) {
202             // return the initial guess if it is a perfect root.
203             return x[1];
204         }
205 
206         // evaluate first  endpoint
207         evaluations.incrementCount();
208         y[0] = f.value(x[0]);
209         if (y[0].isZero()) {
210             // return the first endpoint if it is a perfect root.
211             return x[0];
212         }
213 
214         int nbPoints;
215         int signChangeIndex;
216         if (y[0].multiply(y[1]).negativeOrNull()) {
217 
218             // reduce interval if it brackets the root
219             nbPoints        = 2;
220             signChangeIndex = 1;
221 
222         } else {
223 
224             // evaluate second endpoint
225             evaluations.incrementCount();
226             y[2] = f.value(x[2]);
227             if (y[2].isZero()) {
228                 // return the second endpoint if it is a perfect root.
229                 return x[2];
230             }
231 
232             if (y[1].multiply(y[2]).negativeOrNull()) {
233                 // use all computed point as a start sampling array for solving
234                 nbPoints        = 3;
235                 signChangeIndex = 2;
236             } else {
237                 throw new NoBracketingException(x[0].toDouble(), x[2].toDouble(),
238                                                 y[0].toDouble(), y[2].toDouble());
239             }
240 
241         }
242 
243         // prepare a work array for inverse polynomial interpolation
244         final Dfp[] tmpX = new Dfp[x.length];
245 
246         // current tightest bracketing of the root
247         Dfp xA    = x[signChangeIndex - 1];
248         Dfp yA    = y[signChangeIndex - 1];
249         Dfp absXA = xA.abs();
250         Dfp absYA = yA.abs();
251         int agingA   = 0;
252         Dfp xB    = x[signChangeIndex];
253         Dfp yB    = y[signChangeIndex];
254         Dfp absXB = xB.abs();
255         Dfp absYB = yB.abs();
256         int agingB   = 0;
257 
258         // search loop
259         while (true) {
260 
261             // check convergence of bracketing interval
262             Dfp maxX = absXA.lessThan(absXB) ? absXB : absXA;
263             Dfp maxY = absYA.lessThan(absYB) ? absYB : absYA;
264             final Dfp xTol = absoluteAccuracy.add(relativeAccuracy.multiply(maxX));
265             if (xB.subtract(xA).subtract(xTol).negativeOrNull() ||
266                 maxY.lessThan(functionValueAccuracy)) {
267                 switch (allowedSolution) {
268                 case ANY_SIDE :
269                     return absYA.lessThan(absYB) ? xA : xB;
270                 case LEFT_SIDE :
271                     return xA;
272                 case RIGHT_SIDE :
273                     return xB;
274                 case BELOW_SIDE :
275                     return yA.lessThan(zero) ? xA : xB;
276                 case ABOVE_SIDE :
277                     return yA.lessThan(zero) ? xB : xA;
278                 default :
279                     // this should never happen
280                     throw new MathInternalError(null);
281                 }
282             }
283 
284             // target for the next evaluation point
285             Dfp targetY;
286             if (agingA >= MAXIMAL_AGING) {
287                 // we keep updating the high bracket, try to compensate this
288                 targetY = yB.divide(16).negate();
289             } else if (agingB >= MAXIMAL_AGING) {
290                 // we keep updating the low bracket, try to compensate this
291                 targetY = yA.divide(16).negate();
292             } else {
293                 // bracketing is balanced, try to find the root itself
294                 targetY = zero;
295             }
296 
297             // make a few attempts to guess a root,
298             Dfp nextX;
299             int start = 0;
300             int end   = nbPoints;
301             do {
302 
303                 // guess a value for current target, using inverse polynomial interpolation
304                 System.arraycopy(x, start, tmpX, start, end - start);
305                 nextX = guessX(targetY, tmpX, y, start, end);
306 
307                 if (!(nextX.greaterThan(xA) && nextX.lessThan(xB))) {
308                     // the guessed root is not strictly inside of the tightest bracketing interval
309 
310                     // the guessed root is either not strictly inside the interval or it
311                     // is a NaN (which occurs when some sampling points share the same y)
312                     // we try again with a lower interpolation order
313                     if (signChangeIndex - start >= end - signChangeIndex) {
314                         // we have more points before the sign change, drop the lowest point
315                         ++start;
316                     } else {
317                         // we have more points after sign change, drop the highest point
318                         --end;
319                     }
320 
321                     // we need to do one more attempt
322                     nextX = nan;
323 
324                 }
325 
326             } while (nextX.isNaN() && (end - start > 1));
327 
328             if (nextX.isNaN()) {
329                 // fall back to bisection
330                 nextX = xA.add(xB.subtract(xA).divide(2));
331                 start = signChangeIndex - 1;
332                 end   = signChangeIndex;
333             }
334 
335             // evaluate the function at the guessed root
336             evaluations.incrementCount();
337             final Dfp nextY = f.value(nextX);
338             if (nextY.isZero()) {
339                 // we have found an exact root, since it is not an approximation
340                 // we don't need to bother about the allowed solutions setting
341                 return nextX;
342             }
343 
344             if ((nbPoints > 2) && (end - start != nbPoints)) {
345 
346                 // we have been forced to ignore some points to keep bracketing,
347                 // they are probably too far from the root, drop them from now on
348                 nbPoints = end - start;
349                 System.arraycopy(x, start, x, 0, nbPoints);
350                 System.arraycopy(y, start, y, 0, nbPoints);
351                 signChangeIndex -= start;
352 
353             } else  if (nbPoints == x.length) {
354 
355                 // we have to drop one point in order to insert the new one
356                 nbPoints--;
357 
358                 // keep the tightest bracketing interval as centered as possible
359                 if (signChangeIndex >= (x.length + 1) / 2) {
360                     // we drop the lowest point, we have to shift the arrays and the index
361                     System.arraycopy(x, 1, x, 0, nbPoints);
362                     System.arraycopy(y, 1, y, 0, nbPoints);
363                     --signChangeIndex;
364                 }
365 
366             }
367 
368             // insert the last computed point
369             //(by construction, we know it lies inside the tightest bracketing interval)
370             System.arraycopy(x, signChangeIndex, x, signChangeIndex + 1, nbPoints - signChangeIndex);
371             x[signChangeIndex] = nextX;
372             System.arraycopy(y, signChangeIndex, y, signChangeIndex + 1, nbPoints - signChangeIndex);
373             y[signChangeIndex] = nextY;
374             ++nbPoints;
375 
376             // update the bracketing interval
377             if (nextY.multiply(yA).negativeOrNull()) {
378                 // the sign change occurs before the inserted point
379                 xB = nextX;
380                 yB = nextY;
381                 absYB = yB.abs();
382                 ++agingA;
383                 agingB = 0;
384             } else {
385                 // the sign change occurs after the inserted point
386                 xA = nextX;
387                 yA = nextY;
388                 absYA = yA.abs();
389                 agingA = 0;
390                 ++agingB;
391 
392                 // update the sign change index
393                 signChangeIndex++;
394 
395             }
396 
397         }
398 
399     }
400 
401     /** Guess an x value by n<sup>th</sup> order inverse polynomial interpolation.
402      * <p>
403      * The x value is guessed by evaluating polynomial Q(y) at y = targetY, where Q
404      * is built such that for all considered points (x<sub>i</sub>, y<sub>i</sub>),
405      * Q(y<sub>i</sub>) = x<sub>i</sub>.
406      * </p>
407      * @param targetY target value for y
408      * @param x reference points abscissas for interpolation,
409      * note that this array <em>is</em> modified during computation
410      * @param y reference points ordinates for interpolation
411      * @param start start index of the points to consider (inclusive)
412      * @param end end index of the points to consider (exclusive)
413      * @return guessed root (will be a NaN if two points share the same y)
414      */
415     private Dfp guessX(final Dfp targetY, final Dfp[] x, final Dfp[] y,
416                        final int start, final int end) {
417 
418         // compute Q Newton coefficients by divided differences
419         for (int i = start; i < end - 1; ++i) {
420             final int delta = i + 1 - start;
421             for (int j = end - 1; j > i; --j) {
422                 x[j] = x[j].subtract(x[j-1]).divide(y[j].subtract(y[j - delta]));
423             }
424         }
425 
426         // evaluate Q(targetY)
427         Dfp x0 = targetY.getZero();
428         for (int j = end - 1; j >= start; --j) {
429             x0 = x[j].add(x0.multiply(targetY.subtract(y[j])));
430         }
431 
432         return x0;
433 
434     }
435 
436 }