R/foceiCovAnalytic.R

Defines functions .foceiAnalyticErrFull .omegaVarCovDeriv .foceiAnalyticSolveTol .hasRxSens .foceiCalcRanalytic .foceiAnalyticAssembleR .foceiAnalyticAssembleRFR .foceiAnalyticDirections .foceiEtaOccurrence .foceiNameOccurrence .foceiAnalyticJacLambda .foceiAnalyticDvSensLambda2 .foceiAnalyticDvSensLambda .foceiAnalyticTbsY .foceiCholSECovActive .foceiAnalyticFallback .foceiAnalyticErrWarn .foceiInstallAnalyticCov

# Analytic FOCEI observed-information covariance (covType="analytic"): the exact
# R-matrix from analytic 1st/2nd-order sensitivities (rxode2 .rxSens) with the
# 3rd-order tensor by Shi (2021) finite differences (keeps the augmented ODE at
# O(ndir^2), not O(ndir^3)).  Any solve/scope failure returns NULL and the caller
# falls back to the finite-difference Hessian.
#
# The sensitivity machinery runs over a "direction set": one direction per eta
# (ETA_i_) plus one direction per non-mu-referenced structural theta (THETA_j_,
# e.g. a covariate coefficient).  A mu-referenced theta reuses its eta's
# direction, so a fully mu-referenced model has ndir == neta.

#' Install the stashed analytic covariance as `fit$cov` (the native path fills
#' only the theta block).  No-op on the FD fallback (which foceiCalcR already
#' warned about).  `covFull=TRUE` (default) installs the full theta+sigma+Omega
#' matrix (identical theta SEs); `covFull=FALSE` installs the structural-theta
#' submatrix (NONMEM-matched theta cov, backwards-compatible shape) -- the assembly
#' is always full.
#' @param .ret focei fit environment
#' @noRd
.foceiInstallAnalyticCov <- function(.ret) {
  # only covMethod="r" installs the analytic R^-1; "r,s"/"s" keep the native
  # sandwich / S-matrix cov (which the analytic R already fed via covR).
  if (!identical(as.integer(rxode2::rxGetControl(.ret$ui, "covMethod", 2L)), 2L)) return(invisible())
  if (!exists(".analyticCov", envir = .ret, inherits = FALSE)) return(invisible())
  .cov <- get(".analyticCov", envir = .ret)
  if (!is.matrix(.cov) || !all(is.finite(.cov))) return(invisible())
  .full <- isTRUE(rxode2::rxGetControl(.ret$ui, "covFull", TRUE))
  if (!.full && exists(".analyticThetaNames", envir = .ret, inherits = FALSE)) {
    .th <- get(".analyticThetaNames", envir = .ret)          # structural cov-theta block only
    .th <- .th[.th %in% rownames(.cov)]
    if (length(.th) > 0L) .cov <- .cov[.th, .th, drop = FALSE]
  }
  # PD guard: an indefinite (near-boundary) inverse installs negative variances ->
  # NaN SEs.  Reject and keep the native/FD cov rather than a plausible-looking wrong one.
  .ev <- suppressWarnings(eigen(.cov, symmetric = TRUE, only.values = TRUE)$values)
  if (any(diag(.cov) <= 0) || !all(is.finite(.ev)) || min(.ev) <= 0) {
    warning("analytic covariance is not positive definite; keeping the finite-difference covariance",
            call. = FALSE)
    return(invisible())
  }
  .ret$cov <- .cov                       # analytic-tier cov already carries dimnames
  .ret$covMethod <- "analytic"           # report the analytic observed information (not "r")
  # covFull=TRUE swaps in a larger matrix than C++ foceiFinalizeTables saw, so its
  # condition numbers (computed from the theta-only native cov) are stale -- recompute.
  if (.full) .foceiCovCondition(.ret, .cov, .ev)
  invisible()
}

#' Emit the reason the analytic (covType="analytic") observed-information R-matrix
#' is unavailable for this model and return `NULL` so the caller drops to the
#' finite-difference Hessian.  Only reached in opted-in paths -- the native
#' `.foceiCalcRanalytic` hook checks `covType=="analytic"` first and
#' `foceiCovAnalytic` is an explicit call -- so nothing prints for the default
#' `covType="fd"`.  Kept as a plain `message()` (informational, not a warning):
#' falling back to FD is a valid result, not an error.
#' @param reason short human-readable phrase naming the out-of-scope feature
#' @return `NULL`, invisibly usable as `return(.foceiAnalyticFallback(...))`
#' @noRd
.foceiAnalyticErrWarn <- function(site) {
  # An error here is otherwise indistinguishable from an out-of-scope decline.  Surface
  # it as a warning so it reaches $runInfo (CLAUDE.md: <75 chars, no method prefix), and
  # NUMBER the site so the throwing entry point is identifiable from the fit alone:
  #   err1 = live-fit .foceiCalcRanalytic   err2 = standalone foceiCovAnalytic
  #   err3 = per-subject AGQ node assembly
  function(e) {
    .m <- conditionMessage(e)
    if (nchar(.m) > 54L) .m <- paste0(substr(.m, 1L, 51L), "...")
    warning(sprintf("analytic err%d: %s", site, .m), call. = FALSE)
    NULL
  }
}

.foceiAnalyticFallback <- function(reason) {
  message("covType=\"analytic\": ", reason,
          " is out of analytic-covariance scope; using the finite-difference covariance instead")
  NULL
}

# Shared cholSECov analytic-cov guard, used by both .foceiCalcRanalytic (live hook) and
# .foceiCovAnalyticCalc (standalone entry) so the gate and its reason cannot diverge.
# foceiCalcCov() re-sets every subject's doChol from cholSECov before the cov step (see
# src/inner.cpp, `fInd->doChol = !op_focei.cholSECov`), so the objective's log|Ht| then
# comes from the generalized Cholesky, not chol().  For a non-PD Ht that is a different
# quantity than the analytic R differentiates.  (cholSEOpt is the OPTIMIZATION-phase flag
# and does not affect the cov step.)
.foceiCholSECovReason <- "cholSECov=TRUE (the covariance step re-factors the eta Hessian)"
.foceiCholSECovActive <- function(ui) isTRUE(rxode2::rxGetControl(ui, "cholSECov", FALSE))

#' Transform the observed DV onto the rx_pred_ (transformed) scale for a both-sides
#' transform, using the solved per-observation transform parameters `trans`
#' (lambda/yj/low/hi, each length nobs).  Runs the SAME C++ transform the inner problem
#' uses (`_nlmixr2est_powerD` == `tbs()`), which requires all arguments the same length --
#' the per-observation vectors already are, so a multi-endpoint mix of transforms is
#' handled directly.  Returns `dv` unchanged when `trans` is NULL.
#' @noRd
.foceiAnalyticTbsY <- function(dv, trans) {
  if (is.null(trans)) return(dv)
  .Call(`_nlmixr2est_powerD`, as.double(dv), as.double(trans$lambda),
        as.integer(trans$yj), as.double(trans$low), as.double(trans$hi))
}

#' Per-observation transform value lambda-derivative dy'/dlambda of the DV
#' (`d tbs(DV,lambda)/dlambda`), used to build the residual DV-transform chain
#' when lambda is estimated.  Same per-obs `trans` contract as `.foceiAnalyticTbsY`.
#' @noRd
.foceiAnalyticDvSensLambda <- function(dv, trans) {
  if (is.null(trans)) return(rep(0, length(dv)))
  .Call(`_nlmixr2est_powerDLambda`, as.double(dv), as.double(trans$lambda),
        as.integer(trans$yj), as.double(trans$low), as.double(trans$hi))
}

#' Per-observation d2 tbs(DV,lambda)/dlambda2 (DV second lambda-derivative).
#' @noRd
.foceiAnalyticDvSensLambda2 <- function(dv, trans) {
  if (is.null(trans)) return(rep(0, length(dv)))
  .Call(`_nlmixr2est_powerDLambda2`, as.double(dv), as.double(trans$lambda),
        as.integer(trans$yj), as.double(trans$low), as.double(trans$hi))
}

#' Per-observation d log|dy'/dDV| / dlambda -- the transform Jacobian's
#' lambda-derivative, the extra term the OFV gradient's lambda column carries
#' (-2 * sum over obs); it cancels in the observed-information covariance.
#' @noRd
.foceiAnalyticJacLambda <- function(dv, trans) {
  if (is.null(trans)) return(rep(0, length(dv)))
  .Call(`_nlmixr2est_powerDL`, as.double(dv), as.double(trans$lambda),
        as.integer(trans$yj), as.double(trans$low), as.double(trans$hi))
}

#' Count literal occurrences of each of `.names` on the model-expression RHS.
#'
#' Shared by the eta-occurrence guard (a random effect reused across parameters breaks
#' the covariate/theta direction reuse) and the structural-theta guard (a structural
#' parameter appearing in more than one mu-expression breaks the eta-less covariate
#' reuse -- `df/dtheta` would then gain a path `df/db` lacks).  Returns a named integer
#' vector; `2L` (i.e. ">1", conservative -> no reuse) when `lstExpr` is unavailable.
#' @noRd
.foceiNameOccurrence <- function(ui, .names) {
  .lst <- tryCatch(ui$lstExpr, error = function(e) NULL)
  .count <- function(.e, .nm) {
    if (is.name(.e)) return(as.integer(identical(as.character(.e), .nm)))
    if (is.call(.e)) return(sum(vapply(as.list(.e)[-1L], .count, integer(1), .nm)))
    0L
  }
  .isAssign <- function(.ex) is.call(.ex) && length(.ex) == 3L &&
    (identical(.ex[[1]], as.name("<-")) || identical(.ex[[1]], as.name("=")))
  stats::setNames(vapply(as.character(.names), function(.n)
    if (is.null(.lst)) 2L else sum(vapply(.lst, function(.ex)
      if (.isAssign(.ex)) .count(.ex[[3]], .n) else 0L, integer(1))), integer(1)),
    as.character(.names))
}
#' Literal occurrence count of each diagonal eta across the model RHS.
#'
#' The sensitivity direction-reuse identities -- a mu-ref theta reusing its eta's
#' direction (`df/dtheta = df/deta`) and a mu-ref covariate coefficient scaling it
#' (`df/db = cov*df/deta`) -- are exact only when the eta enters the model through a
#' single additive position.  An eta shared across parameters (e.g. `eta.cl` in both
#' `cl` and `v`) makes `df/dtheta` (one parameter) differ from `df/deta` (all of
#' them), so the reuse would be wrong.  rxode2's mu-reference frames still map such an
#' eta, so its multiplicity is counted here.  Named by the diagonal-eta names in
#' `neta1` order; `2L` (conservative ">1") when `lstExpr` is unavailable.
#' @noRd
.foceiEtaOccurrence <- function(ui) {
  .idf <- ui$iniDf
  .etaRows <- .idf[!is.na(.idf$neta1) & .idf$neta1 == .idf$neta2, , drop = FALSE]
  .etaRows <- .etaRows[order(.etaRows$neta1), , drop = FALSE]
  .foceiNameOccurrence(ui, as.character(.etaRows$name))
}

#' Scope-relevant structural thetas + sensitivity direction map: one direction per
#' eta (ETA_i_) plus one per non-mu-ref structural theta (THETA_j_); a mu-ref theta
#' reuses its eta's direction (so a fully mu-ref model has ndir == neta).  Sigma and
#' IOV SD thetas (`sgName`/`iovVars`) are not structural.  `sharedEta` is a logical
#' flag per diagonal eta (in `neta1` order) marking an eta reused across parameters,
#' for which the theta reuse is invalid.  `NULL` if none remain.
#' @noRd
.foceiAnalyticDirections <- function(ini, thetaForEta, sgName, neta, iovVars = character(0),
                                     sharedEta = logical(0)) {
  thRows <- ini[!is.na(ini$ntheta), , drop = FALSE]
  thRows <- thRows[order(thRows$ntheta), , drop = FALSE]
  thStructRows <- thRows[!.iniIsFixed(ini, thRows$name) & !(thRows$name %in% sgName) &
                           !(thRows$name %in% iovVars), , drop = FALSE]
  thStruct <- thStructRows$name
  nth <- length(thStruct)
  if (nth == 0L) return(.foceiAnalyticFallback("no estimated structural parameter"))
  etaDirs <- paste0("ETA_", seq_len(neta), "_")
  dirTh <- integer(nth); nonMuTheta <- character(0)
  for (p in seq_len(nth)) {
    .mt <- which(thetaForEta == thStruct[p])        # eta(s) this theta mu-references
    if (length(.mt) > 1L)                           # shared by >1 eta -> summing routes NYI, fall back to FD
      return(.foceiAnalyticFallback("a theta shared by two random effects"))
    k <- if (length(.mt) == 1L) .mt else NA_integer_
    # A mu-ref theta reuses its eta's state-sensitivity direction for free (df/dtheta = df/deta,
    # since theta and eta enter the mu identically) -- but ONLY when the eta enters the model in a
    # single additive position.  A shared eta (e.g. eta.cl in both cl and v) makes df/dtheta (this
    # parameter only) differ from df/deta (all parameters the eta appears in), so the reuse is
    # invalid; treat the theta as non-mu and give it its OWN true-sensitivity direction (the eta
    # keeps its own, correct, direction).  Stays analytic -- just one extra direction for that
    # theta, only in the rare shared-eta case.
    .reuse <- !is.na(k) && !(k <= length(sharedEta) && isTRUE(sharedEta[k]))
    if (.reuse) {
      dirTh[p] <- k                                 # reuse eta k's direction (free)
    } else {
      nonMuTheta <- c(nonMuTheta, paste0("THETA_", thStructRows$ntheta[p], "_"))
      dirTh[p] <- neta + length(nonMuTheta)         # its own true-sensitivity direction
    }
  }
  dirs <- c(etaDirs, nonMuTheta)
  # Estimated boxCox/yeoJohnson lambda is an error parameter that (unlike an ordinary
  # sigma with df/dsigma=0) has a NONZERO prediction sensitivity df'/dlambda=rxTBSdL(f,
  # lambda), so it enters as its OWN theta-like direction (appended to the theta block):
  # the kernels then handle its pred + variance sensitivities exactly like a theta, and
  # the DV-transform residual chain (dy'/dlambda) + the -2 log|J| Jacobian are applied on
  # top (in the driver / kernel).  Fixed lambda is filtered out here (static transform).
  .lamRows <- ini[ini$name %in% sgName & !is.na(ini$err) &
                    ini$err %in% c("boxCox", "yeoJohnson") & !.iniIsFixed(ini, ini$name), , drop = FALSE]
  lamNames <- character(0); lamDir <- integer(0)
  if (nrow(.lamRows) > 0L) {
    .lamRows <- .lamRows[order(.lamRows$ntheta), , drop = FALSE]
    for (r in seq_len(nrow(.lamRows))) {
      dirs <- c(dirs, paste0("THETA_", .lamRows$ntheta[r], "_"))
      lamNames <- c(lamNames, .lamRows$name[r]); lamDir <- c(lamDir, length(dirs))
    }
    thStruct <- c(thStruct, lamNames); dirTh <- c(dirTh, lamDir); nth <- nth + length(lamNames)
  }
  lamIdx <- if (length(lamNames)) seq(nth - length(lamNames) + 1L, nth) else integer(0)  # 1-based theta-block positions
  # (f,R) covariance: fold the residual-error parameters (sigmas) into the DIRECTION set.
  # The NONMEM sigma (the unit N(0,1) residual variate) is FIXED to 1, so the estimated
  # error parameters enter ONLY through the weight/variance rx_r_ and NOT through the
  # prediction -- i.e. a sigma direction has df/dsigma=0 (a=A=Ath=0) and only carries the
  # variance sensitivities dR/dsigma (aR/AR/AthR).  Treating each sigma as its own direction
  # therefore lets the SAME 1st/2nd-order solve + Shi-FD 3rd-order machinery supply every
  # sigma derivative (including the 3rd-order mixed sigma-direction terms the determinant
  # block needs) with no special-casing.  NOTE: if the FOCEI methods are ever generalized so
  # the residual sigma itself is estimated (not fixed to 1), that estimated sigma would enter
  # the variance multiplicatively (R -> R*sigma^2) and would need its own chain here.
  .sgRows <- thRows[thRows$name %in% sgName & !(thRows$name %in% lamNames), , drop = FALSE]
  dirSg <- if (nrow(.sgRows) > 0L) seq_len(nrow(.sgRows)) + length(dirs) else integer(0)
  dirsCov <- c(dirs, if (nrow(.sgRows) > 0L) paste0("THETA_", .sgRows$ntheta, "_") else character(0))
  list(thStruct = thStruct, thStructRows = thStructRows, dirs = dirs,
       dirTh = dirTh, ndir = length(dirs), nth = nth,
       dirsCov = dirsCov, ndirCov = length(dirsCov), sgName = .sgRows$name,
       lamNames = lamNames, lamDir = lamDir, lamIdx = lamIdx,
       dirP = c(dirTh, dirSg))                          # every non-Omega param -> a direction
}

#' Sum the per-subject general (f,R) observed-information R over subjects (FOCEI).
#' Builds the augmented model over the covariance direction set (etas + non-mu thetas +
#' SIGMA thetas), Shi-FDs the 3rd-order tensors (Ath + AthR), and assembles via
#' `.foceiAnalyticSubjectRFR`.  Handles any variance structure; `NULL` on failure.
#' @noRd
.foceiAnalyticAssembleRFR <- function(ui, th, ebes, ids, data, Om, ef, neta, ndirP, dirP, omd,
                                      dirsCov, ndirCov, startedEnv = NULL, solveTol = 1e-10,
                                      interaction = 1L, foceType = 0L, lamDir = integer(0)) {
  # Route A (default ON; set FOCEI_NO_RSIG=1 to opt out): build the gradient's `dirs` model
  # (drop the sigma directions) so the (f,R) cov SHARES the gradient's already-compiled model
  # -- one augmented-model compile instead of two, and the sigma directions' wasted (all-zero)
  # state-sensitivity columns vanish.  The sigma tensor slots (aR/AR/AthR) are rebuilt from the
  # model's own rsig outputs in SolveAllFD3 (.foceiAnalyticExpandSigma), which expands the solved
  # E back to ndirCov -- numerically identical to carrying explicit sigma directions.
  .rsigA <- !nzchar(Sys.getenv("FOCEI_NO_RSIG"))
  .erN <- ui$iniDf$ntheta[!is.na(ui$iniDf$err) & !(ui$iniDf$err %in% c("boxCox", "yeoJohnson"))]
  .sigDirs <- if (.rsigA) intersect(dirsCov, paste0("THETA_", .erN, "_")) else character(0)
  .dirsF <- dirsCov[!(dirsCov %in% .sigDirs)]
  # Flag BEFORE building the augmented model: .foceiAnalyticAugModelDirs loads a new
  # compiled model, which calls rxSolveFree() on the fit's global solve.  If the build
  # then declines (is.null(am)/ndir mismatch/no rx_r_), the C++ hook must still restore the
  # freed solve before the FD fallback -- otherwise the FD sandwich solves against freed memory.
  if (!is.null(startedEnv)) assign(".analyticStarted", TRUE, startedEnv)
  am <- .foceiAnalyticAugModelDirs(ui, if (.rsigA) .dirsF else dirsCov)
  # Which rsig columns correspond to the directions Route A dropped?  am$sigTh is the
  # model's sigma-theta list (sorted) and fixes the rx_rsig_ column order; .sigDirs is
  # what we removed.  Any dropped direction the model does not emit rsig for is an
  # inconsistency, not something to paper over -- decline instead of mis-shaping E.
  .sigSel <- NULL
  if (!is.null(am) && length(.sigDirs) > 0L) {
    .sigN <- suppressWarnings(as.integer(sub("^THETA_(\\d+)_$", "\\1", .sigDirs)))
    .sigSel <- match(.sigN, am$sigTh)
    if (anyNA(.sigSel))
      return(.foceiAnalyticFallback("a sigma direction the model omits from rsig"))
  }
  if (is.null(am))
    return(.foceiAnalyticFallback("an augmented model that would not build"))
  if (am$ndir != (ndirCov - length(.sigDirs)))
    return(.foceiAnalyticFallback(paste0("an augmented-model direction mismatch (",
                                         am$ndir, "/", ndirCov - length(.sigDirs), ")")))
  if (!isTRUE(am$hasRvar))
    return(.foceiAnalyticFallback("an augmented model without rx_r_"))
  np <- ndirP + omd$nom; Oi <- solve(Om)
  etav <- paste0("ETA_", seq_len(neta), "_")
  .foce <- identical(as.integer(interaction), 0L)      # FOCE re-solves EBEs to S_FOCE=0
  .fp <- identical(as.integer(foceType), 1L)           # foce+ keeps the live conditional R
  .byId <- split(data, as.character(data$ID))
  .idCode <- if (is.factor(ids)) as.integer(ids) else match(ids, sort(unique(ids)))
  R <- matrix(0, np, np)
  if (.foce) {
    # FOCE cov: the per-subject eta=0 population solve + EBE re-solve + 3rd-order Shi FD3 stay
    # in R (inherently per-subject), then ONE OpenMP C++ call (foceiRAllFoceFR_) sums the
    # observed information over subjects.  The frozen-R0 sensitivities are resolved per subject
    # (nonmem: aRe/ARe=0, aRc/ARc/R0 from E0; foce+: all from the eta-hat solve E).
    nsub <- length(ids); Elist <- vector("list", nsub); E0list <- vector("list", nsub)
    eta0list <- vector("list", nsub); nobsAll <- integer(nsub)
    # BATCHED (f,R) FOCE/foce+ solves: the eta=0 population solve (nonmem frozen R0) is batched,
    # the EBE re-solve stays per-subject (Newton), then ONE batched SolveAllFD3 delivers f/a/A/Ath
    # AND R/aR/AR for all subjects (withR=FALSE: FOCE reads R/aR/AR from the base solve + Ath, but
    # never AthR).  foce+ (foceType=1) keeps the live R (no eta=0 solve).  Per-subject Shi fallback.
    .obsAll <- lapply(seq_len(nsub), function(i) { .s <- .byId[[as.character(.idCode[i])]]
      if (is.null(.s) || nrow(.s) == 0L) NULL else .s[.s$EVID == 0, , drop = FALSE] })
    if (any(vapply(.obsAll, is.null, logical(1L))))
      return(.foceiAnalyticFallback("a subject with no observations"))
    .obsT <- lapply(.obsAll, function(.o) .o$TIME)
    E0all <- if (!.fp) .foceiAnalyticSolveAll(am, th, matrix(0, nsub, neta), .idCode, data, .obsT, solveTol) else NULL
    if (!.fp && is.null(E0all))
      return(.foceiAnalyticFallback("a failed eta=0 population solve"))
    if (.rsigA && !.fp) {                                # Route A: rebuild frozen-R0 sigma slots (aR/AR) from rsig
      .ns0 <- ncol(E0all[[1L]]$Rsig)
      if (!is.null(.ns0) && .ns0 > 0L) E0all <- lapply(E0all, function(.E0) .foceiAnalyticExpandSigma(.E0, .ns0, neta, NULL, NULL, .sigSel))
    }
    eta0Mat <- .foceiAnalyticFoceEbeBatch(am, th, ebes, .idCode, data, .obsAll, .obsT, etav, Oi, neta,
                                          solveTol, foceType = foceType, E0all = E0all)   # batched EBE re-solve
    if (is.null(eta0Mat))
      return(.foceiAnalyticFallback("an EBE re-solve that will not solve"))
    .batch <- !nzchar(Sys.getenv("FOCEI_NO_FD3_BATCH"))
    .EsAll <- if (.batch) .foceiAnalyticSolveAllFD3(am, th, eta0Mat, .idCode, data, .obsT, tol = solveTol, withR = FALSE,
                                                   sigSel = .sigSel) else NULL
    if (.batch && is.null(.EsAll)) .batch <- FALSE
    for (i in seq_len(nsub)) {
      s <- .byId[[as.character(.idCode[i])]]; obs <- .obsAll[[i]]; eta0 <- eta0Mat[i, ]
      E0 <- if (.fp) NULL else E0all[[i]]
      E <- if (.batch) .EsAll[[i]] else .foceiAnalyticSolveSubjectFD3(am, c(th, setNames(eta0, etav)), s, obs$TIME, tol = solveTol, withR = FALSE)
      if (is.null(E))
        return(.foceiAnalyticFallback("a failed subject sensitivity solve"))
      if (isTRUE(ef$canVanish)) { .fa <- abs(E$f)
        if (any(!is.finite(.fa)) || min(.fa) < 1e-6 * max(.fa))
          return(.foceiAnalyticFallback("proportional error near a zero prediction")) }
      E$y <- .foceiAnalyticTbsY(obs$DV, E$trans)
      # E0list[i] <- list(E0) keeps a NULL slot (foce+) without shrinking the list.
      Elist[[i]] <- E; E0list[i] <- list(E0); eta0list[[i]] <- eta0; nobsAll[i] <- length(E$f)
    }
    .fpG <- identical(as.integer(foceType), 1L) || is.null(E0list[[1L]])
    totObs <- sum(nobsAll); off <- c(0L, cumsum(nobsAll)); nd2 <- ndirCov * ndirCov
    aB <- matrix(0, totObs, ndirCov); aReB <- matrix(0, totObs, ndirCov); aRcB <- matrix(0, totObs, ndirCov)
    AB <- array(0, c(totObs, ndirCov, ndirCov)); AReB <- array(0, c(totObs, ndirCov, ndirCov)); ARcB <- array(0, c(totObs, ndirCov, ndirCov))
    AthB <- array(0, c(totObs, neta, nd2))
    fB <- numeric(totObs); yB <- numeric(totObs); R0B <- numeric(totObs); ehatB <- matrix(0, nsub, neta)
    dvSensB <- if (length(lamDir)) matrix(0, totObs, ndirCov) else matrix(0, totObs, 0L)
    dvSens2B <- dvSensB
    .hasCens <- (!is.null(data$CENS) && any(data$CENS != 0, na.rm = TRUE)) ||
      (!is.null(data$LIMIT) && any(is.finite(data$LIMIT)))
    censB <- if (.hasCens) integer(totObs) else integer(0)
    limB <- if (.hasCens) rep(NA_real_, totObs) else numeric(0)
    for (i in seq_len(nsub)) {
      E <- Elist[[i]]; E0 <- E0list[[i]]; no <- nobsAll[i]; rows <- (off[i] + 1L):off[i + 1L]
      aB[rows, ] <- E$a; AB[rows, , ] <- E$A; AthB[rows, , ] <- array(E$Ath, c(no, neta, nd2))
      fB[rows] <- E$f; yB[rows] <- E$y
      if (.fpG) { R0B[rows] <- E$R; aReB[rows, ] <- E$aR; aRcB[rows, ] <- E$aR; AReB[rows, , ] <- E$AR; ARcB[rows, , ] <- E$AR }
      else { R0B[rows] <- E0$R; aRcB[rows, ] <- E0$aR; ARcB[rows, , ] <- E0$AR }   # aReB/AReB stay 0 (frozen)
      if (length(lamDir) || .hasCens) {
        s <- .byId[[as.character(.idCode[i])]]; obs <- s[s$EVID == 0, , drop = FALSE]
        if (length(lamDir)) {                          # DV-transform chain (estimated lambda)
          dvSensB[rows, lamDir] <- .foceiAnalyticDvSensLambda(obs$DV, E$trans)
          dvSens2B[rows, lamDir] <- .foceiAnalyticDvSensLambda2(obs$DV, E$trans)
        }
        if (.hasCens) {
          censB[rows] <- if (is.null(obs$CENS)) 0L else as.integer(obs$CENS)
          .lim <- if (is.null(obs$LIMIT)) rep(NA_real_, length(rows)) else as.numeric(obs$LIMIT)
          limB[rows] <- .foceiAnalyticTbsY(.lim, E$trans)
        }
      }
      ehatB[i, ] <- eta0list[[i]]
    }
    nom <- omd$nom
    dOiC <- array(0, c(neta, neta, max(nom, 1L)))
    if (nom > 0L) for (k in seq_len(nom)) dOiC[, , k] <- omd$dOi[[k]]
    d2OiC <- array(0, c(neta, neta, max(nom * nom, 1L)))
    if (nom > 0L) for (aa in seq_len(nom)) for (bb in seq_len(nom)) d2OiC[, , (aa - 1L) * nom + bb] <- omd$d2Oi[[aa]][[bb]]
    d2LD <- if (nom > 0L) omd$d2LD else matrix(0, 1, 1)
    ncores <- tryCatch(as.integer(rxode2::getRxThreads()), error = function(e) 1L)
    if (length(ncores) != 1L || is.na(ncores) || ncores < 1L) ncores <- 1L
    R <- tryCatch(foceiRAllFoceFR_(aB, AB, AthB, aReB, aRcB, AReB, ARcB, dvSensB, dvSens2B,
                                   as.integer(censB), as.numeric(limB), fB, yB, R0B, ehatB, as.integer(off),
                                   Oi, dOiC, d2OiC, d2LD, neta, ndirCov, ndirP, nom, as.integer(dirP), ncores),
                  error = function(e) NULL)
    if (is.null(R))
      return(.foceiAnalyticFallback("a failed observed-information kernel"))
    if (!all(is.finite(R)))
      return(.foceiAnalyticFallback("a non-finite observed information"))
    return(R)
  }
  # FOCEI: per-subject FD3 (3rd-order Shi) solves collected in R, then ONE OpenMP C++ call
  # (foceiRAllFR_) sums the observed information over subjects -- no per-subject R<->C++ round-trip.
  nsub <- length(ids); Elist <- vector("list", nsub); nobsAll <- integer(nsub)
  # BATCHED (f,R) FOCEI solves: ONE SolveAllFD3 gives f/a/A/Ath AND R/aR/AR/AthR for all subjects
  # (withR=TRUE), instead of the per-subject Shi.  Per-subject Shi is the fallback.
  .obsAll <- lapply(seq_len(nsub), function(i) { .s <- .byId[[as.character(.idCode[i])]]
    if (is.null(.s) || nrow(.s) == 0L) NULL else .s[.s$EVID == 0, , drop = FALSE] })
  if (any(vapply(.obsAll, is.null, logical(1L))))
    return(.foceiAnalyticFallback("a subject with no observations"))
  .obsT <- lapply(.obsAll, function(.o) .o$TIME)
  .batch <- !nzchar(Sys.getenv("FOCEI_NO_FD3_BATCH"))
  .EsAll <- if (.batch) .foceiAnalyticSolveAllFD3(am, th, ebes, .idCode, data, .obsT, tol = solveTol, withR = TRUE,
                                                 sigSel = .sigSel) else NULL
  if (.batch && is.null(.EsAll)) .batch <- FALSE
  for (i in seq_len(nsub)) {
    s <- .byId[[as.character(.idCode[i])]]; obs <- .obsAll[[i]]
    E <- if (.batch) .EsAll[[i]] else .foceiAnalyticSolveSubjectFD3(am, c(th, setNames(ebes[i, ], etav)), s, obs$TIME, tol = solveTol, withR = TRUE)
    if (is.null(E))
      return(.foceiAnalyticFallback("a failed subject sensitivity solve"))
    if (isTRUE(ef$canVanish)) { .fa <- abs(E$f)
      if (any(!is.finite(.fa)) || min(.fa) < 1e-6 * max(.fa))
        return(.foceiAnalyticFallback("proportional error near a zero prediction")) }
    E$y <- .foceiAnalyticTbsY(obs$DV, E$trans)
    Elist[[i]] <- E; nobsAll[i] <- length(E$f)
  }
  totObs <- sum(nobsAll); off <- c(0L, cumsum(nobsAll)); nd2 <- ndirCov * ndirCov
  aB <- matrix(0, totObs, ndirCov); aRB <- matrix(0, totObs, ndirCov)
  AB <- array(0, c(totObs, ndirCov, ndirCov)); ARB <- array(0, c(totObs, ndirCov, ndirCov))
  AthB <- array(0, c(totObs, neta, nd2)); AthRB <- array(0, c(totObs, neta, nd2))
  fB <- numeric(totObs); yB <- numeric(totObs); RB <- numeric(totObs); ehatB <- matrix(0, nsub, neta)
  dvSensB <- if (length(lamDir)) matrix(0, totObs, ndirCov) else matrix(0, totObs, 0L)
  dvSens2B <- dvSensB
  # censored (M2/M3/M4): per-obs CENS + transformed LIMIT drive the censored score partials
  # (the determinant stays Gauss-Newton, matching censOption="gauss"); empty when no censoring.
  .hasCens <- (!is.null(data$CENS) && any(data$CENS != 0, na.rm = TRUE)) ||
    (!is.null(data$LIMIT) && any(is.finite(data$LIMIT)))
  censB <- if (.hasCens) integer(totObs) else integer(0)
  limB <- if (.hasCens) rep(NA_real_, totObs) else numeric(0)
  for (i in seq_len(nsub)) {
    E <- Elist[[i]]; no <- nobsAll[i]; rows <- (off[i] + 1L):off[i + 1L]
    aB[rows, ] <- E$a; aRB[rows, ] <- E$aR; AB[rows, , ] <- E$A; ARB[rows, , ] <- E$AR
    AthB[rows, , ] <- array(E$Ath, c(no, neta, nd2)); AthRB[rows, , ] <- array(E$AthR, c(no, neta, nd2))
    fB[rows] <- E$f; yB[rows] <- E$y; RB[rows] <- E$R; ehatB[i, ] <- ebes[i, ]
    if (length(lamDir) || .hasCens) {
      s <- .byId[[as.character(.idCode[i])]]; obs <- s[s$EVID == 0, , drop = FALSE]
      if (length(lamDir)) {                            # DV-transform chain (estimated lambda)
        dvSensB[rows, lamDir] <- .foceiAnalyticDvSensLambda(obs$DV, E$trans)
        dvSens2B[rows, lamDir] <- .foceiAnalyticDvSensLambda2(obs$DV, E$trans)
      }
      if (.hasCens) {
        censB[rows] <- if (is.null(obs$CENS)) 0L else as.integer(obs$CENS)
        .lim <- if (is.null(obs$LIMIT)) rep(NA_real_, length(rows)) else as.numeric(obs$LIMIT)
        limB[rows] <- .foceiAnalyticTbsY(.lim, E$trans)   # transform the censoring bound like the DV
      }
    }
  }
  nom <- omd$nom
  dOiC <- array(0, c(neta, neta, max(nom, 1L)))
  if (nom > 0L) for (k in seq_len(nom)) dOiC[, , k] <- omd$dOi[[k]]
  d2OiC <- array(0, c(neta, neta, max(nom * nom, 1L)))
  if (nom > 0L) for (aa in seq_len(nom)) for (bb in seq_len(nom)) d2OiC[, , (aa - 1L) * nom + bb] <- omd$d2Oi[[aa]][[bb]]
  d2LD <- if (nom > 0L) omd$d2LD else matrix(0, 1, 1)
  ncores <- tryCatch(as.integer(rxode2::getRxThreads()), error = function(e) 1L)
  if (length(ncores) != 1L || is.na(ncores) || ncores < 1L) ncores <- 1L
  R <- tryCatch(foceiRAllFR_(aB, AB, AthB, aRB, ARB, AthRB, dvSensB, dvSens2B, as.integer(censB), as.numeric(limB),
                             fB, yB, RB, ehatB, as.integer(off),
                             Oi, dOiC, d2OiC, d2LD, neta, ndirCov, ndirP, nom, as.integer(dirP), ncores),
                error = function(e) NULL)
  if (is.null(R))
    return(.foceiAnalyticFallback("a failed observed-information kernel"))
  if (!all(is.finite(R)))
    return(.foceiAnalyticFallback("a non-finite observed information"))
  R
}

