Add parallel processing option

Allow for use of multiple processing via parfor. Also changed progress
output to work with parfor.

demos/fwe_control.m

@@ -90,4 +90,4 @@
         fp(i,j) = sum(pcorr_pos{j}(:)<.05)+sum(pcorr_neg{j}(:)<.05);
     end
 end
-sum(fp)/nsim
+fp/nsim

matlab_tfce.m

@@ -25,10 +25,10 @@
 % corrections, is substantially less conservative due to the fact that
 % it capitalizes on spatial dependency in the data. 
 %
-% [varargout] = matlab_tfce(analysis,tails,imgs,imgs2,covariate,nperm,H,E,C,dh)
-% [pcorr] = matlab_tfce(analysis,1,imgs,imgs2,covariate,nperm,H,E,C,dh)
-% [pcorr_pos,pcorr_neg] = matlab_tfce(analysis,2,imgs,imgs2,covariate,nperm,H,E,C,dh)
-% [pcorr_fac1,pcorr_fac2,pcorr_int] = matlab_tfce(analysis,1,imgs,imgs2,covariate,nperm,H,E,C,dh)
+% [varargout] = matlab_tfce(analysis,tails,imgs,imgs2,covariate,nperm,H,E,C,dh,parworkers)
+% [pcorr] = matlab_tfce(analysis,1,imgs,imgs2,covariate,nperm,H,E,C,dh,parworkers)
+% [pcorr_pos,pcorr_neg] = matlab_tfce(analysis,2,imgs,imgs2,covariate,nperm,H,E,C,dh,parworkers)
+% [pcorr_fac1,pcorr_fac2,pcorr_int] = matlab_tfce(analysis,1,imgs,imgs2,covariate,nperm,H,E,C,dh,parworkers)
 %
 % Arguments:
 %
@@ -42,7 +42,8 @@
 %   -- 'regression' -- multiple linear regression with covariate matrix
 %
 % tails -- specify a 1 or 2 tailed test (unidirectional or bidirectional)
-% that can be combined with t-tests and correlations.
+% that can be combined with t-tests, correlations, and regressions. Ignored
+% for repeated measures ANOVAs.
 %
 % imgs -- a 4D matrix of imaging data for analysis. Dimensions are expected
 % to be x,y,z,subject. Alternatively, for repeated measures ANOVAs, a cell
@@ -81,6 +82,13 @@
 %   The C default setting matches FSL's ramdomise default setting. To
 %   match SPM's default cluster forming, use 18 instead.
 %
+% parworkers -- number of parallel workers to use, default = 0 (no
+% parallelization). If a parallel pool is already open, its size will
+% override this argument. If a new pool is opened for MatlabTFCE, it will
+% be closed at the end of script. Opening and closing parallel pools is
+% time consuming so it may be more efficient to do outside the package if
+% one is running several analyses.
+%
 % Output: 
 % If tails == 1, a single output image with the same xyz dimensions as imgs
 % consisting of corrected p-values with be returned.
@@ -116,6 +124,8 @@
 E = .5;
 C = 26;
 dh = .1;
+global parworkers;
+parworkers = 0;
 
 % adjusting optional arguments based on input
 fixedargn = 3;
@@ -150,6 +160,11 @@
         dh = varargin{7};
     end
 end
+if nargin > (fixedargn + 7)
+    if ~isempty(varargin{8})
+        parworkers = varargin{8};
+    end
+end
 
 % check that tails are appropriate
 if ~(sum(tails==[1 2]))
@@ -241,12 +256,28 @@
     end  
 end
 
+% check if parallel pool exists; use it if it does, create it if not (but
+% is requested).
+pp = gcp('nocreate');
+try 
+    parworkers = pp.NumWorkers;
+    parexist = 1;
+    warning('Existing parallel pool discovered with %d workers. Overriding parworkers argument to use this pool.',parworkers);
+catch
+    parexist = 0;
+    if parworkers > 0
+        parpool(parworkers);
+    end
+end
+
 %% analysis calls
 
 % indicate analysis start
 fprintf('Starting %i-tailed %s analysis...\n',[tails,analysis]);
+fprintf('MatlabTFCE started at %s \n',datestr(datetime('now')));
 fprintf('Permuting data from %i subjects %i times\n',[nsub,nperm])
 
+
 % select appropriate analysis
 switch analysis
 
@@ -317,7 +348,12 @@
     varargout{1} = pcorr_pos;
     varargout{2} = pcorr_neg;
 end
