From b3325e221ebfa253e17007e03032bbafa9b9f761 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Mon, 22 May 2023 16:21:12 +0200 Subject: [PATCH 1/2] Sensitivities for AlgebraicRule models --- src/model_dae.cpp | 42 +++++++++++++++++++++++++++++++++++------- 1 file changed, 35 insertions(+), 7 deletions(-) diff --git a/src/model_dae.cpp b/src/model_dae.cpp index 43a8d81313..10c2823802 100644 --- a/src/model_dae.cpp +++ b/src/model_dae.cpp @@ -164,14 +164,36 @@ void Model_DAE::fdxdotdp( const realtype t, const const_N_Vector x, const const_N_Vector dx ) { auto x_pos = computeX_pos(x); + fdwdp(t, N_VGetArrayPointerConst(x_pos)); if (pythonGenerated) { - // python generated, not yet implemented for DAEs - throw AmiException("Wrapping of DAEs is not yet implemented from Python" - ); + // python generated + derived_state_.dxdotdp_explicit.zero(); + derived_state_.dxdotdp_implicit.zero(); + if (derived_state_.dxdotdp_explicit.capacity()) { + fdxdotdp_explicit_colptrs(derived_state_.dxdotdp_explicit); + fdxdotdp_explicit_rowvals(derived_state_.dxdotdp_explicit); + fdxdotdp_explicit( + derived_state_.dxdotdp_explicit.data(), t, + N_VGetArrayPointerConst(x_pos), + state_.unscaledParameters.data(), state_.fixedParameters.data(), + state_.h.data(), N_VGetArrayPointerConst(dx), derived_state_.w_.data() + ); + } + + fdxdotdw(t, x_pos, dx); + /* Sparse matrix multiplication + dxdotdp_implicit += dxdotdw * dwdp */ + derived_state_.dxdotdw_.sparse_multiply( + derived_state_.dxdotdp_implicit, derived_state_.dwdp_ + ); + + derived_state_.dxdotdp_full.sparse_add( + derived_state_.dxdotdp_explicit, 1.0, + derived_state_.dxdotdp_implicit, 1.0 + ); } else { // matlab generated - fdwdp(t, N_VGetArrayPointerConst(x_pos)); for (int ip = 0; ip < nplist(); ip++) { N_VConst(0.0, derived_state_.dxdotdp.getNVector(ip)); @@ -496,9 +518,15 @@ void Model_DAE::fsxdot( } if (pythonGenerated) { - // python generated, not yet implemented for DAEs - throw AmiException("Wrapping of DAEs is not yet implemented from Python" - ); + /* copy dxdotdp and the implicit version over */ + // initialize + N_VConst(0.0, sxdot); + realtype* sxdot_tmp = N_VGetArrayPointer(sxdot); + + derived_state_.dxdotdp_full.scatter( + plist(ip), 1.0, nullptr, gsl::make_span(sxdot_tmp, nx_solver), 0, + nullptr, 0 + ); } else { /* copy dxdotdp over */ N_VScale(1.0, derived_state_.dxdotdp.getNVector(ip), sxdot); From 74b9d3da4d27d4030a85beb12d7a76faf4d3a743 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Mon, 22 May 2023 20:09:22 +0200 Subject: [PATCH 2/2] test --- tests/testSBMLSuite.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/tests/testSBMLSuite.py b/tests/testSBMLSuite.py index f11870b60d..3e1cb96d0a 100755 --- a/tests/testSBMLSuite.py +++ b/tests/testSBMLSuite.py @@ -76,12 +76,15 @@ def test_sbml_testsuite_case(test_number, result_path, sbml_semantic_cases_dir): current_test_path, test_id, model_dir, - generate_sensitivity_code=test_id in sensitivity_check_cases, + generate_sensitivity_code=True, ) settings = read_settings_file(current_test_path, test_id) atol, rtol = apply_settings(settings, solver, model, test_id) + solver.setSensitivityOrder(amici.SensitivityOrder.first) + solver.setSensitivityMethod(amici.SensitivityMethod.forward) + # simulate model rdata = amici.runAmiciSimulation(model, solver) if rdata["status"] != amici.AMICI_SUCCESS: @@ -98,8 +101,6 @@ def test_sbml_testsuite_case(test_number, result_path, sbml_semantic_cases_dir): # check sensitivities for selected models if epsilon := sensitivity_check_cases.get(test_id): - solver.setSensitivityOrder(amici.SensitivityOrder.first) - solver.setSensitivityMethod(amici.SensitivityMethod.forward) check_derivatives(model, solver, epsilon=epsilon) except amici.sbml_import.SBMLException as err: