% Analyze the information transfer in the eeg data discounting alternative 
% sender effect using cFIT

%% Setup parameters

numWorkers = 12; % number of threads thar are going to be used for computation
maxDelay = 60; % maximal delay of past values against current values of random variables
bins = 2; % number of bins for binning

portionsOfDataset = 1; % sets how many times should be the dataset cut by 2 in order to compute results

doShuffle = 1;
doBootstrap = 1;
doQE = 0;
%numOfShuffling = 100;
numOfShuffling = 10;
numOfBootstrap = 10;

subject = 1;
load algorithmicSecondMax_electrodes
load algorithmicFrontMin_electrodes

conditioningArea = 'emitter'; % can be 'emitter' or 'receiver'
%% Load and prepare data

load metadata

timeSteps = length(time);

leftFIT = zeros(portionsOfDataset, Ns, timeSteps, maxDelay);
leftdfi = zeros(portionsOfDataset, Ns, timeSteps, maxDelay);
leftdi = zeros(portionsOfDataset, Ns, timeSteps, maxDelay);

rightFIT = zeros(portionsOfDataset, Ns, timeSteps, maxDelay);
rightdfi = zeros(portionsOfDataset, Ns, timeSteps, maxDelay);
rightdi = zeros(portionsOfDataset, Ns, timeSteps, maxDelay);

FITRE_complete = rightdfi; cFITRE = rightdfi;
FITRESh =  zeros(portionsOfDataset, Ns, timeSteps, maxDelay, numOfShuffling);
cFITRESh = FITRESh;  
FITREBoot = zeros(portionsOfDataset, Ns, timeSteps, maxDelay, numOfBootstrap);
cFITREBoot = FITREBoot;
%% Start pool with correct number of workers 

if (numWorkers > 1)
    if (isempty(gcp('nocreate')))
        parpool(numWorkers);
    else
        poolTemp = gcp;
        if(poolTemp.NumWorkers < numWorkers)
            delete(gcp)
            parpool(numWorkers);
        end
    end
end

for fraction = 1:portionsOfDataset


%% The cycle over subjects with computations themselves

for subject = 1:Ns
    
    fprintf('Computing subject %d...\n', subject)
    
    fname = sprintf('data/data_%s.mat',subjects{subject});
    dat = matfile(fname);
    % Select data just for the electrode with the highest MI
    LEeeg = dat.feeg(1:(2^(fraction - 1)):size(dat.feeg, 1), :, dat.LEmaxMI);
    REeeg = dat.feeg(1:(2^(fraction - 1)):size(dat.feeg, 1), :, dat.REmaxMI);
    
    % Select RE confounder
    if strcmp(alternative_sender,'FontMin')
        REconfeeg = dat.feeg(1:(2^(fraction - 1)):size(dat.feeg, 1), :, algorithmicFrontMin.RE(subject));
    elseif strcmp(alternative_sender,'SecondMax')
        REconfeeg = dat.feeg(1:(2^(fraction - 1)):size(dat.feeg, 1), :, algorithmicSecondMax.RE(subject));
    end
    % Load eye visibility
    eyesVisibility = dat.eyebubs;
    
    % Bin eye visibility
    [bLeftEyeVisbility, ~] = equi_binning2(eyesVisibility(:,1), bins, 0);
    [bRightEyeVisbility, ~] = equi_binning2(eyesVisibility(:,2), bins, 0);
    bLeftEyeVisbility = bLeftEyeVisbility(1:(2^(fraction - 1)):size(dat.feeg, 1));
    bRightEyeVisbility = bRightEyeVisbility(1:(2^(fraction - 1)):size(dat.feeg, 1));

    % Compute temporal derivatives of eeg
    LEdeeg = zeros(size(LEeeg));
    LEdeeg(:, 2:timeSteps) = LEeeg(:, 2:timeSteps) - LEeeg(:, 1:(timeSteps - 1));

    REdeeg = zeros(size(REeeg));
    REdeeg(:, 2:timeSteps) = REeeg(:, 2:timeSteps) - REeeg(:, 1:(timeSteps - 1));
    
    REconfdeeg = zeros(size(REconfeeg));
    REconfdeeg(:, 2:timeSteps) = REconfdeeg(:, 2:timeSteps) - REconfdeeg(:, 1:(timeSteps - 1));
    
    stimulusLeft = bLeftEyeVisbility;
    stimulusRight = bRightEyeVisbility;
    
    %%% Shuffle the stimulus for establishing p value
     
    if (size(stimulusRight, 2) < size(stimulusRight,1))
        stimulusRight = stimulusRight';
    end
    
    if (size(stimulusLeft, 2) < size(stimulusLeft,1))
        stimulusLeft = stimulusLeft';
    end
    
    
    %% Computation of the information metrics

    parfor time = 151:351
        disp(['Time ', num2str(time), ' subj ', num2str(subject)])
