1 /*****************************************************************************
2 * *
3 * UNURAN -- Universal Non-Uniform Random number generator *
4 * *
5 *****************************************************************************
6 * *
7 * FILE: fmax.c *
8 * *
9 * Find maximum of a function using Brent's algorithm. *
10 * *
11 *****************************************************************************
12 * *
13 * Copyright (c) 2000-2006 Wolfgang Hoermann and Josef Leydold *
14 * Department of Statistics and Mathematics, WU Wien, Austria *
15 * *
16 * This program is free software; you can redistribute it and/or modify *
17 * it under the terms of the GNU General Public License as published by *
18 * the Free Software Foundation; either version 2 of the License, or *
19 * (at your option) any later version. *
20 * *
21 * This program is distributed in the hope that it will be useful, *
22 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
24 * GNU General Public License for more details. *
25 * *
26 * You should have received a copy of the GNU General Public License *
27 * along with this program; if not, write to the *
28 * Free Software Foundation, Inc., *
29 * 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA *
30 * *
31 *****************************************************************************/
32
33 /*---------------------------------------------------------------------------*/
34
35 #include <unur_source.h>
36 #include "fmax_source.h"
37
38 /*---------------------------------------------------------------------------*/
39
40 double
_unur_util_find_max(struct unur_funct_generic fs,double interval_min,double interval_max,double guess_max)41 _unur_util_find_max( struct unur_funct_generic fs, /* function structure */
42 double interval_min, /* lower bound of interval */
43 double interval_max, /* upper bound of interval */
44 double guess_max /* initial guess for maximum */
45 )
46 /*----------------------------------------------------------------------*/
47 /* find max of univariate continuous distribution numerically */
48 /* */
49 /* This is achieved by the following steps: */
50 /* -- Determine a interval containing the max and a third */
51 /* point within ( x0< x1 < x2 ) */
52 /* ++ Determine a region where to search the max; this will be the */
53 /* interval [max-100, max+100] if this is no contrdiction to */
54 /* the given interval;`max' is the best known approx to the max */
55 /* ++ Find a point in the interval with a positve pdf */
56 /* This is done by two geometric sequences of MAX_SRCH elements */
57 /* and find two other points within the given domain */
58 /* ++ Unbounded domains: refine x0, x1, x2 until: */
59 /* f(x0) < f(x1) > f(x2) */
60 /* -- invoke a maximization-routine to determine the max */
61 /* */
62 /* parameters: */
63 /* fs ... function structure */
64 /* interval_min ... left boundary of function definition interval */
65 /* interval_max ... right boundary of function definition interval */
66 /* guess_max ... initial guess for max position */
67 /* */
68 /* return: */
69 /* x ... approximate position of max on success */
70 /* INFINITY ... on error */
71 /*----------------------------------------------------------------------*/
72 {
73 #define MAX_SRCH (100)
74
75 int i;
76
77 double x[3]; /* max (and x[2]) should be between x[0] and x[2] ...*/
78 double fx[3]; /* ... and the respective funtion values */
79 double max; /* (approximative) max of the distribution */
80 double max_l; /* lower bound for max search */
81 double max_u; /* upper bound for max search */
82 double step;
83
84 int unbound_left;
85 int unbound_right;
86
87 /* first guess for max */
88 max = (_unur_FP_is_infinity(guess_max)) ? 0. : guess_max;
89
90 /* determine where to look for the max */
91
92 /* unbounded domain */
93 if ( _unur_FP_is_minus_infinity(interval_min) &&
94 _unur_FP_is_infinity(interval_max) ){
95
96 unbound_left = 1;
97 unbound_right = 1;
98
99 x[1] = max;
100 fx[1] = fs.f(x[1], fs.params);
101 max_l = max - 100.0;
102 max_u = max + 100.0;
103
104 }
105 /* domain unbounded on the right */
106 else if ( ! _unur_FP_is_minus_infinity(interval_min) &&
107 _unur_FP_is_infinity(interval_max) ){
108
109 unbound_left = 0;
110 unbound_right = 1;
111
112 if ( max >= interval_min ){
113 x[1] = max;
114 fx[1] = fs.f(x[1], fs.params);
115 max_l = interval_min;
116 max_u = 2 * max - interval_min;
117 }
118 else{
119 x[1] = interval_min + 100.0;
120 fx[1] = fs.f(x[1], fs.params);
121 max_l = interval_min;
122 max_u = x[1] + 100.0;
123 }
124
125 }
126 /* domain unbounded on the left */
127 else if ( _unur_FP_is_minus_infinity(interval_min) &&
128 ! _unur_FP_is_infinity(interval_max) ){
129
130 unbound_left = 1;
131 unbound_right = 0;
132
133 if ( max <= interval_max ){
134 x[1] = max;
135 fx[1] = fs.f(x[1], fs.params);
136 max_l = interval_max - 2 * max;
137 max_u = interval_max;
138 }
139 else{
140 x[1] = interval_max - 100.0;
141 fx[1] = fs.f(x[1], fs.params);
142 max_l = x[1] - 100.0;
143 max_u = interval_max;
144 }
145
146 }
147 /* domain is bounded */
148 else {
149
150 unbound_left = 0;
151 unbound_right = 0;
152
153 if ( max >= interval_min && max <= interval_max ){
154 x[1] = max;
155 fx[1] = fs.f(x[1], fs.params);
156 }
157 else{
158 x[1] = interval_min/2.0 + interval_max/2.0;
159 fx[1] = fs.f(x[1], fs.params);
160 }
161 max_l = interval_min;
162 max_u = interval_max;
163
164 }
165
166 max = x[1]; /* not exact mode -- best guess */
167
168
169 /* find point with pdf > 0.0 -- max MAX_SRCH trials */
170
171 /* search on the left side */
172 step = pow(x[1]-max_l, 1.0/MAX_SRCH);
173 i = 0;
174 while (i <= MAX_SRCH && _unur_FP_same(0.0, fx[1]) ){
175 x[1] = max - pow(step, (double)i);
176 fx[1] = fs.f(x[1], fs.params);
177 i++;
178 }
179
180 /* search on the right side */
181 if( _unur_FP_same(0.0, fx[1]) ){
182 step = pow(max_u-x[1], 1.0/MAX_SRCH);
183 i = 0;
184 while (i <= MAX_SRCH && _unur_FP_same(0.0, fx[1]) ){
185 x[1] = max + pow(step, (double)i);
186 fx[1] = fs.f(x[1], fs.params);
187 i++;
188 }
189 }
190
191 /* no success -- exit routine */
192 if( _unur_FP_same(fx[1], 0.0) )
193 return INFINITY; /* can't find max in flat region */
194
195 /* x[1] has f > 0 or routines already terminated */
196
197
198
199 /* determine 3 points in the given domain --
200 at least one with pdf > 0 */
201 if ( unbound_left ){
202
203 x[2] = x[1]; fx[2] = fx[1];
204 x[1] = x[2] - 1.0; fx[1] = fs.f(x[1], fs.params);
205 x[0] = x[2] - 2.0; fx[0] = fs.f(x[0], fs.params);
206
207 }
208 else if ( unbound_right ){
209
210 x[0] = x[1]; fx[0] = fx[1];
211 x[1] = x[0] + 1.0; fx[1] = fs.f(x[1], fs.params);
212 x[2] = x[0] + 2.0; fx[2] = fs.f(x[2], fs.params);
213
214 }
215 else{ /* bounded */
216
217 x[0] = interval_min; fx[0] = fs.f(x[0], fs.params);
218 x[2] = interval_max; fx[2] = fs.f(x[2], fs.params);
219
220 if ( _unur_FP_same(x[1], interval_min) ||
221 _unur_FP_same(x[1], interval_max) ){
222 x[1] = interval_min/2.0 + interval_max/2.0;
223 fx[1] = fs.f(x[1], fs.params);
224 }
225
226 }
227 /* points x[i] with their function values determined */
228
229 /* find interval containing the max */
230
231 step = 1.0;
232 if ( unbound_right ){
233 while(fx[0] <= fx[1] && fx[1] <= fx[2]){ /* on the left side of the max */
234
235 step *= 2.0;
236 x[0] = x[1]; fx[0] = fx[1];
237 x[1] = x[2]; fx[1] = fx[2];
238 x[2] += step; fx[2] = fs.f(x[2], fs.params);
239 }
240 }
241
242 step = 1.0; /* reset step size */
243 if ( unbound_left ){
244 while(fx[0] >= fx[1] && fx[1] >= fx[2]){ /* on the right side of the max */
245
246 step *= 2.0;
247 x[2] = x[1]; fx[2] = fx[1];
248 x[1] = x[0]; fx[1] = fx[0];
249 x[0] -= step; fx[0] = fs.f(x[0], fs.params);
250
251 }
252 }
253
254 /* now: the max is between x[0] and x[2] */
255
256 /* printf("x0: %f, fx0: %e\n", x[0], fx[0]); */
257 /* printf("x1: %f, fx1: %e\n", x[1], fx[1]); */
258 /* printf("x2: %f, fx2: %e\n", x[2], fx[2]); */
259
260 /** TODO: FLT_MIN must be much larger than DBL_MIN **/
261
262 max = _unur_util_brent( fs, x[0], x[2], x[1], FLT_MIN );
263 if (!(_unur_FP_is_infinity( max )) ){
264 /* mode successfully computed */
265
266 }
267 else {
268 /* computing max did not work */
269 return INFINITY;
270 }
271
272 /* o.k. */
273 return max;
274
275 #undef MAX_SRCH
276 } /* end of _unur_util_find_max() */
277
278 /*---------------------------------------------------------------------------*/
279 /* */
280 /* Brent algorithm for maximum-calculation of a continous function */
281 /* */
282 /*---------------------------------------------------------------------------*/
283 /*
284 *****************************************************************************
285 * C math library *
286 * function FMINBR - one-dimensional search for a function minimum *
287 * over the given range *
288 * *
289 * Author: Oleg Keselyov. *
290 * *
291 * modified by Josef Leydold (documentation unchanged) *
292 * computed maximum instead of minimum. *
293 * *
294 * Input *
295 * double fminbr(f, a,b,c,tol) *
296 * double a; Minimum will be seeked for over *
297 * double b; a range [a,b], a being < b. *
298 * double c; c within (a,b) is first guess *
299 * double (*f)(double x); Name of the function whose minimum *
300 * will be seeked for *
301 * double tol; Acceptable tolerance for the minimum *
302 * location. It have to be positive *
303 * (e.g. may be specified as EPSILON) *
304 * *
305 * Output *
306 * Fminbr returns an estimate for the minimum location with accuracy *
307 * 3*SQRT_EPSILON*abs(x) + tol. *
308 * The function always obtains a local minimum which coincides with *
309 * the global one only if a function under investigation being *
310 * unimodular. *
311 * If a function being examined possesses no local minimum within *
312 * the given range, Fminbr returns 'a' (if f(a) < f(b)), otherwise *
313 * it returns the right range boundary value b. *
314 * *
315 * Algorithm *
316 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical *
317 * computations. M., Mir, 1980, p.202 of the Russian edition *
318 * *
319 * The function makes use of the "gold section" procedure combined with *
320 * the parabolic interpolation. *
321 * At every step program operates three abscissae - x,v, and w. *
322 * x - the last and the best approximation to the minimum location, *
323 * i.e. f(x) <= f(a) or/and f(x) <= f(b) *
324 * (if the function f has a local minimum in (a,b), then the both *
325 * conditions are fulfiled after one or two steps). *
326 * v,w are previous approximations to the minimum location. They may *
327 * coincide with a, b, or x (although the algorithm tries to make all *
328 * u, v, and w distinct). Points x, v, and w are used to construct *
329 * interpolating parabola whose minimum will be treated as a new *
330 * approximation to the minimum location if the former falls within *
331 * [a,b] and reduces the range enveloping minimum more efficient than *
332 * the gold section procedure. *
333 * When f(x) has a second derivative positive at the minimum location *
334 * (not coinciding with a or b) the procedure converges superlinearly *
335 * at a rate order about 1.324 *
336 * *
337 *****************************************************************************/
338
339 /* in case of any error INFINITY is returned */
340
341 double
_unur_util_brent(struct unur_funct_generic fs,double a,double b,double c,double tol)342 _unur_util_brent( /* An estimate to the min or max location */
343 struct unur_funct_generic fs, /* Function struct */
344 double a, /* Left border | of the range */
345 double b, /* Right border| the min is seeked */
346 double c, /* first guess for the min/max */
347 double tol) /* Acceptable tolerance */
348 {
349 #define SQRT_EPSILON (1.e-7) /* tolerance for relative error */
350 #define MAXIT (1000) /* maximum number of iterations */
351
352 /* compute maximum (instead of minimum as in original implementation ) */
353 #define f(x) ( (-1) * ((fs.f)(x, fs.params)) )
354
355
356 int i;
357
358 double x,v,w; /* Abscissae, descr. see above */
359 double fx; /* f(x) */
360 double fv; /* f(v) */
361 double fw; /* f(w) */
362 const double r = (3.-sqrt(5.0))/2; /* Gold section ratio */
363
364 /* check arguments */
365 CHECK_NULL(fs.f,INFINITY);
366 if ( tol < 0. || b <= a || c <= a || b <= c) {
367 _unur_error("CMAX",UNUR_ERR_SHOULD_NOT_HAPPEN,"");
368 return INFINITY;
369 }
370
371
372 /* Origially the third point was computed by golden section. In the */
373 /* modified version it is given as point `c' by the calling function. */
374 /* v = a + r*(b-a); fv = f(v); First step - always gold section */
375
376 v = c; fv = f(v); /* First step */
377 x = v; w = v;
378 fx=fv; fw=fv;
379
380 for(i=0; i < MAXIT; i++) { /* Main iteration loop */
381 double range = b-a; /* Range over which the minimum is */
382 /* seeked for */
383 double middle_range = (a+b)/2;
384 double tol_act = /* Actual tolerance */
385 SQRT_EPSILON*fabs(x) + tol/3;
386 double new_step; /* Step at this iteration */
387
388 if( fabs(x-middle_range) + range/2 <= 2*tol_act )
389 return x; /* Acceptable approx. is found */
390
391 /* Obtain the gold section step */
392 new_step = r * ( x<middle_range ? b-x : a-x );
393
394
395 /* Decide if the interpolation can be tried */
396 if( fabs(x-w) >= tol_act ) { /* If x and w are distinct */
397 /* interpolatiom may be tried */
398 register double p; /* Interpolation step is calculated */
399 register double q; /* as p/q; division operation */
400 /* is delayed until last moment */
401 register double t;
402
403 t = (x-w) * (fx-fv);
404 q = (x-v) * (fx-fw);
405 p = (x-v)*q - (x-w)*t;
406 q = 2*(q-t);
407
408 if( q>(double)0 ) /* q was calculated with the */
409 p = -p; /* opposite sign; make q positive */
410 else /* and assign possible minus to p */
411 q = -q;
412
413 if( fabs(p) < fabs(new_step*q) && /* If x + p/q falls in [a,b] not too */
414 p > q*(a-x+2*tol_act) && /* close to a and b, and isn't too */
415 p < q*(b-x-2*tol_act) ) /* large, it is accepted. */
416 new_step = p/q; /* If p/q is too large then the */
417 /* gold section procedure can reduce */
418 /* [a,b] range to more extent. */
419 }
420
421 if( fabs(new_step) < tol_act ) { /* Adjust the step to be not less */
422 if( new_step > (double)0 ) /* than tolerance. */
423 new_step = tol_act;
424 else
425 new_step = -tol_act;
426 }
427 /* Obtain the next approximation to min */
428 { /* and reduce the enveloping range. */
429 register double t = x + new_step; /* Tentative point for the min */
430 register double ft = f(t);
431
432 if( ft <= fx ) { /* t is a better approximation */
433 if( t < x )
434 b = x; /* Reduce the range so that */
435 else /* t would fall within it */
436 a = x;
437
438 v = w; w = x; x = t; /* Assign the best approx to x */
439 fv=fw; fw=fx; fx=ft;
440 }
441 else { /* x remains the better approx */
442 if( t < x )
443 a = t; /* Reduce the range enclosing x */
444 else
445 b = t;
446
447 if( ft <= fw || _unur_FP_same(w,x) ) {
448 v = w; w = t;
449 fv=fw; fw=ft;
450 }
451 else if( ft<=fv || _unur_FP_same(v,x) || _unur_FP_same(v,w) ) {
452 v = t;
453 fv=ft;
454 }
455 }
456 } /* ----- end-of-block ----- */
457 } /* ===== End of for loop ===== */
458
459 /* maximal number of iterations exceeded */
460 return INFINITY;
461
462 #undef f
463 #undef MAXIT
464 #undef SQRT_EPSILON
465 } /* end of _unur_util_brent() */
466
467 /*---------------------------------------------------------------------------*/
468