From ea197b16fe5c582dd2a72e81d25c1ebbd5d527b3 Mon Sep 17 00:00:00 2001
From: jranke <jranke@edb9625f-4e0d-4859-8d74-9fd3b1da38cb>
Date: Mon, 23 Apr 2012 23:08:19 +0000
Subject: - Fixed a couple of things - Now the eigenvalue based solutions are
 nicely consistent with the deSolve solutions, if enough   output times are
 specified (100, sometimes more are needed, see test.R) - Workaround for
 invilr not to produce NaN values so often - Still a lot to do (see TODO)

git-svn-id: svn+ssh://svn.r-forge.r-project.org/svnroot/kinfit/pkg/mkin@30 edb9625f-4e0d-4859-8d74-9fd3b1da38cb
---
 R/ilr.R         | 7 +++++++
 R/mkinfit.R     | 6 +++++-
 R/mkinmod.R     | 6 +++++-
 R/mkinpredict.R | 4 +++-
 4 files changed, 20 insertions(+), 3 deletions(-)

(limited to 'R')

diff --git a/R/ilr.R b/R/ilr.R
index b53f8683..389653e7 100644
--- a/R/ilr.R
+++ b/R/ilr.R
@@ -40,5 +40,12 @@ invilr<-function(x) {
   for (i in 1:D) {
     z[i] <- exp(y[i])/sum(exp(y))
   }
+
+  # Work around a numerical problem with NaN values returned by the above
+  # Only works if there is only one NaN value: replace it with 1
+  # if the sum of the other components is < 1e-10
+  if (sum(is.na(z)) == 1 && sum(z[!is.na(z)]) < 1e-10) 
+    z = ifelse(is.na(z), 1, z)
+
   return(z)
 }
diff --git a/R/mkinfit.R b/R/mkinfit.R
index 37eee330..cf58a7d2 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -111,6 +111,9 @@ mkinfit <- function(mkinmod, observed,
   {
     assign("calls", calls+1, inherits=TRUE) # Increase the model solution counter
 
+    # Trace parameter values if quiet is off
+    if(!quiet) cat(P, "\n")
+
     # Time points at which observed data are available
     # Make sure we include time 0, so initial values for state variables are for time 0
     outtimes = sort(unique(c(observed$time,
@@ -305,7 +308,8 @@ mkinfit <- function(mkinmod, observed,
   data$variable <- ordered(data$variable, levels = obs_vars)
   fit$data <- data[order(data$variable, data$time), ]
   fit$atol <- atol
-  fit$parms.all <- parms.all
+  fit$parms.all <- parms.all # Return all backtransformed parameters for summary
+  fit$odeparms.final <- parms.all[mkinmod$parms] # Return ode parameters for further fitting
 
   class(fit) <- c("mkinfit", "modFit") 
   return(fit)
diff --git a/R/mkinmod.R b/R/mkinmod.R
index 54b3b5f4..41dbefa8 100644
--- a/R/mkinmod.R
+++ b/R/mkinmod.R
@@ -92,7 +92,7 @@ mkinmod <- function(..., use_of_ff = "min")
           parms <- c(parms, k_compound_sink)
           decline_term <- paste(k_compound_sink, "*", box_1)
         } else { # otherwise no decline term needed here
-          decline_term = "" 
+          decline_term = "0" 
         }
       } else {
         k_compound <- paste("k", box_1, sep="_")
@@ -128,6 +128,10 @@ mkinmod <- function(..., use_of_ff = "min")
         reversible_binding_term_2 <- paste("+", k_free_bound, "*", box_1, "-",
           k_bound_free, "*", box_2)
       } else { # Use formation fractions also for the free compartment
+	stop("The maximum use of formation fractions is not supported for SFORB models")
+        # The problems were: Calculation of dissipation times did not work in this case
+        # and the coefficient matrix is not generated correctly by the code present 
+        # in this file in this case
         f_free_bound <- paste("f", varname, "free", "bound", sep="_")
         k_bound_free <- paste("k", varname, "bound", "free", sep="_")
         parms <- c(parms, f_free_bound, k_bound_free)
diff --git a/R/mkinpredict.R b/R/mkinpredict.R
index 49166f30..2b7e51d5 100644
--- a/R/mkinpredict.R
+++ b/R/mkinpredict.R
@@ -17,7 +17,9 @@ mkinpredict <- function(mkinmod, odeparms, odeini, outtimes, solution_type = "de
     o <- switch(parent.type,
       SFO = SFO.solution(outtimes, 
           evalparse(parent.name),
-          evalparse(paste("k", parent.name, sep="_"))),
+          ifelse(mkinmod$use_of_ff == "min", 
+	    evalparse(paste("k", parent.name, "sink", sep="_")),
+	    evalparse(paste("k", parent.name, sep="_")))),
       FOMC = FOMC.solution(outtimes,
           evalparse(parent.name),
           evalparse("alpha"), evalparse("beta")),
-- 
cgit v1.2.1