-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathplot_multi_run_from_file.py
196 lines (155 loc) · 8 KB
/
plot_multi_run_from_file.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
# -*- coding: utf-8 -*-
"""
This script provides analysis and plot functions for LFP recordings issued from multiple simulations while varying a single parameter.
Note that a large part of the oscillation analysis (time frequency map, oscillation detection) is based on OpenElectrophy 0.2
The "main" section provides the way to plot results from a simulation output file.
Finally, it also allow to plot few spectograms of LFPs in a simulation output file.
"""
from brian import *
from oscillation_analysis import beta_gamma_detection, spectrum_analysis
import pylab as pl
import pickle
from scipy.stats import sem, nanmean, nanstd
def multi_LFP_analysis(results,osc_th=0.06,freq_cut=40.,burn_time=0.8,distinct_figures=False):
"""
result : list of (param,LFPsig,LFPsr)
From a list of LFPs with repetition of a set of param:
* detect oscillations in each LFP
* compute some oscillation properties
*compute and plot stats of properties across LFPs with the same param
"""
# Analysis of data in results
beta_stats=[]
gamma_stats=[]
params=[]
all_betas=[]
all_gammas=[]
all_LFPlengths=[]
for param,LFPsig,LFPsr in results: # Check each element recorded
params.append(param) # useful to get all simulated parameters
# For that element compute beta and gamma oscillations
list_beta,list_gamma=beta_gamma_detection(LFPsig,LFPsr/Hz,freq_cut=freq_cut,osc_th=osc_th,burn_time=burn_time/second)
all_betas.append(list_beta)
all_gammas.append(list_gamma)
all_LFPlengths.append(1.*LFPsig.size/LFPsr-burn_time)
# Compute oscillation properties from all networks with the same parameter of interest, then compute mean and std of properties
final_list_params=unique(params)
params=array(params)
beta_stats=[]
gamma_stats=[]
for par in final_list_params:
valid_par=where(params==par)[0]
# Prepare to record stats from each recording
list_beta_stats=[]
list_gamma_stats=[]
# Useful function to gather parameters of interest (amp and freq) from all oscillations in a list
def osc_stats(list_osc):
amps,freqs,timeinoscs=[],[],[]
for osc in list_osc:
amps.append(osc.amplitude_max)
freqs.append(osc.freq_max)
dur = (osc.time_stop-osc.time_start)/all_LFPlengths[ind] # proportion of total possible time
timeinoscs.append(dur)
return mean(amps),mean(freqs),sum(timeinoscs),len(amps)
for ind in valid_par:
list_beta_stats.append(osc_stats(all_betas[ind]))
list_gamma_stats.append(osc_stats(all_gammas[ind]))
list_beta_stats=array(list_beta_stats)
list_gamma_stats=array(list_gamma_stats)
# Now compute stats on osc properties for all networks using the same param (use std or sem as wanted)
beta_stats.append([[nanmean(list_beta_stats[:,ind]),nanstd(list_beta_stats[:,ind])] for ind in range(list_beta_stats.shape[1])])
gamma_stats.append([[nanmean(list_gamma_stats[:,ind]),nanstd(list_gamma_stats[:,ind])] for ind in range(list_gamma_stats.shape[1])])
# beta_stats.append([[nanmean(list_beta_stats[:,ind]),sem(list_beta_stats[~isnan(list_beta_stats[:,ind]),ind])] for ind in range(list_beta_stats.shape[1])])
# gamma_stats.append([[nanmean(list_gamma_stats[:,ind]),sem(list_gamma_stats[~isnan(list_gamma_stats[:,ind]),ind])] for ind in range(list_gamma_stats.shape[1])])
# Tranform into arrays for simpler manipulation (shape : param, osc prop, mean/std)
beta_stats=array(beta_stats)
gamma_stats=array(gamma_stats)
lw=2 # linewidth for all curves
if distinct_figures:
figure()
else:
figure(10)
ax=subplot(411) # axis for amp vs param
errorbar(final_list_params,gamma_stats[:,0,0],yerr=gamma_stats[:,0,1],color='k',linewidth=lw)
errorbar(final_list_params,beta_stats[:,0,0],yerr=beta_stats[:,0,1],color='r',linewidth=lw)
ax.set_ylabel("Amplitude")
pl.xlim(final_list_params[0],final_list_params[-1])
pl.ylim(ymin=0.2)
ax=subplot(412) # axis for freq vs param
errorbar(final_list_params,gamma_stats[:,1,0],yerr=gamma_stats[:,1,1],color='k',linewidth=lw)
errorbar(final_list_params,beta_stats[:,1,0],yerr=beta_stats[:,1,1],color='r',linewidth=lw)
pl.ylim(15,90)
pl.xlim(final_list_params[0],final_list_params[-1])
ax.set_ylabel("Frequency")
ax=subplot(413) # total time spent in gamma and in beta
errorbar(final_list_params,gamma_stats[:,2,0],yerr=gamma_stats[:,2,1],color='k',linewidth=lw)
errorbar(final_list_params,beta_stats[:,2,0],yerr=beta_stats[:,2,1],color='r',linewidth=lw)
ax.set_ylabel("Length")
pl.xlim(final_list_params[0],final_list_params[-1])
pl.ylim(ymin=0)
ax=subplot(414) # number of oscillations
errorbar(final_list_params,gamma_stats[:,3,0],yerr=gamma_stats[:,3,1],color='k',linewidth=lw)
errorbar(final_list_params,beta_stats[:,3,0],yerr=beta_stats[:,3,1],color='r',linewidth=lw)
ax.set_ylabel("Number")
pl.xlim(final_list_params[0],final_list_params[-1])
pl.ylim(ymin=-1)
return
def LFP_frequency_content(results,burn_time=0.8*second,freq_min=15.,freq_max=90.,verbose=False):
"""
Compute wavelet map (or FFT ?) of each LFP, to get an average spectral content for the whole LFP
Then compute mean and std (sem ?) of spectral profiles across LFP with the same param
Make a plot with all average profile on the same plot
"""
params=[]
all_spectrums=[]
for param,LFPsig,LFPsr in results: # Check each element recorded
params.append(param) # useful to gather all simulated params
Pxx,freqs=spectrum_analysis(LFPsig,LFPsr/Hz,burn_time=burn_time/second,plot_fig=False,verbose=verbose,return_full=True)
all_spectrums.append(Pxx)
all_spectrums=array(all_spectrums)
params=array(params)
fig=figure()
ax=fig.add_subplot(111)
mask_freq=(freqs>=freq_min)&(freqs<=freq_max)
for ind,param in enumerate(unique(params)):
mask=params==param
ax.plot(freqs[mask_freq],all_spectrums[mask,:][:,mask_freq].mean(axis=0),label=str(param-4.),color=pl.cm.autumn(1.*(ind+1)/(len(unique(params)))),linewidth=2)
ax.legend()
ax.set_xlabel("Frequency (Hz)")
ax.set_ylabel("LFP average spectrum")
def plot_few_spectrograms(results,simul_to_plot=[0],osc_th=0.06,freq_cut=40.,burn_time=0.8*second):
"""
Detect oscillations and plot spectrogram for all simul indices in simul_to_plot
"""
for ind in simul_to_plot:
param,LFPsig,LFPsr=results[ind] # Check each element recorded
print "param: ",param
# For that element compute beta and gamma oscillations
list_beta,list_gamma=beta_gamma_detection(LFPsig,LFPsr/Hz,freq_cut=freq_cut,osc_th=osc_th,burn_time=burn_time/second,plot_fig=True)
return
if __name__=="__main__":
distinct_figures=False
# Files to plot
list_filenames=["put_filename_root_here"]
list_filenames=["test_dict_connectproba"]
# Detection osc params
osc_th=[0.2, 0.2] # oscillation threshold
freq_cut=40.
burn_time=0.8*second
# Spectrum params
freq_min=15.
freq_max=90.
for filename in list_filenames:
# Load data
fid=open(filename+"_multi.out","rb")
param_dict,results=pickle.load(fid)
fid.close()
print "Simulation params of plotted file: ",filename
for key,val in param_dict.items():
print key,': ',val
print
tmp_results=[res[:3] for res in results] # To keep only LFP and param
multi_LFP_analysis(tmp_results,osc_th=osc_th,freq_cut=freq_cut,burn_time=burn_time,distinct_figures=distinct_figures)
LFP_frequency_content(tmp_results,burn_time=burn_time,freq_min=freq_min,freq_max=freq_max)
#~ plot_few_spectrograms(tmp_results,simul_to_plot=range(6),osc_th=osc_th,freq_cut=freq_cut,burn_time=burn_time)
show()