1 /* monte/miser.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Michael Booth
4 * Copyright (C) 2009 Brian Gough
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 3 of the License, or (at
9 * your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful, but
12 * WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
19 */
20
21 /* MISER. Based on the algorithm described in the following article,
22
23 W.H. Press, G.R. Farrar, "Recursive Stratified Sampling for
24 Multidimensional Monte Carlo Integration", Computers in Physics,
25 v4 (1990), pp190-195.
26
27 */
28
29 /* Author: MJB */
30 /* Modified by Brian Gough 12/2000 */
31
32 #include <config.h>
33 #include <math.h>
34 #include <stdlib.h>
35
36 #include <gsl/gsl_math.h>
37 #include <gsl/gsl_errno.h>
38 #include <gsl/gsl_rng.h>
39 #include <gsl/gsl_monte.h>
40 #include <gsl/gsl_monte_miser.h>
41
42 static int
43 estimate_corrmc (gsl_monte_function * f,
44 const double xl[], const double xu[],
45 size_t dim, size_t calls,
46 gsl_rng * r,
47 gsl_monte_miser_state * state,
48 double *result, double *abserr,
49 const double xmid[], double sigma_l[], double sigma_r[]);
50
51
52 int
gsl_monte_miser_integrate(gsl_monte_function * f,const double xl[],const double xu[],size_t dim,size_t calls,gsl_rng * r,gsl_monte_miser_state * state,double * result,double * abserr)53 gsl_monte_miser_integrate (gsl_monte_function * f,
54 const double xl[], const double xu[],
55 size_t dim, size_t calls,
56 gsl_rng * r,
57 gsl_monte_miser_state * state,
58 double *result, double *abserr)
59 {
60 size_t n, estimate_calls, calls_l, calls_r;
61 const size_t min_calls = state->min_calls;
62 size_t i;
63 size_t i_bisect;
64 int found_best;
65
66 double res_est = 0, err_est = 0;
67 double res_r = 0, err_r = 0, res_l = 0, err_l = 0;
68 double xbi_l, xbi_m, xbi_r, s;
69
70 double vol;
71 double weight_l, weight_r;
72
73 double *x = state->x;
74 double *xmid = state->xmid;
75 double *sigma_l = state->sigma_l, *sigma_r = state->sigma_r;
76
77 if (dim != state->dim)
78 {
79 GSL_ERROR ("number of dimensions must match allocated size", GSL_EINVAL);
80 }
81
82 for (i = 0; i < dim; i++)
83 {
84 if (xu[i] <= xl[i])
85 {
86 GSL_ERROR ("xu must be greater than xl", GSL_EINVAL);
87 }
88
89 if (xu[i] - xl[i] > GSL_DBL_MAX)
90 {
91 GSL_ERROR ("Range of integration is too large, please rescale",
92 GSL_EINVAL);
93 }
94 }
95
96 if (state->alpha < 0)
97 {
98 GSL_ERROR ("alpha must be non-negative", GSL_EINVAL);
99 }
100
101 /* Compute volume */
102
103 vol = 1;
104
105 for (i = 0; i < dim; i++)
106 {
107 vol *= xu[i] - xl[i];
108 }
109
110 if (calls < state->min_calls_per_bisection)
111 {
112 double m = 0.0, q = 0.0;
113
114 if (calls < 2)
115 {
116 GSL_ERROR ("insufficient calls for subvolume", GSL_EFAILED);
117 }
118
119 for (n = 0; n < calls; n++)
120 {
121 /* Choose a random point in the integration region */
122
123 for (i = 0; i < dim; i++)
124 {
125 x[i] = xl[i] + gsl_rng_uniform_pos (r) * (xu[i] - xl[i]);
126 }
127
128 {
129 double fval = GSL_MONTE_FN_EVAL (f, x);
130
131 /* recurrence for mean and variance */
132
133 double d = fval - m;
134 m += d / (n + 1.0);
135 q += d * d * (n / (n + 1.0));
136 }
137 }
138
139 *result = vol * m;
140
141 *abserr = vol * sqrt (q / (calls * (calls - 1.0)));
142
143 return GSL_SUCCESS;
144 }
145
146 estimate_calls = GSL_MAX (min_calls, calls * (state->estimate_frac));
147
148 if (estimate_calls < 4 * dim)
149 {
150 GSL_ERROR ("insufficient calls to sample all halfspaces", GSL_ESANITY);
151 }
152
153 /* Flip coins to bisect the integration region with some fuzz */
154
155 for (i = 0; i < dim; i++)
156 {
157 s = (gsl_rng_uniform (r) - 0.5) >= 0.0 ? state->dither : -state->dither;
158 state->xmid[i] = (0.5 + s) * xl[i] + (0.5 - s) * xu[i];
159 }
160
161 /* The idea is to chose the direction to bisect based on which will
162 give the smallest total variance. We could (and may do so later)
163 use MC to compute these variances. But the NR guys simply estimate
164 the variances by finding the min and max function values
165 for each half-region for each bisection. */
166
167 estimate_corrmc (f, xl, xu, dim, estimate_calls,
168 r, state, &res_est, &err_est, xmid, sigma_l, sigma_r);
169
170 /* We have now used up some calls for the estimation */
171
172 calls -= estimate_calls;
173
174 /* Now find direction with the smallest total "variance" */
175
176 {
177 double best_var = GSL_DBL_MAX;
178 double beta = 2.0 / (1.0 + state->alpha);
179 found_best = 0;
180 i_bisect = 0;
181 weight_l = weight_r = 1.0;
182
183 for (i = 0; i < dim; i++)
184 {
185 if (sigma_l[i] >= 0 && sigma_r[i] >= 0)
186 {
187 /* estimates are okay */
188 double var = pow (sigma_l[i], beta) + pow (sigma_r[i], beta);
189
190 if (var <= best_var)
191 {
192 found_best = 1;
193 best_var = var;
194 i_bisect = i;
195 weight_l = pow (sigma_l[i], beta);
196 weight_r = pow (sigma_r[i], beta);
197
198 if (weight_l == 0 && weight_r == 0)
199 {
200 weight_l = 1;
201 weight_r = 1;
202 }
203 }
204 }
205 else
206 {
207 if (sigma_l[i] < 0)
208 {
209 GSL_ERROR ("no points in left-half space!", GSL_ESANITY);
210 }
211 if (sigma_r[i] < 0)
212 {
213 GSL_ERROR ("no points in right-half space!", GSL_ESANITY);
214 }
215 }
216 }
217 }
218
219 if (!found_best)
220 {
221 /* All estimates were the same, so chose a direction at random */
222
223 i_bisect = gsl_rng_uniform_int (r, dim);
224 }
225
226 xbi_l = xl[i_bisect];
227 xbi_m = xmid[i_bisect];
228 xbi_r = xu[i_bisect];
229
230 /* Get the actual fractional sizes of the two "halves", and
231 distribute the remaining calls among them */
232
233 {
234 double fraction_l = fabs ((xbi_m - xbi_l) / (xbi_r - xbi_l));
235 double fraction_r = 1 - fraction_l;
236
237 double a = fraction_l * weight_l;
238 double b = fraction_r * weight_r;
239
240 calls_l = min_calls + (calls - 2 * min_calls) * a / (a + b);
241 calls_r = min_calls + (calls - 2 * min_calls) * b / (a + b);
242 }
243
244 /* Compute the integral for the left hand side of the bisection */
245
246 /* Due to the recursive nature of the algorithm we must allocate
247 some new memory for each recursive call */
248
249 {
250 int status;
251
252 double *xu_tmp = (double *) malloc (dim * sizeof (double));
253
254 if (xu_tmp == 0)
255 {
256 GSL_ERROR_VAL ("out of memory for left workspace", GSL_ENOMEM, 0);
257 }
258
259 for (i = 0; i < dim; i++)
260 {
261 xu_tmp[i] = xu[i];
262 }
263
264 xu_tmp[i_bisect] = xbi_m;
265
266 status = gsl_monte_miser_integrate (f, xl, xu_tmp,
267 dim, calls_l, r, state,
268 &res_l, &err_l);
269 free (xu_tmp);
270
271 if (status != GSL_SUCCESS)
272 {
273 return status;
274 }
275 }
276
277 /* Compute the integral for the right hand side of the bisection */
278
279 {
280 int status;
281
282 double *xl_tmp = (double *) malloc (dim * sizeof (double));
283
284 if (xl_tmp == 0)
285 {
286 GSL_ERROR_VAL ("out of memory for right workspace", GSL_ENOMEM, 0);
287 }
288
289 for (i = 0; i < dim; i++)
290 {
291 xl_tmp[i] = xl[i];
292 }
293
294 xl_tmp[i_bisect] = xbi_m;
295
296 status = gsl_monte_miser_integrate (f, xl_tmp, xu,
297 dim, calls_r, r, state,
298 &res_r, &err_r);
299 free (xl_tmp);
300
301 if (status != GSL_SUCCESS)
302 {
303 return status;
304 }
305 }
306
307 *result = res_l + res_r;
308 *abserr = sqrt (err_l * err_l + err_r * err_r);
309
310 return GSL_SUCCESS;
311 }
312
313 gsl_monte_miser_state *
gsl_monte_miser_alloc(size_t dim)314 gsl_monte_miser_alloc (size_t dim)
315 {
316 gsl_monte_miser_state *s =
317 (gsl_monte_miser_state *) malloc (sizeof (gsl_monte_miser_state));
318
319 if (s == 0)
320 {
321 GSL_ERROR_VAL ("failed to allocate space for miser state struct",
322 GSL_ENOMEM, 0);
323 }
324
325 s->x = (double *) malloc (dim * sizeof (double));
326
327 if (s->x == 0)
328 {
329 free (s);
330 GSL_ERROR_VAL ("failed to allocate space for x", GSL_ENOMEM, 0);
331 }
332
333 s->xmid = (double *) malloc (dim * sizeof (double));
334
335 if (s->xmid == 0)
336 {
337 free (s->x);
338 free (s);
339 GSL_ERROR_VAL ("failed to allocate space for xmid", GSL_ENOMEM, 0);
340 }
341
342 s->sigma_l = (double *) malloc (dim * sizeof (double));
343
344 if (s->sigma_l == 0)
345 {
346 free (s->xmid);
347 free (s->x);
348 free (s);
349 GSL_ERROR_VAL ("failed to allocate space for sigma_l", GSL_ENOMEM, 0);
350 }
351
352 s->sigma_r = (double *) malloc (dim * sizeof (double));
353
354 if (s->sigma_r == 0)
355 {
356 free (s->sigma_l);
357 free (s->xmid);
358 free (s->x);
359 free (s);
360 GSL_ERROR_VAL ("failed to allocate space for sigma_r", GSL_ENOMEM, 0);
361 }
362
363 s->fmax_l = (double *) malloc (dim * sizeof (double));
364
365 if (s->fmax_l == 0)
366 {
367 free (s->sigma_r);
368 free (s->sigma_l);
369 free (s->xmid);
370 free (s->x);
371 free (s);
372 GSL_ERROR_VAL ("failed to allocate space for fmax_l", GSL_ENOMEM, 0);
373 }
374
375 s->fmax_r = (double *) malloc (dim * sizeof (double));
376
377 if (s->fmax_r == 0)
378 {
379 free (s->fmax_l);
380 free (s->sigma_r);
381 free (s->sigma_l);
382 free (s->xmid);
383 free (s->x);
384 free (s);
385 GSL_ERROR_VAL ("failed to allocate space for fmax_r", GSL_ENOMEM, 0);
386 }
387
388 s->fmin_l = (double *) malloc (dim * sizeof (double));
389
390 if (s->fmin_l == 0)
391 {
392 free (s->fmax_r);
393 free (s->fmax_l);
394 free (s->sigma_r);
395 free (s->sigma_l);
396 free (s->xmid);
397 free (s->x);
398 free (s);
399 GSL_ERROR_VAL ("failed to allocate space for fmin_l", GSL_ENOMEM, 0);
400 }
401
402 s->fmin_r = (double *) malloc (dim * sizeof (double));
403
404 if (s->fmin_r == 0)
405 {
406 free (s->fmin_l);
407 free (s->fmax_r);
408 free (s->fmax_l);
409 free (s->sigma_r);
410 free (s->sigma_l);
411 free (s->xmid);
412 free (s->x);
413 free (s);
414 GSL_ERROR_VAL ("failed to allocate space for fmin_r", GSL_ENOMEM, 0);
415 }
416
417 s->fsum_l = (double *) malloc (dim * sizeof (double));
418
419 if (s->fsum_l == 0)
420 {
421 free (s->fmin_r);
422 free (s->fmin_l);
423 free (s->fmax_r);
424 free (s->fmax_l);
425 free (s->sigma_r);
426 free (s->sigma_l);
427 free (s->xmid);
428 free (s->x);
429 free (s);
430 GSL_ERROR_VAL ("failed to allocate space for fsum_l", GSL_ENOMEM, 0);
431 }
432
433 s->fsum_r = (double *) malloc (dim * sizeof (double));
434
435 if (s->fsum_r == 0)
436 {
437 free (s->fsum_l);
438 free (s->fmin_r);
439 free (s->fmin_l);
440 free (s->fmax_r);
441 free (s->fmax_l);
442 free (s->sigma_r);
443 free (s->sigma_l);
444 free (s->xmid);
445 free (s->x);
446 free (s);
447 GSL_ERROR_VAL ("failed to allocate space for fsum_r", GSL_ENOMEM, 0);
448 }
449
450 s->fsum2_l = (double *) malloc (dim * sizeof (double));
451
452 if (s->fsum2_l == 0)
453 {
454 free (s->fsum_r);
455 free (s->fsum_l);
456 free (s->fmin_r);
457 free (s->fmin_l);
458 free (s->fmax_r);
459 free (s->fmax_l);
460 free (s->sigma_r);
461 free (s->sigma_l);
462 free (s->xmid);
463 free (s->x);
464 free (s);
465 GSL_ERROR_VAL ("failed to allocate space for fsum2_l", GSL_ENOMEM, 0);
466 }
467
468 s->fsum2_r = (double *) malloc (dim * sizeof (double));
469
470 if (s->fsum2_r == 0)
471 {
472 free (s->fsum2_l);
473 free (s->fsum_r);
474 free (s->fsum_l);
475 free (s->fmin_r);
476 free (s->fmin_l);
477 free (s->fmax_r);
478 free (s->fmax_l);
479 free (s->sigma_r);
480 free (s->sigma_l);
481 free (s->xmid);
482 free (s->x);
483 free (s);
484 GSL_ERROR_VAL ("failed to allocate space for fsum2_r", GSL_ENOMEM, 0);
485 }
486
487
488 s->hits_r = (size_t *) malloc (dim * sizeof (size_t));
489
490 if (s->hits_r == 0)
491 {
492 free (s->fsum2_r);
493 free (s->fsum2_l);
494 free (s->fsum_r);
495 free (s->fsum_l);
496 free (s->fmin_r);
497 free (s->fmin_l);
498 free (s->fmax_r);
499 free (s->fmax_l);
500 free (s->sigma_r);
501 free (s->sigma_l);
502 free (s->xmid);
503 free (s->x);
504 free (s);
505 GSL_ERROR_VAL ("failed to allocate space for fsum2_r", GSL_ENOMEM, 0);
506 }
507
508 s->hits_l = (size_t *) malloc (dim * sizeof (size_t));
509
510 if (s->hits_l == 0)
511 {
512 free (s->hits_r);
513 free (s->fsum2_r);
514 free (s->fsum2_l);
515 free (s->fsum_r);
516 free (s->fsum_l);
517 free (s->fmin_r);
518 free (s->fmin_l);
519 free (s->fmax_r);
520 free (s->fmax_l);
521 free (s->sigma_r);
522 free (s->sigma_l);
523 free (s->xmid);
524 free (s->x);
525 free (s);
526 GSL_ERROR_VAL ("failed to allocate space for fsum2_r", GSL_ENOMEM, 0);
527 }
528
529 s->dim = dim;
530
531 gsl_monte_miser_init (s);
532
533 return s;
534 }
535
536 int
gsl_monte_miser_init(gsl_monte_miser_state * s)537 gsl_monte_miser_init (gsl_monte_miser_state * s)
538 {
539 /* We use 8 points in each halfspace to estimate the variance. There are
540 2*dim halfspaces. A variance estimate requires a minimum of 2 points. */
541 s->min_calls = 16 * s->dim;
542 s->min_calls_per_bisection = 32 * s->min_calls;
543 s->estimate_frac = 0.1;
544 s->alpha = 2.0;
545 s->dither = 0.0;
546
547 return GSL_SUCCESS;
548 }
549
550 void
gsl_monte_miser_free(gsl_monte_miser_state * s)551 gsl_monte_miser_free (gsl_monte_miser_state * s)
552 {
553 RETURN_IF_NULL (s);
554 free (s->hits_r);
555 free (s->hits_l);
556 free (s->fsum2_r);
557 free (s->fsum2_l);
558 free (s->fsum_r);
559 free (s->fsum_l);
560 free (s->fmin_r);
561 free (s->fmin_l);
562 free (s->fmax_r);
563 free (s->fmax_l);
564 free (s->sigma_r);
565 free (s->sigma_l);
566 free (s->xmid);
567 free (s->x);
568 free (s);
569 }
570
571 void
gsl_monte_miser_params_get(const gsl_monte_miser_state * s,gsl_monte_miser_params * p)572 gsl_monte_miser_params_get (const gsl_monte_miser_state * s, gsl_monte_miser_params * p)
573 {
574 p->estimate_frac = s->estimate_frac;
575 p->min_calls = s->min_calls;
576 p->min_calls_per_bisection = s->min_calls_per_bisection;
577 p->alpha = s->alpha;
578 p->dither = s->dither;
579 }
580
581 void
gsl_monte_miser_params_set(gsl_monte_miser_state * s,const gsl_monte_miser_params * p)582 gsl_monte_miser_params_set (gsl_monte_miser_state * s, const gsl_monte_miser_params * p)
583 {
584 s->estimate_frac = p->estimate_frac;
585 s->min_calls = p->min_calls;
586 s->min_calls_per_bisection = p->min_calls_per_bisection;
587 s->alpha = p->alpha;
588 s->dither = p->dither;
589 }
590
591 static int
estimate_corrmc(gsl_monte_function * f,const double xl[],const double xu[],size_t dim,size_t calls,gsl_rng * r,gsl_monte_miser_state * state,double * result,double * abserr,const double xmid[],double sigma_l[],double sigma_r[])592 estimate_corrmc (gsl_monte_function * f,
593 const double xl[], const double xu[],
594 size_t dim, size_t calls,
595 gsl_rng * r,
596 gsl_monte_miser_state * state,
597 double *result, double *abserr,
598 const double xmid[], double sigma_l[], double sigma_r[])
599 {
600 size_t i, n;
601
602 double *x = state->x;
603 double *fsum_l = state->fsum_l;
604 double *fsum_r = state->fsum_r;
605 double *fsum2_l = state->fsum2_l;
606 double *fsum2_r = state->fsum2_r;
607 size_t *hits_l = state->hits_l;
608 size_t *hits_r = state->hits_r;
609
610 double m = 0.0, q = 0.0;
611 double vol = 1.0;
612
613 for (i = 0; i < dim; i++)
614 {
615 vol *= xu[i] - xl[i];
616 hits_l[i] = hits_r[i] = 0;
617 fsum_l[i] = fsum_r[i] = 0.0;
618 fsum2_l[i] = fsum2_r[i] = 0.0;
619 sigma_l[i] = sigma_r[i] = -1;
620 }
621
622 for (n = 0; n < calls; n++)
623 {
624 double fval;
625
626 unsigned int j = (n/2) % dim;
627 unsigned int side = (n % 2);
628
629 for (i = 0; i < dim; i++)
630 {
631 double z = gsl_rng_uniform_pos (r) ;
632
633 if (i != j)
634 {
635 x[i] = xl[i] + z * (xu[i] - xl[i]);
636 }
637 else
638 {
639 if (side == 0)
640 {
641 x[i] = xmid[i] + z * (xu[i] - xmid[i]);
642 }
643 else
644 {
645 x[i] = xl[i] + z * (xmid[i] - xl[i]);
646 }
647 }
648 }
649
650 fval = GSL_MONTE_FN_EVAL (f, x);
651
652 /* recurrence for mean and variance */
653 {
654 double d = fval - m;
655 m += d / (n + 1.0);
656 q += d * d * (n / (n + 1.0));
657 }
658
659 /* compute the variances on each side of the bisection */
660 for (i = 0; i < dim; i++)
661 {
662 if (x[i] <= xmid[i])
663 {
664 fsum_l[i] += fval;
665 fsum2_l[i] += fval * fval;
666 hits_l[i]++;
667 }
668 else
669 {
670 fsum_r[i] += fval;
671 fsum2_r[i] += fval * fval;
672 hits_r[i]++;
673 }
674 }
675 }
676
677 for (i = 0; i < dim; i++)
678 {
679 double fraction_l = (xmid[i] - xl[i]) / (xu[i] - xl[i]);
680
681 if (hits_l[i] > 0)
682 {
683 fsum_l[i] /= hits_l[i];
684 sigma_l[i] = sqrt (fsum2_l[i] - fsum_l[i] * fsum_l[i] / hits_l[i]);
685 sigma_l[i] *= fraction_l * vol / hits_l[i];
686 }
687
688 if (hits_r[i] > 0)
689 {
690 fsum_r[i] /= hits_r[i];
691 sigma_r[i] = sqrt (fsum2_r[i] - fsum_r[i] * fsum_r[i] / hits_r[i]);
692 sigma_r[i] *= (1 - fraction_l) * vol / hits_r[i];
693 }
694 }
695
696 *result = vol * m;
697
698 if (calls < 2)
699 {
700 *abserr = GSL_POSINF;
701 }
702 else
703 {
704 *abserr = vol * sqrt (q / (calls * (calls - 1.0)));
705 }
706
707 return GSL_SUCCESS;
708 }
709
710