From e2cb0d4668f17f57c65f3ff94a7e17c784eaf4ba Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Mon, 9 Nov 2020 07:31:00 +0100 Subject: Custom analytical solutions for saemix Currently SFO-SFO and DFOP-SFO. Speed increase factor about 60 --- docs/dev/pkgdown.yml | 2 +- docs/dev/reference/Rplot001.png | Bin 27839 -> 1011 bytes docs/dev/reference/saem.html | 80 +++++++++++++++++++++++++++++----------- 3 files changed, 59 insertions(+), 23 deletions(-) (limited to 'docs/dev') diff --git a/docs/dev/pkgdown.yml b/docs/dev/pkgdown.yml index d9720761..8d117fb1 100644 --- a/docs/dev/pkgdown.yml +++ b/docs/dev/pkgdown.yml @@ -10,7 +10,7 @@ articles: web_only/NAFTA_examples: NAFTA_examples.html web_only/benchmarks: benchmarks.html web_only/compiled_models: compiled_models.html -last_built: 2020-11-08T02:20Z +last_built: 2020-11-09T06:03Z urls: reference: https://pkgdown.jrwb.de/mkin/reference article: https://pkgdown.jrwb.de/mkin/articles diff --git a/docs/dev/reference/Rplot001.png b/docs/dev/reference/Rplot001.png index cfc5bc2b..17a35806 100644 Binary files a/docs/dev/reference/Rplot001.png and b/docs/dev/reference/Rplot001.png differ diff --git a/docs/dev/reference/saem.html b/docs/dev/reference/saem.html index 06fcfaa7..25608fc8 100644 --- a/docs/dev/reference/saem.html +++ b/docs/dev/reference/saem.html @@ -229,28 +229,28 @@ using mmkin.

state.ini = c(parent = 100), fixed_initials = "parent", quiet = TRUE) f_saem_p0_fixed <- saem(f_mmkin_parent_p0_fixed)
#> Running main SAEM algorithm -#> [1] "Sun Nov 8 02:44:42 2020" +#> [1] "Mon Nov 9 07:04:09 2020" #> .... #> Minimisation finished -#> [1] "Sun Nov 8 02:44:43 2020"
+#> [1] "Mon Nov 9 07:04:11 2020"
f_mmkin_parent <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE) f_saem_sfo <- saem(f_mmkin_parent["SFO", ])
#> Running main SAEM algorithm -#> [1] "Sun Nov 8 02:44:45 2020" +#> [1] "Mon Nov 9 07:04:12 2020" #> .... #> Minimisation finished -#> [1] "Sun Nov 8 02:44:46 2020"
f_saem_fomc <- saem(f_mmkin_parent["FOMC", ]) +#> [1] "Mon Nov 9 07:04:13 2020"
f_saem_fomc <- saem(f_mmkin_parent["FOMC", ])
#> Running main SAEM algorithm -#> [1] "Sun Nov 8 02:44:47 2020" +#> [1] "Mon Nov 9 07:04:14 2020" #> .... #> Minimisation finished -#> [1] "Sun Nov 8 02:44:49 2020"
f_saem_dfop <- saem(f_mmkin_parent["DFOP", ]) +#> [1] "Mon Nov 9 07:04:16 2020"
f_saem_dfop <- saem(f_mmkin_parent["DFOP", ])
#> Running main SAEM algorithm -#> [1] "Sun Nov 8 02:44:49 2020" +#> [1] "Mon Nov 9 07:04:16 2020" #> .... #> Minimisation finished -#> [1] "Sun Nov 8 02:44:52 2020"
-# The returned saem.mmkin object contains an SaemixObject, we can use +#> [1] "Mon Nov 9 07:04:19 2020"
+# The returned saem.mmkin object contains an SaemixObject, therefore we can use # functions from saemix library(saemix)
#> Package saemix, version 3.1.9000 @@ -262,33 +262,69 @@ using mmkin.

