001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    package org.apache.commons.math3.analysis.solvers;
018    
019    
020    import org.apache.commons.math3.analysis.UnivariateFunction;
021    import org.apache.commons.math3.exception.MathInternalError;
022    import org.apache.commons.math3.exception.NoBracketingException;
023    import org.apache.commons.math3.exception.NumberIsTooSmallException;
024    import org.apache.commons.math3.exception.NumberIsTooLargeException;
025    import org.apache.commons.math3.exception.TooManyEvaluationsException;
026    import org.apache.commons.math3.util.FastMath;
027    import org.apache.commons.math3.util.Precision;
028    
029    /**
030     * This class implements a modification of the <a
031     * href="http://mathworld.wolfram.com/BrentsMethod.html"> Brent algorithm</a>.
032     * <p>
033     * The changes with respect to the original Brent algorithm are:
034     * <ul>
035     *   <li>the returned value is chosen in the current interval according
036     *   to user specified {@link AllowedSolution},</li>
037     *   <li>the maximal order for the invert polynomial root search is
038     *   user-specified instead of being invert quadratic only</li>
039     * </ul>
040     * </p>
041     * The given interval must bracket the root.
042     *
043     * @version $Id: BracketingNthOrderBrentSolver.java 1379560 2012-08-31 19:40:30Z erans $
044     */
045    public class BracketingNthOrderBrentSolver
046        extends AbstractUnivariateSolver
047        implements BracketedUnivariateSolver<UnivariateFunction> {
048    
049        /** Default absolute accuracy. */
050        private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
051    
052        /** Default maximal order. */
053        private static final int DEFAULT_MAXIMAL_ORDER = 5;
054    
055        /** Maximal aging triggering an attempt to balance the bracketing interval. */
056        private static final int MAXIMAL_AGING = 2;
057    
058        /** Reduction factor for attempts to balance the bracketing interval. */
059        private static final double REDUCTION_FACTOR = 1.0 / 16.0;
060    
061        /** Maximal order. */
062        private final int maximalOrder;
063    
064        /** The kinds of solutions that the algorithm may accept. */
065        private AllowedSolution allowed;
066    
067        /**
068         * Construct a solver with default accuracy and maximal order (1e-6 and 5 respectively)
069         */
070        public BracketingNthOrderBrentSolver() {
071            this(DEFAULT_ABSOLUTE_ACCURACY, DEFAULT_MAXIMAL_ORDER);
072        }
073    
074        /**
075         * Construct a solver.
076         *
077         * @param absoluteAccuracy Absolute accuracy.
078         * @param maximalOrder maximal order.
079         * @exception NumberIsTooSmallException if maximal order is lower than 2
080         */
081        public BracketingNthOrderBrentSolver(final double absoluteAccuracy,
082                                             final int maximalOrder)
083            throws NumberIsTooSmallException {
084            super(absoluteAccuracy);
085            if (maximalOrder < 2) {
086                throw new NumberIsTooSmallException(maximalOrder, 2, true);
087            }
088            this.maximalOrder = maximalOrder;
089            this.allowed = AllowedSolution.ANY_SIDE;
090        }
091    
092        /**
093         * Construct a solver.
094         *
095         * @param relativeAccuracy Relative accuracy.
096         * @param absoluteAccuracy Absolute accuracy.
097         * @param maximalOrder maximal order.
098         * @exception NumberIsTooSmallException if maximal order is lower than 2
099         */
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    }