Profiling & Benchmarking

Introduction

This vignette establishes performance baselines for the opal BET model prior to efficiency refactoring. It profiles the full model evaluation, isolates key sub-functions, and records AD tape metrics. After each refactoring change, re-run this vignette and compare against the saved baseline to verify numerical equivalence and measure speedup.

The profiling workflow uses:

• profvis — interactive flame graph profiler for identifying hotspots
• bench — precise microbenchmarking for comparing before/after
• system.time — simple wall-clock timing for quick checks

Load inputs

Load the opal package, RTMB, and profiling dependencies:

library(dplyr)
library(RTMB)
library(profvis)
library(bench)
library(opal)

The baseline data, parameters, and map are loaded from the bundled datasets. These were generated by the baseline.Rmd vignette to ensure consistency across all profiling and comparison tasks.

data(opal_baseline_data)
data <- opal_baseline_data

Length and Weight Data

The composition data and dimensions are already prepared in opal_baseline_data.

Model setup

Parameters

Load initial parameter values from the baseline parameters:

data(opal_baseline_parameters)
parameters <- opal_baseline_parameters

Priors

The priors are already included in the baseline data.

Parameter map

Load the parameter map from the baseline:

data(opal_baseline_map)
map <- opal_baseline_map

Build the AD object

t_build <- system.time({
  obj <- MakeADFun(func = cmb(opal_model, data),
                   parameters = parameters, map = map)
})
cat("MakeADFun build time:", round(t_build["elapsed"], 2), "sec\n")
#> MakeADFun build time: 22.05 sec
cat("Estimated parameters:", length(obj$par), "\n")
#> Estimated parameters: 270
cat("Parameter names:", paste(unique(names(obj$par)), collapse = ", "), "\n")
#> Parameter names: log_B0, log_cpue_q, rdev_y

Model dimensions

These dimensions drive the computational cost. Every unnecessary scalar loop or function call inside the dynamics is amplified by these values:

cat("n_year:    ", data$n_year, "\n")
#> n_year:     268
cat("n_age:     ", data$n_age, "\n")
#> n_age:      40
cat("n_fishery: ", data$n_fishery, "\n")
#> n_fishery:  15
cat("n_season:  ", data$n_season, "\n")
#> n_season:   1
cat("n_len:     ", data$n_len, "\n")
#> n_len:      95
cat("n_lf obs:  ", data$n_lf, "\n")
#> n_lf obs:   210
cat("lf_switch: ", data$lf_switch, "\n")
#> lf_switch:  1
cat("n_wt:      ", data$n_wt, "\n")
#> n_wt:       200
cat("n_wf obs:  ", data$n_wf, "\n")
#> n_wf obs:   112
cat("wf_switch: ", data$wf_switch, "\n")
#> wf_switch:  1

Baseline timing

Warm up the AD tape (the first call constructs it), then time fn() and gr() separately:

invisible(obj$fn(obj$par))
invisible(obj$gr(obj$par))
#> outer mgc:  13629.29

t_fn <- system.time(nll <- obj$fn(obj$par))
t_gr <- system.time(gr  <- obj$gr(obj$par))
#> outer mgc:  13629.29
cat("obj$fn():", round(t_fn["elapsed"], 4), "sec  (NLL =", round(nll, 4), ")\n")
#> obj$fn(): 0.001 sec  (NLL = 852582 )
cat("obj$gr():", round(t_gr["elapsed"], 4), "sec  (max|gr| =", round(max(abs(gr)), 6), ")\n")
#> obj$gr(): 0.027 sec  (max|gr| = 13629.29 )
cat("gr/fn ratio:", round(t_gr["elapsed"] / max(t_fn["elapsed"], 1e-6), 1), "x\n")
#> gr/fn ratio: 27 x

For a well-optimised AD tape, the gradient should cost roughly 2–4× the function evaluation. A much higher ratio suggests redundant operations on the tape that compound in reverse mode.

Repeated timing with bench::mark gives more stable estimates:

