Check condition number when computing sensitivities via newton

include/amici/newton_solver.h

@@ -97,6 +97,14 @@ class NewtonSolver {
      * overwritten by solution to the linear system
      */
     virtual void solveLinearSystem(AmiVector &rhs) = 0;
+
+    /**
+     * @brief Checks whether linear system is singular
+     *
+     * @return boolean indicating whether the linear system is singular
+     * (condition number < 1/machine precision)
+     */
+    virtual bool is_singular() const = 0;
 
     virtual ~NewtonSolver() = default;
 
@@ -185,6 +193,8 @@ class NewtonSolverDense : public NewtonSolver {
      * @param nnewt integer number of current Newton step
      */
     void prepareLinearSystemB(int ntry, int nnewt) override;
+
+    bool is_singular() const override;
 
   private:
     /** temporary storage of Jacobian */
@@ -245,6 +255,8 @@ class NewtonSolverSparse : public NewtonSolver {
      * @param nnewt integer number of current Newton step
      */
     void prepareLinearSystemB(int ntry, int nnewt) override;
+
+    bool is_singular() const override;
 
   private:
     /** temporary storage of Jacobian */

src/newton_solver.cpp

@@ -5,11 +5,19 @@
 
 #include "sunlinsol/sunlinsol_klu.h" // sparse solver
 #include "sunlinsol/sunlinsol_dense.h" // dense solver
+#include <sundials/sundials_math.h> // roundoffs
 
 #include <cstring>
 #include <ctime>
 #include <cmath>
 
+// taken from sundials/src/sunlinsol/klu/sunlinsol_klu.c
+#if defined(SUNDIALS_INT64_T)
+#define KLU_INDEXTYPE long int
+#else
+#define KLU_INDEXTYPE int
+#endif
+
 namespace amici {
 
 NewtonSolver::NewtonSolver(realtype *t, AmiVector *x, Model *model)
@@ -88,6 +96,10 @@ void NewtonSolver::getStep(int ntry, int nnewt, AmiVector &delta) {
 void NewtonSolver::computeNewtonSensis(AmiVectorArray &sx) {
     prepareLinearSystem(0, -1);
     model_->fdxdotdp(*t_, *x_, dx_);
+
+    if (is_singular());
+        model_->app->warningF("AMICI:newton",
+                              "Jacobian is singular at steadystate, sensitivities may be inaccurate");
 
     if (model_->pythonGenerated) {
         for (int ip = 0; ip < model_->nplist(); ip++) {
@@ -157,6 +169,17 @@ void NewtonSolverDense::solveLinearSystem(AmiVector &rhs) {
 
 /* ------------------------------------------------------------------------- */
 
+bool NewtonSolverDense::is_singular() const {
+    // dense solver doesn't have any implementation for rcond/condest, so use
+    // sparse solver interface, not the most efficient solution, but who is
+    // concerned about speed and used the dense solver anyways ¯\_(ツ)_/¯
+    auto sparse_solver = new NewtonSolverSparse(t_, x_, model_);
+    sparse_solver->prepareLinearSystem(0,0);
+    auto is_singular = sparse_solver->is_singular();
+    delete sparse_solver;
+    return is_singular;
+}
+
 NewtonSolverDense::~NewtonSolverDense() {
     if(linsol_)
         SUNLinSolFree_Dense(linsol_);
@@ -212,6 +235,31 @@ void NewtonSolverSparse::solveLinearSystem(AmiVector &rhs) {
 
 /* ------------------------------------------------------------------------- */
 
+bool NewtonSolverSparse::is_singular() const {
+    // adapted from SUNLinSolSetup_KLU in sunlinsol/klu/sunlinsol_klu.c
+    auto content = (SUNLinearSolverContent_KLU)(linsol_->content);
+    // first cheap check via rcond
+    int status = sun_klu_rcond(content->symbolic, content->numeric,
+                               &content->common);
+    if(status != AMICI_SUCCESS)
+        throw NewtonFailure(status, "sun_klu_rcond");
+
+    auto precision = SUNRpowerR(UNIT_ROUNDOFF, 2.0/3.0);
+
+    if (content->common.rcond < precision) {
+        // cheap check indicates singular, expensive check via condest
+        status = sun_klu_condest((KLU_INDEXTYPE*) SM_INDEXPTRS_S(Jtmp_.get()),
+                                 SM_DATA_S(Jtmp_.get()),
+                                 content->symbolic,
+                                 content->numeric,
+                                 &content->common);
+        if(status != AMICI_SUCCESS)
+            throw NewtonFailure(status, "sun_klu_rcond");
+        return content->common.condest > 1.0/precision;
+    }
+    return false;
+}
+
 NewtonSolverSparse::~NewtonSolverSparse() {
     if(linsol_)
         SUNLinSolFree_KLU(linsol_);