Merge pull request #160 from daisybio/development

v1.2.5

Author: JudithBernett

docs/conf.py

@@ -56,9 +56,9 @@
 # the built documents.
 #
 # The short X.Y version.
-version = "1.2.4"
+version = "1.2.5"
 # The full version, including alpha/beta/rc tags.
-release = "1.2.4"
+release = "1.2.5"
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

drevalpy/evaluation.py

@@ -255,6 +255,16 @@ def evaluate(dataset: DrugResponseDataset, metric: list[str] | str):
                     )
                 )
             else:
-                results[m] = float(AVAILABLE_METRICS[m](y_pred=predictions, y_true=response))
+                # check whether the predictions contain NaNs
+                if np.any(np.isnan(predictions)):
+                    # if there are only NaNs in the predictions, the metric is NaN
+                    if np.all(np.isnan(predictions)):
+                        results[m] = float(np.nan)
+                    else:
+                        # remove the rows with NaNs in the predictions and response
+                        mask = ~np.isnan(predictions)
+                        results[m] = float(AVAILABLE_METRICS[m](y_pred=predictions[mask], y_true=response[mask]))
+                else:
+                    results[m] = float(AVAILABLE_METRICS[m](y_pred=predictions, y_true=response))
 
     return results

drevalpy/visualization/utils.py

@@ -111,7 +111,7 @@ def parse_results(path_to_results: str, dataset: str) -> tuple[pd.DataFrame, pd.
 @pipeline_function
 def evaluate_file(
     pred_file: pathlib.Path, test_mode: str, model_name: str
-) -> tuple[pd.DataFrame, pd.DataFrame, pd.DataFrame, pd.DataFrame, str]:
+) -> tuple[pd.DataFrame, pd.DataFrame | None, pd.DataFrame | None, pd.DataFrame, str]:
     """
     Evaluate the predictions from the final models.
 

poetry.lock

@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "drevalpy"
-version = "1.2.4"
+version = "1.2.5"
 description = "Drug response evaluation of cancer cell line drug response models in a fair setting"
 authors = ["DrEvalPy development team"]
 license = "GPL-3.0"

pyproject.toml

@@ -1,4 +1,4 @@
-aiohappyeyeballs==2.4.8 ; python_version >= "3.11" and python_version < "3.13"
+aiohappyeyeballs==2.5.0 ; python_version >= "3.11" and python_version < "3.13"
 aiohttp==3.11.13 ; python_version >= "3.11" and python_version < "3.13"
 aiosignal==1.3.2 ; python_version >= "3.11" and python_version < "3.13"
 anyio==4.8.0 ; python_version >= "3.11" and python_version < "3.13"
@@ -38,18 +38,18 @@ jaraco-classes==3.4.0 ; python_version >= "3.11" and python_version < "3.13"
 jaraco-context==6.0.1 ; python_version >= "3.11" and python_version < "3.13"
 jaraco-functools==4.1.0 ; python_version >= "3.11" and python_version < "3.13"
 jeepney==0.9.0 ; python_version >= "3.11" and python_version < "3.13" and sys_platform == "linux"
-jinja2==3.1.5 ; python_version >= "3.11" and python_version < "3.13"
+jinja2==3.1.6 ; python_version >= "3.11" and python_version < "3.13"
 joblib==1.4.2 ; python_version >= "3.11" and python_version < "3.13"
 keyring==25.6.0 ; python_version >= "3.11" and python_version < "3.13"
 kiwisolver==1.4.8 ; python_version >= "3.11" and python_version < "3.13"
-lightning-utilities==0.13.1 ; python_version >= "3.11" and python_version < "3.13"
+lightning-utilities==0.14.0 ; python_version >= "3.11" and python_version < "3.13"
 markupsafe==3.0.2 ; python_version >= "3.11" and python_version < "3.13"
 matplotlib==3.10.1 ; python_version >= "3.11" and python_version < "3.13"
 more-itertools==10.6.0 ; python_version >= "3.11" and python_version < "3.13"
 mpmath==1.3.0 ; python_version >= "3.11" and python_version < "3.13"
 msgpack==1.1.0 ; python_version >= "3.11" and python_version < "3.13"
 multidict==6.1.0 ; python_version >= "3.11" and python_version < "3.13"
-narwhals==1.29.0 ; python_version >= "3.11" and python_version < "3.13"
+narwhals==1.29.1 ; python_version >= "3.11" and python_version < "3.13"
 networkx==3.4.2 ; python_version >= "3.11" and python_version < "3.13"
 numpy==1.26.4 ; python_version >= "3.11" and python_version < "3.13"
 nvidia-cublas-cu12==12.4.5.8 ; python_version >= "3.11" and python_version < "3.13" and platform_system == "Linux" and platform_machine == "x86_64"
@@ -117,7 +117,7 @@ trove-classifiers==2025.3.3.18 ; python_version >= "3.11" and python_version < "
 typing-extensions==4.12.2 ; python_version >= "3.11" and python_version < "3.13"
 tzdata==2025.1 ; python_version >= "3.11" and python_version < "3.13"
 urllib3==2.3.0 ; python_version >= "3.11" and python_version < "3.13"
-virtualenv==20.29.2 ; python_version >= "3.11" and python_version < "3.13"
+virtualenv==20.29.3 ; python_version >= "3.11" and python_version < "3.13"
 xarray==2025.1.2 ; python_version >= "3.11" and python_version < "3.13"
 xattr==1.1.4 ; python_version >= "3.11" and python_version < "3.13" and sys_platform == "darwin"
 xyzservices==2025.1.0 ; python_version >= "3.11" and python_version < "3.13"

requirements.txt

@@ -1,15 +1,23 @@
 """Tests for the baselines in the models module."""
 
+import pathlib
 import tempfile
 from typing import cast
 
 import numpy as np
+import pandas as pd
 import pytest
 from sklearn.linear_model import ElasticNet, Ridge
 
 from drevalpy.datasets.dataset import DrugResponseDataset, FeatureDataset
-from drevalpy.evaluation import evaluate, pearson
-from drevalpy.experiment import cross_study_prediction
+from drevalpy.evaluation import evaluate
+from drevalpy.experiment import (
+    consolidate_single_drug_model_predictions,
+    cross_study_prediction,
+    generate_data_saving_path,
+    get_datasets_from_cv_split,
+    train_and_predict,
+)
 from drevalpy.models import (
     MODEL_FACTORY,
     NaiveCellLineMeanPredictor,
@@ -19,6 +27,7 @@
 )
 from drevalpy.models.baselines.sklearn_models import SklearnModel
 from drevalpy.models.drp_model import DRPModel
+from drevalpy.visualization.utils import evaluate_file
 
 
 @pytest.mark.parametrize(
@@ -146,61 +155,67 @@ def test_single_drug_baselines(
     :param test_mode: either LPO or LCO
     :param cross_study_dataset: dataset
     """
