From ec39fa12446c163456edd71faf436e2f1a197ee3 Mon Sep 17 00:00:00 2001
From: nicholst <t.e.nichols@warwick.ac.uk>
Date: Wed, 4 Feb 2015 15:02:29 +0100
Subject: [PATCH 01/22] Previously, implicit masking *was* given as an option,
 but in fact was not implimented for integer input data (i.e. data with no NaN
 representation).  Now it is properly implimented.

---
 config/snpm_bch_ui.m |  4 +++-
 snpm_cp.m            | 12 +++++++++---
 snpm_ui.m            |  6 +++---
 3 files changed, 15 insertions(+), 7 deletions(-)

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_cp.m b/snpm_cp.m
index 1a82799..2ca485d 100644
--- a/snpm_cp.m
+++ b/snpm_cp.m
@@ -247,12 +247,14 @@ function snpm_cp(CWD)
 if ~bVolm & pU_ST_Ut>=0
   error('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
@@ -497,8 +499,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);
diff --git a/snpm_ui.m b/snpm_ui.m
index 64b7f42..1cd644a 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
@@ -431,7 +431,7 @@ function snpm_ui(varargin)
 
 %-Implicit Masking - Batch only!
 %-----------------------------------------------------------------------
-ImMASKING=job.masking.im;
+ImMASK=job.masking.im;
 
 %-Get analysis mask
 %-----------------------------------------------------------------------
@@ -577,7 +577,7 @@ 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])

From d17cf584a5b6bedd3914455e7fb527515fd81254 Mon Sep 17 00:00:00 2001
From: nicholst <t.e.nichols@warwick.ac.uk>
Date: Thu, 5 Feb 2015 02:02:26 +0100
Subject: [PATCH 02/22] Added check for NaN from spm_global, making it easier
 to track down pathological problems (e.g. from all and highly negative
 images).

---
 snpm_ui.m | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/snpm_ui.m b/snpm_ui.m
index 1cd644a..f66064f 100644
--- a/snpm_ui.m
+++ b/snpm_ui.m
@@ -479,6 +479,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('Global computation returned NaN! Cannot continue')
+  end
   GX     = rg;
 elseif iGXcalc==1
   rg     = [];

From 46e196ab3189a28b30324760fd1bb4fd00e8f10a Mon Sep 17 00:00:00 2001
From: nicholst <t.e.nichols@warwick.ac.uk>
Date: Sat, 7 Feb 2015 10:09:29 +0000
Subject: [PATCH 03/22] Fixed Matlab R2014b compatibility problem; XlabelTicks,
 formerly a string array, is now a cell array. This broke the FWE report
 *only* for cluster size inference (as the X-ticks are cubed).  Found a fix
 that is compatible with both R2014b and R2013a (the 2 versions I have).

---
 snpm_pp.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/snpm_pp.m b/snpm_pp.m
index 93e42ef..3e603ae 100644
--- a/snpm_pp.m
+++ b/snpm_pp.m
@@ -1619,7 +1619,7 @@ function ShowDist(T,cT,aT,C,cC,aC,Typ)
     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','-');

From b9a1fdbd67d734346854a84aef9077ce8e3eff34 Mon Sep 17 00:00:00 2001
From: nicholst <t.e.nichols@warwick.ac.uk>
Date: Tue, 17 May 2016 17:00:11 +0100
Subject: [PATCH 04/22] Allows MC permutation (with reptitions) for large
 enough nPerms

---
 snpm_pi_TwoSampT.m | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/snpm_pi_TwoSampT.m b/snpm_pi_TwoSampT.m
index b82861a..e14615b 100644
--- a/snpm_pi_TwoSampT.m
+++ b/snpm_pi_TwoSampT.m
@@ -115,6 +115,9 @@
 if isempty(TEST) || ~TEST % When testing the code we need a fixed seed
     rand('seed',sum(100*clock));	% Initialise random number generator
 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 +221,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      

From ba6a2026f3414043083c1732496e74bf8f86ada7 Mon Sep 17 00:00:00 2001
From: cmaumet <c.m.j.maumet@warwick.ac.uk>
Date: Wed, 20 Jul 2016 10:23:55 +0100
Subject: [PATCH 05/22] New syntax for random number generator seed

