2023_Blaschke_Stroke_tDCS_rsfMRI/scripts/addendum/rsfmri_metrics.m

function [mNetX,mNet_Trim,SW,CC,CPL] = rsfmri_metrics(TCfile,GSR_signal,AUC,cut,scratch,want2c,GSR)

[matrix] = load(fullfile(TCfile.folder,TCfile.name));


    funcmat         = double(matrix.matrix);
    ts.ts           = funcmat;
    ts.tr           = 2.48;   
    ts.Nnodes       = size(funcmat,2);
    ts.Ntimepoints  = size(funcmat,1);
    ts.Nsubjects    = 1;
    ts.NnodesOrig   = ts.Nnodes;
    ts.NtimepointsPerSubject = size(funcmat,1);
    ts.DD           = 1:ts.Nnodes;
    ts.UNK          = [];   

    ts_scrub = ts;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        if GSR==1
%             for rix = 1:ts.Nnodes
%                 [~,~,Resid] = regress(funcmat(:,rix),GSR_signal);
%                 funcmat_GSR(:,rix) = Resid';
%             end
            
            yadj = zeros(size(funcmat));
            X = detrend(GSR_signal,'constant');
            for v = 1:ts.Nnodes
%                 b = regression(funcmat(:,v),X);
                    [Q, R]=qr(X,0);
                    b = R\(Q'*funcmat(:,v));
                yvoxmodel = b(1)*X(:,1);

                % Adjust time-series
               yadj(:,v) = funcmat(:,v) - yvoxmodel;
            end
  


            ts_scrub.ts = double(yadj);
        end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        if cut==1
        funcmat_scrub = ts_scrub.ts;
        funcmat_scrub(scratch,:) = [];
        ts_scrub.ts = funcmat_scrub;
        ts_scrub.Ntimepoints  = size(funcmat_scrub,1);
        ts_scrub.NtimepointsPerSubject = size(funcmat_scrub,1);
        end

    

netmats       = nets_netmats_md(ts,1,0,'corr'); % correlation transformed 2 Z-scores
mNet          = reshape(netmats(1,:),[sqrt(length(netmats)) sqrt(length(netmats))]);
 
netmats_scrub       = nets_netmats_md(ts_scrub,1,0,'corr'); % correlation transformed 2 Z-scores
mNet_scrub          = reshape(netmats_scrub(1,:),[sqrt(length(netmats_scrub)) sqrt(length(netmats_scrub))]);
 
mNetX = mNet_scrub;

mNetTresh = zeros(98,98,AUC(end));
CPLx = zeros(1,AUC(end));
CCx = zeros(1,AUC(end));
NCC = zeros(1,AUC(end));
NCPL = zeros(1,AUC(end));
SWx = zeros(1,AUC(end));

f = waitbar(0,strcat('0/100'));
for gi=AUC
        p=gi/100;
        waitbar(1i/100,f,num2str(gi),'/100');
        mNetD              = threshold_proportional(mNetX,p);
        mNetTresh(:,:,gi)  = mNetD;
%         mNetD(mNetD==0)   = nan;
        mNetDWght   = weight_conversion(mNetD,'normalize');
        mNetDBin   = weight_conversion(mNetD,'binarize');

%         Eglob (gi)          = efficiency_wei(mNetD);         % global efficieny
        D                   = distance_wei_floyd(mNetDWght,'inv');
        Dx                   = distance_wei_floyd(mNetDBin,'inv');
        
        lambda              = charpath(D,0,0);
        CPLx(gi)            = lambda;
        cplb                = charpath(Dx,0,0);
        clx                 = clustering_coef_wu(mNetDWght);
        CCx(gi)             = squeeze(mean(clx));
        clb                 = clustering_coef_wu(mNetDBin);

                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    % small worldness density threshold
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                      
                    %   create random network to asses small worldness
                    lambda_rdm = zeros(1,100);
                    clcf_rdm = zeros(1,100);
                    for k=1:100
                        random              = makerandCIJ_und(size(mNetD,1),nnz(mNetD)/2); % compute random network with same number of nodes and edges
                        D_rdm               = distance_bin(random);
                        lambda_rdm(k)       = charpath(D_rdm,0,0);
                        clcf_rdm(k)         = mean(clustering_coef_bu(random));
                    end

                     % small-worldness 
                     norm_clcf      = abs(squeeze(mean(clb)))/abs(mean(clcf_rdm));
                     NCC(gi)         = norm_clcf;
                     norm_lambda    = (cplb/(mean(lambda_rdm)));
                     NCPL(gi)        = norm_lambda;
              SWx(gi) = norm_clcf/norm_lambda;
                     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
        mNetTresh(mNetTresh==0)=nan;
        mNet_Trim = squeeze(mean(mNetTresh,3,'omitnan'));

%         for aucx = 1:length(AUC)
%             auc_value(auc) = 
            CC = trapz (AUC, CCx(AUC))/length(AUC);
            CPL = trapz (AUC, CPLx(AUC))/length(AUC);
            SW = trapz (AUC, SWx(AUC))/length(AUC);
            
    F = findall(0,'type','figure','tag','TMWWaitbar');
    delete(F)
    
    %%
    if want2c ==1
    figure('units','normalized','outerposition',[0 0 1 1])
         
        subplot(131)
        plotmat_values_rsfMRI(mNet,[],[-1 1],0,[]);
        hold on
        axis square;
        plot(get(gca,'xlim'),[mean(get(gca,'ylim')), mean(get(gca,'ylim'))],'k-','linewidth',2);
        plot([mean(get(gca,'xlim')), mean(get(gca,'xlim'))], get(gca,'ylim'),'k-','linewidth',2);
        title('raw mNet')
% %         ntx = threshold_proportional(mNetX,median(AUC)/100);
% %         ntx(ntx==0)=nan;
% %         hx = histogram(ntx,'BinEdges',linspace(min(min(ntx)),max(max(ntx)),20),'Normalization','probability');

        subplot(132)
        plotmat_values_rsfMRI(mNet_scrub,[],[-1 1],0,[]);
        hold on
        axis square;
        plot(get(gca,'xlim'),[mean(get(gca,'ylim')), mean(get(gca,'ylim'))],'k-','linewidth',2);
        plot([mean(get(gca,'xlim')), mean(get(gca,'xlim'))], get(gca,'ylim'),'k-','linewidth',2);
        title('scratched mNet')

        mNetDelta = mNet - mNet_scrub;
        subplot(133)
        plotmat_values_rsfMRI(mNetDelta,[],[-.5 .5],0,[]);
        hold on
        axis square;
        plot(get(gca,'xlim'),[mean(get(gca,'ylim')), mean(get(gca,'ylim'))],'k-','linewidth',2);
        plot([mean(get(gca,'xlim')), mean(get(gca,'xlim'))], get(gca,'ylim'),'k-','linewidth',2);
        title('delta')
% 
%         figure('units','normalized','outerposition',[0 0 1 1])
%         tpx = linspace(10,90,5)
%         for histi = 1:5
%             ntx = threshold_proportional(mNetX,tpx(histi)/100);
%             ntx(ntx==0)=nan;
%             wbin = linspace(min(min(ntx)),max(max(ntx)),20);
%             subplot(1,5,histi)
%             hy = histogram(ntx,'Normalization','probability','BinEdges',wbin);
%             title(num2str(tpx(histi)))
%             ylim([0 0.2])
%         end
%          
% %         
% pause(3)
    end
% close all
%%

end
  

     
%         h = histogram(mNet,'Normalization','probability')