diff --git a/src/model/Block.cpp b/src/model/Block.cpp index 7800fa10d..2d46cd2a0 100644 --- a/src/model/Block.cpp +++ b/src/model/Block.cpp @@ -62,7 +62,6 @@ void Block::update_solution( void Block::post_solve(Eigen::Matrix& y) {} void Block::update_gradient(Eigen::SparseMatrix& jacobian, - Eigen::Matrix& residual, Eigen::Matrix& alpha, std::vector& y, std::vector& dy) { throw std::runtime_error("Gradient calculation not implemented for block " + diff --git a/src/model/Block.h b/src/model/Block.h index 849aeb956..1cd5d6ae0 100644 --- a/src/model/Block.h +++ b/src/model/Block.h @@ -252,20 +252,22 @@ class Block { virtual void post_solve(Eigen::Matrix& 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& jacobian, - Eigen::Matrix& residual, - Eigen::Matrix& alpha, std::vector& y, - std::vector& dy); + virtual void update_gradient(Eigen::SparseMatrix& jacobian, + Eigen::Matrix& alpha, + std::vector& y, std::vector& dy); /** * @brief Number of triplets of element diff --git a/src/model/BloodVessel.cpp b/src/model/BloodVessel.cpp index faa6f25c1..6b7140b5f 100644 --- a/src/model/BloodVessel.cpp +++ b/src/model/BloodVessel.cpp @@ -58,13 +58,9 @@ void BloodVessel::update_solution( void BloodVessel::update_gradient( Eigen::SparseMatrix& jacobian, - Eigen::Matrix& residual, Eigen::Matrix& alpha, std::vector& y, std::vector& 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]]; @@ -72,7 +68,6 @@ void BloodVessel::update_gradient( 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) { @@ -80,25 +75,20 @@ void BloodVessel::update_gradient( } 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; } diff --git a/src/model/BloodVessel.h b/src/model/BloodVessel.h index 77c500417..cdf50dd91 100644 --- a/src/model/BloodVessel.h +++ b/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& jacobian, - Eigen::Matrix& residual, Eigen::Matrix& alpha, std::vector& y, std::vector& dy) override; diff --git a/src/model/BloodVesselCRL.cpp b/src/model/BloodVesselCRL.cpp index 5c264ca35..bd07735bf 100644 --- a/src/model/BloodVesselCRL.cpp +++ b/src/model/BloodVesselCRL.cpp @@ -50,38 +50,19 @@ void BloodVesselCRL::update_solution( void BloodVesselCRL::update_gradient( Eigen::SparseMatrix& jacobian, - Eigen::Matrix& residual, - Eigen::Matrix& alpha, std::vector& y, + Eigen::Matrix& /*alpha*/, std::vector& y, std::vector& 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; } diff --git a/src/model/BloodVesselCRL.h b/src/model/BloodVesselCRL.h index 9bcb8467b..1a17edeb4 100644 --- a/src/model/BloodVesselCRL.h +++ b/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& jacobian, - Eigen::Matrix& residual, Eigen::Matrix& alpha, std::vector& y, std::vector& dy) override; diff --git a/src/model/BloodVesselJunction.cpp b/src/model/BloodVesselJunction.cpp index 56469e31e..0455ad171 100644 --- a/src/model/BloodVesselJunction.cpp +++ b/src/model/BloodVesselJunction.cpp @@ -54,43 +54,24 @@ void BloodVesselJunction::update_solution( void BloodVesselJunction::update_gradient( Eigen::SparseMatrix& jacobian, - Eigen::Matrix& residual, - Eigen::Matrix& alpha, std::vector& y, + Eigen::Matrix& /*alpha*/, std::vector& y, std::vector& 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; } } diff --git a/src/model/BloodVesselJunction.h b/src/model/BloodVesselJunction.h index c77dfae05..aa3b82bbc 100644 --- a/src/model/BloodVesselJunction.h +++ b/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& jacobian, - Eigen::Matrix& residual, Eigen::Matrix& alpha, std::vector& y, std::vector& dy) override; diff --git a/src/model/Junction.cpp b/src/model/Junction.cpp index 4f773df09..64391e39b 100644 --- a/src/model/Junction.cpp +++ b/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& jacobian, - Eigen::Matrix& residual, - Eigen::Matrix& alpha, std::vector& y, - std::vector& 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]]; -} diff --git a/src/model/Junction.h b/src/model/Junction.h index 185c9bd9b..f9e31cf61 100644 --- a/src/model/Junction.h +++ b/src/model/Junction.h @@ -120,22 +120,6 @@ class Junction : public Block { void update_constant(SparseSystem& system, std::vector& 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& jacobian, - Eigen::Matrix& residual, - Eigen::Matrix& alpha, - std::vector& y, - std::vector& dy) override; - /** * @brief Number of triplets of element * diff --git a/src/optimize/LevenbergMarquardtOptimizer.cpp b/src/optimize/LevenbergMarquardtOptimizer.cpp index a91e584db..e46d6c9f6 100644 --- a/src/optimize/LevenbergMarquardtOptimizer.cpp +++ b/src/optimize/LevenbergMarquardtOptimizer.cpp @@ -26,6 +26,11 @@ LevenbergMarquardtOptimizer::LevenbergMarquardtOptimizer( mat = Eigen::Matrix(num_active, num_active); vec = Eigen::Matrix::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 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 y_dpoint(num_vars); + Eigen::Matrix 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; } diff --git a/src/optimize/LevenbergMarquardtOptimizer.h b/src/optimize/LevenbergMarquardtOptimizer.h index 7d5a93d86..6cafbbc07 100644 --- a/src/optimize/LevenbergMarquardtOptimizer.h +++ b/src/optimize/LevenbergMarquardtOptimizer.h @@ -11,6 +11,7 @@ #include #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 vec; std::vector active_param_ids; Model* model; + SparseSystem system; ///< Solver system used to reuse the residual assembly double lambda; int num_obs; diff --git a/src/optimize/calibrate.cpp b/src/optimize/calibrate.cpp index 0d880c3e3..0dc2b5f4b 100644 --- a/src/optimize/calibrate.cpp +++ b/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 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 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,