From 1c559833c191def7013456eeedc89d512952d323 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Wed, 31 Jan 2024 12:42:05 +0100 Subject: [PATCH] Use NrnRandom123 everywhere --- src/gnu/CMakeLists.txt | 1 + src/gnu/NrnRandom123RNG.cpp | 32 +++++++ src/gnu/NrnRandom123RNG.hpp | 20 ++--- src/gnu/Rand.cpp | 5 +- src/gnu/Rand.hpp | 9 -- src/ivoc/ivocrand.cpp | 167 +++++++++--------------------------- 6 files changed, 79 insertions(+), 155 deletions(-) diff --git a/src/gnu/CMakeLists.txt b/src/gnu/CMakeLists.txt index e16ae608b9..bee75b180b 100644 --- a/src/gnu/CMakeLists.txt +++ b/src/gnu/CMakeLists.txt @@ -15,6 +15,7 @@ add_library( NegExp.cpp Normal.cpp nrnisaac.cpp + NrnRandom123RNG.cpp nrnran123.cpp Poisson.cpp Rand.cpp diff --git a/src/gnu/NrnRandom123RNG.cpp b/src/gnu/NrnRandom123RNG.cpp index e69de29bb2..dca4301bcf 100644 --- a/src/gnu/NrnRandom123RNG.cpp +++ b/src/gnu/NrnRandom123RNG.cpp @@ -0,0 +1,32 @@ +#include "NrnRandom123RNG.hpp" + +NrnRandom123::NrnRandom123(std::uint32_t id1, std::uint32_t id2, std::uint32_t id3) { + s_ = nrnran123_newstream3(id1, id2, id3); +} + +NrnRandom123::~NrnRandom123() { + nrnran123_deletestream(s_); +} + +std::uint32_t NrnRandom123::asLong() { + return nrnran123_ipick(s_); +} + +double NrnRandom123::asDouble() { + return nrnran123_dblpick(s_); +} + +void NrnRandom123::reset() { + nrnran123_setseq(s_, 0, 0); +} + +std::uint32_t NrnRandom123::get_seq() const { + std::uint32_t seq; + char which; + nrnran123_getseq(s_, &seq, &which); + return seq; +} + +void NrnRandom123::set_seq(std::uint32_t seq) { + nrnran123_setseq(s_, seq, 0); +} diff --git a/src/gnu/NrnRandom123RNG.hpp b/src/gnu/NrnRandom123RNG.hpp index 5e790c939c..d4484d45f8 100644 --- a/src/gnu/NrnRandom123RNG.hpp +++ b/src/gnu/NrnRandom123RNG.hpp @@ -9,20 +9,10 @@ class NrnRandom123: public RNG { public: NrnRandom123(std::uint32_t id1, std::uint32_t id2, std::uint32_t id3 = 0); ~NrnRandom123(); - std::uint32_t asLong() { - return nrnran123_ipick(s_); - } - double asDouble() { - return nrnran123_dblpick(s_); - } - void reset() { - nrnran123_setseq(s_, 0, 0); - } + std::uint32_t asLong(); + double asDouble(); + void reset(); + std::uint32_t get_seq() const; + void set_seq(std::uint32_t seq); nrnran123_State* s_; }; -NrnRandom123::NrnRandom123(std::uint32_t id1, std::uint32_t id2, std::uint32_t id3) { - s_ = nrnran123_newstream3(id1, id2, id3); -} -NrnRandom123::~NrnRandom123() { - nrnran123_deletestream(s_); -} diff --git a/src/gnu/Rand.cpp b/src/gnu/Rand.cpp index a3b526de5e..195a415a75 100644 --- a/src/gnu/Rand.cpp +++ b/src/gnu/Rand.cpp @@ -1,13 +1,12 @@ #include "Rand.hpp" -#include "ACG.h" +#include "NrnRandom123RNG.hpp" #include "Normal.h" Rand::Rand(unsigned long seed, int size, Object* obj) { // printf("Rand\n"); - gen = new ACG(seed, size); + gen = new NrnRandom123(seed, size); rand = new Normal(0., 1., gen); - type_ = 0; obj_ = obj; } diff --git a/src/gnu/Rand.hpp b/src/gnu/Rand.hpp index f5b47a6a61..010e435898 100644 --- a/src/gnu/Rand.hpp +++ b/src/gnu/Rand.hpp @@ -5,20 +5,11 @@ struct Object; -/* type_: - * 0: ACG - * 1: MLCG - * 2: MCellRan4 - * 3: Isaac64 - * 4: Random123 - */ class Rand { public: Rand(unsigned long seed = 0, int size = 55, Object* obj = nullptr); ~Rand(); RNG* gen; Random* rand; - int type_; // can do special things with some kinds of RNG - // double* looks like random variable that gets changed on every fadvance Object* obj_; }; diff --git a/src/ivoc/ivocrand.cpp b/src/ivoc/ivocrand.cpp index f250049299..9661c00650 100644 --- a/src/ivoc/ivocrand.cpp +++ b/src/ivoc/ivocrand.cpp @@ -17,8 +17,6 @@ #include "ocobserv.h" #include -#include -#include #include #include #include @@ -32,9 +30,7 @@ #include #include #include -#include #include -#include #if HAVE_IV #include "ivoc.h" @@ -115,21 +111,30 @@ static void r_destruct(void* r) { // syntax: // r.ACG([seed],[size]) -static double r_ACG(void* r) { +static double r_nrnran123(void* r) { Rand* x = (Rand*) r; - - unsigned long seed = 0; - int size = 55; - + uint32_t id1 = 0; + uint32_t id2 = 0; + uint32_t id3 = 0; if (ifarg(1)) - seed = long(*getarg(1)); + id1 = (uint32_t) (chkarg(1, 0., dmaxuint)); if (ifarg(2)) - size = int(chkarg(2, 7, 98)); - - x->rand->generator(new ACG(seed, size)); - x->type_ = 0; + id2 = (uint32_t) (chkarg(2, 0., dmaxuint)); + if (ifarg(3)) + id3 = (uint32_t) (chkarg(3, 0., dmaxuint)); + try { + NrnRandom123* r123 = new NrnRandom123(id1, id2, id3); + x->rand->generator(r123); + } catch (const std::bad_alloc& e) { + hoc_execerror("Bad allocation for 'NrnRandom123'", e.what()); + } delete x->gen; x->gen = x->rand->generator(); + return 0.; +} + +static double r_ACG(void* r) { + r_nrnran123(r); return 1.; } @@ -139,83 +144,31 @@ static double r_ACG(void* r) { // r.MLCG([seed1],[seed2]) static double r_MLCG(void* r) { - Rand* x = (Rand*) r; - - unsigned long seed1 = 0; - unsigned long seed2 = 0; - - if (ifarg(1)) - seed1 = long(*getarg(1)); - if (ifarg(2)) - seed2 = long(*getarg(2)); - - x->rand->generator(new MLCG(seed1, seed2)); - delete x->gen; - x->gen = x->rand->generator(); - x->type_ = 1; + r_nrnran123(r); return 1.; } static double r_MCellRan4(void* r) { - Rand* x = (Rand*) r; - - uint32_t seed1 = 0; - uint32_t ilow = 0; + return r_nrnran123(r); +} - if (ifarg(1)) - seed1 = (uint32_t) (chkarg(1, 0., dmaxuint)); - if (ifarg(2)) - ilow = (uint32_t) (chkarg(2, 0., dmaxuint)); - MCellRan4* mcr = new MCellRan4(seed1, ilow); - x->rand->generator(mcr); - delete x->gen; - x->gen = x->rand->generator(); - x->type_ = 2; - return (double) mcr->orig_; +static double r_Isaac64(void* r) { + return r_nrnran123(r); } + long nrn_get_random_sequence(Rand* r) { - assert(r->type_ == 2); - MCellRan4* mcr = (MCellRan4*) r->gen; - return mcr->ihigh_; + auto* gen = static_cast(r->gen); + return gen->get_seq(); } void nrn_set_random_sequence(Rand* r, long seq) { - assert(r->type_ == 2); - MCellRan4* mcr = (MCellRan4*) r->gen; - mcr->ihigh_ = seq; + auto* gen = static_cast(r->gen); + gen->set_seq(seq); } int nrn_random_isran123(Rand* r, uint32_t* id1, uint32_t* id2, uint32_t* id3) { - if (r->type_ == 4) { - NrnRandom123* nr = (NrnRandom123*) r->gen; - nrnran123_getids3(nr->s_, id1, id2, id3); - return 1; - } - return 0; -} - -static double r_nrnran123(void* r) { - Rand* x = (Rand*) r; - uint32_t id1 = 0; - uint32_t id2 = 0; - uint32_t id3 = 0; - if (ifarg(1)) - id1 = (uint32_t) (chkarg(1, 0., dmaxuint)); - if (ifarg(2)) - id2 = (uint32_t) (chkarg(2, 0., dmaxuint)); - if (ifarg(3)) - id3 = (uint32_t) (chkarg(3, 0., dmaxuint)); - try { - NrnRandom123* r123 = new NrnRandom123(id1, id2, id3); - x->rand->generator(r123); - } catch (const std::bad_alloc& e) { - hoc_execerror("Bad allocation for 'NrnRandom123'", e.what()); - } - delete x->gen; - x->gen = x->rand->generator(); - x->type_ = 4; - return 0.; + return 1; } static double r_ran123_globalindex(void* r) { @@ -228,71 +181,29 @@ static double r_ran123_globalindex(void* r) { static double r_sequence(void* r) { Rand* x = (Rand*) r; - if (x->type_ != 2 && x->type_ != 4) { - hoc_execerror( - "Random.seq() can only be used if the random generator was MCellRan4 or Random123", 0); - } - - if (x->type_ == 4) { - uint32_t seq; - char which; - if (ifarg(1)) { - double s = chkarg(1, 0., 17179869183.); /* 2^34 - 1 */ - seq = (uint32_t) (s / 4.); - which = char(s - seq * 4.); - NrnRandom123* nr = (NrnRandom123*) x->gen; - nrnran123_setseq(nr->s_, seq, which); - } - nrnran123_getseq(((NrnRandom123*) x->gen)->s_, &seq, &which); - return double(seq) * 4. + double(which); - } - - MCellRan4* mcr = (MCellRan4*) x->gen; + uint32_t seq; + char which; if (ifarg(1)) { - mcr->ihigh_ = (long) (*getarg(1)); + double s = chkarg(1, 0., 17179869183.); /* 2^34 - 1 */ + seq = (uint32_t) (s / 4.); + which = char(s - seq * 4.); + NrnRandom123* nr = (NrnRandom123*) x->gen; + nrnran123_setseq(nr->s_, seq, which); } - return (double) mcr->ihigh_; + nrnran123_getseq(((NrnRandom123*) x->gen)->s_, &seq, &which); + return double(seq) * 4. + double(which); } int nrn_random123_setseq(Rand* r, uint32_t seq, char which) { - if (r->type_ != 4) { - return 0; - } nrnran123_setseq(((NrnRandom123*) r->gen)->s_, seq, which); return 1; } int nrn_random123_getseq(Rand* r, uint32_t* seq, char* which) { - if (r->type_ != 4) { - return 0; - } nrnran123_getseq(((NrnRandom123*) r->gen)->s_, seq, which); return 1; } -static double r_Isaac64(void* r) { - Rand* x = (Rand*) r; - - uint32_t seed1 = 0; - - if (ifarg(1)) { - seed1 = static_cast(*getarg(1)); - } - - double seed{}; - try { - Isaac64* mcr = new Isaac64(seed1); - x->rand->generator(mcr); - delete x->gen; - x->gen = x->rand->generator(); - x->type_ = 3; - seed = mcr->seed(); - } catch (const std::bad_alloc& e) { - hoc_execerror("Bad allocation for Isaac64 generator", e.what()); - } - return seed; -} - // Pick again from the distribution last used // syntax: // r.repick()