+
+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
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/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/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/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
diff --git a/python/sdist/amici/plotting.py b/python/sdist/amici/plotting.py
index 25607638d7..19dbe05f89 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,12 +100,16 @@ 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:
- 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,
@@ -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):
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/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 f214519f26..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")
@@ -511,6 +517,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):
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
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"},
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 cefdf1ac97..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) {
@@ -1671,7 +1673,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 +1790,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 +1880,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;
@@ -3138,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/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 {
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");
}
/**
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/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
%{
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;
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