#$Id$
# Generates fitted curves so the GUI can plot them
# Based on code in IRLSkinfit
# Modifications developed by Hybrid Intelligence (formerly Tessella), part of
# Capgemini Engineering, for Syngenta, Copyright (C) 2011-2022 Syngenta
# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091
#
#    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/>.

CakeDoPlot <- function(fit, xlim = range(fit$data$time), ...)
{
    original_fit <- fit
    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 <- CakeSinglePlot(fit, xlim)
    
    if (length(metabolites) > 0) {
        for (i in 1:length(metabolites)) {
            metabolite <- metabolites[i]
            
            fit_method = names(fit$map[[metabolite]])[1]
            if (fit_method != "SFO" && fit_method != "DFOP") {
                # Only SFO and DFOP are current compatible with the process below
                # so skip this metabolite.
                next
            }
            
            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]
            
            metabolite0Parameter <- metabolite
            
            if (!(metabolite %in% names(odeini))) {
                metabolite0Parameter <- paste0(metabolite, "_0")
            }
            
            # If the metabolite is DFOP, we should find a k2_ decay variable,
            # otherwise we can fall back to the SFO k_ variant
            decay_var <- paste("k2", metabolite, sep="_")
            if (!(decay_var %in% names(fit$par)) &&
                !(decay_var %in% rownames(fit$fixed))) {
                decay_var <- paste("k", metabolite, sep="_")
            }
            
            # Generate the curve for DT50=1000 and ffm as fitted.
            if (decay_var %in% names(fit$par)) {
                fit$par[decay_var] <- log(2)/1000
            } else {
                # If decay_var was fixed, need to modify it there.
                fit$fixed[decay_var,"value"] <- log(2)/1000
            }
            
            new_col_name <- paste(metabolite, "DT50_1000_FFM_FITTED", sep="_")
            final[, new_col_name] <- CakeSinglePlot(fit, xlim)[metabolite]
            
            ffm1_fit <- SetFormationFractionsToOne(fit, metabolite)
            if (!is.null(ffm1_fit)) {
                # It was possible to set the formation fractions to 1, so generate the curves
                new_col_name <- paste(metabolite, "DT50_1000_FFM_1", sep="_")
                final[, new_col_name] <- CakeSinglePlot(ffm1_fit, xlim)[metabolite]
                
                # Reset the k values back to their originals
                # Then we can generate the curve for DT50 as fitted and ffm=1
                if (decay_var %in% names(ffm1_fit$par)) {
                    ffm1_fit$par[decay_var] <- original_fit$par[decay_var]
                } else {
                    ffm1_fit$fixed[decay_var,"value"] <- original_fit$fixed[decay_var,"value"]
                }
    
                new_col_name <- paste(metabolite, "DT50_FITTED_FFM_1", sep="_")
                final[, new_col_name] <- CakeSinglePlot(ffm1_fit, xlim)[metabolite]
            }
            
            # Finally reset the fit for the next metabolite
            fit <- original_fit
        }
    }
    
    cat("<PLOT MODEL START>\n")
    
    write.table(final, quote=FALSE)
    
    cat("<PLOT MODEL END>\n")

    # View(final)
}

# Get the immediate parents of a metabolite in a topology
GetParents <- function(topology, metabolite)
{
  parents <- NULL
  
  for (topo_metab in names(topology)) {
    if (metabolite %in% topology[[topo_metab]]$to) {
      parents <- append(parents, topo_metab)
    }
  }
  
  return (parents)
}

# Get all the ancestors of a metabolite in a topology
GetAncestors <- function(topology, metabolite, ancestors = NULL)
{
  for (parent in GetParents(topology, metabolite)) {
    if (!(parent %in% ancestors)) {
      ancestors <- GetAncestors(topology, parent, append(ancestors, parent))
    }
  }
  
  return (ancestors)
}

# Set the formation leading to the metabolite to 1 from each parent
# Returns the updated fit or NULL if the metabolite can't be updated for FFM=1
SetFormationFractionsToOne <- function(fit, metabolite)
{
  parents <- GetParents(fit$topology, metabolite)

  for (parent in parents) {
    if (any(GetAncestors(fit$topology, parent) %in% parents)) {
      # Reject metabolites where any parent is also an ancestor of another parent
      return (NULL)
    }
    
    # Try to find the formation fraction in the fitted and fixed parameters
    ffm_var <- paste("f",  parent, "to", metabolite, sep="_")
    updated_fit <- SetFfmVariable(fit, ffm_var, 1)
    if (!is.null(updated_fit)) {
      # Formation fraction was found and set to 1
      fit <- updated_fit
      next
    }
    
    # The formation fraction for this metabolite is implicitly 1 minus all the others from this parent
    # We need to iterate other formation fractions from the parent and set them to zero
    for (other_metabolite in fit$topology[[parent]]$to) {
      if (other_metabolite == metabolite) {
        # Skip the metabolite itself because we know there isn't a formation fraction for it
        next
      }
      
      other_ffm_var <- paste("f", parent, "to", other_metabolite, sep="_")
      updated_fit <- SetFfmVariable(fit, other_ffm_var, 0)
      
      # The updated fit shouldn't be able to be null because otherwise the topology was invalid
      # or didn't match the parameters in the fit.  If it does, we'll throw an error.
      if (is.null(updated_fit)) {
        stop(paste("The updated fit was NULL. The formation fraction '", other_ffm_var, "' doesn't exist in the fit.", sep=""))
      }
      
      fit <- updated_fit
    }
  }
  
  return (fit)
}

# Set the value of ffm_var in the fit, looking in both fitted and fixed parameters
# Returns the updated fit or NULL if the ffm_var couldn't be found
SetFfmVariable <- function(fit, ffm_var, value)
{
  if (ffm_var %in% names(fit$par)) {
    # The ffm_var was found as a fitted parameter
    fit$par[ffm_var] <- value
    return (fit)
  }
  
  if (ffm_var %in% rownames(fit$fixed)) {
    # The ffm_var was found as a fixed parameter
    fit$fixed[ffm_var,"value"] <- value
    return (fit)
  }
  
  # ffm_var was not found in the fit
  return (NULL)
}

CakeSinglePlot <- 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-5
  }

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

  # Solve the system
  odeResult <- CakeOdeSolve(fit, outtimes, solution, atol)
  
  odeResult_transformed <- PostProcessOdeOutput(odeResult, fit, atol)
  
  # Replace values that are incalculably small with 0.
  for (compartment.name in names(fit$map)) {
    if (length(odeResult_transformed[, compartment.name][!is.nan(odeResult_transformed[, compartment.name])]) > 0) {
        odeResult_transformed[, compartment.name][is.nan(odeResult_transformed[, compartment.name])] <- 0
    }
    
    odeResult_transformed[, compartment.name][odeResult_transformed[, compartment.name] < 0] <- 0
  }
  
  return (odeResult_transformed)
}