bugfix: Stop rxd.Region keeping sections alive when they go out of sc…

…ope. (#1072)

* bugfix: Stop rxd.Region keeping sections alive when they go out of scope.

* Update test_region_sections.py

* Prevent species from keeping sections alive.

* Updated following @ramcdougal review.

* Remove deleted section from rxd.

* Additional check and Simplify test.

share/lib/python/neuron/rxd/region.py

@@ -29,7 +29,12 @@ def _sort_secs(secs):
     #for sec in secs:
     #    if sec.orientation():
     #        raise RxDException('still need to deal with backwards sections')
-    return [secs_names[sec.hoc_internal_name()] for sec in all_sorted if sec.hoc_internal_name() in secs_names]
+    secs = [secs_names[sec.hoc_internal_name()] for sec in all_sorted if sec.hoc_internal_name() in secs_names]
+    # return an empty list rather than an empty SectionList because
+    # bool([]) == False
+    # bool(h.SectionList([])) == True
+    # secs are checked in some reactions to determine active regions
+    return [] if secs == [] else h.SectionList(secs)
 
 
 class _c_region:
@@ -41,7 +46,7 @@ class _c_region:
     def __init__(self, regions):
         global _c_region_lookup
         self._regions = [weakref.ref(r) for r in regions]
-        self._overlap = h.SectionList(self._regions[0]()._secs1d)
+        self._overlap = self._regions[0]()._secs1d
         self.num_regions = len(self._regions)
         self.num_species = 0
         self.num_params = 0
@@ -510,15 +515,16 @@ def __init__(self, secs=None, nrn_region=None, geometry=None, dimension=None, dx
         """
         self._allow_setting = True
         if hasattr(secs,'__len__'):
-            self.secs = secs
+            self._secs = secs
         else:
-            self.secs = [secs]
+            self._secs = [secs]
         if secs == [] or secs == None:
             warnings.warn("Warning: No sections. Region 'secs' should be a list of NEURON sections.")
         from nrn import Section
-        for sec in self.secs:
+        for sec in self._secs:
             if not isinstance(sec,Section):
                 raise RxDException("Error: Region 'secs' must be a list of NEURON sections, %r is not a valid NEURON section." % sec)
+        self._secs = h.SectionList(self._secs)
         self.nrn_region = nrn_region
         self.geometry = geometry
 
@@ -620,7 +626,16 @@ def secs(self):
     @secs.setter
     def secs(self, value):
         if hasattr(self, '_allow_setting'):
-            self._secs = value
+            if hasattr(value, '__len__'):
+                secs = value
+            else:
+                secs = [value]
+
+            from nrn import Section
+            for sec in secs:
+                if not isinstance(sec,Section):
+                    raise RxDException("Error: Region 'secs' must be a list of NEURON sections, %r is not a valid NEURON section." % sec)
+            self._secs = h.SectionList(secs)
         else:
             raise RxDException('Cannot set secs now; model already instantiated')
     def volume(self, index):

share/lib/python/neuron/rxd/section1d.py

@@ -18,30 +18,45 @@ def _donothing(): pass
 class _SectionLookup:
     class Lookup:
         def __init__(self):
-            self.rxd_sec_lookup = {}
+            self.rxd_sec_list = {} 
+            self.nrn_sec_list = h.SectionList()
     _instance = None
 
     def __init__(self):
         if _SectionLookup._instance is None:
            _SectionLookup._instance = _SectionLookup.Lookup()
 
     def __getitem__(self, key):
-        return _SectionLookup._instance.rxd_sec_lookup[key]
+        if (key in _SectionLookup._instance.nrn_sec_list and
+            key.hoc_internal_name() in _SectionLookup._instance.rxd_sec_list):
+            return _SectionLookup._instance.rxd_sec_list[key.hoc_internal_name()]
+        return []
 
     def __setitem__(self, key, value):
-        _SectionLookup._instance.rxd_sec_lookup[key] = value
+        _SectionLookup._instance.nrn_sec_list.append(key)
+        _SectionLookup._instance.rxd_sec_list[key.hoc_internal_name()] = value
 
     def __delitem__(self, key):
-        del _SectionLookup._instance.rxd_sec_lookup[key]
+        _SectionLookup._instance.nrn_sec_list.remove(key)
+        if key.hoc_internal_name() in _SectionLookup._instance.rxd_sec_list:
+            del _SectionLookup._instance.rxd_sec_list[key.hoc_internal_name()]
 
     def __iter__(self):
-        return iter(_SectionLookup._instance.rxd_sec_lookup)
+        return iter(_SectionLookup._instance.nrn_sec_list)
 
     def values(self):
-        return _SectionLookup._instance.rxd_sec_lookup.values()
+        return _SectionLookup._instance.rxd_sec_list.values()
 
     def items(self):
-        return _SectionLookup._instance.rxd_sec_lookup.items()
+        res = {}
+        for sec in _SectionLookup._instance.nrn_sec_list:
+            res[sec] = _SectionLookup._instance.rxd_sec_list[sec.hoc_internal_name()]
+        return res.items()
+
+    def remove(self, rxdsec):
+        for key,val in _SectionLookup._instance.rxd_sec_list.items():
+            if val == rxdsec:
+                del _SectionLookup._instance.rxd_sec_list[key]
 
 
 
@@ -58,21 +73,23 @@ def add_values(mat, i, js, vals):
 
 def _parent(sec):
     """Return the parent of seg or None if sec is a root"""
-    seg = sec.trueparentseg()
-    if seg:
-        return seg
-    else:
-        seg = sec.parentseg()
-    if seg:
-        par = seg.sec
-        parseg = par.parentseg()
-        if parseg and seg.x == par.orientation:
-            # connection point belongs to temp's ancestor
-            return _parent(seg.sec)
-        else:
+    if sec:
+        seg = sec.trueparentseg()
+        if seg:
             return seg
-    else:
-        return None
+        else:
+            seg = sec.parentseg()
+        if seg:
+            par = seg.sec
+            parseg = par.parentseg()
+            if parseg and seg.x == par.orientation:
+                # connection point belongs to temp's ancestor
+                return _parent(seg.sec)
+            else:
+                return seg
+        else:
+            return None
+    return None
 
 def _purge_cptrs():
     """purges all cptr information"""
@@ -127,7 +144,7 @@ class Section1D(rxdsection.RxDSection):
     def __init__(self, species, sec, diff, r):
         self._species = weakref.ref(species)
         self._diff = diff
-        self._sec = sec
+        self._secref = h.SectionList([sec])
         self._concentration_ptrs = None
         self._offset = node._allocate(sec.nseg + 1)
         self._nseg = sec.nseg
@@ -138,6 +155,10 @@ def __init__(self, species, sec, diff, r):
             self._rxd_sec_lookup[sec].append(self)
         else:
             self._rxd_sec_lookup[sec] = [self]
+    @property
+    def _sec(self):
+        sl = list(self._secref)
+        return None if not sl else sl[0]
 
     def __req__(self, other):
         if isinstance(other, nrn.Section):
@@ -162,12 +183,15 @@ def _init_diffusion_rates(self):
         node._diffs[self._offset : self._offset + self.nseg] = self._diff
 
     def _update_node_data(self):
-        nseg_changed = 0
+        nseg_changed = False
+        if self._sec is None:
+            self._delete()
+            return True
         if self._nseg != self._sec.nseg:
             num_roots = self._species()._num_roots
             offset = node._allocate(self._sec.nseg + num_roots + 1)
             replace(self, offset, self._sec.nseg)
-            nseg_changed = 1
+            nseg_changed = True
         volumes, surface_area, diffs = node._get_data()
         geo = self._region._geometry
         volumes[self._offset : self._offset + self._nseg] = geo.volumes1d(self)
@@ -183,32 +207,34 @@ def _delete(self):
         _rxd_sec_lookup = _SectionLookup()
 
         # remove ref to this section -- at exit weakref.ref might be none 
-        if self._sec in _rxd_sec_lookup:
-            sec_list = _rxd_sec_lookup[self._sec]
-            sec_list = list(filter(lambda x: not x == self, sec_list))
-            if sec_list == []:
-                del _rxd_sec_lookup[self._sec]
-            else:
-                _rxd_sec_lookup[self._sec] = sec_list
+        if self._sec:
+            if self._sec in _rxd_sec_lookup:
+                sec_list = [s for s in _rxd_sec_lookup[self._sec] if s != self]
+                if sec_list == []:
+                    del _rxd_sec_lookup[self._sec]
+                else:
+                    _rxd_sec_lookup[self._sec] = sec_list
+        else:
+            _rxd_sec_lookup.remove(self)
 
-            # check if memory has already been free'd
-            if self._nseg < 0:
-                return
+        # check if memory has already been free'd
+        if self._nseg < 0:
+             return
 
-            # node data is removed here in case references to sec remains
-            if hasattr(self, "_num_roots"):
-                self._nseg += self._num_roots
-                self._offset -= self._num_roots
-            node._remove(self._offset, self._offset + self._nseg + 1)
+        # node data is removed here in case references to sec remains
+        if hasattr(self, "_num_roots"):
+            self._nseg += self._num_roots
+            self._offset -= self._num_roots
+        node._remove(self._offset, self._offset + self._nseg + 1)
 
-            # shift offset to account for deleted sec
-            for secs in _rxd_sec_lookup.values():
-                for s in secs:
-                    if s._offset > self._offset:
-                        s._offset -= self._nseg + 1
+        # shift offset to account for deleted sec
+        for secs in _rxd_sec_lookup.values():
+            for s in secs:
+                if s._offset > self._offset:
+                    s._offset -= self._nseg + 1
 
-            # mark memory as free'd
-            self._nseg = -1
+        # mark memory as free'd
+        self._nseg = -1
 
     def __del__(self):
         self._delete()
@@ -265,6 +291,10 @@ def _register_cptrs(self):
             self._concentration_ptrs = []
 
     def _setup_diffusion_matrix(self, mat):
+        if not self._sec:
+            if self._nseg > 0:
+                self._delete()
+            return
         _volumes, _surface_area, _diffs = node._get_data()
         offset = self._offset
         dx = self.L / self.nseg

share/lib/python/neuron/rxd/species.py

@@ -1525,6 +1525,8 @@ def _update_node_data(self):
         nsegs_changed = 0
         for sec in self._secs:
             nsegs_changed += sec._update_node_data()
+        # remove deleted sections
+        self._secs = [sec for sec in self._secs if sec and sec._nseg > 0]
         return nsegs_changed
 
     def concentrations(self):

test/rxd/test_region_sections.py

@@ -0,0 +1,28 @@
+import pytest
+
+def test_region_sections(neuron_instance):
+    """A test that regions and species do not keep sections alive"""
+
+    h, rxd, data, save_path = neuron_instance
+    soma = h.Section(name='soma')
+    cyt = rxd.Region(soma.wholetree(), name='cyt')
+    ca = rxd.Species(cyt, name='ca', charge=2)
+    h.finitialize(-65)
+    soma = None
+
+    for sec in h.allsec():
+        assert(False)
+
+def test_section_removal(neuron_instance):
+    """A test that deleted sections are removed from rxd"""
+
+    h, rxd, data, save_path = neuron_instance
+    soma = h.Section(name='soma')
+    cyt = rxd.Region(soma.wholetree(), name='cyt')
+    ca = rxd.Species(cyt, name='ca', charge=2)
+    h.finitialize(-65)
+
+    soma = h.Section(name='soma')
+    h.finitialize(-65)
+
+    assert(not ca.nodes)