#' Sum the per-subject analytic observed-information R (theta+sigma+Omega) over
#' subjects; `NULL` on any build/solve/non-finite failure.  `startedEnv` (when
#' given) is flagged `.analyticStarted` before the first solve so the C++ hook
#' skips its FD fallback rather than run it on the replaced global solve.
#' @noRd
.foceiAnalyticAssembleR <- function(ui, th, ebes, ids, data, Om, ef, neta, nth, nsg, omd,
                                    dirs, dirTh, ndir, startedEnv = NULL, solveTol = 1e-10,
                                    iovDirScale = NULL, etaScale = NULL, interaction = 1L,
                                    foceType = 0L) {
  # Flag BEFORE building the augmented model: .foceiAnalyticAugModelDirs loads a new
  # compiled model, which calls rxSolveFree() on the fit's global solve.  If the build
  # then declines (is.null(am)/ndir mismatch), the C++ hook must still restore the freed
  # solve before the FD fallback -- otherwise the FD sandwich solves against freed memory.
  if (!is.null(startedEnv)) assign(".analyticStarted", TRUE, startedEnv)
  am <- .foceiAnalyticAugModelDirs(ui, dirs)
  if (is.null(am) || am$ndir != ndir) return(NULL)
  # AGQ (nAGQ>1) swaps the objective's data half for -log(sum_k a_k); NULL keeps the
  # FOCEI/Laplace path byte-for-byte.
  .nAGQ <- as.integer(rxode2::rxGetControl(ui, "nAGQ", 1L))
  .ag <- if (.nAGQ > 1L) .agq(neta, .nAGQ) else NULL
  np <- nth + nsg + omd$nom
  etav <- paste0("ETA_", seq_len(neta), "_")
  .foce <- identical(as.integer(interaction), 0L)     # FOCE re-solves EBEs to S_FOCE=0
  # IOV reparameterization to xi = w * eta (unit occasion eta -> variance w^2): the
  # augmented model is solved at the ACTUAL (Param A) EBEs, then the occasion-eta
  # sensitivities are rescaled by 1/w (iovDirScale) and their EBEs by w (etaScale)
  # so the Omega-variance branch sees a genuine variance-w^2 random effect.
  rescale <- !is.null(iovDirScale) && any(iovDirScale != 1)
  R <- matrix(0, np, np)
  Oi <- solve(Om)                                     # constant across subjects (invert once)
  # `data` (dataSav) keys on the etTrans integer code (1..N); `ids` (etaObf$ID /
  # fit$eta$ID) is a factor whose LABELS are the original dataset IDs.  Key BOTH on
  # the integer code so non-1..N IDs (e.g. 101,102,..) and permutations join right.
  .byId <- split(data, as.character(data$ID))         # pre-split once (avoid per-subject rescan)
  .idCode <- if (is.factor(ids)) as.integer(ids) else match(ids, sort(unique(ids)))
  # Batch the 3rd-order solve across ALL subjects (FOCEI *and* FOCE, no IOV rescale) so both take
  # the IDENTICAL method and both get the speedup (1 + 2*neta population solves vs the per-subject
  # Shi's O(nsub*neta)).  CORRECTED FOCE freezes R0 at the eta=0 population solve (batched) and
  # re-solves each EBE (per-subject Newton, censoring-aware) to S_FOCE=0; that eta0 matrix feeds
  # the batched FD3.  foce+ keeps the live R (no eta=0 solve).  Per-subject Shi is the fallback.
  .obsAll <- lapply(seq_along(ids), function(i) { .s <- .byId[[as.character(.idCode[i])]]
    if (is.null(.s) || nrow(.s) == 0L) NULL else .s[.s$EVID == 0, , drop = FALSE] })
  if (any(vapply(.obsAll, is.null, logical(1L)))) return(NULL)   # unmatched subject -> caller FD
  .obsT <- lapply(.obsAll, function(.o) .o$TIME)
  .needF0 <- .foce && identical(as.integer(foceType), 0L) && isTRUE(ef$dependsF0)
  E0List <- vector("list", length(ids))
  if (.needF0) {                                       # eta=0 population solve (frozen R0), batched
    E0List <- .foceiAnalyticSolveAll(am, th, matrix(0, length(ids), neta), .idCode, data, .obsT, solveTol)
    if (is.null(E0List)) return(NULL)
  }
  eta0Mat <- ebes
  if (.foce) for (i in seq_along(ids)) {              # FOCE EBE re-solve (per-subject Newton)
    .o <- .obsAll[[i]]
    .e0 <- .foceiAnalyticFoceEbe(am, th, ebes[i, ], .byId[[as.character(.idCode[i])]], .o$TIME, .o$DV, etav,
                                 if (is.null(E0List[[i]])) NULL else E0List[[i]]$R, Oi, neta, solveTol,
                                 foceType = foceType, cens = .o$CENS, limit = .o$LIMIT)
    if (is.null(.e0)) return(NULL)
    eta0Mat[i, ] <- .e0
  }
  .batch <- !rescale && !nzchar(Sys.getenv("FOCEI_NO_FD3_BATCH"))
  .EsAll <- NULL
  if (.batch) {
    .EsAll <- .foceiAnalyticSolveAllFD3(am, th, eta0Mat, .idCode, data, .obsT, tol = solveTol, withR = FALSE)
    if (is.null(.EsAll)) .batch <- FALSE              # batched solve failed -> per-subject fallback
  }
  for (i in seq_along(ids)) {
    s <- .byId[[as.character(.idCode[i])]]
    obs <- .obsAll[[i]]
    eta0 <- eta0Mat[i, ]
    E0 <- E0List[[i]]
    E <- if (.batch) .EsAll[[i]] else .foceiAnalyticSolveSubjectFD3(am, c(th, setNames(eta0, etav)), s, obs$TIME, tol = solveTol)
    if (is.null(E)) return(NULL)                       # solve failure -> caller falls back to FD
    # near-zero-prediction guard for a vanishing residual variance (pure proportional
    # R = sp^2 f^2 -> 0): the 1/R observed-information terms blow up, so drop to FD.
    if (isTRUE(ef$canVanish)) {
      .fa <- abs(E$f)
      if (any(!is.finite(.fa)) || min(.fa) < 1e-6 * max(.fa))
        return(.foceiAnalyticFallback("proportional error near a zero prediction"))
    }
    E$y <- .foceiAnalyticTbsY(obs$DV, E$trans)
    ehat <- eta0
    if (rescale) {
      E$a <- sweep(E$a, 2, iovDirScale, `*`)           # a_B = a_A / w on occasion directions
      for (d1 in seq_len(ndir)) for (d2 in seq_len(ndir))
        E$A[, d1, d2] <- E$A[, d1, d2] * iovDirScale[d1] * iovDirScale[d2]
      for (d1 in seq_len(neta)) for (d2 in seq_len(ndir)) for (d3 in seq_len(ndir))  # Ath is [obs,neta,ndir,ndir]
        E$Ath[, d1, d2, d3] <- E$Ath[, d1, d2, d3] * iovDirScale[d1] * iovDirScale[d2] * iovDirScale[d3]
      if (!is.null(E0)) {                              # population sensitivities share the rescaling
        E0$a <- sweep(E0$a, 2, iovDirScale, `*`)
        for (d1 in seq_len(ndir)) for (d2 in seq_len(ndir))
          E0$A[, d1, d2] <- E0$A[, d1, d2] * iovDirScale[d1] * iovDirScale[d2]
      }
      ehat <- ehat * etaScale                          # xi_hat = w * eta_hat on occasion etas
    }
    # AGQ: the quadrature nodes are placed per subject (from that subject's Ht), so they are
    # solved inside the subject assembly via this callback.  Order 2 (`SolveFA`) is enough --
    # the nodes need a/A but never the 3rd-order Ath, which is an eta-hat-only quantity.
    .sn <- if (is.null(.ag)) NULL else function(etak)
      .foceiAnalyticSolveFA(am, c(th, setNames(etak, etav)), s, obs$TIME, tol = solveTol)
    Ri <- tryCatch(.foceiAnalyticSubjectR(E, ehat, Om, ef, neta, nth, nsg, ef$sgVar, omd,
                                          ndir = ndir, dirTh = dirTh, Oi = Oi, interaction = interaction,
                                          E0 = E0, foceType = foceType,
                                          qx = .ag$x, qw = .ag$w, solveNode = .sn),
                   error = function(e) NULL)
    if (is.null(Ri) || !all(is.finite(Ri))) return(NULL)
    R <- R + Ri
  }
  R
}

#' C++ `foceiCalcR` hook for covType="analytic": from the live focei env `e`,
#' assemble the full theta+sigma+Omega observed-information R, stash the natural
#' cov (installed by [.foceiInstallAnalyticCov]), and return the theta `R.0` for
#' the native SE path -- or `NULL` to fall back to the FD Hessian.
#' @noRd
.foceiCalcRanalytic <- function(e) {
  tryCatch({
    ui <- get("ui", e)
    # covType="analytic" opt-in; anything else keeps the finite-difference Hessian
    if (!identical(rxode2::rxGetControl(ui, "covType", "fd"), "analytic")) return(NULL)
    if (!.hasRxSens()) return(.foceiAnalyticFallback("an rxode2 without symbolic sensitivities"))
    # bounded thetas are estimated on a transformed scale and the pre-final Jacobian
    # hook corrects env$cov to the natural scale; an analytic install would overwrite
    # that with the internal-scale cov -> bow out to the (Jacobian-correct) FD path.
    if (!is.null(ui$boundedTransforms) && length(ui$boundedTransforms) > 0L)
      return(.foceiAnalyticFallback("a bounded parameter transform"))
    # dose-history functions (tad/podo/tafd/tlast/tfirst/dosenum) are functions of
    # time and the dose record only -- they carry no eta/theta dependence, so rxode2
    # differentiates them to zero (.rxToSEDualVarFunction) and they no longer force
    # the finite-difference fallback.
    # scope: conditional methods only (FOCEI and both FOCE variance modes) with a
    # single Gaussian endpoint; FO/FOI and anything else -> FD.
    if (isTRUE(as.logical(rxode2::rxGetControl(ui, "fo", FALSE))))
      return(.foceiAnalyticFallback("the FO/FOI method"))
    # linCmt() has no symbolic state sensitivities for the augmented model
    if (isTRUE(any(ui$predDf$linCmt)))
      return(.foceiAnalyticFallback("a linCmt() model"))
    if (.foceiCholSECovActive(ui))
      return(.foceiAnalyticFallback(.foceiCholSECovReason))
    interaction <- as.integer(rxode2::rxGetControl(ui, "interaction", 1L))                   # 1 FOCEI / 0 FOCE
    # foceType picks the FOCE variance mode (0 "nonmem" frozen R0, 1 "foce+" live R);
    # it only matters when interaction=0 (FOCEI always uses the live conditional R).
    foceType <- if (interaction == 0L) as.integer(rxode2::rxGetControl(ui, "foceType", 0L)) else 0L
    # AGQ is supported for FOCEI (interaction=1) only -- the node terms are derived on the
    # FOCEI conditional R.  This gate is load-bearing: .foceiAnalyticSubjectR routes
    # interaction=0 to the FOCE assembler, which takes no qx/qw and would silently return the
    # nAGQ=1 FOCE covariance for an AGQ fit.  Matches the analytic gradient, which declines
    # the same combination.
    if (as.integer(rxode2::rxGetControl(ui, "nAGQ", 1L)) > 1L && interaction != 1L)
      return(.foceiAnalyticFallback("adaptive Gaussian quadrature (nAGQ > 1) without interaction"))
    ef <- .foceiAnalyticErrFull(ui)
    if (is.null(ef)) return(NULL)
    # Estimated boxCox/yeoJohnson lambda: both the FOCEI and FOCE (nonmem / foce+)
    # observed-information cov carry the DV-transform 2nd-order chain (dy'/dlambda residual
    # split in the rho DATA terms + d2y'/dlambda2 in the lambda-lambda block).
    # Both FOCEI and FOCE use the general (f,R) cov path for a general (foceiOnly)
    # variance structure; add/prop keeps the fast symbolic assembly.

    ini <- ui$iniDf
    .map <- .foceiEtaThetaMap(ui)
    etaNames <- .map$etaNames; neta <- length(etaNames)
    if (neta == 0L) return(.foceiAnalyticFallback("no random effects"))
    thetaForEta <- .map$thetaForEta

    # IOV (inter-occasion variability): nlmixr rewrites `v ~ .. | occ` BEFORE the fit
    # into a SD-scale THETA `v` (in skipCov) plus k FIXED unit-variance per-occasion
    # etas `rx.v.<occ>` entering the predictor as `abs(v) * sum(rx.v.<occ>*(occ==j))`.
    # We treat the shared variance omega_v = v^2 as ONE Omega-variance parameter over
    # those k occasion etas (reusing the "om" branch), so the occasion etas are the
    # only permitted non-mu-referenced etas.
    iovVars <- .uiIovEnv$iovVars
    # only the SD-scale IOV predictor (abs(v)*occ-eta) is implemented; var/logsd/logvar
    # use a different predictor + chain rule -> bow out to FD.
    if (length(iovVars) > 0L &&
          !identical(rxode2::rxGetControl(ui, "iovXform", "sd"), "sd"))
      return(.foceiAnalyticFallback("IOV with a non-SD parameterization"))
    occEta <- lapply(iovVars, function(v) grep(paste0("^rx\\.", gsub(".", "\\.", v, fixed = TRUE), "\\."), etaNames))
    occEtaAll <- sort(unique(unlist(occEta)))
    naEta <- which(is.na(thetaForEta))
    # a plain non-mu-ref eta keeps its own ETA_i_ direction + eta-named Omega variance;
    # occasion (IOV) etas are the grouped shared-variance case handled below.
    if (length(occEtaAll) > 0L && !all(occEtaAll %in% naEta)) return(NULL)  # sanity: occ etas must be non-mu-ref

    # converged estimates (e$theta, NOT ui$iniDf$est which holds the initials)
    thNames <- get("thetaNames", e)
    thVals  <- get("theta", e)$theta
    names(thVals) <- thNames
    # cov parameters = non-skipped thetas (skipCov drops residual/fixed/IOV)
    skip <- get("skipCov", e)[seq_along(thNames)]
    covParams <- thNames[!skip]

    # converged residual sigma into the error-model evaluator (ef reads initials)
    .valc <- setNames(as.numeric(thVals[ef$sgName]), ef$sgVar)
    ef$ev <- local({ v <- .valc; function(expr, f, y, f0 = f) eval(expr, c(list(f = f, y = y, f0 = f0), as.list(v))) })

    # scope: fixed structural thetas break the eta<->theta indexing
    if (any(.iniIsFixed(ini, thetaForEta)))
      return(.foceiAnalyticFallback("a fixed mu-referenced structural parameter"))
    keep <- !.iniIsFixed(ini, ef$sgName); ef$sgVar <- ef$sgVar[keep]; ef$sgName <- ef$sgName[keep]  # drop fixed sigma

    Om <- get("omega", e)
    # Param B: give each occasion-eta group its shared variance w^2 (the internal
    # model fixes them to unit variance and carries w=|v| in the predictor); the
    # occasion-eta sensitivities/EBEs are rescaled to match in .foceiAnalyticAssembleR.
    iovW <- vapply(iovVars, function(v) abs(as.numeric(thVals[v])), numeric(1))
    etaScale <- rep(1, neta)
    iovGroups <- vector("list", length(iovVars))
    for (g in seq_along(iovVars)) {
      w <- iovW[g]; oe <- occEta[[g]]
      if (length(oe) == 0L || !is.finite(w) || w <= 0)                # unidentified/zero IOV -> FD
        return(.foceiAnalyticFallback("an unidentified or zero IOV variance"))
      for (j in oe) Om[j, j] <- w^2
      etaScale[oe] <- w
      iovGroups[[g]] <- list(idx = oe, w = w)
    }
    pairs <- .foceiOmegaPairs(Om, ini)                # free Omega lower-triangle (declared blocks)
    # occasion-eta diagonals are grouped into the IOV param(s), so drop them from the
    # ordinary element list (they are FIXED and thus already dropped by .omegaFixed).
    if (length(occEtaAll) > 0L && nrow(pairs) > 0L)
      pairs <- pairs[!(pairs[, 1] %in% occEtaAll | pairs[, 2] %in% occEtaAll), , drop = FALSE]
    omd <- .omegaVarCovDeriv(Om, pairs, iov = if (length(iovGroups)) iovGroups else NULL)

    # Uniform direction assembly (shared with the standalone oracle).  IOV variance
    # thetas are NOT structural (predictor) thetas: they never get a THETA_j_ direction
    # (the occasion-eta sensitivities already carry the w factor).
    .dir <- .foceiAnalyticDirections(ini, thetaForEta, ef$sgName, neta, iovVars,
                                     sharedEta = unname(.foceiEtaOccurrence(ui) > 1L))
    if (is.null(.dir)) return(.foceiAnalyticFallback("an unresolvable direction set"))
    thStruct <- .dir$thStruct
    dirs <- .dir$dirs; dirTh <- .dir$dirTh; ndir <- .dir$ndir; nth <- .dir$nth
    # occasion-eta directions carry a[,occ] = w * a_B; rescale by 1/w to the variance-w^2 basis
    iovDirScale <- rep(1, ndir)
    for (g in seq_along(iovVars)) iovDirScale[occEta[[g]]] <- 1 / iovW[g]

    nsg <- length(ef$sgVar)
    # om-param order matches .omegaVarCovDeriv: ordinary Omega elements THEN the IOV
    # shared-variance params (reported on the SD scale as `v`, matching the theta name).
    onm <- etaNames                                            # Omega named by the eta (om.eta.cl)
    # R matrix param order = dirP = [dirTh (structural thetas THEN estimated-lambda) | dirSg
    # (sigmas)] then Omega.  An estimated boxCox/yeoJohnson lambda is a theta-like DIRECTION
    # (already in thStruct/dirTh), so it must NOT be re-listed as a sigma -- use the direction
    # set's sigma names (.dir$sgName, lambda-excluded), not ef$sgName, to avoid a duplicate
    # name that would make fullNm longer than Rfull.
    fullNm <- c(thStruct, .dir$sgName, .foceiOmegaCovNames(pairs, onm), iovVars)  # full natural-scale order

    # op_focei cov-params (non-skipped structural + residual thetas) must all live in
    # the full natural-scale cov (structural thetas + estimated lambda in `thStruct`, residual
    # sigmas in `.dir$sgName`); the R.0 / covFull=FALSE block is taken from it by name below.
    if (!all(covParams %in% fullNm)) return(NULL)

    th <- setNames(as.numeric(thVals[thNames]), paste0("THETA_", seq_along(thNames), "_"))
    etav <- paste0("ETA_", seq_len(neta), "_")
    etaObf <- get("etaObf", e)
    ebes <- as.matrix(etaObf[, paste0("ETA[", seq_len(neta), "]"), drop = FALSE])
    ids  <- etaObf$ID
    data <- get("dataSav", e)
    # Censored data (M2/M3/M4): the analytic observed-information cov carries the censored
    # determinant partials through the general (f,R) path (AssembleRFR), matching the exact
    # censored analytic outer gradient.  Any censored/BLOQ observation (CENS != 0 or a finite
    # LIMIT) routes to that path (and uses the same batched FD3 solve as the uncensored (f,R) cov).
    .hasCensD <- (!is.null(data$CENS) && any(data$CENS != 0, na.rm = TRUE)) ||
      (!is.null(data$LIMIT) && any(is.finite(data$LIMIT)))
    # only the laplace censored determinant stays on the FD cov; gauss (default) is analytic
    if (.hasCensD && as.integer(rxode2::rxGetControl(ui, "censOption", 0L)) == 1L)
      return(.foceiAnalyticFallback("censoring with censOption='laplace'"))

    # AGQ scope.  The node terms live only in .foceiAnalyticSubjectR, so every route that
    # leaves it returns the nAGQ=1 Laplace cov stamped covMethod="analytic" -- finite, no
    # error -- and the gate must sit above the routing, not in the assembler.  foceiOnly and
    # censored route to AssembleRFR (no qx/qw); IOV would hand the node callback xi-space
    # etas and un-rescaled a/A; agqLow/agqHi kinks the objective (inner.cpp clamps each
    # node's lik).  cholSECov is already declined at the top of this hook.
    if (as.integer(rxode2::rxGetControl(ui, "nAGQ", 1L)) > 1L) {
      .agqLo <- suppressWarnings(as.numeric(rxode2::rxGetControl(ui, "agqLow", -Inf)))
      .agqHi <- suppressWarnings(as.numeric(rxode2::rxGetControl(ui, "agqHi", Inf)))
      .why <- if (isTRUE(ef$foceiOnly)) "a general or multi-endpoint residual variance"
        else if (.hasCensD) "censored observations"
        else if (length(iovVars) > 0L) "inter-occasion variability (IOV)"
        else if (isTRUE(is.finite(.agqLo)) || isTRUE(is.finite(.agqHi))) "a finite agqLow/agqHi node clamp"
        else NULL
      if (!is.null(.why))
        return(.foceiAnalyticFallback(paste0("adaptive Gaussian quadrature (nAGQ > 1) with ", .why)))
    }

    # Full natural-scale observed-information R (theta + sigma + Omega), summed over
    # subjects.  startedEnv=e flags `.analyticStarted` before the augmented solve so
    # the C++ hook skips its finite-difference fallback if the assembly then fails.
    # FOCEI with a general (non-add/prop / multi-endpoint) variance uses the (f,R) cov
    # (any structure, sigmas as directions); add/prop FOCEI keeps the fast symbolic
    # assembly until the (f,R) cov is ported to C++ (the R version is correct but slow).
    # FOCE and IOV keep the symbolic add/prop assembly.
    # censored FOCEI must use the general (f,R) path (the fast add/prop assembler has no
    # censored partials); it also carries any general/estimated-lambda variance.
    Rfull <- if (length(iovVars) == 0L && (isTRUE(ef$foceiOnly) || .hasCensD))
      .foceiAnalyticAssembleRFR(ui, th, ebes, ids, data, Om, ef, neta, length(.dir$dirP), .dir$dirP, omd,
                                dirsCov = .dir$dirsCov, ndirCov = .dir$ndirCov,
                                startedEnv = e, solveTol = .foceiAnalyticSolveTol(ui),
                                interaction = interaction, foceType = foceType, lamDir = .dir$lamDir)
    else .foceiAnalyticAssembleR(ui, th, ebes, ids, data, Om, ef, neta, nth, nsg, omd,
                                 dirs = dirs, dirTh = dirTh, ndir = ndir,
                                 startedEnv = e, solveTol = .foceiAnalyticSolveTol(ui),
                                 iovDirScale = iovDirScale, etaScale = etaScale,
                                 interaction = interaction, foceType = foceType)
    if (is.null(Rfull)) return(NULL)
    dimnames(Rfull) <- list(fullNm, fullNm)

    # full natural-scale cov -> fit$cov; theta SEs flow through the native path
    .covNat <- tryCatch(solve(Rfull), error = function(e) NULL)
    if (is.null(.covNat)) return(NULL)
    dimnames(.covNat) <- list(fullNm, fullNm)
    # the op_focei cov-param block (structural + residual thetas, by name) decides
    # success; the native covR = Rinv path (issue #666 fix) reports these SEs directly.
    .ct <- .covNat[covParams, covParams, drop = FALSE]
    .R0 <- tryCatch(solve(.ct), error = function(e) NULL)
    if (is.null(.R0)) return(NULL)   # inversion failed -> FD, and do NOT leave a stale .analyticCov
    assign(".analyticCov", .covNat, envir = e)           # stash only after the deciding inversion
    # covFull=FALSE installs the non-skipped theta block (structural + residual),
    # matching the covType="fd" shape (skipCov drops only fixed/IOV/mixProb thetas).
    assign(".analyticThetaNames", covParams, envir = e)
    .R0
  }, error = .foceiAnalyticErrWarn(1L))
}

