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