Skip to content

Commit 29edd8b

Browse files
fpellegrinifpellegrini
committed
made simulations nsg-cluster ready
1 parent 971df97 commit 29edd8b

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

45 files changed

+2253
-1310
lines changed

‎.DS_Store

-4 KB
Binary file not shown.

‎functions/MI_canolty.m renamed to ‎MI_canolty.m

File renamed without changes.

‎bs_results.mat

5.87 MB
Binary file not shown.

‎eucl.m

+110
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,110 @@
1+
% EUCL: Calculates the euclidean distances among a set of points, or between a
2+
% reference point and a set of points, or among all possible pairs of two
3+
% sets of points, in P dimensions. Returns a single distance for two points.
4+
%
5+
% Syntax: dists = eucl(crds1,crds2)
6+
%
7+
% crds1 = [N1 x P] matrix of point coordinates. If N=1, it is taken to
8+
% be the reference point.
9+
% crds2 = [N2 x P] matrix of point coordinates. If N=1, it is taken to
10+
% be the reference point.
11+
% -----------------------------------------------------------------------
12+
% dists = [N1 x N1] symmetric matrix of pairwise distances (if only crds1
13+
% is specified);
14+
% [N1 x 1] col vector of euclidean distances (if crds1 & ref
15+
% are specified);
16+
% [1 x N2] row vector of euclidean distances (if ref & crds2
17+
% are specified);
18+
% [N1 x N2] rectangular matrix of pairwise distances (if crds1
19+
% & crds2 are specified);
20+
% [1 x 1] scalar (if crds1 is a [2 x P] matrix or ref1 & ref2
21+
% are specified);
22+
%
23+
24+
% RE Strauss, 5/4/94
25+
% 10/28/95 - output row (rather than column) vector for the (reference
26+
% point)-(set of points) case; still outputs column vector for the
27+
% (set of points)-(reference point) case.
28+
% 10/30/95 - for double for-loops, put one matrix-col access in outer loop
29+
% to increase speed.
30+
% 10/12/96 - vectorize inner loop to increase speed.
31+
% 6/12/98 - allow for P=1.
32+
% 11/11/03 - initialize dists to NaN for error return.
33+
34+
function dists = eucl(crds1,crds2)
35+
if nargin == 0, help eucl; return; end;
36+
dists = NaN;
37+
38+
if (nargin < 2) % If only crds1 provided,
39+
[N,P] = size(crds1);
40+
if (N<2)
41+
error(' EUCL: need at least two points');
42+
end;
43+
44+
crds1 = crds1'; % Transpose crds
45+
dists = zeros(N,N); % Calculate pairwise distances
46+
47+
for i = 1:N-1
48+
c1 = crds1(:,i) * ones(1,N-i);
49+
if (P>1)
50+
d = sqrt(sum((c1-crds1(:,(i+1:N))).^2));
51+
else
52+
d = abs(c1-crds1(:,(i+1:N)));
53+
end;
54+
dists(i,(i+1:N)) = d;
55+
dists((i+1:N),i) = d';
56+
end;
57+
if (N==2) % Single distance for two points
58+
dists = dists(1,2);
59+
end;
60+
61+
else % If crds1 & crds2 provided,
62+
[N1,P1] = size(crds1);
63+
[N2,P2] = size(crds2);
64+
if (P1~=P2)
65+
error(' EUCL: sets of coordinates must be of same dimension');
66+
else
67+
P = P1;
68+
end;
69+
70+
crds1 = crds1'; % Transpose crds
71+
crds2 = crds2';
72+
73+
if (N1>1 & N2>1) % If two matrices provided,
74+
dists = zeros(N1,N2); % Calc all pairwise distances between them
75+
for i = 1:N1
76+
c1 = crds1(:,i) * ones(1,N2);
77+
if (P>1)
78+
d = sqrt(sum((c1-crds2).^2));
79+
else
80+
d = abs(c1-crds2);
81+
end;
82+
dists(i,:) = d;
83+
end;
84+
end;
85+
86+
if (N1==1 & N2==1) % If two vectors provided,
87+
dists = sqrt(sum((crds1-crds2).^2)); % Calc scalar
88+
end;
89+
90+
if (N1>1 & N2==1) % If matrix & reference point provided,
91+
crds1 = crds1 - (ones(N1,1)*crds2')'; % Center points on reference point
92+
if (P>1) % Calc euclidean distances in P-space
93+
dists = sqrt(sum(crds1.^2))';
94+
else
95+
dists = abs(crds1)';
96+
end;
97+
end; % Return column vector
98+
99+
if (N1==1 & N2>1) % If reference point & matrix provided,
100+
crds2 = crds2 - (ones(N2,1)*crds1')'; % Center points on reference point
101+
if (P>1) % Calc euclidean distances in P-space
102+
dists = sqrt(sum(crds2.^2));
103+
else
104+
dists = abs(crds2);
105+
end;
106+
end; % Return row vector
107+
end;
108+
109+
return;
110+

‎fp_2chan_sim.m

+98-125
Original file line numberDiff line numberDiff line change
@@ -1,133 +1,106 @@
1-
function fp_2chan_sim
1+
function fp_2chan_sim(seed,iit)
22

3-
fp_addpath_pac
3+
rng('default')
4+
rng(seed)
5+
DIROUT = './';
46

5-
DIRLOG ='/home/bbci/data/haufe/Franziska/log/2chan_sim9/';
6-
if ~exist(DIRLOG); mkdir(DIRLOG); end
7+
logname = num2str(iit);
78

8-
DIROUT = '/home/bbci/data/haufe/Franziska/data/2chan_sim9/';
9-
if ~exist(DIROUT);mkdir(DIROUT); end
9+
%% Parameters
1010

11-
stack_id = str2num(getenv('SGE_TASK_ID'));
12-
rng(stack_id)
11+
N = 120000;
12+
nchan= 2;
1313

14-
iit_ids = (stack_id*6)-5 :(stack_id*6);
14+
%Sampling frequency
15+
fs = 200;
16+
n_trials_s = 60;
17+
n_shuffles = 1000;
1518

19+
%BW
20+
low = [9 11];
21+
high = [58 62];
22+
filt.low = low;
23+
filt.high = high;
1624

17-
for iit = iit_ids
25+
%SNR
26+
snr_v = [0 0.2 0.4 0.6 0.8 1]; %in sim 3
1827

19-
clearvars -except iit iit_ids DIROUT DIRLOG stack_id
20-
21-
logname=num2str(iit);
22-
23-
if ~exist(sprintf('%s%s_work',DIRLOG,logname)) & ~exist(sprintf('%s%s_done',DIRLOG,logname))
24-
eval(sprintf('!touch %s%s_work',DIRLOG,logname))
25-
tic
26-
27-
%% Parameters
28-
29-
% N = 1000000;
30-
N = 120000;
31-
32-
nchan= 2;
33-
34-
%Sampling frequency
35-
fs = 200;
36-
n_trials_s = 60; %before: 500
37-
38-
n_shuffles = 1000;
39-
40-
%BW
41-
low = [9 11];
42-
high = [58 62];
43-
filt.low = low;
44-
filt.high = high;
45-
46-
%SNR
47-
snr_v = [0 0.2 0.4 0.6 0.8 1]; %in sim 3
48-
49-
%% Signal generation
50-
51-
[high_osc, low_osc, pac_0] = syn_sig(N,fs, low, high);
52-
[~,~, pac_1] = syn_sig(N,fs, low, high);
53-
54-
channels_noise = randn(N,2);
55-
channels_noise = channels_noise./ norm(channels_noise(:),'fro');
56-
57-
mixing_matrix = eye(2);
58-
mixing_matrix(1,2)=2*(rand-.5);
59-
mixing_matrix(2,1)=2*(rand-.5);
60-
61-
%% Loop through snrs
62-
63-
for isnr = 1:length(snr_v)
64-
65-
66-
%% Case 1: true bivariate interaction
67-
68-
cse = 1;
69-
fprintf(['SNR ' num2str(isnr) ', case ' num2str(cse) '\n'])
70-
71-
%generate signal
72-
X_1 = [high_osc low_osc];
73-
X_1 = X_1./norm(X_1(:),'fro');
74-
X_1 = snr_v(isnr)*X_1 + (1-snr_v(isnr))*channels_noise;
75-
X_1 = reshape(X_1',nchan,[],n_trials_s);
76-
77-
[pac{cse,isnr}, p{cse,isnr}] = fp_get_all_pac(X_1, fs, filt, n_shuffles);
78-
79-
80-
%% Case 2: true bivariate pac + mixing
81-
82-
cse=2;
83-
fprintf(['SNR ' num2str(isnr) ', case ' num2str(cse) '\n'])
84-
85-
%generate signal with mixing
86-
X_2 = [high_osc low_osc];
87-
X_2 = X_2*mixing_matrix;
88-
X_2 = X_2./norm(X_2(:),'fro');
89-
X_2 = snr_v(isnr)*X_2 + (1-snr_v(isnr))*channels_noise;
90-
X_2 = reshape(X_2',nchan,[],n_trials_s);
91-
92-
[pac{cse,isnr}, p{cse,isnr}] = fp_get_all_pac(X_2, fs, filt, n_shuffles);
93-
94-
%% Case 3: two univariate pac signals
95-
96-
cse = 5;
97-
fprintf(['SNR ' num2str(isnr) ', case ' num2str(cse) '\n'])
98-
99-
%generate signal
100-
X_5 = [pac_0 pac_1];
101-
X_5 = X_5./norm(X_5(:),'fro');
102-
X_5 = snr_v(isnr)*X_5 + (1-snr_v(isnr))*channels_noise;
103-
X_5 = reshape(X_5',nchan,[],n_trials_s);
104-
105-
[pac{cse,isnr}, p{cse,isnr}] = fp_get_all_pac(X_5, fs, filt, n_shuffles);
106-
107-
108-
%% Case 4: two univariate pac signals + mixing
109-
110-
cse = 6;
111-
fprintf(['SNR ' num2str(isnr) ', case ' num2str(cse) '\n'])
112-
113-
%generate signal with mixing
114-
X_6 = [pac_0 pac_1];
115-
X_6 = X_6*mixing_matrix;
116-
X_6 = X_6./norm(X_6(:),'fro');
117-
X_6 = snr_v(isnr)*X_6 + (1-snr_v(isnr))*channels_noise;
118-
X_6 = reshape(X_6',nchan,[],n_trials_s);
119-
120-
[pac{cse,isnr}, p{cse,isnr}] = fp_get_all_pac(X_6, fs, filt, n_shuffles);
121-
122-
end
123-
124-
125-
%% Save
126-
127-
t = toc;
128-
save([DIROUT logname '.mat'],'-v7.3')
129-
save([DIROUT 'pvals_' logname '.mat'],'p','-v7.3')
130-
eval(sprintf('!mv %s%s_work %s%s_done',DIRLOG,logname,DIRLOG,logname))
131-
132-
end
133-
end
28+
%% Signal generation
29+
30+
[high_osc, low_osc, pac_0] = syn_sig(N,fs, low, high);
31+
[~,~, pac_1] = syn_sig(N,fs, low, high);
32+
33+
channels_noise = randn(N,2);
34+
channels_noise = channels_noise./ norm(channels_noise(:),'fro');
35+
36+
mixing_matrix = eye(2);
37+
mixing_matrix(1,2)=2*(rand-.5);
38+
mixing_matrix(2,1)=2*(rand-.5);
39+
40+
%% Loop through snrs
41+
42+
for isnr = 1:length(snr_v)
43+
44+
45+
%% Case 1: true bivariate interaction
46+
47+
cse = 1;
48+
fprintf(['SNR ' num2str(isnr) ', case ' num2str(cse) '\n'])
49+
50+
%generate signal
51+
X_1 = [high_osc low_osc];
52+
X_1 = X_1./norm(X_1(:),'fro');
53+
X_1 = snr_v(isnr)*X_1 + (1-snr_v(isnr))*channels_noise;
54+
X_1 = reshape(X_1',nchan,[],n_trials_s);
55+
56+
[pac{cse,isnr}, p{cse,isnr}] = fp_get_all_pac(X_1, fs, filt, n_shuffles);
57+
58+
59+
%% Case 2: true bivariate pac + mixing
60+
61+
cse=2;
62+
fprintf(['SNR ' num2str(isnr) ', case ' num2str(cse) '\n'])
63+
64+
%generate signal with mixing
65+
X_2 = [high_osc low_osc];
66+
X_2 = X_2*mixing_matrix;
67+
X_2 = X_2./norm(X_2(:),'fro');
68+
X_2 = snr_v(isnr)*X_2 + (1-snr_v(isnr))*channels_noise;
69+
X_2 = reshape(X_2',nchan,[],n_trials_s);
70+
71+
[pac{cse,isnr}, p{cse,isnr}] = fp_get_all_pac(X_2, fs, filt, n_shuffles);
72+
73+
%% Case 3: two univariate pac signals
74+
75+
cse = 5;
76+
fprintf(['SNR ' num2str(isnr) ', case ' num2str(cse) '\n'])
77+
78+
%generate signal
79+
X_5 = [pac_0 pac_1];
80+
X_5 = X_5./norm(X_5(:),'fro');
81+
X_5 = snr_v(isnr)*X_5 + (1-snr_v(isnr))*channels_noise;
82+
X_5 = reshape(X_5',nchan,[],n_trials_s);
83+
84+
[pac{cse,isnr}, p{cse,isnr}] = fp_get_all_pac(X_5, fs, filt, n_shuffles);
85+
86+
87+
%% Case 4: two univariate pac signals + mixing
88+
89+
cse = 6;
90+
fprintf(['SNR ' num2str(isnr) ', case ' num2str(cse) '\n'])
91+
92+
%generate signal with mixing
93+
X_6 = [pac_0 pac_1];
94+
X_6 = X_6*mixing_matrix;
95+
X_6 = X_6./norm(X_6(:),'fro');
96+
X_6 = snr_v(isnr)*X_6 + (1-snr_v(isnr))*channels_noise;
97+
X_6 = reshape(X_6',nchan,[],n_trials_s);
98+
99+
[pac{cse,isnr}, p{cse,isnr}] = fp_get_all_pac(X_6, fs, filt, n_shuffles);
100+
101+
end
102+
103+
104+
%% Save
105+
% save([DIROUT logname '.mat'],'-v7.3')
106+
save([DIROUT 'pvals_' logname '.mat'],'p','-v7.3')

0 commit comments

Comments
 (0)