diff --git a/.github/workflows/deploy_release.yml b/.github/workflows/deploy_release.yml index 57eacbe91d..543ad17caa 100644 --- a/.github/workflows/deploy_release.yml +++ b/.github/workflows/deploy_release.yml @@ -33,3 +33,28 @@ jobs: password: ${{ secrets.pypi_password }} packages_dir: python/sdist/dist + bioSimulatorsUpdateCliAndDockerImage: + name: Release to BioSimulators + needs: pypi + runs-on: ubuntu-latest + env: + # Owner/repository-id for the GitHub repository for the downstream command-line interface and Docker image + DOWNSTREAM_REPOSITORY: biosimulators/Biosimulators_AMICI + + # Username/token to use the GitHub API to trigger an action on the GitHub repository for the downstream + # command-line interface and Docker image. Tokens can be generated at https://github.com/settings/tokens. + # The token should have the scope `repo` + GH_ISSUE_USERNAME: ${{ secrets.BIOSIMULATORS_USERNAME }} + GH_ISSUE_TOKEN: ${{ secrets.BIOSIMULATORS_TOKEN }} + steps: + - name: Trigger GitHub action that will build and release the downstream command-line interface and Docker image + run: | + PACKAGE_VERSION="${GITHUB_REF/refs\/tags\/v/}" + WORKFLOW_FILE=ci.yml + + curl \ + -X POST \ + -u ${GH_ISSUE_USERNAME}:${GH_ISSUE_TOKEN} \ + -H "Accept: application/vnd.github.v3+json" \ + https://api.github.com/repos/${DOWNSTREAM_REPOSITORY}/actions/workflows/${WORKFLOW_FILE}/dispatches \ + -d "{\"ref\": \"dev\", \"inputs\": {\"simulatorVersion\": \"${PACKAGE_VERSION}\", \"simulatorVersionLatest\": \"true\"}}" diff --git a/.github/workflows/release_biosimulators.yml b/.github/workflows/release_biosimulators.yml deleted file mode 100644 index 53d5a55595..0000000000 --- a/.github/workflows/release_biosimulators.yml +++ /dev/null @@ -1,32 +0,0 @@ -name: Release to BioSimulators - -on: - release: - types: - - published - -jobs: - updateCliAndDockerImage: - name: Build and release downstream command-line interface and Docker image - runs-on: ubuntu-latest - env: - # Owner/repository-id for the GitHub repository for the downstream command-line interface and Docker image - DOWNSTREAM_REPOSITORY: biosimulators/Biosimulators_AMICI - - # Username/token to use the GitHub API to trigger an action on the GitHub repository for the downstream - # command-line interface and Docker image. Tokens can be generated at https://github.com/settings/tokens. - # The token should have the scope `repo` - GH_ISSUE_USERNAME: ${{ secrets.BIOSIMULATORS_USERNAME }} - GH_ISSUE_TOKEN: ${{ secrets.BIOSIMULATORS_TOKEN }} - steps: - - name: Trigger GitHub action that will build and release the downstream command-line interface and Docker image - run: | - PACKAGE_VERSION="${GITHUB_REF/refs\/tags\/v/}" - WORKFLOW_FILE=ci.yml - - curl \ - -X POST \ - -u ${GH_ISSUE_USERNAME}:${GH_ISSUE_TOKEN} \ - -H "Accept: application/vnd.github.v3+json" \ - https://api.github.com/repos/${DOWNSTREAM_REPOSITORY}/actions/workflows/${WORKFLOW_FILE}/dispatches \ - -d "{\"ref\": \"dev\", \"inputs\": {\"simulatorVersion\": \"${PACKAGE_VERSION}\", \"simulatorVersionLatest\": \"true\"}}" diff --git a/.github/workflows/test_windows.yml b/.github/workflows/test_windows.yml index 50a2cd1c44..bf6c8d80df 100644 --- a/.github/workflows/test_windows.yml +++ b/.github/workflows/test_windows.yml @@ -33,7 +33,7 @@ jobs: python -m pip install --upgrade pip \ && pip install pytest petab \ && choco install -y ninja \ - && choco install -y swig --version=4.0.1 + && choco install -y swig - name: Install OpenBLAS shell: powershell diff --git a/CHANGELOG.md b/CHANGELOG.md index e39c780dbd..90649a6fc0 100644 --- a/CHANGELOG.md +++ b/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/include/amici/defines.h b/include/amici/defines.h index 0979a72fb2..7e69576179 100644 --- a/include/amici/defines.h +++ b/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/include/amici/exception.h b/include/amici/exception.h index 6556e5fcd4..5d767fd411 100644 --- a/include/amici/exception.h +++ b/include/amici/exception.h @@ -15,6 +15,13 @@ namespace amici { */ class AmiException : public std::exception { public: + /** + * @brief Constructor with printf style interface + * @param fmt error message with printf format + * @param ... printf formatting variables + */ + AmiException(); + /** * @brief Constructor with printf style interface * @param fmt error message with printf format @@ -40,6 +47,14 @@ class AmiException : public std::exception { */ void storeBacktrace(int nMaxFrames); + protected: + /** + * @brief Store the provided message + * @param fmt error message with printf format + * @param argptr pointer to variadic argument list + */ + void storeMessage(const char *fmt, va_list argptr); + private: std::array msg_; std::array trace_; @@ -131,7 +146,13 @@ class IntegrationFailureB : public AmiException { */ class SetupFailure : public AmiException { public: - using AmiException::AmiException; + /** + * @brief Constructor with printf style interface + * @param fmt error message with printf format + * @param ... printf formatting variables + */ + explicit SetupFailure(char const* fmt, ...); + }; diff --git a/include/amici/model_dimensions.h b/include/amici/model_dimensions.h index b624d87118..ca7f0f1af0 100644 --- a/include/amici/model_dimensions.h +++ b/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/include/amici/serialization.h b/include/amici/serialization.h index ebb2ec87dc..6539bc36ca 100644 --- a/include/amici/serialization.h +++ b/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/include/amici/solver.h b/include/amici/solver.h index 74400846d1..4482619540 100644 --- a/include/amici/solver.h +++ b/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/include/amici/solver_cvodes.h b/include/amici/solver_cvodes.h index 7a876fe277..1f457418d5 100644 --- a/include/amici/solver_cvodes.h +++ b/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/include/amici/solver_idas.h b/include/amici/solver_idas.h index e984a3e64a..331a1a6206 100644 --- a/include/amici/solver_idas.h +++ b/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/python/amici/__init__.py b/python/amici/__init__.py index bb474d89d2..6d5420b163 100644 --- a/python/amici/__init__.py +++ b/python/amici/__init__.py @@ -18,17 +18,14 @@ :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 types import ModuleType +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() @@ -134,6 +139,17 @@ def _imported_from_setup() -> bool: from .sbml_import import SbmlImporter, assignmentRules2observables from .ode_export import ODEModel, ODEExporter + try: + # Requires Python>=3.8 + from typing import Protocol + + class ModelModule(Protocol): + """Enable Python static type checking for AMICI-generated model modules""" + def getModel(self) -> amici.Model: + pass + except ImportError: + pass + hdf5_enabled = 'readSolverSettingsFromHDF5' in dir() @@ -177,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) @@ -224,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, @@ -283,7 +298,7 @@ def __exit__(self, exc_type, exc_value, traceback): def import_model_module(module_name: str, - module_path: Optional[str] = None) -> ModuleType: + module_path: Optional[str] = None) -> ModelModule: """ Import Python module of an AMICI model diff --git a/python/amici/custom_commands.py b/python/amici/custom_commands.py index 585a091202..e7f36ebd2d 100644 --- a/python/amici/custom_commands.py +++ b/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/python/amici/parameter_mapping.py b/python/amici/parameter_mapping.py index e7dcc0d2ec..60b50ab63e 100644 --- a/python/amici/parameter_mapping.py +++ b/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/python/amici/petab_import.py b/python/amici/petab_import.py index bfbbf744ad..4f8ad93c9e 100644 --- a/python/amici/petab_import.py +++ b/python/amici/petab_import.py @@ -378,7 +378,7 @@ def import_model_sbml( model_output_dir: Optional[str] = None, verbose: Optional[Union[bool, int]] = True, allow_reinit_fixpar_initcond: bool = True, - **kwargs) -> None: + **kwargs) -> amici.SbmlImporter: """ Create AMICI model from PEtab problem @@ -415,6 +415,9 @@ def import_model_sbml( :param kwargs: Additional keyword arguments to be passed to :meth:`amici.sbml_import.SbmlImporter.sbml2amici`. + + :return: + The created :class:`amici.sbml_import.SbmlImporter` instance. """ set_log_level(logger, verbose) @@ -584,6 +587,8 @@ def import_model_sbml( verbose=verbose, **kwargs) + return sbml_importer + # for backwards compatibility import_model = import_model_sbml diff --git a/python/amici/petab_objective.py b/python/amici/petab_objective.py index 6ae8ab0ca1..88ae41bd5b 100644 --- a/python/amici/petab_objective.py +++ b/python/amici/petab_objective.py @@ -17,6 +17,7 @@ import numpy as np import pandas as pd import petab +import sympy as sp from petab.C import * # noqa: F403 from . import AmiciModel, AmiciExpData @@ -48,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. @@ -77,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 @@ -135,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) @@ -346,11 +357,21 @@ def _set_initial_concentration(condition_id, species_id, init_par_id, value = petab.to_float_if_float( petab_problem.condition_df.loc[condition_id, species_id]) if pd.isna(value): - value = float( - get_species_initial( - petab_problem.sbml_model.getSpecies(species_id) - ) + value = get_species_initial( + petab_problem.sbml_model.getSpecies(species_id) ) + try: + value = float(value) + except (ValueError, TypeError): + if sp.nsimplify(value).is_Atom: + # Get rid of multiplication with one + value = sp.nsimplify(value) + else: + raise NotImplementedError( + "Cannot handle non-trivial expressions for " + f"species initial for {species_id}: {value}") + # this should be a parameter ID + value = str(value) logger.debug(f'The species {species_id} has no initial value ' f'defined for the condition {condition_id} in ' 'the PEtab conditions table. The initial value is ' @@ -728,4 +749,4 @@ def rdatas_to_simulation_df( df = rdatas_to_measurement_df(rdatas=rdatas, model=model, measurement_df=measurement_df) - return df.rename(columns={MEASUREMENT: SIMULATION}) \ No newline at end of file + return df.rename(columns={MEASUREMENT: SIMULATION}) diff --git a/python/amici/pysb_import.py b/python/amici/pysb_import.py index ce1e1ef0b0..603eb3ef23 100644 --- a/python/amici/pysb_import.py +++ b/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/python/amici/swig.py b/python/amici/swig.py index 3bc4fd7f2b..da189fd48d 100644 --- a/python/amici/swig.py +++ b/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/python/tests/test_hdf5.py b/python/tests/test_hdf5.py index 2141045b50..40098f9212 100644 --- a/python/tests/test_hdf5.py +++ b/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/python/tests/valgrind-python.supp b/python/tests/valgrind-python.supp index 6e23c24221..10e33a4de1 100644 --- a/python/tests/valgrind-python.supp +++ b/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/src/CMakeLists.template.cmake b/src/CMakeLists.template.cmake index a89fc82940..6aa13400be 100644 --- a/src/CMakeLists.template.cmake +++ b/src/CMakeLists.template.cmake @@ -18,10 +18,13 @@ set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_POSITION_INDEPENDENT_CODE ON) include(CheckCXXCompilerFlag) -set(MY_CXX_FLAGS -Wall -Wno-unused-function -Wno-unused-variable -Wno-unused-but-set-variable) +set(MY_CXX_FLAGS -Wall -Wno-unused-function -Wno-unused-variable) +if (CMAKE_CXX_COMPILER_ID STREQUAL "GNU") + list(APPEND MY_CXX_FLAGS -Wno-unused-but-set-variable) +endif() foreach(FLAG ${MY_CXX_FLAGS}) unset(CUR_FLAG_SUPPORTED CACHE) - CHECK_CXX_COMPILER_FLAG(${FLAG} CUR_FLAG_SUPPORTED) + CHECK_CXX_COMPILER_FLAG(-Werror ${FLAG} CUR_FLAG_SUPPORTED) if(${CUR_FLAG_SUPPORTED}) set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${FLAG}") endif() diff --git a/src/amici.cpp b/src/amici.cpp index 841fdb98a7..4d6e57b1c8 100644 --- a/src/amici.cpp +++ b/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/src/exception.cpp b/src/exception.cpp index 97676beabd..2d1f72ecdc 100644 --- a/src/exception.cpp +++ b/src/exception.cpp @@ -8,12 +8,18 @@ namespace amici { -AmiException::AmiException(const char *fmt, ...) { +AmiException::AmiException() +{ + storeBacktrace(12); +} + +AmiException::AmiException(const char *fmt, ...) + : AmiException() +{ va_list ap; va_start(ap, fmt); - vsnprintf(msg_.data(), msg_.size(), fmt, ap); + storeMessage(fmt, ap); va_end(ap); - storeBacktrace(12); } const char *AmiException::what() const noexcept { @@ -29,6 +35,11 @@ void AmiException::storeBacktrace(const int nMaxFrames) { backtraceString(nMaxFrames).c_str()); } +void AmiException::storeMessage(const char *fmt, va_list argptr) +{ + vsnprintf(msg_.data(), msg_.size(), fmt, argptr); +} + CvodeException::CvodeException(const int error_code, const char *function) : AmiException("Cvode routine %s failed with error code %i",function,error_code){} @@ -48,4 +59,12 @@ NewtonFailure::NewtonFailure(int code, const char *function) : error_code = code; } +SetupFailure::SetupFailure(const char *fmt, ...) +{ + va_list ap; + va_start(ap, fmt); + storeMessage(fmt, ap); + va_end(ap); +} + } // namespace amici diff --git a/src/hdf5.cpp b/src/hdf5.cpp index b548e4d54f..626f3af249 100644 --- a/src/hdf5.cpp +++ b/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/src/solver.cpp b/src/solver.cpp index 13bdf5ac14..89cb67ae07 100644 --- a/src/solver.cpp +++ b/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/src/solver_cvodes.cpp b/src/solver_cvodes.cpp index 3b96d34af4..fb0c8fd0f6 100644 --- a/src/solver_cvodes.cpp +++ b/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/src/solver_idas.cpp b/src/solver_idas.cpp index 964cdb9efa..18033f42c9 100644 --- a/src/solver_idas.cpp +++ b/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/tests/cpputest/steadystate/tests1.cpp b/tests/cpputest/steadystate/tests1.cpp index 3c9992df89..78be0a2adc 100644 --- a/tests/cpputest/steadystate/tests1.cpp +++ b/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/version.txt b/version.txt index c12244635e..168a8f63a5 100644 --- a/version.txt +++ b/version.txt @@ -1 +1 @@ -0.11.17 +0.11.18