## Not run:
#
# ################################################################
# ## Walk through the most frequently used functions in
# ## connection with ODE models
# ################################################################
#
# library(deSolve)
#
# ## Model definition (text-based, scripting part)
# r <- NULL
# r <- addReaction(r, "A", "B", "k1*A", "translation")
# r <- addReaction(r, "B", "", "k2*B", "degradation")
#
# f <- as.eqnvec(r)
# observables <- eqnvec(Bobs = "s1*B")
#
# innerpars <- getSymbols(c(names(f), f, observables))
# trafo0 <- structure(innerpars, names = innerpars)
# trafo0 <- replaceSymbols(innerpars,
# paste0("exp(", innerpars, ")"),
# trafo0)
#
# conditions <- c("a", "b")
#
# trafo <- list()
# trafo$a <- replaceSymbols("s1", "s1_a", trafo0)
# trafo$b <- replaceSymbols("s1", "s1_b", trafo0)
#
# events <- list()
# events$a <- data.frame(var = "A", time = 5, value = 1, method = "add")
# events$b <- data.frame(var = "A", time = 5, value = 2, method = "add")
#
#
# ## Model definition (generate compiled objects
# ## and functions from above information)
#
# # ODE model
# <<<<<<< HEAD
# model <- odemodel(f, jacobian = "inz.lsodes")
# =======
# model <- odemodel(f)
# >>>>>>> development
#
# # Observation function
# g <- Y(observables, f, compile = TRUE, modelname = "obsfn")
#
# # Parameter transformation for steady states
# pSS <- P(f, method = "implicit", condition = NULL)
#
# # Condition-specific transformation and prediction functions
# p0 <- x <- NULL
# for (C in conditions) {
# p0 <- p0 + P(trafo[[C]], condition = C)
# x <- x + Xs(model, events = events[[C]], condition = C)
# }
#
# ## 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))
# )
#
# ## Evaluate model/data
#
# # Initialize times and parameters
# times <- seq(0, 10, .1)
# outerpars <- attr(p0, "parameters")
# pouter <- structure(rnorm(length(outerpars)),
# names = outerpars)
#
# # Plot prediction
# plot((g*x*p0)(times, pouter))
# plot((g*x*pSS*p0)(times, pouter))
# plotFluxes(pouter, g*x*p0, times, getFluxes(r)$B)
#
# # Fit data with L2 constraint
# constr <- constraintL2(mu = 0*pouter, sigma = 5)
# myfit <- trust(normL2(data, g*x*p0) + constr,
# pouter, rinit = 1, rmax = 10)
# plot((g*x*p0)(times, myfit$argument), data)
#
# # List of fits
# obj <- normL2(data, g*x*p0) + constr
# mylist <- mstrust(normL2(data, g*x*p0) + constr,
# pouter, fits = 10, cores = 4, sd = 4)
# plot(mylist)
# pars <- as.parframe(mylist)
# plotValues(pars)
#
# bestfit <- as.parvec(pars)
#
# # Compute profiles of the fit
# profile <- profile(normL2(data, g*x*p0) + constr,
# bestfit, "k1", limits = c(-5, 5))
# plotProfile(profile)
#
# # Add a validation objective function
# vali <- datapointL2("A", 2, "mypoint", .1, condition = "a")
# # Get validation point parameters
# mypoint <- trust(normL2(data, g*x*p0) + constr + vali,
# c(mypoint = 0), rinit = 1, rmax = 10,
# fixed = bestfit)$argument
# # Compoute profile and plot
# profile <- profile(normL2(data, g*x*p0) + constr + vali,
# c(bestfit, mypoint), "mypoint")
# plotProfile(profile)
#
#
# ## Using the controls() function to manipulate objects ------------------------
#
# # List the available controls of an object
# controls(x)
#
# # List a specific control
# controls(x, condition = "a", name = "optionsSens")
#
# # Change specific control
# controls(x, condition = "a", name = "optionsSens") <-
# list(method = "lsoda", rtol = 1e-8, atol = 1e-8)
#
# # Almost every function (objfn, parfn, prdfn, obsfn) saves the arguments
# # by which it was invoked in a list named "controls", which can be manipulated
# # by the controls() function.
#
#
# ## New example: condition-specific observation parameters ---------------------
#
# # Condition-specific observation parameters
# f <- eqnvec(A = "-k1*A", B = "k1*A - k2*B")
# observables <- eqnvec(Bobs = "s1*B")
# conditions <- c("scale1", "scale2")
#
# dynpars <- getSymbols(c(names(f), f))
# obspars <- setdiff(getSymbols(observables), dynpars)
# trafo0 <- structure(c(dynpars, obspars), names = c(dynpars, obspars))
#
# obspars_all <- outer(obspars, conditions, paste, sep = "_")
# trafo_all <- structure(c(dynpars, obspars_all),
# names = c(dynpars, obspars_all))
# trafo_all <- replaceSymbols(dynpars,
# paste0("exp(", dynpars, ")"),
# trafo_all)
# trafo_all <- replaceSymbols(obspars_all,
# paste0("exp(", obspars_all, ")"),
# trafo_all)
#
# pObs <- NULL
# for (C in conditions) {
# pObs <- pObs + P(replaceSymbols(obspars,
# paste(obspars, C, sep = "_"),
# trafo0),
# condition = C)
# }
#
#
# x <- Xs(model)
# p <- P(trafo_all)
# y <- g*pObs*x*p
#
# times <- seq(0, 10, .1)
# outerpars <- attr(p, "parameters")
# pouter <- structure(rnorm(length(outerpars)), names = outerpars)
#
# # Plot prediction
# plot((g*pObs*x*p)(times, pouter))
#
#
# ## End(Not run)
Run the code above in your browser using DataLab