diff --git a/README.rst b/README.rst index b5322320..6d31205b 100644 --- a/README.rst +++ b/README.rst @@ -12,9 +12,9 @@ A universal workflow for creation, validation, and generalization of detailed ne Introduction --------- -The biophysically detailed electrical neuron model (e-model) is one of the central tools in neuroscience. Here we present a demo of the single neuron e-model creation, validation, and generalization described in `A universal workflow for creation, validation, and generalization of detailed neuronal models (Reva, Rossert et al) `_. +The biophysically detailed electrical neuron model (e-model) is one of the central tools in neuroscience. Here we present a demo of the single neuron e-model creation, validation, and generalization described in `A universal workflow for creation, validation, and generalization of detailed neuronal models (Reva, Rossert et al) `_. -This demo is built on the example of L5PC of the SSCx in juvenile rat. +This demo is built on the example of L5PC of the SSCx in juvenile rat. Structure of the codebase ~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -34,7 +34,7 @@ When you using methods or code from this repository for your research, we ask yo A universal workflow for creation, validation, and generalization of detailed neuronal models" (Reva, Rossert et al). -.. code-block:: +.. code-block:: @article{reva2023universal, title={A universal workflow for creation, validation, and generalization of detailed neuronal models}, @@ -86,7 +86,7 @@ The visualization of the bAP/EPSP validations can be found in `validation.ipynb The morphologies for these validations are located in the `input/morphologies `_ folder. -To run bAP/EPSP validations use:: +To run bAP/EPSP validations use:: python main.py att_conf.json @@ -98,10 +98,10 @@ Somatic validations are located in the `somatic_validation `_ needs to be installed and the mechanisms need to be compiled with the following command before running the notebooks:: +`somatic-val-requirements.txt `_ needs to be installed and the mechanisms need to be compiled with the following command before running the notebooks:: + + nrnivmodl mechanisms - nrnivmodl mechanisms - First, e-features for the validations have to be extracted from the chosen patch clamp protocol. To extract e-features use `feature-extraction.ipynb `_, the results of this extraction can be found in the `somatic_validation/L5TPC `_ folder. To run and visualize results of the somatic validation run `somatic-validation.ipynb `_. 4. Generalization @@ -136,6 +136,7 @@ Instead of manually compiling the mechanisms for each step of the pipeline, the Requirements ------------ +For the exact reproducibility of paper results, we recommend `NEURON>=7.8.0,<8.0.0` and repository code release version `1.0.1 `. Users with other versions may encounter some small discrepancies in results. The `requirements.txt `_ at the main directory should be used for all steps except for the somatic validations. Install `somatic-val-requirements.txt `_ before running the somatic validation notebooks or tests. @@ -147,7 +148,7 @@ This project/research was supported by funding to the Blue Brain Project, a rese License ------- -This work is licensed under `Creative Commons (CC BY) 4.0 `_ +This work is licensed under `Creative Commons (CC BY) 4.0 `_ For MOD files for which the original source is available on ModelDB, any specific licenses on mentioned on ModelDB, or the generic License of ModelDB apply. diff --git a/model_management/mm_run_minimal/emodels_dir/singlecell-optimization/mechanisms/StochKv3.mod b/model_management/mm_run_minimal/emodels_dir/singlecell-optimization/mechanisms/StochKv3.mod index 92d955e9..06203e15 100644 --- a/model_management/mm_run_minimal/emodels_dir/singlecell-optimization/mechanisms/StochKv3.mod +++ b/model_management/mm_run_minimal/emodels_dir/singlecell-optimization/mechanisms/StochKv3.mod @@ -1,4 +1,4 @@ -TITLE StochKv2.mod +TITLE StochKv3.mod COMMENT ---------------------------------------------------------------- @@ -28,11 +28,11 @@ NEURON { USEION k READ ek WRITE ik RANGE N, eta, gk, gamma, deterministic, gkbar, ik RANGE N0, N1, n0_n1, n1_n0 - GLOBAL ninf, linf, ltau, ntau, an, bn, al, bl - GLOBAL P_an, P_bn, P_al, P_bl + RANGE ninf, linf, ltau, ntau, an, bn, al, bl + RANGE P_an, P_bn, P_al, P_bl GLOBAL vmin, vmax - :BBCOREPOINTER rng - POINTER rng + BBCOREPOINTER rng + :POINTER rng } UNITS { @@ -99,6 +99,7 @@ for comparison with Pr to decide whether to activate the synapse or not ENDCOMMENT VERBATIM +#ifndef NRN_VERSION_GTEQ_8_2_0 #include "nrnran123.h" extern int cvode_active_; @@ -106,8 +107,14 @@ extern int cvode_active_; #include #include +#ifndef CORENEURON_BUILD double nrn_random_pick(void* r); void* nrn_random_arg(int argpos); +#endif +#define RANDCAST +#else +#define RANDCAST (Rand*) +#endif ENDVERBATIM : ---------------------------------------------------------------- @@ -117,6 +124,10 @@ INITIAL { if (cvode_active_ && !deterministic) { hoc_execerror("StochKv2 with deterministic=0", "cannot be used with cvode"); } + + if( usingR123 ) { + nrnran123_setseq((nrnran123_State*)_p_rng, 0, 0); + } ENDVERBATIM eta = (gkbar / gamma) * (10000) @@ -248,9 +259,12 @@ PROCEDURE trates(v (mV)) { FUNCTION strap(x) { if (x < 0) { strap = 0 -VERBATIM - fprintf (stderr,"skv.mod:strap: negative state"); -ENDVERBATIM +: This function gets executed in the DERIVATIVE and BREAKPOINT blocks +: which should get vectorized and fprintf doesn't allow this +: Also Intel classic compilers generate invalid code when #pragma omp simd is used +:VERBATIM +: fprintf (stderr,"skv.mod:strap: negative state"); +:ENDVERBATIM } else { strap = x } @@ -273,11 +287,12 @@ PROCEDURE setRNG() { VERBATIM // For compatibility, allow for either MCellRan4 or Random123. Distinguish by the arg types // Object => MCellRan4, seeds (double) => Random123 -#if !NRNBBCORE +#ifndef CORENEURON_BUILD usingR123 = 0; if( ifarg(1) && hoc_is_double_arg(1) ) { nrnran123_State** pv = (nrnran123_State**)(&_p_rng); uint32_t a2 = 0; + uint32_t a3 = 0; if (*pv) { nrnran123_deletestream(*pv); @@ -286,7 +301,10 @@ VERBATIM if (ifarg(2)) { a2 = (uint32_t)*getarg(2); } - *pv = nrnran123_newstream((uint32_t)*getarg(1), a2); + if (ifarg(3)) { + a3 = (uint32_t)*getarg(3); + } + *pv = nrnran123_newstream3((uint32_t)*getarg(1), a2, a3); usingR123 = 1; } else if( ifarg(1) ) { void** pv = (void**)(&_p_rng); @@ -302,46 +320,56 @@ ENDVERBATIM FUNCTION urand() { VERBATIM - double value; + double value = 0.0; if( usingR123 ) { value = nrnran123_dblpick((nrnran123_State*)_p_rng); } else if (_p_rng) { - value = nrn_random_pick(_p_rng); +#ifndef CORENEURON_BUILD + value = nrn_random_pick(RANDCAST _p_rng); +#endif } else { - value = 0.5; + // see BBPBGLIB-972 + value = 0.0; } _lurand = value; ENDVERBATIM } VERBATIM -/* static void bbcore_write(double* x, int* d, int* xx, int* offset, _threadargsproto_) { if (d) { uint32_t* di = ((uint32_t*)d) + *offset; // temporary just enough to see how much space is being used if (!_p_rng) { - di[0] = 0; di[1] = 0; + di[0] = 0; di[1] = 0, di[2] = 0; }else{ nrnran123_State** pv = (nrnran123_State**)(&_p_rng); - nrnran123_getids(*pv, di, di+1); + nrnran123_getids3(*pv, di, di+1, di+2); + // write stream sequence + char which; + nrnran123_getseq(*pv, di+3, &which); + di[4] = (int)which; } -//printf("StochKv2.mod %p: bbcore_write offset=%d %d %d\n", _p, *offset, d?di[0]:-1, d?di[1]:-1); + //printf("StochKv3.mod %p: bbcore_write offset=%d %d %d\n", _p, *offset, d?di[0]:-1, d?di[1]:-1); } - *offset += 2; + *offset += 5; } static void bbcore_read(double* x, int* d, int* xx, int* offset, _threadargsproto_) { - assert(!_p_rng); uint32_t* di = ((uint32_t*)d) + *offset; - if (di[0] != 0 || di[1] != 0) + if (di[0] != 0 || di[1] != 0|| di[2] != 0) { nrnran123_State** pv = (nrnran123_State**)(&_p_rng); - *pv = nrnran123_newstream(di[0], di[1]); +#if !NRNBBCORE + if(*pv) { + nrnran123_deletestream(*pv); + } +#endif + *pv = nrnran123_newstream3(di[0], di[1], di[2]); + nrnran123_setseq(*pv, di[3], (char)di[4]); } -//printf("StochKv2.mod %p: bbcore_read offset=%d %d %d\n", _p, *offset, di[0], di[1]); - *offset += 2; + //printf("StochKv3.mod %p: bbcore_read offset=%d %d %d\n", _p, *offset, di[0], di[1]); + *offset += 5; } -*/ ENDVERBATIM : Returns random numbers drawn from a binomial distribution @@ -471,13 +499,16 @@ VERBATIM FUNCTION bbsavestate() { bbsavestate = 0 VERBATIM -#ifdef ENABLE_SAVE_STATE + #ifndef CORENEURON_BUILD // TODO: since N0,N1 are no longer state variables, they will need to be written using this callback // provided that it is the version that supports multivalue writing /* first arg is direction (-1 get info, 0 save, 1 restore), second is value*/ - double *xdir, *xval, *hoc_pgetarg(); + double *xdir, *xval; + #ifndef NRN_VERSION_GTEQ_8_2_0 + double *hoc_pgetarg(); long nrn_get_random_sequence(void* r); void nrn_set_random_sequence(void* r, int val); + #endif xdir = hoc_pgetarg(1); xval = hoc_pgetarg(2); if (_p_rng) { @@ -498,13 +529,13 @@ VERBATIM xval[0] = (double) seq; xval[1] = (double) which; } else { - xval[0] = (double)nrn_get_random_sequence(_p_rng); + xval[0] = (double)nrn_get_random_sequence(RANDCAST _p_rng); } } else{ if( usingR123 ) { nrnran123_setseq( (nrnran123_State*)_p_rng, (uint32_t)xval[0], (char)xval[1] ); } else { - nrn_set_random_sequence(_p_rng, (long)(xval[0])); + nrn_set_random_sequence(RANDCAST _p_rng, (long)(xval[0])); } } } diff --git a/model_management/mm_run_minimal/emodels_dir/singlecell-optimization/setup/__init__.py b/model_management/mm_run_minimal/emodels_dir/singlecell-optimization/setup/__init__.py index aed8ffa0..a27e8646 100644 --- a/model_management/mm_run_minimal/emodels_dir/singlecell-optimization/setup/__init__.py +++ b/model_management/mm_run_minimal/emodels_dir/singlecell-optimization/setup/__init__.py @@ -17,6 +17,7 @@ To view a copy of this license, visit https://creativecommons.org/licenses/by/4.0/legalcode or send a letter to Creative Commons, 171 Second Street, Suite 300, San Francisco, California, 94105, USA. """ + from . import evaluator from . import template diff --git a/model_management/mm_run_minimal/emodels_dir/singlecell-optimization/setup/evaluator.py b/model_management/mm_run_minimal/emodels_dir/singlecell-optimization/setup/evaluator.py index ce88d15a..335a2918 100644 --- a/model_management/mm_run_minimal/emodels_dir/singlecell-optimization/setup/evaluator.py +++ b/model_management/mm_run_minimal/emodels_dir/singlecell-optimization/setup/evaluator.py @@ -405,14 +405,16 @@ def define_protocols( "type" in protocol_definition and protocol_definition["type"] == "RatSSCxThresholdDetectionProtocol" ): - protocols_dict[ - "ThresholdDetection" - ] = protocols.RatSSCxThresholdDetectionProtocol( - "IDRest", - step_protocol_template=read_step_protocol( - "Threshold", protocol_definition["step_template"], recordings - ), - prefix=prefix, + protocols_dict["ThresholdDetection"] = ( + protocols.RatSSCxThresholdDetectionProtocol( + "IDRest", + step_protocol_template=read_step_protocol( + "Threshold", + protocol_definition["step_template"], + recordings, + ), + prefix=prefix, + ) ) else: stimuli = [] @@ -488,7 +490,6 @@ def calculate_score(self, responses, trace_check=False): class eFELFeatureExtra(eFELFeature): - """eFEL feature extra""" SERIALIZED_FIELDS = ( @@ -647,7 +648,6 @@ def calculate_score(self, responses, trace_check=False): class SingletonWeightObjective(EFeatureObjective): - """Single EPhys feature""" def __init__(self, name, feature, weight): @@ -775,7 +775,6 @@ def define_fitness_calculator( class MultiEvaluator(bpopt.evaluators.Evaluator): - """Multiple cell evaluator""" def __init__( diff --git a/model_management/mm_run_minimal/emodels_dir/singlecell-optimization/setup/protocols.py b/model_management/mm_run_minimal/emodels_dir/singlecell-optimization/setup/protocols.py index 81d19d3a..02e39092 100644 --- a/model_management/mm_run_minimal/emodels_dir/singlecell-optimization/setup/protocols.py +++ b/model_management/mm_run_minimal/emodels_dir/singlecell-optimization/setup/protocols.py @@ -38,7 +38,6 @@ class RatSSCxMainProtocol(ephys.protocols.Protocol): - """Main protocol to fit RatSSCx neuron ephys parameters. Pseudo code: @@ -193,7 +192,6 @@ def run(self, cell_model, param_values, sim=None, isolate=None): class RatSSCxRinHoldcurrentProtocol(ephys.protocols.Protocol): - """IDRest protocol to fit RatSSCx neuron ephys parameters""" def __init__( @@ -385,7 +383,6 @@ def voltage_base(self, current, cell_model, param_values, sim=None, short=False) class RatSSCxThresholdDetectionProtocol(ephys.protocols.Protocol): - """IDRest protocol to fit RatSSCx neuron ephys parameters""" def __init__( @@ -641,7 +638,6 @@ def search_spike_threshold( class StepProtocol(ephys.protocols.SweepProtocol): - """Protocol consisting of step and holding current""" def __init__( @@ -665,9 +661,11 @@ def __init__( super(StepProtocol, self).__init__( name, - stimuli=step_stimuli + [holding_stimulus] - if holding_stimulus is not None - else step_stimuli, + stimuli=( + step_stimuli + [holding_stimulus] + if holding_stimulus is not None + else step_stimuli + ), recordings=recordings, cvode_active=cvode_active, ) @@ -752,7 +750,6 @@ def step_amplitude(self): class StepThresholdProtocol(StepProtocol): - """Step protocol based on threshold""" def __init__( diff --git a/model_management/mm_run_minimal/emodels_dir/singlecell-optimization/setup/template.py b/model_management/mm_run_minimal/emodels_dir/singlecell-optimization/setup/template.py index c1370c4e..8c7044c7 100644 --- a/model_management/mm_run_minimal/emodels_dir/singlecell-optimization/setup/template.py +++ b/model_management/mm_run_minimal/emodels_dir/singlecell-optimization/setup/template.py @@ -122,12 +122,12 @@ def define_parameters(params_filename, stage=None, past_params=None): dist_param_names = definition["parameters"] else: dist_param_names = None - distributions[ - distribution - ] = ephys.parameterscalers.NrnSegmentSomaDistanceScaler( - name=distribution, - distribution=definition["fun"], - dist_param_names=dist_param_names, + distributions[distribution] = ( + ephys.parameterscalers.NrnSegmentSomaDistanceScaler( + name=distribution, + distribution=definition["fun"], + dist_param_names=dist_param_names, + ) ) params_definitions = definitions["parameters"] diff --git a/model_management/mm_run_minimal/tests/test_model_management.py b/model_management/mm_run_minimal/tests/test_model_management.py index 406cebfc..518a5418 100644 --- a/model_management/mm_run_minimal/tests/test_model_management.py +++ b/model_management/mm_run_minimal/tests/test_model_management.py @@ -18,7 +18,6 @@ Second Street, Suite 300, San Francisco, California, 94105, USA. """ - """Test to validate the model management output.""" from pathlib import Path diff --git a/optimization/opt_module/mechanisms/StochKv3.mod b/optimization/opt_module/mechanisms/StochKv3.mod index 92d955e9..06203e15 100644 --- a/optimization/opt_module/mechanisms/StochKv3.mod +++ b/optimization/opt_module/mechanisms/StochKv3.mod @@ -1,4 +1,4 @@ -TITLE StochKv2.mod +TITLE StochKv3.mod COMMENT ---------------------------------------------------------------- @@ -28,11 +28,11 @@ NEURON { USEION k READ ek WRITE ik RANGE N, eta, gk, gamma, deterministic, gkbar, ik RANGE N0, N1, n0_n1, n1_n0 - GLOBAL ninf, linf, ltau, ntau, an, bn, al, bl - GLOBAL P_an, P_bn, P_al, P_bl + RANGE ninf, linf, ltau, ntau, an, bn, al, bl + RANGE P_an, P_bn, P_al, P_bl GLOBAL vmin, vmax - :BBCOREPOINTER rng - POINTER rng + BBCOREPOINTER rng + :POINTER rng } UNITS { @@ -99,6 +99,7 @@ for comparison with Pr to decide whether to activate the synapse or not ENDCOMMENT VERBATIM +#ifndef NRN_VERSION_GTEQ_8_2_0 #include "nrnran123.h" extern int cvode_active_; @@ -106,8 +107,14 @@ extern int cvode_active_; #include #include +#ifndef CORENEURON_BUILD double nrn_random_pick(void* r); void* nrn_random_arg(int argpos); +#endif +#define RANDCAST +#else +#define RANDCAST (Rand*) +#endif ENDVERBATIM : ---------------------------------------------------------------- @@ -117,6 +124,10 @@ INITIAL { if (cvode_active_ && !deterministic) { hoc_execerror("StochKv2 with deterministic=0", "cannot be used with cvode"); } + + if( usingR123 ) { + nrnran123_setseq((nrnran123_State*)_p_rng, 0, 0); + } ENDVERBATIM eta = (gkbar / gamma) * (10000) @@ -248,9 +259,12 @@ PROCEDURE trates(v (mV)) { FUNCTION strap(x) { if (x < 0) { strap = 0 -VERBATIM - fprintf (stderr,"skv.mod:strap: negative state"); -ENDVERBATIM +: This function gets executed in the DERIVATIVE and BREAKPOINT blocks +: which should get vectorized and fprintf doesn't allow this +: Also Intel classic compilers generate invalid code when #pragma omp simd is used +:VERBATIM +: fprintf (stderr,"skv.mod:strap: negative state"); +:ENDVERBATIM } else { strap = x } @@ -273,11 +287,12 @@ PROCEDURE setRNG() { VERBATIM // For compatibility, allow for either MCellRan4 or Random123. Distinguish by the arg types // Object => MCellRan4, seeds (double) => Random123 -#if !NRNBBCORE +#ifndef CORENEURON_BUILD usingR123 = 0; if( ifarg(1) && hoc_is_double_arg(1) ) { nrnran123_State** pv = (nrnran123_State**)(&_p_rng); uint32_t a2 = 0; + uint32_t a3 = 0; if (*pv) { nrnran123_deletestream(*pv); @@ -286,7 +301,10 @@ VERBATIM if (ifarg(2)) { a2 = (uint32_t)*getarg(2); } - *pv = nrnran123_newstream((uint32_t)*getarg(1), a2); + if (ifarg(3)) { + a3 = (uint32_t)*getarg(3); + } + *pv = nrnran123_newstream3((uint32_t)*getarg(1), a2, a3); usingR123 = 1; } else if( ifarg(1) ) { void** pv = (void**)(&_p_rng); @@ -302,46 +320,56 @@ ENDVERBATIM FUNCTION urand() { VERBATIM - double value; + double value = 0.0; if( usingR123 ) { value = nrnran123_dblpick((nrnran123_State*)_p_rng); } else if (_p_rng) { - value = nrn_random_pick(_p_rng); +#ifndef CORENEURON_BUILD + value = nrn_random_pick(RANDCAST _p_rng); +#endif } else { - value = 0.5; + // see BBPBGLIB-972 + value = 0.0; } _lurand = value; ENDVERBATIM } VERBATIM -/* static void bbcore_write(double* x, int* d, int* xx, int* offset, _threadargsproto_) { if (d) { uint32_t* di = ((uint32_t*)d) + *offset; // temporary just enough to see how much space is being used if (!_p_rng) { - di[0] = 0; di[1] = 0; + di[0] = 0; di[1] = 0, di[2] = 0; }else{ nrnran123_State** pv = (nrnran123_State**)(&_p_rng); - nrnran123_getids(*pv, di, di+1); + nrnran123_getids3(*pv, di, di+1, di+2); + // write stream sequence + char which; + nrnran123_getseq(*pv, di+3, &which); + di[4] = (int)which; } -//printf("StochKv2.mod %p: bbcore_write offset=%d %d %d\n", _p, *offset, d?di[0]:-1, d?di[1]:-1); + //printf("StochKv3.mod %p: bbcore_write offset=%d %d %d\n", _p, *offset, d?di[0]:-1, d?di[1]:-1); } - *offset += 2; + *offset += 5; } static void bbcore_read(double* x, int* d, int* xx, int* offset, _threadargsproto_) { - assert(!_p_rng); uint32_t* di = ((uint32_t*)d) + *offset; - if (di[0] != 0 || di[1] != 0) + if (di[0] != 0 || di[1] != 0|| di[2] != 0) { nrnran123_State** pv = (nrnran123_State**)(&_p_rng); - *pv = nrnran123_newstream(di[0], di[1]); +#if !NRNBBCORE + if(*pv) { + nrnran123_deletestream(*pv); + } +#endif + *pv = nrnran123_newstream3(di[0], di[1], di[2]); + nrnran123_setseq(*pv, di[3], (char)di[4]); } -//printf("StochKv2.mod %p: bbcore_read offset=%d %d %d\n", _p, *offset, di[0], di[1]); - *offset += 2; + //printf("StochKv3.mod %p: bbcore_read offset=%d %d %d\n", _p, *offset, di[0], di[1]); + *offset += 5; } -*/ ENDVERBATIM : Returns random numbers drawn from a binomial distribution @@ -471,13 +499,16 @@ VERBATIM FUNCTION bbsavestate() { bbsavestate = 0 VERBATIM -#ifdef ENABLE_SAVE_STATE + #ifndef CORENEURON_BUILD // TODO: since N0,N1 are no longer state variables, they will need to be written using this callback // provided that it is the version that supports multivalue writing /* first arg is direction (-1 get info, 0 save, 1 restore), second is value*/ - double *xdir, *xval, *hoc_pgetarg(); + double *xdir, *xval; + #ifndef NRN_VERSION_GTEQ_8_2_0 + double *hoc_pgetarg(); long nrn_get_random_sequence(void* r); void nrn_set_random_sequence(void* r, int val); + #endif xdir = hoc_pgetarg(1); xval = hoc_pgetarg(2); if (_p_rng) { @@ -498,13 +529,13 @@ VERBATIM xval[0] = (double) seq; xval[1] = (double) which; } else { - xval[0] = (double)nrn_get_random_sequence(_p_rng); + xval[0] = (double)nrn_get_random_sequence(RANDCAST _p_rng); } } else{ if( usingR123 ) { nrnran123_setseq( (nrnran123_State*)_p_rng, (uint32_t)xval[0], (char)xval[1] ); } else { - nrn_set_random_sequence(_p_rng, (long)(xval[0])); + nrn_set_random_sequence(RANDCAST _p_rng, (long)(xval[0])); } } } diff --git a/optimization/opt_module/setup/__init__.py b/optimization/opt_module/setup/__init__.py index 0c9e82a9..3dfb55d6 100644 --- a/optimization/opt_module/setup/__init__.py +++ b/optimization/opt_module/setup/__init__.py @@ -17,5 +17,6 @@ To view a copy of this license, visit https://creativecommons.org/licenses/by/4.0/legalcode or send a letter to Creative Commons, 171 Second Street, Suite 300, San Francisco, California, 94105, USA. """ + from . import evaluator from . import template diff --git a/optimization/opt_module/setup/evaluator.py b/optimization/opt_module/setup/evaluator.py index bef7449c..53bc6155 100644 --- a/optimization/opt_module/setup/evaluator.py +++ b/optimization/opt_module/setup/evaluator.py @@ -408,14 +408,16 @@ def define_protocols( "type" in protocol_definition and protocol_definition["type"] == "RatSSCxThresholdDetectionProtocol" ): - protocols_dict[ - "ThresholdDetection" - ] = protocols.RatSSCxThresholdDetectionProtocol( - "IDRest", - step_protocol_template=read_step_protocol( - "Threshold", protocol_definition["step_template"], recordings - ), - prefix=prefix, + protocols_dict["ThresholdDetection"] = ( + protocols.RatSSCxThresholdDetectionProtocol( + "IDRest", + step_protocol_template=read_step_protocol( + "Threshold", + protocol_definition["step_template"], + recordings, + ), + prefix=prefix, + ) ) else: stimuli = [] @@ -491,7 +493,6 @@ def calculate_score(self, responses, trace_check=False): class eFELFeatureExtra(eFELFeature): - """eFEL feature extra""" SERIALIZED_FIELDS = ( @@ -666,7 +667,6 @@ def calculate_score(self, responses, trace_check=False): class SingletonWeightObjective(EFeatureObjective): - """Single EPhys feature""" def __init__(self, name, feature, weight): @@ -794,7 +794,6 @@ def define_fitness_calculator( class MultiEvaluator(bpopt.evaluators.Evaluator): - """Multiple cell evaluator""" def __init__( diff --git a/optimization/opt_module/setup/protocols.py b/optimization/opt_module/setup/protocols.py index 2111553a..540dc5e2 100644 --- a/optimization/opt_module/setup/protocols.py +++ b/optimization/opt_module/setup/protocols.py @@ -18,7 +18,6 @@ Second Street, Suite 300, San Francisco, California, 94105, USA. """ - """Protocols involved in fitting of RatSSCx neuron ephys parameters""" import numpy @@ -45,7 +44,6 @@ class RatSSCxMainProtocol(ephys.protocols.Protocol): - """Main protocol to fit RatSSCx neuron ephys parameters. Pseudo code: @@ -220,7 +218,6 @@ def run(self, cell_model, param_values, sim=None, isolate=None): class RatSSCxRinHoldcurrentProtocol(ephys.protocols.Protocol): - """IDRest protocol to fit RatSSCx neuron ephys parameters""" def __init__( @@ -424,7 +421,6 @@ def voltage_base(self, current, cell_model, param_values, sim=None, short=False) class RatSSCxThresholdDetectionProtocol(ephys.protocols.Protocol): - """IDRest protocol to fit RatSSCx neuron ephys parameters""" def __init__( @@ -704,7 +700,6 @@ def search_spike_threshold( class StepProtocol(ephys.protocols.SweepProtocol): - """Protocol consisting of step and holding current""" def __init__( @@ -728,9 +723,11 @@ def __init__( super(StepProtocol, self).__init__( name, - stimuli=step_stimuli + [holding_stimulus] - if holding_stimulus is not None - else step_stimuli, + stimuli=( + step_stimuli + [holding_stimulus] + if holding_stimulus is not None + else step_stimuli + ), recordings=recordings, cvode_active=cvode_active, ) @@ -815,7 +812,6 @@ def step_amplitude(self): class StepThresholdProtocol(StepProtocol): - """Step protocol based on threshold""" def __init__( diff --git a/optimization/opt_module/setup/template.py b/optimization/opt_module/setup/template.py index 36da1505..d23a86d0 100644 --- a/optimization/opt_module/setup/template.py +++ b/optimization/opt_module/setup/template.py @@ -130,12 +130,12 @@ def define_parameters(params_filename, stage=None, past_params=None): dist_param_names = definition["parameters"] else: dist_param_names = None - distributions[ - distribution - ] = ephys.parameterscalers.NrnSegmentSomaDistanceScaler( - name=distribution, - distribution=definition["fun"], - dist_param_names=dist_param_names, + distributions[distribution] = ( + ephys.parameterscalers.NrnSegmentSomaDistanceScaler( + name=distribution, + distribution=definition["fun"], + dist_param_names=dist_param_names, + ) ) params_definitions = definitions["parameters"] diff --git a/optimization/opt_module/tools/plot.py b/optimization/opt_module/tools/plot.py index d967981e..344bde21 100644 --- a/optimization/opt_module/tools/plot.py +++ b/optimization/opt_module/tools/plot.py @@ -18,7 +18,6 @@ Second Street, Suite 300, San Francisco, California, 94105, USA. """ - import numpy import seaborn as sns import pandas as pd diff --git a/optimization/opt_module/tools/plottools.py b/optimization/opt_module/tools/plottools.py index b1d8b81e..05d2de3a 100644 --- a/optimization/opt_module/tools/plottools.py +++ b/optimization/opt_module/tools/plottools.py @@ -18,7 +18,6 @@ Second Street, Suite 300, San Francisco, California, 94105, USA. """ - import matplotlib import matplotlib.pyplot as plt diff --git a/optimization/tests/test_optimization.py b/optimization/tests/test_optimization.py index c163c179..fd08e30b 100644 --- a/optimization/tests/test_optimization.py +++ b/optimization/tests/test_optimization.py @@ -282,4 +282,4 @@ def test_responses(self): for c_key, c_value in ca_responses.items(): calcium_response = c_value.response calcium_response_gt = np.loadtxt(f"tests/calcium_responses/{c_key}.dat") - assert np.allclose(calcium_response, calcium_response_gt, atol=1e-6) + assert np.allclose(calcium_response, calcium_response_gt, atol=1e-2) diff --git a/requirements.txt b/requirements.txt index 35bec456..d06b6971 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,7 +3,7 @@ efel==3.1.96 scipy bluepyopt>=1.13,<=1.13.12 numpy<1.24 -NEURON>=7.8.0,<8.0.0 +NEURON>=7.8.0,<=8.2.4 matplotlib seaborn bluepymm>=0.9.7 diff --git a/somatic_validation/mechanisms/StochKv3.mod b/somatic_validation/mechanisms/StochKv3.mod index 92d955e9..06203e15 100644 --- a/somatic_validation/mechanisms/StochKv3.mod +++ b/somatic_validation/mechanisms/StochKv3.mod @@ -1,4 +1,4 @@ -TITLE StochKv2.mod +TITLE StochKv3.mod COMMENT ---------------------------------------------------------------- @@ -28,11 +28,11 @@ NEURON { USEION k READ ek WRITE ik RANGE N, eta, gk, gamma, deterministic, gkbar, ik RANGE N0, N1, n0_n1, n1_n0 - GLOBAL ninf, linf, ltau, ntau, an, bn, al, bl - GLOBAL P_an, P_bn, P_al, P_bl + RANGE ninf, linf, ltau, ntau, an, bn, al, bl + RANGE P_an, P_bn, P_al, P_bl GLOBAL vmin, vmax - :BBCOREPOINTER rng - POINTER rng + BBCOREPOINTER rng + :POINTER rng } UNITS { @@ -99,6 +99,7 @@ for comparison with Pr to decide whether to activate the synapse or not ENDCOMMENT VERBATIM +#ifndef NRN_VERSION_GTEQ_8_2_0 #include "nrnran123.h" extern int cvode_active_; @@ -106,8 +107,14 @@ extern int cvode_active_; #include #include +#ifndef CORENEURON_BUILD double nrn_random_pick(void* r); void* nrn_random_arg(int argpos); +#endif +#define RANDCAST +#else +#define RANDCAST (Rand*) +#endif ENDVERBATIM : ---------------------------------------------------------------- @@ -117,6 +124,10 @@ INITIAL { if (cvode_active_ && !deterministic) { hoc_execerror("StochKv2 with deterministic=0", "cannot be used with cvode"); } + + if( usingR123 ) { + nrnran123_setseq((nrnran123_State*)_p_rng, 0, 0); + } ENDVERBATIM eta = (gkbar / gamma) * (10000) @@ -248,9 +259,12 @@ PROCEDURE trates(v (mV)) { FUNCTION strap(x) { if (x < 0) { strap = 0 -VERBATIM - fprintf (stderr,"skv.mod:strap: negative state"); -ENDVERBATIM +: This function gets executed in the DERIVATIVE and BREAKPOINT blocks +: which should get vectorized and fprintf doesn't allow this +: Also Intel classic compilers generate invalid code when #pragma omp simd is used +:VERBATIM +: fprintf (stderr,"skv.mod:strap: negative state"); +:ENDVERBATIM } else { strap = x } @@ -273,11 +287,12 @@ PROCEDURE setRNG() { VERBATIM // For compatibility, allow for either MCellRan4 or Random123. Distinguish by the arg types // Object => MCellRan4, seeds (double) => Random123 -#if !NRNBBCORE +#ifndef CORENEURON_BUILD usingR123 = 0; if( ifarg(1) && hoc_is_double_arg(1) ) { nrnran123_State** pv = (nrnran123_State**)(&_p_rng); uint32_t a2 = 0; + uint32_t a3 = 0; if (*pv) { nrnran123_deletestream(*pv); @@ -286,7 +301,10 @@ VERBATIM if (ifarg(2)) { a2 = (uint32_t)*getarg(2); } - *pv = nrnran123_newstream((uint32_t)*getarg(1), a2); + if (ifarg(3)) { + a3 = (uint32_t)*getarg(3); + } + *pv = nrnran123_newstream3((uint32_t)*getarg(1), a2, a3); usingR123 = 1; } else if( ifarg(1) ) { void** pv = (void**)(&_p_rng); @@ -302,46 +320,56 @@ ENDVERBATIM FUNCTION urand() { VERBATIM - double value; + double value = 0.0; if( usingR123 ) { value = nrnran123_dblpick((nrnran123_State*)_p_rng); } else if (_p_rng) { - value = nrn_random_pick(_p_rng); +#ifndef CORENEURON_BUILD + value = nrn_random_pick(RANDCAST _p_rng); +#endif } else { - value = 0.5; + // see BBPBGLIB-972 + value = 0.0; } _lurand = value; ENDVERBATIM } VERBATIM -/* static void bbcore_write(double* x, int* d, int* xx, int* offset, _threadargsproto_) { if (d) { uint32_t* di = ((uint32_t*)d) + *offset; // temporary just enough to see how much space is being used if (!_p_rng) { - di[0] = 0; di[1] = 0; + di[0] = 0; di[1] = 0, di[2] = 0; }else{ nrnran123_State** pv = (nrnran123_State**)(&_p_rng); - nrnran123_getids(*pv, di, di+1); + nrnran123_getids3(*pv, di, di+1, di+2); + // write stream sequence + char which; + nrnran123_getseq(*pv, di+3, &which); + di[4] = (int)which; } -//printf("StochKv2.mod %p: bbcore_write offset=%d %d %d\n", _p, *offset, d?di[0]:-1, d?di[1]:-1); + //printf("StochKv3.mod %p: bbcore_write offset=%d %d %d\n", _p, *offset, d?di[0]:-1, d?di[1]:-1); } - *offset += 2; + *offset += 5; } static void bbcore_read(double* x, int* d, int* xx, int* offset, _threadargsproto_) { - assert(!_p_rng); uint32_t* di = ((uint32_t*)d) + *offset; - if (di[0] != 0 || di[1] != 0) + if (di[0] != 0 || di[1] != 0|| di[2] != 0) { nrnran123_State** pv = (nrnran123_State**)(&_p_rng); - *pv = nrnran123_newstream(di[0], di[1]); +#if !NRNBBCORE + if(*pv) { + nrnran123_deletestream(*pv); + } +#endif + *pv = nrnran123_newstream3(di[0], di[1], di[2]); + nrnran123_setseq(*pv, di[3], (char)di[4]); } -//printf("StochKv2.mod %p: bbcore_read offset=%d %d %d\n", _p, *offset, di[0], di[1]); - *offset += 2; + //printf("StochKv3.mod %p: bbcore_read offset=%d %d %d\n", _p, *offset, di[0], di[1]); + *offset += 5; } -*/ ENDVERBATIM : Returns random numbers drawn from a binomial distribution @@ -471,13 +499,16 @@ VERBATIM FUNCTION bbsavestate() { bbsavestate = 0 VERBATIM -#ifdef ENABLE_SAVE_STATE + #ifndef CORENEURON_BUILD // TODO: since N0,N1 are no longer state variables, they will need to be written using this callback // provided that it is the version that supports multivalue writing /* first arg is direction (-1 get info, 0 save, 1 restore), second is value*/ - double *xdir, *xval, *hoc_pgetarg(); + double *xdir, *xval; + #ifndef NRN_VERSION_GTEQ_8_2_0 + double *hoc_pgetarg(); long nrn_get_random_sequence(void* r); void nrn_set_random_sequence(void* r, int val); + #endif xdir = hoc_pgetarg(1); xval = hoc_pgetarg(2); if (_p_rng) { @@ -498,13 +529,13 @@ VERBATIM xval[0] = (double) seq; xval[1] = (double) which; } else { - xval[0] = (double)nrn_get_random_sequence(_p_rng); + xval[0] = (double)nrn_get_random_sequence(RANDCAST _p_rng); } } else{ if( usingR123 ) { nrnran123_setseq( (nrnran123_State*)_p_rng, (uint32_t)xval[0], (char)xval[1] ); } else { - nrn_set_random_sequence(_p_rng, (long)(xval[0])); + nrn_set_random_sequence(RANDCAST _p_rng, (long)(xval[0])); } } } diff --git a/somatic_validation/model.py b/somatic_validation/model.py index 5d341d62..e462953d 100644 --- a/somatic_validation/model.py +++ b/somatic_validation/model.py @@ -18,7 +18,6 @@ Second Street, Suite 300, San Francisco, California, 94105, USA. """ - import os import collections import json @@ -110,12 +109,12 @@ def define_parameters(path_params): dist_param_names = definition["parameters"] else: dist_param_names = None - distributions[ - distribution - ] = ephys.parameterscalers.NrnSegmentSomaDistanceScaler( - name=distribution, - distribution=definition["fun"], - dist_param_names=dist_param_names, + distributions[distribution] = ( + ephys.parameterscalers.NrnSegmentSomaDistanceScaler( + name=distribution, + distribution=definition["fun"], + dist_param_names=dist_param_names, + ) ) params_definitions = definitions["parameters"] diff --git a/somatic_validation/stimuli.py b/somatic_validation/stimuli.py index f3854e93..28b00788 100644 --- a/somatic_validation/stimuli.py +++ b/somatic_validation/stimuli.py @@ -28,14 +28,12 @@ class Stimulus(object): - """Stimulus protocol""" pass class NrnHDPulse(Stimulus): - """Square pulse current clamp injection""" def __init__( diff --git a/validation/lib/attenuation.py b/validation/lib/attenuation.py index 6cb86d64..d75b5285 100644 --- a/validation/lib/attenuation.py +++ b/validation/lib/attenuation.py @@ -826,7 +826,6 @@ def EPSP(evaluator, emodels_path, points): class NrnSynStim(ephys.stimuli.Stimulus): - """Square pulse current clamp injection""" def __init__( diff --git a/validation/lib/plottools.py b/validation/lib/plottools.py index b338da0a..0ebedb48 100644 --- a/validation/lib/plottools.py +++ b/validation/lib/plottools.py @@ -18,7 +18,6 @@ Second Street, Suite 300, San Francisco, California, 94105, USA. """ - import matplotlib import matplotlib.pyplot as plt diff --git a/validation/lib/setup/__init__.py b/validation/lib/setup/__init__.py index 749071d3..5fa4349b 100644 --- a/validation/lib/setup/__init__.py +++ b/validation/lib/setup/__init__.py @@ -17,4 +17,5 @@ To view a copy of this license, visit https://creativecommons.org/licenses/by/4.0/legalcode or send a letter to Creative Commons, 171 Second Street, Suite 300, San Francisco, California, 94105, USA. """ + from . import evaluator diff --git a/validation/lib/setup/evaluator.py b/validation/lib/setup/evaluator.py index dba8e4f1..b4733d27 100644 --- a/validation/lib/setup/evaluator.py +++ b/validation/lib/setup/evaluator.py @@ -260,14 +260,16 @@ def define_protocols( "type" in protocol_definition and protocol_definition["type"] == "RatSSCxThresholdDetectionProtocol" ): - protocols_dict[ - "ThresholdDetection" - ] = protocols.RatSSCxThresholdDetectionProtocol( - "IDRest", - step_protocol_template=read_step_protocol( - "Threshold", protocol_definition["step_template"], recordings - ), - prefix=prefix, + protocols_dict["ThresholdDetection"] = ( + protocols.RatSSCxThresholdDetectionProtocol( + "IDRest", + step_protocol_template=read_step_protocol( + "Threshold", + protocol_definition["step_template"], + recordings, + ), + prefix=prefix, + ) ) else: stimuli = [] @@ -331,7 +333,6 @@ def define_protocols( class eFELFeatureExtra(eFELFeature): - """eFEL feature extra""" SERIALIZED_FIELDS = ( @@ -419,7 +420,6 @@ def __init__( class SingletonWeightObjective(EFeatureObjective): - """Single EPhys feature""" def __init__(self, name, feature, weight): @@ -547,7 +547,6 @@ def define_fitness_calculator( class MultiEvaluator(bpopt.evaluators.Evaluator): - """Multiple cell evaluator""" def __init__( diff --git a/validation/lib/setup/protocols.py b/validation/lib/setup/protocols.py index bc8ac74c..d0c8f6b8 100644 --- a/validation/lib/setup/protocols.py +++ b/validation/lib/setup/protocols.py @@ -18,7 +18,6 @@ Second Street, Suite 300, San Francisco, California, 94105, USA. """ - """Protocols involved in fitting of RatSSCx neuron ephys parameters""" import numpy @@ -34,7 +33,6 @@ class RatSSCxMainProtocol(ephys.protocols.Protocol): - """Main protocol to fit RatSSCx neuron ephys parameters. Pseudo code: @@ -115,7 +113,6 @@ def rin_efeature(self, value): class RatSSCxRinHoldcurrentProtocol(ephys.protocols.Protocol): - """IDRest protocol to fit RatSSCx neuron ephys parameters""" def __init__( @@ -156,7 +153,6 @@ def subprotocols(self): class RatSSCxThresholdDetectionProtocol(ephys.protocols.Protocol): - """IDRest protocol to fit RatSSCx neuron ephys parameters""" def __init__( @@ -387,7 +383,6 @@ def search_spike_threshold( class StepProtocol(ephys.protocols.SweepProtocol): - """Protocol consisting of step and holding current""" def __init__( @@ -411,9 +406,11 @@ def __init__( super(StepProtocol, self).__init__( name, - stimuli=step_stimuli + [holding_stimulus] - if holding_stimulus is not None - else step_stimuli, + stimuli=( + step_stimuli + [holding_stimulus] + if holding_stimulus is not None + else step_stimuli + ), recordings=recordings, cvode_active=cvode_active, ) @@ -479,7 +476,6 @@ def step_amplitude(self): class StepThresholdProtocol(StepProtocol): - """Step protocol based on threshold""" def __init__( diff --git a/validation/lib/setup/template.py b/validation/lib/setup/template.py index 41ed85b4..027f5fab 100644 --- a/validation/lib/setup/template.py +++ b/validation/lib/setup/template.py @@ -122,12 +122,12 @@ def define_parameters(params_filename, stage=None, past_params=None): dist_param_names = definition["parameters"] else: dist_param_names = None - distributions[ - distribution - ] = ephys.parameterscalers.NrnSegmentSomaDistanceScaler( - name=distribution, - distribution=definition["fun"], - dist_param_names=dist_param_names, + distributions[distribution] = ( + ephys.parameterscalers.NrnSegmentSomaDistanceScaler( + name=distribution, + distribution=definition["fun"], + dist_param_names=dist_param_names, + ) ) params_definitions = definitions["parameters"] diff --git a/validation/lib/tools.py b/validation/lib/tools.py index d54c9282..6f5077db 100644 --- a/validation/lib/tools.py +++ b/validation/lib/tools.py @@ -18,7 +18,6 @@ Second Street, Suite 300, San Francisco, California, 94105, USA. """ - import contextlib import os diff --git a/validation/mechanisms/StochKv3.mod b/validation/mechanisms/StochKv3.mod index 92d955e9..06203e15 100644 --- a/validation/mechanisms/StochKv3.mod +++ b/validation/mechanisms/StochKv3.mod @@ -1,4 +1,4 @@ -TITLE StochKv2.mod +TITLE StochKv3.mod COMMENT ---------------------------------------------------------------- @@ -28,11 +28,11 @@ NEURON { USEION k READ ek WRITE ik RANGE N, eta, gk, gamma, deterministic, gkbar, ik RANGE N0, N1, n0_n1, n1_n0 - GLOBAL ninf, linf, ltau, ntau, an, bn, al, bl - GLOBAL P_an, P_bn, P_al, P_bl + RANGE ninf, linf, ltau, ntau, an, bn, al, bl + RANGE P_an, P_bn, P_al, P_bl GLOBAL vmin, vmax - :BBCOREPOINTER rng - POINTER rng + BBCOREPOINTER rng + :POINTER rng } UNITS { @@ -99,6 +99,7 @@ for comparison with Pr to decide whether to activate the synapse or not ENDCOMMENT VERBATIM +#ifndef NRN_VERSION_GTEQ_8_2_0 #include "nrnran123.h" extern int cvode_active_; @@ -106,8 +107,14 @@ extern int cvode_active_; #include #include +#ifndef CORENEURON_BUILD double nrn_random_pick(void* r); void* nrn_random_arg(int argpos); +#endif +#define RANDCAST +#else +#define RANDCAST (Rand*) +#endif ENDVERBATIM : ---------------------------------------------------------------- @@ -117,6 +124,10 @@ INITIAL { if (cvode_active_ && !deterministic) { hoc_execerror("StochKv2 with deterministic=0", "cannot be used with cvode"); } + + if( usingR123 ) { + nrnran123_setseq((nrnran123_State*)_p_rng, 0, 0); + } ENDVERBATIM eta = (gkbar / gamma) * (10000) @@ -248,9 +259,12 @@ PROCEDURE trates(v (mV)) { FUNCTION strap(x) { if (x < 0) { strap = 0 -VERBATIM - fprintf (stderr,"skv.mod:strap: negative state"); -ENDVERBATIM +: This function gets executed in the DERIVATIVE and BREAKPOINT blocks +: which should get vectorized and fprintf doesn't allow this +: Also Intel classic compilers generate invalid code when #pragma omp simd is used +:VERBATIM +: fprintf (stderr,"skv.mod:strap: negative state"); +:ENDVERBATIM } else { strap = x } @@ -273,11 +287,12 @@ PROCEDURE setRNG() { VERBATIM // For compatibility, allow for either MCellRan4 or Random123. Distinguish by the arg types // Object => MCellRan4, seeds (double) => Random123 -#if !NRNBBCORE +#ifndef CORENEURON_BUILD usingR123 = 0; if( ifarg(1) && hoc_is_double_arg(1) ) { nrnran123_State** pv = (nrnran123_State**)(&_p_rng); uint32_t a2 = 0; + uint32_t a3 = 0; if (*pv) { nrnran123_deletestream(*pv); @@ -286,7 +301,10 @@ VERBATIM if (ifarg(2)) { a2 = (uint32_t)*getarg(2); } - *pv = nrnran123_newstream((uint32_t)*getarg(1), a2); + if (ifarg(3)) { + a3 = (uint32_t)*getarg(3); + } + *pv = nrnran123_newstream3((uint32_t)*getarg(1), a2, a3); usingR123 = 1; } else if( ifarg(1) ) { void** pv = (void**)(&_p_rng); @@ -302,46 +320,56 @@ ENDVERBATIM FUNCTION urand() { VERBATIM - double value; + double value = 0.0; if( usingR123 ) { value = nrnran123_dblpick((nrnran123_State*)_p_rng); } else if (_p_rng) { - value = nrn_random_pick(_p_rng); +#ifndef CORENEURON_BUILD + value = nrn_random_pick(RANDCAST _p_rng); +#endif } else { - value = 0.5; + // see BBPBGLIB-972 + value = 0.0; } _lurand = value; ENDVERBATIM } VERBATIM -/* static void bbcore_write(double* x, int* d, int* xx, int* offset, _threadargsproto_) { if (d) { uint32_t* di = ((uint32_t*)d) + *offset; // temporary just enough to see how much space is being used if (!_p_rng) { - di[0] = 0; di[1] = 0; + di[0] = 0; di[1] = 0, di[2] = 0; }else{ nrnran123_State** pv = (nrnran123_State**)(&_p_rng); - nrnran123_getids(*pv, di, di+1); + nrnran123_getids3(*pv, di, di+1, di+2); + // write stream sequence + char which; + nrnran123_getseq(*pv, di+3, &which); + di[4] = (int)which; } -//printf("StochKv2.mod %p: bbcore_write offset=%d %d %d\n", _p, *offset, d?di[0]:-1, d?di[1]:-1); + //printf("StochKv3.mod %p: bbcore_write offset=%d %d %d\n", _p, *offset, d?di[0]:-1, d?di[1]:-1); } - *offset += 2; + *offset += 5; } static void bbcore_read(double* x, int* d, int* xx, int* offset, _threadargsproto_) { - assert(!_p_rng); uint32_t* di = ((uint32_t*)d) + *offset; - if (di[0] != 0 || di[1] != 0) + if (di[0] != 0 || di[1] != 0|| di[2] != 0) { nrnran123_State** pv = (nrnran123_State**)(&_p_rng); - *pv = nrnran123_newstream(di[0], di[1]); +#if !NRNBBCORE + if(*pv) { + nrnran123_deletestream(*pv); + } +#endif + *pv = nrnran123_newstream3(di[0], di[1], di[2]); + nrnran123_setseq(*pv, di[3], (char)di[4]); } -//printf("StochKv2.mod %p: bbcore_read offset=%d %d %d\n", _p, *offset, di[0], di[1]); - *offset += 2; + //printf("StochKv3.mod %p: bbcore_read offset=%d %d %d\n", _p, *offset, di[0], di[1]); + *offset += 5; } -*/ ENDVERBATIM : Returns random numbers drawn from a binomial distribution @@ -471,13 +499,16 @@ VERBATIM FUNCTION bbsavestate() { bbsavestate = 0 VERBATIM -#ifdef ENABLE_SAVE_STATE + #ifndef CORENEURON_BUILD // TODO: since N0,N1 are no longer state variables, they will need to be written using this callback // provided that it is the version that supports multivalue writing /* first arg is direction (-1 get info, 0 save, 1 restore), second is value*/ - double *xdir, *xval, *hoc_pgetarg(); + double *xdir, *xval; + #ifndef NRN_VERSION_GTEQ_8_2_0 + double *hoc_pgetarg(); long nrn_get_random_sequence(void* r); void nrn_set_random_sequence(void* r, int val); + #endif xdir = hoc_pgetarg(1); xval = hoc_pgetarg(2); if (_p_rng) { @@ -498,13 +529,13 @@ VERBATIM xval[0] = (double) seq; xval[1] = (double) which; } else { - xval[0] = (double)nrn_get_random_sequence(_p_rng); + xval[0] = (double)nrn_get_random_sequence(RANDCAST _p_rng); } } else{ if( usingR123 ) { nrnran123_setseq( (nrnran123_State*)_p_rng, (uint32_t)xval[0], (char)xval[1] ); } else { - nrn_set_random_sequence(_p_rng, (long)(xval[0])); + nrn_set_random_sequence(RANDCAST _p_rng, (long)(xval[0])); } } } diff --git a/validation/tests/test_validation.py b/validation/tests/test_validation.py index ec2ff558..ade1ff9e 100644 --- a/validation/tests/test_validation.py +++ b/validation/tests/test_validation.py @@ -18,7 +18,6 @@ Second Street, Suite 300, San Francisco, California, 94105, USA. """ - """Test the validation results on a single sample.""" import json