Changeset 340 for trunk/examples/retina/benchmark_retina.py

Timestamp:

11/14/08 19:17:05 (4 years ago)

Author:

LaurentPerrinet

Message:

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

Files:

• trunk/examples/retina/benchmark_retina.py (modified) (7 diffs)

trunk/examples/retina/benchmark_retina.py

@@ -25 +25 @@
 import os, sys, numpy, pylab, shelve, progressbar
 
+N, N_snr, N_seeds = 1000, 5, 10
+from NeuroTools.parameters import *
+p =  ParameterSpace({
+                'snr' : ParameterRange(list(numpy.linspace(0.1,2.0,N_snr))),
+                'kernelseed' : ParameterRange(list([12345 + k for k in range(N_seeds)]))})
 
 name = sys.argv[0].split('.')[0] # name of the current script withpout the '.py' part
-############## MAKING THE SIMULATIONS ###############
-results = shelve.open('results/mat-' + name)
+
 try:
-    DATA = results['DATA']
-except:
+    ############## MAKING THE SIMULATIONS ###############
+    results = shelve.open('results/mat-' + name)
+    try:
+        DATA = results['DATA']
+        params = results['params']
+    except:
+        from retina import *
+        retina = Retina(N)
+        # calculates the dimension of the parameter space
+        results_dim, results_label = p.parameter_space_dimension_labels()
 
-    from NeuroTools.parameters import *
+        DATA = []
+        pbar=progressbar.ProgressBar(widgets=[name, " ", progressbar.Percentage(), ' ',
+                progressbar.Bar(), ' ', progressbar.ETA()], maxval=N_snr*N_seeds)
+        for i_exp, experiment in enumerate(p.iter_inner()):
+            params = retina.params
+            params.update(experiment) # updates what changed in the dictionary
+            # simulate the experiment and get its data
+            data = retina.run(params)#,verbose=False)
+            # store it
+            DATA.append(data)#
+            pbar.update(i_exp)
 
-    N, N_snr, N_seeds = 1000, 5, 10
-    from retina import *
-    retina = Retina(N)
+        results['DATA'] = DATA
+        results['params'] = retina.params
 
-    t_smooth = 100. # ms.  integration time to show fiber activity
-
-    p =  ParameterSpace({
-                    'snr' : ParameterRange(list(numpy.linspace(0.1,2.0,N_snr))),
-                    'kernelseed' : ParameterRange(list([retina.params['kernelseed'] + k for k in range(N_seeds)]))})
-
-    # calculates the dimension of the parameter space
-    results_dim, results_label = p.parameter_space_dimension_labels()
-
-    DATA = []
-
-    pbar=progressbar.ProgressBar(widgets=[name, " ", progressbar.Percentage(), ' ',
-            progressbar.Bar(), ' ', progressbar.ETA()], maxval=N_snr*N_seeds)
-    for i_exp, experiment in enumerate(p.iter_inner()):
-        params = retina.params
-        params.update(experiment) # updates what changed in the dictionary
-        # simulate the experiment and get its data
-        data = retina.run(params)#,verbose=False)
-        # store it
-        DATA.append(data)#
-        pbar.update(i_exp)
-
-    results['DATA'] = DATA
-
-    pbar.finish()
-
-############## PRE-PROCESSING ###########################
-try:
+        pbar.finish()
+    results.close()
+    results = shelve.open('results/mat-pre-' + name)
+    ############## PRE-PROCESSING ###########################
     #boing # uncomment to force recomputing the pre-processing stage
     lower_edges = results['lower_edges']
@@ -73 +70 @@
     temporal_OFF = results['temporal_OFF']
     map_spatial_ON = results['map_spatial_ON']
+    lower_edges = results['lower_edges']
+    results.close()
 
 except:
-
+    def temporal_mean(spike_list):
+        return numpy.sum(spike_list.firing_rate(t_smooth),axis=0)
+        
+    t_smooth = 100. # ms.  integration time to show fiber activity
     lower_edges = DATA[0]['out_ON_DATA'].time_axis(t_smooth)
-    N_smooth = len(lower_edges)
-    params = retina.params
+    N_smooth = len(lower_edges)-1
 
     #N_snr = len(p.snr)
@@ -86 +87 @@
     #
     N_ret, simtime = params['N_ret'], params['simtime']
-    #
     x = params['position'][0]
     y = params['position'][1]
@@ -92 +92 @@
     r = numpy.sqrt(r2)
     id_center = [int(k) for k in numpy.where( r2 < N_ret**2)[0]]
