Random Numbers

  1. Internal random numbers
  2. External random numbers
  3. MIXMAX random numbers
  4. The methods
  5. Random number debugging
This page describes the random-number generator in PYTHIA and how it can be replaced by an external one.

Internal random numbers

The Rndm class generates random numbers, using the Marsaglia-Zaman-Tsang algorithm [Mar90].

Random numbers R uniformly distributed in 0 < R < 1 are obtained with 

 
   Rndm::flat();
There are also methods to generate according to an exponential, to x * exp(-x), to a Gaussian, or picked among a set of possibilities, which make use of flat().

If the random number generator is not initialized before, it will be so the first time it is asked to generate a random number, and then with the default seed, 19780503. This means that, by default, all runs will use identically the same random number sequence. This is convenient for debugging purposes, but dangerous if you intend to run several "identical" jobs to boost statistics. You can initialize, or reinitialize, with your own choice of seed with a 

 
   Rndm::init(seed);
Here values 0 < seed < 900 000 000 gives so many different random number sequences, while seed = 0 will call the Stdlib time(0) function to provide a "random" seed, and seed < 0 will revert back to the default seed.

The Pythia class defines a flag and a mode, that allows the seed to be set in the Pythia::init call. That would be the standard way for a user to pick the random number sequence in a run.

External random numbers

RndmEngine is a base class for the external handling of random-number generation. The user-written derived class is called if a pointer to it has been handed in with the pythia.rndmEnginePtr() method. Since the default Marsaglia-Zaman-Tsang algorithm is quite good, chances are that any replacement would be a step down, but this may still be required by consistency with other program elements in big experimental frameworks.

There is only one pure virtual method in RndmEngine, to generate one random number flat in the range between 0 and 1: 

 
  virtual double flat() = 0;
Note that methods for initialization are not provided in the base class, in part since input parameters may be specific to the generator used, in part since initialization can as well be taken care of externally to the Pythia code.

An example illustrating how to run with an external random number generator is provided in main23.cc.

MIXMAX random numbers

The MIXMAX class of random number generators utilizes matrix-recursion based on Anosov-Kolmogorov C-K systems, with the ability to create a large number of statistically independent sequences of random numbers based on different initial seeds. This is particularly advantageous in creating statistically independent samples when running a large number of parallel jobs, each with a different initial seed. In the plugin header Pythia8Plugins/MixMax.h an implementation of a MIXMAX random number generator is provided [Sav91,Sav15], courtesy of Konstantin Savvidy, as well as a PYTHIA interface through the MixMaxRndm class. In this implementation a dimensionality of 17 is used, as this has been found to provide faster access to large numbers of independent sequences. A timing comparison between the external MIXMAX random number generator, and the default internal PYTHIA random number generator is provided in the example main23.cc. The MIXMAX random number generator is found to be comparable in speed to the default generator. The primary methods of the MixMaxRndm class are given here.

MixMaxRndm::MixMaxRndm(int seed0, int seed1, int seed2, int seed3)  
for the given four 32-bit seed numbers. The sequence of numbers produced from this set of seeds is guaranteed not to collide with another if at least one bit of the four seeds is different, and, less than 10^100 random numbers are thrown.

The methods

We here collect a more complete and formal overview of the Rndm class methods.

Rndm::Rndm()  
construct a random number generator, but does not initialize it.

Rndm::Rndm(int seed)  
construct a random number generator, and initialize it for the given seed number.

bool Rndm::rndmEnginePtr( RndmEngine* rndmPtr)  
pass in pointer for external random number generation.

void Rndm::init(int seed = 0)  
initialize, or reinitialize, the random number generator for the given seed number. Not necessary if the seed was already set in the constructor.

double Rndm::flat()  
generate next random number uniformly between 0 and 1.

double Rndm::exp()  
generate random numbers according to exp(-x).

double Rndm::xexp()  
generate random numbers according to x exp(-x).

double Rndm::gauss()  
generate random numbers according to exp(-x^2/2).

pair<double, double> Rndm::gauss2()  
generate a pair of random numbers according to exp( -(x^2 + y^2) / 2). Is faster than two calls to gauss().

pair<Vec4, Vec4> Rndm::phaseSpace2(double eCM, double m1, double m2)  
generate a pair of vectors according to the phase space distribution of two particles at the specified eCM and with the specified masses.

int Rndm::pick(const vector<double>& prob)  
pick one option among vector of (positive) probabilities.

void Rndm::shuffle(vector<T>& vec)  
randomly shuffle a vector of objects (type T, i.e. any type, templated method) according to the Fisher-Yates algorithm.

bool Rndm::dumpState(string fileName)  
save the current state of the random number generator to a binary file. This involves two integers and 100 double-precision numbers. Intended for debug purposes. Note that binary files may be platform-dependent and thus not transportable.