bm_core <- bench::mark(
  fn = obj$fn(obj$par),
  gr = obj$gr(obj$par),
  iterations = 20,
  check = FALSE
)
#> outer mgc:  13629.29 
#> outer mgc:  13629.29 
#> outer mgc:  13629.29 
#> outer mgc:  13629.29 
#> outer mgc:  13629.29 
#> outer mgc:  13629.29 
#> outer mgc:  13629.29 
#> outer mgc:  13629.29 
#> outer mgc:  13629.29 
#> outer mgc:  13629.29 
#> outer mgc:  13629.29 
#> outer mgc:  13629.29 
#> outer mgc:  13629.29 
#> outer mgc:  13629.29 
#> outer mgc:  13629.29 
#> outer mgc:  13629.29 
#> outer mgc:  13629.29 
#> outer mgc:  13629.29 
#> outer mgc:  13629.29 
#> outer mgc:  13629.29 
#> outer mgc:  13629.29
bm_core[, c("expression", "min", "median", "itr/sec","total_time")]
#> # A tibble: 2 × 4
#>   expression      min   median `itr/sec`
#>   <bch:expr> <bch:tm> <bch:tm>     <dbl>
#> 1 fn           13.3µs   14.1µs   56284. 
#> 2 gr           24.4ms     25ms      39.9

AD tape diagnostics

The AD tape is the sequence of elementary operations that RTMB records during the first fn() call. Tape size directly determines gradient cost: fewer operations = faster gr().

tryCatch({
  tape_env <- obj$env
  if (!is.null(tape_env$ADFun)) {
    cat("ADFun object found in obj$env\n\n")
    capture.output(print(tape_env$ADFun)) |> head(20) |> cat(sep = "\n")
  }
}, error = function(e) {
  cat("Could not access tape info:", e$message, "\n")
})
#> ADFun object found in obj$env
#> 
#> $ptr
#> <pointer: 0x5562e72161f0>
#> attr(,"par")
#>       log_B0   log_cpue_q       rdev_y       rdev_y       rdev_y       rdev_y 
#> 20.000000000  0.000000000  0.057843051  0.250649244 -0.297200433 -0.332988095 
#>       rdev_y       rdev_y       rdev_y       rdev_y       rdev_y       rdev_y 
#>  0.436765054  0.060352409  0.081547956  0.336249134  0.291607395 -0.828532782 
#>       rdev_y       rdev_y       rdev_y       rdev_y       rdev_y       rdev_y 
#> -0.802114156  0.833572916 -0.139021343  0.533254609  0.453715679 -0.069286219 
#>       rdev_y       rdev_y       rdev_y       rdev_y       rdev_y       rdev_y 
#> -1.140285410  0.476695933  0.794901730 -0.115919747  0.415578293  0.221393600 
#>       rdev_y       rdev_y       rdev_y       rdev_y       rdev_y       rdev_y 
#> -0.105304159 -0.616215018 -0.471341695 -0.064217582  0.201819599  0.501254175 
#>       rdev_y       rdev_y       rdev_y       rdev_y       rdev_y       rdev_y 
#> -0.162970476  0.012973481  0.114400539  0.322098968 -0.988475967  0.453626213 
#>       rdev_y       rdev_y       rdev_y       rdev_y       rdev_y       rdev_y 
#> -0.041498056 -0.415121560  0.049271116 -0.202691687  0.794726577 -0.377310752 
#>       rdev_y       rdev_y       rdev_y       rdev_y       rdev_y       rdev_y 
#> -0.460181950 -0.023642665  0.487765666 -0.562912712 -0.962409020 -0.297170182 
#>       rdev_y       rdev_y       rdev_y       rdev_y       rdev_y       rdev_y

cat("\nobj total size:  ", round(object.size(obj) / 1e6, 1), "MB\n")
#> 
#> obj total size:   0.8 MB
cat("obj$env size:    ", round(object.size(obj$env) / 1e6, 1), "MB\n")
#> obj$env size:     0 MB

profvis — full model evaluation

profvis samples the R call stack at regular intervals and produces an interactive flame graph. Run 10 iterations of fn() + gr() to accumulate enough samples:

prof <- profvis({
  for (i in 1:10) {
    obj$fn(obj$par)
    obj$gr(obj$par)
  }
}, simplify = FALSE)

htmlwidgets::saveWidget(prof, "profile_full.html", selfcontained = TRUE)

Open profile_full.html in a browser to inspect the flame graph. Key things to look for:

• Time in get_harvest_rate vs do_dynamics vs opal_model — if get_harvest_rate is a significant fraction, Change 1 (inlining) will help.
• Time in get_length_like / get_weight_like / get_cpue_like — if the composition likelihoods dominate, Changes 7–9 are important. The weight likelihood adds an extra rebinning matrix multiply per observation.
• Time in getAll() calls — indicates overhead from unpacking data/parameter lists.
• Time in .Call (C++ / RTMB internals) vs R code — if most time is in .Call, the bottleneck is tape complexity rather than R dispatch overhead.

