    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    }