Unify normalization on orthonormal; fix 20+ correctness bugs found by audit - #29
Merged
Merged
Conversation
… audit
Settles a long-standing split where the distributed layer, the QST family and
SHTPlan's scalar path exchanged coefficients in the config's normalization while
serial analysis/synthesis and the Parseval energy diagnostics assumed orthonormal.
Everything is orthonormal now, which makes three previously-untestable properties
exact invariants:
dist_analysis(cfg,f) == analysis(cfg,f) 2.5e-16, every norm
QST: Q == analysis(Vr), Vr == synthesis(Q) 0.0e+00
analysis!(SHTPlan, ...) == analysis 0.0e+00
QST needs care: its S/T come from the still-converting serial sphtor pair, so
analysis_qst converts them back to orthonormal rather than dropping a conversion
(which would re-create the mixed-convention bug). The packed API, serial sphtor,
GPU and LoopVec keep converting -- those are cfg-convention by contract.
Correctness fixes, all verified:
Transforms
- dist_synthesis / dist_synthesis_sphtor / 2d_optimized: real_output=false built
a Float64 destination, selecting the real-output ifft method -- the imaginary
half was discarded and the real part halved.
- dist_synthesis_packed_cplx dropped every m<0 coefficient.
- synthesis_sphtor_ml omitted phi_inv_scale, so the _ml pair was not an inverse.
- analysis/synthesis_sphtor_ml took a parameter named `im`, shadowing Julia's
imaginary unit; `1.0im` is juxtaposition, so the S/T coupling term evaluated
REAL and the two components cross-contaminated.
- shtns_rotation_apply_cplx rotated in the Y_l^m basis while every other LM_cplx
function uses P_l^{|m|}; rotating a real field returned a complex one. Fixed
via an explicit basis change; the real path is bit-identical.
- Batch: in-place ifft fallback restarted at k=1 over already-transformed slices,
double-inverting field 1 (~102% error at odd nlon); irfft threw on ComplexF32.
Plans now use FFTW.UNALIGNED and the fallback resumes past completed slices.
- Batch output eltype followed Float64 regardless of input.
MPI
- theta- and phi-locality predicates were evaluated per-rank while gating
collectives, hanging when a pencil has more partitions than points. Reduced at
all 13 sites; the four that only slice a local window are left per-rank.
- _phi_column_color double-counted an empty-phi rank's slab after a gather; the
helper is removed and its reasoning folded into _keep_one_phi_partner!.
- dist_analysis_packed / dist_synthesis_packed / dist_SH_mul_mx! threw for
mres>1; dist_SH_mul_mx! also read unowned l from uninitialized memory and left
skipped columns unwritten.
Automatic differentiation
- The packed and packed_cplx rrules used each transform's inverse as its adjoint;
the old synthesis_packed adjoint was off by per-(l,m) factors of 2.7-43x.
- _adjoint_synthesis assumed phi_inv_scale == nlon, making every synthesis-family
gradient exactly 2*pi too large under phi_scale=:quad.
- Distributed rrules recomputed first(phi_globals) eagerly, throwing on a
zero-column rank inside a collective region.
- The rotation pullback assigned to a captured variable, so angle gradients were
wrong from the second invocation onward.
- Zero cotangents crashed the sphtor pullbacks; rotation angle gradients read
uninitialized memory (real) or threw (complex) when mmax < lmax.
GPU
- sphtor kernels lacked the pole-limit substitution the CPU kernels have, so the
m=1 contribution vanished from pole rows. NOT runtime-verified: no CUDA device
was available.
Verification: serial suite 67268/67268; 6 MPI suites at 4 ranks; 56 custom MPI
assertions at 8 ranks including empty-partition topologies; finite-difference
checks across norm x phi_scale for the packed, sphtor and rotation gradients.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The orthonormal refactor missed ext/ParallelRotationsPencil.jl, which still converted cfg->internal on input and back on output in both Y-rotation paths (dist_SH_Yrotate_allgatherm! and _yrotate_truncgather_rows!). That does not cancel: the Wigner rotation mixes m, so R(M*a)/M != R(a). With the distributed spectral arrays now orthonormal, the conversion corrupted every non-default-norm rotation. Measured against serial SH_Yrotate, which is itself norm-agnostic: exact at :orthonormal but 0.68 off at :schmidt and 2.86 at :fourpi/cs_phase=false. All three now match serial exactly (0.0e+00). Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
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.
Eight review rounds over this codebase. Two themes: a normalization convention that was split across backends, and a cluster of correctness bugs that the split (and a few silent-failure patterns) had been hiding.
The convention change
The distributed layer, the QST family and
SHTPlan's scalar path exchanged coefficients in the config's normalization; serialanalysis/synthesisand the Parseval energy diagnostics assumed orthonormal. Samealm, two readings.Everything is orthonormal now. That turns three previously-untestable properties into exact invariants:
dist_analysis(cfg,f) == analysis(cfg,f), every normQ == analysis(Vr),Vr == synthesis(Q)analysis!(SHTPlan, …) == analysisQST needed care: its S/T come from the still-converting serial sphtor pair, so
analysis_qstconverts them back to orthonormal. Simply dropping a conversion there would have re-created the mixed-convention bug. The packed API, serial sphtor, GPU and LoopVec keep converting — those are cfg-convention by contract.dist_analysisconverted before this PR too, so downstream code assuming cfg-form distributed coefficients will now see orthonormal ones. Worth a release note.Correctness fixes
Transforms —
real_output=falsebuilt aFloat64destination, so the imaginary half was discarded and the real part halved (three functions).dist_synthesis_packed_cplxdropped everym<0coefficient.synthesis_sphtor_mlomittedphi_inv_scale.shtns_rotation_apply_cplxrotated in the Y_l^m basis while every other LM_cplx function uses P̄_l^{|m|} — rotating a real field returned a complex one.A parameter named
im.analysis/synthesis_sphtor_mltookim::Int, shadowing Julia's imaginary unit.1.0imis juxtaposition (1.0 * im), not a lexer suffix, so the S/T coupling term evaluated real and the components cross-contaminated. Renamed tomvalso it cannot recur.Batch FFT — the in-place ifft fallback restarted at
k=1over already-transformed slices, double-inverting field 1 (~102% error at oddnlon, i.e. this repo's own2·lmax+1convention); evennlonhappened to align, so it was silent and size-dependent.MPI — θ- and φ-locality predicates were evaluated per-rank while gating collectives, hanging when a pencil has more partitions than points. Reduced at all 13 sites.
_phi_column_colordouble-counted an empty-φ rank's slab after a gather and is removed. Several functions threw formres>1.AD — the packed rrules used each transform's inverse as its adjoint (old
synthesis_packedadjoint off by per-(l,m) factors of 2.7–43×)._adjoint_synthesisassumedphi_inv_scale == nlon, making every synthesis-family gradient exactly 2π too large underphi_scale=:quad. The rotation pullback assigned to a captured variable, so angle gradients were wrong from the second call onward.Verification
MPITopologynorm × phi_scalefor the packed, sphtor and rotation gradientsNot verified: the GPU pole-limit fix. No CUDA device was available — only the formulas (identical to the FD-verified CPU originals) and compilation were checked. Please exercise that path before relying on it.
Reviewer's note: many of these were found by reviewing the fixes themselves rather than the original code — several rounds caught regressions introduced by the round before. The diff is worth reading with that in mind.
🤖 Generated with Claude Code