%% This functions has the mission to extract the eye tracker data of each participant and save
% them in txt files for posterior analyses in R.
% Alexis Perez Bellido october-2018


clear
%
overwrite = 0;

%% Select your parameters to select the data where you want to perform the analysis

format.baseline = [-0.2 0];
format.lpfreq = 30;
%format.hpfreq = 2;
format.resamplefs = 200;
format.detrend = 'no';
format.trigger_ref = 'target'; % what is the trigger used as a reference


%% Select your parameters to identify the analysis that you want to perform. Change everytime accordingly


param.trials_sel = 'all' %? correct/all % {'MEG'} what trials go in the analysis?
param.channels_sel = {'allchan'}; %{'MEG', 'UADC007'} ; % all, {'MLO*','MRO*','MZO*'}  % {'MEG', 'UADC007'} what channels go in the analyses ,'UADC007'
param.design = '2mod2vRAW'; %''2block2mod2vRAW','2mod2vRAW' %; '2blocksSDT';  % '2mod2v '2blocks2mod2v'; '2blocksSDT'  '2mod2v'% how are u organizing the data?

% for subject 10 there might some problems (onset is weird). Also for 15.
% some subjects might not have the eye data recorded here.
%

for isubj = [39 40] % 10 11 12 13 14 15 
    
    % isubj = 39%  [12 13 14]
    % Define globals
    global G;
    G = [];
    
     % First load the data file
    G.data_file_name = ['dformat_s',num2str(isubj)];
    G.anal_file_name = ['eye_',param.design,'_',param.trials_sel,'_s',num2str(isubj)];
    
    %% General paths
    exp_main_path = '/project/3018029.06';
    
    SUBJS = {['0',num2str(isubj)]}
    curSubj = SUBJS{1}
    
    addpath(genpath( [exp_main_path filesep 'Matlab' filesep 'MEG2_scripts'])); % general processing scripts in this exp.

    G.subjPath = [exp_main_path filesep 'MEG' filesep 'AVE_MEG2' filesep curSubj filesep]; % main subject folder
    G.dataPath = [exp_main_path filesep 'raw', filesep 'sub-' curSubj  filesep]; % subjecf MEG data
    G.outputPath = [G.subjPath ,'func', filesep]; % subject saving artifact and preprocessing matrixes
    G.cfgbehavPath = [exp_main_path filesep 'Behavioral' filesep 'AVE_MEG2' filesep 'EYE_MOVEMENTS_DATA' filesep] ; % behavioral data of participant
    
    subjscriptsPath = [G.subjPath ,'analysis', filesep]; %  participant specific scripts
    data_out =  G.anal_file_name;
    
