Make StochKv3.mod NEURON 9.0a compatible (#35)

* Make StochKv3.mod NEURON9 compatible * black fix * Update README * limit NEURON to 8.0.2 * relax NEURON upper version until 8.2.4 * update precision of opt test * undo threshold_current ground truth update --------- Co-authored-by: Anil Tuncel <[email protected]>
BlueBrain · Jun 18, 2024 · d30e867 · d30e867
1 parent c24f096
commit d30e867
Show file tree

Hide file tree

Showing 28 changed files with 319 additions and 215 deletions.
diff --git 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) <https://doi.org/10.1016/j.patter.2023.100855>`_. 
+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) <https://doi.org/10.1016/j.patter.2023.100855>`_.
 
-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 <https://github.com/BlueBrain/SSCxEModelExamples/tree/main/validation/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 <https://github.com/B
 
 Note that this is the only step that does not use the ``requirements.txt`` in the main directory.
 
-`somatic-val-requirements.txt <https://github.com/BlueBrain/SSCxEModelExamples/blob/main/somatic_validation/somatic-val-requirements.txt>`_ needs to be installed and the mechanisms need to be compiled with the following command before running the notebooks:: 
+`somatic-val-requirements.txt <https://github.com/BlueBrain/SSCxEModelExamples/blob/main/somatic_validation/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 <https://github.com/BlueBrain/SSCxEModelExamples/blob/main/somatic_validation/feature-extraction.ipynb>`_, the results of this extraction can be found in the `somatic_validation/L5TPC <https://github.com/BlueBrain/SSCxEModelExamples/tree/main/somatic_validation/L5TPC>`_ folder. To run and visualize results of the somatic validation run `somatic-validation.ipynb <https://github.com/BlueBrain/SSCxEModelExamples/blob/main/somatic_validation/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 <https://github.com/BlueBrain/SSCxEModelExamples/releases/tag/1.0.1>`. Users with other versions may encounter some small discrepancies in results.
 The `requirements.txt <https://github.com/BlueBrain/SSCxEModelExamples/blob/main/requirements.txt>`_ at the main directory should be used for all steps except for the somatic validations.
 Install `somatic-val-requirements.txt <https://github.com/BlueBrain/SSCxEModelExamples/blob/main/somatic_validation/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 <https://creativecommons.org/licenses/by/4.0/>`_ 
+This work is licensed under `Creative Commons (CC BY) 4.0 <https://creativecommons.org/licenses/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
@@ -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,15 +99,22 @@ 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_;
 
 #include <stdlib.h>
 #include <stdio.h>
 #include <math.h>
 
+#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
@@ -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
@@ -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
@@ -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
@@ -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"]