Fixes for v2.2.7, mostly parallelization

docs/harmonisation.md

@@ -90,6 +90,7 @@ You can tune this command using additional variables and flags as detailed bello
 |```-demos <demographic_file>```| The name of the csv file containing the demographic information as detailled in the [guidelines](https://meld-graph.readthedocs.io/en/latest/prepare_data.html#prepare-the-demographic-information-required-only-to-compute-the-harmonisation-parameters). An example 'demographics_file.csv' is provided in the <meld_data_folder>.|
 | **Optional variables** |
 |```--parallelise``` | use this flag to speed up the segmentation by running Freesurfer/FastSurfer on multiple subjects in parallel. |
+|```--threads <number_of_threads>``` | use this variable to speed up the segmentation of a single subject by running Freesurfer with several threads (OpenMP). Contrary to ```--parallelise```, the subjects are still processed one after another. Only applies to Freesurfer, and is ignored when ```--parallelise``` is used. |
 |```--fastsurfer``` | use this flag to use FastSurfer instead of Freesurfer. Requires FastSurfer installed. |
 |```--harmo_only``` | Use this flag to do all the processes up to the harmonisation. Usefull if you want to harmonise on some subjects but do not wish to predict on them |
 

docs/interpret_results.md

@@ -114,6 +114,8 @@ After viewing the MELD PDF report, it is then important to visualise the predict
 The predictions are saved as NIFTI files in the folder: 
 /output/predictions_reports/<subject_id>/predictions
 
+The file to open is `prediction.nii.gz` (both hemispheres combined). The subfolders of that directory contain the intermediate files created while registering the prediction back into the native volume, see [the description of script 3](https://meld-graph.readthedocs.io/en/latest/run_prediction_pipeline.html#script-3-lesions-prediction-meld-reports).
+
 **It is important to check that the clusters detected are not due to obvious FreeSurfer reconstruction errors, scan artifacts etc.**
 
 To quality control the predicted lesion and the surfaces, you can use the HTML viewer mentioned above or use Freeview (you will need a standalone version of freeview)

docs/run_prediction_pipeline.md

@@ -84,8 +84,11 @@ You can tune the MELD pipeline command using additional variables and flags as d
 | **Optional variables** |
 | ```-harmo_code <harmo_code>```  | provide the harmonisation code if you want to harmonise your data before prediction. This requires to have [computed the harmonisation parameters](https://meld-graph.readthedocs.io/en/latest/harmonisation.html) beforehand. The harmonisation code should start with H, e.g. H1. | 
 |```--parallelise``` | use this flag to speed up the segmentation by running Freesurfer/FastSurfer on multiple subjects in parallel. |
+|```--threads <number_of_threads>``` | use this variable to speed up the segmentation of a single subject by running Freesurfer with several threads (OpenMP). Contrary to ```--parallelise```, the subjects are still processed one after another. Only applies to Freesurfer, and is ignored when ```--parallelise``` is used. |
 |```--fastsurfer``` | use this flag to use FastSurfer instead of Freesurfer. (Requires FastSurfer installed for native installation). |
 |```--skip_feature_extraction``` | use this flag to skips the segmentation and features extraction (processes from script1). Usefull if you already have these outputs and you just want to run the preprocessing and the predictions (e.g: after harmonisation) |
+|```--skip_prediction``` | use this flag to skip the predictions (processes from script3). The segmentation and the preprocessing are still run. This is different from ```--harmo_only```, which also skips the feature extraction steps that are only needed for the predictions. |
+|```--skip_feature_plotting``` | use this flag to skip the plotting of the features for QC during the preprocessing (processes from script2), as this can take a long time. |
 |**More advanced variables** | 
 |```--no_nifti```| use this flag to run to all the processes up saving the predictions as surface vectors in the hdf5 file. Does not produce produce nifti and pdf outputs.|
 |```--no_report``` | use this flag to do all the processes up to creating the prediction as a nifti file. Does not produce the pdf reports. |
@@ -172,6 +175,30 @@ singularity exec meld_graph.sif /bin/bash -c "cd /app && source \$FREESURFER_HOM
 :::
 ::::
 
+To run a large number of subjects in parallel, e.g. on a multi-core machine or on a cluster, you can call the 3 scripts of the pipeline (see [below](#additional-information-about-the-pipeline)) yourself and parallelise the two steps that work on one subject at a time. Compared to ```--parallelise```, this lets you give several threads to each Freesurfer reconstruction and lets your job scheduler handle the queueing.
+
+**Step 1 - segmentation of each subject, one job per subject.** Each job uses ```--threads``` cores, so keep (number of jobs) x (number of threads) below the number of cores available:
+
+```bash
+python scripts/new_patient_pipeline/run_script_segmentation.py -id <subject_id> --threads <number_of_threads>
+```
+
+**Step 2 - feature extraction and harmonisation, once for all the subjects.** The harmonisation parameters are computed from all the subjects of a scanner/site together, so this step cannot be split per subject:
+
+```bash
+python scripts/new_patient_pipeline/new_pt_pipeline.py -ids input/subjects_list.txt -harmo_code <harmo_code> -demos input/demographics_file.csv --skip_prediction --skip_feature_plotting
+```
+
+The Freesurfer reconstructions of step 1 are detected and skipped, the features of each subject are extracted into the feature matrix of the scanner, and the predictions and the QC feature plots are left out.
+
+**Step 3 - predictions and reports, one job per subject**, parallelised in the same way as step 1:
+
+```bash
+python scripts/new_patient_pipeline/run_script_prediction.py -id <subject_id> -harmo_code <harmo_code> -demos input/demographics_file.csv
+```
+
+Adapt this for containerization (usually apptainer/singularity for HPC) and your job scheduler (e.g., SLURM) as appropriate for your system.
+
 ## Additional information about the pipeline
 
 The pipeline is split into 3 main scripts as illustrated below and detailed in the next sub-sections. 
@@ -316,9 +343,33 @@ singularity exec meld_graph.sif /bin/bash -c "cd /app && source \$FREESURFER_HOM
 This script : 
 1. Run the MELD classifier and predict lesion on new subject
 2. Register the prediction back into the native nifti MRI. Results are stored in output/predictions_reports/<subjec_id>/predictions.
-3. Create MELD reports with predicted lesion location on inflated brain, on native MRI and associated saliencies. Reports are stored in output/predictions_reports/<subjec_id>/predictions/reports.
+3. Create MELD reports with predicted lesion location on inflated brain, on native MRI and associated saliencies. Reports are stored in output/predictions_reports/<subjec_id>/reports.
+
+The predictions folder contains the final prediction in a nifti volume (aligned with the input T1w) as well as intermediate files (prediction on the template fsaverage_sym surface space, the subjects native surface space, etc.):
+
+```
+output/predictions_reports/<subject_id>/
+├── predictions/
+│   ├── prediction.nii.gz         <- final prediction (both hemispheres), in the native T1w space
+│   ├── lh.prediction.nii.gz      <- left hemisphere only
+│   ├── rh.prediction.nii.gz      <- right hemisphere only
+│   ├── fsaverage_sym/            <- intermediate: prediction on the template surface
+│   │   ├── lh.prediction.mgh
+│   │   └── rh.prediction.mgh
+│   ├── surf_native/              <- intermediate: prediction on the native surface
+│   │   ├── lh.prediction.mgh
+│   │   └── rh.prediction.mgh
+│   ├── vol_freesurfer/           <- intermediate: prediction in the Freesurfer conformed volume
+│   │   ├── lh.prediction.mgz
+│   │   └── rh.prediction.mgz
+│   └── vol_native/               <- intermediate: prediction in the native T1 volume, mgz format
+│       ├── lh.prediction.mgz
+│       └── rh.prediction.mgz
+└── reports/                      <- MELD PDF report and the figures it is built from
+```
 
 Notes: 
+- The in previous versions, these intermediate files were stored in the Freesurfer subject directory, so there was a risk of multiple predictions using the same Freesurfer outputs (e.g. with different harmonisation parameters) overwriting each other.
 - Features need to have been processed using script 2 and Freesurfer outputs need to be available for each subject
 
 Example to use it on one patient without harmonisation:

meld_graph/data_preprocessing.py

@@ -23,6 +23,7 @@
 import meld_graph.distributedCombat as dc
 import meld_graph.mesh_tools as mt
 import meld_graph.meld_plotting as mpt
+from meld_graph.hdf5_utils import open_hdf5_file
 
 class Preprocess:
     """
@@ -599,9 +600,9 @@ def load_covars(self, subject_ids=None, demographic_file=DEMOGRAPHIC_FEATURES_FI
     def save_norm_combat_parameters(self, feature, estimates, hdf5_file):
         """Save estimates from combat and normalisation parameters on hdf5"""
         if not os.path.isfile(hdf5_file):
-            hdf5_file_context = h5py.File(hdf5_file, "a")
+            hdf5_file_context = open_hdf5_file(hdf5_file, mode="a", create_parent=True)
         else:
-            hdf5_file_context = h5py.File(hdf5_file, "r+")
+            hdf5_file_context = open_hdf5_file(hdf5_file, mode="r+")
 
         with hdf5_file_context as f:
             list_params = list(set(estimates))
@@ -625,7 +626,7 @@ def save_norm_combat_parameters(self, feature, estimates, hdf5_file):
 
     def read_norm_combat_parameters(self, feature, hdf5_file):
         """reconstruct estimates dictionnary from the combat parameters hdf5 file"""
-        hdf5_file_context = h5py.File(hdf5_file, "r")
+        hdf5_file_context = open_hdf5_file(hdf5_file, mode="r")
         estimates = {}
         with hdf5_file_context as f:
             feat_dir = f[feature]

meld_graph/evaluation.py

@@ -1,5 +1,6 @@
 import logging
 import os
+import tempfile
 import torch
 import torch_geometric.data
 from meld_graph.dataset import GraphDataset
@@ -18,6 +19,7 @@
 import itertools
 import seaborn as sns
 from meld_graph.paths import MELD_DATA_PATH
+from meld_graph.hdf5_utils import open_hdf5_file
 
 # for saliency - do not force people to have this
 try:
@@ -906,27 +908,17 @@ def save_prediction(self, subject, prediction, dataset_str="prediction", dtype=N
             dtype = prediction.dtype
 
         filename = os.path.join(self.results_dir, f"predictions{suffix}.hdf5")
-        if not os.path.isfile(filename):
-            mode = "a"
-        else:
-            mode = "r+"
-        done = False
-        while not done:
-            try:
-                with h5py.File(filename, mode=mode) as f:
-                    self.log.info(f"saving {dataset_str} for {subject}")
-                    for i, hemi in enumerate(["lh", "rh"]):
-                        shape = tuple([nvert_hemi] + list(prediction.shape[1:]))
-                        # create dataset
-                        dset = f.require_dataset(f"{subject}/{hemi}/{dataset_str}", shape=shape, dtype=dtype)
-                        # save prediction in dataset
-                        dset[:] = prediction[i * nvert_hemi : (i + 1) * nvert_hemi]
-                        # if dataset_str == "prediction":
-                        # save threshold as attribute in dataset
-                        # dset.attrs["threshold"] = self.threshold
-                    done = True
-            except OSError:
-                done = False
+        with open_hdf5_file(filename, mode="a", create_parent=True) as f:
+            self.log.info(f"saving {dataset_str} for {subject}")
+            for i, hemi in enumerate(["lh", "rh"]):
+                shape = tuple([nvert_hemi] + list(prediction.shape[1:]))
+                # create dataset
+                dset = f.require_dataset(f"{subject}/{hemi}/{dataset_str}", shape=shape, dtype=dtype)
+                # save prediction in dataset
+                dset[:] = prediction[i * nvert_hemi : (i + 1) * nvert_hemi]
+                # if dataset_str == "prediction":
+                # save threshold as attribute in dataset
+                # dset.attrs["threshold"] = self.threshold
 
     def load_prediction(self, subject, dataset_str="prediction", suffix=""):
         """
@@ -937,7 +929,7 @@ def load_prediction(self, subject, dataset_str="prediction", suffix=""):
             # cannot load data
             self.log.debug(f'file {filename} does not exist')
             return None
-        with h5py.File(filename, mode='r') as f:
+        with open_hdf5_file(filename, mode="r") as f:
             prediction = []
             try:
                 for i, hemi in enumerate(["lh", "rh"]):
@@ -1155,22 +1147,23 @@ def create_surface_plots(coords, faces, overlay, flat_map=True, limits=None):
     else:
         vmin = limits[0]
         vmax = limits[1]
-    tmp_file = os.path.join(MELD_DATA_PATH,'tmp.png')
-    msp.plot_surf(
-        coords,
-        faces,
-        overlay,
-        flat_map=flat_map,
-        rotate=[90, 270],
-        filename=tmp_file,
-        vmin=vmin,
-        vmax=vmax,
-    )
-    im = Image.open(tmp_file)
-    im = trim(im)
-    im = im.convert("RGBA")
-    im1 = np.array(im)
-    os.remove(tmp_file)
+    with tempfile.TemporaryDirectory() as tmpdir:
+        tmp_file = os.path.join(tmpdir,'tmp.png')
+        msp.plot_surf(
+            coords,
+            faces,
+            overlay,
+            flat_map=flat_map,
+            rotate=[90, 270],
+            filename=tmp_file,
+            vmin=vmin,
+            vmax=vmax,
+        )
+        im = Image.open(tmp_file)
+        im = trim(im)
+        im = im.convert("RGBA")
+        im1 = np.array(im)
+        os.remove(tmp_file)
     return im1
 
 
@@ -1195,7 +1188,7 @@ def plot_roc_multiple(roc_dictionary, roc_curves_thresholds):
 def load_prediction(subject, hdf5, dset="prediction"):
     """load network predictions"""
     results = {}
-    with h5py.File(hdf5, "r") as f:
+    with open_hdf5_file(hdf5, mode="r") as f:
         for hemi in ["lh", "rh"]:
             results[hemi] = f[subject][hemi][dset][:]
     return results

meld_graph/hdf5_utils.py

@@ -0,0 +1,57 @@
+import os
+import random
+import time
+
+import h5py
+
+
+# Modes that require the file to already exist. h5py raises FileNotFoundError
+# (a subclass of OSError) for these, which must not be retried.
+_READ_MODES = ("r", "r+")
+
+
+def open_hdf5_file(file_path, mode="r", max_retries=100, sleep_max_seconds=1.0, create_parent=False):
+    """Open an HDF5 file with bounded retries for transient OSErrors.
+
+    The retries exist to survive lock contention when several processes work on
+    the same file. A missing file is not transient, so it fails immediately
+    instead of being retried for the full timeout.
+
+    Args:
+        file_path: HDF5 path.
+        mode: Mode passed to h5py.File.
+        max_retries: Maximum number of open attempts.
+        sleep_max_seconds: Max random sleep between retries.
+        create_parent: If True, create parent directory before opening.
+
+    Returns:
+        Open h5py.File handle.
+
+    Raises:
+        FileNotFoundError: if mode requires an existing file and there is none.
+        OSError: if opening fails after all retries.
+    """
+    if max_retries < 1:
+        raise ValueError("max_retries must be >= 1")
+
+    if mode in _READ_MODES and not os.path.isfile(file_path):
+        raise FileNotFoundError(f"No such HDF5 file: {file_path} (mode={mode})")
+
+    if create_parent:
+        parent_dir = os.path.dirname(file_path)
+        if parent_dir:
+            os.makedirs(parent_dir, exist_ok=True)
+
+    last_error = None
+    for attempt in range(max_retries):
+        try:
+            return h5py.File(file_path, mode)
+        except OSError as err:
+            last_error = err
+            if attempt >= max_retries - 1:
+                break
+            time.sleep(random.random() * sleep_max_seconds)
+
+    raise OSError(
+        f"Could not open HDF5 file after {max_retries} retries: {file_path} (mode={mode})"
+    ) from last_error

meld_graph/meld_cohort.py

@@ -21,6 +21,7 @@
 import logging
 import meld_graph.mesh_tools as mt
 import scipy
+from meld_graph.hdf5_utils import open_hdf5_file
 
 
 class MeldCohort:
@@ -191,13 +192,13 @@ def _site_hdf5(self, site_code, group, write=False, hdf5_file_root=None):
             hdf5_file_root = self.hdf5_file_root
 
         p = os.path.join(self.data_dir, f"MELD_{site_code}", hdf5_file_root.format(site_code=site_code, group=group))
-        # open existing file or create new one
+
         if os.path.isfile(p) and not write:
-            f = h5py.File(p, "r")
+            f = open_hdf5_file(p, mode="r")
         elif os.path.isfile(p) and write:
-            f = h5py.File(p, "r+")
+            f = open_hdf5_file(p, mode="r+")
         elif not os.path.isfile(p) and write:
-            f = h5py.File(p, "a")
+            f = open_hdf5_file(p, mode="a", create_parent=True)
         else:
             f = None
         try:
@@ -648,10 +649,7 @@ def write_feature_values(self, feature, feature_values, hemis=["lh", "rh"], hdf5
         n_vert_cortex = sum(self.cohort.cortex_mask)
         # open hdf5 file
         if hdf5_file is not None:
-            if not os.path.isfile(hdf5_file):
-                hdf5_file_context = h5py.File(hdf5_file, "a")
-            else:
-                hdf5_file_context = h5py.File(hdf5_file, "r+")
+            hdf5_file_context = open_hdf5_file(hdf5_file, mode="a", create_parent=True)
         else:
             hdf5_file_context = self.cohort._site_hdf5(
                 self.site_code, self.group, write=True, hdf5_file_root=hdf5_file_root

meld_graph/mesh_tools.py

@@ -1,5 +1,6 @@
 #Tools for mesh-based operations
 #Smoothing, 
+import tempfile
 import nibabel as nb
 import numpy as np
 import subprocess
@@ -99,36 +100,34 @@ def smooth_array(input_array, neighbours, n_iter=70, cortex_mask=None):
 
 def save_mgh(filename, array, demo):
     """save mgh file using nibabel and imported demo mgh file"""
-    rand_int = np.random.randint(100)
-    mmap = np.memmap("/tmp/tmp" + str(rand_int), dtype="float32", mode="w+", shape=demo.get_data().shape)
-    mmap[:, 0, 0] = array[:]
-    output = nb.MGHImage(mmap, demo.affine, demo.header)
-    nb.save(output, filename)
+    with tempfile.NamedTemporaryFile() as mmap_file:
+        mmap = np.memmap(mmap_file.name, dtype="float32", mode="w+", shape=demo.get_data().shape)
+        mmap[:, 0, 0] = array[:]
+        output = nb.MGHImage(mmap, demo.affine, demo.header)
+        nb.save(output, filename)
 
 
 def smoothing_fs(overlay, fwhm, subject="fsaverage_sym", hemi="lh", subjects_dir=MELD_PARAMS_PATH):
     os.environ["SUBJECTS_DIR"] = subjects_dir
 
     """smooth surface overlay on fsaverage_sym"""
-    tmpdir = "/tmp/" + str(np.random.randint(1000000))
-    os.mkdir(tmpdir)
-    dum = nb.load(os.path.join(subjects_dir, subject, "surf", hemi + ".white.avg.area.mgh"))
-    save_mgh(os.path.join(tmpdir, hemi + ".tmp.mgh"), overlay, dum)
-    subprocess.call(
-        "mris_fwhm --s "
-        + subject
-        + " --hemi "
-        + hemi
-        + " --cortex --smooth-only --fwhm "
-        + str(fwhm)
-        + " --i "
-        + os.path.join(tmpdir, hemi + ".tmp.mgh")
-        + " --o "
-        + os.path.join(tmpdir, hemi + ".sm_tmp.mgh"),
-        shell=True,
-    )
-    overlay_smoothed = load_mgh(os.path.join(tmpdir, hemi + ".sm_tmp.mgh"))
-    subprocess.call("rm -r " + tmpdir, shell=True)
+    with tempfile.TemporaryDirectory() as tmpdir:
+        dum = nb.load(os.path.join(subjects_dir, subject, "surf", hemi + ".white.avg.area.mgh"))
+        save_mgh(os.path.join(tmpdir, hemi + ".tmp.mgh"), overlay, dum)
+        subprocess.call(
+            "mris_fwhm --s "
+            + subject
+            + " --hemi "
+            + hemi
+            + " --cortex --smooth-only --fwhm "
+            + str(fwhm)
+            + " --i "
+            + os.path.join(tmpdir, hemi + ".tmp.mgh")
+            + " --o "
+            + os.path.join(tmpdir, hemi + ".sm_tmp.mgh"),
+            shell=True,
+        )
+        overlay_smoothed = load_mgh(os.path.join(tmpdir, hemi + ".sm_tmp.mgh"))
     return overlay_smoothed