From 30cbb4092f6d2d3beff5800603374a0d009ad770 Mon Sep 17 00:00:00 2001
From: jranke <jranke@edb9625f-4e0d-4859-8d74-9fd3b1da38cb>
Date: Tue, 11 May 2010 23:03:37 +0000
Subject: 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
---
 R/mkinfit.R | 66 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 66 insertions(+)
 create mode 100644 R/mkinfit.R

(limited to 'R/mkinfit.R')

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), ...)
+}
-- 
cgit v1.2.1