clear all;

RandStream.setGlobalStream(RandStream('mt19937ar','seed',sum(100*clock)))

mainfolder=''

addpath(genpath(['subfunctions/']))

%loading list of participants (good perfomance)

load('good_partis.mat')

partis=good_partis;
partis_sub=good_partis;

%% Parameters to be modify per analysis

%% object of interest conf

delay=1; %(1 or 2). This variable set everything for analysing different epochs (relative to delay 1 or 2)

removingalarms=0; % Removing epochs that includes or were preceded by an alarm message (1=yes, 0=no)

% selecting stimulus of interest for each iteration (1 = stim presented 1st. 2 stim presented 2nd )
conf.ooi{1}=1;
conf.ooi{2}=1;
conf.ooi{3}=2;
conf.ooi{4}=2;

conf.ooi{5}=1;
conf.ooi{6}=1;
conf.ooi{7}=2;
conf.ooi{8}=2;

conf.ooi{9}=1;
conf.ooi{10}=1;
conf.ooi{11}=2;
conf.ooi{12}=2;

if delay== 1
    if removingalarms
        foldpath='results/noalarm/mean_recentered_th100_messagesOFF/euclideanmodels';
    else
        foldpath='results/mean_recentered_th100_messagesOFF/euclideanmodels';
    end
elseif delay==2
    if removingalarms
        foldpath='results/noalarm/mean_recentered_th100_messagesOFF/euclideanmodels/2nd_cue';
    else
        foldpath='results/mean_recentered_th100_messagesOFF/euclideanmodels/2nd_cue';
    end
end

%selecting the same of the folder to save the output in each iteration.
conf.folder{1} = [foldpath '/single_trial_with_between_RSA_object_1_cued_cont_Pearson/'];
conf.folder{2} = [foldpath '/single_trial_with_between_RSA_object_1_uncued_cont_Pearson/'];
conf.folder{3} = [foldpath '/single_trial_with_between_RSA_object_2_cued_cont_Pearson/'];
conf.folder{4} = [foldpath '/single_trial_with_between_RSA_object_2_uncued_cont_Pearson/'];

conf.folder{5} = [foldpath '/single_trial_with_between_RSA_object_1_cued_cat_udlr_Pearson/'];
conf.folder{6} = [foldpath '/single_trial_with_between_RSA_object_1_uncued_cat_udlr_Pearson/'];
conf.folder{7} = [foldpath '/single_trial_with_between_RSA_object_2_cued_cat_udlr_Pearson/'];
conf.folder{8} = [foldpath '/single_trial_with_between_RSA_object_2_uncued_cat_udlr_Pearson/'];

conf.folder{9} = [foldpath '/single_trial_with_between_RSA_object_1_cued_cat_udlr45_Pearson/'];
conf.folder{10} = [foldpath '/single_trial_with_between_RSA_object_1_uncued_cat_udlr45_Pearson/'];
conf.folder{11} = [foldpath '/single_trial_with_between_RSA_object_2_cued_cat_udlr45_Pearson/'];
conf.folder{12} = [foldpath '/single_trial_with_between_RSA_object_2_uncued_cat_udlr45_Pearson/'];

%do we need to create this folders? uncomment this:

% for fo = 1:length(conf.folder)
%
%     papa=foldpath%'/LindeDomingoJuan/Memoreye_DATA/RSA/eyedata/results/mean_recentered_th120/euclideanmodels'
%     foi=conf.folder{fo};
%     hijo=foi(length(papa)+1:end)
% mkdir(papa, hijo)
% end

%selecting the cue in each iteration
conf.retrooi{1}=1;
conf.retrooi{2}=2;
conf.retrooi{3}=2;
conf.retrooi{4}=1;

conf.retrooi{5}=1
conf.retrooi{6}=2
conf.retrooi{7}=2
conf.retrooi{8}=1

conf.retrooi{9}=1
conf.retrooi{10}=2
conf.retrooi{11}=2
conf.retrooi{12}=1

%selection the model for each iteration
conf.gmodelmat{1}=creatcontmodel(1);
conf.gmodelmat{2}=creatcontmodel(1);
conf.gmodelmat{3}=creatcontmodel(1);
conf.gmodelmat{4}=creatcontmodel(1);

conf.gmodelmat{5}=creatcatcardeucl;
conf.gmodelmat{6}=creatcatcardeucl;
conf.gmodelmat{7}=creatcatcardeucl;
conf.gmodelmat{8}=creatcatcardeucl;

conf.gmodelmat{9}=creatcatcardeucl45;
conf.gmodelmat{10}=creatcatcardeucl45;
conf.gmodelmat{11}=creatcatcardeucl45;
conf.gmodelmat{12}=creatcatcardeucl45;
%% Starting the loop for all analyses

