https://github.com/macshine/integration
Tip revision: d00534027ab8855f7aec3a32c17f98c6c48431a6 authored by Mac Shine on 14 May 2021, 21:29:02 UTC
Update integration_expt.m
Update integration_expt.m
Tip revision: d005340
hungarian1.m
function [ci_new,F] = hungarian1(ci,beta)
%unique modules per timepoint
[nNodes,nTime] = size(ci);
nDer = nTime - 1;
% dyn_mod = zeros(nMod,nTime);
for t = 1:nTime
temp = tabulate(ci(:,t));
dyn_mod(1:size(temp,1),t) = temp(:,1);
end
%number of modules per window
number_mod = zeros(nTime,1);
for t = 1:nTime
number_mod(t,1) = nnz(dyn_mod(:,t));
end
ignore = double(number_mod>5);
nMod = max(max(ci(:,ignore==0)));
dyn_mod(nMod+1:end,:) = [];
%1-of-k encoding
encode = zeros(nNodes,nMod,nTime);
for t = 1:nTime
if ignore(t,1)==0
C = ci(:,t);
R = 1:numel(C);
A = zeros(numel(C),max(C));
A(sub2ind(size(A),R',C)) = 1;
encode(:,1:size(A,2),t) = A;
else
encode(:,:,t) = 0;
end
end
sprintf('%s','1 of k encoding')
%dice between consecutive time points
%this is usually bimodal
dice_coef_encode = zeros(nMod,nMod,nTime);
for t = 1:nDer
dice_coef_encode(:,:,t) = bsxfun(@corr,encode(:,:,t),encode(:,:,t+1));
end
sprintf('%s','similarity')
%threshold & cost
cost = 1/(double(dice_coef_encode>beta));
%hungarian algorithm
assignment = zeros(nTime,nMod);
for t = 1:nTime
[assignment(t,:),~] = munkres(cost(:,:,t));
end
sprintf('%s','munkres')
%hungarian un-twisting
dyn_mod2 = zeros(nMod,nTime);
for t = 1:nTime
for k = 1:nMod
if number_mod(t,1) < k
dyn_mod2(k,t) = NaN;
end
end
end
dyn_mod2(:,1,:) = dyn_mod(:,1,:); %dyn_mod2 starting with dyn_mod's first assignment
% recoding
tally = max(dyn_mod2(:,1));
for w = 2:nTime-1
for k = 1:nMod
for l = 1:nMod
if dyn_mod2(k,w-1) == 0
dyn_mod2(k,w-1) = tally+1;
tally = tally+1;
end
if assignment(w-1,k) == l
dyn_mod2(l,w) = dyn_mod2(k,w-1);
end
end
end
end
% % number of timepoints that a module exists for
% count_mod = zeros(nMod,1);
%
% for x = 1:max(dyn_mod2(:))
% count_mod(x,1) = countmember(x,dyn_mod2(:));
% end
%
%
% % parcel size of each module over time
% count_mod_size = zeros(nMod,nTime);
%
% for t = 1:nTime
% for y = 1:nMod
% count_mod_size(y,t) = countmember(y,ci(:,t));
% end
% end
%
% for t = 1:nTime
% count_mod_unique(1:size(unique(count_mod_size(:,t))),t) = unique(count_mod_size(:,t));
% end
%
% count_mod_unique(1,:) = [];
% recode the net_thr matrix using new module names (~mucha)
ci_temp = zeros(nNodes,nTime);
for t = 1:nTime
for j = 1:nNodes
for k = 1:nMod
if ci(j,t) == dyn_mod(k,t)
ci_temp(j,t) = dyn_mod2(k,t);
end
end
end
end
sprintf('%s','rename')
%%temporal sorting
%ci_new_name_1_of_k
encode2 = zeros(nNodes,tally,nTime);
for t = 1:nTime
for j = 1:nNodes
for p = 1:tally
if ci_temp(j,t) == p
encode2(j,p,t) = 1;
end
end
end
end
%ci_signature
encode2_sum = nanmean(encode2,3);
encode2_perc = encode2_sum / tally;
%ci_correlation matrix -- can we use a better similarity metric here?
corr_sig = corr(encode2_perc);
%use affinity propogation to cluster
[ap,~,~,~] = apcluster(corr_sig,nanmin(corr_sig(:)));
sprintf('%s','affinity propagation')
%identity of clusters estimated by affinity propogation
unique_ap = unique(ap);
%number of clusters estimated by affinity propogation
count_ap = nnz(unique(ap));
for x = 1:size(unique_ap,1)
freq_ap(x,1) = sum(ap==unique_ap(x));
end
%rename ci_new_name according to ap clusters
dyn_mod3 = dyn_mod2;
for t = 1:nTime
for k = 1:nMod
for p = 1:tally
if dyn_mod2(k,t) == p
dyn_mod3(k,t) = ap(p,1);
end
end
end
end
% rename parcels according to affinity clustering
ci_temp2 = zeros(nNodes,nTime);
ci_new = zeros(nNodes,nTime);
for t = 1:nTime
for j = 1:nNodes
for k = 1:nMod
if ci_temp(j,t) == dyn_mod2(k,t)
ci_temp2(j,t) = dyn_mod3(k,t);
end
end
end
end
%% rename parcels according to order of appearance
for t = 1:nTime
for j = 1:nNodes
for x = 1:size(unique_ap,1)
y = unique_ap(x);
if ci_temp2(j,t) == y
ci_new(j,t) = find(unique_ap==y);
end
end
end
end
F = flexibility(ci_new');
sprintf('%s','flexibility')
end