% This script runs a multi-voxel pattern similarity analysis. It reads in
% the output of a single-trial model, extracts multi-voxel patterns for
% each trial, and summarizes similarity across designated pairs of trials.
%
% The script assumes that a single-trial model has been estimated for each
% subject, and that the files are organized into a hierarchy of folders in
% the following format: [stDir]/s01/Run1/trial_001/, where each
% trial folder contains the output of the model corresponding to that
% trial. Trial information is also coded into the regressor name
% stored in the SPM.mat file, e.g., run2tri50objID352day2conID6 etc.
%
% The script initializes parameters that are fed into a series of
% functions that actually run the analysis. The functions are as follows:
%
% 1. Functions used to initialize file structure & information:
%
% CONCAT_BETAS_ST: This function copies beta_0001.img/hdr files from each
% trial folder to a single "betas" directory for each subject. Files are
% renamed according to their run/trial number, e.g., Run1_trial001.img/hdr.
% A trialinfo.mat file containing the regressor names will also be saved.
% It is unlikely that you will need to run this function more than once.
%
% DECODE_TRIALINFO: This function decodes the trialinfo.mat file into an
% easier to manage format. It also allows you to create new condition codes
% by taking the factorial of any two conditions.
%
% 2. Functions used to compute the pattern similarity matrix for an ROI:
%
% GET_TRIALMATS_INDIVROI: This function extracts the single-trial patterns
% from each trial for each ROI, and computes the full correlation
% matrix linking each trial to every other trial. Trial matrices (pattsim)
% are saved in .mat files for every ROI for every subject.
%
% 3. Functions used to pre-process the pattern similarity matrix:
%
% FIND_BAD_TRIALS: This function identifies "bad" trials based on deviant
% values across all voxels in each volume. Trial IDs are saved for future
% filtering.
%
% CLEAN_UP_PATTSIM: This function excludes trials from the pattern
% similarity analysis, based on deviant signal (ROI-specific or
% whole-brain), and excludes similarity values for pairs of temporally
% adjacent trials.
%
% 4. Functions used to summarize pattern similarity:
%
% SUMMARIZE_TRIALMATS: This function summarizes the average pattern
% similarity for designated trial pairs, which are specified on the basis
% of their shared trial information, separately for each subject and ROI.
% Summary data are saved to the group directory, and the entire group is
% typically run in the same pass.
%
%
% Author: Maureen Ritchey, 2012-2014
%
% For the paper: "Delay-dependent contributions of medial temporal regions 
% to episodic memory retrieval" submitted to eLife Jan 2015

function [filename, ps] = pattern_batch()

curdir = pwd;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Main

% Functions to run
runmeplease = {'concat' 'decode'}; % FIRST - initialize file structure & info
% runmeplease = {'get_indiv'}; % SECOND - get pattern similarity matrix for each ROI
% runmeplease = {'find_bad_trials' 'cleanup' }; % THIRD - pattern similarity 'pre-processing'
% runmeplease = {'summarize'}; % whole-group % FOURTH - summarize pattern similarity

% Single trial directory.
% Assumes that there is a subfolder for each subject.
ps.stDir = '/Users/mritchey/maureen/consim/Analysis/ObjRec_SingleTrial_f/';
ps.subjID = {'s03' 's04' 's05' 's06' 's10' 's12' 's15' 's16' 's17' 's18' 's21' 's22' 's23' 's26' 's28' 's29' 's31' 's32' 's35'};

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Options for concat_betas_st

% Run names in single-trial folder
ps.runs = {'Run1' 'Run2' 'Run3'};

% Which beta to use (usually the first one, beta_0001) & label to assign
ps.betaID = 'beta_0001';
ps.triFilter = 'betamap_';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Options for decode_trialinfo

% Specify which conditions should be coded
ps.allcons = {'run' 'tri' 'objID' 'day' 'conID' 'pairID' 'numresp' 'dm'};

% Use ps.vals and ps.newLabel to generate a new factorial trial identifier
ps.vals = [0 0 0 2 0 0 6 0];
ps.newLabel = 'dayxnumresp';

