1 /* monte/vegas.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Michael Booth
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 /* Author: MJB */
21 /* Modified by: Brian Gough, 12/2000 */
22
23 /* This is an implementation of the adaptive Monte-Carlo algorithm
24 of G. P. Lepage, originally described in J. Comp. Phys. 27, 192(1978).
25 The current version of the algorithm was described in the Cornell
26 preprint CLNS-80/447 of March, 1980.
27
28 This code follows most closely the c version by D.R.Yennie, coded
29 in 1984.
30
31 The input coordinates are x[j], with upper and lower limits xu[j]
32 and xl[j]. The integration length in the j-th direction is
33 delx[j]. Each coordinate x[j] is rescaled to a variable y[j] in
34 the range 0 to 1. The range is divided into bins with boundaries
35 xi[i][j], where i=0 corresponds to y=0 and i=bins to y=1. The grid
36 is refined (ie, bins are adjusted) using d[i][j] which is some
37 variation on the squared sum. A third parameter used in defining
38 the real coordinate using random numbers is called z. It ranges
39 from 0 to bins. Its integer part gives the lower index of the bin
40 into which a call is to be placed, and the remainder gives the
41 location inside the bin.
42
43 When stratified sampling is used the bins are grouped into boxes,
44 and the algorithm allocates an equal number of function calls to
45 each box.
46
47 The variable alpha controls how "stiff" the rebinning algorithm is.
48 alpha = 0 means never change the grid. Alpha is typically set between
49 1 and 2.
50
51 */
52
53 /* configuration headers */
54 #include "gsl__config.h"
55
56 /* standard headers */
57 #include <math.h>
58 #include <stdio.h>
59
60 /* gsl headers */
61 #include "gsl_math.h"
62 #include "gsl_errno.h"
63 #include "gsl_rng.h"
64 #include "gsl_monte_vegas.h"
65
66 /* lib-specific headers */
67 #define BINS_MAX 50 /* even integer, will be divided by two */
68
69 /* A separable grid with coordinates and values */
70 #define COORD(s,i,j) ((s)->xi[(i)*(s)->dim + (j)])
71 #define NEW_COORD(s,i) ((s)->xin[(i)])
72 #define VALUE(s,i,j) ((s)->d[(i)*(s)->dim + (j)])
73
74 /* predeclare functions */
75
76 typedef int coord;
77
78 static void init_grid (gsl_monte_vegas_state * s, double xl[], double xu[],
79 size_t dim);
80 static void reset_grid_values (gsl_monte_vegas_state * s);
81 static void init_box_coord (gsl_monte_vegas_state * s, coord box[]);
82 static int change_box_coord (gsl_monte_vegas_state * s, coord box[]);
83 static void accumulate_distribution (gsl_monte_vegas_state * s, coord bin[],
84 double y);
85 static void random_point (double x[], coord bin[], double *bin_vol,
86 const coord box[],
87 const double xl[], const double xu[],
88 gsl_monte_vegas_state * s, gsl_rng * r);
89 static void resize_grid (gsl_monte_vegas_state * s, unsigned int bins);
90 static void refine_grid (gsl_monte_vegas_state * s);
91
92 static void print_lim (gsl_monte_vegas_state * state,
93 double xl[], double xu[], unsigned long dim);
94 static void print_head (gsl_monte_vegas_state * state,
95 unsigned long num_dim, unsigned long calls,
96 unsigned int it_num,
97 unsigned int bins, unsigned int boxes);
98 static void print_res (gsl_monte_vegas_state * state,
99 unsigned int itr, double res, double err,
100 double cum_res, double cum_err,
101 double chi_sq);
102 static void print_dist (gsl_monte_vegas_state * state, unsigned long dim);
103 static void print_grid (gsl_monte_vegas_state * state, unsigned long dim);
104
105 int
gsl_monte_vegas_integrate(gsl_monte_function * f,double xl[],double xu[],size_t dim,size_t calls,gsl_rng * r,gsl_monte_vegas_state * state,double * result,double * abserr)106 gsl_monte_vegas_integrate (gsl_monte_function * f,
107 double xl[], double xu[],
108 size_t dim, size_t calls,
109 gsl_rng * r,
110 gsl_monte_vegas_state * state,
111 double *result, double *abserr)
112 {
113 double cum_int, cum_sig;
114 size_t i, k, it;
115
116 if (dim != state->dim)
117 {
118 GSL_ERROR ("number of dimensions must match allocated size", GSL_EINVAL);
119 }
120
121 for (i = 0; i < dim; i++)
122 {
123 if (xu[i] <= xl[i])
124 {
125 GSL_ERROR ("xu must be greater than xl", GSL_EINVAL);
126 }
127
128 if (xu[i] - xl[i] > GSL_DBL_MAX)
129 {
130 GSL_ERROR ("Range of integration is too large, please rescale",
131 GSL_EINVAL);
132 }
133 }
134
135 if (state->stage == 0)
136 {
137 init_grid (state, xl, xu, dim);
138
139 if (state->verbose >= 0)
140 {
141 print_lim (state, xl, xu, dim);
142 }
143 }
144
145 if (state->stage <= 1)
146 {
147 state->wtd_int_sum = 0;
148 state->sum_wgts = 0;
149 state->chi_sum = 0;
150 state->it_num = 1;
151 state->samples = 0;
152 }
153
154 if (state->stage <= 2)
155 {
156 unsigned int bins = state->bins_max;
157 unsigned int boxes = 1;
158
159 if (state->mode != GSL_VEGAS_MODE_IMPORTANCE_ONLY)
160 {
161 /* shooting for 2 calls/box */
162
163 boxes = floor (pow (calls / 2.0, 1.0 / dim));
164 state->mode = GSL_VEGAS_MODE_IMPORTANCE;
165
166 if (2 * boxes >= state->bins_max)
167 {
168 /* if bins/box < 2 */
169 int box_per_bin = GSL_MAX (boxes / state->bins_max, 1);
170
171 bins = GSL_MIN(boxes / box_per_bin, state->bins_max);
172 boxes = box_per_bin * bins;
173
174 state->mode = GSL_VEGAS_MODE_STRATIFIED;
175 }
176 }
177
178 {
179 double tot_boxes = pow ((double) boxes, (double) dim);
180 state->calls_per_box = GSL_MAX (calls / tot_boxes, 2);
181 calls = state->calls_per_box * tot_boxes;
182 }
183
184 /* total volume of x-space/(avg num of calls/bin) */
185 state->jac = state->vol * pow ((double) bins, (double) dim) / calls;
186
187 state->boxes = boxes;
188
189 /* If the number of bins changes from the previous invocation, bins
190 are expanded or contracted accordingly, while preserving bin
191 density */
192
193 if (bins != state->bins)
194 {
195 resize_grid (state, bins);
196
197 if (state->verbose > 1)
198 {
199 print_grid (state, dim);
200 }
201 }
202
203 if (state->verbose >= 0)
204 {
205 print_head (state,
206 dim, calls, state->it_num, state->bins, state->boxes);
207 }
208 }
209
210 state->it_start = state->it_num;
211
212 cum_int = 0.0;
213 cum_sig = 0.0;
214
215 state->chisq = 0.0;
216
217 for (it = 0; it < state->iterations; it++)
218 {
219 double intgrl = 0.0, intgrl_sq = 0.0;
220 double sig = 0.0;
221 double wgt;
222 size_t calls_per_box = state->calls_per_box;
223 double jacbin = state->jac;
224 double *x = state->x;
225 coord *bin = state->bin;
226
227 state->it_num = state->it_start + it;
228
229 reset_grid_values (state);
230 init_box_coord (state, state->box);
231
232 do
233 {
234 double m = 0, q = 0;
235 double f_sq_sum = 0.0;
236
237 for (k = 0; k < calls_per_box; k++)
238 {
239 double fval, bin_vol;
240
241 random_point (x, bin, &bin_vol, state->box, xl, xu, state, r);
242
243 fval = jacbin * bin_vol * GSL_MONTE_FN_EVAL (f, x);
244
245 /* recurrence for mean and variance */
246
247 {
248 double d = fval - m;
249 m += d / (k + 1.0);
250 q += d * d * (k / (k + 1.0));
251 }
252
253 if (state->mode != GSL_VEGAS_MODE_STRATIFIED)
254 {
255 double f_sq = fval * fval;
256 accumulate_distribution (state, bin, f_sq);
257 }
258 }
259
260 intgrl += m * calls_per_box;
261
262 f_sq_sum = q * calls_per_box ;
263
264 sig += f_sq_sum ;
265
266 if (state->mode == GSL_VEGAS_MODE_STRATIFIED)
267 {
268 accumulate_distribution (state, bin, f_sq_sum);
269 }
270 }
271 while (change_box_coord (state, state->box));
272
273 /* Compute final results for this iteration */
274
275 sig = sig / (calls_per_box - 1.0) ;
276
277 if (sig > 0)
278 {
279 wgt = 1.0 / sig;
280 }
281 else if (state->sum_wgts > 0)
282 {
283 wgt = state->sum_wgts / state->samples;
284 }
285 else
286 {
287 wgt = 0.0;
288 }
289
290 intgrl_sq = intgrl * intgrl;
291
292 state->result = intgrl;
293 state->sigma = sqrt(sig);
294
295 if (wgt > 0.0)
296 {
297 state->samples++ ;
298 state->sum_wgts += wgt;
299 state->wtd_int_sum += intgrl * wgt;
300 state->chi_sum += intgrl_sq * wgt;
301
302 cum_int = state->wtd_int_sum / state->sum_wgts;
303 cum_sig = sqrt (1 / state->sum_wgts);
304
305 if (state->samples > 1)
306 {
307 state->chisq = (state->chi_sum - state->wtd_int_sum * cum_int) /
308 (state->samples - 1.0);
309 }
310 }
311 else
312 {
313 cum_int += (intgrl - cum_int) / (it + 1.0);
314 cum_sig = 0.0;
315 }
316
317
318 if (state->verbose >= 0)
319 {
320 print_res (state,
321 state->it_num, intgrl, sqrt (sig), cum_int, cum_sig,
322 state->chisq);
323 if (it + 1 == state->iterations && state->verbose > 0)
324 {
325 print_grid (state, dim);
326 }
327 }
328
329 if (state->verbose > 1)
330 {
331 print_dist (state, dim);
332 }
333
334 refine_grid (state);
335
336 if (state->verbose > 1)
337 {
338 print_grid (state, dim);
339 }
340
341 }
342
343 /* By setting stage to 1 further calls will generate independent
344 estimates based on the same grid, although it may be rebinned. */
345
346 state->stage = 1;
347
348 *result = cum_int;
349 *abserr = cum_sig;
350
351 return GSL_SUCCESS;
352 }
353
354
355
356 gsl_monte_vegas_state *
gsl_monte_vegas_alloc(size_t dim)357 gsl_monte_vegas_alloc (size_t dim)
358 {
359 gsl_monte_vegas_state *s =
360 (gsl_monte_vegas_state *) malloc (sizeof (gsl_monte_vegas_state));
361
362 if (s == 0)
363 {
364 GSL_ERROR_VAL ("failed to allocate space for vegas state struct",
365 GSL_ENOMEM, 0);
366 }
367
368 s->delx = (double *) malloc (dim * sizeof (double));
369
370 if (s->delx == 0)
371 {
372 free (s);
373 GSL_ERROR_VAL ("failed to allocate space for delx", GSL_ENOMEM, 0);
374 }
375
376 s->d = (double *) malloc (BINS_MAX * dim * sizeof (double));
377
378 if (s->d == 0)
379 {
380 free (s->delx);
381 free (s);
382 GSL_ERROR_VAL ("failed to allocate space for d", GSL_ENOMEM, 0);
383 }
384
385 s->xi = (double *) malloc ((BINS_MAX + 1) * dim * sizeof (double));
386
387 if (s->xi == 0)
388 {
389 free (s->d);
390 free (s->delx);
391 free (s);
392 GSL_ERROR_VAL ("failed to allocate space for xi", GSL_ENOMEM, 0);
393 }
394
395 s->xin = (double *) malloc ((BINS_MAX + 1) * sizeof (double));
396
397 if (s->xin == 0)
398 {
399 free (s->xi);
400 free (s->d);
401 free (s->delx);
402 free (s);
403 GSL_ERROR_VAL ("failed to allocate space for xin", GSL_ENOMEM, 0);
404 }
405
406 s->weight = (double *) malloc (BINS_MAX * sizeof (double));
407
408 if (s->weight == 0)
409 {
410 free (s->xin);
411 free (s->xi);
412 free (s->d);
413 free (s->delx);
414 free (s);
415 GSL_ERROR_VAL ("failed to allocate space for xin", GSL_ENOMEM, 0);
416 }
417
418 s->box = (coord *) malloc (dim * sizeof (coord));
419
420 if (s->box == 0)
421 {
422 free (s->weight);
423 free (s->xin);
424 free (s->xi);
425 free (s->d);
426 free (s->delx);
427 free (s);
428 GSL_ERROR_VAL ("failed to allocate space for box", GSL_ENOMEM, 0);
429 }
430
431 s->bin = (coord *) malloc (dim * sizeof (coord));
432
433 if (s->bin == 0)
434 {
435 free (s->box);
436 free (s->weight);
437 free (s->xin);
438 free (s->xi);
439 free (s->d);
440 free (s->delx);
441 free (s);
442 GSL_ERROR_VAL ("failed to allocate space for bin", GSL_ENOMEM, 0);
443 }
444
445 s->x = (double *) malloc (dim * sizeof (double));
446
447 if (s->x == 0)
448 {
449 free (s->bin);
450 free (s->box);
451 free (s->weight);
452 free (s->xin);
453 free (s->xi);
454 free (s->d);
455 free (s->delx);
456 free (s);
457 GSL_ERROR_VAL ("failed to allocate space for x", GSL_ENOMEM, 0);
458 }
459
460 s->dim = dim;
461 s->bins_max = BINS_MAX;
462
463 gsl_monte_vegas_init (s);
464
465 return s;
466 }
467
468 /* Set some default values and whatever */
469 int
gsl_monte_vegas_init(gsl_monte_vegas_state * state)470 gsl_monte_vegas_init (gsl_monte_vegas_state * state)
471 {
472 state->stage = 0;
473 state->alpha = 1.5;
474 state->verbose = -1;
475 state->iterations = 5;
476 state->mode = GSL_VEGAS_MODE_IMPORTANCE;
477 state->chisq = 0;
478 state->bins = state->bins_max;
479 state->ostream = stdout;
480
481 return GSL_SUCCESS;
482 }
483
484 void
gsl_monte_vegas_free(gsl_monte_vegas_state * s)485 gsl_monte_vegas_free (gsl_monte_vegas_state * s)
486 {
487 free (s->x);
488 free (s->delx);
489 free (s->d);
490 free (s->xi);
491 free (s->xin);
492 free (s->weight);
493 free (s->box);
494 free (s->bin);
495 free (s);
496 }
497
498 static void
init_box_coord(gsl_monte_vegas_state * s,coord box[])499 init_box_coord (gsl_monte_vegas_state * s, coord box[])
500 {
501 size_t i;
502
503 size_t dim = s->dim;
504
505 for (i = 0; i < dim; i++)
506 {
507 box[i] = 0;
508 }
509 }
510
511 /* change_box_coord steps through the box coord like
512 {0,0}, {0, 1}, {0, 2}, {0, 3}, {1, 0}, {1, 1}, {1, 2}, ...
513 */
514 static int
change_box_coord(gsl_monte_vegas_state * s,coord box[])515 change_box_coord (gsl_monte_vegas_state * s, coord box[])
516 {
517 int j = s->dim - 1;
518
519 int ng = s->boxes;
520
521 while (j >= 0)
522 {
523 box[j] = (box[j] + 1) % ng;
524
525 if (box[j] != 0)
526 {
527 return 1;
528 }
529
530 j--;
531 }
532
533 return 0;
534 }
535
536 static void
init_grid(gsl_monte_vegas_state * s,double xl[],double xu[],size_t dim)537 init_grid (gsl_monte_vegas_state * s, double xl[], double xu[], size_t dim)
538 {
539 size_t j;
540
541 double vol = 1.0;
542
543 s->bins = 1;
544
545 for (j = 0; j < dim; j++)
546 {
547 double dx = xu[j] - xl[j];
548 s->delx[j] = dx;
549 vol *= dx;
550
551 COORD (s, 0, j) = 0.0;
552 COORD (s, 1, j) = 1.0;
553 }
554
555 s->vol = vol;
556 }
557
558
559 static void
reset_grid_values(gsl_monte_vegas_state * s)560 reset_grid_values (gsl_monte_vegas_state * s)
561 {
562 size_t i, j;
563
564 size_t dim = s->dim;
565 size_t bins = s->bins;
566
567 for (i = 0; i < bins; i++)
568 {
569 for (j = 0; j < dim; j++)
570 {
571 VALUE (s, i, j) = 0.0;
572 }
573 }
574 }
575
576 static void
accumulate_distribution(gsl_monte_vegas_state * s,coord bin[],double y)577 accumulate_distribution (gsl_monte_vegas_state * s, coord bin[], double y)
578 {
579 size_t j;
580 size_t dim = s->dim;
581
582 for (j = 0; j < dim; j++)
583 {
584 int i = bin[j];
585 VALUE (s, i, j) += y;
586 }
587 }
588
589 static void
random_point(double x[],coord bin[],double * bin_vol,const coord box[],const double xl[],const double xu[],gsl_monte_vegas_state * s,gsl_rng * r)590 random_point (double x[], coord bin[], double *bin_vol,
591 const coord box[], const double xl[], const double xu[],
592 gsl_monte_vegas_state * s, gsl_rng * r)
593 {
594 /* Use the random number generator r to return a random position x
595 in a given box. The value of bin gives the bin location of the
596 random position (there may be several bins within a given box) */
597
598 double vol = 1.0;
599
600 size_t j;
601
602 size_t dim = s->dim;
603 size_t bins = s->bins;
604 size_t boxes = s->boxes;
605
606 DISCARD_POINTER(xu); /* prevent warning about unused parameter */
607
608 for (j = 0; j < dim; ++j)
609 {
610 /* box[j] + ran gives the position in the box units, while z
611 is the position in bin units. */
612
613 double z = ((box[j] + gsl_rng_uniform_pos (r)) / boxes) * bins;
614
615 int k = z;
616
617 double y, bin_width;
618
619 bin[j] = k;
620
621 if (k == 0)
622 {
623 bin_width = COORD (s, 1, j);
624 y = z * bin_width;
625 }
626 else
627 {
628 bin_width = COORD (s, k + 1, j) - COORD (s, k, j);
629 y = COORD (s, k, j) + (z - k) * bin_width;
630 }
631
632 x[j] = xl[j] + y * s->delx[j];
633
634 vol *= bin_width;
635 }
636
637 *bin_vol = vol;
638 }
639
640
641 static void
resize_grid(gsl_monte_vegas_state * s,unsigned int bins)642 resize_grid (gsl_monte_vegas_state * s, unsigned int bins)
643 {
644 size_t j, k;
645 size_t dim = s->dim;
646
647 /* weight is ratio of bin sizes */
648
649 double pts_per_bin = (double) s->bins / (double) bins;
650
651 for (j = 0; j < dim; j++)
652 {
653 double xold;
654 double xnew = 0;
655 double dw = 0;
656 int i = 1;
657
658 for (k = 1; k <= s->bins; k++)
659 {
660 dw += 1.0;
661 xold = xnew;
662 xnew = COORD (s, k, j);
663
664 for (; dw > pts_per_bin; i++)
665 {
666 dw -= pts_per_bin;
667 NEW_COORD (s, i) = xnew - (xnew - xold) * dw;
668 }
669 }
670
671 for (k = 1 ; k < bins; k++)
672 {
673 COORD(s, k, j) = NEW_COORD(s, k);
674 }
675
676 COORD (s, bins, j) = 1;
677 }
678
679 s->bins = bins;
680 }
681
682 static void
refine_grid(gsl_monte_vegas_state * s)683 refine_grid (gsl_monte_vegas_state * s)
684 {
685 size_t i, j, k;
686 size_t dim = s->dim;
687 size_t bins = s->bins;
688
689 for (j = 0; j < dim; j++)
690 {
691 double grid_tot_j, tot_weight;
692 double * weight = s->weight;
693
694 double oldg = VALUE (s, 0, j);
695 double newg = VALUE (s, 1, j);
696
697 VALUE (s, 0, j) = (oldg + newg) / 2;
698 grid_tot_j = VALUE (s, 0, j);
699
700 /* This implements gs[i][j] = (gs[i-1][j]+gs[i][j]+gs[i+1][j])/3 */
701
702 for (i = 1; i < bins - 1; i++)
703 {
704 double rc = oldg + newg;
705 oldg = newg;
706 newg = VALUE (s, i + 1, j);
707 VALUE (s, i, j) = (rc + newg) / 3;
708 grid_tot_j += VALUE (s, i, j);
709 }
710 VALUE (s, bins - 1, j) = (newg + oldg) / 2;
711
712 grid_tot_j += VALUE (s, bins - 1, j);
713
714 tot_weight = 0;
715
716 for (i = 0; i < bins; i++)
717 {
718 weight[i] = 0;
719
720 if (VALUE (s, i, j) > 0)
721 {
722 oldg = grid_tot_j / VALUE (s, i, j);
723 /* damped change */
724 weight[i] = pow (((oldg - 1) / oldg / log (oldg)), s->alpha);
725 }
726
727 tot_weight += weight[i];
728
729 #ifdef DEBUG
730 printf("weight[%d] = %g\n", i, weight[i]);
731 #endif
732 }
733
734 {
735 double pts_per_bin = tot_weight / bins;
736
737 double xold;
738 double xnew = 0;
739 double dw = 0;
740 i = 1;
741
742 for (k = 0; k < bins; k++)
743 {
744 dw += weight[k];
745 xold = xnew;
746 xnew = COORD (s, k + 1, j);
747
748 for (; dw > pts_per_bin; i++)
749 {
750 dw -= pts_per_bin;
751 NEW_COORD (s, i) = xnew - (xnew - xold) * dw / weight[k];
752 }
753 }
754
755 for (k = 1 ; k < bins ; k++)
756 {
757 COORD(s, k, j) = NEW_COORD(s, k);
758 }
759
760 COORD (s, bins, j) = 1;
761 }
762 }
763 }
764
765
766 static void
print_lim(gsl_monte_vegas_state * state,double xl[],double xu[],unsigned long dim)767 print_lim (gsl_monte_vegas_state * state,
768 double xl[], double xu[], unsigned long dim)
769 {
770 unsigned long j;
771
772 fprintf (state->ostream, "The limits of integration are:\n");
773 for (j = 0; j < dim; ++j)
774 fprintf (state->ostream, "\nxl[%lu]=%f xu[%lu]=%f", j, xl[j], j, xu[j]);
775 fprintf (state->ostream, "\n");
776 fflush (state->ostream);
777 }
778
779 static void
print_head(gsl_monte_vegas_state * state,unsigned long num_dim,unsigned long calls,unsigned int it_num,unsigned int bins,unsigned int boxes)780 print_head (gsl_monte_vegas_state * state,
781 unsigned long num_dim, unsigned long calls,
782 unsigned int it_num, unsigned int bins, unsigned int boxes)
783 {
784 fprintf (state->ostream,
785 "\nnum_dim=%lu, calls=%lu, it_num=%d, max_it_num=%d ",
786 num_dim, calls, it_num, state->iterations);
787 fprintf (state->ostream,
788 "verb=%d, alph=%.2f,\nmode=%d, bins=%d, boxes=%d\n",
789 state->verbose, state->alpha, state->mode, bins, boxes);
790 fprintf (state->ostream,
791 "\n single.......iteration ");
792 fprintf (state->ostream, "accumulated......results \n");
793
794 fprintf (state->ostream,
795 "iteration integral sigma integral ");
796 fprintf (state->ostream, " sigma chi-sq/it\n\n");
797 fflush (state->ostream);
798
799 }
800
801 static void
print_res(gsl_monte_vegas_state * state,unsigned int itr,double res,double err,double cum_res,double cum_err,double chi_sq)802 print_res (gsl_monte_vegas_state * state,
803 unsigned int itr,
804 double res, double err,
805 double cum_res, double cum_err,
806 double chi_sq)
807 {
808 fprintf (state->ostream,
809 "%4d %6.4e %10.2e %6.4e %8.2e %10.2e\n",
810 itr, res, err, cum_res, cum_err, chi_sq);
811 fflush (state->ostream);
812 }
813
814 static void
print_dist(gsl_monte_vegas_state * state,unsigned long dim)815 print_dist (gsl_monte_vegas_state * state, unsigned long dim)
816 {
817 unsigned long i, j;
818 int p = state->verbose;
819 if (p < 1)
820 return;
821
822 for (j = 0; j < dim; ++j)
823 {
824 fprintf (state->ostream, "\n axis %lu \n", j);
825 fprintf (state->ostream, " x g\n");
826 for (i = 0; i < state->bins; i++)
827 {
828 fprintf (state->ostream, "weight [%11.2e , %11.2e] = ",
829 COORD (state, i, j), COORD(state,i+1,j));
830 fprintf (state->ostream, " %11.2e\n", VALUE (state, i, j));
831
832 }
833 fprintf (state->ostream, "\n");
834 }
835 fprintf (state->ostream, "\n");
836 fflush (state->ostream);
837
838 }
839
840 static void
print_grid(gsl_monte_vegas_state * state,unsigned long dim)841 print_grid (gsl_monte_vegas_state * state, unsigned long dim)
842 {
843 unsigned long i, j;
844 int p = state->verbose;
845 if (p < 1)
846 return;
847
848 for (j = 0; j < dim; ++j)
849 {
850 fprintf (state->ostream, "\n axis %lu \n", j);
851 fprintf (state->ostream, " x \n");
852 for (i = 0; i <= state->bins; i++)
853 {
854 fprintf (state->ostream, "%11.2e", COORD (state, i, j));
855 if (i % 5 == 4)
856 fprintf (state->ostream, "\n");
857 }
858 fprintf (state->ostream, "\n");
859 }
860 fprintf (state->ostream, "\n");
861 fflush (state->ostream);
862
863 }
864
865