From 3b13b048990ec4a1a2c5924913794ba8355f5a98 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 13 Jul 2021 16:15:50 +0200 Subject: [PATCH] git subrepo clone (merge) --branch=v0.11.18 --force git@github.com:AMICI-dev/AMICI.git deps/AMICI (#350) subrepo: subdir: "deps/AMICI" merged: "bd3b89d3" upstream: origin: "git@github.com:AMICI-dev/AMICI.git" branch: "v0.11.18" commit: "bd3b89d3" git-subrepo: version: "0.4.3" origin: "https://github.com/ingydotnet/git-subrepo" commit: "be9f02a" --- deps/AMICI/.gitrepo | 8 +- deps/AMICI/CHANGELOG.md | 24 +++ deps/AMICI/include/amici/defines.h | 2 + deps/AMICI/include/amici/model_dimensions.h | 6 +- deps/AMICI/include/amici/serialization.h | 19 +++ deps/AMICI/include/amici/solver.h | 50 +++++-- deps/AMICI/include/amici/solver_cvodes.h | 4 +- deps/AMICI/include/amici/solver_idas.h | 4 +- deps/AMICI/python/amici/__init__.py | 34 +++-- deps/AMICI/python/amici/custom_commands.py | 13 +- deps/AMICI/python/amici/parameter_mapping.py | 4 +- deps/AMICI/python/amici/petab_objective.py | 14 +- deps/AMICI/python/amici/pysb_import.py | 9 +- deps/AMICI/python/amici/swig.py | 66 +++++++++ deps/AMICI/python/tests/test_hdf5.py | 3 + deps/AMICI/python/tests/valgrind-python.supp | 10 ++ deps/AMICI/src/amici.cpp | 14 +- deps/AMICI/src/hdf5.cpp | 13 +- deps/AMICI/src/solver.cpp | 24 ++- deps/AMICI/src/solver_cvodes.cpp | 117 +++++++++++---- deps/AMICI/src/solver_idas.cpp | 140 +++++++++++++----- .../tests/cpputest/steadystate/tests1.cpp | 21 +++ deps/AMICI/version.txt | 2 +- 23 files changed, 479 insertions(+), 122 deletions(-) diff --git a/deps/AMICI/.gitrepo b/deps/AMICI/.gitrepo index 2cefac44d..0fd9144bf 100644 --- a/deps/AMICI/.gitrepo +++ b/deps/AMICI/.gitrepo @@ -5,8 +5,8 @@ ; [subrepo] remote = git@github.com:ICB-DCM/AMICI.git - branch = develop - commit = 8ef53c8808437cafc171d6ab9545a59d61bb83af - parent = eb33a10abde3faf2835ce846482817df06afcfe1 - cmdver = 0.4.1 + branch = v0.11.18 + commit = bd3b89d3839e6518532d547e54e7fe66ef3cfec3 + parent = 98a0746a708be9e43c7daeb54ec38bcb6709f7d9 + cmdver = 0.4.3 method = merge diff --git a/deps/AMICI/CHANGELOG.md b/deps/AMICI/CHANGELOG.md index e39c780db..90649a6fc 100644 --- a/deps/AMICI/CHANGELOG.md +++ b/deps/AMICI/CHANGELOG.md @@ -2,6 +2,30 @@ ## v0.X Series +### v0.11.18 (2021-07-12) + +New: +* Allow specifying maximum integration time via + `amici::Solver::setMaxTime()` (#1530) +* Py: Add `failfast` and `num_threads` argument to `simulate_petab` + (#1528, #1524) +* Enable typehints / static type checking for AMICI-generated model modules + (#1514) (`amici.ModelModule`, available with Python>=3.8) + +Fixes: +* Python: Fix unused argument `generate_sensitivity_code` in `pysb2amici` + (#1526) +* Python: Fix C(++) stdout redirection which could have let to deadlocks in + exotic situations (#1529) +* Py: Fixed deprecation warning (#1525) +* Py: Proper typehints for SWIG interface (#1523), enabling better static type + checking and IDE integration (available with Python>=3.9) +* C++: Fixed clang compiler warning (#1521) +* C++: Fix inherited variadic ctor in exception class (#1520) +* PEtab: More informative output for unhandled species overrides (#1519) +* Return SbmlImporter from PEtab model import (#1517) + + ### v0.11.17 (2021-05-30) Fixes: diff --git a/deps/AMICI/include/amici/defines.h b/deps/AMICI/include/amici/defines.h index 0979a72fb..7e6957617 100644 --- a/deps/AMICI/include/amici/defines.h +++ b/deps/AMICI/include/amici/defines.h @@ -68,12 +68,14 @@ constexpr int AMICI_TOO_MUCH_WORK= -1; constexpr int AMICI_TOO_MUCH_ACC= -2; constexpr int AMICI_ERR_FAILURE= -3; constexpr int AMICI_CONV_FAILURE= -4; +constexpr int AMICI_RHSFUNC_FAIL= -8; constexpr int AMICI_ILL_INPUT= -22; constexpr int AMICI_ERROR= -99; constexpr int AMICI_NO_STEADY_STATE= -81; constexpr int AMICI_DAMPING_FACTOR_ERROR= -86; constexpr int AMICI_SINGULAR_JACOBIAN= -809; constexpr int AMICI_NOT_IMPLEMENTED= -999; +constexpr int AMICI_MAX_TIME_EXCEEDED = -1000; constexpr int AMICI_SUCCESS= 0; constexpr int AMICI_DATA_RETURN= 1; constexpr int AMICI_ROOT_RETURN= 2; diff --git a/deps/AMICI/include/amici/model_dimensions.h b/deps/AMICI/include/amici/model_dimensions.h index b624d8711..ca7f0f1af 100644 --- a/deps/AMICI/include/amici/model_dimensions.h +++ b/deps/AMICI/include/amici/model_dimensions.h @@ -148,12 +148,12 @@ struct ModelDimensions { */ int ndwdw {0}; - /** Number of nonzero elements in the \f$w\f$ derivative of \f$xdot\f$ */ + /** Number of nonzero elements in the \f$ w \f$ derivative of \f$ xdot \f$ */ int ndxdotdw {0}; /** - * Number of nonzero elements in the \f$y\f$ derivative of - * \f$dJy\f$ (dimension `nytrue`) + * Number of nonzero elements in the \f$ y \f$ derivative of + * \f$ dJy \f$ (dimension `nytrue`) */ std::vector ndJydy; diff --git a/deps/AMICI/include/amici/serialization.h b/deps/AMICI/include/amici/serialization.h index ebb2ec87d..6539bc36c 100644 --- a/deps/AMICI/include/amici/serialization.h +++ b/deps/AMICI/include/amici/serialization.h @@ -9,6 +9,7 @@ #include #include #include +#include #include #include @@ -83,6 +84,24 @@ void serialize(Archive &ar, amici::Solver &s, const unsigned int /*version*/) { ar &s.cpu_time_; ar &s.cpu_timeB_; ar &s.rdata_mode_; + ar &s.maxtime_; +} + +/** + * @brief Serialize std::chrono::duration to boost archive + * @param ar Archive + * @param d Duration + */ +template +void serialize(Archive &ar, std::chrono::duration &d, const unsigned int /*version*/) { + Period tmp_period; + if (Archive::is_loading::value) { + ar &tmp_period; + d = std::chrono::duration(tmp_period); + } else { + tmp_period = d.count(); + ar &tmp_period; + } } /** diff --git a/deps/AMICI/include/amici/solver.h b/deps/AMICI/include/amici/solver.h index 74400846d..448261954 100644 --- a/deps/AMICI/include/amici/solver.h +++ b/deps/AMICI/include/amici/solver.h @@ -10,6 +10,7 @@ #include #include #include +#include namespace amici { @@ -47,6 +48,9 @@ namespace amici { */ class Solver { public: + /** Type of what is passed to Sundials solvers as user_data */ + using user_data_type = std::pair; + Solver() = default; /** @@ -468,6 +472,30 @@ class Solver { */ void setMaxSteps(long int maxsteps); + /** + * @brief Returns the maximum time allowed for integration + * @return Time in seconds + */ + double getMaxTime() const; + + /** + * @brief Set the maximum time allowed for integration + * @param maxtime Time in seconds + */ + void setMaxTime(double maxtime); + + /** + * @brief Start timer for tracking integration time + */ + void startTimer() const; + + /** + * @brief Check whether maximum integration time was exceeded + * @return True if the maximum integration time was exceeded, + * false otherwise. + */ + bool timeExceeded() const; + /** * @brief returns the maximum number of solver steps for the backward * problem @@ -1113,21 +1141,16 @@ class Solver { virtual void setErrHandlerFn() const = 0; /** - * @brief Attaches the user data instance (here this is a Model) to the - * forward problem - * - * @param model Model instance + * @brief Attaches the user data to the forward problem */ - virtual void setUserData(Model *model) const = 0; + virtual void setUserData() const = 0; /** - * @brief attaches the user data instance (here this is a Model) to the - * backward problem + * @brief attaches the user data to the backward problem * * @param which identifier of the backwards problem - * @param model Model instance */ - virtual void setUserDataB(int which, Model *model) const = 0; + virtual void setUserDataB(int which) const = 0; /** * @brief specifies the maximum number of steps for the forward @@ -1519,6 +1542,9 @@ class Solver { mutable std::vector>> solver_memory_B_; + /** Sundials user_data */ + mutable user_data_type user_data; + /** internal sensitivity method flag used to select the sensitivity solution * method. Only applies for Forward Sensitivities. */ InternalSensitivityMethod ism_ {InternalSensitivityMethod::simultaneous}; @@ -1540,6 +1566,12 @@ class Solver { /** maximum number of allowed integration steps */ long int maxsteps_ {10000}; + /** Maximum wall-time for integration in seconds */ + std::chrono::duration> maxtime_ {std::chrono::duration::max()}; + + /** Time at which solver timer was started */ + mutable std::chrono::time_point starttime_; + /** linear solver for the forward problem */ mutable std::unique_ptr linear_solver_; diff --git a/deps/AMICI/include/amici/solver_cvodes.h b/deps/AMICI/include/amici/solver_cvodes.h index 7a876fe27..1f457418d 100644 --- a/deps/AMICI/include/amici/solver_cvodes.h +++ b/deps/AMICI/include/amici/solver_cvodes.h @@ -137,9 +137,9 @@ class CVodeSolver : public Solver { void setErrHandlerFn() const override; - void setUserData(Model *model) const override; + void setUserData() const override; - void setUserDataB(int which, Model *model) const override; + void setUserDataB(int which) const override; void setMaxNumSteps(long int mxsteps) const override; diff --git a/deps/AMICI/include/amici/solver_idas.h b/deps/AMICI/include/amici/solver_idas.h index e984a3e64..331a1a620 100644 --- a/deps/AMICI/include/amici/solver_idas.h +++ b/deps/AMICI/include/amici/solver_idas.h @@ -134,9 +134,9 @@ class IDASolver : public Solver { void setErrHandlerFn() const override; - void setUserData(Model *model) const override; + void setUserData() const override; - void setUserDataB(int which, Model *model) const override; + void setUserDataB(int which) const override; void setMaxNumSteps(long int mxsteps) const override; diff --git a/deps/AMICI/python/amici/__init__.py b/deps/AMICI/python/amici/__init__.py index be8ad826c..6d5420b16 100644 --- a/deps/AMICI/python/amici/__init__.py +++ b/deps/AMICI/python/amici/__init__.py @@ -18,16 +18,13 @@ :var has_clibs: boolean indicating if this is the full package with swig interface or the raw package without -:var capture_cstdout: - context to redirect C/C++ stdout to python stdout if python stdout was - redirected (doing nothing if not redirected). """ import importlib import os import re import sys -from contextlib import suppress +from contextlib import suppress, contextmanager from types import ModuleType as ModelModule from typing import Optional, Union, Sequence, List @@ -84,14 +81,22 @@ def _imported_from_setup() -> bool: return False -# redirect C/C++ stdout to python stdout if python stdout is redirected, -# e.g. in ipython notebook -capture_cstdout = suppress -if sys.stdout != sys.__stdout__: - try: - from wurlitzer import sys_pipes as capture_cstdout - except ModuleNotFoundError: - pass +try: + from wurlitzer import sys_pipes +except ModuleNotFoundError: + sys_pipes = suppress + + +@contextmanager +def _capture_cstdout(): + """Redirect C/C++ stdout to python stdout if python stdout is redirected, + e.g. in ipython notebook""" + if sys.stdout == sys.__stdout__: + yield + else: + with sys_pipes(): + yield + # Initialize AMICI paths amici_path = _get_amici_path() @@ -188,8 +193,7 @@ def runAmiciSimulation( :returns: ReturnData object with simulation results """ - - with capture_cstdout(): + with _capture_cstdout(): rdata = amici.runAmiciSimulation(_get_ptr(solver), _get_ptr(edata), _get_ptr(model)) return numpy.ReturnDataView(rdata) @@ -235,7 +239,7 @@ def runAmiciSimulations( :returns: list of simulation results """ - with capture_cstdout(): + with _capture_cstdout(): edata_ptr_vector = amici.ExpDataPtrVector(edata_list) rdata_ptr_list = amici.runAmiciSimulations(_get_ptr(solver), edata_ptr_vector, diff --git a/deps/AMICI/python/amici/custom_commands.py b/deps/AMICI/python/amici/custom_commands.py index 585a09120..e7f36ebd2 100644 --- a/deps/AMICI/python/amici/custom_commands.py +++ b/deps/AMICI/python/amici/custom_commands.py @@ -1,6 +1,5 @@ """Custom setuptools commands for AMICI installation""" -import contextlib import glob import os import subprocess @@ -8,6 +7,7 @@ from shutil import copyfile from typing import Dict, List, Tuple +from amici.swig import fix_typehints from amici.setuptools import generate_swig_interface_files from setuptools.command.build_clib import build_clib from setuptools.command.build_ext import build_ext @@ -142,9 +142,6 @@ def run(self): log.debug("running AmiciDevelop") if not self.no_clibs: - generate_swig_interface_files( - swig_outdir=os.path.join(os.path.abspath(os.getcwd()), - "amici")) self.get_finalized_command('build_clib').run() develop.run(self) @@ -226,8 +223,11 @@ def run(self): log.info(f"copying {src} -> {dest}") copyfile(src, dest) - generate_swig_interface_files( - swig_outdir=os.path.join(build_dir, 'amici')) + swig_outdir = os.path.join(os.path.abspath(build_dir), "amici") + generate_swig_interface_files(swig_outdir=swig_outdir) + swig_py_module_path = os.path.join(swig_outdir, 'amici.py') + log.debug("updating typehints") + fix_typehints(swig_py_module_path, swig_py_module_path) # Always force recompilation. The way setuptools/distutils check for # whether sources require recompilation is not reliable and may lead @@ -326,3 +326,4 @@ def set_compiler_specific_extension_options( except AttributeError: # No compiler-specific options set pass + diff --git a/deps/AMICI/python/amici/parameter_mapping.py b/deps/AMICI/python/amici/parameter_mapping.py index e7dcc0d2e..60b50ab63 100644 --- a/deps/AMICI/python/amici/parameter_mapping.py +++ b/deps/AMICI/python/amici/parameter_mapping.py @@ -93,9 +93,9 @@ def __init__( class ParameterMapping(Sequence): - """Parameter mapping for multiple conditions. + r"""Parameter mapping for multiple conditions. - This can be used like a list of :class:`ParameterMappingForCondition`\\ s. + This can be used like a list of :class:`ParameterMappingForCondition`\ s. :param parameter_mappings: List of parameter mappings for specific conditions. diff --git a/deps/AMICI/python/amici/petab_objective.py b/deps/AMICI/python/amici/petab_objective.py index 13e33d1fb..88ae41bd5 100644 --- a/deps/AMICI/python/amici/petab_objective.py +++ b/deps/AMICI/python/amici/petab_objective.py @@ -49,7 +49,9 @@ def simulate_petab( edatas: List[AmiciExpData] = None, parameter_mapping: ParameterMapping = None, scaled_parameters: Optional[bool] = False, - log_level: int = logging.WARNING + log_level: int = logging.WARNING, + num_threads: int = 1, + failfast: bool = True, ) -> Dict[str, Any]: """Simulate PEtab model. @@ -78,6 +80,12 @@ def simulate_petab( are assumed to be in linear scale. :param log_level: Log level, see :mod:`amici.logging` module. + :param num_threads: + Number of threads to use for simulating multiple conditions + (only used if compiled with OpenMP). + :param failfast: + Returns as soon as an integration failure is encountered, skipping + any remaining simulations. :return: Dictionary of @@ -136,7 +144,9 @@ def simulate_petab( amici_model=amici_model) # Simulate - rdatas = amici.runAmiciSimulations(amici_model, solver, edata_list=edatas) + rdatas = amici.runAmiciSimulations( + amici_model, solver, edata_list=edatas, + num_threads=num_threads, failfast=failfast) # Compute total llh llh = sum(rdata['llh'] for rdata in rdatas) diff --git a/deps/AMICI/python/amici/pysb_import.py b/deps/AMICI/python/amici/pysb_import.py index ce1e1ef0b..603eb3ef2 100644 --- a/deps/AMICI/python/amici/pysb_import.py +++ b/deps/AMICI/python/amici/pysb_import.py @@ -51,7 +51,7 @@ def pysb2amici( simplify: Callable = lambda x: sp.powsimp(x, deep=True), generate_sensitivity_code: bool = True, ): - """ + r""" Generate AMICI C++ files for the provided model. .. warning:: @@ -142,6 +142,7 @@ def pysb2amici( verbose=verbose, assume_pow_positivity=assume_pow_positivity, compiler=compiler, + generate_sensitivity_code=generate_sensitivity_code ) exporter.set_name(model.name) exporter.set_paths(output_dir) @@ -297,7 +298,7 @@ def _process_pysb_expressions( sigmas: Dict[str, str], noise_distributions: Optional[Dict[str, Union[str, Callable]]] = None, ) -> None: - """ + r""" Converts pysb expressions/observables into Observables (with corresponding standard deviation SigmaY and LogLikelihood) or Expressions and adds them to the ODEModel instance @@ -306,8 +307,8 @@ def _process_pysb_expressions( pysb model :param observables: - list of names of :class`pysb.Expression`\\ s or - :class:`pysb.Observable`\\ s that are to be mapped to ODEModel + list of names of :class`pysb.Expression`\ s or + :class:`pysb.Observable`\ s that are to be mapped to ODEModel observables :param sigmas: diff --git a/deps/AMICI/python/amici/swig.py b/deps/AMICI/python/amici/swig.py index 3bc4fd7f2..da189fd48 100644 --- a/deps/AMICI/python/amici/swig.py +++ b/deps/AMICI/python/amici/swig.py @@ -1,6 +1,7 @@ """Functions for downloading/building/finding SWIG""" from typing import Tuple +import ast import os import subprocess import re @@ -77,3 +78,68 @@ def get_swig_version(swig_exe: str) -> Tuple: result.stdout.decode('utf-8')) return tuple(int(x) for x in version.split('.')) + + +class TypeHintFixer(ast.NodeTransformer): + """Replaces SWIG-generated C++ typehints by corresponding Python types""" + + mapping = { + 'void': None, + 'std::unique_ptr< amici::Solver >': 'amici.Solver', + 'amici::InternalSensitivityMethod': 'amici.InternalSensitivityMethod', + 'amici::InterpolationType': 'amici.InterpolationType', + 'amici::LinearMultistepMethod': 'amici.LinearMultistepMethod', + 'amici::LinearSolver': 'amici.LinearSolver', + 'amici::Model *': 'amici.Model', + 'amici::Model const *': 'amici.Model', + 'amici::NewtonDampingFactorMode': 'amici.NewtonDampingFactorMode', + 'amici::NonlinearSolverIteration': 'amici.NonlinearSolverIteration', + 'amici::RDataReporting': 'amici.RDataReporting', + 'amici::SensitivityMethod': 'amici.SensitivityMethod', + 'amici::SensitivityOrder': 'amici.SensitivityOrder', + 'amici::Solver *': 'amici.Solver', + 'amici::SteadyStateSensitivityMode': 'amici.SteadyStateSensitivityMode', + 'amici::realtype': 'float', + 'DoubleVector': 'numpy.ndarray', + 'IntVector': 'List[int]', + 'std::string': 'str', + 'std::string const &': 'str', + 'std::unique_ptr< amici::ExpData >': 'amici.ExpData', + 'std::unique_ptr< amici::ReturnData >': 'amici.ReturnData', + } + + def visit_FunctionDef(self, node): + # Has a return type annotation? + if node.returns: + node.returns.value = self._new_annot(node.returns.value) + + # Has arguments? + if node.args.args: + for arg in node.args.args: + if not arg.annotation: + continue + arg.annotation.value = self._new_annot(arg.annotation.value) + return node + + def _new_annot(self, old_annot): + return self.mapping.get(old_annot, old_annot) + + +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, 'r') as f: + source = f.read() + parsed_source = ast.parse(source) + + # Change AST + fixer = TypeHintFixer() + parsed_source = fixer.visit(parsed_source) + + # AST -> file + with open(outfilename, 'w') as f: + f.write(ast.unparse(parsed_source)) diff --git a/deps/AMICI/python/tests/test_hdf5.py b/deps/AMICI/python/tests/test_hdf5.py index 2141045b5..40098f921 100644 --- a/deps/AMICI/python/tests/test_hdf5.py +++ b/deps/AMICI/python/tests/test_hdf5.py @@ -21,6 +21,9 @@ def _modify_solver_attrs(solver): cval = 0 elif attr == 'setReturnDataReportingMode': cval = amici.RDataReporting.likelihood + elif attr == 'setMaxTime': + # default value is the maximum, must not add to that + cval = random.random() elif isinstance(val, int): cval = val + 1 else: diff --git a/deps/AMICI/python/tests/valgrind-python.supp b/deps/AMICI/python/tests/valgrind-python.supp index 6e23c2422..10e33a4de 100644 --- a/deps/AMICI/python/tests/valgrind-python.supp +++ b/deps/AMICI/python/tests/valgrind-python.supp @@ -79,6 +79,16 @@ fun:__pyx_pw_5numpy_6random_13bit_generator_12BitGenerator_1__init__ } +{ + numpy + Memcheck:Leak + match-leak-kinds: definite + fun:malloc + obj:/usr/bin/python3.? + ... + fun:gentype_generic_method +} + # # scipy # diff --git a/deps/AMICI/src/amici.cpp b/deps/AMICI/src/amici.cpp index 841fdb98a..4d6e57b1c 100644 --- a/deps/AMICI/src/amici.cpp +++ b/deps/AMICI/src/amici.cpp @@ -104,6 +104,8 @@ AmiciApplication::runAmiciSimulation(Solver& solver, Model& model, bool rethrow) { + solver.startTimer(); + /* Applies condition-specific model settings and restores them when going * out of scope */ ConditionContext cc1(&model, edata, FixedParameterContext::simulation); @@ -168,7 +170,11 @@ AmiciApplication::runAmiciSimulation(Solver& solver, rdata->status = AMICI_SUCCESS; } catch (amici::IntegrationFailure const& ex) { - rdata->status = ex.error_code; + if(ex.error_code == AMICI_RHSFUNC_FAIL && solver.timeExceeded()) { + rdata->status = AMICI_MAX_TIME_EXCEEDED; + } else { + rdata->status = ex.error_code; + } if (rethrow) throw; warningF("AMICI:simulation", @@ -176,7 +182,11 @@ AmiciApplication::runAmiciSimulation(Solver& solver, ex.time, ex.what()); } catch (amici::IntegrationFailureB const& ex) { - rdata->status = ex.error_code; + if(ex.error_code == AMICI_RHSFUNC_FAIL && solver.timeExceeded()) { + rdata->status = AMICI_MAX_TIME_EXCEEDED; + } else { + rdata->status = ex.error_code; + } if (rethrow) throw; warningF( diff --git a/deps/AMICI/src/hdf5.cpp b/deps/AMICI/src/hdf5.cpp index b548e4d54..626f3af24 100644 --- a/deps/AMICI/src/hdf5.cpp +++ b/deps/AMICI/src/hdf5.cpp @@ -636,6 +636,10 @@ void writeSolverSettingsToHDF5(Solver const& solver, H5LTset_attribute_double(file.getId(), hdf5Location.c_str(), "ss_rtol_sensi", &dbuffer, 1); + dbuffer = solver.getMaxTime(); + H5LTset_attribute_double(file.getId(), hdf5Location.c_str(), + "maxtime", &dbuffer, 1); + ibuffer = static_cast(solver.getMaxSteps()); H5LTset_attribute_int(file.getId(), hdf5Location.c_str(), "maxsteps", &ibuffer, 1); @@ -775,6 +779,11 @@ void readSolverSettingsFromHDF5(H5::H5File const& file, Solver &solver, "ss_rtol_sensi")); } + if(attributeExists(file, datasetPath, "maxtime")) { + solver.setMaxTime( + getDoubleScalarAttribute(file, datasetPath, "maxtime")); + } + if(attributeExists(file, datasetPath, "maxsteps")) { solver.setMaxSteps( getIntScalarAttribute(file, datasetPath, "maxsteps")); @@ -960,12 +969,12 @@ void readModelDataFromHDF5(const H5::H5File &file, Model &model, model.setUnscaledInitialStateSensitivities(sx0); } } - + if(attributeExists(file, datasetPath, "sigma_res")) { auto sigma_res = getIntScalarAttribute(file, datasetPath, "sigma_res"); model.setAddSigmaResiduals(static_cast(sigma_res)); } - + if(attributeExists(file, datasetPath, "min_sigma")) { auto min_sigma = getDoubleScalarAttribute(file, datasetPath, "min_sigma"); diff --git a/deps/AMICI/src/solver.cpp b/deps/AMICI/src/solver.cpp index 13bdf5ac1..89cb67ae0 100644 --- a/deps/AMICI/src/solver.cpp +++ b/deps/AMICI/src/solver.cpp @@ -20,6 +20,7 @@ Solver::Solver(AmiciApplication *app) : app(app) Solver::Solver(const Solver &other) : ism_(other.ism_), lmm_(other.lmm_), iter_(other.iter_), interp_type_(other.interp_type_), maxsteps_(other.maxsteps_), + maxtime_(other.maxtime_), sensi_meth_(other.sensi_meth_), sensi_meth_preeq_(other.sensi_meth_preeq_), stldet_(other.stldet_), ordering_(other.ordering_), newton_maxsteps_(other.newton_maxsteps_), @@ -133,7 +134,8 @@ void Solver::setup(const realtype t0, Model *model, const AmiVector &x0, /* Set optional inputs */ setErrHandlerFn(); /* Attaches userdata */ - setUserData(model); + user_data = std::make_pair(model, this); + setUserData(); /* activates stability limit detection */ setStabLimDet(stldet_); @@ -184,7 +186,7 @@ void Solver::setupB(int *which, const realtype tf, Model *model, binit(*which, tf, xB0, dxB0); /* Attach user data */ - setUserDataB(*which, model); + setUserDataB(*which); if (nx() == 0) return; @@ -495,6 +497,7 @@ bool operator==(const Solver &a, const Solver &b) { (a.linsol_ == b.linsol_) && (a.atol_ == b.atol_) && (a.rtol_ == b.rtol_) && (a.maxsteps_ == b.maxsteps_) && (a.maxstepsB_ == b.maxstepsB_) && (a.quad_atol_ == b.quad_atol_) && (a.quad_rtol_ == b.quad_rtol_) && + (a.maxtime_ == b.maxtime_) && (a.getAbsoluteToleranceSteadyState() == b.getAbsoluteToleranceSteadyState()) && (a.getRelativeToleranceSteadyState() == @@ -837,6 +840,23 @@ void Solver::setAbsoluteToleranceSteadyStateSensi(const double atol) { long int Solver::getMaxSteps() const { return maxsteps_; } +double Solver::getMaxTime() const { return maxtime_.count(); } + +void Solver::setMaxTime(double maxtime) +{ + maxtime_ = std::chrono::duration(maxtime); +} + +void Solver::startTimer() const +{ + starttime_ = std::chrono::system_clock::now(); +} + +bool Solver::timeExceeded() const +{ + return std::chrono::system_clock::now() - starttime_ > maxtime_; +} + void Solver::setMaxSteps(const long int maxsteps) { if (maxsteps <= 0) throw AmiException("maxsteps must be a positive number"); diff --git a/deps/AMICI/src/solver_cvodes.cpp b/deps/AMICI/src/solver_cvodes.cpp index 3b96d34af..fb0c8fd0f 100644 --- a/deps/AMICI/src/solver_cvodes.cpp +++ b/deps/AMICI/src/solver_cvodes.cpp @@ -349,14 +349,17 @@ void CVodeSolver::setErrHandlerFn() const { throw CvodeException(status, "CVodeSetErrHandlerFn"); } -void CVodeSolver::setUserData(Model *model) const { - int status = CVodeSetUserData(solver_memory_.get(), model); +void CVodeSolver::setUserData() const { + int status = CVodeSetUserData( + solver_memory_.get(), + &user_data + ); if (status != CV_SUCCESS) throw CvodeException(status, "CVodeSetUserData"); } -void CVodeSolver::setUserDataB(const int which, Model *model) const { - int status = CVodeSetUserDataB(solver_memory_.get(), which, model); +void CVodeSolver::setUserDataB(const int which) const { + int status = CVodeSetUserDataB(solver_memory_.get(), which, &user_data); if (status != CV_SUCCESS) throw CvodeException(status, "CVodeSetUserDataB"); } @@ -792,7 +795,10 @@ const Model *CVodeSolver::getModel() const { throw AmiException("Solver has not been allocated, information is not " "available"); auto cv_mem = static_cast(solver_memory_.get()); - return static_cast(cv_mem->cv_user_data); + + auto typed_udata = static_cast(cv_mem->cv_user_data); + Expects(typed_udata); + return typed_udata->first; } @@ -812,7 +818,11 @@ const Model *CVodeSolver::getModel() const { static int fJ(realtype t, N_Vector x, N_Vector xdot, SUNMatrix J, void *user_data, N_Vector /*tmp1*/, N_Vector /*tmp2*/, N_Vector /*tmp3*/) { - auto model = static_cast(user_data); + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); + model->fJ(t, x, xdot, J); return model->checkFinite(gsl::make_span(J), "Jacobian"); } @@ -835,7 +845,11 @@ static int fJ(realtype t, N_Vector x, N_Vector xdot, SUNMatrix J, static int fJB(realtype t, N_Vector x, N_Vector xB, N_Vector xBdot, SUNMatrix JB, void *user_data, N_Vector /*tmp1B*/, N_Vector /*tmp2B*/, N_Vector /*tmp3B*/) { - auto model = static_cast(user_data); + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); + model->fJB(t, x, xB, xBdot, JB); return model->checkFinite(gsl::make_span(JB), "Jacobian"); } @@ -856,7 +870,11 @@ static int fJB(realtype t, N_Vector x, N_Vector xB, N_Vector xBdot, static int fJSparse(realtype t, N_Vector x, N_Vector /*xdot*/, SUNMatrix J, void *user_data, N_Vector /*tmp1*/, N_Vector /*tmp2*/, N_Vector /*tmp3*/) { - auto model = static_cast(user_data); + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); + model->fJSparse(t, x, J); return model->checkFinite(gsl::make_span(J), "Jacobian"); } @@ -878,7 +896,11 @@ static int fJSparse(realtype t, N_Vector x, N_Vector /*xdot*/, static int fJSparseB(realtype t, N_Vector x, N_Vector xB, N_Vector xBdot, SUNMatrix JB, void *user_data, N_Vector /*tmp1B*/, N_Vector /*tmp2B*/, N_Vector /*tmp3B*/) { - auto model = static_cast(user_data); + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); + model->fJSparseB(t, x, xB, xBdot, JB); return model->checkFinite(gsl::make_span(JB), "Jacobian"); } @@ -886,9 +908,6 @@ static int fJSparseB(realtype t, N_Vector x, N_Vector xB, N_Vector xBdot, /** * @brief J in banded form (for banded solvers) - * @param N number of states - * @param mupper upper matrix bandwidth - * @param mlower lower matrix bandwidth * @param t timepoint * @param x Vector with the states * @param xdot Vector with the right hand side @@ -907,9 +926,6 @@ static int fJBand(realtype t, N_Vector x, N_Vector xdot, SUNMatrix J, /** * @brief JB in banded form (for banded solvers) - * @param NeqBdot number of states - * @param mupper upper matrix bandwidth - * @param mlower lower matrix bandwidth * @param t timepoint * @param x Vector with the states * @param xB Vector with the adjoint states @@ -942,7 +958,11 @@ static int fJBandB(realtype t, N_Vector x, N_Vector xB, N_Vector xBdot, **/ static int fJv(N_Vector v, N_Vector Jv, realtype t, N_Vector x, N_Vector /*xdot*/, void *user_data, N_Vector /*tmp*/) { - auto model = static_cast(user_data); + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); + model->fJv(v, Jv, t, x); return model->checkFinite(gsl::make_span(Jv), "Jacobian"); } @@ -964,7 +984,11 @@ static int fJv(N_Vector v, N_Vector Jv, realtype t, N_Vector x, static int fJvB(N_Vector vB, N_Vector JvB, realtype t, N_Vector x, N_Vector xB, N_Vector /*xBdot*/, void *user_data, N_Vector /*tmpB*/) { - auto model = static_cast(user_data); + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); + model->fJvB(vB, JvB, t, x, xB); return model->checkFinite(gsl::make_span(JvB), "Jacobian"); } @@ -980,7 +1004,11 @@ static int fJvB(N_Vector vB, N_Vector JvB, realtype t, N_Vector x, */ static int froot(realtype t, N_Vector x, realtype *root, void *user_data) { - auto model = static_cast(user_data); + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); + model->froot(t, x, gsl::make_span(root, model->ne)); return model->checkFinite(gsl::make_span(root, model->ne), "root function"); @@ -996,7 +1024,16 @@ static int froot(realtype t, N_Vector x, realtype *root, * @return status flag indicating successful execution */ static int fxdot(realtype t, N_Vector x, N_Vector xdot, void *user_data) { - auto model = static_cast(user_data); + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); + auto solver = dynamic_cast(typed_udata->second); + Expects(model); + + if(solver->timeExceeded()) { + return AMICI_MAX_TIME_EXCEEDED; + } if (t > 1e200 && !model->checkFinite(gsl::make_span(x), "fxdot")) { /* when t is large (typically ~1e300), CVODES may pass all NaN x @@ -1022,7 +1059,17 @@ static int fxdot(realtype t, N_Vector x, N_Vector xdot, void *user_data) { */ static int fxBdot(realtype t, N_Vector x, N_Vector xB, N_Vector xBdot, void *user_data) { - auto model = static_cast(user_data); + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); + auto solver = dynamic_cast(typed_udata->second); + Expects(model); + + if(solver->timeExceeded()) { + return AMICI_MAX_TIME_EXCEEDED; + } + model->fxBdot(t, x, xB, xBdot); return model->checkFinite(gsl::make_span(xBdot), "fxBdot"); } @@ -1039,7 +1086,11 @@ static int fxBdot(realtype t, N_Vector x, N_Vector xB, N_Vector xBdot, */ static int fqBdot(realtype t, N_Vector x, N_Vector xB, N_Vector qBdot, void *user_data) { - auto model = static_cast(user_data); + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); + model->fqBdot(t, x, xB, qBdot); return model->checkFinite(gsl::make_span(qBdot), "qBdot"); } @@ -1056,7 +1107,11 @@ static int fqBdot(realtype t, N_Vector x, N_Vector xB, N_Vector qBdot, */ static int fxBdot_ss(realtype t, N_Vector xB, N_Vector xBdot, void *user_data) { - auto model = static_cast(user_data); + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); + model->fxBdot_ss(t, xB, xBdot); return model->checkFinite(gsl::make_span(xBdot), "fxBdot_ss"); } @@ -1073,7 +1128,11 @@ static int fxBdot_ss(realtype t, N_Vector xB, N_Vector xBdot, */ static int fqBdot_ss(realtype t, N_Vector xB, N_Vector qBdot, void *user_data) { - auto model = static_cast(user_data); + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); + model->fqBdot_ss(t, xB, qBdot); return model->checkFinite(gsl::make_span(qBdot), "qBdot_ss"); } @@ -1093,7 +1152,11 @@ static int fqBdot_ss(realtype t, N_Vector xB, N_Vector qBdot, static int fJSparseB_ss(realtype /*t*/, N_Vector /*x*/, N_Vector xBdot, SUNMatrix JB, void *user_data, N_Vector /*tmp1*/, N_Vector /*tmp2*/, N_Vector /*tmp3*/) { - auto model = static_cast(user_data); + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); + model->fJSparseB_ss(JB); return model->checkFinite(gsl::make_span(xBdot), "JSparseB_ss"); } @@ -1117,7 +1180,11 @@ static int fJSparseB_ss(realtype /*t*/, N_Vector /*x*/, N_Vector xBdot, static int fsxdot(int /*Ns*/, realtype t, N_Vector x, N_Vector /*xdot*/, int ip, N_Vector sx, N_Vector sxdot, void *user_data, N_Vector /*tmp1*/, N_Vector /*tmp2*/) { - auto model = static_cast(user_data); + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); + model->fsxdot(t, x, ip, sx, sxdot); return model->checkFinite(gsl::make_span(sxdot), "sxdot"); } diff --git a/deps/AMICI/src/solver_idas.cpp b/deps/AMICI/src/solver_idas.cpp index 964cdb9ef..18033f42c 100644 --- a/deps/AMICI/src/solver_idas.cpp +++ b/deps/AMICI/src/solver_idas.cpp @@ -286,14 +286,14 @@ void IDASolver::setErrHandlerFn() const { throw IDAException(status, "IDASetErrHandlerFn"); } -void IDASolver::setUserData(Model *model) const { - int status = IDASetUserData(solver_memory_.get(), model); +void IDASolver::setUserData() const { + int status = IDASetUserData(solver_memory_.get(), &user_data); if (status != IDA_SUCCESS) throw IDAException(status, "IDASetUserData"); } -void IDASolver::setUserDataB(int which, Model *model) const { - int status = IDASetUserDataB(solver_memory_.get(), which, model); +void IDASolver::setUserDataB(int which) const { + int status = IDASetUserDataB(solver_memory_.get(), which, &user_data); if (status != IDA_SUCCESS) throw IDAException(status, "IDASetUserDataB"); } @@ -729,7 +729,10 @@ const Model *IDASolver::getModel() const { throw AmiException( "Solver has not been allocated, information is not available"); auto ida_mem = static_cast(solver_memory_.get()); - return static_cast(ida_mem->ida_user_data); + auto user_data = static_cast(ida_mem->ida_user_data); + if(user_data) + return user_data->first; + return nullptr; } void IDASolver::setLinearSolver() const { @@ -797,7 +800,7 @@ void IDASolver::setNonLinearSolverB(int which) const { * @param dx Vector with the derivative states * @param xdot Vector with the right hand side * @param J Matrix to which the Jacobian will be written - * @param user_data object with user input @type Model_DAE + * @param user_data object with user input * @param tmp1 temporary storage vector * @param tmp2 temporary storage vector * @param tmp3 temporary storage vector @@ -806,15 +809,16 @@ void IDASolver::setNonLinearSolverB(int which) const { int fJ(realtype t, realtype cj, N_Vector x, N_Vector dx, N_Vector xdot, SUNMatrix J, void *user_data, N_Vector /*tmp1*/, N_Vector /*tmp2*/, N_Vector /*tmp3*/) { - - auto model = static_cast(user_data); + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); model->fJ(t, cj, x, dx, xdot, J); return model->checkFinite(gsl::make_span(J), "Jacobian"); } /** * @brief Jacobian of xBdot with respect to adjoint state xB - * @param NeqBdot number of adjoint state variables * @param t timepoint * @param cj scaling factor, inverse of the step size * @param x Vector with the states @@ -833,8 +837,11 @@ int fJB(realtype t, realtype cj, N_Vector x, N_Vector dx, N_Vector xB, N_Vector dxB, N_Vector /*xBdot*/, SUNMatrix JB, void *user_data, N_Vector /*tmp1B*/, N_Vector /*tmp2B*/, N_Vector /*tmp3B*/) { + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); - auto model = static_cast(user_data); model->fJB(t, cj, x, dx, xB, dxB, JB); return model->checkFinite(gsl::make_span(JB), "Jacobian"); } @@ -847,7 +854,7 @@ int fJB(realtype t, realtype cj, N_Vector x, N_Vector dx, * @param dx Vector with the derivative states * @param xdot Vector with the right hand side * @param J Matrix to which the Jacobian will be written - * @param user_data object with user input @type Model_DAE + * @param user_data object with user input * @param tmp1 temporary storage vector * @param tmp2 temporary storage vector * @param tmp3 temporary storage vector @@ -857,7 +864,11 @@ int fJSparse(realtype t, realtype cj, N_Vector x, N_Vector dx, N_Vector /*xdot*/, SUNMatrix J, void *user_data, N_Vector /*tmp1*/, N_Vector /*tmp2*/, N_Vector /*tmp3*/) { - auto model = static_cast(user_data); + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); + model->fJSparse(t, cj, x, dx, J); return model->checkFinite(gsl::make_span(J), "Jacobian"); } @@ -872,7 +883,7 @@ int fJSparse(realtype t, realtype cj, N_Vector x, N_Vector dx, * @param dxB Vector with the adjoint derivative states * @param xBdot Vector with the adjoint right hand side * @param JB Matrix to which the Jacobian will be written - * @param user_data object with user input @type Model_DAE + * @param user_data object with user input * @param tmp1B temporary storage vector * @param tmp2B temporary storage vector * @param tmp3B temporary storage vector @@ -882,23 +893,24 @@ int fJSparseB(realtype t, realtype cj, N_Vector x, N_Vector dx, N_Vector xB, N_Vector dxB, N_Vector /*xBdot*/, SUNMatrix JB, void *user_data, N_Vector /*tmp1B*/, N_Vector /*tmp2B*/, N_Vector /*tmp3B*/) { - auto model = static_cast(user_data); + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); + model->fJSparseB(t, cj, x, dx, xB, dxB, JB); return model->checkFinite(gsl::make_span(JB), "Jacobian"); } /** * @brief J in banded form (for banded solvers) - * @param N number of states - * @param mupper upper matrix bandwidth - * @param mlower lower matrix bandwidth * @param t timepoint * @param cj scalar in Jacobian (inverse stepsize) * @param x Vector with the states * @param dx Vector with the derivative states * @param xdot Vector with the right hand side * @param J Matrix to which the Jacobian will be written - * @param user_data object with user input @type Model_DAE + * @param user_data object with user input * @param tmp1 temporary storage vector * @param tmp2 temporary storage vector * @param tmp3 temporary storage vector @@ -912,9 +924,6 @@ int fJBand(realtype t, realtype cj, N_Vector x, N_Vector dx, /** * @brief JB in banded form (for banded solvers) - * @param NeqBdot number of states - * @param mupper upper matrix bandwidth - * @param mlower lower matrix bandwidth * @param t timepoint * @param cj scalar in Jacobian (inverse stepsize) * @param x Vector with the states @@ -923,7 +932,7 @@ int fJBand(realtype t, realtype cj, N_Vector x, N_Vector dx, * @param dxB Vector with the adjoint derivative states * @param xBdot Vector with the adjoint right hand side * @param JB Matrix to which the Jacobian will be written - * @param user_data object with user input @type Model_DAE + * @param user_data object with user input * @param tmp1B temporary storage vector * @param tmp2B temporary storage vector * @param tmp3B temporary storage vector @@ -946,7 +955,7 @@ int fJBandB(realtype t, realtype cj, N_Vector x, N_Vector dx, * @param xdot Vector with the right hand side * @param v Vector with which the Jacobian is multiplied * @param Jv Vector to which the Jacobian vector product will be written - * @param user_data object with user input @type Model_DAE + * @param user_data object with user input * @param tmp1 temporary storage vector * @param tmp2 temporary storage vector * @return status flag indicating successful execution @@ -955,7 +964,11 @@ int fJv(realtype t, N_Vector x, N_Vector dx, N_Vector /*xdot*/, N_Vector v, N_Vector Jv, realtype cj, void *user_data, N_Vector /*tmp1*/, N_Vector /*tmp2*/) { - auto model = static_cast(user_data); + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); + model->fJv(t, x, dx, v, Jv, cj); return model->checkFinite(gsl::make_span(Jv), "Jacobian"); } @@ -971,7 +984,7 @@ int fJv(realtype t, N_Vector x, N_Vector dx, N_Vector /*xdot*/, * @param vB Vector with which the Jacobian is multiplied * @param JvB Vector to which the Jacobian vector product will be written * @param cj scalar in Jacobian (inverse stepsize) - * @param user_data object with user input @type Model_DAE + * @param user_data object with user input * @param tmpB1 temporary storage vector * @param tmpB2 temporary storage vector * @return status flag indicating successful execution @@ -981,7 +994,11 @@ int fJvB(realtype t, N_Vector x, N_Vector dx, N_Vector xB, realtype cj, void *user_data, N_Vector /*tmpB1*/, N_Vector /*tmpB2*/) { - auto model = static_cast(user_data); + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); + model->fJvB(t, x, dx, xB, dxB, vB, JvB, cj); return model->checkFinite(gsl::make_span(JvB), "Jacobian"); } @@ -992,12 +1009,16 @@ int fJvB(realtype t, N_Vector x, N_Vector dx, N_Vector xB, * @param x Vector with the states * @param dx Vector with the derivative states * @param root array with root function values - * @param user_data object with user input @type Model_DAE + * @param user_data object with user input * @return status flag indicating successful execution */ int froot(realtype t, N_Vector x, N_Vector dx, realtype *root, void *user_data) { - auto model = static_cast(user_data); + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); + model->froot(t, x, dx, gsl::make_span(root, model->ne)); return model->checkFinite(gsl::make_span(root, model->ne), "root function"); @@ -1009,12 +1030,21 @@ int froot(realtype t, N_Vector x, N_Vector dx, realtype *root, * @param x Vector with the states * @param dx Vector with the derivative states * @param xdot Vector with the right hand side - * @param user_data object with user input @type Model_DAE + * @param user_data object with user input * @return status flag indicating successful execution */ int fxdot(realtype t, N_Vector x, N_Vector dx, N_Vector xdot, void *user_data) { - auto model = static_cast(user_data); + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); + auto solver = dynamic_cast(typed_udata->second); + Expects(model); + + if(solver->timeExceeded()) { + return AMICI_MAX_TIME_EXCEEDED; + } if (t > 1e200 && !model->app->checkFinite(gsl::make_span(x), "fxdot")) { /* when t is large (typically ~1e300), CVODES may pass all NaN x @@ -1036,13 +1066,22 @@ int fxdot(realtype t, N_Vector x, N_Vector dx, N_Vector xdot, * @param xB Vector with the adjoint states * @param dxB Vector with the adjoint derivative states * @param xBdot Vector with the adjoint right hand side - * @param user_data object with user input @type Model_DAE + * @param user_data object with user input * @return status flag indicating successful execution */ int fxBdot(realtype t, N_Vector x, N_Vector dx, N_Vector xB, N_Vector dxB, N_Vector xBdot, void *user_data) { + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); + auto solver = dynamic_cast(typed_udata->second); + Expects(model); + + if(solver->timeExceeded()) { + return AMICI_MAX_TIME_EXCEEDED; + } - auto model = static_cast(user_data); model->fxBdot(t, x, dx, xB, dxB, xBdot); return model->checkFinite(gsl::make_span(xBdot), "xBdot"); } @@ -1055,13 +1094,17 @@ int fxBdot(realtype t, N_Vector x, N_Vector dx, N_Vector xB, * @param xB Vector with the adjoint states * @param dxB Vector with the adjoint derivative states * @param qBdot Vector with the adjoint quadrature right hand side - * @param user_data pointer to temp data object @type Model_DAE + * @param user_data pointer to temp data object * @return status flag indicating successful execution */ int fqBdot(realtype t, N_Vector x, N_Vector dx, N_Vector xB, N_Vector dxB, N_Vector qBdot, void *user_data) { - auto model = static_cast(user_data); + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); + model->fqBdot(t, x, dx, xB, dxB, qBdot); return model->checkFinite(gsl::make_span(qBdot), "qBdot"); @@ -1075,12 +1118,16 @@ int fqBdot(realtype t, N_Vector x, N_Vector dx, N_Vector xB, * @param xB Vector with the adjoint states * @param dxB Vector with the adjoint derivative states * @param xBdot Vector with the adjoint right hand side - * @param user_data object with user input @type Model_DAE + * @param user_data object with user input * @return status flag indicating successful execution */ static int fxBdot_ss(realtype t, N_Vector xB, N_Vector dxB, N_Vector xBdot, void *user_data) { - auto model = static_cast(user_data); + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); + model->fxBdot_ss(t, xB, dxB, xBdot); return model->checkFinite(gsl::make_span(xBdot), "xBdot_ss"); } @@ -1098,7 +1145,11 @@ static int fxBdot_ss(realtype t, N_Vector xB, N_Vector dxB, N_Vector xBdot, */ static int fqBdot_ss(realtype t, N_Vector xB, N_Vector dxB, N_Vector qBdot, void *user_data) { - auto model = static_cast(user_data); + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); + model->fqBdot_ss(t, xB, dxB, qBdot); return model->checkFinite(gsl::make_span(qBdot), "qBdot_ss"); } @@ -1111,7 +1162,7 @@ static int fqBdot_ss(realtype t, N_Vector xB, N_Vector dxB, N_Vector qBdot, * @param dx Vector with the derivative states * @param xdot Vector with the right hand side * @param J Matrix to which the Jacobian will be written - * @param user_data object with user input @type Model_DAE + * @param user_data object with user input * @param tmp1 temporary storage vector * @param tmp2 temporary storage vector * @param tmp3 temporary storage vector @@ -1121,7 +1172,11 @@ static int fqBdot_ss(realtype t, N_Vector xB, N_Vector dxB, N_Vector qBdot, N_Vector /*dx*/, N_Vector xBdot, SUNMatrix JB, void *user_data, N_Vector /*tmp1*/, N_Vector /*tmp2*/, N_Vector /*tmp3*/) { - auto model = static_cast(user_data); + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); + model->fJSparseB_ss(JB); return model->checkFinite(gsl::make_span(xBdot), "JSparseB_ss"); } @@ -1136,7 +1191,7 @@ static int fqBdot_ss(realtype t, N_Vector xB, N_Vector dxB, N_Vector qBdot, * @param sx Vector with the state sensitivities * @param sdx Vector with the derivative state sensitivities * @param sxdot Vector with the sensitivity right hand side - * @param user_data object with user input @type Model_DAE + * @param user_data object with user input * @param tmp1 temporary storage vector * @param tmp2 temporary storage vector * @param tmp3 temporary storage vector @@ -1147,7 +1202,10 @@ int fsxdot(int /*Ns*/, realtype t, N_Vector x, N_Vector dx, N_Vector *sxdot, void *user_data, N_Vector /*tmp1*/, N_Vector /*tmp2*/, N_Vector /*tmp3*/) { - auto model = static_cast(user_data); + auto typed_udata = static_cast(user_data); + Expects(typed_udata); + auto model = dynamic_cast(typed_udata->first); + Expects(model); for (int ip = 0; ip < model->nplist(); ip++) { model->fsxdot(t, x, dx, ip, sx[ip], sdx[ip], sxdot[ip]); diff --git a/deps/AMICI/tests/cpputest/steadystate/tests1.cpp b/deps/AMICI/tests/cpputest/steadystate/tests1.cpp index 3c9992df8..78be0a2ad 100644 --- a/deps/AMICI/tests/cpputest/steadystate/tests1.cpp +++ b/deps/AMICI/tests/cpputest/steadystate/tests1.cpp @@ -113,6 +113,27 @@ TEST(groupSteadystate, testRethrow) runAmiciSimulation(*solver, nullptr, *model, true)); } + +TEST(groupSteadystate, testMaxtime) +{ + auto model = amici::generic_model::getModel(); + auto solver = model->getSolver(); + + amici::hdf5::readModelDataFromHDF5( + NEW_OPTION_FILE, *model, "/model_steadystate/nosensi/options"); + amici::hdf5::readSolverSettingsFromHDF5( + NEW_OPTION_FILE, *solver, "/model_steadystate/nosensi/options"); + + auto rdata = runAmiciSimulation(*solver, nullptr, *model); + CHECK_EQUAL(amici::AMICI_SUCCESS, rdata->status); + + solver->setMaxTime(0.000001); + + // must throw + rdata = runAmiciSimulation(*solver, nullptr, *model); + CHECK_EQUAL(amici::AMICI_MAX_TIME_EXCEEDED, rdata->status); +} + TEST(groupSteadystate, testInitialStatesNonEmpty) { auto model = amici::generic_model::getModel(); diff --git a/deps/AMICI/version.txt b/deps/AMICI/version.txt index c12244635..168a8f63a 100644 --- a/deps/AMICI/version.txt +++ b/deps/AMICI/version.txt @@ -1 +1 @@ -0.11.17 +0.11.18