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
# Some of the CAKE R modules are based on mkin
# Based on code in IRLSkinfit
# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2016 Syngenta
# Tessella Project Reference: 6245, 7247, 8361, 7414

#    The CAKE R modules are 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, atol = 1e-6, ...)
{
  # 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)
  parms.all = c(state.ini, parms.ini)

  # Print disappearance 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(k1 = rep(NA, length(obs_vars)), k2 = rep(NA, length(obs_vars)), row.names = obs_vars)
  
  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.all,sannMaxIter)
    extraDT50[obs_var, ] = CakeExtraDT(type, obs_var, parms.all)
  }
  
  cat("\nEstimated disappearance times:\n")
  print(distimes, digits=digits,...)
   
  cat("\nHalf Lives for Extra k values:\n")
  print(extraDT50, digits=digits,...)
   
  ioreRepDT = CakeIORERepresentativeDT("Parent", parms.all)
  printIoreRepDT <- !is.null(ioreRepDT)
  if (printIoreRepDT){
    cat("\nRepresentative Half Life for IORE:", format(signif(ioreRepDT, digits)))
  }
   
  cat("\n")
  
  fomcRepDT = CakeFOMCBackCalculatedDT(parms.all)
  printFomcRepDT <- !is.null(fomcRepDT)
  
  if (printIoreRepDT){
    cat("\nBack-calculated Half Life for FOMC:", format(signif(fomcRepDT, digits)))
  }
  
  cat("\n")
  
  # Plot the model

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

  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))
  
  odeini <- AdjustOdeInitialValues(odeini, model, odeparms)
  
  # Solve at observation times
  out <- ode(y = odeini, times = obstimes, func = mkindiff, parms = odeparms, atol = atol)
    
  out_transformed <- PostProcessOdeOutput(out, model, atol)
  
  predicted_long <- wide_to_long(out_transformed,time="time")
  
  # Calculate chi2 error levels according to FOCUS (2006)
  errmin<-CakeChi2(model, observed, predicted_long, obs_vars, c(), c(), state.ini, parms.ini, fixed)
  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)
  
  # Solve at output times
  out <- ode(
    y = odeini,
    times = outtimes,
    func = mkindiff, 
    parms = odeparms)
    
  out_transformed <- PostProcessOdeOutput(out, model, atol)
   
  cat("\n<PLOT MODEL START>\n")
  print(out_transformed)
  cat("<PLOT MODEL END>\n")
}