Refactor Stan GP and convolution functions for efficiency

[seabbs-bot]

Description

It picks up the efficiency refactoring started in #1274 by @bob-carpenter, which stalled because the Stan code would not compile under rstan.

The blocker was the new matern_indices() helper calling linspaced_vector(M, factor, factor * M) with factor = pi() / (2 * L). The Stan version shipped with rstan requires the bounds of linspaced_vector() to be data-only, and factor derives from the function argument L, which Stan treats as parameter-dependent. The build succeeds under cmdstanr (the dev backend) because that constraint is relaxed there, so the failure only showed up for contributors using rstan.

This PR applies the same lessons while keeping the code rstan-compatible:

• Added a shared matern_indices() helper for the Matern spectral densities, written so the linspaced_vector() bounds stay data-only (literals) and the scaling is applied afterwards. diagSPD_Matern12/32/52() now use it.
• Built the GP basis functions in PHI() and PHI_periodic() via an outer product (vector * row_vector) instead of rep_matrix() plus diag_post_multiply(), avoiding materialising the replicate matrix.
• Tidied convolve_with_rev_pmf() by allocating z after the validity checks and dropping the redundant if (len > xlen) guard (the loop runs zero times otherwise).
• Clarified the nu reject message to use the decimal values it actually compares against.

The Matern refactors are algebraically identical to the originals, and the existing Stan unit tests pass unchanged. A new test covers matern_indices().

Per maintainer preference, the no-op changes from #1274 (inlining setup_noise()'s return, reformatting diagSPD_Periodic()) were left out to keep the diff focused on genuine efficiency and DRY wins.

Initial submission checklist

• My PR is based on a package issue and I have explicitly linked it.
• I have tested my changes locally (expose_stan_fns() compiles under rstan; test-stan-guassian-process.R and test-stan-convole.R pass).
• I have added or updated unit tests where necessary.
• I have updated the documentation if required and rebuilt docs if yes (Stan-only changes, docstrings updated).
• I have followed the established coding standards.
• I have added a news item linked to this PR.

After the initial Pull Request

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

This was opened by a bot. Please ping @seabbs for any questions.

Summary by CodeRabbit

Release Notes

• Refactor

  • Improved efficiency of Gaussian process and convolution functions through code restructuring and shared utility functions.
  • Simplified mathematical operations in spectral-density computations.
  • Updated parameter validation handling.

• Tests

  • Added comprehensive tests to validate refactored functionality.

• Chores

  • Updated documentation to reflect internal optimisations.

[coderabbitai]

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.

⚙️ Run configuration

Configuration used: Organization UI

Review profile: CHILL

Plan: Pro

Run ID: 79b8a446-d052-494f-b931-2dee9ae96ad0

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

Use the checkbox below for a quick retry:

• 🔍 Trigger review

Walkthrough

This PR refactors Gaussian process and convolution Stan functions for improved performance. It introduces a shared matern_indices() helper to eliminate redundant spectral-index computations, simplifies basis function calculations by replacing rep_matrix/diag_post_multiply patterns with elementwise operations, optimises convolution loop structure, and updates release notes.

Changes

Gaussian Process and Convolution Optimization

Layer / File(s)	Summary
Matern indices helper and spectral density refactors inst/stan/functions/gaussian_process.stan, tests/testthat/test-stan-guassian-process.R	New matern_indices(int M, real L) helper computes the shared squared spectral-indices denominator used across diagSPD_Matern12, diagSPD_Matern32, and diagSPD_Matern52. All three functions rewritten to call the helper instead of manually constructing indices and power expressions. Test validates output length, positivity, and numerical correctness.
Basis function simplifications inst/stan/functions/gaussian_process.stan	PHI and PHI_periodic functions refactored to replace rep_matrix/diag_post_multiply patterns with direct elementwise row-vector operations. update_gp validation error for invalid nu adjusted to use decimal notation (0.5, 1.5, 2.5).
Convolution loop and declaration refactor inst/stan/functions/convolve.stan	Output vector declaration moved after validation checks in convolve_with_rev_pmf. Loop guard condition removed from (xlen + 1):len iteration, relying on Stan to run the loop zero times when the range is empty.
Release notes documentation NEWS.md	Release notes updated to document the efficiency improvements and new matern_indices() helper across Gaussian process and convolution functions.

🚥 Pre-merge checks | ✅ 5

✅ Passed checks (5 passed)

Check name	Status	Explanation
Title check	✅ Passed	The title accurately summarizes the main changes: refactoring Stan GP and convolution functions for efficiency improvements.
Linked Issues check	✅ Passed	The PR successfully addresses issue #1273 by implementing multiple micro-optimizations in Stan code including a shared matern_indices() helper, refactored basis functions, and improved convolution logic.
Out of Scope Changes check	✅ Passed	All changes are directly related to the stated objective of performance optimisation for Stan code. The NEWS.md update, tests, and refactored functions are all within scope.
Docstring Coverage	✅ Passed	No functions found in the changed files to evaluate docstring coverage. Skipping docstring coverage check.
Description Check	✅ Passed	Check skipped - CodeRabbit’s high-level summary is enabled.

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

