%-----------------------------------------------------------------------
%Design specific parameters
%-----------------------------------------------------------------------

clear
spm_defaults;

subs = {'030093' '050637' '050636' '040317' '050714' '040630' '030721' '041031' '050715' '050716' '050717' '050718' '050719' '040700' '040858'};
nsubs = length(subs);
nsess = 2;
ntot = nsubs*nsess;
cnames{1} = {'S-nc' 'S-cs' 'S-cc' 'S-test' 'S-first'};
cnames{2} = {'C-nc' 'C-cs' 'C-cc' 'C-test' 'C-first'};
ncond = length(cnames{1});
TR = 1.1005;
dataroot = '/imaging/russell/expt456';
batchroot = '/imaging/russell/test5';
statsdir = fullfile(batchroot,'stats');
border = [2 1 2 1 2 1 2 1 2 1 1 2 1 2 2];
hpf = 60;
incmoves = 1;
modeldur = 0;
imgfilt = '^snuraW.*\.img$';
movefilt = '^rp_.*\.txt';
evdir = fullfile(dataroot,'stats/e6/change/m2/events');
mnames = {'x_trans' 'y_trans' 'z_trans' 'x_rot' 'y_rot' 'z_rot'};
debug=0;

cname = 'An example contrast';
cons{1} = [repmat([-1 0 1 0 0 0 0 0 0 0 0],1,2) 0 0];



%-----------------------------------------------------------------------
%Design setup
%-----------------------------------------------------------------------

% basis functions and timing parameters
%---------------------------------------------------------------------------
% OPTIONS:'hrf'
%         'hrf (with time derivative)'
%         'hrf (with time and dispersion derivatives)'
%         'Fourier set'
%         'Fourier set (Hanning)'
%         'Gamma functions'
%         'Finite Impulse Response'
%---------------------------------------------------------------------------

%MD No slice timing, therefore using temporal-derivates to shift HRF forward and back in time

xBF.name       = 'hrf';
xBF.length     = 32.2;              % length in seconds
xBF.order      = 1;                 % order of basis set
xBF.T          = 16;                % number of time bins per scan
xBF.T0         = 1;                 % first time bin (see slice timing) - middle of TA
xBF.UNITS      = 'scans';           % OPTIONS: 'scans'|'secs' for onsets
xBF.Volterra   = 1;                 % OPTIONS: 1|2 = order of convolution




tc = 0;
allfiles = '';
for sub = 1:2%nsubs
    clear SPM
    disp(subs{sub})
    SPM.xY.RT = TR;
    SPM.xGX.iGXcalc = 'None';
    SPM.xVi.form = 'AR(1)';
    SPM.xBF = xBF;
    csub = subs{sub};
    subdata = fullfile(dataroot, csub);
    anadir = fullfile(statsdir,csub);
    if exist(anadir)~=7; mkdir(statsdir,csub);end
    cd(anadir);
    
    tc = 0;
    allfiles='';
    condtog = border(sub);
    
    for sess = 1:nsess
        tc = tc+1;
        sessdata = fullfile(subdata,['E6-Sess' num2str(sess)]);
        ons = spm_load(fullfile(evdir,['ev_' csub '_' num2str(sess) '.txt']));
        files = spm_select('List', sessdata, imgfilt);
        if incmoves==1
            mfname = spm_select ('List', sessdata, movefilt);
            moves = load(fullfile(sessdata,mfname));
        end
        SPM.nscan(tc) = size(files,1);
        for f =1:size(files,1)
            ffiles(f,:) = fullfile(sessdata,files(f,:));
        end
        allfiles = strvcat(allfiles,ffiles);
        
        for c = 1:ncond
            cidx = find(ons(:,2)==c);
            
            if modeldur==1
                durs = ons(cidx,3);
            else
                durs = zeros(length(cidx),1);
            end
            
            SPM.Sess(tc).U(c) = struct(...
                'ons',ons(cidx,1),...
                'dur',durs,...
                'name',{cnames{condtog}(c)},...
                'P',struct('name','none'));
        end
        
        SPM.xX.K(tc).HParam = hpf;
        if condtog == 1
            condtog = 2;
        elseif condtog == 2
            condtog = 1;
        end
        
        if incmoves==1
            SPM.Sess(tc).C.C    = moves;     % [n x c double] covariates
            SPM.Sess(tc).C.name = mnames; % [1 x c cell]   names
        else
            SPM.Sess(tc).C.C = [];
            SPM.Sess(tc).C.name = {};
        end
    end
    
    if debug==0
        SPM.xY.P = allfiles;
        SPMdes = spm_fmri_spm_ui(SPM);
        spm_unlink(fullfile('.', 'mask.img')); % avoid overwrite dialog
        SPMest = spm_spm(SPMdes);
        for i = 1:size(cons,2)
            if length(SPMest.xCon)==0
                SPMest.xCon = spm_FcUtil('Set',cname,'T','c',cons{i}',SPMest.xX.xKXs);
            else
                SPMest.xCon(end+1) = spm_FcUtil('Set',cname,'T','c',cons{i}',SPMest.xX.xKXs);
            end
        end
        spm_contrasts(SPMest);
    end
end