From de891bf0cb44aaa6d497086e0fce73dfeed736ef Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20Fr=C3=B6hlich?= Date: Wed, 20 Apr 2022 13:33:15 -0400 Subject: [PATCH 01/10] fix steps saving and apply check frequency to nonFSA check (#1780) --- src/steadystateproblem.cpp | 36 +++++++++++------------------------- 1 file changed, 11 insertions(+), 25 deletions(-) diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp index df5104dd9e..d7959085f0 100644 --- a/src/steadystateproblem.cpp +++ b/src/steadystateproblem.cpp @@ -629,7 +629,8 @@ void SteadystateProblem::runSteadystateSimulation(const Solver &solver, /* If run after Newton's method checks again if it converged */ wrms_ = getWrms(model, sensitivityFlag); - int sim_steps = 0; + + int &sim_steps = backward ? numstepsB_ : numsteps_.at(1); int convergence_check_frequency = 1; @@ -653,22 +654,16 @@ void SteadystateProblem::runSteadystateSimulation(const Solver &solver, flagUpdatedState(); } - try { /* Check for convergence */ - wrms_ = getWrms(model, sensitivityFlag); - /* getWrms needs to be called before getWrmsFSA such that the linear - system is prepared for newton type convergence check */ - if (wrms_ < conv_thresh && check_sensi_conv_ && - sensitivityFlag == SensitivityMethod::forward && - sim_steps % convergence_check_frequency == 0) { - updateSensiSimulation(solver); - wrms_ = getWrmsFSA(model); - } - - } catch (NewtonFailure const &) { - /* linear solves in getWrms failed */ - numsteps_.at(1) = sim_steps; - throw; + if (sim_steps % convergence_check_frequency == 0) { + wrms_ = getWrms(model, sensitivityFlag); + /* getWrms needs to be called before getWrmsFSA such that the linear + system is prepared for newton type convergence check */ + if (wrms_ < conv_thresh && check_sensi_conv_ && + sensitivityFlag == SensitivityMethod::forward) { + updateSensiSimulation(solver); + wrms_ = getWrmsFSA(model); + } } if (wrms_ < conv_thresh) @@ -676,12 +671,10 @@ void SteadystateProblem::runSteadystateSimulation(const Solver &solver, /* increase counter, check for maxsteps */ sim_steps++; if (sim_steps >= solver.getMaxSteps()) { - numsteps_.at(1) = sim_steps; throw NewtonFailure(AMICI_TOO_MUCH_WORK, "exceeded maximum number of steps"); } if (state_.t >= 1e200) { - numsteps_.at(1) = sim_steps; throw NewtonFailure(AMICI_NO_STEADY_STATE, "simulated to late time" " point without convergence of RHS"); @@ -691,13 +684,6 @@ void SteadystateProblem::runSteadystateSimulation(const Solver &solver, // if check_sensi_conv_ is deactivated, we still have to update sensis if (sensitivityFlag == SensitivityMethod::forward) updateSensiSimulation(solver); - - /* store information about steps and sensitivities, if necessary */ - if (backward) { - numstepsB_ = sim_steps; - } else { - numsteps_.at(1) = sim_steps; - } } std::unique_ptr From 3ab64c237636a73bbf4da5b683af6d081ec7cea8 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 21 Apr 2022 13:09:09 +0200 Subject: [PATCH 02/10] Fix CI failures: Invalid directory passed to lcov (#1785) Setuptools [recently changed (>=v62.1.0)](https://github.com/pypa/setuptools/commit/1c23f5e1e4b18b50081cbabb2dea22bf345f5894) the name of the temporary build directory and lcov fails trying to access a non-existent directory (`geninfo: ERROR: cannot read /home/runner/work/AMICI/AMICI/python/sdist/build/temp.linux-x86_64-3.8/amici/src!`). This fixes that. --- .github/workflows/test_python_cplusplus.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test_python_cplusplus.yml b/.github/workflows/test_python_cplusplus.yml index d7fbb7d5ae..e62beea5a8 100644 --- a/.github/workflows/test_python_cplusplus.yml +++ b/.github/workflows/test_python_cplusplus.yml @@ -117,7 +117,7 @@ jobs: -d ${AMICI_DIR}/build/CMakeFiles/amici.dir/src \ -b ${AMICI_DIR} -c -o coverage_cpp.info \ && lcov --compat-libtool --no-external \ - -d ${AMICI_DIR}/python/sdist/build/temp.linux-x86_64-$(python3 --version | sed -E 's/.*([0-9]+\.[0-9]+)\..*/\1/')/amici/src \ + -d ${AMICI_DIR}/python/sdist/build/temp.linux-x86_64-$(python -c "import sys; print(sys.implementation.cache_tag)")/amici/src \ -b ${AMICI_DIR}/python/sdist -c -o coverage_py.info \ && lcov -a coverage_cpp.info -a coverage_py.info -o coverage.info From b951f39f238c4baf59a73a41938a72e47a7251f6 Mon Sep 17 00:00:00 2001 From: Dilan Pathirana <59329744+dilpath@users.noreply.github.com> Date: Thu, 21 Apr 2022 19:39:03 +0200 Subject: [PATCH 03/10] Fix error if no measurements in PEtab problem; fix type handling in PEtab parameter mapping (#1783) * Fix logging error if no measurements in PEtab problem * Fix type handling in PEtab parameter mapping (support supplying simulation conditions as a dictionary) Co-authored-by: Daniel Weindl --- python/amici/petab_objective.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/python/amici/petab_objective.py b/python/amici/petab_objective.py index 1e6dfb0e7f..0cb8b4d9d2 100644 --- a/python/amici/petab_objective.py +++ b/python/amici/petab_objective.py @@ -154,8 +154,11 @@ def simulate_petab( # Log results sim_cond = petab_problem.get_simulation_conditions_from_measurement_df() for i, rdata in enumerate(rdatas): - logger.debug(f"Condition: {sim_cond.iloc[i, :].values}, status: " - f"{rdata['status']}, llh: {rdata['llh']}") + sim_cond_id = "N/A" if sim_cond.empty else sim_cond.iloc[i, :].values + logger.debug( + f"Condition: {sim_cond_id}, status: {rdata['status']}, " + f"llh: {rdata['llh']}" + ) return { LLH: llh, @@ -230,7 +233,7 @@ def create_parameterized_edatas( def create_parameter_mapping( petab_problem: petab.Problem, - simulation_conditions: Union[pd.DataFrame, Dict], + simulation_conditions: Union[pd.DataFrame, List[Dict]], scaled_parameters: bool, amici_model: AmiciModel, ) -> ParameterMapping: @@ -254,6 +257,8 @@ def create_parameter_mapping( if simulation_conditions is None: simulation_conditions = \ petab_problem.get_simulation_conditions_from_measurement_df() + if isinstance(simulation_conditions, list): + simulation_conditions = pd.DataFrame(data=simulation_conditions) # Because AMICI globalizes all local parameters during model import, # we need to do that here as well to prevent parameter mapping errors From fd0db583a003f439b4968d72e07c3b026fd5558d Mon Sep 17 00:00:00 2001 From: Dilan Pathirana <59329744+dilpath@users.noreply.github.com> Date: Thu, 21 Apr 2022 21:43:16 +0200 Subject: [PATCH 04/10] Fix substitution of expressions into event roots (partially reverts #1599) (#1784) --- python/amici/ode_export.py | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/python/amici/ode_export.py b/python/amici/ode_export.py index a1f40409b6..74a5d5abb7 100644 --- a/python/amici/ode_export.py +++ b/python/amici/ode_export.py @@ -1297,7 +1297,7 @@ def parse_events(self) -> None: expr.set_val(self._process_heavisides(expr.get_val(), roots)) # remove all possible Heavisides from roots, which may arise from - # the substitution of `'w'` in `_get_unique_root` + # the substitution of `'w'` in `_collect_heaviside_roots` for root in roots: root.set_val(self._process_heavisides(root.get_val(), roots)) @@ -2066,15 +2066,6 @@ def _get_unique_root( unique identifier for root, or ``None`` if the root is not time-dependent """ - # substitute 'w' expressions into root expressions now, to avoid - # rewriting '{model_name}_root.cpp' and '{model_name}_stau.cpp' headers - # to include 'w.h' - w_sorted = toposort_symbols(dict(zip( - [expr.get_id() for expr in self._expressions], - [expr.get_val() for expr in self._expressions], - ))) - root_found = root_found.subs(w_sorted) - if not self._expr_is_time_dependent(root_found): return None @@ -2115,6 +2106,18 @@ def _collect_heaviside_roots( elif arg.has(sp.Heaviside): root_funs.extend(self._collect_heaviside_roots(arg.args)) + # substitute 'w' expressions into root expressions now, to avoid + # rewriting '{model_name}_root.cpp' and '{model_name}_stau.cpp' headers + # to include 'w.h' + w_sorted = toposort_symbols(dict(zip( + [expr.get_id() for expr in self._expressions], + [expr.get_val() for expr in self._expressions], + ))) + root_funs = [ + r.subs(w_sorted) + for r in root_funs + ] + return root_funs def _process_heavisides( From bac5d025966e7949f088d6eac7a9244389b55b5f Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 28 Apr 2022 15:41:48 +0200 Subject: [PATCH 05/10] More informative NaN/Inf warnings (#1640) Show more informative warnings when encountering `NaN` or `inf` in model quantities (include IDs where available and provide 2D indices instead of flat indices). Closes #1375 --- include/amici/misc.h | 11 + include/amici/model.h | 75 +++++- include/amici/sundials_matrix_wrapper.h | 10 + src/amici.cpp | 6 +- src/misc.cpp | 5 + src/model.cpp | 297 ++++++++++++++++++++++-- src/model_dae.cpp | 2 +- src/model_ode.cpp | 2 +- src/solver_cvodes.cpp | 34 +-- src/solver_idas.cpp | 30 +-- src/sundials_matrix_wrapper.cpp | 40 ++++ swig/amici.i | 1 + tests/cpp/unittests/testMisc.cpp | 59 +++++ 13 files changed, 510 insertions(+), 62 deletions(-) diff --git a/include/amici/misc.h b/include/amici/misc.h index ee069b5224..329ca9a48d 100644 --- a/include/amici/misc.h +++ b/include/amici/misc.h @@ -187,6 +187,17 @@ class ContextManager{ ContextManager(ContextManager &&other) = delete; }; + +/** + * @brief Convert a flat index to a pair of row/column indices, + * assuming row-major order. + * @param flat_idx flat index + * @param num_cols number of columns of referred to matrix + * @return row index, column index + */ +auto unravel_index(size_t flat_idx, size_t num_cols) + -> std::pair; + } // namespace amici #endif // AMICI_MISC_H diff --git a/include/amici/model.h b/include/amici/model.h index fac088af2c..e15fa9e51b 100644 --- a/include/amici/model.h +++ b/include/amici/model.h @@ -34,6 +34,42 @@ void serialize(Archive &ar, amici::Model &m, unsigned int version); namespace amici { +/** + * @brief Describes the various model quantities + */ +enum class ModelQuantity { + J, + JB, + Jv, + JvB, + JDiag, + sx, + sy, + ssigmay, + xdot, + sxdot, + xBdot, + x0_rdata, + x0, + x_rdata, + x, + dwdw, + dwdx, + dwdp, + y, + dydp, + dydx, + w, + root, + qBdot, + qBdot_ss, + xBdot_ss, + JSparseB_ss, +}; + +extern const std::map model_quantity_to_str; + + /** * @brief The Model class represents an AMICI ODE/DAE model. * @@ -1201,18 +1237,43 @@ class Model : public AbstractModel, public ModelDimensions { */ void updateHeavisideB(const int *rootsfound); + + /** + * @brief Check if the given array has only finite elements. + * + * For (1D) spans. + * + * @param array + * @param model_quantity The model quantity `array` corresponds to + * @return + */ + int checkFinite(gsl::span array, + ModelQuantity model_quantity) const; + /** + * @brief Check if the given array has only finite elements. + * + * For flattened 2D arrays. + * + * @param array Flattened matrix + * @param model_quantity The model quantity `array` corresponds to + * @param num_cols Number of columns of the non-flattened matrix + * @return + */ + int checkFinite(gsl::span array, + ModelQuantity model_quantity, + size_t num_cols) const; + /** * @brief Check if the given array has only finite elements. * - * If not, try to give hints by which other fields this could be caused. + * For SUNMatrix. * - * @param array Array to check - * @param fun Name of the function that generated the values (for more - * informative messages). - * @return `amici::AMICI_RECOVERABLE_ERROR` if a NaN/Inf value was found, - * `amici::AMICI_SUCCESS` otherwise + * @param m Matrix to check + * @param model_quantity The model quantity `m` corresponds to + * @param t current timepoint + * @return */ - int checkFinite(gsl::span array, const char *fun) const; + int checkFinite(SUNMatrix m, ModelQuantity model_quantity, realtype t) const; /** * @brief Set whether the result of every call to `Model::f*` should be diff --git a/include/amici/sundials_matrix_wrapper.h b/include/amici/sundials_matrix_wrapper.h index b94667e39e..52d6c3377e 100644 --- a/include/amici/sundials_matrix_wrapper.h +++ b/include/amici/sundials_matrix_wrapper.h @@ -532,6 +532,16 @@ class SUNMatrixWrapper { bool ownmat = true; }; + +/** + * @brief Convert a flat index to a pair of row/column indices. + * @param i flat index + * @param m referred to matrix + * @return row index, column index + */ +auto unravel_index(sunindextype i, SUNMatrix m) + -> std::pair; + } // namespace amici namespace gsl { diff --git a/src/amici.cpp b/src/amici.cpp index 5be90dad00..cc20bd5927 100644 --- a/src/amici.cpp +++ b/src/amici.cpp @@ -179,7 +179,7 @@ AmiciApplication::runAmiciSimulation(Solver& solver, throw; warningF("AMICI:simulation", "AMICI forward simulation failed at t = %f: " - "Maximum time exceeed.\n", + "Maximum time exceeded.\n", ex.time); } else { rdata->status = ex.error_code; @@ -199,7 +199,7 @@ AmiciApplication::runAmiciSimulation(Solver& solver, warningF( "AMICI:simulation", "AMICI backward simulation failed when trying to solve until " - "t = %f: Maximum time exceeed.\n", + "t = %f: Maximum time exceeded.\n", ex.time); } else { @@ -305,6 +305,8 @@ AmiciApplication::errorF(const char* identifier, const char* format, ...) const int AmiciApplication::checkFinite(gsl::span array, const char* fun) { + Expects(array.size() + <= static_cast(std::numeric_limits::max())); for (int idx = 0; idx < (int)array.size(); idx++) { if (isNaN(array[idx])) { diff --git a/src/misc.cpp b/src/misc.cpp index 60e9d3ea18..9b4c650cba 100644 --- a/src/misc.cpp +++ b/src/misc.cpp @@ -178,4 +178,9 @@ std::string printfToString(const char *fmt, va_list ap) { return str; } +std::pair unravel_index(size_t flat_idx, size_t num_cols) +{ + return {flat_idx / num_cols, flat_idx % num_cols}; +} + } // namespace amici diff --git a/src/model.cpp b/src/model.cpp index 1edc211dde..fe6d138b70 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -15,6 +15,39 @@ namespace amici { +/** + * @brief Maps ModelQuantity items to their string value + */ +const std::map model_quantity_to_str { + {ModelQuantity::J, "J"}, + {ModelQuantity::JB, "JB"}, + {ModelQuantity::Jv, "Jv"}, + {ModelQuantity::JvB, "JvB"}, + {ModelQuantity::JDiag, "JDiag"}, + {ModelQuantity::sx, "sx"}, + {ModelQuantity::sy, "sy"}, + {ModelQuantity::ssigmay, "ssigmay"}, + {ModelQuantity::xdot, "xdot"}, + {ModelQuantity::sxdot, "sxdot"}, + {ModelQuantity::xBdot, "xBdot"}, + {ModelQuantity::x0, "x0"}, + {ModelQuantity::x0_rdata, "x0_rdata"}, + {ModelQuantity::x, "x"}, + {ModelQuantity::x_rdata, "x_rdata"}, + {ModelQuantity::dwdw, "dwdw"}, + {ModelQuantity::dwdx, "dwdx"}, + {ModelQuantity::dwdp, "dwdp"}, + {ModelQuantity::y, "y"}, + {ModelQuantity::dydp, "dydp"}, + {ModelQuantity::dydx, "dydx"}, + {ModelQuantity::w, "w"}, + {ModelQuantity::root, "root"}, + {ModelQuantity::qBdot, "qBdot"}, + {ModelQuantity::qBdot_ss, "qBdot_ss"}, + {ModelQuantity::xBdot_ss, "xBdot_ss"}, + {ModelQuantity::JSparseB_ss, "JSparseB_ss"}, +}; + /** * @brief local helper function to get parameters * @param ids vector of name/ids of (fixed)Parameters @@ -843,7 +876,7 @@ void Model::getObservableSensitivity(gsl::span sy, const realtype t, writeSlice(derived_state_.dydp_, sy); if (always_check_finite_) - checkFinite(sy, "sy"); + checkFinite(sy, ModelQuantity::sy, nplist()); } void Model::getObservableSigma(gsl::span sigmay, const int it, @@ -875,7 +908,7 @@ void Model::getObservableSigmaSensitivity(gsl::span ssigmay, } if (always_check_finite_) - checkFinite(ssigmay, "ssigmay"); + checkFinite(ssigmay, ModelQuantity::ssigmay, nplist()); } void Model::addObservableObjective(realtype &Jy, const int it, @@ -1218,18 +1251,242 @@ void Model::updateHeavisideB(const int *rootsfound) { } } +int Model::checkFinite(gsl::span array, + ModelQuantity model_quantity) const +{ + auto it = std::find_if(array.begin(), array.end(), + [](realtype x){return !std::isfinite(x);}); + if(it == array.end()) { + return AMICI_SUCCESS; + } + + // there is some issue - produce a meaningful message + auto flat_index = it - array.begin(); + + std::string msg_id; + std::string non_finite_type; + if (std::isnan(array[flat_index])) { + msg_id = "AMICI:NaN"; + non_finite_type = "NaN"; + } else if (std::isinf(array[flat_index])) { + msg_id = "AMICI:Inf"; + non_finite_type = "Inf"; + } + std::string element_id = std::to_string(flat_index); + + switch (model_quantity) { + case ModelQuantity::xdot: + case ModelQuantity::xBdot: + case ModelQuantity::x0: + case ModelQuantity::x: + case ModelQuantity::x_rdata: + case ModelQuantity::x0_rdata: + case ModelQuantity::Jv: + case ModelQuantity::JvB: + case ModelQuantity::JDiag: + if(hasStateIds()) { + element_id = getStateIdsSolver()[flat_index]; + } + break; + case ModelQuantity::y: + if(hasObservableIds()) { + element_id = getObservableIds()[flat_index]; + } + break; + case ModelQuantity::w: + if(hasExpressionIds()) { + element_id = getExpressionIds()[flat_index]; + } + break; + default: + break; + } -int Model::checkFinite(gsl::span array, const char *fun) const { - auto result = app->checkFinite(array, fun); - - if (result != AMICI_SUCCESS) { - app->checkFinite(state_.fixedParameters, "k"); - app->checkFinite(state_.unscaledParameters, "p"); + std::string model_quantity_str; + try { + model_quantity_str = model_quantity_to_str.at(model_quantity); + } catch (std::out_of_range const&) { + // Missing model quantity string - terminate if this is a debug build, + // but show the quantity number if non-debug. + gsl_ExpectsDebug(false); + model_quantity_str = std::to_string(static_cast(model_quantity)); + } + app->warningF(msg_id.c_str(), + "AMICI encountered a %s value for %s[%i] (%s)", + non_finite_type.c_str(), + model_quantity_str.c_str(), + static_cast(flat_index), + element_id.c_str() + ); + + // check upstream + app->checkFinite(state_.fixedParameters, "k"); + app->checkFinite(state_.unscaledParameters, "p"); + app->checkFinite(simulation_parameters_.ts_, "t"); + if(!always_check_finite_) { + // don't check twice if always_check_finite_ is true app->checkFinite(derived_state_.w_, "w"); - app->checkFinite(simulation_parameters_.ts_, "t"); } - return result; + return AMICI_RECOVERABLE_ERROR; +} + +int Model::checkFinite(gsl::span array, + ModelQuantity model_quantity, size_t num_cols) const +{ + auto it = std::find_if(array.begin(), array.end(), + [](realtype x){return !std::isfinite(x);}); + if(it == array.end()) { + return AMICI_SUCCESS; + } + + // there is some issue - produce a meaningful message + auto flat_index = it - array.begin(); + sunindextype row, col; + std::tie(row, col) = unravel_index(flat_index, num_cols); + std::string msg_id; + std::string non_finite_type; + if (std::isnan(array[flat_index])) { + msg_id = "AMICI:NaN"; + non_finite_type = "NaN"; + } else if (std::isinf(array[flat_index])) { + msg_id = "AMICI:Inf"; + non_finite_type = "Inf"; + } + std::string row_id = std::to_string(row); + std::string col_id = std::to_string(col); + + switch (model_quantity) { + case ModelQuantity::sy: + case ModelQuantity::ssigmay: + case ModelQuantity::dydp: + row_id += " " + getObservableIds()[row]; + col_id += " " + getParameterIds()[plist(col)]; + break; + case ModelQuantity::dydx: + row_id += " " + getObservableIds()[row]; + col_id += " " + getStateIdsSolver()[col]; + break; + default: + break; + } + + std::string model_quantity_str; + try { + model_quantity_str = model_quantity_to_str.at(model_quantity); + } catch (std::out_of_range const&) { + // Missing model quantity string - terminate if this is a debug build, + // but show the quantity number if non-debug. + gsl_ExpectsDebug(false); + model_quantity_str = std::to_string(static_cast(model_quantity)); + } + + app->warningF(msg_id.c_str(), + "AMICI encountered a %s value for %s[%i] (%s, %s)", + non_finite_type.c_str(), + model_quantity_str.c_str(), + static_cast(flat_index), + row_id.c_str(), + col_id.c_str() + ); + + // check upstream + app->checkFinite(state_.fixedParameters, "k"); + app->checkFinite(state_.unscaledParameters, "p"); + app->checkFinite(simulation_parameters_.ts_, "t"); + app->checkFinite(derived_state_.w_, "w"); + + return AMICI_RECOVERABLE_ERROR; +} + +int Model::checkFinite(SUNMatrix m, ModelQuantity model_quantity, realtype t) const +{ + // check flat array, to see if there are any issues + // (faster, in particular for sparse arrays) + auto m_flat = gsl::make_span(m); + auto it = std::find_if(m_flat.begin(), m_flat.end(), + [](realtype x){return !std::isfinite(x);}); + if(it == m_flat.end()) { + return AMICI_SUCCESS; + } + + // there is some issue - produce a meaningful message + auto flat_index = it - m_flat.begin(); + sunindextype row, col; + std::tie(row, col) = unravel_index(flat_index, m); + std::string msg_id; + std::string non_finite_type; + if (std::isnan(m_flat[flat_index])) { + msg_id = "AMICI:NaN"; + non_finite_type = "NaN"; + } else if (std::isinf(m_flat[flat_index])) { + msg_id = "AMICI:Inf"; + non_finite_type = "Inf"; + } else { + throw std::runtime_error( + "Value is not finite, but neither infinite nor NaN."); + } + std::string row_id = std::to_string(row); + std::string col_id = std::to_string(col); + + switch (model_quantity) { + case ModelQuantity::J: + case ModelQuantity::JB: + if(hasStateIds()) { + auto state_ids = getStateIdsSolver(); + row_id += " " + state_ids[row]; + col_id += " " + state_ids[col]; + } + break; + case ModelQuantity::dwdx: + if(hasExpressionIds()) + row_id += " " + getExpressionIds()[row]; + if(hasStateIds()) + col_id += " " + getStateIdsSolver()[col]; + break; + case ModelQuantity::dwdw: + if(hasExpressionIds()) { + auto expr_ids = getExpressionIds(); + row_id += " " + expr_ids[row]; + col_id += " " + expr_ids[col]; + } + break; + case ModelQuantity::dwdp: + if(hasExpressionIds()) + row_id += " " + getExpressionIds()[row]; + if(hasParameterIds()) + col_id += " " + getParameterIds()[plist(col)]; + break; + default: + break; + } + + std::string model_quantity_str; + try { + model_quantity_str = model_quantity_to_str.at(model_quantity); + } catch (std::out_of_range const&) { + // Missing model quantity string - terminate if this is a debug build, + // but show the quantity number if non-debug. + gsl_ExpectsDebug(false); + model_quantity_str = std::to_string(static_cast(model_quantity)); + } + app->warningF(msg_id.c_str(), + "AMICI encountered a %s value for %s[%i] (%s, %s) at t=%g", + non_finite_type.c_str(), + model_quantity_str.c_str(), + static_cast(flat_index), + row_id.c_str(), + col_id.c_str(), + t + ); + + // check upstream + app->checkFinite(state_.fixedParameters, "k"); + app->checkFinite(state_.unscaledParameters, "p"); + app->checkFinite(simulation_parameters_.ts_, "t"); + app->checkFinite(derived_state_.w_, "w"); + + return AMICI_RECOVERABLE_ERROR; } void Model::setAlwaysCheckFinite(bool alwaysCheck) { @@ -1250,8 +1507,8 @@ void Model::fx0(AmiVector &x) { state_.fixedParameters.data()); if (always_check_finite_) { - checkFinite(derived_state_.x_rdata_, "x0 x_rdata"); - checkFinite(x.getVector(), "x0 x"); + checkFinite(derived_state_.x_rdata_, ModelQuantity::x0_rdata); + checkFinite(x.getVector(), ModelQuantity::x0); } } @@ -1333,7 +1590,7 @@ void Model::fx_rdata(AmiVector &x_rdata, const AmiVector &x) { fx_rdata(x_rdata.data(), computeX_pos(x), state_.total_cl.data(), state_.unscaledParameters.data(), state_.fixedParameters.data()); if (always_check_finite_) - checkFinite(x_rdata.getVector(), "x_rdata"); + checkFinite(x_rdata.getVector(), ModelQuantity::x_rdata); } void Model::fsx_rdata(AmiVectorArray &sx_rdata, const AmiVectorArray &sx, @@ -1411,7 +1668,7 @@ void Model::fy(const realtype t, const AmiVector &x) { state_.h.data(), derived_state_.w_.data()); if (always_check_finite_) { - app->checkFinite(gsl::make_span(derived_state_.y_.data(), ny), "y"); + checkFinite(gsl::make_span(derived_state_.y_.data(), ny), ModelQuantity::y); } } @@ -1441,7 +1698,7 @@ void Model::fdydp(const realtype t, const AmiVector &x) { } if (always_check_finite_) { - app->checkFinite(derived_state_.dydp_, "dydp"); + checkFinite(derived_state_.dydp_, ModelQuantity::dydp, nplist()); } } @@ -1461,7 +1718,7 @@ void Model::fdydx(const realtype t, const AmiVector &x) { derived_state_.dwdx_.data()); if (always_check_finite_) { - app->checkFinite(derived_state_.dydx_, "dydx"); + checkFinite(derived_state_.dydx_, ModelQuantity::dydx, ny); } } @@ -2076,7 +2333,7 @@ void Model::fw(const realtype t, const realtype *x) { state_.fixedParameters.data(), state_.h.data(), state_.total_cl.data()); if (always_check_finite_) { - app->checkFinite(derived_state_.w_, "w"); + checkFinite(derived_state_.w_, ModelQuantity::w); } } @@ -2116,7 +2373,7 @@ void Model::fdwdp(const realtype t, const realtype *x) { } if (always_check_finite_) { - app->checkFinite(gsl::make_span(derived_state_.dwdp_.get()), "dwdp"); + checkFinite(derived_state_.dwdp_.get(), ModelQuantity::dwdp, t); } } @@ -2157,7 +2414,7 @@ void Model::fdwdx(const realtype t, const realtype *x) { } if (always_check_finite_) { - app->checkFinite(gsl::make_span(derived_state_.dwdx_.get()), "dwdx"); + checkFinite(derived_state_.dwdx_.get(), ModelQuantity::dwdx, t); } } @@ -2172,7 +2429,7 @@ void Model::fdwdw(const realtype t, const realtype *x) { derived_state_.w_.data(), state_.total_cl.data()); if (always_check_finite_) { - app->checkFinite(gsl::make_span(dwdw_.get()), "dwdw"); + checkFinite(dwdw_.get(), ModelQuantity::dwdw, t); } } diff --git a/src/model_dae.cpp b/src/model_dae.cpp index bed9b4c8c9..2e3d6d06ee 100644 --- a/src/model_dae.cpp +++ b/src/model_dae.cpp @@ -87,7 +87,7 @@ void Model_DAE::fJDiag(const realtype t, AmiVector &JDiag, fJSparse(t, 0.0, x.getNVector(), dx.getNVector(), derived_state_.J_.get()); derived_state_.J_.refresh(); derived_state_.J_.to_diag(JDiag.getNVector()); - if (!checkFinite(JDiag.getVector(), "Jacobian")) + if (checkFinite(JDiag.getVector(), ModelQuantity::JDiag) != AMICI_SUCCESS) throw AmiException("Evaluation of fJDiag failed!"); } diff --git a/src/model_ode.cpp b/src/model_ode.cpp index 02905ef2c0..65552b8579 100644 --- a/src/model_ode.cpp +++ b/src/model_ode.cpp @@ -105,7 +105,7 @@ void Model_ODE::fJDiag(const realtype t, AmiVector &JDiag, const realtype /*cj*/, const AmiVector &x, const AmiVector & /*dx*/) { fJDiag(t, JDiag.getNVector(), x.getNVector()); - if (checkFinite(JDiag.getVector(), "Jacobian") != AMICI_SUCCESS) + if (checkFinite(JDiag.getVector(), ModelQuantity::JDiag) != AMICI_SUCCESS) throw AmiException("Evaluation of fJDiag failed!"); } diff --git a/src/solver_cvodes.cpp b/src/solver_cvodes.cpp index c37b1207cd..74998fdf66 100644 --- a/src/solver_cvodes.cpp +++ b/src/solver_cvodes.cpp @@ -807,7 +807,6 @@ const Model *CVodeSolver::getModel() const { /** * @brief Jacobian of xdot with respect to states x - * @param N number of state variables * @param t timepoint * @param x Vector with the states * @param xdot Vector with the right hand side @@ -827,13 +826,12 @@ static int fJ(realtype t, N_Vector x, N_Vector xdot, SUNMatrix J, Expects(model); model->fJ(t, x, xdot, J); - return model->checkFinite(gsl::make_span(J), "Jacobian"); + return model->checkFinite(J, ModelQuantity::J, t); } /** * @brief Jacobian of xBdot with respect to adjoint state xB - * @param NeqBdot number of adjoint state variables * @param t timepoint * @param x Vector with the states * @param xB Vector with the adjoint states @@ -854,7 +852,7 @@ static int fJB(realtype t, N_Vector x, N_Vector xB, N_Vector xBdot, Expects(model); model->fJB(t, x, xB, xBdot, JB); - return model->checkFinite(gsl::make_span(JB), "Jacobian"); + return model->checkFinite(gsl::make_span(JB), ModelQuantity::JB); } @@ -879,7 +877,7 @@ static int fJSparse(realtype t, N_Vector x, N_Vector /*xdot*/, Expects(model); model->fJSparse(t, x, J); - return model->checkFinite(gsl::make_span(J), "Jacobian"); + return model->checkFinite(J, ModelQuantity::J, t); } @@ -905,7 +903,7 @@ static int fJSparseB(realtype t, N_Vector x, N_Vector xB, N_Vector xBdot, Expects(model); model->fJSparseB(t, x, xB, xBdot, JB); - return model->checkFinite(gsl::make_span(JB), "Jacobian"); + return model->checkFinite(gsl::make_span(JB), ModelQuantity::JB); } @@ -967,7 +965,7 @@ static int fJv(N_Vector v, N_Vector Jv, realtype t, N_Vector x, Expects(model); model->fJv(v, Jv, t, x); - return model->checkFinite(gsl::make_span(Jv), "Jacobian"); + return model->checkFinite(gsl::make_span(Jv), ModelQuantity::Jv); } @@ -993,7 +991,7 @@ static int fJvB(N_Vector vB, N_Vector JvB, realtype t, N_Vector x, Expects(model); model->fJvB(vB, JvB, t, x, xB); - return model->checkFinite(gsl::make_span(JvB), "Jacobian"); + return model->checkFinite(gsl::make_span(JvB), ModelQuantity::JvB); } @@ -1014,7 +1012,7 @@ static int froot(realtype t, N_Vector x, realtype *root, model->froot(t, x, gsl::make_span(root, model->ne)); return model->checkFinite(gsl::make_span(root, model->ne), - "root function"); + ModelQuantity::root); } @@ -1038,7 +1036,8 @@ static int fxdot(realtype t, N_Vector x, N_Vector xdot, void *user_data) { return AMICI_MAX_TIME_EXCEEDED; } - if (t > 1e200 && !model->checkFinite(gsl::make_span(x), "fxdot")) { + if (t > 1e200 + && !model->checkFinite(gsl::make_span(x), ModelQuantity::xdot)) { /* when t is large (typically ~1e300), CVODES may pass all NaN x to fxdot from which we typically cannot recover. To save time on normal execution, we do not always want to check finiteness @@ -1047,7 +1046,7 @@ static int fxdot(realtype t, N_Vector x, N_Vector xdot, void *user_data) { } model->fxdot(t, x, xdot); - return model->checkFinite(gsl::make_span(xdot), "fxdot"); + return model->checkFinite(gsl::make_span(xdot), ModelQuantity::xdot); } @@ -1074,7 +1073,7 @@ static int fxBdot(realtype t, N_Vector x, N_Vector xB, N_Vector xBdot, } model->fxBdot(t, x, xB, xBdot); - return model->checkFinite(gsl::make_span(xBdot), "fxBdot"); + return model->checkFinite(gsl::make_span(xBdot), ModelQuantity::xBdot); } @@ -1095,7 +1094,7 @@ static int fqBdot(realtype t, N_Vector x, N_Vector xB, N_Vector qBdot, Expects(model); model->fqBdot(t, x, xB, qBdot); - return model->checkFinite(gsl::make_span(qBdot), "qBdot"); + return model->checkFinite(gsl::make_span(qBdot), ModelQuantity::qBdot); } @@ -1116,7 +1115,7 @@ static int fxBdot_ss(realtype t, N_Vector xB, N_Vector xBdot, Expects(model); model->fxBdot_ss(t, xB, xBdot); - return model->checkFinite(gsl::make_span(xBdot), "fxBdot_ss"); + return model->checkFinite(gsl::make_span(xBdot), ModelQuantity::xBdot_ss); } @@ -1137,7 +1136,7 @@ static int fqBdot_ss(realtype t, N_Vector xB, N_Vector qBdot, Expects(model); model->fqBdot_ss(t, xB, qBdot); - return model->checkFinite(gsl::make_span(qBdot), "qBdot_ss"); + return model->checkFinite(gsl::make_span(qBdot), ModelQuantity::qBdot_ss); } /** @@ -1161,7 +1160,8 @@ static int fJSparseB_ss(realtype /*t*/, N_Vector /*x*/, N_Vector xBdot, Expects(model); model->fJSparseB_ss(JB); - return model->checkFinite(gsl::make_span(xBdot), "JSparseB_ss"); + return model->checkFinite(gsl::make_span(xBdot), + ModelQuantity::JSparseB_ss); } @@ -1189,7 +1189,7 @@ static int fsxdot(int /*Ns*/, realtype t, N_Vector x, N_Vector /*xdot*/, Expects(model); model->fsxdot(t, x, ip, sx, sxdot); - return model->checkFinite(gsl::make_span(sxdot), "sxdot"); + return model->checkFinite(gsl::make_span(sxdot), ModelQuantity::sxdot); } bool operator==(const CVodeSolver &a, const CVodeSolver &b) { diff --git a/src/solver_idas.cpp b/src/solver_idas.cpp index 0742f9c095..0039f6de0a 100644 --- a/src/solver_idas.cpp +++ b/src/solver_idas.cpp @@ -817,7 +817,7 @@ int fJ(realtype t, realtype cj, N_Vector x, N_Vector dx, 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"); + return model->checkFinite(J, ModelQuantity::J, t); } /** @@ -846,7 +846,8 @@ int fJB(realtype t, realtype cj, N_Vector x, N_Vector dx, Expects(model); model->fJB(t, cj, x, dx, xB, dxB, JB); - return model->checkFinite(gsl::make_span(JB), "Jacobian"); + return model->checkFinite(JB, ModelQuantity::JB, t); + } /** @@ -873,7 +874,7 @@ int fJSparse(realtype t, realtype cj, N_Vector x, N_Vector dx, Expects(model); model->fJSparse(t, cj, x, dx, J); - return model->checkFinite(gsl::make_span(J), "Jacobian"); + return model->checkFinite(J, ModelQuantity::J, t); } /** @@ -902,7 +903,7 @@ int fJSparseB(realtype t, realtype cj, N_Vector x, N_Vector dx, Expects(model); model->fJSparseB(t, cj, x, dx, xB, dxB, JB); - return model->checkFinite(gsl::make_span(JB), "Jacobian"); + return model->checkFinite(JB, ModelQuantity::JB, t); } /** @@ -973,7 +974,7 @@ int fJv(realtype t, N_Vector x, N_Vector dx, N_Vector /*xdot*/, Expects(model); model->fJv(t, x, dx, v, Jv, cj); - return model->checkFinite(gsl::make_span(Jv), "Jacobian"); + return model->checkFinite(gsl::make_span(Jv), ModelQuantity::Jv); } /** @@ -1003,7 +1004,7 @@ int fJvB(realtype t, N_Vector x, N_Vector dx, N_Vector xB, Expects(model); model->fJvB(t, x, dx, xB, dxB, vB, JvB, cj); - return model->checkFinite(gsl::make_span(JvB), "Jacobian"); + return model->checkFinite(gsl::make_span(JvB), ModelQuantity::JvB); } /** @@ -1024,7 +1025,7 @@ int froot(realtype t, N_Vector x, N_Vector dx, realtype *root, model->froot(t, x, dx, gsl::make_span(root, model->ne)); return model->checkFinite(gsl::make_span(root, model->ne), - "root function"); + ModelQuantity::root); } /** @@ -1058,7 +1059,7 @@ int fxdot(realtype t, N_Vector x, N_Vector dx, N_Vector xdot, } model->fxdot(t, x, dx, xdot); - return model->checkFinite(gsl::make_span(xdot), "fxdot"); + return model->checkFinite(gsl::make_span(xdot), ModelQuantity::xdot); } /** @@ -1086,7 +1087,7 @@ int fxBdot(realtype t, N_Vector x, N_Vector dx, N_Vector xB, } model->fxBdot(t, x, dx, xB, dxB, xBdot); - return model->checkFinite(gsl::make_span(xBdot), "xBdot"); + return model->checkFinite(gsl::make_span(xBdot), ModelQuantity::xBdot); } /** @@ -1109,7 +1110,7 @@ int fqBdot(realtype t, N_Vector x, N_Vector dx, N_Vector xB, Expects(model); model->fqBdot(t, x, dx, xB, dxB, qBdot); - return model->checkFinite(gsl::make_span(qBdot), "qBdot"); + return model->checkFinite(gsl::make_span(qBdot), ModelQuantity::qBdot); } @@ -1132,7 +1133,7 @@ static int fxBdot_ss(realtype t, N_Vector xB, N_Vector dxB, N_Vector xBdot, Expects(model); model->fxBdot_ss(t, xB, dxB, xBdot); - return model->checkFinite(gsl::make_span(xBdot), "xBdot_ss"); + return model->checkFinite(gsl::make_span(xBdot), ModelQuantity::xBdot_ss); } @@ -1154,7 +1155,7 @@ static int fqBdot_ss(realtype t, N_Vector xB, N_Vector dxB, N_Vector qBdot, Expects(model); model->fqBdot_ss(t, xB, dxB, qBdot); - return model->checkFinite(gsl::make_span(qBdot), "qBdot_ss"); + return model->checkFinite(gsl::make_span(qBdot), ModelQuantity::qBdot_ss); } /** @@ -1181,7 +1182,8 @@ static int fqBdot_ss(realtype t, N_Vector xB, N_Vector dxB, N_Vector qBdot, Expects(model); model->fJSparseB_ss(JB); - return model->checkFinite(gsl::make_span(xBdot), "JSparseB_ss"); + return model->checkFinite(gsl::make_span(xBdot), + ModelQuantity::JSparseB_ss); } /** @@ -1212,7 +1214,7 @@ int fsxdot(int /*Ns*/, realtype t, N_Vector x, N_Vector dx, for (int ip = 0; ip < model->nplist(); ip++) { model->fsxdot(t, x, dx, ip, sx[ip], sdx[ip], sxdot[ip]); - if (model->checkFinite(gsl::make_span(sxdot[ip]), "sxdot") + if (model->checkFinite(gsl::make_span(sxdot[ip]), ModelQuantity::sxdot) != AMICI_SUCCESS) return AMICI_RECOVERABLE_ERROR; } diff --git a/src/sundials_matrix_wrapper.cpp b/src/sundials_matrix_wrapper.cpp index bd3fda4fb8..f3515c357a 100644 --- a/src/sundials_matrix_wrapper.cpp +++ b/src/sundials_matrix_wrapper.cpp @@ -747,5 +747,45 @@ void SUNMatrixWrapper::refresh() { SUNMatrix SUNMatrixWrapper::get() const { return matrix_; } +std::pair unravel_index(sunindextype i, SUNMatrix m) +{ + gsl_ExpectsDebug(i >= 0); + auto mat_id = SUNMatGetID(m); + if(mat_id == SUNMATRIX_DENSE) { + gsl_ExpectsDebug(i < SM_COLUMNS_D(m) * SM_ROWS_D(m)); + + auto num_rows = SM_ROWS_D(m); + // col-major + sunindextype row = i % num_rows; + sunindextype col = i / num_rows; + + gsl_EnsuresDebug(row >= 0); + gsl_EnsuresDebug(row < SM_ROWS_D(m)); + gsl_EnsuresDebug(col >= 0); + gsl_EnsuresDebug(col < SM_COLUMNS_D(m)); + + return {row, col}; + } + + if(mat_id == SUNMATRIX_SPARSE) { + gsl_ExpectsDebug(i < SM_NNZ_S(m)); + sunindextype row = SM_INDEXVALS_S(m)[i]; + sunindextype i_colptr = 0; + while(SM_INDEXPTRS_S(m)[i_colptr] < SM_NNZ_S(m)) { + if(SM_INDEXPTRS_S(m)[i_colptr + 1] > i) { + sunindextype col = i_colptr; + gsl_EnsuresDebug(row >= 0); + gsl_EnsuresDebug(row < SM_ROWS_S(m)); + gsl_EnsuresDebug(col >= 0); + gsl_EnsuresDebug(col < SM_COLUMNS_S(m)); + return {row, col}; + } + ++i_colptr; + } + } + + throw amici::AmiException("Unimplemented SUNMatrix type for unravel_index"); +} + } // namespace amici diff --git a/swig/amici.i b/swig/amici.i index a754cd61b8..ccfbf65933 100644 --- a/swig/amici.i +++ b/swig/amici.i @@ -120,6 +120,7 @@ wrap_unique_ptr(ExpDataPtr, amici::ExpData) %ignore amici::ContextManager; %ignore amici::ModelState; %ignore amici::ModelStateDerived; +%ignore amici::unravel_index; // Include before any other header which uses enums defined there %include "amici/defines.h" diff --git a/tests/cpp/unittests/testMisc.cpp b/tests/cpp/unittests/testMisc.cpp index 278e93b2bc..c689090c8c 100644 --- a/tests/cpp/unittests/testMisc.cpp +++ b/tests/cpp/unittests/testMisc.cpp @@ -6,6 +6,7 @@ #include #include #include +#include #include #include @@ -646,4 +647,62 @@ TEST_F(SunMatrixWrapperTest, BlockTranspose) SM_INDEXPTRS_S(B_sparse.get())[icol]); } +TEST(UnravelIndex, UnravelIndex) +{ + EXPECT_EQ(unravel_index(0, 2), std::make_pair((size_t) 0, (size_t) 0)); + EXPECT_EQ(unravel_index(1, 2), std::make_pair((size_t) 0, (size_t) 1)); + EXPECT_EQ(unravel_index(2, 2), std::make_pair((size_t) 1, (size_t) 0)); + EXPECT_EQ(unravel_index(3, 2), std::make_pair((size_t) 1, (size_t) 1)); + EXPECT_EQ(unravel_index(4, 2), std::make_pair((size_t) 2, (size_t) 0)); + EXPECT_EQ(unravel_index(5, 2), std::make_pair((size_t) 2, (size_t) 1)); + EXPECT_EQ(unravel_index(6, 2), std::make_pair((size_t) 3, (size_t) 0)); +} + +TEST(UnravelIndex, UnravelIndexSunMatDense) +{ + SUNMatrixWrapper A = SUNMatrixWrapper(3, 2); + + A.set_data(0, 0, 0); + A.set_data(1, 0, 1); + A.set_data(2, 0, 2); + A.set_data(0, 1, 3); + A.set_data(1, 1, 4); + A.set_data(2, 1, 5); + + for(int i = 0; i < 6; ++i) { + auto idx = unravel_index(i, A.get()); + EXPECT_EQ(A.get_data(idx.first, idx.second), i); + } +} + +TEST(UnravelIndex, UnravelIndexSunMatSparse) +{ + SUNMatrixWrapper D = SUNMatrixWrapper(4, 2); + + // [0, 3] + // [0, 0] + // [1, 0] + // [2, 0] + // data [1, 2, 3] + // colptrs [0, 2, 3] + // rowidxs [2, 3, 1] + D.set_data(0, 0, 0); + D.set_data(1, 0, 0); + D.set_data(2, 0, 1); + D.set_data(3, 0, 2); + D.set_data(0, 1, 3); + D.set_data(1, 1, 0); + D.set_data(2, 1, 0); + D.set_data(3, 1, 0); + + auto S = SUNSparseFromDenseMatrix(D.get(), 1e-15, CSC_MAT); + + EXPECT_EQ(unravel_index(0, S), std::make_pair((sunindextype) 2, (sunindextype) 0)); + EXPECT_EQ(unravel_index(1, S), std::make_pair((sunindextype) 3, (sunindextype) 0)); + EXPECT_EQ(unravel_index(2, S), std::make_pair((sunindextype) 0, (sunindextype) 1)); + + SUNMatDestroy(S); +} + + } // namespace From a6f9894d29d2453af67a8817cd7bb8c398a7b8d8 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 3 May 2022 14:46:13 +0200 Subject: [PATCH 06/10] Handle function 'add_global_parameter' removed from PEtab (#1790) Removed in https://github.com/PEtab-dev/libpetab-python/pull/140/ --- python/amici/petab_import.py | 34 +++++++++++++++++++++++++++++++++- 1 file changed, 33 insertions(+), 1 deletion(-) diff --git a/python/amici/petab_import.py b/python/amici/petab_import.py index f1183b451f..e3c8495f05 100644 --- a/python/amici/petab_import.py +++ b/python/amici/petab_import.py @@ -40,6 +40,38 @@ PREEQ_INDICATOR_ID = 'preequilibration_indicator' +def _add_global_parameter(sbml_model: libsbml.Model, + parameter_id: str, + parameter_name: str = None, + constant: bool = False, + units: str = 'dimensionless', + value: float = 0.0) -> libsbml.Parameter: + """Add new global parameter to SBML model + + Arguments: + sbml_model: SBML model + parameter_id: ID of the new parameter + parameter_name: Name of the new parameter + constant: Is parameter constant? + units: SBML unit ID + value: parameter value + + Returns: + The created parameter + """ + + if parameter_name is None: + parameter_name = parameter_id + + p = sbml_model.createParameter() + p.setId(parameter_id) + p.setName(parameter_name) + p.setConstant(constant) + p.setValue(value) + p.setUnits(units) + return p + + def get_fixed_parameters( sbml_model: 'libsbml.Model', condition_df: Optional[pd.DataFrame] = None, @@ -505,7 +537,7 @@ def import_model_sbml( output_parameters[sym] = None logger.debug(f"Adding output parameters to model: {output_parameters}") for par in output_parameters.keys(): - petab.add_global_parameter(sbml_model, par) + _add_global_parameter(sbml_model, par) # # TODO: to parameterize initial states or compartment sizes, we currently From 32c6c42f1b9498702b0fe679fff6c4288577497f Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 4 May 2022 10:26:09 +0200 Subject: [PATCH 07/10] Workaround for PEtab problems with state-dependent noise models (#1791) PEtab does currently not allow observables in noiseFormula and AMICI cannot handle states in sigma expressions. Therefore, where possible, replace species occurring in the error model definition by observableIds. See also #1788 --- python/amici/petab_import.py | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/python/amici/petab_import.py b/python/amici/petab_import.py index e3c8495f05..bc520f9d20 100644 --- a/python/amici/petab_import.py +++ b/python/amici/petab_import.py @@ -533,7 +533,8 @@ def import_model_sbml( key=lambda symbol: symbol.name) for free_sym in free_syms: sym = str(free_sym) - if sbml_model.getElementBySId(sym) is None and sym != 'time': + if sbml_model.getElementBySId(sym) is None and sym != 'time' \ + and sym not in observables: output_parameters[sym] = None logger.debug(f"Adding output parameters to model: {output_parameters}") for par in output_parameters.keys(): @@ -656,10 +657,15 @@ def get_observation_model( observables[oid] = {'name': name, 'formula': formula_obs} sigmas[oid] = formula_noise - # Replace observableIds occurring in error model definition + # PEtab does currently not allow observables in noiseFormula and AMICI + # cannot handle states in sigma expressions. Therefore, where possible, + # replace species occurring in error model definition by observableIds. + replacements = { + sp.sympify(observable['formula']): sp.Symbol(observable_id) + for observable_id, observable in observables.items() + } for observable_id, formula in sigmas.items(): - repl = sp.sympify(formula).subs( - observable_id, observables[observable_id]['formula']) + repl = sp.sympify(formula).subs(replacements) sigmas[observable_id] = str(repl) noise_distrs = petab_noise_distributions_to_amici(observable_df) From 82d0b5fecba11530ab954f0807a50a997b058c96 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20Fr=C3=B6hlich?= Date: Wed, 4 May 2022 10:24:15 -0400 Subject: [PATCH 08/10] fix initial events, fixes #1760 (#1789) * initial implementation, fixes #1760 * fixup * fixup * fixup * fix doc * fix doc * fixup swig * fixup sbml import * add sensitivity check, fixup time-independent triggers * fixup * fixup * fixup * fixup parameter event assignments * update test stats * Apply suggestions from code review Co-authored-by: Daniel Weindl Co-authored-by: Daniel Weindl --- documentation/python_interface.rst | 2 +- include/amici/forwardproblem.h | 4 +- include/amici/model.h | 16 +++++-- python/amici/ode_export.py | 50 +++++++++++++-------- python/amici/ode_model.py | 72 +++++++++++------------------- python/amici/sbml_import.py | 22 ++++----- src/forwardproblem.cpp | 48 +++++++++++++------- src/model.ODE_template.cpp | 4 ++ src/model.cpp | 16 +++++-- src/model_header.ODE_template.h | 7 ++- src/steadystateproblem.cpp | 6 ++- swig/model.i | 2 +- tests/testSBMLSuite.py | 2 + 13 files changed, 148 insertions(+), 103 deletions(-) diff --git a/documentation/python_interface.rst b/documentation/python_interface.rst index b4137c27e3..02cc22af35 100644 --- a/documentation/python_interface.rst +++ b/documentation/python_interface.rst @@ -26,7 +26,7 @@ AMICI can import :term:`SBML` models via the Status of SBML support in Python-AMICI ++++++++++++++++++++++++++++++++++++++ -Python-AMICI currently **passes 1014 out of the 1821 (~56%) test cases** from +Python-AMICI currently **passes 1030 out of the 1821 (~57%) test cases** from the semantic `SBML Test Suite `_ (`current status `_). diff --git a/include/amici/forwardproblem.h b/include/amici/forwardproblem.h index 712377d959..38c3ba1785 100644 --- a/include/amici/forwardproblem.h +++ b/include/amici/forwardproblem.h @@ -290,9 +290,11 @@ class ForwardProblem { * * @param tlastroot pointer to the timepoint of the last event * @param seflag Secondary event flag + * @param initial_event initial event flag */ - void handleEvent(realtype *tlastroot,bool seflag); + void handleEvent(realtype *tlastroot, bool seflag, + bool initial_event); /** * @brief Extract output information for events diff --git a/include/amici/model.h b/include/amici/model.h index e15fa9e51b..1897f9b295 100644 --- a/include/amici/model.h +++ b/include/amici/model.h @@ -97,7 +97,8 @@ class Model : public AbstractModel, public ModelDimensions { SimulationParameters simulation_parameters, amici::SecondOrderMode o2mode, std::vector idlist, - std::vector z2event, bool pythonGenerated = false, + std::vector z2event, + bool pythonGenerated = false, int ndxdotdp_explicit = 0, int ndxdotdx_explicit = 0, int w_recursion_depth = 0); @@ -202,9 +203,11 @@ class Model : public AbstractModel, public ModelDimensions { * @param sdx Reference to time derivative of state sensitivities (DAE only) * @param computeSensitivities Flag indicating whether sensitivities are to * be computed + * @param roots_found boolean indicators indicating whether roots were found at t0 by this fun */ void initialize(AmiVector &x, AmiVector &dx, AmiVectorArray &sx, - AmiVectorArray &sdx, bool computeSensitivities); + AmiVectorArray &sdx, bool computeSensitivities, + std::vector &roots_found); /** * @brief Initialize model properties. @@ -236,8 +239,10 @@ class Model : public AbstractModel, public ModelDimensions { * * @param x Reference to state variables * @param dx Reference to time derivative of states (DAE only) + * @param roots_found boolean indicators indicating whether roots were found at t0 by this fun */ - void initHeaviside(const AmiVector &x, const AmiVector &dx); + void initEvents(const AmiVector &x, const AmiVector &dx, + std::vector &roots_found); /** * @brief Get number of parameters wrt to which sensitivities are computed. @@ -1864,6 +1869,11 @@ class Model : public AbstractModel, public ModelDimensions { /** vector of bools indicating whether state variables are to be assumed to * be positive */ std::vector state_is_non_negative_; + + /** Vector of booleans indicating the initial boolean value for every event trigger function. Events at t0 + * can only trigger if the initial value is set to `false`. Must be specified during model compilation by + * setting the `initialValue` attribute of an event trigger. */ + std::vector root_initial_values_; /** boolean indicating whether any entry in stateIsNonNegative is `true` */ bool any_state_non_negative_ {false}; diff --git a/python/amici/ode_export.py b/python/amici/ode_export.py index 74a5d5abb7..b40e025524 100644 --- a/python/amici/ode_export.py +++ b/python/amici/ode_export.py @@ -822,7 +822,7 @@ def transform_dxdt_to_concentration(species_id, dxdt): else: args += ['value'] if symbol_name == SymbolId.EVENT: - args += ['state_update', 'event_observable'] + args += ['state_update', 'event_observable', 'initial_value'] if symbol_name == SymbolId.OBSERVABLE: args += ['transformation'] @@ -1589,34 +1589,43 @@ def _compute_equation(self, name: str) -> None: elif name == 'stau': self._eqs[name] = [ -self.eq('sroot')[ie, :] / self.eq('drootdt_total')[ie] + if not self.eq('drootdt_total')[ie].is_zero else + sp.zeros(*self.eq('sroot')[ie, :].shape) for ie in range(self.num_events()) ] elif name == 'deltasx': event_eqs = [] for ie, event in enumerate(self._events): - if event._state_update is not None: - # ====== chain rule for the state variables =============== - # get xdot with expressions back-substituted - tmp_eq = smart_multiply( + + tmp_eq = sp.zeros(self.num_states_solver(), self.num_par()) + + # only add stau part if trigger is time-dependent + if not self.eq('drootdt_total')[ie].is_zero: + tmp_eq += smart_multiply( (self.sym('xdot_old') - self.sym('xdot')), self.eq('stau')[ie]) - # construct an enhanced state sensitivity, which accounts - # for the time point sensitivity as well + + # only add deltax part if there is state update + if event._state_update is not None: + # partial derivative for the parameters + tmp_eq += self.eq('ddeltaxdp')[ie] + + # initial part of chain rule state variables tmp_dxdp = self.sym('sx') * sp.ones(1, self.num_par()) - tmp_dxdp += smart_multiply(self.sym('xdot'), - self.eq('stau')[ie]) + + # only add stau part if trigger is time-dependent + if not self.eq('drootdt_total')[ie].is_zero: + # chain rule for the time point + tmp_eq += smart_multiply(self.eq('ddeltaxdt')[ie], + self.eq('stau')[ie]) + + # additional part of chain rule state variables + tmp_dxdp += smart_multiply(self.sym('xdot'), + self.eq('stau')[ie]) + # finish chain rule for the state variables tmp_eq += smart_multiply(self.eq('ddeltaxdx')[ie], tmp_dxdp) - # ====== chain rule for the time point ==================== - tmp_eq += smart_multiply(self.eq('ddeltaxdt')[ie], - self.eq('stau')[ie]) - # ====== partial derivative for the parameters ============ - tmp_eq += self.eq('ddeltaxdp')[ie] - else: - tmp_eq = smart_multiply( - (self.eq('xdot_old') - self.eq('xdot')), - self.eq('stau')[ie]) event_eqs.append(tmp_eq) @@ -2937,6 +2946,11 @@ def _write_model_header_cpp(self) -> None: 'W_RECURSION_DEPTH': self.model._w_recursion_depth, 'QUADRATIC_LLH': 'true' if self.model._has_quadratic_nllh else 'false', + 'ROOT_INITIAL_VALUES': + ', '.join([ + 'true' if event.get_initial_value() else 'false' + for event in self.model._events + ]) } for func_name, func_info in self.functions.items(): diff --git a/python/amici/ode_model.py b/python/amici/ode_model.py index a5885216f3..ed494cf3ea 100644 --- a/python/amici/ode_model.py +++ b/python/amici/ode_model.py @@ -2,17 +2,7 @@ import sympy as sp -import numpy as np -import re -import shutil -import subprocess -import sys -import os -import copy import numbers -import logging -import itertools -import contextlib try: import pysb @@ -20,25 +10,11 @@ pysb = None from typing import ( - Callable, Optional, Union, List, Dict, Tuple, SupportsFloat, Sequence, - Set, Any + Optional, Union, Dict, SupportsFloat, Set ) -from dataclasses import dataclass -from string import Template -from sympy.matrices.immutable import ImmutableDenseMatrix -from sympy.matrices.dense import MutableDenseMatrix -from sympy.logic.boolalg import BooleanAtom -from itertools import chain -from .cxxcodeprinter import AmiciCxxCodePrinter, get_switch_statement - -from . import ( - amiciSwigPath, amiciSrcPath, amiciModulePath, __version__, __commit__, - sbml_import -) -from .logging import get_logger, log_execution_time, set_log_level -from .constants import SymbolId -from .import_utils import smart_subs_dict, toposort_symbols, \ - ObservableTransformation, generate_measurement_symbol, RESERVED_SYMBOLS + +from .import_utils import ObservableTransformation, \ + generate_measurement_symbol, RESERVED_SYMBOLS from .import_utils import cast_to_sym __all__ = [ @@ -46,6 +22,7 @@ 'ModelQuantity', 'Observable', 'Parameter', 'SigmaY', 'State' ] + class ModelQuantity: """ Base class for model components @@ -166,14 +143,6 @@ def __init__(self, self._ncoeff: sp.Expr = coefficients[state_id] super(ConservationLaw, self).__init__(identifier, name, value) - def get_state(self) -> sp.Symbol: - """ - Get the identifier of the state that this conservation law replaces - - :return: identifier of the state - """ - return self._state_id - def get_ncoeff(self, state_id) -> Union[sp.Expr, int, float]: """ Computes the normalized coefficient a_i/a_j where i is the index of @@ -211,10 +180,6 @@ class State(ModelQuantity): algebraic formula that defines the temporal derivative of this state """ - - _dt: Union[sp.Expr, None] = None - _conservation_law: Union[sp.Expr, None] = None - def __init__(self, identifier: sp.Symbol, name: str, @@ -276,7 +241,7 @@ def get_dt(self) -> sp.Expr: """ return self._dt - def get_free_symbols(self) -> Set[sp.Symbol]: + def get_free_symbols(self) -> Set[sp.Basic]: """ Gets the set of free symbols in time derivative and initial conditions @@ -307,8 +272,9 @@ def get_x_rdata(self): def get_dx_rdata_dx_solver(self, state_id): """ - Returns the expression that allows computation of ``dx_rdata_dx_solver`` for this - state, accounting for conservation laws. + Returns the expression that allows computation of + ``dx_rdata_dx_solver`` for this state, accounting for conservation + laws. :return: dx_rdata_dx_solver expression """ @@ -514,7 +480,8 @@ def __init__(self, name: str, value: sp.Expr, state_update: Union[sp.Expr, None], - event_observable: Union[sp.Expr, None]): + event_observable: Union[sp.Expr, None], + initial_value: Optional[bool] = True): """ Create a new Event instance. @@ -534,15 +501,30 @@ def __init__(self, :param event_observable: formula a potential observable linked to the event (None for Heaviside functions, empty events without observable) + + :param initial_value: + initial boolean value of the trigger function at t0. If set to + `False`, events may trigger at ``t==t0``, otherwise not. """ super(Event, self).__init__(identifier, name, value) # add the Event specific components self._state_update = state_update self._observable = event_observable + self._initial_value = initial_value + + def get_initial_value(self) -> bool: + """ + Return the initial value for the root function. + + :return: + initial value formula + """ + return self._initial_value def __eq__(self, other): """ Check equality of events at the level of trigger/root functions, as we need to collect unique root functions for ``roots.cpp`` """ - return self.get_val() == other.get_val() + return self.get_val() == other.get_val() and \ + (self.get_initial_value() == other.get_initial_value()) diff --git a/python/amici/sbml_import.py b/python/amici/sbml_import.py index 616f82e102..8b38e06f87 100644 --- a/python/amici/sbml_import.py +++ b/python/amici/sbml_import.py @@ -484,13 +484,6 @@ def check_event_support(self) -> None: if trigger_sbml.getMath() is None: logger.warning(f'Event {event_id} trigger has no trigger ' 'expression, so a dummy trigger will be set.') - if not trigger_sbml.getInitialValue(): - # True: event not executed if triggered at time == 0 - # (corresponding to AMICI default). Raise if set to False. - raise SBMLException( - f'Event {event_id} has a trigger that has an initial ' - 'value of False. This is currently not supported in AMICI.' - ) if not trigger_sbml.getPersistent(): raise SBMLException( @@ -1001,10 +994,16 @@ def _convert_event_assignment_parameter_targets_to_species(self): parameter_def = \ self.symbols[symbol_id].pop(parameter_target) if parameter_def is None: - raise AssertionError( - 'Unexpected error. The parameter target of an event ' - 'assignment could not be found.' + # this happens for parameters that have initial assignments + # or are assignment rule targets + par = self.sbml.getElementBySId(str(parameter_target)) + ia_init = self._get_element_initial_assignment( + par.getId() ) + parameter_def = { + 'name': par.getName() if par.isSetName() else par.getId(), + 'value': par.getValue() if ia_init is None else ia_init + } # Fixed parameters are added as species such that they can be # targets of events. self.symbols[SymbolId.SPECIES][parameter_target] = { @@ -1140,6 +1139,9 @@ def get_empty_bolus_value() -> sp.Float: 'value': trigger, 'state_update': sp.MutableDenseMatrix(bolus), 'event_observable': None, + 'initial_value': + trigger_sbml.getInitialValue() if trigger_sbml is not None + else True, } @log_execution_time('processing SBML observables', logger) diff --git a/src/forwardproblem.cpp b/src/forwardproblem.cpp index 8147edce05..773ffc5c4b 100644 --- a/src/forwardproblem.cpp +++ b/src/forwardproblem.cpp @@ -52,15 +52,21 @@ void ForwardProblem::workForwardProblem() { if (!preequilibrated_) model->initialize(x_, dx_, sx_, sdx_, solver->getSensitivityOrder() >= - SensitivityOrder::first); - else if (model->ne) - model->initHeaviside(x_, dx_); + SensitivityOrder::first, + roots_found_); + else if (model->ne) { + model->initEvents(x_, dx_, roots_found_); + } /* compute initial time and setup solver for (pre-)simulation */ auto t0 = model->t0(); if (presimulate) t0 -= edata->t_presim; solver->setup(t0, model, x_, dx_, sx_, sdx_); + + if (model->ne && std::any_of(roots_found_.begin(), roots_found_.end(), + [](int rf){return rf==1;})) + handleEvent(&t0, false, true); /* perform presimulation if necessary */ if (presimulate) { @@ -69,9 +75,14 @@ void ForwardProblem::workForwardProblem() { " is currently not implemented."); handlePresimulation(); t_ = model->t0(); - if (model->ne) - model->initHeaviside(x_, dx_); + if (model->ne) { + model->initEvents(x_, dx_, roots_found_); + if (std::any_of(roots_found_.begin(), roots_found_.end(), + [](int rf){return rf==1;})) + handleEvent(&t0, false, true); + } } + /* when computing adjoint sensitivity analysis with presimulation, we need to store sx after the reinitialization after preequilibration but before reinitialization after presimulation. As presimulation with ASA @@ -92,7 +103,6 @@ void ForwardProblem::workForwardProblem() { if (solver->computingFSA() || (solver->computingASA() && !presimulate )) sx_ = solver->getStateSensitivity(model->t0()); - /* store initial state and sensitivity*/ initial_state_ = getSimulationState(); @@ -114,7 +124,7 @@ void ForwardProblem::workForwardProblem() { /* clustering of roots => turn off rootfinding */ solver->turnOffRootFinding(); } else if (status == AMICI_ROOT_RETURN) { - handleEvent(&tlastroot_, false); + handleEvent(&tlastroot_, false, false); } } } @@ -138,7 +148,8 @@ void ForwardProblem::handlePresimulation() } -void ForwardProblem::handleEvent(realtype *tlastroot, const bool seflag) { +void ForwardProblem::handleEvent(realtype *tlastroot, const bool seflag, + const bool initial_event) { /* store Heaviside information at event occurrence */ model->froot(t_, x_, dx_, rootvals_); @@ -146,14 +157,14 @@ void ForwardProblem::handleEvent(realtype *tlastroot, const bool seflag) { discs_.push_back(t_); /* extract and store which events occurred */ - if (!seflag) { + if (!seflag && !initial_event) { solver->getRootInfo(roots_found_.data()); } root_idx_.push_back(roots_found_); rval_tmp_ = rootvals_; - if (!seflag) { + if (!seflag && !initial_event) { /* only check this in the first event fired, otherwise this will always * be true */ if (t_ == *tlastroot) { @@ -165,24 +176,24 @@ void ForwardProblem::handleEvent(realtype *tlastroot, const bool seflag) { *tlastroot = t_; } - if(model->nz>0) + if(model->nz > 0) storeEvent(); /* if we need to do forward sensitivities later on we need to store the old * x and the old xdot */ if (solver->getSensitivityOrder() >= SensitivityOrder::first) { /* store x and xdot to compute jump in sensitivities */ - x_old_ = x_; + x_old_.copy(x_); } if (solver->computingFSA()) { model->fxdot(t_, x_, dx_, xdot_); - xdot_old_ = xdot_; - dx_old_ = dx_; + xdot_old_.copy(xdot_); + dx_old_.copy(dx_); /* compute event-time derivative only for primary events, we get * into trouble with multiple simultaneously firing events here (but * is this really well defined then?), in that case just use the * last ie and hope for the best. */ - if (!seflag) { + if (!seflag && !initial_event) { for (int ie = 0; ie < model->ne; ie++) { if (roots_found_.at(ie) == 1) { /* only consider transitions false -> true */ @@ -190,6 +201,8 @@ void ForwardProblem::handleEvent(realtype *tlastroot, const bool seflag) { } } } + if (initial_event) // t0 has no parameter dependency + std::fill(stau_.begin(), stau_.end(), 0.0); } else if (solver->computingASA()) { /* store x to compute jump in discontinuity */ x_disc_.push_back(x_); @@ -197,7 +210,8 @@ void ForwardProblem::handleEvent(realtype *tlastroot, const bool seflag) { xdot_old_disc_.push_back(xdot_old_); } - model->updateHeaviside(roots_found_); + if (!initial_event) + model->updateHeaviside(roots_found_); applyEventBolus(); @@ -239,7 +253,7 @@ void ForwardProblem::handleEvent(realtype *tlastroot, const bool seflag) { "Secondary event was triggered. Depending on " "the bolus of the secondary event, forward " "sensitivities can be incorrect."); - handleEvent(tlastroot, true); + handleEvent(tlastroot, true, false); } /* only reinitialise in the first event fired */ diff --git a/src/model.ODE_template.cpp b/src/model.ODE_template.cpp index c5fb93198d..049124881d 100644 --- a/src/model.ODE_template.cpp +++ b/src/model.ODE_template.cpp @@ -53,6 +53,10 @@ std::array stateIdxsSolver = { TPL_STATE_IDXS_SOLVER_INITIALIZER_LIST }; +std::array rootInitialValues = { + TPL_ROOT_INITIAL_VALUES +}; + } // namespace model_TPL_MODELNAME } // namespace amici diff --git a/src/model.cpp b/src/model.cpp index fe6d138b70..e6097fa88f 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -130,7 +130,8 @@ static int setValueByIdRegex(std::vector const &ids, Model::Model(ModelDimensions const & model_dimensions, SimulationParameters simulation_parameters, - SecondOrderMode o2mode, std::vector idlist, std::vector z2event, + SecondOrderMode o2mode, std::vector idlist, + std::vector z2event, const bool pythonGenerated, const int ndxdotdp_explicit, const int ndxdotdx_explicit, const int w_recursion_depth) : ModelDimensions(model_dimensions), pythonGenerated(pythonGenerated), @@ -153,6 +154,8 @@ Model::Model(ModelDimensions const & model_dimensions, simulation_parameters_.pscale, state_.unscaledParameters); state_.fixedParameters = simulation_parameters_.fixedParameters; state_.plist = simulation_parameters_.plist; + + root_initial_values_.resize(ne, true); /* If Matlab wrapped: dxdotdp is a full AmiVector, if Python wrapped: dxdotdp_explicit and dxdotdp_implicit are CSC matrices @@ -241,7 +244,8 @@ bool operator==(const ModelDimensions &a, const ModelDimensions &b) { void Model::initialize(AmiVector &x, AmiVector &dx, AmiVectorArray &sx, - AmiVectorArray & /*sdx*/, bool computeSensitivities) { + AmiVectorArray & /*sdx*/, bool computeSensitivities, + std::vector &roots_found) { initializeStates(x); if (computeSensitivities) initializeStateSensitivities(sx, x); @@ -251,7 +255,7 @@ void Model::initialize(AmiVector &x, AmiVector &dx, AmiVectorArray &sx, fsdx0(); if (ne) - initHeaviside(x, dx); + initEvents(x, dx, roots_found); } void Model::initializeB(AmiVector &xB, AmiVector &dxB, AmiVector &xQB, @@ -298,14 +302,18 @@ void Model::initializeStateSensitivities(AmiVectorArray &sx, } } -void Model::initHeaviside(AmiVector const &x, AmiVector const &dx) { +void Model::initEvents(AmiVector const &x, AmiVector const &dx, + std::vector &roots_found) { std::vector rootvals(ne, 0.0); froot(simulation_parameters_.tstart_, x, dx, rootvals); + std::fill(roots_found.begin(), roots_found.end(), 0); for (int ie = 0; ie < ne; ie++) { if (rootvals.at(ie) < 0) { state_.h.at(ie) = 0.0; } else { state_.h.at(ie) = 1.0; + if (!root_initial_values_.at(ie)) // only false->true triggers event + roots_found.at(ie) = 1; } } } diff --git a/src/model_header.ODE_template.h b/src/model_header.ODE_template.h index 27e664bd8f..8a885a61ba 100644 --- a/src/model_header.ODE_template.h +++ b/src/model_header.ODE_template.h @@ -27,6 +27,7 @@ extern std::array stateIds; extern std::array observableIds; extern std::array expressionIds; extern std::array stateIdxsSolver; +extern std::array rootInitialValues; TPL_JY_DEF TPL_DJYDSIGMA_DEF @@ -129,7 +130,11 @@ class Model_TPL_MODELNAME : public amici::Model_ODE { TPL_NDXDOTDP_EXPLICIT, // ndxdotdp_explicit TPL_NDXDOTDX_EXPLICIT, // ndxdotdx_explicit TPL_W_RECURSION_DEPTH // w_recursion_depth - ) {} + ) { + root_initial_values_ = std::vector( + rootInitialValues.begin(), rootInitialValues.end() + ); + } /** * @brief Clone this model instance. diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp index d7959085f0..a2acec9596 100644 --- a/src/steadystateproblem.cpp +++ b/src/steadystateproblem.cpp @@ -188,9 +188,11 @@ void SteadystateProblem::initializeForwardProblem(int it, const Solver &solver, /* process solver handling for pre- or postequilibration */ if (it == -1) { /* solver was not run before, set up everything */ + auto roots_found = std::vector(model.ne, 0); model.initialize(state_.x, state_.dx, state_.sx, sdx_, - solver.getSensitivityOrder() >= - SensitivityOrder::first); + solver.getSensitivityOrder() >= + SensitivityOrder::first, + roots_found); state_.t = model.t0(); solver.setup(state_.t, &model, state_.x, state_.dx, state_.sx, sdx_); } else { diff --git a/swig/model.i b/swig/model.i index 81466a5c3b..d2c9f2eabd 100644 --- a/swig/model.i +++ b/swig/model.i @@ -32,7 +32,7 @@ using namespace amici; %ignore getObservable; %ignore getObservableSensitivity; %ignore getExpression; -%ignore initHeaviside; +%ignore initEvents; %ignore initialize; %ignore initializeB; %ignore initializeStateSensitivities; diff --git a/tests/testSBMLSuite.py b/tests/testSBMLSuite.py index 3187720a8f..767c4287fb 100755 --- a/tests/testSBMLSuite.py +++ b/tests/testSBMLSuite.py @@ -60,6 +60,8 @@ def test_sbml_testsuite_case( sensitivity_check_cases = { # parameter-dependent conservation laws '00783': 1.5e-2, + # initial events + '00995': 1e-3, } try: From c9524bbe748d545e5383f9aefc0e46022970948c Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Fri, 6 May 2022 16:57:07 +0200 Subject: [PATCH 09/10] Update to codecov/codecov-action@v3.1.0 (#1792) --- .github/workflows/test_petab_test_suite.yml | 12 ++++++++---- .github/workflows/test_python_cplusplus.yml | 9 +++++---- .github/workflows/test_sbml_semantic_test_suite.yml | 4 ++-- 3 files changed, 15 insertions(+), 10 deletions(-) diff --git a/.github/workflows/test_petab_test_suite.yml b/.github/workflows/test_petab_test_suite.yml index f089b544b9..4594317582 100644 --- a/.github/workflows/test_petab_test_suite.yml +++ b/.github/workflows/test_petab_test_suite.yml @@ -71,7 +71,8 @@ jobs: - name: Run PEtab-related unit tests run: | source ./build/venv/bin/activate \ - && pytest --cov-report=xml --cov=./ python/tests/test_*petab*.py + && pytest --cov-report=xml:coverage.xml \ + --cov=./ python/tests/test_*petab*.py # run test models - name: Run PEtab test suite @@ -79,12 +80,15 @@ jobs: run: | source ./build/venv/bin/activate \ && AMICI_PARALLEL_COMPILE=2 pytest -v \ - --cov-report=xml --cov-append --cov=amici tests/petab_test_suite/ + --cov-report=xml:coverage.xml \ + --cov-append \ + --cov=amici \ + tests/petab_test_suite/ - name: Codecov - uses: codecov/codecov-action@v2 + uses: codecov/codecov-action@v3.1.0 with: token: ${{ secrets.CODECOV_TOKEN }} - file: ./coverage.xml + file: coverage.xml flags: petab fail_ci_if_error: true diff --git a/.github/workflows/test_python_cplusplus.yml b/.github/workflows/test_python_cplusplus.yml index e62beea5a8..69b09e9dde 100644 --- a/.github/workflows/test_python_cplusplus.yml +++ b/.github/workflows/test_python_cplusplus.yml @@ -104,12 +104,13 @@ jobs: scripts/runNotebook.sh documentation/GettingStarted.ipynb - name: Codecov Python - uses: codecov/codecov-action@v2 + uses: codecov/codecov-action@v3.1.0 with: token: ${{ secrets.CODECOV_TOKEN }} - file: ./build/coverage_py.xml + file: build/coverage_py.xml flags: python fail_ci_if_error: true + verbose: true - name: lcov run: | @@ -122,10 +123,10 @@ jobs: && lcov -a coverage_cpp.info -a coverage_py.info -o coverage.info - name: Codecov CPP - uses: codecov/codecov-action@v2 + uses: codecov/codecov-action@v3.1.0 with: token: ${{ secrets.CODECOV_TOKEN }} - file: ./coverage.info + file: coverage.info flags: cpp fail_ci_if_error: true diff --git a/.github/workflows/test_sbml_semantic_test_suite.yml b/.github/workflows/test_sbml_semantic_test_suite.yml index 458f010fd0..7901ef8640 100644 --- a/.github/workflows/test_sbml_semantic_test_suite.yml +++ b/.github/workflows/test_sbml_semantic_test_suite.yml @@ -53,9 +53,9 @@ jobs: path: tests/amici-semantic-results - name: Codecov SBMLSuite - uses: codecov/codecov-action@v2 + uses: codecov/codecov-action@v3.1.0 with: token: ${{ secrets.CODECOV_TOKEN }} - file: ./coverage_SBMLSuite.xml + file: coverage_SBMLSuite.xml flags: sbmlsuite fail_ci_if_error: true From 19526971d3b007f97623f048832343dfc8f12ded Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Fri, 6 May 2022 14:00:55 +0200 Subject: [PATCH 10/10] Bump version number; update changelog --- CHANGELOG.md | 24 ++++++++++++++++++++++++ version.txt | 2 +- 2 files changed, 25 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 4bde64a0b1..07e7c3c6f8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,30 @@ ## v0.X Series +### v0.11.29 (2022-05-06) + +## What's Changed + +Features: +* Performance: Limit newton step convergence check by @FFroehlich in + https://github.com/AMICI-dev/AMICI/pull/1780 +* More informative NaN/Inf warnings by @dweindl in + https://github.com/AMICI-dev/AMICI/pull/1640 +* SBML import can now handle initial events by @FFroehlich in + https://github.com/AMICI-dev/AMICI/pull/1789 + +Fixes: +* Avoid error if no measurements in PEtab problem; fixed type handling in + PEtab parameter mapping by @dilpath in + https://github.com/AMICI-dev/AMICI/pull/1783 +* Fixed substitution of expressions in root and stau by @dilpath in + https://github.com/AMICI-dev/AMICI/pull/1784 +* Workaround for PEtab problems with state-dependent noise models by @dweindl + in https://github.com/AMICI-dev/AMICI/pull/1791 + +**Full Changelog**: https://github.com/AMICI-dev/AMICI/compare/v0.11.28...v0.11.29 + + ### v0.11.28 (2022-04-08) New features: diff --git a/version.txt b/version.txt index 2b71711fca..36e40e13e0 100644 --- a/version.txt +++ b/version.txt @@ -1 +1 @@ -0.11.28 +0.11.29