From c0d95bbe30020f34c73ffbde3d1d426d846f16ab Mon Sep 17 00:00:00 2001 From: mattsignorelli Date: Thu, 18 Dec 2025 12:14:44 -0500 Subject: [PATCH 1/2] NNF inverter --- src/map.jl | 57 +++++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 56 insertions(+), 1 deletion(-) diff --git a/src/map.jl b/src/map.jl index b3871b3..55d6441 100644 --- a/src/map.jl +++ b/src/map.jl @@ -547,13 +547,68 @@ function inv!( m::TaylorMap, m1::TaylorMap; do_spin::Bool=true, + work_vs::NTuple{3}=m.v isa StaticArray ? + ntuple(t->StaticArrays.sacollect(MVector{nvars(m)}, zero(first(m.v)) for i in 1:nvars(m)), Val{3}()) : + ntuple(t->zero.(view(m.v, 1:nvars(m))), Val{3}()), work_q::Union{Nothing,Quaternion}=isnothing(m.q) ? nothing : zero(m.q) ) checkinplace(m, m1) !(m === m1) || error("Cannot inv!(m, m1) with m === m1") - TI.inv!(m.v, m1.v) + #TI.inv!(m.v, m1.v) + + # Following Dragt book: + # We can use the equation on Wikipedia + # for inverse of block matrices: https://en.wikipedia.org/wiki/Block_matrix#Inversion + A = jacobian(m1, VARS) + B = jacobian(m1, PARAMS) + Ri = hcat(inv(A), -inv(A)*B) + + clear!(m) + setray!(m.v, v_matrix=Ri) + + mo = maxord(m1) + if mo > 1 + nn = ndiffs(m) + nv = nvars(m) + np = nparams(m) + + ri = work_vs[1] + Ntilde = work_vs[2] + tmp = work_vs[3] + + # Put N-not-tilde in tmp + for j in 1:nv + TI.clear!(ri[j]) + TI.clear!(Ntilde[j]) + TI.cutord!(tmp[j], m1.v[j], -1) + end + + # (nv x nn) * (nn x 1) = (nv x 1) + # Use first TPSA in ri as temporary + tt = first(ri) + for j in 1:nv + for i in 1:nv + TI.mul!(tt, -Ri[i,j], tmp[j]) + TI.add!(Ntilde[i], Ntilde[i], tt) + end + end + + # now can set ri, tt dead + TI.clear!(tt) + setray!(ri, v_matrix=Ri) + + # At each step in iteration, m should contain the result up to + # that order + for i in 2:mo + TI.compose!(view(tmp, 1:nv), view(Ntilde, 1:nv), m.v) + for j in 1:nv + TI.cutord!(tmp[j], tmp[j], mo+1) + TI.add!(m.v[j], ri[j], tmp[j]) + end + end + end # Now do quaternion: inverse of q(z0) is q^-1(M^-1(zf)) if !isnothing(m.q) From 66ef809b64b4018a334651bfc423cb1cbedfa42b Mon Sep 17 00:00:00 2001 From: mattsignorelli Date: Fri, 19 Dec 2025 08:21:28 -0500 Subject: [PATCH 2/2] own inv, fix sagan-rubin --- src/map.jl | 2 -- src/sagan_rubin.jl | 40 ++++++++++++++++++++++++++++------------ 2 files changed, 28 insertions(+), 14 deletions(-) diff --git a/src/map.jl b/src/map.jl index 55d6441..22a2621 100644 --- a/src/map.jl +++ b/src/map.jl @@ -556,8 +556,6 @@ function inv!( checkinplace(m, m1) !(m === m1) || error("Cannot inv!(m, m1) with m === m1") - #TI.inv!(m.v, m1.v) - # Following Dragt book: # We can use the equation on Wikipedia # for inverse of block matrices: https://en.wikipedia.org/wiki/Block_matrix#Inversion diff --git a/src/sagan_rubin.jl b/src/sagan_rubin.jl index 465a58a..110cfe6 100644 --- a/src/sagan_rubin.jl +++ b/src/sagan_rubin.jl @@ -12,6 +12,22 @@ function compute_sagan_rubin(a1::DAMap{V0}, ::Val{linear}=Val{false}()) where {V # If we are coasting, then we will not return any dispersion stuff, since # that is all well-defined in the parameter-dependent fixed point. + + # Note: generators (for i in _ for j in _) act in ORDER. So this would be like + # for i in _ + # for j in _ + # ... + + # Then, if a matrix is constructed from this result, it is filled in column-major order + #= Observe: + julia> StaticArrays.sacollect(SMatrix{2,2}, (i, j) for j=1:2 for i=1:2) + 2×2 SMatrix{2, 2, Tuple{Int64, Int64}, 4} with indices SOneTo(2)×SOneTo(2): + (1, 1) (1, 2) + (2, 1) (2, 2) + + The tuple of indices works when COLUMN is first + So the order is basically for col in cols for row in rows + =# nv = nvars(a1) nhv = nhvars(a1) coast = iscoasting(a1) @@ -20,12 +36,12 @@ function compute_sagan_rubin(a1::DAMap{V0}, ::Val{linear}=Val{false}()) where {V let a1_mat = jacobian(a1, HVARS), a1i_mat = inv(a1_mat), nvm = nv-2 H = StaticArrays.sacollect(SVector{Int(nv/2),typeof(a1_mat)}, a1_mat*ip_mat(a1, i)*a1i_mat for i in 1:Int(nv/2)) gamma_c = sqrt(H[1][1,1]) - C = nv > 2 ? StaticArrays.sacollect(SMatrix{2,2,eltype(H[1])}, H[1][i,j] for i in 3:4 for j in 3:4) : 0 + C = nv > 2 ? StaticArrays.sacollect(SMatrix{2,2,eltype(H[1])}, -H[1][row,col]/gamma_c for col in 3:4 for row in 1:2) : 0 Ct = nv > 2 ? SA[C[2,2] -C[1,2]; -C[2,1] C[1,1]] : 0 Vi = gamma_c*I + vcat(hcat(zero(C), -C), hcat(Ct, zero(Ct))) - N = Vi*StaticArrays.sacollect(SMatrix{nvm,nvm,eltype(a1_mat)}, a1_mat[i,j] for i in 1:nvm for j in 1:nvm) + N = Vi*StaticArrays.sacollect(SMatrix{nvm,nvm,eltype(a1_mat)}, a1_mat[row,col] for col in 1:nvm for row in 1:nvm) beta = StaticArrays.sacollect(SVector{Int(nvm/2),eltype(N)}, N[i,i]^2 for i in 1:2:Int(nvm)) - alpha = StaticArrays.sacollect(SVector{Int(nvm/2),eltype(N)}, -N[i,i+1]*N[i,i] for i in 1:2:Int(nvm)) + alpha = StaticArrays.sacollect(SVector{Int(nvm/2),eltype(N)}, -N[i+1,i]*N[i,i] for i in 1:2:Int(nvm)) # Note that eta and zeta here are in x,px,y,py, not in a,pa,b,pb eta = StaticArrays.sacollect(SVector{nvm,eltype(H[1])}, H[end][i,nv] for i in 1:(nvm)) zeta = StaticArrays.sacollect(SVector{nvm,eltype(H[1])}, H[end][i,nv-1] for i in 1:(nvm)) @@ -40,12 +56,12 @@ function compute_sagan_rubin(a1::DAMap{V0}, ::Val{linear}=Val{false}()) where {V let a1_mat = jacobian(a1, HVARS), a1i_mat = inv(a1_mat) H = StaticArrays.sacollect(SVector{Int(nhv/2),typeof(a1_mat)}, a1_mat*ip_mat(a1, i)*a1i_mat for i in 1:Int(nhv/2)) gamma_c = sqrt(H[1][1,1]) - C = nv > 2 ? StaticArrays.sacollect(SMatrix{2,2,eltype(H[1])}, H[1][i,j] for i in 3:4 for j in 3:4) : 0 + C = nv > 2 ? StaticArrays.sacollect(SMatrix{2,2,eltype(H[1])}, -H[1][row,col]/gamma_c for col in 3:4 for row in 1:2) : 0 Ct = nv > 2 ? SA[C[2,2] -C[1,2]; -C[2,1] C[1,1]] : 0 Vi = gamma_c*I + vcat(hcat(zero(C), -C), hcat(Ct, zero(Ct))) - N = Vi*StaticArrays.sacollect(SMatrix{nhv,nhv,eltype(a1_mat)}, a1_mat[i,j] for i in 1:nhv for j in 1:nhv) + N = Vi*StaticArrays.sacollect(SMatrix{nhv,nhv,eltype(a1_mat)}, a1_mat[row,col] for col in 1:nhv for row in 1:nhv) beta = StaticArrays.sacollect(SVector{Int(nhv/2),eltype(N)}, N[i,i]^2 for i in 1:2:nhv) - alpha = StaticArrays.sacollect(SVector{Int(nhv/2),eltype(N)}, -N[i,i+1]*N[i,i] for i in 1:2:nhv) + alpha = StaticArrays.sacollect(SVector{Int(nhv/2),eltype(N)}, -N[i+1,i]*N[i,i] for i in 1:2:nhv) return (; beta=beta, alpha=alpha, gamma_c=gamma_c, C=C) end end @@ -61,12 +77,12 @@ function compute_sagan_rubin(a1::DAMap{V0}, ::Val{linear}=Val{false}()) where {V ) mo = maxord(a1) gamma_c = TI.cutord(sqrt(factor_out(H[1].v[1], 1)), mo) - C = nv > 2 ? StaticArrays.sacollect(SMatrix{2,2,eltype(H[1].v)}, factor_out(H[1].v[i], j) for i in 3:4 for j in 3:4) : 0 + C = nv > 2 ? StaticArrays.sacollect(SMatrix{2,2,eltype(H[1].v)}, TI.cutord(-factor_out(H[1].v[row], col)/gamma_c, mo) for col in 3:4 for row in 1:2) : 0 Ct = nv > 2 ? SA[C[2,2] -C[1,2]; -C[2,1] C[1,1]] : 0 Vi = gamma_c*I + vcat(hcat(zero(C), -C), hcat(Ct, zero(Ct))) - N = TI.cutord.(Vi*StaticArrays.sacollect(SMatrix{nvm,nvm,eltype(a1.v)}, factor_out(a1.v[i], j) for i in 1:nvm for j in 1:nvm), mo) + N = TI.cutord.(Vi*StaticArrays.sacollect(SMatrix{nvm,nvm,eltype(a1.v)}, factor_out(a1.v[row], col) for col in 1:nvm for row in 1:nvm), mo) beta = StaticArrays.sacollect(SVector{Int(nvm/2),eltype(N)}, TI.cutord(N[i,i]^2, mo) for i in 1:2:Int(nvm)) - alpha = StaticArrays.sacollect(SVector{Int(nvm/2),eltype(N)}, TI.cutord(-N[i,i+1]*N[i,i], mo) for i in 1:2:Int(nvm)) + alpha = StaticArrays.sacollect(SVector{Int(nvm/2),eltype(N)}, TI.cutord(-N[i+1,i]*N[i,i], mo) for i in 1:2:Int(nvm)) # Note that eta and zeta here are in x,px,y,py, not in a,pa,b,pb eta = StaticArrays.sacollect(SVector{nvm,eltype(H[1].v)}, factor_out(H[end].v[i], nv) for i in 1:(nvm)) zeta = StaticArrays.sacollect(SVector{nvm,eltype(H[1].v)}, factor_out(H[end].v[i], nv-1) for i in 1:(nvm)) @@ -93,12 +109,12 @@ function compute_sagan_rubin(a1::DAMap{V0}, ::Val{linear}=Val{false}()) where {V ) mo = maxord(a1) gamma_c = TI.cutord(sqrt(factor_out(H[1].v[1], 1)), mo) - C = nv > 2 ? StaticArrays.sacollect(SMatrix{2,2,eltype(H[1].v)}, factor_out(H[1].v[i], j) for i in 3:4 for j in 3:4) : 0 + C = nv > 2 ? StaticArrays.sacollect(SMatrix{2,2,eltype(H[1].v)}, TI.cutord(-factor_out(H[1].v[row], col)/gamma_c, mo) for col in 3:4 for row in 1:2) : 0 Ct = nv > 2 ? SA[C[2,2] -C[1,2]; -C[2,1] C[1,1]] : 0 Vi = gamma_c*I + vcat(hcat(zero(C), -C), hcat(Ct, zero(Ct))) - N = TI.cutord.(Vi*StaticArrays.sacollect(SMatrix{nhv,nhv,eltype(a1.v)}, factor_out(a1.v[i], j) for i in 1:nhv for j in 1:nhv), mo) + N = TI.cutord.(Vi*StaticArrays.sacollect(SMatrix{nhv,nhv,eltype(a1.v)}, factor_out(a1.v[row], col) for col in 1:nhv for row in 1:nhv), mo) beta = StaticArrays.sacollect(SVector{Int(nhv/2),eltype(N)}, TI.cutord(N[i,i]^2, mo) for i in 1:2:Int(nhv)) - alpha = StaticArrays.sacollect(SVector{Int(nhv/2),eltype(N)}, TI.cutord(-N[i,i+1]*N[i,i], mo) for i in 1:2:Int(nhv)) + alpha = StaticArrays.sacollect(SVector{Int(nhv/2),eltype(N)}, TI.cutord(-N[i+1,i]*N[i,i], mo) for i in 1:2:Int(nhv)) return (; beta=beta, alpha=alpha, gamma_c=gamma_c, C=C) end end