R/modelDiagram.R

In nlmixr2plot: Nonlinear Mixed Effects Models in Population PK/PD, Plot Functions

# Automatic model diagrams (#10)
#
# The model's differential equations are parsed (not regex matched) into a
# graph of compartments (nodes) and flows (edges); the graph is then laid out
# following pharmacometric conventions and drawn with one of several engines.

#' Build a compartment graph from a model's differential equations
#'
#' The differential equations are parsed into additive terms.  A term that
#' is subtracted from one compartment and added (identically) to another is
#' mass transfer between the two compartments.  A remaining subtracted term
#' that contains the compartment's own amount is an elimination (output); a
#' remaining term that depends on other compartments is an interaction that
#' does not transport mass (for example an effect compartment or a
#' pharmacodynamic stimulation/inhibition); a remaining added term that does
#' not depend on any other compartment is a (zero-order) input.
#' Dependencies through intermediate variables (like `cp <- central/v`) are
#' followed.  A compartment that drives a transfer between two other
#' compartments (like an enzyme) is drawn as an interaction with the
#' destination.  Production or loss driven only by another compartment
#' (like `ke0*cp` in an effect compartment) is represented by the interaction
#' arrow alone, without a separate input/output arrow.
#'
#' @details
#'
#' Some limitations on how equations must be written:
#'
#' - Mass transfer is only detected when the same term (up to reordering of
#'   the factors of a product) is subtracted from the source and added to the
#'   destination, e.g. `d/dt(depot) <- -ka*depot` and
#'   `d/dt(central) <- ka*depot - ...`.  Scaled transfer (like a
#'   stoichiometric or volume conversion in only one of the equations) is
#'   shown as an elimination plus an interaction.
#'
#' - `linCmt()` models are converted to ODEs with `rxode2::linToOde()`,
#'   which requires a version of 'rxode2' that provides it.
#'
#' @param object model to diagram: a model function, an `rxode2` user
#'   interface (`rxUi`) object, a compiled `rxode2` model or a fitted
#'   `nlmixr2` object.
#' @param dosing optional character vector naming the dosing compartments.
#'   When `NULL` the dosing compartments are detected from the dosing records
#'   in `data` (a dataset without dose records has no dosing compartment);
#'   when there is no data, or it has no `evid`/`amt` columns, the first
#'   compartment (the default `rxode2` dosing compartment) is used.
#' @param data optional dataset used to detect the dosing compartments (from
#'   the dosing records' `cmt`).  For fitted models this defaults to the data
#'   the model was fit with.
#' @return a `nlmixr2ModelGraph` object; a list with:
#'
#' - `nodes`: data frame with the compartment `name`, its `role`
#'   (`"dosing"`, `"central"`, `"peripheral"`, `"transit"`, `"metabolite"`,
#'   `"effect"` or `"other"`), whether it is `dosing` and the layout
#'   coordinates `x` and `y`, and an `annotation` with the compartment's
#'   dosing properties (`lag`, `F`, `rate`, `dur`; `""` when there are none),
#'   which the diagrams show next to the compartment.
#'
#' - `edges`: data frame with `from`, `to` (`NA` for inputs/eliminations),
#'   `type` (`"transfer"`, `"elimination"`, `"input"` or `"interaction"`),
#'   `sign` (`1` when the term is added, `-1` when it is subtracted; for
#'   interactions `1` is stimulation, `-1` inhibition and `0` an effect whose
#'   direction cannot be determined from the equations, e.g. through
#'   `ifelse()` or a conditionally assigned variable),
#'   `bidirectional` (for transfers) and `label` (the model term(s)).
#' @export
#' @author Matthew L. Fidler
#' @family model diagrams
#' @examples
#' \donttest{
#' one.cmt <- function() {
#'   ini({
#'     tka <- 0.45
#'     tcl <- 1
#'     tv <- 3.45
#'     add.sd <- 0.7
#'   })
#'   model({
#'     ka <- exp(tka)
#'     cl <- exp(tcl)
#'     v <- exp(tv)
#'     d/dt(depot) <- -ka * depot
#'     d/dt(central) <- ka * depot - cl / v * central
#'     cp <- central / v
#'     cp ~ add(add.sd)
#'   })
#' }
#' modelGraph(one.cmt)
#' }
modelGraph <- function(object, dosing = NULL, data = NULL) {
  .info <- .mdModelInfo(object)
  if (is.null(data)) {
    data <- .info$data
  }
  .states <- .info$states
  if (length(.states) == 0L) {
    stop("the model has no differential equations to diagram", call. = FALSE)
  }
  if (is.null(dosing)) {
    dosing <- .mdDosingFromData(data, .states, .info$order)
    # without dosing information use rxode2's default dosing compartment; a
    # dataset without dose records has no dosing compartment
    if (is.null(dosing)) dosing <- .states[1]
  } else {
    if (!is.character(dosing)) {
      stop("'dosing' must be a character vector of compartment names", call. = FALSE)
    }
    .bad <- setdiff(dosing, .states)
    if (length(.bad) > 0L) {
      stop("'dosing' compartment(s) not in the model: ", paste(.bad, collapse = ", "), call. = FALSE)
    }
  }
  .edges <- .mdClassifyTerms(.info$terms, .states)
  .nodes <- .mdLayout(.states, .edges, dosing)
  .nodes$annotation <- .info$annotation
  structure(list(nodes = .nodes, edges = .edges), class = "nlmixr2ModelGraph")
}

#' Automatic model diagram
#'
#' Draws a compartment diagram of a model from its differential equations
#' (see [modelGraph()] for how the equations are interpreted).
#'
#' The layout follows common pharmacometric conventions: dosing and
#' absorption/transit compartments are above the compartment they feed; the
#' central compartment is in the middle with the compartments it exchanges
#' mass with (peripheral compartments) to its left; unidirectional transfer
#' (e.g. to a metabolite) and eliminations go below; compartments that
#' interact with the model without mass transfer (e.g. effect compartments or
#' pharmacodynamic models) go to the right, with their own inputs above,
#' outputs below and exchange compartments further right.
#'
#' Mass transfer is drawn with solid arrows; interactions without mass
#' transfer are dashed (with a "tee" arrow head for inhibition and a "dot"
#' arrow head when the direction is undetermined with `DiagrammeR`; dotted
#' for inhibition and dot-dashed when undetermined with `ggplot2`).
#'
#' @inheritParams modelGraph
#' @param object model to diagram (see [modelGraph()]) or a
#'   `nlmixr2ModelGraph` object.
#' @param engine drawing engine: `"DiagrammeR"` (a Graphviz htmlwidget from
#'   the 'DiagrammeR' package), `"ggplot2"` (a `ggplot` object) or `"dot"`
#'   (the Graphviz DOT source as a character string, to customize or render
#'   elsewhere).  The default is `"DiagrammeR"` when that package is
#'   installed and `"ggplot2"` otherwise; it may be changed with
#'   `options(nlmixr2plot.diagram.engine = ...)`.
#' @param labels logical; when `TRUE` label the arrows with the model terms.
#' @param ... ignored.
#' @return the diagram drawn by the requested `engine`.
#' @export
#' @author Matthew L. Fidler
#' @family model diagrams
#' @examples
#' \donttest{
#' pk.turnover.emax <- function() {
#'   ini({
#'     tktr <- log(1)
#'     tka <- log(1)
#'     tcl <- log(0.1)
#'     tv <- log(10)
#'     poplogit <- 2
#'     tec50 <- log(0.5)
#'     tkout <- log(0.05)
#'     te0 <- log(100)
#'     prop.err <- 0.1
#'     pkadd.err <- 0.1
#'     pdadd.err <- 10
#'   })
#'   model({
#'     ktr <- exp(tktr)
#'     ka <- exp(tka)
#'     cl <- exp(tcl)
#'     v <- exp(tv)
#'     emax <- expit(poplogit)
#'     ec50 <- exp(tec50)
#'     kout <- exp(tkout)
#'     e0 <- exp(te0)
#'     DCP <- center / v
#'     PD <- 1 - emax * DCP / (ec50 + DCP)
#'     effect(0) <- e0
#'     kin <- e0 * kout
#'     d/dt(depot) <- -ktr * depot
#'     d/dt(gut) <- ktr * depot - ka * gut
#'     d/dt(center) <- ka * gut - cl / v * center
#'     d/dt(effect) <- kin * PD - kout * effect
#'     cp <- center / v
#'     cp ~ prop(prop.err) + add(pkadd.err)
#'     effect ~ add(pdadd.err)
#'   })
#' }
#' modelDiagram(pk.turnover.emax, engine = "ggplot2")
#' if (requireNamespace("DiagrammeR", quietly = TRUE)) {
#'   modelDiagram(pk.turnover.emax, engine = "DiagrammeR")
#' }
#' }
modelDiagram <- function(
  object,
  dosing = NULL,
  data = NULL,
  engine = getOption("nlmixr2plot.diagram.engine"),
  labels = FALSE,
  ...
) {
  if (inherits(object, "nlmixr2ModelGraph")) {
    .graph <- object
  } else {
    .graph <- modelGraph(object, dosing = dosing, data = data)
  }
  if (is.null(engine)) {
    engine <- if (requireNamespace("DiagrammeR", quietly = TRUE)) {
      "DiagrammeR"
    } else {
      "ggplot2"
    }
  }
  engine <- match.arg(engine, c("DiagrammeR", "ggplot2", "dot"))
  if (!(is.logical(labels) && length(labels) == 1L && !is.na(labels))) {
    stop("'labels' must be TRUE or FALSE", call. = FALSE)
  }
  switch(
    engine,
    DiagrammeR = .mdDiagrammeR(.graph, labels),
    ggplot2 = .mdGgplot(.graph, labels),
    dot = .mdDot(.graph, labels)
  )
}

#' @rdname modelDiagram
#' @param x a `nlmixr2ModelGraph` object, an `rxode2` user interface
#'   (`rxUi`) object or a compiled `rxode2` model
#' @details `plot()` of an `rxode2` user interface (`rxUi`) object, like
#'   `rxode2::rxode2(modelFunction)`, or of a compiled `rxode2` model draws
#'   its model diagram, so `plot(rxode2(model))` is the same as
#'   `modelDiagram(model)`.  (A fitted `nlmixr2` model keeps its
#'   goodness-of-fit `plot()`; use `modelDiagram(fit)` for its diagram.)
#' @export
plot.nlmixr2ModelGraph <- function(x, ..., engine = getOption("nlmixr2plot.diagram.engine"), labels = FALSE) {
  modelDiagram(x, engine = engine, labels = labels)
}

#' @rdname modelDiagram
#' @export
plot.rxUi <- function(
  x,
  ...,
  dosing = NULL,
  data = NULL,
  engine = getOption("nlmixr2plot.diagram.engine"),
  labels = FALSE
) {
  modelDiagram(x, dosing = dosing, data = data, engine = engine, labels = labels)
}

