+role="list">
+
Loos, C., C. Marr, F. J. Theis, and J. Hasenauer. 2015.
“Computational Methods in Systems Biology.” In, edited by
O. Roux and J. Bourdon, 9308:52–63. Lecture Notes in Computer Science.
diff --git a/python/examples/example_errors.ipynb b/python/examples/example_errors.ipynb
index 2b35964d8b..3bf81aff61 100644
--- a/python/examples/example_errors.ipynb
+++ b/python/examples/example_errors.ipynb
@@ -77,7 +77,7 @@
"source": [
"petab_problem = benchmark_models_petab.get_problem(\"Fujita_SciSignal2010\")\n",
"amici_model = import_petab_problem(\n",
- " petab_problem, verbose=False, force_compile=False\n",
+ " petab_problem, verbose=False, compile_=None\n",
")\n",
"\n",
"np.random.seed(2991)\n",
@@ -423,7 +423,7 @@
"source": [
"petab_problem = benchmark_models_petab.get_problem(\"Weber_BMC2015\")\n",
"amici_model = import_petab_problem(\n",
- " petab_problem, verbose=False, force_compile=False\n",
+ " petab_problem, verbose=False, compile_=None\n",
")\n",
"\n",
"np.random.seed(4)\n",
@@ -699,7 +699,7 @@
"amici_model = import_petab_problem(\n",
" petab_problem,\n",
" verbose=False,\n",
- " force_compile=True,\n",
+ " compile_=True,\n",
" model_name=\"Blasi_CellSystems2016_1\",\n",
")\n",
"\n",
@@ -819,7 +819,7 @@
" # we need a different model name if we import the model again\n",
" # we cannot load a model with the same name as an already loaded model\n",
" model_name=\"Blasi_CellSystems2016_2\",\n",
- " force_compile=True,\n",
+ " compile_=True,\n",
")\n",
"del os.environ[\"AMICI_EXPERIMENTAL_SBML_NONCONST_CLS\"]\n",
"\n",
diff --git a/python/examples/example_jax/ExampleJax.ipynb b/python/examples/example_jax/ExampleJax.ipynb
index 225bd13667..62391ce9be 100644
--- a/python/examples/example_jax/ExampleJax.ipynb
+++ b/python/examples/example_jax/ExampleJax.ipynb
@@ -265,7 +265,7 @@
"from amici.petab.petab_import import import_petab_problem\n",
"\n",
"amici_model = import_petab_problem(\n",
- " petab_problem, force_compile=True, verbose=False\n",
+ " petab_problem, compile_=True, verbose=False\n",
")"
]
},
diff --git a/python/sdist/MANIFEST.in b/python/sdist/MANIFEST.in
index 762987665f..5be0b1ebc6 100644
--- a/python/sdist/MANIFEST.in
+++ b/python/sdist/MANIFEST.in
@@ -10,7 +10,7 @@ include version.txt
include LICENSE.md
exclude amici/*.so
exclude amici/*.dll
-prune **/build/
+prune **/build
prune amici/share
prune amici/lib
prune amici/ThirdParty/SuiteSparse/lib
diff --git a/python/sdist/amici/de_export.py b/python/sdist/amici/de_export.py
index 4e7e0999f2..f2badbea76 100644
--- a/python/sdist/amici/de_export.py
+++ b/python/sdist/amici/de_export.py
@@ -30,7 +30,6 @@
Union,
)
from collections.abc import Sequence
-
import numpy as np
import sympy as sp
from sympy.matrices.dense import MutableDenseMatrix
@@ -768,7 +767,7 @@ def __init__(
self._expressions: list[Expression] = []
self._conservation_laws: list[ConservationLaw] = []
self._events: list[Event] = []
- self.splines = []
+ self._splines = []
self._symboldim_funs: dict[str, Callable[[], int]] = {
"sx": self.num_states_solver,
"v": self.num_states_solver,
@@ -968,7 +967,7 @@ def import_from_sbml_importer(
value=spline_expr,
)
)
- self.splines = si.splines
+ self._splines = si.splines
# get symbolic expression from SBML importers
symbols = copy.copy(si.symbols)
@@ -1117,11 +1116,10 @@ def transform_dxdt_to_concentration(species_id, dxdt):
for llh in si.symbols[SymbolId.LLHY].values()
)
- self._process_sbml_rate_of(
- symbols
- ) # substitute SBML-rateOf constructs
+ # substitute SBML-rateOf constructs
+ self._process_sbml_rate_of()
- def _process_sbml_rate_of(self, symbols) -> None:
+ def _process_sbml_rate_of(self) -> None:
"""Substitute any SBML-rateOf constructs in the model equations"""
rate_of_func = sp.core.function.UndefinedFunction("rateOf")
species_sym_to_xdot = dict(zip(self.sym("x"), self.sym("xdot")))
@@ -1129,8 +1127,6 @@ def _process_sbml_rate_of(self, symbols) -> None:
def get_rate(symbol: sp.Symbol):
"""Get rate of change of the given symbol"""
- nonlocal symbols
-
if symbol.find(rate_of_func):
raise SBMLException("Nesting rateOf() is not allowed.")
@@ -1142,6 +1138,7 @@ def get_rate(symbol: sp.Symbol):
return 0
# replace rateOf-instances in xdot by xdot symbols
+ made_substitutions = False
for i_state in range(len(self.eq("xdot"))):
if rate_ofs := self._eqs["xdot"][i_state].find(rate_of_func):
self._eqs["xdot"][i_state] = self._eqs["xdot"][i_state].subs(
@@ -1151,9 +1148,14 @@ def get_rate(symbol: sp.Symbol):
for rate_of in rate_ofs
}
)
- # substitute in topological order
- subs = toposort_symbols(dict(zip(self.sym("xdot"), self.eq("xdot"))))
- self._eqs["xdot"] = smart_subs_dict(self.eq("xdot"), subs)
+ made_substitutions = True
+
+ if made_substitutions:
+ # substitute in topological order
+ subs = toposort_symbols(
+ dict(zip(self.sym("xdot"), self.eq("xdot")))
+ )
+ self._eqs["xdot"] = smart_subs_dict(self.eq("xdot"), subs)
# replace rateOf-instances in x0 by xdot equation
for i_state in range(len(self.eq("x0"))):
@@ -1165,9 +1167,55 @@ def get_rate(symbol: sp.Symbol):
}
)
+ # replace rateOf-instances in w by xdot equation
+ # here we may need toposort, as xdot may depend on w
+ made_substitutions = False
+ for i_expr in range(len(self.eq("w"))):
+ if rate_ofs := self._eqs["w"][i_expr].find(rate_of_func):
+ self._eqs["w"][i_expr] = self._eqs["w"][i_expr].subs(
+ {
+ rate_of: get_rate(rate_of.args[0])
+ for rate_of in rate_ofs
+ }
+ )
+ made_substitutions = True
+
+ if made_substitutions:
+ # Sort expressions in self._expressions, w symbols, and w equations
+ # in topological order. Ideally, this would already happen before
+ # adding the expressions to the model, but at that point, we don't
+ # have access to xdot yet.
+ # NOTE: elsewhere, conservations law expressions are expected to
+ # occur before any other w expressions, so we must maintain their
+ # position
+ # toposort everything but conservation law expressions,
+ # then prepend conservation laws
+ w_sorted = toposort_symbols(
+ dict(
+ zip(
+ self.sym("w")[self.num_cons_law() :, :],
+ self.eq("w")[self.num_cons_law() :, :],
+ )
+ )
+ )
+ w_sorted = (
+ dict(
+ zip(
+ self.sym("w")[: self.num_cons_law(), :],
+ self.eq("w")[: self.num_cons_law(), :],
+ )
+ )
+ | w_sorted
+ )
+ old_syms = tuple(self._syms["w"])
+ topo_expr_syms = tuple(w_sorted.keys())
+ new_order = [old_syms.index(s) for s in topo_expr_syms]
+ self._expressions = [self._expressions[i] for i in new_order]
+ self._syms["w"] = sp.Matrix(topo_expr_syms)
+ self._eqs["w"] = sp.Matrix(list(w_sorted.values()))
+
for component in chain(
self.observables(),
- self.expressions(),
self.events(),
self._algebraic_equations,
):
@@ -1690,7 +1738,7 @@ def _generate_symbol(self, name: str) -> None:
# placeholders for the numeric spline values.
# Need to create symbols
self._syms[name] = sp.Matrix(
- [[f"spl_{isp}" for isp in range(len(self.splines))]]
+ [[f"spl_{isp}" for isp in range(len(self._splines))]]
)
return
elif name == "sspl":
@@ -1698,7 +1746,7 @@ def _generate_symbol(self, name: str) -> None:
self._syms[name] = sp.Matrix(
[
[f"sspl_{isp}_{ip}" for ip in range(len(self._syms["p"]))]
- for isp in range(len(self.splines))
+ for isp in range(len(self._splines))
]
)
return
@@ -2050,7 +2098,7 @@ def _compute_equation(self, name: str) -> None:
elif name == "spline_values":
# force symbols
self._eqs[name] = sp.Matrix(
- [y for spline in self.splines for y in spline.values_at_nodes]
+ [y for spline in self._splines for y in spline.values_at_nodes]
)
elif name == "spline_slopes":
@@ -2058,7 +2106,7 @@ def _compute_equation(self, name: str) -> None:
self._eqs[name] = sp.Matrix(
[
d
- for spline in self.splines
+ for spline in self._splines
for d in (
sp.zeros(len(spline.derivatives_at_nodes), 1)
if spline.derivatives_by_fd
@@ -2210,6 +2258,18 @@ def _compute_equation(self, name: str) -> None:
self._eqs[name] = self.sym(name)
elif name == "dwdx":
+ if (
+ expected := list(
+ map(
+ ConservationLaw.get_x_rdata,
+ reversed(self.conservation_laws()),
+ )
+ )
+ ) != (actual := self.eq("w")[: self.num_cons_law()]):
+ raise AssertionError(
+ "Conservation laws are not at the beginning of 'w'. "
+ f"Got {actual}, expected {expected}."
+ )
x = self.sym("x")
self._eqs[name] = sp.Matrix(
[
@@ -2724,6 +2784,9 @@ def _collect_heaviside_roots(
elif arg.has(sp.Heaviside):
root_funs.extend(self._collect_heaviside_roots(arg.args))
+ if not root_funs:
+ return []
+
# substitute 'w' expressions into root expressions now, to avoid
# rewriting 'root.cpp' and 'stau.cpp' headers
# to include 'w.h'
@@ -2892,7 +2955,7 @@ def __init__(
self.model: DEModel = de_model
self.model._code_printer.known_functions.update(
splines.spline_user_functions(
- self.model.splines, self._get_index("p")
+ self.model._splines, self._get_index("p")
)
)
@@ -2972,6 +3035,7 @@ def _generate_c_code(self) -> None:
self._write_c_make_file()
self._write_swig_files()
self._write_module_setup()
+ _write_gitignore(Path(self.model_path))
shutil.copy(
CXX_MAIN_TEMPLATE_FILE, os.path.join(self.model_path, "main.cpp")
@@ -3553,14 +3617,14 @@ def _get_function_body(
return [line for line in lines if line]
def _get_create_splines_body(self):
- if not self.model.splines:
+ if not self.model._splines:
return [" return {};"]
ind4 = " " * 4
ind8 = " " * 8
body = ["return {"]
- for ispl, spline in enumerate(self.model.splines):
+ for ispl, spline in enumerate(self.model._splines):
if isinstance(spline.nodes, splines.UniformGrid):
nodes = (
f"{ind8}{{{spline.nodes.start}, {spline.nodes.stop}}}, "
@@ -3674,7 +3738,7 @@ def _write_model_header_cpp(self) -> None:
"NEVENT": self.model.num_events(),
"NEVENT_SOLVER": self.model.num_events_solver(),
"NOBJECTIVE": "1",
- "NSPL": len(self.model.splines),
+ "NSPL": len(self.model._splines),
"NW": len(self.model.sym("w")),
"NDWDP": len(
self.model.sparsesym(
@@ -4322,3 +4386,17 @@ def _parallel_applyfunc(obj: sp.Matrix, func: Callable) -> sp.Matrix:
"to a module-level function or disable parallelization by "
"setting `AMICI_IMPORT_NPROCS=1`."
) from e
+
+
+def _write_gitignore(dest_dir: Path) -> None:
+ """Write .gitignore file.
+
+ Generate a `.gitignore` file to ignore a model directory.
+
+ :param dest_dir:
+ Path to the directory to write the `.gitignore` file to.
+ """
+ dest_dir.mkdir(exist_ok=True, parents=True)
+
+ with open(dest_dir / ".gitignore", "w") as f:
+ f.write("**")
diff --git a/python/sdist/amici/numpy.py b/python/sdist/amici/numpy.py
index fdf802147c..6e2966ba4b 100644
--- a/python/sdist/amici/numpy.py
+++ b/python/sdist/amici/numpy.py
@@ -359,7 +359,8 @@ def by_id(
) or self._swigptr.parameter_ids
else:
raise NotImplementedError(
- f"Subsetting {field} by ID is not implemented or not possible."
+ f"Subsetting `{field}` by ID (`{entity_id}`) "
+ "is not implemented or not possible."
)
col_index = ids.index(entity_id)
return getattr(self, field)[:, ..., col_index]
diff --git a/python/sdist/amici/petab/petab_import.py b/python/sdist/amici/petab/petab_import.py
index 558f51ca15..902bf837ef 100644
--- a/python/sdist/amici/petab/petab_import.py
+++ b/python/sdist/amici/petab/petab_import.py
@@ -35,7 +35,7 @@ def import_petab_problem(
petab_problem: petab.Problem,
model_output_dir: Union[str, Path, None] = None,
model_name: str = None,
- force_compile: bool = False,
+ compile_: bool = None,
non_estimated_parameters_as_constants=True,
**kwargs,
) -> "amici.Model":
@@ -53,9 +53,10 @@ def import_petab_problem(
Name of the generated model module. Defaults to the ID of the model
or the model file name without the extension.
- :param force_compile:
- Whether to compile the model even if the target folder is not empty,
- or the model exists already.
+ :param compile_:
+ If ``True``, the model will be compiled. If ``False``, the model will
+ not be compiled. If ``None``, the model will be compiled if it cannot
+ be imported.
:param non_estimated_parameters_as_constants:
Whether parameters marked as non-estimated in PEtab should be
@@ -71,6 +72,17 @@ def import_petab_problem(
:return:
The imported model.
"""
+ if "force_compile" in kwargs:
+ if kwargs["force_compile"]:
+ compile_ = True
+ del kwargs["force_compile"]
+ warn(
+ "The `force_compile` option is deprecated, please use the "
+ "new `compile_` option, which also supports 'do not compile'.",
+ DeprecationWarning,
+ stacklevel=2,
+ )
+
if petab_problem.model.type_id not in (MODEL_TYPE_SBML, MODEL_TYPE_PYSB):
raise NotImplementedError(
"Unsupported model type " + petab_problem.model.type_id
@@ -107,12 +119,15 @@ def import_petab_problem(
os.makedirs(model_output_dir)
# check if compilation necessary
- if force_compile or not _can_import_model(model_name, model_output_dir):
+ if compile_ or (
+ compile_ is None
+ and not _can_import_model(model_name, model_output_dir)
+ ):
# check if folder exists
- if os.listdir(model_output_dir) and not force_compile:
+ if os.listdir(model_output_dir) and not compile_:
raise ValueError(
f"Cannot compile to {model_output_dir}: not empty. "
- "Please assign a different target or set `force_compile`."
+ "Please assign a different target or set `compile_` to `True`."
)
# remove folder if exists
diff --git a/python/tests/test_sbml_import_special_functions.py b/python/tests/test_sbml_import_special_functions.py
index 9d8f447511..3f8383ce94 100644
--- a/python/tests/test_sbml_import_special_functions.py
+++ b/python/tests/test_sbml_import_special_functions.py
@@ -12,7 +12,11 @@
from amici.antimony_import import antimony2amici
from amici.gradient_check import check_derivatives
from amici.testing import TemporaryDirectoryWinSafe, skip_on_valgrind
-from numpy.testing import assert_approx_equal, assert_array_almost_equal_nulp
+from numpy.testing import (
+ assert_approx_equal,
+ assert_array_almost_equal_nulp,
+ assert_allclose,
+)
from scipy.special import loggamma
@@ -222,3 +226,51 @@ def test_rateof():
assert_array_almost_equal_nulp(
rdata.by_id("p2"), 1 + rdata.by_id("S1")
)
+
+
+@skip_on_valgrind
+def test_rateof_with_expression_dependent_rate():
+ """Test rateOf, where the rateOf argument depends on `w` and requires
+ toposorting."""
+ ant_model = """
+ model test_rateof_with_expression_dependent_rate
+ S1 = 0;
+ S2 = 0;
+ S1' = rate;
+ S2' = 2 * rateOf(S1);
+ # the id of the following expression must be alphabetically before
+ # `rate`, so that toposort is required to evaluate the expressions
+ # in the correct order
+ e1 := 2 * rateOf(S1);
+ rate := time
+ end
+ """
+ module_name = "test_rateof_with_expression_dependent_rate"
+ with TemporaryDirectoryWinSafe(prefix=module_name) as outdir:
+ antimony2amici(
+ ant_model,
+ 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()
+ t = np.linspace(0, 10, 11)
+ amici_model.setTimepoints(t)
+ amici_solver = amici_model.getSolver()
+ rdata = amici.runAmiciSimulation(amici_model, amici_solver)
+
+ state_ids_solver = amici_model.getStateIdsSolver()
+
+ assert_array_almost_equal_nulp(rdata.by_id("e1"), 2 * t, 1)
+
+ i_S1 = state_ids_solver.index("S1")
+ i_S2 = state_ids_solver.index("S2")
+ assert_approx_equal(rdata["xdot"][i_S1], t[-1])
+ assert_approx_equal(rdata["xdot"][i_S2], 2 * t[-1])
+
+ assert_allclose(np.diff(rdata.by_id("S1")), t[:-1] + 0.5, atol=1e-9)
+ assert_array_almost_equal_nulp(
+ rdata.by_id("S2"), 2 * rdata.by_id("S1"), 10
+ )
diff --git a/src/CMakeLists.template.cmake b/src/CMakeLists.template.cmake
index 572efede80..ee50a551e8 100644
--- a/src/CMakeLists.template.cmake
+++ b/src/CMakeLists.template.cmake
@@ -41,7 +41,12 @@ message(STATUS "Found AMICI ${Amici_DIR}")
# Debug build?
if("$ENV{ENABLE_AMICI_DEBUGGING}" OR "$ENV{ENABLE_GCOV_COVERAGE}")
- add_compile_options(-UNDEBUG -O0 -g)
+ add_compile_options(-UNDEBUG)
+ if(MSVC)
+ add_compile_options(-DEBUG)
+ else()
+ add_compile_options(-O0 -g)
+ endif()
set(CMAKE_BUILD_TYPE "Debug")
endif()
diff --git a/src/model_ode.cpp b/src/model_ode.cpp
index b3d517b36e..52275008e3 100644
--- a/src/model_ode.cpp
+++ b/src/model_ode.cpp
@@ -14,8 +14,6 @@ void Model_ODE::fJ(
void Model_ODE::fJ(
realtype t, const_N_Vector x, const_N_Vector /*xdot*/, SUNMatrix J
) {
- auto x_pos = computeX_pos(x);
- fdwdx(t, N_VGetArrayPointerConst(x_pos));
fJSparse(t, x, derived_state_.J_.get());
derived_state_.J_.refresh();
auto JDense = SUNMatrixWrapper(J);
diff --git a/tests/petab_test_suite/test_petab_suite.py b/tests/petab_test_suite/test_petab_suite.py
index 0924c09576..6d08cfb693 100755
--- a/tests/petab_test_suite/test_petab_suite.py
+++ b/tests/petab_test_suite/test_petab_suite.py
@@ -62,7 +62,7 @@ def _test_case(case, model_type, version):
petab_problem=problem,
model_output_dir=model_output_dir,
model_name=model_name,
- force_compile=True,
+ compile_=True,
)
solver = model.getSolver()
solver.setSteadyStateToleranceFactor(1.0)
diff --git a/version.txt b/version.txt
index 59dad104b0..2157409059 100644
--- a/version.txt
+++ b/version.txt
@@ -1 +1 @@
-0.21.2
+0.22.0