From 9113be2d604999a768bbbcbda2ee1743b30dca76 Mon Sep 17 00:00:00 2001 From: Jorge Blanco Alonso <41900536+jorblancoa@users.noreply.github.com> Date: Fri, 9 Feb 2024 17:59:14 +0100 Subject: [PATCH] [BBPBGLIB-1124] Remove hoc targets related code (#118) ## Context Addresses BBPBGLIB-1124. Removes code related to hoc targets, as well as the need of using a hoc TargetManager. ## Scope Few functions were implemented in the python Targetmanager in order to remove the hoc version: getMETypes(), getPointList() and location_to_point() ## Review * [x] PR description is complete * [x] Coding style (imports, function length, New functions, classes or files) are good * [ ] Unit/Scientific test added * [ ] Updated Readme, in-code, developer documentation --- core/hoc/SerializedSections.hoc | 82 -- core/hoc/Target.hoc | 964 ------------------ core/hoc/TargetManager.hoc | 720 ------------- core/hoc/TargetParser.hoc | 482 --------- core/hoc/neurodamus.hoc | 2 - examples/test_neurodamus.py | 8 +- neurodamus/cell_distributor.py | 6 - neurodamus/connection.py | 3 +- neurodamus/connection_manager.py | 10 +- neurodamus/modification_manager.py | 13 +- neurodamus/node.py | 72 +- neurodamus/stimulus_manager.py | 2 +- neurodamus/target_manager.py | 371 +++---- tests/integration-e2e/test_multicycle_runs.py | 17 - tests/unit/test_target_intersection.py | 29 +- 15 files changed, 164 insertions(+), 2617 deletions(-) delete mode 100644 core/hoc/SerializedSections.hoc delete mode 100644 core/hoc/Target.hoc delete mode 100644 core/hoc/TargetManager.hoc delete mode 100644 core/hoc/TargetParser.hoc diff --git a/core/hoc/SerializedSections.hoc b/core/hoc/SerializedSections.hoc deleted file mode 100644 index 86abd01c..00000000 --- a/core/hoc/SerializedSections.hoc +++ /dev/null @@ -1,82 +0,0 @@ -/** - * @file SerializedSections.hoc - * @brief Allow morphology sections to be accessed from an array by index. - * @author king - * @date 2009-06-12 - * @remark Copyright © BBP/EPFL 2005-2011; All rights reserved. Do not distribute without further notice. - */ - -//make sure we have some basic things loaded -{load_file( "stdrun.hoc" )} -{load_file( "Map.hoc" )} -{load_file("fileUtils.hoc")} - -//------------------------------------------------------------------------------------------ -// -//------------------------------------------------------------------------------------------ - -_serialized_sections_warned = 0 - -/** - * When getting a location to a compartment by section index, the compartments must be serialized - * and stored for future access. - */ -begintemplate SerializedSections - -external _serialized_sections_warned, terminate - -//----------------------------------------------------------------------------------------------- -// Declare member variables -//----------------------------------------------------------------------------------------------- - -objref isec2sec[1] - -//----------------------------------------------------------------------------------------------- -// public members -//----------------------------------------------------------------------------------------------- - -public isec2sec, init, count - -//----------------------------------------------------------------------------------------------- -// Member function implementations -//----------------------------------------------------------------------------------------------- - -/*! - * Constructor serializes the sections of a cell for easier random access. Note that this is possible because - * the v field in the section has been assigned an integer corresponding to the target index as read from the - * morphology file. - * - * $o1 Cell reference (basic cell, not a CCell) - */ -proc init() { local index - - n = $o1.nSecAll - - objref isec2sec[n] - - index=1 - forsec $o1.all { - if( v(0.0001) >= n ) { - printf("%s v(1)=%d n3d()=%d\n", secname(), v(1), n3d()) - terminate("Error: failure in mk2_isec2sec()") - } - if( v(0.0001) < 0 ) { - if( _serialized_sections_warned == 0 ) { - print "[Warning] SerializedSections: v(0.0001) < 0. index=", index, " v()=", v(0.0001) - _serialized_sections_warned = 1 // Dont report more - } - } else { - isec2sec[v(0.0001)] = new SectionRef() - } - index = index+1 - } -} - -/*! - * @return The number of sections the cell originally had - */ -func count() { - return n -} - -endtemplate SerializedSections diff --git a/core/hoc/Target.hoc b/core/hoc/Target.hoc deleted file mode 100644 index 39f4e63f..00000000 --- a/core/hoc/Target.hoc +++ /dev/null @@ -1,964 +0,0 @@ -/** - * @file Target.hoc - * @brief Encapsulate the operations surrounding a conceptual target used by the Blue Brain Project to identify neurons or components of neurons - * @brief for stimulation, reporting, or other operations - * @author king - * @date 2009-06-12 - * @remark Copyright © BBP/EPFL 2005-2011; All rights reserved. Do not distribute without further notice. - */ - -{load_file("stdlib.hoc")} -{load_file( "TPointList.hoc" )} -{load_file("fileUtils.hoc")} - - -// The whole set of gids local to this node -objref local_gids -local_gids = new Vector() - - -/*! - * TargetUtility contains a collection of functions that are useful in the creation/maintenence/usage of targets - * such as finding a target in a list of target objects sorted by name. Also included is a basic binary search function - * that finds an int (gid) in a vector (gidvec) - */ -begintemplate TargetUtility - -external terminate - -public findTarget, binary_search_contains - -/*! - * Utility function for searching a List of targets objects and determining if any have - * the indicated name, returning either the index or a negative value which can be used to hint - * at where a target with such a name should be inserted into the list - * - * @param $s1 Name to search for - * @param $o2 Reference to List of sorted Targets to be searched - * @return The index of the target with the requested name, or negative value if not found. - * Logic : Using a binary search, compare passed name with those of the targets in the list. - * If target with same name found, return its index; otherwise, return special negative - * value = -(potential_index+1) where potential_index is where a target with the passed - * name would go if it were to be inserted. - */ -func findTarget() { local binsrch_low, binsrch_mid, binsrch_high localobj targetList - targetList = $o2 - - //search through list, using binary search to find if this target exists - binsrch_low = 0 - binsrch_high = targetList.count() - - while ( binsrch_low < binsrch_high ) { - binsrch_mid = int(( binsrch_low+binsrch_high)*0.5 ) - - if( strcmp( targetList.object(binsrch_mid).name, $s1 ) < 0 ) { //guess again, lower - binsrch_low = binsrch_mid+1 - } else { - binsrch_high = binsrch_mid - } - } - - if( binsrch_low $2 ) { //guess again, lower - binsrch_high = binsrch_mid-1 - } else if ( $o1.x[binsrch_mid] < $2 ) { //guess again, higher - binsrch_low = binsrch_mid+1 - } else { - return 1 - } - } - - return 0 -} - -/*! - * - getPoints may be called by TargetManager without a Target object - */ -endtemplate TargetUtility - -//------------------------------------------------------------------------------------------------------------ - -objref targetUtil -targetUtil = new TargetUtility() - - -//------------------------------------------------------------------------------------------------------------ -//------------------------------------------------------------------------------------------------------------ -//------------------------------------------------------------------------------------------------------------ - -/*! - * A Target is a collection of gids and potentially more information that designates points in the circuit - * that may be accessed for things such as stimulus injection or report collection. - */ -begintemplate Target - -//------------------------------------------------------------------------------------------------------------ -// Declare Members Variables -//------------------------------------------------------------------------------------------------------------ - -//! Vector containing indexes into the original parsed list of gids, but only those gids on this local cpu -// This is useful to access extra properties of the original gids -objref localCellsMap - -//! -objref gidMembers, subtargets, flattenedGids, completeFlattenedGids, completeSortedGids -objref cellExtraValues, targetExtraValues, cellSubsets, targetSubsets, nilSecRef - -strdef name, type - -external targetUtil, local_gids - -//------------------------------------------------------------------------------------------------------------ -// Public Accessibility -//------------------------------------------------------------------------------------------------------------ - -public init - -//! general information -public name, type -//! all gids listed directly within a target definition as parsed (i.e no gids from subtargets are included) -public gidMembers -//! list of references to other targets whose contents are to be included for certain operations -public subtargets - -//! For section/cmp targets, there may be a subset descriptor (soma, axon, dend, apic, all) for each member (default=all) -public cellSubsets, targetSubsets -//! For section/cmp/synapse targets, lists of vectors for normalized distances or synapse ids -public cellExtraValues, targetExtraValues - -//! memberfunctions - see below for detailed information -public evalType, isCellTarget, isSectionTarget, isCompartmentTarget, isSynapseTarget -public update, getgid, gids, completegids, getCellExtraValues, getCellSubset, evaluateSynapses, getPointList, getCompletePointList -public isDefined, isEmpty, isUpToDate, completeContains, contains, localCount, getCellCount -public set_offset, get_offset - -//------------------------------------------------------------------------------------------------------------ -// Member functions -//------------------------------------------------------------------------------------------------------------ - -/*! - * Constructor - * - * $s1 Target name - optional argument - * $s2 Target type [Cell, Section, Compartment, Synapse] - optional - */ -proc init() { - localCellsMap = new Vector() //this is not public because it maps to indices of the gidMembers vector which might cause confusion - gidMembers = new Vector() - subtargets = new List() - cellExtraValues = new List() - targetExtraValues = new List() - cellSubsets = new List() - targetSubsets = new List() - - if( argtype(1) == 2 ) { - name = $s1 - } - - // 0=Cell, 1=Section, 2=Compartment, 3=Synapse - typeCode = 0 - if( argtype(2) == 2 ) { - type = $s2 - evalType() - } - - // The offset of gids instantiated. Must be checked to avoid the same target - // to be used in two different population contexts - gid_offset = 0 - - flattenedGids = new Vector() - isFlattened = 0 - - completeFlattenedGids = new Vector() - isCompleteFlattened = 0 - - completeSortedGids = new Vector() - isCompleteSorted = 0 - - isDefined = 0 // the target's contents are defined to prevent multiple targets with same name being *merged* - isEmpty = 0 // target's content is empty - completeMode = 0 // indicates that addressing functions (address*Points) should fill in nils - hasCellsCounted = 0 - isUpToDate = 0 -} - - -/*! - * Sets gid offsetting for this and sub targets - * This enables multiple populations to be instantiated. - * Naturally the same targets cant be used by several circuits - * @arg 1: The gid offset to be applied - */ -proc set_offset() { - if ( $1 == gid_offset) { - return - } - gid_offset = $1 - isUpToDate = 0 - - for targetIndex=0, subtargets.count()-1 { - subtargets.object(targetIndex).set_offset($1) - } -} - -/*! - * Gets gid offset of the target - */ -func get_offset() { - return gid_offset -} - -//------------------------------------------------------------------------------------------------------------ - -/*! - * @return 1 if this target is a Cell type; otherwise, returns 0 - */ -func isCellTarget() { - if( typeCode == 0 ) { - return 1 - } - return 0 -} - -//------------------------------------------------------------------------------------------------------------ - -/*! - * @return 1 if this target is a Section type; otherwise, returns 0 - */ -func isSectionTarget() { - if( typeCode == 1 ) { - return 1 - } - return 0 -} - -//------------------------------------------------------------------------------------------------------------ - -/*! - * @return 1 if this target is a Compartment type; otherwise, returns 0 - */ -func isCompartmentTarget() { - if( typeCode == 2 ) { - return 1 - } - return 0 -} - -//------------------------------------------------------------------------------------------------------------ - -/*! - * @return 1 if this target is a Synapse type; otherwise, returns 0 - */ -func isSynapseTarget() { - if( typeCode == 3 ) { - return 1 - } - return 0 -} - -//------------------------------------------------------------------------------------------------------------ - -/*! - * After putting a string in the type field of this target, have it evaluate - * that string and set the appropriate internal target type code - */ -proc evalType() { - if( strcmp( type, "Cell" ) == 0 ) { - typeCode = 0 - } else if ( strcmp( type, "Section" ) == 0 ) { - typeCode = 1 - } else if ( strcmp( type, "Compartment" ) == 0 ) { - typeCode = 2 - } else if ( strcmp( type, "Synapse" ) == 0 ) { - typeCode = 3 - } else { - typeCode = -1 //invalid target type - } -} - -//------------------------------------------------------------------------------------------------------------ - -/*! - * For each cell member in this target, determine if it exists on the local processor, and put it into localCellsMap (i.e. mapping that gid to an index). - * Logic : Since there should be fewer cells on node than in the target, actually checks each gid on node - * for membership in this target. If the cell is present, the index of the cell in the gidMembers vector - * is kept (to maintain coordination with extra data) - */ -proc ensureUpdated() { local gid, cellIndex, binsrch_low, binsrch_high, binsrch_mid, do_reset localobj allCellIndices - if( isUpToDate ) { - return - } - localCellsMap.resize(0) // clear - isFlattened = 0 - //Get a vector with indexes that tell the sorted order of the gidMembers vector (don't actually sort gidMembers or risk mixing up 'extra' data) - allCellIndices = gidMembers.sortindex() - - //for each cell on node, determine if it exists in this target - for cellIndex=0, local_gids.size()-1 { - gid = local_gids.x[cellIndex] - gid_offset - - //special binary_search - indirect vector access - binsrch_low = 0 - binsrch_high = allCellIndices.size() - - while ( binsrch_low < binsrch_high ) { - binsrch_mid = int(( binsrch_low+binsrch_high)*0.5 ) - - if( gidMembers.x[allCellIndices.x[binsrch_mid]] < gid ) { //guess again, higher - binsrch_low = binsrch_mid+1 - } else { - binsrch_high = binsrch_mid - } - } - - if( binsrch_low SecSize ) { - errorMsg = new String() - sprint( errorMsg.s, "Invalid section access for target %s: requested section %s[%d] (%d) available", name, $s3, $o4.x[count], SecSize ) - terminate( errorMsg.s ) - } - - addressSection( tobj, gid, $s3, $o4.x[count], $o5, $o6 ) - } - } - } -} - -/***************************************************************** - Addresses into a section at the given locations, or all locations if none are told - Required Inputs : $o1 cell object - $2 base gid - $s3 SecName - $4 SecNo - $o5 Vector each element containing Compartment ref. - $o6 TPointList - Logic : -*****************************************************************/ -proc addressSection(){local segcnt, SegNo, locIndex localobj VecObj, tObj - strdef tstr, tstr1 - tObj = $o1 - - if(!section_exists($s3, $4, tObj)) { - if( completeMode == 1 ) { - VecObj = $o5 - if( VecObj.size == 0 ) { - VecObj = new Vector() //the vector was empty, so need to address all points - for segcnt=1, tObj.segCounts.x[$4] { - VecObj.append( segcnt / (tObj.segCounts.x[$4]+1) ) - } - } - - for locIndex=0, VecObj.size()-1 { - $o6.append( nilSecRef, -1 ) - } - } - } else { - sprint(tstr, "%s.%s[%d]", tObj, $s3, $4) - VecObj = $o5 - if(VecObj.size == 0){ //the vector was empty, so need to address all points - VecObj = new Vector() - sprint(tstr1, "access %s", tstr) //can I use the saved section counts like for the off-cpu (complete) case - execute(tstr1) - for segcnt=1, nseg { - VecObj.append(segcnt/(nseg+1)) - } - } - - for locIndex=0, VecObj.size()-1 { - $o6.append( tstr, VecObj.x[locIndex] ) - } - } -} - -/***************************************************************** - This function addresses into the soma of the given cell at the given location(s), or - for all sections (i.e. 1 @ 0.5) if no locations were given - Required Inputs : $o1 cell object - $2 base gid - $o3 Vector with each element containing a Compartment location (typically 0.5). - $o4 TPointList for this gid where section refs will be stored - Logic : Get the cell, determine which locations need addressing, put into CellInfo Obj -*****************************************************************/ -proc addressSoma(){local segcnt, SegNo localobj VecObj, cellObj - strdef tstr, tstr1 - cellObj = $o1 - - //make sure it has a soma - if we are returning a complete TPointList, what do I append? Also, how do I remember I am in - // 'complete' mode? Should I just pass a flag around these functions? - if(!section_exists("soma", cellObj)) { - if( completeMode == 1 ) { - - VecObj = $o3 - - //If no normalized section references were specified, then supply all of them (of course the soma typically only has one) - if(VecObj.size == 0){ - VecObj = new Vector() - for segcnt=1, cellObj.segCounts[0] { //note: the soma should always be in the 0th section - VecObj.append(segcnt/(cellObj.segcnt[0]+1)) - } - } - - //add the address string and the vector of normalized section references - for segcnt=0, VecObj.size-1 { - $o4.append( nilSecRef, -1 ) - } - } - } else { - sprint(tstr, "%s.soma", cellObj ) - VecObj = $o3 - - //If no normalized section references were specified, then supply all of them (of course the soma typically only has one) - if(VecObj.size == 0){ - VecObj = new Vector() - sprint(tstr1, "access %s", tstr) - execute(tstr1) - for segcnt=1, nseg { - VecObj.append(segcnt/(nseg+1)) - } - } - - //add the address string and the vector of normalized section references - for segcnt=0, VecObj.size-1 { - $o4.append( tstr, VecObj.x[segcnt] ) - } - } -} - -//------------------------------------------------------------------------------------------------------------ - -/*! - * Map an actual index for a section from a normalized value - * - * @param $1 Total section count - * @param $2 Normalized section value - * @return The resulting section index - */ -func getSectionNo(){ local SecNo - - if( $2>=0 && $2 <=1){ - SecNo = $2 * ($1-1) - }else { - terminate("Error Normalized Section not within range 0 to 1") - } - - return(int(SecNo)) -} - -endtemplate Target diff --git a/core/hoc/TargetManager.hoc b/core/hoc/TargetManager.hoc deleted file mode 100644 index f10387d6..00000000 --- a/core/hoc/TargetManager.hoc +++ /dev/null @@ -1,720 +0,0 @@ -/** - * @file TargetManager.hoc - * @brief Provides an entity to encapsulate the steps to get cell points for use in Stimulation or Reporting - * @brief This should make certain features of the simulation transparent such as load balancing and target file format - * @author king - * @date 2009-06-12 - * @remark Copyright © BBP/EPFL 2005-2011; All rights reserved. Do not distribute without further notice. - */ - -//make sure we have some basic things loaded -{load_file( "nrngui.hoc" )} -{load_file( "Map.hoc" )} -{load_file( "SerializedSections.hoc" )} -{load_file("fileUtils.hoc")} - -/** - * Target Manager itself, providing the interface for the user to get cell points to use - */ -begintemplate TargetManager - -{load_file("Target.hoc")} - -external terminate - -//----------------------------------------------------------------------------------------------- -// Declare member variables -//----------------------------------------------------------------------------------------------- - -objref targetList, cellList, sectionAccess, isec2sec[1], cellDistributor, deletedSecRef - -///a dummy section which will serve as a nil for populating TPointLists -create dummy - - -//----------------------------------------------------------------------------------------------- -// public members -//----------------------------------------------------------------------------------------------- - -public init -public getTarget, setTargetList, getPointList, locationToPoint, cellDistributor, gids -public getSomaticSections, getApicalSections, getBasalSections, getAxonalSections, registerCell -public getApicalPoints, getBasalPoints, compartmentCast -public selectRandomPoints, getCells, getMETypes - -// Ruben's PFE -public getApicalAtDistance, getSomaPoint, getCell - -//----------------------------------------------------------------------------------------------- -// Member function implementations -//----------------------------------------------------------------------------------------------- - -/*! - * $o1 List with targets most likely parsed from Target files - * $o2 CellDistributor object for accessing sections - */ -proc init() { local lastgid, cellIndex localobj cellObj, secList, pointList - - if( argtype(1) == 1 ) { //need to sort out additional use cases for TargetManager - targetList = $o1 - cellDistributor = $o2 - } - nErrors = 0 - - sectionAccess = new Map() - - //can I create a dummy section, reference it, then delte it to keep a null SectionRef for insertion into pointlists? - access dummy - deletedSecRef = new SectionRef() - delete_section() -} - -//------------------------------------------------------------------------------------------ - -/*! - * Assign a targetList (most likely created by a TargetParser) to this TargetManager - * - * @param $o1 Existing targetList this TargetManger should make a reference to. - */ -proc setTargetList() { - targetList = $o1 -} - -//------------------------------------------------------------------------------------------ - -/*! - * If this TargetManager is not being used for a bluebrain simulaion, the user must register - * any cells so that the TargetManager can find them (under a bluebrain sim, the pnm object - * or something wrapping it will be used) - * - * @param $o1 Reference to a cell object the user has created - */ -proc registerCell() { - if( object_id(cellList) == 0 ) { - cellList = new List() - } - - //should I sort this list? well, yes, but not now :P - cellList.append( $o1 ) -} - -//------------------------------------------------------------------------------------------ - -/*! - * Cell registered with the function registerCell can be obtained by gid - * - * @param $1 gid - */ -obfunc getCell() { - if( object_id(cellList) == 0 ) { - print "getCell: No cell in list." - } - - for listindex=0, cellList.count()-1 { - if( cellList.object(listindex).gid == $1 ) { - return cellList.object(listindex) - } - } - // not found - print "getCell: cell not in list." -} - -//------------------------------------------------------------------------------------------ - -/*! - * Convenience functin for objects like StimuluManager to get access to cell objects without - * having a direct line to the CellDistributor object. Subject to debate, but there must be some way - * for the StimulusManager to access the threshold value for the cells to execute a subThreshols stim (for example). - * - * @param $s1 Target Name to get the gids and collect references to cell objects - * @return List containing cell object references - */ -obfunc getCells() { local cellIndex localobj target, gids, resultList - - resultList = new List() - - target = getTarget( $s1 ) - gids = target.gids() - - for cellIndex=0, gids.size()-1 { - resultList.append( cellDistributor.getCell(gids.x[cellIndex]) ) - } - - return resultList -} - -//------------------------------------------------------------------------------------------ - -/*! - * Convenience functin for objects like StimuluManager to get access to cell objects without - * having a direct line to the CellDistributor object. Subject to debate, but there must be some way - * for the StimulusManager to access the threshold value for the cells to execute a subThreshols stim (for example). - * - * @param $s1 Target Name to get the gids and collect references to cell objects - * @return List containing cell object references - */ -obfunc getMETypes() { local cellIndex localobj target, gids, resultList - - resultList = new List() - - target = getTarget( $s1 ) - gids = target.gids() - - for cellIndex=0, gids.size()-1 { - resultList.append( cellDistributor.getMEType(gids.x[cellIndex]) ) - } - - return resultList -} - -//------------------------------------------------------------------------------------------ - -/*! - * Helper function to aid debugging. Print targets and info about them to stdout - * Required Inputs : - * Logic : For each target, print info about all targets to stdout - */ -proc printTargets() { - //print names of targets found - for targetIndex=0, targetList.count()-1 { - print "Found target ", targetList.object(targetIndex).name - print "\t All Cells: ", targetList.object(targetIndex).gidMembers.size() - print "\t Local Cells: ", targetList.object(targetIndex).gids.size() - print "\t Targets: ", targetList.object(targetIndex).subtargets.count() - for subIndex=0, targetList.object(targetIndex).subtargets.count()-1 { - print "\t\t", targetList.object(targetIndex).subtargets.o[subIndex].name - } - } -} - -//------------------------------------------------------------------------------------------ - -/*! - * Public function which retrieve a target from the list and returns it to the caller - * Required Inputs : $s1 Name of target to find and retrieve - * Logic : Use binary search to locate the target, throwing an error if the target - * is not found - */ -obfunc getTarget() { local targetIndex, potentialGid localobj targetUtil, singleTarget - strdef errorMessage, remake - - targetUtil = new TargetUtility() - targetIndex = targetUtil.findTarget( $s1, targetList ) - if( targetIndex < 0 ) { - // not found -> check if it is just a gid and create a basic target - - //is this entry a cell? : make sure not to confuse some names for gids( i.e. a17sources = target, a17 = gid ) - // test this by trying to extract a gid then putting it back. If either step fails, then can't be a gid - sscanf($s1, "a%d", &potentialGid ) - sprint( remake, "a%d", potentialGid ) - - if( strcmp( $s1, remake ) == 0 ) { //cell detected - singleTarget = new Target( "gid", "Cell" ) - singleTarget.gidMembers.append( potentialGid ) - return singleTarget - } - - // otherwise, at a loss for what to do - sprint( errorMessage, "Error: target %s not found", $s1 ) - terminate( errorMessage ) - } else { - return targetList.object(targetIndex) - } -} - -//------------------------------------------------------------------------------------------ - -/*! - * Retrieve TPointList for the requested item. If a gid is passed, a TPointList is retrieved for the given gid; otherwise, - * a TargetName is expected and a List of TPointLists are retrieved for members of the Target which are local to the cpu - * - * @param $1 or $s1 gid of a cell registered with the TargetManager, or name of Target parsed from a target file - * @return TPointList for single gid, List for Target - */ -obfunc getPointList() { local cellIndex, cmpIndex localobj activeTarget, ptList, activeCell - - //if cells were registered directly, use that cell list to get cell objects for addressing sections; otherwise, - // we will use a pnm object (for now - maybe it will be wrapped in another object?) - if( argtype(1) == 0 && object_id(cellList) ) { - //no target object. Just build pointlist - // Note: this does not allow for split cells - ptList = new TPointList($1) - - //find cell - for cellIndex=0, cellList.count()-1 { - activeCell = cellList.o(cellIndex) - - if( activeCell.gid == $1 ) { - - //to avoid code duplication, should probably create a temporary target containing the cell and use that for getting all compartment points - // does that work without a pnm or CellDistributor object? I don't think so. - - //just mimic all compartment report - //TODO: should put in getCompletePointList instead. Let this function return just soma, I think - forsec activeCell.CellRef.all { - for cmpIndex=1, nseg { - ptList.append( new SectionRef(), cmpIndex/(nseg+1) ) - } - } - - return ptList - } - } - - //if reached this point of the code, the cell was not in the list - print "No such gid registered with TargetManager: ", $1 - return ptList - - } else if( argtype(1) == 2 && object_id(cellDistributor) ) { - - activeTarget = getTarget($s1) - return activeTarget.getPointList( cellDistributor ) - - } else { - print "no cell list? no CellDistributor? where am I supposed to get cell sections from?. I'm leaving" - } -} - -//------------------------------------------------------------------------------------------ - -/*! - * Convenience function to get the gids from a target by name - * - * @param $s1 targetName - * @return Vector containing gids of target members local to this cpu - */ -obfunc gids() { - return getTarget($s1).gids() -} - -//------------------------------------------------------------------------------------------ - -/*! - * Convenience function to encapsulate the selection of random points on a cell. - * - * @param $o1 PointList - * @param $2 nPoints - * @param $3 seed - * @return resultant TPointList with selected points in no particular order - */ -obfunc selectRandomPoints() { local synIndex, cmpIndex, x localobj ptList, allPoints, rng - - //first, get all points - allPoints = $o1 - - print "got ", allPoints.count(), " points" - - rng = new Random( $3 ) - rng.MCellRan4( $3 ) - rng.discunif( 0, allPoints.count()-1 ) - - ptList = new TPointList() - - for synIndex=0, $2-1 { - cmpIndex = rng.repick() - - //Do I need to check for existence? - if( !allPoints.sclst.o(cmpIndex).exists() ) { //section not on this cpu - continue - } - - //allPoints.access_(cmpIndex) - for allPoints.point(cmpIndex, &x) { - ptList.append( new SectionRef(), x ) - } - - } - - return ptList -} - -//------------------------------------------------------------------------------------------ - -/*! - * Convenience function to get points for each apical section from a Cell target - * - * @param $1/$s1 gid or TargetName - * @return TPointList/List where each TPointList contains a point to each apical section of the given target - */ -obfunc getApicalPoints() { local cellIndex localobj tempTarget, ptList, activeCell - - if( argtype(1) != 0 ) { - terminate("getApicalPoints() now only accepts one int arg") - } - - if( object_id(cellList, 1) == -1 ) { - terminate( "No cells registered with my null cellList" ) - } - - ptList = new TPointList($1) - - //find cell - for cellIndex=0, cellList.count()-1 { - activeCell = cellList.o(cellIndex) - - if( activeCell.gid == $1 ) { - - //just mimic all compartment report - //TODO: can we make this code more general and put in a single, utility function. avoid copy/paste in other get(*subset*)Points funcs? - forsec activeCell.CellRef.apical { - for cmpIndex=1, nseg { - ptList.append( new SectionRef(), cmpIndex/(nseg+1) ) - } - } - - return ptList - } - } - - //if reached this point of the code, the cell was not in the list - print "No such gid registered with TargetManager: ", $1 - return ptList -} - -//------------------------------------------------------------------------------------------ - -/*! - * Convenience function to get points for each basal section from a Cell target - * - * @param $s1 TargetName - * @return List where each TPointList contains a point to each basal section of the given target - */ -obfunc getBasalPoints() { local cellIndex localobj tempTarget, ptList, activeCell - - if( argtype(1) != 0 ) { - terminate("getApicalPoints() now only accepts one int arg") - } - if( object_id(cellList, 1) == -1 ) { - terminate( "No cells registered with my null cellList" ) - } - - ptList = new TPointList($1) - - //find cell - for cellIndex=0, cellList.count()-1 { - activeCell = cellList.o(cellIndex) - - if( activeCell.gid == $1 ) { - - //just mimic all compartment report - //TODO: can we make this code more general and put in a single, utility function. avoid copy/paste in other get(*subset*)Points funcs? - forsec activeCell.CellRef.basal { - for cmpIndex=1, nseg { - ptList.append( new SectionRef(), cmpIndex/(nseg+1) ) - } - } - - return ptList - } - } - - //if reached this point of the code, the cell was not in the list - print "No such gid registered with TargetManager: ", $1 - return ptList -} - -//------------------------------------------------------------------------------------------ -// The following functions are supposedly temporary. They are here to allow us to compare -// npoisson from old bglib to new bglib. However, they should be supplanted by the get(*subset*)Points -// functions later -//------------------------------------------------------------------------------------------ - -/*! - * Convenience function to get points for each somatic section from a Cell target - * - * @param $s1 TargetName - * @return List where each TPointList contains a point to each somatic section of the given target - */ -obfunc getSomaticSections() { localobj tempTarget - tempTarget = sectionCast( $s1, "soma" ) - return tempTarget.getCompletePointList( cellDistributor ) -} - -//------------------------------------------------------------------------------------------ - -/*! - * Convenience function to get points for each apical section from a Cell target - * - * @param $s1 TargetName - * @return List where each TPointList contains a point to each apical section of the given target - */ -obfunc getApicalSections() { localobj tempTarget - - tempTarget = sectionCast( $s1, "apic" ) - return tempTarget.getCompletePointList( cellDistributor ) -} - -//------------------------------------------------------------------------------------------ - -/*! - * Convenience function to get points for each basal section from a Cell target - * - * @param $s1 TargetName - * @return List where each TPointList contains a point to each basal section of the given target - */ -obfunc getBasalSections() { localobj tempTarget - tempTarget = sectionCast( $s1, "dend" ) - return tempTarget.getCompletePointList( cellDistributor ) -} - -//------------------------------------------------------------------------------------------ - -/*! - * Convenience function to get points for each axonal section from a Cell target - * - * @param $s1 TargetName - * @return List where each TPointList contains a point to each axonal section of the given target - */ -obfunc getAxonalSections() { localobj tempTarget - tempTarget = sectionCast( $s1, "axon" ) - return tempTarget.getCompletePointList( cellDistributor ) -} - -//------------------------------------------------------------------------------------------ - -/*! - * Utility function to take an existing Cell target and wrap it in a Section target which - * will return 1 compartment in thecenter of each section - * - * @param $o1/$s1 Target object or Target name - * @param $s2 subset name [soma, apic, dend, apic, or "" for all] - */ -obfunc sectionCast() { localobj target, wrapper - - if( argtype(1) == 2 ) { - target = getTarget( $s1 ) - } else { - target = $s1 - } - - if( !target.isCellTarget() ) { - terminate( "Attempting to apply sectionCast to non-Cell target" ) - } - - wrapper = new Target( "temp", "Section" ) - wrapper.subtargets.append( target ) - wrapper.targetSubsets.append( new String($s2) ) - wrapper.targetExtraValues.append( new Vector() ) - - return wrapper -} - -//------------------------------------------------------------------------------------------ - -/*! - * Utility function to take an existing Cell target and wrap it in a Compartment target - * which will return a point for each compartment in each section - * - * @param $o1/$s1 Target object or Target name - * @param $s2 subset name [soma, apic, dend, apic, or "" for all] - */ -obfunc compartmentCast() { localobj target, wrapper - - if( argtype(1) == 2 ) { - target = getTarget( $s1 ) - } else { - target = $o1 - } - - if( !target.isCellTarget() ) { - terminate( "Attempting to apply compartmentCast to non-Cell target" ) - } - - wrapper = new Target( "temp", "Compartment" ) - wrapper.subtargets.append( target ) - wrapper.targetSubsets.append( new String($s2) ) - wrapper.targetExtraValues.append( new Vector() ) - - return wrapper -} - -//------------------------------------------------------------------------------------------ - -/*! - * Given a location for a cell, section id, segment id, and offset into the segment, create a TPointList - * containing a section reference to there - * - * @param $1 gid - * @param $2 section index (isec) - * @param $3 distance to start of segment (ipt) - * @param $4 offset distance beyond the ipt (microns) - * @return TPointList with 1 item, where the synapse should go - */ -obfunc locationToPoint() { local gid, distance, isec, ipt, offset, total localobj resultPoint, cellSections, tmpSection, nil, helpCell - strdef errorMessage - - gid = $1 - isec = $2 - ipt = $3 - offset = $4 - - if( offset < 0 ) { //soma connection, just zero it - offset = 0 - } - - resultPoint = new TPointList( gid ) - - cellSections = gidToSections( gid ) - if( object_id(cellSections) == 0 ) { - terminate( "getting locations for non-bg sims is not implemented yet..." ) - } - - if ( isec >= cellSections.count() ) { - sprint( errorMessage, "Error: section %d out of bounds (%d total). Morphology section count is low, is this a good morphology?", isec, cellSections.count() ) - terminate( errorMessage ) - } - - distance = 0.5 - tmpSection = cellSections.isec2sec[isec] - - if (tmpSection == nil) { //assume we are in LoadBalance mode (until we decide to verify this either here or in a higher layer) - resultPoint.append( deletedSecRef, -1 ) - } else if (ipt == -1) { - // SYN2/Sonata spec have the a pre-calculated distance field. - // In such cases segment (ipt) is -1 and offset is that distance. - if( offset == 0.0 ) { - offset = 0.0000001 - } - if( offset >= 1.0 ) { - offset = 0.9999999 - } - resultPoint.append( tmpSection, offset ) - } else { - tmpSection.sec { - // when a cell is split, the path from the split point - // to the root gets its 3d points reversed because - // the section orientation is reversed. - if (section_orientation() == 1) { - ipt = n3d() - 1 - ipt - offset = -offset - } - - if (ipt < n3d()) { - //distance = (arc3d(ipt))/L - distance = (arc3d(ipt)+offset)/L - - if( distance == 0.0 ) { - distance = 0.0000001 - } - - if( distance >= 1.0 ) { - total = (arc3d(ipt)+offset) - distance = 0.9999999 - } - } - if (section_orientation() == 1) { - distance = 1 - distance - } - } - - resultPoint.append( tmpSection, distance ) - } - - return resultPoint -} - -//------------------------------------------------------------------------------------------ - -/*! - * For a given gid, we want SectionRefs stored in a List for random access. This function will - * see if such a List already exists and return it, or build that list as needed, storing it for the future before - * returning the result. Note that in split cases, the deleted sections will still have their SectionRefs. - * Should this function instead go in the Node object? The CellDistributor object? - * - * @param $1 gid - * @return List object with SectionRefs to every section in the cell - */ -obfunc gidToSections() { local n localobj resultSerial, nil, cellRef - - resultSerial = sectionAccess.get( $1 ) - - if( object_id(resultSerial) == 0 ) { //must build new one - //determine if we should use pnm or cellList based on which is actually allocated - if( object_id(cellDistributor) ) { - cellRef = cellDistributor.getCell($1) - } else if( object_id(cellList) ) { - print "cellList implementation pending" - return nil - } else { - print "no cell list? no CellDistributor? where am I supposed to get cell sections from?. I'm leaving" - return nil - } - - resultSerial = new SerializedSections( cellRef ) - - sectionAccess.put( $1, resultSerial ) - } - - return resultSerial -} - -//------------------------------------------------------------------------------------------ -// -// for Ruben's PFE -// - -/*! - * Given cell object find all Points within the given distance to the soma. - * Return the one where the diameter is maximal. - * - * @param $1 cell gid - * @param $2 distance to soma in micrometers - * @return PointList object - */ -obfunc getApicalAtDistance() { local radius, max_diam, location, i, secnum, d0, d1, stimdist, dd localobj cell, ptList, sl - cell = getCell($1) - cell.CellRef.soma distance() - radius = $2 - i = 0 - max_diam = 0. - forsec cell.CellRef.apical { - d0 = distance(0) - d1 = distance(1) - if(d0 <= radius && d1 >= radius) { - if(diam > max_diam) { - location = (radius - distance(0)) / (distance(1) - distance(0)) - secnum = i - max_diam = diam - } - } - i = i + 1 - } - - ptList = new TPointList() - if(max_diam > 0) { - access cell.CellRef.apic[secnum] - ptList.append( location ) - // debug output - //printf("Found apical with diam %f.\n", max_diam) - } - - // compare to Etay's code: - /* - stimdist = 620/cell.getLongestBranch(cell.CellRef.apical) - printf("stimdist: %f\n", stimdist) - sl = cell.locateSites(cell.CellRef.apical, stimdist) - for(i1=0;i1 1) printCounts = ($2 && isVerbose) - - newTargets = new List() - nErrors = 0 - - if( isVerbose ) { - print "Reading target file ", $s1 - } - - if (tFile.isopen == 0) { - terminate( "Error : Could not open target file", $s1 ) - } - - //Read till end of file - while(tFile.gets(tstr) >= 0){ - sscanf(tstr, "%s",tstr1) - if( 0 == strcmp (tstr1, "Target")) { //have reached a target block - targetInstance = new Target() - sscanf(tstr, "%*s%s%s", typeString, targetInstance.name ) //extract target type and name - - //if this target was seen before (e.g. used as subtarget), a reference already exists in the list and needs to be populated; - // otherwise the new reference we just made will be inserted into the list and given back directly - targetInstance = updateTargetList( targetInstance ) - newTargets.append(targetInstance) - - //make sure that contents for this target haven't been seen before - if( targetInstance.isDefined ) { - if ( isVerbose ) { - print "Error. Multiple definitions for Target ", targetInstance.name - } - nErrors = nErrors + 1 - } - - targetInstance.type = typeString - - targetInstance.evalType() - - //target contents - tFile.gets(tstr) //Grab '{' - sscanf(tstr, "%s",tstr1) - if( targetInstance.isCellTarget() ) { //only need to parse into cellnames and target names - targetInstance.isEmpty = parseCellTargetContents( tFile, targetInstance ) - } else { - targetInstance.isEmpty = parseOtherTargetContents( tFile, targetInstance ) - } - - //note that the closing '}' was grabbed during whichever parse*Contents function call - targetInstance.isDefined = 1 - - //Display msg to default output if target instance is empty - if (targetInstance.isEmpty && isVerbose) { - print "Warning: Following target is empty: ", targetInstance.name - } - } - } - - tFile.close() - if (nErrors > 0) { - terminate ( "Target with multiple definitions" ) - } - if (printCounts) for i=0, newTargets.count()-1 { - subtargets = newTargets.o(i).subtargets.count() - if( subtargets ) { - printf(" * Target %-15.15s %6d cells [%d subtargets]\n", newTargets.o(i).name, newTargets.o(i).getCellCount(), subtargets) - } else { - printf(" * Target %-15.15s %6d cells\n", newTargets.o(i).name, newTargets.o(i).getCellCount()) - } - } -} - -//------------------------------------------------------------------------------------------------------------ - -/*! - * Private function used during parsing to determine if targets already exist with a given name. - * If they do, return a reference to that target; otherwise, add the given reference into the list - * and return the same reference. Perhaps a reference isn't needed, but just the name. Then this - * function can create an instance if necessary. It's not broke now, though, so fix it later as - * needed - * - * @param $o1 Reference to a target whose name needs to be searched for in the list - * @return Reference to an pre-existing with same name, or the parameter reference if it was insterted - * Logic : Using the binary search, determine if a target exists with a particular name, - * if not, insert the given reference into the list at the location as derived by - * the negative return value from the binary search. Return the true reference - */ -obfunc updateTargetList() { local targetIndex localobj targetUtil - targetUtil = new TargetUtility() - - targetIndex = targetUtil.findTarget( $o1.name, targetList ) - if( targetIndex < 0 ) { //not found, so insert it at indicated location - targetIndex = -(targetIndex+1) - if( targetIndex == targetList.count() ) { - targetList.append( $o1 ) - } else { - targetList.insrt( targetIndex, $o1 ) - } - return $o1 - } else { - return targetList.object(targetIndex) - } -} - -//------------------------------------------------------------------------------------------ - -/*! - * For each target in the targetList, update them so that there exist abbreviated lists - * related to cells local to the cpu. - * - * @param $o1 Vector containing gids that are on the local cpu - * @param $2 Append mode. 0 (default) will reset the vector first. - */ -proc updateTargets() { local targetIndex, append_mode localobj activeTarget, tmpvec - append_mode = 0 - if (numarg() > 1) { - append_mode = $2 - } - if (append_mode) { - tmpvec = new Vector($o1) // o1 may be a python obj - local_gids.append(tmpvec) - } else { - local_gids = $o1.c - } - local_gids.sort() - - for targetIndex=0, targetList.count()-1 { - targetList.o(targetIndex).isUpToDate = 0 - } -} - -//------------------------------------------------------------------------------------------------------------ - -/*! - * Cell Targets have a simpler structure where a line can have any number of entries, but limited to gids or other target names - * - * @param $o1 Reference to file where data is coming from - * @param $o2 Reference to Target object where data is stored - * @return Returns error code (1 if target is empty, 0 otherwise) - */ -func parseCellTargetContents() { local potentialGid, numberToken localobj tFile, targetInstance, addition - tFile = $o1 - targetInstance = $o2 - strdef tstr1 - - //the line can have any number of entries - cells or other targets - tFile.scanstr(tstr1) //goto first entry - - while(0 != strcmp(tstr1, "}") && 0 != strcmp( tstr1, "{" ) && !tFile.eof() ) { - potentialGid = addMember(tstr1, targetInstance) - - //count the number of parsed token - numberToken = numberToken + 1 - - //get next string - tFile.scanstr(tstr1) - } - - - validateTargetTermination(targetInstance.name, tstr1) - - //Return error code: 1 if target is empty - if (numberToken == 0) { - return 1 - } else { - return 0 - } -} - -//------------------------------------------------------------------------------------------------------------ - -/*! - * For non-Cell targets, any member (cell or subtarget) must be on separate line. The rest of the line can potentially - * contain section names, section ids (int, but treated as float), normalized distances (float), or synapse ids (int) - * - * @param $o1 Reference to file where data is coming from - * @param $o2 Reference to Target object where data is stored - * @return Returns error code (1 if target is empty, 0 otherwise) - */ -func parseOtherTargetContents() { local offset, extraIndex, potentialGid, val, numberToken localobj tFile, targetInstance, addition, extraData, strobj, extraStrings, extraValues, pendingSubset - tFile = $o1 - targetInstance = $o2 - strdef tstr, tstr1, remake - strobj = new StringFunctions() - - //each line contains a name, then potentially a subset and potentially a list of values - getTrimmedLine( tFile, tstr ) - - sscanf(tstr, "%s",tstr1) - while(0 != strcmp(tstr1, "}") && 0 != strcmp( tstr1, "{" ) && !tFile.eof() ) { - - //count the number of parsed token - numberToken = numberToken + 1 - - potentialGid = addMember(tstr1, targetInstance ) - - //look at remainder of line for [cell|target]Extra field - extraData = new String() - strobj.right( tstr, strobj.len(tstr1) ) - extraData.s = tstr - - extraStrings = new List() - extraValues = new Vector() - pendingSubset = new String("") - - if( targetInstance.isSynapseTarget() ) { - if( split( extraData.s, extraStrings ) > 0 ) { - for extraIndex=0, extraStrings.count()-1 { - sscanf(extraStrings.o(extraIndex).s, "%d", &val) - extraValues.append(val) - } - } else { - //print "no extra data" - } - } else { //section|compartment - if( split( extraData.s, extraStrings ) > 0 ) { - //test if first element is a string : soma, dend, apic, axon - offset = 0 - if( (0==strcmp(extraStrings.o(0).s, "soma")) || (0==strcmp(extraStrings.o(0).s, "axon")) || \ - (0==strcmp(extraStrings.o(0).s, "dend")) || (0==strcmp(extraStrings.o(0).s, "apic"))) { - //targetInstance.cellSections.append( extraStrings.o(0) ) - pendingSubset = extraStrings.o(0) - offset = 1 - } - - for extraIndex=offset, extraStrings.count()-1 { - sscanf(extraStrings.o(extraIndex).s, "%f", &val) - extraValues.append(val) - } - } else { - //print "no extra data" - } - } - - //look at potentialGid again to know where to store this extra data - // note: always append - even empty vector - so that indices correspond with gid/subtarget vecs - if( potentialGid != -1 ) { - targetInstance.cellSubsets.append( pendingSubset ) - targetInstance.cellExtraValues.append( extraValues ) - } else { - targetInstance.targetSubsets.append( pendingSubset ) - targetInstance.targetExtraValues.append( extraValues ) - } - - //goto next line; get first string element - getTrimmedLine( tFile, tstr ) - sscanf(tstr, "%s",tstr1) - } - - validateTargetTermination(targetInstance.name, tstr1) - - //Return error code: 1 if target is empty, 0 if - if (numberToken == 0) { - return 1 - } else { - return 0 - } -} - -//------------------------------------------------------------------------------------------------------------ - -/*! - * Consolidate a validation action that must be performed by both parseCellTargetContents and parseOtherTargetContents. - * Make certain that a closing brace was found as the last string; anything else indicates a parse error. - * Should such a parse error be fatal and compel the program to terminate? - * - * @param $s1 The name of the target parsed - * @param $s2 The last string read while parsing the target - */ -proc validateTargetTermination() { - //we should have found a closing brace at the end (encountering an eof or opening brace instead indicates a parse error) - if( 0 != strcmp($s2, "}") ) { //failed to find closing brace - print "Warning: no closing brace found for Target ", $s1 - //nErrors++? - } -} - -//------------------------------------------------------------------------------------------------------------ - -/*! - * Add a member, either cell or subtarget, to a Target object. The string with the new member's name is checked to - * determine if it a gid (format = "a%d") or not. Note that the expression will be validated to make sure that - * certain target names (e.g. a17sources) are not mistaken for gids - * - * @param $s1 The string getting checked - * @param $o2 Reference to Target object to receive the member - * @return The numeric gid if it validated as a gid (i.e. a%d returns %d); otherwise, returns -1 - */ -func addMember() { local potentialGid localobj targetInstance, addition - strdef remake - targetInstance = $o2 - - //is this entry a cell? : make sure not to confuse some names for gids( i.e. a17sources = target, a17 = gid ) - // test this by trying to extract a gid then putting it back. If either step fails, then can't be a gid - sscanf($s1, "a%d", &potentialGid ) - sprint( remake, "a%d", potentialGid ) - - if( strcmp( $s1, remake ) == 0 ) { //cell detected - targetInstance.gidMembers.append( potentialGid ) - return potentialGid - } else { //must be new or existing target - addition = new Target() - addition.name = $s1 - - //following line either replaces addition with reference to existing target, or addition is used to insert new reference into master list - addition = updateTargetList( addition ) - targetInstance.subtargets.append( addition ) - return -1 - } -} - -//------------------------------------------------------------------------------------------------------------ - -/*! - * Public function which retrieve a target from the list and returns it to the caller - * Required Inputs : $s1 Name of target to find and retrieve - * Logic : Use binary search to locate the target, throwing an error if the target - * is not found - */ -obfunc getTarget() { local targetIndex localobj targetUtil - targetUtil = new TargetUtility() - targetIndex = targetUtil.findTarget( $s1, targetList ) - if( targetIndex < 0 ) { //not found, error - strdef errmsg - sprint( errmsg, "Error: target %s not found. Attempting to terminate", $s1) - terminate( errmsg ) - } else { - return targetList.object(targetIndex) - } -} - -//------------------------------------------------------------------------------------------ - -/*! - * Utility function for handling the retrieval of a line from a file and removing leading - * white space or trailing comments (trailing white space is not checked, but should not be a problem) - * - * @param $o1 Reference to file whose next line is required - * @param $s2 String where the read line should be stored - * Logic : Read the next line, check if the first character is ws, removing all leading ws - * if it is, then check for the comment character '#' and remove it and any other - * characters after it. - */ -proc getTrimmedLine() { localobj targetFile, strobj - strdef tstr, wscheck - strobj = new StringFunctions() - - targetFile = $o1 - targetFile.gets(tstr) - - //check for leading whitespace - wscheck = tstr - strobj.left( wscheck, 1 ) - if( strcmp( wscheck, " " ) == 0 || strcmp( wscheck, "\t" ) == 0 ) { - ret = strobj.tail( tstr, "[ \t]+", tstr ) //strip leading white space - } - - //strip comments - if( strobj.substr(tstr, "#" ) != -1 ) { - ret = strobj.head(tstr, "#", tstr ) - } - - //trailing whitespace shouldn't be a problem, but 'trimmedLine' implies none at front or rear - // this doesn't hurt parsing, though, so I won't deal with it - - $s2 = tstr -} - -//------------------------------------------------------------------------------------------------------------ - -/*! - * Split a string across whitespace into a list of many strings, returning the number of strings found - * maybe this will be moved into global shared hoc file as config file parsing uses a similar function? - * - * @param $s1 String to split - * @param $o1 List which will be populated by string objects - * @return The number of items collected into the final list - */ -func split(){ local RetVal, length localobj strobj, str, Token - strobj = new StringFunctions() - str = new String() - Token = new String() - - str.s = $s1 - i = 0 - while(1){ - RetVal = 0 - RetVal = sscanf(str.s, "%s", Token.s) - if(RetVal <1) { - break - } - i = i+1 - $o2.append(new String(Token.s)) - - RetVal = 0 - RetVal = strobj.substr(str.s,Token.s) - if(RetVal <0){ - break - } - length = 0 - length = strobj.len(Token.s) - strobj.right(str.s, length+RetVal ) - } - return i -} - -endtemplate TargetParser diff --git a/core/hoc/neurodamus.hoc b/core/hoc/neurodamus.hoc index 7ecdc008..f54e02c7 100644 --- a/core/hoc/neurodamus.hoc +++ b/core/hoc/neurodamus.hoc @@ -9,8 +9,6 @@ // Neurodamus - Do not change order {load_file("ConfigParser.hoc")} -{load_file("TargetParser.hoc")} -{load_file("TargetManager.hoc")} {load_file("RNGSettings.hoc")} {load_file("ElectrodeManager.hoc")} {load_file("StimulusManager.hoc")} diff --git a/examples/test_neurodamus.py b/examples/test_neurodamus.py index a8060b0c..bfdfec07 100644 --- a/examples/test_neurodamus.py +++ b/examples/test_neurodamus.py @@ -5,12 +5,13 @@ from __future__ import print_function from neurodamus import Node, Neurodamus from neurodamus.core import NeurodamusCore as Nd -from neurodamus.utils import setup_logging +from neurodamus.utils.logging import setup_logging import sys import logging from os import path as Path -RECIPE_FILE = Path.expanduser("~/dev/TestData/build/circuitBuilding_1000neurons/BlueConfig") +RECIPE_FILE = Path.join(Path.dirname(__file__), + "../tests/simulations/usecase3/simulation_sonata.json") DEFAULT_LOG_LEVEL = 2 TRACE_LOG_LEVEL = 5 @@ -38,7 +39,6 @@ def test_node_run(trace=False): logging.info("Create connections") node.create_synapses() - node.create_gap_junctions() logging.info("Enable Stimulus") node.enable_stimulus() @@ -55,7 +55,7 @@ def test_node_run(trace=False): node.enable_reports() logging.info("Run") - node.run(True) + node.run_all() logging.info("Simulation finished. Gather spikes then clean up.") node.spike2file("out.dat") diff --git a/neurodamus/cell_distributor.py b/neurodamus/cell_distributor.py index 8ce49bee..d586e229 100644 --- a/neurodamus/cell_distributor.py +++ b/neurodamus/cell_distributor.py @@ -331,12 +331,6 @@ def store_metype_stats(metype, n_cells): def _update_targets_local_gids(self): logging.info(" > Updating targets") - cell_offset = self._local_nodes.offset - if cell_offset and self._target_spec.name: - target = self._target_manager.get_target(self._target_spec) - if not hasattr(target, "set_offset"): - raise NotImplementedError("No gid offsetting supported by neurodamus Target.hoc") - target.set_offset(cell_offset) # Add local gids to matching targets self._target_manager.register_local_nodes(self._local_nodes) diff --git a/neurodamus/connection.py b/neurodamus/connection.py index 41a963e6..3f9f6308 100644 --- a/neurodamus/connection.py +++ b/neurodamus/connection.py @@ -265,9 +265,8 @@ def add_synapses(self, target_manager, synapses_params, base_id=0): n_synapses = len(synapses_params) synapse_ids = numpy.arange(base_id, base_id+n_synapses, dtype="uint64") mask = numpy.full(n_synapses, True) # We may need to skip invalid synapses (e.g. on Axon) - for i, syn_params in enumerate(synapses_params): - syn_point = target_manager.hoc.locationToPoint( + syn_point = target_manager.location_to_point( self.tgid, syn_params['isec'], syn_params['ipt'], syn_params['offset']) syn_params['location'] = syn_point.x[0] section = syn_point.sclst[0] diff --git a/neurodamus/connection_manager.py b/neurodamus/connection_manager.py index d917f74c..6a83a3af 100644 --- a/neurodamus/connection_manager.py +++ b/neurodamus/connection_manager.py @@ -603,13 +603,9 @@ def _add_synapses(self, cur_conn, syns_params, syn_type_restrict=None, base_id=0 cur_conn.add_synapses(self._target_manager, syns_params, base_id) # - - def get_updated_population_offsets(self, src_target, dst_target): + def get_population_offsets(self): sgid_offset = self._src_cell_manager.local_nodes.offset tgid_offset = self._cell_manager.local_nodes.offset - if src_target: - src_target.set_offset(sgid_offset) - if dst_target: - dst_target.set_offset(tgid_offset) return sgid_offset, tgid_offset # - @@ -656,7 +652,7 @@ def target_gids(gids): gids = target_gids(gids) created_conns_0 = self._cur_population.count() - sgid_offset, tgid_offset = self.get_updated_population_offsets(src_target, dst_target) + sgid_offset, tgid_offset = self.get_population_offsets() self._synapse_reader.configure_override(mod_override) self._synapse_reader.preload_data(gids) @@ -823,7 +819,7 @@ def get_target_connections(self, src_target_name, if src_target and src_target.is_void() or dst_target.is_void(): return - _, tgid_offset = self.get_updated_population_offsets(src_target, dst_target) + _, tgid_offset = self.get_population_offsets() populations: List[ConnectionSet] = (conn_population,) if conn_population is not None \ else self._populations.values() diff --git a/neurodamus/modification_manager.py b/neurodamus/modification_manager.py index 54e4152a..0e4f6008 100644 --- a/neurodamus/modification_manager.py +++ b/neurodamus/modification_manager.py @@ -46,7 +46,7 @@ def interpret(self, target_spec, mod_info): if not mod_t: raise ConfigurationError("Unknown Modification {}".format(mod_t_name)) target = self._target_manager.get_target(target_spec) - cell_manager = self._target_manager.hoc.cellDistributor + cell_manager = self._target_manager._cell_manager mod = mod_t(target, mod_info, cell_manager) self._modifications.append(mod) @@ -65,10 +65,7 @@ class TTX: Uses TTXDynamicsSwitch as in BGLibPy. Overrides HOC version, which is outdated """ def __init__(self, target, mod_info: dict, cell_manager): - if target.isCellTarget(): # get all compartments for each cell - tpoints = self.compartment_cast(target, "").getPointList(cell_manager) - else: # use compartments present in target - tpoints = target.getPointList(cell_manager) + tpoints = target.getPointList(cell_manager, sections='all') # insert and activate TTX mechanism in all sections of each cell in target for tpoint_list in tpoints: @@ -101,11 +98,7 @@ class ConfigureAllSections: """ def __init__(self, target, mod_info: dict, cell_manager): config, config_attrs = self.parse_section_config(mod_info['SectionConfigure']) - - if target.isCellTarget(): # get all compartments for each cell - tpoints = self.compartment_cast(target, "").getPointList(cell_manager) - else: # use compartments present in target - tpoints = target.getPointList(cell_manager) + tpoints = target.getPointList(cell_manager, sections='all') napply = 0 # number of sections where config applies # change mechanism variable in all sections that have it diff --git a/neurodamus/node.py b/neurodamus/node.py index 3ae49a41..cd5d0019 100644 --- a/neurodamus/node.py +++ b/neurodamus/node.py @@ -310,7 +310,7 @@ def __init__(self, config_file, options=None): # Register the global target and cell manager self._target_manager.register_target(self._circuits.global_target) - self._target_manager.init_hoc_manager(self._circuits.global_manager) + self._target_manager.register_cell_manager(self._circuits.global_manager) # # public 'read-only' properties - object modification on user responsibility @@ -340,8 +340,6 @@ def load_targets(self): log_verbose("Loading targets for circuit %s", circuit._name or "(default)") self._target_manager.load_targets(circuit) - self._target_manager.load_user_target() - # - @mpi_no_errors @timeit(name="Compute LB") @@ -958,13 +956,6 @@ def _report_build_params(self, rep_name, rep_conf, target, pop_offsets_alias_pop or self._default_population) log_verbose("Report on Population: %s, Target: %s", population_name, rep_target.name) - target_population = rep_target.population or self._target_spec.population - pop_offsets, alias_pop = pop_offsets_alias_pop - offset = pop_offsets[target_population] if target_population in pop_offsets \ - else pop_offsets[alias_pop[target_population]] - if offset > 0: # dont reset if not needed (recent hoc API) - target.set_offset(offset) - report_on = rep_conf["ReportOn"] return self.ReportParams( os.path.basename(rep_conf.get("FileName", rep_name)), @@ -984,45 +975,6 @@ def _report_build_params(self, rep_name, rep_conf, target, pop_offsets_alias_pop # def _report_write_coreneuron_config(self, rep_name, rep_conf, target, rep_params): target_spec = TargetSpec(rep_conf["Target"]) - corenrn_target = target - - # 0=Compartment, 1=Cell, Section { 2=Soma, 3=Axon, 4=Dendrite, 5=Apical } - target_type = target.isCellTarget() - compartment_offset = 0 - - # Note: CoreNEURON does not support more than one Section in a Compartment target - if target.isCompartmentTarget(): - for activeTarget in target.subtargets: - if activeTarget.isSectionTarget(): - # Ensure only one Section (i.e., number of subtargets is also 1) - if target.subtargets.count() > 1: - logging.error("Report '%s': only a single Section is allowed in a " - "Compartment target", rep_name) - return False - - # If we reach this point, update original target and offset - corenrn_target = activeTarget - compartment_offset = 4 - - if corenrn_target.isSectionTarget(): - if (corenrn_target.subtargets.count() != 1 - or not corenrn_target.subtargets[0].isCellTarget()): - logging.error("Report '%s': only a single Cell subtarget is allowed in a " - "Section target", rep_name) - return False - - section_type = corenrn_target.targetSubsets[0].s - - if section_type == "soma": - target_type = 2 + compartment_offset - elif section_type == "axon": - target_type = 3 + compartment_offset - elif section_type == "dend": - target_type = 4 + compartment_offset - elif section_type == "apic": - target_type = 5 + compartment_offset - else: - target_type = 0 # for sonata config, compute target_type from user inputs if "Sections" in rep_conf and "Compartments" in rep_conf: @@ -1035,6 +987,7 @@ def _compute_corenrn_target_type(section_type, compartment_type): raise ConfigurationError(f"Report: invalid compartment type {compartment_type}") if section_type == "all": # for "all sections", support only target_type=0 return 0 + # 0=Compartment, Section { 2=Soma, 3=Axon, 4=Dendrite, 5=Apical, 6=SomaAll ... } return sections.index(section_type)+1+4*compartments.index(compartment_type) section_type = rep_conf.get("Sections") @@ -1065,15 +1018,14 @@ def _report_setup(self, report, rep_conf, target, rep_type): # value on the soma. Otherwise, get target points as normal. sections = rep_conf.get("Sections") compartments = rep_conf.get("Compartments") - is_cell_target = target.isCellTarget(sections=sections, - compartments=compartments) + sum_currents_into_soma = (sections == "soma" and compartments == "center") # In case of summation in the soma, we need all points anyway - if is_cell_target and rep_type == "Summation": + if sum_currents_into_soma and rep_type == "Summation": sections = "all" compartments = "all" - points = self._target_manager.get_target_points(target, global_manager, - sections=sections, - compartments=compartments) + points = self._target_manager.getPointList(target, + sections=sections, + compartments=compartments) for point in points: gid = point.gid pop_name, pop_offset = global_manager.getPopulationInfo(gid) @@ -1086,7 +1038,7 @@ def _report_setup(self, report, rep_conf, target, rep_type): cell, point, spgid, SimConfig.use_coreneuron, pop_name, pop_offset) elif rep_type == "Summation": report.addSummationReport( - cell, point, is_cell_target, spgid, SimConfig.use_coreneuron, + cell, point, sum_currents_into_soma, spgid, SimConfig.use_coreneuron, pop_name, pop_offset) elif rep_type == "Synapse": report.addSynapseReport( @@ -1139,8 +1091,7 @@ def execute_neuron_configures(self): log_verbose("Apply configuration \"%s\" on target %s", config.get("Configure").s, target_name) - points = self._target_manager.get_target_points(target_name, - self._circuits.base_cell_manager) + points = self._target_manager.getPointList(target_name) # iterate the pointlist and execute the command on the section for tpoint_list in points: for sec_i, sc in enumerate(tpoint_list.sclst): @@ -1349,7 +1300,6 @@ def _sim_corenrn_configure_datadir(self, corenrn_restore): "Increase the number of nodes to reduce the memory footprint " "(Current use node: %d MB / SHM Limit: %d MB / Mem. Limit: %d MB)", (rss_req >> 20), (shm_avail >> 20), (mem_avail >> 20)) - return corenrn_datadir if not self._shm_enabled else corenrn_datadir_shm # - @@ -1512,10 +1462,6 @@ def clear_model(self, avoid_creating_objs=False, avoid_clearing_queues=True): round robin distribution, we want to clear the cells and synapses in order to have a clean slate on which to instantiate the balanced cells. """ - if not self._target_manager.parser: - # Target parser is the ground block. If not there model is clear - return - logging.info("Clearing model") self._pc.gid_clear() self._target_manager.clear_simulation_data() diff --git a/neurodamus/stimulus_manager.py b/neurodamus/stimulus_manager.py index eec7b188..b89c3cc1 100644 --- a/neurodamus/stimulus_manager.py +++ b/neurodamus/stimulus_manager.py @@ -56,7 +56,7 @@ def interpret(self, target_spec, stim_info): if SimConfig.cli_options.experimental_stims or stim_t and stim_t.IsPythonOnly: # New style Stim, in Python log_verbose("Using new-gen stimulus") - cell_manager = self._target_manager.hoc.cellDistributor + cell_manager = self._target_manager._cell_manager stim = stim_t(target, stim_info, cell_manager) self._stimulus.append(stim) else: diff --git a/neurodamus/target_manager.py b/neurodamus/target_manager.py index 449639f2..b6ee2c17 100644 --- a/neurodamus/target_manager.py +++ b/neurodamus/target_manager.py @@ -1,6 +1,5 @@ import itertools import logging -import os.path from abc import ABCMeta, abstractmethod from functools import lru_cache from typing import List @@ -8,8 +7,8 @@ import libsonata import numpy -from .core import MPI, NeurodamusCore as Nd -from .core.configuration import ConfigurationError, SimConfig, GlobalConfig, find_input_file +from .core import NeurodamusCore as Nd +from .core.configuration import ConfigurationError, find_input_file from .core.nodeset import _NodeSetBase, NodeSet, SelectionNodeSet from .utils import compat from .utils.logging import log_verbose @@ -92,13 +91,10 @@ def __init__(self, run_conf): Initializes a new TargetManager """ self._run_conf = run_conf - self.parser = Nd.TargetParser() - self._has_hoc_targets = False - self.hoc = None # The hoc level target manager + self._cell_manager = None self._targets = {} + self._section_access = {} self._nodeset_reader = self._init_nodesets(run_conf) - if MPI.rank == 0: - self.parser.isVerbose = 1 # A list of the local node sets self.local_nodes = [] @@ -114,8 +110,8 @@ def _init_nodesets(cls, run_conf): NodeSetReader(config_nodeset_file, simulation_nodesets_file) def load_targets(self, circuit): - """Provided that the circuit location is known and whether a user.target file has been - specified, load any target files via a TargetParser. + """Provided that the circuit location is known and whether a nodes file has been + specified, load any target files. Note that these will be moved into a TargetManager after the cells have been distributed, instantiated, and potentially split. """ @@ -123,56 +119,27 @@ def _is_sonata_file(file_name): if file_name.endswith(".h5"): return True return False - if circuit.CircuitPath: - self._try_open_start_target(circuit) nodes_file = circuit.get("CellLibraryFile") if nodes_file and _is_sonata_file(nodes_file) and self._nodeset_reader: self._nodeset_reader.register_node_file(find_input_file(nodes_file)) - def _try_open_start_target(self, circuit): - start_target_file = os.path.join(circuit.CircuitPath, "start.target") - if not os.path.isfile(start_target_file): - log_verbose("Circuit %s start.target not available! Skipping", circuit._name) - else: - self.parser.open(start_target_file, False) - self._has_hoc_targets = True - - def load_user_target(self): - # Old target files. Notice new targets with same should not happen - target_file = self._run_conf.get("TargetFile") - if not target_file or target_file.endswith(".json"): # allow any ext, except nodesets - return - user_target = find_input_file(target_file) - self.parser.open(user_target, True) - self._has_hoc_targets = True - if MPI.rank == 0: - logging.info(" => Loaded %d targets", self.parser.targetList.count()) - if GlobalConfig.verbosity >= 3: - self.parser.printCellCounts() - @classmethod def create_global_target(cls): - # In blueconfig mode the _ALL_ target refers to base single population) - if not SimConfig.is_sonata_config: - return _HocTarget(TargetSpec.GLOBAL_TARGET_NAME, None) # None population -> generic return NodesetTarget(TargetSpec.GLOBAL_TARGET_NAME, []) + def register_cell_manager(self, cell_manager): + self._cell_manager = cell_manager + def register_target(self, target): self._targets[target.name] = target - hoc_target = target.get_hoc_target() - if hoc_target: - self.parser.updateTargetList(target) def register_local_nodes(self, local_nodes): """Registers the local nodes so that targets can be scoped to current rank""" self.local_nodes.append(local_nodes) - self.parser.updateTargets(local_nodes.final_gids(), 1) def clear_simulation_data(self): self.local_nodes.clear() - self.parser.updateTargets(Nd.Vector(), 0) - self.init_hoc_manager(None) # Init/release cell manager def get_target(self, target_spec: TargetSpec, target_pop=None): """Retrieves a target from any .target file or Sonata nodeset files. @@ -206,40 +173,10 @@ def get_concrete_target(target): self.register_target(target) return get_concrete_target(target) - if self._has_hoc_targets: - if self.hoc is not None: - log_verbose("Retrieved `%s` from Hoc TargetManager", target_spec) - hoc_target = self.hoc.getTarget(target_name) - else: - log_verbose("Retrieved `%s` from the Hoc TargetParser", target_spec) - hoc_target = self.parser.getTarget(target_name) - target = _HocTarget(target_name, hoc_target) - self._targets[target_name] = target - return get_concrete_target(target) - raise ConfigurationError( "Target {} can't be loaded. Check target sources".format(target_name) ) - def init_hoc_manager(self, cell_manager): - # give a TargetManager the TargetParser's completed targetList - self.hoc = Nd.TargetManager(self.parser.targetList, cell_manager) - - def get_target_points(self, target, cell_manager, **kw): - """Helper to retrieve the points of a target. - If target is a cell then uses compartmentCast to obtain its points. - Otherwise returns the result of calling getPointList directly on the target. - - Args: - target: The target name or object (faster) - manager: The cell manager to access gids and metype infos - - Returns: The target list of points - """ - if isinstance(target, TargetSpec): - target = self.get_target(target) - return target.getPointList(cell_manager, **kw) - @lru_cache() def intersecting(self, target1, target2): """Checks whether two targets intersect""" @@ -269,9 +206,108 @@ def pathways_overlap(self, conn1, conn2, equal_only=False): return TargetSpec(src1) == TargetSpec(src2) and TargetSpec(dst1) == TargetSpec(dst2) return self.intersecting(src1, src2) and self.intersecting(dst1, dst2) - def __getattr__(self, item): - logging.debug("Compat interface to TargetManager::" + item) - return getattr(self.hoc, item) + def getPointList(self, target, **kw): + """Helper to retrieve the points of a target. + Returns the result of calling getPointList directly on the target. + + Args: + target: The target name or object + manager: The cell manager to access gids and metype infos + + Returns: The target list of points + """ + if not isinstance(target, NodesetTarget): + target = self.get_target(target) + return target.getPointList(self._cell_manager, **kw) + + def getMETypes(self, target_name): + """ + Convenience function for objects like StimulusManager to get access to METypes of cell + objects without having a direct line to the CellDistributor object. + + :param target_name: Target Name to get the GIDs and collect references to cell MEtypes + :return: List containing MEtypes for each cell object associated with the target + """ + result_list = compat.List() + target = self.get_target(target_name) + gids = target.get_local_gids() + + for gid in gids: + metype = self._cell_manager.getMEType(gid) + result_list.append(metype) + + return result_list + + def gid_to_sections(self, gid): + """ + For a given gid, return a list of section references stored for random access. + If the list does not exist, it is built and stored for future use. + + :param gid: GID of the cell + :return: SerializedSections object with SectionRefs to every section in the cell + """ + if gid not in self._section_access: + cell_ref = self._cell_manager.get_cellref(gid) + if cell_ref is None: + logging.warning("No cell found for GID: %d", gid) + return None + result_serial = SerializedSections(cell_ref) + self._section_access[gid] = result_serial + else: + result_serial = self._section_access[gid] + + return result_serial + + def location_to_point(self, gid, isec, ipt, offset): + """ + Given a location for a cell, section id, segment id, and offset into the segment, + create a list containing a section reference to there. + + :param gid: GID of the cell + :param isec: Section index + :param ipt: Distance to start of segment + :param offset: Offset distance beyond the ipt (microns) + :return: List with 1 item, where the synapse should go + """ + # Soma connection, just zero it + if offset < 0: + offset = 0 + + result_point = Nd.TPointList(gid) + cell_sections = self.gid_to_sections(gid) + if not cell_sections: + raise Exception("Getting locations for non-bg sims is not implemented yet...") + + if isec >= cell_sections.num_sections: + raise Exception(f"Error: section {isec} out of bounds ({cell_sections.num_sections} " + "total). Morphology section count is low, is this a good morphology?") + + distance = 0.5 + tmp_section = cell_sections.isec2sec[int(isec)] + + if tmp_section is None: # Assume we are in LoadBalance mode + result_point.append(None, -1) + elif ipt == -1: + # Sonata spec have a pre-calculated distance field. + # In such cases, segment (ipt) is -1 and offset is that distance. + offset = max(min(offset, 0.9999999), 0.0000001) + result_point.append(tmp_section, offset) + else: + # Adjust for section orientation and calculate distance + section = tmp_section.sec + if section.orientation() == 1: + ipt = section.n3d() - 1 - ipt + offset = -offset + + if ipt < section.n3d(): + distance = (section.arc3d(int(ipt)) + offset) / section.L + distance = max(min(distance, 0.9999999), 0.0000001) + + if section.orientation() == 1: + distance = 1 - distance + + result_point.append(tmp_section, distance) + return result_point class NodeSetReader: @@ -414,50 +450,7 @@ def update_local_nodes(self, _local_nodes): pass -class _HocTargetInterface: - """ - Interface of Hoc targets to be respected when we want to use objects - in place of hoc targets. - - This interface provides a default implementation suitable for Nodeset targets - where it will be primarily used - """ - - def isCellTarget(self, **kw): - section_type = kw.get("sections", "soma") - compartment_type = kw.get("compartments", "center" if section_type == "soma" else "all") - return section_type == "soma" and compartment_type == "center" - - def isCompartmentTarget(self, *_): - return 0 - - def isSectionTarget(self, *_): - return 0 - - def isSynapseTarget(self, *_): - return 0 - - def getCellCount(self): - return self.gid_count() - - def completegids(self): - return compat.hoc_vector(self.get_raw_gids()) - - def gids(self): - return compat.hoc_vector(self.get_local_gids()) - - @abstractmethod - def getPointList(self, _cell_manager, **_kw): - return NotImplemented - - def set_offset(self, *_): - pass # Only hoc targets require manually setting offsets - - def get_offset(self, *_): - pass # nodeset targets can span multiple multipopulation -> multiple offsets - - -class NodesetTarget(_TargetInterface, _HocTargetInterface): +class NodesetTarget(_TargetInterface): def __init__(self, name, nodesets: List[_NodeSetBase], local_nodes=None, **_kw): self.name = name self.nodesets = nodesets @@ -601,117 +594,31 @@ def generate_subtargets(self, n_parts): for cycle_i in range(n_parts)] -class _HocTarget(_TargetInterface): - """ - A wrapper around Hoc targets to implement _TargetInterface +class SerializedSections: + """ Serializes the sections of a cell for easier random access. + Note that this is possible because the v field in the section has been assigned + an integer corresponding to the target index as read from the morphology file. """ - GID_DTYPE = numpy.uint32 - - def __init__(self, name, hoc_target, pop_name=None, *, _raw_gids=None): - self.name = name - self.population_name = pop_name - self.hoc_target = hoc_target - self.offset = 0 - self._raw_gids = _raw_gids and numpy.array(_raw_gids, dtype=self.GID_DTYPE) - - @property - def population_names(self): - return {self.population_name} - - @property - def populations(self): - return {self.population_name: NodeSet(self.get_gids())} - - def gid_count(self): - return len(self.get_raw_gids()) - - def get_gids(self): - if not self.offset: - return self.get_raw_gids() - try: - return numpy.add(self.get_raw_gids(), self.offset, dtype=self.GID_DTYPE) - except numpy.core._exceptions.UFuncTypeError as e: - logging.error("Type error: please use type uint32 for the array of raw gids.") - raise e - - def get_raw_gids(self): - if self._raw_gids is None: - assert self.hoc_target - self._raw_gids = self.hoc_target.completegids().as_numpy().astype(self.GID_DTYPE) - self._raw_gids.sort() - return self._raw_gids - - def get_hoc_target(self): - return self.hoc_target - - def gids(self): - """This target gids on this rank, with offset""" - return self.hoc_target.gids() - - def get_local_gids(self): - return self.hoc_target.gids().as_numpy().astype(self.GID_DTYPE) - - def getPointList(self, cell_manager, **kw): - return self.hoc_target.getPointList(cell_manager) - - def make_subtarget(self, pop_name): - if pop_name is not None: - # Old targets have only one population. Ensure one doesn't assign more than once - if self.name == TargetSpec.GLOBAL_TARGET_NAME: # This target is special - return _HocTarget(pop_name + "__ALL__", Nd.Target(self.name), pop_name) - if self.population_name not in (None, pop_name): - raise ConfigurationError("Target %s cannot be reassigned population %s (cur: %s)" - % (self.name, pop_name, self.population_name)) - self.population_name = pop_name - return self - - def append_nodeset(self, nodeset: NodeSet): - # Not very common but we may want to set the nodes later (e.g. the _ALL_ target) - if self.population_name is not None: - logging.warning("[Compat] Skipping adding population %s to HOC target %s", - nodeset.population_name, self.name) - return - self.population_name = nodeset.population_name - self.offset = nodeset.offset - self._raw_gids = numpy.asarray(nodeset.raw_gids()) - hoc_gids = compat.hoc_vector(self._raw_gids) - self.hoc_target = Nd.Target(self.name, hoc_gids, self.population_name) - - def __contains__(self, item): - return self.hoc_target.completeContains(item) - - def is_void(self): - return not bool(self.hoc_target) # old targets could match with any population - - def set_offset(self, offset): - self.hoc_target.set_offset(offset) - self.offset = offset - - def isCellTarget(self, **kw): - return self.hoc_target.isCellTarget() - - def generate_subtargets(self, n_parts): - """generate sub hoc targets for multi-cycle runs - Returns: - list of subtargets - """ - if not n_parts or n_parts == 1: - return False - - allgids = self.get_gids() - new_targets = [] - - for cycle_i in range(n_parts): - target = Nd.Target() - target.name = "{}_{}".format(self.name, cycle_i) - new_targets.append(_HocTarget(target.name, target, self.population_name)) - - target_looper = itertools.cycle(new_targets) - for gid in allgids: - target = next(target_looper) - target.hoc_target.gidMembers.append(gid) - - return new_targets - - def __getattr__(self, item): - return getattr(self.hoc_target, item) + def __init__(self, cell): + self.num_sections = int(cell.nSecAll) + # Initialize list to store SectionRef objects + self.isec2sec = [None] * self.num_sections + # Flag to control warning message display + self._serialized_sections_warned = False + + index = 0 + for sec in cell.all: + # Accessing the 'v' value at location 0.0001 of the section + v_value = sec(0.0001).v + if v_value >= self.num_sections: + logging.debug("{%s} v(1)={sec(1).v} n3d()={%f}", sec.name(), sec.n3d()) + raise Exception("Error: failure in mk2_isec2sec()") + if v_value < 0: + if not self._serialized_sections_warned: + logging.warning("SerializedSections: v(0.0001) < 0. index={%d} v()={%f}", + index, v_value) + self._serialized_sections_warned = True + else: + # Store a SectionRef to the section at the index specified by v_value + self.isec2sec[int(v_value)] = Nd.SectionRef(sec=sec) + index += 1 diff --git a/tests/integration-e2e/test_multicycle_runs.py b/tests/integration-e2e/test_multicycle_runs.py index 25ea92eb..f291d1a6 100644 --- a/tests/integration-e2e/test_multicycle_runs.py +++ b/tests/integration-e2e/test_multicycle_runs.py @@ -41,23 +41,6 @@ def test_nodeset_target_generate_subtargets(): assert np.array_equal(subtargets[2][1].get_gids(), np.array([1002])) -def test_hoc_target_generate_subtargets(): - from neurodamus.target_manager import _HocTarget - - N_PARTS = 3 - gids = list(range(10)) - target = _HocTarget("Column", None, pop_name=None) - target._raw_gids = np.array(gids, dtype="uint32") - - subtargets = target.generate_subtargets(N_PARTS) - assert len(subtargets) == N_PARTS - assert subtargets[0].name == "Column_0" - assert subtargets[0].population_names == {None} - assert np.array_equal(subtargets[0].get_gids(), np.array([0, 3, 6, 9])) - assert np.array_equal(subtargets[1].get_gids(), np.array([1, 4, 7])) - assert np.array_equal(subtargets[2].get_gids(), np.array([2, 5, 8])) - - def _create_tmpconfig_coreneuron(config_file): import fileinput import shutil diff --git a/tests/unit/test_target_intersection.py b/tests/unit/test_target_intersection.py index dc9b9687..819b7093 100644 --- a/tests/unit/test_target_intersection.py +++ b/tests/unit/test_target_intersection.py @@ -1,4 +1,3 @@ -import numpy as np import numpy.testing as npt import pytest @@ -39,26 +38,6 @@ def test_targetspec_overlap(): assert not TargetSpec("pop1:t1").overlap(TargetSpec("pop1:t2")) -@pytest.mark.forked -def test_hoc_target_intersect(): - from neurodamus.target_manager import _HocTarget - ht1 = _HocTarget("t1", None, pop_name=None) - ht2 = _HocTarget("t2", None, pop_name=None) - ht1_pop1 = _HocTarget("t1", None, pop_name="pop1") - ht1_pop2 = _HocTarget("t1", None, pop_name="pop2") - - # different population is short circuited -> gids aren't necessary to see they dont intersect - assert not ht1_pop1.intersects(ht1_pop2) - # We override the internal gid cache to avoid creating hoc targets - ht1._raw_gids = np.array([1, 2], dtype="uint32") - ht2._raw_gids = np.array([1, 2], dtype="uint32") - assert ht1.intersects(ht2) - ht2._raw_gids = np.array([2, 3], dtype="uint32") - assert ht1.intersects(ht2) - ht2._raw_gids = np.array([3, 4], dtype="uint32") - assert not ht1.intersects(ht2) - - @pytest.mark.forked def test_nodeset_target_intersect(): from neurodamus.target_manager import NodesetTarget @@ -92,7 +71,7 @@ def test_nodeset_gids(): [nodes_popA, nodes_popB], local_nodes=[local_nodes_popA, local_nodes_popB] ) - gids = t.gids() + gids = t.get_local_gids() npt.assert_array_equal(gids, [5, 6, 1003, 1004, 1005]) t2 = NodesetTarget( @@ -100,7 +79,7 @@ def test_nodeset_gids(): [nodes_popA, nodes_popB], local_nodes=[local_nodes_popA] ) - gids = t2.gids() + gids = t2.get_local_gids() npt.assert_array_equal(gids, [5, 6]) t3 = NodesetTarget( @@ -108,7 +87,7 @@ def test_nodeset_gids(): [nodes_popB], local_nodes=[local_nodes_popA, local_nodes_popB] ) - gids = t3.gids() + gids = t3.get_local_gids() npt.assert_array_equal(gids, [1003, 1004, 1005]) t4 = NodesetTarget( @@ -116,5 +95,5 @@ def test_nodeset_gids(): [nodes_popB], local_nodes=[local_nodes_popA] ) - gids = t4.gids() + gids = t4.get_local_gids() npt.assert_array_equal(gids, [])