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 * @version $Id: BracketingNthOrderBrentSolver.java 1379560 2012-08-31 19:40:30Z erans $
044 */
045public 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}