% Script used to select the two sets of alternative senders used in the
% cFIT analysis of the EEG dataset (Fig. S9 and S12)

% Find max DI and delays between neighbous electrodes
load metadata
load allEEGchanLocs
load original_ROT_LOT_electrodes
bins = 2;

% Define all electrodes in Front Lobe (for the weak alternative senders
% selection) and in LOT and ROT (for the strong alternative senders
% selection)
allROTelectrodes = {'B2','B3','B4','B5','B6','B7','B8','B9','B10','B11',...
    'B12','B13','B14','B15','B16','B17','B18','B19'};
allLOTelectrodes = {'A6','A7','A8','A9','A10','A11','A12','D16','D17',...
    'D24','D25','D26','D27','D28','D29','D30','D31','D32'};
allFrontElectrodes = {'A1','B1','B20','B21','B22','B23',...
    'B24','B25','B26','B27','B28','B29','B30','B31','B32'...
    'C1','C2','C3','C4','C5','C6','C7','C8','C9','C10','C11',...
    'C12','C13','C14','C15','C16','C17','C18','C19','C20','C21','C22','C23',...
    'C24','C25','C26','C27','C28','C29','C30','C31','C32'...
    'D1','D2','D3','D4','D5','D6','D7','D8','D9','D10','D11',...
    'D12','D13','D14','D15','D18','D19','D20','D21','D22','D23'};

ROTelectrodes = cell(1,Ns);
LOTelectrodes = cell(1,Ns);
FrontElectrodes = cell(1,Ns);

maxElectrodeROT = zeros(1,Ns);
maxDelayROT = zeros(1,Ns);
maxElectrodeLOT = zeros(1,Ns);
maxDelayLOT = zeros(1,Ns);
%% Compute DI over time and delays between the origianl and the ROT
% emisphere, pooling times together / averaging over times

delayMax = 20;
tMin = 151;
tMax = 351;
ROTTEmap = cell(1,Ns);
LOTTEmap = cell(1,Ns);

%% Find electrodes in ROT and LOT that are not the original ones
for subj = 1:Ns
    
    disp(['Compute subject ', num2str(subj)])
    
    opts = [];
    opts.method = 'dr';
    opts.bias = 'naive';
    opts.trperm = 0;
    nElectrodes = size(allChanLocs{subj},2);
    allChanLocs{subj} = allChanLocs{subj};
    ROTelectrodes{subj} = allChanLocs{subj}(1,originalMaxROT(subj)).labels; 
    LOTelectrodes{subj} = allChanLocs{subj}(1,originalMaxLOT(subj)).labels; 
    
    electrodesLabel = cell(1,nElectrodes);
    for eidx = 1:nElectrodes
        electrodesLabel{1,eidx} = cell(1,1);
        electrodesLabel{1,eidx}{1} = allChanLocs{subj}(1,eidx).labels;
    end
    % Find ROT and LOT electrodes
    ROTemitterNeighbours{subj} = [];
    LOTemitterNeighbours{subj} = [];
    FrontElectrodes{subj} = [];

    for nidx = 1:size(allROTelectrodes,2)
        if ~strcmp(allROTelectrodes{nidx} , ROTelectrodes{subj})
            findIdx = find(strcmp([electrodesLabel{:}], allROTelectrodes{nidx}));
            if ~isempty(findIdx)
                ROTemitterNeighbours{subj} = [ROTemitterNeighbours{subj}, findIdx];
            end
        end
    end
    
     for nidx = 1:size(allLOTelectrodes,2)
        if ~strcmp(allLOTelectrodes{nidx} , LOTelectrodes{subj})
            findIdx = find(strcmp([electrodesLabel{:}], allLOTelectrodes{nidx}));
            if ~isempty(findIdx)
                LOTemitterNeighbours{subj} = [LOTemitterNeighbours{subj}, findIdx];
            end
        end
     end
     
     for nidx = 1:size(allFrontElectrodes,2)
        findIdx = find(strcmp([electrodesLabel{:}], allFrontElectrodes{nidx}));
        if ~isempty(findIdx)
            FrontElectrodes{subj} = [FrontElectrodes{subj}, findIdx];
        end
    end
end

%% Maximum informative electrodes
load ROT_LOT_electrodes_per_subject
load allSubjects_stimulus_info_maps_2bins_131120
load original_ROT_LOT_electrodes