Component-level profiling

Benchmark individual functions in plain R (outside the AD tape) to understand their base cost. These timings won’t include AD overhead but show the algorithmic cost of each component.

opal_model

The highest-level entry point bundles all individual components together: growth, selectivity, dynamics, and likelihoods. Timing opal_model() directly outside the AD tape gives a sense of the un-taped work that the tape must eventually record.

bm_opal_model <- bench::mark(
  opal_model = opal_model(parameters = parameters, data = data),
  iterations = 20,
  check = FALSE
)
bm_opal_model[, c("expression", "min", "median", "itr/sec","total_time")]
#> # A tibble: 1 × 4
#>   expression      min   median `itr/sec`
#>   <bch:expr> <bch:tm> <bch:tm>     <dbl>
#> 1 opal_model    153ms    154ms      6.50

prof_opal <- profvis({
  for (i in 1:10) {
    opal_model(parameters = parameters, data = data)
  }
}, simplify = FALSE)

htmlwidgets::saveWidget(prof_opal, "profile_opal.html", selfcontained = TRUE)

get_pla

The probability-of-length-at-age matrix is a n_len × n_age matrix computed via pnorm. Currently uses a scalar inner loop over length bins (see Change 4):

par_list <- obj$env$parList(obj$env$last.par.best)
L1  <- exp(par_list$log_L1)
L2  <- exp(par_list$log_L2)
k   <- exp(par_list$log_k)
CV1 <- exp(par_list$log_CV1)
CV2 <- exp(par_list$log_CV2)

len_lower <- seq(from = data$len_bin_start, by = data$len_bin_width,
                 length.out = data$n_len)
len_upper <- len_lower + data$len_bin_width
len_mid   <- len_lower + data$len_bin_width / 2

mu_a <- get_growth(data$n_age, data$A1, data$A2, L1, L2, k)
sd_a <- get_sd_at_age(mu_a, L1, L2, CV1, CV2)

bm_pla <- bench::mark(
  get_pla = get_pla(len_lower, len_upper, mu_a, sd_a),
  iterations = 50,
  check = FALSE
)
bm_pla[, c("expression", "min", "median", "itr/sec","total_time")]
#> # A tibble: 1 × 4
#>   expression      min   median `itr/sec`
#>   <bch:expr> <bch:tm> <bch:tm>     <dbl>
#> 1 get_pla       543µs    580µs     1683.

# prof_pla <- profvis({
#   for (i in 1:10) {
#     get_pla(len_lower, len_upper, mu_a, sd_a)
#   }
# }, simplify = FALSE)

get_selectivity

Selectivity is computed per fishery via logistic or double-normal curves applied through the PLA matrix:

pla <- get_pla(len_lower, len_upper, mu_a, sd_a)

bm_sel <- bench::mark(
  get_selectivity = get_selectivity(data, par_list$par_sel, pla, len_mid),
  iterations = 50,
  check = FALSE
)
bm_sel[, c("expression", "min", "median", "itr/sec","total_time")]
#> # A tibble: 1 × 4
#>   expression           min   median `itr/sec`
#>   <bch:expr>      <bch:tm> <bch:tm>     <dbl>
#> 1 get_selectivity      5ms   5.26ms      188.

PLA matrix multiply

The core operation inside get_length_like — a n_len × n_age matrix-vector multiply that converts catch-at-age to catch-at-length. This executes once per LF observation:

catch_a <- rep(1, data$n_age)

bm_pla_mult <- bench::mark(
  `pla %*% catch_a` = c(pla %*% catch_a),
  iterations = 200,
  check = FALSE
)
bm_pla_mult[, c("expression", "min", "median", "itr/sec","total_time")]
#> # A tibble: 1 × 4
#>   expression           min   median `itr/sec`
#>   <bch:expr>      <bch:tm> <bch:tm>     <dbl>
#> 1 pla %*% catch_a    2.5µs   3.03µs   322473.

get_length_like

The full length composition likelihood, called with plain-R arguments (not AD-taped). This includes the per-observation matrix multiply, bin aggregation, normalisation, and dmultinom evaluation:

invisible(obj$fn(obj$par))
rep <- obj$report()

