A_SpecifyIndividualRegressorMats_withFeature.m

function A_SpecifyIndividualRegressorMats_withFeature

% Create SPM experimental condition info for each subject for later
% batching with tardis

% v3: create Dim-Att regressors, include cue onset, probe onset regressors
% use FAST modelinstead of AR(1): Corbin, N., Todd, N., Friston, K. J., & Callaghan, M. F. (2018). Accurate modeling of temporal correlations in rapidly sampled fMRI time series. Human Brain Mapping, 26(3), 839?14. http://doi.org/10.1002/hbm.24218
% - AP to change back
% v4: AP remove -1 from design timing, but add durations to cue (1 volume) and
% probe (2).

IDs = {'1117';'1118';'1120';'1124';'1125';'1126';'1131';'1132';'1135';'1136';...
    '1151';'1160';'1164';'1167';'1169';'1172';'1173';'1178';'1182';'1214';...
    '1215';'1216';'1219';'1223';'1227';'1228';'1233';'1234';'1237';'1239';'1240';'1243';...
    '1245';'1247';'1250';'1252';'1257';'1261';'1265';'1266';'1268';'1270';'1276';'1281'};

%pn.root = '/Volumes/LNDG/Projects/StateSwitch/dynamic/data/mri/task/analyses/G_GLM/';
pn.root = '/home/beegfs/kosciessa/StateSwitch/WIP/G_GLM/';

%Path to Base Dir of local output
pn.filePathRoot = '/Users/kosciessa/Desktop/beegfs/StateSwitch/WIP/G_GLM/';
pn.batchPath = [pn.filePathRoot, 'B_data/D_batch1stLevel_v2/'];

addpath([pn.filePathRoot, 'T_tools/spm12/']); % add spm functions
spm fmri

