Switch to PROFUMO postprocessing routines

Rather than replicating these in MATLAB, we now run the Python-based
postprocessing routines for PROFUMO. This makes things considerably
simpler, especially if we want to add features to postprocessing.

IO/loadPROFUMO.m

@@ -2,49 +2,36 @@ function [ Pg, P, A, pcA ] = loadPROFUMO( pfmDir, params )
 
 % SPDX-License-Identifier: Apache-2.0
 
-pfmDir = fullfile(pfmDir, 'FinalModel');
+pfmDir = fullfile(pfmDir, 'Results.ppp');
 
 % Load group maps
-mPg = h5read(fullfile(pfmDir, 'GroupSpatialModel', 'SignalMeans.post', ...
-    'Means.hdf5'), '/dataset');
-pPg = h5read(fullfile(pfmDir, 'GroupSpatialModel', 'Memberships.post', ...
-    'Class_2', 'Probabilities.hdf5'), '/dataset');
-Pg = mPg .* pPg; clear mPg pPg;
+Pg = read_avw(fullfile(pfmDir, 'Maps', 'Group.nii.gz'));
+Pg = reshape(Pg, params.V, params.iN);
 
 % And subject maps / timecourses
 P = cell(params.S, 1);
 A = cell(params.S, 1);
 pcA = cell(params.S, 1);
 for s = 1:params.S
-    subj = sprintf('S%02d',s);
-    subjDir = fullfile(pfmDir, 'Subjects', subj);
+    subj = sprintf('sub-S%02d',s);
     
     % Subject maps
-    mPs = h5read(fullfile(subjDir, 'SpatialMaps.post', 'Signal', ...
-        'Means.hdf5'), '/dataset');
-    pPs = h5read(fullfile(subjDir, 'SpatialMaps.post', 'Signal', ...
-        'MembershipProbabilities.hdf5'), '/dataset');
-    P{s} = mPs .* pPs; clear mPs pPs;
+    P{s} = read_avw(fullfile(pfmDir, 'Maps', [subj '.nii.gz']));
+    P{s} = reshape(P{s}, params.V, params.iN);
     
     % Time courses
     A{s} = cell(params.R(s), 1);
     pcA{s} = cell(params.R(s), 1);
     for r = 1:params.R(s)
-        run = sprintf('R%02d',r);
-        runDir = fullfile(subjDir, 'Runs', run);
+        run = sprintf('run-R%02d',r);
         
         % Timecourses
-        A{s}{r} = h5read(fullfile(runDir, 'TimeCourses.post', ...
-            'CleanTimeCourses', 'Means.hdf5'), '/dataset');
-            %'Means.hdf5'), '/dataset');
+        A{s}{r} = csvread(fullfile(pfmDir, 'TimeCourses', [subj '_' run '.csv']))';
         % Modulate by amplitudes
-        h = h5read(fullfile(runDir, 'ComponentWeightings.post', ...
-            'Means.hdf5'), '/dataset');
+        h = csvread(fullfile(pfmDir, 'Amplitudes', [subj '_' run '.csv']));
         A{s}{r} = A{s}{r} .* h;
         
         % Temporal netmats
-        precmat = h5read(fullfile(runDir, 'TemporalPrecisionMatrix.post', ...
-            'Mean.hdf5'), '/dataset');
-        pcA{s}{r} = - corrcov(precmat);
+        pcA{s}{r} = csvread(fullfile(pfmDir, 'NetMats', [subj '_' run '.csv']));
     end
 end

IO/saveNIfTIs.m

@@ -8,7 +8,7 @@ assert(mod(params.V, 100) == 0, '`params.V` must be divisible by 100!')
 vsize = [1 1 1 params.TR];
 vtype = 'f';
 
-% Loop over subjects, saving NIfTIs and creating MELODIC/PROFUMO sped files
+% Loop over subjects, saving NIfTIs and creating MELODIC/PROFUMO spec files
 data = struct();
 fileNames = {};
 for s = 1:params.S
@@ -26,6 +26,10 @@ for s = 1:params.S
     end
 end
 
+% Save a reference NIfTI
+fileName = fullfile(directory, 'Reference.nii.gz');
+save_avw(ones(10, 10, params.V / 100), fileName, vtype, [1 1 1 1]);
+
 % Save spec files
 fileID = fopen(fullfile(directory, 'PROFUMO_SpecFile.json'), 'w');
 fprintf(fileID, '%s\n', jsonencode(data));

Methods/PROFUMO.sh

@@ -12,12 +12,21 @@ TR=$4
 
 mkdir "${output_dir}"
 
-# Run PROFUMO (need to make sure we're on jalapeno18!)
-#nice -n 20 PROFUMO \
-PROFUMO \
+# Will need adjusting!
+profumo_dir=~/FMRIB/Code/PROFUMO
+
+# Run PROFUMO analysis
+"${profumo_dir}/C++/PROFUMO" \
     "${nifti_dir}/PROFUMO_SpecFile.json" \
-    ${dim} \
-    "${output_dir}" \
-    --useHRF ${TR} --hrfFile ~samh/PROFUMO/HRFs/Default.phrf \
+    "${dim}" \
+    "${output_dir}/Analysis.pfm" \
+    --useHRF "${TR}" --hrfFile "${profumo_dir}/HRFs/Default.phrf" \
     --covModel Run --lowRankData 150 -d 0.5 --globalInit --nThreads 10 \
-    > "${output_dir}/TerminalOutput.txt"
+    > "${output_dir}/TerminalOutput-pfm.txt"
+
+# Run postprocessing
+"${profumo_dir}/Python/postprocess_results.py" \
+    "${output_dir}/Analysis.pfm" \
+    "${output_dir}/Results.ppp" \
+    "${nifti_dir}/Reference.nii.gz" \
+    > "${output_dir}/TerminalOutput-ppp.txt"

RunTests.m

@@ -436,7 +436,7 @@ for n = 1:params.nRepeats
     %% Run PROFUMO %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
     try
-        pfmDir = fullfile(repeatDir, 'PROFUMO.pfm');
+        pfmDir = fullfile(repeatDir, 'PROFUMO');
         system(sprintf( ...
             'sh Methods/PROFUMO.sh %s %s %d %1.2f', ...
             pfmDir, niftiDir, params.iN, params.TR));