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