for indID = 1:numel(IDs)
    
    disp(['Processing ', IDs{indID}]);
    
    matlabbatch = cell(1);
    
    % specify general parameters
    
    matlabbatch{1}.spm.stats.fmri_spec.dir = {[pn.root, 'B_data/F_GLM1stLevelout_v1/' IDs{indID} '/']};
    matlabbatch{1}.spm.stats.fmri_spec.timing.units = 'scans';
    matlabbatch{1}.spm.stats.fmri_spec.timing.RT = 0.645;
    matlabbatch{1}.spm.stats.fmri_spec.timing.fmri_t = 16;
    matlabbatch{1}.spm.stats.fmri_spec.timing.fmri_t0 = 8;
    matlabbatch{1}.spm.stats.fmri_spec.fact = struct('name', {}, 'levels', {});
    matlabbatch{1}.spm.stats.fmri_spec.bases.hrf.derivs = [1 0];
    matlabbatch{1}.spm.stats.fmri_spec.volt = 1;
    matlabbatch{1}.spm.stats.fmri_spec.global = 'None';
    matlabbatch{1}.spm.stats.fmri_spec.mthresh = -Inf;
    matlabbatch{1}.spm.stats.fmri_spec.mask = {''};
    matlabbatch{1}.spm.stats.fmri_spec.cvi = 'FAST';
    
    % subjects with different run number
    
    if strcmp(IDs{indID}, '2131') || strcmp(IDs{indID}, '2237')
        numOfRuns = 2;
    else
        numOfRuns = 4;
    end
    
    for indRun = 1:numOfRuns
        
        % get regressor data
        pn.regressors = [pn.filePathRoot, 'B_data/A_regressors/'];
        load([pn.regressors, IDs{indID}, '_Run',num2str(indRun),'_regressors.mat'])
        
        % change regressor matrix, such that column 4 represents trial-by-trial
        % indices of dimensionality
        blockOnsets=find(Regressors(:,1) == 1);
        for indBlock = 1:numel(blockOnsets)
            if indBlock == numel(blockOnsets)
                Regressors(blockOnsets(indBlock):end,4) = ...
                    Regressors(blockOnsets(indBlock),4);
            else
                Regressors(blockOnsets(indBlock):blockOnsets(indBlock+1)-1,4) = ...
                    Regressors(blockOnsets(indBlock),4);
            end
        end
        
        % For the scans, we need to specify separate 3D images via comma seperator
        % Under the hood: design-specified time points are extracted from 2D voxel*time matrix
        basefile = [pn.root, 'B_data/BOLDin/',IDs{indID},'_run-',num2str(indRun),'_regressed.nii'];
        
        allFiles={};
        if strcmp(IDs{indID}, '2132') && indRun == 2
            disp(['Change #volumes for run2']);
            for indVol = 1:1040
                allFiles{indVol,1} = [basefile, ',', num2str(indVol)];
            end
        else
            for indVol = 1:1054
                allFiles{indVol,1} = [basefile, ',', num2str(indVol)];
            end
        end
        
        matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).scans = allFiles;
        
        % stimulus viewing condition
        % To get relevant-feature-specific conditions, we have to create 4
        % stim conditions: color, direction, size & luminance
        AttributeNames = {'1Target'; '2Targets'; '3Targets'; '4Targets'};
        for indAtt = 1:4 % here: contained in target set
            AttOnsets = find(sum(Regressors(:,7:10),2) == indAtt);
            PM = Regressors(AttOnsets,6); % add sequence position as parametric regressor
            AttDur = repmat(4, numel(AttOnsets),1);
            matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(indAtt).name = AttributeNames{indAtt};
            matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(indAtt).onset = AttOnsets; 
            matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(indAtt).duration = AttDur;
            matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(indAtt).tmod = 0;
            matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(indAtt).pmod.name = 'SequencePosition';
            matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(indAtt).pmod.param = PM;
            matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(indAtt).pmod.poly = 1;
            matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(indAtt).orth = 1;
        end
        
        % add cue regressor
        matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(5).name = 'CueOnset';
        matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(5).onset = find(Regressors(:,2) == 1); % IMPORTANT: SPM starts counting at 0
        onsets = matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(5).onset;
        matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(5).duration = repmat(3, numel(onsets), 1); clear onsets; % duration of 1 (VarToolbox) vs 0 (SPM convention)
        matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(5).tmod = 0;
        matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(5).pmod = struct('name', {}, 'param', {}, 'poly', {});
        matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(5).orth = 0;
        
        % add probe regressor with parametric RT modulation
        matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(6).name = 'ProbeOnset';
        matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(6).onset = find(Regressors(:,11) == 1); % IMPORTANT: SPM starts counting at 0
        onsets = matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(6).onset;
        matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(6).duration = repmat(2, numel(onsets), 1); clear onsets; % duration of 1 (VarToolbox) vs 0 (SPM convention)
        matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(6).tmod = 0;
        matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(6).pmod.name = 'RT';
        matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(6).pmod.param = Regressors((Regressors(:,11) == 1),13);
        matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(6).pmod.poly = 1;
        matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(6).pmod = struct('name', {}, 'param', {}, 'poly', {});
        matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(6).orth = 0;
                
        %% add 24 Friston motion regressors, and DVARS
        
        matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).multi_reg={...
            [pn.root, 'B_data/M_MotionParameters/' IDs{indID} '/' IDs{indID} '_sess-' int2str(indRun) '_motion_24dof.txt'];...
            [pn.root, 'B_data/M_MotionParameters/' IDs{indID} '/' IDs{indID} '_run-' int2str(indRun) '_feat_detrended_highpassed_denoised_dvars.txt']};
                
        %% add general session information
        matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).multi = {''};
        matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).regress = struct('name', {}, 'val', {});
        matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).hpf = 128;
        matlabbatch{1}.spm.stats.fmri_spec.fact = struct('name', {}, 'levels', {});
        matlabbatch{1}.spm.stats.fmri_spec.bases.hrf.derivs = [1 0];
        matlabbatch{1}.spm.stats.fmri_spec.volt = 1;
        matlabbatch{1}.spm.stats.fmri_spec.global = 'None';
        matlabbatch{1}.spm.stats.fmri_spec.mthresh = -Inf;
        matlabbatch{1}.spm.stats.fmri_spec.mask = {[pn.root, 'B_data/E_standards/mni_icbm152_gm_tal_nlin_sym_09c_MNI_3mm.nii']};
        matlabbatch{1}.spm.stats.fmri_spec.cvi = 'FAST';
    end % session loop (i.e. runs)
    
    %% Estimate Model
    
    matlabbatch{2}.spm.stats.fmri_est.spmmat(1) = cfg_dep('fMRI model specification: SPM.mat File', substruct('.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','spmmat'));
    matlabbatch{2}.spm.stats.fmri_est.write_residuals = 0;
    matlabbatch{2}.spm.stats.fmri_est.method.Classical = 1;

    save([pn.batchPath, IDs{indID},'_SPM1stBatchGLM.mat'], 'matlabbatch');
    
end % subject loop

end % finish completely