function[job] = spm_bsm(job)
%
% Copyright J.B.C. Marsman
% NeuroImaging Center 
% University Medical Center Groningnen
%
% September 2015
%
% Cisler 2014
% Rissman 2004
% Mumford 2012

% Camara 2008 (reward processing)

display('== == == Beta Series Method: Running analysis == == ==')

N = numel(job.spms);

warning('OFF', 'MATLAB:MKDIR:DirectoryExists');

for spm_i = 1:N
    spmfile = job.spms{spm_i};
    %% load original SPM file
    load(spmfile);

    %% check for canonical HRF :
    if (SPM.xBF.order ~=1 )
        warning([sprintf('Problem found in File : %s\n', spmfile) ...
                'No canonical HRF was used to model the data\n' ...
                'This model will be skipped']);
        continue;
    end
    
	if (length(SPM.Sess) > 1 )
        warning([sprintf('Problem found in File : %s\n', spmfile) ...
                         sprintf('Multiple sessions are not supported.\n') ...
                         sprintf('This model will be skipped\n')]);
        continue;
    end
    
    SPM_orig = SPM;
    U_orig = SPM.Sess.U;

    target_folder = [fileparts(spmfile) filesep 'BSM'];
 
    %% create output folder
    mkdir(target_folder);
    
    %% specify active columns:
    
    all_columns = 1:length(SPM.Sess(1).U);
    if ~isempty(job.acfl)
        passive_columns = eval(['all_columns( job.acfl );']);      
    else
        passive_columns = [];
    end
    active_columns = setdiff(all_columns, passive_columns);
    
    
    job.active_columns = active_columns;
    job.SPM_orig = SPM_orig;
    job.U_orig = U_orig;
    job.target_folder = target_folder;
    
    %% loop for active columns:
        
    methods = {'Traditional BSM Rissman et al', 'Adapted BSM Rissman et al', 'Fast approach (Mumford et al)', 'Adapted fast approach','none'};
    display(sprintf('== == == Beta Series Method: Analysis Type : %s == == ==', methods{job.type}));

    switch(job.type)
        case 1, %% Traditional BSM Rissman et al
            BSM_type1(job);
            
        case 2, %% Adapted Rissman et al
            BSM_type2(job);
            
        case 3, %% Fast approach (Mumford et al)
            BSM_type3(job);
            
        case 4,  %% Adapted Fast approach
            BSM_type4(job);
            
    end 
    pwd
    keyboard
    save('IntermediateState','job');
    display('== == == Beta Series Method : Concatenating betas == == ==')
    BSM_concat_betas(job)
    return; %for now don't do the last few steps
    display('== == == Beta Series Method : Regressing ROIs== == ==')
    BSM_Regress_Regions(job);
    display('== == == Beta Series Method : Contrasting second level betas == == ==')
    BSM_Contrast_manager(job);
    display('== == == Beta Series Method : Done == == ==')
end

end


function[] = BSM_type1(job)
%% Traditional BSM by Rissman et al 2004

active_columns = job.active_columns;


