    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.dfp;
    018
    019
    020    import org.apache.commons.math3.analysis.solvers.AllowedSolution;
    021    import org.apache.commons.math3.exception.MathInternalError;
    022    import org.apache.commons.math3.exception.NoBracketingException;
    023    import org.apache.commons.math3.exception.NullArgumentException;
    024    import org.apache.commons.math3.exception.NumberIsTooSmallException;
    025    import org.apache.commons.math3.util.Incrementor;
    026    import org.apache.commons.math3.util.MathUtils;
    027
    028    /**
    029     * This class implements a modification of the <a
    030     * href="http://mathworld.wolfram.com/BrentsMethod.html"> Brent algorithm</a>.
    031     * <p>
    032     * The changes with respect to the original Brent algorithm are:
    033     * <ul>
    034     *   <li>the returned value is chosen in the current interval according
    035     *   to user specified {@link AllowedSolution},</li>
    036     *   <li>the maximal order for the invert polynomial root search is
    037     *   user-specified instead of being invert quadratic only</li>
    038     * </ul>
    039     * </p>
    040     * The given interval must bracket the root.
    041     *
    042     * @version $Id: BracketingNthOrderBrentSolverDFP.java 1416643 2012-12-03 19:37:14Z tn $
    043     */
    044    public class BracketingNthOrderBrentSolverDFP {
    045
    046       /** Maximal aging triggering an attempt to balance the bracketing interval. */
    047        private static final int MAXIMAL_AGING = 2;
    048
    049        /** Maximal order. */
    050        private final int maximalOrder;
    051
    052        /** Function value accuracy. */
    053        private final Dfp functionValueAccuracy;
    054
    055        /** Absolute accuracy. */
    056        private final Dfp absoluteAccuracy;
    057
    058        /** Relative accuracy. */
    059        private final Dfp relativeAccuracy;
    060
    061        /** Evaluations counter. */
    062        private final Incrementor evaluations = new Incrementor();
    063
    064        /**
    065         * Construct a solver.
    066         *
    067         * @param relativeAccuracy Relative accuracy.
    068         * @param absoluteAccuracy Absolute accuracy.
    069         * @param functionValueAccuracy Function value accuracy.
    070         * @param maximalOrder maximal order.
    071         * @exception NumberIsTooSmallException if maximal order is lower than 2
    072         */
    073        public BracketingNthOrderBrentSolverDFP(final Dfp relativeAccuracy,
    074                                                final Dfp absoluteAccuracy,
    075                                                final Dfp functionValueAccuracy,
    076                                                final int maximalOrder)
    077            throws NumberIsTooSmallException {
    078            if (maximalOrder < 2) {
    079                throw new NumberIsTooSmallException(maximalOrder, 2, true);
    080            }
    081            this.maximalOrder = maximalOrder;
    082            this.absoluteAccuracy = absoluteAccuracy;
    083            this.relativeAccuracy = relativeAccuracy;
    084            this.functionValueAccuracy = functionValueAccuracy;
    085        }
    086
    087        /** Get the maximal order.
    088         * @return maximal order
    089         */
    090        public int getMaximalOrder() {
    091            return maximalOrder;
    092        }
    093
    094        /**
    095         * Get the maximal number of function evaluations.
    096         *
    097         * @return the maximal number of function evaluations.
    098         */
    099        public int getMaxEvaluations() {
    100            return evaluations.getMaximalCount();
    101        }
    102
    103        /**
    104         * Get the number of evaluations of the objective function.
    105         * The number of evaluations corresponds to the last call to the
    106         * {@code optimize} method. It is 0 if the method has not been
    107         * called yet.
    108         *
    109         * @return the number of evaluations of the objective function.
    110         */
    111        public int getEvaluations() {
    112            return evaluations.getCount();
    113        }
    114
    115        /**
    116         * Get the absolute accuracy.
    117         * @return absolute accuracy
    118         */
    119        public Dfp getAbsoluteAccuracy() {
    120            return absoluteAccuracy;
    121        }
    122
    123        /**
    124         * Get the relative accuracy.
    125         * @return relative accuracy
    126         */
    127        public Dfp getRelativeAccuracy() {
    128            return relativeAccuracy;
    129        }
    130
    131        /**
    132         * Get the function accuracy.
    133         * @return function accuracy
    134         */
    135        public Dfp getFunctionValueAccuracy() {
    136            return functionValueAccuracy;
    137        }
    138
    139        /**
    140         * Solve for a zero in the given interval.
    141         * A solver may require that the interval brackets a single zero root.
    142         * Solvers that do require bracketing should be able to handle the case
    143         * where one of the endpoints is itself a root.
    144         *
    145         * @param maxEval Maximum number of evaluations.
    146         * @param f Function to solve.
    147         * @param min Lower bound for the interval.
    148         * @param max Upper bound for the interval.
    149         * @param allowedSolution The kind of solutions that the root-finding algorithm may
    150         * accept as solutions.
    151         * @return a value where the function is zero.
    152         * @exception NullArgumentException if f is null.
    153         * @exception NoBracketingException if root cannot be bracketed
    154         */
    155        public Dfp solve(final int maxEval, final UnivariateDfpFunction f,
    156                         final Dfp min, final Dfp max, final AllowedSolution allowedSolution)
    157            throws NullArgumentException, NoBracketingException {
    158            return solve(maxEval, f, min, max, min.add(max).divide(2), allowedSolution);
    159        }
    160
    161        /**
    162         * Solve for a zero in the given interval, start at {@code startValue}.
    163         * A solver may require that the interval brackets a single zero root.
    164         * Solvers that do require bracketing should be able to handle the case
    165         * where one of the endpoints is itself a root.
    166         *
    167         * @param maxEval Maximum number of evaluations.
    168         * @param f Function to solve.
    169         * @param min Lower bound for the interval.
    170         * @param max Upper bound for the interval.
    171         * @param startValue Start value to use.
    172         * @param allowedSolution The kind of solutions that the root-finding algorithm may
    173         * accept as solutions.
    174         * @return a value where the function is zero.
    175         * @exception NullArgumentException if f is null.
    176         * @exception NoBracketingException if root cannot be bracketed
    177         */
    178        public Dfp solve(final int maxEval, final UnivariateDfpFunction f,
    179                         final Dfp min, final Dfp max, final Dfp startValue,
    180                         final AllowedSolution allowedSolution)
    181            throws NullArgumentException, NoBracketingException {
    182
    183            // Checks.
    184            MathUtils.checkNotNull(f);
    185
    186            // Reset.
    187            evaluations.setMaximalCount(maxEval);
    188            evaluations.resetCount();
    189            Dfp zero = startValue.getZero();
    190            Dfp nan  = zero.newInstance((byte) 1, Dfp.QNAN);
    191
    192            // prepare arrays with the first points
    193            final Dfp[] x = new Dfp[maximalOrder + 1];
    194            final Dfp[] y = new Dfp[maximalOrder + 1];
    195            x[0] = min;
    196            x[1] = startValue;
    197            x[2] = max;
    198
    199            // evaluate initial guess
    200            evaluations.incrementCount();
    201            y[1] = f.value(x[1]);
    202            if (y[1].isZero()) {
    203                // return the initial guess if it is a perfect root.
    204                return x[1];
    205            }
    206
    207            // evaluate first  endpoint
    208            evaluations.incrementCount();
    209            y[0] = f.value(x[0]);
    210            if (y[0].isZero()) {
    211                // return the first endpoint if it is a perfect root.
    212                return x[0];
    213            }
    214
    215            int nbPoints;
    216            int signChangeIndex;
    217            if (y[0].multiply(y[1]).negativeOrNull()) {
    218
    219                // reduce interval if it brackets the root
    220                nbPoints        = 2;
    221                signChangeIndex = 1;
    222
    223            } else {
    224
    225                // evaluate second endpoint
    226                evaluations.incrementCount();
    227                y[2] = f.value(x[2]);
    228                if (y[2].isZero()) {
    229                    // return the second endpoint if it is a perfect root.
    230                    return x[2];
    231                }
    232
    233                if (y[1].multiply(y[2]).negativeOrNull()) {
    234                    // use all computed point as a start sampling array for solving
    235                    nbPoints        = 3;
    236                    signChangeIndex = 2;
    237                } else {
    238                    throw new NoBracketingException(x[0].toDouble(), x[2].toDouble(),
    239                                                    y[0].toDouble(), y[2].toDouble());
    240                }
    241
    242            }
    243
    244            // prepare a work array for inverse polynomial interpolation
    245            final Dfp[] tmpX = new Dfp[x.length];
    246
    247            // current tightest bracketing of the root
    248            Dfp xA    = x[signChangeIndex - 1];
    249            Dfp yA    = y[signChangeIndex - 1];
    250            Dfp absXA = xA.abs();
    251            Dfp absYA = yA.abs();
    252            int agingA   = 0;
    253            Dfp xB    = x[signChangeIndex];
    254            Dfp yB    = y[signChangeIndex];
    255            Dfp absXB = xB.abs();
    256            Dfp absYB = yB.abs();
    257            int agingB   = 0;
    258
    259            // search loop
    260            while (true) {
    261
    262                // check convergence of bracketing interval
    263                Dfp maxX = absXA.lessThan(absXB) ? absXB : absXA;
    264                Dfp maxY = absYA.lessThan(absYB) ? absYB : absYA;
    265                final Dfp xTol = absoluteAccuracy.add(relativeAccuracy.multiply(maxX));
    266                if (xB.subtract(xA).subtract(xTol).negativeOrNull() ||
    267                    maxY.lessThan(functionValueAccuracy)) {
    268                    switch (allowedSolution) {
    269                    case ANY_SIDE :
    270                        return absYA.lessThan(absYB) ? xA : xB;
    271                    case LEFT_SIDE :
    272                        return xA;
    273                    case RIGHT_SIDE :
    274                        return xB;
    275                    case BELOW_SIDE :
    276                        return yA.lessThan(zero) ? xA : xB;
    277                    case ABOVE_SIDE :
    278                        return yA.lessThan(zero) ? xB : xA;
    279                    default :
    280                        // this should never happen
    281                        throw new MathInternalError(null);
    282                    }
    283                }
    284
    285                // target for the next evaluation point
    286                Dfp targetY;
    287                if (agingA >= MAXIMAL_AGING) {
    288                    // we keep updating the high bracket, try to compensate this
    289                    targetY = yB.divide(16).negate();
    290                } else if (agingB >= MAXIMAL_AGING) {
    291                    // we keep updating the low bracket, try to compensate this
    292                    targetY = yA.divide(16).negate();
    293                } else {
    294                    // bracketing is balanced, try to find the root itself
    295                    targetY = zero;
    296                }
    297
    298                // make a few attempts to guess a root,
    299                Dfp nextX;
    300                int start = 0;
    301                int end   = nbPoints;
    302                do {
    303
    304                    // guess a value for current target, using inverse polynomial interpolation
    305                    System.arraycopy(x, start, tmpX, start, end - start);
    306                    nextX = guessX(targetY, tmpX, y, start, end);
    307
    308                    if (!(nextX.greaterThan(xA) && nextX.lessThan(xB))) {
    309                        // the guessed root is not strictly inside of the tightest bracketing interval
    310
    311                        // the guessed root is either not strictly inside the interval or it
    312                        // is a NaN (which occurs when some sampling points share the same y)
    313                        // we try again with a lower interpolation order
    314                        if (signChangeIndex - start >= end - signChangeIndex) {
    315                            // we have more points before the sign change, drop the lowest point
    316                            ++start;
    317                        } else {
    318                            // we have more points after sign change, drop the highest point
    319                            --end;
    320                        }
    321
    322                        // we need to do one more attempt
    323                        nextX = nan;
    324
    325                    }
    326
    327                } while (nextX.isNaN() && (end - start > 1));
    328
    329                if (nextX.isNaN()) {
    330                    // fall back to bisection
    331                    nextX = xA.add(xB.subtract(xA).divide(2));
    332                    start = signChangeIndex - 1;
    333                    end   = signChangeIndex;
    334                }
    335
    336                // evaluate the function at the guessed root
    337                evaluations.incrementCount();
    338                final Dfp nextY = f.value(nextX);
    339                if (nextY.isZero()) {
    340                    // we have found an exact root, since it is not an approximation
    341                    // we don't need to bother about the allowed solutions setting
    342                    return nextX;
    343                }
    344
    345                if ((nbPoints > 2) && (end - start != nbPoints)) {
    346
    347                    // we have been forced to ignore some points to keep bracketing,
    348                    // they are probably too far from the root, drop them from now on
    349                    nbPoints = end - start;
    350                    System.arraycopy(x, start, x, 0, nbPoints);
    351                    System.arraycopy(y, start, y, 0, nbPoints);
    352                    signChangeIndex -= start;
    353
    354                } else  if (nbPoints == x.length) {
    355
    356                    // we have to drop one point in order to insert the new one
    357                    nbPoints--;
    358
    359                    // keep the tightest bracketing interval as centered as possible
    360                    if (signChangeIndex >= (x.length + 1) / 2) {
    361                        // we drop the lowest point, we have to shift the arrays and the index
    362                        System.arraycopy(x, 1, x, 0, nbPoints);
    363                        System.arraycopy(y, 1, y, 0, nbPoints);
    364                        --signChangeIndex;
    365                    }
    366
    367                }
    368
    369                // insert the last computed point
    370                //(by construction, we know it lies inside the tightest bracketing interval)
    371                System.arraycopy(x, signChangeIndex, x, signChangeIndex + 1, nbPoints - signChangeIndex);
    372                x[signChangeIndex] = nextX;
    373                System.arraycopy(y, signChangeIndex, y, signChangeIndex + 1, nbPoints - signChangeIndex);
    374                y[signChangeIndex] = nextY;
    375                ++nbPoints;
    376
    377                // update the bracketing interval
    378                if (nextY.multiply(yA).negativeOrNull()) {
    379                    // the sign change occurs before the inserted point
    380                    xB = nextX;
    381                    yB = nextY;
    382                    absYB = yB.abs();
    383                    ++agingA;
    384                    agingB = 0;
    385                } else {
    386                    // the sign change occurs after the inserted point
    387                    xA = nextX;
    388                    yA = nextY;
    389                    absYA = yA.abs();
    390                    agingA = 0;
    391                    ++agingB;
    392
    393                    // update the sign change index
    394                    signChangeIndex++;
    395
    396                }
    397
    398            }
    399
    400        }
    401
    402        /** Guess an x value by n<sup>th</sup> order inverse polynomial interpolation.
    403         * <p>
    404         * The x value is guessed by evaluating polynomial Q(y) at y = targetY, where Q
    405         * is built such that for all considered points (x<sub>i</sub>, y<sub>i</sub>),
    406         * Q(y<sub>i</sub>) = x<sub>i</sub>.
    407         * </p>
    408         * @param targetY target value for y
    409         * @param x reference points abscissas for interpolation,
    410         * note that this array <em>is</em> modified during computation
    411         * @param y reference points ordinates for interpolation
    412         * @param start start index of the points to consider (inclusive)
    413         * @param end end index of the points to consider (exclusive)
    414         * @return guessed root (will be a NaN if two points share the same y)
    415         */
    416        private Dfp guessX(final Dfp targetY, final Dfp[] x, final Dfp[] y,
    417                           final int start, final int end) {
    418
    419            // compute Q Newton coefficients by divided differences
    420            for (int i = start; i < end - 1; ++i) {
    421                final int delta = i + 1 - start;
    422                for (int j = end - 1; j > i; --j) {
    423                    x[j] = x[j].subtract(x[j-1]).divide(y[j].subtract(y[j - delta]));
    424                }
    425            }
    426
    427            // evaluate Q(targetY)
    428            Dfp x0 = targetY.getZero();
    429            for (int j = end - 1; j >= start; --j) {
    430                x0 = x[j].add(x0.multiply(targetY.subtract(y[j])));
    431            }
    432
    433            return x0;
    434
    435        }
    436
    437    }