SciCompMod · HenrZu · Feb 6, 2024 · Dec 2, 2022 · Dec 9, 2022 · Dec 9, 2022
diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt
@@ -21,7 +21,6 @@ target_compile_options(ode_seir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNI
 add_executable(ode_sir_example ode_sir.cpp)
 target_link_libraries(ode_sir_example PRIVATE memilio ode_sir)
 
-
 add_executable(seir_flows_example ode_seir_flows.cpp)
 target_link_libraries(seir_flows_example PRIVATE memilio ode_seir)
 

diff --git a/cpp/examples/ode_secir.cpp b/cpp/examples/ode_secir.cpp
@@ -76,12 +76,20 @@ int main()
 
     model.apply_constraints();
 
+    // Using default Integrator
+    mio::TimeSeries<double> secir = simulate(t0, tmax, dt, model);
+
+    /*
+    Example of using a different integrator
+   All available integrators are listed in cpp/memilio/math/README.md
+
     auto integrator = std::make_shared<mio::RKIntegratorCore>();
     integrator->set_dt_min(0.3);
     integrator->set_dt_max(1.0);
     integrator->set_rel_tolerance(1e-4);
     integrator->set_abs_tolerance(1e-1);
     mio::TimeSeries<double> secir = simulate(t0, tmax, dt, model, integrator);
+    */
 
     bool print_to_terminal = true;
 

diff --git a/cpp/examples/ode_secir_ageres.cpp b/cpp/examples/ode_secir_ageres.cpp
@@ -50,16 +50,16 @@ int main()
 
     auto& params = model.parameters;
 
-    params.set<mio::osecir::ICUCapacity>(std::numeric_limits<double>::max());
-    params.set<mio::osecir::StartDay>(0);
-    params.set<mio::osecir::Seasonality>(0);
+    params.set<mio::osecir::StartDay>(60);
+    params.set<mio::osecir::Seasonality>(0.2);
+    params.get<mio::osecir::TestAndTraceCapacity>() = 35;
 
     for (auto i = mio::AgeGroup(0); i < nb_groups; i++) {
         params.get<mio::osecir::IncubationTime>()[i]       = 5.2;
-        params.get<mio::osecir::TimeInfectedSymptoms>()[i] = 6.;
+        params.get<mio::osecir::TimeInfectedSymptoms>()[i] = 5.8;
         params.get<mio::osecir::SerialInterval>()[i]       = 4.2;
-        params.get<mio::osecir::TimeInfectedSevere>()[i]   = 12;
-        params.get<mio::osecir::TimeInfectedCritical>()[i] = 8;
+        params.get<mio::osecir::TimeInfectedSevere>()[i]   = 9.5;
+        params.get<mio::osecir::TimeInfectedCritical>()[i] = 7.1;
 
         model.populations[{i, mio::osecir::InfectionState::Exposed}]                     = fact * nb_exp_t0;
         model.populations[{i, mio::osecir::InfectionState::InfectedNoSymptoms}]          = fact * nb_car_t0;
@@ -73,38 +73,44 @@ int main()
         model.populations.set_difference_from_group_total<mio::AgeGroup>({i, mio::osecir::InfectionState::Susceptible},
                                                                          fact * nb_total_t0);
 
-        params.get<mio::osecir::TransmissionProbabilityOnContact>()[i] = 0.05;
-        params.get<mio::osecir::RelativeTransmissionNoSymptoms>()[i]   = 0.67;
-        params.get<mio::osecir::RecoveredPerInfectedNoSymptoms>()[i]   = 0.09;
-        params.get<mio::osecir::RiskOfInfectionFromSymptomatic>()[i]   = 0.25;
-        params.get<mio::osecir::SeverePerInfectedSymptoms>()[i]        = 0.2;
-        params.get<mio::osecir::CriticalPerSevere>()[i]                = 0.25;
-        params.get<mio::osecir::DeathsPerCritical>()[i]                = 0.3;
+        params.get<mio::osecir::TransmissionProbabilityOnContact>()[i]  = 0.05;
+        params.get<mio::osecir::RelativeTransmissionNoSymptoms>()[i]    = 0.7;
+        params.get<mio::osecir::RecoveredPerInfectedNoSymptoms>()[i]    = 0.09;
+        params.get<mio::osecir::RiskOfInfectionFromSymptomatic>()[i]    = 0.25;
+        params.get<mio::osecir::MaxRiskOfInfectionFromSymptomatic>()[i] = 0.45;
+        params.get<mio::osecir::SeverePerInfectedSymptoms>()[i]         = 0.2;
+        params.get<mio::osecir::CriticalPerSevere>()[i]                 = 0.25;
+        params.get<mio::osecir::DeathsPerCritical>()[i]                 = 0.3;
     }
 
+    model.apply_constraints();
+
     mio::ContactMatrixGroup& contact_matrix = params.get<mio::osecir::ContactPatterns>();
     contact_matrix[0] =
         mio::ContactMatrix(Eigen::MatrixXd::Constant((size_t)nb_groups, (size_t)nb_groups, fact * cont_freq));
     contact_matrix.add_damping(Eigen::MatrixXd::Constant((size_t)nb_groups, (size_t)nb_groups, 0.7),
                                mio::SimulationTime(30.));
 
-    model.apply_constraints();
-
     mio::TimeSeries<double> secir = simulate(t0, tmax, dt, model);
 
-    std::vector<std::string> vars = {"S", "E", "C", "C_confirmed", "I", "I_confirmed", "H", "U", "R", "D"};
-    printf("Number of time points :%d\n", static_cast<int>(secir.get_num_time_points()));
-    printf("People in\n");
+    bool print_to_terminal = true;
 
-    for (size_t k = 0; k < (size_t)mio::osecir::InfectionState::Count; k++) {
-        double dummy = 0;
+    if (print_to_terminal) {
 
-        for (size_t i = 0; i < (size_t)params.get_num_groups(); i++) {
-            printf("\t %s[%d]: %.0f", vars[k].c_str(), (int)i,
-                   secir.get_last_value()[k + (size_t)mio::osecir::InfectionState::Count * (int)i]);
-            dummy += secir.get_last_value()[k + (size_t)mio::osecir::InfectionState::Count * (int)i];
-        }
+        std::vector<std::string> vars = {"S", "E", "C", "C_confirmed", "I", "I_confirmed", "H", "U", "R", "D"};
+        printf("Number of time points :%d\n", static_cast<int>(secir.get_num_time_points()));
+        printf("People in\n");
+
+        for (size_t k = 0; k < (size_t)mio::osecir::InfectionState::Count; k++) {
+            double dummy = 0;
 
-        printf("\t %s_otal: %.0f\n", vars[k].c_str(), dummy);
+            for (size_t i = 0; i < (size_t)params.get_num_groups(); i++) {
+                printf("\t %s[%d]: %.0f", vars[k].c_str(), (int)i,
+                    secir.get_last_value()[k + (size_t)mio::osecir::InfectionState::Count * (int)i]);
+                dummy += secir.get_last_value()[k + (size_t)mio::osecir::InfectionState::Count * (int)i];
+            }
+
+            printf("\t %s_otal: %.0f\n", vars[k].c_str(), dummy);
+        }
     }
 }
diff --git a/cpp/examples/ode_seir.cpp b/cpp/examples/ode_seir.cpp
@@ -28,34 +28,40 @@ int main()
     mio::set_log_level(mio::LogLevel::debug);
 
     double t0   = 0;
-    double tmax = 1;
-    double dt   = 0.001;
+    double tmax = 50.;
+    double dt   = 0.1;
 
-    mio::log_info("Simulating SEIR; t={} ... {} with dt = {}.", t0, tmax, dt);
+    mio::log_info("Simulating ODE SEIR; t={} ... {} with dt = {}.", t0, tmax, dt);
 
     mio::oseir::Model model;
 
-    double total_population                                                                            = 10000;
-    model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Exposed)}]   = 100;
-    model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Infected)}]  = 100;
-    model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Recovered)}] = 100;
+    double total_population                                                                            = 1061000;
+    model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Exposed)}]   = 10000;
+    model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Infected)}]  = 1000;
+    model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Recovered)}] = 1000;
     model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Susceptible)}] =
         total_population -
         model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Exposed)}] -
         model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Infected)}] -
         model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Recovered)}];