-fprintf('Permutation TFCE complete!\n');
+fprintf('Permutation complete!\n');
+fprintf('MatlabTFCE finished at %s \n',datestr(datetime('now')));
 
+%% clean up parallel pool if created in-script
+if ~parexist
+    delete(gcp('nocreate'));
+end
 end
 

matlab_tfce_correlation.m

@@ -59,12 +59,12 @@
 end
 
 % initialize progress indicator
-fprintf('Completed: %3d%%', 0);
-indicatorSteps = round(nperm/100);
+parfor_progress(nperm);
+global parworkers
 
 % cycle through permutations
 exceedances = zeros(nvox,1);
-for p = 1:nperm
+parfor(p = 1:nperm,parworkers)
 
     % permute covariates
     rsel = randperm(nsub);
@@ -84,9 +84,9 @@
     curexceeds = max(rstats) >= cvals;
     exceedances = exceedances + curexceeds;
 
-    % update progress indicator every percentage point
-    if ~mod(p,indicatorSteps) || p==nperm
-        fprintf(sprintf('%s%%3d%%%%', repmat('\b', 1, 4)), round(100*p/nperm));
+    % update progress indicator (only does so 1 in 5 to minimize overhead)
+    if ~randi([0 4]);
+        parfor_progress;
     end
 
 end

matlab_tfce_gui.m

@@ -115,35 +115,39 @@ function anadiag(~,~,analysis)
                       'Toolbar','none', ...
                       'NumberTitle','off', ...
                       'Resize', 'off', ...
-                      'Position',[450,400,220,420]); %X,Y then Width, Height
+                      'Position',[450,450,220,420]); %X,Y then Width, Height
 set(anafig, 'Color', ([128,0,0] ./255)); 
 yp = 20;
 
 uipanel('Title','Contrast name',...
              'BackgroundColor','white',...
              'Units','Pixels',...
-             'Position',[30 340+yp 160 40]);
+             'Position',[30 390+yp 160 40]);
 uipanel('Title','Tails of test',...
              'BackgroundColor','white',...
              'Units','Pixels',...
-             'Position',[30 290+yp 160 40]);
+             'Position',[30 340+yp 160 40]);
 uipanel('Title','Permutation number',...
              'BackgroundColor','white',...
              'Units','Pixels',...
-             'Position',[30 240+yp 160 40]);
+             'Position',[30 290+yp 160 40]);
 uipanel('Title','TFCE H parameter',...
              'BackgroundColor','white',...
              'Units','Pixels',...
-             'Position',[30 190+yp 160 40]);
+             'Position',[30 240+yp 160 40]);
 uipanel('Title','TFCE E parameter',...
              'BackgroundColor','white',...
              'Units','Pixels',...
-             'Position',[30 140+yp 160 40]);
+             'Position',[30 190+yp 160 40]);
 uipanel('Title','Connectivity',...
              'BackgroundColor','white',...
              'Units','Pixels',...
-             'Position',[30 90+yp 160 40]);
+             'Position',[30 140+yp 160 40]);
 uipanel('Title','TFCE dh parameter',...
+             'BackgroundColor','white',...
+             'Units','Pixels',...
+             'Position',[30 90+yp 160 40]);
+uipanel('Title','Size of parallel pool',...
              'BackgroundColor','white',...
              'Units','Pixels',...
              'Position',[30 40+yp 160 40]);
@@ -156,12 +160,13 @@ function anadiag(~,~,analysis)
 else
     u2 = uicontrol('Style','popupmenu','HorizontalAlignment','left','String',{'1'},'Position',[50,300+yp,120,20]);
 end
