v2: Startpoint sampling

petab/v1/distributions.py

@@ -1,4 +1,13 @@
-"""Probability distributions used by PEtab."""
+"""Probability distributions used by PEtab.
+
+This module provides a set of univariate probability distributions
+that can be used for sampling and evaluating the probability density
+function (PDF) and cumulative distribution function (CDF).
+Most of these distributions also support log transformations and truncation.
+
+Not all distributions that can be represented by these classes are valid
+as PEtab parameter prior or noise distributions.
+"""
 
 from __future__ import annotations
 
@@ -19,9 +28,14 @@
 
 __all__ = [
     "Distribution",
+    "Cauchy",
+    "ChiSquare",
+    "Exponential",
+    "Gamma",
+    "Laplace",
     "Normal",
+    "Rayleigh",
     "Uniform",
-    "Laplace",
 ]
 
 

petab/v2/core.py

@@ -25,6 +25,7 @@
 )
 from typing_extensions import Self
 
+from ..v1.distributions import *
 from ..v1.lint import is_valid_identifier
 from ..v1.math import petab_math_str, sympify_petab
 from . import C, get_observable_df
@@ -150,6 +151,26 @@ class PriorDistribution(str, Enum):
     f"{set(C.PRIOR_DISTRIBUTIONS)} vs { {e.value for e in PriorDistribution} }"
 )
 
+_prior_to_cls = {
+    PriorDistribution.CAUCHY: Cauchy,
+    PriorDistribution.CHI_SQUARED: ChiSquare,
+    PriorDistribution.EXPONENTIAL: Exponential,
+    PriorDistribution.GAMMA: Gamma,
+    PriorDistribution.LAPLACE: Laplace,
+    PriorDistribution.LOG10_NORMAL: Normal,
+    PriorDistribution.LOG_LAPLACE: Laplace,
+    PriorDistribution.LOG_NORMAL: Normal,
+    PriorDistribution.LOG_UNIFORM: Uniform,
+    PriorDistribution.NORMAL: Normal,
+    PriorDistribution.RAYLEIGH: Rayleigh,
+    PriorDistribution.UNIFORM: Uniform,
+}
+
+assert not (_mismatch := set(PriorDistribution) ^ set(_prior_to_cls)), (
+    "PriorDistribution enum does not match _prior_to_cls. "
+    f"Mismatches: {_mismatch}"
+)
+
 
 class Observable(BaseModel):
     """Observable definition."""
@@ -929,6 +950,38 @@ def _validate(self) -> Self:
 
         return self
 
+    @property
+    def prior_dist(self) -> Distribution:
+        """Get the pior distribution of the parameter."""
+        if self.estimate is False:
+            raise ValueError(f"Parameter `{self.id}' is not estimated.")
+
+        if self.prior_distribution is None:
+            return Uniform(self.lb, self.ub)
+
+        if not (cls := _prior_to_cls.get(self.prior_distribution)):
+            raise ValueError(
+                f"Prior distribution `{self.prior_distribution}' not "
+                "supported."
+            )
+
+        if str(self.prior_distribution).startswith("log-"):
+            log = True
+        elif str(self.prior_distribution).startswith("log10-"):
+            log = 10
+        else:
+            log = False
+
+        if cls == Exponential:
+            # `Exponential.__init__` does not accept the `log` parameter
+            if log is not False:
+                raise ValueError(
+                    "Exponential distribution does not support log "
+                    "transformation."
+                )
+            return cls(*self.prior_parameters, trunc=[self.lb, self.ub])
+        return cls(*self.prior_parameters, log=log, trunc=[self.lb, self.ub])
+
 
 class ParameterTable(BaseModel):
     """PEtab parameter table."""

petab/v2/problem.py

@@ -12,6 +12,7 @@
 from pathlib import Path
 from typing import TYPE_CHECKING, Any
 
+import numpy as np
 import pandas as pd
 import sympy as sp
 from pydantic import AnyUrl, BaseModel, Field
@@ -22,10 +23,10 @@
     observables,
     parameter_mapping,
     parameters,
-    sampling,
     yaml,
 )
 from ..v1.core import concat_tables, get_visualization_df
+from ..v1.distributions import Distribution
 from ..v1.models.model import Model, model_factory
 from ..v1.yaml import get_path_prefix
 from ..v2.C import *  # noqa: F403
@@ -726,24 +727,29 @@ def get_optimization_to_simulation_parameter_mapping(self, **kwargs):
             )
         )
 
-    def sample_parameter_startpoints(self, n_starts: int = 100, **kwargs):
-        """Create 2D array with starting points for optimization
+    def get_priors(self) -> dict[str, Distribution]:
+        """Get prior distributions.
 
-        See :py:func:`petab.sample_parameter_startpoints`.
+        :returns: The prior distributions for the estimated parameters.
         """
-        return sampling.sample_parameter_startpoints(
-            self.parameter_df, n_starts=n_starts, **kwargs
-        )
+        return {
+            p.id: p.prior_dist
+            for p in self.parameter_table.parameters
+            if p.estimate
+        }
+
+    def sample_parameter_startpoints(self, n_starts: int = 100, **kwargs):
+        """Create 2D array with starting points for optimization"""
+        priors = self.get_priors()
+        return np.vstack([p.sample(n_starts) for p in priors.values()]).T
 
     def sample_parameter_startpoints_dict(
         self, n_starts: int = 100
     ) -> list[dict[str, float]]:
         """Create dictionaries with starting points for optimization
 
-        See also :py:func:`petab.sample_parameter_startpoints`.
-
-        Returns:
-            A list of dictionaries with parameter IDs mapping to samples
+        :returns:
+            A list of dictionaries with parameter IDs mapping to sampled
             parameter values.
         """
         return [

tests/v2/test_problem.py

@@ -23,6 +23,7 @@
     TARGET_VALUE,
     UPPER_BOUND,
 )
+from petab.v2.core import *
 
 
 def test_load_remote():
@@ -170,3 +171,22 @@ def test_modify_problem():
         }
     ).set_index([PETAB_ENTITY_ID])
     assert_frame_equal(problem.mapping_df, exp_mapping_df, check_dtype=False)
+
+
+def test_sample_startpoint_shape():
+    """Test startpoint sampling."""
+    problem = Problem()
+    problem += Parameter(id="p1", estimate=True, lb=1, ub=2)
+    problem += Parameter(
+        id="p2",
+        estimate=True,
+        lb=2,
+        ub=3,
+        prior_distribution="normal",
+        prior_parameters=[2.5, 0.5],
+    )
+    problem += Parameter(id="p3", estimate=False, nominal_value=1)
+
+    n_starts = 10
+    sp = problem.sample_parameter_startpoints(n_starts=n_starts)
+    assert sp.shape == (n_starts, 2)