path: root/CakeNullPlot.R

#$Id$
# Generates model curves so the GUI can plot them. Used when all params are fixed so there was no fit
# The CAKE R modules are based on mkin
# Based on code in IRLSkinfit
# Modifications developed by Tessella Plc for Syngenta: Copyright (C) 2011  Syngenta
# Author: Rob Nelson (Tessella)
# Tessella Project Reference: 6245

#    This program is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
# 
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
# 
#    You should have received a copy of the GNU General Public License
#    along with this program.  If not, see <http://www.gnu.org/licenses/>.”

CakeNullPlot <- function(model, parms.ini, state.ini, observed, xlim = range(observed$time), digits = max(3, getOption("digits") - 3), sannMaxIter = 10000, ...)
{
  # Print the fixed values
  fixed <- data.frame(value = c(state.ini, parms.ini))
  fixed$type = c(rep("state", length(state.ini)), rep("deparm", length(parms.ini)))
  cat("\nFixed parameter values:\n")
  print(fixed)

  # Print disappearence times
  obs_vars = unique(as.character(observed$name))
  distimes <- data.frame(DT50 = rep(NA, length(obs_vars)), 
        DT90 = rep(NA, length(obs_vars)), row.names = obs_vars)
  extraDT50<- data.frame(DT50 = rep(NA, 2), row.names = c("k1", "k2"))
  
  for (obs_var in obs_vars) {
      type = names(model$map[[obs_var]])[1]
      if(exists("DT50")) distimes[obs_var, ] = c(DT50, DT90)
      distimes[obs_var, ] = CakeDT(type,obs_var,parms.ini,sannMaxIter)
   }

   cat("\nEstimated disappearance times:\n")
   print(distimes, digits=digits,...)
   
   extraDT50[ ,c("DT50")] = CakeExtraDT(type,parms.ini)
   cat("\nHalf Lives for Extra k values:\n")
   print(extraDT50, digits=digits,...)

   # Plot the model

  outtimes <- seq(0, xlim[2], length.out=101)

  odeini <- state.ini
  odenames <- model$parms
  odeparms <- parms.ini

  # Solve the system
  evalparse <- function(string)
  {
    eval(parse(text=string), as.list(c(odeparms, odeini)))
  }
  
    mod_vars <- names(model$diffs)
    mkindiff <- function(t, state, parms) {
        time <- t
        diffs <- vector()
        for (box in mod_vars) {
            diffname <- paste("d", box, sep = "_")
            diffs[diffname] <- with(as.list(c(time, state, parms)), 
                eval(parse(text = model$diffs[[box]])))
        }
        return(list(c(diffs)))
    }
  
   # Ensure initial state is at time 0
   obstimes <- unique(c(0,observed$time))
    out <- ode(
      y = odeini,
      times = obstimes,
      func = mkindiff, 
      parms = odeparms,
    )
    
  # Output transformation for models with unobserved compartments like SFORB
  out_transformed <- data.frame(time = out[,"time"])
  for (var in names(model$map)) {
    if(length(model$map[[var]]) == 1) {
      out_transformed[var] <- out[, var]
    } else {
      out_transformed[var] <- rowSums(out[, model$map[[var]]])
    }
  }    
  
   predicted_long <- mkin_wide_to_long(out_transformed,time="time")
  # Calculate chi2 error levels according to FOCUS (2006)
   errmin<-CakeChi2(model, observed, predicted_long, obs_vars, c("auto"), c(), state.ini, c("auto"))
   cat("\nChi2 error levels in percent:\n")
   errmin$err.min <- 100 * errmin$err.min
   print(errmin, digits=digits,...)

  
   # Fit to data
   data<-observed
   data$err<-rep(NA,length(data$time))
   data<-merge(data, predicted_long, by=c("time","name"))
   names(data)<-c("time", "variable", "observed","err-var", "predicted")
   data$residual<-data$observed-data$predicted
   data<-data[order(data$variable,data$time),]
   cat("\nData:\n")
   print(format(data, digits = digits, scientific = FALSE,...), row.names = FALSE)
   
  
    out <- ode(
      y = odeini,
      times = outtimes,
      func = mkindiff, 
      parms = odeparms,
    )
    
  # Output transformation for models with unobserved compartments like SFORB
  out_transformed <- data.frame(time = out[,"time"])
  for (var in names(model$map)) {
    if(length(model$map[[var]]) == 1) {
      out_transformed[var] <- out[, var]
    } else {
      out_transformed[var] <- rowSums(out[, model$map[[var]]])
    }
  }    
   
   cat("\n<PLOT MODEL START>\n")
   print(out_transformed)
   cat("<PLOT MODEL END>\n")
}