GPU-safe constructors, custom metric signatures, and Julia 1.10 LTS - #37
GPU-safe constructors, custom metric signatures, and Julia 1.10 LTS#37AndersEdin wants to merge 5 commits into
Conversation
Make arrays of Clifford numbers usable on the GPU.
- Constructors: replace interpolated `@assert` messages in
`check_element_count`, the `KVector` K-range check, and the
`Z2CliffordNumber` parity check with constant `String` literals. An
interpolated message allocates a `String` in the throw path (via
`print_to_string`/`apply_generic`), which GPU compilers reject
(`jl_alloc_string` is unsupported on-device). The checks are kept; for
concrete types they are statically true and fully eliminated, so the
constructors emit zero throws and zero allocations.
- Adapt: add `ext/CliffordNumbersAdaptExt.jl` registering
`Adapt.adapt_structure`. These types are `isbits` with `NTuple{L,T}`
storage, so the only meaningful adaptation is retyping the scalar
coefficients to match a typed-array destination (e.g.
`adapt(Array{Float32}, ::KVector{...,Float64})` -> `...,Float32`); other
destinations preserve the scalar type. Wired into `[weakdeps]`/
`[extensions]`/`[compat]` and covered by `test/ext/Adapt.jl`.
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
On-device testing (RTX 3070, CUDA.jl) showed that `reverse`, `adjoint`, `conj`, and `grade_involution` did not lower on the GPU for `CliffordNumber` (dense) and `Z2CliffordNumber`, and `reverse` failed for `KVector` too: the generic fallback `T(x[f.(BladeIndices(T))])` resolves blade positions at runtime through `to_index`, which the GPU compiler rejects (gpu_gc_pool_alloc + dynamic dispatch). This blocked the headline `rotate.(cu_pts, Ref(R))` / rotor-sandwich path, since the sandwich uses `R'`. - Replace the generic indexing loop with a `@generated` `_sign_automorphism` that bakes every blade position and sign at compile time, emitting pure tuple arithmetic (allocation-free, GPU-safe). Works for all dense / Z2 / user types. - Add the closed-form `reverse(::KVector)` (identical to `adjoint` for real scalars), the one case the existing KVector closed forms had missed. Verified: the four automorphisms and a rotor-sandwich broadcast now compile and run on-device and match the CPU result. A new CPU regression testset pins all four automorphisms against the original indexing semantics across dense / Even / Odd / KVector over VGA(3), PGA(3), and STA. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
metric_tuple(s) returns a signature's metric as an NTuple of Int8 and is the single, world-age-robust source of metric data: the generated blade kernels and the negative/zero square-bit masks both read it, and it works on a downstream AbstractSignature subtype by reading its interface at the caller's world age rather than inside a generator. This makes user-defined signatures a supported extension path. A custom metric is expressed either as a Signature *value* placed in the Clifford number's type parameter, or as an AbstractSignature subtype bridged to a value with Signature(s) (mirroring how CliffordNumbers' own STAP is built). docs/src/custom_signatures.md walks through both, and test/custom_signature.jl exercises a phase-space (D, D, R) family: the signature interface, metric_tuple as the single metric source, the world-age limitation of a bare subtype, and the Signature(s) bridge driving the generated product/automorphism/inverse kernels end to end. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Raise the compat lower bound to 1.10 in the package and test environments (1.9 support was already nominal: package extensions are the only 1.9 feature in use). dimension of a plain Signature now returns a full-width Int instead of an Int8: on 1.10, binomial with a narrow integer falls back to the Float64 gamma-function method, which made nblades of a KVector in a custom-signature algebra return a Float64 and broke construction. The zero-allocation test gates assert what the current stable compiler constant-folds; the 1.10 compiler misses some of those folds, so ALLOCATION_GATES restricts them to 1.12 and later. Correctness and inference tests are unchanged and run on every version. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
brainandforce
left a comment
There was a problem hiding this comment.
Some overall comments here:
- There are a few comments here about the changes to error messages in the constructors. There are some deeper underlying issues that need to be rectified (long story short, I misunderstood the use of
@assert: it can be disabled at some optimization levels and it shouldn't be used when verifying correctness of inputs.) The full fix might have to wait on a future discussion of the package architecture (specifically because we may need more customExceptionsubtypes, but in the meantime, you can turn those intoerrorstatements. Generally we want to avoid throwingErrorExceptionbecause it lacks detail, but it's better than@assert. - For error messages in constructors, we should look into whether
Base.LazyStringwould allow us to retain the same error messages as before, but if that's not possible the current ones are fine. - Also, per Julia convention, the first character of an error message should be lowercase. Feel free to fix this wherever you see it; it's on me for not following it consistently.
- There may need to be some style changes in the documentation, but that's not something you should worry about with this PR - I'll likely take care of it myself. The comments there are also for my future reference.
- On the other hand, please verify that the information provided in generated documentation is correct. There have been a couple of instances where the documentation may refer to incorrect implementation details, nonexistent packages, or names that don't seem to come from CliffordNumbers.jl or its dependencies (I have commented where I suspect there are problems).
- In the future, I'd prefer PRs to be even more granular than this: I'd rather get 100 PRs that change only a few lines and deal with the merge conflicts they'd cause rather than a few huge ones that change a lot of things. Hopefully it'll also be less overwhelming on your end to get reviews on smaller PRs as well. (For reference, I'd have made separate PRs for the Adapt extension, the fixes to constructors to prevent branching, the changes to how signatures will work to avoid world age issues, the documentation page on defining your own signatures, and the bump to 1.10 LTS as the minimum supported version if I was making these changes.)
| # For concrete types these checks are statically true and get eliminated entirely on device. | ||
| @inline function check_element_count(sz, data) | ||
| @assert length(data) == sz "Expected $sz scalars from input, got $(length(data))" | ||
| @assert length(data) == sz "Number of input scalars does not match the expected count for this type." |
There was a problem hiding this comment.
Two notes here:
- The way I've handled errors here is problematic:
@assertshould not be used here, because@assertis not guaranteed to run at certain optimization levels. The entire way errors are handled across the constructors needs to be fixed. - Can we use
Base.LazyStringhere? This is specifically designed to avoid allocations.
| # The closed-form (KVector) and `@generated` (dense / Z2) implementations must agree with the | ||
| # original generic indexing semantics `T(x[f.(BladeIndices(T))])` for every type and grade. | ||
| # This also locks in `reverse(::KVector)` and the GPU-safe generated paths in src/math/duals.jl. | ||
| refauto(f, x) = (T = typeof(x); T(x[f.(CliffordNumbers.BladeIndices(T))])) |
There was a problem hiding this comment.
What does refauto mean here? This may need a more descriptive name.
| Base.has_offset_axes(::AbstractSignature) = true | ||
|
|
||
| """ | ||
| dimension(s::AbstractSignature) -> Int8 |
There was a problem hiding this comment.
The reason this is an Int8 is because there is no practical reason to use a wider type. Due to how Julia NTuple is implemented, the library works best with algebras that have 6 or fewer dimensions (though I hope this changes in the future...), and on a 64-bit computer it wouldn't be possible to work with algebras with more than 64 dimensions anyway.
| is equal to `2^dimension(s)`. | ||
| """ | ||
| dimension(s::AbstractSignature) = signed(s.dimensions) | ||
| dimension(s::AbstractSignature) = Int(s.dimensions) |
There was a problem hiding this comment.
Either way, the correct solution here is to convert the arguments to the binomial call in nblades to Int64 if needed, rather than the output type of dimension.
| `Int(s.dimensions)`. | ||
|
|
||
| The result is a full-width `Int` so that counts derived from it stay integers: on Julia 1.10, | ||
| `binomial` with a narrower integer argument falls back to the `Float64` gamma-function method, |
There was a problem hiding this comment.
I did some digging in the Julia 1.10 code, and I could not find any reference to an implementation of binomial that calls the gamma function, because Julia Base does not implement this function at all (it is only present in SpecialFunctions.jl).
| ## Generated kernels and world age | ||
|
|
||
| CliffordNumbers.jl compiles its hot paths — the geometric product, the automorphisms, the scalar | ||
| product, and `charpoly_coeffs` — with |
| A dedicated type is useful for a family of related algebras or to attach domain meaning to a | ||
| signature. Subtype [`Metrics.AbstractSignature`](@ref CliffordNumbers.Metrics.AbstractSignature) and | ||
| implement the interface: `dimension`, `firstindex`, and `getindex` are required, while | ||
| `is_degenerate` and `is_positive_definite` are recommended. |
There was a problem hiding this comment.
Also, one of the reasons why you'd want to derive a custom type instead of using Signature is because you can provide custom pretty printing behavior for the basis blades. That's how APS uses sigma and LGA uses gamma in its printing.
This would be useful information to include but is not critical to this PR.
| @assert 0 <= K <= dimension(Q) "K can only range from 0 to $(dimension(Q)) (got $K)." | ||
| # Constant (non-interpolated) assertion message to keep the throw path GPU-safe; an | ||
| # interpolated message allocates a `String`, which GPU compilers cannot lower. | ||
| @assert 0 <= K <= dimension(Q) "K is outside the valid grade range 0:dimension(Q)." |
There was a problem hiding this comment.
As mentioned elsewhere, this shouldn't be an @assert (not sure if ArgumentError, DomainError, or some other exception type would be the best choice).
If needed we can derive our own Exception subtypes (maybe GradeException, but we should discuss this and perhaps fall back on an error call in the meantime).
| value's fields. The package's own projective spacetime algebra | ||
| [`STAP`](@ref CliffordNumbers.Metrics.STAP) is defined this way. | ||
|
|
||
| Pass a fourth argument to set a non-default first index. Projective and Lorentzian algebras |
There was a problem hiding this comment.
probably note that it defaults to 1
|
|
||
| ### Fixed | ||
| - `dimension` of a plain `Signature` returns a full-width `Int` instead of an `Int8`. On Julia | ||
| 1.10 the narrow integer sent `binomial` to its `Float64` gamma-function method, so `nblades` |
There was a problem hiding this comment.
This needs to be reworded since it does not seem there has ever been a gamma function fallback in Julia 1.10.
Here comes the bug fixes and maintenance fixes from my previous PR. The issues were found when implementing a geometric algebra package using CliffordNumbers, and both the GPU fixes and the fix to enable defining metric signatures are needed for being able to use the package.
I have raised the minimum Julia to the 1.10 LTS, which makes sense since it is a lot of work to support older versions. This also aligns with what is already written in the package documentation.
Changes
Julia 1.10 LTS. Compat lower bound raised to
1.10in the package and test environments.dimension(::Signature)now returns a full-widthIntinstead ofInt8: on 1.10 a narrow integer sendsbinomialto itsFloat64gamma-function method, which madenbladesof aKVectorin a custom-signature algebra return aFloat64and broke construction. Zero-allocation test gates are placed behindALLOCATION_GATES = VERSION >= v"1.12"; the 1.10 compiler misses some of the constant-folds those gates assert, so they are informative-only there while correctness and inference tests run on every version.GPU enablement. Adds an
Adaptpackage extension soAbstractCliffordNumbervalues adapt to/from device arrays, plus constructor paths that allocate no boxed intermediates on the device path. The sign-changing automorphisms (reverse,grade_involution,conj) are reimplemented as branch-free bit arithmetic on the blade grade — reversion(-1)^{k(k-1)/2}, grade involution(-1)^k— so they carry no data-dependent branch and codegen uniformly on the GPU. Results are unchanged.Custom metric signatures.
metric_tuple(s) -> NTuple{N,Int8}is the single, world-age-robust source of metric data; the generated kernels and thenegative/zerosquare-bit masks all read it. A user metric is expressed either as aSignaturevalue in the type parameter, or as anAbstractSignaturesubtype reduced to a value withSignature(s)(the built-inSTAPuses the same shape). A bare downstream subtype cannot drive the@generatedkernels — their generators run at CliffordNumbers' load-time world age and cannot see interface methods defined later — soSignature(s), computed at the caller's world age, is the bridge. Tutorial indocs/src/custom_signatures.md.Tests
test/custom_signature.jlexercises a phase-space(D, D, R)family: signature interface,nbladesintegrality on 1.10,metric_tupleas the sole metric source, the world-age failure of a bare subtype, and theSignature(s)bridge driving the generated product/automorphism/inverse kernels. GPU extension covered bytest/ext/Adapt.jl.Pkg.test()passes on the 1.10 LTS.Breaking