-    // suscetible now set with every other update
-    // params.nb_sus_t0   = params.nb_total_t0 - params.nb_exp_t0 - params.nb_inf_t0 - params.nb_rec_t0;
-    model.parameters.set<mio::oseir::TimeExposed>(5.2);
-    model.parameters.set<mio::oseir::TimeInfected>(6);
+
     model.parameters.set<mio::oseir::TransmissionProbabilityOnContact>(0.04);
-    model.parameters.get<mio::oseir::ContactPatterns>().get_baseline()(0, 0) = 10;
+    model.parameters.set<mio::oseir::TimeExposed>(5.2);
+    model.parameters.set<mio::oseir::TimeInfected>(2);
+
+    model.parameters.get<mio::oseir::ContactPatterns>().get_baseline()(0, 0) = 2.7;
+    model.parameters.get<mio::oseir::ContactPatterns>().add_damping(0.6, mio::SimulationTime(12.5));
 
     model.check_constraints();
-    // print_seir_params(model);
 
-    auto seir = simulate(t0, tmax, dt, model);
+    mio::TimeSeries<double> seir = simulate(t0, tmax, dt, model);
+
+    bool print_to_terminal = true;
+
+    if (print_to_terminal) {
+        seir.print_table({"S", "E", "I", "R"});
 
-    printf("\n number total: %f\n",
-           seir.get_last_value()[0] + seir.get_last_value()[1] + seir.get_last_value()[2] + seir.get_last_value()[3]);
+        Eigen::VectorXd res_j = seir.get_last_value();
+        printf("\nnumber total: %f\n", res_j[0] + res_j[1] + res_j[2] + res_j[3]);
+    }
 }
