Skip to content

Commit 3ed1a68

Browse files
jubickerreneSchm
jubicker
and
reneSchm
authored
1520 Operate only on current system state and not on result time series in SMM (#1521)
- Evaluate the current rates in the SMM dependent on the model populations instead of the last result value - If an event occurs, only modify the model populations and not the current value Co-authored-by: reneSchm <[email protected]>
1 parent 7b60e83 commit 3ed1a68

2 files changed

Lines changed: 28 additions & 30 deletions

File tree

‎cpp/models/smm/simulation.h‎

Lines changed: 17 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -83,35 +83,30 @@ class Simulation
8383
update_current_rates_and_waiting_times();
8484
size_t next_event = determine_next_event(); // index of the next event
8585
FP current_time = m_result.get_last_time();
86-
// set current time to add next time point in the future
87-
FP last_result_time = current_time;
86+
// set next result time; system state is saved every dt time step
87+
FP next_result_time = current_time + m_dt;
8888
// iterate over time
8989
while (current_time + m_waiting_times[next_event] < tmax) {
90-
// If the next event happens further in the future than the next stored time point, add a new one.
91-
if (current_time + m_waiting_times[next_event] >= last_result_time) {
92-
auto num_dt = std::ceil((current_time + m_waiting_times[next_event] - last_result_time) / m_dt);
93-
last_result_time = std::min(tmax, last_result_time + num_dt * m_dt);
94-
m_result.add_time_point(last_result_time);
95-
// copy from the previous last value
96-
m_result.get_last_value() = m_result[m_result.get_num_time_points() - 2];
97-
}
9890
// update time
9991
current_time += m_waiting_times[next_event];
92+
// Save results until the next event time point
93+
// do not save current time, as it does not yet include the next event
94+
while (next_result_time < current_time) {
95+
m_result.add_time_point(next_result_time);
96+
m_result.get_last_value() = m_model->populations.get_compartments();
97+
next_result_time += m_dt;
98+
}
10099
// decide event type by index and perform it
101100
if (next_event < adoption_rates().size()) {
102101
// perform adoption event
103102
const auto& rate = adoption_rates()[next_event];
104-
m_result.get_last_value()[m_model->populations.get_flat_index({rate.region, rate.from})] -= 1;
105103
m_model->populations[{rate.region, rate.from}] -= 1.0;
106-
m_result.get_last_value()[m_model->populations.get_flat_index({rate.region, rate.to})] += 1;
107104
m_model->populations[{rate.region, rate.to}] += 1.0;
108105
}
109106
else {
110107
// perform transition event
111108
const auto& rate = transition_rates()[next_event - adoption_rates().size()];
112-
m_result.get_last_value()[m_model->populations.get_flat_index({rate.from, rate.status})] -= 1;
113109
m_model->populations[{rate.from, rate.status}] -= 1.0;
114-
m_result.get_last_value()[m_model->populations.get_flat_index({rate.to, rate.status})] += 1;
115110
m_model->populations[{rate.to, rate.status}] += 1.0;
116111
}
117112
// update internal times
@@ -124,10 +119,13 @@ class Simulation
124119
update_current_rates_and_waiting_times();
125120
next_event = determine_next_event();
126121
}
127-
// copy last result, if no event occurs between last_result_time and tmax
128-
if (last_result_time < tmax) {
129-
m_result.add_time_point(tmax);
130-
m_result.get_last_value() = m_result[m_result.get_num_time_points() - 2];
122+
// copy last result, if no event occurs between last result time point and tmax
123+
if (m_result.get_last_time() < tmax) {
124+
while (next_result_time <= tmax) {
125+
m_result.add_time_point(next_result_time);
126+
m_result.get_last_value() = m_model->populations.get_compartments();
127+
next_result_time += m_dt;
128+
}
131129
// update internal times
132130
for (size_t i = 0; i < m_internal_time.size(); i++) {
133131
m_internal_time[i] += m_current_rates[i] * (tmax - current_time);
@@ -184,7 +182,7 @@ class Simulation
184182
inline void update_current_rates_and_waiting_times()
185183
{
186184
size_t i = 0; // shared index for iterating both rates
187-
const auto last_values = m_result.get_last_value();
185+
const auto last_values = m_model->populations.get_compartments();
188186
for (const auto& rate : adoption_rates()) {
189187
m_current_rates[i] = m_model->evaluate(rate, last_values);
190188
m_waiting_times[i] = (m_current_rates[i] > 0)

‎cpp/tests/test_smm_model.cpp‎

Lines changed: 11 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -214,27 +214,27 @@ TEST(TestSMMSimulation, advance)
214214
.WillOnce(Return(0.0031)) //spatial transition event 1->0, corresponds to global time 0.31
215215
.WillRepeatedly(testing::Return(1.0));
216216

217-
auto sim = mio::smm::Simulation(model, 0.0, 0.1);
217+
auto sim = mio::smm::Simulation(model, 0.0, 0.05);
218218
sim.advance(30.);
219219
//Check whether first value in result time series corresponds to initial values
220220
EXPECT_EQ(sim.get_result().get_value(0)[static_cast<size_t>(InfectionState::S)], 1);
221221
EXPECT_EQ(sim.get_result().get_value(0)[static_cast<size_t>(InfectionState::I)], 1);
222222
EXPECT_EQ(sim.get_result().get_value(
223223
0)[static_cast<size_t>(InfectionState::Count) + static_cast<size_t>(InfectionState::R)],
224224
1);
225-
//no event happens in first two time steps i.e. the first saved result is after t=0.2
226-
EXPECT_GE(sim.get_result().get_time(1), 0.2);
225+
//first event happens at 0.2005 which means there are no events in the first four time steps
226+
EXPECT_GE(sim.get_result().get_time(4), 0.2);
227227
//adoption from I to R is first event
228-
EXPECT_EQ(sim.get_result().get_value(1)[static_cast<size_t>(InfectionState::S)], 1);
229-
EXPECT_EQ(sim.get_result().get_value(1)[static_cast<size_t>(InfectionState::I)], 0);
230-
EXPECT_EQ(sim.get_result().get_value(1)[static_cast<size_t>(InfectionState::R)], 1);
228+
EXPECT_EQ(sim.get_result().get_value(5)[static_cast<size_t>(InfectionState::S)], 1);
229+
EXPECT_EQ(sim.get_result().get_value(5)[static_cast<size_t>(InfectionState::I)], 0);
230+
EXPECT_EQ(sim.get_result().get_value(5)[static_cast<size_t>(InfectionState::R)], 1);
231231
EXPECT_EQ(sim.get_result().get_value(
232232
1)[static_cast<size_t>(InfectionState::Count) + static_cast<size_t>(InfectionState::R)],
233233
1);
234234
//spatial transition is second event
235-
EXPECT_EQ(sim.get_result().get_value(2)[static_cast<size_t>(InfectionState::S)], 1);
236-
EXPECT_EQ(sim.get_result().get_value(2)[static_cast<size_t>(InfectionState::I)], 0);
237-
EXPECT_EQ(sim.get_result().get_value(2)[static_cast<size_t>(InfectionState::R)], 2);
235+
EXPECT_EQ(sim.get_result().get_value(7)[static_cast<size_t>(InfectionState::S)], 1);
236+
EXPECT_EQ(sim.get_result().get_value(7)[static_cast<size_t>(InfectionState::I)], 0);
237+
EXPECT_EQ(sim.get_result().get_value(7)[static_cast<size_t>(InfectionState::R)], 2);
238238
}
239239

240240
TEST(TestSMMSimulation, stopsAtTmax)
@@ -262,9 +262,9 @@ TEST(TestSMMSimulation, stopsAtTmax)
262262

263263
//As populations are not set they have value 0 i.e. no events will happen
264264
//Simulate for 30 days
265-
auto sim = mio::smm::Simulation(model, 0.0, 0.1);
265+
auto sim = mio::smm::Simulation(model, 0.0, 30.);
266266
sim.advance(30.);
267-
//As model populations are all zero only t0 and tmax should be logged
267+
//As dt=30 only t0 and tmax should be logged
268268
EXPECT_EQ(sim.get_result().get_num_time_points(), 2);
269269
EXPECT_EQ(sim.get_result().get_last_time(), 30.);
270270
}

0 commit comments

Comments
 (0)