Two hill activation function and separate activation function class

.gitignore

@@ -38,3 +38,7 @@ build*/
 
 # Node modules (for directed graph visualization)
 node_modules/
+
+# Virtual environments
+venv/
+.venv/

src/model/ActivationFunction.cpp

@@ -0,0 +1,150 @@
+// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the
+// University of California, and others. SPDX-License-Identifier: BSD-3-Clause
+
+#include "ActivationFunction.h"
+
+#include <cmath>
+#include <stdexcept>
+#include <string>
+
+#ifndef M_PI
+#define M_PI 3.14159265358979323846
+#endif
+
+ActivationFunction::ActivationFunction(
+    double cardiac_period,
+    const std::vector<std::pair<std::string, InputParameter>>&
+        input_param_properties)
+    : cardiac_period_(cardiac_period),
+      input_param_properties(input_param_properties) {
+  // Initialize parameters with default values
+  for (const auto& p : this->input_param_properties) {
+    if (p.second.is_number) {
+      double default_val = p.second.is_optional ? p.second.default_val : 0.0;
+      params_[p.first] = default_val;
+    }
+  }
+}
+
+void ActivationFunction::set_param(const std::string& name, double value) {
+  params_[name] = value;
+}
+
+std::unique_ptr<ActivationFunction> ActivationFunction::create_default(
+    const std::string& type_str, double cardiac_period) {
+  if (type_str == "half_cosine") {
+    return std::make_unique<HalfCosineActivation>(cardiac_period);
+  }
+  if (type_str == "piecewise_cosine") {
+    return std::make_unique<PiecewiseCosineActivation>(cardiac_period);
+  }
+  if (type_str == "two_hill") {
+    return std::make_unique<TwoHillActivation>(cardiac_period);
+  }
+  throw std::runtime_error(
+      "Unknown activation_function type '" + type_str +
+      "'. Must be one of: half_cosine, piecewise_cosine, two_hill");
+}
+
+double HalfCosineActivation::compute(double time) {
+  double t_in_cycle = std::fmod(time, cardiac_period_);
+  const double t_active = params_.at("t_active");
+  const double t_twitch = params_.at("t_twitch");
+
+  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;
+}
+
+double PiecewiseCosineActivation::compute(double time) {
+  const double contract_start = params_.at("contract_start");
+  const double relax_start = params_.at("relax_start");
+  const double contract_duration = params_.at("contract_duration");
+  const double relax_duration = params_.at("relax_duration");
+
+  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;
+}
+
+void TwoHillActivation::calculate_normalization_factor() {
+  if (cardiac_period_ <= 0.0) {
+    throw std::runtime_error(
+        "TwoHillActivation::calculate_normalization_factor: cardiac_period "
+        "must be positive (got " +
+        std::to_string(cardiac_period_) + ")");
+  }
+
+  const double tau_1 = params_.at("tau_1");
+  const double tau_2 = params_.at("tau_2");
+  const double m1 = params_.at("m1");
+  const double m2 = params_.at("m2");
+
+  constexpr double NORMALIZATION_DT = 1e-5;
+  double max_value = 0.0;
+
+  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));
+    max_value = std::max(max_value, two_hill_val);
+  }
+
+  if (!(max_value > 0.0) || !std::isfinite(max_value)) {
+    throw std::runtime_error(
+        "TwoHillActivation::calculate_normalization_factor: max activation "
+        "value must be positive and finite (got " +
+        std::to_string(max_value) +
+        "). Check tau_1, tau_2, m1, m2 are valid (e.g., tau_1 > 0, tau_2 > "
+        "0).");
+  }
+
+  normalization_factor_ = 1.0 / max_value;
+  normalization_initialized_ = true;
+}
+
+void TwoHillActivation::finalize() { calculate_normalization_factor(); }
+
+double TwoHillActivation::compute(double time) {
+  if (!normalization_initialized_) {
+    throw std::runtime_error(
+        "TwoHillActivation: call finalize() after setting parameters");
+  }
+
+  const double t_shift = params_.at("t_shift");
+  const double tau_1 = params_.at("tau_1");
+  const double tau_2 = params_.at("tau_2");
+  const double m1 = params_.at("m1");
+  const double m2 = params_.at("m2");
+
+  double t_in_cycle = std::fmod(time, 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_;
+
+  double g1 = std::pow(t_shifted / tau_1, m1);
+  double g2 = std::pow(t_shifted / tau_2, m2);
+
+  return normalization_factor_ * (g1 / (1.0 + g1)) * (1.0 / (1.0 + g2));
+}

src/model/ActivationFunction.h

@@ -0,0 +1,215 @@
+// 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 <map>
+#include <memory>
+#include <string>
+#include <vector>
+
+#include "Parameter.h"
+
+/**
+ * @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 Properties of the input parameters for this activation function
+   * [(name, InputParameter), ...]
+   */
+  const std::vector<std::pair<std::string, InputParameter>>
+      input_param_properties;
+
+  /**
+   * @brief Construct activation function
+   *
+   * @param cardiac_period Cardiac cycle period
+   * @param input_param_properties Properties of the input parameters
+   * [(name, InputParameter), ...] for this activation function
+   */
+  ActivationFunction(double cardiac_period,
+                     const std::vector<std::pair<std::string, InputParameter>>&
+                         input_param_properties);
+
+  /**
+   * @brief Virtual destructor
+   */
+  virtual ~ActivationFunction() = default;
+
+  /**
+   * @brief Compute activation value at given time
+   *
+   * @param time Current time
+   * @return Activation value between 0 and 1
+   */
+  virtual double compute(double time) = 0;
+
+  /**
+   * @brief Create a default activation function from activation function type
+   *
+   * @param type_str One of: "half_cosine", "piecewise_cosine", "two_hill"
+   * @param cardiac_period Cardiac cycle period
+   * @return Unique pointer to the created activation function
+   */
+  static std::unique_ptr<ActivationFunction> create_default(
+      const std::string& type_str, double cardiac_period);
+
+  /**
+   * @brief Set a scalar parameter value by name.
+   *
+   * Calling function must validate the parameter name and value
+   *
+   * @param name Parameter name
+   * @param value Parameter value
+   */
+  void set_param(const std::string& name, double value);
+
+  /**
+   * @brief Called after all parameters are set (e.g. by loader).
+   *
+   * Default no-op. TwoHillActivation overrides to recompute normalization.
+   */
+  virtual void finalize() {}
+
+ protected:
+  /**
+   * @brief Time duration of one cardiac cycle
+   */
+  double cardiac_period_;
+
+  /**
+   * @brief Map of parameter names to their values
+   */
+  std::map<std::string, double> params_;
+};
+
+/**
+ * @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 with default parameter values (loader fills via
+   * set_param).
+   *
+   * @param cardiac_period Cardiac cycle period
+   */
+  explicit HalfCosineActivation(double cardiac_period)
+      : ActivationFunction(cardiac_period, {{"t_active", InputParameter()},
+                                            {"t_twitch", InputParameter()}}) {}
+
+  double compute(double time) override;
+};
+
+/**
+ * @brief Piecewise cosine activation function
+ *
+ * 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.
+ *
+ * \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 with default parameter values (loader fills via
+   * set_param).
+   *
+   * @param cardiac_period Cardiac cycle period
+   */
+  explicit PiecewiseCosineActivation(double cardiac_period)
+      : ActivationFunction(cardiac_period,
+                           {{"contract_start", InputParameter()},
+                            {"relax_start", InputParameter()},
+                            {"contract_duration", InputParameter()},
+                            {"relax_duration", InputParameter()}}) {}
+
+  double compute(double time) override;
+};
+
+/**
+ * @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 with default parameter values (loader fills via
+   * set_param).
+   *
+   * @param cardiac_period Cardiac cycle period
+   */
+  explicit TwoHillActivation(double cardiac_period)
+      : ActivationFunction(cardiac_period, {{"t_shift", InputParameter()},
+                                            {"tau_1", InputParameter()},
+                                            {"tau_2", InputParameter()},
+                                            {"m1", InputParameter()},
+                                            {"m2", InputParameter()}}),
+        normalization_factor_(1.0),
+        normalization_initialized_(false) {}
+
+  double compute(double time) override;
+
+  void finalize() override;
+
+ private:
+  void calculate_normalization_factor();
+
+  double normalization_factor_;
+  bool normalization_initialized_;
+};
+
+#endif  // SVZERODSOLVER_MODEL_ACTIVATIONFUNCTION_HPP_