#' Is the rxode2 symbolic-sensitivity machinery available (scope gate)?
#' @return logical, TRUE when `rxExpandSens2_` and symengine are present
#' @noRd
.hasRxSens <- function() {
  .ns <- asNamespace("rxode2")
  exists("rxExpandSens2_", envir = .ns, inherits = FALSE) &&
    exists("rxOmegaVarCovDeriv", envir = .ns, inherits = FALSE) &&
    requireNamespace("symengine", quietly = TRUE)
}

#' Augmented-solve tolerance for the SEs: `covSolveTol` if the user set it, else
#' tightened from `sigdig` (1e-9 at the default sigdig = 3).
#' @noRd
.foceiAnalyticSolveTol <- function(ui) {
  .user <- tryCatch(rxode2::rxGetControl(ui, "covSolveTol", NULL), error = function(e) NULL)
  if (!is.null(.user) && is.finite(.user) && .user > 0) return(.user)
  # fallback tracks the package default sigdig, not a frozen literal
  .sd <- suppressWarnings(as.numeric(rxode2::rxGetControl(ui, "sigdig", 3)))
  if (!is.finite(.sd)) .sd <- 3
  max(1e-14, min(1e-8, 10^-(.sd + 6)))
}

#' Non-Cholesky Omega derivatives for the analytic Omega block: Omega^{-1} and
#' log|Omega| first/second derivatives w.r.t. the free variance-covariance
#' elements `pairs` (each row `c(a, b)`, `a >= b`), from rxode2's
#' `rxOmegaVarCovDeriv`.  Returns `nom`, `dOi` (list of dOmega^{-1}/dw), `d2Oi`
#' (list-of-lists), and `d2LD` (the log-determinant second derivatives).
#'
#' `iov` (optional) is a list of IOV shared-variance parameters; each element is
#' `list(idx = <occasion-eta Omega positions>, w = <SD-scale value>)`.  The k
#' occasion diagonals all equal ONE parameter (variance omega = w^2), so its
#' Omega derivative is the SUM of the per-element derivatives over those positions
#' (the block is diagonal, so within-group cross second derivatives are summed and
#' cross-terms to the disjoint ordinary blocks vanish).  The group is chain-ruled
#' from the variance omega=w^2 to the reported SD scale w (`dOmega/dw = 2w`,
#' `d2Omega/dw^2 = 2`), so the resulting observed-information entry is directly
#' `var(w)` -- no separate delta method.  IOV params are appended AFTER `pairs`.
#' @noRd
.omegaVarCovDeriv <- function(Om, pairs, iov = NULL) {
  d <- rxode2::rxOmegaVarCovDeriv(Om, order = 2L)
  key <- function(P) paste(P[, 1], P[, 2], sep = "-")
  ek <- key(d$elements)
  ordIdx <- match(key(pairs), ek)
  nOrd <- length(ordIdx)
  nIov <- length(iov)
  nom <- nOrd + nIov
  Z <- matrix(0, nrow(Om), ncol(Om))
  dOi <- vector("list", nom)
  d2Oi <- lapply(seq_len(nom), function(i) lapply(seq_len(nom), function(j) Z))
  d2LD <- matrix(0, nom, nom)
  # dLD (1st derivative of log|Omega|): unused by the FOCEI R (only 2nd derivatives appear
  # there), but Cov_pi(Phi_p, Phi_q) needs Phi's 1st derivative on the Omega block.
  dLD <- numeric(nom)
  # ordinary variance/covariance parameters (one per row of `pairs`)
  if (nOrd > 0L) {
    for (a in seq_len(nOrd)) {
      dOi[[a]] <- d$dOmegaInv[[ordIdx[a]]]
      for (b in seq_len(nOrd)) d2Oi[[a]][[b]] <- d$d2OmegaInv[[ordIdx[a]]][[ordIdx[b]]]
    }
    d2LD[seq_len(nOrd), seq_len(nOrd)] <- d$d2LogDet[ordIdx, ordIdx, drop = FALSE]
    dLD[seq_len(nOrd)] <- d$dLogDet[ordIdx]
  }
  # IOV shared-variance parameters (grouped + chain-ruled to the SD scale w)
  if (nIov > 0L) {
    for (g in seq_len(nIov)) {
      p <- nOrd + g
      w <- iov[[g]]$w
      gi <- match(key(cbind(iov[[g]]$idx, iov[[g]]$idx)), ek)   # (occ, occ) element rows
      sdOi <- Reduce(`+`, d$dOmegaInv[gi])                      # sum_e dOmegaInv_e
      dOi[[p]] <- 2 * w * sdOi                                  # dOmegaInv/dw
      s2 <- Z
      for (a in gi) for (b in gi) s2 <- s2 + d$d2OmegaInv[[a]][[b]]
      d2Oi[[p]][[p]] <- 4 * w^2 * s2 + 2 * sdOi                 # d2OmegaInv/dw^2
      d2LD[p, p] <- 4 * w^2 * sum(d$d2LogDet[gi, gi]) + 2 * sum(d$dLogDet[gi])
      dLD[p] <- 2 * w * sum(d$dLogDet[gi])                      # dlog|Omega|/dw
    }
  }
  list(nom = nom, dOi = dOi, d2Oi = d2Oi, d2LD = d2LD, dLD = dLD)
}

#' Symbolic error machinery (rho/p f-derivatives + sigma partials) for one
#' Gaussian add/prop endpoint; `NULL` for anything else (multi-endpoint, lnorm,
#' propT/propF, a DV transform, combined1, or pure proportional) -> FD fallback.
#' @noRd
.foceiAnalyticErrFull <- function(ui) {
  # Distribution guard (Gaussian-only).  The analytic gradient AND covariance assemble a
  # conditional-GAUSSIAN observed information (rho = 0.5*((y-f)^2/R + log R)).  A non-normal
  # endpoint (t, cauchy, poisson, binomial, ordinal, ...) is built on the SAME norm-form
  # rx_pred_/rx_r_ as a Gaussian fit -- its real (heavier-tailed / discrete) density lives in a
  # SEPARATE ..Llik model the augmented sensitivity model never touches.  Such an endpoint that
  # still carries an add()/prop() location scale would slip past the err-based checks below and
  # be assembled as Gaussian (silently wrong).  The gradient is already shielded upstream by the
  # needOptimHess -> fast=FALSE downgrade (focei.R), but covType="analytic" has no such
  # downgrade, so gate it here.  "norm" and its explicit-Gaussian spelling "dnorm" (used for the
  # Gaussian endpoints of a mixed model) stay in scope; everything else falls back to FD.
  .dist <- tryCatch(as.character(ui$predDfFocei$distribution), error = function(e) NULL)
  if (!length(.dist)) .dist <- tryCatch(as.character(ui$predDf$distribution), error = function(e) character(0))
  if (length(.dist) && !all(.dist %in% c("norm", "dnorm")))
    return(.foceiAnalyticFallback("a non-normal likelihood endpoint"))
  # Multiple modeled endpoints: rx_pred_ and rx_r_ are single dvid-conditional
  # expressions that already select the right endpoint per observation when solved
  # against the dataset, so the (f,R) path handles them -- but the single-endpoint
  # symbolic add/prop machinery (one rx_pred_) does not, so force the general path.
  .multiEndpoint <- !is.null(ui$predDf) && nrow(ui$predDf) > 1L
  ini <- ui$iniDf
  # Both-sides transforms: rx_pred_ is already the transformed prediction and rx_r_ the
  # transformed-scale variance, so the analytic path only transforms the DV to that scale
  # (y' = tbs(DV), via the C++ _powerD in the assemblers).  Exact for a STATIC transform
  # (lnorm, or a FIXED boxCox/yeoJohnson lambda): the transform Jacobian -2 log|dy'/dy| is
  # then constant (or linear in a fixed lambda) and drops out of the observed information.
  # An ESTIMATED lambda additionally makes the DV move with the parameter (a dy'/dlambda
  # chain in the residual): the GRADIENT carries it (lambda is a theta-like direction with
  # df'/dlambda plus the residual DV chain and the -2 log|J| Jacobian term); the cov path
  # gates estimated lambda separately (its 2nd-order lambda terms are not yet ported).
  .trans <- as.character(ui$predDf$transform)
  .estLam <- FALSE
  if (!all(.trans == "untransformed")) {
    .lamRows <- ini[!is.na(ini$err) & ini$err %in% c("boxCox", "yeoJohnson"), , drop = FALSE]
    .estLam <- any(!.lamRows$fix)
    # Multiple estimated per-endpoint lambdas need an endpoint->lambda DV mapping that is not
    # wired (the gradient gates this at foceiGradAnalytic.R; without the same gate the covariance
    # would recycle one endpoint's dy'/dlambda chain across every lambda column -- silently
    # wrong).  A single estimated lambda is the ported case; more than one falls back to FD.
    if (sum(!.lamRows$fix) > 1L)
      return(.foceiAnalyticFallback("multiple estimated boxCox/yeoJohnson lambdas"))
  }
  er <- ini[!is.na(ini$err), , drop = FALSE]
  if (nrow(er) == 0L) {
    # All residual-error parameters are fixed: the fit's working iniDf drops
    # them (they are baked into the compiled rx_r_ as constants), yet the model
    # still HAS a residual (predDf carries its errType).  Rather than fall back
    # to FD -- which re-solves the ODE per outer parameter (all omegas here,
    # which do not even enter the ODE) -- run the general (f,R) path with NO
    # free sigma direction (nsg = 0): rx_r_ is read from the solve, the omega
    # gradient uses the precomputed Omega derivatives, and R's constancy needs
    # no re-solve.  Only the plain conditional-Gaussian error types are ported;
    # anything else stays on FD.
    .et <- as.character(ui$predDf$errType)
    .tok <- unique(trimws(unlist(strsplit(.et, "[^A-Za-z0-9]+"))))
    .tok <- .tok[nzchar(.tok)]
    .known <- c("add", "prop", "pow", "combined1", "combined2")
    if (length(.tok) == 0L || !all(.tok %in% .known) ||
          !all(.trans == "untransformed"))
      return(.foceiAnalyticFallback("no estimated residual error"))
    .hasAddFloor <- any(.tok %in% c("add", "combined1", "combined2"))
    return(list(sgVar = character(0), sgName = character(0), sc = NULL,
                per = NULL, pair = NULL, foce = NULL, focePlus = NULL,
                foceiOnly = TRUE, canVanish = !.hasAddFloor,
                dependsF0 = !all(.tok == "add"), estLam = .estLam,
                ev = function(e, f, y, f0 = f) NULL))
  }
  sgNameAll <- er$name                               # ALL error params (excluded from directions)
  # pure proportional / power error (no additive floor) vanishes as f -> 0, making the
  # 1/R observed-information terms blow up near zero predictions; the assembly guards it.
  # lnorm/logitNorm/probitNorm are additive on their transformed (log / logit / probit) scale
  # (R = sd^2 constant), so the variance never vanishes even though the transformed prediction
  # crosses zero (e.g. logit(cp) = 0 at the bound midpoint) -- the canVanish guard, meant for
  # pure proportional/power error whose R -> 0 with f, must not fire for them.
  canVanish <- !any(er$err %in% c("add", "lnorm", "logitNorm", "probitNorm"))
  # model-declared addProp wins; the control applies only when the model says "default"
  addPr <- as.character(ui$predDf$addProp)
  if (length(addPr) != 1L || is.na(addPr) || addPr == "default") {
    addPr <- tryCatch(rxode2::rxGetControl(ui, "addProp", "combined2"), error = function(e) "combined2")
  }
  # The general (f,R) path reads R and dR/dsigma from the model's own rx_r_ via the solve, so ANY
  # conditional-Gaussian variance structure is in scope -- including combined1 (sa+sp*f)^2, which
  # rx_r_ carries directly ((add + pred*prop)^2, the weight squared).  The symbolic add/prop path
  # below is only the FAST route for the plain untransformed add/prop/combined2 case; everything
  # else returns a minimal `ef` (foceiOnly=TRUE) and runs the general (f,R) assembler
  # (.foceiAnalyticAssembleRFR), which handles BOTH FOCEI and FOCE (interaction 0/1).
  # a both-sides transform is never the plain-scale symbolic add/prop machinery: force the
  # general (f,R) path (rx_pred_/rx_r_ carry the transform; only the DV is retransformed).
  .isAddProp <- !.multiEndpoint && all(.trans == "untransformed") &&
    all(er$err %in% c("add", "prop")) && !identical(addPr, "combined1")
  # R0 (FOCE nonmem frozen variance) needs the eta=0 population solve only when R depends
  # on the prediction (any non-additive error term); pure additive R is constant.
  .dependsF0 <- !all(er$err == "add")
  if (!.isAddProp)
    return(list(sgVar = character(0), sgName = sgNameAll, sc = NULL, per = NULL, pair = NULL,
                foce = NULL, focePlus = NULL, foceiOnly = TRUE, canVanish = canVanish,
                dependsF0 = .dependsF0, estLam = .estLam, ev = function(e, f, y, f0 = f) NULL))
  addN <- er$name[er$err == "add"]; propN <- er$name[er$err == "prop"]
  hasA <- length(addN) == 1L; hasP <- length(propN) == 1L
  if (!hasA && !hasP)
    return(.foceiAnalyticFallback("no additive or proportional residual error"))
  Rstr <- if (hasA && hasP) "sa^2+sp^2*f^2" else if (hasP) "sp^2*f^2" else "sa^2"
  Rq <- parse(text = Rstr)[[1]]
  rhoE <- bquote(0.5 * ((y - f)^2 / .(Rq) + log(.(Rq))))
  pE <- bquote(1 / .(Rq) + 0.5 * (.(stats::D(Rq, "f")) / .(Rq))^2)
  # FOCE (interaction=0) pieces: the inner problem drops dR/deta, so its gradient
  # coefficient is the least-squares part only, qE = drho/df = -(y-f)/R
  # (q'=dq/df, q''=d2q/df2 build the FOCE inner Hessian Hf and its 3-tensor); the FOCE
  # Laplace determinant curvature is pFE = 1/R (no 0.5*(R'/R)^2 interaction term).
  # R here is written in the live prediction f, so D(, "f") carries the dR/df chain:
  # these are the "foce+" (live conditional R) pieces; the "nonmem" frozen-R0 variants
  # are built below with the separate f0 symbol.
  qE <- bquote(-(y - f) / .(Rq))
  pFE <- bquote(1 / .(Rq))
  DD <- function(e, ...) { for (v in c(...)) e <- stats::D(e, v); e }
  sgVar <- c(if (hasA) "sa", if (hasP) "sp"); sgName <- c(if (hasA) addN, if (hasP) propN)
  val <- setNames(c(if (hasA) er$est[er$name == addN], if (hasP) er$est[er$name == propN]), sgVar)
  # r0 is rho itself: unused by the FOCEI R, but pi_k is a softmax over Phi(eta_k) and the
  # cov-path solve carries f/a/A with no $R to rebuild it from.
  sc <- list(r0 = rhoE,
             r1 = DD(rhoE, "f"), r2 = DD(rhoE, "f", "f"), r3 = DD(rhoE, "f", "f", "f"),
             p = pE, p1 = DD(pE, "f"), p2 = DD(pE, "f", "f"),
             q0 = qE, q1 = DD(qE, "f"), q2 = DD(qE, "f", "f"),
             pF = pFE, pF1 = DD(pFE, "f"), pF2 = DD(pFE, "f", "f"))
  per <- list(); for (s in sgVar) per[[s]] <- list(rf = DD(rhoE, "f", s), rff = DD(rhoE, "f", "f", s),
    ps = DD(pE, s), pf = DD(pE, "f", s),
    qs = DD(qE, s), qsf = DD(qE, "f", s), psF = DD(pFE, s), pfF = DD(pFE, "f", s),
    rs = DD(rhoE, s))   # pure d(rho)/d(sigma): first-derivative (outer-gradient) sigma term
  pair <- list(); for (i in seq_along(sgVar)) for (j in i:length(sgVar)) { a <- sgVar[i]; b <- sgVar[j]
    pair[[paste0(a, b)]] <- list(rss = DD(rhoE, a, b), rfss = DD(rhoE, "f", a, b), pss = DD(pE, a, b),
      qss = DD(qE, a, b), pssF = DD(pFE, a, b)) }
  # CORRECTED FOCE variance R0 = R at the eta=0 POPULATION prediction (symbol `f0`),
  # while the residual/numerator keeps the eta-hat prediction `f` (the estimator freezes
  # R at eta=0: getPopR/likInner0).  Because f0 is a separate symbol, D(.,"f") treats R0
  # as constant (so q1=1/R0, q2=0, pF1=0, r3=0), and D(.,f0) gives the theta-chain (R0
  # varies with theta through the population prediction) supplied via a0/A0 in the
  # subject assembly.  Additive R0=sa^2 has no f0 -> all f0 fields vanish (== FOCEI-add).
  R0q <- parse(text = gsub("\\bf\\b", "f0", Rstr))[[1]]
  rho0E <- bquote(0.5 * ((y - f)^2 / .(R0q) + log(.(R0q))))
  q0E <- bquote(-(y - f) / .(R0q)); pF0E <- bquote(1 / .(R0q))
  scF <- list(r1 = DD(rho0E, "f"), r2 = DD(rho0E, "f", "f"), r3 = DD(rho0E, "f", "f", "f"),
              q0 = q0E, q1 = DD(q0E, "f"), q2 = DD(q0E, "f", "f"),
              pF = pF0E, pF1 = DD(pF0E, "f"), pF2 = DD(pF0E, "f", "f"))
  perF <- list(); for (s in sgVar) perF[[s]] <- list(rf = DD(rho0E, "f", s), rff = DD(rho0E, "f", "f", s),
    qs = DD(q0E, s), qsf = DD(q0E, "f", s), psF = DD(pF0E, s), pfF = DD(pF0E, "f", s),
    rs = DD(rho0E, s))   # pure d(rho)/d(sigma) for the FOCE first-derivative sigma term
  pairF <- list(); for (i in seq_along(sgVar)) for (j in i:length(sgVar)) { a <- sgVar[i]; b <- sgVar[j]
    pairF[[paste0(a, b)]] <- list(rss = DD(rho0E, a, b), rfss = DD(rho0E, "f", a, b),
      qss = DD(q0E, a, b), pssF = DD(pF0E, a, b)) }
  f0F <- list(qf0 = DD(q0E, "f0"), qff0 = DD(q0E, "f", "f0"), qf0f0 = DD(q0E, "f0", "f0"),
              pFf0 = DD(pF0E, "f0"), pFf0f0 = DD(pF0E, "f0", "f0"),
              rhof0 = DD(rho0E, "f0"), rhof0f0 = DD(rho0E, "f0", "f0"))
  perf0F <- list(); for (s in sgVar) perf0F[[s]] <- list(qf0s = DD(q0E, "f0", s),
    pFf0s = DD(pF0E, "f0", s), rhof0s = DD(rho0E, "f0", s))
  foce <- list(sc = scF, per = perF, pair = pairF, f0 = f0F, perf0 = perf0F, dependsF0 = hasP)
  # "foce+" (foceType=1, live conditional R): the same truncated inner gradient and
  # determinant, but R follows the eta-hat prediction, so the live-R sc/per/pair pieces
  # apply as-is and there is no f0 population solve/chain.
  focePlus <- list(sc = sc, per = per, pair = pair, dependsF0 = FALSE)
  list(sgVar = sgVar, sgName = sgName, sc = sc, per = per, pair = pair, foce = foce,
       focePlus = focePlus, foceiOnly = FALSE, canVanish = canVanish, dependsF0 = .dependsF0,
       estLam = .estLam,
       ev = function(e, f, y, f0 = f) eval(e, c(list(f = f, y = y, f0 = f0), as.list(val))))
}

