diff --git a/README.md b/README.md index fd5368d..ee1fdde 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,199 @@ -### SnPM: Statistical nonParametric Mapping +## SnPM: Statistical nonParametric Mapping -The Statistical non-Parametric Mapping (SnPM) toolbox provides an extensible framework for voxel level non-parametric permutation/randomisation tests of functional Neuroimaging experiments with independent observations. +The Statistical non-Parametric Mapping (SnPM) toolbox provides an extensible framework for voxel level non-parametric permutation/randomisation tests of functional Neuroimaging experiments with independent observations. The SnPM toolbox provides an alternative to the Statistics section of [SPM](http://www.fil.ion.ucl.ac.uk/spm/). SnPM uses the General Linear Model to construct pseudo t-statistic images, which are then assessed for significance using a standard non-parametric multiple comparisons procedure based on randomisation/permutation testing. It is most suitable for single subject PET/SPECT analyses, or designs with low degrees of freedom available for variance estimation. In these situations the freedom to use weighted locally pooled variance estimates, or variance smoothing, makes the non-parametric approach considerably more powerful than conventional parametric approaches, as are implemented in SPM. Further, the non-parametric approach is always valid, given only minimal assumptions. + +**More information at: www.warwick.ac.uk/snpm** + + - [Getting started](http://www.warwick.ac.uk/snpm/man) + - [fMRI example](http://www.warwick.ac.uk/snpm/man/exnew) + - [PET example](http://www.warwick.ac.uk/snpm/man/ex) + - [Download](http://www.warwick.ac.uk/snpm/snpmreg) + + +### Testing + +#### Initial set up +The first time you run the tests, you will first have to create a set of ground truth data and a file named `snpm_test_config.m` containing the reference SnPM version you used to compute the ground truth and the path to the ground truth folder. For example: +``` +global testDataDir; +testDataDir = '~/snpm_test_data'; +global SnPMrefVersion; +SnPMrefVersion = 'SnPM8'; +``` + +#### Run the test suite +Then, the tests can be started (from the test data folder) with: +``` +import matlab.unittest.TestSuite; +suite = TestSuite.fromFolder(fullfile(spm_str_manip(which('snpm'), 'h'), 'test')); +result = run(suite); +``` + +#### Run a single test +``` +run(test_oneSample, 'test_onesample_1') +``` + +### SnPM13: Bugs & Fixes + +This section describes the bugs that have been reported, along with the appropriate fixes. +Updated versions of appropriate SnPM functions are available at [snpm13_updates](http://warwick.ac.uk/snpm/distribution/snpm13_updates) + +#### SnPM 13.1.04 +* snpm +SnPM version 13.1.04 + +* snpm_bch_ui, snpm_cp, snpm_ui +Fix: implicit masking for non-float datatypes + +* snpm_pi_ANOVAbtwnS, snpm_pi_ANOVAwithinS, snpm_pi_Corr, snpm_pi_Corr1S, +snpm_pi_OneSampT, snpm_pi_TwoSampPairT, snpm_pi_TwoSampT, snpm_pi_TwoSampTss +Update syntax to set the seed of the random number generator (rng) + +* snpm_defaults +New SnPM global parameter shuffle_seed. If true, the rng uses a 'shuffle' seed providing different results for each run. If false Matlab's 'default' seed or any seed defined by the user will be used . + +* snpm_pi_TwoSampT +Allows MC permutation (with reptitions) for large enough nPerms + +* snpm_pp, snpm_ui +Do not display figures in command-line mode + +* snpm, snpm_check_nperm, snpm_combo, snpm_cp, snpm_init, snpm_pi_ANOVAbtwnS, +snpm_pi_ANOVAwithinS, snpm_pi_Corr, snpm_pi_Corr1S, snpm_pi_OneSampT, +snpm_pi_PairT, snpm_pi_PairTrand, snpm_pi_TwoSampPairT, snpm_pi_TwoSampT, +snpm_pi_TwoSampTss, snpm_pp, snpm_t2z, snpm_uc_FDR, snpm_ui, spm_append_96 +Add identifiers to errors and warnings + +* snpm_cp +Fix: Do not assume statistic is T when computing cluster forming threshold from P-value + +* generic_test_snpm, skeleton_class, test_ANOVAbetween, +test_multisub_withinsubanova, test_multisubpaired2cond, +test_multisubsimpleregression, test_oneSample, test_onesub_regression, +test_onesub_twocondrepl, test_twoSample, test_twosample_twocond +Test in command-line mode with no shuffling for the random seed, using new rng + syntax. + +#### Updates from SnPM 13.1.03 +* snpm +SnPM version 13.1.03 + +* snpm_pi_PairT +Now allows more than 52 subjects, as previously that would generate a "Maximum variable size allowed by the program" error message. + +#### Updates from SnPM 13.1.02 +* snpm +SnPM version 13.1.02 + +* snpm_pp, snpm_combo_pp +Update contrast display for compatibility with Matlab R2014b. + +* generic_test_snpm +Updates affecting tests only. + +#### Updates from SnPM 13.1.01 +* snpm +SnPM version 13.1.01 + +* snpm_pp +Check that 'Locs_vox' is integer according to a pre-defined tolerance level (as in snpm_cp). + +#### Updates from SnPM13.1.00 +* snpm_cp, snpm_pp, generic_test_snpm, test_oneSample +Cluster-forming threshold defined by a P-value are now extended to pseudo-T statistic (with variance smoothing). The Z-score corresponding to the selected P-value is used as an approximation of pseudo-T values. + +#### Updates from SnPM13.0.14 +* snpm_bch_ui_Corr, snpm_bch_ui_Corr1S, snpm_pi_Corr, snpm_pi_Corr1S, test_multisubsimpleregression, test_onesub_regression +Introduce nuisance covariates in simple-subject and multi-subject regressions. + +#### Updates from SnPM13.0.13 +* snpm_combo_pp +Fix: Set bNeg to 0 for F-tests. + +#### Updates from SnPM13.0.12 +* snpm_pi_OneSampT +Fix: Added new random mode for sampling of sign flips (Above 53 scans the binary-coding fails as Matlab's double significand runs out of precision). + +#### Updates from SnPM13.0.11 +* snpm_ui +Fix: When relative thresholding is set, then global calculation (user defined or mean) must be computed. + +* snpm_pi_TwoSampPairT, snpm_pi_TwoSampT +Fix: Unbalanced two-sample tests were wrongly aborted. + +* test_multisubpaired2cond, test_oneSample, test_twoSample, test_twosample_twocond +Updates affecting tests only. + +#### Updates from SnPM13.0.10 +* snpm_ui +Fix: ANCOVA normalisation (with an image presenting negative mean using spm_globals). + +#### Updates from SnPM13.09 +* snpm_cp +Increase tolerance in checking that 'Locs_vox' is integer. + +#### Updates from SnPM13.08 +* snpm_cp, spm_append_96 +Display warning message if SnPM_ST file becomes very large. Display more meaningful error message when Locs_vox is not interger. + +* snpm_pp +Display warning message if SnPM_ST can not be loaded. + +#### Updates from SnPM13.07 +* snpm +For backward-compatibility, open SPM batch window and display where to find SnPM tools when 'snpm' in called at the Matlab prompt. + +* snpm_pp +Fix: slow cluster inference with variance smoothing. + +* generic_test_snpm, test_oneSample, test_twoSample +Updates affecting tests only. + +#### Updates from SnPM13.06 +* snpm_ui +Apply user-defined scaling. + +* generic_test_snpm, test_oneSample, test_onesub_twocondrepl, test_twoSample +Updates affecting tests only. + +#### Updates from SnPM13.05 +* snpm_STcalc +Cluster labels from spm_max are not necessarily continuous. + +* snpm_pi_ANOVAbtwnS, snpm_pi_ANOVAwithinS, snpm_pi_Corr, snpm_pi_Corr1S, snpm_pi_OneSampT, snpm_pi_PairT, snpm_pi_TwoSampPairT, snpm_pi_TwoSampT, snpm_pi_TwoSampTss +Fix: if number of requested permutations is equal to the maximum number of possible permutations. + +* snpm_pi_TwoSampTss, snpm_pi_TwoSampPairT, snpm_pi_TwoSampT +Compute possible permutations (PiCond) in a more efficient way to avoid potential memory faults for large groups. + +#### Updates from SnPM13.04 +* snpm_ui +Apply absolute and relative masking. + +#### Updates from SnPM13.03 +* snpm_pi_TwoSampT +Compute possible permutations (PiCond) in a more efficient way to avoid potential memory faults for large groups. + +#### Updates from SnPM13.02 +* snpm_pp +Fix: if the extension is missing in the provided result file name, use .nii by default. +Use NaNs in thresholded map background instead of zeros. + +* [deleted spm_SnPM13] +SnPM is now only accessible through the Batch window (SPM -> Tools -> SnPM). + +* snpm_pi_ANOVAbtwnS, snpm_pi_ANOVAwithinS, snpm_pi_Corr1S, snpm_pi_OneSampT, snpm_pi_PairT, snpm_pi_TwoSampPairT, snpm_pi_TwoSampT +Minor update in the display. + +* test/ +Miror corrections affecting tests procedure only. + +#### Updates from SnPM13.01 +* snpm_cp +Synopsis: When half the permutations are computed (i.e. bhPerms is true) and StartPerm is 2 then the initial number of permutation could be 2 depending on sign of T0. + +* snpm_test_config, snpm_test_gound_truth and generic_test_snpm +Miror corrections affecting tests procedure only. diff --git a/config/snpm_bch_ui.m b/config/snpm_bch_ui.m index 6270aad..9d9d85e 100644 --- a/config/snpm_bch_ui.m +++ b/config/snpm_bch_ui.m @@ -283,7 +283,9 @@ '' 'For image data-types without a representation of NaN, zero is the mask value, and the user can choose whether zero voxels should be masked out or not.' '' - 'By default, an implicit mask is used. ' + 'By default, an implicit mask is used.' + '' + 'It is ill-advised to mix input images, but if data are a mix of integer (no NaN''s) and float/double datatypes, implicit masking will only occur with NaN''s.' '' }'; im.labels = { diff --git a/snpm.m b/snpm.m index 1a9477b..d321c7d 100644 --- a/snpm.m +++ b/snpm.m @@ -49,7 +49,7 @@ %-Parameters %----------------------------------------------------------------------- -SnPMver = 'SnPM13.1.03'; +SnPMver = 'SnPM13.1.04'; %-Format arguments %----------------------------------------------------------------------- @@ -106,7 +106,7 @@ else %======================================================================= -error('Unknown action string') +error('SnPM:UnknownAction', 'Unknown action string') %======================================================================= end diff --git a/snpm_check_nperm.m b/snpm_check_nperm.m index 1686167..e42e745 100644 --- a/snpm_check_nperm.m +++ b/snpm_check_nperm.m @@ -16,17 +16,21 @@ function snpm_check_perm(nPerm,TotPerm) if nPerm<100 & nPerm>=0.9*TotPerm str = sprintf(['Very few (%d) permutations used, nonparametric P-values are very coarse (but exact).\n'... 'Smallest possible P-value is %0.4f.'],nPerm,1/nPerm); + id = 'SnPM:VeryFewPermsCoarseExactPValues'; elseif nPerm<100 & nPerm<0.9*TotPerm str = sprintf(['Very few (%d) permutations used, nonparametric P-values are very coarse will vary \n'... 'substantially over repeated re-analyses.'... 'Smallest possible P-value is %0.4f.'],nPerm,1/nPerm); + id = 'SnPM:VeryFewPermsCoarseVaryingPValues'; elseif nPerm>10000 str = sprintf(['Many (%d) permutations used, analysis may take a very long time\n'... 'to complete. Consiser running fewer permutations.'],nPerm,1/nPerm); + id = 'SnPM:ManyPerms'; elseif nPerm>25000 str = sprintf(['%cInsane nummber of permutations requested (%d)!\n'... 'Are you sure!?'],7,nPerm); + id = 'SnPM:ManyManyPerms'; end if ~isempty(str) - warning(str) + warning(id, str) end diff --git a/snpm_combo.m b/snpm_combo.m index 71200f7..825c55d 100644 --- a/snpm_combo.m +++ b/snpm_combo.m @@ -126,7 +126,7 @@ else %======================================================================= -error('Unknown action string') +error('SnPM:UnknownAction', 'Unknown action string') %======================================================================= end diff --git a/snpm_combo_pp.m b/snpm_combo_pp.m index 2f90554..37078de 100644 --- a/snpm_combo_pp.m +++ b/snpm_combo_pp.m @@ -179,10 +179,10 @@ function snpm_combo_pp(CWD, job) str = 'Voxelwise corrected threshold = %g, which is smaller '; str = [str 'than minimum saved suprathreshold information (%g)']; str = [str '\nAll results significant voxelwise.']; - warning(sprintf(str,C_MaxT,ST_Ut)) + warning('SnPM:VoxelwiseCorrThreshSmaller', sprintf(str,C_MaxT,ST_Ut)) end if bSpatEx~=1 - error('Suprathreshold stats not collected! Cannot do cluster-combining!') + error('SnPM:SupraStatMissing', 'Suprathreshold stats not collected! Cannot do cluster-combining!') end %-Get primary threshold for STC analysis if requested @@ -202,12 +202,14 @@ function snpm_combo_pp(CWD, job) % upper tail p-values to specify primary threshold if alpha == 1 % Not filtering on significance if ~(primaryThresh>=ST_Ut) - error(['Using pseudo-statistics you can''t use (uncorrected)'... + error('SnPM:PseudoStatWithPValueThreshold', ... + ['Using pseudo-statistics you can''t use (uncorrected)'... 'upper tail p-values to specify primary threshold']); end else if ~(primaryThresh>=ST_Ut && primaryThresh=ST_Ut || (primaryThresh>0 && primaryThresh<=pU_ST_Ut_filt)) - error(['Primary threshold must be >=' num2str(ST_Ut) ... + error('SnPM:InvalidPrimaryThresh', ['Primary threshold must be >=' num2str(ST_Ut) ... ' and >0 and <=' num2str(pU_ST_Ut_filt) ]); end else pU_C_MaxT = 1-spm_Tcdf(C_MaxT,df); if ~((primaryThresh>=ST_Ut && primaryThreshpU_C_MaxT && primaryThresh<=pU_ST_Ut_filt)) - error(['Primary threshold must be >=' num2str(ST_Ut) ... + error('SnPM:InvalidPrimaryThresh', ['Primary threshold must be >=' num2str(ST_Ut) ... ' and <' num2str(C_MaxT) ' or >' num2str(pU_C_MaxT) ... ' and <= ' num2str(pU_ST_Ut_filt)]); end @@ -241,7 +243,7 @@ function snpm_combo_pp(CWD, job) % Theta = job.Thr.Clus.ClusMass.Theta; if (Theta<0+tol || Theta>1-tol) - error('Theta should be between 0 and 1'); + error('SnPM:InvalidTheta', 'Theta should be between 0 and 1'); end mTheta = Theta/(1-Theta); %-Weight for mass combining @@ -385,7 +387,7 @@ function snpm_combo_pp(CWD, job) %-Check XYZ for points > ST_Ut in perm 1 matches % XYZ computed above for SnPMt > ST_Ut if ~all(all( SnPM_ST(1:3,SnPM_ST(5,:)==1) == XYZ )) - error('ST XYZ don''t match between STCS & thresh') + error('SnPM:InvalidSTXYZ', 'ST XYZ don''t match between STCS & thresh') end %-- VOXEL CLUSTER COMBINED TEST STATISTICS ---- diff --git a/snpm_cp.m b/snpm_cp.m index 1a82799..e322346 100644 --- a/snpm_cp.m +++ b/snpm_cp.m @@ -214,16 +214,17 @@ function snpm_cp(CWD) %----------------------------------------------------------------------- load(CfgFile); if isempty([H C]) - error('No model specified; [H C] empty'); + error('SnPM:NoModel', 'No model specified; [H C] empty'); end if ~isempty(H) & ~isempty(C) - error('Cannot have both heirachical and covariate effects'); + error('SnPM:HierarchicalAndCov', 'Cannot have both heirachical and covariate effects'); end if size(CONT,2) ~= size([H C B G],2) - error('Contrast problem; wrong number of columns'); + error('SnPM:InvalidContrast','Contrast problem; wrong number of columns'); end if size(CONT,1) > 1 - warning('F contrast! F statistic images are being created.'); + warning('SnPM:FContrast', ... + 'F contrast! F statistic images are being created.'); STAT = 'F'; if (CONT(1,:) == -CONT(2,:)) CONT = CONT(1,:); @@ -236,23 +237,25 @@ function snpm_cp(CWD) CONT = full(u*sqrt(s))'; end if ~bVolm & bVarSm & vFWHM(3) - error('Cannot z-smooth variance in non-volumetric mode'); + error('SnPM:ZSmoothVolume', 'Cannot z-smooth variance in non-volumetric mode'); end if exist('bVarAlph')~=1 bVarAlph=0; end if bVarAlph & ~(~bVarSm & bVolm) - error('No pseudo t or nonvolumetric w/ variable alpha'); + error('SnPM:AlphaVolumePseudo', 'No pseudo t or nonvolumetric w/ variable alpha'); end if ~bVolm & pU_ST_Ut>=0 - error('Must work volumetrically to computer STCS on-the-fly'); + error('SnPM:STCSNotVolume', 'Must work volumetrically to computer STCS on-the-fly'); end -% Re-map files to avoid Endian headaches +% Re-map files to avoid Endian headaches; note if NaN's available +NaNrep=0; for i = 1:length(V) curr_pinfo = V(i).pinfo;% Added to keep scaling V(i) = spm_vol([V(i).fname ',' num2str(V(i).n)]); original_pinfo = V(i).pinfo; V(i).pinfo = curr_pinfo; + NaNrep = NaNrep | spm_type(V(i).dt(1),'nanrep'); end %-Delete files from previous analyses, if they exist @@ -276,10 +279,8 @@ function snpm_cp(CWD) %-Suprathreshold parameters %----------------------------------------------------------------------- -global SnPMdefs -if isempty(SnPMdefs), snpm_defaults; end -STalpha = SnPMdefs.STalpha; -STprop = SnPMdefs.STprop; +STalpha = snpm_get_defaults('STalpha'); +STprop = snpm_get_defaults('STprop'); s_SnPM_save = [s_SnPM_save ' STalpha STprop']; % Save for PP @@ -497,8 +498,12 @@ function snpm_cp(CWD) %-Eliminate background voxels (based on threshold TH), and % eliminate voxels where there are no differences across scans. %--------------------------------------------------------------------- - Q = find(all(X>TH) & any(diff(X)) & Wt); - + if ImMASK & NaNrep==0 + Q = find(all(X>TH) & any(diff(X)) & Wt & all(X~=0)); + else + Q = find(all(X>TH) & any(diff(X)) & Wt); + end + if length(Q) X = X(:,Q); @@ -689,7 +694,7 @@ function snpm_cp(CWD) end % Make an error if actually 'no voxels in brain'. -if perm==0, error('No voxels in brain'); end +if perm==0, error('SnPM:NoVoxelsInBrain', 'No voxels in brain'); end save SnPMt SnPMt @@ -727,7 +732,11 @@ function snpm_cp(CWD) if (pU_ST_Ut>1) ST_Ut=pU_ST_Ut; else - ST_Ut=spm_invTcdf(1-pU_ST_Ut, df); + if STAT == 'T' + ST_Ut = spm_invTcdf(1-pU_ST_Ut, df); + else + ST_Ut = spm_invFcdf(1-pU_ST_Ut, df1, df); + end end end end @@ -973,7 +982,7 @@ function snpm_cp(CWD) if diffWithRounded > tolerance Locs_vox_alter = MAT\Locs_mm; diffWithRounded_alter = max(abs(Locs_vox_alter(:)-round(Locs_vox(:)))); - error(['''Locs_vox'' must be integers (difference is ' num2str(diffWithRounded) ... + error('SnPM:NonIntegerLocs', ['''Locs_vox'' must be integers (difference is ' num2str(diffWithRounded) ... ' or ' num2str(diffWithRounded_alter) ')']); else Locs_vox = round(Locs_vox); diff --git a/snpm_defaults.m b/snpm_defaults.m index 8733315..2c4a69e 100644 --- a/snpm_defaults.m +++ b/snpm_defaults.m @@ -49,3 +49,10 @@ '0' }; +% If true, shuffles the seed of the random number generator to get +% different results every time. Use false, if you want to specify your own +% seed, for instance to insure that results can be replicated or when using +% a high performance cluster. +%------------------------------------------------------------------------ +SnPMdefs.shuffle_seed = true; + diff --git a/snpm_init.m b/snpm_init.m index 5825366..46f53e2 100644 --- a/snpm_init.m +++ b/snpm_init.m @@ -10,7 +10,7 @@ rev = '$Rev: 1716 $'; %#ok if exist('snpm_defaults')~=2 - error('Cannot initialize - SnPM toolbox not found in path') + error('SnPM:MissingSnPM', 'Cannot initialize - SnPM toolbox not found in path') end if isempty(whos('global','SnPMdefs')) snpm_defaults diff --git a/snpm_pi_ANOVAbtwnS.m b/snpm_pi_ANOVAbtwnS.m index cc53e23..8d890f9 100644 --- a/snpm_pi_ANOVAbtwnS.m +++ b/snpm_pi_ANOVAbtwnS.m @@ -102,9 +102,15 @@ %%% nCond = 2; % Number of conditions (groups) iGloNorm = '123'; % Allowable Global norm. codes sDesSave = 'iCond GrpCnt'; % PlugIn variables to save in cfg file -global TEST; -if isempty(TEST) || ~TEST % When testing we need a fixed seed - rand('seed',sum(100*clock)); % Initialise random number generator + +if snpm_get_defaults('shuffle_seed') + % Shuffle seed of random number generator + try + rng('shuffle'); + catch + % Old syntax + rand('seed',sum(100*clock)); + end end %-Get filenames and iCond, the condition labels @@ -129,7 +135,7 @@ % end % end -if nCond>255, error('Can''t support more than 255 groups'); end +if nCond>255, error('SnPM:TooManyGroups', 'Can''t support more than 255 groups'); end tmp0='A/B/...'; @@ -259,7 +265,7 @@ row_total = tmp*GrpCnt'; if ~all(all(double(PiCond)*ones(nScan,1)==row_total)) - error('Invalid PiCond computed!'), end + error('SnPM:InvalidPiCond', 'Invalid PiCond computed!'), end %-Find (maybe) iCond in PiCond, move iCond to 1st; %----------------------------------------------------------------------- @@ -287,7 +293,7 @@ PiCond(1,:) = iCond; perm = 1; else - error(['Bad PiCond (' num2str(perm) ')']) + error('SnPM:InvalidPiCond', ['Bad PiCond (' num2str(perm) ')']) end diff --git a/snpm_pi_ANOVAwithinS.m b/snpm_pi_ANOVAwithinS.m index 642dd44..5b7d8a0 100644 --- a/snpm_pi_ANOVAwithinS.m +++ b/snpm_pi_ANOVAwithinS.m @@ -108,21 +108,27 @@ %-Initialisation %----------------------------------------------------------------------- -global TEST; iGloNorm = '123'; % Allowable Global norm. codes sDesSave = 'iRepl sHCform_Mtx'; % PlugIn variables to save in cfg file -if isempty(TEST) || ~TEST - rand('seed',sum(100*clock)); % Initialise random number generator + +if snpm_get_defaults('shuffle_seed') + % Shuffle seed of random number generator + try + rng('shuffle'); + catch + % Old syntax + rand('seed',sum(100*clock)); + end end %-Get number of subjects nSubj = size(job.fsubject,2);%spm_input('# subjects','+1'); -if (nSubj==1), error('Use single subject plug for single subjects'); end +if (nSubj==1), error('SnPM:SingleSubj', 'Use single subject plug for single subjects'); end %-Get number of scans per subject - nSubj x nRepl design nRepl = unique(arrayfun(@(x) numel(x.scans), job.fsubject));%spm_input('# scans per subject','+1'); if numel(nRepl) > 1 - error('All subjects must have the same number of replications') + error('SnPM:DifferentReplications', 'All subjects must have the same number of replications') end @@ -267,7 +273,7 @@ PiCond(1,:) = iCond; perm = 1; else - error(['Bad PiCond (' num2str(perm) ')']) + error('SnPM:InvalidPiCond', ['Bad PiCond (' num2str(perm) ')']) end diff --git a/snpm_pi_Corr.m b/snpm_pi_Corr.m index e11bc7c..d7d5499 100644 --- a/snpm_pi_Corr.m +++ b/snpm_pi_Corr.m @@ -74,7 +74,15 @@ iGloNorm = '123'; % Allowable Global norm. codes sDesSave = 'CovInt'; % PlugIn variables to save in cfg file -global TEST; +if snpm_get_defaults('shuffle_seed') + % Shuffle seed of random number generator + try + rng('shuffle'); + catch + % Old syntax + rand('seed',sum(100*clock)); + end +end %-Get filenames of scans %----------------------------------------------------------------------- @@ -90,7 +98,7 @@ if BATCH CovInt = job.CovInt(:); if ~all(size(CovInt)==[nSubj,1]) - error(sprintf('Covariate [%d,1] doesn''t match number of subjects [%d]',... + error('SnPM:InvalidCovariate', sprintf('Covariate [%d,1] doesn''t match number of subjects [%d]',... size(CovInt,1),nSubj)) end else @@ -117,7 +125,7 @@ end nGcs = size(Gc,2); if size(d,1) ~= q - error(sprintf('Covariate [%d,1] does not match number of subjects [%d]',... + error('SnPM:InvalidCovariate', sprintf('Covariate [%d,1] does not match number of subjects [%d]',... size(job.cov(i).c,1),nSubj)) else %-Save raw covariates for printing later on @@ -181,9 +189,6 @@ %-Approximate test : % Build up random subset of all (within nSubj) permutations %=============================================================== - if isempty(TEST) || ~TEST % When testing the code we need a fixed seed - rand('seed',sum(100*clock)) %-Initialise random number generator - end PiCond = zeros(nPiCond,nSubj); PiCond(1,:) = 1+rem([0:nSubj-1],nSubj); for i = 2:nPiCond @@ -235,17 +240,14 @@ %----------------------------------------------------------------------- %-Check PiConds sum within Xblks to sum to 1 if ~all(all(sum(PiCond,2)== (nSubj+1)*nSubj/2 )) - error('Invalid PiCond computed!'), end + error('SnPM:InvalidPiCond', 'Invalid PiCond computed!'), end %-Convert to full permutations from permutations within blocks nPiCond = size(PiCond,1); %-Randomise order of PiConds (except first) to allow interim analysis -if isempty(TEST) || ~TEST % When testing the code we need a fixed seed - rand('seed',sum(100*clock)) %-Initialise random number generator -end PiCond=[PiCond(1,:);PiCond(randperm(nPiCond-1)+1,:)]; %-Check first permutation is null permutation if ~all(PiCond(1,:)==[1:nSubj]) - error('PiCond(1,:)~=[1:nSubj]'); end + error('SnPM:InvalidFirstPiCond', 'PiCond(1,:)~=[1:nSubj]'); end %-Form non-null design matrix partitions (Globals handled later) diff --git a/snpm_pi_Corr1S.m b/snpm_pi_Corr1S.m index 7a22466..e527e97 100644 --- a/snpm_pi_Corr1S.m +++ b/snpm_pi_Corr1S.m @@ -79,11 +79,18 @@ %-Initialisation %----------------------------------------------------------------------- -global TEST; iGloNorm = '123'; % Allowable Global norm. codes sDesSave = 'CovInt Xblk'; % PlugIn variables to save in cfg file - +if snpm_get_defaults('shuffle_seed') + % Shuffle seed of random number generator + try + rng('shuffle'); + catch + % Old syntax + rand('seed',sum(100*clock)); + end +end %-Get filenames of scans %----------------------------------------------------------------------- @@ -114,7 +121,7 @@ end nGcs = size(Gc,2); if size(d,1) ~= q - error(sprintf('Covariate [%d,1] does not match number of scans [%d]',... + error('SnPM:InvalidCovariate', sprintf('Covariate [%d,1] does not match number of scans [%d]',... size(job.cov(i).c,1),nScan)) else %-Save raw covariates for printing later on @@ -176,9 +183,6 @@ %-Approximate test : % Build up random subset of all (within Xblk) permutations %=============================================================== - if isempty(TEST) || ~TEST - rand('seed',sum(100*clock)) %-Initialise random number generator - end PiCond = zeros(nPiCond,nScan); PiCond(1,:) = 1+rem([0:Xblk*nXblk-1],Xblk); for i = 2:nPiCond @@ -241,18 +245,15 @@ %----------------------------------------------------------------------- %-Check PiConds sum within Xblks to sum of first nXblk natural numbers if ~all(all(PiCond*spm_DesMtx(iXblk)== (Xblk+1)*Xblk/2 )) - error('Invalid PiCond computed!'), end + error('SnPM:InvalidPiCond', 'Invalid PiCond computed!'), end %-Convert to full permutations from permutations within blocks nPiCond = size(PiCond,1); PiCond = PiCond + meshgrid((iXblk-1)*Xblk,1:nPiCond); %-Randomise order of PiConds (except first) to allow interim analysis -if isempty(TEST) || ~TEST - rand('seed',sum(100*clock)) %-Initialise random number generator -end PiCond=[PiCond(1,:);PiCond(randperm(nPiCond-1)+1,:)]; %-Check first permutation is null permutation if ~all(PiCond(1,:)==[1:nScan]) - error('PiCond(1,:)~=[1:nScan]'); end + error('SnPM:InvalidPiCond', 'PiCond(1,:)~=[1:nScan]'); end %-Form non-null design matrix partitions (Globals handled later) diff --git a/snpm_pi_OneSampT.m b/snpm_pi_OneSampT.m index 2359b2f..53db09d 100644 --- a/snpm_pi_OneSampT.m +++ b/snpm_pi_OneSampT.m @@ -94,9 +94,15 @@ %----------------------------------------------------------------------- iGloNorm = '123'; % Allowable Global norm. codes sDesSave = 'iCond'; % PlugIn variables to save in cfg file -global TEST; -if isempty(TEST) || ~TEST % When testing the code we need a fixed seed - rand('seed',sum(100*clock)); % Initialise random number generator + +if snpm_get_defaults('shuffle_seed') + % Shuffle seed of random number generator + try + rng('shuffle'); + catch + % Old syntax + rand('seed',sum(100*clock)); + end end %-Get filenames and iCond, the condition labels @@ -118,7 +124,7 @@ end nGcs = size(Gc,2); if size(d,1) ~= q - error(sprintf('Covariate [%d,1] does not match number of subjects [%d]',... + error(sprintf('SnPM:InvalidCovariate', 'Covariate [%d,1] does not match number of subjects [%d]',... size(job.cov(i).c,1),nScan)) else %-Save raw covariates for printing later on @@ -244,7 +250,7 @@ PiCond(1,:) = iCond; perm = 1; else - error(['Bad PiCond (' num2str(perm) ')']) + error('SnPM:PiCond', ['Bad PiCond (' num2str(perm) ')']) end diff --git a/snpm_pi_PairT.m b/snpm_pi_PairT.m index 05c5eb6..f18e88b 100644 --- a/snpm_pi_PairT.m +++ b/snpm_pi_PairT.m @@ -112,7 +112,7 @@ %-Get number of subjects % nSubj = spm_input('# subjects','+1'); -% if (nSubj==1), error('Use single subject plug for single subjects'); end +% if (nSubj==1), error('SnPM:SingleSubject', 'Use single subject plug for single subjects'); end nSubj = numel(job.fsubject); %-Only consider one replication -- basically a RFX machine. @@ -144,23 +144,23 @@ %-Check validity of tmpCond if length(tmpCond)==nScan if length(find(diff(sort(tmpCond)))) ~= nCond-1 - error('Exactly ',[int2str(nCond), ... + error('SnPM:InvalidnCond', 'Exactly ',[int2str(nCond), ... ' conditions must be supplied']); elseif sum(tmpCond)~=0 - error(['Exactly ',int2str(nRepl),' As and ', ... + error('SnPM:InvalidnRepl', ['Exactly ',int2str(nRepl),' As and ', ... int2str(nRepl),' Bs must be supplied']); elseif isempty(iCond) sCond=setstr(tmp([1,diff(sort(tmp))]~=0)); Cond = tmpCond; elseif any(iCond(1:nScan)~=tmpCond) & ... any(iCond(1:nScan)~=(-tmpCond)) - error(['Conditions index must be same as', ... + error('SnPM:InvalidiCond', ['Conditions index must be same as', ... 'first subject, or flipped']); else Cond = tmpCond; end else - error(['Enter indicies for ',int2str(nCond*nRepl),' scans']) + error('SnPM:InvalidiIndices', ['Enter indicies for ',int2str(nCond*nRepl),' scans']) end end iCond = [iCond, Cond]; @@ -197,7 +197,7 @@ bAproxTst=1; end if rem(nPiSubj,2) - error(['Number of perms must be even']); + error('SnPM:OddPermutations', ['Number of perms must be even']); nPiSubj = 0; end @@ -208,7 +208,7 @@ % tmp = spm_input(sprintf('# perms. to use? (Max %d)',nPiSubj),'+0'); % tmp = floor(max([0,tmp])); % if rem(tmp,2) -% error(['Number of perms must be even']); +% error('SnPM:OddPermutations', ['Number of perms must be even']); % tmp=0; % end % end @@ -241,7 +241,7 @@ % Look for correct labeling d = find(all((PiSubj==meshgrid(iSubjC,1:size(PiSubj,1)))')); if (length(d)~=1 & ~bAproxTst) - error('Internal error: Correct labeling is not in the perms'); + error('SnPM:CorrectLabelMissing', 'Internal error: Correct labeling is not in the perms'); elseif (length(d)~=1) % Correct labeling randomly removed, insert at top PiSubj(1,:) = iSubjC; diff --git a/snpm_pi_PairTrand.m b/snpm_pi_PairTrand.m index e0d1ced..b5cd1ad 100644 --- a/snpm_pi_PairTrand.m +++ b/snpm_pi_PairTrand.m @@ -101,7 +101,7 @@ %-Get number of subjects nSubj = spm_input('# subjects','+1'); -if (nSubj==1), error('Use single subject plug for single subjects'); end +if (nSubj==1), error('SnPM:SingleSubject','Use single subject plug for single subjects'); end %-Get number of replications per condition - 2 x nRepl design nRepl = spm_input('# replications per condition','+1'); @@ -130,23 +130,23 @@ %-Check validity of tmpCond if length(tmpCond)==nScan if length(find(diff(sort(tmpCond)))) ~= nCond-1 - error('Exactly ',[int2str(nCond), ... + error('SnPM:InvalidnCond', 'Exactly ',[int2str(nCond), ... ' conditions must be supplied']); elseif sum(tmpCond)~=0 - error(['Exactly ',int2str(nRepl),' As and ', ... + error('SnPM:InvalidnRepl',['Exactly ',int2str(nRepl),' As and ', ... int2str(nRepl),' Bs must be supplied']); elseif isempty(iCond) sCond=setstr(tmp([1,diff(sort(tmp))]~=0)); Cond = tmpCond; elseif any(iCond(1:nScan)~=tmpCond) & ... any(iCond(1:nScan)~=(-tmpCond)) - error(['Conditions index must be same as ', ... + error('SnPM:InvalidiCond',['Conditions index must be same as ', ... 'first subject, or flipped']); else Cond = tmpCond; end else - error(['Enter indicies for ',int2str(nCond*nRepl),' scans']) + error('SnPM:InvalidnIndices',['Enter indicies for ',int2str(nCond*nRepl),' scans']) end end iCond = [iCond, Cond]; @@ -171,7 +171,8 @@ % Balanced design... everything is cool else % not balanced - warning('Design not balanced; results may be confounded with time'); + warning('SnPM:Unbalanced', ... + 'Design not balanced; results may be confounded with time'); end %-Compute permutations of subjects @@ -186,7 +187,7 @@ tmp = spm_input(sprintf('# perms. to use? (Max %d)',nPiSubj),'+0'); tmp = floor(max([0,tmp])); if rem(tmp,2) - error(['Number of perms must be even']); + error('SnPM:OddPermutations',['Number of perms must be even']); tmp=0; end end @@ -216,7 +217,7 @@ % Move correct labeling to top d = find(all((PiSubj==meshgrid(iSubjC(2:nSubj),1:size(PiSubj,1)))')); if (length(d)~=1) - error('Internal error: Correct labeling is not in the perms'); end + error('SnPM:CorrectLabelMissing', 'Internal error: Correct labeling is not in the perms'); end PiSubj(d,:) = PiSubj(1,:); PiSubj(1,:) = iSubjC(2:nSubj); % Pick a random sample of PiSubj if bAprxTst diff --git a/snpm_pi_TwoSampPairT.m b/snpm_pi_TwoSampPairT.m index 929f3d8..77063e4 100644 --- a/snpm_pi_TwoSampPairT.m +++ b/snpm_pi_TwoSampPairT.m @@ -110,9 +110,15 @@ nStud = 2; % Number of groups iGloNorm = '123'; % Allowable Global norm. codes sDesSave = 'iStud iCond nSubj'; % PlugIn variables to save in cfg file -global TEST; -if isempty(TEST) || ~TEST - rand('seed',sum(100*clock)); % Initialise random number generator + +if snpm_get_defaults('shuffle_seed') + % Shuffle seed of random number generator + try + rng('shuffle'); + catch + % Old syntax + rand('seed',sum(100*clock)); + end end iStudC = []; % +1/-1 version of iStud @@ -277,7 +283,7 @@ %-Check each perm in PiStuds sums to nStud1-nStud2 if ~all(all(PiStud*ones(nSUBJ,1)==nSubj*[1 -1]')) - error('Invalid PiStud computed!'), end + error('SnPM:InvalidPiStud', 'Invalid PiStud computed!'), end %-Find (maybe) iStudC in PiStud, move iStudC to 1st; negate if neccesary %----------------------------------------------------------------------- @@ -303,7 +309,7 @@ PiStud(1,:) = iStudC; perm = 1; else - error(['Bad PiStud (' num2str(perm) ')']) + error('SnPM:InvalidPiStud', ['Bad PiStud (' num2str(perm) ')']) end %-Turn PiStud into PiCond diff --git a/snpm_pi_TwoSampT.m b/snpm_pi_TwoSampT.m index b82861a..086d685 100644 --- a/snpm_pi_TwoSampT.m +++ b/snpm_pi_TwoSampT.m @@ -111,10 +111,19 @@ nCond = 2; % Number of conditions (groups) iGloNorm = '123'; % Allowable Global norm. codes sDesSave = 'iCond GrpCnt'; % PlugIn variables to save in cfg file -global TEST; -if isempty(TEST) || ~TEST % When testing the code we need a fixed seed - rand('seed',sum(100*clock)); % Initialise random number generator + +if snpm_get_defaults('shuffle_seed') + % Shuffle seed of random number generator + try + rng('shuffle'); + catch + % Old syntax + rand('seed',sum(100*clock)); + end end +nPermMC = 2000; % When using more than this many permutations, use + % conventional MC permutation test, i.e. don't try + % to excluded repeated permutations. %-Get filenames and iCond, the condition labels %======================================================================= @@ -218,8 +227,10 @@ % Fill subsequent rows, checking that we're not repeating for i=2:nPiCond tmp=PiCond(i-1,randperm(nScan)); - while any(all(PiCond(1:(i-1),:)'==meshgrid(tmp,1:(i-1))')) - tmp=PiCond(i-1,randperm(nScan)); + if nPerm<=nPermMC + while any(all(PiCond(1:(i-1),:)'==meshgrid(tmp,1:(i-1))')) + tmp=PiCond(i-1,randperm(nScan)); + end end PiCond(i,:)=tmp; end @@ -229,7 +240,7 @@ %-Check PiConds sum to nGrp1-nGrp2 if ~all(all(PiCond*ones(nScan,1)==nScan-2*nFlip)) - error('Invalid PiCond computed!'), end + error('SnPM:InvalidPiCond', 'Invalid PiCond computed!'), end %-Find (maybe) iCond in PiCond, move iCond to 1st; negate if neccesary %----------------------------------------------------------------------- @@ -255,7 +266,7 @@ PiCond(1,:) = iCond; perm = 1; else - error(['Bad PiCond (' num2str(perm) ')']) + error('SnPM:InvalidPiCond', ['Bad PiCond (' num2str(perm) ')']) end diff --git a/snpm_pi_TwoSampTss.m b/snpm_pi_TwoSampTss.m index 504b1f9..771f557 100644 --- a/snpm_pi_TwoSampTss.m +++ b/snpm_pi_TwoSampTss.m @@ -78,7 +78,6 @@ %-Initialisation %----------------------------------------------------------------------- -global TEST; nCond = 2; % Number of conditions iGloNorm = '123'; % Allowable Global norm. codes sDesSave = 'iCond iRepl Xblk'; % PlugIn variables to save in cfg file @@ -86,6 +85,16 @@ %-Get number of replications per condition - 2 x nRepl design nRepl = job.Tss_repc;%spm_input('# replications per condition','+1'); +if snpm_get_defaults('shuffle_seed') + % Shuffle seed of random number generator + try + rng('shuffle'); + catch + % Old syntax + rand('seed',sum(100*clock)); + end +end + %-Work out exchangability blocks - Assumme Xblks of equal size %----------------------------------------------------------------------- %-Valid sizes of X-blk (Xblk) are integer divisors of nCond*nRepl that @@ -129,7 +138,7 @@ end %-Check PiConds sum to zero within Xblks if ~all(all(PiCond*spm_DesMtx(iXblk)==0)) - error('Invalid PiCond computed!'), end + error('SnPM:InvalidPiCond', 'Invalid PiCond computed!'), end %-Take half of PiCond. This is always possible. %----------------------------------------------------------------------- @@ -156,10 +165,10 @@ perm = find(all((meshgrid(iCond,1:size(PiCond,1))==PiCond)')); if (bhPerms), perm=[perm,... -find(all((meshgrid(iCond,1:size(PiCond,1))==-PiCond)'))]; end - if length(perm)>1, error('Multiple iCond in PiCond'), end - if isempty(perm), error('Invalid iCond for this design'), end + if length(perm)>1, error('SnPM:InvalidPiCond', 'Multiple iCond in PiCond'), end + if isempty(perm), error('SnPM:InvalidiCond', 'Invalid iCond for this design'), end else - error(['Enter indicies for ',int2str(nCond*nRepl),' scans']) + error('SnPM:InvalidIndices', ['Enter indicies for ',int2str(nCond*nRepl),' scans']) end end @@ -172,16 +181,12 @@ %----------------------------------------------------------------------- if (perm<0), PiCond=-PiCond; perm=-perm; end %-Lame last ditch check just to make sure I know what's going on! **** -if ~all(iCond==PiCond(perm,:)), error('iCond~=PiCond(perm,:)'), end +if ~all(iCond==PiCond(perm,:)), error('SnPM:InvalidiCond', 'iCond~=PiCond(perm,:)'), end %-Actual labelling must be at top of PiCond if (perm~=1) PiCond(perm,:)=[]; PiCond=[iCond;PiCond]; end -if isempty(TEST) || ~TEST - %-Randomise order of PiConds (except first) to allow interim analysis - rand('seed',sum(100*clock)) %-Initialise random number generator -end PiCond=[PiCond(1,:);PiCond(randperm(size(PiCond,1)-1)+1,:)]; diff --git a/snpm_pp.m b/snpm_pp.m index 93e42ef..5f572c9 100644 --- a/snpm_pp.m +++ b/snpm_pp.m @@ -512,19 +512,19 @@ function snpm_pp(CWD,varargin) else % Cluster-wise inference if exist(fullfile(CWD,'SnPM_ST.mat'))~=2 & exist(fullfile(CWD,'STCS.mat'))~=2 - error(['ERROR: Cluster-wise inference requested, but no cluster information saved.\n',... + error(['SnPM:NoClusterInfo', 'ERROR: Cluster-wise inference requested, but no cluster information saved.\n',... 'Re-configure analysis changing "Cluster inference" to "Yes" and re-run.\n']) end %%% Sort out the cluster-forming threshold if pU_ST_Ut==-1 % No threshold was set in snpm_ui. if isnan(job.Thr.Clus.ClusSize.CFth) - error('ERROR: Cluster-forming threshold set to NaN in results with "slow" cluster inference method used in compoutation. \nRe-run results specifying a cluster-forming threshold.\n') + error('SnPM:NoClusterFormingThresh', 'ERROR: Cluster-forming threshold set to NaN in results with "slow" cluster inference method used in compoutation. \nRe-run results specifying a cluster-forming threshold.\n') end % Save original ST_Ut ST_Ut_0 = ST_Ut; CFth=job.Thr.Clus.ClusSize.CFth; if (CFth<=0) - error('ERROR: Cluster-forming threshold must be strictly positive.\nRe-run results with a cluster-forming threshold greater than 0.\n') + error('SnPM:InvalidClusterFormingThresh', 'ERROR: Cluster-forming threshold must be strictly positive.\nRe-run results with a cluster-forming threshold greater than 0.\n') end if bVarSm %-If using pseudo-statistics then can't use (uncorrected) @@ -544,9 +544,9 @@ function snpm_pp(CWD,varargin) if (CFth < ST_Ut)%(CFth>=ST_Ut-tol) if isnan(pCFth) - error(sprintf('ERROR: Cluster-forming threshold of %0.2f specified, but statistic image information only saved for %0.2f and greater. \nRe-run results with a cluster-forming threshold of %0.2f or higher. (Alternatively, increase SnPMdefs.STprop in snpm_defaults.m, re-start SnPM, and re-compute analysis.)\n',CFth,ST_Ut,ST_Ut)) + error('SnPM:InvalidClusterFormingThresh', sprintf('ERROR: Cluster-forming threshold of %0.2f specified, but statistic image information only saved for %0.2f and greater. \nRe-run results with a cluster-forming threshold of %0.2f or higher. (Alternatively, increase SnPMdefs.STprop in snpm_defaults.m, re-start SnPM, and re-compute analysis.)\n',CFth,ST_Ut,ST_Ut)) else - error(sprintf('ERROR: Cluster-forming threshold of P=%0.4f (T=%0.2f) specified, but statistic image information only saved for %0.2f and greater. \nRe-run results with a cluster-forming P-value threshold of %0.2f or lower. (Alternatively, increase SnPMdefs.STalpha in snpm_defaults.m, re-start SnPM, and re-compute analysis.)\n',pCFth,CFth,ST_Ut,p_ST_Ut)) + error('SnPM:InvalidClusterFormingThresh', sprintf('ERROR: Cluster-forming threshold of P=%0.4f (T=%0.2f) specified, but statistic image information only saved for %0.2f and greater. \nRe-run results with a cluster-forming P-value threshold of %0.2f or lower. (Alternatively, increase SnPMdefs.STalpha in snpm_defaults.m, re-start SnPM, and re-compute analysis.)\n',pCFth,CFth,ST_Ut,p_ST_Ut)) end end else @@ -564,9 +564,9 @@ function snpm_pp(CWD,varargin) if (CFth < ST_Ut) %(CFth>=ST_Ut-tol) if isnan(pCFth) % statistic-value cluster-forming threshold - error(sprintf('ERROR: Cluster-forming threshold of %0.2f specified, but statistic image information only saved for %0.2f and greater. \nRe-run results with a cluster-forming threshold of %0.2f or higher. (Alternatively, increase SnPMdefs.STalpha in snpm_defaults.m, re-start SnPM, and re-compute analysis.)\n',CFth,ST_Ut,ST_Ut)) + error('SnPM:InvalidClusterFormingThresh', sprintf('ERROR: Cluster-forming threshold of %0.2f specified, but statistic image information only saved for %0.2f and greater. \nRe-run results with a cluster-forming threshold of %0.2f or higher. (Alternatively, increase SnPMdefs.STalpha in snpm_defaults.m, re-start SnPM, and re-compute analysis.)\n',CFth,ST_Ut,ST_Ut)) else - error(sprintf('ERROR: Cluster-forming threshold of P=%0.4f (T=%0.2f) specified, but statistic image information only saved for %0.2f and greater. \nRe-run results with a cluster-forming P-value threshold of %0.2f or lower. (Alternatively, increase SnPMdefs.STalpha in snpm_defaults.m, re-start SnPM, and re-compute analysis.)\n',pCFth,CFth,ST_Ut,p_ST_Ut)) + error('SnPM:InvalidClusterFormingThresh', sprintf('ERROR: Cluster-forming threshold of P=%0.4f (T=%0.2f) specified, but statistic image information only saved for %0.2f and greater. \nRe-run results with a cluster-forming P-value threshold of %0.2f or lower. (Alternatively, increase SnPMdefs.STalpha in snpm_defaults.m, re-start SnPM, and re-compute analysis.)\n',pCFth,CFth,ST_Ut,p_ST_Ut)) end end end @@ -576,7 +576,7 @@ function snpm_pp(CWD,varargin) ST_Ut = CFth; else % Threshold *was* set in snpm_ui. if ~isnan(job.Thr.Clus.ClusSize.CFth) - error(sprintf('ERROR: Cluster-forming threshold of T=%0.2f was already set during analysis configuration; hence, in results, cluster-forming threshold must be left as "NaN".\nRe-run results with cluster-forming threshold set to NaN.\n',ST_Ut)) + error('SnPM:InvalidClusterFormingThresh', sprintf('ERROR: Cluster-forming threshold of T=%0.2f was already set during analysis configuration; hence, in results, cluster-forming threshold must be left as "NaN".\nRe-run results with cluster-forming threshold set to NaN.\n',ST_Ut)) end end u=ST_Ut; % Flag use of a statistic-value threshold @@ -823,7 +823,8 @@ function snpm_pp(CWD,varargin) load(fullfile(CWD,'SnPM_ST')) catch exception if strcmp(exception.message, 'File may be corrupt.') - warning(['SnPM_ST file can not be loaded. Consider using' ... + warning('SnPM:SnPMSTFileNotLOaded', ... + ['SnPM_ST file can not be loaded. Consider using' ... ' ''set cluster-forming threshold now (fast)'' option' ... ' in SnPM ''Specify''.']); end @@ -883,7 +884,7 @@ function snpm_pp(CWD,varargin) if diffWithRounded > tolerance Locs_vox_alter = MAT\Locs_mm; diffWithRounded_alter = max(abs(Locs_vox_alter(:)-round(Locs_vox(:)))); - error(['''Locs_vox'' must be integers (difference is ' num2str(diffWithRounded) ... + error('SnPM:NonIntegerLocsvox', ['''Locs_vox'' must be integers (difference is ' num2str(diffWithRounded) ... ' or ' num2str(diffWithRounded_alter) ')']); else Locs_vox = round(Locs_vox); @@ -969,11 +970,11 @@ function snpm_pp(CWD,varargin) % XYZ computed above for SnPMt > ST_Ut %if pU_ST_Ut==-1 % if ~all(all( SnPM_ST(1:3,SnPM_ST(5,:)==1) == XYZ )) - % error('ST XYZ don''t match between STCS & thresh') + % error('SnPM:InvalidSTXYZ', 'ST XYZ don''t match between STCS & thresh') % end %else if ~all(all( STCstats(1:3,:) == XYZ )) - error('ST XYZ don''t match between STCS & thresh') + error('SnPM:InvalidSTXYZ', 'ST XYZ don''t match between STCS & thresh') end %end end @@ -1025,7 +1026,7 @@ function snpm_pp(CWD,varargin) % Thresholding based on uncorrected P-values Q = find(SnPMucp < alpha_ucp); else - error('Coding error') + error('SnPM:CodingError', 'Coding error') end if ~isnan(u) fprintf('Filtering statistic image, voxelwise...'); @@ -1140,8 +1141,9 @@ function snpm_pp(CWD,varargin) end end - - +% Display only if *not* in command line mode +if ~spm_get_defaults('cmdline') + %======================================================================= %-D I S P L A Y : Max report %======================================================================= @@ -1438,6 +1440,8 @@ function snpm_pp(CWD,varargin) end +end + %- Image output? %======================================================================= %-Write out filtered SnPMt? @@ -1525,7 +1529,7 @@ function ShowDist(T,cT,aT,C,cC,aC,Typ) elseif strcmp(Typ,'uncor') TitStr = 'Distribution: Observed Uncorrected P-values'; else - error('Unknown type string') + error('SnPM:UnknownString', 'Unknown type string') end pos1 = [0.125 0.50 0.75 0.3]; @@ -1614,12 +1618,12 @@ function ShowDist(T,cT,aT,C,cC,aC,Typ) title('Permutation Distribution: Maximum Cluster Size (cuberoot plot)','FontSize',14) xlabel(sprintf('Max Cluster Size Observed = %g Cluster Threshold = %g',C(1),cC)) else - warning('Uncorrected cluster size not supported yet!') + warning('SnPM:UncorrectedClustSize', 'Uncorrected cluster size not supported yet!') title('Permutation Distribution: (Uncorrected) Cuberoot Cluster Size','FontSize',14) xlabel(sprintf('Max Cluster Size Observed = %g Cluster Threshold = %g',C(1),cC)) end Ylim = get(gca,'ylim'); Xlim = get(gca,'xlim'); - set(gca,'Xticklabel',num2str(str2num(get(gca,'Xticklabel')).^3)) + set(gca,'Xticklabel',num2str(str2double(cellstr(get(gca,'Xticklabel'))).^3)) line(rC(1)*[1 1],Ylim.*[1 0.95],'LineStyle','--','color',[1 0 0]); text(rC(1)+diff(Xlim)*0.01,Ylim(2)*0.95,'Observed','FontSize',10) line(cC.^(1/3)*[1 1],Ylim.*[1 0.85],'LineStyle','-'); @@ -1643,7 +1647,7 @@ function ShowDist(T,cT,aT,C,cC,aC,Typ) function v = BoundCheck(val,Range,ErrMsg) if valRange(2) - error([ErrMsg ': ' num2str(val)]) + error('SnPM:RangeError', [ErrMsg ': ' num2str(val)]) else v=val; end diff --git a/snpm_t2z.m b/snpm_t2z.m index 0cb6b18..530c101 100644 --- a/snpm_t2z.m +++ b/snpm_t2z.m @@ -66,9 +66,9 @@ %-Argument range and size checks %--------------------------------------------------------------------------- -if nargin<2, error('insufficient arguments'), end -if length(df)~=1, error('df must be a scalar'), end -if df<=0, error('df must be strictly positive'), end +if nargin<2, error('SnPM:InsufficientArgs', 'insufficient arguments'), end +if length(df)~=1, error('SnPM:Invaliddf', 'df must be a scalar'), end +if df<=0, error('SnPM:Invaliddf', 'df must be strictly positive'), end %-Computation %=========================================================================== diff --git a/snpm_uc_FDR.m b/snpm_uc_FDR.m index 6cbf6f4..92d8fe0 100644 --- a/snpm_uc_FDR.m +++ b/snpm_uc_FDR.m @@ -103,7 +103,7 @@ %----------------------------------------------------------------------- if isstruct(Vs) if (abs(n) ~= length(Vs)) - error(sprintf('n & number of mapped images doesn''t match (%d,%d)',... + error('SnPM:Invalidn', sprintf('n & number of mapped images doesn''t match (%d,%d)',... abs(n),length(Vs))); end Ts = spm_read_vols(Vs(1)); diff --git a/snpm_ui.m b/snpm_ui.m index 64b7f42..221a485 100644 --- a/snpm_ui.m +++ b/snpm_ui.m @@ -122,7 +122,7 @@ function snpm_ui(varargin) % sDesign Description of PlugIn design % V Memory mapping handles % MASK Filename of explicit mask image -% ImMASKING Implicit masking; 0=none; 1=zeros are equivalent to NaN +% ImMASK Implicit masking; 0=none; 1=zeros are equivalent to NaN % % df degrees of freedom due to error % sDesSave String of PlugIn variables to save to cfg file @@ -282,7 +282,7 @@ function snpm_ui(varargin) vFWHM=[0,0,0]; % FWHM for variance smoothing sVarSm=''; % String describing Variance Smoothing bVolm=0; % Flag for volumetric computation -nMax4DefVol=SnPMdefs.nMax4DefVol; +nMax4DefVol=snpm_get_defaults('nMax4DefVol'); % Default to volumetric if less than nMax4DefVol scans sPiCond=''; % String describing permutations in PiCond bhPerms=0; % Flag for half permutations. Rest are then their inverses @@ -322,7 +322,8 @@ function snpm_ui(varargin) %-Decide upon volumetric operation bVolm = job.bVolm; if ~bVolm & (vFWHM(3)~=0) -warning(sprintf(['Working volumetrically because of smoothing in z (%g).\n'... +warning('SnPM:MayRunOutOfMemory', ... + sprintf(['Working volumetrically because of smoothing in z (%g).\n'... 'May run out of memory.'],vFWHM(3))); bVolm=1; end @@ -334,7 +335,8 @@ function snpm_ui(varargin) % Add: get primary threshold for STC analysis if requested if bST if ~bVolm - warning(sprintf('Note: Cannot define threshold now, because not working volumetrically\n')); + warning('SnPM:CannotDefineThreshVolumetrically', ... + sprintf('Note: Cannot define threshold now, because not working volumetrically\n')); pU_ST_Ut=-1; % Define the threshold later else pU_ST_Ut = ~isfield(job.ST,'ST_later'); @@ -398,7 +400,7 @@ function snpm_ui(varargin) GX = job.globalc.g_user.global_uval; rg = GX; if length(GX) ~= nScan - error(['User-specified globals length [%d] doesn''t match number of' ... + error('SnPM:InvalidGlobals', ['User-specified globals length [%d] doesn''t match number of' ... ' scans [%d]'],length(GX),nScan); end elseif isfield(job.globalc,'g_mean') @@ -431,7 +433,7 @@ function snpm_ui(varargin) %-Implicit Masking - Batch only! %----------------------------------------------------------------------- -ImMASKING=job.masking.im; +ImMASK=job.masking.im; %-Get analysis mask %----------------------------------------------------------------------- @@ -457,10 +459,10 @@ function snpm_ui(varargin) %-Check compatability of images (Bombs for single image) %----------------------------------------------------------------------- if any(any(diff(cat(1,V(:).dim),1,1),1)&[1,1,1]) - error('images do not all have the same dimensions') + error('SnPM:ImageDirections', 'images do not all have the same dimensions') end if any(any(any(diff(cat(3,V(:).mat),1,3),3))) - error('images do not all have same orientation & voxel size') + error('SnPM:ImageOrientationsResolutions','images do not all have same orientation & voxel size') end %-Get ORIGIN, etc @@ -479,6 +481,10 @@ function snpm_ui(varargin) %-Compute global values rg = zeros(nScan,1); for i = 1:nScan, rg(i)=spm_global(V(i)); end + if any(~isfinite(rg)) + disp(rg) + error('SnPM:NaNGlobal', 'Global computation returned NaN! Cannot continue') + end GX = rg; elseif iGXcalc==1 rg = []; @@ -508,7 +514,7 @@ function snpm_ui(varargin) % Relative threshold TH = THRESH * GX; else - error(['Wrong value for iTHRESH: ', num2str(iTHRESH)]) + error('SnPM:InvalidiTHRESH', ['Wrong value for iTHRESH: ', num2str(iTHRESH)]) end @@ -550,7 +556,7 @@ function snpm_ui(varargin) else Gnames = str2mat(Gnames,GLnames); end else %----------------------------------------------------------------------- - error(sprintf('%cError: invalid iGloNorm option\n',7)) + error('SnPM:InvalidiGloNorm', sprintf('%cError: invalid iGloNorm option\n',7)) end % (if) @@ -577,148 +583,150 @@ function snpm_ui(varargin) %----------------------------------------------------------------------- s_SnPMcfg_save = ['s_SnPMcfg_save H C B G HCBGnames P PiCond ',... 'sPiCond bhPerms sHCform iGloNorm sGloNorm GM rg GX GMscale CONT ',... - 'THRESH MASK TH bVarSm vFWHM sVarSm bVolm bST sDesFile sDesign ',... + 'THRESH MASK ImMASK TH bVarSm vFWHM sVarSm bVolm bST sDesFile sDesign ',... 'V pU_ST_Ut df1 ', ... 'sDesSave ',sDesSave]; eval(['save SnPMcfg ',s_SnPMcfg_save]) +if ~spm_get_defaults('cmdline') -%======================================================================= -%-Display parameters -%======================================================================= - -if BATCH - % OK, only _now_ activate the graphics... so if we crash here, there's - % nothing lost - Finter = spm_figure('FindWin','Interactive'); - Fgraph = spm_figure('FindWin','Graphics'); - if isempty(Fgraph), Fgraph=spm_figure('Create','Graphics'); end - spm_clf(Finter), spm_clf(Fgraph) -end - + %======================================================================= + %-Display parameters + %======================================================================= -%-Muck about a bit to set flags for various indicators - handy for later -bMStud=~isempty(iStud); -bMSubj=~isempty(iSubj); -bMCond=~isempty(iCond); -bMRepl=~isempty(iRepl); -bMXblk=~isempty(iXblk); + if BATCH + % OK, only _now_ activate the graphics... so if we crash here, there's + % nothing lost + Finter = spm_figure('FindWin','Interactive'); + Fgraph = spm_figure('FindWin','Graphics'); + if isempty(Fgraph), Fgraph=spm_figure('Create','Graphics'); end + spm_clf(Finter), spm_clf(Fgraph) + end -%-Compute common path components - all paths will begin with file separator -%----------------------------------------------------------------------- -d = max(find(P(1,1:find(~all(P == ones(nScan,1)*P(1,:)), 1 )-1)==filesep)) - 1; -CPath = P(1,1:d); -Q = P(:,d+1:size(P,2)); -%-Display data parameters -%======================================================================= -figure(Fgraph); spm_clf; axis off -text(0.30,1.02,'Statistical analysis','Fontsize',16,'Fontweight','Bold'); -text(-0.10,0.85,'Scan Index','Rotation',90) -if bMStud, text(-0.05,0.85,'Study', 'Rotation',90); end -if bMSubj, text(+0.00,0.85,'Subject', 'Rotation',90); end -if bMCond, text(+0.05,0.85,'Condition', 'Rotation',90); end -if bMRepl, text(+0.10,0.85,'Replication','Rotation',90); end -if bMXblk, text(+0.15,0.85,'Exchange Blk','Rotation',90); end -x0 = 0.20; y0 = 0.83; -dx = 0.15; dy = 0.02; -x = x0; -for i = 1:size(Cc,2) - text(x + 0.02,0.85,Ccnames(i,:),'Rotation',90); + %-Muck about a bit to set flags for various indicators - handy for later + bMStud=~isempty(iStud); + bMSubj=~isempty(iSubj); + bMCond=~isempty(iCond); + bMRepl=~isempty(iRepl); + bMXblk=~isempty(iXblk); + + %-Compute common path components - all paths will begin with file separator + %----------------------------------------------------------------------- + d = max(find(P(1,1:find(~all(P == ones(nScan,1)*P(1,:)), 1 )-1)==filesep)) - 1; + CPath = P(1,1:d); + Q = P(:,d+1:size(P,2)); + + %-Display data parameters + %======================================================================= + figure(Fgraph); spm_clf; axis off + text(0.30,1.02,'Statistical analysis','Fontsize',16,'Fontweight','Bold'); + text(-0.10,0.85,'Scan Index','Rotation',90) + if bMStud, text(-0.05,0.85,'Study', 'Rotation',90); end + if bMSubj, text(+0.00,0.85,'Subject', 'Rotation',90); end + if bMCond, text(+0.05,0.85,'Condition', 'Rotation',90); end + if bMRepl, text(+0.10,0.85,'Replication','Rotation',90); end + if bMXblk, text(+0.15,0.85,'Exchange Blk','Rotation',90); end + x0 = 0.20; y0 = 0.83; + dx = 0.15; dy = 0.02; + x = x0; + for i = 1:size(Cc,2) + text(x + 0.02,0.85,Ccnames(i,:),'Rotation',90); + x = x + dx; end + for i = 1:size(Gc,2) + text(x + 0.02,0.85,Gcnames(i,:),'Rotation',90); x = x + dx; end -for i = 1:size(Gc,2) - text(x + 0.02,0.85,Gcnames(i,:),'Rotation',90); -x = x + dx; end -text(x,0.92,'Base directory:','FontSize',10,'Fontweight','Bold'); -text(x,0.90,CPath,'FontSize',10,'interpreter','none'); -text(x,0.87,'Filename Tails'); -y = y0; - -for i = 1:nScan - text(-0.12,y,sprintf('%02d :',i)); - if bMStud, text(-0.06,y,sprintf('%2d',iStud(i))); end - if bMSubj, text(-0.01,y,sprintf('%2d',iSubj(i))); end - if bMCond, text(+0.04,y,sprintf('%2d',iCond(i))); end - if bMRepl, text(+0.09,y,sprintf('%2d',iRepl(i))); end - if bMXblk, text(+0.14,y,sprintf('%2d',iXblk(i))); end - x = x0; - for j = 1:size(Cc,2) - text(x,y,sprintf('%-8.6g',Cc(i,j)),'FontSize',10) - x = x + dx; end - for j = 1:size(Gc,2) - text(x,y,sprintf('%-8.6g',Gc(i,j)),'FontSize',10) - x = x + dx; end - text(x,y,Q(i,:),'FontSize',10,'interpreter','none'); - y = y - dy; - if y < 0; - spm_print - spm_clf; axis off - y = y0; - text(0.16,1.02,['Statistical analysis (continued)'],... - 'Fontsize',16,'Fontweight','Bold'); - end -end - -%-Print miscellaneous data parameters -%----------------------------------------------------------------------- -y = y - dy; -dy = dy*1.2; -if (GM~=0) - text(0,y,sprintf(['Images scaled to a grand mean of %g'],GM)) - y = y - dy; -end -text(0,y,sprintf(... - 'Analysis threshold is %3.0f%% of the whole brain mean',THRESH*100)) -spm_print - - -%-Display design parameters -%======================================================================= -figure(Fgraph); spm_clf(Fgraph); axis off -text(0.30,1.02,'Design Matrix','Fontsize',16,'Fontweight','Bold'); + text(x,0.92,'Base directory:','FontSize',10,'Fontweight','Bold'); + text(x,0.90,CPath,'FontSize',10,'interpreter','none'); + text(x,0.87,'Filename Tails'); + y = y0; + + for i = 1:nScan + text(-0.12,y,sprintf('%02d :',i)); + if bMStud, text(-0.06,y,sprintf('%2d',iStud(i))); end + if bMSubj, text(-0.01,y,sprintf('%2d',iSubj(i))); end + if bMCond, text(+0.04,y,sprintf('%2d',iCond(i))); end + if bMRepl, text(+0.09,y,sprintf('%2d',iRepl(i))); end + if bMXblk, text(+0.14,y,sprintf('%2d',iXblk(i))); end + x = x0; + for j = 1:size(Cc,2) + text(x,y,sprintf('%-8.6g',Cc(i,j)),'FontSize',10) + x = x + dx; end + for j = 1:size(Gc,2) + text(x,y,sprintf('%-8.6g',Gc(i,j)),'FontSize',10) + x = x + dx; end + text(x,y,Q(i,:),'FontSize',10,'interpreter','none'); + y = y - dy; + if y < 0; + spm_print + spm_clf; axis off + y = y0; + text(0.16,1.02,['Statistical analysis (continued)'],... + 'Fontsize',16,'Fontweight','Bold'); + end + end -%-Label the effects -%----------------------------------------------------------------------- -hDesMtx = axes('Position',[0.2 0.3 0.6 0.5]); -image((nHCBG + 1)*32); -ylabel('Observations') -set(hDesMtx,'XTick',[],'XTickLabel','') -hEfLabs = axes('Position',[0.2 0.82 0.6 0.1],'Visible','off'); -y = 0.1; -dx = 1/size(nHCBG,2); -for i = 1:size(nHCBG,2) - text((i - 0.5)*dx,y,deblank(HCBGnames(i,:)),... - 'Fontsize',8,'Rotation',90) -end + %-Print miscellaneous data parameters + %----------------------------------------------------------------------- + y = y - dy; + dy = dy*1.2; + if (GM~=0) + text(0,y,sprintf(['Images scaled to a grand mean of %g'],GM)) + y = y - dy; + end + text(0,y,sprintf(... + 'Analysis threshold is %3.0f%% of the whole brain mean',THRESH*100)) + spm_print + + + %-Display design parameters + %======================================================================= + figure(Fgraph); spm_clf(Fgraph); axis off + text(0.30,1.02,'Design Matrix','Fontsize',16,'Fontweight','Bold'); + + %-Label the effects + %----------------------------------------------------------------------- + hDesMtx = axes('Position',[0.2 0.3 0.6 0.5]); + image((nHCBG + 1)*32); + ylabel('Observations') + set(hDesMtx,'XTick',[],'XTickLabel','') + hEfLabs = axes('Position',[0.2 0.82 0.6 0.1],'Visible','off'); + y = 0.1; + dx = 1/size(nHCBG,2); + for i = 1:size(nHCBG,2) + text((i - 0.5)*dx,y,deblank(HCBGnames(i,:)),... + 'Fontsize',8,'Rotation',90) + end -%-Display non-parametric analysis summary -%----------------------------------------------------------------------- -hPramAxes=axes('Position',[0.05 0.08 0.8 0.20],'Visible','off'); -text(0,1.00,sDesign,'Fontsize',10); -text(0,0.90,['SnPM design flie: ',sDesFile],'Fontsize',10); -text(0,0.80,sPiCond,'Fontsize',10); -text(0,0.70,['Global normalisation: ',deblank(sGloNorm)],'Fontsize',10); -text(0,0.60,['Threshold masking: ',deblank(sThresh)],'Fontsize',10); - -%-Display parameter summary -%----------------------------------------------------------------------- -text(0,.5,'Parameters:','Fontsize',10,'Fontweight','Bold'); -text(0,.4,sprintf(['%d Condition + %d Covariate ',... - '+ %d Block + %d Confound'],... - size(H,2),size(C,2),size(B,2),size(G,2)),... - 'Fontsize',10); -text(0,.3,sprintf(['= %d parameters, having %d degrees of freedom, ',... - 'giving %d residual df (%d scans).'],... - size([H C B G],2),rank([H C B G]),nScan-rank([H C B G]),nScan),... - 'Fontsize',10); -if (bVarSm) text(0,0.2,sVarSm,'Fontsize',10); -end + %-Display non-parametric analysis summary + %----------------------------------------------------------------------- + hPramAxes=axes('Position',[0.05 0.08 0.8 0.20],'Visible','off'); + text(0,1.00,sDesign,'Fontsize',10); + text(0,0.90,['SnPM design flie: ',sDesFile],'Fontsize',10); + text(0,0.80,sPiCond,'Fontsize',10); + text(0,0.70,['Global normalisation: ',deblank(sGloNorm)],'Fontsize',10); + text(0,0.60,['Threshold masking: ',deblank(sThresh)],'Fontsize',10); + + %-Display parameter summary + %----------------------------------------------------------------------- + text(0,.5,'Parameters:','Fontsize',10,'Fontweight','Bold'); + text(0,.4,sprintf(['%d Condition + %d Covariate ',... + '+ %d Block + %d Confound'],... + size(H,2),size(C,2),size(B,2),size(G,2)),... + 'Fontsize',10); + text(0,.3,sprintf(['= %d parameters, having %d degrees of freedom, ',... + 'giving %d residual df (%d scans).'],... + size([H C B G],2),rank([H C B G]),nScan-rank([H C B G]),nScan),... + 'Fontsize',10); + if (bVarSm) text(0,0.2,sVarSm,'Fontsize',10); + end -spm_print + spm_print -%-Clear interactive window -%----------------------------------------------------------------------- -spm_clf(Finter) + %-Clear interactive window + %----------------------------------------------------------------------- + spm_clf(Finter) +end \ No newline at end of file diff --git a/spm_append_96.m b/spm_append_96.m index 218fb48..52d4ceb 100644 --- a/spm_append_96.m +++ b/spm_append_96.m @@ -47,7 +47,7 @@ function spm_append_96(MAT,X,SizeWarn) %---------------------------------------------------------------------------- [m n] = size(X); HDR = fread(fid,5,'int32'); -if m ~= HDR(2); error(' Incompatible sizes'); end +if m ~= HDR(2); error('SnPM:InvalidSizes', ' Incompatible sizes'); end fclose(fid); % append and update nunmber of columns @@ -57,7 +57,7 @@ function spm_append_96(MAT,X,SizeWarn) Len = ftell(fid); if Len>2^(10*3)*0.99 % 99% of 1GB if ~strcmp(LAST_MAT_SIZE_WARNING,MAT) - warning(['Very large file! (' MAT ')\n' SizeWarn]) + warning('SnPM:VeryLargeFile', ['Very large file! (' MAT ')\n' SizeWarn]) LAST_MAT_SIZE_WARNING = MAT; end end diff --git a/spm_xyz2e.m b/spm_xyz2e.m index 900f927..39424e4 100644 --- a/spm_xyz2e.m +++ b/spm_xyz2e.m @@ -27,7 +27,7 @@ %-Condition arguments %----------------------------------------------------------------------- -if (nargin<2), error('Insufficient arguments'), end +if (nargin<2), error('SnPM:InsufficientArguments', 'Insufficient arguments'), end ImDim = V(1).dim'; M = V(1).mat(1:3, 1:3); diff --git a/test/common/generic_test_snpm.m b/test/common/generic_test_snpm.m index 7803193..19e0797 100644 --- a/test/common/generic_test_snpm.m +++ b/test/common/generic_test_snpm.m @@ -35,8 +35,17 @@ methods (TestMethodSetup) function setGlobals(testCase) - global TEST; - TEST = true; + % Random number generator should not be initialised with a + % shuffled seed + global SnPMdefs + SnPMdefs.shuffle_seed = false; + + % Run the tests in command line mode (no display) + global defaults; + defaults.cmdline = true; + + % Disable warning on very small number of permutations + warning('off','SnPM:VeryFewPermsCoarseExactPValues') snpm_test_config; cd(spm_str_manip(which('snpm_test_config'), 'h')) @@ -45,7 +54,7 @@ function setGlobals(testCase) testCase.SnPMrefVersion = SnPMrefVersion; if isempty(testDataDir) - error('Test data directory not set, please update snpm_test_config'); + error('SnPM:NotTestDataDir', 'Test data directory not set, please update snpm_test_config'); end testCase.parentDataDir = spm_str_manip(testDataDir, 'h'); @@ -186,7 +195,12 @@ function assert_check(testCase) end testCase.compare_batch_with_inter(zeroingNaNs); - clear global TEST; + % Reinitialize SnPM & SPM defaults + snpm_defaults; + spm_defaults; + + % Reanable all warnings + warning('on','all'); end function complete_batch_and_run(testCase) @@ -245,9 +259,12 @@ function complete_batch_and_run(testCase) if ~exist(testCase.spmDir, 'dir') mkdir(testCase.spmDir); else - % Remove SPM.mat to avoid interactive window asking if - % model can be overwritten - delete(fullfile(testCase.spmDir, 'SPM.mat')); + spmmatfile = fullfile(testCase.spmDir, 'SPM.mat'); + if exist(spmmatfile, 'file') + % Remove SPM.mat to avoid interactive window asking if + % model can be overwritten + delete(spmmatfile); + end end % If grand mean scaling then we should calculate mean (otherwise error?) if ( isfield(testCase.spmBatch{1}.spm.stats.factorial_design, 'globalm') &&... @@ -405,7 +422,7 @@ function compare_batch_with_inter(testCase, zeroingNaNs) end if testCase.checks if numel(testCase.inter_map) ~= numel(testCase.batch_map) - error(['Number of ' testCase.mapName ' maps are not equal between batch (',... + error('SnPM:UnequalTestCases', ['Number of ' testCase.mapName ' maps are not equal between batch (',... num2str(numel(testCase.batch_map)),... ') and interactive mode ',... num2str(numel(testCase.inter_map)) ]); diff --git a/test/common/skeleton_class.m b/test/common/skeleton_class.m index 8c73b31..b434194 100644 --- a/test/common/skeleton_class.m +++ b/test/common/skeleton_class.m @@ -28,9 +28,14 @@ function test_testclassname_var(testCase) % With approximate test function test_testclassname_approx(testCase) - testCase.testName = 'testclassname_approx'; - - rand('seed',200); + testCase.testName = 'testclassname_approx'; + try + % Syntax for newest Matlab versions + rng(200); + catch + % Old syntax + rand('seed',200); + end end end diff --git a/test/ground_truth/snpm_test_ground_truth.m b/test/ground_truth/snpm_test_ground_truth.m index 473be83..708f5e4 100644 --- a/test/ground_truth/snpm_test_ground_truth.m +++ b/test/ground_truth/snpm_test_ground_truth.m @@ -12,7 +12,7 @@ function snpm_test_ground_truth() disp('Recomputing ground thruth for SnPM testing...') if isempty(testDataDir) - error('Test data directory not set, please update snpm_test_config'); + error('SnPM:NotTestDataDir', 'Test data directory not set, please update snpm_test_config'); end disp('Current version of SnPM') @@ -23,7 +23,7 @@ function snpm_test_ground_truth() switch buttonName, case 'no', - error('Re-computation of ground truth stopped!'); + error('SnPM:GroundTruthHalted', 'Re-computation of ground truth stopped!'); case 'yes', disp('Start re-computation of ground truth') diff --git a/test/test_ANOVAbetween.m b/test/test_ANOVAbetween.m index b3a7426..095ecb1 100644 --- a/test/test_ANOVAbetween.m +++ b/test/test_ANOVAbetween.m @@ -59,8 +59,14 @@ function test_ANOVAbetween_var(testCase) % With approximate test function test_ANOVAbetween_approx(testCase) - testCase.testName = 'ANOVAbetween_approx'; - rand('seed',200); + testCase.testName = 'ANOVAbetween_approx'; + try + % Syntax for newest Matlab versions + rng(200); + catch + % Old syntax + rand('seed',200); + end testCase.matlabbatch{1}.spm.tools.snpm.des.ANOVAbtwnS.nPerm = 15; end end diff --git a/test/test_multisub_withinsubanova.m b/test/test_multisub_withinsubanova.m index 8f904b6..3088ea4 100644 --- a/test/test_multisub_withinsubanova.m +++ b/test/test_multisub_withinsubanova.m @@ -23,7 +23,7 @@ function create_basis_matlabbatch(testCase) % Fill the design part in the batch for i = 1:nSubjects for j = 1:nScansPerSub - testCase.matlabbatch{1}.spm.tools.snpm.des.ANOVAwithinS.fsubject(i).scans{j} = ... + testCase.matlabbatch{1}.spm.tools.snpm.des.ANOVAwithinS.fsubject(i).scans{j,1} = ... fullfile(testCase.testDataDir, ['test_data_', num2str((i-1)*nScansPerSub+j, '%02.0f'), '.nii,1']); end end @@ -48,7 +48,13 @@ function test_multisub_withinsubanova_var(testCase) function test_multisub_withinsubanova_approx(testCase) testCase.testName = 'multisub_withinsubanova_approx'; - rand('seed',200); + try + % Syntax for newest Matlab versions + rng(200); + catch + % Old syntax + rand('seed',200); + end testCase.matlabbatch{1}.spm.tools.snpm.des.ANOVAwithinS.nPerm = 13; end end diff --git a/test/test_multisubpaired2cond.m b/test/test_multisubpaired2cond.m index 297d509..df79be4 100644 --- a/test/test_multisubpaired2cond.m +++ b/test/test_multisubpaired2cond.m @@ -50,7 +50,13 @@ function test_multisubpaired2cond_var(testCase) % With approximate test function test_multisubpaired2cond_approx(testCase) testCase.testName = 'multisubpaired2cond_approx'; - rand('seed',200); + try + % Syntax for newest Matlab versions + rng(200); + catch + % Old syntax + rand('seed',200); + end testCase.matlabbatch{1}.spm.tools.snpm.des.PairT.nPerm = 14; end diff --git a/test/test_multisubsimpleregression.m b/test/test_multisubsimpleregression.m index 1bfbdbe..57b2a51 100644 --- a/test/test_multisubsimpleregression.m +++ b/test/test_multisubsimpleregression.m @@ -15,7 +15,7 @@ function create_basis_matlabbatch(testCase) nSubjects = 4; for i = 1:nSubjects - testCase.matlabbatch{1}.spm.tools.snpm.des.Corr.P{i} = ... + testCase.matlabbatch{1}.spm.tools.snpm.des.Corr.P{i,1} = ... fullfile(testCase.testDataDir, ['test_data_', num2str(i, '%02.0f'), '.nii,1']); end testCase.matlabbatch{1}.spm.tools.snpm.des.Corr.CovInt = [1 3 5 0]; @@ -69,7 +69,13 @@ function test_multisubsimpleregression_var(testCase) function test_multisubsimpleregression_approx(testCase) testCase.testName = 'multisubsimpleregression_approx'; - rand('seed',200); + try + % Syntax for newest Matlab versions + rng(200); + catch + % Old syntax + rand('seed',200); + end testCase.matlabbatch{1}.spm.tools.snpm.des.Corr.nPerm = 14; end end diff --git a/test/test_oneSample.m b/test/test_oneSample.m index 1301503..1cb6482 100644 --- a/test/test_oneSample.m +++ b/test/test_oneSample.m @@ -198,7 +198,13 @@ function test_onesample_mask(testCase) function test_onesample_approx(testCase) testCase.testName = 'onesample_approx'; - rand('seed',200); + try + % Syntax for newest Matlab versions + rng(200); + catch + % Old syntax + rand('seed',200); + end testCase.matlabbatch{1}.spm.tools.snpm.des.OneSampT.nPerm = 15; end @@ -246,7 +252,11 @@ function test_onesample_slice(testCase) testCase.testName = 'onesample_slice'; - rand('seed',200); + try + rng(200); + catch + rand('seed',200); + end for i = 6:17 testCase.matlabbatch{1}.spm.tools.snpm.des.OneSampT.P{end+1} = ... fullfile(testCase.testDataDir, ['test_data_' num2str(i, '%02.0f') '.nii']); diff --git a/test/test_onesub_regression.m b/test/test_onesub_regression.m index c1c1b9c..a199ad8 100644 --- a/test/test_onesub_regression.m +++ b/test/test_onesub_regression.m @@ -70,7 +70,13 @@ function test_onesub_regression_var(testCase) function test_onesub_regression_approx(testCase) testCase.testName = 'onesub_regression_approx'; - rand('seed',200); + try + % Syntax for newest Matlab versions + rng(200); + catch + % Old syntax + rand('seed',200); + end testCase.matlabbatch{1}.spm.tools.snpm.des.Corr1S.nPerm = 25; end diff --git a/test/test_onesub_twocondrepl.m b/test/test_onesub_twocondrepl.m index e4c936a..69b21ab 100644 --- a/test/test_onesub_twocondrepl.m +++ b/test/test_onesub_twocondrepl.m @@ -15,7 +15,7 @@ function create_basis_matlabbatch(testCase) % Fill the design part in the batch for i = 1:12 - testCase.matlabbatch{1}.spm.tools.snpm.des.TwoSampTss.P{i} = ... + testCase.matlabbatch{1}.spm.tools.snpm.des.TwoSampTss.P{i,1} = ... fullfile(testCase.testDataDir, ['test_data_', num2str(i, '%02.0f') '.nii,1']); end testCase.matlabbatch{1}.spm.tools.snpm.des.TwoSampTss.condidx = [1 2 1 2 1 2 1 2 1 2 1 2]; @@ -58,7 +58,13 @@ function test_onesub_twocondrepl_approx(testCase) testCase.checks = false; testCase.testName = 'onesub_twocondrepl_approx'; - rand('seed',200); + try + % Syntax for newest Matlab versions + rng(200); + catch + % Old syntax + rand('seed',200); + end testCase.matlabbatch{1}.spm.tools.snpm.des.TwoSampTss.nPerm = 15; end end diff --git a/test/test_twoSample.m b/test/test_twoSample.m index 806d284..c0a44d4 100644 --- a/test/test_twoSample.m +++ b/test/test_twoSample.m @@ -122,7 +122,13 @@ function test_twosample_var(testCase) function test_twosample_approx(testCase) testCase.testName = 'twosample_approx'; - rand('seed',200); + try + % Syntax for newest Matlab versions + rng(200); + catch + % Old syntax + rand('seed',200); + end testCase.matlabbatch{1}.spm.tools.snpm.des.TwoSampT.nPerm = 15; end diff --git a/test/test_twosample_twocond.m b/test/test_twosample_twocond.m index 9c4a2d0..be4d46d 100644 --- a/test/test_twosample_twocond.m +++ b/test/test_twosample_twocond.m @@ -89,7 +89,13 @@ function test_twosample_twocond_var(testCase) % With approximate test function test_twosample_twocond_approx(testCase) testCase.testName = 'twosample_twocond_approx'; - rand('seed',200); + try + % Syntax for newest Matlab versions + rng(200); + catch + % Old syntax + rand('seed',200); + end testCase.matlabbatch{1}.spm.tools.snpm.des.TwoSampPairT.nPerm = 15; end end