Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

1133 Change member functions in Model class from public to private in IDE model #1166

Merged
merged 2 commits into from
Jan 6, 2025
Merged

1133 Change member functions in Model class from public to private in IDE model #1166

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ed31e7a
make simulation functions private and define friend classes and funct…
annawendler Dec 19, 2024
1abd2b4
review
annawendler Jan 6, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
174 changes: 115 additions & 59 deletions cpp/models/ide_secir/model.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -67,6 +67,99 @@ Model::Model(TimeSeries<ScalarType>&& init, CustomIndexArray<ScalarType, AgeGrou
}
}

bool Model::check_constraints(ScalarType dt) const
{

if (!((size_t)m_transitions.get_num_elements() == (size_t)InfectionTransition::Count * m_num_agegroups)) {
log_error("A variable given for model construction is not valid. Number of elements in transition vector "
"does not match the required number.");
return true;
}

for (AgeGroup group = AgeGroup(0); group < AgeGroup(m_num_agegroups); ++group) {

for (int i = 0; i < (int)InfectionState::Count; i++) {
int index = get_state_flat_index(i, group);
if (m_populations[0][index] < 0) {
log_error("Initialization failed. Initial values for populations are less than zero.");
return true;
}
}
}

// It may be possible to run the simulation with fewer time points, but this number ensures that it is possible.
if (m_transitions.get_num_time_points() < (Eigen::Index)std::ceil(get_global_support_max(dt) / dt)) {
log_error("Initialization failed. Not enough time points for transitions given before start of "
"simulation.");
return true;
}

for (AgeGroup group = AgeGroup(0); group < AgeGroup(m_num_agegroups); ++group) {

for (int i = 0; i < m_transitions.get_num_time_points(); i++) {
for (int j = 0; j < (int)InfectionTransition::Count; j++) {
int index = get_transition_flat_index(j, group);
if (m_transitions[i][index] < 0) {
log_error("Initialization failed. One or more initial value for transitions is less than zero.");
return true;
}
}
}
}
if (m_transitions.get_last_time() != m_populations.get_last_time()) {
log_error("Last time point of TimeSeries for transitions does not match last time point of "
"TimeSeries for "
"compartments. Both of these time points have to agree for a sensible simulation.");
return true;
}

if (m_populations.get_num_time_points() != 1) {
log_error("The TimeSeries for the compartments contains more than one time point. It is unclear how to "
"initialize.");
return true;
}

return parameters.check_constraints();
}

// Note that this function computes the global_support_max via the get_support_max() method and does not make use
// of the vector m_transitiondistributions_support_max. This is because the global_support_max is already used in
// check_constraints and we cannot ensure that the vector has already been computed when checking for constraints
// (which usually happens before setting the initial flows and simulating).
ScalarType Model::get_global_support_max(ScalarType dt) const
{
ScalarType global_support_max = 0.;
ScalarType global_support_max_new = 0.;
for (AgeGroup group = AgeGroup(0); group < AgeGroup(m_num_agegroups); ++group) {
global_support_max_new = std::max(
{parameters.get<TransitionDistributions>()[group][(int)InfectionTransition::ExposedToInfectedNoSymptoms]
.get_support_max(dt, m_tol),
parameters
.get<TransitionDistributions>()[group][(int)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms]
.get_support_max(dt, m_tol),
parameters.get<TransitionDistributions>()[group][(int)InfectionTransition::InfectedNoSymptomsToRecovered]
.get_support_max(dt, m_tol),
parameters
.get<TransitionDistributions>()[group][(int)InfectionTransition::InfectedSymptomsToInfectedSevere]
.get_support_max(dt, m_tol),
parameters.get<TransitionDistributions>()[group][(int)InfectionTransition::InfectedSymptomsToRecovered]
.get_support_max(dt, m_tol),
parameters
.get<TransitionDistributions>()[group][(int)InfectionTransition::InfectedSevereToInfectedCritical]
.get_support_max(dt, m_tol),
parameters.get<TransitionDistributions>()[group][(int)InfectionTransition::InfectedSevereToRecovered]
.get_support_max(dt, m_tol),
parameters.get<TransitionDistributions>()[group][(int)InfectionTransition::InfectedCriticalToDead]
.get_support_max(dt, m_tol),
parameters.get<TransitionDistributions>()[group][(int)InfectionTransition::InfectedCriticalToRecovered]
.get_support_max(dt, m_tol)});
if (global_support_max_new > global_support_max) {
global_support_max = global_support_max_new;
}
}
return global_support_max;
}

// ---- Functionality to calculate the sizes of the compartments for time t0. ----
void Model::compute_compartment_from_flows(ScalarType dt, Eigen::Index idx_InfectionState, AgeGroup group,
Eigen::Index idx_IncomingFlow, int idx_TransitionDistribution1,
Expand Down Expand Up @@ -384,7 +477,26 @@ void Model::flows_current_timestep(ScalarType dt)
compute_flow(Eigen::Index(InfectionTransition::InfectedCriticalToRecovered),
Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), dt, group);
}
} // namespace isecir
}

