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

2 files changed, 142 insertions, 0 deletions

diff --git a/R/mkinfit.R b/R/mkinfit.R
new file mode 100644
index 00000000..9651fd66
--- /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), ...)

+}

diff --git a/R/mkinmod.R b/R/mkinmod.R
new file mode 100644
index 00000000..7e64f879
--- /dev/null
+++ b/R/mkinmod.R
@@ -0,0 +1,76 @@
+mkinmod <- function(spec = list(parent = list(type = "SFO", to = NA, sink = TRUE)))

+{

+  if(!is.list(spec)) stop("spec must be a list of model specifications for each observed variable")

+

+  # The returned model will be a list of two character vectors, containing

+  # differential equations and parameter names

+  parms <- vector()

+  diffs <- vector()

+  map <- list()

+

+  # Establish list of differential equations

+  for (varname in names(spec))

+  {

+    if(!spec[[varname]]$type %in% c("SFO", "SFORB")) stop("Available types are SFO and SFORB only")

+    new_parms <- vector()

+

+    # New (sub)compartments (boxes) needed for the model type

+    new_boxes <- switch(spec[[varname]]$type,

+      SFO = varname,

+      SFORB = paste(varname, c("free", "bound"), sep="_")

+    )

+    map[[varname]] <- new_boxes

+

+    # Start a new differential equation for each new box

+    new_diffs <- paste("d_", new_boxes, " =", sep="")

+

+    # Construct terms for transfer to sink and add if appropriate

+    if(spec[[varname]]$sink) {

+      k_compound_sink <- paste("k", new_boxes[[1]], "sink", sep="_")

+      sink_term <- paste("-", k_compound_sink, "*", new_boxes[[1]])

+      # Only add sink term to first (or only) box for SFO and SFORB models

+      if(spec[[varname]]$type %in% c("SFO", "SFORB")) {

+        new_diffs[[1]] <- paste(new_diffs[[1]], sink_term)

+      }

+      new_parms <- k_compound_sink

+    }

+   

+    # Add reversible binding if appropriate

+    if(spec[[varname]]$type == "SFORB") {

+      k_free_bound <- paste("k", varname, "free", "bound", sep="_")      

+      k_bound_free <- paste("k", varname, "bound", "free", sep="_")      

+      reversible_binding_terms <- c(

+        paste("-", k_free_bound, "*", new_boxes[[1]], "+", k_bound_free, "*", new_boxes[[2]]),

+        paste("+", k_free_bound, "*", new_boxes[[1]], "-", k_bound_free, "*", new_boxes[[2]]))

+      new_diffs <- paste(new_diffs, reversible_binding_terms)

+      new_parms <- c(new_parms, k_free_bound, k_bound_free)

+    } 

+

+    # Add observed variable to model

+    parms <- c(parms, new_parms)

+    names(new_diffs) <- new_boxes

+    diffs <- c(diffs, new_diffs)

+  }

+

+  # Transfer between compartments

+  for (varname in names(spec)) {

+    to <- spec[[varname]]$to

+    if(!is.na(to)) {

+      origin_box <- switch(spec[[varname]]$type,

+        SFO = varname,

+        SFORB = paste(varname, "free", sep="_"))

+      for (target in to) {

+        target_box <- switch(spec[[target]]$type,

+          SFO = target,

+          SFORB = paste(target, "free", sep="_"))

+      }

+      k_from_to <- paste("k", origin_box, target_box, sep="_")

+      diffs[[origin_box]] <- paste(diffs[[origin_box]], "-", k_from_to, "*", origin_box)

+      diffs[[target_box]] <- paste(diffs[[target_box]], "+", k_from_to, "*", origin_box)

+      parms <- c(parms, k_from_to)

+    }

+  }

+  model <- list(diffs = diffs, parms = parms, map = map)

+  class(model) <- "mkinmod"

+  return(model)

+}