Skip to content

Commit 24ff34a

Browse files
aabrown100-gitaabrown100-git
committed
Initial commit. Implements ActivationFunction with half cosine, piecewise cosine, and two hill activation functions. ActivationFunction is created in SimulationParameters.cpp and associated with a chamber. Renamed PiecewiseCosineChamber to LinearElastanceChamber. All tests pass
1 parent d40189e commit 24ff34a

17 files changed

Lines changed: 578 additions & 86 deletions

‎.gitignore‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -38,3 +38,7 @@ build*/
3838

3939
# Node modules (for directed graph visualization)
4040
node_modules/
41+
42+
# Virtual environments
43+
venv/
44+
.venv/

‎src/model/ActivationFunction.cpp‎

Lines changed: 149 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,149 @@
1+
// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the
2+
// University of California, and others. SPDX-License-Identifier: BSD-3-Clause
3+
4+
#include "ActivationFunction.h"
5+
6+
#include <cmath>
7+
#include <stdexcept>
8+
9+
ActivationFunction::ActivationFunction(
10+
double cardiac_period,
11+
std::vector<std::pair<std::string, InputParameter>> params)
12+
: cardiac_period_(cardiac_period) {
13+
for (auto& p : params) {
14+
if (p.second.is_number) {
15+
double default_val =
16+
p.second.is_optional ? p.second.default_val : 0.0;
17+
params_[p.first] = {std::move(p.second), default_val};
18+
}
19+
}
20+
}
21+
22+
void ActivationFunction::set_param(const std::string& name, double value) {
23+
auto it = params_.find(name);
24+
if (it == params_.end()) {
25+
throw std::runtime_error(
26+
"ActivationFunction::set_param: unknown parameter '" + name + "'");
27+
}
28+
if (!it->second.first.is_number) {
29+
throw std::runtime_error(
30+
"ActivationFunction::set_param: parameter '" + name +
31+
"' is not a scalar number");
32+
}
33+
it->second.second = value;
34+
}
35+
36+
std::vector<std::pair<std::string, InputParameter>>
37+
ActivationFunction::get_input_params() const {
38+
std::vector<std::pair<std::string, InputParameter>> out;
39+
for (const auto& p : params_) {
40+
out.push_back({p.first, p.second.first});
41+
}
42+
return out;
43+
}
44+
45+
std::unique_ptr<ActivationFunction> ActivationFunction::create_default(
46+
const std::string& type_str, double cardiac_period) {
47+
if (type_str == "half_cosine") {
48+
return std::make_unique<HalfCosineActivation>(cardiac_period);
49+
}
50+
if (type_str == "piecewise_cosine") {
51+
return std::make_unique<PiecewiseCosineActivation>(cardiac_period);
52+
}
53+
if (type_str == "two_hill") {
54+
return std::make_unique<TwoHillActivation>(cardiac_period);
55+
}
56+
throw std::runtime_error(
57+
"Unknown activation_function type '" + type_str +
58+
"'. Must be one of: half_cosine, piecewise_cosine, two_hill");
59+
}
60+
61+
double HalfCosineActivation::compute(double time) {
62+
double t_in_cycle = std::fmod(time, cardiac_period_);
63+
const double t_active = params_.at("t_active").second;
64+
const double t_twitch = params_.at("t_twitch").second;
65+
66+
double t_contract = 0.0;
67+
if (t_in_cycle >= t_active) {
68+
t_contract = t_in_cycle - t_active;
69+
}
70+
71+
double act = 0.0;
72+
if (t_contract <= t_twitch) {
73+
act = -0.5 * std::cos(2.0 * M_PI * t_contract / t_twitch) + 0.5;
74+
}
75+
76+
return act;
77+
}
78+
79+
double PiecewiseCosineActivation::compute(double time) {
80+
const double contract_start = params_.at("contract_start").second;
81+
const double relax_start = params_.at("relax_start").second;
82+
const double contract_duration = params_.at("contract_duration").second;
83+
const double relax_duration = params_.at("relax_duration").second;
84+
85+
double phi = 0.0;
86+
double piecewise_condition =
87+
std::fmod(time - contract_start, cardiac_period_);
88+
89+
if (0.0 <= piecewise_condition &&
90+
piecewise_condition < contract_duration) {
91+
phi = 0.5 * (1.0 - std::cos((M_PI * piecewise_condition) /
92+
contract_duration));
93+
} else {
94+
piecewise_condition = std::fmod(time - relax_start, cardiac_period_);
95+
if (0.0 <= piecewise_condition &&
96+
piecewise_condition < relax_duration) {
97+
phi = 0.5 * (1.0 + std::cos((M_PI * piecewise_condition) /
98+
relax_duration));
99+
}
100+
}
101+
102+
return phi;
103+
}
104+
105+
void TwoHillActivation::calculate_normalization_factor() {
106+
const double tau_1 = params_.at("tau_1").second;
107+
const double tau_2 = params_.at("tau_2").second;
108+
const double m1 = params_.at("m1").second;
109+
const double m2 = params_.at("m2").second;
110+
111+
constexpr double NORMALIZATION_DT = 1e-5;
112+
double max_value = 0.0;
113+
114+
for (double t_temp = 0.0; t_temp < cardiac_period_; t_temp += NORMALIZATION_DT) {
115+
double g1 = std::pow(t_temp / tau_1, m1);
116+
double g2 = std::pow(t_temp / tau_2, m2);
117+
double two_hill_val = (g1 / (1.0 + g1)) * (1.0 / (1.0 + g2));
118+
max_value = std::max(max_value, two_hill_val);
119+
}
120+
121+
normalization_factor_ = 1.0 / max_value;
122+
normalization_initialized_ = true;
123+
}
124+
125+
void TwoHillActivation::finalize() {
126+
calculate_normalization_factor();
127+
}
128+
129+
double TwoHillActivation::compute(double time) {
130+
if (!normalization_initialized_) {
131+
throw std::runtime_error(
132+
"TwoHillActivation: call finalize() after setting parameters");
133+
}
134+
135+
const double t_shift = params_.at("t_shift").second;
136+
const double tau_1 = params_.at("tau_1").second;
137+
const double tau_2 = params_.at("tau_2").second;
138+
const double m1 = params_.at("m1").second;
139+
const double m2 = params_.at("m2").second;
140+
141+
double t_in_cycle = std::fmod(time, cardiac_period_);
142+
double t_shifted = std::fmod(t_in_cycle - t_shift, cardiac_period_);
143+
t_shifted = (t_shifted >= 0.0) ? t_shifted : t_shifted + cardiac_period_;
144+
145+
double g1 = std::pow(t_shifted / tau_1, m1);
146+
double g2 = std::pow(t_shifted / tau_2, m2);
147+
148+
return normalization_factor_ * (g1 / (1.0 + g1)) * (1.0 / (1.0 + g2));
149+
}