lf_year_fi <- split(data$lf_year, data$lf_fishery)
lf_n_fi    <- split(data$lf_n, data$lf_fishery)

bm_lf <- bench::mark(
  get_length_like = get_length_like(
    lf_obs_flat      = data$lf_obs_flat,
    lf_obs_ints      = data$lf_obs_ints,
    lf_obs_prop      = data$lf_obs_prop,
    catch_pred_fya   = rep$catch_pred_fya,
    pla              = pla,
    lf_n_f           = data$lf_n_f,
    lf_fishery_f     = data$lf_fishery_f,
    lf_year_fi       = lf_year_fi,
    lf_n_fi          = lf_n_fi,
    lf_minbin        = data$lf_minbin,
    lf_maxbin        = data$lf_maxbin,
    removal_switch_f = data$removal_switch_f,
    lf_switch        = data$lf_switch,
    n_len            = data$n_len,
    n_lf             = data$n_lf,
    log_lf_tau       = par_list$log_lf_tau
  ),
  iterations = 50,
  check = FALSE
)
bm_lf[, c("expression", "min", "median", "itr/sec","total_time")]
#> # A tibble: 1 × 4
#>   expression           min   median `itr/sec`
#>   <bch:expr>      <bch:tm> <bch:tm>     <dbl>
#> 1 get_length_like   5.48ms   5.58ms      175.

split() overhead

These split() calls currently execute inside the AD-taped function on every model evaluation (see Change 8):

bm_split <- bench::mark(
  split_year = split(data$lf_year, data$lf_fishery),
  split_n    = split(data$lf_n, data$lf_fishery),
  iterations = 500,
  check = FALSE
)
bm_split[, c("expression", "min", "median", "itr/sec","total_time")]
#> # A tibble: 2 × 4
#>   expression      min   median `itr/sec`
#>   <bch:expr> <bch:tm> <bch:tm>     <dbl>
#> 1 split_year   15.7µs   17.4µs    55069.
#> 2 split_n      15.7µs   16.9µs    56802.

Rebinning matrix multiply

The weight likelihood adds an extra matrix multiply per WF observation: wf_rebin_matrix %*% pred_at_length converts the n_len-vector of predicted length composition to an n_wt-vector of predicted weight composition. This is a n_wt × n_len matrix-vector product:

pred_at_length <- rep(1, data$n_len)

bm_rebin_mult <- bench::mark(
  `rebin %*% pred_l` = c(data$wf_rebin_matrix %*% pred_at_length),
  iterations = 200,
  check = FALSE
)
bm_rebin_mult[, c("expression", "min", "median", "itr/sec","total_time")]
#> # A tibble: 1 × 4
#>   expression            min   median `itr/sec`
#>   <bch:expr>       <bch:tm> <bch:tm>     <dbl>
#> 1 rebin %*% pred_l   14.3µs   16.9µs    37043.

get_weight_like

The full weight composition likelihood, called with plain-R arguments (not AD-taped). This includes the per-observation PLA multiply, rebinning to weight bins, normalisation, and likelihood evaluation:

invisible(obj$fn(obj$par))
rep <- obj$report()

wf_year_fi <- split(data$wf_year, data$wf_fishery)
wf_n_fi    <- split(data$wf_n, data$wf_fishery)

bm_wf <- bench::mark(
  get_weight_like = get_weight_like(
    wf_obs_flat      = data$wf_obs_flat,
    wf_obs_ints      = data$wf_obs_ints,
    wf_obs_prop      = data$wf_obs_prop,
    catch_pred_fya   = rep$catch_pred_fya,
    pla              = pla,
    wf_rebin_matrix  = data$wf_rebin_matrix,
    wf_n_f           = data$wf_n_f,
    wf_fishery_f     = data$wf_fishery_f,
    wf_year_fi       = wf_year_fi,
    wf_n_fi          = wf_n_fi,
    wf_minbin        = data$wf_minbin,
    wf_maxbin        = data$wf_maxbin,
    removal_switch_f = data$removal_switch_f,
    wf_switch        = data$wf_switch,
    n_wt             = data$n_wt,
    n_wf             = data$n_wf,
    log_wf_tau       = par_list$log_wf_tau
  ),
  iterations = 50,
  check = FALSE
)
bm_wf[, c("expression", "min", "median", "itr/sec","total_time")]
#> # A tibble: 1 × 4
#>   expression           min   median `itr/sec`
#>   <bch:expr>      <bch:tm> <bch:tm>     <dbl>
#> 1 get_weight_like   8.15ms   8.34ms      119.

