aboutsummaryrefslogtreecommitdiff
path: root/R/mkinfit.R
diff options
context:
space:
mode:
authorjranke <jranke@edb9625f-4e0d-4859-8d74-9fd3b1da38cb>2010-05-11 23:03:37 +0000
committerjranke <jranke@edb9625f-4e0d-4859-8d74-9fd3b1da38cb>2010-05-11 23:03:37 +0000
commit30cbb4092f6d2d3beff5800603374a0d009ad770 (patch)
treeef75421d92823b5b7add1b2d5da9c7499dfd59e5 /R/mkinfit.R
Initial upload of the upcoming multicompartment version of kinfit.
Some functionality is still missing (chi2), some may never be implemented (FOMC model), but in general it is much more powerful than kinfit, owing to the powerful FME package. git-svn-id: svn+ssh://svn.r-forge.r-project.org/svnroot/kinfit/pkg/mkin@8 edb9625f-4e0d-4859-8d74-9fd3b1da38cb
Diffstat (limited to 'R/mkinfit.R')
-rw-r--r--R/mkinfit.R66
1 files changed, 66 insertions, 0 deletions
diff --git a/R/mkinfit.R b/R/mkinfit.R
new file mode 100644
index 0000000..9651fd6
--- /dev/null
+++ b/R/mkinfit.R
@@ -0,0 +1,66 @@
+mkinfit <- function(mkinmod, observed,
+ parms.ini = rep(0.1, length(mkinmod$parms)),
+ state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)),
+ fixed_parms = rep(FALSE, length(mkinmod$parms)),
+ fixed_initials = c(FALSE, rep(TRUE, length(mkinmod$diffs) - 1)),
+ plot = NULL,
+ err = NULL, weight = "none", scaleVar = FALSE,
+ ...)
+{
+ # Name the parameters if they are not named yet
+ if(is.null(names(parms.ini))) names(parms.ini) <- mkinmod$parms
+ # Create a function calculating the differentials specified by the model
+ mkindiff <- function(t, state, parms) {
+ diffs <- vector()
+ for (box in names(mkinmod$diffs))
+ {
+ diffname <- paste("d", box, sep="_")
+ diffs[diffname] <- with(as.list(c(state, parms)),
+ eval(parse(text=mkinmod$diffs[[box]])))
+ }
+ return(list(c(diffs)))
+ }
+
+ # Name the inital parameter values if they are not named yet
+ if(is.null(names(state.ini))) names(state.ini) <- names(mkinmod$diffs)
+
+ # TODO: Collect parameters to be optimised
+ parms.optim <- parms.ini[!fixed_parms]
+ parms.fixed <- parms.ini[fixed_parms]
+
+ state.ini.optim <- state.ini[!fixed_initials]
+ state.ini.optim.boxnames <- names(state.ini.optim)
+ names(state.ini.optim) <- paste(names(state.ini.optim), "0", sep="_")
+ state.ini.fixed <- state.ini[fixed_initials]
+
+ # Define the model cost function
+ cost <- function(P)
+ {
+ if(length(state.ini.optim) > 0) {
+ odeini <- c(P[1:length(state.ini.optim)], state.ini.fixed)
+ names(odeini) <- c(state.ini.optim.boxnames, names(state.ini.fixed))
+ } else odeini <- state.ini.fixed
+
+ odeparms <- c(P[(length(state.ini.optim) + 1):length(P)], parms.fixed)
+ # Solve the ODE
+ out <- ode(
+ y = odeini,
+ times = unique(observed$time),
+ func = mkindiff,
+ parms = odeparms)
+
+ # Output transformation for models with ghost compartments like SFORB
+ out_transformed <- data.frame(time = out[,"time"])
+ for (var in names(mkinmod$map)) {
+ if(length(mkinmod$map[[var]]) == 1) {
+ out_transformed[var] <- out[, var]
+ } else {
+ out_transformed[var] <- rowSums(out[, mkinmod$map[[var]]])
+ }
+ }
+
+ return(modCost(out_transformed, observed, y = "value",
+ err = err, weight = weight, scaleVar = scaleVar))
+ }
+ modFit(cost, c(state.ini.optim, parms.optim), ...)
+}

Contact - Imprint