#$Id$
# Generates fitted curves so the GUI can plot them
# Based on code in IRLSkinfit
# Author: Rob Nelson (Tessella)
# Modifications developed by Tessella Plc for Syngenta: Copyright (C) 2011  Syngenta
# 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/>.”

CakePlotInit <- function(fit, xlim = range(fit$data$time), ...)
{
    t.map.names <- names(fit$map)
    metabolites <- grep("[A-Z]\\d", t.map.names, value = TRUE)
    t.map.rest  <- setdiff(t.map.names, metabolites)
    
    # Generate the normal graphs.
    final <- CakeOlsPlot(fit, xlim)
    final_scaled <- final
    
    if(length(metabolites) > 0){
        for(i in 1:length(metabolites))
        {
            metabolite <- metabolites[i]
            decay_var  <- paste("k", metabolite, sep="_")
            
            # calculate the new ffm and generate the two ffm scale charts
            regex <- paste("f_(.+)_to", metabolite, sep="_")
            decays = grep(regex, names(fit$par), value = TRUE)
            ffm_fitted <- sum(fit$par[decays])
            normal <- final
            ffm_scale <- NULL
                
            # Generate the DT50=1000d and ffm as fitted.
            k_new <- fit$par[decay_var]*fit$distimes[metabolite,]["DT50"]/1000;
            fit$par[decay_var]<- k_new[metabolite,]
            dt50_1000_ffm_fitted <- CakeOlsPlot(fit, xlim)[metabolite]
            
            naming <- c(names(final), paste(metabolite, "DT50_1000_FFM_FITTED", sep="_"))
            normal <- c(final, dt50_1000_ffm_fitted)
            names(normal) <- naming
            final <- normal
            
            # Generate the scaled FFM
            if(ffm_fitted != 0)
            {
                ffm_scale <- 1 / ffm_fitted
                final_scaled <- final[c("time", metabolite, paste(metabolite, "DT50_1000_FFM_FITTED", sep="_"))]
                final_scaled[t.map.rest] <- NULL;
                final_frame <- as.data.frame(final_scaled)
                new_names <- c(paste(metabolite, "DT50_FITTED_FFM_1", sep="_"), paste(metabolite, "DT50_1000_FFM_1", sep="_"))
                names(final_frame) <- c("time", new_names)
                final_frame[new_names]<-final_frame[new_names]*ffm_scale;
                
                cat("<PLOT MODEL START>\n")
            
                write.table(final_frame, quote=FALSE)
                
                cat("<PLOT MODEL END>\n")
            }
        }
    }
    
    cat("<PLOT MODEL START>\n")
    
    write.table(final, quote=FALSE)
    
    cat("<PLOT MODEL END>\n")

    # View(final)
}

CakeOlsPlot <- function(fit, xlim = range(fit$data$time), scale_x = 1.0, ...)
{
   solution = fit$solution
   if ( is.null(solution) ) {
      solution <- "deSolve"
   }
   atol = fit$atol
   if ( is.null(atol) ) {
      atol = 1.0e-6
   }
   
  fixed <- fit$fixed$value
  names(fixed) <- rownames(fit$fixed)
  parms.all <- c(fit$par, fixed)
  ininames <- c(
    rownames(subset(fit$start, type == "state")),
    rownames(subset(fit$fixed, type == "state")))
  odeini <- parms.all[ininames]
  names(odeini) <- names(fit$diffs)

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

  odenames <- c(
    rownames(subset(fit$start, type == "deparm")),
    rownames(subset(fit$fixed, type == "deparm")))
  odeparms <- parms.all[odenames]

  # Solve the system
  evalparse <- function(string)
  {
    eval(parse(text=string), as.list(c(odeparms, odeini)))
  }
  if (solution == "analytical") {
    parent.type = names(fit$map[[1]])[1]  
    parent.name = names(fit$diffs)[[1]]
    o <- switch(parent.type,
      SFO = SFO.solution(outtimes, 
          evalparse(parent.name),
          evalparse(paste("k", parent.name, sep="_"))),
      FOMC = FOMC.solution(outtimes,
          evalparse(parent.name),
          evalparse("alpha"), evalparse("beta")),
      DFOP = DFOP.solution(outtimes,
          evalparse(parent.name),
          evalparse("k1"), evalparse("k2"),
          evalparse("g")),
      HS = HS.solution(outtimes,
          evalparse(parent.name),
          evalparse("k1"), evalparse("k2"),
          evalparse("tb")),
      SFORB = SFORB.solution(outtimes,
          evalparse(parent.name),
          evalparse(paste("k", parent.name, "free_bound", sep="_")),
          evalparse(paste("k", parent.name, "bound_free", sep="_")),
          evalparse(paste("k", parent.name, "free_sink", sep="_")))
    )
    out <- cbind(outtimes, o)
    dimnames(out) <- list(outtimes, c("time", parent.name))
  }
  if (solution == "eigen") {
    coefmat.num <- matrix(sapply(as.vector(fit$coefmat), evalparse), 
      nrow = length(odeini))
    e <- eigen(coefmat.num)
    c <- solve(e$vectors, odeini)
    f.out <- function(t) {
      e$vectors %*% diag(exp(e$values * t), nrow=length(odeini)) %*% c
    }
    o <- matrix(mapply(f.out, outtimes), 
      nrow = length(odeini), ncol = length(outtimes))
    dimnames(o) <- list(names(odeini), NULL)
    out <- cbind(time = outtimes, t(o))
  } 
  if (solution == "deSolve") {
    out <- ode(
      y = odeini,
      times = outtimes,
      func = fit$mkindiff, 
      parms = odeparms,
      atol = atol
    )
  }
    
  # Output transformation for models with unobserved compartments like SFORB
  out_transformed <- data.frame(time = out[,"time"])
  for (var in names(fit$map)) {
    if(length(fit$map[[var]]) == 1) {
      out_transformed[var] <- out[, var]
    } else {
      out_transformed[var] <- rowSums(out[, fit$map[[var]]])
    }
  }
  return(out_transformed)
}