From 26a3680f56d45d79d9e3cf9976453885acb80f06 Mon Sep 17 00:00:00 2001 From: MahdiAll99 Date: Tue, 24 Feb 2026 08:54:05 -0500 Subject: [PATCH 01/15] Random comment deleted --- MEDiml/learning/FSR.py | 4 +--- MEDiml/wrangling/DataManager.py | 11 +++++++---- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/MEDiml/learning/FSR.py b/MEDiml/learning/FSR.py index a3631d8..5460790 100644 --- a/MEDiml/learning/FSR.py +++ b/MEDiml/learning/FSR.py @@ -60,9 +60,6 @@ def __get_fda_corr_table( # Keep only variables that are in both tables _, outcome_table_binary = intersect_var_tables(variable_table, outcome_table_binary) - # Under-sample the outcome table to equalize the number of positive and negative outcomes - #outcome_table_binary_balanced = under_sample(outcome_table_binary) - # Get the patient teach split patients_teach_splits = get_stratified_splits( outcome_table_binary, @@ -648,6 +645,7 @@ def apply_fsr(self, ml: Dict, variable_table: List, outcome_table_binary: pd.Dat ml (dict): Machine learning dictionary containing the learning options. variable_table (List): Table of variables. outcome_table_binary (pd.DataFrame): Table of binary outcomes. + path_save_logging (Path, optional): Path to save logging information. Defaults to None. Returns: List: Table of variables after feature set reduction. diff --git a/MEDiml/wrangling/DataManager.py b/MEDiml/wrangling/DataManager.py index a8f0035..efa6b13 100644 --- a/MEDiml/wrangling/DataManager.py +++ b/MEDiml/wrangling/DataManager.py @@ -456,7 +456,7 @@ def load_mask(_id, file, medscan): path_roi_data = self.paths._path_to_niftis for file in self.__nifti.stack_path_roi: - _id = image_file.name.split("(")[0] if ("(") in image_file.name else image_file.name # id is PatientID__ImagingScanName + _id = file.name.split("(")[0] if ("(") in file.name else file.name # id is PatientID__ImagingScanName load_mask(_id, file, medscan) roi_index += 1 else: @@ -831,11 +831,13 @@ def __pre_radiomics_checks_dimensions( for f in tqdm(range(len(file_paths))): try: if file_paths[f].name.endswith("nii.gz") or file_paths[f].name.endswith("nii"): - medscan = nib.load(file_paths[f]) + with open(file_paths[f], 'rb') as file: + medscan = pickle.load(file) xy_dim["data"][f] = medscan.header.get_zooms()[0] z_dim["data"][f] = medscan.header.get_zooms()[2] else: - medscan = np.load(file_paths[f], allow_pickle=True) + with open(file_paths[f], 'rb') as file: + medscan = pickle.load(file) xy_dim["data"][f] = medscan.data.volume.spatialRef.PixelExtentInWorldX z_dim["data"][f] = medscan.data.volume.spatialRef.PixelExtentInWorldZ except Exception as e: @@ -1007,7 +1009,8 @@ def __pre_radiomics_checks_window( if file.name.endswith('nii.gz') or file.name.endswith('nii'): medscan = self.__process_one_nifti(file, path_data) else: - medscan = np.load(file, allow_pickle=True) + with open(file, 'rb') as file: + medscan = pickle.load(file) if re.search('PTscan', wildcard) and medscan.format != 'nifti': medscan.data.volume.array = compute_suv_map( np.double(medscan.data.volume.array), From bba239a3ed7faeef21ae4de12a60e051d8984d5b Mon Sep 17 00:00:00 2001 From: MahdiAll99 Date: Tue, 24 Feb 2026 08:55:16 -0500 Subject: [PATCH 02/15] Rad names fix --- MEDiml/utils/get_full_rad_names.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/MEDiml/utils/get_full_rad_names.py b/MEDiml/utils/get_full_rad_names.py index 6d94634..c255e46 100644 --- a/MEDiml/utils/get_full_rad_names.py +++ b/MEDiml/utils/get_full_rad_names.py @@ -16,6 +16,6 @@ def get_full_rad_names(str_user_data: str, rad_var_ids: List): full_rad_names = np.array([]) for rad_var in rad_var_ids: ind_var = int(rad_var[6:]) - full_rad_names = np.append(full_rad_names, str_user_data.split('||')[ind_var].split(':')[1]) + full_rad_names = np.append(full_rad_names, str_user_data.split('||')[ind_var].split(f"{rad_var}:")[1]) return full_rad_names From c4a5812e963f2307b9bc35bab22b8bd7a5dbb25a Mon Sep 17 00:00:00 2001 From: MahdiAll99 Date: Tue, 24 Feb 2026 08:56:02 -0500 Subject: [PATCH 03/15] pydicom version downgraded --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index e3cbf27..69c61b5 100644 --- a/requirements.txt +++ b/requirements.txt @@ -13,7 +13,7 @@ pandas<2.0.0 Pillow protobuf pycaret -pydicom +pydicom<2.0.0 PyWavelets ray[default] scikit_image From 96e76147b7c4ccfed0c8b35a23f0a0c53c78d04b Mon Sep 17 00:00:00 2001 From: MahdiAll99 Date: Tue, 24 Feb 2026 09:00:03 -0500 Subject: [PATCH 04/15] badge fix --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 5ee6e70..920043b 100644 --- a/README.md +++ b/README.md @@ -3,7 +3,7 @@ [![PyPI - Python Version](https://img.shields.io/badge/python-3.8%20|%203.9%20|%203.10-blue)](https://www.python.org/downloads/release/python-380/) -[![PyPI - version](https://img.shields.io/badge/pypi-v0.9.10-blue)](https://pypi.org/project/medimage-pkg/) +[![PyPI - version](https://img.shields.io/badge/pypi-v0.9.11-blue)](https://pypi.org/project/mediml/) [![Continuous Integration](https://github.com/MEDomicsLab/MEDiml/actions/workflows/python-app.yml/badge.svg)](https://github.com/MEDomicsLab/MEDiml/actions/workflows/python-app.yml) [![Documentation Status](https://readthedocs.org/projects/mediml/badge/?version=latest)](https://mediml.readthedocs.io/en/latest/?badge=latest) [![License: GPL-3](https://img.shields.io/badge/license-GPLv3-blue)](LICENSE) From 8a17c771ae8d4be9454f2505a36ffe6bb2e796d5 Mon Sep 17 00:00:00 2001 From: MahdiAll99 Date: Tue, 24 Feb 2026 09:00:29 -0500 Subject: [PATCH 05/15] new version 0.9.11 --- MEDiml/__init__.py | 2 +- pyproject.toml | 2 +- setup.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/MEDiml/__init__.py b/MEDiml/__init__.py index 06e15db..7f75fc3 100644 --- a/MEDiml/__init__.py +++ b/MEDiml/__init__.py @@ -14,7 +14,7 @@ logging.getLogger(__name__).addHandler(stream_handler) __author__ = "MEDomicsLab consortium" -__version__ = "0.9.10" +__version__ = "0.9.11" __copyright__ = "Copyright (C) MEDomicsLab consortium" __license__ = "GNU General Public License 3.0" __maintainer__ = "MAHDI AIT LHAJ LOUTFI" diff --git a/pyproject.toml b/pyproject.toml index a889645..baa9417 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [tool.poetry] name = "mediml" -version = "0.9.10" +version = "0.9.11" description = "MEDiml is a Python package for processing and extracting features from medical images" authors = ["MEDomics Consortium "] license = "GPL-3.0" diff --git a/setup.py b/setup.py index 110a76a..dd103c9 100644 --- a/setup.py +++ b/setup.py @@ -14,7 +14,7 @@ setup( name="MEDiml", - version="0.9.10", + version="0.9.11", author="MEDomics consortium", author_email="medomics.info@gmail.com", description="Python Open-source package for medical images processing and radiomic features extraction", From ba73f3377b492cf6942ac2aeeb18c38aae43de39 Mon Sep 17 00:00:00 2001 From: MahdiAll99 Date: Wed, 25 Feb 2026 14:57:28 -0500 Subject: [PATCH 06/15] Added a safe check --- MEDiml/biomarkers/BatchExtractor.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/MEDiml/biomarkers/BatchExtractor.py b/MEDiml/biomarkers/BatchExtractor.py index c08f699..06aeeeb 100644 --- a/MEDiml/biomarkers/BatchExtractor.py +++ b/MEDiml/biomarkers/BatchExtractor.py @@ -593,6 +593,14 @@ def __batch_all_patients(self, im_params: Dict) -> None: # READING CSV EXPERIMENT TABLE tabel_roi = pd.read_csv(self._path_csv / ('roiNames_' + roi_type_label + '.csv')) + + # Check if all the requires columns are present + for col in ['PatientID', 'ImagingScanName', 'ImagingModality', 'ROIname']: + if col not in list(tabel_roi.columns): + raise ValueError(f'Missing column "{col}" in the ROI CSV file for roi type "{roi_type_label}". \ + Please check that the CSV file contains all the required columns: "PatientID", "ImagingScanName", " \ + "ImagingModality" and "ROIname".') + tabel_roi['under'] = '_' tabel_roi['dot'] = '.' tabel_roi['npy'] = '.npy' From f95ddb33e58d873bfecb75e2d6ff313867a08756 Mon Sep 17 00:00:00 2001 From: MahdiAll99 Date: Fri, 27 Feb 2026 11:36:46 -0500 Subject: [PATCH 07/15] removed MEDimage mentions in all notebooks --- notebooks/ibsi/ibsi2p1.ipynb | 2 +- .../tutorial/BatchExtractor-Tutorial.ipynb | 20 +++++------ notebooks/tutorial/DataManager-Tutorial.ipynb | 12 +++---- notebooks/tutorial/Learning-Tutorial.ipynb | 34 +++++++++---------- notebooks/tutorial/MEDscan-Tutorial.ipynb | 20 +++++------ 5 files changed, 44 insertions(+), 44 deletions(-) diff --git a/notebooks/ibsi/ibsi2p1.ipynb b/notebooks/ibsi/ibsi2p1.ipynb index 2512bcd..1cf018b 100644 --- a/notebooks/ibsi/ibsi2p1.ipynb +++ b/notebooks/ibsi/ibsi2p1.ipynb @@ -67,7 +67,7 @@ "source": [ "### Initialization\n", "\n", - "In this chapter and phase we will not need any class and we will use the *MEDimage* package modules directly, especially the module ``filter`` which offers methods for image filtering. \n", + "In this chapter and phase we will not need any class and we will use the *MEDiml* package modules directly, especially the module ``filter`` which offers methods for image filtering. \n", "\n", "#### ⚠️ DOWNLOADING THE DATA IS REQUIRED ⚠️\n", "\n", diff --git a/notebooks/tutorial/BatchExtractor-Tutorial.ipynb b/notebooks/tutorial/BatchExtractor-Tutorial.ipynb index 5b0cee3..dceac52 100644 --- a/notebooks/tutorial/BatchExtractor-Tutorial.ipynb +++ b/notebooks/tutorial/BatchExtractor-Tutorial.ipynb @@ -5,7 +5,7 @@ "id": "ecdbaa79", "metadata": {}, "source": [ - "# Batch extraction Tutorial − Radiomics batch extraction using MEDimage package\n", + "# Batch extraction Tutorial − Radiomics batch extraction using MEDiml package\n", "\n", "@Author : [MEDomics consortium](https://github.com/medomics/)\n", "\n", @@ -32,10 +32,10 @@ "\n", "Running this notebook requires running the [DataManager-tutorial notebook](https://colab.research.google.com/github/MEDomicsLab/MEDiml/blob/dev/notebooks/tutorial/DataManager-Tutorial.ipynb). We also recommend that you take a look at [MEDscan-Tutorial notebook](https://colab.research.google.com/github/MEDomicsLab/MEDiml/blob/dev/notebooks/tutorial/MEDscan-Tutorial.ipynb) as well.\n", "\n", - "This notebook is a tutorial of radiomics batch extraction using the *MEDimage* package and specifically the ``BatchExtractor`` class. For this task, the ``BatchExtractor`` class is the main object used to order scans and prepare batches and run processing and features extraction. The class extracts all type of family features and organizes the results in json files and csv tables.\n", + "This notebook is a tutorial of radiomics batch extraction using the *MEDiml* package and specifically the ``BatchExtractor`` class. For this task, the ``BatchExtractor`` class is the main object used to order scans and prepare batches and run processing and features extraction. The class extracts all type of family features and organizes the results in json files and csv tables.\n", "\n", "\n", - "In a nutshell, This tutorial will help you learn everything you need about batch extraction in the *MEDimage package*. We also advise you to read the [class documentation](https://mediml.readthedocs.io/en/documentation/biomarkers.html#module-MEDimage.biomarkers.BatchExtractor) before starting to test it." + "In a nutshell, This tutorial will help you learn everything you need about batch extraction in the *MEDiml package*. We also advise you to read the [class documentation](https://mediml.readthedocs.io/en/latest/tutorials.html#batchextractor) before starting to test it." ] }, { @@ -55,7 +55,7 @@ "source": [ "### DICOM data\n", "\n", - "In this tutorial we will use data from STS study (soft-tissue-sorcoma) processed by McGill institute, containing 204 scans with different scan types (PTscan, CTscan...). We assume that you have already processed these scans in the [DataManager-tutorial notebook]()." + "In this tutorial we will use data from STS study (soft-tissue-sorcoma) processed by McGill institute, containing 204 scans with different scan types (PTscan, CTscan...). We assume that you have already processed these scans in the [DataManager-tutorial notebook](https://colab.research.google.com/github/MEDomicsLab/MEDiml/blob/dev/notebooks/tutorial/DataManager-Tutorial.ipynb)." ] }, { @@ -68,7 +68,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "id": "af05f061", "metadata": {}, "outputs": [ @@ -84,7 +84,7 @@ "import os\n", "import sys\n", "\n", - "import MEDimage" + "import MEDiml" ] }, { @@ -124,7 +124,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "id": "5ff7be32", "metadata": {}, "outputs": [], @@ -134,7 +134,7 @@ "\n", "path_read = Path(os.getcwd()) / \"data\" / \"npy\"\n", "path_csv = Path(os.getcwd()) / \"CSV\"\n", - "path_to_params = Path(os.getcwd()) / \"settings\" / \"MEDimage-Tutorial.json\"\n", + "path_to_params = Path(os.getcwd()) / \"settings\" / \"MEDiml-Tutorial.json\"\n", "path_save = Path(os.getcwd())" ] }, @@ -148,12 +148,12 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": null, "id": "1afd1975", "metadata": {}, "outputs": [], "source": [ - "batch_extractor = MEDimage.biomarkers.BatchExtractor(\n", + "batch_extractor = MEDiml.biomarkers.BatchExtractor(\n", " path_read=path_read,\n", " path_csv=path_csv,\n", " path_params=path_to_params,\n", diff --git a/notebooks/tutorial/DataManager-Tutorial.ipynb b/notebooks/tutorial/DataManager-Tutorial.ipynb index 817e0fc..5162fdc 100644 --- a/notebooks/tutorial/DataManager-Tutorial.ipynb +++ b/notebooks/tutorial/DataManager-Tutorial.ipynb @@ -31,7 +31,7 @@ "## Introduction\n", "\n", "\n", - "This notebook is a tutorial for the *DataManager* class to give a detailed introduction & explanation on how to use this Python class. The *DataManager* class is the main object used in the *MEDimage* package when it comes to processing raw data in NIfTI and DICOM formats. This class can:\n", + "This notebook is a tutorial for the *DataManager* class to give a detailed introduction & explanation on how to use this Python class. The *DataManager* class is the main object used in the *MEDiml* package when it comes to processing raw data in NIfTI and DICOM formats. This class can:\n", " - Create ``MEDscan`` class objects from the raw data and makes the manipulation of these objects easy.\n", " - Help find the proper dimension and re-segmentation ranges options for radiomics analysis by running some pre-computation checks.\n", "\n", @@ -43,7 +43,7 @@ "id": "95b2b171", "metadata": {}, "source": [ - "The ``DataManager`` class is one of the first operation done in the radiomics analysis workflow, because it helps create the ``MEDscan`` class objects which is the main asset used in the *MEDimage* package.\n", + "The ``DataManager`` class is one of the first operation done in the radiomics analysis workflow, because it helps create the ``MEDscan`` class objects which is the main asset used in the *MEDiml* package.\n", "\n", "\n", "\n", @@ -82,10 +82,10 @@ "import os\n", "import sys\n", "\n", - "MODULE_DIR = os.path.dirname(os.path.abspath('../MEDimage/'))\n", + "MODULE_DIR = os.path.dirname(os.path.abspath('../MEDiml/'))\n", "sys.path.append(os.path.dirname(MODULE_DIR))\n", "\n", - "import MEDimage" + "import MEDiml" ] }, { @@ -566,7 +566,7 @@ "\n", "\n", "\n", - "We will go through all the functionalities of the ``DataManager`` class. For more details about the class please refer to the [DataManager documentation](https://mediml.readthedocs.io/en/documentation/wrangling.html#module-MEDimage.wrangling.DataManager)" + "We will go through all the functionalities of the ``DataManager`` class. For more details about the class please refer to the [DataManager documentation](https://mediml.readthedocs.io/en/documentation/wrangling.html#module-MEDiml.wrangling.DataManager)" ] }, { @@ -599,7 +599,7 @@ "path_dicoms_data = Path(os.getcwd()) / \"data\" / \"DICOM-STS\"\n", "path_save = Path(os.getcwd()) / \"data\" / \"npy\"\n", "path_save.mkdir() if not path_save.exists() else path_save\n", - "dm = MEDimage.wrangling.DataManager(path_to_dicoms=path_dicoms_data,\n", + "dm = MEDiml.wrangling.DataManager(path_to_dicoms=path_dicoms_data,\n", " path_save=path_save,\n", " path_csv=path_csv,\n", " n_batch=2)" diff --git a/notebooks/tutorial/Learning-Tutorial.ipynb b/notebooks/tutorial/Learning-Tutorial.ipynb index 615b3d6..303d123 100644 --- a/notebooks/tutorial/Learning-Tutorial.ipynb +++ b/notebooks/tutorial/Learning-Tutorial.ipynb @@ -58,7 +58,7 @@ "from numpyencoder import NumpyEncoder\n", "\n", "\n", - "import MEDimage" + "import MEDiml" ] }, { @@ -80,7 +80,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "Move the features folder to /home/local/USHERBROOKE/aitm2302/Bureau/MEDOMICS-ALL/MEDimagePackage/MEDimageUds/MEDimage/notebooks/tutorial/learning/experiments/IDH/holdOut__all_learn__001\n" + "Move the features folder to /home/local/USHERBROOKE/aitm2302/Bureau/MEDOMICS-ALL/MEDimagePackage/MEDimageUds/MEDiml/notebooks/tutorial/learning/experiments/IDH/holdOut__all_learn__001\n" ] } ], @@ -98,7 +98,7 @@ "path_save_experiment = path_study / 'experiments'\n", "\n", "# Seperate data (using the outcomes file)\n", - "path_data, _ = MEDimage.learning.ml_utils.create_holdout_set(\n", + "path_data, _ = MEDiml.learning.ml_utils.create_holdout_set(\n", " path_outcome_file=path_outcome_file, \n", " outcome_name='IDH', \n", " path_save_experiments=path_save_experiment,\n", @@ -18570,14 +18570,14 @@ "experiment_label = \"GliomaIDH_Image_All\"\n", "\n", "# Load json and adapt the settings to the test (make sure all modalities are included)\n", - "ml_variables = MEDimage.utils.json_utils.load_json(path_settings / \"ml_variables.json\")\n", + "ml_variables = MEDiml.utils.json_utils.load_json(path_settings / \"ml_variables.json\")\n", "ml_variables[\"var1\"][\"nameType\"] = \"RadiomicsFull\"\n", "ml_variables[\"var1\"][\"imSpaces\"] = [\"image\"] # image refers to the csv where all feature types are included\n", "ml_variables[\"var1\"][\"scans\"] = [\"T1\", \"T1CE\", \"T2\",\"T2F\"]\n", - "MEDimage.utils.json_utils.save_json(path_settings / \"ml_variables.json\", ml_variables, cls=NumpyEncoder)\n", + "MEDiml.utils.json_utils.save_json(path_settings / \"ml_variables.json\", ml_variables, cls=NumpyEncoder)\n", "\n", "# Initialize the radiomics learner class (Main machine learning class)\n", - "learner = MEDimage.learning.RadiomicsLearner(\n", + "learner = MEDiml.learning.RadiomicsLearner(\n", " path_study=path_data, \n", " path_settings=path_settings, \n", " experiment_label=experiment_label\n", @@ -18700,7 +18700,7 @@ } ], "source": [ - "MEDimage.utils.load_json(Path(path_data) / f'learn__{experiment_label}' /'results_avg.json')['test']" + "MEDiml.utils.load_json(Path(path_data) / f'learn__{experiment_label}' /'results_avg.json')['test']" ] }, { @@ -18733,11 +18733,11 @@ "outputs": [], "source": [ "# Load json and adapt the settings to the new test\n", - "ml_variables = MEDimage.utils.json_utils.load_json(path_settings / \"ml_variables.json\")\n", + "ml_variables = MEDiml.utils.json_utils.load_json(path_settings / \"ml_variables.json\")\n", "ml_variables[\"var1\"][\"nameType\"] = \"RadiomicsMorph\"\n", "ml_variables[\"var1\"][\"imSpaces\"] = [\"morph\"]\n", "ml_variables[\"var1\"][\"scans\"] = [\"T2F\"]\n", - "MEDimage.utils.json_utils.save_json(path_settings / \"ml_variables.json\", ml_variables, cls=NumpyEncoder)" + "MEDiml.utils.json_utils.save_json(path_settings / \"ml_variables.json\", ml_variables, cls=NumpyEncoder)" ] }, { @@ -38697,7 +38697,7 @@ "path_settings = Path.cwd() / \"learning\" / \"settings\"\n", "\n", "# Initialize the radiomics learner class (Main machine learning class)\n", - "learner = MEDimage.learning.RadiomicsLearner(\n", + "learner = MEDiml.learning.RadiomicsLearner(\n", " path_study=path_data, \n", " path_settings=path_settings, \n", " experiment_label=experiment_label\n", @@ -38725,10 +38725,10 @@ "outputs": [], "source": [ "# Load json and adapt the settings to the new test\n", - "ml_variables = MEDimage.utils.json_utils.load_json(path_settings / \"ml_variables.json\")\n", + "ml_variables = MEDiml.utils.json_utils.load_json(path_settings / \"ml_variables.json\")\n", "ml_variables[\"var1\"][\"nameType\"] = \"RadiomicsInt\"\n", "ml_variables[\"var1\"][\"imSpaces\"] = [\"intensity\"]\n", - "MEDimage.utils.json_utils.save_json(path_settings / \"ml_variables.json\", ml_variables, cls=NumpyEncoder)" + "MEDiml.utils.json_utils.save_json(path_settings / \"ml_variables.json\", ml_variables, cls=NumpyEncoder)" ] }, { @@ -59814,7 +59814,7 @@ "path_settings = Path.cwd() / \"learning\" / \"settings\"\n", "\n", "# Initialize the radiomics learner class (Main machine learning class)\n", - "learner = MEDimage.learning.RadiomicsLearner(\n", + "learner = MEDiml.learning.RadiomicsLearner(\n", " path_study=path_data, \n", " path_settings=path_settings, \n", " experiment_label=experiment_label\n", @@ -59842,10 +59842,10 @@ "outputs": [], "source": [ "# Load json and adapt the settings to the new test\n", - "ml_variables = MEDimage.utils.json_utils.load_json(path_settings / \"ml_variables.json\")\n", + "ml_variables = MEDiml.utils.json_utils.load_json(path_settings / \"ml_variables.json\")\n", "ml_variables[\"var1\"][\"nameType\"] = \"RadiomicsText\"\n", "ml_variables[\"var1\"][\"imSpaces\"] = [\"texture\"]\n", - "MEDimage.utils.json_utils.save_json(path_settings / \"ml_variables.json\", ml_variables, cls=NumpyEncoder)" + "MEDiml.utils.json_utils.save_json(path_settings / \"ml_variables.json\", ml_variables, cls=NumpyEncoder)" ] }, { @@ -80807,7 +80807,7 @@ "path_settings = Path.cwd() / \"learning\" / \"settings\"\n", "\n", "# Initialize the radiomics learner class (Main machine learning class)\n", - "learner = MEDimage.learning.RadiomicsLearner(\n", + "learner = MEDiml.learning.RadiomicsLearner(\n", " path_study=path_data, \n", " path_settings=path_settings, \n", " experiment_label=experiment_label\n", @@ -80847,7 +80847,7 @@ ], "source": [ "# First instantiate the results class\n", - "result = MEDimage.learning.Results()\n", + "result = MEDiml.learning.Results()\n", "\n", "# Feel free to change plot options (metrics, title, figure size...)\n", "result.plot_heatmap(\n", diff --git a/notebooks/tutorial/MEDscan-Tutorial.ipynb b/notebooks/tutorial/MEDscan-Tutorial.ipynb index 8210ebe..5bd3c2e 100644 --- a/notebooks/tutorial/MEDscan-Tutorial.ipynb +++ b/notebooks/tutorial/MEDscan-Tutorial.ipynb @@ -45,7 +45,7 @@ "\n", "We recommed you run the [DataManager-Tutorial](https://colab.research.google.com/github/MEDomicsLab/MEDiml/blob/dev/notebooks/tutorial/DataManager-Tutorial.ipynb) before going through this one.\n", "\n", - "This notebook is a tutorial for the ``MEDscan`` class to give a detailed introduction & explanation on how the Python class is created, used and saved. The ``MEDscan`` class is the main object used in the *MEDimage* package either when it comes to processing, features extraction or any other type of other image analysis. It contains many attributes, child classes and many methods that holds information about the imaging data we're using, the processing and computation parameters etc. This makes the *MEDimage* package an excellent tool for radiomics studies and the ``MEDscan`` class a main tool to the use of this package. \n", + "This notebook is a tutorial for the ``MEDscan`` class to give a detailed introduction & explanation on how the Python class is created, used and saved. The ``MEDscan`` class is the main object used in the *MEDiml* package either when it comes to processing, features extraction or any other type of other image analysis. It contains many attributes, child classes and many methods that holds information about the imaging data we're using, the processing and computation parameters etc. This makes the *MEDiml* package an excellent tool for radiomics studies and the ``MEDscan`` class a main tool to the use of this package. \n", "\n", "In a nutshell, This tutorial will help you learn everything you need about the ``MEDscan`` class.\n" ] @@ -59,7 +59,7 @@ "\n", "\n", "\n", - "So using the *MEDimage* package and class, we get the following flowchart\n", + "So using the *MEDiml* package and class, we get the following flowchart\n", "\n", "\n" ] @@ -74,7 +74,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "id": "af05f061", "metadata": {}, "outputs": [ @@ -90,7 +90,7 @@ "import os\n", "import sys\n", "\n", - "import MEDimage" + "import MEDiml" ] }, { @@ -143,7 +143,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "id": "e2bb8157", "metadata": {}, "outputs": [], @@ -151,7 +151,7 @@ "from pathlib import Path\n", "\n", "path_nifti_data = Path(os.getcwd()) / \"data\" / \"NIfTI\"\n", - "dm = MEDimage.wrangling.DataManager(path_to_niftis=path_nifti_data)" + "dm = MEDiml.wrangling.DataManager(path_to_niftis=path_nifti_data)" ] }, { @@ -287,7 +287,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": null, "id": "71a71edd", "metadata": {}, "outputs": [], @@ -295,7 +295,7 @@ "from pathlib import Path\n", "\n", "path_dicom_data = Path(os.getcwd()) / \"data\" / \"DICOM\"\n", - "dm = MEDimage.wrangling.DataManager(path_to_dicoms=path_dicom_data, path_save=path_dicom_data)" + "dm = MEDiml.wrangling.DataManager(path_to_dicoms=path_dicom_data, path_save=path_dicom_data)" ] }, { @@ -417,7 +417,7 @@ "\n", "\n", "\n", - "For more details about the class please refer to the [*MEDimage* documentation](https://mediml.readthedocs.io/en/latest/)" + "For more details about the class please refer to the [*MEDiml* documentation](https://mediml.readthedocs.io/en/latest/)" ] }, { @@ -942,7 +942,7 @@ "id": "93d9e8da", "metadata": {}, "source": [ - "You can update every class attribute value using the right class methods for that (check the class diagram above or the [*MEDimage* documentation](https://mediml.readthedocs.io/en/latest/))" + "You can update every class attribute value using the right class methods for that (check the class diagram above or the [*MEDiml* documentation](https://mediml.readthedocs.io/en/latest/))" ] } ], From 1e7a979a89fbc38ee1583652a7bdcfa073a64e37 Mon Sep 17 00:00:00 2001 From: MahdiAll99 Date: Tue, 3 Mar 2026 14:35:38 -0500 Subject: [PATCH 08/15] bug fix for ray --- MEDiml/biomarkers/BatchExtractor.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/MEDiml/biomarkers/BatchExtractor.py b/MEDiml/biomarkers/BatchExtractor.py index 06aeeeb..850b88e 100644 --- a/MEDiml/biomarkers/BatchExtractor.py +++ b/MEDiml/biomarkers/BatchExtractor.py @@ -804,7 +804,7 @@ def compute_radiomics(self, create_tables: bool = True) -> None: if ray.is_initialized(): ray.shutdown() - ray.init(local_mode=True, include_dashboard=True, num_cpus=self.n_bacth) + ray.init(local_mode=True, include_dashboard=False, num_cpus=self.n_bacth) # Batch all scans from CSV file and compute radiomics for each scan self.__batch_all_patients(im_params) From ef4c1ddc48a7f8662ebe87ff0e9986c2d03489d1 Mon Sep 17 00:00:00 2001 From: MahdiAll99 Date: Tue, 3 Mar 2026 14:36:07 -0500 Subject: [PATCH 09/15] very minor bug fix --- MEDiml/learning/DesignExperiment.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/MEDiml/learning/DesignExperiment.py b/MEDiml/learning/DesignExperiment.py index 428b56a..03f1b8f 100644 --- a/MEDiml/learning/DesignExperiment.py +++ b/MEDiml/learning/DesignExperiment.py @@ -200,7 +200,7 @@ def fill_learner_dict(self, path_ml_options: Path) -> Path: scans = var_struct['scans'] # list of imaging sequences rois = var_struct['rois'] # list of roi labels im_spaces = var_struct['imSpaces'] # list of image spaces (filterd and original) - use_combinations = var_struct['use_combinations'] # boolean to use combinations of scans and im_spaces + use_combinations = var_struct['use_combinations'] if 'use_combinations' in list(var_struct.keys()) else False # boolean to use combinations of scans and im_spaces if use_combinations: all_combinations = [] scans = list(var_struct['combinations'].keys()) From 8b192d55bbcc3c7e712ad65c872433e808ad8d68 Mon Sep 17 00:00:00 2001 From: MahdiAll99 Date: Tue, 3 Mar 2026 14:38:09 -0500 Subject: [PATCH 10/15] Huge commit: updated existing classes to 'sklearn' style classes and improved code --- MEDiml/learning/DataCleaner.py | 288 +++++++++------------ MEDiml/learning/Estimator.py | 36 +++ MEDiml/learning/Normalization.py | 194 +++++++------- MEDiml/learning/RadiomicsLearner.py | 383 +++++----------------------- MEDiml/learning/Results.py | 292 ++++++++++++--------- MEDiml/learning/__init__.py | 2 +- MEDiml/learning/ml_utils.py | 208 ++++++++------- MEDiml/utils/rf_learner.py | 137 ++++++++++ MEDiml/utils/xgboost_learner.py | 161 ++++++++++++ 9 files changed, 902 insertions(+), 799 deletions(-) create mode 100644 MEDiml/learning/Estimator.py create mode 100644 MEDiml/utils/rf_learner.py create mode 100644 MEDiml/utils/xgboost_learner.py diff --git a/MEDiml/learning/DataCleaner.py b/MEDiml/learning/DataCleaner.py index d3e2c92..a705c65 100644 --- a/MEDiml/learning/DataCleaner.py +++ b/MEDiml/learning/DataCleaner.py @@ -1,198 +1,152 @@ -import random -from typing import Dict, List - import numpy as np import pandas as pd +from sklearn.base import BaseEstimator, TransformerMixin +from sklearn.impute import SimpleImputer +from sklearn.utils import check_random_state -class DataCleaner: +class DataCleaner(BaseEstimator, TransformerMixin): """ - Class that will clean features of the csv by removing features with too many missing values, - too little variation, too many missing values per sample, too little variation per sample, - and imputing missing values. + A scikit-learn compatible transformer that cleans features by removing those + with too many missing values or too little variation, removes samples with + too many missing features, and imputes missing values. """ - def __init__(self, df_features: pd.DataFrame, type: str = "continuous"): - """ - Constructor of the class DataCleaner - - Args: - df_features (pd.DataFrame): Table of features. - type (str): Type of variable: "continuous", "hcategorical" or "icategorical". Defaults to "continuous". - """ - self.df_features = df_features - self.type = type - - def __update_df_features(self, var_of_type: List[str], flag_var_out: List[bool]) -> List[str]: + def __init__( + self, + var_type: str = "continuous", + imputation: str = "mean", + missingCutoffpf: float = 0.1, + missingCutoffps: float = 0.25, + covCutoff: float = 0.1, + random_state=None + ): """ - Updates the variable table by deleting the features that are not in the variable of type + Initializes the DataCleaner with specified parameters for feature and sample filtering and imputation. Args: - var_of_type (List[str]): List of variable names. - flag_var_out (List[bool]): List of variables to flag out. + var_type (str): Type of variable ("continuous", "hcategorical", "icategorical"). + imputation_method (str): Method of imputation ("mean", "median", "mode", "random"). + missing_cutoff_pf (float): Max % of missing values allowed per feature (column). + missing_cutoff_ps (float): Max % of missing values allowed per sample (row). + cov_cutoff (float): Min coefficient of variation allowed per feature. + random_state (int, RandomState instance or None): Seed for reproducibility. Returns: - List[str]: List of variable names that were not flagged out. + None """ - var_to_delete = np.delete(var_of_type, [i for i, v in enumerate(flag_var_out) if not v]) - var_of_type = np.delete(var_of_type, [i for i, v in enumerate(flag_var_out) if v]) - self.df_features = self.df_features.drop(var_to_delete, axis=1) - return var_of_type + self.var_type = var_type + self.imputation_method = imputation + self.missing_cutoff_pf = missingCutoffpf + self.missing_cutoff_ps = missingCutoffps + self.cov_cutoff = covCutoff + self.random_state = random_state + + # Attributes learned during fit + self.features_to_keep_ = None + self.imputer_ = None - def cut_off_missing_per_sample(self, var_of_type: List[str], missing_cutoff : float = 0.25) -> None: + def fit(self, X: pd.DataFrame, y: pd.DataFrame=None): """ - Removes observations/samples with more than ``missing_cutoff`` missing features. + Learns which features to keep based on missingness and variation thresholds. Args: - var_of_type (List[str]): List of variable names. - missing_cutoff (float): Maximum percentage cut-offs of missing features per sample. Defaults to 25%. + X (pd.DataFrame): Input feature data. + y (pd.DataFrame, optional): Ignored, present for API consistency by convention. Returns: - None. + DataCleaner: Returns self. """ - # Initialization - n_observation, n_features = self.df_features.shape - empty_vec = np.zeros(n_observation, dtype=int) - data = self.df_features[var_of_type] - empty_vec += data.isna().sum(axis=1).values + # Ensure input is a DataFrame + X = self._validate_input(X) - # Gathering results - ind_obs_out = np.where(((empty_vec/n_features) > missing_cutoff) == True) - self.df_features = self.df_features.drop(self.df_features.index[ind_obs_out]) - - def cut_off_missing_per_feature(self, var_of_type: List[str], missing_cutoff : float = 0.1) -> List[str]: - """ - Removes features with more than ``missing_cutoff`` missing patients. + # 1. Identify features to keep based on missingness (per feature) + missing_frac = X.isna().mean() + features_by_missing = missing_frac[missing_frac <= self.missing_cutoff_pf].index.tolist() - Args: - var_of_type (list): List of variable names. - missing_cutoff (float): maximal percentage cut-offs of missing patient samples per variable. + # 2. Identify features to keep based on Coefficient of Variation (CV) + # We calculate CV only on the features that passed the missingness check + X = X[features_by_missing] - Returns: - List[str]: List of variable names that were not flagged out. - """ - flag_var_out = (((self.df_features[var_of_type].isna().sum()) / self.df_features.shape[0]) > missing_cutoff) - return self.__update_df_features(var_of_type, flag_var_out) - - def cut_off_variation(self, var_of_type: List[str], cov_cutoff : float = 0.1) -> List[str]: - """ - Removes features with a coefficient of variation (cov) less than ``cov_cutoff``. - - Args: - var_of_type (list): List of variable names. - cov_cutoff (float): minimal coefficient of variation cut-offs over samples per variable. Defaults to 10%. - - Returns: - List[str]: List of variable names that were not flagged out. - """ + # Handle division by zero or near-zero means by adding epsilon eps = np.finfo(np.float32).eps - cov_df_features = (self.df_features[var_of_type].std(skipna=True) / self.df_features[var_of_type].mean(skipna=True)) - flag_var_out = cov_df_features.abs().add(eps) < cov_cutoff - return self.__update_df_features(var_of_type, flag_var_out) - - def impute_missing(self, var_of_type: List[str], imputation_method : str = "mean") -> None: - """ - Imputes missing values of the features of type. + std = X.std(skipna=True) + mean = X.mean(skipna=True).abs() + eps + cv = std / mean + + self.features_to_keep_ = cv[cv >= self.cov_cutoff].index.tolist() + X = X[self.features_to_keep_] - Args: - var_of_type (list): List of variable names. - imputation_method (str): Method of imputation. Can be "mean", "median", "mode" or "random". - For "random" imputation, a seed can be provided by adding the seed value after the method - name, for example "random42". + # 3. Fit the Imputer on the selected features + self._fit_imputer(X) - Returns: - None. + return self + + def transform(self, X: pd.DataFrame): """ - if self.type in ['continuous', 'hcategorical']: - # random imputation - if 'random' in imputation_method: - if len(imputation_method) > 6: - try: - seed = int(imputation_method[7:]) - random.seed(seed) - except Exception as e: - print(f"Warning: Seed must be an integer. Random seed will be set to None. str({e})") - random.seed(a=None) - else: - random.seed(a=None) - self.df_features[var_of_type] = self.df_features[var_of_type].apply(lambda x: x.fillna(random.choice(list(x.dropna(axis=0))))) - - # Imputation with median - elif 'median' in imputation_method: - self.df_features[var_of_type] = self.df_features[var_of_type].fillna(self.df_features[var_of_type].median()) - - # Imputation with mean - elif 'mean' in imputation_method: - self.df_features[var_of_type] = self.df_features[var_of_type].fillna(self.df_features[var_of_type].mean()) + Applies feature selection, sample filtering, and imputation. + """ + # check is fitted + if self.features_to_keep_ is None: + raise RuntimeError("You must fit the transformer before transforming data.") - else: - raise ValueError("Imputation method for continuous and hcategorical features must be 'random', 'median' or 'mean'.") + X = self._validate_input(X) - elif self.type in ['icategorical']: - if 'random' in imputation_method: - if len(imputation_method) > 6: - seed = int(imputation_method[7:]) - random.seed(seed) - else: - random.seed(a=None) - - self.df_features[var_of_type] = self.df_features[var_of_type].apply(lambda x: x.fillna(random.choice(list(x.dropna(axis=0))))) - - if 'mode' in imputation_method: - self.df_features[var_of_type] = self.df_features[var_of_type].fillna(self.df_features[var_of_type].mode().max()) - else: - raise ValueError("Variable type must be 'continuous', 'hcategorical' or 'icategorical'.") + # 1. Filter Features (Columns) + # Only keep columns learned during fit + X_transformed = X[self.features_to_keep_].copy() - def __call__(self, cleaning_dict: Dict, imputation_method: str = "mean", - missing_cutoff_ps: float = 0.25, missing_cutoff_pf: float = 0.1, - cov_cutoff:float = 0.1) -> pd.DataFrame: - """ - Applies data cleaning to the features of type. - - Args: - cleaning_dict (dict): Dictionary of cleaning parameters (missing cutoffs and coefficient of variation cutoffs etc.). - var_of_type (list, optional): List of variable names. - imputation_method (str): Method of imputation. Can be "mean", "median", "mode" or "random". - For "random" imputation, a seed can be provided by adding the seed value after the method - name, for example "random42". - missing_cutoff_ps (float, optional): maximal percentage cut-offs of missing features per sample. - missing_cutoff_pf (float, optional): maximal percentage cut-offs of missing samples per variable. - cov_cutoff (float, optional): minimal coefficient of variation cut-offs over samples per variable. + # 2. Filter Samples (Rows) based on missingness + missing_frac_rows = X_transformed.isna().mean(axis=1) + mask_rows_keep = missing_frac_rows <= self.missing_cutoff_ps + X_transformed = X_transformed.loc[mask_rows_keep] - Returns: - pd.DataFrame: Cleaned table of features. - """ - - # Initialization - var_of_type = self.df_features.Properties['userData']['variables']['continuous'] - - # Retrieve thresholds from cleaning_dict if not None - if cleaning_dict is not None: - missing_cutoff_pf = cleaning_dict['missingCutoffpf'] - missing_cutoff_ps = cleaning_dict['missingCutoffps'] - cov_cutoff = cleaning_dict['covCutoff'] - imputation_method = cleaning_dict['imputation'] + # 3. Impute Missing Values + X_imputed = self._apply_imputation(X_transformed) - # Replace infinite values with NaNs - self.df_features = self.df_features.replace([np.inf, -np.inf], np.nan) - - # Remove features with more than missing_cutoff_pf missing samples (NaNs) - var_of_type = self.cut_off_missing_per_feature(var_of_type, missing_cutoff_pf) - - # Check - if len(var_of_type) == 0: - return None - - # Remove features with a coefficient of variation less than cov_cutoff - var_of_type = self.cut_off_variation(var_of_type, cov_cutoff) - - # Check - if len(var_of_type) == 0: - return None - - # Remove scans with more than missing_cutoff_ps missing features - self.cut_off_missing_per_sample(var_of_type, missing_cutoff_ps) - - # Impute missing values - self.impute_missing(var_of_type, imputation_method) - - return self.df_features + # Return as DataFrame to maintain column names + return pd.DataFrame(X_imputed, columns=self.features_to_keep_, index=X_transformed.index) + + def _fit_imputer(self, X): + """Helper to initialize and fit the correct imputer logic.""" + # Handle 'random' manually as SimpleImputer doesn't support it + if "random" in self.imputation_method: + self.imputer_ = "random" # Marker logic + return + + # Map methods to SimpleImputer strategies + strategy_map = { + "mean": "mean", + "median": "median", + "mode": "most_frequent" + } + + # Default logic for icategorical (mode) vs continuous (mean/median) + if self.imputation_method not in strategy_map: + # Fallback logic from original class + if self.var_type == "icategorical": + strategy = "most_frequent" + else: + strategy = "mean" + else: + strategy = strategy_map[self.imputation_method] + + self.imputer_ = SimpleImputer(strategy=strategy) + self.imputer_.fit(X) + + def _apply_imputation(self, X): + """Helper to apply the imputation.""" + if self.imputer_ == "random": + rng = check_random_state(self.random_state) + # Custom random imputation logic: fill NaNs with random choice from valid values in that column + return X.apply(lambda col: col.fillna( + np.random.choice(col.dropna().values) if not col.dropna().empty else col.mean() # Fallback if empty + )) + else: + return self.imputer_.transform(X) + + def _validate_input(self, X): + """Ensures X is a DataFrame and handles infinite values.""" + if not isinstance(X, pd.DataFrame): + X = pd.DataFrame(X) + # Replace infs with NaNs (as per original class) + return X.replace([np.inf, -np.inf], np.nan) diff --git a/MEDiml/learning/Estimator.py b/MEDiml/learning/Estimator.py new file mode 100644 index 0000000..ba31b65 --- /dev/null +++ b/MEDiml/learning/Estimator.py @@ -0,0 +1,36 @@ +from sklearn.base import BaseEstimator, ClassifierMixin + +from ..utils.rf_learner import RandomForestEstimator +from ..utils.xgboost_learner import XGBoostEstimator + + +class Estimator(BaseEstimator, ClassifierMixin): + def __init__(self, algorithm: str, ml_config: dict): + self.ml_config = ml_config + self.algorithm = algorithm + self.estimator_ = None + + def _initialize_estimator(self): + """Factory algorithm to select the right internal class.""" + if self.algorithm == 'xgboost': + return XGBoostEstimator(**self.ml_config) + elif self.algorithm == 'rf': + return RandomForestEstimator(**self.ml_config) + else: + raise ValueError(f"Method {self.algorithm} not supported.") + + def fit(self, X, y): + self.estimator_ = self._initialize_estimator() + self.estimator_.fit(X, y) + self.classes_ = self.estimator_.classes_ + return self + + def predict(self, X): + return self.estimator_.predict(X) + + def predict_proba(self, X): + return self.estimator_.predict_proba(X) + + def save(self, filepath): + import joblib + joblib.dump(self, filepath) diff --git a/MEDiml/learning/Normalization.py b/MEDiml/learning/Normalization.py index ea1dfc2..88cd918 100644 --- a/MEDiml/learning/Normalization.py +++ b/MEDiml/learning/Normalization.py @@ -1,112 +1,122 @@ import numpy as np import pandas as pd from neuroCombat import neuroCombat +from sklearn.base import BaseEstimator, TransformerMixin from ..utils.get_institutions_from_ids import get_institutions_from_ids -class Normalization: +class CombatNormalization(BaseEstimator, TransformerMixin): + """ + Sklearn-compatible Transformer for ComBat Normalization. + + This transformer assumes the input X (DataFrame) contains both the features + to be normalized and the column identifying the institution/batch. + """ def __init__( - self, - method: str = 'combat', - variable_table: pd.DataFrame = None, - covariates_df: pd.DataFrame = None, - institutions: list = None - ) -> None: - """ - Constructor of the Normalization class. + self, + institution_col: str = None, + covariates: list = None, + drop_institution: bool = True + ): """ - self.method = method - self.variable_table = variable_table - self.covariates_df = covariates_df - self.institutions = institutions - - def apply_combat( - self, - variable_table: pd.DataFrame, - covariate_df: pd.DataFrame = None, - institutions: list = None - ) -> pd.DataFrame: + Args: + institution_col (str): Name of the column in X containing the institution/batch IDs. + If None, tries to derive from Index using util. + covariates (list): List of column names in X to treat as covariates (biological retention). + drop_institution (bool): If True, removes the institution column from output. """ - Applys ComBat Normalization method to the data. - More details :ref:`this link `. + self.institution_col = institution_col + self.covariates = covariates if covariates is not None else [] + self.drop_institution = drop_institution - Args: - variable_table (pd.DataFrame): pandas data frame on which Combat harmonization will be applied. - This table is of size N X F (Observations X Features) and has the IDs as index. - Requirements for this table + def fit(self, X, y=None): + """ + ComBat calculates parameters on the current batch data provided in transform. + Standard fit does nothing but validate input exists. + """ + return self - - Does not contain NaNs. - - No feature has 0 variance. - - All variables are continuous (For example: Radiomics variables). - covariate_df (pd.DataFrame, optional): N X M pandas data frame, where N must equal the number of - observations in variable_table. M is the number of covariates to include in the algorithm. - institutions (list, optional): List of size n_observations X 1 with the different institutions. - - Returns: - pd.DataFrame: variable_table after Combat harmonization. + def transform(self, X): """ - # Initializing the class attributes from the arguments - if variable_table is None: - if self.variable_table is None: - raise ValueError('variable_table must be given.') - else: - self.variable_table = variable_table - if covariate_df is not None: - self.covariates_df = covariate_df - if institutions: - self.institutions = institutions - - # Intializing the institutions if not given - if self.institutions is None: - patient_ids = pd.Series(self.variable_table.index) - self.institutions = get_institutions_from_ids(patient_ids) - all_institutions = self.institutions.unique() - for n in range(all_institutions.size): - self.institutions[self.institutions == all_institutions[n]] = n+1 - self.institutions = self.institutions.to_numpy(dtype=int) - self.institutions = np.reshape(self.institutions, (-1, 1)) + Applies ComBat Normalization. + """ + # Validate Input + if not isinstance(X, pd.DataFrame): + raise ValueError("Input X must be a pandas DataFrame.") - # No harmonization will be applied if there is only one institution - if np.unique(self.institutions).size < 2: - return self.variable_table + # Avoid modifying the original input + X_df = X.copy() - # Initializing the covariates if not given - if self.covariates_df is not None: - self.covariates_df['institution'] = self.institutions + # 1. Identify Institutions + if self.institution_col and self.institution_col in X_df.columns: + institutions = X_df[self.institution_col] + # If we plan to drop it later, we don't include it in features matrix + if self.drop_institution: + X_df = X_df.drop(columns=[self.institution_col]) else: - # the covars matrix is only a row with the institution - self.covariates_df = pd.DataFrame( - self.institutions, - columns=['institution'], - index=self.variable_table.index.values - ) + # Fallback to index-based logic from original code + institutions = get_institutions_from_ids(pd.Series(X_df.index)) - # Apply combat - n_features = self.variable_table.shape[1] - batch_col = 'institution' - if n_features == 1: - # combat does not work with a single feature so a temporary one is added, - # then removed later (this has no effect on the algorithm). - self.variable_table['temp'] = pd.Series( - np.ones(self.variable_table.shape[0]), - index=self.variable_table.index - ) - data_combat = neuroCombat( - self.variable_table.transpose(), - self.covariates_df, - batch_col - ) - self.variable_table = pd.DataFrame(self.variable_table.drop('temp', axis=1)) - vt_combat = pd.DataFrame(data_combat[:][0].transpose()) - else: - data_combat = neuroCombat( - self.variable_table.transpose(), - self.covariates_df, - batch_col + # Encode institutions to integers (1, 2, 3...) required by logic + institutions = self._process_institutions(institutions) + + # Check: If < 2 institutions, ComBat fails. Return original. + if len(np.unique(institutions)) < 2: + print("Warning: Less than 2 institutions detected. Skipping ComBat.") + return X_df + + # 2. Prepare Covariates + # We need to extract covariate data from X if specified + covars_df = pd.DataFrame({'institution': institutions.flatten()}, index=X_df.index) + + for cov in self.covariates: + if cov in X_df.columns: + covars_df[cov] = X_df[cov] + # Remove covariates from the feature matrix to be harmonized? + # Usually ComBat harmonizes features *adjusting* for covariates. + # We keep covariates in X_df usually, but standard combat expects data + # to ONLY be the features. + # Let's separate them: + X_df = X_df.drop(columns=[cov]) + + # 3. Handle Single Feature Edge Case (from original code) + cols_to_restore = [] + if X_df.shape[1] == 1: + X_df['temp_ones'] = 1 + cols_to_restore.append('temp_ones') + + # 4. Run NeuroCombat + # neuroCombat expects: data (Features x Samples), covars (Samples x Covars) + # Note: We transpose X_df because neuroCombat expects Features as Rows + try: + results = neuroCombat( + dat=X_df.T, + covars=covars_df, + batch_col='institution' ) - vt_combat = pd.DataFrame(data_combat['data']).transpose() + + # 5. Reconstruct DataFrame + # Result 'data' is Features x Samples + harmonized_data = pd.DataFrame(results['data'].T, index=X_df.index, columns=X_df.columns) + + # Remove temp columns if any + if cols_to_restore: + harmonized_data = harmonized_data.drop(columns=cols_to_restore) + + # Re-attach Covariates if they were stripped + for cov in self.covariates: + harmonized_data[cov] = covars_df[cov] + + return harmonized_data - self.variable_table[:] = vt_combat.values + except Exception as e: + print(f"ComBat failed: {e}. Returning original data.") + return X - return self.variable_table + def _process_institutions(self, institutions): + """Helper to map institution strings to integers.""" + institutions = pd.Series(institutions) + unique_inst = institutions.unique() + mapping = {inst: i+1 for i, inst in enumerate(unique_inst)} + return institutions.map(mapping).values.reshape(-1, 1) diff --git a/MEDiml/learning/RadiomicsLearner.py b/MEDiml/learning/RadiomicsLearner.py index a731e57..bcfb58b 100644 --- a/MEDiml/learning/RadiomicsLearner.py +++ b/MEDiml/learning/RadiomicsLearner.py @@ -5,23 +5,19 @@ from pathlib import Path from typing import Dict, List, Tuple -import numpy as np import pandas as pd from numpyencoder import NumpyEncoder from pycaret.classification import * -from sklearn import metrics -from sklearn.model_selection import GridSearchCV, RandomizedSearchCV -from xgboost import XGBClassifier from MEDiml.learning.DataCleaner import DataCleaner from MEDiml.learning.DesignExperiment import DesignExperiment +from MEDiml.learning.Estimator import Estimator from MEDiml.learning.FSR import FSR from MEDiml.learning.ml_utils import (average_results, combine_rad_tables, - feature_imporance_analysis, - finalize_rad_table, get_ml_test_table, - get_radiomics_table, intersect, - intersect_var_tables, save_model) -from MEDiml.learning.Normalization import Normalization + feature_importance_analysis, + get_ml_test_table, get_radiomics_table, + intersect) +from MEDiml.learning.Normalization import CombatNormalization from MEDiml.learning.Results import Results from ..utils.json_utils import load_json, save_json @@ -76,41 +72,6 @@ def __load_ml_info(self, ml_dict_paths: Dict) -> Dict: ml_dict['path_results'] = ml_dict_paths['results'] return ml_dict - - def __find_balanced_threshold( - self, - model: XGBClassifier, - variable_table: pd.DataFrame, - outcome_table_binary: pd.DataFrame - ) -> float: - """ - Finds the balanced threshold for the given machine learning test. - - Args: - model (XGBClassifier): Trained XGBoost classifier for the given machine learning run. - variable_table (pd.DataFrame): Radiomics table. - outcome_table_binary (pd.DataFrame): Outcome table with binary labels. - - Returns: - float: Balanced threshold for the given machine learning test. - """ - # Check is there is a feature mismatch - if model.feature_names_in_.shape[0] != variable_table.columns.shape[0]: - variable_table = variable_table.loc[:, model.feature_names_in_] - - # Getting the probability responses for each patient - prob_xgb = np.zeros((variable_table.index.shape[0], 1)) * np.nan - patient_ids = list(variable_table.index.values) - for p in range(variable_table.index.shape[0]): - prob_xgb[p] = self.predict_xgb(model, variable_table.loc[[patient_ids[p]], :]) - - # Calculating the ROC curve - fpr, tpr, thresholds = metrics.roc_curve(outcome_table_binary.iloc[:, 0], prob_xgb) - - # Calculating the optimal threshold by minizing fpr (false positive rate) and maximizing tpr (true positive rate) - minimum = np.argmin(np.power(fpr, 2) + np.power(1-tpr, 2)) - - return thresholds[minimum] def get_hold_out_set_table(self, ml: Dict, var_id: str, patients_id: List): """ @@ -219,9 +180,18 @@ def pre_process_radiomics_table( # Data cleaning if flags_preprocessing['var_datacleaning']: - cleaning_dict = ml['datacleaning'][ml['variables'][var_id]['var_datacleaning']]['feature']['continuous'] - data_cleaner = DataCleaner(rad_table_learning) - rad_table_learning = data_cleaner(cleaning_dict) + cleaning_dict = ml['datacleaning'][ml['variables'][var_id]['var_datacleaning']]['continuous'] + data_cleaner = DataCleaner(**cleaning_dict) + + # Temp save of properties + temp_properties = deepcopy(rad_table_learning.Properties) + + # Apply data cleaning + rad_table_learning = data_cleaner.fit_transform(rad_table_learning) + + # Re-assign properties + rad_table_learning.Properties = temp_properties + if rad_table_learning is None: continue @@ -242,8 +212,8 @@ def pre_process_radiomics_table( rad_table_learning.Properties['userData']['normalization']['original_data']['datacleaning_method'] = data_cln_method # Apply ComBat - normalization = Normalization('combat') - rad_table_learning = normalization.apply_combat(variable_table=rad_table_learning) # Training data + normalization = CombatNormalization() + rad_table_learning = normalization.fit_transform(rad_table_learning) # Training data else: raise NotImplementedError(f'Normalization method: {normalization_method} not recognized.') @@ -287,223 +257,6 @@ def pre_process_radiomics_table( return rad_tables_training, rad_tables_testing - def train_xgboost_model( - self, - var_table_train: pd.DataFrame, - outcome_table_binary_train: pd.DataFrame, - var_importance_threshold: float = 0.05, - optimal_threshold: float = None, - optimization_metric: str = 'MCC', - method : str = "pycaret", - use_gpu: bool = False, - seed: int = None, - ) -> Dict: - """ - Trains an XGBoost model for the given machine learning test. - - Args: - var_table_train (pd.DataFrame): Radiomics table for the training/learning set. - outcome_table_binary_train (pd.DataFrame): Outcome table with binary labels for the training/learning set. - var_importance_threshold (float): Threshold for the variable importance. Variables with importance below - this threshold will be removed from the model. - optimal_threshold (float, optional): Optimal threshold for the XGBoost model. If not given, it will be - computed using the training set. - optimization_metric (str, optional): String specifying the metric to use to optimize the ml model. - method (str, optional): String specifying the method to use to train the XGBoost model. - - "pycaret": Use PyCaret to train the model (automatic). - - "grid_search": Grid search with cross-validation to find the best parameters. - - "random_search": Random search with cross-validation to find the best parameters. - use_gpu (bool, optional): Boolean specifying if the GPU should be used to train the model. Default is True. - seed (int, optional): Integer specifying the seed to use for the random number generator. - - Returns: - Dict: Dictionary containing info about the trained XGBoost model. - """ - - # Safety check (make sure that the outcome table and the variable table have the same patients) - var_table_train, outcome_table_binary_train = intersect_var_tables(var_table_train, outcome_table_binary_train) - - # Finalize the new radiomics table with the remaining variables - var_table_train = finalize_rad_table(var_table_train) - - if method.lower() == "pycaret": - # Set up data for PyCaret - temp_data = pd.merge(var_table_train, outcome_table_binary_train, left_index=True, right_index=True) - - # PyCaret setup - setup( - data=temp_data, - feature_selection=True, - n_features_to_select=1-var_importance_threshold, - fold=5, - target=temp_data.columns[-1], - use_gpu=use_gpu, - feature_selection_estimator="xgboost", - session_id=seed - ) - - # Set seed - if seed is not None: - set_config('seed', seed) - - # Creating XGBoost model using PyCaret - classifier = create_model('xgboost', verbose=False) - - # Tuning XGBoost model using PyCaret - classifier = tune_model(classifier, optimize=optimization_metric) - - else: - # Initial training to filter features using variable importance - # XGB Classifier - classifier = XGBClassifier() - classifier.fit(var_table_train, outcome_table_binary_train) - var_importance = classifier.feature_importances_ - - # Normalize var_importance if necessary - if np.sum(var_importance) != 1: - var_importance_threshold = var_importance_threshold / np.sum(var_importance) - var_importance = var_importance / np.sum(var_importance) - - # Filter variables - var_table_train = var_table_train.iloc[:, var_importance >= var_importance_threshold] - - # Check if variable table is empty after filtering - if var_table_train.shape[1] == 0: - raise ValueError('Variable table is empty after variable importance filtering. Use a smaller threshold.') - - # Suggested scale_pos_weight - scale_pos_weight = 1 - (outcome_table_binary_train == 0).sum().values[0] \ - / (outcome_table_binary_train == 1).sum().values[0] - - # XGB Classifier - classifier = XGBClassifier(scale_pos_weight=scale_pos_weight) - - # Tune XGBoost parameters - params = { - 'max_depth': [3, 4, 5], - 'learning_rate': [0.1 , 0.01, 0.001], - 'n_estimators': [50, 100, 200] - } - - if method.lower() == "grid_search": - # Set up grid search with cross-validation - grid_search = GridSearchCV( - estimator=classifier, - param_grid=params, - cv=5, - n_jobs=-1, - verbose=3, - scoring='matthews_corrcoef' - ) - elif method.lower() == "random_search": - # Set up random search with cross-validation - grid_search = RandomizedSearchCV( - estimator=classifier, - param_distributions=params, - cv=5, - n_jobs=-1, - verbose=3, - scoring='matthews_corrcoef' - ) - else: - raise NotImplementedError(f'Method: {method} not recognized. Use "grid_search", "random_search", "auto" or "pycaret".') - - # Fit the grid search - grid_search.fit(var_table_train, outcome_table_binary_train) - - # Get the best parameters - best_params = grid_search.best_params_ - - # Fit the XGB Classifier with the best parameters - classifier = XGBClassifier(**best_params) - classifier.fit(var_table_train, outcome_table_binary_train) - - # Saving the information of the model in a dictionary - model_xgb = dict() - model_xgb['algo'] = 'xgb' - model_xgb['type'] = 'binary' - model_xgb['method'] = method - if optimal_threshold: - model_xgb['threshold'] = optimal_threshold - else: - try: - model_xgb['threshold'] = self.__find_balanced_threshold(classifier, var_table_train, outcome_table_binary_train) - except Exception as e: - print('Error in finding optimal threshold, it will be set to 0.5:' + str(e)) - model_xgb['threshold'] = 0.5 - model_xgb['model'] = classifier - model_xgb['var_names'] = list(classifier.feature_names_in_) - model_xgb['var_info'] = deepcopy(var_table_train.Properties['userData']) - if method == "auto": - model_xgb['optimization'] = "auto" - elif method == "pycaret": - model_xgb['optimization'] = classifier.get_params() - else: - model_xgb['optimization'] = best_params - - return model_xgb - - def test_xgb_model(self, model_dict: Dict, variable_table: pd.DataFrame, patient_list: List) -> List: - """ - Tests the XGBoost model for the given dataset patients. - - Args: - model_dict (Dict): Dictionary containing info about the trained XGBoost model. - variable_table (pd.DataFrame): Radiomics table for the test set (should not be normalized). - patient_list (List): List of patients to test. - - Returns: - List: List the model response for the training and test sets. - """ - # Initialization - n_test = len(patient_list) - var_names = model_dict['var_names'] - var_def = model_dict['var_info']['variables']['var_def'] - model_response = list() - - # Preparing the variable table - variable_table = get_ml_test_table(variable_table, var_names, var_def) - - # Test the model - for i in range(n_test): - # Get the patient IDs - patient_ids = patient_list[i] - - # Getting predictions for each patient - n_patients = len(patient_ids) - varargout = np.zeros((n_patients, 1)) * np.nan # NaN if the computation fails - for p in range(n_patients): - try: - varargout[p] = self.predict_xgb(model_dict['model'], variable_table.loc[[patient_ids[p]], :]) - except Exception as e: - print('Error in computing prediction for patient ' + str(patient_ids[p]) + ': ' + str(e)) - varargout[p] = np.nan - - # Save the predictions - model_response.append(varargout) - - return model_response - - def predict_xgb(self, xgb_model: XGBClassifier, variable_table: pd.DataFrame) -> float: - """ - Computes the prediction of the XGBoost model for the given variable table. - - Args: - xgb_model (XGBClassifier): XGBClassifier model. - variable_table (pd.DataFrame): Variable table for the prediction. - - Returns: - float: Prediction of the XGBoost model. - """ - - # Predictions - predictions = xgb_model.predict_proba(variable_table) - - # Get the probability of the positive class - predictions = predictions[:, 1][0] - - return predictions - def ml_run(self, path_ml: Path, holdout_test: bool = True, method: str = 'auto') -> None: """ This function runs the machine learning test for the ceated experiment. @@ -570,48 +323,50 @@ def ml_run(self, path_ml: Path, holdout_test: bool = True, method: str = 'auto') # Serperate variable table for training sets (repetitive but double-checking) var_table_train = processed_training_table.loc[patients_train, :] - # Initializing XGBoost model settings - var_importance_threshold = ml['algorithms']['XGBoost']['varImportanceThreshold'] - optimal_threshold = ml['algorithms']['XGBoost']['optimalThreshold'] - optimization_metric = ml['algorithms']['XGBoost']['optimizationMetric'] - method = ml['algorithms']['XGBoost']['method'] if 'method' in ml['algorithms']['XGBoost'].keys() else method - use_gpu = ml['algorithms']['XGBoost']['useGPU'] if 'useGPU' in ml['algorithms']['XGBoost'].keys() else True - seed = ml['algorithms']['XGBoost']['seed'] if 'seed' in ml['algorithms']['XGBoost'].keys() else None + # Initializing the model settings + algorithm = ml['settings']['algorithm'] + var_importance_threshold = ml['algorithms'][algorithm]['varImportanceThreshold'] + optimal_threshold = ml['algorithms'][algorithm]['optimalThreshold'] + optimization_metric = ml['algorithms'][algorithm]['optimizationMetric'] + method = ml['algorithms'][algorithm]['method'] if 'method' in ml['algorithms'][algorithm].keys() else method + use_gpu = ml['algorithms'][algorithm]['useGPU'] if 'useGPU' in ml['algorithms'][algorithm].keys() else True + seed = ml['algorithms'][algorithm]['seed'] if 'seed' in ml['algorithms'][algorithm].keys() else None - # B.2. Training the XGBoost model + # B.2. Training the model tstart = time.time() - logging.info(f"\n\n--> TRAINING XGBOOST MODEL FOR VARIABLE {var_id}") + logging.info(f"\n\n--> TRAINING {algorithm.upper()} MODEL FOR VARIABLE {var_id}") # Training the model - model = self.train_xgboost_model( - var_table_train, - outcome_table_binary_train, - var_importance_threshold, - optimal_threshold, - method=method, - use_gpu=use_gpu, - optimization_metric=optimization_metric, - seed=seed - ) + estimator = Estimator( + algorithm=algorithm, + ml_config={ + 'var_importance_threshold': var_importance_threshold, + 'optimal_threshold': optimal_threshold, + 'optimization_metric': optimization_metric, + 'use_gpu': use_gpu, + 'seed': seed + }) + estimator.fit(var_table_train, outcome_table_binary_train) # Saving the trained model using pickle - name_save_model = ml['algorithms']['XGBoost']['nameSave'] + name_save_model = ml['algorithms'][algorithm]['nameSave'] model_id = name_save_model + '_' + str(ml['variables']['varStudy']) path_model = os.path.dirname(path_results) + '/' + (model_id + '.pickle') - model_dict = save_model(model, str(ml['variables']['varStudy']), path_model, ml=ml) + estimator.save(path_model) logging.info("{}--> DONE. TOTAL TIME OF LEARNING PROCESS: {:.2f} min".format(" " * 4, (time.time()-tstart) / 60)) # --> C. Testing phase - # C.1. Testing the XGBoost model and computing model response + # C.1. Testing the model and computing model response tstart = time.time() - logging.info(f"\n\n--> TESTING XGBOOST MODEL FOR VARIABLE {var_id}") + logging.info(f"\n\n--> TESTING {algorithm.upper()} MODEL FOR VARIABLE {var_id}") - response_train, response_test = self.test_xgb_model( - model, - processed_testing_table, - [patients_train, patients_test] - ) + # Preparing the variable table + var_table_test = get_ml_test_table(estimator, processed_testing_table) + + # Getting the model response for training and test sets + response_train = estimator.predict_proba(var_table_test.loc[patients_train, :]) + response_test = estimator.predict_proba(var_table_test.loc[patients_test, :]) logging.info('{}--> DONE. TOTAL TIME OF LEARNING PROCESS: {:.2f}'.format(" " * 4, (time.time() - tstart)/60)) @@ -620,11 +375,11 @@ def ml_run(self, path_ml: Path, holdout_test: bool = True, method: str = 'auto') # D.1. Prepare holdout test data var_table_all_holdout = self.get_hold_out_set_table(ml, var_id, patients_holdout) - # D.2. Testing the XGBoost model and computing model response on the holdout set + # D.2. Testing the model and computing model response on the holdout set tstart = time.time() - logging.info(f"\n\n--> TESTING XGBOOST MODEL FOR VARIABLE {var_id} ON THE HOLDOUT SET") + logging.info(f"\n\n--> TESTING {algorithm.upper()} MODEL FOR VARIABLE {var_id} ON THE HOLDOUT SET") - response_holdout = self.test_xgb_model(model, var_table_all_holdout, [patients_holdout])[0] + response_holdout = estimator.predict_proba(var_table_all_holdout.loc[patients_holdout, :]) logging.info('{}--> DONE. TOTAL TIME OF LEARNING PROCESS: {:.2f}'.format(" " * 4, (time.time() - tstart)/60)) @@ -632,7 +387,7 @@ def ml_run(self, path_ml: Path, holdout_test: bool = True, method: str = 'auto') tstart = time.time() # Initialize the Results class - result = Results(model_dict, model_id) + result = Results(estimator.estimator_.model_info_, model_id) if holdout_test: run_results = result.to_json( response_train=response_train, @@ -640,35 +395,19 @@ def ml_run(self, path_ml: Path, holdout_test: bool = True, method: str = 'auto') response_holdout=response_holdout, patients_train=patients_train, patients_test=patients_test, - patients_holdout=patients_holdout + patients_holdout=patients_holdout, + outcome_table_binary_train=outcome_table_binary_train, + outcome_table_binary_test=outcome_table_binary_test, + outcome_table_binary_holdout=outcome_table_binary_holdout ) else: run_results = result.to_json( response_train=response_train, response_test=response_test, - response_holdout=None, patients_train=patients_train, patients_test=patients_test, - patients_holdout=None - ) - - # Calculating performance metrics for training phase and saving the ROC curve - run_results[model_id]['train']['metrics'] = result.get_model_performance( - response_train, - outcome_table_binary_train, - ) - - # Calculating performance metrics for testing phase and saving the ROC curve - run_results[model_id]['test']['metrics'] = result.get_model_performance( - response_test, - outcome_table_binary_test, - ) - - if holdout_test: - # Calculating performance metrics for holdout phase and saving the ROC curve - run_results[model_id]['holdout']['metrics'] = result.get_model_performance( - response_holdout, - outcome_table_binary_holdout, + outcome_table_binary_train=outcome_table_binary_train, + outcome_table_binary_test=outcome_table_binary_test, ) logging.info('\n\n--> COMPUTING PERFORMANCE METRICS ... Done in {:.2f} sec'.format(time.time()-tstart)) @@ -687,7 +426,7 @@ def run_experiment(self, holdout_test: bool = True, method: str = "pycaret") -> Args: holdout_test (bool, optional): Boolean specifying if the hold-out test should be performed. - method (str, optional): String specifying the method to use to train the XGBoost model. + method (str, optional): String specifying the method to use to train the model. - "pycaret": Use PyCaret to train the model (automatic). - "grid_search": Grid search with cross-validation to find the best parameters. - "random_search": Random search with cross-validation to find the best parameters. @@ -710,5 +449,5 @@ def run_experiment(self, holdout_test: bool = True, method: str = "pycaret") -> average_results(self.path_study / f'learn__{self.experiment_label}', save=True) # Analyze the features importance for all the runs - feature_imporance_analysis(self.path_study / f'learn__{self.experiment_label}') + feature_importance_analysis(self.path_study / f'learn__{self.experiment_label}') \ No newline at end of file diff --git a/MEDiml/learning/Results.py b/MEDiml/learning/Results.py index 839b728..a6014fd 100644 --- a/MEDiml/learning/Results.py +++ b/MEDiml/learning/Results.py @@ -17,7 +17,7 @@ from numpyencoder import NumpyEncoder from sklearn import metrics -from MEDiml.learning.ml_utils import feature_imporance_analysis, list_metrics +from MEDiml.learning.ml_utils import feature_importance_analysis, METRICS_LIST from MEDiml.learning.Stats import Stats from MEDiml.utils.json_utils import load_json, save_json from MEDiml.utils.texture_features_names import * @@ -36,7 +36,7 @@ class Results: model_id (str): ID of the model. results_dict (dict): Dictionary containing the results of the model's performance. """ - def __init__(self, model_dict: dict = {}, model_id: str = "") -> None: + def __init__(self, model_dict: dict = None, model_id: str = "") -> None: """ Constructor of the class Results """ @@ -46,7 +46,7 @@ def __init__(self, model_dict: dict = {}, model_id: str = "") -> None: def __calculate_performance( self, - response: list, + response: pd.Series, labels: pd.DataFrame, thresh: float ) -> dict: @@ -61,70 +61,48 @@ def __calculate_performance( Returns: Dict: Dictionary containing the performance metrics. """ - # Recording results - results_dict = dict() - - # Removing Nans - df = labels.copy() - outcome_name = labels.columns.values[0] - df['response'] = response - df.dropna(axis=0, how='any', inplace=True) - - # Confusion matrix elements: - results_dict['TP'] = ((df['response'] >= thresh) & (df[outcome_name] == 1)).sum() - results_dict['TN'] = ((df['response'] < thresh) & (df[outcome_name] == 0)).sum() - results_dict['FP'] = ((df['response'] >= thresh) & (df[outcome_name] == 0)).sum() - results_dict['FN'] = ((df['response'] < thresh) & (df[outcome_name] == 1)).sum() + outcome_name = labels.columns[0] + # Align data and drop NaNs once + df = pd.concat([labels, response.rename('response')], axis=1).dropna() - # Copying confusion matrix elements - TP = results_dict['TP'] - TN = results_dict['TN'] - FP = results_dict['FP'] - FN = results_dict['FN'] - - # AUC - results_dict['AUC'] = metrics.roc_auc_score(df[outcome_name], df['response']) - - # AUPRC - results_dict['AUPRC'] = metrics.average_precision_score(df[outcome_name], df['response']) - - # Sensitivity - try: - results_dict['Sensitivity'] = TP / (TP + FN) - except: - print('TP + FN = 0, Division by 0, replacing sensitivity by 0.0') - results_dict['Sensitivity'] = 0.0 - - # Specificity - try: - results_dict['Specificity'] = TN / (TN + FP) - except: - print('TN + FP= 0, Division by 0, replacing specificity by 0.0') - results_dict['Specificity'] = 0.0 - - # Balanced accuracy - results_dict['BAC'] = (results_dict['Sensitivity'] + results_dict['Specificity']) / 2 - - # Precision - results_dict['Precision'] = TP / (TP + FP) - - # NPV (Negative Predictive Value) - results_dict['NPV'] = TN / (TN + FN) - - # Accuracy - results_dict['Accuracy'] = (TP + TN) / (TP + TN + FP + FN) - - # F1 score - results_dict['F1_score'] = 2 * TP / (2 * TP + FP + FN) - - # mcc (mathews correlation coefficient) - results_dict['MCC'] = (TP * TN - FP * FN) / np.sqrt((TP + FP) * (TP + FN) * (TN + FP) * (TN + FN)) + y_true = df[outcome_name] + y_prob = df['response'] + y_pred = (y_prob >= thresh).astype(int) + + # Vectorized Confusion Matrix + tn, fp, fn, tp = metrics.confusion_matrix(y_true, y_pred, labels=[0, 1]).ravel() + + # Pre-calculate sums to avoid repeated addition + actual_pos = tp + fn + actual_neg = tn + fp + pred_pos = tp + fp + pred_neg = tn + fn + total = actual_pos + actual_neg + + # Calculate metrics using vectorized scalars + res = { + 'TP': int(tp), 'TN': int(tn), 'FP': int(fp), 'FN': int(fn), + 'AUC': metrics.roc_auc_score(y_true, y_prob), + 'AUPRC': metrics.average_precision_score(y_true, y_prob), + 'Sensitivity': tp / actual_pos if actual_pos > 0 else 0.0, + 'Specificity': tn / actual_neg if actual_neg > 0 else 0.0, + 'Precision': tp / pred_pos if pred_pos > 0 else 0.0, + 'NPV': tn / pred_neg if pred_neg > 0 else 0.0, + 'Accuracy': (tp + tn) / total if total > 0 else 0.0, + } + + res['BAC'] = (res['Sensitivity'] + res['Specificity']) / 2 + res['F1_score'] = 2 * tp / (2 * tp + fp + fn) if (2 * tp + fp + fn) > 0 else 0.0 + + # Matthews Correlation Coefficient + mcc_denom = np.sqrt(pred_pos * actual_pos * actual_neg * pred_neg) + res['MCC'] = (tp * tn - fp * fn) / mcc_denom if mcc_denom > 0 else 0.0 - return results_dict + return res def __get_metrics_failure_dict( self, - metrics: list = list_metrics + metrics: list = METRICS_LIST ) -> dict: """ This function fills the metrics with NaNs in case of failure. @@ -218,7 +196,7 @@ def __count_percentage_radiomics(self, results_dict: dict) -> list: for key in list(radiomics_tables_dict.keys()): if key.lower().startswith('radtab'): table_path = radiomics_tables_dict[key]['original_data']['path_radiomics_csv'] - table_name = table_path.split('/')[-1] + table_name = Path(table_path).name table = pd.read_csv(table_path, index_col=0) # Morph if 'morph' in table_name.lower(): @@ -320,8 +298,8 @@ def __count_patients(self, path_results: Path) -> dict: break # The number of patients is the same for all the runs return patients_count - - def average_results(self, path_results: Path, save: bool = False) -> None: + + def average_results(self, path_results: Path, save: bool = False) -> dict: """ Averages the results (AUC, BAC, Sensitivity and Specifity) of all the runs of the same experiment, for training, testing and holdout sets. @@ -331,49 +309,97 @@ def average_results(self, path_results: Path, save: bool = False) -> None: save (bool, optional): If True, saves the results in the same folder as the model. Returns: - None. + dict: Averaged results for each dataset. """ - # Get all tests paths - list_path_tests = [path for path in path_results.iterdir() if path.is_dir()] - - # Initialize dictionaries - results_avg = { - 'train': {}, - 'test': {}, - 'holdout': {} - } + list_path_tests = [path / 'run_results.json' for path in path_results.iterdir() if path.is_dir()] - # Retrieve metrics - for dataset in ['train', 'test', 'holdout']: - dataset_dict = results_avg[dataset] - for metric in list_metrics: - metric_values = [] - for path_test in list_path_tests: - results_dict = load_json(path_test / 'run_results.json') - if dataset in results_dict[list(results_dict.keys())[0]].keys(): - if 'metrics' in results_dict[list(results_dict.keys())[0]][dataset].keys(): - metric_values.append(results_dict[list(results_dict.keys())[0]][dataset]['metrics'][metric]) - else: - continue - else: - continue - - # Fill the dictionary - if metric_values: - dataset_dict[f'{metric}_mean'] = np.nanmean(metric_values) - dataset_dict[f'{metric}_std'] = np.nanstd(metric_values) - dataset_dict[f'{metric}_max'] = np.nanmax(metric_values) - dataset_dict[f'{metric}_min'] = np.nanmin(metric_values) - dataset_dict[f'{metric}_2.5%'] = np.nanpercentile(metric_values, 2.5) - dataset_dict[f'{metric}_97.5%'] = np.nanpercentile(metric_values, 97.5) + all_metrics = [] + for p in list_path_tests: + data = load_json(p) + model_key = list(data.keys())[0] + + # Flatten metrics into a list of dicts with 'dataset' as a key + for ds in ['train', 'test', 'holdout']: + if ds in data[model_key] and 'metrics' in data[model_key][ds]: + m = data[model_key][ds]['metrics'].copy() + m['dataset'] = ds + all_metrics.append(m) + + if not all_metrics: + return {} + + df_all = pd.DataFrame(all_metrics) + results_avg = {} + + for ds in ['train', 'test', 'holdout']: + ds_df = df_all[df_all['dataset'] == ds].drop(columns='dataset') + if ds_df.empty: + results_avg[ds] = {} + continue + + # Vectorized aggregation for all metrics at once + stats = ds_df.agg(['mean', 'std', 'max', 'min']).to_dict() + + # Flatten the nested stats into your required format + results_avg[ds] = { + f"{met}_{stat}": val + for met, s_dict in stats.items() + for stat, val in s_dict.items() + } - # Save the results if save: save_json(path_results / 'results_avg.json', results_avg, cls=NumpyEncoder) - return path_results / 'results_avg.json' - + return results_avg + def bootstrap_metrics( + self, + response: np.ndarray, + labels: pd.DataFrame, + thresh: float, + n_bootstraps: int = 100 + ) -> dict: + """ + Computes 95% Confidence Intervals using bootstrap resampling. + + Args: + response (np.ndarray): Array of the probabilities of class "1" for all instances (prediction). + labels (pd.DataFrame): Column vector specifying the outcome status (1 or 0) for all instances. + thresh (float): Optimal threshold selected from the ROC curve. + n_bootstraps (int, optional): Number of bootstrap samples. Defaults to 100. + + Returns: + dict: Dictionary containing the 95% confidence intervals for each metric. + """ + bootstrapped_stats = [] + rng = np.random.default_rng() + + # Ensure input is numpy for fast indexing + y_true = labels.iloc[:, 0].values + y_prob = np.array(response) + + for _ in range(n_bootstraps): + indices = rng.integers(0, len(y_true), len(y_true)) + if len(np.unique(y_true[indices])) < 2: + continue + + # Reuse the optimized calculation logic + res = self.__calculate_performance( + pd.Series(y_prob[indices]), + pd.DataFrame(y_true[indices]), + thresh + ) + bootstrapped_stats.append(res) + + df_boot = pd.DataFrame(bootstrapped_stats) + ci_results = {} + + for metric in df_boot.columns: + ci_results[f"{metric}_95ci_low"] = np.percentile(df_boot[metric], 2.5) + ci_results[f"{metric}_95ci_high"] = np.percentile(df_boot[metric], 97.5) + + return ci_results + def get_model_performance( self, response: list, @@ -443,7 +469,7 @@ def get_optimal_level( Returns: None. """ - assert metric.split('_')[0] in list_metrics, f'Given metric {list_metrics} is not in the list of metrics. Please choose from {list_metrics}' + assert metric.split('_')[0] in METRICS_LIST, f'Given metric {METRICS_LIST} is not in the list of metrics. Please choose from {METRICS_LIST}' # Extract modalities and initialize the dictionary if type(experiments_labels[0]) == str: @@ -680,7 +706,7 @@ def plot_heatmap( Returns: None. """ - assert metric.split('_')[0] in list_metrics, f'Given metric {list_metrics} is not in the list of metrics. Please choose from {list_metrics}' + assert metric.split('_')[0] in METRICS_LIST, f'Given metric {METRICS_LIST} is not in the list of metrics. Please choose from {METRICS_LIST}' # Extract modalities and initialize the dictionary if type(experiments_labels[0]) == str: @@ -1200,7 +1226,7 @@ def plot_feature_analysis( if 'feature_importance_analysis.json' in os.listdir(path_experiments / exp_full_name): fa_dict = load_json(path_experiments / exp_full_name / 'feature_importance_analysis.json') else: - fa_dict = feature_imporance_analysis(path_experiments / exp_full_name) + fa_dict = feature_importance_analysis(path_experiments / exp_full_name) # Extract percentage of features per level perc_levels = np.round(self.__count_percentage_levels(fa_dict), 2) @@ -1304,7 +1330,7 @@ def plot_original_level_tree( if 'feature_importance_analysis.json' in os.listdir(path_experiments / exp_full_name): fa_dict = load_json(path_experiments / exp_full_name / 'feature_importance_analysis.json') else: - fa_dict = feature_imporance_analysis(path_experiments / exp_full_name) + fa_dict = feature_importance_analysis(path_experiments / exp_full_name) # Organize data feature_data = { @@ -1617,7 +1643,7 @@ def plot_lf_level_tree( if 'feature_importance_analysis.json' in os.listdir(path_experiments / exp_full_name): fa_dict = load_json(path_experiments / exp_full_name / 'feature_importance_analysis.json') else: - fa_dict = feature_imporance_analysis(path_experiments / exp_full_name) + fa_dict = feature_importance_analysis(path_experiments / exp_full_name) # Organize data feature_data = { @@ -1944,7 +1970,7 @@ def plot_tf_level_tree( if 'feature_importance_analysis.json' in os.listdir(path_experiments / exp_full_name): fa_dict = load_json(path_experiments / exp_full_name / 'feature_importance_analysis.json') else: - fa_dict = feature_imporance_analysis(path_experiments / exp_full_name) + fa_dict = feature_importance_analysis(path_experiments / exp_full_name) # Organize data feature_data = { @@ -2190,7 +2216,7 @@ def plot_tf_level_tree( # Save the plot (Mandatory, since the plot is not well displayed on matplotlib) fig.savefig(path_experiments / f'TF_{experiment}_{level}_{modality}_explanation_tree.png', dpi=300) - + def to_json( self, response_train: list = None, @@ -2198,7 +2224,10 @@ def to_json( response_holdout: list = None, patients_train: list = None, patients_test: list = None, - patients_holdout: list = None + patients_holdout: list = None, + outcome_table_binary_train: pd.DataFrame = None, + outcome_table_binary_test: pd.DataFrame = None, + outcome_table_binary_holdout: pd.DataFrame = None ) -> dict: """ Creates a dictionary with the results of the model using the class attributes. @@ -2206,32 +2235,45 @@ def to_json( Args: response_train (list): List of machine learning model predictions for the training set. response_test (list): List of machine learning model predictions for the test set. + response_holdout (list): List of machine learning model predictions for the holdout set. patients_train (list): List of patients in the training set. patients_test (list): List of patients in the test set. patients_holdout (list): List of patients in the holdout set. + outcome_table_binary_train (pd.DataFrame): Binary outcome table for the training set. + outcome_table_binary_test (pd.DataFrame): Binary outcome table for the test set. + outcome_table_binary_holdout (pd.DataFrame): Binary outcome table for the holdout set. Returns: Dict: Dictionary with the the responses of the model and the patients used for training, testing and holdout. """ - run_results = dict() - run_results[self.model_id] = self.model_dict - - # Training results info - run_results[self.model_id]['train'] = dict() - run_results[self.model_id]['train']['patients'] = patients_train - run_results[self.model_id]['train']['response'] = response_train.tolist() if response_train is not None else [] + # Initialization + run_results = {self.model_id: self.model_dict} + threshold = self.model_dict.get('threshold', 0.5) + + # Map datasets for cleaner iteration + datasets = { + 'train': (response_train, patients_train, outcome_table_binary_train), + 'test': (response_test, patients_test, outcome_table_binary_test), + 'holdout': (response_holdout, patients_holdout, outcome_table_binary_holdout) + } - # Testing results info - run_results[self.model_id]['test'] = dict() - run_results[self.model_id]['test']['patients'] = patients_test - run_results[self.model_id]['test']['response'] = response_test.tolist() if response_test is not None else [] + for name, (resp, pts, outcome) in datasets.items(): + run_results[self.model_id][name] = { + 'patients': pts, + 'response': resp.tolist() if hasattr(resp, 'tolist') else (resp or []) + } - # Holdout results info - run_results[self.model_id]['holdout'] = dict() - run_results[self.model_id]['holdout']['patients'] = patients_holdout - run_results[self.model_id]['holdout']['response'] = response_holdout.tolist() if response_holdout is not None else [] + # Only calculate if we have both predictions and ground truth + if resp is not None and outcome is not None: + # 1. Vectorized Point Estimates + metrics_dict = self.get_model_performance(resp, outcome) + + # 2. Bootstrap Confidence Intervals (95% CI) + ci_dict = self.bootstrap_metrics(resp, outcome, threshold) + + # Merge both into the metrics entry + run_results[self.model_id][name]['metrics'] = {**metrics_dict, **ci_dict} - # keep a copy of the results self.results_dict = run_results return run_results diff --git a/MEDiml/learning/__init__.py b/MEDiml/learning/__init__.py index 9c355f8..8e274ec 100644 --- a/MEDiml/learning/__init__.py +++ b/MEDiml/learning/__init__.py @@ -4,7 +4,7 @@ from .DesignExperiment import DesignExperiment from .FSR import FSR from .ml_utils import * -from .Normalization import Normalization +from .Normalization import CombatNormalization from .RadiomicsLearner import RadiomicsLearner from .Results import Results from .Stats import Stats diff --git a/MEDiml/learning/ml_utils.py b/MEDiml/learning/ml_utils.py index cca73e0..a43d1c3 100644 --- a/MEDiml/learning/ml_utils.py +++ b/MEDiml/learning/ml_utils.py @@ -8,82 +8,94 @@ from pathlib import Path from typing import Dict, List, Tuple, Union -import matplotlib.pyplot as plt import numpy as np import pandas import pandas as pd -import seaborn as sns from numpyencoder import NumpyEncoder +from sklearn.base import BaseEstimator from sklearn.model_selection import StratifiedKFold from MEDiml.utils import get_institutions_from_ids from MEDiml.utils.get_full_rad_names import get_full_rad_names from MEDiml.utils.json_utils import load_json, save_json - # Define useful constants # Metrics to process -list_metrics = [ +METRICS_LIST = [ 'AUC', 'AUPRC', 'BAC', 'Sensitivity', 'Specificity', 'Precision', 'NPV', 'F1_score', 'Accuracy', 'MCC', 'TN', 'FP', 'FN', 'TP' ] -def average_results(path_results: Path, save: bool = False) -> None: +def average_results(path_results: Path, save: bool = False) -> dict: """ - Averages the results (AUC, BAC, Sensitivity and Specifity) of all the runs of the same experiment, - for training, testing and holdout sets. + Averages the results (including mean, std, and percentiles) of all the runs + of the same experiment for training, testing, and holdout sets. Args: path_results(Path): path to the folder containing the results of the experiment. save (bool, optional): If True, saves the results in the same folder as the model. Returns: - None. + Dict: Aggregated results for all datasets. """ - # Get all tests paths - list_path_tests = [path for path in path_results.iterdir() if path.is_dir()] + # 1. Optimized File I/O: Get all file paths first + list_path_tests = [path / 'run_results.json' for path in path_results.iterdir() if path.is_dir()] + + # 2. Extract data into a flat list of records (One pass through files) + records = [] + for file_path in list_path_tests: + if not file_path.exists(): + continue + + data = load_json(file_path) + model_id = list(data.keys())[0] + model_data = data[model_id] + + for dataset in ['train', 'test', 'holdout']: + metrics_data = model_data.get(dataset, {}).get('metrics') + if metrics_data: + # Flatten metrics and tag with dataset name + row = {**metrics_data, 'dataset_name': dataset} + records.append(row) + + if not records: + return {} + + # 3. Vectorized Aggregation using Pandas + df = pd.DataFrame(records) + results_avg = {} - # Initialize dictionaries - results_avg = { - 'train': {}, - 'test': {}, - 'holdout': {} - } + for dataset in ['train', 'test', 'holdout']: + ds_df = df[df['dataset_name'] == dataset].drop(columns='dataset_name') + + # Skipt if empty + if ds_df.empty: + results_avg[dataset] = {} + continue - # Metrics to process - metrics = ['AUC', 'AUPRC', 'BAC', 'Sensitivity', 'Specificity', - 'Precision', 'NPV', 'F1_score', 'Accuracy', 'MCC', - 'TN', 'FP', 'FN', 'TP'] + # Calculate all stats at once across all metric columns + summary = ds_df.apply(lambda x: pd.Series({ + 'mean': np.nanmean(x), + 'std': np.nanstd(x), + 'max': np.nanmax(x), + 'min': np.nanmin(x), + '2.5%': np.nanpercentile(x, 2.5), + '97.5%': np.nanpercentile(x, 97.5) + })) + + # Reshape into the requested {metric}_{stat} format + results_avg[dataset] = { + f"{col}_{stat}": val + for col in summary.columns + for stat, val in summary[col].items() + } - # Process metrics - for dataset in ['train', 'test', 'holdout']: - dataset_dict = results_avg[dataset] - for metric in metrics: - metric_values = [] - for path_test in list_path_tests: - results_dict = load_json(path_test / 'run_results.json') - if dataset in results_dict[list(results_dict.keys())[0]].keys(): - if 'metrics' in results_dict[list(results_dict.keys())[0]][dataset].keys(): - metric_values.append(results_dict[list(results_dict.keys())[0]][dataset]['metrics'][metric]) - else: - continue - else: - continue - - # Fill the dictionary - if metric_values: - dataset_dict[f'{metric}_mean'] = np.nanmean(metric_values) - dataset_dict[f'{metric}_std'] = np.nanstd(metric_values) - dataset_dict[f'{metric}_max'] = np.nanmax(metric_values) - dataset_dict[f'{metric}_min'] = np.nanmin(metric_values) - dataset_dict[f'{metric}_2.5%'] = np.nanpercentile(metric_values, 2.5) - dataset_dict[f'{metric}_97.5%'] = np.nanpercentile(metric_values, 97.5) - - # Save the results + # 4. Save and Return if save: - save_json(path_results / 'results_avg.json', results_avg, cls=NumpyEncoder) - return path_results / 'results_avg.json' + save_path = path_results / 'results_avg.json' + save_json(save_path, results_avg, cls=NumpyEncoder) + return save_path return results_avg @@ -483,13 +495,13 @@ def find_best_model(path_results: Path, metric: str = 'AUC', second_metric: str Returns: Tuple[Dict, Path]: Tuple containing the best model result dict and the path to the best model. """ - list_metrics = [ + METRICS_LIST = [ 'AUC', 'Sensitivity', 'Specificity', 'BAC', 'AUPRC', 'Precision', 'NPV', 'Accuracy', 'F1_score', 'MCC', 'TP', 'TN', 'FP', 'FN' ] - assert metric in list_metrics, f'Given metric {metric} is not in the list of metrics. Please choose from {list_metrics}' + assert metric in METRICS_LIST, f'Given metric {metric} is not in the list of metrics. Please choose from {METRICS_LIST}' # Get all tests paths list_path_tests = [path for path in path_results.iterdir() if path.is_dir()] @@ -524,58 +536,66 @@ def find_best_model(path_results: Path, metric: str = 'AUC', second_metric: str return model, results_dict_best -def feature_imporance_analysis(path_results: Path): +def feature_importance_analysis(path_results: Path): """ - Averages the results (AUC, BAC, Sensitivity and Specifity) of all the runs of the same experiment, - for training, testing and holdout sets. + Analyzes and averages feature importance across all experimental runs. + Calculates the mean importance (among selections) and selection frequency. Args: - path_results(Path): path to the folder containing the results of the experiment. - save (bool, optional): If True, saves the results in the same folder as the model. - - Returns: - None. + path_results (Path): Path to the folder containing the experiment run directories. """ - # Get all tests paths - list_path_tests = [path for path in path_results.iterdir() if path.is_dir()] + list_path_tests = [path for path in path_results.iterdir() if path.is_dir()] + importance_accumulator = {} - # Initialization - results_avg_temp = {} - results_avg = {} - - # Process metrics for path_test in list_path_tests: - variables = [] list_models = list(path_test.glob('*.pickle')) - if len(list_models) == 0 or len(list_models) > 1: - raise ValueError(f'Path {path_test} does not contain a single model.') - model_obj = list_models[0] - with open(model_obj, "rb") as f: - model_dict = pickle.load(f) - if model_dict["var_names"]: - variables = get_full_rad_names(model_dict['var_info']['variables']['var_def'], model_dict["var_names"]) - for index, var in enumerate(variables): - var = var.split("\\")[-1] # Remove the path for windows - var = var.split("/")[-1] # Remove the path for linux - if var not in results_avg_temp: - results_avg_temp[var] = { - 'importance_mean': [], - 'times_selected': 0 - } - - results_avg_temp[var]['importance_mean'].append(model_dict['model'].feature_importances_[index]) - results_avg_temp[var]['times_selected'] += 1 - for var in results_avg_temp: - results_avg[var] = { - 'importance_mean': np.sum(results_avg_temp[var]['importance_mean']) / len(list_path_tests), - 'times_selected': results_avg_temp[var]['times_selected'] + + if len(list_models) != 1: + print(f"Skipping {path_test}: Expected 1 pickle model, found {len(list_models)}.") + continue + + import joblib + pipeline = joblib.load(list_models[0]) + + # Extract feature names and importances + if hasattr(pipeline, 'estimator_') \ + and hasattr(pipeline.estimator_, 'model_info_') \ + and hasattr(pipeline.estimator_, 'classifier_') \ + and hasattr(pipeline.estimator_.classifier_, 'feature_importances_'): + variables = get_full_rad_names( + pipeline.estimator_.model_info_['var_info']['variables']['var_def'], + pipeline.estimator_.model_info_['var_names'] + ) + importances = pipeline.estimator_.classifier_.feature_importances_ + + # Accumulate importance values for each variable + for index, var_path in enumerate(variables): + var_name = Path(var_path).name + + if var_name not in importance_accumulator: + importance_accumulator[var_name] = [] + + importance_accumulator[var_name].append(importances[index]) + + # Aggregate results + total_runs = len(list_path_tests) + final_analysis = { + var: { + 'importance_mean': np.sum(vals) / total_runs if total_runs > 0 else 0, # Average over all runs (including zeros for non-selections) + 'importance_std': np.std(vals), + 'times_selected': len(vals), + 'selection_frequency': (len(vals) / total_runs) * 100 if total_runs > 0 else 0 } - - del results_avg_temp - - save_json(path_results / 'feature_importance_analysis.json', results_avg, cls=NumpyEncoder) + for var, vals in importance_accumulator.items() + } + + # Sort by importance_mean descending for better readability + final_analysis = dict(sorted(final_analysis.items(), key=lambda item: item[1]['importance_mean'], reverse=True)) -def get_ml_test_table(variable_table: pd.DataFrame, var_names: List, var_def: str) -> pd.DataFrame: + save_json(path_results / 'feature_importance_analysis.json', final_analysis, cls=NumpyEncoder) + return final_analysis + +def get_ml_test_table(estimator: BaseEstimator, variable_table: pd.DataFrame) -> pd.DataFrame: """ Gets the test table with the variables that are present in the training table. @@ -589,6 +609,10 @@ def get_ml_test_table(variable_table: pd.DataFrame, var_names: List, var_def: st pd.DataFrame: Table with the variables that are present in the training table. """ + # retrieve the necessary information from the variable table + var_names = estimator.estimator_.model_info_['var_names'] + var_def = estimator.estimator_.model_info_['var_info']['variables']['var_def'] + # Get the full variable names for training full_radvar_names_trained = get_full_rad_names(var_def, var_names).tolist() diff --git a/MEDiml/utils/rf_learner.py b/MEDiml/utils/rf_learner.py new file mode 100644 index 0000000..23faf7b --- /dev/null +++ b/MEDiml/utils/rf_learner.py @@ -0,0 +1,137 @@ +from copy import deepcopy + +import numpy as np +import pandas as pd +from pycaret.classification import * +from sklearn import metrics +from sklearn.base import BaseEstimator, ClassifierMixin + +from ..learning.ml_utils import finalize_rad_table, intersect_var_tables + + +class RandomForestEstimator(BaseEstimator, ClassifierMixin): + def __init__( + self, + optimization_metric='MCC', + var_importance_threshold=0.05, + internal_cv_folds=5, + optimal_threshold=None, + use_gpu=False, + seed=None + ): + self.optimization_metric = optimization_metric + self.var_importance_threshold = var_importance_threshold + self.internal_cv_folds = internal_cv_folds + self.optimal_threshold = optimal_threshold + self.use_gpu = use_gpu + self.seed = seed + + self.model_info_ = None + self.classifier_ = None + self.selected_features_ = None + self.selected_features_definitions_ = None + + def fit(self, X, y): + if not isinstance(X, pd.DataFrame): + X = pd.DataFrame(X) + + # Ensure y is a DataFrame for merging in PyCaret logic + if not isinstance(y, pd.DataFrame): + y = pd.DataFrame(y) + + results, self.classifier_ = self._train_logic(X, y) + + self.model_info_ = results + self.selected_features_ = results['var_names'] + self.selected_features_definitions_ = results.get('var_def') + self.classes_ = np.unique(y) + + return self + + def predict(self, X): + if self.classifier_ is None: + raise ValueError("Model has not been fitted yet.") + + probas = self.predict_proba(X) + threshold = self.model_info_.get('threshold', 0.5) + return (probas >= threshold).astype(int) + + def predict_proba(self, X): + if not isinstance(X, pd.DataFrame): + X = pd.DataFrame(X) + + # Filter X to include only features selected during fit + X_filtered = X[self.selected_features_] + return self.classifier_.predict_proba(X_filtered)[:, 1] + + def _train_logic(self, var_table_train, outcome_table_binary_train): + # Align tables + var_table_train, outcome_table_binary_train = intersect_var_tables(var_table_train, outcome_table_binary_train) + var_table_train = finalize_rad_table(var_table_train) + + # Merge for PyCaret + temp_data = pd.merge(var_table_train, outcome_table_binary_train, left_index=True, right_index=True) + target_col = outcome_table_binary_train.columns[0] + + # PyCaret setup + setup( + data=temp_data, + target=target_col, + feature_selection=True, + n_features_to_select=1-self.var_importance_threshold, + fold=self.internal_cv_folds, + use_gpu=self.use_gpu, + feature_selection_estimator="rf", + session_id=self.seed, + html=False, + verbose=False + ) + + if self.seed is not None: + set_config('seed', self.seed) + + # Create RF model. 'balanced' is crucial for your 5/143 imbalance. + # This penalizes mistakes on the 5 progressors more heavily. + classifier = create_model('rf', class_weight='balanced', verbose=False) + + # Tune model + classifier = tune_model(classifier, optimize=self.optimization_metric, verbose=False) + + # Assemble dictionary + model_rf = dict() + model_rf['algo'] = 'rf' + model_rf['type'] = 'binary' + + # Find threshold + try: + model_rf['threshold'] = self.__find_balanced_threshold(classifier, var_table_train, outcome_table_binary_train) + except Exception as e: + print(f'Threshold calculation failed: {e}. Defaulting to 0.5') + model_rf['threshold'] = 0.5 + + # Store metadata safely + user_data = var_table_train.Properties.get('userData', {}) if hasattr(var_table_train, 'Properties') else {} + model_rf['var_info'] = deepcopy(user_data) + model_rf['var_def'] = deepcopy(user_data.get('variables', {}).get('var_def')) + model_rf['var_names'] = list(classifier.feature_names_in_) + model_rf['optimization'] = classifier.get_params() + + return model_rf, classifier + + def __find_balanced_threshold(self, model, variable_table, outcome_table_binary) -> float: + # Align features + if hasattr(model, 'feature_names_in_'): + variable_table = variable_table[list(model.feature_names_in_)] + + # Get probabilities + y_probs = model.predict_proba(variable_table)[:, 1] + + # ROC Calculation + fpr, tpr, thresholds = metrics.roc_curve(outcome_table_binary.iloc[:, 0], y_probs) + + # Geometric optimization (closest to top-left corner) + # Distance = sqrt( fpr^2 + (1-tpr)^2 ) + dist = np.sqrt(np.power(fpr, 2) + np.power(1 - tpr, 2)) + best_idx = np.argmin(dist) + + return thresholds[best_idx] diff --git a/MEDiml/utils/xgboost_learner.py b/MEDiml/utils/xgboost_learner.py new file mode 100644 index 0000000..eba3e1f --- /dev/null +++ b/MEDiml/utils/xgboost_learner.py @@ -0,0 +1,161 @@ +from copy import deepcopy + +import numpy as np +import pandas as pd +from pycaret.classification import * +from sklearn import metrics +from sklearn.base import BaseEstimator, ClassifierMixin + +from ..learning.ml_utils import finalize_rad_table, intersect_var_tables + + +class XGBoostEstimator(BaseEstimator, ClassifierMixin): + def __init__( + self, + optimization_metric='MCC', + var_importance_threshold=0.05, + internal_cv_folds=5, + optimal_threshold=None, + use_gpu=False, + seed=None + ): + # Store all parameters as attributes + self.optimization_metric = optimization_metric + self.var_importance_threshold = var_importance_threshold + self.internal_cv_folds = internal_cv_folds + self.optimal_threshold = optimal_threshold + self.use_gpu = use_gpu + self.seed = seed + + # This will hold the "model_xgb" dictionary result + self.model_info_ = None + self.classifier_ = None + self.selected_features_ = None + + def fit(self, X, y): + # 1. Standardize input format (ensure DataFrame) + if not isinstance(X, pd.DataFrame): + X = pd.DataFrame(X) + + # 2. Call your existing logic + # Note: I am assuming 'intersect_var_tables' and 'finalize_rad_table' + # are available in your namespace. + results, self.classifier_ = self._train_logic(X, y) + + # 3. Store results for sklearn + self.model_info_ = results + self.selected_features_ = results['var_names'] + self.selected_features_definitions_ = results['var_def'] + self.classes_ = np.unique(y) + + return self + + def predict(self, X): + # Safety check + if self.selected_features_ is None or self.selected_features_definitions_ is None or self.features_names_in_ is None: + raise ValueError("Model has no selected features or definitions. " \ + "Ensure that fit() has been called successfully before predict().") + # Apply the threshold stored in model_info_ + probas = self.predict_proba(X) + threshold = self.model_info_.get('threshold', 0.5) + return (probas >= threshold).astype(int) + + def predict_proba(self, X): + if not isinstance(X, pd.DataFrame): + X = pd.DataFrame(X) + + # Important: Filter X to only include features selected during fit + features_names = self.selected_features_ or self.selected_features_definitions_ or self.features_names_in_ + X_filtered = X[features_names] + return self.classifier_.predict_proba(X_filtered)[:, 1] + + def _train_logic(self, var_table_train, outcome_table_binary_train): + """ + Trains an XGBoost model for the given machine learning test. + + Args: + var_table_train (pd.DataFrame): Radiomics table for the training/learning set. + outcome_table_binary_train (pd.DataFrame): Outcome table with binary labels for the training/learning set. + + Returns: + Dict: Dictionary containing info about the trained XGBoost model. + """ + # Safety check (make sure that the outcome table and the variable table have the same patients) + var_table_train, outcome_table_binary_train = intersect_var_tables(var_table_train, outcome_table_binary_train) + + # Finalize the new radiomics table with the remaining variables + var_table_train = finalize_rad_table(var_table_train) + + # Set up data for PyCaret + temp_data = pd.merge(var_table_train, outcome_table_binary_train, left_index=True, right_index=True) + + # PyCaret setup + setup( + data=temp_data, + feature_selection=True, + n_features_to_select=1-self.var_importance_threshold, + fold=self.internal_cv_folds, + target=temp_data.columns[-1], + use_gpu=self.use_gpu, + feature_selection_estimator="xgboost", + session_id=self.seed + ) + + # Set seed + if self.seed is not None: + set_config('seed', self.seed) + + # Creating XGBoost model using PyCaret + classifier = create_model('xgboost', verbose=False) + + # Tuning XGBoost model using PyCaret + classifier = tune_model(classifier, optimize=self.optimization_metric) + + # Saving the information of the model in a dictionary + model_xgb = dict() + model_xgb['algo'] = 'xgb' + model_xgb['type'] = 'binary' + try: + model_xgb['threshold'] = self.__find_balanced_threshold(classifier, var_table_train, outcome_table_binary_train) + except Exception as e: + print('Error in finding optimal threshold, it will be set to 0.5:' + str(e)) + model_xgb['threshold'] = 0.5 + model_xgb['var_info'] = deepcopy(var_table_train.Properties['userData']) + model_xgb['var_def'] = deepcopy(var_table_train.Properties['userData']['variables']['var_def']) + model_xgb['var_names'] = list(classifier.feature_names_in_) + model_xgb['optimization'] = classifier.get_params() + + return model_xgb, classifier + + def __find_balanced_threshold( + self, + model: object, + variable_table: pd.DataFrame, + outcome_table_binary: pd.DataFrame + ) -> float: + """ + Finds the balanced threshold for the given machine learning test. + + Args: + model (XGBClassifier): Trained XGBoost classifier for the given machine learning run. + variable_table (pd.DataFrame): Radiomics table. + outcome_table_binary (pd.DataFrame): Outcome table with binary labels. + + Returns: + float: Balanced threshold for the given machine learning test. + """ + # Check is there is a feature mismatch + if model.feature_names_in_.shape[0] != variable_table.columns.shape[0]: + variable_table = variable_table.loc[:, model.feature_names_in_] + + # Getting the probability responses for each patient + patient_ids = list(variable_table.index.values) + prob_xgb = self.predict(variable_table.loc[patient_ids, :]) + + # Calculating the ROC curve + fpr, tpr, thresholds = metrics.roc_curve(outcome_table_binary.iloc[:, 0], prob_xgb) + + # Calculating the optimal threshold by minizing fpr (false positive rate) and maximizing tpr (true positive rate) + minimum = np.argmin(np.power(fpr, 2) + np.power(1-tpr, 2)) + + return thresholds[minimum] \ No newline at end of file From 23961b1a79f7a078bfa6a9c7403fce92e306652f Mon Sep 17 00:00:00 2001 From: MahdiAll99 Date: Tue, 3 Mar 2026 18:16:57 -0500 Subject: [PATCH 11/15] Merged all json settings files into one single yml file --- MEDiml/learning/DesignExperiment.py | 211 +++++++++------------------- MEDiml/learning/RadiomicsLearner.py | 40 +++--- MEDiml/utils/rf_learner.py | 15 +- MEDiml/utils/xgboost_learner.py | 18 ++- 4 files changed, 106 insertions(+), 178 deletions(-) diff --git a/MEDiml/learning/DesignExperiment.py b/MEDiml/learning/DesignExperiment.py index 03f1b8f..b1d7ac5 100644 --- a/MEDiml/learning/DesignExperiment.py +++ b/MEDiml/learning/DesignExperiment.py @@ -5,6 +5,7 @@ from typing import Dict, List import pandas as pd +import yaml from ..utils.get_institutions_from_ids import get_institutions_from_ids from ..utils.json_utils import load_json, posix_to_string, save_json @@ -19,7 +20,7 @@ def __init__(self, path_study: Path, path_settings: Path, experiment_label: str) Args: path_study (Path): Path to the main study folder where the outcomes, learning patients and holdout patients dictionaries are found. - path_settings (Path): Path to the settings folder. + path_settings (Path): Path to the settings file. experiment_label (str): String specifying the label to attach to a given learning experiment in "path_experiments". This label will be attached to the ml__$experiments_label$.json file as well as the learn__$experiment_label$ folder. This label is used to keep track of different experiments @@ -72,84 +73,56 @@ def __create_folder_and_content( save_json(path_file, paths_ml) return ml_path + + def __load_config(self) -> Dict: + """Loads the YAML master configuration file.""" + with open(self.path_settings, 'r') as file: + return yaml.safe_load(file) - def generate_learner_dict(self) -> dict: + def __get_learning_dict(self) -> Path: """ - Generates a dictionary containing all the settings for the learning experiment. + Generates a dictionary containing all settings for the learning experiment + using a single YAML master configuration. Returns: - dict: Dictionary containing all the settings for the learning experiment. + Path: Path to the saved experiment-specific YAML options. """ - ml_options = dict() - - # operating system - ml_options['os'] = platform.system() - - # design experiment settings - ml_options['design'] = self.path_settings / 'ml_design.json' - # check if file exist: - if not ml_options['design'].exists(): - raise FileNotFoundError(f"File {ml_options['design']} does not exist.") - - # ML run settings - run = dict() - ml_options['run'] = run - - # Machine learning settings - ml_options['settings'] = self.path_settings / 'ml_settings.json' - # check if file exist: - if not ml_options['settings'].exists(): - raise FileNotFoundError(f"File {ml_options['settings']} does not exist.") - - # variables settings - ml_options['variables'] = self.path_settings / 'ml_variables.json' - # check if file exist: - if not ml_options['variables'].exists(): - raise FileNotFoundError(f"File {ml_options['variables']} does not exist.") - - # ML algorithms settings - ml_options['algorithms'] = self.path_settings / 'ml_algorithms.json' - # check if file exist: - if not ml_options['algorithms'].exists(): - raise FileNotFoundError(f"File {ml_options['algorithms']} does not exist.") - - # Data cleaning settings - ml_options['datacleaning'] = self.path_settings / 'ml_datacleaning.json' - # check if file exist: - if not ml_options['datacleaning'].exists(): - raise FileNotFoundError(f"File {ml_options['datacleaning']} does not exist.") - - # Normalization settings - ml_options['normalization'] = self.path_settings / 'ml_normalization.json' - # check if file exist: - if not ml_options['normalization'].exists(): - raise FileNotFoundError(f"File {ml_options['normalization']} does not exist.") - - # Feature set reduction settings - ml_options['fSetReduction'] = self.path_settings / 'ml_fset_reduction.json' - # check if file exist: - if not ml_options['fSetReduction'].exists(): - raise FileNotFoundError(f"File {ml_options['fSetReduction']} does not exist.") - - # Experiment label check - if self.experiment_label == "": + # Safety check: Verify master config exists + if not self.path_settings.exists(): + raise FileNotFoundError( + f"Master configuration file not found at: {self.path_settings}. " + "Please ensure the consolidated YAML is in the settings folder." + ) + + # Load the Master Config + config = self.__load_config() + + # Assemble Experiment Metadata + # We maintain the structure while adding run-specific info + ml_options = { + 'os': platform.system(), + 'experiment_label': self.experiment_label, + 'config_source': str(self.path_settings), + # Directly map the sections from the YAML for downstream use + 'design': config.get('design'), + 'variables': config.get('variables'), + 'datacleaning': config.get('data_cleaning'), + 'fSetReduction': config.get('feature_reduction'), + 'normalization': config.get('normalization'), + 'modeling': config.get('modeling'), + 'study_metadata': config.get('study_metadata') + } + + # Experiment Label Safety Check + if not self.experiment_label: raise ValueError("Experiment label is empty. Class was not initialized properly.") - - # save all the ml options and return the path to the saved file - name_save_options = 'ml_options_' + self.experiment_label + '.json' - path_ml_options = self.path_settings / name_save_options - ml_options = posix_to_string(ml_options) - save_json(path_ml_options, ml_options) - - return path_ml_options - def fill_learner_dict(self, path_ml_options: Path) -> Path: + return ml_options + + def __fill_learner_dict(self) -> Path: """ Fills the main expirement dictionary from the settings in the different json files. This main dictionary will hold all the settings for the data processing and learning experiment. - - Args: - path_ml_options (Path): Path to the ml_options json file for the experiment. Returns: Path: Path to the learner object. @@ -158,35 +131,19 @@ def fill_learner_dict(self, path_ml_options: Path) -> Path: all_datacleaning = list() all_normalization = list() all_fset_reduction = list() - - # Load ml options dict - ml_options = load_json(path_ml_options) - options = ml_options.keys() - - # Design options - ml = dict() - ml['design'] = load_json(ml_options['design']) - - # ML run options - ml['run'] = ml_options['run'] - - # Machine learning options - if 'settings' in options: - ml['settings'] = load_json(ml_options['settings']) + ml = self.__get_learning_dict() # Machine learning variables - if 'variables' in options: - ml['variables'] = dict() - var_options = load_json(ml_options['variables']) + if 'variables' in list(ml.keys()): + var_options = ml['variables'] fields = list(var_options.keys()) vars = [(idx, s) for idx, s in enumerate(fields) if re.match(r"^var[0-9]{1,}$", s)] var_names = [var[1] for var in vars] # list of var names - + # For each variable, organize the option in the ML dictionary for (idx, var) in vars: - vars_dict = dict() - vars_dict[var] = var_options[var] - var_struct = var_options[var] + vars_dict = ml['variables'] + var_struct = vars_dict[var] # Radiomics variables if 'radiomics' in var_struct['nameType'].lower(): @@ -248,12 +205,12 @@ def fill_learner_dict(self, path_ml_options: Path) -> Path: ml['variables'].update(vars_dict) # Initialize data processing methods - if 'var_datacleaning' in var_struct.keys(): - all_datacleaning.append(var_struct['var_datacleaning']) - if 'var_normalization' in var_struct.keys(): - all_normalization.append((var_struct['var_normalization'])) - if 'var_fSetReduction' in var_struct.keys(): - all_fset_reduction.append(var_struct['var_fSetReduction']['method']) + if 'cleaning_profile' in var_struct.keys(): + all_datacleaning.append(var_struct['cleaning_profile']) + if 'normalization' in var_struct.keys(): + all_normalization.append((var_struct['normalization'])) + if 'reduction_method' in var_struct.keys(): + all_fset_reduction.append(var_struct['reduction_method']) # Combinations of variables if 'combinations' in var_options.keys(): @@ -261,53 +218,18 @@ def fill_learner_dict(self, path_ml_options: Path) -> Path: combs = [comb for i in range(len(vars)) for comb in combinations(var_names, i+1)] combstrings = ['_'.join(elt) for elt in combs] ml['variables']['combinations'] = combstrings - else: - ml['variables']['combinations'] = var_options['combinations'] - - # Varibles to use for ML - ml['variables']['varStudy'] = var_options['varStudy'] - - # ML algorithms - if 'algorithms' in options: - algorithm = ml['settings']['algorithm'] - algorithms = load_json(ml_options['algorithms']) - ml['algorithms'] = {} - ml['algorithms'][algorithm] = algorithms[algorithm] - - # ML data processing methods and its options - for (method, method_list) in [ - ('datacleaning', all_datacleaning), - ('normalization', all_normalization), - ('fSetReduction', all_fset_reduction) - ]: - # Skip if no method is selected - if all(v == "" for v in method_list): - continue - if method in options: - # Add algorithm specific methods - if method in ml['settings'].keys(): - method_list.append(ml['settings'][method]) - method_list = list(set(method_list)) # to only get unique values of all_datacleaning - method_options = load_json(ml_options[method]) # load json file of each method - if method == 'normalization' and 'combat' in method_list: - ml[method] = 'combat' - continue - ml[method] = dict() - for name in list(set(method_list)): - if name != "": - ml[method][name] = method_options[name] # Save the ML dictionary if self.experiment_label == "": raise ValueError("Experiment label is empty. Class was not initialized properly.") - path_ml_object = self.path_study / f'ml__{self.experiment_label}.json' + path_ml_object = self.path_study / f'ml_test__{self.experiment_label}.json' ml = posix_to_string(ml) # Convert all paths to string save_json(path_ml_object, ml) # return ml return path_ml_object - def create_experiment(self, ml: dict = None) -> Dict: + def create_experiment(self) -> Dict: """ Create the machine learning experiment dictionary, organizes each test/split information in a seperate folder. @@ -319,7 +241,7 @@ def create_experiment(self, ml: dict = None) -> Dict: """ # Initialization ml_path = list() - ml = load_json(self.path_ml_object) if ml is None else ml + ml = load_json(self.path_ml_object) # Learning set patients_learn = load_json(self.path_study / 'patientsLearn.json') @@ -343,7 +265,7 @@ def create_experiment(self, ml: dict = None) -> Dict: Path.mkdir(path_learn, exist_ok=True) # Getting the type of test_sets - test_sets_types = ml['design']['testSets'] + test_sets_types = ml['design']['active_method'] # Creating the sets for the different machine learning runs for type_set in test_sets_types: @@ -427,13 +349,13 @@ def create_experiment(self, ml: dict = None) -> Dict: elif type_set.lower() == 'cv': # Get the experiment options for the sets cv_info = ml['design'][type_set] - n_splits = cv_info['nSplits'] + n_folds = cv_info['nFolds'] seed = cv_info['seed'] - + # Get the training and testing sets patients_train, patients_test = cross_validation_split( - outcomes_table, - n_splits, + outcomes_table, + n_folds, seed=seed ) @@ -441,8 +363,8 @@ def create_experiment(self, ml: dict = None) -> Dict: if type(patients_train) != list and not hasattr((patients_train), "__len__"): patients_train = [patients_train] patients_test = [patients_test] - - for i in range(n_splits): + + for i in range(n_folds): # Create a folder for each split/run run_name = "test__{0:03}".format(i+1) ml_path = self.__create_folder_and_content( @@ -463,11 +385,8 @@ def generate_experiment(self): Generate the json files containing all the options the experiment. The json files will then be used in machine learning. """ - # Generate the ml options dictionary - path_ml_options = self.generate_learner_dict() - # Fill the ml options dictionary - self.path_ml_object = self.fill_learner_dict(path_ml_options) + self.path_ml_object = self.__fill_learner_dict() # Generate the experiment dictionary experiment_dict = self.create_experiment() diff --git a/MEDiml/learning/RadiomicsLearner.py b/MEDiml/learning/RadiomicsLearner.py index bcfb58b..5e386b4 100644 --- a/MEDiml/learning/RadiomicsLearner.py +++ b/MEDiml/learning/RadiomicsLearner.py @@ -163,10 +163,10 @@ def pre_process_radiomics_table( # Initialization patient_ids = list(outcome_table_binary.index) outcome_table_binary_training = outcome_table_binary.loc[patients_train] - var_names = ['var_datacleaning', 'var_normalization', 'var_fSetReduction'] + var_names = ['cleaning_profile', 'normalization', 'reduction_method'] flags_preprocessing = {key: key in ml['variables'][var_id].keys() for key in var_names} flags_preprocessing_test = flags_preprocessing.copy() - flags_preprocessing_test['var_fSetReduction'] = False + flags_preprocessing_test['reduction_method'] = False # Pre-processing rad_var_struct = ml['variables'][var_id] @@ -179,8 +179,8 @@ def pre_process_radiomics_table( rad_table_learning = get_radiomics_table(path_radiomics_csv, path_radiomics_txt, image_type, patient_ids) # Data cleaning - if flags_preprocessing['var_datacleaning']: - cleaning_dict = ml['datacleaning'][ml['variables'][var_id]['var_datacleaning']]['continuous'] + if flags_preprocessing['cleaning_profile']: + cleaning_dict = ml['datacleaning'][ml['variables'][var_id]['cleaning_profile']]['continuous'] data_cleaner = DataCleaner(**cleaning_dict) # Temp save of properties @@ -196,8 +196,8 @@ def pre_process_radiomics_table( continue # Normalization (ComBat) - if flags_preprocessing['var_normalization']: - normalization_method = ml['variables'][var_id]['var_normalization'] + if flags_preprocessing['normalization']: + normalization_method = ml['variables'][var_id]['normalization'] # Some information must be stored to re-apply combat for testing data if 'combat' in normalization_method.lower(): # Training data @@ -207,10 +207,10 @@ def pre_process_radiomics_table( rad_table_learning.Properties['userData']['normalization']['original_data']['path_radiomics_txt'] = path_radiomics_txt rad_table_learning.Properties['userData']['normalization']['original_data']['image_type'] = image_type rad_table_learning.Properties['userData']['normalization']['original_data']['patient_ids'] = patient_ids - if flags_preprocessing['var_datacleaning']: - data_cln_method = ml['variables'][var_id]['var_datacleaning'] + if flags_preprocessing['cleaning_profile']: + data_cln_method = ml['variables'][var_id]['cleaning_profile'] rad_table_learning.Properties['userData']['normalization']['original_data']['datacleaning_method'] = data_cln_method - + # Apply ComBat normalization = CombatNormalization() rad_table_learning = normalization.fit_transform(rad_table_learning) # Training data @@ -233,8 +233,8 @@ def pre_process_radiomics_table( temp_properties.append(deepcopy(rad_tab.Properties)) # Feature set reduction (for training data only) - if flags_preprocessing['var_fSetReduction']: - f_set_reduction_method = ml['variables'][var_id]['var_fSetReduction']['method'] + if flags_preprocessing['reduction_method']: + f_set_reduction_method = ml['variables'][var_id]['reduction_method'] fsr = FSR(f_set_reduction_method) # Apply FDA @@ -324,13 +324,13 @@ def ml_run(self, path_ml: Path, holdout_test: bool = True, method: str = 'auto') var_table_train = processed_training_table.loc[patients_train, :] # Initializing the model settings - algorithm = ml['settings']['algorithm'] - var_importance_threshold = ml['algorithms'][algorithm]['varImportanceThreshold'] - optimal_threshold = ml['algorithms'][algorithm]['optimalThreshold'] - optimization_metric = ml['algorithms'][algorithm]['optimizationMetric'] - method = ml['algorithms'][algorithm]['method'] if 'method' in ml['algorithms'][algorithm].keys() else method - use_gpu = ml['algorithms'][algorithm]['useGPU'] if 'useGPU' in ml['algorithms'][algorithm].keys() else True - seed = ml['algorithms'][algorithm]['seed'] if 'seed' in ml['algorithms'][algorithm].keys() else None + algorithm = ml['modeling']['method'] if 'method' in ml['modeling'].keys() else method + var_importance_threshold = ml['modeling']['var_importance_threshold'] + optimize_threshold = ml['modeling']['optimize_threshold'] + optimization_metric = ml['modeling']['optimization_metric'] + method = ml['modeling']['method'] if 'method' in ml['modeling'].keys() else method + use_gpu = ml['modeling']['useGPU'] if 'useGPU' in ml['modeling'].keys() else True + seed = ml['modeling']['seed'] if 'seed' in ml['modeling'].keys() else None # B.2. Training the model tstart = time.time() @@ -341,7 +341,7 @@ def ml_run(self, path_ml: Path, holdout_test: bool = True, method: str = 'auto') algorithm=algorithm, ml_config={ 'var_importance_threshold': var_importance_threshold, - 'optimal_threshold': optimal_threshold, + 'optimize_threshold': optimize_threshold, 'optimization_metric': optimization_metric, 'use_gpu': use_gpu, 'seed': seed @@ -349,7 +349,7 @@ def ml_run(self, path_ml: Path, holdout_test: bool = True, method: str = 'auto') estimator.fit(var_table_train, outcome_table_binary_train) # Saving the trained model using pickle - name_save_model = ml['algorithms'][algorithm]['nameSave'] + name_save_model = ml['modeling']['nameSave'] if 'nameSave' in ml['modeling'].keys() else None model_id = name_save_model + '_' + str(ml['variables']['varStudy']) path_model = os.path.dirname(path_results) + '/' + (model_id + '.pickle') estimator.save(path_model) diff --git a/MEDiml/utils/rf_learner.py b/MEDiml/utils/rf_learner.py index 23faf7b..42db3d3 100644 --- a/MEDiml/utils/rf_learner.py +++ b/MEDiml/utils/rf_learner.py @@ -15,14 +15,14 @@ def __init__( optimization_metric='MCC', var_importance_threshold=0.05, internal_cv_folds=5, - optimal_threshold=None, + optimize_threshold=None, use_gpu=False, seed=None ): self.optimization_metric = optimization_metric self.var_importance_threshold = var_importance_threshold self.internal_cv_folds = internal_cv_folds - self.optimal_threshold = optimal_threshold + self.optimize_threshold = optimize_threshold self.use_gpu = use_gpu self.seed = seed @@ -103,10 +103,13 @@ def _train_logic(self, var_table_train, outcome_table_binary_train): model_rf['type'] = 'binary' # Find threshold - try: - model_rf['threshold'] = self.__find_balanced_threshold(classifier, var_table_train, outcome_table_binary_train) - except Exception as e: - print(f'Threshold calculation failed: {e}. Defaulting to 0.5') + if self.optimize_threshold: + try: + model_rf['threshold'] = self.__find_balanced_threshold(classifier, var_table_train, outcome_table_binary_train) + except Exception as e: + print(f'Threshold calculation failed: {e}. Defaulting to 0.5') + model_rf['threshold'] = 0.5 + else: model_rf['threshold'] = 0.5 # Store metadata safely diff --git a/MEDiml/utils/xgboost_learner.py b/MEDiml/utils/xgboost_learner.py index eba3e1f..31160ae 100644 --- a/MEDiml/utils/xgboost_learner.py +++ b/MEDiml/utils/xgboost_learner.py @@ -15,7 +15,7 @@ def __init__( optimization_metric='MCC', var_importance_threshold=0.05, internal_cv_folds=5, - optimal_threshold=None, + optimize_threshold=None, use_gpu=False, seed=None ): @@ -23,7 +23,7 @@ def __init__( self.optimization_metric = optimization_metric self.var_importance_threshold = var_importance_threshold self.internal_cv_folds = internal_cv_folds - self.optimal_threshold = optimal_threshold + self.optimize_threshold = optimize_threshold self.use_gpu = use_gpu self.seed = seed @@ -115,11 +115,17 @@ def _train_logic(self, var_table_train, outcome_table_binary_train): model_xgb = dict() model_xgb['algo'] = 'xgb' model_xgb['type'] = 'binary' - try: - model_xgb['threshold'] = self.__find_balanced_threshold(classifier, var_table_train, outcome_table_binary_train) - except Exception as e: - print('Error in finding optimal threshold, it will be set to 0.5:' + str(e)) + + # Find threshold + if self.optimize_threshold: + try: + model_xgb['threshold'] = self.__find_balanced_threshold(classifier, var_table_train, outcome_table_binary_train) + except Exception as e: + print('Error in finding optimal threshold, it will be set to 0.5:' + str(e)) + model_xgb['threshold'] = 0.5 + else: model_xgb['threshold'] = 0.5 + model_xgb['var_info'] = deepcopy(var_table_train.Properties['userData']) model_xgb['var_def'] = deepcopy(var_table_train.Properties['userData']['variables']['var_def']) model_xgb['var_names'] = list(classifier.feature_names_in_) From 316d35bee936796741df7f7ade74ea7cf7c96ca0 Mon Sep 17 00:00:00 2001 From: MahdiAll99 Date: Wed, 18 Mar 2026 15:59:26 -0400 Subject: [PATCH 12/15] New standardized process to compute SUV values --- MEDiml/MEDscan.py | 5 +- MEDiml/processing/__init__.py | 3 +- MEDiml/processing/compute_suv_map.py | 121 -------------------------- MEDiml/wrangling/DataManager.py | 103 +++++++++++----------- MEDiml/wrangling/ProcessDICOM.py | 122 ++++++++++++++++++++++++--- 5 files changed, 167 insertions(+), 187 deletions(-) delete mode 100644 MEDiml/processing/compute_suv_map.py diff --git a/MEDiml/MEDscan.py b/MEDiml/MEDscan.py index 857a185..9af9f2a 100644 --- a/MEDiml/MEDscan.py +++ b/MEDiml/MEDscan.py @@ -293,8 +293,9 @@ def init_params(self, im_param_scan: Dict) -> None: if self.type == 'PTscan' and _compute_suv_map and self.format != 'nifti': try: - from .processing.compute_suv_map import compute_suv_map - self.data.volume.array = compute_suv_map(self.data.volume.array, self.dicomH[0]) + from .processing.PETSUVConverter import PETSUVConverter + suv_converter = PETSUVConverter(self.dicomH) + self.data.volume.array = suv_converter.compute(self.data.volume.array) except Exception as e : message = f"\n ERROR COMPUTING SUV MAP - SOME FEATURES WILL BE INVALID: \n {e}" logging.error(message) diff --git a/MEDiml/processing/__init__.py b/MEDiml/processing/__init__.py index 32514ea..a7d72ce 100644 --- a/MEDiml/processing/__init__.py +++ b/MEDiml/processing/__init__.py @@ -1,6 +1,7 @@ from . import * -from .compute_suv_map import * from .discretisation import * from .interpolation import * +from .PETSUVConverter import * from .resegmentation import * from .segmentation import * +from .SUVHeaderProxy import * diff --git a/MEDiml/processing/compute_suv_map.py b/MEDiml/processing/compute_suv_map.py deleted file mode 100644 index 48e7783..0000000 --- a/MEDiml/processing/compute_suv_map.py +++ /dev/null @@ -1,121 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - - -import numpy as np -import pydicom - - -def compute_suv_map(raw_pet: np.ndarray, - dicom_h: pydicom.Dataset) -> np.ndarray: - """Computes the suv_map of a raw input PET volume. It is assumed that - the calibration factor was applied beforehand to the PET volume. - **E.g: raw_pet = raw_pet*RescaleSlope + RescaleIntercept.** - - Args: - raw_pet (ndarray):3D array representing the PET volume in raw format. - dicom_h (pydicom.dataset.FileDataset): DICOM header of one of the - corresponding slice of ``raw_pet``. - - Returns: - ndarray: ``raw_pet`` converted to SUVs (standard uptake values). - """ - def dcm_hhmmss(date_str: str) -> float: - """"Converts to seconds - - Args: - date_str (str): date string - - Returns: - float: total seconds - """ - # Converts to seconds - if not isinstance(date_str, str): - date_str = str(date_str) - hh = float(date_str[0:2]) - mm = float(date_str[2:4]) - ss = float(date_str[4:6]) - tot_sec = hh*60.0*60.0 + mm*60.0 + ss - return tot_sec - - def pydicom_has_tag(dcm_seq, tag): - # Checks if tag exists - return get_pydicom_meta_tag(dcm_seq, tag, test_tag=True) - - def get_pydicom_meta_tag(dcm_seq, tag, tag_type=None, default=None, - test_tag=False): - # Reads dicom tag - # Initialise with default - tag_value = default - # Read from header using simple itk - try: - tag_value = dcm_seq[tag].value - except KeyError: - if test_tag: - return False - if test_tag: - return True - # Find empty entries - if tag_value is not None: - if tag_value == "": - tag_value = default - # Cast to correct type (meta tags are usually passed as strings) - if tag_value is not None: - # String - if tag_type == "str": - tag_value = str(tag_value) - # Float - elif tag_type == "float": - tag_value = float(tag_value) - # Multiple floats - elif tag_type == "mult_float": - tag_value = [float(str_num) for str_num in tag_value] - # Integer - elif tag_type == "int": - tag_value = int(tag_value) - # Multiple floats - elif tag_type == "mult_int": - tag_value = [int(str_num) for str_num in tag_value] - # Boolean - elif tag_type == "bool": - tag_value = bool(tag_value) - - return tag_value - - # Get patient weight - if pydicom_has_tag(dcm_seq=dicom_h, tag=(0x0010, 0x1030)): - weight = get_pydicom_meta_tag(dcm_seq=dicom_h, tag=(0x0010, 0x1030), - tag_type="float") * 1000.0 # in grams - else: - weight = None - if weight is None: - weight = 75000.0 # estimation - try: - # Get Scan time - scantime = dcm_hhmmss(date_str=get_pydicom_meta_tag( - dcm_seq=dicom_h, tag=(0x0008, 0x0032), tag_type="str")) - # Start Time for the Radiopharmaceutical Injection - injection_time = dcm_hhmmss(date_str=get_pydicom_meta_tag( - dcm_seq=dicom_h[0x0054, 0x0016][0], - tag=(0x0018, 0x1072), tag_type="str")) - # Half Life for Radionuclide - half_life = get_pydicom_meta_tag( - dcm_seq=dicom_h[0x0054, 0x0016][0], - tag=(0x0018, 0x1075), tag_type="float") - # Total dose injected for Radionuclide - injected_dose = get_pydicom_meta_tag( - dcm_seq=dicom_h[0x0054, 0x0016][0], - tag=(0x0018, 0x1074), tag_type="float") - # Calculate decay - decay = np.exp(-np.log(2)*(scantime-injection_time)/half_life) - # Calculate the dose decayed during procedure - injected_dose_decay = injected_dose*decay # in Bq - except KeyError: - # 90 min waiting time, 15 min preparation - decay = np.exp(-np.log(2)*(1.75*3600)/6588) - injected_dose_decay = 420000000 * decay # 420 MBq - - # Calculate SUV - suv_map = raw_pet * weight / injected_dose_decay - - return suv_map diff --git a/MEDiml/wrangling/DataManager.py b/MEDiml/wrangling/DataManager.py index efa6b13..37eedeb 100644 --- a/MEDiml/wrangling/DataManager.py +++ b/MEDiml/wrangling/DataManager.py @@ -13,7 +13,6 @@ import numpy as np import pandas as pd import pydicom -import pydicom.errors import pydicom.misc import ray from nilearn import image @@ -21,7 +20,7 @@ from tqdm import tqdm, trange from ..MEDscan import MEDscan -from ..processing.compute_suv_map import compute_suv_map +from ..processing.PETSUVConverter import PETSUVConverter from ..processing.segmentation import get_roi_from_indexes from ..utils.get_file_paths import get_file_paths from ..utils.get_patient_names import get_patient_names @@ -294,7 +293,10 @@ def process_all_dicoms(self) -> Union[List[MEDscan], None]: Returns: List[MEDscan]: List of MEDscan instances. """ - ray.init(local_mode=True, include_dashboard=True) + # Initialize ray + if ray.is_initialized(): + ray.shutdown() + ray.init(local_mode=True, include_dashboard=False) print('--> Reading all DICOM objects to create MEDscan classes') self.__read_all_dicoms() @@ -319,50 +321,10 @@ def process_all_dicoms(self) -> Union[List[MEDscan], None]: ids = [pd.process_files() for pd in pds] # Update the path to the created instances - for name_save in ray.get(ids): - if self.paths._path_save: - self.path_to_objects.append(str(self.paths._path_save / name_save)) - # Update processing summary - if name_save.split('_')[0].count('-') >= 2: - scan_type = name_save[name_save.find('__')+2 : name_save.find('.')] - if name_save.split('-')[0] not in self.__studies: - self.__studies.append(name_save.split('-')[0]) # add new study - if name_save.split('-')[1] not in self.__institutions: - self.__institutions.append(name_save.split('-')[1]) # add new study - if name_save.split('-')[0] not in self.summary: - self.summary[name_save.split('-')[0]] = {} - if name_save.split('-')[1] not in self.summary[name_save.split('-')[0]]: - self.summary[name_save.split('-')[0]][name_save.split('-')[1]] = {} # add new institution - if scan_type not in self.__scans: - self.__scans.append(scan_type) - if scan_type not in self.summary[name_save.split('-')[0]][name_save.split('-')[1]]: - self.summary[name_save.split('-')[0]][name_save.split('-')[1]][scan_type] = [] - if name_save not in self.summary[name_save.split('-')[0]][name_save.split('-')[1]][scan_type]: - self.summary[name_save.split('-')[0]][name_save.split('-')[1]][scan_type].append(name_save) - else: - if self.save: - logging.warning(f"The patient ID of the following file: {name_save} does not respect the MEDiml "\ - "naming convention 'study-institution-id' (Ex: Glioma-TCGA-001)") - - nb_job_left = n_scans - n_batch - - # Distribute the remaining tasks - for _ in trange(n_scans): - _, ids = ray.wait(ids, num_returns=1) - if nb_job_left > 0: - idx = n_scans - nb_job_left - pd = ProcessDICOM( - self.__dicom.cell_path_images[idx], - self.__dicom.cell_path_rs[idx], - self.paths._path_save, - self.save) - ids.extend([pd.process_files()]) - nb_job_left -= 1 - - # Update the path to the created instances + if self.save: for name_save in ray.get(ids): if self.paths._path_save: - self.path_to_objects.extend(str(self.paths._path_save / name_save)) + self.path_to_objects.append(str(self.paths._path_save / name_save)) # Update processing summary if name_save.split('_')[0].count('-') >= 2: scan_type = name_save[name_save.find('__')+2 : name_save.find('.')] @@ -381,9 +343,53 @@ def process_all_dicoms(self) -> Union[List[MEDscan], None]: if name_save not in self.summary[name_save.split('-')[0]][name_save.split('-')[1]][scan_type]: self.summary[name_save.split('-')[0]][name_save.split('-')[1]][scan_type].append(name_save) else: - if self.save: + logging.warning(f"The patient ID of the following file: {name_save} does not respect the MEDiml "\ + "naming convention 'study-institution-id' (Ex: Glioma-TCGA-001)") + + nb_job_left = n_scans - n_batch + + return ray.get(ids) if not self.save else None + + # Distribute the remaining tasks + for _ in trange(n_scans): + _, ids = ray.wait(ids, num_returns=1) + if nb_job_left > 0: + idx = n_scans - nb_job_left + pd = ProcessDICOM( + self.__dicom.cell_path_images[idx], + self.__dicom.cell_path_rs[idx], + self.paths._path_save, + self.save) + ids.extend([pd.process_files()]) + nb_job_left -= 1 + + # Update the path to the created instances + if self.save: + for name_save in ray.get(ids): + if self.paths._path_save: + self.path_to_objects.extend(str(self.paths._path_save / name_save)) + # Update processing summary + if name_save.split('_')[0].count('-') >= 2: + scan_type = name_save[name_save.find('__')+2 : name_save.find('.')] + if name_save.split('-')[0] not in self.__studies: + self.__studies.append(name_save.split('-')[0]) # add new study + if name_save.split('-')[1] not in self.__institutions: + self.__institutions.append(name_save.split('-')[1]) # add new study + if name_save.split('-')[0] not in self.summary: + self.summary[name_save.split('-')[0]] = {} + if name_save.split('-')[1] not in self.summary[name_save.split('-')[0]]: + self.summary[name_save.split('-')[0]][name_save.split('-')[1]] = {} # add new institution + if scan_type not in self.__scans: + self.__scans.append(scan_type) + if scan_type not in self.summary[name_save.split('-')[0]][name_save.split('-')[1]]: + self.summary[name_save.split('-')[0]][name_save.split('-')[1]][scan_type] = [] + if name_save not in self.summary[name_save.split('-')[0]][name_save.split('-')[1]][scan_type]: + self.summary[name_save.split('-')[0]][name_save.split('-')[1]][scan_type].append(name_save) + else: logging.warning(f"The patient ID of the following file: {name_save} does not respect the MEDiml "\ "naming convention 'study-institution-id' (Ex: Glioma-TCGA-001)") + else: + return ray.get(ids) print('DONE') def __read_all_niftis(self) -> None: @@ -1012,9 +1018,8 @@ def __pre_radiomics_checks_window( with open(file, 'rb') as file: medscan = pickle.load(file) if re.search('PTscan', wildcard) and medscan.format != 'nifti': - medscan.data.volume.array = compute_suv_map( - np.double(medscan.data.volume.array), - medscan.dicomH[2]) + suv_converter = PETSUVConverter(medscan.dicomH) + medscan.data.volume.array = suv_converter.compute(np.double(medscan.data.volume.array)) patient_names = pd.Index(patient_names) ind_roi = patient_names.get_loc(patient_name) name_roi = roi_table.loc[ind_roi][3] diff --git a/MEDiml/wrangling/ProcessDICOM.py b/MEDiml/wrangling/ProcessDICOM.py index afd7269..368f8da 100644 --- a/MEDiml/wrangling/ProcessDICOM.py +++ b/MEDiml/wrangling/ProcessDICOM.py @@ -76,6 +76,81 @@ def __get_dicom_scan_orientation(self, dicom_header: List[pydicom.dataset.FileDa return orientation + def __get_minimal_suv_header(self, dcm: pydicom.Dataset) -> dict: + """ + Extracts only the tags required for SUV conversion and PET scaling. + This dict is Ray-serializable and free of weakrefs. + """ + def find_philips_private_tags(ds): + # Find which block Philips reserved + offset = None + for i in range(0x10, 0x100, 0x01): + tag = (0x7053, i) + if tag in ds and ds[tag].value.lower().startswith("philips"): + # If (7053, 0011) is the creator, the offset is 0x1100 + offset = i << 8 + break + + if offset: + suv_tag = (0x7053, offset + 0x00) + act_tag = (0x7053, offset + 0x09) + print(f"Philips Tags Found at: SUV={hex(suv_tag[1])}, Act={hex(act_tag[1])}") + return ds.get(suv_tag), ds.get(act_tag) + + print("Creator 'Philips PET Private Group' not found in group 0x7053.") + return None, None + + suv_elem, act_elem = find_philips_private_tags(dcm) + # Map the tags to their values (storing as hex strings for keys) + tags = [ + 0x00101030, 0x00100040, 0x00080031, 0x00080032, 0x00080021, 0x00541102, + 0x00181072, 0x00181078, 0x00281052, 0x00281053, 0x00080070, 0x00541001, + 0x00541001, 0x00541006, 0x00101020, 0x00101040, 0x00280030, 0x00180050, + ] + + suv_data = {} + for t in tags: + if t in dcm: + suv_data[t] = dcm[t].value + elif t == 0x00541006 and 0x00541001 in dcm and dcm[0x00541001].value == 'GML': + suv_data[t] = 'BW' # If absent, and the Units are GML, then the type of SUV shall be assumed to be BW. + + # Handle the Radiopharmaceutical Sequence specially + radio_tag = 0x00540016 + if radio_tag in dcm and dcm[radio_tag]: + item = dcm[radio_tag][0] + sub_tags = [0x00181072, 0x00181074, 0x00181075, 0x00181078] + radio_dict = {st: item[st].value for st in sub_tags if st in item} + suv_data[radio_tag] = [radio_dict] # List of dicts + + # Extra tags depending on the unit type + if 0x00541001 in dcm: + unit = str(dcm[0x00541001].value).lower() + if unit == 'cnts': + # SD SUV scale factor + if 0x70531000 in dcm: + suv_data[0x70531000] = dcm[0x70531000].value + # If not found, try DS Activity Concentration Scale Factor + elif 0x70531009 in dcm: + suv_data[0x70531000] = dcm[0x70531009].value + # If still not found, try Frame Duration (for dynamic PET) + elif 0x00181242 in dcm: + suv_data[0x00181242] = dcm[0x00181242].value + # Dose Calibration Factor (Needed to convert to CPS then to BQML) + if 0x00541322 in dcm: + suv_data[0x00541322] = dcm[0x00541322].value + # Corrected image tag + if 0x00280051 in dcm: + suv_data[0x00280051] = dcm[0x00280051].value + elif unit == 'cps': + # Doe Calibration Factor + if 0x00541322 in dcm: + suv_data[0x00541322] = dcm[0x00541322].value + # Corrected image tag + if 0x00280051 in dcm: + suv_data[0x00280051] = dcm[0x00280051].value + return suv_data + def __merge_slice_pixel_arrays(self, slice_datasets): first_dataset = slice_datasets[0] num_rows = first_dataset.Rows @@ -84,14 +159,16 @@ def __merge_slice_pixel_arrays(self, slice_datasets): sorted_slice_datasets = self.__sort_by_slice_spacing(slice_datasets) - if any(self.__requires_rescaling(d) for d in sorted_slice_datasets): - voxels = np.empty( - (num_columns, num_rows, num_slices), dtype=np.float32) + if self.__requires_rescaling(sorted_slice_datasets): + if not self.__rescaling_is_the_same(sorted_slice_datasets): + if not self.__intercept_is_zero(sorted_slice_datasets): + # The scan is skipped is the RescaleSlope attribute has multiple values across slices and the RescaleIntercept is not zero. + return None + voxels = np.empty((num_columns, num_rows, num_slices), dtype=np.float32) for k, dataset in enumerate(sorted_slice_datasets): slope = float(getattr(dataset, 'RescaleSlope', 1)) intercept = float(getattr(dataset, 'RescaleIntercept', 0)) - voxels[:, :, k] = dataset.pixel_array.T.astype( - np.float32)*slope + intercept + voxels[:, :, k] = dataset.pixel_array.T.astype(np.float32) * slope + intercept else: dtype = first_dataset.pixel_array.dtype voxels = np.empty((num_columns, num_rows, num_slices), dtype=dtype) @@ -100,8 +177,23 @@ def __merge_slice_pixel_arrays(self, slice_datasets): return voxels - def __requires_rescaling(self, dataset): - return hasattr(dataset, 'RescaleSlope') or hasattr(dataset, 'RescaleIntercept') + def __requires_rescaling(self, slice_datasets): + return any(hasattr(dataset, 'RescaleSlope') or hasattr(dataset, 'RescaleIntercept') for dataset in slice_datasets) + + def __rescaling_is_the_same(self, slice_datasets): + first_slope = float(getattr(slice_datasets[0], 'RescaleSlope', 1)) + for dataset in slice_datasets[1:]: + slope = float(getattr(dataset, 'RescaleSlope', 1)) + if slope != first_slope: + return False + return True + + def __intercept_is_zero(self, slice_datasets): + for dataset in slice_datasets: + intercept = float(getattr(dataset, 'RescaleIntercept', 0)) + if intercept != 0: + return False + return True def __ijk_to_patient_xyz_transform_matrix(self, slice_datasets): first_dataset = self.__sort_by_slice_spacing(slice_datasets)[0] @@ -299,6 +391,8 @@ def combine_slices(self, slice_datasets: List[pydicom.dataset.FileDataset]) -> L self.__validate_slices_form_uniform_grid(slice_datasets) voxels = self.__merge_slice_pixel_arrays(slice_datasets) + if voxels is None: + return None, None, None, None transform, rotation, scaling = self.__ijk_to_patient_xyz_transform_matrix( slice_datasets) @@ -345,6 +439,8 @@ def process_files_wrapper(self) -> MEDscan: # https://dicom-numpy.readthedocs.io/en/latest/index.html#dicom_numpy.combine_slices try: voxel_ndarray, ijk_to_xyz, rotation_m, scaling_m = self.combine_slices(dicom_hi) + if voxel_ndarray is None: + return None except ValueError as e: raise ValueError(f'Invalid DICOM data for combine_slices(). Error: {e}') @@ -387,11 +483,13 @@ def process_files_wrapper(self) -> MEDscan: # DICOM HEADERS OF IMAGING DATA dicom_h = [ - pydicom.dcmread(str(dicom_file),stop_before_pixels=True,force=True) for dicom_file in self.path_images + pydicom.dcmread(str(dicom_file),stop_before_pixels=True) for dicom_file in self.path_images ] for i in range(0, len(dicom_h)): dicom_h[i].remove_private_tags() - medscan.dicomH = dicom_h + + # Save the minimal header required for SUV conversion and PET scaling in the MEDscan class + medscan.dicomH = self.__get_minimal_suv_header(dicom_h[0]) # DICOM RTstruct (if applicable) if self.path_rs is not None and len(self.path_rs) > 0: @@ -496,11 +594,7 @@ def process_files_wrapper(self) -> MEDscan: name_complete = save_MEDscan(medscan, self.path_save) del medscan else: - series_description = medscan.series_description.translate({ord(ch): '-' for ch in '/\\ ()&:*'}) - name_id = medscan.patientID.translate({ord(ch): '-' for ch in '/\\ ()&:*'}) - - # final saving name - name_complete = name_id + '__' + series_description + '.' + medscan.type + '.npy' + return medscan except Exception as e: if 'SeriesDescription' in dicom_hi[0]: From 775af691cacc7c7c1511c909e3835d720b322d8a Mon Sep 17 00:00:00 2001 From: MahdiAll99 Date: Wed, 18 Mar 2026 16:13:11 -0400 Subject: [PATCH 13/15] New standardized process to compute SUV values (init commit) --- MEDiml/processing/PETSUVConverter.py | 229 +++++++++++++++++++++++++++ 1 file changed, 229 insertions(+) create mode 100644 MEDiml/processing/PETSUVConverter.py diff --git a/MEDiml/processing/PETSUVConverter.py b/MEDiml/processing/PETSUVConverter.py new file mode 100644 index 0000000..77e8667 --- /dev/null +++ b/MEDiml/processing/PETSUVConverter.py @@ -0,0 +1,229 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import logging + +import dateutil.parser +import numpy as np + + +class PETSUVConverter: + """ + A class for converting raw PET volumes into Standardized Uptake Value (SUV) maps. + """ + + def __init__(self, dicom_proxy): + """ + Initializes the converter with a DICOM proxy object. + """ + self.dcm = dicom_proxy + self.logger = logging.getLogger(self.__class__.__name__) + + # Strategy pattern for dynamic computation routing + self._strategies = { + 'gml': self._compute_gml, + 'bqml': self._compute_bqml, + 'cm2ml': self._compute_cm2ml, + 'cnts': self._compute_cnts, + 'cps': self._compute_cps + } + + # ========================================== # + # PROPERTIES # + # ========================================== # + + @property + def unit(self) -> str: + return str(self.dcm.get(0x00541001, 'unknown')).lower() + + @property + def patient_weight_g(self) -> float: + """Returns patient weight_kg in grams.""" + return float(self.dcm[0x0010, 0x1030].value) * 1000.0 if (0x0010, 0x1030) in self.dcm else 75000.0 + + @property + def patient_height_cm(self) -> float: + """Returns patient height_m in cm.""" + return float(self.dcm[0x0010, 0x1020].value) * 100.0 if (0x0010, 0x1020) in self.dcm else 170.0 + + @property + def patient_sex(self) -> str: + return str(self.dcm[0x0010, 0x0040].value).upper() if (0x0010, 0x0040) in self.dcm else 'O' + + # ========================================== # + # PUBLIC METHODS # + # ========================================== # + + def get_conversion_factors(self) -> tuple: + """ + Calculates conversion factors (suv_factor, rescale_slope, rescale_intercept). + """ + nuclide_dose = self.dcm[0x0054, 0x0016][0][0x0018, 0x1074].value + weight_kg = self.patient_weight_g / 1000.0 + half_life = float(self.dcm[0x0054, 0x0016][0][0x0018, 0x1075].value) + + series_time = str(self.dcm[0x0008, 0x0031].value) + series_date = str(self.dcm[0x0008, 0x0021].value) + series_dt = dateutil.parser.parse(f"{series_date} {series_time}") + + nuclide_time = str(self.dcm[0x0054, 0x0016][0][0x0018, 0x1072].value) + nuclide_dt = dateutil.parser.parse(f"{series_date} {nuclide_time}") + + delta_time = (series_dt - nuclide_dt).total_seconds() + decay_correction = 2 ** (-1 * delta_time / half_life) + suv_factor = (weight_kg * 1000) / (decay_correction * nuclide_dose) + + rescale_slope = self.dcm[0x0028, 0x1053].value + rescale_intercept = self.dcm[0x0028, 0x1052].value + + manufacturer = str(self.dcm[0x0008, 0x0070].value).lower() + + # Philips private tag logic + if "philips" in manufacturer and "bqml" not in self.unit: + if 0x70531000 in self.dcm: + suv_factor = float(self.dcm[0x7053, 0x1000].value) + + return (suv_factor, rescale_slope, rescale_intercept) + + def compute(self, raw_pet: np.ndarray) -> np.ndarray: + """ + Main entry point. Routes the raw PET array to the correct computation strategy based on unit. + """ + strategy = self._strategies.get(self.unit) + if not strategy: + raise ValueError(f"Unsupported unit '{self.unit}' for SUV computation.") + + return strategy(raw_pet) + + # ========================================== # + # COMPUTATION STRATEGIES # + # ========================================== # + + def _compute_gml(self, raw_pet: np.ndarray) -> np.ndarray: + weight_kg = self.patient_weight_g / 1000.0 + sex = self.patient_sex + suv_type = self.dcm.get(0x00541006, 'unknown').lower() + + lbm = None + bmi = weight_kg / (self.patient_height_cm ** 2) if self.patient_height_cm > 0 else 0 + + if suv_type == 'lbm': + if sex == 'M': + lbm = 1.10 * weight_kg - 120 * (weight_kg / self.patient_height_cm) ** 2 + elif sex == 'F' or sex == 'O': + lbm = 1.07 * weight_kg - 148 * (weight_kg / self.patient_height_cm) ** 2 + + elif suv_type == 'lbmjames128': + if sex == 'M': + lbm = 1.10 * weight_kg - 128 * (weight_kg / self.patient_height_cm) ** 2 + elif sex == 'F' or sex == 'O': + lbm = 1.07 * weight_kg - 148 * (weight_kg / self.patient_height_cm) ** 2 + + elif suv_type == 'lbmjamna': + if sex == 'M': + lbm = 9270 * weight_kg / (6680 + 216 * bmi) + else: + lbm = 9270 * weight_kg / (8780 + 244 * bmi) + + elif suv_type == 'ibw': + if sex == 'M': + lbm = 48 + 1.06 * (self.patient_height_cm - 152) + else: + lbm = 45.5 + 0.91 * (self.patient_height_cm - 152) + + elif suv_type == 'bw': + return raw_pet + + else: + raise ValueError(f"Unsupported SUV type '{suv_type}'.") + + return (raw_pet / lbm) * weight_kg + + def _compute_cm2ml(self, raw_pet: np.ndarray) -> np.ndarray: + weight_kg = self.patient_weight_g / 1000.0 + bsa = 0.007184 * (weight_kg ** 0.425) * (self.patient_height_cm ** 0.725) + return (raw_pet / bsa) * weight_kg / 10 + + def _compute_cnts(self, raw_pet: np.ndarray) -> np.ndarray: + if 0x70531000 in self.dcm and float(self.dcm.get(0x70531000)) != 0: + return raw_pet * float(self.dcm.get(0x70531000)) + + if 0x70531009 in self.dcm and float(self.dcm.get(0x70531009)) != 0: + act_scale = float(self.dcm.get(0x70531009)) + return self._compute_bqml(raw_pet * act_scale) + + if 0x00181242 in self.dcm and float(self.dcm.get(0x00181242)) != 0: + frame_duration_sec = float(self.dcm.get(0x00181242)) / 1000.0 + return self._compute_cps(raw_pet / frame_duration_sec) + + raise ValueError("No valid scale factor found for 'cnts' unit in DICOM header.") + + def _compute_cps(self, cps_map: np.ndarray) -> np.ndarray: + corrected_image_tags = self.dcm.get(0x00280051) if 0x00280051 in self.dcm else [] + is_dcal = "DCAL" in corrected_image_tags + + pixel_spacing = self.dcm.get(0x00280030) + slice_thickness = self.dcm.get(0x00180050) + + if not pixel_spacing or not slice_thickness: + raise KeyError("Voxel dimensions (0028,0030 or 0018,0050) missing.") + + voxel_vol_ml = (float(pixel_spacing[0]) * float(pixel_spacing[1]) * float(slice_thickness)) / 1000.0 + + if is_dcal: + bqml_map = cps_map / voxel_vol_ml + else: + cal_factor = self.dcm.get(0x00541322) + if cal_factor is None: + raise ValueError("Image is not DCAL and Dose Calibration Factor (0054,1322) is unknown.") + bqml_map = (cps_map * float(cal_factor)) / voxel_vol_ml + + return self._compute_bqml(bqml_map) + + def _compute_bqml(self, raw_pet: np.ndarray) -> np.ndarray: + try: + scantime = self._parse_time(str(self.dcm[0x0008, 0x0032].value)) + radio_item = self.dcm[0x0054, 0x0016][0] + + if (0x0018, 0x1072) not in radio_item and (0x0018, 0x1078) in radio_item: + injection_time = self._parse_time(str(radio_item[0x0018, 0x1078].value)[8:]) + elif (0x0018, 0x1072) in radio_item: + injection_time = self._parse_time(str(radio_item[0x0018, 0x1072].value)) + else: + raise KeyError("Radiopharmaceutical Start Time tags missing.") + + half_life = float(radio_item[0x0018, 0x1075].value) + injected_dose = float(radio_item[0x0018, 0x1074].value) + decay_correction = str(self.dcm.get(0x00541102, '')).upper() + + if decay_correction == 'ADMIN': + injected_dose_decay = injected_dose + elif decay_correction in ['START', 'NONE']: + decay = np.exp(-np.log(2) * (scantime - injection_time) / half_life) + injected_dose_decay = injected_dose * decay + else: + raise ValueError(f"Unrecognized decay correction status: {decay_correction}") + + raw_pet = raw_pet * self.patient_weight_g / injected_dose_decay + + # Convert MBq to Bq if necessary + if (np.any(raw_pet > 0) and np.nanmean(raw_pet[raw_pet > 0]) > 100): + raw_pet = raw_pet / 1_000_000.0 + + except Exception as e: + self.logger.warning(f"Error computing BQML ({e}). Using standard 1.75h fallback decay.") + decay = np.exp(-np.log(2) * (1.75 * 3600) / 6588) + raw_pet = raw_pet * self.patient_weight_g / (420000000 * decay) + + return raw_pet + + # ========================================== # + # STATIC HELPERS # + # ========================================== # + + @staticmethod + def _parse_time(time_str: str) -> float: + """Helper to convert HHMMSS string to total seconds.""" + time_str = str(time_str).zfill(6) + hh, mm, ss = float(time_str[0:2]), float(time_str[2:4]), float(time_str[4:6]) + return hh * 3600.0 + mm * 60.0 + ss From 26817213d0434c47b50c49f892b31909a403f4d2 Mon Sep 17 00:00:00 2001 From: MahdiAll99 Date: Wed, 18 Mar 2026 16:13:32 -0400 Subject: [PATCH 14/15] docs update: yaml config documentation added to readthedocs --- docs/learning_config.rst | 692 +++++++++++---------------------------- docs/processing.rst | 4 +- 2 files changed, 194 insertions(+), 502 deletions(-) diff --git a/docs/learning_config.rst b/docs/learning_config.rst index b65b6e6..8744ff5 100644 --- a/docs/learning_config.rst +++ b/docs/learning_config.rst @@ -1,528 +1,220 @@ Learning -------- -This section will walk you through the details on how to set up the configuration file for the machine learning part of the pipeline. -It will be separated to the following subdivisions: +This section walks you through setting up the consolidated master configuration file (``config.yaml``) for the machine learning pipeline. +Instead of multiple JSON files in the previous versions, all parameters are now managed in a single YAML structure, +allowing for easier maintenance and the use of YAML anchors for consistency (e.g., shared seeds). -- :ref:`Design` -- :ref:`Data Cleaning` -- :ref:`Data Normalization` -- :ref:`Feature Set Reduction` -- :ref:`Machine Learning` -- :ref:`Variables Definition` +The configuration is separated into the following subdivisions: + +* :ref:`Study Metadata` +* :ref:`Experiment Design Parameters` +* :ref:`Variables Definition` +* :ref:`Data Cleaning Parameters` +* :ref:`Data Normalization Parameters` +* :ref:`Feature Set Reduction Parameters` +* :ref:`Machine Learning Parameters` + +Study Metadata +^^^^^^^^^^^^^^ + +Defines high-level experiment identifiers and global variables. + +.. code-block:: yaml + + study_metadata: + var_study: "var1" + combinations: ["var1"] + seed: &global_seed 54288 # YAML anchor used to sync seeds across the pipeline Experiment Design Parameters ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -This set of parameters is used to define the experiment design (data splitting, splitting proportion...), it is organized as follows: - -.. code-block:: JSON - - { - "testSets": ["Define method here"], - "method name": "Define method here" - - } - -Now let's specify the parameters for the selected method; for instance, in the case of the ``Random`` and ``CV`` methods: - -.. jsonschema:: - - { - "$schema": "http://json-schema.org/draft-04/schema#", - "title": "Splitting methods", - "description": "Type of sets to create.", - "type": "object", - "properties": { - "Random": { - "description": "Random splitting method.", - "type": "object", - "properties": { - "method": { - "description": "Method of splitting the data.", - "type": "string", - "options": { - "SubSampling": { - "description": "The data will be randomly split", - "type": "string" - }, - "Institutions": { - "description": "The data will be split based on institutions", - "type": "string" - } - } - }, - "nSplits": { - "description": "Number of splits to create.", - "type": "int" - }, - "stratifyInstitutions": { - "description": "If ``True``, the data will be stratified based on institutions.", - "type": "bool" - }, - "testProportion": { - "description": "Proportion of the test set.", - "type": "float" - }, - "seed": { - "description": "Seed for the random number generator.", - "type": "int" - } - } - }, - "CV" : { - "description": "Cross-validation splitting method.", - "type": "object", - "properties": { - "nFolds": { - "description": "Number of folds to use.", - "type": "int" - }, - "seed": { - "description": "Seed for the random number generator.", - "type": "int" - } - } - } - } - } - -- **Example** - -.. code-block:: JSON - - { - "Random": { - "method": "SubSampling", - "nSplits": 10, - "stratifyInstitutions": 1, - "testProportion": 0.33, - "seed": 54288 - } - } +Used to define the data splitting and validation strategy. You can define multiple profiles and select the active one. + +.. code-block:: yaml + + design: + active_method: "CrossValidation" # Options: "Random", "CrossValidation", "Bootstrapping" + + Random: + method: "SubSampling" + nSplits: 10 + stratifyInstitutions: 1 + testProportion: 0.33 + seed: *global_seed + + CrossValidation: + method: "StratifiedKFold" + nFolds: 5 + nRepeats: 10 + seed: *global_seed + + Bootstrapping: + method: "Out-of-Bag" + nIterations: 1000 + seed: *global_seed + +Variables Definition +^^^^^^^^^^^^^^^^^^^^ + +Defines the data sources and maps them to specific cleaning and reduction profiles defined later in the file. + +.. code-block:: yaml + variables: + var1: + nameType: "RadiomicsFull" + path: "setToFeaturesinWorkspace" + scans: ["CECT"] + rois: ["tumor"] + imSpaces: ["image"] + cleaning_profile: "default" + normalization: "combat" + reduction_method: "FDA" Data Cleaning Parameters ^^^^^^^^^^^^^^^^^^^^^^^^ -This set of parameters is used to define the data cleaning process Parameters, it is organized as follows: - -.. code-block:: JSON - - { - "method name": { - "define parameters here" - }, - "another method": { - "define parameters here" - } - } - -.. jsonschema:: - - { - "$schema": "http://json-schema.org/draft-04/schema#", - "title": "Cleaning methods", - "description": "Feature cleaning method name.", - "type": "object", - "properties": { - "default": { - "description": "Default cleaning method.", - "type": "string" - } - } - } - -Now let's specify the parameters for the selected cleaning method; for instance, in the case of the ``default`` method: - -.. jsonschema:: - - { - "$schema": "http://json-schema.org/draft-04/schema#", - "title": "Chosen method's parameters", - "description": "Feature cleaning parameters.", - "type": "object", - "properties": { - "continuous": { - "description": "Continuous feature cleaning parameters.", - "type": "object", - "properties": { - "missingCutoffps": { - "description": "Maximum percentage cut-offs of missing features per sample. Samples with more missing features than this cut-off will be removed.", - "type": "float" - }, - "covCutoff": { - "description": "Minimal coefficient of variation cut-offs over samples per variable. Variables with less coefficient of variation than this cut-off will be removed.", - "type": "float" - }, - "missingCutoffpf": { - "description": "Maximal percentage cut-offs of missing samples per variable. Features with more missing samples than this cut-off will be removed.", - "type": "float" - }, - "imputation": { - "description": "Imputation method for missing values. Default is ``mean``.", - "type": "string", - "options": { - "mean": { - "description": "Impute missing values with the mean of the feature.", - "type": "string" - }, - "median": { - "description": "Impute missing values with the median of the feature.", - "type": "string" - }, - "random": { - "description": "Impute missing values with the a random value from the feature set.", - "type": "string" - } - } - } - } - } - } - } - -- **Example** - -.. code-block:: JSON - - { - "default": - { - "feature": { - "continuous": { - "missingCutoffps": 0.25, - "covCutoff": 0.1, - "missingCutoffpf": 0.1, - "imputation": "mean" - } - } - } -.. note:: - Note that you can add as many methods as you want, for other feature types (categorical, ordinal, etc.) and for other cleaning methods (e.g. ``PCA``). +Defines how missing values and low-variance features are handled. + +.. code-block:: yaml + + data_cleaning: + default: + continuous: + missingCutoffps: 0.25 # Max % missing features per sample + covCutoff: 0.1 # Min coefficient of variation + missingCutoffpf: 0.1 # Max % missing samples per feature + imputation: "mean" # Options: "mean", "median", "random" Data Normalization Parameters ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -Data normalization aims to remove batch effects from the data. This set of parameters is used to define the data normalization process Parameters, it is organized as follows: - -.. code-block:: JSON - - { - "standardCombat": { - "define parameters here" - } - } - -.. jsonschema:: - - { - "$schema": "http://json-schema.org/draft-04/schema#", - "title": "Chosen method parameters", - "description": "Normalization method name.", - "type": "string", - "options": { - "standardCombat": { - "description": "Standard Combat normalization method.", - "type": "string" - } - } - } -.. note:: - For now only the ``standardCombat`` method is available and it does not require any parameters. +Aims to remove batch effects (e.g., multicenter differences). + +.. code-block:: yaml + + normalization: + standardCombat: "RUN" + standardization: + perClass: 0 + perInstitution: 0 + minmax: + min: 0 + max: 1 Feature Set Reduction Parameters ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -Feature set reduction consists of reducing the number of features in the data by removing correlated features, selecting important features, etc. This set of parameters is used to define the feature set reduction process Parameters, it is organized as follows: - -.. code-block:: JSON - - { - "selected method": { - "define parameters here" - } - } - -.. jsonschema:: - - { - "$schema": "http://json-schema.org/draft-04/schema#", - "title": "method name", - "description": "Feature set reduction method name.", - "type": "string", - "options": { - "FDA": { - "description": "False discovery avoidance method. `Read the paper. `__", - "type": "string" - }, - "FDAbalanced": { - "description": "Balanced version of the False discovery avoidance method, where the selected number of features is the same for each table.", - "type": "string" - } - } - } - -Now let's specify the parameters for the selected feature set reduction method; for instance, in the case of the ``FDA`` method: - -.. jsonschema:: - - { - "$schema": "http://json-schema.org/draft-04/schema#", - "title": "FDA method", - "description": "Feature set reduction parameters.", - "type": "object", - "properties": { - "FDA": { - "description": "FDA method's parameters.", - "type": "object", - "properties": { - "nSplits": { - "description": "Number of splits to use for the FDA algorithm.", - "type": "int" - }, - "corrType": { - "description": "Type of correlation to use for the FDA algorithm. Default is ``Spearman``.", - "type": "string", - "options": { - "Spearman": { - "description": "Spearman correlation.", - "type": "string" - }, - "Pearson": { - "description": "Pearson correlation.", - "type": "string" - } - } - }, - "threshStableStart": { - "description": "Stability threshold to cut-off the unstable features at the beginning of the FDA algorithm.", - "type": "float" - }, - "threshInterCorr": { - "description": "Threshold to cut-off the inter-correlated features.", - "type": "float" - }, - "minNfeatStable": { - "description": "Minimum number of stable features to keep before inter-correlation step.", - "type": "int" - }, - "minNfeatInterCorr": { - "description": "Minimum number of inter-correlated features to keep.", - "type": "int" - }, - "minNfeat": { - "description": "Minimum number of features to keep at the end of the FDA algorithm.", - "type": "int" - }, - "seed": { - "description": "Seed for the random number generator.", - "type": "int" - } - } - } - } - } - -- **Example** - -.. code-block:: JSON - - { - "FDA": { - "nSplits": 100, - "corrType": "Spearman", - "threshStableStart": 0.5, - "threshInterCorr": 0.7, - "minNfeatStable": 100, - "minNfeatInterCorr": 60, - "minNfeat": 5, - "seed": 54288 - } - } -.. note:: - Only ``FDA`` and ``FDAbalanced`` methods are available for now and they share the same parameters. +Parameters for reducing high-dimensional feature sets (like Radiomics) to a stable subset. + +.. code-block:: yaml + + feature_reduction: + FDA: + nSplits: 100 + corrType: "Spearman" # Options: "Spearman", "Pearson" + threshStableStart: 0.5 + threshInterCorr: 0.7 + minNfeatStable: 100 + minNfeat: 10 + seed: *global_seed Machine Learning Parameters ^^^^^^^^^^^^^^^^^^^^^^^^^^^ -This set of parameters is used to define the machine learning process, algorithm, and parameters, it is organized as follows: - -.. code-block:: JSON - - { - "selected algorithm": { - "define parameters here" - } - } - -Now let's specify the parameters for the selected machine learning algorithm; for instance, in the case of the ``XGBoost`` algorithm: - -.. jsonschema:: - - { - "$schema": "http://json-schema.org/draft-04/schema#", - "title": "ML Algorithm", - "description": "Machine learning algorithm name.", - "type": "object", - "properties": { - "XGBoost": { - "description": "`XGBoost `__ algorithm.", - "type": "object", - "properties": { - "varImportanceThreshold": { - "description": "Variable importance threshold. Default is ``0.3``. Variables with importance below this threshold will be removed.", - "type": "float" - }, - "optimalThreshold": { - "description": "If ``null``, the optimal threshold will be computed. Default is ``0.5``.", - "type": "float" - }, - "optimizationMetric": { - "description": "Model's optimization metric. Default is ``AUC``. Only used if ``method`` is ``pycaret``.", - "type": "string" - }, - "method": { - "description": "Method to use for the XGBoost algorithm. Default is ``pycaret``.", - "type": "string", - "options": { - "pycaret": { - "description": "Automated using `PyCaret `__.", - "type": "string" - }, - "random_search": { - "description": "Random search using a pre-defined grid of parameters.", - "type": "string" - }, - "grid_search": { - "description": "Grid search using a pre-defined grid of parameters.", - "type": "string" - } - } - }, - "nameSave" : { - "description": "Name of the file to save the model.", - "type": "string" - }, - "seed" : { - "description": "Seed for the random number generator.", - "type": "int" - } - } - } - } - } - -- **Example** - -.. code-block:: JSON - - { - "XGBoost": { - "varImportanceThreshold": 0.3, - "optimalThreshold": null, - "optimizationMetric": "AUC", - "method": "pycaret", - "nameSave": "XGBoost03AUC", - "seed": 54288 - } - } +Defines the algorithm and hyperparameter optimization settings. -.. note:: - Only the ``XGBoost`` algorithm is available for now. +.. code-block:: yaml -Variables Definition -^^^^^^^^^^^^^^^^^^^^ + modeling: + method: "firth" # Options: "firth", "rf", "xgboost" + optimization_metric: "MCC" + cv_folds: 5 + var_importance_threshold: 0.05 -This set of parameters is used to define the variables to use for the machine learning process, it is organized as follows: - -.. code-block:: JSON - - { - "selected variable": { - "define parameters here" - }, - "combinations": [ - "Insert combinations of variables here" - ] - } - -.. jsonschema:: - - { - "$schema": "http://json-schema.org/draft-04/schema#", - "title": "Variables", - "description": "Variables to use for the machine learning process.", - "type": "object", - "properties": { - "combinations": { - "description": "List of variables combinations to use for the study.", - "type": "List[str]" - } - } - } - -For the selected variable, you can specify the following parameters: - -.. jsonschema:: - - { - "$schema": "http://json-schema.org/draft-04/schema#", - "title": "selected variable", - "description": "Variable name to use for the machine learning process.", - "type": "object", - "properties": { - "nameType": { - "description": "Type of variable to use. Must contain ``Radiomics`` for radiomics features.", - "type": "string" - }, - "path": { - "description": "Path to the variable file. Use ``\"setToFolderNameinWorkspace\"`` to set the features folder to ``FolderName`` in the workspace.", - "type": "string" - }, - "scans": { - "description": "List of scans to use for the variable. For example is ``T1C``.", - "type": "List[str]" - }, - "rois": { - "description": "List of ROIs to include in the study (will be used to identify the features fie). For example is ``GTV``.", - "type": "List[str]" - }, - "imSpaces": { - "description": "Radiomics level, the features file must end with this level. For example is ``morph``.", - "type": "List[str]" - }, - "var_datacleaning": { - "description": "Data cleaning method to use for the variable. Default is ``default``.", - "type": "string" - }, - "var_normalization": { - "description": "Data normalization method to use for the variable. Default is ``combat``.", - "type": "string" - }, - "var_fSetReduction": { - "description": "Feature set reduction method to use for the variable. Default is ``FDA``.", - "type": "string" - } - } - } - -- **Example** - -.. code-block:: JSON - - { - "var1": { - "nameType": "RadiomicsMorph", - "path": "setToMyFeaturesInWorkspace", - "scans": ["T1CE"], - "rois": ["GTV"], - "imSpaces": ["morph"], - "var_datacleaning": "default", - "var_normalization": "combat", - "var_fSetReduction": "FDA" - }, - "combinations": [ - "var1" - ] - } +.. note:: + For rare-event studies (e.g., only 5 positive cases), it is highly recommended to use the ``firth`` method or a ``Random Forest`` with ``class_weight='balanced'``. + +Full Configuration Example +^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Below is a complete example of a ``config.yaml`` file incorporating all the sections discussed above. You can copy this into your project workspace to get started. + +.. code-block:: yaml + + # ============================================================================== + # Master Configuration for Machine Learning Pipeline + # ============================================================================== + # Note: All lines below are indented by 3 spaces to satisfy the readthedocs directive + + study_metadata: + var_study: "var1" + combinations: ["var1"] + seed: &global_seed 54288 + + design: + active_method: "CrossValidation" + + Random: + method: "SubSampling" + nSplits: 10 + stratifyInstitutions: 1 + testProportion: 0.33 + seed: *global_seed + + CrossValidation: + method: "StratifiedKFold" + nFolds: 5 + nRepeats: 10 + seed: *global_seed + + Bootstrapping: + method: "Out-of-Bag" + nIterations: 1000 + seed: *global_seed + + variables: + var1: + nameType: "RadiomicsFull" + path: "path/to/features/workspace" + scans: ["CECT"] + rois: ["tumor"] + imSpaces: ["image"] + cleaning_profile: "default" + normalization: "combat" + reduction_method: "FDA" + + data_cleaning: + default: + continuous: + missingCutoffps: 0.25 + covCutoff: 0.1 + missingCutoffpf: 0.1 + imputation: "mean" + + normalization: + standardCombat: "RUN" + standardization: + perClass: 0 + perInstitution: 0 + minmax: + min: 0 + max: 1 + + feature_reduction: + FDA: + nSplits: 100 + corrType: "Spearman" + threshStableStart: 0.5 + threshInterCorr: 0.7 + minNfeatStable: 100 + minNfeat: 10 + seed: *global_seed + + modeling: + method: "firth" + optimization_metric: "MCC" + cv_folds: 5 + var_importance_threshold: 0.05 \ No newline at end of file diff --git a/docs/processing.rst b/docs/processing.rst index 6f56feb..b91fff2 100644 --- a/docs/processing.rst +++ b/docs/processing.rst @@ -2,10 +2,10 @@ Processing =========================== -compute\_suv\_map +PET SUV map computation -------------------------------------------- -.. automodule:: MEDiml.processing.compute_suv_map +.. automodule:: MEDiml.processing.PETSUVConverter :members: :undoc-members: :show-inheritance: From d15cce6823a4d70d70c45aa1c1f17cb7dab80e2f Mon Sep 17 00:00:00 2001 From: MahdiAll99 Date: Thu, 19 Mar 2026 08:45:18 -0400 Subject: [PATCH 15/15] forgot this important file - init commit --- MEDiml/processing/SUVHeaderProxy.py | 46 +++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) create mode 100644 MEDiml/processing/SUVHeaderProxy.py diff --git a/MEDiml/processing/SUVHeaderProxy.py b/MEDiml/processing/SUVHeaderProxy.py new file mode 100644 index 0000000..63d9559 --- /dev/null +++ b/MEDiml/processing/SUVHeaderProxy.py @@ -0,0 +1,46 @@ +class SUVHeaderProxy: + """ + Mimics pydicom's dcm[tag].value and dcm[seq][0][tag].value syntax. + """ + class Element: + def __init__(self, value): + self.value = value + + def __getitem__(self, idx): + # Allows dcm[seq][0] + if isinstance(self.value, list): + return SUVHeaderProxy(self.value[idx]) + return self.value + + def __len__(self): + return len(self.value) if isinstance(self.value, list) else 1 + + def __init__(self, data_dict): + self.data = data_dict if data_dict is not None else {} + + def __getitem__(self, tag): + # Convert (0xGGGG, 0xEEEE) tuple to 0xGGGGEEEE integer + if isinstance(tag, tuple): + tag = (tag[0] << 16) | tag[1] + + val = self.data.get(tag) + if val is None: + raise KeyError(f"Tag {hex(tag) if isinstance(tag, int) else tag} not found") + return self.Element(val) + + def __getattribute__(self, name): + # Fallback for .value calls on the proxy itself (during recursion) + if name == "value": + return self.data + return super().__getattribute__(name) + + def __contains__(self, tag): + if isinstance(tag, tuple): + tag = (tag[0] << 16) | tag[1] + return tag in self.data + + def contains(self, tag): + return self.__contains__(tag) + + def get(self, tag, default=None): + return self.data.get(tag, default)