rng(x) instead of rand('seed', x)
---
 test/common/skeleton_class.m         | 11 ++++++++---
 test/test_ANOVAbetween.m             | 10 ++++++++--
 test/test_multisub_withinsubanova.m  |  8 +++++++-
 test/test_multisubpaired2cond.m      |  8 +++++++-
 test/test_multisubsimpleregression.m |  8 +++++++-
 test/test_oneSample.m                |  8 +++++++-
 test/test_onesub_regression.m        |  8 +++++++-
 test/test_onesub_twocondrepl.m       |  8 +++++++-
 test/test_twoSample.m                |  8 +++++++-
 test/test_twosample_twocond.m        |  8 +++++++-
 10 files changed, 72 insertions(+), 13 deletions(-)

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/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..069761f 100644
--- a/test/test_multisub_withinsubanova.m
+++ b/test/test_multisub_withinsubanova.m
@@ -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..e53d055 100644
--- a/test/test_multisubsimpleregression.m
+++ b/test/test_multisubsimpleregression.m
@@ -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..c300b34 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
         
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..1130c91 100644
--- a/test/test_onesub_twocondrepl.m
+++ b/test/test_onesub_twocondrepl.m
@@ -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

From 4058c9b0cd9b8abccff1ce67d35859821e9c14f8 Mon Sep 17 00:00:00 2001
From: cmaumet <c.m.j.maumet@warwick.ac.uk>
Date: Wed, 20 Jul 2016 10:32:26 +0100
Subject: [PATCH 06/22] New SnPM parameter to disable seed shuffling

fixes #3
---
 snpm_defaults.m                 |  7 +++++++
 snpm_pi_ANOVAbtwnS.m            | 13 ++++++++++---
 snpm_pi_ANOVAwithinS.m          | 13 ++++++++++---
 snpm_pi_Corr.m                  | 17 ++++++++++-------
 snpm_pi_Corr1S.m                | 18 ++++++++++--------
 snpm_pi_OneSampT.m              | 12 +++++++++---
 snpm_pi_TwoSampPairT.m          | 12 +++++++++---
 snpm_pi_TwoSampT.m              | 12 +++++++++---
 snpm_pi_TwoSampTss.m            | 16 +++++++++++-----
 test/common/generic_test_snpm.m |  9 ++++++---
 10 files changed, 91 insertions(+), 38 deletions(-)

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_pi_ANOVAbtwnS.m b/snpm_pi_ANOVAbtwnS.m
index cc53e23..9fd9458 100644
--- a/snpm_pi_ANOVAbtwnS.m
+++ b/snpm_pi_ANOVAbtwnS.m
@@ -102,9 +102,16 @@
 %%% 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
+global SnPMdefs;
+
+if SnPMdefs.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
diff --git a/snpm_pi_ANOVAwithinS.m b/snpm_pi_ANOVAwithinS.m
index 642dd44..4e1c2c5 100644
--- a/snpm_pi_ANOVAwithinS.m
+++ b/snpm_pi_ANOVAwithinS.m
@@ -108,11 +108,18 @@
 
 %-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
+
+global SnPMdefs
+if SnPMdefs.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
diff --git a/snpm_pi_Corr.m b/snpm_pi_Corr.m
index e11bc7c..0eaab5b 100644
--- a/snpm_pi_Corr.m
+++ b/snpm_pi_Corr.m
@@ -74,7 +74,16 @@
 iGloNorm = '123';		% Allowable Global norm. codes
 sDesSave = 'CovInt';		% PlugIn variables to save in cfg file
 
-global TEST;
+global SnPMdefs
+if SnPMdefs.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
 %-----------------------------------------------------------------------
@@ -181,9 +190,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
@@ -239,9 +245,6 @@
 %-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])
diff --git a/snpm_pi_Corr1S.m b/snpm_pi_Corr1S.m
index 7a22466..a3da34d 100644
--- a/snpm_pi_Corr1S.m
+++ b/snpm_pi_Corr1S.m
@@ -79,11 +79,19 @@
 
 %-Initialisation
 %-----------------------------------------------------------------------
-global TEST;
 iGloNorm = '123';		% Allowable Global norm. codes
 sDesSave = 'CovInt Xblk';	% PlugIn variables to save in cfg file
 
-
+global SnPMdefs
+if SnPMdefs.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
 %-----------------------------------------------------------------------
@@ -176,9 +184,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
@@ -246,9 +251,6 @@
 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])
diff --git a/snpm_pi_OneSampT.m b/snpm_pi_OneSampT.m
index 2359b2f..13d2540 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
+global SnPMdefs
+if SnPMdefs.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
diff --git a/snpm_pi_TwoSampPairT.m b/snpm_pi_TwoSampPairT.m
index 929f3d8..e88f08c 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
+global SnPMdefs
+if SnPMdefs.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
 
diff --git a/snpm_pi_TwoSampT.m b/snpm_pi_TwoSampT.m
index b82861a..206ba9e 100644
--- a/snpm_pi_TwoSampT.m
+++ b/snpm_pi_TwoSampT.m
@@ -111,9 +111,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 the code we need a fixed seed
-    rand('seed',sum(100*clock));	% Initialise random number generator
+global SnPMdefs
+if SnPMdefs.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
diff --git a/snpm_pi_TwoSampTss.m b/snpm_pi_TwoSampTss.m
index 504b1f9..fa3e82e 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,17 @@
 %-Get number of replications per condition - 2 x nRepl design
 nRepl    = job.Tss_repc;%spm_input('# replications per condition','+1');
 