for c = 1:length(active_columns) %% split out all columns but baseline
    column = active_columns(c);
    other_columns = setdiff(1:length(active_columns), column);
    
    display(sprintf('*) Exploding column %d...', column));
    
    target_folder_c = [job.target_folder filesep sprintf('col_%04d', column)];
    mkdir(target_folder_c);
    
    pos_index =1;
    
    U_new = struct('name',[],'ons', [], 'dur',[], 'P',[], 'dt', [], 'u', [],'pst',[]);
    
    for i = 1:length(job.U_orig(column).ons) %% walk over all onsets within 'column'
        
        SPM = job.SPM_orig; % copy settings from original SPM model
        target_folder_i = [target_folder_c filesep sprintf('trial_%04d', i)];
        mkdir(target_folder_i);
        SPM.swd = target_folder_i;
        
        %% add event i from 'column'
        U_new(pos_index).name{1} = sprintf('%s_%d', SPM.Sess.U(column).name{1} , i);
        U_new(pos_index).ons  =  job.SPM_orig.Sess.U(column).ons(i) ;
        U_new(pos_index).dur  = job.SPM_orig.Sess.U(column).dur(i);
        U_new(pos_index).P = job.SPM_orig.Sess.U(column).P;
        U_new(pos_index).dt = job.SPM_orig.Sess.U(column).dt;
        
        %% add all other events but i from 'column' to next column
        pos_index = pos_index +1;
        other_trials = setdiff(1:length(job.U_orig(column).ons), i);
        
        U_new(pos_index).name{1} = sprintf('%s_oc', SPM.Sess.U(column).name{1});
        U_new(pos_index).ons  =  job.SPM_orig.Sess.U(column).ons(other_trials) ;
        U_new(pos_index).dur  = job.SPM_orig.Sess.U(column).dur(other_trials);
        U_new(pos_index).P = job.SPM_orig.Sess.U(column).P;
        U_new(pos_index).dt = job.SPM_orig.Sess.U(column).dt;
        
        pos_index = pos_index + 1;
        
        for oc = 1:length(other_columns)
            other_column = other_columns(oc);
            
            U_new(pos_index).name{1} = 'other conds';
            U_new(pos_index).ons  =  [U_new(pos_index).ons job.SPM_orig.Sess.U(other_column).ons];
            
            if length(job.SPM_orig.Sess.U(other_column).dur) == 1
                U_new(pos_index).dur  = [U_new(pos_index).dur repmat(job.SPM_orig.Sess.U(other_column).dur, 1, length(job.SPM_orig.Sess.U(other_column).ons))];
            else
                U_new(pos_index).dur  = [U_new(pos_index).dur job.SPM_orig.Sess.U(other_column).dur];
            end
            U_new(pos_index).P = job.SPM_orig.Sess.U(other_column).P;
            U_new(pos_index).dt = job.SPM_orig.Sess.U(other_column).dt;
        end
        
        %% store the new U structure
        SPM.Sess.U = U_new;
        
        BSM_finalize_model(job, SPM);
        
    end %% end loop for trial i    
end %% end loop for active columns
end

function BSM_type2(job)
%% Adapted Rissman et al:
active_columns = job.active_columns;
for c = 1:length(active_columns) %% split out all columns but baseline
    column = active_columns(c);
    other_columns = setdiff(1:length(active_columns), column);
    
    display(sprintf('*) Exploding column %d...', column));
    
    target_folder_c = [job.target_folder filesep sprintf('col_%04d', column)];
    mkdir(target_folder_c);
    
    pos_index =1;
    
    U_new = struct('name',[],'ons', [], 'dur',[], 'P',[], 'dt', [], 'u', [],'pst',[]);
    
    for i = 1:length(job.U_orig(column).ons) %% walk over all onsets within 'column'
        
        SPM = job.SPM_orig; % copy settings from original SPM model
        target_folder_i = [target_folder_c filesep sprintf('trial_%04d', i)];
        mkdir(target_folder_i);
        SPM.swd = target_folder_i;
        
        %% add event i from 'column'
        U_new(pos_index).name{1} = sprintf('%s_%d', SPM.Sess.U(column).name{1} , i);
        U_new(pos_index).ons  =  job.SPM_orig.Sess.U(column).ons(i) ;
        U_new(pos_index).dur  = job.SPM_orig.Sess.U(column).dur(i);
        U_new(pos_index).P = job.SPM_orig.Sess.U(column).P;
        U_new(pos_index).dt = job.SPM_orig.Sess.U(column).dt;
        
        %% add all other events but i from 'column' to next column
        pos_index = pos_index +1;
        other_trials = setdiff(1:length(job.U_orig(column).ons), i);
        
        U_new(pos_index).name{1} = sprintf('%s_oc', SPM.Sess.U(column).name{1});
        U_new(pos_index).ons  =  job.SPM_orig.Sess.U(column).ons(other_trials) ;
        U_new(pos_index).dur  = job.SPM_orig.Sess.U(column).dur(other_trials);
        U_new(pos_index).P = job.SPM_orig.Sess.U(column).P;
        U_new(pos_index).dt = job.SPM_orig.Sess.U(column).dt;
        
        pos_index = pos_index + 1;
        
        for oc = 1:length(other_columns)
            other_column = other_columns(oc);
            
            U_new(pos_index).name{1} = sprintf('%s_%d', SPM.Sess.U(other_column).name{1} , i);
            U_new(pos_index).ons  =  job.SPM_orig.Sess.U(other_column).ons;
            U_new(pos_index).dur  = job.SPM_orig.Sess.U(other_column).dur;
            U_new(pos_index).P = job.SPM_orig.Sess.U(other_column).P;
            U_new(pos_index).dt = job.SPM_orig.Sess.U(other_column).dt;
            
            pos_index = pos_index+ 1;
        end
        
        %% store the new U structure
        SPM.Sess.U = U_new;
        
        BSM_finalize_model(job, SPM);
    end %% end loop for trial i
    
