Segfault when running tests with PyNN

[JCGoran]

Context

https://github.com/NeuralEnsemble/PyNN/ uses NEURON. It is also part of the EBRAINS 2.0 software distribution (ESD).

❗ NOTE: in order to make NEURON 9 part of the next ESD release, this needs to be resolved by 17-10-2025.

Overview of the issue

I've added fixes for NEURON 9 compatibility in PyNN here: NeuralEnsemble/PyNN#822
Unfortunately, the test here segfaults: https://github.com/neuronsimulator/PyNN/blob/766233794eef65f79de52b368d780becd393c224/test/system/scenarios/test_connectors.py#L51

Note that running just that test works, so I'm guessing NEURON does not clean-up or set-up certain objects properly since the tests all share the same NEURON instance.

Expected result/behavior

No segfault.

NEURON setup

• Version: 9.0.0 and newer
• Installation method: any

Minimal working example - MWE

git clone --branch jelic/neuron_9_compat https://github.com/neuronsimulator/PyNN
cd PyNN
python3 -m venv env && source env/bin/activate
pip install pytest matplotlib 'neuron>=9' .
python -m pytest -vvv

Output:

test/system/scenarios/test_connectors.py::test_fixed_number_pre_no_replacement[neuron] Fatal Python error: Segmentation fault

Current thread 0x000075e3a35f9080 (most recent call first):
  File "PyNN/pyNN/neuron/simulator.py", line 266 in reset
  File "PyNN/pyNN/neuron/simulator.py", line 257 in clear
  File "PyNN/pyNN/neuron/control.py", line 44 in setup
  File "PyNN/test/system/scenarios/test_connectors.py", line 52 in test_fixed_number_pre_no_replacement

Compiling NEURON with Debug type gives the source of the problem the following mod file:

in _net_init (_pnt=0x2dbb350, _args=0x2dbb550, _lflag=0) at PyNN/pyNN/neuron/nmodl/x86_64/tsodyksmarkram.cpp:316
316      t_last = t ;

which comes from here: https://github.com/neuronsimulator/PyNN/blob/766233794eef65f79de52b368d780becd393c224/pyNN/neuron/nmodl/tsodyksmarkram.mod#L37-L40

[JCGoran]

The following test seems to fail on the same line as well:

python -m pytest -vv -k 'test_scenario4[neuron]'

with:

Current thread 0x00007835e3a9f080 (most recent call first):
  File "PyNN/pyNN/neuron/simulator.py", line 275 in _pre_run
  File "PyNN/pyNN/neuron/simulator.py", line 299 in run_until
  File "PyNN/pyNN/common/control.py", line 96 in run_until
  File "PyNN/pyNN/common/control.py", line 114 in run
  File "PyNN/test/system/scenarios/test_scenario4.py", line 59 in test_scenario4

Note that the previous test (test_fixed_number_pre_no_replacement[neuron]) can be fixed by using pytest-forked (from here), but this one always seems to fail.

[JCGoran]

I've narrowed down the source of the problem to the SoA mega-commit: #2027
Unfortunately, even after restoring the full branch, due to build problems I couldn't pinpoint if there was a change there that caused the issue.

[JCGoran]

After a discussion with @nrnhines, it was discovered that adding THREADSAFE to the mod file in question allows us to avoid the segfault. The difference in generated code is as follows:

diff --git a/tsodyksmarkram_unsafe.cpp b/tsodyksmarkram_safe.cpp
index fb12f514..c6db3fe2 100644
--- a/tsodyksmarkram_unsafe.cpp
+++ b/tsodyksmarkram_safe.cpp
@@ -1,6 +1,6 @@
 /* Created by Language version: 7.7.0 */
-/* NOT VECTORIZED */
-#define NRN_VECTORIZED 0
+/* VECTORIZED */
+#define NRN_VECTORIZED 1
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
@@ -16,11 +16,8 @@
 #include "nrnconf.h"
 // clang-format on
 #include "neuron/cache/mechanism_range.hpp"
