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