% Analyze the information transfer in the eeg data

% This is mcOriginakl + perm new null

%% 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;
shuffType = 'simple'; % either simple or conditioned (run for both, the maximum is compute in the plot script)
doQE = 0;
numOfShuffling = 10;

subject = 1;

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

    if (doShuffle)
        for s = 1:numOfShuffling
            index = randperm(numel(stimulusLeft));
            stimShLeft(s, :) = stimulusLeft(index);
            index = randperm(numel(stimulusRight));
            stimShRight(s, :) = stimulusRight(index);
        end
    end

    %% Computation of the information metrics

    tic

    for time = 151:351
        if mod(time-1,10) == 0
            disp(['Time ', num2str(time)])
        end
%         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);                
                
                %%% 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;
                
                
                % 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);
                
                if (doShuffle)
                    XLSh = zeros(numOfShuffling, numel(pastX));
                    XRSh = zeros(numOfShuffling, numel(pastXR));
                    SvalL = unique(stimulusLeft);
                    SvalR = unique(stimulusRight);
                    
                    for s = 1:numOfShuffling
                        stimLsh = stimShLeft(s,:);
                        stimRsh = stimShRight(s,:);
                        
                        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

                    if strcmp(shuffType,'conditioned')
                        parfor s = 1:numOfShuffling
                            [~, leftdfiSh(fraction, subject, time, delay, s), leftFITSh(fraction, subject, time, delay, s)] = compute_FIT_TE(stimulusLeft, XLSh(s,:), Y, pastY);
                            [~, rightdfiSh(fraction, subject, time, delay, s), rightFITSh(fraction, subject, time, delay, s)] = compute_FIT_TE(stimulusRight, XRSh(s,:), YR, pastYR);
                        end
                        
                    elseif strcmp(shuffType,'simple')
                        parfor s = 1:numOfShuffling
                            [~, leftdfiSh(fraction, subject, time, delay, s), leftFITSh(fraction, subject, time, delay, s)] = compute_FIT_TE(stimLsh, pastX, Y, pastY);
                            [~, rightdfiSh(fraction, subject, time, delay, s), rightFITSh(fraction, subject, time, delay, s)] = compute_FIT_TE(stimRsh, pastXR, YR, pastYR);
                        end
                    end
                        
                    end
                end                 
            end
        end
    end
    
    toc
end

%% Save the results

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

leftEyeInf.FITSh = leftFITSh;
leftEyeInf.dfiSh = leftdfiSh;

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

rightEyeInf.FITSh = rightFITSh;
rightEyeInf.dfiSh = rightdfiSh;

fname = ['FIT_contralateral_',shuffType];
save(fname,'leftEyeInf','rightEyeInf','bins','-v7.3')