Fix test_priors_to_measurements (#336)

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.

---------

Co-authored-by: Dilan Pathirana <[email protected]>

Author: dweindl

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,31 @@
 )
 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)
+    # All priors in this test case are defined on parameter scale, hence
+    # the dummy measurements will take the scaled nominal values.
+    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 +64,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 +79,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 +88,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 +137,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,
@@ -117,9 +150,7 @@ def test_priors_to_measurements(problem_id):
         prior_type, prior_pars, par_scale, par_bounds = prior
         if prior_type == petab.v1.PARAMETER_SCALE_NORMAL:
             prior_contrib += norm.logpdf(
-                petab_problem_priors.x_nominal_free_scaled[
-                    petab_problem_priors.x_free_ids.index(parameter_id)
-                ],
+                x_scaled_dict[parameter_id],
                 loc=prior_pars[0],
                 scale=prior_pars[1],
             )
@@ -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)