FEAT: variable number of parfor workers dependning on template resolu…

…tion

Not enabled by default (-> uses Inf for all resolutions). Enable by specifying run.nworker as an array with the number of workers per resolution level, e.g., run.nworker = [Inf, 8] means that Inf number of workers will be used for all levels except the last (finest) resolution. This is useful when the template is large and too many workers can cause memory issues.

spm_mb_fit.m

@@ -17,15 +17,6 @@
 %--------------------------------------------------------------------------
 rng('default'); rng(1);
 
-% If SPM has been compiled with OpenMP support then the number of threads
-% are here set to speed up the algorithm
-%--------------------------------------------------------------------------
-if sett.nworker > 1
-    setenv('SPM_NUM_THREADS',sprintf('%d',0));
-else
-    setenv('SPM_NUM_THREADS',sprintf('%d',-1));
-end
-
 % Get template size and orientation
 %--------------------------------------------------------------------------
 if isfield(sett.mu,'exist')
@@ -40,6 +31,24 @@
 sz      = spm_mb_shape('zoom_settings',sett.v_settings,sett.mu,nz);
 sett.ms = sz(end);
 
+% Configure number of parfor workers
+%--------------------------------------------------------------------------
+nworker_zm = sett.nworker;
+if isscalar(nworker_zm), nworker_zm = [nworker_zm]; end
+if numel(nworker_zm) > numel(sz), nworker_zm = nworker_zm(1:numel(sz)); end
+nworker_zm = padarray(nworker_zm,[0, numel(sz) - numel(nworker_zm)],nworker_zm(1),'pre');
+nworker_zm = fliplr(nworker_zm);   
+sett.nworker = nworker_zm(end);    
+
+% If SPM has been compiled with OpenMP support then the number of threads
+% are here set to speed up the algorithm
+%--------------------------------------------------------------------------
+if sett.nworker > 1
+    setenv('SPM_NUM_THREADS',sprintf('%d',0));
+else
+    setenv('SPM_NUM_THREADS',sprintf('%d',-1));
+end
+
 % Init shape model parameters
 %--------------------------------------------------------------------------
 dat = spm_mb_shape('init_def',dat,sett.ms);
@@ -78,6 +87,7 @@
 spm_plot_convergence('Init','Rigid Alignment','Objective','Iteration');
 E      = Inf;
 for it0=1:nit_aff
+
     if it0>1
         oE  = E/nvox(dat);
     else
@@ -125,6 +135,9 @@
 for zm=numel(sz):-1:1 % loop over zoom levels
     fprintf('\nzoom=1/%d: %d x %d x %d\n', 2^(zm-1), sett.ms.d);
 
+    sett.nworker = nworker_zm(zm);
+    disp(['sett.nworker = ' num2str(sett.nworker)])
+
     if updt_mu
         dat = spm_mb_appearance('restart',dat,sett);
     end
@@ -202,6 +215,8 @@
     end
     do_save(mu,sett,dat);
 end
+% reset number of workers
+sett.nworker = nworker_zm(end);
 %spm_plot_convergence('Clear');
 %==========================================================================
 

spm_mb_output.m

@@ -24,6 +24,8 @@
 % If SPM has been compiled with OpenMP support then the number of threads
 % are here set to speed up the algorithm
 %--------------------------------------------------------------------------
+if isscalar(sett.nworker), sett.nworker = [sett.nworker]; end
+sett.nworker = sett.nworker(end);
 nw   = spm_mb_shape('get_num_workers',sett,max(27,sett.K*5+17));
 if sett.nworker > 1
     setenv('SPM_NUM_THREADS',sprintf('%d',0));