Wrap ground truth dual regression in try/catch

As per the convention for external calls to PROFUMO/MELODIC too.

RunTests.m

@@ -327,51 +327,63 @@ for n = 1:params.nRepeats
     
     % No obvious way to combine group maps, so skip if fewer inferred
     if params.iN >= params.N
-        
-        %Extract mean group map
-        Pg = 0;
-        for s = 1:params.S
-            Pg = Pg + P{s};
+        try
+            %Extract mean group map
+            Pg = 0;
+            for s = 1:params.S
+                Pg = Pg + P{s};
+            end
+            Pg = Pg / params.S;
+            % And save
+            mapFile = fullfile(repeatDir, 'GT_Pg.nii.gz');
+            save_avw( ...
+                reshape(Pg, 10, 10, params.V / 100, params.N), ...
+                mapFile, 'f', [1 1 1 1]);
+            
+            % Record actual size - don't include any extra maps asked for
+            % This is basically equivalent to setting extras to zero, but more
+            % numerically stable
+            gt_params = params; gt_params.iN = gt_params.N;
+            
+            % Dual-reg from ground truth group maps
+            dualregDir = fullfile(repeatDir, 'GT_Pg.dr');
+            system(sprintf( ...
+                'sh Methods/DualRegression.sh %s %s %s', ...
+                dualregDir, niftiDir, mapFile));
+            % `t` means thresholded as per Bijsterbosch et al., eLife, 2019
+            [gtdr_P, gtdr_A, gtdr_Pt, gtdr_At] ...
+                = loadDualRegression(dualregDir, gt_params);
+            gtdr_pcA  = calculateNetmats(gtdr_A, gt_params);
+            gtdr_pcAt = calculateNetmats(gtdr_At, gt_params);
+            
+            scores.GTg_DR(n) = calculateDecompositionAccuracy( ...
+                    P, gtdr_P, A, gtdr_A, pcA, gtdr_pcA, gt_params);
+            scores.GTg_DRt(n) = calculateDecompositionAccuracy( ...
+                    P, gtdr_Pt, A, gtdr_At, pcA, gtdr_pcAt, gt_params);
+            
+            % Add dummy maps for any extras inferred
+            %if params.iN > params.N
+            %    Pg = [Pg, zeros([params.V, params.iN - params.N])];
+            %end
+            
+            % Dual-reg from ground truth group maps
+            %[gtdr_P, gtdr_A] = runDR(D, Pg, params);
+            %gtdr_pcA = calculateNetmats(gtdr_A, gt_params);
+            %scores.GTg_DR(n) = calculateDecompositionAccuracy( ...
+            %        P, gtdr_P, A, gtdr_A, pcA, gtdr_pcA, gt_params);
+            
+            % Repeat, but now 'clean' dual regression (i.e. from PA, not D)
+            %[gtdr_P, gtdr_A] = runDR(PA, Pg, params);
+            %gtdr_pcA = calculateNetmats(gtdr_A, params);
+            %scores.GTg_PADR(n) = calculateDecompositionAccuracy( ...
+            %        P, gtdr_P, A, gtdr_A, pcA, gtdr_pcA, params);
+            
+            %Ground truth linear model
+            %scores.GT.PA(n) = calculateBoldRecovery(PA, makePA(P,A,params), params);
+            % N.B. This is essentially deprecated now there is noise in the PA subspace
+        catch
+            warning('Ground-truth dual regression failed.');
         end
-        Pg = Pg / params.S;
-        % Add dummy maps for any extras inferred
-        %if params.iN > params.N
-        %    Pg = [Pg, zeros([params.V, params.iN - params.N])];
-        %end
-        % Save
-        mapFile = fullfile(repeatDir, 'GT_Pg.nii.gz');
-        save_avw( ...
-            reshape(Pg, 10, 10, params.V / 100, params.N), ...
-            mapFile, 'f', [1 1 1 1]);
-        gt_params = params; gt_params.iN = gt_params.N;
-        
-        % Dual-reg from ground truth group maps
-        %[gtdr_P, gtdr_A] = runDR(D, Pg, params);
-        dualregDir = fullfile(repeatDir, 'GT_Pg.dr');
-        system(sprintf( ...
-            'sh Methods/DualRegression.sh %s %s %s', ...
-            dualregDir, niftiDir, mapFile));
-        % `t` means thresholded as per Bijsterbosch et al., eLife, 2019
-        [gtdr_P, gtdr_A, gtdr_Pt, gtdr_At] ...
-            = loadDualRegression(dualregDir, gt_params);
-        gtdr_pcA  = calculateNetmats(gtdr_A, gt_params);
-        gtdr_pcAt = calculateNetmats(gtdr_At, gt_params);
-        
-        scores.GTg_DR(n) = calculateDecompositionAccuracy( ...
-                P, gtdr_P, A, gtdr_A, pcA, gtdr_pcA, gt_params);
-        scores.GTg_DRt(n) = calculateDecompositionAccuracy( ...
-                P, gtdr_Pt, A, gtdr_At, pcA, gtdr_pcAt, gt_params);
-        
-        % Repeat, but now 'clean' dual regression (i.e. from PA, not D)
-        %[gtdr_P, gtdr_A] = runDR(PA, Pg, params);
-        %gtdr_pcA = calculateNetmats(gtdr_A, params);
-        
-        %scores.GTg_PADR(n) = calculateDecompositionAccuracy( ...
-        %        P, gtdr_P, A, gtdr_A, pcA, gtdr_pcA, params);
-        
-        %Ground truth linear model
-        %scores.GT.PA(n) = calculateBoldRecovery(PA, makePA(P,A,params), params);
-        % N.B. This is essentially deprecated now there is noise in the PA subspace
         
     end