bool Rndm::readState(string fileName)  
set the state of the random number generator by reading in a binary file saved by the above command. Comments as above.

RndmState Rndm::getState()  
save the current state of the random number generator as a struct RndmState. This can then later be used within the same run, e.g. as input to another Pythia instance. It circumvents the intermediate file of dumpState, but cannot be saved for a later run.

void Rndm::setState(RndmState& state)  
set the state of the random number generator by reading in a struct RndmState saved by the above command. Comments as above.

virtual double RndmEngine::flat()  
if you want to construct an external random number generator (or generator interface) then you must implement this method in your class derived from the RndmEningen base class, to give a random number between 0 and 1.

Random number debugging

In some cases, when trying to determine where two different versions of the PYTHIA code might diverge, it can be useful to trace all random numbers that are used. In this way, it is possible to see exactly where a random sequence first diverges. An experimental random number number debugging tool is available for advanced debugging. This tool is only available when using the internal PYTHIA random number generation class, Rndm, and can be enabled by specifying the following during the configuration of PYTHIA. 
 
./configure --obj-common=-DRNGDEBUG
Then, when building any of the examples, the random number debugging will be enabled. The output, for example when running main01, will look something like the following. 
 
ParticleDecays::twoBody:flat        3.5433e-01 src/ParticleDecays.cc:546 
Rndm::exp:flat                      2.7883e-01 src/Basics.cc:96 
Rndm::exp:flat                      7.8534e-01 src/Basics.cc:96 
ParticleDataEntry::pickChannel:flat 8.4012e-01 src/ParticleData.cc:542
The printout is specified by the following format. 
 
Class::ClassMethod:RndmMethod RandomNumber Location:LineNumber
Here, Class is the class where the random number is being called, ClassMethod is the method being called, and RndmMethod is the method of Rndm that is being called. This is followed by RandomNumber which is the random number which was returned by the method, and then Location which is the source file where the call was made from, and LineNumber the line number of the call. In some methods, more than just a single number is returned, for example Rndm::phaseSpace2. In this case, all the elements of the returned value are printed (here, four elements for two vectors).

It is possible to control the output of the debugging with a number of statically defined variables. These variables can be modified in whatever code is being run.

bool Rndm::debugNow  
This can be used to switch the debug statements on and off. If set to false then the debug statements will not be printed. If set to true (the default) they will be printed. In this way, printing of the debugging can be turned on only after a certain event has passed, for example.

bool Rndm::debugLocation  
By default, set to true. If set to false the location of the random number generator call within the source code will no longer be printed. This is useful when running a difference on two sets of output.

bool Rndm::debugIndex  
By default, set to false. If set to true the number of RNG calls will be counted, and the index of the RNG call will be printed. This can be useful when comparing two different outputs where other print statements have been inserted in the code.

int Rndm::debugCall  
This integer tracks the number of random number generation calls, and can be used to reset the index.

int Rndm::debugPrecision  
Sets the precision of the random number values being printed. By default, this is 4 but could be set to 17, for example, to print the full precision of a double.

Because there are a large number of random number calls within PYTHIA, one may wish to filter the output. This can be done by inserting filter strings into the following static variables. These filter strings are applied to the Class::ClassMethod portion of the output, and if any of the defined sets are found, then the RNG call is printed (an OR operation).

set<string> Rndm::debugStarts  
Print the call if the class name and method starts with this string.

set<string> Rndm::debugEnds  
Print the call if the class name and method ends with this string.

set<string> Rndm::debugContains  
Print the call if the class name and method contains this string.

set<string> Rndm::debugContains  
Print the call if the class name and method exactly matches this string. For example, the following code only keeps calls if they begin with SimpleTimeShower or SimpleSpaceShower. 

 
Rndm::debugStarts.insert("SimpleTimeShower"); 
Rndm::debugStarts.insert("SimpleSpaceShower");

There are some limitations to the tool. The biggest is that this level of debugging is only possible when the Rndm method name is not used elsewhere in the PYTHIA code. For example, the exp method is used also as just the standard exponential. This means that whenever Rndm::exp is called, the location from where it was called can no longer be reported, only that this method was called, and what value it produced (this can be seen in the example above). In the PYTHIA code, this also currently applies to the pick and shuffle methods. This similarly applies to user code. If flat is used as a method name in user code, than compiling with the -DRNGDEBUG flag will result in the code failing to compile. This can be overcome by including the include/Pythia8/RngDebug.h header after any PYTHIA headers, but before any additional headers. As an example, main93 will fail to compile with RNG debugging on. However, the following modification in main93.cc resolves the issue. 

 
#include "Pythia8/Pythia.h" 
#include "Pythia8/HeavyIons.h" 
#include "Pythia8/RngDebug.h" 
#include "Pythia8Plugins/Pythia8Rivet.h"
Here, RngDebug.h has been included after the last PYTHIA header, but before additional headers.