#' @rdname modelDiagram
#' @export
plot.rxode2 <- plot.rxUi

#' @export
print.nlmixr2ModelGraph <- function(x, ...) {
  cat("nlmixr2 model graph\n\ncompartments:\n")
  .n <- x$nodes[, c("name", "role", "dosing")]
  .a <- x$nodes$annotation
  if (!is.null(.a) && any(!is.na(.a) & nzchar(.a))) {
    .a[is.na(.a)] <- ""
    .n$annotation <- gsub("\n", "; ", .a, fixed = TRUE)
  }
  print(.n, row.names = FALSE)
  cat("\nflows:\n")
  .e <- x$edges
  if (nrow(.e) == 0L) {
    cat("  (none)\n")
  } else {
    .e$from[is.na(.e$from)] <- "(input)"
    .e$to[is.na(.e$to)] <- "(output)"
    print(.e[, c("from", "to", "type", "sign", "label")], row.names = FALSE)
  }
  invisible(x)
}

# ---------------------------------------------------------------------------
# Model extraction
# ---------------------------------------------------------------------------

#' Get the model lines, states and data from a supported object
#'
#' @param object model object
#' @return list(states, terms, data)
#' @noRd
.mdModelInfo <- function(object) {
  .data <- NULL
  if (inherits(object, "nlmixr2FitCore") || inherits(object, "nlmixr2FitData")) {
    .data <- tryCatch(object$origData, error = function(e) NULL)
    object <- object$ui
  }
  if (inherits(object, "rxode2")) {
    .mv <- rxode2::rxModelVars(object)
    # braces let a normalized `}\nelse` parse
    .lines <- as.list(str2lang(paste0("{\n", .mv$model["normModel"], "\n}")))[-1]
    .order <- .mv$state
  } else {
    if (is.function(object)) {
      object <- rxode2::rxode2(object)
    }
    if (!inherits(object, "rxUi")) {
      stop(
        "cannot create a model diagram from an object of class '",
        paste(class(object), collapse = "', '"),
        "'",
        call. = FALSE
      )
    }
    if (.mdIsLinCmt(object)) {
      object <- .mdLinToOde(object)
    }
    .lines <- object$lstExpr
    # residual error lines (`cp ~ add(sd)`) are not assignments
    .err <- object$predDf$line
    if (length(.err) > 0L) {
      .lines <- .lines[-.err]
    }
    .order <- object$mv0$state
  }
  # substitute definitions into the equations; when that makes an equation
  # too large (QSP/PBPK models), redo the whole model with smaller
  # substitutions so every equation is written the same way
  .parsed <- .mdParseLines(.lines)
  for (.level in c("small", "none")) {
    if (!isTRUE(.parsed$overflow)) {
      break
    }
    .parsed <- .mdParseLines(.lines, .level)
  }
  .states <- .parsed$states
  # keep rxode2's compartment order (used to map numeric `cmt` values)
  .states <- c(intersect(.order, .states), setdiff(.states, .order))
  .foldIndex <- .mdFoldIndex(.parsed$defs)
  .terms <- do.call(
    rbind,
    lapply(.states, function(.s) {
      .t <- .parsed$ode[[.s]]
      if (length(.t) == 0L) {
        return(NULL)
      }
      data.frame(
        state = .s,
        sign = vapply(.t, function(x) x$sign, numeric(1)),
        key = vapply(.t, function(x) .mdTermKey(x$expr), character(1)),
        label = vapply(
          .t,
          function(x) {
            # `flux#2` (one name per assignment) is shown as `flux`
            gsub("#[0-9]+", "", .mdDeparse(.mdFold(x$expr, .parsed$defs, .foldIndex)), perl = TRUE)
          },
          character(1)
        ),
        stringsAsFactors = FALSE
      )
    })
  )
  if (!is.null(.terms)) {
    .terms$states <- unlist(
      lapply(.states, function(.s) {
        lapply(.parsed$ode[[.s]], function(x) {
          .mdExprStates(x$expr, .states, x$snap)
        })
      }),
      recursive = FALSE
    )
    # direction (1 increasing, -1 decreasing, NA unknown) of each term in
    # each compartment it depends on
    .terms$dir <- unlist(
      lapply(.states, function(.s) {
        lapply(.parsed$ode[[.s]], function(x) {
          .st <- .mdExprStates(x$expr, .states, x$snap)
          .dep <- lapply(.st, .mdDependents, states = .states, closure = x$snap)
          stats::setNames(.mdDirections(x$expr, .st, .dep), .st)
        })
      }),
      recursive = FALSE
    )
  }
  # rxode2's full compartment order, including compartments without ODEs
  # (e.g. from `cmt()`), maps numeric `cmt` values
  list(
    states = .states,
    terms = .terms,
    data = .data,
    annotation = .mdAnnotation(.parsed$props, .states),
    order = c(.order, setdiff(.states, .order))
  )
}

#' Dosing compartments from the dosing records of a dataset
#'
#' @param data dataset (or NULL)
#' @param states compartment names that are diagrammed
#' @param order rxode2's compartment order (numeric `cmt` values index it)
#' @return character vector of dosed compartments (empty when the data has
#'   no dose records) or `NULL` when the data has no dosing information
#' @noRd
.mdDosingFromData <- function(data, states, order = states) {
  if (!is.data.frame(data) || nrow(data) == 0L) {
    return(NULL)
  }
  .nm <- tolower(names(data))
  .col <- function(n) {
    .w <- which(.nm == n)
    if (length(.w) == 0L) {
      return(NULL)
    }
    data[[.w[1]]]
  }
  .evid <- .col("evid")
  .amt <- .col("amt")
  if (is.null(.evid) && is.null(.amt)) {
    return(NULL)
  }
  .dose <- rep(TRUE, nrow(data))
  # factors are converted through their labels, not their level codes
  .num <- function(v) suppressWarnings(as.numeric(as.character(v)))
  if (!is.null(.evid)) {
    .evid <- .num(.evid)
    # 0 = observation, 2 = other event, 3 = reset
    .dose <- .dose & !is.na(.evid) & !(.evid %in% c(0, 2, 3))
  }
  if (!is.null(.amt)) {
    .amt <- .num(.amt)
    .dose <- .dose & !is.na(.amt) & .amt != 0
  }
  if (!any(.dose)) {
    return(character(0))
  }
  .cmt <- .col("cmt")
  if (is.null(.cmt)) {
    return(intersect(states, order[1]))
  }
  .cmt <- .cmt[.dose]
  if (is.factor(.cmt)) {
    .cmt <- as.character(.cmt)
  }
  # a missing compartment doses the default (first) compartment
  .ret <- if (anyNA(.cmt)) order[1] else character(0)
  .cmt <- .cmt[!is.na(.cmt)]
  if (is.character(.cmt)) {
    .num <- suppressWarnings(as.numeric(.cmt))
    .chr <- .cmt[is.na(.num)]
    .chr[.chr %in% c("(default)", "")] <- order[1]
    .ret <- c(.ret, .chr)
    .cmt <- .num[!is.na(.num)]
  } else {
    .cmt <- as.numeric(.cmt)
  }
  # negative compartment numbers turn compartments off; they are not doses
  .cmt <- .cmt[!is.na(.cmt) & .cmt > 0 & .cmt <= length(order)]
  .ret <- c(.ret, order[.cmt])
  intersect(states, .ret)
}

#' @noRd
.mdIsLinCmt <- function(ui) {
  any(vapply(
    ui$lstExpr,
    function(x) {
      "linCmt" %in% all.names(x)
    },
    logical(1)
  ))
}

#' @noRd
.mdLinToOde <- function(ui) {
  if (!("linToOde" %in% getNamespaceExports("rxode2"))) {
    stop("diagramming 'linCmt()' models requires a version of 'rxode2' with 'linToOde()'", call. = FALSE)
  }
  .fun <- getExportedValue("rxode2", "linToOde")
  .fun(ui)
}