time_range = 151:351; %i.e. 100-250ms post-stimulus (100-200ms in Roussellet papers)
x_time = (100:2:200);
K = 15; % number of maxima
selected_max_electrode.LE = zeros(K,Ns); valsLE = zeros(K,Ns);
selected_max_electrode.RE = zeros(K,Ns); valsRE = zeros(K,Ns);

selected_min_electrode.LE = zeros(K,Ns); 
selected_min_electrode.RE = zeros(K,Ns); 

for subj = 1:Ns   
    [valsLE(:,subj),selected_max_electrode.LE(:,subj)]=maxk(mean(spatiotemp_info.LE{subj}(:,time_range),2),K);
    [~,selected_min_electrode.LE(:,subj)]=mink(mean(spatiotemp_info.LE{subj}(:,time_range),2),K);

    [valsRE(:,subj),selected_max_electrode.RE(:,subj)]=maxk(mean(spatiotemp_info.RE{subj}(:,time_range),2),K);
    [~, selected_min_electrode.RE(:,subj)]=mink(mean(spatiotemp_info.RE{subj}(:,time_range),2),K);
end

%% Add original electrode to the LOT and ROTs
for subj = 1:Ns
    ROTemitterNeighbours{subj} = [ROTemitterNeighbours{subj} originalMaxROT(subj)];
    LOTemitterNeighbours{subj} = [LOTemitterNeighbours{subj} originalMaxLOT(subj)];
end
%%
for subj = 1:Ns
    for k = 1:K
        findIdxLE = find(ROTemitterNeighbours{subj}==selected_max_electrode.LE(k,subj));
        findIdxRE = find(LOTemitterNeighbours{subj}==selected_max_electrode.RE(k,subj));
        findIdxFrontLE = find(FrontElectrodes{subj}==selected_min_electrode.LE(k,subj));
        findIdxFrontRE = find(FrontElectrodes{subj}==selected_min_electrode.RE(k,subj));
    
        if isempty(findIdxLE)
            selected_max_electrode.LE(k,subj) = nan;
        end
        if isempty(findIdxRE)
            selected_max_electrode.RE(k,subj) = nan;
        end
        if isempty(findIdxFrontLE)
            selected_min_electrode.LE(k,subj) = nan;
        end
        if isempty(findIdxFrontRE)
            selected_min_electrode.RE(k,subj) = nan;
        end
    end
end
%% Find second max and minFrontal electrodes
J = 1; % Rank second, third,... J+1th maximum in one hemisphere
algorithmicSecondMax.LE = zeros(J,Ns);
algorithmicSecondMax.RE = zeros(J,Ns);
algorithmicFrontMin.LE = zeros(J,Ns);
algorithmicFrontMin.RE = zeros(J,Ns);


for subj = 1:Ns
    findIdxLE = find(selected_max_electrode.LE(:,subj)==originalMaxROT(subj));
    findIdxRE = find(selected_max_electrode.RE(:,subj)==originalMaxLOT(subj));
    
    selected_max_electrode.LE(findIdxLE,subj) = nan;
    selected_max_electrode.RE(findIdxRE,subj) = nan;
    
    % SecondMax LE
    validIndices = ~isnan(selected_max_electrode.LE(:,subj));
    validElectrodes = selected_max_electrode.LE(validIndices,subj);
    algorithmicSecondMax.LE(:,subj) = validElectrodes(1:J);
    
    % SecondMax RE
    validIndices = ~isnan(selected_max_electrode.RE(:,subj));
    validElectrodes = selected_max_electrode.RE(validIndices,subj);
    algorithmicSecondMax.RE(:,subj) = validElectrodes(1:J);
    
    % FrontMin LE
    validIndices = ~isnan(selected_min_electrode.LE(:,subj));
    validElectrodes = selected_min_electrode.LE(validIndices,subj);
    algorithmicFrontMin.LE(:,subj) = validElectrodes(1:J);
    
    % FrontMin RE
    validIndices = ~isnan(selected_min_electrode.RE(:,subj));
    validElectrodes = selected_min_electrode.RE(validIndices,subj);
    algorithmicFrontMin.RE(:,subj) = validElectrodes(1:J);
end

