vegas.c

该文件为c++的数学函数库!是一个非常有用的编程工具.它含有各种数学函数,为科学计算、工程应用等程序编写提供方便！

/* monte/vegas.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Michael Booth
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author: MJB */
/* Modified by: Brian Gough, 12/2000 */

/* This is an implementation of the adaptive Monte-Carlo algorithm
   of G. P. Lepage, originally described in J. Comp. Phys. 27, 192(1978).
   The current version of the algorithm was described in the Cornell
   preprint CLNS-80/447 of March, 1980.

   This code follows most closely the c version by D.R.Yennie, coded
   in 1984.

   The input coordinates are x[j], with upper and lower limits xu[j]
   and xl[j].  The integration length in the j-th direction is
   delx[j].  Each coordinate x[j] is rescaled to a variable y[j] in
   the range 0 to 1.  The range is divided into bins with boundaries
   xi[i][j], where i=0 corresponds to y=0 and i=bins to y=1.  The grid
   is refined (ie, bins are adjusted) using d[i][j] which is some
   variation on the squared sum.  A third parameter used in defining
   the real coordinate using random numbers is called z.  It ranges
   from 0 to bins.  Its integer part gives the lower index of the bin
   into which a call is to be placed, and the remainder gives the
   location inside the bin.

   When stratified sampling is used the bins are grouped into boxes,
   and the algorithm allocates an equal number of function calls to
   each box.

   The variable alpha controls how "stiff" the rebinning algorithm is.  
   alpha = 0 means never change the grid.  Alpha is typically set between
   1 and 2.

   */

/* configuration headers */
#include <config.h>

/* standard headers */
#include <math.h>
#include <stdio.h>

/* gsl headers */
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_monte_vegas.h>

/* lib-specific headers */
#define BINS_MAX 50             /* even integer, will be divided by two */

/* A separable grid with coordinates and values */
#define COORD(s,i,j) ((s)->xi[(i)*(s)->dim + (j)])
#define NEW_COORD(s,i) ((s)->xin[(i)])
#define VALUE(s,i,j) ((s)->d[(i)*(s)->dim + (j)])

/* predeclare functions */

typedef int coord;

static void init_grid (gsl_monte_vegas_state * s, double xl[], double xu[],
                size_t dim);
static void reset_grid_values (gsl_monte_vegas_state * s);
static void init_box_coord (gsl_monte_vegas_state * s, coord box[]);
static int change_box_coord (gsl_monte_vegas_state * s, coord box[]);
static void accumulate_distribution (gsl_monte_vegas_state * s, coord bin[],
                                     double y);
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);
static void resize_grid (gsl_monte_vegas_state * s, unsigned int bins);
static void refine_grid (gsl_monte_vegas_state * s);

