% do within-stage PLI surrogates after geting bump locations and define tempral windows for connectivity

clear all
close all
clc

nSurro  = 200;
buOff   = 2; % FC stges everything between bumps except for bumps
withBL  = true;
dimBL   = [10, 6]; % form 40 - 10 to onset - 6 samples -> BL with 25 samples, 250 millisecons
dataHSMM = '/hmmBu4RroiAveMap3StagBL.mat';
savename = '/surroCirchmmBu4RroiAveMap3StagBL.mat';


blOff   = 40; % 400 milliseconds baseline
nSj     = 18;
fs      = 100;
nRoi    = 68;
nFC     = (nRoi^2)/2 - nRoi/2;

savePath = './dataSave';
dataPath = './data';
dataMEG  = './data/thetaROIs_R.mat';
load([dataPath dataHSMM], 'proBump', 'x','y')
load(dataMEG, 'roiaveR','xR', 'yR', 'subjtR')

nStage  = size(proBump, 3)+1;
nStageBL= nStage;
buOff_i = buOff * ones(1, nStage-1);
buOff_i = [0, buOff_i];
buOff_e = buOff * ones(1, nStage-1);
buOff_e = [buOff_e, 0];
if withBL, 
    nStageBL = nStage + 1; 
    buOff_i  = [dimBL(1), buOff_i];
    buOff_e  = [dimBL(2), buOff_e];
end

roiR = single(roiaveR);
xd   = xR;
yd  = yR;
subj = subjtR;

clear roiaveR xR yR subjtR 

pairs =  combnk([1:nRoi],2);
tooShort = [];
pliSurro   = zeros(nStageBL, nFC, nSurro, 'single');

for sj = 1:nSj
    xj  = x(subj==sj); % bumps
    yj  = y(subj==sj);
    xdj = xd(subj==sj); % data 
    ydj = yd(subj==sj);
    nTr     = length(xj);
    lens    = yj - xj + 1;
    pliTr   = zeros(nStageBL, nFC, nSurro, nTr,'single');
    maxN    = max(lens);
    buSampl = ones(nTr, nStage,'single');
    % location of bumps given by HsMM
    for bu = 1:nStage-1
        buSampl(:,bu+1)   = round( [1:size(proBump,1)] * squeeze(proBump(:, subj==sj, bu)));
    end
    buSampl = horzcat(buSampl, lens);
    buSampl = buSampl + blOff;
    if withBL,
        buSampl = horzcat(ones(size(buSampl,1), 1), buSampl);
    end
    for tr = 1:nTr
        % pli: corss spectral power
        data    = hilbert(roiR(xdj(tr):ydj(tr),:));
        data    = repmat(data,1,1,nSurro);
        dataSu  = zeros(size(data));
        raIx    = uint32(randi(lens(tr)-1,1,nSurro));
        for si = 1:nSurro
            dataSu(:,:,si) = circshift(data(:,:,1), raIx(si));
        end
        % corss spectral density
        csd = data(:, pairs(:,1),:) .* conj(dataSu(:, pairs(:,2),:));
        % sample phase differences sign
        siCsd = sign( imag( csd )); % [samples, FC, surros]
         % stage mean PLI
        for s = 1:size(buSampl,2)-1  
            if (buSampl(tr,s+1)-buOff_e(s)) - (buSampl(tr,s)+buOff_i(s)) <= 0,
                pliTr(s,:,:,tr) = 0; % [stages, FC, surro, trials]
                tooShort = vertcat(tooShort, [s,tr]);
            end    
            pliTr(s,:,:,tr) = mean( siCsd( buSampl(tr,s)+buOff_i(s):buSampl(tr,s+1)-buOff_e(s) ,:,:), 1);
        end
    end
    if any(any(any(any(isnan(pliTr)))))
         aNaN = 1;
    end
    pliSurro(:,:,:) =  squeeze(pliSurro(:,:,:)) + squeeze( mean(pliTr, 4) ) ./ nSj;  % mean over subjects
end
p95a = prctile(abs(pliSurro),95,3); % minimum FC respects to Surrogates
p98a = prctile(abs(pliSurro),98,3); % minimum FC respects to Surrogates

%save([savePath savename], 'pliSurro','p95a', 'p98a', 'tooShort')