From ea197b16fe5c582dd2a72e81d25c1ebbd5d527b3 Mon Sep 17 00:00:00 2001 From: jranke 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 +++- TODO | 1 + 5 files changed, 21 insertions(+), 3 deletions(-) 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")), diff --git a/TODO b/TODO index 10ccf914..13bc6fce 100644 --- a/TODO +++ b/TODO @@ -1,4 +1,5 @@ - Fix transformation and backtransformation of formation fractions +- Fix problems with the schaefer complex test data (test.R) - Adapt mkinplot function - Calculate confidence intervals for parameters assuming normal distribution - Calculate confidence intervals for DT50 and DT90 values when only one parameter is involved -- cgit v1.2.1