#' Augmented rxode2 model with state + 1st/2nd-order sensitivities and the
#' prediction chain f1/f2 over an arbitrary direction set `dirs` (each ETA_i_ or
#' THETA_j_).  State sensitivities use rxode2's `.rxSens`; the higher-order
#' prediction chain f1/f2 is built here.  O(ndir^2), so it compiles far more
#' cheaply than a 3rd-order model.
#' @return list(augMod, dirs, ndir, st, P2) or `NULL` on failure
#' @noRd
.foceiAnalyticAugCache <- new.env(parent = emptyenv())     # session cache: (model digest | dirs) -> aug model
#' Solve-output column names + index maps for the augmented model (from dirs, P2, sigTh).
#' Precomputed once on `am$cols` at build time so the matrix-extraction readers avoid per-call
#' paste0; recomputed on the fly for a reconstructed `am` (from foceiModel$outer) that predates it.
#' @noRd
.foceiAnalyticCols <- function(dirs, fDirs, P2, P2r, sigTh) {
  sigP2 <- if (length(sigTh)) do.call(rbind, lapply(seq_along(sigTh), function(.a)
    data.frame(a = .a, b = seq_along(sigTh)[seq_along(sigTh) >= .a]))) else NULL
  # f-part (prediction a/A) spans f-directions only (sigma slots are zero-filled by the reader);
  # R-part (aR/AR) spans every direction.  iiF/jjF index the f2 pairs into the FULL dirs vector.
  list(f1 = paste0("rx_f1_", fDirs), fDirIdx = match(fDirs, dirs),
       # the nrow() guards are load-bearing: paste0 treats a zero-length argument as "", so an
       # empty P2 (order = 1) gives "rx_f2__" rather than character(0) -- a bogus column name
       # that makes the solve's column check fail and every node solve return NULL.
       f2 = if (nrow(P2)) paste0("rx_f2_", P2$i, "_", P2$j) else character(0),
       iiF = match(P2$i, dirs), jjF = match(P2$j, dirs),
       rvar1 = paste0("rx_rvar1_", dirs),
       rvar2 = if (nrow(P2r)) paste0("rx_rvar2_", P2r$i, "_", P2r$j) else character(0),
       ii = match(P2r$i, dirs), jj = match(P2r$j, dirs),
       rsig = if (length(sigTh)) paste0("rx_rsig_", sigTh, "_") else character(0),
       rsig1 = lapply(sigTh, function(.n) paste0("rx_rsig1_", .n, "_", dirs)),
       rsig2 = if (is.null(sigP2)) character(0) else paste0("rx_rsig2_", sigTh[sigP2$a], "_", sigTh[sigP2$b]),
       sigP2 = sigP2)
}
.foceiAnalyticEtCache <- new.env(parent = emptyenv())      # per-fit translated event table (etTrans reuse)
#' The augmented model's events are IDENTICAL across every solve of a fit (only the theta/eta
#' params vary), so translate them once with etTrans and reuse -- ~40% off each population solve
#' (validated: identical predictions).  Keyed by the model key + a content hash of the data
#' (`digest`, already a dependency), so distinct datasets that happen to share summary statistics
#' cannot collide onto the same translated event table; falls back to raw data.  The hash is far
#' cheaper than the etTrans it guards, so the per-solve lookup stays a net win.
#' @noRd
.foceiAnalyticEvents <- function(am, data) {
  if (is.null(am$key)) return(data)
  .fp <- tryCatch(digest::digest(data), error = function(e) NULL)
  if (is.null(.fp)) return(data)                          # cannot fingerprint -> translate fresh (no cache)
  .ek <- paste0(am$key, "|et|", nrow(data), "|", .fp)
  .et <- get0(.ek, envir = .foceiAnalyticEtCache, inherits = FALSE)
  if (is.null(.et)) {
    .et <- tryCatch(rxode2::etTrans(data, am$augMod), error = function(e) NULL)
    if (is.null(.et)) return(data)
    if (length(ls(.foceiAnalyticEtCache, all.names = TRUE)) >= 256L)    # bound memory in long sessions
      rm(list = ls(.foceiAnalyticEtCache, all.names = TRUE), envir = .foceiAnalyticEtCache)
    assign(.ek, .et, envir = .foceiAnalyticEtCache)
  }
  .et
}
#' Covariate-coefficient directions that can be reconstructed by eta-scaling.
#'
#' A mu-referenced covariate coefficient `b` enters linearly through an eta's parameter
#' (`cl = exp(tcl + eta.cl + b*cov)`), so for a subject-constant covariate every sensitivity
#' of `b` equals the linked eta's scaled by the covariate value: df/db = cov*df/deta,
#' d2f/(db dx) = cov*d2f/(deta dx), d2f/db2 = cov^2*d2f/deta2 (verified exact).  We can therefore
#' skip integrating `b`'s state-sensitivity ODEs and emit its f1/f2/rvar/rsig columns as scaled
#' copies of the linked eta's columns.  Returns a named list `covDir -> list(etaDir, scale)`
#' (scale = the covariate expression in rxode2 syntax) for such thetas present in `fDirs`.
#'
#' Detection reuses rxode2's OWN mu-reference classification rather than re-parsing the model --
#' the same machinery the mu-referenced (muModel="lin"/"irls") family consumes (see #711/#712):
#'   * ui$muRefDataFrame                  : the theta<->eta link.  A well-formed mu-reference has
#'       the eta as a BARE `+eta` in the parameter; a scaled eta (`0.5*eta`) is simply absent, so
#'       such coefficients are auto-excluded with no extra check.
#'   * ui$muRefCovariateDataFrame         : BARE data-column covariate coefficients.
#'   * ui$mu2RefCovariateReplaceDataFrame : ALGEBRAIC covariate-expression coefficients (allometric
#'       `log(WT/70)`, centered `(WT-70)`, ...) -- exactly the set `muRefCovAlg` rewrites; its
#'       `covariate` column is the (equivalent) covariate expression used as the scale.
#' Both covariate frames carry (theta, covariateParameter, covariate); the theta is mapped to its
#' eta via muRefDataFrame.  The ONE structural condition rxode2's frames do not encode is an eta
#' REUSED in another parameter (then df/deta gains a path df/db lacks), so an eta-occurrence guard
#' is applied.  The scale being covariate-and-constant only is enforced downstream by the caller's
#' subject-constant gate.  Empty when the info is unavailable.
#' @noRd
.foceiAnalyticCovDirMap <- function(ui, fDirs) {
  .ini <- ui$iniDf
  # eta name -> ETA_k_ (k = order among diagonal etas, matching .foceiAnalyticDirections' etaDirs)
  .etaRows <- .ini[!is.na(.ini$neta1) & .ini$neta1 == .ini$neta2, , drop = FALSE]
  .etaRows <- .etaRows[order(.etaRows$neta1), , drop = FALSE]
  .eta2Dir <- stats::setNames(paste0("ETA_", seq_len(nrow(.etaRows)), "_"), .etaRows$name)
  # id-level theta -> eta (occasion/IOV etas excluded -- IOV is out of analytic scope anyway).
  # May be empty: a covariate can also scale an ETA-LESS structural parameter (handled below), so
  # do NOT bail on an empty eta map -- only bail when there are no covariate coefficients at all.
  .mrd <- tryCatch(ui$muRefDataFrame, error = function(e) NULL)
  .theta2eta <- if (!is.null(.mrd) && is.data.frame(.mrd) && nrow(.mrd) > 0L) {
    .idRef <- .mrd[is.na(.mrd$level) | .mrd$level == "id", , drop = FALSE]
    stats::setNames(as.character(.idRef$eta), as.character(.idRef$theta))
  } else character(0)
  # (theta, coefficient, scale) rows from BOTH rxode2 covariate frames (bare + algebraic)
  .grab <- function(df, scaleCol) {
    if (is.null(df) || !is.data.frame(df) || nrow(df) == 0L ||
        !all(c("theta", "covariateParameter", scaleCol) %in% names(df))) return(NULL)
    data.frame(theta = as.character(df$theta), coef = as.character(df$covariateParameter),
               scale = as.character(df[[scaleCol]]), stringsAsFactors = FALSE)
  }
  .cov <- do.call(rbind, list(
    .grab(tryCatch(ui$muRefCovariateDataFrame, error = function(e) NULL), "covariate"),
    .grab(tryCatch(ui$mu2RefCovariateReplaceDataFrame, error = function(e) NULL), "covariate")))
  if (is.null(.cov) || nrow(.cov) == 0L) return(list())
  # A covariate coefficient's sensitivities are those of the parameter position it scales, times the
  # covariate value.  That position is either a random effect (`cl = exp(tcl + eta.cl + b*cov)` ->
  # df/db = cov*df/deta) or, for a covariate on an ETA-LESS parameter (`v = exp(tv + b*cov)`, no
  # eta.v), the structural theta itself (df/db = cov*df/dtheta, since theta and b enter the mu
  # identically -- tv's own direction is already integrated, so b reuses it for free).  Guard the
  # ONE structural condition the rxode2 frames don't encode: the reused position must occur exactly
  # once; a shared eta/theta gains a state path df/db lacks (counts are conservative -> `2L` when
  # lstExpr is unavailable, so uncertain cases fall back to the full symbolic build).
  .eo <- .foceiEtaOccurrence(ui)
  .to <- .foceiNameOccurrence(ui, unique(.cov$theta))
  .res <- list(); .dup <- character(0)
  for (.r in seq_len(nrow(.cov))) {
    .thN <- .cov$theta[.r]
    .etN <- if (.thN %in% names(.theta2eta)) .theta2eta[[.thN]] else NA_character_
    if (!is.na(.etN)) {
      .ed <- .eta2Dir[[.etN]]; if (is.null(.ed) || !(.ed %in% fDirs)) next
      if (!identical(.eo[[.etN]], 1L)) next                      # eta shared/reused -> scaling not exact
    } else {
      # eta-less parameter: reuse the structural theta's own (already integrated) direction
      .snt <- .ini$ntheta[.ini$name == .thN]; if (length(.snt) != 1L || is.na(.snt)) next
      .ed <- paste0("THETA_", .snt, "_"); if (!(.ed %in% fDirs)) next
      if (!identical(.to[[.thN]], 1L)) next                      # structural theta shared -> scaling not exact
    }
    .nt <- .ini$ntheta[.ini$name == .cov$coef[.r]]; if (length(.nt) != 1L || is.na(.nt)) next
    .cd <- paste0("THETA_", .nt, "_"); if (!(.cd %in% fDirs) || .cd %in% .dup) next
    # `etaDir` names the source direction to scale (an eta direction, or a structural-theta
    # direction for an eta-less parameter); the downstream emitter treats both identically.
    .new <- list(etaDir = .ed, scale = .cov$scale[.r])
    if (!is.null(.res[[.cd]])) {
      # same coefficient re-listed identically (bare frame + algebraic frame overlap) -> no-op; a
      # different (etaDir, scale) is a genuine ambiguity (b*cov1 + b*cov2) -> drop, fall back.
      if (identical(.res[[.cd]], .new)) next
      .res[[.cd]] <- NULL; .dup <- c(.dup, .cd); next
    }
    .res[[.cd]] <- .new
  }
  .res
}
#' @param order sensitivity tier: 2 (default) emits the 1st and 2nd-order chains; 1 emits
#'   the 1st-order chain only, for a consumer needing just a/aR (the AGQ quadrature nodes
#'   never read A/AR/RsigDir).  With `order = 1` the reader zero-fills A/AR naturally:
#'   `P2`/`P2r` are empty, so `cols$f2`/`cols$rvar2` are `character(0)` and the A/AR fill
#'   loops never execute.
#' @noRd
.foceiAnalyticAugModelDirs <- function(ui, dirs, order = 2L) {
  order <- as.integer(order)
  .key <- tryCatch(paste0(rxUiGet.foceiModelDigest(list(ui)), "|", paste(dirs, collapse = ","),
                          "|o", order,                                  # sensitivity tier -> distinct cached model
                          "|sk", Sys.getenv("FOCEI_NO_SIGMA_SKIP"),     # skip flag -> distinct cached model
                          # subject-constant covariate set -> which covariate directions get eta-scaling reuse
                          "|cc", paste(sort(tryCatch(rxode2::rxGetControl(ui, "foceiConstCovs", NULL),
                                                     error = function(e) NULL)), collapse = ",")),
                   error = function(e) NULL)
  if (!is.null(.key)) { .hit <- get0(.key, envir = .foceiAnalyticAugCache, inherits = FALSE)
    if (!is.null(.hit)) return(.hit) }
  .res <- tryCatch({
    .s <- ui$loadPruneSens
    .st <- rxode2::rxStateOde(.s)
    # matExp()/indLin(): the base ODEs are materialized into the pruned env (via
    # .rxJacobian -> .rxInjectMatExpOdes) with correct RHS, but rxStateOde() can return
    # them reversed (an indLin() state parses as compartment 1).  The augmented model
    # numbers compartments by first d/dt() appearance, so emitting .baseOde in that
    # reversed order misplaces default dosing (the dose lands in the wrong compartment)
    # and the eta sensitivities collapse.  Reorder source-first, matching the inner model
    # (.rxInjectMatExpDdt / rxUiGet.foceiCmtPreModel).
    .mvS <- rxode2::rxModelVars(.s)
    if (is.list(.mvS$indLin) && length(.mvS$indLin) == 4L) {
      .st <- .rxMatExpStateOrder(.st, ls(envir = .s, all.names = TRUE))
    }
    # Sigma-skip (default ON; set FOCEI_NO_SIGMA_SKIP=1 to opt out): the residual-error (sigma) directions have
    # df/dsigma=0, so their f-sensitivity columns (f1/f2) and state-sensitivities are all
    # ZERO -- pure wasted compile (each sigma direction costs a full structural parameter's
    # worth of gcc + states).  Build the prediction chain f1/f2 + state sensitivities over
    # the model-affecting directions only (.fDirs); the variance chain still spans EVERY
    # direction (sigma R-derivatives are algebraic, no state chain).  The reader zero-fills
    # the sigma a/A slots (the assembly already treats a=A=Ath=0 for a sigma direction).
    .erN <- ui$iniDf$ntheta[!is.na(ui$iniDf$err) & !(ui$iniDf$err %in% c("boxCox", "yeoJohnson"))]
    .sigDirs <- if (!nzchar(Sys.getenv("FOCEI_NO_SIGMA_SKIP"))) intersect(dirs, paste0("THETA_", .erN, "_")) else character(0)
    .fDirs <- dirs[!(dirs %in% .sigDirs)]
    # Covariate-coefficient reuse (always on; the subject-constant gate below keeps it exact): a
    # mu-ref covariate direction's sensitivities are the linked eta's scaled by the covariate value,
    # so we build state sensitivities over the model directions .mfDirs only (no covariate-direction
    # ODEs) and emit the covariate columns algebraically below.  The covariate directions stay in the
    # full direction set (.fDirs/dirs/.P2/.P2r/am$cols) so the assembly and FD3 are unchanged -- they
    # read the algebraically-emitted covariate columns like any other.
    .covMap <- .foceiAnalyticCovDirMap(ui, .fDirs)
    # Subject-constant gate: the eta-scaling identity holds only for covariates that are constant
    # within each subject (a time-varying covariate makes d(state)/db != cov*d(state)/deta).  The
    # constant-covariate set is computed from the data and stashed on the control by the fit setup;
    # reuse only a covariate whose scale expression uses exclusively constant covariates (and drop
    # all reuse when the set is unavailable, e.g. standalone use -- the full build is always correct).
    if (length(.covMap)) {
      .constCovs <- tryCatch(rxode2::rxGetControl(ui, "foceiConstCovs", NULL), error = function(e) NULL)
      .covMap <- .covMap[vapply(.covMap, function(.m) {
        .v <- all.vars(parse(text = .m$scale)[[1]]); length(.v) > 0L && all(.v %in% .constCovs)
      }, logical(1))]
    }
    .covDirs <- names(.covMap)
    .mfDirs <- .fDirs[!(.fDirs %in% .covDirs)]      # directions that get integrated state sensitivities
    rxode2::.rxJacobian(.s, c(.st, .mfDirs))
    # 1st-order sensitivities are already expanded for the gradient (free); if
    # unavailable the model is not differentiable, so bail before the 2nd-order build.
    # Gated on there being STATES to expand: a purely algebraic model (no d/dt(), e.g. a
    # generalized ll() endpoint) legitimately has no state sensitivities, and needs none --
    # .g1/.g2 below reduce to plain symbolic derivatives when .st is empty, which is exactly
    # right.  Bailing on the empty expansion made every such model fall back to finite
    # differences.
    .s1 <- rxode2::.rxSens(.s, .mfDirs)
    if (length(.st) > 0L && length(.s1) == 0L) return(NULL)
    # 2nd order (model f-directions only): the expensive expansion.  order = 1 skips it
    # entirely -- no 2nd-order state-sensitivity compartments, no f2/rvar2 chains.
    .s2 <- if (order >= 2L) rxode2::.rxSens(.s, .mfDirs, .mfDirs) else character(0)
    .pred <- get("rx_pred_", .s)
    # residual variance rx_r_ (any structure) with its own 1st/2nd sensitivity chains,
    # exactly like the prediction -- so R and dR/ddir, d2R/ddir2 come from the SOLVE
    # (matching the inner ODE model's d(R)/d(eta)); the assembly then treats the
    # transformed prediction f and the variance R as independent solved quantities.
    .rvar <- tryCatch(get("rx_r_", .s), error = function(e) NULL)
    .Dn <- function(.e, .v) symengine::D(.e, symengine::S(.v))
    .sn1 <- function(.j, ...) symengine::S(paste0("rx__sens_", .j, "_BY_", paste(c(...), collapse = "_BY_"), "__"))
    .toRx <- function(.l) rxode2::rxFromSE(.l)
    # state-sensitivity chains only for model directions .mfDirs (a sigma or covariate-reuse
    # direction contributes no integrated state chain -- sigma has d(state)/dsigma=0, and a
    # covariate direction's columns are emitted by scaling below, not via .g1/.g2).
    .g1 <- function(.ex, .p) { .e <- .Dn(.ex, .p)
      if (.p %in% .mfDirs) for (.j in .st) .e <- .e + .Dn(.ex, .j) * .sn1(.j, .p); .e }
    .g2 <- function(.ex, .p, .q) { .gq <- .g1(.ex, .q); .e <- .Dn(.gq, .p)
      if (.p %in% .mfDirs) for (.k in .st) .e <- .e + .Dn(.gq, .k) * .sn1(.k, .p)
      if (.p %in% .mfDirs && .q %in% .mfDirs) for (.j in .st) .e <- .e + .Dn(.ex, .j) * .sn1(.j, .p, .q); .e }
    # Covariate-direction reuse: map a direction to its model substitute (the linked eta for a
    # covariate direction) + scale; emit a 1st/2nd-order line symbolically when it involves only
    # model directions, else algebraically from the (eta-substituted, reordered) model column
    # scaled by the covariate value(s).  `.emitScaled1` does the 1st-order (single-direction)
    # columns; `.emitPair` the 2nd-order (pair) columns.
    .covSub <- function(.d) { .m <- .covMap[[.d]]; if (is.null(.m)) list(d = .d, s = NULL) else list(d = .m$etaDir, s = .m$scale) }
    .emitScaled1 <- function(.pre, .suf, .cd) paste0(.pre, .cd, .suf, "=(", .covMap[[.cd]]$scale, ")*", .pre, .covMap[[.cd]]$etaDir, .suf)
    .emitPair <- function(.pre, .i, .j, .symFn, .ord) {
      .si <- .covSub(.i); .sj <- .covSub(.j)
      if (is.null(.si$s) && is.null(.sj$s)) return(paste0(.pre, .i, "_", .j, "=", .symFn(.i, .j)))
      .a <- .si$d; .b <- .sj$d
      if (match(.a, .ord) > match(.b, .ord)) { .t <- .a; .a <- .b; .b <- .t }
      paste0(.pre, .i, "_", .j, "=", paste(sprintf("(%s)", c(.si$s, .sj$s)), collapse = "*"), "*", .pre, .a, "_", .b)
    }
    # f2_i_j == f2_j_i, so only the i<=j triangle is differentiated/compiled (~2x off
    # the dominant model-build cost); the reader mirrors A[,i,j]=A[,j,i].
    .P2 <- expand.grid(i = .fDirs, j = .fDirs, stringsAsFactors = FALSE)     # prediction f2: f-directions only
    .P2 <- .P2[match(.P2$i, .fDirs) <= match(.P2$j, .fDirs), , drop = FALSE]
    .P2r <- expand.grid(i = dirs, j = dirs, stringsAsFactors = FALSE)        # variance rvar2: every direction
    .P2r <- .P2r[match(.P2r$i, dirs) <= match(.P2r$j, dirs), , drop = FALSE]
    # order = 1: no 2nd-order pairs at all -- empties .fL2/.rL2, the 2nd-order IC block,
    # and cols$f2/cols$rvar2, so the reader zero-fills A/AR.
    if (order < 2L) { .P2 <- .P2[0L, , drop = FALSE]; .P2r <- .P2r[0L, , drop = FALSE] }
    # model directions first (they define the columns the covariate columns scale), covariate
    # directions after -- rxode2 evaluates assignments in order.
    .fL1 <- c(vapply(.mfDirs, function(.p) paste0("rx_f1_", .p, "=", .toRx(.g1(.pred, .p))), character(1)),
              vapply(.covDirs, function(.cd) .emitScaled1("rx_f1_", "", .cd), character(1)))
    .p2mod <- !(.P2$i %in% .covDirs) & !(.P2$j %in% .covDirs)
    .fL2 <- c(vapply(which(.p2mod), function(.r) .emitPair("rx_f2_", .P2$i[.r], .P2$j[.r], function(.a, .b) .toRx(.g2(.pred, .a, .b)), .fDirs), character(1)),
              vapply(which(!.p2mod), function(.r) .emitPair("rx_f2_", .P2$i[.r], .P2$j[.r], NULL, .fDirs), character(1)))
    # Variance rx_r_ and its 1st/2nd direction sensitivities.  The `rx_..._` naming
    # matters: rxode2 keeps a constant assignment to an `rx_<name>_` variable as a
    # real lhs output column (like the inner model's rx__sens_*), whereas a plainly
    # named `x = <constant>` would be parsed as an initial value and drop out of the
    # solve.  So a constant sensitivity (e.g. d(R)/d(eta)=0 for additive error) still
    # comes back as a column of that constant -- no special-casing needed.
    .rvarL <- character(0)
    if (!is.null(.rvar)) {
      .rL1 <- c(vapply(dirs[!(dirs %in% .covDirs)], function(.p) paste0("rx_rvar1_", .p, "=", .toRx(.g1(.rvar, .p))), character(1)),
                vapply(.covDirs, function(.cd) .emitScaled1("rx_rvar1_", "", .cd), character(1)))
      .pr2mod <- !(.P2r$i %in% .covDirs) & !(.P2r$j %in% .covDirs)
      .rL2 <- c(vapply(which(.pr2mod), function(.r) .emitPair("rx_rvar2_", .P2r$i[.r], .P2r$j[.r], function(.a, .b) .toRx(.g2(.rvar, .a, .b)), dirs), character(1)),
                vapply(which(!.pr2mod), function(.r) .emitPair("rx_rvar2_", .P2r$i[.r], .P2r$j[.r], NULL, dirs), character(1)))
      .rvarL <- c(paste0("rx_rvarf_=", .toRx(.rvar)), .rL1, .rL2)
    }
    # Residual-variance sigma sensitivities: for each error parameter the variance R
    # depends on (an error-model theta appearing in rx_r_), emit dR/dsigma
    # (rx_rsig_<n>_), its direction sensitivity d2R/(dsigma ddir) (rx_rsig1_<n>_<dir>)
    # and the sigma-pair 2nd derivative d2R/(dsigma dsigma') (rx_rsig2_<n>_<n'>).  In
    # the (f,R) assembly a sigma is a pseudo-direction with df/dsigma=0 and dR/dsigma
    # from here, so every parameter contracts through the same rho(f,R) machinery.
    .sigTh <- integer(0); .sigL <- character(0)
    if (!is.null(.rvar)) {
      .rvarStr <- as.character(.rvar)
      .erTh <- ui$iniDf$ntheta[!is.na(ui$iniDf$err)]
      .frTh <- suppressWarnings(unique(as.integer(sub("THETA_([0-9]+)_", "\\1",
        regmatches(.rvarStr, gregexpr("THETA_[0-9]+_", .rvarStr))[[1]]))))
      # An ESTIMATED boxCox/yeoJohnson lambda is an error theta that occupies a sigma
      # slot in the parameter accounting, but rx_r_ never mentions it (dR/dlambda = 0 --
      # the transform moves f and the DV, not the variance).  The intersection above
      # therefore dropped it, leaving nsg = 2 sigma directions against ONE emitted
      # rx_rsig_ column, and the pooled kernel threw
      #   copy into submatrix: incompatible matrix dimensions: 11x2 and 11x1
      # which the bare catch(...) in analyticOuterGradDirect swallowed into a silent
      # "grad: fd".  Emit the full sigma set instead: dR/dlambda is identically zero, and
      # an `rx_rsig_<n>_ = 0` still comes back as a real solve column (the rx_<name>_
      # naming rule -- see CLAUDE.md), so the column count matches the direction count and
      # the zero contributes nothing.  Lambda's real gradient still arrives through its
      # theta direction and the -2*jacSum Jacobian term, unchanged.
      .li <- ui$iniDf
      .lamTh <- .li$ntheta[!is.na(.li$ntheta) & !is.na(.li$err) & !.li$fix &
                             .li$err %in% c("boxCox", "yeoJohnson")]
      .sigTh <- sort(union(intersect(.erTh, .frTh), .lamTh))
      for (.n in .sigTh) {
        .sg <- paste0("THETA_", .n, "_"); .dRs <- .Dn(.rvar, .sg)
        .sigL <- c(.sigL, paste0("rx_rsig_", .n, "_=", .toRx(.dRs)),
          vapply(dirs[!(dirs %in% .covDirs)], function(.p) paste0("rx_rsig1_", .n, "_", .p, "=", .toRx(.g1(.dRs, .p))), character(1)),
          vapply(.covDirs, function(.cd) .emitScaled1(paste0("rx_rsig1_", .n, "_"), "", .cd), character(1)),
          vapply(.sigTh[.sigTh >= .n], function(.n2)
            paste0("rx_rsig2_", .n, "_", .n2, "=", .toRx(.Dn(.dRs, paste0("THETA_", .n2, "_")))), character(1)))
      }
    }
    # Transform-of-both-sides parameters (rx_yj_/rx_lambda_/rx_low_/rx_hi_): emit copies
    # so the DV can be transformed in R (y = tbs(DV)) the same way the inner model does.
    .tvarL <- character(0); .hasTrans <- FALSE
    for (.tv in c("yj", "lambda", "low", "hi")) {
      .tval <- tryCatch(get(paste0("rx_", .tv, "_"), .s), error = function(e) NULL)
      if (!is.null(.tval)) .tvarL <- c(.tvarL, paste0("rx_t", .tv, "_=", .toRx(.tval)))
    }
    .hasTrans <- length(.tvarL) > 0L && !identical(as.character(ui$predDf$transform), "untransformed")
    .baseOde <- vapply(.st, function(.x) paste0("d/dt(", .x, ")=", .toRx(get(paste0("rx__d_dt_", .x, "__"), .s))), character(1))
    # State initial conditions (state(0) <- expr).  Emit the base IC AND its
    # per-direction 1st/2nd-order sensitivity-compartment ICs -- otherwise the
    # augmented model starts every state (and sensitivity compartment) at 0,
    # which for a model whose prediction is driven by a parameter-dependent IC
    # (e.g. y(0) <- exp(ta), or a DDE history amplitude) makes f -- and every
    # residual-weighted covariance term -- badly wrong.  The IC is a function of
    # parameters/covariates evaluated at t=0 (before integration), so its
    # direction derivatives are direct partials (no state-sensitivity chain).
    # Compartments whose IC .rxSens already emitted (e.g. the DDE delay-sensitivity
    # augmentation writes the sensitivity-compartment histories/ICs itself); skip
    # those to avoid a duplicate `cmt(0)=` assignment.
    # as.character(): with no ODE states .rxSens returns an empty LIST, not an empty
    # character vector, and strsplit() rejects it ("non-character argument").
    .icDone <- trimws(sub("\\(0\\)=.*$", "", grep("\\(0\\)=", unlist(strsplit(as.character(c(.s1, .s2)), "\n")), value = TRUE)))
    .emitIc <- function(.cmt, .expr) if (identical(.toRx(.expr), "0") || .cmt %in% .icDone) character(0)
      else paste0(.cmt, "(0)=", .toRx(.expr))
    .icL <- character(0)
    for (.x in .st) {
      .ic <- tryCatch(get(paste0("rx_", .x, "_ini_0__"), .s), error = function(e) NULL)
      if (is.null(.ic)) next
      if (!(.x %in% .icDone)) .icL <- c(.icL, paste0(.x, "(0)=", .toRx(.ic)))  # base state IC
      # only model directions have sensitivity compartments (covariate directions are emitted by
      # scaling, so their IC contribution is already in the scaled eta columns -- emitting a
      # phantom rx__sens_..._BY_<covDir>_(0)= for a non-existent compartment would break the model).
      for (.p in .mfDirs) .icL <- c(.icL, .emitIc(paste0("rx__sens_", .x, "_BY_", .p, "__"), .Dn(.ic, .p)))
      .p2m <- which(!(.P2$i %in% .covDirs) & !(.P2$j %in% .covDirs))
      for (.r in .p2m)
        .icL <- c(.icL, .emitIc(paste0("rx__sens_", .x, "_BY_", .P2$i[.r], "_BY_", .P2$j[.r], "__"),
                                .Dn(.Dn(.ic, .P2$i[.r]), .P2$j[.r])))
    }
    # Dosing modifiers (bioavailability f(), lag()/alag(), rate(), dur()) live in the
    # pruned env as rx_<mod>_<state>_ and are NOT part of rx__d_dt_*.  Emit them so the
    # augmented model doses correctly AND so eventSens="jump" (below) fills the analytic
    # dose-parameter ("jump") sensitivities for the rx__sens_* compartments -- otherwise
    # a modeled dosing parameter's sensitivity is silently zero.
    .dosVars <- grep("^rx_(f|lag|alag|rate|dur)_.+_$", ls(envir = .s, all.names = TRUE), value = TRUE)
    .dose <- vapply(.dosVars, function(.v) {
      .m <- regmatches(.v, regexec("^rx_(f|lag|alag|rate|dur)_(.+)_$", .v))[[1]]
      .fun <- if (.m[2L] == "lag") "alag" else .m[2L]     # rxode2 stores lag() as alag()
      paste0(.fun, "(", .m[3L], ")=", .toRx(get(.v, envir = .s)))
    }, character(1))
    # DDE non-constant delay() pre-history: base past(state,tau)<-expr + the
    # per-sensitivity-compartment 1st/2nd-order histories that .rxSens()
    # accumulated as a side effect of the .s1/.s2 calls above (rxode2's
    # .rxDelaySensAugment/.rxDelaySensAugment2).  Spliced after all d/dt +
    # sensitivity-ODE declarations so the referenced state/sens compartments
    # are already defined -- same convention as .rxFinalizeInner (focei.R).
    # NULL (no delay()) for ordinary models.
    .pastLines <- .s$..pastLines
    if (is.null(.pastLines)) .pastLines <- character(0)
    .modTxt <- paste(c(.baseOde, .dose, .s1, .s2, .icL, .pastLines, paste0("rx_predf_=", .toRx(.pred)), .fL1, .fL2, .rvarL, .sigL, .tvarL), collapse = "\n")
    .modTxt <- gsub("ETA\\[([0-9]+)\\]", "ETA_\\1_", .modTxt); .modTxt <- gsub("THETA\\[([0-9]+)\\]", "THETA_\\1_", .modTxt)
    # Optimize common subexpressions (as the inner model does): the augmented model
    # has heavy shared subexpressions across the sensitivity ODEs and the f1/f2
    # prediction chains, so rxOptExpr materially shrinks the per-solve work.
    # rxode2's rxOptExpr() itself optimizes a large model in cost-balanced chunks by
    # default (it disguises the compartment-scoped constructs -- state ICs and
    # f/alag/lag/rate/dur dosing modifiers -- a chunk would otherwise orphan, and
    # falls back to a single whole-model pass if a chunk cannot be optimized), so the
    # ~200-350 line augmented model no longer needs hand-rolled chunking here.  The
    # fit's rxControl(cores=) is passed through so the chunks are optimized in
    # parallel with the same thread setting the solves use (`.cores` also rides along
    # in the result so the batched solves run parallel).
    .cores <- .optExprCores(ui)
    if (isTRUE(rxode2::rxGetControl(ui, "optExpression", TRUE))) {
      .modTxt <- tryCatch(
        rxode2::rxOptExpr(.modTxt, "FOCEi outer gradient model", parallel = .cores),
        error = function(e) .modTxt)
    }
    # Declare the theta/eta inputs AND the model covariates up front (param()) so the
    # solve parameter order is fixed and positional.  Reuse .uiGetThetaEtaParams -- the
    # SAME theta/eta/covariate ordering the inner sensitivity model uses -- so the
    # augmented outer/analytic model's parameter order matches the inner problem
    # exactly (thetas in ntheta order, etas in neta order, then covariates).
    .param <- .uiGetThetaEtaParams(ui, TRUE)          # params(THETA[1], .., ETA[1], .., covs)
    .param <- gsub("ETA\\[([0-9]+)\\]", "ETA_\\1_", gsub("THETA\\[([0-9]+)\\]", "THETA_\\1_", .param))
    # Same compartment/endpoint prologue and epilogue the inner, predOnly and predNoLhs
    # models get from .toRx (which pastes toRxParam + body + toRxDvidCmt).  This model
    # builds its own text and calls .nlmixr2estRxode2 directly, so without this it was
    # the only peer with no cmt() pins and an EMPTY dvid, i.e. a different endpoint
    # contract from every model it shares the solve pool with.  Aligning it is right on
    # its own terms and is a prerequisite for ever pooling a multi-endpoint model.
    #
    # It is NOT sufficient, and multi-endpoint models still do not pool (see the
    # outerPoolOk gate in focei.R).  The endpoint cmt() is only legal AFTER the d/dt
    # block, so it is numbered after that model's generated sensitivity states -- 5/6
    # for the 4-state inner model, 63/64 for this 62-state one.  CMT reaches a model as
    # that solve-compartment index (measured: par_ptr last slot, 5 vs 63), so peers whose
    # state counts differ cannot share one translated event table and both resolve their
    # `CMT ==` endpoint switch.  Single-endpoint models have no switch, which is the only
    # reason they pool safely today.
    #
    # Verified neutral where it does apply: single-endpoint theo_sd fast=TRUE still pools
    # (pool = outer), still takes the analytic gradient, objf 133.654382798 against
    # fast=FALSE's 133.654382066.
    .cmtPre <- ui$foceiCmtPreModel
    .interp <- ui$interpLinesStr
    if (!is.null(.interp) && .interp != "") .cmtPre <- paste0(.cmtPre, "\n", .interp)
    .modTxt <- paste(c(.param, .cmtPre, .modTxt, .foceiToCmtLinesAndDvid(ui)), collapse = "\n")
    # no splitBolus() in the augmented model -- it translates the already-split
    # dataSav, so declaring it would split the doses twice (.foceiPreProcessData)
    # eventSens="jump" attaches rxode2's analytic event/dosing-parameter sensitivities
    # (forward variational jumps at dose times) for the sensitivity compartments.  `cols`
    # precomputes solve-output column names/index maps; `cores` carries the fit's rxControl thread
    # count so the batched solves run parallel; `key` seeds the per-fit event-table reuse cache.
    list(augMod = .nlmixr2estRxode2(.modTxt, "rxOuter", eventSens = "jump"),
         dirs = dirs, ndir = length(dirs), fDirs = .fDirs,
         st = .st, P2 = .P2, P2r = .P2r, hasRvar = !is.null(.rvar), sigTh = .sigTh, hasTrans = .hasTrans,
         cols = .foceiAnalyticCols(dirs, .fDirs, .P2, .P2r, .sigTh), cores = .cores, key = .key)
  }, error = function(e) NULL)
  if (!is.null(.key) && !is.null(.res)) {
    if (length(ls(.foceiAnalyticAugCache, all.names = TRUE)) >= 64L)    # bound retained compiled models
      rm(list = ls(.foceiAnalyticAugCache, all.names = TRUE), envir = .foceiAnalyticAugCache)
    assign(.key, .res, envir = .foceiAnalyticAugCache)
  }
  .res
}

