test_rig

 runs tests on MarsBaR using specified designs
 FORMAT res = test(design_paths, params)
 
 Inputs
 design_paths     - path(s) to SPM design files
 params           - structure giving params to pass to estimate method,
                    see help for do_estimate methods for details
                    Default is
                    params = struct('redo_covar', 0, ...
                                'redo_whitening', 0);
 
 Outputs
 res              - 1 if all tests passed, 0 otherwise
 
 The function depends on the SPM design having estimated contrasts to
 play with.  It uses these to:
 Get the maximum voxel in the first F and first T contrast
 Records the T/F statistic value
 Makes an ROI out of this voxel
 Estimates in MarsBaR
 Checks the statistic value is that same.
 
 Along the way, it uses much of the MarsBaR machinery
 
 $Id$ 

0001 function res = test_rig(design_paths, params)
0002 % runs tests on MarsBaR using specified designs
0003 % FORMAT res = test(design_paths, params)
0004 %
0005 % Inputs
0006 % design_paths     - path(s) to SPM design files
0007 % params           - structure giving params to pass to estimate method,
0008 %                    see help for do_estimate methods for details
0009 %                    Default is
0010 %                    params = struct('redo_covar', 0, ...
0011 %                                'redo_whitening', 0);
0012 %
0013 % Outputs
0014 % res              - 1 if all tests passed, 0 otherwise
0015 %
0016 % The function depends on the SPM design having estimated contrasts to
0017 % play with.  It uses these to:
0018 % Get the maximum voxel in the first F and first T contrast
0019 % Records the T/F statistic value
0020 % Makes an ROI out of this voxel
0021 % Estimates in MarsBaR
0022 % Checks the statistic value is that same.
0023 %
0024 % Along the way, it uses much of the MarsBaR machinery
0025 %
0026 % $Id$
0027   
0028 if nargin < 1
0029   design_paths = spm_get([0 Inf], 'SPM*.mat', 'Select SPM designs');
0030 end
0031 if nargin < 2
0032   params = struct('redo_covar', 0, ...
0033           'redo_whitening', 0);
0034 end
0035 
0036 n_designs = size(design_paths, 1);
0037 res = zeros(n_designs, 1);
0038 for d = 1:n_designs
0039   d_path = deblank(design_paths(d,:));
0040   res(d) = sf_test_design(d_path, params);
0041 end
0042 return
0043 
0044 function res = sf_test_design(d_path, params)
0045 % tests one design
0046   
0047 % Check for SPM estimated design, with estimated contrasts
0048 D = mardo(d_path);
0049 if ~is_spm_estimated(D)
0050   error('Need an SPM estimated design');
0051 end
0052 if ~has_contrasts(D)
0053   error(['Design ' d_path ' does not contain contrasts']);
0054 end
0055 if ~has_images(D)
0056   error(['Design ' d_path ' does not contain images']);
0057 end
0058 
0059 % try to get one F and one T contrast
0060 Swd = fileparts(d_path);
0061 xCon = get_contrasts(D);
0062 stats = [xCon(:).STAT];
0063 Ic    = []; fnames = {};
0064 for t = 'TF'
0065   for c = fliplr(find(stats == t))
0066     F = xCon(c).Vspm;
0067     if ~isempty(F)
0068       % SPM99 = filename, SPM2 = vol_struct
0069       if isstruct(F), F = F.fname; end
0070       % SPM5 has full paths for the contrast images
0071       con_pth = fileparts(F);
0072       if isempty(con_pth)
0073           F = fullfile(Swd, F);
0074       end
0075       if exist(F, 'file'), Ic = [Ic c]; fnames{end+1} = F; break, end
0076     end
0077   end
0078 end
0079 if isempty(Ic)
0080   error(['Could not find any contrast images for ' d_path]);
0081 end
0082 
0083 % find maximum voxel coordinate for contrasts and test
0084 res = 1;
0085 for c = 1:length(Ic)
0086   V = spm_vol(fnames{c});
0087   img = spm_read_vols(V);
0088   [mx(c) i] = max(img(:));
0089   xyz(:, c) = mars_utils('e2xyz', i, V.dim(1:3));
0090   mx_roi(c) = maroi_pointlist(struct('XYZ', xyz(:, c), ...
0091                 'mat', V.mat), 'vox');
0092   Y = get_marsy(mx_roi(c), D, 'mean');
0093   E = estimate(D, Y, params);
0094   [E n_Ic] = add_contrasts(E, D, Ic(c));
0095   marsS = compute_contrasts(E, n_Ic);
0096   fprintf('SPM statistic %7.4f; MarsBaR statistic %7.4f\n',...
0097             mx(c), marsS.stat(1));
0098   st_spm = mx(c);
0099   st_mars = marsS.stat(1);
0100   bad_test = abs(st_mars - st_spm) > 1e-5;
0101   if bad_test % Statistics are different - SPM8 fudge?
0102     spmV = lower(mars_utils('spm_version'));
0103     if any(strcmp(spmV, {'spm8', 'spm12b', 'spm12'}))
0104       xCon = get_contrasts(E);
0105       this_con = xCon(n_Ic);
0106       if this_con.STAT == 'T' & (st_mars > st_spm)
0107         % Probably the fudge factor
0108         fudge = (st_mars / st_spm - 1) / exp(-8);
0109         fprintf('Fudge value was %f\n', fudge);
0110         % fudge should now be max SE across all voxels in the
0111         % estimation block (e.g. slice) / SE for this voxel
0112         if fudge < 10
0113           bad_test = 0;
0114         end
0115       end
0116     end
0117   end
0118   if bad_test
0119     disp('MarsBaR gives a different result for contrast');
0120     res = 0;
0121   end
0122 end
0123 return