From 1ac4668804a931c1e909a4aa2510f1d92a669710 Mon Sep 17 00:00:00 2001 From: Jorge Blanco Alonso Date: Tue, 27 Feb 2024 15:45:26 +0100 Subject: [PATCH] Move CompartmentMapping to python --- core/hoc/CompartmentMapping.hoc | 285 -------------------- core/hoc/mapping.py | 1 - core/python/mapping.py | 20 -- neurodamus/core/_neurodamus.py | 1 - neurodamus/core/coreneuron_configuration.py | 98 +++++++ neurodamus/node.py | 6 +- 6 files changed, 101 insertions(+), 310 deletions(-) delete mode 100644 core/hoc/CompartmentMapping.hoc delete mode 120000 core/hoc/mapping.py delete mode 100644 core/python/mapping.py diff --git a/core/hoc/CompartmentMapping.hoc b/core/hoc/CompartmentMapping.hoc deleted file mode 100644 index 17a73dba..00000000 --- a/core/hoc/CompartmentMapping.hoc +++ /dev/null @@ -1,285 +0,0 @@ -/** - * @file CompartmentMapping.hoc - * @brief Implement interface to register section segment mapping with NEURON. - * @remark Copyright © BBP/EPFL 2005-2017; All rights reserved. Do not distribute without further notice. - */ - -{nrnpython( "import mapping" )} - -// global object which python will manipulate -objref secvec -objref segvec -objref lfp_factors -objref all_lfp_factors - -/*! - * Get the number of sections of certain type - * @param $o1 Cell object who owns the currently accessed section - * @param $s2 Section name - * @return number of sections of certain type - */ -func get_num_sec() { local num_sections localobj strobj, cellref - strdef section_name - strobj = new StringFunctions() - cellref = $o1 - section_name = $s2 - num_sections = 0 - if(strcmp( section_name, "ais" ) == 0) { - result = execute1( "{nSecLastAIS}", cellref, 0 ) - if( result != 0 ) { - num_sections = cellref.nSecLastAIS - } - } else if(strcmp( section_name, "node" ) == 0) { - result = execute1( "{nSecNodal}", cellref, 0 ) - if( result != 0 ) { - num_sections = cellref.nSecNodal - } - } - return num_sections -} - -/*! - * Determine the absolute index of a section within the cell - * @param $o1 Cell object who owns the currently accessed section - * @return index which can be used to map back onto the neuron - */ -func sectionNo() { local baseOffset, secIndex localobj str, strobj - strdef subset, section_name - str = new String() - strobj = new StringFunctions() - - if( strobj.substr( secname(), "soma" ) > 0 ) { - subset = "soma" - baseOffset = 0 - } else if (strobj.substr( secname(), "axon") >0) { - subset = "axon" - baseOffset = $o1.nSecSoma - }else if (strobj.substr( secname(), "dend") >0) { - subset = "dend" - baseOffset = $o1.nSecSoma + $o1.nSecAxonalOrig - }else if (strobj.substr( secname(), "apic") > 0) { - subset = "apic" - baseOffset = $o1.nSecSoma + $o1.nSecAxonalOrig + $o1.nSecBasal - } else if (strobj.substr( secname(), "ais") >0) { - subset = "ais" - baseOffset = $o1.nSecSoma + $o1.nSecAxonalOrig + $o1.nSecBasal + $o1.nSecApical - } else if (strobj.substr( secname(), "node") >0) { - subset = "node" - baseOffset = $o1.nSecSoma + $o1.nSecAxonalOrig + $o1.nSecBasal + $o1.nSecApical + get_num_sec($o1, "ais") - } else if (strobj.substr(secname(), "myelin") >0) { - subset = "myelin" - baseOffset = $o1.nSecSoma + $o1.nSecAxonalOrig + $o1.nSecBasal + $o1.nSecApical + get_num_sec($o1, "ais") + get_num_sec($o1, "node") - } - - // parse out the index from the secname - strobj.tail( secname(), subset, str.s ) - if (strobj.len(str.s) > 0) { - strobj.left(str.s, strobj.len(str.s)-1) // result should look like [n - strobj.right(str.s, 1) - sscanf( str.s, "%d", &secIndex ) - } else { - // todo : abort on error here - print "Error while getting section number (should this ever happen?) !" - } - - return baseOffset + secIndex -} - -/*! - * Register section/segment mapping to NEURON for report generation - * todo : removing logic of target selection as its incomplete - * @param $o1 cellDistributor object - * @param $o2 target object - */ - -proc registerMapping () { local cellIndex, activegid, nSoma, nAxon, nDend, nApic, nSec, indexSec, nSeg, val, num_electrodes \ - localobj gidvec, cellref, pyobj, pop_info, commstr, target, cellDistributor, pc, python_obj - strdef tstr - - cellDistributor = $o1 - //target = $o2 - - pyobj = new PythonObject() - commstr = new String() - pc = new ParallelContext() - - // list of gids on current rank - gidvec = cellDistributor.getGidListForProcessor() - pop_info = new PythonObject() - - for cellIndex=0, gidvec.size()-1 { - activegid = gidvec.x[cellIndex] - cellref = cellDistributor.getCell(activegid) - pop_info = cellDistributor.getPopulationInfo(activegid) - - // skip gids that we don't want to report, instead of skipping - // we should add a flag in mapping - //if(!target.contains(activegid)) { - //continue - //} - - nSoma = 0 - nAxon = 0 - nDend = 0 - nApic = 0 - nAis = 0 - nNode = 0 - nMyel = 0 - nSec = 0 - indexSec = 0 - - // vectors holding section and segments - secvec = new Vector() - segvec = new Vector() - lfp_factors = new Vector() - all_lfp_factors = new Vector() - - // gather all somas and register somas with neuron - forsec cellref.somatic { - nSec = nSec + 1 - val = sectionNo( cellref ) - sprint( commstr.s, "mappinginfo = mapping.getNodeID(%d, %d)", val, nSoma ) - nrnpython( commstr.s ) - nSoma = pyobj.mappinginfo - } - - num_electrodes = 0 - python_obj = new PythonObject() - has_lfp_manager = python_obj.hasattr(cellDistributor, "_lfp_manager") - if( has_lfp_manager != 0 ) { - num_electrodes = cellDistributor._lfp_manager.get_number_electrodes(activegid, pop_info) - all_lfp_factors = cellDistributor._lfp_manager.read_lfp_factors(activegid, pop_info) - if ( all_lfp_factors.size() > 0 && nSoma != 0 ) { - // Filter the factors of each section - lfp_factors = all_lfp_factors.c(indexSec*num_electrodes, (indexSec+nSoma)*num_electrodes-1) - } - } - pc.nrnbbcore_register_mapping(activegid, "soma", secvec, segvec, lfp_factors, num_electrodes ) - indexSec = indexSec + nSoma - // gather axonal compartments and register with neuron - - nSec = 0 - secvec = new Vector() - segvec = new Vector() - lfp_factors = new Vector() - - forsec cellref.axonal { - nSec = nSec + 1 - val = sectionNo( cellref ) - sprint( commstr.s, "mappinginfo = mapping.getNodeID( %d, %d )", val, nAxon ) - nrnpython( commstr.s ) - nAxon = pyobj.mappinginfo - } - if( has_lfp_manager != 0 && all_lfp_factors.size() > 0 && nAxon != 0 ) { - lfp_factors = all_lfp_factors.c(indexSec*num_electrodes, (indexSec+nAxon)*num_electrodes-1) - } - pc.nrnbbcore_register_mapping(activegid, "axon", secvec, segvec, lfp_factors, num_electrodes ) - indexSec = indexSec + nAxon - // gather basal compartments and register with neuron - - nSec = 0 - secvec = new Vector() - segvec = new Vector() - lfp_factors = new Vector() - - forsec cellref.basal { - nSec = nSec + 1 - val = sectionNo( cellref ) - sprint( commstr.s, "mappinginfo = mapping.getNodeID( %d, %d )", val, nDend ) - nrnpython( commstr.s ) - nDend = pyobj.mappinginfo - } - if( has_lfp_manager != 0 && all_lfp_factors.size() > 0 && nDend != 0 ) { - lfp_factors = all_lfp_factors.c(indexSec*num_electrodes, (indexSec+nDend)*num_electrodes-1) - } - pc.nrnbbcore_register_mapping(activegid, "dend", secvec, segvec, lfp_factors, num_electrodes ) - indexSec = indexSec + nDend - - nSec = 0 - secvec = new Vector() - segvec = new Vector() - lfp_factors = new Vector() - - // gather apical compartments and register with neuron - - forsec cellref.apical { - nSec = nSec + 1 - val = sectionNo( cellref ) - sprint( commstr.s, "mappinginfo = mapping.getNodeID( %d, %d )", val, nApic ) - nrnpython( commstr.s ) - nApic = pyobj.mappinginfo - } - - if( has_lfp_manager != 0 && all_lfp_factors.size() > 0 && nApic != 0 ) { - lfp_factors = all_lfp_factors.c(indexSec*num_electrodes, (indexSec+nApic)*num_electrodes-1) - } - pc.nrnbbcore_register_mapping(activegid, "apic", secvec, segvec, lfp_factors, num_electrodes ) - indexSec = indexSec + nApic - - nSec = 0 - secvec = new Vector() - segvec = new Vector() - lfp_factors = new Vector() - - result = execute1( "{AIS}", cellref, 0 ) - if( result != 0 ) { - forsec cellref.AIS { - nSec = nSec + 1 - val = sectionNo( cellref ) - sprint( commstr.s, "mappinginfo = mapping.getNodeID( %d, %d )", val, nAis ) - nrnpython( commstr.s ) - nAis = pyobj.mappinginfo - } - if( has_lfp_manager != 0 && all_lfp_factors.size() > 0 && nAis != 0 ) { - lfp_factors = all_lfp_factors.c(indexSec*num_electrodes, (indexSec+nAis)*num_electrodes-1) - } - pc.nrnbbcore_register_mapping(activegid, "ais", secvec, segvec, lfp_factors, num_electrodes ) - indexSec = indexSec + nAis - } - - // gather nodal compartments and register with neuron - nSec = 0 - secvec = new Vector() - segvec = new Vector() - lfp_factors = new Vector() - - result = execute1( "{nodal}", cellref, 0 ) - if( result != 0 ) { - forsec cellref.nodal { - nSec = nSec + 1 - val = sectionNo( cellref ) - sprint( commstr.s, "mappinginfo = mapping.getNodeID( %d, %d )", val, nNode ) - nrnpython( commstr.s ) - nNode = pyobj.mappinginfo - } - if( has_lfp_manager != 0 && all_lfp_factors.size() > 0 && nNode != 0 ) { - lfp_factors = all_lfp_factors.c(indexSec*num_electrodes, (indexSec+nNode)*num_electrodes-1) - } - pc.nrnbbcore_register_mapping(activegid, "node", secvec, segvec, lfp_factors, num_electrodes ) - indexSec = indexSec + nNode - } - - // gather myelinated compartments and register with neuron - nSec = 0 - secvec = new Vector() - segvec = new Vector() - lfp_factors = new Vector() - - result = execute1( "{myelinated}", cellref, 0 ) - if( result != 0 ) { - forsec cellref.myelinated { - nSec = nSec + 1 - val = sectionNo( cellref ) - sprint( commstr.s, "mappinginfo = mapping.getNodeID( %d, %d )", val, nMyel ) - nrnpython( commstr.s ) - nMyel = pyobj.mappinginfo - } - if( has_lfp_manager != 0 && all_lfp_factors.size() > 0 && nMyel != 0 ) { - lfp_factors = all_lfp_factors.c(indexSec*num_electrodes, (indexSec+nMyel)*num_electrodes-1) - } - pc.nrnbbcore_register_mapping(activegid, "myelin", secvec, segvec, lfp_factors, num_electrodes ) - indexSec = indexSec + nMyel - } - - } // end for gidIndex -} diff --git a/core/hoc/mapping.py b/core/hoc/mapping.py deleted file mode 120000 index 6ff3e6ec..00000000 --- a/core/hoc/mapping.py +++ /dev/null @@ -1 +0,0 @@ -../python/mapping.py \ No newline at end of file diff --git a/core/python/mapping.py b/core/python/mapping.py deleted file mode 100644 index df4de299..00000000 --- a/core/python/mapping.py +++ /dev/null @@ -1,20 +0,0 @@ -from neuron import h - -pc = h.ParallelContext() -nhost = int(pc.nhost()) -rank = int(pc.id()) - - -def getNodeID(sectionID, nSeg): - """ - Since hoc has no interface to access section node info (node_index), - provide a python interface. - """ - sec = h.cas() - - for seg in sec: - h.secvec.append(sectionID) - h.segvec.append(seg.node_index()) - nSeg += 1 - - return nSeg diff --git a/neurodamus/core/_neurodamus.py b/neurodamus/core/_neurodamus.py index 0679a824..a66f79f8 100644 --- a/neurodamus/core/_neurodamus.py +++ b/neurodamus/core/_neurodamus.py @@ -58,7 +58,6 @@ def _init(cls, **kwargs): # Additional libraries introduced in saveUpdate sys.path.append(os.environ['HOC_LIBRARY_PATH']) - cls.load_hoc("CompartmentMapping") # Attempt to instantiate BBSaveState to early detect errors cls._h.BBSaveState() diff --git a/neurodamus/core/coreneuron_configuration.py b/neurodamus/core/coreneuron_configuration.py index e987c685..59a9283f 100644 --- a/neurodamus/core/coreneuron_configuration.py +++ b/neurodamus/core/coreneuron_configuration.py @@ -5,6 +5,104 @@ from . import NeurodamusCore as Nd +class CompartmentMapping: + """ Interface to register section segment mapping with NEURON. + """ + def __init__(self, cell_distributor): + self.cell_distributor = cell_distributor + self.pc = Nd.ParallelContext() + + def get_section_index(self, cell, section): + """ + Calculate the global index of a given section within its cell. + :param cell: The cell instance containing the section of interest + :param section: The specific section for which the index is required + :return: The global index of the section, applicable for neuron mapping + """ + section_name = str(section) + base_offset = 0 + section_index = 0 + if "soma" in section_name: + pass # base_offset is 0 + elif "axon" in section_name: + base_offset = cell.nSecSoma + elif "dend" in section_name: + base_offset = cell.nSecSoma + cell.nSecAxonalOrig + elif "apic" in section_name: + base_offset = cell.nSecSoma + cell.nSecAxonalOrig + cell.nSecBasal + elif "ais" in section_name: + base_offset = cell.nSecSoma + cell.nSecAxonalOrig + cell.nSecBasal + cell.nSecApical + elif "node" in section_name: + base_offset = cell.nSecSoma + cell.nSecAxonalOrig + cell.nSecBasal + cell.nSecApical \ + + getattr(cell, 'nSecLastAIS', 0) + elif "myelin" in section_name: + base_offset = cell.nSecSoma + cell.nSecAxonalOrig + cell.nSecBasal + cell.nSecApical \ + + getattr(cell, 'nSecLastAIS', 0) + getattr(cell, 'nSecNodal', 0) + + # Extract the index from the section name + try: + index_str = section_name.split('[')[-1].rstrip(']') + section_index = int(index_str) + except ValueError: + logging.warning(f"Error while getting section index {index_str}") + + return int(base_offset + section_index) + + def create_section_vectors(self, section_id, section, secvec, segvec): + num_segments = 0 + for seg in section: + secvec.append(section_id) + segvec.append(seg.node_index()) + num_segments += 1 + + return num_segments + + def process_section(self, cell, sections, num_electrodes, all_lfp_factors, section_offset): + secvec, segvec, lfp_factors = Nd.Vector(), Nd.Vector(), Nd.Vector() + num_segments = 0 + section_attr = getattr(cell, sections[0], None) + if section_attr: + for sec in section_attr: + section_index = self.get_section_index(cell, sec) + num_segments += self.create_section_vectors(section_index, sec, secvec, segvec) + + if num_electrodes > 0 and all_lfp_factors.size() > 0 and num_segments > 0: + start_idx = section_offset * num_electrodes + end_idx = (section_offset + num_segments) * num_electrodes - 1 + lfp_factors.copy(all_lfp_factors, start_idx, end_idx) + + self.pc.nrnbbcore_register_mapping(cell.gid, sections[1], secvec, segvec, + lfp_factors, num_electrodes) + return num_segments + + def register_mapping(self): + sections = [ + ('somatic', 'soma'), + ('axonal', 'axon'), + ('basal', 'dend'), + ('apical', 'apic'), + ('AIS', 'ais'), + ('nodal', 'node'), + ('myelinated', 'myelin') + ] + gidvec = self.cell_distributor.getGidListForProcessor() + for activegid in gidvec: + cellref = self.cell_distributor.getCell(activegid) + all_lfp_factors = Nd.Vector() + num_electrodes = 0 + lfp_manager = getattr(self.cell_distributor, "_lfp_manager", None) + if lfp_manager: + pop_info = self.cell_distributor.getPopulationInfo(activegid) + num_electrodes = lfp_manager.get_number_electrodes(activegid, pop_info) + all_lfp_factors = lfp_manager.read_lfp_factors(activegid, pop_info) + + section_offset = 0 + for section in sections: + processed_segments = self.process_section(cellref, section, num_electrodes, + all_lfp_factors, section_offset) + section_offset += processed_segments + + class _CoreNEURONConfig(object): """ The CoreConfig class is responsible for managing the configuration of the CoreNEURON simulation. diff --git a/neurodamus/node.py b/neurodamus/node.py index 3d28e620..1174bd98 100644 --- a/neurodamus/node.py +++ b/neurodamus/node.py @@ -33,7 +33,7 @@ from .utils.logging import log_stage, log_verbose, log_all from .utils.memory import DryRunStats, trim_memory, pool_shrink, free_event_queues, print_mem_usage from .utils.timeit import TimerManager, timeit -from .core.coreneuron_configuration import CoreConfig +from .core.coreneuron_configuration import CoreConfig, CompartmentMapping from .io.sonata_config import ConnectionTypes # Internal Plugins from . import ngv as _ngv # NOQA @@ -1181,7 +1181,7 @@ def _coreneuron_ensure_all_ranks_have_gids(self, corenrn_data): base_manager.load_artificial_cell(fake_gid, CoreConfig.artificial_cell_object) yield - # Nd.registerMapping doesn't work for this artificial cell as somatic attr is + # register_mapping() doesn't work for this artificial cell as somatic attr is # missing, so create a dummy mapping file manually, required for reporting cur_files = glob.glob("%s/*_3.dat" % corenrn_data) example_mapfile = cur_files[0] @@ -1259,7 +1259,7 @@ def _sim_corenrn_write_config(self, corenrn_restore=False): fwd_skip = self._run_conf.get("ForwardSkip", 0) if not corenrn_restore else 0 if not corenrn_restore: - Nd.registerMapping(self._circuits.global_manager) + CompartmentMapping(self._circuits.global_manager).register_mapping() with self._coreneuron_ensure_all_ranks_have_gids(CoreConfig.datadir): self._pc.nrnbbcore_write(CoreConfig.datadir) MPI.barrier() # wait for all ranks to finish corenrn data generation