The SKIRT project
advanced radiative transfer for astrophysics
Public Member Functions | Protected Member Functions | Private Types | Private Attributes | Friends | List of all members
Random Class Reference

#include <Random.hpp>

Inheritance diagram for Random:
Inheritance graph
[legend]

Public Member Functions

double cdfLinLin (const Array &xv, const Array &Pv)
 
double cdfLogLog (const Array &xv, const Array &pv, const Array &Pv)
 
Direction direction ()
 
Direction direction (Direction bfk, double costheta)
 
double expon ()
 
double exponCutoff (double xmax)
 
double gauss ()
 
Vec maxwell ()
 
void pop ()
 
Position position (const Box &box)
 
void push (int seed)
 
int seed () const
 
double uniform ()
 
- Public Member Functions inherited from SimulationItem
template<class T >
T * find (bool setup=true) const
 
template<class T >
T * interface (int levels=-999999, bool setup=true) const
 
virtual string itemName () const
 
void setup ()
 
string typeAndName () const
 
- Public Member Functions inherited from Item
 Item (const Item &)=delete
 
virtual ~Item ()
 
void addChild (Item *child)
 
const vector< Item * > & children () const
 
virtual void clearItemListProperty (const PropertyDef *property)
 
void destroyChild (Item *child)
 
virtual bool getBoolProperty (const PropertyDef *property) const
 
virtual vector< double > getDoubleListProperty (const PropertyDef *property) const
 
virtual double getDoubleProperty (const PropertyDef *property) const
 
virtual string getEnumProperty (const PropertyDef *property) const
 
virtual int getIntProperty (const PropertyDef *property) const
 
virtual vector< Item * > getItemListProperty (const PropertyDef *property) const
 
virtual ItemgetItemProperty (const PropertyDef *property) const
 
virtual string getStringProperty (const PropertyDef *property) const
 
int getUtilityProperty (string name) const
 
virtual void insertIntoItemListProperty (const PropertyDef *property, int index, Item *item)
 
Itemoperator= (const Item &)=delete
 
Itemparent () const
 
virtual void removeFromItemListProperty (const PropertyDef *property, int index)
 
virtual void setBoolProperty (const PropertyDef *property, bool value)
 
virtual void setDoubleListProperty (const PropertyDef *property, vector< double > value)
 
virtual void setDoubleProperty (const PropertyDef *property, double value)
 
virtual void setEnumProperty (const PropertyDef *property, string value)
 
virtual void setIntProperty (const PropertyDef *property, int value)
 
virtual void setItemProperty (const PropertyDef *property, Item *item)
 
virtual void setStringProperty (const PropertyDef *property, string value)
 
void setUtilityProperty (string name, int value)
 
virtual string type () const
 

Protected Member Functions

 Random ()
 
void setupSelfBefore () override
 
- Protected Member Functions inherited from SimulationItem
 SimulationItem ()
 
virtual bool offersInterface (const std::type_info &interfaceTypeInfo) const
 
virtual void setupSelfAfter ()
 
virtual void setupSelfBefore ()
 
- Protected Member Functions inherited from Item
 Item ()
 

Private Types

using BaseType = SimulationItem
 
using ItemType = Random
 

Private Attributes

int _seed
 

Friends

class ItemRegistry
 

Detailed Description

The Random class offers pseudo-random number generation facilities in a multi-threaded environment. At the core of its operation is the function to generate a uniform deviate, i.e. a pseudo-random number drawn from the uniform probability distribution over the unit interval. Additional functions allow drawing pseudo-random numbers from other frequently-used probability distributions.

To avoid the need for synchronization between multiple execution threads, each thread receives its own thread-local pseudo-random number generator instance. The parent thread in each process (i.e. the thread that calls the setup() function) receives a predictable generator. This generator is initialized with a fixed state depending only on the value of the user-configurable seed property, so that it delivers exactly the same pseudo-random sequence for every parent thread in every process and for every execution of the program. All other (child) threads receive an arbitrary generator. These generators are seeded with a truly random state obtained from the operating system, so that they deliver a unique and unpredictable pseudo-random sequence for each thread, and for each execution of the program.

