% Analyze the information transfer in the eeg data

%% 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;
numOfShuffling = 1; % Shuff ipsilateral not used 

%% 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);

leftFITSh = zeros(portionsOfDataset, Ns, timeSteps, maxDelay, numOfShuffling);
leftdfiSh = zeros(portionsOfDataset, Ns, timeSteps, maxDelay, numOfShuffling);

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

rightFITSh = zeros(portionsOfDataset, Ns, timeSteps, maxDelay, numOfShuffling);
rightdfiSh = zeros(portionsOfDataset, Ns, timeSteps, maxDelay, numOfShuffling);

%% 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);

    % 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));
    
    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
    
    stimShLeft = zeros(numOfShuffling, length(stimulusLeft));
    stimShRight = zeros(numOfShuffling, length(stimulusRight));

    %% Computation of the information metrics

    parfor time = 151:351
        disp(['Time ', num2str(time), ', subj ', num2str(subject)])
        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);                
                
                %%% Left Eye visibility, Left to Right transfer
                pastX = (bLEeegPast - 1) .* bins + bLEdeegPast; 
                Y = (bREeegCurrent - 1) .* bins + bREdeegCurrent;
                pastY = (bREeegPast - 1) .* bins + bREdeegPast;
                

                % 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, Right to Left transfer
                pastXR = (bREeegPast - 1) .* bins + bREdeegPast; 
                YR = (bLEeegCurrent - 1) .* bins + bLEdeegCurrent;
                pastYR = (bLEeegPast - 1) .* bins + bLEdeegPast;
                
                
                % 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';
                end
                
                [rightdi(fraction, subject, time, delay), rightdfi(fraction, subject, time, delay), rightFIT(fraction, subject, time, delay)] = compute_FIT_TE(stimulusRight, pastXR, YR, pastYR);
                
                [leftdi(fraction, subject, time, delay), leftdfi(fraction, subject, time, delay), leftFIT(fraction, subject, time, delay)] = compute_FIT_TE(stimulusLeft, pastX, Y, pastY);
                
            end
        end
    end
    end
end

%% Save the results

leftEyeInf.rDFI = leftFIT;
leftEyeInf.dfi = leftdfi;
leftEyeInf.di = leftdi;

leftEyeInf.rDFISh = leftFITSh;
leftEyeInf.dfiSh = leftdfiSh;

rightEyeInf.rDFI = rightFIT;
rightEyeInf.dfi = rightdfi;
rightEyeInf.di = rightdi;

rightEyeInf.rDFISh = rightFITSh;
rightEyeInf.dfiSh = rightdfiSh;
%%
fname = ['FIT_ipsilateral'];
save(fname,'leftEyeInf','rightEyeInf','bins','-v7.3')