Skip to content

refactoring Stan code for efficiency [WIP]#1274

Open
bob-carpenter wants to merge 3 commits intoepiforecasts:mainfrom
bob-carpenter:1273-stan-performance
Open

refactoring Stan code for efficiency [WIP]#1274
bob-carpenter wants to merge 3 commits intoepiforecasts:mainfrom
bob-carpenter:1273-stan-performance

Conversation

@bob-carpenter
Copy link

@bob-carpenter bob-carpenter commented Jan 29, 2026

Work in Progress: do not merge

I am creating this intermediate PR to see if the EpiNow2 devs want this kind of PR and to ask what kind of testing you want before proceeding to refactor the rest of the code base for efficiency.

Description

This PR closes #1273

Initial submission checklist

  • My PR is based on a package issue and I have explicitly linked it.
  • I have tested my changes locally (using devtools::test() and devtools::check()).
  • I have added or updated unit tests where necessary.
  • I have updated the documentation if required and rebuilt docs if yes (using devtools::document()).
  • I have followed the established coding standards (and checked using lintr::lint_package()).
  • I have added a news item linked to this PR.

I'm not an R developer. I installed devtools::test(), but it's failing with an error I don't understand that looks like it's Stan related as it mentions a Stan signature. I made sure that all of my Stan changes compile in Stan, so I'm not sure where it's coming from or how to track this down. I'll append the error dump below.

I don't know what you're expecting from the linter run. The linter dumps out dozens of pages of lint errors when I run it in R on files I haven't touched. Is there a way to see if the current PR introduced new errors?

After the initial Pull Request

  • I have reviewed Checks for this PR and addressed any issues as far as I am able.

ERROR I'M GETTING WITH devtools::test()

> devtools::test()
ℹ Testing EpiNow2
Error in `(function (command = NULL, args = character(), error_on_status = TRUE, …`:
! System command 'R' failed
---
Exit status: 1
Stdout & stderr (last 10 lines, see `$stderr` for more):
(int, real, real)
Available signatures:
(data int, data real, data real) => vector
  The second argument must be data-only. (Local variables are assumed to
  depend on parameters; same goes for function inputs unless they are marked
  with the keyword 'data'.)
Calls: <Anonymous> -> sapply -> lapply -> FUN -> <Anonymous>
Execution halted
ERROR: configuration failed for package ‘EpiNow2’
* removing ‘/private/var/folders/8f/5nrkg18127x2b8rq18tshp5r0000gr/T/RtmpazDR3K/devtools_install_128b95e5bd56/EpiNow2’
---
Type .Last.error to see the more details
> .Last.error
<system_command_status_error/rlib_error_3_0/rlib_error/error>
Error in `(function (command = NULL, args = character(), error_on_status = TRUE, …`:
! System command 'R' failed
---
Exit status: 1
Stdout & stderr (last 10 lines, see `$stderr` for more):
(int, real, real)
Available signatures:
(data int, data real, data real) => vector
  The second argument must be data-only. (Local variables are assumed to
  depend on parameters; same goes for function inputs unless they are marked
  with the keyword 'data'.)
