From 4d37031a026651ad06db8f806b266450ad586eca Mon Sep 17 00:00:00 2001 From: Polina Lakrisenko Date: Tue, 15 Mar 2022 19:33:43 +0100 Subject: [PATCH 01/15] Steady State Sensitivity modes --- include/amici/defines.h | 3 ++- src/steadystateproblem.cpp | 21 +++++++++++++++------ 2 files changed, 17 insertions(+), 7 deletions(-) diff --git a/include/amici/defines.h b/include/amici/defines.h index 852b0c4370..8efeacf491 100644 --- a/include/amici/defines.h +++ b/include/amici/defines.h @@ -180,7 +180,8 @@ enum class NonlinearSolverIteration { /** Sensitivity computation mode in steadyStateProblem */ enum class SteadyStateSensitivityMode { newtonOnly, - simulationFSA + integrationOnly, + integrateIfNewtonFails }; /** State in which the steady state computation finished */ diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp index d7c3850e64..c058bcc185 100644 --- a/src/steadystateproblem.cpp +++ b/src/steadystateproblem.cpp @@ -250,11 +250,16 @@ void SteadystateProblem::computeSteadyStateQuadrature(NewtonSolver *newtonSolver We therefore compute the integral over xB first and then do a matrix-vector multiplication */ + auto sensitivityMode = model->getSteadyStateSensitivityMode(); + /* Try to compute the analytical solution for quadrature algebraically */ - getQuadratureByLinSolve(newtonSolver, model); + if (sensitivityMode == SteadyStateSensitivityMode::newtonOnly + || sensitivityMode == SteadyStateSensitivityMode::integrateIfNewtonFails) + getQuadratureByLinSolve(newtonSolver, model); - /* Analytical solution didn't work, perform simulation instead */ - if (!hasQuadrature()) + /* Perform simulation */ + if (sensitivityMode == SteadyStateSensitivityMode::integrationOnly || + (sensitivityMode == SteadyStateSensitivityMode::integrateIfNewtonFails && !hasQuadrature())) getQuadratureBySimulation(solver, model); /* If analytic solution and integration did not work, throw an Exception */ @@ -377,7 +382,7 @@ bool SteadystateProblem::getSensitivityFlag(const Model *model, bool forwardSensisAlreadyComputed = solver->getSensitivityOrder() >= SensitivityOrder::first && steady_state_status_[1] == SteadyStateStatus::success && - model->getSteadyStateSensitivityMode() == SteadyStateSensitivityMode::simulationFSA; + model->getSteadyStateSensitivityMode() == SteadyStateSensitivityMode::integrationOnly; bool simulationStartedInSteadystate = steady_state_status_[0] == SteadyStateStatus::success && @@ -402,7 +407,11 @@ bool SteadystateProblem::getSensitivityFlag(const Model *model, /* When we're creating a new solver object */ bool needForwardSensiAtCreation = needForwardSensisPreeq && - model->getSteadyStateSensitivityMode() == SteadyStateSensitivityMode::simulationFSA; + model->getSteadyStateSensitivityMode() == SteadyStateSensitivityMode::integrationOnly; + + if (model->getSteadyStateSensitivityMode() == SteadyStateSensitivityMode::integrateIfNewtonFails) + throw AmiException(" SteadyStateSensitivityMode is not compatible with " + "forward sensitivities approach"); /* Check if we need to store sensis */ switch (context) { @@ -651,7 +660,7 @@ void SteadystateProblem::runSteadystateSimulation(const Solver *solver, numsteps_.at(1) = sim_steps; if (solver->getSensitivityOrder() > SensitivityOrder::none && model->getSteadyStateSensitivityMode() == - SteadyStateSensitivityMode::simulationFSA) + SteadyStateSensitivityMode::integrationOnly) sx_ = solver->getStateSensitivity(t_); } } From 3e9df5bc9b35f087fe6baa634ab6f9a0be118c64 Mon Sep 17 00:00:00 2001 From: Polina Lakrisenko Date: Mon, 21 Mar 2022 15:37:00 +0100 Subject: [PATCH 02/15] integrationOnly and integrateIfNewtonFails in forward sensitivity case --- src/steadystateproblem.cpp | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp index c058bcc185..aaf38e51f9 100644 --- a/src/steadystateproblem.cpp +++ b/src/steadystateproblem.cpp @@ -26,6 +26,14 @@ SteadystateProblem::SteadystateProblem(const Solver &solver, const Model &model) xB_(model.nJ * model.nx_solver), xQ_(model.nJ * model.nx_solver), xQB_(model.nplist()), xQBdot_(model.nplist()), dJydx_(model.nJ * model.nx_solver * model.nt(), 0.0) { + + /* Turn off Newton's method if SteadyStateSensitivityMode is integrationOnly + in forward sensitivity case */ + if (solver.getSensitivityMethod() == SensitivityMethod::forward && + model.getSteadyStateSensitivityMode == SteadyStateSensitivityMode::integrationOnly) { + solver.setNewtonMaxSteps(0); + } + /* maxSteps must be adapted if iterative linear solvers are used */ if (solver.getLinearSolver() == LinearSolver::SPBCG) { max_steps_ = solver.getNewtonMaxSteps(); @@ -382,7 +390,8 @@ bool SteadystateProblem::getSensitivityFlag(const Model *model, bool forwardSensisAlreadyComputed = solver->getSensitivityOrder() >= SensitivityOrder::first && steady_state_status_[1] == SteadyStateStatus::success && - model->getSteadyStateSensitivityMode() == SteadyStateSensitivityMode::integrationOnly; + (model->getSteadyStateSensitivityMode() == SteadyStateSensitivityMode::integrationOnly || + model->getSteadyStateSensitivityMode() == SteadyStateSensitivityMode::integrateIfNewtonFails); bool simulationStartedInSteadystate = steady_state_status_[0] == SteadyStateStatus::success && @@ -407,11 +416,9 @@ bool SteadystateProblem::getSensitivityFlag(const Model *model, /* When we're creating a new solver object */ bool needForwardSensiAtCreation = needForwardSensisPreeq && - model->getSteadyStateSensitivityMode() == SteadyStateSensitivityMode::integrationOnly; - - if (model->getSteadyStateSensitivityMode() == SteadyStateSensitivityMode::integrateIfNewtonFails) - throw AmiException(" SteadyStateSensitivityMode is not compatible with " - "forward sensitivities approach"); + (model->getSteadyStateSensitivityMode() == SteadyStateSensitivityMode::integrationOnly || + model->getSteadyStateSensitivityMode() == SteadyStateSensitivityMode::integrateIfNewtonFails + ); /* Check if we need to store sensis */ switch (context) { @@ -659,8 +666,10 @@ void SteadystateProblem::runSteadystateSimulation(const Solver *solver, } else { numsteps_.at(1) = sim_steps; if (solver->getSensitivityOrder() > SensitivityOrder::none && + (model->getSteadyStateSensitivityMode() == + SteadyStateSensitivityMode::integrationOnly || model->getSteadyStateSensitivityMode() == - SteadyStateSensitivityMode::integrationOnly) + SteadyStateSensitivityMode::integrateIfNewtonFails)) sx_ = solver->getStateSensitivity(t_); } } From 12eede5da0b22004fd72e217459128d5d07a9434 Mon Sep 17 00:00:00 2001 From: Polina Lakrisenko Date: Mon, 21 Mar 2022 16:25:22 +0100 Subject: [PATCH 03/15] turn off newton method in case integrationOnly and forward sensitivity --- include/amici/steadystateproblem.h | 2 +- src/steadystateproblem.cpp | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/include/amici/steadystateproblem.h b/include/amici/steadystateproblem.h index 420d91b9b3..350046a935 100644 --- a/include/amici/steadystateproblem.h +++ b/include/amici/steadystateproblem.h @@ -29,7 +29,7 @@ class SteadystateProblem { * @param solver Solver instance * @param model Model instance */ - explicit SteadystateProblem(const Solver &solver, Model &model); + explicit SteadystateProblem(Solver &solver, Model &model); /** * @brief Handles steady state computation in the forward case: diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp index 32b0f18857..18640fa73f 100644 --- a/src/steadystateproblem.cpp +++ b/src/steadystateproblem.cpp @@ -20,7 +20,7 @@ constexpr realtype conv_thresh = 1.0; namespace amici { -SteadystateProblem::SteadystateProblem(const Solver &solver, Model &model) +SteadystateProblem::SteadystateProblem(Solver &solver, Model &model) : delta_(model.nx_solver), ewt_(model.nx_solver), ewtQB_(model.nplist()), x_old_(model.nx_solver), xdot_(model.nx_solver), sdx_(model.nx_solver, model.nplist()), xB_(model.nJ * model.nx_solver), @@ -53,7 +53,7 @@ SteadystateProblem::SteadystateProblem(const Solver &solver, Model &model) /* Turn off Newton's method if SteadyStateSensitivityMode is integrationOnly in forward sensitivity case */ if (solver.getSensitivityMethod() == SensitivityMethod::forward && - model.getSteadyStateSensitivityMode == + model.getSteadyStateSensitivityMode() == SteadyStateSensitivityMode::integrationOnly) { solver.setNewtonMaxSteps(0); } From 2709d51f4cb8131decd067c4f6c6507a3a4e8f93 Mon Sep 17 00:00:00 2001 From: Polina Lakrisenko Date: Mon, 21 Mar 2022 17:20:59 +0100 Subject: [PATCH 04/15] python tests --- python/tests/test_compare_conservation_laws_sbml.py | 2 +- python/tests/test_preequilibration.py | 7 +++---- python/tests/test_pysb.py | 4 ++-- python/tests/test_sbml_import.py | 3 +-- tests/petab_test_suite/test_petab_suite.py | 4 ++-- 5 files changed, 9 insertions(+), 11 deletions(-) diff --git a/python/tests/test_compare_conservation_laws_sbml.py b/python/tests/test_compare_conservation_laws_sbml.py index 2c760ae014..d0758020f2 100644 --- a/python/tests/test_compare_conservation_laws_sbml.py +++ b/python/tests/test_compare_conservation_laws_sbml.py @@ -135,7 +135,7 @@ def test_compare_conservation_laws_sbml(models, edata_fixture): # ----- compare simulations wo edata, sensi = 0, states and sensis ------- model_without_cl.setSteadyStateSensitivityMode( - amici.SteadyStateSensitivityMode.simulationFSA + amici.SteadyStateSensitivityMode.integrationOnly ) # run simulations diff --git a/python/tests/test_preequilibration.py b/python/tests/test_preequilibration.py index 5cd67a39d3..d3a01d7c5a 100644 --- a/python/tests/test_preequilibration.py +++ b/python/tests/test_preequilibration.py @@ -286,7 +286,8 @@ def test_equilibration_methods_with_adjoints(preeq_fixture): rdatas = {} equil_meths = [amici.SteadyStateSensitivityMode.newtonOnly, - amici.SteadyStateSensitivityMode.simulationFSA] + amici.SteadyStateSensitivityMode.integrationOnly, + amici.SteadyStateSensitivityMode.integrateIfNewtonFails] sensi_meths = [amici.SensitivityMethod.forward, amici.SensitivityMethod.adjoint] settings = itertools.product(equil_meths, sensi_meths) @@ -333,7 +334,7 @@ def test_newton_solver_equilibration(preeq_fixture): edata.setObservedDataStdDev(np.hstack([stdy, stdy[0]])) rdatas = {} - settings = [amici.SteadyStateSensitivityMode.simulationFSA, + settings = [amici.SteadyStateSensitivityMode.integrationOnly, amici.SteadyStateSensitivityMode.newtonOnly] for equil_meth in settings: @@ -343,8 +344,6 @@ def test_newton_solver_equilibration(preeq_fixture): model.setSteadyStateSensitivityMode(equil_meth) if equil_meth == amici.SteadyStateSensitivityMode.newtonOnly: solver.setNewtonMaxSteps(10) - else: - solver.setNewtonMaxSteps(0) # add rdatas rdatas[equil_meth] = amici.runAmiciSimulation(model, solver, edata) diff --git a/python/tests/test_pysb.py b/python/tests/test_pysb.py index 0fd050b596..5129b6021d 100644 --- a/python/tests/test_pysb.py +++ b/python/tests/test_pysb.py @@ -201,7 +201,7 @@ def get_data(model): solver = model.getSolver() model.setTimepoints(np.linspace(0, 60, 61)) model.setSteadyStateSensitivityMode( - amici.SteadyStateSensitivityMode.simulationFSA + amici.SteadyStateSensitivityMode.integrateIfNewtonFails ) rdata = amici.runAmiciSimulation(model, solver) @@ -220,7 +220,7 @@ def get_results(model, edata): edata.reinitializeFixedParameterInitialStates = True model.setTimepoints(np.linspace(0, 60, 61)) model.setSteadyStateSensitivityMode( - amici.SteadyStateSensitivityMode.simulationFSA + amici.SteadyStateSensitivityMode.integrateIfNewtonFails ) return amici.runAmiciSimulation(model, solver, edata) diff --git a/python/tests/test_sbml_import.py b/python/tests/test_sbml_import.py index 4637e149ad..64bf2610f1 100644 --- a/python/tests/test_sbml_import.py +++ b/python/tests/test_sbml_import.py @@ -205,10 +205,9 @@ def test_presimulation(sbml_example_presimulation_module): """Test 'presimulation' test model""" model = sbml_example_presimulation_module.getModel() solver = model.getSolver() - solver.setNewtonMaxSteps(0) model.setTimepoints(np.linspace(0, 60, 61)) model.setSteadyStateSensitivityMode( - amici.SteadyStateSensitivityMode.simulationFSA + amici.SteadyStateSensitivityMode.integrationOnly ) solver.setSensitivityOrder(amici.SensitivityOrder.first) model.setReinitializeFixedParameterInitialStates(True) diff --git a/tests/petab_test_suite/test_petab_suite.py b/tests/petab_test_suite/test_petab_suite.py index 59c13c72b4..16cfa9136b 100755 --- a/tests/petab_test_suite/test_petab_suite.py +++ b/tests/petab_test_suite/test_petab_suite.py @@ -16,7 +16,7 @@ from amici.petab_import import import_petab_problem, PysbPetabProblem from amici.petab_objective import ( simulate_petab, rdatas_to_measurement_df, create_parameterized_edatas) -from amici import SteadyStateSensitivityMode_simulationFSA +from amici import SteadyStateSensitivityMode_integrateIfNewtonFails, SteadyStateSensitivityMode_newtonOnly logger = get_logger(__name__, logging.DEBUG) set_log_level(get_logger("amici.petab_import"), logging.DEBUG) @@ -134,7 +134,7 @@ def check_derivatives(problem: petab.Problem, model: amici.Model) -> None: # Required for case 9 to not fail in # amici::NewtonSolver::computeNewtonSensis model.setSteadyStateSensitivityMode( - SteadyStateSensitivityMode_simulationFSA) + SteadyStateSensitivityMode_integrateIfNewtonFails) for edata in create_parameterized_edatas( amici_model=model, petab_problem=problem, From e4d922255fecc044723af4effd63be62f907a9c6 Mon Sep 17 00:00:00 2001 From: Polina Lakrisenko Date: Mon, 21 Mar 2022 17:24:31 +0100 Subject: [PATCH 05/15] fix --- tests/petab_test_suite/test_petab_suite.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/petab_test_suite/test_petab_suite.py b/tests/petab_test_suite/test_petab_suite.py index 16cfa9136b..131de9f918 100755 --- a/tests/petab_test_suite/test_petab_suite.py +++ b/tests/petab_test_suite/test_petab_suite.py @@ -16,7 +16,7 @@ from amici.petab_import import import_petab_problem, PysbPetabProblem from amici.petab_objective import ( simulate_petab, rdatas_to_measurement_df, create_parameterized_edatas) -from amici import SteadyStateSensitivityMode_integrateIfNewtonFails, SteadyStateSensitivityMode_newtonOnly +from amici import SteadyStateSensitivityMode_integrateIfNewtonFails logger = get_logger(__name__, logging.DEBUG) set_log_level(get_logger("amici.petab_import"), logging.DEBUG) From a515db381261d172cd73f89c8932ce0ba8ae63e7 Mon Sep 17 00:00:00 2001 From: Polina Lakrisenko Date: Mon, 21 Mar 2022 18:24:23 +0100 Subject: [PATCH 06/15] Apply suggestions from code review MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Fabian Fröhlich --- tests/petab_test_suite/test_petab_suite.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/petab_test_suite/test_petab_suite.py b/tests/petab_test_suite/test_petab_suite.py index 131de9f918..97d4bc6866 100755 --- a/tests/petab_test_suite/test_petab_suite.py +++ b/tests/petab_test_suite/test_petab_suite.py @@ -16,7 +16,7 @@ from amici.petab_import import import_petab_problem, PysbPetabProblem from amici.petab_objective import ( simulate_petab, rdatas_to_measurement_df, create_parameterized_edatas) -from amici import SteadyStateSensitivityMode_integrateIfNewtonFails +from amici import SteadyStateSensitivityMode logger = get_logger(__name__, logging.DEBUG) set_log_level(get_logger("amici.petab_import"), logging.DEBUG) @@ -134,7 +134,7 @@ def check_derivatives(problem: petab.Problem, model: amici.Model) -> None: # Required for case 9 to not fail in # amici::NewtonSolver::computeNewtonSensis model.setSteadyStateSensitivityMode( - SteadyStateSensitivityMode_integrateIfNewtonFails) + SteadyStateSensitivityMode.integrateIfNewtonFails) for edata in create_parameterized_edatas( amici_model=model, petab_problem=problem, From f974f7fba4d36928b554eb4ee9813f38282f84f6 Mon Sep 17 00:00:00 2001 From: Polina Lakrisenko Date: Mon, 21 Mar 2022 20:07:25 +0100 Subject: [PATCH 07/15] turn off newtin in findSteadyState --- include/amici/steadystateproblem.h | 2 +- src/steadystateproblem.cpp | 24 +++++++++++------------- 2 files changed, 12 insertions(+), 14 deletions(-) diff --git a/include/amici/steadystateproblem.h b/include/amici/steadystateproblem.h index 350046a935..420d91b9b3 100644 --- a/include/amici/steadystateproblem.h +++ b/include/amici/steadystateproblem.h @@ -29,7 +29,7 @@ class SteadystateProblem { * @param solver Solver instance * @param model Model instance */ - explicit SteadystateProblem(Solver &solver, Model &model); + explicit SteadystateProblem(const Solver &solver, Model &model); /** * @brief Handles steady state computation in the forward case: diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp index 18640fa73f..bd017b2a3b 100644 --- a/src/steadystateproblem.cpp +++ b/src/steadystateproblem.cpp @@ -20,7 +20,7 @@ constexpr realtype conv_thresh = 1.0; namespace amici { -SteadystateProblem::SteadystateProblem(Solver &solver, Model &model) +SteadystateProblem::SteadystateProblem(const Solver &solver, Model &model) : delta_(model.nx_solver), ewt_(model.nx_solver), ewtQB_(model.nplist()), x_old_(model.nx_solver), xdot_(model.nx_solver), sdx_(model.nx_solver, model.nplist()), xB_(model.nJ * model.nx_solver), @@ -48,15 +48,7 @@ SteadystateProblem::SteadystateProblem(Solver &solver, Model &model) solver.getSensitivityOrder() > SensitivityOrder::none) throw AmiException("Preequilibration using adjoint sensitivities " "is not compatible with using forward " - "sensitivities during simulation"); - - /* Turn off Newton's method if SteadyStateSensitivityMode is integrationOnly - in forward sensitivity case */ - if (solver.getSensitivityMethod() == SensitivityMethod::forward && - model.getSteadyStateSensitivityMode() == - SteadyStateSensitivityMode::integrationOnly) { - solver.setNewtonMaxSteps(0); - } + "sensitivities during simulation"); } void SteadystateProblem::workSteadyStateProblem(Solver *solver, Model *model, @@ -98,16 +90,21 @@ void SteadystateProblem::workSteadyStateBackwardProblem( void SteadystateProblem::findSteadyState(const Solver *solver, Model *model, int it) { steady_state_status_.resize(3, SteadyStateStatus::not_run); + bool turnOffNewton = solver->getSensitivityMethod() == + SensitivityMethod::forward && model->getSteadyStateSensitivityMode() == + SteadyStateSensitivityMode::integrationOnly; + /* First, try to run the Newton solver */ - findSteadyStateByNewtonsMethod(model, false); + if (!turnOffNewton) + findSteadyStateByNewtonsMethod(model, false); /* Newton solver didn't work, so try to simulate to steady state */ if (!checkSteadyStateSuccess()) findSteadyStateBySimulation(solver, model, it); /* Simulation didn't work, retry the Newton solver from last sim state. */ - if (!checkSteadyStateSuccess()) + if (!turnOffNewton && !checkSteadyStateSuccess()) findSteadyStateByNewtonsMethod(model, true); /* Nothing worked, throw an as informative error as possible */ @@ -254,7 +251,8 @@ void SteadystateProblem::computeSteadyStateQuadrature(const Solver *solver, /* Perform simulation */ if (sensitivityMode == SteadyStateSensitivityMode::integrationOnly || - (sensitivityMode == SteadyStateSensitivityMode::integrateIfNewtonFails && !hasQuadrature())) + (sensitivityMode == SteadyStateSensitivityMode::integrateIfNewtonFails + && !hasQuadrature())) getQuadratureBySimulation(solver, model); /* If analytic solution and integration did not work, throw an Exception */ From f5d5e08efd5a0fdd564912c395e145d184996af5 Mon Sep 17 00:00:00 2001 From: Polina Lakrisenko Date: Wed, 23 Mar 2022 15:49:57 +0100 Subject: [PATCH 08/15] fix test_compare_conservation_laws_sbml --- python/tests/test_compare_conservation_laws_sbml.py | 2 +- src/steadystateproblem.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/python/tests/test_compare_conservation_laws_sbml.py b/python/tests/test_compare_conservation_laws_sbml.py index d0758020f2..1cabe90c3b 100644 --- a/python/tests/test_compare_conservation_laws_sbml.py +++ b/python/tests/test_compare_conservation_laws_sbml.py @@ -135,7 +135,7 @@ def test_compare_conservation_laws_sbml(models, edata_fixture): # ----- compare simulations wo edata, sensi = 0, states and sensis ------- model_without_cl.setSteadyStateSensitivityMode( - amici.SteadyStateSensitivityMode.integrationOnly + amici.SteadyStateSensitivityMode.integrateIfNewtonFails ) # run simulations diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp index bd017b2a3b..9dd16ee86f 100644 --- a/src/steadystateproblem.cpp +++ b/src/steadystateproblem.cpp @@ -48,7 +48,7 @@ SteadystateProblem::SteadystateProblem(const Solver &solver, Model &model) solver.getSensitivityOrder() > SensitivityOrder::none) throw AmiException("Preequilibration using adjoint sensitivities " "is not compatible with using forward " - "sensitivities during simulation"); + "sensitivities during simulation"); } void SteadystateProblem::workSteadyStateProblem(Solver *solver, Model *model, From ee0ba61fb71560c627f26b44f9bcf010eab69313 Mon Sep 17 00:00:00 2001 From: Polina Lakrisenko Date: Tue, 29 Mar 2022 16:39:35 +0200 Subject: [PATCH 09/15] fix --- src/steadystateproblem.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp index f368adf432..0f682e66f4 100644 --- a/src/steadystateproblem.cpp +++ b/src/steadystateproblem.cpp @@ -94,8 +94,8 @@ void SteadystateProblem::workSteadyStateBackwardProblem( void SteadystateProblem::findSteadyState(const Solver &solver, Model &model, int it) { steady_state_status_.resize(3, SteadyStateStatus::not_run); - bool turnOffNewton = solver->getSensitivityMethod() == - SensitivityMethod::forward && model->getSteadyStateSensitivityMode() == + bool turnOffNewton = solver.getSensitivityMethod() == + SensitivityMethod::forward && model.getSteadyStateSensitivityMode() == SteadyStateSensitivityMode::integrationOnly; @@ -377,9 +377,9 @@ bool SteadystateProblem::getSensitivityFlag(const Model &model, bool forwardSensisAlreadyComputed = solver.getSensitivityOrder() >= SensitivityOrder::first && steady_state_status_[1] == SteadyStateStatus::success && - (model->getSteadyStateSensitivityMode() == + (model.getSteadyStateSensitivityMode() == SteadyStateSensitivityMode::integrationOnly || - model->getSteadyStateSensitivityMode() == + model.getSteadyStateSensitivityMode() == SteadyStateSensitivityMode::integrateIfNewtonFails); bool simulationStartedInSteadystate = @@ -407,9 +407,9 @@ bool SteadystateProblem::getSensitivityFlag(const Model &model, /* When we're creating a new solver object */ bool needForwardSensiAtCreation = needForwardSensisPreeq && - (model->getSteadyStateSensitivityMode() == + (model.getSteadyStateSensitivityMode() == SteadyStateSensitivityMode::integrationOnly || - model->getSteadyStateSensitivityMode() == + model.getSteadyStateSensitivityMode() == SteadyStateSensitivityMode::integrateIfNewtonFails ); From 57f7fc41b6a9b4cc5432aff86b2c741fb5c9c16a Mon Sep 17 00:00:00 2001 From: Polina Lakrisenko Date: Tue, 29 Mar 2022 16:49:01 +0200 Subject: [PATCH 10/15] fix --- src/steadystateproblem.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp index 0f682e66f4..a2a27ef1aa 100644 --- a/src/steadystateproblem.cpp +++ b/src/steadystateproblem.cpp @@ -248,7 +248,7 @@ void SteadystateProblem::computeSteadyStateQuadrature(const Solver &solver, We therefore compute the integral over xB first and then do a matrix-vector multiplication */ - auto sensitivityMode = model->getSteadyStateSensitivityMode(); + auto sensitivityMode = model.getSteadyStateSensitivityMode(); /* Try to compute the analytical solution for quadrature algebraically */ if (sensitivityMode == SteadyStateSensitivityMode::newtonOnly From 5f8dce132732190396b5d1681d8693b612dab11d Mon Sep 17 00:00:00 2001 From: Polina Lakrisenko Date: Wed, 30 Mar 2022 00:45:32 +0200 Subject: [PATCH 11/15] take preequilibration into account --- src/steadystateproblem.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp index a2a27ef1aa..f528ddaabf 100644 --- a/src/steadystateproblem.cpp +++ b/src/steadystateproblem.cpp @@ -94,10 +94,11 @@ void SteadystateProblem::workSteadyStateBackwardProblem( void SteadystateProblem::findSteadyState(const Solver &solver, Model &model, int it) { steady_state_status_.resize(3, SteadyStateStatus::not_run); - bool turnOffNewton = solver.getSensitivityMethod() == - SensitivityMethod::forward && model.getSteadyStateSensitivityMode() == - SteadyStateSensitivityMode::integrationOnly; - + bool turnOffNewton = model.getSteadyStateSensitivityMode() == + SteadyStateSensitivityMode::integrationOnly && + (it == -1 && solver.getSensitivityMethodPreequilibration() == + SensitivityMethod::forward || solver.getSensitivityMethod() == + SensitivityMethod::forward); /* First, try to run the Newton solver */ if (!turnOffNewton) From 2c197a4358d2b1f3b3e665d15dde6c9720fabc0e Mon Sep 17 00:00:00 2001 From: Polina Lakrisenko Date: Wed, 30 Mar 2022 12:10:30 +0200 Subject: [PATCH 12/15] add extra parentheses --- src/steadystateproblem.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp index f528ddaabf..be16d56981 100644 --- a/src/steadystateproblem.cpp +++ b/src/steadystateproblem.cpp @@ -96,8 +96,8 @@ void SteadystateProblem::findSteadyState(const Solver &solver, Model &model, steady_state_status_.resize(3, SteadyStateStatus::not_run); bool turnOffNewton = model.getSteadyStateSensitivityMode() == SteadyStateSensitivityMode::integrationOnly && - (it == -1 && solver.getSensitivityMethodPreequilibration() == - SensitivityMethod::forward || solver.getSensitivityMethod() == + ((it == -1 && solver.getSensitivityMethodPreequilibration() == + SensitivityMethod::forward) || solver.getSensitivityMethod() == SensitivityMethod::forward); /* First, try to run the Newton solver */ From c6001d3093f317993c9fdccea649774dc0ba26a4 Mon Sep 17 00:00:00 2001 From: Polina Lakrisenko Date: Mon, 4 Apr 2022 18:15:12 +0200 Subject: [PATCH 13/15] fix test_compare_conservation_laws_sbml --- .../test_compare_conservation_laws_sbml.py | 66 +++++++++++-------- 1 file changed, 39 insertions(+), 27 deletions(-) diff --git a/python/tests/test_compare_conservation_laws_sbml.py b/python/tests/test_compare_conservation_laws_sbml.py index 1cabe90c3b..52c8a5621b 100644 --- a/python/tests/test_compare_conservation_laws_sbml.py +++ b/python/tests/test_compare_conservation_laws_sbml.py @@ -10,7 +10,7 @@ def edata_fixture(): """edata is generated to test pre- and postequilibration""" edata_pre = amici.ExpData(2, 0, 0, - np.array([0., 0.1, 0.2, 0.5, 1., 2., 5., 10.])) + np.array([0., 0.1, 0.2, 0.5, 1., 2., 5., 10.])) edata_pre.setObservedData([1.5] * 16) edata_pre.fixedParameters = np.array([5., 20.]) edata_pre.fixedParametersPreequilibration = np.array([0., 10.]) @@ -18,13 +18,15 @@ def edata_fixture(): # edata for postequilibration edata_post = amici.ExpData(2, 0, 0, - np.array([float('inf')] * 3)) + np.array([float('inf')] * 3)) edata_post.setObservedData([0.75] * 6) edata_post.fixedParameters = np.array([7.5, 30.]) # edata with both equilibrations edata_full = amici.ExpData(2, 0, 0, - np.array([0., 0., 0., 1., 2., 2., 4., float('inf'), float('inf')])) + np.array( + [0., 0., 0., 1., 2., 2., 4., float('inf'), + float('inf')])) edata_full.setObservedData([3.14] * 18) edata_full.fixedParameters = np.array([1., 2.]) edata_full.fixedParametersPreequilibration = np.array([3., 4.]) @@ -42,7 +44,7 @@ def models(): sbml_importer = amici.SbmlImporter(sbml_file) # Name of the model that will also be the name of the python module - model_name = model_output_dir ='model_constant_species' + model_name = model_output_dir = 'model_constant_species' model_name_cl = model_output_dir_cl = 'model_constant_species_cl' # Define constants, observables, sigmas @@ -67,9 +69,11 @@ def models(): compute_conservation_laws=False) # load both models - model_without_cl_module = amici.import_model_module(model_name, + model_without_cl_module = amici.import_model_module( + model_name, module_path=os.path.abspath(model_name)) - model_with_cl_module = amici.import_model_module(model_name_cl, + model_with_cl_module = amici.import_model_module( + model_name_cl, module_path=os.path.abspath(model_name_cl)) # get the models and return @@ -81,8 +85,8 @@ def models(): def get_results(model, edata=None, sensi_order=0, sensi_meth=amici.SensitivityMethod.forward, sensi_meth_preeq=amici.SensitivityMethod.forward, + stst_sensi_mode=amici.SteadyStateSensitivityMode.newtonOnly, reinitialize_states=False): - # set model and data properties model.setReinitializeFixedParameterInitialStates(reinitialize_states) @@ -92,6 +96,7 @@ def get_results(model, edata=None, sensi_order=0, solver.setSensitivityOrder(sensi_order) solver.setSensitivityMethodPreequilibration(sensi_meth_preeq) solver.setSensitivityMethod(sensi_meth) + model.setSteadyStateSensitivityMode(stst_sensi_mode) if edata is None: model.setTimepoints(np.linspace(0, 5, 101)) else: @@ -178,43 +183,50 @@ def test_compare_conservation_laws_sbml(models, edata_fixture): def test_adjoint_pre_and_post_equilibration(edata_fixture): # get both models - model_module = amici.import_model_module('model_constant_species', + model_module = amici.import_model_module( + 'model_constant_species', module_path=os.path.abspath('model_constant_species')) model = model_module.getModel() - model_module_cl = amici.import_model_module('model_constant_species_cl', + model_module_cl = amici.import_model_module( + 'model_constant_species_cl', module_path=os.path.abspath('model_constant_species_cl')) model_cl = model_module_cl.getModel() # check gradient with and without state reinitialization for edata in edata_fixture: for reinit in [False, True]: - # --- compare different ways of preequilibration, full rank Jacobian --------- + # --- compare different ways of preequilibration, full rank Jacobian ------ # forward preequilibration, forward simulation - rff_cl = get_results(model_cl, edata=edata, sensi_order=1, - sensi_meth=amici.SensitivityMethod.forward, - sensi_meth_preeq=amici.SensitivityMethod.forward, - reinitialize_states=reinit) + rff_cl = get_results( + model_cl, edata=edata, sensi_order=1, + sensi_meth=amici.SensitivityMethod.forward, + sensi_meth_preeq=amici.SensitivityMethod.forward, + reinitialize_states=reinit) # forward preequilibration, adjoint simulation - rfa_cl = get_results(model_cl, edata=edata, sensi_order=1, - sensi_meth=amici.SensitivityMethod.adjoint, - sensi_meth_preeq=amici.SensitivityMethod.forward, - reinitialize_states=reinit) + rfa_cl = get_results( + model_cl, edata=edata, sensi_order=1, + sensi_meth=amici.SensitivityMethod.adjoint, + sensi_meth_preeq=amici.SensitivityMethod.forward, + reinitialize_states=reinit) # adjoint preequilibration, adjoint simulation - raa_cl = get_results(model_cl, edata=edata, sensi_order=1, - sensi_meth=amici.SensitivityMethod.adjoint, - sensi_meth_preeq=amici.SensitivityMethod.adjoint, - reinitialize_states=reinit) + raa_cl = get_results( + model_cl, edata=edata, sensi_order=1, + sensi_meth=amici.SensitivityMethod.adjoint, + sensi_meth_preeq=amici.SensitivityMethod.adjoint, + reinitialize_states=reinit) # assert all are close assert np.isclose(rff_cl['sllh'], rfa_cl['sllh']).all() assert np.isclose(rfa_cl['sllh'], raa_cl['sllh']).all() assert np.isclose(raa_cl['sllh'], rff_cl['sllh']).all() - # --- compare fully adjoint approach to simulation with singular Jacobian ---- - raa = get_results(model, edata=edata, sensi_order=1, - sensi_meth=amici.SensitivityMethod.adjoint, - sensi_meth_preeq=amici.SensitivityMethod.adjoint, - reinitialize_states=reinit) + # --- compare fully adjoint approach to simulation with singular Jacobian ---- + raa = get_results( + model, edata=edata, sensi_order=1, + sensi_meth=amici.SensitivityMethod.adjoint, + sensi_meth_preeq=amici.SensitivityMethod.adjoint, + stst_sensi_mode=amici.SteadyStateSensitivityMode.integrateIfNewtonFails, + reinitialize_states=reinit) # assert gradients are close (quadrature tolerances are laxer) assert np.isclose(raa_cl['sllh'], raa['sllh'], 1e-5, 1e-5).all() From 249f1a7bcc9b940d9121171d39a77ea11cf24b0a Mon Sep 17 00:00:00 2001 From: Polina Lakrisenko Date: Mon, 4 Apr 2022 18:49:59 +0200 Subject: [PATCH 14/15] fix test_compare_conservation_laws_sbml --- python/tests/test_compare_conservation_laws_sbml.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/python/tests/test_compare_conservation_laws_sbml.py b/python/tests/test_compare_conservation_laws_sbml.py index 52c8a5621b..fc0644ef0a 100644 --- a/python/tests/test_compare_conservation_laws_sbml.py +++ b/python/tests/test_compare_conservation_laws_sbml.py @@ -139,9 +139,6 @@ def test_compare_conservation_laws_sbml(models, edata_fixture): assert np.isclose(rdata[field], rdata_cl[field]).all(), field # ----- compare simulations wo edata, sensi = 0, states and sensis ------- - model_without_cl.setSteadyStateSensitivityMode( - amici.SteadyStateSensitivityMode.integrateIfNewtonFails - ) # run simulations edata, _, _ = edata_fixture @@ -159,7 +156,10 @@ def test_compare_conservation_laws_sbml(models, edata_fixture): # run simulations rdata_cl = get_results(model_with_cl, edata=edata, sensi_order=1) assert rdata_cl['status'] == amici.AMICI_SUCCESS - rdata = get_results(model_without_cl, edata=edata, sensi_order=1) + rdata = get_results( + model_without_cl, edata=edata, sensi_order=1, + stst_sensi_mode=amici.SteadyStateSensitivityMode.integrateIfNewtonFails + ) assert rdata['status'] == amici.AMICI_SUCCESS # check that steady state computation succeeded only by sim in full model assert (rdata['preeq_status'] == np.array([-3, 1, 0])).all() From b39b9254ab407a2e0c9fb9edd807cab9d0f0c0fe Mon Sep 17 00:00:00 2001 From: Polina Lakrisenko Date: Tue, 5 Apr 2022 01:27:08 +0200 Subject: [PATCH 15/15] quick notebook fix --- .../ExampleEquilibrationLogic.ipynb | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/python/examples/example_constant_species/ExampleEquilibrationLogic.ipynb b/python/examples/example_constant_species/ExampleEquilibrationLogic.ipynb index f65005b1d0..1e24e62a56 100644 --- a/python/examples/example_constant_species/ExampleEquilibrationLogic.ipynb +++ b/python/examples/example_constant_species/ExampleEquilibrationLogic.ipynb @@ -496,9 +496,9 @@ } ], "source": [ - "# Call simulation with singular Jacobian and simulationFSA mode\n", + "# Call simulation with singular Jacobian and integrateIfNewtonFails mode\n", "model.setTimepoints(np.full(1, np.inf))\n", - "model.setSteadyStateSensitivityMode(amici.SteadyStateSensitivityMode.simulationFSA)\n", + "model.setSteadyStateSensitivityMode(amici.SteadyStateSensitivityMode.integrateIfNewtonFails)\n", "solver = model.getSolver()\n", "solver.setNewtonMaxSteps(10)\n", "solver.setSensitivityMethod(amici.SensitivityMethod.forward)\n", @@ -883,7 +883,7 @@ ], "source": [ "# Singluar Jacobian, use simulation\n", - "model.setSteadyStateSensitivityMode(amici.SteadyStateSensitivityMode.simulationFSA)\n", + "model.setSteadyStateSensitivityMode(amici.SteadyStateSensitivityMode.integrateIfNewtonFails)\n", "solver = model.getSolver()\n", "solver.setNewtonMaxSteps(10)\n", "solver.setSensitivityMethod(amici.SensitivityMethod.forward)\n", @@ -1135,7 +1135,7 @@ ], "source": [ "# Non-singular Jacobian, use simulaiton\n", - "model_reduced.setSteadyStateSensitivityMode(amici.SteadyStateSensitivityMode.simulationFSA)\n", + "model_reduced.setSteadyStateSensitivityMode(amici.SteadyStateSensitivityMode.integrateIfNewtonFails)\n", "solver_reduced = model_reduced.getSolver()\n", "solver_reduced.setNewtonMaxSteps(0)\n", "solver_reduced.setSensitivityMethod(amici.SensitivityMethod.forward)\n", @@ -1186,7 +1186,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.2" + "version": "3.8.10" }, "toc": { "base_numbering": 1,