Bench suite + Stable mode - #8
Merged
Merged
Conversation
DoubleFloats v1.9.x cbrt/∛ on a Double64 returns the raw (hi, lo) component tuple instead of a Double64, breaking every ∛ call reached with a Double64 argument (volume_equivalent_radius, Riccati-Bessel term estimators) and making the EBCM path unusable for Double64. Add a behavior-gated, Double64-specific method that recovers full precision via Newton refinement and disables itself once upstream returns a Double64.
Add a BenchmarkTools SUITE (benchmark/benchmarks.jl) compatible with both PkgBenchmark and AirspeedVelocity, a local runner with tuned params, and a 'just bench' recipe. Covers special functions, EBCM, IITM, post-processing, linearization, and Float64-vs-Double64 precision. Raw result JSONs are gitignored; tuned params.json is committed for reproducible runs.
Run the benchmark SUITE on pull requests touching src/, benchmark/, or Project.toml, comparing the PR head against main back-to-back on the same runner. Uses the pull_request event with a read-only token and job-summary output, so fork PR code never runs with a write token.
…ian slices 𝐓_from_𝐏_and_𝐔 now solves -𝐏/𝐐 via an LU factorization instead of forming inv(𝐐) (~1.3-1.5x faster, numerically cleaner). ∂𝐓_from_𝐏_and_𝐔 uses ∂𝐓 = -(∂𝐏 + 𝐓 ∂𝐐) 𝐐⁻¹ with a single factorization. linearize_transition_matrix factors each m-block's 𝐐 once and reuses it for the value block and every parameter slice instead of recomputing inv per parameter, ~2x faster for the analytic Jacobian. Arblib (Arb/Acb) matrices keep their dedicated inv via a type guard.
…sion The IITM forward solvers (axisymmetric, n-fold, arbitrary) computed inv(𝐈 - wri·𝐔·𝐆)·𝐔 and inv(𝐈 - 𝐓·𝐐ₕₕ)·𝐓 with an explicit inverse plus a matrix product on every radial layer. Replace each with a left division via _iitm_ldiv, which factorizes once and avoids forming the inverse (~1.3-1.5x faster per inner solve, numerically cleaner). Arb/Acb matrices keep their dedicated inv via a type guard. The IITM linearization path already used lu/\.
Add src/EBCM/stabilization.jl implementing F⁺_{nk}(s,x), the non-negative-power part of x·χ_n(x)·ψ_k(sx) (Somerville, Auguié & Le Ru, JQSRT 123:153 (2013), Eq. 45-46), the cancellation-free replacement for χ_n(x)·ψ_k(sx) in the spheroid EBCM U-matrix integrand. Coefficients come from the DLMF 10.53 power series, accumulated in BigFloat to defeat the (2n-1)!!-scale cancellation, then evaluated in the target precision. Implemented from the published equations only (SMARTIES is CC BY-NC; this package is MIT). Validated against the package Bessel product (n≤k+2), the Somerville recurrence Eq.51, and Float64-vs-BigFloat evaluation. Not yet wired into the solver.
ricattibessely! unconditionally wrote χ[2], so ricattibessely(1, x) threw a BoundsError (whereas ricattibesselj(1, …) worked). Guard the χ[2] write behind nₘₐₓ ≥ 2; the χ′ loop is already empty for nₘₐₓ=1. Adds a regression testitem.
… U integrand (S2) Add the F⁺-based primitives the spheroid EBCM U-matrix integrand needs (Somerville 2013, Eqs. 59-62): [xχ_nψ′_k]⁺, [xχ′_nψ_k]⁺, and the L⁷/L⁸ radial factors [x(χ′_nψ′_k + n(n+1)χ_nψ_k/(sx²))]⁺ and its k(k+1) partner. Each is built from F⁺ at shifted indices via the Riccati-Bessel recurrences. Validated against the direct Bessel products for n ≤ k (rel ≤ 1e-22). Next: assemble the m=0 U blocks from these and check against the Arb ground truth.
…(S2b)
The bare F⁺ power series (S1) is only well conditioned where the multipole
order exceeds the argument, so it cannot fill the small-n/large-x corner of
the matrix. Following Somerville, Auguié & Le Ru, JQSRT 123 (2013) §4.2:
- _eval_F⁺ now truncates the series once it has started to decay and a term
is negligible, which also prevents x^{2·nterms} from overflowing to NaN for
large x (previously NaN for x ≳ 25).
- _F⁺_matrix builds the whole F⁺_{nk} matrix (0 ≤ n,k ≤ N+1, n+k even) for a
single x by combining: direct Riccati-Bessel products for n ≤ k+2, the power
series for the last row n=N+1 (accurate since N+1 ≫ x), and the stable Eq.51
scheme-(c) recursion for the n ≥ k+4 cancellation block.
Validated against high-precision BigFloat series over the entire bottom-left
block for x = 0.1 … 35: max relative error ≤ 3e-14, no entry above 1e-8.
…0 (S3) transition_matrix gains a stable kwarg (default false). stable=true is gated to Spheroid (the cancellation-free F⁺ integrands rely on the spheroid surface making the divergent Laurent terms integrate to zero) and is threaded through transition_matrix_m₀/_m → ebcm_matrices_m₀/_m → _transition_matrix_m_core. - m>0: the K¹/K² integrands (Somerville Eqs 53-54, n+k odd) use _xχψ′⁺_mat / _xχ′ψ⁺_mat; the L blocks (n+k even) reuse _Fp / _L⁷⁺_mat / _L⁸⁺_mat as for m=0. Diagonal and 𝐏 keep the standard direct products. - The per-quadrature-point F⁺ matrices are m-independent, so they are built once in the m=0 block and carried through the reuse cache to all m>0 blocks. Literature anchor (Somerville, Auguié & Le Ru, JQSRT 123 (2013), Table 2, Model 1 of their Ref. [35]): prolate s=1.55+0.01i, h=4, xmax=k·c=10.079368. With stable=true (N=32, Ng=400) Qsca=3.2129055420 and Qext=3.3672129262 match the published New DP/AP values to ~1e-11; the standard double-precision path degrades to 3.2298 at the same N. Added as a testitem. stable=false leaves all existing results unchanged (full suite green).
…table (S4) - _F⁺_lastrow precomputes the x-independent BigFloat last-row series coefficients once per block; _F⁺_matrix(…; lastrow=) reuses them across every quadrature point instead of recomputing the (expensive) coefficients per point. The stabilized spheroid path is now ≈2.4× the standard assembly (was far worse). - benchmark/benchmarks.jl: add highaspect_std_n24 / highaspect_stable_n24 to the ebcm group so the stabilization cost is tracked. - docs: document the stable kwarg in the transition_matrix docstring and add a 'High-aspect-ratio spheroids' section to docs/src/usage.md (spheroid-only, 2-3× cost, requires nₘₐₓ+1 ≳ k·c + 15). Case C unchanged after the refactor (Qsca/Qext match Somerville to ~1e-11); full suite green.
|
Note Reviews pausedIt looks like this branch is under active development. To avoid overwhelming you with review comments due to an influx of new commits, CodeRabbit has automatically paused this review. You can configure this behavior by changing the Use the following commands to manage reviews:
Use the checkboxes below for quick actions:
📝 WalkthroughWalkthroughAdds benchmark CI and a local benchmark suite; integrates stabilized EBCM (F⁺) and a Sh-matrix precompute/assembly path; replaces explicit matrix inverses in IITM with ldiv helpers; introduces solver abstractions and an iterative driver; refactors linearization internals; updates docs, examples, tests, and compat shims. ChangesBenchmark Infrastructure
Numerical Stability and Linear Algebra
Linearization and Performance
Code Quality, Compat, Tests, Docs, Examples
Estimated code review effort:
✨ Finishing Touches🧪 Generate unit tests (beta)
|
|
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #8 +/- ##
==========================================
+ Coverage 95.20% 95.63% +0.43%
==========================================
Files 29 32 +3
Lines 3732 4403 +671
==========================================
+ Hits 3553 4211 +658
- Misses 179 192 +13
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
![coderabbitai[bot]](https://avatars.githubusercontent.com/in/347564?s=60&v=4){kind=small}
There was a problem hiding this comment.
Actionable comments posted: 1
🧹 Nitpick comments (1)
benchmark/runbenchmarks.jl (1)
15-22: 💤 Low valueClarify the parameter-loading message.
The
@infoat line 16 says "Loaded tuned parameters", but the actual loading happens inbenchmarks.jl(lines 148–152) when that file is included. This check is just confirming the file exists before deciding whether to tune. The message could mislead readers trying to understand the control flow.💡 Suggested clarification
if isfile(PARAMS_PATH) - `@info` "Loaded tuned parameters" path=PARAMS_PATH + `@info` "Using cached tuned parameters" path=PARAMS_PATH else `@info` "Tuning benchmark parameters (first run; this is slow)…"or add a comment:
+# Note: benchmarks.jl already loaded params.json if present; this check +# determines whether we need to tune and save them. if isfile(PARAMS_PATH) `@info` "Loaded tuned parameters" path=PARAMS_PATH🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@benchmark/runbenchmarks.jl` around lines 15 - 22, The log message "Loaded tuned parameters" is misleading because this block only checks PARAMS_PATH existence and triggers tuning (tune!) or saving (BenchmarkTools.save(params(SUITE))); update the `@info` text to reflect that it only detected an existing parameters file (e.g., "Found parameters file; will load when benchmarks.jl is included") or add a concise comment explaining that actual loading happens later in benchmarks.jl, referencing PARAMS_PATH, SUITE, tune!, BenchmarkTools.save, and params to locate the logic.
🤖 Prompt for all review comments with AI agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
Inline comments:
In @.github/workflows/Benchmark.yml:
- Around line 45-48: Replace the mutable tag usage of the GitHub Action
MilesCranmer/AirspeedVelocity.jl@action-v1 with the specific commit SHA
315c11b51ceee8ebd6063d70cff6ae499a040d28 to pin the action; locate the line
containing "uses: MilesCranmer/AirspeedVelocity.jl@action-v1" in the
Benchmark.yml workflow and update it to reference the commit SHA instead of the
tag (keeping a comment like "# action-v1" optional for context).
---
Nitpick comments:
In `@benchmark/runbenchmarks.jl`:
- Around line 15-22: The log message "Loaded tuned parameters" is misleading
because this block only checks PARAMS_PATH existence and triggers tuning (tune!)
or saving (BenchmarkTools.save(params(SUITE))); update the `@info` text to reflect
that it only detected an existing parameters file (e.g., "Found parameters file;
will load when benchmarks.jl is included") or add a concise comment explaining
that actual loading happens later in benchmarks.jl, referencing PARAMS_PATH,
SUITE, tune!, BenchmarkTools.save, and params to locate the logic.
🪄 Autofix (Beta)
Fix all unresolved CodeRabbit comments on this PR:
- Push a commit to this branch (recommended)
- Create a new PR with the fixes
ℹ️ Review info
⚙️ Run configuration
Configuration used: defaults
Review profile: CHILL
Plan: Pro Plus
Run ID: 636b2c90-cf9b-4731-a13a-4def589f68d5
📒 Files selected for processing (40)
.github/workflows/Benchmark.yml.github/workflows/README.md.gitignorebenchmark/Project.tomlbenchmark/README.mdbenchmark/benchmarks.jlbenchmark/params.jsonbenchmark/runbenchmarks.jldocs/make.jldocs/src/usage.mdjustfilepackages/EBCMPrecisionLossEstimators/src/EBCMPrecisionLossEstimators.jlsrc/EBCM/axisymmetric.jlsrc/EBCM/index.jlsrc/EBCM/linearization.jlsrc/EBCM/routines.jlsrc/EBCM/stabilization.jlsrc/IITM/arbitrary.jlsrc/IITM/axisymmetric.jlsrc/IITM/fourier.jlsrc/IITM/linearization.jlsrc/IITM/nfold.jlsrc/Mie/MieTransitionMatrix.jlsrc/Mie/bhcoat.jlsrc/Mie/linearization.jlsrc/common/AbstractTransitionMatrix.jlsrc/common/AxisymmetricTransitionMatrix.jlsrc/common/RandomOrientationTransitionMatrix.jlsrc/compat/index.jlsrc/linearization.jlsrc/shapes/chebyshev.jlsrc/shapes/cylinder.jlsrc/shapes/prism.jlsrc/shapes/spheroid.jlsrc/special_functions/bessel.jlsrc/special_functions/quadrature.jlsrc/special_functions/wignerd.jltest/compat.jltest/linearization.jltest/runtests.jl
Comment on lines
+45
to
+48
| - uses: MilesCranmer/AirspeedVelocity.jl@action-v1 | ||
| with: | ||
| julia-version: "1" | ||
| job-summary: "true" |
There was a problem hiding this comment.
🧩 Analysis chain
🏁 Script executed:
#!/bin/bash
# Fetch the commit SHA behind the action-v1 tag
gh api repos/MilesCranmer/AirspeedVelocity.jl/git/ref/tags/action-v1 --jq '.object.sha' 2>/dev/null || \
gh api repos/MilesCranmer/AirspeedVelocity.jl/commits/action-v1 --jq '.sha'Repository: JuliaRemoteSensing/TransitionMatrices.jl
Length of output: 125
🏁 Script executed:
#!/bin/bash
set -euo pipefail
echo "== Benchmark.yml uses lines around 40-60 =="
nl -ba .github/workflows/Benchmark.yml | sed -n '35,70p'
echo
echo "== Any other unpinned uses: lines (tags, no sha) in the workflow =="
rg -n '^\s*uses:\s*.*@[^#\s]+' .github/workflows/Benchmark.yml || true
echo
echo "== Specifically check for action-v1 in the workflow =="
rg -n 'action-v1' .github/workflows/Benchmark.yml || trueRepository: JuliaRemoteSensing/TransitionMatrices.jl
Length of output: 189
🏁 Script executed:
#!/bin/bash
set -euo pipefail
echo "== Benchmark.yml uses lines around 40-60 =="
cat -n .github/workflows/Benchmark.yml | sed -n '35,70p'
echo
echo "== Any other unpinned uses: lines (tags, no sha) in the workflow =="
rg -n '^\s*uses:\s*.*@[^#\s]+' .github/workflows/Benchmark.yml || true
echo
echo "== Specifically check for action-v1 in the workflow =="
rg -n 'action-v1' .github/workflows/Benchmark.yml || trueRepository: JuliaRemoteSensing/TransitionMatrices.jl
Length of output: 929
Pin the AirspeedVelocity action to a commit SHA
MilesCranmer/AirspeedVelocity.jl@action-v1 is a mutable tag; pin it to commit 315c11b51ceee8ebd6063d70cff6ae499a040d28 to prevent supply-chain mutation.
Replace:
- uses: MilesCranmer/AirspeedVelocity.jl@action-v1With:
- uses: MilesCranmer/AirspeedVelocity.jl@315c11b51ceee8ebd6063d70cff6ae499a040d28 # action-v1🧰 Tools
🪛 zizmor (1.25.2)
[error] 45-45: unpinned action reference (unpinned-uses): action is not pinned to a hash (required by blanket policy)
(unpinned-uses)
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
In @.github/workflows/Benchmark.yml around lines 45 - 48, Replace the mutable
tag usage of the GitHub Action MilesCranmer/AirspeedVelocity.jl@action-v1 with
the specific commit SHA 315c11b51ceee8ebd6063d70cff6ae499a040d28 to pin the
action; locate the line containing "uses:
MilesCranmer/AirspeedVelocity.jl@action-v1" in the Benchmark.yml workflow and
update it to reference the commit SHA instead of the tag (keeping a comment like
"# action-v1" optional for context).
Add src/EBCM/shmatrix.jl: the Sh-matrix (Farafonov-family) moment-separation
machinery. Every EBCM surface integral factors as
∫ = Σ_terms [size/material coeff(k,s)] × [shape-only moment M(q)],
so the geometry moments are computed once and each (λ,s) evaluation is a cheap
coefficient×moment sum — the basis for fast parameter sweeps. The same term-by-
term expand-and-integrate also recovers the spheroid U stabilization: the
negative-r-power moments vanish analytically for a spheroid, so computed at
high precision they contribute ≈0 and the χ_n·ψ_{n′} cancellation never forms.
This commit lands the m=0 block:
- _sh_series: radial power-series term tables (reuse _β,_γ; DLMF 10.53)
- _sh_recon: coefficient×moment reconstruction (incremental powers)
- _sh_moments_m: shape-only moment families (τd, dτ, πd, ππττ), BigFloat-
accumulated, stored Float64; high-precision shape widening for spheroids
- _sh_matrices_m₀: m=0 P,U assembly mirroring ebcm_matrices_m₀ term by term
Validated vs ebcm_matrices_m₀ (aspect-2 spheroid, λ=2π): P 1.6e-14, U 2.3e-12.
On an aspect-4 prolate at nmax=24 the std Float64 U loses to 6.4e-4 while the
Sh reconstruction holds 7.2e-8 vs BigFloat ground truth.
Add _sh_matrices_m: the m>0 P,U assembly from moments, mirroring
_transition_matrix_m_core term by term. K-blocks (P₁₂,P₂₁,U₁₂,U₂₁) use the
M^{πd} moment family; L-blocks reuse the same reconstruction as m=0 with the
m-dependent Wigner functions baked into the moment tables. Fix _sh_moments_m to
index the Wigner arrays over max(0,m):nmax as the recursion requires.
Validated vs ebcm_matrices_m (aspect-2 spheroid, λ=2π, m=1..4): P ≲2e-14,
U ≲6e-12.
… spectrum Add the public moment-separation API: - ShPreparation: precomputed shape-only moment tables (m=0:nmax). The moments depend only on geometry, not on λ or mᵣ, so one preparation feeds an arbitrarily long sweep cheaply. Records `stable` (spheroid-only U). - prepare_sh(shape, nmax, Ng; B, momtype, store): one-time geometry quadrature. - transition_matrix(prep, λ, mᵣ): assemble the full T from the moments — the cheap inner call of a sweep. 2-arg form reuses the prepared shape's mᵣ. - transition_matrix_spectrum(prep, λs, mᵣs): sweep helper (scalar mᵣ held fixed, or a per-λ dispersion table). Export prepare_sh, transition_matrix_spectrum, ShPreparation. Validated: - aspect-2 spheroid, full T vs package: Qsca/Qext rel err ~1e-15 - aspect-4 prolate vs package stable=true: Csca rel err 3.7e-12 - 3-point (λ,mᵣ) sweep reusing one preparation: Csca rel err ~1e-15 each
- 3 @testitem blocks in shmatrix.jl: block-level reconstruction vs
ebcm_matrices_m₀/ebcm_matrices_m (m=0 and m=1..4); full transition_matrix
Qsca/Qext vs the package + high-aspect vs stable=true; spectrum sweep vs
per-point assembly (and scalar-index broadcast). Full suite: 47676 pass.
- benchmark/benchmarks.jl: SUITE["sh_matrix"] = {prepare_n8, sweep_point_n8,
fresh_assembly_n8}. Measured 6.6× per-point speedup at nmax=8,Ng=32 (the
margin grows with Ng since prepare_sh absorbs the quadrature).
- docs/src/usage.md: "Parameter sweeps with the Sh-matrix method" section
(prepare_sh / transition_matrix(prep,λ,mᵣ) / transition_matrix_spectrum,
stability note, B / size-parameter caveat, spheroid-only U gating). The new
exported docstrings are auto-listed by api.md's @autodocs.
Two size/material-independent costs were being repeated at every sweep point. Hoist them and fold the reconstruction: - _sh_coeff_tables: the radial series coefficients (BigInt _β/_γ rational arithmetic → R) depend only on order and B, so build them once in prepare_sh and store in ShPreparation. This was the dominant per-point cost. - _sh_weighted: per (λ,mᵣ) point, fold z^(base+2b) into the coefficients once per order (incremental powers, no ^), reused across all integrands and m. - _sh_conv: replace _sh_recon's double loop with an anti-diagonal fold by j=b+c. The r-power q₀+2j depends only on j, so the moment is fetched O(B) times instead of O(B²). Direct convolution (not FFT): FFT's global roundoff mixing would recontaminate the well-conditioned terms with the ill-conditioned (huge-coefficient, vanishing-moment) ones — the very cancellation this method avoids. _sh_matrices_m₀/_sh_matrices_m now take the coefficient tables and drop the B argument (derived from the tables). Correctness unchanged (Csca rel err ~1e-15; high-aspect vs stable ~6e-10). Full suite 47676 pass. Per-point cost 10861µs → 610µs; sweep speedup vs fresh assembly 6.6× → 119× at nmax=8,Ng=32 (grows with Ng).
Reduce the one-time prepare_sh precompute cost with three structural changes that do not touch precision: - _sh_geometry: the Gauss–Legendre nodes/weights, ϑ, and the r^q weight tables (WR = w·r′·r^q, Wq = w·r^q) are m-independent, so compute them once instead of rebuilding inside every one of the nmax+1 per-m moment passes. Splitting the weight into WR (with r′) and Wq (no r′) also removes the per-term r′ multiply. - Drop the Mdτ family: Mdτ[n,n′,q] = Mτd[n′,n,q] exactly, so store only Mτd and read it with the orders swapped in the assembly. ~1/3 less pair-moment work and storage. - Hoist the angular product τ[i,n]·d[i,n′] (and π·d) out of the q loop — formed once per (n,n′) instead of once per q. Correctness unchanged (Csca rel err ~1e-15; high-aspect vs stable ~6e-10). Full suite 47676 pass (2m33s, down from 3m05s). Warm prepare_sh: nmax=12 0.65s, nmax=16 2.1s, nmax=24 12.6s.
The catastrophic cancellation that forces high-precision accumulation lives entirely in the negative-r-power (q<0) moments — for a spheroid those vanish analytically and must reach ≈0, not round-off. The q≥0 bulk is well conditioned and never needed the extra precision. Split the contraction by q-range: _sh_geometry also returns the weight tables cast to the fast type Tf = store<:IEEEFloat ? store : momtype (WRf/Wqf), and _sh_moments_m contracts q≥0 in Tf (hardware float, @simd, no heap allocation) while only the q<0 slice keeps the slow momtype accumulation. The angular functions get a one-time Tf cast per m for the bulk. Results unchanged (Csca rel err ~4e-15; high-aspect aspect-4 vs stable 5.0e-10, i.e. still the ~1e-9 stable level). Full suite 47676 pass. Warm prepare_sh ~4.6× faster: nmax=12 0.65→0.148s, nmax=16 2.1→0.455s, nmax=24 12.6→2.7s.
The remaining precompute cost was ~99% the BigFloat q<0 contraction. For a spheroid the negative-r-power moments vanish analytically (the F⁺ result), so they can be set to zero instead of integrated, and the whole precompute then runs in the shape's hardware precision with no BigFloat at all. prepare_sh now auto-selects: spheroid → momtype = shape's real type and zero the q≤-2 moments (compute only the well-conditioned q≥-1 in fast precision); other shapes (negative moments don't vanish) → BigFloat as before. _sh_moments_m gains a vanish_negative flag that skips the q<0 integration. momtype=BigFloat can be passed to force the high-precision path for a spheroid. Validated against the BigFloat full-EBCM ground truth across aspect 2–6 and s→1: observables match the package to 1e-9..1e-16, and at aspect 6 the T-block accuracy (4.1e-7) equals the existing stable=true path (2.8e-7) — the limit is the shared reconstruction-series truncation, not the zeroing. Full suite 47676 pass. Warm prepare_sh: nmax=16 2.1s→6.9ms, nmax=24 12.6s→19ms (~300–650×).
BREAKING (0.5): reorganize the axisymmetric T-matrix solvers into two clearly
separated layers, selected by a solver object passed to transition_matrix(s, λ, solver):
- Fixed solution (explicit discretization, one build):
EBCM(nₘₐₓ, Ng; stable), IITM(nₘₐₓ, Nr, Nϑ[, Nφ]; rₘᵢₙ), ShMatrix(nₘₐₓ, Ng; B, …)
- Iterative solver (sweeps discretization to Qsca/Qext convergence):
Iterative(EBCM; stable, threshold, ndgs, maxiter, …)
New file src/solvers.jl: AbstractSolver/AbstractFixedSolver, the fixed structs,
ConvergencePolicy, Iterative{EBCM}, a method-generic convergence driver with a
four-hook interface (initial/probe/step/assemble), the EBCM cheap m=0 probe
specialization, and a full-T cross-section fallback (seam for Iterative(IITM)).
Fixed solvers are thin shims over the existing functions (unchanged).
Closes the main gap: `Iterative(EBCM; stable=true)` gives an auto-converged AND
stabilized high-aspect spheroid — previously impossible (stable was fixed-order
only). The iterative path guards stable to Spheroid (the per-m functions, called
directly here, lack the 4-arg's guard).
The bare transition_matrix(s, λ) / calc_T(s, λ) is unchanged behaviorally
(= Iterative(EBCM)) but now routes through the unified driver. The old kwargs
auto-converge transition_matrix(s, λ; threshold, ndgs, …) is REMOVED; its kwargs
move to Iterative(EBCM; …) / ConvergencePolicy. routine_mishchenko and the m=0
efficiency probes are reused unchanged.
Export AbstractSolver, AbstractFixedSolver, EBCM, IITM, ShMatrix, Iterative,
ConvergencePolicy. Docs (usage.md) rewritten to the solver-type API. Full suite
47686 pass.
Deferred: Iterative(IITM) (axisymmetric needs the full-T probe, stubbed),
arbitrary/nfold IITM iterative (needs cross sections for TransitionMatrix),
Sh-matrix convergence.
Add a @testitem confirming the Sh-matrix backend works for non-spheroid axisymmetric shapes (Cylinder, Chebyshev), not just spheroids: m=0 P and U blocks match ebcm_matrices_m₀ to <1e-10, and the full-T Csca/Cext match the classic EBCM build to ~1e-9, at moderate size. Documents that prep.stable=false for non-spheroids (the analytic-zeroing fast precompute is spheroid-only) while the moment machinery, P, U, and sweep reuse are fully general. Full suite 47696 pass.
Add an examples/ folder with its own environment (Project.toml: Pluto, Plots, ForwardDiff, BenchmarkTools, dev'd TransitionMatrices) and a Pluto notebook showcasing fast wavelength/refractive-index sensitivities: prepare_sh (geometry quadrature once) → ForwardDiff over the cheap reconstruction → exact ∂Csca/∂λ and ∂Csca/∂mᵣ across a spectrum, at a large speedup over differentiating the from-scratch assembly. The notebook plots the Csca/Cext spectra and their sensitivities and times the Sh path vs classic AD. Generated through Pluto's own serializer (genuine cell IDs) and verified end-to-end headless (all 11 cells run, 0 errors). examples/ Manifest.toml stays gitignored; the notebook's first cell activates the env and develops the parent so it runs from a fresh checkout. README documents usage.
Add verified Pluto notebooks (each generated via Pluto's serializer with genuine cell IDs, run headless with 0 errors): - shapes_gallery.jl — define Spheroid/Cylinder/Chebyshev/Prism, compute T and Qsca/Qext/g/ω, bar chart, solver-choice notes (Prism via IITM). - angular_scattering.jl — orientation-averaged scattering matrix vs angle (phase function F₁₁, polarization -F₁₂/F₁₁, F₂₂/F₁₁, F₄₄/F₁₁) + fixed-geometry amplitude/phase matrices. - orientation_averaging.jl — analytic RandomOrientationTransitionMatrix vs numerical orientation_average, with a convergence plot. - solver_landscape.jl — the two-layer API (EBCM/IITM/Iterative/stable); classic EBCM diverges with nₘₐₓ on a high-aspect spheroid while stable holds. README lists all notebooks. wip: post showcase notebooks
There was a problem hiding this comment.
Actionable comments posted: 2
🧹 Nitpick comments (1)
examples/spectral_sensitivity.jl (1)
110-115: API compatibility looks intentional; optionally align with the solver-object form.
transition_matrix(s, λ, nₘₐₓ, Ng; stable=false)is still defined (src/EBCM/axisymmetric.jl), andtransition_matrix(s, λ, EBCM(nₘₐₓ, Ng))is a thin wrapper over the same 4-arg implementation (src/solvers.jl), with tests assertingtransition_matrix(s, λ, EBCM(8, 64)).𝐓 == transition_matrix(s, λ, 8, 64).𝐓(and forstable=true). So the old signature is intentionally maintained; updating the example call totransition_matrix(s, λ, EBCM(nmax, Ng))would be for consistency only.🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@examples/spectral_sensitivity.jl` around lines 110 - 115, The example uses the 4-argument transition_matrix call but the project also exposes the solver-object wrapper; to make the example consistent with the solver-based API, change the call inside Csca_classic to pass an EBCM solver object (use transition_matrix(s, λ, EBCM(nmax, Ng))) instead of transition_matrix(s, λ, nmax, Ng); update the Csca_classic function definition accordingly so calc_Csca receives transition_matrix(s, λ, EBCM(nmax, Ng)) and leave the rest (s construction and calc_Csca call) unchanged.
🤖 Prompt for all review comments with AI agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
Inline comments:
In `@src/EBCM/shmatrix.jl`:
- Around line 621-625: The assertion that Ng is even in prepare_sh is too
strict; change it to only enforce evenness when the shape uses the symmetry
shortcut by calling has_symmetric_plane(shape) — i.e., replace the unconditional
`@assert` iseven(Ng) with a conditional check that runs `@assert` iseven(Ng) "Ng
must be even!" only if has_symmetric_plane(shape) is true (this aligns
prepare_sh with _sh_geometry which folds quadrature to Ng ÷ 2 only for symmetric
shapes).
- Around line 681-684: The function transition_matrix_spectrum currently uses
zip(λs, mᵣs) which silently truncates to the shorter iterable; instead, validate
inputs and reject mismatched lengths: if mᵣs is a Number behave as before,
otherwise check that length(mᵣs) == length(λs) and throw an ArgumentError (or
similar) if they differ; keep the rest of the implementation (the comprehension
calling transition_matrix(prep, λ, m)) unchanged so you still construct the full
spectrum when lengths match.
---
Nitpick comments:
In `@examples/spectral_sensitivity.jl`:
- Around line 110-115: The example uses the 4-argument transition_matrix call
but the project also exposes the solver-object wrapper; to make the example
consistent with the solver-based API, change the call inside Csca_classic to
pass an EBCM solver object (use transition_matrix(s, λ, EBCM(nmax, Ng))) instead
of transition_matrix(s, λ, nmax, Ng); update the Csca_classic function
definition accordingly so calc_Csca receives transition_matrix(s, λ, EBCM(nmax,
Ng)) and leave the rest (s construction and calc_Csca call) unchanged.
🪄 Autofix (Beta)
Fix all unresolved CodeRabbit comments on this PR:
- Push a commit to this branch (recommended)
- Create a new PR with the fixes
ℹ️ Review info
⚙️ Run configuration
Configuration used: defaults
Review profile: CHILL
Plan: Pro Plus
Run ID: 50b01f38-063c-4bf4-9989-cf444099f2bc
📒 Files selected for processing (42)
benchmark/benchmarks.jldocs/make.jldocs/src/usage.mdexamples/Project.tomlexamples/README.mdexamples/angular_scattering.jlexamples/orientation_averaging.jlexamples/shapes_gallery.jlexamples/solver_landscape.jlexamples/spectral_sensitivity.jlpackages/EBCMPrecisionLossEstimators/src/EBCMPrecisionLossEstimators.jlsrc/EBCM/axisymmetric.jlsrc/EBCM/index.jlsrc/EBCM/linearization.jlsrc/EBCM/routines.jlsrc/EBCM/shmatrix.jlsrc/EBCM/stabilization.jlsrc/IITM/arbitrary.jlsrc/IITM/axisymmetric.jlsrc/IITM/fourier.jlsrc/IITM/linearization.jlsrc/IITM/nfold.jlsrc/Mie/MieTransitionMatrix.jlsrc/Mie/bhcoat.jlsrc/Mie/linearization.jlsrc/TransitionMatrices.jlsrc/common/AbstractTransitionMatrix.jlsrc/common/AxisymmetricTransitionMatrix.jlsrc/common/RandomOrientationTransitionMatrix.jlsrc/compat/index.jlsrc/linearization.jlsrc/shapes/chebyshev.jlsrc/shapes/cylinder.jlsrc/shapes/prism.jlsrc/shapes/spheroid.jlsrc/solvers.jlsrc/special_functions/bessel.jlsrc/special_functions/quadrature.jlsrc/special_functions/wignerd.jltest/compat.jltest/linearization.jltest/runtests.jl
💤 Files with no reviewable changes (1)
- src/EBCM/routines.jl
✅ Files skipped from review due to trivial changes (15)
- src/shapes/cylinder.jl
- examples/README.md
- test/runtests.jl
- src/shapes/prism.jl
- src/Mie/bhcoat.jl
- src/special_functions/quadrature.jl
- docs/make.jl
- src/EBCM/index.jl
- src/common/AxisymmetricTransitionMatrix.jl
- src/common/RandomOrientationTransitionMatrix.jl
- src/common/AbstractTransitionMatrix.jl
- src/special_functions/wignerd.jl
- src/Mie/linearization.jl
- src/Mie/MieTransitionMatrix.jl
- test/linearization.jl
🚧 Files skipped from review as they are similar to previous changes (13)
- src/shapes/spheroid.jl
- packages/EBCMPrecisionLossEstimators/src/EBCMPrecisionLossEstimators.jl
- src/compat/index.jl
- src/IITM/fourier.jl
- src/special_functions/bessel.jl
- src/IITM/nfold.jl
- src/linearization.jl
- src/IITM/axisymmetric.jl
- src/IITM/arbitrary.jl
- src/IITM/linearization.jl
- src/EBCM/linearization.jl
- src/EBCM/stabilization.jl
- src/EBCM/axisymmetric.jl
Two things, both about making the Documenter build clean and the API page useful: 1. Resolve the CrossReferences failure: several docstrings linked [`name`](@ref) to internal helpers that are not rendered — `routine_mishchenko` and the five `_sh_*` helpers. Demote those to plain code formatting. 2. Filter internals out of the API reference: api.md's @autodocs now excludes underscore-prefixed bindings (Filter = t -> !startswith(string(nameof(t)), "_")), and makedocs uses checkdocs = :exported so the now-omitted internals don't trip the missing_docs check. Docs build is green; api.html drops 150.8 → 119.8 KiB (the remainder is the legitimate public API; over the 100 KiB warn limit but under the 200 KiB hard limit).
Integrate the examples/ Pluto notebooks into the Documenter build so users see the executed output (plots, tables, timings) on the docs site, under a new "Examples" section. - docs/make.jl: before makedocs, run PlutoStaticHTML.build_notebooks over examples/ (documenter_output) and copy the generated .md under docs/src/examples/. Each notebook self-activates the examples env in its own Pluto worker, so it runs without disturbing the docs process. ENV["GKSwstype"] = "100" for headless GR; SKIP_EXAMPLES=true bypasses rendering for fast local builds; previous_dir caching reuses unchanged renders. - Add the rendered pages to the nav and exempt them from Documenter's HTML size_threshold (embedded plots are large). - Add PlutoStaticHTML to docs/Project.toml. - spectral_sensitivity.jl: drop the @belapsed timing cell to samples=1 so the render isn't dominated by benchmarking. - .gitignore the generated example .md (examples/*.md except README, and docs/src/examples/). Verified: full docs build is green; all five notebooks render with their plots embedded as base64 SVG (e.g. angular_scattering.html ≈ 415 KB), Examples section present, no cross-reference / missing_docs / size errors.
Two review fixes in src/EBCM/shmatrix.jl plus an example tidy: - prepare_sh: the unconditional `@assert iseven(Ng)` was too strict. Only the equatorial-symmetry shortcut folds the quadrature to Ng÷2 (see _sh_geometry), so require even Ng only when has_symmetric_plane(shape); non-symmetric shapes (e.g. odd-order Chebyshev) now accept odd Ng. - transition_matrix_spectrum: zip(λs, mᵣs) silently truncated to the shorter iterable. Now, when mᵣs is not a Number, require length(mᵣs) == length(λs) and throw ArgumentError otherwise; scalar mᵣs still broadcasts. - examples/spectral_sensitivity.jl: use the solver-object form transition_matrix(s, λ, EBCM(nmax, Ng)) in Csca_classic for API consistency. Validated: odd Ng on a non-symmetric Chebyshev works (finite Csca), symmetric+odd still throws; mismatched spectrum lengths throw, matched/scalar work; full suite 47696 pass; the notebook re-renders headless with 0 errors.
…ased Add CHANGELOG.md (Keep a Changelog format) with an entry and compare/release link for every release tag (v0.1.0, v0.2.0, v0.3.0, v0.3.1, v0.4.0). The post-v0.4.0 work (F⁺ stable mode, Sh-matrix method, two-layer solver API [breaking], example notebooks, docs, perf and fixes) is recorded under [Unreleased] on top of v0.4.0, flagged as targeting v0.5.0 due to the breaking solver-API change.
The docs CI failed loading Plots/GR on the Linux runner (GR_jll → GLFW → Pango, undefined HarfBuzz symbol). Plots 1.41.6 hardcodes GR as the default backend, so `using Plots` always loads it — there is no in-Plots way to dodge it. Instead, stop rendering the notebooks in CI: - Commit the rendered pages (docs/src/examples/*.md, plots embedded as base64 SVG) and have makedocs use them as-is; track that directory in git. - make.jl renders only when BUILD_EXAMPLES=true; the PlutoStaticHTML/Pluto load and the headless-GR setting move inside that branch, so the default (CI) build loads neither — verified: a `just docs` build is green with PlutoStaticHTML and Pluto not loaded and all five example pages built from the committed sources. - Justfile: add `docs-examples` to re-render + rebuild locally (commit the regenerated pages afterward). `docs`/`ci` leave BUILD_EXAMPLES unset. PlutoStaticHTML stays in docs/Project.toml (needed by `docs-examples`) but is no longer loaded by the CI docs build.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary by CodeRabbit
New Features
Bug Fixes
Tests
Documentation