#' Parse model lines into ODE terms and variable -> state dependencies
#'
#' Variables that are assigned once (outside of `if` blocks) and depend on a
#' compartment amount are substituted into the differential equations before
#' they are split into terms, so that the direction of an effect expressed
#' through an intermediate variable (e.g. `PD <- 1 - emax*cp/(ec50 + cp)`) is
#' seen.
#'
#' @param lines list of model expressions
#' @return list(ode = named list of term lists, deps = named list of
#'   variables to the (possibly empty) set of variables they depend on)
#' @noRd
.mdParseLines <- function(lines, level = c("full", "small", "none")) {
  level <- match.arg(level)
  .env <- new.env(parent = emptyenv())
  .env$ode <- list()
  # compartments each variable stands for, as the model is read (a later
  # assignment must not change what an earlier use meant)
  .env$stateOf <- list()
  .env$defs <- list()
  .env$count <- list()
  .env$inIf <- character(0)
  .env$states <- character(0)
  .env$props <- list()
  .env$defsSmall <- list()
  .env$overflow <- FALSE
  .env$version <- list()
  .env$versioned <- list()
  # the definition to substitute for variable `n` (NULL: keep the variable)
  .mdDefinition <- function(n, rhs, maxSize) {
    if (n %in% .env$states || n %in% .env$inIf || length(all.names(rhs)) > maxSize) {
      # values from `if` branches cannot be substituted, and very large
      # definitions (QSP/PBPK models) are kept as variables so the
      # equations do not grow combinatorially.  A reassigned variable that
      # is not substituted keeps one name per assignment, so that two
      # different values (`flux = k*A` ... `flux = h*B`) are not mistaken
      # for one flow
      return(.env$versioned[[n]])
    }
    # a constant (`x = -1`) is substituted so that its value and sign are
    # known where it is used
    if (length(all.vars(rhs)) == 0L) {
      return(rhs)
    }
    if (!identical(.env$count[[n]], 1)) {
      # a reassigned variable is substituted with its current value so
      # that a reused name (like `flux`) is not mistaken for one flow
      return(rhs)
    }
    if (any(all.vars(rhs) %in% .env$states)) {
      return(rhs)
    }
    NULL
  }
  .isAssign <- function(x) {
    length(x) == 3L &&
      (identical(x[[1]], quote(`<-`)) || identical(x[[1]], quote(`=`)) || identical(x[[1]], quote(`~`)))
  }
  # first pass: compartment names and how often (and where) variables are
  # assigned
  .count <- function(x, inIf) {
    if (!is.call(x)) {
      return(invisible())
    }
    .f <- x[[1]]
    if (identical(.f, quote(`{`))) {
      for (.i in seq_along(x)[-1]) {
        .count(x[[.i]], inIf)
      }
    } else if (identical(.f, quote(`if`))) {
      .count(x[[3]], TRUE)
      if (length(x) == 4L) .count(x[[4]], TRUE)
    } else if (.isAssign(x)) {
      .state <- .mdDdtState(x[[2]])
      if (!is.null(.state)) {
        .env$states <- union(.env$states, .state)
      } else if (is.name(x[[2]])) {
        .n <- as.character(x[[2]])
        .env$count[[.n]] <- (if (is.null(.env$count[[.n]])) 0 else .env$count[[.n]]) + 1
        if (inIf) .env$inIf <- union(.env$inIf, .n)
      }
    }
    invisible()
  }
  for (.l in lines) {
    .count(.l, FALSE)
  }
  # add terms to a compartment's equation; a repeated term (`-k*A - k*A`)
  # is kept, since each copy moves mass (if/else branches are merged
  # before they get here)
  .addTerms <- function(state, terms) {
    .keep <- Filter(
      function(t) {
        # a literal zero (e.g. `d/dt(x) <- 0`) is no flow
        !(is.numeric(t$expr) && all(t$expr == 0))
      },
      terms
    )
    # what the variables mean here, for the compartments of this term
    .snap <- .env$stateOf
    .keep <- lapply(.keep, function(t) {
      if (is.null(t$snap)) {
        t$snap <- .snap
      }
      t
    })
    .env$ode[[state]] <- c(.env$ode[[state]], .keep)
  }
  # second pass; `cond` holds the variables of the enclosing `if` conditions
  .walk <- function(x, cond) {
    if (!is.call(x)) {
      return(invisible())
    }
    .f <- x[[1]]
    if (identical(.f, quote(`{`))) {
      for (.i in seq_along(x)[-1]) {
        .walk(x[[.i]], cond)
      }
    } else if (identical(.f, quote(`if`))) {
      .cond <- union(cond, all.vars(x[[2]]))
      # walk each branch separately, then merge: terms in both branches
      # apply either way, the others become `ifelse(cond, term, 0)` /
      # `ifelse(cond, 0, term)` so that the condition is kept
      .base <- .env$ode
      .baseState <- .env$stateOf
      .env$ode <- list()
      .walk(x[[3]], .cond)
      .yes <- .env$ode
      .yesState <- .env$stateOf
      .env$ode <- list()
      .env$stateOf <- .baseState
      if (length(x) == 4L) {
        .walk(x[[4]], .cond)
      }
      .no <- .env$ode
      .noState <- .env$stateOf
      # a variable assigned in a branch may hold either branch's value
      .merged <- .noState
      for (.v in union(names(.yesState), names(.noState))) {
        .merged[[.v]] <- union(.yesState[[.v]], .noState[[.v]])
      }
      .env$stateOf <- .merged
      .env$ode <- .base
      for (.s in union(names(.yes), names(.no))) {
        .addTerms(.s, .mdMergeBranches(x[[2]], .yes[[.s]], .no[[.s]]))
      }
    } else if (.isAssign(x)) {
      .lhs <- x[[2]]
      # full substitution, and substitution of small definitions only (used
      # when full substitution makes an equation too large)
      .rhs <- switch(
        level,
        full = .mdSubstitute(x[[3]], .env$defs),
        small = .mdSubstitute(x[[3]], .env$defsSmall),
        none = .mdSubstitute(x[[3]], .env$versioned)
      )
      .rhsSmall <- .mdSubstitute(x[[3]], .env$defsSmall)
      .state <- .mdDdtState(.lhs)
      .prop <- .mdDoseProperty(.lhs)
      if (!is.null(.prop)) {
        # dosing properties (lag, F, rate, dur) are annotations only
        .old <- .env$props[[.prop$state]]
        .val <- .mdDeparse(x[[3]])
        .cur <- .old[.prop$name]
        .old[.prop$name] <- if (is.null(.old) || is.na(.cur)) {
          .val
        } else {
          paste(.cur, .val, sep = " / ")
        }
        .env$props[[.prop$state]] <- .old
      } else if (!is.null(.state)) {
        # an equation that is too large after substitution makes the whole
        # model fall back to a lower substitution level, so that the same
        # flow is written the same way in every equation
        if (length(all.names(.rhs)) > .mdMaxOdeSize) {
          .env$overflow <- TRUE
        }
        .addTerms(.state, .mdSplitTerms(.rhs, .env$states, .state))
      } else if (is.name(.lhs)) {
        .n <- as.character(.lhs)
        # the compartments this value stands for, resolved now (a
        # conditional assignment also depends on its condition)
        .dep <- union(all.vars(x[[3]]), cond)
        .set <- unique(unlist(c(
          intersect(.dep, .env$states),
          lapply(setdiff(.dep, .env$states), function(.v) .env$stateOf[[.v]])
        )))
        if (is.null(.set)) {
          .set <- character(0)
        }
        .env$stateOf[[.n]] <- .set
        if (!identical(.env$count[[.n]], 1) && !(.n %in% .env$inIf)) {
          # one name per assignment of a reassigned variable; a variable
          # assigned in `if` branches keeps one name, since its value is one
          # of the branches (and depends on all of them)
          .v <- (if (is.null(.env$version[[.n]])) 0L else .env$version[[.n]]) + 1L
          .env$version[[.n]] <- .v
          .vn <- paste0(.n, "#", .v)
          .env$versioned[[.n]] <- as.name(.vn)
          .env$stateOf[[.vn]] <- .set
        }
        # a NULL definition removes the variable from the substitutions
        .env$defs[[.n]] <- .mdDefinition(.n, .rhs, .mdMaxDefSize)
        .env$defsSmall[[.n]] <- .mdDefinition(.n, .rhsSmall, .mdMaxSmallDefSize)
      }
    }
    invisible()
  }
  for (.l in lines) {
    .walk(.l, character(0))
  }
  # only single-assignment definitions can be folded back into labels
  .defs <- c(.env$defs, .env$defsSmall)
  .defs <- .defs[!vapply(.defs, is.null, logical(1))]
  .fold <- .defs[vapply(
    names(.defs),
    function(.n) {
      identical(.env$count[[.n]], 1)
    },
    logical(1)
  )]
  list(ode = .env$ode, defs = .fold, states = .env$states, props = .env$props, overflow = .env$overflow)
}

#' Replace substituted definitions by their variable names (for labels)
#' @noRd
.mdFold <- function(expr, defs, index = .mdFoldIndex(defs)) {
  if (length(defs) == 0L) {
    return(expr)
  }
  .fold <- function(e) {
    if (is.call(e)) {
      # only deparse sub-expressions the size of some definition
      if (length(all.names(e)) %in% index$size) {
        .w <- which(index$key == .mdDeparse(e))
        if (length(.w) > 0L) return(as.name(names(defs)[.w[length(.w)]]))
      }
      for (.i in seq_along(e)[-1]) {
        .v <- .fold(e[[.i]])
        if (!is.null(.v)) e[[.i]] <- .v
      }
    }
    e
  }
  .fold(expr)
}

#' Deparsed keys and sizes of definitions, computed once per model
#' @noRd
.mdFoldIndex <- function(defs) {
  list(key = vapply(defs, .mdDeparse, character(1)), size = vapply(defs, function(d) length(all.names(d)), numeric(1)))
}

#' Substitute variable definitions into an expression
#' @noRd
.mdSubstitute <- function(expr, defs) {
  if (length(defs) == 0L) {
    return(expr)
  }
  do.call(substitute, list(expr, defs))
}

#' Dosing property set by an assignment like `alag(depot) <- tlag`
#'
#' @param lhs left hand side of an assignment
#' @return list(state, name) with the display name (`lag`, `F`, `rate` or
#'   `dur`), or NULL
#' @noRd
.mdDoseProperty <- function(lhs) {
  if (!is.call(lhs) || length(lhs) != 2L || !is.name(lhs[[1]]) || !is.name(lhs[[2]])) {
    return(NULL)
  }
  .n <- c(lag = "lag", alag = "lag", f = "F", F = "F", rate = "rate", dur = "dur")[as.character(lhs[[1]])]
  if (is.na(.n)) {
    return(NULL)
  }
  list(state = as.character(lhs[[2]]), name = unname(.n))
}

#' Annotation text for each compartment's dosing properties
#' @param props named list (by compartment) of named character vectors
#' @param states compartment names
#' @return character vector (`""` without properties)
#' @noRd
.mdAnnotation <- function(props, states) {
  vapply(
    states,
    function(.s) {
      .p <- props[[.s]]
      if (is.null(.p) || length(.p) == 0L) {
        return("")
      }
      .p <- .p[intersect(c("lag", "F", "rate", "dur"), names(.p))]
      paste0(names(.p), " = ", .p, collapse = "\n")
    },
    character(1),
    USE.NAMES = FALSE
  )
}

#' Return the state name for a `d/dt(state)` expression or NULL
#' @noRd
.mdDdtState <- function(lhs) {
  if (
    is.call(lhs) &&
      identical(lhs[[1]], quote(`/`)) &&
      identical(lhs[[2]], quote(d)) &&
      is.call(lhs[[3]]) &&
      identical(lhs[[3]][[1]], quote(dt)) &&
      length(lhs[[3]]) == 2L
  ) {
    return(as.character(lhs[[3]][[2]]))
  }
  NULL
}

#' States an expression depends on (directly or through variables)
#' @param closure named list of the compartments each variable stands for
#' @noRd
.mdExprStates <- function(expr, states, closure) {
  .v <- all.vars(expr)
  # a state's own amount is not expanded through assignments
  .direct <- .v[.v %in% states]
  .via <- unlist(closure[setdiff(.v, states)], use.names = FALSE)
  intersect(states, c(.direct, .via))
}

# limits that keep large (QSP/PBPK) models fast: the size (number of names)
# of a definition substituted into the equations, and the number of terms a
# single product may be distributed into
.mdMaxDefSize <- 2000L
.mdMaxTerms <- 64L
# beyond this size (after substitution) an equation only gets definitions up
# to `.mdMaxSmallDefSize` substituted
.mdMaxOdeSize <- 5000L
# a summand is only expanded when its terms stay within this multiple of its
# size
.mdMaxTermGrowth <- 50L
# more neighbors than this in one direction are fanned out on an arc
.mdMaxLine <- 4L
# models with more compartments than this use a bounded placement search
.mdMaxWideSearch <- 100L
# graphs with more edges than this are drawn with straight edges in DOT
.mdMaxSplineEdges <- 200L
.mdMaxSmallDefSize <- 200L

#' Split an equation into terms, bounding the growth of each summand
#'
#' The top-level sum is split first (which never grows); each summand is then
#' expanded with `.mdTerms()` unless that would make it much larger (e.g. a
#' long numerator split over a large denominator), in which case it is kept
#' as one term.
#' @param x expression
#' @param states compartment names
#' @return list of list(sign, expr)
#' @noRd
.mdSplitTerms <- function(x, states = character(0), own = NULL) {
  .top <- .mdTopSum(x)
  # expand within a size budget; NULL when the budget runs out
  .try <- function(expr, own, limit) {
    .budget <- new.env(parent = emptyenv())
    .budget$left <- limit
    tryCatch(.mdTerms(expr, states, own, .budget), mdBudget = function(e) NULL)
  }
  unlist(
    lapply(.top, function(.t) {
      .limit <- .mdMaxTermGrowth * length(all.names(.t$expr)) + 200
      .e <- .try(.t$expr, NULL, .limit)
      if (is.null(.e) && !is.null(own)) {
        # coarse expansion (own compartment / other compartments / constants)
        .e <- .try(.t$expr, own, .limit)
      }
      if (is.null(.e)) {
        return(list(.t))
      }
      if (.t$sign < 0) {
        .e <- .mdNeg(.e)
      }
      .e
    }),
    recursive = FALSE
  )
}