% Set ps.startOver to 1 if you wish to completely regenerate trialinfo_num
ps.startOver = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% ROI options

% options for using individually-defined ROIs (get_indiv)
ps.Pdir = '/Users/mritchey/maureen/consim/Data/mprage/'; %assumes there are subject directories inside, and that masks are named ['rs01_' ps.Pname] etc
ps.Pname = {'HIPP_ANT_L_filtered_1SD.nii.gz' 'PRC_L_filtered_1SD.nii.gz' 'PHC_L_filtered_1SD.nii.gz' 'HIPP_POST_L_filtered_1SD.nii.gz'  ...
    'HIPP_ANT_R_filtered_1SD.nii.gz' 'PRC_R_filtered_1SD.nii.gz' 'PHC_R_filtered_1SD.nii.gz' 'HIPP_POST_R_filtered_1SD.nii.gz'
    };

% Summary label
ps.sum_label = 'all_mtl_anat_betas'; % append "corrected" for fourth stage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Options for get_trialmats and get_trialmats_indivroi

% Flag to zscore patterns within a trial (1) or across trials (2)
ps.do_zscore = 1;

% Specify which distance metric to use: 1 for Pearson's correlation and 2
% for Euclidean distance
ps.dist_metric = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Options for find_bad_trials & clean_up_pattsim

% Define mask to use for find_bad_trials (uses wholebrain mask from
% single-trial model if left empty)
ps.wbmaskname =  ''; %'gm_func_mask.nii';

% Use histograms to determine threshold for determining outlying signal
% values - generates exclusions fiel
ps.do_hist = 0;
ps.badthresh = 1.2; %seems to best capture 'tail'

% Choose whether to filter out outlying trials
%1 for filter based on current ROI only, 2 for exclusions file only, 3 for both, 0 for no-filter
ps.do_filter = 2; 
ps.exclusionsFile = [ps.triFilter 'bad_trials_1.2.mat'];

% Choose whether to exclude comparisons between immediately adjacent trials
ps.exc_adj = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Options for summarize_trialmats (pattern similarity analysis)

% Specify which factor to limit by (leave empty if do not want to limit trials)
ps.limitcons = {'dayxnumresp' [6 12]}; % limit to day 1 recollection and day 2 recollection
%  ps.limitcons = {'dayxnumresp' [4 5 10 11]}; % limit to day 1 familiarity and day 2 familiarity

% Label for output
ps.out_label = 'RFX_REC';

% Specify which factor to summarize by (group pairs of trials)
ps.cons = 'conID'; 

% Choose whether to generate and print figures
ps.printimg = 0;
ps.write_csv = 1;

% Pair names
%%% same=same room, rel=grouped rooms, diff=ungrouped rooms, dif_day=between day
ps.allset_names = {'del_same' 'del_rel' 'del_diff' 'imm_same' 'imm_rel' 'imm_diff' 'dif_day'};

% Symmetric pairs of interest
%%% cons 1-4 correspond to day 1 rooms; cons 5-8 correspond to day 2 rooms
ps.allsets = {[1 1; 2 2; 3 3; 4 4] ...
    [1 2; 2 1; 3 4; 4 3] ...
    [1 3; 3 1; 1 4; 4 1; 2 3; 3 2; 2 4; 4 2] ...
    [5 5;6 6;7 7;8 8] ...
    [5 6; 6 5;7 8; 8 7] ...
    [5 7; 7 5; 5 8; 8 5; 6 7; 7 6; 6 8; 8 6] ...
    [1 5; 5 1; 1 6; 6 1; 1 7; 7 1; 1 8; 8 1; 2 5; 5 2; 2 6; 6 2; 2 7; ...
    7 2; 2 8; 8 2; 3 5; 5 3; 3 6; 6 3; 3 7; 7 3; 3 8; 8 3; 4 5; 5 4; ...
    4 6; 6 4; 4 7; 7 4; 4 8; 8 4]};

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Save the workspace
if ~exist([ps.stDir 'batch_spec'],'dir')
    mkdir([ps.stDir 'batch_spec']);
