Commit 694e3f8f authored by Lynn Garren's avatar Lynn Garren
Browse files

use std::shared_ptr and add MixMax Random engine

parent 43d0774c
......@@ -7,6 +7,7 @@ set( pkginclude_HEADERS DoubConv.hh
Hurd160Engine.h
Hurd288Engine.h
JamesRandom.h
MixMaxRng.h
MTwistEngine.h
NonRandomEngine.h
RandBinomial.h
......@@ -59,7 +60,9 @@ set( pkginclude_HEADERS DoubConv.hh
SeedTable.h
Stat.h
StaticRandomStates.h
TripleRand.h )
TripleRand.h
mixmax.h
mixmax_skip_N256.icc )
# notice that defs.h is not referenced here
INSTALL (FILES ${pkginclude_HEADERS}
......
// $Id:$
// -*- C++ -*-
//
// -----------------------------------------------------------------------
// HEP Random
// --- MixMaxRng ---
// class header file
// -----------------------------------------------------------------------
//
// This file interfaces the PseudoRandom Number Generator
// proposed by:
// N.Z. Akopov, G.K.Saviddy & N.G. Ter-Arutyunian
// "Matrix Generator of Pseudorandom Numbers",
// J.Compt.Phy. 97, 573 (1991)
// Preprint: EPI-867(18)-86, Yerevan June 1986.
// G. Savvidy & N. Savvidy
// "On the Monte Carlo Simulation of Physical Systems",
// J.Comput.Phys. 97 (1991) 566
// =======================================================================
// Implementation by Konstantin Savvidy - 2004-2015
// Release 0.99 and later: released under the LGPL license version 3.0
// =======================================================================
// CLHEP interface implemented by
// J. Apostolakis, G. Cosmo & K. Savvidy - Created: 6th July 2015
// CLHEP interface released under the LGPL license version 3.0
// =======================================================================
#ifndef MixMaxRng_h
#define MixMaxRng_h 1
#include "CLHEP/Random/defs.h"
#include "CLHEP/Random/RandomEngine.h"
#include "CLHEP/Random/mixmax.h"
namespace CLHEP {
/**
* @author K. Savvidy
* @ingroup random
*/
class MixMaxRng: public HepRandomEngine {
public:
MixMaxRng(std::istream& is);
MixMaxRng();
MixMaxRng(long seed);
MixMaxRng(int rowIndex, int colIndex);
virtual ~MixMaxRng();
// Constructor and destructor.
MixMaxRng(const MixMaxRng& rng);
MixMaxRng& operator=(const MixMaxRng& rng);
// Copy constructor and assignment operator.
double flat();
// Returns a pseudo random number between 0 and 1
// (excluding the zero: in (0,1] )
// smallest number which it will give is approximately 10^-19
void flatArray (const int size, double* vect);
// Fills the array "vect" of specified size with flat random values.
void setSeed(long seed, int dum=0);
// Sets the state of the algorithm according to seed.
void setSeeds(const long * seeds, int seedNum=0);
// Sets the initial state of the engine according to the array of between one and four 32-bit seeds.
// If the size of long is greater on the platform, only the lower 32-bits are used.
// Streams created from seeds differing by at least one bit somewhere are guaranteed absolutely
// to be independent and non-colliding for at least the next 10^100 random numbers
void saveStatus( const char filename[] = "MixMaxRngState.conf" ) const;
// Saves the the current engine state in the file given, by default MixMaxRngState.conf
void restoreStatus( const char filename[] = "MixMaxRngState.conf" );
// Reads a valid engine state from a given file, by default MixMaxRngState.conf
// and restores it.
void showStatus() const;
// Dumps the engine status on the screen.
operator unsigned int();
// 32-bit flat
virtual std::ostream & put (std::ostream & os) const;
virtual std::istream & get (std::istream & is);
static std::string beginTag ( );
virtual std::istream & getState ( std::istream & is );
std::string name() const;
static std::string engineName() {return "MixMaxRng";}
std::vector<unsigned long> put () const;
bool get (const std::vector<unsigned long> & v);
bool getState (const std::vector<unsigned long> & v);
static const unsigned int VECTOR_STATE_SIZE = 516; // 2N+4 for MIXMAX
private:
// Members defining the current status of the generator.
rng_state_st* fRngState;
};
} // namespace CLHEP
#endif
......@@ -103,7 +103,7 @@ private:
static double genBinomial( HepRandomEngine *anEngine, long n, double p );
shared_ptr<HepRandomEngine> localEngine;
std::shared_ptr<HepRandomEngine> localEngine;
long defaultN;
double defaultP;
......
......@@ -134,7 +134,7 @@ public:
private:
shared_ptr<HepRandomEngine> localEngine;
std::shared_ptr<HepRandomEngine> localEngine;
double defaultA;
double defaultB;
......
......@@ -101,7 +101,7 @@ private:
static double genChiSquare( HepRandomEngine *anEngine, double a );
shared_ptr<HepRandomEngine> localEngine;
std::shared_ptr<HepRandomEngine> localEngine;
double defaultA;
};
......
/*
Code adapted from:
http://www.jstatsoft.org/v05/i08/
Original disclaimer:
The ziggurat method for RNOR and REXP
Combine the code below with the main program in which you want
normal or exponential variates. Then use of RNOR in any expression
will provide a standard normal variate with mean zero, variance 1,
while use of REXP in any expression will provide an exponential variate
with density exp(-x),x>0.
Before using RNOR or REXP in your main, insert a command such as
zigset(86947731 );
with your own choice of seed value>0, rather than 86947731.
(If you do not invoke zigset(...) you will get all zeros for RNOR and REXP.)
For details of the method, see Marsaglia and Tsang, "The ziggurat method
for generating random variables", Journ. Statistical Software.
*/
// $Id:$
// -*- C++ -*-
//
// -----------------------------------------------------------------------
// HEP Random
// --- RandExpZiggurat ---
// class header file
// -----------------------------------------------------------------------
// This file is part of Geant4 (simulation toolkit for HEP).
// Class defining methods for shooting or firing Landau distributed
// random values.
// =======================================================================
//
// Code adapted from:
// http://www.jstatsoft.org/v05/i08/
//
//
// Original disclaimer:
// The ziggurat method for RNOR and REXP
// Combine the code below with the main program in which you want
// normal or exponential variates. Then use of RNOR in any expression
// will provide a standard normal variate with mean zero, variance 1,
// while use of REXP in any expression will provide an exponential variate
// with density exp(-x),x>0.
// Before using RNOR or REXP in your main, insert a command such as
// zigset(86947731 );
// with your own choice of seed value>0, rather than 86947731.
// (If you do not invoke zigset(...) you will get all zeros for RNOR and REXP.)
// For details of the method, see Marsaglia and Tsang, "The ziggurat method
// for generating random variables", Journ. Statistical Software.
// =======================================================================
#ifndef RandExpZiggurat_h
#define RandExpZiggurat_h 1
#include "CLHEP/Random/defs.h"
#include "CLHEP/Random/Random.h"
#include "CLHEP/Utility/memory.h"
#include "CLHEP/Utility/thread_local.h"
namespace CLHEP {
......@@ -83,11 +98,11 @@ public:
// the static generator.
inline float fire() {return fire(defaultMean);};
inline float fire( float mean ) {return ziggurat_REXP(localEngine)*mean;};
inline float fire( float mean ) {return ziggurat_REXP(localEngine.get())*mean;};
/* ENGINE IS INTRINSIC FLOAT
inline double fire() {return fire(defaultMean);};
inline double fire( double mean ) {return ziggurat_REXP(localEngine)*mean;};
inline double fire( double mean ) {return ziggurat_REXP(localEngine.get())*mean;};
*/
void fireArray ( const int size, float* vect );
......@@ -147,8 +162,7 @@ private:
// Private copy constructor. Defining it here disallows use.
RandExpZiggurat(const RandExpZiggurat& d);
HepRandomEngine* localEngine;
bool deleteEngine;
std::shared_ptr<HepRandomEngine> localEngine;
double defaultMean;
};
......@@ -161,11 +175,11 @@ using namespace CLHEP;
namespace CLHEP {
inline RandExpZiggurat::RandExpZiggurat(HepRandomEngine & anEngine, double mean ) : localEngine(&anEngine), deleteEngine(false), defaultMean(mean)
inline RandExpZiggurat::RandExpZiggurat(HepRandomEngine & anEngine, double mean ) : localEngine(&anEngine, do_nothing_deleter()), defaultMean(mean)
{
}
inline RandExpZiggurat::RandExpZiggurat(HepRandomEngine * anEngine, double mean ) : localEngine(anEngine), deleteEngine(true), defaultMean(mean)
inline RandExpZiggurat::RandExpZiggurat(HepRandomEngine * anEngine, double mean ) : localEngine(anEngine), defaultMean(mean)
{
}
......
......@@ -96,7 +96,7 @@ public:
private:
shared_ptr<HepRandomEngine> localEngine;
std::shared_ptr<HepRandomEngine> localEngine;
double defaultMean;
};
......
......@@ -198,7 +198,7 @@ private:
static CLHEP_THREAD_LOCAL unsigned long staticRandomInt;
static CLHEP_THREAD_LOCAL unsigned long staticFirstUnusedBit;
shared_ptr<HepRandomEngine> localEngine;
std::shared_ptr<HepRandomEngine> localEngine;
double defaultWidth;
double defaultA;
double defaultB;
......
......@@ -106,7 +106,7 @@ private:
static double genGamma( HepRandomEngine *anEngine, double k,
double lambda );
shared_ptr<HepRandomEngine> localEngine;
std::shared_ptr<HepRandomEngine> localEngine;
double defaultK;
double defaultLambda;
......
......@@ -153,7 +153,7 @@ protected:
double defaultMean;
double defaultStdDev;
shared_ptr<HepRandomEngine> localEngine;
std::shared_ptr<HepRandomEngine> localEngine;
private:
......
......@@ -129,7 +129,7 @@ public:
private:
shared_ptr<HepRandomEngine> localEngine;
std::shared_ptr<HepRandomEngine> localEngine;
std::vector<double> theIntegralPdf;
int nBins;
double oneOverNbins;
......
......@@ -108,7 +108,7 @@ protected:
private:
shared_ptr<HepRandomEngine> localEngine;
std::shared_ptr<HepRandomEngine> localEngine;
};
......
......@@ -116,7 +116,7 @@ protected:
private:
shared_ptr<HepRandomEngine> localEngine;
std::shared_ptr<HepRandomEngine> localEngine;
double status[3], oldm;
// static data
......
......@@ -104,7 +104,7 @@ public:
private:
shared_ptr<HepRandomEngine> localEngine;
std::shared_ptr<HepRandomEngine> localEngine;
double defaultA;
};
......
......@@ -31,6 +31,7 @@
#include "CLHEP/Random/Hurd160Engine.h"
#include "CLHEP/Random/Hurd288Engine.h"
#include "CLHEP/Random/JamesRandom.h"
#include "CLHEP/Random/MixMaxRng.h"
#include "CLHEP/Random/MTwistEngine.h"
#include "CLHEP/Random/RandEngine.h"
#include "CLHEP/Random/RanecuEngine.h"
......
// $Id:$
// -*- C++ -*-
//
// -----------------------------------------------------------------------
// MixMax Matrix PseudoRandom Number Generator
// --- MixMax ---
// class header file
// -----------------------------------------------------------------------
//
//
// Created by Konstantin Savvidy on Sun Feb 22 2004.
// The code is released under
// GNU Lesser General Public License v3
//
// Generator described in
// N.Z.Akopov, G.K.Savvidy and N.G.Ter-Arutyunian, Matrix Generator of Pseudorandom Numbers,
// J.Comput.Phys. 97, 573 (1991);
// Preprint EPI-867(18)-86, Yerevan Jun.1986;
//
// and
//
// K.Savvidy
// The MIXMAX random number generator
// Comp. Phys. Commun. (2015)
// http://dx.doi.org/10.1016/j.cpc.2015.06.003
//
// -----------------------------------------------------------------------
#ifndef CLHEP_MIXMAX_H_
#define CLHEP_MIXMAX_H_ 1
#include <stdio.h>
#include <stdint.h>
#define USE_INLINE_ASM YES
namespace CLHEP {
#ifdef __cplusplus
extern "C" {
#endif
const int N = 256;
/* The currently recommended N are 3150, 1260, 508, 256, 240, 88
Since the algorithm is linear in N, the cost per number is almost independent of N.
*/
#ifndef __LP64__
typedef uint64_t myuint;
//#warning but no problem, 'myuint' is 'uint64_t'
#else
typedef unsigned long long int myuint;
//#warning but no problem, 'myuint' is 'unsigned long long int'
#endif
struct rng_state_st
{
myuint V[N];
myuint sumtot;
int counter;
FILE* fh;
};
typedef struct rng_state_st rng_state_t; // C struct alias
int rng_get_N(void); // get the N programmatically, useful for checking the value for which the library was compiled
rng_state_t *rng_alloc(); /* allocate the state */
int rng_free(rng_state_t* X); /* free memory occupied by the state */
rng_state_t *rng_copy(myuint *Y); /* init from vector, takes the vector Y,
returns pointer to the newly allocated and initialized state */
void read_state(rng_state_t* X, const char filename[] );
void print_state(rng_state_t* X);
int iterate(rng_state_t* X);
myuint iterate_raw_vec(myuint* Y, myuint sumtotOld);
// FUNCTIONS FOR SEEDING
typedef uint32_t myID_t;
void seed_uniquestream(rng_state_t* X, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID );
/*
best choice: will make a state vector from which you can get at least 10^100 numbers
guaranteed mathematically to be non-colliding with any other stream prepared from another set of 32bit IDs,
so long as it is different by at least one bit in at least one of the four IDs
-- useful if you are running a parallel simulation with many clusters, many CPUs each
*/
void seed_spbox(rng_state_t* X, myuint seed); // non-linear method, makes certified unique vectors, probability for streams to collide is < 1/10^4600
void seed_vielbein(rng_state_t* X, unsigned int i); // seeds with the i-th unit vector, i = 0..N-1, for testing only
// FUNCTIONS FOR GETTING RANDOM NUMBERS
#ifdef __MIXMAX_C
myuint get_next(rng_state_t* X); // returns 64-bit int, which is between 1 and 2^61-1 inclusive
double get_next_float(rng_state_t* X); // returns double precision floating point number in (0,1]
#endif //__MIXMAX_C
void fill_array(rng_state_t* X, unsigned int n, double *array); // fastest method: set n to a multiple of N (e.g. n=256)
void iterate_and_fill_array(rng_state_t* X, double *array); // fills the array with N numbers
myuint precalc(rng_state_t* X);
/* needed if the state has been changed by something other than iterate, but no worries, seeding functions call this for you when necessary */
myuint apply_bigskip(myuint* Vout, myuint* Vin, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID );
// applies a skip of some number of steps calculated from the four IDs
void branch_inplace( rng_state_t* Xin, myID_t* ID ); // almost the same as apply_bigskip, but in-place and from a vector of IDs
#define BITS 61
/* magic with Mersenne Numbers */
#define M61 2305843009213693951ULL
myuint modadd(myuint foo, myuint bar);
myuint modmulM61(myuint s, myuint a);
myuint fmodmulM61(myuint cum, myuint s, myuint a);
#define MERSBASE M61 //xSUFF(M61)
#define MOD_PAYNE(k) ((((k)) & MERSBASE) + (((k)) >> BITS) ) // slightly faster than my old way, ok for addition
#define MOD_REM(k) ((k) % MERSBASE ) // latest Intel CPU is supposed to do this in one CPU cycle, but on my machines it seems to be 20% slower than the best tricks
#define MOD_MERSENNE(k) MOD_PAYNE(k)
#define INV_MERSBASE (0.43368086899420177360298E-18L)
// the charpoly is irreducible for the combinations of N and SPECIAL and has maximal period for N=508, 256, half period for 1260, and 1/12 period for 3150
// #if (N==256)
#define SPECIALMUL 0
#define SPECIAL 487013230256099064ULL // s=487013230256099064, m=1 -- good old MIXMAX
#define MOD_MULSPEC(k) fmodmulM61( 0, SPECIAL , (k) );
/*
#elif (N==17)
#define SPECIALMUL 36 // m=2^37+1
#elif (N==8)
#define SPECIALMUL 53 // m=2^53+1
#elif (N==40)
#define SPECIALMUL 42 // m=2^42+1
#elif (N==96)
#define SPECIALMUL 55 // m=2^55+1
#elif (N==64)
#define SPECIALMUL 55 // m=2^55 (!!!) and m=2^37+2
#elif (N==120)
#define SPECIALMUL 51 // m=2^51+1 and a SPECIAL=+1 (!!!)
#define SPECIAL 1
#define MOD_MULSPEC(k) (k);
#else
#warning Not a verified N, you are on your own!
#define SPECIALMUL 58
#endif // list of interesting N for modulus M61 ends here
*/
#ifndef __MIXMAX_C // c++ can put code into header files, why cant we? (with the inline declaration, should be safe from duplicate-symbol error)
#define get_next(X) GET_BY_MACRO(X)
#define get_next_float(X) get_next_float_BY_MACRO(X)
#endif // __MIXMAX_C
inline myuint GET_BY_MACRO(rng_state_t* X) {
int i;
i=X->counter;
if (i<=(N-1) ){
X->counter++;
return X->V[i];
}else{
X->sumtot = iterate_raw_vec(X->V, X->sumtot);
X->counter=2;
return X->V[1];
}
}
inline double get_next_float_BY_MACRO(rng_state_t* X){
int64_t Z=(int64_t)get_next(X);
#if defined(__x86_64__) && defined(__SSE__) && defined(__AVX__) && defined(USE_INLINE_ASM)
double F;
__asm__ __volatile__( "pxor %0, %0;"
//"cvtsi2sdq %1, %0;"
:"=x"(F)
//:"r"(Z)
);
F=Z;
return F*INV_MERSBASE;
#else
return Z*INV_MERSBASE;
#endif
}
// ERROR CODES - exit() is called with these
#define ARRAY_INDEX_OUT_OF_BOUNDS 0xFF01
#define SEED_WAS_ZERO 0xFF02
#define ERROR_READING_STATE_FILE 0xFF03
#define ERROR_READING_STATE_COUNTER 0xFF04
#define ERROR_READING_STATE_CHECKSUM 0xFF05
#ifdef __cplusplus
}
#endif
//#define HOOKUP_GSL 1
#ifdef HOOKUP_GSL // if you need to use mixmax through GSL, pass -DHOOKUP_GSL=1 to the compiler
#include <gsl/gsl_rng.h>
unsigned long gsl_get_next(void *vstate);
double gsl_get_next_float(void *vstate);
void seed_for_gsl(void *vstate, unsigned long seed);
static const gsl_rng_type mixmax_type =
{"MIXMAX", /* name */
MERSBASE, /* RAND_MAX */
1, /* RAND_MIN */
sizeof (rng_state_t),
&seed_for_gsl,
&gsl_get_next,
&gsl_get_next_float
};
unsigned long gsl_get_next(void *vstate) {
rng_state_t* X= (rng_state_t*)vstate;
return (unsigned long)get_next(X);
}
double gsl_get_next_float(void *vstate) {
rng_state_t* X= (rng_state_t*)vstate;
return ( (double)get_next(X)) * INV_MERSBASE;
}
void seed_for_gsl(void *vstate, unsigned long seed){
rng_state_t* X= (rng_state_t*)vstate;
seed_spbox(X,(myuint)seed);
}
const gsl_rng_type *gsl_rng_ran3 = &mixmax_type;
#endif // HOOKUP_GSL
} // namespace CLHEP
#endif // closing CLHEP_MIXMAX_H_
This diff is collapsed.
......@@ -7,9 +7,11 @@ clhep_build_library( Random DoubConv.cc
erfQ.cc
flatToGaussian.cc
gammln.cc
mixmax.cc
Hurd160Engine.cc
Hurd288Engine.cc
JamesRandom.cc
MixMaxRng.cc
MTwistEngine.cc
NonRandomEngine.cc
RandBinomial.cc
......