end %% end loop for active columns
end


    
function BSM_type3(job)
%% Fast approach

active_columns = job.active_columns;

for c = 1:length(active_columns) %% split out all columns but baseline
    column = active_columns(c);
    other_columns = setdiff(1:length(active_columns), column);
    
    display(sprintf('*) Exploding column %d...', column));
    
    target_folder_c = [job.target_folder filesep sprintf('col_%04d', column)];
    mkdir(target_folder_c);
    
    pos_index =1;
	
    SPM = job.SPM_orig;
    SPM.swd = target_folder_c;
    U_new = struct('name',[],'ons', [], 'dur',[], 'P',[], 'dt', [], 'u', [],'pst',[]);
    
    for i = 1:length(job.U_orig(column).ons) %% walk over all onsets within 'column'
                   
        %% add event i from 'column'
        U_new(pos_index).name{i} = sprintf('%s_%d', SPM.Sess.U(column).name{1} , i);
        U_new(pos_index).ons  =  job.SPM_orig.Sess.U(column).ons(i) ;
        U_new(pos_index).dur  = job.SPM_orig.Sess.U(column).dur(i);
        U_new(pos_index).P = job.SPM_orig.Sess.U(column).P;
        U_new(pos_index).dt = job.SPM_orig.Sess.U(column).dt;
        
        %% add all other events but i from 'column' to next column
        pos_index = pos_index +1;
    end
       
    for oc = 1:length(other_columns)
    	other_column = other_columns(oc);
            
        U_new(pos_index).name{1} = sprintf('%s_%d', SPM.Sess.U(other_column).name{1} , i);
        U_new(pos_index).ons  =  job.SPM_orig.Sess.U(other_column).ons;
        U_new(pos_index).dur  = job.SPM_orig.Sess.U(other_column).dur;
        U_new(pos_index).P = job.SPM_orig.Sess.U(other_column).P;
        U_new(pos_index).dt = job.SPM_orig.Sess.U(other_column).dt;
            
        pos_index = pos_index+ 1;
    end
        
    %% store the new U structure
    SPM.Sess.U = U_new;
        
    BSM_finalize_model(job, SPM);
        
    %% end loop for active columns
end
end

function BSM_type4(job)
active_columns = job.active_columns;

