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.Incrementor;
20  import org.apache.commons.math3.exception.NotStrictlyPositiveException;
21  import org.apache.commons.math3.exception.TooManyEvaluationsException;
22  import org.apache.commons.math3.exception.MaxCountExceededException;
23  import org.apache.commons.math3.analysis.UnivariateFunction;
24  import org.apache.commons.math3.optim.nonlinear.scalar.GoalType;
25  
26  /**
27   * Provide an interval that brackets a local optimum of a function.
28   * This code is based on a Python implementation (from <em>SciPy</em>,
29   * module {@code optimize.py} v0.5).
30   *
31   * @version $Id: BracketFinder.java 1435539 2013-01-19 13:27:24Z tn $
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, GoalType goal, double xA, double xB) {
113         evaluations.resetCount();
114         final boolean isMinim = goal == GoalType.MINIMIZE;
115 
116         double fA = eval(func, xA);
117         double fB = eval(func, xB);
118         if (isMinim ?
119             fA < fB :
120             fA > fB) {
121 
122             double tmp = xA;
123             xA = xB;
124             xB = tmp;
125 
126             tmp = fA;
127             fA = fB;
128             fB = tmp;
129         }
130 
131         double xC = xB + GOLD * (xB - xA);
132         double fC = eval(func, xC);
133 
134         while (isMinim ? fC < fB : fC > fB) {
135             double tmp1 = (xB - xA) * (fB - fC);
136             double tmp2 = (xB - xC) * (fB - fA);
137 
138             double val = tmp2 - tmp1;
139             double denom = Math.abs(val) < EPS_MIN ? 2 * EPS_MIN : 2 * val;
140 
141             double w = xB - ((xB - xC) * tmp2 - (xB - xA) * tmp1) / denom;
142             double wLim = xB + growLimit * (xC - xB);
143 
144             double fW;
145             if ((w - xC) * (xB - w) > 0) {
146                 fW = eval(func, w);
147                 if (isMinim ?
148                     fW < fC :
149                     fW > fC) {
150                     xA = xB;
151                     xB = w;
152                     fA = fB;
153                     fB = fW;
154                     break;
155                 } else if (isMinim ?
156                            fW > fB :
157                            fW < fB) {
158                     xC = w;
159                     fC = fW;
160                     break;
161                 }
162                 w = xC + GOLD * (xC - xB);
163                 fW = eval(func, w);
164             } else if ((w - wLim) * (wLim - xC) >= 0) {
165                 w = wLim;
166                 fW = eval(func, w);
167             } else if ((w - wLim) * (xC - w) > 0) {
168                 fW = eval(func, w);
169                 if (isMinim ?
170                     fW < fC :
171                     fW > fC) {
172                     xB = xC;
173                     xC = w;
174                     w = xC + GOLD * (xC - xB);
175                     fB = fC;
176                     fC =fW;
177                     fW = eval(func, w);
178                 }
179             } else {
180                 w = xC + GOLD * (xC - xB);
181                 fW = eval(func, w);
182             }
183 
184             xA = xB;
185             fA = fB;
186             xB = xC;
187             fB = fC;
188             xC = w;
189             fC = fW;
190         }
191 
192         lo = xA;
193         fLo = fA;
194         mid = xB;
195         fMid = fB;
196         hi = xC;
197         fHi = fC;
198 
199         if (lo > hi) {
200             double tmp = lo;
201             lo = hi;
202             hi = tmp;
203 
204             tmp = fLo;
205             fLo = fHi;
206             fHi = tmp;
207         }
208     }
209 
210     /**
211      * @return the number of evalutations.
212      */
213     public int getMaxEvaluations() {
214         return evaluations.getMaximalCount();
215     }
216 
217     /**
218      * @return the number of evalutations.
219      */
220     public int getEvaluations() {
221         return evaluations.getCount();
222     }
223 
224     /**
225      * @return the lower bound of the bracket.
226      * @see #getFLo()
227      */
228     public double getLo() {
229         return lo;
230     }
231 
232     /**
233      * Get function value at {@link #getLo()}.
234      * @return function value at {@link #getLo()}
235      */
236     public double getFLo() {
237         return fLo;
238     }
239 
240     /**
241      * @return the higher bound of the bracket.
242      * @see #getFHi()
243      */
244     public double getHi() {
245         return hi;
246     }
247 
248     /**
249      * Get function value at {@link #getHi()}.
250      * @return function value at {@link #getHi()}
251      */
252     public double getFHi() {
253         return fHi;
254     }
255 
256     /**
257      * @return a point in the middle of the bracket.
258      * @see #getFMid()
259      */
260     public double getMid() {
261         return mid;
262     }
263 
264     /**
265      * Get function value at {@link #getMid()}.
266      * @return function value at {@link #getMid()}
267      */
268     public double getFMid() {
269         return fMid;
270     }
271 
272     /**
273      * @param f Function.
274      * @param x Argument.
275      * @return {@code f(x)}
276      * @throws TooManyEvaluationsException if the maximal number of evaluations is
277      * exceeded.
278      */
279     private double eval(UnivariateFunction f, double x) {
280         try {
281             evaluations.incrementCount();
282         } catch (MaxCountExceededException e) {
283             throw new TooManyEvaluationsException(e.getMax());
284         }
285         return f.value(x);
286     }
287 }