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.optim.univariate;
018
019import org.apache.commons.math3.util.FastMath;
020import org.apache.commons.math3.util.IntegerSequence;
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.optim.nonlinear.scalar.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 * @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     * Counter for function evaluations.
047     */
048    private IntegerSequence.Incrementor evaluations;
049    /**
050     * Lower bound of the bracket.
051     */
052    private double lo;
053    /**
054     * Higher bound of the bracket.
055     */
056    private double hi;
057    /**
058     * Point inside the bracket.
059     */
060    private double mid;
061    /**
062     * Function value at {@link #lo}.
063     */
064    private double fLo;
065    /**
066     * Function value at {@link #hi}.
067     */
068    private double fHi;
069    /**
070     * Function value at {@link #mid}.
071     */
072    private double fMid;
073
074    /**
075     * Constructor with default values {@code 100, 500} (see the
076     * {@link #BracketFinder(double,int) other constructor}).
077     */
078    public BracketFinder() {
079        this(100, 500);
080    }
081
082    /**
083     * Create a bracketing interval finder.
084     *
085     * @param growLimit Expanding factor.
086     * @param maxEvaluations Maximum number of evaluations allowed for finding
087     * a bracketing interval.
088     */
089    public BracketFinder(double growLimit,
090                         int maxEvaluations) {
091        if (growLimit <= 0) {
092            throw new NotStrictlyPositiveException(growLimit);
093        }
094        if (maxEvaluations <= 0) {
095            throw new NotStrictlyPositiveException(maxEvaluations);
096        }
097
098        this.growLimit = growLimit;
099        evaluations = IntegerSequence.Incrementor.create().withMaximalCount(maxEvaluations);
100    }
101
102    /**
103     * Search new points that bracket a local optimum of the function.
104     *
105     * @param func Function whose optimum should be bracketed.
106     * @param goal {@link GoalType Goal type}.
107     * @param xA Initial point.
108     * @param xB Initial point.
109     * @throws TooManyEvaluationsException if the maximum number of evaluations
110     * is exceeded.
111     */
112    public void search(UnivariateFunction func,
113                       GoalType goal,
114                       double xA,
115                       double xB) {
116        evaluations = evaluations.withStart(0);
117        final boolean isMinim = goal == GoalType.MINIMIZE;
118
119        double fA = eval(func, xA);
120        double fB = eval(func, xB);
121        if (isMinim ?
122            fA < fB :
123            fA > fB) {
124
125            double tmp = xA;
126            xA = xB;
127            xB = tmp;
128
129            tmp = fA;
130            fA = fB;
131            fB = tmp;
132        }
133
134        double xC = xB + GOLD * (xB - xA);
135        double fC = eval(func, xC);
136
137        while (isMinim ? fC < fB : fC > fB) {
138            double tmp1 = (xB - xA) * (fB - fC);
139            double tmp2 = (xB - xC) * (fB - fA);
140
141            double val = tmp2 - tmp1;
142            double denom = FastMath.abs(val) < EPS_MIN ? 2 * EPS_MIN : 2 * val;
143
144            double w = xB - ((xB - xC) * tmp2 - (xB - xA) * tmp1) / denom;
145            double wLim = xB + growLimit * (xC - xB);
146
147            double fW;
148            if ((w - xC) * (xB - w) > 0) {
149                fW = eval(func, w);
150                if (isMinim ?
151                    fW < fC :
152                    fW > fC) {
153                    xA = xB;
154                    xB = w;
155                    fA = fB;
156                    fB = fW;
157                    break;
158                } else if (isMinim ?
159                           fW > fB :
160                           fW < fB) {
161                    xC = w;
162                    fC = fW;
163                    break;
164                }
165                w = xC + GOLD * (xC - xB);
166                fW = eval(func, w);
167            } else if ((w - wLim) * (wLim - xC) >= 0) {
168                w = wLim;
169                fW = eval(func, w);
170            } else if ((w - wLim) * (xC - w) > 0) {
171                fW = eval(func, w);
172                if (isMinim ?
173                    fW < fC :
174                    fW > fC) {
175                    xB = xC;
176                    xC = w;
177                    w = xC + GOLD * (xC - xB);
178                    fB = fC;
179                    fC =fW;
180                    fW = eval(func, w);
181                }
182            } else {
183                w = xC + GOLD * (xC - xB);
184                fW = eval(func, w);
185            }
186
187            xA = xB;
188            fA = fB;
189            xB = xC;
190            fB = fC;
191            xC = w;
192            fC = fW;
193        }
194
195        lo = xA;
196        fLo = fA;
197        mid = xB;
198        fMid = fB;
199        hi = xC;
200        fHi = fC;
201
202        if (lo > hi) {
203            double tmp = lo;
204            lo = hi;
205            hi = tmp;
206
207            tmp = fLo;
208            fLo = fHi;
209            fHi = tmp;
210        }
211    }
212
213    /**
214     * @return the number of evalutations.
215     */
216    public int getMaxEvaluations() {
217        return evaluations.getMaximalCount();
218    }
219
220    /**
221     * @return the number of evalutations.
222     */
223    public int getEvaluations() {
224        return evaluations.getCount();
225    }
226
227    /**
228     * @return the lower bound of the bracket.
229     * @see #getFLo()
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 getFLo() {
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 TooManyEvaluationsException if the maximal number of evaluations is
280     * exceeded.
281     */
282    private double eval(UnivariateFunction f, double x) {
283        try {
284            evaluations.increment();
285        } catch (MaxCountExceededException e) {
286            throw new TooManyEvaluationsException(e.getMax());
287        }
288        return f.value(x);
289    }
290}