Fit a grouped call's levels under the plan, and stop the formula exporting the session (#338, #329)

[beckyfisher]

What

Two changes to what reaches a future worker, which land together because the
first multiplies the second.

bnec_group() fitted its levels one after another whatever plan was set. It now
dispatches them through the dispatcher the model loop already uses, and decides
which of the two loops the plan drives (#338).

A formula records the environment it was created in, and serialising it writes
that environment out in full. Every fit therefore sent the calling session to
every worker, and kept it in the saved fit. The formula's environment is now
rebuilt to hold only the names the formula mentions and the data does not
supply (#329).

A grouped call also gives different numbers than it did. Every level is now
given a seed realised in the calling session, so that the arrangement cannot
change the estimates; the reason is under Reproducibility below.

Why it matters

A grouped call multiplies the model set by the number of levels, so it is where
a model-averaged analysis stops being practical. Seven herbicides at eleven
equations is 77 fits.

The export made the parallel path unusable from any session that already holds a
fit, which is what precompiling a vignette does: it failed at
future.globals.maxSize rather than merely slowing. A grouped call would have
sent that environment once per level as well as once per model, so #338 without
#329 multiplies a failure. Two further consequences of the same fix: a saved fit
no longer stores the session it was fitted in, and amend() no longer exports
it either, by way of the .Environment of the model frame's terms attribute.

One premise of #338 does not survive measurement, and it changes what the level
loop buys. One plan drives one loop, because future evaluates a nested future
sequentially unless the plan is a list, so dispatching the levels takes the
workers away from the model loop rather than adding to them. Counting the two in
rounds of one fit, the levels take ceiling(L / W) * M and the models
L * ceiling(M / W). The level arrangement cannot go below M rounds however
many workers there are, and the model arrangement reaches L, so with fewer
levels than equations the levels are the smaller count only in a band. Seven
levels of eleven equations at the four cores #338 measured is 22 rounds against
21, so that call is left exactly as it is; at eight workers it is 11 against 14,
and the levels take the plan. bnec_group() makes that comparison and reports
it. A nested plan drives both loops and is honoured without being counted, and
plan(list(sequential, ...)) asks for the release arrangement explicitly.

Reproducibility

A grouped call's estimates no longer depend on how it was arranged, and this
needed making true rather than being inherited. bnec() realises the seed for
its weighted posterior draw from whatever stream it is running in, so a level
fitted in a worker drew from that worker's stream and a level fitted in the
parent from the session's. Since the arrangement is read off the worker count,
one script at one seed would have reported different model-averaged estimates on
a four-core and an eight-core machine, which nothing in the released package
does.

Every level is therefore given a seed realised in the calling session before the
levels are dispatched, derived from seed where one was passed and from the
session's stream otherwise, and the sampler kind is pinned so that the same
seed means one thing in a session set to the pre-3.6.0 sampler. The realised
seeds are stored on the returned object. The calling session's own stream and
generator kind are restored afterwards, on both arrangements.

Two consequences. A grouped call now repeats under seed alone with no
set.seed() in the session, which is more than bnec() offers for a single set.
And a grouped call reports different numbers than earlier versions did, because
the levels are seeded rather than drawing from the stream as each is reached;
they are another realisation of the same weighting.

Evidence

Nested futures

R 4.6.1, future 1.70.0. Inside a worker of plan(multisession, workers = 2)
the strategy is sequential, nbrOfWorkers() is 1, and tempdir() differs
between workers. Inside a level worker of plan(list(tweak(multisession, workers = 2), tweak(multisession, workers = I(3)))), nbrOfWorkers() is 3 and
the model loop parallelises over them. Under plan(list(sequential, tweak(multisession, workers = I(3)))) the parent reports one worker and the
level future reports 3, which is what makes that plan the release arrangement
rather than a fully sequential run.

The I() is needed rather than decorative: future pins mc.cores to 1 inside
a worker, parallelly reads that as the core budget, and a bare workers = 4
inner count returns Attempting to set up 4 localhost parallel workers with only
1 CPU cores available ... The hard limit is set to 300% where I(4) returns 4.

The size of the export

R 4.6.1, a 76 MiB vector bound beside the formula in the environment that
created it, measured with serialize():

object	released	here
the bayesnecformula	76.29 MiB	under 0.01 MiB
the closure narrow_environment() builds around it	76.30 MiB	under 0.01 MiB
the model frame built from it	76.30 MiB	0.006 MiB

The model frame is the route amend() exports, through the .Environment of
its terms attribute. A formula whose environment is the global environment is
left alone, because R serialises that one by reference.

Separately, on a 200,000-row frame over seven levels: what each level future
is sent fell from 5.53 MiB per future, 38.71 MiB over the seven, to 5.53 MiB for
the seven subsets together, because the subsets are now the dispatched elements
rather than indices into a data frame held in the exported environment.

Dispatch packing

R 4.6.1, plan(multisession, workers = 4), a body that sleeps two seconds
standing in for one fit, after one call to pay the namespace load: four elements
took 2.2 s under future_lapply() directly and 4.6 s under
bnec_parallel_lapply() on its first call; seven elements took 4.1 s, which is
the two rounds group_loop_rounds() counts and not the one a chunk-per-worker
dispatch would have given.

Wall clock

No wall-clock measurement of a real grouped call was taken. Three attempts were
made and none completed: the development host was running other work throughout,
and repeated multisession cluster creation in one session stalls on it often
enough that two sleep-based harnesses and one four-worker grouped fit all
exceeded twenty minutes, while the same dispatcher measured on its own returned
in seconds. The arrangement claim therefore rests on the arithmetic, which is
unit tested, on the packing measurement above, and on the future semantics
measured at the top of this section. The shapes where the level loop is the
smaller count need eight or more workers, so the measurement that would settle
it belongs on the HPC and not here.

Tests

tests/testthat/test-parallel_models.R gains the round arithmetic including
both shapes #338 quotes and the unbounded-worker case, the five decisions
plan_group_levels() can take, the compile directories asserted through a real
future, four properties of the formula narrowing including the case where it
saves nothing and the NULL binding that used to be dropped, the level seeding
and its sampler pinning, that a seeded body gives one answer on both
arrangements, that a grouped call leaves the caller's stream where it found it,
and one grouped fit, with no plan set, asserting the per-level rows, the
ordering and naming of the returned fits, the recorded level seeds and the
narrowed formula.

Nothing in the suite fits under a plan, deliberately. Cross-arrangement equality
is a property of the per-level seeding and is pinned without fitting; what a
grouped call under a plan would add is a Stan compilation per worker. It is also
unreliable on at least one machine: on a WSL2 host a single bnec() fit inside a
multisession worker did not return in 300 s, with the workers asleep in poll
having used 9 s of processor time, while the same fit with no plan returns in
173 s and bnec_parallel_lapply() over a body that loads brms, raises, or
calls bnec() without fitting returns in under 3 s. That path predates this
branch --- it is the model loop #184 added --- and the interaction was not
identified.

Implementation detail

#338, the level loop

R/bnec_group.R. The for loop is replaced by bnec_parallel_lapply(), which
is bnec_model_lapply() renamed: it serves both loops now and its guarantees --
one element per chunk, future.seed = TRUE, the RNG kind restored inside the
worker, the caller's stream restored after -- are properties of the dispatch
rather than of what is dispatched.

plan_group_levels() in R/parallel_models.R is the counterpart of
plan_model_set() and reads the decision from the plan for the same reason: an
argument beside it would duplicate state the plan already holds. It returns
dispatch and concurrent separately. dispatch says whether the levels go
through future at all; concurrent says whether two of them can be fitted at
one time, and it is what decides whether the compile cache has to be separated
and whether a per-level message can still be printed in order. A list plan whose
first strategy is sequential has dispatch without concurrent, which is the
release arrangement asked for explicitly -- entering a future is what advances
the plan to the inner strategy, so the model loop still gets all of it. A list
plan with nothing parallel after the first strategy is not read as a request at
all.

group_loop_rounds() is the pure function the decision rests on and is unit
tested on its own. Ties go to the models, which is what the release does; an
unbounded worker count is counted as one round of whichever loop is dispatched,
and NA is reserved for a count that genuinely could not be read.

The comparison assumes every fit takes the same time, and that is wrong in a
known direction: the 23 equations differ in fitting time by an order of
magnitude, a round is set by its slowest member, and the level arrangement puts
M fits in a round against the model arrangement's W. Where M is the larger
the count therefore favours the levels more than the clock would, which is the
second reason ties go to the models. n_models is the set the formula asks for,
not the set that will be fitted: check_models() runs inside each level's
bnec() call and may drop an equation, so this is a scheduling estimate and the
message says so.

The compile cache

Every level fits the same equations, so parallel levels
build the same Stan programs at once, and cmdstanr does not lock its cache --
the hazard hpc/precompile-hpc.sh:255 serialises its job array to avoid.
worker_stan_cache_dir() gives each worker process a directory of its own,
under cmdstanr_write_stan_file_dir where the user has set one and tempdir()
otherwise. Keyed on the process rather than on the level: a level key would make
a worker that draws three levels compile every equation three times, which for
the #338 shape is 77 compilations against at most one set per worker. It is
applied on both backends although only multicore shares a tempdir(), because
a directory the user chose is shared by every worker whatever the backend. What
this adds is compilation on a cold cache, and group_loop_rounds() counts
sampling only, so on short fits the arrangement the count prefers can be the
slower one in wall clock; vignette("example2") says so under The levels of a
grouped call, and says that the per-process directories accumulate under a root
the user set and are safe to delete. The worked treatment of the level loop
lives there rather than in ?bnec_group, which keeps the rule, the
reproducibility guarantee and level_seeds and defers the rest, as ?bnec
defers to the same vignette. The
root is read in the parent and passed into the worker: future exports globals
and not options, so read inside the worker a persistent cache would have been
honoured under a forking plan and silently ignored under every other one.
rstan under rstan_options(auto_write = TRUE) caches in one place a forked
worker shares with its parent and offers no argument that changes it, so that is
documented rather than worked around.

Peak memory

Peak memory is reported rather than capped. A worker fitting a level holds
that level's whole model-averaged set, so W level workers hold up to W of
them where fitting in sequence holds one. The count above already declines to
dispatch the levels at the small worker counts where memory is tightest, which
is the case #190 hits.

Argument forwarding

do.call("bnec", ...) replaces the forwarding of
..., because ... cannot be put in the environment narrow_environment()
builds. By name rather than by value, so that an error raised inside a level
still names bnec(). Two consequences were checked rather than assumed. The
formula reaches bnec() with its environment intact, which #319 needs for
crf() to resolve a symbol -- pinned by a test. And the family arrives as a
value rather than as the symbol it was, which family_link_source() already
treats as "symbol"; R/validate_family.R names do.call() as that case, and
the family has in any event been validated and marked by then.

The elements dispatched

The dispatched elements are the subsets
themselves, built once with split(data, grp), so neither the data frame nor
the grouping vector is in the exported environment. split() orders its pieces
by levels(grp), which is the order levs is in, and an empty level cannot
arise because a level under four observations is refused earlier.

#329, what the formula carries

narrow_formula_environment() in R/parallel_models.R. It walks up from
environment(formula), copying each name the formula mentions and data does
not supply, and stops at the first environment R serialises by reference --
tested with environmentName(), which is non-empty for exactly those. The
replacement is parented at the environment the walk stopped at, not at the
package namespace, so the lookup chain is the one the formula had: a formula
written in a frame whose parent is baseenv() still resolves through
baseenv(), which is what the trials() test added by #334 rests on.

Blanking the environment would have been wrong: model.frame() resolves a term
against data first and the formula's environment second, so a formula naming
anything outside data needs it, and that includes the model argument of
crf() where it is a variable. Since #319 merged into dev that resolution is
deliberate rather than incidental -- formula_eval_env() evaluates in a frame
whose parent is formula_env(). This branch was rebased onto dev after that
merge and the formula test files were run against it. A test asserts that a
mods variable holding the model set survives, and another that a name bound to
NULL is kept rather than dropped -- vals[[nm]] <- NULL deletes, which left
such a name to resolve against the parent chain, where a decoy of the same name
answered for it.

Applied in bnec(), amend() and bnec_group() before the model frame is
built, so that the frame, the brms formula and the stored fit are narrowed by
one call, and on the sequential path as well as the parallel one, for the reason
narrow_environment() gives: a name left behind then fails on the first
ordinary call rather than only for whoever sets a plan.

Two limits, both measured and both documented:

• A name bound in the global environment, an attached package or a namespace is
  left where it is and resolved through the parent chain, because copying it
  would reintroduce the size this removes. A worker's global environment is not
  the caller's, so such a name was already out of reach under a plan and still
  is.
• A function defined beside the formula is a closure over that same environment,
  so the name is copied with everything it closed over. Measured at 76.294 MiB
  before and after: unchanged, and there is nothing to be saved, because the
  object is needed. A helper written in a knitr chunk and used in a formula
  therefore still gives a large fit.

One pre-existing behaviour was found and left alone: model.frame() on a
bayesnecformula does not resolve a user function named inside crf() from the
formula's environment at all, failing with "could not find function" both before
and after this change.

Not done here

• The two loops are still nested rather than flattened. Dispatching all
  L x M fits as one set would pack the workers better than any nesting and
  needs no comparison at all, but it means separating bnec()'s fitting loop
  from its preparation and its assembly, which is a change to the package's
  central function and is not this PR.
• expand_manec() still realises w_draw_seed from the ambient stream. The
  level seeds make that harmless for a grouped call; deriving it from
  brm_args$seed would close it for bnec() as well, and is the change Model-averaged output is not reproducible: predict() on a bayesmanecfit differs between calls #216
  decided against.
• brm(file = ) passed through ... names one file for every model of a set
  and every level of a group, so the second fit reads back the first. That is a
  defect of bnec() rather than of the level loop, and predates it.

🤖 Generated with Claude Code

https://claude.ai/code/session_01B7T9mkqxuRAKcdXmN5pctD

[beckyfisher]

The change was reviewed over four rounds by an independent reviewer working from
the branch rather than the diff. Every finding but one is fixed; that one and one
other proposed remedy were declined, with reasons below. The body above describes
the current state of the change, which is not the state first pushed.

The findings that mattered

Wrong answers

A grouped call's model-averaged estimates
depended on the worker count, because expand_manec() realises the seed for its
weighted posterior draw from whatever stream it is running in and the arrangement
is chosen from that count: one script at one seed reported different estimates on
a four-core and an eight-core machine. And the formula narrowing deleted a name
bound to NULL rather than keeping it, so the name resolved against the parent
chain instead — a mods <- NULL beside the formula, with an unrelated mods in
the global environment, turned a refusal from get_model_from_formula() into a
fit of something else.

Wasted transfer, wasted compilation, a disturbed RNG stream

Each level future was sent the whole data frame: 38.71 MiB over seven levels
against the 5.53 MiB they need between them. A compile directory per level made a worker that
drew three levels compile every equation three times, 77 compilations against at
most one set per worker for the shape #338 measures. And a grouped call left the
caller's stream at the last level's seed, and at a constant where seed was
supplied, so set.seed(i); bnec_group(..., seed = 1); rnorm(1) returned one
number for every i.

Documented calls that did not do what they said

The nested plan
?bnec_group printed errored, because future pins mc.cores to 1 inside a
worker and parallelly refuses four workers against one core; it now uses
workers = I(4). plan(list(sequential, ...)), the natural way to ask for the
arrangement earlier versions used, fitted both loops in sequence and said
nothing; it is now that request. And nbrOfWorkers() reports Inf on a
scheduler backend, which was read as a count that could not be read and declined
— the setting where the level loop pays best.

Every finding, and what was done with it

First round

1. The nested-plan recipe in ?bnec_group, NEWS.md and the message errored.
   Fixed with workers = I(4) and one sentence saying why it is needed.
2. expand_manec()'s weighted draw was realised in a worker under the level
   arrangement, so the worker count changed the model-averaged estimates. Closed
   rather than documented: group_level_seeds() realises one seed per level in
   the calling session, derived from seed where one was passed.
3. Each level future was sent the whole data frame. The dispatched elements are
   now the subsets, built once with split(data, grp).
4. A compile directory per level forced M x L compilations, and two claims
   about cache reuse were false. Keyed on the process instead; the false claims
   are withdrawn and what the arrangement adds on a cold cache is now stated.
5. There was no way to ask for the release arrangement, and
   plan(list(sequential, ...)) silently fitted both loops in sequence.
   plan_group_levels() now returns dispatch and concurrent separately, and
   that plan is the explicit request.
6. An unbounded worker count was treated as unknown. Counted as one round of
   whichever loop is dispatched; NA is reserved for a count that could not be
   read.
7. do.call(bnec, ...) recorded the function object, so every error raised in a
   level read Error in (function (formula, data, .... do.call("bnec", ...).
8. vals[[nm]] <- NULL deletes; a name bound to NULL was left to resolve
   against the parent chain. vals[nm] <- list(...), with a test using a decoy
   of the same name.
9. The end-to-end grouped test would have passed with the level loop removed. It
   now asserts the dispatch message from the real call and that the compile
   option is restored in the parent.
10. The message said "models" where the repository says "equations", and stated
    the equation count as fact where it is the set the formula asks for before
    any is ruled out.
11. House style: one idiom, and two metaphors standing in for the mechanism.

Second round, on the revision

12. The per-level set.seed() left the caller's stream pinned on the sequential
    path. The stream and the generator kind are now restored around the whole
    loop, including on an error exit.
13. On the release arrangement every per-level message arrived at the end,
    because a sequential future relays its conditions at collection. The levels
    are now dispatched one at a time and announced before each starts, which
    still reaches the inner plan.
14. The realised level seeds were neither recorded nor stable across
    sample.kind. They are stored on the returned object, and the sampler is
    pinned.
15. plan(list(sequential, sequential)) was accepted and promised parallelism a
    sequential inner strategy does not provide.
16. n_workers(Inf) printed "Inf workers".
17. dots$seed partial-matches; dots[["seed"]].

Third round

18. level_seeds was on an exported object and documented nowhere; the class
    documentation was also missing weights_method. Both are @slots now, and
    ?bnec_group says the realisation is kept and why.
19. Only the sampler was pinned, not the generator, so the same seed realised
    different level seeds in a session set to L'Ecuyer-CMRG — which is what a
    user doing their own parallel work sets.
20. The concurrent message named the next strategy in a list plan where that
    strategy is sequential.

Raised and not done

A test that builds the recipe's nested cluster

The reviewer measured that
plan at 1.7 s in a fresh process. The same test inside this file did not complete
in 280 s on two runs, and repeated multisession cluster creation in one session
stalls on the development host often enough that it bit a pre-existing test three
times today. Both measurements are true; the second describes the suite. The test
asserts instead that the advice the package prints keeps the I(), with both
measurements recorded in the comment. The stall itself is worth its own issue,
with connection exhaustion as the hypothesis to start from.

Deleting the per-worker compile directories after use

That gives back the
whole benefit of keying on the process for any worker that draws a second level,
and it has a modelling function removing files under a directory the user
nominated. ?bnec_group says instead that the directories accumulate, are never
read again, and are safe to delete.

[beckyfisher]

CI runs were failing repeatedly following the Claude solution. Codex sol (high) fixed this.

Changes:

• Replaced the unreliable model-frame size ratio with direct checks that big is unreachable through both the narrowed formula and the model frame’s
  terms environment.

• Added a planned bnec_group() test using two forked workers and a mocked bnec().

• The new test covers concurrent level dispatch, subsets, ordering, worker processes, compile-cache separation and formula narrowing without additional
  Stan compilation.

• The existing real sequential grouped fit remains as the fitting contract.

Verification:

• Exact edited tests: 18 assertions passed in about 8 seconds.

• R parse and git diff --check: passed.

• A full test-file run was stopped after ten minutes inside the unchanged real Stan fit; no failure was emitted. This confirms that restoring the former
  four-real-fit test would substantially increase runtime.

• The new fork test runs on macOS/Linux; Windows retains the existing socket-dispatch coverage.