rxSolve(params = <multi-row data.frame>, omega = NA) reads etas out of bounds -- silently wrong, non-deterministic results with >= 8 subjects

[billdenney]

Summary

When rxSolve() is given a multi-row params data.frame (one parameter set per
id) together with omega = NA, it returns different results for
byte-identical inputs. Some of those results are silently wrong, and some are
non-finite. omega = NA is documented to set the between-subject etas to zero;
instead the etas are filled with whatever happens to sit next to a length-1
vector on R's heap.

This is a silent-corruption bug: no warning, no error, and most values in each
solve are correct, so it is easy to miss in a simulation.

Confirmed on current main (v5.1.7, 7d74e29). The offending line was
introduced on 2022-04-01 in cd69147 ("pred-only sim") and first released in
v2.0.6, so this is long-standing rather than a regression.

Reproducible example

Deliberately trivial so the correct answer is known exactly, with no solver
tolerance involved: base <- tbase + eta.base, and with omega = NA the eta
must be zero, so base must be exactly 1..16 on every solve.

library(rxode2)

mod <- function() {
  ini({
    tbase <- 1
    eta.base ~ 0.1
    addSd <- 1
  })
  model({
    base <- tbase + eta.base
    base ~ add(addSd)
  })
}

n <- 16L
pars <- data.frame(id = seq_len(n), tbase = as.numeric(seq_len(n)))
ev <- data.frame(id = seq_len(n), time = 0, evid = 0L, amt = 0)

# `omega = NA` is documented to set the between-subject etas to zero, so
# `base` must equal `tbase` -- exactly 1..16 -- on every one of these solves.
for (i in 1:10) {
  print(rxSolve(mod, ev, params = pars, omega = NA, returnType = "data.frame")$base)
}

Observed (rxode2 5.1.7)

 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16
 [1]  1  2  3  4  5  6  7 10  9 10 11 12 13 14 15 16     <- element 8 wrong
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 16 16     <- element 15 wrong
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 14 16     <- element 15 wrong
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16

Expected

All ten lines identical and exactly 1 2 ... 16.

Which runs go wrong, and which elements within a run, changes from run to run
and from session to session. At 200 subjects it is far worse -- in one measured
batch, 50/100 solves were wrong, 12/100 contained non-finite values, and the
largest deviation from the correct value was 1.06e+11.

Root cause

R/rxsolve.R:2073-2079,
in the omega = NA branch of rxSolve.rxUi:

} else if (is.logical(.rxControl$omega)) {
    if (is.na(.rxControl$omega)) {
      .omega <- object$omega
      params <- c(params, setNames(rep(0, dim(.omega)[1]), dimnames(.omega)[[2]]))
      .rxControl$omega <- NULL
    }
  }

c() on a data.frame drops the data.frame class and returns a plain ragged
list: the per-id parameter columns keep length n_id, while the appended eta
zeros have length 1.

pars <- data.frame(id = 1:4, lbase = c(.1,.2,.3,.4), lic50 = c(-1,-2,-3,-4))
lengths(c(pars, setNames(rep(0, 2), c("etalbase", "etalic50"))))
#>       id    lbase    lic50 etalbase etalic50
#>        4        4        4        1        1

Downstream, the solver reads n_id values out of each length-1 eta vector, i.e.
out of bounds. I isolated this by bypassing omega = NA entirely and handing
rxSolve() the list that line 2076 builds, versus the same list with a
full-length eta:

params handed to rxSolve()	eta element	result (50 solves, 16 ids)
ragged list, as built by line 2076	length 1	1/50 wrong, 2 distinct
same list, full-length eta	length 16	0/50 wrong, 1 distinct
data.frame with an eta column	length 16	0/50 wrong, 1 distinct

The only thing that differs between the broken and working cases is the length
of the eta element -- the second case is still a plain list, and it is fine. So
the list class is not what matters; the length-1 vector being read n_id times
is.

The >= 8 threshold is consistent with R's small-vector size classes: a length-1
REALSXP is allocated in a node with room for several doubles, so reads just
past the end stay inside that node (and find zeros); further reads cross into
adjacent live heap.

