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