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