# Was Id: mkinmod.R 71 2010-09-12 01:13:36Z jranke 

# Some of the CAKE R modules are based on mkin
# Portions Johannes Ranke, 2010
# Contact: mkin-devel@lists.berlios.de

# 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/>.

CakeModel <- function(..., use_of_ff = "max")
{
  # While spec is modified throughout this function, topology is added as-is to the returned model
  topology <- list(...)
  spec <- list(...)
  obs_vars <- names(spec)

  if (!use_of_ff %in% c("min", "max")) stop("The use of formation fractions 'use_of_ff' can only be 'min' or 'max'")
  
  # differential equations, parameter names and a mapping from model variables
  # to observed variables. If possible, a matrix representation of the 
  # differential equations is included
  parms <- vector()
  diffs <- vector()
  map <- list()
  nlt <- list()

  # Establish list of differential equations
  for (varname in obs_vars)
  {
    if(is.null(spec[[varname]]$type)) stop("Every argument to CakeModel must be a list containing a type component")
    if(!spec[[varname]]$type %in% c("SFO", "FOMC", "DFOP", "HS", "IORE")) stop("Available types are SFO, FOMC, DFOP, HS and IORE only")
    
    new_parms <- vector()

    # New (sub)compartments (boxes) needed for the model type
    if (spec[[varname]]$type == "DFOP" && varname != "Parent"){
      # If this is a DFOP metabolite, we need to form a system of two differential equations
      # (see FOCUS page 137)
      new_boxes <- paste(varname, c("1", "2"), sep="_")
    } else{
      new_boxes <- varname
    }
    
    map[[varname]] <- new_boxes
    names(map[[varname]]) <- rep(spec[[varname]]$type, length(new_boxes))

    # Start a new differential equation for each new box
    new_diffs <- paste("d_", new_boxes, " =", sep="")
    
    # Turn on sink if not specified otherwise
    if(is.null(spec[[varname]]$sink)) spec[[varname]]$sink <- TRUE
    
    # Construct and add SFO term and add SFO parameters if needed
    if(spec[[varname]]$type == "SFO") {
      if (use_of_ff == "min") { # Minimum use of formation fractions
        if(spec[[varname]]$sink) {
          # From p. 53 of the FOCUS kinetics report
          k_term <- paste("k", new_boxes[[1]], "sink", sep="_")
          nonlinear_term <- paste(k_term, new_boxes[[1]], sep=" * ")
          spec[[varname]]$nlt<-nonlinear_term
          new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term)
          new_parms <- k_term
          ff <- vector()
          decline_term <- paste(nonlinear_term, "*", new_boxes[[1]])
        } else { # otherwise no decline term needed here
          decline_term <- "0" 
        }
      } else {
        k_term <- paste("k", new_boxes[[1]], sep="_")
        nonlinear_term <- paste(k_term, new_boxes[[1]], sep=" * ")
        spec[[varname]]$nlt<-nonlinear_term
        new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term)
        new_parms <- k_term
        ff <- vector()
        decline_term <- paste(nonlinear_term, "*", new_boxes[[1]])
      }
    } 

    # Construct and add FOMC term and add FOMC parameters if needed
    if(spec[[varname]]$type == "FOMC") {
      if(match(varname, obs_vars) != 1) {
        stop("Type FOMC is only possible for the first compartment, which is assumed to be the source compartment")
      }
      # From p. 53 of the FOCUS kinetics report
      nonlinear_term <- paste("(alpha/beta) * ((time/beta) + 1)^-1 *", new_boxes[[1]])
      spec[[varname]]$nlt<-nonlinear_term
      new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term)
      new_parms <- c("alpha", "beta")
      ff <- vector()
    } 

    # Construct and add DFOP term and add DFOP parameters if needed
    if(spec[[varname]]$type == "DFOP") {
      k1_term <- paste("k1", varname, sep="_")
      k2_term <- paste("k2", varname, sep="_")
      g_term <- paste("g", varname, sep="_")
      
      new_parms <- c(k1_term, k2_term, g_term)
      ff <- vector()
      
      if (varname == "Parent"){
        # From p. 57 of the FOCUS kinetics report
        # Looks like this: paste("((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) *", new_boxes[[1]])
        nonlinear_term <- paste("((", k1_term, " * ", g_term, " * exp(-", k1_term, " * time) + ", k2_term, " * (1 - ", g_term, ") * exp(-", k2_term, " * time)) / (", g_term, " * exp(-", k1_term, " * time) + (1 - ", g_term, ") * exp(-", k2_term, " * time))) *", new_boxes[[1]])
        
        spec[[varname]]$nlt<-nonlinear_term
        new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term)
      } else {
        # From p. 137 of the FOCUS kinetics report - two sub-compartments with g applied to initial conditions (in ODE pre-processing) and in formation terms.
        first_nonlinear_term <- paste(k1_term, new_boxes[[1]], sep=" * ")
        second_nonlinear_term <- paste(k2_term, new_boxes[[2]], sep=" * ")
        
        # CU-150: formation has already been accounted for on the formation into the two sub-compartments; it does not need to be factored in again on the way out
        overall_term_for_formation <- paste("(", first_nonlinear_term, " + ", second_nonlinear_term, ")")
        spec[[varname]]$nlt <- overall_term_for_formation
        
        new_diffs[[1]] <- paste(new_diffs[[1]], "-", first_nonlinear_term)
        new_diffs[[2]] <- paste(new_diffs[[2]], "-", second_nonlinear_term)
      }
    } 

    # Construct and add HS term and add HS parameters if needed
    if(spec[[varname]]$type == "HS") {
      if(match(varname, obs_vars) != 1) {
        stop("Type HS is only possible for the first compartment, which is assumed to be the source compartment")
      }
      # From p. 55 of the FOCUS kinetics report: nonlinear_term <- paste("ifelse(time <= tb, k1, k2)", "*", new_boxes[[1]]) 
      # Replaced this with the smoothed version on the next line to ease convergence.
      nonlinear_term <- paste("((k1 - k2) / ( 1 + exp( 300*(time - tb)/(tb + 0.00001) ) ) + k2)", "*", new_boxes[[1]])
      spec[[varname]]$nlt<-nonlinear_term
      new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term)
      new_parms <- c("k1", "k2", "tb")
      ff <- vector()
    } 

    # Construct and add IORE term and add IORE parameters if needed
    if(spec[[varname]]$type == "IORE") {
      k_term <- paste("k", new_boxes[[1]], sep="_")
      nonlinear_term <- paste(k_term, new_boxes[[1]], sep=" * ")
      nonlinear_term <- paste(nonlinear_term, "^N")
      spec[[varname]]$nlt<-nonlinear_term
      new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term)
      new_parms <- c(k_term, "N")
      ff <- vector()
    }

    # Add observed variable to model
    parms <- c(parms, new_parms)
    names(new_diffs) <- new_boxes
    diffs <- c(diffs, new_diffs)
  }
  
  # Transfer between compartments
  for (varname in obs_vars) {
    to <- spec[[varname]]$to
    
    if(!is.null(to)) {
      fraction_left <- NULL
      
      for (target in to) {
        target_boxes <- switch(spec[[target]]$type,
          SFO = target,
          DFOP = paste(target, c("1", "2"), sep="_"))
        
        if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS", "SFO", "IORE")) {
          if ( length(to)==1 && !spec[[varname]]$sink ) {
            # There is only one output, so no need for any flow fractions. Just add the whole flow from the parent
            formation_term <- spec[[varname]]$nlt
          } else {
            fraction_to_target <- paste("f", varname, "to", target, sep="_")
            fraction_not_to_target <- paste("(1 - ", fraction_to_target, ")", sep="")
            
            if(is.null(fraction_left)) {
              fraction_really_to_target <- fraction_to_target
              fraction_left <- fraction_not_to_target
            } else {
              # If this is the last output and there is no sink, it gets what's left
              if ( target==tail(to,1) && !spec[[varname]]$sink ) {
                 fraction_really_to_target <- fraction_left
              } else {
                 fraction_really_to_target <- fraction_to_target
                 fraction_left <- paste(fraction_left, " - ", fraction_to_target, sep="")
              }
            }
            
            ff[target] <- fraction_really_to_target
            formation_term <- paste(ff[target], "*", spec[[varname]]$nlt)
            
            # Add the flow fraction parameter (if it exists)
            if ( target!=tail(to,1) || spec[[varname]]$sink ) {
              parms <- c(parms, fraction_to_target)
            }
          }
          
          if (spec[[target]]$type == "DFOP") {
            g_term <- paste("g", target, sep="_")
            
            first_formation_term <- paste(g_term, "*", formation_term)
            second_formation_term <- paste("(1 -", g_term, ") *", formation_term)
            
            diffs[[target_boxes[[1]]]] <- paste(diffs[[target_boxes[[1]]]], "+", first_formation_term)
            diffs[[target_boxes[[2]]]] <- paste(diffs[[target_boxes[[2]]]], "+", second_formation_term)
          } else{
            diffs[[target_boxes[[1]]]] <- paste(diffs[[target_boxes[[1]]]], "+", formation_term)
          }
        }
      }
    }
  }
  
  model <- list(diffs = diffs, parms = parms, map = map, use_of_ff = use_of_ff, topology = topology)

  if (exists("ff")) {
    model$ff <- ff
  }
  
  class(model) <- "mkinmod"
  invisible(model)
}