    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.optim.univariate;
    018
    019    import org.apache.commons.math3.util.Incrementor;
    020    import org.apache.commons.math3.exception.NotStrictlyPositiveException;
    021    import org.apache.commons.math3.exception.TooManyEvaluationsException;
    022    import org.apache.commons.math3.exception.MaxCountExceededException;
    023    import org.apache.commons.math3.analysis.UnivariateFunction;
    024    import org.apache.commons.math3.optim.nonlinear.scalar.GoalType;
    025
    026    /**
    027     * Provide an interval that brackets a local optimum of a function.
    028     * This code is based on a Python implementation (from <em>SciPy</em>,
    029     * module {@code optimize.py} v0.5).
    030     *
    031     * @version $Id: BracketFinder.java 1413186 2012-11-24 13:47:59Z erans $
    032     * @since 2.2
    033     */
    034    public class BracketFinder {
    035        /** Tolerance to avoid division by zero. */
    036        private static final double EPS_MIN = 1e-21;
    037        /**
    038         * Golden section.
    039         */
    040        private static final double GOLD = 1.618034;
    041        /**
    042         * Factor for expanding the interval.
    043         */
    044        private final double growLimit;
    045        /**
    046         * Counter for function evaluations.
    047         */
    048        private final Incrementor evaluations = new Incrementor();
    049        /**
    050         * Lower bound of the bracket.
    051         */
    052        private double lo;
    053        /**
    054         * Higher bound of the bracket.
    055         */
    056        private double hi;
    057        /**
    058         * Point inside the bracket.
    059         */
    060        private double mid;
    061        /**
    062         * Function value at {@link #lo}.
    063         */
    064        private double fLo;
    065        /**
    066         * Function value at {@link #hi}.
    067         */
    068        private double fHi;
    069        /**
    070         * Function value at {@link #mid}.
    071         */
    072        private double fMid;
    073
    074        /**
    075         * Constructor with default values {@code 100, 50} (see the
    076         * {@link #BracketFinder(double,int) other constructor}).
    077         */
    078        public BracketFinder() {
    079            this(100, 50);
    080        }
    081
    082        /**
    083         * Create a bracketing interval finder.
    084         *
    085         * @param growLimit Expanding factor.
    086         * @param maxEvaluations Maximum number of evaluations allowed for finding
    087         * a bracketing interval.
    088         */
    089        public BracketFinder(double growLimit,
    090                             int maxEvaluations) {
    091            if (growLimit <= 0) {
    092                throw new NotStrictlyPositiveException(growLimit);
    093            }
    094            if (maxEvaluations <= 0) {
    095                throw new NotStrictlyPositiveException(maxEvaluations);
    096            }
    097
    098            this.growLimit = growLimit;
    099            evaluations.setMaximalCount(maxEvaluations);
    100        }
    101
    102        /**
    103         * Search new points that bracket a local optimum of the function.
    104         *
    105         * @param func Function whose optimum should be bracketed.
    106         * @param goal {@link GoalType Goal type}.
    107         * @param xA Initial point.
    108         * @param xB Initial point.
    109         * @throws TooManyEvaluationsException if the maximum number of evaluations
    110         * is exceeded.
    111         */
    112        public void search(UnivariateFunction func, GoalType goal, double xA, double xB) {
    113            evaluations.resetCount();
    114            final boolean isMinim = goal == GoalType.MINIMIZE;
    115
    116            double fA = eval(func, xA);
    117            double fB = eval(func, xB);
    118            if (isMinim ?
    119                fA < fB :
    120                fA > fB) {
    121
    122                double tmp = xA;
    123                xA = xB;
    124                xB = tmp;
    125
    126                tmp = fA;
    127                fA = fB;
    128                fB = tmp;
    129            }
    130
    131            double xC = xB + GOLD * (xB - xA);
    132            double fC = eval(func, xC);
    133
    134            while (isMinim ? fC < fB : fC > fB) {
    135                double tmp1 = (xB - xA) * (fB - fC);
    136                double tmp2 = (xB - xC) * (fB - fA);
    137
    138                double val = tmp2 - tmp1;
    139                double denom = Math.abs(val) < EPS_MIN ? 2 * EPS_MIN : 2 * val;
    140
    141                double w = xB - ((xB - xC) * tmp2 - (xB - xA) * tmp1) / denom;
    142                double wLim = xB + growLimit * (xC - xB);
    143
    144                double fW;
    145                if ((w - xC) * (xB - w) > 0) {
    146                    fW = eval(func, w);
    147                    if (isMinim ?
    148                        fW < fC :
    149                        fW > fC) {
    150                        xA = xB;
    151                        xB = w;
    152                        fA = fB;
    153                        fB = fW;
    154                        break;
    155                    } else if (isMinim ?
    156                               fW > fB :
    157                               fW < fB) {
    158                        xC = w;
    159                        fC = fW;
    160                        break;
    161                    }
    162                    w = xC + GOLD * (xC - xB);
    163                    fW = eval(func, w);
    164                } else if ((w - wLim) * (wLim - xC) >= 0) {
    165                    w = wLim;
    166                    fW = eval(func, w);
    167                } else if ((w - wLim) * (xC - w) > 0) {
    168                    fW = eval(func, w);
    169                    if (isMinim ?
    170                        fW < fC :
    171                        fW > fC) {
    172                        xB = xC;
    173                        xC = w;
    174                        w = xC + GOLD * (xC - xB);
    175                        fB = fC;
    176                        fC =fW;
    177                        fW = eval(func, w);
    178                    }
    179                } else {
    180                    w = xC + GOLD * (xC - xB);
    181                    fW = eval(func, w);
    182                }
    183
    184                xA = xB;
    185                fA = fB;
    186                xB = xC;
    187                fB = fC;
    188                xC = w;
    189                fC = fW;
    190            }
    191
    192            lo = xA;
    193            fLo = fA;
    194            mid = xB;
    195            fMid = fB;
    196            hi = xC;
    197            fHi = fC;
    198
    199            if (lo > hi) {
    200                double tmp = lo;
    201                lo = hi;
    202                hi = tmp;
    203
    204                tmp = fLo;
    205                fLo = fHi;
    206                fHi = tmp;
    207            }
    208        }
    209
    210        /**
    211         * @return the number of evalutations.
    212         */
    213        public int getMaxEvaluations() {
    214            return evaluations.getMaximalCount();
    215        }
    216
    217        /**
    218         * @return the number of evalutations.
    219         */
    220        public int getEvaluations() {
    221            return evaluations.getCount();
    222        }
    223
    224        /**
    225         * @return the lower bound of the bracket.
    226         * @see #getFLo()
    227         */
    228        public double getLo() {
    229            return lo;
    230        }
    231
    232        /**
    233         * Get function value at {@link #getLo()}.
    234         * @return function value at {@link #getLo()}
    235         */
    236        public double getFLo() {
    237            return fLo;
    238        }
    239
    240        /**
    241         * @return the higher bound of the bracket.
    242         * @see #getFHi()
    243         */
    244        public double getHi() {
    245            return hi;
    246        }
    247
    248        /**
    249         * Get function value at {@link #getHi()}.
    250         * @return function value at {@link #getHi()}
    251         */
    252        public double getFHi() {
    253            return fHi;
    254        }
    255
    256        /**
    257         * @return a point in the middle of the bracket.
    258         * @see #getFMid()
    259         */
    260        public double getMid() {
    261            return mid;
    262        }
    263
    264        /**
    265         * Get function value at {@link #getMid()}.
    266         * @return function value at {@link #getMid()}
    267         */
    268        public double getFMid() {
    269            return fMid;
    270        }
    271
    272        /**
    273         * @param f Function.
    274         * @param x Argument.
    275         * @return {@code f(x)}
    276         * @throws TooManyEvaluationsException if the maximal number of evaluations is
    277         * exceeded.
    278         */
    279        private double eval(UnivariateFunction f, double x) {
    280            try {
    281                evaluations.incrementCount();
    282            } catch (MaxCountExceededException e) {
    283                throw new TooManyEvaluationsException(e.getMax());
    284            }
    285            return f.value(x);
    286        }
    287    }