From c2ead08b9aec48bc2050b38e9f0cf7eb8a1203c5 Mon Sep 17 00:00:00 2001 From: FFroehlich Date: Wed, 16 Mar 2022 11:11:01 -0400 Subject: [PATCH 1/6] remove SPBCG solver --- include/amici/newton_solver.h | 90 -------------------- src/newton_solver.cpp | 152 +--------------------------------- src/steadystateproblem.cpp | 5 -- 3 files changed, 1 insertion(+), 246 deletions(-) diff --git a/include/amici/newton_solver.h b/include/amici/newton_solver.h index 01c033aeb6..9dd3c8a92f 100644 --- a/include/amici/newton_solver.h +++ b/include/amici/newton_solver.h @@ -254,96 +254,6 @@ class NewtonSolverSparse : public NewtonSolver { SUNLinearSolver linsol_ {nullptr}; }; -/** - * @brief The NewtonSolverIterative provides access to the iterative linear - * solver for the Newton method. - */ - -class NewtonSolverIterative : public NewtonSolver { - - public: - /** - * @brief Constructor, initializes all members with the provided objects - * @param t pointer to time variable - * @param x pointer to state variables - * @param model pointer to the model object - */ - NewtonSolverIterative(realtype *t, AmiVector *x, Model *model); - - ~NewtonSolverIterative() override = default; - - /** - * @brief Solves the linear system for the Newton step by passing it to - * linsolveSPBCG - * - * @param rhs containing the RHS of the linear system, will be - * overwritten by solution to the linear system - */ - void solveLinearSystem(AmiVector &rhs) override; - - /** - * Writes the Jacobian (J) for the Newton iteration and passes it to the linear - * solver. - * Also wraps around getSensis for iterative linear solver. - * - * @param ntry integer newton_try integer start number of Newton solver - * (1 or 2) - * @param nnewt integer number of current Newton step - */ - void prepareLinearSystem(int ntry, int nnewt) override; - - /** - * Writes the Jacobian (JB) for the Newton iteration and passes it to the linear - * solver. - * Also wraps around getSensis for iterative linear solver. - * - * @param ntry integer newton_try integer start number of Newton solver - * (1 or 2) - * @param nnewt integer number of current Newton step - */ - void prepareLinearSystemB(int ntry, int nnewt) override; - - /** - * Iterative linear solver created from SPILS BiCG-Stab. - * Solves the linear system within each Newton step if iterative solver is - * chosen. - * - * @param ntry integer newton_try integer start number of Newton solver - * (1 or 2) - * @param nnewt integer number of current Newton step - * @param ns_delta Newton step - */ - void linsolveSPBCG(int ntry, int nnewt, AmiVector &ns_delta); - - private: - /** number of tries */ - int newton_try_ {0}; - /** number of iterations */ - int i_newton_ {0}; - /** ??? */ - AmiVector ns_p_; - /** ??? */ - AmiVector ns_h_; - /** ??? */ - AmiVector ns_t_; - /** ??? */ - AmiVector ns_s_; - /** ??? */ - AmiVector ns_r_; - /** ??? */ - AmiVector ns_rt_; - /** ??? */ - AmiVector ns_v_; - /** ??? */ - AmiVector ns_Jv_; - /** ??? */ - AmiVector ns_tmp_; - /** ??? */ - AmiVector ns_Jdiag_; - /** temporary storage of Jacobian */ - SUNMatrixWrapper ns_J_; -}; - } // namespace amici diff --git a/src/newton_solver.cpp b/src/newton_solver.cpp index b770f04d8b..59b9cd6763 100644 --- a/src/newton_solver.cpp +++ b/src/newton_solver.cpp @@ -49,8 +49,7 @@ std::unique_ptr NewtonSolver::getSolver(realtype *t, AmiVector *x, throw NewtonFailure(AMICI_NOT_IMPLEMENTED, "getSolver"); case LinearSolver::SPBCG: - solver.reset(new NewtonSolverIterative(t, x, model)); - break; + throw NewtonFailure(AMICI_NOT_IMPLEMENTED, "getSolver"); case LinearSolver::SPTFQMR: throw NewtonFailure(AMICI_NOT_IMPLEMENTED, "getSolver"); @@ -71,8 +70,6 @@ std::unique_ptr NewtonSolver::getSolver(realtype *t, AmiVector *x, solver->damping_factor_mode_ = simulationSolver.getNewtonDampingFactorMode(); solver->damping_factor_lower_bound = simulationSolver.getNewtonDampingFactorLowerBound(); - if (simulationSolver.getLinearSolver() == LinearSolver::SPBCG) - solver->num_lin_steps_.resize(simulationSolver.getNewtonMaxSteps(), 0); return solver; } @@ -220,152 +217,5 @@ NewtonSolverSparse::~NewtonSolverSparse() { SUNLinSolFree_KLU(linsol_); } -/* ------------------------------------------------------------------------- */ -/* - Iterative linear solver------------------------------------------------ */ -/* ------------------------------------------------------------------------- */ - -NewtonSolverIterative::NewtonSolverIterative(realtype *t, AmiVector *x, - Model *model) - : NewtonSolver(t, x, model), ns_p_(model->nx_solver), - ns_h_(model->nx_solver), ns_t_(model->nx_solver), ns_s_(model->nx_solver), - ns_r_(model->nx_solver), ns_rt_(model->nx_solver), ns_v_(model->nx_solver), - ns_Jv_(model->nx_solver), ns_tmp_(model->nx_solver), - ns_Jdiag_(model->nx_solver), ns_J_(model->nx_solver, model->nx_solver) - { -} - -/* ------------------------------------------------------------------------- */ - -void NewtonSolverIterative::prepareLinearSystem(int ntry, int nnewt) { - newton_try_ = ntry; - i_newton_ = nnewt; - if (nnewt == -1) { - throw NewtonFailure(AMICI_NOT_IMPLEMENTED, - "Linear solver SPBCG does not support sensitivity " - "computation for steady state problems."); - } - - // Get the Jacobian and its diagonal for preconditioning - model_->fJ(*t_, 0.0, *x_, dx_, xdot_, ns_J_.get()); - ns_J_.refresh(); - model_->fJDiag(*t_, ns_Jdiag_, 0.0, *x_, dx_); - - // Ensure positivity of entries in ns_Jdiag - ns_p_.set(1.0); - ns_Jdiag_.abs(); - N_VCompare(1e-15, ns_Jdiag_.getNVector(), ns_tmp_.getNVector()); - linearSum(-1.0, ns_tmp_, 1.0, ns_p_, ns_tmp_); - linearSum(1.0, ns_Jdiag_, 1.0, ns_tmp_, ns_Jdiag_); -} - -/* ------------------------------------------------------------------------- */ - -void NewtonSolverIterative::prepareLinearSystemB(int ntry, int nnewt) { - newton_try_ = ntry; - i_newton_ = nnewt; - if (nnewt == -1) { - throw AmiException("Linear solver SPBCG does not support sensitivity " - "computation for steady state problems."); - } - - // Get the Jacobian and its diagonal for preconditioning - model_->fJB(*t_, 0.0, *x_, dx_, xB_, dxB_, xdot_, ns_J_.get()); - ns_J_.refresh(); - // Get the diagonal and ensure negativity of entries is ns_J. Note that diag(JB) = -diag(J). - model_->fJDiag(*t_, ns_Jdiag_, 0.0, *x_, dx_); - - ns_p_.set(1.0); - ns_Jdiag_.abs(); - N_VCompare(1e-15, ns_Jdiag_.getNVector(), ns_tmp_.getNVector()); - linearSum(-1.0, ns_tmp_, 1.0, ns_p_, ns_tmp_); - linearSum(1.0, ns_Jdiag_, 1.0, ns_tmp_, ns_Jdiag_); - ns_Jdiag_.minus(); -} - -/* ------------------------------------------------------------------------- */ - -void NewtonSolverIterative::solveLinearSystem(AmiVector &rhs) { - linsolveSPBCG(newton_try_, i_newton_, rhs); - rhs.minus(); -} - -/* ------------------------------------------------------------------------- */ - -void NewtonSolverIterative::linsolveSPBCG(int /*ntry*/, int nnewt, - AmiVector &ns_delta) { - xdot_ = ns_delta; - xdot_.minus(); - - // Initialize for linear solve - ns_p_.zero(); - ns_v_.zero(); - ns_delta.zero(); - ns_tmp_.zero(); - double rho = 1.0; - double omega = 1.0; - double alpha = 1.0; - - ns_J_.multiply(ns_Jv_, ns_delta); - - // ns_r = xdot - ns_Jv; - linearSum(-1.0, ns_Jv_, 1.0, xdot_, ns_r_); - ns_r_ /= ns_Jdiag_; - ns_rt_ = ns_r_; - - for (int i_linstep = 0; i_linstep < max_lin_steps_; i_linstep++) { - // Compute factors - double rho1 = rho; - rho = dotProd(ns_rt_, ns_r_); - double beta = rho * alpha / (rho1 * omega); - - // ns_p = ns_r + beta * (ns_p - omega * ns_v); - linearSum(1.0, ns_p_, -omega, ns_v_, ns_p_); - linearSum(1.0, ns_r_, beta, ns_p_, ns_p_); - - // ns_v = J * ns_p - ns_v_.zero(); - ns_J_.multiply(ns_v_, ns_p_); - ns_v_ /= ns_Jdiag_; - - // Compute factor - alpha = rho / dotProd(ns_rt_, ns_v_); - - // ns_h = ns_delta + alpha * ns_p; - linearSum(1.0, ns_delta, alpha, ns_p_, ns_h_); - // ns_s = ns_r - alpha * ns_v; - linearSum(1.0, ns_r_, -alpha, ns_v_, ns_s_); - - // ns_t = J * ns_s - ns_t_.zero(); - ns_J_.multiply(ns_t_, ns_s_); - ns_t_ /= ns_Jdiag_; - - // Compute factor - omega = dotProd(ns_t_, ns_s_) / dotProd(ns_t_, ns_t_); - - // ns_delta = ns_h + omega * ns_s; - linearSum(1.0, ns_h_, omega, ns_s_, ns_delta); - // ns_r = ns_s - omega * ns_t; - linearSum(1.0, ns_s_, -omega, ns_t_, ns_r_); - - // Compute the (unscaled) residual - ns_r_ *= ns_Jdiag_; - double res = sqrt(dotProd(ns_r_, ns_r_)); - - // Test convergence - if (res < atol_) { - // Write number of steps needed - num_lin_steps_.at(nnewt) = i_linstep + 1; - - // Return success - return; - } - - // Scale back - ns_r_ /= ns_Jdiag_; - } - throw NewtonFailure(AMICI_CONV_FAILURE, "linsolveSPBCG"); -} - } // namespace amici diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp index d7c3850e64..1bdae812e8 100644 --- a/src/steadystateproblem.cpp +++ b/src/steadystateproblem.cpp @@ -26,11 +26,6 @@ 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) { - /* maxSteps must be adapted if iterative linear solvers are used */ - if (solver.getLinearSolver() == LinearSolver::SPBCG) { - max_steps_ = solver.getNewtonMaxSteps(); - numlinsteps_.resize(2 * max_steps_, 0); - } /* Check for compatibility of options */ if (solver.getSensitivityMethod() == SensitivityMethod::forward && solver.getSensitivityMethodPreequilibration() == SensitivityMethod::adjoint && From 54a0b2feef188ba0c120b9e8994c1b64de362725 Mon Sep 17 00:00:00 2001 From: FFroehlich Date: Wed, 16 Mar 2022 13:31:03 -0400 Subject: [PATCH 2/6] remove test --- tests/cpp/expectedResults.h5 | Bin 4198368 -> 4198368 bytes tests/cpp/testOptions.h5 | Bin 345728 -> 345728 bytes .../generateTestConfig/example_steadystate.py | 11 ----------- 3 files changed, 11 deletions(-) diff --git a/tests/cpp/expectedResults.h5 b/tests/cpp/expectedResults.h5 index a24be851566e794a6dffd0f0564bed2660e06996..c88ecca4b2b6dd690197a5b5c27cd2a10d3cb848 100644 GIT binary patch delta 208 zcmWN|Sz18>0D#f&8Vecf{z;jYu|ia2T$&#KcRR2OtFQ@6;o&@<^I;b`U+u#n_(uYX zNK|4HmxLtcMs6h~Y01c)WF;qgDagG%NGR9-C`Eaal03_cl%*n7sYzWL(v+6G%A2&M rBVFl9Uk36nAMzcd%9ni0S~jvhuseqTI7c+- diff --git a/tests/cpp/testOptions.h5 b/tests/cpp/testOptions.h5 index e2bf3f658e83688bab6aa03e5fe7a117d1360861..a48d8ac1d9cae37d7869d4f867110da936e23b39 100644 GIT binary patch delta 29 kcmZqZ6>aDhogl)<*eKd6%GfH()GEr{D$24|lr`H80D1ifa{vGU delta 29 kcmZqZ6>aDhogl)<&?wp}%GfH()GEr{D$24|lr`H80D0*LasU7T diff --git a/tests/generateTestConfig/example_steadystate.py b/tests/generateTestConfig/example_steadystate.py index b9f2885bdd..fae863b3fc 100755 --- a/tests/generateTestConfig/example_steadystate.py +++ b/tests/generateTestConfig/example_steadystate.py @@ -54,16 +54,6 @@ def writeSensiForwardDense(filename): ex.writeToFile(filename, '/model_steadystate/sensiforwarddense/') -def writeNosensiSPBCG(filename): - ex = ExampleSteadystate() - - ex.modelOptions['ts'] = np.append(np.linspace(0, 100, 50), np.inf) - ex.solverOptions['sensi'] = 0 - ex.solverOptions['linsol'] = 7 - - ex.writeToFile(filename, '/model_steadystate/nosensiSPBCG/') - - def writeSensiForwardErrorInt(filename): ex = ExampleSteadystate() @@ -379,7 +369,6 @@ def main(): writeSensiForward(filename) writeSensiForwardPlist(filename) writeSensiForwardDense(filename) - writeNosensiSPBCG(filename) writeSensiForwardErrorInt(filename) writeSensiForwardErrorNewt(filename) writeSensiFwdNewtonPreeq(filename) From 2f3ca8b84b01925b4f4c79c78499246e386bd602 Mon Sep 17 00:00:00 2001 From: FFroehlich Date: Wed, 16 Mar 2022 14:19:40 -0400 Subject: [PATCH 3/6] fixup --- tests/cpp/steadystate/tests1.cpp | 6 ------ 1 file changed, 6 deletions(-) diff --git a/tests/cpp/steadystate/tests1.cpp b/tests/cpp/steadystate/tests1.cpp index 505f8f15ed..93f48ceb27 100644 --- a/tests/cpp/steadystate/tests1.cpp +++ b/tests/cpp/steadystate/tests1.cpp @@ -174,12 +174,6 @@ TEST(ExampleSteadystate, SensitivityForwardDense) amici::simulateVerifyWrite("/model_steadystate/sensiforwarddense/"); } -TEST(ExampleSteadystate, SensitivityForwardSPBCG) -{ - amici::simulateVerifyWrite( - "/model_steadystate/nosensiSPBCG/", 10 * TEST_ATOL, 10 * TEST_RTOL); -} - TEST(ExampleSteadystate, SensiFwdNewtonPreeq) { amici::simulateVerifyWrite("/model_steadystate/sensifwdnewtonpreeq/"); From 1f96d3d44c56d83128658c31e467d6bc61ec5f74 Mon Sep 17 00:00:00 2001 From: FFroehlich Date: Wed, 16 Mar 2022 14:26:17 -0400 Subject: [PATCH 4/6] retry deleting from hdf5 --- tests/cpp/expectedResults.h5 | Bin 4198368 -> 4198368 bytes 1 file changed, 0 insertions(+), 0 deletions(-) diff --git a/tests/cpp/expectedResults.h5 b/tests/cpp/expectedResults.h5 index c88ecca4b2b6dd690197a5b5c27cd2a10d3cb848..a68de6ea7eda08275a442ed7b639639b8429aacd 100644 GIT binary patch delta 304 zcmW;GNlwCG06^ik_)D3TNmNjfxmCadGN@2E010pe_AEI-W14jBs!J|m(k)@d)^LH| zK#t(kmwfx@Z!N`#Wy`i=2kf&iv70Fe$=lCFBOnX+XD)Pix)7(&o#q%zfpgDT2Gx(y zLYDrEyQ^>W$ypzRP!U2H5kwKg1uhXs0!gHh#uYNS#tm+fMGkkkM;-+{pokL6sGy1( o>d>I$5r0n&G||E{47AZf7cc0cj{$~w#T(u+!WgEvK2A2HUwcu71ONa4 delta 335 zcmXZPxlRIM07YSD@P|PWLWD*(@hwuy(Hc;>a#7shKM?=C(kPM}v zqNkzc4XEfjL2mLT_pHW!d4;|q8M0wD4kqHSA?uPkYlwtZvJKZ?A{v=n{`Z?xc1POx zbl=N@GW6Vxufd{AX;tW3c zagGaI;tBx-5yCZY5XLR;5Wzhj@Q5coBZ?T}NFa$6(ope&3|^7N8{Uya9t9NfffCB7 KnA)N`Ds2Iqn19v) From 30383b2c4ea09331df7634baf2dd2c20233f1317 Mon Sep 17 00:00:00 2001 From: FFroehlich Date: Wed, 16 Mar 2022 14:30:13 -0400 Subject: [PATCH 5/6] fixup remove options --- tests/cpp/testOptions.h5 | Bin 345728 -> 345728 bytes 1 file changed, 0 insertions(+), 0 deletions(-) diff --git a/tests/cpp/testOptions.h5 b/tests/cpp/testOptions.h5 index a48d8ac1d9cae37d7869d4f867110da936e23b39..25280e9cec208ed2aaecb54b27a2106c61723607 100644 GIT binary patch delta 81 zcmZqZ6>aDhogl)<&?wp}%GfH(v{jV3CX|t3`=n529cD)M>H8y?HKz+iGJlx3O=H?2 jCduuRQOxVta|SRnfPujDg*ME>?N2yawm;!yonr$4voaS| delta 106 zcmZqZ6>aDhogl)<*eKd6%GfH(v{jV3CX|tJ`=n529cD(3>H8y?HI)JwAmGGqW-x;h t%Krf63+$WzD1!OhwkYNk>v{6>i&OK8GlK)1oZZ{6aI$Q_!pXYE1^`ez95(;} From 400a6025946e71d04ddb9fdfa207f56e01a5a378 Mon Sep 17 00:00:00 2001 From: FFroehlich Date: Wed, 16 Mar 2022 15:30:14 -0400 Subject: [PATCH 6/6] remove python test --- python/tests/test_preequilibration.py | 16 ---------------- 1 file changed, 16 deletions(-) diff --git a/python/tests/test_preequilibration.py b/python/tests/test_preequilibration.py index 2af6da888b..5cd67a39d3 100644 --- a/python/tests/test_preequilibration.py +++ b/python/tests/test_preequilibration.py @@ -359,19 +359,3 @@ def test_newton_solver_equilibration(preeq_fixture): rdatas[settings[1]][variable], 1e-6, 1e-6 ).all(), variable - - # test failure for iterative linear solver with sensitivities - edata.fixedParametersPreequilibration = () - edata.t_presim = 0.0 - edata.fixedParametersPresimulation = () - - solver.setLinearSolver(amici.LinearSolver.SPBCG) - solver.setSensitivityMethod(amici.SensitivityMethod.adjoint) - solver.setSensitivityOrder(amici.SensitivityOrder.first) - model.setSteadyStateSensitivityMode( - amici.SteadyStateSensitivityMode.newtonOnly) - solver.setNewtonMaxSteps(10) - solver.setNewtonMaxLinearSteps(10) - rdata_spbcg = amici.runAmiciSimulation(model, solver, edata) - - assert rdata_spbcg['status'] == amici.AMICI_ERROR