diff --git a/CepGen/Core/Generator.cpp b/CepGen/Core/Generator.cpp
index d14629e..3fc35cd 100644
--- a/CepGen/Core/Generator.cpp
+++ b/CepGen/Core/Generator.cpp
@@ -1,184 +1,179 @@
#include "CepGen/Generator.h"
#include "CepGen/Parameters.h"
#include "CepGen/Version.h"
#include "CepGen/Core/Integrator.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Core/Timer.h"
+#include "CepGen/Core/utils.h"
#include "CepGen/Physics/PDG.h"
#include "CepGen/Processes/ProcessesHandler.h"
#include "CepGen/Event/Event.h"
#include <fstream>
#include <chrono>
#include <atomic>
namespace cepgen
{
namespace utils { std::atomic<int> gSignal; }
Generator::Generator() :
parameters_( new Parameters ), result_( -1. ), result_error_( -1. )
{
CG_DEBUG( "Generator:init" ) << "Generator initialized";
try {
printHeader();
} catch ( const Exception& e ) {
e.dump();
}
//--- random number initialization
std::chrono::system_clock::time_point time = std::chrono::system_clock::now();
srandom( time.time_since_epoch().count() );
}
Generator::Generator( Parameters* ip ) :
parameters_( ip ), result_( -1. ), result_error_( -1. )
{}
Generator::~Generator()
{
if ( parameters_->generation().enabled
&& parameters_->process()
&& parameters_->numGeneratedEvents() > 0 )
CG_INFO( "Generator" )
<< "Mean generation time / event: "
<< parameters_->totalGenerationTime()*1.e3/parameters_->numGeneratedEvents()
<< " ms.";
}
size_t
Generator::numDimensions() const
{
if ( !parameters_->process() )
return 0;
return parameters_->process()->numDimensions();
}
void
Generator::clearRun()
{
if ( parameters_->process() ) {
parameters_->process()->first_run = true;
parameters_->process()->addEventContent();
parameters_->process()->setKinematics( parameters_->kinematics );
}
result_ = result_error_ = -1.;
{
std::ostringstream os;
for ( const auto& pr : cepgen::proc::ProcessesHandler::get().modules() )
os << " " << pr;
CG_DEBUG( "Generator:clearRun" ) << "Processes handled:" << os.str() << ".";
}
}
void
Generator::setParameters( const Parameters& ip )
{
parameters_.reset( new Parameters( (Parameters&)ip ) ); // copy constructor
}
void
Generator::printHeader()
{
std::string tmp;
std::ostringstream os; os << "version " << version() << std::endl;
std::ifstream hf( "README" );
if ( !hf.good() )
throw CG_WARNING( "Generator" ) << "Failed to open README file.";
while ( true ) {
if ( !hf.good() ) break;
getline( hf, tmp );
os << "\n " << tmp;
}
hf.close();
CG_INFO( "Generator" ) << os.str();
}
double
Generator::computePoint( double* x )
{
double res = integrand::eval( x, numDimensions(), (void*)parameters_.get() );
std::ostringstream os;
- for ( unsigned int i = 0; i < numDimensions(); ++i )
+ for ( size_t i = 0; i < numDimensions(); ++i )
os << x[i] << " ";
CG_DEBUG( "Generator:computePoint" )
<< "Result for x[" << numDimensions() << "] = { " << os.str() << "}:\n\t"
<< res << ".";
return res;
}
void
Generator::computeXsection( double& xsec, double& err )
{
CG_INFO( "Generator" ) << "Starting the computation of the process cross-section.";
integrate();
xsec = result_;
err = result_error_;
- //CG_INFO( "Generator" ) << "Total cross section: "
if ( xsec < 1.e-2 )
CG_INFO( "Generator" )
<< "Total cross section: " << xsec*1.e3 << " +/- " << err*1.e3 << " fb.";
else if ( xsec < 0.5e3 )
CG_INFO( "Generator" )
<< "Total cross section: " << xsec << " +/- " << err << " pb.";
else if ( xsec < 0.5e6 )
CG_INFO( "Generator" )
<< "Total cross section: " << xsec*1.e-3 << " +/- " << err*1.e-3 << " nb.";
else if ( xsec < 0.5e9 )
CG_INFO( "Generator" )
<< "Total cross section: " << xsec*1.e-6 << " +/- " << err*1.e-6 << " µb.";
else
CG_INFO( "Generator" )
<< "Total cross section: " << xsec*1.e-9 << " +/- " << err*1.e-9 << " mb.";
}
void
Generator::integrate()
{
clearRun();
// first destroy and recreate the integrator instance
- if ( !integrator_ )
- integrator_ = std::unique_ptr<Integrator>( new Integrator( numDimensions(), integrand::eval, parameters_.get() ) );
- else if ( integrator_->dimensions() != numDimensions() )
- integrator_.reset( new Integrator( numDimensions(), integrand::eval, parameters_.get() ) );
+ if ( !integrator_ || integrator_->dimensions() != numDimensions() )
+ integrator_.reset( new Integrator( numDimensions(), integrand::eval, *parameters_ ) );
CG_DEBUG( "Generator:newInstance" )
<< "New integrator instance created\n\t"
<< "Considered topology: " << parameters_->kinematics.mode << " case\n\t"
<< "Will proceed with " << numDimensions() << "-dimensional integration.";
- const int res = integrator_->integrate( result_, result_error_ );
- if ( res != 0 )
- throw CG_FATAL( "Generator" )
- << "Error while computing the cross-section!\n\t"
- << "GSL error: " << gsl_strerror( res ) << ".";
+ integrator_->integrate( result_, result_error_ );
}
std::shared_ptr<Event>
Generator::generateOneEvent()
{
integrator_->generateOne();
return parameters_->process()->last_event;
}
void
Generator::generate( std::function<void( const Event&, unsigned long )> callback )
{
const utils::Timer tmr;
CG_INFO( "Generator" )
<< parameters_->generation().maxgen << " events will be generated.";
integrator_->generate( parameters_->generation().maxgen, callback );
const double gen_time_s = tmr.elapsed();
CG_INFO( "Generator" )
- << parameters_->numGeneratedEvents() << " events generated "
+ << parameters_->numGeneratedEvents() << " event" << utils::s( parameters_->numGeneratedEvents() )
+ << " generated "
<< "in " << gen_time_s << " s "
<< "(" << gen_time_s/parameters_->numGeneratedEvents()*1.e3 << " ms/event).";
}
}
diff --git a/CepGen/Core/Integrator.cpp b/CepGen/Core/Integrator.cpp
index 36b0bb3..1acd9d9 100644
--- a/CepGen/Core/Integrator.cpp
+++ b/CepGen/Core/Integrator.cpp
@@ -1,491 +1,494 @@
#include "CepGen/Core/Integrator.h"
#include "CepGen/Core/GridParameters.h"
#include "CepGen/Core/utils.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Parameters.h"
#include "CepGen/Processes/GenericProcess.h"
#include "CepGen/Hadronisers/GenericHadroniser.h"
#include "CepGen/Event/Event.h"
#include <thread>
#include <math.h>
#include <gsl/gsl_monte_miser.h>
#define COORD(s,i,j) ((s)->xi[(i)*(s)->dim + (j)])
namespace cepgen
{
- Integrator::Integrator( unsigned int ndim, double integrand( double*, size_t, void* ), Parameters* params ) :
- ps_bin_( INVALID_BIN ), input_params_( params, []( Parameters* ){} ),
- function_( new gsl_monte_function{ integrand, ndim, (void*)input_params_.get() } ),
- rng_( gsl_rng_alloc( input_params_->integration().rng_engine ), gsl_rng_free ),
+ Integrator::Integrator( unsigned int ndim, double integrand( double*, size_t, void* ), Parameters& params ) :
+ ps_bin_( INVALID_BIN ), input_params_( params ),
+ function_( new gsl_monte_function{ integrand, ndim, (void*)&input_params_ } ),
+ rng_( gsl_rng_alloc( input_params_.integration().rng_engine ) ),
grid_( new GridParameters( ndim ) )
{
//--- initialise the random number generator
- unsigned long seed = ( input_params_->integration().rng_seed > 0 )
- ? input_params_->integration().rng_seed
+ unsigned long seed = ( input_params_.integration().rng_seed > 0 )
+ ? input_params_.integration().rng_seed
: time( nullptr ); // seed with time
gsl_rng_set( rng_.get(), seed );
//--- a bit of printout for debugging
CG_DEBUG( "Integrator:build" )
<< "Number of integration dimensions: " << function_->dim << ",\n\t"
- << "Number of function calls: " << input_params_->integration().ncvg << ",\n\t"
+ << "Number of function calls: " << input_params_.integration().ncvg << ",\n\t"
<< "Random numbers generator: " << gsl_rng_name( rng_.get() ) << ".";
- switch ( input_params_->integration().type ) {
+ switch ( input_params_.integration().type ) {
case IntegratorType::Vegas:
CG_DEBUG( "Integrator:build" ) << "Vegas parameters:\n\t"
- << "Number of iterations in Vegas: " << input_params_->integration().vegas.iterations << ",\n\t"
- << "α-value: " << input_params_->integration().vegas.alpha << ",\n\t"
- << "Verbosity: " << input_params_->integration().vegas.verbose << ",\n\t"
- << "Grid interpolation mode: " << (Integrator::VegasMode)input_params_->integration().vegas.mode << ".";
+ << "Number of iterations in Vegas: " << input_params_.integration().vegas.iterations << ",\n\t"
+ << "α-value: " << input_params_.integration().vegas.alpha << ",\n\t"
+ << "Verbosity: " << input_params_.integration().vegas.verbose << ",\n\t"
+ << "Grid interpolation mode: " << (Integrator::VegasMode)input_params_.integration().vegas.mode << ".";
break;
case IntegratorType::MISER:
CG_DEBUG( "Integrator:build" ) << "MISER parameters:\n\t"
- << "Number of calls: " << input_params_->integration().miser.min_calls << ", "
- << "per bisection: " << input_params_->integration().miser.min_calls_per_bisection << ",\n\t"
- << "Estimate fraction: " << input_params_->integration().miser.estimate_frac << ",\n\t"
- << "α-value: " << input_params_->integration().miser.alpha << ",\n\t"
- << "Dither: " << input_params_->integration().miser.dither << ".";
+ << "Number of calls: " << input_params_.integration().miser.min_calls << ", "
+ << "per bisection: " << input_params_.integration().miser.min_calls_per_bisection << ",\n\t"
+ << "Estimate fraction: " << input_params_.integration().miser.estimate_frac << ",\n\t"
+ << "α-value: " << input_params_.integration().miser.alpha << ",\n\t"
+ << "Dither: " << input_params_.integration().miser.dither << ".";
break;
case IntegratorType::plain:
break;
}
}
Integrator::~Integrator()
{}
//-----------------------------------------------------------------------------------------------
// cross section computation part
//-----------------------------------------------------------------------------------------------
- int
+ void
Integrator::integrate( double& result, double& abserr )
{
int res = -1;
//--- integration bounds
std::vector<double> x_low( function_->dim, 0. ), x_up( function_->dim, 1. );
//--- launch integration
- switch ( input_params_->integration().type ) {
+ switch ( input_params_.integration().type ) {
case IntegratorType::plain: {
std::unique_ptr<gsl_monte_plain_state,void(*)( gsl_monte_plain_state* )>
pln_state( gsl_monte_plain_alloc( function_->dim ), gsl_monte_plain_free );
res = gsl_monte_plain_integrate( function_.get(),
&x_low[0], &x_up[0],
- function_->dim, input_params_->integration().ncvg,
+ function_->dim, input_params_.integration().ncvg,
rng_.get(), pln_state.get(),
&result, &abserr );
} break;
case IntegratorType::Vegas: {
//----- warmup (prepare the grid)
- res = warmupVegas( x_low, x_up, 25000 );
+ warmupVegas( x_low, x_up, 25000 );
//----- integration
unsigned short it_chisq = 0;
do {
res = gsl_monte_vegas_integrate( function_.get(),
&x_low[0], &x_up[0],
- function_->dim, 0.2 * input_params_->integration().ncvg,
+ function_->dim, 0.2 * input_params_.integration().ncvg,
rng_.get(), veg_state_.get(),
&result, &abserr );
CG_LOG( "Integrator:integrate" )
<< "\t>> at call " << ( ++it_chisq ) << ": "
<< Form( "average = %10.6f "
"sigma = %10.6f chi2 = %4.3f.",
result, abserr,
gsl_monte_vegas_chisq( veg_state_.get() ) );
} while ( fabs( gsl_monte_vegas_chisq( veg_state_.get() )-1. )
- > input_params_->integration().vegas_chisq_cut-1. );
+ > input_params_.integration().vegas_chisq_cut-1. );
CG_DEBUG( "Integrator:integrate" )
<< "Vegas grid information:\n\t"
<< "ran for " << veg_state_->dim << " dimensions, and generated " << veg_state_->bins_max << " bins.\n\t"
<< "Integration volume: " << veg_state_->vol << ".";
grid_->r_boxes = std::pow( veg_state_->bins, function_->dim );
} break;
case IntegratorType::MISER: {
std::unique_ptr<gsl_monte_miser_state,void(*)( gsl_monte_miser_state* )>
mis_state( gsl_monte_miser_alloc( function_->dim ), gsl_monte_miser_free );
- gsl_monte_miser_params_set( mis_state.get(), &input_params_->integration().miser );
+ gsl_monte_miser_params_set( mis_state.get(), &input_params_.integration().miser );
res = gsl_monte_miser_integrate( function_.get(),
&x_low[0], &x_up[0],
- function_->dim, input_params_->integration().ncvg,
+ function_->dim, input_params_.integration().ncvg,
rng_.get(), mis_state.get(),
&result, &abserr );
} break;
}
- input_params_->integration().result = result;
- input_params_->integration().err_result = abserr;
+ input_params_.integration().result = result;
+ input_params_.integration().err_result = abserr;
- if ( input_params_->hadroniser() )
- input_params_->hadroniser()->setCrossSection( result, abserr );
+ if ( input_params_.hadroniser() )
+ input_params_.hadroniser()->setCrossSection( result, abserr );
- return res;
+ if ( res != GSL_SUCCESS )
+ throw CG_FATAL( "Integrator:integrate" )
+ << "Error while computing the cross-section!\n\t"
+ << "GSL error: " << gsl_strerror( res ) << ".";
}
- int
+ void
Integrator::warmupVegas( std::vector<double>& x_low, std::vector<double>& x_up, unsigned int ncall )
{
// start by preparing the grid/state
veg_state_.reset( gsl_monte_vegas_alloc( function_->dim ) );
- gsl_monte_vegas_params_set( veg_state_.get(), &input_params_->integration().vegas );
+ gsl_monte_vegas_params_set( veg_state_.get(), &input_params_.integration().vegas );
// then perform a first integration with the given calls count
double result = 0., abserr = 0.;
- int res = gsl_monte_vegas_integrate( function_.get(),
+ const int res = gsl_monte_vegas_integrate( function_.get(),
&x_low[0], &x_up[0],
- function_->dim, ncall,
- rng_.get(), veg_state_.get(),
+ function_->dim, ncall, rng_.get(), veg_state_.get(),
&result, &abserr );
// ensure the operation was successful
if ( res != GSL_SUCCESS )
throw CG_ERROR( "Integrator:vegas" )
<< "Failed to warm-up the Vegas grid.\n\t"
<< "GSL error: " << gsl_strerror( res ) << ".";
CG_INFO( "Integrator:vegas" )
<< "Finished the Vegas warm-up.";
- return res;
}
//-----------------------------------------------------------------------------------------------
// events generation part
//-----------------------------------------------------------------------------------------------
void
Integrator::generateOne( std::function<void( const Event&, unsigned long )> callback )
{
if ( !grid_->gen_prepared )
computeGenerationParameters();
std::vector<double> x( function_->dim, 0. );
//--- correction cycles
if ( ps_bin_ != INVALID_BIN ) {
bool has_correction = false;
while ( !correctionCycle( x, has_correction ) ) {}
if ( has_correction ) {
storeEvent( x, callback );
return;
}
}
double weight = 0.;
//--- normal generation cycle
while ( true ) {
double y = -1.;
//----- select a and reject if fmax is too small
while ( true ) {
// ...
ps_bin_ = uniform() * grid_->max;
y = uniform() * grid_->f_max_global;
grid_->num[ps_bin_] += 1;
if ( y <= grid_->f_max[ps_bin_] )
break;
}
// shoot a point x in this bin
std::vector<unsigned short> grid_n = grid_->n_map.at( ps_bin_ );
for ( unsigned int i = 0; i < function_->dim; ++i )
x[i] = ( uniform() + grid_n[i] ) * GridParameters::INV_M_BIN;
// get weight for selected x value
weight = eval( x );
if ( weight <= 0. )
continue;
if ( weight > y )
break;
}
if ( weight <= grid_->f_max[ps_bin_] )
ps_bin_ = INVALID_BIN;
else {
//--- if weight is higher than local or global maximum,
// init correction cycle
grid_->f_max_old = grid_->f_max[ps_bin_];
grid_->f_max[ps_bin_] = weight;
grid_->f_max_diff = weight-grid_->f_max_old;
if ( weight > grid_->f_max_global )
grid_->f_max_global = weight;
grid_->correc = ( grid_->num[ps_bin_]-1. ) * grid_->f_max_diff / grid_->f_max_global - 1.;
CG_DEBUG("Integrator::generateOne")
<< "Correction " << grid_->correc << " will be applied for phase space bin " << ps_bin_ << ".";
}
// return with an accepted event
if ( weight > 0. )
storeEvent( x, callback );
}
void
Integrator::generate( unsigned long num_events, std::function<void( const Event&, unsigned long )> callback )
{
if ( num_events < 1 )
- num_events = input_params_->generation().maxgen;
+ num_events = input_params_.generation().maxgen;
try {
- while ( input_params_->numGeneratedEvents() < num_events )
+ while ( input_params_.numGeneratedEvents() < num_events )
generateOne( callback );
} catch ( const Exception& e ) { return; }
}
bool
Integrator::correctionCycle( std::vector<double>& x, bool& has_correction )
{
CG_DEBUG_LOOP( "Integrator:correction" )
<< "Correction cycles are started.\n\t"
<< "bin = " << ps_bin_ << "\t"
<< "correc = " << grid_->correc << "\t"
<< "corre2 = " << grid_->correc2 << ".";
if ( grid_->correc >= 1. )
grid_->correc -= 1.;
if ( uniform() < grid_->correc ) {
grid_->correc = -1.;
std::vector<double> xtmp( function_->dim );
// Select x values in phase space bin
const std::vector<unsigned short> grid_n = grid_->n_map.at( ps_bin_ );
for ( unsigned int k = 0; k < function_->dim; ++k )
xtmp[k] = ( uniform() + grid_n[k] ) * GridParameters::INV_M_BIN;
const double weight = eval( xtmp );
// Parameter for correction of correction
if ( weight > grid_->f_max[ps_bin_] ) {
grid_->f_max2 = std::max( grid_->f_max2, weight );
grid_->correc += 1.;
grid_->correc2 -= 1.;
}
// Accept event
if ( weight >= grid_->f_max_diff*uniform() + grid_->f_max_old ) {
x = xtmp;
has_correction = true;
return true;
}
return false;
}
// Correction if too big weight is found while correction
// (All your bases are belong to us...)
if ( grid_->f_max2 > grid_->f_max[ps_bin_] ) {
grid_->f_max_old = grid_->f_max[ps_bin_];
grid_->f_max[ps_bin_] = grid_->f_max2;
grid_->f_max_diff = grid_->f_max2-grid_->f_max_old;
grid_->correc = ( grid_->num[ps_bin_]-1. ) * grid_->f_max_diff / grid_->f_max_global;
if ( grid_->f_max2 >= grid_->f_max_global ) {
grid_->correc *= grid_->f_max2 / grid_->f_max_global;
grid_->f_max_global = grid_->f_max2;
}
grid_->correc -= grid_->correc2;
grid_->correc2 = 0.;
grid_->f_max2 = 0.;
return false;
}
return true;
}
bool
Integrator::storeEvent( const std::vector<double>& x, std::function<void( const Event&, unsigned long )> callback )
{
//--- start by computing the matrix element for that point
const double weight = eval( x );
//--- reject if unphysical
if ( weight <= 0. )
return false;
{
- if ( input_params_->numGeneratedEvents() % input_params_->generation().gen_print_every == 0 ) {
+ if ( input_params_.numGeneratedEvents() % input_params_.generation().gen_print_every == 0 ) {
CG_INFO( "Integrator:store" )
- << "Generated events: " << input_params_->numGeneratedEvents();
- input_params_->process()->last_event->dump();
+ << "Generated events: " << input_params_.numGeneratedEvents();
+ input_params_.process()->last_event->dump();
}
- const Event& last_event = *input_params_->process()->last_event;
+ const Event& last_event = *input_params_.process()->last_event;
if ( callback )
- callback( last_event, input_params_->numGeneratedEvents() );
- input_params_->addGenerationTime( last_event.time_total );
+ callback( last_event, input_params_.numGeneratedEvents() );
+ input_params_.addGenerationTime( last_event.time_total );
}
return true;
}
//-----------------------------------------------------------------------------------------------
// initial preparation run before the generation of unweighted events
//-----------------------------------------------------------------------------------------------
void
Integrator::computeGenerationParameters()
{
- input_params_->setStorage( false );
+ input_params_.setStorage( false );
- if ( input_params_->generation().treat
- && input_params_->integration().type != IntegratorType::Vegas ) {
+ if ( input_params_.generation().treat
+ && input_params_.integration().type != IntegratorType::Vegas ) {
CG_INFO( "Integrator:setGen" )
<< "Treat switched on without a proper Vegas grid; running a warm-up beforehand.";
std::vector<double> x_low( function_->dim, 0. ), x_up( function_->dim, 1. );
- int res = warmupVegas( x_low, x_up, 25000 );
- if ( res != GSL_SUCCESS )
+ try {
+ warmupVegas( x_low, x_up, 25000 );
+ } catch ( const Exception& ) {
throw CG_FATAL( "Integrator::setGen" )
<< "Failed to perform a Vegas warm-up.\n\t"
<< "Try to re-run while disabling integrand treatment...";
+ }
}
CG_INFO( "Integrator:setGen" )
- << "Preparing the grid (" << input_params_->generation().num_points << " points/bin) "
+ << "Preparing the grid (" << input_params_.generation().num_points << " points/bin) "
<< "for the generation of unweighted events.";
- const double inv_num_points = 1./input_params_->generation().num_points;
+ const double inv_num_points = 1./input_params_.generation().num_points;
std::vector<double> x( function_->dim, 0. );
std::vector<unsigned short> n( function_->dim, 0 );;
// ...
double sum = 0., sum2 = 0., sum2p = 0.;
//--- main loop
for ( unsigned int i = 0; i < grid_->max; ++i ) {
unsigned int jj = i;
for ( unsigned int j = 0; j < function_->dim; ++j ) {
unsigned int tmp = jj*GridParameters::INV_M_BIN;
n[j] = jj-tmp*GridParameters::M_BIN;
jj = tmp;
}
grid_->n_map[i] = n;
if ( CG_EXCEPT_MATCH( "Integrator:setGen", debugInsideLoop ) ) {
std::ostringstream os;
for ( const auto& ni : n )
os << ni << " ";
CG_DEBUG_LOOP( "Integrator:setGen" )
<< "n-vector for bin " << i << ": " << os.str();
}
double fsum = 0., fsum2 = 0.;
- for ( unsigned int j = 0; j < input_params_->generation().num_points; ++j ) {
+ for ( unsigned int j = 0; j < input_params_.generation().num_points; ++j ) {
for ( unsigned int k = 0; k < function_->dim; ++k )
x[k] = ( uniform()+n[k] ) * GridParameters::INV_M_BIN;
const double weight = eval( x );
grid_->f_max[i] = std::max( grid_->f_max[i], weight );
fsum += weight;
fsum2 += weight*weight;
}
const double av = fsum*inv_num_points, av2 = fsum2*inv_num_points, sig2 = av2-av*av;
sum += av;
sum2 += av2;
sum2p += sig2;
grid_->f_max_global = std::max( grid_->f_max_global, grid_->f_max[i] );
// per-bin debugging loop
if ( CG_EXCEPT_MATCH( "Integrator:setGen", debugInsideLoop ) ) {
const double sig = sqrt( sig2 );
const double eff = ( grid_->f_max[i] != 0. )
? grid_->f_max[i]/av
: 1.e4;
std::ostringstream os;
for ( unsigned int j = 0; j < function_->dim; ++j )
os << ( j != 0 ? ", " : "" ) << n[j];
CG_DEBUG_LOOP( "Integrator:setGen" )
<< "In iteration #" << i << ":\n\t"
<< "av = " << av << "\n\t"
<< "sig = " << sig << "\n\t"
<< "fmax = " << grid_->f_max[i] << "\n\t"
<< "eff = " << eff << "\n\t"
<< "n = (" << os.str() << ")";
}
} // end of main loop
const double inv_max = 1./grid_->max;
sum *= inv_max;
sum2 *= inv_max;
sum2p *= inv_max;
const double sig = sqrt( sum2-sum*sum ), sigp = sqrt( sum2p );
double eff1 = 0.;
for ( unsigned int i = 0; i < grid_->max; ++i )
eff1 += sum/grid_->max*grid_->f_max[i];
const double eff2 = sum/grid_->f_max_global;
CG_DEBUG( "Integrator:setGen" )
<< "Average function value = sum = " << sum << "\n\t"
<< "Average squared function value = sum2 = " << sum2 << "\n\t"
<< "Overall standard deviation = sig = " << sig << "\n\t"
<< "Average standard deviation = sigp = " << sigp << "\n\t"
<< "Maximum function value = f_max = " << grid_->f_max_global << "\n\t"
<< "Average inefficiency = eff1 = " << eff1 << "\n\t"
<< "Overall inefficiency = eff2 = " << eff2;
grid_->gen_prepared = true;
- input_params_->setStorage( true );
+ input_params_.setStorage( true );
CG_INFO( "Integrator:setGen" ) << "Grid prepared! Now launching the production.";
}
//------------------------------------------------------------------------------------------------
// helper / alias methods
//------------------------------------------------------------------------------------------------
unsigned short
Integrator::dimensions() const
{
if ( !function_ )
return 0;
return function_->dim;
}
double
Integrator::eval( const std::vector<double>& x )
{
- if ( !input_params_->generation().treat )
- return function_->f( (double*)&x[0], function_->dim, (void*)input_params_.get() );
+ if ( !input_params_.generation().treat )
+ return function_->f( (double*)&x[0], function_->dim, (void*)&input_params_ );
//--- treatment of the integration grid
double w = grid_->r_boxes;
std::vector<double> x_new( x.size() );
for ( unsigned short j = 0; j < function_->dim; ++j ) {
//--- find surrounding coordinates and interpolate
const double z = x[j]*veg_state_->bins;
const unsigned int id = z; // coordinate of point before
const double rel_pos = z-id; // position between coordinates (norm.)
const double bin_width = ( id == 0 )
? COORD( veg_state_, 1, j )
: COORD( veg_state_, id+1, j )-COORD( veg_state_, id, j );
//--- build new coordinate from linear interpolation
x_new[j] = COORD( veg_state_, id+1, j )-bin_width*( 1.-rel_pos );
w *= bin_width;
}
- return w*function_->f( (double*)&x_new[0], function_->dim, (void*)input_params_.get() );
+ return w*function_->f( (double*)&x_new[0], function_->dim, (void*)&input_params_ );
}
double
Integrator::uniform() const
{
return gsl_rng_uniform( rng_.get() );
}
//------------------------------------------------------------------------------------------------
std::ostream&
operator<<( std::ostream& os, const IntegratorType& type )
{
switch ( type ) {
case IntegratorType::plain:
return os << "plain";
case IntegratorType::Vegas:
return os << "Vegas";
case IntegratorType::MISER:
return os << "MISER";
}
return os;
}
std::ostream&
operator<<( std::ostream& os, const Integrator::VegasMode& mode )
{
switch ( mode ) {
case Integrator::VegasMode::importance:
return os << "importance";
case Integrator::VegasMode::importanceOnly:
return os << "importance-only";
case Integrator::VegasMode::stratified:
return os << "stratified";
}
return os;
}
}
diff --git a/CepGen/Core/Integrator.h b/CepGen/Core/Integrator.h
index 0f5e408..40259c3 100644
--- a/CepGen/Core/Integrator.h
+++ b/CepGen/Core/Integrator.h
@@ -1,112 +1,115 @@
#ifndef CepGen_Core_Integrator_h
#define CepGen_Core_Integrator_h
#include <vector>
#include <memory>
#include <functional>
#include <string.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_vegas.h>
#include <gsl/gsl_rng.h>
namespace cepgen
{
class Parameters;
class Event;
class GridParameters;
namespace utils { class Timer; }
/// Flavour of integration algorithm
enum class IntegratorType {
plain = 0, ///< Simple trial-and-error algorithm
Vegas = 1, ///< Vegas algorithm (G.P. Lepage, 1977 \cite Lepage:1977sw)
MISER = 2 ///< MISER stratified sampling algorithm
};
/// Monte-Carlo integrator instance
class Integrator
{
public:
enum class VegasMode { importance = 1, importanceOnly = 0, stratified = -1 };
/**
* Book the memory slots and structures for the integrator
* \note Three integration algorithms are currently supported:
* * the plain algorithm randomly sampling points in the phase space
* * the Vegas algorithm developed by P. Lepage, as documented in \cite Lepage:1977sw
* * the MISER algorithm developed by W.H. Press and G.R. Farrar, as documented in \cite Press:1989vk.
* \param[in] ndim Number of dimensions on which the function will be integrated
* \param[in] integrand Function to be integrated
* \param[inout] params Run parameters to define the phase space on which this integration is performed (embedded in an Parameters object)
*/
- Integrator( unsigned int ndim, double integrand(double*,size_t,void*), Parameters* params );
+ Integrator( unsigned int ndim, double integrand(double*,size_t,void*), Parameters& params );
/// Class destructor
~Integrator();
/**
* Algorithm to perform the n-dimensional Monte Carlo integration of a given function.
* \author This C++ implementation: GSL
* \param[out] result_ The cross section as integrated for the given phase space restrictions
* \param[out] abserr_ The error associated to the computed cross section
- * \return 0 if the integration was performed successfully
*/
- int integrate( double& result_, double& abserr_ );
+ void integrate( double& result_, double& abserr_ );
/// Dimensional size of the phase space
unsigned short dimensions() const;
/// Generate a single event
void generateOne( std::function<void( const Event&, unsigned long )> callback = nullptr );
/// Launch the event generation for a given number of events
void generate( unsigned long num_events = 0, std::function<void( const Event&, unsigned long )> callback = nullptr );
private:
/**
* Store the event characterized by its _ndim-dimensional point in
* the phase space to the output file
* \brief Store the event in the output file
* \param[in] x The d-dimensional point in the phase space defining the unique event to store
* \return A boolean stating whether or not the event could be saved
*/
bool storeEvent( const std::vector<double>& x, std::function<void( const Event&, unsigned long )> callback = nullptr );
/// Start the correction cycle on the grid
/// \param x Point in the phase space considered
/// \param has_correction Correction cycle started?
bool correctionCycle( std::vector<double>& x, bool& has_correction );
/// Prepare Vegas for an integration/event generation cycle
- int warmupVegas( std::vector<double>& x_low, std::vector<double>& x_up, unsigned int ncall );
+ void warmupVegas( std::vector<double>& x_low, std::vector<double>& x_up, unsigned int ncall );
/**
* Set all the generation mode variables and align them to the
* integration grid set while computing the cross-section
* \brief Prepare the class for events generation
*/
void computeGenerationParameters();
/// Generate a uniformly distributed (between 0 and 1) random number
double uniform() const;
/// Compute the function value at the given phase space point
double eval( const std::vector<double>& x );
/// Selected bin at which the function will be evaluated
int ps_bin_;
static constexpr int INVALID_BIN = -1;
/// List of parameters to specify the integration range and the
/// physics determining the phase space
- std::shared_ptr<Parameters> input_params_;
+ Parameters& input_params_;
/// GSL structure storing the function to be integrated by this
/// integrator instance (along with its parameters)
std::unique_ptr<gsl_monte_function> function_;
+ struct gsl_rng_deleter
+ {
+ inline void operator()( gsl_rng* rng ) { gsl_rng_free( rng ); }
+ };
/// Instance of random number generator service
- std::unique_ptr<gsl_rng,void(*)( gsl_rng* )> rng_;
+ std::unique_ptr<gsl_rng,gsl_rng_deleter> rng_;
/// Set of parameters for the integration/event generation grid
std::unique_ptr<GridParameters> grid_;
/// A trivial deleter for the Vegas integrator
struct gsl_monte_vegas_deleter
{
void operator()( gsl_monte_vegas_state* state ) {
gsl_monte_vegas_free( state );
}
};
/// A Vegas integrator state for integration (optional) and/or
/// "treated" event generation
std::unique_ptr<gsl_monte_vegas_state,gsl_monte_vegas_deleter> veg_state_;
};
std::ostream& operator<<( std::ostream&, const IntegratorType& );
std::ostream& operator<<( std::ostream&, const Integrator::VegasMode& );
}
#endif
diff --git a/CepGen/Core/utils.h b/CepGen/Core/utils.h
index 4728216..bb63dee 100644
--- a/CepGen/Core/utils.h
+++ b/CepGen/Core/utils.h
@@ -1,38 +1,50 @@
#ifndef CepGen_Core_utils_h
#define CepGen_Core_utils_h
#include <string>
+#include <vector>
+#include <numeric>
namespace cepgen
{
- /// Add a closing "s" when needed
- inline const char* s( unsigned short num ) { return ( num > 1 ) ? "s" : ""; }
/// Format a string using a printf style format descriptor.
std::string Form( const std::string fmt, ... );
/// Human-readable boolean printout
inline const char* yesno( const bool& test ) { return ( test ) ? "\033[32;1myes\033[0m" : "\033[31;1mno\033[0m"; }
//inline const char* boldify( const char* str ) { const std::string out = std::string( "\033[33;1m" ) + std::string( str ) + std::string( "\033[0m" ); return out.c_str(); }
/// Boldify a string for TTY-type output streams
inline std::string boldify( const std::string& str ) { return Form( "\033[1m%s\033[0m", str.c_str() ); }
/// Boldify a string for TTY-type output streams
inline std::string boldify( const char* str ) { return boldify( std::string( str ) ); }
/// Boldify a double floating point number for TTY-type output streams
inline std::string boldify( const double& dbl ) { return boldify( Form("%.2f", dbl ) ); }
/// Boldify an integer for TTY-type output streams
inline std::string boldify( const int& i ) { return boldify( Form("% d", i ) ); }
/// Boldify an unsigned integer for TTY-type output streams
inline std::string boldify( const unsigned int& ui ) { return boldify( Form("%d", ui ) ); }
/// Boldify an unsigned long integer for TTY-type output streams
inline std::string boldify( const unsigned long& ui ) { return boldify( Form("%lu", ui ) ); }
/// TTY-type enumeration of colours
enum class Colour { gray = 30, red = 31, green = 32, yellow = 33, blue = 34, purple = 35 };
/// Colourise a string for TTY-type output streams
inline std::string colourise( const std::string& str, const Colour& col ) { return Form( "\033[%d%s\033[0m", (int)col, str.c_str() ); }
/// Replace all occurences of a text by another
size_t replace_all( std::string& str, const std::string& from, const std::string& to );
+ namespace utils
+ {
+ /// Add a trailing "s" when needed
+ inline const char* s( unsigned short num ) { return ( num > 1 ) ? "s" : ""; }
+ /// Helper to print a vector
+ template<class T> std::string repr( const std::vector<T>& vec, const std::string& sep = "," ) {
+ return std::accumulate( std::next( vec.begin() ), vec.end(),
+ std::to_string( *vec.begin() ), [&sep]( std::string str, T xv ) {
+ return std::move( str )+sep+std::to_string( xv );
+ } );
+ }
+ }
}
/// Provide a random number generated along a uniform distribution between 0 and 1
#define drand() static_cast<double>( rand()/RAND_MAX )
#endif