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 */
017package org.apache.commons.math3.analysis.solvers;
018
019
020import org.apache.commons.math3.analysis.UnivariateFunction;
021import org.apache.commons.math3.exception.MathInternalError;
022import org.apache.commons.math3.exception.NoBracketingException;
023import org.apache.commons.math3.exception.NumberIsTooSmallException;
024import org.apache.commons.math3.exception.NumberIsTooLargeException;
025import org.apache.commons.math3.exception.TooManyEvaluationsException;
026import org.apache.commons.math3.util.FastMath;
027import 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 */
044public class BracketingNthOrderBrentSolver
045    extends AbstractUnivariateSolver
046    implements BracketedUnivariateSolver<UnivariateFunction> {
047
048    /** Default absolute accuracy. */
049    private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
050
051    /** Default maximal order. */
052    private static final int DEFAULT_MAXIMAL_ORDER = 5;
053
054    /** Maximal aging triggering an attempt to balance the bracketing interval. */
055    private static final int MAXIMAL_AGING = 2;
056
057    /** Reduction factor for attempts to balance the bracketing interval. */
058    private static final double REDUCTION_FACTOR = 1.0 / 16.0;
059
060    /** Maximal order. */
061    private final int maximalOrder;
062
063    /** The kinds of solutions that the algorithm may accept. */
064    private AllowedSolution allowed;
065
066    /**
067     * Construct a solver with default accuracy and maximal order (1e-6 and 5 respectively)
068     */
069    public BracketingNthOrderBrentSolver() {
070        this(DEFAULT_ABSOLUTE_ACCURACY, DEFAULT_MAXIMAL_ORDER);
071    }
072
073    /**
074     * Construct a solver.
075     *
076     * @param absoluteAccuracy Absolute accuracy.
077     * @param maximalOrder maximal order.
078     * @exception NumberIsTooSmallException if maximal order is lower than 2
079     */
080    public BracketingNthOrderBrentSolver(final double absoluteAccuracy,
081                                         final int maximalOrder)
082        throws NumberIsTooSmallException {
083        super(absoluteAccuracy);
084        if (maximalOrder < 2) {
085            throw new NumberIsTooSmallException(maximalOrder, 2, true);
086        }
087        this.maximalOrder = maximalOrder;
088        this.allowed = AllowedSolution.ANY_SIDE;
089    }
090
091    /**
092     * Construct a solver.
093     *
094     * @param relativeAccuracy Relative accuracy.
095     * @param absoluteAccuracy Absolute accuracy.
096     * @param maximalOrder maximal order.
097     * @exception NumberIsTooSmallException if maximal order is lower than 2
098     */
099    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}