From d4d4a8069a959f1e0fabef566f81b24cbb90a40f Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 2 Oct 2025 08:43:03 +0200 Subject: [PATCH] Combine observables, sigmas and noise distributions parameters for model import Combine observables, sigmas and noise distributions parameters for model import. Closes #2923. --- CHANGELOG.md | 12 + .../ExampleExperimentalConditions.ipynb | 18 +- .../ExampleEquilibrationLogic.ipynb | 169 +++--- .../GettingStartedExtended.ipynb | 57 +-- pytest.ini | 4 +- python/benchmark/benchmark_pysb.py | 6 +- python/sdist/amici/__init__.py | 3 +- python/sdist/amici/de_model.py | 41 +- python/sdist/amici/import_utils.py | 134 ++++- python/sdist/amici/petab/import_helpers.py | 85 ++-- python/sdist/amici/petab/pysb_import.py | 40 +- python/sdist/amici/petab/sbml_import.py | 50 +- python/sdist/amici/pysb_import.py | 215 +++----- python/sdist/amici/sbml_import.py | 479 ++++++------------ python/sdist/amici/testing/models.py | 114 ++--- python/tests/conftest.py | 8 +- python/tests/test_bngl.py | 4 +- .../test_compare_conservation_laws_sbml.py | 23 +- python/tests/test_events.py | 13 +- python/tests/test_jax.py | 9 +- python/tests/test_pysb.py | 14 +- python/tests/test_sbml_import.py | 101 ++-- .../test_sbml_import_special_functions.py | 27 +- 23 files changed, 773 insertions(+), 853 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8165a8187d..c5c1efeffa 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,6 +14,18 @@ BREAKING CHANGES * `ReturnDataView.posteq_status` and `ReturnDataView.preeq_status` now return `list[SteadyStateStatus]` instead of an `ndarray[int]` of shape `(1, 3)`. +* The way the observation model is passed to the different import functions + (`sbml2amici`, `pysb2amici`, ...) has changed: + The information previously split across `observables`, `sigmas`, + `noise_distributions`, `event_observables`, `event_sigmas`, and + `event_noise_distributions` dicts is now passed as a single + `observation_model: list[MeasurementChannel]` argument. +* `assignmentRules2observables` has been renamed to + `assignment_rules_to_observables` and now returns `list[MeasurementChannel]` + to be passed to the import functions. +* `assignmentRules2observables` has been renamed to + `assignment_rules_to_observables` and now returns `list[MeasurementChannel]` + to be passed to the import functions. * `AmiVector::getLength` and `AmiVectorArray::getLength` have been renamed to `size` to be more consistent with STL containers. * The following deprecated functionality has been removed: diff --git a/doc/examples/example_presimulation/ExampleExperimentalConditions.ipynb b/doc/examples/example_presimulation/ExampleExperimentalConditions.ipynb index ee06015152..572c07085d 100644 --- a/doc/examples/example_presimulation/ExampleExperimentalConditions.ipynb +++ b/doc/examples/example_presimulation/ExampleExperimentalConditions.ipynb @@ -189,9 +189,7 @@ { "cell_type": "markdown", "metadata": {}, - "source": [ - "The SBML model specifies a single observable named `pPROT` which describes the fraction of phosphorylated Protein. We load this observable using [amici.assignmentRules2observables](https://amici.readthedocs.io/en/latest/generated/amici.sbml_import.html#amici.sbml_import.assignmentRules2observables)." - ] + "source": "The SBML model specifies a single observable named `pPROT` which describes the fraction of phosphorylated Protein. We load this observable using [amici.assignment_rules_to_observables](https://amici.readthedocs.io/en/latest/generated/amici.sbml_import.html#amici.sbml_import.assignment_rules_to_observables)." }, { "cell_type": "code", @@ -211,7 +209,7 @@ ], "source": [ "# Retrieve model output names and formulae from AssignmentRules and remove the respective rules\n", - "observables = amici.assignmentRules2observables(\n", + "observables = amici.assignment_rules_to_observables(\n", " sbml_importer.sbml, # the libsbml model object\n", " filter_function=lambda variable: variable.getName() == \"pPROT\",\n", ")\n", @@ -236,7 +234,7 @@ " model_name,\n", " model_output_dir,\n", " verbose=False,\n", - " observables=observables,\n", + " observation_model=observables,\n", " constant_parameters=fixedParameters,\n", ")\n", "# load the generated module\n", @@ -292,7 +290,7 @@ "# Run simulation using default model parameters and solver options\n", "model.setTimepoints(np.linspace(0, 60, 60))\n", "rdata = amici.runAmiciSimulation(model, solver)\n", - "amici.plotting.plotObservableTrajectories(rdata)" + "amici.plotting.plot_observable_trajectories(rdata)" ] }, { @@ -324,7 +322,7 @@ "edata = amici.ExpData(rdata, 0.1, 0.0)\n", "edata.fixedParameters = [0, 2]\n", "rdata = amici.runAmiciSimulation(model, solver, edata)\n", - "amici.plotting.plotObservableTrajectories(rdata)" + "amici.plotting.plot_observable_trajectories(rdata)" ] }, { @@ -356,7 +354,7 @@ "source": [ "edata.fixedParametersPreequilibration = [3, 0]\n", "rdata = amici.runAmiciSimulation(model, solver, edata)\n", - "amici.plotting.plotObservableTrajectories(rdata)" + "amici.plotting.plot_observable_trajectories(rdata)" ] }, { @@ -404,7 +402,7 @@ ], "source": [ "rdata = amici.runAmiciSimulation(model, solver, edata)\n", - "amici.plotting.plotObservableTrajectories(rdata)" + "amici.plotting.plot_observable_trajectories(rdata)" ] }, { @@ -449,7 +447,7 @@ "print(edata.fixedParametersPresimulation)\n", "print(edata.fixedParameters)\n", "rdata = amici.runAmiciSimulation(model, solver, edata)\n", - "amici.plotting.plotObservableTrajectories(rdata)" + "amici.plotting.plot_observable_trajectories(rdata)" ] } ], diff --git a/doc/examples/example_steady_states/ExampleEquilibrationLogic.ipynb b/doc/examples/example_steady_states/ExampleEquilibrationLogic.ipynb index 8d104ccf8b..a7840e1127 100644 --- a/doc/examples/example_steady_states/ExampleEquilibrationLogic.ipynb +++ b/doc/examples/example_steady_states/ExampleEquilibrationLogic.ipynb @@ -80,7 +80,9 @@ }, { "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ "# Import necessary libraries and define the model\n", "import amici\n", @@ -107,6 +109,7 @@ " AMICI_SUCCESS,\n", " AMICI_ERROR,\n", " ExpData,\n", + " MeasurementChannel as MC,\n", ")\n", "\n", "# We encode the model in Antimony format (https://tellurium.readthedocs.io/en/latest/antimony.html),\n", @@ -156,23 +159,23 @@ "temp_dir = Path(tempfile.mkdtemp())\n", "model_output_dir = temp_dir / model_name\n", "model_reduced_output_dir = temp_dir / model_reduced_name" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ "# Import the model\n", "sbml_importer = amici.SbmlImporter(antimony2sbml(ant_model), from_file=False)\n", "\n", "# specify observables and constant parameters\n", "constant_parameters = [\"synthesis_substrate\", \"init_enzyme\"]\n", - "observables = {\n", - " \"observable_product\": {\"name\": \"\", \"formula\": \"product\"},\n", - " \"observable_substrate\": {\"name\": \"\", \"formula\": \"substrate\"},\n", - "}\n", + "observation_model = [\n", + " MC(id_=\"observable_product\", name=\"\", formula=\"product\"),\n", + " MC(id_=\"observable_substrate\", name=\"\", formula=\"substrate\"),\n", + "]\n", "sigmas = {\"observable_product\": 1.0, \"observable_substrate\": 1.0}\n", "\n", "# import the model without handled the conserved quantity\n", @@ -181,9 +184,8 @@ "sbml_importer.sbml2amici(\n", " model_name,\n", " model_output_dir,\n", - " observables=observables,\n", + " observation_model=observation_model,\n", " constant_parameters=constant_parameters,\n", - " sigmas=sigmas,\n", " compute_conservation_laws=False,\n", ")\n", "\n", @@ -192,17 +194,16 @@ "sbml_importer.sbml2amici(\n", " model_reduced_name,\n", " model_reduced_output_dir,\n", - " observables=observables,\n", + " observation_model=observation_model,\n", " constant_parameters=constant_parameters,\n", - " sigmas=sigmas,\n", ")" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ "# import the models and run some test simulations\n", "model_reduced_module = import_model_module(\n", @@ -232,22 +233,20 @@ "plot_state_trajectories(rdata_reduced, model=model_reduced, ax=axes[0, 1])\n", "plot_observable_trajectories(rdata_reduced, model=model_reduced, ax=axes[1, 1])\n", "fig.tight_layout()" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ "# the enzyme state was removed from the ODEs of the reduced model\n", "print(model.getStateIdsSolver())\n", "print(model_reduced.getStateIdsSolver())\n", "print(model.getStateIds())\n", "print(model_reduced.getStateIds())" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -291,9 +290,11 @@ }, { "cell_type": "code", + "execution_count": null, "metadata": { "scrolled": false }, + "outputs": [], "source": [ "# Call post-equilibration by setting an infinity timepoint\n", "model.setTimepoints([np.inf])\n", @@ -322,9 +323,7 @@ " \"posteq_cpu_timeB\",\n", "]:\n", " print(f\"{key:>16}: {rdata[key]}\")" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -350,15 +349,17 @@ ] }, { - "metadata": {}, "cell_type": "code", - "source": "list(SteadyStateStatus)", + "execution_count": null, + "metadata": {}, "outputs": [], - "execution_count": null + "source": [ + "list(SteadyStateStatus)" + ] }, { - "metadata": {}, "cell_type": "markdown", + "metadata": {}, "source": [ "In our case, only the second entry of `posteq_status` contains a positive integer:\n", "The first run of Newton's method failed due to a Jacobian, which could not be factorized, but the second run (simulation) contains the entry `1` (success).\n", @@ -372,7 +373,9 @@ }, { "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ "# steady state value found by pre-equilibration\n", "steady_state = rdata[\"x\"][0]\n", @@ -384,9 +387,7 @@ "\n", "for stst_value in steady_state:\n", " plt.axhline(y=stst_value, color=\"gray\", linestyle=\"--\", linewidth=1)" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -397,7 +398,9 @@ }, { "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ "# reduce maxsteps for integration\n", "model.setTimepoints([np.inf])\n", @@ -409,9 +412,7 @@ "print(\n", " \"Number of steps employed in post-equilibration:\", rdata[\"posteq_numsteps\"]\n", ")" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -423,9 +424,11 @@ }, { "cell_type": "code", + "execution_count": null, "metadata": { "scrolled": false }, + "outputs": [], "source": [ "model_reduced.setTimepoints([np.inf])\n", "model_reduced.setSteadyStateComputationMode(\n", @@ -449,9 +452,7 @@ " \"Number of steps employed in post-equilibration:\",\n", " rdata_reduced[\"posteq_numsteps\"],\n", ")" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -464,7 +465,9 @@ }, { "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ "# create edata, with 3 timepoints and 2 observables:\n", "edata = ExpData(2, 0, 0, np.array([0.0, 0.1, 1.0]))\n", @@ -473,15 +476,15 @@ "# set parameters for pre-equilibration\n", "edata.fixedParametersPreequilibration = np.array([0.0, 2.0])\n", "edata.reinitializeFixedParameterInitialStates = True" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "code", + "execution_count": null, "metadata": { "scrolled": true }, + "outputs": [], "source": [ "# create the solver object and run the simulation\n", "solver_reduced = model_reduced.getSolver()\n", @@ -503,9 +506,7 @@ "\n", "plot_state_trajectories(rdata_reduced, model=model_reduced)\n", "plot_observable_trajectories(rdata_reduced, model=model_reduced)" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -519,7 +520,9 @@ }, { "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ "# Change the last timepoint to an infinity timepoint.\n", "edata.setTimepoints(np.array([0.0, 0.1, float(\"inf\")]))\n", @@ -527,9 +530,7 @@ "# run the simulation\n", "rdata_reduced = runAmiciSimulation(model_reduced, solver_reduced, edata)\n", "assert rdata_reduced.status == AMICI_SUCCESS" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -581,9 +582,11 @@ }, { "cell_type": "code", + "execution_count": null, "metadata": { "scrolled": false }, + "outputs": [], "source": [ "# Call simulation with singular Jacobian and `integrateIfNewtonFails` mode\n", "model.setTimepoints([np.inf])\n", @@ -605,9 +608,7 @@ ")\n", "print(\"Computed state sensitivities:\")\n", "print(rdata[\"sx\"][0, :, :])" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -619,7 +620,9 @@ }, { "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ "# Call simulation with singular Jacobian and newtonOnly mode (will fail)\n", "model.setTimepoints([np.inf])\n", @@ -640,9 +643,7 @@ ")\n", "print(\"Computed state sensitivities:\")\n", "print(rdata[\"sx\"][0, :, :])" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -655,9 +656,11 @@ }, { "cell_type": "code", + "execution_count": null, "metadata": { "scrolled": true }, + "outputs": [], "source": [ "# Try `newtonOnly` option with reduced model\n", "model_reduced.setTimepoints([np.inf])\n", @@ -689,9 +692,7 @@ ")\n", "print(\"Computed state sensitivities:\")\n", "print(rdata_reduced[\"sx\"][0, :, :])" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -765,7 +766,9 @@ }, { "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ "# Call post-equilibration and sensitivities computation using adjoint sensitivity analysis\n", "# by setting an infinity timepoint\n", @@ -798,9 +801,7 @@ " rdata_reduced[\"posteq_numstepsB\"],\n", ")\n", "print(\"Computed gradient:\", rdata_reduced[\"sllh\"])" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -812,7 +813,9 @@ }, { "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ "# Call adjoint post-equilibration with model with singular Jacobian\n", "model.setSteadyStateSensitivityMode(SteadyStateSensitivityMode.newtonOnly)\n", @@ -833,9 +836,7 @@ " rdata[\"posteq_numstepsB\"],\n", ")\n", "print(\"Computed gradient:\", rdata[\"sllh\"])" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -856,9 +857,11 @@ }, { "cell_type": "code", + "execution_count": null, "metadata": { "scrolled": false }, + "outputs": [], "source": [ "# No post-equilibration this time.\n", "# create edata, with 3 timepoints and 2 observables:\n", @@ -881,15 +884,15 @@ "for key, value in rdata.items():\n", " if key[0:6] == \"preeq_\":\n", " print(f\"{key:20s}:\", value)" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "code", + "execution_count": null, "metadata": { "scrolled": false }, + "outputs": [], "source": [ "# Singular Jacobian, use simulation\n", "model.setSteadyStateSensitivityMode(\n", @@ -910,13 +913,13 @@ "for key, value in rdata.items():\n", " if key[0:6] == \"preeq_\":\n", " print(f\"{key:20s}:\", value)" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ "# Non-singular Jacobian, use Newton solver\n", "solver_reduced = model_reduced.getSolver()\n", @@ -928,9 +931,7 @@ "for key, value in rdata_reduced.items():\n", " if key[0:6] == \"preeq_\":\n", " print(f\"{key:20s}:\", value)" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -958,7 +959,9 @@ }, { "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ "# Non-singular Jacobian, use Newton solver and adjoints with initial state sensitivities\n", "solver_reduced = model_reduced.getSolver()\n", @@ -977,13 +980,13 @@ " if key[0:6] == \"preeq_\":\n", " print(f\"{key:20s}:\", value)\n", "print(\"Gradient:\", rdata_reduced[\"sllh\"])" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ "# Non-singular Jacobian, use simulation solver and adjoints with initial state sensitivities\n", "solver_reduced = model_reduced.getSolver()\n", @@ -1004,13 +1007,13 @@ " if key[0:6] == \"preeq_\":\n", " print(f\"{key:20s}:\", value)\n", "print(\"Gradient:\", rdata_reduced[\"sllh\"])" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ "# Non-singular Jacobian, use Newton solver and adjoints with fully adjoint pre-equilibration\n", "solver_reduced = model_reduced.getSolver()\n", @@ -1031,9 +1034,7 @@ " if key[0:6] == \"preeq_\":\n", " print(f\"{key:20s}:\", value)\n", "print(\"Gradient:\", rdata_reduced[\"sllh\"])" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -1042,7 +1043,9 @@ }, { "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ "# Singular Jacobian, use try Newton solver and adjoints with fully adjoint pre-equilibration\n", "solver = model.getSolver()\n", @@ -1063,9 +1066,7 @@ " if key[0:6] == \"preeq_\":\n", " print(f\"{key:20s}:\", value)\n", "print(\"Gradient:\", rdata[\"sllh\"])" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -1087,11 +1088,13 @@ }, { "cell_type": "code", + "execution_count": null, "metadata": { "pycharm": { "name": "#%%\n" } }, + "outputs": [], "source": [ "# Non-singular Jacobian, use simulation\n", "model_reduced.setSteadyStateSensitivityMode(\n", @@ -1128,9 +1131,7 @@ "print(\"\\nCPU time to reach steady state (ms):\")\n", "print(\" lax tolerances: \", rdata_reduced_lax[\"preeq_cpu_time\"])\n", "print(\" strict tolerances:\", rdata_reduced_strict[\"preeq_cpu_time\"])" - ], - "outputs": [], - "execution_count": null + ] } ], "metadata": { diff --git a/doc/examples/getting_started_extended/GettingStartedExtended.ipynb b/doc/examples/getting_started_extended/GettingStartedExtended.ipynb index aba2d94f06..e95bf54801 100644 --- a/doc/examples/getting_started_extended/GettingStartedExtended.ipynb +++ b/doc/examples/getting_started_extended/GettingStartedExtended.ipynb @@ -152,46 +152,40 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### Observables\n", + "### Observation model\n", "\n", - "Specifying observables is beyond the scope of SBML. Here we define them manually.\n", + "Specifying the observation model (i.e., the quantities that are observed, as well as the respective error models) is beyond the scope of SBML. Here we define that manually.\n", + "\n", + "If you are looking for a more scalable way of defining observables, then checkout [PEtab](https://github.com/PEtab-dev/PEtab). Another possibility is using SBML's [`AssignmentRule`s](https://sbml.org/software/libsbml/5.18.0/docs/formatted/python-api/classlibsbml_1_1_assignment_rule.html) to specify model outputs within the SBML file.\n", "\n", - "If you are looking for a more scalable way for defining observables, then checkout [PEtab](https://github.com/PEtab-dev/PEtab). Another possibility is using SBML's [`AssignmentRule`s](https://sbml.org/software/libsbml/5.18.0/docs/formatted/python-api/classlibsbml_1_1_assignment_rule.html) to specify model outputs within the SBML file." - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [], - "source": [ - "# Define observables\n", - "observables = {\n", - " \"observable_x1\": {\"name\": \"\", \"formula\": \"x1\"},\n", - " \"observable_x2\": {\"name\": \"\", \"formula\": \"x2\"},\n", - " \"observable_x3\": {\"name\": \"\", \"formula\": \"x3\"},\n", - " \"observable_x1_scaled\": {\"name\": \"\", \"formula\": \"scaling_x1 * x1\"},\n", - " \"observable_x2_offsetted\": {\"name\": \"\", \"formula\": \"offset_x2 + x2\"},\n", - " \"observable_x1withsigma\": {\"name\": \"\", \"formula\": \"x1\"},\n", - "}" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### $\\sigma$ parameters\n", "\n", - "To specify measurement noise as a parameter, we simply provide a dictionary with (preexisting) parameter names as keys and a list of observable names as values to indicate which sigma parameter is to be used for which observable." + "\n", + "For model import in AMICI, the different types of measurements are represented as `MeasurementChannels`.\n", + "The measurement channel is characterized by an ID, an optional name, the observation function (`MeasurementChannels.formula`), the type of noise distribution (`MeasurementChannels.noise_distribution`, defaults to a normal distribution), and the scale parameter of that distribution (`MeasurementChannels.sigma`).\n", + "The symbols used in the observation function and for the scale parameter must already be defined in the model." ] }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ - "sigmas = {\"observable_x1withsigma\": \"observable_x1withsigma_sigma\"}" + "# Define observation model\n", + "from amici import MeasurementChannel as MC\n", + "\n", + "observation_model = [\n", + " MC(id_=\"observable_x1\", formula=\"x1\"),\n", + " MC(id_=\"observable_x2\", formula=\"x2\"),\n", + " MC(id_=\"observable_x3\", formula=\"x3\"),\n", + " MC(id_=\"observable_x1_scaled\", formula=\"scaling_x1 * x1\"),\n", + " MC(id_=\"observable_x2_offsetted\", formula=\"offset_x2 + x2\"),\n", + " MC(\n", + " id_=\"observable_x1withsigma\",\n", + " formula=\"x1\",\n", + " sigma=\"observable_x1withsigma_sigma\",\n", + " ),\n", + "]" ] }, { @@ -340,9 +334,8 @@ " model_name,\n", " model_output_dir,\n", " verbose=logging.INFO,\n", - " observables=observables,\n", + " observation_model=observation_model,\n", " constant_parameters=constant_parameters,\n", - " sigmas=sigmas,\n", ")" ] }, diff --git a/pytest.ini b/pytest.ini index 5a478bc1f9..a32b3badf3 100644 --- a/pytest.ini +++ b/pytest.ini @@ -1,6 +1,6 @@ [pytest] -addopts = -vv --strict-markers +addopts = -vv --strict-markers --doctest-modules filterwarnings = # warnings are errors @@ -30,4 +30,4 @@ filterwarnings = ignore:jax.* is deprecated:DeprecationWarning -norecursedirs = .git amici_models build doc documentation matlab models ThirdParty amici sdist examples *build* +norecursedirs = .git amici_models build doc documentation models ThirdParty amici sdist examples *build* petab_test_problems pysb_test_models diff --git a/python/benchmark/benchmark_pysb.py b/python/benchmark/benchmark_pysb.py index 079041ed4d..10b3d424ae 100644 --- a/python/benchmark/benchmark_pysb.py +++ b/python/benchmark/benchmark_pysb.py @@ -75,7 +75,11 @@ pysb_model, outdir, compute_conservation_laws=compute_conservation_laws, - observables=list(pysb_model.observables.keys()), + observation_model=list( + map( + amici.MeasurementChannel, pysb_model.observables.keys() + ) + ), ) amici_model_module = amici.import_model_module( diff --git a/python/sdist/amici/__init__.py b/python/sdist/amici/__init__.py index ebab7efd22..f5beae52d5 100644 --- a/python/sdist/amici/__init__.py +++ b/python/sdist/amici/__init__.py @@ -131,8 +131,9 @@ def _imported_from_setup() -> bool: from .de_export import DEExporter # noqa: F401 from .sbml_import import ( # noqa: F401 SbmlImporter, - assignmentRules2observables, + assignment_rules_to_observables, ) + from .import_utils import MeasurementChannel # noqa: F401 try: from .jax import JAXModel diff --git a/python/sdist/amici/de_model.py b/python/sdist/amici/de_model.py index 9666847204..b4e1f043ca 100644 --- a/python/sdist/amici/de_model.py +++ b/python/sdist/amici/de_model.py @@ -1934,9 +1934,46 @@ def _compute_equation(self, name: str) -> None: syms_x = self.sym("x") syms_yz = self.sym(name.removeprefix("sigma")) xs_in_sigma = {} - for sym_yz, eq_yz in zip(syms_yz, self._eqs[name], strict=True): - yz_free_syms = eq_yz.free_symbols + for i, (sym_yz, eq_sigma_yz) in enumerate( + zip( + syms_yz, + self._eqs[name], + strict=True, + ) + ): + yz_free_syms = eq_sigma_yz.free_symbols if tmp := {x for x in syms_x if x in yz_free_syms}: + # Can we replace x symbols by an equivalent observable? + # (currently, only the matching observable is supported) + x_to_eliminate = next(iter(tmp)) + eq_yz = ( + self.eq("y")[i] + if name == "sigmay" + else self.eq("z")[self._z2event[i]][i] + ) + + try: + # solve for the next best x symbol and substitute + # (maybe try another one if we fail?) + replacement = sp.solve( + sp.Eq(sym_yz, eq_yz), x_to_eliminate + ) + except NotImplementedError: + # can't solve + replacement = [] + + if len(replacement) == 1: + self._eqs[name][i] = self._eqs[name][i].subs( + x_to_eliminate, replacement[0] + ) + if not any( + x in self._eqs[name][i].free_symbols + for x in syms_x + ): + # successfully eliminated x symbols + continue + + # Report all x symbols that cannot be replaced xs_in_sigma[sym_yz] = tmp if xs_in_sigma: msg = ", ".join( diff --git a/python/sdist/amici/import_utils.py b/python/sdist/amici/import_utils.py index 6ca7b1a10b..3a8e2ae61b 100644 --- a/python/sdist/amici/import_utils.py +++ b/python/sdist/amici/import_utils.py @@ -67,6 +67,138 @@ class ObservableTransformation(str, enum.Enum): LIN = "lin" +class MeasurementChannel: + """ + A measurement channel (observable) definition. + + Measurement channels define how model state and parameters are mapped to + observables, including any associated noise models. + + Measurement channels can be time-resolved (i.e., defined over the course of + a simulation) or event-resolved (i.e., defined at specific events during + a simulation). Event-resolved observables are associated with a specific + event via the `event_id` attribute. + """ + + def __init__( + self, + id_: str, + name: str | None = None, + formula: str | sp.Expr | None = None, + noise_distribution: str | Callable[[str], str] | None = None, + sigma: str | sp.Expr | float | None = None, + event_id: str | None = None, + ): + """ + Initialize a measurement channel. + + .. note:: + + When providing expressions for (event) observables and their sigmas + as strings (see below), those will be passed to + :func:`sympy.sympify`. The supported grammar is not well-defined. + Note there can be issues with, for example, ``==`` or n-ary (n>2) + comparison operators. + Also note that operator precedence and function names may differ + from SBML L3 formulas or PEtab math expressions. + Passing a sympy expression directly will + be the safer option for more complex expressions. + + .. note:: + + In any math expressions passed to this function, ``time`` will + be interpreted as the time symbol. + + :param id_: + Unique identifier for the measurement channel. + :param name: + Human-readable name for the measurement channel. + :param formula: + Expression defining how the observable is computed from model state + and parameters. + :param noise_distribution: + Noise distribution associated with the observable. + This is usually a string identifier (e.g., 'normal', 'log-normal'; + see + :func:`amici.import_utils.noise_distribution_to_cost_function`). + If ``None``, a normal distribution is assumed. + + Alternatively, a callable can be passed to account for a + custom noise model. The callable must take a single argument + ``str_symbol``, and return a string defining the negative + log-likelihood contribution for a single data point, using + variables ``{str_symbol}``, ``m{str_symbol}``, and + ``sigma{str_symbol}`` for the simulation, measurement, and + scale parameter, respectively. + :param sigma: + Expression representing the scale parameter of the noise + distribution. This can be a numeric value, a sympy expression, + or an expression string that will be passed to + :func:`sympy.sympify`. + :param event_id: + Identifier of the associated event for event-resolved observables. + `None` for time-resolved observables. + + Example usage: + >>> # Time-resolved observable + >>> mc1 = MeasurementChannel( + ... id_="obs1", + ... name="Observable 1", + ... formula="k1 * x1 + k2", + ... noise_distribution="log-normal", + ... sigma="sigma1" + ... ) + >>> mc1 # doctest: +NORMALIZE_WHITESPACE + MeasurementChannel(id_='obs1', name='Observable 1', \ +formula='k1 * x1 + k2', noise_distribution='log-normal', \ +sigma='sigma1', event_id=None) + >>> mc1.is_time_resolved + True + >>> mc1.is_event_resolved + False + >>> # Event-resolved observable + >>> mc2 = MeasurementChannel( + ... id_="obs2", + ... name="Observable 2", + ... formula="x3", + ... noise_distribution="log-normal", + ... sigma="sigma1", + ... event_id="event1" + ... ) + >>> mc2 # doctest: +NORMALIZE_WHITESPACE + MeasurementChannel(id_='obs2', name='Observable 2', formula='x3', \ + noise_distribution='log-normal', sigma='sigma1', event_id='event1') + >>> mc2.is_event_resolved + True + >>> mc2.is_time_resolved + False + """ + self.id = id_ + self.name = name + self.formula = formula + self.noise_distribution = noise_distribution or "normal" + self.sigma = sigma + self.event_id = event_id + + def __repr__(self): + return ( + f"{self.__class__.__name__}(id_={self.id!r}, name={self.name!r}, " + f"formula={self.formula!r}, noise_distribution=" + f"{self.noise_distribution!r}, sigma={self.sigma!r}, " + f"event_id={self.event_id!r})" + ) + + @property + def is_time_resolved(self): + """Whether this is a time-resolved observable.""" + return self.event_id is None + + @property + def is_event_resolved(self) -> bool: + """Whether this is an event-resolved observable.""" + return self.event_id is not None + + def noise_distribution_to_observable_transformation( noise_distribution: str | Callable, ) -> ObservableTransformation: @@ -174,7 +306,7 @@ def noise_distribution_to_cost_function( AMICI uses the logarithm :math:`\\log(\\pi(m|y,\\sigma)`. - In addition to the above mentioned distributions, it is also possible to + In addition to the above-mentioned distributions, it is also possible to pass a function taking a symbol string and returning a log-distribution string with variables '{str_symbol}', 'm{str_symbol}', 'sigma{str_symbol}' for y, m, sigma, respectively. diff --git a/python/sdist/amici/petab/import_helpers.py b/python/sdist/amici/petab/import_helpers.py index d42e99b1e3..1d248c1239 100644 --- a/python/sdist/amici/petab/import_helpers.py +++ b/python/sdist/amici/petab/import_helpers.py @@ -23,38 +23,36 @@ ) from petab.v1.parameters import get_valid_parameters_for_parameter_table from sympy.abc import _clash +from amici.import_utils import MeasurementChannel logger = logging.getLogger(__name__) def get_observation_model( - observable_df: pd.DataFrame, -) -> tuple[dict[str, dict[str, str]], dict[str, str], dict[str, str | float]]: + observable_df: pd.DataFrame | None, +) -> list[MeasurementChannel]: """ Get observables, sigmas, and noise distributions from PEtab observation table in a format suitable for :meth:`amici.sbml_import.SbmlImporter.sbml2amici`. :param observable_df: - PEtab observables table + PEtab observable table - :return: - Tuple of dicts with observables, noise distributions, and sigmas. + :return: The observation model. """ if observable_df is None: - return {}, {}, {} + return [] - observables = {} - sigmas = {} nan_pat = r"^[nN]a[nN]$" - + observation_model = [] for _, observable in observable_df.iterrows(): oid = str(observable.name) # need to sanitize due to https://github.com/PEtab-dev/PEtab/issues/447 name = re.sub(nan_pat, "", str(observable.get(OBSERVABLE_NAME, ""))) formula_obs = re.sub(nan_pat, "", str(observable[OBSERVABLE_FORMULA])) formula_noise = re.sub(nan_pat, "", str(observable[NOISE_FORMULA])) - observables[oid] = {"name": name, "formula": formula_obs} + # FIXME: use PEtab's formula parser here formula_noise_sym = sp.sympify(formula_noise, locals=_clash) formula_obs_sym = sp.sympify(formula_obs, locals=_clash) # PEtab does currently not allow observables in noiseFormula, and @@ -66,48 +64,51 @@ def get_observation_model( formula_noise_sym = formula_noise_sym.subs( formula_obs_sym, sp.Symbol(oid) ) - sigmas[oid] = str(formula_noise_sym) - - noise_distrs = petab_noise_distributions_to_amici(observable_df) + observation_model.append( + MeasurementChannel( + id_=oid, + name=name, + formula=formula_obs, + # FIXME: get rid of the string conversion here + sigma=str(formula_noise_sym), + noise_distribution=petab_noise_distribution_to_amici( + observable + ), + ) + ) - return observables, noise_distrs, sigmas + return observation_model -def petab_noise_distributions_to_amici( - observable_df: pd.DataFrame, -) -> dict[str, str]: +def petab_noise_distribution_to_amici(observable_row) -> str: """ Map from the petab to the amici format of noise distribution identifiers. - :param observable_df: - PEtab observable table + :param observable_row: + A row of the PEtab observable table :return: - dictionary of observable_id => AMICI noise-distributions + The corresponding AMICI noise-distribution """ - amici_distrs = {} - for _, observable in observable_df.iterrows(): - amici_val = "" - - if ( - OBSERVABLE_TRANSFORMATION in observable - and isinstance(observable[OBSERVABLE_TRANSFORMATION], str) - and observable[OBSERVABLE_TRANSFORMATION] - ): - amici_val += observable[OBSERVABLE_TRANSFORMATION] + "-" - - if ( - NOISE_DISTRIBUTION in observable - and isinstance(observable[NOISE_DISTRIBUTION], str) - and observable[NOISE_DISTRIBUTION] - ): - amici_val += observable[NOISE_DISTRIBUTION] - else: - amici_val += "normal" - amici_distrs[observable.name] = amici_val - - return amici_distrs + amici_val = "" + + if ( + OBSERVABLE_TRANSFORMATION in observable_row + and isinstance(observable_row[OBSERVABLE_TRANSFORMATION], str) + and observable_row[OBSERVABLE_TRANSFORMATION] + ): + amici_val += observable_row[OBSERVABLE_TRANSFORMATION] + "-" + + if ( + NOISE_DISTRIBUTION in observable_row + and isinstance(observable_row[NOISE_DISTRIBUTION], str) + and observable_row[NOISE_DISTRIBUTION] + ): + amici_val += observable_row[NOISE_DISTRIBUTION] + else: + amici_val += "normal" + return amici_val def petab_scale_to_amici_scale(scale_str: str) -> int: diff --git a/python/sdist/amici/petab/pysb_import.py b/python/sdist/amici/petab/pysb_import.py index 0c6933c205..a46aae0fdd 100644 --- a/python/sdist/amici/petab/pysb_import.py +++ b/python/sdist/amici/petab/pysb_import.py @@ -13,7 +13,12 @@ import pysb import pysb.bng import sympy as sp -from petab.v1.C import CONDITION_NAME, NOISE_FORMULA, OBSERVABLE_FORMULA +from amici import MeasurementChannel +from petab.v1.C import ( + CONDITION_NAME, + NOISE_FORMULA, + OBSERVABLE_FORMULA, +) from petab.v1.models.pysb_model import PySBModel from ..import_utils import strip_pysb @@ -21,7 +26,7 @@ from . import PREEQ_INDICATOR_ID from .import_helpers import ( get_fixed_parameters, - petab_noise_distributions_to_amici, + petab_noise_distribution_to_amici, ) from .util import get_states_in_condition_table @@ -270,22 +275,19 @@ def import_model_pysb( ) if petab_problem.observable_df is None: - observables = None - sigmas = None - noise_distrs = None + observation_model = [] else: - observables = [ - expr.name - for expr in pysb_model.expressions - if expr.name in petab_problem.observable_df.index + observation_model = [ + MeasurementChannel( + id_=observable.name, + sigma=f"{observable.name}_sigma", + noise_distribution=petab_noise_distribution_to_amici( + observable + ), + ) + for _, observable in petab_problem.observable_df.iterrows() ] - sigmas = {obs_id: f"{obs_id}_sigma" for obs_id in observables} - - noise_distrs = petab_noise_distributions_to_amici( - petab_problem.observable_df - ) - if jax: from amici.pysb_import import pysb2jax @@ -294,9 +296,7 @@ def import_model_pysb( output_dir=model_output_dir, model_name=model_name, verbose=True, - observables=observables, - sigmas=sigmas, - noise_distributions=noise_distrs, + observation_model=observation_model, pysb_model_has_obs_and_noise=True, **kwargs, ) @@ -309,10 +309,8 @@ def import_model_pysb( output_dir=model_output_dir, model_name=model_name, verbose=True, - observables=observables, - sigmas=sigmas, constant_parameters=constant_parameters, - noise_distributions=noise_distrs, + observation_model=observation_model, pysb_model_has_obs_and_noise=True, **kwargs, ) diff --git a/python/sdist/amici/petab/sbml_import.py b/python/sdist/amici/petab/sbml_import.py index 4bcd105657..eb8872e05c 100644 --- a/python/sdist/amici/petab/sbml_import.py +++ b/python/sdist/amici/petab/sbml_import.py @@ -12,6 +12,8 @@ import petab.v1 as petab import sympy as sp from _collections import OrderedDict + +from amici import MeasurementChannel from amici.logging import log_execution_time, set_log_level from petab.v1.models import MODEL_TYPE_SBML from sympy.abc import _clash @@ -152,8 +154,7 @@ def _workaround_initial_states( def _workaround_observable_parameters( - observables: dict[str, dict[str, str]], - sigmas: dict[str, str | float], + observation_model: list[MeasurementChannel], sbml_model: libsbml.Model, output_parameter_defaults: dict[str, float] | None, jax: bool = False, @@ -166,8 +167,10 @@ def _workaround_observable_parameters( the actual SBML import. """ formulas = chain( - (val["formula"] for val in observables.values()), sigmas.values() + (obs.formula for obs in observation_model), + (obs.sigma for obs in observation_model), ) + id_to_observable = {obs.id: obs for obs in observation_model} output_parameters = OrderedDict() for formula in formulas: # we want reproducible parameter ordering upon repeated import @@ -180,25 +183,26 @@ def _workaround_observable_parameters( if jax and (m := re.match(r"(noiseParameter\d+)_(\w+)", sym)): # group1 is the noise parameter, group2 is the observable, don't add to sbml but replace with generic # noise parameter - sigmas[m.group(2)] = str( - sp.sympify(sigmas[m.group(2)], locals=_clash).subs( - free_sym, sp.Symbol(m.group(1)) - ) + # FIXME: get rid of those str(sympify(...)) here and below + id_to_observable[m.group(2)].sigma = str( + sp.sympify( + id_to_observable[m.group(2)].sigma, locals=_clash + ).subs(free_sym, sp.Symbol(m.group(1))) ) elif jax and ( m := re.match(r"(observableParameter\d+)_(\w+)", sym) ): # group1 is the noise parameter, group2 is the observable, don't add to sbml but replace with generic # observable parameter - observables[m.group(2)]["formula"] = str( + id_to_observable[m.group(2)].formula = str( sp.sympify( - observables[m.group(2)]["formula"], locals=_clash + id_to_observable[m.group(2)].formula, locals=_clash ).subs(free_sym, sp.Symbol(m.group(1))) ) elif ( sbml_model.getElementBySId(sym) is None and sym != "time" - and sym not in observables + and sym not in id_to_observable ): output_parameters[sym] = None logger.debug( @@ -351,24 +355,12 @@ def import_model_sbml( "the problem and try again." ) - if petab_problem.observable_df is not None: - observables, noise_distrs, sigmas = get_observation_model( - petab_problem.observable_df - ) - else: - observables = noise_distrs = sigmas = None - - logger.info(f"Observables: {len(observables)}") - logger.info(f"Sigmas: {len(sigmas)}") + observation_model = get_observation_model(petab_problem.observable_df) - if len(sigmas) != len(observables): - raise AssertionError( - f"Number of provided observables ({len(observables)}) and sigmas " - f"({len(sigmas)}) do not match." - ) + logger.info(f"Number of observables: {len(observation_model)}") _workaround_observable_parameters( - observables, sigmas, sbml_model, output_parameter_defaults, jax=jax + observation_model, sbml_model, output_parameter_defaults, jax=jax ) if not jax: @@ -399,9 +391,7 @@ def import_model_sbml( sbml_importer.sbml2jax( model_name=model_name, output_dir=model_output_dir, - observables=observables, - sigmas=sigmas, - noise_distributions=noise_distrs, + observation_model=observation_model, verbose=verbose, **kwargs, ) @@ -410,11 +400,9 @@ def import_model_sbml( sbml_importer.sbml2amici( model_name=model_name, output_dir=model_output_dir, - observables=observables, + observation_model=observation_model, constant_parameters=fixed_parameters, - sigmas=sigmas, allow_reinit_fixpar_initcond=allow_reinit_fixpar_initcond, - noise_distributions=noise_distrs, verbose=verbose, **kwargs, ) diff --git a/python/sdist/amici/pysb_import.py b/python/sdist/amici/pysb_import.py index fb54dd745c..3bc2156981 100644 --- a/python/sdist/amici/pysb_import.py +++ b/python/sdist/amici/pysb_import.py @@ -42,6 +42,7 @@ noise_distribution_to_cost_function, noise_distribution_to_observable_transformation, _default_simplify, + MeasurementChannel, ) from .logging import get_logger, log_execution_time, set_log_level @@ -54,9 +55,7 @@ def pysb2jax( model: pysb.Model, output_dir: str | Path | None = None, - observables: list[str] = None, - sigmas: dict[str, str] = None, - noise_distributions: dict[str, str | Callable] | None = None, + observation_model: list[MeasurementChannel] = None, verbose: int | bool = False, compute_conservation_laws: bool = True, simplify: Callable = _default_simplify, @@ -86,20 +85,16 @@ def pysb2jax( :param output_dir: see :meth:`amici.de_export.ODEExporter.set_paths` - :param observables: - list of :class:`pysb.core.Expression` or :class:`pysb.core.Observable` - names in the provided model that should be mapped to observables - - :param sigmas: - dict of :class:`pysb.core.Expression` names that should be mapped to - sigmas - - :param noise_distributions: - dict with names of observable Expressions as keys and a noise type - identifier, or a callable generating a custom noise formula string - (see :py:func:`amici.import_utils.noise_distribution_to_cost_function` - ). If nothing is passed for some observable id, a normal model is - assumed as default. + :param observation_model: + The different measurement channels that make up the observation + model, see also :class:`amici.import_utils.MeasurementChannel`. + The ID is expected to be the name of a :class:`pysb.Expression` or + :class:`pysb.Observable` in the provided model that should be mapped to + an observable. + ``sigma`` is expected to be the name of a :class:`pysb.Expression` to + be mapped to the scale parameter of the noise distribution. + ``MeasurementChannel.formula`` is expected to be + ``None``. Event-observables are not supported. :param verbose: verbosity level for logging, True/False default to :attr:`logging.DEBUG`/:attr:`logging.ERROR` @@ -125,20 +120,12 @@ def pysb2jax( :param pysb_model_has_obs_and_noise: if set to ``True``, the pysb model is expected to have extra observables and noise variables added """ - if observables is None: - observables = [] - - if sigmas is None: - sigmas = {} - model_name = model_name or model.name set_log_level(logger, verbose) ode_model = ode_model_from_pysb_importer( model, - observables=observables, - sigmas=sigmas, - noise_distributions=noise_distributions, + observation_model=observation_model, compute_conservation_laws=compute_conservation_laws, simplify=simplify, cache_simplify=cache_simplify, @@ -161,10 +148,8 @@ def pysb2jax( def pysb2amici( model: pysb.Model, output_dir: str | Path | None = None, - observables: list[str] = None, + observation_model: list[MeasurementChannel] = None, constant_parameters: list[str] = None, - sigmas: dict[str, str] = None, - noise_distributions: dict[str, str | Callable] | None = None, verbose: int | bool = False, assume_pow_positivity: bool = False, compiler: str = None, @@ -198,20 +183,16 @@ def pysb2amici( :param output_dir: see :meth:`amici.de_export.ODEExporter.set_paths` - :param observables: - list of :class:`pysb.core.Expression` or :class:`pysb.core.Observable` - names in the provided model that should be mapped to observables - - :param sigmas: - dict of :class:`pysb.core.Expression` names that should be mapped to - sigmas - - :param noise_distributions: - dict with names of observable Expressions as keys and a noise type - identifier, or a callable generating a custom noise formula string - (see :py:func:`amici.import_utils.noise_distribution_to_cost_function` - ). If nothing is passed for some observable id, a normal model is - assumed as default. + :param observation_model: + The different measurement channels that make up the observation + model, see also :class:`amici.import_utils.MeasurementChannel`. + The ID is expected to be the name of a :class:`pysb.Expression` or + :class:`pysb.Observable` in the provided model that should be mapped to + an observable. + ``sigma`` is expected to be the name of a :class:`pysb.Expression` to + be mapped to the scale parameter of the noise distribution. + ``MeasurementChannel.formula`` is expected to be + ``None``. Event-observables are not supported. :param constant_parameters: list of :class:`pysb.core.Parameter` names that should be mapped as @@ -233,7 +214,7 @@ def pysb2amici( if set to ``True``, conservation laws are automatically computed and applied such that the state-jacobian of the ODE right-hand-side has full rank. This option should be set to ``True`` when using the Newton - algorithm to compute steadystates + algorithm to compute steady states :param compile: If ``True``, build the python module for the generated model. If false, @@ -256,25 +237,21 @@ def pysb2amici( will be used. :param pysb_model_has_obs_and_noise: - if set to ``True``, the pysb model is expected to have extra observables and noise variables added + if set to ``True``, the pysb model is expected to have extra + observables and noise variables added """ - if observables is None: - observables = [] + if observation_model is None: + observation_model = [] if constant_parameters is None: constant_parameters = [] - if sigmas is None: - sigmas = {} - model_name = model_name or model.name set_log_level(logger, verbose) ode_model = ode_model_from_pysb_importer( model, constant_parameters=constant_parameters, - observables=observables, - sigmas=sigmas, - noise_distributions=noise_distributions, + observation_model=observation_model, compute_conservation_laws=compute_conservation_laws, simplify=simplify, cache_simplify=cache_simplify, @@ -304,9 +281,7 @@ def pysb2amici( def ode_model_from_pysb_importer( model: pysb.Model, constant_parameters: list[str] = None, - observables: list[str] = None, - sigmas: dict[str, str] = None, - noise_distributions: dict[str, str | Callable] | None = None, + observation_model: list[MeasurementChannel] = None, compute_conservation_laws: bool = True, simplify: Callable = sp.powsimp, # Do not enable by default without testing. @@ -326,14 +301,7 @@ def ode_model_from_pysb_importer( :param constant_parameters: see :func:`amici.pysb_import.pysb2amici` - :param observables: - see :func:`amici.pysb_import.pysb2amici` - - :param sigmas: - dict with names of observable Expressions as keys and names of sigma - Expressions as value sigma - - :param noise_distributions: + :param observation_model: see :func:`amici.pysb_import.pysb2amici` :param compute_conservation_laws: @@ -354,7 +322,8 @@ def ode_model_from_pysb_importer( if set to ``True``, the generated model will be compatible with JAX export :param pysb_model_has_obs_and_noise: - if set to ``True``, the pysb model is expected to have extra observables and noise variables added + if set to ``True``, the pysb model is expected to have extra + observables and noise variables added :return: New DEModel instance according to pysbModel @@ -369,11 +338,9 @@ def ode_model_from_pysb_importer( if constant_parameters is None: constant_parameters = [] - if observables is None: - observables = [] - - if sigmas is None: - sigmas = {} + if observation_model is None: + observation_model = [] + observation_model = {channel.id: channel for channel in observation_model} pysb.bng.generate_equations(model, verbose=verbose) @@ -384,22 +351,19 @@ def ode_model_from_pysb_importer( _process_pysb_observables( model, ode, - observables, - sigmas, - noise_distributions, + observation_model, pysb_model_has_obs_and_noise, ) _process_pysb_expressions( model, ode, - observables, - sigmas, - noise_distributions, + observation_model, pysb_model_has_obs_and_noise, ) - ode._has_quadratic_nllh = not noise_distributions or all( - noise_distr in ["normal", "lin-normal", "log-normal", "log10-normal"] - for noise_distr in noise_distributions.values() + ode._has_quadratic_nllh = all( + channel.noise_distribution + in ["normal", "lin-normal", "log-normal", "log10-normal"] + for channel in observation_model.values() ) _process_stoichiometric_matrix(model, ode, constant_parameters) @@ -582,9 +546,7 @@ def _process_pysb_parameters( def _process_pysb_expressions( pysb_model: pysb.Model, ode_model: DEModel, - observables: list[str], - sigmas: dict[str, str], - noise_distributions: dict[str, str | Callable] | None = None, + observation_model: dict[str, MeasurementChannel], pysb_model_has_obs_and_noise: bool = False, ) -> None: r""" @@ -595,24 +557,17 @@ def _process_pysb_expressions( :param pysb_model: pysb model - :param observables: - list of names of :class`pysb.Expression`\ s or - :class:`pysb.Observable`\ s that are to be mapped to DEModel - observables - - :param sigmas: - dict with names of observable :class:`pysb.Expression` / - :class:`pysb.Observable` names as keys and names of sigma - :class:`pysb.Expressions` as values - - :param noise_distributions: - see :func:`amici.pysb_import.pysb2amici` + :param observation_model: + A map of observable names to + :class:`amici.import_utils.MeasurementChannel`. + see also :func:`amici.pysb_import.pysb2amici` :param ode_model: DEModel instance :param pysb_model_has_obs_and_noise: - if set to ``True``, the pysb model is expected to have extra observables and noise variables added + if set to ``True``, the pysb model is expected to have extra + observables and noise variables added """ # we no longer expand expressions here. pysb/bng guarantees that # they are ordered according to their dependency and we can @@ -639,9 +594,7 @@ def _process_pysb_expressions( expr.expr, pysb_model, ode_model, - observables, - sigmas, - noise_distributions, + observation_model, pysb_model_has_obs_and_noise, ) @@ -652,9 +605,7 @@ def _add_expression( expr: sp.Basic, pysb_model: pysb.Model, ode_model: DEModel, - observables: list[str], - sigmas: dict[str, str], - noise_distributions: dict[str, str | Callable] | None = None, + observation_model: dict[str, MeasurementChannel], pysb_model_has_obs_and_noise: bool = False, ): """ @@ -673,23 +624,20 @@ def _add_expression( :param pysb_model: see :py:func:`_process_pysb_expressions` - :param observables: + :param observation_model: see :py:func:`_process_pysb_expressions` - :param sigmas: - see :py:func:`_process_pysb_expressions` - - :param noise_distributions: - see :py:func:`amici.pysb_import.pysb2amici` - :param ode_model: see :py:func:`_process_pysb_expressions` :param pysb_model_has_obs_and_noise: - if set to ``True``, the pysb model is expected to have extra observables and noise variables added + if set to ``True``, the pysb model is expected to have extra + observables and noise variables added """ - if not pysb_model_has_obs_and_noise or name not in observables: - if name in list(sigmas.values()): + if not pysb_model_has_obs_and_noise or name not in observation_model: + if any( + name == channel.sigma for channel in observation_model.values() + ): component = SigmaY else: component = Expression @@ -697,12 +645,8 @@ def _add_expression( component(sym, name, _parse_special_functions(expr)) ) - if name in observables: - noise_dist = ( - noise_distributions.get(name, "normal") - if noise_distributions - else "normal" - ) + if name in observation_model: + noise_dist = observation_model[name].noise_distribution y = sp.Symbol(name) trafo = noise_distribution_to_observable_transformation(noise_dist) @@ -715,7 +659,7 @@ def _add_expression( ) ode_model.add_component(obs) - sigma = _get_sigma(pysb_model, name, sigmas) + sigma = _get_sigma(pysb_model, name, observation_model[name].sigma) if not pysb_model_has_obs_and_noise: ode_model.add_component( SigmaY(sigma, f"sigma_{name}", sp.Float(1.0)) @@ -741,7 +685,7 @@ def _add_expression( def _get_sigma( - pysb_model: pysb.Model, obs_name: str, sigmas: dict[str, str] + pysb_model: pysb.Model, obs_name: str, sigma_name: str | None ) -> sp.Symbol: """ Tries to extract standard deviation symbolic identifier and formula @@ -754,15 +698,14 @@ def _get_sigma( :param obs_name: name of the observable - :param sigmas: - dict of :class:`pysb.core.Expression` names that should be mapped to - sigmas + :param sigma_name: + Name of a :class:`pysb.core.Expression` that should be mapped to a + sigma or ``None``. :return: symbolic variable representing the standard deviation of the observable """ - if obs_name in sigmas: - sigma_name = sigmas[obs_name] + if sigma_name is not None: if sigma_name in pysb_model.expressions.keys(): return pysb_model.expressions[sigma_name] raise ValueError( @@ -776,9 +719,7 @@ def _get_sigma( def _process_pysb_observables( pysb_model: pysb.Model, ode_model: DEModel, - observables: list[str], - sigmas: dict[str, str], - noise_distributions: dict[str, str | Callable] | None = None, + observation_model: dict[str, MeasurementChannel], pysb_model_has_obs_and_noise: bool = False, ) -> None: """ @@ -791,19 +732,13 @@ def _process_pysb_observables( :param ode_model: DEModel instance - :param observables: - list of names of pysb.Expressions or pysb.Observables that are to be - mapped to DEModel observables - - :param sigmas: - dict with names of observable pysb.Expressions/pysb.Observables - names as keys and names of sigma pysb.Expressions as values - - :param noise_distributions: - see :func:`amici.pysb_import.pysb2amici` + :param observation_model: + Mapping from observable name to + :class:`amici.import_utils.MeasurementChannel`. :param pysb_model_has_obs_and_noise: - if set to ``True``, the pysb model is expected to have extra observables and noise variables added + if set to ``True``, the pysb model is expected to have extra + observables and noise variables added """ # only add those pysb observables that occur in the added # Observables as expressions @@ -814,9 +749,7 @@ def _process_pysb_observables( obs.expand_obs(), pysb_model, ode_model, - observables, - sigmas, - noise_distributions, + observation_model, pysb_model_has_obs_and_noise, ) diff --git a/python/sdist/amici/sbml_import.py b/python/sdist/amici/sbml_import.py index 5139f830cc..efce53a601 100644 --- a/python/sdist/amici/sbml_import.py +++ b/python/sdist/amici/sbml_import.py @@ -55,6 +55,7 @@ _xor_to_or, _eq_to_and, _ne_to_or, + MeasurementChannel, ) from .logging import get_logger, log_execution_time, set_log_level from .sbml_utils import SBMLException @@ -278,13 +279,8 @@ def sbml2amici( self, model_name: str, output_dir: str | Path = None, - observables: dict[str, dict[str, str | sp.Expr]] = None, - event_observables: dict[str, dict[str, str | sp.Expr]] = None, constant_parameters: Iterable[str] = None, - sigmas: dict[str, str | float | sp.Expr] = None, - event_sigmas: dict[str, str | float | sp.Expr] = None, - noise_distributions: dict[str, str | Callable] = None, - event_noise_distributions: dict[str, str | Callable] = None, + observation_model: list[MeasurementChannel] = None, verbose: int | bool = logging.ERROR, assume_pow_positivity: bool = False, compiler: str = None, @@ -311,23 +307,6 @@ def sbml2amici( Sensitivity analysis for local parameters is enabled by creating global parameters ``_{reactionId}_{localParameterName}``. - .. note:: - - When providing expressions for (event) observables and their sigmas - as strings (see below), those will be passed to - :func:`sympy.sympify`. The supported grammar is not well defined. - Note there can be issues with, for example, ``==`` or n-ary (n>2) - comparison operators. - Also note that operator precedence and function names may differ - from SBML L3 formulas or PEtab math expressions. - Passing a sympy expression directly will - be the safer option for more complex expressions. - - .. note:: - - In any math expressions passed to this function, ``time`` will - be interpreted as the time symbol. - :param model_name: Name of the generated model package. Note that in a given Python session, only one model with a given @@ -338,72 +317,17 @@ def sbml2amici( :param output_dir: Directory where the generated model package will be stored. - :param observables: - Observables to be added to the model: - - .. code-block:: - - dict( - observableId: { - 'name': observableName, # optional - 'formula': formulaString or sympy expression, - } - ) - - If the observation function is passed as a string, - it will be passed to :func:`sympy.sympify` (see note above). - - :param event_observables: - Event observables to be added to the model: - - .. code-block:: - - dict( - eventObservableId: { - 'name': eventObservableName, # optional - 'event':eventId, - 'formula': formulaString or sympy expression, - } - ) - - If the formula is passed as a string, it will be passed to - :func:`sympy.sympify` (see note above). - :param constant_parameters: list of SBML Ids identifying constant parameters - :param sigmas: - Expression for the scale parameter of the noise distribution for - each observable. This can be a numeric value, a sympy expression, - or an expression string that will be passed to - :func:`sympy.sympify`. - - ``{observableId: sigma}`` - - :param event_sigmas: - Expression for the scale parameter of the noise distribution for - each observable. This can be a numeric value, a sympy expression, - or an expression string that will be passed to - :func:`sympy.sympify`. - - ``{eventObservableId: sigma}`` - - :param noise_distributions: - dictionary(observableId: noise type). - If nothing is passed for some observable id, a normal model is - assumed as default. Either pass a noise type identifier, or a - callable generating a custom noise string. - For noise identifiers, see - :func:`amici.import_utils.noise_distribution_to_cost_function`. - - :param event_noise_distributions: - dictionary(eventObservableId: noise type). - If nothing is passed for some observable id, a normal model is - assumed as default. Either pass a noise type identifier, or a - callable generating a custom noise string. - For noise identifiers, see - :func:`amici.import_utils.noise_distribution_to_cost_function`. - + :param observation_model: + The different measurement channels that make up the observation + model, see :class:`amici.import_utils.MeasurementChannel`. + If ``None``, default observables will be added (for all + state variables of the model and all species, compartments, + and assignment rule targets) and normally distributed + measurement noise is assumed. + If ``[]``, no observables will be added. :param verbose: Verbosity level for logging, ``True``/``False`` defaults to ``logging.Error``/``logging.DEBUG``. @@ -460,13 +384,8 @@ def sbml2amici( set_log_level(logger, verbose) ode_model = self._build_ode_model( - observables=observables, - event_observables=event_observables, constant_parameters=constant_parameters, - sigmas=sigmas, - event_sigmas=event_sigmas, - noise_distributions=noise_distributions, - event_noise_distributions=event_noise_distributions, + observation_model=observation_model, verbose=verbose, compute_conservation_laws=compute_conservation_laws, simplify=simplify, @@ -499,9 +418,7 @@ def sbml2jax( self, model_name: str, output_dir: str | Path = None, - observables: dict[str, dict[str, str]] = None, - sigmas: dict[str, str | float] = None, - noise_distributions: dict[str, str | Callable] = None, + observation_model: list[MeasurementChannel] = None, verbose: int | bool = logging.ERROR, compute_conservation_laws: bool = True, simplify: Callable | None = _default_simplify, @@ -527,21 +444,15 @@ def sbml2jax( :param output_dir: Directory where the generated model package will be stored. - :param observables: - Observables to be added to the model: - ``dictionary( observableId:{'name':observableName - (optional), 'formula':formulaString)})``. - - :param sigmas: - dictionary(observableId: sigma value or (existing) parameter name) - - :param noise_distributions: - dictionary(observableId: noise type). - If nothing is passed for some observable id, a normal model is - assumed as default. Either pass a noise type identifier, or a - callable generating a custom noise string. - For noise identifiers, see - :func:`amici.import_utils.noise_distribution_to_cost_function`. + :param observation_model: + The different measurement channels that make up the observation + model, see :class:`amici.import_utils.MeasurementChannel`. + Only time-resolved observables are supported here. + If ``None``, default observables will be added (for all + state variables of the model and all species, compartments, + and assignment rule targets) and normally distributed + measurement noise is assumed. + If ``[]``, no observables will be added. :param verbose: verbosity level for logging, ``True``/``False`` default to @@ -573,9 +484,7 @@ def sbml2jax( set_log_level(logger, verbose) ode_model = self._build_ode_model( - observables=observables, - sigmas=sigmas, - noise_distributions=noise_distributions, + observation_model=observation_model, verbose=verbose, compute_conservation_laws=compute_conservation_laws, simplify=simplify, @@ -594,13 +503,8 @@ def sbml2jax( def _build_ode_model( self, - observables: dict[str, dict[str, str]] = None, - event_observables: dict[str, dict[str, str]] = None, constant_parameters: Iterable[str] = None, - sigmas: dict[str, str | float] = None, - event_sigmas: dict[str, str | float] = None, - noise_distributions: dict[str, str | Callable] = None, - event_noise_distributions: dict[str, str | Callable] = None, + observation_model: list[MeasurementChannel] = None, verbose: int | bool = logging.ERROR, compute_conservation_laws: bool = True, simplify: Callable | None = _default_simplify, @@ -622,18 +526,6 @@ def _build_ode_model( f"and hard-coded which is not allowed: {invalid}" ) - if sigmas is None: - sigmas = {} - - if event_sigmas is None: - event_sigmas = {} - - if noise_distributions is None: - noise_distributions = {} - - if event_noise_distributions is None: - event_noise_distributions = {} - self._reset_symbols() self._process_sbml( constant_parameters=constant_parameters, @@ -656,10 +548,7 @@ def _build_ode_model( ) compute_conservation_laws = False - self._process_observables(observables, sigmas, noise_distributions) - self._process_event_observables( - event_observables, event_sigmas, event_noise_distributions - ) + self._process_observation_model(observation_model) self._replace_compartments_with_volumes() self._clean_reserved_symbols() @@ -1948,182 +1837,149 @@ def _process_events(self) -> None: "priority": self._sympify(event.getPriority()), } - @log_execution_time("processing SBML observables", logger) - def _process_observables( + @log_execution_time("processing observation model", logger) + def _process_observation_model( self, - observables: dict[str, dict[str, str | sp.Expr]] | None, - sigmas: dict[str, str | float | sp.Expr], - noise_distributions: dict[str, str], + observation_model: list[MeasurementChannel] = None, ) -> None: """ Perform symbolic computations required for observable and objective function evaluation. - :param observables: - dictionary(observableId: {'name':observableName - (optional), 'formula':formulaString)}) - to be added to the model - - :param sigmas: - dictionary(observableId: sigma value or (existing) - parameter name) - - :param noise_distributions: - dictionary(observableId: noise type) - See :py:func:`sbml2amici`. + Add user-provided observables or make all species, and compartments + with assignment rules, observable. """ + if observation_model: + # prepare and validate + for channel in observation_model: + # gather local symbols before parsing observable and sigma + # formulas + self.add_local_symbol( + channel.id, symbol_with_assumptions(channel.id) + ) - _validate_observables( - observables, sigmas, noise_distributions, events=False - ) - - # add user-provided observables or make all species, and compartments - # with assignment rules, observable - if observables: - # gather local symbols before parsing observable and sigma formulas - for obs in observables.keys(): - self.add_local_symbol(obs, symbol_with_assumptions(obs)) + if ( + channel.event_id is not None + and sp.Symbol(channel.event_id) + not in self.symbols[SymbolId.EVENT] + ): + raise ValueError( + "Could not find an event with the event identifier " + f"{channel.event_id} for the event observable with " + f"name {channel.name or channel.id}." + ) + # Add time-resolved observables self.symbols[SymbolId.OBSERVABLE] = { - symbol_with_assumptions(obs): { - "name": definition.get("name", f"y{iobs}"), - "value": self._sympify(definition["formula"]), + symbol_with_assumptions(channel.id): { + "name": channel.name or f"y{iobs}", + "value": self._sympify(channel.formula), "transformation": noise_distribution_to_observable_transformation( - noise_distributions.get(obs, "normal") + channel.noise_distribution ), } - for iobs, (obs, definition) in enumerate(observables.items()) + for iobs, channel in enumerate( + filter(lambda c: c.is_time_resolved, observation_model) + ) } # check for nesting of observables (unsupported) - observable_syms = set(self.symbols[SymbolId.OBSERVABLE].keys()) - for obs in self.symbols[SymbolId.OBSERVABLE].values(): - if any( - sym in observable_syms for sym in obs["value"].free_symbols - ): - raise ValueError( - "Nested observables are not supported, " - f"but observable `{obs['name']} = {obs['value']}` " - "references another observable." - ) - elif observables is None: - self._generate_default_observables() + _check_symbol_nesting( + self.symbols[SymbolId.OBSERVABLE], "observable" + ) - _check_symbol_nesting( - self.symbols[SymbolId.OBSERVABLE], "eventObservable" - ) + # Add event observables + self.symbols[SymbolId.EVENT_OBSERVABLE] = { + symbol_with_assumptions(channel.id): { + "name": channel.name or f"z{iobs}", + "value": self._sympify(channel.formula), + "event": sp.Symbol(channel.event_id), + "transformation": noise_distribution_to_observable_transformation( + channel.noise_distribution + ), + } + for iobs, channel in enumerate( + filter(lambda c: c.is_event_resolved, observation_model) + ) + } + # Check for time errors + wrong_t = sp.Symbol("t") + for eo in self.symbols[SymbolId.EVENT_OBSERVABLE].values(): + if eo["value"].has(wrong_t): + warnings.warn( + f"Event observable {eo['name']} uses `t` in " + "it's formula which is not the time variable. " + "For the time variable, please use `time` " + "instead!", + stacklevel=1, + ) + # check for nesting of observables (unsupported) + _check_symbol_nesting( + self.symbols[SymbolId.EVENT_OBSERVABLE], "eventObservable" + ) - if sigmas: + # FIXME: that name-based approach is a bit fragile noise_pars = list( { name - for sigma in sigmas.values() - for symbol in self._sympify(sigma).free_symbols + for channel in observation_model + if channel.is_time_resolved and channel.sigma is not None + for symbol in self._sympify(channel.sigma).free_symbols if re.match(r"noiseParameter\d+$", (name := str(symbol))) } ) - else: - noise_pars = [] - self.symbols[SymbolId.NOISE_PARAMETER] = { - symbol_with_assumptions(np): {"name": np} for np in noise_pars - } - if observables: + self.symbols[SymbolId.NOISE_PARAMETER] = { + symbol_with_assumptions(noise_par): {"name": noise_par} + for noise_par in noise_pars + } + observable_pars = list( { name - for obs in observables.values() - for symbol in self._sympify(obs["formula"]).free_symbols + for channel in observation_model + if channel.is_time_resolved and channel.sigma is not None + for symbol in self._sympify(channel.formula).free_symbols if re.match( r"observableParameter\d+$", (name := str(symbol)) ) } ) - else: - observable_pars = [] - self.symbols[SymbolId.OBSERVABLE_PARAMETER] = { - symbol_with_assumptions(op): {"name": op} for op in observable_pars - } - - self._process_log_likelihood(sigmas, noise_distributions) - - @log_execution_time("processing SBML event observables", logger) - def _process_event_observables( - self, - event_observables: dict[str, dict[str, str | sp.Expr]], - event_sigmas: dict[str, str | float | sp.Float], - event_noise_distributions: dict[str, str], - ) -> None: - """ - Perform symbolic computations required for observable and objective - function evaluation. - - :param event_observables: - See :py:func:`sbml2amici`. - - :param event_sigmas: - See :py:func:`sbml2amici`. - - :param event_noise_distributions: - See :py:func:`sbml2amici`. - """ - if event_observables is None: - return - - _validate_observables( - event_observables, - event_sigmas, - event_noise_distributions, - events=True, - ) - - # gather local symbols before parsing observable and sigma formulas - for obs, definition in event_observables.items(): - self.add_local_symbol(obs, symbol_with_assumptions(obs)) - # check corresponding event exists - if ( - sp.Symbol(definition["event"]) - not in self.symbols[SymbolId.EVENT] - ): - raise ValueError( - "Could not find an event with the event identifier " - f"{definition['event']} for the event observable with name" - f"{definition['name']}." - ) - - self.symbols[SymbolId.EVENT_OBSERVABLE] = { - symbol_with_assumptions(obs): { - "name": definition.get("name", f"z{iobs}"), - "value": self._sympify(definition["formula"]), - "event": sp.Symbol(definition.get("event")), - "transformation": noise_distribution_to_observable_transformation( - event_noise_distributions.get(obs, "normal") - ), + self.symbols[SymbolId.OBSERVABLE_PARAMETER] = { + symbol_with_assumptions(obs_par): {"name": obs_par} + for obs_par in observable_pars } - for iobs, (obs, definition) in enumerate(event_observables.items()) - } - wrong_t = sp.Symbol("t") - for eo in self.symbols[SymbolId.EVENT_OBSERVABLE].values(): - if eo["value"].has(wrong_t): - warnings.warn( - f"Event observable {eo['name']} uses `t` in " - "it's formula which is not the time variable. " - "For the time variable, please use `time` " - "instead!", - stacklevel=1, - ) + elif observation_model is None: + self._generate_default_observables() - # check for nesting of observables (unsupported) - _check_symbol_nesting( - self.symbols[SymbolId.EVENT_OBSERVABLE], "eventObservable" - ) + # TODO: _process_log_likelihood everything below and + # needs some refactoring + sigmas = { + channel.id: channel.sigma + for channel in observation_model or [] + if channel.is_time_resolved and channel.sigma is not None + } + noise_distributions = { + channel.id: channel.noise_distribution + for channel in observation_model or [] + if channel.is_time_resolved + } + self._process_log_likelihood(sigmas, noise_distributions) + sigmas = { + channel.id: channel.sigma + for channel in observation_model or [] + if channel.is_event_resolved and channel.sigma is not None + } + noise_distributions = { + channel.id: channel.noise_distribution + for channel in observation_model or [] + if channel.is_event_resolved + } + self._process_log_likelihood(sigmas, noise_distributions, events=True) self._process_log_likelihood( - event_sigmas, event_noise_distributions, events=True - ) - self._process_log_likelihood( - event_sigmas, - event_noise_distributions, + sigmas, + noise_distributions, events=True, event_reg=True, ) @@ -3213,12 +3069,16 @@ def _parse_event_trigger(trigger: sp.Expr) -> sp.Expr: ) -def assignmentRules2observables( - sbml_model: libsbml.Model, filter_function: Callable = lambda *_: True -): +def assignment_rules_to_observables( + sbml_model: libsbml.Model, + filter_function: Callable = lambda *_: True, + as_dict=False, +) -> list[MeasurementChannel] | dict[str, MeasurementChannel]: """ Turn assignment rules into observables. + The matching assignment rules and their targets are removed from the model. + :param sbml_model: Model to operate on @@ -3227,30 +3087,35 @@ def assignmentRules2observables( ``True``/``False`` to indicate if the respective rule should be turned into an observable. + :param as_dict: + If ``True``, return a dict instead of a list, using the observable IDs + as keys. + :return: - A dictionary(observableId:{ - 'name': observableName, - 'formula': formulaString - }) + A list or dict of measurement channels. """ observables = {} for rule in sbml_model.getListOfRules(): if rule.getTypeCode() != libsbml.SBML_ASSIGNMENT_RULE: continue - parameter_id = rule.getVariable() - if (p := sbml_model.getParameter(parameter_id)) and filter_function(p): - observables[parameter_id] = { - "name": p.getName() if p.isSetName() else parameter_id, - "formula": sbml_model.getAssignmentRuleByVariable( - parameter_id + target_id = rule.getVariable() + if (p := sbml_model.getParameter(target_id)) and filter_function(p): + observables[target_id] = MeasurementChannel( + id_=target_id, + name=p.getName() if p.isSetName() else target_id, + formula=sbml_model.getAssignmentRuleByVariable( + target_id ).getFormula(), - } + ) + + for target_id in observables: + sbml_model.removeRuleByVariable(target_id) + sbml_model.removeParameter(target_id) - for parameter_id in observables: - sbml_model.removeRuleByVariable(parameter_id) - sbml_model.removeParameter(parameter_id) + if as_dict: + return observables - return observables + return list(observables.values()) def _add_conservation_for_constant_species( @@ -3401,36 +3266,6 @@ def _check_unsupported_functions_sbml( raise SBMLException(str(err)) -def _validate_observables( - observables: dict[str, dict[str, str | sp.Expr]] | None, - sigmas: dict[str, str | float | sp.Expr], - noise_distributions: dict[str, str], - events: bool = False, -) -> None: - if observables is None or not observables: - return - - # Ensure no non-existing observableIds have been specified - # (no problem here, but usually an upstream bug) - unknown_ids = set(sigmas.keys()) - set(observables.keys()) - if unknown_ids: - raise ValueError( - f"Sigma provided for unknown " - f"{'eventO' if events else 'o'}bservableIds: " - f"{unknown_ids}." - ) - - # Ensure no non-existing observableIds have been specified - # (no problem here, but usually an upstream bug) - unknown_ids = set(noise_distributions.keys()) - set(observables.keys()) - if unknown_ids: - raise ValueError( - f"Noise distribution provided for unknown " - f"{'eventO' if events else 'o'}bservableIds: " - f"{unknown_ids}." - ) - - def _check_symbol_nesting( symbols: dict[sp.Symbol, dict[str, sp.Expr]], symbol_type: str ): diff --git a/python/sdist/amici/testing/models.py b/python/sdist/amici/testing/models.py index 542d6bb8ae..a92c9daac1 100644 --- a/python/sdist/amici/testing/models.py +++ b/python/sdist/amici/testing/models.py @@ -4,7 +4,7 @@ from ..antimony_import import antimony2amici, antimony2sbml from pathlib import Path import libsbml -from amici import SbmlImporter, AmiciModel +from amici import SbmlImporter, AmiciModel, MeasurementChannel import sys import tempfile from amici.sbml_utils import amici_time_symbol @@ -113,11 +113,11 @@ def import_model_robertson(outdir: Path = None) -> Model: antimony2amici( robertson_ant, constant_parameters=["k1"], - observables={ - "obs_x1": {"formula": "x1"}, - "obs_x2": {"formula": "1e4 * x2"}, - "obs_x3": {"formula": "x3"}, - }, + observation_model=[ + MeasurementChannel(id_="obs_x1", formula="x1"), + MeasurementChannel(id_="obs_x2", formula="1e4 * x2"), + MeasurementChannel(id_="obs_x3", formula="x3"), + ], model_name=model_name, output_dir=outdir, ) @@ -136,14 +136,14 @@ def import_model_calvetti(outdir: Path = None) -> Model: antimony2amici( calvetti_ant, constant_parameters=["V1ss", "R1ss", "V2ss", "R2ss", "V3ss", "R3ss"], - observables={ - "obs_V1": {"formula": "V1"}, - "obs_V2": {"formula": "V2"}, - "obs_V3": {"formula": "V3"}, - "obs_f0": {"formula": "f0"}, - "obs_f1": {"formula": "f1"}, - "obs_f2": {"formula": "f2"}, - }, + observation_model=[ + MeasurementChannel(id_="obs_V1", formula="V1"), + MeasurementChannel(id_="obs_V2", formula="V2"), + MeasurementChannel(id_="obs_V3", formula="V3"), + MeasurementChannel(id_="obs_f0", formula="f0"), + MeasurementChannel(id_="obs_f1", formula="f1"), + MeasurementChannel(id_="obs_f2", formula="f2"), + ], model_name=model_name, output_dir=outdir, hardcode_symbols=["p1"], @@ -172,9 +172,7 @@ def import_model_dirac(outdir: Path = None) -> Model: antimony2amici( model_dirac_ant, - observables={ - "obs_x2": {"formula": "x2"}, - }, + observation_model=[MeasurementChannel(id_="obs_x2", formula="x2")], model_name=model_name, output_dir=outdir, ) @@ -228,17 +226,6 @@ def import_model_neuron(outdir: Path = None) -> AmiciModel: }, } - observables = { - "y1": { - "name": "v", - "formula": "v", - } - } - - event_observables = { - "z1": {"name": "z1", "event": "event_1", "formula": "time"} - } - sbml_document, sbml_model = create_sbml_model( initial_assignments=initial_assignments, parameters=parameters, @@ -254,9 +241,11 @@ def import_model_neuron(outdir: Path = None) -> AmiciModel: model = create_amici_model( sbml_model, model_name=model_name, - observables=observables, + observation_model=[ + MeasurementChannel(id_="y1", name="v", formula="v"), + MeasurementChannel(id_="z1", event_id="event_1", formula="time"), + ], constant_parameters=constants, - event_observables=event_observables, output_dir=outdir, ) return model @@ -312,18 +301,6 @@ def import_model_events(outdir: Path = None) -> AmiciModel: "event_2": {"trigger": "x1 > x3", "target": [], "assignment": []}, } - observables = { - "y1": { - "name": "y1", - "formula": "p4*(x1+x2+x3)", - } - } - - event_observables = { - "z1": {"name": "z1", "event": "event_1", "formula": "time"}, - "z2": {"name": "z2", "event": "event_2", "formula": "time"}, - } - sbml_document, sbml_model = create_sbml_model( initial_assignments=initial_assignments, parameters=parameters, @@ -336,12 +313,16 @@ def import_model_events(outdir: Path = None) -> AmiciModel: model_name = "model_events_py" constants = ["k1", "k2", "k3", "k4"] + model = create_amici_model( sbml_model, model_name=model_name, - observables=observables, + observation_model=[ + MeasurementChannel(id_="y1", formula="p4*(x1+x2+x3)", name="y1"), + MeasurementChannel(id_="z1", event_id="event_1", formula="time"), + MeasurementChannel(id_="z2", event_id="event_2", formula="time"), + ], constant_parameters=constants, - event_observables=event_observables, output_dir=outdir, ) return model @@ -432,20 +413,21 @@ def import_model_jakstat(outdir: Path = None) -> AmiciModel: SbmlImporter(sbml_model).sbml2amici( constant_parameters=["Omega_cyt", "Omega_nuc"], - observables={ - "obs_pSTAT": { - "formula": "offset_pSTAT + scale_pSTAT/init_STAT*(pSTAT + 2*pSTAT_pSTAT)" - }, - "obs_tSTAT": { - "formula": "offset_tSTAT + scale_tSTAT/init_STAT*(STAT + pSTAT + 2*(pSTAT_pSTAT))" - }, - "obs_spline": {"formula": "u"}, - }, - sigmas={ - "obs_pSTAT": "sigma_pSTAT", - "obs_tSTAT": "sigma_tSTAT", - "obs_spline": "sigma_pEpoR", - }, + observation_model=[ + MeasurementChannel( + id_="obs_pSTAT", + formula="offset_pSTAT + scale_pSTAT/init_STAT*(pSTAT + 2*pSTAT_pSTAT)", + sigma="sigma_pSTAT", + ), + MeasurementChannel( + id_="obs_tSTAT", + formula="offset_tSTAT + scale_tSTAT/init_STAT*(STAT + pSTAT + 2*(pSTAT_pSTAT))", + sigma="sigma_tSTAT", + ), + MeasurementChannel( + id_="obs_spline", formula="u", sigma="sigma_pEpoR" + ), + ], model_name=model_name, output_dir=outdir, ) @@ -480,9 +462,9 @@ def import_model_nested_events(outdir: Path = None) -> AmiciModel: antimony2amici( ant_str, - observables={ - "obs_Virus": {"formula": "Virus"}, - }, + observation_model=[ + MeasurementChannel(id_="obs_Virus", formula="Virus"), + ], model_name=model_name, output_dir=outdir, ) @@ -523,11 +505,11 @@ def import_model_steadystate(outdir: Path = None) -> AmiciModel: antimony2amici( ant_str, constant_parameters=["k1", "k2", "k3", "k4"], - observables={ - "obs_x1": {"formula": "x1"}, - "obs_x2": {"formula": "x2"}, - "obs_x3": {"formula": "x3"}, - }, + observation_model=[ + MeasurementChannel(id_="obs_x1", formula="x1"), + MeasurementChannel(id_="obs_x2", formula="x2"), + MeasurementChannel(id_="obs_x3", formula="x3"), + ], model_name=model_name, output_dir=outdir, ) diff --git a/python/tests/conftest.py b/python/tests/conftest.py index 557d9187a5..59fcadaa9b 100644 --- a/python/tests/conftest.py +++ b/python/tests/conftest.py @@ -7,6 +7,8 @@ import amici import pytest from pathlib import Path + +from amici import MeasurementChannel from amici.testing import TemporaryDirectoryWinSafe as TemporaryDirectory @@ -30,7 +32,7 @@ def sbml_example_presimulation_module(): constant_parameters = ["DRUG_0", "KIN_0"] - observables = amici.assignmentRules2observables( + observables = amici.assignment_rules_to_observables( sbml_importer.sbml, # the libsbml model object filter_function=lambda variable: variable.getName() == "pPROT_obs", ) @@ -41,7 +43,7 @@ def sbml_example_presimulation_module(): model_name=module_name, output_dir=outdir, verbose=False, - observables=observables, + observation_model=observables, constant_parameters=constant_parameters, ) @@ -78,7 +80,7 @@ def pysb_example_presimulation_module(): model, outdir, verbose=True, - observables=["pPROT_obs"], + observation_model=[MeasurementChannel("pPROT_obs")], constant_parameters=constant_parameters, ) diff --git a/python/tests/test_bngl.py b/python/tests/test_bngl.py index 42926e379a..9dbdb68ec1 100644 --- a/python/tests/test_bngl.py +++ b/python/tests/test_bngl.py @@ -69,7 +69,9 @@ def test_compare_to_pysb_simulation(example): kwargs = { "compute_conservation_laws": cl, - "observables": list(pysb_model.observables.keys()), + "observation_model": list( + map(amici.MeasurementChannel, pysb_model.observables.keys()) + ), } with TemporaryDirectoryWinSafe(prefix=pysb_model.name) as outdir: diff --git a/python/tests/test_compare_conservation_laws_sbml.py b/python/tests/test_compare_conservation_laws_sbml.py index 44701f70d5..897aee40c4 100644 --- a/python/tests/test_compare_conservation_laws_sbml.py +++ b/python/tests/test_compare_conservation_laws_sbml.py @@ -4,7 +4,7 @@ import amici import numpy as np import pytest -from amici import SteadyStateStatus +from amici import SteadyStateStatus, ExpData, MeasurementChannel as MC from numpy.testing import assert_allclose from amici.testing import skip_on_valgrind @@ -12,7 +12,7 @@ @pytest.fixture def edata_fixture(): """edata is generated to test pre- and postequilibration""" - edata_pre = amici.ExpData( + edata_pre = ExpData( 2, 0, 0, np.array([0.0, 0.1, 0.2, 0.5, 1.0, 2.0, 5.0, 10.0]) ) edata_pre.setObservedData([1.5] * 16) @@ -21,12 +21,12 @@ def edata_fixture(): edata_pre.reinitializeFixedParameterInitialStates = True # edata for postequilibration - edata_post = amici.ExpData(2, 0, 0, np.array([float("inf")] * 3)) + edata_post = ExpData(2, 0, 0, np.array([float("inf")] * 3)) edata_post.setObservedData([0.75] * 6) edata_post.fixedParameters = np.array([7.5, 30.0]) # edata with both equilibrations - edata_full = amici.ExpData( + edata_full = ExpData( 2, 0, 0, @@ -85,26 +85,23 @@ def models(): # Define constants, observables, sigmas constant_parameters = ["synthesis_substrate", "init_enzyme"] - observables = { - "observable_product": {"name": "", "formula": "product"}, - "observable_substrate": {"name": "", "formula": "substrate"}, - } - sigmas = {"observable_product": 1.0, "observable_substrate": 1.0} + observation_model = [ + MC("observable_product", formula="product", sigma=1.0, name=""), + MC("observable_substrate", formula="substrate", sigma=1.0, name=""), + ] # wrap models with and without conservations laws sbml_importer.sbml2amici( model_name_cl, model_output_dir_cl, - observables=observables, constant_parameters=constant_parameters, - sigmas=sigmas, + observation_model=observation_model, ) sbml_importer.sbml2amici( model_name, model_output_dir, - observables=observables, constant_parameters=constant_parameters, - sigmas=sigmas, + observation_model=observation_model, compute_conservation_laws=False, ) diff --git a/python/tests/test_events.py b/python/tests/test_events.py index 8a16b7cba4..4a9096aaee 100644 --- a/python/tests/test_events.py +++ b/python/tests/test_events.py @@ -5,7 +5,12 @@ import amici import numpy as np import pytest -from amici import import_model_module, SensitivityMethod, SensitivityOrder +from amici import ( + import_model_module, + SensitivityMethod, + SensitivityOrder, + MeasurementChannel as MC, +) from amici.antimony_import import antimony2amici from amici.gradient_check import check_derivatives from amici.testing import skip_on_valgrind @@ -1127,11 +1132,7 @@ def test_posteq_events_are_handled(tempdir): E1: at time > 1: target = target + bolus E2: at some_time >= 2: target = target + bolus """, - observables={ - "obs_target": { - "formula": "target", - } - }, + observation_model=[MC("obs_target", formula="target")], model_name=model_name, output_dir=tempdir, verbose=True, diff --git a/python/tests/test_jax.py b/python/tests/test_jax.py index a7fe0ac41a..3e5d0ee184 100644 --- a/python/tests/test_jax.py +++ b/python/tests/test_jax.py @@ -20,6 +20,7 @@ from amici.jax import JAXProblem, ReturnValue, run_simulations from numpy.testing import assert_allclose from test_petab_objective import lotka_volterra # noqa: F401 +from amici import MeasurementChannel as MC pysb = pytest.importorskip("pysb") @@ -42,8 +43,8 @@ def test_conversion(): pysb.Observable("ab", a(s="b")) with TemporaryDirectoryWinSafe() as outdir: - pysb2amici(model, outdir, verbose=True, observables=["ab"]) - pysb2jax(model, outdir, verbose=True, observables=["ab"]) + pysb2amici(model, outdir, verbose=True, observation_model=[MC("ab")]) + pysb2jax(model, outdir, verbose=True, observation_model=[MC("ab")]) amici_module = amici.import_model_module( module_name=model.name, module_path=outdir @@ -97,13 +98,13 @@ def test_dimerization(): model, outdir, verbose=True, - observables=["a_obs", "b_obs"], + observation_model=[MC("a_obs"), MC("b_obs")], constant_parameters=["ksyn_a", "ksyn_b"], ) pysb2jax( model, outdir, - observables=["a_obs", "b_obs"], + observation_model=[MC("a_obs"), MC("b_obs")], ) amici_module = amici.import_model_module( diff --git a/python/tests/test_pysb.py b/python/tests/test_pysb.py index 1fd7647481..fb2d7d5b9c 100644 --- a/python/tests/test_pysb.py +++ b/python/tests/test_pysb.py @@ -185,7 +185,12 @@ def test_compare_to_pysb_simulation(example): outdir, verbose=logging.INFO, compute_conservation_laws=compute_conservation_laws, - observables=list(pysb_model.observables.keys()), + observation_model=list( + map( + amici.MeasurementChannel, + pysb_model.observables.keys(), + ) + ), ) amici_model_module = amici.import_model_module( @@ -338,7 +343,12 @@ def test_heavyside_and_special_symbols(): ) with TemporaryDirectoryWinSafe(prefix=model.name) as outdir: - pysb2amici(model, outdir, verbose=True, observables=["a"]) + pysb2amici( + model, + outdir, + verbose=True, + observation_model=[amici.MeasurementChannel("a")], + ) model_module = amici.import_model_module( module_name=model.name, module_path=outdir diff --git a/python/tests/test_sbml_import.py b/python/tests/test_sbml_import.py index f72773fc86..2dafd4acbe 100644 --- a/python/tests/test_sbml_import.py +++ b/python/tests/test_sbml_import.py @@ -11,7 +11,10 @@ import pytest from amici.gradient_check import check_derivatives from amici.sbml_import import SbmlImporter, SymbolId -from amici.import_utils import symbol_with_assumptions +from amici.import_utils import ( + symbol_with_assumptions, + MeasurementChannel as MC, +) from numpy.testing import assert_allclose, assert_array_equal from amici import import_model_module from amici.testing import skip_on_valgrind @@ -66,7 +69,7 @@ def test_sbml2amici_no_observables(tempdir): sbml_importer.sbml2amici( model_name=model_name, output_dir=tempdir, - observables=None, + observation_model=None, compute_conservation_laws=False, ) @@ -86,10 +89,10 @@ def test_sbml2amici_nested_observables_fail(tempdir): sbml_importer.sbml2amici( model_name=model_name, output_dir=tempdir, - observables={ - "outer": {"formula": "inner"}, - "inner": {"formula": "S1"}, - }, + observation_model=[ + MC("outer", formula="inner"), + MC("inner", formula="S1"), + ], compute_conservation_laws=False, generate_sensitivity_code=False, compile=False, @@ -104,7 +107,7 @@ def test_nosensi(tempdir): sbml_importer.sbml2amici( model_name=model_name, output_dir=tempdir, - observables=None, + observation_model=None, compute_conservation_laws=False, generate_sensitivity_code=False, ) @@ -142,14 +145,18 @@ def observable_dependent_error_model(): sbml_importer.sbml2amici( model_name=model_name, output_dir=tmpdir, - observables={ - "observable_s1": {"formula": "S1"}, - "observable_s1_scaled": {"formula": "0.5 * S1"}, - }, - sigmas={ - "observable_s1": "0.1 + relative_sigma * observable_s1", - "observable_s1_scaled": "0.02 * observable_s1_scaled", - }, + observation_model=[ + MC( + "observable_s1", + formula="S1", + sigma="0.1 + relative_sigma * S1", + ), + MC( + "observable_s1_scaled", + formula="0.5 * S1", + sigma="0.02 * observable_s1_scaled", + ), + ], ) yield amici.import_model_module( module_name=model_name, module_path=tmpdir @@ -220,22 +227,25 @@ def model_steadystate_module(): sbml_file = MODEL_STEADYSTATE_SCALED_XML sbml_importer = amici.SbmlImporter(sbml_file) - observables = amici.assignmentRules2observables( + observables = amici.assignment_rules_to_observables( sbml_importer.sbml, filter_function=lambda variable: variable.getId().startswith( "observable_" ) and not variable.getId().endswith("_sigma"), + as_dict=True, ) + observables[ + "observable_x1withsigma" + ].sigma = "observable_x1withsigma_sigma" module_name = "test_model_steadystate_scaled" with TemporaryDirectory(prefix=module_name) as outdir: sbml_importer.sbml2amici( model_name=module_name, output_dir=outdir, - observables=observables, constant_parameters=["k0"], - sigmas={"observable_x1withsigma": "observable_x1withsigma_sigma"}, + observation_model=list(observables.values()), ) yield amici.import_model_module( @@ -548,36 +558,27 @@ def model_test_likelihoods(tempdir): sbml_file = MODEL_STEADYSTATE_SCALED_XML sbml_importer = amici.SbmlImporter(sbml_file) - # define observables - observables = { - "o1": {"formula": "x1"}, - "o2": {"formula": "10^x1"}, - "o3": {"formula": "10^x1"}, - "o4": {"formula": "x1"}, - "o5": {"formula": "10^x1"}, - "o6": {"formula": "10^x1"}, - "o7": {"formula": "x1"}, - } - - # define different noise models - noise_distributions = { - "o1": "normal", - "o2": "log-normal", - "o3": "log10-normal", - "o4": "laplace", - "o5": "log-laplace", - "o6": "log10-laplace", - "o7": lambda str_symbol: f"Abs({str_symbol} - m{str_symbol}) " - f"/ sigma{str_symbol}", - } + # define observation model + observation_model = [ + MC("o1", formula="x1", noise_distribution="normal"), + MC("o2", formula="10^x1", noise_distribution="log-normal"), + MC("o3", formula="10^x1", noise_distribution="log10-normal"), + MC("o4", formula="x1", noise_distribution="laplace"), + MC("o5", formula="10^x1", noise_distribution="log-laplace"), + MC("o6", formula="10^x1", noise_distribution="log10-laplace"), + MC( + "o7", + formula="x1", + noise_distribution=lambda str_symbol: f"Abs({str_symbol} - m{str_symbol}) / sigma{str_symbol}", + ), + ] module_name = "model_test_likelihoods" sbml_importer.sbml2amici( model_name=module_name, output_dir=tempdir, - observables=observables, constant_parameters=["k0"], - noise_distributions=noise_distributions, + observation_model=observation_model, ) yield amici.import_model_module( @@ -656,21 +657,16 @@ def test_likelihoods_error(): sbml_file = MODEL_STEADYSTATE_SCALED_XML sbml_importer = amici.SbmlImporter(sbml_file) - # define observables - observables = {"o1": {"formula": "x1"}} - - # define different noise models - noise_distributions = {"o1": "nörmal"} - module_name = "test_likelihoods_error" outdir = "test_likelihoods_error" with pytest.raises(ValueError): sbml_importer.sbml2amici( model_name=module_name, output_dir=outdir, - observables=observables, constant_parameters=["k0"], - noise_distributions=noise_distributions, + observation_model=[ + MC("o1", formula="x1", noise_distribution="nörmal") + ], ) @@ -1178,7 +1174,10 @@ def test_time_dependent_initial_assignment(compute_conservation_laws: bool): si = SbmlImporter(sbml_model, from_file=False) de_model = si._build_ode_model( - observables={"obs_p1": {"formula": "p1"}, "obs_p2": {"formula": "p2"}}, + observation_model=[ + MC("obs_p1", formula="p1"), + MC("obs_p2", formula="p2"), + ], compute_conservation_laws=compute_conservation_laws, ) # "species", because the initial assignment expression is time-dependent diff --git a/python/tests/test_sbml_import_special_functions.py b/python/tests/test_sbml_import_special_functions.py index 0a528c5eb0..70f7bfc187 100644 --- a/python/tests/test_sbml_import_special_functions.py +++ b/python/tests/test_sbml_import_special_functions.py @@ -10,6 +10,7 @@ from amici.antimony_import import antimony2amici from amici.gradient_check import check_derivatives from amici.testing import skip_on_valgrind, TemporaryDirectoryWinSafe +from amici import SbmlImporter, MeasurementChannel as MC from numpy.testing import ( assert_approx_equal, assert_array_almost_equal_nulp, @@ -22,29 +23,21 @@ @pytest.fixture(scope="session") def model_special_likelihoods(): """Test model for special likelihood functions.""" - # load sbml model - sbml_importer = amici.SbmlImporter(MODEL_STEADYSTATE_SCALED_XML) - - # define observables - observables = { - "o1": {"formula": "100*10^x1"}, - "o2": {"formula": "100*10^x1"}, - } - - # define different noise models - noise_distributions = { - "o1": "binomial", - "o2": "negative-binomial", - } - + sbml_importer = SbmlImporter(MODEL_STEADYSTATE_SCALED_XML) module_name = "test_special_likelihoods" with TemporaryDirectoryWinSafe(prefix=module_name) as outdir: sbml_importer.sbml2amici( model_name=module_name, output_dir=outdir, - observables=observables, constant_parameters=["k0"], - noise_distributions=noise_distributions, + observation_model=[ + MC("o1", formula="100*10^x1", noise_distribution="binomial"), + MC( + "o2", + formula="100*10^x1", + noise_distribution="negative-binomial", + ), + ], ) yield amici.import_model_module(