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