% Batch for within-subject ANOVAs in SPM5       Rik Henson  Jan 06

spm('defaults', 'FMRI')
global defaults
global UFp
UFp=0.001;

bwd = '/myhome';
cd(bwd)

subnam = [060001 060002 060003 060004 060005 060013 060014 060015 060016 060017 060018 060019 060023 060024 060025 060035 060036];
nsub = length(subnam);

ses = 'Ses1'; cons = [1:12]; ana = 'Mono12'; 

ncon = length(cons);
fnm = 'con';
%fnm = 'beta';
nsc = 1;	

nscan = nsub*ncon;

nwd = sprintf('%s/Group/%s/%s',bwd,ses,ana)
eval(sprintf('!mkdir %s',nwd));
cd(nwd);


% get image files names
P={};
for c=1:ncon
	con = cons(c);
	for s=1:nsub
	     	sub = subnam(s);
      	     	P{(c-1)*nsub+s} = sprintf('%s/CBU%06d/%s/Stats/%s_%04d.img',bwd,subnam(s),ses,fnm,con);
	end 
end

%-Assemble SPM structure
%=======================================================================
clear SPM

SPM.nscan = nscan;
SPM.xY.P = P;

for i=1:SPM.nscan
 SPM.xY.VY(i) = spm_vol(SPM.xY.P{i});
end

cname = cell(ncon,1);
for c=1:ncon
	cname{c} = sprintf('con %d',c);
end
for c=1:nsub
    cname{c+ncon} = sprintf('sub %d',c);
end

os = ones(nsub,1); oc = ones(ncon,1);
SPM.xX = struct(...
        'X',[kron(eye(ncon),os) kron(oc,eye(nsub))] ,...
        'iH',[1:(ncon+nsub)],'iC',zeros(1,0),'iB',zeros(1,0),'iG',zeros(1,0),...
        'name',{cname},'I',[ones(nscan,1) kron([1:ncon]',os) kron(oc,[1:nsub]') ones(nscan,1)],...
        'sF',{{'repl'  'cond'  'subj'  ''}});

SPM.xC = [];	

SPM.xGX = struct(...
        'iGXcalc',1,    'sGXcalc','omit',                               'rg',[],...
        'iGMsca',9,     'sGMsca','<no grand Mean scaling>',...
        'GM',0,         'gSF',ones(nscan,1),...
        'iGC',  12,     'sGC',  '(redundant: not doing AnCova)',        'gc',[],...
        'iGloNorm',9,   'sGloNorm','<no global normalisation>');

x=zeros(ncon);
s=eye(nsub);
nv=0;
% first do the leading diagonal elements
for i=1:ncon
        nv=nv+1;
        v=x;
        v(i,i)=1;
        vi{nv}=sparse(kron(v,s));
end
% then do the off-diagonals
for i=1:(ncon-1)
        for j=(i+1):ncon
                nv=nv+1;
                v=x;
                v(j,i)=1; v(i,j)=1;
                vi{nv}=sparse(kron(v,s));
       end
end

if ncon==1 | nsc==0
 SPM.xVi = struct('iid',1, 'V',speye(nsub*ncon) );       
else
 SPM.xVi = struct('iid',0, 'I',SPM.xX.I, 'Vi',{vi} );       
end

Mdes    = struct(...	
        'Analysis_threshold',   {'None (-Inf)'},...
        'Implicit_masking',     {'Yes: NaNs treated as missing'},...
        'Explicit_masking',     {'No'});

SPM.xM  = struct(...
        'T',-Inf,'TH',ones(nscan,1)*-Inf,...
        'I',1,'VM',[],'xs',Mdes);

Pdes    = {{sprintf('%d condition, +0 covariate, +%d block, +0 nuisance',ncon,nsub); sprintf('%d total, having %d degrees of freedom',ncon+nsub,ncon+nsub-1); sprintf('leaving %d degrees of freedom from %d images',nscan-ncon-nsub+1,nscan)}};

Pdes{1}

SPM.xsDes = struct(...
        'Design',               {'1-way ANOVA (within-subjects)'},...
        'Global_calculation',   {'omit'},...
        'Grand_mean_scaling',   {'<no grand Mean scaling>'},...
        'Global_normalisation', {'<no global normalisation>'},...
        'Parameters',           Pdes);

save SPM SPM


% Estimate parameters
%===========================================================================
SPM = spm_spm(SPM);


SPM = rmfield(SPM,'xCon');

c = eye(ncon);
c = [c zeros(ncon,nsub)];
cname       = 'Unwhitened EOI';
SPM.xCon(1) = spm_FcUtil('Set',cname,'F','c',c',SPM.xX.xKXs);


% F-contrasts 
%-------------

hemi = [1 1; 1 -1]; hemin = {'L+R';'L-R'};
stim = [1 1; 1 -1]; stimn = {'F+H';'F-H'};
reps = [1 1 1; 2 -1 -1]; repsn = {'1+2';'1-2'};

cc=1
for nh=1:size(hemi,1)
 for ns=1:size(stim,1)
  for nr=1:size(reps,1)
   c = kron(hemi(nh,:),stim(ns,:));
   c = kron(c,reps(nr,:));
   c = norm_cons(c);
   c = [c zeros(1,nsub)];
   cc=cc+1;
   cname = sprintf('%s X %s X %s',hemin{nh},stimn{ns},repsn{nr});
   SPM.xCon(cc) = spm_FcUtil('Set',cname,'F','c',c',SPM.xX.xKXs);
  end
 end
end

 
% T-contrasts 
%-------------

hemi = [1 1; 1 -1; -1 1]; hemin = {'L+R';'L-R';'R-L'};
stim = [1 1; 1 -1; -1 1; 1 0; 0 1]; stimn = {'F+H';'F-H';'H-F';'F';'H'};
reps = [1 1 1; 1 -1 -1; -1 1 1; 1 -1 0; -1 1 0];
repsn = {'All';'1-2';'2-1';'1-S';'S-1'};

%cc=1
for nh=1:size(hemi,1)
 for ns=1:size(stim,1)
  for nr=1:size(reps,1)
   c = kron(hemi(nh,:),stim(ns,:));
   c = kron(c,reps(nr,:));
   c = norm_cons(c);
   c = [c zeros(1,nsub)];
   cc=cc+1;
   cname = sprintf('%s X %s X %s',hemin{nh},stimn{ns},repsn{nr});
   SPM.xCon(cc) = spm_FcUtil('Set',cname,'T','c',c',SPM.xX.xKXs);
  end
 end
end

spm_contrasts(SPM);