From 1c559833c191def7013456eeedc89d512952d323 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Wed, 31 Jan 2024 12:42:05 +0100 Subject: [PATCH 01/12] 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() From 764070a22f6580a2f6a0dc8fe3967ff7c79c377f Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Wed, 31 Jan 2024 12:46:31 +0100 Subject: [PATCH 02/12] Removal of isaac --- cmake/NeuronFileLists.cmake | 1 - src/gnu/ACG.cpp | 292 ---------------------------------- src/gnu/ACG.h | 68 -------- src/gnu/CMakeLists.txt | 6 - src/gnu/Isaac64RNG.cpp | 20 --- src/gnu/Isaac64RNG.hpp | 33 ---- src/gnu/MCellRan4RNG.cpp | 15 -- src/gnu/MCellRan4RNG.hpp | 34 ---- src/gnu/MLCG.cpp | 110 ------------- src/gnu/MLCG.h | 87 ---------- src/gnu/isaac64.cpp | 145 ----------------- src/gnu/isaac64.h | 88 ---------- src/gnu/nrnisaac.cpp | 25 --- src/gnu/nrnisaac.h | 9 -- src/ivoc/ivocmain.cpp | 9 -- src/ivoc/ivocrand.cpp | 1 - src/nrniv/multisend_setup.cpp | 10 +- 17 files changed, 5 insertions(+), 948 deletions(-) delete mode 100755 src/gnu/ACG.cpp delete mode 100755 src/gnu/ACG.h delete mode 100644 src/gnu/Isaac64RNG.cpp delete mode 100644 src/gnu/Isaac64RNG.hpp delete mode 100644 src/gnu/MCellRan4RNG.cpp delete mode 100644 src/gnu/MCellRan4RNG.hpp delete mode 100755 src/gnu/MLCG.cpp delete mode 100755 src/gnu/MLCG.h delete mode 100644 src/gnu/isaac64.cpp delete mode 100644 src/gnu/isaac64.h delete mode 100644 src/gnu/nrnisaac.cpp delete mode 100644 src/gnu/nrnisaac.h diff --git a/cmake/NeuronFileLists.cmake b/cmake/NeuronFileLists.cmake index fb6ca521e7..84cd4786ac 100644 --- a/cmake/NeuronFileLists.cmake +++ b/cmake/NeuronFileLists.cmake @@ -10,7 +10,6 @@ set(STRUCTURED_HEADER_FILES_TO_INSTALL neuron/model_data_fwd.hpp) set(HEADER_FILES_TO_INSTALL gnu/mcran4.h - gnu/nrnisaac.h gnu/nrnran123.h nrniv/backtrace_utils.h nrniv/bbsavestate.h diff --git a/src/gnu/ACG.cpp b/src/gnu/ACG.cpp deleted file mode 100755 index 0f9136b347..0000000000 --- a/src/gnu/ACG.cpp +++ /dev/null @@ -1,292 +0,0 @@ -// This may look like C code, but it is really -*- C++ -*- -/* -Copyright (C) 1989 Free Software Foundation - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library 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 Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. -*/ - -#include -#ifdef __GNUG__ -#pragma implementation -#endif -#include "ACG.h" -#include - -// -// This is an extension of the older implementation of Algorithm M -// which I previously supplied. The main difference between this -// version and the old code are: -// -// + Andres searched high & low for good constants for -// the LCG. -// -// + theres more bit chopping going on. -// -// The following contains his comments. -// -// agn@UNH.CS.CMU.EDU sez.. -// -// The generator below is based on 2 well known -// methods: Linear Congruential (LCGs) and Additive -// Congruential generators (ACGs). -// -// The LCG produces the longest possible sequence -// of 32 bit random numbers, each being unique in -// that sequence (it has only 32 bits of state). -// It suffers from 2 problems: a) Independence -// isnt great, that is the (n+1)th number is -// somewhat related to the preceding one, unlike -// flipping a coin where knowing the past outcomes -// dont help to predict the next result. b) -// Taking parts of a LCG generated number can be -// quite non-random: for example, looking at only -// the least significant byte gives a permuted -// 8-bit counter (that has a period length of only -// 256). The advantage of an LCA is that it is -// perfectly uniform when run for the entire period -// length (and very uniform for smaller sequences -// too, if the parameters are chosen carefully). -// -// ACGs have extremly long period lengths and -// provide good independence. Unfortunately, -// uniformity isnt not too great. Furthermore, I -// didnt find any theoretically analysis of ACGs -// that addresses uniformity. -// -// The RNG given below will return numbers -// generated by an LCA that are permuted under -// control of a ACG. 2 permutations take place: the -// 4 bytes of one LCG generated number are -// subjected to one of 16 permutations selected by -// 4 bits of the ACG. The permutation a such that -// byte of the result may come from each byte of -// the LCG number. This effectively destroys the -// structure within a word. Finally, the sequence -// of such numbers is permuted within a range of -// 256 numbers. This greatly improves independence. -// -// -// Algorithm M as describes in Knuths "Art of Computer Programming", -// Vol 2. 1969 -// is used with a linear congruential generator (to get a good uniform -// distribution) that is permuted with a Fibonacci additive congruential -// generator to get good independence. -// -// Bit, byte, and word distributions were extensively tested and pass -// Chi-squared test near perfect scores (>7E8 numbers tested, Uniformity -// assumption holds with probability > 0.999) -// -// Run-up tests for on 7E8 numbers confirm independence with -// probability > 0.97. -// -// Plotting random points in 2d reveals no apparent structure. -// -// Autocorrelation on sequences of 5E5 numbers (A(i) = SUM X(n)*X(n-i), -// i=1..512) -// results in no obvious structure (A(i) ~ const). -// -// Except for speed and memory requirements, this generator outperforms -// random() for all tests. (random() scored rather low on uniformity tests, -// while independence test differences were less dramatic). -// -// AGN would like to.. -// thanks to M.Mauldin, H.Walker, J.Saxe and M.Molloy for inspiration & help. -// -// And I would (DGC) would like to thank Donald Kunth for AGN for letting me -// use his extensions in this implementation. -// - -// -// Part of the table on page 28 of Knuth, vol II. This allows us -// to adjust the size of the table at the expense of shorter sequences. -// - -static int randomStateTable[][3] = { -{3,7,16}, {4,9, 32}, {3,10, 32}, {1,11, 32}, {1,15,64}, {3,17,128}, -{7,18,128}, {3,20,128}, {2,21, 128}, {1,22, 128}, {5,23, 128}, {3,25, 128}, -{2,29, 128}, {3,31, 128}, {13,33, 256}, {2,35, 256}, {11,36, 256}, -{14,39,256}, {3,41,256}, {9,49,256}, {3,52,256}, {24,55,256}, {7,57, 256}, -{19,58,256}, {38,89,512}, {17,95,512}, {6,97,512}, {11,98,512}, {-1,-1,-1} }; - -// -// spatial permutation table -// RANDOM_PERM_SIZE must be a power of two -// - -#define RANDOM_PERM_SIZE 64 -uint32_t randomPermutations[RANDOM_PERM_SIZE] = { -0xffffffff, 0x00000000, 0x00000000, 0x00000000, // 3210 -0x0000ffff, 0x00ff0000, 0x00000000, 0xff000000, // 2310 -0xff0000ff, 0x0000ff00, 0x00000000, 0x00ff0000, // 3120 -0x00ff00ff, 0x00000000, 0xff00ff00, 0x00000000, // 1230 - -0xffff0000, 0x000000ff, 0x00000000, 0x0000ff00, // 3201 -0x00000000, 0x00ff00ff, 0x00000000, 0xff00ff00, // 2301 -0xff000000, 0x00000000, 0x000000ff, 0x00ffff00, // 3102 -0x00000000, 0x00000000, 0x00000000, 0xffffffff, // 2103 - -0xff00ff00, 0x00000000, 0x00ff00ff, 0x00000000, // 3012 -0x0000ff00, 0x00000000, 0x00ff0000, 0xff0000ff, // 2013 -0x00000000, 0x00000000, 0xffffffff, 0x00000000, // 1032 -0x00000000, 0x0000ff00, 0xffff0000, 0x000000ff, // 1023 - -0x00000000, 0xffffffff, 0x00000000, 0x00000000, // 0321 -0x00ffff00, 0xff000000, 0x00000000, 0x000000ff, // 0213 -0x00000000, 0xff000000, 0x0000ffff, 0x00ff0000, // 0132 -0x00000000, 0xff00ff00, 0x00000000, 0x00ff00ff // 0123 -}; - -// -// SEED_TABLE_SIZE must be a power of 2 -// -#define SEED_TABLE_SIZE 32 -static uint32_t seedTable[SEED_TABLE_SIZE] = { -0xbdcc47e5, 0x54aea45d, 0xec0df859, 0xda84637b, -0xc8c6cb4f, 0x35574b01, 0x28260b7d, 0x0d07fdbf, -0x9faaeeb0, 0x613dd169, 0x5ce2d818, 0x85b9e706, -0xab2469db, 0xda02b0dc, 0x45c60d6e, 0xffe49d10, -0x7224fea3, 0xf9684fc9, 0xfc7ee074, 0x326ce92a, -0x366d13b5, 0x17aaa731, 0xeb83a675, 0x7781cb32, -0x4ec7c92d, 0x7f187521, 0x2cf346b4, 0xad13310f, -0xb89cff2b, 0x12164de1, 0xa865168d, 0x32b56cdf -}; - -// -// The LCG used to scramble the ACG -// -// -// LC-parameter selection follows recommendations in -// "Handbook of Mathematical Functions" by Abramowitz & Stegun 10th, edi. -// -// LC_A = 251^2, ~= sqrt(2^32) = 66049 -// LC_C = result of a long trial & error series = 3907864577 -// - -static const uint32_t LC_A = 66049; -static const uint32_t LC_C = 3907864577U; -static inline uint32_t LCG(uint32_t x) -{ - return( x * LC_A + LC_C ); -} - - -ACG::ACG(uint32_t seed, int size) -{ - int l; - initialSeed = seed; - // - // Determine the size of the state table - // - - for (l = 0; - randomStateTable[l][0] != -1 && randomStateTable[l][1] < size; - l++); - - if (randomStateTable[l][1] == -1) { - l--; - } - - initialTableEntry = l; - - stateSize = randomStateTable[ initialTableEntry ][ 1 ]; - auxSize = randomStateTable[ initialTableEntry ][ 2 ]; - - // - // Allocate the state table & the auxillary table in a single malloc - // - - state = new uint32_t[stateSize + auxSize]; - auxState = &state[stateSize]; - - reset(); -} - -// -// Initialize the state -// -void -ACG::reset() -{ - uint32_t u; - - if (initialSeed < SEED_TABLE_SIZE) { - u = seedTable[ initialSeed ]; - } else { - u = initialSeed ^ seedTable[ initialSeed & (SEED_TABLE_SIZE-1) ]; - } - - - j = randomStateTable[ initialTableEntry ][ 0 ] - 1; - k = randomStateTable[ initialTableEntry ][ 1 ] - 1; - - int i; - for(i = 0; i < stateSize; i++) { - state[i] = u = LCG(u); - } - - for (i = 0; i < auxSize; i++) { - auxState[i] = u = LCG(u); - } - - k = u % stateSize; - int tailBehind = (stateSize - randomStateTable[ initialTableEntry ][ 0 ]); - j = k - tailBehind; - if (j < 0) { - j += stateSize; - } - - lcgRecurr = u; - - assert(sizeof(double) == 2 * sizeof(int32_t)); -} - -ACG::~ACG() -{ - if (state) delete [] state; - state = 0; - // don't delete auxState, it's really an alias for state. -} - -// -// Returns 32 bits of random information. -// - -uint32_t -ACG::asLong() -{ - uint32_t result = state[k] + state[j]; - state[k] = result; - j = (j <= 0) ? (stateSize-1) : (j-1); - k = (k <= 0) ? (stateSize-1) : (k-1); - - short int auxIndex = (result >> 24) & (auxSize - 1); - uint32_t auxACG = auxState[auxIndex]; - auxState[auxIndex] = lcgRecurr = LCG(lcgRecurr); - - // - // 3c is a magic number. We are doing four masks here, so we - // do not want to run off the end of the permutation table. - // This insures that we have always got four entries left. - // - uint32_t *perm = & randomPermutations[result & 0x3c]; - - result = *(perm++) & auxACG; - result |= *(perm++) & ((auxACG << 24) - | ((auxACG >> 8)& 0xffffff)); - result |= *(perm++) & ((auxACG << 16) - | ((auxACG >> 16) & 0xffff)); - result |= *(perm++) & ((auxACG << 8) - | ((auxACG >> 24) & 0xff)); - - return(result); -} diff --git a/src/gnu/ACG.h b/src/gnu/ACG.h deleted file mode 100755 index 63230e0031..0000000000 --- a/src/gnu/ACG.h +++ /dev/null @@ -1,68 +0,0 @@ -// This may look like C code, but it is really -*- C++ -*- -/* -Copyright (C) 1988 Free Software Foundation - written by Dirk Grunwald (grunwald@cs.uiuc.edu) - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library 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 Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ -#ifndef _ACG_h -#define _ACG_h 1 - -#include "RNG.h" -#include -#ifdef __GNUG__ -//#pragma interface -#endif - -// -// Additive number generator. This method is presented in Volume II -// of The Art of Computer Programming by Knuth. I've coded the algorithm -// and have added the extensions by Andres Nowatzyk of CMU to randomize -// the result of algorithm M a bit by using an LCG & a spatial -// permutation table. -// -// The version presented uses the same constants for the LCG that Andres -// uses (chosen by trial & error). The spatial permutation table is -// the same size (it's based on word size). This is for 32-bit words. -// -// The ``auxillary table'' used by the LCG table varies in size, and -// is chosen to be the the smallest power of two which is larger than -// twice the size of the state table. -// - -class ACG : public RNG { - - uint32_t initialSeed; // used to reset generator - int initialTableEntry; - - uint32_t *state; - uint32_t *auxState; - short stateSize; - short auxSize; - uint32_t lcgRecurr; - short j; - short k; - -protected: - -public: - ACG(uint32_t seed = 0, int size = 55); - virtual ~ACG(); - // - // Return a long-words word of random bits - // - virtual uint32_t asLong(); - virtual void reset(); -}; - -#endif diff --git a/src/gnu/CMakeLists.txt b/src/gnu/CMakeLists.txt index bee75b180b..1d8783a621 100644 --- a/src/gnu/CMakeLists.txt +++ b/src/gnu/CMakeLists.txt @@ -1,20 +1,14 @@ add_library( nrngnu STATIC - ACG.cpp Binomial.cpp DiscUnif.cpp Erlang.cpp Geom.cpp HypGeom.cpp - isaac64.cpp - Isaac64RNG.cpp LogNorm.cpp - MCellRan4RNG.cpp mcran4.cpp - MLCG.cpp NegExp.cpp Normal.cpp - nrnisaac.cpp NrnRandom123RNG.cpp nrnran123.cpp Poisson.cpp diff --git a/src/gnu/Isaac64RNG.cpp b/src/gnu/Isaac64RNG.cpp deleted file mode 100644 index 559c307d4d..0000000000 --- a/src/gnu/Isaac64RNG.cpp +++ /dev/null @@ -1,20 +0,0 @@ -#include "Isaac64RNG.hpp" - -uint32_t Isaac64::cnt_ = 0; - -Isaac64::Isaac64(std::uint32_t seed) { - if (cnt_ == 0) { - cnt_ = 0xffffffff; - } - --cnt_; - seed_ = seed; - if (seed_ == 0) { - seed_ = cnt_; - } - rng_ = nrnisaac_new(); - reset(); -} - -Isaac64::~Isaac64() { - nrnisaac_delete(rng_); -} diff --git a/src/gnu/Isaac64RNG.hpp b/src/gnu/Isaac64RNG.hpp deleted file mode 100644 index c825c741ca..0000000000 --- a/src/gnu/Isaac64RNG.hpp +++ /dev/null @@ -1,33 +0,0 @@ -#pragma once - -#include - -#include "RNG.h" -#include "nrnisaac.h" - -class Isaac64: public RNG { - public: - Isaac64(std::uint32_t seed = 0); - ~Isaac64(); - std::uint32_t asLong() { - return nrnisaac_uint32_pick(rng_); - } - void reset() { - nrnisaac_init(rng_, seed_); - } - double asDouble() { - return nrnisaac_dbl_pick(rng_); - } - std::uint32_t seed() { - return seed_; - } - void seed(std::uint32_t s) { - seed_ = s; - reset(); - } - - private: - std::uint32_t seed_; - void* rng_; - static std::uint32_t cnt_; -}; diff --git a/src/gnu/MCellRan4RNG.cpp b/src/gnu/MCellRan4RNG.cpp deleted file mode 100644 index 43b732cbe1..0000000000 --- a/src/gnu/MCellRan4RNG.cpp +++ /dev/null @@ -1,15 +0,0 @@ -#include "MCellRan4RNG.hpp" - -MCellRan4::MCellRan4(std::uint32_t ihigh, std::uint32_t ilow) { - ++cnt_; - ilow_ = ilow; - ihigh_ = ihigh; - if (ihigh_ == 0) { - ihigh_ = cnt_; - ihigh_ = (std::uint32_t) asLong(); - } - orig_ = ihigh_; -} -MCellRan4::~MCellRan4() {} - -std::uint32_t MCellRan4::cnt_ = 0; diff --git a/src/gnu/MCellRan4RNG.hpp b/src/gnu/MCellRan4RNG.hpp deleted file mode 100644 index a4356562a9..0000000000 --- a/src/gnu/MCellRan4RNG.hpp +++ /dev/null @@ -1,34 +0,0 @@ -#pragma once - -#include - -#include "RNG.h" -#include "mcran4.h" - -// The decision that has to be made is whether each generator instance -// should have its own seed or only one seed for all. We choose separate -// seed for each but if arg not present or 0 then seed chosen by system. - -// the addition of ilow > 0 means that value is used for the lowindex -// instead of the mcell_ran4_init global 32 bit lowindex. - -class MCellRan4: public RNG { - public: - MCellRan4(std::uint32_t ihigh = 0, std::uint32_t ilow = 0); - virtual ~MCellRan4(); - virtual std::uint32_t asLong() { - return (std::uint32_t) (ilow_ == 0 ? mcell_iran4(&ihigh_) : nrnRan4int(&ihigh_, ilow_)); - } - virtual void reset() { - ihigh_ = orig_; - } - virtual double asDouble() { - return (ilow_ == 0 ? mcell_ran4a(&ihigh_) : nrnRan4dbl(&ihigh_, ilow_)); - } - std::uint32_t ihigh_; - std::uint32_t orig_; - std::uint32_t ilow_; - - private: - static std::uint32_t cnt_; -}; diff --git a/src/gnu/MLCG.cpp b/src/gnu/MLCG.cpp deleted file mode 100755 index dce8358a43..0000000000 --- a/src/gnu/MLCG.cpp +++ /dev/null @@ -1,110 +0,0 @@ -// This may look like C code, but it is really -*- C++ -*- -/* -Copyright (C) 1989 Free Software Foundation - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library 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 Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ -#ifdef __GNUG__ -#pragma implementation -#endif -#include "MLCG.h" -// -// SEED_TABLE_SIZE must be a power of 2 -// - - -#define SEED_TABLE_SIZE 32 - -static int32_t seedTable[SEED_TABLE_SIZE] = { - /* 0xbdcc47e5 */ -1110685723, /* 0x54aea45d */ 1420731485, - /* 0xec0df859 */ -334628775, /* 0xda84637b */ -628857989, - /* 0xc8c6cb4f */ -926495921, /* 0x35574b01 */ 894913281, - /* 0x28260b7d */ 673581949, /* 0x0d07fdbf */ 218627519, - /* 0x9faaeeb0 */ -1616187728, /* 0x613dd169 */ 1631441257, - /* 0x5ce2d818 */ 1558370328, /* 0x85b9e706 */ -2051414266, - /* 0xab2469db */ -1423676965, /* 0xda02b0dc */ -637357860, - /* 0x45c60d6e */ 1170607470, /* 0xffe49d10 */ -1794800, - /* 0x7224fea3 */ 1915027107, /* 0xf9684fc9 */ -110604343, - /* 0xfc7ee074 */ -58793868, /* 0x326ce92a */ 845998378, - /* 0x366d13b5 */ 913118133, /* 0x17aaa731 */ 397059889, - /* 0xeb83a675 */ -343693707, /* 0x7781cb32 */ 2004994866, - /* 0x4ec7c92d */ 1321716013, /* 0x7f187521 */ 2132309281, - /* 0x2cf346b4 */ 754140852, /* 0xad13310f */ -1391251185, - /* 0xb89cff2b */ -1197670613, /* 0x12164de1 */ 303451617, - /* 0xa865168d */ -1469770099, /* 0x32b56cdf */ 850750687 -}; - -MLCG::MLCG(int32_t seed1, int32_t seed2) -{ - initialSeedOne = seed1; - initialSeedTwo = seed2; - reset(); -} - -void -MLCG::reset() -{ - int32_t seed1 = initialSeedOne; - int32_t seed2 = initialSeedTwo; - - // Most people pick stupid seed numbers that do not have enough - // bits. In this case, if they pick a small seed number, we - // map that to a specific seed. - - if (seed1 < 0) { - seed1 = (seed1 + 2147483561); - seed1 = (seed1 < 0) ? -seed1 : seed1; - } - - if (seed2 < 0) { - seed2 = (seed2 + 2147483561); - seed2 = (seed2 < 0) ? -seed2 : seed2; - } - - if (seed1 > -1 && seed1 < SEED_TABLE_SIZE) { - seedOne = seedTable[seed1]; - } else { - seedOne = seed1 ^ seedTable[seed1 & (SEED_TABLE_SIZE-1)]; - } - - if (seed2 > -1 && seed2 < SEED_TABLE_SIZE) { - seedTwo = seedTable[seed2]; - } else { - seedTwo = seed2 ^ seedTable[ seed2 & (SEED_TABLE_SIZE-1) ]; - } - seedOne = (seedOne % 2147483561) + 1; - seedTwo = (seedTwo % 2147483397) + 1; -} - -uint32_t MLCG::asLong() -{ - int32_t k = seedOne % 53668; - - seedOne = 40014 * (seedOne-k * 53668) - k * 12211; - if (seedOne < 0) { - seedOne += 2147483563; - } - - k = seedTwo % 52774; - seedTwo = 40692 * (seedTwo - k * 52774) - k * 3791; - if (seedTwo < 0) { - seedTwo += 2147483399; - } - - int32_t z = seedOne - seedTwo; - if (z < 1) { - z += 2147483562; - } - return( (unsigned long) z); -} - diff --git a/src/gnu/MLCG.h b/src/gnu/MLCG.h deleted file mode 100755 index 2b566acb91..0000000000 --- a/src/gnu/MLCG.h +++ /dev/null @@ -1,87 +0,0 @@ -// This may look like C code, but it is really -*- C++ -*- -/* -Copyright (C) 1988 Free Software Foundation - written by Dirk Grunwald (grunwald@cs.uiuc.edu) - -This file is part of the GNU C++ Library. This library is free -software; you can redistribute it and/or modify it under the terms of -the GNU Library General Public License as published by the Free -Software Foundation; either version 2 of the License, or (at your -option) any later version. This library 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 Library General Public License for more details. -You should have received a copy of the GNU Library General Public -License along with this library; if not, write to the Free Software -Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -*/ -#ifndef _MLCG_h -#define _MLCG_h 1 -#ifdef __GNUG__ -//#pragma interface -#endif - -#include "RNG.h" -#include - -// -// Multiplicative Linear Conguential Generator -// - -class MLCG : public RNG { - int32_t initialSeedOne; - int32_t initialSeedTwo; - int32_t seedOne; - int32_t seedTwo; - -protected: - -public: - MLCG(int32_t seed1 = 0, int32_t seed2 = 1); - // - // Return a long-words word of random bits - // - virtual uint32_t asLong(); - virtual void reset(); - int32_t seed1(); - void seed1(int32_t); - int32_t seed2(); - void seed2(int32_t); - void reseed(int32_t, int32_t); -}; - -inline int32_t -MLCG::seed1() -{ - return(seedOne); -} - -inline void -MLCG::seed1(int32_t s) -{ - initialSeedOne = s; - reset(); -} - -inline int32_t -MLCG::seed2() -{ - return(seedTwo); -} - -inline void -MLCG::seed2(int32_t s) -{ - initialSeedTwo = s; - reset(); -} - -inline void -MLCG::reseed(int32_t s1, int32_t s2) -{ - initialSeedOne = s1; - initialSeedTwo = s2; - reset(); -} - -#endif diff --git a/src/gnu/isaac64.cpp b/src/gnu/isaac64.cpp deleted file mode 100644 index b0c6c40317..0000000000 --- a/src/gnu/isaac64.cpp +++ /dev/null @@ -1,145 +0,0 @@ -/* ------------------------------------------------------------------------------- -isaac64.cpp: A Fast cryptographic random number generator - for 32-bit and 64-bit machines. -By Bob Jenkins, 1996. Public Domain. - -Modified for modularity by Tom Bartol and Rex Kerr ------------------------------------------------------------------------------- -*/ -#include "isaac64.h" - - -#define ind(mm, x) (*(ub8*) ((ub1*) (mm) + ((x) & ((RANDSIZ - 1) << 3)))) - -#define rngstep(mix, a, b, mm, m, m2, r, x) \ - { \ - x = *m; \ - a = (mix) + *(m2++); \ - *(m++) = y = ind(mm, x) + a + b; \ - *(r++) = b = ind(mm, y >> RANDSIZL) + x; \ - } - -#define mix(a, b, c, d, e, f, g, h) \ - { \ - a -= e; \ - f ^= h >> 9; \ - h += a; \ - b -= f; \ - g ^= a << 9; \ - a += b; \ - c -= g; \ - h ^= b >> 23; \ - b += c; \ - d -= h; \ - a ^= c << 15; \ - c += d; \ - e -= a; \ - b ^= d >> 14; \ - d += e; \ - f -= b; \ - c ^= e << 20; \ - e += f; \ - g -= c; \ - d ^= f >> 17; \ - f += g; \ - h -= d; \ - e ^= g << 14; \ - g += h; \ - } - - -void isaac64_generate(struct isaac64_state* rng) { - ub8 a, b, x, y, *m, *m2, *r, *mend; - - m = rng->mm; - r = rng->randrsl; - a = rng->aa; - b = rng->bb + (++rng->cc); - for (m = rng->mm, mend = m2 = m + (RANDSIZ / 2); m < mend;) { - rngstep(~(a ^ (a << 21)), a, b, rng->mm, m, m2, r, x); - rngstep(a ^ (a >> 5), a, b, rng->mm, m, m2, r, x); - rngstep(a ^ (a << 12), a, b, rng->mm, m, m2, r, x); - rngstep(a ^ (a >> 33), a, b, rng->mm, m, m2, r, x); - } - for (m2 = rng->mm; m2 < mend;) { - rngstep(~(a ^ (a << 21)), a, b, rng->mm, m, m2, r, x); - rngstep(a ^ (a >> 5), a, b, rng->mm, m, m2, r, x); - rngstep(a ^ (a << 12), a, b, rng->mm, m, m2, r, x); - rngstep(a ^ (a >> 33), a, b, rng->mm, m, m2, r, x); - } - rng->bb = b; - rng->aa = a; -} - - -void isaac64_init(struct isaac64_state* rng, ub4 seed) { - ub8 *r, *m; - ub8 a, b, c, d, e, f, g, h; - ub4 i; - - rng->aa = (ub8) 0; - rng->bb = (ub8) 0; - rng->cc = (ub8) 0; - - a = b = c = d = e = f = g = h = 0x9e3779b97f4a7c13LL; /* the golden ratio */ - - r = rng->randrsl; - m = rng->mm; - - for (i = 0; i < RANDSIZ; ++i) - r[i] = (ub8) 0; - - r[0] = seed; - - for (i = 0; i < 4; ++i) /* scramble it */ - { - mix(a, b, c, d, e, f, g, h); - } - - for (i = 0; i < RANDSIZ; i += 8) /* fill in m[] with messy stuff */ - { - /* use all the information in the seed */ - a += r[i]; - b += r[i + 1]; - c += r[i + 2]; - d += r[i + 3]; - e += r[i + 4]; - f += r[i + 5]; - g += r[i + 6]; - h += r[i + 7]; - mix(a, b, c, d, e, f, g, h); - m[i] = a; - m[i + 1] = b; - m[i + 2] = c; - m[i + 3] = d; - m[i + 4] = e; - m[i + 5] = f; - m[i + 6] = g; - m[i + 7] = h; - } - - /* do a second pass to make all of the seed affect all of m[] */ - for (i = 0; i < RANDSIZ; i += 8) { - a += m[i]; - b += m[i + 1]; - c += m[i + 2]; - d += m[i + 3]; - e += m[i + 4]; - f += m[i + 5]; - g += m[i + 6]; - h += m[i + 7]; - mix(a, b, c, d, e, f, g, h); - m[i] = a; - m[i + 1] = b; - m[i + 2] = c; - m[i + 3] = d; - m[i + 4] = e; - m[i + 5] = f; - m[i + 6] = g; - m[i + 7] = h; - } - - isaac64_generate(rng); /* fill in the first set of results */ - rng->randcnt = RANDMAX; /* prepare to use the first set of results */ -} diff --git a/src/gnu/isaac64.h b/src/gnu/isaac64.h deleted file mode 100644 index 4af3b7586e..0000000000 --- a/src/gnu/isaac64.h +++ /dev/null @@ -1,88 +0,0 @@ -#pragma once - -/* ------------------------------------------------------------------------------- -isaac64.h: Definitions for a fast cryptographic random number generator -Bob Jenkins, 1996, Public Domain - -Modified for modularity by Tom Bartol and Rex Kerr - -Reference: -http://burtleburtle.net/bob/rand/isaacafa.html -Jenkins, R.J. (1996) ISAAC, in Fast Software Encryption, vol. 1039, - ed. Gollmann, D. Spinger-Verlag, Cambridge - ------------------------------------------------------------------------------- -*/ - -#include - -#define RANDSIZL (4) /* I recommend 8 for crypto, 4 for simulations */ -#define RANDSIZ (1 << RANDSIZL) -#define RANDMAX (2 * RANDSIZ) - -typedef unsigned long long ub8; -#if defined(uint32_t) -typedef uint32_t ub4; -#else -typedef unsigned int ub4; -#endif -typedef unsigned short int ub2; -typedef unsigned char ub1; - -#define DBL32 (2.3283064365386962890625e-10) -#define DBL53 (1.1102230246251565404236316680908203125e-16) -#define DBL64 (5.42101086242752217003726400434970855712890625e-20) -#define MSK53 0x001FFFFFFFFFFFFFLL - -struct isaac64_state { - int randcnt; - ub8 aa; - ub8 bb; - ub8 cc; - ub8 randrsl[RANDSIZ]; - ub8 mm[RANDSIZ]; -}; - - -void isaac64_init(struct isaac64_state* rng, ub4 seed); - -void isaac64_generate(struct isaac64_state* rng); - -/* ------------------------------------------------------------------------------- -Macros to get individual random numbers ------------------------------------------------------------------------------- -*/ - -#define isaac64_uint32(rng) \ - (rng->randcnt > 0 ? (*(((ub4*) (rng->randrsl)) + (rng->randcnt -= 1))) \ - : (isaac64_generate(rng), \ - rng->randcnt = RANDMAX - 1, \ - *(((ub4*) (rng->randrsl)) + rng->randcnt))) - -#define isaac64_uint64(rng) \ - (rng->randcnt > 1 ? (*((ub8*) (((ub4*) (rng->randrsl)) + (rng->randcnt -= 2)))) \ - : (isaac64_generate(rng), \ - rng->randcnt = RANDMAX - 2, \ - *((ub8*) (((ub4*) (rng->randrsl)) + rng->randcnt)))) - -#define isaac64_dbl32(rng) \ - (rng->randcnt > 0 ? (DBL32 * (*(((ub4*) (rng->randrsl)) + (rng->randcnt -= 1)))) \ - : (isaac64_generate(rng), \ - rng->randcnt = RANDMAX - 1, \ - DBL32 * (*(((ub4*) (rng->randrsl)) + rng->randcnt)))) - -#define isaac64_dbl53(rng) \ - (rng->randcnt > 1 \ - ? (DBL53 * ((*((ub8*) (((ub4*) (rng->randrsl)) + (rng->randcnt -= 2)))) >> 11)) \ - : (isaac64_generate(rng), \ - rng->randcnt = RANDMAX - 2, \ - DBL64 * ((*((ub8*) (((ub4*) (rng->randrsl)) + rng->randcnt))) >> 11))) - -#define isaac64_dbl64(rng) \ - (rng->randcnt > 1 ? (DBL64 * (*((ub8*) (((ub4*) (rng->randrsl)) + (rng->randcnt -= 2))))) \ - : (isaac64_generate(rng), \ - rng->randcnt = RANDMAX - 2, \ - DBL64 * (*((ub8*) (((ub4*) (rng->randrsl)) + rng->randcnt))))) - diff --git a/src/gnu/nrnisaac.cpp b/src/gnu/nrnisaac.cpp deleted file mode 100644 index f0549a4d9a..0000000000 --- a/src/gnu/nrnisaac.cpp +++ /dev/null @@ -1,25 +0,0 @@ -#include -#include "nrnisaac.h" -#include "isaac64.h" - -using RNG = struct isaac64_state; - -void* nrnisaac_new() { - return new RNG; -} - -void nrnisaac_delete(void* v) { - delete static_cast(v); -} - -void nrnisaac_init(void* v, unsigned long int seed) { - isaac64_init(static_cast(v), seed); -} - -double nrnisaac_dbl_pick(void* v) { - return isaac64_dbl32(static_cast(v)); -} - -std::uint32_t nrnisaac_uint32_pick(void* v) { - return isaac64_uint32(static_cast(v)); -} diff --git a/src/gnu/nrnisaac.h b/src/gnu/nrnisaac.h deleted file mode 100644 index 5322febafe..0000000000 --- a/src/gnu/nrnisaac.h +++ /dev/null @@ -1,9 +0,0 @@ -#pragma once - -#include - -void* nrnisaac_new(); -void nrnisaac_delete(void* rng); -void nrnisaac_init(void* rng, unsigned long int seed); -double nrnisaac_dbl_pick(void* rng); -std::uint32_t nrnisaac_uint32_pick(void* rng); diff --git a/src/ivoc/ivocmain.cpp b/src/ivoc/ivocmain.cpp index 315a375eda..ac0f4fc599 100644 --- a/src/ivoc/ivocmain.cpp +++ b/src/ivoc/ivocmain.cpp @@ -233,14 +233,6 @@ extern void nrnmpi_stubs(); extern std::string nrnmpi_load(); #endif -// some things are defined in libraries earlier than they are used so... -#include -static void force_load() { - if (always_false) { - nrnisaac_new(); - } -} - #ifdef MINGW // see iv/src/OS/directory.cpp #include @@ -381,7 +373,6 @@ int ivocmain_session(int argc, const char** argv, const char** env, int start_se // extern char** environ; int i; // prargs("at beginning", argc, argv); - force_load(); nrn_global_argc = argc; // https://en.cppreference.com/w/cpp/language/main_function, note that argv is // of length argc + 1 and argv[argc] is null. diff --git a/src/ivoc/ivocrand.cpp b/src/ivoc/ivocrand.cpp index 9661c00650..c44455864b 100644 --- a/src/ivoc/ivocrand.cpp +++ b/src/ivoc/ivocrand.cpp @@ -9,7 +9,6 @@ #include #include "classreg.h" #include "oc2iv.h" -#include "nrnisaac.h" #include "utils/enumerate.h" #include diff --git a/src/nrniv/multisend_setup.cpp b/src/nrniv/multisend_setup.cpp index 479d30bf78..b57c07ddc7 100644 --- a/src/nrniv/multisend_setup.cpp +++ b/src/nrniv/multisend_setup.cpp @@ -94,18 +94,18 @@ void TarList::alloc() { } } -#include -static void* ranstate; +#include +static nrnran123_State* ranstate; static void random_init(int i) { if (!ranstate) { - ranstate = nrnisaac_new(); + ranstate = nrnran123_newstream(0, 0); } - nrnisaac_init(ranstate, nrnmpi_myid + 1); + nrnran123_setseq(ranstate, nrnmpi_myid + 1, 0); } static unsigned int get_random() { - return nrnisaac_uint32_pick(ranstate); + return nrnran123_ipick(ranstate); } static int iran(int i1, int i2) { From 2fc44a49150767db5279db5b01ef33c67c54e5e8 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Wed, 31 Jan 2024 18:54:07 +0100 Subject: [PATCH 03/12] Remove usage of MCellRan4 --- cmake/NeuronFileLists.cmake | 1 - src/gnu/CMakeLists.txt | 1 - src/gnu/mcran4.cpp | 139 ------------------------------------ src/gnu/mcran4.h | 40 ----------- src/nrniv/kschan.cpp | 9 ++- src/nrniv/kssingle.cpp | 2 +- src/nrniv/kssingle.h | 9 ++- src/oc/hoc_init.cpp | 1 - src/oc/oc_mcran4.cpp | 30 +++----- src/oc/oc_mcran4.hpp | 2 - src/oc/scoprand.cpp | 16 +---- src/scopmath/praxis.cpp | 64 +++-------------- 12 files changed, 33 insertions(+), 281 deletions(-) delete mode 100644 src/gnu/mcran4.cpp delete mode 100644 src/gnu/mcran4.h diff --git a/cmake/NeuronFileLists.cmake b/cmake/NeuronFileLists.cmake index 84cd4786ac..1a18db8807 100644 --- a/cmake/NeuronFileLists.cmake +++ b/cmake/NeuronFileLists.cmake @@ -9,7 +9,6 @@ set(STRUCTURED_HEADER_FILES_TO_INSTALL neuron/container/generic_data_handle.hpp neuron/container/non_owning_soa_identifier.hpp neuron/model_data_fwd.hpp) set(HEADER_FILES_TO_INSTALL - gnu/mcran4.h gnu/nrnran123.h nrniv/backtrace_utils.h nrniv/bbsavestate.h diff --git a/src/gnu/CMakeLists.txt b/src/gnu/CMakeLists.txt index 1d8783a621..e2f31b85cf 100644 --- a/src/gnu/CMakeLists.txt +++ b/src/gnu/CMakeLists.txt @@ -6,7 +6,6 @@ add_library( Geom.cpp HypGeom.cpp LogNorm.cpp - mcran4.cpp NegExp.cpp Normal.cpp NrnRandom123RNG.cpp diff --git a/src/gnu/mcran4.cpp b/src/gnu/mcran4.cpp deleted file mode 100644 index a9553687f6..0000000000 --- a/src/gnu/mcran4.cpp +++ /dev/null @@ -1,139 +0,0 @@ -/* -The copyrighted code from Numerical Recipes Software has been removed -and replaced by an independent implementation found in the random.cpp file -in function Ranint4 -from http://www.inference.phy.cam.ac.uk/bayesys/BayesSys.tar.gz -by John Skilling -http://www.inference.phy.cam.ac.uk/bayesys/ -The code fragment was further modified by Michael Hines to change prototypes -and produce output identical to the old version. This code is now -placed under the General GNU Public License, version 2. The random.cpp file -contained the header: -//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -// Filename: random.cpp -// -// Purpose: Random utility procedures for BayeSys3. -// -// History: Random.cpp 17 Nov 1994 - 13 Sep 2003 -// -// Acknowledgments: -// "Numerical Recipes", Press et al, for ideas -// "Handbook of Mathematical Functions", Abramowitz and Stegun, for formulas -// Peter J Acklam website, for inverse normal code -//----------------------------------------------------------------------------- - Copyright (c) 1994-2003 Maximum Entropy Data Consultants Ltd, - 114c Milton Road, Cambridge CB4 1XE, England - - This library is free software; you can redistribute it and/or - modify it under the terms of the GNU Lesser General Public - License as published by the Free Software Foundation; either - version 2.1 of the License, or (at your option) any later version. - - This library 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 - Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public - License along with this library; if not, write to the Free Software - Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -*/ - -#include -#include -#include -#include -#include -#include "mcran4.h" - -static uint32_t lowindex = 0; - -double mcell_lowindex() { - return static_cast(lowindex); -} - -void mcell_ran4_init(uint32_t low) { - lowindex = low; -} - -double mcell_ran4(uint32_t* high, double* x, unsigned int n, double range) { - int i; - for (i = 0; i < n; i++) { - x[i] = range * nrnRan4dbl(high, lowindex); - } - return x[0]; -} - -double mcell_ran4a(uint32_t* high) { - return nrnRan4dbl(high, lowindex); -} - -uint32_t mcell_iran4(uint32_t* high) { - return nrnRan4int(high, lowindex); -} - -uint32_t nrnRan4int(uint32_t* idx1, uint32_t idx2) { - uint32_t u, v, w, m, n; - /* 64-bit hash */ - n = (*idx1)++; - m = idx2; - - w = n ^ 0xbaa96887; - v = w >> 16; - w &= 0xffff; - u = ~((v - w) * (v + w)); - /*m ^= (((u >> 16) | (u << 16)) ^ 0xb4f0c4a7) + w * v;*/ - m ^= (((u >> 16) | (u << 16)) ^ 0x4b0f3b58) + w * v; - - w = m ^ 0x1e17d32c; - v = w >> 16; - w &= 0xffff; - u = ~((v - w) * (v + w)); - /*n ^= (((u >> 16) | (u << 16)) ^ 0x178b0f3c) + w * v;*/ - n ^= (((u >> 16) | (u << 16)) ^ 0xe874f0c3) + w * v; - return n; - - w = n ^ 0x03bcdc3c; - v = w >> 16; - w &= 0xffff; - u = (v - w) * (v + w); - m ^= (((u >> 16) | (u << 16)) ^ 0x96aa3a59) + w * v; - - w = m ^ 0x0f33d1b2; - v = w >> 16; - w &= 0xffff; - u = (v - w) * (v + w); - n ^= (((u >> 16) | (u << 16)) ^ 0xaa5835b9) + w * v; - return n; -} - -/* -//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -// Function: Randouble (hines now nrnRan4dbl -// -// Purpose: Random double-precision floating-point sample. -// The 2^52 allowed values are odd multiples of 2^-53, -// symmetrically placed in strict interior of (0,1). -// -// Notes: (1) Tuned to 52-bit mantissa in "double" format. -// (2) Uses one call to Ranint to get 64 random bits, with extra -// random integer available in Rand[3]. -// (3) All floating-point random calls are directed through this code, -// except Rangauss which uses the extra random integer in Rand[3]. -// -// History: John Skilling 6 May 1995, 3 Dec 1995, 24 Aug 1996 -// 20 Oct 2002, 17 Dec 2002 -//----------------------------------------------------------------------------- -// -*/ -static const double SHIFT32 = 1.0 / 4294967296.0; /* 2^-32 */ -double nrnRan4dbl(uint32_t* idx1, uint32_t idx2) { - uint32_t hi, lo, extra; - hi = (uint32_t) nrnRan4int(idx1, idx2); /*top 32 bits*/ - /* - // lo = (extra // low bits - // & 0xfffff000) ^ 0x00000800; // switch lowest (2^-53) bit ON - // return ((double)hi + (double)lo * SHIFT32) * SHIFT32; - */ - return ((double) hi) * SHIFT32; -} diff --git a/src/gnu/mcran4.h b/src/gnu/mcran4.h deleted file mode 100644 index 9d7ff79dd9..0000000000 --- a/src/gnu/mcran4.h +++ /dev/null @@ -1,40 +0,0 @@ -#pragma once -#include - -double mcell_lowindex(); -void mcell_ran4_init(uint32_t); -double mcell_ran4(uint32_t* idx1, double* x, unsigned int n, double range); -double mcell_ran4a(uint32_t* idx1); -uint32_t mcell_iran4(uint32_t* idx1); -double nrnRan4dbl(uint32_t* idx1, uint32_t idx2); -uint32_t nrnRan4int(uint32_t* idx1, uint32_t idx2); - -/* - The original ran4 generator was copyrighted by "Numerical Recipes in C" - and therefore has been removed from the NEURON sources and replaced by code - fragments obtained from http://www.inference.phy.cam.ac.uk/bayesys/ - by John Skilling - The function mcell_ran4a only returns - a single uniform random number in the distribution 0.0 to 1.0 . - The prototype for mcell_ran4 is the original. -*/ - -/* -Michael Hines added the prefix mcell to the global names. -These functions were obtained from Tom Bartol -who uses them in his mcell program. -He comments: -For MCell, Ran4 has the distinct advantage of generating -streams of random bits not just random numbers. This means that you can -break the 32 bits of a single returned random number into several smaller -chunks without fear of correlations between the chunks. Ran4 is not the -fastest generator in the universe but it's pretty fast (16 million -floating point random numbers per second on my 1GHz Intel PIII and 20 -million integer random numbers per second) and of near cryptographic -quality. I've modified it so that a given seed will generate the same -sequence on Intel, Alpha, RS6000, PowerPC, and MIPS architectures (and -probably anything else out there). It's also been modified to generate -arbitrary length vectors of random numbers. This makes generating numbers -more efficient because you can generate many numbers per function call. -MCell generates them in chunks of 10000 at a time. -*/ diff --git a/src/nrniv/kschan.cpp b/src/nrniv/kschan.cpp index b3ab59ddc1..05a6bbf5ff 100644 --- a/src/nrniv/kschan.cpp +++ b/src/nrniv/kschan.cpp @@ -279,9 +279,12 @@ static double ks_vres(void* v) { static double ks_rseed(void* v) { if (ifarg(1)) { - KSSingle::idum_ = (unsigned int) chkarg(1, 0, 1e9); + nrnran123_setseq(KSSingle::rand_state, static_cast(chkarg(1, 0, 1e9)), 0); } - return (double) KSSingle::idum_; + std::uint32_t seq; + char which; + nrnran123_getseq(KSSingle::rand_state, &seq, &which); + return (double) seq; } static double ks_usetable(void* v) { @@ -804,7 +807,7 @@ void KSChan_reg() { ksgate_sym = hoc_lookup("KSGate"); kstrans_sym = hoc_lookup("KSTrans"); KSSingle::vres_ = 0.1; - KSSingle::idum_ = 0; + KSSingle::rand_state = nrnran123_newstream(0, 0); } // param is gmax, g, i --- if change then change numbers below diff --git a/src/nrniv/kssingle.cpp b/src/nrniv/kssingle.cpp index f23267ee33..32db40e581 100644 --- a/src/nrniv/kssingle.cpp +++ b/src/nrniv/kssingle.cpp @@ -130,7 +130,7 @@ extern NetCvode* net_cvode_instance; double KSSingle::vres_; -unsigned int KSSingle::idum_; +nrnran123_State* KSSingle::rand_state; KSSingle::KSSingle(KSChan* c) { // implemenation assumes one ks gate complex with power 1 diff --git a/src/nrniv/kssingle.h b/src/nrniv/kssingle.h index 22034ee696..dbd2139aae 100644 --- a/src/nrniv/kssingle.h +++ b/src/nrniv/kssingle.h @@ -5,8 +5,7 @@ #include #include #include - -#include +#include class KSSingleTrans; class KSSingleState; @@ -71,10 +70,10 @@ class KSSingle { return MyMath::eq(x, y, vres_); } double exprand() { - return -log(mcell_ran4a(&idum_)); + return nrnran123_negexp(rand_state); } double unifrand(double range) { - return mcell_ran4a(&idum_) * range; + return nrnran123_dblpick(rand_state) * range; } int rvalrand(int); @@ -84,7 +83,7 @@ class KSSingle { double* rval_; bool uses_ligands_; static double vres_; - static unsigned int idum_; + static nrnran123_State* rand_state; static unsigned long singleevent_deliver_; static unsigned long singleevent_move_; diff --git a/src/oc/hoc_init.cpp b/src/oc/hoc_init.cpp index 17675d7753..455630a653 100644 --- a/src/oc/hoc_init.cpp +++ b/src/oc/hoc_init.cpp @@ -255,7 +255,6 @@ void hoc_init(void) /* install constants and built-ins table */ } } - set_use_mcran4(false); nrn_xopen_broadcast_ = 255; extern void hoc_init_space(void); hoc_init_space(); diff --git a/src/oc/oc_mcran4.cpp b/src/oc/oc_mcran4.cpp index 128ccd7d04..2bd1ef043d 100644 --- a/src/oc/oc_mcran4.cpp +++ b/src/oc/oc_mcran4.cpp @@ -1,15 +1,7 @@ #include "hocdec.h" -#include "mcran4.h" +#include "nrnran123.h" -int use_mcell_ran4_; - -void set_use_mcran4(bool value) { - use_mcell_ran4_ = value ? 1 : 0; -} - -bool use_mcran4() { - return use_mcell_ran4_ != 0; -} +static const nrnran123_State* state = nrnran123_newstream(0, 0); void hoc_mcran4() { uint32_t idx; @@ -17,25 +9,25 @@ void hoc_mcran4() { double x; xidx = hoc_pgetarg(1); idx = (uint32_t) (*xidx); - x = mcell_ran4a(&idx); - *xidx = idx; + nrnran123_setseq(state, idx, 0); + x = nrnran123_dblpick(state); + *xidx = idx + 1; hoc_ret(); hoc_pushx(x); } void hoc_mcran4init() { - double prev = mcell_lowindex(); + std::uint32_t seq; + char which; + nrnran123_getseq(state, &seq, &which); + double prev = static_cast(seq); if (ifarg(1)) { uint32_t idx = (uint32_t) chkarg(1, 0., 4294967295.); - mcell_ran4_init(idx); + nrnran123_setseq(state, idx, 0); } hoc_ret(); hoc_pushx(prev); } void hoc_usemcran4() { - double prev = (double) use_mcell_ran4_; - if (ifarg(1)) { - use_mcell_ran4_ = (int) chkarg(1, 0., 1.); - } hoc_ret(); - hoc_pushx(prev); + hoc_pushx(0); } diff --git a/src/oc/oc_mcran4.hpp b/src/oc/oc_mcran4.hpp index d275e3e9bf..1dca3b6feb 100644 --- a/src/oc/oc_mcran4.hpp +++ b/src/oc/oc_mcran4.hpp @@ -1,5 +1,3 @@ -void set_use_mcran4(bool value); -bool use_mcran4(); void hoc_mcran4(); void hoc_mcran4init(); void hoc_usemcran4(); diff --git a/src/oc/scoprand.cpp b/src/oc/scoprand.cpp index 70651f5608..37e55ee789 100644 --- a/src/oc/scoprand.cpp +++ b/src/oc/scoprand.cpp @@ -23,8 +23,7 @@ static char RCSid[] = "random.cpp,v 1.4 1999/01/04 12:46:49 hines Exp"; #endif #include -#include "oc_mcran4.hpp" -#include "mcran4.h" +#include "nrnran123.h" #include "scoplib.h" static uint32_t value = 1; @@ -59,17 +58,8 @@ static uint32_t value = 1; *--------------------------------------------------------------------------- */ double scop_random(void) { - if (use_mcran4()) { - /*perhaps 4 times slower but much higher quality*/ - return mcell_ran4a(&value); - } else { - uint32_t a = 2147437301, c = 453816981, - /* m = 2^32 - 1, the largest long int value that can be represented */ - /*m = 0xFFFFFFFF;*/ /* limited to 32 bit integers*/ - m = ~0; - value = a * value + c; - return (fabs((double) value / (double) m)); - } + auto *rng = nrnran123_newstream(value, 0); + return nrnran123_dblpick(rng); } /*----------------------------------------------------------------------------- diff --git a/src/scopmath/praxis.cpp b/src/scopmath/praxis.cpp index effb3f8844..daea1160bd 100644 --- a/src/scopmath/praxis.cpp +++ b/src/scopmath/praxis.cpp @@ -1,13 +1,13 @@ #include <../../nrnconf.h> -#include "mcran4.h" -#include "oc_ansi.h" -#include "scoplib.h" - #include #include #include #include +#include "oc_ansi.h" +#include "scoplib.h" +#include + #undef small typedef long int integer; @@ -1730,55 +1730,7 @@ static /* Subroutine */ int maprnt_(integer* option, doublereal* v, integer* m, uint32_t nrn_praxis_ran_index; -static doublereal random_(integer* naught) { - double x; - return mcell_ran4(&nrn_praxis_ran_index, &x, 1, 1.); -#if 0 - /* Initialized data */ - - static logical init = FALSE_; - - /* System generated locals */ - doublereal ret_val; - - /* Local variables */ - static doublereal half; - static integer i, j, q, r; - static doublereal ran1; - static integer ran2; - static doublereal ran3[127]; - - if (init) { - goto L3; - } - r = *naught % 8190 + 1; - ran2 = 128; - for (i = 1; i <= 127; ++i) { - --ran2; - ran1 = -36028797018963968.; - for (j = 1; j <= 7; ++j) { - r = r * 1756 % 8191; - q = r / 32; -/* L1: */ - ran1 = (ran1 + q) * .00390625; - } -/* L2: */ - ran3[ran2 - 1] = ran1; - } - init = TRUE_; -L3: - if (ran2 == 1) { - ran2 = 128; - } - --ran2; - ran1 += ran3[ran2 - 1]; - half = .5; - if (ran1 >= 0.) { - half = -half; - } - ran1 += half; - ran3[ran2 - 1] = ran1; - ret_val = ran1 + .5; - return ret_val; -#endif -} /* random_ */ +static doublereal random_(integer* /* naught */) { + auto *rng = nrnran123_newstream(nrn_praxis_ran_index, 0); + return nrnran123_dblpick(rng); +} From 6a6ab8581b091818be067b32d3ca07ae00436a52 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Wed, 31 Jan 2024 19:47:44 +0100 Subject: [PATCH 04/12] Fix mcran --- src/oc/oc_mcran4.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/oc/oc_mcran4.cpp b/src/oc/oc_mcran4.cpp index 2bd1ef043d..18f8bc8113 100644 --- a/src/oc/oc_mcran4.cpp +++ b/src/oc/oc_mcran4.cpp @@ -1,7 +1,7 @@ #include "hocdec.h" #include "nrnran123.h" -static const nrnran123_State* state = nrnran123_newstream(0, 0); +static nrnran123_State* state = nrnran123_newstream(0, 0); void hoc_mcran4() { uint32_t idx; From 8ee380ae3932a0abd236e68c9fb3f63623cdd9fa Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Wed, 31 Jan 2024 19:48:19 +0100 Subject: [PATCH 05/12] Format --- src/nrniv/kschan.cpp | 2 +- src/oc/scoprand.cpp | 2 +- src/scopmath/praxis.cpp | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/nrniv/kschan.cpp b/src/nrniv/kschan.cpp index 05a6bbf5ff..8f9356bed6 100644 --- a/src/nrniv/kschan.cpp +++ b/src/nrniv/kschan.cpp @@ -279,7 +279,7 @@ static double ks_vres(void* v) { static double ks_rseed(void* v) { if (ifarg(1)) { - nrnran123_setseq(KSSingle::rand_state, static_cast(chkarg(1, 0, 1e9)), 0); + nrnran123_setseq(KSSingle::rand_state, static_cast(chkarg(1, 0, 1e9)), 0); } std::uint32_t seq; char which; diff --git a/src/oc/scoprand.cpp b/src/oc/scoprand.cpp index 37e55ee789..f6fd48dc3c 100644 --- a/src/oc/scoprand.cpp +++ b/src/oc/scoprand.cpp @@ -58,7 +58,7 @@ static uint32_t value = 1; *--------------------------------------------------------------------------- */ double scop_random(void) { - auto *rng = nrnran123_newstream(value, 0); + auto* rng = nrnran123_newstream(value, 0); return nrnran123_dblpick(rng); } diff --git a/src/scopmath/praxis.cpp b/src/scopmath/praxis.cpp index daea1160bd..67ae03eb96 100644 --- a/src/scopmath/praxis.cpp +++ b/src/scopmath/praxis.cpp @@ -1731,6 +1731,6 @@ static /* Subroutine */ int maprnt_(integer* option, doublereal* v, integer* m, uint32_t nrn_praxis_ran_index; static doublereal random_(integer* /* naught */) { - auto *rng = nrnran123_newstream(nrn_praxis_ran_index, 0); + auto* rng = nrnran123_newstream(nrn_praxis_ran_index, 0); return nrnran123_dblpick(rng); } From c03f943950a73068d256e06510341d485e9844f7 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Wed, 31 Jan 2024 20:24:57 +0100 Subject: [PATCH 06/12] Take care of mech_api too --- src/oc/mech_api.h | 1 - 1 file changed, 1 deletion(-) diff --git a/src/oc/mech_api.h b/src/oc/mech_api.h index 8d3ed11ecd..f7ff66ea0c 100644 --- a/src/oc/mech_api.h +++ b/src/oc/mech_api.h @@ -10,7 +10,6 @@ * we leave to nocmodl. */ #include "bbsavestate.h" -#include "mcran4.h" #include "nrncvode.h" #include "nrnran123.h" #include "nrnrandom.h" From 301edcfafd5b3c45bf9f9eed339fc50b64d79609 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 1 Feb 2024 10:00:11 +0100 Subject: [PATCH 07/12] Revert "Take care of mech_api too" This reverts commit c03f943950a73068d256e06510341d485e9844f7. --- src/oc/mech_api.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/oc/mech_api.h b/src/oc/mech_api.h index f7ff66ea0c..8d3ed11ecd 100644 --- a/src/oc/mech_api.h +++ b/src/oc/mech_api.h @@ -10,6 +10,7 @@ * we leave to nocmodl. */ #include "bbsavestate.h" +#include "mcran4.h" #include "nrncvode.h" #include "nrnran123.h" #include "nrnrandom.h" From f676951e2945428286fc1152e9c58c05b8b6fb65 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 1 Feb 2024 10:00:27 +0100 Subject: [PATCH 08/12] Revert "Format" This reverts commit 8ee380ae3932a0abd236e68c9fb3f63623cdd9fa. --- src/nrniv/kschan.cpp | 2 +- src/oc/scoprand.cpp | 2 +- src/scopmath/praxis.cpp | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/nrniv/kschan.cpp b/src/nrniv/kschan.cpp index 8f9356bed6..05a6bbf5ff 100644 --- a/src/nrniv/kschan.cpp +++ b/src/nrniv/kschan.cpp @@ -279,7 +279,7 @@ static double ks_vres(void* v) { static double ks_rseed(void* v) { if (ifarg(1)) { - nrnran123_setseq(KSSingle::rand_state, static_cast(chkarg(1, 0, 1e9)), 0); + nrnran123_setseq(KSSingle::rand_state, static_cast(chkarg(1, 0, 1e9)), 0); } std::uint32_t seq; char which; diff --git a/src/oc/scoprand.cpp b/src/oc/scoprand.cpp index f6fd48dc3c..37e55ee789 100644 --- a/src/oc/scoprand.cpp +++ b/src/oc/scoprand.cpp @@ -58,7 +58,7 @@ static uint32_t value = 1; *--------------------------------------------------------------------------- */ double scop_random(void) { - auto* rng = nrnran123_newstream(value, 0); + auto *rng = nrnran123_newstream(value, 0); return nrnran123_dblpick(rng); } diff --git a/src/scopmath/praxis.cpp b/src/scopmath/praxis.cpp index 67ae03eb96..daea1160bd 100644 --- a/src/scopmath/praxis.cpp +++ b/src/scopmath/praxis.cpp @@ -1731,6 +1731,6 @@ static /* Subroutine */ int maprnt_(integer* option, doublereal* v, integer* m, uint32_t nrn_praxis_ran_index; static doublereal random_(integer* /* naught */) { - auto* rng = nrnran123_newstream(nrn_praxis_ran_index, 0); + auto *rng = nrnran123_newstream(nrn_praxis_ran_index, 0); return nrnran123_dblpick(rng); } From eca7536321969cbca81d174192254426d7ff169c Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 1 Feb 2024 10:00:41 +0100 Subject: [PATCH 09/12] Revert "Fix mcran" This reverts commit 6a6ab8581b091818be067b32d3ca07ae00436a52. --- src/oc/oc_mcran4.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/oc/oc_mcran4.cpp b/src/oc/oc_mcran4.cpp index 18f8bc8113..2bd1ef043d 100644 --- a/src/oc/oc_mcran4.cpp +++ b/src/oc/oc_mcran4.cpp @@ -1,7 +1,7 @@ #include "hocdec.h" #include "nrnran123.h" -static nrnran123_State* state = nrnran123_newstream(0, 0); +static const nrnran123_State* state = nrnran123_newstream(0, 0); void hoc_mcran4() { uint32_t idx; From fee9977ed1fbdbb5a0a8e2de766549d568cabbe2 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 1 Feb 2024 10:00:52 +0100 Subject: [PATCH 10/12] Revert "Remove usage of MCellRan4" This reverts commit 2fc44a49150767db5279db5b01ef33c67c54e5e8. --- cmake/NeuronFileLists.cmake | 1 + src/gnu/CMakeLists.txt | 1 + src/gnu/mcran4.cpp | 139 ++++++++++++++++++++++++++++++++++++ src/gnu/mcran4.h | 40 +++++++++++ src/nrniv/kschan.cpp | 9 +-- src/nrniv/kssingle.cpp | 2 +- src/nrniv/kssingle.h | 9 +-- src/oc/hoc_init.cpp | 1 + src/oc/oc_mcran4.cpp | 30 +++++--- src/oc/oc_mcran4.hpp | 2 + src/oc/scoprand.cpp | 16 ++++- src/scopmath/praxis.cpp | 64 ++++++++++++++--- 12 files changed, 281 insertions(+), 33 deletions(-) create mode 100644 src/gnu/mcran4.cpp create mode 100644 src/gnu/mcran4.h diff --git a/cmake/NeuronFileLists.cmake b/cmake/NeuronFileLists.cmake index 1a18db8807..84cd4786ac 100644 --- a/cmake/NeuronFileLists.cmake +++ b/cmake/NeuronFileLists.cmake @@ -9,6 +9,7 @@ set(STRUCTURED_HEADER_FILES_TO_INSTALL neuron/container/generic_data_handle.hpp neuron/container/non_owning_soa_identifier.hpp neuron/model_data_fwd.hpp) set(HEADER_FILES_TO_INSTALL + gnu/mcran4.h gnu/nrnran123.h nrniv/backtrace_utils.h nrniv/bbsavestate.h diff --git a/src/gnu/CMakeLists.txt b/src/gnu/CMakeLists.txt index e2f31b85cf..1d8783a621 100644 --- a/src/gnu/CMakeLists.txt +++ b/src/gnu/CMakeLists.txt @@ -6,6 +6,7 @@ add_library( Geom.cpp HypGeom.cpp LogNorm.cpp + mcran4.cpp NegExp.cpp Normal.cpp NrnRandom123RNG.cpp diff --git a/src/gnu/mcran4.cpp b/src/gnu/mcran4.cpp new file mode 100644 index 0000000000..a9553687f6 --- /dev/null +++ b/src/gnu/mcran4.cpp @@ -0,0 +1,139 @@ +/* +The copyrighted code from Numerical Recipes Software has been removed +and replaced by an independent implementation found in the random.cpp file +in function Ranint4 +from http://www.inference.phy.cam.ac.uk/bayesys/BayesSys.tar.gz +by John Skilling +http://www.inference.phy.cam.ac.uk/bayesys/ +The code fragment was further modified by Michael Hines to change prototypes +and produce output identical to the old version. This code is now +placed under the General GNU Public License, version 2. The random.cpp file +contained the header: +//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +// Filename: random.cpp +// +// Purpose: Random utility procedures for BayeSys3. +// +// History: Random.cpp 17 Nov 1994 - 13 Sep 2003 +// +// Acknowledgments: +// "Numerical Recipes", Press et al, for ideas +// "Handbook of Mathematical Functions", Abramowitz and Stegun, for formulas +// Peter J Acklam website, for inverse normal code +//----------------------------------------------------------------------------- + Copyright (c) 1994-2003 Maximum Entropy Data Consultants Ltd, + 114c Milton Road, Cambridge CB4 1XE, England + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library 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 + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + +#include +#include +#include +#include +#include +#include "mcran4.h" + +static uint32_t lowindex = 0; + +double mcell_lowindex() { + return static_cast(lowindex); +} + +void mcell_ran4_init(uint32_t low) { + lowindex = low; +} + +double mcell_ran4(uint32_t* high, double* x, unsigned int n, double range) { + int i; + for (i = 0; i < n; i++) { + x[i] = range * nrnRan4dbl(high, lowindex); + } + return x[0]; +} + +double mcell_ran4a(uint32_t* high) { + return nrnRan4dbl(high, lowindex); +} + +uint32_t mcell_iran4(uint32_t* high) { + return nrnRan4int(high, lowindex); +} + +uint32_t nrnRan4int(uint32_t* idx1, uint32_t idx2) { + uint32_t u, v, w, m, n; + /* 64-bit hash */ + n = (*idx1)++; + m = idx2; + + w = n ^ 0xbaa96887; + v = w >> 16; + w &= 0xffff; + u = ~((v - w) * (v + w)); + /*m ^= (((u >> 16) | (u << 16)) ^ 0xb4f0c4a7) + w * v;*/ + m ^= (((u >> 16) | (u << 16)) ^ 0x4b0f3b58) + w * v; + + w = m ^ 0x1e17d32c; + v = w >> 16; + w &= 0xffff; + u = ~((v - w) * (v + w)); + /*n ^= (((u >> 16) | (u << 16)) ^ 0x178b0f3c) + w * v;*/ + n ^= (((u >> 16) | (u << 16)) ^ 0xe874f0c3) + w * v; + return n; + + w = n ^ 0x03bcdc3c; + v = w >> 16; + w &= 0xffff; + u = (v - w) * (v + w); + m ^= (((u >> 16) | (u << 16)) ^ 0x96aa3a59) + w * v; + + w = m ^ 0x0f33d1b2; + v = w >> 16; + w &= 0xffff; + u = (v - w) * (v + w); + n ^= (((u >> 16) | (u << 16)) ^ 0xaa5835b9) + w * v; + return n; +} + +/* +//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +// Function: Randouble (hines now nrnRan4dbl +// +// Purpose: Random double-precision floating-point sample. +// The 2^52 allowed values are odd multiples of 2^-53, +// symmetrically placed in strict interior of (0,1). +// +// Notes: (1) Tuned to 52-bit mantissa in "double" format. +// (2) Uses one call to Ranint to get 64 random bits, with extra +// random integer available in Rand[3]. +// (3) All floating-point random calls are directed through this code, +// except Rangauss which uses the extra random integer in Rand[3]. +// +// History: John Skilling 6 May 1995, 3 Dec 1995, 24 Aug 1996 +// 20 Oct 2002, 17 Dec 2002 +//----------------------------------------------------------------------------- +// +*/ +static const double SHIFT32 = 1.0 / 4294967296.0; /* 2^-32 */ +double nrnRan4dbl(uint32_t* idx1, uint32_t idx2) { + uint32_t hi, lo, extra; + hi = (uint32_t) nrnRan4int(idx1, idx2); /*top 32 bits*/ + /* + // lo = (extra // low bits + // & 0xfffff000) ^ 0x00000800; // switch lowest (2^-53) bit ON + // return ((double)hi + (double)lo * SHIFT32) * SHIFT32; + */ + return ((double) hi) * SHIFT32; +} diff --git a/src/gnu/mcran4.h b/src/gnu/mcran4.h new file mode 100644 index 0000000000..9d7ff79dd9 --- /dev/null +++ b/src/gnu/mcran4.h @@ -0,0 +1,40 @@ +#pragma once +#include + +double mcell_lowindex(); +void mcell_ran4_init(uint32_t); +double mcell_ran4(uint32_t* idx1, double* x, unsigned int n, double range); +double mcell_ran4a(uint32_t* idx1); +uint32_t mcell_iran4(uint32_t* idx1); +double nrnRan4dbl(uint32_t* idx1, uint32_t idx2); +uint32_t nrnRan4int(uint32_t* idx1, uint32_t idx2); + +/* + The original ran4 generator was copyrighted by "Numerical Recipes in C" + and therefore has been removed from the NEURON sources and replaced by code + fragments obtained from http://www.inference.phy.cam.ac.uk/bayesys/ + by John Skilling + The function mcell_ran4a only returns + a single uniform random number in the distribution 0.0 to 1.0 . + The prototype for mcell_ran4 is the original. +*/ + +/* +Michael Hines added the prefix mcell to the global names. +These functions were obtained from Tom Bartol +who uses them in his mcell program. +He comments: +For MCell, Ran4 has the distinct advantage of generating +streams of random bits not just random numbers. This means that you can +break the 32 bits of a single returned random number into several smaller +chunks without fear of correlations between the chunks. Ran4 is not the +fastest generator in the universe but it's pretty fast (16 million +floating point random numbers per second on my 1GHz Intel PIII and 20 +million integer random numbers per second) and of near cryptographic +quality. I've modified it so that a given seed will generate the same +sequence on Intel, Alpha, RS6000, PowerPC, and MIPS architectures (and +probably anything else out there). It's also been modified to generate +arbitrary length vectors of random numbers. This makes generating numbers +more efficient because you can generate many numbers per function call. +MCell generates them in chunks of 10000 at a time. +*/ diff --git a/src/nrniv/kschan.cpp b/src/nrniv/kschan.cpp index 05a6bbf5ff..b3ab59ddc1 100644 --- a/src/nrniv/kschan.cpp +++ b/src/nrniv/kschan.cpp @@ -279,12 +279,9 @@ static double ks_vres(void* v) { static double ks_rseed(void* v) { if (ifarg(1)) { - nrnran123_setseq(KSSingle::rand_state, static_cast(chkarg(1, 0, 1e9)), 0); + KSSingle::idum_ = (unsigned int) chkarg(1, 0, 1e9); } - std::uint32_t seq; - char which; - nrnran123_getseq(KSSingle::rand_state, &seq, &which); - return (double) seq; + return (double) KSSingle::idum_; } static double ks_usetable(void* v) { @@ -807,7 +804,7 @@ void KSChan_reg() { ksgate_sym = hoc_lookup("KSGate"); kstrans_sym = hoc_lookup("KSTrans"); KSSingle::vres_ = 0.1; - KSSingle::rand_state = nrnran123_newstream(0, 0); + KSSingle::idum_ = 0; } // param is gmax, g, i --- if change then change numbers below diff --git a/src/nrniv/kssingle.cpp b/src/nrniv/kssingle.cpp index 32db40e581..f23267ee33 100644 --- a/src/nrniv/kssingle.cpp +++ b/src/nrniv/kssingle.cpp @@ -130,7 +130,7 @@ extern NetCvode* net_cvode_instance; double KSSingle::vres_; -nrnran123_State* KSSingle::rand_state; +unsigned int KSSingle::idum_; KSSingle::KSSingle(KSChan* c) { // implemenation assumes one ks gate complex with power 1 diff --git a/src/nrniv/kssingle.h b/src/nrniv/kssingle.h index dbd2139aae..22034ee696 100644 --- a/src/nrniv/kssingle.h +++ b/src/nrniv/kssingle.h @@ -5,7 +5,8 @@ #include #include #include -#include + +#include class KSSingleTrans; class KSSingleState; @@ -70,10 +71,10 @@ class KSSingle { return MyMath::eq(x, y, vres_); } double exprand() { - return nrnran123_negexp(rand_state); + return -log(mcell_ran4a(&idum_)); } double unifrand(double range) { - return nrnran123_dblpick(rand_state) * range; + return mcell_ran4a(&idum_) * range; } int rvalrand(int); @@ -83,7 +84,7 @@ class KSSingle { double* rval_; bool uses_ligands_; static double vres_; - static nrnran123_State* rand_state; + static unsigned int idum_; static unsigned long singleevent_deliver_; static unsigned long singleevent_move_; diff --git a/src/oc/hoc_init.cpp b/src/oc/hoc_init.cpp index 455630a653..17675d7753 100644 --- a/src/oc/hoc_init.cpp +++ b/src/oc/hoc_init.cpp @@ -255,6 +255,7 @@ void hoc_init(void) /* install constants and built-ins table */ } } + set_use_mcran4(false); nrn_xopen_broadcast_ = 255; extern void hoc_init_space(void); hoc_init_space(); diff --git a/src/oc/oc_mcran4.cpp b/src/oc/oc_mcran4.cpp index 2bd1ef043d..128ccd7d04 100644 --- a/src/oc/oc_mcran4.cpp +++ b/src/oc/oc_mcran4.cpp @@ -1,7 +1,15 @@ #include "hocdec.h" -#include "nrnran123.h" +#include "mcran4.h" -static const nrnran123_State* state = nrnran123_newstream(0, 0); +int use_mcell_ran4_; + +void set_use_mcran4(bool value) { + use_mcell_ran4_ = value ? 1 : 0; +} + +bool use_mcran4() { + return use_mcell_ran4_ != 0; +} void hoc_mcran4() { uint32_t idx; @@ -9,25 +17,25 @@ void hoc_mcran4() { double x; xidx = hoc_pgetarg(1); idx = (uint32_t) (*xidx); - nrnran123_setseq(state, idx, 0); - x = nrnran123_dblpick(state); - *xidx = idx + 1; + x = mcell_ran4a(&idx); + *xidx = idx; hoc_ret(); hoc_pushx(x); } void hoc_mcran4init() { - std::uint32_t seq; - char which; - nrnran123_getseq(state, &seq, &which); - double prev = static_cast(seq); + double prev = mcell_lowindex(); if (ifarg(1)) { uint32_t idx = (uint32_t) chkarg(1, 0., 4294967295.); - nrnran123_setseq(state, idx, 0); + mcell_ran4_init(idx); } hoc_ret(); hoc_pushx(prev); } void hoc_usemcran4() { + double prev = (double) use_mcell_ran4_; + if (ifarg(1)) { + use_mcell_ran4_ = (int) chkarg(1, 0., 1.); + } hoc_ret(); - hoc_pushx(0); + hoc_pushx(prev); } diff --git a/src/oc/oc_mcran4.hpp b/src/oc/oc_mcran4.hpp index 1dca3b6feb..d275e3e9bf 100644 --- a/src/oc/oc_mcran4.hpp +++ b/src/oc/oc_mcran4.hpp @@ -1,3 +1,5 @@ +void set_use_mcran4(bool value); +bool use_mcran4(); void hoc_mcran4(); void hoc_mcran4init(); void hoc_usemcran4(); diff --git a/src/oc/scoprand.cpp b/src/oc/scoprand.cpp index 37e55ee789..70651f5608 100644 --- a/src/oc/scoprand.cpp +++ b/src/oc/scoprand.cpp @@ -23,7 +23,8 @@ static char RCSid[] = "random.cpp,v 1.4 1999/01/04 12:46:49 hines Exp"; #endif #include -#include "nrnran123.h" +#include "oc_mcran4.hpp" +#include "mcran4.h" #include "scoplib.h" static uint32_t value = 1; @@ -58,8 +59,17 @@ static uint32_t value = 1; *--------------------------------------------------------------------------- */ double scop_random(void) { - auto *rng = nrnran123_newstream(value, 0); - return nrnran123_dblpick(rng); + if (use_mcran4()) { + /*perhaps 4 times slower but much higher quality*/ + return mcell_ran4a(&value); + } else { + uint32_t a = 2147437301, c = 453816981, + /* m = 2^32 - 1, the largest long int value that can be represented */ + /*m = 0xFFFFFFFF;*/ /* limited to 32 bit integers*/ + m = ~0; + value = a * value + c; + return (fabs((double) value / (double) m)); + } } /*----------------------------------------------------------------------------- diff --git a/src/scopmath/praxis.cpp b/src/scopmath/praxis.cpp index daea1160bd..effb3f8844 100644 --- a/src/scopmath/praxis.cpp +++ b/src/scopmath/praxis.cpp @@ -1,13 +1,13 @@ #include <../../nrnconf.h> +#include "mcran4.h" +#include "oc_ansi.h" +#include "scoplib.h" + #include #include #include #include -#include "oc_ansi.h" -#include "scoplib.h" -#include - #undef small typedef long int integer; @@ -1730,7 +1730,55 @@ static /* Subroutine */ int maprnt_(integer* option, doublereal* v, integer* m, uint32_t nrn_praxis_ran_index; -static doublereal random_(integer* /* naught */) { - auto *rng = nrnran123_newstream(nrn_praxis_ran_index, 0); - return nrnran123_dblpick(rng); -} +static doublereal random_(integer* naught) { + double x; + return mcell_ran4(&nrn_praxis_ran_index, &x, 1, 1.); +#if 0 + /* Initialized data */ + + static logical init = FALSE_; + + /* System generated locals */ + doublereal ret_val; + + /* Local variables */ + static doublereal half; + static integer i, j, q, r; + static doublereal ran1; + static integer ran2; + static doublereal ran3[127]; + + if (init) { + goto L3; + } + r = *naught % 8190 + 1; + ran2 = 128; + for (i = 1; i <= 127; ++i) { + --ran2; + ran1 = -36028797018963968.; + for (j = 1; j <= 7; ++j) { + r = r * 1756 % 8191; + q = r / 32; +/* L1: */ + ran1 = (ran1 + q) * .00390625; + } +/* L2: */ + ran3[ran2 - 1] = ran1; + } + init = TRUE_; +L3: + if (ran2 == 1) { + ran2 = 128; + } + --ran2; + ran1 += ran3[ran2 - 1]; + half = .5; + if (ran1 >= 0.) { + half = -half; + } + ran1 += half; + ran3[ran2 - 1] = ran1; + ret_val = ran1 + .5; + return ret_val; +#endif +} /* random_ */ From 84c7da1fcde37c18c3609230687ca826f821a2ad Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 1 Feb 2024 11:26:51 +0100 Subject: [PATCH 11/12] Fix --- src/gnu/CMakeLists.txt | 1 + src/gnu/MCellRan4RNG.cpp | 15 ++++++++ src/gnu/MCellRan4RNG.hpp | 34 ++++++++++++++++++ src/gnu/NrnRandom123RNG.hpp | 6 ++-- src/gnu/RNG.h | 8 ++--- src/gnu/Rand.cpp | 1 + src/gnu/Rand.hpp | 8 +++++ src/ivoc/ivocrand.cpp | 70 +++++++++++++++++++++++++++++-------- 8 files changed, 120 insertions(+), 23 deletions(-) create mode 100644 src/gnu/MCellRan4RNG.cpp create mode 100644 src/gnu/MCellRan4RNG.hpp diff --git a/src/gnu/CMakeLists.txt b/src/gnu/CMakeLists.txt index 1d8783a621..c8b5f59cdb 100644 --- a/src/gnu/CMakeLists.txt +++ b/src/gnu/CMakeLists.txt @@ -6,6 +6,7 @@ add_library( Geom.cpp HypGeom.cpp LogNorm.cpp + MCellRan4RNG.cpp mcran4.cpp NegExp.cpp Normal.cpp diff --git a/src/gnu/MCellRan4RNG.cpp b/src/gnu/MCellRan4RNG.cpp new file mode 100644 index 0000000000..43b732cbe1 --- /dev/null +++ b/src/gnu/MCellRan4RNG.cpp @@ -0,0 +1,15 @@ +#include "MCellRan4RNG.hpp" + +MCellRan4::MCellRan4(std::uint32_t ihigh, std::uint32_t ilow) { + ++cnt_; + ilow_ = ilow; + ihigh_ = ihigh; + if (ihigh_ == 0) { + ihigh_ = cnt_; + ihigh_ = (std::uint32_t) asLong(); + } + orig_ = ihigh_; +} +MCellRan4::~MCellRan4() {} + +std::uint32_t MCellRan4::cnt_ = 0; diff --git a/src/gnu/MCellRan4RNG.hpp b/src/gnu/MCellRan4RNG.hpp new file mode 100644 index 0000000000..a4356562a9 --- /dev/null +++ b/src/gnu/MCellRan4RNG.hpp @@ -0,0 +1,34 @@ +#pragma once + +#include + +#include "RNG.h" +#include "mcran4.h" + +// The decision that has to be made is whether each generator instance +// should have its own seed or only one seed for all. We choose separate +// seed for each but if arg not present or 0 then seed chosen by system. + +// the addition of ilow > 0 means that value is used for the lowindex +// instead of the mcell_ran4_init global 32 bit lowindex. + +class MCellRan4: public RNG { + public: + MCellRan4(std::uint32_t ihigh = 0, std::uint32_t ilow = 0); + virtual ~MCellRan4(); + virtual std::uint32_t asLong() { + return (std::uint32_t) (ilow_ == 0 ? mcell_iran4(&ihigh_) : nrnRan4int(&ihigh_, ilow_)); + } + virtual void reset() { + ihigh_ = orig_; + } + virtual double asDouble() { + return (ilow_ == 0 ? mcell_ran4a(&ihigh_) : nrnRan4dbl(&ihigh_, ilow_)); + } + std::uint32_t ihigh_; + std::uint32_t orig_; + std::uint32_t ilow_; + + private: + static std::uint32_t cnt_; +}; diff --git a/src/gnu/NrnRandom123RNG.hpp b/src/gnu/NrnRandom123RNG.hpp index d4484d45f8..d68dae4832 100644 --- a/src/gnu/NrnRandom123RNG.hpp +++ b/src/gnu/NrnRandom123RNG.hpp @@ -9,9 +9,9 @@ class NrnRandom123: public RNG { public: NrnRandom123(std::uint32_t id1, std::uint32_t id2, std::uint32_t id3 = 0); ~NrnRandom123(); - std::uint32_t asLong(); - double asDouble(); - void reset(); + std::uint32_t asLong() override; + double asDouble() override; + void reset() override; std::uint32_t get_seq() const; void set_seq(std::uint32_t seq); nrnran123_State* s_; diff --git a/src/gnu/RNG.h b/src/gnu/RNG.h index 578be4bd85..38b6eebb86 100755 --- a/src/gnu/RNG.h +++ b/src/gnu/RNG.h @@ -27,17 +27,15 @@ Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. * use these in place of the ones that libg++ used to provide. */ #include -#include -#include union PrivateRNGSingleType { // used to access floats as unsigneds float s; - uint32_t u; + std::uint32_t u; }; union PrivateRNGDoubleType { // used to access doubles as unsigneds double d; - uint32_t u[2]; + std::uint32_t u[2]; }; // @@ -52,7 +50,7 @@ class RNG { // // Return a long-words word of random bits // - virtual uint32_t asLong() = 0; + virtual std::uint32_t asLong() = 0; virtual void reset() = 0; // // Return random bits converted to either a float or a double diff --git a/src/gnu/Rand.cpp b/src/gnu/Rand.cpp index 195a415a75..f0289c0f37 100644 --- a/src/gnu/Rand.cpp +++ b/src/gnu/Rand.cpp @@ -7,6 +7,7 @@ Rand::Rand(unsigned long seed, int size, Object* obj) { // printf("Rand\n"); 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 010e435898..008ab39d1d 100644 --- a/src/gnu/Rand.hpp +++ b/src/gnu/Rand.hpp @@ -5,11 +5,19 @@ struct Object; +/* type_: + * 0: ACG (unused) + * 1: MLCG (unused) + * 2: MCellRan4 + * 3: Isaac64 (unused) + * 4: Random123 + */ class Rand { public: Rand(unsigned long seed = 0, int size = 55, Object* obj = nullptr); ~Rand(); RNG* gen; Random* rand; + int type_; Object* obj_; }; diff --git a/src/ivoc/ivocrand.cpp b/src/ivoc/ivocrand.cpp index c44455864b..f30d4e3b0d 100644 --- a/src/ivoc/ivocrand.cpp +++ b/src/ivoc/ivocrand.cpp @@ -30,6 +30,7 @@ #include #include #include +#include #if HAVE_IV #include "ivoc.h" @@ -129,6 +130,7 @@ static double r_nrnran123(void* r) { } delete x->gen; x->gen = x->rand->generator(); + x->type_ = 4; return 0.; } @@ -148,7 +150,23 @@ static double r_MLCG(void* r) { } static double r_MCellRan4(void* r) { - return r_nrnran123(r); + auto* x = static_cast(r); + + std::uint32_t seed1 = 0; + std::uint32_t ilow = 0; + + if (ifarg(1)) { + seed1 = static_cast(chkarg(1, 0., dmaxuint)); + } + if (ifarg(2)) { + ilow = static_cast(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 static_cast(mcr->orig_); } static double r_Isaac64(void* r) { @@ -157,16 +175,25 @@ static double r_Isaac64(void* r) { long nrn_get_random_sequence(Rand* r) { - auto* gen = static_cast(r->gen); - return gen->get_seq(); + nrn_assert(r->type_ == 2); + auto* mcr = static_cast(r->gen); + return mcr->ihigh_; } + void nrn_set_random_sequence(Rand* r, long seq) { - auto* gen = static_cast(r->gen); - gen->set_seq(seq); + nrn_assert(r->type_ == 2); + auto* mcr = static_cast(r->gen); + mcr->ihigh_ = seq; } int nrn_random_isran123(Rand* r, uint32_t* id1, uint32_t* id2, uint32_t* id3) { + if (r->type_ != 4) { + return 0; + } + + auto* nr = static_cast(r->gen); + nrnran123_getids3(nr->s_, id1, id2, id3); return 1; } @@ -179,26 +206,39 @@ static double r_ran123_globalindex(void* r) { } static double r_sequence(void* r) { - Rand* x = (Rand*) r; - uint32_t seq; - char which; + auto* x = static_cast(r); + if (x->type_ == 4) { + std::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); + } + auto* mcr = static_cast(x->gen); 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); + mcr->ihigh_ = static_cast(*getarg(1)); } - nrnran123_getseq(((NrnRandom123*) x->gen)->s_, &seq, &which); - return double(seq) * 4. + double(which); + return static_cast(mcr->ihigh_); } 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; } From efbf7c570fc852add83497e1eb93dd9551b1ab20 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 1 Feb 2024 11:46:02 +0100 Subject: [PATCH 12/12] Fix tests --- test/pytest_coreneuron/test_nlayer.py | 10 +++++++--- test/pytest_coreneuron/test_vector_api.py | 8 ++++---- 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/test/pytest_coreneuron/test_nlayer.py b/test/pytest_coreneuron/test_nlayer.py index 5a7670da58..fed680f7e5 100644 --- a/test/pytest_coreneuron/test_nlayer.py +++ b/test/pytest_coreneuron/test_nlayer.py @@ -36,7 +36,7 @@ def run(self): return res -def test_nlayer(): +def get_result(): res = [] for nlayer in [1, 2, 3]: h.nlayer_extracellular(nlayer) @@ -44,16 +44,20 @@ def test_nlayer(): m = Model() res.append(m.run()) del m + return res + + +def test_nlayer(): + res = get_result() for r in res[1:]: for i in range(2): x = r[i].c().sub(res[0][i]).abs().sum() print(x) assert x < 0.01 # At 10ms. At 5ms would be 2e-5 - return res if __name__ == "__main__": - res = test_nlayer() + res = get_result() for i in range(len(res[0][0])): for j in range(2): for k in range(3): diff --git a/test/pytest_coreneuron/test_vector_api.py b/test/pytest_coreneuron/test_vector_api.py index ff90bafb2d..faaac86636 100644 --- a/test/pytest_coreneuron/test_vector_api.py +++ b/test/pytest_coreneuron/test_vector_api.py @@ -514,10 +514,10 @@ def test_vector_api(): vrand = h.Vector((1, 2, 3)) r = h.Random() r.poisson(12) - assert vrand.cl().setrand(r).to_python() == [10.0, 16.0, 11.0] - assert vrand.cl().setrand(r, 1, 2).to_python() == [1.0, 9.0, 18.0] - assert vrand.cl().addrand(r).to_python() == [9.0, 9.0, 16.0] - assert vrand.cl().addrand(r, 0, 1).to_python() == [13.0, 16.0, 3.0] + assert vrand.cl().setrand(r).to_python() == [22.0, 16.0, 16.0] + assert vrand.cl().setrand(r, 1, 2).to_python() == [1.0, 8.0, 12.0] + assert vrand.cl().addrand(r).to_python() == [7.0, 13.0, 12.0] + assert vrand.cl().addrand(r, 0, 1).to_python() == [9.0, 16.0, 3.0] """ Misc