-u3 = uicontrol('Style','edit','HorizontalAlignment','left','String','5000','Position',[50,250+yp,120,20]);
-u4 = uicontrol('Style','edit','HorizontalAlignment','left','String','2','Position',[50,200+yp,120,20]);
-u5 = uicontrol('Style','edit','HorizontalAlignment','left','String','0.5','Position',[50,150+yp,120,20]);
-u6 = uicontrol('Style','popupmenu','HorizontalAlignment','left','String',{'26','18','6'},'Position',[50,100+yp,120,20]);
-u7 = uicontrol('Style','edit','HorizontalAlignment','left','String','0.1','Position',[50,50+yp,120,20]);
-uicontrol('Style','PushButton','HorizontalAlignment','center','String','Begin','Position',[30,10,160,40],'Callback',{@run_tfce,[u1 u2 u3 u4 u5 u6 u7],analysis,niiout,imgs,covariate});
+u3 = uicontrol('Style','edit','HorizontalAlignment','left','String','5000','Position',[50,300+yp,120,20]);
+u4 = uicontrol('Style','edit','HorizontalAlignment','left','String','2','Position',[50,250+yp,120,20]);
+u5 = uicontrol('Style','edit','HorizontalAlignment','left','String','0.5','Position',[50,200+yp,120,20]);
+u6 = uicontrol('Style','popupmenu','HorizontalAlignment','left','String',{'26','18','6'},'Position',[50,150+yp,120,20]);
+u7 = uicontrol('Style','edit','HorizontalAlignment','left','String','0.1','Position',[50,100+yp,120,20]);
+u8 = uicontrol('Style','edit','HorizontalAlignment','left','String','0','Position',[50,50+yp,120,20]);
+uicontrol('Style','PushButton','HorizontalAlignment','center','String','Begin','Position',[30,10,160,40],'Callback',{@run_tfce,[u1 u2 u3 u4 u5 u6 u7,u8],analysis,niiout,imgs,covariate});
 
 end
 
@@ -231,6 +236,7 @@ function run_tfce(~,~,uis,analysis,niiout,imgs,covariate)
 E = str2double(strings{5});
 C = str2double(strings{6}{values(6)});
 dh = str2double(strings{7});
+parworkers = str2double(strings{8});
 if ~(strcmp(analysis,'rm_anova1') || strcmp(analysis,'rm_anova2')) && length(imgs) == 2 
     imgs2 = imgs{2};
     imgs = imgs{1};
@@ -242,11 +248,11 @@ function run_tfce(~,~,uis,analysis,niiout,imgs,covariate)
 % run analysis
 disp(['Analysis ' ananame ' ' analysis ' started at ' datestr(now)]);
 if strcmp(analysis,'rm_anova2')
-    [pcorr_fac1,pcorr_fac2,pcorr_int] = matlab_tfce(analysis,tails,imgs,imgs2,covariate,nperm,H,E,C,dh);
+    [pcorr_fac1,pcorr_fac2,pcorr_int] = matlab_tfce(analysis,tails,imgs,imgs2,covariate,nperm,H,E,C,dh,parworkers);
 elseif tails == 1
-    pcorr = matlab_tfce(analysis,tails,imgs,imgs2,covariate,nperm,H,E,C,dh);
+    pcorr = matlab_tfce(analysis,tails,imgs,imgs2,covariate,nperm,H,E,C,dh,parworkers);
 else
-    [pcorr_pos,pcorr_neg] = matlab_tfce(analysis,tails,imgs,imgs2,covariate,nperm,H,E,C,dh);
+    [pcorr_pos,pcorr_neg] = matlab_tfce(analysis,tails,imgs,imgs2,covariate,nperm,H,E,C,dh,parworkers);
 end
 disp(['Analysis ' ananame ' ' analysis ' completed at ' datestr(now)]);
 

matlab_tfce_regression.m

@@ -82,12 +82,12 @@
 end
 
 % initialize progress indicator
-fprintf('Completed: %3d%%', 0);
-indicatorSteps = round(nperm/100);
+parfor_progress(nperm);
+global parworkers
 
 % cycle through permutations
 exceedances = zeros(nvox,npred);
-for p = 1:nperm,parworkers
+parfor(p = 1:nperm,parworkers)
 
     % permute predictors
     rsel = randperm(nsub);
@@ -110,6 +110,7 @@
     end
 
     rtstats = bs./serr;
+    curexceeds = zeros(nvox,npred);
     for i = 1:npred
         rbrain=NaN(bsize);
         rbrain(implicitmask) = rtstats(i,:);
@@ -119,13 +120,13 @@
             rstats = abs(rstats);
         end  
         % compare maxima to t-values and increment as appropriate
-        curexceeds = max(rstats) >= tvals(i,:);
-        exceedances(:,i) = exceedances(:,i) + curexceeds';
+        curexceeds(:,i) = max(rstats) >= tvals(i,:);
     end
+    exceedances = exceedances + curexceeds;
 
-    % update progress indicator every percentage point
-    if ~mod(p,indicatorSteps) || p==nperm
-        fprintf(sprintf('%s%%3d%%%%', repmat('\b', 1, 4)), round(100*p/nperm));
+    % update progress indicator (only does so 1 in 5 to minimize overhead)
+    if ~randi([0 4]);
+        parfor_progress;
     end
 
 end

matlab_tfce_rm_anova1.m

@@ -73,27 +73,28 @@
 tfcestat = tfcestat(implicitmask);
 
 % initialize progress indicator