end
filename = [ps.stDir 'batch_spec/pattern_batch_' datestr(now,'yyyy-mm-dd_HH-MM')];
fprintf('\n--------------------\n\nSaving batch spec to: %s\n',filename);
save(filename);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Run functions

% 1. Functions used to initialize file structure & information:

if max(strcmp('concat',runmeplease))
    fprintf('\n--------------------\n\nRunning concat_betas_st...\n');
    concat_betas_st(ps);
    cd(curdir);
end

if max(strcmp('decode',runmeplease))
    fprintf('\n--------------------\n\nRunning decode_trialinfo...\n');
    decode_trialinfo(ps);
    cd(curdir);
end

% 2. Functions used to compute the pattern similarity matrix for an ROI:

if max(strcmp('get_indiv',runmeplease))
    fprintf('\n--------------------\n\nRunning get_trialmats...\n');
    get_trialmats_indivroi(ps);
    cd(curdir);
end

% 3. Functions used to pre-process the pattern similarity matrix:

if max(strcmp('find_bad_trials',runmeplease))
    fprintf('\n--------------------\n\nRunning find_bad_trials...\n');
    find_bad_trials(ps);
    cd(curdir);
end

if max(strcmp('cleanup',runmeplease))
    fprintf('\n--------------------\n\nRunning clean_up_pattsim...\n');
    clean_up_pattsim(ps);
    cd(curdir);
end

% 4. Functions used to summarize pattern similarity:

if max(strcmp('summarize',runmeplease))
    fprintf('\n--------------------\n\nRunning summarize_trialmats...\n');
    summarize_trialmats(ps);
    cd(curdir);
end

end %main function


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [] = concat_betas_st(ps)

%%%Copy and rename betas to new folder
for iSub=1:length(ps.subjID)
    fprintf('Working on Subject %s\n',ps.subjID{iSub})
    
    %initialize variables and directories
    counter = 1;
    trialinfo = {};
    if ~exist([ps.stDir ps.subjID{iSub} '/betas/'])
        mkdir([ps.stDir ps.subjID{iSub} '/betas/']);
    end
    
    %loop over ps.runs
    for iRun = 1:length(ps.runs)
        fprintf('Working on Run %d\n',iRun)
        curDir = [ps.stDir ps.subjID{iSub} '/' ps.runs{iRun} '/'];
        cd(curDir)
        
        %loop over trials in run
        tri = dir('trial_*');
        for iTri = 1:length(tri)
            
            try
                %get & store the current regressor name, which contains info about the trial
                regs = load([curDir tri(iTri).name '/regs.mat']);
                trialinfo = [trialinfo; regs.names{1}];
                
                %get current beta image, rename with run and trial info, & copy to single directory
                curImage = [curDir tri(iTri).name '/' ps.betaID '.img'];
                if exist(curImage,'file')
                    destImage = [ps.stDir ps.subjID{iSub} '/betas/' ps.triFilter ps.runs{iRun} '_' tri(iTri).name '.img'];
                    copyfile(curImage,destImage);
                    curImage = [curDir tri(iTri).name '/' ps.betaID '.hdr'];
                    destImage = [ps.stDir ps.subjID{iSub} '/betas/' ps.triFilter ps.runs{iRun} '_' tri(iTri).name '.hdr'];
                    copyfile(curImage,destImage);
                end
                clear regs
                
            catch err
                fprintf('Error on Subject %s, Run %d, Tri %d\n',ps.subjID{iSub},iRun,iTri);
                continue;
            end
            
        end
    end
    
    %save trial information for future use
    save([ps.stDir ps.subjID{iSub} '/trialinfo.mat'],'trialinfo');
end

end %concat_betas_st

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [] = decode_trialinfo(ps)