This mechanism is intended to support the following use cases. In single thread/single process mode, repeated execution (with the same seed value) must produce identical results. This is important, for example, for automated testing. The user can force a different randomization by specifying another seed value. In multi-threading and/or multi-processing mode, for tasks employing parallelization, each of the threads in each of the processes requires a different random sequence. On the other hand, in this mode, for serialized tasks, the (single) thread employed in each of the processes must receive the same random sequence. This is important for functions that rely on randomness but need to produce the same result in every parallel process.

The recommended use of the Random class is to include a single instance in each simulation run-time hierarchy. The use cases discussed above can be achieved as follows. In single thread/single process mode, perform all tasks in the parent thread. In multi-threading and/or multi-processing mode, perform serial tasks in the parent thread of each process, and perform all parallized tasks in a child thread.

As an additional service, this class allows temporarily installing a predictable random number generator with a given seed in the current execution thread through the push() and pop() functions. This supports the use case where, usually during setup, the same pseudo-random sequence is required in multiple places.

All random number generators used in this class are based on the 64-bit Mersenne twister, which offers a sufficiently long period and acceptable spectral properties for most purposes.

Constructor & Destructor Documentation

◆ Random()

Random::Random ( )
inlineprotected

Default constructor for concrete Item subclass Random : "the default random generator" .

Member Function Documentation

◆ cdfLinLin()

double Random::cdfLinLin ( const Array xv,
const Array Pv 
)