%    addpath(genpath('/project/3018029.06/Matlab/MEG_scripts/'));
    
    G.inputPath = [G.dataPath  'ses-meg01' filesep ];
    
    

    if exist([G.outputPath G.data_file_name,'.mat'],'file')
        load([G.outputPath G.data_file_name,'.mat']);
        % file exist; checking whether the analyses is new...
        
        number_of_datasets = size(dformat,1);
        disp(['there are ' num2str(number_of_datasets) ' datasets in this participant file'])
        
        for i = 1 : number_of_datasets % lets find where are the data that we want to perform the analysis
            data_exist = 0;
            
            if strcmp(num2str(dformat(i).format.baseline), num2str(format.baseline)) & ...
                    strcmp(num2str(dformat(i).format.lpfreq) ,num2str(format.lpfreq)) & ...
                    strcmp(num2str(dformat(i).format.resamplefs) ,num2str(format.resamplefs)) & ...
                    strcmp(num2str(dformat(i).format.detrend) ,num2str(format.detrend)) & ...
                    strcmp(num2str(dformat(i).format.trigger_ref) ,num2str(format.trigger_ref));  data_exist = 1; end;
            
            if (data_exist > 0); n = i;
                disp(['the dataset is in field number ' num2str(n),' now preparing to validate the analysis type  ' ]); end; % overwrite in case that this is the dataset repeated
            if (i == number_of_datasets & data_exist == 0); error('This particular data preproc doesnt exist for this participant...'); end; % if we did not find this analyses already performed, save data in the next folder.
        end
        
    else
        error(['No datasets for this participant. You need to create the datase using prep_data_format function'])
        
    end
    
    eye_link_channels =  {'UADC007' 'UADC005' 'UADC006' };
    channels_names = {'pupil', 'x','y'}
    
    cfg = [];
    cfg.channel =  eye_link_channels;
    
    data = ft_selectdata(cfg, dformat(n).data);
    
   % cfg = [];
   % cfg.resamplefs  = 100;
   % [data] =  ft_resampledata(cfg, data)
    
    %% To reduce the size of the files cutting time and downsampling at
    % 100hz
    % cutting some time before the stimulus is present
    cfg = [];
    cfg.toilim  = [-0.5 1.25]
    data = ft_redefinetrial(cfg, data);

    
    % adding a new column to datainfo
    data.trialinfo =  [data.trialinfo zeros(size(data.trialinfo,1),1)];
    
    % re-coding yes / no responses based on signal noise and correct
    % responses. If is correct and there was signal is because the
    % participant pressed YES
    data.trialinfo(find(data.trialinfo(:,4) == 1 & data.trialinfo(:,10) == 1),12) = 1;


    %%  Conditions labels
    condition = {'noise','signal'};
    modalities = {'visual','audiovisual'};
    visualfields = {'center','periphery'};
    response = {'no','yes'};
    session = [1:6]
    
    conds = unique(data.trialinfo(:,4)); % signal noise
    modals = unique(data.trialinfo(:,5)); % visual audiovisual
    vfields = unique(data.trialinfo(:,6)); % center periphery
    resps = unique(data.trialinfo(:,12)); % yes or no
    sessions = unique(data.trialinfo(:,2)); % session?
    corrects = unique(data.trialinfo(:,10)); % yes or no
    
    time = data.time;
    numT = length(data.time);
    numN = length(data.trial);
    
    
    % instantiating a matrix to create data
    eye_data = []
 
    for icond = 1 : length(condition)
        for imod = 1 : length(modals)
            for ivfield = 1 : length(vfields)
                for icorr = 1 :  length(corrects)
                    % for isession = 1 :  length(sessions)
                  
                    eyes_vectors = [];
                    
                    trials = find( ... % data.trialinfo(:, 2) == session(isession)
                        data.trialinfo(:, 4) == conds(icond) ...
                        & data.trialinfo(:, 5) == modals(imod) ...
                        & data.trialinfo(:, 6) == vfields(ivfield) ...
                        & data.trialinfo(:, 10) == corrects(icorr) );
            
                    if ~isempty(trials) % if there are no trials for this, just create a fake matrix to record data. then erase the matrix in the analyses                   
                        for i_channel =  1 : length(eye_link_channels)
                            
                              i_channel = 1                          
                            %%  Blinks during stimulus
                            % I only read from the eyetracker (no EOG data)
                            % stimOn = [-0.1, 0.1];
                     %      data.trialinfo(trials, :)   
                      %     data_preproc_clean.trialinfo(trials, :)    
                             cfg = [];
                             cfg.channel = eye_link_channels{i_channel};
                             tmp_data = ft_selectdata(cfg, data);
                             X_eye =   tmp_data.trial; 
                          
                             %zscore all the trials values for this channel
                             %and attach to dataset
                             tmp_data.trial = (X_eye - nanmean(X_eye(:))) ./ nanstd(X_eye(:));
                             figure;  imagesc( mean(squeeze( tmp_data.trial),1))
                             % select the zscored trials for this condition
                             % combination
                             cfg = [];
                             cfg.trials  = trials;
                             tmp_data = ft_selectdata(cfg, tmp_data);
                            
                            RT_avg = mean(tmp_data.trialinfo(:,11)) ; % calculate the average RT
                            nobs = size(tmp_data.trialinfo,1); % number of observed trials in this manifold
                            
                            disp('done')
                            % not sure what this reshape is doing
                            % X = reshape((tmp_data.trial), [1, length(tmp_data.time), length(tmp_data.trial)]);
                            X_eye = squeeze(tmp_data.trial);
                            % mask data in all trials below this value.
                         
                            %imagesc(time,1 , X_eye)
                
                            %% generate mask with eye blinks.
                            if i_channel == 1
                                cutoff_z_pupDil = 5;
                                blinkMask =  (abs(X_eye) > cutoff_z_pupDil);
                                
                              
                                if nobs > 1
                                blinks_time = mean(blinkMask,1);
                                else
                                blinks_time =  blinkMask'; 
                                end
                            end
                            % replace eye blinks with NaN
                            X_eye(blinkMask) = NaN;
                            
                            % imagesc( time, 1,  X_eye)
                            
                            
                            
                            % stimOnMask = repmat(inrange(tmp_data.time{1}, stimOn + [-1 1]*stimOnBuffer), [length(tmp_data.trial), 1]);
                            % {Visualize blinks
                            % figure;
                            % subplot(4, 1, 1); imagesc(1:numN, time, X_EOGv'); colorbar;
                            % subplot(3, 1, 1); imagesc(1:numN, time, X_eye'); colorbar;
                            % subplot(3, 1, 2); imagesc(1:numN, time, blinkMask'); colorbar;
                            % subplot(3, 1, 3); imagesc(1:numN, time, stimO
                            
                            % ntrials = size( X_eye, 1);
                            % ntimepoints = size( X_eye, 2);
                          
                            if nobs > 1
                                aux_xyp = [nanmedian(X_eye,1)]; % saving median vectors of x, y and pupil vectors (i use the median because it should be less affected by extreme values)
                                eyes_vectors = [eyes_vectors,aux_xyp'];
                                % imagesc( time, 1,  aux_xyp)
                            else
                                aux_xyp = X_eye; % if there is only one trial it doesnt make sense to average.
                                eyes_vectors = [eyes_vectors,aux_xyp];
                                % imagesc(time, 1, eyes_vectors')
                            end
                        end
                        
                        % attaching eyeblinks per unit of time
                        eyes_vectors = [eyes_vectors, blinks_time'];
                    end
                    
                    
                        if isempty(trials) % if there are no trials for this, just create a fake matrix to record data. then erase the matrix in the analyses
                            aux = num2cell(zeros( length(time),20));
                        else 
                            
                            aux = [num2cell(repmat(isubj,length(time),1)), num2cell(repmat([ tmp_data.trialinfo(1, :),nobs, RT_avg], length(time),1) ), num2cell(eyes_vectors),  num2cell(time')];
                        end
                        
                           % aux = [num2cell(repmat(tmp_data.trialinfo(3, 3 : 12), ntimepoints,1)), num2cell( mean(X_eye))', repmat(channels_names( i_channel),ntimepoints,1)];
                           eye_data = [eye_data ; aux]; 
                  
                end             
            end
        end
    end
    
    % save here
    
    save(['/project/3018029.06/Behavioral/EYE_MOVEMENTS_DATA/' data_out ] ,'eye_data','-v7.3')
       
    fileID = fopen(['/project/3018029.06/Behavioral/EYE_MOVEMENTS_DATA/' data_out '.txt'],'wt'); % try to write on top
    formatSpec = '%i %i %i %i %i %i %i %f %i %i %i %f %i %i %f %f %f %f %f %f\r\n';
   
    for row = 1 : length(eye_data)
        fprintf(fileID, formatSpec, eye_data{row,:}) ;
    end
    fclose(fileID);


end






figure;
plot(time,mean(X_eye))