diff --git a/cpp/examples/ode_sir.cpp b/cpp/examples/ode_sir.cpp
@@ -30,7 +30,7 @@ int main()
 
     double t0   = 0.;
     double tmax = 50.;
-    double dt   = 0.1002004008016032;
+    double dt   = 0.1;
 
     double total_population = 1061000;
 
@@ -45,7 +45,7 @@ int main()
         model.populations[{mio::Index<mio::osir::InfectionState>(mio::osir::InfectionState::Infected)}] -
         model.populations[{mio::Index<mio::osir::InfectionState>(mio::osir::InfectionState::Recovered)}];
     model.parameters.set<mio::osir::TimeInfected>(2);
-    model.parameters.set<mio::osir::TransmissionProbabilityOnContact>(1);
+    model.parameters.set<mio::osir::TransmissionProbabilityOnContact>(0.04);
     model.parameters.get<mio::osir::ContactPatterns>().get_baseline()(0, 0) = 2.7;
     model.parameters.get<mio::osir::ContactPatterns>().add_damping(0.6, mio::SimulationTime(12.5));
 

diff --git a/cpp/tests/data/ode-secihurd-ageres-compare.csv b/cpp/tests/data/ode-secihurd-ageres-compare.csv
@@ -0,0 +1,37 @@
+ # t S E C I H U R D
+0.00000000000000  9760.00000000000000 99.99999999999999 49.99999999999999 49.99999999999999 20.00000000000000 10.00000000000000 10.00000000000000 0.00000000000000
+0.10000000000000  9757.42416799591410 99.45942737057143 50.60117953460856 51.41459534037426 19.96451694367820 9.91235721061906 11.18168764569956 0.04206795853284
+0.55000000000000  9745.36796156551645 97.66667068776195 52.80317079869250 57.77716750828912 19.86873991622822 9.53203585185988 16.75736344521052 0.22689022644118
+1.42115996873272  9720.21765090973895 96.45016368613835 55.52314466738508 69.71226872067037 19.96186679690890 8.86043024423564 28.70933468271379 0.56514029220753
+2.33084745755960  9691.92497869655745 97.24907014366903 57.32956232056675 81.23625092962342 20.40215502749755 8.24788282351597 42.71644458897381 0.89365546959368
+3.36568976301034  9657.67036545653718 99.69315091310736 59.03242578073674 93.06520105587174 21.25033245683768 7.65656176296593 60.39098866328460 1.24097391065658
+4.54471398597572  9616.26545410955441 103.66470935699886 61.07602911913177 105.10496285204607 22.56760663135830 7.11114693206236 82.60181454739303 1.60827645145411
+5.91864691579154  9565.01269242785384 109.23983456621171 63.86573782398520 117.71690271212263 24.45445337661126 6.63300432235815 111.07092792625593 2.00644684459844
+7.57080600739402  9499.19410816582058 116.74870764248212 67.82411588230345 131.63822819160845 27.08530103801698 6.25645008704367 148.79791363713463 2.45517535558465
+9.69771495700291  9407.74307882261564 127.17291599480350 73.67468778267948 148.64104161943021 30.88485223841903 6.04179128982807 202.83591199883398 3.00572025338731
+12.08190364281839  9296.22452465780225 139.40684417028925 80.88564407650992 167.46159621257237 35.54957306972640 6.09401637298649 270.76310550157035 3.61469593853811
+14.61257442846575  9167.47575404596682 152.67226650290124 88.94997334310170 187.73811597049439 40.87331291140326 6.41469191936146 351.59455070711243 4.28133459965344
+17.14324521411311  9028.18674950193599 165.96117773034928 97.19250101362846 208.47519095166331 46.52382925600715 6.94737206775424 441.71910300301539 4.99407647564145
+19.82441300844518  8869.48653938637290 179.78757071614550 105.91416208766753 230.84318467044426 52.81822443851458 7.69140829181014 547.63716741492431 5.82174299411559
+22.50558080277724  8699.97241830218991 193.05536641291368 114.43150055055563 253.35949078445978 59.37503625763978 8.57981152931013 664.48415894668506 6.74221721624076
+25.30653394080261  8512.34585123302168 205.96572333385399 122.89715587156614 276.65157340485166 66.43207932010536 9.62746749807299 798.26156830851005 7.81858103001262
+28.10748707882798  8315.35625857618652 217.50715139376189 130.68467876930453 299.22996543796091 73.60777346949087 10.76780833299896 943.82167042326114 9.02469359702818
+28.66767770643306  8274.98768218145415 219.61356715329242 132.13821335488299 303.61273232383508 75.04747944156549 11.00462905592607 974.31333081450384 9.28236567453371
+28.97740381328433  8252.54405679252704 220.74593564993978 132.92505833098369 306.01255424024566 75.84323126146357 11.13662505255828 991.36529228931340 9.42724638296136
+29.25615730945047  8232.95277181661731 221.07772161590265 133.60950480142310 308.15700261585653 76.55906996753524 11.25603004959417 1006.82877941063930 9.55911972242535
+29.66087406925001  8212.50000875450860 213.90199553300408 134.13116596969613 311.21408310190498 77.59742185917916 11.43036501733785 1029.47186492954302 9.75309483482066
+29.93136513659281  8205.03364858298846 203.70048022164997 133.67982437482837 313.14312260953017 78.29003556932764 11.54748801370110 1044.72099699982095 9.88440362814734
+30.42441236974887  8194.26130160793400 184.59299320224341 130.90060414734285 316.11077596114535 79.54460397670003 11.76211076966458 1072.70040487430151 10.12720546066069
+30.93570782641099  8183.32209682123721 167.44119761873691 126.10763937673261 318.07266977989138 80.82145727143748 11.98594009135113 1101.86526865666474 10.38373038394032
+31.92630883416595  8163.04766297332026 140.26922500517372 113.99478252566954 317.88530652259146 83.14770734522330 12.42104537026428 1158.33974581994107 10.89452443780815
+32.91690984192091  8144.24910893506058 119.05792346639640 100.80810246705657 312.39863449130485 85.14890038732679 12.85234284243618 1214.06150890051595 11.42347850989474
+34.05259196215597  8124.60934757494215 99.97111781834724 86.43806709438758 300.61658062527658 86.87960504925786 13.33030939621634 1276.10318731695952 12.05178512460471
+35.29523465270979  8105.36080427264824 83.69071548778604 72.65846842329395 282.90656771083434 87.95159039706978 13.81512728730399 1340.85207374139804 12.76465267965806
+36.66527228700479  8086.64261846955378 69.69706716310475 60.06281915658900 259.97122123657249 88.07832712118997 14.27969956798211 1407.68999419840611 13.57825308659419
+38.20375523583810  8068.45816165088490 57.44138015005768 48.80677283699718 232.57612301190954 86.95586648364471 14.68651811735760 1476.55470621060886 14.52047153853275
+39.96139207227289  8050.89681498883419 46.56647995949557 38.90245168439026 201.70099549638979 84.24995875981710 14.97567696666532 1547.08447951536982 15.62314262903155
+42.02550896030483  8033.97169999245671 36.76078768796594 30.20133202570899 168.22978683986778 79.57096172332076 15.05859374455560 1619.27193168032113 16.93490630579746
+44.32412090137599  8018.96672096901966 28.47471409726841 23.08423003818334 135.82892549275056 73.09572632375743 14.83362586895918 1687.32693365609134 18.38912355396305
+46.81974882514553  8006.31077052649653 21.69394281539434 17.42822627377479 106.70406823908149 65.29956431717308 14.25862505451351 1748.37907826304240 19.92572451051782
+49.51135437715445  7995.94087278360803 16.22872156383918 12.96730471030573 81.69575580275861 56.71311101101572 13.33966860961966 1801.61700443673908 21.49756108210813
+50.00000000000000  7994.35698178958410 15.39841472610330 12.29605286381966 77.78664353600600 55.17888164591182 13.14756251713682 1810.06444973708744 21.77101318434512
diff --git a/cpp/tests/data/ode-sir-js-compare.csv → cpp/tests/data/ode-sir-compare.csv b/cpp/tests/data/ode-sir-js-compare.csv → cpp/tests/data/ode-sir-compare.csv
diff --git a/cpp/tests/data/secihurd-compare.csv b/cpp/tests/data/secihurd-compare.csv
@@ -50,4 +50,4 @@ Time             S                E                C                C+
 47.10000000000000 8003.92894509990219 21.08747690227999 16.92856585059912 0.00000000000000 103.98945790666005 0.00000000000000 64.47480658783135 17.02081150115228 1749.05198133817930 23.51795481340114
 48.10000000000000 7999.84084473854455 18.92921094525236 15.16120071824191 0.00000000000000 94.20859326541573 0.00000000000000 61.27126538818517 16.63555801635884 1769.72415687649345 24.22917005151288
 49.10000000000000 7996.16540124720086 16.99674972514239 13.58911131653709 0.00000000000000 85.28098792803146 0.00000000000000 58.08203838912079 16.20670178089632 1788.75584706136237 24.92316255171389
-50.00000000000000 7993.17638252225152 15.42943711043520 12.32062487449738 0.00000000000000 77.92692837057375 0.00000000000000 55.24515753967167 15.78886448206867 1804.58098849539601 25.53161660511103
+50.00000000000000 7993.17638252225152 15.42943711043520 12.32062487449738 0.00000000000000 77.92692837057375 0.00000000000000 55.24515753967167 15.78886448206867 1804.58098849539601 25.53161660511103