#' Solve the direction-set 2nd-order model for one subject and recover the
#' 3rd-order tensor `Ath` by Shi (2021) central differences of the analytic
#' 2nd-order sensitivities `A` (C++ `shi21CentralWrap`, perturbing one full-param
#' coordinate per direction), then symmetrize.  Returns `list(f, a, A, Ath)` or
#' `NULL` on failure.  (`rxode2::rxExpandSens3_` would give `Ath` analytically but
#' at O(ndir^3) augmented-model compile cost; the Shi FD keeps it at O(ndir^2).)
#' @noRd
.foceiAnalyticSolveSubjectFD3 <- function(aug, params, ev, times, .fdEps = 7e-7, tol = 1e-10, withR = FALSE) {
  dirs <- aug$dirs; nd <- length(dirs)
  # base solve (f, a, A) via the shared helper; this tier adds the 3rd-order Ath by
  # Shi-differencing A.
  E0 <- .foceiAnalyticSolveFA(aug, params, ev, times, tol = tol); if (is.null(E0)) return(NULL)
  nobs <- length(E0$f)
  # 3rd-order sensitivities Ath = d3f/(deta ddir ddir') and (f,R cov) AthR the same for R,
  # by Shi-differencing the analytic 2nd-order A / AR.  The cov assembly only ever reads
  # these with the FIRST index an ETA (Tn[l,s,t], d2HtDD both have l,m in the eta set), so
  # only the eta directions need to be Shi-differenced -- the theta/sigma "slices" are never
  # used.  This is O(neta) solves instead of O(ndir): the returned tensors are
  # [obs, neta, ndir, ndir] with the eta (differencing) axis first.  A is symmetric in its
  # last two axes, so d(A)/deta is too -- no cross-axis symmetrization is needed.
  # withR=TRUE (the (f,R) cov) also differences AR; the symbolic add/prop cov (R=R(f)) skips it.
  .etaDir <- which(grepl("^ETA_[0-9]+_$", dirs)); neta <- length(.etaDir)
  if (neta == 0L) return(NULL)
  .hasR <- isTRUE(withR) && !is.null(E0$AR); nA <- nobs * nd * nd
  f0 <- if (.hasR) c(as.vector(E0$A), as.vector(E0$AR)) else as.vector(E0$A)
  # shi21CentralWrap differences the closure at a perturbed full-param vector; it
  # strips names, so re-attach names(params) before mapping back through the solve.
  Aflat <- function(.tt) { E <- .foceiAnalyticSolveFA(aug, setNames(.tt, names(params)), ev, times, tol = tol)
    if (is.null(E) || !all(is.finite(E$A))) return(NULL)
    if (.hasR) { if (!all(is.finite(E$AR))) return(NULL); c(as.vector(E$A), as.vector(E$AR)) } else as.vector(E$A) }
  Ath <- array(0, c(nobs, neta, nd, nd)); AthR <- if (.hasR) array(0, c(nobs, neta, nd, nd)) else NULL
  for (li in seq_len(neta)) {
    idx <- match(dirs[.etaDir[li]], names(params))       # this eta's coordinate in params
    if (is.na(idx)) return(NULL)
    sc <- shi21CentralWrap(Aflat, params, f0, idx, .fdEps)  # C++ shi21Central
    if (is.null(sc$gr) || !all(is.finite(sc$gr))) return(NULL)
    Ath[, li, , ] <- array(sc$gr[seq_len(nA)], c(nobs, nd, nd))
    if (.hasR) AthR[, li, , ] <- array(sc$gr[nA + seq_len(nA)], c(nobs, nd, nd))
  }
  .out <- c(list(f = E0$f, a = E0$a, A = E0$A, Ath = Ath),
            if (.hasR) list(R = E0$R, aR = E0$aR, AR = E0$AR, AthR = AthR,
                            Rsig = E0$Rsig, RsigDir = E0$RsigDir, Rsig2 = E0$Rsig2) else NULL,
            if (!is.null(E0$trans)) list(trans = E0$trans) else NULL)
  if (!all(is.finite(.out$f)) || !all(is.finite(.out$a)) || !all(is.finite(.out$A)) || !all(is.finite(.out$Ath))) return(NULL)
  .out
}

#' Route-A (default on; FOCEI_NO_RSIG opts out) sigma reconstruction: the (f,R) cov shares the gradient's `dirs` model
#' (no sigma directions) and rebuilds the sigma tensor slots from the model's own residual-sigma
#' outputs (Rsig=dR/dsigma, RsigDir=d2R/(dsigma ddir), Rsig2=d2R/(dsigma dsigma')) plus their
#' eta-derivatives (RsigDirEta, Rsig2Eta from the FD3).  A sigma has df/dsigma=0, so the prediction
#' slots (a/A/Ath) are zero; the variance slots (aR/AR/AthR) carry the rsig values.  The result is
#' the SAME ndirCov tensor the sigma-as-direction model produced -- but from a model shared with the
#' gradient.
#' @noRd
.foceiAnalyticExpandSigma <- function(E, nsig, neta, RsigDirEta, Rsig2Eta, sigSel = NULL) {
  # `sigSel` selects WHICH rsig columns to append: exactly the sigma directions the
  # caller dropped.  The model emits rx_rsig_ for every sigTh (lambda included), which
  # is a SUPERSET of what Route A drops (.erN excludes boxCox/yeoJohnson) -- appending
  # all of them widens E past the caller's ndirCov buffers and `aB[rows, ] <- E$a`
  # throws.  NULL keeps the old all-columns behaviour.
  if (!is.null(sigSel)) {
    E$Rsig <- E$Rsig[, sigSel, drop = FALSE]
    if (!is.null(E$RsigDir)) E$RsigDir <- E$RsigDir[, , sigSel, drop = FALSE]
    if (!is.null(E$Rsig2)) E$Rsig2 <- E$Rsig2[, sigSel, sigSel, drop = FALSE]
    if (!is.null(RsigDirEta)) RsigDirEta <- RsigDirEta[, , , sigSel, drop = FALSE]
    if (!is.null(Rsig2Eta)) Rsig2Eta <- Rsig2Eta[, , sigSel, sigSel, drop = FALSE]
    nsig <- length(sigSel)
  }
  nd <- ncol(E$a); no <- nrow(E$a); ndc <- nd + nsig; sg <- nd + seq_len(nsig)
  a <- cbind(E$a, matrix(0, no, nsig))
  A <- array(0, c(no, ndc, ndc)); A[, seq_len(nd), seq_len(nd)] <- E$A
  .hasAth <- !is.null(E$Ath)                              # E0 (frozen-R0 solve) has no Ath
  Ath <- if (.hasAth) { .Y <- array(0, c(no, neta, ndc, ndc)); .Y[, , seq_len(nd), seq_len(nd)] <- E$Ath; .Y } else NULL
  aR <- cbind(E$aR, E$Rsig)
  AR <- array(0, c(no, ndc, ndc)); AR[, seq_len(nd), seq_len(nd)] <- E$AR
  .hasAthR <- !is.null(E$AthR) && !is.null(RsigDirEta)     # withR: 3rd-order sigma slots too
  AthR <- if (.hasAthR) { .X <- array(0, c(no, neta, ndc, ndc)); .X[, , seq_len(nd), seq_len(nd)] <- E$AthR; .X } else NULL
  for (k in seq_len(nsig)) {
    AR[, seq_len(nd), sg[k]] <- E$RsigDir[, , k]; AR[, sg[k], seq_len(nd)] <- E$RsigDir[, , k]
    if (.hasAthR) { AthR[, , seq_len(nd), sg[k]] <- RsigDirEta[, , , k]; AthR[, , sg[k], seq_len(nd)] <- RsigDirEta[, , , k] }
    for (l in seq_len(nsig)) { AR[, sg[k], sg[l]] <- E$Rsig2[, k, l]; if (.hasAthR) AthR[, , sg[k], sg[l]] <- Rsig2Eta[, , k, l] }
  }
  E$a <- a; E$A <- A; if (.hasAth) E$Ath <- Ath; E$aR <- aR; E$AR <- AR; if (.hasAthR) E$AthR <- AthR; E
}

#' Batched analogue of [.foceiAnalyticSolveSubjectFD3]: recover the 3rd-order tensor `Ath` (and
#' `AthR` when `withR`) for ALL subjects at once by central-differencing the analytic 2nd-order
#' sensitivities `A` w.r.t. each ETA coordinate via BATCHED population solves ([.foceiAnalyticSolveAll],
#' 1 base + 2*neta perturbed) instead of the per-subject Shi (~O(nsub*neta) solves).  RICHARDSON-
#' extrapolated over steps (h, h/2) to 4th order, with the perturbed solves at a tighter tol, so the
#' batched Ath reproduces the per-subject adaptive-Shi Ath (the FOCEI-vs-FOCE analytic-R agreement,
#' not just the SEs).  Returns the per-subject E-list with `Ath` attached, or NULL (-> per-subject).
#' @noRd
.foceiAnalyticSolveAllFD3 <- function(am, thv, ebes, ids, data, obsTimes, tol = 1e-10,
                                      fdEps = 1e-3, withR = FALSE, sigSel = NULL) {
  dirs <- am$dirs; nd <- length(dirs); neta <- ncol(ebes)
  E0 <- .foceiAnalyticSolveAll(am, thv, ebes, ids, data, obsTimes, tol)
  if (is.null(E0)) return(NULL)
  nsub <- length(E0)
  Ath  <- lapply(E0, function(E) array(0, c(nrow(E$a), neta, nd, nd)))
  AthR <- if (withR) lapply(E0, function(E) array(0, c(nrow(E$a), neta, nd, nd))) else NULL
  # Route A (default ON; FOCEI_NO_RSIG=1 opts out): am is the gradient's `dirs` model (no sigma
  # directions); rebuild the sigma tensor slots from the rsig outputs + their eta-derivatives
  # (differenced here alongside A/AR).  Guarded on Rsig actually being present in the solve.
  # expand only when the caller actually dropped sigma directions (sigSel); expanding on
  # Rsig presence alone widens E past the caller's buffers when the two disagree
  .rsig <- !nzchar(Sys.getenv("FOCEI_NO_RSIG")) && !is.null(E0[[1L]]$Rsig) &&
    length(E0[[1L]]$Rsig) > 0L && (is.null(sigSel) || length(sigSel) > 0L)
  nsig <- if (.rsig) ncol(E0[[1L]]$Rsig) else 0L
  RsigDirEta <- if (.rsig && withR) lapply(E0, function(E) array(0, c(nrow(E$a), neta, nd, nsig))) else NULL
  Rsig2Eta   <- if (.rsig && withR) lapply(E0, function(E) array(0, c(nrow(E$a), neta, nsig, nsig))) else NULL
  # differencing A carries the ODE solve's ~tol error floor; solve the PERTURBED models tighter
  # than the base so the 3rd-order Ath stays as clean as the per-subject adaptive Shi.
  .ptol <- min(tol, 1e-12)
  .cdiff <- function(li, h) {
    ep <- ebes; ep[, li] <- ebes[, li] + h; em <- ebes; em[, li] <- ebes[, li] - h
    Ep <- .foceiAnalyticSolveAll(am, thv, ep, ids, data, obsTimes, .ptol)
    Em <- .foceiAnalyticSolveAll(am, thv, em, ids, data, obsTimes, .ptol)
    if (is.null(Ep) || is.null(Em) || length(Ep) != nsub || length(Em) != nsub) return(NULL)
    list(Ep = Ep, Em = Em, h = h)
  }
  for (li in seq_len(neta)) {                             # ETA_li coordinate == ebes column li
    h <- fdEps * max(abs(ebes[, li]), 1)
    s1 <- .cdiff(li, h); s2 <- .cdiff(li, h / 2)
    if (is.null(s1) || is.null(s2)) return(NULL)
    for (i in seq_len(nsub)) {
      d1 <- (s1$Ep[[i]]$A - s1$Em[[i]]$A) / (2 * s1$h); d2 <- (s2$Ep[[i]]$A - s2$Em[[i]]$A) / (2 * s2$h)
      Ath[[i]][, li, , ] <- (4 * d2 - d1) / 3
      if (withR) {
        e1 <- (s1$Ep[[i]]$AR - s1$Em[[i]]$AR) / (2 * s1$h); e2 <- (s2$Ep[[i]]$AR - s2$Em[[i]]$AR) / (2 * s2$h)
        AthR[[i]][, li, , ] <- (4 * e2 - e1) / 3
        if (.rsig) {
          g1 <- (s1$Ep[[i]]$RsigDir - s1$Em[[i]]$RsigDir) / (2 * s1$h); g2 <- (s2$Ep[[i]]$RsigDir - s2$Em[[i]]$RsigDir) / (2 * s2$h)
          RsigDirEta[[i]][, li, , ] <- (4 * g2 - g1) / 3
          k1 <- (s1$Ep[[i]]$Rsig2 - s1$Em[[i]]$Rsig2) / (2 * s1$h); k2 <- (s2$Ep[[i]]$Rsig2 - s2$Em[[i]]$Rsig2) / (2 * s2$h)
          Rsig2Eta[[i]][, li, , ] <- (4 * k2 - k1) / 3
        }
      }
    }
  }
  for (i in seq_len(nsub)) {
    E0[[i]]$Ath <- Ath[[i]]
    if (withR) E0[[i]]$AthR <- AthR[[i]]
    if (.rsig) E0[[i]] <- .foceiAnalyticExpandSigma(E0[[i]], nsig, neta,
                            if (is.null(RsigDirEta)) NULL else RsigDirEta[[i]],
                            if (is.null(Rsig2Eta)) NULL else Rsig2Eta[[i]], sigSel)
    if (!all(is.finite(E0[[i]]$A)) || !all(is.finite(E0[[i]]$Ath))) return(NULL)
  }
  E0
}

#' Per-subject observed-information R in the general (f,R) form: the prediction f and
#' the variance R are independent solved quantities with sensitivities a/A/Ath and
#' aR/AR/AthR (Ath/AthR by Shi-FD, as the add/prop cov FD-differences A -> Ath).  Every
#' non-Omega parameter (structural theta AND residual sigma) is a DIRECTION; a sigma
#' direction has a=A=Ath=0 (f is sigma-independent) so only its variance sensitivities
#' contribute, and no separate sigma machinery is needed.  The rho(f,R,y) partials up to
#' 3rd order are model-independent closed forms.  `dirP` maps each non-Omega param to its
#' direction; Omega params follow, using `omd`.  Reduces to `.foceiAnalyticSubjectR` when
#' R=R(f).  FOCE (interaction=0) delegates to the frozen-R0 variant.
#' @noRd
.foceiAnalyticSubjectRFR <- function(E, ehat, Om, neta, ndirP, dirP, omd,
                                     ndir = neta, Oi = solve(Om), interaction = 1L,
                                     E0 = NULL, foceType = 0L) {
  if (identical(as.integer(interaction), 0L))
    return(.foceiAnalyticSubjectRfoceFR(E, ehat, Om, neta, ndirP, dirP, omd,
                                        ndir = ndir, Oi = Oi, E0 = E0, foceType = foceType))
  tr <- function(M) sum(diag(M))
  f <- E$f; y <- E$y; R <- E$R; a <- E$a; A <- E$A; Ath <- E$Ath
  aR <- E$aR; AR <- E$AR; AthR <- E$AthR
  res <- y - f
  rf <- -res / R; rR <- 0.5 * (1 / R - res^2 / R^2)
  rff <- 1 / R; rfR <- res / R^2; rRR <- 0.5 * (-1 / R^2 + 2 * res^2 / R^3)
  rffR <- -1 / R^2; rfRR <- -2 * res / R^3; rRRR <- 0.5 * (2 / R^3 - 6 * res^2 / R^4)  # rfff = 0
  iR <- 1 / R; iR2 <- iR^2; iR3 <- iR^3; iR4 <- iR^4
  nom <- omd$nom; np <- ndirP + nom
  ei <- seq_len(neta); di <- seq_len(ndir)
  # (f,R) 2nd total derivative of the density w.r.t. two directions da, db
  Gdd <- function(da, db) sum(rff * a[, da] * a[, db] + rfR * (a[, da] * aR[, db] + aR[, da] * a[, db]) +
                                rRR * aR[, da] * aR[, db] + rf * A[, da, db] + rR * AR[, da, db])
  H <- Oi; for (l in ei) for (m in ei) H[l, m] <- H[l, m] + Gdd(l, m)
  HiM <- solve(H)
  N <- matrix(0, neta, ndir); for (l in ei) for (d in di) N[l, d] <- Gdd(l, d)
  # exact 3rd total derivative Tn[l,s,t] = d2(Phi_l)/ddir_s ddir_t, Phi_l = rf a_l + rR aR_l
  Tn <- array(0, c(neta, ndir, ndir))
  for (l in ei) for (s in di) for (t in di) {
    us  <- rff * a[, s] + rfR * aR[, s]; ut <- rff * a[, t] + rfR * aR[, t]
    ust <- rffR * (a[, s] * aR[, t] + aR[, s] * a[, t]) + rfRR * aR[, s] * aR[, t] + rff * A[, s, t] + rfR * AR[, s, t]
    ws  <- rfR * a[, s] + rRR * aR[, s]; wt <- rfR * a[, t] + rRR * aR[, t]
    wst <- rffR * a[, s] * a[, t] + rfRR * (a[, s] * aR[, t] + aR[, s] * a[, t]) + rRRR * aR[, s] * aR[, t] +
      rfR * A[, s, t] + rRR * AR[, s, t]
    Tn[l, s, t] <- sum(ust * a[, l] + us * A[, l, t] + ut * A[, l, s] + rf * Ath[, l, s, t] +
                         wst * aR[, l] + ws * AR[, l, t] + wt * AR[, l, s] + rR * AthR[, l, s, t])
  }
  # Laplace determinant Ht = Oi + sum(a a / R + 0.5 aR aR / R^2) and its 1st/2nd
  # direction derivatives (the interaction-free E[rho_fR]=0 split into f- and R-quadratics)
  Ht <- Oi; for (l in ei) for (m in ei) Ht[l, m] <- Ht[l, m] + sum(a[, l] * a[, m] * iR + 0.5 * aR[, l] * aR[, m] * iR2)
  Hti <- solve(Ht)
  dHtD <- lapply(di, function(s) { D <- matrix(0, neta, neta); for (l in ei) for (m in ei)
    D[l, m] <- sum((A[, l, s] * a[, m] + a[, l] * A[, m, s]) * iR - a[, l] * a[, m] * aR[, s] * iR2 +
                     0.5 * (AR[, l, s] * aR[, m] + aR[, l] * AR[, m, s]) * iR2 - aR[, l] * aR[, m] * aR[, s] * iR3); D })
  d2HtDD <- lapply(di, function(s) lapply(di, function(t) { D <- matrix(0, neta, neta); for (l in ei) for (m in ei)
    D[l, m] <- sum(
      (Ath[, l, s, t] * a[, m] + A[, l, s] * A[, m, t] + A[, l, t] * A[, m, s] + a[, l] * Ath[, m, s, t]) * iR -
        (A[, l, s] * a[, m] + a[, l] * A[, m, s]) * aR[, t] * iR2 -
        (A[, l, t] * a[, m] + a[, l] * A[, m, t]) * aR[, s] * iR2 - a[, l] * a[, m] * AR[, s, t] * iR2 +
        2 * a[, l] * a[, m] * aR[, s] * aR[, t] * iR3 +
        0.5 * (AthR[, l, s, t] * aR[, m] + AR[, l, s] * AR[, m, t] + AR[, l, t] * AR[, m, s] + aR[, l] * AthR[, m, s, t]) * iR2 -
        (AR[, l, s] * aR[, m] + aR[, l] * AR[, m, s]) * aR[, t] * iR3 -
        (AR[, l, t] * aR[, m] + aR[, l] * AR[, m, t]) * aR[, s] * iR3 - aR[, l] * aR[, m] * AR[, s, t] * iR3 +
        3 * aR[, l] * aR[, m] * aR[, s] * aR[, t] * iR4); D }))
  Cen <- vapply(ei, function(l) 0.5 * tr(Hti %*% dHtD[[l]]), numeric(1))
  Cee <- matrix(0, neta, neta); for (s in ei) for (t in ei)
    Cee[s, t] <- 0.5 * (tr(Hti %*% d2HtDD[[s]][[t]]) - tr(Hti %*% dHtD[[s]] %*% Hti %*% dHtD[[t]]))
  typ <- function(p) if (p <= ndirP) "dir" else "om"
  dOf <- function(p) dirP[p]; omc <- function(p) p - ndirP
  Mcol <- function(p) if (typ(p) == "dir") N[, dOf(p)] else as.numeric(omd$dOi[[omc(p)]] %*% ehat)
  dHt_p <- function(p) if (typ(p) == "dir") dHtD[[dOf(p)]] else omd$dOi[[omc(p)]]
  d2HtEtaP <- function(p, l) if (typ(p) == "dir") d2HtDD[[dOf(p)]][[l]] else matrix(0, neta, neta)
  d2Ht_pp <- function(aa, bb) { ta <- typ(aa); tb <- typ(bb)
    if (ta == "dir" && tb == "dir") return(d2HtDD[[dOf(aa)]][[dOf(bb)]])
    if (ta == "om" && tb == "om") return(omd$d2Oi[[omc(aa)]][[omc(bb)]])
    matrix(0, neta, neta) }
  Smat <- function(p) { if (typ(p) == "om") return(omd$dOi[[omc(p)]])
    M <- matrix(0, neta, ndir); for (l in ei) for (s in di) M[l, s] <- Tn[l, dOf(p), s]; M }
  Svec <- function(aa, bb) { ta <- typ(aa); tb <- typ(bb)
    if (ta == "dir" && tb == "dir") return(Tn[, dOf(aa), dOf(bb)])
    if (ta == "om" && tb == "om") return(as.numeric(omd$d2Oi[[omc(aa)]][[omc(bb)]] %*% ehat))
    rep(0, neta) }
  d2Phi <- function(aa, bb) { ta <- typ(aa); tb <- typ(bb)
    if (ta == "dir" && tb == "dir") return(Gdd(dOf(aa), dOf(bb)))
    if (ta == "om" && tb == "om") return(0.5 * as.numeric(t(ehat) %*% omd$d2Oi[[omc(aa)]][[omc(bb)]] %*% ehat) + 0.5 * omd$d2LD[omc(aa), omc(bb)])
    0 }
  .Mcols <- lapply(1:np, Mcol)
  etaP <- matrix(vapply(1:np, function(p) as.numeric(-HiM %*% .Mcols[[p]]), numeric(neta)), nrow = neta)
  eta2 <- function(aa, bb) { b <- Svec(aa, bb) + Smat(aa)[, ei, drop = FALSE] %*% etaP[, bb] + Smat(bb)[, ei, drop = FALSE] %*% etaP[, aa]
    for (l in ei) b[l] <- b[l] + as.numeric(t(etaP[, aa]) %*% Tn[l, ei, ei] %*% etaP[, bb]); as.numeric(-HiM %*% b) }
  Cpe <- function(p, l) 0.5 * (tr(Hti %*% d2HtEtaP(p, l)) - tr(Hti %*% dHt_p(p) %*% Hti %*% dHtD[[l]]))
  Cpp <- function(aa, bb) 0.5 * (tr(Hti %*% d2Ht_pp(aa, bb)) - tr(Hti %*% dHt_p(aa) %*% Hti %*% dHt_p(bb)))
  .CpeRow <- lapply(1:np, function(p) vapply(ei, function(l) Cpe(p, l), numeric(1)))
  R <- matrix(0, np, np)
  for (aa in 1:np) for (bb in aa:np) {
    dat <- d2Phi(aa, bb) - as.numeric(t(.Mcols[[aa]]) %*% HiM %*% .Mcols[[bb]])
    ld <- Cpp(aa, bb) + sum(.CpeRow[[aa]] * etaP[, bb]) + sum(.CpeRow[[bb]] * etaP[, aa]) +
          as.numeric(t(etaP[, aa]) %*% Cee %*% etaP[, bb]) + sum(Cen * eta2(aa, bb))
    R[aa, bb] <- R[bb, aa] <- dat + ld
  }
  R
}

