Add option to select activation function for ChamberSphere

src/model/ActivationFunction.cpp

@@ -41,6 +41,9 @@ std::unique_ptr<ActivationFunction> ActivationFunction::create_default(
   if (type_str == "two_hill") {
     return std::make_unique<TwoHillActivation>(cardiac_period);
   }
+  if (type_str == "double_tanh") {
+    return std::make_unique<DoubleTanhActivation>(cardiac_period);
+  }
   if (type_str == "fourier") {
     return std::make_unique<FourierActivation>(cardiac_period);
   }
@@ -49,8 +52,8 @@ std::unique_ptr<ActivationFunction> 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)

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_+ (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]
+ */
+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
  *

src/model/ChamberSphere.cpp

@@ -110,31 +110,15 @@ void ChamberSphere::update_solution(
 void ChamberSphere::get_elastance_values(std::vector<double>& 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);
-
-  auto warp_signed = [T_cardiac](double dt) {
-    return fmod(dt + 1.5 * T_cardiac, T_cardiac) - 0.5 * T_cardiac;
-  };
-
-  const double phase_tsys = warp_signed(t_in_cycle - tsys);
-  const double phase_tdias = warp_signed(t_in_cycle - tdias);
-
-  const double S_plus = 0.5 * (1.0 + tanh(phase_tsys / steepness));
-  const double S_minus = 0.5 * (1.0 - tanh(phase_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<ActivationFunction> af) {
+  activation_function_ = std::move(af);
 }

src/model/ChamberSphere.h

@@ -8,7 +8,11 @@
 #define SVZERODSOLVER_MODEL_ChamberSphere_HPP_
 
 #include <math.h>
+#include <map>              
+#include <memory>           
+#include <string> 
 
+#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
@@ -225,9 +228,12 @@ class ChamberSphere : public Block {
    */
   void get_elastance_values(std::vector<double>& parameters);
 
+  void set_activation_function(std::unique_ptr<ActivationFunction> af) override;
+
  private:
   double act = 0.0;       // activation function
   double act_plus = 0.0;  // act_plus = max(act, 0)
+  std::unique_ptr<ActivationFunction> activation_function_;
 
   /**
    * @brief Number of triplets of element

src/solve/SimulationParameters.cpp

@@ -145,14 +145,14 @@ std::unique_ptr<ActivationFunction> generate_activation_function(
     throw std::runtime_error(
         "Missing 'activation_function' for chamber " + chamber_name +
         ". Required with structure: {\"type\": \"half_cosine\", \"t_active\": "
-        "0.2, \"t_twitch\": 0.3} (or type piecewise_cosine / two_hill with "
+        "0.2, \"t_twitch\": 0.3} (or type piecewise_cosine / two_hill / double_tanh with "
         "their parameters).");
   }
   if (!j.contains("type") || !j["type"].is_string()) {
     throw std::runtime_error(
         "Missing or invalid 'type' in activation_function for chamber " +
         chamber_name +
-        ". Must be one of: half_cosine, piecewise_cosine, two_hill");
+        ". Must be one of: half_cosine, piecewise_cosine, two_hill, double_tanh");
   }
 
   // Extract activation function type
@@ -254,6 +254,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<std::tuple<std::string, std::string>> connections;
 
@@ -337,10 +350,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")) {
@@ -601,16 +622,6 @@ void create_chambers(
     Model& model,
     std::vector<std::tuple<std::string, std::string>>& 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"];

tests/cases/chamber_sphere.json

@@ -27,6 +27,7 @@
     "output_variable_based": true,
     "absolute_tolerance": 1e-9,
     "output_all_cycles": true,
+    "cardiac_period": 1.0,
     "max_iter_error_to_warning": false
   },
   "vessels": [
@@ -58,7 +59,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