-fprintf('Completed: %3d%%', 0);
-indicatorSteps = round(nperm/100);
+parfor_progress(nperm);
+global parworkers
 
 % cycle through permutations
 exceedances = zeros(nvox,1);
-rbrain = zeros(bsize);
-for p = 1:nperm
+parfor(p = 1:nperm,parworkers)
 
     % permute labels
+    roccimgs = NaN(size(occimgs));
     for s = 1:nsub
         relabeling = randperm(levels)';
-        occimgs(s,relabeling,:) = occimgs(s,:,:);
+        roccimgs(s,relabeling,:) = occimgs(s,:,:);
     end
 
     % calculate permutation statistic
-    rcmeans = mean(occimgs,1);
+    rcmeans = mean(roccimgs,1);
     rSSfac = sum((rcmeans-gmeans_c).^2,2)*nsub;
     rSSerr = SSw-rSSfac;
     rMSfac = rSSfac./(levels-1);
     rMSerr = rSSerr./((levels-1)*(nsub-1));
     rfvals = rMSfac./rMSerr;
+    rbrain = zeros(bsize);
     rbrain(implicitmask) = rfvals;
     rbrain = transform(rbrain,H,E,C,dh);
     rstats = rbrain(implicitmask);
@@ -102,10 +103,11 @@
     curexceeds = max(rstats) >= tfcestat;
     exceedances = exceedances + curexceeds;
 
-    % update progress indicator every percentage point
-    if ~mod(p,indicatorSteps) || p==nperm
-        fprintf(sprintf('%s%%3d%%%%', repmat('\b', 1, 4)), round(100*p/nperm));
+    % update progress indicator (only does so 1 in 5 to minimize overhead)
+    if ~randi([0 4]);
+        parfor_progress;
     end
+
 end
 
 % create corrected p-value image

matlab_tfce_rm_anova2.m

@@ -107,31 +107,29 @@
 tfcestati = tfcestati(implicitmask);
 
 % initialize progress indicator
-fprintf('Completed: %3d%%', 0);
-indicatorSteps = round(nperm/100);
+parfor_progress(nperm);
+global parworkers
 
 % cycle through permutations
 exceedances1 = zeros(nvox,1);
 exceedances2 = zeros(nvox,1);
 exceedancesi = zeros(nvox,1);
-rbrain1 = zeros(bsize);
-rbrain2 = zeros(bsize);
-rbraini = zeros(bsize);
-for p = 1:nperm
+parfor(p = 1:nperm,parworkers)
 
     % permute labels
+    roccimgs = NaN(size(occimgs));
     for s = 1:nsub
         relabeling1 = randperm(levels(1))';
         relabeling2 = randperm(levels(2))';
-        occimgs(s,relabeling1,relabeling2,:) = occimgs(s,:,:,:);
+        roccimgs(s,relabeling1,relabeling2,:) = occimgs(s,:,:,:);
     end
 
     % calculate permuted means
-    rc1means = mean(mean(occimgs,1),3);
-    rc2means = mean(mean(occimgs,1),2);
-    rc1smeans = mean(occimgs,3);
-    rc2smeans = mean(occimgs,2);
-    rccmeans = mean(occimgs,1);
+    rc1means = mean(mean(roccimgs,1),3);
+    rc2means = mean(mean(roccimgs,1),2);
+    rc1smeans = mean(roccimgs,3);
+    rc2smeans = mean(roccimgs,2);
+    rccmeans = mean(roccimgs,1);
 
     % calculate SS
     rSSfac1 = sum((rc1means-gmeans_c1).^2,2)*levels(2)*nsub;
@@ -158,6 +156,9 @@
     rfvalsi = rMSint./rMSerr;
 
     % perform TFCE
+    rbrain1 = zeros(bsize);
+    rbrain2 = zeros(bsize);
+    rbraini = zeros(bsize);
     rbrain1(implicitmask) = rfvals1;
     rbrain2(implicitmask) = rfvals2;
     rbraini(implicitmask) = rfvalsi;
@@ -176,9 +177,9 @@
     exceedances2 = exceedances2 + curexceeds2;
     exceedancesi = exceedancesi + curexceedsi;
 
-    % update progress indicator every percentage point
-    if ~mod(p,indicatorSteps) || p==nperm
-        fprintf(sprintf('%s%%3d%%%%', repmat('\b', 1, 4)), round(100*p/nperm));
+    % update progress indicator (only does so 1 in 5 to minimize overhead)
+    if ~randi([0 4]);
+        parfor_progress;
     end
 
 end