#' (f,R) FOCE per-subject observed-information R (interaction=0).  The inner problem is
#' interaction-free (S_FOCE = sum(q0 a) + Omega^-1 eta, q0=-(y-f)/R0, q1=1/R0), so the
#' EBE does not stationarize the full Laplace objective and the general non-envelope form
#' R_ab = F_ab + F_aeta eta_b + F_beta eta_a + eta_a' F_etaeta eta_b + F_eta eta_ab is used.
#' R0 is frozen at the eta=0 population variance (nonmem, from E0) or the live conditional
#' variance (foce+, from E); its theta-chain enters the parameter columns via aRc/ARc
#' (dR0/ddir, d2R0/ddir2 from E0), with the eta-block frozen (aRe=0 for nonmem).  Ath is
#' the eta-hat prediction 3rd-order tensor [obs, neta, ndir, ndir].  Sigmas are directions
#' (a=A=Ath=0).  Reduces to `.foceiAnalyticSubjectRfoce` on add/prop.
#' @noRd
.foceiAnalyticSubjectRfoceFR <- function(E, ehat, Om, neta, ndirP, dirP, omd,
                                         ndir = neta, Oi = solve(Om), E0 = NULL, foceType = 0L) {
  tr <- function(M) sum(diag(M))
  f <- E$f; y <- E$y; a <- E$a; A <- E$A; Ath <- E$Ath
  ei <- seq_len(neta); di <- seq_len(ndir); nobs <- length(f)
  .fp <- identical(as.integer(foceType), 1L) || is.null(E0)
  # frozen variance R0 and its sensitivities.  aRe drives the eta-block (0 for nonmem --
  # R0 is eta-independent -- live E$aR for foce+); aRc drives the parameter columns
  # (E0's full dR0/ddir for nonmem, live E$aR for foce+), so a mu-referenced theta (whose
  # direction is an eta) keeps its frozen-R0 theta sensitivity.
  if (.fp) { R0 <- E$R; aRe <- E$aR; aRc <- E$aR; ARc <- E$AR }
  else { R0 <- E0$R; aRe <- matrix(0, nobs, ndir); aRc <- E0$aR; ARc <- E0$AR }
  res <- y - f
  rf <- -res / R0; rR <- 0.5 * (1 / R0 - res^2 / R0^2)
  rff <- 1 / R0; rfR <- res / R0^2; rRR <- 0.5 * (-1 / R0^2 + 2 * res^2 / R0^3)
  q0 <- rf; q1 <- 1 / R0; iR <- 1 / R0; iR2 <- iR^2; iR3 <- iR^3
  nom <- omd$nom; np <- ndirP + nom
  ae <- a[, ei, drop = FALSE]
  isD <- function(p) p <= ndirP; dOf <- function(p) dirP[p]; omc <- function(p) p - ndirP
  # ---- Phi (data) tensors: H=Phi_etaeta, gPhi=Phi_eta (aRe eta-block) ----
  gPhi <- as.numeric(Oi %*% ehat); for (l in ei) gPhi[l] <- gPhi[l] + sum(rf * a[, l] + rR * aRe[, l])
  H <- Oi; for (l in ei) for (m in ei)
    H[l, m] <- H[l, m] + sum(rff * a[, l] * a[, m] + rfR * (a[, l] * aRe[, m] + aRe[, l] * a[, m]) +
                              rRR * aRe[, l] * aRe[, m] + rf * A[, l, m] + rR * E_ARelm(E, l, m, .fp))
  # ---- FOCE inner (EBE) tensors: interaction-free q-based Hf/Nf/Tnf ----
  Hf <- Oi; Nf <- matrix(0, neta, ndir)
  for (l in ei) { for (m in ei) Hf[l, m] <- Hf[l, m] + sum(q1 * a[, l] * a[, m] + q0 * A[, l, m])
    for (d in di) Nf[l, d] <- sum(q1 * a[, l] * a[, d] + q0 * A[, l, d]) }
  HfInv <- solve(Hf)
  Tnf <- array(0, c(neta, ndir, ndir)); for (l in ei) for (s in di) for (t in di)
    Tnf[l, s, t] <- sum(q1 * (A[, l, s] * a[, t] + A[, l, t] * a[, s] + A[, s, t] * a[, l]) + q0 * Ath[, l, s, t])
  # ---- determinant Ht = Oi + sum(a a / R0) (interaction-free) + its derivatives ----
  # dHtDir/d2HtDir are parameterized by the R0-sensitivity of each direction: the eta-block
  # uses aRe (frozen, 0 for nonmem), the parameter columns use aRc (E0's dR0/ddir), so a
  # mu-referenced theta keeps its frozen-R0 chain.  ARv is d2R0 for the pair (0 unless both
  # indices carry a live R0 dependence).
  Ht <- Oi; for (l in ei) for (m in ei) Ht[l, m] <- Ht[l, m] + sum(a[, l] * a[, m] * iR); Hti <- solve(Ht)
  ARblk <- function(s, t) if (.fp) E$AR[, s, t] else rep(0, nobs)   # d2R0/(dir_s dir_t), eta-block
  dHtDir <- function(s, aRvS) { D <- matrix(0, neta, neta); for (l in ei) for (m in ei)
    D[l, m] <- sum((A[, l, s] * a[, m] + a[, l] * A[, m, s]) * iR - a[, l] * a[, m] * aRvS[, s] * iR2); D }
  d2HtDir <- function(s, t, aRvS, aRvT, ARv) { D <- matrix(0, neta, neta); for (l in ei) for (m in ei)
    D[l, m] <- sum((Ath[, l, s, t] * a[, m] + A[, l, s] * A[, m, t] + A[, l, t] * A[, m, s] + a[, l] * Ath[, m, s, t]) * iR -
      (A[, l, s] * a[, m] + a[, l] * A[, m, s]) * aRvT[, t] * iR2 -
      (A[, l, t] * a[, m] + a[, l] * A[, m, t]) * aRvS[, s] * iR2 - a[, l] * a[, m] * ARv * iR2 +
      2 * a[, l] * a[, m] * aRvS[, s] * aRvT[, t] * iR3); D }
  dHtE <- lapply(ei, function(l) dHtDir(l, aRe))                    # eta-block dHt/deta_l
  d2HtEE <- lapply(ei, function(s) lapply(ei, function(t) d2HtDir(s, t, aRe, aRe, ARblk(s, t))))
  Cen <- vapply(ei, function(l) 0.5 * tr(Hti %*% dHtE[[l]]), numeric(1))
  Cee <- matrix(0, neta, neta); for (s in ei) for (t in ei)
    Cee[s, t] <- 0.5 * (tr(Hti %*% d2HtEE[[s]][[t]]) - tr(Hti %*% dHtE[[s]] %*% Hti %*% dHtE[[t]]))
  dHtP <- function(p) if (isD(p)) dHtDir(dOf(p), aRc) else omd$dOi[[omc(p)]]  # parameter dHt/dp
  # ---- R0 theta-chains for the parameter accessors (aRc/ARc); sigma is a direction ----
  chQ <- function(d) (res / R0^2) * aRc[, d]                       # d(q0)/dtheta via R0
  # Phi_(eta,p) and S_p share the (res/R0^2) aRc chain (q0 = Phi_f for FOCE)
  McolData <- function(p) { if (!isD(p)) return(as.numeric(omd$dOi[[omc(p)]] %*% ehat))
    d <- dOf(p); Nf[, d] + as.numeric(crossprod(ae, chQ(d))) }
  McolEBE <- McolData
  d2Phi <- function(aa, bb) { if (!isD(aa) && !isD(bb))
      return(0.5 * as.numeric(t(ehat) %*% omd$d2Oi[[omc(aa)]][[omc(bb)]] %*% ehat) + 0.5 * omd$d2LD[omc(aa), omc(bb)])
    if (!isD(aa) || !isD(bb)) return(0)
    da <- dOf(aa); db <- dOf(bb)
    sum(rff * a[, da] * a[, db] + rf * A[, da, db] + rfR * (a[, da] * aRc[, db] + aRc[, da] * a[, db]) +
          rRR * aRc[, da] * aRc[, db] + rR * ARc[, da, db]) }
  # S_(p,eta) row (SmatEBE): Tnf plus the theta chain -aRc/R0^2 a a + (res/R0^2) aRc A
  SmatEBE <- function(p) { if (!isD(p)) return(omd$dOi[[omc(p)]])
    d <- dOf(p); M <- matrix(0, neta, ndir); for (l in ei) for (s in di)
      M[l, s] <- Tnf[l, d, s] + sum(-aRc[, d] * iR2 * a[, s] * a[, l] + (res * iR2) * aRc[, d] * A[, l, s]); M }
  # S_(p,p') vector (SvecEBE): Tnf + the combined 2nd-order R0 chain (R0'A0 cancels)
  SvecEBE <- function(aa, bb) { ta <- isD(aa); tb <- isD(bb)
    if (!ta && !tb) return(as.numeric(omd$d2Oi[[omc(aa)]][[omc(bb)]] %*% ehat))
    if (!ta || !tb) return(rep(0, neta))
    da <- dOf(aa); db <- dOf(bb); v <- Tnf[, da, db]
    w <- -(a[, da] * aRc[, db] + aRc[, da] * a[, db]) * iR2 + (res * iR2) * ARc[, da, db] - 2 * res * aRc[, da] * aRc[, db] * iR3
    for (l in ei) v[l] <- v[l] + sum(w * a[, l] + (res * iR2) * (aRc[, da] * A[, l, db] + aRc[, db] * A[, l, da])); v }
  # determinant d2Ht/(deta_l dp) uses aRc for the theta-direction, aRe for the eta; the mixed
  # d2R0/(deta dtheta) is 0 for nonmem (R0 frozen w.r.t. eta) and E$AR for foce+.
  d2HtEtaP <- function(p, l) { if (!isD(p)) return(matrix(0, neta, neta))
    d <- dOf(p); d2HtDir(d, l, aRc, aRe, if (.fp) E$AR[, d, l] else rep(0, nobs)) }
  d2Ht_pp <- function(aa, bb) { ta <- isD(aa); tb <- isD(bb)
    if (ta && tb) return(d2HtDir(dOf(aa), dOf(bb), aRc, aRc, ARc[, dOf(aa), dOf(bb)]))
    if (!ta && !tb) return(omd$d2Oi[[omc(aa)]][[omc(bb)]])
    matrix(0, neta, neta) }
  Cpe <- function(p, l) 0.5 * (tr(Hti %*% d2HtEtaP(p, l)) - tr(Hti %*% dHtP(p) %*% Hti %*% dHtE[[l]]))
  Cpp <- function(aa, bb) 0.5 * (tr(Hti %*% d2Ht_pp(aa, bb)) -
    tr(Hti %*% dHtP(aa) %*% Hti %*% dHtP(bb)))
  .McD <- lapply(1:np, McolData)
  etaP <- matrix(vapply(1:np, function(p) as.numeric(-HfInv %*% McolEBE(p)), numeric(neta)), nrow = neta)
  eta2 <- function(aa, bb) { b <- SvecEBE(aa, bb) + SmatEBE(aa)[, ei, drop = FALSE] %*% etaP[, bb] +
      SmatEBE(bb)[, ei, drop = FALSE] %*% etaP[, aa]
    for (l in ei) b[l] <- b[l] + as.numeric(t(etaP[, aa]) %*% Tnf[l, ei, ei] %*% etaP[, bb]); as.numeric(-HfInv %*% b) }
  .CpeRow <- lapply(1:np, function(p) vapply(ei, function(l) Cpe(p, l), numeric(1)))
  R <- matrix(0, np, np)
  for (aa in 1:np) for (bb in aa:np) {
    e_ab <- eta2(aa, bb)
    dat <- d2Phi(aa, bb) + sum(.McD[[aa]] * etaP[, bb]) + sum(.McD[[bb]] * etaP[, aa]) +
           as.numeric(t(etaP[, aa]) %*% H %*% etaP[, bb]) + sum(gPhi * e_ab)
    ld <- Cpp(aa, bb) + sum(.CpeRow[[aa]] * etaP[, bb]) + sum(.CpeRow[[bb]] * etaP[, aa]) +
          as.numeric(t(etaP[, aa]) %*% Cee %*% etaP[, bb]) + sum(Cen * e_ab)
    R[aa, bb] <- R[bb, aa] <- dat + ld
  }
  R
}

#' d2R0/(deta_l deta_m) for the FOCE eta-block H: 0 (nonmem, R0 frozen) or live E$AR (foce+).
#' @noRd
E_ARelm <- function(E, l, m, fp) if (fp) E$AR[, l, m] else 0

#' C++/Armadillo port of `.foceiAnalyticSubjectRfoceFR` (FOCE (f,R) observed information).
#' Resolves the frozen-R0 sensitivities (eta-block aRe/ARe: 0 for nonmem, live E for foce+;
#' parameter columns aRc/ARc: E0 for nonmem, live E for foce+), reshapes Ath, and calls the
#' kernel.  Matches `.foceiAnalyticSubjectRfoceFR` exactly.
#' @noRd
.foceiAnalyticSubjectRfoceFRCpp <- function(E, ehat, Om, neta, ndirP, dirP, omd, ndir,
                                            Oi = solve(Om), E0 = NULL, foceType = 0L,
                                            dvSens = matrix(0, length(E$f), 0L),
                                            dvSens2 = matrix(0, length(E$f), 0L),
                                            censv = integer(0), limv = numeric(0)) {
  nobs <- length(E$f); nom <- omd$nom
  .fp <- identical(as.integer(foceType), 1L) || is.null(E0)
  if (.fp) { R0 <- E$R; aRe <- E$aR; aRc <- E$aR; ARc <- E$AR; ARe <- E$AR }
  else { R0 <- E0$R; aRe <- matrix(0, nobs, ndir); aRc <- E0$aR; ARc <- E0$AR; ARe <- array(0, c(nobs, ndir, ndir)) }
  AthC <- array(E$Ath, c(nobs, neta, ndir * ndir))
  dOiC <- array(0, c(neta, neta, max(nom, 1L)))
  if (nom > 0L) for (k in seq_len(nom)) dOiC[, , k] <- omd$dOi[[k]]
  d2OiC <- array(0, c(neta, neta, max(nom * nom, 1L)))
  if (nom > 0L) for (aa in seq_len(nom)) for (bb in seq_len(nom)) d2OiC[, , (aa - 1L) * nom + bb] <- omd$d2Oi[[aa]][[bb]]
  d2LD <- if (nom > 0L) omd$d2LD else matrix(0, 1, 1)
  foceiSubjectRfoceFR_(E$a, E$A, AthC, aRe, aRc, ARe, ARc, dvSens, dvSens2, as.integer(censv), as.numeric(limv),
                       E$f, E$y, R0, as.numeric(ehat), Oi,
                       dOiC, d2OiC, d2LD, neta, ndir, ndirP, nom, as.integer(dirP))
}

#' C++/Armadillo port of `.foceiAnalyticSubjectRFR` (FOCEI (f,R) observed information).
#' Reshapes the 3rd-order Ath/AthR tensors and the Omega derivative lists for the kernel.
#' @noRd
.foceiAnalyticSubjectRFRCpp <- function(E, ehat, Om, neta, ndirP, dirP, omd, ndir, Oi = solve(Om),
                                        dvSens = matrix(0, length(E$f), 0L),
                                        dvSens2 = matrix(0, length(E$f), 0L),
                                        censv = integer(0), limv = numeric(0)) {
  nobs <- length(E$f); nom <- omd$nom
  # Ath/AthR are [obs, neta, ndir, ndir] (eta axis first); reshape to (obs, neta, ndir^2)
  # so the kernel reads Ath(o, l, s + t*ndir) with l over the etas.
  AthC <- array(E$Ath, c(nobs, neta, ndir * ndir)); AthRC <- array(E$AthR, c(nobs, neta, ndir * ndir))
  dOiC <- array(0, c(neta, neta, max(nom, 1L)))
  if (nom > 0L) for (k in seq_len(nom)) dOiC[, , k] <- omd$dOi[[k]]
  d2OiC <- array(0, c(neta, neta, max(nom * nom, 1L)))
  if (nom > 0L) for (aa in seq_len(nom)) for (bb in seq_len(nom)) d2OiC[, , (aa - 1L) * nom + bb] <- omd$d2Oi[[aa]][[bb]]
  d2LD <- if (nom > 0L) omd$d2LD else matrix(0, 1, 1)
  foceiSubjectRFR_(E$a, E$A, AthC, E$aR, E$AR, AthRC, dvSens, dvSens2, as.integer(censv), as.numeric(limv),
                   E$f, E$y, E$R, as.numeric(ehat), Oi,
                   dOiC, d2OiC, d2LD, neta, ndir, ndirP, nom, as.integer(dirP))
}