-#include <vector>
-using std::size_t;
-static auto& std_cerr_stream = std::cerr;
 static constexpr auto number_of_datum_variables = 3;
-static constexpr auto number_of_floating_point_variables = 11;
+static constexpr auto number_of_floating_point_variables = 12;
 namespace {
 template <typename T>
 using _nrn_mechanism_std_vector = std::vector<T>;
@@ -53,20 +50,20 @@ void _nrn_mechanism_register_data_fields(Args&&... args) {
 #define nrn_state _nrn_state__TsodyksMarkramWA
 #define _net_receive _net_receive__TsodyksMarkramWA 
  
-#define _threadargscomma_ /**/
-#define _threadargsprotocomma_ /**/
-#define _internalthreadargsprotocomma_ /**/
-#define _threadargs_ /**/
-#define _threadargsproto_ /**/
-#define _internalthreadargsproto_ /**/
+#define _threadargscomma_ _ml, _iml, _ppvar, _thread, _globals, _nt,
+#define _threadargsprotocomma_ Memb_list* _ml, size_t _iml, Datum* _ppvar, Datum* _thread, double* _globals, NrnThread* _nt,
+#define _internalthreadargsprotocomma_ _nrn_mechanism_cache_range* _ml, size_t _iml, Datum* _ppvar, Datum* _thread, double* _globals, NrnThread* _nt,
+#define _threadargs_ _ml, _iml, _ppvar, _thread, _globals, _nt
+#define _threadargsproto_ Memb_list* _ml, size_t _iml, Datum* _ppvar, Datum* _thread, double* _globals, NrnThread* _nt
+#define _internalthreadargsproto_ _nrn_mechanism_cache_range* _ml, size_t _iml, Datum* _ppvar, Datum* _thread, double* _globals, NrnThread* _nt
  	/*SUPPRESS 761*/
 	/*SUPPRESS 762*/
 	/*SUPPRESS 763*/
 	/*SUPPRESS 765*/
 	 extern double *hoc_getarg(int);
  
-#define t nrn_threads->_t
-#define dt nrn_threads->_dt
+#define t _nt->_t
+#define dt _nt->_dt
 #define tau_rec _ml->template fpfield<0>(_iml)
 #define tau_rec_columnindex 0
 #define tau_facil _ml->template fpfield<1>(_iml)
@@ -87,16 +84,16 @@ void _nrn_mechanism_register_data_fields(Args&&... args) {
 #define u_columnindex 8
 #define t_last _ml->template fpfield<9>(_iml)
 #define t_last_columnindex 9
-#define _tsav _ml->template fpfield<10>(_iml)
-#define _tsav_columnindex 10
+#define v _ml->template fpfield<10>(_iml)
+#define v_columnindex 10
+#define _tsav _ml->template fpfield<11>(_iml)
+#define _tsav_columnindex 11
 #define _nd_area *_ml->dptr_field<0>(_iml)
 #define wsyn	*_ppvar[2].get<double*>()
 #define _p_wsyn _ppvar[2].literal_value<void*>()
- static _nrn_mechanism_cache_instance _ml_real{nullptr};
-static _nrn_mechanism_cache_range *_ml{&_ml_real};
-static size_t _iml{0};
-static Datum *_ppvar;
+ /* Thread safe. No static _ml, _iml or _ppvar. */
  static int hoc_nrnpointerindex =  2;
+ static _nrn_mechanism_std_vector<Datum> _extcall_thread;
  /* external NEURON variables */
  /* declaration of user functions */
  static int _mechtype;
@@ -155,7 +152,6 @@ static void register_nmodl_text_and_filename(int mechtype);
  {"tau_syn", "ms"},
  {0, 0}
 };
- static double v = 0;
  /* connect global user variables to hoc */
  static DoubScal hoc_scdoub[] = {
  {0, 0}
@@ -166,10 +162,6 @@ static void register_nmodl_text_and_filename(int mechtype);
  static double _sav_indep;
  extern void _nrn_setdata_reg(int, void(*)(Prop*));
  static void _setdata(Prop* _prop) {
- neuron::legacy::set_globals_from_prop(_prop, _ml_real, _ml, _iml);
-_ppvar = _nrn_mechanism_access_dparam(_prop);
- Node * _node = _nrn_mechanism_access_node(_prop);
-v = _nrn_mechanism_access_voltage(_node);
  }
  static void _hoc_setdata(void* _vptr) { Prop* _prop;
  _prop = ((Point_process*)_vptr)->_prop;
@@ -219,7 +211,7 @@ static void nrn_alloc(Prop* _prop) {
      _nrn_mechanism_cache_instance _ml_real{_prop};
     auto* const _ml = &_ml_real;
     size_t const _iml{};
-    assert(_nrn_mechanism_get_num_vars(_prop) == 11);
+    assert(_nrn_mechanism_get_num_vars(_prop) == 12);
  	/*initialize range parameters*/
  	tau_rec = _parm_default[0]; /* 100 */
  	tau_facil = _parm_default[1]; /* 1000 */
@@ -227,7 +219,7 @@ static void nrn_alloc(Prop* _prop) {
  	u0 = _parm_default[3]; /* 0 */
  	tau_syn = _parm_default[4]; /* 2 */
   }
- 	 assert(_nrn_mechanism_get_num_vars(_prop) == 11);
+ 	 assert(_nrn_mechanism_get_num_vars(_prop) == 12);
  	_nrn_mechanism_access_dparam(_prop) = _ppvar;
  	/*connect ionic variables to this model*/
  
@@ -242,11 +234,11 @@ extern void hoc_register_tolerance(int, HocStateTolerance*, Symbol***);
 extern void _cvode_abstol( Symbol**, double*, int);
 
  extern "C" void _tsodyksmarkram_reg() {
-	int _vectorized = 0;
+	int _vectorized = 1;
   _initlists();
  	_pointtype = point_register_mech(_mechanism,
 	 nrn_alloc,nullptr, nullptr, nullptr, nrn_init,
-	 hoc_nrnpointerindex, 0,
+	 hoc_nrnpointerindex, 1,
 	 _hoc_create_pnt, _hoc_destroy_pnt, _member_func);
  _mechtype = nrn_get_mechtype(_mechanism[1]);
  hoc_register_parm_default(_mechtype, &_parm_default);
@@ -265,11 +257,12 @@ extern void _cvode_abstol( Symbol**, double*, int);
                                        _nrn_mechanism_field<double>{"z"} /* 7 */,
                                        _nrn_mechanism_field<double>{"u"} /* 8 */,
                                        _nrn_mechanism_field<double>{"t_last"} /* 9 */,
-                                       _nrn_mechanism_field<double>{"_tsav"} /* 10 */,
+                                       _nrn_mechanism_field<double>{"v"} /* 10 */,
+                                       _nrn_mechanism_field<double>{"_tsav"} /* 11 */,
                                        _nrn_mechanism_field<double*>{"_nd_area", "area"} /* 0 */,
                                        _nrn_mechanism_field<Point_process*>{"_pntproc", "pntproc"} /* 1 */,
                                        _nrn_mechanism_field<double*>{"wsyn", "pointer"} /* 2 */);
-  hoc_register_prop_size(_mechtype, 11, 3);
+  hoc_register_prop_size(_mechtype, 12, 3);
   hoc_register_dparam_semantics(_mechtype, 0, "area");
   hoc_register_dparam_semantics(_mechtype, 1, "pntproc");
   hoc_register_dparam_semantics(_mechtype, 2, "pointer");
@@ -291,8 +284,11 @@ static int _match_recurse=1;
 static void _modl_cleanup(){ _match_recurse=1;}
  
 static void _net_receive (Point_process* _pnt, double* _args, double _lflag) 
-{   neuron::legacy::set_globals_from_prop(_pnt->_prop, _ml_real, _ml, _iml);
-    _ppvar = _nrn_mechanism_access_dparam(_pnt->_prop);
+{  Prop* _p; Datum* _ppvar; Datum* _thread; NrnThread* _nt;
+   _nrn_mechanism_cache_instance _ml_real{_pnt->_prop};
+  auto* const _ml = &_ml_real;
+  size_t const _iml{};
+   _thread = nullptr; double* _globals = nullptr; _nt = (NrnThread*)_pnt->_vnt;   _ppvar = _nrn_mechanism_access_dparam(_pnt->_prop);
   if (_tsav > t){ hoc_execerror(hoc_object_name(_pnt->ob), ":Event arrived out of order. Must call ParallelContext.set_maxstep AFTER assigning minimum NetCon.delay");}
  _tsav = t; {
    z = z * exp ( - ( t - t_last ) / tau_rec ) ;
@@ -312,13 +308,18 @@ static void _net_receive (Point_process* _pnt, double* _args, double _lflag)
    } }
  
 static void _net_init(Point_process* _pnt, double* _args, double _lflag) {
-       _ppvar = _nrn_mechanism_access_dparam(_pnt->_prop);
+     _nrn_mechanism_cache_instance _ml_real{_pnt->_prop};
+  auto* const _ml = &_ml_real;
+  size_t const _iml{};
+  Datum* _ppvar = _nrn_mechanism_access_dparam(_pnt->_prop);
+  Datum* _thread = nullptr;
+  double* _globals = nullptr;
+  NrnThread* _nt = (NrnThread*)_pnt->_vnt;
  t_last = t ;
    }
 
-static void initmodel() {
-  int _i; double _save;_ninits++;
-{
+static void initmodel(_internalthreadargsproto_) {
+  int _i; double _save;{
  {
    y = 0.0 ;
    z = 0.0 ;
@@ -330,32 +331,43 @@ static void initmodel() {
 }
 
 static void nrn_init(_nrn_model_sorted_token const& _sorted_token, NrnThread* _nt, Memb_list* _ml_arg, int _type){
-Node *_nd; double _v; int* _ni; int _cntml;
 _nrn_mechanism_cache_range _lmr{_sorted_token, *_nt, *_ml_arg, _type};
 auto* const _vec_v = _nt->node_voltage_storage();
-_ml = &_lmr;
+auto* const _ml = &_lmr;
+Datum* _ppvar; Datum* _thread;
+Node *_nd; double _v; int* _ni; int _iml, _cntml;
 _ni = _ml_arg->_nodeindices;
 _cntml = _ml_arg->_nodecount;
+_thread = _ml_arg->_thread;
+double* _globals = nullptr;
+if (gind != 0 && _thread != nullptr) { _globals = _thread[_gth].get<double*>(); }
 for (_iml = 0; _iml < _cntml; ++_iml) {
  _ppvar = _ml_arg->_pdata[_iml];
  _tsav = -1e20;
    _v = _vec_v[_ni[_iml]];
  v = _v;
- initmodel();
-}}
+ initmodel(_threadargs_);
+}
+}
 
-static double _nrn_current(double _v){double _current=0.;v=_v;{
+static double _nrn_current(_internalthreadargsprotocomma_ double _v) {
+double _current=0.; v=_v;
+{
 } return _current;
 }
 
-static void nrn_state(_nrn_model_sorted_token const& _sorted_token, NrnThread* _nt, Memb_list* _ml_arg, int _type){
-Node *_nd; double _v = 0.0; int* _ni; int _cntml;
+static void nrn_state(_nrn_model_sorted_token const& _sorted_token, NrnThread* _nt, Memb_list* _ml_arg, int _type) {
 _nrn_mechanism_cache_range _lmr{_sorted_token, *_nt, *_ml_arg, _type};
 auto* const _vec_v = _nt->node_voltage_storage();
-_ml = &_lmr;
+auto* const _ml = &_lmr;
+Datum* _ppvar; Datum* _thread;
+Node *_nd; double _v = 0.0; int* _ni;
 _ni = _ml_arg->_nodeindices;
-_cntml = _ml_arg->_nodecount;
-for (_iml = 0; _iml < _cntml; ++_iml) {
+size_t _cntml = _ml_arg->_nodecount;
+_thread = _ml_arg->_thread;
+double* _globals = nullptr;
+if (gind != 0 && _thread != nullptr) { _globals = _thread[_gth].get<double*>(); }
+for (size_t _iml = 0; _iml < _cntml; ++_iml) {
  _ppvar = _ml_arg->_pdata[_iml];
  _nd = _ml_arg->_nodelist[_iml];
    _v = _vec_v[_ni[_iml]];
@@ -367,7 +379,7 @@ for (_iml = 0; _iml < _cntml; ++_iml) {
 
 static void terminal(){}
 
-static void _initlists() {
+static void _initlists(){
  int _i; static int _first = 1;
   if (!_first) return;
 _first = 0;
@@ -384,6 +396,7 @@ static void register_nmodl_text_and_filename(int mech_type) {
   "ENDCOMMENT\n"
   "\n"
   "NEURON {\n"
+  "    THREADSAFE\n"
   "    POINT_PROCESS TsodyksMarkramWA\n"
   "    RANGE tau_rec, tau_facil, U, u0, tau_syn\n"
   "    POINTER wsyn\n"

The key difference being:

 static void _net_init(Point_process* _pnt, double* _args, double _lflag) {
-       _ppvar = _nrn_mechanism_access_dparam(_pnt->_prop);
+     _nrn_mechanism_cache_instance _ml_real{_pnt->_prop};
+  auto* const _ml = &_ml_real;
+  size_t const _iml{};
+  Datum* _ppvar = _nrn_mechanism_access_dparam(_pnt->_prop);
+  Datum* _thread = nullptr;
+  double* _globals = nullptr;
+  NrnThread* _nt = (NrnThread*)_pnt->_vnt;
  t_last = t ;
    }

Since t_last is a macro equal to _ml->template fpfield<9>(_iml), I'm guessing the presence of THREADSAFE properly initializes _ml, and hence accessing the fpfield<9>(_iml) is well-behaved.

I consider this to be a workaround, a proper solution would be to fix NOCMODL, but for purposes of the ESD release, it should be satisfactory.

[nrnhines]

I've reviewed in nocmodl all the cases where the internal variable vectorize is set to 0 (the default value is 1) and in almost all cases it is being used in association with a notice about thread safety. I'm speculating that the fix is to keep the notice and merely remove the vectorize = 0; statements. All translated version 9 mod files are in fact vectorized in the sense of data layout (particularly RANGE) and the bug was due to legacy code paths when vectorize = 0 which failed to create or pass to functions the index of the mechanism instance.

To remind me of vectorize = 0 usage, here are the notes from assessing src/nocmodl$ grep -r -n -w vectorize and looking at the assignment statements

vectorize = 0
all of the following also print a thread safe notice

VERBATIM
call to EXTDEF5 function
DERIVATIVE and not cnexp, derivimplicit, euler
NONLINEAR
LINEAR  
DISCRETE
SOLVE with ERROR
SUFFIX nothing
artificial_cellin some circumstances
state_discontinuity outside NET_RECEIVE

calls to threadsafe(const char* s)
EXTERNAL
POINTER
BBCOREPOINTER
Assignment to global variable

[nrnhines]

Ironically, assignment to the t_last RANGE variable in the NET_RECEIVE INITIAL block abuses a bit the purpose of that INITIAL block which is intended to initialize only the arguments (after the first) of the NET_RECEIVE block. Really, t_last should be initialized in the top level mod file's INITIAL block. As is, the t_last instance will be initialized for however many NetCon's target this mechanism instance.