+global SnPMdefs
+if SnPMdefs.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
@@ -178,10 +188,6 @@
 	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/test/common/generic_test_snpm.m b/test/common/generic_test_snpm.m
index 7803193..ed5adaa 100644
--- a/test/common/generic_test_snpm.m
+++ b/test/common/generic_test_snpm.m
@@ -35,8 +35,10 @@
     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;
 
             snpm_test_config;
             cd(spm_str_manip(which('snpm_test_config'), 'h'))
@@ -186,7 +188,8 @@ function assert_check(testCase)
             end
             testCase.compare_batch_with_inter(zeroingNaNs); 
             
-            clear global TEST;
+            % Reinitialize SnPM defaults
+            snpm_defaults;
         end
         
         function complete_batch_and_run(testCase)

From 16255d3ebf1e299d984825a6e0c8b7a2c82cacd4 Mon Sep 17 00:00:00 2001
From: cmaumet <c.m.j.maumet@warwick.ac.uk>
Date: Wed, 20 Jul 2016 11:21:26 +0100
Subject: [PATCH 07/22] Do not show figures in command-line mode

Re-use SPM's defaults.cmdline
Fixes #6
---
 snpm_ui.m | 262 +++++++++++++++++++++++++++---------------------------
 1 file changed, 132 insertions(+), 130 deletions(-)

diff --git a/snpm_ui.m b/snpm_ui.m
index f66064f..15ded20 100644
--- a/snpm_ui.m
+++ b/snpm_ui.m
@@ -588,141 +588,143 @@ function snpm_ui(varargin)
 
 
 
+if ~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
-%=======================================================================
-
-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
-
-
-%-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);
+    %-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

From 82a9469a0d2b85d13b6daf0bef3a90bb4f16b661 Mon Sep 17 00:00:00 2001
From: cmaumet <c.m.j.maumet@warwick.ac.uk>
Date: Wed, 20 Jul 2016 11:54:25 +0100
Subject: [PATCH 08/22] Add identifiers to errors and warnings

---
 snpm.m                                     |  2 +-
 snpm_check_nperm.m                         |  6 +++-
 snpm_combo.m                               |  2 +-
 snpm_combo_pp.m                            | 18 ++++++------
 snpm_cp.m                                  | 19 +++++++------
 snpm_init.m                                |  2 +-
 snpm_pi_ANOVAbtwnS.m                       |  6 ++--
 snpm_pi_ANOVAwithinS.m                     |  6 ++--
 snpm_pi_Corr.m                             |  8 +++---
 snpm_pi_Corr1S.m                           |  6 ++--
 snpm_pi_OneSampT.m                         |  4 +--
 snpm_pi_PairT.m                            | 16 +++++------
 snpm_pi_PairTrand.m                        | 17 +++++------
 snpm_pi_TwoSampPairT.m                     |  4 +--
 snpm_pi_TwoSampT.m                         |  4 +--
 snpm_pi_TwoSampTss.m                       | 10 +++----
 snpm_pp.m                                  | 33 +++++++++++-----------
 snpm_t2z.m                                 |  6 ++--
 snpm_uc_FDR.m                              |  2 +-
 snpm_ui.m                                  | 18 ++++++------
 spm_append_96.m                            |  4 +--
 spm_xyz2e.m                                |  2 +-
 test/common/generic_test_snpm.m            |  4 +--
 test/ground_truth/snpm_test_ground_truth.m |  4 +--
 24 files changed, 107 insertions(+), 96 deletions(-)

diff --git a/snpm.m b/snpm.m
index 1a9477b..6877320 100644
--- a/snpm.m
+++ b/snpm.m
@@ -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<C_MaxT)
-                        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  
         end
