-
Notifications
You must be signed in to change notification settings - Fork 0
/
simulations_toys_sc1_grad_grid.py
executable file
·420 lines (353 loc) · 16.2 KB
/
simulations_toys_sc1_grad_grid.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
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
"""
-----------------------------------------------------------------------
Harmoni: a Novel Method for Eliminating Spurious Neuronal Interactions due to the Harmonic Components in Neuronal Data
Mina Jamshidi Idaji, Juanli Zhang, Tilman Stephani, Guido Nolte, Klaus-Robert Mueller, Arno Villringer, Vadim V. Nikulin
https://doi.org/10.1101/2021.10.06.463319
-----------------------------------------------------------------------
script for:
** simulation of toy examples**
In the manuscript figure:
panel A: scenario 1
panel B: scenario 3
panel C: scenario 4
panel D: scenario 5
In this script more other scenarios are included than the 4 that are presented in teh ms.
-----------------------------------------------------------------------
(c) Mina Jamshidi ([email protected]) @ Neurolgy Dept, MPI CBS, 2021
https://github.com/minajamshidi
(c) please cite the above paper in case of using this code for your research
License: MIT License
-----------------------------------------------------------------------
last modified: 20210929 \Mina
-----------------------------------------------------------------------
-----------------------------------------------------------------------
the two-signal scenario without noise
±±±±±±±±±±± ±±±±±±±±±±± ±±±±±±±±±±±
Scenario 1:
------------
x1 -R- x3
| |
y1 -S- y3
scenario 2:
--------------
x1 -R- x3
| |
y1 -S- y3
y2 y4
scenario 3:
--------------
x1 -R- x3
| |
y1 -S- y3
y2 -R- y4
scenario 4:
--------------
x1 x3
| |
y1 y3
y2 -R- y4
scenario 5:
--------------
x1 -R- y4
|
y1 x3
|
y2 y3
scenario 6:
--------------
x1 <--> x3
x1 <--> y4
x1 -R- y4
| ..
y1 ...x3
|
y2 y3
±±±±±±±±±±± ±±±±±±±±±±± ±±±±±±±±±±±
"""
import numpy as np
from numpy import pi
import os.path as op
from matplotlib import pyplot as plt
from scipy.signal import filtfilt, butter
from tools_signal import *
from tools_simulations import data_fun_pink_noise, filtered_randn, produce_nm_phase_locked_sig, adjust_snr
from tools_general import *
from tools_connectivity import *
from scipy import stats
# from harmoni.harmonitools import harmonic_removal_simple
from harmoni.extratools import *
from tools_harmonic_removal import harmonic_removal_simple
# --------------------
# Scenario
# --------------------
scenario = 1 # the scenario to be simulated - pls check the header for the scenario descriptions
# in the following we encode the scenario in the parameters identifying which components exist in the signals
if scenario == 1:
x1_x3_coupling = 1
y2_y4_exist = 0
elif scenario == 2:
x1_x3_coupling = 1
y2_y4_exist = 1
y2_y4_coupling = 0
elif scenario == 3:
x1_x3_coupling = 1
y2_y4_exist = 1
y2_y4_coupling = 1
elif scenario == 4:
x1_x3_coupling = 0
y2_y4_exist = 1
y2_y4_coupling = 1
elif scenario == 5:
x1_x3_coupling = 0
y2_y4_exist = 1
y2_y4_coupling = 0
elif scenario == 6:
x1_x3_coupling = 1
y2_y4_exist = 1
y2_y4_coupling = 0
elif scenario == 7:
x1_x3_coupling = 1
y2_y4_exist = 1
y2_y4_coupling = 0
# --------------------
# general settings
# --------------------
path_save_results = '/data/p_02076/CODES/Codes_CurrentlyWorking/EEG_Networks/build_nets_python36/harmoni/harmoni-supplementary-data/' # fill this in, if you wanna save the results. Otherwise leave it as ''
path_save_fig = '' # fill this in, if you wanna save the figures. Otherwise leave it as ''
# in case you have the seeds for the simulations, fill this in. Otherwise leave it as ''
# path_seeds = ''
path_seeds = ''
# --------------------
# parameters
# --------------------
fs = 256 # sampling frequency
duration = 60 # seconds
n_samples = int(duration * fs) # number of time samples
times = np.arange(0, n_samples)/fs # the time points - used for plotting purpose
max_iter = 50 # number of interactions
c_y2 = 1 # the weight of y2 in the signal
c_y4 = 1 # the weight of y4
noisy = 1 # if the additive noise should be added to the signals. noisy = 1 --> noisy signals
SNR_alpha = dBinv(5) # SNR of the alpha band
SNR_beta = dBinv(-5) # SNR of the beta band
coh = True # to use coherence or PLV as the connectivity measure
# the filter coefficients
b10, a10 = butter(N=2, Wn=np.array([8, 12])/fs*2, btype='bandpass')
b20, a20 = butter(N=2, Wn=np.array([18, 22])/fs*2, btype='bandpass')
# the containers for the optimum values of c and phi
c_abs_opt_1 = np.empty((max_iter,))
c_phi_opt_1 = np.empty((max_iter,))
c_abs_opt_2 = np.empty((max_iter,))
c_phi_opt_2 = np.empty((max_iter,))
# the containers for the synchronization values
synch_sig1x_sig1y = np.empty((max_iter,))
synch_sig1x_yres1_grid = np.empty((max_iter,))
synch_sig1x_yres1_grad = np.empty((max_iter,))
synch_sig2x_sig2y = np.empty((max_iter,))
synch_sig2x_yres2_grid = np.empty((max_iter,))
synch_sig2x_yres2_grad = np.empty((max_iter,))
synch_sig1x_sig2y = np.empty((max_iter,))
synch_sig1x_yres2_grid = np.empty((max_iter,))
synch_sig1x_yres2_grad = np.empty((max_iter,))
synch_sig2x_sig1y = np.empty((max_iter,))
synch_sig2x_yres1_grid = np.empty((max_iter,))
synch_sig2x_yres1_grad = np.empty((max_iter,))
synch_sig1y_sig2y = np.empty((max_iter,))
synch_yres1_yres2_grid = np.empty((max_iter,))
synch_yres1_yres2_grad = np.empty((max_iter,))
if path_seeds == '':
seed = np.random.randint(low=0, high=2 ** 32, size=(max_iter,))
else:
seed = load_pickle(path_seeds)
for n_iter in range(max_iter):
# n_iter = 0
print(n_iter)
np.random.seed(int(seed[n_iter]))
"""
dphi_y1 = 0
dphi_y3 = 0
dphi_x3 = 0
dphi_y4 = 0
"""
dphi_y1 = pi / 2 * np.random.random(1) + pi / 4 # phase-shift of y1 comparing to the phase warped x1
dphi_y3 = pi / 2 * np.random.random(1) + pi / 4 # phase-shift of y3 comparing to the phase of warped x3
dphi_x3 = pi / 2 * np.random.random(1) + pi / 4 # phase-shift of x3 comparing to x1(in case of coupling of x1 & x3)
dphi_y4 = pi / 2 * np.random.random(1) + pi / 4 # phase-shift of y4 comparing to y2
# --------------------------------------------------------------
# generate narrow-band components of sig1 and sig2
# --------------------------------------------------------------
# x1 is the alpha component of sig1 - produced by band-pass filtering random noise
x1 = filtered_randn(8, 12, fs, n_samples)
if x1_x3_coupling: # if sig1 and sig2 are coupled, generate x3 by shifting the phase of x1
x3 = produce_nm_phase_locked_sig(x1, dphi_x3, 1, 1, [8, 12], fs, nonsin_mode=1)
else: # otherwise, also generate x3 by band-pass filtering random noise
x3 = filtered_randn(8, 12, fs, n_samples)
# generate y1 and y3 by phase-warping of x1 and x3, and then adding a phase-shift
y1 = produce_nm_phase_locked_sig(sig=x1, phase_lag=dphi_y1, n=1, m=2, wn_base=[8, 12], sfreq=fs)
y3 = produce_nm_phase_locked_sig(x3, dphi_y3, 1, 2, [8, 12], fs)
# generate a band-pass filtering random noise, it will be saved as y2
y2 = filtered_randn(16, 24, fs, n_samples)
if y2_y4_exist: # if y2 and y4 are contained in sig1 and sig2:
if y2_y4_coupling: # if y2 and y4 are coupled, generate y4 by phase-shifting y2
y4 = produce_nm_phase_locked_sig(y2, dphi_y4, 1, 1, [18, 22], fs, nonsin_mode=1)
else: # otherwise, if y2 and y4 are not coupled:
if scenario == 5 or scenario == 6 or scenario == 7: # if there is a geneuine CFS:
# use phase warping on x1, to generate y4 cross-frequency coupled to x1
y4 = produce_nm_phase_locked_sig(sig=x1, phase_lag=dphi_y4, n=1, m=2, wn_base=[8, 12], sfreq=fs, nonsin_mode=1)
else: # if non of the above cases, generate y4 by band-pass filtering random noise
y4 = filtered_randn(16, 24, fs, n_samples)
# the alpha components of sig1 and sig2 ------------
x_sig1 = x1
x_sig2 = x3
# the beta components of sig1 and sig2 ---------------
if scenario == 7:
y_sig1 = y1
y_sig2 = 0
else:
y_sig1 = y1
y_sig2 = y3
if y2_y4_exist:
if scenario == 7:
y_sig2 = y_sig2 + c_y4 * y4
else:
y_sig1 = y_sig1 + c_y2 * y2
y_sig2 = y_sig2 + c_y4 * y4
# --------------------------------------------------------------
# generate and add the pink noise - SNR is also tuned here
# --------------------------------------------------------------
if noisy:
# generate the noise components ---------
pink_noise_1 = data_fun_pink_noise(times)[np.newaxis, :]
pink_noise_2 = data_fun_pink_noise(times)[np.newaxis, :]
# SNR adjustment ------------
factor_x_sig1 = adjust_snr(np.real(x_sig1), pink_noise_1, SNR_alpha, np.array([8, 12]) / fs * 2)
x_sig1 = x_sig1 / factor_x_sig1
factor_x_sig2 = adjust_snr(np.real(x_sig2), pink_noise_2, SNR_alpha, np.array([8, 12]) / fs * 2)
x_sig2 = x_sig2 / factor_x_sig2
factor_y_sig1 = adjust_snr(np.real(y_sig1), pink_noise_1, SNR_beta, np.array([16, 24]) / fs * 2)
y_sig1 = y_sig1 / factor_y_sig1
factor_y_sig2 = adjust_snr(np.real(y_sig2), pink_noise_2, SNR_beta, np.array([16, 24]) / fs * 2)
y_sig2 = y_sig2 / factor_y_sig2
# final sig1 and sig1 ---------------------------------------
sig1 = np.real(x_sig1 + y_sig1)
sig2 = np.real(x_sig2 + y_sig2)
if noisy: # if noisy add teh pink noise
sig1 += pink_noise_1
sig2 += pink_noise_2
"""
from here on, we pretend that we have the noisy non-sin signal and we wanna use Harmoni to suppress the
harmonic info
"""
# --------------------------------------------------------------
# HARMONI
# --------------------------------------------------------------
# filter sig1 and sig2 in narrow-band
sig1_x = filtfilt(b10, a10, sig1)
sig1_y = filtfilt(b20, a20, sig1)
sig2_x = filtfilt(b10, a10, sig2)
sig2_y = filtfilt(b20, a20, sig2)
# optimization for sig1 and sig2 -------------
y_sig1_res_grad = harmonic_removal_simple(sig1_x, sig1_y, fs)
y_sig1_res_grid = harmonic_removal_simple(sig1_x, sig1_y, fs, method='grid')
y_sig2_res_grad = harmonic_removal_simple(sig2_x, sig2_y, fs)
y_sig2_res_grid = harmonic_removal_simple(sig2_x, sig2_y, fs, method='grid')
# compute the synchronization indices
# we use the absolute coherency as the metric
synch_sig1x_sig1y[n_iter] = compute_phaseconn_with_permtest(sig1_x, sig1_y, 1, 2, fs, plv_type='abs', coh=coh)
synch_sig1x_yres1_grad[n_iter] = compute_phaseconn_with_permtest(sig1_x, y_sig1_res_grad, 1, 2, fs, plv_type='abs', coh=coh)
synch_sig1x_yres1_grid[n_iter] = compute_phaseconn_with_permtest(sig1_x, y_sig1_res_grid, 1, 2, fs, plv_type='abs', coh=coh)
synch_sig2x_sig2y[n_iter] = compute_phaseconn_with_permtest(sig2_x, sig2_y, 1, 2, fs, plv_type='abs', coh=coh)
synch_sig2x_yres2_grad[n_iter] = compute_phaseconn_with_permtest(sig2_x, y_sig2_res_grad, 1, 2, fs, plv_type='abs', coh=coh)
synch_sig2x_yres2_grid[n_iter] = compute_phaseconn_with_permtest(sig2_x, y_sig2_res_grid, 1, 2, fs, plv_type='abs', coh=coh)
synch_sig1x_sig2y[n_iter] = compute_phaseconn_with_permtest(sig1_x, sig2_y, 1, 2, fs, plv_type='abs', coh=coh)
synch_sig1x_yres2_grad[n_iter] = compute_phaseconn_with_permtest(sig1_x, y_sig2_res_grad, 1, 2, fs, plv_type='abs', coh=coh)
synch_sig1x_yres2_grid[n_iter] = compute_phaseconn_with_permtest(sig1_x, y_sig2_res_grid, 1, 2, fs, plv_type='abs', coh=coh)
synch_sig2x_sig1y[n_iter] = compute_phaseconn_with_permtest(sig2_x, sig1_y, 1, 2, fs, plv_type='abs', coh=coh)
synch_sig2x_yres1_grad[n_iter] = compute_phaseconn_with_permtest(sig2_x, y_sig1_res_grad, 1, 2, fs, plv_type='abs', coh=coh)
synch_sig2x_yres1_grid[n_iter] = compute_phaseconn_with_permtest(sig2_x, y_sig1_res_grid, 1, 2, fs, plv_type='abs', coh=coh)
synch_sig1y_sig2y[n_iter] = compute_phaseconn_with_permtest(sig1_y, sig2_y, 1, 1, fs, plv_type='abs', coh=coh)
synch_yres1_yres2_grad[n_iter] = compute_phaseconn_with_permtest(y_sig1_res_grad, y_sig2_res_grad, 1, 1, fs, plv_type='abs', coh=coh)
synch_yres1_yres2_grid[n_iter] = compute_phaseconn_with_permtest(y_sig1_res_grid, y_sig2_res_grid, 1, 1, fs, plv_type='abs', coh=coh)
dict1 = {'seed': seed,
'synch_sig1x_sig1y': synch_sig1x_sig1y, 'synch_sig1x_yres1_grid': synch_sig1x_yres1_grid,
'synch_sig1x_yres1_grad': synch_sig1x_yres1_grad,
'synch_sig2x_sig2y': synch_sig2x_sig2y, 'synch_sig2x_yres2_grad': synch_sig2x_yres2_grad,
'synch_sig2x_yres2_grid': synch_sig2x_yres2_grid,
'synch_sig1x_sig2y': synch_sig1x_sig2y, 'synch_sig1x_yres2_grad': synch_sig1x_yres2_grad,
'synch_sig1x_yres2_grid': synch_sig1x_yres2_grid,
'synch_sig2x_sig1y': synch_sig2x_sig1y, 'synch_sig2x_yres1_grad': synch_sig2x_yres1_grad,
'synch_sig2x_yres1_grid': synch_sig2x_yres1_grid,
'synch_sig1y_sig2y': synch_sig1y_sig2y, 'synch_yres1_yres2_grad': synch_yres1_yres2_grad,
'synch_yres1_yres2_grid': synch_yres1_yres2_grid}
if len(path_save_results):
save_pickle(path_save_results + '/toys_grad_grid_' + 'scenario' + str(scenario), dict1)
#
# # ------------------------------------
# # plotting
# # ------------------------------------
#
# # fig = plt.figure()
# #
# # ax = plt.subplot(231)
# # plot_boxplot_paired(ax, dict1['plv_sig1x_sig1y'], dic1t['plv_sig1x_yres1'], datapoints=True,
# # labels=['plv(s1_x, s1_y)', 'plv(s1_x, s1_y_res)'])
# #
# # ax = plt.subplot(232)
# # plot_boxplot_paired(ax, dict1['plv_sig2x_sig2y'], dict1['plv_sig2x_yres2'], datapoints=True,
# # labels=['plv(s2_x, s2_y)', 'plv(s2_x, s2_y_res)'])
# #
# # ax = plt.subplot(233)
# # plot_boxplot_paired(ax, dict1['plv_sig1x_sig2y'], dict1['plv_sig1x_yres2'], datapoints=True,
# # labels=['plv(s1_x, s2_y)', 'plv(s1_x, s2_y_res)'])
# #
# # ax = plt.subplot(234)
# # plot_boxplot_paired(ax, dict1['plv_sig2x_sig1y'], dict1['plv_sig2x_yres1'], datapoints=True,
# # labels=['plv(s2_x, s1_y)', 'plv(s2_x, s1_y_res)'])
# #
# # ax = plt.subplot(235)
# # plot_boxplot_paired(ax, dict1['plv_sig1y_sig2y'], dict1['plv_yres1_yres2'], datapoints=True,
# # labels=['plv(s1_y, s2_y)', 'plv(s1_y_res, s2_y_res)'])
# #
# # fname_fig = op.join(path_save_fig, 'sc' + str(scenario) + '.eps')
# # fig.savefig(fname_fig, facecolor='white')
#
#
# # ------------------------------------
# # plot by loading your saved results
# # ------------------------------------
# # if you wanna load your saved results. uncomment the follwoing line:
# # dict1 = load_pickle(path_save_results + 'toys_scenario' + str(scenario))
#
# data = (dict1['synch_sig1x_sig1y'][:, np.newaxis], dict1['synch_sig1x_yres1'][:, np.newaxis],
# dict1['synch_sig2x_sig2y'][:, np.newaxis], dict1['synch_sig2x_yres2'][:, np.newaxis],
# dict1['synch_sig1x_sig2y'][:, np.newaxis], dict1['synch_sig1x_yres2'][:, np.newaxis],
# dict1['synch_sig2x_sig1y'][:, np.newaxis], dict1['synch_sig2x_yres1'][:, np.newaxis],
# dict1['synch_sig1y_sig2y'][:, np.newaxis], dict1['synch_yres1_yres2'][:, np.newaxis])
#
# random_coh = random_synchronization_dist(1, 2, duration, f0=10, fs=fs, maxiter=5000)
# perc95, perc99 = np.percentile(random_coh, [95, 99])
#
# fig = plt.figure()
# plt.hlines([perc95, perc99], 0, 11, linestyle='--', color='lightgray')
#
# plt.boxplot(np.concatenate(data, axis=1), notch=True)
# violin_plot([random_coh], positions=[11])
#
# for k in range(max_iter):
# for i1 in range(0, 9, 2):
# plt.plot(np.ones((1, 1)) * (i1+1) + np.random.randn(1, 1) * 0.02, data[i1][k],
# marker='.', color='lightskyblue', markersize=3)
# plt.plot(np.ones((1, 1)) * (i1+2) + np.random.randn(1, 1) * 0.02, data[i1+1][k],
# marker='.', color='lightskyblue', markersize=3)
# x = np.array([i1+1, i1+2])
# y = np.array([data[i1][k], data[i1+1][k]])
# plt.plot(x, y, '-', linewidth=.05)
#
# if len(path_save_fig):
# fname_fig = op.join(path_save_fig, 'sc' + str(scenario) + '.eps')
# fig.savefig(fname_fig, facecolor='white')
#
# # do the stats
# for ii in range(0, 9, 2):
# res = stats.wilcoxon(data[ii].ravel(), data[ii+1].ravel())
# print(res[1], res[1]*5)