%% Adapted Fast approach
for c = 1:length(active_columns) %% split out all columns but baseline
    column = active_columns(c);
    other_columns = setdiff(1:length(active_columns), column);
    
    display(sprintf('*) Exploding column %d...', column));
    
    target_folder_c = [job.target_folder filesep sprintf('col_%04d', column)];
    mkdir(target_folder_c);
    
    pos_index =1;
	
    SPM = job.SPM_orig;
    SPM.swd = target_folder_c;
    U_new = struct('name',[],'ons', [], 'dur',[], 'P',[], 'dt', [], 'u', [],'pst',[]);
    
    for i = 1:length(job.U_orig(column).ons) %% walk over all onsets within 'column'
                   
        %% add event i from 'column'
        U_new(pos_index).name{i} = sprintf('%s_%d', SPM.Sess.U(column).name{1} , i);
        U_new(pos_index).ons  =  job.SPM_orig.Sess.U(column).ons(i) ;
        U_new(pos_index).dur  = job.SPM_orig.Sess.U(column).dur(i);
        U_new(pos_index).P = job.SPM_orig.Sess.U(column).P;
        U_new(pos_index).dt = job.SPM_orig.Sess.U(column).dt;
        
        %% add all other events but i from 'column' to next column
        pos_index = pos_index +1;
    end
    
    for oc = 1:length(other_columns)
        other_column = other_columns(oc);
        
        U_new(pos_index).name{1} = 'other conds';
        U_new(pos_index).ons  =  [U_new(pos_index).ons job.SPM_orig.Sess.U(other_column).ons];
        
        if length(job.SPM_orig.Sess.U(other_column).dur) == 1
            U_new(pos_index).dur  = [U_new(pos_index).dur repmat(job.SPM_orig.Sess.U(other_column).dur, 1, length(job.SPM_orig.Sess.U(other_column).ons))];
        else
            U_new(pos_index).dur  = [U_new(pos_index).dur M_orig.Sess.U(other_column).dur];
        end
        U_new(pos_index).P = job.SPM_orig.Sess.U(other_column).P;
        U_new(pos_index).dt = job.SPM_orig.Sess.U(other_column).dt;
    end
    
    
    %% store the new U structure        
    SPM.Sess.U = U_new;
    
    BSM_finalize_model(job, SPM);
        
    %% end loop for active columns
end
end

function BSM_finalize_model(job, SPM)

U = spm_get_ons(SPM,1);
SPM.Sess(1).U = U;
[X,Xn,Fc] = spm_Volterra(SPM.Sess.U,SPM.xBF.bf,SPM.xBF.Volterra);

k = SPM.nscan(1);
fMRI_T     = SPM.xBF.T;
fMRI_T0    = SPM.xBF.T0;

X = X((0:(k - 1))*fMRI_T + fMRI_T0 + 32,:);
% and orthogonalise (within trial type)
for i = 1:length(Fc)
    X(:,Fc(i).i) = spm_orth(X(:,Fc(i).i));
end

Xname={};
for i = 1:size(X,2)
    Xname{end + 1} = [sprintf('Sn(%i) ',1) SPM.Sess(1).U(i).name{1}];
end

%% reinstate motion parameters
if (job.keepr)
    C = job.SPM_orig.Sess(1).C.C;
    Cname = job.SPM_orig.Sess(1).C.name;
    
    X     = [X spm_detrend(C)];
    Xname    = {Xname{:} Cname{:}};
end    


%% reinstate baseline column
X(:,end+1) = 1;
SPM.xX.X = X;

SPM.xX.iB = size(X,2);
SPM.xX.iC = 1:size(X,2)-1;
Xname{end+1} = [sprintf('Sn(%i) ',1) 'constant'];
SPM.xX.name = Xname;

try
    SPM.xX = rmfield(SPM.xX, 'W');
catch
end
save([SPM.swd filesep 'SPM.mat'], 'SPM');

%% estimate model:
matlabbatch{1}.spm.stats.fmri_est.spmmat = {[SPM.swd filesep 'SPM.mat']};
matlabbatch{1}.spm.stats.fmri_est.method.Classical = 1;
spm_jobman('run', matlabbatch);
end

%% Concatenate betas per condition in a Beta series 4D file
function BSM_concat_betas(job)
active_columns = job.active_columns;
target_folder = job.target_folder;
if (job.type < 3)
    
    for c = 1:length(active_columns) %% split out all columns but baseline
        
        column = active_columns(c);        
        target_folder_c = [target_folder filesep sprintf('col_%04d', column)];
        beta_files = spm_select('FPListRec', target_folder_c, 'beta_0001.nii');
    
        matlabbatch{1}.spm.util.cat.vols = cellstr(beta_files);
        matlabbatch{1}.spm.util.cat.name = [target_folder_c filesep 'betas.nii'];
        matlabbatch{1}.spm.util.cat.dtype = 4;
        spm_jobman('run', matlabbatch);
        
        if (job.keepfl == 0)
            
                to_delete = spm_select('FPListRec', target_folder_c, '^*.*');
                items = ~regexp(to_delete, 'betas.nii');
                to_delete = to_delete(items);
                spm_unlink(to_delete);  %% throw away all files
        end
    end
    
