From 7836a6f37fee4566a5c8a4962e3e6448c0548124 Mon Sep 17 00:00:00 2001 From: "charo.heijnen" Date: Tue, 21 Jul 2026 10:27:40 -0400 Subject: [PATCH 1/3] changed Chamber Sphere so that it is possible to select from a variety of activation funcitons. Added double_tanh activation funciton to activation function file. Changed chamver_sphere.json test case to be able to choose activaiton function. when double_tanh activation funciton was selected test case runs as previously ran. --- src/model/ActivationFunction.cpp | 21 ++++++++++++-- src/model/ActivationFunction.h | 32 +++++++++++++++++++++- src/model/ChamberSphere.cpp | 21 ++++---------- src/model/ChamberSphere.h | 44 +++++++++++++++++------------- src/solve/SimulationParameters.cpp | 35 ++++++++++++++++-------- tests/cases/chamber_sphere.json | 8 ++++-- 6 files changed, 110 insertions(+), 51 deletions(-) diff --git a/src/model/ActivationFunction.cpp b/src/model/ActivationFunction.cpp index abc121838..adf3cc74a 100644 --- a/src/model/ActivationFunction.cpp +++ b/src/model/ActivationFunction.cpp @@ -41,6 +41,9 @@ std::unique_ptr ActivationFunction::create_default( if (type_str == "two_hill") { return std::make_unique(cardiac_period); } + if (type_str == "double_tanh") { + return std::make_unique(cardiac_period); + } if (type_str == "fourier") { return std::make_unique(cardiac_period); } @@ -49,8 +52,8 @@ std::unique_ptr ActivationFunction::create_default( } throw std::runtime_error( "Unknown activation_function type '" + type_str + - "'. Must be one of: half_cosine, piecewise_cosine, two_hill, fourier, " - "wrapping_cosine"); + "'. Must be one of: half_cosine, piecewise_cosine, two_hill, " + "double_tanh, fourier, wrapping_cosine"); } double HalfCosineActivation::compute(double time) { @@ -156,6 +159,20 @@ double TwoHillActivation::compute(double time) { return normalization_factor_ * (g1 / (1.0 + g1)) * (1.0 / (1.0 + g2)); } +double DoubleTanhActivation::compute(double time) { + const double tsys = params_.at("tsys"); + const double tdias = params_.at("tdias"); + const double steepness = params_.at("steepness"); + + const double t_in_cycle = std::fmod(time, cardiac_period_); + + const double S_plus = 0.5 * (1.0 + std::tanh((t_in_cycle - tsys) / steepness)); + const double S_minus = + 0.5 * (1.0 - std::tanh((t_in_cycle - tdias) / steepness)); + + return S_plus * S_minus; +} + // ============================================================ // WrappingCosineActivation — atrial activation that wraps // across the cycle boundary (Sankaran 2012, Menon 2023) diff --git a/src/model/ActivationFunction.h b/src/model/ActivationFunction.h index 2c136368b..360a1f64f 100644 --- a/src/model/ActivationFunction.h +++ b/src/model/ActivationFunction.h @@ -59,7 +59,8 @@ class ActivationFunction { /** * @brief Create a default activation function from activation function type * - * @param type_str One of: "half_cosine", "piecewise_cosine", "two_hill" + * @param type_str One of: "half_cosine", "piecewise_cosine", "two_hill", + * "double_tanh", "wrapping_cosine", "fourier" * @param cardiac_period Cardiac cycle period * @return Unique pointer to the created activation function */ @@ -212,6 +213,35 @@ class TwoHillActivation : public ActivationFunction { bool normalization_initialized_; }; +/** + * @brief Double tanh (systole/diastole sigmoid product) activation function + * + * This implements the original ChamberSphere activation: a smooth indicator + * function built from the product of two tanh sigmoids, one rising at + * systole and one falling at diastole. + * + * \f[ + * f(t) = S_+ \cdot S_-, \quad S_\pm = \frac{1}{2} \left(1.0 \pm + * \text{tanh}\left( \frac{t_{in\_cycle} - t_\text{sys/dias}} {\gamma} + * \right) \right) + * \f] + */ +class DoubleTanhActivation : public ActivationFunction { + public: + /** + * @brief Construct with default parameter values (loader fills via + * set_param). + * + * @param cardiac_period Cardiac cycle period + */ + explicit DoubleTanhActivation(double cardiac_period) + : ActivationFunction(cardiac_period, {{"tsys", InputParameter()}, + {"tdias", InputParameter()}, + {"steepness", InputParameter()}}) {} + + double compute(double time) override; +}; + /** * @brief Wrapping cosine activation function * diff --git a/src/model/ChamberSphere.cpp b/src/model/ChamberSphere.cpp index 6971e7713..4d2cb3658 100644 --- a/src/model/ChamberSphere.cpp +++ b/src/model/ChamberSphere.cpp @@ -110,24 +110,15 @@ void ChamberSphere::update_solution( void ChamberSphere::get_elastance_values(std::vector& parameters) { const double alpha_max = parameters[global_param_ids[ParamId::alpha_max]]; const double alpha_min = parameters[global_param_ids[ParamId::alpha_min]]; - const double tsys = parameters[global_param_ids[ParamId::tsys]]; - const double tdias = parameters[global_param_ids[ParamId::tdias]]; - const double steepness = parameters[global_param_ids[ParamId::steepness]]; - const double t = model->time; - - const auto T_cardiac = model->cardiac_cycle_period; - const auto t_in_cycle = fmod(model->time, T_cardiac); - - const double S_plus = 0.5 * (1.0 + tanh((t_in_cycle - tsys) / steepness)); - const double S_minus = 0.5 * (1.0 - tanh((t_in_cycle - tdias) / steepness)); - - // indicator function - const double f = S_plus * S_minus; - - // activation rates + const double f = activation_function_->compute(model->time); const double act_t = alpha_max * f + alpha_min * (1 - f); act = std::abs(act_t); act_plus = std::max(act_t, 0.0); +} + +void ChamberSphere::set_activation_function( + std::unique_ptr af) { + activation_function_ = std::move(af); } \ No newline at end of file diff --git a/src/model/ChamberSphere.h b/src/model/ChamberSphere.h index a0c60cc05..f1173a88c 100644 --- a/src/model/ChamberSphere.h +++ b/src/model/ChamberSphere.h @@ -8,7 +8,11 @@ #define SVZERODSOLVER_MODEL_ChamberSphere_HPP_ #include +#include +#include +#include +#include "ActivationFunction.h" #include "Block.h" #include "SparseSystem.h" @@ -55,11 +59,9 @@ * \dot{\tau} + a \tau - \sigma_\text{max} a_+ = 0, \quad a_+ = \max(a, 0), \quad a = f\alpha_\text{max} + (1 - f)\alpha_\text{min} * \f] - * with indicator function - * \f[ - * f = S_+ \cdot S_-, \quad S_\pm = \frac{1}{2} \left(1.0 \pm \text{tanh}\left( - \frac{t - t_\text{sys/dias}} {\gamma} \right) \right) - * \f] + * where \f$f \in [0, 1]\f$ is the activation function, evaluated by a + * separate \ref ActivationFunction object (e.g. two_hill, half_cosine, + * piecewise_cosine) selected in the JSON configuration. * * 5. Acceleration: * \f[ @@ -89,9 +91,11 @@ * * `sigma_max` - Maximum active stress \f$\sigma_\text{max}\f$ * * `alpha_max` - Maximum activation parameter \f$\alpha_\text{max}\f$ * * `alpha_min` - Minimum activation parameter \f$\alpha_\text{min}\f$ - * * `tsys` - Systole timing parameter \f$t_\text{sys}\f$ - * * `tdias` - Diastole timing parameter \f$t_\text{dias}\f$ - * * `steepness` - Activation steepness parameter \f$\gamma\f$ + * + * An `activation_function` object is also required alongside + * `zero_d_element_values` to select and parameterize the activation function + * \f$f(t)\f$ (see \ref ActivationFunction, e.g. `two_hill`, `half_cosine`, + * `piecewise_cosine`, `wrapping_cosine`, `fourier`, `double_tanh`). * * ### Usage in json configuration file * @@ -111,10 +115,15 @@ * "eta" : 10.0, * "sigma_max" : 185e3, * "alpha_max": 30.0, - * "alpha_min": -30.0, - * "tsys": 0.170, - * "tdias": 0.484, - * "steepness": 0.005 + * "alpha_min": -30.0 + * }, + * "activation_function": { + * "type": "two_hill", + * "t_shift": 0.0, + * "tau_1": 0.25, + * "tau_2": 0.45, + * "m1": 1.5, + * "m2": 8.0 * } * } * ] @@ -146,9 +155,6 @@ class ChamberSphere : public Block { sigma_max = 6, alpha_max = 7, alpha_min = 8, - tsys = 9, - tdias = 10, - steepness = 11 }; /** @@ -167,10 +173,7 @@ class ChamberSphere : public Block { {"eta", InputParameter()}, {"sigma_max", InputParameter()}, {"alpha_max", InputParameter()}, - {"alpha_min", InputParameter()}, - {"tsys", InputParameter()}, - {"tdias", InputParameter()}, - {"steepness", InputParameter()}}) {} + {"alpha_min", InputParameter()}}) {} /** * @brief Set up the degrees of freedom (DOF) of the block @@ -222,9 +225,12 @@ class ChamberSphere : public Block { */ void get_elastance_values(std::vector& parameters); + void set_activation_function(std::unique_ptr af) override; + private: double act = 0.0; // activation function double act_plus = 0.0; // act_plus = max(act, 0) + std::unique_ptr activation_function_; /** * @brief Number of triplets of element diff --git a/src/solve/SimulationParameters.cpp b/src/solve/SimulationParameters.cpp index 036609994..4bcc8e5f9 100644 --- a/src/solve/SimulationParameters.cpp +++ b/src/solve/SimulationParameters.cpp @@ -252,6 +252,19 @@ SimulationParameters load_simulation_params(const nlohmann::json& config) { void load_simulation_model(const nlohmann::json& config, Model& model) { DEBUG_MSG("Loading model"); + + // Set cardiac period from simulation_parameters so activation functions + // have it available while blocks are created below. May already be set by + // closed_loop_blocks. + if (model.cardiac_cycle_period < 0.0 && + config.contains("simulation_parameters") && + config["simulation_parameters"].contains("cardiac_period")) { + double period = config["simulation_parameters"]["cardiac_period"]; + if (period > 0.0) { + model.cardiac_cycle_period = period; + } + } + // Create list to store block connections while generating blocks std::vector> connections; @@ -335,10 +348,18 @@ void create_vessels( JsonWrapper(config, component, "vessel_name", i); const auto& vessel_values = vessel_config["zero_d_element_values"]; const std::string vessel_name = vessel_config["vessel_name"]; + const std::string vessel_type = vessel_config["zero_d_element_type"]; vessel_id_map.insert({vessel_config["vessel_id"], vessel_name}); - generate_block(model, vessel_values, vessel_config["zero_d_element_type"], - vessel_name); + generate_block(model, vessel_values, vessel_type, vessel_name); + + // Create and set activation_function for vessel types that use one + if (vessel_type == "ChamberSphere") { + auto act_func = generate_activation_function( + model, vessel_config["activation_function"], vessel_name); + model.get_block(vessel_name) + ->set_activation_function(std::move(act_func)); + } // Read connected boundary conditions if (vessel_config.contains("boundary_conditions")) { @@ -599,16 +620,6 @@ void create_chambers( Model& model, std::vector>& connections, const nlohmann::json& config, const std::string& component) { - // Set cardiac period from simulation_parameters so activation functions have - // it. May already be set by closed_loop_blocks. - if (model.cardiac_cycle_period < 0.0 && - config.contains("simulation_parameters") && - config["simulation_parameters"].contains("cardiac_period")) { - double period = config["simulation_parameters"]["cardiac_period"]; - if (period > 0.0) { - model.cardiac_cycle_period = period; - } - } for (size_t i = 0; i < config[component].size(); i++) { const auto& chamber_config = JsonWrapper(config, component, "name", i); std::string chamber_type = chamber_config["type"]; diff --git a/tests/cases/chamber_sphere.json b/tests/cases/chamber_sphere.json index 4936ccf09..b287c738d 100644 --- a/tests/cases/chamber_sphere.json +++ b/tests/cases/chamber_sphere.json @@ -26,7 +26,8 @@ "steady_initial": false, "output_variable_based": true, "absolute_tolerance": 1e-9, - "output_all_cycles": true + "output_all_cycles": true, + "cardiac_period": 1.0 }, "vessels": [ { @@ -57,7 +58,10 @@ "eta": 10.0, "sigma_max": 185e3, "alpha_max": 30.0, - "alpha_min": -30.0, + "alpha_min": -30.0 + }, + "activation_function": { + "type": "double_tanh", "tsys": 0.17, "tdias": 0.484, "steepness": 0.005 From ec94a3911fdfc2263a436f7108ef65d561ea6c2a Mon Sep 17 00:00:00 2001 From: "charo.heijnen" Date: Fri, 24 Jul 2026 11:27:26 -0400 Subject: [PATCH 2/3] Change in documentation activationfunction.h --- src/model/ActivationFunction.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/model/ActivationFunction.h b/src/model/ActivationFunction.h index 360a1f64f..2f4e1b3b6 100644 --- a/src/model/ActivationFunction.h +++ b/src/model/ActivationFunction.h @@ -221,8 +221,8 @@ class TwoHillActivation : public ActivationFunction { * systole and one falling at diastole. * * \f[ - * f(t) = S_+ \cdot S_-, \quad S_\pm = \frac{1}{2} \left(1.0 \pm - * \text{tanh}\left( \frac{t_{in\_cycle} - t_\text{sys/dias}} {\gamma} + * f(t) = S_+ (t_{in\_cycle} - t_sys) \cdot S_- (t_{in\_cycle} - t_dias), \quad S_\pm (\Delta t) = \frac{1}{2} \left(1.0 \pm + * \text{tanh}\left( \frac{\Delta t} {\gamma} * \right) \right) * \f] */ From 11a484193455b08ed10060a668954af562082ffb Mon Sep 17 00:00:00 2001 From: "charo.heijnen" Date: Fri, 24 Jul 2026 14:23:06 -0400 Subject: [PATCH 3/3] node version change --- .github/workflows/gui.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/gui.yml b/.github/workflows/gui.yml index 9309c2eb2..7ac0e36f8 100644 --- a/.github/workflows/gui.yml +++ b/.github/workflows/gui.yml @@ -19,7 +19,7 @@ jobs: - name: Set up Node.js uses: actions/setup-node@v3 with: - node-version: '20' + node-version: '24' - name: Install dependencies working-directory: tests/cypress