TTB_modellingvortexinteraction/Energy_vortex_vs_bold.m at main · decolab/TTB_modellingvortexinteraction · GitHub

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
clear all
path2=[ '../../Nonequilibrium/'];
addpath(genpath(path2));
path3=[ '../../Tenet/TENET/'];
addpath(genpath(path3));
path4=[ '../../TaskReservoire/DataHCP100ordered'];
addpath(genpath(path4));
path5=[ '../TaskReservoire/DMF_Receptors'];
addpath(genpath(path5));
path6=[ '../../QuantumDiffMap/CHARM'];
addpath(genpath(path6));
path7=[ '../../TaskReservoir/Behaviour'];
addpath(genpath(path7));

load results_turbu_fdiffs.mat;

N=100;
indexN=1:N;

sigma=0.02;
Tmax=10;
dt=0.1;

Isubdiag = find(tril(ones(N),-1));

load results_Ceff_Rest_sub_bold.mat;
Ceff_bold=squeeze(mean(Ceff_rest_sub_bold));

load results_Ceff_Rest_sub_vortex.mat;
Ceff_vortex=squeeze(mean(Ceff_rest_sub_vortex));

%%%  Energy  vortex

[FC,CV,Theta0,J]=hopf_int(Ceff_vortex,f_diff,sigma);
A=-J;
D = 0.5*(sigma^2)*eye(2*N);
iD=inv(D);
Theta=Theta0;
nn=1;
for t=0:dt:Tmax
    Theta=(Theta+Theta')/2;
    Theta=Theta+dt*(-(A*Theta+Theta*A')+2*D);
    DTA=D*inv(Theta)-A;
    [P L]=eig(Theta);
    AP=A*P;
    TT=P'*Theta*P;
    XTX=AP*TT*AP';

    for node=1:2*N
        EntroProd2(nn,node)=DTA(node,node)+ iD(node,node)*XTX(node,node)-A(node,node);
        InfoProcess2(nn,node)=DTA(node,node);
        EntroFlow2(nn,node)=iD(node,node)*XTX(node,node)-A(node,node);
    end
    nn=nn+1;
end

for node=1:2*N
    EntroProd(node)=dt*sum(EntroProd2(:,node));
    InfoProcess(node)=dt*sum(InfoProcess2(:,node));
    EntroFlow(node)=dt*sum(EntroFlow2(:,node));
end
EntroProdTotal_vortex=squeeze(mean(EntroProd));
InfoProcessTotal_vortex=squeeze(mean(InfoProcess));
EntroFlowTotal_vortex=squeeze(mean(EntroFlow));


%%% Energy bold

[FC,CV,Theta0,J]=hopf_int(Ceff_bold,f_diff,sigma);
A=-J;
D = 0.5*(sigma^2)*eye(2*N);
iD=inv(D);
Theta=Theta0;
nn=1;
for t=0:dt:Tmax
    Theta=(Theta+Theta')/2;
    Theta=Theta+dt*(-(A*Theta+Theta*A')+2*D);
    DTA=D*inv(Theta)-A;
    [P L]=eig(Theta);
    AP=A*P;
    TT=P'*Theta*P;
    XTX=AP*TT*AP';

    for node=1:2*N
        EntroProd2(nn,node)=DTA(node,node)+ iD(node,node)*XTX(node,node)-A(node,node);
        InfoProcess2(nn,node)=DTA(node,node);
        EntroFlow2(nn,node)=iD(node,node)*XTX(node,node)-A(node,node);
    end
    nn=nn+1;
end

for node=1:2*N
    EntroProd(node)=dt*sum(EntroProd2(:,node));
    InfoProcess(node)=dt*sum(InfoProcess2(:,node));
    EntroFlow(node)=dt*sum(EntroFlow2(:,node));
end
EntroProdTotal_bold=squeeze(mean(EntroProd));
InfoProcessTotal_bold=squeeze(mean(InfoProcess));
EntroFlowTotal_bold=squeeze(mean(EntroFlow));