vegas.c

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

    {
      free (s->box);
      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 bin", GSL_ENOMEM, 0);
    }

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

  if (s->x == 0)
    {
      free (s->bin);
      free (s->box);
      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 x", GSL_ENOMEM, 0);
    }

  s->dim = dim;
  s->bins_max = BINS_MAX;

  gsl_monte_vegas_init (s);

  return s;
}

/* Set some default values and whatever */
int
gsl_monte_vegas_init (gsl_monte_vegas_state * state)
{
  state->stage = 0;
  state->alpha = 1.5;
  state->verbose = -1;
  state->iterations = 5;
  state->mode = GSL_VEGAS_MODE_IMPORTANCE;
  state->chisq = 0;
  state->bins = state->bins_max;
  state->ostream = stdout;

  return GSL_SUCCESS;
}

void
gsl_monte_vegas_free (gsl_monte_vegas_state * s)
{
  free (s->x);
  free (s->delx);
  free (s->d);
  free (s->xi);
  free (s->xin);
  free (s->weight);
  free (s->box);
  free (s->bin);
  free (s);
}

static void
init_box_coord (gsl_monte_vegas_state * s, coord box[])
{
  size_t i;

  size_t dim = s->dim;

  for (i = 0; i < dim; i++)
    {
      box[i] = 0;
    }
}

/* change_box_coord steps through the box coord like
   {0,0}, {0, 1}, {0, 2}, {0, 3}, {1, 0}, {1, 1}, {1, 2}, ...
*/
static int
change_box_coord (gsl_monte_vegas_state * s, coord box[])
{
  int j = s->dim - 1;

  int ng = s->boxes;

  while (j >= 0)
    {
      box[j] = (box[j] + 1) % ng;

      if (box[j] != 0)
        {
          return 1;
        }

      j--;
    }

  return 0;
}

static void
init_grid (gsl_monte_vegas_state * s, double xl[], double xu[], size_t dim)
{
  size_t j;

  double vol = 1.0;

  s->bins = 1;

  for (j = 0; j < dim; j++)
    {
      double dx = xu[j] - xl[j];
      s->delx[j] = dx;
      vol *= dx;

      COORD (s, 0, j) = 0.0;
      COORD (s, 1, j) = 1.0;
    }

  s->vol = vol;
}


static void
reset_grid_values (gsl_monte_vegas_state * s)
{
  size_t i, j;

  size_t dim = s->dim;
  size_t bins = s->bins;

  for (i = 0; i < bins; i++)
    {
      for (j = 0; j < dim; j++)
        {
          VALUE (s, i, j) = 0.0;
        }
    }
}

