R/parse_pml_observes.R

In Certara.Xpose.NLME: Enhances 'xpose' Diagnostics for Pharmacometric Models from 'Certara.RsNLME' and Phoenix NLME

# Lightweight PML parsing scoped to the two needs of get_summaryNlme():
#  (1) recover the ObsName -> sigma map from xp$code so per-sigma IWRES
#      pooling can be reproduced for shrinkage = "sd" / "var";
#  (2) classify each sigma's role in its observe() expression (additive /
#      proportional / other / mixed / unknown) for the residual-transform
#      advisories.
#
# Neither of these can assume anything about which tool produced the PML --
# it may have been generated by a model-building helper, hand-written from
# scratch, or edited afterwards. Everything here works only on the text
# actually stored in `xp$code`.
#
# Block/token extraction (1) uses a recursive PCRE for balanced parens and
# possessive quantifiers to avoid catastrophic backtracking. Deliberately
# scoped: no full PML statement classifier here.

# Strip line and block comments and normalise newlines without flattening
# the source -- preserving newlines avoids accidental token glue across
# adjacent statements.
.strip_pml_comments <- function(src) {
  pattern <- paste(
    "(?:\\/\\/(?:\\\\\\n|[^\\n])*(?=$|\\n))",   # //...
    "(?:#(?:\\\\\\n|[^\\n])*(?=$|\\n))",         # #...
    "(?:\\/\\*[\\s\\S]*?\\*\\/)",                # /* ... */
    sep = "|"
  )
  gsub(pattern, "", src, perl = TRUE)
}

# Return the inner content (sans the outer parens) of every `keyword(...)`
# block in `src`. The recursive subroutine `(?-1)` matches arbitrarily nested
# parentheses; possessive quantifiers prevent catastrophic backtracking.
.match_keyword_blocks <- function(src, keyword) {
  pattern <- paste0("(?<=\\b", keyword, ")(\\((?:[^()]++|(?-1))*+\\))")
  hits <- regmatches(src, gregexpr(pattern, src, perl = TRUE))[[1]]
  if (!length(hits)) return(character())
  substr(hits, 2L, nchar(hits) - 1L)
}

# Identifier tokens (words starting with a letter or underscore), preserving
# the order of first appearance.
.identifier_tokens <- function(s) {
  toks <- regmatches(s, gregexpr("[A-Za-z_][A-Za-z0-9_]*", s, perl = TRUE))[[1]]
  toks[!duplicated(toks)]
}

# --------------------------------------------------------------------------
# Symbolic sigma-role classification
# --------------------------------------------------------------------------

# Pull the right-hand side out of an `observe(...)` block, which may arrive
# either as the bare expression or as the full `ObsName = expression` /
# `ObsName(arg) = expression` text (the latter is how the engine dumps a
# table-column observe such as `EObs(C) = E * (1 + EEps)`). Parsing once and
# checking the AST root -- rather than regex-stripping a leading prefix --
# means only a genuine *top-level* assignment is removed: a named function
# argument (`f(x = 1) + CEps`) has a top-level `+`, and a comparison
# (`C == CEps`) has a top-level `==`, so both are left untouched regardless
# of the LHS shape or surrounding whitespace.
.parse_observe_rhs <- function(rhs) {
  expr <- tryCatch(str2lang(rhs), error = function(e) NULL)
  if (is.null(expr)) return(NULL)
  if (is.call(expr) &&
      identical(expr[[1L]], as.name("=")) &&
      length(expr) == 3L) {
    return(expr[[3L]])
  }
  expr
}

# Replace every zero-argument call `name()` in `expr` with `replacement`,
# recursing into every sub-call. Used to turn `sigma()` -- which reports the
# fitted residual SD and never depends on the error variable -- into a
# plain symbol so it can be treated as an ordinary (constant-in-epsilon)
# free variable during differentiation.
.subst_call <- function(expr, name, replacement) {
  if (is.call(expr)) {
    if (identical(expr[[1L]], as.name(name)) && length(expr) == 1L) {
      return(replacement)
    }
    return(as.call(lapply(as.list(expr), .subst_call,
                          name = name, replacement = replacement)))
  }
  expr
}

# A handful of well-separated numeric probe points for the free variables
# (everything in the expression except the error variable itself). Using
# more than one point lets the classifier check whether a quantity actually
# stays constant as the surrounding model variables change, rather than
# relying on how the expression happens to be written.
.numeric_probes <- function(free_vars) {
  if (!length(free_vars)) return(list(list()))
  bases <- c(2, 7, 15)
  lapply(bases, function(b) {
    vals <- b + seq_along(free_vars) * 1.3
    stats::setNames(as.list(vals), free_vars)
  })
}

# Evaluate `expr` against a named list of numeric bindings, returning
# `NA_real_` on any error or non-scalar-numeric result instead of raising --
# callers treat evaluation failure as "can't tell" (role "unknown"), not as
# a hard error.
.safe_eval <- function(expr, env_list) {
  env <- list2env(env_list, parent = baseenv())
  tryCatch({
    v <- eval(expr, env)
    if (!is.numeric(v) || length(v) != 1L || is.na(v)) NA_real_ else as.numeric(v)
  }, error = function(e) NA_real_)
}