else
    
    for c = 1:length(active_columns) %% split out all columns but baseline
        column = active_columns(c);
        beta_for_cond_c = 1:length(job.U_orig(column).ons);
        target_folder_c = [target_folder filesep sprintf('col_%04d', column)];
        beta_files = {}; %% empty beta_files
        for i =1:length(job.U_orig(column).ons) 
            beta_files{i} = spm_select('FPList', target_folder_c, sprintf('^beta_%04d.nii', i));
        end
        
        matlabbatch{1}.spm.util.cat.vols = beta_files';
        matlabbatch{1}.spm.util.cat.name = [target_folder_c filesep 'betas.nii'];
        matlabbatch{1}.spm.util.cat.dtype = 4;
        spm_jobman('run', matlabbatch);
        
        if (job.keepfl == 0)
            display('*) Cleaning up ...');
            to_delete = cellstr(spm_select('FPListRec', target_folder_c, '^*.*'));
            items_to_retain = ~cellfun(@isempty,regexp(to_delete, 'betas'));
            to_delete(find(items_to_retain)) = [];
            
            for f = 1:length(to_delete)
                spm_unlink(to_delete{f});  %% throw away all files
            end
        end
        
    end    
end
end

%% Step 3: Use Betas series files for GLM on a ROI timecourse (average Beta series extracted for a ROI)
%% Organized now as follows  : / BSM / <ROI_filename> / col_000n
function[] = BSM_Regress_Regions(job)
target_folder = job.target_folder;
active_columns = job.active_columns;

for r = 1:length(job.rois)
    ROI(r).volume = spm_read_vols(spm_vol(job.rois{r}));
    ROI(r).ind = find(ROI(r).volume);          
    [ROI(r).x ROI(r).y ROI(r).z] = ind2sub(size(ROI(r).volume), ROI(r).ind);
end
for r = 1:length(job.rois)
    
    if (isempty(ROI(r).ind))
        ROI(r) = [];
        warning([sprintf('  %d ) Problem found for ROI : ''%s'', no active voxels.\n',r, job.rois{r}) ...
                         sprintf('This ROI will be skipped\n')]);
    else
        display(sprintf('  %d ) ROI ''%s'' contains %d voxels', r, job.rois{r}, length(ROI(r).ind)));
    end
    
end