-    drug_response = sample_dataset
-    drug_response.split_dataset(
+    sample_dataset.split_dataset(
         n_cv_splits=5,
         mode=test_mode,
     )
-    assert drug_response.cv_splits is not None
-    split = drug_response.cv_splits[0]
-    train_dataset = split["train"]
-    val_dataset = split["validation"]
-
+    assert sample_dataset.cv_splits is not None
+    split = sample_dataset.cv_splits[0]
     model = MODEL_FACTORY[model_name]()
-    cell_line_input = model.load_cell_line_features(data_path="../data", dataset_name="TOYv1")
-    cell_lines_to_keep = cell_line_input.identifiers
 
-    len_train_before = len(train_dataset)
-    len_pred_before = len(val_dataset)
-    train_dataset.reduce_to(cell_line_ids=cell_lines_to_keep, drug_ids=None)
-    val_dataset.reduce_to(cell_line_ids=cell_lines_to_keep, drug_ids=None)
-    print(f"Reduced training dataset from {len_train_before} to {len(train_dataset)}")
-    print(f"Reduced val dataset from {len_pred_before} to {len(val_dataset)}")
-
-    all_unique_drugs = np.unique(train_dataset.drug_ids)
+    # test what happens if a drug is only in the original dataset, not in the cross-study dataset
+    exclusive_drugs = list(set(sample_dataset.drug_ids).difference(set(cross_study_dataset.drug_ids)))
+    all_unique_drugs = list(set(sample_dataset.drug_ids).intersection(set(cross_study_dataset.drug_ids)))
+    all_unique_drugs.sort()
+    exclusive_drugs.sort()
+    all_unique_drugs_arr = np.array(all_unique_drugs)
+    exclusive_drugs_arr = np.array(exclusive_drugs)
     # randomly sample a drug to speed up testing
     np.random.seed(123)
-    np.random.shuffle(all_unique_drugs)
-    random_drug = all_unique_drugs[:1]
-
-    all_predictions = np.zeros_like(val_dataset.drug_ids, dtype=float)
+    np.random.shuffle(all_unique_drugs_arr)
+    np.random.shuffle(exclusive_drugs_arr)
+    random_drugs = all_unique_drugs_arr[:1]
+    random_drugs = np.concatenate([random_drugs, exclusive_drugs_arr[:1]])
+    # test what happens if the training and validation dataset is empty for a drug but the test set is not
+    drug_to_remove = all_unique_drugs_arr[2]
+    random_drugs = np.concatenate([random_drugs, [drug_to_remove]])
 
     hpam_combi = model.get_hyperparameter_set()[0]
+    result_path = tempfile.TemporaryDirectory()
     if model_name == "SingleDrugRandomForest":
         hpam_combi["n_estimators"] = 2  # reduce test time
         hpam_combi["max_depth"] = 2  # reduce test time
-
-    model.build_model(hpam_combi)
-    output_mask = train_dataset.drug_ids == random_drug
-    drug_train = train_dataset.copy()
-    drug_train.mask(output_mask)
-    model.train(output=drug_train, cell_line_input=cell_line_input, drug_input=None)
-
-    val_mask = val_dataset.drug_ids == random_drug
-    all_predictions[val_mask] = model.predict(
-        drug_ids=random_drug,
-        cell_line_ids=val_dataset.cell_line_ids[val_mask],
-        cell_line_input=cell_line_input,
-    )
-    # check whether predictions are constant
-    if np.all(all_predictions[val_mask] == all_predictions[val_mask][0]):
-        print("Predictions are constant")
-    else:
-        pcc_drug = pearson(val_dataset.response[val_mask], all_predictions[val_mask])
-        print(f"{test_mode}: Performance of {model_name} for drug {random_drug}: PCC = {pcc_drug}")
-        assert pcc_drug >= -1.0
-    with tempfile.TemporaryDirectory() as temp_dir:
-        print(f"Running cross-study prediction for {model_name}")
+    for random_drug in random_drugs:
+        predictions_path = generate_data_saving_path(
+            model_name=model_name,
+            drug_id=str(random_drug),
+            result_path=result_path.name,
+            suffix="predictions",
+        )
+        prediction_file = pathlib.Path(predictions_path, "predictions_split_0.csv")
+        (
+            train_dataset,
+            validation_dataset,
+            early_stopping_dataset,
+            test_dataset,
+        ) = get_datasets_from_cv_split(split, MODEL_FACTORY[model_name], model_name, random_drug)
+        train_dataset.add_rows(validation_dataset)
+        if random_drug == drug_to_remove:
+            reduce_to_drugs = np.array(list(set(train_dataset.drug_ids) - {random_drug}))
+            train_dataset.reduce_to(cell_line_ids=None, drug_ids=reduce_to_drugs)
+        train_dataset.shuffle(random_state=42)
+        test_dataset = train_and_predict(
+            model=model,
+            hpams=hpam_combi,
+            path_data="../data",
+            train_dataset=train_dataset,
+            prediction_dataset=test_dataset,
+            early_stopping_dataset=None,
+            response_transformation=None,
+            model_checkpoint_dir="TEMPORARY",
+        )
+        cross_study_dataset.remove_nan_responses()
+        parent_dir = str(pathlib.Path(predictions_path).parent)
         cross_study_prediction(
             dataset=cross_study_dataset,
             model=model,
@@ -209,10 +224,38 @@ def test_single_drug_baselines(
             path_data="../data",
             early_stopping_dataset=None,
             response_transformation=None,
-            path_out=temp_dir,
+            path_out=parent_dir,
             split_index=0,
-            single_drug_id=str(random_drug[0]),
+            single_drug_id=str(random_drug),
         )
+        test_dataset.to_csv(prediction_file)
+    consolidate_single_drug_model_predictions(
+        models=[MODEL_FACTORY[model_name]],
+        n_cv_splits=1,
+        results_path=result_path.name,
+        cross_study_datasets=[cross_study_dataset.dataset_name],
+        randomization_mode=None,
+        n_trials_robustness=0,
+        out_path=result_path.name,
+    )
+    # get cross-study predictions and assert that each drug-cell line combination only occurs once
+    cross_study_predictions = pd.read_csv(
+        pathlib.Path(result_path.name, model_name, "cross_study", "cross_study_TOYv2_split_0.csv")
+    )
+    assert len(cross_study_predictions) == len(cross_study_predictions.drop_duplicates(["drug_id", "cell_line_id"]))
+    predictions_file = pathlib.Path(result_path.name, model_name, "predictions", "predictions_split_0.csv")
+    cross_study_file = pathlib.Path(result_path.name, model_name, "cross_study", "cross_study_TOYv2_split_0.csv")
+    for file in [predictions_file, cross_study_file]:
+        (
+            overall_eval,
+            eval_results_per_drug,
+            eval_results_per_cl,
+            t_vs_p,
+            model_name,
+        ) = evaluate_file(pred_file=file, test_mode=test_mode, model_name=model_name)
+        assert len(overall_eval) == 1
+        print(f"Performance of {model_name}: PCC = {overall_eval['Pearson'][0]}")
+        assert overall_eval["Pearson"][0] >= -1.0
 
 
 def _call_naive_predictor(