From b15d57e6e7b91768cbe2e21337140a940dc76e5e Mon Sep 17 00:00:00 2001
From: Daniel Weindl <daniel.weindl@uni-bonn.de>
Date: Sun, 8 Dec 2024 16:40:02 +0100
Subject: [PATCH] Fix test_priors_to_measurements

Fixes an issue in `test_priors_to_measurements` which led to evaluating the prior at the wrong parameter values
(using the location parameter of the prior instead of the actually estimated parameters).
The problem was in the test code, not the tested code.
---
 tests/v1/test_priors.py | 67 ++++++++++++++++++++++++++++++-----------
 1 file changed, 49 insertions(+), 18 deletions(-)

diff --git a/tests/v1/test_priors.py b/tests/v1/test_priors.py
index ea47e54f..37fd6590 100644
--- a/tests/v1/test_priors.py
+++ b/tests/v1/test_priors.py
@@ -8,7 +8,15 @@
 from scipy.stats import norm
 
 import petab.v1
-from petab.v1 import get_simulation_conditions
+from petab.v1 import (
+    ESTIMATE,
+    MEASUREMENT,
+    OBJECTIVE_PRIOR_TYPE,
+    OBSERVABLE_ID,
+    SIMULATION,
+    get_simulation_conditions,
+    get_simulation_df,
+)
 from petab.v1.priors import priors_to_measurements
 
 
@@ -17,20 +25,29 @@
 )
 def test_priors_to_measurements(problem_id):
     """Test the conversion of priors to measurements."""
+    # setup
     petab_problem_priors: petab.v1.Problem = (
         benchmark_models_petab.get_problem(problem_id)
     )
     petab_problem_priors.visualization_df = None
     assert petab.v1.lint_problem(petab_problem_priors) is False
-
     if problem_id == "Isensee_JCB2018":
         # required to match the stored simulation results below
         petab.v1.flatten_timepoint_specific_output_overrides(
             petab_problem_priors
         )
         assert petab.v1.lint_problem(petab_problem_priors) is False
+
     original_problem = deepcopy(petab_problem_priors)
+    x_scaled_dict = dict(
+        zip(
+            original_problem.x_free_ids,
+            original_problem.x_nominal_free_scaled,
+            strict=True,
+        )
+    )
 
+    # convert priors to measurements
     petab_problem_measurements = priors_to_measurements(petab_problem_priors)
 
     # check that the original problem is not modified
@@ -45,6 +62,7 @@ def test_priors_to_measurements(problem_id):
                 getattr(original_problem, attr)
             )
         ).empty, diff
+
     # check that measurements and observables were added
     assert petab.v1.lint_problem(petab_problem_measurements) is False
     assert (
@@ -59,6 +77,7 @@ def test_priors_to_measurements(problem_id):
         petab_problem_measurements.measurement_df.shape[0]
         > petab_problem_priors.measurement_df.shape[0]
     )
+
     # ensure we didn't introduce any new conditions
     assert len(
         get_simulation_conditions(petab_problem_measurements.measurement_df)
@@ -67,26 +86,40 @@ def test_priors_to_measurements(problem_id):
     # verify that the objective function value is the same
 
     # load/construct the simulation results
-    simulation_df_priors = petab.v1.get_simulation_df(
+    simulation_df_priors = get_simulation_df(
         Path(
             benchmark_models_petab.MODELS_DIR,
             problem_id,
             f"simulatedData_{problem_id}.tsv",
         )
     )
-    simulation_df_measurements = pd.concat(
-        [
-            petab_problem_measurements.measurement_df.rename(
-                columns={petab.v1.MEASUREMENT: petab.v1.SIMULATION}
-            )[
-                petab_problem_measurements.measurement_df[
-                    petab.v1.C.OBSERVABLE_ID
-                ].str.startswith("prior_")
-            ],
-            simulation_df_priors,
+    # for the prior observables, we need to "simulate" the model with the
+    #  nominal parameter values
+    simulated_prior_observables = (
+        petab_problem_measurements.measurement_df.rename(
+            columns={MEASUREMENT: SIMULATION}
+        )[
+            petab_problem_measurements.measurement_df[
+                OBSERVABLE_ID
+            ].str.startswith("prior_")
         ]
     )
 
+    def apply_parameter_values(row):
+        # apply the parameter values to the observable formula for the prior
+        if row[OBSERVABLE_ID].startswith("prior_"):
+            row[SIMULATION] = x_scaled_dict[
+                row[OBSERVABLE_ID].removeprefix("prior_")
+            ]
+        return row
+
+    simulated_prior_observables = simulated_prior_observables.apply(
+        apply_parameter_values, axis=1
+    )
+    simulation_df_measurements = pd.concat(
+        [simulation_df_priors, simulated_prior_observables]
+    )
+
     llh_priors = petab.v1.calculate_llh_for_table(
         petab_problem_priors.measurement_df,
         simulation_df_priors,
@@ -102,10 +135,8 @@ def test_priors_to_measurements(problem_id):
 
     # get prior objective function contribution
     parameter_ids = petab_problem_priors.parameter_df.index.values[
-        (petab_problem_priors.parameter_df[petab.v1.ESTIMATE] == 1)
-        & petab_problem_priors.parameter_df[
-            petab.v1.OBJECTIVE_PRIOR_TYPE
-        ].notna()
+        (petab_problem_priors.parameter_df[ESTIMATE] == 1)
+        & petab_problem_priors.parameter_df[OBJECTIVE_PRIOR_TYPE].notna()
     ]
     priors = petab.v1.get_priors_from_df(
         petab_problem_priors.parameter_df,
@@ -134,4 +165,4 @@ def test_priors_to_measurements(problem_id):
         llh_priors + prior_contrib, llh_measurements, rtol=1e-3, atol=1e-16
     ), (llh_priors + prior_contrib, llh_measurements)
     # check that the tolerance is not too high
-    assert np.abs(prior_contrib) > 1e-3 * np.abs(llh_priors)
+    assert np.abs(prior_contrib) > 1e-8 * np.abs(llh_priors)