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 */
017package org.apache.commons.math3.optimization.univariate;
018
019import org.apache.commons.math3.util.FastMath;
020import org.apache.commons.math3.util.Incrementor;
021import org.apache.commons.math3.exception.NotStrictlyPositiveException;
022import org.apache.commons.math3.exception.TooManyEvaluationsException;
023import org.apache.commons.math3.exception.MaxCountExceededException;
024import org.apache.commons.math3.analysis.UnivariateFunction;
025import org.apache.commons.math3.optimization.GoalType;
026
027/**
028 * Provide an interval that brackets a local optimum of a function.
029 * This code is based on a Python implementation (from <em>SciPy</em>,
030 * module {@code optimize.py} v0.5).
031 *
032 * @deprecated As of 3.1 (to be removed in 4.0).
033 * @since 2.2
034 */
035@Deprecated
036public 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 = FastMath.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}