for r = 1:length(job.rois)
    for c = 1:length(job.active_columns)
        column = active_columns(c);        
        target_folder_c = [target_folder filesep sprintf('col_%04d', column)];
        
        betas = cellstr(spm_select('ExtFPList', target_folder_c, 'betas.nii', 1:999));
        B = zeros(1, length(betas));
        
        for b = 1:length(betas)
            H = spm_vol(betas{b});
            B(b) = mean(spm_sample_vol(H, ROI(r).x, ROI(r).y, ROI(r).z,0));
        end
        
        [base roiname ext]= fileparts(job.rois{r});
        target_folder_r = [target_folder filesep roiname filesep sprintf('col_%04d', column)];
        mkdir(target_folder_r);
        
        matlabbatch{1}.spm.stats.factorial_design.dir = { target_folder_r};
        matlabbatch{1}.spm.stats.factorial_design.des.t1.scans = betas;
        matlabbatch{1}.spm.stats.factorial_design.cov.c = B;
        matlabbatch{1}.spm.stats.factorial_design.cov.cname = roiname;
        matlabbatch{1}.spm.stats.factorial_design.cov.iCFI = 1;
        matlabbatch{1}.spm.stats.factorial_design.cov.iCC = 1;
        matlabbatch{1}.spm.stats.factorial_design.multi_cov = struct('files', {}, 'iCFI', {}, 'iCC', {});
        matlabbatch{1}.spm.stats.factorial_design.masking.tm.tm_none = 1;
        matlabbatch{1}.spm.stats.factorial_design.masking.im = 1;
        matlabbatch{1}.spm.stats.factorial_design.masking.em = {''};
        matlabbatch{1}.spm.stats.factorial_design.globalc.g_omit = 1;
        matlabbatch{1}.spm.stats.factorial_design.globalm.gmsca.gmsca_no = 1;
        matlabbatch{1}.spm.stats.factorial_design.globalm.glonorm = 1;
        
        matlabbatch{2}.spm.stats.fmri_est.spmmat(1) = cfg_dep('Factorial design 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;
        
        matlabbatch{3}.spm.stats.con.spmmat(1) = cfg_dep('Model estimation: SPM.mat File', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','spmmat'));
        matlabbatch{3}.spm.stats.con.consess{1}.tcon.name = 'r';
        matlabbatch{3}.spm.stats.con.consess{1}.tcon.weights = [0 1];
        matlabbatch{3}.spm.stats.con.consess{1}.tcon.sessrep = 'none';
        matlabbatch{3}.spm.stats.con.consess{2}.tcon.name = '-r';
        matlabbatch{3}.spm.stats.con.consess{2}.tcon.weights = [0 -1];
        matlabbatch{3}.spm.stats.con.consess{2}.tcon.sessrep = 'none';
        matlabbatch{3}.spm.stats.con.delete = 0;
        
        spm_jobman('run', matlabbatch);
    end
end

end

%% Step 4: Contrast Second level Betas for each ROI
function[] = BSM_Contrast_manager(job)
target_folder = job.target_folder;
for r = 1:length(job.rois)
    [base roiname ext]= fileparts(job.rois{r});
    target_folder_final = [target_folder filesep roiname filesep 'results'];
    mkdir(target_folder_final);
    
    for c1 = 1:length(job.active_columns)
        for c2 = setdiff(1:length(job.active_columns),c1)
            
            betafile_c1 = spm_select('FPList', [target_folder filesep roiname filesep sprintf('col_%04d', c1)], 'beta_0002.nii');
            betafile_c2 = spm_select('FPList', [target_folder filesep roiname filesep sprintf('col_%04d', c2)], 'beta_0002.nii');
            
            col1 = job.active_columns(c1);
            col2 = job.active_columns(c2);
            
            matlabbatch{1}.spm.util.imcalc.input = { betafile_c1; betafile_c2 };
            matlabbatch{1}.spm.util.imcalc.output = sprintf('con_%04d_minus_%04d.nii', col1, col2);
            matlabbatch{1}.spm.util.imcalc.outdir = { target_folder_final};
            matlabbatch{1}.spm.util.imcalc.expression = 'i1 - i2';
            matlabbatch{1}.spm.util.imcalc.var = struct('name', {}, 'value', {});
            matlabbatch{1}.spm.util.imcalc.options.dmtx = 0;
            matlabbatch{1}.spm.util.imcalc.options.mask = 0;
            matlabbatch{1}.spm.util.imcalc.options.interp = 1;
            matlabbatch{1}.spm.util.imcalc.options.dtype = 4;
            
            matlabbatch{2}.spm.util.imcalc.input = { betafile_c1; betafile_c2 };
            matlabbatch{2}.spm.util.imcalc.output = sprintf('con_%04d_minus_%04d.nii', col2, col1);
            matlabbatch{2}.spm.util.imcalc.outdir = { target_folder_final};
            matlabbatch{2}.spm.util.imcalc.expression = 'i2 - i1';
            matlabbatch{2}.spm.util.imcalc.var = struct('name', {}, 'value', {});
            matlabbatch{2}.spm.util.imcalc.options.dmtx = 0;
            matlabbatch{2}.spm.util.imcalc.options.mask = 0;
            matlabbatch{2}.spm.util.imcalc.options.interp = 1;
            matlabbatch{2}.spm.util.imcalc.options.dtype = 4;
            
            spm_jobman('run', matlabbatch);
        end
    end
     display(sprintf(' *) Results for roi %s are stored in : %s', roiname, job.target_folder_final));
   end
   display(' *) These are difference in beta images and can be viewed via mricron or xjview');

end