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   * @since 2.2
33   */
34  public class BracketFinder {
35      /** Tolerance to avoid division by zero. */
36      private static final double EPS_MIN = 1e-21;
37      /**
38       * Golden section.
39       */
40      private static final double GOLD = 1.618034;
41      /**
42       * Factor for expanding the interval.
43       */
44      private final double growLimit;
45      /**
46       * Counter for function evaluations.
47       */
48      private final Incrementor evaluations = new Incrementor();
49      /**
50       * Lower bound of the bracket.
51       */
52      private double lo;
53      /**
54       * Higher bound of the bracket.
55       */
56      private double hi;
57      /**
58       * Point inside the bracket.
59       */
60      private double mid;
61      /**
62       * Function value at {@link #lo}.
63       */
64      private double fLo;
65      /**
66       * Function value at {@link #hi}.
67       */
68      private double fHi;
69      /**
70       * Function value at {@link #mid}.
71       */
72      private double fMid;
73  
74      /**
75       * Constructor with default values {@code 100, 50} (see the
76       * {@link #BracketFinder(double,int) other constructor}).
77       */
78      public BracketFinder() {
79          this(100, 50);
80      }
81  
82      /**
83       * Create a bracketing interval finder.
84       *
85       * @param growLimit Expanding factor.
86       * @param maxEvaluations Maximum number of evaluations allowed for finding
87       * a bracketing interval.
88       */
89      public BracketFinder(double growLimit,
90                           int maxEvaluations) {
91          if (growLimit <= 0) {
92              throw new NotStrictlyPositiveException(growLimit);
93          }
94          if (maxEvaluations <= 0) {
95              throw new NotStrictlyPositiveException(maxEvaluations);
96          }
97  
98          this.growLimit = growLimit;
99          evaluations.setMaximalCount(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.resetCount();
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.incrementCount();
285         } catch (MaxCountExceededException e) {
286             throw new TooManyEvaluationsException(e.getMax());
287         }
288         return f.value(x);
289     }
290 }