From 12634eff25ad3c66d39af4d567c6f6c67b6dfcc6 Mon Sep 17 00:00:00 2001 From: WeinaJi Date: Thu, 29 Feb 2024 13:28:04 +0100 Subject: [PATCH] [BBPBGLIB-1132]Replace StimulusManager.hoc with the python classes (#132) ## Context Deprecate the hoc implementation of the stimulus manager. ## Scope Replace StimulusManager.hoc by the exsiting stimulus_manager.py where all the stimuli have already been defined in the dedicated python classes, except the type `extracellular_stimulation` (pending for BBPBGLIB-890). ## Testing New unit test in `test_sonata_config.py` for parsing "subthreshold" stimulus. ## Review * [x] PR description is complete * [x] Coding style (imports, function length, New functions, classes or files) are good * [x] Unit/Scientific test added * [ ] Updated Readme, in-code, developer documentation --- core/hoc/EStim.hoc | 197 ---- core/hoc/StimulusManager.hoc | 1263 ----------------------- core/hoc/TDistFunc.hoc | 82 -- core/hoc/TStim.hoc | 248 ----- core/hoc/neurodamus.hoc | 3 +- neurodamus/io/sonata_config.py | 8 +- neurodamus/node.py | 3 +- neurodamus/stimulus_manager.py | 55 +- tests/integration/test_sonata_config.py | 11 + 9 files changed, 37 insertions(+), 1833 deletions(-) delete mode 100644 core/hoc/EStim.hoc delete mode 100644 core/hoc/StimulusManager.hoc delete mode 100644 core/hoc/TDistFunc.hoc delete mode 100644 core/hoc/TStim.hoc diff --git a/core/hoc/EStim.hoc b/core/hoc/EStim.hoc deleted file mode 100644 index bbba114e..00000000 --- a/core/hoc/EStim.hoc +++ /dev/null @@ -1,197 +0,0 @@ -/** - * @file EStim.hoc - * @brief - * @author reimann - * @date 2010-12-30 - * @remark Copyright © BBP/EPFL 2005-2011; All rights reserved. Do not distribute without further notice. - */ - -D__ESTIMHOC = 1 - - - /** - * EStim class is a derivative of TStim for stimuli with an extracelular electrode. The main difference is that it - * collects all elementary stimuli pulses and converts them using a VirtualElectrode object before it injects anything - * - * The stimulus is defined on the hoc level by using the addpoint function for every (step) change in extracellular - * electrode voltage. At this stage only step changes can be used. Gradual, i.e. sinusoidal changes will be - * implemented in the future - * After every step has been defined, you have to call initElec() to perform the frequency dependent transformation. - * This transformation turns e_electrode into e_extracellular at distance d=1 micron from the electrode. - * After the transformation is complete, NO MORE STEPS CAN BE ADDED! You can then use inject() to first scale - * e_extracellular down by a distance dependent factor and then vector.play() it into the currently accessed compartment - * - * TODO: 1. more stimulus primitives than step. 2. a dt of 0.1 ms is hardcoded. make this flexible! - */ - -begintemplate EStim -objref tvec, stimVec, elecObj, this -objref stim,rand,null -public ic,pulse,ramp,stim,tvec,delay,segment,sinspec,noise,x,c,loc,mark,get_loc, train, initElec, initRawElec, xPos, inject - -/** - * - */ -proc init(/*x[,Random]*/){ - stim = new Vector() //where the new voltage of every step change is stored - tvec = new Vector() //times of the steps. Formal: e_electrode = stim[i] for t in [tvec[i], tvec[i+1]) - stimVec = new List() // this is where the transformed and injected Vector objects are stored - if(numarg()>0 && argtype(1)!=0) rand = $o1 //At this stage I use no random stimuli... reserved for future versions -} - -/** - * initElec(): Performs frequency dependent transform using a VirtualElectrode object. - * The VirtualElectrode takes the time/amplitude of the step changes and returns the signal in an unfortunately - * hardcoded frequency of 10000Hz. - * - * - */ -proc initElec(){ local i, numElem, eValue, startValue - //Instantiation of the VirtualElectrode object wrapped in an execute statement because for some machines - //(linsrv1 and sg1 at the moment) the VirtualElectrode.mod cannot be compiled due to missing fftw library. - //For these machines we just remove the .mod file (disabling extracellular stimulation!) - execute("elecObj = new VirtualElectrode()",this) - //elecObj = new VirtualElectrode() - elecObj.setLgth(tvec.size()+2) //number of step changes. Add 2 to the end of it to make sure that it really is 0 at the end - if(stim.size()==0)return //no steps: not much to do. - startValue = tvec.x[0] - for i=0, tvec.size()-1 { - elecObj.addPoint(tvec.x[i],stim.x[i]) // fill values into the VirtualElectrode - } - elecObj.addPoint(tvec.x[tvec.size()-1]+0.1,0) //Make sure it goes down to 0 - elecObj.addPoint(tvec.x[tvec.size()-1]+20,0) //arbitrary 20 is arbitrary. Maybe use Duration instead? - numElem = tvec.x[tvec.size()-1]+20 //We need to know until which time we have to extract samples - tvec.resize(0) //empty tvec and stim. These will be re-used to hold the 10000 Hz sampled signal instead - stim.resize(0) - elecObj.convert() - i = 0.0 - while(i<=numElem) { //now start filling stim with the converted signal - eValue = elecObj.getValue() - if(i >= startValue){ //To make sure that nothing happens before the stimulus starts. elecObj.getValue() can be >0 due to numerical issues - stim.append(eValue) - } else { - stim.append(0) - } - tvec.append(i) - i = i+0.1 - } -} - -proc initRawElec(){ local j, tm, vl, nextTime localobj stimPts, tvecPts - stimPts = stim.c - tvecPts = tvec.c - stim.resize(0) - tvec.resize(0) - tm = 0.0 - vl = 0 - j = 0 - nextTime = tvecPts.x[j] - while( tm < tvecPts.x[tvecPts.size()-1]+20) { - if((nextTime > 0) && (nextTime <= tm)){ - vl = stimPts.x[j] - j = j+1 - if(j < tvecPts.size()){ - nextTime = tvecPts.x[j] - } else { - nextTime = -1 - } - } - tvec.append(tm) - tm = tm+0.1 - stim.append(vl) - } -} - -/* - * inject() CALL ONLY AFTER initElect(). TODO: actually check if it has been initiated - * injects a scaled down (scale depends on distance) version of the converted signal into the active compartment - * - * */ - -proc inject() { local i localobj vecCopy - vecCopy = stim.c //copy... - //scale down... - vecCopy.mul($2) - vecCopy.play(&e_extracellular($1),tvec,1) //and inject - stimVec.append(vecCopy) -} - -//basically used to start and end a stimulus -proc addpoint(/*amp,tstart*/){ local i - i=0 - while(i < tvec.size()){ - if($2 < tvec.x[i]) break - if($2 == tvec.x[i]){ - if(abs($1)>abs(stim.x[i])) stim.x[i] = $1 - return - } - i=i+1 - } - tvec.insrt(i,$2) - stim.insrt(i,$1) -} - - -proc pulse(/*start,dur,max_amp[,min_amp]*/){local min_amp -min_amp = 0 -if(numarg()>3) min_amp = $4 -addpoint($3,$1) -addpoint(min_amp,$1+$2) -} - - -/** - * Stimulus with repeated pulse injections at a specified frquency - * - * @param $1 del - time when stimulation starts - * @param $2 dur - time when stimulation ends - * @param $3 amp - the amplitude of a each pulse - * @param $4 frequency - determines the number of pulses per second (hz) - * @param $5 width - the duration of a single pulse - */ -proc train() { local tpulse, remtime, width - - tpulse = 1000/$4 - width = $5 - - if( width > tpulse ) { //Duration of a single pulse is longer than the time between pulses -> just use a standard pulse - pulse( $1, $2, $3 ) - return - } - - remtime = $2 - nextpulse = $1 - -//TODO: for better clarity, don't try to shift pulse onsets by one dt (only done to match old bglib) - - while( remtime > 0 ) { - if( width < remtime ) { //can fit whole pulse in remaining time - addpoint( $3, nextpulse ) - addpoint( 0, nextpulse+width ) - } else { //final pulse, just paritial - addpoint( $3, nextpulse ) - addpoint( 0, $1+$2 ) - } - - //TODO: determine if we should stick with old bglib way, where time spent in pulse doesn't count against tpulse - //nextpulse = nextpulse + tpulse - nextpulse = nextpulse + tpulse + width //by adding time of width, tpulse is kept full - remtime = remtime - tpulse // width doesn't count against remtime -} - - -/* -for remtime=0, tvec.size()-1 { - print tvec.x[remtime], " ", stim.x[remtime] -} - */ - -} - - -proc delay(/*dur*/){ - tvec.add($1) -} - - -endtemplate EStim diff --git a/core/hoc/StimulusManager.hoc b/core/hoc/StimulusManager.hoc deleted file mode 100644 index dd803cee..00000000 --- a/core/hoc/StimulusManager.hoc +++ /dev/null @@ -1,1263 +0,0 @@ -/** - * @file StimulusManager.hoc - * @brief Process Synapse connections Rules of a BlueConfig (if any) and then instantiate the synapses for a bluebrain sim - * @author king - * @date 2009-10-23 - * @remark Copyright © BBP/EPFL 2005-2011; All rights reserved. Do not distribute without further notice. - */ - -{load_file("TStim.hoc")} //note that the TStim functions are used in a way where our naming schemes don't match 100% -{load_file("EStim.hoc")} -{load_file("fileUtils.hoc")} -{load_file("Map.hoc")} - -begintemplate StimulusManager - -external terminate - -//----------------------------------------------------------------------------------------------- -// Declare member variables -//----------------------------------------------------------------------------------------------- - -objref stimList, elecManager, targetManager, cellDistributor, this - -//! For use in poisson stims, when creating via the execute stmt, need to have class level object (not just local to function) -objref stimHelper, poissonNetconList, handleList, poissonSynList, vecStimList, voltageReplayList, timeReplayList, synapseEventGenerator - -//----------------------------------------------------------------------------------------------- -// public members -//----------------------------------------------------------------------------------------------- - -public init, linear, interpret, FillSpikeTrain, interpretExtracellulars, saveStatePreparation, reevent - -//note: for some reason, if we do not make 'noise' function public, we get a syntax error. -public noise - -//----------------------------------------------------------------------------------------------- -// Member function implementations -//----------------------------------------------------------------------------------------------- - -/*! - * Typical constructor, create empty list where stimulus are stored as they get generated, keep reference - * to TargetManager for accessing Targets - * - * @param $o1 TargetManager - * @param $o2 Electrode Manager for extracellular stimulation (optional) - * @param $3 [deprecated] optional globalBaseSeed to apply to random numbers (default 0) - */ -proc init() { local argIndex - - stimList = new List() - targetManager = $o1 - - //keep track of how many poisson and noise stims we create so that the seeds can be offset (otherwise we would have same injections or syns/spike trains in a given cell) - poissonStimCount = 0 - noiseStimCount = 0 - poissonNetconList = new List() - poissonSynList = new List() - handleList = new List() - vecStimList = new List() - voltageReplayList = new List() - timeReplayList = new List() - synapseEventGenerator = new List() - - // if we have been passed an elec manager or base seed argument, need to determine now by arg type. May want to clean this up in the future to - // remove excess freedom for the users - if( numarg() > 1 ) { - if( argtype(2) == 0 ) { // 0 indicates value type arg, which goes to baseSeed - globalBaseSeed = $2 - } else if ( argtype(2) == 1 ) { // 1 indicates objref, which goes to elecManager - elecManager = $o2 - } - } - if( numarg() > 2 ) { - if( argtype(3) == 0 ) { // 0 indicates value type arg, which goes to baseSeed - globalBaseSeed = $3 - } else if ( argtype(3) == 1 ) { // 1 indicates objref, which goes to elecManager - elecManager = $o3 - } - } -} - -//----------------------------------------------------------------------------------------------- - -/*! - * Helper function to get an overview over the requested extracellular stimuli. interpretExtracellulars() needs - * them in the following format: A map(String,List), where the key is the target name and the object is - * a List of Maps containing Stimulus parameters. Because it is a rather complicated structure to set up, I do - * it in a separate function. - * - * @param $o1 injRequests - the list of maps containing info on requested StimulusInjects - * @param $o2 parsedStimuli - Maps containing the information on defined stimuli - * - */ -obfunc getExtracellularParams() { local injIndex localobj listForTarget, newListObj, injRequests, stimName, targetName, stim, paramValue, extraList, nil - if($o1 == nil){ - print "injRequests are NULL!" - } - if($o2 == nil){ - print "parsedStimuli are NULL!" - } - injRequests = $o1 - extraList = new Map() - for injIndex=0, injRequests.count()-1 { - stimName = injRequests.o(injIndex).get( "Stimulus" ) - targetName = injRequests.o(injIndex).get( "Target" ) - stim = $o2.get( stimName ) - paramValue = stim.get( "Mode" ) - if( strcmp( paramValue.s, "Extracellular" ) == 0 ) { - listForTarget = extraList.get(targetName) - if( listForTarget == nil ){ - newListObj = new List() - newListObj.append(stim) - extraList.put( targetName, newListObj ) - } else { - listForTarget.append(stim) - } - } - } - return extraList -} - -//----------------------------------------------------------------------------------------------- - -/*! - * This function handles all the requested extracellular stimuli. It uses helper function getExtracellularParams - * to get an overview over requested extracellular stimuli - * - * - * @param $o1 injRequests - the list of maps containing info on requested StimulusInjects - * @param $o2 parsedStimuli - Maps containing the information on defined stimuli - * - */ -proc interpretExtracellulars() { local cellIndex, targetIndex, stimIndex, elecIndex, scaleFac, secType, secNum, offset, tx, ty, tz, x localobj stimlist, eName, pts, gidObj, stim, estim, paramValue, activeElectrode, nil - strdef tstr - //verify that electrode manager was passed to StimulusManager - if( elecManager == nil ) { - terminate("Error: attempting to use extracellular stimulation without ElectrodeManager" ) - } - - stimlist = getExtracellularParams($o1,$o2) - for targetIndex=0, stimlist.count()-1{ //over all unique targets - //print "Treating extracellulars for target: ", stimlist.key(targetIndex).s - pts = targetManager.getPointList(stimlist.key(targetIndex).s) //get that target - gidObj = targetManager.gids(stimlist.key(targetIndex).s) - - for elecIndex = 0, elecManager.count()-1{ //over all electrodes - //Now we have every combination of electrode and target. We collect all Stimuli using a given combination, - //Fill them into one EStim object and do the frequency dependent conversion. Then we Vector.play() it into - //Every compartment of the target - eName = elecManager.getName(elecIndex) - activeElectrode = elecManager.getElectrode(elecIndex) - estim = new EStim() //Create one EStim per electrode! - for stimIndex=0, stimlist.o(targetIndex).count()-1 { //Now gather every stimulus using this electrode - stim = stimlist.o(targetIndex).o(stimIndex) - // No longer checking for type of stimulus at this stage! - paramValue = stim.get( "Electrode" ) - //Compare Electrode name. TODO: Need to find a solution for lookup by location. - //Also: Maybe do this with a Map object? - if( paramValue != nil && (strcmp(paramValue.s,eName.s)==0) ){ - estim.train( stim.valueOf( "Delay" ), stim.valueOf( "Duration" ), stim.valueOf( "AmpStart" ), stim.valueOf( "Frequency" ), stim.valueOf( "Width" ) ) - } - } - estim.initRawElec() //now that we treated all stimuli using this electrode we can initiate it - for cellIndex=0, pts.count()-1 { //Now over all compartments in the target - scaleFac = 0 - for pts.o(cellIndex).each_point(&x) { - sprint(tstr,"{%s insert extracellular}",secname()) - execute1(tstr,0) - if(!getSecTypeAndIndex(&secType,&secNum)){ //Now get the scale factor based on distance - if(activeElectrode.vInfo()>4){ //version 4 has one value per compartment, <4 only one per Section - scaleFac = activeElectrode.getValueForSection(gidObj.x[cellIndex],secType,secNum,x) - } else { - scaleFac = activeElectrode.getValueForSection(gidObj.x[cellIndex],secType,secNum) - } - if(scaleFac == 0){ //No need to inject anything. 0 is default value for e_extracellular - continue - } - } else { - terminate( "Problem identifying section type and number for ", secname() ) - } - estim.inject(x,scaleFac) - } - } - stimList.append(estim) - } - } -} - -//----------------------------------------------------------------------------------------------- -//----------------------------------------------------------------------------------------------- - -/*! - * For handling stimulus which has been parsed from a file - look at the parameters collected and - * call the appropriate stimulus functions - * - * @param $s1 targetName - which target where stimulus should be injected - * @param $o2 stim - map containing stimulus parameters - * - * alternatively for Ruben's PFE - * @param $o1 TPointList where stimulus will be injected - * @param $o2 stim - map containing stimulus parameters - * @param $3 gid of cell - */ -proc interpret() { local tempVal, del, dur, relamp, threshold, targetVm localobj stim, target, paramValue, paramVect, patternName, synType, synConfig, fillFunc, nil - - stim = $o2 - - paramVect = new Vector() - paramValue = stim.get( "Mode" ) - - //Extracellular stimuli handled separately; Current and Voltage Modes only - if( strcmp( paramValue.s, "Extracellular" ) == 0 ) { - return - } - - //all will have Delay and Duration, so parse immediately - del = stim.valueOf( "Delay" ) - dur = stim.valueOf( "Duration" ) - - if( t > 0 ) { - del = del + t - } - - //more specific parameters based on pattern - patternName = stim.get( "Pattern" ) - - // Ruben's pfe chunk - // - // check for TPointList argument - if( argtype(1) == 1 ) { - // SquarePulse for TPointList arg - if( strcmp( patternName.s, "SquarePulse" ) == 0 ) { - relamp = stim.valueOf( "relAmp" ) - //threshold = targetManager.getCell($3).getCCell().getThreshold() - threshold = targetManager.getCell($3).getThreshold() - squarepulse( $o1, relamp * threshold, del, dur ) - } - // - // check for target name argument - } else if( argtype(1) == 2 ) { - // SquarePulse for target name arg - if( strcmp( patternName.s, "SEClamp" ) == 0 ) { - targetVm = stim.valueOf( "Voltage" ) - seclamp( $s1, targetVm, del, dur ) - } else if( strcmp( patternName.s, "Hyperpolarizing" ) == 0 ) { - hypamp( $s1, del, dur ) - } else if( strcmp( patternName.s, "SquarePulse" ) == 0 ) { - /* - relamp = stim.valueOf( "relAmp" ) - threshold = targetManager.getCells($s1).o(0).getCCell().getThreshold() - linear( $s1, relamp * threshold, relamp * threshold, del, dur ) - targetManager.getCells($s1) - */ - } else if( strcmp( patternName.s, "Linear" ) == 0 ) { - paramVect.append( stim.valueOf("AmpStart") ) - - //optional option - AmpEnd. If not there, use same value as AmpStart - if( stim.exists( "AmpEnd" ) ) { - paramVect.append( stim.valueOf("AmpEnd") ) - } else { - paramVect.append( paramVect.x[paramVect.size()-1] ) - } - - linear( $s1, paramVect.x[0], paramVect.x[1], del, dur ) - } else if( strcmp( patternName.s, "RelativeLinear" ) == 0 ) { - paramVect.append( stim.valueOf("PercentStart") ) - - //optional option - PercentEnd. If not there, use same value as PercentStart - if( stim.get( "PercentEnd" ) != nil ) { - paramVect.append( stim.valueOf("PercentEnd") ) - } else { - paramVect.append( paramVect.x[paramVect.size()-1] ) - } - - relativeLinear( $s1, paramVect.x[0], paramVect.x[1], del, dur ) - } else if ( strcmp( patternName.s, "Pulse" ) == 0 ) { - //looking at bglib, it expects potentially multiple delay/duration pairs (and an AmpEnd) - // TODO: rename bglib's stimulus pattern names to be inline with TStim, creating 'train' as replacement for bglib's 'pulse' - - paramVect.append( stim.valueOf( "AmpStart" ) ) - paramVect.append( stim.valueOf( "Frequency" ) ) - paramVect.append( stim.valueOf( "Width" ) ) - - //optional argument - Offset; default = 0. - // TODO: re-evaluate offset for Pulse stimulus and debug because it looks like it might not be multisplit safe - if( stim.exists( "Offset" ) ) { - paramVect.append( stim.valueOf( "Offset" ) ) - } else { - paramVect.append( 0 ) - } - - train( $s1, paramVect.x[0], paramVect.x[1], paramVect.x[2], paramVect.x[3], del, dur ) - //linear( $s1, paramVect.x[0], paramVect.x[0], del, dur ) - } else if(strcmp(patternName.s, "Sinusoidal")==0){ - paramVect.append( stim.valueOf("AmpStart")) - paramVect.append( stim.valueOf("Frequency")) - if( stim.exists( "Offset" ) ) { - paramVect.append( stim.valueOf( "Offset" ) ) - } else { - paramVect.append( 0 ) - } - sinusoid($s1, paramVect.x[0], paramVect.x[1], paramVect.x[2], del, dur) - } else if( strcmp( patternName.s, "NPoisson" ) == 0 ) { - paramVect.append( stim.valueOf( "Lambda" ) ) - paramVect.append( stim.valueOf( "Weight" ) ) - paramVect.append( stim.valueOf( "NumOfSynapses" ) ) - - //Note: SynapseType is optional and may return nil (handled by npoisson function) - //Note: SynapseConfigure is optional and may return nil (handled by npoisson function) - //Note: FillFunctionOverride is optional and under development, may return nil (handled by npoisson function) - //Note: requires special Target handling, so just pass targetName through - if( stim.exists( "SynapseType" ) ) { synType = stim.get( "SynapseType" ) } - if( stim.exists( "SynapseConfigure" ) ) { synConfig = stim.get( "SynapseConfigure" ) } - if( stim.exists( "FillFunctionOverride" ) ) { fillFunc = stim.get( "FillFunctionOverride" ) } - npoisson( $s1, paramVect.x[0], paramVect.x[1], paramVect.x[2], synType, synConfig, fillFunc, del, dur ) - } else if( strcmp( patternName.s, "NPoissonInhomogeneous" ) == 0 ) { - //Note: SynapseType is optional and may return nil (handled by npoisson function) - //Note: SynapseConfigure is optional and may return nil (handled by npoisson function) - - //Note: even though we dont use delay/duration, they must still be present in the BlueConfig until I can refactor the code better - //Note: requires special Target handling, so just pass targetName through - if( stim.exists( "SynapseType" ) ) { synType = stim.get( "SynapseType" ) } - if( stim.exists( "SynapseConfigure" ) ) { synConfig = stim.get( "SynapseConfigure" ) } - npoissonInhomogeneous( $s1, stim.valueOf( "Weight" ), stim.valueOf( "NumOfSynapses" ), synType, \ - synConfig, stim.get( "H5File" ), stim.get( "TBins" ), stim.get( "Rate" ) ) - } else if( strcmp( patternName.s, "SubThreshold" ) == 0 ) { - subthreshold( $s1, stim.valueOf( "PercentLess" ), del, dur ) - } else if( strcmp( patternName.s, "Noise" ) == 0 ) { - if( stim.exists( "Mean" ) ) { - noise( $s1, stim.valueOf( "Mean" ), stim.valueOf( "Variance" ), del, dur ) - } else { - relativeNoise( $s1, stim.valueOf( "MeanPercent" ), stim.valueOf( "Variance" ), del, dur ) - } - } else if ( strcmp( patternName.s, "ReplayVoltageTrace" ) == 0 ) { - //will replay a voltage trace into a cell's soma via SEClamp - if( stim.exists( "File" ) ) { - //expecting path to hdf5 files - // TODO: add suport for other file formats - replayVoltageTrace( $s1, stim.get( "File" ).s, stim.get( "DataSetLabel" ).s, del, dur ) - } - } else { - terminate("Unknown stimulus type", patternName.s) - } - } -} - -//----------------------------------------------------------------------------------------------- - -/*! - * Injects a squarepulse stimulus into a cell. Takes a TPointList, not a target, as first argument - * - * @param $o1 TPointList of points where iclamps will be placed - * @param $2 amp amplitude of the current injected (nA) - * @param $3 delay from start of simulation to when current is injected (ms) - * @param $4 duration of stimulus (ms) - */ -proc squarepulse() { local cellIndex, x, amp, del, dur localobj pts, tstim - - amp = $2 - del = $3 - dur = $4 - - //since there are no random elements, we can just evaluate the pointlist and proceed as normal - pts = $o1 - - for pts.each_point(&x) { - - tstim = new TStim(x) - - tstim.pulse( del, dur, amp) - stimList.append( tstim ) - } -} - -//----------------------------------------------------------------------------------------------- - -/** - * Apply a cell's hyperpolarizing current for the given time span. The current value (I) to apple is - * part of the CCell template (provided it has been calculated) - * - * @param $s1 Target name - * @param $2 delay - * @param $3 duration - */ -proc hypamp() { local hypampI, ret, x, del, dur localobj pts, cellList, tstim - pts = targetManager.getPointList($s1) - cellList = targetManager.getMETypes($s1) - //cellList = targetManager.getCells($s1) - del = $2 - dur = $3 - - for cellIndex=0, pts.count()-1 { - for pts.o(cellIndex).each_point(&x) { - // Now with MType instances this always works - hypampI = cellList.o(cellIndex).getHypAmp() - tstim = new TStim(x) - tstim.pulse( del, dur, hypampI ) - stimList.append(tstim) - } - } -} - -//----------------------------------------------------------------------------------------------- - -/** - * Uses TStim noise functionality to give each TPoint a random current injection. Like noise stimulus, but uses each cell's - * threshold firing current and scales it according to a percentage. Only tested for Cell targets - * TODO: There seems to be an issue where it continues to inject some stimulus after the duration should expire. - * For now, we use it for background stimulation over the whole runtime, so it is low priority to investigate - * - * @param $s1 targetName whose points will receive the stimulation - * @param $2 mean percentage of a cell's threshold current - the noise injection will center around this value - * @param $3 variance of current - * @param $4 delay from start of simulation to when current is injected (ms) - * @param $5 duration of stimulus (ms) - */ -proc relativeNoise() { local x, threshold localobj pts, cellList, rand, tstim, rngInfo - pts = targetManager.getPointList($s1) - cellList = targetManager.getMETypes($s1) - - rngInfo = new RNGSettings() - - for cellIndex=0, pts.count()-1 { - for pts.o(cellIndex).each_point(&x) { - threshold = cellList.o(cellIndex).getThreshold() - - if( rngInfo.getRNGMode() == rngInfo.COMPATIBILITY ) { - rand = new Random( pts.o(cellIndex).gid + noiseStimCount ) - } else if( rngInfo.getRNGMode() == rngInfo.UPMCELLRAN4 ) { - rand = new Random() - rand.MCellRan4( noiseStimCount*10000+100, rngInfo.getGlobalSeed() + rngInfo.getStimulusSeed() + pts.o(cellIndex).gid*1000 ) - } else if( rngInfo.getRNGMode() == rngInfo.RANDOM123 ) { - rand = new Random() - rand.Random123( noiseStimCount+100, rngInfo.getStimulusSeed() + 500, pts.o(cellIndex).gid + 300 ) - } - tstim = new TStim(x, rand) - tstim.noise( $4, $5, threshold*$2/100, threshold*$3/100 ) - stimList.append(tstim) - } - } - - noiseStimCount = noiseStimCount+1 -} - -//----------------------------------------------------------------------------------------------- - -/** - * Uses TStim noise functionality to give each TPoint a random current injection. Only tested for Cell targets - * TODO: There seems to be an issue where it continues to inject some stimulus after the duration should expire. - * For now, we use it for background stimulation over the whole runtime, so it is low priority to investigate - * - * @param $s1 targetName whose points will receive the stimulation - * @param $2 mean current where the noise should be centered - * @param $3 variance of current - * @param $4 delay from start of simulation to when current is injected (ms) - * @param $5 duration of stimulus (ms) - */ -proc noise() { local x localobj pts, rand, tstim, rngInfo - pts = targetManager.getPointList($s1) - - rngInfo = new RNGSettings() - - for cellIndex=0, pts.count()-1 { - for pts.o(cellIndex).each_point(&x) { - if( rngInfo.getRNGMode() == rngInfo.COMPATIBILITY ) { - rand = new Random( pts.o(cellIndex).gid + noiseStimCount ) - } else if( rngInfo.getRNGMode() == rngInfo.UPMCELLRAN4 ) { - rand = new Random() - rand.MCellRan4( noiseStimCount*10000+100, rngInfo.getGlobalSeed() + rngInfo.getStimulusSeed() + pts.o(cellIndex).gid*1000 ) - } else if( rngInfo.getRNGMode() == rngInfo.RANDOM123 ) { - rand = new Random() - rand.Random123( noiseStimCount+100, rngInfo.getStimulusSeed() + 500, pts.o(cellIndex).gid + 300 ) - } - tstim = new TStim(x, rand) - tstim.noise( $4, $5, $2, $3 ) - stimList.append(tstim) - } - } - - noiseStimCount = noiseStimCount+1 -} - -//----------------------------------------------------------------------------------------------- - -/** - * Injects a pulse stimulus into the targeted locations scaled to some percentage of that cell's threshold - * current (the current needed to cause 1 AP) - * - * @param $s1 Target name which should receive the stimulus - * @param $2 Percent less than full threshold current (e.g. val of 5 indicates that 95% of threshold is to be injected - * @param $3 del - * @param $4 dur - */ -proc subthreshold() { local cellIndex, x localobj pts, cellList, tstim - - pts = targetManager.getPointList($s1) - cellList = targetManager.getMETypes($s1) - //cellList = targetManager.getCells($s1) - - for cellIndex=0, pts.count()-1 { - for pts.o(cellIndex).each_point(&x) { - tstim = new TStim(x) - - tstim.pulse( $3, $4, cellList.o(cellIndex).getThreshold()*(100-$2)/100 ) - stimList.append(tstim) - } - } -} - -//----------------------------------------------------------------------------------------------- - -/*! - * A modified version of the liner stimulus where the injected current is relative to a cell's threshold - * current. Note, unlike the subthreshold stimulus, this uses the percent value directly rather than - * subtracting from 100%. - * - * @param $s1 Name of Target where iclamps will be placed - * @param $2 percentStart Starting scalar of the cell's threshold current when starting injection - * @param $3 ampEnd Final scalar of the cell's thrshold current when done injecting - * @param $4 delay from start of simulation to when current is injected (ms) - * @param $5 duration of stimulus (ms) - */ -proc relativeLinear() { local cellIndex, x, percentStart, percentEnd, cellThresh, del, dur localobj pts, tstim, cellList - - percentStart = $2 - percentEnd = $3 - del = $4 - dur = $5 - - pts = targetManager.getPointList($s1) - cellList = targetManager.getMETypes($s1) - //cellList = targetManager.getCells($s1) - - for cellIndex=0, pts.count()-1 { - for pts.o(cellIndex).each_point(&x) { - cellThresh = cellList.o(cellIndex).getThreshold() - tstim = new TStim(x) - - if( percentStart == percentEnd ) { - tstim.pulse( del, dur, cellThresh*(percentStart/100) ) - } else { - tstim.ramp( 0, del, cellThresh*(percentStart/100), cellThresh*(percentEnd/100), dur, 0, 0 ) - } - - stimList.append( tstim ) - } - } -} - -//----------------------------------------------------------------------------------------------- - -/*! - * Injects a train stimulus into cells wih regular frequency. - * - * @param $s1 Name of Target where stims will be placed - * @param $2 amp - * @param $3 frequency - * @param $4 width - * @param $5 offset std dev - * @param $6 del - * @param $7 dur - */ -proc train() { local cellIndex, x, offset, secIndex, foo localobj pts, tstim, jitter - - pts = targetManager.getPointList($s1) - - for cellIndex=0, pts.count()-1 { - offset = 0 - if( $5 != 0 ) { - jitter = new Random( pts.o(cellIndex).gid ) - jitter.normal(0,$5) - } - for pts.o(cellIndex).each_point(&x) { - tstim = new TStim(x) - if( $5 != 0 ) { - if($5<0){ - offset = offset-$5 - } else { - offset = jitter.repick() - } - } - tstim.train( $6+offset, $7, $2, $3, $4 ) - - stimList.append(tstim) - } - } -} - -//----------------------------------------------------------------------------------------------- - -/*! - * Injects a sinusoidal stimulus into cells wih regular frequency. - * - * @param $s1 Name of Target where stims will be placed - * @param $2 amp - * @param $3 frequency - * @param $4 offset std dev - * @param $5 del - * @param $6 dur - */ -proc sinusoid() { local cellIndex, x, offset, secIndex, foo localobj pts, tstim, jitter - pts = targetManager.getPointList($s1) - offset = 0 - for cellIndex=0, pts.count()-1 { - offset = 0 - if( $4 != 0 ) { - jitter = new Random( pts.o(cellIndex).gid ) - jitter.normal(0,$4) - } - for pts.o(cellIndex).each_point(&x) { - tstim = new TStim(x) - if( $4 != 0 ) { - if($4<0){ - offset = offset-$4 - } else { - offset = jitter.repick() - } - } - tstim.sin($2,$5+offset,$6,$3) - stimList.append(tstim) - } - } -} - -//----------------------------------------------------------------------------------------------- - -/*! - * Injects a square pulse into a cell. Note that if the starting and ending amplitude are the same, - * a simple - * - * @param $s1 Name of Target where iclamps will be placed - * @param $2 ampStart Starting amplitude of the current injected (nA) - * @param $3 ampEnd Final amplitude of the current when done injecting (nA) - * @param $4 delay from start of simulation to when current is injected (ms) - * @param $5 duration of stimulus (ms) - */ -proc linear() { local cellIndex, x, ampStart, ampEnd, del, dur localobj pts, tstim - - ampStart = $2 - ampEnd = $3 - del = $4 - dur = $5 - - //since there are no random elements, we can just evaluate the pointlist and proceed as normal - pts = targetManager.getPointList($s1) - - for cellIndex=0, pts.count()-1 { - for pts.o(cellIndex).each_point(&x) { - - tstim = new TStim(x) - - if( ampStart == ampEnd ) { - tstim.pulse( del, dur, ampStart ) - } else { - tstim.ramp( 0, del, ampStart, ampEnd, dur, 0, 0 ) - } - - stimList.append( tstim ) - } - } -} - -//----------------------------------------------------------------------------------------------- - -/*! - * Apply a poisson stimulus over n artificial synapses distributed randomly on a cell's basal and apical sections - * - * @param $s1 TargetName - * @param $2 lambda - * @param $3 weight - * @param $4 number of synapses - * @param $o5 synapseType as String Name of object to be instantiated for Synapses (nil for default of ExpSyn) - * @param $o6 synapseConfigure as String: executes code to modify synapse parameters (nil for default values) - * @param $o7 fillFunctionOverride as String: determines method used to generate the npoison spike timings (nil for default) - * @param $8 delay from start of simulation to when stimulus starts (ms) - * @param $9 duration of stimulus (ms) - */ -proc npoisson() { local cellIndex, del, dur localobj apicalPoints, basalPoints, gids - - del = $8 - dur = $9 - - //The synapses generated should only go on the apical or basal sections of any target cells - // difficulty for doing exactly same thing as old bglib: that could operate easily on the section - // of a cell, but now we work with compartment points only(?) If so, random selection won't be the same - // maybe the utility function of targetManager can return only the middle of all sections? - // e.g Target section helper - // { - // base_target apic //note that base_target must be a cell target, if section target, that is fine - // base_target dend - // } - - // open for interpretation - should I get the target, check they type and take action based on that? - // or should I send two requests to the targetManager asking for the apic & basal sections, letting - // the manager do checks and balances - gids = targetManager.gids( $s1 ) - apicalPoints = targetManager.getApicalSections( $s1 ) - basalPoints = targetManager.getBasalSections( $s1 ) - - //for each gid/TPointList - place synapses randomly and supply a spike train - for cellIndex=0, gids.size()-1 { - - //so all cells should be represented and have the same point count regardless of splitting - distributePoissonSynapses( gids.x[cellIndex], $2, $3, $4, $o5, $o6, $o7, del, dur, apicalPoints.o(cellIndex), basalPoints.o(cellIndex) ) - } - - poissonStimCount = poissonStimCount + 1 -} - -//----------------------------------------------------------------------------------------------- - -/*! - * Apply a poisson stimulus over n artificial synapses distributed randomly on a cell's basal and apical sections - * - * @param $s1 TargetName - * @param $2 weight - * @param $3 number of synapses - * @param $o4 synapseType as String Name of object to be instantiated for Synapses (nil for default of ExpSyn) - * @param $o5 synapseConfigure as String: executes code to modify synapse parameters (nil for default values) - * @param $o6 h5file - * @param $o7 tbins - * @param $o8 rate - */ -proc npoissonInhomogeneous() { local cellIndex, del, dur localobj apicalPoints, basalPoints, gids - - //TODO: refactor with npoisson - - //The synapses generated should only go on the apical or basal sections of any target cells - // difficulty for doing exactly same thing as old bglib: that could operate easily on the section - // of a cell, but now we work with compartment points only(?) If so, random selection won't be the same - // maybe the utility function of targetManager can return only the middle of all sections? - // e.g Target section helper - // { - // base_target apic //note that base_target must be a cell target, if section target, that is fine - // base_target dend - // } - - // open for interpretation - should I get the target, check they type and take action based on that? - // or should I send two requests to the targetManager asking for the apic & basal sections, letting - // the manager do checks and balances - gids = targetManager.gids( $s1 ) - apicalPoints = targetManager.getApicalSections( $s1 ) - basalPoints = targetManager.getBasalSections( $s1 ) - - //for each gid/TPointList - place synapses randomly and supply a spike train - for cellIndex=0, gids.size()-1 { - - //so all cells should be represented and have the same point count regardless of splitting - distributePoissonSynapsesInhomogeneous( gids.x[cellIndex], $2, $3, $o4, $o5, $o6, $o7, $o8, apicalPoints.o(cellIndex), basalPoints.o(cellIndex) ) - } - - poissonStimCount = poissonStimCount + 1 -} - -//----------------------------------------------------------------------------------------------- - -func getSecTypeAndIndex() { local baseOffset, secIndex localobj str, strobj - strdef subset - str = new String() -strobj = new StringFunctions() -if( strobj.substr( secname(), "soma" ) > 0 ) { - subset = "soma" - $&1 = 1 -} else if (strobj.substr( secname(), "axon") >0) { - subset = "axon" - $&1 = 2 -}else if (strobj.substr( secname(), "dend") >0) { - subset = "dend" - $&1=3 -}else if (strobj.substr( secname(), "apic") > 0) { - subset = "apic" - $&1=4 -} else { - return 1 -} -// 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 ) - $&2 = secIndex - } else { - return 1 - } -$&2=secIndex - return 0 -} - -//----------------------------------------------------------------------------------------------- - -/*! - * Given some points on a targeted cell, randomly place synapses and generate their spike trains - * - * @param $1 gid - * @param $2 lambda - * @param $3 weight - * @param $4 number of synapses - * @param $o5 synapseType as String Name of object to be instantiated for Synapses (nil for default of ExpSyn) - * @param $o6 synapseConfigure as String: executes code to modify synapse parameters (nil for default values) - * @param $o7 fillFunctionOverride as String: determines method used to generate the npoison spike timings (nil for default) - * @param $8 delay from start of simulation to when simulation starts - * @param $9 duration of stimulus - * @param $o10 TPointList with apical points - * @param $o11 TPointList with basal points - */ -proc distributePoissonSynapses() { local nCreated, nNeeded, iteration, baseSeed, sectionIndex, lambda, weight, del, dur, substituteCount, vIndex \ - localobj ran1, ran2, apicalPoints, basalPoints, activeSubset, synapseType, synapseConfigure, nc, fih, sf, fillFunctionOverride, vecstim, myvec, nil, rngInfo - strdef tstr, tstr1 - - iteration = 0 - nCreated = 0 - - rngInfo = new RNGSettings() - - baseSeed = rngInfo.getGlobalSeed() + $1 + poissonStimCount * 10000 //-> should baseSeed be an arg so that it doesn't need knowledge of # of previous npoisson stims? - lambda = $2 - weight = $3 - nNeeded = int($4) - synapseConfigure = $o6 - fillFunctionOverride = $o7 - del = $8 - dur = $9 - apicalPoints = $o10 - basalPoints = $o11 - - //determine now if the synapseType is given or we should use the default - if( object_id($o5,1) < 0 ) { - synapseType = new String("ExpSyn") - } else { - synapseType = $o5 - } - - while ( nCreated < nNeeded ) { - - //count the iterations and use that to offset the seed - iteration = iteration + 1 - - ran2 = new Random(baseSeed + iteration) - ran2.MCellRan4(baseSeed + iteration, $1+100 ) - ran2.discunif(0,1) - - ran1 = new Random(baseSeed + iteration) - ran1.MCellRan4(baseSeed + iteration, $1+110 ) - - //first choose to target apical of basal sections - if(ran2.repick()) { - activeSubset = basalPoints - } else { - activeSubset = apicalPoints - } - - if( activeSubset.count() == 0 ) { continue } - - ran1.discunif( 0, activeSubset.count()-1 ) - sectionIndex = ran1.repick() - - //do we have random access to TPointList SectionRef members? how should such a utility function look/behave? - if( !activeSubset.sclst.o(sectionIndex).exists() ) { //section not on this cpu - nCreated = nCreated + 1 //assume another split piece handles it instead - continue - } - - //the section is on this cpu, so create the synapse and fill the spike train - activeSubset.access_(sectionIndex) - sprint(tstr,"stimHelper = new %s(0.5)", synapseType.s ) - execute1(tstr, this) - - //do we need to execute some statements on the synapse? - if( object_id( synapseConfigure, 1 ) != -1 ) { - sf = new StringFunctions() - - //need to replace all '%s' with stimHelper object - substituteCount = 0 - tstr = synapseConfigure.s - - // keep checking the string for '%s'; as long as one is there, rebuild the string around it - while( sf.substr( tstr, "%s" ) != -1 ) { - sf.head( tstr, "%s", tstr1 ) - sf.tail( tstr, "%s", tstr ) - sprint( tstr, "%s%s%s", tstr1, stimHelper, tstr ) - } - - //surround with curly braces then execute - sprint(tstr, "{%s}", tstr ) - //print "execute ", tstr - execute1(tstr,this) - - //debug - only works for ExpSyn - //print stimHelper.e, " ", stimHelper.tau - } - - if( object_id( fillFunctionOverride, 1 ) == -1 ) { //no override - so use default spike timing fill method - - nc = new NetCon( nil, stimHelper,10,1,weight ) //has nil source - - // register FInitializeHandler for our FillSpikeTrain - sprint(tstr,"FillSpikeTrain(%s,%f,%f,%f,%d)", nc, lambda, del, dur, baseSeed+iteration ) - synapseEventGenerator.append( new String(tstr) ) - fih = new FInitializeHandler(tstr, this) - //print tstr - } else { - tstr = "print \"fillFunctionOverride not implemented.\"" - } - - //keep variables so they don't get deleted - poissonNetconList.append(nc) - poissonSynList.append(stimHelper) - handleList.append(fih) - - nCreated = nCreated + 1 - } -} - -//----------------------------------------------------------------------------------------------- - -/*! - * Given some points on a targeted cell, randomly place synapses and generate their spike trains - * - * @param $1 gid - * @param $2 weight - * @param $3 number of synapses - * @param $o4 synapseType as String Name of object to be instantiated for Synapses (nil for default of ExpSyn) - * @param $o5 synapseConfigure as String: executes code to modify synapse parameters (nil for default values) - * @param $o6 h5file - * @param $o7 tbins - * @param $o8 rate - * @param $o9 apical points - * @param $o10 basal points - */ -proc distributePoissonSynapsesInhomogeneous() { local nCreated, nNeeded, iteration, baseSeed, sectionIndex, weight, substituteCount, tIndex \ - localobj ran1, ran2, apicalPoints, basalPoints, activeSubset, synapseType, synapseConfigure, nc, sf, h5file, tbins, rate, h5reader, vecstim, tbins_vec, rate_vec, inhpoiss, exprng, uniformrng, nil, rngInfo - strdef tstr, tstr1 - - //TODO: refactor this code with distributePoissonSynapses - - iteration = 0 - nCreated = 0 - - rngInfo = new RNGSettings() - - baseSeed = rngInfo.getGlobalSeed() + $1 + poissonStimCount * 10000 //-> should baseSeed be an arg so that it doesn't need knowledge of # of previous npoisson stims? - weight = $2 - nNeeded = int($3) - synapseConfigure = $o5 - h5file = $o6 - tbins = $o7 - rate = $o8 - apicalPoints = $o9 - basalPoints = $o10 - - //determine now if the synapseType is given or we should use the default - if( object_id($o4,1) < 0 ) { - synapseType = new String("ExpSyn") - } else { - synapseType = $o4 - } - - - //open h5file ... only need to read it once for all synapses - // they all share the same tbins_vec and rate_vec object - h5reader = new HDF5Reader( h5file.s ) - - //get tbins from file - put in vector for InhPoissonStim - h5reader.loadData( tbins.s ) - nrows = h5reader.numberofrows( tbins.s ) - - tbins_vec = new Vector() - for tIndex=0,nrows-1 { - //print h5reader.getData( tbins.s, tIndex, 0 ) - tbins_vec.append( h5reader.getData( tbins.s, tIndex, 0 ) ) - } - - //get rate from file - put in vector for InhPoissonStim - h5reader.loadData( rate.s ) - nrows = h5reader.numberofrows( rate.s ) - - rate_vec = new Vector() - for tIndex=0,nrows-1 { - //print h5reader.getData( rate.s, tIndex, 0 ) - rate_vec.append( h5reader.getData( rate.s, tIndex, 0 ) ) - } - - vecStimList.append( tbins_vec ) - vecStimList.append( rate_vec ) - - - while ( nCreated < nNeeded ) { - - //count the iterations and use that to offset the seed - iteration = iteration + 1 - - ran2 = new Random(baseSeed + iteration) - ran2.MCellRan4(baseSeed + iteration, $1+100 ) - ran2.discunif(0,1) - - ran1 = new Random(baseSeed + iteration) - ran1.MCellRan4(baseSeed + iteration, $1+110 ) - - //first choose to target apical of basal sections - if(ran2.repick()) { - activeSubset = basalPoints - } else { - activeSubset = apicalPoints - } - - if( activeSubset.count() == 0 ) { continue } - - ran1.discunif( 0, activeSubset.count()-1 ) - sectionIndex = ran1.repick() - - //do we have random access to TPointList SectionRef members? how should such a utility function look/behave? - if( !activeSubset.sclst.o(sectionIndex).exists() ) { //section not on this cpu - nCreated = nCreated + 1 //assume another split piece handles it instead - continue - } - - //the section is on this cpu, so create the synapse and fill the spike train - activeSubset.access_(sectionIndex) - sprint(tstr,"stimHelper = new %s(0.5)", synapseType.s ) - execute1(tstr, this) - - //do we need to execute some statements on the synapse? - if( object_id( synapseConfigure, 1 ) != -1 ) { - sf = new StringFunctions() - - //need to replace all '%s' with stimHelper object - substituteCount = 0 - tstr = synapseConfigure.s - - // keep checking the string for '%s'; as long as one is there, rebuild the string around it - while( sf.substr( tstr, "%s" ) != -1 ) { - sf.head( tstr, "%s", tstr1 ) - sf.tail( tstr, "%s", tstr ) - sprint( tstr, "%s%s%s", tstr1, stimHelper, tstr ) - } - - //surround with curly braces then execute - sprint(tstr, "{%s}", tstr ) - //print "execute ", tstr - execute1(tstr,this) - - //debug - only works for ExpSyn - //print stimHelper.e, " ", stimHelper.tau - } - - - //ignore rate dataset for now - //h5reader.loadData( rate.s ) - //nrows = h5reader.numberofrows( rate.s ) - - inhpoiss = new InhPoissonStim(0.5) - nc = new NetCon(inhpoiss, stimHelper) - nc.weight[0] = weight - nc.delay = 0.1 - inhpoiss.setTbins( tbins_vec ) - inhpoiss.setRate( rate_vec ) - - // The InhPoissonStim needs random #s - exprng = new Random() - exprng.ACG(baseSeed+iteration) - exprng.negexp(1) - uniformrng = new Random() - uniformrng.ACG(baseSeed+iteration) - uniformrng.uniform(0.0,1.0) - inhpoiss.setRNGs(exprng, uniformrng) - - //keep variables so they don't get deleted - poissonNetconList.append(nc) - poissonSynList.append(stimHelper) - vecStimList.append( inhpoiss ) - vecStimList.append( exprng ) - vecStimList.append( uniformrng ) - - nCreated = nCreated + 1 - } -} - -//----------------------------------------------------------------------------------------------- - -/*! - * This function creates spike events distributed as the intervals between events in a poisson distribution. - * It populates a vector with randomly generated time values and assigns those to the given netcon. - * - * @param $o1 Netcon Object - * @param $2 mean firing rate. - * @param $3 Delay - * @param $4 TimeDuration - * @param $5 seed - */ -proc FillSpikeTrain() { local i,tmpt, time, index, SpikeNeeded, nSpikeNeeded, SpikeDiff localobj ran, tVec, ran2 - //print "Seed = ", $5 - tVec = new Vector() - ran = new Random() - ran.ACG($5) - //ran.MCellRan4() - ran.negexp(1000/$2) - //ran.negexp($2) - time = 0 - SpikeNeeded = (($4 * $2)/1000) - nSpikeNeeded = int (SpikeNeeded) - if((SpikeNeeded - nSpikeNeeded) >= 0.5){ - nSpikeNeeded = nSpikeNeeded + 1 - } - - while (time < $4){ - tmpt = ran.repick() - //time = time + ran.repick() - if(tmpt+time > $4){ - break - } - time = time + tmpt - tVec.append(time) - } - if(SpikeNeeded < 1) { - time = 0 - } - - for i =0, tVec.size-1 { - $o1.event(tVec.x[i] + $3) - } -} - -//----------------------------------------------------------------------------------------------- - -/*! - * Take voltage values from a file source and play to an SEClamp on a cell's soma. This assumes that the files are - * HDF5 reports produced by previous runs of the simulator using the fixed time step method. TODO: support other file - * layouts. Difficult to say if the Bin report bbp format will be among those since it is laid out more appropriately - * for accessing a single timestep of the circuit at once rather than a single cell's states all at one - * - * @param $s1 Target name to receive the SEClamp - * @param $s2 path where data files can be found - * @param $s3 label of dataset to use from h5 files - * @param $4 Stim delay - * @param $5 Stim duration - */ -proc replayVoltageTrace() { local gidIndex, x, dtReport, del, dur \ - localobj gids, pts, dataset, filename, voltageReader, vvec, dvec, tvec, stim, pc - - del = $4 - dur = $5 - - dataset = new String() - filename = new String() - gids = targetManager.gids( $s1 ) - pts = targetManager.getPointList( $s1 ) - - for gidIndex=0, gids.size()-1 { - - gid = gids.x[gidIndex] - sprint( filename.s, "%s/a%d.h5", $s2, gid ) - voltageReader = new HDF5Reader( filename.s ) - - sprint( dataset.s, "/a%d/%s/data", gid, $s3 ) - voltageReader.loadData( dataset.s ) - - dtReport = voltageReader.getAttributeValue( dataset.s, "Dt" ) - - vvec = new Vector() - voltageReader.getColumnData( dataset.s, 0, vvec ) - - // we need to give a tvec to SEClamp, use the stimulation timing values to generate - tvec = new Vector() - tvec = tvec.indgen( del, del+dur, dtReport ) - - // if we have more voltage data than time or the reverse, balance the vectors - if( vvec.size() > tvec.size() ) { - vvec.resize( tvec.size() ) - } else if ( tvec.size() < vvec.size() ) { - tvec.resize( vvec.size() ) - } - - for pts.o(gidIndex).each_point(&x) { - stim = new SEClamp(x) - stim.rs = 0.01 - stim.dur1 = dur - - dvec = vvec.c - dvec.play_remove() - dvec.play( &stim.amp1, tvec ) - voltageReplayList.append( dvec ) - - stimList.append(stim) - } - - timeReplayList.append( tvec ) - - } - -} - -//----------------------------------------------------------------------------------------------- - -/** - * For save state, we need to indicate that point processes made by this manager object will be ignored during the save action. - * - * @param $o1 savestate object - */ -proc saveStatePreparation() { local ppIndex - for ppIndex=0, poissonSynList.count()-1 { - //print poissonSynList.o(ppIndex) - $o1.ignore( poissonSynList.o(ppIndex) ) - } - - for ppIndex=0, stimList.count()-1 { - $o1.ignore( stimList.o(ppIndex).ic ) - } -} - -proc reevent() { local eIndex, vIndex, res - for eIndex=0, synapseEventGenerator.count()-1 { - //print synapseEventGenerator.o(eIndex).s - res = execute1( synapseEventGenerator.o(eIndex).s, this ) - } - - //while I am here, let me check all other stims - /* - print "replay as needed ", stimList.count() - for eIndex=0, stimList.count()-1 { - //print "stim info ", stimList.o(eIndex).ic.amp, " ", stimList.o(eIndex).ic.del, " ", stimList.o(eIndex).ic.dur - //print stimList.o(eIndex).tvec.x[0], " ", stimList.o(eIndex).stim.size(), " ", stimList.o(eIndex).stim.x[0] - for vIndex=0, stimList.o(eIndex).tvec.size()-1 { - print "time ", stimList.o(eIndex).tvec.x[vIndex] - } - for vIndex=0, stimList.o(eIndex).stim.size()-1 { - print "amp ", stimList.o(eIndex).stim.x[vIndex] - } - stimList.o(eIndex).stim.play(&stimList.o(eIndex).ic.amp,stimList.o(eIndex).tvec,1) - }*/ -} - -//----------------------------------------------------------------------------------------------- - -/* - * Places an SEClamp on the targeted locations which will hold the voltage to the indicated value - * - * @param $s1 Target name - * @param $2 target membrane voltage Vm - * @param $3 delay for stimulus onset (not used) - * @param $4 duration of stimulus - */ -proc seclamp() { local del, dur, gidIndex, x localobj pts, stim - del = $3 - dur = $4 - - //gids = targetManager.gids( $s1 ) - pts = targetManager.getPointList( $s1 ) - - for gidIndex=0, pts.count()-1 { - for pts.o(gidIndex).each_point(&x) { - stim = new SEClamp(x) - stim.rs = 0.01 - stim.dur1 = dur - stim.amp1 = $2 - - stimList.append(stim) - } - } - -} - -endtemplate StimulusManager diff --git a/core/hoc/TDistFunc.hoc b/core/hoc/TDistFunc.hoc deleted file mode 100644 index d4813552..00000000 --- a/core/hoc/TDistFunc.hoc +++ /dev/null @@ -1,82 +0,0 @@ -/** - * @file TDistFunc.hoc - * @brief - * @author rajnish - * @date 2010-11-23 - * @remark Copyright © BBP/EPFL 2005-2011; All rights reserved. Do not distribute without further notice. - */ - -D__TDISTFUNCHOC = 1 - -begintemplate TDistFunc - strdef stmp,distfunc,mech - public distribute, ndistribute, distribute_section_variable - - proc distribute(){local x localobj sl - sl = $o1 - mech = $s2 - distfunc = $s3 - sprint(distfunc,"%%s %s(%%f) = %s",mech, distfunc) - //print distfunc - forsec sl for(x,0) { - sprint (stmp,distfunc,secname(),x,distance(x)) - //print stmp - execute(stmp) - } - } - - proc distribute_section_variable(){ localobj sl - sl = $o1 - mech = $s2 - distfunc = $s3 - sprint(distfunc,"%%s %s = %s",mech, distfunc) - //print distfunc - forsec sl { - sprint (stmp,distfunc,secname(),distance(0.5)) - //print stmp - execute(stmp) - } - } - - proc ndistribute(){local x,xroot,norm_factor localobj sl - //print "Inside distribute .. " - sl = $o1 - mech = $s2 - distfunc = $s3 - norm_factor = 1 - if(numarg()>3) norm_factor = getLongestBranch(sl) - sprint(distfunc,"%%s %s(%%f) = %s",mech, distfunc) - //print distfunc - forsec sl for(x,0) { - sprint (stmp,distfunc,secname(),x,distance(x)/norm_factor) - // print stmp - execute(stmp) - } - } - - - func getLongestBranch(){local maxL,d localobj distallist,sref - forsec $o1 { - distance() - break - } - maxL = 0 - d = 0 - distallist = new SectionList() - forsec $o1 { - sref = new SectionRef() - if (sref.nchild==0) distallist.append() - } - forsec distallist{ - d = distance(1) - if(maxL1) rand = $o2 - //the most important argument is continuos. this essentily - // allowed me to draw the input current with points. - stim.play(&ic.amp,tvec,1) - } - -/** -* stim and tvec can not be destroyed for any reason or the pointer to the -* ic.amp (play function) will be lost. Therefore for each new stimulus -* decleration, vecs are initialized but not deleted. -*/ - proc initinject(){ - tvec.resize(0) - stim.resize(0) - } - - //basically used to start and end a stimulus - proc addpoint(/*amp,tstart*/){tvec.append($2) stim.append($1)} - - proc appendpulse(/*amp,dur*/){tvec.append(0,$2) stim.append($1,$1)} - - - proc pulse(/*start,dur,max_amp[,min_amp]*/){local min_amp - min_amp = 0 - initinject() - if(numarg()>3) min_amp = $4 - addpoint(min_amp,$1) - appendpulse($3,$2) - addpoint(min_amp,0) - addpoint(min_amp,0) - tvec.integral //make the set of durations a timeline - } - - proc ramp(/*start_amp,dur,amp1,amp2,dur,end_amp,dur*/){local min_amp - initinject() - appendpulse($1,$2) - tvec.append(0,$5) stim.append($3,$4) //the ramp - appendpulse($6,$7) - tvec.integral - } - - /** - * Stimulus with repeated pulse injections at a specified frquency - * - * @param $1 del - time when stimulation starts - * @param $2 dur - time when stimulation ends - * @param $3 amp - the amplitude of a each pulse - * @param $4 frequency - determines the number of pulses per second (hz) - * @param $5 width - the duration of a single pulse - */ - proc train() { local tpulse, remtime, width - - tpulse = 1000/$4 - width = $5 - - if( width > tpulse ) { //just use a standard pulse - pulse( $1, $2, $3 ) - return - } - - initinject() - //addpoint( 0, $1 ) - remtime = $2 - nextpulse = $1 - -//TODO: for better clarity, don't try to shift pulse onsets by one dt (only done to match old bglib) - - while( remtime > 0 ) { - if( width < remtime ) { //can fit whole pulse in remaining time - addpoint( 0, nextpulse ) - addpoint( $3, nextpulse ) - addpoint( $3, nextpulse+width ) - addpoint( 0, nextpulse+width ) - } else { //final pulse, just paritial - addpoint( 0, nextpulse ) - addpoint( $3, nextpulse ) - addpoint( $3, $1+$2 ) - addpoint( 0, $1+$2 ) - } - - nextpulse = nextpulse + tpulse - remtime = remtime - tpulse - } - - addpoint( 0, $1+$2 ) - } - - proc sin(/*amp,start,dur,freq*/){ - initinject() - tvec.indgen($2,$2+$3,0.025) - stim.resize(tvec.size()) - stim.sin($4,0,0.025) - stim.mul($1) - tvec.append(tvec.x[tvec.size()-1]+0.025) - stim.append(0) - } - - func zap(){return sin(exp(exp($1/5000)))}//original was in seconds (fron PulseQ) - - proc sinspec(/*start,dur*/){ - initinject() - tvec.indgen(0,10000,0.5) - stim.copy(tvec)//must preserve the original object - stim.apply("zap") - } - - proc segment(/*amp1,dur1,amp2,dur2,...*/){local i,k - initinject() - for(i=1;i<=numarg();i+=1) { - k = $i i+=1 appendpulse(k,$i) - } - addpoint(0,0)//putting it to zero after the last injection - tvec.integral - } - - /** - * Note that the easiest way to write the noise function is: - * rand.play(&iclamp.amp) - * However, this does not allow you the control of when to start and to end the - * stimulus, so there is a memory overhead of preparing the noise vector in advanced - * and put it in memory. I do now know of any differences in performance between - * the vector.play and the rand.play functions. - */ - proc noise(/*start,dur,mean,var[,dt]*/){local dT, tstim_end, prev_t, next_t, tstep localobj rngInfo - if(numarg()>4) dT = $5 else dT = 0.5 - initinject() - rand.normal($3,$4) - tstim_end = $1+$2 - - rngInfo = new RNGSettings() - - // To help with save state, we want Noise stimulus to continue where it left off - // so first draw the numbers that were already used (not very efficient, but can be improved later) - if( rngInfo.getRNGMode() != rngInfo.COMPATIBILITY ) if( $1 > 0 ) { - - // this can use up a lot of memory as simulation times increase. Should do a proper rewrite, but for now - // make sure no chunk is larger than what we will use for the actual simulation - prev_t = 0 - tstep = tstim_end-$1-dt - - while( prev_t < $1-dt ) { - if( prev_t + tstep < $1-dt ) { - next_t = prev_t + tstep - } else { - next_t = $1-dt - } - - tvec.indgen(prev_t,next_t,dT) - stim.copy(tvec.c) - stim.setrand(rand) - - prev_t = next_t+dt - } - } - - tvec.indgen($1,tstim_end,dT) - stim.copy(tvec.c) - stim.setrand(rand) - - // If the stimulus starts with a delay from simulation start time, we will want 0 amps until the next Dt. - // Note that for Resume, we will be at t > 0 and delay was offset by previous save time. Ergo, we must check that delay > t - // Note that for legacy versions, we always start with 0 amps; when we stop support for legacy, we should selectively set initial slot to 0 - if( rngInfo.getRNGMode() == rngInfo.COMPATIBILITY || $1 > t ) { - stim.x[0] = 0 - } - - // If the stimulus should end before the simulation, make sure the final value will set the amp to 0 - if(tstim_end < tstop) { - stim.x[stim.size()-1] = 0 - } - - addpoint(0,0) - } - - proc delay(/*dur*/){ - tvec.add($1) - } - - -/** -* Meant to give a copy option (like vector.c) but it is very different if the concept. -* Note that the stim and tvec vectors of the newly created TStim are pointing to the -* old one. So when changing stimulation (pulse, ramp, noise ets') it will be altered -* in all the objects created by c. -*/ - obfunc c(/*[x]*/){localobj ts - if(rand!=null) { - ts = new TStim(x,rand) - }else{ - ts = new TStim(x,rand) - } - if(numarg()>0) ts.loc($1) - ts.stim = stim //set stim of the new object to reference the old one - ts.tvec = tvec //set stim of the new object to reference the old one - return ts - } - - -/** -* It is better to thing about TStim class as a poitprocess object. So I exposed some -* of the functionallity of pointprocess object here. More is always wellcome. -*/ - proc loc(){ ic.loc(x=$1) } - func get_loc(){ return get_loc() }//don't forget the pop_section() - proc mark(/*Shape,color_index*/){$o1.point_mark(ic,$2)} - - -endtemplate TStim diff --git a/core/hoc/neurodamus.hoc b/core/hoc/neurodamus.hoc index 92a83527..a6e9f2dc 100644 --- a/core/hoc/neurodamus.hoc +++ b/core/hoc/neurodamus.hoc @@ -9,9 +9,8 @@ // Neurodamus - Do not change order {load_file("RNGSettings.hoc")} -{load_file("StimulusManager.hoc")} {load_file("Cell.hoc")} -{load_file("TDistFunc.hoc")} {load_file("Report.hoc")} {load_file("ConnectionUtils.hoc")} {load_file("MorphIO.hoc")} +{load_file("Map.hoc")} diff --git a/neurodamus/io/sonata_config.py b/neurodamus/io/sonata_config.py index fc20274f..b87dd1dd 100644 --- a/neurodamus/io/sonata_config.py +++ b/neurodamus/io/sonata_config.py @@ -364,13 +364,17 @@ def parsedStimuli(self): "extracellular_stimulation": "Extracellular", "conductance": "Conductance" } + _module_translation = { + "seclamp": "SEClamp", + "subthreshold": "SubThreshold" + } stimuli = {} for name in self._sim_conf.list_input_names: stimulus = self._translate_dict("inputs", self._sim_conf.input(name)) self._adapt_libsonata_fields(stimulus) - stimulus["Pattern"] = "SEClamp" if stimulus["Pattern"] == "seclamp" \ - else snake_to_camel(stimulus["Pattern"]) + stimulus["Pattern"] = _module_translation.get(stimulus["Pattern"], + snake_to_camel(stimulus["Pattern"])) stimulus["Mode"] = _input_type_translation.get(stimulus["Mode"], stimulus["Mode"]) stimuli[name] = stimulus return stimuli diff --git a/neurodamus/node.py b/neurodamus/node.py index ea8260aa..0e064b64 100644 --- a/neurodamus/node.py +++ b/neurodamus/node.py @@ -643,7 +643,7 @@ def enable_stimulus(self): log_stage("Stimulus Apply.") # for each stimulus defined in the config file, request the stimmanager to instantiate - self._stim_manager = StimulusManager(self._target_manager, None) + self._stim_manager = StimulusManager(self._target_manager) # build a dictionary of stims for faster lookup : useful when applying 10k+ stims # while we are at it, check if any stims are using extracellular @@ -1151,7 +1151,6 @@ def _sim_init_neuron(self): with timeit(name="restoretime"): logging.info("Restoring state...") self._stim_manager.saveStatePreparation(self._bbss) - self._stim_manager.reevent() self._bbss.vector_play_init() self._restart_events() # On restore the event queue is cleared return # Upon restore sim is ready diff --git a/neurodamus/stimulus_manager.py b/neurodamus/stimulus_manager.py index b89c3cc1..169e7a86 100644 --- a/neurodamus/stimulus_manager.py +++ b/neurodamus/stimulus_manager.py @@ -21,7 +21,7 @@ import logging from .core import NeurodamusCore as Nd from .utils.logging import log_verbose -from .core.configuration import SimConfig +from .core.configuration import SimConfig, ConfigurationError from .core.stimuli import CurrentSource, ConductanceSource from .core import random @@ -35,10 +35,7 @@ class StimulusManager: _stim_types = {} # stimulus handled in Python - def __init__(self, target_manager, elec_manager=None): - base_seed = SimConfig.run_conf.get("BaseSeed") - xargs = () if base_seed is None else (base_seed,) - self._hoc = Nd.StimulusManager(target_manager, elec_manager, *xargs) + def __init__(self, target_manager): self._target_manager = target_manager self._stim_seed = SimConfig.run_conf.get("StimulusSeed") self._stimulus = [] @@ -46,30 +43,21 @@ def __init__(self, target_manager, elec_manager=None): def interpret(self, target_spec, stim_info): stim_t = self._stim_types.get(stim_info["Pattern"]) + if not stim_t: + raise ConfigurationError("No implementation for Stimulus %s " % stim_info["Pattern"]) if self._stim_seed is None and getattr(stim_t, 'IsNoise', False): logging.warning("StimulusSeed unset (default %d), " "set explicitly to vary noisy stimuli across runs", SimConfig.rng_info.getStimulusSeed()) - # Get either hoc target or sonata node_set, needed for python and hoc interpret - # If sonata node_set, internally register the target and add to hoc TargetList target = self._target_manager.get_target(target_spec) - 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._cell_manager - stim = stim_t(target, stim_info, cell_manager) - self._stimulus.append(stim) - else: - # Fallback to hoc stim manager - self._hoc.interpret(target_spec.name, stim_info.hoc_map) + log_verbose("Interpret stimulus") + cell_manager = self._target_manager._cell_manager + stim = stim_t(target, stim_info, cell_manager) + self._stimulus.append(stim) def interpret_extracellulars(self, injects, stimuli): - """Hoc only implementation for extra-cellulars""" - self._hoc.interpretExtracellulars(injects.hoc_map, stimuli.hoc_map) - - def __getattr__(self, item): - logging.debug("Pass unknown method request to Hoc") - return getattr(self._hoc, item) + """Pending for BBPBGLIB-890""" + raise ConfigurationError("input_type extracellular_stimulation is not implemented") def reset_helpers(self): ShotNoise.stimCount = 0 @@ -82,13 +70,16 @@ def register_type(cls, stim_class): cls._stim_types[stim_class.__name__] = stim_class return stim_class + def saveStatePreparation(self, ss_obj): + for stim in self._stimulus: + ss_obj.ignore(stim) + class BaseStim: """ Barebones stimulus class """ - IsPythonOnly = False IsNoise = False def __init__(self, _target, stim_info: dict, _cell_manager): @@ -101,7 +92,6 @@ class OrnsteinUhlenbeck(BaseStim): """ Ornstein-Uhlenbeck process, injected as current or conductance """ - IsPythonOnly = True IsNoise = True stimCount = 0 # global count for seeding @@ -193,7 +183,6 @@ class RelativeOrnsteinUhlenbeck(OrnsteinUhlenbeck): Ornstein-Uhlenbeck process, injected as current or conductance, relative to cell threshold current or inverse input resistance """ - IsPythonOnly = True IsNoise = True def __init__(self, target, stim_info: dict, cell_manager): @@ -231,7 +220,6 @@ class ShotNoise(BaseStim): ShotNoise stimulus handler implementing Poisson shot noise with bi-exponential response and gamma-distributed amplitudes """ - IsPythonOnly = True IsNoise = True stimCount = 0 # global count for seeding @@ -361,7 +349,6 @@ class RelativeShotNoise(ShotNoise): RelativeShotNoise stimulus handler, same as ShotNoise but parameters relative to cell threshold current or inverse input resistance """ - IsPythonOnly = True IsNoise = True def __init__(self, target, stim_info: dict, cell_manager): @@ -411,7 +398,6 @@ class AbsoluteShotNoise(ShotNoise): AbsoluteShotNoise stimulus handler, same as ShotNoise but parameters from given mean and std. dev. """ - IsPythonOnly = True IsNoise = True def __init__(self, target, stim_info: dict, cell_manager): @@ -568,17 +554,12 @@ def __init__(self, target, stim_info: dict, cell_manager): rng_mode = SimConfig.rng_info.getRNGMode() # simulation RNGMode # setup RNG - if rng_mode == SimConfig.rng_info.COMPATIBILITY: - rand = lambda gid: random.RNG(seed=gid + Noise.stimCount) - elif rng_mode == SimConfig.rng_info.UPMCELLRAN4: - rand = lambda gid: random.MCellRan4(Noise.stimCount * 10000 + 100, - SimConfig.rng_info.getGlobalSeed() + - SimConfig.rng_info.getStimulusSeed() + - gid * 1000) - elif rng_mode == SimConfig.rng_info.RANDOM123: + if rng_mode == SimConfig.rng_info.RANDOM123: rand = lambda gid: random.Random123(Noise.stimCount + 100, SimConfig.rng_info.getStimulusSeed() + 500, gid + 300) + else: + raise Exception("rng_mod is not RANDOM123") # apply stim to each point in target tpoints = target.getPointList(cell_manager) @@ -590,7 +571,7 @@ def __init__(self, target, stim_info: dict, cell_manager): rng = rand(gid) # setup RNG # draw already used numbers - if rng_mode != SimConfig.rng_info.COMPATIBILITY and self.delay > 0: + if self.delay > 0: self.draw_already_used_numbers(rng, sim_dt) for sec_id, sc in enumerate(tpoint_list.sclst): diff --git a/tests/integration/test_sonata_config.py b/tests/integration/test_sonata_config.py index a0793245..6243cabd 100644 --- a/tests/integration/test_sonata_config.py +++ b/tests/integration/test_sonata_config.py @@ -130,6 +130,14 @@ def test_SimConfig_from_sonata(): "module": "hyperpolarizing", "delay": 0.0, "duration": 10000.0 + }, + "subthreshould_mosaic": { + "module": "subthreshold", + "input_type": "current_clamp", + "delay": 0.0, + "duration": 30000.0, + "node_set": "Mosaic", + "percent_less": 50.0 } } } @@ -210,3 +218,6 @@ def test_parse_inputs(sonataconfig): assert input_RSN["SDPercent"] == 40. assert input_RSN["Dt"] == 0.25 assert input_RSN.get("Seed") is None + input_subthreshold = SimConfig.stimuli["subthreshould_mosaic"] + assert input_subthreshold["Pattern"] == "SubThreshold" + assert input_subthreshold["PercentLess"] == 50.0