@@ -216,14 +218,14 @@ function snpm_combo_pp(CWD, job)
         pU_ST_Ut  = 1-spm_Tcdf(ST_Ut,df);
         if alpha==1	% Not filtering on significance
             if ~( 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 && primaryThresh<C_MaxT) || ...
                     (primaryThresh>pU_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 2ca485d..6d20161 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,16 +237,16 @@ 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; note if NaN's available
 NaNrep=0;
@@ -695,7 +696,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
 
@@ -979,7 +980,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_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 9fd9458..7473bf0 100644
--- a/snpm_pi_ANOVAbtwnS.m
+++ b/snpm_pi_ANOVAbtwnS.m
@@ -136,7 +136,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/...';
 
@@ -266,7 +266,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; 
 %-----------------------------------------------------------------------
@@ -294,7 +294,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 4e1c2c5..7b696ef 100644
--- a/snpm_pi_ANOVAwithinS.m
+++ b/snpm_pi_ANOVAwithinS.m
@@ -124,12 +124,12 @@
 
 %-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
 
 
@@ -274,7 +274,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 0eaab5b..2a81cd9 100644
--- a/snpm_pi_Corr.m
+++ b/snpm_pi_Corr.m
@@ -99,7 +99,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
@@ -126,7 +126,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
@@ -241,14 +241,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
 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 a3da34d..2ec3811 100644
--- a/snpm_pi_Corr1S.m
+++ b/snpm_pi_Corr1S.m
@@ -122,7 +122,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
@@ -246,7 +246,7 @@
 %-----------------------------------------------------------------------
 %-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);
@@ -254,7 +254,7 @@
 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 13d2540..e5cbd13 100644
--- a/snpm_pi_OneSampT.m
+++ b/snpm_pi_OneSampT.m
@@ -124,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
@@ -250,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 e88f08c..eb3f355 100644
--- a/snpm_pi_TwoSampPairT.m
+++ b/snpm_pi_TwoSampPairT.m
@@ -283,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
 %-----------------------------------------------------------------------
@@ -309,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 4fc4f8f..c8210a2 100644
--- a/snpm_pi_TwoSampT.m
+++ b/snpm_pi_TwoSampT.m
@@ -240,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
 %-----------------------------------------------------------------------
@@ -266,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 fa3e82e..0b826fb 100644
--- a/snpm_pi_TwoSampTss.m
+++ b/snpm_pi_TwoSampTss.m
@@ -139,7 +139,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.
 %-----------------------------------------------------------------------
@@ -166,10 +166,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
 
@@ -182,7 +182,7 @@
 %-----------------------------------------------------------------------
 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,:)=[];
diff --git a/snpm_pp.m b/snpm_pp.m
index 3e603ae..b76661b 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...');
@@ -1525,7 +1526,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,7 +1615,7 @@ 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
@@ -1643,7 +1644,7 @@ function ShowDist(T,cT,aT,C,cC,aC,Typ)
 
 function v = BoundCheck(val,Range,ErrMsg)
 if val<Range(1) | val>Range(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 15ded20..f8d3c06 100644
--- a/snpm_ui.m
+++ b/snpm_ui.m
@@ -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')
@@ -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
@@ -481,7 +483,7 @@ function snpm_ui(varargin)
   for i  = 1:nScan, rg(i)=spm_global(V(i)); end
   if any(~isfinite(rg))
     disp(rg)
-    error('Global computation returned NaN! Cannot continue')
+    error('SnPM:NaNGlobal', 'Global computation returned NaN! Cannot continue')
   end
   GX     = rg;
 elseif iGXcalc==1
@@ -512,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
 
 
@@ -554,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)
 
 
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 ed5adaa..ffb9e40 100644
--- a/test/common/generic_test_snpm.m
+++ b/test/common/generic_test_snpm.m
@@ -47,7 +47,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');
@@ -408,7 +408,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/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')
       

From 55d22b6038d9027f96cb73689d78710a088e9641 Mon Sep 17 00:00:00 2001
From: cmaumet <c.m.j.maumet@warwick.ac.uk>
Date: Wed, 20 Jul 2016 14:15:25 +0100
Subject: [PATCH 09/22] Statistic value from p-value depend on stat type

fixes #14
---
 snpm_cp.m | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/snpm_cp.m b/snpm_cp.m
index 6d20161..cc76b21 100644
--- a/snpm_cp.m
+++ b/snpm_cp.m
@@ -734,7 +734,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  

From e279f6fd725dae88e920ff405055252584f8c934 Mon Sep 17 00:00:00 2001
From: cmaumet <c.m.j.maumet@warwick.ac.uk>
Date: Wed, 20 Jul 2016 14:21:53 +0100
Subject: [PATCH 10/22] Add how to run the tests in README

---
 README.md | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/README.md b/README.md