Calls: <Anonymous> -> sapply -> lapply -> FUN -> <Anonymous>
Execution halted
ERROR: configuration failed for package ‘EpiNow2’
* removing ‘/private/var/folders/8f/5nrkg18127x2b8rq18tshp5r0000gr/T/RtmpazDR3K/devtools_install_128b95e5bd56/EpiNow2’
---
Backtrace:
 1. devtools::test()
 2. testthat::test_local(pkg$path, filter = filter, stop_on_failure = stop_on_failure, …
 3. testthat::test_dir(test_path, package = package, reporter = reporter, ..., …
 4. testthat:::test_files(test_dir = path, test_paths = test_paths, test_package = package, …
 5. testthat:::test_files_serial(test_dir = test_dir, test_package = test_package, …
 6. testthat:::test_files_setup_env(test_package, test_dir, load_package, env)
 7. pkgload::load_all(test_dir, export_all = args[["export_all"]], …
 8. pkgbuild::compile_dll(path, quiet = quiet, debug = debug)
 9. pkgbuild:::withr_with_makevars(compiler_flags(debug), { …
10. pkgbuild:::withr_with_envvar(c(R_MAKEVARS_USER = makevars_file), { …
11. base::force(code)
12. base::force(code)
13. pkgbuild:::withr_with_envvar(c(DEBUG = "true"), build())
14. base::force(code)
15. pkgbuild::build()
16. pkgbuild:::install_min(path, dest = install_dir, components = "libs", args = if (needs_clean(path)) "--preclean", …
17. pkgbuild::rcmd_build_tools("INSTALL", c(path, paste("--library=", dest, …
18. pkgbuild::with_build_tools({ …
19. base::withCallingHandlers(callr::rcmd_safe(..., env = env, spinner = FALSE, …
20. callr::rcmd_safe(..., env = env, spinner = FALSE, show = FALSE, …
21. callr:::run_r(options)
22. base::with(options, with_envvar(env, do.call(processx::run, c(list(bin, …
23. base::with.default(options, with_envvar(env, do.call(processx::run, …
24. base::eval(substitute(expr), data, enclos = parent.frame())
25. base::eval(substitute(expr), data, enclos = parent.frame())
26. callr:::with_envvar(env, do.call(processx::run, c(list(bin, args = real_cmdargs, …
27. base::force(code)
28. base::do.call(processx::run, c(list(bin, args = real_cmdargs, stdout_line_callback = real_callback(stdout), …
29. (function (command = NULL, args = character(), error_on_status = TRUE, …
30. base::throw(new_process_error(res, call = sys.call(), echo = echo, …
31. | base::signalCondition(cond)
32. (function (e) …
33. asNamespace("callr")$err$throw(e)

Summary by CodeRabbit

  • New Features

    • Exported new Matern indices helper function for public use
  • Refactor

    • Optimised Gaussian Process spectral density computations for Matern kernels
    • Streamlined periodic spectral density calculations
    • Simplified noise setup process
  • Bug Fixes

    • Improved error messages for unsupported Matern parameter values

✏️ Tip: You can customize this high-level summary in your review settings.

@coderabbitai
Copy link
Contributor

coderabbitai bot commented Jan 29, 2026

Important

Review skipped

Auto incremental reviews are disabled on this repository.

Please check the settings in the CodeRabbit UI or the .coderabbit.yaml file in this repository. To trigger a single review, invoke the @coderabbitai review command.

You can disable this status message by setting the reviews.review_status to false in the CodeRabbit configuration file.

  • 🔍 Trigger a full review

Walkthrough

The pull request modifies two Stan function files. In convolve.stan, the loop bounds logic is adjusted and allocation timing changed. In gaussian_process.stan, a new matern_indices helper function is introduced and Matern spectral density computations are refactored to use it.

Changes

Cohort / File(s) Summary
Convolution function
inst/stan/functions/convolve.stan
Loop bound changed to use max(len, xlen) instead of xlen; vector z allocation moved into loop; post-loop convolution block for len > xlen case removed.
Gaussian process functions
inst/stan/functions/gaussian_process.stan
New exported helper function matern_indices(int M, int L) added; diagSPD_Matern12, diagSPD_Matern32, and diagSPD_Matern52 refactored to use matern_indices; PHI_periodic rewritten to compute w0xk via row vector; setup_noise simplified; error message updated for Matern nu values.
🚥 Pre-merge checks | ✅ 4 | ❌ 1
❌ Failed checks (1 warning)
Check name Status Explanation Resolution
Out of Scope Changes check ⚠️ Warning All changes are within scope: convolve.stan and gaussian_process.stan refactoring directly support the stated objective of Stan code performance improvements; however, the convolve.stan changes introduce potential out-of-bounds issues. Review and fix the convolve.stan loop logic to prevent out-of-bounds writes to z when len > xlen; ensure changes maintain algorithmic correctness without introducing safety issues.
✅ Passed checks (4 passed)
Check name Status Explanation
Description Check ✅ Passed Check skipped - CodeRabbit’s high-level summary is enabled.
Title check ✅ Passed The title 'refactoring Stan code for efficiency [WIP]' accurately describes the main changes, which involve refactoring convolve.stan and gaussian_process.stan functions for performance improvements.
Linked Issues check ✅ Passed The PR addresses the objective from issue #1273 to identify and apply micro-optimizations in Stan code; changes to convolve.stan and gaussian_process.stan represent optimisation efforts.
Docstring Coverage ✅ Passed No functions found in the changed files to evaluate docstring coverage. Skipping docstring coverage check.

✏️ Tip: You can configure your own custom pre-merge checks in the settings.

✨ Finishing touches
🧪 Generate unit tests (beta)
  • Create PR with unit tests
  • Post copyable unit tests in a comment

Tip

🧪 Unit Test Generation v2 is now available!

We have significantly improved our unit test generation capabilities.

To enable: Add this to your .coderabbit.yaml configuration:

reviews:
  finishing_touches:
    unit_tests:
      enabled: true

Try it out by using the @coderabbitai generate unit tests command on your code files or under ✨ Finishing Touches on the walkthrough!

Have feedback? Share your thoughts on our Discord thread!


Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out.

❤️ Share

Comment @coderabbitai help to get the list of available commands and usage tips.

Copy link
Contributor

@coderabbitai coderabbitai bot left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actionable comments posted: 3

🤖 Fix all issues with AI agents
In `@inst/stan/functions/convolve.stan`:
- Around line 74-79: The loop currently calls calc_conv_indices_xlen for all s
up to max(len, xlen), which causes end_x = s to exceed xlen when s > xlen and
leads to out-of-bounds access; modify the loop in the convolve implementation so
that for s <= xlen it uses calc_conv_indices_xlen but for s > xlen it uses
calc_conv_indices_len (or otherwise ensure end_x is capped at xlen), updating
z[s] = dot_product(...) to use the indices returned by the correct helper (refer
to calc_conv_indices_xlen, calc_conv_indices_len and the loop variable s) so no
slice of x exceeds xlen.

In `@inst/stan/functions/gaussian_process.stan`:
- Around line 69-71: The code declares real factor1 but returns factor which is
undefined; change the return to use factor1 (i.e., return factor1 ./ denom) so
the function references the declared variable factor1 instead of the
non-existent factor in gaussian_process.stan (look for the declarations of
factor1 and denom and the return statement).
- Around line 36-38: The function matern_indices currently declares parameter L
as int but callers (diagSPD_Matern12, diagSPD_Matern32, diagSPD_Matern52) pass a
real, causing a Stan type mismatch; change the signature of matern_indices to
accept real L instead of int, and ensure uses inside (e.g., factor = pi() / (2 *
L) and linspaced_vector(M, factor, factor * M)) remain consistent with real
arithmetic so the function compiles and matches its callers.
🧹 Nitpick comments (1)
inst/stan/functions/gaussian_process.stan (1)

86-89: Unused variable indices.

The indices vector on line 86 is computed but never used — the function now relies on matern_indices(M, L) instead. This is dead code that should be removed.

♻️ Proposed fix to remove dead code
 vector diagSPD_Matern52(real alpha, real rho, real L, int M) {
-  vector[M] indices = linspaced_vector(M, 1, M);
   real factor = 16 * pow(sqrt(5) / rho, 5);
   vector[M] denom = 3 * pow(5 / square(rho) + matern_indices(M, L), 3);
   return alpha * sqrt(factor ./ denom);
 }

@seabbs
Copy link
Contributor

seabbs commented Jan 30, 2026

Hi @bob-carpenter,

Thank you for this, and yes we are very keen to have some help with efficiency tuning. Sorry its a bit annoying to work with.

Stan error: This seems like an odd one. I would imagine that the problem here is due to a mismatch between stan versions in rstan and cmdstanr. More recent work from us is using the latter but this is still on rstan so we can provide an easy to use package to Public health users (who generally can't install cmdstan). If it isn't that than I am not clear what it could be.

Linting: That isn't ideal. We have the lint-changed-files CI action that shows lints for just PR diff (unfortunately, that is failing due to the Stan issue).

In terms of testing, our CI has correctness checks and benchmarking (it will post as a comment if not failing!) so it should be enough to bounce off, I would hope. When I do local development I use our "cmdstanr" backend for ease (https://epiforecasts.io/EpiNow2/reference/stan_sampling_opts.html?q=cmdstan#arg-backend) and then use devtools::test() as you are doing.

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.

improve performance of Stan code

2 participants