#' Total size (number of names) of a list of terms
#' @noRd
.mdSize <- function(terms) {
  sum(vapply(terms, function(t) length(all.names(t$expr)), numeric(1)))
}

#' Charge an expansion budget, signalling `mdBudget` when it runs out
#' @noRd
.mdCharge <- function(budget, size) {
  if (is.null(budget)) {
    return(invisible())
  }
  budget$left <- budget$left - size
  if (budget$left < 0) {
    stop(structure(
      class = c("mdBudget", "error", "condition"),
      list(message = "expansion budget exceeded", call = NULL)
    ))
  }
  invisible()
}

#' Split only the top-level sum of an expression into signed summands
#' @noRd
.mdTopSum <- function(x, sign = 1) {
  if (is.call(x)) {
    .f <- x[[1]]
    if (identical(.f, quote(`(`))) {
      return(.mdTopSum(x[[2]], sign))
    }
    if (identical(.f, quote(`+`))) {
      if (length(x) == 2L) {
        return(.mdTopSum(x[[2]], sign))
      }
      return(c(.mdTopSum(x[[2]], sign), .mdTopSum(x[[3]], sign)))
    }
    if (identical(.f, quote(`-`))) {
      if (length(x) == 2L) {
        return(.mdTopSum(x[[2]], -sign))
      }
      return(c(.mdTopSum(x[[2]], sign), .mdTopSum(x[[3]], -sign)))
    }
  }
  list(list(sign = sign, expr = x))
}

#' Split an expression into signed additive terms
#'
#' Sums and differences are split, products distribute over sums and a sum
#' in a numerator is split over its denominator.
#' @param x expression
#' @param states compartment names, used to group long sums (see
#'   `.mdCollapse()`)
#' @param own when given, the equation's own compartment: products are only
#'   expanded into terms with and without it, with other compartments and
#'   without any (a coarse expansion that stays small)
#' @return list of list(sign, expr)
#' @noRd
.mdTerms <- function(x, states = character(0), own = NULL, budget = NULL) {
  if (is.call(x)) {
    .f <- x[[1]]
    if (identical(.f, quote(`(`))) {
      return(.mdTerms(x[[2]], states, own, budget))
    }
    if (identical(.f, quote(`+`))) {
      if (length(x) == 2L) {
        return(.mdTerms(x[[2]], states, own, budget))
      }
      return(c(.mdTerms(x[[2]], states, own, budget), .mdTerms(x[[3]], states, own, budget)))
    }
    if (identical(.f, quote(`-`))) {
      if (length(x) == 2L) {
        return(.mdNeg(.mdTerms(x[[2]], states, own, budget)))
      }
      return(c(.mdTerms(x[[2]], states, own, budget), .mdNeg(.mdTerms(x[[3]], states, own, budget))))
    }
    if (identical(.f, quote(`*`)) && length(x) == 3L) {
      .a <- .mdTerms(x[[2]], states, own, budget)
      .b <- .mdTerms(x[[3]], states, own, budget)
      # distributing a product of long sums grows combinatorially; beyond a
      # limit first group like terms (same sign and variables) of the longer
      # factor, then of the other, and only then keep the product whole
      if (!is.null(own)) {
        # coarse expansion: only keep apart what the diagram needs (own
        # compartment / other compartments / constants)
        .a <- .mdCollapse(.a, states, own = own)
        .b <- .mdCollapse(.b, states, own = own)
      }
      for (.coarse in c(FALSE, TRUE)) {
        if (length(.a) * length(.b) > .mdMaxTerms) {
          .a <- .mdCollapse(.a, states, .coarse)
          .b <- .mdCollapse(.b, states, .coarse)
        }
      }
      if (length(.a) * length(.b) > .mdMaxTerms) {
        return(list(list(sign = 1, expr = x)))
      }
      .mdCharge(budget, length(.b) * .mdSize(.a) + length(.a) * .mdSize(.b))
      .ret <- list()
      for (.i in .a) {
        for (.j in .b) {
          .ret[[length(.ret) + 1L]] <-
            list(sign = .i$sign * .j$sign, expr = .mdMult(.i$expr, .j$expr))
        }
      }
      return(.ret)
    }
    if (identical(.f, quote(`/`)) && length(x) == 3L) {
      .den <- .mdTerms(x[[3]], states, own, budget)
      .denSign <- 1
      if (length(.den) == 1L) {
        # keep a single signed denominator's sign on the term
        .denSign <- .den[[1]]$sign
        .denExpr <- .den[[1]]$expr
      } else {
        .denExpr <- x[[3]]
      }
      .num <- .mdTerms(x[[2]], states, own, budget)
      # each numerator term carries a copy of the denominator; when that is
      # too much, group like numerator terms first so that a production and
      # a loss over the same denominator stay apart
      .denSize <- length(all.names(.denExpr))
      if (!is.null(budget) && length(.num) * .denSize > budget$left) {
        .num <- .mdCollapse(.num, states, own = own)
      }
      if (!is.null(budget) && length(.num) * .denSize > budget$left) {
        .num <- .mdCollapse(.num, states, coarse = TRUE)
      }
      .mdCharge(budget, length(.num) * .denSize)
      return(lapply(.num, function(.t) {
        list(sign = .t$sign * .denSign, expr = as.call(list(quote(`/`), .t$expr, .denExpr)))
      }))
    }
  }
  if (is.call(x) && identical(x[[1]], quote(ifelse)) && length(x) == 4L) {
    # distribute `ifelse(cond, a, b)` over the terms of each branch:
    # `ifelse(cond, a1, 0) + ... + ifelse(cond, 0, b1) + ...`; a term in both
    # branches applies either way
    .isZero <- function(t) is.numeric(t$expr) && all(t$expr == 0)
    return(.mdMergeBranches(
      x[[2]],
      Filter(Negate(.isZero), .mdTerms(x[[3]], states, own, budget)),
      Filter(Negate(.isZero), .mdTerms(x[[4]], states, own, budget))
    ))
  }
  if (is.numeric(x) && length(x) == 1L && !is.na(x) && x < 0) {
    return(list(list(sign = -1, expr = -x)))
  }
  list(list(sign = 1, expr = x))
}

#' Merge the terms of the two branches of a condition
#'
#' A term in both branches applies either way; the others become
#' `ifelse(cond, term, 0)` or `ifelse(cond, 0, term)`.  Repeated terms are
#' paired one-to-one, so `-k*A` against `-k*A - k*A` leaves one conditional
#' `-k*A`.
#' @param cond condition expression
#' @param yes,no lists of terms of each branch
#' @return list of terms
#' @noRd
.mdMergeBranches <- function(cond, yes, no) {
  .id <- function(t) paste(t$sign, .mdTermKey(t$expr))
  .wrap <- function(t, isYes) {
    list(sign = t$sign, expr = as.call(list(quote(ifelse), cond, if (isYes) t$expr else 0, if (isYes) 0 else t$expr)))
  }
  .noId <- vapply(no, .id, character(1))
  .noUsed <- rep(FALSE, length(no))
  .ret <- list()
  for (.t in yes) {
    .w <- which(!.noUsed & .noId == .id(.t))
    if (length(.w) > 0L) {
      .noUsed[.w[1]] <- TRUE
      .ret[[length(.ret) + 1L]] <- .t
    } else {
      .ret[[length(.ret) + 1L]] <- .wrap(.t, TRUE)
    }
  }
  c(.ret, lapply(no[!.noUsed], .wrap, isYes = FALSE))
}

#' Canonical form of an expression: operands of sums and products sorted
#' @noRd
.mdCanon <- function(x) {
  if (!is.call(x)) {
    return(x)
  }
  if (identical(x[[1]], quote(`(`))) {
    return(.mdCanon(x[[2]]))
  }
  # `delay(x, tau)` moves the same mass as `x`, only later: a delayed
  # transfer `-ka*depot` / `+ka*delay(depot, tlag)` is still one flow
  if (identical(x[[1]], quote(delay)) && length(x) >= 2L) {
    return(.mdCanon(x[[2]]))
  }
  for (.op in list(quote(`+`), quote(`*`))) {
    if (identical(x[[1]], .op) && length(x) == 3L) {
      .ops <- list()
      .flat <- function(e) {
        if (is.call(e) && identical(e[[1]], quote(`(`))) {
          return(.flat(e[[2]]))
        }
        if (is.call(e) && identical(e[[1]], .op) && length(e) == 3L) {
          .flat(e[[2]])
          .flat(e[[3]])
        } else {
          .ops[[length(.ops) + 1L]] <<- .mdCanon(e)
        }
      }
      .flat(x)
      .ops <- .ops[order(vapply(.ops, .mdDeparse, character(1)))]
      return(Reduce(function(a, b) as.call(list(.op, a, b)), .ops))
    }
  }
  for (.i in seq_along(x)[-1]) {
    .v <- .mdCanon(x[[.i]])
    if (!is.null(.v)) x[[.i]] <- .v
  }
  x
}

#' Group like terms of a long sum into one term
#'
#' Terms are grouped by sign and the compartments they contain (or, when
#' `coarse`, only whether they contain any compartment), which keeps the
#' constant, self-dependent and compartment-driven parts of a long sum apart
#' while bounding the number of terms.
#' @param terms list of terms
#' @param states compartment names (without them, the variables are used)
#' @param coarse group only by sign and whether a compartment appears
#' @param own group by sign, whether `own` appears and whether other
#'   compartments appear
#' @noRd
.mdCollapse <- function(terms, states = character(0), coarse = FALSE, own = NULL) {
  .key <- vapply(
    terms,
    function(t) {
      .v <- all.vars(t$expr)
      if (length(states) > 0L) {
        .v <- intersect(.v, states)
      }
      .v <- if (!is.null(own)) {
        paste(own %in% .v, length(setdiff(.v, own)) > 0L)
      } else if (coarse) {
        as.character(length(.v) > 0L)
      } else {
        paste(sort(.v), collapse = ",")
      }
      paste(t$sign, .v)
    },
    character(1)
  )
  lapply(split(terms, factor(.key, levels = unique(.key))), function(g) {
    list(
      sign = g[[1]]$sign,
      expr = Reduce(function(a, b) as.call(list(quote(`+`), a, b)), lapply(g, function(t) t$expr))
    )
  })
}

#' @noRd
.mdNeg <- function(terms) {
  lapply(terms, function(.t) {
    .t$sign <- -.t$sign
    .t
  })
}

