From b38a15067b361f86e42fa4e08c41dad39b4e61da Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 10 Feb 2026 19:00:10 +0000 Subject: [PATCH 01/23] Initial plan From 7bf169e2d468c32c9c446e37ee6820e38ec88aa7 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 10 Feb 2026 19:05:15 +0000 Subject: [PATCH 02/23] Add ActivationFunction class and rename PiecewiseCosineChamber to LinearElastanceChamber Co-authored-by: aabrown100-git <71424733+aabrown100-git@users.noreply.github.com> --- src/model/ActivationFunction.h | 312 +++++++++++++++++++++++++++ src/model/BlockType.h | 3 +- src/model/CMakeLists.txt | 3 + src/model/LinearElastanceChamber.cpp | 62 ++++++ src/model/LinearElastanceChamber.h | 200 +++++++++++++++++ src/model/Model.cpp | 3 +- src/model/Model.h | 1 + 7 files changed, 582 insertions(+), 2 deletions(-) create mode 100644 src/model/ActivationFunction.h create mode 100644 src/model/LinearElastanceChamber.cpp create mode 100644 src/model/LinearElastanceChamber.h diff --git a/src/model/ActivationFunction.h b/src/model/ActivationFunction.h new file mode 100644 index 000000000..2646a3af0 --- /dev/null +++ b/src/model/ActivationFunction.h @@ -0,0 +1,312 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +/** + * @file ActivationFunction.h + * @brief Activation function classes for cardiac chamber models + */ + +#ifndef SVZERODSOLVER_MODEL_ACTIVATIONFUNCTION_HPP_ +#define SVZERODSOLVER_MODEL_ACTIVATIONFUNCTION_HPP_ + +#include +#include +#include +#include +#include + +/** + * @brief Enum for activation function types + */ +enum class ActivationType { + HALF_COSINE = 0, + PIECEWISE_COSINE = 1, + TWO_HILL = 2 +}; + +/** + * @brief Base class for activation functions + * + * Activation functions compute the activation value (between 0 and 1) at a + * given time point within a cardiac cycle. These are used to modulate + * chamber elastance over time. + */ +class ActivationFunction { + public: + /** + * @brief Virtual destructor + */ + virtual ~ActivationFunction() = default; + + /** + * @brief Compute activation value at given time + * + * @param time Current time + * @param cardiac_period Period of cardiac cycle + * @return Activation value between 0 and 1 + */ + virtual double compute(double time, double cardiac_period) = 0; + + /** + * @brief Get the type of activation function + * + * @return Type of activation function + */ + virtual ActivationType get_type() const = 0; +}; + +/** + * @brief Half cosine activation function + * + * This implements the activation function used in the original + * ChamberElastanceInductor. The activation follows a half cosine wave + * during the contraction period. + * + * \f[ + * A(t) = \begin{cases} + * -\frac{1}{2}\cos(2\pi t_{contract}/t_{twitch}) + \frac{1}{2}, & \text{if } t_{contract} \le t_{twitch} \\ + * 0, & \text{otherwise} + * \end{cases} + * \f] + * + * where \f$t_{contract} = \max(0, t_{in\_cycle} - t_{active})\f$ + */ +class HalfCosineActivation : public ActivationFunction { + public: + /** + * @brief Construct a new Half Cosine Activation object + * + * @param t_active Time when activation begins within cardiac cycle + * @param t_twitch Duration of the contraction twitch + */ + HalfCosineActivation(double t_active, double t_twitch) + : t_active_(t_active), t_twitch_(t_twitch) {} + + /** + * @brief Compute activation value + * + * @param time Current time + * @param cardiac_period Period of cardiac cycle + * @return Activation value between 0 and 1 + */ + double compute(double time, double cardiac_period) override { + double t_in_cycle = std::fmod(time, cardiac_period); + + double t_contract = 0.0; + if (t_in_cycle >= t_active_) { + t_contract = t_in_cycle - t_active_; + } + + double act = 0.0; + if (t_contract <= t_twitch_) { + act = -0.5 * std::cos(2.0 * M_PI * t_contract / t_twitch_) + 0.5; + } + + return act; + } + + /** + * @brief Get the type of activation function + * + * @return ActivationType::HALF_COSINE + */ + ActivationType get_type() const override { + return ActivationType::HALF_COSINE; + } + + private: + double t_active_; ///< Time when activation begins + double t_twitch_; ///< Duration of contraction twitch +}; + +/** + * @brief Piecewise cosine activation function + * + * This implements the activation function from the PiecewiseCosineChamber + * (Regazzoni chamber model). The activation consists of separate contraction + * and relaxation phases, each following a cosine curve. + * + * \f[ + * \phi(t, t_C, t_R, T_C, T_R) = \begin{cases} + * \frac{1}{2}\left[1 - \cos\left(\frac{\pi}{T_C} \bmod(t - t_C, T_{HB})\right)\right], + * & \text{if } 0 \le \bmod(t - t_C, T_{HB}) < T_C \\ + * \frac{1}{2}\left[1 + \cos\left(\frac{\pi}{T_R} \bmod(t - t_R, T_{HB})\right)\right], + * & \text{if } 0 \le \bmod(t - t_R, T_{HB}) < T_R \\ + * 0, & \text{otherwise} + * \end{cases} + * \f] + */ +class PiecewiseCosineActivation : public ActivationFunction { + public: + /** + * @brief Construct a new Piecewise Cosine Activation object + * + * @param contract_start Time when contraction starts + * @param relax_start Time when relaxation starts + * @param contract_duration Duration of contraction phase + * @param relax_duration Duration of relaxation phase + */ + PiecewiseCosineActivation(double contract_start, double relax_start, + double contract_duration, double relax_duration) + : contract_start_(contract_start), + relax_start_(relax_start), + contract_duration_(contract_duration), + relax_duration_(relax_duration) {} + + /** + * @brief Compute activation value + * + * @param time Current time + * @param cardiac_period Period of cardiac cycle + * @return Activation value between 0 and 1 + */ + double compute(double time, double cardiac_period) override { + double phi = 0.0; + + double piecewise_condition = std::fmod(time - contract_start_, cardiac_period); + + if (0.0 <= piecewise_condition && piecewise_condition < contract_duration_) { + phi = 0.5 * (1.0 - std::cos((M_PI * piecewise_condition) / contract_duration_)); + } else { + piecewise_condition = std::fmod(time - relax_start_, cardiac_period); + if (0.0 <= piecewise_condition && piecewise_condition < relax_duration_) { + phi = 0.5 * (1.0 + std::cos((M_PI * piecewise_condition) / relax_duration_)); + } + } + + return phi; + } + + /** + * @brief Get the type of activation function + * + * @return ActivationType::PIECEWISE_COSINE + */ + ActivationType get_type() const override { + return ActivationType::PIECEWISE_COSINE; + } + + private: + double contract_start_; ///< Start time of contraction + double relax_start_; ///< Start time of relaxation + double contract_duration_; ///< Duration of contraction + double relax_duration_; ///< Duration of relaxation +}; + +/** + * @brief Two hill activation function + * + * This implements the two-hill activation function which provides more + * flexible and physiologically realistic waveforms. See + * https://link.springer.com/article/10.1007/s10439-022-03047-3 + * + * The activation is computed as: + * \f[ + * A(t) = C \cdot \frac{g_1(t)}{1 + g_1(t)} \cdot \frac{1}{1 + g_2(t)} + * \f] + * + * where: + * \f[ + * g_1(t) = \left(\frac{t_{shifted}}{\tau_1}\right)^{m_1}, \quad + * g_2(t) = \left(\frac{t_{shifted}}{\tau_2}\right)^{m_2} + * \f] + * + * and \f$t_{shifted} = (t - t_{shift}) \bmod T_{cardiac}\f$, and \f$C\f$ is a + * normalization constant to ensure max activation is 1. + */ +class TwoHillActivation : public ActivationFunction { + public: + /** + * @brief Construct a new Two Hill Activation object + * + * @param t_shift Time shift parameter + * @param tau_1 Time constant for first hill + * @param tau_2 Time constant for second hill + * @param m1 Exponent for first hill + * @param m2 Exponent for second hill + * @param cardiac_period Cardiac cycle period (needed for normalization) + */ + TwoHillActivation(double t_shift, double tau_1, double tau_2, double m1, + double m2, double cardiac_period) + : t_shift_(t_shift), + tau_1_(tau_1), + tau_2_(tau_2), + m1_(m1), + m2_(m2), + normalization_factor_(1.0), + normalization_initialized_(false) { + initialize_normalization(cardiac_period); + } + + /** + * @brief Compute activation value + * + * @param time Current time + * @param cardiac_period Period of cardiac cycle + * @return Activation value between 0 and 1 + */ + double compute(double time, double cardiac_period) override { + if (!normalization_initialized_) { + throw std::runtime_error("TwoHillActivation: Normalization not initialized"); + } + + double t_in_cycle = std::fmod(time, cardiac_period); + + // Compute shifted time (handle negative modulo) + double t_shifted = std::fmod(t_in_cycle - t_shift_, cardiac_period); + t_shifted = (t_shifted >= 0.0) ? t_shifted : t_shifted + cardiac_period; + + // Compute hill functions + double g1 = std::pow(t_shifted / tau_1_, m1_); + double g2 = std::pow(t_shifted / tau_2_, m2_); + + // Compute activation with normalization + double act = normalization_factor_ * (g1 / (1.0 + g1)) * (1.0 / (1.0 + g2)); + + return act; + } + + /** + * @brief Get the type of activation function + * + * @return ActivationType::TWO_HILL + */ + ActivationType get_type() const override { + return ActivationType::TWO_HILL; + } + + private: + /** + * @brief Initialize normalization factor + * + * Computes the maximum value of the two-hill function over one cardiac + * cycle to normalize the activation to [0, 1]. + * + * @param cardiac_period Period of cardiac cycle + */ + void initialize_normalization(double cardiac_period) { + double max_value = 0.0; + double dt = 1e-5; + + for (double t_temp = 0.0; t_temp < cardiac_period; t_temp += dt) { + double g1 = std::pow(t_temp / tau_1_, m1_); + double g2 = std::pow(t_temp / tau_2_, m2_); + double two_hill_val = (g1 / (1.0 + g1)) * (1.0 / (1.0 + g2)); + + max_value = std::max(max_value, two_hill_val); + } + + normalization_factor_ = 1.0 / max_value; + normalization_initialized_ = true; + } + + double t_shift_; ///< Time shift parameter + double tau_1_; ///< Time constant for first hill + double tau_2_; ///< Time constant for second hill + double m1_; ///< Exponent for first hill + double m2_; ///< Exponent for second hill + double normalization_factor_; ///< Normalization constant + bool normalization_initialized_; ///< Flag for normalization +}; + +#endif // SVZERODSOLVER_MODEL_ACTIVATIONFUNCTION_HPP_ diff --git a/src/model/BlockType.h b/src/model/BlockType.h index a57972647..76d862c96 100644 --- a/src/model/BlockType.h +++ b/src/model/BlockType.h @@ -31,7 +31,8 @@ enum class BlockType { chamber_sphere = 15, blood_vessel_CRL = 16, piecewise_cosine_chamber = 17, - piecewise_valve = 18 + piecewise_valve = 18, + linear_elastance_chamber = 19 }; /** diff --git a/src/model/CMakeLists.txt b/src/model/CMakeLists.txt index a76fc23e3..62e1fe414 100644 --- a/src/model/CMakeLists.txt +++ b/src/model/CMakeLists.txt @@ -30,10 +30,12 @@ set(CXXSRCS ValveTanh.cpp WindkesselBC.cpp PiecewiseCosineChamber.cpp + LinearElastanceChamber.cpp PiecewiseValve.cpp ) set(HDRS + ActivationFunction.h Block.h BlockType.h BloodVessel.h @@ -59,6 +61,7 @@ set(HDRS ValveTanh.h WindkesselBC.h PiecewiseCosineChamber.h + LinearElastanceChamber.h PiecewiseValve.h ) diff --git a/src/model/LinearElastanceChamber.cpp b/src/model/LinearElastanceChamber.cpp new file mode 100644 index 000000000..b91dc6546 --- /dev/null +++ b/src/model/LinearElastanceChamber.cpp @@ -0,0 +1,62 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +#include "LinearElastanceChamber.h" + +void LinearElastanceChamber::setup_dofs(DOFHandler& dofhandler) { + // Internal variable is chamber volume + Block::setup_dofs_(dofhandler, 3, {"Vc"}); +} + +void LinearElastanceChamber::update_constant(SparseSystem& system, + std::vector& parameters) { + // Eq 0: P_in - E(t)(Vc - Vrest) = 0 + system.F.coeffRef(global_eqn_ids[0], global_var_ids[0]) = 1.0; + + // Eq 1: P_in - P_out = 0 + system.F.coeffRef(global_eqn_ids[1], global_var_ids[0]) = 1.0; + system.F.coeffRef(global_eqn_ids[1], global_var_ids[2]) = -1.0; + + // Eq 2: Q_in - Q_out - dVc = 0 + system.F.coeffRef(global_eqn_ids[2], global_var_ids[1]) = 1.0; + system.F.coeffRef(global_eqn_ids[2], global_var_ids[3]) = -1.0; + system.E.coeffRef(global_eqn_ids[2], global_var_ids[4]) = -1.0; +} + +void LinearElastanceChamber::update_time(SparseSystem& system, + std::vector& parameters) { + get_elastance_values(parameters); + + // Eq 0: P_in - E(t)(Vc - Vrest) = P_in - E(t)*Vc + E(t)*Vrest = 0 + system.F.coeffRef(global_eqn_ids[0], global_var_ids[4]) = -1 * Elas; + system.C.coeffRef(global_eqn_ids[0]) = + Elas * parameters[global_param_ids[ParamId::VREST]]; +} + +void LinearElastanceChamber::get_elastance_values( + std::vector& parameters) { + double Emax = parameters[global_param_ids[ParamId::EMAX]]; + double Epass = parameters[global_param_ids[ParamId::EPASS]]; + double Vrest = parameters[global_param_ids[ParamId::VREST]]; + double contract_start = parameters[global_param_ids[ParamId::CSTART]]; + double relax_start = parameters[global_param_ids[ParamId::RSTART]]; + double contract_duration = parameters[global_param_ids[ParamId::CDUR]]; + double relax_duration = parameters[global_param_ids[ParamId::RDUR]]; + + auto T_HB = model->cardiac_cycle_period; + + double phi = 0; + + auto piecewise_condition = fmod(model->time - contract_start, T_HB); + + if (0 <= piecewise_condition && piecewise_condition < contract_duration) { + phi = 0.5 * (1 - cos((M_PI * piecewise_condition) / contract_duration)); + } else { + piecewise_condition = fmod(model->time - relax_start, T_HB); + if (0 <= piecewise_condition && piecewise_condition < relax_duration) { + phi = 0.5 * (1 + cos((M_PI * piecewise_condition) / relax_duration)); + } + } + + Elas = Epass + Emax * phi; +} diff --git a/src/model/LinearElastanceChamber.h b/src/model/LinearElastanceChamber.h new file mode 100644 index 000000000..e668817eb --- /dev/null +++ b/src/model/LinearElastanceChamber.h @@ -0,0 +1,200 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +/** + * @file LinearElastanceChamber.h + * @brief model::LinearElastanceChamber source file + */ + +#ifndef SVZERODSOLVER_MODEL_LINEARELASTANCECHAMBER_HPP_ +#define SVZERODSOLVER_MODEL_LINEARELASTANCECHAMBER_HPP_ + +#include + +#include "Block.h" +#include "Model.h" +#include "SparseSystem.h" +#include "debug.h" + +/** + * @brief Cardiac chamber with linear elastance (no inductor). + * + * Models a cardiac chamber as a time-varying capacitor (elastance with + * specified resting volumes) without inductance. See \cite Regazzoni2022 + * + * This chamber block can be connected to other blocks using junctions. + * + * \f[ + * \begin{circuitikz} + * \draw node[left] {$Q_{in}$} [-latex] (0,0) -- (0.8,0); + * \draw (1,0) node[anchor=south]{$P_{in}$} + * to[short, *-*] (3,0) node[anchor=south]{$P_{out}$} + * (1,0) to [vC, l=$E$, *-] (1,-1.5) + * node[ground]{}; + * \draw (3.2,0) -- (4.0,0) node[right] {$Q_{out}$} [-latex]; + * \end{circuitikz} + * \f] + * + * ### Governing equations + * + * \f[ + * P_{in}-E(t)(V_c-V_{rest})=0 + * \f] + * + * \f[ + * P_{in}-P_{out}=0 + * \f] + * + * \f[ + * Q_{in}-Q_{out}-\dot{V}_c=0 + * \f] + * + * ### Local contributions + * + * \f[ + * \mathbf{y}^{e}=\left[\begin{array}{lllll}P_{in} & Q_{in} & + * P_{out} & Q_{out} & V_c\end{array}\right]^{T} \f] + * + * \f[ + * \mathbf{E}^{e}=\left[\begin{array}{ccccc} + * 0 & 0 & 0 & 0 & 0\\ + * 0 & 0 & 0 & 0 & 0\\ + * 0 & 0 & 0 & 0 & -1 + * \end{array}\right] + * \f] + * + * \f[ + * \mathbf{F}^{e}=\left[\begin{array}{ccccc} + * 1 & 0 & 0 & 0 & E(t) \\ + * 1 & 0 & -1 & 0 & 0 \\ + * 0 & 1 & 0 & -1 & 0 + * \end{array}\right] + * \f] + * + * \f[ + * \mathbf{c}^{e}=\left[\begin{array}{c} + * E(t)V_{rest} \\ + * 0 \\ + * 0 + * \end{array}\right] + * \f] + * + * In the above equations, + * + * \f[ + * E_i(t) = E_i^{\text{pass}} + E_i^{\text{act,max}} \, + * \phi\!\left(t, t_C^i, t_R^i, T_C^i, T_R^i\right), + * \f] + * + * \f[ + * \phi(t, t_C, t_R, T_C, T_R) = + * \begin{cases} + * \frac{1}{2}\left[1 - \cos\!\left(\frac{\pi}{T_C} \operatorname{mod}(t - t_C, + * T_{\mathrm{HB}})\right)\right], & \text{if } 0 \le \operatorname{mod}(t - + * t_C, T_{\mathrm{HB}}) < T_C, \\[1.2em] \frac{1}{2}\left[1 + + * \cos\!\left(\frac{\pi}{T_R} \operatorname{mod}(t - t_R, + * T_{\mathrm{HB}})\right)\right], & \text{if } 0 \le \operatorname{mod}(t - + * t_R, T_{\mathrm{HB}}) < T_R, \\[1.2em] 0, & \text{otherwise.} \end{cases} \f] + * + * ### Parameters + * + * Parameter sequence for constructing this block + * + * * `0` Emax: Maximum elastance + * * `1` Epass: Passive elastance + * * `2` Vrest: Rest diastolic volume + * * `3` contract_start: Contract start time + * * `4` relax_start: Relax start time + * * `5` contract_duration: Contract duration + * * `6` relax_duration: Relaxation duration + * + * ### Internal variables + * + * Names of internal variables in this block's output: + * + * * `Vc`: Chamber volume + * + */ +class LinearElastanceChamber : public Block { + public: + /** + * @brief Construct a new LinearElastanceChamber object + * + * @param id Global ID of the block + * @param model The model to which the block belongs + */ + LinearElastanceChamber(int id, Model* model) + : Block(id, model, BlockType::linear_elastance_chamber, + BlockClass::chamber, + {{"Emax", InputParameter()}, + {"Epass", InputParameter()}, + {"Vrest", InputParameter()}, + {"contract_start", InputParameter()}, + {"relax_start", InputParameter()}, + {"contract_duration", InputParameter()}, + {"relax_duration", InputParameter()}}) {} + + /** + * @brief Local IDs of the parameters + * + */ + enum ParamId { + EMAX = 0, + EPASS = 1, + VREST = 2, + CSTART = 3, + RSTART = 4, + CDUR = 5, + RDUR = 6 + }; + + /** + * @brief Set up the degrees of freedom (DOF) of the block + * + * Set global_var_ids and global_eqn_ids of the element based on the + * number of equations and the number of internal variables of the + * element. + * + * @param dofhandler Degree-of-freedom handler to register variables and + * equations at + */ + void setup_dofs(DOFHandler& dofhandler); + + /** + * @brief Update the constant contributions of the element in a sparse + system + * + * @param system System to update contributions at + * @param parameters Parameters of the model + */ + void update_constant(SparseSystem& system, std::vector& parameters); + + /** + * @brief Update the time-dependent contributions of the element in a sparse + * system + * + * @param system System to update contributions at + * @param parameters Parameters of the model + */ + void update_time(SparseSystem& system, std::vector& parameters); + + /** + * @brief Number of triplets of element + * + * Number of triplets that the element contributes to the global system + * (relevant for sparse memory reservation) + */ + TripletsContributions num_triplets{6, 2, 0}; + + private: + double Elas; // Chamber Elastance + + /** + * @brief Update the elastance functions which depend on time + * + * @param parameters Parameters of the model + */ + void get_elastance_values(std::vector& parameters); +}; + +#endif // SVZERODSOLVER_MODEL_LINEARELASTANCECHAMBER_HPP_ diff --git a/src/model/Model.cpp b/src/model/Model.cpp index a527161a5..f4d2143eb 100644 --- a/src/model/Model.cpp +++ b/src/model/Model.cpp @@ -31,7 +31,8 @@ Model::Model() { {"ChamberElastanceInductor", block_factory()}, {"BloodVesselCRL", block_factory()}, {"PiecewiseValve", block_factory()}, - {"PiecewiseCosineChamber", block_factory()}}; + {"PiecewiseCosineChamber", block_factory()}, + {"LinearElastanceChamber", block_factory()}}; } Model::~Model() {} diff --git a/src/model/Model.h b/src/model/Model.h index 6a444e93a..23b673e58 100644 --- a/src/model/Model.h +++ b/src/model/Model.h @@ -29,6 +29,7 @@ #include "DOFHandler.h" #include "FlowReferenceBC.h" #include "Junction.h" +#include "LinearElastanceChamber.h" #include "Node.h" #include "OpenLoopCoronaryBC.h" #include "Parameter.h" From 980cec7c2cad03d347c32e07a87b12528d32414c Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 10 Feb 2026 19:07:32 +0000 Subject: [PATCH 03/23] Refactor ChamberElastanceInductor and LinearElastanceChamber to use ActivationFunction Co-authored-by: aabrown100-git <71424733+aabrown100-git@users.noreply.github.com> --- src/model/ChamberElastanceInductor.cpp | 66 +++++++++++++++++++----- src/model/ChamberElastanceInductor.h | 38 ++++++++++++-- src/model/LinearElastanceChamber.cpp | 70 ++++++++++++++++++++------ src/model/LinearElastanceChamber.h | 36 +++++++++++-- 4 files changed, 173 insertions(+), 37 deletions(-) diff --git a/src/model/ChamberElastanceInductor.cpp b/src/model/ChamberElastanceInductor.cpp index 440835e8a..143d502c8 100644 --- a/src/model/ChamberElastanceInductor.cpp +++ b/src/model/ChamberElastanceInductor.cpp @@ -9,6 +9,11 @@ void ChamberElastanceInductor::setup_dofs(DOFHandler& dofhandler) { void ChamberElastanceInductor::update_constant( SparseSystem& system, std::vector& parameters) { + // Initialize activation function on first call + if (!activation_func_) { + initialize_activation_function(parameters); + } + double L = parameters[global_param_ids[ParamId::IMPEDANCE]]; // Eq 0: P_in - E(t)(Vc - Vrest) = 0 @@ -40,22 +45,59 @@ void ChamberElastanceInductor::get_elastance_values( double Emin = parameters[global_param_ids[ParamId::EMIN]]; double Vrd = parameters[global_param_ids[ParamId::VRD]]; double Vrs = parameters[global_param_ids[ParamId::VRS]]; - double t_active = parameters[global_param_ids[ParamId::TACTIVE]]; - double t_twitch = parameters[global_param_ids[ParamId::TTWITCH]]; auto T_cardiac = model->cardiac_cycle_period; - auto t_in_cycle = fmod(model->time, T_cardiac); + + // Compute activation using the activation function + double act = activation_func_->compute(model->time, T_cardiac); - double t_contract = 0; - if (t_in_cycle >= t_active) { - t_contract = t_in_cycle - t_active; - } + Vrest = (1.0 - act) * (Vrd - Vrs) + Vrs; + Elas = (Emax - Emin) * act + Emin; +} - double act = 0; - if (t_contract <= t_twitch) { - act = -0.5 * cos(2 * M_PI * t_contract / t_twitch) + 0.5; +void ChamberElastanceInductor::initialize_activation_function( + std::vector& parameters) { + // Check if activation_type parameter is provided (optional parameter) + int activation_type_int = 0; // Default to HALF_COSINE + if (global_param_ids.count(ParamId::ACTIVATION_TYPE) > 0) { + activation_type_int = static_cast( + parameters[global_param_ids[ParamId::ACTIVATION_TYPE]]); } - Vrest = (1.0 - act) * (Vrd - Vrs) + Vrs; - Elas = (Emax - Emin) * act + Emin; + auto T_cardiac = model->cardiac_cycle_period; + + switch (activation_type_int) { + case 0: // HALF_COSINE (default) + { + double t_active = parameters[global_param_ids[ParamId::TACTIVE]]; + double t_twitch = parameters[global_param_ids[ParamId::TTWITCH]]; + activation_func_ = std::make_unique(t_active, t_twitch); + break; + } + case 1: // PIECEWISE_COSINE + { + double contract_start = parameters[global_param_ids[ParamId::CSTART]]; + double relax_start = parameters[global_param_ids[ParamId::RSTART]]; + double contract_duration = parameters[global_param_ids[ParamId::CDUR]]; + double relax_duration = parameters[global_param_ids[ParamId::RDUR]]; + activation_func_ = std::make_unique( + contract_start, relax_start, contract_duration, relax_duration); + break; + } + case 2: // TWO_HILL + { + double t_shift = parameters[global_param_ids[ParamId::TSHIFT]]; + double tau_1 = parameters[global_param_ids[ParamId::TAU_1]]; + double tau_2 = parameters[global_param_ids[ParamId::TAU_2]]; + double m1 = parameters[global_param_ids[ParamId::M1]]; + double m2 = parameters[global_param_ids[ParamId::M2]]; + activation_func_ = std::make_unique( + t_shift, tau_1, tau_2, m1, m2, T_cardiac); + break; + } + default: + throw std::runtime_error( + "ChamberElastanceInductor: Invalid activation_type " + + std::to_string(activation_type_int)); + } } diff --git a/src/model/ChamberElastanceInductor.h b/src/model/ChamberElastanceInductor.h index 527947cc7..2476afe78 100644 --- a/src/model/ChamberElastanceInductor.h +++ b/src/model/ChamberElastanceInductor.h @@ -8,7 +8,9 @@ #define SVZERODSOLVER_MODEL_CHAMBERELASTANCEINDUCTOR_HPP_ #include +#include +#include "ActivationFunction.h" #include "Block.h" #include "Model.h" #include "SparseSystem.h" @@ -152,9 +154,19 @@ class ChamberElastanceInductor : public Block { {"Emin", InputParameter()}, {"Vrd", InputParameter()}, {"Vrs", InputParameter()}, - {"t_active", InputParameter()}, - {"t_twitch", InputParameter()}, - {"Impedance", InputParameter()}}) {} + {"t_active", InputParameter(false)}, + {"t_twitch", InputParameter(false)}, + {"Impedance", InputParameter()}, + {"activation_type", InputParameter(false)}, + {"contract_start", InputParameter(false)}, + {"relax_start", InputParameter(false)}, + {"contract_duration", InputParameter(false)}, + {"relax_duration", InputParameter(false)}, + {"t_shift", InputParameter(false)}, + {"tau_1", InputParameter(false)}, + {"tau_2", InputParameter(false)}, + {"m1", InputParameter(false)}, + {"m2", InputParameter(false)}}) {} /** * @brief Local IDs of the parameters @@ -167,7 +179,17 @@ class ChamberElastanceInductor : public Block { VRS = 3, TACTIVE = 4, TTWITCH = 5, - IMPEDANCE = 6 + IMPEDANCE = 6, + ACTIVATION_TYPE = 7, + CSTART = 8, + RSTART = 9, + CDUR = 10, + RDUR = 11, + TSHIFT = 12, + TAU_1 = 13, + TAU_2 = 14, + M1 = 15, + M2 = 16 }; /** @@ -211,6 +233,14 @@ class ChamberElastanceInductor : public Block { private: double Elas; // Chamber Elastance double Vrest; // Rest Volume + std::unique_ptr activation_func_; // Activation function + + /** + * @brief Initialize the activation function based on parameters + * + * @param parameters Parameters of the model + */ + void initialize_activation_function(std::vector& parameters); /** * @brief Update the elastance functions which depend on time diff --git a/src/model/LinearElastanceChamber.cpp b/src/model/LinearElastanceChamber.cpp index b91dc6546..464e75332 100644 --- a/src/model/LinearElastanceChamber.cpp +++ b/src/model/LinearElastanceChamber.cpp @@ -10,6 +10,11 @@ void LinearElastanceChamber::setup_dofs(DOFHandler& dofhandler) { void LinearElastanceChamber::update_constant(SparseSystem& system, std::vector& parameters) { + // Initialize activation function on first call + if (!activation_func_) { + initialize_activation_function(parameters); + } + // Eq 0: P_in - E(t)(Vc - Vrest) = 0 system.F.coeffRef(global_eqn_ids[0], global_var_ids[0]) = 1.0; @@ -37,26 +42,59 @@ void LinearElastanceChamber::get_elastance_values( std::vector& parameters) { double Emax = parameters[global_param_ids[ParamId::EMAX]]; double Epass = parameters[global_param_ids[ParamId::EPASS]]; - double Vrest = parameters[global_param_ids[ParamId::VREST]]; - double contract_start = parameters[global_param_ids[ParamId::CSTART]]; - double relax_start = parameters[global_param_ids[ParamId::RSTART]]; - double contract_duration = parameters[global_param_ids[ParamId::CDUR]]; - double relax_duration = parameters[global_param_ids[ParamId::RDUR]]; - auto T_HB = model->cardiac_cycle_period; + auto T_cardiac = model->cardiac_cycle_period; + + // Compute activation using the activation function + double phi = activation_func_->compute(model->time, T_cardiac); - double phi = 0; + Elas = Epass + Emax * phi; +} - auto piecewise_condition = fmod(model->time - contract_start, T_HB); +void LinearElastanceChamber::initialize_activation_function( + std::vector& parameters) { + // Check if activation_type parameter is provided (optional parameter) + // Default to PIECEWISE_COSINE (1) for backward compatibility + int activation_type_int = 1; + if (global_param_ids.count(ParamId::ACTIVATION_TYPE) > 0) { + activation_type_int = static_cast( + parameters[global_param_ids[ParamId::ACTIVATION_TYPE]]); + } - if (0 <= piecewise_condition && piecewise_condition < contract_duration) { - phi = 0.5 * (1 - cos((M_PI * piecewise_condition) / contract_duration)); - } else { - piecewise_condition = fmod(model->time - relax_start, T_HB); - if (0 <= piecewise_condition && piecewise_condition < relax_duration) { - phi = 0.5 * (1 + cos((M_PI * piecewise_condition) / relax_duration)); + auto T_cardiac = model->cardiac_cycle_period; + + switch (activation_type_int) { + case 0: // HALF_COSINE + { + double t_active = parameters[global_param_ids[ParamId::TACTIVE]]; + double t_twitch = parameters[global_param_ids[ParamId::TTWITCH]]; + activation_func_ = std::make_unique(t_active, t_twitch); + break; + } + case 1: // PIECEWISE_COSINE (default for backward compatibility) + { + double contract_start = parameters[global_param_ids[ParamId::CSTART]]; + double relax_start = parameters[global_param_ids[ParamId::RSTART]]; + double contract_duration = parameters[global_param_ids[ParamId::CDUR]]; + double relax_duration = parameters[global_param_ids[ParamId::RDUR]]; + activation_func_ = std::make_unique( + contract_start, relax_start, contract_duration, relax_duration); + break; + } + case 2: // TWO_HILL + { + double t_shift = parameters[global_param_ids[ParamId::TSHIFT]]; + double tau_1 = parameters[global_param_ids[ParamId::TAU_1]]; + double tau_2 = parameters[global_param_ids[ParamId::TAU_2]]; + double m1 = parameters[global_param_ids[ParamId::M1]]; + double m2 = parameters[global_param_ids[ParamId::M2]]; + activation_func_ = std::make_unique( + t_shift, tau_1, tau_2, m1, m2, T_cardiac); + break; } + default: + throw std::runtime_error( + "LinearElastanceChamber: Invalid activation_type " + + std::to_string(activation_type_int)); } - - Elas = Epass + Emax * phi; } diff --git a/src/model/LinearElastanceChamber.h b/src/model/LinearElastanceChamber.h index e668817eb..2c9b78d15 100644 --- a/src/model/LinearElastanceChamber.h +++ b/src/model/LinearElastanceChamber.h @@ -10,7 +10,9 @@ #define SVZERODSOLVER_MODEL_LINEARELASTANCECHAMBER_HPP_ #include +#include +#include "ActivationFunction.h" #include "Block.h" #include "Model.h" #include "SparseSystem.h" @@ -129,10 +131,18 @@ class LinearElastanceChamber : public Block { {{"Emax", InputParameter()}, {"Epass", InputParameter()}, {"Vrest", InputParameter()}, - {"contract_start", InputParameter()}, - {"relax_start", InputParameter()}, - {"contract_duration", InputParameter()}, - {"relax_duration", InputParameter()}}) {} + {"contract_start", InputParameter(false)}, + {"relax_start", InputParameter(false)}, + {"contract_duration", InputParameter(false)}, + {"relax_duration", InputParameter(false)}, + {"activation_type", InputParameter(false)}, + {"t_active", InputParameter(false)}, + {"t_twitch", InputParameter(false)}, + {"t_shift", InputParameter(false)}, + {"tau_1", InputParameter(false)}, + {"tau_2", InputParameter(false)}, + {"m1", InputParameter(false)}, + {"m2", InputParameter(false)}}) {} /** * @brief Local IDs of the parameters @@ -145,7 +155,15 @@ class LinearElastanceChamber : public Block { CSTART = 3, RSTART = 4, CDUR = 5, - RDUR = 6 + RDUR = 6, + ACTIVATION_TYPE = 7, + TACTIVE = 8, + TTWITCH = 9, + TSHIFT = 10, + TAU_1 = 11, + TAU_2 = 12, + M1 = 13, + M2 = 14 }; /** @@ -188,6 +206,14 @@ class LinearElastanceChamber : public Block { private: double Elas; // Chamber Elastance + std::unique_ptr activation_func_; // Activation function + + /** + * @brief Initialize the activation function based on parameters + * + * @param parameters Parameters of the model + */ + void initialize_activation_function(std::vector& parameters); /** * @brief Update the elastance functions which depend on time From 2f7c593d1333e3b9f5fde9f5e72629a9aec6af04 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 10 Feb 2026 19:08:33 +0000 Subject: [PATCH 04/23] Add test cases for two hill activation function Co-authored-by: aabrown100-git <71424733+aabrown100-git@users.noreply.github.com> --- tests/cases/chamber_two_hill_activation.json | 117 +++++++++++++++ .../linear_chamber_two_hill_activation.json | 138 ++++++++++++++++++ 2 files changed, 255 insertions(+) create mode 100644 tests/cases/chamber_two_hill_activation.json create mode 100644 tests/cases/linear_chamber_two_hill_activation.json diff --git a/tests/cases/chamber_two_hill_activation.json b/tests/cases/chamber_two_hill_activation.json new file mode 100644 index 000000000..e8968f3d9 --- /dev/null +++ b/tests/cases/chamber_two_hill_activation.json @@ -0,0 +1,117 @@ +{ + "description": { + "description of test case": "pulsatile pressure -> valve -> chamber (two hill activation) -> valve -> RC -> R", + "analytical results": [ + "Boundary conditions:", + "inlet:", + "pressure = sinusoid from 0 to 10 over 0.17s, then 0 for the rest of cycle", + "outlet:", + "Resistance: Rd = 1000, Pd = 64" + ] + }, + "simulation_parameters": { + "number_of_cardiac_cycles": 30, + "number_of_time_pts_per_cardiac_cycle": 101, + "absolute_tolerance": 1e-9, + "output_all_cycles": false, + "output_variable_based": true + }, + "boundary_conditions": [ + { + "bc_name": "INLET", + "bc_type": "PRESSURE", + "bc_values": { + "P": [ + 0.0, 0.003414699, 0.4075651, 1.431949, 2.938218, 4.722943, 6.545085, + 8.158555, 9.345445, 9.945458, 9.87756, 9.15092, 7.863676, 6.189676, + 4.355004, 2.607442, 1.183009, 0.274081, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0 + ], + "t": [ + 0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1, 0.11, + 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2, 0.21, 0.22, 0.23, + 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.3, 0.31, 0.32, 0.33, 0.34, 0.35, + 0.36, 0.37, 0.38, 0.39, 0.4, 0.41, 0.42, 0.43, 0.44, 0.45, 0.46, 0.47, + 0.48, 0.49, 0.5, 0.51, 0.52, 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, + 0.6, 0.61, 0.62, 0.63, 0.64, 0.65, 0.66, 0.67, 0.68, 0.69, 0.7, 0.71, + 0.72, 0.73, 0.74, 0.75, 0.76, 0.77, 0.78, 0.79, 0.8, 0.81, 0.82, 0.83, + 0.84, 0.85, 0.86, 0.87, 0.88, 0.89, 0.9, 0.91, 0.92, 0.93, 0.94, 0.95, + 0.96, 0.97, 0.98, 0.99, 1.0 + ] + } + }, + { + "bc_name": "OUTLET", + "bc_type": "RESISTANCE", + "bc_values": { + "Pd": 64.0, + "R": 1000.0 + } + } + ], + "junctions": [], + "vessels": [ + { + "boundary_conditions": { + "outlet": "OUTLET" + }, + "vessel_id": 0, + "vessel_length": 10.0, + "vessel_name": "vessel", + "zero_d_element_type": "BloodVessel", + "zero_d_element_values": { + "C": 0.0001, + "R_poiseuille": 100.0 + } + } + ], + "valves": [ + { + "type": "ValveTanh", + "name": "valve0", + "params": { + "Rmax": 100000.0, + "Rmin": 100.0, + "Steepness": 100.0, + "upstream_block": "INLET", + "downstream_block": "ventricle" + } + }, + { + "type": "ValveTanh", + "name": "valve1", + "params": { + "Rmax": 100000.0, + "Rmin": 100.0, + "Steepness": 100.0, + "upstream_block": "ventricle", + "downstream_block": "vessel" + } + } + ], + "chambers": [ + { + "type": "ChamberElastanceInductor", + "name": "ventricle", + "values": { + "Emax": 1.057, + "Emin": 0.091, + "Vrd": 26.1, + "Vrs": 18.0, + "Impedance": 0.000351787, + "activation_type": 2, + "t_shift": 0.15, + "tau_1": 0.25, + "tau_2": 0.35, + "m1": 1.5, + "m2": 10.0 + } + } + ], + "initial_condition": { + "Vc:ventricle": 96.07 + } +} diff --git a/tests/cases/linear_chamber_two_hill_activation.json b/tests/cases/linear_chamber_two_hill_activation.json new file mode 100644 index 000000000..70d241acf --- /dev/null +++ b/tests/cases/linear_chamber_two_hill_activation.json @@ -0,0 +1,138 @@ +{ + "description": { + "description of test case": "LinearElastanceChamber with two hill activation function", + "analytical results": [ + "Testing LinearElastanceChamber (formerly PiecewiseCosineChamber) with two hill activation" + ] + }, + "simulation_parameters": { + "number_of_cardiac_cycles": 30, + "number_of_time_pts_per_cardiac_cycle": 689, + "output_variable_based": true, + "output_all_cycles": false, + "cardiac_period": 0.68900516753 + }, + "vessels": [ + { + "vessel_id": 1, + "vessel_length": 10.0, + "vessel_name": "pul_artery", + "zero_d_element_type": "BloodVessel", + "zero_d_element_values": { + "C": 0.0, + "L": 0.0, + "R_poiseuille": 0.0 + } + }, + { + "vessel_id": 5, + "vessel_length": 10.0, + "vessel_name": "sys_artery", + "zero_d_element_type": "BloodVessel", + "zero_d_element_values": { + "C": 0.001138164597527442, + "R_poiseuille": 596.82, + "L": 6.665 + } + }, + { + "vessel_id": 7, + "vessel_length": 10.0, + "vessel_name": "sys_vein", + "zero_d_element_type": "BloodVessel", + "zero_d_element_values": { + "C": 0.045011252813952737, + "R_poiseuille": 249.42, + "L": 0.6665 + } + } + ], + "junctions": [ + { + "inlet_blocks": ["sys_vein"], + "junction_name": "J1", + "junction_type": "NORMAL_JUNCTION", + "outlet_blocks": ["right_atrium"] + }, + { + "inlet_vessels": [5], + "junction_name": "J2", + "junction_type": "NORMAL_JUNCTION", + "outlet_vessels": [7] + } + ], + "boundary_conditions": [], + "chambers": [ + { + "type": "LinearElastanceChamber", + "name": "right_atrium", + "values": { + "Emax": 199.95, + "Epass": 89.80375933024839, + "Vrest": 41.680015938842274, + "activation_type": 2, + "t_shift": 0.025, + "tau_1": 0.06, + "tau_2": 0.15, + "m1": 2.0, + "m2": 8.0 + } + }, + { + "type": "LinearElastanceChamber", + "name": "right_ventricle", + "values": { + "Emax": 1662.0158240056528, + "Epass": 40.85565535747109, + "Vrest": 72.05452710344869, + "activation_type": 2, + "t_shift": 0.207, + "tau_1": 0.08, + "tau_2": 0.22, + "m1": 2.5, + "m2": 12.0 + } + } + ], + "valves": [ + { + "type": "PiecewiseValve", + "name": "tricuspid", + "params": { + "Rmax": 6665, + "Rmin": 6.665, + "upstream_block": "right_atrium", + "downstream_block": "right_ventricle" + } + }, + { + "type": "PiecewiseValve", + "name": "pulmonary", + "params": { + "Rmax": 6665, + "Rmin": 6.665, + "upstream_block": "right_ventricle", + "downstream_block": "pul_artery" + } + }, + { + "type": "PiecewiseValve", + "name": "aortic", + "params": { + "Rmax": 6665, + "Rmin": 6.665, + "upstream_block": "right_ventricle", + "downstream_block": "sys_artery" + } + } + ], + "initial_condition": { + "Vc:right_ventricle": 128.58981029386334, + "Vc:right_atrium": 76.8340776729488, + "pressure:aortic:sys_artery": 84604.18388331511, + "pressure:J2:sys_vein": 31311.64989129829, + "pressure:pulmonary:pul_artery": 20525.08438550143, + "flow:sys_artery:J2": 91.00177508885831, + "flow:sys_vein:J1": 112.86832799795421 + } +} From 88e0ad126e711aced3a62522e0abfe3a4965c5c5 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 10 Feb 2026 19:10:01 +0000 Subject: [PATCH 05/23] Address code review comments: improve readability and documentation Co-authored-by: aabrown100-git <71424733+aabrown100-git@users.noreply.github.com> --- src/model/ActivationFunction.h | 7 +++++-- src/model/LinearElastanceChamber.cpp | 2 +- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/model/ActivationFunction.h b/src/model/ActivationFunction.h index 2646a3af0..a554f371d 100644 --- a/src/model/ActivationFunction.h +++ b/src/model/ActivationFunction.h @@ -285,10 +285,13 @@ class TwoHillActivation : public ActivationFunction { * @param cardiac_period Period of cardiac cycle */ void initialize_normalization(double cardiac_period) { + // Time step for numerical integration when finding maximum activation + // Value chosen to balance accuracy and computational cost + constexpr double NORMALIZATION_DT = 1e-5; + double max_value = 0.0; - double dt = 1e-5; - for (double t_temp = 0.0; t_temp < cardiac_period; t_temp += dt) { + for (double t_temp = 0.0; t_temp < cardiac_period; t_temp += NORMALIZATION_DT) { double g1 = std::pow(t_temp / tau_1_, m1_); double g2 = std::pow(t_temp / tau_2_, m2_); double two_hill_val = (g1 / (1.0 + g1)) * (1.0 / (1.0 + g2)); diff --git a/src/model/LinearElastanceChamber.cpp b/src/model/LinearElastanceChamber.cpp index 464e75332..942b8174d 100644 --- a/src/model/LinearElastanceChamber.cpp +++ b/src/model/LinearElastanceChamber.cpp @@ -33,7 +33,7 @@ void LinearElastanceChamber::update_time(SparseSystem& system, get_elastance_values(parameters); // Eq 0: P_in - E(t)(Vc - Vrest) = P_in - E(t)*Vc + E(t)*Vrest = 0 - system.F.coeffRef(global_eqn_ids[0], global_var_ids[4]) = -1 * Elas; + system.F.coeffRef(global_eqn_ids[0], global_var_ids[4]) = -Elas; system.C.coeffRef(global_eqn_ids[0]) = Elas * parameters[global_param_ids[ParamId::VREST]]; } From 0b2ecf036b5ef6dfef4c6603cf12da93f8d2070a Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 10 Feb 2026 19:11:06 +0000 Subject: [PATCH 06/23] Add documentation for activation functions Co-authored-by: aabrown100-git <71424733+aabrown100-git@users.noreply.github.com> --- _codeql_detected_source_root | 1 + docs/pages/activation_functions.md | 153 +++++++++++++++++++++++++++++ 2 files changed, 154 insertions(+) create mode 120000 _codeql_detected_source_root create mode 100644 docs/pages/activation_functions.md diff --git a/_codeql_detected_source_root b/_codeql_detected_source_root new file mode 120000 index 000000000..945c9b46d --- /dev/null +++ b/_codeql_detected_source_root @@ -0,0 +1 @@ +. \ No newline at end of file diff --git a/docs/pages/activation_functions.md b/docs/pages/activation_functions.md new file mode 100644 index 000000000..2d91eff11 --- /dev/null +++ b/docs/pages/activation_functions.md @@ -0,0 +1,153 @@ +# Activation Functions for Chamber Models + +## Overview + +Chamber models in svZeroDSolver support different activation functions to model the time-varying elastance of cardiac chambers. The activation function determines how the chamber contracts and relaxes over the cardiac cycle. + +## Available Activation Functions + +### 1. Half Cosine Activation (Type 0) + +The default activation function for `ChamberElastanceInductor`. It uses a half cosine wave to model contraction: + +``` +A(t) = -0.5 * cos(2π * t_contract / t_twitch) + 0.5 if t_contract ≤ t_twitch +A(t) = 0 otherwise +``` + +where `t_contract = max(0, t - t_active)` + +**Parameters:** +- `t_active`: Time when activation begins within the cardiac cycle +- `t_twitch`: Duration of the contraction twitch + +**Example usage in JSON:** +```json +{ + "type": "ChamberElastanceInductor", + "name": "ventricle", + "values": { + "Emax": 1.057, + "Emin": 0.091, + "Vrd": 26.1, + "Vrs": 18.0, + "t_active": 0.2, + "t_twitch": 0.3, + "Impedance": 0.000351787 + } +} +``` + +### 2. Piecewise Cosine Activation (Type 1) + +The default activation function for `LinearElastanceChamber` (formerly `PiecewiseCosineChamber`). It models separate contraction and relaxation phases: + +``` +φ(t) = 0.5 * [1 - cos(π * (t - t_C) / T_C)] during contraction +φ(t) = 0.5 * [1 + cos(π * (t - t_R) / T_R)] during relaxation +φ(t) = 0 otherwise +``` + +**Parameters:** +- `contract_start`: Time when contraction starts +- `relax_start`: Time when relaxation starts +- `contract_duration`: Duration of contraction phase +- `relax_duration`: Duration of relaxation phase + +**Example usage in JSON:** +```json +{ + "type": "LinearElastanceChamber", + "name": "left_atrium", + "values": { + "Emax": 199.95, + "Epass": 260.59, + "Vrest": 26.24, + "contract_start": 0.025, + "relax_start": 0.08625, + "contract_duration": 0.06125, + "relax_duration": 0.18375 + } +} +``` + +### 3. Two Hill Activation (Type 2) + +A more flexible and physiologically realistic activation function based on the two-hill model. This allows for more precise control over the activation waveform shape. + +The activation is computed as: +``` +A(t) = C * [g₁(t) / (1 + g₁(t))] * [1 / (1 + g₂(t))] +``` + +where: +``` +g₁(t) = (t_shifted / τ₁)^m₁ +g₂(t) = (t_shifted / τ₂)^m₂ +t_shifted = (t - t_shift) mod T_cardiac +``` + +C is a normalization constant ensuring maximum activation is 1. + +**Parameters:** +- `t_shift`: Time shift parameter to align activation with cardiac cycle +- `tau_1`: Time constant for first hill (controls rise) +- `tau_2`: Time constant for second hill (controls fall) +- `m1`: Exponent for first hill (controls steepness of rise) +- `m2`: Exponent for second hill (controls steepness of fall) + +**Example usage in JSON:** +```json +{ + "type": "ChamberElastanceInductor", + "name": "ventricle", + "values": { + "Emax": 1.057, + "Emin": 0.091, + "Vrd": 26.1, + "Vrs": 18.0, + "Impedance": 0.000351787, + "activation_type": 2, + "t_shift": 0.15, + "tau_1": 0.25, + "tau_2": 0.35, + "m1": 1.5, + "m2": 10.0 + } +} +``` + +## Switching Between Activation Functions + +To use a specific activation function, add the `activation_type` parameter to your chamber configuration: + +- `activation_type: 0` - Half Cosine (default for ChamberElastanceInductor) +- `activation_type: 1` - Piecewise Cosine (default for LinearElastanceChamber) +- `activation_type: 2` - Two Hill + +## Elastance Calculation + +For all activation functions, the time-varying elastance E(t) is computed as: + +**For ChamberElastanceInductor:** +``` +E(t) = (Emax - Emin) * A(t) + Emin +Vrest(t) = (1 - A(t)) * (Vrd - Vrs) + Vrs +``` + +**For LinearElastanceChamber:** +``` +E(t) = Epass + Emax * φ(t) +``` + +## References + +The two-hill activation function is described in: +- Kaiser, A. D., et al. (2022). "A design-based model of the aortic valve for fluid-structure interaction." Biomechanics and Modeling in Mechanobiology. https://link.springer.com/article/10.1007/s10439-022-03047-3 + +## Backward Compatibility + +Existing configuration files will continue to work without modification: +- `ChamberElastanceInductor` defaults to half cosine activation (type 0) +- `LinearElastanceChamber` defaults to piecewise cosine activation (type 1) +- `PiecewiseCosineChamber` is maintained as an alias to `LinearElastanceChamber` for backward compatibility From b68ffd5d2d68bb6395a99a68231bbfc2debd3113 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 10 Feb 2026 19:25:22 +0000 Subject: [PATCH 07/23] Support string activation_function_type and nested activation_function_values in JSON Co-authored-by: aabrown100-git <71424733+aabrown100-git@users.noreply.github.com> --- docs/pages/activation_functions.md | 84 +++++++++++++++---- src/solve/SimulationParameters.cpp | 43 +++++++++- tests/cases/chamber_two_hill_activation.json | 14 ++-- .../linear_chamber_two_hill_activation.json | 28 ++++--- 4 files changed, 133 insertions(+), 36 deletions(-) diff --git a/docs/pages/activation_functions.md b/docs/pages/activation_functions.md index 2d91eff11..e9b18208d 100644 --- a/docs/pages/activation_functions.md +++ b/docs/pages/activation_functions.md @@ -38,6 +38,26 @@ where `t_contract = max(0, t - t_active)` } ``` +Or with explicit activation function specification: +```json +{ + "type": "ChamberElastanceInductor", + "name": "ventricle", + "values": { + "Emax": 1.057, + "Emin": 0.091, + "Vrd": 26.1, + "Vrs": 18.0, + "Impedance": 0.000351787, + "activation_function_type": "half_cosine", + "activation_function_values": { + "t_active": 0.2, + "t_twitch": 0.3 + } + } +} +``` + ### 2. Piecewise Cosine Activation (Type 1) The default activation function for `LinearElastanceChamber` (formerly `PiecewiseCosineChamber`). It models separate contraction and relaxation phases: @@ -63,10 +83,13 @@ The default activation function for `LinearElastanceChamber` (formerly `Piecewis "Emax": 199.95, "Epass": 260.59, "Vrest": 26.24, - "contract_start": 0.025, - "relax_start": 0.08625, - "contract_duration": 0.06125, - "relax_duration": 0.18375 + "activation_function_type": "piecewise_cosine", + "activation_function_values": { + "contract_start": 0.025, + "relax_start": 0.08625, + "contract_duration": 0.06125, + "relax_duration": 0.18375 + } } } ``` @@ -107,23 +130,49 @@ C is a normalization constant ensuring maximum activation is 1. "Vrd": 26.1, "Vrs": 18.0, "Impedance": 0.000351787, - "activation_type": 2, - "t_shift": 0.15, - "tau_1": 0.25, - "tau_2": 0.35, - "m1": 1.5, - "m2": 10.0 + "activation_function_type": "two_hill", + "activation_function_values": { + "t_shift": 0.15, + "tau_1": 0.25, + "tau_2": 0.35, + "m1": 1.5, + "m2": 10.0 + } } } ``` ## Switching Between Activation Functions -To use a specific activation function, add the `activation_type` parameter to your chamber configuration: +To use a specific activation function, add the `activation_function_type` parameter to your chamber configuration: + +- `activation_function_type: "half_cosine"` - Half Cosine (default for ChamberElastanceInductor) +- `activation_function_type: "piecewise_cosine"` - Piecewise Cosine (default for LinearElastanceChamber) +- `activation_function_type: "two_hill"` - Two Hill + +Group the activation function-specific parameters under `activation_function_values` as a nested dictionary. + +## Backward Compatibility + +For backward compatibility, the old flat format is also supported: + +```json +{ + "type": "ChamberElastanceInductor", + "name": "ventricle", + "values": { + "Emax": 1.057, + "Emin": 0.091, + "Vrd": 26.1, + "Vrs": 18.0, + "t_active": 0.2, + "t_twitch": 0.3, + "Impedance": 0.000351787 + } +} +``` -- `activation_type: 0` - Half Cosine (default for ChamberElastanceInductor) -- `activation_type: 1` - Piecewise Cosine (default for LinearElastanceChamber) -- `activation_type: 2` - Two Hill +This format will continue to work, with the activation function determined by which parameters are provided. ## Elastance Calculation @@ -148,6 +197,7 @@ The two-hill activation function is described in: ## Backward Compatibility Existing configuration files will continue to work without modification: -- `ChamberElastanceInductor` defaults to half cosine activation (type 0) -- `LinearElastanceChamber` defaults to piecewise cosine activation (type 1) -- `PiecewiseCosineChamber` is maintained as an alias to `LinearElastanceChamber` for backward compatibility +- `ChamberElastanceInductor` defaults to half cosine activation when `t_active` and `t_twitch` are provided +- `LinearElastanceChamber` defaults to piecewise cosine activation when contraction/relaxation parameters are provided +- The old numeric `activation_type` parameter (0, 1, 2) is still supported but deprecated in favor of the string format +- Parameters can be specified at the top level without using `activation_function_values` for backward compatibility diff --git a/src/solve/SimulationParameters.cpp b/src/solve/SimulationParameters.cpp index b914c401d..9a80579a3 100644 --- a/src/solve/SimulationParameters.cpp +++ b/src/solve/SimulationParameters.cpp @@ -544,7 +544,48 @@ void create_chambers( const auto& chamber_config = JsonWrapper(config, component, "name", i); std::string chamber_type = chamber_config["type"]; std::string chamber_name = chamber_config["name"]; - generate_block(model, chamber_config["values"], chamber_type, chamber_name); + + // Create a mutable copy of the values to handle activation function parameters + nlohmann::json chamber_values = chamber_config["values"]; + + // Check if activation function type is specified as a string + if (chamber_values.contains("activation_function_type") && + chamber_values["activation_function_type"].is_string()) { + std::string act_type_str = chamber_values["activation_function_type"]; + + // Map string to numeric value + int act_type_num = 0; // Default to half_cosine + if (act_type_str == "half_cosine") { + act_type_num = 0; + } else if (act_type_str == "piecewise_cosine") { + act_type_num = 1; + } else if (act_type_str == "two_hill") { + act_type_num = 2; + } else { + throw std::runtime_error("Unknown activation_function_type: " + act_type_str + + " in chamber " + chamber_name); + } + + // Replace string with numeric value + chamber_values["activation_type"] = act_type_num; + chamber_values.erase("activation_function_type"); + } + + // Check if activation function values are specified in nested dict + if (chamber_values.contains("activation_function_values") && + chamber_values["activation_function_values"].is_object()) { + const auto& act_values = chamber_values["activation_function_values"]; + + // Flatten the nested dictionary - move all parameters to top level + for (auto& [key, value] : act_values.items()) { + chamber_values[key] = value; + } + + // Remove the nested dict + chamber_values.erase("activation_function_values"); + } + + generate_block(model, chamber_values, chamber_type, chamber_name); DEBUG_MSG("Created chamber " << chamber_name); } } diff --git a/tests/cases/chamber_two_hill_activation.json b/tests/cases/chamber_two_hill_activation.json index e8968f3d9..d4c2fbb54 100644 --- a/tests/cases/chamber_two_hill_activation.json +++ b/tests/cases/chamber_two_hill_activation.json @@ -102,12 +102,14 @@ "Vrd": 26.1, "Vrs": 18.0, "Impedance": 0.000351787, - "activation_type": 2, - "t_shift": 0.15, - "tau_1": 0.25, - "tau_2": 0.35, - "m1": 1.5, - "m2": 10.0 + "activation_function_type": "two_hill", + "activation_function_values": { + "t_shift": 0.15, + "tau_1": 0.25, + "tau_2": 0.35, + "m1": 1.5, + "m2": 10.0 + } } } ], diff --git a/tests/cases/linear_chamber_two_hill_activation.json b/tests/cases/linear_chamber_two_hill_activation.json index 70d241acf..f2a03469e 100644 --- a/tests/cases/linear_chamber_two_hill_activation.json +++ b/tests/cases/linear_chamber_two_hill_activation.json @@ -70,12 +70,14 @@ "Emax": 199.95, "Epass": 89.80375933024839, "Vrest": 41.680015938842274, - "activation_type": 2, - "t_shift": 0.025, - "tau_1": 0.06, - "tau_2": 0.15, - "m1": 2.0, - "m2": 8.0 + "activation_function_type": "two_hill", + "activation_function_values": { + "t_shift": 0.025, + "tau_1": 0.06, + "tau_2": 0.15, + "m1": 2.0, + "m2": 8.0 + } } }, { @@ -85,12 +87,14 @@ "Emax": 1662.0158240056528, "Epass": 40.85565535747109, "Vrest": 72.05452710344869, - "activation_type": 2, - "t_shift": 0.207, - "tau_1": 0.08, - "tau_2": 0.22, - "m1": 2.5, - "m2": 12.0 + "activation_function_type": "two_hill", + "activation_function_values": { + "t_shift": 0.207, + "tau_1": 0.08, + "tau_2": 0.22, + "m1": 2.5, + "m2": 12.0 + } } } ], From 2a7a9d271c75f23b006204215d73a48134f1dd4e Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 10 Feb 2026 20:13:17 +0000 Subject: [PATCH 08/23] Refactor initialize_activation_function to use factory method in ActivationFunction Co-authored-by: aabrown100-git <71424733+aabrown100-git@users.noreply.github.com> --- src/model/ActivationFunction.h | 61 ++++++++++++++++++++++++++ src/model/ChamberElastanceInductor.cpp | 55 ++++++++++++----------- src/model/LinearElastanceChamber.cpp | 55 ++++++++++++----------- 3 files changed, 115 insertions(+), 56 deletions(-) diff --git a/src/model/ActivationFunction.h b/src/model/ActivationFunction.h index a554f371d..9d435ccdf 100644 --- a/src/model/ActivationFunction.h +++ b/src/model/ActivationFunction.h @@ -24,6 +24,31 @@ enum class ActivationType { TWO_HILL = 2 }; +/** + * @brief Parameters for creating activation functions + */ +struct ActivationFunctionParams { + ActivationType type; + double cardiac_period; + + // Half cosine parameters + double t_active = 0.0; + double t_twitch = 0.0; + + // Piecewise cosine parameters + double contract_start = 0.0; + double relax_start = 0.0; + double contract_duration = 0.0; + double relax_duration = 0.0; + + // Two hill parameters + double t_shift = 0.0; + double tau_1 = 0.0; + double tau_2 = 0.0; + double m1 = 0.0; + double m2 = 0.0; +}; + /** * @brief Base class for activation functions * @@ -53,6 +78,15 @@ class ActivationFunction { * @return Type of activation function */ virtual ActivationType get_type() const = 0; + + /** + * @brief Factory method to create activation function from parameters + * + * @param params Parameters for creating the activation function + * @return Unique pointer to the created activation function + */ + static std::unique_ptr create( + const ActivationFunctionParams& params); }; /** @@ -312,4 +346,31 @@ class TwoHillActivation : public ActivationFunction { bool normalization_initialized_; ///< Flag for normalization }; +/** + * @brief Factory method implementation for creating activation functions + */ +inline std::unique_ptr ActivationFunction::create( + const ActivationFunctionParams& params) { + switch (params.type) { + case ActivationType::HALF_COSINE: + return std::make_unique( + params.t_active, params.t_twitch); + + case ActivationType::PIECEWISE_COSINE: + return std::make_unique( + params.contract_start, params.relax_start, + params.contract_duration, params.relax_duration); + + case ActivationType::TWO_HILL: + return std::make_unique( + params.t_shift, params.tau_1, params.tau_2, + params.m1, params.m2, params.cardiac_period); + + default: + throw std::runtime_error( + "ActivationFunction::create: Invalid activation type " + + std::to_string(static_cast(params.type))); + } +} + #endif // SVZERODSOLVER_MODEL_ACTIVATIONFUNCTION_HPP_ diff --git a/src/model/ChamberElastanceInductor.cpp b/src/model/ChamberElastanceInductor.cpp index 143d502c8..8fce01733 100644 --- a/src/model/ChamberElastanceInductor.cpp +++ b/src/model/ChamberElastanceInductor.cpp @@ -64,40 +64,39 @@ void ChamberElastanceInductor::initialize_activation_function( parameters[global_param_ids[ParamId::ACTIVATION_TYPE]]); } - auto T_cardiac = model->cardiac_cycle_period; - - switch (activation_type_int) { - case 0: // HALF_COSINE (default) - { - double t_active = parameters[global_param_ids[ParamId::TACTIVE]]; - double t_twitch = parameters[global_param_ids[ParamId::TTWITCH]]; - activation_func_ = std::make_unique(t_active, t_twitch); + // Prepare parameters for factory method + ActivationFunctionParams params; + params.type = static_cast(activation_type_int); + params.cardiac_period = model->cardiac_cycle_period; + + // Set parameters based on type + switch (params.type) { + case ActivationType::HALF_COSINE: + params.t_active = parameters[global_param_ids[ParamId::TACTIVE]]; + params.t_twitch = parameters[global_param_ids[ParamId::TTWITCH]]; break; - } - case 1: // PIECEWISE_COSINE - { - double contract_start = parameters[global_param_ids[ParamId::CSTART]]; - double relax_start = parameters[global_param_ids[ParamId::RSTART]]; - double contract_duration = parameters[global_param_ids[ParamId::CDUR]]; - double relax_duration = parameters[global_param_ids[ParamId::RDUR]]; - activation_func_ = std::make_unique( - contract_start, relax_start, contract_duration, relax_duration); + + case ActivationType::PIECEWISE_COSINE: + params.contract_start = parameters[global_param_ids[ParamId::CSTART]]; + params.relax_start = parameters[global_param_ids[ParamId::RSTART]]; + params.contract_duration = parameters[global_param_ids[ParamId::CDUR]]; + params.relax_duration = parameters[global_param_ids[ParamId::RDUR]]; break; - } - case 2: // TWO_HILL - { - double t_shift = parameters[global_param_ids[ParamId::TSHIFT]]; - double tau_1 = parameters[global_param_ids[ParamId::TAU_1]]; - double tau_2 = parameters[global_param_ids[ParamId::TAU_2]]; - double m1 = parameters[global_param_ids[ParamId::M1]]; - double m2 = parameters[global_param_ids[ParamId::M2]]; - activation_func_ = std::make_unique( - t_shift, tau_1, tau_2, m1, m2, T_cardiac); + + case ActivationType::TWO_HILL: + params.t_shift = parameters[global_param_ids[ParamId::TSHIFT]]; + params.tau_1 = parameters[global_param_ids[ParamId::TAU_1]]; + params.tau_2 = parameters[global_param_ids[ParamId::TAU_2]]; + params.m1 = parameters[global_param_ids[ParamId::M1]]; + params.m2 = parameters[global_param_ids[ParamId::M2]]; break; - } + default: throw std::runtime_error( "ChamberElastanceInductor: Invalid activation_type " + std::to_string(activation_type_int)); } + + // Use factory method to create activation function + activation_func_ = ActivationFunction::create(params); } diff --git a/src/model/LinearElastanceChamber.cpp b/src/model/LinearElastanceChamber.cpp index 942b8174d..0353d0ebe 100644 --- a/src/model/LinearElastanceChamber.cpp +++ b/src/model/LinearElastanceChamber.cpp @@ -61,40 +61,39 @@ void LinearElastanceChamber::initialize_activation_function( parameters[global_param_ids[ParamId::ACTIVATION_TYPE]]); } - auto T_cardiac = model->cardiac_cycle_period; - - switch (activation_type_int) { - case 0: // HALF_COSINE - { - double t_active = parameters[global_param_ids[ParamId::TACTIVE]]; - double t_twitch = parameters[global_param_ids[ParamId::TTWITCH]]; - activation_func_ = std::make_unique(t_active, t_twitch); + // Prepare parameters for factory method + ActivationFunctionParams params; + params.type = static_cast(activation_type_int); + params.cardiac_period = model->cardiac_cycle_period; + + // Set parameters based on type + switch (params.type) { + case ActivationType::HALF_COSINE: + params.t_active = parameters[global_param_ids[ParamId::TACTIVE]]; + params.t_twitch = parameters[global_param_ids[ParamId::TTWITCH]]; break; - } - case 1: // PIECEWISE_COSINE (default for backward compatibility) - { - double contract_start = parameters[global_param_ids[ParamId::CSTART]]; - double relax_start = parameters[global_param_ids[ParamId::RSTART]]; - double contract_duration = parameters[global_param_ids[ParamId::CDUR]]; - double relax_duration = parameters[global_param_ids[ParamId::RDUR]]; - activation_func_ = std::make_unique( - contract_start, relax_start, contract_duration, relax_duration); + + case ActivationType::PIECEWISE_COSINE: + params.contract_start = parameters[global_param_ids[ParamId::CSTART]]; + params.relax_start = parameters[global_param_ids[ParamId::RSTART]]; + params.contract_duration = parameters[global_param_ids[ParamId::CDUR]]; + params.relax_duration = parameters[global_param_ids[ParamId::RDUR]]; break; - } - case 2: // TWO_HILL - { - double t_shift = parameters[global_param_ids[ParamId::TSHIFT]]; - double tau_1 = parameters[global_param_ids[ParamId::TAU_1]]; - double tau_2 = parameters[global_param_ids[ParamId::TAU_2]]; - double m1 = parameters[global_param_ids[ParamId::M1]]; - double m2 = parameters[global_param_ids[ParamId::M2]]; - activation_func_ = std::make_unique( - t_shift, tau_1, tau_2, m1, m2, T_cardiac); + + case ActivationType::TWO_HILL: + params.t_shift = parameters[global_param_ids[ParamId::TSHIFT]]; + params.tau_1 = parameters[global_param_ids[ParamId::TAU_1]]; + params.tau_2 = parameters[global_param_ids[ParamId::TAU_2]]; + params.m1 = parameters[global_param_ids[ParamId::M1]]; + params.m2 = parameters[global_param_ids[ParamId::M2]]; break; - } + default: throw std::runtime_error( "LinearElastanceChamber: Invalid activation_type " + std::to_string(activation_type_int)); } + + // Use factory method to create activation function + activation_func_ = ActivationFunction::create(params); } From 6c1835771d2cf2f00f67bf5f97ba8cbb04095f80 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 10 Feb 2026 20:38:10 +0000 Subject: [PATCH 09/23] Make activation function parameters members of ActivationFunction instances Co-authored-by: aabrown100-git <71424733+aabrown100-git@users.noreply.github.com> --- src/model/ChamberElastanceInductor.cpp | 44 +++++---------- src/model/ChamberElastanceInductor.h | 34 +++++------- src/model/LinearElastanceChamber.cpp | 50 +++++++---------- src/model/LinearElastanceChamber.h | 30 +++++----- src/solve/SimulationParameters.cpp | 77 ++++++++++++++++++++++---- src/solve/SimulationParameters.h | 3 + 6 files changed, 128 insertions(+), 110 deletions(-) diff --git a/src/model/ChamberElastanceInductor.cpp b/src/model/ChamberElastanceInductor.cpp index 8fce01733..5888caf4b 100644 --- a/src/model/ChamberElastanceInductor.cpp +++ b/src/model/ChamberElastanceInductor.cpp @@ -62,41 +62,23 @@ void ChamberElastanceInductor::initialize_activation_function( if (global_param_ids.count(ParamId::ACTIVATION_TYPE) > 0) { activation_type_int = static_cast( parameters[global_param_ids[ParamId::ACTIVATION_TYPE]]); + activation_params_.type = static_cast(activation_type_int); } - // Prepare parameters for factory method - ActivationFunctionParams params; - params.type = static_cast(activation_type_int); - params.cardiac_period = model->cardiac_cycle_period; + // Set cardiac period + activation_params_.cardiac_period = model->cardiac_cycle_period; - // Set parameters based on type - switch (params.type) { - case ActivationType::HALF_COSINE: - params.t_active = parameters[global_param_ids[ParamId::TACTIVE]]; - params.t_twitch = parameters[global_param_ids[ParamId::TTWITCH]]; - break; - - case ActivationType::PIECEWISE_COSINE: - params.contract_start = parameters[global_param_ids[ParamId::CSTART]]; - params.relax_start = parameters[global_param_ids[ParamId::RSTART]]; - params.contract_duration = parameters[global_param_ids[ParamId::CDUR]]; - params.relax_duration = parameters[global_param_ids[ParamId::RDUR]]; - break; - - case ActivationType::TWO_HILL: - params.t_shift = parameters[global_param_ids[ParamId::TSHIFT]]; - params.tau_1 = parameters[global_param_ids[ParamId::TAU_1]]; - params.tau_2 = parameters[global_param_ids[ParamId::TAU_2]]; - params.m1 = parameters[global_param_ids[ParamId::M1]]; - params.m2 = parameters[global_param_ids[ParamId::M2]]; - break; - - default: - throw std::runtime_error( - "ChamberElastanceInductor: Invalid activation_type " + - std::to_string(activation_type_int)); + // For backward compatibility, if activation parameters weren't set from JSON, + // try to get them from the parameter vector (for old flat format) + if (activation_params_.type == ActivationType::HALF_COSINE) { + if (global_param_ids.count(ParamId::TACTIVE) > 0) { + activation_params_.t_active = parameters[global_param_ids[ParamId::TACTIVE]]; + } + if (global_param_ids.count(ParamId::TTWITCH) > 0) { + activation_params_.t_twitch = parameters[global_param_ids[ParamId::TTWITCH]]; + } } // Use factory method to create activation function - activation_func_ = ActivationFunction::create(params); + activation_func_ = ActivationFunction::create(activation_params_); } diff --git a/src/model/ChamberElastanceInductor.h b/src/model/ChamberElastanceInductor.h index 2476afe78..dd8077264 100644 --- a/src/model/ChamberElastanceInductor.h +++ b/src/model/ChamberElastanceInductor.h @@ -157,16 +157,7 @@ class ChamberElastanceInductor : public Block { {"t_active", InputParameter(false)}, {"t_twitch", InputParameter(false)}, {"Impedance", InputParameter()}, - {"activation_type", InputParameter(false)}, - {"contract_start", InputParameter(false)}, - {"relax_start", InputParameter(false)}, - {"contract_duration", InputParameter(false)}, - {"relax_duration", InputParameter(false)}, - {"t_shift", InputParameter(false)}, - {"tau_1", InputParameter(false)}, - {"tau_2", InputParameter(false)}, - {"m1", InputParameter(false)}, - {"m2", InputParameter(false)}}) {} + {"activation_type", InputParameter(false)}}) {} /** * @brief Local IDs of the parameters @@ -180,16 +171,7 @@ class ChamberElastanceInductor : public Block { TACTIVE = 4, TTWITCH = 5, IMPEDANCE = 6, - ACTIVATION_TYPE = 7, - CSTART = 8, - RSTART = 9, - CDUR = 10, - RDUR = 11, - TSHIFT = 12, - TAU_1 = 13, - TAU_2 = 14, - M1 = 15, - M2 = 16 + ACTIVATION_TYPE = 7 }; /** @@ -234,7 +216,19 @@ class ChamberElastanceInductor : public Block { double Elas; // Chamber Elastance double Vrest; // Rest Volume std::unique_ptr activation_func_; // Activation function + ActivationFunctionParams activation_params_; // Activation function parameters + public: + /** + * @brief Set activation function parameters from JSON configuration + * + * @param params Activation function parameters + */ + void set_activation_params(const ActivationFunctionParams& params) { + activation_params_ = params; + } + + private: /** * @brief Initialize the activation function based on parameters * diff --git a/src/model/LinearElastanceChamber.cpp b/src/model/LinearElastanceChamber.cpp index 0353d0ebe..60f5ca065 100644 --- a/src/model/LinearElastanceChamber.cpp +++ b/src/model/LinearElastanceChamber.cpp @@ -59,41 +59,29 @@ void LinearElastanceChamber::initialize_activation_function( if (global_param_ids.count(ParamId::ACTIVATION_TYPE) > 0) { activation_type_int = static_cast( parameters[global_param_ids[ParamId::ACTIVATION_TYPE]]); + activation_params_.type = static_cast(activation_type_int); } - // Prepare parameters for factory method - ActivationFunctionParams params; - params.type = static_cast(activation_type_int); - params.cardiac_period = model->cardiac_cycle_period; + // Set cardiac period + activation_params_.cardiac_period = model->cardiac_cycle_period; - // Set parameters based on type - switch (params.type) { - case ActivationType::HALF_COSINE: - params.t_active = parameters[global_param_ids[ParamId::TACTIVE]]; - params.t_twitch = parameters[global_param_ids[ParamId::TTWITCH]]; - break; - - case ActivationType::PIECEWISE_COSINE: - params.contract_start = parameters[global_param_ids[ParamId::CSTART]]; - params.relax_start = parameters[global_param_ids[ParamId::RSTART]]; - params.contract_duration = parameters[global_param_ids[ParamId::CDUR]]; - params.relax_duration = parameters[global_param_ids[ParamId::RDUR]]; - break; - - case ActivationType::TWO_HILL: - params.t_shift = parameters[global_param_ids[ParamId::TSHIFT]]; - params.tau_1 = parameters[global_param_ids[ParamId::TAU_1]]; - params.tau_2 = parameters[global_param_ids[ParamId::TAU_2]]; - params.m1 = parameters[global_param_ids[ParamId::M1]]; - params.m2 = parameters[global_param_ids[ParamId::M2]]; - break; - - default: - throw std::runtime_error( - "LinearElastanceChamber: Invalid activation_type " + - std::to_string(activation_type_int)); + // For backward compatibility, if activation parameters weren't set from JSON, + // try to get them from the parameter vector (for old flat format) + if (activation_params_.type == ActivationType::PIECEWISE_COSINE) { + if (global_param_ids.count(ParamId::CSTART) > 0) { + activation_params_.contract_start = parameters[global_param_ids[ParamId::CSTART]]; + } + if (global_param_ids.count(ParamId::RSTART) > 0) { + activation_params_.relax_start = parameters[global_param_ids[ParamId::RSTART]]; + } + if (global_param_ids.count(ParamId::CDUR) > 0) { + activation_params_.contract_duration = parameters[global_param_ids[ParamId::CDUR]]; + } + if (global_param_ids.count(ParamId::RDUR) > 0) { + activation_params_.relax_duration = parameters[global_param_ids[ParamId::RDUR]]; + } } // Use factory method to create activation function - activation_func_ = ActivationFunction::create(params); + activation_func_ = ActivationFunction::create(activation_params_); } diff --git a/src/model/LinearElastanceChamber.h b/src/model/LinearElastanceChamber.h index 2c9b78d15..8b744abcb 100644 --- a/src/model/LinearElastanceChamber.h +++ b/src/model/LinearElastanceChamber.h @@ -135,14 +135,7 @@ class LinearElastanceChamber : public Block { {"relax_start", InputParameter(false)}, {"contract_duration", InputParameter(false)}, {"relax_duration", InputParameter(false)}, - {"activation_type", InputParameter(false)}, - {"t_active", InputParameter(false)}, - {"t_twitch", InputParameter(false)}, - {"t_shift", InputParameter(false)}, - {"tau_1", InputParameter(false)}, - {"tau_2", InputParameter(false)}, - {"m1", InputParameter(false)}, - {"m2", InputParameter(false)}}) {} + {"activation_type", InputParameter(false)}}) {} /** * @brief Local IDs of the parameters @@ -156,14 +149,7 @@ class LinearElastanceChamber : public Block { RSTART = 4, CDUR = 5, RDUR = 6, - ACTIVATION_TYPE = 7, - TACTIVE = 8, - TTWITCH = 9, - TSHIFT = 10, - TAU_1 = 11, - TAU_2 = 12, - M1 = 13, - M2 = 14 + ACTIVATION_TYPE = 7 }; /** @@ -207,7 +193,19 @@ class LinearElastanceChamber : public Block { private: double Elas; // Chamber Elastance std::unique_ptr activation_func_; // Activation function + ActivationFunctionParams activation_params_; // Activation function parameters + public: + /** + * @brief Set activation function parameters from JSON configuration + * + * @param params Activation function parameters + */ + void set_activation_params(const ActivationFunctionParams& params) { + activation_params_ = params; + } + + private: /** * @brief Initialize the activation function based on parameters * diff --git a/src/solve/SimulationParameters.cpp b/src/solve/SimulationParameters.cpp index 9a80579a3..a1045f604 100644 --- a/src/solve/SimulationParameters.cpp +++ b/src/solve/SimulationParameters.cpp @@ -548,27 +548,35 @@ void create_chambers( // Create a mutable copy of the values to handle activation function parameters nlohmann::json chamber_values = chamber_config["values"]; + // Prepare activation function parameters + ActivationFunctionParams act_params; + act_params.type = ActivationType::HALF_COSINE; // Default + act_params.cardiac_period = 0.0; // Will be set later from model + // Check if activation function type is specified as a string if (chamber_values.contains("activation_function_type") && chamber_values["activation_function_type"].is_string()) { std::string act_type_str = chamber_values["activation_function_type"]; - // Map string to numeric value - int act_type_num = 0; // Default to half_cosine + // Map string to activation type if (act_type_str == "half_cosine") { - act_type_num = 0; + act_params.type = ActivationType::HALF_COSINE; } else if (act_type_str == "piecewise_cosine") { - act_type_num = 1; + act_params.type = ActivationType::PIECEWISE_COSINE; } else if (act_type_str == "two_hill") { - act_type_num = 2; + act_params.type = ActivationType::TWO_HILL; } else { throw std::runtime_error("Unknown activation_function_type: " + act_type_str + " in chamber " + chamber_name); } - // Replace string with numeric value - chamber_values["activation_type"] = act_type_num; + // Replace string with numeric value for backward compatibility + chamber_values["activation_type"] = static_cast(act_params.type); chamber_values.erase("activation_function_type"); + } else if (chamber_values.contains("activation_type")) { + // Handle numeric activation_type for backward compatibility + int act_type_num = chamber_values["activation_type"]; + act_params.type = static_cast(act_type_num); } // Check if activation function values are specified in nested dict @@ -576,16 +584,61 @@ void create_chambers( chamber_values["activation_function_values"].is_object()) { const auto& act_values = chamber_values["activation_function_values"]; - // Flatten the nested dictionary - move all parameters to top level - for (auto& [key, value] : act_values.items()) { - chamber_values[key] = value; + // Extract activation function parameters based on type + if (act_params.type == ActivationType::HALF_COSINE) { + if (act_values.contains("t_active")) act_params.t_active = act_values["t_active"]; + if (act_values.contains("t_twitch")) act_params.t_twitch = act_values["t_twitch"]; + } else if (act_params.type == ActivationType::PIECEWISE_COSINE) { + if (act_values.contains("contract_start")) act_params.contract_start = act_values["contract_start"]; + if (act_values.contains("relax_start")) act_params.relax_start = act_values["relax_start"]; + if (act_values.contains("contract_duration")) act_params.contract_duration = act_values["contract_duration"]; + if (act_values.contains("relax_duration")) act_params.relax_duration = act_values["relax_duration"]; + } else if (act_params.type == ActivationType::TWO_HILL) { + if (act_values.contains("t_shift")) act_params.t_shift = act_values["t_shift"]; + if (act_values.contains("tau_1")) act_params.tau_1 = act_values["tau_1"]; + if (act_values.contains("tau_2")) act_params.tau_2 = act_values["tau_2"]; + if (act_values.contains("m1")) act_params.m1 = act_values["m1"]; + if (act_values.contains("m2")) act_params.m2 = act_values["m2"]; } - // Remove the nested dict + // Remove the nested dict (don't flatten to avoid parameter conflicts) chamber_values.erase("activation_function_values"); + } else { + // Handle flat parameter structure for backward compatibility + if (act_params.type == ActivationType::HALF_COSINE) { + if (chamber_values.contains("t_active")) act_params.t_active = chamber_values["t_active"]; + if (chamber_values.contains("t_twitch")) act_params.t_twitch = chamber_values["t_twitch"]; + } else if (act_params.type == ActivationType::PIECEWISE_COSINE) { + if (chamber_values.contains("contract_start")) act_params.contract_start = chamber_values["contract_start"]; + if (chamber_values.contains("relax_start")) act_params.relax_start = chamber_values["relax_start"]; + if (chamber_values.contains("contract_duration")) act_params.contract_duration = chamber_values["contract_duration"]; + if (chamber_values.contains("relax_duration")) act_params.relax_duration = chamber_values["relax_duration"]; + } else if (act_params.type == ActivationType::TWO_HILL) { + if (chamber_values.contains("t_shift")) act_params.t_shift = chamber_values["t_shift"]; + if (chamber_values.contains("tau_1")) act_params.tau_1 = chamber_values["tau_1"]; + if (chamber_values.contains("tau_2")) act_params.tau_2 = chamber_values["tau_2"]; + if (chamber_values.contains("m1")) act_params.m1 = chamber_values["m1"]; + if (chamber_values.contains("m2")) act_params.m2 = chamber_values["m2"]; + } + } + + // Create the chamber block + int block_id = generate_block(model, chamber_values, chamber_type, chamber_name); + + // Set activation parameters on the chamber + auto chamber_block = model.get_block(chamber_name); + if (chamber_type == "ChamberElastanceInductor") { + auto chamber = dynamic_cast(chamber_block); + if (chamber) { + chamber->set_activation_params(act_params); + } + } else if (chamber_type == "LinearElastanceChamber" || chamber_type == "PiecewiseCosineChamber") { + auto chamber = dynamic_cast(chamber_block); + if (chamber) { + chamber->set_activation_params(act_params); + } } - generate_block(model, chamber_values, chamber_type, chamber_name); DEBUG_MSG("Created chamber " << chamber_name); } } diff --git a/src/solve/SimulationParameters.h b/src/solve/SimulationParameters.h index 512a5aaac..34e09d1d1 100644 --- a/src/solve/SimulationParameters.h +++ b/src/solve/SimulationParameters.h @@ -12,6 +12,9 @@ #include #include +#include "ActivationFunction.h" +#include "ChamberElastanceInductor.h" +#include "LinearElastanceChamber.h" #include "Model.h" #include "State.h" #include "debug.h" From fbc1dc88204f2ee3fdf8dc6d7299c9c372e1e6c6 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 10 Feb 2026 21:12:45 +0000 Subject: [PATCH 10/23] Store cardiac_period as member of ActivationFunction base class Co-authored-by: aabrown100-git <71424733+aabrown100-git@users.noreply.github.com> --- src/model/ActivationFunction.h | 47 +++++++++++++++----------- src/model/ChamberElastanceInductor.cpp | 4 +-- src/model/LinearElastanceChamber.cpp | 4 +-- 3 files changed, 29 insertions(+), 26 deletions(-) diff --git a/src/model/ActivationFunction.h b/src/model/ActivationFunction.h index 9d435ccdf..524538bb6 100644 --- a/src/model/ActivationFunction.h +++ b/src/model/ActivationFunction.h @@ -67,10 +67,9 @@ class ActivationFunction { * @brief Compute activation value at given time * * @param time Current time - * @param cardiac_period Period of cardiac cycle * @return Activation value between 0 and 1 */ - virtual double compute(double time, double cardiac_period) = 0; + virtual double compute(double time) = 0; /** * @brief Get the type of activation function @@ -87,6 +86,9 @@ class ActivationFunction { */ static std::unique_ptr create( const ActivationFunctionParams& params); + + protected: + double cardiac_period_; ///< Cardiac cycle period }; /** @@ -112,19 +114,21 @@ class HalfCosineActivation : public ActivationFunction { * * @param t_active Time when activation begins within cardiac cycle * @param t_twitch Duration of the contraction twitch + * @param cardiac_period Cardiac cycle period */ - HalfCosineActivation(double t_active, double t_twitch) - : t_active_(t_active), t_twitch_(t_twitch) {} + HalfCosineActivation(double t_active, double t_twitch, double cardiac_period) + : t_active_(t_active), t_twitch_(t_twitch) { + cardiac_period_ = cardiac_period; + } /** * @brief Compute activation value * * @param time Current time - * @param cardiac_period Period of cardiac cycle * @return Activation value between 0 and 1 */ - double compute(double time, double cardiac_period) override { - double t_in_cycle = std::fmod(time, cardiac_period); + double compute(double time) override { + double t_in_cycle = std::fmod(time, cardiac_period_); double t_contract = 0.0; if (t_in_cycle >= t_active_) { @@ -179,30 +183,33 @@ class PiecewiseCosineActivation : public ActivationFunction { * @param relax_start Time when relaxation starts * @param contract_duration Duration of contraction phase * @param relax_duration Duration of relaxation phase + * @param cardiac_period Cardiac cycle period */ PiecewiseCosineActivation(double contract_start, double relax_start, - double contract_duration, double relax_duration) + double contract_duration, double relax_duration, + double cardiac_period) : contract_start_(contract_start), relax_start_(relax_start), contract_duration_(contract_duration), - relax_duration_(relax_duration) {} + relax_duration_(relax_duration) { + cardiac_period_ = cardiac_period; + } /** * @brief Compute activation value * * @param time Current time - * @param cardiac_period Period of cardiac cycle * @return Activation value between 0 and 1 */ - double compute(double time, double cardiac_period) override { + double compute(double time) override { double phi = 0.0; - double piecewise_condition = std::fmod(time - contract_start_, cardiac_period); + double piecewise_condition = std::fmod(time - contract_start_, cardiac_period_); if (0.0 <= piecewise_condition && piecewise_condition < contract_duration_) { phi = 0.5 * (1.0 - std::cos((M_PI * piecewise_condition) / contract_duration_)); } else { - piecewise_condition = std::fmod(time - relax_start_, cardiac_period); + piecewise_condition = std::fmod(time - relax_start_, cardiac_period_); if (0.0 <= piecewise_condition && piecewise_condition < relax_duration_) { phi = 0.5 * (1.0 + std::cos((M_PI * piecewise_condition) / relax_duration_)); } @@ -269,6 +276,7 @@ class TwoHillActivation : public ActivationFunction { m2_(m2), normalization_factor_(1.0), normalization_initialized_(false) { + cardiac_period_ = cardiac_period; initialize_normalization(cardiac_period); } @@ -276,19 +284,18 @@ class TwoHillActivation : public ActivationFunction { * @brief Compute activation value * * @param time Current time - * @param cardiac_period Period of cardiac cycle * @return Activation value between 0 and 1 */ - double compute(double time, double cardiac_period) override { + double compute(double time) override { if (!normalization_initialized_) { throw std::runtime_error("TwoHillActivation: Normalization not initialized"); } - double t_in_cycle = std::fmod(time, cardiac_period); + double t_in_cycle = std::fmod(time, cardiac_period_); // Compute shifted time (handle negative modulo) - double t_shifted = std::fmod(t_in_cycle - t_shift_, cardiac_period); - t_shifted = (t_shifted >= 0.0) ? t_shifted : t_shifted + cardiac_period; + double t_shifted = std::fmod(t_in_cycle - t_shift_, cardiac_period_); + t_shifted = (t_shifted >= 0.0) ? t_shifted : t_shifted + cardiac_period_; // Compute hill functions double g1 = std::pow(t_shifted / tau_1_, m1_); @@ -354,12 +361,12 @@ inline std::unique_ptr ActivationFunction::create( switch (params.type) { case ActivationType::HALF_COSINE: return std::make_unique( - params.t_active, params.t_twitch); + params.t_active, params.t_twitch, params.cardiac_period); case ActivationType::PIECEWISE_COSINE: return std::make_unique( params.contract_start, params.relax_start, - params.contract_duration, params.relax_duration); + params.contract_duration, params.relax_duration, params.cardiac_period); case ActivationType::TWO_HILL: return std::make_unique( diff --git a/src/model/ChamberElastanceInductor.cpp b/src/model/ChamberElastanceInductor.cpp index 5888caf4b..9ec9d1b74 100644 --- a/src/model/ChamberElastanceInductor.cpp +++ b/src/model/ChamberElastanceInductor.cpp @@ -45,11 +45,9 @@ void ChamberElastanceInductor::get_elastance_values( double Emin = parameters[global_param_ids[ParamId::EMIN]]; double Vrd = parameters[global_param_ids[ParamId::VRD]]; double Vrs = parameters[global_param_ids[ParamId::VRS]]; - - auto T_cardiac = model->cardiac_cycle_period; // Compute activation using the activation function - double act = activation_func_->compute(model->time, T_cardiac); + double act = activation_func_->compute(model->time); Vrest = (1.0 - act) * (Vrd - Vrs) + Vrs; Elas = (Emax - Emin) * act + Emin; diff --git a/src/model/LinearElastanceChamber.cpp b/src/model/LinearElastanceChamber.cpp index 60f5ca065..1d3bd33c4 100644 --- a/src/model/LinearElastanceChamber.cpp +++ b/src/model/LinearElastanceChamber.cpp @@ -42,11 +42,9 @@ void LinearElastanceChamber::get_elastance_values( std::vector& parameters) { double Emax = parameters[global_param_ids[ParamId::EMAX]]; double Epass = parameters[global_param_ids[ParamId::EPASS]]; - - auto T_cardiac = model->cardiac_cycle_period; // Compute activation using the activation function - double phi = activation_func_->compute(model->time, T_cardiac); + double phi = activation_func_->compute(model->time); Elas = Epass + Emax * phi; } From 0f53d8c22f511215b6bba9da1feb3c384dac1f30 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 10 Feb 2026 21:16:40 +0000 Subject: [PATCH 11/23] Make PiecewiseCosineChamber an alias to LinearElastanceChamber for backward compatibility Co-authored-by: aabrown100-git <71424733+aabrown100-git@users.noreply.github.com> --- src/model/BlockType.h | 2 +- src/model/PiecewiseCosineChamber.cpp | 62 +-------- src/model/PiecewiseCosineChamber.h | 192 +-------------------------- 3 files changed, 11 insertions(+), 245 deletions(-) diff --git a/src/model/BlockType.h b/src/model/BlockType.h index 76d862c96..00af9d041 100644 --- a/src/model/BlockType.h +++ b/src/model/BlockType.h @@ -30,7 +30,7 @@ enum class BlockType { chamber_elastance_inductor = 14, chamber_sphere = 15, blood_vessel_CRL = 16, - piecewise_cosine_chamber = 17, + piecewise_cosine_chamber = 17, ///< Deprecated: Use linear_elastance_chamber piecewise_valve = 18, linear_elastance_chamber = 19 }; diff --git a/src/model/PiecewiseCosineChamber.cpp b/src/model/PiecewiseCosineChamber.cpp index 010016576..bb540c765 100644 --- a/src/model/PiecewiseCosineChamber.cpp +++ b/src/model/PiecewiseCosineChamber.cpp @@ -1,62 +1,6 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "PiecewiseCosineChamber.h" - -void PiecewiseCosineChamber::setup_dofs(DOFHandler& dofhandler) { - // Internal variable is chamber volume - Block::setup_dofs_(dofhandler, 3, {"Vc"}); -} - -void PiecewiseCosineChamber::update_constant(SparseSystem& system, - std::vector& parameters) { - // Eq 0: P_in - E(t)(Vc - Vrest) = 0 - system.F.coeffRef(global_eqn_ids[0], global_var_ids[0]) = 1.0; - - // Eq 1: P_in - P_out = 0 - system.F.coeffRef(global_eqn_ids[1], global_var_ids[0]) = 1.0; - system.F.coeffRef(global_eqn_ids[1], global_var_ids[2]) = -1.0; - - // Eq 2: Q_in - Q_out - dVc = 0 - system.F.coeffRef(global_eqn_ids[2], global_var_ids[1]) = 1.0; - system.F.coeffRef(global_eqn_ids[2], global_var_ids[3]) = -1.0; - system.E.coeffRef(global_eqn_ids[2], global_var_ids[4]) = -1.0; -} - -void PiecewiseCosineChamber::update_time(SparseSystem& system, - std::vector& parameters) { - get_elastance_values(parameters); - - // Eq 0: P_in - E(t)(Vc - Vrest) = P_in - E(t)*Vc + E(t)*Vrest = 0 - system.F.coeffRef(global_eqn_ids[0], global_var_ids[4]) = -1 * Elas; - system.C.coeffRef(global_eqn_ids[0]) = - Elas * parameters[global_param_ids[ParamId::VREST]]; -} - -void PiecewiseCosineChamber::get_elastance_values( - std::vector& parameters) { - double Emax = parameters[global_param_ids[ParamId::EMAX]]; - double Epass = parameters[global_param_ids[ParamId::EPASS]]; - double Vrest = parameters[global_param_ids[ParamId::VREST]]; - double contract_start = parameters[global_param_ids[ParamId::CSTART]]; - double relax_start = parameters[global_param_ids[ParamId::RSTART]]; - double contract_duration = parameters[global_param_ids[ParamId::CDUR]]; - double relax_duration = parameters[global_param_ids[ParamId::RDUR]]; - - auto T_HB = model->cardiac_cycle_period; - - double phi = 0; - - auto piecewise_condition = fmod(model->time - contract_start, T_HB); - - if (0 <= piecewise_condition && piecewise_condition < contract_duration) { - phi = 0.5 * (1 - cos((M_PI * piecewise_condition) / contract_duration)); - } else { - piecewise_condition = fmod(model->time - relax_start, T_HB); - if (0 <= piecewise_condition && piecewise_condition < relax_duration) { - phi = 0.5 * (1 + cos((M_PI * piecewise_condition) / relax_duration)); - } - } - - Elas = Epass + Emax * phi; -} +// This file intentionally left empty. +// PiecewiseCosineChamber is now an alias to LinearElastanceChamber. +// See PiecewiseCosineChamber.h for the alias definition. diff --git a/src/model/PiecewiseCosineChamber.h b/src/model/PiecewiseCosineChamber.h index b9a461503..ba37e11f6 100644 --- a/src/model/PiecewiseCosineChamber.h +++ b/src/model/PiecewiseCosineChamber.h @@ -3,198 +3,20 @@ /** * @file PiecewiseCosineChamber.h - * @brief model::PiecewiseCosineChamber source file + * @brief Backward compatibility alias for LinearElastanceChamber */ #ifndef SVZERODSOLVER_MODEL_PIECEWISECOSINECHAMBER_HPP_ #define SVZERODSOLVER_MODEL_PIECEWISECOSINECHAMBER_HPP_ -#include - -#include "Block.h" -#include "Model.h" -#include "SparseSystem.h" -#include "debug.h" +#include "LinearElastanceChamber.h" /** - * @brief Cardiac chamber with piecewise elastance and inductor. - * - * Models a cardiac chamber as a time-varying capacitor (elastance with - * specified resting volumes) and an inductor. See \cite Regazzoni2022 - * - * This chamber block can be connected to other blocks using junctions. - * - * \f[ - * \begin{circuitikz} - * \draw node[left] {$Q_{in}$} [-latex] (0,0) -- (0.8,0); - * \draw (1,0) node[anchor=south]{$P_{in}$} - * to[short, *-*] (3,0) node[anchor=south]{$P_{out}$} - * (1,0) to [vC, l=$E$, *-] (1,-1.5) - * node[ground]{}; - * \draw (3.2,0) -- (4.0,0) node[right] {$Q_{out}$} [-latex]; - * \end{circuitikz} - * \f] - * - * ### Governing equations - * - * \f[ - * P_{in}-E(t)(V_c-V_{rest})=0 - * \f] - * - * \f[ - * P_{in}-P_{out}=0 - * \f] - * - * \f[ - * Q_{in}-Q_{out}-\dot{V}_c=0 - * \f] - * - * ### Local contributions - * - * \f[ - * \mathbf{y}^{e}=\left[\begin{array}{lllll}P_{in} & Q_{in} & - * P_{out} & Q_{out} & V_c\end{array}\right]^{T} \f] - * - * \f[ - * \mathbf{E}^{e}=\left[\begin{array}{ccccc} - * 0 & 0 & 0 & 0 & 0\\ - * 0 & 0 & 0 & 0 & 0\\ - * 0 & 0 & 0 & 0 & -1 - * \end{array}\right] - * \f] - * - * \f[ - * \mathbf{F}^{e}=\left[\begin{array}{ccccc} - * 1 & 0 & 0 & 0 & E(t) \\ - * 1 & 0 & -1 & 0 & 0 \\ - * 0 & 1 & 0 & -1 & 0 - * \end{array}\right] - * \f] - * - * \f[ - * \mathbf{c}^{e}=\left[\begin{array}{c} - * E(t)V_{rest} \\ - * 0 \\ - * 0 - * \end{array}\right] - * \f] - * - * In the above equations, - * - * \f[ - * E_i(t) = E_i^{\text{pass}} + E_i^{\text{act,max}} \, - * \phi\!\left(t, t_C^i, t_R^i, T_C^i, T_R^i\right), - * \f] - * - * \f[ - * \phi(t, t_C, t_R, T_C, T_R) = - * \begin{cases} - * \frac{1}{2}\left[1 - \cos\!\left(\frac{\pi}{T_C} \operatorname{mod}(t - t_C, - * T_{\mathrm{HB}})\right)\right], & \text{if } 0 \le \operatorname{mod}(t - - * t_C, T_{\mathrm{HB}}) < T_C, \\[1.2em] \frac{1}{2}\left[1 + - * \cos\!\left(\frac{\pi}{T_R} \operatorname{mod}(t - t_R, - * T_{\mathrm{HB}})\right)\right], & \text{if } 0 \le \operatorname{mod}(t - - * t_R, T_{\mathrm{HB}}) < T_R, \\[1.2em] 0, & \text{otherwise.} \end{cases} \f] - * - * ### Parameters - * - * Parameter sequence for constructing this block - * - * * `0` Emax: Maximum elastance - * * `1` Epass: Passive elastance - * * `2` Vrest: Rest diastolic volume - * * `3` contract_start: Contract start time - * * `4` relax_start: Relax start time - * * `5` contract_duration: Contract duration - * * `6` relax_duration: Relaxation duration - * - * ### Internal variables - * - * Names of internal variables in this block's output: - * - * * `Vc`: Chamber volume - * + * @brief Backward compatibility alias for LinearElastanceChamber + * + * PiecewiseCosineChamber has been renamed to LinearElastanceChamber. + * This alias is maintained for backward compatibility with existing configurations. */ -class PiecewiseCosineChamber : public Block { - public: - /** - * @brief Construct a new BloodVessel object - * - * @param id Global ID of the block - * @param model The model to which the block belongs - */ - PiecewiseCosineChamber(int id, Model* model) - : Block(id, model, BlockType::piecewise_cosine_chamber, - BlockClass::chamber, - {{"Emax", InputParameter()}, - {"Epass", InputParameter()}, - {"Vrest", InputParameter()}, - {"contract_start", InputParameter()}, - {"relax_start", InputParameter()}, - {"contract_duration", InputParameter()}, - {"relax_duration", InputParameter()}}) {} - - /** - * @brief Local IDs of the parameters - * - */ - enum ParamId { - EMAX = 0, - EPASS = 1, - VREST = 2, - CSTART = 3, - RSTART = 4, - CDUR = 5, - RDUR = 6 - }; - - /** - * @brief Set up the degrees of freedom (DOF) of the block - * - * Set global_var_ids and global_eqn_ids of the element based on the - * number of equations and the number of internal variables of the - * element. - * - * @param dofhandler Degree-of-freedom handler to register variables and - * equations at - */ - void setup_dofs(DOFHandler& dofhandler); - - /** - * @brief Update the constant contributions of the element in a sparse - system - * - * @param system System to update contributions at - * @param parameters Parameters of the model - */ - void update_constant(SparseSystem& system, std::vector& parameters); - - /** - * @brief Update the time-dependent contributions of the element in a sparse - * system - * - * @param system System to update contributions at - * @param parameters Parameters of the model - */ - void update_time(SparseSystem& system, std::vector& parameters); - - /** - * @brief Number of triplets of element - * - * Number of triplets that the element contributes to the global system - * (relevant for sparse memory reservation) - */ - TripletsContributions num_triplets{6, 2, 0}; - - private: - double Elas; // Chamber Elastance - - /** - * @brief Update the elastance functions which depend on time - * - * @param parameters Parameters of the model - */ - void get_elastance_values(std::vector& parameters); -}; +using PiecewiseCosineChamber = LinearElastanceChamber; #endif // SVZERODSOLVER_MODEL_PIECEWISECOSINECHAMBER_HPP_ From 7038e2d50a1e7246d77a5a7718c5454f869ec0e1 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 10 Feb 2026 21:19:23 +0000 Subject: [PATCH 12/23] Remove backward compatibility code and require activation_type parameter Co-authored-by: aabrown100-git <71424733+aabrown100-git@users.noreply.github.com> --- src/model/ChamberElastanceInductor.cpp | 11 --------- src/model/ChamberElastanceInductor.h | 8 ++----- src/model/LinearElastanceChamber.cpp | 32 +++++++------------------- 3 files changed, 10 insertions(+), 41 deletions(-) diff --git a/src/model/ChamberElastanceInductor.cpp b/src/model/ChamberElastanceInductor.cpp index 9ec9d1b74..d710dbdef 100644 --- a/src/model/ChamberElastanceInductor.cpp +++ b/src/model/ChamberElastanceInductor.cpp @@ -66,17 +66,6 @@ void ChamberElastanceInductor::initialize_activation_function( // Set cardiac period activation_params_.cardiac_period = model->cardiac_cycle_period; - // For backward compatibility, if activation parameters weren't set from JSON, - // try to get them from the parameter vector (for old flat format) - if (activation_params_.type == ActivationType::HALF_COSINE) { - if (global_param_ids.count(ParamId::TACTIVE) > 0) { - activation_params_.t_active = parameters[global_param_ids[ParamId::TACTIVE]]; - } - if (global_param_ids.count(ParamId::TTWITCH) > 0) { - activation_params_.t_twitch = parameters[global_param_ids[ParamId::TTWITCH]]; - } - } - // Use factory method to create activation function activation_func_ = ActivationFunction::create(activation_params_); } diff --git a/src/model/ChamberElastanceInductor.h b/src/model/ChamberElastanceInductor.h index dd8077264..2357a9b88 100644 --- a/src/model/ChamberElastanceInductor.h +++ b/src/model/ChamberElastanceInductor.h @@ -154,8 +154,6 @@ class ChamberElastanceInductor : public Block { {"Emin", InputParameter()}, {"Vrd", InputParameter()}, {"Vrs", InputParameter()}, - {"t_active", InputParameter(false)}, - {"t_twitch", InputParameter(false)}, {"Impedance", InputParameter()}, {"activation_type", InputParameter(false)}}) {} @@ -168,10 +166,8 @@ class ChamberElastanceInductor : public Block { EMIN = 1, VRD = 2, VRS = 3, - TACTIVE = 4, - TTWITCH = 5, - IMPEDANCE = 6, - ACTIVATION_TYPE = 7 + IMPEDANCE = 4, + ACTIVATION_TYPE = 5 }; /** diff --git a/src/model/LinearElastanceChamber.cpp b/src/model/LinearElastanceChamber.cpp index 1d3bd33c4..42a4d5128 100644 --- a/src/model/LinearElastanceChamber.cpp +++ b/src/model/LinearElastanceChamber.cpp @@ -51,35 +51,19 @@ void LinearElastanceChamber::get_elastance_values( void LinearElastanceChamber::initialize_activation_function( std::vector& parameters) { - // Check if activation_type parameter is provided (optional parameter) - // Default to PIECEWISE_COSINE (1) for backward compatibility - int activation_type_int = 1; - if (global_param_ids.count(ParamId::ACTIVATION_TYPE) > 0) { - activation_type_int = static_cast( - parameters[global_param_ids[ParamId::ACTIVATION_TYPE]]); - activation_params_.type = static_cast(activation_type_int); + // Check if activation_type parameter is provided (required parameter) + if (global_param_ids.count(ParamId::ACTIVATION_TYPE) == 0) { + throw std::runtime_error( + "LinearElastanceChamber: activation_type parameter is required"); } + + int activation_type_int = static_cast( + parameters[global_param_ids[ParamId::ACTIVATION_TYPE]]); + activation_params_.type = static_cast(activation_type_int); // Set cardiac period activation_params_.cardiac_period = model->cardiac_cycle_period; - // For backward compatibility, if activation parameters weren't set from JSON, - // try to get them from the parameter vector (for old flat format) - if (activation_params_.type == ActivationType::PIECEWISE_COSINE) { - if (global_param_ids.count(ParamId::CSTART) > 0) { - activation_params_.contract_start = parameters[global_param_ids[ParamId::CSTART]]; - } - if (global_param_ids.count(ParamId::RSTART) > 0) { - activation_params_.relax_start = parameters[global_param_ids[ParamId::RSTART]]; - } - if (global_param_ids.count(ParamId::CDUR) > 0) { - activation_params_.contract_duration = parameters[global_param_ids[ParamId::CDUR]]; - } - if (global_param_ids.count(ParamId::RDUR) > 0) { - activation_params_.relax_duration = parameters[global_param_ids[ParamId::RDUR]]; - } - } - // Use factory method to create activation function activation_func_ = ActivationFunction::create(activation_params_); } From aa4787df98e27a44833b08e2a30d6cc10faa3c9a Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 10 Feb 2026 21:22:41 +0000 Subject: [PATCH 13/23] Remove activation function parameters from LinearElastanceChamber and flat parameter backward compatibility Co-authored-by: aabrown100-git <71424733+aabrown100-git@users.noreply.github.com> --- src/model/LinearElastanceChamber.h | 10 +--------- src/solve/SimulationParameters.cpp | 17 ----------------- 2 files changed, 1 insertion(+), 26 deletions(-) diff --git a/src/model/LinearElastanceChamber.h b/src/model/LinearElastanceChamber.h index 8b744abcb..0c1000546 100644 --- a/src/model/LinearElastanceChamber.h +++ b/src/model/LinearElastanceChamber.h @@ -131,10 +131,6 @@ class LinearElastanceChamber : public Block { {{"Emax", InputParameter()}, {"Epass", InputParameter()}, {"Vrest", InputParameter()}, - {"contract_start", InputParameter(false)}, - {"relax_start", InputParameter(false)}, - {"contract_duration", InputParameter(false)}, - {"relax_duration", InputParameter(false)}, {"activation_type", InputParameter(false)}}) {} /** @@ -145,11 +141,7 @@ class LinearElastanceChamber : public Block { EMAX = 0, EPASS = 1, VREST = 2, - CSTART = 3, - RSTART = 4, - CDUR = 5, - RDUR = 6, - ACTIVATION_TYPE = 7 + ACTIVATION_TYPE = 3 }; /** diff --git a/src/solve/SimulationParameters.cpp b/src/solve/SimulationParameters.cpp index a1045f604..bb8c7312d 100644 --- a/src/solve/SimulationParameters.cpp +++ b/src/solve/SimulationParameters.cpp @@ -603,23 +603,6 @@ void create_chambers( // Remove the nested dict (don't flatten to avoid parameter conflicts) chamber_values.erase("activation_function_values"); - } else { - // Handle flat parameter structure for backward compatibility - if (act_params.type == ActivationType::HALF_COSINE) { - if (chamber_values.contains("t_active")) act_params.t_active = chamber_values["t_active"]; - if (chamber_values.contains("t_twitch")) act_params.t_twitch = chamber_values["t_twitch"]; - } else if (act_params.type == ActivationType::PIECEWISE_COSINE) { - if (chamber_values.contains("contract_start")) act_params.contract_start = chamber_values["contract_start"]; - if (chamber_values.contains("relax_start")) act_params.relax_start = chamber_values["relax_start"]; - if (chamber_values.contains("contract_duration")) act_params.contract_duration = chamber_values["contract_duration"]; - if (chamber_values.contains("relax_duration")) act_params.relax_duration = chamber_values["relax_duration"]; - } else if (act_params.type == ActivationType::TWO_HILL) { - if (chamber_values.contains("t_shift")) act_params.t_shift = chamber_values["t_shift"]; - if (chamber_values.contains("tau_1")) act_params.tau_1 = chamber_values["tau_1"]; - if (chamber_values.contains("tau_2")) act_params.tau_2 = chamber_values["tau_2"]; - if (chamber_values.contains("m1")) act_params.m1 = chamber_values["m1"]; - if (chamber_values.contains("m2")) act_params.m2 = chamber_values["m2"]; - } } // Create the chamber block From 5d80d12da4e51ea24ce07864fc35138c946572ba Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 10 Feb 2026 21:25:25 +0000 Subject: [PATCH 14/23] Remove PiecewiseCosineChamber from factory registration and type checks Co-authored-by: aabrown100-git <71424733+aabrown100-git@users.noreply.github.com> --- src/model/Model.cpp | 1 - src/solve/SimulationParameters.cpp | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/model/Model.cpp b/src/model/Model.cpp index f4d2143eb..2344a6b60 100644 --- a/src/model/Model.cpp +++ b/src/model/Model.cpp @@ -31,7 +31,6 @@ Model::Model() { {"ChamberElastanceInductor", block_factory()}, {"BloodVesselCRL", block_factory()}, {"PiecewiseValve", block_factory()}, - {"PiecewiseCosineChamber", block_factory()}, {"LinearElastanceChamber", block_factory()}}; } diff --git a/src/solve/SimulationParameters.cpp b/src/solve/SimulationParameters.cpp index bb8c7312d..ac0787b60 100644 --- a/src/solve/SimulationParameters.cpp +++ b/src/solve/SimulationParameters.cpp @@ -615,7 +615,7 @@ void create_chambers( if (chamber) { chamber->set_activation_params(act_params); } - } else if (chamber_type == "LinearElastanceChamber" || chamber_type == "PiecewiseCosineChamber") { + } else if (chamber_type == "LinearElastanceChamber") { auto chamber = dynamic_cast(chamber_block); if (chamber) { chamber->set_activation_params(act_params); From 8031791122b76a8f5b2b06dfdea19feb4573124d Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 10 Feb 2026 21:34:58 +0000 Subject: [PATCH 15/23] Remove PiecewiseCosineChamber completely - no longer needed with activation function separation Co-authored-by: aabrown100-git <71424733+aabrown100-git@users.noreply.github.com> --- docs/pages/activation_functions.md | 2 +- src/model/ActivationFunction.h | 2 +- src/model/BlockType.h | 1 - src/model/CMakeLists.txt | 2 - src/model/Model.h | 1 - src/model/PiecewiseCosineChamber.cpp | 6 --- src/model/PiecewiseCosineChamber.h | 22 -------- src/solve/csv_writer.cpp | 2 +- .../linear_chamber_two_hill_activation.json | 2 +- tests/cases/piecewise_Chamber_and_Valve.json | 52 ++++++++++++------- 10 files changed, 36 insertions(+), 56 deletions(-) delete mode 100644 src/model/PiecewiseCosineChamber.cpp delete mode 100644 src/model/PiecewiseCosineChamber.h diff --git a/docs/pages/activation_functions.md b/docs/pages/activation_functions.md index e9b18208d..b83c83170 100644 --- a/docs/pages/activation_functions.md +++ b/docs/pages/activation_functions.md @@ -60,7 +60,7 @@ Or with explicit activation function specification: ### 2. Piecewise Cosine Activation (Type 1) -The default activation function for `LinearElastanceChamber` (formerly `PiecewiseCosineChamber`). It models separate contraction and relaxation phases: +The default activation function for `LinearElastanceChamber`. It models separate contraction and relaxation phases: ``` φ(t) = 0.5 * [1 - cos(π * (t - t_C) / T_C)] during contraction diff --git a/src/model/ActivationFunction.h b/src/model/ActivationFunction.h index 524538bb6..ff1a3950b 100644 --- a/src/model/ActivationFunction.h +++ b/src/model/ActivationFunction.h @@ -160,7 +160,7 @@ class HalfCosineActivation : public ActivationFunction { /** * @brief Piecewise cosine activation function * - * This implements the activation function from the PiecewiseCosineChamber + * This implements the activation function from the LinearElastanceChamber * (Regazzoni chamber model). The activation consists of separate contraction * and relaxation phases, each following a cosine curve. * diff --git a/src/model/BlockType.h b/src/model/BlockType.h index 00af9d041..aff0b7ce6 100644 --- a/src/model/BlockType.h +++ b/src/model/BlockType.h @@ -30,7 +30,6 @@ enum class BlockType { chamber_elastance_inductor = 14, chamber_sphere = 15, blood_vessel_CRL = 16, - piecewise_cosine_chamber = 17, ///< Deprecated: Use linear_elastance_chamber piecewise_valve = 18, linear_elastance_chamber = 19 }; diff --git a/src/model/CMakeLists.txt b/src/model/CMakeLists.txt index 62e1fe414..c6128e356 100644 --- a/src/model/CMakeLists.txt +++ b/src/model/CMakeLists.txt @@ -29,7 +29,6 @@ set(CXXSRCS ResistiveJunction.cpp ValveTanh.cpp WindkesselBC.cpp - PiecewiseCosineChamber.cpp LinearElastanceChamber.cpp PiecewiseValve.cpp ) @@ -60,7 +59,6 @@ set(HDRS ResistiveJunction.h ValveTanh.h WindkesselBC.h - PiecewiseCosineChamber.h LinearElastanceChamber.h PiecewiseValve.h ) diff --git a/src/model/Model.h b/src/model/Model.h index 23b673e58..c837af5e6 100644 --- a/src/model/Model.h +++ b/src/model/Model.h @@ -33,7 +33,6 @@ #include "Node.h" #include "OpenLoopCoronaryBC.h" #include "Parameter.h" -#include "PiecewiseCosineChamber.h" #include "PiecewiseValve.h" #include "PressureReferenceBC.h" #include "ResistanceBC.h" diff --git a/src/model/PiecewiseCosineChamber.cpp b/src/model/PiecewiseCosineChamber.cpp deleted file mode 100644 index bb540c765..000000000 --- a/src/model/PiecewiseCosineChamber.cpp +++ /dev/null @@ -1,6 +0,0 @@ -// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the -// University of California, and others. SPDX-License-Identifier: BSD-3-Clause - -// This file intentionally left empty. -// PiecewiseCosineChamber is now an alias to LinearElastanceChamber. -// See PiecewiseCosineChamber.h for the alias definition. diff --git a/src/model/PiecewiseCosineChamber.h b/src/model/PiecewiseCosineChamber.h deleted file mode 100644 index ba37e11f6..000000000 --- a/src/model/PiecewiseCosineChamber.h +++ /dev/null @@ -1,22 +0,0 @@ -// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the -// University of California, and others. SPDX-License-Identifier: BSD-3-Clause - -/** - * @file PiecewiseCosineChamber.h - * @brief Backward compatibility alias for LinearElastanceChamber - */ - -#ifndef SVZERODSOLVER_MODEL_PIECEWISECOSINECHAMBER_HPP_ -#define SVZERODSOLVER_MODEL_PIECEWISECOSINECHAMBER_HPP_ - -#include "LinearElastanceChamber.h" - -/** - * @brief Backward compatibility alias for LinearElastanceChamber - * - * PiecewiseCosineChamber has been renamed to LinearElastanceChamber. - * This alias is maintained for backward compatibility with existing configurations. - */ -using PiecewiseCosineChamber = LinearElastanceChamber; - -#endif // SVZERODSOLVER_MODEL_PIECEWISECOSINECHAMBER_HPP_ diff --git a/src/solve/csv_writer.cpp b/src/solve/csv_writer.cpp index 4e521597a..8dbdda49b 100644 --- a/src/solve/csv_writer.cpp +++ b/src/solve/csv_writer.cpp @@ -45,7 +45,7 @@ std::string to_vessel_csv(const std::vector& times, if (dynamic_cast(block) == nullptr && dynamic_cast(block) == nullptr && - dynamic_cast(block) == nullptr && + dynamic_cast(block) == nullptr && dynamic_cast(block) == nullptr && dynamic_cast(block) == nullptr) { continue; diff --git a/tests/cases/linear_chamber_two_hill_activation.json b/tests/cases/linear_chamber_two_hill_activation.json index f2a03469e..5b0f6bf8f 100644 --- a/tests/cases/linear_chamber_two_hill_activation.json +++ b/tests/cases/linear_chamber_two_hill_activation.json @@ -2,7 +2,7 @@ "description": { "description of test case": "LinearElastanceChamber with two hill activation function", "analytical results": [ - "Testing LinearElastanceChamber (formerly PiecewiseCosineChamber) with two hill activation" + "Testing LinearElastanceChamber with two hill activation" ] }, "simulation_parameters": { diff --git a/tests/cases/piecewise_Chamber_and_Valve.json b/tests/cases/piecewise_Chamber_and_Valve.json index e4ec38b93..1efe3d2a0 100644 --- a/tests/cases/piecewise_Chamber_and_Valve.json +++ b/tests/cases/piecewise_Chamber_and_Valve.json @@ -126,55 +126,67 @@ "boundary_conditions": [], "chambers": [ { - "type": "PiecewiseCosineChamber", + "type": "LinearElastanceChamber", "name": "right_atrium", "values": { "Emax": 199.95, "Epass": 89.80375933024839, "Vrest": 41.680015938842274, - "contract_start": 0.025, - "relax_start": 0.08625, - "contract_duration": 0.06125, - "relax_duration": 0.18375 + "activation_function_type": "piecewise_cosine", + "activation_function_values": { + "contract_start": 0.025, + "relax_start": 0.08625, + "contract_duration": 0.06125, + "relax_duration": 0.18375 + } } }, { - "type": "PiecewiseCosineChamber", + "type": "LinearElastanceChamber", "name": "right_ventricle", "values": { "Emax": 1662.0158240056528, "Epass": 40.85565535747109, "Vrest": 72.05452710344869, - "contract_start": 0.207, - "relax_start": 0.29625, - "contract_duration": 0.08925, - "relax_duration": 0.26975 + "activation_function_type": "piecewise_cosine", + "activation_function_values": { + "contract_start": 0.207, + "relax_start": 0.29625, + "contract_duration": 0.08925, + "relax_duration": 0.26975 + } } }, { - "type": "PiecewiseCosineChamber", + "type": "LinearElastanceChamber", "name": "left_atrium", "values": { "Emax": 199.95, "Epass": 260.59064494602634, "Vrest": 26.23534455334443, - "contract_start": 0.025, - "relax_start": 0.08625, - "contract_duration": 0.06125, - "relax_duration": 0.18375 + "activation_function_type": "piecewise_cosine", + "activation_function_values": { + "contract_start": 0.025, + "relax_start": 0.08625, + "contract_duration": 0.06125, + "relax_duration": 0.18375 + } } }, { - "type": "PiecewiseCosineChamber", + "type": "LinearElastanceChamber", "name": "left_ventricle", "values": { "Emax": 17837.499965252825, "Epass": 122.82901158252674, "Vrest": 32.41857776349184, - "contract_start": 0.207, - "relax_start": 0.29625, - "contract_duration": 0.08925, - "relax_duration": 0.26975 + "activation_function_type": "piecewise_cosine", + "activation_function_values": { + "contract_start": 0.207, + "relax_start": 0.29625, + "contract_duration": 0.08925, + "relax_duration": 0.26975 + } } } ], From 0ce3c3247e9fac50344639d757a5867329e5b674 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Wed, 11 Feb 2026 01:52:33 +0000 Subject: [PATCH 16/23] Improve error reporting to list all missing parameters at once Co-authored-by: ncdorn <126125269+ncdorn@users.noreply.github.com> --- src/solve/SimulationParameters.cpp | 78 ++++++++++++++++++++++++------ 1 file changed, 62 insertions(+), 16 deletions(-) diff --git a/src/solve/SimulationParameters.cpp b/src/solve/SimulationParameters.cpp index ac0787b60..a91a611c3 100644 --- a/src/solve/SimulationParameters.cpp +++ b/src/solve/SimulationParameters.cpp @@ -79,6 +79,9 @@ int generate_block(Model& model, const nlohmann::json& block_params_json, } } } else { + // First pass: collect all missing parameters + std::vector missing_params; + for (const auto& block_param : block->input_params) { // Time parameter is read at the same time as time-dependent value if (block_param.first.compare("t") == 0) { @@ -90,42 +93,85 @@ int generate_block(Model& model, const nlohmann::json& block_params_json, continue; } - // Get vector parameter + // Check vector parameter if (block_param.second.is_array) { - // Get parameter vector + // Check parameter vector std::vector val; err = get_param_vector(block_params_json, block_param.first, block_param.second, val); if (err) { - throw std::runtime_error("Array parameter " + block_param.first + - " is mandatory in " + block_type + - " block " + static_cast(name)); + missing_params.push_back(block_param.first + " (array)"); } - // Get time vector + // Check time vector InputParameter t_param{false, true}; std::vector time; err = get_param_vector(block_params_json, "t", t_param, time); if (err) { - throw std::runtime_error("Array parameter t is mandatory in " + - block_type + " block " + - static_cast(name)); + // Only add "t" if not already in the list + if (std::find(missing_params.begin(), missing_params.end(), "t (array)") == missing_params.end()) { + missing_params.push_back("t (array)"); + } } - - // Add parameters to model - new_id = model.add_parameter(time, val, periodic); } - // Get scalar parameter + // Check scalar parameter else { double val; err = get_param_scalar(block_params_json, block_param.first, block_param.second, val); if (err) { - throw std::runtime_error("Scalar parameter " + block_param.first + - " is mandatory in " + block_type + - " block " + static_cast(name)); + missing_params.push_back(block_param.first); + } + } + } + + // If any parameters are missing, throw error with complete list + if (!missing_params.empty()) { + std::string error_msg = "Missing required parameter(s) in " + block_type + + " block " + static_cast(name) + ": "; + for (size_t i = 0; i < missing_params.size(); ++i) { + error_msg += missing_params[i]; + if (i < missing_params.size() - 1) { + error_msg += ", "; } + } + throw std::runtime_error(error_msg); + } + + // Second pass: actually read and store the parameters + for (const auto& block_param : block->input_params) { + // Time parameter is read at the same time as time-dependent value + if (block_param.first.compare("t") == 0) { + continue; + } + + // Skip reading parameters that are not a number + if (!block_param.second.is_number) { + continue; + } + + // Get vector parameter + if (block_param.second.is_array) { + // Get parameter vector + std::vector val; + get_param_vector(block_params_json, block_param.first, + block_param.second, val); + + // Get time vector + InputParameter t_param{false, true}; + std::vector time; + get_param_vector(block_params_json, "t", t_param, time); + + // Add parameters to model + new_id = model.add_parameter(time, val, periodic); + } + + // Get scalar parameter + else { + double val; + get_param_scalar(block_params_json, block_param.first, + block_param.second, val); // Add parameter to model new_id = model.add_parameter(val); From 3c417709fa6b0ff540acf7c50466d4867df33d47 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Wed, 11 Feb 2026 01:59:59 +0000 Subject: [PATCH 17/23] Fix BlockType enum values to be sequential (piecewise_valve = 17, linear_elastance_chamber = 18) Co-authored-by: aabrown100-git <71424733+aabrown100-git@users.noreply.github.com> --- src/model/BlockType.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/model/BlockType.h b/src/model/BlockType.h index aff0b7ce6..2284b33b5 100644 --- a/src/model/BlockType.h +++ b/src/model/BlockType.h @@ -30,8 +30,8 @@ enum class BlockType { chamber_elastance_inductor = 14, chamber_sphere = 15, blood_vessel_CRL = 16, - piecewise_valve = 18, - linear_elastance_chamber = 19 + piecewise_valve = 17, + linear_elastance_chamber = 18 }; /** From 6d09f571cc6b81a0186f3b86bc40a32b5cdf8536 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 12 Feb 2026 01:01:06 +0000 Subject: [PATCH 18/23] Fix compilation error: use vector size check instead of non-existent count() method Co-authored-by: ncdorn <126125269+ncdorn@users.noreply.github.com> --- src/model/ChamberElastanceInductor.cpp | 6 +++--- src/model/ChamberElastanceInductor.h | 2 +- src/model/LinearElastanceChamber.cpp | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/model/ChamberElastanceInductor.cpp b/src/model/ChamberElastanceInductor.cpp index d710dbdef..ca85147a0 100644 --- a/src/model/ChamberElastanceInductor.cpp +++ b/src/model/ChamberElastanceInductor.cpp @@ -56,9 +56,9 @@ void ChamberElastanceInductor::get_elastance_values( void ChamberElastanceInductor::initialize_activation_function( std::vector& parameters) { // Check if activation_type parameter is provided (optional parameter) - int activation_type_int = 0; // Default to HALF_COSINE - if (global_param_ids.count(ParamId::ACTIVATION_TYPE) > 0) { - activation_type_int = static_cast( + // Defaults to HALF_COSINE if not provided + if (global_param_ids.size() > ParamId::ACTIVATION_TYPE) { + int activation_type_int = static_cast( parameters[global_param_ids[ParamId::ACTIVATION_TYPE]]); activation_params_.type = static_cast(activation_type_int); } diff --git a/src/model/ChamberElastanceInductor.h b/src/model/ChamberElastanceInductor.h index 2357a9b88..53b56799c 100644 --- a/src/model/ChamberElastanceInductor.h +++ b/src/model/ChamberElastanceInductor.h @@ -155,7 +155,7 @@ class ChamberElastanceInductor : public Block { {"Vrd", InputParameter()}, {"Vrs", InputParameter()}, {"Impedance", InputParameter()}, - {"activation_type", InputParameter(false)}}) {} + {"activation_type", InputParameter(true, false, true, 0.0)}}) {} /** * @brief Local IDs of the parameters diff --git a/src/model/LinearElastanceChamber.cpp b/src/model/LinearElastanceChamber.cpp index 42a4d5128..b1f73cf46 100644 --- a/src/model/LinearElastanceChamber.cpp +++ b/src/model/LinearElastanceChamber.cpp @@ -52,7 +52,7 @@ void LinearElastanceChamber::get_elastance_values( void LinearElastanceChamber::initialize_activation_function( std::vector& parameters) { // Check if activation_type parameter is provided (required parameter) - if (global_param_ids.count(ParamId::ACTIVATION_TYPE) == 0) { + if (global_param_ids.size() <= ParamId::ACTIVATION_TYPE) { throw std::runtime_error( "LinearElastanceChamber: activation_type parameter is required"); } From b0aa76c4ceef99862df71022b0c3fbf350a3a18b Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 12 Feb 2026 05:19:40 +0000 Subject: [PATCH 19/23] Make activation_type a required parameter for ChamberElastanceInductor Co-authored-by: aabrown100-git <71424733+aabrown100-git@users.noreply.github.com> --- src/model/ChamberElastanceInductor.cpp | 11 ++++------- src/model/ChamberElastanceInductor.h | 2 +- 2 files changed, 5 insertions(+), 8 deletions(-) diff --git a/src/model/ChamberElastanceInductor.cpp b/src/model/ChamberElastanceInductor.cpp index ca85147a0..a1ff39bc1 100644 --- a/src/model/ChamberElastanceInductor.cpp +++ b/src/model/ChamberElastanceInductor.cpp @@ -55,13 +55,10 @@ void ChamberElastanceInductor::get_elastance_values( void ChamberElastanceInductor::initialize_activation_function( std::vector& parameters) { - // Check if activation_type parameter is provided (optional parameter) - // Defaults to HALF_COSINE if not provided - if (global_param_ids.size() > ParamId::ACTIVATION_TYPE) { - int activation_type_int = static_cast( - parameters[global_param_ids[ParamId::ACTIVATION_TYPE]]); - activation_params_.type = static_cast(activation_type_int); - } + // Get required activation_type parameter + int activation_type_int = static_cast( + parameters[global_param_ids[ParamId::ACTIVATION_TYPE]]); + activation_params_.type = static_cast(activation_type_int); // Set cardiac period activation_params_.cardiac_period = model->cardiac_cycle_period; diff --git a/src/model/ChamberElastanceInductor.h b/src/model/ChamberElastanceInductor.h index 53b56799c..bf69e6ed6 100644 --- a/src/model/ChamberElastanceInductor.h +++ b/src/model/ChamberElastanceInductor.h @@ -155,7 +155,7 @@ class ChamberElastanceInductor : public Block { {"Vrd", InputParameter()}, {"Vrs", InputParameter()}, {"Impedance", InputParameter()}, - {"activation_type", InputParameter(true, false, true, 0.0)}}) {} + {"activation_type", InputParameter()}}) {} /** * @brief Local IDs of the parameters From 0a9dc481babf4e70ba9d9cf48d1a171f60c8d45b Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 12 Feb 2026 05:33:40 +0000 Subject: [PATCH 20/23] Remove backward compatibility for numeric activation_type parameter - only support activation_function_type string Co-authored-by: aabrown100-git <71424733+aabrown100-git@users.noreply.github.com> --- src/solve/SimulationParameters.cpp | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/solve/SimulationParameters.cpp b/src/solve/SimulationParameters.cpp index a91a611c3..ebfbf815c 100644 --- a/src/solve/SimulationParameters.cpp +++ b/src/solve/SimulationParameters.cpp @@ -616,13 +616,6 @@ void create_chambers( " in chamber " + chamber_name); } - // Replace string with numeric value for backward compatibility - chamber_values["activation_type"] = static_cast(act_params.type); - chamber_values.erase("activation_function_type"); - } else if (chamber_values.contains("activation_type")) { - // Handle numeric activation_type for backward compatibility - int act_type_num = chamber_values["activation_type"]; - act_params.type = static_cast(act_type_num); } // Check if activation function values are specified in nested dict From fe33a8000ed8e69b2f80d075e4006b64776eb109 Mon Sep 17 00:00:00 2001 From: aabrown100-git Date: Wed, 11 Feb 2026 21:39:32 -0800 Subject: [PATCH 21/23] Update chamber_elastance_inductor test case with activation function --- tests/cases/chamber_elastance_inductor.json | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/tests/cases/chamber_elastance_inductor.json b/tests/cases/chamber_elastance_inductor.json index 9b756054e..bac17b8ad 100644 --- a/tests/cases/chamber_elastance_inductor.json +++ b/tests/cases/chamber_elastance_inductor.json @@ -101,9 +101,12 @@ "Emin": 0.091, "Vrd": 26.1, "Vrs": 18.0, - "t_active": 0.2, - "t_twitch": 0.3, - "Impedance": 0.000351787 + "Impedance": 0.000351787, + "activation_function_type": "half_cosine", + "activation_function_values": { + "t_active": 0.2, + "t_twitch": 0.3 + } } } ], From ffddd68db91178e2d9a8a6656d07b90b3fae717d Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 12 Feb 2026 05:40:58 +0000 Subject: [PATCH 22/23] Update documentation to reflect mandatory string-based activation_function_type parameter Co-authored-by: aabrown100-git <71424733+aabrown100-git@users.noreply.github.com> --- docs/pages/activation_functions.md | 69 ++++++------------------------ 1 file changed, 14 insertions(+), 55 deletions(-) diff --git a/docs/pages/activation_functions.md b/docs/pages/activation_functions.md index b83c83170..5ddb0c91f 100644 --- a/docs/pages/activation_functions.md +++ b/docs/pages/activation_functions.md @@ -6,9 +6,9 @@ Chamber models in svZeroDSolver support different activation functions to model ## Available Activation Functions -### 1. Half Cosine Activation (Type 0) +### 1. Half Cosine Activation -The default activation function for `ChamberElastanceInductor`. It uses a half cosine wave to model contraction: +A half cosine wave activation function to model contraction: ``` A(t) = -0.5 * cos(2π * t_contract / t_twitch) + 0.5 if t_contract ≤ t_twitch @@ -23,23 +23,6 @@ where `t_contract = max(0, t - t_active)` **Example usage in JSON:** ```json -{ - "type": "ChamberElastanceInductor", - "name": "ventricle", - "values": { - "Emax": 1.057, - "Emin": 0.091, - "Vrd": 26.1, - "Vrs": 18.0, - "t_active": 0.2, - "t_twitch": 0.3, - "Impedance": 0.000351787 - } -} -``` - -Or with explicit activation function specification: -```json { "type": "ChamberElastanceInductor", "name": "ventricle", @@ -58,9 +41,9 @@ Or with explicit activation function specification: } ``` -### 2. Piecewise Cosine Activation (Type 1) +### 2. Piecewise Cosine Activation -The default activation function for `LinearElastanceChamber`. It models separate contraction and relaxation phases: +A piecewise activation function that models separate contraction and relaxation phases: ``` φ(t) = 0.5 * [1 - cos(π * (t - t_C) / T_C)] during contraction @@ -94,7 +77,7 @@ The default activation function for `LinearElastanceChamber`. It models separate } ``` -### 3. Two Hill Activation (Type 2) +### 3. Two Hill Activation A more flexible and physiologically realistic activation function based on the two-hill model. This allows for more precise control over the activation waveform shape. @@ -144,35 +127,19 @@ C is a normalization constant ensuring maximum activation is 1. ## Switching Between Activation Functions -To use a specific activation function, add the `activation_function_type` parameter to your chamber configuration: +**Both `ChamberElastanceInductor` and `LinearElastanceChamber` require the `activation_function_type` parameter to be explicitly specified as a string:** -- `activation_function_type: "half_cosine"` - Half Cosine (default for ChamberElastanceInductor) -- `activation_function_type: "piecewise_cosine"` - Piecewise Cosine (default for LinearElastanceChamber) -- `activation_function_type: "two_hill"` - Two Hill +- `activation_function_type: "half_cosine"` - Half Cosine activation +- `activation_function_type: "piecewise_cosine"` - Piecewise Cosine activation +- `activation_function_type: "two_hill"` - Two Hill activation -Group the activation function-specific parameters under `activation_function_values` as a nested dictionary. +All activation function-specific parameters must be grouped under `activation_function_values` as a nested dictionary. -## Backward Compatibility +### Important Notes -For backward compatibility, the old flat format is also supported: - -```json -{ - "type": "ChamberElastanceInductor", - "name": "ventricle", - "values": { - "Emax": 1.057, - "Emin": 0.091, - "Vrd": 26.1, - "Vrs": 18.0, - "t_active": 0.2, - "t_twitch": 0.3, - "Impedance": 0.000351787 - } -} -``` - -This format will continue to work, with the activation function determined by which parameters are provided. +- The `activation_function_type` parameter is **required** for both chamber types +- Only string values are accepted (e.g., "half_cosine", not numeric values) +- Parameters must be specified within the `activation_function_values` nested dictionary ## Elastance Calculation @@ -193,11 +160,3 @@ E(t) = Epass + Emax * φ(t) The two-hill activation function is described in: - Kaiser, A. D., et al. (2022). "A design-based model of the aortic valve for fluid-structure interaction." Biomechanics and Modeling in Mechanobiology. https://link.springer.com/article/10.1007/s10439-022-03047-3 - -## Backward Compatibility - -Existing configuration files will continue to work without modification: -- `ChamberElastanceInductor` defaults to half cosine activation when `t_active` and `t_twitch` are provided -- `LinearElastanceChamber` defaults to piecewise cosine activation when contraction/relaxation parameters are provided -- The old numeric `activation_type` parameter (0, 1, 2) is still supported but deprecated in favor of the string format -- Parameters can be specified at the top level without using `activation_function_values` for backward compatibility From 28d2faae1c9ea633f64d851e7469ec656a4b17d2 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 12 Feb 2026 05:50:23 +0000 Subject: [PATCH 23/23] Add comprehensive validation for activation_function_type and activation_function_values parameters Co-authored-by: aabrown100-git <71424733+aabrown100-git@users.noreply.github.com> --- src/solve/SimulationParameters.cpp | 140 +++++++++++++++++++++-------- 1 file changed, 102 insertions(+), 38 deletions(-) diff --git a/src/solve/SimulationParameters.cpp b/src/solve/SimulationParameters.cpp index ebfbf815c..dac52629d 100644 --- a/src/solve/SimulationParameters.cpp +++ b/src/solve/SimulationParameters.cpp @@ -596,54 +596,118 @@ void create_chambers( // Prepare activation function parameters ActivationFunctionParams act_params; - act_params.type = ActivationType::HALF_COSINE; // Default act_params.cardiac_period = 0.0; // Will be set later from model - // Check if activation function type is specified as a string - if (chamber_values.contains("activation_function_type") && - chamber_values["activation_function_type"].is_string()) { - std::string act_type_str = chamber_values["activation_function_type"]; - - // Map string to activation type - if (act_type_str == "half_cosine") { - act_params.type = ActivationType::HALF_COSINE; - } else if (act_type_str == "piecewise_cosine") { - act_params.type = ActivationType::PIECEWISE_COSINE; - } else if (act_type_str == "two_hill") { - act_params.type = ActivationType::TWO_HILL; - } else { - throw std::runtime_error("Unknown activation_function_type: " + act_type_str + - " in chamber " + chamber_name); - } - + // Check if activation function type is specified - it's required + if (!chamber_values.contains("activation_function_type")) { + throw std::runtime_error("Missing required parameter 'activation_function_type' in chamber " + + chamber_name + ". Must be one of: 'half_cosine', 'piecewise_cosine', 'two_hill'"); + } + + if (!chamber_values["activation_function_type"].is_string()) { + throw std::runtime_error("Parameter 'activation_function_type' must be a string in chamber " + + chamber_name + ". Must be one of: 'half_cosine', 'piecewise_cosine', 'two_hill'"); + } + + std::string act_type_str = chamber_values["activation_function_type"]; + + // Map string to activation type + if (act_type_str == "half_cosine") { + act_params.type = ActivationType::HALF_COSINE; + } else if (act_type_str == "piecewise_cosine") { + act_params.type = ActivationType::PIECEWISE_COSINE; + } else if (act_type_str == "two_hill") { + act_params.type = ActivationType::TWO_HILL; + } else { + throw std::runtime_error("Unknown activation_function_type: '" + act_type_str + + "' in chamber " + chamber_name + + ". Must be one of: 'half_cosine', 'piecewise_cosine', 'two_hill'"); } // Check if activation function values are specified in nested dict - if (chamber_values.contains("activation_function_values") && - chamber_values["activation_function_values"].is_object()) { - const auto& act_values = chamber_values["activation_function_values"]; + if (!chamber_values.contains("activation_function_values")) { + throw std::runtime_error("Missing required parameter 'activation_function_values' in chamber " + + chamber_name + ". This nested dictionary must contain parameters for the " + + act_type_str + " activation function."); + } + + if (!chamber_values["activation_function_values"].is_object()) { + throw std::runtime_error("Parameter 'activation_function_values' must be an object/dictionary in chamber " + + chamber_name); + } + + const auto& act_values = chamber_values["activation_function_values"]; + std::vector missing_params; + + // Extract and validate activation function parameters based on type + if (act_params.type == ActivationType::HALF_COSINE) { + // Required parameters for half_cosine + if (!act_values.contains("t_active")) missing_params.push_back("t_active"); + if (!act_values.contains("t_twitch")) missing_params.push_back("t_twitch"); + + if (!missing_params.empty()) { + std::string error_msg = "Missing required activation_function_values for 'half_cosine' in chamber " + + chamber_name + ": "; + for (size_t i = 0; i < missing_params.size(); i++) { + error_msg += missing_params[i]; + if (i < missing_params.size() - 1) error_msg += ", "; + } + throw std::runtime_error(error_msg); + } + + act_params.t_active = act_values["t_active"]; + act_params.t_twitch = act_values["t_twitch"]; + + } else if (act_params.type == ActivationType::PIECEWISE_COSINE) { + // Required parameters for piecewise_cosine + if (!act_values.contains("contract_start")) missing_params.push_back("contract_start"); + if (!act_values.contains("relax_start")) missing_params.push_back("relax_start"); + if (!act_values.contains("contract_duration")) missing_params.push_back("contract_duration"); + if (!act_values.contains("relax_duration")) missing_params.push_back("relax_duration"); - // Extract activation function parameters based on type - if (act_params.type == ActivationType::HALF_COSINE) { - if (act_values.contains("t_active")) act_params.t_active = act_values["t_active"]; - if (act_values.contains("t_twitch")) act_params.t_twitch = act_values["t_twitch"]; - } else if (act_params.type == ActivationType::PIECEWISE_COSINE) { - if (act_values.contains("contract_start")) act_params.contract_start = act_values["contract_start"]; - if (act_values.contains("relax_start")) act_params.relax_start = act_values["relax_start"]; - if (act_values.contains("contract_duration")) act_params.contract_duration = act_values["contract_duration"]; - if (act_values.contains("relax_duration")) act_params.relax_duration = act_values["relax_duration"]; - } else if (act_params.type == ActivationType::TWO_HILL) { - if (act_values.contains("t_shift")) act_params.t_shift = act_values["t_shift"]; - if (act_values.contains("tau_1")) act_params.tau_1 = act_values["tau_1"]; - if (act_values.contains("tau_2")) act_params.tau_2 = act_values["tau_2"]; - if (act_values.contains("m1")) act_params.m1 = act_values["m1"]; - if (act_values.contains("m2")) act_params.m2 = act_values["m2"]; + if (!missing_params.empty()) { + std::string error_msg = "Missing required activation_function_values for 'piecewise_cosine' in chamber " + + chamber_name + ": "; + for (size_t i = 0; i < missing_params.size(); i++) { + error_msg += missing_params[i]; + if (i < missing_params.size() - 1) error_msg += ", "; + } + throw std::runtime_error(error_msg); } - // Remove the nested dict (don't flatten to avoid parameter conflicts) - chamber_values.erase("activation_function_values"); + act_params.contract_start = act_values["contract_start"]; + act_params.relax_start = act_values["relax_start"]; + act_params.contract_duration = act_values["contract_duration"]; + act_params.relax_duration = act_values["relax_duration"]; + + } else if (act_params.type == ActivationType::TWO_HILL) { + // Required parameters for two_hill + if (!act_values.contains("t_shift")) missing_params.push_back("t_shift"); + if (!act_values.contains("tau_1")) missing_params.push_back("tau_1"); + if (!act_values.contains("tau_2")) missing_params.push_back("tau_2"); + if (!act_values.contains("m1")) missing_params.push_back("m1"); + if (!act_values.contains("m2")) missing_params.push_back("m2"); + + if (!missing_params.empty()) { + std::string error_msg = "Missing required activation_function_values for 'two_hill' in chamber " + + chamber_name + ": "; + for (size_t i = 0; i < missing_params.size(); i++) { + error_msg += missing_params[i]; + if (i < missing_params.size() - 1) error_msg += ", "; + } + throw std::runtime_error(error_msg); + } + + act_params.t_shift = act_values["t_shift"]; + act_params.tau_1 = act_values["tau_1"]; + act_params.tau_2 = act_values["tau_2"]; + act_params.m1 = act_values["m1"]; + act_params.m2 = act_values["m2"]; } + // Remove the nested dict (don't flatten to avoid parameter conflicts) + chamber_values.erase("activation_function_values"); + // Create the chamber block int block_id = generate_block(model, chamber_values, chamber_type, chamber_name);