#' Per-subject observed-information R over structural theta, sigma and Omega, from
#' already-evaluated sensitivities `E$a/A/Ath` + the error model `ef` + the Omega
#' derivatives `omd`.  The inner Hessian / sigma / Omega machinery is eta-indexed
#' (`ei`); the sensitivity slots (N, Tn, dHtD, d2HtDD, and the theta accessors) use
#' the direction index (`di`), and each structural theta differentiates in its own
#' direction-slot via `dirTh` (a mu-ref theta reuses its eta's direction, so a
#' fully mu-referenced model has ndir == neta).  Param order: `nth` theta, `nsg`
#' sigma, then Omega.
#' @noRd
.foceiAnalyticSubjectR <- function(E, ehat, Om, ef, neta, nth, nsg, sgVar, omd,
                                   ndir = neta, dirTh = seq_len(nth), Oi = solve(Om),
                                   interaction = 1L, E0 = NULL, foceType = 0L,
                                   qx = NULL, qw = NULL, solveNode = NULL) {
  if (identical(as.integer(interaction), 0L))
    return(.foceiAnalyticSubjectRfoce(E, ehat, Om, ef, neta, nth, nsg, sgVar, omd,
                                      ndir = ndir, dirTh = dirTh, Oi = Oi, E0 = E0,
                                      foceType = foceType))
  tr <- function(M) sum(diag(M))
  a <- E$a; A <- E$A; Ath <- E$Ath; f <- E$f; y <- E$y
  evf <- function(e) ef$ev(e, f, y)
  rd <- list(r1 = evf(ef$sc$r1), r2 = evf(ef$sc$r2), r3 = evf(ef$sc$r3))
  pf <- list(p = evf(ef$sc$p), p1 = evf(ef$sc$p1), p2 = evf(ef$sc$p2))
  np <- nth + nsg + omd$nom
  ei <- seq_len(neta); di <- seq_len(ndir)
  ae <- a[, ei, drop = FALSE]                          # eta-cols for sigma/Omega
  .dirOf <- function(p) dirTh[p]
  H <- Oi; for (l in ei) for (m in ei) H[l, m] <- H[l, m] + sum(rd$r2 * a[, l] * a[, m] + rd$r1 * A[, l, m])
  N <- matrix(0, neta, ndir); for (l in ei) for (d in di) N[l, d] <- sum(rd$r2 * a[, l] * a[, d] + rd$r1 * A[, l, d])
  HiM <- solve(H)
  Ht <- Oi; for (l in ei) for (m in ei) Ht[l, m] <- Ht[l, m] + sum(pf$p * a[, l] * a[, m]); Hti <- solve(Ht)
  ouAA <- function(v) { M <- matrix(0, neta, neta); for (l in ei) for (m in ei) M[l, m] <- sum(v * a[, l] * a[, m]); M }
  dHtD <- lapply(di, function(s) { D <- matrix(0, neta, neta); for (l in ei) for (m in ei)
    D[l, m] <- sum(pf$p1 * a[, s] * a[, l] * a[, m] + pf$p * A[, l, s] * a[, m] + pf$p * a[, l] * A[, m, s]); D })
  d2HtDD <- lapply(di, function(s) lapply(di, function(t) { D <- matrix(0, neta, neta); for (l in ei) for (m in ei)
    D[l, m] <- sum(pf$p2 * a[, s] * a[, t] * a[, l] * a[, m] +
      pf$p1 * (A[, s, t] * a[, l] * a[, m] + a[, s] * A[, l, t] * a[, m] + a[, s] * a[, l] * A[, m, t] + a[, t] * A[, l, s] * a[, m] + a[, t] * a[, l] * A[, m, s]) +
      pf$p * (Ath[, l, s, t] * a[, m] + A[, l, s] * A[, m, t] + A[, l, t] * A[, m, s] + a[, l] * Ath[, m, s, t])); D }))
  Cen <- vapply(ei, function(l) 0.5 * tr(Hti %*% dHtD[[l]]), numeric(1))
  Cee <- matrix(0, neta, neta); for (s in ei) for (t in ei)
    Cee[s, t] <- 0.5 * (tr(Hti %*% d2HtDD[[s]][[t]]) - tr(Hti %*% dHtD[[s]] %*% Hti %*% dHtD[[t]]))
  Tn <- array(0, c(neta, ndir, ndir)); for (l in ei) for (s in di) for (t in di)
    Tn[l, s, t] <- sum(rd$r3 * a[, l] * a[, s] * a[, t] + rd$r2 * (A[, l, s] * a[, t] + A[, l, t] * a[, s] + A[, s, t] * a[, l]) + rd$r1 * Ath[, l, s, t])
  Ndd <- function(da, db) sum(rd$r2 * a[, da] * a[, db] + rd$r1 * A[, da, db])
  typ <- function(p) if (p <= nth) "th" else if (p <= nth + nsg) "sg" else "om"
  sgi <- function(p) sgVar[p - nth]
  omc <- function(p) p - nth - nsg
  PVper <- function(p) lapply(ef$per[[sgi(p)]], evf)
  PVpair <- function(aa, bb) { s1 <- sgi(aa); s2 <- sgi(bb); key <- if (paste0(s1, s2) %in% names(ef$pair)) paste0(s1, s2) else paste0(s2, s1)
    lapply(ef$pair[[key]], evf) }
  # Omega enters only the prior (1/2 eta' Omega^-1 eta + 1/2 ln|Omega|) and H~'s
  # +Omega^-1 term, so every Omega-block quantity is an E-basis contraction from
  # `omd` (non-Cholesky variance-covariance derivatives) -- diagonal or block.
  Mcol <- function(p) { t <- typ(p)
    if (t == "th") return(N[, .dirOf(p)]); if (t == "sg") return(as.numeric(crossprod(ae, PVper(p)$rf)))
    as.numeric(omd$dOi[[omc(p)]] %*% ehat) }
  dHt_p <- function(p) { t <- typ(p)
    if (t == "th") return(dHtD[[.dirOf(p)]]); if (t == "sg") return(ouAA(PVper(p)$ps))
    omd$dOi[[omc(p)]] }
  d2HtEtaP <- function(p, l) { t <- typ(p)
    if (t == "th") return(d2HtDD[[.dirOf(p)]][[l]]); if (t == "om") return(matrix(0, neta, neta))
    P <- PVper(p); D <- matrix(0, neta, neta); for (s in ei) for (m in ei)
      D[s, m] <- sum(P$pf * a[, l] * a[, s] * a[, m] + P$ps * (A[, s, l] * a[, m] + a[, s] * A[, m, l])); D }
  d2Ht_pp <- function(aa, bb) { ta <- typ(aa); tb <- typ(bb)
    if (ta == "th" && tb == "th") return(d2HtDD[[.dirOf(aa)]][[.dirOf(bb)]])
    if (ta == "om" && tb == "om") return(omd$d2Oi[[omc(aa)]][[omc(bb)]])
    if (ta == "om" || tb == "om") return(matrix(0, neta, neta))
    if (ta == "sg" && tb == "sg") return(ouAA(PVpair(aa, bb)$pss))
    thp <- if (ta == "th") aa else bb; sg <- if (ta == "th") bb else aa; d2HtEtaP(sg, .dirOf(thp)) }
  Smat <- function(p) { t <- typ(p)
    if (t == "th") { M <- matrix(0, neta, ndir); for (l in ei) for (s in di) M[l, s] <- Tn[l, .dirOf(p), s]; return(M) }
    if (t == "om") return(omd$dOi[[omc(p)]])
    P <- PVper(p); M <- matrix(0, neta, ndir); for (l in ei) for (s in di) M[l, s] <- sum(P$rff * a[, s] * a[, l] + P$rf * A[, l, s]); M }
  Svec <- function(aa, bb) { ta <- typ(aa); tb <- typ(bb)
    if (ta == "th" && tb == "th") return(Tn[, .dirOf(aa), .dirOf(bb)])
    if (ta == "om" && tb == "om") return(as.numeric(omd$d2Oi[[omc(aa)]][[omc(bb)]] %*% ehat))
    if (ta == "om" || tb == "om") return(rep(0, neta))
    if (ta == "sg" && tb == "sg") return(as.numeric(crossprod(ae, PVpair(aa, bb)$rfss)))
    thp <- if (ta == "th") aa else bb; sg <- if (ta == "th") bb else aa; Smat(sg)[, .dirOf(thp)] }
  d2Phi <- function(aa, bb) { ta <- typ(aa); tb <- typ(bb)
    if (ta == "th" && tb == "th") return(Ndd(.dirOf(aa), .dirOf(bb)))
    if (ta == "om" && tb == "om") return(0.5 * as.numeric(t(ehat) %*% omd$d2Oi[[omc(aa)]][[omc(bb)]] %*% ehat) + 0.5 * omd$d2LD[omc(aa), omc(bb)])
    if (ta == "om" || tb == "om") return(0)
    if (ta == "sg" && tb == "sg") return(sum(PVpair(aa, bb)$rss))
    thp <- if (ta == "th") aa else bb; sg <- if (ta == "th") bb else aa; as.numeric(crossprod(a[, .dirOf(thp)], PVper(sg)$rf)) }
  npE <- np                                        # full theta + sigma + Omega block
  .Mcols <- lapply(1:npE, Mcol)                     # Mcol(p) depends on ONE index -- cache
  etaP <- matrix(vapply(1:npE, function(p) as.numeric(-HiM %*% .Mcols[[p]]), numeric(neta)), nrow = neta)  # neta x npE (neta==1 safe)
  eta2 <- function(aa, bb) { b <- Svec(aa, bb) + Smat(aa)[, ei, drop = FALSE] %*% etaP[, bb] + Smat(bb)[, ei, drop = FALSE] %*% etaP[, aa]
    for (l in ei) b[l] <- b[l] + as.numeric(t(etaP[, aa]) %*% Tn[l, ei, ei] %*% etaP[, bb]); as.numeric(-HiM %*% b) }
  Cpe <- function(p, l) 0.5 * (tr(Hti %*% d2HtEtaP(p, l)) - tr(Hti %*% dHt_p(p) %*% Hti %*% dHtD[[l]]))
  Cpp <- function(aa, bb) 0.5 * (tr(Hti %*% d2Ht_pp(aa, bb)) - tr(Hti %*% dHt_p(aa) %*% Hti %*% dHt_p(bb)))
  .CpeRow <- lapply(1:npE, function(p) vapply(ei, function(l) Cpe(p, l), numeric(1)))  # length-neta row per p
  .ldOf <- function(aa, bb, e2) Cpp(aa, bb) + sum(.CpeRow[[aa]] * etaP[, bb]) +
    sum(.CpeRow[[bb]] * etaP[, aa]) + as.numeric(t(etaP[, aa]) %*% Cee %*% etaP[, bb]) + sum(Cen * e2)
  if (is.null(qx)) {                                  # ---- FOCEI/Laplace (nAGQ == 1) ----
    R <- matrix(0, npE, npE)
    for (aa in 1:npE) for (bb in aa:npE) {              # R is symmetric -- fill upper, mirror
      dat <- d2Phi(aa, bb) - as.numeric(t(.Mcols[[aa]]) %*% HiM %*% .Mcols[[bb]])
      R[aa, bb] <- R[bb, aa] <- dat + .ldOf(aa, bb, eta2(aa, bb))
    }
    return(R)
  }
  # ---- AGQ (nAGQ > 1) ------------------------------------------------------------------
  # L_i = -log(sum_k a_k) + 0.5*log|Ht|, so only the DATA half changes; `ld` above is the
  # log-det half and is reused verbatim:
  #     R_i = ld + E_pi[Phi_pq^total(k)] - Cov_pi(Phi_p^total, Phi_q^total)
  # At nAGQ=1 there is one node at etahat with pi=1: Cov vanishes and E_pi[.] collapses to
  # the total d2Phi(etahat) = d2Phi - M'H^-1 M, i.e. exactly the `dat` branch above.
  .phiU <- function(M) { U <- M; U[lower.tri(U)] <- 0; diag(U) <- diag(U) / 2; U }
  .ch <- tryCatch(chol(Ht), error = function(e) NULL)   # non-PD Ht -> caller falls back to FD
  if (is.null(.ch)) return(NULL)
  # PD margin, as in the gradient path: chol() succeeding is not proof the objective placed
  # the nodes via chol(Ht).  arma's is_sympd() and R's chol() disagree near the PD boundary,
  # and at fit time calcEtaHessian switches to nmNearPD/cholSE when !is_sympd() -- so a
  # near-boundary Ht would mean the cov differentiates a different node placement than was
  # optimized.  Require the margin and fall back to FD otherwise.
  if (!is.finite(rcond(Ht)) || rcond(Ht) < 1e-10) return(NULL)
  Ginv <- backsolve(.ch, diag(neta))                    # chol(Ht)^-1 places the nodes
  # TOTAL derivatives of Ht through etahat(theta) -- the node placement moves with the mode
  dHtStar <- lapply(1:npE, function(p) { D <- dHt_p(p)
    for (l in ei) D <- D + dHtD[[l]] * etaP[l, p]; D })
  .e2 <- lapply(1:npE, function(aa) vector("list", npE))
  for (aa in 1:npE) for (bb in aa:npE) { v <- eta2(aa, bb); .e2[[aa]][[bb]] <- v; .e2[[bb]][[aa]] <- v }
  d2HtStar <- function(aa, bb) { D <- d2Ht_pp(aa, bb)
    for (l in ei) D <- D + d2HtEtaP(aa, l) * etaP[l, bb] + d2HtEtaP(bb, l) * etaP[l, aa]
    for (l in ei) for (m in ei) D <- D + d2HtDD[[l]][[m]] * etaP[l, aa] * etaP[m, bb]
    for (l in ei) D <- D + dHtD[[l]] * .e2[[aa]][[bb]][l]; D }
  # dGinv_p = -Ginv*PhiU(Ginv' dHtStar_p Ginv); differentiate once more for d2Ginv
  UL <- lapply(1:npE, function(p) .phiU(t(Ginv) %*% dHtStar[[p]] %*% Ginv))
  dGinvL <- lapply(1:npE, function(p) -Ginv %*% UL[[p]])
  d2GinvF <- function(aa, bb) Ginv %*% UL[[bb]] %*% UL[[aa]] -
    Ginv %*% .phiU(t(dGinvL[[bb]]) %*% dHtStar[[aa]] %*% Ginv + t(Ginv) %*% d2HtStar(aa, bb) %*% Ginv +
                   t(Ginv) %*% dHtStar[[aa]] %*% dGinvL[[bb]])
  # Phi and its p/eta derivatives at a node.  Phi_eta(node) != 0 (unlike at etahat), which
  # is what makes the d2etaCur (w_kpq) term survive here and vanish in the FOCEI branch.
  .nodeParts <- function(Ek, etak) {
    ak <- Ek$a; Ak <- Ek$A; fk <- Ek$f
    evk <- function(e) ef$ev(e, fk, y)
    r1k <- evk(ef$sc$r1); r2k <- evk(ef$sc$r2)
    aek <- ak[, ei, drop = FALSE]
    Hk <- Oi; for (l in ei) for (m in ei) Hk[l, m] <- Hk[l, m] + sum(r2k * ak[, l] * ak[, m] + r1k * Ak[, l, m])
    Nk <- matrix(0, neta, ndir); for (l in ei) for (d in di) Nk[l, d] <- sum(r2k * ak[, l] * ak[, d] + r1k * Ak[, l, d])
    gk <- as.numeric(Oi %*% etak) + vapply(ei, function(l) sum(r1k * ak[, l]), numeric(1))
    PVperk <- function(p) lapply(ef$per[[sgi(p)]], evk)
    PVpairk <- function(aa, bb) { s1 <- sgi(aa); s2 <- sgi(bb)
      key <- if (paste0(s1, s2) %in% names(ef$pair)) paste0(s1, s2) else paste0(s2, s1)
      lapply(ef$pair[[key]], evk) }
    Mk <- function(p) { t <- typ(p)
      if (t == "th") return(Nk[, .dirOf(p)])
      if (t == "sg") return(as.numeric(crossprod(aek, PVperk(p)$rf)))
      as.numeric(omd$dOi[[omc(p)]] %*% etak) }
    dPk <- function(p) { t <- typ(p)
      if (t == "th") return(sum(r1k * ak[, .dirOf(p)]))
      if (t == "sg") return(sum(PVperk(p)$rs))
      0.5 * as.numeric(t(etak) %*% omd$dOi[[omc(p)]] %*% etak) + 0.5 * omd$dLD[omc(p)] }
    d2Pk <- function(aa, bb) { ta <- typ(aa); tb <- typ(bb)
      if (ta == "th" && tb == "th")
        return(sum(r2k * ak[, .dirOf(aa)] * ak[, .dirOf(bb)] + r1k * Ak[, .dirOf(aa), .dirOf(bb)]))
      if (ta == "om" && tb == "om")
        return(0.5 * as.numeric(t(etak) %*% omd$d2Oi[[omc(aa)]][[omc(bb)]] %*% etak) +
                 0.5 * omd$d2LD[omc(aa), omc(bb)])
      if (ta == "om" || tb == "om") return(0)
      if (ta == "sg" && tb == "sg") return(sum(PVpairk(aa, bb)$rss))
      thp <- if (ta == "th") aa else bb; sg <- if (ta == "th") bb else aa
      as.numeric(crossprod(ak[, .dirOf(thp)], PVperk(sg)$rf)) }
    list(H = Hk, g = gk, M = Mk, dP = dPk, d2P = d2Pk,
         Phi = sum(evk(ef$sc$r0)) + 0.5 * as.numeric(t(etak) %*% Oi %*% etak))
  }
  nn <- nrow(qx)
  P <- vector("list", nn); lg <- numeric(nn); vK <- vector("list", nn)
  for (k in seq_len(nn)) {
    x <- qx[k, ]
    etak <- ehat + sqrt(2) * as.numeric(Ginv %*% x)   # sqrt(2)*Ginv*x, matching inner.cpp
    Ek <- solveNode(etak)
    if (is.null(Ek)) return(NULL)                       # node solve failed -> caller drops to FD
    P[[k]] <- .nodeParts(Ek, etak)
    lg[k] <- sum(log(qw[k, ])) + sum(x^2) - P[[k]]$Phi   # u_k = log a_k (exp(x'x) untilt)
    # matrix(, nrow=neta) is load-bearing: at neta==1 vapply drops the dim, vK[[k]][, p]
    # throws and the tryCatch silently falls back to FD.  Same guard as etaP above.
    vK[[k]] <- matrix(vapply(1:npE, function(p) etaP[, p] + sqrt(2) * as.numeric(dGinvL[[p]] %*% x),
                             numeric(neta)), nrow = neta)
  }
  lg <- lg - max(lg); pk <- exp(lg); pk <- pk / sum(pk)  # pi_k (softmax, shifted for stability)
  gp <- matrix(0, nn, npE)                               # Phi_p^total at each node
  for (k in seq_len(nn)) for (p in 1:npE)
    gp[k, p] <- P[[k]]$dP(p) + sum(P[[k]]$g * vK[[k]][, p])
  Eg <- as.numeric(crossprod(pk, gp))
  R <- matrix(0, npE, npE)
  for (aa in 1:npE) for (bb in aa:npE) {
    d2G <- d2GinvF(aa, bb); e2 <- .e2[[aa]][[bb]]; s <- 0
    for (k in seq_len(nn)) {
      va <- vK[[k]][, aa]; vb <- vK[[k]][, bb]
      wab <- e2 + sqrt(2) * as.numeric(d2G %*% qx[k, ])   # d2etaCur_k/dp dq (sqrt(2) on d2Ginv*x)
      Pk <- P[[k]]
      s <- s + pk[k] * (Pk$d2P(aa, bb) + sum(Pk$M(aa) * vb) + sum(Pk$M(bb) * va) +
                          as.numeric(t(va) %*% Pk$H %*% vb) + sum(Pk$g * wab))
    }
    cv <- sum(pk * gp[, aa] * gp[, bb]) - Eg[aa] * Eg[bb]
    R[aa, bb] <- R[bb, aa] <- .ldOf(aa, bb, e2) + s - cv
  }
  R
}

#' Per-subject observed-information R for FOCE (interaction=0).
#'
#' Unlike FOCEI, the FOCE EBE eta-hat does NOT stationarize the full Laplace
#' objective Phi (the inner problem drops the dR/deta interaction term), so the
#' envelope/Schur shortcut used in [.foceiAnalyticSubjectR] fails.  This uses the
#' GENERAL total-derivative Hessian of F = Phi_full + 0.5 log|H~_FOCE|:
#'
#'   R_ab = F_ab + F_aeta eta_b + F_beta eta_a + eta_a' F_etaeta eta_b + F_eta eta_ab
#'
#' with eta_b, eta_ab the 1st/2nd parameter-derivatives of the EBE obtained by
#' implicit differentiation of the FOCE inner stationarity S_FOCE = sum(q a) +
#' Omega^-1 eta = 0 (q = -eps/R, the interaction-free gradient coefficient; its
#' Jacobian is the FOCE inner Hessian Hf = sum(q' a a' + q A) + Omega^-1 and its
#' 3-tensor Tnf).  Phi keeps the FULL rho derivatives (r1/r2 = Phi_eta / Phi_etaeta
#' and Phi_eta is nonzero here); the log-determinant uses p = 1/R (no interaction).
#' Splitting F into the Phi (data) and 0.5 log|H~| (log-det) parts, the data term is
#' the general (non-envelope) form and the log-det term is the same moving-mode
#' assembly as FOCEI but built from p=1/R.  Same parameter order as
#' [.foceiAnalyticSubjectR].
#' @noRd
.foceiAnalyticSubjectRfoce <- function(E, ehat, Om, ef, neta, nth, nsg, sgVar, omd,
                                       ndir = neta, dirTh = seq_len(nth), Oi = solve(Om), E0 = NULL,
                                       foceType = 0L) {
  tr <- function(M) sum(diag(M))
  a <- E$a; A <- E$A; Ath <- E$Ath; f <- E$f; y <- E$y
  # CORRECTED FOCE (foceType=0): the variance R0 is frozen at the eta=0 POPULATION
  # prediction f0 (E0$f), with population sensitivities a0/A0 (E0$a/E0$A) supplying R0's
  # theta-chain.  The residual/numerator/inner sensitivities (y-f, f, a, A, Ath) stay at
  # eta-hat.  R0 is eta-independent (dR0/deta=0), so q1=1/R0, q2=0, pF1=0, r3=0; its
  # theta/sigma derivatives (through f0) enter as a0-chain corrections on the `th`
  # accessors only.  "foce+" (foceType=1) keeps the live conditional R: the live-R
  # pieces carry the dR/df chain through f's own sensitivities, so no f0 solve/chain.
  .fc <- if (identical(as.integer(foceType), 1L)) ef$focePlus else ef$foce
  f0 <- if (!is.null(E0)) E0$f else f
  evf <- function(e) ef$ev(e, f, y, f0)
  rd <- list(r1 = evf(.fc$sc$r1), r2 = evf(.fc$sc$r2), r3 = evf(.fc$sc$r3))   # full rho (Phi), R0
  qd <- list(q0 = evf(.fc$sc$q0), q1 = evf(.fc$sc$q1), q2 = evf(.fc$sc$q2))   # FOCE inner gradient coef, R0
  pf <- list(p = evf(.fc$sc$pF), p1 = evf(.fc$sc$pF1), p2 = evf(.fc$sc$pF2))  # FOCE determinant p = 1/R0
  np <- nth + nsg + omd$nom
  ei <- seq_len(neta); di <- seq_len(ndir); ae <- a[, ei, drop = FALSE]
  .dirOf <- function(p) dirTh[p]
  # a0-chain (only when the variance depends on the population prediction, i.e. prop part)
  .cf0 <- isTRUE(.fc$dependsF0) && !is.null(E0)
  a0 <- if (.cf0) E0$a else NULL; A0 <- if (.cf0) E0$A else NULL
  fq <- if (.cf0) list(qf0 = evf(.fc$f0$qf0), qff0 = evf(.fc$f0$qff0), qf0f0 = evf(.fc$f0$qf0f0),
                       pFf0 = evf(.fc$f0$pFf0), pFf0f0 = evf(.fc$f0$pFf0f0),
                       rhof0 = evf(.fc$f0$rhof0), rhof0f0 = evf(.fc$f0$rhof0f0)) else NULL

  # ---- Phi (data) pieces: FULL rho derivatives ----
  H <- Oi; for (l in ei) for (m in ei) H[l, m] <- H[l, m] + sum(rd$r2 * a[, l] * a[, m] + rd$r1 * A[, l, m])  # Phi_etaeta
  Ndat <- matrix(0, neta, ndir); for (l in ei) for (d in di) Ndat[l, d] <- sum(rd$r2 * a[, l] * a[, d] + rd$r1 * A[, l, d])  # Phi_(eta,theta)
  gPhi <- as.numeric(Oi %*% ehat); for (l in ei) gPhi[l] <- gPhi[l] + sum(rd$r1 * a[, l])  # Phi_eta (nonzero at eta-hat_FOCE)

  # ---- FOCE inner (EBE) pieces: q-based Jacobian S_eta = Hf and its 3-tensor ----
  Hf <- Oi; for (l in ei) for (m in ei) Hf[l, m] <- Hf[l, m] + sum(qd$q1 * a[, l] * a[, m] + qd$q0 * A[, l, m])  # S_eta = Hf
  HfInv <- solve(Hf)
  Nf <- matrix(0, neta, ndir); for (l in ei) for (d in di) Nf[l, d] <- sum(qd$q1 * a[, l] * a[, d] + qd$q0 * A[, l, d])  # S_(eta,theta)
  Tnf <- array(0, c(neta, ndir, ndir)); for (l in ei) for (s in di) for (t in di)
    Tnf[l, s, t] <- sum(qd$q2 * a[, l] * a[, s] * a[, t] + qd$q1 * (A[, l, s] * a[, t] + A[, l, t] * a[, s] + A[, s, t] * a[, l]) + qd$q0 * Ath[, l, s, t])  # S_etaeta

  # ---- determinant H~ (p = 1/R) pieces (same moving-mode assembly as FOCEI) ----
  Ht <- Oi; for (l in ei) for (m in ei) Ht[l, m] <- Ht[l, m] + sum(pf$p * a[, l] * a[, m]); Hti <- solve(Ht)
  ouAA <- function(v) { M <- matrix(0, neta, neta); for (l in ei) for (m in ei) M[l, m] <- sum(v * a[, l] * a[, m]); M }
  dHtD <- lapply(di, function(s) { D <- matrix(0, neta, neta); for (l in ei) for (m in ei)
    D[l, m] <- sum(pf$p1 * a[, s] * a[, l] * a[, m] + pf$p * A[, l, s] * a[, m] + pf$p * a[, l] * A[, m, s]); D })
  d2HtDD <- lapply(di, function(s) lapply(di, function(t) { D <- matrix(0, neta, neta); for (l in ei) for (m in ei)
    D[l, m] <- sum(pf$p2 * a[, s] * a[, t] * a[, l] * a[, m] +
      pf$p1 * (A[, s, t] * a[, l] * a[, m] + a[, s] * A[, l, t] * a[, m] + a[, s] * a[, l] * A[, m, t] + a[, t] * A[, l, s] * a[, m] + a[, t] * a[, l] * A[, m, s]) +
      pf$p * (Ath[, l, s, t] * a[, m] + A[, l, s] * A[, m, t] + A[, l, t] * A[, m, s] + a[, l] * Ath[, m, s, t])); D }))
  Cen <- vapply(ei, function(l) 0.5 * tr(Hti %*% dHtD[[l]]), numeric(1))
  Cee <- matrix(0, neta, neta); for (s in ei) for (t in ei)
    Cee[s, t] <- 0.5 * (tr(Hti %*% d2HtDD[[s]][[t]]) - tr(Hti %*% dHtD[[s]] %*% Hti %*% dHtD[[t]]))

  typ <- function(p) if (p <= nth) "th" else if (p <= nth + nsg) "sg" else "om"
  sgi <- function(p) sgVar[p - nth]; omc <- function(p) p - nth - nsg
  perR <- function(p) lapply(.fc$per[[sgi(p)]][c("rf", "rff")], evf)                 # rho sigma partials (Phi), R0
  perQ <- function(p) list(qs = evf(.fc$per[[sgi(p)]]$qs), qsf = evf(.fc$per[[sgi(p)]]$qsf))  # FOCE inner sigma, R0
  perP <- function(p) list(ps = evf(.fc$per[[sgi(p)]]$psF), pf = evf(.fc$per[[sgi(p)]]$pfF))   # FOCE det sigma, R0
  # sigma x population-prediction cross fields (a0-chain of the sigma x theta blocks)
  PVf0 <- function(p) lapply(.fc$perf0[[sgi(p)]], evf)                               # qf0s, pFf0s, rhof0s
  pairKey <- function(aa, bb) { s1 <- sgi(aa); s2 <- sgi(bb); if (paste0(s1, s2) %in% names(.fc$pair)) paste0(s1, s2) else paste0(s2, s1) }
  pairR <- function(aa, bb) lapply(.fc$pair[[pairKey(aa, bb)]][c("rss", "rfss")], evf)
  pairQ <- function(aa, bb) list(qss = evf(.fc$pair[[pairKey(aa, bb)]]$qss))
  pairP <- function(aa, bb) list(pss = evf(.fc$pair[[pairKey(aa, bb)]]$pssF))

  # ---- data-side accessors (Phi partials) ----
  # a0-correction to a theta-derivative of q0=Phi_f: d(q0)/d(theta) picks up qf0*a0 (the
  # population prediction moves with theta); rho_ff0 == qf0.
  McolData <- function(p) { t <- typ(p)                                            # Phi_(eta,p)
    if (t == "th") { v <- Ndat[, .dirOf(p)]
      if (.cf0) v <- v + as.numeric(crossprod(ae, fq$qf0 * a0[, .dirOf(p)])); return(v) }
    if (t == "sg") return(as.numeric(crossprod(ae, perR(p)$rf)))
    as.numeric(omd$dOi[[omc(p)]] %*% ehat) }
  d2Phi <- function(aa, bb) { ta <- typ(aa); tb <- typ(bb)                          # Phi_(p,p) at fixed eta
    if (ta == "th" && tb == "th") { da <- .dirOf(aa); db <- .dirOf(bb)
      v <- sum(rd$r2 * a[, da] * a[, db] + rd$r1 * A[, da, db])
      if (.cf0) v <- v + sum(fq$qf0 * (a[, da] * a0[, db] + a0[, da] * a[, db]) +
                             fq$rhof0f0 * a0[, da] * a0[, db] + fq$rhof0 * A0[, da, db])
      return(v) }
    if (ta == "om" && tb == "om") return(0.5 * as.numeric(t(ehat) %*% omd$d2Oi[[omc(aa)]][[omc(bb)]] %*% ehat) + 0.5 * omd$d2LD[omc(aa), omc(bb)])
    if (ta == "om" || tb == "om") return(0)
    if (ta == "sg" && tb == "sg") return(sum(pairR(aa, bb)$rss))
    thp <- if (ta == "th") aa else bb; sg <- if (ta == "th") bb else aa
    v <- as.numeric(crossprod(a[, .dirOf(thp)], perR(sg)$rf))
    if (.cf0) v <- v + as.numeric(crossprod(a0[, .dirOf(thp)], PVf0(sg)$rhof0s)); v }

  # ---- EBE-side accessors (S_FOCE partials); S=Phi_eta, so S_p carries the same qf0*a0
  # theta-chain as McolData; the eta-role indices (contracted with etaP) get no a0 term ----
  McolEBE <- function(p) { t <- typ(p)                                             # S_p
    if (t == "th") { v <- Nf[, .dirOf(p)]
      if (.cf0) v <- v + as.numeric(crossprod(ae, fq$qf0 * a0[, .dirOf(p)])); return(v) }
    if (t == "sg") return(as.numeric(crossprod(ae, perQ(p)$qs)))
    as.numeric(omd$dOi[[omc(p)]] %*% ehat) }
  SmatEBE <- function(p) { t <- typ(p)                                             # S_(p,eta), eta-role s
    if (t == "th") { d <- .dirOf(p); M <- matrix(0, neta, ndir)
      for (l in ei) for (s in di) M[l, s] <- Tnf[l, d, s]
      if (.cf0) for (l in ei) for (s in ei)
        M[l, s] <- M[l, s] + sum(fq$qff0 * a[, s] * a0[, d] * a[, l] + fq$qf0 * a0[, d] * A[, l, s])
      return(M) }
    if (t == "om") return(omd$dOi[[omc(p)]])
    P <- perQ(p); M <- matrix(0, neta, ndir); for (l in ei) for (s in di) M[l, s] <- sum(P$qsf * a[, s] * a[, l] + P$qs * A[, l, s]); M }
  SvecEBE <- function(aa, bb) { ta <- typ(aa); tb <- typ(bb)                       # S_(p,p)
    if (ta == "th" && tb == "th") { da <- .dirOf(aa); db <- .dirOf(bb); v <- Tnf[, da, db]
      if (.cf0) { w <- fq$qff0 * (a[, da] * a0[, db] + a0[, da] * a[, db]) + fq$qf0f0 * a0[, da] * a0[, db] + fq$qf0 * A0[, da, db]
        for (l in ei) v[l] <- v[l] + sum(w * a[, l] + fq$qf0 * (a0[, da] * A[, l, db] + a0[, db] * A[, l, da])) }
      return(v) }
    if (ta == "om" && tb == "om") return(as.numeric(omd$d2Oi[[omc(aa)]][[omc(bb)]] %*% ehat))
    if (ta == "om" || tb == "om") return(rep(0, neta))
    if (ta == "sg" && tb == "sg") return(as.numeric(crossprod(ae, pairQ(aa, bb)$qss)))
    thp <- if (ta == "th") aa else bb; sg <- if (ta == "th") bb else aa
    v <- SmatEBE(sg)[, .dirOf(thp)]
    if (.cf0) v <- v + as.numeric(crossprod(ae, PVf0(sg)$qf0s * a0[, .dirOf(thp)])); v }

  # ---- determinant-side accessors (H~ partials, p = 1/R0).  dHtD/d2HtDD are eta-role
  # (pF1=pF2=0), so a theta-derivative of pF=1/R0 adds a pFf0*a0 chain on the `th` blocks. ----
  dHt_p <- function(p) { t <- typ(p)
    if (t == "th") { D <- dHtD[[.dirOf(p)]]
      if (.cf0) { d <- .dirOf(p); for (m in ei) for (n in ei) D[m, n] <- D[m, n] + sum(fq$pFf0 * a0[, d] * a[, m] * a[, n]) }
      return(D) }
    if (t == "sg") return(ouAA(perP(p)$ps)); omd$dOi[[omc(p)]] }
  d2HtEtaP <- function(p, l) { t <- typ(p)
    if (t == "th") { D <- d2HtDD[[.dirOf(p)]][[l]]
      if (.cf0) { d <- .dirOf(p); for (m in ei) for (n in ei)
        D[m, n] <- D[m, n] + sum(fq$pFf0 * a0[, d] * (A[, m, l] * a[, n] + a[, m] * A[, n, l])) }
      return(D) }
    if (t == "om") return(matrix(0, neta, neta))
    P <- perP(p); D <- matrix(0, neta, neta); for (s in ei) for (m in ei)
      D[s, m] <- sum(P$pf * a[, l] * a[, s] * a[, m] + P$ps * (A[, s, l] * a[, m] + a[, s] * A[, m, l])); D }
  d2Ht_pp <- function(aa, bb) { ta <- typ(aa); tb <- typ(bb)
    if (ta == "th" && tb == "th") { da <- .dirOf(aa); db <- .dirOf(bb); D <- d2HtDD[[da]][[db]]
      if (.cf0) for (m in ei) for (n in ei)
        D[m, n] <- D[m, n] + sum(fq$pFf0f0 * a0[, da] * a0[, db] * a[, m] * a[, n] + fq$pFf0 * A0[, da, db] * a[, m] * a[, n] +
          fq$pFf0 * (a0[, da] * A[, m, db] + a0[, db] * A[, m, da]) * a[, n] +
          fq$pFf0 * (a0[, da] * A[, n, db] + a0[, db] * A[, n, da]) * a[, m])
      return(D) }
    if (ta == "om" && tb == "om") return(omd$d2Oi[[omc(aa)]][[omc(bb)]])
    if (ta == "om" || tb == "om") return(matrix(0, neta, neta))
    if (ta == "sg" && tb == "sg") return(ouAA(pairP(aa, bb)$pss))
    thp <- if (ta == "th") aa else bb; sg <- if (ta == "th") bb else aa; d <- .dirOf(thp)
    D <- d2HtEtaP(sg, d)
    if (.cf0) { P <- PVf0(sg); for (m in ei) for (n in ei) D[m, n] <- D[m, n] + sum(P$pFf0s * a0[, d] * a[, m] * a[, n]) }
    D }
  Cpe <- function(p, l) 0.5 * (tr(Hti %*% d2HtEtaP(p, l)) - tr(Hti %*% dHt_p(p) %*% Hti %*% dHtD[[l]]))
  Cpp <- function(aa, bb) 0.5 * (tr(Hti %*% d2Ht_pp(aa, bb)) - tr(Hti %*% dHt_p(aa) %*% Hti %*% dHt_p(bb)))

  npE <- np
  .McD <- lapply(1:npE, McolData)                                                  # Phi_(eta,p), cached
  etaP <- matrix(vapply(1:npE, function(p) as.numeric(-HfInv %*% McolEBE(p)), numeric(neta)), nrow = neta)  # eta_p = -Hf^-1 S_p
  eta2 <- function(aa, bb) {                                                        # eta_ab (2nd EBE deriv)
    b <- SvecEBE(aa, bb) + SmatEBE(aa)[, ei, drop = FALSE] %*% etaP[, bb] + SmatEBE(bb)[, ei, drop = FALSE] %*% etaP[, aa]
    for (l in ei) b[l] <- b[l] + as.numeric(t(etaP[, aa]) %*% Tnf[l, ei, ei] %*% etaP[, bb]); as.numeric(-HfInv %*% b) }
  .CpeRow <- lapply(1:npE, function(p) vapply(ei, function(l) Cpe(p, l), numeric(1)))
  R <- matrix(0, npE, npE)
  for (aa in 1:npE) for (bb in aa:npE) {                # R is symmetric -- fill upper, mirror
    e_ab <- eta2(aa, bb)
    # data term (general, non-envelope): F_ab + F_aeta eta_b + F_beta eta_a + eta_a' F_etaeta eta_b + F_eta eta_ab, Phi part
    dat <- d2Phi(aa, bb) + sum(.McD[[aa]] * etaP[, bb]) + sum(.McD[[bb]] * etaP[, aa]) +
           as.numeric(t(etaP[, aa]) %*% H %*% etaP[, bb]) + sum(gPhi * e_ab)
    # log-determinant term (moving mode), 0.5 log|H~| part
    ld <- Cpp(aa, bb) + sum(.CpeRow[[aa]] * etaP[, bb]) + sum(.CpeRow[[bb]] * etaP[, aa]) +
          as.numeric(t(etaP[, aa]) %*% Cee %*% etaP[, bb]) + sum(Cen * e_ab)
    R[aa, bb] <- R[bb, aa] <- dat + ld
  }
  R
}