This function generates a random number drawn from an arbitrary probability distribution \(p(x)\,{\text{d}}x\) with corresponding cumulative distribution function \(P(x)\). The function accepts a discretized version \(P_i\) of the cdf sampled at a set of \(N\) points \(x_i\). A uniform deviate \(\cal{X}\) is generated, and the equation \({\cal{X}}=P(x)\) is solved using linear interpolation (i.e. assuming piece-wise linear behavior of the cdf (and equivalently, of the underlying pdf).

◆ cdfLogLog()

double Random::cdfLogLog ( const Array xv,
const Array pv,
const Array Pv 
)

This function generates a random number drawn from an arbitrary probability distribution \(p(x)\,{\text{d}}x\) with corresponding cumulative distribution function \(P(x)\). The function accepts discretized versions \(p_i\) and \(P_i\) of the pdf and cdf sampled at a set of \(N\) points \(x_i\). A uniform deviate \(\cal{X}\) is generated, and the equation \({\cal{X}}=P(x)\) is solved using a interpolation mechanism that assumes that the pdf is linear in log-log space between any two grid points (equivalent to power-law behavior), as described below.

Consider the pdf values \(p_i\) and \(p_{i+1}\) at two consecutive grid points \(x_i\) and \(x_{i+1}\). Assuming power-law behavior, the pdf between these two grid points can be written as

\[ p(x) = p_i \left(\frac{x}{x_i}\right)^{\alpha_i}, \quad\mathrm{with}\; \alpha_i = \frac{\ln p_{i+1}/\ln p_i}{\ln x_{i+1}/\ln x_i} \]

With \(\mathcal{X}\) a random deviate for which \(x_i\leq\mathcal{X}\leq x_{i+1}\) happens to be true, we thus need to invert the relation

\[ \mathcal{X}-x_i = \int_{x_i}^x p(x')\,\mathrm{d}x' = \int_{x_i}^x p_i \left(\frac{x'}{x_i}\right)^{\alpha_i}\mathrm{d}x' = p_i x_i \;\mathrm{gln}\left(-\alpha_i, \frac{x}{x_i}\right) \]

which leads to

\[x = x_i \;\mathrm{gexp}\left(-\alpha_i, \frac{\mathcal{X}-x_i}{p_i x_i}\right)\]

where \(\mathrm{gln}(a,x)\) is the generalized logarithm and \(\mathrm{gexp}(a,x)\) the generalized exponential, defined in the description of respectively the SpecialFunctions::gln() and SpecialFunctions::gexp() functions.

◆ direction() [1/2]

Direction Random::direction ( )

This function generates a random direction on the unit sphere, i.e. a couple \((\theta,\phi)\) from the two-dimensional probability density

\[ p(\theta,\phi)\,d\theta\,d\phi = \left(\frac{\sin\theta}{2}\,d\theta\right) \left(\frac{1}{2\pi}\,d\varphi\right).\]

A random direction on the unit sphere can thus be constructed by taking two uniform deviates \({\cal{X}}_1\) and \({\cal{X}}_2\), and solving the two equations

\[ \begin{split} {\cal{X}}_1 &= \int_0^\theta \frac{\sin\theta'\,d\theta'}{2} \\ {\cal{X}}_2 &= \int_0^\varphi \frac{d\varphi'}{2\pi} \end{split} \]

for \(\theta\) and \(\varphi\). The solution is readily found,

\[ \begin{split} \theta &= \arccos\left(2{\cal{X}}_1-1\right) \\ \varphi &= 2\pi\,{\cal {X}}_2. \end{split} \]

Once these spherical coordinates are calculated, a Direction object can be constructed by calling the constructor Direction::Direction(double theta, double phi).

◆ direction() [2/2]

Direction Random::direction ( Direction  bfk,
double  costheta 
)

This function generates a new direction on the unit sphere deviating from a given original direction \(\bf{k}\) by a given polar angle \(\theta\) (specified through its cosine) and a uniformly random azimuth angle \(\phi\). The function can use an arbitrary reference point for the azimuth angle since it is distributed uniformly. We use equation (A37) of Bianchi et al. 1996 (ApJ 465,127), which uses the projection of the z-axis as the reference point.

◆ expon()

double Random::expon ( )

This function generates a random number from an exponential distribution function, defined by the probability distribution

\[ p(x)\,{\rm d}x = {\rm e}^{-x}\,{\rm d}x.\]

A simple inversion technique is used.

◆ exponCutoff()

double Random::exponCutoff ( double  xmax)

This function generates a random number from an exponential distribution function with a cut-off, defined by the probability distribution

\[ p(x)\,{\rm d}x \propto \begin{cases} \;{\rm e}^{-x}\,{\rm d}x &\qquad \text{for $0<x<x_{\text{max}}$,} \\ \;0 &\qquad \text{for $x>x_{\text{max}}$.} \end{cases} \]

A simple inversion technique is used.

◆ gauss()

double Random::gauss ( )

This function generates a random number from a Gaussian distribution function with mean 0 and standard deviation 1, i.e. defined by the probability distribution

\[ p(x)\,{\rm d}x = \frac{1}{\sqrt{2\pi}}\, {\rm e}^{-\frac12\,x^2}\,{\rm d}x.\]

The algorithm used and the implementation are taken from Press et al. (2002).

◆ maxwell()

Vec Random::maxwell ( )

This function generates a random velocity from a three-dimensional Maxwell-Boltzmann distribution with velocity dispersion 1, which is equivalent to a Gaussian distribution with mean 0 and standard deviation 1 for each of the Cartesian velocity components. The implementation simply calls the function gauss() to obtain each of these three velocity components.

◆ pop()

void Random::pop ( )

This function pops the most recently pushed random number generator from the stack for the current thread and establishes it as the active random number generator for the current thread. If the stack does not contain a random number generator, the behavior of this function is undefined.

◆ position()

Position Random::position ( const Box box)

This function generates a uniformly distributed random position in a given box (i.e. a cuboid lined up with the coordinate axes).

◆ push()

void Random::push ( int  seed)

This function pushes the active random number generator for the current thread onto the stack for the current thread and establishes a new predictable random number generator, initialized with the specified seed, as the active random number generator for the current thread.

◆ seed()

Random::seed ( ) const
inline

This function returns the value of the discoverable integer property seed : "the seed for the random generator" .

The minimum value for this property is "0" .

The maximum value for this property is "1000000" .

The default value for this property is given by the conditional value expression "0" .

This property is displayed only if the Boolean expression "Level3" evaluates to true after replacing the names by true or false depending on their presence.

◆ setupSelfBefore()

void Random::setupSelfBefore ( )
overrideprotectedvirtual

This function initializes the random generator for the calling thread (deemed the parent thread) to its fixed initial state depending on the value of the user-configurable seed property.

Reimplemented from SimulationItem.

◆ uniform()

double Random::uniform ( )

This function generates a uniform deviate, i.e. a random double precision number in the open interval (0,1). The interval borders zero and one are never returned.


The documentation for this class was generated from the following file: