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++) {
302 const ubvector_int &lbox,
const vec_t &xl,
309 for (j=0;j<
dim;++j) {
319 double z=((lbox[j]+rdn)/boxes)*
bins;
331 bin_width=xi[(k+1)*dim+j]-xi[k*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];
405 d[(bins-1)*dim+j]=(newg+oldg)/2;
407 grid_tot_j+=d[(bins-1)*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];
417 weight[i]=pow (((oldg-1)/oldg/log (oldg)), alpha);
420 tot_weight+=weight[i];
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++) {
438 (xin[i])=xnew-(xnew-xold)*dw/weight[k];
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=" 475 << iterations <<
", verb=" << this->
verbose <<
", alph=" << alpha
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;
574 d.resize(bins_max*ldim);
575 xi.resize((bins_max+1)*ldim);
576 xin.resize(bins_max+1);
577 weight.resize(bins_max);
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;
664 if (2*lboxes >= bins_max) {
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);
679 calls_per_box=((
unsigned int)(GSL_MAX(calls/tot_boxes,2)));
680 calls=((size_t)( calls_per_box*tot_boxes));
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);
801 chisq*=(samples-2.0);
802 chisq+=(wgt/(1+(wgt/lsum_wgts)))*q*q;
803 chisq/=(samples-1.0);
807 cum_int+=(intgrl-cum_int)/(it+1.0);
812 print_res(it_num,intgrl,sig,cum_int,cum_sig,chisq);
813 if (it+1 == iterations && this->
verbose > 1) {
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 ubvector_int box
The boxes for each direction.
rng_gsl_uniform_real rng_dist
The random number distribution.
virtual void print_head(unsigned long num_dim, unsigned long calls, unsigned int lit_num, unsigned int lbins, unsigned int lboxes)
Print header.
std::function< double(size_t, const boost::numeric::ublas::vector< double > &)> multi_funct
Multi-dimensional function typedef.
ubvector delx
The iteration length in each direction.
double chisq
The chi-squared per degree of freedom for the weighted estimate of the integral.
ubvector_int bin
The bins for each direction.
int change_box_coord(ubvector_int &boxt)
Change box coordinates.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
ubvector xin
Storage for grid refinement.
unsigned int it_start
The starting iteration number.
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.
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.
virtual void print_grid(unsigned long ldim)
Print grid.
virtual void refine_grid()
Refine the grid.
invalid argument supplied by user
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 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.
void accumulate_distribution(ubvector_int &lbin, double y)
Add the most recently generated result to the distribution.
virtual void print_dist(unsigned long ldim)
Print distribution.
unsigned int iterations
Set the number of iterations (default 5)
size_t dim
Number of dimensions.
unsigned int it_num
The total number of iterations.
virtual void reset_grid_values()
Reset grid values.
virtual void init_grid(const vec_t &xl, const vec_t &xu, size_t ldim)
Initialize grid.
ubvector weight
The weight for each bin.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
std::string dtos(double x, int prec=6, bool auto_prec=false)
Convert a double to a string.
virtual const char * type()
Return string denoting type ("mcarlo_vegas")
virtual void resize_grid(unsigned int lbins)
Resize the grid.
unsigned long n_points
Number of integration points (default 1000)
virtual int allocate(size_t ldim)
Allocate memory.
double result
Result from last iteration.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
unsigned int calls_per_box
Number of function calls per box.
Monte-Carlo integration [abstract base].
double vol
The volume of the current bin.
virtual void init_box_coord(ubvector_int &boxt)
Initialize box coordinates.
static const size_t bins_max
Maximum number of bins.
virtual void print_res(unsigned int itr, double res, double err, double cum_res, double cum_err, double chi_sq)
Print results.
vec_t x
Point for function evaluation.
double interror
The uncertainty for the last integration computation.
Multidimensional integration using Vegas Monte Carlo (GSL)
virtual void print_lim(const vec_t &xl, const vec_t &xu, unsigned long ldim)
Print limits of integration.
double sigma
Uncertainty from last iteration.
unsigned int boxes
Number of boxes.
std::ostream * outs
The output stream to send output information (default std::cout)
unsigned int bins
Number of bins.
unsigned int samples
Number of samples for computing chi squared.
ubvector xi
Boundaries for each bin.
double alpha
The stiffness of the rebinning algorithm (default 1.5)
rng_t rng
The random number generator.