Quickstart Guide: Runge Kutta and model products
Source:vignettes/quickstart_products.Rmd
quickstart_products.Rmd
epi_dir = system.file("starter_models", "product_example", "Epi_model", package = "macpan2")
age_dir = system.file("starter_models", "product_example", "Age_model", package = "macpan2")Model definition directories can be converted directly into simulators using the “SimulatorConstructor” function. When using this function users must input the usual arguments for creating a simulator (e.g. state, flow, parameters) as well as “per_capita_transmission” a matrix with a number of columns equal to the number of infectious compartments and a number of rows equal to the number of infectious flows. The “per_capita_transmission” matrix can either be filled directly with numerical values or it can be created as an “empty_matrix” and computed via expressions in the “derivations.R” file.
The SimulatorConstructor function can create simulators that use one of two different integration methods, the specific method is set using the “integration_method” argument. The defult option is “RK4” which denotes the Runge-Kutta 4 algorithm, to use Eulers method set “integration_method = Euler”.
epi_rk4 = SimulatorConstructor(epi_dir,
time_steps = 25L,
state = c(S = 999, E = 1, I = 0, R = 0),
flow = c(total_foi = NA, progression = 0.1, recovery = 0.05),
N = empty_matrix,
transmissability = 0.75,
per_capita_transmission = empty_matrix,
.mats_to_return = c("state"))
epi_rk4$report()
#> matrix time row col value
#> 1 state 0 S 9.990000e+02
#> 2 state 0 E 1.000000e+00
#> 3 state 0 I 0.000000e+00
#> 4 state 0 R 0.000000e+00
#> 5 state 1 S 9.989848e+02
#> 6 state 1 E 9.502965e-01
#> 7 state 1 I 6.390140e-02
#> 8 state 1 R 1.015230e-03
#> 9 state 2 S 9.989388e+02
#> 10 state 2 E 9.336567e-01
#> 11 state 2 I 1.234922e-01
#> 12 state 2 R 4.086569e-03
#> 13 state 3 S 9.988636e+02
#> 14 state 3 E 9.463694e-01
#> 15 state 3 I 1.809617e-01
#> 16 state 3 R 9.105480e-03
#> 17 state 4 S 9.987598e+02
#> 18 state 4 E 9.859539e-01
#> 19 state 4 I 2.382230e-01
#> 20 state 4 R 1.603183e-02
#> 21 state 5 S 9.986272e+02
#> 22 state 5 E 1.050966e+00
#> 23 state 5 I 2.969966e-01
#> 24 state 5 R 2.488603e-02
#> 25 state 6 S 9.984645e+02
#> 26 state 6 E 1.140854e+00
#> 27 state 6 I 3.588803e-01
#> 28 state 6 R 3.574375e-02
#> 29 state 7 S 9.982700e+02
#> 30 state 7 E 1.255856e+00
#> 31 state 7 I 4.254080e-01
#> 32 state 7 R 4.873277e-02
#> 33 state 8 S 9.980409e+02
#> 34 state 8 E 1.396925e+00
#> 35 state 8 I 4.981013e-01
#> 36 state 8 R 6.403167e-02
#> 37 state 9 S 9.977739e+02
#> 38 state 9 E 1.565692e+00
#> 39 state 9 I 5.785159e-01
#> 40 state 9 R 8.187016e-02
#> 41 state 10 S 9.974647e+02
#> 42 state 10 E 1.764440e+00
#> 43 state 10 I 6.682835e-01
#> 44 state 10 R 1.025309e-01
#> 45 state 11 S 9.971084e+02
#> 46 state 11 E 1.996109e+00
#> 47 state 11 I 7.691523e-01
#> 48 state 11 R 1.263526e-01
#> 49 state 12 S 9.966989e+02
#> 50 state 12 E 2.264316e+00
#> 51 state 12 I 8.830273e-01
#> 52 state 12 R 1.537347e-01
#> 53 state 13 S 9.962295e+02
#> 54 state 13 E 2.573390e+00
#> 55 state 13 I 1.012010e+00
#> 56 state 13 R 1.851430e-01
#> 57 state 14 S 9.956920e+02
#> 58 state 14 E 2.928424e+00
#> 59 state 14 I 1.158440e+00
#> 60 state 14 R 2.211169e-01
#> 61 state 15 S 9.950774e+02
#> 62 state 15 E 3.335345e+00
#> 63 state 15 I 1.324941e+00
#> 64 state 15 R 2.622784e-01
#> 65 state 16 S 9.943752e+02
#> 66 state 16 E 3.800997e+00
#> 67 state 16 I 1.514466e+00
#> 68 state 16 R 3.093415e-01
#> 69 state 17 S 9.935733e+02
#> 70 state 17 E 4.333237e+00
#> 71 state 17 I 1.730352e+00
#> 72 state 17 R 3.631250e-01
#> 73 state 18 S 9.926580e+02
#> 74 state 18 E 4.941057e+00
#> 75 state 18 I 1.976377e+00
#> 76 state 18 R 4.245654e-01
#> 77 state 19 S 9.916137e+02
#> 78 state 19 E 5.634709e+00
#> 79 state 19 I 2.256823e+00
#> 80 state 19 R 4.947327e-01
#> 81 state 20 S 9.904227e+02
#> 82 state 20 E 6.425857e+00
#> 83 state 20 I 2.576547e+00
#> 84 state 20 R 5.748487e-01
#> 85 state 21 S 9.890649e+02
#> 86 state 21 E 7.327745e+00
#> 87 state 21 I 2.941060e+00
#> 88 state 21 R 6.663070e-01
#> 89 state 22 S 9.875173e+02
#> 90 state 22 E 8.355375e+00
#> 91 state 22 I 3.356616e+00
#> 92 state 22 R 7.706962e-01
#> 93 state 23 S 9.857541e+02
#> 94 state 23 E 9.525717e+00
#> 95 state 23 I 3.830308e+00
#> 96 state 23 R 8.898264e-01
#> 97 state 24 S 9.837461e+02
#> 98 state 24 E 1.085792e+01
#> 99 state 24 I 4.370175e+00
#> 100 state 24 R 1.025758e+00
#> 101 state 25 S 9.814603e+02
#> 102 state 25 E 1.237357e+01
#> 103 state 25 I 4.985324e+00
#> 104 state 25 R 1.180837e+00
#> 105 state 26 S 9.814603e+02
#> 106 state 26 E 1.237357e+01
#> 107 state 26 I 4.985324e+00
#> 108 state 26 R 1.180837e+00
age_rk4 = SimulatorConstructor(age_dir,
time_steps = 25L,
state = c(young = 333, medium = 333, old = 333),
flow = c(ageing_rate = 0.03, birth_rate = 5, death_rate = 0.01),
per_capita_transmission = empty_matrix,
.mats_to_return = c("state"))
age_rk4$report()
#> matrix time row col value
#> 1 state 0 young 333.0000
#> 2 state 0 medium 333.0000
#> 3 state 0 old 333.0000
#> 4 state 1 young 329.7043
#> 5 state 1 medium 332.9692
#> 6 state 1 old 337.4260
#> 7 state 2 young 326.4740
#> 8 state 2 medium 332.8743
#> 9 state 2 old 341.8213
#> 10 state 3 young 323.3076
#> 11 state 3 medium 332.7179
#> 12 state 3 old 346.1850
#> 13 state 4 young 320.2040
#> 14 state 4 medium 332.5024
#> 15 state 4 old 350.5159
#> 16 state 5 young 317.1619
#> 17 state 5 medium 332.2303
#> 18 state 5 old 354.8133
#> 19 state 6 young 314.1800
#> 20 state 6 medium 331.9038
#> 21 state 6 old 359.0762
#> 22 state 7 young 311.2572
#> 23 state 7 medium 331.5253
#> 24 state 7 old 363.3037
#> 25 state 8 young 308.3924
#> 26 state 8 medium 331.0970
#> 27 state 8 old 367.4951
#> 28 state 9 young 305.5843
#> 29 state 9 medium 330.6208
#> 30 state 9 old 371.6497
#> 31 state 10 young 302.8318
#> 32 state 10 medium 330.0990
#> 33 state 10 old 375.7668
#> 34 state 11 young 300.1339
#> 35 state 11 medium 329.5335
#> 36 state 11 old 379.8456
#> 37 state 12 young 297.4894
#> 38 state 12 medium 328.9263
#> 39 state 12 old 383.8857
#> 40 state 13 young 294.8973
#> 41 state 13 medium 328.2792
#> 42 state 13 old 387.8865
#> 43 state 14 young 292.3566
#> 44 state 14 medium 327.5940
#> 45 state 14 old 391.8474
#> 46 state 15 young 289.8663
#> 47 state 15 medium 326.8725
#> 48 state 15 old 395.7680
#> 49 state 16 young 287.4252
#> 50 state 16 medium 326.1165
#> 51 state 16 old 399.6478
#> 52 state 17 young 285.0326
#> 53 state 17 medium 325.3275
#> 54 state 17 old 403.4864
#> 55 state 18 young 282.6873
#> 56 state 18 medium 324.5072
#> 57 state 18 old 407.2836
#> 58 state 19 young 280.3885
#> 59 state 19 medium 323.6571
#> 60 state 19 old 411.0388
#> 61 state 20 young 278.1353
#> 62 state 20 medium 322.7787
#> 63 state 20 old 414.7519
#> 64 state 21 young 275.9267
#> 65 state 21 medium 321.8735
#> 66 state 21 old 418.4225
#> 67 state 22 young 273.7618
#> 68 state 22 medium 320.9429
#> 69 state 22 old 422.0504
#> 70 state 23 young 271.6399
#> 71 state 23 medium 319.9882
#> 72 state 23 old 425.6355
#> 73 state 24 young 269.5600
#> 74 state 24 medium 319.0108
#> 75 state 24 old 429.1775
#> 76 state 25 young 267.5213
#> 77 state 25 medium 318.0120
#> 78 state 25 old 432.6762
#> 79 state 26 young 267.5213
#> 80 state 26 medium 318.0120
#> 81 state 26 old 432.6762
epi_euler = SimulatorConstructor(epi_dir,
integration_method = Euler,
time_steps = 25L,
state = c(S = 999, E = 1, I = 0, R = 0),
flow = c(total_foi = NA, progression = 0.1, recovery = 0.05),
N = empty_matrix,
transmissability = 0.75,
per_capita_transmission = empty_matrix,
.mats_to_return = c("state"))
epi_euler$report()
#> matrix time row col value
#> 1 state 0 S 999.00000000
#> 2 state 0 E 1.00000000
#> 3 state 0 I 0.00000000
#> 4 state 0 R 0.00000000
#> 5 state 1 S 999.00000000
#> 6 state 1 E 0.90000000
#> 7 state 1 I 0.10000000
#> 8 state 1 R 0.00000000
#> 9 state 2 S 998.92507500
#> 10 state 2 E 0.88492500
#> 11 state 2 I 0.18500000
#> 12 state 2 R 0.00500000
#> 13 state 3 S 998.78647415
#> 14 state 3 E 0.93503335
#> 15 state 3 I 0.26424250
#> 16 state 3 R 0.01425000
#> 17 state 4 S 998.58853277
#> 18 state 4 E 1.03947139
#> 19 state 4 I 0.34453371
#> 20 state 4 R 0.02746213
#> 21 state 5 S 998.33049721
#> 22 state 5 E 1.19355981
#> 23 state 5 I 0.43125416
#> 24 state 5 R 0.04468881
#> 25 state 6 S 998.00759657
#> 26 state 6 E 1.39710447
#> 27 state 6 I 0.52904744
#> 28 state 6 R 0.06625152
#> 29 state 7 S 997.61160155
#> 30 state 7 E 1.65338905
#> 31 state 7 I 0.64230551
#> 32 state 7 R 0.09270389
#> 33 state 8 S 997.13102298
#> 34 state 8 E 1.96862871
#> 35 state 8 I 0.77552914
#> 36 state 8 R 0.12481917
#> 37 state 9 S 996.55104485
#> 38 state 9 E 2.35174397
#> 39 state 9 I 0.93361556
#> 40 state 9 R 0.16359562
#> 41 state 10 S 995.85324818
#> 42 state 10 E 2.81436624
#> 43 state 10 I 1.12210918
#> 44 state 10 R 0.21027640
#> 45 state 11 S 995.01515613
#> 46 state 11 E 3.37102167
#> 47 state 11 I 1.34744034
#> 48 state 11 R 0.26638186
#> 49 state 12 S 994.00961346
#> 50 state 12 E 4.03946217
#> 51 state 12 I 1.61717049
#> 52 state 12 R 0.33375388
#> 53 state 13 S 992.80400120
#> 54 state 13 E 4.84112821
#> 55 state 13 I 1.94025818
#> 56 state 13 R 0.41461240
#> 57 state 14 S 991.35927914
#> 58 state 14 E 5.80173746
#> 59 state 14 I 2.32735809
#> 60 state 14 R 0.51162531
#> 61 state 15 S 989.62884311
#> 62 state 15 E 6.95199974
#> 63 state 15 I 2.79116394
#> 64 state 15 R 0.62799322
#> 65 state 16 S 987.55718085
#> 66 state 16 E 8.32846202
#> 67 state 16 I 3.34680571
#> 68 state 16 R 0.76755141
#> 69 state 17 S 985.07830934
#> 70 state 17 E 9.97448733
#> 71 state 17 I 4.01231163
#> 72 state 17 R 0.93489170
#> 73 state 18 S 982.11397847
#> 74 state 18 E 11.94136947
#> 75 state 18 I 4.80914478
#> 76 state 18 R 1.13550728
#> 77 state 19 S 978.57163224
#> 78 state 19 E 14.28957875
#> 79 state 19 I 5.76282449
#> 80 state 19 R 1.37596452
#> 81 state 20 S 974.34212981
#> 82 state 20 E 17.09012330
#> 83 state 20 I 6.90364114
#> 84 state 20 R 1.66410574
#> 85 state 21 S 969.29724851
#> 86 state 21 E 20.42599228
#> 87 state 21 I 8.26747141
#> 88 state 21 R 2.00928780
#> 89 state 22 S 963.28702054
#> 90 state 22 E 24.39362102
#> 91 state 22 I 9.89669707
#> 92 state 22 R 2.42266137
#> 93 state 23 S 956.13700066
#> 94 state 23 E 29.10427880
#> 95 state 23 I 11.84122432
#> 96 state 23 R 2.91749622
#> 97 state 24 S 947.64562613
#> 98 state 24 E 34.68522545
#> 99 state 24 I 14.15959098
#> 100 state 24 R 3.50955744
#> 101 state 25 S 937.58192028
#> 102 state 25 E 41.28040875
#> 103 state 25 I 16.92013398
#> 104 state 25 R 4.21753699
#> 105 state 26 S 937.58192028
#> 106 state 26 E 41.28040875
#> 107 state 26 I 16.92013398
#> 108 state 26 R 4.21753699
age_euler = SimulatorConstructor(age_dir,
integration_method = Euler,
time_steps = 25L,
state = c(young = 333, medium = 333, old = 333),
flow = c(ageing_rate = 0.03, birth_rate = 5, death_rate = 0.01),
per_capita_transmission = empty_matrix,
.mats_to_return = c("state"))
age_euler$report()
#> matrix time row col value
#> 1 state 0 young 333.0000
#> 2 state 0 medium 333.0000
#> 3 state 0 old 333.0000
#> 4 state 1 young 328.0100
#> 5 state 1 medium 333.0000
#> 6 state 1 old 339.6600
#> 7 state 2 young 323.1697
#> 8 state 2 medium 332.8503
#> 9 state 2 old 346.2534
#> 10 state 3 young 318.4746
#> 11 state 3 medium 332.5599
#> 12 state 3 old 352.7764
#> 13 state 4 young 313.9204
#> 14 state 4 medium 332.1373
#> 15 state 4 old 359.2254
#> 16 state 5 young 309.5028
#> 17 state 5 medium 331.5908
#> 18 state 5 old 365.5973
#> 19 state 6 young 305.2177
#> 20 state 6 medium 330.9282
#> 21 state 6 old 371.8890
#> 22 state 7 young 301.0611
#> 23 state 7 medium 330.1569
#> 24 state 7 old 378.0980
#> 25 state 8 young 297.0293
#> 26 state 8 medium 329.2840
#> 27 state 8 old 384.2217
#> 28 state 9 young 293.1184
#> 29 state 9 medium 328.3163
#> 30 state 9 old 390.2580
#> 31 state 10 young 289.3249
#> 32 state 10 medium 327.2604
#> 33 state 10 old 396.2049
#> 34 state 11 young 285.6451
#> 35 state 11 medium 326.1223
#> 36 state 11 old 402.0607
#> 37 state 12 young 282.0758
#> 38 state 12 medium 324.9080
#> 39 state 12 old 407.8237
#> 40 state 13 young 278.6135
#> 41 state 13 medium 323.6231
#> 42 state 13 old 413.4927
#> 43 state 14 young 275.2551
#> 44 state 14 medium 322.2728
#> 45 state 14 old 419.0665
#> 46 state 15 young 271.9974
#> 47 state 15 medium 320.8622
#> 48 state 15 old 424.5440
#> 49 state 16 young 268.8375
#> 50 state 16 medium 319.3963
#> 51 state 16 old 429.9245
#> 52 state 17 young 265.7724
#> 53 state 17 medium 317.8795
#> 54 state 17 old 435.2071
#> 55 state 18 young 262.7992
#> 56 state 18 medium 316.3163
#> 57 state 18 old 440.3914
#> 58 state 19 young 259.9152
#> 59 state 19 medium 314.7108
#> 60 state 19 old 445.4770
#> 61 state 20 young 257.1178
#> 62 state 20 medium 313.0669
#> 63 state 20 old 450.4635
#> 64 state 21 young 254.4043
#> 65 state 21 medium 311.3885
#> 66 state 21 old 455.3509
#> 67 state 22 young 251.7721
#> 68 state 22 medium 309.6789
#> 69 state 22 old 460.1391
#> 70 state 23 young 249.2190
#> 71 state 23 medium 307.9417
#> 72 state 23 old 464.8280
#> 73 state 24 young 246.7424
#> 74 state 24 medium 306.1801
#> 75 state 24 old 469.4180
#> 76 state 25 young 244.3401
#> 77 state 25 medium 304.3969
#> 78 state 25 old 473.9092
#> 79 state 26 young 244.3401
#> 80 state 26 medium 304.3969
#> 81 state 26 old 473.9092Eventually it will be possible to create a new model directory that corresponds to the product of two models using the “ModelProduct” function. This function accepts as arguments two model simulators and will write a model definition corresponding to the product of those models. Currently the file writing capability is unavailable so this function is not complete. Note that every model directory created this way will have a blank “derivations.json” file, the use must either give numerical values for the “per_capita_transmission” matrix directly or add expressions to “derivations.json” that compute the “per_capita_transmission” matrix from some other supplied parameters.