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.math4.analysis.solvers;
018
019
020import org.apache.commons.math4.analysis.UnivariateFunction;
021import org.apache.commons.math4.exception.MathInternalError;
022import org.apache.commons.math4.exception.NoBracketingException;
023import org.apache.commons.math4.exception.NumberIsTooLargeException;
024import org.apache.commons.math4.exception.NumberIsTooSmallException;
025import org.apache.commons.math4.exception.TooManyEvaluationsException;
026import org.apache.commons.math4.util.FastMath;
027import org.apache.commons.numbers.core.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><p>
040 * The given interval must bracket the root.</p>
041 *
042 */
043public class BracketingNthOrderBrentSolver
044    extends AbstractUnivariateSolver
045    implements BracketedUnivariateSolver<UnivariateFunction> {
046
047    /** Default absolute accuracy. */
048    private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
049
050    /** Default maximal order. */
051    private static final int DEFAULT_MAXIMAL_ORDER = 5;
052
053    /** Maximal aging triggering an attempt to balance the bracketing interval. */
054    private static final int MAXIMAL_AGING = 2;
055
056    /** Reduction factor for attempts to balance the bracketing interval. */
057    private static final double REDUCTION_FACTOR = 1.0 / 16.0;
058
059    /** Maximal order. */
060    private final int maximalOrder;
061
062    /** The kinds of solutions that the algorithm may accept. */
063    private AllowedSolution allowed;
064
065    /**
066     * Construct a solver with default accuracy and maximal order (1e-6 and 5 respectively)
067     */
068    public BracketingNthOrderBrentSolver() {
069        this(DEFAULT_ABSOLUTE_ACCURACY, DEFAULT_MAXIMAL_ORDER);
070    }
071
072    /**
073     * Construct a solver.
074     *
075     * @param absoluteAccuracy Absolute accuracy.
076     * @param maximalOrder maximal order.
077     * @exception NumberIsTooSmallException if maximal order is lower than 2
078     */
079    public BracketingNthOrderBrentSolver(final double absoluteAccuracy,
080                                         final int maximalOrder)
081        throws NumberIsTooSmallException {
082        super(absoluteAccuracy);
083        if (maximalOrder < 2) {
084            throw new NumberIsTooSmallException(maximalOrder, 2, true);
085        }
086        this.maximalOrder = maximalOrder;
087        this.allowed = AllowedSolution.ANY_SIDE;
088    }
089
090    /**
091     * Construct a solver.
092     *
093     * @param relativeAccuracy Relative accuracy.
094     * @param absoluteAccuracy Absolute accuracy.
095     * @param maximalOrder maximal order.
096     * @exception NumberIsTooSmallException if maximal order is lower than 2
097     */
098    public BracketingNthOrderBrentSolver(final double relativeAccuracy,
099                                         final double absoluteAccuracy,
100                                         final int maximalOrder)
101        throws NumberIsTooSmallException {
102        super(relativeAccuracy, absoluteAccuracy);
103        if (maximalOrder < 2) {
104            throw new NumberIsTooSmallException(maximalOrder, 2, true);
105        }
106        this.maximalOrder = maximalOrder;
107        this.allowed = AllowedSolution.ANY_SIDE;
108    }
109
110    /**
111     * Construct a solver.
112     *
113     * @param relativeAccuracy Relative accuracy.
114     * @param absoluteAccuracy Absolute accuracy.
115     * @param functionValueAccuracy Function value accuracy.
116     * @param maximalOrder maximal order.
117     * @exception NumberIsTooSmallException if maximal order is lower than 2
118     */
119    public BracketingNthOrderBrentSolver(final double relativeAccuracy,
120                                         final double absoluteAccuracy,
121                                         final double functionValueAccuracy,
122                                         final int maximalOrder)
123        throws NumberIsTooSmallException {
124        super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
125        if (maximalOrder < 2) {
126            throw new NumberIsTooSmallException(maximalOrder, 2, true);
127        }
128        this.maximalOrder = maximalOrder;
129        this.allowed = AllowedSolution.ANY_SIDE;
130    }
131
132    /** Get the maximal order.
133     * @return maximal order
134     */
135    public int getMaximalOrder() {
136        return maximalOrder;
137    }
138
139    /**
140     * {@inheritDoc}
141     */
142    @Override
143    protected double doSolve()
144        throws TooManyEvaluationsException,
145               NumberIsTooLargeException,
146               NoBracketingException {
147        // prepare arrays with the first points
148        final double[] x = new double[maximalOrder + 1];
149        final double[] y = new double[maximalOrder + 1];
150        x[0] = getMin();
151        x[1] = getStartValue();
152        x[2] = getMax();
153        verifySequence(x[0], x[1], x[2]);
154
155        // evaluate initial guess
156        y[1] = computeObjectiveValue(x[1]);
157        if (Precision.equals(y[1], 0.0, 1)) {
158            // return the initial guess if it is a perfect root.
159            return x[1];
160        }
161
162        // evaluate first  endpoint
163        y[0] = computeObjectiveValue(x[0]);
164        if (Precision.equals(y[0], 0.0, 1)) {
165            // return the first endpoint if it is a perfect root.
166            return x[0];
167        }
168
169        int nbPoints;
170        int signChangeIndex;
171        if (y[0] * y[1] < 0) {
172
173            // reduce interval if it brackets the root
174            nbPoints        = 2;
175            signChangeIndex = 1;
176
177        } else {
178
179            // evaluate second endpoint
180            y[2] = computeObjectiveValue(x[2]);
181            if (Precision.equals(y[2], 0.0, 1)) {
182                // return the second endpoint if it is a perfect root.
183                return x[2];
184            }
185
186            if (y[1] * y[2] < 0) {
187                // use all computed point as a start sampling array for solving
188                nbPoints        = 3;
189                signChangeIndex = 2;
190            } else {
191                throw new NoBracketingException(x[0], x[2], y[0], y[2]);
192            }
193
194        }
195
196        // prepare a work array for inverse polynomial interpolation
197        final double[] tmpX = new double[x.length];
198
199        // current tightest bracketing of the root
200        double xA    = x[signChangeIndex - 1];
201        double yA    = y[signChangeIndex - 1];
202        double absYA = FastMath.abs(yA);
203        int agingA   = 0;
204        double xB    = x[signChangeIndex];
205        double yB    = y[signChangeIndex];
206        double absYB = FastMath.abs(yB);
207        int agingB   = 0;
208
209        // search loop
210        while (true) {
211
212            // check convergence of bracketing interval
213            final double xTol = getAbsoluteAccuracy() +
214                                getRelativeAccuracy() * FastMath.max(FastMath.abs(xA), FastMath.abs(xB));
215            if (((xB - xA) <= xTol) || (FastMath.max(absYA, absYB) < getFunctionValueAccuracy())) {
216                switch (allowed) {
217                case ANY_SIDE :
218                    return absYA < absYB ? xA : xB;
219                case LEFT_SIDE :
220                    return xA;
221                case RIGHT_SIDE :
222                    return xB;
223                case BELOW_SIDE :
224                    return (yA <= 0) ? xA : xB;
225                case ABOVE_SIDE :
226                    return (yA <  0) ? xB : xA;
227                default :
228                    // this should never happen
229                    throw new MathInternalError();
230                }
231            }
232
233            // target for the next evaluation point
234            double targetY;
235            if (agingA >= MAXIMAL_AGING) {
236                // we keep updating the high bracket, try to compensate this
237                final int p = agingA - MAXIMAL_AGING;
238                final double weightA = (1 << p) - 1;
239                final double weightB = p + 1;
240                targetY = (weightA * yA - weightB * REDUCTION_FACTOR * yB) / (weightA + weightB);
241            } else if (agingB >= MAXIMAL_AGING) {
242                // we keep updating the low bracket, try to compensate this
243                final int p = agingB - MAXIMAL_AGING;
244                final double weightA = p + 1;
245                final double weightB = (1 << p) - 1;
246                targetY = (weightB * yB - weightA * REDUCTION_FACTOR * yA) / (weightA + weightB);
247            } else {
248                // bracketing is balanced, try to find the root itself
249                targetY = 0;
250            }
251
252            // make a few attempts to guess a root,
253            double nextX;
254            int start = 0;
255            int end   = nbPoints;
256            do {
257
258                // guess a value for current target, using inverse polynomial interpolation
259                System.arraycopy(x, start, tmpX, start, end - start);
260                nextX = guessX(targetY, tmpX, y, start, end);
261
262                if (!((nextX > xA) && (nextX < xB))) {
263                    // the guessed root is not strictly inside of the tightest bracketing interval
264
265                    // the guessed root is either not strictly inside the interval or it
266                    // is a NaN (which occurs when some sampling points share the same y)
267                    // we try again with a lower interpolation order
268                    if (signChangeIndex - start >= end - signChangeIndex) {
269                        // we have more points before the sign change, drop the lowest point
270                        ++start;
271                    } else {
272                        // we have more points after sign change, drop the highest point
273                        --end;
274                    }
275
276                    // we need to do one more attempt
277                    nextX = Double.NaN;
278
279                }
280
281            } while (Double.isNaN(nextX) && (end - start > 1));
282
283            if (Double.isNaN(nextX)) {
284                // fall back to bisection
285                nextX = xA + 0.5 * (xB - xA);
286                start = signChangeIndex - 1;
287                end   = signChangeIndex;
288            }
289
290            // evaluate the function at the guessed root
291            final double nextY = computeObjectiveValue(nextX);
292            if (Precision.equals(nextY, 0.0, 1)) {
293                // we have found an exact root, since it is not an approximation
294                // we don't need to bother about the allowed solutions setting
295                return nextX;
296            }
297
298            if ((nbPoints > 2) && (end - start != nbPoints)) {
299
300                // we have been forced to ignore some points to keep bracketing,
301                // they are probably too far from the root, drop them from now on
302                nbPoints = end - start;
303                System.arraycopy(x, start, x, 0, nbPoints);
304                System.arraycopy(y, start, y, 0, nbPoints);
305                signChangeIndex -= start;
306
307            } else  if (nbPoints == x.length) {
308
309                // we have to drop one point in order to insert the new one
310                nbPoints--;
311
312                // keep the tightest bracketing interval as centered as possible
313                if (signChangeIndex >= (x.length + 1) / 2) {
314                    // we drop the lowest point, we have to shift the arrays and the index
315                    System.arraycopy(x, 1, x, 0, nbPoints);
316                    System.arraycopy(y, 1, y, 0, nbPoints);
317                    --signChangeIndex;
318                }
319
320            }
321
322            // insert the last computed point
323            //(by construction, we know it lies inside the tightest bracketing interval)
324            System.arraycopy(x, signChangeIndex, x, signChangeIndex + 1, nbPoints - signChangeIndex);
325            x[signChangeIndex] = nextX;
326            System.arraycopy(y, signChangeIndex, y, signChangeIndex + 1, nbPoints - signChangeIndex);
327            y[signChangeIndex] = nextY;
328            ++nbPoints;
329
330            // update the bracketing interval
331            if (nextY * yA <= 0) {
332                // the sign change occurs before the inserted point
333                xB = nextX;
334                yB = nextY;
335                absYB = FastMath.abs(yB);
336                ++agingA;
337                agingB = 0;
338            } else {
339                // the sign change occurs after the inserted point
340                xA = nextX;
341                yA = nextY;
342                absYA = FastMath.abs(yA);
343                agingA = 0;
344                ++agingB;
345
346                // update the sign change index
347                signChangeIndex++;
348
349            }
350
351        }
352
353    }
354
355    /** Guess an x value by n<sup>th</sup> order inverse polynomial interpolation.
356     * <p>
357     * The x value is guessed by evaluating polynomial Q(y) at y = targetY, where Q
358     * is built such that for all considered points (x<sub>i</sub>, y<sub>i</sub>),
359     * Q(y<sub>i</sub>) = x<sub>i</sub>.
360     * </p>
361     * @param targetY target value for y
362     * @param x reference points abscissas for interpolation,
363     * note that this array <em>is</em> modified during computation
364     * @param y reference points ordinates for interpolation
365     * @param start start index of the points to consider (inclusive)
366     * @param end end index of the points to consider (exclusive)
367     * @return guessed root (will be a NaN if two points share the same y)
368     */
369    private double guessX(final double targetY, final double[] x, final double[] y,
370                          final int start, final int end) {
371
372        // compute Q Newton coefficients by divided differences
373        for (int i = start; i < end - 1; ++i) {
374            final int delta = i + 1 - start;
375            for (int j = end - 1; j > i; --j) {
376                x[j] = (x[j] - x[j-1]) / (y[j] - y[j - delta]);
377            }
378        }
379
380        // evaluate Q(targetY)
381        double x0 = 0;
382        for (int j = end - 1; j >= start; --j) {
383            x0 = x[j] + x0 * (targetY - y[j]);
384        }
385
386        return x0;
387
388    }
389
390    /** {@inheritDoc} */
391    @Override
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    @Override
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}