R/mread_yaml.R

In kylebaron/mread.yaml: Model Specification Using 'yaml'

default_rxn_width <- function(x,plus=1) ceiling(log10(length(x))) + plus
eq_pad <- function(x) formatC(x,width = max(nchar(x)),flag='-')

block_sep <- "\n//-------------------------------\n"
valid_names <- c("value", "label", "name", "unit", "reference")

mread_yaml <- function(file, yaml_dir = '.', project = yaml_dir, ...) {
  cpp <- yaml_to_cfile(file=file,yaml_dir=yaml_dir,project=project)
  stopifnot(file_writeable(project))
  mread_cache(cpp, project = project, ...)
}

mread_yaml_ex <- function(file,project=tempdir(),...) {
  y_dir <- system.file("examples", package="mread.yaml")
  mread_yaml(file,yaml_dir=y_dir,project=project,...)
}

yaml_to_cfile <- function(file,yaml_dir,project,ext=".model") {
  mod <- yaml.load_file(file.path(yaml_dir,file))
  ode <- make_ode(mod)
  main <- make_main(mod)
  table <- make_table(mod)
  capture <- make_capture(mod)
  global <- make_global(mod)
  param <- make_param(mod)
  init <- make_init(mod)
  prob <- make_prob(mod)
  stem <- tools::file_path_sans_ext(basename(file))
  cpp <- paste0(stem,ext)
  out <- file.path(project,cpp)
  message("Writing model to ",cpp)
  mesg <- "// Generated by yaml_to_cfile; do not edit by hand"
  src <- paste0("// Source: ", file)
  writeLines(
    c(mesg,src,block_sep,prob,global,param,init,main,ode,table,capture),
    con = out
  )
  return(cpp)
}

conf_input <- function(x,name,context) {
  if(is.numeric(x)) return(list(value = x))
  if(!all(names(x) %in% valid_names)) {
    diff <- setdiff(names(x),valid_names)
    msg <- paste0("- ", diff)
    message("context: ", context)
    message(" - name: ", name)
    message(" - invalid field name: ", diff)
    message(" - input must be either numeric or a named list")
    stop("invalid input", call.=FALSE)
  }
  if(is.null(x$value)) {
    message("context: ", context)
    message("  - the 'value' field must not be NULL")
    stop("invalid input", call.=FALSE)
  }
  x
}

make_param <- function(x) {
  if(is.null(x$param)) return(NULL)
  pnames <- names(x$param)
  p <- imap(x$param, conf_input, context = "validating param block")
  pvalues <- map_dbl(p, ~eval(parse(text = .x$value)))
  ans <- paste0(pnames, " = ", pvalues)
  c("[ param ]", ans, block_sep)
}

make_init <- function(x) {
  if(is.null(x$init)) return(NULL)
  inames <- names(x$init)
  i <- imap(x$init, conf_input, context = "validating init block")
  ivalues <- map_dbl(i,  ~eval(parse(text = .x$value)))
  ans <- paste0(inames," = ", ivalues)
  c("[ init ]", ans, block_sep)
}

make_fluxes <- function(x) {
  x <- as.data.frame(x)
  x$fluxes <- paste0("double ",x$J, " = ",x$formula, ";")
  x
}

make_dadt <- function(lhs,rhs,j_names) {
  a <- imap_dfr(lhs, expand.grid, stringsAsFactors=FALSE, reactant = TRUE)
  b <- imap_dfr(rhs, expand.grid, stringsAsFactors=FALSE, reactant = FALSE)
  dx <- bind_rows(a,b)
  names(dx) <- c("cmt", "j_number", "reactant")
  dx <- filter(dx, !(cmt %in%  c("null", "void", "NULL")))
  dx <- arrange(dx,j_number)
  dx <- mutate(dx, J = j_names[j_number])
  dx <- mutate(dx, J = if_else(reactant, paste0(" -",J), J))
  dx <- mutate(dx, sep = if_else(reactant, "", " +"))
  dx <- mutate(dx, cmtf = factor(cmt, levels = unique(cmt)))
  ode <-
    dx %>%
    group_by(cmtf) %>%
    summarise(values = paste0(sep,J, collapse = "")) %>%
    ungroup() %>%
    mutate(values = sub("^ \\+", "  ", values))
  ode <- mutate(ode, dxdt = paste0("dxdt_", eq_pad(as.character(cmtf)), " =", values, ";"))
  ode
}


reactions_to_ode <- function(eq) {
  ans <- parse_reactions(eq)
  fluxes <- make_fluxes(ans$species)
  dadt <- make_dadt(ans$lhs,ans$rhs,ans$j_names)
  list(code = c(fluxes$fluxes, " ", dadt$dxdt), ans = ans)
}

make_main <- function(x) {
  if(is.null(x$main)) {
    return(NULL)
  }
  c("[ main ]", x$main, block_sep )
}

make_table <- function(x) {
  if(is.null(x$table)) {
    return(NULL)
  }
  c("[ table ]", x$table, block_sep)
}

make_capture <- function(x) {
  if(is.null(x$capture)) {
    return(NULL)
  }
  c("[ capture ]", x$capture, block_sep)
}

make_global <- function(x) {
  if(is.null(x$global)) {
    return(NULL)
  }
  c("[ global ]", x$global, block_sep)
}

make_ode <- function(mod) {
  if(is.null(mod$reactions)) return(NULL)
  ode <- reactions_to_ode(mod$reactions)
  c("[ ode ]", mod$ode_assignments, ode$code, block_sep)
}

make_prob <- function(x) {
  if(is.null(x$prob)) return(NULL)
  c("[ prob ]", x$prob)
}