htswanalysis/MACS/lib/gsl/gsl-1.11/monte/vegas.c

   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 <config.h>
  55
  56 /* standard headers */
  57 #include <math.h>
  58 #include <stdio.h>
  59
  60 /* gsl headers */
  61 #include <gsl/gsl_math.h>
  62 #include <gsl/gsl_errno.h>
  63 #include <gsl/gsl_rng.h>
  64 #include <gsl/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
 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 *
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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