1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
|
#' Create degradation functions for known analytical solutions
#'
#' @param spec List of model specifications as contained in mkinmod objects
#' @param use_of_ff Minimum or maximum use of formation fractions
#' @return Degradation function to be attached to mkinmod objects
#' @examples
#'
#' SFO_SFO <- mkinmod(
#' parent = mkinsub("SFO", "m1"),
#' m1 = mkinsub("SFO"))
#' FOCUS_D <- subset(FOCUS_2006_D, value != 0) # to avoid warnings
#' fit_1 <- mkinfit(SFO_SFO, FOCUS_D, solution_type = "analytical", quiet = TRUE)
#' \dontrun{
#' fit_2 <- mkinfit(SFO_SFO, FOCUS_D, solution_type = "deSolve", quiet = TRUE)
#' if (require(rbenchmark))
#' benchmark(
#' analytical = mkinfit(SFO_SFO, FOCUS_D, solution_type = "analytical", quiet = TRUE),
#' deSolve = mkinfit(SFO_SFO, FOCUS_D, solution_type = "deSolve", quiet = TRUE),
#' replications = 2)
#' DFOP_SFO <- mkinmod(
#' parent = mkinsub("DFOP", "m1"),
#' m1 = mkinsub("SFO"))
#' benchmark(
#' analytical = mkinfit(DFOP_SFO, FOCUS_D, solution_type = "analytical", quiet = TRUE),
#' deSolve = mkinfit(DFOP_SFO, FOCUS_D, solution_type = "deSolve", quiet = TRUE),
#' replications = 2)
#' }
#' @export
create_deg_func <- function(spec, use_of_ff = c("min", "max")) {
use_of_ff <- match.arg(use_of_ff)
min_ff <- use_of_ff == "min"
obs_vars <- names(spec)
parent <- obs_vars[1]
parent_type <- spec[[1]]$type
supported <- TRUE # This may be modified below
predicted_text <- character(0)
if (parent_type == "SFO") {
if (min_ff) {
targets <- c(spec[[1]]$to, if (spec[[1]]$sink) "sink" else NULL)
k_parent <- paste(paste0("k_", parent, "_", targets), collapse = " + ")
} else {
k_parent <- paste0("k_", parent)
}
}
predicted_text[parent] <- paste0(parent_type, ".solution(t, odeini['", parent,
if (parent_type == "SFORB") "_free", "'], ",
switch(parent_type,
SFO = k_parent,
FOMC = "alpha, beta",
IORE = paste0("k__iore_", parent, if (min_ff) "_sink" else "", ", N_", parent),
DFOP = "k1, k2, g",
SFORB = paste0("k_", parent, "_free_bound, k_", parent, "_bound_free, k_", parent, "_free", if (min_ff) "_sink" else ""),
HS = "k1, k2, tb",
logistic = "kmax, k0, r"
),
")")
if (length(obs_vars) >= 2) {
supported <- FALSE # except for the following cases
n1 <- obs_vars[1]
n2 <- obs_vars[2]
n10 <- paste0("odeini['", parent, "']")
n20 <- paste0("odeini['", n2, "']")
# sfo_sfo
if (all(spec[[1]]$sink == FALSE, length(obs_vars) == 2,
spec[[1]]$type == "SFO", spec[[2]]$type == "SFO")) {
supported <- TRUE
k1 <- k_parent
k2 <- paste0("k_", n2, if(min_ff) "_sink" else "")
predicted_text[n2] <- paste0(
"(((", k2, "-", k1, ")*", n20, "-", k1, "*", n10, ")*exp(-", k2, "*t)+",
k1, "*", n10, "*exp(-", k1, "*t))/(", k2, "-", k1, ")")
}
# sfo_f12_sfo
if (all(use_of_ff == "max", spec[[1]]$sink == TRUE, length(obs_vars) == 2,
spec[[1]]$type == "SFO", spec[[2]]$type == "SFO")) {
supported <- TRUE
k1 <- k_parent
k2 <- paste0("k_", n2)
f12 <- paste0("f_", n1, "_to_", n2)
predicted_text[n2] <- paste0(
"(((", k2, "-", k1, ")*", n20, "-", f12, "*", k1, "*", n10, ")*exp(-", k2, "*t)+",
f12, "*", k1, "*", n10, "*exp(-", k1, "*t))/(", k2, "-", k1, ")")
}
# sfo_k120_sfo
if (all(use_of_ff == "min", spec[[1]]$sink == TRUE, length(obs_vars) == 2,
spec[[1]]$type == "SFO", spec[[2]]$type == "SFO")) {
supported <- TRUE
k12 <- paste0("k_", n1, "_", n2)
k10 <- paste0("k_", n1, "_sink")
k2 <- paste0("k_", n2, "_sink")
predicted_text[n2] <- paste0(
"(((", k2, "-", k12, "-", k10, ")*", n20, "-", k12, "*", n10, ")*exp(-", k2, "*t)+",
k12, "*", n10, "*exp(-(", k_parent, ")*t))/(", k2, "-(", k_parent, "))")
}
# dfop_f12_sfo
if (all(use_of_ff == "max", spec[[1]]$sink == TRUE, length(obs_vars) == 2,
spec[[1]]$type == "DFOP", spec[[2]]$type == "SFO")) {
supported <- TRUE
f12 <- paste0("f_", n1, "_to_", n2)
k2 <- paste0("k_", n2)
predicted_text[n2] <- paste0(
"((", f12, "* g - ", f12, ") * k2 * ", n10, " * exp(- k2 * t))/(k2 - ", k2, ") - ",
"((", f12, "* g) * k1 * ", n10, " * exp(- k1 * t))/(k1 - ", k2, ") + ",
"(((k1 - ", k2, ") * k2 - ", k2, "* k1 + ", k2, "^2) * ", n20, "+",
"((", f12, "* k1 + (", f12, "*g - ", f12, ") * ", k2, ") * k2 - ", f12, " * g * ", k2, " * k1) * ", n10, ") * ",
"exp( - ", k2, " * t)/((k1 - ", k2, ") * k2 - ", k2, " * k1 + ", k2, "^2)")
}
}
if (supported) {
deg_func <- function(observed, odeini, odeparms) {}
f_body <- paste0("{\n",
"predicted <- numeric(0)\n",
"with(as.list(odeparms), {\n")
for (obs_var in obs_vars) {
f_body <- paste0(f_body,
"t <- observed[observed$name == '", obs_var, "', 'time']\n",
"predicted <<- c(predicted, ", predicted_text[obs_var], ")\n")
}
f_body <- paste0(f_body,
"})\n",
"return(predicted)\n}\n")
body(deg_func) <- parse(text = f_body)
return(deg_func)
} else {
return(NULL)
}
}
|