index fd5368d..b9a1391 100644
--- a/README.md
+++ b/README.md
@@ -3,3 +3,11 @@
 The <b>S</b>tatistical <b>n</b>on-<b>P</b>arametric <b>M</b>apping (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.
+
+##### Testing
+To run the test suite you will first have to create a set of ground truth data. Then, the tests can be started with:
+```
+import matlab.unittest.TestSuite;
+suite = TestSuite.fromFolder(fullfile(spm_str_manip(which('snpm'), 'h'), 'test'));
+result = run(suite);
+```
\ No newline at end of file

From 300303e8f86015f9a9a33a127f287fc8e840c57a Mon Sep 17 00:00:00 2001
From: cmaumet <c.m.j.maumet@warwick.ac.uk>
Date: Wed, 20 Jul 2016 14:30:26 +0100
Subject: [PATCH 11/22] First set up of the tests

---
 README.md | 15 +++++++++++++--
 1 file changed, 13 insertions(+), 2 deletions(-)

diff --git a/README.md b/README.md
index b9a1391..0b6eeda 100644
--- a/README.md
+++ b/README.md
@@ -4,8 +4,19 @@ The <b>S</b>tatistical <b>n</b>on-<b>P</b>arametric <b>M</b>apping (SnPM) toolbo
 
 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.
 
-##### Testing
-To run the test suite you will first have to create a set of ground truth data. Then, the tests can be started with:
+#### 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 tests
+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'));

From caf8f7aeeb167c85e099be0565398d895ac9cd38 Mon Sep 17 00:00:00 2001
From: cmaumet <c.m.j.maumet@warwick.ac.uk>
Date: Wed, 20 Jul 2016 14:36:27 +0100
Subject: [PATCH 12/22] Example to run a  single test

---
 README.md | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/README.md b/README.md
index 0b6eeda..0febebf 100644
--- a/README.md
+++ b/README.md
@@ -15,10 +15,15 @@ global SnPMrefVersion;
 SnPMrefVersion = 'SnPM8';
 ```
 
-##### Run the tests
+##### 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')
 ```
\ No newline at end of file

From b7d899a272e904cd88b6cdd791b4376c1c621358 Mon Sep 17 00:00:00 2001
From: cmaumet <c.m.j.maumet@warwick.ac.uk>
Date: Wed, 20 Jul 2016 14:38:35 +0100
Subject: [PATCH 13/22] Test in command line mode

---
 test/common/generic_test_snpm.m | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/test/common/generic_test_snpm.m b/test/common/generic_test_snpm.m
index ffb9e40..84bb36f 100644
--- a/test/common/generic_test_snpm.m
+++ b/test/common/generic_test_snpm.m
@@ -39,6 +39,10 @@ function setGlobals(testCase)
             % shuffled seed
             global SnPMdefs
             SnPMdefs.shuffle_seed = false;
+            
+            % Run the tests in command line mode (no display)
+            global defaults;
+            defaults.cmdline = true;
 
             snpm_test_config;
             cd(spm_str_manip(which('snpm_test_config'), 'h'))
@@ -188,8 +192,9 @@ function assert_check(testCase)
             end
             testCase.compare_batch_with_inter(zeroingNaNs); 
             
-            % Reinitialize SnPM defaults
+            % Reinitialize SnPM & SPM defaults
             snpm_defaults;
+            spm_defaults;
         end
         
         function complete_batch_and_run(testCase)

From 6aebb79a47c98f52a17479427a71f488e511a644 Mon Sep 17 00:00:00 2001
From: cmaumet <c.m.j.maumet@warwick.ac.uk>
Date: Wed, 20 Jul 2016 14:49:06 +0100
Subject: [PATCH 14/22] Get snpm default using snpm_get_defaults function

(rather than direct call to global variable)
---
 snpm_cp.m              | 6 ++----
 snpm_pi_ANOVAbtwnS.m   | 3 +--
 snpm_pi_ANOVAwithinS.m | 3 +--
 snpm_pi_Corr.m         | 3 +--
 snpm_pi_Corr1S.m       | 3 +--
 snpm_pi_OneSampT.m     | 4 ++--
 snpm_pi_TwoSampPairT.m | 4 ++--
 snpm_pi_TwoSampT.m     | 4 ++--
 snpm_pi_TwoSampTss.m   | 3 +--
 snpm_ui.m              | 2 +-
 10 files changed, 14 insertions(+), 21 deletions(-)

diff --git a/snpm_cp.m b/snpm_cp.m
index cc76b21..e322346 100644
--- a/snpm_cp.m
+++ b/snpm_cp.m
@@ -279,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
 
diff --git a/snpm_pi_ANOVAbtwnS.m b/snpm_pi_ANOVAbtwnS.m
index 7473bf0..8d890f9 100644
--- a/snpm_pi_ANOVAbtwnS.m
+++ b/snpm_pi_ANOVAbtwnS.m
@@ -102,9 +102,8 @@
 %%% nCond    = 2;			% Number of conditions (groups)
 iGloNorm = '123';		% Allowable Global norm. codes
 sDesSave = 'iCond GrpCnt';	% PlugIn variables to save in cfg file
-global SnPMdefs;
 
-if SnPMdefs.shuffle_seed
+if snpm_get_defaults('shuffle_seed')
     % Shuffle seed of random number generator
     try
         rng('shuffle');
diff --git a/snpm_pi_ANOVAwithinS.m b/snpm_pi_ANOVAwithinS.m
index 7b696ef..5b7d8a0 100644
--- a/snpm_pi_ANOVAwithinS.m
+++ b/snpm_pi_ANOVAwithinS.m
@@ -111,8 +111,7 @@
 iGloNorm = '123';		% Allowable Global norm. codes
 sDesSave = 'iRepl sHCform_Mtx';		        % PlugIn variables to save in cfg file
 
-global SnPMdefs
-if SnPMdefs.shuffle_seed
+if snpm_get_defaults('shuffle_seed')
     % Shuffle seed of random number generator
     try
         rng('shuffle');
diff --git a/snpm_pi_Corr.m b/snpm_pi_Corr.m
index 2a81cd9..d7d5499 100644
--- a/snpm_pi_Corr.m
+++ b/snpm_pi_Corr.m
@@ -74,8 +74,7 @@
 iGloNorm = '123';		% Allowable Global norm. codes
 sDesSave = 'CovInt';		% PlugIn variables to save in cfg file
 
-global SnPMdefs
-if SnPMdefs.shuffle_seed
+if snpm_get_defaults('shuffle_seed')
     % Shuffle seed of random number generator
     try
         rng('shuffle');
diff --git a/snpm_pi_Corr1S.m b/snpm_pi_Corr1S.m
index 2ec3811..e527e97 100644
--- a/snpm_pi_Corr1S.m
+++ b/snpm_pi_Corr1S.m
@@ -82,8 +82,7 @@
 iGloNorm = '123';		% Allowable Global norm. codes
 sDesSave = 'CovInt Xblk';	% PlugIn variables to save in cfg file
 
-global SnPMdefs
-if SnPMdefs.shuffle_seed
+if snpm_get_defaults('shuffle_seed')
     % Shuffle seed of random number generator
     try
         rng('shuffle');
diff --git a/snpm_pi_OneSampT.m b/snpm_pi_OneSampT.m
index e5cbd13..53db09d 100644
--- a/snpm_pi_OneSampT.m
+++ b/snpm_pi_OneSampT.m
@@ -94,8 +94,8 @@
 %-----------------------------------------------------------------------
 iGloNorm = '123';		% Allowable Global norm. codes
 sDesSave = 'iCond';		% PlugIn variables to save in cfg file
-global SnPMdefs
-if SnPMdefs.shuffle_seed
+
+if snpm_get_defaults('shuffle_seed')
     % Shuffle seed of random number generator
     try
         rng('shuffle');
diff --git a/snpm_pi_TwoSampPairT.m b/snpm_pi_TwoSampPairT.m
index eb3f355..77063e4 100644
--- a/snpm_pi_TwoSampPairT.m
+++ b/snpm_pi_TwoSampPairT.m
@@ -110,8 +110,8 @@
 nStud    = 2;			% Number of groups
 iGloNorm = '123';		% Allowable Global norm. codes
 sDesSave = 'iStud iCond nSubj'; % PlugIn variables to save in cfg file
-global SnPMdefs
-if SnPMdefs.shuffle_seed
+
+if snpm_get_defaults('shuffle_seed')
     % Shuffle seed of random number generator
     try
         rng('shuffle');
diff --git a/snpm_pi_TwoSampT.m b/snpm_pi_TwoSampT.m
index c8210a2..086d685 100644
--- a/snpm_pi_TwoSampT.m
+++ b/snpm_pi_TwoSampT.m
@@ -111,8 +111,8 @@
 nCond    = 2;			% Number of conditions (groups)
 iGloNorm = '123';		% Allowable Global norm. codes
 sDesSave = 'iCond GrpCnt';	% PlugIn variables to save in cfg file
-global SnPMdefs
-if SnPMdefs.shuffle_seed
+
+if snpm_get_defaults('shuffle_seed')
     % Shuffle seed of random number generator
     try
         rng('shuffle');
diff --git a/snpm_pi_TwoSampTss.m b/snpm_pi_TwoSampTss.m
index 0b826fb..771f557 100644
--- a/snpm_pi_TwoSampTss.m
+++ b/snpm_pi_TwoSampTss.m
@@ -85,8 +85,7 @@
 %-Get number of replications per condition - 2 x nRepl design
 nRepl    = job.Tss_repc;%spm_input('# replications per condition','+1');
 
-global SnPMdefs
-if SnPMdefs.shuffle_seed
+if snpm_get_defaults('shuffle_seed')
     % Shuffle seed of random number generator
     try
         rng('shuffle');
diff --git a/snpm_ui.m b/snpm_ui.m
index f8d3c06..d3aa45d 100644
--- a/snpm_ui.m
+++ b/snpm_ui.m
@@ -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

From a99a5033ca4f51a800dc02df02d9ebac28878a3d Mon Sep 17 00:00:00 2001
From: cmaumet <c.m.j.maumet@warwick.ac.uk>
Date: Wed, 20 Jul 2016 14:54:49 +0100
Subject: [PATCH 15/22] Disable warning on few permutations for testing

---
 test/common/generic_test_snpm.m | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/test/common/generic_test_snpm.m b/test/common/generic_test_snpm.m
index 84bb36f..51222c4 100644
--- a/test/common/generic_test_snpm.m
+++ b/test/common/generic_test_snpm.m
@@ -43,6 +43,9 @@ function setGlobals(testCase)
             % 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'))
@@ -195,6 +198,9 @@ function assert_check(testCase)
             % Reinitialize SnPM & SPM defaults
             snpm_defaults;
             spm_defaults;
+            
+            % Reanable all warnings
+            warning('on','all');
         end
         
         function complete_batch_and_run(testCase)

From 98f0902ac6838f1472f32644af40b19c7c43bd00 Mon Sep 17 00:00:00 2001
From: cmaumet <c.m.j.maumet@warwick.ac.uk>
Date: Wed, 20 Jul 2016 15:04:09 +0100
Subject: [PATCH 16/22] Even less display in command-line mode

---
 snpm_pp.m | 7 +++++--
 snpm_ui.m | 2 +-
 2 files changed, 6 insertions(+), 3 deletions(-)

diff --git a/snpm_pp.m b/snpm_pp.m
index b76661b..5f572c9 100644
--- a/snpm_pp.m
+++ b/snpm_pp.m
@@ -1141,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
 %=======================================================================
@@ -1439,6 +1440,8 @@ function snpm_pp(CWD,varargin)
 
 end
 
+end
+
 %- Image output?
 %=======================================================================
 %-Write out filtered SnPMt?
diff --git a/snpm_ui.m b/snpm_ui.m
index d3aa45d..221a485 100644
--- a/snpm_ui.m
+++ b/snpm_ui.m
@@ -590,7 +590,7 @@ function snpm_ui(varargin)
 
 
 
-if ~defaults.cmdline
+if ~spm_get_defaults('cmdline')
 
     %=======================================================================
     %-Display parameters

From 6aa9511f544468a6eaa029433c18a4e7cf0e4b7b Mon Sep 17 00:00:00 2001
From: cmaumet <c.m.j.maumet@warwick.ac.uk>
Date: Wed, 20 Jul 2016 15:08:48 +0100
Subject: [PATCH 17/22] Delete previous SPM.mat only if it exists

---
 test/common/generic_test_snpm.m | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/test/common/generic_test_snpm.m b/test/common/generic_test_snpm.m
index 51222c4..19e0797 100644
--- a/test/common/generic_test_snpm.m
+++ b/test/common/generic_test_snpm.m
@@ -259,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') &&...

From fe57ff3bc6c005a5ecb7584050b27ffa1855b7f2 Mon Sep 17 00:00:00 2001
From: cmaumet <c.m.j.maumet@warwick.ac.uk>
Date: Wed, 20 Jul 2016 15:10:03 +0100
Subject: [PATCH 18/22] Add missing rng

---
 test/test_oneSample.m | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/test/test_oneSample.m b/test/test_oneSample.m
index c300b34..1cb6482 100644
--- a/test/test_oneSample.m
+++ b/test/test_oneSample.m
@@ -252,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']);

From f35e58a3de97085a4627a347439fa5e6056a89ab Mon Sep 17 00:00:00 2001
From: cmaumet <c.m.j.maumet@warwick.ac.uk>
Date: Wed, 20 Jul 2016 16:05:19 +0100
Subject: [PATCH 19/22] One scan per row

---
 test/test_multisub_withinsubanova.m  | 2 +-
 test/test_multisubsimpleregression.m | 2 +-
 test/test_onesub_twocondrepl.m       | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/test/test_multisub_withinsubanova.m b/test/test_multisub_withinsubanova.m
index 069761f..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
diff --git a/test/test_multisubsimpleregression.m b/test/test_multisubsimpleregression.m
index e53d055..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];
diff --git a/test/test_onesub_twocondrepl.m b/test/test_onesub_twocondrepl.m
index 1130c91..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];