WF split() overhead

These split() calls for WF data also execute inside the AD-taped function on every model evaluation, mirroring the LF split overhead:

bm_wf_split <- bench::mark(
  split_wf_year = split(data$wf_year, data$wf_fishery),
  split_wf_n    = split(data$wf_n, data$wf_fishery),
  iterations = 500,
  check = FALSE
)
bm_wf_split[, c("expression", "min", "median", "itr/sec","total_time")]
#> # A tibble: 2 × 4
#>   expression         min   median `itr/sec`
#>   <bch:expr>    <bch:tm> <bch:tm>     <dbl>
#> 1 split_wf_year   14.4µs   15.2µs    64128.
#> 2 split_wf_n      14.1µs   15.1µs    64208.

Likelihood component breakdown

Inspect the contribution of each likelihood component to the total NLL:

cat("lp_prior:  ", round(rep$lp_prior, 4), "\n")
#> lp_prior:   0
cat("lp_penalty:", round(rep$lp_penalty, 4), "\n")
#> lp_penalty: 0
cat("lp_rec:    ", round(rep$lp_rec, 4), "\n")
#> lp_rec:     205.8579
cat("lp_cpue:   ", round(sum(rep$lp_cpue), 4),
    " (sum of", length(rep$lp_cpue), "obs)\n")
#> lp_cpue:    490.8378  (sum of 268 obs)
cat("lp_lf:     ", round(sum(rep$lp_lf), 4),
    " (sum of", length(rep$lp_lf), "obs)\n")
#> lp_lf:      133665.5  (sum of 210 obs)
cat("lp_wf:     ", round(sum(rep$lp_wf), 4),
    " (sum of", length(rep$lp_wf), "obs)\n")
#> lp_wf:      718157.9  (sum of 112 obs)
cat("Total NLL: ", round(nll, 4), "\n")
#> Total NLL:  852582

Estimated tape operation counts

Before profiling, we can estimate the theoretical number of scalar AD tape operations from the model dimensions. These estimates help prioritise which refactoring changes will have the largest impact:

n_yr  <- data$n_year
n_a   <- data$n_age
n_f   <- data$n_fishery
n_s   <- data$n_season
n_l   <- data$n_len
n_obs <- data$n_lf
n_w   <- data$n_wt
n_wobs <- data$n_wf

cat("Dynamics loop:\n")
#> Dynamics loop:
cat("  Harvest rate function calls: ", n_yr, "×", n_s, "=", n_yr * n_s, "\n")
#>   Harvest rate function calls:  268 × 1 = 268
cat("  Scalar catch pred ops:       ", n_f, "×", n_a, "×", n_yr * n_s,
    "=", n_f * n_a * n_yr * n_s, "\n")
#>   Scalar catch pred ops:        15 × 40 × 268 = 160800

cat("\nTime-invariant array overhead:\n")
#> 
#> Time-invariant array overhead:
cat("  sel_fya year replication:  ", n_f, "×", n_yr, "=", n_f * n_yr, "\n")
#>   sel_fya year replication:   15 × 268 = 4020
cat("  weight_fya replication:    ", n_f, "×", n_yr, "=", n_f * n_yr, "\n")
#>   weight_fya replication:     15 × 268 = 4020

cat("\nPLA construction:\n")
#> 
#> PLA construction:
cat("  Current (scalar pnorm):    ", n_a, "×", n_l, "=", n_a * n_l, "\n")
#>   Current (scalar pnorm):     40 × 95 = 3800
cat("  After vectorisation:       ", n_a, "× 2 =", n_a * 2, "\n")
#>   After vectorisation:        40 × 2 = 80

cat("\nLength likelihood:\n")
#> 
#> Length likelihood:
cat("  Matrix multiplies:         ", n_obs, "×", n_l, "×", n_a, "=",
    n_obs * n_l * n_a, "multiply-adds\n")
#>   Matrix multiplies:          210 × 95 × 40 = 798000 multiply-adds
cat("  if-branch evaluations:     ", n_obs, "× 3 =", n_obs * 3, "\n")
#>   if-branch evaluations:      210 × 3 = 630

cat("\nWeight likelihood:\n")
#> 
#> Weight likelihood:
cat("  PLA multiplies:            ", n_wobs, "×", n_l, "×", n_a, "=",
    n_wobs * n_l * n_a, "multiply-adds\n")
