Skip to content

Add run_rceattle(): Rceattle as a sixth estimation model - #16

Open
grantdadams wants to merge 3 commits into
update-emsfrom
update-ems-w-Rceattle
Open

Add run_rceattle(): Rceattle as a sixth estimation model#16
grantdadams wants to merge 3 commits into
update-emsfrom
update-ems-w-Rceattle

Conversation

@grantdadams

Copy link
Copy Markdown

Adds run_rceattle() so Rceattle can join the model comparison as a sixth estimation model. Also adds a standalone smoke test.

What it adds

  • R/run_rceattle.R — reads each OM replicate, converts it to an Rceattle data_list, and fits three recruitment scenarios, writing one set of RDS files per replicate and scenario.
  • R/run_em.R — one dispatch line.
  • NAMESPACE, man/run_rceattle.Rd — the man page is generated from roxygen, so document() regenerates it byte-identically.
  • FIMS_integration_test_data/Rceattle_integration_test.R — end-to-end smoke test.

Scenarios: random_effects, random_effects_sigmaR_constant, fixed_effects. The figures use the first and last, matching how WHAM and FIMS each contribute two series.

Estimation-model configuration

All three fit Beverton-Holt (R = alpha*S / (1 + beta*S)) with alpha fixed at the value implied by the OM's steepness, and year-1 numbers-at-age estimated as free parameters.

Both choices follow the existing EMs. WHAM fixes log_SR_a from OM steepness and estimates 12 free initial numbers-at-age; FIMS estimates log_init_naa and holds steepness at its default 0.75. Alpha and beta cannot be estimated jointly on these data — the Hessian is singular — so every platform here fixes one of them.

Fits are hindcast-only, so no reference points, which Bai confirmed are not compared.

Non-regression: run_fims.R is untouched

run_fims.R already defines cv_2_sd, count_na_standard_errors(sdreport), and get_condition_number(obj, opt, sdreport) at package scope. DESCRIPTION has no Collate:, so R collates alphabetically and run_rceattle.R would win — a naive port would silently redefine get_condition_number with a different arity and break run_fims() at runtime.

Handled by sharing cv_2_sd (identical bodies; sharing it keeps observation error specified the same way across EMs) and prefixing the other two as rceattle_*. Each rename carries a comment explaining why, so it does not get tidied away later.

R/run_fims.R has a zero-line diff in this PR.

Verification

Smoke test on the bundled FIMS_C1 OM, against OM truth: SSB correlation 0.9999, biomass 0.9999, F 0.9986, recruitment 0.985, survey 0.9993, max gradient 6.4e-08, positive-definite Hessian, sigmaR 0.358 against a true 0.4. Catchability estimates 3.258e-07 against an OM truth of 3.315e-07.

The comparison-script readers were also run end-to-end against this output: all five metrics present, year counts matching the OM, no NA values.

A real run_fims() run is not possible on the bundled OMs (see below), so FIMS non-regression rests on the zero-line diff plus arity assertions.

Worth flagging

  1. FIMS fits the biomass survey; every other EM fits numbers. run_fims.R uses surveyB.obs in mt, while the other EMs fit survey.obs and read_output_data() compares against the OM's survey_q. The two series differ in scale and shape. It also means FIMS cannot run on the bundled FIMS_integration_test_data OMs at all — none of them contain surveyB.obs.
  2. Adding a sixth EM changes the other five's numbers. check_convergence() returns the intersection of converged simulations across all models, and the global outlier filter is shared, so the figures for every existing model shift. Worth knowing before a 500-replicate run.
  3. Rceattle uses analytical (closed-form MLE) catchability; FIMS estimates log_q freely. Comparable, but not an identical treatment. We can change that in Rceattle.
  4. On a deterministic OM (logR_sd == 0) the two random-effects scenarios are ill-posed — the Laplace inner problem for a zero-variance random effect is singular — and will hit time_limit and record as non-converged. Use a stochastic case.
  5. DESCRIPTION: Depends: R (>= 3.5.0) is already inaccurate on this branch, since run_fims.R uses the native pipe throughout. Not changed here.

The output-file manifest was hand-written into man/run_rceattle.Rd, so the
next devtools::document() would have silently dropped it. Moved into the
roxygen @return; the man page now regenerates byte-identically.

run_fims.R never sets logit_steep -- FIMS holds steepness constant at its
own package default of 0.75, which matches these operating models only by
coincidence. Say that rather than claiming FIMS pins it at the OM value,
since it is the justification for fixing alpha here.

Drop the uppercase output ignore rule; nothing writes that directory.
Three faults found in adversarial review, the first two demonstrated by
execution:

  * The socket cluster used the default outfile, which is /dev/null, so every
    message() and warning() raised inside the foreach body was discarded.
    run_rceattle(survey_units = "biomass") on an OM without surveyB.obs ran to
    completion, printed nothing and wrote no files -- the explanatory stop()
    never reached anyone. The cluster now keeps worker output on the console,
    and each replicate returns a status that the master reports as a warning.
  * The replicate-level tryCatch wrapped the whole scenario loop, but only
    fit_rceattle_scenario() handled its own errors. A stop() in
    rceattle_estimates() or a failed saveRDS in the first scenario therefore
    abandoned the remaining two. Each scenario is now fit and written under its
    own handler.
  * A non-converged fit wrote no run_time_rceattle_<scenario>.RDS, so its file
    set differed from a converged one. run_fims.R writes its run time
    unconditionally; this now does too. That state is expected on a
    deterministic case, where the random-effects scenarios time out.

Also corrects the fit_control() comment: newtonsteps = 1 is more than the
default 0, and phase = FALSE is the default, so "tighter and cheaper" was
wrong on two of three counts. The loopnum = 1 trade-off is now stated.

Smoke test green: 21 files, all three scenarios pdHess with max|grad| < 1e-06,
SSB correlation 0.9995.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant