1/*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements.  See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License.  You may obtain a copy of the License at
8 *
9 *      http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17package org.apache.commons.math.optimization.univariate;
18
19import org.apache.commons.math.FunctionEvaluationException;
20import org.apache.commons.math.exception.NotStrictlyPositiveException;
21import org.apache.commons.math.MaxIterationsExceededException;
22import org.apache.commons.math.analysis.UnivariateRealFunction;
23import org.apache.commons.math.optimization.GoalType;
24
25/**
26 * Provide an interval that brackets a local optimum of a function.
27 * This code is based on a Python implementation (from <em>SciPy</em>,
28 * module {@code optimize.py} v0.5).
29 * @version $Revision$ $Date$
30 * @since 2.2
31 */
32public class BracketFinder {
33    /** Tolerance to avoid division by zero. */
34    private static final double EPS_MIN = 1e-21;
35    /**
36     * Golden section.
37     */
38    private static final double GOLD = 1.618034;
39    /**
40     * Factor for expanding the interval.
41     */
42    private final double growLimit;
43    /**
44     * Maximum number of iterations.
45     */
46    private final int maxIterations;
47    /**
48     * Number of iterations.
49     */
50    private int iterations;
51    /**
52     * Number of function evaluations.
53     */
54    private int evaluations;
55    /**
56     * Lower bound of the bracket.
57     */
58    private double lo;
59    /**
60     * Higher bound of the bracket.
61     */
62    private double hi;
63    /**
64     * Point inside the bracket.
65     */
66    private double mid;
67    /**
68     * Function value at {@link #lo}.
69     */
70    private double fLo;
71    /**
72     * Function value at {@link #hi}.
73     */
74    private double fHi;
75    /**
76     * Function value at {@link #mid}.
77     */
78    private double fMid;
79
80    /**
81     * Constructor with default values {@code 100, 50} (see the
82     * {@link #BracketFinder(double,int) other constructor}).
83     */
84    public BracketFinder() {
85        this(100, 50);
86    }
87
88    /**
89     * Create a bracketing interval finder.
90     *
91     * @param growLimit Expanding factor.
92     * @param maxIterations Maximum number of iterations allowed for finding
93     * a bracketing interval.
94     */
95    public BracketFinder(double growLimit,
96                         int maxIterations) {
97        if (growLimit <= 0) {
98            throw new NotStrictlyPositiveException(growLimit);
99        }
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}
296