%% Compute stimulus information profiles shown in Figure 4B,C

%% Setup parameters

numWorkers = 30; % number of threads thar are going to be used for computation

bins = 2; % number of bins for binning

%% Load and prepare data

load metadata

timeSteps = length(time);

leftEyeLeftSensorMI = zeros(Ns, timeSteps);
rightEyeLeftSensorMI = zeros(Ns, timeSteps);
leftEyeRightSensorMI = zeros(Ns, timeSteps);
rightEyeRightSensorMI = zeros(Ns, timeSteps);

for subject = 1:Ns
    
    fname = sprintf('data/data_%s.mat',subjects{subject});
    dat = matfile(fname);

    % Select data just for the electrode with the highest MI
    LEfeeg = dat.feeg(:, :, dat.LEmaxMI);
    REfeeg = dat.feeg(:, :, dat.REmaxMI);

    % Compute temporal derivatives of eeg
    LEdfeeg = zeros(size(LEfeeg));
    LEdfeeg(:, 2:timeSteps) = LEfeeg(:, 2:timeSteps) - LEfeeg(:, 1:(timeSteps - 1));
    
    REdfeeg = zeros(size(REfeeg));
    REdfeeg(:, 2:timeSteps) = REfeeg(:, 2:timeSteps) - REfeeg(:, 1:(timeSteps - 1));

    % Load eye visibility
    eyesVisibility = dat.eyebubs;
    
    % Bin eye visibility
    [bLeftEyeVisbility, ~] = equi_binning2(eyesVisibility(:,1), bins, 0);
    [bRightEyeVisbility, ~] = equi_binning2(eyesVisibility(:,2), bins, 0);
    % Compute prob. dist. of the visibility
    
    numberOfTrials = size(eyesVisibility, 1);
    idx = 1:numberOfTrials;
    mlev = sparse(idx, bLeftEyeVisbility, 1, numberOfTrials, bins * bins, numberOfTrials);
    mrev = sparse(idx, bRightEyeVisbility, 1, numberOfTrials, bins * bins, numberOfTrials);
    
    plev = sum(mlev) / numberOfTrials;
    prev = sum(mrev) / numberOfTrials;
    
    hlev = -dot(plev, log2(plev + eps));
    hrev = -dot(prev, log2(prev + eps));

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

    %% Computation of the information metrics

    tic

    parfor time = 1:timeSteps
        
        sprintf('Time %d\n', time);
        
         
            % binning and construction of new variables from the
            % combination of absolute value and the derivative of eeg

            % binning to values = 1 .. bins

            [bLEfeeg, ~] = equi_binning2(LEfeeg(:, time), bins, 0);
            [bLEfdeeg, ~] = equi_binning2(LEdfeeg(:, time), bins, 0);
            
            [bREfeeg, ~] = equi_binning2(REfeeg(:, time), bins, 0);
            [bREfdeeg, ~] = equi_binning2(REdfeeg(:, time), bins, 0);

            % Compute prob. dist. and MI
            
            bLEf2Deeg = (bLEfeeg - 1) .* bins + bLEfdeeg;
            bREf2Deeg = (bREfeeg - 1) .* bins + bREfdeeg;
            
            mleeg = sparse(idx, bLEf2Deeg, 1, numberOfTrials, bins * bins, numberOfTrials);
            pleeg = sum(mleeg) / numberOfTrials;
            mreeg = sparse(idx, bREf2Deeg, 1, numberOfTrials, bins * bins, numberOfTrials);
            preeg = sum(mreeg) / numberOfTrials;
            
                                 
            plevLEeg = mlev' * mleeg / numberOfTrials;
            prevLEeg = mrev' * mleeg / numberOfTrials;
            plevREeg = mlev' * mreeg / numberOfTrials;
            prevREeg = mrev' * mreeg / numberOfTrials;
            
            hleeg = -dot(pleeg, log2(pleeg + eps));
            hreeg = -dot(preeg, log2(preeg + eps));
            
            hlevLEeg = -dot(nonzeros(plevLEeg), log2(nonzeros(plevLEeg) + eps));
            hrevLEeg = -dot(nonzeros(prevLEeg), log2(nonzeros(prevLEeg) + eps));
            hlevREeg = -dot(nonzeros(plevREeg), log2(nonzeros(plevREeg) + eps));
            hrevREeg = -dot(nonzeros(prevREeg), log2(nonzeros(prevREeg) + eps));
            
            leftEyeLeftSensorMI(subject, time) = hleeg + hlev - hlevLEeg;
            rightEyeLeftSensorMI(subject, time) = hleeg + hrev - hrevLEeg;
            leftEyeRightSensorMI(subject, time) = hreeg + hlev - hlevREeg;
            rightEyeRightSensorMI(subject, time) = hreeg + hrev - hrevREeg;
            
    end

    toc
    
end

%% Save the results

save(['./output/miEyeVisXEEG' num2str(bins) '.mat'], 'leftEyeLeftSensorMI', 'rightEyeLeftSensorMI', 'leftEyeRightSensorMI', 'rightEyeRightSensorMI');



