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