function [timeDelay_per_group,timeDelay_per_groupSh] = time_delay_groupsPairs_visNet_v4(timeDelayMap, timeDelayMapSh, iLab, bandLab, computedROIs, sort_idxs, params)
% Average time delay maps over groups of emitter and receiver ROIs
% v3 + bootstap shuff

corrErrLabels = {'all','corr','err'};

for ceLab = corrErrLabels
    tPoints = size(timeDelayMap.(ceLab{1}).(iLab).(bandLab),3);
    dPoints = numel(params.minSelDelay:params.maxSelDelay);

    % sort time-delay maps
    tdMap = timeDelayMap.(ceLab{1}).(iLab).(bandLab)(sort_idxs,sort_idxs,:,:); % not taking params.minSelDelay:params.maxSelDelay because TD maps are already precomputed in that interval
    if strcmp(ceLab{1},'all')
        tdMapSh = timeDelayMapSh.(ceLab{1}).(iLab).(bandLab)(sort_idxs,sort_idxs,:,:,:);
    end
    
    % All to all ROIs
    roisPairLab = 'all_all';
    all_idxs = logical(ones(1,numel(sort_idxs)));
    
    % Pick maps for selected ROIs
    tmpMap = tdMap(all_idxs,all_idxs,:,:);
    if strcmp(ceLab{1},'all')
        tmpMapSh = tdMapSh(all_idxs,all_idxs,:,:,:);
    end
    
    nTmpPairs = (sum(all_idxs)*(sum(all_idxs)-1))/2;
    timeDelay_per_group.(ceLab{1}).(roisPairLab).forward = nan(nTmpPairs,tPoints,dPoints);
    timeDelay_per_group.(ceLab{1}).(roisPairLab).feedback = nan(nTmpPairs,tPoints,dPoints);
    if strcmp(ceLab{1},'all')
        timeDelay_per_groupSh.(ceLab{1}).(roisPairLab).forward = nan(nTmpPairs,tPoints,dPoints,params.nBoot);
        timeDelay_per_groupSh.(ceLab{1}).(roisPairLab).feedback = nan(nTmpPairs,tPoints,dPoints,params.nBoot);
    end
    
    pairCount = 1;
    for relatEmitIdx = 1:sum(all_idxs)
        for relatRecIdx = relatEmitIdx+1:sum(all_idxs)
            timeDelay_per_group.(ceLab{1}).(roisPairLab).forward(pairCount,:,:) = tmpMap(relatEmitIdx,relatRecIdx,:,:);
            timeDelay_per_group.(ceLab{1}).(roisPairLab).feedback(pairCount,:,:) = tmpMap(relatRecIdx,relatEmitIdx,:,:);
            if strcmp(ceLab{1},'all')
                for shIdx = 1:params.nBoot
                   timeDelay_per_groupSh.(ceLab{1}).(roisPairLab).forward(pairCount,:,:,shIdx) = tmpMapSh(relatEmitIdx,relatRecIdx,:,:,shIdx);
                   timeDelay_per_groupSh.(ceLab{1}).(roisPairLab).feedback(pairCount,:,:,shIdx) = tmpMapSh(relatRecIdx,relatEmitIdx,:,:,shIdx);
                end
            end
            pairCount = pairCount+1;
        end
    end
    timeDelay_per_group.(ceLab{1}).(roisPairLab).forward = squeeze(mean(timeDelay_per_group.(ceLab{1}).(roisPairLab).forward,1));
    timeDelay_per_group.(ceLab{1}).(roisPairLab).feedback = squeeze(mean(timeDelay_per_group.(ceLab{1}).(roisPairLab).feedback,1));
    if strcmp(ceLab{1},'all')
        timeDelay_per_groupSh.(ceLab{1}).(roisPairLab).forward = squeeze(mean(timeDelay_per_groupSh.(ceLab{1}).(roisPairLab).forward,1));
        timeDelay_per_groupSh.(ceLab{1}).(roisPairLab).feedback = squeeze(mean(timeDelay_per_groupSh.(ceLab{1}).(roisPairLab).feedback,1));
    end
    
    % groups pairs
    for geIdx = 1:numel(computedROIs) % emitter group
        emit_idxs = strcmp(computedROIs(sort_idxs),computedROIs{geIdx});
        for grIdx = geIdx:numel(computedROIs) % receiver group
            roisPairLab = [computedROIs{geIdx},'_',computedROIs{grIdx}];
            rec_idxs = strcmp(computedROIs(sort_idxs),computedROIs{grIdx});

            all_idxs = emit_idxs | rec_idxs;
            % Squared array with first 2 dimensions equal to communication only
            % within the selected subgroups pair
            tmpMap = tdMap(all_idxs,all_idxs,:,:);
            if strcmp(ceLab{1},'all')
                tmpMapSh = tdMapSh(all_idxs,all_idxs,:,:,:);
            end
            
            nTmpPairs = (sum(all_idxs)*(sum(all_idxs)-1))/2;
            timeDelay_per_group.(ceLab{1}).(roisPairLab).forward = nan(nTmpPairs,tPoints,dPoints);
            timeDelay_per_group.(ceLab{1}).(roisPairLab).feedback = nan(nTmpPairs,tPoints,dPoints);
            if strcmp(ceLab{1},'all')
                timeDelay_per_groupSh.(ceLab{1}).(roisPairLab).forward = nan(nTmpPairs,tPoints,dPoints,params.nBoot);
                timeDelay_per_groupSh.(ceLab{1}).(roisPairLab).feedback = nan(nTmpPairs,tPoints,dPoints,params.nBoot);
            end
            
            pairCount = 1;
            for relatEmitIdx = 1:sum(all_idxs)
                for relatRecIdx = relatEmitIdx+1:sum(all_idxs)
                    timeDelay_per_group.(ceLab{1}).(roisPairLab).forward(pairCount,:,:) = tmpMap(relatEmitIdx,relatRecIdx,:,:);
                    timeDelay_per_group.(ceLab{1}).(roisPairLab).feedback(pairCount,:,:) = tmpMap(relatRecIdx,relatEmitIdx,:,:);
                    if strcmp(ceLab{1},'all')
                        for shIdx = 1:params.nBoot
                           timeDelay_per_groupSh.(ceLab{1}).(roisPairLab).forward(pairCount,:,:,shIdx) = tmpMapSh(relatEmitIdx,relatRecIdx,:,:,shIdx);
                           timeDelay_per_groupSh.(ceLab{1}).(roisPairLab).feedback(pairCount,:,:,shIdx) = tmpMapSh(relatRecIdx,relatEmitIdx,:,:,shIdx);
                        end
                    end
                    pairCount = pairCount+1;
                end
            end
            timeDelay_per_group.(ceLab{1}).(roisPairLab).forward = squeeze(mean(timeDelay_per_group.(ceLab{1}).(roisPairLab).forward,1));
            timeDelay_per_group.(ceLab{1}).(roisPairLab).feedback = squeeze(mean(timeDelay_per_group.(ceLab{1}).(roisPairLab).feedback,1));
            if strcmp(ceLab{1},'all')
                timeDelay_per_groupSh.(ceLab{1}).(roisPairLab).forward = squeeze(mean(timeDelay_per_groupSh.(ceLab{1}).(roisPairLab).forward,1));
                timeDelay_per_groupSh.(ceLab{1}).(roisPairLab).feedback = squeeze(mean(timeDelay_per_groupSh.(ceLab{1}).(roisPairLab).feedback,1));
            end
        end
    end
end

end