From c20946cb31e6f63e47398af50761a69d0acf2106 Mon Sep 17 00:00:00 2001
From: cmaumet <c.m.j.maumet@warwick.ac.uk>
Date: Thu, 21 Jul 2016 16:24:06 +0100
Subject: [PATCH 20/22] Version 13.1.04

---
 snpm.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/snpm.m b/snpm.m
index 6877320..d321c7d 100644
--- a/snpm.m
+++ b/snpm.m
@@ -49,7 +49,7 @@
 
 %-Parameters
 %-----------------------------------------------------------------------
-SnPMver    = 'SnPM13.1.03';
+SnPMver    = 'SnPM13.1.04';
 
 %-Format arguments
 %-----------------------------------------------------------------------

From 20f7f37eb8c0e33528036f1155cecc5d2d35a1cf Mon Sep 17 00:00:00 2001
From: cmaumet <c.m.j.maumet@warwick.ac.uk>
Date: Fri, 22 Jul 2016 08:39:29 +0100
Subject: [PATCH 21/22] Add release notes to README

---
 README.md | 164 +++++++++++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 163 insertions(+), 1 deletion(-)

diff --git a/README.md b/README.md
index 0febebf..dca7e15 100644
--- a/README.md
+++ b/README.md
@@ -26,4 +26,166 @@ result = run(suite);
 ##### Run a single test
 ```
 run(test_oneSample, 'test_onesample_1')
-```
\ No newline at end of file
+```
+
+#### 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.

From 6e2c4e3a83a6571dac8b9eaaf71821412146cf5b Mon Sep 17 00:00:00 2001
From: cmaumet <c.m.j.maumet@warwick.ac.uk>
Date: Fri, 22 Jul 2016 08:49:00 +0100
Subject: [PATCH 22/22] Links to main documentation page

---
 README.md | 60 +++++++++++++++++++++++++++++++------------------------
 1 file changed, 34 insertions(+), 26 deletions(-)

diff --git a/README.md b/README.md
index dca7e15..ee1fdde 100644
--- a/README.md
+++ b/README.md
@@ -1,12 +1,20 @@
-### SnPM: Statistical nonParametric Mapping
+## SnPM: Statistical nonParametric Mapping
 
-The <b>S</b>tatistical <b>n</b>on-<b>P</b>arametric <b>M</b>apping (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.
 
-#### Testing
+**More information at: www.warwick.ac.uk/snpm**
 
-##### Initial set up 
+ - [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;
@@ -15,7 +23,7 @@ global SnPMrefVersion;
 SnPMrefVersion = 'SnPM8';
 ```
 
