    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.optimization.univariate;
    018
    019    import org.apache.commons.math3.util.Precision;
    020    import org.apache.commons.math3.util.FastMath;
    021    import org.apache.commons.math3.exception.NumberIsTooSmallException;
    022    import org.apache.commons.math3.exception.NotStrictlyPositiveException;
    023    import org.apache.commons.math3.optimization.ConvergenceChecker;
    024    import org.apache.commons.math3.optimization.GoalType;
    025
    026    /**
    027     * For a function defined on some interval {@code (lo, hi)}, this class
    028     * finds an approximation {@code x} to the point at which the function
    029     * attains its minimum.
    030     * It implements Richard Brent's algorithm (from his book "Algorithms for
    031     * Minimization without Derivatives", p. 79) for finding minima of real
    032     * univariate functions.
    033     * <br/>
    034     * This code is an adaptation, partly based on the Python code from SciPy
    035     * (module "optimize.py" v0.5); the original algorithm is also modified
    036     * <ul>
    037     *  <li>to use an initial guess provided by the user,</li>
    038     *  <li>to ensure that the best point encountered is the one returned.</li>
    039     * </ul>
    040     *
    041     * @version $Id: BrentOptimizer.java 1422230 2012-12-15 12:11:13Z erans $
    042     * @deprecated As of 3.1 (to be removed in 4.0).
    043     * @since 2.0
    044     */
    045    @Deprecated
    046    public class BrentOptimizer extends BaseAbstractUnivariateOptimizer {
    047        /**
    048         * Golden section.
    049         */
    050        private static final double GOLDEN_SECTION = 0.5 * (3 - FastMath.sqrt(5));
    051        /**
    052         * Minimum relative tolerance.
    053         */
    054        private static final double MIN_RELATIVE_TOLERANCE = 2 * FastMath.ulp(1d);
    055        /**
    056         * Relative threshold.
    057         */
    058        private final double relativeThreshold;
    059        /**
    060         * Absolute threshold.
    061         */
    062        private final double absoluteThreshold;
    063
    064        /**
    065         * The arguments are used implement the original stopping criterion
    066         * of Brent's algorithm.
    067         * {@code abs} and {@code rel} define a tolerance
    068         * {@code tol = rel |x| + abs}. {@code rel} should be no smaller than
    069         * <em>2 macheps</em> and preferably not much less than <em>sqrt(macheps)</em>,
    070         * where <em>macheps</em> is the relative machine precision. {@code abs} must
    071         * be positive.
    072         *
    073         * @param rel Relative threshold.
    074         * @param abs Absolute threshold.
    075         * @param checker Additional, user-defined, convergence checking
    076         * procedure.
    077         * @throws NotStrictlyPositiveException if {@code abs <= 0}.
    078         * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
    079         */
    080        public BrentOptimizer(double rel,
    081                              double abs,
    082                              ConvergenceChecker<UnivariatePointValuePair> checker) {
    083            super(checker);
    084
    085            if (rel < MIN_RELATIVE_TOLERANCE) {
    086                throw new NumberIsTooSmallException(rel, MIN_RELATIVE_TOLERANCE, true);
    087            }
    088            if (abs <= 0) {
    089                throw new NotStrictlyPositiveException(abs);
    090            }
    091
    092            relativeThreshold = rel;
    093            absoluteThreshold = abs;
    094        }
    095
    096        /**
    097         * The arguments are used for implementing the original stopping criterion
    098         * of Brent's algorithm.
    099         * {@code abs} and {@code rel} define a tolerance
    100         * {@code tol = rel |x| + abs}. {@code rel} should be no smaller than
    101         * <em>2 macheps</em> and preferably not much less than <em>sqrt(macheps)</em>,
    102         * where <em>macheps</em> is the relative machine precision. {@code abs} must
    103         * be positive.
    104         *
    105         * @param rel Relative threshold.
    106         * @param abs Absolute threshold.
    107         * @throws NotStrictlyPositiveException if {@code abs <= 0}.
    108         * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
    109         */
    110        public BrentOptimizer(double rel,
    111                              double abs) {
    112            this(rel, abs, null);
    113        }
    114
    115        /** {@inheritDoc} */
    116        @Override
    117        protected UnivariatePointValuePair doOptimize() {
    118            final boolean isMinim = getGoalType() == GoalType.MINIMIZE;
    119            final double lo = getMin();
    120            final double mid = getStartValue();
    121            final double hi = getMax();
    122
    123            // Optional additional convergence criteria.
    124            final ConvergenceChecker<UnivariatePointValuePair> checker
    125                = getConvergenceChecker();
    126
    127            double a;
    128            double b;
    129            if (lo < hi) {
    130                a = lo;
    131                b = hi;
    132            } else {
    133                a = hi;
    134                b = lo;
    135            }
    136
    137            double x = mid;
    138            double v = x;
    139            double w = x;
    140            double d = 0;
    141            double e = 0;
    142            double fx = computeObjectiveValue(x);
    143            if (!isMinim) {
    144                fx = -fx;
    145            }
    146            double fv = fx;
    147            double fw = fx;
    148
    149            UnivariatePointValuePair previous = null;
    150            UnivariatePointValuePair current
    151                = new UnivariatePointValuePair(x, isMinim ? fx : -fx);
    152            // Best point encountered so far (which is the initial guess).
    153            UnivariatePointValuePair best = current;
    154
    155            int iter = 0;
    156            while (true) {
    157                final double m = 0.5 * (a + b);
    158                final double tol1 = relativeThreshold * FastMath.abs(x) + absoluteThreshold;
    159                final double tol2 = 2 * tol1;
    160
    161                // Default stopping criterion.
    162                final boolean stop = FastMath.abs(x - m) <= tol2 - 0.5 * (b - a);
    163                if (!stop) {
    164                    double p = 0;
    165                    double q = 0;
    166                    double r = 0;
    167                    double u = 0;
    168
    169                    if (FastMath.abs(e) > tol1) { // Fit parabola.
    170                        r = (x - w) * (fx - fv);
    171                        q = (x - v) * (fx - fw);
    172                        p = (x - v) * q - (x - w) * r;
    173                        q = 2 * (q - r);
    174
    175                        if (q > 0) {
    176                            p = -p;
    177                        } else {
    178                            q = -q;
    179                        }
    180
    181                        r = e;
    182                        e = d;
    183
    184                        if (p > q * (a - x) &&
    185                            p < q * (b - x) &&
    186                            FastMath.abs(p) < FastMath.abs(0.5 * q * r)) {
    187                            // Parabolic interpolation step.
    188                            d = p / q;
    189                            u = x + d;
    190
    191                            // f must not be evaluated too close to a or b.
    192                            if (u - a < tol2 || b - u < tol2) {
    193                                if (x <= m) {
    194                                    d = tol1;
    195                                } else {
    196                                    d = -tol1;
    197                                }
    198                            }
    199                        } else {
    200                            // Golden section step.
    201                            if (x < m) {
    202                                e = b - x;
    203                            } else {
    204                                e = a - x;
    205                            }
    206                            d = GOLDEN_SECTION * e;
    207                        }
    208                    } else {
    209                        // Golden section step.
    210                        if (x < m) {
    211                            e = b - x;
    212                        } else {
    213                            e = a - x;
    214                        }
    215                        d = GOLDEN_SECTION * e;
    216                    }
    217
    218                    // Update by at least "tol1".
    219                    if (FastMath.abs(d) < tol1) {
    220                        if (d >= 0) {
    221                            u = x + tol1;
    222                        } else {
    223                            u = x - tol1;
    224                        }
    225                    } else {
    226                        u = x + d;
    227                    }
    228
    229                    double fu = computeObjectiveValue(u);
    230                    if (!isMinim) {
    231                        fu = -fu;
    232                    }
    233
    234                    // User-defined convergence checker.
    235                    previous = current;
    236                    current = new UnivariatePointValuePair(u, isMinim ? fu : -fu);
    237                    best = best(best,
    238                                best(previous,
    239                                     current,
    240                                     isMinim),
    241                                isMinim);
    242
    243                    if (checker != null) {
    244                        if (checker.converged(iter, previous, current)) {
    245                            return best;
    246                        }
    247                    }
    248
    249                    // Update a, b, v, w and x.
    250                    if (fu <= fx) {
    251                        if (u < x) {
    252                            b = x;
    253                        } else {
    254                            a = x;
    255                        }
    256                        v = w;
    257                        fv = fw;
    258                        w = x;
    259                        fw = fx;
    260                        x = u;
    261                        fx = fu;
    262                    } else {
    263                        if (u < x) {
    264                            a = u;
    265                        } else {
    266                            b = u;
    267                        }
    268                        if (fu <= fw ||
    269                            Precision.equals(w, x)) {
    270                            v = w;
    271                            fv = fw;
    272                            w = u;
    273                            fw = fu;
    274                        } else if (fu <= fv ||
    275                                   Precision.equals(v, x) ||
    276                                   Precision.equals(v, w)) {
    277                            v = u;
    278                            fv = fu;
    279                        }
    280                    }
    281                } else { // Default termination (Brent's criterion).
    282                    return best(best,
    283                                best(previous,
    284                                     current,
    285                                     isMinim),
    286                                isMinim);
    287                }
    288                ++iter;
    289            }
    290        }
    291
    292        /**
    293         * Selects the best of two points.
    294         *
    295         * @param a Point and value.
    296         * @param b Point and value.
    297         * @param isMinim {@code true} if the selected point must be the one with
    298         * the lowest value.
    299         * @return the best point, or {@code null} if {@code a} and {@code b} are
    300         * both {@code null}. When {@code a} and {@code b} have the same function
    301         * value, {@code a} is returned.
    302         */
    303        private UnivariatePointValuePair best(UnivariatePointValuePair a,
    304                                              UnivariatePointValuePair b,
    305                                              boolean isMinim) {
    306            if (a == null) {
    307                return b;
    308            }
    309            if (b == null) {
    310                return a;
    311            }
    312
    313            if (isMinim) {
    314                return a.getValue() <= b.getValue() ? a : b;
    315            } else {
    316                return a.getValue() >= b.getValue() ? a : b;
    317            }
    318        }
    319    }