|
| 1 | +/* |
| 2 | +* Copyright (C) 2020-2025 MEmilio |
| 3 | +* |
| 4 | +* Authors: Henrik Zunker |
| 5 | +* |
| 6 | +* Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de> |
| 7 | +* |
| 8 | +* Licensed under the Apache License, Version 2.0 (the "License"); |
| 9 | +* you may not use this file except in compliance with the License. |
| 10 | +* You may obtain a copy of the License at |
| 11 | +* |
| 12 | +* http://www.apache.org/licenses/LICENSE-2.0 |
| 13 | +* |
| 14 | +* Unless required by applicable law or agreed to in writing, software |
| 15 | +* distributed under the License is distributed on an "AS IS" BASIS, |
| 16 | +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| 17 | +* See the License for the specific language governing permissions and |
| 18 | +* limitations under the License. |
| 19 | +*/ |
| 20 | +#include "ode_mseirs4/model.h" |
| 21 | +#include "memilio/compartments/simulation.h" |
| 22 | +#include <iostream> |
| 23 | + |
| 24 | +int main() |
| 25 | +{ |
| 26 | + using FP = double; |
| 27 | + mio::omseirs4::Model<FP> model; |
| 28 | + auto& params = model.parameters; |
| 29 | + |
| 30 | + // Example parameter values for day-based time unit (t in days) |
| 31 | + params.get<mio::omseirs4::BaseTransmissionRate<FP>>() = 0.4; // b0 per day |
| 32 | + params.get<mio::omseirs4::SeasonalAmplitude<FP>>() = 0.15; // b1 |
| 33 | + params.get<mio::omseirs4::SeasonalPhase<FP>>() = 0.0; // phi; for phase shift use 2*pi*offsetDays/365 |
| 34 | + params.get<mio::omseirs4::NaturalBirthDeathRate<FP>>() = 1.0 / (70.0 * 365.0); // mu per day |
| 35 | + params.get<mio::omseirs4::LossMaternalImmunityRate<FP>>() = 1.0 / 90.0; // xi per day (~3 months) |
| 36 | + params.get<mio::omseirs4::ProgressionRate<FP>>() = 1.0 / 7.0; // sigma per day (≈7 days latent) |
| 37 | + params.get<mio::omseirs4::RecoveryRate<FP>>() = 1.0 / 14.0; // nu per day (≈14 days infectious) |
| 38 | + params.get<mio::omseirs4::ImmunityWaningRate<FP>>() = 1.0 / (5.0 * 365.0); // gamma per day (5 years) |
| 39 | + // factors default already set (0.5, 0.35, 0.25) as in the paper https://doi.org/10.1016/S0025-5564(01)00066-9. |
| 40 | + |
| 41 | + // Initial population (absolute counts), set all compartments explicitly. |
| 42 | + double N = 1e6; |
| 43 | + |
| 44 | + // Initial populations. |
| 45 | + model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::MaternalImmune)}] = |
| 46 | + 5000.0; |
| 47 | + |
| 48 | + model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::E1)}] = 300.0; |
| 49 | + model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::E2)}] = 150.0; |
| 50 | + model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::E3)}] = 80.0; |
| 51 | + model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::E4)}] = 70.0; |
| 52 | + |
| 53 | + model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::I1)}] = 200.0; |
| 54 | + model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::I2)}] = 100.0; |
| 55 | + model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::I3)}] = 50.0; |
| 56 | + model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::I4)}] = 50.0; |
| 57 | + |
| 58 | + model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::R1)}] = 40000.0; |
| 59 | + model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::R2)}] = 30000.0; |
| 60 | + model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::R3)}] = 20000.0; |
| 61 | + model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::R4)}] = 10000.0; |
| 62 | + |
| 63 | + model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::S2)}] = 100000.0; |
| 64 | + model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::S3)}] = 50000.0; |
| 65 | + model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::S4)}] = 50000.0; |
| 66 | + |
| 67 | + // Compute S1 as residual to match N |
| 68 | + double assigned = 5000.0 + (300.0 + 150.0 + 80.0 + 70.0) + (200.0 + 100.0 + 50.0 + 50.0) + |
| 69 | + (40000.0 + 30000.0 + 20000.0 + 10000.0) + (100000.0 + 50000.0 + 50000.0); |
| 70 | + double S1 = N - assigned; |
| 71 | + if (S1 < 0) |
| 72 | + S1 = 0; |
| 73 | + model.populations[{mio::Index<mio::omseirs4::InfectionState>(mio::omseirs4::InfectionState::S1)}] = S1; |
| 74 | + |
| 75 | + model.check_constraints(); |
| 76 | + |
| 77 | + // simulate |
| 78 | + double t0 = 0.0; |
| 79 | + double tmax = 20.0; // days |
| 80 | + double dt = 1.0; // daily output |
| 81 | + auto result = mio::simulate(t0, tmax, dt, model); |
| 82 | + |
| 83 | + // print header |
| 84 | + std::cout << "t M S1 S2 S3 S4 E1 E2 E3 E4 I1 I2 I3 I4 R1 R2 R3 R4\n"; |
| 85 | + for (size_t i = 0; i < (size_t)result.get_num_time_points(); ++i) { |
| 86 | + std::cout << result.get_time(i); |
| 87 | + const auto& y = result.get_value(i); |
| 88 | + for (size_t k = 0; k < (size_t)mio::omseirs4::InfectionState::Count; ++k) { |
| 89 | + std::cout << ' ' << y[(Eigen::Index)k]; |
| 90 | + } |
| 91 | + std::cout << '\n'; |
| 92 | + } |
| 93 | + return 0; |
| 94 | +} |
0 commit comments