Reuse solver residual in calibrator

src/model/Block.cpp

@@ -62,7 +62,6 @@ void Block::update_solution(
 void Block::post_solve(Eigen::Matrix<double, Eigen::Dynamic, 1>& y) {}
 
 void Block::update_gradient(Eigen::SparseMatrix<double>& jacobian,
-                            Eigen::Matrix<double, Eigen::Dynamic, 1>& residual,
                             Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha,
                             std::vector<double>& y, std::vector<double>& dy) {
   throw std::runtime_error("Gradient calculation not implemented for block " +

src/model/Block.h

@@ -252,20 +252,22 @@ class Block {
   virtual void post_solve(Eigen::Matrix<double, Eigen::Dynamic, 1>& y);
 
   /**
-   * @brief Set the gradient of the block contributions with respect to the
+   * @brief Set the gradient of the block residual with respect to the
    * parameters
    *
+   * Only the Jacobian of the residual with respect to the parameters is
+   * assembled here. The residual itself is identical to the one assembled by
+   * the solver (see SparseSystem::update_residual) and is therefore reused
+   * during calibration instead of being redefined.
+   *
    * @param jacobian Jacobian with respect to the parameters
    * @param alpha Current parameter vector
-   * @param residual Residual with respect to the parameters
    * @param y Current solution
    * @param dy Time-derivative of the current solution
    */
-  virtual void update_gradient(
-      Eigen::SparseMatrix<double>& jacobian,
-      Eigen::Matrix<double, Eigen::Dynamic, 1>& residual,
-      Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha, std::vector<double>& y,
-      std::vector<double>& dy);
+  virtual void update_gradient(Eigen::SparseMatrix<double>& jacobian,
+                               Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha,
+                               std::vector<double>& y, std::vector<double>& dy);
 
   /**
    * @brief Number of triplets of element

src/model/BloodVessel.cpp

@@ -58,47 +58,37 @@ void BloodVessel::update_solution(
 
 void BloodVessel::update_gradient(
     Eigen::SparseMatrix<double>& jacobian,
-    Eigen::Matrix<double, Eigen::Dynamic, 1>& residual,
     Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha, std::vector<double>& y,
     std::vector<double>& dy) {
-  auto y0 = y[global_var_ids[0]];
   auto y1 = y[global_var_ids[1]];
-  auto y2 = y[global_var_ids[2]];
-  auto y3 = y[global_var_ids[3]];
 
   auto dy0 = dy[global_var_ids[0]];
   auto dy1 = dy[global_var_ids[1]];
   auto dy3 = dy[global_var_ids[3]];
 
   auto resistance = alpha[global_param_ids[ParamId::RESISTANCE]];
   auto capacitance = alpha[global_param_ids[ParamId::CAPACITANCE]];
-  auto inductance = alpha[global_param_ids[ParamId::INDUCTANCE]];
   double stenosis_coeff = 0.0;
 
   if (global_param_ids.size() > 3) {
     stenosis_coeff = alpha[global_param_ids[ParamId::STENOSIS_COEFFICIENT]];
   }
   auto stenosis_resistance = stenosis_coeff * fabs(y1);
 
-  jacobian.coeffRef(global_eqn_ids[0], global_param_ids[0]) = -y1;
-  jacobian.coeffRef(global_eqn_ids[0], global_param_ids[2]) = -dy3;
+  jacobian.coeffRef(global_eqn_ids[0], global_param_ids[0]) = y1;
+  jacobian.coeffRef(global_eqn_ids[0], global_param_ids[2]) = dy3;
 
   if (global_param_ids.size() > 3) {
-    jacobian.coeffRef(global_eqn_ids[0], global_param_ids[3]) = -fabs(y1) * y1;
+    jacobian.coeffRef(global_eqn_ids[0], global_param_ids[3]) = fabs(y1) * y1;
   }
 
-  jacobian.coeffRef(global_eqn_ids[1], global_param_ids[0]) = capacitance * dy1;
+  jacobian.coeffRef(global_eqn_ids[1], global_param_ids[0]) =
+      -capacitance * dy1;
   jacobian.coeffRef(global_eqn_ids[1], global_param_ids[1]) =
-      -dy0 + (resistance + 2 * stenosis_resistance) * dy1;
+      dy0 - (resistance + 2 * stenosis_resistance) * dy1;
 
   if (global_param_ids.size() > 3) {
     jacobian.coeffRef(global_eqn_ids[1], global_param_ids[3]) =
-        2.0 * capacitance * fabs(y1) * dy1;
+        -2.0 * capacitance * fabs(y1) * dy1;
   }
-
-  residual(global_eqn_ids[0]) =
-      y0 - (resistance + stenosis_resistance) * y1 - y2 - inductance * dy3;
-  residual(global_eqn_ids[1]) =
-      y1 - y3 - capacitance * dy0 +
-      capacitance * (resistance + 2.0 * stenosis_resistance) * dy1;
 }

src/model/BloodVessel.h

@@ -202,12 +202,10 @@ class BloodVessel : public Block {
    *
    * @param jacobian Jacobian with respect to the parameters
    * @param alpha Current parameter vector
-   * @param residual Residual with respect to the parameters
    * @param y Current solution
    * @param dy Time-derivative of the current solution
    */
   void update_gradient(Eigen::SparseMatrix<double>& jacobian,
-                       Eigen::Matrix<double, Eigen::Dynamic, 1>& residual,
                        Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha,
                        std::vector<double>& y,
                        std::vector<double>& dy) override;

src/model/BloodVesselCRL.cpp

@@ -50,38 +50,19 @@ void BloodVesselCRL::update_solution(
 
 void BloodVesselCRL::update_gradient(
     Eigen::SparseMatrix<double>& jacobian,
-    Eigen::Matrix<double, Eigen::Dynamic, 1>& residual,
-    Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha, std::vector<double>& y,
+    Eigen::Matrix<double, Eigen::Dynamic, 1>& /*alpha*/, std::vector<double>& y,
     std::vector<double>& dy) {
-  auto y0 = y[global_var_ids[0]];
-  auto y1 = y[global_var_ids[1]];
-  auto y2 = y[global_var_ids[2]];
   auto y3 = y[global_var_ids[3]];
 
   auto dy0 = dy[global_var_ids[0]];
-  auto dy1 = dy[global_var_ids[1]];
   auto dy3 = dy[global_var_ids[3]];
 
-  auto resistance = alpha[global_param_ids[ParamId::RESISTANCE]];
-  auto capacitance = alpha[global_param_ids[ParamId::CAPACITANCE]];
-  auto inductance = alpha[global_param_ids[ParamId::INDUCTANCE]];
-  double stenosis_coeff = 0.0;
+  jacobian.coeffRef(global_eqn_ids[0], global_param_ids[0]) = y3;
+  jacobian.coeffRef(global_eqn_ids[0], global_param_ids[2]) = dy3;
 
   if (global_param_ids.size() > 3) {
-    stenosis_coeff = alpha[global_param_ids[ParamId::STENOSIS_COEFFICIENT]];
+    jacobian.coeffRef(global_eqn_ids[0], global_param_ids[3]) = fabs(y3) * y3;
   }
-  auto stenosis_resistance = stenosis_coeff * fabs(y3);
 
-  jacobian.coeffRef(global_eqn_ids[0], global_param_ids[0]) = -y3;
-  jacobian.coeffRef(global_eqn_ids[0], global_param_ids[2]) = -dy3;
-
-  if (global_param_ids.size() > 3) {
-    jacobian.coeffRef(global_eqn_ids[0], global_param_ids[3]) = -fabs(y3) * y3;
-  }
-
-  jacobian.coeffRef(global_eqn_ids[1], global_param_ids[1]) = -dy0;
-
-  residual(global_eqn_ids[0]) =
-      y0 - (resistance + stenosis_resistance) * y3 - y2 - inductance * dy3;
-  residual(global_eqn_ids[1]) = y1 - y3 - capacitance * dy0;
+  jacobian.coeffRef(global_eqn_ids[1], global_param_ids[1]) = dy0;
 }

src/model/BloodVesselCRL.h

@@ -181,12 +181,10 @@ class BloodVesselCRL : public Block {
    *
    * @param jacobian Jacobian with respect to the parameters
    * @param alpha Current parameter vector
-   * @param residual Residual with respect to the parameters
    * @param y Current solution
    * @param dy Time-derivative of the current solution
    */
   void update_gradient(Eigen::SparseMatrix<double>& jacobian,
-                       Eigen::Matrix<double, Eigen::Dynamic, 1>& residual,
                        Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha,
                        std::vector<double>& y,
                        std::vector<double>& dy) override;

src/model/BloodVesselJunction.cpp

@@ -54,43 +54,24 @@ void BloodVesselJunction::update_solution(
 
 void BloodVesselJunction::update_gradient(
     Eigen::SparseMatrix<double>& jacobian,
-    Eigen::Matrix<double, Eigen::Dynamic, 1>& residual,
-    Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha, std::vector<double>& y,
+    Eigen::Matrix<double, Eigen::Dynamic, 1>& /*alpha*/, std::vector<double>& y,
     std::vector<double>& dy) {
-  auto p_in = y[global_var_ids[0]];
-  auto q_in = y[global_var_ids[1]];
-
-  residual(global_eqn_ids[0]) = q_in;
   for (size_t i = 0; i < num_outlets; i++) {
-    // Get parameters
-    auto resistance = alpha[global_param_ids[i]];
-    auto inductance = alpha[global_param_ids[num_outlets + i]];
-    double stenosis_coeff = 0.0;
-    if (global_param_ids.size() / num_outlets > 2) {
-      stenosis_coeff = alpha[global_param_ids[2 * num_outlets + i]];
-    }
     auto q_out = y[global_var_ids[3 + 2 * i]];
-    auto p_out = y[global_var_ids[2 + 2 * i]];
     auto dq_out = dy[global_var_ids[3 + 2 * i]];
-    auto stenosis_resistance = stenosis_coeff * fabs(q_out);
 
     // Resistance
-    jacobian.coeffRef(global_eqn_ids[i + 1], global_param_ids[i]) = -q_out;
+    jacobian.coeffRef(global_eqn_ids[i + 1], global_param_ids[i]) = q_out;
 
     // Inductance
     jacobian.coeffRef(global_eqn_ids[i + 1],
-                      global_param_ids[num_outlets + i]) = -dq_out;
+                      global_param_ids[num_outlets + i]) = dq_out;
 
     // Stenosis Coefficent
     if (global_param_ids.size() / num_outlets > 2) {
       jacobian.coeffRef(global_eqn_ids[i + 1],
                         global_param_ids[2 * num_outlets + i]) =
-          -fabs(q_out) * q_out;
+          fabs(q_out) * q_out;
     }
-
-    residual(global_eqn_ids[0]) -= q_out;
-    residual(global_eqn_ids[i + 1]) =
-        p_in - p_out - (resistance + stenosis_resistance) * q_out -
-        inductance * dq_out;
   }
 }

src/model/BloodVesselJunction.h

@@ -210,12 +210,10 @@ class BloodVesselJunction : public Block {
    *
    * @param jacobian Jacobian with respect to the parameters
    * @param alpha Current parameter vector
-   * @param residual Residual with respect to the parameters
    * @param y Current solution
    * @param dy Time-derivative of the current solution
    */
   void update_gradient(Eigen::SparseMatrix<double>& jacobian,
-                       Eigen::Matrix<double, Eigen::Dynamic, 1>& residual,
                        Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha,
                        std::vector<double>& y,
                        std::vector<double>& dy) override;

src/model/Junction.cpp

@@ -31,14 +31,3 @@ void Junction::update_constant(SparseSystem& system,
                       global_var_ids[i]) = -1.0;
   }
 }
-
-void Junction::update_gradient(
-    Eigen::SparseMatrix<double>& jacobian,
-    Eigen::Matrix<double, Eigen::Dynamic, 1>& residual,
-    Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha, std::vector<double>& y,
-    std::vector<double>& dy) {
-  // Pressure conservation
-  residual(global_eqn_ids[0]) = y[global_var_ids[0]] - y[global_var_ids[2]];
-
-  residual(global_eqn_ids[1]) = y[global_var_ids[1]] - y[global_var_ids[3]];
-}

src/model/Junction.h

@@ -120,22 +120,6 @@ class Junction : public Block {
   void update_constant(SparseSystem& system,
                        std::vector<double>& parameters) override;
 
-  /**
-   * @brief Set the gradient of the block contributions with respect to the
-   * parameters
-   *
-   * @param jacobian Jacobian with respect to the parameters
-   * @param alpha Current parameter vector
-   * @param residual Residual with respect to the parameters
-   * @param y Current solution
-   * @param dy Time-derivative of the current solution
-   */
-  void update_gradient(Eigen::SparseMatrix<double>& jacobian,
-                       Eigen::Matrix<double, Eigen::Dynamic, 1>& residual,
-                       Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha,
-                       std::vector<double>& y,
-                       std::vector<double>& dy) override;
-
   /**
    * @brief Number of triplets of element
    *

src/optimize/LevenbergMarquardtOptimizer.cpp

@@ -26,6 +26,11 @@ LevenbergMarquardtOptimizer::LevenbergMarquardtOptimizer(
   mat = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>(num_active,
                                                               num_active);
   vec = Eigen::Matrix<double, Eigen::Dynamic, 1>::Zero(num_active);
+
+  // Set up the solver system so that its residual assembly can be reused for a
+  // single observation at a time.
+  system = SparseSystem(num_vars);
+  system.reserve(model);
 }
 
 Eigen::Matrix<double, Eigen::Dynamic, 1> LevenbergMarquardtOptimizer::run(
@@ -68,14 +73,46 @@ void LevenbergMarquardtOptimizer::update_gradient(
   jacobian.setZero();
   residual.setZero();
 
+  // Push the current parameter vector into the model so that the solver's
+  // residual assembly uses it.
+  for (size_t k = 0; k < num_params; k++) {
+    model->update_parameter_value(k, alpha[k]);
+  }
+
+  // Assemble the parameter-dependent constant system contributions (E and F)
+  // once for the current parameters.
+  model->update_constant(system);
+
+  Eigen::Matrix<double, Eigen::Dynamic, 1> y_dpoint(num_vars);
+  Eigen::Matrix<double, Eigen::Dynamic, 1> dy_dpoint(num_vars);
+
   // Assemble gradient and residual
   for (size_t i = 0; i < num_obs; i++) {
+    // Copy the observation for this datapoint into Eigen vectors.
+    for (size_t k = 0; k < num_vars; k++) {
+      y_dpoint[k] = y_obs[i][k];
+      dy_dpoint[k] = dy_obs[i][k];
+    }
+
+    // Reuse the residual of the solver: r = -C - E*ydot - F*y. This is the
+    // same residual the solver assembles, so it does not need to be redefined
+    // here.
+    model->update_solution(system, y_dpoint, dy_dpoint);
+    system.update_residual(y_dpoint, dy_dpoint);
+    residual.segment(num_eqns * i, num_eqns) = system.residual;
+
+    // Assemble the Jacobian of the residual with respect to the parameters.
     for (size_t j = 0; j < model->get_num_blocks(true); j++) {
       auto block = model->get_block(j);
+      // Blocks without calibration parameters do not contribute to the
+      // Jacobian (e.g. a normal junction).
+      if (block->global_param_ids.empty()) {
+        continue;
+      }
       for (size_t l = 0; l < block->global_eqn_ids.size(); l++) {
         block->global_eqn_ids[l] += num_eqns * i;
       }
-      block->update_gradient(jacobian, residual, alpha, y_obs[i], dy_obs[i]);
+      block->update_gradient(jacobian, alpha, y_obs[i], dy_obs[i]);
       for (size_t l = 0; l < block->global_eqn_ids.size(); l++) {
         block->global_eqn_ids[l] -= num_eqns * i;
       }

src/optimize/LevenbergMarquardtOptimizer.h

@@ -11,6 +11,7 @@
 #include <Eigen/Sparse>
 
 #include "Model.h"
+#include "SparseSystem.h"
 
 /**
  * @brief Levenberg-Marquardt optimization class
@@ -29,6 +30,13 @@
  * \mathbb{R}^{NxN}\f$, and system vector \f$\boldsymbol{c} \in
  * \mathbb{R}^{N}\f$.
  *
+ * This residual is identical (up to an overall sign) to the one the solver
+ * assembles in SparseSystem::update_residual. It is therefore reused here
+ * instead of being redefined, and only its Jacobian with respect to the
+ * parameters is assembled by the blocks (see Block::update_gradient). The
+ * overall sign cancels in the normal equations below, so it does not affect
+ * the parameter increment.
+ *
  * The least squares problem can be formulated as
  *
  * \f[
@@ -115,6 +123,7 @@ class LevenbergMarquardtOptimizer {
   Eigen::Matrix<double, Eigen::Dynamic, 1> vec;
   std::vector<int> active_param_ids;
   Model* model;
+  SparseSystem system;  ///< Solver system used to reuse the residual assembly
   double lambda;
 
   int num_obs;

src/optimize/calibrate.cpp

@@ -93,7 +93,12 @@ nlohmann::json calibrate(const nlohmann::json& config) {
     auto* block = model.create_block(block_type);
     int num_slots = num_param_slots(*block, /*stride=*/1);
     std::vector<int> param_ids;
-    for (int k = 0; k < num_slots; k++) param_ids.push_back(param_counter++);
+    // Register a model parameter for each slot so the solver's residual
+    // assembly can read the current parameter values during calibration.
+    for (int k = 0; k < num_slots; k++) {
+      param_ids.push_back(param_counter++);
+      model.add_parameter(0.0);
+    }
     model.add_block(block, vessel_name, param_ids);
     vessel_id_map.insert({vessel_config["vessel_id"], vessel_name});
     DEBUG_MSG("Created vessel " << vessel_name);
@@ -124,7 +129,12 @@ nlohmann::json calibrate(const nlohmann::json& config) {
       auto* block = model.create_block("BloodVesselJunction");
       int num_slots = num_param_slots(*block, num_outlets);
       std::vector<int> param_ids;
-      for (int k = 0; k < num_slots; k++) param_ids.push_back(param_counter++);
+      // Register a model parameter for each slot so the solver's residual
+      // assembly can read the current parameter values during calibration.
+      for (int k = 0; k < num_slots; k++) {
+        param_ids.push_back(param_counter++);
+        model.add_parameter(0.0);
+      }
       model.add_block(block, junction_name, param_ids);
 
       register_active(*block, param_ids,