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