for iSub = 1:length(ps.subjID)
    curDir = [ps.stDir ps.subjID{iSub} '/'];
    
    if ~exist([curDir 'trialinfo_num.mat']) || ps.startOver
        %load trial info
        curInfo = [curDir 'trialinfo.mat'];
        load(curInfo)
        labels = ps.cons;
        
        %extract category information stored in trialinfo (originally embedded in regressor names)
        for iTri = 1:length(trialinfo)
            for iCon = 1:length(ps.cons)
                i = strfind(trialinfo{iTri},ps.cons{iCon}) + length(ps.cons{iCon});
                if iCon < length(ps.cons)
                    j = strfind(trialinfo{iTri},ps.cons{iCon+1}) - 1;
                    cat(iTri,iCon) = str2num(trialinfo{iTri}(i:j));
                else
                    cat(iTri,iCon) = str2num(trialinfo{iTri}(i:end));
                end
            end
        end
    else
        %if this function has already been run, can run again to define new column
        curInfo = [curDir 'trialinfo_num.mat'];
        load(curInfo)
        cat = trialinfo;
        if ~isempty(ps.newLabel)
            if max(cell2mat(strfind(labels,ps.newLabel)))>0
                error('This label has been used already.');
            end
        end
    end
    
    %create new columns-- NOTE: only works for 2 factors so far
    new = [];
    if ~isempty(ps.vals)
        new = zeros(size(cat,1),1);
        counter = 1;
        tmp = find(ps.vals>0);
        j = tmp(1); k = tmp(2);
        
        for jx = 1:ps.vals(j)
            for kx = 1:ps.vals(k)
                new(cat(:,j)==jx & cat(:,k)==kx) = counter;
                counter = counter + 1;
            end
        end
    end
    
    %add in new column and save
    trialinfo = [];
    trialinfo = [cat new];
    labels = [labels ps.newLabel];
    filename =[curDir 'trialinfo_num.mat'];
    save(filename,'trialinfo','labels');
    
end
end %decode_trialinfo

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [] = get_trialmats_indivroi(ps)