for ii=1:12
    
    %folder to save results
    
    folder=conf.folder{ii}
    
    %object of interest (stim first of second)
    ooi=conf.ooi{ii};
    
    %taking only trials with this retrocue
    retrooi=conf.retrooi{ii};
    
    selectingretrocue=1;
    
    %info about the analysis
    analysis_info = ' '
    
    %selecting the model that we will use
    
    modelmat=conf.gmodelmat{ii}; %
    
    %%
    
    for ppp=1:length(partis) %loop of participants
        
        %% Loading dataset
        
        if delay== 1
            toload =[mainfolder 'eyelink_preprocessed/mean_recentered_th100_messagesOFF/object_1_onset_11_seconds/' partis{ppp,:} '.mat']; %for 1st delay
        elseif delay==2
            toload =[mainfolder 'eyelink_preprocessed/mean_recentered_th100_messagesOFF/2nd_cue_onset_17_5__seconds/' partis{ppp,:} '.mat']; %for 2nd delay
        end
        load(toload)
        
        %Loading logfiles
        
        subtable=[];
        
        if strcmp( partis_sub{ppp,:},'p11')
            load([mainfolder 'logfiles/resultfile_p11_table_repaired.mat'])
            subtable=sub
        else          
            sub = tdfread([mainfolder 'logfiles/resultfile_' partis_sub{ppp,:} '.txt'],'tab'); %Logfile for this participant
            subtable=struct2table(sub);
        end

        %% marking whether trials contains a FIXATE! message before the retention period of interest
        
        subtable.delay1_invalid(1:size(subtable,1))=0;
        subtable.delay2_invalid(1:size(subtable,1))=0;
        
        for tr= 1:size(subtable,1)
            if subtable.object_null_1(tr)==1|subtable.object_null_2(tr)==1|subtable.delay_1_null(tr)==1
                subtable.delay1_invalid(tr)=1;
            end
            if subtable.object_null_1(tr)==1|subtable.object_null_2(tr)==1|subtable.delay_1_null(tr)==1|subtable.delay_2_null(tr)==1
                subtable.delay2_invalid(tr)=1;
            end
            
            if subtable.testing_uncued(tr)==0
                subtable.delay2_invalid(tr)=NaN;
            end
        end
        
        %% Adding this information in the trialinfo
        % To make sure that repeated trials (if any) has a different trial number
        
        %% Logfile
        
        foundrep=[];
        for tr =2:size(subtable.trial,1)
            if abs(subtable.trial(tr)-subtable.trial(tr-1))>400
                foundrep=tr
            end
        end
        subtable.trial(foundrep:end)=subtable.trial(foundrep:end)+1000
        
        %% Eye data
        
        foundrep=[]
        
        for tr =2:size(data_epoch_thr.trialinfo.trial,1)
            if abs(data_epoch_thr.trialinfo.trial(tr)-data_epoch_thr.trialinfo.trial(tr-1))>400
                foundrep=tr
            end
        end
        data_epoch_thr.trialinfo.trial(foundrep:end)=data_epoch_thr.trialinfo.trial(foundrep:end)+1000
        
        %%
        
        for tr = 1:size(data_epoch_thr.trialinfo,1)
            
            toi=data_epoch_thr.trialinfo.trial(tr);
            found=find(subtable.trial==toi);
            
            data_epoch_thr.trialinfo.delay1_invalid(tr)=subtable.delay1_invalid(found);
            data_epoch_thr.trialinfo.delay2_invalid(tr)=subtable.delay2_invalid(found);
        end
        
        %% If the position is out of the threshold, it will be a NaN
        
        %     thr=100 %threshold in pixels
        %     thr=70 %threshold in pixels
        %
        %     data_art_reject.eyedat(abs(data_epoch_recenter.eyedat(:,:,:))>thr)=NaN;
        
        data_art_reject=data_epoch_thr; % this threshold step has been already implementes in the 2022 preprocessing pipeline
        
        %% Keeping trials based on retrocue and (optional) based on FIXATE! messages
        
        if selectingretrocue
            
            if removingalarms&delay==1
                toi =data_art_reject.trialinfo.retrocue==retrooi&data_epoch_thr.trialinfo.delay1_invalid~=1;
            elseif removingalarms&delay==2
                toi =data_art_reject.trialinfo.retrocue==retrooi&data_epoch_thr.trialinfo.delay2_invalid~=1;
            elseif removingalarms==0
                toi =data_art_reject.trialinfo.retrocue==retrooi;
            end
            data_art_reject.eyedat=data_art_reject.eyedat(toi,:,:);
            data_art_reject.trialinfo=data_art_reject.trialinfo(toi,:);
            
        else
            
            if removingalarms&delay==1
                toi =data_epoch_thr.trialinfo.delay1_invalid~=1;
            elseif removingalarms&delay==2
                toi=data_epoch_thr.trialinfo.delay2_invalid~=1;
            elseif removingalarms==0
                toi= 1:size(data_epoch_thr.trialinfo,1);
            end
            data_art_reject.eyedat=data_art_reject.eyedat(toi,:,:);
            data_art_reject.trialinfo=data_art_reject.trialinfo(toi,:);
            
        end
        
        %%
        objlist=unique(subtable.object_cued_rot);
        objlistid=unique(subtable.object_cued_id);
        
        rho=nan(size(data_art_reject.eyedat,1),3,size(data_art_reject.eyedat,3));
        
        parfor ttt=1:size(data_art_reject.trialinfo,1) %for each trial
            
            t_s_ave=[];
            item_id_rsa=nan(length(t_s_ave),3,size(data_art_reject.eyedat,3));
            D=nan(16,1);
            RRR=nan(3,size(data_art_reject.eyedat,3));
            trialid=[];
            
            %% Making averages, excluding the present trial of interest Selecting trials per item ID and rot            

            counter=1;
            for id = 1:length(objlistid)
                for k =1:length(objlist)
                    
                    actual_trial=data_art_reject.trialinfo.trial(ttt);
                    
                    if ooi==1                        
                        triales = find(data_art_reject.trialinfo.object_1_rot==objlist(k)&data_art_reject.trialinfo.object_1_id==objlistid(id)&...
                            data_art_reject.trialinfo.trial~=actual_trial); % here we group different objects                        
                    elseif ooi==2                        
                        triales = find(data_art_reject.trialinfo.object_2_rot==objlist(k)&data_art_reject.trialinfo.object_2_id==objlistid(id)&...
                            data_art_reject.trialinfo.trial~=actual_trial); % here we group different objects                        
                    end
                    [t_s_ave{counter}] = squeeze(nanmean(data_art_reject.eyedat(triales,:,:),1));                   
                    counter=counter+1;                    
                end               
            end
            %%  Just changing the format
            
            for k = 1:length(t_s_ave);
                
                if ~isempty(t_s_ave{k})
                    
                    rr=(t_s_ave{k});
                    item_id_rsa(k,:,:)=[rr];
                    
                elseif isnan(t_s_ave{k})
                    item_id_rsa(k,1:3,1:size(data_art_reject.eyedat,3))=NaN;
                end
            end
            
            %%
            
            if ooi==1
                obj1rotnum= data_art_reject.trialinfo.object_1_rot(ttt); %this is the number of the rotation for the single trial
            elseif ooi==2;
                obj1rotnum= data_art_reject.trialinfo.object_2_rot(ttt); %this is the number of the rotation for the single trial
            end

            o1=1:length(objlist);
            o2=length(objlist)+1:(length(objlist)*2);
            o3=length(objlist)*2+1:length(objlist)*3;
                        
            for p=1:(size(item_id_rsa,3))
                
                to_rsa_id = item_id_rsa(:,1:2,p);
                to_rsa_id_reshape=reshape(to_rsa_id, size(to_rsa_id,1), []);
                
                %We take the euclidean distance between the actual trial and the
                %averages
                for a1=1:size(to_rsa_id_reshape,1)
                    D(a1,1) = norm(to_rsa_id_reshape(a1,:)- data_art_reject.eyedat(ttt,1:2,p,:));
                end
                
                rot_pos=find(objlist==obj1rotnum);
                RRR(1,p)=corr(modelmat(rot_pos,:)',D(o1,1),'type','Pearson','rows','pairwise');
                RRR(2,p)=corr(modelmat(rot_pos,:)',D(o2,1),'type','Pearson','rows','pairwise');
                RRR(3,p)=corr(modelmat(rot_pos,:)',D(o3,1),'type','Pearson','rows','pairwise');                
            end
            
            rho(ttt,:,:)=RRR;
                        
        end
        
        %%
        times=data_art_reject.time;
        scriptname= mfilename('fullpath');
        trialinfo=data_art_reject.trialinfo;
        RSA_obj_rot = struct('rho',{rho},'time',{times},'trialinfo',{trialinfo},'script',{scriptname},'analysis_info',{analysis_info});
              
        %% Saving results
        fold_tosave =[mainfolder folder partis{ppp,:}];
        save (fold_tosave,'RSA_obj_rot','-v7.3')
        clearvars -except ppp partis partis_sub mainfolder fold_tosave tobefound tobefound_tr analysis_info modelmat folder retrooi toload_folder ooi conf ii delay removingalarms selectingretrocue
        
    end
end