View Javadoc
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   */
17  package org.apache.commons.math3.optim.univariate;
18  
19  import org.apache.commons.math3.util.FastMath;
20  import org.apache.commons.math3.util.Incrementor;
21  import org.apache.commons.math3.exception.NotStrictlyPositiveException;
22  import org.apache.commons.math3.exception.TooManyEvaluationsException;
23  import org.apache.commons.math3.exception.MaxCountExceededException;
24  import org.apache.commons.math3.analysis.UnivariateFunction;
25  import org.apache.commons.math3.optim.nonlinear.scalar.GoalType;
26  
27  /**
28   * Provide an interval that brackets a local optimum of a function.
29   * This code is based on a Python implementation (from <em>SciPy</em>,
30   * module {@code optimize.py} v0.5).
31   *
32   * @version $Id: BracketFinder.java 1573316 2014-03-02 14:54:37Z erans $
33   * @since 2.2
34   */
35  public class BracketFinder {
36      /** Tolerance to avoid division by zero. */
37      private static final double EPS_MIN = 1e-21;
38      /**
39       * Golden section.
40       */
41      private static final double GOLD = 1.618034;
42      /**
43       * Factor for expanding the interval.
44       */
45      private final double growLimit;
46      /**
47       * Counter for function evaluations.
48       */
49      private final Incrementor evaluations = new Incrementor();
50      /**
51       * Lower bound of the bracket.
52       */
53      private double lo;
54      /**
55       * Higher bound of the bracket.
56       */
57      private double hi;
58      /**
59       * Point inside the bracket.
60       */
61      private double mid;
62      /**
63       * Function value at {@link #lo}.
64       */
65      private double fLo;
66      /**
67       * Function value at {@link #hi}.
68       */
69      private double fHi;
70      /**
71       * Function value at {@link #mid}.
72       */
73      private double fMid;
74  
75      /**
76       * Constructor with default values {@code 100, 50} (see the
77       * {@link #BracketFinder(double,int) other constructor}).
78       */
79      public BracketFinder() {
80          this(100, 50);
81      }
82  
83      /**
84       * Create a bracketing interval finder.
85       *
86       * @param growLimit Expanding factor.
87       * @param maxEvaluations Maximum number of evaluations allowed for finding
88       * a bracketing interval.
89       */
90      public BracketFinder(double growLimit,
91                           int maxEvaluations) {
92          if (growLimit <= 0) {
93              throw new NotStrictlyPositiveException(growLimit);
94          }
95          if (maxEvaluations <= 0) {
96              throw new NotStrictlyPositiveException(maxEvaluations);
97          }
98  
99          this.growLimit = growLimit;
100         evaluations.setMaximalCount(maxEvaluations);
101     }
102 
103     /**
104      * Search new points that bracket a local optimum of the function.
105      *
106      * @param func Function whose optimum should be bracketed.
107      * @param goal {@link GoalType Goal type}.
108      * @param xA Initial point.
109      * @param xB Initial point.
110      * @throws TooManyEvaluationsException if the maximum number of evaluations
111      * is exceeded.
112      */
113     public void search(UnivariateFunction func,
114                        GoalType goal,
115                        double xA,
116                        double xB) {
117         evaluations.resetCount();
118         final boolean isMinim = goal == GoalType.MINIMIZE;
119 
120         double fA = eval(func, xA);
121         double fB = eval(func, xB);
122         if (isMinim ?
123             fA < fB :
124             fA > fB) {
125 
126             double tmp = xA;
127             xA = xB;
128             xB = tmp;
129 
130             tmp = fA;
131             fA = fB;
132             fB = tmp;
133         }
134 
135         double xC = xB + GOLD * (xB - xA);
136         double fC = eval(func, xC);
137 
138         while (isMinim ? fC < fB : fC > fB) {
139             double tmp1 = (xB - xA) * (fB - fC);
140             double tmp2 = (xB - xC) * (fB - fA);
141 
142             double val = tmp2 - tmp1;
143             double denom = FastMath.abs(val) < EPS_MIN ? 2 * EPS_MIN : 2 * val;
144 
145             double w = xB - ((xB - xC) * tmp2 - (xB - xA) * tmp1) / denom;
146             double wLim = xB + growLimit * (xC - xB);
147 
148             double fW;
149             if ((w - xC) * (xB - w) > 0) {
150                 fW = eval(func, w);
151                 if (isMinim ?
152                     fW < fC :
153                     fW > fC) {
154                     xA = xB;
155                     xB = w;
156                     fA = fB;
157                     fB = fW;
158                     break;
159                 } else if (isMinim ?
160                            fW > fB :
161                            fW < fB) {
162                     xC = w;
163                     fC = fW;
164                     break;
165                 }
166                 w = xC + GOLD * (xC - xB);
167                 fW = eval(func, w);
168             } else if ((w - wLim) * (wLim - xC) >= 0) {
169                 w = wLim;
170                 fW = eval(func, w);
171             } else if ((w - wLim) * (xC - w) > 0) {
172                 fW = eval(func, w);
173                 if (isMinim ?
174                     fW < fC :
175                     fW > fC) {
176                     xB = xC;
177                     xC = w;
178                     w = xC + GOLD * (xC - xB);
179                     fB = fC;
180                     fC =fW;
181                     fW = eval(func, w);
182                 }
183             } else {
184                 w = xC + GOLD * (xC - xB);
185                 fW = eval(func, w);
186             }
187 
188             xA = xB;
189             fA = fB;
190             xB = xC;
191             fB = fC;
192             xC = w;
193             fC = fW;
194         }
195 
196         lo = xA;
197         fLo = fA;
198         mid = xB;
199         fMid = fB;
200         hi = xC;
201         fHi = fC;
202 
203         if (lo > hi) {
204             double tmp = lo;
205             lo = hi;
206             hi = tmp;
207 
208             tmp = fLo;
209             fLo = fHi;
210             fHi = tmp;
211         }
212     }
213 
214     /**
215      * @return the number of evalutations.
216      */
217     public int getMaxEvaluations() {
218         return evaluations.getMaximalCount();
219     }
220 
221     /**
222      * @return the number of evalutations.
223      */
224     public int getEvaluations() {
225         return evaluations.getCount();
226     }
227 
228     /**
229      * @return the lower bound of the bracket.
230      * @see #getFLo()
231      */
232     public double getLo() {
233         return lo;
234     }
235 
236     /**
237      * Get function value at {@link #getLo()}.
238      * @return function value at {@link #getLo()}
239      */
240     public double getFLo() {
241         return fLo;
242     }
243 
244     /**
245      * @return the higher bound of the bracket.
246      * @see #getFHi()
247      */
248     public double getHi() {
249         return hi;
250     }
251 
252     /**
253      * Get function value at {@link #getHi()}.
254      * @return function value at {@link #getHi()}
255      */
256     public double getFHi() {
257         return fHi;
258     }
259 
260     /**
261      * @return a point in the middle of the bracket.
262      * @see #getFMid()
263      */
264     public double getMid() {
265         return mid;
266     }
267 
268     /**
269      * Get function value at {@link #getMid()}.
270      * @return function value at {@link #getMid()}
271      */
272     public double getFMid() {
273         return fMid;
274     }
275 
276     /**
277      * @param f Function.
278      * @param x Argument.
279      * @return {@code f(x)}
280      * @throws TooManyEvaluationsException if the maximal number of evaluations is
281      * exceeded.
282      */
283     private double eval(UnivariateFunction f, double x) {
284         try {
285             evaluations.incrementCount();
286         } catch (MaxCountExceededException e) {
287             throw new TooManyEvaluationsException(e.getMax());
288         }
289         return f.value(x);
290     }
291 }