egs_rndm.h

Go to the documentation of this file.

/*
###############################################################################
#
#  EGSnrc egs++ random number generators headers
#  Copyright (C) 2015 National Research Council Canada
#
#  This file is part of EGSnrc.
#
#  EGSnrc is free software: you can redistribute it and/or modify it under
#  the terms of the GNU Affero General Public License as published by the
#  Free Software Foundation, either version 3 of the License, or (at your
#  option) any later version.
#
#  EGSnrc is distributed in the hope that it will be useful, but WITHOUT ANY
#  WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
#  FOR A PARTICULAR PURPOSE.  See the GNU Affero General Public License for
#  more details.
#
#  You should have received a copy of the GNU Affero General Public License
#  along with EGSnrc. If not, see <http://www.gnu.org/licenses/>.
#
###############################################################################
#
#  Author:          Iwan Kawrakow, 2005
#
#  Contributors:    Reid Townson
#
###############################################################################
*/
 
 
#ifndef EGS_RANDOM_GENERATOR_
#define EGS_RANDOM_GENERATOR_
 
#include "egs_libconfig.h"
#include "egs_math.h"
#include "egs_functions.h"
#include "egs_input_struct.h"
 
class EGS_Input;
 
inline void addRngDefinitionBlockImpl(EGS_BlockInput &rngBlock) {
    rngBlock.addSingleInput("type", true, "Algorithm used for RNG. Defaults to ranmar", {"ranmar"});
    rngBlock.addSingleInput("initial seeds", true, "Two integers that represent the inital seed. The same input file running on the same computer will give idential results for the same seeds, and statistically independent for different seeds. Parallel runs automatically vary the seeds for each job.");
    rngBlock.addSingleInput("high resolution", false, "Defaults to no. Microscale simulations may require high resolution RNG.", {"No", "Yes"});
}
 
inline void addRngDefinitionBlock(const std::shared_ptr<EGS_InputStruct> &blockPtr) {
    auto rngBlock = blockPtr->addBlockInput("rng definition");
    addRngDefinitionBlockImpl(*rngBlock);
}
 
inline void addRngDefinitionBlock(const std::shared_ptr<EGS_BlockInput> &blockPtr) {
    auto rngBlock = blockPtr->addBlockInput("rng definition");
    addRngDefinitionBlockImpl(*blockPtr);
}
 
inline string addRngDefinitionExample() {
    string example = {
        R"(
:start rng definition:
    type = ranmar
    initial seeds = 33 97
    high resolution = no
:stop rng definition:
)"};
    return example;
}
 
class EGS_EXPORT EGS_RandomGenerator {
 
public:
 
    EGS_RandomGenerator(int n=128);
 
    EGS_RandomGenerator(const EGS_RandomGenerator &r) : np(0) {
        copyBaseState(r);
    };
 
    virtual ~EGS_RandomGenerator() {
        delete [] rarray;
    };
 
    inline EGS_Float getUniform() {
        if (ip >= np) {
            fillArray(np,rarray);
            ip = 0;
        }
        return rarray[ip++];
    };
 
    EGS_I64 numbersGenerated() const {
        return count;
    };
 
    EGS_I64 numbersUsed() const {
        return ip<np ? count - np + ip : count;
    };
 
    inline void getAzimuth(EGS_Float &cphi, EGS_Float &sphi) {
#ifndef FAST_SINCOS
        EGS_Float xphi,xphi2,yphi,yphi2,rhophi;
        do {
            xphi = 2*getUniform() - 1;
            xphi2 = xphi*xphi;
            yphi = getUniform();
            yphi2 = yphi*yphi;
            rhophi = xphi2 + yphi2;
        }
        while (rhophi > 1);
        cphi = (xphi2 - yphi2)/rhophi;
        sphi = 2*xphi*yphi/rhophi;
#else
        EGS_Float phi = 2*M_PI*getUniform();
        cphi = cos(phi);
        sphi = sin(phi);
#endif
    };
 
    inline EGS_Float getGaussian() {
        if (have_x) {
            have_x = false;
            return the_x;
        }
        EGS_Float r = sqrt(-2*log(1-getUniform()));
        EGS_Float cphi, sphi;
        getAzimuth(cphi,sphi);
        have_x = true;
        the_x = r*sphi;
        return r*cphi;
    };
 
    static EGS_RandomGenerator *createRNG(EGS_Input *inp, int sequence=0);
 
    static EGS_RandomGenerator *defaultRNG(int sequence=0);
 
    virtual void fillArray(int n, EGS_Float *array) = 0;
 
 
    bool storeState(ostream &data);
    bool setState(istream &data);
    bool addState(istream &data);
    void resetCounter() {
        count = 0;
    };
 
    virtual EGS_RandomGenerator *getCopy() = 0;
 
    virtual void setState(EGS_RandomGenerator *r) = 0;
 
    virtual void saveState() = 0;
 
    virtual void resetState() = 0;
 
    virtual int  rngSize() const = 0;
 
    virtual void describeRNG() const {};
 
protected:
 
    EGS_I64   count;  
    int       np;     
    int       ip;     
    EGS_Float *rarray;
 
    virtual bool storePrivateState(ostream &data) =  0;
 
    virtual bool setPrivateState(istream &data) =  0;
 
    void copyBaseState(const EGS_RandomGenerator &r);
 
    void allocate(int n);
 
    int  baseSize() const {
        return 2*sizeof(EGS_I32) + sizeof(EGS_I64) + sizeof(bool) +
               sizeof(EGS_Float) + np*sizeof(EGS_Float);
    };
 
 
private:
 
    bool      have_x; 
    EGS_Float the_x;  
 
};
 
#endif