#' Base subject solve: `f` plus the 1st/2nd analytic sensitivities (`a`, `A`), no
#' 3rd-order Shi tensor.  Shared by [.foceiAnalyticSolveSubjectFD3] (which adds `Ath`)
#' and the FOCE EBE re-solve (which needs only `a`/`A`).  Muffles benign solver
#' warnings (a real error returns `NULL` -> FD fallback); the nrow guard bails when an
#' EVID==2/covariate-update row shares an obs timestamp (would misalign f against y).
#' @noRd
.foceiAnalyticSolveFA <- function(aug, params, ev, times, tol = 1e-10) {
  dirs <- aug$dirs; nd <- length(dirs)
  .ev <- .foceiAnalyticEvents(aug, ev)                    # reuse the pre-translated events (FOCE Newton reuse)
  .nc <- if (is.null(aug$cores)) 0L else aug$cores
  # DDE augmented solve: force pure dop853 (dense, no Jacobian).  Its 8th-order dense history
  # reproduces the delayed sensitivity solve exactly and needs no Jacobian, so it sidesteps the
  # composite/ros4 on-the-fly Jacobian generation for this THETA/ETA-named augmented model.
  .ddeArgs <- if (isTRUE(rxode2::rxModelVars(aug$augMod)$flags[["hasDelay"]] == 1L))
    list(method = "dop853", stiff2 = 0L, dense = TRUE) else list()
  .d <- tryCatch(withCallingHandlers(
      as.data.frame(do.call(rxode2::rxSolve, c(list(aug$augMod, params = params, events = .ev, cores = .nc,
          returnType = "data.frame", atol = tol, rtol = tol), .ddeArgs))),
      warning = function(w) invokeRestart("muffleWarning")),
    error = function(e) NULL)
  if (is.null(.d)) return(NULL)
  .d <- .d[.d$time %in% times, , drop = FALSE]
  .fD <- if (is.null(aug$fDirs)) dirs else aug$fDirs   # sigma-skip: f-sensitivities only over f-directions
  if (nrow(.d) == 0L || nrow(.d) != length(times) ||
        !all(c("rx_predf_", paste0("rx_f1_", .fD)) %in% names(.d))) return(NULL)
  no <- nrow(.d)
  a <- matrix(0, no, nd); a[, match(.fD, dirs)] <- vapply(.fD, function(q) .d[[paste0("rx_f1_", q)]], numeric(no))
  A <- array(0, c(no, nd, nd))
  for (r in seq_len(nrow(aug$P2))) {                   # only i<=j emitted -> mirror to i>j (f-dir pairs)
    .ii <- match(aug$P2$i[r], dirs); .jj <- match(aug$P2$j[r], dirs)
    .v2 <- .d[[paste0("rx_f2_", aug$P2$i[r], "_", aug$P2$j[r])]]
    A[, .ii, .jj] <- .v2; A[, .jj, .ii] <- .v2
  }
  .out <- list(f = .d$rx_predf_, a = a, A = A)
  if (isTRUE(aug$hasRvar)) .out <- c(.out, .foceiReadRvar(.d, aug, no))
  .out
}

#' Read the residual variance R and its 1st/2nd direction sensitivities (aR, AR)
#' from an augmented solve (rx_rvarf_ / rx_rvar1_<dir> / rx_rvar2_<i>_<j>).  The
#' `rx_..._` naming keeps even constant sensitivities as real output columns, so
#' every piece is read directly from the solve.
#' @noRd
.foceiReadRvar <- function(.d, aug, no) {
  dirs <- aug$dirs; nd <- length(dirs); P2 <- if (is.null(aug$P2r)) aug$P2 else aug$P2r   # rvar2 spans every direction
  aR <- matrix(vapply(dirs, function(q) .d[[paste0("rx_rvar1_", q)]], numeric(no)), no, nd)
  AR <- array(0, c(no, nd, nd))
  for (r in seq_len(nrow(P2))) {
    .ii <- match(P2$i[r], dirs); .jj <- match(P2$j[r], dirs)
    .v2 <- .d[[paste0("rx_rvar2_", P2$i[r], "_", P2$j[r])]]
    AR[, .ii, .jj] <- .v2; AR[, .jj, .ii] <- .v2
  }
  .out <- list(R = .d$rx_rvarf_, aR = aR, AR = AR)
  # sigma pseudo-directions: dR/dsigma (Rsig), d2R/(dsigma ddir) (RsigDir), and
  # d2R/(dsigma dsigma') (Rsig2), keyed by the error-parameter theta indices sigTh.
  .sig <- aug$sigTh
  if (length(.sig) > 0L) {
    .out$Rsig <- matrix(vapply(.sig, function(n) .d[[paste0("rx_rsig_", n, "_")]], numeric(no)), no, length(.sig))
    .out$RsigDir <- array(vapply(.sig, function(n)
      vapply(dirs, function(q) .d[[paste0("rx_rsig1_", n, "_", q)]], numeric(no)), matrix(0, no, nd)),
      c(no, nd, length(.sig)))
    .Rs2 <- array(0, c(no, length(.sig), length(.sig)))
    for (a in seq_along(.sig)) for (b in seq_along(.sig)[seq_along(.sig) >= a]) {
      .v <- .d[[paste0("rx_rsig2_", .sig[a], "_", .sig[b])]]; .Rs2[, a, b] <- .v; .Rs2[, b, a] <- .v
    }
    .out$Rsig2 <- .Rs2
  }
  if (isTRUE(aug$hasTrans))
    .out$trans <- list(yj = .d$rx_tyj_, lambda = .d$rx_tlambda_, low = .d$rx_tlow_, hi = .d$rx_thi_)
  .out
}

#' Re-solve one subject's EBE to the FOCE inner stationary point S_FOCE = sum(q a)
#' + Omega^-1 eta = 0 (q = -eps/R) via Newton on the FOCE inner Hessian
#' Hf = sum(q' a a' + q A) + Omega^-1, starting from the stored eta `eta0`.
#' nlmixr's stored FOCE-combined EBEs do NOT satisfy S_FOCE=0 (an estimation-side
#' inconsistency), so R must be formed at the re-solved eta.  For additive/FOCEI the
#' stored eta is already stationary (|S_FOCE| < `skip`) -> returns `eta0` unchanged
#' (byte no-op).  `NULL` on a solve/Newton failure -> caller falls back to FD.
#' @noRd
.foceiAnalyticFoceEbe <- function(aug, th, eta0, s, times, y, etav, R0, Oi, neta, tol,
                                  maxit = 30L, skip = 1e-3, conv = 1e-9,
                                  foceType = 0L, cens = NULL, limit = NULL) {
  ei <- seq_len(neta)
  # interaction-free FOCE inner gradient/curvature from (f,R0): q0 = -(y-f)/R0 = rho_f,
  # q1 = 1/R0 = rho_ff.  For censored (M2/M3/M4) observations q0/q1 are the EXACT censored
  # rho_f/rho_ff at the frozen R0 (censNormalPartials_) so the re-solved eta* is the censored
  # FOCE stationary point.  foce+ (foceType=1) uses the live conditional R at the trial eta;
  # nonmem freezes R0 at the eta=0 population value passed in.
  .fp <- identical(as.integer(foceType), 1L) || is.null(R0)
  # censored (M2/M3/M4) per-obs CENS + LIMIT (NA/NULL -> uncensored) and the censored-obs index
  .cv <- if (is.null(cens)) integer(length(y)) else as.integer(ifelse(is.na(cens), 0L, cens))
  .lv <- if (is.null(limit)) rep(NA_real_, length(y)) else as.numeric(limit)
  .cw <- which(.cv != 0 | is.finite(.lv))              # censored observations
  .SH <- function(eta) {                               # FOCE S_FOCE and its Jacobian Hf at eta
    E <- .foceiAnalyticSolveFA(aug, c(th, setNames(eta, etav)), s, times, tol = tol)
    if (is.null(E)) return(NULL)
    yt <- .foceiAnalyticTbsY(y, E$trans)               # DV -> rx_pred_ (transformed) scale; no-op if untransformed
    R0e <- if (.fp) E$R else R0
    q0 <- -(yt - E$f) / R0e; q1 <- 1 / R0e
    if (length(.cw)) {                                 # censored: exact rho_f/rho_ff at frozen R0
      .limt <- .foceiAnalyticTbsY(.lv, E$trans)        # transform the censoring bound like the DV
      .cp <- censNormalPartials_(.cv, yt, .limt, E$f, R0e, 2L)
      q0[.cw] <- .cp[.cw, 1]; q1[.cw] <- .cp[.cw, 3]   # cp[,1]=rho_f, cp[,3]=rho_ff
    }
    S <- as.numeric(Oi %*% eta); for (l in ei) S[l] <- S[l] + sum(q0 * E$a[, l])
    Hf <- Oi; for (l in ei) for (m in ei) Hf[l, m] <- Hf[l, m] + sum(q1 * E$a[, l] * E$a[, m] + q0 * E$A[, l, m])
    list(S = S, Hf = Hf)
  }
  eta <- eta0
  sh <- .SH(eta); if (is.null(sh)) return(NULL)
  if (!all(is.finite(sh$S))) return(NULL)              # unsolvable subject -> non-finite score
  if (max(abs(sh$S)) < skip) return(eta0)              # already FOCE-stationary (additive/FOCEI) -> no-op
  for (it in seq_len(maxit)) {
    step <- tryCatch(solve(sh$Hf, sh$S), error = function(e) NULL)
    if (is.null(step)) return(NULL)
    eta <- eta - step
    sh <- .SH(eta); if (is.null(sh)) return(NULL)
    if (!all(is.finite(sh$S))) return(NULL)            # unsolvable subject -> non-finite score
    if (max(abs(sh$S)) < conv) break
  }
  if (max(abs(sh$S)) >= conv) return(NULL)               # Newton did not converge -> FD fallback
  eta
}

#' Batched FOCE/foce+ EBE re-solve: the same interaction-free Newton as
#' [.foceiAnalyticFoceEbe] but over ALL subjects at once via [.foceiAnalyticSolveAll]
#' (one batched solve per Newton iteration instead of per-subject SolveFA).  Bit-identical
#' to the per-subject Newton; avoids the per-subject solve entirely (needed for the shared
#' `dirs` model, which solves batched but not per-subject in the fit's cov-hook context) and
#' is faster.  Returns the nsub x neta eta-hat matrix, or NULL if any subject fails to converge.
#' @noRd
.foceiAnalyticFoceEbeBatch <- function(am, th, ebes, ids, data, obsAll, obsTimes, etav, Oi, neta, tol,
                                       foceType = 0L, E0all = NULL, maxit = 30L, skip = 1e-3, conv = 1e-9) {
  nsub <- nrow(ebes); ei <- seq_len(neta)
  .fp <- identical(as.integer(foceType), 1L) || is.null(E0all)
  Y  <- lapply(obsAll, function(.o) .o$DV)
  CV <- lapply(obsAll, function(.o) if (is.null(.o$CENS)) integer(length(.o$DV)) else as.integer(ifelse(is.na(.o$CENS), 0L, .o$CENS)))
  LV <- lapply(obsAll, function(.o) if (is.null(.o$LIMIT)) rep(NA_real_, length(.o$DV)) else as.numeric(.o$LIMIT))
  .SHi <- function(E, eta_i, i) {                        # S_FOCE + Hf for subject i (censored-aware)
    yt <- .foceiAnalyticTbsY(Y[[i]], E$trans)
    R0e <- if (.fp) E$R else E0all[[i]]$R
    q0 <- -(yt - E$f) / R0e; q1 <- 1 / R0e
    .cw <- which(CV[[i]] != 0 | is.finite(LV[[i]]))
    if (length(.cw)) { .limt <- .foceiAnalyticTbsY(LV[[i]], E$trans)
      .cp <- censNormalPartials_(CV[[i]], yt, .limt, E$f, R0e, 2L); q0[.cw] <- .cp[.cw, 1]; q1[.cw] <- .cp[.cw, 3] }
    S <- as.numeric(Oi %*% eta_i); for (l in ei) S[l] <- S[l] + sum(q0 * E$a[, l])
    Hf <- Oi; for (l in ei) for (m in ei) Hf[l, m] <- Hf[l, m] + sum(q1 * E$a[, l] * E$a[, m] + q0 * E$A[, l, m])
    list(S = S, Hf = Hf)
  }
  eta <- ebes; active <- rep(TRUE, nsub)
  for (it in seq_len(maxit + 1L)) {                      # it=1 evaluates at eta0 (skip test), then Newton steps
    Es <- .foceiAnalyticSolveAll(am, th, eta, ids, data, obsTimes, tol)
    if (is.null(Es)) return(NULL)
    for (i in which(active)) {
      sh <- .SHi(Es[[i]], eta[i, ], i)
      # a subject that cannot be solved gives a non-finite score; max(abs(S)) is then NA
      # and the test below would error rather than fall back
      if (!all(is.finite(sh$S))) return(NULL)
      if (max(abs(sh$S)) < (if (it == 1L) skip else conv)) { active[i] <- FALSE; next }
      if (it == maxit + 1L) return(NULL)                 # did not converge -> FD fallback
      step <- tryCatch(solve(sh$Hf, sh$S), error = function(e) NULL); if (is.null(step)) return(NULL)
      eta[i, ] <- eta[i, ] - step
    }
    if (!any(active)) break
  }
  if (any(active)) return(NULL)
  eta
}

#' Compute the full analytic FOCEI covariance (theta + sigma + Omega) for a fitted
#' object, or `NULL` when out of scope / the augmented solve fails.  The cached,
#' env-installing entry point is [foceiCovAnalytic]; this is the raw compute.
#' @param fit a fitted nlmixr2 focei object
#' @return list(cov, se, R, params, method) or `NULL`
#' @noRd
.foceiCovAnalyticCalc <- function(fit) {
  ui <- fit$finalUi
  if (!.hasRxSens())
    return(.foceiAnalyticFallback("an rxode2 without symbolic sensitivities"))
  # FO/FOI is out of scope (the analytic (f,R) path is a FOCEI/FOCE observed information).  The
  # runtime `fo` flag is not persisted to fit$finalUi, so this standalone entry also keys on the
  # persisted estimation method (ui$control$est) -- otherwise an FO fit would be silently assembled
  # as FOCE and mislabelled "analytic".  (The live covType="analytic" hook reads the in-fit ui
  # where `fo` is set, and FO forces covMethod=0, so the production path is already safe.)
  if (isTRUE(as.logical(rxode2::rxGetControl(ui, "fo", FALSE))) ||
      isTRUE(rxode2::rxGetControl(ui, "est", "") %in% c("fo", "foi")))
    return(.foceiAnalyticFallback("the FO/FOI method"))
  # linCmt() has no symbolic state sensitivities for the augmented model
  if (isTRUE(any(ui$predDf$linCmt)))
    return(.foceiAnalyticFallback("a linCmt() model"))
  if (.foceiCholSECovActive(ui))
    return(.foceiAnalyticFallback(.foceiCholSECovReason))
  interaction <- as.integer(rxode2::rxGetControl(ui, "interaction", 1L))                   # 1 FOCEI / 0 FOCE
  # FOCE variance mode (0 "nonmem" frozen R0, 1 "foce+" live R); FOCEI ignores it
  foceType <- if (interaction == 0L) as.integer(rxode2::rxGetControl(ui, "foceType", 0L)) else 0L
  # AGQ needs interaction=1 -- see .foceiCalcRanalytic for why this gate is load-bearing.
  if (as.integer(rxode2::rxGetControl(ui, "nAGQ", 1L)) > 1L && interaction != 1L)
    return(.foceiAnalyticFallback("adaptive Gaussian quadrature (nAGQ > 1) without interaction"))
  # IOV is out of scope; derive the flag from THIS fit, not the process-global
  # .uiIovEnv (which reflects the LAST-preprocessed model).  An occasion (IOV) eta
  # carries a non-"id" `condition` -- the same signal .uiApplyIov keys on.
  .idf0 <- ui$iniDf
  if (any(!is.na(.idf0$condition) & .idf0$condition != "id" & is.na(.idf0$err)))            # IOV
    return(.foceiAnalyticFallback("inter-occasion variability (IOV)"))
  # censored (M2/M3/M4): FOCEI and FOCE with censOption="gauss" are in scope (censored score
  # partials + Gauss-Newton determinant); only the laplace censored determinant uses FD.
  .hasCens <- (!is.null(fit$dataSav$CENS) && any(fit$dataSav$CENS != 0, na.rm = TRUE)) ||
    (!is.null(fit$dataSav$LIMIT) && any(is.finite(fit$dataSav$LIMIT)))
  if (.hasCens && as.integer(rxode2::rxGetControl(ui, "censOption", 0L)) == 1L)
    return(.foceiAnalyticFallback("censoring with censOption='laplace'"))
  ef <- .foceiAnalyticErrFull(ui)
  if (is.null(ef)) return(NULL)                     # unsupported error model -> errFull already messaged

  # AGQ scope -- see .foceiCalcRanalytic for why this must gate above the routing (IOV is
  # already declined outright on this path).
  if (as.integer(rxode2::rxGetControl(ui, "nAGQ", 1L)) > 1L) {
    .agqLo <- suppressWarnings(as.numeric(rxode2::rxGetControl(ui, "agqLow", -Inf)))
    .agqHi <- suppressWarnings(as.numeric(rxode2::rxGetControl(ui, "agqHi", Inf)))
    .why <- if (isTRUE(ef$foceiOnly)) "a general or multi-endpoint residual variance"
      else if (.hasCens) "censored observations"
      else if (isTRUE(is.finite(.agqLo)) || isTRUE(is.finite(.agqHi))) "a finite agqLow/agqHi node clamp"
      else NULL
    if (!is.null(.why))
      return(.foceiAnalyticFallback(paste0("adaptive Gaussian quadrature (nAGQ > 1) with ", .why)))
  }

  ini <- ui$iniDf
  .map <- .foceiEtaThetaMap(ui)                    # theta <-> eta pairing
  etaNames <- .map$etaNames
  neta <- length(etaNames)
  if (neta == 0L) return(.foceiAnalyticFallback("no random effects"))
  thetaForEta <- .map$thetaForEta
  # a non-mu-ref (orphan) eta is in scope: it keeps its own ETA_i_ sensitivity direction
  # and an eta-named Omega variance (matching the production hook); no theta maps to it.
  if (any(.iniIsFixed(ini, thetaForEta)))          # fixed structural theta breaks eta indexing
    return(.foceiAnalyticFallback("a fixed mu-referenced structural parameter"))
  keep <- !.iniIsFixed(ini, ef$sgName); ef$sgVar <- ef$sgVar[keep]; ef$sgName <- ef$sgName[keep]  # drop fixed sigma
  Om <- fit$omega
  pairs <- .foceiOmegaPairs(Om, ini)               # free Omega lower-triangle (declared blocks)
  omd <- .omegaVarCovDeriv(Om, pairs)

  # Uniform direction assembly (shared with the production covType="analytic" hook).
  .dir <- .foceiAnalyticDirections(ini, thetaForEta, ef$sgName, neta,
                                   sharedEta = unname(.foceiEtaOccurrence(ui) > 1L))
  if (is.null(.dir)) return(.foceiAnalyticFallback("an unresolvable direction set"))
  thStruct <- .dir$thStruct; dirs <- .dir$dirs; dirTh <- .dir$dirTh
  ndir <- .dir$ndir; nth <- .dir$nth

  thRows <- ini[!is.na(ini$ntheta), , drop = FALSE]
  thRows <- thRows[order(thRows$ntheta), , drop = FALSE]
  # converged estimates (fit$theta), not ui$iniDf$est which can hold initials
  .thv <- fit$theta[thRows$name]
  if (anyNA(.thv)) .thv <- thRows$est
  th <- setNames(as.numeric(.thv), paste0("THETA_", seq_len(nrow(thRows)), "_"))
  ebes <- as.matrix(fit$eta[, etaNames, drop = FALSE])
  nsg <- length(ef$sgVar)

  # FOCEI with a general (non-add/prop, multi-endpoint, or estimated-lambda) variance uses the
  # (f,R) cov (sigmas as directions), matching the live covType="analytic" hook; add/prop keeps
  # the fast symbolic assembly.  (IOV already bowed out above.)  An estimated boxCox/yeoJohnson
  # lambda sits in thStruct as a theta-like direction, so name the sigma block with the
  # lambda-excluded .dir$sgName (matching the live hook's fullNm).  Censored FOCEI also routes
  # here (the fast add/prop assembler has no censored partials).
  R <- if (isTRUE(ef$foceiOnly) || .hasCens)
    .foceiAnalyticAssembleRFR(ui, th, ebes, fit$eta$ID, fit$dataSav, Om, ef, neta,
                              length(.dir$dirP), .dir$dirP, omd,
                              dirsCov = .dir$dirsCov, ndirCov = .dir$ndirCov,
                              solveTol = .foceiAnalyticSolveTol(ui),
                              interaction = interaction, foceType = foceType, lamDir = .dir$lamDir)
  else .foceiAnalyticAssembleR(ui, th, ebes, fit$eta$ID, fit$dataSav, Om, ef, neta, nth, nsg, omd,
                               dirs = dirs, dirTh = dirTh, ndir = ndir,
                               solveTol = .foceiAnalyticSolveTol(ui), interaction = interaction,
                               foceType = foceType)
  if (is.null(R)) return(.foceiAnalyticFallback("an observed information that would not assemble"))
  cov <- tryCatch(solve(R), error = function(e) NULL)
  if (is.null(cov)) return(.foceiAnalyticFallback("an observed information that would not invert"))
  onm <- etaNames                                            # Omega named by the eta (om.eta.cl)
  nm <- c(thStruct, .dir$sgName, .foceiOmegaCovNames(pairs, onm))
  dimnames(R) <- dimnames(cov) <- list(nm, nm)
  list(cov = cov, se = setNames(suppressWarnings(sqrt(diag(cov))), nm),  # NaN flags non-PD
       R = R, params = nm, method = "analytic")
}

#' Full analytic FOCEI covariance (theta + sigma + Omega) for a fitted object, or
#' `NULL` when out of scope / the augmented solve fails.
#'
#' The result is cached on the fit environment (`.covAnalytic`) and its `$cov` is
#' installed as the fit's `$cov` on the first call, so repeated calls -- and
#' `getVarCov()` -- return the stored covariance instead of recomputing the
#' augmented sensitivity solve every time.
#' @param fit a fitted nlmixr2 focei object
#' @return list(cov, se, R, params, method) or `NULL`
#' @noRd
foceiCovAnalytic <- function(fit) {
  .env <- fit
  if (rxode2::rxIs(fit, "nlmixr2FitData")) .env <- fit$env
  if (exists(".covAnalytic", envir = .env, inherits = FALSE)) {
    return(get(".covAnalytic", envir = .env))
  }
  # Match the live covType="analytic" hook (.foceiCalcRanalytic), which wraps the whole assembly
  # in tryCatch and returns NULL on any error -> FD fallback.  A direct foceiCovAnalytic()/
  # getVarCov() call must fall back just as gracefully (e.g. a pure-proportional FOCE fit whose
  # near-zero-prediction branch can hit an NA), never throw.
  .ret <- tryCatch(.foceiCovAnalyticCalc(fit), error = .foceiAnalyticErrWarn(2L))
  assign(".covAnalytic", .ret, envir = .env)   # cache (incl. NULL) -- do not recompute
  if (!is.null(.ret) && is.matrix(.ret$cov)) {
    .env$cov <- .ret$cov                        # install so getVarCov()/$cov reuse it
    .env$covMethod <- "analytic"                # report the analytic observed information
  }
  .ret
}

Try the nlmixr2est package in your browser

Any scripts or data that you put into this service are public.

nlmixr2est documentation built on Aug. 5, 2026, 1:11 a.m.