#' Multiply two expressions, dropping multiplications by one
#' @noRd
.mdMult <- function(a, b) {
  if (is.numeric(a) && length(a) == 1L && a == 1) {
    return(b)
  }
  if (is.numeric(b) && length(b) == 1L && b == 1) {
    return(a)
  }
  as.call(list(quote(`*`), a, b))
}

#' Remove all parentheses; deparse re-adds the ones that are needed
#' @noRd
.mdStripParen <- function(x) {
  if (is.call(x)) {
    if (identical(x[[1]], quote(`(`))) {
      return(.mdStripParen(x[[2]]))
    }
    for (.i in seq_along(x)[-1]) {
      .v <- .mdStripParen(x[[.i]])
      if (!is.null(.v)) x[[.i]] <- .v
    }
  }
  x
}

#' @noRd
.mdDeparse <- function(x) {
  paste(deparse(.mdStripParen(x), width.cutoff = 500L), collapse = " ")
}

#' Direction of the dependence of an expression on a compartment
#'
#' Rates and parameters (including symbolic exponents) are assumed positive.
#' Monotone functions (`exp`, `log`, `sqrt`, `expit`, positive powers) and
#' `delay(x, tau)` keep the direction of their argument, and a quotient whose numerator and
#' denominator both increase is taken as a saturating (Emax/Hill) increase.
#' @param expr expression
#' @param o compartment name
#' @return 1 (increasing), -1 (decreasing), 0 (independent) or NA (unknown)
#' @noRd
.mdDirection <- function(expr, o, states, deps, dependents = .mdDependents(o, states, deps)) {
  .mdDirections(expr, o, list(dependents))[[1]]
}

#' Directions of an expression in several compartments at once
#'
#' The expression is walked once; each node's direction is a vector with
#' one entry per compartment in `os` (1 increasing, -1 decreasing, 0
#' independent, NA unknown).
#' @param expr expression
#' @param os compartment names
#' @param dependents list (one per compartment) of the compartment and the
#'   variables that depend on it
#' @return numeric vector, one direction per compartment
#' @noRd
.mdDirections <- function(expr, os, dependents) {
  .k <- length(os)
  if (.k == 0L) {
    return(numeric(0))
  }
  .zero <- rep(0, .k)
  # does `e` depend on each compartment?
  .dep <- function(e) {
    .v <- all.vars(e)
    if (length(.v) == 0L) {
      return(rep(FALSE, .k))
    }
    vapply(dependents, function(d) any(.v %in% d), logical(1))
  }
  # combine the directions of the parts of a sum/product, per compartment
  .comb <- function(m) {
    if (!is.matrix(m)) {
      m <- matrix(m, nrow = .k)
    }
    apply(m, 1L, function(d) {
      d <- d[is.na(d) | d != 0]
      if (length(d) == 0L) {
        return(0)
      }
      if (anyNA(d) || length(unique(d)) > 1L) {
        return(NA_real_)
      }
      d[1]
    })
  }
  # the part of a product that depends on compartment `i` (NULL when none)
  .depPart <- function(e, i) {
    e <- .mdStripParen(e)
    if (is.call(e) && identical(e[[1]], quote(`*`)) && length(e) == 3L) {
      .a <- .depPart(e[[2]], i)
      .b <- .depPart(e[[3]], i)
      if (is.null(.a)) {
        return(.b)
      }
      if (is.null(.b)) {
        return(.a)
      }
      return(as.call(list(quote(`*`), .a, .b)))
    }
    if (.dep(e)[i]) e else NULL
  }
  # is `num/den` of the form `c*N/(K + N)` in compartment `i`, with K
  # positive and independent of it?
  .saturating <- function(e, i) {
    .n <- .depPart(e[[2]], i)
    .den <- .mdStripParen(e[[3]])
    if (is.null(.n) || !(is.call(.den) && identical(.den[[1]], quote(`+`)))) {
      return(FALSE)
    }
    .s <- list()
    .flat <- function(x) {
      x <- .mdStripParen(x)
      if (is.call(x) && identical(x[[1]], quote(`+`)) && length(x) == 3L) {
        .flat(x[[2]])
        .flat(x[[3]])
      } else {
        .s[[length(.s) + 1L]] <<- x
      }
    }
    .flat(.den)
    .isDep <- vapply(.s, function(x) .dep(x)[i], logical(1))
    .ds <- .s[.isDep]
    .ind <- .s[!.isDep]
    # K must be positive (`C/(C - 1)` decreases)
    length(.ds) == 1L &&
      length(.ind) > 0L &&
      isTRUE(all(vapply(.ind, .sgn, numeric(1)) > 0)) &&
      identical(.mdDeparse(.mdCanon(.ds[[1]])), .mdDeparse(.mdCanon(.n)))
  }
  # sign of an expression that does not depend on the compartment: parameters and
  # `exp()`/`sqrt()` are positive, constant arithmetic is evaluated, and
  # anything else is unknown (NA)
  .sgn <- function(e) {
    if (is.numeric(e) && length(e) == 1L && !is.na(e)) {
      return(sign(e))
    }
    if (is.name(e)) {
      return(1)
    }
    if (!is.call(e)) {
      return(NA_real_)
    }
    if (length(all.vars(e)) == 0L) {
      .v <- tryCatch(eval(e, baseenv()), error = function(err) NA_real_)
      if (is.numeric(.v) && length(.v) == 1L && !is.na(.v)) {
        return(sign(.v))
      }
      return(NA_real_)
    }
    .f <- e[[1]]
    if (identical(.f, quote(`(`))) {
      return(.sgn(e[[2]]))
    }
    if (identical(.f, quote(`-`)) && length(e) == 2L) {
      return(-.sgn(e[[2]]))
    }
    if (
      (identical(.f, quote(`*`)) || identical(.f, quote(`/`))) &&
        length(e) == 3L
    ) {
      return(.sgn(e[[2]]) * .sgn(e[[3]]))
    }
    if (identical(.f, quote(`+`)) && length(e) == 3L) {
      .a <- .sgn(e[[2]])
      .b <- .sgn(e[[3]])
      if (identical(.a, 1) && identical(.b, 1)) {
        return(1)
      }
      if (identical(.a, -1) && identical(.b, -1)) {
        return(-1)
      }
      return(NA_real_)
    }
    if (
      (identical(.f, quote(exp)) || identical(.f, quote(sqrt))) &&
        length(e) == 2L
    ) {
      return(1)
    }
    # a power of a positive base is positive
    if (
      (identical(.f, quote(`^`)) || identical(.f, quote(`**`))) &&
        length(e) == 3L &&
        identical(.sgn(e[[2]]), 1)
    ) {
      return(1)
    }
    NA_real_
  }
  # `a` scaled by the sign of an independent factor, where `a` is non-zero
  .scale <- function(a, e) {
    .s <- .sgn(e)
    ifelse(!is.na(a) & a == 0, 0, a * .s)
  }
  .d <- function(e) {
    if (is.name(e)) {
      .n <- as.character(e)
      .r <- ifelse(os == .n, 1, 0)
      # a variable that was not substituted (conditional/reassigned)
      .r[.r == 0 & .dep(e)] <- NA_real_
      return(.r)
    }
    if (!is.call(e)) {
      return(.zero)
    }
    .f <- e[[1]]
    .fn <- if (is.name(.f)) as.character(.f) else ""
    if (.fn == "(") {
      return(.d(e[[2]]))
    }
    if (.fn == "+") {
      return(.comb(vapply(as.list(e)[-1], .d, numeric(.k))))
    }
    if (.fn == "-") {
      if (length(e) == 2L) {
        return(-.d(e[[2]]))
      }
      return(.comb(cbind(.d(e[[2]]), -.d(e[[3]]))))
    }
    if (.fn == "*") {
      .a <- .d(e[[2]])
      .b <- .d(e[[3]])
      # a factor that does not depend on the compartment may still carry a
      # sign (R parses `-k*C` as `(-k)*C`)
      .a0 <- !is.na(.a) & .a == 0
      .b0 <- !is.na(.b) & .b == 0
      .r <- .comb(cbind(.a, .b))
      .r[.a0 & !.b0] <- .scale(.b, e[[2]])[.a0 & !.b0]
      .r[.b0 & !.a0] <- .scale(.a, e[[3]])[.b0 & !.a0]
      .r[.a0 & .b0] <- 0
      return(.r)
    }
    if (.fn == "/") {
      .n <- .d(e[[2]])
      .m <- .d(e[[3]])
      .n0 <- !is.na(.n) & .n == 0
      .m0 <- !is.na(.m) & .m == 0
      .r <- .comb(cbind(.n, -.m))
      .r[.m0 & !.n0] <- .scale(.n, e[[3]])[.m0 & !.n0]
      .r[.n0 & !.m0] <- .scale(-.m, e[[2]])[.n0 & !.m0]
      .r[.n0 & .m0] <- 0
      # saturating forms `N/(K + N)` (Emax/Hill: `C^g/(ec50^g + C^g)`)
      # increase with C; other quotients may not be monotone
      for (.i in which(!is.na(.n) & .n == 1 & !is.na(.m) & .m == 1)) {
        if (.saturating(e, .i)) .r[.i] <- 1
      }
      return(.r)
    }
    if (.fn %in% c("exp", "log", "sqrt", "expit", "log1p", "log10", "log2") && length(e) == 2L) {
      return(.d(e[[2]]))
    }
    # a delayed value moves with the value itself (the delay time is a
    # parameter)
    if (.fn == "delay" && length(e) >= 2L) {
      .r <- .d(e[[2]])
      for (.a in as.list(e)[-(1:2)]) {
        .r[.dep(.a)] <- NA_real_
      }
      return(.r)
    }
    if (.fn %in% c("^", "**") && length(e) == 3L) {
      .b <- .d(e[[2]])
      # `(A - 1)^2` is not monotone in A: the base must be non-negative
      if (!.mdNonNeg(e[[2]])) {
        .b[is.na(.b) | .b != 0] <- NA_real_
      }
      .p <- .mdStripParen(e[[3]])
      .r <- if (is.numeric(.p)) {
        sign(.p) * .b
      } else if (is.call(.p) && identical(.p[[1]], quote(`-`)) && length(.p) == 2L) {
        # `C^(-gamma)`: a negated exponent decreases
        -.sgn(.p[[2]]) * .b
      } else {
        # a symbolic exponent (e.g. a Hill coefficient) is assumed positive
        .b
      }
      .r[.dep(e[[3]])] <- NA_real_
      return(.r)
    }
    ifelse(.dep(e), NA_real_, 0)
  }
  .d(expr)
}

#' A compartment and the variables that depend on it
#' @noRd
.mdDependents <- function(o, states, closure) {
  .v <- names(closure)[vapply(closure, function(s) o %in% s, logical(1))]
  c(o, setdiff(.v, states))
}