‎src/model/ActivationFunction.h‎

Lines changed: 208 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,208 @@
1+
// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the
2+
// University of California, and others. SPDX-License-Identifier: BSD-3-Clause
3+
4+
/**
5+
* @file ActivationFunction.h
6+
* @brief Activation function classes for cardiac chamber models
7+
*/
8+
9+
#ifndef SVZERODSOLVER_MODEL_ACTIVATIONFUNCTION_HPP_
10+
#define SVZERODSOLVER_MODEL_ACTIVATIONFUNCTION_HPP_
11+
12+
#include <map>
13+
#include <memory>
14+
#include <string>
15+
#include <vector>
16+
17+
#include "Parameter.h"
18+
19+
/**
20+
* @brief Base class for activation functions
21+
*
22+
* Activation functions compute the activation value (between 0 and 1) at a
23+
* given time point within a cardiac cycle. These are used to modulate
24+
* chamber elastance over time.
25+
*/
26+
class ActivationFunction {
27+
public:
28+
/**
29+
* @brief Construct activation function
30+
*
31+
* @param cardiac_period Cardiac cycle period
32+
* @param params Parameter definitions (name, InputParameter) for this activation function
33+
*/
34+
ActivationFunction(
35+
double cardiac_period,
36+
std::vector<std::pair<std::string, InputParameter>> params);
37+
38+
/**
39+
* @brief Virtual destructor
40+
*/
41+
virtual ~ActivationFunction() = default;
42+
43+
/**
44+
* @brief Compute activation value at given time
45+
*
46+
* @param time Current time
47+
* @return Activation value between 0 and 1
48+
*/
49+
virtual double compute(double time) = 0;
50+
51+
/**
52+
* @brief Create a default activation function from activation function type
53+
*
54+
* @param type_str One of: "half_cosine", "piecewise_cosine", "two_hill"
55+
* @param cardiac_period Cardiac cycle period
56+
* @return Unique pointer to the created activation function
57+
*/
58+
static std::unique_ptr<ActivationFunction> create_default(
59+
const std::string& type_str, double cardiac_period);
60+
61+
/**
62+
* @brief Set a scalar parameter value by name.
63+
*
64+
* Validates that name is in params_ and has a number schema, then stores
65+
* the value. No per-class logic needed.
66+
*
67+
* @param name Parameter name (must be a key in params_)
68+
* @param value Parameter value
69+
*/
70+
void set_param(const std::string& name, double value);
71+
72+
/**
73+
* @brief Parameter definitions for validation/loading (analogous to Block::input_params).
74+
*
75+
* Returns (name, InputParameter) for each parameter. Built from params_.
76+
*/
77+
std::vector<std::pair<std::string, InputParameter>> get_input_params() const;
78+
79+
/**
80+
* @brief Called after all parameters are set (e.g. by loader).
81+
*
82+
* Default no-op. TwoHillActivation overrides to recompute normalization.
83+
*/
84+
virtual void finalize() {}
85+
86+
protected:
87+
double cardiac_period_;
88+
std::map<std::string, std::pair<InputParameter, double>> params_;
89+
};
90+
91+
/**
92+
* @brief Half cosine activation function
93+
*
94+
* This implements the activation function used in the original
95+
* ChamberElastanceInductor. The activation follows a half cosine wave
96+
* during the contraction period.
97+
*
98+
* \f[
99+
* A(t) = \begin{cases}
100+
* -\frac{1}{2}\cos(2\pi t_{contract}/t_{twitch}) + \frac{1}{2}, & \text{if } t_{contract} \le t_{twitch} \\
101+
* 0, & \text{otherwise}
102+
* \end{cases}
103+
* \f]
104+
*
105+
* where \f$t_{contract} = \max(0, t_{in\_cycle} - t_{active})\f$
106+
*/
107+
class HalfCosineActivation : public ActivationFunction {
108+
public:
109+
/**
110+
* @brief Construct with default parameter values (loader fills via set_param).
111+
*
112+
* @param cardiac_period Cardiac cycle period
113+
*/
114+
explicit HalfCosineActivation(double cardiac_period)
115+
: ActivationFunction(cardiac_period, {{"t_active", InputParameter()},
116+
{"t_twitch", InputParameter()}}) {}
117+
118+
double compute(double time) override;
119+
};
120+
121+
122+
/**
123+
* @brief Piecewise cosine activation function
124+
*
125+
* This implements the activation function from the LinearElastanceChamber
126+
* (Regazzoni chamber model). The activation consists of separate contraction
127+
* and relaxation phases, each following a cosine curve.
128+
*
129+
* \f[
130+
* \phi(t, t_C, t_R, T_C, T_R) = \begin{cases}
131+
* \frac{1}{2}\left[1 - \cos\left(\frac{\pi}{T_C} \bmod(t - t_C, T_{HB})\right)\right],
132+
* & \text{if } 0 \le \bmod(t - t_C, T_{HB}) < T_C \\
133+
* \frac{1}{2}\left[1 + \cos\left(\frac{\pi}{T_R} \bmod(t - t_R, T_{HB})\right)\right],
134+
* & \text{if } 0 \le \bmod(t - t_R, T_{HB}) < T_R \\
135+
* 0, & \text{otherwise}
136+
* \end{cases}
137+
* \f]
138+
*/
139+
class PiecewiseCosineActivation : public ActivationFunction {
140+
public:
141+
/**
142+
* @brief Construct with default parameter values (loader fills via set_param).
143+
*
144+
* @param cardiac_period Cardiac cycle period
145+
*/
146+
explicit PiecewiseCosineActivation(double cardiac_period)
147+
: ActivationFunction(cardiac_period,
148+
{{"contract_start", InputParameter()},
149+
{"relax_start", InputParameter()},
150+
{"contract_duration", InputParameter()},
151+
{"relax_duration", InputParameter()}}) {}
152+
153+
double compute(double time) override;
154+
};
155+
156+
/**
157+
* @brief Two hill activation function
158+
*
159+
* This implements the two-hill activation function which provides more
160+
* flexible and physiologically realistic waveforms. See
161+
* https://link.springer.com/article/10.1007/s10439-022-03047-3
162+
*
163+
* The activation is computed as:
164+
* \f[
165+
* A(t) = C \cdot \frac{g_1(t)}{1 + g_1(t)} \cdot \frac{1}{1 + g_2(t)}
166+
* \f]
167+
*
168+
* where:
169+
* \f[
170+
* g_1(t) = \left(\frac{t_{shifted}}{\tau_1}\right)^{m_1}, \quad
171+
* g_2(t) = \left(\frac{t_{shifted}}{\tau_2}\right)^{m_2}
172+
* \f]
173+
*
174+
* and \f$t_{shifted} = (t - t_{shift}) \bmod T_{cardiac}\f$, and \f$C\f$ is a
175+
* normalization constant to ensure max activation is 1.
176+
*/
177+
class TwoHillActivation : public ActivationFunction {
178+
public:
179+
/**
180+
* @brief Construct with default parameter values (loader fills via set_param).
181+
*
182+
* Defaults for tau_1, tau_2, m1, m2 are 1.0 to avoid division by zero.
183+
* Call finalize() after all set_param to recompute normalization.
184+
*
185+
* @param cardiac_period Cardiac cycle period
186+
*/
187+
explicit TwoHillActivation(double cardiac_period)
188+
: ActivationFunction(cardiac_period,
189+
{{"t_shift", InputParameter()},
190+
{"tau_1", InputParameter(false, false, true, 1.0)},
191+
{"tau_2", InputParameter(false, false, true, 1.0)},
192+
{"m1", InputParameter(false, false, true, 1.0)},
193+
{"m2", InputParameter(false, false, true, 1.0)}}),
194+
normalization_factor_(1.0),
195+
normalization_initialized_(false) {}
196+
197+
double compute(double time) override;
198+
199+
void finalize() override;
200+
201+
private:
202+
void calculate_normalization_factor();
203+
204+
double normalization_factor_;
205+
bool normalization_initialized_;
206+
};
207+
208+
#endif // SVZERODSOLVER_MODEL_ACTIVATIONFUNCTION_HPP_

‎src/model/Block.h‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,7 @@
1313
#include <map>
1414
#include <vector>
1515

16+
#include "ActivationFunction.h"
1617
#include "BlockType.h"
1718
#include "DOFHandler.h"
1819
#include "Parameter.h"
@@ -283,6 +284,17 @@ class Block {
283284
* @return TripletsContributions Number of triplets of element
284285
*/
285286
virtual TripletsContributions get_num_triplets();
287+
288+
/**
289+
* @brief Set activation function (for chamber blocks that use one).
290+
*
291+
* Default no-op. Overridden by ChamberElastanceInductor and
292+
* LinearElastanceChamber to take ownership of the activation function.
293+
*
294+
* @param af Unique pointer to the activation function (caller transfers ownership)
295+
*/
296+
virtual void set_activation_function(
297+
std::unique_ptr<ActivationFunction> /*af*/) {}
286298
};
287299

288300
#endif

0 commit comments

Comments
 (0)