void Model::update_compartment_from_flow(InfectionState infectionState,
std::vector<InfectionTransition> const& IncomingFlows,
std::vector<InfectionTransition> const& OutgoingFlows, AgeGroup group)
{
int state_idx = get_state_flat_index(Eigen::Index(infectionState), group);

Eigen::Index num_time_points = m_populations.get_num_time_points();
ScalarType updated_compartment = m_populations[num_time_points - 2][state_idx];
for (const InfectionTransition& inflow : IncomingFlows) {
int inflow_idx = get_transition_flat_index(Eigen::Index(inflow), group);
updated_compartment += m_transitions.get_last_value()[inflow_idx];
}
for (const InfectionTransition& outflow : OutgoingFlows) {
int outflow_idx = get_transition_flat_index(Eigen::Index(outflow), group);
updated_compartment -= m_transitions.get_last_value()[outflow_idx];
}
m_populations.get_last_value()[state_idx] = updated_compartment;
}

void Model::update_compartments()
{
Expand Down Expand Up @@ -434,25 +546,6 @@ void Model::update_compartments()
}
}

void Model::update_compartment_from_flow(InfectionState infectionState,
std::vector<InfectionTransition> const& IncomingFlows,
std::vector<InfectionTransition> const& OutgoingFlows, AgeGroup group)
{
int state_idx = get_state_flat_index(Eigen::Index(infectionState), group);

Eigen::Index num_time_points = m_populations.get_num_time_points();
ScalarType updated_compartment = m_populations[num_time_points - 2][state_idx];
for (const InfectionTransition& inflow : IncomingFlows) {
int inflow_idx = get_transition_flat_index(Eigen::Index(inflow), group);
updated_compartment += m_transitions.get_last_value()[inflow_idx];
}
for (const InfectionTransition& outflow : OutgoingFlows) {
int outflow_idx = get_transition_flat_index(Eigen::Index(outflow), group);
updated_compartment -= m_transitions.get_last_value()[outflow_idx];
}
m_populations.get_last_value()[state_idx] = updated_compartment;
}

void Model::compute_forceofinfection(ScalarType dt, bool initialization)
{

Expand Down Expand Up @@ -532,8 +625,9 @@ void Model::compute_forceofinfection(ScalarType dt, bool initialization)
m_forceofinfection[i] += divNj * sum;
}
}
} // namespace mio
}

// ---- Functionality to set vectors with necessary information regarding TransitionDistributions. ----
void Model::set_transitiondistributions_support_max(ScalarType dt)
{
m_transitiondistributions_support_max = CustomIndexArray<std::vector<ScalarType>, AgeGroup>(
Expand Down Expand Up @@ -618,43 +712,5 @@ void Model::set_transitiondistributions_in_forceofinfection(ScalarType dt)
}
}

// Note that this function computes the global_support_max via the get_support_max() method and does not make use
// of the vector m_transitiondistributions_support_max. This is because the global_support_max is already used in
// check_constraints and we cannot ensure that the vector has already been computed when checking for constraints
// (which usually happens before setting the initial flows and simulating).
ScalarType Model::get_global_support_max(ScalarType dt) const
{
ScalarType global_support_max = 0.;
ScalarType global_support_max_new = 0.;
for (AgeGroup group = AgeGroup(0); group < AgeGroup(m_num_agegroups); ++group) {
global_support_max_new = std::max(
{parameters.get<TransitionDistributions>()[group][(int)InfectionTransition::ExposedToInfectedNoSymptoms]
.get_support_max(dt, m_tol),
parameters
.get<TransitionDistributions>()[group][(int)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms]
.get_support_max(dt, m_tol),
parameters.get<TransitionDistributions>()[group][(int)InfectionTransition::InfectedNoSymptomsToRecovered]
.get_support_max(dt, m_tol),
parameters
.get<TransitionDistributions>()[group][(int)InfectionTransition::InfectedSymptomsToInfectedSevere]
.get_support_max(dt, m_tol),
parameters.get<TransitionDistributions>()[group][(int)InfectionTransition::InfectedSymptomsToRecovered]
.get_support_max(dt, m_tol),
parameters
.get<TransitionDistributions>()[group][(int)InfectionTransition::InfectedSevereToInfectedCritical]
.get_support_max(dt, m_tol),
parameters.get<TransitionDistributions>()[group][(int)InfectionTransition::InfectedSevereToRecovered]
.get_support_max(dt, m_tol),
parameters.get<TransitionDistributions>()[group][(int)InfectionTransition::InfectedCriticalToDead]
.get_support_max(dt, m_tol),
parameters.get<TransitionDistributions>()[group][(int)InfectionTransition::InfectedCriticalToRecovered]
.get_support_max(dt, m_tol)});
if (global_support_max_new > global_support_max) {
global_support_max = global_support_max_new;
}
}
return global_support_max;
}

} // namespace isecir
} // namespace mio
Loading
Loading