static void
accumulate_distribution (gsl_monte_vegas_state * s, coord bin[], double y)
{
  size_t j;
  size_t dim = s->dim;

  for (j = 0; j < dim; j++)
    {
      int i = bin[j];
      VALUE (s, i, j) += 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)
{
  /* Use the random number generator r to return a random position x
     in a given box.  The value of bin gives the bin location of the
     random position (there may be several bins within a given box) */

  double vol = 1.0;

  size_t j;

  size_t dim = s->dim;
  size_t bins = s->bins;
  size_t boxes = s->boxes;

  DISCARD_POINTER(xu); /* prevent warning about unused parameter */

  for (j = 0; j < dim; ++j)
    {
      /* box[j] + ran gives the position in the box units, while z
         is the position in bin units.  */

      double z = ((box[j] + gsl_rng_uniform_pos (r)) / boxes) * bins;

      int k = z;

      double y, bin_width;

      bin[j] = k;

      if (k == 0)
        {
          bin_width = COORD (s, 1, j);
          y = z * bin_width;
        }
      else
        {
          bin_width = COORD (s, k + 1, j) - COORD (s, k, j);
          y = COORD (s, k, j) + (z - k) * bin_width;
        }

      x[j] = xl[j] + y * s->delx[j];

      vol *= bin_width;
    }

  *bin_vol = vol;
}


static void
resize_grid (gsl_monte_vegas_state * s, unsigned int bins)
{
  size_t j, k;
  size_t dim = s->dim;

  /* weight is ratio of bin sizes */

  double pts_per_bin = (double) s->bins / (double) bins;

  for (j = 0; j < dim; j++)
    {
      double xold;
      double xnew = 0;
      double dw = 0;
      int i = 1;

      for (k = 1; k <= s->bins; k++)
        {
          dw += 1.0;
          xold = xnew;
          xnew = COORD (s, k, j);

          for (; dw > pts_per_bin; i++)
            {
              dw -= pts_per_bin;
              NEW_COORD (s, i) = xnew - (xnew - xold) * dw;
            }
        }

      for (k = 1 ; k < bins; k++)
        {
          COORD(s, k, j) = NEW_COORD(s, k);
        }

      COORD (s, bins, j) = 1;
    }

  s->bins = bins;
}

static void
refine_grid (gsl_monte_vegas_state * s)
{
  size_t i, j, k;
  size_t dim = s->dim;
  size_t bins = s->bins;

  for (j = 0; j < dim; j++)
    {
      double grid_tot_j, tot_weight;
      double * weight = s->weight;

      double oldg = VALUE (s, 0, j);
      double newg = VALUE (s, 1, j);

      VALUE (s, 0, j) = (oldg + newg) / 2;
      grid_tot_j = VALUE (s, 0, j);

      /* This implements gs[i][j] = (gs[i-1][j]+gs[i][j]+gs[i+1][j])/3 */

      for (i = 1; i < bins - 1; i++)
        {
          double rc = oldg + newg;
          oldg = newg;
          newg = VALUE (s, i + 1, j);
          VALUE (s, i, j) = (rc + newg) / 3;
          grid_tot_j += VALUE (s, i, j);
        }
      VALUE (s, bins - 1, j) = (newg + oldg) / 2;

      grid_tot_j += VALUE (s, bins - 1, j);

      tot_weight = 0;

      for (i = 0; i < bins; i++)
        {
          weight[i] = 0;

          if (VALUE (s, i, j) > 0)
            {
              oldg = grid_tot_j / VALUE (s, i, j);
              /* damped change */
              weight[i] = pow (((oldg - 1) / oldg / log (oldg)), s->alpha);
            }

          tot_weight += weight[i];

#ifdef DEBUG
          printf("weight[%d] = %g\n", i, weight[i]);
#endif
        }

      {
        double pts_per_bin = tot_weight / bins;

        double xold;
        double xnew = 0;
        double dw = 0;
        i = 1;

        for (k = 0; k < bins; k++)
          {
            dw += weight[k];
            xold = xnew;
            xnew = COORD (s, k + 1, j);

            for (; dw > pts_per_bin; i++)
              {
                dw -= pts_per_bin;
                NEW_COORD (s, i) = xnew - (xnew - xold) * dw / weight[k];
              }
          }

        for (k = 1 ; k < bins ; k++)
          {
            COORD(s, k, j) = NEW_COORD(s, k);
          }

        COORD (s, bins, j) = 1;
      }
    }
}


static void
print_lim (gsl_monte_vegas_state * state,
           double xl[], double xu[], unsigned long dim)
{
  unsigned long j;

  fprintf (state->ostream, "The limits of integration are:\n");
  for (j = 0; j < dim; ++j)
    fprintf (state->ostream, "\nxl[%lu]=%f    xu[%lu]=%f", j, xl[j], j, xu[j]);
  fprintf (state->ostream, "\n");
  fflush (state->ostream);
}

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)
{
  fprintf (state->ostream,
           "\nnum_dim=%lu, calls=%lu, it_num=%d, max_it_num=%d ",
           num_dim, calls, it_num, state->iterations);
  fprintf (state->ostream,
           "verb=%d, alph=%.2f,\nmode=%d, bins=%d, boxes=%d\n",
           state->verbose, state->alpha, state->mode, bins, boxes);
  fprintf (state->ostream,
           "\n       single.......iteration                   ");
  fprintf (state->ostream, "accumulated......results   \n");

  fprintf (state->ostream,
           "iteration     integral    sigma             integral   ");
  fprintf (state->ostream, "      sigma     chi-sq/it\n\n");
  fflush (state->ostream);

}

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)
{
  fprintf (state->ostream,
           "%4d        %6.4e %10.2e          %6.4e      %8.2e  %10.2e\n",
           itr, res, err, cum_res, cum_err, chi_sq);
  fflush (state->ostream);
}

static void
print_dist (gsl_monte_vegas_state * state, unsigned long dim)
{
  unsigned long i, j;
  int p = state->verbose;
  if (p < 1)
    return;

  for (j = 0; j < dim; ++j)
    {
      fprintf (state->ostream, "\n axis %lu \n", j);
      fprintf (state->ostream, "      x   g\n");
      for (i = 0; i < state->bins; i++)
        {
          fprintf (state->ostream, "weight [%11.2e , %11.2e] = ", 
                   COORD (state, i, j), COORD(state,i+1,j));
          fprintf (state->ostream, " %11.2e\n", VALUE (state, i, j));

        }
      fprintf (state->ostream, "\n");
    }
  fprintf (state->ostream, "\n");
  fflush (state->ostream);

}

static void
print_grid (gsl_monte_vegas_state * state, unsigned long dim)
{
  unsigned long i, j;
  int p = state->verbose;
  if (p < 1)
    return;

  for (j = 0; j < dim; ++j)
    {
      fprintf (state->ostream, "\n axis %lu \n", j);
      fprintf (state->ostream, "      x   \n");
      for (i = 0; i <= state->bins; i++)
        {
          fprintf (state->ostream, "%11.2e", COORD (state, i, j));
          if (i % 5 == 4)
            fprintf (state->ostream, "\n");
        }
      fprintf (state->ostream, "\n");
    }
  fprintf (state->ostream, "\n");
  fflush (state->ostream);

}