PhaseClustering/FindCom_permute.m at main · TemporalBrain/PhaseClustering · 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
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
function [finaldiv,com, iter]=FindCom_permute(temp,com,iter,N,finaldiv,mini,phasf, perm_nr,opp)

% Author: Marika Strindberg 2025-2026, marika.strindberg@ki.se

% What the code does: optimizes for communities with strong integration with in communities and segregation between commuities.

% temp = network matrix with weights cos(delta0) (at time t)
% phasf = vector with smallest instantaneous phase value, range [-pi,pi]
% iter = which iterative step
% N = number of nodes in the network
% finaldiv = constant which updates to one when all nodes are assigned to a community
% mini = thr1= maximum phase difference (cos) within community, set by user

% permn_nr = number of permutations in order to find the largest and most
% phase consitent community


assigned=find(com~= 0); % assigned nodes
un_assigned=find(com==0); % not assigned nodes


phasfzero=phasf; %THIS VERSION uses smallest negative
phasfzero(assigned)=10; %% ADdded so that already assigned areas are not reused.

temp_remain=temp;
temp_remain(assigned,:)=10;
temp_remain(:,assigned)=10;

% r1 is the seed for the first community, it is the node with most
% negative instantaneous phase value, range [-pi,pi]
% r2 is the seed for the second community that is most segregated to r1

r1=find(phasfzero==min(phasfzero(:)));

tempR1=temp_remain(r1,:); %Vector with the weights between r1 and rest of unassigned nodes

if min(tempR1)< opp %


    r2=find(tempR1==min(tempR1)); % the node most segregated with r1

    r3=setdiff(un_assigned,[r1,r2]); % remaining nodes other than r1 and r2

    r1=AddCom_permute(r1,r3,temp,mini,opp,r2,phasf,perm_nr); % first grow the r1 community
                            %r1,r3,temp,mini,opp,r2,phasf,perm_nr
    r3=setdiff(un_assigned,[r1,r2]); % remaning nodes not in r1 community or r2

    r2=AddCom_permute(r2,r3,temp,mini,opp,r1,phasf,perm_nr); % then grow the r2 community


else   %IF NO SegREATION ACCORDING TO thr opp no com 2


    r2=[];
    r3=setdiff(un_assigned,[r1,r2]);
    r1=AddComSingle_permute(r1,r3,temp,mini,phasf,perm_nr);


end


% Check that r1 and r2 are same dim
[r1s1, r1s2]=size(r1);
if r1s1>r1s2

    r1=r1';
end


[r2s1, r2s2]=size(r2);
if r2s1>r2s2


    r2=r2';
end

% Find potential third community
allN=[1:N]';
rtemp = [r1,r2,assigned'];
r3=setdiff(allN,rtemp); % the areas that are not part of temporary com 1 and 2

% %% If a third community is present
% % see if can be added to an existing community
%
if numel(r3) > 0

    %First r1
    if isempty(r1)==0

        r1=AddCom_permute(r1,r3,temp,mini,opp,r2,phasf,perm_nr);
        rtemp = [r1,r2,assigned'];
        r3=setdiff(allN,rtemp); % the areas that are not part of temporary com 1 and 2
    end
    if isempty(r2)==0

        r2=AddCom_permute(r2,r3,temp,mini,opp,r1,phasf,perm_nr); % returns the rows of the cluster
    end

else % if all areas are assigned then stop

end


com1=iter+1;
com2=iter+2;
com(r1)=iter+1; %iter changes with 2 each step

if isempty(r2)==0

    com(r2)=iter+2;
else

    com2=0;
end


if sum(com==0)==0
    finaldiv=0;

end


r1=[]; r2=[]; r3=[];
com1=[]; com2=[];


end