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    }