-
     # mean activity accross kernelseeds as a function of SNR
     for i_exp, experiment in enumerate(p.iter_inner()):
@@ -98 +97 @@
         index = p.parameter_space_index(experiment)
         # getting SpikeLists corresponding to the interesting parts (within the center)
-        temporal_ON[:,index[1]] += sum(DATA[i_exp]['out_ON_DATA'][id_center].firing_rate(t_smooth))/N_seeds
-        temporal_OFF[:,index[1]] += sum(DATA[i_exp]['out_OFF_DATA'][id_center].firing_rate(t_smooth))/N_seeds
-        map_spatial_ON[:,index[1]] += DATA[i_exp]['out_ON_DATA'].mean_rates(t_start=simtime/4.,t_stop=3*simtime/4.)/N_seeds
-        map_spatial_OFF[:,index[1]] += DATA[i_exp]['out_OFF_DATA'].mean_rates(t_start=simtime/4.,t_stop=3*simtime/4.)/N_seeds
+        temporal_ON[:,index[1]] += temporal_mean(DATA[i_exp]['out_ON_DATA'].id_slice(id_center))/N_seeds
+        temporal_OFF[:,index[1]] += temporal_mean(DATA[i_exp]['out_ON_DATA'].id_slice(id_center))/N_seeds
+        map_spatial_ON[:,index[1]] += DATA[i_exp]['out_ON_DATA'].mean_rates(t_start=simtime/4.,t_stop=3*simtime/4.)#/N_seeds
+        map_spatial_OFF[:,index[1]] += DATA[i_exp]['out_OFF_DATA'].mean_rates(t_start=simtime/4.,t_stop=3*simtime/4.)#/N_seeds
 
+    results = shelve.open('results/mat-pre-' + name)
     results['temporal_ON'] = temporal_ON
     results['map_spatial_OFF'] = map_spatial_OFF
     results['temporal_OFF'] = temporal_OFF
     results['map_spatial_ON'] = map_spatial_ON
+    results['lower_edges'] = lower_edges
+    results.close()
 
 results.close()
@@ -133 +135 @@
 #pylab.axes([Lmargin, dmargin , 1.0 - Rmargin- Lmargin,1.0-umargin-dmargin]) # [left, bottom, width, height]
 pylab.subplot(131)
-pylab.scatter(x,y,c=params['amplitude'],faceted=False)
+pylab.scatter(x,y,c=params['amplitude'], edgecolors='none')
 
 pylab.subplot(232)
-pylab.plot(lower_edges,temporal_ON)
+pylab.plot(lower_edges[:-1],temporal_ON)
 pylab.title('time course (ROI) ',fontsize = 'small')
 #pylab.title('time course ON',fontsize = 'small')
-pylab.xticks( numpy.linspace(0, simtime, 5) )
+pylab.xticks( numpy.linspace(0, params.simtime, 5) )
 pylab.ylabel('ON activity (spike / s / neuron)')
 #pylab.axis('tight')
 pylab.subplot(235)
-pylab.plot(lower_edges,temporal_OFF)
+pylab.plot(lower_edges[:-1],temporal_OFF)
 #pylab.title('time course OFF',fontsize = 'small')
-pylab.xticks( numpy.linspace(0, simtime, 5) )
+pylab.xticks( numpy.linspace(0, params.simtime, 5) )
 pylab.ylabel('OFF activity (spike / s / neuron)')
 #pylab.axis('tight')
 pylab.xlabel('time (ms)')
 pylab.subplot(233)
-pylab.scatter(x, y, c= map_spatial_ON[:,-1], faceted=False)
+pylab.scatter(x, y, c= map_spatial_ON[:,-1], edgecolors='none')
 #pylab.title('spatial distribution ON',fontsize ='small')
 #pylab.title('spatial distribution',fontsize ='small')
 pylab.subplot(236)
-pylab.scatter(x, y, c= map_spatial_OFF[:,-1], faceted=False)
+pylab.scatter(x, y, c= map_spatial_OFF[:,-1], edgecolors='none')
 #pylab.title('spatial distribution OFF',fontsize ='small')
 #pylab.xlabel('distance from center')
@@ -161 +163 @@
 
 
-pylab.savefig('results/fig-' + name + '.pdf')
-pylab.savefig('results/fig-' + name + '.png', dpi = 300)
+if 1:
+    pylab.ion()
+    pylab.show()
+else:
+    pylab.savefig('results/fig-' + name + '.pdf')
+    pylab.savefig('results/fig-' + name + '.png', dpi = 300)
+
+