Changeset 340

Show
Ignore:
Timestamp:
11/14/08 19:17:05 (2 months ago)
Author:
LaurentPerrinet
Message:

some polishing of the wiki:examples, the retina needs a bit more work...

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • trunk/examples/retina/benchmark_linear.py

    r333 r340  
    1616import os, sys, numpy, pylab, shelve 
    1717 
     18N, N_exp = 1000, 6 
     19t_smooth = 100. # width (in ms) of the integration window 
     20from NeuroTools.parameters import ParameterSpace, ParameterRange 
     21snr  = 2.0 * numpy.linspace(0.1,2.0,N_exp) 
     22p = ParameterSpace({'snr' : ParameterRange(list(snr))}) 
     23 
    1824 
    1925name = sys.argv[0].split('.')[0] # name of the current script withpout the '.py' part 
    …… 
    2834 
    2935except: 
    30     N, N_exp = 1000, 6 
    31     t_smooth = 100. # width (in ms) of the integration window 
    3236    from retina import * 
    3337    retina = Retina(N) 
    3438    retina.params['amplitude'] = numpy.ones(retina.params['amplitude'].shape) 
    35     snr  = retina.params['snr']* numpy.linspace(0.1,2.0,N_exp) 
    36     from NeuroTools.parameters import ParameterSpace 
    37     p = ParameterSpace({'snr' : ParameterRange(list(snr))}) 
     39     
    3840 
    3941    # calculates the dimension of the parameter space 
    …… 
    6567    results['temporal_ON'] = temporal_ON 
    6668    results['temporal_OFF'] = temporal_OFF 
    67     results['params'] = params 
     69    results['params'] = retina.params 
    6870 
    6971    pbar.finish() 
    …… 
    8486#pylab.rcParams.update(pylab_params(fig_width_pt = 497.9/2., ratio = 1.)) 
    8587pylab.figure(1) 
    86 fmax = numpy.max([numpy.max(temporal_OFF[:]),numpy.max(temporal_ON[:])]) 
     88#fmax = numpy.max([numpy.max(temporal_OFF[:]),numpy.max(temporal_ON[:])]) 
    8789 
    8890pylab.subplot(211) 
    8991for i_exp in range(N_exp): 
    90     pylab.plot(lower_edges[:-1] + t_smooth/2,  # TODO add a half bin 
    91                     temporal_ON[i_exp])# 
    92 pylab.xticks( numpy.round(numpy.linspace(0, retina.params.simtime, 5),0) ) 
    93 pylab.ylabel('ON Firing frequency (Hz/neuron)') 
    94 pylab.axis([0, retina.params.simtime, 0.0, fmax]) 
     92    pylab.plot(lower_edges[:-1] + t_smooth/2, temporal_ON[i_exp], 
     93                        label= '%5.2f' % p.snr._values[i_exp]) 
     94pylab.xticks( numpy.round(numpy.linspace(0, params.simtime, 5),0) ) 
     95pylab.ylabel('ON Firing frequency (Hz)') 
     96pylab.axis([0, params.simtime, 0.0, numpy.max(temporal_ON[:])]) 
     97pylab.legend(loc='upper right') 
    9598pylab.subplot(212) 
    9699for i_exp in range(N_exp): 
    97     pylab.plot(lower_edges[:-1] + t_smooth/2,  # TODO add a half bin 
    98                     temporal_OFF[i_exp], 
    99                         label= 'snr= %5.3f' % p.snr._values[i_exp] ) 
    100 pylab.xticks( numpy.round(numpy.linspace(0, retina.params.simtime, 5),0) ) 
    101 pylab.ylabel('OFF Firing frequency (Hz/neuron)') 
     100    pylab.plot(lower_edges[:-1] + t_smooth/2, temporal_OFF[i_exp]) 
     101pylab.xticks( numpy.round(numpy.linspace(0, params.simtime, 5),0) ) 
     102pylab.ylabel('OFF Firing frequency (Hz)') 
    102103pylab.xlabel('time (ms)') 
    103 pylab.axis([0, retina.params.simtime, 0.0, fmax ]) 
    104  
    105 pylab.legend(loc='upper right') 
     104pylab.axis([0, params.simtime, 0.0, numpy.max(temporal_OFF[:]) ]) 
    106105 
    107106 

  • trunk/examples/retina/benchmark_noise.py

    r333 r340  
    1919import os, sys, numpy, shelve 
    2020 
    21 retina.params['snr'] = 0 # no input 
    22  
     21N, N_exp_noise = 1000, 22 
     22from NeuroTools.parameters import * 
     23p =  ParameterSpace({'noise_std' : ParameterRange(list(10.**(numpy.linspace(-.50,1.,N_exp_noise))))}) 
    2324 
    2425name = sys.argv[0].split('.')[0] # name of the current script withpout the '.py' part 
    …… 
    2728    CRF = results['CRF'] 
    2829except: 
    29          
    30     from NeuroTools.parameters import * 
    31  
    3230    # this is not mandatory but just a "easy_install progressbar" away 
    3331    # else remove all corresponding lines in this code... 
    3432    import progressbar # see http://projects.scipy.org/pipermail/scipy-dev/2008-January/008200.html 
    35  
    36     N, N_exp_noise = 1000, 22 
    37  
    38     p =  ParameterSpace({'noise_std' : ParameterRange(list(10.**(numpy.linspace(-.50,1.,N_exp_noise))))}) 
    39  
    4033    import retina as model 
    4134    retina = model.Retina(N)  
     35    retina.params['snr'] = 0 # no input 
     36     
    4237    # calculates the dimension of the parameter space 
    4338    results_dim, results_label = p.parameter_space_dimension_labels() 
    …… 
    7570""" 
    7671 
    77 N, simtime = retina.params['N'], retina.params['simtime'] 
    78  
    7972import pylab 
    8073 

  • trunk/examples/retina/benchmark_retina.py

    r333 r340  
    2525import os, sys, numpy, pylab, shelve, progressbar 
    2626 
     27N, N_snr, N_seeds = 1000, 5, 10 
     28from NeuroTools.parameters import * 
     29p =  ParameterSpace({ 
     30                'snr' : ParameterRange(list(numpy.linspace(0.1,2.0,N_snr))), 
     31                'kernelseed' : ParameterRange(list([12345 + k for k in range(N_seeds)]))}) 
    2732 
    2833name = sys.argv[0].split('.')[0] # name of the current script withpout the '.py' part 
    29 ############## MAKING THE SIMULATIONS ############### 
    30 results = shelve.open('results/mat-' + name) 
     34 
    3135try: 
    32     DATA = results['DATA'] 
    33 except: 
     36    ############## MAKING THE SIMULATIONS ############### 
     37    results = shelve.open('results/mat-' + name) 
     38    try: 
     39        DATA = results['DATA'] 
     40        params = results['params'] 
     41    except: 
     42        from retina import * 
     43        retina = Retina(N) 
     44        # calculates the dimension of the parameter space 
     45        results_dim, results_label = p.parameter_space_dimension_labels() 
    3446 
    35     from NeuroTools.parameters import * 
     47        DATA = [] 
     48        pbar=progressbar.ProgressBar(widgets=[name, " ", progressbar.Percentage(), ' ', 
     49                progressbar.Bar(), ' ', progressbar.ETA()], maxval=N_snr*N_seeds) 
     50        for i_exp, experiment in enumerate(p.iter_inner()): 
     51            params = retina.params 
     52            params.update(experiment) # updates what changed in the dictionary 
     53            # simulate the experiment and get its data 
     54            data = retina.run(params)#,verbose=False) 
     55            # store it 
     56            DATA.append(data)# 
     57            pbar.update(i_exp) 
    3658 
    37     N, N_snr, N_seeds = 1000, 5, 10 
    38     from retina import * 
    39     retina = Retina(N) 
     59        results['DATA'] = DATA 
     60        results['params'] = retina.params 
    4061 
    41     t_smooth = 100. # ms.  integration time to show fiber activity 
    42  
    43     p =  ParameterSpace({ 
    44                     'snr' : ParameterRange(list(numpy.linspace(0.1,2.0,N_snr))), 
    45                     'kernelseed' : ParameterRange(list([retina.params['kernelseed'] + k for k in range(N_seeds)]))}) 
    46  
    47     # calculates the dimension of the parameter space 
    48     results_dim, results_label = p.parameter_space_dimension_labels() 
    49  
    50     DATA = [] 
    51  
    52     pbar=progressbar.ProgressBar(widgets=[name, " ", progressbar.Percentage(), ' ', 
    53             progressbar.Bar(), ' ', progressbar.ETA()], maxval=N_snr*N_seeds) 
    54     for i_exp, experiment in enumerate(p.iter_inner()): 
    55         params = retina.params 
    56         params.update(experiment) # updates what changed in the dictionary 
    57         # simulate the experiment and get its data 
    58         data = retina.run(params)#,verbose=False) 
    59         # store it 
    60         DATA.append(data)# 
    61         pbar.update(i_exp) 
    62  
    63     results['DATA'] = DATA 
    64  
    65     pbar.finish() 
    66  
    67 ############## PRE-PROCESSING ########################### 
    68 try: 
     62        pbar.finish() 
     63    results.close() 
     64    results = shelve.open('results/mat-pre-' + name) 
     65    ############## PRE-PROCESSING ########################### 
    6966    #boing # uncomment to force recomputing the pre-processing stage 
    7067    lower_edges = results['lower_edges'] 
    …… 
    7370    temporal_OFF = results['temporal_OFF'] 
    7471    map_spatial_ON = results['map_spatial_ON'] 
     72    lower_edges = results['lower_edges'] 
     73    results.close() 
    7574 
    7675except: 
    77  
     76    def temporal_mean(spike_list): 
     77        return numpy.sum(spike_list.firing_rate(t_smooth),axis=0) 
     78         
     79    t_smooth = 100. # ms.  integration time to show fiber activity 
    7880    lower_edges = DATA[0]['out_ON_DATA'].time_axis(t_smooth) 
    79     N_smooth = len(lower_edges) 
    80     params = retina.params 
     81    N_smooth = len(lower_edges)-1 
    8182 
    8283    #N_snr = len(p.snr) 
    …… 
    8687    # 
    8788    N_ret, simtime = params['N_ret'], params['simtime'] 
    88     # 
    8989    x = params['position'][0] 
    9090    y = params['position'][1] 
    …… 
    9292    r = numpy.sqrt(r2) 
    9393    id_center = [int(k) for k in numpy.where( r2 < N_ret**2)[0]] 
    94  
    9594    # mean activity accross kernelseeds as a function of SNR 
    9695    for i_exp, experiment in enumerate(p.iter_inner()): 
    …… 
    9897        index = p.parameter_space_index(experiment) 
    9998        # getting SpikeLists corresponding to the interesting parts (within the center) 
    100         temporal_ON[:,index[1]] += sum(DATA[i_exp]['out_ON_DATA'][id_center].firing_rate(t_smooth))/N_seeds 
    101         temporal_OFF[:,index[1]] += sum(DATA[i_exp]['out_OFF_DATA'][id_center].firing_rate(t_smooth))/N_seeds 
    102         map_spatial_ON[:,index[1]] += DATA[i_exp]['out_ON_DATA'].mean_rates(t_start=simtime/4.,t_stop=3*simtime/4.)/N_seeds 
    103         map_spatial_OFF[:,index[1]] += DATA[i_exp]['out_OFF_DATA'].mean_rates(t_start=simtime/4.,t_stop=3*simtime/4.)/N_seeds 
     99        temporal_ON[:,index[1]] += temporal_mean(DATA[i_exp]['out_ON_DATA'].id_slice(id_center))/N_seeds 
     100        temporal_OFF[:,index[1]] += temporal_mean(DATA[i_exp]['out_ON_DATA'].id_slice(id_center))/N_seeds 
     101        map_spatial_ON[:,index[1]] += DATA[i_exp]['out_ON_DATA'].mean_rates(t_start=simtime/4.,t_stop=3*simtime/4.)#/N_seeds 
     102        map_spatial_OFF[:,index[1]] += DATA[i_exp]['out_OFF_DATA'].mean_rates(t_start=simtime/4.,t_stop=3*simtime/4.)#/N_seeds 
    104103 
     104    results = shelve.open('results/mat-pre-' + name) 
    105105    results['temporal_ON'] = temporal_ON 
    106106    results['map_spatial_OFF'] = map_spatial_OFF 
    107107    results['temporal_OFF'] = temporal_OFF 
    108108    results['map_spatial_ON'] = map_spatial_ON 
     109    results['lower_edges'] = lower_edges 
     110    results.close() 
    109111 
    110112results.close() 
    …… 
    133135#pylab.axes([Lmargin, dmargin , 1.0 - Rmargin- Lmargin,1.0-umargin-dmargin]) # [left, bottom, width, height] 
    134136pylab.subplot(131) 
    135 pylab.scatter(x,y,c=params['amplitude'],faceted=False) 
     137pylab.scatter(x,y,c=params['amplitude'], edgecolors='none') 
    136138 
    137139pylab.subplot(232) 
    138 pylab.plot(lower_edges,temporal_ON) 
     140pylab.plot(lower_edges[:-1],temporal_ON) 
    139141pylab.title('time course (ROI) ',fontsize = 'small') 
    140142#pylab.title('time course ON',fontsize = 'small') 
    141 pylab.xticks( numpy.linspace(0, simtime, 5) ) 
     143pylab.xticks( numpy.linspace(0, params.simtime, 5) ) 
    142144pylab.ylabel('ON activity (spike / s / neuron)') 
    143145#pylab.axis('tight') 
    144146pylab.subplot(235) 
    145 pylab.plot(lower_edges,temporal_OFF) 
     147pylab.plot(lower_edges[:-1],temporal_OFF) 
    146148#pylab.title('time course OFF',fontsize = 'small') 
    147 pylab.xticks( numpy.linspace(0, simtime, 5) ) 
     149pylab.xticks( numpy.linspace(0, params.simtime, 5) ) 
    148150pylab.ylabel('OFF activity (spike / s / neuron)') 
    149151#pylab.axis('tight') 
    150152pylab.xlabel('time (ms)') 
    151153pylab.subplot(233) 
    152 pylab.scatter(x, y, c= map_spatial_ON[:,-1], faceted=False) 
     154pylab.scatter(x, y, c= map_spatial_ON[:,-1], edgecolors='none') 
    153155#pylab.title('spatial distribution ON',fontsize ='small') 
    154156#pylab.title('spatial distribution',fontsize ='small') 
    155157pylab.subplot(236) 
    156 pylab.scatter(x, y, c= map_spatial_OFF[:,-1], faceted=False) 
     158pylab.scatter(x, y, c= map_spatial_OFF[:,-1], edgecolors='none') 
    157159#pylab.title('spatial distribution OFF',fontsize ='small') 
    158160#pylab.xlabel('distance from center') 
    …… 
    161163 
    162164 
    163 pylab.savefig('results/fig-' + name + '.pdf') 
    164 pylab.savefig('results/fig-' + name + '.png', dpi = 300) 
     165if 1: 
     166    pylab.ion() 
     167    pylab.show() 
     168else: 
     169    pylab.savefig('results/fig-' + name + '.pdf') 
     170    pylab.savefig('results/fig-' + name + '.png', dpi = 300) 
     171 
     172 

  • trunk/examples/retina/retina.py

    r326 r340  
    117117        params = {'description' : 'default retina', 
    118118            'N': N, # integer;  square of total number of Ganglion Cells LUP: how do we include types and units in parameters? (or by default it is a float in ISO standards) TODO: go rectangular 
    119             'N_ret': .2, # float;  diameter of Ganglion Cell's RF 
     119            'N_ret': .2, # float;  diameter of Ganglion Cell's RF (max: 1) 
    120120            'K_ret': 4.0, # float; ratio of center vs. surround in DOG 
    121121            'dt'         : 0.1,# discretization step in simulations (ms)