#' Is an expression non-negative wherever it is defined?
#'
#' Parameters and compartment amounts are assumed non-negative; a difference
#' (or a negated value) may be of either sign.
#' @noRd
.mdNonNeg <- function(e) {
  e <- .mdStripParen(e)
  if (is.numeric(e)) {
    return(all(!is.na(e) & e >= 0))
  }
  if (is.name(e)) {
    return(TRUE)
  }
  if (!is.call(e)) {
    return(FALSE)
  }
  .f <- e[[1]]
  if (identical(.f, quote(`+`)) || identical(.f, quote(`*`)) || identical(.f, quote(`/`))) {
    return(all(vapply(as.list(e)[-1], .mdNonNeg, logical(1))))
  }
  if (
    (identical(.f, quote(exp)) || identical(.f, quote(sqrt))) &&
      length(e) == 2L
  ) {
    return(TRUE)
  }
  if (
    (identical(.f, quote(`^`)) || identical(.f, quote(`**`))) &&
      length(e) == 3L
  ) {
    return(.mdNonNeg(e[[2]]))
  }
  FALSE
}

#' Interaction sign from a direction (unknown directions are 0)
#' @noRd
.mdSign <- function(d) {
  d <- unname(d)
  if (length(d) == 0L || is.na(d)) 0 else d
}

#' Numerator and denominator factors of a product/quotient term
#' @return list(num = list of expressions, den = list of expressions)
#' @noRd
.mdTermFactors <- function(x) {
  .env <- new.env(parent = emptyenv())
  .env$num <- list()
  .env$den <- list()
  .flat <- function(e, num) {
    if (is.call(e) && identical(e[[1]], quote(`(`))) {
      return(.flat(e[[2]], num))
    }
    if (is.call(e) && length(e) == 3L && identical(e[[1]], quote(`*`))) {
      .flat(e[[2]], num)
      .flat(e[[3]], num)
    } else if (is.call(e) && length(e) == 3L && identical(e[[1]], quote(`/`))) {
      .flat(e[[2]], num)
      .flat(e[[3]], !num)
    } else if (num) {
      .env$num[[length(.env$num) + 1L]] <- e
    } else {
      .env$den[[length(.env$den) + 1L]] <- e
    }
  }
  .flat(x, TRUE)
  list(num = .env$num, den = .env$den)
}

#' Canonical key of a product/quotient term (factor order ignored)
#' @noRd
.mdTermKey <- function(x) {
  .f <- .mdTermFactors(x)
  .num <- vapply(.f$num, function(e) .mdDeparse(.mdCanon(e)), character(1))
  .den <- vapply(.f$den, function(e) .mdDeparse(.mdCanon(e)), character(1))
  paste0(paste(sort(.num), collapse = "*"), "/", paste(sort(.den), collapse = "*"))
}

# ---------------------------------------------------------------------------
# Classification
# ---------------------------------------------------------------------------

#' Classify ODE terms into graph edges
#' @param terms data frame of terms from `.mdModelInfo()`
#' @param states compartment names
#' @return edge data frame
#' @noRd
.mdClassifyTerms <- function(terms, states) {
  .empty <- data.frame(
    from = character(0),
    to = character(0),
    type = character(0),
    sign = numeric(0),
    label = character(0),
    stringsAsFactors = FALSE
  )
  if (is.null(terms) || nrow(terms) == 0L) {
    .empty$bidirectional <- logical(0)
    return(.empty)
  }
  .n <- nrow(terms)
  .used <- rep(FALSE, .n)
  .matchedFrom <- vector("list", .n)
  .rows <- list()
  .add <- function(from, to, type, sign, label) {
    .rows[[length(.rows) + 1L]] <<-
      data.frame(from = from, to = to, type = type, sign = sign, label = label, stringsAsFactors = FALSE)
  }
  # mass transfer: -term in the source (containing the source amount) and
  # +term in another compartment.  One +term may receive mass from several
  # sources (e.g. binding `kon*C*R` into a complex from both C and R) ...
  for (.i in seq_len(.n)) {
    if (.used[.i] || terms$sign[.i] > 0) {
      next
    }
    .src <- terms$state[.i]
    # a term subtracted from one compartment and added to another is
    # conserved flow, whatever drives it (first-order `k*A`, zero-order
    # `rate`, or another compartment like an enzyme `Vmax*E`)
    .j <- which(
      terms$sign > 0 &
        terms$state != .src &
        terms$key == terms$key[.i] &
        !vapply(.matchedFrom, function(m) .src %in% m, logical(1))
    )
    if (length(.j) == 0L) {
      next
    }
    # ... and one -term may go to several destinations (e.g. dissociation
    # `koff*RC` back to both C and R); keep one destination per compartment
    .j <- .j[!duplicated(terms$state[.j])]
    if (!(.src %in% terms$states[[.i]])) {
      # a rate that does not depend on the source (zero-order or driven by
      # another compartment) is one flow: pair it with one unused +term
      .j <- .j[!.used[.j]][1]
      if (is.na(.j)) next
    }
    .used[c(.i, .j)] <- TRUE
    for (.k in .j) {
      .matchedFrom[[.k]] <- c(.matchedFrom[[.k]], .src)
      .add(.src, terms$state[.k], "transfer", 1, terms$label[.k])
    }
  }
  # compartments that drive a transfer without being its source (e.g. an
  # enzyme `E` in `Vmax*E*A`) stimulate/inhibit the destination
  for (.j in which(!vapply(.matchedFrom, is.null, logical(1)))) {
    .drivers <- setdiff(terms$states[[.j]], c(.matchedFrom[[.j]], terms$state[.j]))
    for (.o in .drivers) {
      .add(.o, terms$state[.j], "interaction", .mdSign(terms$dir[[.j]][.o]), terms$label[.j])
    }
  }
  for (.i in which(!.used)) {
    .s <- terms$state[.i]
    .st <- terms$states[[.i]]
    .others <- setdiff(.st, .s)
    # direction of the effect of each other compartment on d/dt(.s)
    .dir <- function(o) {
      terms$sign[.i] * .mdSign(terms$dir[[.i]][o])
    }
    if (terms$sign[.i] < 0) {
      if (.s %in% .st || length(.others) == 0L) {
        .add(.s, NA_character_, "elimination", -1, terms$label[.i])
      }
    } else if (length(.others) == 0L) {
      .add(NA_character_, .s, "input", 1, terms$label[.i])
    }
    for (.o in .others) {
      .add(.o, .s, "interaction", .dir(.o), terms$label[.i])
    }
  }
  .e <- do.call(rbind, .rows)
  if (is.null(.e)) {
    .e <- .empty
  }
  # combine duplicated flows
  if (nrow(.e) > 0L) {
    .id <- paste(.e$from, .e$to, .e$type, .e$sign, sep = "\r")
    .e <- do.call(
      rbind,
      lapply(unique(.id), function(.k) {
        .w <- .e[.id == .k, , drop = FALSE]
        # keep every contribution (`k*A + k*A` is twice `k*A`)
        .w$label[1] <- paste(.w$label, collapse = " + ")
        .w[1, , drop = FALSE]
      })
    )
  }
  .tr <- .e$type == "transfer"
  .pairs <- paste(.e$from, .e$to, sep = "\r")
  .rev <- paste(.e$to, .e$from, sep = "\r")
  .e$bidirectional <- .tr & .rev %in% .pairs[.tr]
  rownames(.e) <- NULL
  .e
}

# ---------------------------------------------------------------------------
# Layout
# ---------------------------------------------------------------------------

#' Choose the central compartment
#' @noRd
.mdCentral <- function(states, edges) {
  .cn <- states[tolower(states) %in% c("central", "center", "centr", "cent")]
  if (length(.cn) > 0L) {
    return(.cn[1])
  }
  .tr <- edges[edges$type == "transfer", , drop = FALSE]
  .deg <- vapply(states, function(.s) sum(.tr$from == .s | .tr$to == .s), numeric(1))
  .elim <- states %in% edges$from[edges$type == "elimination"]
  .score <- .deg + 0.5 * .elim
  if (all(.deg == 0)) {
    return(states[1])
  }
  states[which.max(.score)]
}

