From f64b609dca40ee1c84e7e059ea45ca2d4560085e Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 25 Jan 2024 23:35:48 +0100 Subject: [PATCH 01/15] Modernize nrnran123 (#2685) --- src/ivoc/ivocrand.cpp | 8 ++- src/oc/nrnran123.cpp | 117 ++++++++++++++++++++---------------------- src/oc/nrnran123.h | 49 ++++++++++-------- 3 files changed, 92 insertions(+), 82 deletions(-) diff --git a/src/ivoc/ivocrand.cpp b/src/ivoc/ivocrand.cpp index d3501913db..fa45403dcc 100644 --- a/src/ivoc/ivocrand.cpp +++ b/src/ivoc/ivocrand.cpp @@ -337,8 +337,12 @@ static double r_nrnran123(void* r) { id2 = (uint32_t) (chkarg(2, 0., dmaxuint)); if (ifarg(3)) id3 = (uint32_t) (chkarg(3, 0., dmaxuint)); - NrnRandom123* r123 = new NrnRandom123(id1, id2, id3); - x->rand->generator(r123); + 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; diff --git a/src/oc/nrnran123.cpp b/src/oc/nrnran123.cpp index d9db49351c..eeb2768115 100644 --- a/src/oc/nrnran123.cpp +++ b/src/oc/nrnran123.cpp @@ -1,81 +1,79 @@ -#include <../../nrnconf.h> -#include -#include -#include -#include -#include +#include +#include "nrnran123.h" +#include +#include #include -static const double SHIFT32 = 1.0 / 4294967297.0; /* 1/(2^32 + 1) */ +using RNG = r123::Philox4x32; -static philox4x32_key_t k = {{0}}; +static RNG::key_type k = {{0}}; struct nrnran123_State { - philox4x32_ctr_t c; - philox4x32_ctr_t r; + RNG::ctr_type c; + RNG::ctr_type r; char which_; }; -void nrnran123_set_globalindex(uint32_t gix) { - k.v[0] = gix; +void nrnran123_set_globalindex(std::uint32_t gix) { + k[0] = gix; } /* if one sets the global, one should reset all the stream sequences. */ -uint32_t nrnran123_get_globalindex() { - return k.v[0]; +std::uint32_t nrnran123_get_globalindex() { + return k[0]; } -nrnran123_State* nrnran123_newstream(uint32_t id1, uint32_t id2) { +nrnran123_State* nrnran123_newstream(std::uint32_t id1, std::uint32_t id2) { return nrnran123_newstream3(id1, id2, 0); } -nrnran123_State* nrnran123_newstream3(uint32_t id1, uint32_t id2, uint32_t id3) { - nrnran123_State* s; - s = (nrnran123_State*) ecalloc(sizeof(nrnran123_State), 1); - s->c.v[1] = id3; - s->c.v[2] = id1; - s->c.v[3] = id2; +nrnran123_State* nrnran123_newstream3(std::uint32_t id1, std::uint32_t id2, std::uint32_t id3) { + auto* s = new nrnran123_State; + s->c[1] = id3; + s->c[2] = id1; + s->c[3] = id2; nrnran123_setseq(s, 0, 0); return s; } void nrnran123_deletestream(nrnran123_State* s) { - free(s); + delete s; } -void nrnran123_getseq(nrnran123_State* s, uint32_t* seq, char* which) { - *seq = s->c.v[0]; +void nrnran123_getseq(nrnran123_State* s, std::uint32_t* seq, char* which) { + *seq = s->c[0]; *which = s->which_; } -void nrnran123_setseq(nrnran123_State* s, uint32_t seq, char which) { +void nrnran123_setseq(nrnran123_State* s, std::uint32_t seq, char which) { if (which > 3 || which < 0) { s->which_ = 0; } else { s->which_ = which; } - s->c.v[0] = seq; + s->c[0] = seq; s->r = philox4x32(s->c, k); } -void nrnran123_getids(nrnran123_State* s, uint32_t* id1, uint32_t* id2) { - *id1 = s->c.v[2]; - *id2 = s->c.v[3]; +void nrnran123_getids(nrnran123_State* s, std::uint32_t* id1, std::uint32_t* id2) { + *id1 = s->c[2]; + *id2 = s->c[3]; } -void nrnran123_getids3(nrnran123_State* s, uint32_t* id1, uint32_t* id2, uint32_t* id3) { - *id3 = s->c.v[1]; - *id1 = s->c.v[2]; - *id2 = s->c.v[3]; +void nrnran123_getids3(nrnran123_State* s, + std::uint32_t* id1, + std::uint32_t* id2, + std::uint32_t* id3) { + *id3 = s->c[1]; + *id1 = s->c[2]; + *id2 = s->c[3]; } -uint32_t nrnran123_ipick(nrnran123_State* s) { - uint32_t rval; +std::uint32_t nrnran123_ipick(nrnran123_State* s) { char which = s->which_; - assert(which < 4); - rval = s->r.v[which++]; + std::uint32_t rval = s->r[which++]; if (which > 3) { which = 0; - s->c.v[0]++; + s->c.incr(); s->r = philox4x32(s->c, k); } s->which_ = which; @@ -83,15 +81,17 @@ uint32_t nrnran123_ipick(nrnran123_State* s) { } double nrnran123_dblpick(nrnran123_State* s) { - return nrnran123_uint2dbl(nrnran123_ipick(s)); + static const double SHIFT32 = 1.0 / 4294967297.0; /* 1/(2^32 + 1) */ + auto u = nrnran123_ipick(s); + return ((double) u + 1.0) * SHIFT32; } double nrnran123_negexp(nrnran123_State* s) { /* min 2.3283064e-10 to max 22.18071 */ - return -log(nrnran123_dblpick(s)); + return -std::log(nrnran123_dblpick(s)); } -/* At cost of a cached value we could compute two at a time. */ +/* At cost of a cached value we could compute two at a time. */ /* But that would make it difficult to transfer to coreneuron for t > 0 */ double nrnran123_normal(nrnran123_State* s) { double w, x, y; @@ -104,31 +104,28 @@ double nrnran123_normal(nrnran123_State* s) { w = (u1 * u1) + (u2 * u2); } while (w > 1); - y = sqrt((-2. * log(w)) / w); + y = std::sqrt((-2. * log(w)) / w); x = u1 * y; return x; } -nrnran123_array4x32 nrnran123_iran(uint32_t seq, uint32_t id1, uint32_t id2) { +nrnran123_array4x32 nrnran123_iran(std::uint32_t seq, std::uint32_t id1, std::uint32_t id2) { return nrnran123_iran3(seq, id1, id2, 0); } -nrnran123_array4x32 nrnran123_iran3(uint32_t seq, uint32_t id1, uint32_t id2, uint32_t id3) { +nrnran123_array4x32 nrnran123_iran3(std::uint32_t seq, + std::uint32_t id1, + std::uint32_t id2, + std::uint32_t id3) { nrnran123_array4x32 a; - philox4x32_ctr_t c; - c.v[0] = seq; - c.v[1] = id3; - c.v[2] = id1; - c.v[3] = id2; - philox4x32_ctr_t r = philox4x32(c, k); - a.v[0] = r.v[0]; - a.v[1] = r.v[1]; - a.v[2] = r.v[2]; - a.v[3] = r.v[3]; + RNG::ctr_type c; + c[0] = seq; + c[1] = id3; + c[2] = id1; + c[3] = id2; + RNG::ctr_type r = philox4x32(c, k); + a.v[0] = r[0]; + a.v[1] = r[1]; + a.v[2] = r[2]; + a.v[3] = r[3]; return a; } - -double nrnran123_uint2dbl(uint32_t u) { - /* 0 to 2^32-1 transforms to double value in open (0,1) interval */ - /* min 2.3283064e-10 to max (1 - 2.3283064e-10) */ - return ((double) u + 1.0) * SHIFT32; -} diff --git a/src/oc/nrnran123.h b/src/oc/nrnran123.h index 0a0acb6a71..a4e4844bf6 100644 --- a/src/oc/nrnran123.h +++ b/src/oc/nrnran123.h @@ -23,40 +23,49 @@ of the full distribution available from http://www.deshawresearch.com/resources_random123.html */ -#include +#include -typedef struct nrnran123_State nrnran123_State; +struct nrnran123_State; -typedef struct nrnran123_array4x32 { - uint32_t v[4]; -} nrnran123_array4x32; +struct nrnran123_array4x32 { + std::uint32_t v[4]; +}; /* global index. eg. run number */ /* all generator instances share this global index */ -extern void nrnran123_set_globalindex(uint32_t gix); -extern uint32_t nrnran123_get_globalindex(); +extern void nrnran123_set_globalindex(std::uint32_t gix); +extern std::uint32_t nrnran123_get_globalindex(); /* minimal data stream */ -extern nrnran123_State* nrnran123_newstream(uint32_t id1, uint32_t id2); -extern nrnran123_State* nrnran123_newstream3(uint32_t id1, uint32_t id2, uint32_t id3); +extern nrnran123_State* nrnran123_newstream(std::uint32_t id1, std::uint32_t id2); +extern nrnran123_State* nrnran123_newstream3(std::uint32_t id1, + std::uint32_t id2, + std::uint32_t id3); extern void nrnran123_deletestream(nrnran123_State*); -extern void nrnran123_getseq(nrnran123_State*, uint32_t* seq, char* which); -extern void nrnran123_setseq(nrnran123_State*, uint32_t seq, char which); -extern void nrnran123_getids(nrnran123_State*, uint32_t* id1, uint32_t* id2); -extern void nrnran123_getids3(nrnran123_State*, uint32_t* id1, uint32_t* id2, uint32_t* id3); +extern void nrnran123_getseq(nrnran123_State*, std::uint32_t* seq, char* which); +extern void nrnran123_setseq(nrnran123_State*, std::uint32_t seq, char which); +extern void nrnran123_getids(nrnran123_State*, std::uint32_t* id1, std::uint32_t* id2); +extern void nrnran123_getids3(nrnran123_State*, + std::uint32_t* id1, + std::uint32_t* id2, + std::uint32_t* id3); -extern double nrnran123_negexp(nrnran123_State*); /* mean 1.0 */ -extern uint32_t nrnran123_ipick(nrnran123_State*); /* uniform 0 to 2^32-1 */ -extern double nrnran123_dblpick(nrnran123_State*); /* uniform open interval (0,1)*/ -/* nrnran123_dblpick minimum value is 2.3283064e-10 and max value is 1-min */ +// Get a random uint32_t in [0, 2^32-1] +extern std::uint32_t nrnran123_ipick(nrnran123_State*); +// Get a random double on [0, 1] +// nrnran123_dblpick minimum value is 2.3283064e-10 and max value is 1-min +extern double nrnran123_dblpick(nrnran123_State*); /* nrnran123_negexp min value is 2.3283064e-10, max is 22.18071 */ +extern double nrnran123_negexp(nrnran123_State*); /* mean 1.0 */ extern double nrnran123_normal(nrnran123_State*); /* mean 0.0, std 1.0 */ /* more fundamental (stateless) (though the global index is still used) */ -extern nrnran123_array4x32 nrnran123_iran(uint32_t seq, uint32_t id1, uint32_t id2); -extern nrnran123_array4x32 nrnran123_iran3(uint32_t seq, uint32_t id1, uint32_t id2, uint32_t id3); -extern double nrnran123_uint2dbl(uint32_t); +extern nrnran123_array4x32 nrnran123_iran(std::uint32_t seq, std::uint32_t id1, std::uint32_t id2); +extern nrnran123_array4x32 nrnran123_iran3(std::uint32_t seq, + std::uint32_t id1, + std::uint32_t id2, + std::uint32_t id3); #endif From 08920e75e7b65f1c3731b94356f7277bb1b12264 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Sat, 27 Jan 2024 09:22:30 +0100 Subject: [PATCH 02/15] Move random number related code to src/gnu (#2689) * This patch moves most of the code that is not HOC to the directory gnu. * `Isaac`, `MCellRan4`, `Random123`, `Rand` * In the future this directory will be cleaned to remove gnu stuff and renamed something like 'src/random'. * This way nrn will have an external random library, and a clean hoc_interface. --- cmake/CompilerHelper.cmake | 2 +- cmake/NeuronFileLists.cmake | 13 +-- src/gnu/CMakeLists.txt | 9 ++ src/gnu/Isaac64RNG.cpp | 20 ++++ src/gnu/Isaac64RNG.hpp | 33 ++++++ src/gnu/MCellRan4RNG.cpp | 15 +++ src/gnu/MCellRan4RNG.hpp | 34 ++++++ src/gnu/NrnRandom123RNG.cpp | 0 src/gnu/NrnRandom123RNG.hpp | 28 +++++ src/gnu/Rand.cpp | 19 ++++ src/{ivoc/random1.h => gnu/Rand.hpp} | 14 ++- src/{oc => gnu}/isaac64.cpp | 5 - src/{oc => gnu}/isaac64.h | 7 +- src/{oc => gnu}/mcran4.cpp | 42 +------ src/{oc => gnu}/mcran4.h | 1 + src/gnu/nrnisaac.cpp | 25 ++++ src/gnu/nrnisaac.h | 9 ++ src/{oc => gnu}/nrnran123.cpp | 0 src/{oc => gnu}/nrnran123.h | 0 src/ivoc/ivocrand.cpp | 163 ++++----------------------- src/ivoc/ivocvect.cpp | 2 +- src/nrniv/CMakeLists.txt | 1 - src/oc/hoc_init.cpp | 4 +- src/oc/nrnisaac.cpp | 39 ------- src/oc/nrnisaac.h | 15 --- src/oc/oc_mcran4.cpp | 41 +++++++ src/oc/oc_mcran4.hpp | 5 + src/oc/ocfunc.h | 1 - src/oc/scoprand.cpp | 6 +- 29 files changed, 286 insertions(+), 267 deletions(-) create mode 100644 src/gnu/Isaac64RNG.cpp create mode 100644 src/gnu/Isaac64RNG.hpp create mode 100644 src/gnu/MCellRan4RNG.cpp create mode 100644 src/gnu/MCellRan4RNG.hpp create mode 100644 src/gnu/NrnRandom123RNG.cpp create mode 100644 src/gnu/NrnRandom123RNG.hpp create mode 100644 src/gnu/Rand.cpp rename src/{ivoc/random1.h => gnu/Rand.hpp} (63%) rename src/{oc => gnu}/isaac64.cpp (98%) rename src/{oc => gnu}/isaac64.h (97%) rename src/{oc => gnu}/mcran4.cpp (86%) rename src/{oc => gnu}/mcran4.h (98%) create mode 100644 src/gnu/nrnisaac.cpp create mode 100644 src/gnu/nrnisaac.h rename src/{oc => gnu}/nrnran123.cpp (100%) rename src/{oc => gnu}/nrnran123.h (100%) delete mode 100644 src/oc/nrnisaac.cpp delete mode 100644 src/oc/nrnisaac.h create mode 100644 src/oc/oc_mcran4.cpp create mode 100644 src/oc/oc_mcran4.hpp diff --git a/cmake/CompilerHelper.cmake b/cmake/CompilerHelper.cmake index 54adf1ce31..e50b120351 100644 --- a/cmake/CompilerHelper.cmake +++ b/cmake/CompilerHelper.cmake @@ -51,7 +51,7 @@ if(CMAKE_C_COMPILER_ID MATCHES "PGI" OR CMAKE_C_COMPILER_ID MATCHES "NVHPC") # Random123 does not play nicely with NVHPC 21.11+'s detection of ABM features, see: # https://github.com/BlueBrain/CoreNeuron/issues/724 and # https://github.com/DEShawResearch/random123/issues/6. - list(APPEND NRN_COMPILE_DEFS R123_USE_INTRIN_H=0) + list(APPEND NRN_R123_COMPILE_DEFS R123_USE_INTRIN_H=0) endif() else() set(NRN_HAVE_NVHPC_COMPILER OFF) diff --git a/cmake/NeuronFileLists.cmake b/cmake/NeuronFileLists.cmake index 252f3e9ff6..c1beffb9fd 100644 --- a/cmake/NeuronFileLists.cmake +++ b/cmake/NeuronFileLists.cmake @@ -9,6 +9,9 @@ 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/nrnisaac.h + gnu/nrnran123.h nrniv/backtrace_utils.h nrniv/bbsavestate.h nrnmpi/nrnmpidec.h @@ -37,15 +40,12 @@ set(HEADER_FILES_TO_INSTALL oc/hocgetsym.h oc/hoclist.h oc/hocparse.h - oc/mcran4.h oc/mech_api.h oc/memory.hpp oc/nrnapi.h oc/nrnassrt.h - oc/nrnisaac.h oc/nrnmpi.h oc/nrnrandom.h - oc/nrnran123.h oc/oc_ansi.h oc/ocfunc.h oc/ocmisc.h @@ -77,16 +77,10 @@ set(HEADER_FILES_TO_INSTALL # ============================================================================= set(NRN_HEADERS_INCLUDE_LIST) -# ============================================================================= -# Lists of random number related files -# ============================================================================= -set(RAN_FILE_LIST isaac64.cpp mcran4.cpp nrnisaac.cpp nrnran123.cpp) - # ============================================================================= # Files in oc directory # ============================================================================= set(OC_FILE_LIST - ${RAN_FILE_LIST} audit.cpp axis.cpp code.cpp @@ -102,6 +96,7 @@ set(OC_FILE_LIST hoc_oop.cpp list.cpp math.cpp + oc_mcran4.cpp memory.cpp mswinprt.cpp nonlin.cpp diff --git a/src/gnu/CMakeLists.txt b/src/gnu/CMakeLists.txt index 0e9562297c..e16ae608b9 100644 --- a/src/gnu/CMakeLists.txt +++ b/src/gnu/CMakeLists.txt @@ -6,14 +6,23 @@ add_library( Erlang.cpp Geom.cpp HypGeom.cpp + isaac64.cpp + Isaac64RNG.cpp LogNorm.cpp + MCellRan4RNG.cpp + mcran4.cpp MLCG.cpp NegExp.cpp Normal.cpp + nrnisaac.cpp + nrnran123.cpp Poisson.cpp + Rand.cpp Random.cpp RndInt.cpp RNG.cpp Uniform.cpp Weibull.cpp) set_property(TARGET nrngnu PROPERTY POSITION_INDEPENDENT_CODE ON) +target_include_directories(nrngnu SYSTEM PRIVATE "${PROJECT_SOURCE_DIR}/external/Random123/include") +target_compile_definitions(nrngnu PRIVATE ${NRN_R123_COMPILE_DEFS}) diff --git a/src/gnu/Isaac64RNG.cpp b/src/gnu/Isaac64RNG.cpp new file mode 100644 index 0000000000..559c307d4d --- /dev/null +++ b/src/gnu/Isaac64RNG.cpp @@ -0,0 +1,20 @@ +#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 new file mode 100644 index 0000000000..c825c741ca --- /dev/null +++ b/src/gnu/Isaac64RNG.hpp @@ -0,0 +1,33 @@ +#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 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.cpp b/src/gnu/NrnRandom123RNG.cpp new file mode 100644 index 0000000000..e69de29bb2 diff --git a/src/gnu/NrnRandom123RNG.hpp b/src/gnu/NrnRandom123RNG.hpp new file mode 100644 index 0000000000..5e790c939c --- /dev/null +++ b/src/gnu/NrnRandom123RNG.hpp @@ -0,0 +1,28 @@ +#pragma once + +#include + +#include "nrnran123.h" +#include "RNG.h" + +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); + } + 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 new file mode 100644 index 0000000000..a3b526de5e --- /dev/null +++ b/src/gnu/Rand.cpp @@ -0,0 +1,19 @@ +#include "Rand.hpp" + +#include "ACG.h" +#include "Normal.h" + +Rand::Rand(unsigned long seed, int size, Object* obj) { + // printf("Rand\n"); + gen = new ACG(seed, size); + rand = new Normal(0., 1., gen); + type_ = 0; + obj_ = obj; +} + +Rand::~Rand() { + // printf("~Rand\n"); + delete gen; + delete rand; +} + diff --git a/src/ivoc/random1.h b/src/gnu/Rand.hpp similarity index 63% rename from src/ivoc/random1.h rename to src/gnu/Rand.hpp index 8fedd8f5b0..f5b47a6a61 100644 --- a/src/ivoc/random1.h +++ b/src/gnu/Rand.hpp @@ -1,14 +1,20 @@ -#ifndef random1_h -#define random1_h +#pragma once #include "RNG.h" #include "Random.h" 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 = NULL); + Rand(unsigned long seed = 0, int size = 55, Object* obj = nullptr); ~Rand(); RNG* gen; Random* rand; @@ -16,5 +22,3 @@ class Rand { // double* looks like random variable that gets changed on every fadvance Object* obj_; }; - -#endif diff --git a/src/oc/isaac64.cpp b/src/gnu/isaac64.cpp similarity index 98% rename from src/oc/isaac64.cpp rename to src/gnu/isaac64.cpp index 579f52d4f0..b0c6c40317 100644 --- a/src/oc/isaac64.cpp +++ b/src/gnu/isaac64.cpp @@ -1,8 +1,3 @@ -#include <../../nrnconf.h> -#if HAVE_SYS_TYPES_H -#include -#endif - /* ------------------------------------------------------------------------------ isaac64.cpp: A Fast cryptographic random number generator diff --git a/src/oc/isaac64.h b/src/gnu/isaac64.h similarity index 97% rename from src/oc/isaac64.h rename to src/gnu/isaac64.h index 6b13dd72f1..4af3b7586e 100644 --- a/src/oc/isaac64.h +++ b/src/gnu/isaac64.h @@ -1,3 +1,5 @@ +#pragma once + /* ------------------------------------------------------------------------------ isaac64.h: Definitions for a fast cryptographic random number generator @@ -12,10 +14,7 @@ Jenkins, R.J. (1996) ISAAC, in Fast Software Encryption, vol. 1039, ------------------------------------------------------------------------------ */ -#ifndef ISAAC64_H -#define ISAAC64_H -#include #include #define RANDSIZL (4) /* I recommend 8 for crypto, 4 for simulations */ @@ -87,5 +86,3 @@ Macros to get individual random numbers rng->randcnt = RANDMAX - 2, \ DBL64 * (*((ub8*) (((ub4*) (rng->randrsl)) + rng->randcnt))))) - -#endif /* ISAAC64_H */ diff --git a/src/oc/mcran4.cpp b/src/gnu/mcran4.cpp similarity index 86% rename from src/oc/mcran4.cpp rename to src/gnu/mcran4.cpp index 1b303164d6..a9553687f6 100644 --- a/src/oc/mcran4.cpp +++ b/src/gnu/mcran4.cpp @@ -39,17 +39,19 @@ contained the header: Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ -#include <../../nrnconf.h> #include #include #include #include #include -#include -#include "hocdec.h" +#include "mcran4.h" static uint32_t lowindex = 0; +double mcell_lowindex() { + return static_cast(lowindex); +} + void mcell_ran4_init(uint32_t low) { lowindex = low; } @@ -70,40 +72,6 @@ uint32_t mcell_iran4(uint32_t* high) { return nrnRan4int(high, lowindex); } -/* Hoc interface */ -extern double chkarg(); -extern int use_mcell_ran4_; - - -void hoc_mcran4() { - uint32_t idx; - double* xidx; - double x; - xidx = hoc_pgetarg(1); - idx = (uint32_t) (*xidx); - x = mcell_ran4a(&idx); - *xidx = idx; - hoc_ret(); - hoc_pushx(x); -} -void hoc_mcran4init() { - double prev = (double) lowindex; - if (ifarg(1)) { - uint32_t idx = (uint32_t) chkarg(1, 0., 4294967295.); - 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(prev); -} - uint32_t nrnRan4int(uint32_t* idx1, uint32_t idx2) { uint32_t u, v, w, m, n; /* 64-bit hash */ diff --git a/src/oc/mcran4.h b/src/gnu/mcran4.h similarity index 98% rename from src/oc/mcran4.h rename to src/gnu/mcran4.h index 7d251b4cc1..9d7ff79dd9 100644 --- a/src/oc/mcran4.h +++ b/src/gnu/mcran4.h @@ -1,6 +1,7 @@ #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); diff --git a/src/gnu/nrnisaac.cpp b/src/gnu/nrnisaac.cpp new file mode 100644 index 0000000000..f0549a4d9a --- /dev/null +++ b/src/gnu/nrnisaac.cpp @@ -0,0 +1,25 @@ +#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 new file mode 100644 index 0000000000..5322febafe --- /dev/null +++ b/src/gnu/nrnisaac.h @@ -0,0 +1,9 @@ +#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/oc/nrnran123.cpp b/src/gnu/nrnran123.cpp similarity index 100% rename from src/oc/nrnran123.cpp rename to src/gnu/nrnran123.cpp diff --git a/src/oc/nrnran123.h b/src/gnu/nrnran123.h similarity index 100% rename from src/oc/nrnran123.h rename to src/gnu/nrnran123.h diff --git a/src/ivoc/ivocrand.cpp b/src/ivoc/ivocrand.cpp index fa45403dcc..f250049299 100644 --- a/src/ivoc/ivocrand.cpp +++ b/src/ivoc/ivocrand.cpp @@ -4,7 +4,7 @@ #include #include -#include "random1.h" +#include "Rand.hpp" #include #include "classreg.h" @@ -17,7 +17,6 @@ #include "ocobserv.h" #include -#include #include #include #include @@ -33,6 +32,9 @@ #include #include #include +#include +#include +#include #if HAVE_IV #include "ivoc.h" @@ -59,118 +61,6 @@ class RandomPlay: public Observer, public Resource { using RandomPlayList = std::vector; static RandomPlayList* random_play_list_; -#include - -class NrnRandom123: public RNG { - public: - NrnRandom123(uint32_t id1, uint32_t id2, uint32_t id3 = 0); - virtual ~NrnRandom123(); - virtual uint32_t asLong() { - return nrnran123_ipick(s_); - } - virtual double asDouble() { - return nrnran123_dblpick(s_); - } - virtual void reset() { - nrnran123_setseq(s_, 0, 0); - } - nrnran123_State* s_; -}; -NrnRandom123::NrnRandom123(uint32_t id1, uint32_t id2, uint32_t id3) { - s_ = nrnran123_newstream3(id1, id2, id3); -} -NrnRandom123::~NrnRandom123() { - nrnran123_deletestream(s_); -} - - -// 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(uint32_t ihigh = 0, uint32_t ilow = 0); - virtual ~MCellRan4(); - virtual uint32_t asLong() { - return (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_)); - } - uint32_t ihigh_; - uint32_t orig_; - uint32_t ilow_; - - private: - static uint32_t cnt_; -}; - -MCellRan4::MCellRan4(uint32_t ihigh, uint32_t ilow) { - ++cnt_; - ilow_ = ilow; - ihigh_ = ihigh; - if (ihigh_ == 0) { - ihigh_ = cnt_; - ihigh_ = (uint32_t) asLong(); - } - orig_ = ihigh_; -} -MCellRan4::~MCellRan4() {} - -uint32_t MCellRan4::cnt_ = 0; - -class Isaac64: public RNG { - public: - Isaac64(uint32_t seed = 0); - virtual ~Isaac64(); - virtual uint32_t asLong() { - return (uint32_t) nrnisaac_uint32_pick(rng_); - } - virtual void reset() { - nrnisaac_init(rng_, seed_); - } - virtual double asDouble() { - return nrnisaac_dbl_pick(rng_); - } - uint32_t seed() { - return seed_; - } - void seed(uint32_t s) { - seed_ = s; - reset(); - } - - private: - uint32_t seed_; - void* rng_; - static uint32_t cnt_; -}; - -Isaac64::Isaac64(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_); -} - -uint32_t Isaac64::cnt_ = 0; - RandomPlay::RandomPlay(Rand* r, neuron::container::data_handle px) : r_{r} , px_{std::move(px)} { @@ -197,27 +87,6 @@ void RandomPlay::update(Observable*) { list_remove(); } -Rand::Rand(unsigned long seed, int size, Object* obj) { - // printf("Rand\n"); - gen = new ACG(seed, size); - rand = new Normal(0., 1., gen); - type_ = 0; - obj_ = obj; -} - -Rand::~Rand() { - // printf("~Rand\n"); - delete gen; - delete rand; -} - -// constructor for a random number generator based on the RNG class -// from the gnu c++ class library -// defaults to the ACG generator (see below) - -// syntax: -// a = new Rand([seed],[size]) - static void* r_cons(Object* obj) { unsigned long seed = 0; int size = 55; @@ -406,14 +275,22 @@ static double r_Isaac64(void* r) { uint32_t seed1 = 0; - if (ifarg(1)) - seed1 = (uint32_t) (*getarg(1)); - Isaac64* mcr = new Isaac64(seed1); - x->rand->generator(mcr); - delete x->gen; - x->gen = x->rand->generator(); - x->type_ = 3; - return (double) mcr->seed(); + 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 diff --git a/src/ivoc/ivocvect.cpp b/src/ivoc/ivocvect.cpp index cc6565bea3..7b9087d328 100644 --- a/src/ivoc/ivocvect.cpp +++ b/src/ivoc/ivocvect.cpp @@ -111,7 +111,7 @@ static double dmaxint_ = 9007199254740992; #include "ivocvect.h" // definition of random numer generator -#include "random1.h" +#include "Rand.hpp" #include #if HAVE_IV diff --git a/src/nrniv/CMakeLists.txt b/src/nrniv/CMakeLists.txt index babfe56cc4..48f1c16531 100644 --- a/src/nrniv/CMakeLists.txt +++ b/src/nrniv/CMakeLists.txt @@ -211,7 +211,6 @@ set(NRN_INCLUDE_DIRS ${PROJECT_BINARY_DIR}/src/parallel ${PROJECT_BINARY_DIR}/src/sundials ${PROJECT_BINARY_DIR}/src/sundials/shared - ${PROJECT_SOURCE_DIR}/external/Random123/include ${PROJECT_SOURCE_DIR}/src ${PROJECT_SOURCE_DIR}/src/gnu ${PROJECT_SOURCE_DIR}/src/nrncvode diff --git a/src/oc/hoc_init.cpp b/src/oc/hoc_init.cpp index 6e022f7047..17675d7753 100644 --- a/src/oc/hoc_init.cpp +++ b/src/oc/hoc_init.cpp @@ -10,6 +10,7 @@ #include "nrn_ansi.h" #include "ocfunc.h" +#include "oc_mcran4.hpp" extern void hoc_nrnmpi_init(); @@ -231,7 +232,6 @@ double hoc_default_dll_loaded_; char* neuron_home; const char* nrn_mech_dll; /* but actually only for NEURON mswin and linux */ int nrn_noauto_dlopen_nrnmech; /* 0 except when binary special. */ -int use_mcell_ran4_; int nrn_xopen_broadcast_; void hoc_init(void) /* install constants and built-ins table */ @@ -255,7 +255,7 @@ void hoc_init(void) /* install constants and built-ins table */ } } - use_mcell_ran4_ = 0; + set_use_mcran4(false); nrn_xopen_broadcast_ = 255; extern void hoc_init_space(void); hoc_init_space(); diff --git a/src/oc/nrnisaac.cpp b/src/oc/nrnisaac.cpp deleted file mode 100644 index 1aa0050339..0000000000 --- a/src/oc/nrnisaac.cpp +++ /dev/null @@ -1,39 +0,0 @@ -#include <../../nrnconf.h> -#include -#if HAVE_SYS_TYPES_H -#include -#endif -#include -#include -#include "hocdec.h" - -typedef struct isaac64_state Rng; - -void* nrnisaac_new(void) { - Rng* rng; - rng = (Rng*) hoc_Emalloc(sizeof(Rng)); - hoc_malchk(); - return (void*) rng; -} - -void nrnisaac_delete(void* v) { - free(v); -} - -void nrnisaac_init(void* v, unsigned long int seed) { - isaac64_init((Rng*) v, seed); -} - -double nrnisaac_dbl_pick(void* v) { - Rng* rng = (Rng*) v; - double x = isaac64_dbl32(rng); - /*printf("dbl %d %d %d %d %g\n", sizeof(ub8), sizeof(ub4), sizeof(ub2), sizeof(ub1), x);*/ - return x; -} - -uint32_t nrnisaac_uint32_pick(void* v) { - Rng* rng = (Rng*) v; - double x = isaac64_uint32(rng); - /*printf("uint32 %g\n", x);*/ - return x; -} diff --git a/src/oc/nrnisaac.h b/src/oc/nrnisaac.h deleted file mode 100644 index 053366684c..0000000000 --- a/src/oc/nrnisaac.h +++ /dev/null @@ -1,15 +0,0 @@ -#ifndef nrnisaac_h -#define nrnisaac_h - -#include <../../nrnconf.h> -#include - - -void* nrnisaac_new(void); -void nrnisaac_delete(void* rng); -void nrnisaac_init(void* rng, unsigned long int seed); -double nrnisaac_dbl_pick(void* rng); -uint32_t nrnisaac_uint32_pick(void* rng); - - -#endif diff --git a/src/oc/oc_mcran4.cpp b/src/oc/oc_mcran4.cpp new file mode 100644 index 0000000000..128ccd7d04 --- /dev/null +++ b/src/oc/oc_mcran4.cpp @@ -0,0 +1,41 @@ +#include "hocdec.h" +#include "mcran4.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; +} + +void hoc_mcran4() { + uint32_t idx; + double* xidx; + double x; + xidx = hoc_pgetarg(1); + idx = (uint32_t) (*xidx); + x = mcell_ran4a(&idx); + *xidx = idx; + hoc_ret(); + hoc_pushx(x); +} +void hoc_mcran4init() { + double prev = mcell_lowindex(); + if (ifarg(1)) { + uint32_t idx = (uint32_t) chkarg(1, 0., 4294967295.); + 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(prev); +} diff --git a/src/oc/oc_mcran4.hpp b/src/oc/oc_mcran4.hpp new file mode 100644 index 0000000000..d275e3e9bf --- /dev/null +++ b/src/oc/oc_mcran4.hpp @@ -0,0 +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/ocfunc.h b/src/oc/ocfunc.h index fca598c0fa..61a620aa44 100644 --- a/src/oc/ocfunc.h +++ b/src/oc/ocfunc.h @@ -32,7 +32,6 @@ extern void hoc_neuronhome(void), hoc_Execerror(void); extern void hoc_sscanf(void), hoc_save_session(void), hoc_print_session(void); extern void hoc_Chdir(void), hoc_getcwd(void), hoc_Symbol_units(void), hoc_stdout(void); extern void hoc_name_declared(void), hoc_unix_mac_pc(void), hoc_show_winio(void); -extern void hoc_usemcran4(void), hoc_mcran4(void), hoc_mcran4init(void); extern void hoc_nrn_load_dll(void), hoc_nrnversion(void), hoc_object_pushed(void); extern void hoc_mallinfo(void); extern void hoc_Setcolor(void); diff --git a/src/oc/scoprand.cpp b/src/oc/scoprand.cpp index 7f08ae8814..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 +#include "oc_mcran4.hpp" +#include "mcran4.h" #include "scoplib.h" static uint32_t value = 1; @@ -58,8 +59,7 @@ static uint32_t value = 1; *--------------------------------------------------------------------------- */ double scop_random(void) { - extern int use_mcell_ran4_; - if (use_mcell_ran4_) { + if (use_mcran4()) { /*perhaps 4 times slower but much higher quality*/ return mcell_ran4a(&value); } else { From 638c30c7f023fd9445282ab4374feb64325a0b69 Mon Sep 17 00:00:00 2001 From: Luc Grosheintz Date: Sat, 27 Jan 2024 09:24:28 +0100 Subject: [PATCH 03/15] Only use `nrn_prop_datum_alloc`. (#2683) --- src/nrniv/hocmech.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/nrniv/hocmech.cpp b/src/nrniv/hocmech.cpp index 0235647100..e72027c07e 100644 --- a/src/nrniv/hocmech.cpp +++ b/src/nrniv/hocmech.cpp @@ -134,7 +134,7 @@ static void alloc_pnt(Prop* p) { p->ob = nrn_point_prop_->ob; // printf("p->ob comes from nrn_point_prop_ %s\n", hoc_object_name(p->ob)); } else { - p->dparam = (Datum*) hoc_Ecalloc(2, sizeof(Datum)); + nrn_prop_datum_alloc(p->_type, 2, p); if (last_created_pp_ob_) { p->ob = last_created_pp_ob_; // printf("p->ob comes from last_created %s\n", hoc_object_name(p->ob)); From c1e245d04780cb19a925d2432c659435582564aa Mon Sep 17 00:00:00 2001 From: Pramod Kumbhar Date: Sat, 27 Jan 2024 13:41:52 +0100 Subject: [PATCH 04/15] Remove usage of -DDISABLE_OPENACC (GPU build simplification) (#2653) * Remove usage of -DDISABLE_OPENACC (GPU build simplification) - always use Random123 via unified memory support - some of the cpp files were avoided it via `-DDISABLE_OPENACC` but with the addition of RANDOM construct, the Random123 streams are allocated during model setup. So it's better to do this in one way unless there is performance degradation. * update nmodl as well : BlueBrain/nmodl/pull/1133 --- external/nmodl | 2 +- src/coreneuron/CMakeLists.txt | 16 ---------------- src/coreneuron/utils/randoms/nrnran123.h | 5 +---- 3 files changed, 2 insertions(+), 21 deletions(-) diff --git a/external/nmodl b/external/nmodl index 6c89dc639c..bb4dfd2fbc 160000 --- a/external/nmodl +++ b/external/nmodl @@ -1 +1 @@ -Subproject commit 6c89dc639c1ee59e8cfdd69719780e0185f1ec08 +Subproject commit bb4dfd2fbc5209bb1893d3c681157db743ca1dfc diff --git a/src/coreneuron/CMakeLists.txt b/src/coreneuron/CMakeLists.txt index a18bb38d7d..2552c80abd 100644 --- a/src/coreneuron/CMakeLists.txt +++ b/src/coreneuron/CMakeLists.txt @@ -406,22 +406,6 @@ set(CORENEURON_BUILTIN_MODFILES # coreneuron GPU library # ============================================================================= if(CORENRN_ENABLE_GPU) - # ~~~ - # artificial cells and some other cpp files (using Random123) should be compiled - # without OpenACC to avoid use of GPU Random123 streams - # OL210813: this shouldn't be needed anymore, but it may have a small performance benefit - # ~~~ - set(OPENACC_EXCLUDED_FILES - ${CMAKE_CURRENT_BINARY_DIR}/netstim.cpp - ${CMAKE_CURRENT_BINARY_DIR}/netstim_inhpoisson.cpp - ${CMAKE_CURRENT_BINARY_DIR}/pattern.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/io/nrn_setup.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/io/setup_fornetcon.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/io/corenrn_data_return.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/io/global_vars.cpp) - - set_source_files_properties(${OPENACC_EXCLUDED_FILES} PROPERTIES COMPILE_FLAGS - "-DDISABLE_OPENACC") # Only compile the explicit CUDA implementation of the Hines solver in GPU builds. Because of # https://forums.developer.nvidia.com/t/cannot-dynamically-load-a-shared-library-containing-both-openacc-and-cuda-code/210972 # this cannot be included in the same shared library as the rest of the OpenACC code. diff --git a/src/coreneuron/utils/randoms/nrnran123.h b/src/coreneuron/utils/randoms/nrnran123.h index d4108612d0..57b30b2d93 100644 --- a/src/coreneuron/utils/randoms/nrnran123.h +++ b/src/coreneuron/utils/randoms/nrnran123.h @@ -41,10 +41,7 @@ of the full distribution available from #include -// Some files are compiled with DISABLE_OPENACC, and some builds have no GPU -// support at all. In these two cases, request that the random123 state is -// allocated using new/delete instead of CUDA unified memory. -#if defined(CORENEURON_ENABLE_GPU) && !defined(DISABLE_OPENACC) +#if defined(CORENEURON_ENABLE_GPU) #define CORENRN_RAN123_USE_UNIFIED_MEMORY true #else #define CORENRN_RAN123_USE_UNIFIED_MEMORY false From 128c684d4da0781372d56f797f544a9dda1eef1b Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Sat, 27 Jan 2024 17:55:38 +0100 Subject: [PATCH 05/15] Remove usage of IVOS/regexp (#2204) Replace with `std::regex` Co-authored-by: Michael Hines Co-authored-by: Pramod Kumbhar --- cmake/NeuronFileLists.cmake | 2 +- share/lib/hoc/import3d/read_nlcda3.hoc | 2 +- src/ivoc/strfun.cpp | 57 +- src/ivos/InterViews/_defines.h | 1 - src/ivos/InterViews/_undefs.h | 1 - src/ivos/InterViews/regexp.h | 80 -- src/ivos/regexp.cpp | 1219 ------------------------ src/nrncvode/netcvode.cpp | 129 +-- 8 files changed, 83 insertions(+), 1408 deletions(-) delete mode 100755 src/ivos/InterViews/regexp.h delete mode 100644 src/ivos/regexp.cpp diff --git a/cmake/NeuronFileLists.cmake b/cmake/NeuronFileLists.cmake index c1beffb9fd..e42d769355 100644 --- a/cmake/NeuronFileLists.cmake +++ b/cmake/NeuronFileLists.cmake @@ -421,7 +421,7 @@ set(NMODL_FILES_LIST units.cpp version.cpp) -set(IVOS_FILES_LIST observe.cpp regexp.cpp resource.cpp) +set(IVOS_FILES_LIST observe.cpp resource.cpp) set(MPI_DYNAMIC_INCLUDE nrnmpi_dynam.h nrnmpi_dynam_cinc nrnmpi_dynam_wrappers.inc) diff --git a/share/lib/hoc/import3d/read_nlcda3.hoc b/share/lib/hoc/import3d/read_nlcda3.hoc index 142d78cd34..f1a2602ab5 100755 --- a/share/lib/hoc/import3d/read_nlcda3.hoc +++ b/share/lib/hoc/import3d/read_nlcda3.hoc @@ -545,7 +545,7 @@ proc b2soption_split() {local i, n, id, ip localobj p, newsec, tobj } proc remove_trailspace() { // yuck - hoc_sf_.head(line, " *$", line) + hoc_sf_.rtrim(line, line) sprint(line, "%s\n", line) } diff --git a/src/ivoc/strfun.cpp b/src/ivoc/strfun.cpp index f21820b2d8..cb5375eb82 100644 --- a/src/ivoc/strfun.cpp +++ b/src/ivoc/strfun.cpp @@ -1,10 +1,10 @@ #include <../../nrnconf.h> -#include #include #include #include "classreg.h" #include "oc2iv.h" #include +#include // for alias #include #include @@ -39,31 +39,52 @@ static double l_len(void*) { static double l_head(void*) { std::string text(gargstr(1)); - Regexp r(gargstr(2)); - r.Search(text.c_str(), text.size(), 0, text.size()); - int i = r.BeginningOfMatch(); - // text.set_to_left(i); doesnt work - char** head = hoc_pgargstr(3); - if (i > 0) { - hoc_assign_str(head, text.substr(0, i).c_str()); - } else { - hoc_assign_str(head, ""); + { // Clean the text so we keep only the first line + // Imitation of std::multiline in our case + std::regex r("^(.*)(\n|$)"); + std::smatch sm; + std::regex_search(text, sm, r); + text = sm[1]; } + int i = -1; + std::string result{}; + try { + std::regex r(gargstr(2), std::regex::egrep); + if (std::smatch sm; std::regex_search(text, sm, r)) { + i = sm.position(); + result = sm.prefix().str(); + } + } catch (const std::regex_error& e) { + std::cerr << e.what() << std::endl; + } + char** head = hoc_pgargstr(3); + hoc_assign_str(head, result.c_str()); hoc_return_type_code = 1; // integer return double(i); } static double l_tail(void*) { std::string text(gargstr(1)); - Regexp r(gargstr(2)); - r.Search(text.c_str(), text.size(), 0, text.size()); - int i = r.EndOfMatch(); - char** tail = hoc_pgargstr(3); - if (i >= 0) { - hoc_assign_str(tail, text.c_str() + i); - } else { - hoc_assign_str(tail, ""); + { // Clean the text so we keep only the first line + // Imitation of std::multiline in our case + std::regex r("^(.*)(\n|$)"); + std::smatch sm; + std::regex_search(text, sm, r); + text = sm[1]; } + int i = -1; + std::string result{}; + try { + std::regex r(gargstr(2), std::regex::egrep); + if (std::smatch sm; std::regex_search(text, sm, r)) { + i = sm.position() + sm.length(); + result = sm.suffix().str(); + } + } catch (const std::regex_error& e) { + std::cerr << e.what() << std::endl; + } + char** tail = hoc_pgargstr(3); + hoc_assign_str(tail, result.c_str()); hoc_return_type_code = 1; // integer return double(i); } diff --git a/src/ivos/InterViews/_defines.h b/src/ivos/InterViews/_defines.h index 7bb9234e5c..9b2a95fa6f 100755 --- a/src/ivos/InterViews/_defines.h +++ b/src/ivos/InterViews/_defines.h @@ -217,7 +217,6 @@ #define Raster _lib_iv(Raster) #define RasterRep _lib_iv(RasterRep) #define Reducer _lib_iv(Reducer) -#define Regexp _lib_iv(Regexp) #define ReqErr _lib_iv(ReqErr) #define Requirement _lib_iv(Requirement) #define Requisition _lib_iv(Requisition) diff --git a/src/ivos/InterViews/_undefs.h b/src/ivos/InterViews/_undefs.h index a93e162890..6891baa495 100755 --- a/src/ivos/InterViews/_undefs.h +++ b/src/ivos/InterViews/_undefs.h @@ -217,7 +217,6 @@ #undef Raster #undef RasterRep #undef Reducer -#undef Regexp #undef ReqErr #undef Requirement #undef Requisition diff --git a/src/ivos/InterViews/regexp.h b/src/ivos/InterViews/regexp.h deleted file mode 100755 index c1466fd7dc..0000000000 --- a/src/ivos/InterViews/regexp.h +++ /dev/null @@ -1,80 +0,0 @@ -/* - * Copyright (c) 1987, 1988, 1989, 1990, 1991 Stanford University - * Copyright (c) 1991 Silicon Graphics, Inc. - * - * Permission to use, copy, modify, distribute, and sell this software and - * its documentation for any purpose is hereby granted without fee, provided - * that (i) the above copyright notices and this permission notice appear in - * all copies of the software and related documentation, and (ii) the names of - * Stanford and Silicon Graphics may not be used in any advertising or - * publicity relating to the software without the specific, prior written - * permission of Stanford and Silicon Graphics. - * - * THE SOFTWARE IS PROVIDED "AS-IS" AND WITHOUT WARRANTY OF ANY KIND, - * EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY - * WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. - * - * IN NO EVENT SHALL STANFORD OR SILICON GRAPHICS BE LIABLE FOR - * ANY SPECIAL, INCIDENTAL, INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY KIND, - * OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, - * WHETHER OR NOT ADVISED OF THE POSSIBILITY OF DAMAGE, AND ON ANY THEORY OF - * LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE - * OF THIS SOFTWARE. - */ - -/* - * Regexp - regular expression searching - */ - -#ifndef iv_regexp_h -#define iv_regexp_h - -#include - -/* - * These definitions are from Henry Spencers public-domain regular - * expression matching routines. - * - * Definitions etc. for regexp(3) routines. - * - * Caveat: this is V8 regexp(3) [actually, a reimplementation thereof], - * not the System V one. - */ -#define NSUBEXP 10 -struct regexp { - char *startp[NSUBEXP]; - char *endp[NSUBEXP]; - char *textStart; - char regstart; /* Internal use only. */ - char reganch; /* Internal use only. */ - char *regmust; /* Internal use only. */ - int regmlen; /* Internal use only. */ - char program[1]; /* Unwarranted chumminess with compiler. */ -}; - -/* - * The first byte of the regexp internal "program" is actually this magic - * number; the start node begins in the second byte. This used to be the octal - * integer literal 0234 = 156, which would be implicitly converted to -100 when - * narrowing to signed 8 bit char. This conversion was implementation defined - * before C++20. - */ -#define REGEXP_MAGIC static_cast(-100) - -class Regexp { -public: - Regexp(const char*); - Regexp(const char*, int length); - ~Regexp(); - - const char* pattern() const; - int Search(const char* text, int length, int index, int range); - int Match(const char* text, int length, int index); - int BeginningOfMatch(int subexp = 0); - int EndOfMatch(int subexp = 0); -private: - char* pattern_; - regexp* c_pattern; -}; - -#endif diff --git a/src/ivos/regexp.cpp b/src/ivos/regexp.cpp deleted file mode 100644 index ed96c43bcf..0000000000 --- a/src/ivos/regexp.cpp +++ /dev/null @@ -1,1219 +0,0 @@ -#ifdef HAVE_CONFIG_H -#include <../../nrnconf.h> -#endif -/* - * Copyright (c) 1987, 1988, 1989, 1990, 1991 Stanford University - * Copyright (c) 1991 Silicon Graphics, Inc. - * - * Permission to use, copy, modify, distribute, and sell this software and - * its documentation for any purpose is hereby granted without fee, provided - * that (i) the above copyright notices and this permission notice appear in - * all copies of the software and related documentation, and (ii) the names of - * Stanford and Silicon Graphics may not be used in any advertising or - * publicity relating to the software without the specific, prior written - * permission of Stanford and Silicon Graphics. - * - * THE SOFTWARE IS PROVIDED "AS-IS" AND WITHOUT WARRANTY OF ANY KIND, - * EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY - * WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. - * - * IN NO EVENT SHALL STANFORD OR SILICON GRAPHICS BE LIABLE FOR - * ANY SPECIAL, INCIDENTAL, INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY KIND, - * OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, - * WHETHER OR NOT ADVISED OF THE POSSIBILITY OF DAMAGE, AND ON ANY THEORY OF - * LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE - * OF THIS SOFTWARE. - */ - -/* - * Regexp - regular expression searching - */ - -#include -#include -#include - -/* - * This version is based on the Henry Spencers public domain reimplementation - * of the regular expression matching subroutines. They are included as - * static subroutines after the externally accessible routines. - */ - -/* - * Forward declarations for regcomp()'s friends. - */ -static regexp* regcomp(const char* exp); -static char* reg(int paren, int* flagp); -static char* regbranch(int* flagp); -static char* regpiece(int* flagp); -static char* regatom(int* flagp); -static char* regnode(char op); -static char* regnext(char* p); -static void regc(char b); -static void reginsert(char op, char* opnd); -static void regtail(char* p, char* val); -static void regoptail(char* p, char* val); -static void regerror(const char* s); -static int regexec(regexp* prog, char* string); -static int regtry(regexp* prog, char* string); -static int regmatch(char* prog); -static int regrepeat(char* p); - - -inline char * -FindNewline(char* s) { - return strchr(s, '\n'); -} - -inline char * -NextLine(char* s) { - char* newstart; - - if ((newstart = FindNewline(s)) != nil) - newstart++; - return newstart; -} - -Regexp::Regexp (const char* pat) { - int length = strlen(pat); - pattern_ = new char[length+1]; - strncpy(pattern_, pat, length); - pattern_[length] = '\0'; - c_pattern = regcomp(pattern_); - if (!c_pattern) { - delete [] pattern_; - pattern_ = nil; - } -} - -Regexp::Regexp (const char* pat, int length) { - pattern_ = new char[length+1]; - strncpy(pattern_, pat, length); - pattern_[length] = '\0'; - c_pattern = regcomp(pattern_); - if (!c_pattern) { - delete [] pattern_; - pattern_ = nil; - } -} - -Regexp::~Regexp () { - if (pattern_) { - delete [] pattern_; - } - if (c_pattern) { - delete [] c_pattern; - } -} - -const char* Regexp::pattern() const { return pattern_; } - -int Regexp::Search (const char* text, int length, int index, int range) { - bool forwardSearch; - bool frontAnchored; - bool endAnchored; - char* searchStart; - char* searchLimit; - char* endOfLine = nil; - char* lastMatch = nil; - char csave; - - /* - * A small sanity check. Otherwise length is unused in this function. - * This is really what the logic embedded in the old version of this - * routine enforced. - */ - if (index + range > length) { - range = length - index; - if (range < 0) - return -1; - } - - if (c_pattern == nil) { - return -1; - } - - c_pattern->startp[0] = nil; - - if (range < 0) { - forwardSearch = false; - searchLimit = (char *) text + index; - searchStart = (char *) searchLimit + range; /* range is negative */ - } else { - forwardSearch = true; - searchStart = (char *) text + index; - searchLimit = (char *) searchStart + range; - } - - /* Mark end of text string so search will stop */ - char save = *searchLimit; - *searchLimit = '\0'; - - frontAnchored = pattern_[0] == '^'; - endAnchored = pattern_[strlen(pattern_)-1] == '$'; - if (frontAnchored && (searchStart != text && searchStart[-1] != '\n')) { - searchStart = NextLine(searchStart); - } - - while (searchStart && searchStart < searchLimit) { - int result; - - if (endAnchored && (endOfLine = FindNewline(searchStart)) != nil) { - csave = *endOfLine; - *endOfLine = '\0'; - } - - result = regexec(c_pattern, searchStart); - - if (endOfLine) - *endOfLine = csave; - - if (result) { - /* Found a match */ - if (forwardSearch) - break; /* Done */ - else { - lastMatch = c_pattern->startp[0]; - searchStart = c_pattern->endp[0]; - if (frontAnchored) - searchStart = NextLine(searchStart); - continue; - } - } - /* Did not find a match */ - if (frontAnchored || endAnchored) - searchStart = NextLine(searchStart); - else - break; - } - - if (!forwardSearch && lastMatch) { - if (endAnchored && (endOfLine = FindNewline(lastMatch)) != nil) { - csave = *endOfLine; - *endOfLine = '\0'; - } - (void) regexec(c_pattern, lastMatch); /* Refill startp and endp */ - if (endOfLine) - *endOfLine = csave; - } - - *searchLimit = save; - c_pattern->textStart = (char *) text; - - return c_pattern->startp[0] - c_pattern->textStart; -} - -int Regexp::Match (const char* text, int length, int index) { - - if (c_pattern == nil) - return -1; - - c_pattern->startp[0] = nil; - - char save = *(text+length); - *(char*)(text+length) = '\0'; - - c_pattern->textStart = (char *) text; - (void) regexec(c_pattern, (char *) text + index); - - *(char*)(text+length) = save; - - if (c_pattern->startp[0] != nil) - return c_pattern->endp[0] - c_pattern->startp[0]; - else - return -1; -} - -int Regexp::BeginningOfMatch (int subexp) { - if (subexp < 0 || subexp > NSUBEXP || - c_pattern == nil || c_pattern->startp[0] == nil) - return -1; - return c_pattern->startp[subexp] - c_pattern->textStart; -} - -int Regexp::EndOfMatch (int subexp) { - if (subexp < 0 || subexp > NSUBEXP || - c_pattern == nil || c_pattern->startp[0] == nil) - return -1; - return c_pattern->endp[subexp] - c_pattern->textStart; -} - -/* - * regcomp and regexec - * - * Copyright (c) 1986 by University of Toronto. - * Written by Henry Spencer. Not derived from licensed software. - * - * Permission is granted to anyone to use this software for any - * purpose on any computer system, and to redistribute it freely, - * subject to the following restrictions: - * - * 1. The author is not responsible for the consequences of use of - * this software, no matter how awful, even if they arise - * from defects in it. - * - * 2. The origin of this software must not be misrepresented, either - * by explicit claim or by omission. - * - * 3. Altered versions must be plainly marked as such, and must not - * be misrepresented as being the original software. - * - * Beware that some of this code is subtly aware of the way operator - * precedence is structured in regular expressions. Serious changes in - * regular-expression syntax might require a total rethink. - */ - -/* - * The "internal use only" fields in regexp.h are present to pass info from - * compile to execute that permits the execute phase to run lots faster on - * simple cases. They are: - * - * regstart char that must begin a match; '\0' if none obvious - * reganch is the match anchored (at beginning-of-line only)? - * regmust string (pointer into program) that match must include, or nil - * regmlen length of regmust string - * - * Regstart and reganch permit very fast decisions on suitable starting points - * for a match, cutting down the work a lot. Regmust permits fast rejection - * of lines that cannot possibly match. The regmust tests are costly enough - * that regcomp() supplies a regmust only if the r.e. contains something - * potentially expensive (at present, the only such thing detected is * or + - * at the start of the r.e., which can involve a lot of backup). Regmlen is - * supplied because the test in regexec() needs it and regcomp() is computing - * it anyway. - */ - -/* - * Structure for regexp "program". This is essentially a linear encoding - * of a nondeterministic finite-state machine (aka syntax charts or - * "railroad normal form" in parsing technology). Each node is an opcode - * plus a "next" pointer, possibly plus an operand. "Next" pointers of - * all nodes except BRANCH implement concatenation; a "next" pointer with - * a BRANCH on both ends of it is connecting two alternatives. (Here we - * have one of the subtle syntax dependencies: an individual BRANCH (as - * opposed to a collection of them) is never concatenated with anything - * because of operator precedence.) The operand of some types of node is - * a literal string; for others, it is a node leading into a sub-FSM. In - * particular, the operand of a BRANCH node is the first node of the branch. - * (NB this is *not* a tree structure: the tail of the branch connects - * to the thing following the set of BRANCHes.) The opcodes are: - */ - -/* definition number opnd? meaning */ -#define END 0 /* no End of program. */ -#define BOL 1 /* no Match "" at beginning of line. */ -#define EOL 2 /* no Match "" at end of line. */ -#define ANY 3 /* no Match any one character. */ -#define ANYOF 4 /* str Match any character in this string. */ -#define ANYBUT 5 /* str Match any character not in this string. */ -#define BRANCH 6 /* node Match this alternative, or the next... */ -#define BACK 7 /* no Match "", "next" ptr points backward. */ -#define EXACTLY 8 /* str Match this string. */ -#define NOTHING 9 /* no Match empty string. */ -#define STAR 10 /* node Match this (simple) thing 0 or more times. */ -#define PLUS 11 /* node Match this (simple) thing 1 or more times. */ -#define OPEN 20 /* no Mark this point in input as start of #n. */ - /* OPEN+1 is number 1, etc. */ -#define CLOSE 30 /* no Analogous to OPEN. */ - -/* - * Opcode notes: - * - * BRANCH The set of branches constituting a single choice are hooked - * together with their "next" pointers, since precedence prevents - * anything being concatenated to any individual branch. The - * "next" pointer of the last BRANCH in a choice points to the - * thing following the whole choice. This is also where the - * final "next" pointer of each individual branch points; each - * branch starts with the operand node of a BRANCH node. - * - * BACK Normal "next" pointers all implicitly point forward; BACK - * exists to make loop structures possible. - * - * STAR,PLUS '?', and complex '*' and '+', are implemented as circular - * BRANCH structures using BACK. Simple cases (one character - * per match) are implemented with STAR and PLUS for speed - * and to minimize recursive plunges. - * - * OPEN,CLOSE ...are numbered at compile time. - */ - -/* - * A node is one char of opcode followed by two chars of "next" pointer. - * "Next" pointers are stored as two 8-bit pieces, high order first. The - * value is a positive offset from the opcode of the node containing it. - * An operand, if any, simply follows the node. (Note that much of the - * code generation knows about this implicit relationship.) - * - * Using two bytes for the "next" pointer is vast overkill for most things, - * but allows patterns to get big without disasters. - */ -#define OP(p) (*(p)) -#define NEXT(p) (((*((p)+1)&0377)<<8) + (*((p)+2)&0377)) -#define OPERAND(p) ((p) + 3) - -/* - * Utility definitions. - */ - -/** - * This replaces a macro of the same name with some bit manipulation magic in - * it. The does not seem well-suited now, but it's not clear that it was before - * either. - */ -inline int UCHARAT(const char* p) { - return *p; -} - -#define FAIL(m) { regerror(m); return(nil); } -#define ISMULT(c) ((c) == '*' || (c) == '+' || (c) == '?') -#define META "^$.[()|?+*\\" - -/* - * Flags to be passed up and down. - */ -#define HASWIDTH 01 /* Known never to match null string. */ -#define SIMPLE 02 /* Simple enough to be STAR/PLUS operand. */ -#define SPSTART 04 /* Starts with * or +. */ -#define WORST 0 /* Worst case. */ - -/* - * Global work variables for regcomp(). - */ -static const char *regparse; /* Input-scan pointer. */ -static int regnpar; /* () count. */ -static char regdummy; -static char *regcode; /* Code-emit pointer; ®dummy = don't. */ -static long regsize; /* Code size. */ - -/* - - regcomp - compile a regular expression into internal code - * - * We can't allocate space until we know how big the compiled form will be, - * but we can't compile it (and thus know how big it is) until we've got a - * place to put the code. So we cheat: we compile it twice, once with code - * generation turned off and size counting turned on, and once "for real". - * This also means that we don't allocate space until we are sure that the - * thing really will compile successfully, and we never have to move the - * code and thus invalidate pointers into it. (Note that it has to be in - * one piece because free() must be able to free it all.) - * - * Beware that the optimization-preparation code in here knows about some - * of the structure of the compiled regexp. - */ -static regexp * -regcomp(const char* exp) { - regexp *r; - char *scan; - char *longest; - int len; - int flags; - - if (exp == nil) - FAIL("nil argument"); - - /* First pass: determine size, legality. */ - regparse = exp; - regnpar = 1; - regsize = 0L; - regcode = ®dummy; - regc(REGEXP_MAGIC); - if (reg(0, &flags) == nil) - return(nil); - - /* Small enough for pointer-storage convention? */ - if (regsize >= 32767L) /* Probably could be 65535L. */ - FAIL("regexp too big"); - - /* Allocate space. */ - r = (regexp *) new char[sizeof(regexp) + (unsigned)regsize]; - - /* Second pass: emit code. */ - regparse = exp; - regnpar = 1; - regcode = r->program; - regc(REGEXP_MAGIC); - if (reg(0, &flags) == nil) { - delete [] r; - return(nil); - } - - /* Dig out information for optimizations. */ - r->regstart = '\0'; /* Worst-case defaults. */ - r->reganch = 0; - r->regmust = nil; - r->regmlen = 0; - scan = r->program+1; /* First BRANCH. */ - if (OP(regnext(scan)) == END) { /* Only one top-level choice. */ - scan = OPERAND(scan); - - /* Starting-point info. */ - if (OP(scan) == EXACTLY) - r->regstart = *OPERAND(scan); - else if (OP(scan) == BOL) - r->reganch++; - - /* - * If there's something expensive in the r.e., find the - * longest literal string that must appear and make it the - * regmust. Resolve ties in favor of later strings, since - * the regstart check works with the beginning of the r.e. - * and avoiding duplication strengthens checking. Not a - * strong reason, but sufficient in the absence of others. - */ - if (flags&SPSTART) { - longest = nil; - len = 0; - for (; scan != nil; scan = regnext(scan)) - if (OP(scan) == EXACTLY && strlen(OPERAND(scan)) >= len) { - longest = OPERAND(scan); - len = strlen(OPERAND(scan)); - } - r->regmust = longest; - r->regmlen = len; - } - } - - return(r); -} - -/* - - reg - regular expression, i.e. main body or parenthesized thing - * - * Caller must absorb opening parenthesis. - * - * Combining parenthesis handling with the base level of regular expression - * is a trifle forced, but the need to tie the tails of the branches to what - * follows makes it hard to avoid. - */ -static char * -reg(int paren, int* flagp) { - char *ret; - char *br; - char *ender; - int parno; - int flags; - - *flagp = HASWIDTH; /* Tentatively. */ - - /* Make an OPEN node, if parenthesized. */ - if (paren) { - if (regnpar >= NSUBEXP) - FAIL("too many ()"); - parno = regnpar; - regnpar++; - ret = regnode(OPEN+parno); - } else - ret = nil; - - /* Pick up the branches, linking them together. */ - br = regbranch(&flags); - if (br == nil) - return(nil); - if (ret != nil) - regtail(ret, br); /* OPEN -> first. */ - else - ret = br; - if (!(flags&HASWIDTH)) - *flagp &= ~HASWIDTH; - *flagp |= flags&SPSTART; - while (*regparse == '|') { - regparse++; - br = regbranch(&flags); - if (br == nil) - return(nil); - regtail(ret, br); /* BRANCH -> BRANCH. */ - if (!(flags&HASWIDTH)) - *flagp &= ~HASWIDTH; - *flagp |= flags&SPSTART; - } - - /* Make a closing node, and hook it on the end. */ - ender = regnode((paren) ? CLOSE+parno : END); - regtail(ret, ender); - - /* Hook the tails of the branches to the closing node. */ - for (br = ret; br != nil; br = regnext(br)) - regoptail(br, ender); - - /* Check for proper termination. */ - if (paren && *regparse++ != ')') { - FAIL("unmatched ()"); - } else if (!paren && *regparse != '\0') { - if (*regparse == ')') { - FAIL("unmatched ()"); - } else - FAIL("junk on end"); /* "Can't happen". */ - /* NOTREACHED */ - } - - return(ret); -} - -/* - - regbranch - one alternative of an | operator - * - * Implements the concatenation operator. - */ -static char * -regbranch(int* flagp) { - char *ret; - char *chain; - char *latest; - int flags; - - *flagp = WORST; /* Tentatively. */ - - ret = regnode(BRANCH); - chain = nil; - while (*regparse != '\0' && *regparse != '|') { - if (*regparse == '\\' && regparse[1] == ')') { - regparse++; - break; - } - latest = regpiece(&flags); - if (latest == nil) - return(nil); - *flagp |= flags&HASWIDTH; - if (chain == nil) /* First piece. */ - *flagp |= flags&SPSTART; - else - regtail(chain, latest); - chain = latest; - } - if (chain == nil) /* Loop ran zero times. */ - (void) regnode(NOTHING); - - return(ret); -} - -/* - - regpiece - something followed by possible [*+?] - * - * Note that the branching code sequences used for ? and the general cases - * of * and + are somewhat optimized: they use the same NOTHING node as - * both the endmarker for their branch list and the body of the last branch. - * It might seem that this node could be dispensed with entirely, but the - * endmarker role is not redundant. - */ -static char * -regpiece(int* flagp) { - char *ret; - char op; - char *next; - int flags; - - ret = regatom(&flags); - if (ret == nil) - return(nil); - - op = *regparse; - if (!ISMULT(op)) { - *flagp = flags; - return(ret); - } - - if (!(flags&HASWIDTH) && op != '?') - FAIL("*+ operand could be empty"); - *flagp = (op != '+') ? (WORST|SPSTART) : (WORST|HASWIDTH); - - if (op == '*' && (flags&SIMPLE)) - reginsert(STAR, ret); - else if (op == '*') { - /* Emit x* as (x&|), where & means "self". */ - reginsert(BRANCH, ret); /* Either x */ - regoptail(ret, regnode(BACK)); /* and loop */ - regoptail(ret, ret); /* back */ - regtail(ret, regnode(BRANCH)); /* or */ - regtail(ret, regnode(NOTHING)); /* null. */ - } else if (op == '+' && (flags&SIMPLE)) - reginsert(PLUS, ret); - else if (op == '+') { - /* Emit x+ as x(&|), where & means "self". */ - next = regnode(BRANCH); /* Either */ - regtail(ret, next); - regtail(regnode(BACK), ret); /* loop back */ - regtail(next, regnode(BRANCH)); /* or */ - regtail(ret, regnode(NOTHING)); /* null. */ - } else if (op == '?') { - /* Emit x? as (x|) */ - reginsert(BRANCH, ret); /* Either x */ - regtail(ret, regnode(BRANCH)); /* or */ - next = regnode(NOTHING); /* null. */ - regtail(ret, next); - regoptail(ret, next); - } - regparse++; - if (ISMULT(*regparse)) - FAIL("nested *?+"); - - return(ret); -} - -/* - - regatom - the lowest level - * - * Optimization: gobbles an entire sequence of ordinary characters so that - * it can turn them into a single node, which is smaller to store and - * faster to run. Backslashed characters are exceptions, each becoming a - * separate node; the code is simpler that way and it's not worth fixing. - */ -static char * -regatom(int* flagp) { - char *ret; - int flags; - - *flagp = WORST; /* Tentatively. */ - - switch (*regparse++) { - case '^': - ret = regnode(BOL); - break; - case '$': - ret = regnode(EOL); - break; - case '.': - ret = regnode(ANY); - *flagp |= HASWIDTH|SIMPLE; - break; - case '[': { - int classbeg; - int classend; - - if (*regparse == '^') { /* Complement of range. */ - ret = regnode(ANYBUT); - regparse++; - } else - ret = regnode(ANYOF); - if (*regparse == ']' || *regparse == '-') - regc(*regparse++); - while (*regparse != '\0' && *regparse != ']') { - if (*regparse == '-') { - regparse++; - if (*regparse == ']' || *regparse == '\0') - regc('-'); - else { - classbeg = UCHARAT(regparse-2)+1; - classend = UCHARAT(regparse); - if (classbeg > classend+1) - FAIL("invalid [] range"); - for (; classbeg <= classend; classbeg++) - regc(classbeg); - regparse++; - } - } else - regc(*regparse++); - } - regc('\0'); - if (*regparse != ']') - FAIL("unmatched []"); - regparse++; - *flagp |= HASWIDTH|SIMPLE; - } - break; - case '\0': - case '|': - FAIL("internal urp"); /* Supposed to be caught earlier. */ - break; - case '?': - case '+': - case '*': - FAIL("?+* follows nothing"); - break; - case '\\': - if (*regparse == '\0') - FAIL("trailing \\"); - if (*regparse == '(') { - regparse++; - ret = reg(1, &flags); - if (ret == nil) - return(nil); - *flagp |= flags&(HASWIDTH|SPSTART); - } else { - ret = regnode(EXACTLY); - regc(*regparse++); - regc('\0'); - *flagp |= HASWIDTH|SIMPLE; - } - break; - default: { - int len; - char ender; - - regparse--; - len = strcspn(regparse, META); - if (len <= 0) - FAIL("internal disaster"); - ender = *(regparse+len); - if (len > 1 && ISMULT(ender)) - len--; /* Back off clear of ?+* operand. */ - *flagp |= HASWIDTH; - if (len == 1) - *flagp |= SIMPLE; - ret = regnode(EXACTLY); - while (len > 0) { - regc(*regparse++); - len--; - } - regc('\0'); - } - break; - } - - return(ret); -} - -/* - - regnode - emit a node - */ -static char * /* Location. */ -regnode(char op) { - char *ret; - char *ptr; - - ret = regcode; - if (ret == ®dummy) { - regsize += 3; - return(ret); - } - - ptr = ret; - *ptr++ = op; - *ptr++ = '\0'; /* Null "next" pointer. */ - *ptr++ = '\0'; - regcode = ptr; - - return(ret); -} - -/* - - regc - emit (if appropriate) a byte of code - */ -static void -regc(char b) { - if (regcode != ®dummy) - *regcode++ = b; - else - regsize++; -} - -/* - - reginsert - insert an operator in front of already-emitted operand - * - * Means relocating the operand. - */ -static void -reginsert(char op, char* opnd) { - char *src; - char *dst; - char *place; - - if (regcode == ®dummy) { - regsize += 3; - return; - } - - src = regcode; - regcode += 3; - dst = regcode; - while (src > opnd) - *--dst = *--src; - - place = opnd; /* Op node, where operand used to be. */ - *place++ = op; - *place++ = '\0'; - *place++ = '\0'; -} - -/* - - regtail - set the next-pointer at the end of a node chain - */ -static void -regtail(char* p, char* val) { - char *scan; - char *temp; - int offset; - - if (p == ®dummy) - return; - - /* Find last node. */ - scan = p; - for (;;) { - temp = regnext(scan); - if (temp == nil) - break; - scan = temp; - } - - if (OP(scan) == BACK) - offset = scan - val; - else - offset = val - scan; - *(scan+1) = (offset>>8)&0377; - *(scan+2) = offset&0377; -} - -/* - - regoptail - regtail on operand of first argument; nop if operandless - */ -static void -regoptail(char* p, char* val) { - /* "Operandless" and "op != BRANCH" are synonymous in practice. */ - if (p == nil || p == ®dummy || OP(p) != BRANCH) - return; - regtail(OPERAND(p), val); -} - -/* - * regexec and friends - */ - -/* - * Global work variables for regexec(). - */ -static char *reginput; /* String-input pointer. */ -static char *regbol; /* Beginning of input, for ^ check. */ -static char **regstartp; /* Pointer to startp array. */ -static char **regendp; /* Ditto for endp. */ - -/* - - regexec - match a regexp against a string - */ -static int -regexec(regexp* prog, char* string) { - char *s; - - /* Be paranoid... */ - if (prog == nil || string == nil) { - regerror("nil parameter"); - return(0); - } - - /* Check validity of program. */ - if (UCHARAT(prog->program) != REGEXP_MAGIC) { - regerror("corrupted program"); - return(0); - } - - /* If there is a "must appear" string, look for it. */ - if (prog->regmust != nil) { - s = string; - while ((s = strchr(s, prog->regmust[0])) != nil) { - if (strncmp(s, prog->regmust, prog->regmlen) == 0) - break; /* Found it. */ - s++; - } - if (s == nil) /* Not present. */ - return(0); - } - - /* Mark beginning of line for ^ . */ - regbol = string; - - /* Simplest case: anchored match need be tried only once. */ - if (prog->reganch) - return(regtry(prog, string)); - - /* Messy cases: unanchored match. */ - s = string; - if (prog->regstart != '\0') - /* We know what char it must start with. */ - while ((s = strchr(s, prog->regstart)) != nil) { - if (regtry(prog, s)) - return(1); - s++; - } - else - /* We don't -- general case. */ - do { - if (regtry(prog, s)) - return(1); - } while (*s++ != '\0'); - - /* Failure. */ - return(0); -} - -/* - - regtry - try match at specific point - */ -static int /* 0 failure, 1 success */ -regtry(regexp* prog, char* string) { - int i; - char **sp; - char **ep; - - reginput = string; - regstartp = prog->startp; - regendp = prog->endp; - - sp = prog->startp; - ep = prog->endp; - for (i = NSUBEXP; i > 0; i--) { - *sp++ = nil; - *ep++ = nil; - } - if (regmatch(prog->program + 1)) { - prog->startp[0] = string; - prog->endp[0] = reginput; - return(1); - } else - return(0); -} - -/* - - regmatch - main matching routine - * - * Conceptually the strategy is simple: check to see whether the current - * node matches, call self recursively to see whether the rest matches, - * and then act accordingly. In practice we make some effort to avoid - * recursion, in particular by going through "ordinary" nodes (that don't - * need to know whether the rest of the match failed) by a loop instead of - * by recursion. - */ -static int /* 0 failure, 1 success */ -regmatch(char* prog) { - char *scan; /* Current node. */ - char *next; /* Next node. */ - - scan = prog; - while (scan != nil) { - next = regnext(scan); - - switch (OP(scan)) { - case BOL: - if (reginput != regbol) - return(0); - break; - case EOL: - if (*reginput != '\0') - return(0); - break; - case ANY: - if (*reginput == '\0') - return(0); - reginput++; - break; - case EXACTLY: { - int len; - char *opnd; - - opnd = OPERAND(scan); - /* Inline the first character, for speed. */ - if (*opnd != *reginput) - return(0); - len = strlen(opnd); - if (len > 1 && strncmp(opnd, reginput, len) != 0) - return(0); - reginput += len; - } - break; - case ANYOF: - if (*reginput == '\0') - return(0); - if (strchr(OPERAND(scan), *reginput) == nil) - return(0); - reginput++; - break; - case ANYBUT: - if (*reginput == '\0') - return(0); - if (strchr(OPERAND(scan), *reginput) != nil) - return(0); - reginput++; - break; - case NOTHING: - break; - case BACK: - break; - case OPEN+1: - case OPEN+2: - case OPEN+3: - case OPEN+4: - case OPEN+5: - case OPEN+6: - case OPEN+7: - case OPEN+8: - case OPEN+9: { - int no; - char *save; - - no = OP(scan) - OPEN; - save = reginput; - - if (regmatch(next)) { - /* - * Don't set startp if some later - * invocation of the same parentheses - * already has. - */ - if (regstartp[no] == nil) - regstartp[no] = save; - return(1); - } else - return(0); - } - break; - case CLOSE+1: - case CLOSE+2: - case CLOSE+3: - case CLOSE+4: - case CLOSE+5: - case CLOSE+6: - case CLOSE+7: - case CLOSE+8: - case CLOSE+9: { - int no; - char *save; - - no = OP(scan) - CLOSE; - save = reginput; - - if (regmatch(next)) { - /* - * Don't set endp if some later - * invocation of the same parentheses - * already has. - */ - if (regendp[no] == nil) - regendp[no] = save; - return(1); - } else - return(0); - } - break; - case BRANCH: { - char *save; - - if (OP(next) != BRANCH) /* No choice. */ - next = OPERAND(scan); /* Avoid recursion. */ - else { - do { - save = reginput; - if (regmatch(OPERAND(scan))) - return(1); - reginput = save; - scan = regnext(scan); - } while (scan != nil && OP(scan) == BRANCH); - return(0); - /* NOTREACHED */ - } - } - break; - case STAR: - case PLUS: { - char nextch; - int no; - char *save; - int min; - - /* - * Lookahead to avoid useless match attempts - * when we know what character comes next. - */ - nextch = '\0'; - if (OP(next) == EXACTLY) - nextch = *OPERAND(next); - min = (OP(scan) == STAR) ? 0 : 1; - save = reginput; - no = regrepeat(OPERAND(scan)); - while (no >= min) { - /* If it could work, try it. */ - if (nextch == '\0' || *reginput == nextch) - if (regmatch(next)) - return(1); - /* Couldn't or didn't -- back up. */ - no--; - reginput = save + no; - } - return(0); - } - break; - case END: - return(1); /* Success! */ - default: - regerror("memory corruption"); - return(0); - } - - scan = next; - } - - /* - * We get here only if there's trouble -- normally "case END" is - * the terminating point. - */ - regerror("corrupted pointers"); - return(0); -} - -/* - - regrepeat - repeatedly match something simple, report how many - */ -static int -regrepeat(char* p) { - int count = 0; - char *scan; - char *opnd; - - scan = reginput; - opnd = OPERAND(p); - switch (OP(p)) { - case ANY: - count = strlen(scan); - scan += count; - break; - case EXACTLY: - while (*opnd == *scan) { - count++; - scan++; - } - break; - case ANYOF: - while (*scan != '\0' && strchr(opnd, *scan) != nil) { - count++; - scan++; - } - break; - case ANYBUT: - while (*scan != '\0' && strchr(opnd, *scan) == nil) { - count++; - scan++; - } - break; - default: /* Oh dear. Called inappropriately. */ - regerror("internal foulup"); - count = 0; /* Best compromise. */ - break; - } - reginput = scan; - - return(count); -} - -/* - - regnext - dig the "next" pointer out of a node - */ -static char * -regnext(char* p) { - int offset; - - if (p == ®dummy) - return(nil); - - offset = NEXT(p); - if (offset == 0) - return(nil); - - if (OP(p) == BACK) - return(p-offset); - else - return(p+offset); -} - -static void -regerror(const char* s) { - std::cerr << "regexp: " << s << "\n"; -} - diff --git a/src/nrncvode/netcvode.cpp b/src/nrncvode/netcvode.cpp index d2dd5ffcb5..f86afcfdce 100644 --- a/src/nrncvode/netcvode.cpp +++ b/src/nrncvode/netcvode.cpp @@ -7,7 +7,7 @@ #include #include #include -#include +#include #include "classreg.h" #include "nrnoc2iv.h" #include "parse.hpp" @@ -930,98 +930,74 @@ Object** NetCvode::netconlist() { Object** po = newoclist(4, o); Object *opre = nullptr, *opost = nullptr, *otar = nullptr; - Regexp *spre = nullptr, *spost = nullptr, *star = nullptr; - char* s; - int n; + std::regex spre, spost, star; if (hoc_is_object_arg(1)) { opre = *hoc_objgetarg(1); } else { - s = gargstr(1); - if (s[0] == '\0') { - spre = new Regexp(".*"); + std::string s(gargstr(1)); + if (s.empty()) { + spre = std::regex(".*"); } else { - spre = new Regexp(escape_bracket(s)); - } - if (!spre->pattern()) { - delete std::exchange(spre, nullptr); - hoc_execerror(gargstr(1), "not a valid regular expression"); + try { + spre = std::regex(escape_bracket(s.data())); + } catch (std::regex_error&) { + hoc_execerror(gargstr(1), "not a valid regular expression"); + } } } if (hoc_is_object_arg(2)) { opost = *hoc_objgetarg(2); } else { - s = gargstr(2); - if (s[0] == '\0') { - spost = new Regexp(".*"); + std::string s(gargstr(2)); + if (s.empty()) { + spost = std::regex(".*"); } else { - spost = new Regexp(escape_bracket(s)); - } - if (!spost->pattern()) { - delete std::exchange(spost, nullptr); - delete std::exchange(spre, nullptr); - hoc_execerror(gargstr(2), "not a valid regular expression"); + try { + spost = std::regex(escape_bracket(s.data())); + } catch (std::regex_error&) { + hoc_execerror(gargstr(2), "not a valid regular expression"); + } } } if (hoc_is_object_arg(3)) { otar = *hoc_objgetarg(3); } else { - s = gargstr(3); - if (s[0] == '\0') { - star = new Regexp(".*"); + std::string s(gargstr(3)); + if (s.empty()) { + star = std::regex(".*"); } else { - star = new Regexp(escape_bracket(s)); - } - if (!star->pattern()) { - delete std::exchange(star, nullptr); - delete std::exchange(spre, nullptr); - delete std::exchange(spost, nullptr); - hoc_execerror(gargstr(3), "not a valid regular expression"); + try { + star = std::regex(escape_bracket(s.data())); + } catch (std::regex_error&) { + hoc_execerror(gargstr(3), "not a valid regular expression"); + } } } - bool b; hoc_Item* q; if (psl_) { ITERATE(q, psl_) { PreSyn* ps = (PreSyn*) VOIDITM(q); - b = false; + bool b = false; if (ps->ssrc_) { Object* precell = nrn_sec2cell(ps->ssrc_); if (opre) { - if (precell == opre) { - b = true; - } else { - b = false; - } + b = precell == opre; } else { - s = hoc_object_name(precell); - n = strlen(s); - if (spre->Match(s, n, 0) > 0) { - b = true; - } else { - b = false; - } + std::string s(hoc_object_name(precell)); + b = std::regex_search(s, spre); } } else if (ps->osrc_) { Object* presyn = ps->osrc_; if (opre) { - if (presyn == opre) { - b = true; - } else { - b = false; - } + b = presyn == opre; } else { - s = hoc_object_name(presyn); - n = strlen(s); - if (spre->Match(s, n, 0) > 0) { - b = true; - } else { - b = false; - } + std::string s(hoc_object_name(presyn)); + b = std::regex_search(s, spre); } } - if (b == true) { + if (b) { for (const auto& d: ps->dil_) { Object* postcell = nullptr; Object* target = nullptr; @@ -1033,37 +1009,19 @@ Object** NetCvode::netconlist() { } } if (opost) { - if (postcell == opost) { - b = true; - } else { - b = false; - } + b = postcell == opost; } else { - s = hoc_object_name(postcell); - n = strlen(s); - if (spost->Match(s, n, 0) > 0) { - b = true; - } else { - b = false; - } + std::string s(hoc_object_name(postcell)); + b = std::regex_search(s, spost); } - if (b == true) { + if (b) { if (otar) { - if (target == otar) { - b = true; - } else { - b = false; - } + b = target == otar; } else { - s = hoc_object_name(target); - n = strlen(s); - if (star->Match(s, n, 0) > 0) { - b = true; - } else { - b = false; - } + std::string s(hoc_object_name(target)); + b = std::regex_search(s, star); } - if (b == true) { + if (b) { o->append(d->obj_); } } @@ -1071,9 +1029,6 @@ Object** NetCvode::netconlist() { } } } - delete std::exchange(spre, nullptr); - delete std::exchange(spost, nullptr); - delete std::exchange(star, nullptr); return po; } From ed0b46126332e5f817e90368aff5ac14114a57dc Mon Sep 17 00:00:00 2001 From: JCGoran Date: Mon, 29 Jan 2024 19:12:45 +0100 Subject: [PATCH 06/15] Refactor `nrn_mlh_gsort` with `std::sort` (#2692) Co-authored-by: Nicolas Cornu --- src/ivoc/ivocvect.cpp | 196 +------------------------------------- test/CMakeLists.txt | 1 + test/unit_tests/iovec.cpp | 45 +++++++++ 3 files changed, 50 insertions(+), 192 deletions(-) create mode 100644 test/unit_tests/iovec.cpp diff --git a/src/ivoc/ivocvect.cpp b/src/ivoc/ivocvect.cpp index 7b9087d328..866fab6836 100644 --- a/src/ivoc/ivocvect.cpp +++ b/src/ivoc/ivocvect.cpp @@ -1,6 +1,7 @@ #include <../../nrnconf.h> //#include +#include #include #include #include @@ -3875,198 +3876,9 @@ void Vector_reg() { #endif } -// hacked version of gsort from ../gnu/d_vec.cpp -// the transformation is that everything that used to be a double* becomes -// an int* and cmp(*arg1, *arg2) becomes cmp(vec[*arg1], vec[*arg2]) -// I am not sure what to do about the BYTES_PER_WORD - -// An adaptation of Schmidt's new quicksort - -static inline void SWAP(int* A, int* B) { - int tmp = *A; - *A = *B; - *B = tmp; -} - -/* This should be replaced by a standard ANSI macro. */ -#define BYTES_PER_WORD 8 -#define BYTES_PER_LONG 4 - -/* The next 4 #defines implement a very fast in-line stack abstraction. */ - -#define STACK_SIZE (BYTES_PER_WORD * BYTES_PER_LONG) -#define PUSH(LOW, HIGH) \ - do { \ - top->lo = LOW; \ - top++->hi = HIGH; \ - } while (0) -#define POP(LOW, HIGH) \ - do { \ - LOW = (--top)->lo; \ - HIGH = top->hi; \ - } while (0) -#define STACK_NOT_EMPTY (stack < top) - -/* Discontinue quicksort algorithm when partition gets below this size. - This particular magic number was chosen to work best on a Sun 4/260. */ -#define MAX_THRESH 4 - - -/* Order size using quicksort. This implementation incorporates - four optimizations discussed in Sedgewick: - - 1. Non-recursive, using an explicit stack of pointer that - store the next array partition to sort. To save time, this - maximum amount of space required to store an array of - MAX_INT is allocated on the stack. Assuming a 32-bit integer, - this needs only 32 * sizeof (stack_node) == 136 bits. Pretty - cheap, actually. - - 2. Chose the pivot element using a median-of-three decision tree. - This reduces the probability of selecting a bad pivot value and - eliminates certain extraneous comparisons. - - 3. Only quicksorts TOTAL_ELEMS / MAX_THRESH partitions, leaving - insertion sort to order the MAX_THRESH items within each partition. - This is a big win, since insertion sort is faster for small, mostly - sorted array segements. - - 4. The larger of the two sub-partitions is always pushed onto the - stack first, with the algorithm then concentrating on the - smaller partition. This *guarantees* no more than log (n) - stack size is needed! */ - int nrn_mlh_gsort(double* vec, int* base_ptr, int total_elems, int (*cmp)(double, double)) { - /* Stack node declarations used to store unfulfilled partition obligations. */ - struct stack_node { - int* lo; - int* hi; - }; - int pivot_buffer; - int max_thresh = MAX_THRESH; - - if (total_elems > MAX_THRESH) { - int* lo = base_ptr; - int* hi = lo + (total_elems - 1); - int* left_ptr; - int* right_ptr; - stack_node stack[STACK_SIZE]; /* Largest size needed for 32-bit int!!! */ - stack_node* top = stack + 1; - - while (STACK_NOT_EMPTY) { - { - int* pivot = &pivot_buffer; - { - /* Select median value from among LO, MID, and HI. Rearrange - LO and HI so the three values are sorted. This lowers the - probability of picking a pathological pivot value and - skips a comparison for both the LEFT_PTR and RIGHT_PTR. */ - - int* mid = lo + ((hi - lo) >> 1); - - if (cmp(vec[*mid], vec[*lo]) < 0) - SWAP(mid, lo); - if (cmp(vec[*hi], vec[*mid]) < 0) { - SWAP(mid, hi); - if (cmp(vec[*mid], vec[*lo]) < 0) - SWAP(mid, lo); - } - *pivot = *mid; - pivot = &pivot_buffer; - } - left_ptr = lo + 1; - right_ptr = hi - 1; - - /* Here's the famous ``collapse the walls'' section of quicksort. - Gotta like those tight inner loops! They are the main reason - that this algorithm runs much faster than others. */ - do { - while (cmp(vec[*left_ptr], vec[*pivot]) < 0) - left_ptr += 1; - - while (cmp(vec[*pivot], vec[*right_ptr]) < 0) - right_ptr -= 1; - - if (left_ptr < right_ptr) { - SWAP(left_ptr, right_ptr); - left_ptr += 1; - right_ptr -= 1; - } else if (left_ptr == right_ptr) { - left_ptr += 1; - right_ptr -= 1; - break; - } - } while (left_ptr <= right_ptr); - } - - /* Set up pointers for next iteration. First determine whether - left and right partitions are below the threshold size. If so, - ignore one or both. Otherwise, push the larger partition's - bounds on the stack and continue sorting the smaller one. */ - - if ((right_ptr - lo) <= max_thresh) { - if ((hi - left_ptr) <= max_thresh) /* Ignore both small partitions. */ - POP(lo, hi); - else /* Ignore small left partition. */ - lo = left_ptr; - } else if ((hi - left_ptr) <= max_thresh) /* Ignore small right partition. */ - hi = right_ptr; - else if ((right_ptr - lo) > (hi - left_ptr)) /* Push larger left partition indices. */ - { - PUSH(lo, right_ptr); - lo = left_ptr; - } else /* Push larger right partition indices. */ - { - PUSH(left_ptr, hi); - hi = right_ptr; - } - } - } - - /* Once the BASE_PTR array is partially sorted by quicksort the rest - is completely sorted using insertion sort, since this is efficient - for partitions below MAX_THRESH size. BASE_PTR points to the beginning - of the array to sort, and END_PTR points at the very last element in - the array (*not* one beyond it!). */ - - - { - int* end_ptr = base_ptr + 1 * (total_elems - 1); - int* run_ptr; - int* tmp_ptr = base_ptr; - int* thresh = (end_ptr < (base_ptr + max_thresh)) ? end_ptr : (base_ptr + max_thresh); - - /* Find smallest element in first threshold and place it at the - array's beginning. This is the smallest array element, - and the operation speeds up insertion sort's inner loop. */ - - for (run_ptr = tmp_ptr + 1; run_ptr <= thresh; run_ptr += 1) - if (cmp(vec[*run_ptr], vec[*tmp_ptr]) < 0) - tmp_ptr = run_ptr; - - if (tmp_ptr != base_ptr) - SWAP(tmp_ptr, base_ptr); - - /* Insertion sort, running from left-hand-side up to `right-hand-side.' - Pretty much straight out of the original GNU qsort routine. */ - - for (run_ptr = base_ptr + 1; (tmp_ptr = run_ptr += 1) <= end_ptr;) { - while (cmp(vec[*run_ptr], vec[*(tmp_ptr -= 1)]) < 0) - ; - - if ((tmp_ptr += 1) != run_ptr) { - int* trav; - - for (trav = run_ptr + 1; --trav >= run_ptr;) { - int c = *trav; - int *hi, *lo; - - for (hi = lo = trav; (lo -= 1) >= tmp_ptr; hi = lo) - *hi = *lo; - *hi = c; - } - } - } - } + std::sort(base_ptr, base_ptr + total_elems, [&](int a, int b) { + return cmp(vec[a], vec[b]) < 0; + }); return 1; } diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index d5e9762cd7..e5731f55fb 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -17,6 +17,7 @@ add_executable( testneuron common/catch2_main.cpp unit_tests/basic.cpp + unit_tests/iovec.cpp unit_tests/container/container.cpp unit_tests/container/generic_data_handle.cpp unit_tests/container/mechanism.cpp diff --git a/test/unit_tests/iovec.cpp b/test/unit_tests/iovec.cpp new file mode 100644 index 0000000000..946ad7c7fa --- /dev/null +++ b/test/unit_tests/iovec.cpp @@ -0,0 +1,45 @@ +#include +#include + +#include "oc_ansi.h" + +#include + +// This function is the one that is used in all nrn-modeldb-ci +// Keep as is +int cmpdfn(double a, double b) { + return ((a) <= (b)) ? (((a) == (b)) ? 0 : -1) : 1; +} + +TEST_CASE("Test nrn_mlh_gsort output", "[nrn_gsort]") { + std::vector input{1.2, -2.5, 5.1}; + + { + std::vector indices(input.size()); + // all values from 0 to size - 1 + std::iota(indices.begin(), indices.end(), 0); + + // for comparison + auto sorted_input = input; + std::sort(sorted_input.begin(), sorted_input.end()); + + SECTION("Test sorting") { + nrn_mlh_gsort(input.data(), indices.data(), input.size(), cmpdfn); + for (auto i = 0; i < input.size(); ++i) { + REQUIRE(sorted_input[i] == input[indices[i]]); + } + } + } + + { + std::vector indices{2, 1, 1}; + std::vector expected_result{1, 1, 2}; // as -2,5 < 5.1 + + SECTION("Test sorting with repeated indices") { + nrn_mlh_gsort(input.data(), indices.data(), input.size(), cmpdfn); + for (auto i = 0; i < input.size(); ++i) { + REQUIRE(indices[i] == expected_result[i]); + } + } + } +} From a71036671ab736c810f2428eec31698fa8816a6e Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 30 Jan 2024 10:30:52 +0100 Subject: [PATCH 07/15] Replace hoc_regexp_* by std::regex (#2694) --- cmake/NeuronFileLists.cmake | 1 - src/nrnoc/cabcode.cpp | 46 +++- src/oc/oc_ansi.h | 2 - src/oc/regexp.cpp | 466 ------------------------------------ 4 files changed, 39 insertions(+), 476 deletions(-) delete mode 100644 src/oc/regexp.cpp diff --git a/cmake/NeuronFileLists.cmake b/cmake/NeuronFileLists.cmake index e42d769355..fb6ca521e7 100644 --- a/cmake/NeuronFileLists.cmake +++ b/cmake/NeuronFileLists.cmake @@ -103,7 +103,6 @@ set(OC_FILE_LIST ocerf.cpp plot.cpp plt.cpp - regexp.cpp scoprand.cpp settext.cpp symbol.cpp diff --git a/src/nrnoc/cabcode.cpp b/src/nrnoc/cabcode.cpp index cfd0faa111..5fdc65fd56 100644 --- a/src/nrnoc/cabcode.cpp +++ b/src/nrnoc/cabcode.cpp @@ -1,10 +1,12 @@ #include <../../nrnconf.h> /* /local/src/master/nrn/src/nrnoc/cabcode.cpp,v 1.37 1999/07/08 14:24:59 hines Exp */ -#define HOC_L_LIST 1 +#include #include #include #include + +#define HOC_L_LIST 1 #include "section.h" #include "nrn_ansi.h" #include "nrniv_mf.h" @@ -13,6 +15,36 @@ #include "hocparse.h" #include "membdef.h" +static char* escape_bracket(const char* s) { + static char* b; + const char* p1; + char* p2; + if (!b) { + b = new char[256]; + } + for (p1 = s, p2 = b; *p1; ++p1, ++p2) { + switch (*p1) { + case '<': + *p2 = '['; + break; + case '>': + *p2 = ']'; + break; + case '[': + case ']': + *p2 = '\\'; + *(++p2) = *p1; + break; + default: + *p2 = *p1; + break; + } + } + *p2 = '\0'; + return b; +} + + extern int hoc_execerror_messages; #define symlist hoc_symlist @@ -1979,8 +2011,8 @@ void forall_section(void) { Section* sec = hocSEC(qsec); qsec = qsec->next; if (buf[0]) { - hoc_regexp_compile(buf); - if (!hoc_regexp_search(secname(sec))) { + std::regex pattern(escape_bracket(buf)); + if (!std::regex_match(secname(sec), pattern)) { continue; } } @@ -2011,8 +2043,8 @@ void hoc_ifsec(void) { s = hoc_strpop(); Sprintf(buf, ".*%s.*", *s); - hoc_regexp_compile(buf); - if (hoc_regexp_search(secname(chk_access()))) { + std::regex pattern(escape_bracket(buf)); + if (std::regex_match(secname(chk_access()), pattern)) { hoc_execute(relative(savepc)); } if (!hoc_returning) @@ -2020,8 +2052,8 @@ void hoc_ifsec(void) { } void issection(void) { /* returns true if string is the access section */ - hoc_regexp_compile(gargstr(1)); - if (hoc_regexp_search(secname(chk_access()))) { + std::regex pattern(escape_bracket(gargstr(1))); + if (std::regex_match(secname(chk_access()), pattern)) { hoc_retpushx(1.); } else { hoc_retpushx(0.); diff --git a/src/oc/oc_ansi.h b/src/oc/oc_ansi.h index 34cc37062f..c83f9300fd 100644 --- a/src/oc/oc_ansi.h +++ b/src/oc/oc_ansi.h @@ -358,8 +358,6 @@ int is_vector_arg(int); char* vector_get_label(IvocVect*); void vector_set_label(IvocVect*, char*); -void hoc_regexp_compile(const char*); -int hoc_regexp_search(const char*); Symbol* hoc_install_var(const char*, double*); void hoc_class_registration(); void hoc_spinit(); diff --git a/src/oc/regexp.cpp b/src/oc/regexp.cpp deleted file mode 100644 index 4387fa1613..0000000000 --- a/src/oc/regexp.cpp +++ /dev/null @@ -1,466 +0,0 @@ -#include <../../nrnconf.h> -/* /local/src/master/nrn/src/oc/regexp.cpp,v 1.1.1.1 1994/10/12 17:22:13 hines Exp */ -/* -regexp.cpp,v - * Revision 1.1.1.1 1994/10/12 17:22:13 hines - * NEURON 3.0 distribution - * - * Revision 2.19 93/02/02 10:34:37 hines - * static functions declared before used - * - * Revision 1.3 92/07/31 12:11:31 hines - * following merged from hoc - * The regular expression has been augmented with - * {istart-iend} where istart and iend are integers. The expression matches - * any integer that falls in this range. - * - * Revision 1.2 92/01/30 08:17:19 hines - * bug fixes found in hoc incorporated. if()return, no else, objectcenter - * warnings. - * - * Revision 1.1 91/10/11 11:12:16 hines - * Initial revision - * - * Revision 3.108 90/10/24 09:44:14 hines - * saber warnings gone - * - * Revision 3.58 90/05/17 16:30:52 jamie - * changed global functions to start with hoc_ - * moved regexp.cpp from project 'neuron' to 'hoc' - * - * Revision 1.25 89/08/31 10:28:46 mlh - * regular expressions for issection() - * differences between standard regular expressions are: - * allways match from beginning to end of target (implicit ^ and $) - * change [] to <> - * eliminate \( - * - * Revision 1.2 89/08/31 09:22:17 mlh - * works as in e.cpp and lint free - * - * Revision 1.1 89/08/31 08:24:59 mlh - * Initial revision - * -*/ - -/* regular expression match for section names - grabbed prototype from e.cpp - Use by first compiling the search string with hoc_regexp_compile(pattern) - Then checking target strings one at a time with hoc_regexp_search(target) -*/ - -#include -#include "hocdec.h" -#define CABLESECTION 1 -/* Always match from beginning of string (implicit ^), - Always match end of string (implicit $), - change [] to <>, - eliminate \( -*/ - -#define STAR 01 -#define SUFF '.' -#define TILDE '~' - -#define EREGEXP 24 -#define error(enum) hoc_execerror("search string format error", pattern) -#define CBRA 1 -#define CCHR 2 -#define CDOT 4 -#define CCL 6 -#define NCCL 8 -#define CDOL 10 -#define CEOF 11 -#define CKET 12 -#if CABLESECTION -#define INTRANGE 14 -#endif -#define NBRA 5 -#define ESIZE 256 -#define eof '\0' -static char expbuf[ESIZE + 4]; -static const char* pattern = ""; -static char* loc1; -static char* loc2; -static char* locs; -static char* braslist[NBRA]; -static char* braelist[NBRA]; -static int circfl; -#if CABLESECTION -static int int_range_start[NBRA]; -static int int_range_stop[NBRA]; -#endif - -static int advance(char* lp, char* ep); -static int hoc_cclass(char* set, char c, int af); - -void hoc_regexp_compile(const char* pat) { - char* cp = (char*) pat; - int c; - char* ep; - char* lastep = 0; -#if (!CABLESECTION) - char bracket[NBRA], *bracketp; - int nbra; -#else - int int_range_index = 0; -#endif - int cclcnt; - int tempc; - - - if (!cp) { - pattern = ""; - error(EREGEXP); - } - if (pattern == cp && strcmp(pattern, cp)) { - /* if previous pattern != cp then may have been freed */ - return; - } - pattern = cp; - ep = expbuf; -#if (!CABLESECTION) - bracketp = bracket; - nbra = 0; -#endif - if ((c = *cp++) == '\n') { - cp--; - c = eof; - } - if (c == eof) { - if (*ep == 0) - error(EREGEXP); - return; - } -#if CABLESECTION - circfl = 1; -#else - circfl = 0; - if (c == '^') { - c = *cp++; - circfl++; - } -#endif - if (c == '*') - goto cerror; - cp--; - for (;;) { - if (ep >= &expbuf[ESIZE]) - goto cerror; - c = *cp++; - if (c == '\n') { - cp--; - c = eof; - } - if (c == eof) { -#if CABLESECTION - *ep++ = CDOL; -#endif - *ep++ = CEOF; - return; - } - if (c != '*') - lastep = ep; - switch (c) { - case '\\': -#if (!CABLESECTION) - if ((c = *cp++) == '(') { - if (nbra >= NBRA) - goto cerror; - *bracketp++ = nbra; - *ep++ = CBRA; - *ep++ = nbra++; - continue; - } - if (c == ')') { - if (bracketp <= bracket) - goto cerror; - *ep++ = CKET; - *ep++ = *--bracketp; - continue; - } -#endif - *ep++ = CCHR; - if (c == '\n') - goto cerror; - *ep++ = c; - continue; - - case '.': - *ep++ = CDOT; - continue; - - case '\n': - goto cerror; - - case '*': - if (*lastep == CBRA || *lastep == CKET) - error(EREGEXP); - *lastep |= STAR; - continue; - -#if (!CABLESECTION) - case '$': - tempc = *cp; - if (tempc != eof && tempc != '\n') - goto defchar; - *ep++ = CDOL; - continue; -#endif - -#if CABLESECTION - case '{': { - char* cp1 = cp; - if (int_range_index >= NBRA) - goto cerror; - *ep++ = INTRANGE; - do { - if (!(*cp >= '0' && *cp <= '9') && *cp != '-') { - error(EREGEXP); - } - } while (*(++cp) != '}'); - cp++; - if (2 != sscanf(cp1, - "%d-%d", - int_range_start + int_range_index, - int_range_stop + int_range_index)) { - error(EREGEXP); - } - *ep++ = int_range_index++; - } - continue; -#endif -#if CABLESECTION - case '<': -#else - case '[': -#endif - *ep++ = CCL; - *ep++ = 0; - cclcnt = 1; - if ((c = *cp++) == '^') { - c = *cp++; - ep[-2] = NCCL; - } - do { - if (c == '\n') - goto cerror; - /* - * Handle the escaped '-' - */ - if (c == '-' && *(ep - 1) == '\\') - *(ep - 1) = '-'; - /* - * Handle ranges of characters (e.g. a-z) - */ - else if ((tempc = *cp++) != ']' && c == '-' && cclcnt > 1 && tempc != '\n' && - (c = *(ep - 1)) <= tempc) { - while (++c <= tempc) { - *ep++ = c; - cclcnt++; - if (ep >= &expbuf[ESIZE]) - goto cerror; - } - } - /* - * Normal case. Add character to buffer - */ - else { - cp--; - *ep++ = c; - cclcnt++; - if (ep >= &expbuf[ESIZE]) - goto cerror; - } -#if CABLESECTION - } while ((c = *cp++) != '>'); -#else - - } while ((c = *cp++) != ']'); -#endif - lastep[1] = cclcnt; - continue; - -#if (!CABLESECTION) - defchar: -#endif - default: - *ep++ = CCHR; - *ep++ = c; - } - } -cerror: - expbuf[0] = 0; - error(EREGEXP); -} - -int hoc_regexp_search(const char* tar) { /*return true if target matches pattern*/ - char* target = (char*) tar; - char *p1, *p2, c; - -#if 1 - if (target == (char*) 0) { - return (0); - } - p1 = target; - locs = (char*) 0; -#else /* in e, apparently for searches within or at begining of string */ - if (gf) { - if (circfl) - return (0); - p1 = linebuf; - p2 = genbuf; - while (*p1++ = *p2++) - ; - locs = p1 = loc2; - } else { - if (addr == zero) - return (0); - p1 = getline(*addr); - locs = NULL; - } -#endif - p2 = expbuf; - if (circfl) { - loc1 = p1; - return (advance(p1, p2)); - } - /* fast check for first character */ - if (*p2 == CCHR) { - c = p2[1]; - do { - if (*p1 != c) - continue; - if (advance(p1, p2)) { - loc1 = p1; - return (1); - } - } while (*p1++); - return (0); - } - /* regular algorithm */ - do { - if (advance(p1, p2)) { - loc1 = p1; - return (1); - } - } while (*p1++); - return (0); -} - -static int advance(char* lp, char* ep) { - char* curlp; - - for (;;) - switch (*ep++) { - case CCHR: - if (*ep++ == *lp++) - continue; - return (0); - - case CDOT: - if (*lp++) - continue; - return (0); - - case CDOL: - if (*lp == 0) - continue; - return (0); - - case CEOF: - loc2 = lp; - return (1); - -#if CABLESECTION - case INTRANGE: { - int start, stop, num; - start = int_range_start[*ep]; - stop = int_range_stop[*ep++]; - num = *lp++ - '0'; - if (num < 0 || num > 9) { - return (0); - } - while (*lp >= '0' && *lp <= '9') { - num = 10 * num + *lp - '0'; - ++lp; - } - if (num >= start && num <= stop) { - continue; - } - } - return (0); -#endif - - case CCL: - if (hoc_cclass(ep, *lp++, 1)) { - ep += *ep; - continue; - } - return (0); - - case NCCL: - if (hoc_cclass(ep, *lp++, 0)) { - ep += *ep; - continue; - } - return (0); - - case CBRA: - braslist[*ep++] = lp; - continue; - - case CKET: - braelist[*ep++] = lp; - continue; - - case CDOT | STAR: - curlp = lp; - /*EMPTY*/ - while (*lp++) - ; - goto star; - - case CCHR | STAR: - curlp = lp; - /*EMPTY*/ - while (*lp++ == *ep) - ; - ep++; - goto star; - - case CCL | STAR: - case NCCL | STAR: - curlp = lp; - /*EMPTY*/ - while (hoc_cclass(ep, *lp++, ep[-1] == (CCL | STAR))) - ; - ep += *ep; - goto star; - - star: - do { - lp--; - if (lp == locs) - break; - if (advance(lp, ep)) - return (1); - } while (lp > curlp); - return (0); - - default: - error(EREGEXP); - } -} - -static int hoc_cclass(char* set, char c, int af) { - int n; - - if (c == 0) - return (0); - n = *set++; - while (--n) - if (*set++ == c) - return (af); - return (!af); -} From 94ce2eb3eb717c54941a8832eda858d4749f0b44 Mon Sep 17 00:00:00 2001 From: Erik Heeren Date: Fri, 2 Feb 2024 12:29:24 +0100 Subject: [PATCH 08/15] readthedocs: update config (#2703) --- .readthedocs.yml | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/.readthedocs.yml b/.readthedocs.yml index ea276443ce..dc4223186f 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -1,4 +1,9 @@ version: 2 +build: + os: "ubuntu-22.04" + tools: + python: "mambaforge-22.9" + conda: environment: docs/conda_environment.yml From ab5d4e513c7f9d388fdadd218ad6abb9a8b061c5 Mon Sep 17 00:00:00 2001 From: Erik Heeren Date: Mon, 5 Feb 2024 17:26:14 +0100 Subject: [PATCH 09/15] pin packaging to keep LegacyVersion support (#2706) * readthedocs: update config * pin packaging==21.3 - 22.0 dropped LegacyVersion: https://github.com/pypa/packaging/pull/407 - This makes it so our branches can no longer be parsed as versions --- docs/conda_environment.yml | 1 + docs/docs_requirements.txt | 1 + 2 files changed, 2 insertions(+) diff --git a/docs/conda_environment.yml b/docs/conda_environment.yml index f22bf80040..d13bdc8775 100644 --- a/docs/conda_environment.yml +++ b/docs/conda_environment.yml @@ -3,6 +3,7 @@ channels: - conda-forge - defaults dependencies: + - packaging=21.3 - python=3.11 - bison - cmake diff --git a/docs/docs_requirements.txt b/docs/docs_requirements.txt index b6639c9680..7a4ea6dd3f 100644 --- a/docs/docs_requirements.txt +++ b/docs/docs_requirements.txt @@ -15,3 +15,4 @@ plotly nbsphinx jinja2 sphinx-design +packaging==21.3 From 77130a86ae2fced55de8c47159ad3ba916c4a161 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 6 Feb 2024 09:24:47 +0100 Subject: [PATCH 10/15] Use new action for removing support of node16 to node20 (#2705) --- .github/workflows/coverage.yml | 10 +++++----- .github/workflows/docs.yml | 4 ++-- .github/workflows/formatting.yml | 2 +- .github/workflows/neuron-ci.yml | 8 ++++---- .github/workflows/release.yml | 2 +- .github/workflows/windows.yml | 6 +++--- 6 files changed, 16 insertions(+), 16 deletions(-) diff --git a/.github/workflows/coverage.yml b/.github/workflows/coverage.yml index 5c8a2dc9ed..d7c273564c 100644 --- a/.github/workflows/coverage.yml +++ b/.github/workflows/coverage.yml @@ -73,7 +73,7 @@ jobs: sudo apt-get install xvfb sudo /usr/bin/Xvfb $DISPLAY -screen 0 1600x1200x24 -noreset -nolock -shmem & # run in bg - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 with: fetch-depth: 2 @@ -83,7 +83,7 @@ jobs: git submodule update --init --recursive --force --depth 1 -- external/nmodl - name: Set up Python@${{ env.PY_MIN_VERSION }} - uses: actions/setup-python@v4 + uses: actions/setup-python@v5 with: python-version: ${{ env.PY_MIN_VERSION }} @@ -95,12 +95,12 @@ jobs: python -m pip install --upgrade -r ci_requirements.txt - name: Set up Python@${{ env.PY_MID_VERSION }} - uses: actions/setup-python@v4 + uses: actions/setup-python@v5 with: python-version: ${{ env.PY_MID_VERSION }} - name: Set up Python@${{ env.PY_MAX_VERSION }} - uses: actions/setup-python@v4 + uses: actions/setup-python@v5 with: python-version: ${{ env.PY_MAX_VERSION }} @@ -184,7 +184,7 @@ jobs: if: failure() && contains(github.event.pull_request.title, 'live-debug-coverage') uses: mxschmitt/action-tmate@v3 - - uses: codecov/codecov-action@v3 + - uses: codecov/codecov-action@v4 with: directory: ./build fail_ci_if_error: true diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index d4206d7345..ffc24629e8 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -36,11 +36,11 @@ jobs: shell: bash - name: Set up Python@${{ env.DEFAULT_PY_VERSION }} - uses: actions/setup-python@v4 + uses: actions/setup-python@v5 with: python-version: ${{ env.DEFAULT_PY_VERSION }} - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - name: Install Python dependencies working-directory: ${{runner.workspace}}/nrn diff --git a/.github/workflows/formatting.yml b/.github/workflows/formatting.yml index 1b490ed4b0..ab3e59723d 100644 --- a/.github/workflows/formatting.yml +++ b/.github/workflows/formatting.yml @@ -20,7 +20,7 @@ jobs: runs-on: ubuntu-22.04 timeout-minutes: 5 steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - name: Update submodule working-directory: ${{runner.workspace}}/nrn run: git submodule update --init external/coding-conventions diff --git a/.github/workflows/neuron-ci.yml b/.github/workflows/neuron-ci.yml index 0d0d05cc46..0c3ecc14a5 100644 --- a/.github/workflows/neuron-ci.yml +++ b/.github/workflows/neuron-ci.yml @@ -127,7 +127,7 @@ jobs: echo CI_OS_NAME=linux >> $GITHUB_ENV shell: bash - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 with: fetch-depth: 2 @@ -138,7 +138,7 @@ jobs: - name: Set up Python@${{ env.PY_MIN_VERSION }} if: ${{matrix.config.python_dynamic == 'ON'}} - uses: actions/setup-python@v4 + uses: actions/setup-python@v5 with: python-version: ${{ env.PY_MIN_VERSION }} @@ -151,7 +151,7 @@ jobs: python -m pip install --upgrade -r ci_requirements.txt - name: Set up Python@${{ env.PY_MAX_VERSION }} - uses: actions/setup-python@v4 + uses: actions/setup-python@v5 with: python-version: ${{ env.PY_MAX_VERSION }} @@ -198,7 +198,7 @@ jobs: # Workaround for https://github.com/actions/cache/issues/92 - name: Checkout cache action - uses: actions/checkout@v3 + uses: actions/checkout@v4 with: repository: actions/cache ref: v3 diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index 994c3c0028..3ba15345b3 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -25,7 +25,7 @@ jobs: release_url: ${{ steps.create_release.outputs.upload_url }} rel_tag: ${{ env.REL_TAG }} steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 name: Checkout branch ${{ env.REL_BRANCH }} with: ref: ${{ env.REL_BRANCH }} diff --git a/.github/workflows/windows.yml b/.github/workflows/windows.yml index fd9dcc0771..052a0d6679 100644 --- a/.github/workflows/windows.yml +++ b/.github/workflows/windows.yml @@ -35,7 +35,7 @@ jobs: timeout-minutes: 45 steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 with: ref: ${{ inputs.tag }} @@ -46,7 +46,7 @@ jobs: working-directory: ${{runner.workspace}}\nrn - name: Set up Python3 - uses: actions/setup-python@v4 + uses: actions/setup-python@v5 with: python-version: '3.11' @@ -79,7 +79,7 @@ jobs: uses: mxschmitt/action-tmate@v3 - name: Upload build artifact - uses: actions/upload-artifact@v3 + uses: actions/upload-artifact@v4 with: name: nrn-nightly-AMD64.exe path: ${{runner.workspace}}\nrn\nrn-nightly-AMD64.exe From c2566eac21dcc0fce1102da051c6bcaf7149d01c Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 6 Feb 2024 19:49:27 +0100 Subject: [PATCH 11/15] Cmake: Create 'get_link_libraries' to make it recursive (#2709) The issue was spotted when building with Caliper master. It returns MPI as a target and hence we need this fix. --- cmake/CMakeListsNrnMech.cmake | 94 ++++++++++++++++++++--------------- 1 file changed, 53 insertions(+), 41 deletions(-) diff --git a/cmake/CMakeListsNrnMech.cmake b/cmake/CMakeListsNrnMech.cmake index 38ea5084c8..6d53842273 100644 --- a/cmake/CMakeListsNrnMech.cmake +++ b/cmake/CMakeListsNrnMech.cmake @@ -23,48 +23,60 @@ endif() # Interview might have linked to libnrniv but we don't want to link to special list(REMOVE_ITEM NRN_LINK_LIBS "interviews") -# CMake does some magic to transform sys libs to -l. We replicate it -foreach(link_lib ${NRN_LINK_LIBS}) - # skip static readline library as it will be linked to nrniv (e.g. with wheel) also stub libraries - # from OSX can be skipped - if("${link_lib}" MATCHES "(libreadline.a|/*.tbd)") - continue() - endif() +function(get_link_libraries libs) + # CMake does some magic to transform sys libs to -l. We replicate it + foreach(link_lib ${libs}) + # skip static readline library as it will be linked to nrniv (e.g. with wheel) also stub + # libraries from OSX can be skipped + if("${link_lib}" MATCHES "(libreadline.a|/*.tbd)") + continue() + endif() - get_filename_component(dir_path ${link_lib} DIRECTORY) - if(TARGET ${link_lib}) - get_property( - link_flag - TARGET ${link_lib} - PROPERTY INTERFACE_LINK_LIBRARIES) - set(description - "Extracting link flags from target '${link_lib}', beware that this can be fragile.") - # Not use it yet because it can be generator expressions get_property(compile_flag TARGET - # ${link_lib} PROPERTY INTERFACE_COMPILE_OPTIONS) string(APPEND NRN_COMPILE_DEFS - # ${compile_flag}) - elseif(NOT dir_path) - set(link_flag "-l${link_lib}") - set(description - "Generating link flags from name '${link_lib}', beware that this can be fragile.") - # avoid library paths from special directory /nrnwheel which used to build wheels under docker - # container - elseif("${dir_path}" MATCHES "^/nrnwheel") - continue() - elseif("${dir_path}" MATCHES "^(/lib|/lib64|/usr/lib|/usr/lib64)$") - # NAME_WLE not avaialble with CMake version < 3.14 - get_filename_component(libname ${link_lib} NAME) - string(REGEX REPLACE "\\.[^.]*$" "" libname_wle ${libname}) - string(REGEX REPLACE "^lib" "" libname_wle ${libname_wle}) - set(link_flag "-l${libname_wle}") - set(description - "Extracting link flags from path '${link_lib}', beware that this can be fragile.") - else() - set(link_flag "${link_lib} -Wl,-rpath,${dir_path}") - set(description "Generating link flags from path ${link_lib}") - endif() - message(NOTICE "${description} Got: ${link_flag}") - string(APPEND NRN_LINK_DEFS " ${link_flag}") -endforeach() + get_filename_component(dir_path ${link_lib} DIRECTORY) + if(TARGET ${link_lib}) + get_property( + sublink_flag + TARGET ${link_lib} + PROPERTY INTERFACE_LINK_LIBRARIES) + set(description + "Extracting link flags from target '${link_lib}', beware that this can be fragile.") + # Not use it yet because it can be generator expressions get_property(compile_flag TARGET + # ${link_lib} PROPERTY INTERFACE_COMPILE_OPTIONS) string(APPEND NRN_COMPILE_DEFS + # ${compile_flag}) + foreach(sublink_lib ${sublink_flag}) + if(TARGET ${sublink_lib}) + message(NOTICE "For '${link_lib}' going to see TARGET '${sublink_lib}' recursively.") + get_link_libraries(${sublink_lib}) + else() + set(link_flag "${link_flag} ${sublink_flag}") + endif() + endforeach() + elseif(NOT dir_path) + set(link_flag "-l${link_lib}") + set(description + "Generating link flags from name '${link_lib}', beware that this can be fragile.") + # avoid library paths from special directory /nrnwheel which used to build wheels under docker + # container + elseif("${dir_path}" MATCHES "^/nrnwheel") + continue() + elseif("${dir_path}" MATCHES "^(/lib|/lib64|/usr/lib|/usr/lib64)$") + # NAME_WLE not avaialble with CMake version < 3.14 + get_filename_component(libname ${link_lib} NAME) + string(REGEX REPLACE "\\.[^.]*$" "" libname_wle ${libname}) + string(REGEX REPLACE "^lib" "" libname_wle ${libname_wle}) + set(link_flag "-l${libname_wle}") + set(description + "Extracting link flags from path '${link_lib}', beware that this can be fragile.") + else() + set(link_flag "${link_lib} -Wl,-rpath,${dir_path}") + set(description "Generating link flags from path ${link_lib}") + endif() + message(NOTICE "${description} Got: ${link_flag}") + string(APPEND NRN_LINK_DEFS " ${link_flag}") + endforeach() +endfunction(get_link_libraries) + +get_link_libraries("${NRN_LINK_LIBS}") # Compiler flags depending on cmake build type from BUILD_TYPE__FLAGS string(TOUPPER "${CMAKE_BUILD_TYPE}" _BUILD_TYPE) From 758d15c8e39e5dc18d8f5ba91ff027ad232dc0e1 Mon Sep 17 00:00:00 2001 From: nrnhines Date: Thu, 8 Feb 2024 05:57:58 -0500 Subject: [PATCH 12/15] Support for RANDOM language construct in NMODL (#2627) * Added support for RANDOM construct parsing * Updated netstim.mod to use new RANDOM construct * Add tests for usage of RANDOM in HOC and Python * Documentation for NEURON { RANDOM123 ranvar} and NMODLRandom * Update NMODL to support CoreNEURON side * Update NetStim with RANDOM but allow legacy API. * Update CoreNEURON <-> NEURON for NMODL RANDOM Co-authored-by: Pramod S Kumbhar Co-authored-by: Omar Awile --- cmake/NeuronFileLists.cmake | 1 + docs/cmake_doc/options.rst | 14 + .../programmatic/mechanisms/nmodl2.rst | 50 ++ docs/python/programming/math/random.rst | 89 ++++ docs/python/simctrl/bbsavestate.rst | 79 ++- external/nmodl | 2 +- src/coreneuron/io/core2nrn_data_return.cpp | 65 +++ src/coreneuron/io/nrn2core_direct.h | 1 + src/coreneuron/io/nrn_checkpoint.cpp | 33 +- src/coreneuron/io/nrn_setup.cpp | 18 +- src/coreneuron/io/phase2.cpp | 44 ++ src/coreneuron/io/phase2.hpp | 1 + .../mechanism/mech/modfile/netstim.mod | 477 +++++------------- src/coreneuron/mechanism/register_mech.cpp | 2 + src/coreneuron/nrniv/nrniv_decl.h | 1 + src/coreneuron/utils/nrnoc_aux.cpp | 2 +- src/coreneuron/utils/randoms/nrnran123.h | 25 + src/gnu/nrnran123.cpp | 66 ++- src/gnu/nrnran123.h | 79 ++- src/modlunit/extdef.h | 3 +- src/modlunit/init.cpp | 1 + src/modlunit/model.h | 1 + src/modlunit/nrnunit.cpp | 6 + src/modlunit/parse1.ypp | 3 + src/nmodl/init.cpp | 16 + src/nmodl/modl.h | 21 +- src/nmodl/nocpout.cpp | 84 ++- src/nmodl/parsact.cpp | 6 +- src/nmodl/parse1.ypp | 19 + src/nrniv/bbsavestate.cpp | 33 +- src/nrniv/ndatclas.cpp | 2 +- src/nrniv/nmodlrandom.cpp | 141 ++++++ src/nrniv/nrnclass.h | 6 +- .../callbacks/nrncore_callbacks.cpp | 62 ++- .../callbacks/nrncore_callbacks.h | 8 + src/nrniv/nrncore_write/data/cell_group.cpp | 5 +- src/nrniv/nrncore_write/io/nrncore_io.cpp | 17 +- src/nrniv/nrncore_write/io/nrncore_io.h | 2 + src/nrniv/nrnmenu.cpp | 6 +- src/nrnoc/cabcode.cpp | 18 + src/nrnoc/init.cpp | 131 +++-- src/nrnoc/membfunc.h | 10 +- src/nrnoc/netstim.mod | 310 +++--------- src/nrnoc/point.cpp | 3 + src/nrnoc/treeset.cpp | 3 + src/nrnpython/nrnpy_hoc.cpp | 2 +- src/nrnpython/nrnpy_nrn.cpp | 15 +- src/oc/code.h | 1 + src/oc/hoc_oop.cpp | 27 + src/oc/oc_ansi.h | 3 + src/oc/parse.ypp | 19 +- test/CMakeLists.txt | 25 + test/coreneuron/mod files/noisychan.mod | 59 +++ test/coreneuron/mod files/watchrange.mod | 17 + test/coreneuron/mod files/watchrange2.mod | 76 +++ test/coreneuron/test_nmodlrandom.py | 169 +++++++ test/coreneuron/test_nmodlrandom_syntax.py | 96 ++++ test/coreneuron/test_watchrange.py | 27 +- test/external/CMakeLists.txt | 2 +- test/hoctests/rand.mod | 43 ++ test/hoctests/rand_art.mod | 50 ++ test/hoctests/rand_pp.mod | 43 ++ test/hoctests/tests/test_random.hoc | 79 +++ test/hoctests/tests/test_random.json | 175 +++++++ test/hoctests/tests/test_random.py | 178 +++++++ test/nmodl/test_random.py | 129 +++++ test/rxd/test_currents.py | 7 +- test/rxd/test_pure_diffusion.py | 7 +- 68 files changed, 2491 insertions(+), 724 deletions(-) create mode 100644 src/nrniv/nmodlrandom.cpp create mode 100644 test/coreneuron/mod files/noisychan.mod create mode 100644 test/coreneuron/mod files/watchrange2.mod create mode 100644 test/coreneuron/test_nmodlrandom.py create mode 100644 test/coreneuron/test_nmodlrandom_syntax.py create mode 100644 test/hoctests/rand.mod create mode 100644 test/hoctests/rand_art.mod create mode 100644 test/hoctests/rand_pp.mod create mode 100644 test/hoctests/tests/test_random.hoc create mode 100644 test/hoctests/tests/test_random.json create mode 100644 test/hoctests/tests/test_random.py create mode 100644 test/nmodl/test_random.py diff --git a/cmake/NeuronFileLists.cmake b/cmake/NeuronFileLists.cmake index fb6ca521e7..9141aa4562 100644 --- a/cmake/NeuronFileLists.cmake +++ b/cmake/NeuronFileLists.cmake @@ -228,6 +228,7 @@ set(NRNIV_FILE_LIST multisplit.cpp ndatclas.cpp netpar.cpp + nmodlrandom.cpp nonlinz.cpp nrncore_write.cpp nrncore_write/callbacks/nrncore_callbacks.cpp diff --git a/docs/cmake_doc/options.rst b/docs/cmake_doc/options.rst index 83497d5c0b..1d9bfba9e6 100644 --- a/docs/cmake_doc/options.rst +++ b/docs/cmake_doc/options.rst @@ -579,6 +579,15 @@ NRN_COVERAGE_FILES:STRING= ``-DNRN_COVERAGE_FILES="src/nrniv/partrans.cpp;src/nmodl/parsact.cpp;src/nrnpython/nrnpy_hoc.cpp"`` + For a list of all the cpp files changed in a pull request, consider + copy/pasting the ``;`` separated list obtained with + + .. code-block:: shell + + a=`git diff --name-only master | grep '\.cpp'` + echo $a | sed 's/ /;/g' + + NRN_SANITIZERS:STRING= ---------------------- Enable some combination of AddressSanitizer, LeakSanitizer, ThreadSanitizer @@ -589,6 +598,11 @@ NRN_SANITIZERS:STRING= ``-DNRN_SANITIZERS=address`` with the use of Python virtual environments; if you attempt this then the CMake code should recommend a solution. + Note: the ``address`` sanitizer also prints leak infornation when a + launch exits. That can be avoided with + + ``export ASAN_OPTIONS=detect_leaks=0`` + Miscellaneous Rarely used options specific to NEURON: ===================================================== diff --git a/docs/python/modelspec/programmatic/mechanisms/nmodl2.rst b/docs/python/modelspec/programmatic/mechanisms/nmodl2.rst index dc949f32ec..e5ad80727d 100644 --- a/docs/python/modelspec/programmatic/mechanisms/nmodl2.rst +++ b/docs/python/modelspec/programmatic/mechanisms/nmodl2.rst @@ -44,6 +44,7 @@ Description: POINT_PROCESS ... POINTER pointer1, ... BBCOREPOINTER bbcore1, ... + RANDOM ranvar1, ... EXTERNAL external1, ... THREADSAFE REPRESENTS ontology_id @@ -410,6 +411,55 @@ Description: ``TODO``: Add description (?) and existing example mod file (provided by link) +.. _nmodlrandom: + +RANDOM +###### + +Description: + .. code-block:: + + NEURON { + RANDOM ranvar1, ... + } + + These names refer to random variable streams that are automatically + associated with nrnran123 generators. Such nrnran123 generators are also used, for example to implement + :meth:`Random.Random123` + These names are analogous to range variables in that the streams are distinct for every mechanism instance + of a POINT_PROCESS, ARTIFICIAL_CELL, or instance of a density mechanism in a segment of a cable section. + Each stream exists for the lifetime of the mechanism instance. While a stream exists, its properties can + be changed from the interpreter. + + Prior to the introduction of this keyword, random streams required a POINTER variable and + fairly elaborate VERBATIM blocks + to setup the streams and manage the stream properties from HOC or Python so that each stream was + statistically independent of all other streams. + + From the interpreter, the ranvar1 stream properties are assigned and evaluated using standard + range variable syntax where mention of ranvar1 returns a :class:`~NMODLRandom` object that wraps the stream + and provides method calls to get and set the three stream ids and the starting sequence number. + + When a stream is instantiated, its identifier triplet is default initialized to + (1, :meth:`mpiworldrank `, ++internal_id3) + so all streams are statistically independent (at launch time, internal_id3 = 0). + However since the identifier triplet depends on the order of + construction, it is recommended for parallel simulation reproducibility that triplets be algorithmically specified + at the interpreter level. And see :meth:`Random.Random123_globalindex`. + + At present, the list of random_... methods available for use within mod files (outside of VERBATIM blocks) are: + + * random_setseq(ranvar1, uint34_value) + * random_setids(ranvar1, id1_uint32, id2_uint32, id3_uint32) + * x = random_uniform(ranvar1) : uniform 0 to 1 -- minimum value is 2.3283064e-10 and max value is 1-min + * x = random_uniform(ranvar1, min, max) + * x = random_negexp(ranvar1) : mean 1.0 -- min value is 2.3283064e-10, max is 22.18071 + * x = random_negexp(ranvar1, mean) + * x = random_normal(ranvar1) : mean 1.0, std 1.0 + * x = random_normal(ranvar1, mean, std) + * x = random_ipick(ranvar1) : range 0 to 2^32-1 + * x = random_dpick(ranvar1) + EXTERNAL ######## diff --git a/docs/python/programming/math/random.rst b/docs/python/programming/math/random.rst index dfaea09a04..210b43500d 100755 --- a/docs/python/programming/math/random.rst +++ b/docs/python/programming/math/random.rst @@ -768,3 +768,92 @@ Random Class +---- + +NMODLRandom Class +================= + +.. class:: NMODLRandom + + Syntax: + ``r = point_process.ranvar`` + + ``r = section(x).mech.ranvar`` + + ``r = section(x).ranvar_mech`` + + + Description: + Returns an NMODLRandom wrapper for the nrnran123_State associated with the mechanism + :ref:`RANDOM ranvar ` variable. + Note that an attempt to assign a value to ranvar will raise an error. + At present, all mentions of ranvar in the context of a specific mechanism instance return a wrapper for + the same nrnran123_State (though the NMODLRandom instances are different). + +---- + +.. method:: NMODLRandom.get_ids + + Syntax: + ``vector = r.get_ids()`` + + Description: + Returns a HOC Vector of size 3 containing the 32 bit id1, id2, id3 of the nrnran123_State + +---- + +.. method:: NMODLRandom.set_ids + + Syntax: + ``r = r.set_ids(id1, id2, id3)`` + + Description: + Sets the 32 bit id1, id2, id3 of the nrnran123_State and returns the same NModlRandom instance. + + +---- + +.. method:: NMODLRandom.get_seq + + Syntax: + ``x = r.get_seq()`` + + Description: + Returns as a float, the 34 bit sequence position of the nrnran123_State + +---- + +.. method:: NMODLRandom.set_seq + + Syntax: + ``r = r.set_seq(x)`` + + Description: + Sets the 34 bit sequence position of the nrnran123_State. Returns the same NMODLRandom instance. + +---- + +.. method:: NMODLRandom.uniform + + Syntax: + ``x = r.uniform()`` + + Description: + Returns as a float, the uniform random value in the open interval 0 to 1 at the current sequence + position of the nrnran123_State (the current sequence position is then incremented by 1) + This is, for testing purposes, the only distribution exposed to + the interpreter. We don't forsee any practical use of + NMODLRandom within the interpreter in regard to sampling. The purpose + of NMODLRandom is to allow setting of stream properties for a + mod file RANDOM variable. Indeed, if one explicitly constructs an NMODLRandom + from the interpreter, then + + .. code-block:: + python + + from neuron import h + r = h.NMODLRandom() + print(r.uniform()) + + NEURON: NMODLRandom wrapped handle is not valid + diff --git a/docs/python/simctrl/bbsavestate.rst b/docs/python/simctrl/bbsavestate.rst index 8d4b627f82..cb22f24995 100644 --- a/docs/python/simctrl/bbsavestate.rst +++ b/docs/python/simctrl/bbsavestate.rst @@ -30,18 +30,24 @@ BBSaveState h.stdinit() bbss = h.BBSaveState() if restore: - bbss.restore_test() + bbss.restore("temp.dat") print(f'after restore t={h.t}') else: pc.psolve(tstop/2) - bbss.save_test() + bbss.save("temp.dat") pc.psolve(tstop) - Note that files are saved in a subdirectory called "out" and restored - from a subdirectory called "in". An empty "out" folder should be created by - the user prior to calling save_test(). A script filter + In this case, the entire model state for this MPI rank is in the filename for save and restore. + + If multisplit is involved, or it is desired to reassemble the model cells on a different set of MPI ranks, + one instead must use :meth:`BBSaveState.save_test` + and :meth:`BBSaveState.restore_test`. This allows reassembly of the multisplit subtrees back into + their complete cells to allow different multisplitting and different cell distribution on different ranks. + In this case + files are saved in a subdirectory called "bbss_out" and restored + from a subdirectory called "bbss_in". A script filter (see :meth:`BBSaveState.save_test`) is needed to copy and sometimes - concatenate files from the out to the in subfolders. These files have + concatenate files from the bbss_out to the bbss_in subfolders. These files have an ascii format. BBSaveState has a c++ API that allows one to replace the file reader and @@ -59,7 +65,7 @@ BBSaveState Because a restore clears the event queue and because one cannot call finitialize from hoc without vitiating the restore, :meth:`Vector.play` will not work unless one calls :meth:`BBSaveState.vector_play_init` after a - restore (similarly :func:`frecord` must be called for :meth:`Vector.record` to work. + restore (similarly :func:`frecord_init` must be called for :meth:`Vector.record` to work. Note that it is necessary that Vector.play use a tvec argument with a first element greater than or equal to the restore time. @@ -72,9 +78,10 @@ BBSaveState with a base gid. 4. NetCon.event in Hoc can be used only with NetCon's with a None source. - - To allow extra state, such as Random sequence, to be saved for - POINT_PROCESS or SUFFIX density nmodl mechanisms, + RANDOM variables declared in a mod file NEURON block have their sequence value saved + automatically. + + To allow extra state to be saved, eg. when POINTER and BBCOREPOINTER are used to manage objects, declare FUNCTION bbsavestate() within the mechanism. That function is called when the mechanism instance is saved and restored. @@ -119,16 +126,16 @@ BBSaveState Description: - State of the model is saved in files within the subdirectory, `out`. - The file `out/tmp` contains the value of t. Other files have the + State of the model is saved in files within the subdirectory, `bbss_out`. + The file `bbss_out/tmp` contains the value of t. Other files have the filename format tmp.. . Only in the case of multisplit is it possible to have the same gid in more than one filename. Note that the out folder needs to be created by the user prior to a call to save_test(). To prepare for a restore, the tmp.. files should be copied - from the `out` subfolder to a subfolder called `in`, with the filename - in/tmp. . Each file should begin with a first line that specifies + from the `bbss_out` subfolder to a subfolder called `bbss_in`, with the filename + bbss_in/tmp. . Each file should begin with a first line that specifies the number of files in the `out` folder that had the same gid. The following out2in.sh script shows how to do this (not particularly @@ -138,16 +145,16 @@ BBSaveState bash #!/usr/bin/env bash - rm -f in/* - cat out/tmp > in/tmp - for f in out/tmp.*.* ; do + rm -f bbss_in/* + cat bbss_out/tmp > bbss_in/tmp + for f in bbss_out/tmp.*.* ; do echo $f i=`echo "$f" | sed 's/.*tmp\.\([0-9]*\)\..*/\1/'` echo $i - if test ! -f in/tmp.$i ; then - cnt=`ls out/tmp.$i.* | wc -l` - echo $cnt > in/tmp.$i - cat out/tmp.$i.* >> in/tmp.$i + if test ! -f bbss_in/tmp.$i ; then + cnt=`ls bbss_out/tmp.$i.* | wc -l` + echo $cnt > bbss_in/tmp.$i + cat bbss_out/tmp.$i.* >> bbss_in/tmp.$i fi done @@ -166,16 +173,40 @@ BBSaveState Description: State of the model is restored from files within the - subdirectory, "in". The file "in/tmp" supplies the value of t. - Other files have the filename format tmp. and are read when + subdirectory, "bbss_in". The file "bbss_in/tmp" supplies the value of t. + Other files have the filename format tmp. and are read when that gid is restored. Note that in a multisplit context, the same - "in/tmp." file will be read by multiple ranks, but only the state + "bbss_in/tmp." file will be read by multiple ranks, but only the state assocated with sections that exist on a rank will be restored. ---- +.. method:: BBSaveState.save + + Syntax: + ``.save("filename")`` + Description: + Saves the state of the entire model (on this rank). This is simpler to use than the ``save_test``, ``restore_test`` + pattern but does not work if one has multisplit cells or desires a different distribution of cells on a different + number of ranks. + + +---- + + +.. method:: BBSaveState.restore + + Syntax: + ``.restore("filename")`` + + Description: + Restores the state of the entire model (on this rank). This is simpler to use than the ``save_test``, ``restore_test`` + pattern but does not work if one has multisplit cells or desires a different distribution of cells on a different + number of ranks. + +---- .. method:: BBSaveState.ignore diff --git a/external/nmodl b/external/nmodl index bb4dfd2fbc..8f7eb99fd3 160000 --- a/external/nmodl +++ b/external/nmodl @@ -1 +1 @@ -Subproject commit bb4dfd2fbc5209bb1893d3c681157db743ca1dfc +Subproject commit 8f7eb99fd36ab886eac5c1ab050272fd2c46fa04 diff --git a/src/coreneuron/io/core2nrn_data_return.cpp b/src/coreneuron/io/core2nrn_data_return.cpp index 623dc1f477..5b13d4880d 100644 --- a/src/coreneuron/io/core2nrn_data_return.cpp +++ b/src/coreneuron/io/core2nrn_data_return.cpp @@ -7,6 +7,7 @@ */ #include +#include #include "coreneuron/coreneuron.hpp" #include "coreneuron/io/nrn2core_direct.h" @@ -40,6 +41,11 @@ int (*core2nrn_corepointer_mech_)(int tid, int dcnt, int* iArray, double* dArray); + +int (*core2nrn_nmodlrandom_)(int tid, + int type, + const std::vector& indices, + const std::vector& nmodlrandom); } namespace coreneuron { @@ -179,6 +185,64 @@ static void core2nrn_corepointer(int tid, NrnThreadMembList* tml) { (*core2nrn_corepointer_mech_)(tid, type, icnt, dcnt, iArray.get(), dArray.get()); } +// based on code from nrncore_callbacks.cpp +std::vector& nrn_mech_random_indices(int type) { + static std::unordered_map> mech_random_indices{}; + static std::mutex mx; + std::unique_lock lock(mx); + if (mech_random_indices.count(type) == 0) { + // if no element, create empty one and search dparam_semantics to fill + auto& mri = mech_random_indices[type]; + int* semantics = corenrn.get_memb_func(type).dparam_semantics; + int dparam_size = corenrn.get_prop_dparam_size()[type]; + for (int i = 0; i < dparam_size; ++i) { + if (semantics[i] == -11) { + mri.push_back(i); + } + } + } + lock.unlock(); + return mech_random_indices[type]; +} + +/** @brief Copy back NMODL RANDOM sequence to NEURON + */ +static void c2n_nmodlrandom(int tid, NrnThreadMembList* tml) { + // Started out as a copy of corenrn_corepointer above. + // overall algorithm for nmodlrandom is similar to nrnthread_dat2_mech. + int type = tml->index; + auto& indices = nrn_mech_random_indices(type); + if (indices.size() == 0) { + return; + } + NrnThread& nt = nrn_threads[tid]; + Memb_list* ml = tml->ml; + int layout = corenrn.get_mech_data_layout()[type]; + int pdsz = corenrn.get_prop_dparam_size()[type]; + int aln_cntml = nrn_soa_padded_size(ml->nodecount, layout); + int n = ml->nodecount; + + // will send back vector of 34 bit uints (aka double) + std::vector nmodlrandom{}; + nmodlrandom.reserve(n * indices.size()); + for (int ix: indices) { + for (int j = 0; j < n; ++j) { + int jp = j; + if (ml->_permute) { + jp = ml->_permute[j]; + } + int pv = ml->pdata[nrn_i_layout(jp, n, ix, pdsz, layout)]; + nrnran123_State* state = (nrnran123_State*) nt._vdata[pv]; + uint32_t seq; + char which; + nrnran123_getseq(state, &seq, &which); + nmodlrandom.push_back(double(seq) * 4 + which); + } + } + + (*core2nrn_nmodlrandom_)(tid, type, indices, nmodlrandom); +} + /** @brief Copy event queue and related state back to NEURON. */ static void core2nrn_tqueue(NrnThread&); @@ -275,6 +339,7 @@ void core2nrn_data_return() { } core2nrn_corepointer(tid, tml); + c2n_nmodlrandom(tid, tml); } // Copy the event queue and related state. diff --git a/src/coreneuron/io/nrn2core_direct.h b/src/coreneuron/io/nrn2core_direct.h index 64c6930a53..5790f2197d 100644 --- a/src/coreneuron/io/nrn2core_direct.h +++ b/src/coreneuron/io/nrn2core_direct.h @@ -59,6 +59,7 @@ extern int (*nrn2core_get_dat2_mech_)(int tid, int*& nodeindices, double*& data, int*& pdata, + std::vector& nmodlrandom, std::vector& pointer2type); extern int (*nrn2core_get_dat2_3_)(int tid, diff --git a/src/coreneuron/io/nrn_checkpoint.cpp b/src/coreneuron/io/nrn_checkpoint.cpp index 77e301b4e5..ed91d7de81 100644 --- a/src/coreneuron/io/nrn_checkpoint.cpp +++ b/src/coreneuron/io/nrn_checkpoint.cpp @@ -73,6 +73,8 @@ void CheckPoints::write_checkpoint(NrnThread* nt, int nb_threads) const { /** * if openmp threading needed: * #pragma omp parallel for private(i) shared(nt, nb_threads) schedule(runtime) + * but note that nrn_mech_random_indices(type) is not threadsafe on first + * call for each type. */ for (int i = 0; i < nb_threads; i++) { if (nt[i].ncell || nt[i].tml) { @@ -307,12 +309,41 @@ void CheckPoints::write_phase2(NrnThread& nt) const { } } fh.write_array(d, cnt * sz); - delete[] d; size_t s = pointer2type.size(); fh << s << " npointer\n"; if (s) { fh.write_array(pointer2type.data(), s); } + + // nmodlrandom + auto& indices = nrn_mech_random_indices(type); + s = indices.size() ? (1 + indices.size() + 5 * cnt * indices.size()) : 0; + fh << s << " nmodlrandom\n"; + if (s) { + std::vector nmodlrandom{}; + nmodlrandom.reserve(s); + nmodlrandom.push_back((uint32_t) indices.size()); + for (auto ix: indices) { + nmodlrandom.push_back((uint32_t) ix); + } + for (auto ix: indices) { + uint32_t data[5]; + char which; + for (int i = 0; i < cnt; ++i) { + void* v = nt._vdata[d[i * sz + ix]]; + nrnran123_State* r = (nrnran123_State*) v; + nrnran123_getids3(r, &data[0], &data[1], &data[2]); + nrnran123_getseq(r, &data[3], &which); + data[4] = uint32_t(which); + for (auto j: data) { + nmodlrandom.push_back(j); + } + } + } + fh.write_array(nmodlrandom.data(), nmodlrandom.size()); + } + + delete[] d; } } diff --git a/src/coreneuron/io/nrn_setup.cpp b/src/coreneuron/io/nrn_setup.cpp index c346f84bf5..8697a9015d 100644 --- a/src/coreneuron/io/nrn_setup.cpp +++ b/src/coreneuron/io/nrn_setup.cpp @@ -738,6 +738,16 @@ void nrn_cleanup() { (*s)(nt, ml, tml->index); } + // Moved from below as priv_dtor is now deleting the RANDOM streams, + // and at this moment need an undeleted pdata. + // Destroy the global variables struct allocated in nrn_init + if (auto* const priv_dtor = corenrn.get_memb_func(tml->index).private_destructor) { + (*priv_dtor)(nt, ml, tml->index); + assert(!ml->instance); + assert(!ml->global_variables); + assert(ml->global_variables_size == 0); + } + ml->data = nullptr; // this was pointing into memory owned by nt free_memory(ml->pdata); ml->pdata = nullptr; @@ -753,14 +763,6 @@ void nrn_cleanup() { ml->_thread = nullptr; } - // Destroy the global variables struct allocated in nrn_init - if (auto* const priv_dtor = corenrn.get_memb_func(tml->index).private_destructor) { - (*priv_dtor)(nt, ml, tml->index); - assert(!ml->instance); - assert(!ml->global_variables); - assert(ml->global_variables_size == 0); - } - NetReceiveBuffer_t* nrb = ml->_net_receive_buffer; if (nrb) { if (nrb->_size) { diff --git a/src/coreneuron/io/phase2.cpp b/src/coreneuron/io/phase2.cpp index 54e59ee37b..0502129068 100644 --- a/src/coreneuron/io/phase2.cpp +++ b/src/coreneuron/io/phase2.cpp @@ -50,6 +50,7 @@ int (*nrn2core_get_dat2_mech_)(int tid, int*& nodeindices, double*& data, int*& pdata, + std::vector& nmodlrandom, std::vector& pointer2type); int (*nrn2core_get_dat2_3_)(int tid, @@ -182,6 +183,13 @@ void Phase2::read_file(FileHandler& F, const NrnThread& nt) { auto& p2t = tmls.back().pointer2type; p2t = F.read_vector(sz); } + + // nmodlrandom + sz = F.read_int(); + if (sz) { + auto& nmodlrandom = tmls.back().nmodlrandom; + nmodlrandom = F.read_vector(sz); + } } } output_vindex = F.read_vector(nt.n_presyn); @@ -325,13 +333,24 @@ void Phase2::read_direct(int thread_id, const NrnThread& nt) { int* nodeindices_ = nullptr; double* data_ = _data + offset; int* pdata_ = const_cast(tml.pdata.data()); + + // nmodlrandom, receives: + // number of random variables + // dparam index (neuron side) of each random variable + // 5 uint32 for each var of each instance + // id1, id2, id3, seq, uint32_t(which) + // all instances of ranvar1 first, then all instances of ranvar2, etc. + auto& nmodlrandom = tml.nmodlrandom; + (*nrn2core_get_dat2_mech_)(thread_id, i, dparam_sizes[type] > 0 ? dsz_inst : 0, nodeindices_, data_, pdata_, + nmodlrandom, tml.pointer2type); + if (dparam_sizes[type] > 0) dsz_inst++; offset += nrn_soa_padded_size(nodecounts[i], layout) * param_sizes[type]; @@ -1111,6 +1130,31 @@ void Phase2::populate(NrnThread& nt, const UserParams& userParams) { pp->_tid = nt.id; } } + + auto& r = tmls[itml].nmodlrandom; + if (r.size()) { + size_t ix{}; + uint32_t n_randomvar = r[ix++]; + assert(r.size() == 1 + n_randomvar + 5 * n_randomvar * n); + std::vector indices(n_randomvar); + for (uint32_t i = 0; i < n_randomvar; ++i) { + indices[i] = r[ix++]; + } + int cnt = ml->nodecount; + for (auto index: indices) { + // should we also verify that index on corenrn side same as on nrn side? + // sonarcloud thinks ml_pdata can be nullptr, so ... + assert(index >= 0 && index < szdp); + for (int i = 0; i < n; ++i) { + nrnran123_State* state = nrnran123_newstream3(r[ix], r[ix + 1], r[ix + 2]); + nrnran123_setseq(state, r[ix + 3], char(r[ix + 4])); + ix += 5; + int ipd = ml->pdata[nrn_i_layout(i, cnt, index, szdp, layout)]; + assert(ipd >= 0 && ipd < n_vdata + extra_nv); + nt._vdata[ipd] = state; + } + } + } } // pnt_offset needed for SelfEvent transfer from NEURON. Not needed on GPU. diff --git a/src/coreneuron/io/phase2.hpp b/src/coreneuron/io/phase2.hpp index 8b84305e43..9071e7ee38 100644 --- a/src/coreneuron/io/phase2.hpp +++ b/src/coreneuron/io/phase2.hpp @@ -115,6 +115,7 @@ class Phase2 { std::vector iArray; std::vector dArray; std::vector pointer2type; + std::vector nmodlrandom{}; }; std::vector tmls; std::vector output_vindex; diff --git a/src/coreneuron/mechanism/mech/modfile/netstim.mod b/src/coreneuron/mechanism/mech/modfile/netstim.mod index 1ec52940ac..86bb7562fe 100644 --- a/src/coreneuron/mechanism/mech/modfile/netstim.mod +++ b/src/coreneuron/mechanism/mech/modfile/netstim.mod @@ -1,377 +1,134 @@ : $Id: netstim.mod 2212 2008-09-08 14:32:26Z hines $ : comments at end -: the Random idiom has been extended to support CoreNEURON. - -: For backward compatibility, noiseFromRandom(hocRandom) can still be used -: as well as the default low-quality scop_exprand generator. -: However, CoreNEURON will not accept usage of the low-quality generator, -: and, if noiseFromRandom is used to specify the random stream, that stream -: must be using the Random123 generator. - -: The recommended idiom for specfication of the random stream is to use -: noiseFromRandom123(id1, id2[, id3]) - -: If any instance uses noiseFromRandom123, then no instance can use noiseFromRandom -: and vice versa. - NEURON { - ARTIFICIAL_CELL NetStim - RANGE interval, number, start - RANGE noise - THREADSAFE : only true if every instance has its own distinct Random - BBCOREPOINTER donotuse + ARTIFICIAL_CELL NetStim + THREADSAFE + RANGE interval, number, start + RANGE noise + RANDOM ranvar } PARAMETER { - interval = 10 (ms) <1e-9,1e9>: time between spikes (msec) - number = 10 <0,1e9> : number of spikes (independent of noise) - start = 50 (ms) : start of first spike - noise = 0 <0,1> : amount of randomness (0.0 - 1.0) + interval = 10 (ms) <1e-9,1e9> : time between spikes (msec) + number = 10 <0,1e9> : number of spikes (independent of noise) + start = 50 (ms) : start of first spike + noise = 0 <0,1> : amount of randomness (0.0 - 1.0) } ASSIGNED { - event (ms) - on - ispike - donotuse -} - -VERBATIM -#if NRNBBCORE /* running in CoreNEURON */ - -#define IFNEWSTYLE(arg) arg - -#else /* running in NEURON */ - -/* - 1 means noiseFromRandom was called when _ran_compat was previously 0 . - 2 means noiseFromRandom123 was called when _ran_compat was previously 0. -*/ -static int _ran_compat; /* specifies the noise style for all instances */ -#define IFNEWSTYLE(arg) if(_ran_compat == 2) { arg } - -#endif /* running in NEURON */ -ENDVERBATIM - -:backward compatibility -PROCEDURE seed(x) { -VERBATIM -#if !NRNBBCORE -ENDVERBATIM - set_seed(x) -VERBATIM -#endif -ENDVERBATIM + event (ms) + on + ispike } INITIAL { - - VERBATIM - if (_p_donotuse) { - /* only this style initializes the stream on finitialize */ - IFNEWSTYLE(nrnran123_setseq((nrnran123_State*)_p_donotuse, 0, 0);) - } - ENDVERBATIM - - on = 0 : off - ispike = 0 - if (noise < 0) { - noise = 0 - } - if (noise > 1) { - noise = 1 - } - if (start >= 0 && number > 0) { - on = 1 - : randomize the first spike so on average it occurs at - : start + noise*interval - event = start + invl(interval) - interval*(1. - noise) - : but not earlier than 0 - if (event < 0) { - event = 0 - } - net_send(event, 3) - } + seed(0) + on = 0 : off + ispike = 0 + if (noise < 0) { + noise = 0 + } + if (noise > 1) { + noise = 1 + } + if (start >= 0 && number > 0) { + on = 1 + : randomize the first spike so on average it occurs at + : start + noise*interval + event = start + invl(interval) - interval*(1. - noise) + : but not earlier than 0 + if (event < 0) { + event = 0 + } + net_send(event, 3) + } } PROCEDURE init_sequence(t(ms)) { - if (number > 0) { - on = 1 - event = 0 - ispike = 0 - } + if (number > 0) { + on = 1 + event = 0 + ispike = 0 + } } FUNCTION invl(mean (ms)) (ms) { - if (mean <= 0.) { - mean = .01 (ms) : I would worry if it were 0. - } - if (noise == 0) { - invl = mean - }else{ - invl = (1. - noise)*mean + noise*mean*erand() - } + if (mean <= 0.) { + mean = .01 (ms) : I would worry if it were 0. + } + if (noise == 0) { + invl = mean + } else { + invl = (1. - noise)*mean + noise*mean*erand() + } } -VERBATIM -#include "nrnran123.h" - -#if !NRNBBCORE -/* backward compatibility */ -double nrn_random_pick(void* r); -void* nrn_random_arg(int argpos); -int nrn_random_isran123(void* r, uint32_t* id1, uint32_t* id2, uint32_t* id3); -int nrn_random123_setseq(void* r, uint32_t seq, char which); -int nrn_random123_getseq(void* r, uint32_t* seq, char* which); -#endif -ENDVERBATIM FUNCTION erand() { -VERBATIM - if (_p_donotuse) { - /* - :Supports separate independent but reproducible streams for - : each instance. However, the corresponding hoc Random - : distribution MUST be set to Random.negexp(1) - */ -#if !NRNBBCORE - if (_ran_compat == 2) { - _lerand = nrnran123_negexp((nrnran123_State*)_p_donotuse); - }else{ - _lerand = nrn_random_pick(_p_donotuse); - } -#else - _lerand = nrnran123_negexp((nrnran123_State*)_p_donotuse); -#endif - return _lerand; - }else{ -#if NRNBBCORE - assert(0); -#else - /* - : the old standby. Cannot use if reproducible parallel sim - : independent of nhost or which host this instance is on - : is desired, since each instance on this cpu draws from - : the same stream - */ -#endif - } -#if !NRNBBCORE -ENDVERBATIM - erand = exprand(1) -VERBATIM -#endif -ENDVERBATIM -} - -PROCEDURE noiseFromRandom() { -VERBATIM -#if !NRNBBCORE - { - void** pv = (void**)(&_p_donotuse); - if (_ran_compat == 2) { - fprintf(stderr, "NetStim.noiseFromRandom123 was previously called\n"); - assert(0); - } - _ran_compat = 1; - if (ifarg(1)) { - *pv = nrn_random_arg(1); - }else{ - *pv = (void*)0; - } - } -#endif -ENDVERBATIM + erand = random_negexp(ranvar) } - -PROCEDURE noiseFromRandom123() { -VERBATIM -#if !NRNBBCORE - { - nrnran123_State** pv = (nrnran123_State**)(&_p_donotuse); - if (_ran_compat == 1) { - fprintf(stderr, "NetStim.noiseFromRandom was previously called\n"); - assert(0); - } - _ran_compat = 2; - if (*pv) { - nrnran123_deletestream(*pv); - *pv = (nrnran123_State*)0; - } - if (ifarg(3)) { - *pv = nrnran123_newstream3((uint32_t)*getarg(1), (uint32_t)*getarg(2), (uint32_t)*getarg(3)); - }else if (ifarg(2)) { - *pv = nrnran123_newstream((uint32_t)*getarg(1), (uint32_t)*getarg(2)); - } - } -#endif -ENDVERBATIM -} - -DESTRUCTOR { -VERBATIM - if (!noise) { return; } - if (_p_donotuse) { -#if NRNBBCORE - { /* but note that mod2c does not translate DESTRUCTOR */ -#else - if (_ran_compat == 2) { -#endif - nrnran123_State** pv = (nrnran123_State**)(&_p_donotuse); - nrnran123_deletestream(*pv); - *pv = (nrnran123_State*)0; - } - } -ENDVERBATIM +PROCEDURE next_invl() { + if (number > 0) { + event = invl(interval) + } + if (ispike >= number) { + on = 0 + } } -VERBATIM -static void bbcore_write(double* x, int* d, int* xx, int *offset, _threadargsproto_) { - if (!noise) { return; } - /* error if using the legacy scop_exprand */ - if (!_p_donotuse) { - fprintf(stderr, "NetStim: cannot use the legacy scop_negexp generator for the random stream.\n"); - assert(0); - } - if (d) { - char which; - uint32_t* di = ((uint32_t*)d) + *offset; -#if !NRNBBCORE - if (_ran_compat == 1) { - void** pv = (void**)(&_p_donotuse); - /* error if not using Random123 generator */ - if (!nrn_random_isran123(*pv, di, di+1, di+2)) { - fprintf(stderr, "NetStim: Random123 generator is required\n"); - assert(0); - } - nrn_random123_getseq(*pv, di+3, &which); - di[4] = (int)which; - }else{ -#else - { -#endif - nrnran123_State** pv = (nrnran123_State**)(&_p_donotuse); - nrnran123_getids3(*pv, di, di+1, di+2); - nrnran123_getseq(*pv, di+3, &which); - di[4] = (int)which; -#if NRNBBCORE - /* CORENeuron does not call DESTRUCTOR so... */ - nrnran123_deletestream(*pv); - *pv = (nrnran123_State*)0; -#endif - } - /*printf("Netstim bbcore_write %d %d %d\n", di[0], di[1], di[3]);*/ - } - *offset += 5; +NET_RECEIVE (w) { + if (flag == 0) { : external event + if (w > 0 && on == 0) { : turn on spike sequence + : but not if a netsend is on the queue + init_sequence(t) + : randomize the first spike so on average it occurs at + : noise*interval (most likely interval is always 0) + next_invl() + event = event - interval*(1. - noise) + net_send(event, 1) + }else if (w < 0) { : turn off spiking definitively + on = 0 + } + } + if (flag == 3) { : from INITIAL + if (on == 1) { : but ignore if turned off by external event + init_sequence(t) + net_send(0, 1) + } + } + if (flag == 1 && on == 1) { + ispike = ispike + 1 + net_event(t) + next_invl() + if (on == 1) { + net_send(event, 1) + } + } } -static void bbcore_read(double* x, int* d, int* xx, int* offset, _threadargsproto_) { - if (!noise) { return; } - /* Generally, CoreNEURON, in the context of psolve, begins with - an empty model so this call takes place in the context of a freshly - created instance and _p_donotuse is not NULL. - However, this function - is also now called from NEURON at the end of coreneuron psolve - in order to transfer back the nrnran123 sequence state. That - allows continuation with a subsequent psolve within NEURON or - properly transfer back to CoreNEURON if we continue the psolve - there. So now, extra logic is needed for this call to work in - a NEURON context. - */ +:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +: Legacy API +: +: Difference: seed(x) merely sets ranvar sequence to ((uint32_t)x, 0) +: noiseFromRandom HOC Random object must use Random123 +: generator. The ids and sequence are merely copied +: into ranvar. +:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - uint32_t* di = ((uint32_t*)d) + *offset; -#if NRNBBCORE - nrnran123_State** pv = (nrnran123_State**)(&_p_donotuse); - assert(!_p_donotuse); - *pv = nrnran123_newstream3(di[0], di[1], di[2]); - nrnran123_setseq(*pv, di[3], (char)di[4]); -#else - uint32_t id1, id2, id3; - assert(_p_donotuse); - if (_ran_compat == 1) { /* Hoc Random.Random123 */ - void** pv = (void**)(&_p_donotuse); - int b = nrn_random_isran123(*pv, &id1, &id2, &id3); - assert(b); - nrn_random123_setseq(*pv, di[3], (char)di[4]); - }else{ - assert(_ran_compat == 2); - nrnran123_State** pv = (nrnran123_State**)(&_p_donotuse); - nrnran123_getids3(*pv, &id1, &id2, &id3); - nrnran123_setseq(*pv, di[3], (char)di[4]); - } - /* Random123 on NEURON side has same ids as on CoreNEURON side */ - assert(di[0] == id1 && di[1] == id2 && di[2] == id3); -#endif - *offset += 5; -} -ENDVERBATIM +: the Random idiom has been extended to support CoreNEURON. -PROCEDURE next_invl() { - if (number > 0) { - event = invl(interval) - } - if (ispike >= number) { - on = 0 - } -} +: For backward compatibility, noiseFromRandom(hocRandom) can still be used +: as well as the default low-quality scop_exprand generator. +: However, CoreNEURON will not accept usage of the low-quality generator, +: and, if noiseFromRandom is used to specify the random stream, that stream +: must be using the Random123 generator. -NET_RECEIVE (w) { - if (flag == 0) { : external event - if (w > 0 && on == 0) { : turn on spike sequence - : but not if a netsend is on the queue - init_sequence(t) - : randomize the first spike so on average it occurs at - : noise*interval (most likely interval is always 0) - next_invl() - event = event - interval*(1. - noise) - net_send(event, 1) - }else if (w < 0) { : turn off spiking definitively - on = 0 - } - } - if (flag == 3) { : from INITIAL - if (on == 1) { : but ignore if turned off by external event - init_sequence(t) - net_send(0, 1) - } - } - if (flag == 1 && on == 1) { - ispike = ispike + 1 - net_event(t) - next_invl() - if (on == 1) { - net_send(event, 1) - } - } -} +: The recommended idiom for specfication of the random stream is to use +: noiseFromRandom123(id1, id2[, id3]) -FUNCTION bbsavestate() { - bbsavestate = 0 - : limited to noiseFromRandom123 -VERBATIM -#if !NRNBBCORE - if (_ran_compat == 2) { - nrnran123_State** pv = (nrnran123_State**)(&_p_donotuse); - if (!*pv) { return 0.0; } - char which; - uint32_t seq; - double *xdir, *xval; - xdir = hoc_pgetarg(1); - if (*xdir == -1.) { *xdir = 2; return 0.0; } - xval = hoc_pgetarg(2); - if (*xdir == 0.) { - nrnran123_getseq(*pv, &seq, &which); - xval[0] = (double)seq; - xval[1] = (double)which; - } - if (*xdir == 1) { - nrnran123_setseq(*pv, (uint32_t)xval[0], (char)xval[1]); - } - } /* else do nothing */ -#endif -ENDVERBATIM -} +: If any instance uses noiseFromRandom123, then no instance can use noiseFromRandom +: and vice versa. COMMENT @@ -409,3 +166,39 @@ its sequence. ENDCOMMENT +PROCEDURE seed(x) { + random_setseq(ranvar, x) +} + +PROCEDURE noiseFromRandom() { +VERBATIM +#if !NRNBBCORE + { + if (ifarg(1)) { + Rand* r = nrn_random_arg(1); + uint32_t id[3]; + if (!nrn_random_isran123(r, &id[0], &id[1], &id[2])) { + hoc_execerr_ext("NetStim: Random.Random123 generator is required."); + } + nrnran123_setids(ranvar, id[0], id[1], id[2]); + char which; + nrn_random123_getseq(r, &id[0], &which); + nrnran123_setseq(ranvar, id[0], which); + } + } +#endif +ENDVERBATIM +} + +PROCEDURE noiseFromRandom123() { +VERBATIM +#if !NRNBBCORE + if (ifarg(3)) { + nrnran123_setids(ranvar, static_cast(*getarg(1)), static_cast(*getarg(2)), static_cast(*getarg(3))); + } else if (ifarg(2)) { + nrnran123_setids(ranvar, static_cast(*getarg(1)), static_cast(*getarg(2)), 0); + } + nrnran123_setseq(ranvar, 0, 0); +#endif +ENDVERBATIM +} diff --git a/src/coreneuron/mechanism/register_mech.cpp b/src/coreneuron/mechanism/register_mech.cpp index 6e3741043d..9853c1243a 100644 --- a/src/coreneuron/mechanism/register_mech.cpp +++ b/src/coreneuron/mechanism/register_mech.cpp @@ -234,6 +234,8 @@ void hoc_register_dparam_semantics(int type, int ix, const char* name) { memb_func[type].dparam_semantics[ix] = -9; } else if (strcmp(name, "fornetcon") == 0) { memb_func[type].dparam_semantics[ix] = -10; + } else if (strcmp(name, "random") == 0) { + memb_func[type].dparam_semantics[ix] = -11; } else { int i = name[0] == '#' ? 1 : 0; int etype = nrn_get_mechtype(name + i); diff --git a/src/coreneuron/nrniv/nrniv_decl.h b/src/coreneuron/nrniv/nrniv_decl.h index 3cc20db667..00cdaef16d 100644 --- a/src/coreneuron/nrniv/nrniv_decl.h +++ b/src/coreneuron/nrniv/nrniv_decl.h @@ -48,6 +48,7 @@ extern void nrn_set_extra_thread0_vdata(void); extern Point_process* nrn_artcell_instantiate(const char* mechname); extern int nrnmpi_spike_compress(int nspike, bool gidcompress, int xchng); extern bool nrn_use_bin_queue_; +extern std::vector& nrn_mech_random_indices(int type); extern void nrn_outputevent(unsigned char, double); extern void ncs2nrn_integrate(double tstop); diff --git a/src/coreneuron/utils/nrnoc_aux.cpp b/src/coreneuron/utils/nrnoc_aux.cpp index 0ff43f3d2d..6d51dbe1f6 100644 --- a/src/coreneuron/utils/nrnoc_aux.cpp +++ b/src/coreneuron/utils/nrnoc_aux.cpp @@ -21,7 +21,7 @@ int v_structure_change; int diam_changed; #define MAXERRCOUNT 5 int hoc_errno_count; -const char* bbcore_write_version = "1.6"; // Allow multiple gid and PreSyn per real cell. +const char* bbcore_write_version = "1.7"; // NMODLRandom char* pnt_name(Point_process* pnt) { return corenrn.get_memb_func(pnt->_type).sym; diff --git a/src/coreneuron/utils/randoms/nrnran123.h b/src/coreneuron/utils/randoms/nrnran123.h index 57b30b2d93..efd4760691 100644 --- a/src/coreneuron/utils/randoms/nrnran123.h +++ b/src/coreneuron/utils/randoms/nrnran123.h @@ -134,6 +134,15 @@ inline double nrnran123_dblpick(nrnran123_State* s) { return nrnran123_uint2dbl(nrnran123_ipick(s)); } +// same as dblpick +inline double nrnran123_uniform(nrnran123_State* s) { + return nrnran123_uint2dbl(nrnran123_ipick(s)); +} + +inline double nrnran123_uniform(nrnran123_State* s, double low, double high) { + return low + nrnran123_uint2dbl(nrnran123_ipick(s)) * (high - low); +} + /* this could be called from openacc parallel construct (in INITIAL block) */ inline void nrnran123_setseq(nrnran123_State* s, uint32_t seq, char which) { if (which > 3) { @@ -145,6 +154,22 @@ inline void nrnran123_setseq(nrnran123_State* s, uint32_t seq, char which) { s->r = coreneuron_random123_philox4x32_helper(s); } +/* this could be called from openacc parallel construct (in INITIAL block) */ +inline void nrnran123_setseq(nrnran123_State* s, double seq34) { + if (seq34 < 0.0) { + seq34 = 0.0; + } + if (seq34 > double(0XffffffffffLL)) { + seq34 = 0.0; + } + + // at least 64 bits even on 32 bit machine (could be more) + unsigned long long x = ((unsigned long long) seq34) & 0X3ffffffffLL; + char which = x & 0X3; + uint32_t seq = x >> 2; + nrnran123_setseq(s, seq, which); +} + // nrnran123_negexp min value is 2.3283064e-10, max is 22.18071, mean 1.0 inline double nrnran123_negexp(nrnran123_State* s) { return -std::log(nrnran123_dblpick(s)); diff --git a/src/gnu/nrnran123.cpp b/src/gnu/nrnran123.cpp index eeb2768115..6180140c09 100644 --- a/src/gnu/nrnran123.cpp +++ b/src/gnu/nrnran123.cpp @@ -23,10 +23,18 @@ std::uint32_t nrnran123_get_globalindex() { return k[0]; } -nrnran123_State* nrnran123_newstream(std::uint32_t id1, std::uint32_t id2) { - return nrnran123_newstream3(id1, id2, 0); -} +/* deprecated */ nrnran123_State* nrnran123_newstream3(std::uint32_t id1, std::uint32_t id2, std::uint32_t id3) { + return nrnran123_newstream(id1, id2, id3); +} + +nrnran123_State* nrnran123_newstream() { + extern int nrnmpi_myid; + static std::uint32_t id3{}; + return nrnran123_newstream(1, nrnmpi_myid, ++id3); +} + +nrnran123_State* nrnran123_newstream(std::uint32_t id1, std::uint32_t id2, std::uint32_t id3) { auto* s = new nrnran123_State; s->c[1] = id3; s->c[2] = id1; @@ -54,12 +62,29 @@ void nrnran123_setseq(nrnran123_State* s, std::uint32_t seq, char which) { s->r = philox4x32(s->c, k); } +/** @brief seq4which is 34 bit uint encoded as double(seq)*4 + which + * More convenient to get and set from interpreter +*/ +void nrnran123_setseq(nrnran123_State* s, double seq4which) { + if (seq4which < 0.0) { + seq4which = 0.0; + } + if (seq4which > double(0XffffffffffLL)) { + seq4which = 0.0; + } + // at least 64 bits even on 32 bit machine (could be more) + unsigned long long x = ((unsigned long long) seq4which) & 0X3ffffffffLL; + char which = x & 0X3; + uint32_t seq = x >> 2; + nrnran123_setseq(s, seq, which); +} + void nrnran123_getids(nrnran123_State* s, std::uint32_t* id1, std::uint32_t* id2) { *id1 = s->c[2]; *id2 = s->c[3]; } -void nrnran123_getids3(nrnran123_State* s, +void nrnran123_getids(nrnran123_State* s, std::uint32_t* id1, std::uint32_t* id2, std::uint32_t* id3) { @@ -68,6 +93,20 @@ void nrnran123_getids3(nrnran123_State* s, *id2 = s->c[3]; } +/* Deprecated */ +void nrnran123_getids3(nrnran123_State* s, + std::uint32_t* id1, + std::uint32_t* id2, + std::uint32_t* id3) { + nrnran123_getids(s, id1, id2, id3); +} + +void nrnran123_setids(nrnran123_State* s, std::uint32_t id1, std::uint32_t id2, std::uint32_t id3) { + s->c[1] = id3; + s->c[2] = id1; + s->c[3] = id2; +} + std::uint32_t nrnran123_ipick(nrnran123_State* s) { char which = s->which_; std::uint32_t rval = s->r[which++]; @@ -86,6 +125,19 @@ double nrnran123_dblpick(nrnran123_State* s) { return ((double) u + 1.0) * SHIFT32; } +double nrnran123_uniform(nrnran123_State* s) { + return nrnran123_dblpick(s); +} + +double nrnran123_uniform(nrnran123_State* s, double a, double b) { + return a + nrnran123_dblpick(s) * (b - a); +} + +double nrnran123_negexp(nrnran123_State* s, double mean) { + /* min 2.3283064e-10 to max 22.18071 (if mean is 1) */ + return -std::log(nrnran123_dblpick(s)) * mean; +} + double nrnran123_negexp(nrnran123_State* s) { /* min 2.3283064e-10 to max 22.18071 */ return -std::log(nrnran123_dblpick(s)); @@ -104,11 +156,15 @@ double nrnran123_normal(nrnran123_State* s) { w = (u1 * u1) + (u2 * u2); } while (w > 1); - y = std::sqrt((-2. * log(w)) / w); + y = std::sqrt((-2. * std::log(w)) / w); x = u1 * y; return x; } +double nrnran123_normal(nrnran123_State* s, double mu, double sigma) { + return mu + nrnran123_normal(s) * sigma; +} + nrnran123_array4x32 nrnran123_iran(std::uint32_t seq, std::uint32_t id1, std::uint32_t id2) { return nrnran123_iran3(seq, id1, id2, 0); } diff --git a/src/gnu/nrnran123.h b/src/gnu/nrnran123.h index a4e4844bf6..e087a1607e 100644 --- a/src/gnu/nrnran123.h +++ b/src/gnu/nrnran123.h @@ -37,32 +37,91 @@ struct nrnran123_array4x32 { extern void nrnran123_set_globalindex(std::uint32_t gix); extern std::uint32_t nrnran123_get_globalindex(); -/* minimal data stream */ -extern nrnran123_State* nrnran123_newstream(std::uint32_t id1, std::uint32_t id2); +/** Construct a new Random123 stream based on the philox4x32 generator. + * + * @param id1 stream ID + * @param id2 optional defaults to 0 + * @param id3 optional defaults to 0 + * @return an nrnran123_State object representing this stream + */ +extern nrnran123_State* nrnran123_newstream(std::uint32_t id1, + std::uint32_t id2 = 0, + std::uint32_t id3 = 0); + +/** @deprecated use nrnran123_newstream instead **/ extern nrnran123_State* nrnran123_newstream3(std::uint32_t id1, std::uint32_t id2, std::uint32_t id3); -extern void nrnran123_deletestream(nrnran123_State*); -extern void nrnran123_getseq(nrnran123_State*, std::uint32_t* seq, char* which); -extern void nrnran123_setseq(nrnran123_State*, std::uint32_t seq, char which); -extern void nrnran123_getids(nrnran123_State*, std::uint32_t* id1, std::uint32_t* id2); + +/** Construct a new Random123 stream based on the philox4x32 generator. + * + * @note This overload constructs each stream instance as an independent stream. + * Independence is derived by using id1=1, id2=nrnmpi_myid, + * id3 = ++internal_static_uint32_initialized_to_0 + */ +extern nrnran123_State* nrnran123_newstream(); + +/** Destroys the given Random123 stream. */ +extern void nrnran123_deletestream(nrnran123_State* s); + +/** Get sequence number and selector from an nrnran123_State object */ +extern void nrnran123_getseq(nrnran123_State* s, std::uint32_t* seq, char* which); + +/** Set a Random123 sequence for a sequnece ID and which selector. + * + * @param s an Random123 state object + * @param seq the sequence ID for which to initialize the random number sequence + * @param which the selector (0 <= which < 4) of the sequence + */ +extern void nrnran123_setseq(nrnran123_State* s, std::uint32_t seq, char which); + +/** Set a Random123 sequence for a sequnece ID and which selector. + * + * This overload encodes the sequence ID and which in one double. This is done specifically to be + * able to expose the Random123 API in HOC, which only supports real numbers. + * + * @param s an Random123 state object + * @param seq4which encodes both seq and which as seq*4+which + */ +extern void nrnran123_setseq(nrnran123_State* s, double seq4which); // seq*4+which); + +/** Get stream IDs from Random123 State object */ +extern void nrnran123_getids(nrnran123_State* s, std::uint32_t* id1, std::uint32_t* id2); + + +/** Get stream IDs from Random123 State object */ +extern void nrnran123_getids(nrnran123_State* s, + std::uint32_t* id1, + std::uint32_t* id2, + std::uint32_t* id3); + +/** @brief. Deprecated, use nrnran123_getids **/ extern void nrnran123_getids3(nrnran123_State*, std::uint32_t* id1, std::uint32_t* id2, std::uint32_t* id3); +extern void nrnran123_setids(nrnran123_State*, std::uint32_t id1, std::uint32_t id2, std::uint32_t id3); + // Get a random uint32_t in [0, 2^32-1] extern std::uint32_t nrnran123_ipick(nrnran123_State*); // Get a random double on [0, 1] // nrnran123_dblpick minimum value is 2.3283064e-10 and max value is 1-min extern double nrnran123_dblpick(nrnran123_State*); -/* nrnran123_negexp min value is 2.3283064e-10, max is 22.18071 */ -extern double nrnran123_negexp(nrnran123_State*); /* mean 1.0 */ -extern double nrnran123_normal(nrnran123_State*); /* mean 0.0, std 1.0 */ +/* nrnran123_negexp min value is 2.3283064e-10, max is 22.18071 if mean is 1.0 */ +extern double nrnran123_negexp(nrnran123_State*); // mean = 1.0 +extern double nrnran123_negexp(nrnran123_State*, double mean); +extern double nrnran123_normal(nrnran123_State*); // mean = 0.0, std = 1.0 +extern double nrnran123_normal(nrnran123_State*, double mean, double std); + +extern double nrnran123_uniform(nrnran123_State*); // same as dblpick +extern double nrnran123_uniform(nrnran123_State*, double min, double max); /* more fundamental (stateless) (though the global index is still used) */ -extern nrnran123_array4x32 nrnran123_iran(std::uint32_t seq, std::uint32_t id1, std::uint32_t id2); +extern nrnran123_array4x32 nrnran123_iran(std::uint32_t seq, std::uint32_t id1, std::uint32_t id2 = 0, std::uint32_t id3 = 0); + +/** @brief. Deprecated, use nrnran123_iran **/ extern nrnran123_array4x32 nrnran123_iran3(std::uint32_t seq, std::uint32_t id1, std::uint32_t id2, diff --git a/src/modlunit/extdef.h b/src/modlunit/extdef.h index aa7c8a51e7..36da115736 100644 --- a/src/modlunit/extdef.h +++ b/src/modlunit/extdef.h @@ -8,4 +8,5 @@ "gauss", "normrand", "poisrand", "poisson", "setseed", "scop_random", "boundary", "romberg", "legendre", "invert", "stepforce", "schedule", "set_seed", "nrn_pointing", "state_discontinuity", "net_send", "net_move", "net_event", "nrn_random_play", "at_time", - "nrn_ghk", + "nrn_ghk", "random_negexp", "random_normal", "random_uniform", "random_setseq", "random_setids", + "random_ipick", diff --git a/src/modlunit/init.cpp b/src/modlunit/init.cpp index 9e92331a97..80250d47ab 100644 --- a/src/modlunit/init.cpp +++ b/src/modlunit/init.cpp @@ -79,6 +79,7 @@ static struct { /* Keywords */ {"READ", READ}, {"WRITE", WRITE}, {"RANGE", RANGE}, + {"RANDOM", RANDOM}, {"VALENCE", VALENCE}, {"CHARGE", VALENCE}, {"GLOBAL", GLOBAL}, diff --git a/src/modlunit/model.h b/src/modlunit/model.h index 1118488599..2433436c53 100644 --- a/src/modlunit/model.h +++ b/src/modlunit/model.h @@ -128,6 +128,7 @@ extern List* _LST(Item* q, char* file, int line); #define LOCL 0400000L #define CNVFAC 01000000L #define UFACTOR 02000000L +#define RANGEOBJ 04000000L #define EXPLICIT_DECL 01 /* usage field, variable occurs in input file */ diff --git a/src/modlunit/nrnunit.cpp b/src/modlunit/nrnunit.cpp index 6af59eb7fd..fa71749b68 100644 --- a/src/modlunit/nrnunit.cpp +++ b/src/modlunit/nrnunit.cpp @@ -85,6 +85,12 @@ void nrn_list(Item* qtype, Item* qlist) { point_process = 1; } break; + case RANDOM: + plist = (List**) 0; + ITERATE(q, qlist) { + declare(RANGEOBJ, q, nullptr); + } + break; default: plist = (List**) 0; break; diff --git a/src/modlunit/parse1.ypp b/src/modlunit/parse1.ypp index af77480730..6e0e4778b3 100755 --- a/src/modlunit/parse1.ypp +++ b/src/modlunit/parse1.ypp @@ -102,6 +102,7 @@ extern int lexcontext; %type initstmt bablk %token CONDUCTANCE %type conducthint +%token RANDOM /* precedence in expressions--- low to high */ %left OR @@ -752,6 +753,8 @@ nrnstmt: /*nothing*/ { P1{nrn_list($2,$3);}} | nrnstmt RANGE nrnlist { P1{nrn_list($2, $3);}} + | nrnstmt RANDOM nrnlist + { P1{nrn_list($2, $3);}} | nrnstmt GLOBAL nrnlist { P1{nrn_list($2, $3);}} | nrnstmt POINTER nrnlist diff --git a/src/nmodl/init.cpp b/src/nmodl/init.cpp index 35b0cf6715..b5c89c3c33 100644 --- a/src/nmodl/init.cpp +++ b/src/nmodl/init.cpp @@ -91,6 +91,7 @@ static struct { /* Keywords */ {"MUTEXLOCK", NRNMUTEXLOCK}, {"MUTEXUNLOCK", NRNMUTEXUNLOCK}, {"REPRESENTS", REPRESENTS}, + {"RANDOM", RANDOM}, {0, 0}}; /* @@ -153,6 +154,17 @@ static const char* extdef5[] = {/* the extdef names that are not threadsafe */ #include "extdef5.h" 0}; +/* random to nrnran123 functions */ +std::map extdef_rand = { + {"random_setseq", "nrnran123_setseq"}, + {"random_setids", "nrnran123_setids"}, + {"random_uniform", "nrnran123_uniform"}, + {"random_negexp", "nrnran123_negexp"}, + {"random_normal", "nrnran123_normal"}, + {"random_ipick", "nrnran123_ipick"}, + {"random_dpick", "nrnran123_dblpick"}, +}; + List *constructorfunc, *destructorfunc; void init() { @@ -197,6 +209,10 @@ void init() { assert(s); s->subtype |= EXTDEF5; } + for (auto it: extdef_rand) { + s = install(it.first.c_str(), NAME); + s->subtype = EXTDEF_RANDOM; + } intoken = newlist(); initfunc = newlist(); modelfunc = newlist(); diff --git a/src/nmodl/modl.h b/src/nmodl/modl.h index 858134e57a..0da34de44f 100644 --- a/src/nmodl/modl.h +++ b/src/nmodl/modl.h @@ -4,6 +4,8 @@ #include #include #include +#include +#include /** * \dir @@ -155,6 +157,8 @@ typedef struct Symbol { } Symbol; #define SYM0 (Symbol*) 0 +extern std::map extdef_rand; + /* * this is convenient way to get the element pointer if you know what type * the item is @@ -195,13 +199,14 @@ typedef struct Symbol { #define EXTDEF 0100000 #define LINF 0200000 #define UNITDEF 0400000L -#define EXTDEF2 01000000L /* functions that can take array or function name arguments */ -#define nmodlCONST 02000000L /* constants that do not appear in .var file */ -#define EXTDEF3 04000000L /* get two extra reset arguments at beginning */ -#define INTGER 010000000L /* must be cast to double in expr */ -#define EXTDEF4 020000000L /* get extra NrnThread* arg at beginning */ -#define EXTDEF5 040000000L /* not threadsafe from the extdef list */ -#define EXPLICIT_DECL 01 /* usage field, variable occurs in input file */ +#define EXTDEF2 01000000L /* functions that can take array or function name arguments */ +#define nmodlCONST 02000000L /* constants that do not appear in .var file */ +#define EXTDEF3 04000000L /* get two extra reset arguments at beginning */ +#define INTGER 010000000L /* must be cast to double in expr */ +#define EXTDEF4 020000000L /* get extra NrnThread* arg at beginning */ +#define EXTDEF5 040000000L /* not threadsafe from the extdef list */ +#define EXTDEF_RANDOM 0600000000L /* functions that can be used with RANDOM type */ +#define EXPLICIT_DECL 01 /* usage field, variable occurs in input file */ #define NRNEXTRN 01 /* t, dt, celsius, etc. */ @@ -219,6 +224,7 @@ typedef struct Symbol { #define NRNPOINTER 04000 #define IONCONC 010000 #define NRNBBCOREPOINTER 020000 +#define NMODLRANDOM 040000 // Implicit ion concentration variable that has been added so we can call nrn_wrote_conc, but which // is not used in the MOD file #define IONCONC_IMPLICIT 040000 @@ -325,6 +331,5 @@ extern Item* qlint; #endif using neuron::Sprintf; - void verbatim_adjust(char* q); /** @} */ // end of hoc_functions diff --git a/src/nmodl/nocpout.cpp b/src/nmodl/nocpout.cpp index 743d224166..67a2d85672 100644 --- a/src/nmodl/nocpout.cpp +++ b/src/nmodl/nocpout.cpp @@ -1,4 +1,5 @@ #include <../../nmodlconf.h> + /* /local/src/master/nrn/src/nmodl/nocpout.c,v 4.1 1997/08/30 20:45:28 hines Exp */ /* @@ -125,6 +126,9 @@ static List* rangeparm; static List* rangedep; static List* rangestate; static List* nrnpointers; +static List* nmodlrandoms; +static List* nrn_mech_inst_destruct_list; +static int num_random_vars = 0; static char suffix[256]; static const char* rsuffix; /* point process range and functions don't have suffix*/ static char* mechname; @@ -193,6 +197,7 @@ void nrninit() { debugging_ = 1; thread_cleanup_list = newlist(); thread_mem_init_list = newlist(); + nmodlrandoms = newlist(); } void parout() { @@ -759,7 +764,11 @@ extern void nrn_promote(Prop*, int, int);\n\ "Memb_list*, int);\n"); } /* count the number of pointers needed */ - ppvar_cnt = ioncount + diamdec + pointercount + areadec; + num_random_vars = 0; + ITERATE(q, nmodlrandoms) { + num_random_vars++; + } + ppvar_cnt = ioncount + diamdec + pointercount + num_random_vars + areadec; if (net_send_seen_) { tqitem_index = ppvar_cnt; ppvar_semantics( @@ -901,6 +910,8 @@ static const char *_mechanism[] = {\n\ q = q->next->next->next; } + Item* before_nrn_alloc = lappendstr(defs_list, "\n"); + Lappendstr(defs_list, "\n" "extern Prop* need_memb(Symbol*);\n" @@ -1047,6 +1058,30 @@ static const char *_mechanism[] = {\n\ } } + + // I've put all the nrn_mech_inst_destruct here with nmodlrandoms allocation. + // Refactor if ever things other than nmodlrandoms need it. + nrn_mech_inst_destruct_list = newlist(); + ITERATE(q, nmodlrandoms) { + Sprintf(buf, "_p_%s = (void*)nrnran123_newstream();\n", SYM(q)->name); + Lappendstr(defs_list, buf); + Sprintf(buf, "nrnran123_deletestream(%s);\n", SYM(q)->name); + Lappendstr(nrn_mech_inst_destruct_list, buf); + } + if (nrn_mech_inst_destruct_list != nrn_mech_inst_destruct_list->next) { + auto& list = nrn_mech_inst_destruct_list; + // registration just means adding to nrn_mech_inst_destruct + Lappendstr(defs_list, "nrn_mech_inst_destruct[_mechtype] = _mech_inst_destruct;\n"); + // boilerplate for _mech_inst_destruct + Linsertstr(list, + "\nstatic void _mech_inst_destruct(Prop* _prop) {\n" + " Datum* _ppvar = _nrn_mechanism_access_dparam(_prop);\n"); + Lappendstr(list, "}\n"); + movelist(list->next, list->prev, procfunc); + // need a forward declaration before nrn_alloc. + insertstr(before_nrn_alloc, "\nstatic void _mech_inst_destruct(Prop* _prop);\n"); + } + if (constructorfunc->next != constructorfunc) { Lappendstr(defs_list, "if (!nrn_point_prop_) {_constructor(_prop);}\n"); if (vectorize) { @@ -1840,6 +1875,17 @@ void nrn_list(Item* q1, Item* q2) { } use_bbcorepointer = 1; break; + case RANDOM: + for (q = q1->next; q != q2->next; q = q->next) { + Symbol* s = SYM(q); + if (s->type != NAME || s->subtype || s->nrntype) { + diag(s->name, " cannot be redeclared as RANDOM"); + } + s->nrntype |= NRNNOTP | EXTDEF_RANDOM; + s->type = RANDOMVAR; + } + plist = &nmodlrandoms; + break; } if (plist) { if (!*plist) { @@ -2397,8 +2443,38 @@ int iondef(int* p_pointercount) { (*p_pointercount)++; } + // print all RANDOM variables + num_random_vars = 0; + ITERATE(q, nmodlrandoms) { + num_random_vars++; + } + if (num_random_vars) { + Sprintf(buf, "\n //RANDOM variables \n"); + lappendstr(defs_list, buf); + + int index = 0; + ITERATE(q, nmodlrandoms) { + Symbol* s = SYM(q); + Sprintf(buf, + "#define %s (nrnran123_State*)_ppvar[%d].get()\n", + s->name, + ioncount + *p_pointercount + index); + lappendstr(defs_list, buf); + Sprintf(buf, + "#define _p_%s _ppvar[%d].literal_value()\n", + s->name, + ioncount + *p_pointercount + index); + lappendstr(defs_list, buf); + ppvar_semantics(ioncount + *p_pointercount + index, "random", s->name, "void*"); + index++; + } + lappendstr(defs_list, "\n"); + } + if (diamdec) { /* must be last */ - Sprintf(buf, "#define diam *_ppvar[%d].get()\n", ioncount + *p_pointercount); + Sprintf(buf, + "#define diam *_ppvar[%d].get()\n", + ioncount + *p_pointercount + num_random_vars); q2 = lappendstr(defs_list, buf); q2->itemtype = VERBATIM; } /* notice that ioncount is not incremented */ @@ -2406,13 +2482,13 @@ int iondef(int* p_pointercount) { procedures must be redone */ Sprintf(buf, "#define area *_ppvar[%d].get()\n", - ioncount + *p_pointercount + diamdec); + ioncount + *p_pointercount + num_random_vars + diamdec); q2 = lappendstr(defs_list, buf); q2->itemtype = VERBATIM; } /* notice that ioncount is not incremented */ Sprintf(buf, "static constexpr auto number_of_datum_variables = %d;\n", - ioncount + *p_pointercount + diamdec + areadec); + ioncount + *p_pointercount + num_random_vars + diamdec + areadec); linsertstr(defs_list, buf)->itemtype = VERBATIM; return ioncount; } diff --git a/src/nmodl/parsact.cpp b/src/nmodl/parsact.cpp index 1bfcea0664..79e8b96584 100644 --- a/src/nmodl/parsact.cpp +++ b/src/nmodl/parsact.cpp @@ -222,7 +222,7 @@ void vectorize_scan_for_func(Item* q1, Item* q2) { for (q = q1; q != q2; q = q->next) { if (q->itemtype == SYMBOL) { Symbol* s = SYM(q); - if ((s->usage & FUNCT) && !(s->subtype & (EXTDEF))) { + if ((s->usage & FUNCT) && !(s->subtype & (EXTDEF | EXTDEF_RANDOM))) { if (q->next->itemtype == SYMBOL && strcmp(SYM(q->next)->name, "(") == 0) { int b = func_arg_examine(q->next, q2); if (b == 0) { /* no args */ @@ -845,7 +845,7 @@ void hocfunc(Symbol* n, Item* qpar1, Item* qpar2) /*interface between modl and h /* ARGSUSED */ void vectorize_use_func(Item* qname, Item* qpar1, Item* qexpr, Item* qpar2, int blocktype) { Item* q; - if (SYM(qname)->subtype & EXTDEF) { + if (SYM(qname)->subtype & (EXTDEF | EXTDEF_RANDOM)) { if (strcmp(SYM(qname)->name, "nrn_pointing") == 0) { // TODO: this relies on undefined behaviour in C++. &*foo is not // guaranteed to be equivalent to foo if foo is null. See @@ -914,6 +914,8 @@ void vectorize_use_func(Item* qname, Item* qpar1, Item* qexpr, Item* qpar2, int } else { diag("net_move", "only allowed in NET_RECEIVE block"); } + } else if (SYM(qname)->subtype & EXTDEF_RANDOM) { + replacstr(qname, extdef_rand[SYM(qname)->name]); } return; } diff --git a/src/nmodl/parse1.ypp b/src/nmodl/parse1.ypp index a697d50578..03c7108c46 100755 --- a/src/nmodl/parse1.ypp +++ b/src/nmodl/parse1.ypp @@ -111,6 +111,7 @@ static int nr_argcnt_, argcnt_; /* for matching number of args in NET_RECEIVE %type neuronblk nrnuse nrnlist valence initstmt bablk optontology %token CONDUCTANCE %type conducthint +%token RANDOM RANDOMVAR /* precedence in expressions--- low to high */ %left OR @@ -186,6 +187,11 @@ Name: NAME SYM($1) = checklocal(SYM($1)); /* it was a bug when this was done to the lookahead token in lex */ } + | RANDOMVAR error { + std::string s{SYM($1)->name}; + s += " RANDOM var can only be used as the first arg of a random_ function"; + myerr(s.c_str()); + } ; declare: parmblk | indepblk | depblk | stateblk | neuronblk | unitblk | constblk @@ -622,10 +628,21 @@ fprintf(stderr, "Notice: %s is not thread safe\n", SYM($1)->name); vectorize = 0; } } + Item* arg = $2->next; + if (SYM($1)->subtype & EXTDEF_RANDOM) { + if (arg == $2 || arg->itemtype != SYMBOL || SYM(arg)->type != RANDOMVAR) { + diag(SYM($1)->name, " must have RANDOM var as first argument"); + } + }else{ + if (arg != $2 && arg->itemtype == SYMBOL && SYM(arg)->type == RANDOMVAR) { + diag(SYM($1)->name, " cannot have RANDOM var as an argument"); + } + } vectorize_use_func($1,$2,$4,$5,blocktype); } ; exprlist: /*nothing*/{$$ = ITEM0;} + | RANDOMVAR | expr | STRING | exprlist ',' expr @@ -1010,6 +1027,8 @@ nrnstmt: /*nothing*/ { nrn_list($2, $3);} | nrnstmt THREADSAFE { assert_threadsafe = 1; } + | nrnstmt RANDOM nrnlist + { nrn_list($2, $3);} ; nrnuse: USEION NAME READ nrnlist valence optontology {nrn_use($2, $4, ITEM0, $5);} diff --git a/src/nrniv/bbsavestate.cpp b/src/nrniv/bbsavestate.cpp index 3fd6cc30c0..7633c4a1b5 100644 --- a/src/nrniv/bbsavestate.cpp +++ b/src/nrniv/bbsavestate.cpp @@ -172,6 +172,7 @@ callback to bbss_early when needed. #include "ndatclas.h" #include "nrncvode.h" #include "nrnoc2iv.h" +#include "nrnran123.h" #include "ocfile.h" #include #include @@ -2056,12 +2057,40 @@ void BBSaveState::mech(Prop* p) { f->s(buf, 1); { auto const size = ssi[p->_type].size; // sum over array dimensions for range variables + auto& random_indices = nrn_mech_random_indices(p->_type); + auto size_random = random_indices.size(); std::vector tmp{}; - tmp.reserve(size); + tmp.reserve(size + size_random); for (auto i = 0; i < size; ++i) { tmp.push_back(static_cast(p->param_handle_legacy(ssi[p->_type].offset + i))); } - f->d(size, tmp.data()); + + // read or write the RANDOM 34 sequence values by pointing last + // size_random tmp elements to seq34 double slots. + std::vector seq34(size_random, 0); + for (auto i = 0; i < size_random; ++i) { + tmp.push_back(static_cast(&seq34[i])); + } + // if writing, nrnran123_getseq into seq34 + if (f->type() == BBSS_IO::OUT) { // save + for (auto i = 0; i < size_random; ++i) { + uint32_t seq{}; + char which{}; + auto& datum = p->dparam[random_indices[i]]; + nrnran123_State* n123s = (nrnran123_State*) datum.get(); + nrnran123_getseq(n123s, &seq, &which); + seq34[i] = 4.0 * double(seq) + double(which); + } + } + f->d(size + size_random, tmp.data()); + // if reading, seq34 into nrnran123_setseq + if (f->type() == BBSS_IO::IN) { // restore + for (auto i = 0; i < size_random; ++i) { + auto& datum = p->dparam[random_indices[i]]; + nrnran123_State* n123s = (nrnran123_State*) datum.get(); + nrnran123_setseq(n123s, seq34[i]); + } + } } Point_process* pp{}; if (memb_func[p->_type].is_point) { diff --git a/src/nrniv/ndatclas.cpp b/src/nrniv/ndatclas.cpp index c6d4defdd6..d336e459ab 100644 --- a/src/nrniv/ndatclas.cpp +++ b/src/nrniv/ndatclas.cpp @@ -217,7 +217,7 @@ Symbol* NrnProperty::find(const char* name) { } int NrnProperty::prop_index(const Symbol* s) const { assert(s); - if (s->type != RANGEVAR) { + if (s->type != RANGEVAR && s->type != RANGEOBJ) { hoc_execerror(s->name, "not a range variable"); } return s->u.rng.index; diff --git a/src/nrniv/nmodlrandom.cpp b/src/nrniv/nmodlrandom.cpp new file mode 100644 index 0000000000..8432c64e78 --- /dev/null +++ b/src/nrniv/nmodlrandom.cpp @@ -0,0 +1,141 @@ +#include <../../nrnconf.h> + +/* +HOC wrapper object for NMODL NEURON block RANDOM variables to give a HOC +pointprocess.ranvar.method(...) +sec.ranvar_mech(x).method(...) +and Python +poinprocess.ranvar.method(...) +sec(x).mech.ranvar.method(...) +syntax +*/ + +#include +#include +#include +#include +#include + +struct NMODLRandom { + NMODLRandom(Object*) {} + ~NMODLRandom() {} + nrnran123_State* r() { + return (nrnran123_State*) hr_.get(); + } + void chk() { + if (!prop_id_) { + hoc_execerr_ext("NMODLRandom wrapped handle is not valid"); + } + } + neuron::container::generic_data_handle hr_{}; + neuron::container::non_owning_identifier_without_container prop_id_{}; +}; + +static Symbol* nmodlrandom_sym{}; +#undef dmaxuint +#define dmaxuint 4294967295. + +static Object** set_ids(void* v) { // return this NMODLRandom instance + NMODLRandom* r = (NMODLRandom*) v; + r->chk(); + uint32_t id[3]; + for (int i = 0; i < 3; ++i) { + id[i] = (uint32_t) (chkarg(i + 1, 0., dmaxuint)); + } + nrnran123_setids(r->r(), id[0], id[1], id[2]); + return hoc_temp_objptr(nrn_get_gui_redirect_obj()); +} + +static Object** get_ids(void* v) { // return a Vector of size 3. + NMODLRandom* r = (NMODLRandom*) v; + r->chk(); + uint32_t id[3]{}; + nrnran123_getids3(r->r(), id, id + 1, id + 2); + IvocVect* vec = vector_new1(3); + double* data = vector_vec(vec); + for (int i = 0; i < 3; ++i) { + data[i] = double(id[i]); + } + return vector_temp_objvar(vec); +} + +static Object** set_seq(void* v) { // return this NModlRandom instance + NMODLRandom* r = (NMODLRandom*) v; + r->chk(); + double seq = *getarg(1); + nrnran123_setseq(r->r(), seq); + return hoc_temp_objptr(nrn_get_gui_redirect_obj()); +} + +static double get_seq(void* v) { // return the 34 bits (seq*4 + which) as double + NMODLRandom* r = (NMODLRandom*) v; + r->chk(); + uint32_t seq; + char which; + nrnran123_getseq(r->r(), &seq, &which); + return double(seq) * 4.0 + which; +} + +static double uniform(void* v) { + NMODLRandom* r = (NMODLRandom*) v; + r->chk(); + return nrnran123_uniform(r->r()); +} + +static Member_func members[] = {{"get_seq", get_seq}, {"uniform", uniform}, {nullptr, nullptr}}; + +static Member_ret_obj_func retobj_members[] = {{"set_ids", set_ids}, + {"get_ids", get_ids}, + {"set_seq", set_seq}, + {nullptr, nullptr}}; + +static void* nmodlrandom_cons(Object*) { + NMODLRandom* r = new NMODLRandom(nullptr); + return r; +} + +static void nmodlrandom_destruct(void* v) { + NMODLRandom* r = (NMODLRandom*) v; + delete r; +} + +void NMODLRandom_reg() { + class2oc("NMODLRandom", + nmodlrandom_cons, + nmodlrandom_destruct, + members, + nullptr, + retobj_members, + nullptr); + if (!nmodlrandom_sym) { + nmodlrandom_sym = hoc_lookup("NMODLRandom"); + assert(nmodlrandom_sym); + } +} + +Object* nrn_pntproc_nmodlrandom_wrap(void* v, Symbol* sym) { + auto* const pnt = static_cast(v); + if (!pnt->prop) { + if (nrn_inpython_ == 1) { /* python will handle the error */ + hoc_warning("point process not located in a section", nullptr); + nrn_inpython_ = 2; + return {}; + } else { + hoc_execerror("point process not located in a section", nullptr); + } + } + + return nrn_nmodlrandom_wrap(pnt->prop, sym); +} + +Object* nrn_nmodlrandom_wrap(Prop* prop, Symbol* sym) { + assert(sym->type == RANGEOBJ && sym->subtype == NMODLRANDOM); + auto& datum = prop->dparam[sym->u.rng.index]; + assert(datum.holds()); + + NMODLRandom* r = new NMODLRandom(nullptr); + r->hr_ = datum; + r->prop_id_ = prop->id(); + Object* wrap = hoc_new_object(nmodlrandom_sym, r); + return wrap; +} diff --git a/src/nrniv/nrnclass.h b/src/nrniv/nrnclass.h index 73e89cec23..02659d7993 100644 --- a/src/nrniv/nrnclass.h +++ b/src/nrniv/nrnclass.h @@ -3,9 +3,9 @@ , Shape_reg(), PlotShape_reg(), PPShape_reg(), RangeVarPlot_reg(), SectionBrowser_reg(), MechanismStandard_reg(), MechanismType_reg(), NetCon_reg(), LinearMechanism_reg(), KSChan_reg(), Impedance_reg(), SaveState_reg(), BBSaveState_reg(), FInitializeHandler_reg(), - StateTransitionEvent_reg(), nrnpython_reg() + StateTransitionEvent_reg(), nrnpython_reg(), NMODLRandom_reg() #if USEDASPK - , + , Daspk_reg() #endif #if USECVODE @@ -26,7 +26,7 @@ , Shape_reg, PlotShape_reg, PPShape_reg, RangeVarPlot_reg, SectionBrowser_reg, MechanismStandard_reg, MechanismType_reg, NetCon_reg, LinearMechanism_reg, KSChan_reg, Impedance_reg, SaveState_reg, BBSaveState_reg, FInitializeHandler_reg, StateTransitionEvent_reg, - nrnpython_reg + nrnpython_reg, NMODLRandom_reg #if USEDASPK , Daspk_reg diff --git a/src/nrniv/nrncore_write/callbacks/nrncore_callbacks.cpp b/src/nrniv/nrncore_write/callbacks/nrncore_callbacks.cpp index 49a1bf6806..ec88e782a7 100644 --- a/src/nrniv/nrncore_write/callbacks/nrncore_callbacks.cpp +++ b/src/nrniv/nrncore_write/callbacks/nrncore_callbacks.cpp @@ -258,7 +258,7 @@ int nrnthread_dat2_1(int tid, cg.ml_vdata_offset[j] = vdata_offset; int* ds = memb_func[type].dparam_semantics; for (int psz = 0; psz < bbcore_dparam_size[type]; ++psz) { - if (ds[psz] == -4 || ds[psz] == -6 || ds[psz] == -7 || ds[psz] == 0) { + if (ds[psz] == -4 || ds[psz] == -6 || ds[psz] == -7 || ds[psz] == -11 || ds[psz] == 0) { // printf("%s ds[%d]=%d vdata_offset=%d\n", memb_func[type].sym->name, psz, ds[psz], // vdata_offset); vdata_offset += ml->nodecount; @@ -332,6 +332,7 @@ int nrnthread_dat2_mech(int tid, int*& nodeindices, double*& data, int*& pdata, + std::vector& nmodlrandom, // 5 uint32_t per var per instance std::vector& pointer2type) { if (tid >= nrn_nthread) { return 0; @@ -392,6 +393,32 @@ int nrnthread_dat2_mech(int tid, pdata = NULL; } + // nmodlrandom: reserve 5 uint32 for each var of each instance + // id1, id2, id3, seq, uint32_t(which) + // Header is number of random variables followed by dparam indices + // if no destructor, skip. There are no random variables. + if (nrn_mech_inst_destruct.count(type)) { + auto& indices = nrn_mech_random_indices(type); + nmodlrandom.reserve(1 + indices.size() + 5 * n * indices.size()); + nmodlrandom.push_back(indices.size()); + for (int ix: indices) { + nmodlrandom.push_back((uint32_t) ix); + } + for (int ix: indices) { + uint32_t data[5]; + char which; + for (int i = 0; i < n; ++i) { + auto& datum = ml->pdata[i][ix]; + nrnran123_State* r = (nrnran123_State*) datum.get(); + nrnran123_getids3(r, &data[0], &data[1], &data[2]); + nrnran123_getseq(r, &data[3], &which); + data[4] = uint32_t(which); + for (auto j: data) { + nmodlrandom.push_back(j); + } + } + } + } return 1; } @@ -523,6 +550,37 @@ int core2nrn_corepointer_mech(int tid, int type, int icnt, int dcnt, int* iArray return 1; } +// NMODL RANDOM seq34 data return from coreneuron +int core2nrn_nmodlrandom(int tid, + int type, + const std::vector& indices, + const std::vector& nmodlrandom) { + if (tid >= nrn_nthread) { + return 0; + } + NrnThread& nt = nrn_threads[tid]; + Memb_list* ml = nt._ml_list[type]; + // ARTIFICIAL_CELL are not in nt. + if (!ml) { + ml = CellGroup::deferred_type2artml_[tid][type]; + assert(ml); + } + + auto& nrnindices = nrn_mech_random_indices(type); // for sanity checking + assert(nrnindices == indices); + assert(nmodlrandom.size() == indices.size() * ml->nodecount); + + int ir = 0; // into nmodlrandom + for (const auto ix: nrnindices) { + for (int i = 0; i < ml->nodecount; ++i) { + auto& datum = ml->pdata[i][ix]; + nrnran123_State* state = (nrnran123_State*) datum.get(); + nrnran123_setseq(state, nmodlrandom[ir++]); + } + } + return 1; +} + int* datum2int(int type, Memb_list* ml, NrnThread& nt, @@ -577,6 +635,8 @@ int* datum2int(int type, } else if (etype == -7) { // bbcorepointer pdata[jj] = ml_vdata_offset + eindex; // printf("etype %d jj=%d eindex=%d pdata=%d\n", etype, jj, eindex, pdata[jj]); + } else if (etype == -11) { // random + pdata[jj] = ml_vdata_offset + eindex; } else { // uninterpreted assert(eindex != -3); // avoided if last pdata[jj] = 0; diff --git a/src/nrniv/nrncore_write/callbacks/nrncore_callbacks.h b/src/nrniv/nrncore_write/callbacks/nrncore_callbacks.h index ddc07a1520..efce1359fb 100644 --- a/src/nrniv/nrncore_write/callbacks/nrncore_callbacks.h +++ b/src/nrniv/nrncore_write/callbacks/nrncore_callbacks.h @@ -6,6 +6,8 @@ #include #include #include +#include + // includers need several pieces of info for nrn_get_partrans_setup_info #include "partrans.h" @@ -66,6 +68,7 @@ int nrnthread_dat2_mech(int tid, int*& nodeindices, double*& data, int*& pdata, + std::vector& nmodlrandom, std::vector& pointer2type); int nrnthread_dat2_3(int tid, int nweight, @@ -123,6 +126,10 @@ int nrnthread_all_spike_vectors_return(std::vector& spiketvec, void nrnthreads_all_weights_return(std::vector& weights); size_t nrnthreads_type_return(int type, int tid, double*& data, std::vector& mdata); int core2nrn_corepointer_mech(int tid, int type, int icnt, int dcnt, int* iarray, double* darray); +int core2nrn_nmodlrandom(int tid, + int type, + const std::vector& indices, + const std::vector& nmodlrandom); } // For direct transfer of event queue information from CoreNEURON @@ -227,6 +234,7 @@ static core2nrn_callback_t cnbs[] = { {"core2nrn_clear_queues_", (CNB) core2nrn_clear_queues}, {"core2nrn_corepointer_mech_", (CNB) core2nrn_corepointer_mech}, + {"core2nrn_nmodlrandom_", (CNB) core2nrn_nmodlrandom}, {"core2nrn_NetCon_event_", (CNB) core2nrn_NetCon_event}, {"core2nrn_SelfEvent_event_", (CNB) core2nrn_SelfEvent_event}, {"core2nrn_SelfEvent_event_noweight_", (CNB) core2nrn_SelfEvent_event_noweight}, diff --git a/src/nrniv/nrncore_write/data/cell_group.cpp b/src/nrniv/nrncore_write/data/cell_group.cpp index fae073634a..fe25bc1793 100644 --- a/src/nrniv/nrncore_write/data/cell_group.cpp +++ b/src/nrniv/nrncore_write/data/cell_group.cpp @@ -258,7 +258,7 @@ void CellGroup::datumindex_fill(int ith, CellGroup& cg, DatumIndices& di, Memb_l int vdata_size = 0; for (int i = 0; i < dsize; ++i) { int* ds = memb_func[di.type].dparam_semantics; - if (ds[i] == -4 || ds[i] == -6 || ds[i] == -7 || ds[i] == 0) { + if (ds[i] == -4 || ds[i] == -6 || ds[i] == -7 || ds[i] == -11 || ds[i] == 0) { ++vdata_size; } } @@ -311,6 +311,9 @@ void CellGroup::datumindex_fill(int ith, CellGroup& cg, DatumIndices& di, Memb_l } else if (dmap[j] == -10) { // fornetcon etype = -10; eindex = 0; + } else if (dmap[j] == -11) { // random + etype = -11; + eindex = vdata_offset++; } else if (dmap[j] == -9) { // diam cg.ndiam = nt.end; etype = -9; diff --git a/src/nrniv/nrncore_write/io/nrncore_io.cpp b/src/nrniv/nrncore_write/io/nrncore_io.cpp index 094e727093..4b3fa7f709 100644 --- a/src/nrniv/nrncore_write/io/nrncore_io.cpp +++ b/src/nrniv/nrncore_write/io/nrncore_io.cpp @@ -22,7 +22,7 @@ extern NetCvode* net_cvode_instance; extern void (*nrnthread_v_transfer_)(NrnThread*); int chkpnt; -const char* bbcore_write_version = "1.6"; // Allow muliple gid and PreSyn per real cell. +const char* bbcore_write_version = "1.7"; // NMODLRandom /// create directory with given path void create_dir_path(const std::string& path) { @@ -206,7 +206,9 @@ void write_nrnthread(const char* path, NrnThread& nt, CellGroup& cg) { int *nodeindices = NULL, *pdata = NULL; double* data = NULL; std::vector pointer2type; - nrnthread_dat2_mech(nt.id, i, dsz_inst, nodeindices, data, pdata, pointer2type); + std::vector nmodlrandom; + nrnthread_dat2_mech( + nt.id, i, dsz_inst, nodeindices, data, pdata, nmodlrandom, pointer2type); Memb_list* ml = mla[i].second; int n = ml->nodecount; int sz = nrn_prop_param_size_[type]; @@ -227,6 +229,11 @@ void write_nrnthread(const char* path, NrnThread& nt, CellGroup& cg) { if (sz > 0) { writeint(pointer2type.data(), sz); } + + fprintf(f, "%d nmodlrandom\n", int(nmodlrandom.size())); + if (nmodlrandom.size()) { + write_uint32vec(nmodlrandom, f); + } } } @@ -300,6 +307,12 @@ void writedbl_(double* p, size_t size, FILE* f) { assert(n == size); } +void write_uint32vec(std::vector& vec, FILE* f) { + fprintf(f, "chkpnt %d\n", chkpnt++); + size_t n = fwrite(vec.data(), sizeof(uint32_t), vec.size(), f); + assert(n == vec.size()); +} + #define writeint(p, size) writeint_(p, size, f) #define writedbl(p, size) writedbl_(p, size, f) diff --git a/src/nrniv/nrncore_write/io/nrncore_io.h b/src/nrniv/nrncore_write/io/nrncore_io.h index ee690b694d..d191de5f89 100644 --- a/src/nrniv/nrncore_write/io/nrncore_io.h +++ b/src/nrniv/nrncore_write/io/nrncore_io.h @@ -31,6 +31,8 @@ void write_nrnthread(const char* fname, NrnThread& nt, CellGroup& cg); void writeint_(int* p, size_t size, FILE* f); void writedbl_(double* p, size_t size, FILE* f); +void write_uint32vec(std::vector& vec, FILE* f); + #define writeint(p, size) writeint_(p, size, f) #define writedbl(p, size) writedbl_(p, size, f) // also for read diff --git a/src/nrniv/nrnmenu.cpp b/src/nrniv/nrnmenu.cpp index cfe20eea35..463ab45280 100644 --- a/src/nrniv/nrnmenu.cpp +++ b/src/nrniv/nrnmenu.cpp @@ -479,7 +479,11 @@ static void point_menu(Object* ob, int make_label) { if (psym->s_varn) { for (k = 0; k < psym->s_varn; k++) { vsym = psym->u.ppsym[k]; - if (nrn_vartype(vsym) == nrnocCONST) { + int vartype = nrn_vartype(vsym); + if (vartype == NMODLRANDOM) { // skip + continue; + } + if (vartype == nrnocCONST) { deflt = true; #if defined(MikeNeubig) diff --git a/src/nrnoc/cabcode.cpp b/src/nrnoc/cabcode.cpp index 5fdc65fd56..8cc9ac7a07 100644 --- a/src/nrnoc/cabcode.cpp +++ b/src/nrnoc/cabcode.cpp @@ -1396,6 +1396,24 @@ void rangepoint(void) /* symbol at pc, return value on stack */ rangevareval(); } +void rangeobjeval(void) /* symbol at pc, section location on stack, return object on stack*/ +{ + Symbol* s{(pc++)->sym}; + assert(s->subtype == NMODLRANDOM); // the only possibility at the moment + double d = xpop(); + Section* sec{nrn_sec_pop()}; + auto const i = node_index(sec, d); + Prop* m = nrn_mechanism_check(s->u.rng.type, sec, i); + Object* ob = nrn_nmodlrandom_wrap(m, s); + hoc_push_object(ob); +} + +void rangeobjevalmiddle(void) /* symbol at pc, return object on stack*/ +{ + hoc_pushx(0.5); + rangeobjeval(); +} + int node_index(Section* sec, double x) /* returns nearest index to x */ { int i; diff --git a/src/nrnoc/init.cpp b/src/nrnoc/init.cpp index 4feddaa3df..0b09f07af4 100644 --- a/src/nrnoc/init.cpp +++ b/src/nrnoc/init.cpp @@ -21,6 +21,7 @@ #include "nrnmpi.h" #include +#include /* change this to correspond to the ../nmodl/nocpout nmodl_version_ string*/ static char nmodl_version_[] = "7.7.0"; @@ -159,6 +160,7 @@ int* nrn_prop_dparam_size_; int* nrn_dparam_ptr_start_; int* nrn_dparam_ptr_end_; NrnWatchAllocateFunc_t* nrn_watch_allocate_; +std::unordered_map nrn_mech_inst_destruct; void hoc_reg_watch_allocate(int type, NrnWatchAllocateFunc_t waf) { nrn_watch_allocate_[type] = waf; @@ -463,7 +465,7 @@ void initialize_memb_func(int mechtype, nrn_init_t initialize, int vectorized); void check_mech_version(const char** m); -std::pair count_variables_in_mechanism(const char** m2, int modltypemax); +int count_variables_in_mechanism(const char** m2, int modltypemax); void register_mech_vars(const char** var_buffers, int modltypemax, Symbol* mech_symbol, @@ -503,9 +505,9 @@ void nrn_register_mech_common(const char** m, } else { modltypemax = NRNPOINTER; } - auto [nvartypes, nvars] = count_variables_in_mechanism(m2, modltypemax); + int nvars = count_variables_in_mechanism(m2, modltypemax); mech_symbol->s_varn = nvars; - mech_symbol->u.ppsym = (Symbol**) emalloc((unsigned) (nvartypes * sizeof(Symbol*))); + mech_symbol->u.ppsym = (Symbol**) emalloc((unsigned) (nvars * sizeof(Symbol*))); register_mech_vars(m2, modltypemax, mech_symbol, mechtype, nrnpointerindex); ++mechtype; @@ -552,16 +554,8 @@ void register_mech_vars(const char** var_buffers, nsub = 1; varname.erase(subscript); } - /*SUPPRESS 624*/ if ((var_symbol = hoc_lookup(varname.c_str()))) { -#if 0 - if (var_symbol->subtype != RANGEVAR) { - IGNORE(fprintf(stderr, CHKmes, - varname.c_str())); - } -#else // not 0 IGNORE(fprintf(stderr, CHKmes, varname.c_str())); -#endif // not 0 } else { var_symbol = hoc_install(varname.c_str(), RANGEVAR, 0.0, &hoc_symlist); var_symbol->subtype = modltype; @@ -595,16 +589,18 @@ void register_mech_vars(const char** var_buffers, } } -std::pair count_variables_in_mechanism(const char** m2, int modltypemax) { - int j, k, modltype; +int count_variables_in_mechanism(const char** m2, int modltypemax) { + int j; + int modltype; + int nvars; // count the number of variables registered in this mechanism - for (j = 0, k = 0, modltype = nrnocCONST; modltype <= modltypemax; modltype++) { + for (j = 0, nvars = 0, modltype = nrnocCONST; modltype <= modltypemax; modltype++) { // while we have not encountered a 0 (sentinel for variable type) while (m2[j++]) { - k++; + nvars++; } } - return std::make_pair(j, k); + return nvars; } void reallocate_mech_data(int mechtype) { @@ -772,30 +768,24 @@ namespace { * * This logic used to live inside hoc_register_dparam_semantics. */ + +// name to int map for the negative types +// xx_ion and #xx_ion will get values of type and type+1000 respectively +static std::unordered_map name_to_negint = {{"area", -1}, + {"iontype", -2}, + {"cvodeieq", -3}, + {"netsend", -4}, + {"pointer", -5}, + {"pntproc", -6}, + {"bbcorepointer", -7}, + {"watch", -8}, + {"diam", -9}, + {"fornetcon", -10}, + {"random", -11}}; + int dparam_semantics_to_int(std::string_view name) { - // only interested in area, iontype, cvode_ieq, netsend, pointer, pntproc, bbcorepointer, watch, - // diam, fornetcon, xx_ion and #xx_ion which will get a semantics value of -1, -2, -3, -4, -5, - // -6, -7, -8, -9, -10 type, and type+1000 respectively - if (name == "area") { - return -1; - } else if (name == "iontype") { - return -2; - } else if (name == "cvodeieq") { - return -3; - } else if (name == "netsend") { - return -4; - } else if (name == "pointer") { - return -5; - } else if (name == "pntproc") { - return -6; - } else if (name == "bbcorepointer") { - return -7; - } else if (name == "watch") { - return -8; - } else if (name == "diam") { - return -9; - } else if (name == "fornetcon") { - return -10; + if (auto got = name_to_negint.find(std::string{name}); got != name_to_negint.end()) { + return got->second; } else { bool const i{name[0] == '#'}; Symbol* s = hoc_lookup(std::string{name.substr(i)}.c_str()); @@ -805,6 +795,57 @@ int dparam_semantics_to_int(std::string_view name) { throw std::runtime_error("unknown dparam semantics: " + std::string{name}); } } + +std::vector indices_of_type( + const char* semantic_type, + std::vector> const& dparam_info) { + std::vector indices{}; + int inttype = dparam_semantics_to_int(std::string{semantic_type}); + for (auto i = 0; i < dparam_info.size(); ++i) { + if (dparam_semantics_to_int(dparam_info[i].second) == inttype) { + indices.push_back(i); + } + } + return indices; +} + +static std::unordered_map> mech_random_indices{}; + +void update_mech_ppsym_for_modlrandom( + int mechtype, + std::vector> const& dparam_info) { + std::vector indices = indices_of_type("random", dparam_info); + mech_random_indices[mechtype] = indices; + if (indices.empty()) { + return; + } + Symbol* mechsym = memb_func[mechtype].sym; + int is_point = memb_func[mechtype].is_point; + + int k = mechsym->s_varn; + mechsym->s_varn += int(indices.size()); + mechsym->u.ppsym = (Symbol**) erealloc(mechsym->u.ppsym, mechsym->s_varn * sizeof(Symbol*)); + + + for (auto i: indices) { + auto& p = dparam_info[i]; + Symbol* ransym{}; + if (is_point) { + ransym = hoc_install(p.first, RANGEOBJ, 0.0, &(nrn_pnt_template_[mechtype]->symtable)); + } else { + std::string s{p.first}; + s += "_"; + s += mechsym->name; + ransym = hoc_install(s.c_str(), RANGEOBJ, 0.0, &hoc_symlist); + } + ransym->subtype = NMODLRANDOM; + ransym->u.rng.type = mechtype; + ransym->cpublic = 1; + ransym->u.rng.index = i; + mechsym->u.ppsym[k++] = ransym; + } +} + } // namespace namespace neuron::mechanism::detail { @@ -842,6 +883,7 @@ void register_data_fields(int mechtype, // of double-valued // per-instance variables memb_list[mechtype].set_storage_pointer(&mech_data); + update_mech_ppsym_for_modlrandom(mechtype, dparam_info); } } // namespace neuron::mechanism::detail namespace neuron::mechanism { @@ -909,6 +951,17 @@ void hoc_register_dparam_semantics(int mechtype, int ix, const char* name) { assert(memb_func[mechtype].dparam_semantics[ix] == dparam_semantics_to_int(name)); } +int nrn_dparam_semantics_to_int(const char* name) { + return dparam_semantics_to_int(name); +} + +/** + * @brief dparam indices with random semantics for mechtype + */ +std::vector& nrn_mech_random_indices(int type) { + return mech_random_indices[type]; +} + void hoc_register_cvode(int i, nrn_ode_count_t cnt, nrn_ode_map_t map, diff --git a/src/nrnoc/membfunc.h b/src/nrnoc/membfunc.h index 916d871f7f..5cf6ca55bf 100644 --- a/src/nrnoc/membfunc.h +++ b/src/nrnoc/membfunc.h @@ -160,8 +160,8 @@ inline std::size_t vext_pseudoindex() { pointers which connect variables from other mechanisms via the _ppval array. \ */ -#define _AMBIGUOUS 5 // for Ions -#define RAND 6 // RANDOM +#define _AMBIGUOUS 5 // for Ions +#define NMODLRANDOM 6 // RANDOM variable in NEURON block #define BEFORE_INITIAL 0 #define AFTER_INITIAL 1 @@ -180,6 +180,8 @@ extern std::vector memb_func; extern int n_memb_func; extern int* nrn_prop_param_size_; extern int* nrn_prop_dparam_size_; +extern int nrn_dparam_semantics_to_int(const char*); +extern std::vector& nrn_mech_random_indices(int type); extern std::vector memb_list; /* for finitialize, order is same up through extracellular, then ions, @@ -312,3 +314,7 @@ _nrn_mechanism_get_param_handle(Prop* prop, int field, int array_index = 0) { [[nodiscard]] NrnThread* _nrn_mechanism_get_thread(Node*); [[nodiscard]] int _nrn_mechanism_get_type(Prop*); [[nodiscard]] int _nrn_mechanism_get_v_node_index(Node*); + +// Rarely (e.g. NEURON {RANDOM123 ranvar}) instances of a mod file +// need to deallocate owning objects at end of their life. +extern std::unordered_map nrn_mech_inst_destruct; diff --git a/src/nrnoc/netstim.mod b/src/nrnoc/netstim.mod index 5db68ca1f1..86bb7562fe 100755 --- a/src/nrnoc/netstim.mod +++ b/src/nrnoc/netstim.mod @@ -1,26 +1,12 @@ : $Id: netstim.mod 2212 2008-09-08 14:32:26Z hines $ : comments at end -: the Random idiom has been extended to support CoreNEURON. - -: For backward compatibility, noiseFromRandom(hocRandom) can still be used -: as well as the default low-quality scop_exprand generator. -: However, CoreNEURON will not accept usage of the low-quality generator, -: and, if noiseFromRandom is used to specify the random stream, that stream -: must be using the Random123 generator. - -: The recommended idiom for specfication of the random stream is to use -: noiseFromRandom123(id1, id2[, id3]) - -: If any instance uses noiseFromRandom123, then no instance can use noiseFromRandom -: and vice versa. - NEURON { ARTIFICIAL_CELL NetStim + THREADSAFE RANGE interval, number, start RANGE noise - THREADSAFE : only true if every instance has its own distinct Random - BBCOREPOINTER donotuse + RANDOM ranvar } PARAMETER { @@ -34,41 +20,10 @@ ASSIGNED { event (ms) on ispike - donotuse -} - -VERBATIM -#if !NRNBBCORE -/** If we're running in NEURON, specify the noise style for all instances. - * 1 means noiseFromRandom was called when _ran_compat was previously 0. - * 2 means noiseFromRandom123 was called when _ran_compat was previously 0. - */ -static int _ran_compat; -#endif -ENDVERBATIM - -:backward compatibility -PROCEDURE seed(x) { -VERBATIM -#if !NRNBBCORE -ENDVERBATIM - set_seed(x) -VERBATIM -#endif -ENDVERBATIM } INITIAL { - VERBATIM -#if NRNBBCORE - if(_p_donotuse) { -#else - if(_p_donotuse && _ran_compat == 2) { -#endif - /* only this style initializes the stream on finitialize */ - nrnran123_setseq(reinterpret_cast(_p_donotuse), 0, 0); - } - ENDVERBATIM + seed(0) on = 0 : off ispike = 0 if (noise < 0) { @@ -110,180 +65,9 @@ FUNCTION invl(mean (ms)) (ms) { } FUNCTION erand() { -VERBATIM - if (_p_donotuse) { - /* - :Supports separate independent but reproducible streams for - : each instance. However, the corresponding hoc Random - : distribution MUST be set to Random.negexp(1) - */ -#if !NRNBBCORE - if (_ran_compat == 2) { - _lerand = nrnran123_negexp(reinterpret_cast(_p_donotuse)); - } else { - _lerand = nrn_random_pick(reinterpret_cast(_p_donotuse)); - } -#else - _lerand = nrnran123_negexp(reinterpret_cast(_p_donotuse)); -#endif - return _lerand; - } else { -#if NRNBBCORE - assert(0); -#else - /* - : the old standby. Cannot use if reproducible parallel sim - : independent of nhost or which host this instance is on - : is desired, since each instance on this cpu draws from - : the same stream - */ -#endif - } -#if !NRNBBCORE -ENDVERBATIM - erand = exprand(1) -VERBATIM -#endif -ENDVERBATIM -} - -PROCEDURE noiseFromRandom() { -VERBATIM -#if !NRNBBCORE - { - if (_ran_compat == 2) { - fprintf(stderr, "NetStim.noiseFromRandom123 was previously called\n"); - assert(0); - } - _ran_compat = 1; - auto& randstate = reinterpret_cast(_p_donotuse); - if (ifarg(1)) { - randstate = nrn_random_arg(1); - } else { - randstate = nullptr; - } - } -#endif -ENDVERBATIM -} - -PROCEDURE noiseFromRandom123() { -VERBATIM -#if !NRNBBCORE - if (_ran_compat == 1) { - fprintf(stderr, "NetStim.noiseFromRandom was previously called\n"); - assert(0); - } - _ran_compat = 2; - auto& r123state = reinterpret_cast(_p_donotuse); - if (r123state) { - nrnran123_deletestream(r123state); - r123state = nullptr; - } - if (ifarg(3)) { - r123state = nrnran123_newstream3(static_cast(*getarg(1)), static_cast(*getarg(2)), static_cast(*getarg(3))); - } else if (ifarg(2)) { - r123state = nrnran123_newstream(static_cast(*getarg(1)), static_cast(*getarg(2))); - } -#endif -ENDVERBATIM -} - -DESTRUCTOR { -VERBATIM - if (!noise) { return; } - if (_p_donotuse) { -#if NRNBBCORE - { /* but note that mod2c does not translate DESTRUCTOR */ -#else - if (_ran_compat == 2) { -#endif - auto& r123state = reinterpret_cast(_p_donotuse); - nrnran123_deletestream(r123state); - r123state = nullptr; - } - } -ENDVERBATIM -} - -VERBATIM -static void bbcore_write(double* x, int* d, int* xx, int *offset, _threadargsproto_) { - if (!noise) { return; } - /* error if using the legacy scop_exprand */ - if (!_p_donotuse) { - fprintf(stderr, "NetStim: cannot use the legacy scop_negexp generator for the random stream.\n"); - assert(0); - } - if (d) { - char which; - uint32_t* di = reinterpret_cast(d) + *offset; -#if !NRNBBCORE - if (_ran_compat == 1) { - auto* rand = reinterpret_cast(_p_donotuse); - /* error if not using Random123 generator */ - if (!nrn_random_isran123(rand, di, di+1, di+2)) { - fprintf(stderr, "NetStim: Random123 generator is required\n"); - assert(0); - } - nrn_random123_getseq(rand, di+3, &which); - di[4] = which; - } else { -#else - { -#endif - auto& r123state = reinterpret_cast(_p_donotuse); - nrnran123_getids3(r123state, di, di+1, di+2); - nrnran123_getseq(r123state, di+3, &which); - di[4] = which; -#if NRNBBCORE - /* CoreNEURON does not call DESTRUCTOR so... */ - nrnran123_deletestream(r123state); - r123state = nullptr; -#endif - } - /*printf("Netstim bbcore_write %d %d %d\n", di[0], di[1], di[3]);*/ - } - *offset += 5; + erand = random_negexp(ranvar) } -static void bbcore_read(double* x, int* d, int* xx, int* offset, _threadargsproto_) { - if (!noise) { return; } - /* Generally, CoreNEURON, in the context of psolve, begins with an empty - * model, so this call takes place in the context of a freshly created - * instance and _p_donotuse is not NULL. - * However, this function is also now called from NEURON at the end of - * coreneuron psolve in order to transfer back the nrnran123 sequence state. - * That allows continuation with a subsequent psolve within NEURON or - * properly transfer back to CoreNEURON if we continue the psolve there. - * So now, extra logic is needed for this call to work in a NEURON context. - */ - uint32_t* di = reinterpret_cast(d) + *offset; -#if NRNBBCORE - auto& r123state = reinterpret_cast(_p_donotuse); - assert(!r123state); - r123state = nrnran123_newstream3(di[0], di[1], di[2]); - nrnran123_setseq(r123state, di[3], di[4]); -#else - uint32_t id1, id2, id3; - assert(_p_donotuse); - if (_ran_compat == 1) { /* Hoc Random.Random123 */ - auto* pv = reinterpret_cast(_p_donotuse); - int b = nrn_random_isran123(pv, &id1, &id2, &id3); - assert(b); - nrn_random123_setseq(pv, di[3], (char)di[4]); - } else { - assert(_ran_compat == 2); - auto* r123state = reinterpret_cast(_p_donotuse); - nrnran123_getids3(r123state, &id1, &id2, &id3); - nrnran123_setseq(r123state, di[3], di[4]); - } - /* Random123 on NEURON side has same ids as on CoreNEURON side */ - assert(di[0] == id1 && di[1] == id2 && di[2] == id3); -#endif - *offset += 5; -} -ENDVERBATIM - PROCEDURE next_invl() { if (number > 0) { event = invl(interval) @@ -323,34 +107,28 @@ NET_RECEIVE (w) { } } -FUNCTION bbsavestate() { - bbsavestate = 0 - : limited to noiseFromRandom123 -VERBATIM -#if !NRNBBCORE - if (_ran_compat == 2) { - auto r123state = reinterpret_cast(_p_donotuse); - if (!r123state) { return 0.0; } - double* xdir = hoc_pgetarg(1); - if (*xdir == -1.) { - *xdir = 2; - return 0.0; - } - double* xval = hoc_pgetarg(2); - if (*xdir == 0.) { - char which; - uint32_t seq; - nrnran123_getseq(r123state, &seq, &which); - xval[0] = seq; - xval[1] = which; - } - if (*xdir == 1) { - nrnran123_setseq(r123state, xval[0], xval[1]); - } - } -#endif -ENDVERBATIM -} +:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +: Legacy API +: +: Difference: seed(x) merely sets ranvar sequence to ((uint32_t)x, 0) +: noiseFromRandom HOC Random object must use Random123 +: generator. The ids and sequence are merely copied +: into ranvar. +:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + +: the Random idiom has been extended to support CoreNEURON. + +: For backward compatibility, noiseFromRandom(hocRandom) can still be used +: as well as the default low-quality scop_exprand generator. +: However, CoreNEURON will not accept usage of the low-quality generator, +: and, if noiseFromRandom is used to specify the random stream, that stream +: must be using the Random123 generator. + +: The recommended idiom for specfication of the random stream is to use +: noiseFromRandom123(id1, id2[, id3]) + +: If any instance uses noiseFromRandom123, then no instance can use noiseFromRandom +: and vice versa. COMMENT @@ -388,3 +166,39 @@ its sequence. ENDCOMMENT +PROCEDURE seed(x) { + random_setseq(ranvar, x) +} + +PROCEDURE noiseFromRandom() { +VERBATIM +#if !NRNBBCORE + { + if (ifarg(1)) { + Rand* r = nrn_random_arg(1); + uint32_t id[3]; + if (!nrn_random_isran123(r, &id[0], &id[1], &id[2])) { + hoc_execerr_ext("NetStim: Random.Random123 generator is required."); + } + nrnran123_setids(ranvar, id[0], id[1], id[2]); + char which; + nrn_random123_getseq(r, &id[0], &which); + nrnran123_setseq(ranvar, id[0], which); + } + } +#endif +ENDVERBATIM +} + +PROCEDURE noiseFromRandom123() { +VERBATIM +#if !NRNBBCORE + if (ifarg(3)) { + nrnran123_setids(ranvar, static_cast(*getarg(1)), static_cast(*getarg(2)), static_cast(*getarg(3))); + } else if (ifarg(2)) { + nrnran123_setids(ranvar, static_cast(*getarg(1)), static_cast(*getarg(2)), 0); + } + nrnran123_setseq(ranvar, 0, 0); +#endif +ENDVERBATIM +} diff --git a/src/nrnoc/point.cpp b/src/nrnoc/point.cpp index e9a7b95e30..4159358d6b 100644 --- a/src/nrnoc/point.cpp +++ b/src/nrnoc/point.cpp @@ -328,6 +328,9 @@ static void free_one_point(Point_process* pnt) { if (memb_func[p->_type].destructor) { memb_func[p->_type].destructor(p); } + if (auto got = nrn_mech_inst_destruct.find(p->_type); got != nrn_mech_inst_destruct.end()) { + (got->second)(p); + } if (p->dparam) { nrn_prop_datum_free(p->_type, p->dparam); } diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index b3b698851c..0849d1cdb5 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -721,6 +721,9 @@ void single_prop_free(Prop* p) { clear_point_process_struct(p); return; } + if (auto got = nrn_mech_inst_destruct.find(p->_type); got != nrn_mech_inst_destruct.end()) { + (got->second)(p); + } if (p->dparam) { if (p->_type == CABLESECTION) { notify_freed_val_array(&(p->dparam[2].literal_value()), 6); diff --git a/src/nrnpython/nrnpy_hoc.cpp b/src/nrnpython/nrnpy_hoc.cpp index 678a4977fa..19e8531a1b 100644 --- a/src/nrnpython/nrnpy_hoc.cpp +++ b/src/nrnpython/nrnpy_hoc.cpp @@ -1185,7 +1185,7 @@ static PyObject* hocobj_getattr(PyObject* subself, PyObject* pyname) { // an array int t = sym->type; if (t == VAR || t == STRING || t == OBJECTVAR || t == RANGEVAR || t == SECTION || - t == SECTIONREF || t == VARALIAS || t == OBJECTALIAS) { + t == SECTIONREF || t == VARALIAS || t == OBJECTALIAS || t == RANGEOBJ) { if (sym != nrn_child_sym && !ISARRAY(sym)) { hoc_push_object(po->ho_); nrn_inpython_ = 1; diff --git a/src/nrnpython/nrnpy_nrn.cpp b/src/nrnpython/nrnpy_nrn.cpp index e1f82848bf..c875edc10e 100644 --- a/src/nrnpython/nrnpy_nrn.cpp +++ b/src/nrnpython/nrnpy_nrn.cpp @@ -1931,7 +1931,13 @@ static PyObject* segment_getattro(NPySegObj* self, PyObject* pyname) { } } else if ((rv = PyDict_GetItemString(rangevars_, n)) != NULL) { sym = ((NPyRangeVar*) rv)->sym_; - if (ISARRAY(sym)) { + if (sym->type == RANGEOBJ) { + int mtype = sym->u.rng.type; + Node* nd = node_exact(sec, self->x_); + Prop* p = nrn_mechanism(mtype, nd); + Object* ob = nrn_nmodlrandom_wrap(p, sym); + result = nrnpy_ho2po(ob); + } else if (ISARRAY(sym)) { NPyRangeVar* r = PyObject_New(NPyRangeVar, range_type); r->pymech_ = new_pymechobj(); r->pymech_->pyseg_ = self; @@ -2150,7 +2156,7 @@ static PyObject* mech_getattro(NPyMechObj* self, PyObject* pyname) { std::snprintf(buf, bufsz, "%s_%s", isptr ? n + 5 : n, mname); } Symbol* sym = np.find(buf); - if (sym) { + if (sym && sym->type == RANGEVAR) { // printf("mech_getattro sym %s\n", sym->name); if (ISARRAY(sym)) { NPyRangeVar* r = PyObject_New(NPyRangeVar, range_type); @@ -2170,6 +2176,9 @@ static PyObject* mech_getattro(NPyMechObj* self, PyObject* pyname) { result = Py_BuildValue("d", *px); } } + } else if (sym && sym->type == RANGEOBJ) { + Object* ob = nrn_nmodlrandom_wrap(self->prop_, sym); + result = nrnpy_ho2po(ob); } else if (strcmp(n, "__dict__") == 0) { result = PyDict_New(); for (Symbol* s = np.first_var(); np.more_var(); s = np.next_var()) { @@ -2614,7 +2623,7 @@ static PyMethodDef nrnpy_methods[] = { static PyObject* nrnmodule_; static void rangevars_add(Symbol* sym) { - assert(sym && sym->type == RANGEVAR); + assert(sym && (sym->type == RANGEVAR || sym->type == RANGEOBJ)); NPyRangeVar* r = PyObject_New(NPyRangeVar, range_type); // printf("%s\n", sym->name); r->sym_ = sym; diff --git a/src/oc/code.h b/src/oc/code.h index cf280b4ee8..7ecfbb0c82 100644 --- a/src/oc/code.h +++ b/src/oc/code.h @@ -38,6 +38,7 @@ extern void hoc_autoobject(void), hocobjret(void), hoc_newobj_ret(void); extern void connectsection(void), add_section(void), range_const(void), range_interpolate(void); extern void clear_sectionlist(void), install_sectionlist(void); extern void rangevareval(void), sec_access(void), mech_access(void); +extern void rangeobjeval(void), rangeobjevalmiddle(void); extern void for_segment(void), for_segment1(void); extern void sec_access_temp(void), sec_access_push(void), sec_access_pop(void); extern void rangepoint(void), forall_section(void), hoc_ifsec(void); diff --git a/src/oc/hoc_oop.cpp b/src/oc/hoc_oop.cpp index af1c88e7d2..c56f418654 100644 --- a/src/oc/hoc_oop.cpp +++ b/src/oc/hoc_oop.cpp @@ -948,6 +948,22 @@ static void range_suffix(Symbol* sym, int nindex, int narg) { hoc_execerror(sym->name, "section property can't have argument"); } hoc_pushs(sym); + } else if (sym->type == RANGEOBJ) { + // must return NMODLObject on stack + assert(sym->subtype == NMODLRANDOM); // the only possibility at the moment + double x{0.5}; + if (narg) { + if (narg > 1) { + hoc_execerr_ext("%s range object can have only one arg length parameter", + sym->name); + } + x = xpop(); + } + Section* sec{nrn_sec_pop()}; + auto const i = node_index(sec, x); + Prop* m = nrn_mechanism_check(sym->u.rng.type, sec, i); + Object* ob = nrn_nmodlrandom_wrap(m, sym); + hoc_push_object(ob); } else { hoc_execerror(sym->name, "suffix not a range variable or section property"); } @@ -1263,6 +1279,17 @@ void hoc_object_component() { hoc_nopop(); /* get rid of iterator statement context */ break; } + case RANGEOBJ: { + assert(sym->subtype == NMODLRANDOM); + if (sym->subtype == NMODLRANDOM) { // NMODL NEURON block RANDOM var + // RANGE type. The void* is a nrnran123_State*. Wrap in a + // NMODLRandom and push_object + Object* o = nrn_pntproc_nmodlrandom_wrap(obp->u.this_pointer, sym); + hoc_pop_defer(); + hoc_push_object(o); + } + break; + } default: if (cplus) { if (nindex) { diff --git a/src/oc/oc_ansi.h b/src/oc/oc_ansi.h index c83f9300fd..1499f77512 100644 --- a/src/oc/oc_ansi.h +++ b/src/oc/oc_ansi.h @@ -39,6 +39,7 @@ union Objectdata; struct Symbol; struct Symlist; struct VoidFunc; +struct Prop; namespace neuron { struct model_sorted_token; @@ -312,6 +313,8 @@ Object** hoc_temp_objvar(Symbol* template_symbol, void* cpp_object); Object** hoc_temp_objptr(Object*); Object* hoc_new_object(Symbol* symtemp, void* v); void hoc_new_object_asgn(Object** obp, Symbol* template_symbol, void* cpp_object); +Object* nrn_pntproc_nmodlrandom_wrap(void* pnt, Symbol* sym); +Object* nrn_nmodlrandom_wrap(Prop* prop, Symbol* sym); HocSymExtension* hoc_var_extra(const char*); double check_domain_limits(float*, double); Object* hoc_obj_get(int i); diff --git a/src/oc/parse.ypp b/src/oc/parse.ypp index 8b38a4629b..70cccbbce0 100755 --- a/src/oc/parse.ypp +++ b/src/oc/parse.ypp @@ -112,9 +112,9 @@ static void hoc_opasgn_invalid(int op); /* NEWCABLE */ %token SECTIONKEYWORD SECTION CONNECTKEYWORD ACCESSKEYWORD %token RANGEVAR MECHANISM INSERTKEYWORD FORALL NRNPNTVAR FORSEC IFSEC -%token UNINSERTKEYWORD SETPOINTERKEYWORD SECTIONREF +%token UNINSERTKEYWORD SETPOINTERKEYWORD SECTIONREF RANGEOBJ %type sectiondecl sectionname -%type rangevar rangevar1 section section_or_ob +%type rangevar rangevar1 section section_or_ob rangeobj rangeobj1 /* END NEWCABLE */ /* OOP */ @@ -202,6 +202,15 @@ asgn: varname ROP expr ; /* OOP */ + +rangeobj: rangeobj1 + { code(sec_access_push); codesym((Symbol *)0);} + | section '.' rangeobj1 + ; +rangeobj1: RANGEOBJ {pushs($1); pushi(CHECK);} wholearray + { $$ = $3;} + ; + object: OBJECTVAR {pushi(OBJECTVAR);pushs($1); pushi(CHECK);} wholearray {$$ = $3; code(hoc_objectvar); spop(); codesym($1);} | OBJECTARG @@ -216,6 +225,10 @@ object: OBJECTVAR {pushi(OBJECTVAR);pushs($1); pushi(CHECK);} wholearray | HOCOBJFUNCTION begin '(' arglist ')' { $$ = $2; code(call); codesym($1); codei($4); code(hoc_known_type); codei(OBJECTVAR); pushi(OBJECTVAR);} + | rangeobj + { code(rangeobjevalmiddle); codesym(spop()); pushi(OBJECTVAR);} + | rangeobj '(' expr ')' + {TPD; code(rangeobjeval); codesym(spop()); pushi(OBJECTVAR);} ; ob: ob1 { spop(); } @@ -975,7 +988,7 @@ ckvar: VAR anyname: STRING|VAR|UNDEF|FUNCTION|PROCEDURE|FUN_BLTIN|SECTION|RANGEVAR |NRNPNTVAR|OBJECTVAR|TEMPLATE|OBFUNCTION|AUTO|AUTOOBJ|SECTIONREF |MECHANISM|BLTIN|STRFUNCTION|HOCOBJFUNCTION|ITERATOR|STRINGFUNC - |OBJECTFUNC + |OBJECTFUNC|RANGEOBJ ; %% /* end of grammar */ diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index e5731f55fb..c1e5411f03 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -503,6 +503,14 @@ if(NRN_ENABLE_PYTHON) SCRIPT_PATTERNS test/nmodl/test_kinetic.py ENVIRONMENT ${sonata_zero_gid_env} ${nrnpython_mpi_env} COMMAND ${modtests_launch_py} test/nmodl/test_kinetic.py) + nrn_add_test( + GROUP nmodl_tests + NAME test_random + REQUIRES ${modtests_preload_sanitizer} + SCRIPT_PATTERNS test/nmodl/test_random.py + ENVIRONMENT ${sonata_zero_gid_env} ${nrnpython_mpi_env} + COMMAND ${modtests_launch_py} test/nmodl/test_random.py) + nrn_add_test_group( CORENEURON NAME coreneuron_modtests @@ -634,6 +642,22 @@ if(NRN_ENABLE_PYTHON) ENVIRONMENT ${modtests_processor_env} ${nrnpython_mpi_env} COVERAGE_FILE=.coverage.coreneuron_test_ba_py COMMAND ${modtests_launch_py} test/coreneuron/test_ba.py) + nrn_add_test( + GROUP coreneuron_modtests + NAME test_nmodlrandom_py_${processor} + REQUIRES coreneuron ${processor} ${modtests_preload_sanitizer} + SCRIPT_PATTERNS test/coreneuron/test_nmodlrandom.py + ENVIRONMENT ${modtests_processor_env} ${nrnpython_mpi_env} + COVERAGE_FILE=.coverage.coreneuron_test_nmodlrandom_py + COMMAND ${modtests_launch_py} test/coreneuron/test_nmodlrandom.py) + nrn_add_test( + GROUP coreneuron_modtests + NAME test_nmodlrandom_syntax_py_${processor} + REQUIRES coreneuron ${processor} ${modtests_preload_sanitizer} + SCRIPT_PATTERNS test/coreneuron/test_nmodlrandom_syntax.py + ENVIRONMENT ${modtests_processor_env} ${nrnpython_mpi_env} + COVERAGE_FILE=.coverage.coreneuron_test_nmodlrandom_syntax_py + COMMAND ${modtests_launch_py} test/coreneuron/test_nmodlrandom_syntax.py) nrn_add_test( GROUP coreneuron_modtests NAME test_natrans_py_${processor} @@ -760,6 +784,7 @@ add_modlunit_test("${PROJECT_SOURCE_DIR}/test/pytest_coreneuron/unitstest.mod") add_modlunit_test("${PROJECT_SOURCE_DIR}/src/nrnoc/hh.mod") add_modlunit_test("${PROJECT_SOURCE_DIR}/src/nrnoc/stim.mod") add_modlunit_test("${PROJECT_SOURCE_DIR}/src/nrnoc/pattern.mod") +add_modlunit_test("${PROJECT_SOURCE_DIR}/test/hoctests/rand.mod") set_property(TEST modlunit_pattern PROPERTY WILL_FAIL ON) set_tests_properties(${TESTS} PROPERTIES ENVIRONMENT "${NRN_RUN_FROM_BUILD_DIR_ENV}") diff --git a/test/coreneuron/mod files/noisychan.mod b/test/coreneuron/mod files/noisychan.mod new file mode 100644 index 0000000000..780feea691 --- /dev/null +++ b/test/coreneuron/mod files/noisychan.mod @@ -0,0 +1,59 @@ +: The idea is that the voltage follows the change in e because of high +: conductance. So every negexp event causes voltage to pass threshold + +NEURON { + SUFFIX noisychan + NONSPECIFIC_CURRENT i + RANGE tau, invl + RANDOM ran +} + +UNITS { + (mA) = (milliamp) + (mV) = (millivolt) + (S) = (siemens) +} + +PARAMETER { + g = .1 (S/cm2) + tau = .1 (ms) + invl = 1 (ms) + emax = 50 (mV) +} + +ASSIGNED { + v (mV) + i (mA/cm2) + tnext (ms) +} + +STATE { + e (mV) +} + +INITIAL { + random_setseq(ran, 0) + e = -65(mV) + tnext = negexp() +} + +BEFORE STEP { + if (t > tnext) { + tnext = tnext + negexp() + e = emax + } +} + +BREAKPOINT { + SOLVE conduct METHOD cnexp + i = g*(v - e) +} + +DERIVATIVE conduct { + e' = (-65(mV) - e)/tau +} + +FUNCTION negexp()(ms) { + negexp = invl*random_negexp(ran) +} + diff --git a/test/coreneuron/mod files/watchrange.mod b/test/coreneuron/mod files/watchrange.mod index a710a67e20..dd051ac23a 100644 --- a/test/coreneuron/mod files/watchrange.mod +++ b/test/coreneuron/mod files/watchrange.mod @@ -106,6 +106,18 @@ VERBATIM ENDVERBATIM } +DESTRUCTOR { +VERBATIM + { + nrnran123_State** pv = (nrnran123_State**)(&_p_ran); + if (*pv) { + nrnran123_deletestream(*pv); + *pv = nullptr; + } + } +ENDVERBATIM +} + VERBATIM static void bbcore_write(double* z, int* d, int* zz, int* offset, _threadargsproto_) { if (d) { @@ -115,6 +127,11 @@ static void bbcore_write(double* z, int* d, int* zz, int* offset, _threadargspro nrnran123_getids3(*pv, di, di+1, di+2); nrnran123_getseq(*pv, di+3, &which); di[4] = (int)which; +#if NRNBBCORE + /* CoreNEURON does not call DESTRUCTOR so ... */ + nrnran123_deletestream(*pv); + *pv = nullptr; +#endif } *offset += 5; } diff --git a/test/coreneuron/mod files/watchrange2.mod b/test/coreneuron/mod files/watchrange2.mod new file mode 100644 index 0000000000..13100e45a4 --- /dev/null +++ b/test/coreneuron/mod files/watchrange2.mod @@ -0,0 +1,76 @@ +: for testing multiple WATCH statements activated as same time +: high, low, and mid regions watch a random uniform variable. +: The random variable ranges from 0 to 1 and changes at random times in +: the neighborhood of interval tick. + +NEURON { + THREADSAFE + POINT_PROCESS Bounce2 + RANGE r, result, n_high, n_low, n_mid, tick, x, t1 + RANDOM ran +} + +PARAMETER { + tick = 0.25 (ms) + LowThresh = 0.3 + HighThresh = 0.7 +} + +ASSIGNED { + x (1) + t1 (ms) + r (1) + n_high (1) + n_mid (1) + n_low (1) + result (1) +} + +DEFINE Low 1 +DEFINE Mid 2 +DEFINE High 3 +DEFINE Clock 4 + +INITIAL { + random_setseq(ran, 0) + n_high = 0 + n_mid = 0 + n_low = 0 + r = uniform() + t1 = t + x = 0 + net_send(0, Mid) + net_send(0, Clock) +} + +:AFTER SOLVE { +: result = t1*100/1(ms) + x +:} + +NET_RECEIVE(w) { + t1 = t + if (flag == Clock) { + r = uniform() + net_send(tick*(uniform() + .5), Clock) + } + if (flag == High) { + x = High + n_high = n_high + 1 + WATCH (r < LowThresh) Low + WATCH (r < HighThresh ) Mid + }else if (flag == Mid) { + x = Mid + n_mid = n_mid + 1 + WATCH (r < LowThresh) Low + WATCH (r > HighThresh) High + }else if (flag == Low) { + x = Low + n_low = n_low + 1 + WATCH (r > HighThresh) High + WATCH (r > LowThresh) Mid + } +} + +FUNCTION uniform() { + uniform = random_uniform(ran) +} diff --git a/test/coreneuron/test_nmodlrandom.py b/test/coreneuron/test_nmodlrandom.py new file mode 100644 index 0000000000..067f26979e --- /dev/null +++ b/test/coreneuron/test_nmodlrandom.py @@ -0,0 +1,169 @@ +# augmented to also checkpoint verify when RANDOM is permuted +from neuron.tests.utils.strtobool import strtobool +import os + +from neuron import h, coreneuron + +pc = h.ParallelContext() + + +class Cell: + def __init__(self, gid): + self.soma = h.Section(name="soma", cell=self) + self.soma.insert("noisychan") + if gid % 2 == 0: + # CoreNEURON permutation not the identity if cell topology not homogeneous + self.dend = h.Section(name="dend", cell=self) + self.dend.connect(self.soma(0.5)) + self.dend.insert("noisychan") + self.gid = gid + pc.set_gid2node(gid, pc.id()) + pc.cell(gid, h.NetCon(self.soma(0.5)._ref_v, None, sec=self.soma)) + + +def model(): + nslist = [h.NetStim() for _ in range(3)] + cells = [Cell(gid) for gid in range(10, 15)] + for gid, ns in enumerate(nslist): + ns.start = 0 + ns.number = 1e9 + ns.interval = 1 + ns.noise = 1 + pc.set_gid2node(gid, pc.id()) + pc.cell(gid, h.NetCon(ns, None)) + spiketime = h.Vector() + spikegid = h.Vector() + pc.spike_record(-1, spiketime, spikegid) + return nslist, spiketime, spikegid, cells + + +def pspike(m): + print("spike raster") + for i, t in enumerate(m[1]): + print("%.5f %g" % (t, m[2][i])) + + +def run(tstop, m): + pc.set_maxstep(10) + h.finitialize(-65) + pc.psolve(tstop) + + +def chk(std, m): + assert std[0].eq(m[1]) + assert std[1].eq(m[2]) + + +def test_embeded_run(): + m = model() + run(5, m) + std = [m[1].c(), m[2].c()] + pc.psolve(7) + std2 = [m[1].c(), m[2].c()] + + coreneuron.enable = True + coreneuron.verbose = 0 + coreneuron.gpu = bool(strtobool(os.environ.get("CORENRN_ENABLE_GPU", "false"))) + run(5, m) + chk(std, m) + + coreneuron.enable = False + pc.psolve(7) + chk(std2, m) + + +def sortspikes(spiketime, gidvec): + return sorted(zip(spiketime, gidvec)) + + +def test_chkpnt(): + import shutil, os, platform, subprocess + + # clear out the old if any exist + shutil.rmtree("coredat", ignore_errors=True) + + m = model() + + # std spikes from 0-5 and 5-10 + run(5, m) + std1 = [m[1].c(), m[2].c()] + m[1].resize(0) + m[2].resize(0) + pc.psolve(10) + std2 = [m[1].c(), m[2].c()] + pspike(m) + + # Files for coreneuron runs + h.finitialize(-65) + pc.nrncore_write("coredat") + + def runcn(tstop, perm, args): + exe = os.path.join(os.getcwd(), platform.machine(), "special-core") + common = [ + "-d", + "coredat", + "--voltage", + "1000", + "--verbose", + "0", + "--cell-permute", + str(perm), + ] + + gpu_run = bool(strtobool(os.environ.get("CORENRN_ENABLE_GPU", "false"))) + if gpu_run: + common.append("--gpu") + + cmd = [exe] + ["--tstop", "{:g}".format(tstop)] + common + args + print(" ".join(cmd)) + subprocess.run( + cmd, + check=True, + shell=False, + ) + + runcn(5, 1, ["--checkpoint", "coredat/chkpnt", "-o", "coredat"]) + runcn(10, 1, ["--restore", "coredat/chkpnt", "-o", "coredat/chkpnt"]) + cmp_spks(sortspikes(std2[0], std2[1]), "coredat") + + # cleanup + shutil.rmtree("coredat", ignore_errors=True) + + +def cmp_spks(spikes, dir): # modified from test_pointer.py + import os, subprocess, shutil + + # sorted nrn standard spikes into dir/out.spk + with open(os.path.join(dir, "temp"), "w") as f: + for spike in spikes: + f.write("{:.8g}\t{}\n".format(spike[0], int(spike[1]))) + + # sometimes roundoff to %.8g gives different sort. + def help(cmd, name_in, name_out): + # `cmd` is some generic utility, which does not need to have a + # sanitizer runtime pre-loaded. LD_PRELOAD=/path/to/libtsan.so can + # cause problems for *nix utilities, so drop it if it was present. + env = os.environ.copy() + try: + del env["LD_PRELOAD"] + except KeyError: + pass + subprocess.run( + [ + shutil.which(cmd), + os.path.join(dir, name_in), + os.path.join(dir, name_out), + ], + check=True, + env=env, + shell=False, + ) + + help("sortspike", "temp", "nrn.spk") + help("sortspike", "chkpnt/out.dat", "chkpnt/out.spk") + help("cmp", "chkpnt/out.spk", "nrn.spk") + + +if __name__ == "__main__": + test_embeded_run() + test_chkpnt() diff --git a/test/coreneuron/test_nmodlrandom_syntax.py b/test/coreneuron/test_nmodlrandom_syntax.py new file mode 100644 index 0000000000..1dbe14ec61 --- /dev/null +++ b/test/coreneuron/test_nmodlrandom_syntax.py @@ -0,0 +1,96 @@ +from neuron import h +import subprocess +from pathlib import Path + + +# default args generate accepted nmodl string +def modfile( + s0=":s0", + s1="RANDOM rv1, rv2", + s2=":s2", + s3=":s3", + s4="x1 = random_negexp(rv1, 1.0)", + s5=":s5", + s6="foo = random_negexp(rv1, 1.0)", +): + txt = """ +: temp.mod file with format elements to test for RANDOM syntax errors + +%s : 0 error if randomvar is mentioned + +NEURON { + SUFFIX temp + RANGE x1 + %s : 1 declare randomvars +} + +%s : 2 error if randomvar is mentioned + +ASSIGNED { + x1 +} + +BEFORE STEP { + %s : 3 error if assign or eval a randomvar + %s : 4 random_function accepted but ranvar must be first arg +} + +FUNCTION foo(arg) { + %s : 5 LOCAL ranvar makes it a double in this scope + %s : 6 random_function accepted but ranvar must be first arg + foo = arg +} +""" % ( + s0, + s1, + s2, + s3, + s4, + s5, + s6, + ) + return txt + + +def run(cmd): + result = subprocess.run(cmd, capture_output=True) + return result + + +def chk_nmodl(txt, program="nocmodl", rcode=False): + f = open("temp.mod", "w") + f.write(txt) + f.close() + result = run([program, "temp.mod"]) + ret = (result.returncode == 0) == rcode + if ret: + Path("temp.mod").unlink(missing_ok=True) + Path("temp.cpp").unlink(missing_ok=True) + else: + print("chk_nmodl ", program, " return code ", result.returncode) + print(txt) + print(result.stderr.decode()) + print(result.stdout.decode()) + return (result.returncode == 0) == rcode + + +def test_syntax(): + for program in ["nocmodl", "nmodl"]: + foo = False + assert chk_nmodl(modfile(), program, rcode=True) + assert chk_nmodl(modfile(s0="ASSIGNED{rv1}"), program) + foo = True if program == "nocmodl" else False + assert chk_nmodl(modfile(s0="LOCAL rv1"), program, rcode=foo) + assert chk_nmodl(modfile(s2="ASSIGNED{rv1}"), program) + assert chk_nmodl(modfile(s2="LOCAL rv1"), program) + assert chk_nmodl(modfile(s3="rv1 = 1"), program) + assert chk_nmodl(modfile(s3="x1 = rv1"), program) + assert chk_nmodl(modfile(s4="foo(rv1)"), program) + assert chk_nmodl(modfile(s4="random_negexp()"), program) + assert chk_nmodl(modfile(s4="random_negexp(1.0)"), program) + assert chk_nmodl(modfile(s5="LOCAL rv1"), program) + assert chk_nmodl(modfile(s4="random_negexp(rv1, rv2)"), program) + + +if __name__ == "__main__": + test_syntax() diff --git a/test/coreneuron/test_watchrange.py b/test/coreneuron/test_watchrange.py index bc7b255bde..d121e87a3c 100644 --- a/test/coreneuron/test_watchrange.py +++ b/test/coreneuron/test_watchrange.py @@ -1,5 +1,6 @@ # Basically want to test that net_move statement doesn't get # mixed up with other instances. +# Augmented to also test RANDOM (Bounce2) with coreneuron permutation. from neuron.tests.utils.strtobool import strtobool import os @@ -13,7 +14,7 @@ class Cell: - def __init__(self, gid): + def __init__(self, gid, bounce=0): self.soma = h.Section(name="soma", cell=self) if gid % 2 == 0: # CoreNEURON permutation not the identity if cell topology not homogeneous @@ -21,11 +22,13 @@ def __init__(self, gid): self.dend.connect(self.soma(0.5)) self.gid = gid pc.set_gid2node(gid, pc.id()) - self.r = h.Random() - self.r.Random123(gid, 0, 0) - self.syn = h.Bounce(self.soma(0.5)) + if bounce == 0: + self.syn = h.Bounce(self.soma(0.5)) + self.syn.noiseFromRandom123(gid, 0, 1) + else: + self.syn = h.Bounce2(self.soma(0.5)) + self.syn.ran.set_ids(gid, 0, 1) pc.cell(gid, h.NetCon(self.soma(0.5)._ref_v, None, sec=self.soma)) - self.syn.noiseFromRandom123(gid, 0, 1) self.t1vec = h.Vector() self.t1vec.record(self.syn._ref_t1, sec=self.soma) self.xvec = h.Vector() @@ -44,7 +47,7 @@ def result(self): ) -def test_watchrange(): +def watchrange(): from neuron import coreneuron coreneuron.enable = False @@ -145,16 +148,26 @@ def runassert(mode): for mode in [0, 1, 2]: runassert(mode) + # replace Bounce with Bounce2 (Uses RANDOM declaration) + pc.gid_clear() + cells = [Cell(gid, bounce=2) for gid in gids] + for mode in [0, 1, 2]: + runassert(mode) + coreneuron.enable = False # teardown pc.gid_clear() return stdlist, tvec +def test_watchrange(): + watchrange() + + if __name__ == "__main__": from neuron import gui - stdlist, tvec = test_watchrange() + stdlist, tvec = watchrange() g = h.Graph() print("n_high n_mid n_low") for i, result in enumerate(stdlist): diff --git a/test/external/CMakeLists.txt b/test/external/CMakeLists.txt index 6026122a62..7a9c8f3e63 100644 --- a/test/external/CMakeLists.txt +++ b/test/external/CMakeLists.txt @@ -54,7 +54,7 @@ if("channel-benchmark" IN_LIST NRN_ENABLE_MODEL_TESTS) FetchContent_Declare( channel-benchmark GIT_REPOSITORY https://github.com/bluebrain/nmodlbench.git - GIT_TAG 85b63eafb914f61632ff138e0d787af71e10e527 + GIT_TAG 2313db91599bcdd83e4291ab508d1e4474e87f25 SOURCE_DIR ${PROJECT_SOURCE_DIR}/external/tests/channel-benchmark) FetchContent_MakeAvailable(channel-benchmark) add_subdirectory(channel-benchmark) diff --git a/test/hoctests/rand.mod b/test/hoctests/rand.mod new file mode 100644 index 0000000000..a0147a2291 --- /dev/null +++ b/test/hoctests/rand.mod @@ -0,0 +1,43 @@ +NEURON { + SUFFIX rantst + RANGE x1, x2 + RANDOM ran1, ran2 +} + +ASSIGNED { + x1 x2 +} + +INITIAL { + random_setseq(ran1, 5) + random_setseq(ran2, 5) + x1 = random_uniform(ran1) +} + +BEFORE STEP { + x2 = random_normal(ran2) +} + +FUNCTION uniform0() { + uniform0 = random_uniform(ran1) +} + +FUNCTION uniform2(min, max) { + uniform2 = random_uniform(ran1, min, max) +} + +FUNCTION negexp0() { + negexp0 = random_negexp(ran1) +} + +FUNCTION negexp1(mean) { + negexp1 = random_negexp(ran1, mean) +} + +FUNCTION normal0() { + normal0 = random_normal(ran1) +} + +FUNCTION normal2(mean, std) { + normal2 = random_normal(ran1, mean, std) +} diff --git a/test/hoctests/rand_art.mod b/test/hoctests/rand_art.mod new file mode 100644 index 0000000000..bdeac87022 --- /dev/null +++ b/test/hoctests/rand_art.mod @@ -0,0 +1,50 @@ +NEURON { + ARTIFICIAL_CELL RanArt + RANGE x1, x2, mean + RANDOM ran1, ran2 +} + +PARAMETER { mean = .1 (ms) } + +ASSIGNED { + x1 x2 +} + +INITIAL { + random_setseq(ran1, 5) + random_setseq(ran2, 5) + x1 = random_uniform(ran1) + net_send(mean*negexp0(), 1) +} + +NET_RECEIVE(w){ + if (flag == 1) { + net_send(mean*negexp0(), 1) + net_event(t) + x2 = t + } +} + +FUNCTION uniform0() { + uniform0 = random_uniform(ran1) +} + +FUNCTION uniform2(min, max) { + uniform2 = random_uniform(ran1, min, max) +} + +FUNCTION negexp0() { + negexp0 = random_negexp(ran1) +} + +FUNCTION negexp1(mean) { + negexp1 = random_negexp(ran1, mean) +} + +FUNCTION normal0() { + normal0 = random_normal(ran1) +} + +FUNCTION normal2(mean, std) { + normal2 = random_normal(ran1, mean, std) +} diff --git a/test/hoctests/rand_pp.mod b/test/hoctests/rand_pp.mod new file mode 100644 index 0000000000..046f8bcf96 --- /dev/null +++ b/test/hoctests/rand_pp.mod @@ -0,0 +1,43 @@ +NEURON { + POINT_PROCESS RanPP + RANGE x1, x2 + RANDOM ran1, ran2 +} + +ASSIGNED { + x1 x2 +} + +INITIAL { + random_setseq(ran1, 5) + random_setseq(ran2, 5) + x1 = random_uniform(ran1) +} + +BEFORE STEP { + x2 = random_normal(ran2) +} + +FUNCTION uniform0() { + uniform0 = random_uniform(ran1) +} + +FUNCTION uniform2(min, max) { + uniform2 = random_uniform(ran1, min, max) +} + +FUNCTION negexp0() { + negexp0 = random_negexp(ran1) +} + +FUNCTION negexp1(mean) { + negexp1 = random_negexp(ran1, mean) +} + +FUNCTION normal0() { + normal0 = random_normal(ran1) +} + +FUNCTION normal2(mean, std) { + normal2 = random_normal(ran1, mean, std) +} diff --git a/test/hoctests/tests/test_random.hoc b/test/hoctests/tests/test_random.hoc new file mode 100644 index 0000000000..61f08a0dab --- /dev/null +++ b/test/hoctests/tests/test_random.hoc @@ -0,0 +1,79 @@ +load_file("expect_err.hoc") + +objref rt, z +z = new List("NMODLRandom") +proc assert() { + if ($1 == 0) { + hoc_execerror("assert", "") + } +} + +// ARTIFICIAL_CELL syntax tests (same as POINT_PROCESS wrt RANDOM) +rt = new RanArt() + +print rt.ran1 +assert(z.count == 0) + +assert(rt.ran1.set_seq(5).set_ids(7,8,9).get_seq() == 5) +assert(rt.ran1.get_ids().x[2] == 9) +assert(z.count == 0) + +objref r +r = rt.ran1 +r.set_seq(25) +assert(rt.ran1.get_seq() == 25) +assert(z.count == 1) +assert(r.get_ids().x[1] == 8) +objref r +assert(z.count == 0) + +create cable +access cable +nseg=3 +insert rantst +finitialize(-70) + +r = cable.ran1_rantst +assert(cable.ran1_rantst.set_seq(10).get_seq() == 10) +assert(cable.ran1_rantst(.5).get_ids().eq(r.get_ids())) + +for (x,0) { + r = cable.ran1_rantst(x) + assert(ran1_rantst(x).get_ids().eq(r.get_ids())) +} +objref r +assert(z.count == 0) + +// test syntax in a cell type +begintemplate Cell +public cable, rt +create cable +objref rt +proc init() { + create cable + cable { + nseg = 3 + insert rantst + rt = new RanPP(.1) + } +} +endtemplate Cell + +objref cell +cell = new Cell() +cell.cable print ran2_rantst +cell.cable print ran2_rantst(.1) +expect_err("cell.cable print ran2_rantst(.1, .2)") +print cell.cable.ran2_rantst +print cell.cable.ran2_rantst(.1) +expect_err("print cell.cable.ran2_rantst(.1, .2)") +cell.cable.ran2_rantst(.1).set_seq(10).set_ids(8,9,10) +cell.cable.ran2_rantst(.5).set_seq(11).set_ids(1,2,3) +assert( cell.cable.ran2_rantst(.1).get_seq() == 10) + +objref rp +rp = cell.rt +objref cell // but rp is still a reference to unlocated point process. +expect_err("rp.ran1") + + diff --git a/test/hoctests/tests/test_random.json b/test/hoctests/tests/test_random.json new file mode 100644 index 0000000000..4ad254d50c --- /dev/null +++ b/test/hoctests/tests/test_random.json @@ -0,0 +1,175 @@ +{ + "results": [ + 0.17893146859041148, + 2.5858652899862125, + 0.5900787127623134, + 13.879620017511858, + 2.5285665260297696, + -3.8255380765067546, + 0.6195943787182694, + 0.8640101598892337, + 0.8705624195582786, + 0.5816326112923024, + 0.6824968245405592, + 12.935409682585064, + 5.5629342829617965, + 15.621557479479236, + 0.00018058368000653952, + 0.7241277485796884, + 0.5950188712228511, + 1.4151762421390146, + -1.0273020350279256, + 11.370265127306276, + 1.5457504368518737, + 7.40405624730305, + 0.7815047526309489, + 0.9767377004081529, + 0.9199947437923414, + 1.770347276557777, + 1.1348050752784042, + 10.21118748695329, + 0.8290698041860384, + 3.2620640633572284, + 0.7935756776031163, + 0.6215100738635496, + 13.0, + 1.0, + 9.0, + 1.0, + 9.0, + 1.0, + 9.0, + 1.0 + ], + "bbsavestate 1 to 2": [ + [ + -0.7247871714935129, + 0.7219008192457033, + 0.5391448928750727, + -0.22338622131375607, + -0.5209245315875433, + -0.20192660100822468, + -0.40971682383225255, + 0.40125116545729983, + -1.018708027424587, + -0.6435348201983898, + -1.110505571561422, + -0.5456431266293069, + -0.6431811174215588, + -0.2050222092130618, + -0.41871095517167045, + 2.1327022497045847, + 1.2475561990837372, + -0.8148916287356835, + 1.3637834172361567, + 0.44911001868255424, + -0.8038191925528606, + -1.0182609292625484, + 0.814562018050233, + -0.2001845683092376, + -0.7036837840067589, + -0.1749096487126645, + -0.1384874349183572, + -0.27802269131343627, + 0.5024084209351181, + 0.2857886341240531, + 1.1578619121493123, + 0.8320672365682753, + 0.8388016796407406, + -2.204403101188663, + -0.41677220453431346, + -0.3493109618344672, + 0.6394206556529147, + 0.433099968696616, + -1.6289203897217996, + -0.5961253470352493, + -0.028642516485346325 + ], + [ + -0.007237015674221396, + 1.0485038145975802, + 0.5666028312087761, + 1.5478221474628917, + -1.9262214185465736, + -1.1342009341579276, + 0.7839089952205353, + -2.1453640858801464, + 0.5459020601695509, + 0.5702794219336828, + -2.7987171853712787, + -0.5441341689538746, + 0.16498316711023667, + 0.5332911939974267, + 1.5636051898318526, + -0.6390453305292507, + -0.14494230468386957, + -0.7174754220603327, + 0.0539835666341803, + 0.9773045846743287, + 0.36026791233149735, + -1.1293229505809468, + -0.6829744883194595, + -0.6966328933631747, + -1.5001099642041043, + -0.07403611736759062, + -0.8287345152039473, + -1.7481478211594699, + -0.0019827309900973702, + -0.2831062597817745, + -1.3860153441895837, + -0.7725712206688502, + 0.9644474198597336, + 0.038712287163401804, + -0.49387615024354914, + -0.5149700079858756, + -0.06700819948825379, + -1.1356059525253281, + -0.4846497787175894, + -0.05689192371120626, + 0.31197352531127626 + ], + [ + 0.952070199416881, + 0.952070199416881, + 0.952070199416881, + 0.952070199416881, + 0.952070199416881, + 0.952070199416881, + 0.952070199416881, + 1.1554634134011021, + 1.1554634134011021, + 1.1554634134011021, + 1.1554634134011021, + 1.1554634134011021, + 1.1554634134011021, + 1.1554634134011021, + 1.1554634134011021, + 1.1554634134011021, + 1.1554634134011021, + 1.1554634134011021, + 1.1554634134011021, + 1.4540666926704917, + 1.4540666926704917, + 1.4540666926704917, + 1.4540666926704917, + 1.4540666926704917, + 1.4540666926704917, + 1.4540666926704917, + 1.633567197773812, + 1.633567197773812, + 1.633567197773812, + 1.633567197773812, + 1.7246756477971599, + 1.7246756477971599, + 1.7246756477971599, + 1.7246756477971599, + 1.7246756477971599, + 1.7246756477971599, + 1.7246756477971599, + 1.8967211458457505, + 1.8967211458457505, + 1.8967211458457505, + 1.8967211458457505 + ] + ] +} diff --git a/test/hoctests/tests/test_random.py b/test/hoctests/tests/test_random.py new file mode 100644 index 0000000000..90d60389dd --- /dev/null +++ b/test/hoctests/tests/test_random.py @@ -0,0 +1,178 @@ +from neuron import h +from neuron.expect_hocerr import expect_err, set_quiet +from neuron.tests.utils.checkresult import Chk +import os +import math + +pc = h.ParallelContext() + +dir_path = os.path.dirname(os.path.realpath(__file__)) +chk = Chk(os.path.join(dir_path, "test_random.json")) + +set_quiet(False) + +z = h.List("NMODLRandom") + +# ARTIFICIAL_CELL syntax tests +rt = h.RanArt() + +print(rt.ran1) +assert z.count() == 0 +assert rt.ran1.set_seq(5).set_ids(7, 8, 9).get_seq() == 5 +assert rt.ran1.get_ids().x[2] == 9 +assert z.count() == 0 + +x = rt.ran1 +x.set_seq(25) +assert rt.ran1.get_seq() == 25 +assert z.count() == 1 +assert x.get_ids().x[1] == 8 + +y = rt.ran1.set_seq +y(50) +assert x.get_seq() == 50 + +del x, y +assert z.count() == 0 + +x = rt.ran1 +expect_err("rt.ran1 = rt.ran2") # cannot assign +del rt +expect_err("x.get_seq()") +del x +assert z.count() == 0 + +# density mechanism tests +cable = h.Section(name="cable") +cable.nseg = 3 +cable.insert("rantst") +nr = [(seg.x, seg.rantst.ran1) for seg in cable] +for r in nr: + assert r[1].get_ids().eq(cable(r[0]).ran1_rantst.get_ids()) +del nr, r +assert z.count() == 0 + +rt = h.RanArt() +r1 = rt.ran1 +# wrap around at end +assert r1.set_seq(2**34 - 1).get_seq() == (2**34 - 1) +r1.uniform() +assert r1.get_seq() == 0 +assert r1.set_seq(2**34 + 1).get_seq() == 1 +assert r1.set_seq(-10).get_seq() == 0 # all neg setseq to 0 +assert r1.set_seq(2**40 - 1).get_seq() == 17179869183.0 # up to 2**40 wrap +assert r1.set_seq(2**40 + 1).get_seq() == 0 # all above 2**40 setseq to 0 + +# negexp(mean) has proper scale +r1.set_seq(0) +x = rt.negexp0() +r1.set_seq(0) +assert math.isclose(rt.negexp1(5), 5 * x) + +r1 = h.NMODLRandom() +expect_err("r1.uniform()") +del r1 + +# test the mod random_... functions +rt = h.RanArt() +cable = h.Section(name="cable") +cable.nseg = 3 +cable.insert("rantst") +mlist = [rt] +mlist.extend(seg.rantst for seg in cable) +for i, m in enumerate(mlist): + m.ran1.set_ids(i, 1, 0).set_seq(0) + m.ran2.set_ids(i, 2, 0).set_seq(0) + +results = [] +for m in mlist: + results.extend([m.uniform0(), m.negexp0(), m.normal0()]) + results.extend( + [ + m.uniform2(10, 20), + m.negexp1(5), + m.normal2(5, 8), + ] + ) + results.extend([m.ran1.uniform(), m.ran2.uniform()]) +for m in mlist: + results.extend([m.ran1.get_seq(), m.ran2.get_seq()]) + +assert z.count() == 0 + +chk("results", results, tol=1e-12) + +# test access to hoc cell class +h( + """ +begintemplate Cell +public cable, rt +create cable +objref rpp, rart +proc init() { + create cable + cable { + nseg = 3 + insert rantst + rpp = new RanPP(0.1) + } + rart = new RanArt() + rart.mean = 0.1 +} +endtemplate Cell +""" +) +cell = h.Cell() +cell.cable(0.1).ran2_rantst.set_seq(10).set_ids(8, 9, 10) +cell.cable(0.5).ran2_rantst.set_seq(11).set_ids(1, 2, 3) +assert cell.cable(0.1).rantst.ran2.get_seq() == 10 + +rp = cell.rpp +del cell +expect_err("rp.ran2") + + +def set_gids(cells): + gid = 0 + for cell in cells: + for port in [cell.cable(0.5)._ref_v, cell.rart]: + gid += 1 + pc.set_gid2node(gid, pc.id()) + pc.cell(gid, h.NetCon(port, None, sec=cell.cable)) + + +def test_bbsavestate(): + cell = h.Cell() + set_gids([cell]) ## BBSaveState requires that real cells have gids + mechs = [cell.cable(0.5).rantst, cell.rpp, cell.rart] + rec = [] + for i, m in enumerate(mechs): + m.ran2.set_ids(i, 2, 3) + rec.append(h.Vector().record(m._ref_x2, sec=cell.cable)) + + def run(tstop, init=True): + if init: + pc.set_maxstep(10) + h.finitialize(-65) + else: + h.frecord_init() + pc.psolve(tstop) + + # rerunning in two parts gives same result. + run(1) + bbss = h.BBSaveState() + bbss.save("bbss_random.txt") + + # second half is our standard + run(2, False) + chk("bbsavestate 1 to 2", [v.to_python() for v in rec], tol=1e-12) + + # savestate restore and redo second half. + bbss.restore("bbss_random.txt") + run(2, False) + chk("bbsavestate 1 to 2", [v.to_python() for v in rec], tol=1e-12) + + +test_bbsavestate() + +chk.save() diff --git a/test/nmodl/test_random.py b/test/nmodl/test_random.py new file mode 100644 index 0000000000..75a6ca093a --- /dev/null +++ b/test/nmodl/test_random.py @@ -0,0 +1,129 @@ +from math import isclose + +from neuron import h + +pc = h.ParallelContext() + + +def model(): + pc.gid_clear() + for s in h.allsec(): + h.delete_section(sec=s) + s = h.Section() + s.L = 10 + s.diam = 10 + s.insert("hh") + ic = h.IClamp(s(0.5)) + ic.delay = 0.1 + ic.dur = 0.1 + ic.amp = 0.5 * 0 + syn = h.ExpSyn(s(0.5)) + nc = h.NetCon(None, syn) + nc.weight[0] = 0.001 + return {"s": s, "ic": ic, "syn": syn, "nc": nc} + + +def test_netstim_noise(): + cells = {gid: (h.NetStim()) for gid in range(pc.id(), 5, pc.nhost())} + for gid, cell in cells.items(): + pc.set_gid2node(gid, pc.id()) + pc.cell(gid, h.NetCon(cell, None)) + + cell.interval = gid + 1 + cell.number = 100 + cell.start = 0 + cell.noise = 1 + + # Initialize RANDOM variable ids and initial sequence. + cell.ranvar.set_ids(gid, 2, 3).set_seq(0) + + spiketime = h.Vector() + spikegid = h.Vector() + pc.spike_record(-1, spiketime, spikegid) + + pc.set_maxstep(10) + tstop = 5 + + h.finitialize() + pc.psolve(tstop) + + spiketime_result = spiketime.c() + spikegid_result = spikegid.c() + + spikegid_ref = [ + 1.0, + 0.0, + 2.0, + 0.0, + 0.0, + 0.0, + 4.0, + 3.0, + 0.0, + 4.0, + 1.0, + 0.0, + 0.0, + 2.0, + 2.0, + ] + spiketime_ref = [ + 0.038647213491710304, + 0.08268113588796304, + 0.5931985927363619, + 0.7687313066056471, + 0.867367543646173, + 1.1822988033563793, + 1.3476448598895432, + 1.748395215773899, + 1.9382702939333631, + 2.3381219031177376, + 2.9858151911753863, + 3.3721447007688603, + 3.4714402733585277, + 4.130940076465841, + 4.406639959683753, + ] + + # check if gid and spike time matches + for i in range(int(spiketime_result.size())): + assert spikegid_ref[i] == spikegid_result[i] and isclose( + spiketime_ref[i], spiketime_result[i] + ) + + +def test_random(): + pc.gid_clear() + ncell = 10 + cells = [] + gids = range(pc.id(), ncell, pc.nhost()) + rng = [] + for gid in gids: + pc.set_gid2node(gid, pc.id()) + cell = h.NetStim() + + cell.ranvar.set_ids(gid, 2, 3).set_seq(0) + + r = cell.erand() + rng.append(r) + + rng_ref = [ + 0.08268113588796304, + 0.019323606745855152, + 0.19773286424545394, + 0.4370988039434747, + 0.26952897197790865, + 0.27785183823847076, + 1.4834024918566038, + 1.1439786359830195, + 0.13094521398166833, + 0.3743746759204156, + ] + + for x, y in zip(rng, rng_ref): + assert isclose(x, y) + + +if __name__ == "__main__": + test_netstim_noise() + test_random() diff --git a/test/rxd/test_currents.py b/test/rxd/test_currents.py index ae69b82804..32c052fbbe 100644 --- a/test/rxd/test_currents.py +++ b/test/rxd/test_currents.py @@ -1,5 +1,10 @@ import pytest from .testutils import compare_data, tol +from platform import platform + + +def applearm(): + return "macOS-" in platform() and "-arm64-" in platform() @pytest.fixture @@ -80,7 +85,7 @@ def test_currents_cvode(model_pump): h.continuerun(10) if not save_path: max_err = compare_data(data) - assert max_err < tol + assert max_err < (1e-8 if applearm() else tol) def test_currents_stucture_change(model_pump): diff --git a/test/rxd/test_pure_diffusion.py b/test/rxd/test_pure_diffusion.py index 1e627d67ee..423a21e366 100644 --- a/test/rxd/test_pure_diffusion.py +++ b/test/rxd/test_pure_diffusion.py @@ -1,4 +1,9 @@ from .testutils import compare_data, tol +from platform import platform + + +def applearm(): + return "macOS-" in platform() and "-arm64-" in platform() def test_pure_diffusion(neuron_instance): @@ -52,4 +57,4 @@ def test_pure_diffusion_cvode(neuron_instance): h.continuerun(t) if not save_path: max_err = compare_data(data) - assert max_err < tol + assert max_err < (2e-10 if applearm() else tol) From 0ced58524b26f7e37a93e89acf858d6e936bfd68 Mon Sep 17 00:00:00 2001 From: Luc Grosheintz Date: Thu, 8 Feb 2024 20:17:29 +0100 Subject: [PATCH 13/15] Call `nrnivmodl` directly during building (#2716) * Prior to this commit `neurondemo` as called to generate `special`. During development there can be bugs that cause `neurondemo` to crash without building `special`. This causes the build to fail, which makes it harder to debug the failure. * This commit calls `nrnivmodl` directly to generate `special`. --- CMakeLists.txt | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 44983a771b..c35b8507fd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -827,10 +827,11 @@ if(NRN_ENABLE_TESTS) "${neurondemo_prefix}/${CMAKE_SHARED_LIBRARY_PREFIX}nrnmech${CMAKE_SHARED_LIBRARY_SUFFIX}") add_custom_command( OUTPUT ${neurondemo_files} - COMMAND - ${CMAKE_COMMAND} -E env ${NRN_RUN_FROM_BUILD_DIR_ENV} ${NRN_SANITIZER_ENABLE_ENVIRONMENT} - "${PROJECT_BINARY_DIR}/bin/neurondemo" "-nobanner" "-nogui" "-c" "quit()" + COMMAND ${CMAKE_COMMAND} -E env ${NRN_RUN_FROM_BUILD_DIR_ENV} + ${NRN_SANITIZER_ENABLE_ENVIRONMENT} "${PROJECT_BINARY_DIR}/bin/nrnivmodl" + WORKING_DIRECTORY ${PROJECT_BINARY_DIR}/share/nrn/demo/release VERBATIM) + add_custom_target(generate-neurondemo-mechanism-library ALL DEPENDS ${neurondemo_files} nrniv copy_share_lib_to_build) # Initialize the submodule *before* including the test/CMakeLists.txt that uses it. This ensures From 1de44de0bf9ac4ea9e0ac55f990f9388730ef56f Mon Sep 17 00:00:00 2001 From: Luc Grosheintz Date: Thu, 8 Feb 2024 20:18:58 +0100 Subject: [PATCH 14/15] Accessors for cache pointers (#2712) In NMODL we have a struct: struct ionic_Instance { double* ena{}; double* ina{}; double* v_unused{}; double* g_unused{}; const double* const* ion_ina{}; double* const* ion_ena{}; ionic_Store* global{&ionic_global}; }; The pointers are views into contiguous arrays of doubles. While we can the address for parameters like `ena`, we can't get the address of "double pointer" variables, with the previous existing API. * fixup: missing template arguments. --- src/neuron/cache/mechanism_range.hpp | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/neuron/cache/mechanism_range.hpp b/src/neuron/cache/mechanism_range.hpp index 4a133b0b44..751639c6dc 100644 --- a/src/neuron/cache/mechanism_range.hpp +++ b/src/neuron/cache/mechanism_range.hpp @@ -87,6 +87,11 @@ struct MechanismRange { return std::next(m_data_ptrs[variable], array_size * (m_offset + instance)); } + template + [[nodiscard]] double* data_array_ptr() { + return data_array(0); + } + /** * @brief Get a RANGE variable value. * @tparam variable The index of the RANGE variable in the mechanism. @@ -99,6 +104,11 @@ struct MechanismRange { return *data_array(instance); } + template + [[nodiscard]] double* fpfield_ptr() { + return data_array(0); + } + /** * @brief Get a RANGE variable value. * @param instance Which mechanism instance to access inside this MechanismRange. @@ -124,6 +134,12 @@ struct MechanismRange { return m_pdata_ptrs[variable][m_offset + instance]; } + template + [[nodiscard]] double* const* dptr_field_ptr() { + static_assert(variable < NumDatumFields); + return m_pdata_ptrs[variable] + m_offset; + } + protected: /** * @brief Pointer to a range of pointers to the start of RANGE variable storage. From f7aba0abe56224f1f84d0461833555a597270a5c Mon Sep 17 00:00:00 2001 From: Erik Heeren Date: Thu, 8 Feb 2024 20:21:45 +0100 Subject: [PATCH 15/15] Issue #2678: add container build job on gitlab (#2688) * Issue #2678: add container build job on gitlab * Update documentation for the automated docker build * A manual action to build the wheels container for x86_64 * Make specifying the tag through a variable mandatory As a consequence, only run on manually-triggered pipelines * document which variables need to be set when starting the podman machine on the m1 runner --- .gitlab-ci.yml | 58 ++++++++++++++++++++++++++++++++++- docs/install/python_wheels.md | 20 ++++++++++-- 2 files changed, 75 insertions(+), 3 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index bd944ee463..fa9dc5cf0e 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -7,7 +7,6 @@ include: - project: hpc/gitlab-upload-logs file: enable-upload.yml - # see https://gitlab.com/gitlab-org/gitlab/-/issues/263401 for why we specify the flags like this now # 130 characters seems to be the point at which jobs refuse to run .matrix: @@ -263,3 +262,60 @@ test:neuron:nmodl:nvhpc:omp:legacy: test:neuron:nmodl:nvhpc:omp: extends: [.test_neuron, .gpu_node] needs: ["build:neuron:nmodl:nvhpc:omp"] + + +# Container building +mac_m1_container_build: + stage: .pre + tags: + - macos-arm64 + script: + - if [ -z "${ARM64_IMAGE_TAG}" ]; then + - echo "Please set the ARM64_IMAGE_TAG variable" + - exit 1 + - fi + - cd packaging/python + - echo "Replacing symlinks with their targets to keep podman happy" + - find . -type l -exec cp $(realpath {}) ./TEMP \; -exec rm {} \; -exec mv TEMP {} \; + - ls -l + - export BUILDAH_FORMAT=docker # enables ONBUILD instructions which are not OCI compatible + - machine_status=$(podman machine inspect | awk '/State/ {print $2}' | tr -d '",') + # If you start the machine yourself, make sure BUILDAH_FORMAT and the http proxy variables are set in the shell before doing so! + - if [[ "${machine_status}" != "running" ]]; then + - echo "Machine is in ${machine_status} status - starting" + - podman machine start + - fi + - podman build -t neuronsimulator/neuron_wheel:${ARM64_IMAGE_TAG} --build-arg MANYLINUX_IMAGE=manylinux2014_aarch64 -f Dockerfile . + - podman login -u ${DOCKER_HUB_USER} -p ${DOCKER_HUB_AUTH_TOKEN} docker.io + - podman push neuronsimulator/neuron_wheel:${ARM64_IMAGE_TAG} + - podman rmi localhost/neuronsimulator/neuron_wheel:${ARM64_IMAGE_TAG} + rules: + - if: $CI_PIPELINE_SOURCE == "web" + when: manual + +x86_64_container_build: + stage: .pre + image: + name: quay.io/buildah/stable + entrypoint: [""] + variables: + KUBERNETES_CPU_LIMIT: 4 + KUBERNETES_CPU_REQUEST: 2 + KUBERNETES_MEMORY_LIMIT: 8Gi + KUBERNETES_MEMORY_REQUEST: 4Gi + tags: + - kubernetes + rules: + - if: $CI_PIPELINE_SOURCE == "web" + when: manual + script: + - if [ -z "${X86_IMAGE_TAG}" ]; then + - echo "Please set the X86_IMAGE_TAG variable" + - exit 1 + - fi + - export STORAGE_DRIVER=vfs # allows to build inside containers without additional mounts + - export BUILDAH_FORMAT=docker # enables ONBUILD instructions which are not OCI compatible + - cd packaging/python + - buildah bud --iidfile image_id -t neuronsimulator/neuron_wheel:${X86_IMAGE_TAG} -f Dockerfile . + - buildah login -u ${DOCKER_HUB_USER} -p ${DOCKER_HUB_AUTH_TOKEN} docker.io + - buildah push $(