-
Notifications
You must be signed in to change notification settings - Fork 1
/
p2_perturb_hopf_psilodep.m
76 lines (63 loc) · 2.11 KB
/
p2_perturb_hopf_psilodep.m
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
clear all;
load empiricalLEiDA.mat;
load optimizedhopfawake.mat
P1emp=mean(P1emp);
P2emp=mean(P2emp);
PERTURB=0:0.01:0.1; %% This is synchronisation protocol
%PERTURB=0:-0.025:-0.4; %% This is for the noise Protocol
we=0.09;
ITER=30;
N=90;
a=zeros(N,2);
C=squeeze(Coptim(find(abs(WE-we)<0.0001),:,:));
TSmax=1000;
NSUB=n_Subjects;
TR=2.08; % Repetition Time (seconds)
NumClusters=Number_Clusters;
%%%%%%%%%%%%%%%%%%
omega = repmat(2*pi*f_diff',1,2); omega(:,1) = -omega(:,1);
dt=0.1*TR/2;
Tmax=NSUB*TSmax;
sig=0.02;
dsig = sqrt(dt)*sig; % to avoid sqrt(dt) at each time step
%%%%%%%%%%%%
%% Optimize
%%
wC = we*C;
sumC = repmat(sum(wC,2),1,2); % for sum Cij*xj
for node=1:N/2
node
iwe=1;
for perturb=PERTURB
a=zeros(N,2);
a(node,:)=a(node,:)+perturb;
a(N+1-node,:)=a(N+1-node,:)+perturb;
for iter=1:ITER
xs=zeros(Tmax,N);
z = 0.1*ones(N,2); % --> x = z(:,1), y = z(:,2)
nn=0;
% discard first 3000 time steps
for t=0:dt:3000
suma = wC*z - sumC.*z; % sum(Cij*xi) - sum(Cij)*xj
zz = z(:,end:-1:1); % flipped z, because (x.*x + y.*y)
z = z + dt*(a.*z + zz.*omega - z.*(z.*z+zz.*zz) + suma) + dsig*randn(N,2);
end
% actual modeling (x=BOLD signal (Interpretation), y some other oscillation)
for t=0:dt:((Tmax-1)*TR)
suma = wC*z - sumC.*z; % sum(Cij*xi) - sum(Cij)*xj
zz = z(:,end:-1:1); % flipped z, because (x.*x + y.*y)
z = z + dt*(a.*z + zz.*omega - z.*(z.*z+zz.*zz) + suma) + dsig*randn(N,2);
if abs(mod(t,TR))<0.01
nn=nn+1;
xs(nn,:)=z(:,1)';
end
end
%%%% KL dist between PTR2emp
[PTRsim,Pstates]=LEiDA_fix_cluster(xs',NumClusters,Vemp,TR);
KLps_p(iter)=0.5*(sum(Pstates.*log(Pstates./P2emp))+sum(P2emp.*log(P2emp./Pstates)));
end
KLpstatessleep_perturbed(node,iwe)=mean(KLps_p)
iwe=iwe+1;
end
end
imagesc(KLpstatessleep_perturbed);