Feat: SVector as common currency

Project.toml

@@ -1,7 +1,7 @@
 name = "Gradus"
 uuid = "c5b7b928-ea15-451c-ad6f-a331a0f3e4b7"
 authors = ["fjebaker <[email protected]>"]
-version = "0.2.3"
+version = "0.2.4"
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"

src/Gradus.jl

@@ -81,6 +81,9 @@ export AbstractMetricParams,
     AbstractIntegrationParameters,
     IntegrationParameters
 
+# export static arrays too
+export SVector, @SVector
+
 """
     abstract type AbstractPointFunction
 

src/tracing/constraints.jl

@@ -1,22 +1,12 @@
 # vectors of vectors 
-function constrain_all(m::AbstractMetricParams, us, vs, μ)
+function constrain_all(m::AbstractMetricParams, us::AbstractArray{<:StaticVector}, vs::AbstractArray{<:StaticVector}, μ)
     @inbounds for i in eachindex(vs)
         vs[i] = constrain_all(m, us[i], vs[i], μ)
     end
     vs
 end
 
-function constrain_all(
-    m::AbstractMetricParams{T},
-    u::AbstractVector{T},
-    v::AbstractVector{T},
-    μ,
-) where {T<:Number}
-    v[1] = constrain(m, u, v, μ = μ)
-    v
-end
-
-function constrain_all(
+@inline function constrain_all(
     m::AbstractMetricParams{T},
     u::StaticVector{S,T},
     v::StaticVector{S,T},
@@ -25,7 +15,6 @@ function constrain_all(
     # mut = MVector{S,T}(v)
     # mut[1] = constrain(m, u, v, μ = μ)
     # SVector{S,T}(mut)
-    #
     # this results in no performance gain over the above, but i think it
     # reads nicer
     @inbounds SVector{S,T}(constrain(m, u, v, μ = μ), v[2], v[3], v[4])

src/tracing/method-implementations/auto-diff.jl

@@ -74,7 +74,7 @@ inv(metric)
         g4 = g_comp[4],
         g5 = g_comp[5]
 
-        (
+        SVector{5}(
             (g2 * g3 * g4) / (g1 * g2 * g3 * g4 - (g5^2) * g2 * g3),
             (g1 * g3 * g4 - (g5^2) * g3) / (g1 * g2 * g3 * g4 - (g5^2) * g2 * g3),
             (g1 * g2 * g4 - (g5^2) * g2) / (g1 * g2 * g3 * g4 - (g5^2) * g2 * g3),
@@ -176,7 +176,7 @@ jacobian = (j0, j1_mat, j2_mat, j0)
         v3 = v[3],
         v4 = v[4]
 
-        (
+        SVector{4}(
             -2(0.5gi5 * j14 + 0.5gi1 * j15) * v2 * v4 -
             2(0.5gi1 * j11 + 0.5gi5 * j15) * v1 * v2 -
             2(0.5gi1 * j21 + 0.5gi5 * j25) * v1 * v3 -
@@ -229,7 +229,7 @@ Limitations:
     end
 end
 
-function constrain(
+@inline function constrain(
     m::AbstractAutoDiffStaticAxisSymmetricParams{T},
     u,
     v;

src/tracing/tracing.jl

@@ -81,26 +81,6 @@ function integrator_problem(
     ODEProblem{false}(f, u_init, time_domain, IntegrationParameters(StatusCodes.NoStatus))
 end
 
-function integrator_problem(
-    m::AbstractMetricParams{T},
-    pos::AbstractVector{T},
-    vel::AbstractVector{T},
-    time_domain,
-) where {T}
-    u_init = vcat(pos, vel)
-
-    function f!(du, u::AbstractVector{T}, p, λ) where {T}
-        @inbounds let x = @view(u[1:4]), v = @view(u[5:8])
-            dv = SVector{4,T}(geodesic_eq(m, x, v))
-
-            du[1:4] .= v
-            du[5:8] .= dv
-        end
-    end
-
-    ODEProblem{true}(f!, u_init, time_domain, IntegrationParameters(StatusCodes.NoStatus))
-end
-
 # single position and single velocity
 function __tracegeodesics(
     m::AbstractMetricParams{T},

src/tracing/utility.jl

@@ -3,15 +3,6 @@
 
 Map impact parameters `α` and `β` to a velocity vector at some position `u` in the given metric `m`.
 """
-function map_impact_parameters(
-    m::AbstractMetricParams{T},
-    u::AbstractVector{T},
-    α::N,
-    β::N,
-) where {T,N<:Number}
-    [T(0.0), impact_parameters_to_vel(m, u, α, β)...]
-end
-
 function map_impact_parameters(
     m::AbstractMetricParams{T},
     u::SVector{S,T},
@@ -31,9 +22,9 @@ function map_impact_parameters(
     curried.(α, β)
 end
 
-function impact_parameters_to_vel(m::AbstractMetricParams{T}, u, α, β) where {T}
+@inline function impact_parameters_to_vel(m::AbstractMetricParams{T}, u, α, β) where {T}
     mcomp = metric_components(m, @view(u[2:3]))
     T(-1.0), -β / mcomp[3], -α / √(mcomp[3] * mcomp[4])
 end
 
-export map_impact_parameters, impact_parameters_to_vel
+export map_impact_parameters

test/smoke-tests/tracegeodesics.jl

@@ -38,22 +38,6 @@ using StaticArrays
         end
     end
 
-    @testset "regular-arrays" begin
-        u = Float64[0.0, 100.0, deg2rad(85), 0.0]
-        # arbitrary single velocity vector
-        v = Float64[0.0, -1.0, -0.0, -3.5e-6]
-        for m in [
-            KerrMetric(),
-            JohannsenMetric(),
-            KerrSpacetimeFirstOrder(),
-            MorrisThorneWormhole(),
-        ]
-            test_single(m, u, v)
-            test_many(m, u, v)
-            @test true
-        end
-    end
-
     @testset "corona-models" begin
         # only implemented for KerrMetric at the moment
         # because of the vector to local sky methods