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    }