| Title: | Dynamic Modeling and Parameter Estimation in ODE Models |
|---|---|
| Description: | The framework provides functions to generate ODEs of reaction networks, parameter transformations, observation functions, residual functions, etc. The framework follows the paradigm that derivative information should be used for optimization whenever possible. Therefore, all major functions produce and can handle expressions for symbolic derivatives. |
| Authors: | Daniel Kaschek |
| Maintainer: | Daniel Kaschek <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 1.1.0 |
| Built: | 2025-02-03 15:24:32 UTC |
| Source: | https://github.com/dkaschek/dmod |
Used to concatenate observation functions, prediction functions and parameter transformation functions.
## S3 method for class 'fn' p1 * p2## S3 method for class 'fn' p1 * p2
p1 |
function of class |
p2 |
function of class |
Object of the same class as x1 and x2.
# Define a time grid on which to make a prediction by peace-wise linear function. # Then define a (generic) prediction function based on thid grid. times <- 0:5 grid <- data.frame(name = "A", time = times, row.names = paste0("p", times)) x <- Xd(grid) # Define an observable and an observation function observables <- eqnvec(Aobs = "s*A") g <- Y(g = observables, f = NULL, states = "A", parameters = "s") # Collect parameters and define an overarching parameter transformation # for two "experimental condtions". dynpars <- attr(x, "parameters") obspars <- attr(g, "parameters") innerpars <- c(dynpars, obspars) trafo <- structure(innerpars, names = innerpars) trafo_C1 <- replaceSymbols(innerpars, paste(innerpars, "C1", sep = "_"), trafo) trafo_C2 <- replaceSymbols(innerpars, paste(innerpars, "C2", sep = "_"), trafo) p <- NULL p <- p + P(trafo = trafo_C1, condition = "C1") p <- p + P(trafo = trafo_C2, condition = "C2") # Collect outer (overarching) parameters and # initialize with random values outerpars <- attr(p, "parameters") pars <- structure(runif(length(outerpars), 0, 1), names = outerpars) # Predict internal/unobserved states out1 <- (x*p)(times, pars) plot(out1) # Predict observed states in addition to unobserved out2 <- (g*x*p)(times, pars) plot(out2)# Define a time grid on which to make a prediction by peace-wise linear function. # Then define a (generic) prediction function based on thid grid. times <- 0:5 grid <- data.frame(name = "A", time = times, row.names = paste0("p", times)) x <- Xd(grid) # Define an observable and an observation function observables <- eqnvec(Aobs = "s*A") g <- Y(g = observables, f = NULL, states = "A", parameters = "s") # Collect parameters and define an overarching parameter transformation # for two "experimental condtions". dynpars <- attr(x, "parameters") obspars <- attr(g, "parameters") innerpars <- c(dynpars, obspars) trafo <- structure(innerpars, names = innerpars) trafo_C1 <- replaceSymbols(innerpars, paste(innerpars, "C1", sep = "_"), trafo) trafo_C2 <- replaceSymbols(innerpars, paste(innerpars, "C2", sep = "_"), trafo) p <- NULL p <- p + P(trafo = trafo_C1, condition = "C1") p <- p + P(trafo = trafo_C2, condition = "C2") # Collect outer (overarching) parameters and # initialize with random values outerpars <- attr(p, "parameters") pars <- structure(runif(length(outerpars), 0, 1), names = outerpars) # Predict internal/unobserved states out1 <- (x*p)(times, pars) plot(out1) # Predict observed states in addition to unobserved out2 <- (g*x*p)(times, pars) plot(out2)
The %.*% operator allows to multiply objects of class objlist or objfn with
a scalar.
x1 %.*% x2x1 %.*% x2
x1 |
object of class objfn or objlist. |
x2 |
numeric of length one. |
An objective function or objlist object.
Used to merge datasets with overlapping conditions.
## S3 method for class 'datalist' data1 + data2## S3 method for class 'datalist' data1 + data2
data1 |
dataset of class |
data2 |
dataset of class |
Each data list contains data frames for a number of conditions. The direct sum of datalist is meant as merging the two data lists and returning the overarching datalist.
Object of class datalist for the
union of conditions.
# Start with two data frames mydata1 <- data.frame( name = "A", time = 0:1, value = 1:2, sigma = .1, compound = c("DEM", "APAP"), dose = "0.1" ) mydata2 <- data.frame( name = "A", time = 0:1, value = 3:4, sigma = .1, compound = c("APAP", "DCF"), dose = "0.1" ) # Create datalists from dataframes data1 <- as.datalist(mydata1, split.by = c("compound", "dose")) data2 <- as.datalist(mydata2, split.by = c("compound", "dose")) # Direct sum of datalists data <- data1 + data2 print(data) # Check the condition.grid (if available) condition.grid <- attr(data, "condition.grid") print(condition.grid)# Start with two data frames mydata1 <- data.frame( name = "A", time = 0:1, value = 1:2, sigma = .1, compound = c("DEM", "APAP"), dose = "0.1" ) mydata2 <- data.frame( name = "A", time = 0:1, value = 3:4, sigma = .1, compound = c("APAP", "DCF"), dose = "0.1" ) # Create datalists from dataframes data1 <- as.datalist(mydata1, split.by = c("compound", "dose")) data2 <- as.datalist(mydata2, split.by = c("compound", "dose")) # Direct sum of datalists data <- data1 + data2 print(data) # Check the condition.grid (if available) condition.grid <- attr(data, "condition.grid") print(condition.grid)
Used to add prediction function, parameter transformation functions or observation functions.
## S3 method for class 'fn' x1 + x2## S3 method for class 'fn' x1 + x2
x1 |
function of class |
x2 |
function of class |
Each prediction function is associated to a number of conditions. Adding functions means merging or overwriting the set of conditions.
Object of the same class as x1 and x2 which returns results for the
union of conditions.
# Define a time grid on which to make a prediction by peace-wise linear function. # Then define a (generic) prediction function based on thid grid. times <- 0:5 grid <- data.frame(name = "A", time = times, row.names = paste0("p", times)) x <- Xd(grid) # Define an observable and an observation function observables <- eqnvec(Aobs = "s*A") g <- Y(g = observables, f = NULL, states = "A", parameters = "s") # Collect parameters and define an overarching parameter transformation # for two "experimental condtions". dynpars <- attr(x, "parameters") obspars <- attr(g, "parameters") innerpars <- c(dynpars, obspars) trafo <- structure(innerpars, names = innerpars) trafo_C1 <- replaceSymbols(innerpars, paste(innerpars, "C1", sep = "_"), trafo) trafo_C2 <- replaceSymbols(innerpars, paste(innerpars, "C2", sep = "_"), trafo) p <- NULL p <- p + P(trafo = trafo_C1, condition = "C1") p <- p + P(trafo = trafo_C2, condition = "C2") # Collect outer (overarching) parameters and # initialize with random values outerpars <- attr(p, "parameters") pars <- structure(runif(length(outerpars), 0, 1), names = outerpars) # Predict internal/unobserved states out1 <- (x*p)(times, pars) plot(out1) # Predict observed states in addition to unobserved out2 <- (g*x*p)(times, pars) plot(out2)# Define a time grid on which to make a prediction by peace-wise linear function. # Then define a (generic) prediction function based on thid grid. times <- 0:5 grid <- data.frame(name = "A", time = times, row.names = paste0("p", times)) x <- Xd(grid) # Define an observable and an observation function observables <- eqnvec(Aobs = "s*A") g <- Y(g = observables, f = NULL, states = "A", parameters = "s") # Collect parameters and define an overarching parameter transformation # for two "experimental condtions". dynpars <- attr(x, "parameters") obspars <- attr(g, "parameters") innerpars <- c(dynpars, obspars) trafo <- structure(innerpars, names = innerpars) trafo_C1 <- replaceSymbols(innerpars, paste(innerpars, "C1", sep = "_"), trafo) trafo_C2 <- replaceSymbols(innerpars, paste(innerpars, "C2", sep = "_"), trafo) p <- NULL p <- p + P(trafo = trafo_C1, condition = "C1") p <- p + P(trafo = trafo_C2, condition = "C2") # Collect outer (overarching) parameters and # initialize with random values outerpars <- attr(p, "parameters") pars <- structure(runif(length(outerpars), 0, 1), names = outerpars) # Predict internal/unobserved states out1 <- (x*p)(times, pars) plot(out1) # Predict observed states in addition to unobserved out2 <- (g*x*p)(times, pars) plot(out2)
Direct sum of objective functions
## S3 method for class 'objfn' x1 + x2## S3 method for class 'objfn' x1 + x2
x1 |
function of class |
x2 |
function of class |
The objective functions are evaluated and their results as added. Sometimes, the evaluation of an objective function depends on results that have been computed internally in a preceding objective function. Therefore, environments are forwarded and all evaluations take place in the same environment. The first objective function in a sum of functions generates a new environment.
Object of class objfn.
normL2, constraintL2, priorL2, datapointL2
## Generate three objective functions prior <- structure(rep(0, 5), names = letters[1:5]) obj1 <- constraintL2(mu = prior, attr.name = "center") obj2 <- constraintL2(mu = prior + 1, attr.name = "right") obj3 <- constraintL2(mu = prior - 1, attr.name = "left") ## Evaluate first objective function on a random vector pouter <- prior + rnorm(length(prior)) print(obj1(pouter)) ## Split into fixed and non-fixed part fixed <- pouter[4:5] pouter <- pouter[1:3] print(obj1(pouter, fixed = fixed)) ## Visualize the result by a parameter profile myfit <- trust(obj1, pouter, rinit = 1, rmax = 10, fixed = fixed) myprof <- profile(obj1, myfit$argument, "a", fixed = fixed) plotProfile(myprof) ## Create new objective function by adding the single ones, ## then evalue the random vector again pouter <- prior + rnorm(length(prior)) obj <- obj1 + obj2 + obj3 print(obj(pouter))## Generate three objective functions prior <- structure(rep(0, 5), names = letters[1:5]) obj1 <- constraintL2(mu = prior, attr.name = "center") obj2 <- constraintL2(mu = prior + 1, attr.name = "right") obj3 <- constraintL2(mu = prior - 1, attr.name = "left") ## Evaluate first objective function on a random vector pouter <- prior + rnorm(length(prior)) print(obj1(pouter)) ## Split into fixed and non-fixed part fixed <- pouter[4:5] pouter <- pouter[1:3] print(obj1(pouter, fixed = fixed)) ## Visualize the result by a parameter profile myfit <- trust(obj1, pouter, rinit = 1, rmax = 10, fixed = fixed) myprof <- profile(obj1, myfit$argument, "a", fixed = fixed) plotProfile(myprof) ## Create new objective function by adding the single ones, ## then evalue the random vector again pouter <- prior + rnorm(length(prior)) obj <- obj1 + obj2 + obj3 print(obj(pouter))
Add two lists element by element
## S3 method for class 'objlist' out1 + out2## S3 method for class 'objlist' out1 + out2
out1 |
List of numerics or matrices |
out2 |
List with the same structure as out1 (there will be no warning when mismatching) |
If out1 has names, out2 is assumed to share these names. Each element of the list out1
is inspected. If it has a names attributed, it is used to do a matching between out1 and out2.
The same holds for the attributed dimnames. In all other cases, the "+" operator is applied
the corresponding elements of out1 and out2 as they are.
List of length of out1.
objlist1 <- dMod:::init_empty_objlist(c(a = 1, b = 2)) objlist1$value = 1; objlist1$gradient[1:2] <- 1; objlist1$hessian[1:4] <- 1 objlist2 <- dMod:::init_empty_objlist(c(a = 1, d = 2)) objlist2$value = 1; objlist2$gradient[1:2] <- 1; objlist2$hessian[1:4] <- 1 objlist1 + objlist2objlist1 <- dMod:::init_empty_objlist(c(a = 1, b = 2)) objlist1$value = 1; objlist1$gradient[1:2] <- 1; objlist1$hessian[1:4] <- 1 objlist2 <- dMod:::init_empty_objlist(c(a = 1, d = 2)) objlist2$value = 1; objlist2$gradient[1:2] <- 1; objlist2$hessian[1:4] <- 1 objlist1 + objlist2
Add reaction to reaction table
addReaction(eqnlist, from, to, rate, description = names(rate))addReaction(eqnlist, from, to, rate, description = names(rate))
eqnlist |
equation list, see eqnlist |
from |
character with the left hand side of the reaction, e.g. "2*A + B" |
to |
character with the right hand side of the reaction, e.g. "C + 2*D" |
rate |
character. The rate associated with the reaction. The name is employed as a description of the reaction. |
description |
Optional description instead of |
An object of class eqnlist.
f <- eqnlist() f <- addReaction(f, "2*A+B", "C + 2*D", "k1*B*A^2") f <- addReaction(f, "C + A", "B + A", "k2*C*A") # Write your example here. You can also add more Start..End blocks if needed. # Please mask all output such as print() with the special tag # # such that the test is not littered. Statements guarded by are enabled # in the example file which is extracted from this test file. To extract the # example run # extractExamples() # on the R command line. ## Generate another equation list eq <- eqnlist() eq <- addReaction(eq, "A", "pA", "act_A * A * stimulus", "Phosphorylation of A") eq <- addReaction(eq, "pA", "A", "deact_A * pA", "Deposphorylation of pA") eq <- addReaction(eq, "2*pA", "pA_pA", "form_complex_pA * pA^2", "Complex formation of pA") eq <- addReaction(eq, "B", "pB", "act_B * B * pA_pA", "Phosphorylation of B") eq <- addReaction(eq, "pB", "B", "deact_B * pB", "Deposphorylation of pB") ## Extract data.frame of reactions reactions <- getReactions(eq) print(reactions) ## Get conserved quantities cq <- conservedQuantities(eq$smatrix) print(cq) ## Get fluxes fluxes <- getFluxes(eq) print(fluxes) ## Subsetting of equation list subeq1 <- subset(eq, "pB" %in% Product) print(subeq1) subeq2 <- subset(eq, grepl("not_available", Description)) print(subeq2) ## Time derivatives of observables observables <- eqnvec(pA_obs = "s1*pA", tA_obs = "s2*(A + pA)") dobs <- dot(observables, eq) ## Combined equation vector for ODE and observables f <- c(as.eqnvec(eq), dobs) print(f)f <- eqnlist() f <- addReaction(f, "2*A+B", "C + 2*D", "k1*B*A^2") f <- addReaction(f, "C + A", "B + A", "k2*C*A") # Write your example here. You can also add more Start..End blocks if needed. # Please mask all output such as print() with the special tag # # such that the test is not littered. Statements guarded by are enabled # in the example file which is extracted from this test file. To extract the # example run # extractExamples() # on the R command line. ## Generate another equation list eq <- eqnlist() eq <- addReaction(eq, "A", "pA", "act_A * A * stimulus", "Phosphorylation of A") eq <- addReaction(eq, "pA", "A", "deact_A * pA", "Deposphorylation of pA") eq <- addReaction(eq, "2*pA", "pA_pA", "form_complex_pA * pA^2", "Complex formation of pA") eq <- addReaction(eq, "B", "pB", "act_B * B * pA_pA", "Phosphorylation of B") eq <- addReaction(eq, "pB", "B", "deact_B * pB", "Deposphorylation of pB") ## Extract data.frame of reactions reactions <- getReactions(eq) print(reactions) ## Get conserved quantities cq <- conservedQuantities(eq$smatrix) print(cq) ## Get fluxes fluxes <- getFluxes(eq) print(fluxes) ## Subsetting of equation list subeq1 <- subset(eq, "pB" %in% Product) print(subeq1) subeq2 <- subset(eq, grepl("not_available", Description)) print(subeq2) ## Time derivatives of observables observables <- eqnvec(pA_obs = "s1*pA", tA_obs = "s2*(A + pA)") dobs <- dot(observables, eq) ## Combined equation vector for ODE and observables f <- c(as.eqnvec(eq), dobs) print(f)
Append an objective function to a basic dMod.frame
appendObj( dMod.frame, prd = list(g * (x * p)), obj_data = list(normL2(data, prd, e)), obj = list(obj_data), pars = list(structure(rnorm(length(getParameters(obj))), names = getParameters(obj))), times = list(seq(min(as.data.frame(data)[["time"]]), max(as.data.frame(data)[["time"]]) * 1.1, length.out = 200)), ..., keepCalls = F )appendObj( dMod.frame, prd = list(g * (x * p)), obj_data = list(normL2(data, prd, e)), obj = list(obj_data), pars = list(structure(rnorm(length(getParameters(obj))), names = getParameters(obj))), times = list(seq(min(as.data.frame(data)[["time"]]), max(as.data.frame(data)[["time"]]) * 1.1, length.out = 200)), ..., keepCalls = F )
dMod.frame |
A dMod.frame |
prd |
Expression after which the concatenated prediction function is formed. Has to wrapped in list() |
obj_data |
Expression after which the objective function which describes the data is formed. Has to wrapped in list() |
obj |
This object is taken by the standard fitting functions. At typical expression would be |
pars |
A named vector of parameters to run e.g. test simulations of the model. Defaults to random parameters |
times |
A vector of times to run e.g. test simulations of the model. Defaults to |
... |
Other columns which are mutations of existing ones or new columns. |
keepCalls |
Store a record of the calls in a new colun? See mutatedMod.frame. |
The dMod.frame augmented by standardized columns
Most plotting functions rely on a column "parframes" to be existent in the dMod.frame
appendParframes( dMod.frame, parframes = list(as.parframe(fits)), keepFits = F, ..., keepCalls = F )appendParframes( dMod.frame, parframes = list(as.parframe(fits)), keepFits = F, ..., keepCalls = F )
dMod.frame |
A dmod.frame, preferably with a column |
parframes |
Expression to turn a column containing a parlist (e.g. fits) into a column of parframes |
keepFits |
|
... |
Other columns you want to mutate |
keepCalls |
Store a record of the calls in a new colun? See mutatedMod.frame. |
The dMod.frame containing the column "parframes"
Coerce to a Data Frame
## S3 method for class 'datalist' as.data.frame(x, ...) ## S3 method for class 'prdlist' as.data.frame(x, ..., data = NULL, errfn = NULL)## S3 method for class 'datalist' as.data.frame(x, ...) ## S3 method for class 'prdlist' as.data.frame(x, ..., data = NULL, errfn = NULL)
x |
any R object |
... |
not used right now |
data |
data list oject |
errfn |
obsfn object, the error model function to predict sigma |
a data frame
Coerce equation list into a data frame
## S3 method for class 'eqnlist' as.data.frame(x, ...)## S3 method for class 'eqnlist' as.data.frame(x, ...)
x |
object of class eqnlist |
... |
other arguments |
a data.frame with columns "Description" (character),
"Rate" (character), and one column per ODE state with the state names.
The state columns correspond to the stoichiometric matrix.
An equation list stores an ODE in a list format. The function translates this list into the right-hand sides of the ODE.
as.eqnvec(x, ...) ## S3 method for class 'character' as.eqnvec(x = NULL, names = NULL, ...) ## S3 method for class 'eqnlist' as.eqnvec(x, ...)as.eqnvec(x, ...) ## S3 method for class 'character' as.eqnvec(x = NULL, names = NULL, ...) ## S3 method for class 'eqnlist' as.eqnvec(x, ...)
x |
object of class |
... |
arguments going to the corresponding methods |
names |
character, the left-hand sides of the equation |
If x is of class eqnlist, getFluxes is called and coerced
into a vector of equations.
object of class eqnvec.
Coerce to eventlist
as.eventlist(x, ...) ## S3 method for class 'list' as.eventlist(x, ...) ## S3 method for class 'data.frame' as.eventlist(x, ...)as.eventlist(x, ...) ## S3 method for class 'list' as.eventlist(x, ...) ## S3 method for class 'data.frame' as.eventlist(x, ...)
x |
list, data.frame |
... |
not used |
Generate objective list from numeric vector
as.objlist(p)as.objlist(p)
p |
Named numeric vector |
list with entries value (0),
gradient (rep(0, length(p))) and
hessian (matrix(0, length(p), length(p))) of class obj.
p <- c(A = 1, B = 2) as.objlist(p)p <- c(A = 1, B = 2) as.objlist(p)
Coerce object to a parameter frame
## S3 method for class 'parlist' as.parframe(x, sort.by = "value", ...) as.parframe(x, ...)## S3 method for class 'parlist' as.parframe(x, sort.by = "value", ...) as.parframe(x, ...)
x |
object to be coerced |
sort.by |
character indicating by which colum the returned parameter frame
should be sorted. Defaults to |
... |
other arguments |
object of class parframe.
## Generate a prediction function regfn <- c(y = "sin(a*time)") g <- Y(regfn, parameters = "a") x <- Xt(condition = "C1") ## Generate data data <- datalist( C1 = data.frame( name = "y", time = 1:5, value = sin(1:5) + rnorm(5, 0, .1), sigma = .1 ) ) ## Initialize parameters and time pars <- c(a = 1) times <- seq(0, 5, .1) plot((g*x)(times, pars), data) ## Do many fits from random positions and store them into parlist out <- as.parlist(lapply(1:50, function(i) { trust(normL2(data, g*x), pars + rnorm(length(pars), 0, 1), rinit = 1, rmax = 10) })) summary(out) ## Reduce parlist to parframe parframe <- as.parframe(out) plotValues(parframe) ## Reduce parframe to best fit bestfit <- as.parvec(parframe) plot((g*x)(times, bestfit), data)## Generate a prediction function regfn <- c(y = "sin(a*time)") g <- Y(regfn, parameters = "a") x <- Xt(condition = "C1") ## Generate data data <- datalist( C1 = data.frame( name = "y", time = 1:5, value = sin(1:5) + rnorm(5, 0, .1), sigma = .1 ) ) ## Initialize parameters and time pars <- c(a = 1) times <- seq(0, 5, .1) plot((g*x)(times, pars), data) ## Do many fits from random positions and store them into parlist out <- as.parlist(lapply(1:50, function(i) { trust(normL2(data, g*x), pars + rnorm(length(pars), 0, 1), rinit = 1, rmax = 10) })) summary(out) ## Reduce parlist to parframe parframe <- as.parframe(out) plotValues(parframe) ## Reduce parframe to best fit bestfit <- as.parvec(parframe) plot((g*x)(times, bestfit), data)
Obtain a parameter vector from a parameter frame.
## S3 method for class 'parframe' as.parvec(x, index = 1, ...)## S3 method for class 'parframe' as.parvec(x, index = 1, ...)
x |
A parameter frame, e.g., the output of
|
index |
Integer, the parameter vector with the |
... |
not used right now |
With this command, additional information included in the parameter frame as the objective value and the convergence state are removed and a parameter vector is returned. This parameter vector can be used to e.g., evaluate an objective function.
On selection, the parameters in the parameter frame are ordered such, that the parameter vector with the lowest objective value is at index 1. Thus, the parameter vector with the index-th lowest objective value is easily obtained.
The parameter vector with the index-th lowest objective value.
Wolfgang Mader, [email protected]
Select or discard attributes from an object.
attrs(x, atr = NULL, keep = TRUE)attrs(x, atr = NULL, keep = TRUE)
x |
The object to work on |
atr |
An optional list of attributes which are either kept or removed. This parameter defaults to dim, dimnames, names, col.names, and row.names. |
keep |
For keep = TRUE, atr is a positive list on attributes which are kept, for keep = FALSE, atr are removed. |
x with selected attributes.
Wolfgang Mader, [email protected]
Mirjam Fehling-Kaschek, [email protected]
Embed two matrices into one blockdiagonal matrix
blockdiagSymb(M, N)blockdiagSymb(M, N)
M |
matrix of type character |
N |
matrix of type character |
Matrix of type character containing M and N as upper left and lower right block
M <- matrix(1:9, 3, 3, dimnames = list(letters[1:3], letters[1:3])) N <- matrix(1:4, 2, 2, dimnames = list(LETTERS[1:2], LETTERS[1:2])) blockdiagSymb(M, N)M <- matrix(1:9, 3, 3, dimnames = list(letters[1:3], letters[1:3])) N <- matrix(1:4, 2, 2, dimnames = list(LETTERS[1:2], LETTERS[1:2])) blockdiagSymb(M, N)
If prediction is NA for observables which are not observed in a condition, they don't matter. In this case, replace NA by 0, such that the error model can be evaluated.
check_and_sanitize_prediction(prediction, data, cn, FLAGNaNInfwarnings)check_and_sanitize_prediction(prediction, data, cn, FLAGNaNInfwarnings)
prediction |
prediction for condition cn |
data |
datalist |
cn |
condition name for which the prediction was made |
FLAGNaNInfwarnings |
print warnings? |
the prediction with harmless NA's replaced by 0
Daniel Lill ([email protected])
Load one row of a dMod.frame into the .GlobalEnv
checkout_hypothesis(dMod.frame, hypothesis, prefix = "", suffix = "")checkout_hypothesis(dMod.frame, hypothesis, prefix = "", suffix = "")
dMod.frame |
A dMod.frame |
hypothesis |
character or numeric. specifying the name or the index of the hypothesis |
prefix |
Prefix appended to the object names in .GlobalEnv |
suffix |
Suffix appended to the object names in .GlobalEnv |
testframe <- dplyr::tibble(hypothesis = c("linear", "quadratic"), plots = list(plot(1:10,1:10), plot(1:10,(1:10)^2)), myfun = list(function(x,a) {a * x}, function(x,a) {a * x^2}), a = c(1:2)) checkout_hypothesis(testframe, "quadratic", prefix = "quad") quadplots quadmyfun(1:10, quada)testframe <- dplyr::tibble(hypothesis = c("linear", "quadratic"), plots = list(plot(1:10,1:10), plot(1:10,(1:10)^2)), myfun = list(function(x,a) {a * x}, function(x,a) {a * x^2}), a = c(1:2)) checkout_hypothesis(testframe, "quadratic", prefix = "quad") quadplots quadmyfun(1:10, quada)
Combine several data.frames by rowbind
combine(...)combine(...)
... |
data.frames or matrices with not necessarily overlapping colnames |
This function is useful when separating models into independent csv model files,
e.g.~a receptor model and several downstream pathways. Then, the models can be recombined
into one model by combine().
A data.frame
data1 <- data.frame(Description = "reaction 1", Rate = "k1*A", A = -1, B = 1) data2 <- data.frame(Description = "reaction 2", Rate = "k2*B", B = -1, C = 1) combine(data1, data2)data1 <- data.frame(Description = "reaction 1", Rate = "k1*A", A = -1, B = 1) data2 <- data.frame(Description = "reaction 2", Rate = "k2*B", B = -1, C = 1) combine(data1, data2)
Works eigher on a list or on two arguments. In case of a list, comparison is done with respect to a reference entry. Besides the objects themselves also some of their attributes are compared, i.e. "equations", "parameters" and "events" and "forcings".
compare(vec1, ...) ## S3 method for class 'list' compare(vec1, vec2 = NULL, reference = 1, ...) ## S3 method for class 'character' compare(vec1, vec2 = NULL, ...) ## S3 method for class 'eqnvec' compare(vec1, vec2 = NULL, ...) ## S3 method for class 'data.frame' compare(vec1, vec2 = NULL, ...)compare(vec1, ...) ## S3 method for class 'list' compare(vec1, vec2 = NULL, reference = 1, ...) ## S3 method for class 'character' compare(vec1, vec2 = NULL, ...) ## S3 method for class 'eqnvec' compare(vec1, vec2 = NULL, ...) ## S3 method for class 'data.frame' compare(vec1, vec2 = NULL, ...)
vec1 |
object of class eqnvec, |
... |
arguments going to the corresponding methods |
vec2 |
same as vec1. Not used if vec1 is a list. |
reference |
numeric of length one, the reference entry. |
data.frame or list of data.frames with the differences.
## Compare equation vectors eq1 <- eqnvec(a = "-k1*a + k2*b", b = "k2*a - k2*b") eq2 <- eqnvec(a = "-k1*a", b = "k2*a - k2*b", c = "k2*b") compare(eq1, eq2) ## Compare character vectors c1 <- c("a", "b") c2 <- c("b", "c") compare(c1, c2) ## Compare data.frames d1 <- data.frame(var = "a", time = 1, value = 1:3, method = "replace") d2 <- data.frame(var = "a", time = 1, value = 2:4, method = "replace") compare(d1, d2) ## Compare structures like prediction functions fn1 <- function(x) x^2 attr(fn1, "equations") <- eq1 attr(fn1, "parameters") <- c1 attr(fn1, "events") <- d1 fn2 <- function(x) x^3 attr(fn2, "equations") <- eq2 attr(fn2, "parameters") <- c2 attr(fn2, "events") <- d2 mylist <- list(f1 = fn1, f2 = fn2) compare(mylist)## Compare equation vectors eq1 <- eqnvec(a = "-k1*a + k2*b", b = "k2*a - k2*b") eq2 <- eqnvec(a = "-k1*a", b = "k2*a - k2*b", c = "k2*b") compare(eq1, eq2) ## Compare character vectors c1 <- c("a", "b") c2 <- c("b", "c") compare(c1, c2) ## Compare data.frames d1 <- data.frame(var = "a", time = 1, value = 1:3, method = "replace") d2 <- data.frame(var = "a", time = 1, value = 2:4, method = "replace") compare(d1, d2) ## Compare structures like prediction functions fn1 <- function(x) x^2 attr(fn1, "equations") <- eq1 attr(fn1, "parameters") <- c1 attr(fn1, "events") <- d1 fn2 <- function(x) x^3 attr(fn2, "equations") <- eq2 attr(fn2, "parameters") <- c2 attr(fn2, "events") <- d2 mylist <- list(f1 = fn1, f2 = fn2) compare(mylist)
Compile one or more prdfn, obsfn or parfn objects
compile(..., output = NULL, args = NULL, cores = 1, verbose = F)compile(..., output = NULL, args = NULL, cores = 1, verbose = F)
... |
Objects of class parfn, obsfn or prdfn |
output |
Optional character of the file to be produced. If several objects were passed, the different C files are all compiled into one shared object file. |
args |
Additional arguments for the R CMD SHLIB call, e.g. |
cores |
Number of cores used for compilation when several files are compiled. |
verbose |
Print compiler output to R command line. |
extract parameter uncertainties from profiles
## S3 method for class 'parframe' confint(object, parm = NULL, level = 0.95, ..., val.column = "data")## S3 method for class 'parframe' confint(object, parm = NULL, level = 0.95, ..., val.column = "data")
object |
object of class |
parm |
a specification of which parameters are to be given confidence intervals, either a vector of numbers or a vector of names. If missing, all parameters are considered. |
level |
the confidence level required. |
... |
not used right now. |
val.column |
the value column used in the parframe, usually 'data'. |
Determine conserved quantites by finding the kernel of the stoichiometric matrix
conservedQuantities(S)conservedQuantities(S)
S |
Stoichiometric matrix |
Data frame with conserved quantities carrying an attribute with the number of conserved quantities.
Malenke Mader, [email protected]
# Write your example here. You can also add more Start..End blocks if needed. # Please mask all output such as print() with the special tag # # such that the test is not littered. Statements guarded by are enabled # in the example file which is extracted from this test file. To extract the # example run # extractExamples() # on the R command line. ## Generate another equation list eq <- eqnlist() eq <- addReaction(eq, "A", "pA", "act_A * A * stimulus", "Phosphorylation of A") eq <- addReaction(eq, "pA", "A", "deact_A * pA", "Deposphorylation of pA") eq <- addReaction(eq, "2*pA", "pA_pA", "form_complex_pA * pA^2", "Complex formation of pA") eq <- addReaction(eq, "B", "pB", "act_B * B * pA_pA", "Phosphorylation of B") eq <- addReaction(eq, "pB", "B", "deact_B * pB", "Deposphorylation of pB") ## Extract data.frame of reactions reactions <- getReactions(eq) print(reactions) ## Get conserved quantities cq <- conservedQuantities(eq$smatrix) print(cq) ## Get fluxes fluxes <- getFluxes(eq) print(fluxes) ## Subsetting of equation list subeq1 <- subset(eq, "pB" %in% Product) print(subeq1) subeq2 <- subset(eq, grepl("not_available", Description)) print(subeq2) ## Time derivatives of observables observables <- eqnvec(pA_obs = "s1*pA", tA_obs = "s2*(A + pA)") dobs <- dot(observables, eq) ## Combined equation vector for ODE and observables f <- c(as.eqnvec(eq), dobs) print(f)# Write your example here. You can also add more Start..End blocks if needed. # Please mask all output such as print() with the special tag # # such that the test is not littered. Statements guarded by are enabled # in the example file which is extracted from this test file. To extract the # example run # extractExamples() # on the R command line. ## Generate another equation list eq <- eqnlist() eq <- addReaction(eq, "A", "pA", "act_A * A * stimulus", "Phosphorylation of A") eq <- addReaction(eq, "pA", "A", "deact_A * pA", "Deposphorylation of pA") eq <- addReaction(eq, "2*pA", "pA_pA", "form_complex_pA * pA^2", "Complex formation of pA") eq <- addReaction(eq, "B", "pB", "act_B * B * pA_pA", "Phosphorylation of B") eq <- addReaction(eq, "pB", "B", "deact_B * pB", "Deposphorylation of pB") ## Extract data.frame of reactions reactions <- getReactions(eq) print(reactions) ## Get conserved quantities cq <- conservedQuantities(eq$smatrix) print(cq) ## Get fluxes fluxes <- getFluxes(eq) print(fluxes) ## Subsetting of equation list subeq1 <- subset(eq, "pB" %in% Product) print(subeq1) subeq2 <- subset(eq, grepl("not_available", Description)) print(subeq2) ## Time derivatives of observables observables <- eqnvec(pA_obs = "s1*pA", tA_obs = "s2*(A + pA)") dobs <- dot(observables, eq) ## Combined equation vector for ODE and observables f <- c(as.eqnvec(eq), dobs) print(f)
Compute a differentiable box prior
constraintExp2(p, mu, sigma = 1, k = 0.05, fixed = NULL)constraintExp2(p, mu, sigma = 1, k = 0.05, fixed = NULL)
p |
Named numeric, the parameter value |
mu |
Named numeric, the prior values, means of boxes |
sigma |
Named numeric, half box width |
k |
Named numeric, shape of box; if 0 a quadratic prior is obtained, the higher k the more box shape, gradient at border of the box (-sigma, sigma) is equal to sigma*k |
fixed |
Named numeric with fixed parameter values (contribute to the prior value but not to gradient and Hessian) |
list with entries: value (numeric, the weighted residual sum of squares),
gradient (numeric, gradient) and
hessian (matrix of type numeric). Object of class objlist.
Soft L2 constraint on parameters
constraintL2(mu, sigma = 1, attr.name = "prior", condition = NULL)constraintL2(mu, sigma = 1, attr.name = "prior", condition = NULL)
mu |
named numeric, the prior values |
sigma |
named numeric of length of mu or numeric of length one or character of length of mu or character of length one |
attr.name |
character. The constraint value is additionally returned in an attributed with this name |
condition |
character, the condition for which the constraint should apply. If
|
If sigma is numeric, the function computes the constraint value
and its derivatives with respect to p. If sigma is a character, the function computes
and its derivatives with respect to p and sigma. Sigma parameters being
passed to the function are ALWAYS assumed to be on a log scale, i.e. internally
sigma parameters are converted by exp().
object of class objfn
wrss
mu <- c(A = 0, B = 0) sigma <- c(A = 0.1, B = 1) myfn <- constraintL2(mu, sigma) myfn(pars = c(A = 1, B = -1)) # Introduce sigma parameter but fix them (sigma parameters # are assumed to be passed on log scale) mu <- c(A = 0, B = 0) sigma <- paste("sigma", names(mu), sep = "_") myfn <- constraintL2(mu, sigma) pars <- c(A = .8, B = -.3, sigma_A = -1, sigma_B = 1) myfn(pars = pars[c(1, 3)], fixed = pars[c(2, 4)]) # Assume same sigma parameter for both A and B # sigma is assumed to be passed on log scale mu <- c(A = 0, B = 0) myfn <- constraintL2(mu, sigma = "sigma") pars <- c(A = .8, B = -.3, sigma = 0) myfn(pars = pars)mu <- c(A = 0, B = 0) sigma <- c(A = 0.1, B = 1) myfn <- constraintL2(mu, sigma) myfn(pars = c(A = 1, B = -1)) # Introduce sigma parameter but fix them (sigma parameters # are assumed to be passed on log scale) mu <- c(A = 0, B = 0) sigma <- paste("sigma", names(mu), sep = "_") myfn <- constraintL2(mu, sigma) pars <- c(A = .8, B = -.3, sigma_A = -1, sigma_B = 1) myfn(pars = pars[c(1, 3)], fixed = pars[c(2, 4)]) # Assume same sigma parameter for both A and B # sigma is assumed to be passed on log scale mu <- c(A = 0, B = 0) myfn <- constraintL2(mu, sigma = "sigma") pars <- c(A = .8, B = -.3, sigma = 0) myfn(pars = pars)
Applies to objects of class objfn,
parfn, prdfn and obsfn. Allows to manipulate
different arguments that have been set when creating the
objects.
controls(x, ...) ## S3 method for class 'objfn' controls(x, name = NULL, ...) ## S3 method for class 'fn' controls(x, condition = NULL, name = NULL, ...) controls(x, ...) <- value ## S3 replacement method for class 'objfn' controls(x, name, ...) <- value ## S3 replacement method for class 'fn' controls(x, condition = NULL, name, ...) <- valuecontrols(x, ...) ## S3 method for class 'objfn' controls(x, name = NULL, ...) ## S3 method for class 'fn' controls(x, condition = NULL, name = NULL, ...) controls(x, ...) <- value ## S3 replacement method for class 'objfn' controls(x, name, ...) <- value ## S3 replacement method for class 'fn' controls(x, condition = NULL, name, ...) <- value
x |
function |
... |
arguments going to the appropriate S3 methods |
name |
character, the name of the control |
condition |
character, the condition name |
value |
the new value |
If called without further arguments, controls(x) lists the
available controls within an object. Calling controls() with name
and condition returns the control value. The value can be overwritten. If
a list or data.frame ist returned, elements of those can be manipulated by the
$- or []-operator.
Either a print-out or the values of the control.
## parfn with condition p <- P(eqnvec(x = "-a*x"), method = "implicit", condition = "C1") controls(p) controls(p, "C1", "keep.root") controls(p, "C1", "keep.root") <- FALSE ## obsfn with NULL condition g <- Y(g = eqnvec(y = "s*x"), f = NULL, states = "x", parameters = "s") controls(g) controls(g, NULL, "attach.input") controls(g, NULL, "attach.input") <- FALSE## parfn with condition p <- P(eqnvec(x = "-a*x"), method = "implicit", condition = "C1") controls(p) controls(p, "C1", "keep.root") controls(p, "C1", "keep.root") <- FALSE ## obsfn with NULL condition g <- Y(g = eqnvec(y = "s*x"), f = NULL, states = "x", parameters = "s") controls(g) controls(g, NULL, "attach.input") controls(g, NULL, "attach.input") <- FALSE
Applies a symbolically defined transformation to the value
column of a data frame. Additionally, if a sigma column is
present, those values are transformed according to Gaussian error
propagation.
coordTransform(data, transformations)coordTransform(data, transformations)
data |
data frame with at least columns "name" (character) and "value" (numeric). Can optionally contain a column "sigma" (numeric). |
transformations |
character (the transformation) or named list of characters. In this case, the list names must be a subset of those contained in the "name" column. |
The data frame with the transformed values and sigma uncertainties.
mydata1 <- data.frame(name = c("A", "B"), time = 0:5, value = 0:5, sigma = .1) coordTransform(mydata1, "log(value)") coordTransform(mydata1, list(A = "exp(value)", B = "sqrt(value)"))mydata1 <- data.frame(name = c("A", "B"), time = 0:5, value = 0:5, sigma = .1) coordTransform(mydata1, "log(value)") coordTransform(mydata1, list(A = "exp(value)", B = "sqrt(value)"))
Access the covariates in the data
## S3 method for class 'tbl_df' covariates(x, hypothesis = 1) covariates(x) ## S3 method for class 'datalist' covariates(x) ## S3 method for class 'data.frame' covariates(x)## S3 method for class 'tbl_df' covariates(x, hypothesis = 1) covariates(x) ## S3 method for class 'datalist' covariates(x) ## S3 method for class 'data.frame' covariates(x)
x |
Either a datalist or a |
hypothesis |
The hypothesis in the dMod.frame |
The condition.grid of the data
Open a unit test template.
createExample(test, testPath = NULL)createExample(test, testPath = NULL)
test |
Name of the unit test used as file name with "test-" prefixed. |
testPath |
Unit test folder. Defaults to "inst/tests". |
Wolfgang Mader, [email protected]
The datalist object stores time-course data in a list of data.frames. The names of the list serve as identifiers, e.g. of an experimental condition, etc.
datalist(...) as.datalist(x, ...) ## S3 method for class 'data.frame' as.datalist(x, split.by = NULL, keep.covariates = NULL, ...) ## S3 method for class 'list' as.datalist(x, names = NULL, ..., condition.grid = attr(x, "condition.grid")) ## S3 replacement method for class 'datalist' names(x) <- value is.datalist(x) ## S3 method for class 'datalist' c(...)datalist(...) as.datalist(x, ...) ## S3 method for class 'data.frame' as.datalist(x, split.by = NULL, keep.covariates = NULL, ...) ## S3 method for class 'list' as.datalist(x, names = NULL, ..., condition.grid = attr(x, "condition.grid")) ## S3 replacement method for class 'datalist' names(x) <- value is.datalist(x) ## S3 method for class 'datalist' c(...)
... |
data.frame objects to be coerced into a list and additional arguments |
x |
object of class |
split.by |
vector of columns names which yield a unique identifier (conditions). If NULL, all columns except for the expected standard columns "name", "time", "value", "sigma" and "lloq" will be selected. |
keep.covariates |
vector of additional column names which should be kept in the condition.grid. |
names |
optional names vector, otherwise names are taken from |
condition.grid |
Optionally, to manually specify a condition.grid |
value |
The new condition names of the datalist and its condition.grid |
Datalists can be plotted, see plotData and merged, see sumdatalist. They are the basic structure when combining model prediction and data via the normL2 objective function.
The standard columns of the datalist data frames are "name" (observable name),
"time" (time points), "value" (data value), "sigma" (uncertainty, can be NA), and
"lloq" (lower limit of quantification, -Inf by default).
Datalists carry the attribute condition.grid which contains additional information about different
conditions, such as dosing information for the experiment. It can be conveniently accessed by the covariates-function.
Reassigning names to a datalist also renames the rows of the condition.grid.
Object of class datalist.
Object of class datalist
## Generate datalist from scratch mydata1 <- data.frame(name = "A", time = 0:5, value = 0:5, sigma = .1, lloq = -0.5) mydata2 <- data.frame(name = "A", time = 0:5, value = sin(0:5), sigma = .1) data <- datalist(C1 = mydata1, C2 = mydata2) print(data) plot(data) ## Generate datalist from singla data.frame times <- seq(0, 2*pi, length.out = 20) mydata <- data.frame(name = "A", time = times, value = c(sin(times), 1.5 * sin(times)), sigma = .1, stage = rep(c("upper", "lower"), each = 10), phase = rep(c("first", "second"), each = 20), amplitude = rep(c(1,1.5), each = 20)) data <- as.datalist(mydata, split.by = c("stage", "phase"), keep.covariates = "amplitude") print(data) plot(data) condition.grid <- attr(data, "condition.grid") print(condition.grid)## Generate datalist from scratch mydata1 <- data.frame(name = "A", time = 0:5, value = 0:5, sigma = .1, lloq = -0.5) mydata2 <- data.frame(name = "A", time = 0:5, value = sin(0:5), sigma = .1) data <- datalist(C1 = mydata1, C2 = mydata2) print(data) plot(data) ## Generate datalist from singla data.frame times <- seq(0, 2*pi, length.out = 20) mydata <- data.frame(name = "A", time = times, value = c(sin(times), 1.5 * sin(times)), sigma = .1, stage = rep(c("upper", "lower"), each = 10), phase = rep(c("first", "second"), each = 20), amplitude = rep(c(1,1.5), each = 20)) data <- as.datalist(mydata, split.by = c("stage", "phase"), keep.covariates = "amplitude") print(data) plot(data) condition.grid <- attr(data, "condition.grid") print(condition.grid)
L2 objective function for validation data point
datapointL2(name, time, value, sigma = 1, attr.name = "validation", condition)datapointL2(name, time, value, sigma = 1, attr.name = "validation", condition)
name |
character, the name of the prediction, e.g. a state name. |
time |
numeric, the time-point associated to the prediction |
value |
character, the name of the parameter which contains the prediction value. |
sigma |
numeric, the uncertainty of the introduced test data point |
attr.name |
character. The constraint value is additionally returned in an attributed with this name |
condition |
character, the condition for which the prediction is made. |
Computes the constraint value
and its derivatives with respect to p.
List of class objlist, i.e. objective value, gradient and Hessian as list.
wrss, constraintL2
prediction <- list(a = matrix(c(0, 1), nrow = 1, dimnames = list(NULL, c("time", "A")))) derivs <- matrix(c(0, 1, 0.1), nrow = 1, dimnames = list(NULL, c("time", "A.A", "A.k1"))) attr(prediction$a, "deriv") <- derivs p0 <- c(A = 1, k1 = 2) vali <- datapointL2(name = "A", time = 0, value = "newpoint", sigma = 1, condition = "a") vali(pars = c(p0, newpoint = 1), env = .GlobalEnv)prediction <- list(a = matrix(c(0, 1), nrow = 1, dimnames = list(NULL, c("time", "A")))) derivs <- matrix(c(0, 1, 0.1), nrow = 1, dimnames = list(NULL, c("time", "A.A", "A.k1"))) attr(prediction$a, "deriv") <- derivs p0 <- c(A = 1, k1 = 2) vali <- datapointL2(name = "A", time = 0, value = "newpoint", sigma = 1, condition = "a") vali(pars = c(p0, newpoint = 1), env = .GlobalEnv)
DatapointL2 without access to stored predictions
datapointL2_indiv( name, time, value, sigma = 1, attr.name = "validation", condition, prd_indiv )datapointL2_indiv( name, time, value, sigma = 1, attr.name = "validation", condition, prd_indiv )
name |
character, the name of the prediction, e.g. a state name. |
time |
numeric, the time-point associated to the prediction |
value |
character, the name of the parameter which contains the prediction value. |
sigma |
numeric, the uncertainty of the introduced test data point |
attr.name |
character. The constraint value is additionally returned in an attributed with this name |
condition |
character, the condition for which the prediction is made. |
prd_indiv |
a prediction function |
Daniel Lill ([email protected])
define(), branch() and insert()
Define parameter transformations by define(), branch() and insert()
define(trafo, expr, ..., conditionMatch = NULL) insert(trafo, expr, ..., conditionMatch = NULL) branch(trafo, table = NULL, conditions = rownames(table))define(trafo, expr, ..., conditionMatch = NULL) insert(trafo, expr, ..., conditionMatch = NULL) branch(trafo, table = NULL, conditions = rownames(table))
trafo |
named character vector of parametric expressions or object
of class |
expr |
character of the form |
... |
used to pass values for symbols as named arguments |
conditionMatch |
optional character, Use as regular expression to apply the reparameterization only to conditions containing conditionMatch |
table |
table of covariates as data frame. Rownames are used as unique identifier, usually called "conditions", and columns represent covariates associated with these conditions. |
conditions |
character vector with condition names. Overwrites the rownames of table. |
object of the same class as trafo or list thereof, if branch() has been
used.
# Define some parameter names parameters <- c("A", "B", "k1", "k2") # Define a covariate table covtable <- data.frame(dose = c(1, 1, 10), inhibitor = c("no", "inh", "no"), row.names = c("Low_noInh", "Low_Inh", "High_noInh")) # Start with an empty transformation trans <- NULL # Generate the identity transformation for parameters trans <- define(trans, "x ~ x", x = parameters); print(trans) # Insert exp(x) wherever you find x trans <- insert(trans, "x ~ exp(x)", x = parameters); print(trans) # Some new expressions instead of k1 and k2 trans <- insert(trans, "x ~ y", x = c("k1", "k2"), y = c("q1 + q2", "q1 - q2")); print(trans) # Define some parameters as 0 trans <- define(trans, "x ~ 0", x = "B"); print(trans) # The parameter name can also be directly used in the formula trans <- insert(trans, "q1 ~ Q"); print(trans) # Replicate the transformation 3 times with the rownames of covtable as list names trans <- branch(trans, table = covtable); print(trans) # Insert the rhs wherever the lhs is found in the transformation # column names of covtable can be used to perform specific replacements # for each transformation trans <- insert(trans, "x ~ x_inh", x = c("Q", "q2"), inh = inhibitor); print(trans) # Also numbers can be inserted trans <- define(trans, "A ~ dose", dose = dose); print(trans) # Turn that into a parameter transformation function p <- P(trans) parnames <- getParameters(p) pars <- rnorm(length(parnames)) names(pars) <- parnames p(pars) # Advanced tricks exploiting the quoting-mechanism when capturing "..." mydataframe <- data.frame( name = rep(letters[1:2], each = 3), value = 1:6, time = rep(1:3, 2), sigma = 0.1, par1 = rep(0:1, each = 3), par2 = rep(9:10, each = 3), par3 = rep(1:3, each = 2), stringsAsFactors = FALSE ) parameters <- c("a", "b", "par1", "par2", "par3") pars_to_insert <- c("par1", "par2") # this would be the usual way when setting up a model # pars_to_insert <- intersect(getParameters(g*x), names(data)) trafo <- define(NULL, "x~x", x = parameters) trafo <- branch(trafo, covariates(as.datalist(mydataframe))) # Trick 1: Access values from covariates()-Table with get/mget. # The names of the parameters which are supplied in the covariates()-table # have to be supplied manually. trafo <- insert(trafo, "name ~ value", value = unlist(mget(pars_to_insert)), name = pars_to_insert) # Trick 2: Access symbols from current condition-specific trafo with .currentSymbols, access # current condition-specific trafo by .currentTrafo # The input passed by the dots is "quoted" (substituted) and eval()'ed in the environment # of the lapply(1:length(conditions), function(i) {}) trafo <- insert(trafo, "x~exp(X)", x = .currentSymbols, X = toupper(.currentSymbols)) # Trick 3: Condition specificity. There are two ways to do this # 1. Apply reparametrization only for specific conditions using Regular Expressions for the # conditionMatch argument. This matches the condition name agains a regex trafo <- define(NULL, "x~x", x = parameters) trafo <- branch(trafo, covariates(as.datalist(mydataframe))) # Conditions starting with 0_9 insert(trafo, "x~x_par3", x = "a", conditionMatch = "^0_9", par3 = par3) # Conditions NOT starting with 0_9 insert(trafo, "x~0", x = "a", conditionMatch = "^(?!0_9)") # 2. Specify conditions by boolean arguments # Conditions which satisfy par1 == 0 insert(trafo, "x~x_par2", par1 == 0, x = parameters, par2 = par2) # Special case: Pass two arguments with the same name. This is only possible if one of them # is logical and the other is not. # Conditions which satisfy par2 == 9 insert(trafo, "x~x_par2", par2 == 9, x = .currentSymbols, par2 = par2)# Define some parameter names parameters <- c("A", "B", "k1", "k2") # Define a covariate table covtable <- data.frame(dose = c(1, 1, 10), inhibitor = c("no", "inh", "no"), row.names = c("Low_noInh", "Low_Inh", "High_noInh")) # Start with an empty transformation trans <- NULL # Generate the identity transformation for parameters trans <- define(trans, "x ~ x", x = parameters); print(trans) # Insert exp(x) wherever you find x trans <- insert(trans, "x ~ exp(x)", x = parameters); print(trans) # Some new expressions instead of k1 and k2 trans <- insert(trans, "x ~ y", x = c("k1", "k2"), y = c("q1 + q2", "q1 - q2")); print(trans) # Define some parameters as 0 trans <- define(trans, "x ~ 0", x = "B"); print(trans) # The parameter name can also be directly used in the formula trans <- insert(trans, "q1 ~ Q"); print(trans) # Replicate the transformation 3 times with the rownames of covtable as list names trans <- branch(trans, table = covtable); print(trans) # Insert the rhs wherever the lhs is found in the transformation # column names of covtable can be used to perform specific replacements # for each transformation trans <- insert(trans, "x ~ x_inh", x = c("Q", "q2"), inh = inhibitor); print(trans) # Also numbers can be inserted trans <- define(trans, "A ~ dose", dose = dose); print(trans) # Turn that into a parameter transformation function p <- P(trans) parnames <- getParameters(p) pars <- rnorm(length(parnames)) names(pars) <- parnames p(pars) # Advanced tricks exploiting the quoting-mechanism when capturing "..." mydataframe <- data.frame( name = rep(letters[1:2], each = 3), value = 1:6, time = rep(1:3, 2), sigma = 0.1, par1 = rep(0:1, each = 3), par2 = rep(9:10, each = 3), par3 = rep(1:3, each = 2), stringsAsFactors = FALSE ) parameters <- c("a", "b", "par1", "par2", "par3") pars_to_insert <- c("par1", "par2") # this would be the usual way when setting up a model # pars_to_insert <- intersect(getParameters(g*x), names(data)) trafo <- define(NULL, "x~x", x = parameters) trafo <- branch(trafo, covariates(as.datalist(mydataframe))) # Trick 1: Access values from covariates()-Table with get/mget. # The names of the parameters which are supplied in the covariates()-table # have to be supplied manually. trafo <- insert(trafo, "name ~ value", value = unlist(mget(pars_to_insert)), name = pars_to_insert) # Trick 2: Access symbols from current condition-specific trafo with .currentSymbols, access # current condition-specific trafo by .currentTrafo # The input passed by the dots is "quoted" (substituted) and eval()'ed in the environment # of the lapply(1:length(conditions), function(i) {}) trafo <- insert(trafo, "x~exp(X)", x = .currentSymbols, X = toupper(.currentSymbols)) # Trick 3: Condition specificity. There are two ways to do this # 1. Apply reparametrization only for specific conditions using Regular Expressions for the # conditionMatch argument. This matches the condition name agains a regex trafo <- define(NULL, "x~x", x = parameters) trafo <- branch(trafo, covariates(as.datalist(mydataframe))) # Conditions starting with 0_9 insert(trafo, "x~x_par3", x = "a", conditionMatch = "^0_9", par3 = par3) # Conditions NOT starting with 0_9 insert(trafo, "x~0", x = "a", conditionMatch = "^(?!0_9)") # 2. Specify conditions by boolean arguments # Conditions which satisfy par1 == 0 insert(trafo, "x~x_par2", par1 == 0, x = parameters, par2 = par2) # Special case: Pass two arguments with the same name. This is only possible if one of them # is logical and the other is not. # Conditions which satisfy par2 == 9 insert(trafo, "x~x_par2", par2 == 9, x = .currentSymbols, par2 = par2)
Read /proc/loadavg and subtract from the number of cores
detectFreeCores(machine = NULL)detectFreeCores(machine = NULL)
machine |
character, e.g. "user@localhost". |
Generates R and bash scrips, transfers them to remote via ssh. 'sshpass' needs to be installed on your local machine to circumvent password entry.
distributed_computing( ..., jobname, partition = "single", cores = 16, nodes = 1, walltime = "01:00:00", ssh_passwd = NULL, machine = "cluster", var_values = NULL, no_rep = NULL, recover = T, purge_local = F, compile = F, custom_folders = NULL, resetSeeds = TRUE, returnAll = TRUE )distributed_computing( ..., jobname, partition = "single", cores = 16, nodes = 1, walltime = "01:00:00", ssh_passwd = NULL, machine = "cluster", var_values = NULL, no_rep = NULL, recover = T, purge_local = F, compile = F, custom_folders = NULL, resetSeeds = TRUE, returnAll = TRUE )
... |
R code to be remotely executed, parameters to be changed for runs
must be named |
jobname |
Name in characters for the run must be unique, existing will be overwritten. Must not contain the string "Minus". |
partition |
Define the partition used for the SLURM manager. Default is "single". Should only be changed if explicitly wanted. |
cores |
Number of cores per node to be used. If a value is set to 16 < value < 25, the number of possible nodes isseverely limited. 16 or less is possible on all nodes, more on none. |
nodes |
Nodes per task. Default is 1, should not be changed since
|
walltime |
Estimated runtime in the format |
ssh_passwd |
To be set when the sshpass should be used to automatically authenticate on remote via passphrase. This is an obvious security nightmare... |
machine |
SSH address in the form |
var_values |
List of parameter arrays. The number of arrays (i.e. list
entries) must correspond to the number of parameters in the function passed
to |
no_rep |
Number of repetitions. Usage of this parameter and |
recover |
Logical parameter, if set to |
purge_local |
Logical, if set to |
compile |
Logical, if set to |
custom_folders |
named vector with exact three entries named 'compiled',
'output' and 'tmp'. The values are strings with relative paths from the current
working directory to the respective directory of the compiled files, the temporary
folder from which files will be copied to the cluster and the output folder in
which the calculated result from the cluster will be saved.
The default is |
resetSeeds |
logical, if set to |
returnAll |
logical if set to |
distribute_computing() generates R and bash scrips designed to run
on a HPC system managed by the SLURM batch manager. The current workspace
together with the scripts are exported and transferred to remote via SSH. If
ssh-key authentication is not possible, the ssh-password can be passed and
is used by sshpass (which has to be installed on the local machine).
The code to be executed remotely is passed to the ... argument, its final
output is stored in cluster_result, which is loaded in the local
workspace by the get() function.
It is possible to either run repetitions of the same program realization (by
use of the no_rep parameter), or to pass a list of parameter arrays
via var_values. The parameters to be changed for every run *must* be
named var_i where i corresponds to the i-th array in the var_values
parameter.
List of functions check(), get() and purge().
check() checks, if the result is ready.
get() copies all files from the remote working directory to local and
loads all present results (even if not all nodes where done) in the current
active workspace as the object cluster_result which is a list with the
results of each node as entries.
purge() deletes the temporary folder on remote and if purge_local
is set to TRUE.
## Not run: out_distributed_computing <- distributed_computing( { mstrust( objfun=objective_function, center=outer_pars, studyname = "study", rinit = 1, rmax = 10, fits = 48, cores = 16, iterlim = 700, sd = 4 ) }, jobname = "my_name", partition = "single", cores = 16, nodes = 1, walltime = "02:00:00", ssh_passwd = "password", machine = "cluster", var_values = NULL, no_rep = 20, recover = F, compile = F ) out_distributed_computing$check() out_distributed_computing$get() out_distributed_computing$purge() result <- cluster_result print(result) # calculate profiles var_list <- profile_pars_per_node(best_fit, 4) profile_jobname <- paste0(fit_filename,"_profiles_opt") method <- "optimize" profiles_distributed_computing <- distributed_computing( { profile( obj = obj, pars = best_fit, whichPar = (as.numeric(var_1):as.numeric(var_2)), limits = c(-5, 5), cores = 16, method = method, stepControl = list( stepsize = 1e-6, min = 1e-4, max = Inf, atol = 1e-2, rtol = 1e-2, limit = 100 ), optControl = list(iterlim = 20) ) }, jobname = profile_jobname, partition = "single", cores = 16, nodes = 1, walltime = "02:00:00", ssh_passwd = "password", machine = "cluster", var_values = var_list, no_rep = NULL, recover = F, compile = F ) profiles_distributed_computing$check() profiles_distributed_computing$get() profiles_distributed_computing$purge() profiles <- NULL for (i in cluster_result) { profiles <- rbind(profiles, i) } ## End(Not run)## Not run: out_distributed_computing <- distributed_computing( { mstrust( objfun=objective_function, center=outer_pars, studyname = "study", rinit = 1, rmax = 10, fits = 48, cores = 16, iterlim = 700, sd = 4 ) }, jobname = "my_name", partition = "single", cores = 16, nodes = 1, walltime = "02:00:00", ssh_passwd = "password", machine = "cluster", var_values = NULL, no_rep = 20, recover = F, compile = F ) out_distributed_computing$check() out_distributed_computing$get() out_distributed_computing$purge() result <- cluster_result print(result) # calculate profiles var_list <- profile_pars_per_node(best_fit, 4) profile_jobname <- paste0(fit_filename,"_profiles_opt") method <- "optimize" profiles_distributed_computing <- distributed_computing( { profile( obj = obj, pars = best_fit, whichPar = (as.numeric(var_1):as.numeric(var_2)), limits = c(-5, 5), cores = 16, method = method, stepControl = list( stepsize = 1e-6, min = 1e-4, max = Inf, atol = 1e-2, rtol = 1e-2, limit = 100 ), optControl = list(iterlim = 20) ) }, jobname = profile_jobname, partition = "single", cores = 16, nodes = 1, walltime = "02:00:00", ssh_passwd = "password", machine = "cluster", var_values = var_list, no_rep = NULL, recover = F, compile = F ) profiles_distributed_computing$check() profiles_distributed_computing$get() profiles_distributed_computing$purge() profiles <- NULL for (i in cluster_result) { profiles <- rbind(profiles, i) } ## End(Not run)
Basically, a dMod.frame is a tibble, which is grouped rowwise.
Since the dMod.frame is also designed for interactive use, the class will not be called "dMod.frame" as I initially planned, but will be c("tbl_df"). This way, I don't have to overwrite all the dplyr-verbs.
The dMod.frame object stores all objects that are needed to reproducibly produce a result, such as a plot or profiles, in one row of a tibble. Each row corresponds to a distinct hypothesis, e.g. one would have two distinct rows for fitting with and without a prior.
dMod.frame(hypothesis, g, x, p, data, e = NULL, ...)dMod.frame(hypothesis, g, x, p, data, e = NULL, ...)
hypothesis |
Character. Description of the hypothesis |
g |
fn |
x |
fn |
p |
fn |
data |
data.frame or datalist, will be coerced to datalist |
e |
fn |
... |
other columns, have to be named. They will be placed in a list. |
Object of class tbl_df, which is grouped rowwise.
## Not run: # ---- Example library(dMod) library(dplyr) # devtools::install_github("dlill/conveniencefunctions") setwd(tempdir()) ## Model definition (text-based, scripting part) f <- NULL %>% addReaction("A", "B", "k1*A", "translation") %>% addReaction("B", "", "k2*B", "degradation") %>% as.eqnvec() events <- eventlist(var = "A", time = 5, value = "A_add", method = "add") x <- odemodel(f, events = events) %>% Xs g <- Y(c(Bobs = "s1*B"), x, compile = T, modelname = "obsfn") conditions <- c("a", "b") # there is a bug in # getParameters(g*x) parameters <- union(getParameters(g), getParameters(x)) trafo <- NULL %>% define("x~x", x = parameters) %>% branch(conditions = conditions) %>% insert("x~x_cond", x = "s1", cond = condition) %>% insert("x~1", x = "added", conditionMatch = "a") %>% insert("x~5", x = "added", conditionMatch = "b") %>% insert("x~exp(x)", x = getSymbols(mytrafo[[i]])) %>% {.} p <- P(trafo) # Parameter transformation for steady states pSS <- P(f, method = "implicit", condition = NULL, compile = T) ## Process data # Data data <- datalist( a = data.frame(time = c(0, 2, 7), name = "Bobs", value = c(.3, .3, .3), sigma = c(.03, .03, .03)), b = data.frame(time = c(0, 2, 7), name = "Bobs", value = c(.1, .1, .2), sigma = c(.01, .01, .02)) ) # construct dMod.frame myframe1 <- dMod.frame("no steady states", g, x, p, data) print(myframe1) # Augment by derived objects: prd, obj_data, obj, times, pars set.seed(4) myframe2 <- appendObj(myframe1) print(myframe2) # Plot the model with random pars plotCombined(myframe2) # Fit with prior myframe3 <- myframe2 %>% mutate(constr = list(constraintL2(mu = 0*pars, sigma = 5)), obj = list(obj_data + constr)) myframe4 <- myframe3 %>% mutate(fits = list(mstrust(obj, pars, studyname = "Fits", fits = 20, cores = 4, blather = T))) # Inspect Fits myframe5 <- myframe4 %>% appendParframes # Visualize # Little Bug: If you want to use the dots that go to subset(), you need to specify the others plotCombined(myframe5, 1, 1, str_detect(name, "B")) plotPars(myframe5) plotValues(myframe5) plotValues(myframe5, 1, tol =0.0000000001, value >1) # Profiles myframe6 <- myframe5 %>% mutate(profiles = list(profile(obj, as.parvec(parframes), whichPar = "k1"))) myframe6$profiles %>% plotProfile() # Validation myframe7 <- myframe6 %>% mutate(vali = list(datapointL2("A", 2, "mypoint", .1, condition = "a")), obj_vali = list(obj_data + constr + vali), par_vali = list(c(dMod:::sanitizePars(as.parvec(parframes))$pars, "mypoint" = 0.1 )), fits_vali = list(mstrust(obj_vali, par_vali)), profile_vali = list(profile(obj_vali, fits_vali %>% as.parframe %>% as.parvec, "mypoint"))) myframe7$profile_vali %>% plotProfile() # Saving ---- # Save dMod.frame as rds, stage the rds and the .so's it depends on for commit system("git status") git_add_dMod.frame(myframe7) system("git status") # Alternatetively: zip_dMod.frame(myframe7) # Look into you current wd to find the dMod.frame and the .so's zipped into a file. # # Short version of above ---- # # myframe <- dMod.frame("no steady states", g, x, p, data) %>% # appendObj() %>% # mutate(constr = list(constraintL2(mu = 0*pars, sigma = 5)), # obj = list(obj_data + constr), # fits = list(mstrust(obj, pars, studyname = "Fits", fits = 20, cores = 4, blather = T))) %>% # appendParframes() %>% # mutate(profiles = list(profile(obj, as.parvec(parframes), whichPar = "k1"))) %>% # mutate(vali = list(datapointL2("A", 2, "mypoint", .1, condition = "a")), # obj_vali = list(obj_data + constr + vali), # par_vali = list(c(dMod:::sanitizePars(as.parvec(parframes))$pars, "mypoint" = 0.1 )), # fits_vali = list(mstrust(obj_vali, par_vali)), # profile_vali = list(profile(obj_vali, fits_vali %>% as.parframe %>% as.parvec, "mypoint"))) # Easily test several hypotheses ---- # Fit with various prior strengths multiframe <- dMod.frame("no steady states", g, x, p, data) %>% appendObj() %>% rbind(.,.,.,.) %>% # replicate four times ungroup() %>% # If you don't ungroup and run a mutate with an "lapply(1:nrow(), function(i) ...", the index i always gets restored to 1, as mutate() runs the expressions (...) independently for each group. mutate(constr = map(seq_along(x), function(i) constraintL2(mu = 0*pars[[i]], sigma = 10^(i-3))), hypothesis = map_chr(seq_along(x), function(i) paste0(hypothesis[[i]], ", prior sigma = ", 10^(i-3)))) %>% rowwise() %>% #regroup by row for convenient interface to mutate() mutate(obj = list(obj_data + constr), fits = list(mstrust(obj, pars, studyname = "Fits", fits = 10, cores = 4, blather = T))) %>% appendParframes() ## Plot map(1:4, function(i) multiframe %>% plotValues(hypothesis = i)) map(1:4, function(i) multiframe %>% plotCombined(hypothesis = i)) ## Profiles multiframe <- multiframe %>% mutate(profiles = list(profile(obj, parframes %>% as.parvec, names(pars), cores = 4))) ## Plot profiles: The "profiles"-column is already a proflist :) multiframe$profiles %>% plotProfile() + coord_cartesian(ylim = c(-0.5, 4), xlim = c(-2,2)) # Quick and dirty analysis of one single hypothesis ---- checkout_hypothesis(multiframe, 4, suffix = "_weak_prior") parframes_weak_prior ## End(Not run)## Not run: # ---- Example library(dMod) library(dplyr) # devtools::install_github("dlill/conveniencefunctions") setwd(tempdir()) ## Model definition (text-based, scripting part) f <- NULL %>% addReaction("A", "B", "k1*A", "translation") %>% addReaction("B", "", "k2*B", "degradation") %>% as.eqnvec() events <- eventlist(var = "A", time = 5, value = "A_add", method = "add") x <- odemodel(f, events = events) %>% Xs g <- Y(c(Bobs = "s1*B"), x, compile = T, modelname = "obsfn") conditions <- c("a", "b") # there is a bug in # getParameters(g*x) parameters <- union(getParameters(g), getParameters(x)) trafo <- NULL %>% define("x~x", x = parameters) %>% branch(conditions = conditions) %>% insert("x~x_cond", x = "s1", cond = condition) %>% insert("x~1", x = "added", conditionMatch = "a") %>% insert("x~5", x = "added", conditionMatch = "b") %>% insert("x~exp(x)", x = getSymbols(mytrafo[[i]])) %>% {.} p <- P(trafo) # Parameter transformation for steady states pSS <- P(f, method = "implicit", condition = NULL, compile = T) ## Process data # Data data <- datalist( a = data.frame(time = c(0, 2, 7), name = "Bobs", value = c(.3, .3, .3), sigma = c(.03, .03, .03)), b = data.frame(time = c(0, 2, 7), name = "Bobs", value = c(.1, .1, .2), sigma = c(.01, .01, .02)) ) # construct dMod.frame myframe1 <- dMod.frame("no steady states", g, x, p, data) print(myframe1) # Augment by derived objects: prd, obj_data, obj, times, pars set.seed(4) myframe2 <- appendObj(myframe1) print(myframe2) # Plot the model with random pars plotCombined(myframe2) # Fit with prior myframe3 <- myframe2 %>% mutate(constr = list(constraintL2(mu = 0*pars, sigma = 5)), obj = list(obj_data + constr)) myframe4 <- myframe3 %>% mutate(fits = list(mstrust(obj, pars, studyname = "Fits", fits = 20, cores = 4, blather = T))) # Inspect Fits myframe5 <- myframe4 %>% appendParframes # Visualize # Little Bug: If you want to use the dots that go to subset(), you need to specify the others plotCombined(myframe5, 1, 1, str_detect(name, "B")) plotPars(myframe5) plotValues(myframe5) plotValues(myframe5, 1, tol =0.0000000001, value >1) # Profiles myframe6 <- myframe5 %>% mutate(profiles = list(profile(obj, as.parvec(parframes), whichPar = "k1"))) myframe6$profiles %>% plotProfile() # Validation myframe7 <- myframe6 %>% mutate(vali = list(datapointL2("A", 2, "mypoint", .1, condition = "a")), obj_vali = list(obj_data + constr + vali), par_vali = list(c(dMod:::sanitizePars(as.parvec(parframes))$pars, "mypoint" = 0.1 )), fits_vali = list(mstrust(obj_vali, par_vali)), profile_vali = list(profile(obj_vali, fits_vali %>% as.parframe %>% as.parvec, "mypoint"))) myframe7$profile_vali %>% plotProfile() # Saving ---- # Save dMod.frame as rds, stage the rds and the .so's it depends on for commit system("git status") git_add_dMod.frame(myframe7) system("git status") # Alternatetively: zip_dMod.frame(myframe7) # Look into you current wd to find the dMod.frame and the .so's zipped into a file. # # Short version of above ---- # # myframe <- dMod.frame("no steady states", g, x, p, data) %>% # appendObj() %>% # mutate(constr = list(constraintL2(mu = 0*pars, sigma = 5)), # obj = list(obj_data + constr), # fits = list(mstrust(obj, pars, studyname = "Fits", fits = 20, cores = 4, blather = T))) %>% # appendParframes() %>% # mutate(profiles = list(profile(obj, as.parvec(parframes), whichPar = "k1"))) %>% # mutate(vali = list(datapointL2("A", 2, "mypoint", .1, condition = "a")), # obj_vali = list(obj_data + constr + vali), # par_vali = list(c(dMod:::sanitizePars(as.parvec(parframes))$pars, "mypoint" = 0.1 )), # fits_vali = list(mstrust(obj_vali, par_vali)), # profile_vali = list(profile(obj_vali, fits_vali %>% as.parframe %>% as.parvec, "mypoint"))) # Easily test several hypotheses ---- # Fit with various prior strengths multiframe <- dMod.frame("no steady states", g, x, p, data) %>% appendObj() %>% rbind(.,.,.,.) %>% # replicate four times ungroup() %>% # If you don't ungroup and run a mutate with an "lapply(1:nrow(), function(i) ...", the index i always gets restored to 1, as mutate() runs the expressions (...) independently for each group. mutate(constr = map(seq_along(x), function(i) constraintL2(mu = 0*pars[[i]], sigma = 10^(i-3))), hypothesis = map_chr(seq_along(x), function(i) paste0(hypothesis[[i]], ", prior sigma = ", 10^(i-3)))) %>% rowwise() %>% #regroup by row for convenient interface to mutate() mutate(obj = list(obj_data + constr), fits = list(mstrust(obj, pars, studyname = "Fits", fits = 10, cores = 4, blather = T))) %>% appendParframes() ## Plot map(1:4, function(i) multiframe %>% plotValues(hypothesis = i)) map(1:4, function(i) multiframe %>% plotCombined(hypothesis = i)) ## Profiles multiframe <- multiframe %>% mutate(profiles = list(profile(obj, parframes %>% as.parvec, names(pars), cores = 4))) ## Plot profiles: The "profiles"-column is already a proflist :) multiframe$profiles %>% plotProfile() + coord_cartesian(ylim = c(-0.5, 4), xlim = c(-2,2)) # Quick and dirty analysis of one single hypothesis ---- checkout_hypothesis(multiframe, 4, suffix = "_weak_prior") parframes_weak_prior ## End(Not run)
The time evolution of the internal states is defined in the equation list. Time derivatives of observation functions are expressed in terms of the rates of the internal states.
dot(observable, eqnlist)dot(observable, eqnlist)
observable |
named character vector or object of type eqnvec |
eqnlist |
equation list |
Observables are translated into an ODE
An object of class eqnvec
# Write your example here. You can also add more Start..End blocks if needed. # Please mask all output such as print() with the special tag # # such that the test is not littered. Statements guarded by are enabled # in the example file which is extracted from this test file. To extract the # example run # extractExamples() # on the R command line. ## Generate another equation list eq <- eqnlist() eq <- addReaction(eq, "A", "pA", "act_A * A * stimulus", "Phosphorylation of A") eq <- addReaction(eq, "pA", "A", "deact_A * pA", "Deposphorylation of pA") eq <- addReaction(eq, "2*pA", "pA_pA", "form_complex_pA * pA^2", "Complex formation of pA") eq <- addReaction(eq, "B", "pB", "act_B * B * pA_pA", "Phosphorylation of B") eq <- addReaction(eq, "pB", "B", "deact_B * pB", "Deposphorylation of pB") ## Extract data.frame of reactions reactions <- getReactions(eq) print(reactions) ## Get conserved quantities cq <- conservedQuantities(eq$smatrix) print(cq) ## Get fluxes fluxes <- getFluxes(eq) print(fluxes) ## Subsetting of equation list subeq1 <- subset(eq, "pB" %in% Product) print(subeq1) subeq2 <- subset(eq, grepl("not_available", Description)) print(subeq2) ## Time derivatives of observables observables <- eqnvec(pA_obs = "s1*pA", tA_obs = "s2*(A + pA)") dobs <- dot(observables, eq) ## Combined equation vector for ODE and observables f <- c(as.eqnvec(eq), dobs) print(f)# Write your example here. You can also add more Start..End blocks if needed. # Please mask all output such as print() with the special tag # # such that the test is not littered. Statements guarded by are enabled # in the example file which is extracted from this test file. To extract the # example run # extractExamples() # on the R command line. ## Generate another equation list eq <- eqnlist() eq <- addReaction(eq, "A", "pA", "act_A * A * stimulus", "Phosphorylation of A") eq <- addReaction(eq, "pA", "A", "deact_A * pA", "Deposphorylation of pA") eq <- addReaction(eq, "2*pA", "pA_pA", "form_complex_pA * pA^2", "Complex formation of pA") eq <- addReaction(eq, "B", "pB", "act_B * B * pA_pA", "Phosphorylation of B") eq <- addReaction(eq, "pB", "B", "deact_B * pB", "Deposphorylation of pB") ## Extract data.frame of reactions reactions <- getReactions(eq) print(reactions) ## Get conserved quantities cq <- conservedQuantities(eq$smatrix) print(cq) ## Get fluxes fluxes <- getFluxes(eq) print(fluxes) ## Subsetting of equation list subeq1 <- subset(eq, "pB" %in% Product) print(subeq1) subeq2 <- subset(eq, grepl("not_available", Description)) print(subeq2) ## Time derivatives of observables observables <- eqnvec(pA_obs = "s1*pA", tA_obs = "s2*(A + pA)") dobs <- dot(observables, eq) ## Combined equation vector for ODE and observables f <- c(as.eqnvec(eq), dobs) print(f)
The eqnlist object stores an ODE as a list of stoichiometric matrix, rate expressions, state names and compartment volumes.
Translates a reaction network, e.g. defined by a data.frame, into an equation list object.
eqnlist( smatrix = NULL, states = colnames(smatrix), rates = NULL, volumes = NULL, description = NULL ) as.eqnlist(data, volumes) ## S3 method for class 'data.frame' as.eqnlist(data, volumes = NULL) is.eqnlist(x)eqnlist( smatrix = NULL, states = colnames(smatrix), rates = NULL, volumes = NULL, description = NULL ) as.eqnlist(data, volumes) ## S3 method for class 'data.frame' as.eqnlist(data, volumes = NULL) is.eqnlist(x)
smatrix |
Matrix of class numeric. The stoichiometric matrix, one row per reaction/process and one column per state. |
states |
Character vector. Names of the states. |
rates |
Character vector. The rate expressions. |
volumes |
Named character, volume parameters for states. Names must be a subset of the states.
Values can be either characters, e.g. "V1", or numeric values for the volume. If |
description |
Character vector. Description of the single processes. |
data |
data.frame with columns Description, Rate, and one colum for each state reflecting the stoichiometric matrix |
x |
object of class |
... |
additional arguments to be passed to or from methods. |
If data is a data.frame, it must contain columns "Description" (character),
"Rate" (character), and one column per ODE state with the state names.
The state columns correspond to the stoichiometric matrix.
An object of class eqnlist, basically a list.
Object of class eqnlist
# Generate eqnlist from the constructor S <- matrix(c(-1, 1, 1, -1), nrow = 2, ncol = 2, dimnames = list(NULL, c("A", "B"))) rates <- c("k1*A", "k2*B") description <- c("forward", "backward") f <- eqnlist(smatrix = S, rates = rates, description = description) print(f) # Convert to data.frame fdata <- as.data.frame(f) print(fdata) # Generate eqnlist from data.frame and add volume parameter f <- as.eqnlist(fdata, volumes = c(A = "Vcyt", B = "Vnuc")) print(f) print(as.eqnvec(f)) print(as.eqnvec(f, type = "amount"))# Generate eqnlist from the constructor S <- matrix(c(-1, 1, 1, -1), nrow = 2, ncol = 2, dimnames = list(NULL, c("A", "B"))) rates <- c("k1*A", "k2*B") description <- c("forward", "backward") f <- eqnlist(smatrix = S, rates = rates, description = description) print(f) # Convert to data.frame fdata <- as.data.frame(f) print(fdata) # Generate eqnlist from data.frame and add volume parameter f <- as.eqnlist(fdata, volumes = c(A = "Vcyt", B = "Vnuc")) print(f) print(as.eqnvec(f)) print(as.eqnvec(f, type = "amount"))
The eqnvec object stores explicit algebraic equations, like the right-hand sides of an ODE, observation functions or parameter transformations as named character vectors.
eqnvec(...) is.eqnvec(x)eqnvec(...) is.eqnvec(x)
... |
mathematical expressions as characters to be coerced, the right-hand sides of the equations |
x |
obect of any class |
object of class eqnvec, basically a named character.
v <- eqnvec(y = "2*x + offset") print(v) is.eqnvec(v)v <- eqnvec(y = "2*x + offset") print(v) is.eqnvec(v)
An eventlist is a data.frame with the necessary parameters to define an event as columns and specific events as rows. Event time and value can be passed as parameters, which can also be estimated.
eventlist(var = NULL, time = NULL, value = NULL, root = NULL, method = NULL) addEvent(event, var, time = 0, value = 0, root = NA, method = "replace", ...)eventlist(var = NULL, time = NULL, value = NULL, root = NULL, method = NULL) addEvent(event, var, time = 0, value = 0, root = NA, method = "replace", ...)
var |
Character, the state to which the event is applied |
time |
Character or Numeric, the time at which the event happens |
value |
Character or Numeric, the value of the event |
root |
Character or NA, condition to trigger the event instead of time |
method |
Character, options are "replace", "add" or "multiply" |
event |
object of class |
... |
not used |
The function addEvent is pipe-friendly
data.frame with class eventlist
eventlist(var = "A", time = "5", value = 1, method = "add") events <- addEvent(NULL, var = "A", time = "5", value = 1, method = "add") events <- addEvent(events, var = "A", time = "10", value = 1, method = "add") # With symbols events <- eventlist() # Set "A" to "value_switch" at time "time_Switch" events <- addEvent(events, var = "A", time = "time_switch", value = "value_switch", method = "replace") # Set "B" to 2 when "A" reaches "A_target". The time-parameter for internal use will be "time_root". events <- addEvent(events, var = "B", time = "time_root", value = 2, root = "A - A_target", method = "replace")eventlist(var = "A", time = "5", value = 1, method = "add") events <- addEvent(NULL, var = "A", time = "5", value = 1, method = "add") events <- addEvent(events, var = "A", time = "10", value = 1, method = "add") # With symbols events <- eventlist() # Set "A" to "value_switch" at time "time_Switch" events <- addEvent(events, var = "A", time = "time_switch", value = "value_switch", method = "replace") # Set "B" to 2 when "A" reaches "A_target". The time-parameter for internal use will be "time_root". events <- addEvent(events, var = "B", time = "time_root", value = 2, root = "A - A_target", method = "replace")
Extract an example from a unit test file.
exmpextr(test, testPath = NULL, examplePath = NULL)exmpextr(test, testPath = NULL, examplePath = NULL)
test |
File name of the test file, source. |
testPath |
Path to test files. Defaults to "inst/tests". |
examplePath |
Path under which the example is stored. Defaults to "inst/examples". |
From test, an example is extracted and written to a file of the same name, but with removed "test-". The file is saved under examplePath. The special tag #-! can be used to hide code in the test which is enabled again in the example.
Wolfgang Mader, [email protected]
Alternative version of expand.grid
expand.grid.alt(seq1, seq2)expand.grid.alt(seq1, seq2)
seq1 |
Vector, numeric or character |
seq2 |
Vector, numeric or character |
Matrix ob combinations of elemens of seq1 and seq2
Extract example from unit tests.
extractExamples(testPath = NULL, examplePath = NULL)extractExamples(testPath = NULL, examplePath = NULL)
testPath |
character, path to test folder |
examplePath |
character path to example folder |
Calls exmpextr for all unit test files.
Wolfgang Mader, [email protected]
Fit an error model to reduced replicate data, see
reduceReplicates.
fitErrorModel( data, factors, errorModel = "exp(s0)+exp(srel)*x^2", par = c(s0 = 1, srel = 0.1), plotting = TRUE, blather = FALSE, ... )fitErrorModel( data, factors, errorModel = "exp(s0)+exp(srel)*x^2", par = c(s0 = 1, srel = 0.1), plotting = TRUE, blather = FALSE, ... )
data |
Reduced replicate data, see |
factors |
data is pooled with respect to the columns named here, see Details. |
errorModel |
Character vector defining the error model in terms of the variance. Use x to reference the independend variable, see Details. |
par |
Inital values for the parameters of the error model. |
plotting |
If TRUE, a plot of the pooled variance together with the fit of the error model is shown. |
blather |
If TRUE, additional information is returned, such as fit parameters and sigmaLS (original sigma given in input data). |
... |
Parameters handed to the optimizer |
The variance estimator using data points is
distributed with degrees of freedom. Given replicates for
consecutive time points, the sample variance can be assumed a function of
the sample mean. By defining an error model which must hold for all time
points, a maximum likelihood estimator for the parameters of the error
model can be derived. The parameter errorModel takes the error
model as a character vector, where the mean (independent variable) is
refered to as x.
It is desireable to estimate the variance from many replicates. The parameter data must provide one or more columns which define the pooling of data. In case more than one column is announced by factors, all combinations are constructed. If, e.g., factors = c("condition", "name") is used, where "condition" is "a", "b", "c" and repeating and "name" is "d", "e" and repeating, the effective conditions used for pooling are "a d", "b e", "c d", "a e", "b d", and "c e".
By default, a plot of the pooled data, sigma and its confidence bound at 68% and 95% is shown.
Returned by default is a data frame with columns as in data, but with the sigma values replaced by the derived values, obtained by evaluating the error model with the fit parameters.
If the blather = TRUE option is chosen, fit values of the parameters of the error model are appended, with the column names equal to the parameter names. The error model is appended as the attribute "errorModel". Confidence bounds for sigma at confidence level 68% and 95% are calculated, their values come next in the returned data frame. Finally, the effective conditions are appended to easily check how the pooling was done.
Wolfgang Mader, [email protected]
Return some useful forcing functions as strings
forcingsSymb( type = c("Gauss", "Fermi", "1-Fermi", "MM", "Signal", "Dose"), parameters = NULL )forcingsSymb( type = c("Gauss", "Fermi", "1-Fermi", "MM", "Signal", "Dose"), parameters = NULL )
type |
Which function to be returned |
parameters |
Named vector, character or numeric. Replace parameters by the corresponding valus
in |
String with the function
Encode equation vector in format with sufficient spaces
## S3 method for class 'eqnvec' format(x, ...)## S3 method for class 'eqnvec' format(x, ...)
x |
object of class eqnvec. Alternatively, a named parsable character vector. |
... |
additional arguments |
named character
Evaluation of algebraic expressions defined by characters
funC0( x, variables = getSymbols(x, exclude = parameters), parameters = NULL, compile = FALSE, modelname = NULL, verbose = FALSE, convenient = TRUE, warnings = TRUE )funC0( x, variables = getSymbols(x, exclude = parameters), parameters = NULL, compile = FALSE, modelname = NULL, verbose = FALSE, convenient = TRUE, warnings = TRUE )
x |
Object of class |
variables |
character vector, the symbols that should be treated as variables |
parameters |
character vector, the symbols that should be treated as parameters |
compile |
Logical. Directly compile the file. If |
modelname |
file name of the generated C file. See description of parameter |
verbose |
Print compiler output to R command line. |
convenient |
logical, if TRUE return a function with argument |
warnings |
logical. Suppress warnings about missing variables/parameters that are automatically replaced by zero values. |
Either a prediction function f(..., attach.input = FALSE) where the
variables/parameters are passed as named arguments or a prediction function
f(M, p, attach.input = FALSE) where M is the matrix of variable values
(colums with colnames correspond to different variables) and p is the vector of
parameter values.
The argument attach.input determines whether M is attached to the output.
The function f returns a matrix.
library(ggplot2) myfun <- funC0(c(y = "a*x^4 + b*x^2 + c")) out <- myfun(a = -1, b = 2, c = 3, x = seq(-2, 2, .1), attach.input = TRUE) qplot(x = x, y = y, data = as.data.frame(out), geom = "line")library(ggplot2) myfun <- funC0(c(y = "a*x^4 + b*x^2 + c")) out <- myfun(a = -1, b = 2, c = 3, x = seq(-2, 2, .1), attach.input = TRUE) qplot(x = x, y = y, data = as.data.frame(out), geom = "line")
Get coefficients from a character
getCoefficients(char, symbol)getCoefficients(char, symbol)
char |
character, e.g. "2*x + y" |
symbol |
single character, e.g. "x" or "y" |
numeric vector with the coefficients
Extract the conditions of an object
getConditions(x, ...) ## S3 method for class 'list' getConditions(x, ...) ## S3 method for class 'fn' getConditions(x, ...) ## S3 method for class 'tbl_df' getConditions(x, hypothesis = 1, ...)getConditions(x, ...) ## S3 method for class 'list' getConditions(x, ...) ## S3 method for class 'fn' getConditions(x, ...) ## S3 method for class 'tbl_df' getConditions(x, hypothesis = 1, ...)
x |
object from which the conditions should be extracted |
... |
additional arguments (not used right now) |
hypothesis |
The hypothesis in the dMod.frame |
The conditions in a format that depends on the class of x.
Extract the derivatives of an object
getDerivs(x, ...) ## S3 method for class 'parvec' getDerivs(x, ...) ## S3 method for class 'prdframe' getDerivs(x, ...) ## S3 method for class 'prdlist' getDerivs(x, ...) ## S3 method for class 'list' getDerivs(x, ...) ## S3 method for class 'objlist' getDerivs(x, ...)getDerivs(x, ...) ## S3 method for class 'parvec' getDerivs(x, ...) ## S3 method for class 'prdframe' getDerivs(x, ...) ## S3 method for class 'prdlist' getDerivs(x, ...) ## S3 method for class 'list' getDerivs(x, ...) ## S3 method for class 'objlist' getDerivs(x, ...)
x |
object from which the derivatives should be extracted |
... |
additional arguments (not used right now) |
The derivatives in a format that depends on the class of x.
This is
parvec -> matrix,
prdframe -> prdframe,
prdlist -> prdlist,
objlist -> named numeric.
Get Symbols and Numeric constants from a character
getElements(char, exclude = NULL)getElements(char, exclude = NULL)
char |
Character vector (e.g. equation) |
exclude |
Character vector, the symbols to be excluded from the return value |
getElements(c("A*AB+B^2"))getElements(c("A*AB+B^2"))
Extract the equations of an object
getEquations(x, conditions = NULL) ## S3 method for class 'odemodel' getEquations(x, conditions = NULL) ## S3 method for class 'prdfn' getEquations(x, conditions = NULL) ## S3 method for class 'fn' getEquations(x, conditions = NULL)getEquations(x, conditions = NULL) ## S3 method for class 'odemodel' getEquations(x, conditions = NULL) ## S3 method for class 'prdfn' getEquations(x, conditions = NULL) ## S3 method for class 'fn' getEquations(x, conditions = NULL)
x |
object from which the equations should be extracted |
conditions |
character or numeric vector specifying the conditions to
which |
The equations as list of eqnvec objects.
Get Parameter mappings outer.estgrid - inner.estgrid
getEstGridParameterMapping(x)getEstGridParameterMapping(x)
x |
est.grid |
named character
Daniel Lill ([email protected])
est.grid <- data.frame(ID = 1:2, condition = c("A", "B"), k1 = c("k1_A", "k1_B"), k2 = c("k2_A", "k2_B"), k3 = c("k3", NA), k4 = c("k4", "k4"), stringsAsFactors = FALSE) getEstGridParameterMapping(est.grid)est.grid <- data.frame(ID = 1:2, condition = c("A", "B"), k1 = c("k1_A", "k1_B"), k2 = c("k2_A", "k2_B"), k3 = c("k3", NA), k4 = c("k4", "k4"), stringsAsFactors = FALSE) getEstGridParameterMapping(est.grid)
Generate list of fluxes from equation list
getFluxes(eqnlist, type = c("conc", "amount"))getFluxes(eqnlist, type = c("conc", "amount"))
eqnlist |
object of class eqnlist. |
type |
"conc." or "amount" for fluxes in units of concentrations or number of molecules. |
list of named characters, the in- and out-fluxes for each state.
# Write your example here. You can also add more Start..End blocks if needed. # Please mask all output such as print() with the special tag # # such that the test is not littered. Statements guarded by are enabled # in the example file which is extracted from this test file. To extract the # example run # extractExamples() # on the R command line. ## Generate another equation list eq <- eqnlist() eq <- addReaction(eq, "A", "pA", "act_A * A * stimulus", "Phosphorylation of A") eq <- addReaction(eq, "pA", "A", "deact_A * pA", "Deposphorylation of pA") eq <- addReaction(eq, "2*pA", "pA_pA", "form_complex_pA * pA^2", "Complex formation of pA") eq <- addReaction(eq, "B", "pB", "act_B * B * pA_pA", "Phosphorylation of B") eq <- addReaction(eq, "pB", "B", "deact_B * pB", "Deposphorylation of pB") ## Extract data.frame of reactions reactions <- getReactions(eq) print(reactions) ## Get conserved quantities cq <- conservedQuantities(eq$smatrix) print(cq) ## Get fluxes fluxes <- getFluxes(eq) print(fluxes) ## Subsetting of equation list subeq1 <- subset(eq, "pB" %in% Product) print(subeq1) subeq2 <- subset(eq, grepl("not_available", Description)) print(subeq2) ## Time derivatives of observables observables <- eqnvec(pA_obs = "s1*pA", tA_obs = "s2*(A + pA)") dobs <- dot(observables, eq) ## Combined equation vector for ODE and observables f <- c(as.eqnvec(eq), dobs) print(f)# Write your example here. You can also add more Start..End blocks if needed. # Please mask all output such as print() with the special tag # # such that the test is not littered. Statements guarded by are enabled # in the example file which is extracted from this test file. To extract the # example run # extractExamples() # on the R command line. ## Generate another equation list eq <- eqnlist() eq <- addReaction(eq, "A", "pA", "act_A * A * stimulus", "Phosphorylation of A") eq <- addReaction(eq, "pA", "A", "deact_A * pA", "Deposphorylation of pA") eq <- addReaction(eq, "2*pA", "pA_pA", "form_complex_pA * pA^2", "Complex formation of pA") eq <- addReaction(eq, "B", "pB", "act_B * B * pA_pA", "Phosphorylation of B") eq <- addReaction(eq, "pB", "B", "deact_B * pB", "Deposphorylation of pB") ## Extract data.frame of reactions reactions <- getReactions(eq) print(reactions) ## Get conserved quantities cq <- conservedQuantities(eq$smatrix) print(cq) ## Get fluxes fluxes <- getFluxes(eq) print(fluxes) ## Subsetting of equation list subeq1 <- subset(eq, "pB" %in% Product) print(subeq1) subeq2 <- subset(eq, grepl("not_available", Description)) print(subeq2) ## Time derivatives of observables observables <- eqnvec(pA_obs = "s1*pA", tA_obs = "s2*(A + pA)") dobs <- dot(observables, eq) ## Combined equation vector for ODE and observables f <- c(as.eqnvec(eq), dobs) print(f)
Determine loaded DLLs available in working directory
getLocalDLLs()getLocalDLLs()
Character vector with the names of the loaded DLLs available in the working directory
Extract the observables of an object
getObservables(x, ...) ## S3 method for class 'tbl_df' getObservables(x, hypothesis = 1, ...)getObservables(x, ...) ## S3 method for class 'tbl_df' getObservables(x, hypothesis = 1, ...)
x |
object from which the equations should be extracted |
... |
not used |
hypothesis |
The hypothesis in the dMod.frame |
The equations as a character.
Extract the parameters of an object
getParameters(..., conditions = NULL) ## S3 method for class 'odemodel' getParameters(x, conditions = NULL) ## S3 method for class 'fn' getParameters(x, conditions = NULL) ## S3 method for class 'parvec' getParameters(x, conditions = NULL) ## S3 method for class 'prdframe' getParameters(x, conditions = NULL) ## S3 method for class 'prdlist' getParameters(x, conditions = NULL) ## S3 method for class 'eqnlist' getParameters(x) ## S3 method for class 'eventlist' getParameters(x)getParameters(..., conditions = NULL) ## S3 method for class 'odemodel' getParameters(x, conditions = NULL) ## S3 method for class 'fn' getParameters(x, conditions = NULL) ## S3 method for class 'parvec' getParameters(x, conditions = NULL) ## S3 method for class 'prdframe' getParameters(x, conditions = NULL) ## S3 method for class 'prdlist' getParameters(x, conditions = NULL) ## S3 method for class 'eqnlist' getParameters(x) ## S3 method for class 'eventlist' getParameters(x)
... |
objects from which the parameters should be extracted |
conditions |
character vector specifying the conditions to
which |
x |
object from which the parameters are extracted |
The parameters in a format that depends on the class of x.
Title
## S3 method for class 'data.table' getParameters(x, ...)## S3 method for class 'data.table' getParameters(x, ...)
x |
|
... |
Daniel Lill ([email protected])
Determine which parameters need sensitivities
getParametersToEstimate(est.grid, trafo, reactions)getParametersToEstimate(est.grid, trafo, reactions)
est.grid |
est.grid = data.table(condition, par1,...,parN) |
trafo |
symbolic base trafo |
reactions |
Object of class |
character, to go into the estimate-argument of odemodel
Daniel Lill ([email protected])
Get est.grid and fixed.grid
getParGrids(mytrafo, mytrafoL, mycondition.grid, SS_pars = NULL)getParGrids(mytrafo, mytrafoL, mycondition.grid, SS_pars = NULL)
mytrafo |
base trafo |
mytrafoL |
condition specific branched trafo list |
mycondition.grid |
condition.grid with condition names as rownames e.g. as output from attr(datalist, "condition.grid") |
SS_pars |
parameters determined by the steady state |
Svenja Kemmer
Generate a table of reactions (data.frame) from an equation list
getReactions(eqnlist)getReactions(eqnlist)
eqnlist |
object of class eqnlist |
data.frame with educts, products, rate and description. The first
column is a check if the reactions comply with reaction kinetics.
# Write your example here. You can also add more Start..End blocks if needed. # Please mask all output such as print() with the special tag # # such that the test is not littered. Statements guarded by are enabled # in the example file which is extracted from this test file. To extract the # example run # extractExamples() # on the R command line. ## Generate another equation list eq <- eqnlist() eq <- addReaction(eq, "A", "pA", "act_A * A * stimulus", "Phosphorylation of A") eq <- addReaction(eq, "pA", "A", "deact_A * pA", "Deposphorylation of pA") eq <- addReaction(eq, "2*pA", "pA_pA", "form_complex_pA * pA^2", "Complex formation of pA") eq <- addReaction(eq, "B", "pB", "act_B * B * pA_pA", "Phosphorylation of B") eq <- addReaction(eq, "pB", "B", "deact_B * pB", "Deposphorylation of pB") ## Extract data.frame of reactions reactions <- getReactions(eq) print(reactions) ## Get conserved quantities cq <- conservedQuantities(eq$smatrix) print(cq) ## Get fluxes fluxes <- getFluxes(eq) print(fluxes) ## Subsetting of equation list subeq1 <- subset(eq, "pB" %in% Product) print(subeq1) subeq2 <- subset(eq, grepl("not_available", Description)) print(subeq2) ## Time derivatives of observables observables <- eqnvec(pA_obs = "s1*pA", tA_obs = "s2*(A + pA)") dobs <- dot(observables, eq) ## Combined equation vector for ODE and observables f <- c(as.eqnvec(eq), dobs) print(f)# Write your example here. You can also add more Start..End blocks if needed. # Please mask all output such as print() with the special tag # # such that the test is not littered. Statements guarded by are enabled # in the example file which is extracted from this test file. To extract the # example run # extractExamples() # on the R command line. ## Generate another equation list eq <- eqnlist() eq <- addReaction(eq, "A", "pA", "act_A * A * stimulus", "Phosphorylation of A") eq <- addReaction(eq, "pA", "A", "deact_A * pA", "Deposphorylation of pA") eq <- addReaction(eq, "2*pA", "pA_pA", "form_complex_pA * pA^2", "Complex formation of pA") eq <- addReaction(eq, "B", "pB", "act_B * B * pA_pA", "Phosphorylation of B") eq <- addReaction(eq, "pB", "B", "deact_B * pB", "Deposphorylation of pB") ## Extract data.frame of reactions reactions <- getReactions(eq) print(reactions) ## Get conserved quantities cq <- conservedQuantities(eq$smatrix) print(cq) ## Get fluxes fluxes <- getFluxes(eq) print(fluxes) ## Subsetting of equation list subeq1 <- subset(eq, "pB" %in% Product) print(subeq1) subeq2 <- subset(eq, grepl("not_available", Description)) print(subeq2) ## Time derivatives of observables observables <- eqnvec(pA_obs = "s1*pA", tA_obs = "s2*(A + pA)") dobs <- dot(observables, eq) ## Combined equation vector for ODE and observables f <- c(as.eqnvec(eq), dobs) print(f)
Get the indices of the n largest (not necessarily best) steps of a parframe
getStepIndices(myparframe, nsteps = 5, tol = 1)getStepIndices(myparframe, nsteps = 5, tol = 1)
myparframe |
parframe, result from mstrust |
nsteps |
number of steps |
tol |
tolerance for stepdetection |
indices of the largest steps
## Generate a prediction function regfn <- c(y = "sin(a*time)") g <- Y(regfn, parameters = "a") x <- Xt(condition = "C1") ## Generate data data <- datalist( C1 = data.frame( name = "y", time = 1:5, value = sin(1:5) + rnorm(5, 0, .1), sigma = .1 ) ) ## Initialize parameters and time pars <- c(a = 1) times <- seq(0, 5, .1) ## Do many fits from random positions and store them into parlist out <- as.parlist(lapply(1:50, function(i) { trust(normL2(data, g*x), pars + rnorm(length(pars), 0, 1), rinit = 1, rmax = 10) })) ## Reduce parlist to parframe parframe <- as.parframe(out) plotValues(parframe) ## Get steps getStepIndices(parframe, nsteps = 2, tol = 1) getSteps(parframe, nsteps = 2) plotValues(getSteps(parframe, nsteps = 2))## Generate a prediction function regfn <- c(y = "sin(a*time)") g <- Y(regfn, parameters = "a") x <- Xt(condition = "C1") ## Generate data data <- datalist( C1 = data.frame( name = "y", time = 1:5, value = sin(1:5) + rnorm(5, 0, .1), sigma = .1 ) ) ## Initialize parameters and time pars <- c(a = 1) times <- seq(0, 5, .1) ## Do many fits from random positions and store them into parlist out <- as.parlist(lapply(1:50, function(i) { trust(normL2(data, g*x), pars + rnorm(length(pars), 0, 1), rinit = 1, rmax = 10) })) ## Reduce parlist to parframe parframe <- as.parframe(out) plotValues(parframe) ## Get steps getStepIndices(parframe, nsteps = 2, tol = 1) getSteps(parframe, nsteps = 2) plotValues(getSteps(parframe, nsteps = 2))
Get the rows of the n largest steps of a parframe
getSteps(myparframe, nsteps = 5, tol = 1)getSteps(myparframe, nsteps = 5, tol = 1)
myparframe |
parframe, result from mstrust |
nsteps |
number of steps |
tol |
tolerance for stepdetection |
parframe subsetted to the n largest steps
## Generate a prediction function regfn <- c(y = "sin(a*time)") g <- Y(regfn, parameters = "a") x <- Xt(condition = "C1") ## Generate data data <- datalist( C1 = data.frame( name = "y", time = 1:5, value = sin(1:5) + rnorm(5, 0, .1), sigma = .1 ) ) ## Initialize parameters and time pars <- c(a = 1) times <- seq(0, 5, .1) ## Do many fits from random positions and store them into parlist out <- as.parlist(lapply(1:50, function(i) { trust(normL2(data, g*x), pars + rnorm(length(pars), 0, 1), rinit = 1, rmax = 10) })) ## Reduce parlist to parframe parframe <- as.parframe(out) plotValues(parframe) ## Get steps getStepIndices(parframe, nsteps = 2, tol = 1) getSteps(parframe, nsteps = 2) plotValues(getSteps(parframe, nsteps = 2))## Generate a prediction function regfn <- c(y = "sin(a*time)") g <- Y(regfn, parameters = "a") x <- Xt(condition = "C1") ## Generate data data <- datalist( C1 = data.frame( name = "y", time = 1:5, value = sin(1:5) + rnorm(5, 0, .1), sigma = .1 ) ) ## Initialize parameters and time pars <- c(a = 1) times <- seq(0, 5, .1) ## Do many fits from random positions and store them into parlist out <- as.parlist(lapply(1:50, function(i) { trust(normL2(data, g*x), pars + rnorm(length(pars), 0, 1), rinit = 1, rmax = 10) })) ## Reduce parlist to parframe parframe <- as.parframe(out) plotValues(parframe) ## Get steps getStepIndices(parframe, nsteps = 2, tol = 1) getSteps(parframe, nsteps = 2) plotValues(getSteps(parframe, nsteps = 2))
Convenience function to show last plot in an external viewer.
ggopen(plot = last_plot(), command = "xdg-open", ...)ggopen(plot = last_plot(), command = "xdg-open", ...)
plot |
|
command |
character, indicatig which pdf viewer is started. |
... |
arguments going to |
Ensure all tables are data.tables
gridlist(est.grid, fix.grid)gridlist(est.grid, fix.grid)
est.grid, fix.grid
|
data.tables, will be coerced to one |
list of the grid
An identity function which vanishes upon concatenation of fns
Id()Id()
fn of class idfn
x <- Xt() id <- Id() (id*x)(1:10, pars = c(a = 1)) (x*id)(1:10, pars = c(a = 1)) str(id*x) str(x*id)x <- Xt() id <- Id() (id*x)(1:10, pars = c(a = 1)) (x*id)(1:10, pars = c(a = 1)) str(id*x) str(x*id)
Requires AMICI https://github.com/ICB-DCM/AMICI/ and dependencies to be installed on your system Big thanks go to Daniel Weindl!
import_sbml(modelpath, amicipath = NULL)import_sbml(modelpath, amicipath = NULL)
modelpath |
Path to the sbml file |
amicipath |
Path to your amici-python-installation, e.g.: AMICIPATH/python |
list of eqnlist, parameters and inits
Add single-valued parameters to the pargrids
indiv_addGlobalParsToGridlist(pars, gridlist, FLAGoverwrite = FALSE)indiv_addGlobalParsToGridlist(pars, gridlist, FLAGoverwrite = FALSE)
pars |
character or numeric |
gridlist |
list(fix.grid, est.grid) |
FLAGoverwrite |
Overwrite existing parameters in fix.grid or est.grid? |
gridlist
Daniel Lill ([email protected])
pars <- c(NewParSymbolic = "NewParSymbolic", NewParFixed = 1) est.grid <- data.table(ID = 1:2, condition = c("A", "B"), k1 = c("k1_A", "k1_B"), k2 = c("k2_A", "k2_B"), k3 = c("k3", NA), k4 = c("k4", "k4"), stringsAsFactors = FALSE) fix.grid <- data.table(ID = 1:2, condition = c("A", "B"), k3 = c(NA, 3.5), k5 = c(5.1,5.2), k6 = c(6,6), stringsAsFactors = FALSE) indiv_addGlobalParsToGridlist(c(NewParSymbolic = "NewParSymbolic", NewParFixed = 1), gridlist(est.grid = est.grid, fix.grid = fix.grid)) indiv_addGlobalParsToGridlist(c(k1 = 1), gridlist(est.grid = est.grid, fix.grid = fix.grid), FLAGoverwrite = FALSE) # nothing happens indiv_addGlobalParsToGridlist(c(k1 = 1), gridlist(est.grid = est.grid, fix.grid = fix.grid), FLAGoverwrite = TRUE) # k1 is replaced and moved to fix.gridpars <- c(NewParSymbolic = "NewParSymbolic", NewParFixed = 1) est.grid <- data.table(ID = 1:2, condition = c("A", "B"), k1 = c("k1_A", "k1_B"), k2 = c("k2_A", "k2_B"), k3 = c("k3", NA), k4 = c("k4", "k4"), stringsAsFactors = FALSE) fix.grid <- data.table(ID = 1:2, condition = c("A", "B"), k3 = c(NA, 3.5), k5 = c(5.1,5.2), k6 = c(6,6), stringsAsFactors = FALSE) indiv_addGlobalParsToGridlist(c(NewParSymbolic = "NewParSymbolic", NewParFixed = 1), gridlist(est.grid = est.grid, fix.grid = fix.grid)) indiv_addGlobalParsToGridlist(c(k1 = 1), gridlist(est.grid = est.grid, fix.grid = fix.grid), FLAGoverwrite = FALSE) # nothing happens indiv_addGlobalParsToGridlist(c(k1 = 1), gridlist(est.grid = est.grid, fix.grid = fix.grid), FLAGoverwrite = TRUE) # k1 is replaced and moved to fix.grid
Can only add parameters, cannot update existing parameters
indiv_addLocalParsToGridList(pars, gridlist, FLAGoverwrite = FALSE)indiv_addLocalParsToGridList(pars, gridlist, FLAGoverwrite = FALSE)
pars |
data.table(ID, condition, pars...) Only one of ID, condition must be present |
gridlist |
modified gridlist
Daniel Lill ([email protected])
pars <- data.table(ID = 1:2, newFix = c(1,2), newEst = c("a", "b"), newMix = c("a", 1)) est.grid <- data.table(ID = 1:2, condition = c("A", "B"), k1 = c("k1_A", "k1_B"), k2 = c("k2_A", "k2_B"), k3 = c("k3", NA), k4 = c("k4", "k4"), stringsAsFactors = FALSE) fix.grid <- data.table(ID = 1:2, condition = c("A", "B"), k3 = c(NA, 3.5), k5 = c(5.1,5.2), k6 = c(6,6), stringsAsFactors = FALSE) gl <- gridlist(est.grid = est.grid, fix.grid = fix.grid) indiv_addLocalParsToGridList(pars, gl)pars <- data.table(ID = 1:2, newFix = c(1,2), newEst = c("a", "b"), newMix = c("a", 1)) est.grid <- data.table(ID = 1:2, condition = c("A", "B"), k1 = c("k1_A", "k1_B"), k2 = c("k2_A", "k2_B"), k3 = c("k3", NA), k4 = c("k4", "k4"), stringsAsFactors = FALSE) fix.grid <- data.table(ID = 1:2, condition = c("A", "B"), k3 = c(NA, 3.5), k5 = c(5.1,5.2), k6 = c(6,6), stringsAsFactors = FALSE) gl <- gridlist(est.grid = est.grid, fix.grid = fix.grid) indiv_addLocalParsToGridList(pars, gl)
Create an objlist with zeros as entries
init_empty_objlist(pars, deriv = TRUE, FLAGchisquare = FALSE)init_empty_objlist(pars, deriv = TRUE, FLAGchisquare = FALSE)
pars |
named vector. Only names and length are used |
deriv |
TRUE or FALSE |
Daniel Lill ([email protected])
init_empty_objlist(setNames(rnorm(5), letters[1:5])) init_empty_objlist(setNames(rnorm(5), letters[1:5]), FLAGchisquare = TRUE)init_empty_objlist(setNames(rnorm(5), letters[1:5])) init_empty_objlist(setNames(rnorm(5), letters[1:5]), FLAGchisquare = TRUE)
Installs Julia and the necessary julia packages
installJuliaForSteadyStates(installJulia = FALSE, installJuliaPackages = TRUE)installJuliaForSteadyStates(installJulia = FALSE, installJuliaPackages = TRUE)
installJulia |
boolean, default |
installJuliaPackages |
boolean, default |
nothing
Phosphorylated Epo receptor (pEpoR), phosphorylated STAT in the cytoplasm (tpSTAT) and total STAT (tSTAT) in the cytoplasmhave been measured at times 0, ..., 60.
Bind named list of data.frames into one data.frame
lbind(mylist)lbind(mylist)
mylist |
A named list of data.frame. The data.frames are expected to have the same structure. |
Each data.frame ist augented by a "condition" column containing the name attributed of
the list entry. Subsequently, the augmented data.frames are bound together by rbind.
data.frame with the originial columns augmented by a "condition" column.
An aborted mstrust
leaves behind results of already completed fits. This command loads these
fits into a fitlist.
load.parlist(folder)load.parlist(folder)
folder |
Path to the folder where the fit has left its results. |
The command mstrust saves
each completed fit along the multi-start sequence such that the results can
be resurected on abortion. This command loads a fitlist from these
intermediate results.
An object of class parlist.
Wolfgang Mader, [email protected]
Usually when restarting the R session, although all objects are saved in
the workspace, the dynamic libraries are not linked any more. loadDLL
is a wrapper for dyn.load that uses the "modelname" attribute of
dMod objects like prediction functions, observation functions, etc. to
load the corresponding shared object.
loadDLL(...)loadDLL(...)
... |
objects of class prdfn, obsfn, parfn, objfn, ... |
Translate long to wide format (inverse of wide2long.matrix)
long2wide(out)long2wide(out)
out |
data.frame in long format |
data.frame in wide format
Lists the objects for a set of classes.
lsdMod( classlist = c("odemodel", "parfn", "prdfn", "obsfn", "objfn", "datalist"), envir = .GlobalEnv )lsdMod( classlist = c("odemodel", "parfn", "prdfn", "obsfn", "objfn", "datalist"), envir = .GlobalEnv )
classlist |
List of object classes to print. |
envir |
Alternative environment to search for objects. |
## Not run: lsdMod() lsdMod(classlist = "prdfn", envir = environment(obj)) ## End(Not run)## Not run: lsdMod() lsdMod(classlist = "prdfn", envir = environment(obj)) ## End(Not run)
Extract pars, fixed and parnames from grids for a given condition
make_pars(pars, fixed = NULL, est.grid, fix.grid, ID)make_pars(pars, fixed = NULL, est.grid, fix.grid, ID)
est.grid |
needs condition and ID |
fix.grid |
|
ID |
|
est.vec |
pars <- c("k1_A" = 1, "k1_B" = 1, "k2_A" = 2, "k3" = 3, "k4" = 4) fixed <- c("k2_B" = 2.5) est.grid <- data.frame(ID = 1:2, condition = c("A", "B"), k1 = c("k1_A", "k1_B"), k2 = c("k2_A", "k2_B"), k3 = c("k3", NA), k4 = c("k4", "k4"), stringsAsFactors = FALSE) fix.grid <- data.frame(ID = 1:2, condition = c("A", "B"), k3 = c(NA, 3.5), k5 = c(5.1,5.2), k6 = c(6,6), stringsAsFactors = FALSE) make_pars(pars, fixed, est.grid, fix.grid, 1) make_pars(pars, fixed, est.grid, fix.grid, 2)pars <- c("k1_A" = 1, "k1_B" = 1, "k2_A" = 2, "k3" = 3, "k4" = 4) fixed <- c("k2_B" = 2.5) est.grid <- data.frame(ID = 1:2, condition = c("A", "B"), k1 = c("k1_A", "k1_B"), k2 = c("k2_A", "k2_B"), k3 = c("k3", NA), k4 = c("k4", "k4"), stringsAsFactors = FALSE) fix.grid <- data.frame(ID = 1:2, condition = c("A", "B"), k3 = c(NA, 3.5), k5 = c(5.1,5.2), k6 = c(6,6), stringsAsFactors = FALSE) make_pars(pars, fixed, est.grid, fix.grid, 1) make_pars(pars, fixed, est.grid, fix.grid, 2)
The function is exported for dependency reasons
match.fnargs(arglist, choices)match.fnargs(arglist, choices)
arglist |
list |
choices |
character |
Get modelname from single object (used internally)
mname(x, conditions = NULL) ## S3 method for class ''NULL'' mname(x, conditions = NULL) ## S3 method for class 'character' mname(x, conditions = NULL) ## S3 method for class 'objfn' mname(x, conditions = NULL) ## S3 method for class 'fn' mname(x, conditions = NULL)mname(x, conditions = NULL) ## S3 method for class ''NULL'' mname(x, conditions = NULL) ## S3 method for class 'character' mname(x, conditions = NULL) ## S3 method for class 'objfn' mname(x, conditions = NULL) ## S3 method for class 'fn' mname(x, conditions = NULL)
x |
dMod object |
conditions |
character vector of conditions |
The modelname attribute refers to the name of a C file associated with a dMod function object like prediction-, parameter transformation- or objective functions.
modelname(..., conditions = NULL) modelname(x, ...) <- value ## S3 replacement method for class 'fn' modelname(x, conditions = NULL, ...) <- value ## S3 replacement method for class 'objfn' modelname(x, conditions = NULL, ...) <- valuemodelname(..., conditions = NULL) modelname(x, ...) <- value ## S3 replacement method for class 'fn' modelname(x, conditions = NULL, ...) <- value ## S3 replacement method for class 'objfn' modelname(x, conditions = NULL, ...) <- value
... |
objects of type |
conditions |
character vector of conditions |
x |
dMod object for which the model name should be set |
value |
character, the new modelname (does not change the C file) |
character vector of model names, corresponding to C files in the local directory.
The output of this function can be used for the center - argument of mstrust
msParframe( pars, n = 20, seed = 12345, samplefun = stats::rnorm, keepfirst = TRUE, ... )msParframe( pars, n = 20, seed = 12345, samplefun = stats::rnorm, keepfirst = TRUE, ... )
pars |
Named vector. If |
n |
Integer how many lines should the parframe have |
seed |
Seed for the random number generator |
samplefun |
|
keepfirst |
boolean, if set to |
... |
arguments going to samplefun |
parframe (without metanames)
msParframe(c(a = 0, b = 100000), 5) # Parameter specific sigma msParframe(c(a = 0, b = 100000), 5, samplefun = rnorm, sd = c(100, 0.5))msParframe(c(a = 0, b = 100000), 5) # Parameter specific sigma msParframe(c(a = 0, b = 100000), 5, samplefun = rnorm, sd = c(100, 0.5))
Wrapper around trust allowing for multiple fits
from randomly chosen initial values.
mstrust( objfun, center, studyname, rinit = 0.1, rmax = 10, fits = 20, cores = 1, optmethod = "trust", samplefun = "rnorm", resultPath = ".", stats = FALSE, output = FALSE, cautiousMode = FALSE, start1stfromCenter = FALSE, ... )mstrust( objfun, center, studyname, rinit = 0.1, rmax = 10, fits = 20, cores = 1, optmethod = "trust", samplefun = "rnorm", resultPath = ".", stats = FALSE, output = FALSE, cautiousMode = FALSE, start1stfromCenter = FALSE, ... )
objfun |
Objective function, see |
center |
Parameter values around which the initial values for each fit
are randomly sampled. The initial values handed to trust are the sum
of center and the output of samplefun, center +
samplefun. See |
studyname |
The names of the study or fit. This name is used to determine filenames for interim and final results. See Details. |
rinit |
Starting trust region radius, see |
rmax |
Maximum allowed trust region radius, see |
fits |
Number of fits (jobs). |
cores |
Number of cores for job parallelization. |
samplefun |
Function to sample random initial values. It is assumed,
that samplefun has a named parameter "n" which defines how many
random numbers are to be returned, such as for |
resultPath |
character indicating the folder where the results should be stored. Defaults to ".". |
stats |
If true, the same summary statistic as written to the logfile is printed to command line on mstrust completion. |
output |
logical. If true, writes output to the disc. |
cautiousMode |
Logical, write every fit to disk in deparsed form (avoids the RDA incompatibility trap) and don't delete intermediate results |
... |
Additional parameters handed to trust(), samplefun(), or the objective function by matching parameter names. All unmatched parameters are handed to the objective function objfun(). The log file starts with a table telling which parameter was assigend to which function. |
By running multiple fits starting at randomly chosen inital
parameters, the chisquare landscape can be explored using a deterministic
optimizer. Here, trust is used for optimization. The standard
procedure to obtain random initial values is to sample random variables
from a uniform distribution (rnorm) and adding these to
center. It is, however, possible, to employ any other sampling
strategy by handing the respective function to mstrust(),
samplefun.
In case a special sampling is required, a customized sampling function can be used. If, e.g., inital values leading to a non-physical systems are to be discarded upfront, the objective function can be addapted accordingly.
All started fits either lead to an error or complete converged or unconverged. A statistics about the return status of fits can be shown by setting stats to TRUE.
Fit final and intermediat results are stored under studyname. For
each run of mstrust for the same study name, a folder under
studyname of the form "trial-x-date" is created. "x" is the number
of the trial, date is the current time stamp. In this folder, the
intermediate results are stored. These intermediate results can be loaded
by load.parlist. These are removed on successfull completion
of mstrust. In this case, the final list of fit parameters
(parameterList.Rda) and the fit log (mstrust.log) are found instead.
A parlist holding errored and converged fits.
Wolfgang Mader, [email protected]
1. trust, for the used optimizer,
2. rnorm, runif for two common sampling functions,
3. msParframe for passing a reproducible set of random initial
guesses to mstrust,
4. as.parframe for formatting the output to a handy table
This is basically dplyr::mutate with two differences
1. It groups the tibble rowwise before mutating
2. It allows to store the calls. This is intended for use when your objects are
not standard transformations, such as prd = g*x*p but more complicated and
you want to keep a record of what you did.
If the result of your ... expressions is not atomic, make sure to wrap your
expression in list()
mutatedMod.frame(dMod.frame, ..., keepCalls = F)mutatedMod.frame(dMod.frame, ..., keepCalls = F)
dMod.frame |
A dMod.frame |
... |
Expressions going to mutate() |
keepCalls |
Should the dots ... be recorded? |
## Not run: mytbl <- tibble(a = 1:2, b = letters[1:2]) %>% mutatedMod.frame(e = paste0(a,b), keepCalls = T) %>% mutatedMod.frame(f = paste0(e, "=", a, "*", b), keepCalls = T) %>% mutatedMod.frame(e = paste0(a,"*", b), keepCalls = T) mytbl$.calls ## End(Not run)## Not run: mytbl <- tibble(a = 1:2, b = letters[1:2]) %>% mutatedMod.frame(e = paste0(a,b), keepCalls = T) %>% mutatedMod.frame(f = paste0(e, "=", a, "*", b), keepCalls = T) %>% mutatedMod.frame(e = paste0(a,"*", b), keepCalls = T) mytbl$.calls ## End(Not run)
Gaussian Log-likelihood. Supports NONMEM-like BLOQ handling methods M1, M3 and M4 and estimation of error models
nll( nout, pars, deriv, opt.BLOQ = "M3", opt.hessian = c(ALOQ_part1 = TRUE, ALOQ_part2 = TRUE, ALOQ_part3 = TRUE, BLOQ_part1 = TRUE, BLOQ_part2 = TRUE, BLOQ_part3 = TRUE, PD = TRUE) )nll( nout, pars, deriv, opt.BLOQ = "M3", opt.hessian = c(ALOQ_part1 = TRUE, ALOQ_part2 = TRUE, ALOQ_part3 = TRUE, BLOQ_part1 = TRUE, BLOQ_part2 = TRUE, BLOQ_part3 = TRUE, PD = TRUE) )
nout |
|
pars |
Example named vector of outer parameters to construct the objlist |
deriv |
TRUE or FALSE |
opt.BLOQ |
see normIndiv |
opt.hessian |
see normIndiv |
list with entries value (numeric, the weighted residual sum of squares), gradient (numeric, gradient) and hessian (matrix of type numeric).
Non-linear log likelihood for the ALOQ part of the data
nll_ALOQ( nout, derivs, derivs.err, opt.BLOQ = c("M3", "M4NM", "M4BEAL", "M1"), opt.hessian = c(ALOQ_part1 = TRUE, ALOQ_part2 = TRUE, ALOQ_part3 = TRUE) )nll_ALOQ( nout, derivs, derivs.err, opt.BLOQ = c("M3", "M4NM", "M4BEAL", "M1"), opt.hessian = c(ALOQ_part1 = TRUE, ALOQ_part2 = TRUE, ALOQ_part3 = TRUE) )
nout |
output of |
derivs, derivs.err
|
attributes of output of |
opt.BLOQ |
Character denoting the method to deal with BLOQ data |
opt.hessian |
Named logical vector to include or exclude various non-convex summands of the hessian matrix |
Non-linear log likelihood for the BLOQ part of the data
nll_BLOQ( nout.bloq, derivs.bloq, derivs.err.bloq, opt.BLOQ = c("M3", "M4NM", "M4BEAL", "M1"), opt.hessian = c(BLOQ_part1 = TRUE, BLOQ_part2 = TRUE, BLOQ_part3 = TRUE) )nll_BLOQ( nout.bloq, derivs.bloq, derivs.err.bloq, opt.BLOQ = c("M3", "M4NM", "M4BEAL", "M1"), opt.hessian = c(BLOQ_part1 = TRUE, BLOQ_part2 = TRUE, BLOQ_part3 = TRUE) )
nout.bloq |
The bloq output of |
derivs.bloq, derivs.err.bloq
|
attributes of output of |
opt.BLOQ |
Character denoting the method to deal with BLOQ data |
opt.hessian |
Named logical vector to include or exclude various summands of the hessian matrix |
For parameter estimation and optimization, an objective function
is needed. normL2 returns an objective function for the L2 norm of
data and model prediction. The resulting objective function can be used for
optimization with the trust optimizer, see mstrust.
normL2(data, x, errmodel = NULL, times = NULL, attr.name = "data")normL2(data, x, errmodel = NULL, times = NULL, attr.name = "data")
data |
object of class datalist |
x |
object of class prdfn |
errmodel |
object of class obsfn. |
times |
numeric vector, additional time points where the prediction function is
evaluated. If NULL, time points are extacted from the datalist solely. If the prediction
function makes use of events, hand over event |
attr.name |
character. The constraint value is additionally returned in an attributed with this name |
forcings |
TO BE FILLED BY DANIEL K |
iiv |
Example: c("ka", "ETA_EMAX"). |
conditional |
Example: data.frame(parname = "GR", covname = "SEX", covvalue = "1", stringsAsFactors = FALSE). |
fixed.grid |
data.frame(parname, partask, ids...) Lookup table for fixed parameters |
nauxtimes |
additional simulation times |
cores |
to parallelize over conditions not over fits |
Objective functions can be combined by the "+" operator, see sumobjfn.
Object of class obsfn, i.e. a function
obj(..., fixed, deriv, conditions, env) that returns an objective list,
objlist.
## Generate a prediction function times <- 0:5 grid <- data.frame(name = "A", time = times, row.names = paste0("p", times)) x <- Xd(grid, condition = "C1") pars <- structure(rep(0, nrow(grid)), names = row.names(grid)) ## Simulate data data.list <- lapply(1:3, function(i) { prediction <- x(times, pars + rnorm(length(pars), 0, 1)) cbind(wide2long(prediction), sigma = 1) }) data <- as.datalist(do.call(rbind, data.list)) ## Generate objective function based on data and model ## Then fit the data and plot the result obj <- normL2(data, x) myfit <- trust(obj, pars, rinit = 1, rmax = 10) plot(x(times, myfit$argument), data)## Generate a prediction function times <- 0:5 grid <- data.frame(name = "A", time = times, row.names = paste0("p", times)) x <- Xd(grid, condition = "C1") pars <- structure(rep(0, nrow(grid)), names = row.names(grid)) ## Simulate data data.list <- lapply(1:3, function(i) { prediction <- x(times, pars + rnorm(length(pars), 0, 1)) cbind(wide2long(prediction), sigma = 1) }) data <- as.datalist(do.call(rbind, data.list)) ## Generate objective function based on data and model ## Then fit the data and plot the result obj <- normL2(data, x) myfit <- trust(obj, pars, rinit = 1, rmax = 10) plot(x(times, myfit$argument), data)
Fast normL2
normL2_indiv( data, prd0, errmodel = NULL, est.grid, fix.grid, times = NULL, attr.name = "data", mycores = 1 )normL2_indiv( data, prd0, errmodel = NULL, est.grid, fix.grid, times = NULL, attr.name = "data", mycores = 1 )
data |
|
prd0 |
|
errmodel |
|
est.grid |
|
fix.grid |
|
times |
|
attr.name |
|
mycores |
You can use mycores > 1 for parallel computation of conditions to speed up single fits or profiles. Use carefully when also performing multi-start optimization via mstrust. |
objective function
Daniel Lill ([email protected])
Find integer-null space of matrix A
nullZ(A, tol = sqrt(.Machine$double.eps))nullZ(A, tol = sqrt(.Machine$double.eps))
A |
matrix for which the null space is searched |
tol |
tolerance to find pivots in rref-function below |
null space of A with only integers in it
Malenka Mader, [email protected]
An objective frame is supposed to store the residuals of a model prediction with respect to a data frame.
objframe(mydata, deriv = NULL, deriv.err = NULL)objframe(mydata, deriv = NULL, deriv.err = NULL)
mydata |
data.frame as being generated by res. |
deriv |
matrix of the derivatives of the residuals with respect to parameters. |
deriv.err |
matrix of the derivatives of the error model. |
An object of class objframe, i.e. a data frame with attribute "deriv".
An objective list contains an objective value, a gradient, and a Hessian matrix.
Objective lists can contain additional numeric attributes that are preserved or combined with the corresponding attributes of another objective list when both are added by the "+" operator, see sumobjlist.
Objective lists are returned by objective functions as being generated by normL2, constraintL2, priorL2 and datapointL2.
objlist(value, gradient, hessian)objlist(value, gradient, hessian)
value |
numeric of length 1 |
gradient |
named numeric |
hessian |
matrix with rownames and colnames according to gradient names |
Object of class objlist
# objlist(1, c(a = 1, b = 2), matrix(2, nrow = 2, ncol = 2, dimnames = list(c("a", "b"),c("a", "b"))))# objlist(1, c(a = 1, b = 2), matrix(2, nrow = 2, ncol = 2, dimnames = list(c("a", "b"),c("a", "b"))))
to help the objective function integrator
objtimes(datatimes, eventtimes = NULL, Nobjtimes = 25)objtimes(datatimes, eventtimes = NULL, Nobjtimes = 25)
datatimes |
times present in data |
eventtimes |
times present in events (not yet implemented) |
vector of times, including datatimes
objtimes(c(0,30,60,90,600)) objtimes(c(30,60,90,600)) objtimes(c(-30,60,90,600))objtimes(c(0,30,60,90,600)) objtimes(c(30,60,90,600)) objtimes(c(-30,60,90,600))
An observation function is a function is that is concatenated
with a prediction function via prodfn to yield a new prediction function,
see prdfn. Observation functions are generated by Y. Handling
of the conditions is then organized by the obsfn object.
obsfn(X2Y, parameters = NULL, condition = NULL)obsfn(X2Y, parameters = NULL, condition = NULL)
X2Y |
the low-level observation function generated e.g. by Y. |
parameters |
character vector with parameter names |
condition |
character, the condition name |
Observation functions can be "added" by the "+" operator, see sumfn. Thereby, observations for different conditions are merged or, overwritten. Observation functions can also be concatenated with other functions, e.g. observation functions (obsfn) or prediction functions (prdfn) by the "*" operator, see prodfn.
Object of class obsfn, i.e. a function x(..., fixed, deriv, conditions, env)
which returns a prdlist. The arguments out (prediction) and pars (parameter values)
should be passed via the ... argument.
# Define a time grid on which to make a prediction by peace-wise linear function. # Then define a (generic) prediction function based on thid grid. times <- 0:5 grid <- data.frame(name = "A", time = times, row.names = paste0("p", times)) x <- Xd(grid) # Define an observable and an observation function observables <- eqnvec(Aobs = "s*A") g <- Y(g = observables, f = NULL, states = "A", parameters = "s") # Collect parameters and define an overarching parameter transformation # for two "experimental condtions". dynpars <- attr(x, "parameters") obspars <- attr(g, "parameters") innerpars <- c(dynpars, obspars) trafo <- structure(innerpars, names = innerpars) trafo_C1 <- replaceSymbols(innerpars, paste(innerpars, "C1", sep = "_"), trafo) trafo_C2 <- replaceSymbols(innerpars, paste(innerpars, "C2", sep = "_"), trafo) p <- NULL p <- p + P(trafo = trafo_C1, condition = "C1") p <- p + P(trafo = trafo_C2, condition = "C2") # Collect outer (overarching) parameters and # initialize with random values outerpars <- attr(p, "parameters") pars <- structure(runif(length(outerpars), 0, 1), names = outerpars) # Predict internal/unobserved states out1 <- (x*p)(times, pars) plot(out1) # Predict observed states in addition to unobserved out2 <- (g*x*p)(times, pars) plot(out2)# Define a time grid on which to make a prediction by peace-wise linear function. # Then define a (generic) prediction function based on thid grid. times <- 0:5 grid <- data.frame(name = "A", time = times, row.names = paste0("p", times)) x <- Xd(grid) # Define an observable and an observation function observables <- eqnvec(Aobs = "s*A") g <- Y(g = observables, f = NULL, states = "A", parameters = "s") # Collect parameters and define an overarching parameter transformation # for two "experimental condtions". dynpars <- attr(x, "parameters") obspars <- attr(g, "parameters") innerpars <- c(dynpars, obspars) trafo <- structure(innerpars, names = innerpars) trafo_C1 <- replaceSymbols(innerpars, paste(innerpars, "C1", sep = "_"), trafo) trafo_C2 <- replaceSymbols(innerpars, paste(innerpars, "C2", sep = "_"), trafo) p <- NULL p <- p + P(trafo = trafo_C1, condition = "C1") p <- p + P(trafo = trafo_C2, condition = "C2") # Collect outer (overarching) parameters and # initialize with random values outerpars <- attr(p, "parameters") pars <- structure(runif(length(outerpars), 0, 1), names = outerpars) # Predict internal/unobserved states out1 <- (x*p)(times, pars) plot(out1) # Predict observed states in addition to unobserved out2 <- (g*x*p)(times, pars) plot(out2)
Generate the model objects for use in Xs (models with sensitivities)
odemodel( f, deriv = TRUE, forcings = NULL, events = NULL, outputs = NULL, fixed = NULL, estimate = NULL, modelname = "odemodel", solver = c("deSolve", "Sundials"), gridpoints = NULL, verbose = FALSE, ... )odemodel( f, deriv = TRUE, forcings = NULL, events = NULL, outputs = NULL, fixed = NULL, estimate = NULL, modelname = "odemodel", solver = c("deSolve", "Sundials"), gridpoints = NULL, verbose = FALSE, ... )
f |
Something that can be converted to eqnvec, e.g. a named character vector with the ODE |
deriv |
logical, generate sensitivities or not |
forcings |
Character vector with the names of the forcings |
events |
data.frame of events with columns "var" (character, the name of the state to be
affected), "time" (character or numeric, time point), "value" (character or numeric, value),
"method" (character, either
"replace" or "add"). See events. Events need to be defined here if they contain
parameters, like the event time or value. If both, time and value are purely numeric, they
can be specified in |
outputs |
Named character vector for additional output variables. |
fixed |
Character vector with the names of parameters (initial values and dynamic) for which no sensitivities are required (will speed up the integration). |
estimate |
Character vector specifying parameters (initial values and dynamic) for which sensitivities are returned. If estimate is specified, it overwrites 'fixed'. |
modelname |
Character, the name of the C file being generated. |
solver |
Solver for which the equations are prepared. |
gridpoints |
Integer, the minimum number of time points where the ODE is evaluated internally |
verbose |
Print compiler output to R command line. |
... |
Further arguments being passed to funC. |
list with func (ODE object) and extended (ODE+Sensitivities object)
## Not run: ## Generate a compiled ODE model from an equation vector ## The model will not return sensitivities for "switch" ## Files will be generated in your working directory! f <- eqnvec(A = "-k*A + switch*F") model <- odemodel(f, forcings = "F", fixed = "switch") print(model) ## Generate the same model from an equation list f <- addReaction(NULL, from = "", to = "A", rate = "switch*F", description = "production") f <- addReaction(f , from = "A", to = "", rate = "k*A", description = "degradation") print(f) model <- odemodel(f, forcings = "F", fixed = "switch") print(model) # create forcings forc1 <- data.frame(name = "F", time = seq(0,5, 0.1), value = sin(seq(0,5,0.1))) forc2 <- data.frame(name = "F", time = seq(0,5, 0.1), value = exp(-seq(0,5,0.1))) forc3 <- data.frame(name = "F", time= 0, value = 0.1) x <- Xs(model, forc1, condition = "forc1") + Xs(model, forc2, condition = "forc2") + Xs(model, forc3, condition = "forc3") g <- Y(c(out1 = "F * A", out2 = "F"), x) times <- seq(0,5, 0.001) pars <- setNames(runif(length(getParameters(x))), getParameters(x)) pred <- (g*x)(times, pars) plot(pred) ## End(Not run)## Not run: ## Generate a compiled ODE model from an equation vector ## The model will not return sensitivities for "switch" ## Files will be generated in your working directory! f <- eqnvec(A = "-k*A + switch*F") model <- odemodel(f, forcings = "F", fixed = "switch") print(model) ## Generate the same model from an equation list f <- addReaction(NULL, from = "", to = "A", rate = "switch*F", description = "production") f <- addReaction(f , from = "A", to = "", rate = "k*A", description = "degradation") print(f) model <- odemodel(f, forcings = "F", fixed = "switch") print(model) # create forcings forc1 <- data.frame(name = "F", time = seq(0,5, 0.1), value = sin(seq(0,5,0.1))) forc2 <- data.frame(name = "F", time = seq(0,5, 0.1), value = exp(-seq(0,5,0.1))) forc3 <- data.frame(name = "F", time= 0, value = 0.1) x <- Xs(model, forc1, condition = "forc1") + Xs(model, forc2, condition = "forc2") + Xs(model, forc3, condition = "forc3") g <- Y(c(out1 = "F * A", out2 = "F"), x) times <- seq(0,5, 0.001) pars <- setNames(runif(length(getParameters(x))), getParameters(x)) pred <- (g*x)(times, pars) plot(pred) ## End(Not run)
Get default arguments for integrators
optionsLSODES( rtol = 1e-06, atol = 1e-06, jacvec = NULL, sparsetype = "sparseint", nnz = NULL, inz = NULL, rootfunc = NULL, verbose = FALSE, nroot = 0, tcrit = NULL, hmin = 0, hmax = NULL, hini = 0, ynames = TRUE, maxord = NULL, maxsteps = 5000, lrw = NULL, liw = NULL ) optionsLSODE( rtol = 1e-06, atol = 1e-06, jacfunc = NULL, jactype = "fullint", mf = NULL, rootfunc = NULL, verbose = FALSE, nroot = 0, tcrit = NULL, hmin = 0, hmax = NULL, hini = 0, ynames = TRUE, maxord = NULL, bandup = NULL, banddown = NULL, maxsteps = 5000, rpar = NULL, ipar = NULL, nout = 0, outnames = NULL ) optionsLSODA( rtol = 1e-06, atol = 1e-06, jacfunc = NULL, jactype = "fullint", rootfunc = NULL, verbose = FALSE, nroot = 0, tcrit = NULL, hmin = 0, hmax = NULL, hini = 0, ynames = TRUE, maxordn = 12, maxords = 5, bandup = NULL, banddown = NULL, maxsteps = 5000 )optionsLSODES( rtol = 1e-06, atol = 1e-06, jacvec = NULL, sparsetype = "sparseint", nnz = NULL, inz = NULL, rootfunc = NULL, verbose = FALSE, nroot = 0, tcrit = NULL, hmin = 0, hmax = NULL, hini = 0, ynames = TRUE, maxord = NULL, maxsteps = 5000, lrw = NULL, liw = NULL ) optionsLSODE( rtol = 1e-06, atol = 1e-06, jacfunc = NULL, jactype = "fullint", mf = NULL, rootfunc = NULL, verbose = FALSE, nroot = 0, tcrit = NULL, hmin = 0, hmax = NULL, hini = 0, ynames = TRUE, maxord = NULL, bandup = NULL, banddown = NULL, maxsteps = 5000, rpar = NULL, ipar = NULL, nout = 0, outnames = NULL ) optionsLSODA( rtol = 1e-06, atol = 1e-06, jacfunc = NULL, jactype = "fullint", rootfunc = NULL, verbose = FALSE, nroot = 0, tcrit = NULL, hmin = 0, hmax = NULL, hini = 0, ynames = TRUE, maxordn = 12, maxords = 5, bandup = NULL, banddown = NULL, maxsteps = 5000 )
List of arguments
Generate parameter transformation function from a named character vector or object of class eqnvec. This is a wrapper function for Pexpl and Pimpl. See for more details there.
P( trafo = NULL, parameters = NULL, condition = NULL, attach.input = FALSE, keep.root = TRUE, compile = FALSE, modelname = NULL, method = c("explicit", "implicit"), verbose = FALSE )P( trafo = NULL, parameters = NULL, condition = NULL, attach.input = FALSE, keep.root = TRUE, compile = FALSE, modelname = NULL, method = c("explicit", "implicit"), verbose = FALSE )
trafo |
object of class |
parameters |
character vector |
condition |
character, the condition for which the transformation is generated |
attach.input |
attach those incoming parameters to output which are not overwritten by the parameter transformation. |
keep.root |
logical, applies for |
compile |
logical, compile the function (see funC0) |
modelname |
character, see funC0 |
method |
character, either |
verbose |
Print out information during compilation |
An object of class parfn.
Title
P_indiv(p0, est.grid, fix.grid)P_indiv(p0, est.grid, fix.grid)
p0 |
|
est.grid |
|
fix.grid |
Daniel Lill ([email protected])
Generate functions that transform one parameter vector into another by means of a transformation, pushing forward the jacobian matrix of the original parameter. Usually, this function is called internally, e.g. by P. However, you can use it to add your own specialized parameter transformations to the general framework.
parfn(p2p, parameters = NULL, condition = NULL)parfn(p2p, parameters = NULL, condition = NULL)
p2p |
a transformation function for one condition, i.e. a function
|
parameters |
character vector, the parameters accepted by the function |
condition |
character, the condition for which the transformation is defined |
object of class parfn, i.e. a function p(..., fixed, deriv,
conditions, env). The argument pars should be passed via the ...
argument.
Contains attributes "mappings", a list of p2p
functions, "parameters", the union of parameters acceted by the mappings and
"conditions", the total set of conditions.
# Define a time grid on which to make a prediction by peace-wise linear function. # Then define a (generic) prediction function based on thid grid. times <- 0:5 grid <- data.frame(name = "A", time = times, row.names = paste0("p", times)) x <- Xd(grid) # Define an observable and an observation function observables <- eqnvec(Aobs = "s*A") g <- Y(g = observables, f = NULL, states = "A", parameters = "s") # Collect parameters and define an overarching parameter transformation # for two "experimental condtions". dynpars <- attr(x, "parameters") obspars <- attr(g, "parameters") innerpars <- c(dynpars, obspars) trafo <- structure(innerpars, names = innerpars) trafo_C1 <- replaceSymbols(innerpars, paste(innerpars, "C1", sep = "_"), trafo) trafo_C2 <- replaceSymbols(innerpars, paste(innerpars, "C2", sep = "_"), trafo) p <- NULL p <- p + P(trafo = trafo_C1, condition = "C1") p <- p + P(trafo = trafo_C2, condition = "C2") # Collect outer (overarching) parameters and # initialize with random values outerpars <- attr(p, "parameters") pars <- structure(runif(length(outerpars), 0, 1), names = outerpars) # Predict internal/unobserved states out1 <- (x*p)(times, pars) plot(out1) # Predict observed states in addition to unobserved out2 <- (g*x*p)(times, pars) plot(out2)# Define a time grid on which to make a prediction by peace-wise linear function. # Then define a (generic) prediction function based on thid grid. times <- 0:5 grid <- data.frame(name = "A", time = times, row.names = paste0("p", times)) x <- Xd(grid) # Define an observable and an observation function observables <- eqnvec(Aobs = "s*A") g <- Y(g = observables, f = NULL, states = "A", parameters = "s") # Collect parameters and define an overarching parameter transformation # for two "experimental condtions". dynpars <- attr(x, "parameters") obspars <- attr(g, "parameters") innerpars <- c(dynpars, obspars) trafo <- structure(innerpars, names = innerpars) trafo_C1 <- replaceSymbols(innerpars, paste(innerpars, "C1", sep = "_"), trafo) trafo_C2 <- replaceSymbols(innerpars, paste(innerpars, "C2", sep = "_"), trafo) p <- NULL p <- p + P(trafo = trafo_C1, condition = "C1") p <- p + P(trafo = trafo_C2, condition = "C2") # Collect outer (overarching) parameters and # initialize with random values outerpars <- attr(p, "parameters") pars <- structure(runif(length(outerpars), 0, 1), names = outerpars) # Predict internal/unobserved states out1 <- (x*p)(times, pars) plot(out1) # Predict observed states in addition to unobserved out2 <- (g*x*p)(times, pars) plot(out2)
A parameter frame is a data.frame where the rows correspond to different parameter specifications. The columns are divided into three parts. (1) the meta-information columns (e.g. index, value, constraint, etc.), (2) the attributes of an objective function (e.g. data contribution and prior contribution) and (3) the parameters.
parframe( x = NULL, parameters = colnames(x), metanames = NULL, obj.attributes = NULL ) is.parframe(x) ## S3 method for class 'parframe' x[i = NULL, j = NULL, drop = FALSE] ## S3 method for class 'parframe' subset(x, ...)parframe( x = NULL, parameters = colnames(x), metanames = NULL, obj.attributes = NULL ) is.parframe(x) ## S3 method for class 'parframe' x[i = NULL, j = NULL, drop = FALSE] ## S3 method for class 'parframe' subset(x, ...)
x |
data.frame. |
parameters |
character vector, the names of the parameter columns. |
metanames |
character vector, the names of the meta-information columns. |
obj.attributes |
character vector, the names of the objective function attributes. |
i |
row index in any format |
j |
column index in any format |
drop |
logical. If TRUE the result is coerced to the lowest possible dimension |
... |
additional arguments |
Parameter frames can be subsetted either by [ , ] or by subset. If
[ , index] is used, the names of the removed columns will also be removed from
the corresponding attributes, i.e. metanames, obj.attributes and parameters.
An object of class parframe, i.e. a data.frame with attributes for the
different names. Inherits from data.frame.
## Generate a prediction function regfn <- c(y = "sin(a*time)") g <- Y(regfn, parameters = "a") x <- Xt(condition = "C1") ## Generate data data <- datalist( C1 = data.frame( name = "y", time = 1:5, value = sin(1:5) + rnorm(5, 0, .1), sigma = .1 ) ) ## Initialize parameters and time pars <- c(a = 1) times <- seq(0, 5, .1) plot((g*x)(times, pars), data) ## Do many fits from random positions and store them into parlist out <- as.parlist(lapply(1:50, function(i) { trust(normL2(data, g*x), pars + rnorm(length(pars), 0, 1), rinit = 1, rmax = 10) })) summary(out) ## Reduce parlist to parframe parframe <- as.parframe(out) plotValues(parframe) ## Reduce parframe to best fit bestfit <- as.parvec(parframe) plot((g*x)(times, bestfit), data)## Generate a prediction function regfn <- c(y = "sin(a*time)") g <- Y(regfn, parameters = "a") x <- Xt(condition = "C1") ## Generate data data <- datalist( C1 = data.frame( name = "y", time = 1:5, value = sin(1:5) + rnorm(5, 0, .1), sigma = .1 ) ) ## Initialize parameters and time pars <- c(a = 1) times <- seq(0, 5, .1) plot((g*x)(times, pars), data) ## Do many fits from random positions and store them into parlist out <- as.parlist(lapply(1:50, function(i) { trust(normL2(data, g*x), pars + rnorm(length(pars), 0, 1), rinit = 1, rmax = 10) })) summary(out) ## Reduce parlist to parframe parframe <- as.parframe(out) plotValues(parframe) ## Reduce parframe to best fit bestfit <- as.parvec(parframe) plot((g*x)(times, bestfit), data)
The special use of a parameter list is to save the outcome of multiple optimization runs provided by mstrust, into one list.
Fitlists carry an fit index which must be held unique on merging multiple fitlists.
parlist(...) as.parlist(x = NULL) ## S3 method for class 'parlist' summary(object, ...) ## S3 method for class 'parlist' c(...)parlist(...) as.parlist(x = NULL) ## S3 method for class 'parlist' summary(object, ...) ## S3 method for class 'parlist' c(...)
... |
Objects to be coerced to parameter list. |
x |
list of lists, as returned by |
object |
a parlist |
Wolfgang Mader, [email protected]
## Generate a prediction function regfn <- c(y = "sin(a*time)") g <- Y(regfn, parameters = "a") x <- Xt(condition = "C1") ## Generate data data <- datalist( C1 = data.frame( name = "y", time = 1:5, value = sin(1:5) + rnorm(5, 0, .1), sigma = .1 ) ) ## Initialize parameters and time pars <- c(a = 1) times <- seq(0, 5, .1) plot((g*x)(times, pars), data) ## Do many fits from random positions and store them into parlist out <- as.parlist(lapply(1:50, function(i) { trust(normL2(data, g*x), pars + rnorm(length(pars), 0, 1), rinit = 1, rmax = 10) })) summary(out) ## Reduce parlist to parframe parframe <- as.parframe(out) plotValues(parframe) ## Reduce parframe to best fit bestfit <- as.parvec(parframe) plot((g*x)(times, bestfit), data)## Generate a prediction function regfn <- c(y = "sin(a*time)") g <- Y(regfn, parameters = "a") x <- Xt(condition = "C1") ## Generate data data <- datalist( C1 = data.frame( name = "y", time = 1:5, value = sin(1:5) + rnorm(5, 0, .1), sigma = .1 ) ) ## Initialize parameters and time pars <- c(a = 1) times <- seq(0, 5, .1) plot((g*x)(times, pars), data) ## Do many fits from random positions and store them into parlist out <- as.parlist(lapply(1:50, function(i) { trust(normL2(data, g*x), pars + rnorm(length(pars), 0, 1), rinit = 1, rmax = 10) })) summary(out) ## Reduce parlist to parframe parframe <- as.parframe(out) plotValues(parframe) ## Reduce parframe to best fit bestfit <- as.parvec(parframe) plot((g*x)(times, bestfit), data)
A parameter vector is a named numeric vector (the parameter values) together with a "deriv" attribute (the Jacobian of a parameter transformation by which the parameter vector was generated).
parvec(..., deriv = NULL) as.parvec(x, ...) ## S3 method for class 'numeric' as.parvec(x, names = NULL, deriv = NULL, ...) ## S3 method for class 'parvec' x[..., drop = FALSE] ## S3 method for class 'parvec' c(...)parvec(..., deriv = NULL) as.parvec(x, ...) ## S3 method for class 'numeric' as.parvec(x, names = NULL, deriv = NULL, ...) ## S3 method for class 'parvec' x[..., drop = FALSE] ## S3 method for class 'parvec' c(...)
... |
objects to be concatenated |
deriv |
matrix with rownames (according to names of |
x |
numeric or named numeric, the parameter values |
names |
optional character vector, the parameter names. Otherwise, names
are taken from |
drop |
logical, drop empty columns in Jacobian after subsetting. ATTENTION: Be careful with this option. The default behavior is to keep the columns in the Jacobian. This can lead to unintended results when subsetting the parvec and using it e.g. in another parameter transformation. |
An object of class parvec, i.e. a named numeric vector with attribute "deriv".
# Generate a parameter vector v <- parvec(a = 2, b = 3) print(v) print(getDerivs(v)) # Parameter vector from a named numeric M <- matrix(c(1, 1, 0, 1), nrow = 2, ncol = 2, dimnames = list(c("a", "b"), c("A", "B")) ) v <- as.parvec(x = c(a = 2, b = 3), deriv = M) print(v) print(getDerivs(v)) # Subsetting of parameter vectors # Case 1: Dependencies in the Jacobian are maintained w <- v[1] print(w) print(getDerivs(w)) # Case 2: Dependencies are dropped w <- v[1, drop = TRUE] print(w) print(getDerivs(w)) # Concatenating parameter vectors w <- parvec(c = 4, d = 5) print(c(v, w)) print(getDerivs(c(v, w)))# Generate a parameter vector v <- parvec(a = 2, b = 3) print(v) print(getDerivs(v)) # Parameter vector from a named numeric M <- matrix(c(1, 1, 0, 1), nrow = 2, ncol = 2, dimnames = list(c("a", "b"), c("A", "B")) ) v <- as.parvec(x = c(a = 2, b = 3), deriv = M) print(v) print(getDerivs(v)) # Subsetting of parameter vectors # Case 1: Dependencies in the Jacobian are maintained w <- v[1] print(w) print(getDerivs(w)) # Case 2: Dependencies are dropped w <- v[1, drop = TRUE] print(w) print(getDerivs(w)) # Concatenating parameter vectors w <- parvec(c = 4, d = 5) print(c(v, w)) print(getDerivs(c(v, w)))
Parameter transformation
Pexpl( trafo, parameters = NULL, attach.input = FALSE, condition = NULL, compile = FALSE, modelname = NULL, verbose = FALSE )Pexpl( trafo, parameters = NULL, attach.input = FALSE, condition = NULL, compile = FALSE, modelname = NULL, verbose = FALSE )
trafo |
Named character vector. Names correspond to the parameters being fed into the model (the inner parameters). The elements of tafo are equations that express the inner parameters in terms of other parameters (the outer parameters) |
parameters |
Character vector. Optional. If given, the generated parameter
transformation returns values for each element in |
attach.input |
attach those incoming parameters to output which are not overwritten by the parameter transformation. |
condition |
character, the condition for which the transformation is generated |
compile |
Logical, compile the function (see funC0) |
modelname |
Character, used if |
verbose |
Print compiler output to R command line. |
a function p2p(p, fixed = NULL, deriv = TRUE) representing the parameter
transformation. Here, p is a named numeric vector with the values of the outer parameters,
fixed is a named numeric vector with values of the outer parameters being considered
as fixed (no derivatives returned) and deriv is a logical determining whether the Jacobian
of the parameter transformation is returned as attribute "deriv".
Pimpl for implicit parameter transformations
logtrafo <- c(k1 = "exp(logk1)", k2 = "exp(logk2)", A = "exp(logA)", B = "exp(logB)") p_log <- P(logtrafo) pars <- c(logk1 = 1, logk2 = -1, logA = 0, logB = 0) out <- p_log(pars) getDerivs(out)logtrafo <- c(k1 = "exp(logk1)", k2 = "exp(logk2)", A = "exp(logA)", B = "exp(logB)") p_log <- P(logtrafo) pars <- c(logk1 = 1, logk2 = -1, logA = 0, logB = 0) out <- p_log(pars) getDerivs(out)
Parameter transformation (implicit)
Pimpl( trafo, parameters = NULL, condition = NULL, keep.root = TRUE, positive = TRUE, compile = FALSE, modelname = NULL, verbose = FALSE )Pimpl( trafo, parameters = NULL, condition = NULL, keep.root = TRUE, positive = TRUE, compile = FALSE, modelname = NULL, verbose = FALSE )
trafo |
Named character vector defining the equations to be set to zero. Names correspond to dependent variables. |
parameters |
Character vector, the independent variables. |
condition |
character, the condition for which the transformation is generated |
keep.root |
logical, applies for |
positive |
logical, returns projection to the (semi)positive range. Comes with a warning if the steady state has been found to be negative. |
compile |
Logical, compile the function (see funC0) |
modelname |
Character, used if |
verbose |
Print compiler output to R command line. |
Usually, the equations contain the dependent variables, the independent variables and
other parameters. The argument p of p2p must provide values for the independent
variables and the parameters but ALSO FOR THE DEPENDENT VARIABLES. Those serve as initial guess
for the dependent variables. The dependent variables are then numerically computed by
multiroot. The Jacobian of the solution with respect to dependent variables
and parameters is computed by the implicit function theorem. The function p2p returns
all parameters as they are with corresponding 1-entries in the Jacobian.
a function p2p(p, fixed = NULL, deriv = TRUE) representing the parameter
transformation. Here, p is a named numeric vector with the values of the outer parameters,
fixed is a named numeric vector with values of the outer parameters being considered
as fixed (no derivatives returned) and deriv is a logical determining whether the Jacobian
of the parameter transformation is returned as attribute "deriv".
Pexpl for explicit parameter transformations
######################################################################## ## Example 1: Steady-state trafo ######################################################################## f <- c(A = "-k1*A + k2*B", B = "k1*A - k2*B") P.steadyState <- Pimpl(f, "A") p.outerValues <- c(k1 = 1, k2 = 0.1, A = 10, B = 1) P.steadyState(p.outerValues) ######################################################################## ## Example 2: Steady-state trafo combined with log-transform ######################################################################## f <- c(A = "-k1*A + k2*B", B = "k1*A - k2*B") P.steadyState <- Pimpl(f, "A") logtrafo <- c(k1 = "exp(logk1)", k2 = "exp(logk2)", A = "exp(logA)", B = "exp(logB)") P.log <- P(logtrafo) p.outerValue <- c(logk1 = 1, logk2 = -1, logA = 0, logB = 0) (P.log)(p.outerValue) (P.steadyState * P.log)(p.outerValue) ######################################################################## ## Example 3: Steady-states with conserved quantitites ######################################################################## f <- c(A = "-k1*A + k2*B", B = "k1*A - k2*B") replacement <- c(B = "A + B - total") f[names(replacement)] <- replacement pSS <- Pimpl(f, "total") pSS(c(k1 = 1, k2 = 2, A = 5, B = 5, total = 3))######################################################################## ## Example 1: Steady-state trafo ######################################################################## f <- c(A = "-k1*A + k2*B", B = "k1*A - k2*B") P.steadyState <- Pimpl(f, "A") p.outerValues <- c(k1 = 1, k2 = 0.1, A = 10, B = 1) P.steadyState(p.outerValues) ######################################################################## ## Example 2: Steady-state trafo combined with log-transform ######################################################################## f <- c(A = "-k1*A + k2*B", B = "k1*A - k2*B") P.steadyState <- Pimpl(f, "A") logtrafo <- c(k1 = "exp(logk1)", k2 = "exp(logk2)", A = "exp(logA)", B = "exp(logB)") P.log <- P(logtrafo) p.outerValue <- c(logk1 = 1, logk2 = -1, logA = 0, logB = 0) (P.log)(p.outerValue) (P.steadyState * P.log)(p.outerValue) ######################################################################## ## Example 3: Steady-states with conserved quantitites ######################################################################## f <- c(A = "-k1*A + k2*B", B = "k1*A - k2*B") replacement <- c(B = "A + B - total") f[names(replacement)] <- replacement pSS <- Pimpl(f, "total") pSS(c(k1 = 1, k2 = 2, A = 5, B = 5, total = 3))
Plot a list data points
## S3 method for class 'datalist' plot(x, ..., scales = "free", facet = "wrap")## S3 method for class 'datalist' plot(x, ..., scales = "free", facet = "wrap")
x |
Named list of data.frames as being used in res, i.e. with columns |
... |
Further arguments going to |
scales |
The scales argument of |
facet |
Either |
The data.frame being plotted has columns time, value, sigma,
name and condition.
A plot object of class ggplot.
Plot a parameter list.
## S3 method for class 'parlist' plot(x, path = FALSE, ...)## S3 method for class 'parlist' plot(x, path = FALSE, ...)
x |
fitlist obtained from mstrust |
path |
print path of parameters from initials to convergence. For this
option to be TRUE |
... |
additional arguments |
If path=TRUE:
Malenka Mader, [email protected]
Plot an array of trajectories along the profile of a parameter
plotArray( par, profs, prd, times, direction = c("up", "down"), covtable, ..., nsimus = 4 )plotArray( par, profs, prd, times, direction = c("up", "down"), covtable, ..., nsimus = 4 )
par |
Character of parameter name for which the array should be generated. |
profs |
Lists of profiles as being returned by profile. |
prd |
Named list of matrices or data.frames, usually the output of a prediction function as generated by Xs. |
times |
Numeric vector of time points for the model prediction. |
direction |
Character "up" or "down" indicating the direction the value should be traced along the profile starting at the bestfit value. |
covtable |
Optional covariate table or condition.grid necessary if subsetting is required. |
... |
Further arguments for subsetting the plot. |
nsimus |
Number of trajectories/ simulation to be calculated. |
A plot object of class ggplot.
Svenja Kemmer, [email protected]
## Not run: plotArray("myparameter", myprofiles, g*x*p, seq(0, 250, 1), "up", condition.grid, name == "ProteinA" & condition == "c1") ## End(Not run)## Not run: plotArray("myparameter", myprofiles, g*x*p, seq(0, 250, 1), "up", condition.grid, name == "ProteinA" & condition == "c1") ## End(Not run)
Plot a list of model predictions and a list of data points in a combined plot
## S3 method for class 'tbl_df' plotCombined(model, hypothesis = 1, index = 1, ..., plotErrorBands = F) plotCombined(prediction, ...) ## S3 method for class 'prdlist' plot(x, data = NULL, ..., scales = "free", facet = "wrap", transform = NULL) ## S3 method for class 'prdlist' plotCombined( prediction, data = NULL, ..., scales = "free", facet = "wrap", transform = NULL, aesthetics = NULL ) ## S3 method for class 'prdframe' plot(x, data = NULL, ..., scales = "free", facet = "wrap", transform = NULL)## S3 method for class 'tbl_df' plotCombined(model, hypothesis = 1, index = 1, ..., plotErrorBands = F) plotCombined(prediction, ...) ## S3 method for class 'prdlist' plot(x, data = NULL, ..., scales = "free", facet = "wrap", transform = NULL) ## S3 method for class 'prdlist' plotCombined( prediction, data = NULL, ..., scales = "free", facet = "wrap", transform = NULL, aesthetics = NULL ) ## S3 method for class 'prdframe' plot(x, data = NULL, ..., scales = "free", facet = "wrap", transform = NULL)
... |
Further arguments going to |
plotErrorBands |
If the dMod.frame contains an observation function for the error model, use it to show error bands. |
prediction |
Named list of matrices or data.frames, usually the output of a prediction function as generated by Xs. |
x |
prediction |
data |
Named list of data.frames as being used in res, i.e. with columns |
scales |
The scales argument of |
facet |
|
transform |
list of transformation for the states, see coordTransform. |
aesthetics |
Named list of aesthetic mappings, specified as character, e.g. |
The data.frame being plotted has columns time, value, sigma,
name and condition.
A plot object of class ggplot.
## Observation function fn <- eqnvec( sine = "1 + sin(6.28*omega*time)", cosine = "cos(6.28*omega*time)" ) g <- Y(fn, parameters = "omega") ## Prediction function for time x <- Xt() ## Parameter transformations to split conditions p <- NULL for (i in 1:3) { p <- p + P(trafo = c(omega = paste0("omega_", i)), condition = paste0("frequency_", i)) } ## Evaluate prediction times <- seq(0, 1, .01) pars <- structure(seq(1, 2, length.out = 3), names = attr(p, "parameters")) prediction <- (g*x*p)(times, pars) ## Plotting prediction # plot(prediction) plotPrediction(prediction) plotPrediction(prediction, scales = "fixed") plotPrediction(prediction, facet = "grid") plotPrediction(prediction, scales = "fixed", transform = list(sine = "x^2", cosine = "x - 1")) ## Simulate data dataset <- wide2long(prediction) dataset <- dataset[seq(1, nrow(dataset), 5),] set.seed(1) dataset$value <- dataset$value + rnorm(nrow(dataset), 0, .1) dataset$sigma <- 0.1 data <- as.datalist(dataset, split.by = "condition") ## Plotting data # plot(data) plot1 <- plotData(data) plot1 ## Plotting data and prediction with subsetting # plot(prediction, data) plot2 <- plotCombined(prediction, data) plot2 plot3 <- plotCombined(prediction, data, time <= 0.5 & condition == "frequency_1") plot3 plot4 <- plotCombined(prediction, data, time <= 0.5 & condition != "frequency_1", facet = "grid") plot4 plot5 <- plotCombined(prediction, data, aesthetics = list(linetype = "condition")) plot5## Observation function fn <- eqnvec( sine = "1 + sin(6.28*omega*time)", cosine = "cos(6.28*omega*time)" ) g <- Y(fn, parameters = "omega") ## Prediction function for time x <- Xt() ## Parameter transformations to split conditions p <- NULL for (i in 1:3) { p <- p + P(trafo = c(omega = paste0("omega_", i)), condition = paste0("frequency_", i)) } ## Evaluate prediction times <- seq(0, 1, .01) pars <- structure(seq(1, 2, length.out = 3), names = attr(p, "parameters")) prediction <- (g*x*p)(times, pars) ## Plotting prediction # plot(prediction) plotPrediction(prediction) plotPrediction(prediction, scales = "fixed") plotPrediction(prediction, facet = "grid") plotPrediction(prediction, scales = "fixed", transform = list(sine = "x^2", cosine = "x - 1")) ## Simulate data dataset <- wide2long(prediction) dataset <- dataset[seq(1, nrow(dataset), 5),] set.seed(1) dataset$value <- dataset$value + rnorm(nrow(dataset), 0, .1) dataset$sigma <- 0.1 data <- as.datalist(dataset, split.by = "condition") ## Plotting data # plot(data) plot1 <- plotData(data) plot1 ## Plotting data and prediction with subsetting # plot(prediction, data) plot2 <- plotCombined(prediction, data) plot2 plot3 <- plotCombined(prediction, data, time <= 0.5 & condition == "frequency_1") plot3 plot4 <- plotCombined(prediction, data, time <= 0.5 & condition != "frequency_1", facet = "grid") plot4 plot5 <- plotCombined(prediction, data, aesthetics = list(linetype = "condition")) plot5
Plot a list data points
## S3 method for class 'tbl_df' plotData(dMod.frame, hypothesis = 1, ...) ## S3 method for class 'datalist' plotData(data, ..., scales = "free", facet = "wrap", transform = NULL) plotData(data, ...) ## S3 method for class 'data.frame' plotData(data, ...)## S3 method for class 'tbl_df' plotData(dMod.frame, hypothesis = 1, ...) ## S3 method for class 'datalist' plotData(data, ..., scales = "free", facet = "wrap", transform = NULL) plotData(data, ...) ## S3 method for class 'data.frame' plotData(data, ...)
... |
Further arguments going to |
data |
Named list of data.frames as being used in res, i.e. with columns |
scales |
The scales argument of |
facet |
Either |
transform |
list of transformation for the states, see coordTransform. |
The data.frame being plotted has columns time, value, sigma,
name and condition.
A plot object of class ggplot.
## Observation function fn <- eqnvec( sine = "1 + sin(6.28*omega*time)", cosine = "cos(6.28*omega*time)" ) g <- Y(fn, parameters = "omega") ## Prediction function for time x <- Xt() ## Parameter transformations to split conditions p <- NULL for (i in 1:3) { p <- p + P(trafo = c(omega = paste0("omega_", i)), condition = paste0("frequency_", i)) } ## Evaluate prediction times <- seq(0, 1, .01) pars <- structure(seq(1, 2, length.out = 3), names = attr(p, "parameters")) prediction <- (g*x*p)(times, pars) ## Plotting prediction # plot(prediction) plotPrediction(prediction) plotPrediction(prediction, scales = "fixed") plotPrediction(prediction, facet = "grid") plotPrediction(prediction, scales = "fixed", transform = list(sine = "x^2", cosine = "x - 1")) ## Simulate data dataset <- wide2long(prediction) dataset <- dataset[seq(1, nrow(dataset), 5),] set.seed(1) dataset$value <- dataset$value + rnorm(nrow(dataset), 0, .1) dataset$sigma <- 0.1 data <- as.datalist(dataset, split.by = "condition") ## Plotting data # plot(data) plot1 <- plotData(data) plot1 ## Plotting data and prediction with subsetting # plot(prediction, data) plot2 <- plotCombined(prediction, data) plot2 plot3 <- plotCombined(prediction, data, time <= 0.5 & condition == "frequency_1") plot3 plot4 <- plotCombined(prediction, data, time <= 0.5 & condition != "frequency_1", facet = "grid") plot4 plot5 <- plotCombined(prediction, data, aesthetics = list(linetype = "condition")) plot5## Observation function fn <- eqnvec( sine = "1 + sin(6.28*omega*time)", cosine = "cos(6.28*omega*time)" ) g <- Y(fn, parameters = "omega") ## Prediction function for time x <- Xt() ## Parameter transformations to split conditions p <- NULL for (i in 1:3) { p <- p + P(trafo = c(omega = paste0("omega_", i)), condition = paste0("frequency_", i)) } ## Evaluate prediction times <- seq(0, 1, .01) pars <- structure(seq(1, 2, length.out = 3), names = attr(p, "parameters")) prediction <- (g*x*p)(times, pars) ## Plotting prediction # plot(prediction) plotPrediction(prediction) plotPrediction(prediction, scales = "fixed") plotPrediction(prediction, facet = "grid") plotPrediction(prediction, scales = "fixed", transform = list(sine = "x^2", cosine = "x - 1")) ## Simulate data dataset <- wide2long(prediction) dataset <- dataset[seq(1, nrow(dataset), 5),] set.seed(1) dataset$value <- dataset$value + rnorm(nrow(dataset), 0, .1) dataset$sigma <- 0.1 data <- as.datalist(dataset, split.by = "condition") ## Plotting data # plot(data) plot1 <- plotData(data) plot1 ## Plotting data and prediction with subsetting # plot(prediction, data) plot2 <- plotCombined(prediction, data) plot2 plot3 <- plotCombined(prediction, data, time <= 0.5 & condition == "frequency_1") plot3 plot4 <- plotCombined(prediction, data, time <= 0.5 & condition != "frequency_1", facet = "grid") plot4 plot5 <- plotCombined(prediction, data, aesthetics = list(linetype = "condition")) plot5
Plot Fluxes given a list of flux Equations
plotFluxes(pouter, x, times, fluxEquations, nameFlux = "Fluxes:", ...)plotFluxes(pouter, x, times, fluxEquations, nameFlux = "Fluxes:", ...)
pouter |
parameters |
x |
The model prediction function |
times |
Numeric vector of time points for the model prediction |
fluxEquations |
list of chars containing expressions for the fluxes, if names are given, they are shown in the legend. Easy to obtain via subset.eqnlist, see Examples. |
nameFlux |
character, name of the legend. |
... |
Further arguments going to x, such as |
A plot object of class ggplot.
## Not run: plotFluxes(bestfit, x, times, subset(f, "B"%in%Product)$rates, nameFlux = "B production") ## End(Not run)## Not run: plotFluxes(bestfit, x, times, subset(f, "B"%in%Product)$rates, nameFlux = "B production") ## End(Not run)
Plot parameter values for a fitlist
## S3 method for class 'tbl_df' plotPars(dMod.frame, hypothesis = 1, ..., nsteps = 3, tol = 1) ## S3 method for class 'parframe' plotPars(x, tol = 1, ...) plotPars(x, ...)## S3 method for class 'tbl_df' plotPars(dMod.frame, hypothesis = 1, ..., nsteps = 3, tol = 1) ## S3 method for class 'parframe' plotPars(x, tol = 1, ...) plotPars(x, ...)
... |
arguments for subsetting of x |
nsteps |
number of steps from the waterfall plot |
tol |
maximal allowed difference between neighboring objective values to be recognized as one. |
x |
parameter frame as obtained by as.parframe(mstrust) |
Profile likelihood: plot of the parameter paths.
plotPaths.tbl_df(dMod.frame, hypothesis = 1, ...) plotPaths( profs, ..., whichPar = NULL, sort = FALSE, relative = TRUE, scales = "fixed" )plotPaths.tbl_df(dMod.frame, hypothesis = 1, ...) plotPaths( profs, ..., whichPar = NULL, sort = FALSE, relative = TRUE, scales = "fixed" )
hypothesis |
numeric, can be longer than 1 |
... |
arguments going to subset |
profs |
profile or list of profiles as being returned by profile |
whichPar |
Character or index vector, indicating the parameters that are taken as possible reference (x-axis) |
sort |
Logical. If paths from different parameter profiles are plotted together, possible combinations are either sorted or all combinations are taken as they are. |
relative |
logical indicating whether the origin should be shifted. |
scales |
character, either |
See profile for examples.
A plot object of class ggplot.
Profile likelihood: plot all parameter paths belonging to one profile in one plot
plotPathsMulti(profs, whichpars, npars = 5)plotPathsMulti(profs, whichpars, npars = 5)
profs |
Lists of profiles as being returned by profile. |
whichpars |
Character vector of parameter names for which the profile paths should be generated. |
npars |
Numeric vector of number of colored and named parameter paths. |
A plot object of class ggplot for length(whichpars) = 1 and otherwise an object of class cowplot.
Svenja Kemmer, [email protected]
## Not run: plotPathsMulti(myprofiles, c("mypar1", "mypar2"), npars = 5) ## End(Not run)## Not run: plotPathsMulti(myprofiles, c("mypar1", "mypar2"), npars = 5) ## End(Not run)
Plot a list of model predictions
## S3 method for class 'tbl_df' plotPrediction(dMod.frame, hypothesis = 1, index = 1, ...) plotPrediction(prediction, ...) ## S3 method for class 'prdlist' plotPrediction( prediction, ..., errfn = NULL, scales = "free", facet = "wrap", transform = NULL )## S3 method for class 'tbl_df' plotPrediction(dMod.frame, hypothesis = 1, index = 1, ...) plotPrediction(prediction, ...) ## S3 method for class 'prdlist' plotPrediction( prediction, ..., errfn = NULL, scales = "free", facet = "wrap", transform = NULL )
... |
Further arguments going to |
prediction |
Named list of matrices or data.frames, usually the output of a prediction function as generated by Xs. |
errfn |
error model function |
scales |
The scales argument of |
facet |
Either |
transform |
list of transformation for the states, see coordTransform. |
The data.frame being plotted has columns time, value, name and condition.
A plot object of class ggplot.
## Observation function fn <- eqnvec( sine = "1 + sin(6.28*omega*time)", cosine = "cos(6.28*omega*time)" ) g <- Y(fn, parameters = "omega") ## Prediction function for time x <- Xt() ## Parameter transformations to split conditions p <- NULL for (i in 1:3) { p <- p + P(trafo = c(omega = paste0("omega_", i)), condition = paste0("frequency_", i)) } ## Evaluate prediction times <- seq(0, 1, .01) pars <- structure(seq(1, 2, length.out = 3), names = attr(p, "parameters")) prediction <- (g*x*p)(times, pars) ## Plotting prediction # plot(prediction) plotPrediction(prediction) plotPrediction(prediction, scales = "fixed") plotPrediction(prediction, facet = "grid") plotPrediction(prediction, scales = "fixed", transform = list(sine = "x^2", cosine = "x - 1")) ## Simulate data dataset <- wide2long(prediction) dataset <- dataset[seq(1, nrow(dataset), 5),] set.seed(1) dataset$value <- dataset$value + rnorm(nrow(dataset), 0, .1) dataset$sigma <- 0.1 data <- as.datalist(dataset, split.by = "condition") ## Plotting data # plot(data) plot1 <- plotData(data) plot1 ## Plotting data and prediction with subsetting # plot(prediction, data) plot2 <- plotCombined(prediction, data) plot2 plot3 <- plotCombined(prediction, data, time <= 0.5 & condition == "frequency_1") plot3 plot4 <- plotCombined(prediction, data, time <= 0.5 & condition != "frequency_1", facet = "grid") plot4 plot5 <- plotCombined(prediction, data, aesthetics = list(linetype = "condition")) plot5## Observation function fn <- eqnvec( sine = "1 + sin(6.28*omega*time)", cosine = "cos(6.28*omega*time)" ) g <- Y(fn, parameters = "omega") ## Prediction function for time x <- Xt() ## Parameter transformations to split conditions p <- NULL for (i in 1:3) { p <- p + P(trafo = c(omega = paste0("omega_", i)), condition = paste0("frequency_", i)) } ## Evaluate prediction times <- seq(0, 1, .01) pars <- structure(seq(1, 2, length.out = 3), names = attr(p, "parameters")) prediction <- (g*x*p)(times, pars) ## Plotting prediction # plot(prediction) plotPrediction(prediction) plotPrediction(prediction, scales = "fixed") plotPrediction(prediction, facet = "grid") plotPrediction(prediction, scales = "fixed", transform = list(sine = "x^2", cosine = "x - 1")) ## Simulate data dataset <- wide2long(prediction) dataset <- dataset[seq(1, nrow(dataset), 5),] set.seed(1) dataset$value <- dataset$value + rnorm(nrow(dataset), 0, .1) dataset$sigma <- 0.1 data <- as.datalist(dataset, split.by = "condition") ## Plotting data # plot(data) plot1 <- plotData(data) plot1 ## Plotting data and prediction with subsetting # plot(prediction, data) plot2 <- plotCombined(prediction, data) plot2 plot3 <- plotCombined(prediction, data, time <= 0.5 & condition == "frequency_1") plot3 plot4 <- plotCombined(prediction, data, time <= 0.5 & condition != "frequency_1", facet = "grid") plot4 plot5 <- plotCombined(prediction, data, aesthetics = list(linetype = "condition")) plot5
Profile likelihood plot
## S3 method for class 'tbl_df' plotProfile(dMod.frame, hypothesis = 1, ...) ## S3 method for class 'parframe' plotProfile(profs, ..., maxvalue = 5, parlist = NULL) ## S3 method for class 'list' plotProfile(profs, ..., maxvalue = 5, parlist = NULL) plotProfile(profs, ...)## S3 method for class 'tbl_df' plotProfile(dMod.frame, hypothesis = 1, ...) ## S3 method for class 'parframe' plotProfile(profs, ..., maxvalue = 5, parlist = NULL) ## S3 method for class 'list' plotProfile(profs, ..., maxvalue = 5, parlist = NULL) plotProfile(profs, ...)
hypothesis |
numeric, can be longer than 1 |
... |
logical going to subset before plotting. |
profs |
Lists of profiles as being returned by profile. |
maxvalue |
Numeric, the value where profiles are cut off. |
parlist |
Matrix or data.frame with columns for the parameters to be added to the plot as points. If a "value" column is contained, deltas are calculated with respect to lowest chisquare of profiles. |
See profile for examples.
A plot object of class ggplot.
Profile likelihood: plot profiles along with their parameter paths
plotProfilesAndPaths(profs, whichpars, npars = 5)plotProfilesAndPaths(profs, whichpars, npars = 5)
profs |
Lists of profiles as being returned by profile. |
whichpars |
Character vector of parameter names for which the profile paths should be generated. |
npars |
Numeric vector of colored and named parameter paths. |
A plot object of class ggplot.
Svenja Kemmer, [email protected]
## Not run: plotProfilesAndPaths(myprofiles, c("mypar1", "mypar2"), npars = 5) ## End(Not run)## Not run: plotProfilesAndPaths(myprofiles, c("mypar1", "mypar2"), npars = 5) ## End(Not run)
Plot residuals for a fitlist
plotResiduals(parframe, x, data, split = "condition", errmodel = NULL, ...)plotResiduals(parframe, x, data, split = "condition", errmodel = NULL, ...)
parframe |
Object of class |
x |
Prediction function returning named list of data.frames with names as |
data |
Named list of data.frames, i.e. with columns |
split |
List of characters specifying how to summarise the residuals by |
errmodel |
object of type prdfn, the error model function. |
... |
Additional arguments for x |
A plot object of class ggplot with data.frame as attribute attr(P,"out").
## Not run: # time axis: plotResiduals(myfitlist, g*x*p, data, c("time","index","condition","name"), conditions = myconditions[1:4]) # condition axis (residuals summed over time for each observable and condition): plotResiduals(myfitlist, g*x*p, data, c("condition","name","index")) ## End(Not run)## Not run: # time axis: plotResiduals(myfitlist, g*x*p, data, c("time","index","condition","name"), conditions = myconditions[1:4]) # condition axis (residuals summed over time for each observable and condition): plotResiduals(myfitlist, g*x*p, data, c("condition","name","index")) ## End(Not run)
Plotting objective values of a collection of fits
## S3 method for class 'tbl_df' plotValues(dMod.frame, hypothesis = 1, ..., tol = 1) ## S3 method for class 'parframe' plotValues(x, tol = 1, ...) plotValues(x, ...)## S3 method for class 'tbl_df' plotValues(dMod.frame, hypothesis = 1, ..., tol = 1) ## S3 method for class 'parframe' plotValues(x, tol = 1, ...) plotValues(x, ...)
... |
arguments for subsetting of x |
tol |
maximal allowed difference between neighboring objective values to be recognized as one. |
x |
data.frame with columns "value", "converged" and "iterations", e.g. a parframe. |
Title
PRD_indiv(prd0, est.grid, fix.grid)PRD_indiv(prd0, est.grid, fix.grid)
prd0 |
|
est.grid |
|
fix.grid |
Daniel Lill ([email protected])
A prediction function is a function x(..., fixed, deriv, conditions).
Prediction functions are generated by Xs, Xf or Xd. For an example
see the last one.
prdfn(P2X, parameters = NULL, condition = NULL)prdfn(P2X, parameters = NULL, condition = NULL)
P2X |
transformation function as being produced by Xs. |
parameters |
character vector with parameter names |
condition |
character, the condition name |
Prediction functions can be "added" by the "+" operator, see sumfn. Thereby, predictions for different conditions are merged or overwritten. Prediction functions can also be concatenated with other functions, e.g. observation functions (obsfn) or parameter transformation functions (parfn) by the "*" operator, see prodfn.
Object of class prdfn, i.e. a function x(..., fixed, deriv, conditions, env)
which returns a prdlist. The arguments times and
pars (parameter values) should be passed via the ... argument, in this order.
# Define a time grid on which to make a prediction by peace-wise linear function. # Then define a (generic) prediction function based on thid grid. times <- 0:5 grid <- data.frame(name = "A", time = times, row.names = paste0("p", times)) x <- Xd(grid) # Define an observable and an observation function observables <- eqnvec(Aobs = "s*A") g <- Y(g = observables, f = NULL, states = "A", parameters = "s") # Collect parameters and define an overarching parameter transformation # for two "experimental condtions". dynpars <- attr(x, "parameters") obspars <- attr(g, "parameters") innerpars <- c(dynpars, obspars) trafo <- structure(innerpars, names = innerpars) trafo_C1 <- replaceSymbols(innerpars, paste(innerpars, "C1", sep = "_"), trafo) trafo_C2 <- replaceSymbols(innerpars, paste(innerpars, "C2", sep = "_"), trafo) p <- NULL p <- p + P(trafo = trafo_C1, condition = "C1") p <- p + P(trafo = trafo_C2, condition = "C2") # Collect outer (overarching) parameters and # initialize with random values outerpars <- attr(p, "parameters") pars <- structure(runif(length(outerpars), 0, 1), names = outerpars) # Predict internal/unobserved states out1 <- (x*p)(times, pars) plot(out1) # Predict observed states in addition to unobserved out2 <- (g*x*p)(times, pars) plot(out2)# Define a time grid on which to make a prediction by peace-wise linear function. # Then define a (generic) prediction function based on thid grid. times <- 0:5 grid <- data.frame(name = "A", time = times, row.names = paste0("p", times)) x <- Xd(grid) # Define an observable and an observation function observables <- eqnvec(Aobs = "s*A") g <- Y(g = observables, f = NULL, states = "A", parameters = "s") # Collect parameters and define an overarching parameter transformation # for two "experimental condtions". dynpars <- attr(x, "parameters") obspars <- attr(g, "parameters") innerpars <- c(dynpars, obspars) trafo <- structure(innerpars, names = innerpars) trafo_C1 <- replaceSymbols(innerpars, paste(innerpars, "C1", sep = "_"), trafo) trafo_C2 <- replaceSymbols(innerpars, paste(innerpars, "C2", sep = "_"), trafo) p <- NULL p <- p + P(trafo = trafo_C1, condition = "C1") p <- p + P(trafo = trafo_C2, condition = "C2") # Collect outer (overarching) parameters and # initialize with random values outerpars <- attr(p, "parameters") pars <- structure(runif(length(outerpars), 0, 1), names = outerpars) # Predict internal/unobserved states out1 <- (x*p)(times, pars) plot(out1) # Predict observed states in addition to unobserved out2 <- (g*x*p)(times, pars) plot(out2)
A prediction frame is used to store a model prediction in a matrix. The columns of the matrix are "time" and one column per state. The prediction frame has attributes "deriv", the matrix of sensitivities with respect to "outer parameters" (see P), an attribute "sensitivities", the matrix of sensitivities with respect to the "inner parameters" (the model parameters, left-hand-side of the parameter transformation) and an attributes "parameters", the parameter vector of inner parameters to produce the prediction frame.
Prediction frames are usually the constituents of prediction lists (prdlist). They are
produced by Xs, Xd or Xf. When you define your own prediction functions,
see P2X in prdfn, the result should be returned as a prediction frame.
prdframe( prediction = NULL, deriv = NULL, sensitivities = NULL, parameters = NULL )prdframe( prediction = NULL, deriv = NULL, sensitivities = NULL, parameters = NULL )
prediction |
matrix of model prediction |
deriv |
matrix of sensitivities wrt outer parameters |
sensitivities |
matrix of sensitivitie wrt inner parameters |
parameters |
names of the outer paramters |
Object of class prdframe, i.e. a matrix with other matrices and vectors as attributes.
A prediction list is used to store a list of model predictions from different prediction functions or the same prediction function with different parameter specifications. Each entry of the list is a prdframe.
prdlist(...) as.prdlist(x, ...) ## S3 method for class 'list' as.prdlist(x = NULL, names = NULL, ...)prdlist(...) as.prdlist(x, ...) ## S3 method for class 'list' as.prdlist(x = NULL, names = NULL, ...)
... |
objects of class prdframe conditions. |
x |
list of prediction frames |
names |
character vector, the list names, e.g. the names of the experimental |
Make a model prediction for times and a parameter frame. The function is a generalization of the standard prediction by a prediction function object in that it allows to pass a parameter frame instead of a single parameter vector.
## S3 method for class 'prdfn' predict(object, ..., times, pars, data = NULL)## S3 method for class 'prdfn' predict(object, ..., times, pars, data = NULL)
object |
prediction function |
... |
Further arguments goint to the prediction function |
times |
numeric vector of time points |
pars |
|
data |
data list object. If data is passed, its condition.grid
attribute is used to augment the output dataframe by additional
columns. |
A data frame
for nice plots
predtimes(datatimes, eventtimes = NULL, Nobjtimes = 100)predtimes(datatimes, eventtimes = NULL, Nobjtimes = 100)
datatimes |
|
eventtimes |
|
Nobjtimes |
number of time points in total |
vector of time points
predtimes(c(30,60,60,90)) predtimes(c(-30,60,60,90))predtimes(c(30,60,60,90)) predtimes(c(-30,60,60,90))
Print or pander equation list
## S3 method for class 'eqnlist' print(x, pander = FALSE, ...)## S3 method for class 'eqnlist' print(x, pander = FALSE, ...)
x |
object of class eqnlist |
pander |
logical, use pander for output (used with R markdown) |
... |
additional arguments |
Wolfgang Mader, [email protected]
Daniel Kaschek, [email protected]
Print equation vector
## S3 method for class 'eqnvec' print(x, width = 140, pander = FALSE, ...)## S3 method for class 'eqnvec' print(x, width = 140, pander = FALSE, ...)
x |
object of class eqnvec. |
width |
numeric, width of the print-out |
pander |
logical, use pander for output (used with R markdown) |
... |
not used right now |
Wolfgang Mader, [email protected]
Pretty printing parameter transformations
## S3 method for class 'parfn' print(x, ...)## S3 method for class 'parfn' print(x, ...)
x |
prediction function |
... |
additional arguments |
Wolfgang Mader, [email protected]
Pretty printing for a parameter vector
## S3 method for class 'parvec' print(x, ...)## S3 method for class 'parvec' print(x, ...)
x |
object of class |
... |
not used yet. |
Wolfgang Mader, [email protected]
Print object and its "default" attributes only.
print0(x, list_attributes = TRUE)print0(x, list_attributes = TRUE)
x |
Object to be printed |
list_attributes |
Prints the names of all attribute of x, defaults to TRUE |
Before the x is printed by print.default, all its arguments
not in the default list of attrs are removed.
Wolfgang Mader, [email protected]
Mirjam Fehling-Kaschek, [email protected]
As a prior function, it returns derivatives with respect to the penalty parameter in addition to parameter derivatives.
priorL2(mu, lambda = "lambda", attr.name = "prior", condition = NULL)priorL2(mu, lambda = "lambda", attr.name = "prior", condition = NULL)
mu |
Named numeric, the prior values |
lambda |
Character of length one. The name of the penalty paramter in |
attr.name |
character. The constraint value is additionally returned in an attributed with this name |
condition |
character, the condition for which the constraint should apply. If
|
Computes the constraint value
and its derivatives with respect to p and lambda.
List of class objlist, i.e. objective value, gradient and Hessian as list.
wrss
p <- c(A = 1, B = 2, C = 3, lambda = 0) mu <- c(A = 0, B = 0) obj <- priorL2(mu = mu, lambda = "lambda") obj(pars = p + rnorm(length(p), 0, .1))p <- c(A = 1, B = 2, C = 3, lambda = 0) mu <- c(A = 0, B = 0) obj <- priorL2(mu = mu, lambda = "lambda") obj(pars = p + rnorm(length(p), 0, .1))
Profile-likelihood (PL) computation
profile( obj, pars, whichPar, alpha = 0.05, limits = c(lower = -Inf, upper = Inf), method = c("integrate", "optimize"), stepControl = NULL, algoControl = NULL, optControl = NULL, verbose = FALSE, cores = 1, cautiousMode = FALSE, ... )profile( obj, pars, whichPar, alpha = 0.05, limits = c(lower = -Inf, upper = Inf), method = c("integrate", "optimize"), stepControl = NULL, algoControl = NULL, optControl = NULL, verbose = FALSE, cores = 1, cautiousMode = FALSE, ... )
obj |
Objective function |
pars |
Parameter vector corresponding to the log-liklihood optimum. |
whichPar |
Numeric or character vector. The parameters for which the profile is computed. |
alpha |
Numeric, the significance level based on the chisquare distribution with df=1 |
limits |
Numeric vector of length 2, the lower and upper deviance from the original
value of |
method |
Character, either |
stepControl |
List of arguments controlling the step adaption. Defaults to integration set-up, i.e.
|
algoControl |
List of arguments controlling the fast PL algorithm. defaults to
|
optControl |
List of arguments controlling the |
verbose |
Logical, print verbose messages. |
cores |
number of cores used when computing profiles for several parameters. |
cautiousMode |
Logical, write every step to disk and don't delete intermediate results |
... |
Arguments going to obj() |
Computation of the profile likelihood is based on the method of Lagrangian multipliers and Euler integration of the corresponding differential equation of the profile likelihood paths.
algoControl: Since the Hessian which is needed for the differential equation is frequently misspecified,
the error in integration needs to be compensated by a correction factor gamma. Instead of the
Hessian, an identity matrix can be used. To guarantee that the profile likelihood path stays on
the true path, each point proposed by the differential equation can be used as starting point for
an optimization run when reoptimize = TRUE. The correction factor gamma is adapted
based on the amount of actual correction. If this exceeds the value correction, gamma is
reduced. In some cases, the Hessian becomes singular. This leads to problems when inverting the
Hessian. To avoid this problem, the pseudoinverse is computed by removing all singular values lower
than reg.
stepControl: The Euler integration starts with stepsize. In each step the predicted change
of the objective function is compared with the actual change. If this is larger than atol, the
stepsize is reduced. For small deviations, either compared the abolute tolerance atol or the
relative tolerance rtol, the stepsize may be increased. max and min are upper and lower
bounds for stepsize. limit is the maximum number of steps that are take for the profile computation.
stop is a character, usually "value" or "data", for which the significance level alpha
is evaluated.
Named list of length one. The name is the parameter name. The list enty is a matrix with columns "value" (the objective value), "constraint" (deviation of the profiled paramter from the original value), "stepsize" (the stepsize take for the iteration), "gamma" (the gamma value employed for the iteration), "valueData" and "valuePrior" (if specified in obj), one column per parameter (the profile paths).
## Not run: ## Parameter transformation trafo <- eqnvec(a = "exp(loga)", b = "exp(logb)", c = "exp(loga)*exp(logb)*exp(logc)") p <- P(trafo) ## Objective function obj1 <- constraintL2(mu = c(a = .1, b = 1, c = 10), sigma = .6) obj2 <- constraintL2(mu = c(loga = 0, logb = 0), sigma = 10) obj <- obj1*p + obj2 ## Initialize parameters and obtain fit pars <- c(loga = 1, logb = 1, logc = 1) myfit <- trust(obj, pars, rinit = 1, rmax = 10) myfit.fixed <- trust(obj, pars[-1], rinit = 1, rmax = 10, fixed = pars[1]) ## Compute profiles by integration method profiles.approx <- do.call( rbind, lapply(1:3, function(i) { profile(obj, myfit$argument, whichPar = i, limits = c(-3, 3), method = "integrate") }) ) ## Compute profiles by repeated optimization profiles.exact <- do.call( rbind, lapply(1:3, function(i) { profile(obj, myfit$argument, whichPar = i, limits = c(-3, 3), method = "optimize") }) ) ## Compute profiles for fit with fixed element by integration method profiles.approx.fixed <- do.call( rbind, lapply(1:2, function(i) { profile(obj, myfit.fixed$argument, whichPar = i, limits = c(-3, 3), method = "integrate", fixed = pars[1]) }) ) ## Plotting plotProfile(profiles.approx) plotProfile(list(profiles.approx, profiles.exact)) plotProfile(list(profiles.approx, profiles.approx.fixed)) plotPaths(profiles.approx, sort = TRUE) plotPaths(profiles.approx, whichPar = "logc") plotPaths(list(profiles.approx, profiles.approx.fixed), whichPar = "logc") ## Confidence Intervals confint(profiles.approx, val.column = "value") ## End(Not run)## Not run: ## Parameter transformation trafo <- eqnvec(a = "exp(loga)", b = "exp(logb)", c = "exp(loga)*exp(logb)*exp(logc)") p <- P(trafo) ## Objective function obj1 <- constraintL2(mu = c(a = .1, b = 1, c = 10), sigma = .6) obj2 <- constraintL2(mu = c(loga = 0, logb = 0), sigma = 10) obj <- obj1*p + obj2 ## Initialize parameters and obtain fit pars <- c(loga = 1, logb = 1, logc = 1) myfit <- trust(obj, pars, rinit = 1, rmax = 10) myfit.fixed <- trust(obj, pars[-1], rinit = 1, rmax = 10, fixed = pars[1]) ## Compute profiles by integration method profiles.approx <- do.call( rbind, lapply(1:3, function(i) { profile(obj, myfit$argument, whichPar = i, limits = c(-3, 3), method = "integrate") }) ) ## Compute profiles by repeated optimization profiles.exact <- do.call( rbind, lapply(1:3, function(i) { profile(obj, myfit$argument, whichPar = i, limits = c(-3, 3), method = "optimize") }) ) ## Compute profiles for fit with fixed element by integration method profiles.approx.fixed <- do.call( rbind, lapply(1:2, function(i) { profile(obj, myfit.fixed$argument, whichPar = i, limits = c(-3, 3), method = "integrate", fixed = pars[1]) }) ) ## Plotting plotProfile(profiles.approx) plotProfile(list(profiles.approx, profiles.exact)) plotProfile(list(profiles.approx, profiles.approx.fixed)) plotPaths(profiles.approx, sort = TRUE) plotPaths(profiles.approx, whichPar = "logc") plotPaths(list(profiles.approx, profiles.approx.fixed), whichPar = "logc") ## Confidence Intervals confint(profiles.approx, val.column = "value") ## End(Not run)
Generates list of WhichPar entries to facillitate distribute
profile calculation.
profile_pars_per_node(parameters, fits_per_node)profile_pars_per_node(parameters, fits_per_node)
parameters |
list of parameters |
fits_per_node |
numerical, number of parameters that will be send to each node. |
Lists to split the parameters for which the profiles are calculated on the different nodes.
List with two arrays: from contains the number of the starting
parameter, while to stores the respective upper end of the parameter list
per node.
## Not run: parameter_list <- setNames(1:10, letters[1:10]) var_list <- profile_pars_per_node(parameter_list, 4) ## End(Not run)## Not run: parameter_list <- setNames(1:10, letters[1:10]) var_list <- profile_pars_per_node(parameter_list, 4) ## End(Not run)
Progress bar
progressBar(percentage, size = 50, number = TRUE)progressBar(percentage, size = 50, number = TRUE)
percentage |
Numeric between 0 and 100 |
size |
Integer, the size of the bar print-out |
number |
Logical, Indicates whether the percentage should be printed out. |
rPython is liked against a certain Python version found on the system. If Python code called from R requires a specific Python version, the rPython package needs to be reinstalled. This functions helps to do this in one line.
python_version_request(version)python_version_request(version)
version |
character indicating the requested Python version |
TRUE if rPython is linked against the requested version. Otherwise, the user is asked if rPython should be reinstalled with the correctly linked Python version.
Get the Python version to which rPython is linked
python_version_rpython()python_version_rpython()
The Python version and additional information
Check which Python versions are installed on the system
python_version_sys(version = NULL)python_version_sys(version = NULL)
version |
NULL or character. Check for specific version |
Character vector with the python versions and where they are located.
Obtain the mean and standard deviation from replicates per condition.
reduceReplicates(data, select = "condition", datatrans = NULL)reduceReplicates(data, select = "condition", datatrans = NULL)
data |
A data frame containing the measurements. See Format for details. |
select |
Names of the columns in the data frame used to define conditions, see Details. |
datatrans |
Character vector describing a function to transform data. Use x to refer to data. |
The following columns are mandatory for the data frame:
Name of the observed species.
Measurement time point.
Measurement value.
The condition under which the observation was made.
In addition to these columns, any number of columns can follow to allow a
fine-grained definition of conditions. The values of all columns named in
select are then merged to get the set of conditions.
Experiments are usually repeated multiple times possibly under different
conditions leading to replicated measurements. The column "condition" in the
data allows grouping the data by their condition. However, sometimes, a more
fine-grained grouping is desirable. In this case, any number of additional
columns can be appended to the data. These columns are referred to as
"condition identifiers". Which of the condition identifiers are used for
grouping is user-defined by specifying their names in select. The mandatory
column "condition" is always used. The total set of different conditions is
thus defined by all combinations of values occurring in the selected condition
identifiers. The replicates of each condition are then reduced to mean and
standard deviation. New condition names are derived by merging all conditions
which were used in mean and standard deviation.
A data frame of the following variables:
Measurement time point.
Name of the observed species.
Mean of replicates.
Standard error of the mean, NA for single measurements.
The number of replicates reduced.
The condition for which the value and sigma were calculated. If more than one column was used to define the condition, this variable holds the effective condition which is the combination of all applied single conditions.
Wolfgang Mader, [email protected]
Simon Beyer, [email protected]
Method for files (character)
## S3 method for class 'character' reduceReplicates(data, select = "condition", datatrans = NULL)## S3 method for class 'character' reduceReplicates(data, select = "condition", datatrans = NULL)
Method for data frames
## S3 method for class 'data.frame' reduceReplicates(data, select = "condition", datatrans = NULL)## S3 method for class 'data.frame' reduceReplicates(data, select = "condition", datatrans = NULL)
Update attr(prediction, "deriv") to correct par.grid.outer names
renameDerivPars(pred0, pars, est.grid, cn)renameDerivPars(pred0, pars, est.grid, cn)
pred0 |
prd0(times,pars)[1] |
pars |
pars |
est.grid |
est.grid |
cn |
name of condition |
pred0 with updated deriv argument
Daniel Lill ([email protected])
Rename the gradient and hessian of an objlist
renameDerivParsInObjlist(objlist, parnames)renameDerivParsInObjlist(objlist, parnames)
objlist |
an objlist with grid.inner parnames |
parnames |
vector(grid.inner = grid.outer), e.g. as the |
objlist with new names
Daniel Lill ([email protected])
ol <- dMod:::init_empty_objlist(c(a = 1, b = 2)) parnames <- c(a = "c", b = "d") renameDerivParsInObjlist(ol, parnames)ol <- dMod:::init_empty_objlist(c(a = 1, b = 2)) parnames <- c(a = "c", b = "d") renameDerivParsInObjlist(ol, parnames)
Reparameterization
repar(expr, trafo = NULL, ..., reset = FALSE)repar(expr, trafo = NULL, ..., reset = FALSE)
expr |
character of the form |
trafo |
character or equation vector or list thereof. The object where the replacement takes place in |
... |
pass symbols as named arguments |
reset |
logical. If true, the trafo element corresponding to lhs is reset according to rhs. If false, lhs wherever it occurs in the rhs of trafo is replaced by rhs of the formula. |
Left and right-hand side of expr are searched for symbols. If separated by
"_", symbols are recognized as such, e.g. in Delta_x where the symbols are
"Delta" and "x". Each symbol for which values (character or numbers) are passed by the
... argument is replaced.
an equation vector with the reparameterization.
innerpars <- letters[1:3] constraints <- c(a = "b + c") mycondition <- "cond1" trafo <- repar("x ~ x", x = innerpars) trafo <- repar("x ~ y", trafo, x = names(constraints), y = constraints) trafo <- repar("x ~ exp(x)", trafo, x = innerpars) trafo <- repar("x ~ x + Delta_x_condition", trafo, x = innerpars, condition = mycondition)innerpars <- letters[1:3] constraints <- c(a = "b + c") mycondition <- "cond1" trafo <- repar("x ~ x", x = innerpars) trafo <- repar("x ~ y", trafo, x = names(constraints), y = constraints) trafo <- repar("x ~ exp(x)", trafo, x = innerpars) trafo <- repar("x ~ x + Delta_x_condition", trafo, x = innerpars, condition = mycondition)
Wrapper on rep() to input names instead of length.
repWithNames(x, names)repWithNames(x, names)
x |
Value to be repeated. |
names |
List of names. |
When sigma is given in the data, it always has priority over the sigma in the error model
res(data, out, err = NULL)res(data, out, err = NULL)
data |
data.frame with name (factor), time (numeric), value (numeric) and sigma (numeric) |
out |
output of ode(), optionally augmented with attributes "deriv" (output of ode() for the sensitivity equations) and "parameters" (character vector of parameter names, a subsest of those contained in the sensitivity equations). If "deriv" is given, also "parameters" needs to be given. |
err |
output of the error model function |
data.frame with the original data augmented by columns "prediction" ( numeric, the model prediction), "residual" (numeric, difference between prediction and data value), "weighted.residual" (numeric, residual devided by sigma). If "deriv" was given, the returned data.frame has an attribute "deriv" (data.frame with the derivatives of the residuals with respect to the parameters).
Place top elements into bottom elemens
resolveRecurrence(variables)resolveRecurrence(variables)
variables |
named character vector |
If the names of top vector elements occur in the bottom of the vector, they are replaced by the character of the top entry. Useful for steady state conditions.
named character vector of the same length as variables
resolveRecurrence(c(A = "k1*B/k2", C = "A*k3+k4", D="A*C*k5"))resolveRecurrence(c(A = "k1*B/k2", C = "A*k3+k4", D="A*C*k5"))
Transform matrix A into reduced row echelon form this function is written along the lines of the rref-matlab function.
rref(A, tol = sqrt(.Machine$double.eps), verbose = FALSE, fractions = FALSE)rref(A, tol = sqrt(.Machine$double.eps), verbose = FALSE, fractions = FALSE)
A |
matrix for which the reduced row echelon form is searched |
tol |
tolerance to find pivots |
verbose |
logical, print verbose information |
fractions |
logical, not used right now. |
a list of two entries is returned; ret[[1]] is the reduced row echelon form of A, ret[[2]] is the index of columns in which a pivot was found
Malenka Mader, [email protected]
Generate an R code of the expression that is copied via scp
to any machine (ssh-key needed). Then collect the results.
runbg( ..., machine = "localhost", filename = NULL, input = ls(.GlobalEnv), compile = FALSE, wait = FALSE, recover = F )runbg( ..., machine = "localhost", filename = NULL, input = ls(.GlobalEnv), compile = FALSE, wait = FALSE, recover = F )
... |
Some R code |
machine |
Character vector, e.g. |
filename |
Character, defining the filename of the temporary file. Random file name ist chosen if NULL. |
input |
Character vector, the objects in the workspace that are stored into an R data file and copied to the remove machine. |
compile |
Logical. If |
wait |
Logical. Wait until executed. If |
recover |
Logical. This option is useful to recover the three functions check(), get() and purge(), e.g. when a session has crashed. Then, the three functions are recreated without restarting the job. They can then be used to get the results of a job wihtout having to do it manually. Requires the correct filename, so if the previous runbg was run with filename = NULL, you have to specify the tmp_filename manually. |
runbg() generates a workspace from the input argument
and copies the workspace and all C files or .so files to the remote machines via
scp. This will only work if *an ssh-key had been generated and added
to the authorized keys on the remote machine*. The
code snippet, i.e. the ... argument, can include several intermediate results
but only the last call which is not redirected into a variable is returned via the
variable .runbgOutput, see example below.
List of functions check, get() and purge().
check() checks, if the result is ready.
get() copies the result file
to the working directory and loads it into the workspace as an object called .runbgOutput.
This object is a list named according to the machines that contains the results returned by each
machine.
purge() deletes the temporary folder
from the working directory and the remote machines.
## Not run: out_job1 <- runbg({ M <- matrix(rnorm(1e2), 10, 10) solve(M) }, machine = c("localhost", "localhost"), filename = "job1") out_job1$check() out_job1$get() result <- .runbgOutput print(result) out_job1$purge() ## End(Not run) ## Not run: #' Recover a runbg job with the option "recover" out_job1 <- runbg({ M <- matrix(rnorm(1e2), 10, 10) solve(M) }, machine = c("localhost", "localhost"), filename = "job1") Sys.sleep(1) remove(out_job1) try(out_job1$check()) out_job1 <- runbg({ "This code is not run" }, machine = c("localhost", "localhost"), filename = "job1", recover = T) out_job1$get() result <- .runbgOutput print(result) out_job1$purge() ## End(Not run)## Not run: out_job1 <- runbg({ M <- matrix(rnorm(1e2), 10, 10) solve(M) }, machine = c("localhost", "localhost"), filename = "job1") out_job1$check() out_job1$get() result <- .runbgOutput print(result) out_job1$purge() ## End(Not run) ## Not run: #' Recover a runbg job with the option "recover" out_job1 <- runbg({ M <- matrix(rnorm(1e2), 10, 10) solve(M) }, machine = c("localhost", "localhost"), filename = "job1") Sys.sleep(1) remove(out_job1) try(out_job1$check()) out_job1 <- runbg({ "This code is not run" }, machine = c("localhost", "localhost"), filename = "job1", recover = T) out_job1$get() result <- .runbgOutput print(result) out_job1$purge() ## End(Not run)
Generate an R code of the expression that is copied via scp
to the bwForCluster (ssh-key needed). Then collect the results.
runbg_bwfor( ..., machine, filename = NULL, nodes = 1, cores = 1, walltime = "01:00:00", input = ls(.GlobalEnv), compile = TRUE, recover = F )runbg_bwfor( ..., machine, filename = NULL, nodes = 1, cores = 1, walltime = "01:00:00", input = ls(.GlobalEnv), compile = TRUE, recover = F )
... |
Some R code |
machine |
e.g. |
filename |
Character, defining the filename of the temporary file. Random file name ist chosen if NULL. |
nodes |
Number of nodes, e.g. 10 |
cores |
Number of cores, e.g. 16 |
walltime |
estimated runtime in the format |
input |
Character vector, the objects in the workspace that are stored into an R data file and copied to the remove machine. |
compile |
Logical. If |
recover |
Logical, If |
runbg() generates a workspace from the input argument
and copies the workspace and all C files or .so files to the remote machines via
scp. This will only work if *an ssh-key had been generated and added
to the authorized keys on the remote machine*. The
code snippet, i.e. the ... argument, can include several intermediate results
but only the last call which is not redirected into a variable is returned via the
variable .runbgOutput, see example below.
List of functions check(), get() and purge().
check() checks, if the result is ready.
get() copies the result file
to the working directory and loads it into the workspace as an object called .runbgOutput.
This object is a list named according to the machines that contains the results returned by each
machine.
purge() deletes the temporary folder
from the working directory and the remote machines.
## Not run: out_job1 <- runbg({ mstrust(obj, center, fits = 10, cores = 2) }, machine = "bwfor", nodes = 2, cores = "2:best", walltime = "00:01:00", filename = "job1") out_job1$check() out_job1$get() out_job1$purge() result <- .runbgOutput print(result) ## End(Not run)## Not run: out_job1 <- runbg({ mstrust(obj, center, fits = 10, cores = 2) }, machine = "bwfor", nodes = 2, cores = "2:best", walltime = "00:01:00", filename = "job1") out_job1$check() out_job1$get() out_job1$purge() result <- .runbgOutput print(result) ## End(Not run)
Generate an R code of the expression that is copied via scp
to the bwForCluster. Then collect the results. ssh-key not needed. Password can be provided via an additional argument.
sshpass needs to be installed on your local machine.
runbg_bwfor_slurm( ..., machine, filename = NULL, nodes = 1, cores = 1, partition = "single", walltime = "01:00:00", input = ls(.GlobalEnv), compile = TRUE, recover = F, password = "'begin__end'" )runbg_bwfor_slurm( ..., machine, filename = NULL, nodes = 1, cores = 1, partition = "single", walltime = "01:00:00", input = ls(.GlobalEnv), compile = TRUE, recover = F, password = "'begin__end'" )
... |
Some R code |
machine |
e.g. |
filename |
Character, defining the filename of the temporary file. Random file name ist chosen if NULL. Must not contain the string "Minus". |
nodes |
Number of nodes, e.g. 10 |
cores |
Number of cores, e.g. 16 |
partition |
character, the partition where to start the job |
walltime |
estimated runtime in the format |
input |
Character vector, the objects in the workspace that are stored into an R data file and copied to the remove machine. |
compile |
Logical. If |
recover |
Logical, If |
password |
Your ssh password in plain text (yes, no joke unfortunately), the password is handed over to sshpass for automatic login on the cluster. If NULL, the standard ssh/scp is used and you will be asked for your password multiple times while uploading the scripts. |
runbg() generates a workspace from the input argument
and copies the workspace and all C files or .so files to the remote machines via
scp. This will only work if *an ssh-key had been generated and added
to the authorized keys on the remote machine*. The
code snippet, i.e. the ... argument, can include several intermediate results
but only the last call which is not redirected into a variable is returned via the
variable .runbgOutput, see example below.
List of functions check(), get() and purge().
check() checks, if the result is ready.
get() copies the result file
to the working directory and loads it into the workspace as an object called .runbgOutput.
This object is a list named according to the machines that contains the results returned by each
machine.
purge() deletes the temporary folder
from the working directory and the remote machines.
## Not run: out_job1 <- runbg({ mstrust(obj, center, fits = 10, cores = 2) }, machine = "bwfor", nodes = 2, cores = "2:best", walltime = "00:01:00", filename = "job1") out_job1$check() out_job1$get() out_job1$purge() result <- .runbgOutput print(result) ## End(Not run)## Not run: out_job1 <- runbg({ mstrust(obj, center, fits = 10, cores = 2) }, machine = "bwfor", nodes = 2, cores = "2:best", walltime = "00:01:00", filename = "job1") out_job1$check() out_job1$get() out_job1$purge() result <- .runbgOutput print(result) ## End(Not run)
Generate an R code of the expression that is copied via scp
to the bwForCluster. Then collect the results. ssh-key not needed. Password can be provided via an additional argument.
sshpass needs to be installed on your local machine.
runbg_bwfor_sshpass( ..., machine, filename = NULL, nodes = 1, cores = 1, walltime = "01:00:00", input = ls(.GlobalEnv), compile = TRUE, recover = F, password = "'begin__end'" )runbg_bwfor_sshpass( ..., machine, filename = NULL, nodes = 1, cores = 1, walltime = "01:00:00", input = ls(.GlobalEnv), compile = TRUE, recover = F, password = "'begin__end'" )
... |
Some R code |
machine |
e.g. |
filename |
Character, defining the filename of the temporary file. Random file name ist chosen if NULL. |
nodes |
Number of nodes, e.g. 10 |
cores |
Number of cores, e.g. 16 |
walltime |
estimated runtime in the format |
input |
Character vector, the objects in the workspace that are stored into an R data file and copied to the remove machine. |
compile |
Logical. If |
recover |
Logical, If |
password |
Your ssh password in plain text (yes, no joke unfortunately), the password is handed over to sshpass for automatic login on the cluster. |
runbg() generates a workspace from the input argument
and copies the workspace and all C files or .so files to the remote machines via
scp. This will only work if *an ssh-key had been generated and added
to the authorized keys on the remote machine*. The
code snippet, i.e. the ... argument, can include several intermediate results
but only the last call which is not redirected into a variable is returned via the
variable .runbgOutput, see example below.
List of functions check(), get() and purge().
check() checks, if the result is ready.
get() copies the result file
to the working directory and loads it into the workspace as an object called .runbgOutput.
This object is a list named according to the machines that contains the results returned by each
machine.
purge() deletes the temporary folder
from the working directory and the remote machines.
## Not run: out_job1 <- runbg({ mstrust(obj, center, fits = 10, cores = 2) }, machine = "bwfor", nodes = 2, cores = "2:best", walltime = "00:01:00", filename = "job1") out_job1$check() out_job1$get() out_job1$purge() result <- .runbgOutput print(result) ## End(Not run)## Not run: out_job1 <- runbg({ mstrust(obj, center, fits = 10, cores = 2) }, machine = "bwfor", nodes = 2, cores = "2:best", walltime = "00:01:00", filename = "job1") out_job1$check() out_job1$get() out_job1$purge() result <- .runbgOutput print(result) ## End(Not run)
Install your local dMod version to a remote host via ssh.
runbgInstall(sshtarget, source = NULL, type = "local")runbgInstall(sshtarget, source = NULL, type = "local")
sshtarget |
The ssh host url. |
source |
If type = local, source must point to the source directory of your dMod version. This is most probably you local dMod git repository. |
type |
Which dMod to install. At the moment, only your local version is supported. |
Wolfgang Mader, [email protected]
Save the necessary objects for dModtoShiny with the correct names
saveShiny( reactions = myeqnlist, x, parameters, fixed = NULL, data, profiles, pubref = "none", errmodel = NULL )saveShiny( reactions = myeqnlist, x, parameters, fixed = NULL, data, profiles, pubref = "none", errmodel = NULL )
reactions |
eqnlist-object |
x |
prediction funciton |
parameters |
Parframe, result of mstrust() |
fixed |
Character? fixed parameters? |
data |
datalist |
profiles |
Parframe, result of profile(). Is a proflist also possible? |
pubref |
Link to publication |
errmodel |
obsfn |
SaveShiny for dMod.frames
saveShiny_dMod.frame( dMod.frame, hypothesis = 1, reactions = dMod.frame$reactions[[hypothesis]], pubref = "none", fixed = dMod.frame$fixed[[hypothesis]], projectname, appfolder = "/home/mfehling/ShinyApps/dModtoShiny/" )saveShiny_dMod.frame( dMod.frame, hypothesis = 1, reactions = dMod.frame$reactions[[hypothesis]], pubref = "none", fixed = dMod.frame$fixed[[hypothesis]], projectname, appfolder = "/home/mfehling/ShinyApps/dModtoShiny/" )
dMod.frame |
dmod.frame |
hypothesis |
hypothesis |
reactions |
eqnlist which are the basis for prd |
pubref |
character. pubmed link or similar |
fixed |
|
projectname |
Name of the folder being created on fermi |
## Not run: library(dMod) # library(dplyr) # devtools::install_github("dlill/conveniencefunctions") library(conveniencefunctions) setwd(tempdir()) ## Model definition (text-based, scripting part) reactions <- NULL %>% addReaction("A", "B", "k1*A", "translation") %>% addReaction("B", "", "k2*B", "degradation") f <- reactions %>% as.eqnvec() events <- eventlist(var = "A", time = 5, value = "A_add", method = "add") x <- odemodel(f, events = events) %>% Xs g <- Y(c(Bobs = "s1*B"), x, compile = T, modelname = "obsfn") conditions <- c("a", "b") # there is a bug in # getParameters(g*x) parameters <- union(getParameters(g), getParameters(x)) trafo <- NULL %>% define("x~x", x = parameters) %>% branch(conditions = conditions) %>% insert("x~x_cond", x = "s1", cond = condition) %>% insert("x~1", x = "added", conditionMatch = "a") %>% insert("x~5", x = "added", conditionMatch = "b") %>% insert("x~exp(x)", x = getSymbols(mytrafo[[i]])) %>% {.} p <- P(trafo) # Parameter transformation for steady states pSS <- P(f, method = "implicit", condition = NULL, compile = T) ## Process data # Data data <- datalist( a = data.frame(time = c(0, 2, 7), name = "Bobs", value = c(.3, .3, .3), sigma = c(.03, .03, .03)), b = data.frame(time = c(0, 2, 7), name = "Bobs", value = c(.1, .1, .2), sigma = c(.01, .01, .02)) ) myframe <- dMod.frame("no steady states", g, x, p, data) %>% appendObj() %>% mutate(constr = list(constraintL2(mu = 0*pars, sigma = 5)), obj = list(obj_data + constr), fits = list(mstrust(obj, pars, studyname = "Fits", fits = 20, cores = 4, blather = T))) %>% appendParframes() %>% mutate(profiles = list(profile(obj, as.parvec(parframes), whichPar = "k1"))) %>% mutate(vali = list(datapointL2("A", 2, "mypoint", .1, condition = "a")), obj_vali = list(obj_data + constr + vali), par_vali = list(c(dMod:::sanitizePars(as.parvec(parframes))$pars, "mypoint" = 0.1 )), fits_vali = list(mstrust(obj_vali, par_vali)), profile_vali = list(profile(obj_vali, fits_vali %>% as.parframe %>% as.parvec, "mypoint"))) %>% {.} myframe <- myframe %>% tibble::add_column(reactions = list(reactions), fixed = list(NULL)) saveShiny_dMod.frame(myframe, projectname = "dModFrameTest") ## End(Not run)## Not run: library(dMod) # library(dplyr) # devtools::install_github("dlill/conveniencefunctions") library(conveniencefunctions) setwd(tempdir()) ## Model definition (text-based, scripting part) reactions <- NULL %>% addReaction("A", "B", "k1*A", "translation") %>% addReaction("B", "", "k2*B", "degradation") f <- reactions %>% as.eqnvec() events <- eventlist(var = "A", time = 5, value = "A_add", method = "add") x <- odemodel(f, events = events) %>% Xs g <- Y(c(Bobs = "s1*B"), x, compile = T, modelname = "obsfn") conditions <- c("a", "b") # there is a bug in # getParameters(g*x) parameters <- union(getParameters(g), getParameters(x)) trafo <- NULL %>% define("x~x", x = parameters) %>% branch(conditions = conditions) %>% insert("x~x_cond", x = "s1", cond = condition) %>% insert("x~1", x = "added", conditionMatch = "a") %>% insert("x~5", x = "added", conditionMatch = "b") %>% insert("x~exp(x)", x = getSymbols(mytrafo[[i]])) %>% {.} p <- P(trafo) # Parameter transformation for steady states pSS <- P(f, method = "implicit", condition = NULL, compile = T) ## Process data # Data data <- datalist( a = data.frame(time = c(0, 2, 7), name = "Bobs", value = c(.3, .3, .3), sigma = c(.03, .03, .03)), b = data.frame(time = c(0, 2, 7), name = "Bobs", value = c(.1, .1, .2), sigma = c(.01, .01, .02)) ) myframe <- dMod.frame("no steady states", g, x, p, data) %>% appendObj() %>% mutate(constr = list(constraintL2(mu = 0*pars, sigma = 5)), obj = list(obj_data + constr), fits = list(mstrust(obj, pars, studyname = "Fits", fits = 20, cores = 4, blather = T))) %>% appendParframes() %>% mutate(profiles = list(profile(obj, as.parvec(parframes), whichPar = "k1"))) %>% mutate(vali = list(datapointL2("A", 2, "mypoint", .1, condition = "a")), obj_vali = list(obj_data + constr + vali), par_vali = list(c(dMod:::sanitizePars(as.parvec(parframes))$pars, "mypoint" = 0.1 )), fits_vali = list(mstrust(obj_vali, par_vali)), profile_vali = list(profile(obj_vali, fits_vali %>% as.parframe %>% as.parvec, "mypoint"))) %>% {.} myframe <- myframe %>% tibble::add_column(reactions = list(reactions), fixed = list(NULL)) saveShiny_dMod.frame(myframe, projectname = "dModFrameTest") ## End(Not run)
Standard dMod color palette
scale_color_dMod(...)scale_color_dMod(...)
... |
arguments goint to codescale_color_manual() |
library(ggplot2) times <- seq(0, 2*pi, 0.1) values <- sin(times) data <- data.frame( time = times, value = c(values, 1.2*values, 1.4*values, 1.6*values), group = rep(c("C1", "C2", "C3", "C4"), each = length(times)) ) qplot(time, value, data = data, color = group, geom = "line") + theme_dMod() + scale_color_dMod()library(ggplot2) times <- seq(0, 2*pi, 0.1) values <- sin(times) data <- data.frame( time = times, value = c(values, 1.2*values, 1.4*values, 1.6*values), group = rep(c("C1", "C2", "C3", "C4"), each = length(times)) ) qplot(time, value, data = data, color = group, geom = "line") + theme_dMod() + scale_color_dMod()
Standard dMod color scheme
scale_fill_dMod(...)scale_fill_dMod(...)
... |
arguments goint to codescale_color_manual() |
This function performs the scale dependent flux sensitivity analysis for a simulated system. It prints a pdf-file to the working directory for each condition of the system.
scaledependent( frame, x, obstates, bestfit, time, p, bound = 0.001, figurename = "scaledependent", timeunit = "h", width = 16, height = 9, linewidth = 2 )scaledependent( frame, x, obstates, bestfit, time, p, bound = 0.001, figurename = "scaledependent", timeunit = "h", width = 16, height = 9, linewidth = 2 )
frame |
Reactionlist of the model |
x |
Xs function of the model |
obstates |
The states under observation whose changes are analysed |
bestfit |
The parameter values which are used for the simulation |
time |
Array of observed timepoints |
p |
P function of the system |
bound |
Absolute value which a flux sensitive component has to reach at any point in time to be shown in the outputplots |
figurename |
Name of the saved pdf |
timeunit |
Unit of the time at the x-axis of the plots |
linewidth |
Linewidth of the coloured fluxes in the plots |
time, width
|
Dimension of the plots |
#Scale dependent flux sensitivity Analysis for the SIR-Model with two different values for beta frame <- NULL frame <-addReaction(frame, from = "S", to = "I", rate = "beta*S*I") frame <-addReaction(frame, from = "I", to = "R", rate = "gamma*I") x <-Xs(odemodel(frame)) obstates = c("S", "I", "R") par <- getParameters(x) p <- NULL for (i in c(0.0026, 0.0027)){ trafo <- repar("x~x", x=par) trafo <- repar("x~y", x = "beta", y= i, trafo) parameter <- structure(trafo, names = par) p <- p+P(parameter, condition= paste0(i))} time <- seq(0, 14, length.out=1000) bestfit <- c(S = 7.637603e+02, I = 6.184001e-01, R = 8.251977e-16, gamma = 4.582934e-01) scaledependent(frame, x, obstates, bestfit, time, p)#Scale dependent flux sensitivity Analysis for the SIR-Model with two different values for beta frame <- NULL frame <-addReaction(frame, from = "S", to = "I", rate = "beta*S*I") frame <-addReaction(frame, from = "I", to = "R", rate = "gamma*I") x <-Xs(odemodel(frame)) obstates = c("S", "I", "R") par <- getParameters(x) p <- NULL for (i in c(0.0026, 0.0027)){ trafo <- repar("x~x", x=par) trafo <- repar("x~y", x = "beta", y= i, trafo) parameter <- structure(trafo, names = par) p <- p+P(parameter, condition= paste0(i))} time <- seq(0, 14, length.out=1000) bestfit <- c(S = 7.637603e+02, I = 6.184001e-01, R = 8.251977e-16, gamma = 4.582934e-01) scaledependent(frame, x, obstates, bestfit, time, p)
This function performs the scale invariant flux sensitivity analysis for a simulated system. It prints a pdf-file to the working directory for each condition of the system.
scaleinvariant( frame, x, obstates, bestfit, time, p, bound = 0.001, figurename = "scaleinvariant", timeunit = "h", width = 16, height = 9, linewidth = 2 )scaleinvariant( frame, x, obstates, bestfit, time, p, bound = 0.001, figurename = "scaleinvariant", timeunit = "h", width = 16, height = 9, linewidth = 2 )
frame |
Reactionlist of the model |
x |
Xs function of the Model |
obstates |
The states under observation whose changes are analysed |
bestfit |
The parameter values which are used for the simulation |
time |
Array of observed timepoints |
p |
P function of the system |
bound |
Absolute value which a flux sensitive component has to reach at any point in time to be shown in the outputplots |
figurename |
Name of the saved pdf |
timeunit |
Unit of the time at the x-axis of the plots |
linewidth |
Linewidth of the coloured fluxes in the plots |
time, width
|
Dimensions of the plots |
#Scale invariant flux sensitivity Analysis for the SIR-Model with two different values for beta frame <- NULL frame <-addReaction(frame, from = "S", to = "I", rate = "beta*S*I") frame <-addReaction(frame, from = "I", to = "R", rate = "gamma*I") x <-Xs(odemodel(frame)) obstates = c("S", "I", "R") par <- getParameters(x) p <- NULL for (i in c(0.0026, 0.0027)){ trafo <- repar("x~x", x=par) trafo <- repar("x~y", x = "beta", y= i, trafo) parameter <- structure(trafo, names = par) p <- p+P(parameter, condition= paste0(i))} time <- seq(0, 14, length.out=1000) bestfit <- c(S = 7.637603e+02, I = 6.184001e-01, R = 8.251977e-16, gamma = 4.582934e-01) scaleinvariant(frame, x, obstates, bestfit, time, p)#Scale invariant flux sensitivity Analysis for the SIR-Model with two different values for beta frame <- NULL frame <-addReaction(frame, from = "S", to = "I", rate = "beta*S*I") frame <-addReaction(frame, from = "I", to = "R", rate = "gamma*I") x <-Xs(odemodel(frame)) obstates = c("S", "I", "R") par <- getParameters(x) p <- NULL for (i in c(0.0026, 0.0027)){ trafo <- repar("x~x", x=par) trafo <- repar("x~y", x = "beta", y= i, trafo) parameter <- structure(trafo, names = par) p <- p+P(parameter, condition= paste0(i))} time <- seq(0, 14, length.out=1000) bestfit <- c(S = 7.637603e+02, I = 6.184001e-01, R = 8.251977e-16, gamma = 4.582934e-01) scaleinvariant(frame, x, obstates, bestfit, time, p)
Gather statistics of a fitlist
stat.parlist(x)stat.parlist(x)
x |
The fitlist |
This function follows the method published in [1]. Find the latest version of the tool and some examples under [2]. The determined steady-state solution is tailored to parameter estimation. Please note that kinetic parameters might be fixed for solution of steady-state equations. Note that additional parameters might be introduced to ensure positivity of the solution.
The function calls a python script via the reticulate package. Use python3.x
steadyStates( model, file = NULL, rates = NULL, forcings = NULL, givenCQs = NULL, neglect = NULL, sparsifyLevel = 2, outputFormat = "R", testSteady = "T" )steadyStates( model, file = NULL, rates = NULL, forcings = NULL, givenCQs = NULL, neglect = NULL, sparsifyLevel = 2, outputFormat = "R", testSteady = "T" )
model |
Either name of the csv-file or the eqnlist of the model. |
file |
Name of the file to which the steady-state equations are saved. |
rates |
Character vector, flux vector of the system |
forcings |
Character vector with the names of the forcings |
givenCQs |
(Unnamed) Character vector with conserved quantities. Use the format c("A + pA = totA", "B + pB = totB"). The format c("A + pA", "B + pB") works also. If NULL, conserved quantities are automatically calculated. |
neglect |
Character vector with names of states and parameters that must not be used for solving the steady-state equations |
sparsifyLevel |
numeric, Upper bound for length of linear combinations used for simplifying the stoichiometric matrix |
outputFormat |
Define the output format. By default "R" generating dMod compatible output. To obtain an output appropriate for d2d [2] "M" must be selected. |
testSteady |
Boolean, if "T" the correctness of the obtained steady states is numerically checked (this can be very time intensive). If "F" this is skipped. |
Character vector of steady-state equations.
Marcus Rosenblatt, [email protected]
[1] https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4863410/
[2] https://github.com/marcusrosenblatt/AlyssaPetit
[3] https://github.com/Data2Dynamics/d2d
## Not run: library(dMod) setwd(tempdir()) reactions <- eqnlist() reactions <- addReaction(reactions, "Tca_buffer", "Tca_cyto", "import_Tca*Tca_buffer", "Basolateral uptake") reactions <- addReaction(reactions, "Tca_cyto", "Tca_buffer", "export_Tca_baso*Tca_cyto", "Basolateral efflux") reactions <- addReaction(reactions, "Tca_cyto", "Tca_canalicular", "export_Tca_cana*Tca_cyto", "Canalicular efflux") reactions <- addReaction(reactions, "Tca_canalicular", "Tca_buffer", "transport_Tca*Tca_canalicular", "Transport bile") mysteadies <- steadyStates(reactions) print(mysteadies) x <- Xs(odemodel(reactions)) parameters <- getParameters(x) trafo <- `names<-`(parameters, parameters) trafo <- repar("inner~steadyEqn", trafo, inner = names(mysteadies), steadyEqn = mysteadies) pSS <- P(trafo, condition = "steady") set.seed(2) pars <- structure(runif( length(getParameters(pSS)), 0,1), names = getParameters(pSS)) prediction <- (x*pSS)(seq(0,10, 0.1), pars, deriv = F) plot(prediction) ## End(Not run)## Not run: library(dMod) setwd(tempdir()) reactions <- eqnlist() reactions <- addReaction(reactions, "Tca_buffer", "Tca_cyto", "import_Tca*Tca_buffer", "Basolateral uptake") reactions <- addReaction(reactions, "Tca_cyto", "Tca_buffer", "export_Tca_baso*Tca_cyto", "Basolateral efflux") reactions <- addReaction(reactions, "Tca_cyto", "Tca_canalicular", "export_Tca_cana*Tca_cyto", "Canalicular efflux") reactions <- addReaction(reactions, "Tca_canalicular", "Tca_buffer", "transport_Tca*Tca_canalicular", "Transport bile") mysteadies <- steadyStates(reactions) print(mysteadies) x <- Xs(odemodel(reactions)) parameters <- getParameters(x) trafo <- `names<-`(parameters, parameters) trafo <- repar("inner~steadyEqn", trafo, inner = names(mysteadies), steadyEqn = mysteadies) pSS <- P(trafo, condition = "steady") set.seed(2) pars <- structure(runif( length(getParameters(pSS)), 0,1), names = getParameters(pSS)) prediction <- (x*pSS)(seq(0,10, 0.1), pars, deriv = F) plot(prediction) ## End(Not run)
Uses julia to calculate the steady state transformations
steadyStateToolJulia(el, forcings = NULL, neglect = NULL, verboseLevel = 1)steadyStateToolJulia(el, forcings = NULL, neglect = NULL, verboseLevel = 1)
el |
the equation list |
forcings |
vector of strings, default |
neglect |
vector of strings, default |
verboseLevel |
integer, default |
named vector with the steady state transformations. The names are the inner, the values are the outer parameters
Elide character vector
strelide(string, width, where = "right", force = FALSE)strelide(string, width, where = "right", force = FALSE)
string |
String subject to eliding |
width |
Width including eliding ... of return string |
where |
Eliding can happen at 'left', 'middel', or 'right'. Defaults to 'right'. |
force |
Elide, even is <string> is shorter than <width>. Default to 'FALSE'. |
Elide a string to <width>. Eliding can happen at 'left', 'middle', or 'right'. #' If forcing = FALSE, which is the default, strings shorten than <width> are returend unaltered; forcing = TRUE inserts eliding symbols (...) in any case.
Elided string of length <width>.
Wolfgang Mader, [email protected]
Pad string to desired width
strpad(string, width, where = "right", padding = " ", autoelide = FALSE)strpad(string, width, where = "right", padding = " ", autoelide = FALSE)
string |
String to pad |
width |
Desired width of padded string |
where |
Padding can be inserted to the right or left of <string>. Default to 'right'. |
padding |
A single character with with the padding space is filled. Defaults to blank ' ' yielding invisible padding. |
autoelide |
If TRUE, <string> is elided if it is wider than <width>. The position of eliding follows <where>. Defaults to FALSE. |
Padded string of length <width>.
Wolfgang Mader, [email protected]
Submatrix of a matrix returning ALWAYS a matrix
submatrix(M, rows = 1:nrow(M), cols = 1:ncol(M))submatrix(M, rows = 1:nrow(M), cols = 1:ncol(M))
M |
matrix |
rows |
Index vector |
cols |
Index vector |
The matrix M[rows, cols], keeping/adjusting attributes like ncol nrow and dimnames.
subset of an equation list
## S3 method for class 'eqnlist' subset(x, ...)## S3 method for class 'eqnlist' subset(x, ...)
x |
the equation list |
... |
logical expression for subsetting |
The argument ... can contain "Educt", "Product", "Rate" and "Description".
The "
An object of class eqnlist
reactions <- data.frame(Description = c("Activation", "Deactivation"), Rate = c("act*A", "deact*pA"), A=c(-1,1), pA=c(1, -1) ) f <- as.eqnlist(reactions) subset(f, "A" %in% Educt) subset(f, "pA" %in% Product) subset(f, grepl("act", Rate))reactions <- data.frame(Description = c("Activation", "Deactivation"), Rate = c("act*A", "deact*pA"), A=c(-1,1), pA=c(1, -1) ) f <- as.eqnlist(reactions) subset(f, "A" %in% Educt) subset(f, "pA" %in% Product) subset(f, grepl("act", Rate))
Remove redundancies (happens when a parameter is duplicated and mapped to the same outer parameter). For example in est.grid = data.table(ID = 1, init_A = "Atot", Atot= "Atot", ...)
sumDuplicatedParsInDeriv(der)sumDuplicatedParsInDeriv(der)
der |
deriv of a single condition, as used in [renameDerivPars] |
der with redundant columns summed and duplicates removed
Remove redundancies (happens when a parameter is duplicated and mapped to the same outer parameter). For example in est.grid = data.table(ID = 1, init_A = "Atot", Atot= "Atot", ...)
sumDuplicatedParsInObjlist(ol)sumDuplicatedParsInObjlist(ol)
objlist |
objlist with potentially duplicated names |
objlist with duplicated gradient and hessian elements summed and redundancies removed
Daniel Lill ([email protected])
ol <- dMod:::init_empty_objlist(c("S2" = 2, "S3" = 3, S2 = 3)) ol$gradient <- ol$gradient + 1:3 ol$hessian <- ol$hessian + 1:9 sumDuplicatedParsInObjlist(ol)ol <- dMod:::init_empty_objlist(c("S2" = 2, "S3" = 3, S2 = 3)) ol$gradient <- ol$gradient + 1:3 ol$hessian <- ol$hessian + 1:9 sumDuplicatedParsInObjlist(ol)
Summary of an equation vector
## S3 method for class 'eqnvec' summary(object, ...)## S3 method for class 'eqnvec' summary(object, ...)
object |
of class eqnvec. |
... |
additional arguments |
Wolfgang Mader, [email protected]
This function follows the method published in [1].
The function calls a python script via rPython. Usage problems might occur when different python versions are used. The script was written and tested for python 2.7.12, sympy 0.7.6.
Recently, users went into problems with RJSONIO when rPython was used. Unless a sound solution is available, please try to reinstall RJSONIO in these cases.
symmetryDetection( f, obsvect = NULL, prediction = NULL, initial = NULL, ansatz = "uni", pMax = 2, inputs = NULL, fixed = NULL, cores = 1, allTrafos = FALSE )symmetryDetection( f, obsvect = NULL, prediction = NULL, initial = NULL, ansatz = "uni", pMax = 2, inputs = NULL, fixed = NULL, cores = 1, allTrafos = FALSE )
f |
object containing the ODE for which |
obsvect |
vector of observation functions |
prediction |
vector containing prediction to be tested |
initial |
vector containing initial values |
ansatz |
type of infinitesimal ansatz used for the analysis (uni, par, multi) |
pMax |
maximal degree of infinitesimal ansatz |
inputs |
specify the input variables |
fixed |
variables to concider fixed |
cores |
maximal number of cores used for the analysis |
allTrafos |
do not remove transformations with a common parameter factor |
[1] https://journals.aps.org/pre/abstract/10.1103/PhysRevE.92.012920
## Not run: eq <- NULL eq <- addReaction(eq, "A", "B", "k1*A") eq <- addReaction(eq, "B", "A", "k2*B") observables <- eqnvec(Aobs = "alpha * A") symmetryDetection(eq, observables) ## End(Not run)## Not run: eq <- NULL eq <- addReaction(eq, "A", "B", "k1*A") eq <- addReaction(eq, "B", "A", "k2*B") observables <- eqnvec(Aobs = "alpha * A") symmetryDetection(eq, observables) ## End(Not run)
Standard plotting theme of dMod
theme_dMod(base_size = 12, base_family = "")theme_dMod(base_size = 12, base_family = "")
base_size |
numeric, font-size |
base_family |
character, font-name |
DatapointL2 without env access
timepointL2_indiv( name, time, value, sigma = 1, attr.name = "timepointL2", condition, prd_indiv )timepointL2_indiv( name, time, value, sigma = 1, attr.name = "timepointL2", condition, prd_indiv )
name |
|
time |
|
value |
|
sigma |
|
attr.name |
|
condition |
|
prd_indiv |
Daniel Lill ([email protected])
This function carries out a minimization or maximization of a function using a trust region algorithm. See the references for details.
trust( objfun, parinit, rinit, rmax, parscale, iterlim = 100, fterm = sqrt(.Machine$double.eps), mterm = sqrt(.Machine$double.eps), minimize = TRUE, blather = FALSE, parupper = Inf, parlower = -Inf, printIter = FALSE, traceFile = NULL, ... ) trustL1( objfun, parinit, mu = 0 * parinit, one.sided = FALSE, lambda = 1, rinit, rmax, parscale, iterlim = 100, fterm = sqrt(.Machine$double.eps), mterm = sqrt(.Machine$double.eps), minimize = TRUE, blather = FALSE, blather2 = FALSE, parupper = Inf, parlower = -Inf, printIter = FALSE, ... )trust( objfun, parinit, rinit, rmax, parscale, iterlim = 100, fterm = sqrt(.Machine$double.eps), mterm = sqrt(.Machine$double.eps), minimize = TRUE, blather = FALSE, parupper = Inf, parlower = -Inf, printIter = FALSE, traceFile = NULL, ... ) trustL1( objfun, parinit, mu = 0 * parinit, one.sided = FALSE, lambda = 1, rinit, rmax, parscale, iterlim = 100, fterm = sqrt(.Machine$double.eps), mterm = sqrt(.Machine$double.eps), minimize = TRUE, blather = FALSE, blather2 = FALSE, parupper = Inf, parlower = -Inf, printIter = FALSE, ... )
objfun |
an R function that computes value, gradient, and Hessian of the
function to be minimized or maximized and returns them as a list with
components value, gradient, and hessian. Its first argument should be a
vector of the length of parinit followed by any other arguments specified
by the |
parinit |
starting parameter values for the optimization. Must be feasible (in the domain). |
rinit |
starting trust region radius. The trust region radius (see details below) is adjusted as the algorithm proceeds. A bad initial value wastes a few steps while the radius is adjusted, but does not keep the algorithm from working properly. |
rmax |
maximum allowed trust region radius. This may be set very large. If set small, the algorithm traces a steepest descent path (steepest ascent, when minimize = FALSE). |
parscale |
an estimate of the size of each parameter at the minimum. The algorithm operates as if optimizing function(x, ...) objfun(x / parscale, ...). May be missing in which case no rescaling is done. See also the details section below. |
iterlim |
a positive integer specifying the maximum number of iterations to be performed before the program is terminated. |
fterm |
a positive scalar giving the tolerance at which the difference in objective function values in a step is considered close enough to zero to terminate the algorithm. |
mterm |
a positive scalar giving the tolerance at which the two-term Taylor-series approximation to the difference in objective function values in a step is considered close enough to zero to terminate the algorithm. |
minimize |
If TRUE minimize. If FALSE maximize. |
blather |
If TRUE return extra info. |
parupper |
named numeric vector of upper bounds. If not named, first value will be used for all parameters. |
parlower |
named numeric vector of lower bounds. If not named, first value will be used for all parameters. |
printIter |
print iteration information to R console |
... |
additional argument to objfun |
mu |
named numeric value. The reference value for L1 penalized parameters. |
one.sided |
logical. One-sided penalization. |
lambda |
strength of the L1 penalty |
blather2 |
even more information |
See Fletcher (1987, Section 5.1) or Nocedal and Wright (1999, Section 4.2) for detailed expositions.
A list containing the following components:
value: the value returned by objfun at the final iterate.
gradient: the gradient returned by objfun at the final iterate.
hessian: the Hessian returned by objfun at the final iterate.
argument: the final iterate
converged: if TRUE the final iterate was deemed optimal by the specified termination criteria.
iterations: number of trust region subproblems done (including those whose solutions are not accepted).
argpath: (if blather == TRUE) the sequence of iterates, not including the final iterate.
argtry: (if blather == TRUE) the sequence of solutions of the trust region subproblem.
steptype: (if blather == TRUE) the sequence of cases that arise in solutions of the trust region subproblem. "Newton" means the Newton step solves the subproblem (lies within the trust region). Other values mean the subproblem solution is constrained. "easy-easy" means the eigenvectors corresponding to the minimal eigenvalue of the rescaled Hessian are not all orthogonal to the gradient. The other cases are rarely seen. "hard-hard" means the Lagrange multiplier for the trust region constraint is minus the minimal eigenvalue of the rescaled Hessian; "hard-easy" means it isn't.
accept: (if blather == TRUE) indicates which of the sequence of solutions of the trust region subproblem were accepted as the next iterate. (When not accepted the trust region radius is reduced, and the previous iterate is kept.)
r: (if blather == TRUE) the sequence of trust region radii.
rho: (if blather == TRUE) the sequence of ratios of actual over predicted decrease in the objective function in the trust region subproblem, where predicted means the predicted decrease in the two-term Taylor series model used in the subproblem.
valpath: (if blather == TRUE) the sequence of objective function values at the iterates.
valtry: (if blather == TRUE) the sequence of objective function values at the solutions of the trust region subproblem.
preddiff: (if blather == TRUE) the sequence of predicted differences using the two-term Taylor-series model between the function values at the current iterate and at the solution of the trust region subproblem.
stepnorm: (if blather == TRUE) the sequence of norms of steps, that is distance between current iterate and proposed new iterate found in the trust region subproblem.
Extract those lines of a parameter frame with unique elements in the value column
## S3 method for class 'parframe' unique(x, incomparables = FALSE, tol = 1, ...)## S3 method for class 'parframe' unique(x, incomparables = FALSE, tol = 1, ...)
x |
parameter frame |
incomparables |
not used. Argument exists for compatibility with S3 generic. |
tol |
tolerance to decide when values are assumed to be equal, see |
... |
additional arguments being passed to |
A subset of the parameter frame x.
Get variance-covariance matrix from trust result
vcov(fit, parupper = NULL, parlower = NULL)vcov(fit, parupper = NULL, parlower = NULL)
fit |
return from trust or selected fit from parlist |
parupper |
upper limit for parameter values. Parameters are considered fixed when exceeding these limits. |
parlower |
lower limit for parameter values. Parameters are considered fixed when exceeding these limits. |
Translate wide output format (e.g. from ode) into long format
wide2long(out, keep = 1, na.rm = FALSE)wide2long(out, keep = 1, na.rm = FALSE)
out |
data.frame or matrix or list of matrices in wide format |
keep |
Index vector, the columns to keep |
na.rm |
Logical, if |
The function assumes that out[,1] represents a time-like vector whereas out[,-1] represents the values. Useful for plotting with ggplot. If a list is supplied, the names of the list are added as extra column names "condition"
data.frame in long format, i.e. columns "time" (out[,1]), "name" (colnames(out[,-1])), "value" (out[,-1]) and, if out was a list, "condition" (names(out))
Translate wide output format (e.g. from ode) into long format
## S3 method for class 'data.frame' wide2long(out, keep = 1, na.rm = FALSE)## S3 method for class 'data.frame' wide2long(out, keep = 1, na.rm = FALSE)
out |
data.frame or matrix or list of matrices in wide format |
keep |
Index vector, the columns to keep |
na.rm |
Logical, if |
The function assumes that out[,1] represents a time-like vector whereas out[,-1] represents the values. Useful for plotting with ggplot. If a list is supplied, the names of the list are added as extra column names "condition"
data.frame in long format, i.e. columns "time" (out[,1]), "name" (colnames(out[,-1])), "value" (out[,-1]) and, if out was a list, "condition" (names(out))
Translate wide output format (e.g. from ode) into long format
## S3 method for class 'list' wide2long(out, keep = 1, na.rm = FALSE)## S3 method for class 'list' wide2long(out, keep = 1, na.rm = FALSE)
out |
list of matrices in wide format |
keep |
Index vector, the columns to keep |
na.rm |
Logical, if |
The function assumes that out[,1] represents a time-like vector whereas out[,-1] represents the values. Useful for plotting with ggplot. If a list is supplied, the names of the list are added as extra column names "condition"
data.frame in long format, i.e. columns "time" (out[,1]), "name" (colnames(out[,-1])), "value" (out[,-1]) and, if out was a list, "condition" (names(out))
Translate wide output format (e.g. from ode) into long format
## S3 method for class 'matrix' wide2long(out, keep = 1, na.rm = FALSE)## S3 method for class 'matrix' wide2long(out, keep = 1, na.rm = FALSE)
out |
data.frame or matrix or list of matrices in wide format |
keep |
Index vector, the columns to keep |
na.rm |
Logical, if |
The function assumes that out[,1] represents a time-like vector whereas out[,-1] represents the values. Useful for plotting with ggplot. If a list is supplied, the names of the list are added as extra column names "condition"
data.frame in long format, i.e. columns "time" (out[,1]), "name" (colnames(out[,-1])), "value" (out[,-1]) and, if out was a list, "condition" (names(out))
Write equation list into a csv file
write.eqnlist(eqnlist, ...)write.eqnlist(eqnlist, ...)
eqnlist |
object of class eqnlist |
... |
Arguments going to write.table |
Model prediction function from data.frame
Xd(data, condition = NULL)Xd(data, condition = NULL)
data |
data.frame with columns "name", "time", and row names that are taken as parameter names. The data frame can contain a column "value" to initialize the parameters. |
condition |
either NULL (generic prediction for any condition) or a character, denoting the condition for which the function makes a prediction. |
Object of class prdfn, i.e.
a function x(times pars, deriv = TRUE, conditions = NULL),
see also Xs. Attributes are "parameters", the parameter names (row names of
the data frame), and possibly "pouter", a named numeric vector which is generated
from data$value.
# Generate a data.frame and corresponding prediction function timesD <- seq(0, 2*pi, 0.5) mydata <- data.frame(name = "A", time = timesD, value = sin(timesD), row.names = paste0("par", 1:length(timesD))) x <- Xd(mydata) # Evaluate the prediction function at different time points times <- seq(0, 2*pi, 0.01) pouter <- structure(mydata$value, names = rownames(mydata)) prediction <- x(times, pouter) plot(prediction)# Generate a data.frame and corresponding prediction function timesD <- seq(0, 2*pi, 0.5) mydata <- data.frame(name = "A", time = timesD, value = sin(timesD), row.names = paste0("par", 1:length(timesD))) x <- Xd(mydata) # Evaluate the prediction function at different time points times <- seq(0, 2*pi, 0.01) pouter <- structure(mydata$value, names = rownames(mydata)) prediction <- x(times, pouter) plot(prediction)
Interface to get an ODE
into a model function x(times, pars, forcings, events) returning ODE output.
It is a reduced version of Xs, missing the sensitivities.
Xf( odemodel, forcings = NULL, events = NULL, condition = NULL, optionsOde = list(method = "lsoda"), fcontrol = NULL )Xf( odemodel, forcings = NULL, events = NULL, condition = NULL, optionsOde = list(method = "lsoda"), fcontrol = NULL )
odemodel |
Object of class odemodel. |
forcings |
see Xs |
events |
see Xs |
condition |
either NULL (generic prediction for any condition) or a character, denoting the condition for which the function makes a prediction. |
optionsOde |
list with arguments to be passed to odeC() for the ODE integration. |
fcontrol |
list with additional fine-tuning arguments for the forcing interpolation. See approxfun for possible arguments. |
Can be used to integrate additional quantities, e.g. fluxes, by adding them to f.
All quantities that are not initialised by pars
in x(..., forcings, events) are initialized with 0. For more details and
the return value see Xs.
Interface to combine an ODE and its sensitivity equations
into one model function x(times, pars, deriv = TRUE) returning ODE output and sensitivities.
Xs( odemodel, forcings = NULL, events = NULL, names = NULL, condition = NULL, optionsOde = list(method = "lsoda"), optionsSens = list(method = "lsodes"), fcontrol = NULL )Xs( odemodel, forcings = NULL, events = NULL, names = NULL, condition = NULL, optionsOde = list(method = "lsoda"), optionsSens = list(method = "lsodes"), fcontrol = NULL )
odemodel |
object of class odemodel |
forcings |
data.frame with columns name (factor), time (numeric) and value (numeric). The ODE forcings. |
events |
data.frame of events with columns "var" (character, the name of the state to be
affected), "time" (numeric, time point), "value" (numeric, value), "method" (character, either
"replace", "add" or "multiply"). See events.
ATTENTION: Sensitivities for event states will only be correctly computed if defined within
|
names |
character vector with the states to be returned. If NULL, all states are returned. |
condition |
either NULL (generic prediction for any condition) or a character, denoting the condition for which the function makes a prediction. |
optionsOde |
list with arguments to be passed to odeC() for the ODE integration. |
optionsSens |
list with arguments to be passed to odeC() for integration of the extended system |
fcontrol |
list with additional fine-tuning arguments for the forcing interpolation. See approxfun for possible arguments. |
Object of class prdfn. If the function is called with parameters that result from a parameter transformation (see P), the Jacobian of the parameter transformation and the sensitivities of the ODE are multiplied according to the chain rule for differentiation. The result is saved in the attributed "deriv", i.e. in this case the attibutes "deriv" and "sensitivities" do not coincide.
Function to deal with non-ODE models within the framework of dMod. See example.
Xt(condition = NULL)Xt(condition = NULL)
condition |
either NULL (generic prediction for any condition) or a character, denoting the condition for which the function makes a prediction. |
Object of class prdfn.
x <- Xt() g <- Y(c(y = "a*time^2+b"), f = NULL, parameters = c("a", "b")) times <- seq(-1, 1, by = .05) pars <- c(a = .1, b = 1) plot((g*x)(times, pars))x <- Xt() g <- Y(c(y = "a*time^2+b"), f = NULL, parameters = c("a", "b")) times <- seq(-1, 1, by = .05) pars <- c(a = .1, b = 1) plot((g*x)(times, pars))
Creates an object of type obsfn that evaluates an observation function and its derivatives based on the output of a model prediction function, see prdfn, as e.g. produced by Xs.
Y( g, f = NULL, states = NULL, parameters = NULL, condition = NULL, attach.input = TRUE, deriv = TRUE, compile = FALSE, modelname = NULL, verbose = FALSE )Y( g, f = NULL, states = NULL, parameters = NULL, condition = NULL, attach.input = TRUE, deriv = TRUE, compile = FALSE, modelname = NULL, verbose = FALSE )
g |
Named character vector or equation vector defining the observation function |
f |
Named character of equations or object that can be converted to eqnvec or object of class fn. If f is provided, states and parameters are guessed from f. |
states |
character vector, alternative definition of "states", usually the names of |
parameters |
character vector, alternative definition of the "parameters",
usually the symbols contained in "g" and "f" except for |
condition |
either NULL (generic prediction for any condition) or a character, denoting the condition for which the function makes a prediction. |
attach.input |
logical, indiating whether the original input should be returned with the output. |
deriv |
logical, generate function to evaluate derivatives of observables. Necessary for parameter estimation. |
compile |
Logical, compile the function (see funC0) |
modelname |
Character, used if |
verbose |
Print compiler output to R command line. |
For odemodels with forcings, it is best, to pass the prediction function x to the "f"-argument
instead of the equations themselves. If an eqnvec is passed to "f" in this case, the forcings and states
have to be specified manually via the "states"-argument.
Object of class obsfn, i.e.
a function y(..., deriv = TRUE, conditions = NULL) representing the evaluation of the
observation function. Arguments out (model prediction) and pars (parameter values)
shoudl be passed by the ... argument.
If out has the attribute "sensitivities", the result of
y(out, pars), will have an attributed "deriv" which reflecs the sensitivities of
the observation with respect to the parameters.
If pars is the result of a parameter transformation p(pars) (see P),
the Jacobian
of the parameter transformation and the sensitivities of the observation function
are multiplied according to the chain rule for differentiation.
# Define a time grid on which to make a prediction by peace-wise linear function. # Then define a (generic) prediction function based on thid grid. times <- 0:5 grid <- data.frame(name = "A", time = times, row.names = paste0("p", times)) x <- Xd(grid) # Define an observable and an observation function observables <- eqnvec(Aobs = "s*A") g <- Y(g = observables, f = NULL, states = "A", parameters = "s") # Collect parameters and define an overarching parameter transformation # for two "experimental condtions". dynpars <- attr(x, "parameters") obspars <- attr(g, "parameters") innerpars <- c(dynpars, obspars) trafo <- structure(innerpars, names = innerpars) trafo_C1 <- replaceSymbols(innerpars, paste(innerpars, "C1", sep = "_"), trafo) trafo_C2 <- replaceSymbols(innerpars, paste(innerpars, "C2", sep = "_"), trafo) p <- NULL p <- p + P(trafo = trafo_C1, condition = "C1") p <- p + P(trafo = trafo_C2, condition = "C2") # Collect outer (overarching) parameters and # initialize with random values outerpars <- attr(p, "parameters") pars <- structure(runif(length(outerpars), 0, 1), names = outerpars) # Predict internal/unobserved states out1 <- (x*p)(times, pars) plot(out1) # Predict observed states in addition to unobserved out2 <- (g*x*p)(times, pars) plot(out2)# Define a time grid on which to make a prediction by peace-wise linear function. # Then define a (generic) prediction function based on thid grid. times <- 0:5 grid <- data.frame(name = "A", time = times, row.names = paste0("p", times)) x <- Xd(grid) # Define an observable and an observation function observables <- eqnvec(Aobs = "s*A") g <- Y(g = observables, f = NULL, states = "A", parameters = "s") # Collect parameters and define an overarching parameter transformation # for two "experimental condtions". dynpars <- attr(x, "parameters") obspars <- attr(g, "parameters") innerpars <- c(dynpars, obspars) trafo <- structure(innerpars, names = innerpars) trafo_C1 <- replaceSymbols(innerpars, paste(innerpars, "C1", sep = "_"), trafo) trafo_C2 <- replaceSymbols(innerpars, paste(innerpars, "C2", sep = "_"), trafo) p <- NULL p <- p + P(trafo = trafo_C1, condition = "C1") p <- p + P(trafo = trafo_C2, condition = "C2") # Collect outer (overarching) parameters and # initialize with random values outerpars <- attr(p, "parameters") pars <- structure(runif(length(outerpars), 0, 1), names = outerpars) # Predict internal/unobserved states out1 <- (x*p)(times, pars) plot(out1) # Predict observed states in addition to unobserved out2 <- (g*x*p)(times, pars) plot(out2)