Go to the documentation of this file.
42 #ifndef O2SCL_MCARLO_VEGAS_H
43 #define O2SCL_MCARLO_VEGAS_H
51 #include <gsl/gsl_math.h>
52 #include <gsl/gsl_monte.h>
53 #include <gsl/gsl_machine.h>
54 #include <gsl/gsl_monte_vegas.h>
56 #include <o2scl/mcarlo.h>
58 #ifndef DOXYGEN_NO_O2NS
126 static const int mode_importance=1;
127 static const int mode_importance_only=0;
128 static const int mode_stratified=-1;
160 #ifndef DOXYGEN_INTERNAL
214 for (i=0;i<
dim;i++) {
236 boxt[j]=(boxt[j]+1) % ng;
248 virtual void init_grid(
const vec_t &xl,
const vec_t &xu,
size_t ldim) {
253 for (j=0;j<ldim;j++) {
254 double dx=xu[j]-xl[j];
269 for (i=0;i<
bins;i++) {
270 for (j=0;j<
dim;j++) {
285 for (j=0;j<
dim;j++) {
309 for (j=0;j<
dim;++j) {
332 y=
xi[k*
dim+j]+(z-k)*bin_width;
335 lx[j]=xl[j]+y*
delx[j];
351 double pts_per_bin=(double)
bins/(
double) lbins;
353 for (j=0;j<
dim;j++) {
359 for (k=1;k <=
bins;k++) {
364 for (;dw > pts_per_bin;i++) {
366 (
xin[i])=xnew-(xnew-xold)*dw;
370 for (k=1 ;k<lbins;k++) {
387 for (j=0;j<
dim;j++) {
388 double grid_tot_j, tot_weight;
391 double newg=
d[
dim+j];
398 for (i=1;i<
bins-1;i++) {
402 d[i*
dim+j]=(rc+newg)/3;
403 grid_tot_j+=
d[i*
dim+j];
411 for (i=0;i<
bins;i++) {
414 if (
d[i*
dim+j] > 0) {
415 oldg=grid_tot_j/
d[i*
dim+j];
424 double pts_per_bin=tot_weight/
bins;
431 for (k=0;k<
bins;k++) {
434 xnew=
xi[(k+1)*
dim+j];
436 for (;dw > pts_per_bin;i++) {
442 for (k=1 ;k<
bins ;k++) {
453 virtual void print_lim(
const vec_t &xl,
const vec_t &xu,
454 unsigned long ldim) {
458 (*outs) <<
"The limits of integration are:\n" << std::endl;
459 for (j=0;j<ldim;++j) {
460 (*outs) <<
"xl[" << j <<
"]=" << xl[j] <<
" xu[" << j
461 <<
"]=" << xu[j] << std::endl;
463 (*outs) << std::endl;
469 virtual void print_head(
unsigned long num_dim,
unsigned long calls,
470 unsigned int lit_num,
unsigned int lbins,
471 unsigned int lboxes) {
473 (*outs) <<
"num_dim=" << num_dim <<
", calls=" << calls
474 <<
", it_num=" << lit_num <<
", max_it_num="
476 <<
",\n mode=" << mode <<
", boxes=" << lboxes
477 <<
"\n\n single.......iteration "
478 <<
"accumulated......results\n"
479 <<
"iter integral sigma integral "
480 <<
" sigma chi-sq/it" << std::endl;
487 double err,
double cum_res,
double cum_err,
489 (*outs).precision(5);
490 (*outs) << itr <<
" ";
491 outs->setf(std::ios::showpos);
492 (*outs) << res <<
" ";
493 outs->unsetf(std::ios::showpos);
494 (*outs) << err <<
" ";
495 outs->setf(std::ios::showpos);
496 (*outs) << cum_res <<
" ";
497 outs->unsetf(std::ios::showpos);
498 (*outs) << cum_err <<
" " << chi_sq << std::endl;
499 (*outs).precision(6);
511 for (j=0;j<ldim;++j) {
512 (*outs) <<
"\n Axis " << j << std::endl;
513 (*outs) <<
" x g" << std::endl;
514 outs->setf(std::ios::showpos);
515 for (i=0;i<
bins;i++) {
516 (*outs) <<
"weight [ " << (
xi[(i)*
dim+(j)]) <<
" , "
517 <<
xi[(i+1)*
dim+j] <<
" ] = ";
518 (*outs) <<
" " << (
d[(i)*
dim+(j)]) << std::endl;
520 outs->unsetf(std::ios::showpos);
521 (*outs) << std::endl;
523 (*outs) << std::endl;
535 for (j=0;j<ldim;++j) {
536 (*outs) <<
"\n Axis " << j << std::endl;
537 (*outs) <<
" x " << std::endl;
538 outs->setf(std::ios::showpos);
539 for (i=0;i <=
bins;i++) {
540 (*outs) << (
xi[(i)*
dim+(j)]) <<
" ";
542 (*outs) << std::endl;
545 (*outs) << std::endl;
546 outs->unsetf(std::ios::showpos);
548 (*outs) << std::endl;
564 mode=mode_importance;
609 const vec_t &xl,
const vec_t &xu,
610 double &res,
double &err) {
614 double cum_int, cum_sig;
617 for (i=0;i<
dim;i++) {
618 if (xu[i] <= xl[i]) {
619 std::string serr=
"Upper limit, "+
dtos(xu[i])+
", must be greater "+
620 "than lower limit, "+
dtos(xl[i])+
", in mcarlo_vegas::"+
621 "vegas_minteg_err().";
625 if (xu[i]-xl[i] > GSL_DBL_MAX) {
626 O2SCL_ERR2(
"Range of integration is too large, please rescale ",
627 "in mcarlo_vegas::vegas_minteg_err().",
exc_einval);
650 unsigned int lboxes=1;
652 if (mode != mode_importance_only) {
661 lboxes=((
unsigned int)(floor(pow(calls/2.0,1.0/
dim))));
662 mode=mode_importance;
666 int box_per_bin=GSL_MAX(lboxes/
bins_max,1);
668 if (lboxes/box_per_bin<
bins_max) lbins=lboxes/box_per_bin;
670 lboxes=box_per_bin*lbins;
672 mode=mode_stratified;
678 double tot_boxes=pow((
double)lboxes,(
double)
dim);
683 jac=
vol*pow((
double) lbins, (
double)
dim)/calls;
707 double intgrl=0.0, intgrl_sq=0.0;
709 double wgt, var, sig;
719 volatile double m=0, q=0;
722 for (k=0;k<lcalls_per_box;k++) {
723 double fval, bin_vol;
728 fval*=jacbin*bin_vol;
735 q+=dt*dt*(k/(k+1.0));
738 if (mode != mode_stratified) {
739 double f_sq=fval*fval;
744 intgrl+=m*lcalls_per_box;
746 f_sq_sum=q*lcalls_per_box;
750 if (mode == mode_stratified) {
758 var=tss/(lcalls_per_box-1.0);
762 }
else if (sum_wgts>0) {
768 intgrl_sq=intgrl*intgrl;
776 double lsum_wgts=sum_wgts;
777 double m=(sum_wgts > 0) ? (wtd_int_sum/sum_wgts) : 0;
782 wtd_int_sum+=intgrl*wgt;
783 chi_sum+=intgrl_sq*wgt;
785 cum_int= wtd_int_sum/sum_wgts;
786 cum_sig=sqrt(1/sum_wgts);
802 chisq+=(wgt/(1+(wgt/lsum_wgts)))*q*q;
807 cum_int+=(intgrl-cum_int)/(it+1.0);
845 virtual int minteg_err(func_t &func,
size_t ndim,
const vec_t &a,
846 const vec_t &b,
double &res,
double &err) {
858 virtual double minteg(func_t &func,
size_t ndim,
const vec_t &a,
866 virtual const char *
type() {
return "mcarlo_vegas"; }
870 #ifndef DOXYGEN_NO_O2NS
double interror
The uncertainty for the last integration computation.
std::string dtos(const fp_t &x, int prec=6, bool auto_prec=false)
Convert a double to a string.
unsigned long n_points
Number of integration points (default 1000)
virtual const char * type()
Return string denoting type ("mcarlo_vegas")
virtual void reset_grid_values()
Reset grid values.
rng_gsl_uniform_real rng_dist
The random number distribution.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
virtual void init_box_coord(ubvector_int &boxt)
Initialize box coordinates.
virtual void print_lim(const vec_t &xl, const vec_t &xu, unsigned long ldim)
Print limits of integration.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
double sigma
Uncertainty from last iteration.
vec_t x
Point for function evaluation.
unsigned int calls_per_box
Number of function calls per box.
virtual void print_head(unsigned long num_dim, unsigned long calls, unsigned int lit_num, unsigned int lbins, unsigned int lboxes)
Print header.
ubvector xi
Boundaries for each bin.
unsigned int bins
Number of bins.
static const size_t bins_max
Maximum number of bins.
unsigned int boxes
Number of boxes.
double vol
The volume of the current bin.
virtual void print_res(unsigned int itr, double res, double err, double cum_res, double cum_err, double chi_sq)
Print results.
rng_gsl rng
The random number generator.
Multidimensional integration using Vegas Monte Carlo (GSL)
std::function< double(size_t, const boost::numeric::ublas::vector< double > &)> multi_funct
Multi-dimensional function typedef in src/base/multi_funct.h.
ubvector delx
The iteration length in each direction.
@ exc_einval
invalid argument supplied by user
std::ostream * outs
The output stream to send output information (default std::cout)
virtual int minteg_err(func_t &func, size_t ndim, const vec_t &a, const vec_t &b, double &res, double &err)
Integrate function func from x=a to x=b.
int change_box_coord(ubvector_int &boxt)
Change box coordinates.
Monte-Carlo integration [abstract base].
double alpha
The stiffness of the rebinning algorithm (default 1.5)
virtual double minteg(func_t &func, size_t ndim, const vec_t &a, const vec_t &b)
Integrate function func over the hypercube from to for ndim-1.
virtual void print_grid(unsigned long ldim)
Print grid.
ubvector_int box
The boxes for each direction.
unsigned int samples
Number of samples for computing chi squared.
double chisq
The chi-squared per degree of freedom for the weighted estimate of the integral.
virtual void print_dist(unsigned long ldim)
Print distribution.
ubvector xin
Storage for grid refinement.
unsigned int iterations
Set the number of iterations (default 5)
ubvector_int bin
The bins for each direction.
void random_point(vec_t &lx, ubvector_int &lbin, double &bin_vol, const ubvector_int &lbox, const vec_t &xl, const vec_t &xu)
Generate a random position in a given box.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
unsigned int it_start
The starting iteration number.
virtual void resize_grid(unsigned int lbins)
Resize the grid.
virtual void refine_grid()
Refine the grid.
void accumulate_distribution(ubvector_int &lbin, double y)
Add the most recently generated result to the distribution.
ubvector weight
The weight for each bin.
virtual int vegas_minteg_err(int stage, func_t &func, size_t ndim, const vec_t &xl, const vec_t &xu, double &res, double &err)
Integrate function func from x=a to x=b.
virtual void init_grid(const vec_t &xl, const vec_t &xu, size_t ldim)
Initialize grid.
size_t dim
Number of dimensions.
virtual int allocate(size_t ldim)
Allocate memory.
double result
Result from last iteration.
unsigned int it_num
The total number of iterations.
Documentation generated with Doxygen. Provided under the
GNU Free Documentation License (see License Information).