Merge branch 'main' into empty_columns_fine

Author: dweindl

README.md

@@ -27,9 +27,9 @@ It will require Python>=3.10 to run. (We are following the
 
 Development versions of the PEtab library can be installed using
 
-    pip3 install https://github.com/PEtab-dev/libpetab-python/archive/develop.zip
+    pip3 install https://github.com/PEtab-dev/libpetab-python/archive/main.zip
 
-(replace `develop` by the branch or commit you would like to install).
+(replace `main` by the branch or commit you would like to install).
 
 When setting up a new parameter estimation problem, the most useful tools will
 be:

doc/conf.py

@@ -114,7 +114,7 @@
     "display_github": True,
     "github_user": "petab-dev",
     "github_repo": "libpetab-python",
-    "github_version": "develop",
+    "github_version": "main",
     "conf_py_path": "/doc",
 }
 

petab/petablint.py

@@ -6,6 +6,7 @@
 import logging
 import sys
 
+import pydantic
 from colorama import Fore
 from colorama import init as init_colorama
 from jsonschema.exceptions import ValidationError as SchemaValidationError
@@ -160,10 +161,23 @@ def main():
         try:
             validate(args.yaml_file_name)
         except SchemaValidationError as e:
+            path = ""
+            if e.absolute_path:
+                # construct a path to the error location inside the YAML file
+                path = list(e.absolute_path)
+                path = (
+                    f" at {path[0]}"
+                    + "".join(f"[{str(p)}]" for p in path[1:])
+                    + ": "
+                )
             logger.error(
-                f"Provided YAML file does not adhere to PEtab schema: {e}"
+                "Provided YAML file does not adhere to the PEtab schema"
+                f"{path}: {e.args[0]}"
             )
             sys.exit(1)
+        except ValueError as e:
+            logger.error(e)
+            sys.exit(1)
 
         if petab.is_composite_problem(args.yaml_file_name):
             # TODO: further checking:
@@ -179,12 +193,24 @@ def main():
             case 2:
                 from petab.v2.lint import lint_problem
 
-                validation_issues = lint_problem(args.yaml_file_name)
-                if validation_issues:
-                    validation_issues.log(logger=logger)
+                try:
+                    validation_issues = lint_problem(args.yaml_file_name)
+                    if validation_issues:
+                        # Handle petab.v2.lint.ValidationTask issues
+                        validation_issues.log(logger=logger)
+                        sys.exit(1)
+                    logger.info("PEtab format check completed successfully.")
+                    sys.exit(0)
+                except pydantic.ValidationError as e:
+                    # Handle Pydantic validation errors
+                    for err in e.errors():
+                        loc = ", ".join(str(loc) for loc in err["loc"])
+                        msg = err["msg"]
+                        # TODO: include model info here once available
+                        #  https://github.com/pydantic/pydantic/issues/7224
+                        logger.error(f"Error in field(s) `{loc}`: {msg}")
                     sys.exit(1)