static void print_lim (gsl_monte_vegas_state * state,
                       double xl[], double xu[], unsigned long dim);
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);
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);
static void print_dist (gsl_monte_vegas_state * state, unsigned long dim);
static void print_grid (gsl_monte_vegas_state * state, unsigned long dim);

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)
{
  double cum_int, cum_sig;
  size_t i, k, it;

  if (dim != state->dim)
    {
      GSL_ERROR ("number of dimensions must match allocated size", GSL_EINVAL);
    }

  for (i = 0; i < dim; i++)
    {
      if (xu[i] <= xl[i])
        {
          GSL_ERROR ("xu must be greater than xl", GSL_EINVAL);
        }

      if (xu[i] - xl[i] > GSL_DBL_MAX)
        {
          GSL_ERROR ("Range of integration is too large, please rescale",
                     GSL_EINVAL);
        }
    }

  if (state->stage == 0)
    {
      init_grid (state, xl, xu, dim);

      if (state->verbose >= 0)
        {
          print_lim (state, xl, xu, dim);
        }
    }

  if (state->stage <= 1)
    {
      state->wtd_int_sum = 0;
      state->sum_wgts = 0;
      state->chi_sum = 0;
      state->it_num = 1;
      state->samples = 0;
    }

  if (state->stage <= 2)
    {
      unsigned int bins = state->bins_max;
      unsigned int boxes = 1;

      if (state->mode != GSL_VEGAS_MODE_IMPORTANCE_ONLY)
        {
          /* shooting for 2 calls/box */

          boxes = floor (pow (calls / 2.0, 1.0 / dim));
          state->mode = GSL_VEGAS_MODE_IMPORTANCE;

          if (2 * boxes >= state->bins_max)
            {
              /* if bins/box < 2 */
              int box_per_bin = GSL_MAX (boxes / state->bins_max, 1);

              bins = GSL_MIN(boxes / box_per_bin, state->bins_max);
              boxes = box_per_bin * bins;

              state->mode = GSL_VEGAS_MODE_STRATIFIED;
            }
        }

      {
        double tot_boxes = pow ((double) boxes, (double) dim);
        state->calls_per_box = GSL_MAX (calls / tot_boxes, 2);
        calls = state->calls_per_box * tot_boxes;
      }

      /* total volume of x-space/(avg num of calls/bin) */
      state->jac = state->vol * pow ((double) bins, (double) dim) / calls;

      state->boxes = boxes;

      /* If the number of bins changes from the previous invocation, bins
         are expanded or contracted accordingly, while preserving bin
         density */

      if (bins != state->bins)
        {
          resize_grid (state, bins);

          if (state->verbose > 1)
            {
              print_grid (state, dim);
            }
        }

      if (state->verbose >= 0)
        {
          print_head (state,
                      dim, calls, state->it_num, state->bins, state->boxes);
        }
    }

  state->it_start = state->it_num;

  cum_int = 0.0;
  cum_sig = 0.0;

  state->chisq = 0.0;

  for (it = 0; it < state->iterations; it++)
    {
      double intgrl = 0.0, intgrl_sq = 0.0;
      double sig = 0.0;
      double wgt;
      size_t calls_per_box = state->calls_per_box;
      double jacbin = state->jac;
      double *x = state->x;
      coord *bin = state->bin;

      state->it_num = state->it_start + it;

      reset_grid_values (state);
      init_box_coord (state, state->box);
      
      do
        {
          double m = 0, q = 0;
          double f_sq_sum = 0.0;

          for (k = 0; k < calls_per_box; k++)
            {
              double fval, bin_vol;

              random_point (x, bin, &bin_vol, state->box, xl, xu, state, r);

              fval = jacbin * bin_vol * GSL_MONTE_FN_EVAL (f, x);

              /* recurrence for mean and variance */

              {
                double d = fval - m;
                m += d / (k + 1.0);
                q += d * d * (k / (k + 1.0));
              }

              if (state->mode != GSL_VEGAS_MODE_STRATIFIED)
                {
                  double f_sq = fval * fval;
                  accumulate_distribution (state, bin, f_sq);
                }
            }

          intgrl += m * calls_per_box;

          f_sq_sum = q * calls_per_box ;

          sig += f_sq_sum ;

          if (state->mode == GSL_VEGAS_MODE_STRATIFIED)
            {
              accumulate_distribution (state, bin, f_sq_sum);
            }
        }
      while (change_box_coord (state, state->box));

      /* Compute final results for this iteration   */

      sig = sig / (calls_per_box - 1.0)  ;

      if (sig > 0) 
        {
          wgt = 1.0 / sig;
        }
      else if (state->sum_wgts > 0) 
        {
          wgt = state->sum_wgts / state->samples;
        }
      else 
        {
          wgt = 0.0;
        }
        
     intgrl_sq = intgrl * intgrl;

     state->result = intgrl;
     state->sigma  = sqrt(sig);

     if (wgt > 0.0)
       {
         state->samples++ ;
         state->sum_wgts += wgt;
         state->wtd_int_sum += intgrl * wgt;
         state->chi_sum += intgrl_sq * wgt;

         cum_int = state->wtd_int_sum / state->sum_wgts;
         cum_sig = sqrt (1 / state->sum_wgts);

         if (state->samples > 1)
           {
             state->chisq = (state->chi_sum - state->wtd_int_sum * cum_int) /
               (state->samples - 1.0);
           }
       }
     else
       {
         cum_int += (intgrl - cum_int) / (it + 1.0);
         cum_sig = 0.0;
       }         


      if (state->verbose >= 0)
        {
          print_res (state,
                     state->it_num, intgrl, sqrt (sig), cum_int, cum_sig,
                     state->chisq);
          if (it + 1 == state->iterations && state->verbose > 0)
            {
              print_grid (state, dim);
            }
        }

      if (state->verbose > 1)
        {
          print_dist (state, dim);
        }

      refine_grid (state);

      if (state->verbose > 1)
        {
          print_grid (state, dim);
        }

    }

  /* By setting stage to 1 further calls will generate independent
     estimates based on the same grid, although it may be rebinned. */

  state->stage = 1;  

  *result = cum_int;
  *abserr = cum_sig;

  return GSL_SUCCESS;
}



gsl_monte_vegas_state *
gsl_monte_vegas_alloc (size_t dim)
{
  gsl_monte_vegas_state *s =
    (gsl_monte_vegas_state *) malloc (sizeof (gsl_monte_vegas_state));

  if (s == 0)
    {
      GSL_ERROR_VAL ("failed to allocate space for vegas state struct",
                     GSL_ENOMEM, 0);
    }

  s->delx = (double *) malloc (dim * sizeof (double));

  if (s->delx == 0)
    {
      free (s);
      GSL_ERROR_VAL ("failed to allocate space for delx", GSL_ENOMEM, 0);
    }

  s->d = (double *) malloc (BINS_MAX * dim * sizeof (double));

  if (s->d == 0)
    {
      free (s->delx);
      free (s);
      GSL_ERROR_VAL ("failed to allocate space for d", GSL_ENOMEM, 0);
    }

  s->xi = (double *) malloc ((BINS_MAX + 1) * dim * sizeof (double));

  if (s->xi == 0)
    {
      free (s->d);
      free (s->delx);
      free (s);
      GSL_ERROR_VAL ("failed to allocate space for xi", GSL_ENOMEM, 0);
    }

  s->xin = (double *) malloc ((BINS_MAX + 1) * sizeof (double));

  if (s->xin == 0)
    {
      free (s->xi);
      free (s->d);
      free (s->delx);
      free (s);
      GSL_ERROR_VAL ("failed to allocate space for xin", GSL_ENOMEM, 0);
    }

  s->weight = (double *) malloc (BINS_MAX * sizeof (double));

  if (s->weight == 0)
    {
      free (s->xin);
      free (s->xi);
      free (s->d);
      free (s->delx);
      free (s);
      GSL_ERROR_VAL ("failed to allocate space for xin", GSL_ENOMEM, 0);
    }

  s->box = (coord *) malloc (dim * sizeof (coord));

  if (s->box == 0)
    {
      free (s->weight);
      free (s->xin);
      free (s->xi);
      free (s->d);
      free (s->delx);
      free (s);
      GSL_ERROR_VAL ("failed to allocate space for box", GSL_ENOMEM, 0);
    }

  s->bin = (coord *) malloc (dim * sizeof (coord));

  if (s->bin == 0)