# Decide a role from three probe points' worth of noise-free prediction
# (`preds`), first derivative of the observe expression with respect to the
# error variable at that prediction (`g1`), and second derivative (`g2`).
#
# `g1` is the local error-magnitude function: how much one unit of the
# error variable perturbs the observation. `g1` constant across probes
# (independent of the prediction and every other model variable) means the
# noise is added on the raw scale -- additive error, for which
# `multiplicative_cv` is simply the wrong scale. `g1` proportional to the
# prediction, combined with a vanishing second derivative (the expression
# is exactly linear in the error variable), is the one shape
# `multiplicative_cv` (100 * sigma) is exact for. Everything else --
# `g1` depending on the prediction in a non-proportional way (combined,
# mix-ratio, power, ...), or a non-zero second derivative (the error enters
# non-linearly, e.g. multiplicatively on a log scale, where 100 * sigma is
# only a small-sigma approximation) -- is reported as "other" so the caller
# advises rather than silently assuming proportional error.
.classify_from_probes <- function(preds, g1, g2) {
  rel_tol <- function(x) 1e-6 * max(1, abs(x[1]))

  if (diff(range(g1)) < rel_tol(g1)) {
    if (abs(g1[1]) < 1e-9) return("unknown")
    return("additive")
  }

  if (any(abs(preds) < 1e-9)) return("other")
  ratio <- g1 / preds
  if (!all(is.finite(ratio)) || diff(range(ratio)) >= rel_tol(ratio)) {
    return("other")
  }

  if (all(abs(g2) < 1e-6 * max(1, abs(g1[1])))) "proportional" else "other"
}

# Classify the role of `sigma` inside a single observe expression by
# differentiating it symbolically (base R's `stats::D()`) rather than
# pattern-matching the source text. This works the same way regardless of
# how the expression happens to be written or which tool produced it --
# only the actual dependence of the expression on `sigma` matters.
#
# Expressions that `D()` cannot differentiate (a function outside its
# built-in derivative table, or anything that fails to parse as a plain R
# expression) fall back to "unknown"; callers treat "unknown" the same way
# as a detected non-proportional shape -- both mean the `multiplicative_cv`
# default should not be assumed correct.
.classify_sigma_role <- function(rhs, sigma) {
  expr <- .parse_observe_rhs(rhs)
  if (is.null(expr)) return("unknown")

  expr <- .subst_call(expr, "sigma", as.name(".sigma_ref"))

  g1 <- tryCatch(stats::D(expr, sigma), error = function(e) NULL)
  if (is.null(g1)) return("unknown")
  g2 <- tryCatch(stats::D(g1, sigma), error = function(e) NULL)
  if (is.null(g2)) return("unknown")

  # `.sigma_ref` must stay out of the probed free variables: it stands in
  # for `sigma()`, the model's own fitted residual SD, which is a single
  # fixed number at classification time -- not a quantity that varies
  # from one probe point to the next like the surrounding model
  # variables do. Pinning it to a constant keeps the probes isolating
  # only genuine free-variable dependence.
  free_vars <- setdiff(all.vars(expr), c(sigma, ".sigma_ref"))
  probes <- lapply(.numeric_probes(free_vars), function(env) {
    env[[sigma]] <- 0
    env[[".sigma_ref"]] <- 1
    list(pred = .safe_eval(expr, env),
        g1   = .safe_eval(g1, env),
        g2   = .safe_eval(g2, env))
  })

  vals <- unlist(probes, use.names = FALSE)
  if (anyNA(vals) || !all(is.finite(vals))) return("unknown")

  preds <- vapply(probes, `[[`, numeric(1), "pred")
  g1s   <- vapply(probes, `[[`, numeric(1), "g1")
  g2s   <- vapply(probes, `[[`, numeric(1), "g2")
  .classify_from_probes(preds, g1s, g2s)
}

# Walk every observe(...) block, extract its first identifier (the ObsName)
# and the unique sigma drawn from the closed `sigma_names` set. Every
# observe's right-hand side is expected to reference exactly one error
# variable, so `intersect` must yield length 1; anything else means the PML
# source can't be mapped and callers are told to fall back to
# engine-reported shrinkage instead. Also returns a per-sigma role
# classification used by the residual-transform advisories in
# `get_summaryNlme()`.
.map_obs_to_sigma <- function(xpdb_code, sigma_names) {
  src <- .strip_pml_comments(paste(xpdb_code, collapse = "\n"))
  blocks <- .match_keyword_blocks(src, "observe")
  if (!length(blocks)) {
    return(list(map = stats::setNames(character(), character()),
                roles = stats::setNames(character(), character())))
  }

  toks <- lapply(blocks, .identifier_tokens)
  obs_names <- vapply(toks, `[`, character(1), 1L)
  sig <- vapply(toks, function(t) {
    hit <- intersect(t[-1L], sigma_names)
    if (length(hit) == 1L) hit else NA_character_
  }, character(1))

  if (anyNA(sig)) {
    bad <- obs_names[is.na(sig)]
    stop("Cannot identify a unique sigma in observe(",
         paste(bad, collapse = ", "), "). ",
         "Use shrinkage = \"engine\".",
         call. = FALSE)
  }

  per_observe <- vapply(seq_along(blocks),
                        function(i) .classify_sigma_role(blocks[[i]], sig[[i]]),
                        character(1))
  roles <- vapply(unique(sig), function(s) {
    these <- per_observe[sig == s]
    if (length(unique(these)) == 1L) these[[1]] else "mixed"
  }, character(1))

  list(map = stats::setNames(sig, obs_names),
       roles = stats::setNames(unname(roles), unique(sig)))
}