From e3162e617bc268d9de92640311e2fbe650aa636a Mon Sep 17 00:00:00 2001
From: jranke <jranke@edb9625f-4e0d-4859-8d74-9fd3b1da38cb>
Date: Mon, 23 Apr 2012 06:36:24 +0000
Subject: - Some fixes to mkinfit, to account for the changes in mkinmod - See
 TODO for what remains to be done (not complete)

git-svn-id: svn+ssh://svn.r-forge.r-project.org/svnroot/kinfit/pkg/mkin@29 edb9625f-4e0d-4859-8d74-9fd3b1da38cb
---
 R/mkinfit.R | 43 +++++++++++++++++++++++++------------------
 1 file changed, 25 insertions(+), 18 deletions(-)

(limited to 'R')

diff --git a/R/mkinfit.R b/R/mkinfit.R
index 64edaba..37eee33 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -26,25 +26,22 @@ mkinfit <- function(mkinmod, observed,
   state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)), 
   fixed_parms = NULL,
   fixed_initials = names(mkinmod$diffs)[-1],
-  eigen = FALSE,
+  solution_type = "auto",
   plot = FALSE, quiet = FALSE,
   err = NULL, weight = "none", scaleVar = FALSE,
-  atol = 1e-6,
+  atol = 1e-6, n.outtimes = 100,
   ...)
 {
   # Get the names of the state variables in the model
   mod_vars <- names(mkinmod$diffs)
 
-  # See which variant of the model specification was used
-  use_of_ff <- mkinmod$use_of_ff
-
   # Subset observed data with names of observed data in the model
   observed <- subset(observed, name %in% names(mkinmod$map))
 
   # Get names of observed variables
   obs_vars = unique(as.character(observed$name))
 
-  # Define initial values for parameters where not specified by the user
+  # Define starting values for parameters where not specified by the user
   if (parms.ini[[1]] == "auto") parms.ini = vector()
   defaultpar.names <- setdiff(mkinmod$parms, names(parms.ini))
   for (parmname in defaultpar.names) {
@@ -54,7 +51,7 @@ mkinfit <- function(mkinmod, observed,
     if (grepl("free_bound$", parmname)) parms.ini[parmname] = 0.1 
     if (grepl("bound_free$", parmname)) parms.ini[parmname] = 0.02
     # Default values for formation fractions
-    if (substr(parmname, 1, 2) == "f_") parms.ini[parmname] = 0.1
+    if (substr(parmname, 1, 2) == "f_") parms.ini[parmname] = 0.2
     # Default values for the FOMC, DFOP and HS models
     if (parmname == "alpha") parms.ini[parmname] = 1
     if (parmname == "beta") parms.ini[parmname] = 10
@@ -88,13 +85,21 @@ mkinfit <- function(mkinmod, observed,
   # Decide if the solution of the model can be based on a simple analytical
   # formula, the spectral decomposition of the matrix (fundamental system)
   # or a numeric ode solver from the deSolve package
-  if (length(mkinmod$map) == 1) {
-    solution = "analytical"
-  } else {
-    if (is.matrix(mkinmod$coefmat) && eigen) {
-      solution = "eigen"
+  if (!solution_type %in% c("auto", "analytical", "eigen", "deSolve"))
+      stop("solution_type must be auto, analytical, eigen or de Solve")
+  if (solution_type == "analytical" && length(mkinmod$map) > 1)
+      stop("Analytical solution not implemented for models with metabolites.")
+  if (solution_type == "eigen" && !is.matrix(mkinmod$coefmat))
+      stop("Eigenvalue based solution not possible, coefficient matrix not present.")
+  if (solution_type == "auto") {
+    if (length(mkinmod$map) == 1) {
+      solution_type = "analytical"
     } else {
-      solution = "deSolve"
+      if (is.matrix(mkinmod$coefmat)) {
+	solution_type = "eigen"
+      } else {
+	solution_type = "deSolve"
+      }
     }
   }
 
@@ -107,7 +112,9 @@ mkinfit <- function(mkinmod, observed,
     assign("calls", calls+1, inherits=TRUE) # Increase the model solution counter
 
     # Time points at which observed data are available
-    outtimes = unique(observed$time)
+    # Make sure we include time 0, so initial values for state variables are for time 0
+    outtimes = sort(unique(c(observed$time,
+		      seq(min(observed$time), max(observed$time), length.out=n.outtimes))))
 
     if(length(state.ini.optim) > 0) {
       odeini <- c(P[1:length(state.ini.optim)], state.ini.fixed)
@@ -119,7 +126,7 @@ mkinfit <- function(mkinmod, observed,
     parms <- backtransform_odeparms(odeparms, mod_vars)
 
     # Solve the system with current transformed parameter values
-    out <- mkinpredict(mkinmod, parms, odeini, outtimes, solution_type = solution, ...)
+    out <- mkinpredict(mkinmod, parms, odeini, outtimes, solution_type = solution_type, ...)
 
     assign("out_predicted", out, inherits=TRUE)
 
@@ -135,7 +142,7 @@ mkinfit <- function(mkinmod, observed,
         outtimes_plot = seq(min(observed$time), max(observed$time), length.out=100)
 
         out_plot <- mkinpredict(mkinmod, parms, odeini, outtimes_plot, 
-          solution_type = solution, ...)
+          solution_type = solution_type, ...)
 
         plot(0, type="n", 
           xlim = range(observed$time), ylim = range(observed$value, na.rm=TRUE),
@@ -158,8 +165,8 @@ mkinfit <- function(mkinmod, observed,
   fit <- modFit(cost, c(state.ini.optim, parms.optim), ...)
 
   # We need to return some more data for summary and plotting
-  fit$solution <- solution
-  if (solution == "eigen") {
+  fit$solution_type <- solution_type
+  if (solution_type == "eigen") {
     fit$coefmat <- mkinmod$coefmat
   } 
 
-- 
cgit v1.2.1