Changeset 141

Show
Ignore:
Timestamp:
10/05/07 11:09:44 (1 year ago)
Author:
apdavison
Message:

* Added new standard cell: AdaptiveExponentialIF_alpha - only available in nest2 for now.
* Refactored recording and printing in nest2 to reduce/eliminate duplicated code. This has not yet been well tested, and may be broken.
* Fixed several of the standard cells in nest2 so the IF_curr_alpha* and IF_curr_exp tests now pass.
* Fixed the connect() function in nest2 so it returns something useful rather than just None. It now returns a Connection object, which has a (partially implemented) weight property, and will have a delay property. Should consider making this object available in all modules, since the low-level API currently lacks a way of accessing/changing weights/delays of individual connections.

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • trunk/common.py

    r137 r141  
    1919 
    2020dt = 0.1 
    21  
    2221 
    2322# ============================================================================== 
    …… 
    9897                return self._position 
    9998            except (AttributeError, KeyError): 
    100                 self._position = (int(self), 0, 0) 
     99                self._position = (float(self), 0.0, 0.0) 
    101100                return self._position 
    102101 
    …… 
    283282        'tau_syn_I' : 2.0, 
    284283        'i_offset'  : 0.0, 
    285     } # v_init??? 
     284        'v_init'    : -65.0, 
     285    } 
     286 
     287class AdaptiveExponentialIF_alpha(StandardCellType): 
     288    """adaptive exponential integrate and fire neuron according to  
     289    Brette R and Gerstner W (2005) Adaptive Exponential Integrate-and-Fire Model as 
     290            an Effective Description of Neuronal Activity. J Neurophysiol 94:3637-3642 
     291    """ 
     292     
     293    default_parameters = { 
     294        'v_init'    : -70.6, #'V_m'        # Initial membrane potential in mV 
     295        'w_init'    : 0.0,   #'w'          # Spike-adaptation current in nA (nest:pA) 
     296        'cm'        : 0.281, #'C_m'        # Capacity of the membrane in nF 
     297        'tau_refrac': 0.0,   #'t_ref'      # Duration of refractory period in ms. 
     298        'v_spike'   : 0.0,   #'V_peak'     # Spike detection threshold in mV. 
     299        'v_reset'   : -70.6, #'V_reset'    # Reset value for V_m after a spike. In mV. 
     300        'v_rest'    : -70.6, #'E_L'        # Resting membrane potential (Leak reversal potential) in mV. 
     301        'tau_m'     : 9.3667,  #'g_L'        # Membrane time constant in ms (nest:Leak conductance in nS.) 
     302        'i_offset'  : 0.0,   #'I_e'        # Offset current in nA (nest:pA) 
     303        'a'         : 4.0,                # Subthreshold adaptation conductance in nS. 
     304        'b'         : 0.0805,           # Spike-triggered adaptation in nA (nest:pA). 
     305        'delta_T'   : 2.0,   # Delta_T    # Slope factor in mV 
     306        'tau_w'     : 144.0, #'tau_w'      # Adaptation time constant in ms 
     307        'v_thresh'  : -50.4,  #'V_t'        # Spike initiation threshold in mV (V_th can also be used for compatibility). 
     308        'e_rev_E'   : 0.0,    #'E_ex'       # Excitatory reversal potential in mV. 
     309        'tau_syn_E' : 5.0,  #'tau_ex'     # Rise time of excitatory synaptic conductance in ms (alpha function). 
     310        'e_rev_I'   : -80.0,  #'E_in'       # Inhibitory reversal potential in mV. 
     311        'tau_syn_I' : 5.0, #'tau_in'     # Rise time of the inhibitory synaptic conductance in ms (alpha function). 
     312    } 
    286313 
    287314class SpikeSourcePoisson(StandardCellType): 

  • trunk/nest2.py

    r140 r141  
    22""" 
    33PyNEST implementation of the PyNN API. 
    4 $Id:nest.py 5 2007-04-16 15:01:24Z davison $ 
     4$Id$ 
    55""" 
    66__version__ = "$Revision:5 $" 
    …… 
    1515from math import * 
    1616 
    17 ll_spike_files = [] 
    18 ll_v_files     = [] 
    19 hl_spike_files = {} 
    20 hl_v_files     = {} 
    21 hl_c_files     = {} 
     17#ll_spike_files = [] 
     18#ll_v_files     = {} 
     19#hl_spike_files = {} 
     20#hl_v_files     = {} 
     21#hl_c_files     = {} 
     22recorder_dict = {} 
    2223tempdirs       = [] 
    2324dt             = 0.1 
    24  
    25  
     25recording_device_names = {'spikes': 'spike_detector', 
     26                          'v': 'voltmeter', 
     27                          'conductance': 'conductancemeter'} 
    2628 
    2729# ============================================================================== 
    …… 
    6668             
    6769 
     70class Connection(object): 
     71     
     72    def __init__(self, pre, post): 
     73        self.pre = pre 
     74        self.post = post 
     75        conn_dict = nest.GetConnections([pre], 'static_synapse')[0] 
     76        if conn_dict: 
     77            self.port = len(conn_dict['targets'])-1 
     78        else: 
     79            raise Exception("Could not get port number for connection between %s and %s" % (pre, post)) 
     80 
     81    def _set_weight(self, w): 
     82        pass 
     83 
     84    def _get_weight(self): 
     85        # this needs to be modified to take account of threads 
     86        # also see nest.GetConnection (was nest.GetSynapseStatus) 
     87        conn_dict = nest.GetConnections([self.pre],'static_synapse')[0] 
     88        if conn_dict: 
     89            return conn_dict['weights'][self.port] 
     90        else: 
     91            return None 
     92         
     93    weight = property(_get_weight, _set_weight) 
     94 
    6895# ============================================================================== 
    6996#   Standard cells 
    …… 
    75102     
    76103    translations = { 
    77             'v_rest'    : ('U0'    , "parameters['v_rest']"), 
    78             'v_reset'   : ('Vreset', "parameters['v_reset']"), 
    79             'cm'        : ('C'     ,  "parameters['cm']*1000.0"), # C is in pF, cm in nF 
    80             'tau_m'     : ('Tau'   , "parameters['tau_m']"), 
    81             'tau_refrac': ('TauR'  , "max(dt,parameters['tau_refrac'])"), 
    82             'tau_syn_E' : ('TauSynE', "parameters['tau_syn_E']"), 
    83             'tau_syn_I' : ('TauSynI', "parameters['tau_syn_I']"), 
    84             'v_thresh'  : ('Theta' , "parameters['v_thresh']"), 
    85             'i_offset'  : ('I0'    ,  "parameters['i_offset']*1000.0"), # I0 is in pA, i_offset in nA 
    86             'v_init'    : ('u'     , "parameters['v_init']"), 
     104            'v_rest'    : ('E_L'    , "parameters['v_rest']"), 
     105            'v_reset'   : ('V_reset', "parameters['v_reset']"), 
     106            'cm'        : ('C_m'    , "parameters['cm']*1000.0"), # C_m is in pF, cm in nF 
     107            'tau_m'     : ('tau_m'  , "parameters['tau_m']"), 
     108            'tau_refrac': ('tau_ref', "max(dt,parameters['tau_refrac'])"), 
     109            'tau_syn_E' : ('tau_ex' , "parameters['tau_syn_E']"), 
     110            'tau_syn_I' : ('tau_in' , "parameters['tau_syn_I']"), 
     111            'v_thresh'  : ('V_th'   , "parameters['v_thresh']"), 
     112            'i_offset'  : ('I_e'    , "parameters['i_offset']*1000.0"), # I_e is in pA, i_offset in nA 
     113            'v_init'    : ('V_m'    , "parameters['v_init']"), 
    87114    } 
    88     nest_name = "iaf_neuron2" 
     115    nest_name = "iaf_psc_alpha" 
    89116     
    90117    def __init__(self,parameters): 
    …… 
    121148     
    122149    translations = { 
    123             'v_rest'    : ('U0'          , "parameters['v_rest']"), 
    124             'v_reset'   : ('Vreset'      , "parameters['v_reset']"), 
    125             'cm'        : ('C'           , "parameters['cm']*1000.0"), # C is in pF, cm in nF 
    126             'tau_m'     : ('Tau'         , "parameters['tau_m']"), 
    127             'tau_refrac': ('TauR'        , "max(dt,parameters['tau_refrac'])"), 
    128             'tau_syn_E' : ('TauSyn_E'    , "parameters['tau_syn_E']"), 
    129             'tau_syn_I' : ('TauSyn_I'    , "parameters['tau_syn_I']"), 
    130             'v_thresh'  : ('Theta'       , "parameters['v_thresh']"), 
    131             'i_offset'  : ('Istim'       , "parameters['i_offset']*1000.0"), # I0 is in pA, i_offset in nA 
    132             'e_rev_E'   : ('V_reversal_E', "parameters['e_rev_E']"), 
    133             'e_rev_I'   : ('V_reversal_I', "parameters['e_rev_I']"), 
    134             'v_init'    : ('u'           , "parameters['v_init']"), 
     150            'v_rest'    : ('E_L'       , "parameters['v_rest']"), 
     151            'v_reset'   : ('V_reset'   , "parameters['v_reset']"), 
     152            'cm'        : ('C_m'       , "parameters['cm']*1000.0"), # C is in pF, cm in nF 
     153            'tau_m'     : ('g_L'       , "parameters['cm']/parameters['tau_m']*1000.0"), 
     154            'tau_refrac': ('t_ref'     , "max(dt,parameters['tau_refrac'])"), 
     155            'tau_syn_E' : ('tau_syn_ex', "parameters['tau_syn_E']"), 
     156            'tau_syn_I' : ('tau_syn_in', "parameters['tau_syn_I']"), 
     157            'v_thresh'  : ('V_th'      , "parameters['v_thresh']"), 
     158            #'i_offset'  : ('Istim'    , "parameters['i_offset']*1000.0"), # I0 is in pA, i_offset in nA 
     159            'e_rev_E'   : ('E_ex'      , "parameters['e_rev_E']"), 
     160            'e_rev_I'   : ('E_in'      , "parameters['e_rev_I']"), 
     161            'v_init'    : ('V_m'       , "parameters['v_init']"), 
    135162    } 
    136163    nest_name = "iaf_cond_alpha" 
    …… 
    140167                                                       # values for not-specified parameters. 
    141168        self.parameters = self.translate(self.parameters) 
    142         self.parameters['gL'] = self.parameters['C']/self.parameters['Tau'] # Trick to fix the leak conductance 
    143  
     169         
    144170 
    145171class IF_cond_exp(common.IF_cond_exp): 
    …… 
    148174     
    149175    translations = { 
    150             'v_rest'    : ('U0'          , "parameters['v_rest']"), 
    151             'v_reset'   : ('Vreset'      , "parameters['v_reset']"), 
    152             'cm'        : ('C'           , "parameters['cm']*1000.0"), # C is in pF, cm in nF 
    153             'tau_m'     : ('Tau'         , "parameters['tau_m']"), 
    154             'tau_refrac': ('TauR'        , "max(dt,parameters['tau_refrac'])"), 
    155             'tau_syn_E' : ('TauSyn_E'    , "parameters['tau_syn_E']"), 
    156             'tau_syn_I' : ('TauSyn_I'    , "parameters['tau_syn_I']"), 
    157             'v_thresh'  : ('Theta'       , "parameters['v_thresh']"), 
    158             'i_offset'  : ('Istim'       , "parameters['i_offset']*1000.0"), # I0 is in pA, i_offset in nA 
    159             'e_rev_E'   : ('V_reversal_E', "parameters['e_rev_E']"), 
    160             'e_rev_I'   : ('V_reversal_I', "parameters['e_rev_I']"), 
    161             'v_init'    : ('u'           , "parameters['v_init']"), 
     176            'v_rest'    : ('E_L'          , "parameters['v_rest']"), 
     177            'v_reset'   : ('V_reset'      , "parameters['v_reset']"), 
     178            'cm'        : ('C_m'           , "parameters['cm']*1000.0"), # C is in pF, cm in nF 
     179            'tau_m'     : ('g_L'         , "parameters['cm']/parameters['tau_m']*1000.0"), 
     180            'tau_refrac': ('t_ref'        , "max(dt,parameters['tau_refrac'])"), 
     181            'tau_syn_E' : ('tau_syn_ex'    , "parameters['tau_syn_E']"), 
     182            'tau_syn_I' : ('tau_syn_in'    , "parameters['tau_syn_I']"), 
     183            'v_thresh'  : ('V_th'       , "parameters['v_thresh']"), 
     184            #'i_offset'  : ('Istim'       , "parameters['i_offset']*1000.0"), # I0 is in pA, i_offset in nA 
     185            'e_rev_E'   : ('E_ex', "parameters['e_rev_E']"), 
     186            'e_rev_I'   : ('E_in', "parameters['e_rev_I']"), 
     187            'v_init'    : ('V_m'           , "parameters['v_init']"), 
    162188    } 
    163189    nest_name = "iaf_cond_exp" 
    …… 
    167193                                                       # values for not-specified parameters. 
    168194        self.parameters = self.translate(self.parameters) 
    169         self.parameters['gL'] = self.parameters['C']/self.parameters['Tau'] # Trick to fix the leak conductance 
    170  
     195         
    171196 
    172197class HH_cond_exp(common.HH_cond_exp): 
    …… 
    187212        'tau_syn_I' : ('tau_in', "parameters['tau_syn_I']"), 
    188213        'i_offset'  : ('I_stim', "parameters['i_offset']*1000.0"), 
     214        'v_init'    : ('V_m',    "parameters['v_init']"), 
    189215    } 
    190216    nest_name = "hh_cond_exp_traub" 
    …… 
    195221        self.parameters = self.translate(self.parameters) 
    196222         
     223class AdaptiveExponentialIF_alpha(common.AdaptiveExponentialIF_alpha): 
     224    """adaptive exponential integrate and fire neuron according to Brette and Gerstner (2005)""" 
     225     
     226    translations = { 
     227        'v_init'    : ('V_m',        "parameters['v_init']"), 
     228        'w_init'    : ('w',          "parameters['w_init']*1000.0"), # nA -> pA 
     229        'cm'        : ('C_m',        "parameters['cm']*1000.0"),     # nF -> pF 
     230        'tau_refrac': ('t_ref',      "parameters['tau_refrac']"),  
     231        'v_spike'   : ('V_peak',     "parameters['v_spike']"), 
     232        'v_reset'   : ('V_reset',    "parameters['v_reset']"), 
     233        'v_rest'    : ('E_L',        "parameters['v_rest']"), 
     234        'tau_m'     : ('g_L',        "parameters['cm']/parameters['tau_m']*1000.0"), 
     235        'i_offset'  : ('I_e',        "parameters['i_offset']*1000.0"), # nA -> pA 
     236        'a'         : ('a',          "parameters['a']"),        
     237        'b'         : ('b',          "parameters['b']*1000.0"),  # nA -> pA. 
     238        'delta_T'   : ('Delta_T',    "parameters['delta_T']"),  
     239        'tau_w'     : ('tau_w',      "parameters['tau_w']"),  
     240        'v_thresh'  : ('V_th',       "parameters['v_thresh']"),  
     241        'e_rev_E'   : ('E_ex',       "parameters['e_rev_E']"), 
     242        'tau_syn_E' : ('tau_syn_ex', "parameters['tau_syn_E']"),  
     243        'e_rev_I'   : ('E_in',       "parameters['e_rev_I']"),  
     244        'tau_syn_I' : ('tau_syn_in', "parameters['tau_syn_I']"), 
     245    } 
     246    nest_name = "aeif_cond_alpha" 
     247     
     248    def __init__(self,parameters): 
     249        common.AdaptiveExponentialIF_alpha.__init__(self,parameters) 
     250        self.parameters = self.translate(self.parameters) 
    197251         
    198252class SpikeSourcePoisson(common.SpikeSourcePoisson): 
    …… 
    229283# ============================================================================== 
    230284 
    231 def setup(timestep=0.1,debug=False,**extra_params): 
     285def setup(timestep=0.1,min_delay=0.1,max_delay=0.1,debug=False,**extra_params): 
    232286    """ 
    233287    Should be called at the very beginning of a script. 
    …… 
    239293    global dt 
    240294    global tempdir 
    241     global hl_spike_files, hl_v_files 
     295    global _min_delay 
     296    #global hl_spike_files, hl_v_files 
    242297    dt = timestep 
    243  
     298    _min_delay = min_delay 
     299     
    244300    # reset the simulation kernel 
    245301    nest.ResetKernel() 
    …… 
    247303    nest.sr('clear') 
    248304 
    249     # check if hl_spike_files , hl_v_files are empty 
    250     if not len(hl_spike_files) == 0: 
    251         print 'hl_spike_files still contained files, please close all open files before setup' 
    252         hl_spike_files = {} 
    253     if not len(hl_v_files) == 0: 
    254         print 'hl_v_files still contained files, please close all open files before setup' 
    255         hl_v_files = {} 
     305#    # check if hl_spike_files , hl_v_files are empty 
     306#    if not len(hl_spike_files) == 0: 
     307#        print 'hl_spike_files still contained files, please close all open files before setup' 
     308#        hl_spike_files = {} 
     309#    if not len(hl_v_files) == 0: 
     310#        print 'hl_v_files still contained files, please close all open files before setup' 
     311#        hl_v_files = {} 
    256312     
    257313    tempdir = tempfile.mkdtemp() 
    258314    tempdirs.append(tempdir) # append tempdir to tempdirs list 
     315     
    259316    # set tempdir 
    260     nest.SetStatus([0],[{'device_prefix':tempdir}]) 
     317    nest.SetStatus([0], {'device_prefix':tempdir}) 
    261318    # set resolution 
    262     nest.SetStatus([0],[{'resolution': dt}]) 
     319    nest.SetStatus([0], {'resolution': dt}) 
    263320     
    264321    if extra_params.has_key('threads'): 
    …… 
    295352 
    296353    logging.info("Initialization of Nest")     
    297     return 0 
     354    return nest.Rank() 
    298355 
    299356def end(compatible_output=True): 
    …… 
    303360 
    304361    # NEST will soon close all its output files after the simulate function is over, therefore this step is not necessary 
    305     global tempdir 
     362    global tempdirs 
    306363     
    307364    # And we postprocess the low level files opened by record() 
    308365    # and record_v() method 
    309     for file in ll_spike_files: 
    310         _printSpikes(file, compatible_output) 
    311     for file in ll_v_files: 
    312         _print_v(file, compatible_output) 
     366    #for file in ll_spike_files: 
     367    #    _printSpikes(file, compatible_output) 
     368    #for file, nest_file in ll_v_files: 
     369    #    _print_v(file, nest_file, compatible_output) 
     370    print "Saving the following files:", recorder_dict.keys() 
     371    for filename in recorder_dict.keys(): 
     372        _print(filename, gather=False, compatible_output=compatible_output) 
     373     
    313374    for tempdir in tempdirs: 
    314375        os.system("rm -rf %s" %tempdir) 
    315     nest.end() 
    316376 
    317377def run(simtime): 
    …… 
    364424        weight = 0.0 
    365425    if delay is None: 
    366         delay = nest.getLimits()['min_delay'] 
     426        delay = _min_delay 
    367427    weight = weight*1000 # weights should be in nA or uS, but iaf_neuron uses pA and iaf_cond_neuron uses nS. 
    368428                         # Using convention in this way is not ideal. We should be able to look up the units used by each model somewhere. 
    …… 
    371431    try: 
    372432        if type(source) != types.ListType and type(target) != types.ListType: 
    373             connect_id = nest.ConnectWD([source],[target],[weight],[delay]) 
     433            nest.ConnectWD([source],[target],[weight],[delay]) 
     434            connect_id = Connection(source, target) 
    374435        else: 
    375436            connect_id = [] 
    …… 
    386447                for j,tgt in enumerate(target): 
    387448                    if p >= 1 or rarr[j] < p: 
    388                         connect_id += [nest.ConnectWD(src,tgt,[weight],[delay])] 
     449                        nest.ConnectWD([src],[tgt],[weight],[delay]) 
     450                        connect_id += [Connection(src,tgt)] 
    389451    except nest.SLIError: 
    390452        raise common.ConnectionError 
    …… 
    410472    nest.SetStatus(cells,[param]) 
    411473 
    412 def record(source,filename): 
     474def _connect_recording_device(recorder, record_from=None): 
     475    device = nest.GetStatus(recorder, "model")[0] 
     476    if device == "spike_detector": 
     477        nest.ConvergentConnect(record_from, recorder) 
     478    elif device in ('voltmeter', 'conductancemeter'):         
     479        nest.DivergentConnect(recorder, record_from) 
     480    else: 
     481        raise Exception("Not a valid recording device") 
     482 
     483def _record(variable, source, filename): 
    413484    """Record spikes to a file. source can be an individual cell or a list of 
    414485    cells.""" 
    415486    # would actually like to be able to record to an array and choose later 
    416487    # whether to write to a file. 
    417     spike_detector = nest.Create('spike_detector') 
    418     nest.setDict(spike_detector,{'withtime':True,  # record time of spikes 
    419                                    'withpath':True}) # record which neuron spiked 
    420     if type(source) == types.ListType: 
    421         source = [nest.getAddress(src) for src in source] 
    422     else: 
    423         source = [nest.getAddress(source)] 
    424     for src in source: 
    425         nest.connect(src,spike_detector[0]) 
    426         nest.sr('/%s (%s/%s) (w) file def' % (filename, tempdir, filename)) 
    427         nest.sr('%s << /output_stream %s >> SetStatus' % (nest.getGID(spike_detector[0]),filename)) 
    428     ll_spike_files.append(filename) 
    429  
    430  
    431 def record_v(source,filename_user): 
     488    device_name = recording_device_names[variable] 
     489    print "Trying to record %s from cell %s using a %s" % (variable, source, device_name) 
     490    recording_device = nest.Create(device_name) # does it need to be an ID? 
     491    nest.SetStatus(recording_device, 
     492                   {"to_file" : True, "withgid" : True, "withtime" : True, 
     493                    "interval": nest.GetStatus([0], "resolution")[0]}) 
     494             
     495    if type(source) != types.ListType: 
     496        source = [source] 
     497    _connect_recording_device(recording_device, record_from=source) 
     498    recorder_dict[filename] = recording_device  
     499 
     500def record(source, filename): 
     501    """Record spikes to a file. source can be an individual cell or a list of 
     502    cells.""" 
     503    # would actually like to be able to record to an array and choose later 
     504    # whether to write to a file. 
     505    _record('spikes', source, filename)  
     506     
     507def record_v(source, filename): 
    432508    """ 
    433509    Record membrane potential to a file. source can be an individual cell or 
    …… 
    435511    # would actually like to be able to record to an array and 
    436512    # choose later whether to write to a file. 
    437     if type(source) != types.ListType: 
    438         source = [source] 
    439     record_file = filename_user 
    440     #ll_v_files.append(filename) 
    441     filename = nest.record_v(source,record_file) 
    442     ll_v_files.append(filename) 
    443  
    444 def _printSpikes(filename, compatible_output=True): 
    445     """ Print spikes into a file, and postprocessed them if 
    446     needed and asked to produce a compatible output for all the simulator 
    447     Should actually work with record() and allow to dissociate the recording of the 
    448     writing process, which is not the case for the moment""" 
    449     tempfilename = "%s/%s" %(tempdir, filename) 
    450     nest.sr('%s close' %filename)  
    451     if (compatible_output): 
    452         # Here we postprocess the file to have effectively the 
    453         # desired format : 
    454         # First line: # dimensions of the population 
    455         # Then spiketime (in ms) cell_id-min(cell_id) 
    456         result = open(filename,'w',1000) 
    457         # Writing # such that Population.printSpikes and this have same output format 
    458         result.write("# "+"\n") 
    459         # Writing spiketimes, cell_id-min(cell_id) 
    460         # Pylab has a great load() function, but it is not necessary to import 
    461         # it into pyNN. The fromfile() function of numpy has trouble on several 
    462         # machine with Python 2.5, so that's why a dedicated _readArray function 
    463         # has been created to load from file the raster or the membrane potentials 
    464         # saved by NEST 
    465         try: 
    466             raster = _readArray(tempfilename, sepchar=None) 
    467             raster = raster[:,1:3] 
    468             raster[:,1] = raster[:,1]*dt 
    469             for idx in xrange(len(raster)): 
    470                 result.write("%g\t%d\n" %(raster[idx][1], raster[idx][0])) 
    471         except Exception: 
    472             print "Error while writing data into a compatible mode" 
    473         result.close() 
    474         os.system("rm %s" %tempfilename) 
     513    _record('v', source, filename)  
     514 
     515def _print(user_filename, gather=True, compatible_output=True, population=None, variable=None): 
     516    global recorder_dict 
     517     
     518    if population is None: 
     519        recorder = recorder_dict[user_filename] 
    475520    else: 
    476         shutil.move(tempfilename, filename) 
    477  
    478  
    479 def _print_v(filename, compatible_output=True): 
    480     """ Print membrane potentials in a file, and postprocessed them if 
    481     needed and asked to produce a compatible output for all the simulator 
    482     Should actually work with record_v() and allow to dissociate the recording of the 
    483     writing process, which is not the case for the moment""" 
    484     tempfilename = tempdir+'/'+filename 
    485     nest.sr('%s close' %tempfilename.replace('/','_'))  
    486     result = open(filename,'w',1000) 
    487     dt = nest.getNESTStatus()['resolution'] 
    488     n = int(nest.getNESTStatus()['time']/dt) 
    489     result.write("# dt = %f\n# n = %d\n" % (dt,n)) 
    490     if (compatible_output): 
    491         # Here we postprocess the file to have effectively the 
    492         # desired format : 
    493         # First line: dimensions of the population 
    494         # Then spiketime cell_id-min(cell_id) 
    495  
    496         # Pylab has a great load() function, but it is not necessary to import 
    497         # it into pyNN. The fromfile() function of numpy has trouble on several 
    498         # machine with Python 2.5, so that's why a dedicated _readArray function 
    499         # has been created to load from file the raster or the membrane potentials 
    500         # saved by NEST 
    501         try: 
    502             raster = _readArray(tempfilename.replace('/','_'), sepchar=None) 
    503             for idx in xrange(len(raster)): 
    504                 result.write("%g\t%d\n" %(raster[idx][1], raster[idx][0])) 
    505         except Exception: 
    506             print "Error while writing data into a compatible mode" 
    507     else: 
    508         f = open(tempfilename.replace('/','_'),'r',1000) 
    509         lines = f.readlines() 
    510         f.close() 
    511         for line in lines: 
    512             result.write(line) 
    513     result.close() 
    514     os.system("rm %s" %tempfilename.replace('/','_')) 
    515  
    516 def _readArray(filename, sepchar = None, skipchar = '#'): 
     521        assert variable in ['spikes', 'v', 'conductance'] 
     522        recorder = population.recorders[variable] 
     523     
     524    nest.FlushDevice(recorder)  
     525    status = nest.GetStatus([0])[0] 
     526    local_num_threads = status['local_num_threads'] 
     527    node_list = range(nest.GetStatus([0], "num_processes")[0]) 
     528     
     529    # First combine data from different threads 
     530    os.system("rm -f %s" % user_filename) 
     531    for nest_thread in range(local_num_threads): 
     532        nest.sps(recorder[0]) 
     533        nest.sr("%i GetAddress %i append" % (recorder[0], nest_thread)) 
     534        nest.sr("GetStatus /filename get") 
     535        nest_filename = nest.spp() #nest.GetStatus(recorder, "filename") 
     536        ##system_line = 'cat %s >> %s' % (nest_filename, "%s_%d" % (user_filename, nest.Rank())) 
     537        merged_filename = "%s/%s" % (os.path.dirname(nest_filename), user_filename) 
     538        system_line = 'cat %s >> %s' % (nest_filename, merged_filename) # will fail if writing to a common directory, e.g. using NFS 
     539        print system_line 
     540        os.system(system_line) 
     541    if gather: 
     542        raise Exception("'gather' not currently supported.") 
     543        if nest.Rank() == 0: # only on the master node (?) 
     544            for node in node_list: 
     545                pass # not a good way to do it at the moment 
     546     
     547    if compatible_output: 
     548        if gather == False or nest.Rank() == 0: # if we gather, only do this on the master node 
     549            logging.info("Writing %s in compatible format." % user_filename) 
     550             
     551            # Here we postprocess the file to have effectively the 
     552            # desired format: spiketime (in ms) cell_id-min(cell_id) 
     553            #if not os.path.exists(user_filename): 
     554            result = open(user_filename,'w',1000) 
     555            #else: 
     556            #    result = open(user_filename,'a',1000) 
     557                 
     558            ## Writing header info (e.g., dimensions of the population) 
     559            if population is not None: 
     560                result.write("# " + "\t".join([str(d) for d in self.dim]) + "\n") 
     561                padding = population.cell.flatten()[0] 
     562            else: 
     563                padding = 0 
     564            result.write("# dt = %g\n" % nest.GetStatus([0], "resolution")[0]) 
     565                             
     566            # Writing spiketimes, cell_id-min(cell_id) 
     567                  
     568            # (Pylab has a great load() function, but it is not necessary to import 
     569            # it into pyNN. The fromfile() function of numpy has trouble on several 
     570            # machine with Python 2.5, so that's why a dedicated _readArray function 
     571            # has been created to load from file the raster or the membrane potentials 
     572            # saved by NEST). 
     573                     
     574            # open file 
     575            if int(os.path.getsize(merged_filename)) > 0: 
     576                raster = _readArray(merged_filename, sepchar=None) 
     577                result.write("# n = %d\n" % len(raster)) # shouldn't be called raster 
     578                raster[:,0] = raster[:,0] - padding 
     579                for idx in xrange(len(raster)): 
     580                    result.write("%g\t%d\n" % (raster[idx][2], raster[idx][0])) # v id 
     581            else: 
     582                logging.info("%s_tmp is empty" % user_filename) 
     583            result.close() 
     584     
     585    recorder_dict.pop(user_filename)     
     586 
     587def _readArray(filename, sepchar=None, skipchar='#'): 
    517588    logging.debug(filename) 
    518589    myfile = open(filename, "r") 
    …… 
    526597            if (stripped_line[0] != skipchar): 
    527598                items = stripped_line.split(sepchar) 
    528                 # Here we have to deal with the fact that quite often, NEST 
    529                 # does not write correctly the last line of Vm recordings. 
    530                 # More precisely, it is often not complete 
    531                 #try : 
     599                if len(items) != 3: 
     600                    print stripped_line 
     601                    print items 
     602                    raise Exception() 
    532603                data.append(map(float, items)) 
    533                 #except Exception: 
    534                 #    # The last line has a gid and just a "-" sign... 
    535                 #    pass 
    536604    try : 
    537605        a = numpy.array(data) 
    538606    except Exception: 
     607        raise 
    539608        # The last line has just a gid, so we has to remove it 
    540         a = numpy.array(data[0:len(data)-2]) 
     609        #a = numpy.array(data[0:len(data)-2]) 
    541610    (Nrow,Ncol) = a.shape 
    542611    logging.debug(str(a.shape)) 
    …… 
    596665        if not self.label: 
    597666            self.label = 'population%d' % Population.nPop 
     667        self.recorders = {'spikes': None, 'v': None, 'conductance': None} 
    598668        Population.nPop += 1 
    599669     
    …… 
    756826        raise Exception("Method not yet implemented") 
    757827 
    758     def record(self,record_from=None,rng=None): 
    759         """ 
    760         If record_from is not given, record spikes from all cells in the Population. 
    761         record_from can be an integer - the number of cells to record from, chosen 
    762         at random (in this case a random number generator can also be supplied) 
    763         - or a list containing the ids 
    764         of the cells to record. 
    765         """ 
    766         global hl_spike_files 
    767          
     828    def _record(self, variable, record_from=None, rng=None): 
     829        assert variable in ('spikes', 'v', 'conductance') 
     830     
    768831        # create device 
    769         self.spike_detector = ID(nest.Create('spike_detector')[0]) 
    770         params = {"to_file" : True, "withgid" : True, "withtime" : True}#,'withpath':True} 
    771         nest.SetStatus([self.spike_detector], [params]) 
    772          
    773         filename = nest.GetStatus([self.spike_detector], "filename") 
    774         hl_spike_files[self.label] = filename 
    775          
     832        device_name = recording_device_names[variable] 
     833        if self.recorders[variable] is None: 
     834            self.recorders[variable] = ID(nest.Create(device_name)[0]) # does it need to be an ID? 
     835            nest.SetStatus([self.recorders[variable]], 
     836                           {"to_file" : True, "withgid" : True, "withtime" : True})       
     837                 
    776838        # create list of neurons         
    777839        fixed_list = False 
    …… 
    787849            n_rec = self.size 
    788850         
    789  
    790851        tmp_list = [] 
    791852        if (fixed_list == True): 
    792853            for neuron in record_from: 
    793854                tmp_list = [neuron for neuron in record_from] 
    794                 #nest.Connect([neuron],[self.spike_detector[0]]) 
    795855        else: 
    796856            for neuron in numpy.random.permutation(numpy.reshape(self.cell,(self.cell.size,)))[0:n_rec]: 
    797857                tmp_list.append(neuron) 
    798                 #nest.Connect([neuron],[self.spike_detector[0]]) 
    799  
     858                 
    800859        # connect device to neurons 
    801         nest.ConvergentConnect(tmp_list,[self.spike_detector]) 
    802          
    803          
    804         #hl_spike_files.append('%s.spikes' % self.label) 
    805         self.n_rec = n_rec 
    806  
     860        _record(variable, tmp_list) 
     861 
     862    def record(self, record_from=None, rng=None): 
     863        """ 
     864        If record_from is not given, record spikes from all cells in the Population. 
     865        record_from can be an integer - the number of cells to record from, chosen 
     866        at random (in this case a random number generator can also be supplied) 
     867        - or a list containing the ids 
     868        of the cells to record. 
     869        """ 
     870        self._record('spikes', record_from, rng) 
     871         
    807872    def record_v(self,record_from=None,rng=None): 
    808873        """ 
    …… 
    813878        - or a list containing the ids of the cells to record. 
    814879        """ 
    815         global hl_v_files 
    816          
    817         # create device 
    818          
    819         self.voltmeter = nest.Create("voltmeter") 
    820         params = {"to_file" : True, "withgid" : True, "withtime" : True} 
    821         nest.SetStatus(self.voltmeter, [params]) 
    822          
    823         filename = nest.GetStatus(self.voltmeter, "filename") 
    824         hl_v_files[self.label] = filename 
    825          
    826          
    827         # create list of neurons 
    828         fixed_list = False 
    829         if record_from: 
    830             if type(record_from) == types.ListType: 
    831                 fixed_list = True 
    832                 n_rec = len(record_from) 
    833             elif type(record_from) == types.IntType: 
    834                 n_rec = record_from 
    835             else: 
    836                 raise "record_from must be a list or an integer" 
    837         else: 
    838             n_rec = self.size 
    839  
    840         tmp_list = [] 
    841         if (fixed_list == True): 
    842             tmp_list = [neuron for neuron in record_from] 
    843         else: 
    844             for neuron in numpy.random.permutation(numpy.reshape(self.cell,(self.cell.size,)))[0:n_rec]: 
    845                 tmp_list.append(neuron) 
    846          
    847         # connect device to neurons 
    848         nest.DivergentConnect(self.voltmeter, tmp_list) 
    849          
     880        self._record('v', record_from, rng) 
    850881     
    851882    def record_c(self,record_from=None,rng=None): 
    …… 
    857888        - or a list containing the ids of the cells to record. 
    858889        """ 
    859         global hl_c_files 
    860          
    861         # create device 
    862          
    863         self.conductancemeter = nest.Create("conductancemeter") 
    864         params = {"to_file" : True, "withgid" : True, "withtime" : True} 
    865         nest.SetStatus(self.conductancemeter, [params]) 
    866          
    867         filename = nest.GetStatus(self.conductancemeter, "filename") 
    868         hl_c_files[self.label] = filename 
    869          
    870          
    871         # create list of neurons 
    872         fixed_list = False 
    873         if record_from: 
    874             if type(record_from) == types.ListType: 
    875                 fixed_list = True 
    876                 n_rec = len(record_from) 
    877             elif type(record_from) == types.IntType: 
    878                 n_rec = record_from 
    879             else: 
    880                 raise "record_from must be a list or an integer" 
    881         else: 
    882             n_rec = self.size 
    883  
    884         tmp_list = [] 
    885         if (fixed_list == True): 
    886             tmp_list = [neuron for neuron in record_from] 
    887         else: 
    888             for neuron in numpy.random.permutation(numpy.reshape(self.cell,(self.cell.size,)))[0:n_rec]: 
    889                 tmp_list.append(neuron) 
    890          
    891         # connect device to neurons 
    892         nest.DivergentConnect(self.conductancemeter, tmp_list) 
    893      
    894     def printSpikes(self,filename,gather=True, compatible_output=False): 
     890        self._record('conductance', record_from, rng) 
     891     
     892    def printSpikes(self, filename, gather=True, compatible_output=False): 
    895893        """ 
    896894        Writes spike times to file. 
    …… 
    904902         
    905903        If compatible_output is False, the raw format produced by the simulator 
    906         is used. This may be faster, sinc# Writing spiketimes, cell_id-min(cell_id) 
     904        is used. This may be faster. 
    907905        If gather is True, the file will only be created on the master node, 
    908906        otherwise, a file will be written on each node. 
    909         """         
    910         global hl_spike_files 
    911         global tempdir 
    912         # closing of the file 
    913         # just a workaround, nest will do that automatically soon 
    914         if hl_spike_files.has_key(self.label):#   __contains__(tempfilename): 
    915             nest.sps(self.spike_detector) 
    916             nest.sr("FlushDevice") 
    917          
    918         status = nest.GetStatus([0])[0] 
    919         np = status['num_processes'] 
    920         vp = status['vp'] 
    921         local_num_threads = status['local_num_threads'] 
    922         rank = numpy.mod(vp,np) 
    923  
    924          
    925         filename = filename + '-%d' % rank 
    926         if (compatible_output): 
    927             if rank == 0: 
    928                 logging.info("Writing %s in compatible format." % filename) 
    929                 for nest_thread in range(local_num_threads): 
    930                     label = '-%d-%d-%d.gdf' %(rank,nest_thread,self.spike_detector) 
    931                     # Here we postprocess the file to have effectively the 
    932                     # desired format: spiketime (in ms) cell_id-min(cell_id) 
    933                     if not os.path.exists(filename): 
    934                         result = open(filename,'w',1000) 
    935                         # Writing dimensions of the population: 
    936                         result.write("# " + "\t".join([str(d) for d in self.dim]) + "\n") 
    937                         # Writing spiketimes, cell_id-min(cell_id) 
    938                         ###padding = numpy.reshape(self.cell,self.cell.size)[0] # moved below 
    939                     else: 
    940                         result = open(filename,'a',1000) 
    941                     padding = numpy.reshape(self.cell,self.cell.size)[0] 
    942                      
    943                     # Pylab has a great load() function, but it is not necessary to import 
    944                     # it into pyNN. The fromfile() function of numpy has trouble on several 
    945                     # machine with Python 2.5, so that's why a dedicated _readArray function 
    946                     # has been created to load from file the raster or the membrane potentials 
    947                     # saved by NEST 
    948                     # open file 
    949                     simfile = tempdir + '/spike_detector'+ label 
    950                     # if int(os.path.getsize(hl_spike_files[self.label][0])) > 0: 
    951                     if int(os.path.getsize(simfile)) > 0: 
    952 #                        try: 
    953                         raster = _readArray(simfile ,sepchar=None) 
    954                         # Sometimes, nest doesn't write the last line entirely, so we need 
    955                         # to trunk it to avoid errors 
    956                         ###raster = raster[:,1:3] 
    957                         raster[:,0] = raster[:,0] - padding 
    958                         raster[:,1] = raster[:,1]*dt 
    959                         for idx in xrange(len(raster)): 
    960                             result.write("%g\t%d\n" %(raster[idx][1], raster[idx][0])) 
    961                                  
    962 #                        except Exception, inst: 
    963 #                            print "Error while writing data into a compatible mode" 
    964 #                            logging.error("Error while writing data into a compatible mode: %s" % str(inst)) 
    965                     else: 
    966                         logging.info("%s is empty" % simfile) 
    967                 result.close() 
    968                 # os.system("rm %s" % hl_spike_files[self.label][0]) 
    969         else: 
    970             if gather: 
    971                 if os.path.exists(filename): 
    972                     # 'target file exists, will be removed before copying the new data.' 
    973                     os.remove(filename) 
    974                 for nest_thread in range(local_num_threads): 
    975                     label = '-%d-%d-%d.gdf' %(rank,nest_thread,self.spike_detector) 
    976                     simfile = tempdir + '/spike_detector'+ label 
    977                     system_line = 'cat %s >> %s' %(simfile,filename) 
    978                     if os.system(system_line) == 0: # cat was successful 
    979                         os.remove(simfile) 
    980             else: 
    981                 print 'spike data is not gathered and located in: ',tempdir 
    982          
    983         hl_spike_files.pop(self.label) 
    984  
    985     def collectdata(self,filename): 
    986         """ 
    987         merges all mpi data into one final file, should only be called once, after all spikes have been printed 
    988         """ 
    989         status = nest.GetStatus([0])[0] 
    990         np = status['num_processes'] 
    991         vp = status['vp'] 
    992         local_num_threads = status['local_num_threads'] 
    993         rank = numpy.mod(vp,np) 
    994         if os.path.exists(filename): 
    995             os.remove(filename) 
    996          
    997         for processor in range(np): 
    998             filename_p = filename+'-'+str(processor) 
    999             system_line = 'cat %s >> %s' %(filename_p,filename) 
    1000             if os.system(system_line) == 0: # cat was successful 
    1001                 os.remove(filename_p)     
    1002          
    1003     def meanSpikeCount(self,gather=True): 
     907        """ 
     908        _print(filename, gather=gather, compatible_output=compatible_output, 
     909               population=self, variable="spikes") 
     910        
     911    def meanSpikeCount(self, gather=True): 
    1004912        """ 
    1005913        Returns the mean number of spikes per neuron. 
    …… 
    1010918        return float(n_spikes)/self.n_rec 
    1011919 
    1012     def randomInit(self,rand_distr): 
     920    def randomInit(self, rand_distr): 
    1013921        """ 
    1014922        Sets initial membrane potentials for all the cells in the population to 
    …… 
    1016924        """ 
    1017925        self.rset('v_init',rand_distr) 
    1018         #cells = numpy.reshape(self.cell,self.cell.size) 
    1019         #rvals = rand_distr.next(n=self.cell.size) 
    1020         #for node, v_init in zip(cells,rvals): 
    1021         #    nest.setDict([node],{'u': v_init}) 
    1022      
    1023     def print_v(self,filename,gather=True, compatible_output=False): 
     926 
     927    def print_v(self, filename, gather=True, compatible_output=False): 
    1024928        """ 
    1025929        Write membrane potential traces to file. 
    1026         If compatible_output is True, the format is "v cell_id", 
     930        If compatible_output is True, the format is "t v cell_id", 
    1027931        where cell_id is the index of the cell counting along rows and down 
    1028932        columns (and the extension of that for 3-D). 
    1029         This allows easy plotting of a `raster' plot of spiketimes, with one 
    1030         line for each cell. 
    1031933        The timestep and number of data points per cell is written as a header, 
    1032934        indicated by a '#' at the beginning of the line. 
    …… 
    1036938        voltage files. 
    1037939        """ 
    1038         global hl_v_files 
    1039         &nb