#' Lay out the compartments on a grid
#'
#' @return node data frame with name, role, dosing, x, y
#' @noRd
.mdLayout <- function(states, edges, dosing) {
  .central <- .mdCentral(states, edges)
  .tr <- edges[edges$type == "transfer", , drop = FALSE]
  .bi <- .tr[.tr$bidirectional, , drop = FALSE]
  .uni <- .tr[!.tr$bidirectional, , drop = FALSE]
  .int <- edges[edges$type == "interaction", , drop = FALSE]
  # every arrow drawn between two compartments
  .links <- edges[!is.na(edges$from) & !is.na(edges$to), , drop = FALSE]
  .x <- stats::setNames(rep(NA_real_, length(states)), states)
  .y <- .x
  .role <- stats::setNames(rep("other", length(states)), states)
  .free <- function(x, y) {
    !any(!is.na(.x) & abs(.x - x) < 0.9 & abs(.y - y) < 0.9)
  }
  # is (x, y) free for `s`, with the arrows between `s` and the placed
  # compartments clear of other compartments, and no placed arrow running
  # through the new box?
  # arrow partners of each compartment and edge end point indices, computed
  # once (large QSP/PBPK models have thousands of edges)
  .fromI <- match(.links$from, states)
  .toI <- match(.links$to, states)
  .adj <- lapply(seq_along(states), function(.i) {
    unique(c(.fromI[.toI == .i], .toI[.fromI == .i]))
  })
  names(.adj) <- states
  .ok <- function(s, x, y) {
    if (!.free(x, y)) {
      return(FALSE)
    }
    .isPlaced <- !is.na(.x)
    # arrows between `s` and placed partners must be clear of the other
    # placed compartments
    .p <- .adj[[s]]
    .p <- .p[.isPlaced[.p]]
    if (length(.p) > 0L) {
      .others <- which(.isPlaced)
      .g <- expand.grid(p = .p, o = .others)
      .g <- .g[.g$p != .g$o, , drop = FALSE]
      if (
        nrow(.g) > 0L &&
          any(.mdSegRect(.x[.g$p], .y[.g$p], x, y, .x[.g$o], .y[.g$o]))
      ) {
        return(FALSE)
      }
    }
    # no placed arrow may run through the new box
    .pe <- .isPlaced[.fromI] & .isPlaced[.toI]
    if (
      any(.pe) &&
        any(.mdSegRect(.x[.fromI[.pe]], .y[.fromI[.pe]], .x[.toI[.pe]], .y[.toI[.pe]], x, y))
    ) {
      return(FALSE)
    }
    TRUE
  }
  .place <- function(s, x, y, step) {
    # move away (in the `step` direction) until the position is clear,
    # trying neighboring columns in each row
    # large models search fewer rows before falling back
    # a wide search gives cleaner diagrams; very large models use a
    # bounded one so the layout stays fast
    .wide <- length(states) <= .mdMaxWideSearch
    .kmax <- if (.wide) 2L * length(states) + 4L else 12L
    for (.k in 0:.kmax) {
      for (.dx in if (.wide) c(0, 1, -1, 2, -2, 3, -3, 4, -4) else c(0, 1, -1, 2, -2)) {
        if (.ok(s, x + .dx, y + .k * step)) {
          .x[s] <<- x + .dx
          .y[s] <<- y + .k * step
          return(invisible())
        }
      }
    }
    # give up on arrow clearance rather than fail
    while (!.free(x, y)) {
      y <- y + step
    }
    .x[s] <<- x
    .y[s] <<- y
  }
  # breadth-first placement of the (mass transfer) neighbors of placed nodes
  # `side` is where exchange (bidirectional) partners go: left of the PK
  # model, right of PD models so they stay off the arrows coming in from the
  # left
  # place a group of neighbors of `s` in direction `dir` (radians: up is
  # pi/2, down -pi/2, left pi, right 0).  A few neighbors are placed in the
  # usual line; many (a PBPK hub with its tissues) are fanned out on an arc
  # whose radius grows with their number, so the arrows radiate out instead
  # of grazing the boxes in between
  # crossings involving the compartments in `new` (the rest of the diagram
  # is the same whichever arrangement is chosen, so they need not be
  # counted)
  .nCrossings <- function(new) {
    .p <- !is.na(.x)
    .e <- which(.p[.fromI] & .p[.toI])
    .o <- which(.p)
    if (length(.e) == 0L || length(.o) == 0L) {
      return(0)
    }
    .f <- .fromI[.e]
    .t <- .toI[.e]
    .g <- expand.grid(e = seq_along(.e), o = .o)
    .keep <- .f[.g$e] != .g$o & .t[.g$e] != .g$o & (.f[.g$e] %in% new | .t[.g$e] %in% new | .g$o %in% new)
    .g <- .g[.keep, , drop = FALSE]
    if (nrow(.g) == 0L) {
      return(0)
    }
    sum(.mdSegRect(.x[.f[.g$e]], .y[.f[.g$e]], .x[.t[.g$e]], .y[.t[.g$e]], .x[.g$o], .y[.g$o]))
  }
  # place a group of neighbors of `s` in direction `dir` (radians: up is
  # pi/2, down -pi/2, left pi, right 0), either in a line or -- for a hub
  # with many neighbors, like a PBPK central compartment -- fanned out on an
  # arc whose radius grows with their number.  Both are tried and the one
  # with fewer crossings is kept.
  .placeGroup <- function(nodes, s, dir, step) {
    .m <- length(nodes)
    if (.m == 0L) {
      return(invisible())
    }
    .line <- function() {
      .off <- (seq_len(.m) - 1) - if (abs(cos(dir)) > 0.5) 0 else (.m - 1) / 2
      for (.k in seq_len(.m)) {
        .dx <- if (abs(cos(dir)) > 0.5) round(cos(dir)) else .off[.k]
        .dy <- if (abs(cos(dir)) > 0.5) .off[.k] else round(sin(dir))
        .place(nodes[.k], .x[s] + .dx, .y[s] + .dy, step)
      }
    }
    .fan <- function() {
      .sector <- pi * 2 / 3
      .d <- .sector / (.m - 1)
      .r <- max(1.5, 1.3 / .d)
      .a <- dir - .sector / 2 + (seq_len(.m) - 1) * .d
      for (.k in seq_len(.m)) {
        .place(nodes[.k], round((.x[s] + .r * cos(.a[.k])) * 2) / 2, round((.y[s] + .r * sin(.a[.k])) * 2) / 2, step)
      }
    }
    if (.m <= 2L) {
      .line()
      return(invisible())
    }
    .new <- match(nodes, states)
    .x0 <- .x
    .y0 <- .y
    .line()
    .nLine <- .nCrossings(.new)
    .xLine <- .x
    .yLine <- .y
    .x <<- .x0
    .y <<- .y0
    .fan()
    if (.nCrossings(.new) > .nLine) {
      .x <<- .xLine
      .y <<- .yLine
    }
    invisible()
  }

  .spread <- function(start, side = -1) {
    .queue <- start
    while (length(.queue) > 0L) {
      .s <- .queue[1]
      .queue <- .queue[-1]
      .isCentral <- .s == .central
      # upstream (unidirectional into .s): above
      .up <- setdiff(unique(.uni$from[.uni$to == .s]), names(.x)[!is.na(.x)])
      .placeGroup(.up, .s, pi / 2, 1)
      for (.n in .up) {
        if (.role[.n] == "other") .role[.n] <<- "transit"
      }
      # bidirectional exchange: to the left (PD models: right), fanned
      .lr <- unique(c(.bi$to[.bi$from == .s], .bi$from[.bi$to == .s]))
      .lr <- setdiff(.lr, names(.x)[!is.na(.x)])
      .placeGroup(.lr, .s, if (side < 0) pi else 0, -1)
      for (.n in .lr) {
        if (.role[.n] == "other") .role[.n] <<- "peripheral"
      }
      # downstream (unidirectional out of .s): below
      .dn <- setdiff(unique(.uni$to[.uni$from == .s]), names(.x)[!is.na(.x)])
      .placeGroup(.dn, .s, -pi / 2, -1)
      for (.n in .dn) {
        if (.role[.n] == "other") {
          .role[.n] <<- if (.isCentral) "metabolite" else "other"
        }
      }
      .queue <- c(.queue, .up, .lr, .dn)
    }
  }

  .role[.central] <- "central"
  .place(.central, 0, 0, -1)
  .spread(.central)
  # choose a row in column `nx` for `s` where the straight interaction
  # arrows between `s` and the already placed compartments (either
  # direction) are clear of the other compartments
  .pickRow <- function(s, nx, y0) {
    .m <- if (length(states) <= .mdMaxWideSearch) length(states) else 12L
    for (.d in c(0, rbind(-seq_len(.m), seq_len(.m)))) {
      .cy <- y0 + .d
      if (.ok(s, nx, .cy)) return(.cy)
    }
    y0
  }
  # compartments interacting without mass transfer: to the right
  repeat {
    .placed <- names(.x)[!is.na(.x)]
    .cand <- .int[.int$from %in% .placed & !(.int$to %in% .placed), , drop = FALSE]
    if (nrow(.cand) > 0L) {
      .s <- .cand$to[1]
      .nx <- max(.x, na.rm = TRUE) + 1
      .place(.s, .nx, .pickRow(.s, .nx, .y[.cand$from[1]]), -1)
      .role[.s] <- "effect"
      .spread(.s, side = 1)
      next
    }
    .left <- states[is.na(.x)]
    if (length(.left) == 0L) {
      break
    }
    # compartments acting on the placed ones first, then dosing compartments
    .acting <- intersect(.left, .int$from[.int$to %in% .placed])
    .s <- c(.acting, intersect(dosing, .left), .left)[1]
    .nx <- max(.x, na.rm = TRUE) + 1
    .place(.s, .nx, .pickRow(.s, .nx, 0), -1)
    .spread(.s, side = if (length(.acting) > 0L) 1 else -1)
  }
  .role[states %in% dosing & .role != "central"] <- "dosing"
  data.frame(
    name = states,
    role = unname(.role[states]),
    dosing = states %in% dosing,
    x = unname(.x[states]),
    y = unname(.y[states]),
    stringsAsFactors = FALSE
  )
}

#' Does the segment (x0, y0)-(x1, y1) pass through any of the node boxes?
#'
#' Boxes are centered on (x, y) (grid units) with half width `hw` and half
#' height `hh`; the segment is sampled finely enough for unit-grid layouts.
#' @noRd
.mdSegmentCrosses <- function(x0, y0, x1, y1, x, y, hw = 0.35, hh = 0.25) {
  if (length(x) == 0L) {
    return(FALSE)
  }
  any(.mdSegRect(x0, y0, x1, y1, x, y, hw, hh))
}

#' Exact segment / open rectangle intersection (Liang-Barsky), vectorized
#'
#' @param x0,y0,x1,y1 segment end points (recycled)
#' @param cx,cy rectangle centers (recycled)
#' @param hw,hh rectangle half width and half height
#' @return logical vector: does each segment pass through the interior of
#'   its rectangle?
#' @noRd
.mdSegRect <- function(x0, y0, x1, y1, cx, cy, hw = 0.35, hh = 0.25) {
  .n <- max(length(x0), length(cx))
  x0 <- rep_len(x0, .n)
  y0 <- rep_len(y0, .n)
  x1 <- rep_len(x1, .n)
  y1 <- rep_len(y1, .n)
  cx <- rep_len(cx, .n)
  cy <- rep_len(cy, .n)
  .lo <- rep(0, .n)
  .hi <- rep(1, .n)
  .ok <- rep(TRUE, .n)
  .slab <- function(p0, d, c, h) {
    .z <- abs(d) < 1e-12
    # parallel to the slab: inside only when strictly within it
    .ok <<- .ok & (!.z | abs(p0 - c) < h)
    .t1 <- ifelse(.z, -Inf, (c - h - p0) / d)
    .t2 <- ifelse(.z, Inf, (c + h - p0) / d)
    .lo <<- pmax(.lo, pmin(.t1, .t2))
    .hi <<- pmin(.hi, pmax(.t1, .t2))
  }
  .slab(x0, x1 - x0, cx, hw)
  .slab(y0, y1 - y0, cy, hh)
  .ok & .lo < .hi
}

#' Positions of the invisible input/output end points
#'
#' @return edge data frame with added columns x0, y0, x1, y1 (in grid units)
#' @noRd
.mdEdgeCoords <- function(graph) {
  .n <- graph$nodes
  .e <- graph$edges
  .px <- stats::setNames(.n$x, .n$name)
  .py <- stats::setNames(.n$y, .n$name)
  .occupied <- function(x, y) any(abs(.n$x - x) < 0.5 & abs(.n$y - y) < 0.6)
  .e$x0 <- .px[.e$from]
  .e$y0 <- .py[.e$from]
  .e$x1 <- .px[.e$to]
  .e$y1 <- .py[.e$to]
  for (.i in seq_len(nrow(.e))) {
    if (.e$type[.i] == "elimination") {
      .x <- .e$x0[.i]
      .y <- .e$y0[.i]
      .e$x1[.i] <- .x
      .e$y1[.i] <- .y - 0.7
      if (.occupied(.x, .y - 1)) {
        .e$x1[.i] <- .x + 0.45
        .e$y1[.i] <- .y - 0.6
      }
    } else if (.e$type[.i] == "input") {
      .x <- .e$x1[.i]
      .y <- .e$y1[.i]
      .e$x0[.i] <- .x
      .e$y0[.i] <- .y + 0.7
      if (.occupied(.x, .y + 1)) {
        .e$x0[.i] <- .x + 0.45
        .e$y0[.i] <- .y + 0.6
      }
    }
  }
  .e
}

.mdRoleColors <- c(
  dosing = "#F2C57C",
  central = "#7FB3D5",
  peripheral = "#A9CCE3",
  transit = "#FAD7A0",
  metabolite = "#D2B4DE",
  effect = "#A9DFBF",
  other = "#E5E7E9"
)