-##### Run the test suite
+#### Run the test suite
 Then, the tests can be started (from the test data folder) with:
 ```
 import matlab.unittest.TestSuite;
@@ -23,17 +31,17 @@ suite = TestSuite.fromFolder(fullfile(spm_str_manip(which('snpm'), 'h'), 'test')
 result = run(suite);
 ```
 
-##### Run a single test
+#### Run a single test
 ```
 run(test_oneSample, 'test_onesample_1')
 ```
 
-#### SnPM13: Bugs & Fixes
+### 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 13.1.04
 * snpm
 SnPM version 13.1.04
 
@@ -69,14 +77,14 @@ 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
+####  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
+####  Updates from SnPM 13.1.02
 * snpm
 SnPM version 13.1.02
 
@@ -86,30 +94,30 @@ Update contrast display for compatibility with Matlab R2014b.
 * generic_test_snpm
 Updates affecting tests only.
 
-#####  Updates from SnPM 13.1.01
+####  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
+####  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
+####  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
+####  Updates from SnPM13.0.13
 * snpm_combo_pp
 Fix: Set bNeg to 0 for F-tests.
 
-#####  Updates from SnPM13.0.12
+####  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
+####  Updates from SnPM13.0.11
 * snpm_ui
 Fix: When relative thresholding is set, then global calculation (user defined or mean) must be computed.
 
@@ -119,22 +127,22 @@ 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 
+####  Updates from SnPM13.0.10 
 * snpm_ui
 Fix: ANCOVA normalisation (with an image presenting negative mean using spm_globals).
 
-#####  Updates from SnPM13.09
+####  Updates from SnPM13.09
 * snpm_cp
 Increase tolerance in checking that 'Locs_vox' is integer.
 
-#####  Updates from SnPM13.08
+####  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
+####  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.
 
@@ -144,14 +152,14 @@ Fix: slow cluster inference with variance smoothing.
 * generic_test_snpm, test_oneSample, test_twoSample
 Updates affecting tests only.
 
-#####  Updates from SnPM13.06
+####  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
+####  Updates from SnPM13.05
 * snpm_STcalc
 Cluster labels from spm_max are not necessarily continuous.
 
@@ -161,15 +169,15 @@ Fix: if number of requested permutations is equal to the maximum number of possi
 * 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
+####  Updates from SnPM13.04
 * snpm_ui
 Apply absolute and relative masking.
 
-#####  Updates from SnPM13.03
+####  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
+####  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.
@@ -183,7 +191,7 @@ Minor update in the display.
 * test/
 Miror corrections affecting tests procedure only.
 
-#####  Updates from SnPM13.01
+####  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.