%% Function to compute Between Item Similarity Matrices for Encoding and Recognition separately

function BetweenItemSimilarity(subject,channels,atlasType,atlasIdx,rsawidth,rsastep,maxT,rsa_tag,roi_tag)
    %% set base path
    fs = '\';
    dataDir = '..\Data';
    ftpath = 'C:\Program Files\MATLAB\R2020b\toolbox\fieldtrip';
    save_to_disk = 1;    
        
    %load folder names of subjects
    load('..\subject_list.mat'); %subNames for easier folder accessing
    
    fileName = subNames{subject}; %set folderName of current subject
    disp(fileName)

    %% load RSA-prepped tf data for specified RSA settings
    fileTag = ['tfRSA_' num2str(rsawidth) 'width_' num2str(rsastep) 'step_' num2str(maxT) 'ms'];
    load([dataDir fs fileName fs 'tf_data' fs fileTag '_enc.mat']); %encRSA
    load([dataDir fs fileName fs 'tf_data' fs fileTag '_ret.mat']); %retRSA
    
    %% Fetch indices of channels over which to calculate the RSA
    load([dataDir fs fileName fs 'coordinates' fs 'electrodes_bp_atlas.mat'])

    chanIdx = channels;
    chanLabels = electrodes.label(chanIdx);
    
    %%% code replaced by previous channel index selection in
    %%% FindSubjOfROI_xy
    
    %% Now compute pairwise similarity matrices for encoding and retrieval separately
    encoding = zeros(length(encRSA.time),size(encRSA.trialinfo,1),size(encRSA.trialinfo,1)); % TIME X STIM X STIM
    recognition = zeros(size(encoding));

    % compute indices for pairwise stimuli 
    %get linear indices of half matrix
    confMat = ones(size(encRSA.powspctrm,1),size(encRSA.powspctrm,1)); %STIM X STIM
    linMat = tril(confMat,-1); %get only half matrix but keep diagonal
    linIdx = find(linMat); %stores linear indices for stimXstim combinations

    [r,c] = ind2sub(size(confMat,1),linIdx);
    linIdx(:,2) = r;
    linIdx(:,3) = c;

    for phase = 1:2 %1 encoding, 2 retrieval
        if phase == 1
            disp('ENCODING')
            data = encRSA;
            %temp_r_mem = data;
        elseif phase == 2
            disp('RETRIEVAL')
            data = retRSA;
            %temp_r_mem = data;
        end
         
        for win = 1:size(encoding,1)
            for ix = 1:length(linIdx)
                stim1_data = permute(squeeze(data.powspctrm(linIdx(ix,2),chanIdx,:,win)),[3 1 2]); %win, chan freq
                stim1 = squeeze(reshape(stim1_data,size(stim1_data,1),[]))';
                stim2_data = permute(squeeze(data.powspctrm(linIdx(ix,3),chanIdx,:,win)),[3 1 2]); %win, chan freq
                stim2 = squeeze(reshape(stim2_data,size(stim2_data,1),[]))';
                
                
                %[r,p]= corr(stim1,stim2,'rows','complete','type','spearman'); %win*win for each stimpair
                if phase == 1
                    encoding(win,linIdx(ix,2),linIdx(ix,3)) = corr(stim1, stim2, 'rows','complete','type','spearman');
                elseif phase == 2
                    recognition(win,linIdx(ix,2),linIdx(ix,3)) = corr(stim1, stim2, 'rows','complete','type','spearman');
                end
            end

            %normalize and make full matrix!                   
            if phase == 1
                mat = squeeze(encoding(win,:,:)); %half matrix w/o diagonal
                lower_vals = mat(mat~=0); %vector of half matrix
                FishZ_vals = 0.5.*log((1+lower_vals) ./ (1-lower_vals)); %FisherZ transform half matrix
                mat(mat~=0) = FishZ_vals;% set as new half matrix
                encoding(win,:,:) = mat;
            elseif phase == 2
                mat = squeeze(recognition(win,:,:)); %half matrix with diagonal 1
                lower_vals = mat(mat~=0); %vector of half matrix
                FishZ_vals = 0.5.*log((1+lower_vals) ./ (1-lower_vals)); %FisherZ transform half matrix
                mat(mat~=0) = FishZ_vals;% set as new half matrix
                recognition(win,:,:) = mat;
            end
        end
    end
    
    %%BtwInfo file
    info = {};
    info.atlasType = atlasType;
    info.areaIdx = atlasIdx;
    info.subject = subject;
    info.channels = chanIdx;
    info.chanLabels = chanLabels;
    info.trialinfo = encRSA.trialinfo;

    %% set saving directory and save
    saveDir = ['stats' fs 'BtwItem' fs rsa_tag fs roi_tag];
    
    status = mkdir([dataDir fs fileName fs saveDir]);
    save([dataDir fs fileName fs saveDir fs 'BtwItem_encoding.mat'],'encoding')
    save([dataDir fs fileName fs saveDir fs 'BtwItem_recognition.mat'],'recognition')
    save([dataDir fs fileName fs saveDir fs 'BtwInfo.mat'],'info')
                
end
