diff --git a/cpp/models/smm/simulation.h b/cpp/models/smm/simulation.h index 5acaae3d51..14930ac50c 100644 --- a/cpp/models/smm/simulation.h +++ b/cpp/models/smm/simulation.h @@ -83,19 +83,20 @@ class Simulation update_current_rates_and_waiting_times(); size_t next_event = determine_next_event(); // index of the next event FP current_time = m_result.get_last_time(); - // set in the past to add a new time point immediately - FP last_result_time = current_time - m_dt; + // set current time to add next time point in the future + FP last_result_time = current_time; // iterate over time while (current_time + m_waiting_times[next_event] < tmax) { - // update time - current_time += m_waiting_times[next_event]; - // regularily save current state in m_results - if (current_time > last_result_time + m_dt) { - last_result_time = current_time; - m_result.add_time_point(current_time); + // If the next event happens further in the future than the next stored time point, add a new one. + if (current_time + m_waiting_times[next_event] >= last_result_time) { + auto num_dt = std::ceil((current_time + m_waiting_times[next_event] - last_result_time) / m_dt); + last_result_time = std::min(tmax, last_result_time + num_dt * m_dt); + m_result.add_time_point(last_result_time); // copy from the previous last value m_result.get_last_value() = m_result[m_result.get_num_time_points() - 2]; } + // update time + current_time += m_waiting_times[next_event]; // decide event type by index and perform it if (next_event < adoption_rates().size()) { // perform adoption event