for iSub = 1:length(ps.subjID)
    fprintf('\nWorking on %s ...\n', ps.subjID{iSub})
    curDir = [ps.stDir ps.subjID{iSub} '/'];
    
    %loop over ROIs
    for nROI=1:length(ps.Pname)
        
        %define P based on ps.Pdir and ps.Pname for current ROI
        P = [];
        P = [ps.Pdir,ps.subjID{iSub},'/r',ps.subjID{iSub},'_',ps.Pname{nROI}];
        [a b c] = fileparts(ps.Pname{nROI});
        
        %if it's zipped, unzip it & correct filename
        if strfind(c,'.gz')
            gunzip(P);
            P = strrep(P,'.gz','');
        end
        
        %use as ROI mask
        M=spm_vol(P);
        mask = spm_read_vols(M);
        mask = round(mask); % make sure it's binary
        
        %load beta patterns for each trial
        fprintf('- Extracting ROI beta series for %s.\n', b)
        [t B] = getAllFiles_filtered([curDir 'betas/'],'img',ps.triFilter);
        numvox = length(find(mask>0));
        % B = B(1:20); %NOTE: I USE THIS FOR TESTING ONLY
        V=spm_vol(char(B));
        for j=1:size(V,1)
            beta=spm_read_vols(V(j));
            beta=beta(mask>0);
            tmp(:,j) = reshape(beta,numel(beta),1);
            dispcounter(j);
        end
        
        %get rid of any rows (voxels) containing a NaN
        patt = []; filt_v = [];
        for n=1:size(tmp,1)
            if max(isnan(tmp(n,:)))==0
                patt = [patt; tmp(n,:)];
            else
                filt_v = [filt_v n];
            end
        end
        fprintf('\n - Filtered %i voxels.',length(filt_v));
        clear tmp;
        
        %identify columns (trials) where mean activation is deviant but
        %don't filter til clean_up_pattsim if requested
        pattmean = mean(patt,1);
        zpattmean = zscore(pattmean');
        filt_t = find((abs(zpattmean)>3))';
        fprintf('\n - Identified %i timepoints for future filtering.',length(filt_t));
        
        %if requested, normalize along each column (trial) - doesn't make a
        %difference if using correlation but would for Euc dist
        if ps.do_zscore
            for n=1:size(patt,2)
                patt(:,n) = (patt(:,n)-mean(patt(:,n)))./std(patt(:,n));
            end
        end
        
        %calculate similarity matrix
        fprintf('\n - Calculating similarity matrices.')
        if ps.dist_metric==1
            pattsim = corrcoef(patt);
        elseif ps.dist_metric==2
            pattsim = pdist(patt);
        else
            error('Specify which distance metric you prefer.');
        end
        
        %Save out distance values
        if ~exist(strcat(curDir,ps.sum_label),'dir')
            mkdir(strcat(curDir,ps.sum_label));
        end
        filename = strcat(curDir,ps.sum_label,'/',b,'_trialmat.mat');
        fprintf('\n - Writing file %s\n\n',filename);
        save(filename,'pattsim','filt_t','filt_v','ps.do_zscore','ps.dist_metric','numvox');
        
        %if mask was zipped before, rezip it
        if strfind(c,'.gz')
            gzip(P);
            delete(P);
        end
        
        clear pattsim patt filt_t pattmean zpattmean
        
    end %nROI
    
end %subjects

end %get_trialmats_indivroi

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [] = find_bad_trials(ps)

% This function identifies bad trials based on deviant values across
% within-brain voxels. Deviant values are defined by absolute z-scores
% across the trial series, then averaged across all voxels. Trials with an
% average absolute z-score greater than 'thresh' are identified for future
% filtering and saved in bad_trials*.mat.

for iSub = 1:length(ps.subjID)
    
    fprintf('\nWorking on %s.\n', ps.subjID{iSub})
    curDir = [ps.stDir ps.subjID{iSub} '/'];
    
    % if maskname is empty, it will automatically grab single-trial model mask
    % if maskname doesn't include a path, it will automatically look in each
    % subject's folder in ps.stDir
    if ~isempty(ps.wbmaskname)
        newmaskname = strcat(ps.stDir,ps.subjID{iSub},'/',ps.wbmaskname);
    else
        newmaskname = ps.wbmaskname;
    end
    if isempty(newmaskname)
        brainmask = [ps.stDir ps.subjID{iSub} '/Run1/trial_001/mask.img'];
    else
        brainmask = newmaskname;
    end
    fprintf('Using mask: %s\n',brainmask);
    M=spm_vol(brainmask);
    mask = spm_read_vols(M);
    
    %load beta patterns for each trial
    [t B] = getAllFiles_filtered([curDir 'betas'],'img',ps.triFilter);
    V=spm_vol(char(B));
    for iTri = 1:size(V,1)
        dispcounter(iTri);
        tmp=spm_read_vols(V(iTri)); % read in data and apply mask
        tmp=tmp(mask>0);
        allbeta(iTri,:) = tmp';
        clear tmp
    end
    
    %get zscores & average across voxels
    zallbeta = abs(zscore(allbeta));
    zavgbeta = nanmean(zallbeta,2);
    
    %plot histogram of zscores
    %visual inspection used to determine threshold for 'tail'
    if ps.do_hist
        figure;
        hist(zavgbeta);
        title(ps.subjID{iSub});
    end
    
    %label trials that exceed threshold
    bad_trials = find(zavgbeta>ps.badthresh);
    fprintf('\nIdentified %0.0f bad trials with %0.2f threshold\n\n',length(bad_trials),ps.badthresh);
    
    info_filename = [curDir ps.triFilter 'bad_trials_' num2str(ps.badthresh) '.mat'];
    save(info_filename,'bad_trials');
    
    clear mask M allbeta V B brainmask zallbeta zavgbeta bad_trials
    
end

end %find_bad_trials

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [] = clean_up_pattsim(ps)

for iSub = 1:length(ps.subjID)
    
    fprintf('\nWorking on %s\n',ps.subjID{iSub})
    
    %load matrix coding adjacent voxels - could also use this to exclude other trials a priori
    corrections = {};
    curInfo = [ps.stDir 'adjacentvox2.mat'];
    load(curInfo)
    corrections = [corrections curInfo];
    
    P = ps.Pname;
    for nROI=1:length(P)
        
        dispcounter(nROI);
        
        %load pattsim
        curDir = [ps.stDir ps.subjID{iSub} '/'];
        [a b c] = fileparts(P{nROI});
        b = strrep(b,'.gz','''');
        if ps.do_zscore == 1
            curMat = [curDir ps.sum_label '/' b '_trialmat.mat'];
        elseif ps.do_zscore == 2
            curMat = [curDir ps.sum_label '/' b '_trialmat_z.mat'];
        end
        fprintf('\nROI: %s\n',b);
        
        load(curMat)
        
        %make sure pattsim has the same dimensions as the adjacent voxels matrix
        if size(pattsim,1) ~= size(adjvox,1)
            error('pattsim does not have the same dimensions as adjvox. Check get_trialmats.');
        end
        
        %exclude adjacent trials
        if ps.exc_adj==1
            pattsim(adjvox==1)=NaN;
            fprintf('Excluded adjacent +- 2 trials\n');
        end
        
        %filter out bad trials
        %find exclusions, if they exist
        if exist([curDir ps.exclusionsFile],'file') %for user-defined exclusions (e.g., motion ps.runs)
            exclusions = load([curDir ps.exclusionsFile]);
        else
            exclusions.bad_trials = [];
        end
        
        if ps.do_filter==1 %filter out filt_t (based on ROI outliers) only
            pattsim(filt_t,:) = NaN;
            pattsim(:,filt_t) = NaN;
            fprintf('Filtered %0.0f trials based on ROI outliers\n',length(filt_t));
        elseif ps.do_filter==2 %filter out based on bad_trials.mat only
            pattsim(exclusions.bad_trials,:) = NaN;
            pattsim(:,exclusions.bad_trials) = NaN;
            fprintf('Filtered %0.0f trials based on %s\n',length(exclusions.bad_trials),ps.exclusionsFile);
        elseif ps.do_filter==3 %filter out both
            filt_t2 = unique([filt_t exclusions.bad_trials]);
            pattsim(filt_t2,:) = NaN;
            pattsim(:,filt_t2) = NaN;
            fprintf('Filtered %0.0f trials based on ROI + %s\n',length(filt_t2),ps.exclusionsFile);
        else %don't do any filtering
            fprintf('No trial filtering applied.\n');
            
        end
        
        %save new trialmat file appended with "corrected" label
        [aa bb cc] = fileparts(curMat);
        if ~exist([aa '_corrected' num2str(ps.do_filter)],'dir')
            mkdir([aa '_corrected' num2str(ps.do_filter)]);
        end
        
        newMat = [aa '_corrected' num2str(ps.do_filter) '/' bb cc];
        save(newMat,'pattsim','filt_t','filt_v','ps.do_zscore','ps.dist_metric','corrections');
        
        clear pattsim filt_t filt_v
        
    end
end

end %clean_up_pattsim

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [] = summarize_trialmats(ps)

%initizalize cell array for csv output
outputdata = {'subject' 'ROI' 'similarity' 'set.type' 'sum.label' 'out.label' 'limit.bin'};
counter = 2;

%loop over subjects
for iSub = 1:length(ps.subjID)
    
    fprintf('\nWorking on %s\n',ps.subjID{iSub})
    
    %load trialinfo for current subject
    curDir = [ps.stDir ps.subjID{iSub} '/'];
    curInfo = [curDir 'trialinfo_num.mat'];
    load(curInfo)
    
    % loop over ROIs
    P = ps.Pname;
    for nROI=1:length(P)
        
        dispcounter(nROI);
        
        %load pattsim - the trial similarity matrix for the current ROI
        [a b c] = fileparts(P{nROI});
        b = strrep(b,'.gz','''');
        if ps.do_zscore == 1
            curMat = [curDir ps.sum_label '/' b '_trialmat.mat'];
        elseif ps.do_zscore == 2
            curMat = [curDir ps.sum_label '/' b '_trialmat_z.mat'];
        end
        load(curMat)
        
        %make sure pattsim has the same dimensions as trialinfo
        if size(pattsim,1) ~= size(trialinfo,1)
            error('pattsim does not have the same dimensions as trialinfo. Check get_trialmats.');
        end
        
        %calculate similarity within each run
        for iRun = unique(trialinfo(:,1))'; %run number stored in first column
            
            %find con label for summarizing pairs
            iCon = [];
            tmp = strcmp(labels,ps.cons);
            iCon = find(tmp);
            if isempty(iCon)
                error('cons label not found. If necessary, run decode_trialinfo.');
            end
            
            %find column number for limiting contrast (if specified)
            if ~isempty(ps.limitcons)
                iLimCon = [];
                tmp = strcmp(labels,ps.limitcons{1,1});
                iLimCon = find(tmp);
                if isempty(iLimCon)
                    error('cons label not found for limiting contrast. If necessary, run decode_trialinfo.');
                end
            end
            
            %limit pattsim to only lower half (no double counting, but
            %pairs must be specified symmetrically)
            pattsim(tril(pattsim)==0) = NaN;
            pattsim(eye(size(pattsim))==1) = NaN;
            
            %find means of each set: first, find all pairwise combinations
            %that contribute to the set. then, find the mean pattern
            %similarity across all of those pairs.
            for iset = 1:length(ps.allsets)
                
                allxy = []; % refresh for current set
                
                % get the current trial combination for the current set
                for icombo = 1:size(ps.allsets{iset},1)
                    x = []; y = []; % refresh for current combination
                    
                    i = ps.allsets{iset}(icombo,1);
                    j = ps.allsets{iset}(icombo,2);
                    
                    if ~isempty(ps.limitcons)
                        x = find(trialinfo(:,iCon)==i & ismember(trialinfo(:,iLimCon),ps.limitcons{1,2}) & trialinfo(:,1)==iRun);
                        y = find(trialinfo(:,iCon)==j & ismember(trialinfo(:,iLimCon),ps.limitcons{1,2}) & trialinfo(:,1)==iRun);
                    else
                        x = find(trialinfo(:,iCon)==i & trialinfo(:,1)==iRun);
                        y = find(trialinfo(:,iCon)==j & trialinfo(:,1)==iRun);
                    end
                    
                    % get all pattern similarity values for this
                    % combination
                    curxy = [];
                    for ix = 1:length(x)
                        for iy = 1:length(y)
                            curxy = [curxy; pattsim(x(ix),y(iy))];
                        end
                    end
                    
                    % tack on these values to entire set
                    allxy = [allxy; curxy];
                end %icombo
                
                % take mean of all values for this set
                if ~isempty(allxy) & find(~isnan(allxy))
                    %make sure both trial types are found
                    tmpsum(iset,iRun) = nanmean(allxy);
                else
                    tmpsum(iset,iRun) = NaN;
                end
                
            end %iset
            
        end %run
        
        % average within-run similarity across ps.runs & store in roi variable
        curtmpsum = nanmean(tmpsum,2);
        roi(nROI).summary(:,iSub) = curtmpsum;
        clear tmpsum
        
        % store labels & analysis information in roi variable
        roi(nROI).roiname = b;
        roi(nROI).roilabel = b; % could have shorter label
        roi(nROI).cons = ps.cons;
        roi(nROI).subjects = ps.subjID;
        
        limitbin = [ps.limitcons{1,1} num2str(ps.limitcons{1,2})];
        limitbin = strrep(limitbin,' ','');
        
        % write to cell array that will be exported for summary stats
        for setid = 1:size(curtmpsum,1)
            
            outputdata{counter,1} = ps.subjID{iSub};
            outputdata{counter,2} = b;
            outputdata{counter,3} = curtmpsum(setid);
            outputdata{counter,4} = ps.allset_names{setid};
            outputdata{counter,5} = ps.sum_label;
            outputdata{counter,6} = ps.out_label;
            outputdata{counter,7} = limitbin;
            
            counter = counter + 1;
        end
        
    end
end

% create output directory
if ~exist([ps.stDir 'group/' ps.sum_label],'dir')
    mkdir([ps.stDir 'group/' ps.sum_label]);
end

filename = [ps.stDir 'group/' ps.sum_label '/' ps.sum_label '_' ps.out_label '_summary'];

% write out csv file for summary stats
if ps.write_csv
    cell2csv([filename '.csv'],outputdata);
end

% save .mat file
save(filename, 'roi');

end %summarize_trialmats
