Fix correctness and performance defects found by codebase review - #32
Merged
Merged
Conversation
Two rounds of review over the branch diff and then over the files the branch never touched. Every finding below was reproduced before fixing and re-verified after. Breaking (see CHANGELOG.md for porting instructions): - analysis_axisym / analysis_axisym_l omitted the phi quadrature factor cphi*nlon = 2pi, so every coefficient was 1/2pi too small and they inverted neither synthesis_axisym nor the m=0 column of analysis. - dist_synthesis! / dist_synthesis_sphtor! keep rejecting real_output=false with real-backed output arrays. real_output=false now computes a different function (m>=0 half, no Hermitian mirror), so no tolerance can accept the old usage; the error messages now name the porting path. Correctness: - Batch QST and sphtor transforms derived their output eltype from a single input instead of promoting across all of them, silently truncating double-precision components (2.05e-8 error, now 1.39e-17). - ForwardDiff could not flow through any plan-based batch transform. SHTPlan is FFTW-backed and cannot hold Dual, so the batch entry points now route non-FFTW eltypes through the plan-free cfg-form transforms. - load_config silently rebuilt :driscoll_healy grids as :regular, changing theta nodes and weights (8.7e-16 -> 3e-3 error). save_config now records use_dh_weights. - shtns_set_grid returned success for unrecognized grid codes while producing a mirrored grid; the fallback branch was missing reverse!. - Pole-inclusive grids with nlat == 1 built an all-NaN config from pi/0 and returned NaN from every later transform without raising. - dist_SH_mul_mx! crashed on every mres > 1 config (missing m % mres guard). - dist_SH_Yrotate crashed on mres > 1; a Y-rotation mixes orders and is not representable in an mres-strided layout, so it now says so up front. - device_transfer_arrays rejected :cuda/:amdgpu, the config vocabulary. - Eleven comments left by the orthonormal refactor prescribed conversions the code no longer performs, including the AD extension's only statement of the rrule/primal normalization contract, which stated its inverse. Performance: - Distributed analysis went from three collectives per call to one on first use and none after: the phi_is_local_all / theta_is_distributed predicates are reduced together in a single Allreduce and cached per (pencil, communicator), with a rank-symmetric eviction policy. - dist_synthesis_packed_cplx is single-pass, down from two full distributed syntheses; it now matches the serial reference exactly. - prepare_plm_tables! built NP/NdP tables bit-for-bit identical to plm/dplm. They now alias, halving table memory and build time, and estimate_table_memory finally matches reality (it reported half). - Cached FFT plans are built UNALIGNED, removing a silent fallback to an O(n^2) DFT; batch FFT helpers reuse the cache instead of re-planning. Internal and tests: - pack_lm!/pack_lm/unpack_lm!/unpack_lm in src/layout.jl replace six open-coded copies of the packed<->dense (l,m) mapping. The m % mres guard had needed fixing three separate times across those copies. - The axisym round-trip test fitted and divided out its own scale factor, so it passed while every coefficient was 1/2pi off. It now asserts the absolute identity and pins axisym analysis to the m=0 column of analysis. Verified: serial suite 67279/67279; MPI comprehensive and audit-fix suites pass at 2 and 4 ranks; distributed analysis matches serial to 2.2e-16 on 1/2/4 ranks for both decompositions. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…actly CI failed on x86_64 at test_plan.jl:158 with a 2e-16 relative difference, while the same assertion passed on arm64. The cause is this branch's own UNALIGNED change: `synthesis` reaches its phi transform through the shared plan cache, which is now built UNALIGNED because it reuses plans across arbitrary caller arrays, while `SHTPlan` plans its own stably-aligned buffers with default flags. FFTW can pick different codelets for the two, so results may differ in the last ulp. That is legitimate floating point, not a convention bug, and it varies by CPU and FFTW build — so bit-exactness was never a property to assert. Keeping the plan's own flags as they are: SHTPlan owns buffers of stable alignment, and forcing UNALIGNED there would permanently give up the aligned SIMD codelets that are the entire point of planning, purely to satisfy an over-specified test. The assertions become isapprox at rtol=1e-12, which is ~4 orders above roundoff and ~10 orders below what this test exists to catch. Verified that a simulated normalization revert still fails it: the relative error would be 0.845 for :schmidt and 4.29 for :fourpi. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Codecov Report❌ Patch coverage is
🚀 New features to boost your workflow:
|
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
1 participant
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.
Two rounds of review: first over this branch's diff, then over the ~9k lines of
src//ext/the branch never touched. Every finding was reproduced before fixing and re-verified after.Breaking changes
Both are documented with porting instructions in the new
CHANGELOG.md.analysis_axisym/analysis_axisym_lreturn different values. They omitted the φ quadrature factorcphi*nlon = 2π, so every coefficient was1/2πtoo small — they inverted neithersynthesis_axisymnor the m=0 column ofanalysis. Anything compensating for the old scale downstream must drop that compensation.dist_synthesis!/dist_synthesis_sphtor!still rejectreal_output=falsewith real-backed output arrays. I built the softer option (accept when the imaginary part is negligible) and measured that it cannot work:real_output=falseno longer means "the real field, typed complex" — it sums only the m≥0 half without the Hermitian mirror. On a typical config the result has|imag|up to 1.34 and a real part differing from the real field by 1.31, against field magnitude 2.96, so no tolerance bridges the two. Callers who wanted the real field should passreal_output=true; the error messages now say exactly that.Correctness
SHTPlanis FFTW-backed and cannot holdDual, so batch entry points now route non-FFTW eltypes through the plan-freecfg-form transforms. Batch gradients match the non-batch reference to 0.0.load_configsilently rebuilt:driscoll_healygrids as:regular, changing θ nodes and weights (8.7e-16 → 3e-3).save_confignow recordsuse_dh_weights.shtns_set_gridreturned success for unrecognized grid codes while producing a mirrored grid — the fallback branch was missing thereverse!.nlat == 1built an all-NaN config fromπ/0and returned NaN from every later transform without raising.dist_SH_mul_mx!crashed on everymres > 1config (missingm % mresguard).dist_SH_Yrotatecrashed onmres > 1. Worth noting the obvious fix is wrong: a Y-rotation mixes orders and is not representable in an mres-strided layout at all, so adding the stride just moves the failure deeper. It now rejects up front with the reason.device_transfer_arraysrejected:cuda/:amdgpu— the config's own vocabulary.M[l,m].Performance
φ_is_local_allandθ_is_distributedare reduced together in a singleAllreduceand cached per(pencil, communicator). The cache is identity-keyed with wholesale eviction so the hit/miss pattern stays identical on every rank — a rank-asymmetric cache would deadlock the collective.dist_synthesis_packed_cplxis single-pass, down from two full distributed syntheses. The negative-m φ bins are filled in the same θ/m traversal as the positive ones, reusing one Legendre row per(m, θ). It now matches the serial reference exactly.prepare_plm_tables!was buildingNP_tables/NdP_tablesbit-for-bit identical toplm_tables/dplm_tables; they now alias.estimate_table_memoryreported half the true figure, so jobs sized by it allocated twice their budget and were OOM-killed.UNALIGNED, removing a silent fallback to an O(n²) DFT; batch FFT helpers reuse the cache instead of re-planning per call.Internal and tests
pack_lm!/pack_lm/unpack_lm!/unpack_lminsrc/layout.jlreplace six open-coded copies of the packed↔dense(l,m)mapping. Them % mresguard had needed fixing three separate times across those copies.1/2πoff — and passed again after the fix. It now asserts the absolute identity and pins axisym analysis to the m=0 column ofanalysis.Verification
🤖 Generated with Claude Code