f_mmkin_parent_tc <- update(f_mmkin_parent, error_model = "tc") f_saem_fomc_tc <- saem(f_mmkin_parent_tc["FOMC", ])
#> Running main SAEM algorithm -#> [1] "Sun Nov 8 02:44:54 2020" +#> [1] "Mon Nov 9 07:04:21 2020" #> .... #> Minimisation finished -#> [1] "Sun Nov 8 02:44:59 2020"
compare.saemix(list(f_saem_fomc$so, f_saem_fomc_tc$so)) +#> [1] "Mon Nov 9 07:04:26 2020"
compare.saemix(list(f_saem_fomc$so, f_saem_fomc_tc$so))
#> Likelihoods computed by importance sampling
#> AIC BIC #> 1 467.7644 465.0305 #> 2 469.4862 466.3617
-dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"), +sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"), A1 = mkinsub("SFO")) -
#> Successfully compiled differential equation model from auto-generated C code.
f_mmkin <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, solution_type = "analytical") -# This takes about 4 minutes on my system -f_saem <- saem(f_mmkin) +
#> Successfully compiled differential equation model from auto-generated C code.
fomc_sfo <- mkinmod(parent = mkinsub("FOMC", "A1"), + A1 = mkinsub("SFO")) +
#> Successfully compiled differential equation model from auto-generated C code.
dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"), + A1 = mkinsub("SFO")) +
#> Successfully compiled differential equation model from auto-generated C code.
# The following fit uses analytical solutions for SFO-SFO and DFOP-SFO, +# and compiled ODEs for FOMC, both are fast +f_mmkin <- mmkin(list( + "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo), + ds, quiet = TRUE) +# These take about five seconds each on this system, as we use +# analytical solutions written for saemix. When using the analytical +# solutions written for mkin this took around four minutes +f_saem_sfo_sfo <- saem(f_mmkin["SFO-SFO", ]) +
#> Running main SAEM algorithm +#> [1] "Mon Nov 9 07:04:28 2020" +#> .... +#> Minimisation finished +#> [1] "Mon Nov 9 07:04:33 2020"
f_saem_dfop_sfo <- saem(f_mmkin["SFO-SFO", ])
#> Running main SAEM algorithm -#> [1] "Sun Nov 8 02:45:00 2020" +#> [1] "Mon Nov 9 07:04:33 2020" #> .... #> Minimisation finished -#> [1] "Sun Nov 8 02:49:01 2020"
#> Error in exp(trans_k): non-numeric argument to mathematical function
#> Timing stopped at: 317.7 169.9 258.9
-f_mmkin_des <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, solution_type = "deSolve") +#> [1] "Mon Nov 9 07:04:39 2020"
# Using a single core, the following takes about 6 minutes, using 10 cores # it is slower instead of faster -f_saem_des <- saem(f_mmkin_des, cores = 1) +f_saem_fomc <- saem(f_mmkin["FOMC-SFO", ], cores = 1)
#> Running main SAEM algorithm -#> [1] "Sun Nov 8 02:49:20 2020" +#> [1] "Mon Nov 9 07:04:39 2020" +#> DLSODA- At current T (=R1), MXSTEP (=I1) steps +#> taken on this call before reaching TOUT +#> In above message, I1 = 5000 +#> +#> In above message, R1 = 0.00156238 +#> +#> DLSODA- At T (=R1) and step size H (=R2), the +#> corrector convergence failed repeatedly +#> or with ABS(H) = HMIN +#> In above message, R1 = 0, R2 = 1.1373e-10 +#> +#> DLSODA- At current T (=R1), MXSTEP (=I1) steps +#> taken on this call before reaching TOUT +#> In above message, I1 = 5000 +#> +#> In above message, R1 = 2.24752e-06 +#> +#> DLSODA- At current T (=R1), MXSTEP (=I1) steps +#> taken on this call before reaching TOUT +#> In above message, I1 = 5000 +#> +#> In above message, R1 = 0.000585935 +#> #> .... #> Minimisation finished -#> [1] "Sun Nov 8 02:57:30 2020"
#> Error in exp(trans_k): non-numeric argument to mathematical function
#> Timing stopped at: 591.7 205.6 521.4
compare.saemix(list(f_saem$so, f_saem_des$so)) -
#> Error in compare.saemix(list(f_saem$so, f_saem_des$so)): object 'f_saem' not found
# } +#> [1] "Mon Nov 9 07:11:24 2020"
#> Warning: Creating predictions from the saemix model failed
# }