From 758d15c8e39e5dc18d8f5ba91ff027ad232dc0e1 Mon Sep 17 00:00:00 2001 From: nrnhines Date: Thu, 8 Feb 2024 05:57:58 -0500 Subject: [PATCH] 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)