From 032fb6c174723080ff33d7665e325bd0f6bc3c8e Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 13 Mar 2024 17:59:11 +0100 Subject: [PATCH 01/17] Add status code AMICI_CONSTR_FAIL (#2379) --- include/amici/defines.h | 1 + src/amici.cpp | 1 + 2 files changed, 2 insertions(+) diff --git a/include/amici/defines.h b/include/amici/defines.h index bd97d9a0b2..6ccba89923 100644 --- a/include/amici/defines.h +++ b/include/amici/defines.h @@ -71,6 +71,7 @@ constexpr int AMICI_CONV_FAILURE= -4; constexpr int AMICI_LSETUP_FAIL= -6; constexpr int AMICI_RHSFUNC_FAIL= -8; constexpr int AMICI_FIRST_RHSFUNC_ERR= -9; +constexpr int AMICI_CONSTR_FAIL= -15; constexpr int AMICI_ILL_INPUT= -22; constexpr int AMICI_ERROR= -99; constexpr int AMICI_NO_STEADY_STATE= -81; diff --git a/src/amici.cpp b/src/amici.cpp index 2f60a8e12a..c74c88e3a5 100644 --- a/src/amici.cpp +++ b/src/amici.cpp @@ -59,6 +59,7 @@ std::map simulation_status_to_str_map = { {AMICI_ERR_FAILURE, "AMICI_ERR_FAILURE"}, {AMICI_CONV_FAILURE, "AMICI_CONV_FAILURE"}, {AMICI_FIRST_RHSFUNC_ERR, "AMICI_FIRST_RHSFUNC_ERR"}, + {AMICI_CONSTR_FAIL, "AMICI_CONSTR_FAIL"}, {AMICI_RHSFUNC_FAIL, "AMICI_RHSFUNC_FAIL"}, {AMICI_ILL_INPUT, "AMICI_ILL_INPUT"}, {AMICI_ERROR, "AMICI_ERROR"}, From 41b4ce4c8c6340bfe599446732593054dd7358be Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 13 Mar 2024 19:22:38 +0100 Subject: [PATCH 02/17] Optionally include measurements in `plot_observable_trajectories` (#2381) If some `ExpData` is provided, `plot_observable_trajectories` will now also visualize the measurements. --- python/sdist/amici/plotting.py | 46 ++++++++++++++++++++++++++-------- 1 file changed, 36 insertions(+), 10 deletions(-) diff --git a/python/sdist/amici/plotting.py b/python/sdist/amici/plotting.py index 25607638d7..d27f2994ce 100644 --- a/python/sdist/amici/plotting.py +++ b/python/sdist/amici/plotting.py @@ -12,6 +12,7 @@ import seaborn as sns from matplotlib.axes import Axes +import amici from . import Model, ReturnDataView from .numpy import StrOrExpr, evaluate @@ -66,10 +67,11 @@ def plot_state_trajectories( for ix, label in zip(state_indices, labels): ax.plot(rdata["t"], rdata["x"][:, ix], marker=marker, label=label) - ax.set_xlabel("$t$") - ax.set_ylabel("$x(t)$") - ax.legend() - ax.set_title("State trajectories") + + ax.set_xlabel("$t$") + ax.set_ylabel("$x(t)$") + ax.legend() + ax.set_title("State trajectories") def plot_observable_trajectories( @@ -79,6 +81,7 @@ def plot_observable_trajectories( model: Model = None, prefer_names: bool = True, marker=None, + edata: Union[amici.ExpData, amici.ExpDataView] = None, ) -> None: """ Plot observable trajectories. @@ -97,8 +100,12 @@ def plot_observable_trajectories( :param marker: Point marker for plotting (see `matplotlib documentation `_). - + :param edata: + Experimental data to be plotted (no event observables yet). """ + if isinstance(edata, amici.amici.ExpData): + edata = amici.ExpDataView(edata) + if not ax: fig, ax = plt.subplots() if not observable_indices: @@ -125,11 +132,30 @@ def plot_observable_trajectories( labels = np.asarray(rdata.ptr.observable_ids)[list(observable_indices)] for iy, label in zip(observable_indices, labels): - ax.plot(rdata["t"], rdata["y"][:, iy], marker=marker, label=label) - ax.set_xlabel("$t$") - ax.set_ylabel("$y(t)$") - ax.legend() - ax.set_title("Observable trajectories") + (l,) = ax.plot( + rdata["t"], rdata["y"][:, iy], marker=marker, label=label + ) + + if edata is not None: + ax.plot( + edata.ts, + edata.observedData[:, iy], + "x", + label=f"exp. {label}", + color=l.get_color(), + ) + ax.errorbar( + edata.ts, + edata.observedData[:, iy], + yerr=rdata.sigmay[:, iy], + fmt="none", + color=l.get_color(), + ) + + ax.set_xlabel("$t$") + ax.set_ylabel("$y(t)$") + ax.set_title("Observable trajectories") + ax.legend() def plot_jacobian(rdata: ReturnDataView): From 790ab4427508de57ad430a6cc156e1869458bb31 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Mon, 18 Mar 2024 12:23:17 +0100 Subject: [PATCH 03/17] Fix initial state issues with PEtab (#2382) * Ensure initial state parameters are always fixed parameters * Ensure initial concentration for preequilibration and simulation are both set in both phases to avoid NaN issues * Fix initialAssignment handling during PEtab import (missed in #2359) --- python/sdist/amici/petab/parameter_mapping.py | 15 +++++++++------ python/sdist/amici/petab/sbml_import.py | 4 +++- 2 files changed, 12 insertions(+), 7 deletions(-) diff --git a/python/sdist/amici/petab/parameter_mapping.py b/python/sdist/amici/petab/parameter_mapping.py index 9c527f0395..369ad13a96 100644 --- a/python/sdist/amici/petab/parameter_mapping.py +++ b/python/sdist/amici/petab/parameter_mapping.py @@ -486,12 +486,11 @@ def create_parameter_mapping_for_condition( condition_scale_map_preeq, preeq_value, ) - else: - # need to set dummy value for preeq parameter anyways, as it - # is expected below (set to 0, not nan, because will be - # multiplied with indicator variable in initial assignment) - condition_map_sim[init_par_id] = 0.0 - condition_scale_map_sim[init_par_id] = LIN + # need to set dummy value for preeq parameter anyways, as it + # is expected below (set to 0, not nan, because will be + # multiplied with indicator variable in initial assignment) + condition_map_sim[init_par_id] = 0.0 + condition_scale_map_sim[init_par_id] = LIN # for simulation condition_id = condition[SIMULATION_CONDITION_ID] @@ -505,6 +504,10 @@ def create_parameter_mapping_for_condition( condition_scale_map_sim, value, ) + # set dummy value as above + if condition_map_preeq: + condition_map_preeq[init_par_id] = 0.0 + condition_scale_map_preeq[init_par_id] = LIN ########################################################################## # separate fixed and variable AMICI parameters, because we may have diff --git a/python/sdist/amici/petab/sbml_import.py b/python/sdist/amici/petab/sbml_import.py index e349683505..2484d57a7a 100644 --- a/python/sdist/amici/petab/sbml_import.py +++ b/python/sdist/amici/petab/sbml_import.py @@ -300,6 +300,8 @@ def import_model_sbml( init_par = sbml_model.createParameter() init_par.setId(init_par_id) init_par.setName(init_par_id) + # must be a fixed parameter in any case to allow reinitialization + fixed_parameters.append(init_par_id) assignment = sbml_model.getInitialAssignment(assignee_id) if assignment is None: assignment = sbml_model.createInitialAssignment() @@ -535,7 +537,7 @@ def _get_fixed_parameters_sbml( ia.getMath(), parser_settings ) ) - if not sym_math.is_Number: + if not sym_math.evalf().is_Number: fixed_parameters.remove(fixed_parameter) continue From ebf5eeffb12d75ad8499d03102aab056958ee10d Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 28 Mar 2024 09:46:58 +0100 Subject: [PATCH 04/17] Fix cmake: `cannot create directory: /cmake/Amici` (#2389) Fixes cmake-install failures during model import in cases where BLAS was not found via FindBLAS. Needs a prettier solution at some point, but it does the job for now. --- cmake/AmiciFindBLAS.cmake | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/cmake/AmiciFindBLAS.cmake b/cmake/AmiciFindBLAS.cmake index 6b2f9cf5f9..51d9c9b257 100644 --- a/cmake/AmiciFindBLAS.cmake +++ b/cmake/AmiciFindBLAS.cmake @@ -83,13 +83,14 @@ if(NOT TARGET BLAS::BLAS) INTERFACE_INCLUDE_DIRECTORIES "${BLAS_INCLUDE_DIRS}" INTERFACE_LINK_LIBRARIES "${BLAS_LIBRARIES}") add_library(BLAS::BLAS ALIAS BLAS) - install(TARGETS BLAS EXPORT BLAS) - export(EXPORT BLAS NAMESPACE BLAS::) - install( - EXPORT BLAS - DESTINATION "${CMAKE_INSTALL_LIBDIR}/cmake/Amici" - NAMESPACE BLAS::) - + if("${PROJECT_NAME}" STREQUAL "amici") + install(TARGETS BLAS EXPORT BLAS) + export(EXPORT BLAS NAMESPACE BLAS::) + install( + EXPORT BLAS + DESTINATION "${CMAKE_INSTALL_LIBDIR}/cmake/Amici" + NAMESPACE BLAS::) + endif() # legacy python package environment variables: if(DEFINED ENV{BLAS_CFLAGS}) From b28b83b904504f7aa6c12da1254be3f335756644 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 28 Mar 2024 10:39:36 +0100 Subject: [PATCH 05/17] Fix Solver operator== and copyctor (#2388) A few recently added members were missing. Constraints were not copied when cloning a Solver. --- src/solver.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/solver.cpp b/src/solver.cpp index 6fc1132c4d..b9e497f309 100644 --- a/src/solver.cpp +++ b/src/solver.cpp @@ -17,6 +17,7 @@ Solver::Solver(Solver const& other) , maxsteps_(other.maxsteps_) , maxtime_(other.maxtime_) , simulation_timer_(other.simulation_timer_) + , constraints_(other.constraints_) , sensi_meth_(other.sensi_meth_) , sensi_meth_preeq_(other.sensi_meth_preeq_) , stldet_(other.stldet_) @@ -564,7 +565,9 @@ bool operator==(Solver const& a, Solver const& b) { == b.check_sensi_steadystate_conv_) && (a.rdata_mode_ == b.rdata_mode_) && (a.max_conv_fails_ == b.max_conv_fails_) - && (a.max_nonlin_iters_ == b.max_nonlin_iters_); + && (a.max_nonlin_iters_ == b.max_nonlin_iters_) + && (a.max_step_size_ == b.max_step_size_) + && (a.constraints_.getVector() == b.constraints_.getVector()); } void Solver::applyTolerances() const { From d9933622f2672360fe50c3a59d75bc14a0fcdaef Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 4 Apr 2024 16:41:16 +0200 Subject: [PATCH 06/17] Update valgrind suppressions (#2390) * Update valgrind suppressions --- python/tests/valgrind-python.supp | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/python/tests/valgrind-python.supp b/python/tests/valgrind-python.supp index b4e0177033..01bc776aec 100644 --- a/python/tests/valgrind-python.supp +++ b/python/tests/valgrind-python.supp @@ -213,6 +213,19 @@ fun:__Pyx__PyObject_CallOneArg } +{ + other + Memcheck:Leak + match-leak-kinds: definite + fun:realloc + fun:resize_compact + fun:_PyUnicodeWriter_Finish + fun:PyUnicode_FromFormatV + fun:PyUnicode_FromFormat + fun:PyFortranObject_NewAsAttr + ... +} + { other Memcheck:Value8 From 6bbcde2cb12abe4286a0bc383fe3224b71fdc0db Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 17 Apr 2024 11:54:13 +0200 Subject: [PATCH 07/17] Require sphinx<7.3.0 (#2404) Closes #2403 --- documentation/rtd_requirements.txt | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/documentation/rtd_requirements.txt b/documentation/rtd_requirements.txt index 8d2c2100f9..c5b1450dc2 100644 --- a/documentation/rtd_requirements.txt +++ b/documentation/rtd_requirements.txt @@ -1,5 +1,7 @@ # NOTE: relative paths are expected to be relative to the repository root -sphinx + +# sphinx<7.3.0: https://github.com/AMICI-dev/AMICI/issues/2403 +sphinx<7.3.0 mock>=5.0.2 setuptools>=67.7.2 pysb>=1.11.0 From 51c8cad7cb6130928091dd2a22e310433f351144 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 18 Apr 2024 08:37:41 +0200 Subject: [PATCH 08/17] Revert "Require sphinx<7.3.0 (#2404)" This reverts commit 6bbcde2cb12abe4286a0bc383fe3224b71fdc0db. Fixed in sphinx 7.3.6 (https://www.sphinx-doc.org/en/master/changes.html) --- documentation/rtd_requirements.txt | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/documentation/rtd_requirements.txt b/documentation/rtd_requirements.txt index c5b1450dc2..8d2c2100f9 100644 --- a/documentation/rtd_requirements.txt +++ b/documentation/rtd_requirements.txt @@ -1,7 +1,5 @@ # NOTE: relative paths are expected to be relative to the repository root - -# sphinx<7.3.0: https://github.com/AMICI-dev/AMICI/issues/2403 -sphinx<7.3.0 +sphinx mock>=5.0.2 setuptools>=67.7.2 pysb>=1.11.0 From 18d4af316818208a1bd3ab316420b8e8c6ed4af7 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Fri, 19 Apr 2024 22:21:04 +0200 Subject: [PATCH 09/17] Model::checkFinite: exclude timepoints (#2395) Removes the check for non-finite timepoints. This rather confusing (#2370) than helpful. The current timepoint is already included in all other nan/inf messages, and the other timepoints don't matter in this context. Closes #2370. --- src/model.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/model.cpp b/src/model.cpp index cefdf1ac97..452e484f58 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -1671,7 +1671,6 @@ int Model::checkFinite( && model_quantity != ModelQuantity::ts) { checkFinite(state_.fixedParameters, ModelQuantity::k, t); checkFinite(state_.unscaledParameters, ModelQuantity::p, t); - checkFinite(simulation_parameters_.ts_, ModelQuantity::ts, t); if (!always_check_finite_ && model_quantity != ModelQuantity::w) { // don't check twice if always_check_finite_ is true checkFinite(derived_state_.w_, ModelQuantity::w, t); @@ -1789,7 +1788,6 @@ int Model::checkFinite( // check upstream checkFinite(state_.fixedParameters, ModelQuantity::k, t); checkFinite(state_.unscaledParameters, ModelQuantity::p, t); - checkFinite(simulation_parameters_.ts_, ModelQuantity::ts, t); checkFinite(derived_state_.w_, ModelQuantity::w, t); return AMICI_RECOVERABLE_ERROR; @@ -1880,7 +1878,6 @@ int Model::checkFinite(SUNMatrix m, ModelQuantity model_quantity, realtype t) // check upstream checkFinite(state_.fixedParameters, ModelQuantity::k, t); checkFinite(state_.unscaledParameters, ModelQuantity::p, t); - checkFinite(simulation_parameters_.ts_, ModelQuantity::ts, t); checkFinite(derived_state_.w_, ModelQuantity::w, t); return AMICI_RECOVERABLE_ERROR; From 0bba2361124aca5e8510fe331a2a1623cf8067a9 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Fri, 19 Apr 2024 22:22:16 +0200 Subject: [PATCH 10/17] =?UTF-8?q?More=20comprehensive/robust=20conversion?= =?UTF-8?q?=20of=20docstring=20type/rtype=20to=20type=20=E2=80=A6=20(#2401?= =?UTF-8?q?)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit …annotations Should cover everything currently present in `amici.py`, except for overloaded methods. Exclude two methods from swig interface that aren't easily usable from Python. --- python/sdist/amici/swig.py | 75 ++++++++++++++++++++++++++++---------- swig/misc.i | 2 + 2 files changed, 57 insertions(+), 20 deletions(-) diff --git a/python/sdist/amici/swig.py b/python/sdist/amici/swig.py index fbc486c301..5ba8017005 100644 --- a/python/sdist/amici/swig.py +++ b/python/sdist/amici/swig.py @@ -1,11 +1,12 @@ """Functions related to SWIG or SWIG-generated code""" +from __future__ import annotations import ast import contextlib import re class TypeHintFixer(ast.NodeTransformer): - """Replaces SWIG-generated C++ typehints by corresponding Python types""" + """Replaces SWIG-generated C++ typehints by corresponding Python types.""" mapping = { "void": None, @@ -53,9 +54,13 @@ class TypeHintFixer(ast.NodeTransformer): "std::allocator< amici::ParameterScaling > > const &": ast.Constant( "ParameterScalingVector" ), + "H5::H5File": None, } def visit_FunctionDef(self, node): + # convert type/rtype from docstring to annotation, if possible. + # those may be c++ types, not valid in python, that need to be + # converted to python types below. self._annotation_from_docstring(node) # Has a return type annotation? @@ -67,14 +72,17 @@ def visit_FunctionDef(self, node): for arg in node.args.args: if not arg.annotation: continue - if isinstance(arg.annotation, ast.Name): + if not isinstance(arg.annotation, ast.Constant): # there is already proper annotation continue arg.annotation = self._new_annot(arg.annotation.value) return node - def _new_annot(self, old_annot: str): + def _new_annot(self, old_annot: str | ast.Name): + if isinstance(old_annot, ast.Name): + old_annot = old_annot.id + with contextlib.suppress(KeyError): return self.mapping[old_annot] @@ -117,6 +125,8 @@ def _annotation_from_docstring(self, node: ast.FunctionDef): Swig sometimes generates ``:type solver: :py:class:`Solver`` instead of ``:type solver: Solver``. Those need special treatment. + + Overloaded functions are skipped. """ docstring = ast.get_docstring(node, clean=False) if not docstring or "*Overload 1:*" in docstring: @@ -127,22 +137,18 @@ def _annotation_from_docstring(self, node: ast.FunctionDef): lines_to_remove = set() for line_no, line in enumerate(docstring): - if ( - match := re.match( - r"\s*:rtype:\s*(?::py:class:`)?(\w+)`?\s+$", line - ) - ) and not match.group(1).startswith(":"): - node.returns = ast.Constant(match.group(1)) + if type_str := self.extract_rtype(line): + # handle `:rtype:` + node.returns = ast.Constant(type_str) lines_to_remove.add(line_no) + continue - if ( - match := re.match( - r"\s*:type\s*(\w+):\W*(?::py:class:`)?(\w+)`?\s+$", line - ) - ) and not match.group(1).startswith(":"): + arg_name, type_str = self.extract_type(line) + if arg_name is not None: + # handle `:type ...:` for arg in node.args.args: - if arg.arg == match.group(1): - arg.annotation = ast.Constant(match.group(2)) + if arg.arg == arg_name: + arg.annotation = ast.Constant(type_str) lines_to_remove.add(line_no) if lines_to_remove: @@ -155,13 +161,42 @@ def _annotation_from_docstring(self, node: ast.FunctionDef): ) node.body[0].value = ast.Str(new_docstring) + @staticmethod + def extract_type(line: str) -> tuple[str, str] | tuple[None, None]: + """Extract argument name and type string from ``:type:`` docstring + line.""" + match = re.match(r"\s*:type\s+(\w+):\s+(.+?)(?:, optional)?\s*$", line) + if not match: + return None, None + + arg_name = match.group(1) + + # get rid of any :py:class`...` in the type string if necessary + if not match.group(2).startswith(":py:"): + return arg_name, match.group(2) + + match = re.match(r":py:\w+:`(.+)`", match.group(2)) + assert match + return arg_name, match.group(1) + + @staticmethod + def extract_rtype(line: str) -> str | None: + """Extract type string from ``:rtype:`` docstring line.""" + match = re.match(r"\s*:rtype:\s+(.+)\s*$", line) + if not match: + return None + + # get rid of any :py:class`...` in the type string if necessary + if not match.group(1).startswith(":py:"): + return match.group(1) + + match = re.match(r":py:\w+:`(.+)`", match.group(1)) + assert match + return match.group(1) + def fix_typehints(infilename, outfilename): """Change SWIG-generated C++ typehints to Python typehints""" - # Only available from Python3.9 - if not getattr(ast, "unparse", None): - return - # file -> AST with open(infilename) as f: source = f.read() diff --git a/swig/misc.i b/swig/misc.i index 8015e28bfe..af166b48ed 100644 --- a/swig/misc.i +++ b/swig/misc.i @@ -4,6 +4,8 @@ %ignore amici::regexErrorToString; %ignore amici::writeSlice; %ignore ContextManager; +%ignore amici::scaleParameters; +%ignore amici::unscaleParameters; // Add necessary symbols to generated header %{ From 206df037fc5ab2a5804982f87e5a17a080401cf5 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Fri, 19 Apr 2024 22:23:07 +0200 Subject: [PATCH 11/17] Fix ModuleNotFoundError in tox -e doc (#2400) `import exhale_multiproject_monkeypatch` fails on some systems due to unclear sys.path issues. Therefore, we need to adjust sys.path manually. --- documentation/conf.py | 28 +++++++++++++++++++++------- 1 file changed, 21 insertions(+), 7 deletions(-) diff --git a/documentation/conf.py b/documentation/conf.py index 8b2379a299..25c6dab647 100644 --- a/documentation/conf.py +++ b/documentation/conf.py @@ -9,19 +9,33 @@ import subprocess import sys from enum import EnumType - -# need to import before setting typing.TYPE_CHECKING=True, fails otherwise -import amici import exhale.deploy -import exhale_multiproject_monkeypatch from unittest import mock -import pandas as pd import sphinx -import sympy as sp from exhale import configs as exhale_configs from sphinx.transforms.post_transforms import ReferencesResolver -exhale_multiproject_monkeypatch, pd, sp # to avoid removal of unused import +try: + import exhale_multiproject_monkeypatch # noqa: F401 +except ModuleNotFoundError: + # for unclear reasons, the import of exhale_multiproject_monkeypatch + # fails on some systems, because the the location of the editable install + # is not automatically added to sys.path ¯\_(ツ)_/¯ + from importlib.metadata import Distribution + import json + from urllib.parse import unquote_plus, urlparse + + dist = Distribution.from_name("sphinx-contrib-exhale-multiproject") + url = json.loads(dist.read_text("direct_url.json"))["url"] + package_dir = unquote_plus(urlparse(url).path) + sys.path.append(package_dir) + import exhale_multiproject_monkeypatch # noqa: F401 + +# need to import before setting typing.TYPE_CHECKING=True, fails otherwise +import amici +import pandas as pd # noqa: F401 +import sympy as sp # noqa: F401 + # BEGIN Monkeypatch exhale from exhale.deploy import _generate_doxygen as exhale_generate_doxygen From 09e0581ddc55d46e8ace6938cce33b50823c4fe0 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Fri, 19 Apr 2024 22:24:19 +0200 Subject: [PATCH 12/17] Fix incorrect exception types / messages for IDASolver (#2398) --- src/solver_idas.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/solver_idas.cpp b/src/solver_idas.cpp index 4f96c95f86..8093b336e5 100644 --- a/src/solver_idas.cpp +++ b/src/solver_idas.cpp @@ -17,7 +17,7 @@ namespace amici { /* - * The following static members are callback function to CVODES. + * The following static members are callback function to IDAS. * Their signatures must not be changes. */ @@ -437,7 +437,7 @@ void IDASolver::reInitPostProcess( auto status = IDASetStopTime(ida_mem, tout); if (status != IDA_SUCCESS) - throw IDAException(status, "CVodeSetStopTime"); + throw IDAException(status, "IDASetStopTime"); status = IDASolve( ami_mem, tout, t, yout->getNVector(), ypout->getNVector(), IDA_ONE_STEP @@ -853,7 +853,7 @@ void IDASolver::setNonLinearSolver() const { solver_memory_.get(), non_linear_solver_->get() ); if (status != IDA_SUCCESS) - throw CvodeException(status, "CVodeSetNonlinearSolver"); + throw IDAException(status, "IDASetNonlinearSolver"); } void IDASolver::setNonLinearSolverSens() const { @@ -883,7 +883,7 @@ void IDASolver::setNonLinearSolverSens() const { } if (status != IDA_SUCCESS) - throw CvodeException(status, "CVodeSolver::setNonLinearSolverSens"); + throw IDAException(status, "IDASolver::setNonLinearSolverSens"); } void IDASolver::setNonLinearSolverB(int which) const { @@ -891,7 +891,7 @@ void IDASolver::setNonLinearSolverB(int which) const { solver_memory_.get(), which, non_linear_solver_B_->get() ); if (status != IDA_SUCCESS) - throw CvodeException(status, "CVodeSetNonlinearSolverB"); + throw IDAException(status, "IDASetNonlinearSolverB"); } /** From 28b2611cc10c91f2d781944e9166b7d1e65341cf Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Fri, 19 Apr 2024 22:26:09 +0200 Subject: [PATCH 13/17] SwigPtrView: look up unhandled attributes in _swigptr (#2405) So far, only the explicitly listed attributes are accessible directly via SwigPtrView. Things like ReturnData.ny are only available through ReturnDataView.ptr.ny, which is inconvenient. Let's just look all non-private attributes that are not already explicitly handled by SwigPtrView on _swigptr and return them as is. --- python/sdist/amici/numpy.py | 17 ++++++++++------- python/sdist/amici/plotting.py | 2 +- python/tests/test_swig_interface.py | 3 +++ 3 files changed, 14 insertions(+), 8 deletions(-) diff --git a/python/sdist/amici/numpy.py b/python/sdist/amici/numpy.py index c1aef949c6..f40d0f4c6e 100644 --- a/python/sdist/amici/numpy.py +++ b/python/sdist/amici/numpy.py @@ -55,15 +55,18 @@ def __getitem__(self, item: str) -> Union[np.ndarray, float]: if item in self._cache: return self._cache[item] - if item == "id": - return getattr(self._swigptr, item) + if item in self._field_names: + value = _field_as_numpy( + self._field_dimensions, item, self._swigptr + ) + self._cache[item] = value - if item not in self._field_names: - self.__missing__(item) + return value + + if not item.startswith("_") and hasattr(self._swigptr, item): + return getattr(self._swigptr, item) - value = _field_as_numpy(self._field_dimensions, item, self._swigptr) - self._cache[item] = value - return value + self.__missing__(item) def __missing__(self, key: str) -> None: """ diff --git a/python/sdist/amici/plotting.py b/python/sdist/amici/plotting.py index d27f2994ce..19dbe05f89 100644 --- a/python/sdist/amici/plotting.py +++ b/python/sdist/amici/plotting.py @@ -109,7 +109,7 @@ def plot_observable_trajectories( if not ax: fig, ax = plt.subplots() if not observable_indices: - observable_indices = range(rdata["y"].shape[1]) + observable_indices = range(rdata.ny) if marker is None: # Show marker if only one time point is available, diff --git a/python/tests/test_swig_interface.py b/python/tests/test_swig_interface.py index f214519f26..b5063ca3cc 100644 --- a/python/tests/test_swig_interface.py +++ b/python/tests/test_swig_interface.py @@ -511,6 +511,9 @@ def test_rdataview(sbml_example_presimulation_module): rdata = amici.runAmiciSimulation(model, model.getSolver()) assert isinstance(rdata, amici.ReturnDataView) + # check that non-array attributes are looked up in the wrapped object + assert rdata.ptr.ny == rdata.ny + # fields are accessible via dot notation and [] operator, # __contains__ and __getattr__ are implemented correctly with pytest.raises(AttributeError): From e8493cf32773847078151954754b8a33eaa5047d Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Sat, 20 Apr 2024 00:49:20 +0200 Subject: [PATCH 14/17] cmake: set sundials path hint for python package (#2397) So far, CMakeLists.txt only provided a path hint for sundials for the standalone c++ library, but not for the python package. The latter should be added to ensure the right sundials version is found. There have been issues with system installations of sundials. Closes #2394 --- CMakeLists.txt | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 744847930d..d8a1109d90 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -133,9 +133,16 @@ elseif(AMICI_TRY_ENABLE_HDF5) endif() set(VENDORED_SUNDIALS_DIR ${CMAKE_CURRENT_SOURCE_DIR}/ThirdParty/sundials) -set(VENDORED_SUNDIALS_BUILD_DIR ${VENDORED_SUNDIALS_DIR}/build) -set(VENDORED_SUNDIALS_INSTALL_DIR ${VENDORED_SUNDIALS_BUILD_DIR}) set(SUNDIALS_PRIVATE_INCLUDE_DIRS "${VENDORED_SUNDIALS_DIR}/src") +# Handle different sundials build/install dirs, depending on whether we are +# building the Python extension only or the full C++ interface +if(AMICI_PYTHON_BUILD_EXT_ONLY) + set(VENDORED_SUNDIALS_BUILD_DIR ${CMAKE_CURRENT_SOURCE_DIR}) + set(VENDORED_SUNDIALS_INSTALL_DIR ${VENDORED_SUNDIALS_BUILD_DIR}) +else() + set(VENDORED_SUNDIALS_BUILD_DIR ${VENDORED_SUNDIALS_DIR}/build) + set(VENDORED_SUNDIALS_INSTALL_DIR ${VENDORED_SUNDIALS_BUILD_DIR}) +endif() find_package( SUNDIALS REQUIRED PATHS "${VENDORED_SUNDIALS_INSTALL_DIR}/${CMAKE_INSTALL_LIBDIR}/cmake/sundials/") From 01cc4e9cb00bb74a0f2f87e65331b1c2b783f885 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Sat, 20 Apr 2024 08:34:23 +0200 Subject: [PATCH 15/17] Allow subselection of state variables for steady-state simulations (#2387) * Allow subselection of state variables for steady-state simulations Closes #2368 * .. --- include/amici/model.h | 43 ++++++++++++++++++ include/amici/serialization.h | 1 + include/amici/steadystateproblem.h | 6 ++- include/amici/vector.h | 4 +- python/sdist/amici/swig_wrappers.py | 1 + python/tests/test_preequilibration.py | 63 +++++++++++++++++++++++++++ python/tests/test_swig_interface.py | 14 ++++-- src/hdf5.cpp | 6 +++ src/model.cpp | 21 ++++++++- src/steadystateproblem.cpp | 27 +++++++++--- swig/model.i | 1 + 11 files changed, 171 insertions(+), 16 deletions(-) diff --git a/include/amici/model.h b/include/amici/model.h index 29b98aa913..1c020400ac 100644 --- a/include/amici/model.h +++ b/include/amici/model.h @@ -1481,6 +1481,40 @@ class Model : public AbstractModel, public ModelDimensions { */ virtual std::vector get_trigger_timepoints() const; + /** + * @brief Get steady-state mask as std::vector. + * + * See `set_steadystate_mask` for details. + * + * @return Steady-state mask + */ + std::vector get_steadystate_mask() const { + return steadystate_mask_.getVector(); + }; + + /** + * @brief Get steady-state mask as AmiVector. + * + * See `set_steadystate_mask` for details. + * @return Steady-state mask + */ + AmiVector const& get_steadystate_mask_av() const { + return steadystate_mask_; + }; + + /** + * @brief Set steady-state mask. + * + * The mask is used to exclude certain state variables from the steady-state + * convergence check. Positive values indicate that the corresponding state + * variable should be included in the convergence check, while non-positive + * values indicate that the corresponding state variable should be excluded. + * An empty mask is interpreted as including all state variables. + * + * @param mask Mask of length `nx_solver`. + */ + void set_steadystate_mask(std::vector const& mask); + /** * Flag indicating whether for * `amici::Solver::sensi_` == `amici::SensitivityOrder::second` @@ -2087,6 +2121,15 @@ class Model : public AbstractModel, public ModelDimensions { /** Simulation parameters, initial state, etc. */ SimulationParameters simulation_parameters_; + + /** + * Mask for state variables that should be checked for steady state + * during pre-/post-equilibration. Positive values indicate that the + * corresponding state variable should be checked for steady state. + * Negative values indicate that the corresponding state variable should + * be ignored. + */ + AmiVector steadystate_mask_; }; bool operator==(Model const& a, Model const& b); diff --git a/include/amici/serialization.h b/include/amici/serialization.h index e4b1e3afe5..a03540104f 100644 --- a/include/amici/serialization.h +++ b/include/amici/serialization.h @@ -150,6 +150,7 @@ void serialize(Archive& ar, amici::Model& m, unsigned int const /*version*/) { ar & m.steadystate_computation_mode_; ar & m.steadystate_sensitivity_mode_; ar & m.state_independent_events_; + ar & m.steadystate_mask_; } /** diff --git a/include/amici/steadystateproblem.h b/include/amici/steadystateproblem.h index 55c9aaca77..72d248f8bf 100644 --- a/include/amici/steadystateproblem.h +++ b/include/amici/steadystateproblem.h @@ -246,14 +246,16 @@ class SteadystateProblem { * w_i = 1 / ( rtol * x_i + atol ) * @param x current state (sx[ip] for sensitivities) * @param xdot current rhs (sxdot[ip] for sensitivities) + * @param mask mask for state variables to include in WRMS norm. + * Positive value: include; non-positive value: exclude; empty: include all. * @param atol absolute tolerance * @param rtol relative tolerance * @param ewt error weight vector * @return root-mean-square norm */ realtype getWrmsNorm( - AmiVector const& x, AmiVector const& xdot, realtype atol, realtype rtol, - AmiVector& ewt + AmiVector const& x, AmiVector const& xdot, AmiVector const& mask, + realtype atol, realtype rtol, AmiVector& ewt ) const; /** diff --git a/include/amici/vector.h b/include/amici/vector.h index 486ea77682..27f1d81a81 100644 --- a/include/amici/vector.h +++ b/include/amici/vector.h @@ -430,8 +430,8 @@ inline span make_span(N_Vector nv) { } /** - * @brief Create span from N_Vector - * @param nv + * @brief Create span from AmiVector + * @param av * */ inline span make_span(amici::AmiVector const& av) { diff --git a/python/sdist/amici/swig_wrappers.py b/python/sdist/amici/swig_wrappers.py index e942cc76dd..72f850118c 100644 --- a/python/sdist/amici/swig_wrappers.py +++ b/python/sdist/amici/swig_wrappers.py @@ -255,6 +255,7 @@ def writeSolverSettingsToHDF5( "SteadyStateSensitivityMode", ("t0", "setT0"), "Timepoints", + "_steadystate_mask", ] diff --git a/python/tests/test_preequilibration.py b/python/tests/test_preequilibration.py index d003507199..d9a2335459 100644 --- a/python/tests/test_preequilibration.py +++ b/python/tests/test_preequilibration.py @@ -8,6 +8,10 @@ from amici.debugging import get_model_for_preeq from numpy.testing import assert_allclose, assert_equal from test_pysb import get_data +from amici.testing import ( + TemporaryDirectoryWinSafe as TemporaryDirectory, + skip_on_valgrind, +) @pytest.fixture @@ -658,3 +662,62 @@ def test_get_model_for_preeq(preeq_fixture): rdata1.sx, rdata2.sx, ) + + +@skip_on_valgrind +def test_partial_eq(): + """Check that partial equilibration is possible.""" + from amici.antimony_import import antimony2amici + + ant_str = """ + model test_partial_eq + explodes = 1 + explodes' = explodes + A = 1 + B = 0 + R: A -> B; k*A - k*B + k = 1 + end + """ + module_name = "test_partial_eq" + with TemporaryDirectory(prefix=module_name) as outdir: + antimony2amici( + ant_str, + model_name=module_name, + output_dir=outdir, + ) + model_module = amici.import_model_module( + module_name=module_name, module_path=outdir + ) + amici_model = model_module.getModel() + amici_model.setTimepoints([np.inf]) + amici_solver = amici_model.getSolver() + amici_solver.setRelativeToleranceSteadyState(1e-12) + + # equilibration of `explodes` will fail + rdata = amici.runAmiciSimulation(amici_model, amici_solver) + assert rdata.status == amici.AMICI_ERROR + assert rdata.messages[0].identifier == "EQUILIBRATION_FAILURE" + + # excluding `explodes` should enable equilibration + amici_model.set_steadystate_mask( + [ + 0 if state_id == "explodes" else 1 + for state_id in amici_model.getStateIdsSolver() + ] + ) + rdata = amici.runAmiciSimulation(amici_model, amici_solver) + assert rdata.status == amici.AMICI_SUCCESS + assert_allclose( + rdata.by_id("A"), + 0.5, + atol=amici_solver.getAbsoluteToleranceSteadyState(), + rtol=amici_solver.getRelativeToleranceSteadyState(), + ) + assert_allclose( + rdata.by_id("B"), + 0.5, + atol=amici_solver.getAbsoluteToleranceSteadyState(), + rtol=amici_solver.getRelativeToleranceSteadyState(), + ) + assert rdata.t_last < 100 diff --git a/python/tests/test_swig_interface.py b/python/tests/test_swig_interface.py index b5063ca3cc..2dfb46e0a7 100644 --- a/python/tests/test_swig_interface.py +++ b/python/tests/test_swig_interface.py @@ -39,7 +39,7 @@ def test_copy_constructors(pysb_example_presimulation_module): val = get_val(obj, attr) try: - modval = get_mod_val(val, attr) + modval = get_mod_val(val, attr, obj) except ValueError: # happens for everything that is not bool or scalar continue @@ -79,6 +79,10 @@ def test_copy_constructors(pysb_example_presimulation_module): tuple([1.0] + [0.0] * 35), tuple([0.1] * 36), ], + "_steadystate_mask": [ + (), + tuple([0] * 3), + ], "MinimumSigmaResiduals": [ 50.0, 60.0, @@ -361,19 +365,21 @@ def get_val(obj, attr): return getattr(obj, attr) -def get_mod_val(val, attr): +def get_mod_val(val, attr, obj): if attr == "getReturnDataReportingMode": return amici.RDataReporting.likelihood elif attr == "getParameterList": - return tuple(get_mod_val(val[0], "") for _ in val) + return tuple(get_mod_val(val[0], "", obj) for _ in val) elif attr == "getStateIsNonNegative": raise ValueError("Cannot modify value") + elif attr == "get_steadystate_mask": + return [0 for _ in range(obj.nx_solver)] elif isinstance(val, bool): return not val elif isinstance(val, numbers.Number): return val + 1 elif isinstance(val, tuple): - return tuple(get_mod_val(v, attr) for v in val) + return tuple(get_mod_val(v, attr, obj) for v in val) raise ValueError("Cannot modify value") diff --git a/src/hdf5.cpp b/src/hdf5.cpp index 0454b634f1..f9914452eb 100644 --- a/src/hdf5.cpp +++ b/src/hdf5.cpp @@ -1244,6 +1244,12 @@ void readModelDataFromHDF5( model.setInitialStates(x0); } + if (locationExists(file, datasetPath + "/steadystate_mask")) { + auto mask = getDoubleDataset1D(file, datasetPath + "/steadystate_mask"); + if (!mask.empty()) + model.set_steadystate_mask(mask); + } + if (locationExists(file, datasetPath + "/sx0")) { hsize_t length0 = 0; hsize_t length1 = 0; diff --git a/src/model.cpp b/src/model.cpp index 452e484f58..c27a5310d9 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -287,7 +287,9 @@ bool operator==(Model const& a, Model const& b) { && (a.nmaxevent_ == b.nmaxevent_) && (a.state_is_non_negative_ == b.state_is_non_negative_) && (a.sigma_res_ == b.sigma_res_) && (a.min_sigma_ == b.min_sigma_) - && a.state_ == b.state_; + && (a.state_ == b.state_) + && (a.steadystate_mask_.getVector() + == b.steadystate_mask_.getVector()); } bool operator==(ModelDimensions const& a, ModelDimensions const& b) { @@ -3135,6 +3137,23 @@ std::vector Model::get_trigger_timepoints() const { return trigger_timepoints; } +void Model::set_steadystate_mask(std::vector const& mask) { + if (mask.size() == 0) { + if (steadystate_mask_.getLength() != 0) { + steadystate_mask_ = AmiVector(); + } + return; + } + + if (gsl::narrow(mask.size()) != nx_solver) + throw AmiException( + "Steadystate mask has wrong size: %d, expected %d", + gsl::narrow(mask.size()), nx_solver + ); + + steadystate_mask_ = AmiVector(mask); +} + const_N_Vector Model::computeX_pos(const_N_Vector x) { if (any_state_non_negative_) { for (int ix = 0; ix < derived_state_.x_pos_tmp_.getLength(); ++ix) { diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp index 98c36589f7..d78f9a8705 100644 --- a/src/steadystateproblem.cpp +++ b/src/steadystateproblem.cpp @@ -509,8 +509,8 @@ bool SteadystateProblem::getSensitivityFlag( } realtype SteadystateProblem::getWrmsNorm( - AmiVector const& x, AmiVector const& xdot, realtype atol, realtype rtol, - AmiVector& ewt + AmiVector const& x, AmiVector const& xdot, AmiVector const& mask, + realtype atol, realtype rtol, AmiVector& ewt ) const { /* Depending on what convergence we want to check (xdot, sxdot, xQBdot) we need to pass ewt[QB], as xdot and xQBdot have different sizes */ @@ -522,7 +522,14 @@ realtype SteadystateProblem::getWrmsNorm( N_VAddConst(ewt.getNVector(), atol, ewt.getNVector()); /* ewt = 1/ewt (ewt = 1/(rtol*x+atol)) */ N_VInv(ewt.getNVector(), ewt.getNVector()); - /* wrms = sqrt(sum((xdot/ewt)**2)/n) where n = size of state vector */ + + // wrms = sqrt(sum((xdot/ewt)**2)/n) where n = size of state vector + if (mask.getLength()) { + return N_VWrmsNormMask( + const_cast(xdot.getNVector()), ewt.getNVector(), + const_cast(mask.getNVector()) + ); + } return N_VWrmsNorm( const_cast(xdot.getNVector()), ewt.getNVector() ); @@ -543,7 +550,10 @@ SteadystateProblem::getWrms(Model& model, SensitivityMethod sensi_method) { "Newton type convergence check is not implemented for adjoint " "steady state computations. Stopping." ); - wrms = getWrmsNorm(xQB_, xQBdot_, atol_quad_, rtol_quad_, ewtQB_); + wrms = getWrmsNorm( + xQB_, xQBdot_, model.get_steadystate_mask_av(), atol_quad_, + rtol_quad_, ewtQB_ + ); } else { /* If we're doing a forward simulation (with or without sensitivities: Get RHS and compute weighted error norm */ @@ -552,7 +562,8 @@ SteadystateProblem::getWrms(Model& model, SensitivityMethod sensi_method) { else updateRightHandSide(model); wrms = getWrmsNorm( - state_.x, newton_step_conv_ ? delta_ : xdot_, atol_, rtol_, ewt_ + state_.x, newton_step_conv_ ? delta_ : xdot_, + model.get_steadystate_mask_av(), atol_, rtol_, ewt_ ); } return wrms; @@ -573,8 +584,10 @@ realtype SteadystateProblem::getWrmsFSA(Model& model) { ); if (newton_step_conv_) newton_solver_->solveLinearSystem(xdot_); - wrms - = getWrmsNorm(state_.sx[ip], xdot_, atol_sensi_, rtol_sensi_, ewt_); + wrms = getWrmsNorm( + state_.sx[ip], xdot_, model.get_steadystate_mask_av(), atol_sensi_, + rtol_sensi_, ewt_ + ); /* ideally this function would report the maximum of all wrms over all ip, but for practical purposes we can just report the wrms for the first ip where we know that the convergence threshold is not diff --git a/swig/model.i b/swig/model.i index fdb8ad9c85..93a21662c7 100644 --- a/swig/model.i +++ b/swig/model.i @@ -94,6 +94,7 @@ using namespace amici; %ignore fdx_rdatadtcl; %ignore fdx_rdatadx_solver; %ignore fdsigmaydy; +%ignore get_steadystate_mask_av; %newobject amici::Model::clone; From 97642613109d3fda86e7c51a0b4dd25e71a83576 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Mon, 22 Apr 2024 09:47:50 +0200 Subject: [PATCH 16/17] Update references (#2330) * KissVen2024 * .. --- documentation/amici_refs.bib | 31 +++++++++++++++++++++++++++++++ documentation/references.md | 17 ++++++++++++++++- 2 files changed, 47 insertions(+), 1 deletion(-) diff --git a/documentation/amici_refs.bib b/documentation/amici_refs.bib index ef283eaf52..550a58ecc6 100644 --- a/documentation/amici_refs.bib +++ b/documentation/amici_refs.bib @@ -1275,6 +1275,37 @@ @Article{SluijsZho2024 publisher = {Springer Science and Business Media LLC}, } +@Article{KissVen2024, + author = {Kiss, Anna E and Venkatasubramani, Anuroop V and Pathirana, Dilan and Krause, Silke and Sparr, Aline Campos and Hasenauer, Jan and Imhof, Axel and Müller, Marisa and Becker, Peter B}, + journal = {Nucleic Acids Research}, + title = {{Processivity and specificity of histone acetylation by the male-specific lethal complex}}, + year = {2024}, + issn = {0305-1048}, + month = {02}, + pages = {gkae123}, + abstract = {{Acetylation of lysine 16 of histone H4 (H4K16ac) stands out among the histone modifications, because it decompacts the chromatin fiber. The metazoan acetyltransferase MOF (KAT8) regulates transcription through H4K16 acetylation. Antibody-based studies had yielded inconclusive results about the selectivity of MOF to acetylate the H4 N-terminus. We used targeted mass spectrometry to examine the activity of MOF in the male-specific lethal core (4-MSL) complex on nucleosome array substrates. This complex is part of the Dosage Compensation Complex (DCC) that activates X-chromosomal genes in male Drosophila. During short reaction times, MOF acetylated H4K16 efficiently and with excellent selectivity. Upon longer incubation, the enzyme progressively acetylated lysines 12, 8 and 5, leading to a mixture of oligo-acetylated H4. Mathematical modeling suggests that MOF recognizes and acetylates H4K16 with high selectivity, but remains substrate-bound and continues to acetylate more N-terminal H4 lysines in a processive manner. The 4-MSL complex lacks non-coding roX RNA, a critical component of the DCC. Remarkably, addition of RNA to the reaction non-specifically suppressed H4 oligo-acetylation in favor of specific H4K16 acetylation. Because RNA destabilizes the MSL-nucleosome interaction in vitro we speculate that RNA accelerates enzyme-substrate turn-over in vivo, thus limiting the processivity of MOF, thereby increasing specific H4K16 acetylation.}}, + creationdate = {2024-02-28T18:25:06}, + doi = {10.1093/nar/gkae123}, + eprint = {https://academic.oup.com/nar/advance-article-pdf/doi/10.1093/nar/gkae123/56756494/gkae123.pdf}, + modificationdate = {2024-02-28T18:25:06}, + url = {https://doi.org/10.1093/nar/gkae123}, +} + +@Article{DoresicGre2024, + author = {Domagoj Dore{\v s}i{\'c} and Stephan Grein and Jan Hasenauer}, + journal = {bioRxiv}, + title = {Efficient parameter estimation for ODE models of cellular processes using semi-quantitative data}, + year = {2024}, + abstract = {Quantitative dynamical models facilitate the understanding of biological processes and the prediction of their dynamics. The parameters of these models are commonly estimated from experimental data. Yet, experimental data generated from different techniques do not provide direct information about the state of the system but a non-linear (monotonic) transformation of it. For such semi-quantitative data, when this transformation is unknown, it is not apparent how the model simulations and the experimental data can be compared. Here, we propose a versatile spline-based approach for the integration of a broad spectrum of semi-quantitative data into parameter estimation. We derive analytical formulas for the gradients of the hierarchical objective function and show that this substantially increases the estimation efficiency. Subsequently, we demonstrate that the method allows for the reliable discovery of unknown measurement transformations. Furthermore, we show that this approach can significantly improve the parameter inference based on semi-quantitative data in comparison to available methods. Modelers can easily apply our method by using our implementation in the open-source Python Parameter EStimation TOolbox (pyPESTO).Competing Interest StatementThe authors have declared no competing interest.}, + creationdate = {2024-04-20T13:05:06}, + doi = {10.1101/2024.01.26.577371}, + elocation-id = {2024.01.26.577371}, + eprint = {https://www.biorxiv.org/content/early/2024/01/30/2024.01.26.577371.full.pdf}, + modificationdate = {2024-04-20T13:05:06}, + publisher = {Cold Spring Harbor Laboratory}, + url = {https://www.biorxiv.org/content/early/2024/01/30/2024.01.26.577371}, +} + @Comment{jabref-meta: databaseType:bibtex;} @Comment{jabref-meta: grouping: diff --git a/documentation/references.md b/documentation/references.md index 2164037aaf..0a2759f4eb 100644 --- a/documentation/references.md +++ b/documentation/references.md @@ -1,6 +1,6 @@ # References -List of publications using AMICI. Total number is 83. +List of publications using AMICI. Total number is 85. If you applied AMICI in your work and your publication is missing, please let us know via a new GitHub issue. @@ -14,6 +14,21 @@ If you applied AMICI in your work and your publication is missing, please let us

2024

+
+Dorešić, Domagoj, Stephan Grein, and Jan Hasenauer. 2024. +“Efficient Parameter Estimation for ODE Models of Cellular +Processes Using Semi-Quantitative Data.” bioRxiv. https://doi.org/10.1101/2024.01.26.577371. +
+
+Kiss, Anna E, Anuroop V Venkatasubramani, Dilan Pathirana, Silke Krause, +Aline Campos Sparr, Jan Hasenauer, Axel Imhof, Marisa Müller, and Peter +B Becker. 2024. “Processivity and specificity +of histone acetylation by the male-specific lethal +complex.” Nucleic Acids Research, February, +gkae123. https://doi.org/10.1093/nar/gkae123. +
Lang, Paul F., David R. Penas, Julio R. Banga, Daniel Weindl, and Bela Novak. 2024. “Reusable Rule-Based Cell Cycle Model Explains From ca9d46799e5ceaa39c74e7647eb5fba6db46c6fb Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Sat, 20 Apr 2024 08:27:19 +0200 Subject: [PATCH 17/17] Bump version, update release notes --- CHANGELOG.md | 42 ++++++++++++++++++++++++++++++++++++++++++ version.txt | 2 +- 2 files changed, 43 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7052c85b08..e658b12593 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,48 @@ ## v0.X Series +### v0.24.0 (2024-04-22) + +This will be the last release supporting Python 3.9. +Future releases will require Python 3.10. + +**Fixes** + +* Fix cmake error `cannot create directory: /cmake/Amici` + during model import in cases where BLAS was not found via `FindBLAS` + by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2389 +* Added status code `AMICI_CONSTR_FAIL` + by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2379 +* Fixed certain initial state issues with PEtab + by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2382 +* Fixed Solver `operator==` and copyctor + (constraints were not copied correctly) + by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2388 +* Avoid confusing warnings about non-finite timepoints + by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2395 +* Fixed incorrect exception types / messages for `IDASolver` + by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2398 +* cmake: set SUNDIALS path hint for python package to help CMake find + the correct SUNDIALS installation + by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2397 + +* **Features** + +* Optionally include measurements in `plot_observable_trajectories` + by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2381 +* Improved type annotations in swig-wrappers + by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2401 +* Additional attributes are accessible directly via `ReturnDataView` and + `ExpDataView`, e.g. `ReturnDataView.ny`, `ReturnDataView.nx` + by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2405 +* Allow subselection of state variables for convergence check during + steady-state simulations via `Model.set_steadystate_mask([1, 0, ..., 1])` + (positive value: check; non-positive: don't check). + by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2387 + +**Full Changelog**: https://github.com/AMICI-dev/AMICI/compare/v0.23.1...v0.24.0 + + ### v0.23.1 (2024-03-11) * Fixes installation issues related to building SuiteSparse on some systems diff --git a/version.txt b/version.txt index 610e28725b..2094a100ca 100644 --- a/version.txt +++ b/version.txt @@ -1 +1 @@ -0.23.1 +0.24.0