#>   PLA multiplies:             112 × 95 × 40 = 425600 multiply-adds
cat("  Rebin multiplies:          ", n_wobs, "×", n_w, "×", n_l, "=",
    n_wobs * n_w * n_l, "multiply-adds\n")
#>   Rebin multiplies:           112 × 200 × 95 = 2128000 multiply-adds
cat("  if-branch evaluations:     ", n_wobs, "× 3 =", n_wobs * 3, "\n")
#>   if-branch evaluations:      112 × 3 = 336

scalar_total <- n_f * n_a * n_yr * n_s + n_a * n_l + n_f * n_yr * 2 +
                n_wobs * n_w * n_l
cat("\nTotal estimated scalar tape ops: ~", scalar_total, "\n")
#> 
#> Total estimated scalar tape ops: ~ 2300640

Optimisation timing

Uncomment the block below to time a full nlminb optimisation run. This is useful for measuring end-to-end improvement but takes several minutes:

Lwr <- rep(-Inf, length(obj$par))
Upr <- rep(Inf, length(obj$par))
Lwr[grep("log_B0", names(obj$par))] <- log(1)
Upr[grep("log_B0", names(obj$par))] <- log(exp(22))
Lwr[grep("log_cpue_q", names(obj$par))] <- log(0.1)
Upr[grep("log_cpue_q", names(obj$par))] <- log(10)
Lwr[grep("rdev_y", names(obj$par))] <- -5
Upr[grep("rdev_y", names(obj$par))] <- 5

control <- list(eval.max = 10000, iter.max = 10000)

t_opt <- system.time({
  opt <- nlminb(start = obj$par, objective = obj$fn, gradient = obj$gr,
                hessian = obj$he, lower = Lwr, upper = Upr, control = control)
})
cat("Optimization time: ", round(t_opt["elapsed"], 1), "sec\n")
cat("Convergence:       ", ifelse(opt$convergence == 0, "TRUE", "FALSE"), "\n")
cat("Final NLL:         ", round(opt$objective, 4), "\n")
cat("Function evals:    ", opt$evaluations["function"], "\n")
cat("Gradient evals:    ", opt$evaluations["gradient"], "\n")
cat("Sec per fn eval:   ", round(t_opt["elapsed"] / opt$evaluations["function"], 4), "\n")
cat("Sec per gr eval:   ", round(t_opt["elapsed"] / opt$evaluations["gradient"], 4), "\n")

Save baseline

Save the current results to an RDS file for later comparison. Note that for long-term consistency, the reference baseline targets are also available in the opal_baseline dataset.

baseline <- list(
  timestamp   = Sys.time(),
  description = "Pre-refactoring baseline",
  dimensions  = list(
    n_year    = data$n_year,
    n_age     = data$n_age,
    n_fishery = data$n_fishery,
    n_season  = data$n_season,
    n_len     = data$n_len,
    n_lf      = data$n_lf,
    n_wt      = data$n_wt,
    n_wf      = data$n_wf
  ),
  nll         = nll,
  max_gr      = max(abs(gr)),
  timing      = list(
    fn_median = as.numeric(bm_core$median[1]),
    gr_median = as.numeric(bm_core$median[2])
  ),
  report      = list(
    spawning_biomass_y = rep$spawning_biomass_y,
    catch_pred_fya     = rep$catch_pred_fya,
    lp_cpue            = rep$lp_cpue,
    lp_lf              = rep$lp_lf,
    lp_wf              = rep$lp_wf,
    number_ysa         = rep$number_ysa
  ),
  par_values  = obj$par
)

saveRDS(baseline, "opal_baseline.rds")

To compare after refactoring using the bundled opal_baseline dataset:

data(opal_baseline)
old <- opal_baseline

# Numerical equivalence (should be ~1e-14)
max(abs(rep$number_ysa - old$report$number_ysa))
max(abs(rep$catch_pred_fya - old$report$catch_pred_fya))
max(abs(rep$spawning_biomass_y - old$report$spawning_biomass_y))
abs(nll - old$nll)

# Speedup
cat("fn speedup:", round(old$timing$fn_median / as.numeric(bm_core$median[1]), 1), "x\n")
cat("gr speedup:", round(old$timing$gr_median / as.numeric(bm_core$median[2]), 1), "x\n")