-                logger.info("PEtab format check completed successfully.")
-                sys.exit(0)
+
             case _:
                 logger.error(
                     "The provided PEtab files are of unsupported version "

petab/v1/calculate.py

@@ -1,6 +1,7 @@
 """Functions performing various calculations."""
 
 import numbers
+import operator
 from functools import reduce
 
 import numpy as np
@@ -139,19 +140,17 @@ def calculate_residuals_for_table(
             # apply scaling
             observable = observable_df.loc[row[OBSERVABLE_ID]]
             trafo = observable.get(OBSERVABLE_TRANSFORMATION, LIN)
-            simulation = petab.scale(simulation, trafo)
-            measurement = petab.scale(measurement, trafo)
+            scaled_simulation = petab.scale(simulation, trafo)
+            scaled_measurement = petab.scale(measurement, trafo)
 
         # non-normalized residual is just the difference
-        residual = simulation - measurement
+        residual = scaled_measurement - scaled_simulation
 
-        noise_value = 1
         if normalize:
-            # look up noise standard deviation
-            noise_value = evaluate_noise_formula(
+            # divide by standard deviation
+            residual /= evaluate_noise_formula(
                 row, noise_formulas, parameter_df, simulation
             )
-        residual /= noise_value
 
         # fill in value
         residual_df.loc[irow, RESIDUAL] = residual
@@ -169,13 +168,10 @@ def get_symbolic_noise_formulas(observable_df) -> dict[str, sp.Expr]:
     """
     noise_formulas = {}
     # iterate over observables
-    for row in observable_df.itertuples():
-        observable_id = row.Index
-        if NOISE_FORMULA not in observable_df.columns:
-            noise_formula = None
-        else:
-            noise_formula = sympify_petab(row.noiseFormula)
-        noise_formulas[observable_id] = noise_formula
+    for observable_id, row in observable_df.iterrows():
+        noise_formulas[observable_id] = (
+            sympify_petab(row.noiseFormula) if NOISE_FORMULA in row else None
+        )
     return noise_formulas
 
 
@@ -364,7 +360,7 @@ def calculate_llh_for_table(
             (simulation_df[col] == row[col]) | petab.is_empty(row[col])
             for col in compared_cols
         ]
-        mask = reduce(lambda x, y: x & y, masks)
+        mask = reduce(operator.and_, masks)
 
         simulation = simulation_df.loc[mask][SIMULATION].iloc[0]
 
@@ -375,7 +371,7 @@ def calculate_llh_for_table(
 
         # get noise standard deviation
         noise_value = evaluate_noise_formula(
-            row, noise_formulas, parameter_df, petab.scale(simulation, scale)
+            row, noise_formulas, parameter_df, simulation
         )
 
         # get noise distribution

petab/v1/measurements.py

@@ -289,7 +289,7 @@ def assert_overrides_match_parameter_count(
             )
         }
         if NOISE_FORMULA in observable_df.columns
-        else {obs_id: 0 for obs_id in observable_df.index.values}
+        else dict.fromkeys(observable_df.index.values, 0)
     )
 
     for _, row in measurement_df.iterrows():

petab/v1/sampling.py

@@ -33,7 +33,7 @@ def sample_from_prior(
         bounds=tuple(bounds),
         transformation=scaling,
     )
-    return prior.sample(shape=(n_starts,))
+    return prior.sample(shape=(n_starts,), x_scaled=True)
 
 
 def sample_parameter_startpoints(

petab/v1/simulate.py

@@ -241,20 +241,39 @@ def sample_noise(
         simulated_value,
     )
 
-    # default noise distribution is petab.C.NORMAL
-    noise_distribution = petab_problem.observable_df.loc[
+    observable_row = petab_problem.observable_df.loc[
         measurement_row[petab.C.OBSERVABLE_ID]
-    ].get(petab.C.NOISE_DISTRIBUTION, petab.C.NORMAL)
+    ]
+    # default noise distribution is petab.C.NORMAL
+    noise_distribution = observable_row.get(
+        petab.C.NOISE_DISTRIBUTION, petab.C.NORMAL
+    )
     # an empty noise distribution column in an observables table can result in
     # `noise_distribution == float('nan')`
     if pd.isna(noise_distribution):
         noise_distribution = petab.C.NORMAL
 
+    observable_transformation = observable_row.get(
+        petab.C.OBSERVABLE_TRANSFORMATION, petab.C.LIN
+    )
+    transform = lambda x: x  # noqa: E731
+    # observableTransformation=log -> the log of the simulated value is
+    #  distributed according to `noise_distribution`
+    if observable_transformation == petab.C.LOG:
+        simulated_value = np.log(simulated_value)
+        transform = np.exp
+    elif observable_transformation == petab.C.LOG10:
+        simulated_value = np.log10(simulated_value)
+        transform = lambda x: np.power(10, x)  # noqa: E731
+
     # below is e.g.: `np.random.normal(loc=simulation, scale=noise_value)`
     simulated_value_with_noise = getattr(rng, noise_distribution)(
         loc=simulated_value, scale=noise_value * noise_scaling_factor
     )
 
+    # apply observable transformation, ensure `float` type
+    simulated_value_with_noise = float(transform(simulated_value_with_noise))
+
     if zero_bounded and np.sign(simulated_value) != np.sign(
         simulated_value_with_noise
     ):

petab/v1/visualize/plot_residuals.py

@@ -134,6 +134,7 @@ def plot_goodness_of_fit(
     petab_problem: Problem,
     simulations_df: str | Path | pd.DataFrame,
     size: tuple = (10, 7),
+    color=None,
     ax: plt.Axes | None = None,
 ) -> matplotlib.axes.Axes:
     """
@@ -148,6 +149,9 @@ def plot_goodness_of_fit(
         output data file.
     size:
         Figure size.
+    color:
+        The marker colors, matches the `c` parameter of
+        `matplotlib.pyplot.scatter`.
     ax:
         Axis object.
 
@@ -171,8 +175,8 @@ def plot_goodness_of_fit(
         parameter_dfs=petab_problem.parameter_df,
     )[0]
     slope, intercept, r_value, p_value, std_err = stats.linregress(
-        petab_problem.measurement_df["measurement"],
         simulations_df["simulation"],
+        petab_problem.measurement_df["measurement"],
     )  # x, y
 
     if ax is None:
@@ -182,6 +186,7 @@ def plot_goodness_of_fit(
     ax.scatter(
         petab_problem.measurement_df["measurement"],
         simulations_df["simulation"],
+        c=color,
     )
 
     ax.axis("square")

petab/v1/visualize/plotting.py

@@ -861,9 +861,15 @@ def parse_from_id_list(
         """
         if ids_per_plot is None:
             # this is the default case. If no grouping is specified,
-            # all observables are plotted. One observable per plot.
-            unique_obs_list = self._data_df[OBSERVABLE_ID].unique()
-            ids_per_plot = [[obs_id] for obs_id in unique_obs_list]
+            # each group_by category will be plotted on a separate plot
+            unique_ids_list = self._data_df[
+                {
+                    "dataset": DATASET_ID,
+                    "observable": OBSERVABLE_ID,
+                    "simulation": SIMULATION_CONDITION_ID,
+                }[group_by]
+            ].unique()
+            ids_per_plot = [[id_] for id_ in unique_ids_list]
 
         if group_by == "dataset" and DATASET_ID not in self._data_df:
             raise ValueError(

petab/v1/yaml.py

@@ -77,16 +77,14 @@ def validate_yaml_syntax(
         #  but let's still use the latest PEtab schema for full validation
         version = yaml_config.get(FORMAT_VERSION, None)
         version = (
-            parse_version(version)[:2]
-            if version
-            else list(SCHEMAS.values())[-1]
+            parse_version(version)[:2] if version else list(SCHEMAS.keys())[-1]
         )
 
         try:
             schema = SCHEMAS[version]
         except KeyError as e:
             raise ValueError(
-                "Unknown PEtab version given in problem "
+                "No or unknown PEtab version given in problem "
                 f"specification: {version}"
             ) from e
     schema = load_yaml(schema)