✨ Finishing Touches

🧪 Generate unit tests (beta)

• Create PR with unit tests
• Commit unit tests in branch stan-efficiency-refactor

---

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

• X
• Mastodon
• Reddit
• LinkedIn

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

[seabbs-bot]

Automated review (quality gate)

Reviewed the diff against main. No Critical or Important issues found. Summary:

Correctness — verified algebraically identical to the originals:

• matern_indices(M, L) = (pi/(2L) * i)^2, matching the squared term in every Matern denominator.
• diagSPD_Matern12: new 1/rho + rho*matern_indices equals original rho*((1/rho)^2 + (pi/2L*i)^2). ✓
• diagSPD_Matern32: new 3/rho^2 + matern_indices equals original (sqrt(3)/rho)^2 + (pi/2L*i)^2; factor unchanged. ✓
• diagSPD_Matern52: new 3*(5/rho^2 + matern_indices)^3 equals original; factor unchanged. ✓
• PHI/PHI_periodic: the outer product vector * row_vector reproduces diag_post_multiply(rep_matrix(...), linspaced_vector(...)) entry-for-entry, while avoiding materialising the replicate matrix. ✓
• convolve_with_rev_pmf: dropping the if (len > xlen) guard is safe because a Stan for (s in (xlen+1):len) loop runs zero times when xlen+1 > len; the inline comment documents this. Moving vector[len] z below the reject checks is a minor allocation tidy-up. ✓

rstan constraint: matern_indices keeps the linspaced_vector(M, 1, M) bounds as data-only literals and applies the pi()/(2*L) scaling afterwards, which is the fix for the data-only-argument constraint that blocked #1274. PHI/PHI_periodic likewise use literal bounds via linspaced_row_vector(M, 1, M).

Tests: new matern_indices test independently reconstructs the expected vector and checks dimensions/positivity.

Style: new lines all within the 80-char limit; no trailing whitespace. The over-80 lines flagged by grep are all pre-existing.

Minor suggestions (non-blocking):

• The clarified nu reject message (0.5, 1.5, 2.5) now matches the literal values compared against — good.
• Consider a brief note that the convolve zero-iteration loop relies on Stan range semantics; the comment already covers this adequately.

This is an automated review from the PR quality gate. CI is in progress; will relay results.

[seabbs-bot]

Automated re-review — commit 3bd55a5

Reviewed the new commit. No Critical or Important issues.

reporting_phi equivalence — verified exact:

• report_lp: reached only inside if (model_type), so the helper returns inv_square(reporting_overdispersion) — identical to the original inline call.
• report_log_lik: the NB branch is the else (Poisson is the if), so model_type is truthy there → inv_square(...). Identical.
• report_rng: original defaulted phi = 1e5 then overwrote with inv_square(...) iff model_type. The helper model_type ? inv_square(...) : 1e5 reproduces this exactly, including the Poisson fallback.
• model_type is int; nonzero-as-true in the ternary matches the original if (model_type) semantics.

Other:

• convolve.stan: stray ;; → ; in calc_conv_indices_len — pure cleanup.
• New reporting_phi test covers both branches (NB 1/od^2, Poisson 1e5) independently of the helper internals.
• Doc block follows the file's roxygen/@ingroup observation_model convention; new lines within 80 chars.

This is an automated review from the PR quality gate. CI re-running on the new commit; will relay results.

[github-actions]

Thank you for your contribution seabbs-bot 🚀! Your synthetic_recovery markdown is ready for download 👉 here 👈!
_{(The artifact expires on 2026-06-08T17:58:42Z. You can re-generate it by re-running the workflow here.)}

[github-actions]

This is how benchmark results would change (along with a 95% confidence interval in relative change) if 3bd55a5 is merged into main:

• ✔️default: 22.8s -> 22.6s [-12.61%, +10.1%]
• ✔️no_delays: 23.7s -> 27.6s [-14.32%, +47.89%]
• ✔️random_walk: 15.3s -> 9.14s [-106.02%, +25.73%]
• ✔️stationary: 12.6s -> 12.4s [-15.44%, +11.24%]
• ✔️uncertain: 1.03m -> 1.56m [-43.13%, +146.82%]
  Further explanation regarding interpretation and methodology can be found in the documentation.

[seabbs]

Doesn't seem to help runtimes much but is slightly cleaner I think. @sbfnk thoughts?

[sbfnk]

I think the simplification of code itself is enough of a positive.