Note that 25 lines earlier, the same function already branches correctly on
data.frame-vs-numeric params
(R/rxsolve.R:2049-2061),
using params[[.t]] <- .theta[.t] for the data.frame case and c() only for the
numeric case. The omega = NA branch just uses c() unconditionally.

Boundary map

Measured on 5.1.7 with the reprex model, 100 solves per cell, counting solves
whose output differs from the known-correct 1..n. Rates are stochastic (they
depend on heap layout), so a 0/100 near the threshold is not proof of safety --
but <= 7 subjects was solid across every repetition I ran.

subjects	cores = 16	cores = 1
2	0/100	0/100
4	0/100	0/100
6	0/100	0/100
7	0/100	0/100
8	0/100	2/100
9	1/100	0/100
16	4/100	5/100
64	19/100	10/100
200	50/100	40/100

It fails identically at cores = 1, so this is not a threading race.

Ruled out

• R's RNG -- the per-id parameter draws are generated once, before the solve
  loop, and hash identically across runs. The reprex above uses no random draws
  at all.
• Threading -- reproduces at cores = 1.
• The residual error family -- reproduces with add(); originally found with
  logitNorm().
• ODE integration -- the reprex model is purely algebraic, with no d/dt()
  states.
• Number of observation rows -- reproduces with a single observation per
  subject.
• The specific model -- reproduces with a one-parameter, one-eta model.

Suggested fix

Make the omega = NA branch data.frame-aware. The codebase already does exactly
this for the same "treat these omega/sigma items as zero" operation at
R/rxsolve.R:3148-3155:

if (inherits(params, "data.frame")) {
  for (v in .ctl$.zeros) {
    params[[v]] <- 0.0
  }
} else if (inherits(params, "numeric") ||
             inherits(params, "integer")) {
  params <- c(params, setNames(rep(0.0, length(.ctl$.zeros)), .ctl$.zeros))
}

Applying that same shape at line 2076 fixes every case above -- 0/100 wrong at
every subject count, at both cores = 1 and cores = 16, and the original
200-subject model I first hit this with becomes deterministic:

.etaZero <- setNames(rep(0, dim(.omega)[1]), dimnames(.omega)[[2]])
if (inherits(params, "data.frame")) {
  for (.e in names(.etaZero)) params[[.e]] <- 0.0
} else {
  params <- c(params, .etaZero)
}

(Offered as a diagnosis rather than a polished patch -- I have not checked it
against the full test suite, and you may prefer to factor the two sites into one
helper.)

Two related things at the same code site

1. sigma = NA has the identical unguarded c() at
   R/rxsolve.R:2109.
   It did not silently corrupt in my testing -- with a multi-row params
   data.frame it errors with The following parameter(s) are required for solving: rxerr.base, addSd, tbase -- but it looks like the same latent
   defect and would be worth fixing together.

2. omega = NA on a model with no etas gives a confusing error. object$omega
   is NULL, so dim(NULL)[1] is NA and rep(0, NA) raises
   invalid 'times' argument from inside rxSolve. Arguably omega = NA
   should be a no-op there, or say what is actually wrong.

For contrast, passing params as a matrix fails loudly and clearly
(The following parameter(s) are required for solving: ...) both before and
after the fix, so that path is not affected by this.

Environment

rxode2      5.1.7 (origin/main @ 7d74e29e4, built from source)
rxode2ll    2.0.16
R           4.6.1 (2026-06-24), x86_64-pc-linux-gnu
OS          Ubuntu 24.04.4 LTS
rxCores()   16

How this surfaced

An intermittently-failing pkgdown build in
nlmixr2/nlmixr2lib#488 -- a vignette that solved an
uncertainty band over 200 parameter draws with
rxSolve(..., params = <data.frame>, omega = NA) and got a different band, and
occasionally Inf/NaN, on each render. Worked around there by evaluating the
band in closed form; the workaround stays until this is fixed upstream.