# ---------------------------------------------------------------------------
# Engines
# ---------------------------------------------------------------------------

#' Graphviz DOT source for a model graph
#' @noRd
.mdDot <- function(graph, labels = FALSE) {
  # column spacing (inches) wide enough for the longest compartment name
  # (14pt Helvetica is about 0.11 inch per character plus margins)
  .xs <- max(1.6, 0.11 * max(nchar(graph$nodes$name), 0L) + 0.9)
  .ys <- 1.1
  # dosing properties shown as an external label next to the compartment
  .ann <- graph$nodes$annotation
  if (is.null(.ann)) {
    .ann <- rep("", nrow(graph$nodes))
  }
  .ann[is.na(.ann)] <- ""
  # quote a DOT string; "\r" (from combined labels) becomes a DOT line break
  .q <- function(x) {
    paste0("\"", gsub("\r", "\\n", gsub("\"", "\\\\\"", x), fixed = TRUE), "\"")
  }
  .n <- graph$nodes
  .e <- .mdEdgeCoords(graph)
  # routing curved edges around the boxes is slow for large graphs, which
  # do not gain much from it (the positions are pinned)
  .splines <- if (nrow(graph$edges) > .mdMaxSplineEdges) "line" else "true"
  .lines <- c(
    "digraph model {",
    sprintf(
      paste0("  graph [layout = neato, splines = %s, ", "outputorder = edgesfirst, forcelabels = true];"),
      .splines
    ),
    "  node [shape = box, style = \"rounded,filled\", fontname = Helvetica];",
    "  edge [fontname = Helvetica, fontsize = 10];"
  )
  for (.i in seq_len(nrow(.n))) {
    .lines <- c(
      .lines,
      sprintf(
        "  %s [pos = \"%g,%g!\", fillcolor = %s%s%s];",
        .q(.n$name[.i]),
        .n$x[.i] * .xs,
        .n$y[.i] * .ys,
        .q(.mdRoleColors[[.n$role[.i]]]),
        if (.n$dosing[.i]) ", penwidth = 2" else "",
        if (nzchar(.ann[.i])) paste0(", xlabel = ", .q(gsub("\n", "\r", .ann[.i], fixed = TRUE))) else ""
      )
    )
  }
  .done <- rep(FALSE, nrow(.e))
  for (.i in seq_len(nrow(.e))) {
    if (.done[.i]) {
      next
    }
    .t <- .e$type[.i]
    .from <- .e$from[.i]
    .to <- .e$to[.i]
    .attr <- character(0)
    .lab <- .e$label[.i]
    if (.t %in% c("elimination", "input")) {
      .pt <- paste0(".", .t, .i)
      .lines <- c(
        .lines,
        sprintf(
          "  %s [shape = point, style = invis, width = 0.01, pos = \"%g,%g!\"];",
          .q(.pt),
          (if (.t == "elimination") .e$x1[.i] else .e$x0[.i]) * .xs,
          (if (.t == "elimination") .e$y1[.i] else .e$y0[.i]) * .ys
        )
      )
      if (.t == "elimination") .to <- .pt else .from <- .pt
    } else if (.t == "transfer" && .e$bidirectional[.i]) {
      .j <- which(.e$type == "transfer" & .e$from == .to & .e$to == .from)
      .done[.j] <- TRUE
      .attr <- c(.attr, "dir = both")
      .lab <- paste(c(.lab, .e$label[.j]), collapse = "\r")
    } else if (.t == "interaction") {
      .attr <- c(
        .attr,
        "style = dashed",
        "color = gray40",
        if (.e$sign[.i] < 0) "arrowhead = tee",
        if (.e$sign[.i] == 0) "arrowhead = dot"
      )
    }
    if (labels) {
      .attr <- c(.attr, paste0("label = ", .q(.lab)))
    }
    .lines <- c(
      .lines,
      sprintf(
        "  %s -> %s%s;",
        .q(.from),
        .q(.to),
        if (length(.attr) > 0L) paste0(" [", paste(.attr, collapse = ", "), "]") else ""
      )
    )
    .done[.i] <- TRUE
  }
  paste(c(.lines, "}"), collapse = "\n")
}

#' @noRd
.mdDiagrammeR <- function(graph, labels = FALSE) {
  if (!requireNamespace("DiagrammeR", quietly = TRUE)) {
    stop(
      "the 'DiagrammeR' engine requires the 'DiagrammeR' package; install it or use engine = \"ggplot2\"",
      call. = FALSE
    )
  }
  DiagrammeR::grViz(.mdDot(graph, labels))
}

#' Move a segment's end points to the boundaries of the node boxes
#' @noRd
.mdClip <- function(x0, y0, x1, y1, hw, hh, clip0, clip1) {
  .dx <- x1 - x0
  .dy <- y1 - y0
  .t <- function(dx, dy) {
    .tx <- if (abs(dx) > 1e-8) hw / abs(dx) else Inf
    .ty <- if (abs(dy) > 1e-8) hh / abs(dy) else Inf
    min(.tx, .ty, 0.45)
  }
  .t0 <- if (clip0) .t(.dx, .dy) else 0
  .t1 <- if (clip1) .t(.dx, .dy) else 0
  c(x0 + .t0 * .dx, y0 + .t0 * .dy, x1 - .t1 * .dx, y1 - .t1 * .dy)
}

#' @noRd
.mdGgplot <- function(graph, labels = FALSE) {
  .xs <- 1.6
  .hw <- 0.5
  .hh <- 0.2
  .n <- graph$nodes
  .n$x <- .n$x * .xs
  .e <- .mdEdgeCoords(graph)
  .e$x0 <- .e$x0 * .xs
  .e$x1 <- .e$x1 * .xs
  .px <- stats::setNames(.n$x, .n$name)
  .py <- stats::setNames(.n$y, .n$name)
  # arrows between the same two compartments (in either direction)
  .pair <- ifelse(
    is.na(.e$from) | is.na(.e$to),
    paste0(".", seq_len(nrow(.e))),
    paste(pmin(.e$from, .e$to), pmax(.e$from, .e$to), sep = "\r")
  )
  .pairN <- as.integer(stats::ave(seq_along(.pair), .pair, FUN = length))
  .pairK <- as.integer(stats::ave(seq_along(.pair), .pair, FUN = seq_along))
  .seg <- do.call(
    rbind,
    lapply(seq_len(nrow(.e)), function(.i) {
      .x0 <- .e$x0[.i]
      .y0 <- .e$y0[.i]
      .x1 <- .e$x1[.i]
      .y1 <- .e$y1[.i]
      if (.e$type[.i] %in% c("transfer", "interaction") && .e$from[.i] == .e$to[.i]) {
        return(NULL)
      }
      if (.pairN[.i] > 1L) {
        # spread several arrows between the same two compartments (both
        # directions of an exchange, or a stimulation and an inhibition) to
        # either side of the center line
        .a <- c(pmin(.e$from[.i], .e$to[.i]), pmax(.e$from[.i], .e$to[.i]))
        .dx <- .px[.a[2]] - .px[.a[1]]
        .dy <- .py[.a[2]] - .py[.a[1]]
        .len <- sqrt(.dx^2 + .dy^2)
        .off <- (.pairK[.i] - (.pairN[.i] + 1) / 2) * 0.12
        .ox <- -.dy / .len * .off
        .oy <- .dx / .len * .off
        .x0 <- .x0 + .ox
        .x1 <- .x1 + .ox
        .y0 <- .y0 + .oy
        .y1 <- .y1 + .oy
      }
      .c <- .mdClip(.x0, .y0, .x1, .y1, .hw, .hh, clip0 = .e$type[.i] != "input", clip1 = .e$type[.i] != "elimination")
      data.frame(
        x = .c[1],
        y = .c[2],
        xend = .c[3],
        yend = .c[4],
        flow = ifelse(
          .e$type[.i] == "interaction",
          c("inhibition", "modulation", "stimulation")[sign(.e$sign[.i]) + 2],
          "mass transfer"
        ),
        label = .e$label[.i],
        stringsAsFactors = FALSE
      )
    })
  )
  .n$role <- factor(.n$role, levels = names(.mdRoleColors))
  if (is.null(.n$annotation)) {
    .n$annotation <- ""
  }
  .n$annotation[is.na(.n$annotation)] <- ""
  .ann <- .n[nzchar(.n$annotation), , drop = FALSE]
  .ann$x <- .ann$x + .hw * 0.9
  .ann$y <- .ann$y + .hh * 1.1
  .annLines <- strsplit(.ann$annotation, "\n", fixed = TRUE)
  .annExtent <- data.frame(
    x = .ann$x + 0.075 * vapply(.annLines, function(l) max(nchar(l)), numeric(1)),
    y = .ann$y + 0.13 * lengths(.annLines)
  )
  .p <- ggplot2::ggplot() +
    ggplot2::geom_tile(
      data = .n,
      ggplot2::aes(x = .data$x, y = .data$y, fill = .data$role),
      width = 2 * .hw,
      height = 2 * .hh,
      color = "gray30",
      linewidth = ifelse(.n$dosing, 1, 0.4)
    ) +
    ggplot2::geom_text(data = .n, ggplot2::aes(x = .data$x, y = .data$y, label = .data$name)) +
    # dosing properties (lag, F, rate, dur) as an annotation at the upper
    # right corner of the compartment
    ggplot2::geom_text(
      data = .ann,
      ggplot2::aes(x = .data$x, y = .data$y, label = .data$annotation),
      hjust = 0,
      vjust = 0,
      size = 2.6,
      fontface = "italic",
      lineheight = 0.9
    ) +
    # keep the annotations inside the plot (away from the legend)
    ggplot2::geom_blank(data = .annExtent, ggplot2::aes(x = .data$x, y = .data$y)) +
    ggplot2::scale_fill_manual(values = .mdRoleColors, drop = TRUE, name = "compartment") +
    ggplot2::coord_equal(clip = "off") +
    ggplot2::theme_void() +
    ggplot2::theme(plot.margin = ggplot2::margin(10, 10, 10, 10))
  if (!is.null(.seg) && nrow(.seg) > 0L) {
    .p <- .p +
      ggplot2::geom_segment(
        data = .seg,
        ggplot2::aes(x = .data$x, y = .data$y, xend = .data$xend, yend = .data$yend, linetype = .data$flow),
        arrow = ggplot2::arrow(length = ggplot2::unit(0.08, "inches"), type = "closed")
      ) +
      ggplot2::scale_linetype_manual(
        values = c("mass transfer" = "solid", stimulation = "dashed", inhibition = "dotted", modulation = "dotdash"),
        name = "flow"
      )
    if (labels) {
      .p <- .p +
        ggplot2::geom_label(
          data = .seg,
          ggplot2::aes(x = (.data$x + .data$xend) / 2, y = (.data$y + .data$yend) / 2, label = .data$label),
          size = 2.5
        )
    }
  }
  .p
}