%         sprintf('Time %d...\n', time)
        for delay = 1:maxDelay
            if ((time - delay) > 0) 
                % binning and construction of new variables from the
                % combination of absolute value and the derivative of eeg

                % binning to values = 1 .. bins
                
                [bLEeegCurrent, ~] = equi_binning2(LEeeg(:, time), bins, 0);
                [bLEdeegCurrent, ~] = equi_binning2(LEdeeg(:, time), bins, 0);
                [bLEeegPast, ~] = equi_binning2(LEeeg(:, time - delay), bins, 0);
                [bLEdeegPast, ~] = equi_binning2(LEdeeg(:, time - delay), bins, 0);
                
                [bREeegCurrent, ~] = equi_binning2(REeeg(:, time), bins, 0);
                [bREdeegCurrent, ~] = equi_binning2(REdeeg(:, time), bins, 0);
                [bREeegPast, ~] = equi_binning2(REeeg(:, time - delay), bins, 0);
                [bREdeegPast, ~] = equi_binning2(REdeeg(:, time - delay), bins, 0);                
                
                
                if strcmp(conditioningArea, 'emitter')
                    addDelay = 0;
                    [bREconfeegPast, ~] = equi_binning2(REconfeeg(:, time - delay- addDelay), bins, 0);
                    [bREconfdeegPast, ~] = equi_binning2(REconfdeeg(:, time - delay - addDelay), bins, 0);
                elseif strcmp(conditioningArea, 'receiver')
                    addDelay = maxDelayROT(subject);
                    [bREconfeegPast, ~] = equi_binning2(REconfeeg(:, time - addDelay), bins, 0);
                    [bREconfdeegPast, ~] = equi_binning2(REconfdeeg(:, time - addDelay), bins, 0);
                end
                %%% Left Eye visibility, Right to Left transfer
                pastX = (bREeegPast - 1) .* bins + bREdeegPast; 
                Y = (bLEeegCurrent - 1) .* bins + bLEdeegCurrent;
                pastY = (bLEeegPast - 1) .* bins + bLEdeegPast;

                % compute probability distribution
                
                if (size(pastX, 2) < size(pastX,1))
                    pastX = pastX';
                end
                
                if (size(Y, 2) < size(Y,1))
                    Y = Y';
                end
                
                if (size(pastY, 2) < size(pastY,1))
                    pastY = pastY';
                end
                
                %%% Right Eye visibility, Left to Right transfer
                pastXR = (bLEeegPast - 1) .* bins + bLEdeegPast; 
                YR = (bREeegCurrent - 1) .* bins + bREdeegCurrent;
                pastYR = (bREeegPast - 1) .* bins + bREdeegPast;
                pastZR = (bREconfeegPast - 1) .* bins + bREconfdeegPast;
                
                % compute probability distribution
                
                if (size(pastXR, 2) < size(pastXR,1))
                    pastXR = pastXR';
                end
                
                if (size(YR, 2) < size(YR,1))
                    YR = YR';
                end
                
                if (size(pastYR, 2) < size(pastYR,1))
                    pastYR = pastYR';
                    pastZR = pastZR';
                end
                
                [rightdi(fraction, subject, time, delay), rightdfi(fraction, subject, time, delay), rightFIT(fraction, subject, time, delay), cFITRE(fraction, subject, time, delay)] = ...
                    compute_FIT_TE_cFIT(stimulusRight, pastXR, YR, pastYR);
                
                [leftdi(fraction, subject, time, delay), leftdfi(fraction, subject, time, delay), leftFIT(fraction, subject, time, delay), cFITLE(fraction, subject, time, delay)] = ...
                    compute_FIT_TE_cFIT(stimulusLeft, pastX, Y, pastY);
                
                if (doShuffle)
                    XLSh = zeros(numOfShuffling, numel(pastX));
                    XRSh = zeros(numOfShuffling, numel(pastXR));
                    SvalL = unique(stimulusLeft);
                    SvalR = unique(stimulusRight);
                    
                    for s = 1:numOfShuffling
                        for sS = 1:numel(SvalL)
                            idxL = (stimulusLeft == SvalL(sS));
                            idxR = (stimulusRight == SvalR(sS));
                            hXL = pastX(idxL);
                            hXR = pastXR(idxR);
                            ridxL = randperm(sum(idxL));
                            ridxR = randperm(sum(idxR));
                            XLSh(s,idxL) = hXL(ridxL);
                            XRSh(s,idxR) = hXR(ridxR);
                        end
                    end
                            
                    for s = 1:numOfShuffling

                        [ ~, ~, FITRESh(fraction, subject, time, delay,s), cFITRESh(fraction, subject, time, delay,s)] = compute_FIT_TE_cFIT(stimulusRight, XRSh(s,:), YR, pastYR, pastZR)

                    end
                end  
               
            end
        end
    end
end

if strcmp(alternative_sender,'FontMin')
    save allSubjects_EEG_FrontMinConfounders_onRE
elseif strcmp(alternative_sender,'SecondMax')
    save allSubjects_EEG_SecondMaxConfounders_onRE
end