1 files changed, 66 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), ...)
+}