Skip to content

Commit

Permalink
extended support for dyadics
Browse files Browse the repository at this point in the history
  • Loading branch information
@chakravala
chakravala committed Apr 26, 2024
1 parent b36c93b commit 90ba554
Show file tree
Hide file tree
Showing 4 changed files with 151 additions and 20 deletions.
15 changes: 7 additions & 8 deletions src/Grassmann.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,12 +26,12 @@ import AbstractTensors: Values, Variables, FixedVector, clifford, hodge, wedge,
export ⊕, ℝ, @V_str, @S_str, @D_str, Manifold, Submanifold, Signature, DiagonalForm, value
export @basis, @basis_str, @dualbasis, @dualbasis_str, @mixedbasis, @mixedbasis_str, Λ
export ℝ0, ℝ1, ℝ2, ℝ3, ℝ4, ℝ5, ℝ6, ℝ7, ℝ8, ℝ9, mdims, tangent, metric, antimetric
export hodge, wedge, vee, complement, dot, antidot
export hodge, wedge, vee, complement, dot, antidot, istangent

import Base: @pure, ==, isapprox
import Base: print, show, getindex, setindex!, promote_rule, convert, adjoint
import DirectSum: V0, ⊕, generate, basis, getalgebra, getbasis, dual, Zero, One, Zero, One
import Leibniz: hasinf, hasorigin, dyadmode, value, pre, vsn, metric, mdims
import Leibniz: hasinf, hasorigin, dyadmode, value, pre, vsn, metric, mdims, gdims
import Leibniz: bit2int, indexbits, indices, diffvars, diffmask
import Leibniz: symmetricmask, indexstring, indexsymbol, combo, digits_fast
import DirectSum: antimetric
Expand Down Expand Up @@ -117,8 +117,8 @@ end
d(ω::T) where T<:TensorAlgebra = Manifold(ω)(∇)∧ω
δ(ω::T) where T<:TensorAlgebra = -∂(ω)

function boundary_rank(t,d=gdims(t))
out = gdims(∂(t))
function boundary_rank(t,d=count_gdims(t))
out = count_gdims(∂(t))
out[1] = 0
for k ∈ 2:length(out)-1
@inbounds out[k] = min(out[k],d[k+1])
Expand All @@ -127,7 +127,7 @@ function boundary_rank(t,d=gdims(t))
end

function boundary_null(t)
d = gdims(t)
d = count_gdims(t)
r = boundary_rank(t,d)
l = length(d)
out = zeros(Variables{l,Int})
Expand All @@ -143,7 +143,7 @@ end
Compute the Betti numbers.
"""
function betti(t::T) where T<:TensorAlgebra
d = gdims(t)
d = count_gdims(t)
r = boundary_rank(t,d)
l = length(d)-1
out = zeros(Variables{l,Int})
Expand Down Expand Up @@ -569,8 +569,7 @@ function __init__()
DyadicChain(m::StaticArrays.SMatrix{N,N}) where N = Chain{Submanifold(N),1}(m)
Chain{V,G}(m::StaticArrays.SMatrix{N,N}) where {V,G,N} = Chain{V,G}(Chain{V,G}.(getindex.(Ref(m),:,StaticArrays.SVector{N}(1:N))))
Chain{V,G,Chain{W,G}}(m::StaticArrays.SMatrix{M,N}) where {V,W,G,M,N} = Chain{V,G}(Chain{W,G}.(getindex.(Ref(m),:,StaticArrays.SVector{N}(1:N))))
Base.exp(A::Chain{V,G,<:Chain{V,G}}) where {V,G} = Chain{V,G}(exp(StaticArrays.SMatrix(A)))
Base.log(A::Chain{V,G,<:Chain{V,G}}) where {V,G} = Chain{V,G}(log(StaticArrays.SMatrix(A)))
#Base.log(A::Chain{V,G,<:Chain{V,G}}) where {V,G} = Chain{V,G}(log(StaticArrays.SMatrix(A)))
LinearAlgebra.eigvals(A::Chain{V,G,<:Chain{V,G}}) where {V,G} = Chain(Values{binomial(mdims(V),G)}(LinearAlgebra.eigvals(StaticArrays.SMatrix(A))))
LinearAlgebra.eigvecs(A::Chain{V,G,<:Chain{V,G}}) where {V,G} = Chain(Chain.(Values{binomial(mdims(A),G)}.(getindex.(Ref(LinearAlgebra.eigvecs(StaticArrays.SMatrix(A))),:,list(1,binomial(mdims(A),G))))))
function LinearAlgebra.eigen(A::Chain{V,G,<:Chain{V,G}}) where {V,G}
Expand Down
104 changes: 104 additions & 0 deletions src/composite.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -216,6 +216,109 @@ function Base.exp(t::T,::Val{hint}) where T<:TensorGraded{V} where {V,hint}
end
end

@inline Base.expm1(A::DyadicChain) = exp(A)-I
@inline Base.exp(A::DyadicChain{V,G,T,1}) where {V,G,T} = Chain{V,G}(Values(Chain{V,G}(exp(A[1][1]))))

@inline function Base.exp(A::DyadicChain{V,G,<:Real,2}) where {V,G}
T = typeof(exp(zero(valuetype(A))))
@inbounds a = A[1][1]
@inbounds c = A[1][2]
@inbounds b = A[2][1]
@inbounds d = A[2][2]
v = (a-d)^2 + 4*b*c
if v > 0
z = sqrt(v)
z1 = cosh(z / 2)
z2 = sinh(z / 2) / z
elseif v < 0
z = sqrt(-v)
z1 = cos(z / 2)
z2 = sin(z / 2) / z
else # if v == 0
z1 = T(1.0)
z2 = T(0.5)
end
r = exp((a + d) / 2)
m11 = r * (z1 + (a - d) * z2)
m12 = r * 2b * z2
m21 = r * 2c * z2
m22 = r * (z1 - (a - d) * z2)
Chain{V,G}(Chain{V,G}(m11, m21), Chain{V,G}(m12, m22))
end

@inline function Base.exp(A::DyadicChain{V,G,<:Complex,2}) where {V,G}
T = typeof(exp(zero(valuetype(A))))
@inbounds a = A[1][1]
@inbounds c = A[1][2]
@inbounds b = A[2][1]
@inbounds d = A[2][2]
z = sqrt((a - d)*(a - d) + 4*b*c )
e = expm1((a + d - z) / 2)
f = expm1((a + d + z) / 2)
ϵ = eps()
g = abs2(z) < ϵ^2 ? exp((a + d) / 2) * (1 + z^2 / 24) : (f - e) / z
m11 = (g * (a - d) + f + e) / 2 + 1
m12 = g * b
m21 = g * c
m22 = (-g * (a - d) + f + e) / 2 + 1
Chain{V,G}(Chain{V,G}(m11, m21), Chain{V,G}(m12, m22))
end

# Adapted from implementation in Base; algorithm from
# Higham, "Functions of Matrices: Theory and Computation", SIAM, 2008
function Base.exp(_A::DyadicChain{W,G,T,N}) where {W,G,T,N}
S = typeof((zero(T)*zero(T) + zero(T)*zero(T))/one(T))
A = Chain{W,G}(map.(S,value(_A)))
# omitted: matrix balancing, i.e., LAPACK.gebal!
nA = maximum(sum.(value.(map.(abs,value(A)))))
# marginally more performant than norm(A, 1)
## For sufficiently small nA, use lower order Padé-Approximations
if (nA <= 2.1)
A2 = A*A
if nA > 0.95
U = S(8821612800)*I+A2*(S(302702400)*I+A2*(S(2162160)*I+A2*(S(3960)*I+A2)))
U = A*U
V = S(17643225600)*I+A2*(S(2075673600)*I+A2*(S(30270240)*I+A2*(S(110880)*I+S(90)*A2)))
elseif nA > 0.25
U = S(8648640)*I+A2*(S(277200)*I+A2*(S(1512)*I+A2))
U = A*U
V = S(17297280)*I+A2*(S(1995840)*I+A2*(S(25200)*I+S(56)*A2))
elseif nA > 0.015
U = S(15120)*I+A2*(S(420)*I+A2)
U = A*U
V = S(30240)*I+A2*(S(3360)*I+S(30)*A2)
else
U = S(60)*I+A2
U = A*U
V = S(120)*I+S(12)*A2
end
expA = (V - U) \ (V + U)
else
s = log2(nA/5.4) # power of 2 later reversed by squaring
if s > 0
si = ceil(Int,s)
A = A / S(2^si)
end
A2 = A*A
A4 = A2*A2
A6 = A2*A4
U = A6*(A6 + S(16380)*A4 + S(40840800)*A2) +
(S(33522128640)*A6 + S(10559470521600)*A4 + S(1187353796428800)*A2) +
S(32382376266240000)*I
U = A*U
V = A6*(S(182)*A6 + S(960960)*A4 + S(1323241920)*A2) +
(S(670442572800)*A6 + S(129060195264000)*A4 + S(7771770303897600)*A2) +
S(64764752532480000)*I
expA = (V - U) \ (V + U)
if s > 0 # squaring to reverse dividing by power of 2
for t=1:si
expA = expA*expA
end
end
end
expA
end

qlog(b::PseudoCouple,x::Int=10000) = qlog(multispin(b),x)
qlog(b::AntiSpinor,x::Int=10000) = qlog(Multivector(b),x)
function qlog(w::T,x::Int=10000) where T<:TensorAlgebra
Expand Down Expand Up @@ -804,6 +907,7 @@ function inv_approx(t::Chain{M,1,<:Chain{V,1}}) where {M,V}
mdims(M) < mdims(V) ? (inv(tt⋅t))⋅tt : tt⋅inv(t⋅tt)
end

Base.:\(t::Chain{M,1,<:Chain{W,1}},v::Chain{V,1,<:Chain}) where {M,W,V} = Chain{V,1}(t.\value(v))
Base.:\(t::Chain{M,1,<:Chain{W,1}},v::Chain{V,1}) where {M,W,V} = value(t)\v
Base.in(v::Chain{V,1},t::Chain{W,1,<:Chain{V,1}}) where {V,W} = v ∈ value(t)
Base.inv(t::Chain{V,1,<:Chain{W,1}}) where {W,V} = inv(value(t))
Expand Down
45 changes: 36 additions & 9 deletions src/multivectors.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,18 @@ end
import Leibniz: Fields, parval, mixed, mvecs, svecs, spinsum, spinsum_set
parsym = (Symbol,parval...)

"""
compactio(::Bool)
Toggles compact numbers for extended `Grassmann` elements (enabled by default)
"""
compact = ( () -> begin
io=true
return (tf=io)->(io≠tf && (io=tf); return io)
end)()

compactio(io::IO) = compact() ? IOContext(io, :compact => true) : io

function showterm(io::IO,V,B::UInt,i::T,compact=get(io,:compact,false)) where T
if !(|(broadcast(<:,T,parsym)...)) && signbit(i) && !isnan(i)
print(io, compact ? "-" : " - ")
Expand Down Expand Up @@ -84,17 +96,24 @@ Chain(v::Chain{V,G,𝕂}) where {V,G,𝕂} = v

DyadicProduct{V,W,G,T,N} = Chain{V,G,Chain{W,G,T,N},N}
DyadicChain{V,G,T,N} = DyadicProduct{V,V,G,T,N}
#TriadicProduct{V,W,G,T,N} = Chain{V,G,DyadicChain{W,G,T,N},N}
#DyadicChain{V,G,T,N} = Chain{V,G,Chain{V,G,T,N},N}
#TriadicChain{V,G,T,N,M} = Chain{V,G,DyadicChain{V,1,T,M},N}

Base.Matrix(m::Chain{V,G,<:TensorGraded{W,G}}) where {V,W,G} = hcat(value.(Chain.(value(m)))...)
Base.Matrix(m::Chain{V,G,<:Chain{W,G}}) where {V,W,G} = hcat(value.(value(m))...)
DyadicChain(m::Matrix) = Chain{Submanifold(size(m)[1]),1}(m)
Chain{V}(m::Matrix) where V = Chain{V,1}(m)
function Chain{V,G}(m::Matrix) where {V,G}
N = size(m)[2]
Chain{V,G,Chain{N≠mdims(V) ? Submanifold(N) : V,G}}(m)
end
Chain{V,G,Chain{W,G}}(m::Matrix) where {V,W,G} = Chain{V,G}(Chain{W,G}.(getindex.(Ref(m),:,list(1,size(m)[2]))))
DyadicChain{V,G}(m::Matrix) where {V,G} = Chain{V,G}(Chain{V,G}.(getindex.(Ref(m),:,list(1,size(m)[2]))))
DyadicChain{V}(m::Matrix) where V = DyadicChain{V,1}(m)
DyadicChain(m::Matrix) = DyadicChain{Submanifold(size(m)[1])}(m)
Base.log(A::DyadicChain{V,G}) where {V,G} = DyadicChain{V,G}(log(Matrix(A)))

export Chain, DyadicProduct, DyadicChain
export Chain, DyadicChain, TriadicChain, DyadicProduct
getindex(m::Chain,i::Int) = m.v[i]
getindex(m::Chain,i::UnitRange{Int}) = m.v[i]
getindex(m::Chain,i::T) where T<:AbstractVector = m.v[i]
Expand All @@ -118,6 +137,7 @@ Base.transpose(t::Chain{V,1,<:Chain{W,1}}) where {V,W} = _transpose(value(t),V)

function show(io::IO, m::Chain{V,G,T}) where {V,G,T}
ib,compact = indexbasis(mdims(V),G),get(io,:compact,false)
io = compactio(io)
@inbounds Leibniz.showvalue(io,V,ib[1],m.v[1])
for k ∈ list(2,binomial(mdims(V),G))
@inbounds showterm(io,V,ib[k],m.v[k],compact)
Expand Down Expand Up @@ -370,6 +390,7 @@ getindex(m::Multivector,i::T) where T<:AbstractVector = m.v[i]
setindex!(m::Multivector{V,T} where V,k::T,i::Int,j::Int) where T = (m[i][j] = k)
Base.firstindex(m::Multivector) = 0
Base.lastindex(m::Multivector{V,T} where T) where V = mdims(V)
@pure Base.length(m::Multivector{V}) where V = 1<<mdims(V)

(m::Multivector{V,T})(g::Val{G}) where {V,T,G} = Chain{V,G,T}(m[g])
(m::Multivector{V,T})(g::Int,i) where {V,T} = m(Val(g),i)
Expand All @@ -388,6 +409,7 @@ end

function show(io::IO, m::Multivector{V,T}) where {V,T}
N,compact,bases = mdims(V),get(io,:compact,false),true
io = compactio(io)
bs = binomsum_set(N)
@inbounds print(io,m.v[1])
for i ∈ list(2,N+1)
Expand Down Expand Up @@ -492,6 +514,7 @@ for pinor ∈ (:Spinor,:AntiSpinor)
#setindex!(m::$pinor{V,T} where V,k::T,i::Int,j::Int) where T = (m[i][j] = k)
Base.firstindex(m::$pinor) = 0
Base.lastindex(m::$pinor{V,T} where T) where V = mdims(V)
@pure Base.length(m::$pinor{V}) where V = 1<<(mdims(V)-1)
grade(m::$pinor,g::Val) = m(g)
(m::$pinor{V,T})(g::Int,i::Int) where {T,V} = m(Val(g),i)
Multivector(t::$pinor{V}) where V = Multivector{V}(t)
Expand Down Expand Up @@ -643,6 +666,7 @@ end

function Base.show(io::IO, m::Spinor{V,T}) where {V,T}
N,compact = mdims(V),get(io,:compact,false)
io = compactio(io)
bs = spinsum_set(N)
@inbounds print(io,m.v[1])
for i ∈ evens(2,N)
Expand All @@ -654,6 +678,7 @@ function Base.show(io::IO, m::Spinor{V,T}) where {V,T}
end
function Base.show(io::IO, m::AntiSpinor{V,T}) where {V,T}
N,compact = mdims(V),get(io,:compact,false)
io = compactio(io)
bs = antisum_set(N)
for i ∈ evens(1,N)
ib = indexbasis(N,i)
Expand Down Expand Up @@ -866,6 +891,7 @@ for couple ∈ (:Couple,:PseudoCouple)
Spinor{V}(val::$couple{V,B,𝕂}) where {V,B,𝕂} = Spinor{V,𝕂}(val)
Base.zero(::$couple{V,B,T}) where {V,B,T} = $couple{V,B}(zero(Complex{T}))
Base.zero(::Type{$couple{V,B,T}}) where {V,B,T} = $couple{V,B}(zero(Complex{T}))
@pure Base.length(m::$couple) = 2
end
end

Expand Down Expand Up @@ -1007,6 +1033,7 @@ getindex(P::Dyadic,i::Int,j::Int) = P.x[i]*P.y[j]
show(io::IO,P::Dyadic) = print(io,"(",P.x,")⊗(",P.y,")")

DyadicChain(P::Dyadic{V}) where V = DyadicProduct{V}(P)
#DyadicChain{V}(P::Dyadic{V}) where V = outer(P.x,P.y)
DyadicChain{V}(P::Dyadic{V}) where V = DyadicProduct{V}(p)
DyadicProduct(P::Dyadic{V}) where V = DyadicProduct{V}(P)
DyadicProduct{V}(P::Dyadic{V}) where V = outer(P.x,P.y)
Expand Down Expand Up @@ -1184,27 +1211,27 @@ end
@pure maxpseudograde(t::Type{<:AntiSpinor{V}}) where V = mdims(V)-1
@pure nextmaxpseudograde(t) = maxpseudograde(t)-nextgrade(t)

Leibniz.gdims(t::Tuple{Vector{<:Chain},Vector{Int}}) = gdims(t[1][findall(x->!iszero(x),t[2])])
function Leibniz.gdims(t::Vector{<:Chain})
Leibniz.count_gdims(t::Tuple{Vector{<:Chain},Vector{Int}}) = count_gdims(t[1][findall(x->!iszero(x),t[2])])
function Leibniz.count_gdims(t::Vector{<:Chain})
out = @inbounds zeros(Variables{mdims(Manifold(points(t)))+1,Int})
@inbounds out[mdims(Manifold(t))+1] = length(t)
return out
end
function Leibniz.gdims(t::Values{N,<:Vector}) where N
function Leibniz.count_gdims(t::Values{N,<:Vector}) where N
out = @inbounds zeros(Variables{mdims(points(t[1]))+1,Int})
for i ∈ list(1,N)
@inbounds out[mdims(Manifold(t[i]))+1] = length(t[i])
end
return out
end
function Leibniz.gdims(t::Values{N,<:Tuple}) where N
function Leibniz.count_gdims(t::Values{N,<:Tuple}) where N
out = @inbounds zeros(Variables{mdims(points(t[1][1]))+1,Int})
for i ∈ list(1,N)
@inbounds out[mdims(Manifold(t[i][1]))+1] = length(t[i][1])
end
return out
end
function Leibniz.gdims(t::Multivector{V}) where V
function Leibniz.count_gdims(t::Multivector{V}) where V
N = mdims(V)
out = zeros(Variables{N+1,Int})
bs = binomsum_set(N)
Expand All @@ -1217,8 +1244,8 @@ function Leibniz.gdims(t::Multivector{V}) where V
return out
end

Leibniz.χ(t::Values{N,<:Vector}) where N = (B=gdims(t);sum([B[t]*(-1)^t for t ∈ list(1,length(B))]))
Leibniz.χ(t::Values{N,<:Tuple}) where N = (B=gdims(t);sum([B[t]*(-1)^t for t ∈ list(1,length(B))]))
Leibniz.χ(t::Values{N,<:Vector}) where N = (B=count_gdims(t);sum([B[t]*(-1)^t for t ∈ list(1,length(B))]))
Leibniz.χ(t::Values{N,<:Tuple}) where N = (B=count_gdims(t);sum([B[t]*(-1)^t for t ∈ list(1,length(B))]))

## Adjoint

Expand Down
7 changes: 4 additions & 3 deletions src/products.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -714,8 +714,8 @@ contraction(a::Proj{V,<:Chain{V,1,<:TensorNested}},b::TensorGraded{V,0}) where V
#contraction(a::Chain{W,1,<:Proj{V}},b::Chain{V,1}) where {W,V} = Chain{W,1}(value(a).⋅b)
contraction(a::Chain{W,1,<:Dyadic{V}},b::Chain{V,1}) where {W,V} = Chain{W,1}(value(a).⋅Ref(b))
contraction(a::Proj{W,<:Chain{W,1,<:TensorNested{V}}},b::Chain{V,1}) where {W,V} = a.v:b
contraction(a::Chain{W},b::Chain{V,G,<:Chain}) where {W,G,V} = Chain{V,G}(value.(Ref(a).⋅value(b)))
contraction(a::Chain{W,L,<:Chain},b::Chain{V,G,<:Chain{W,L}}) where {W,L,G,V} = Chain{V,G}(value.(Ref(a).⋅value(b)))
contraction(a::Chain{W},b::Chain{V,G,<:Chain}) where {W,G,V} = Chain{V,G}(value.(Single.(Ref(a).⋅value(b))))
contraction(a::Chain{W,L,<:Chain,N},b::Chain{V,G,<:Chain{W,L},M}) where {W,L,G,V,N,M} = Chain{V,G}(value.(Ref(a).⋅value(b)))
contraction(a::Multivector{W,<:Multivector},b::Multivector{V,<:Multivector{W}}) where {W,V} = Multivector{V}(column(Ref(a).⋅value(b)))
Base.:(:)(a::Chain{V,1,<:Chain},b::Chain{V,1,<:Chain}) where V = sum(value(a).⋅value(b))
Base.:(:)(a::Chain{W,1,<:Dyadic{V}},b::Chain{V,1}) where {W,V} = sum(value(a).⋅Ref(b))
Expand All @@ -725,7 +725,8 @@ contraction(a::Submanifold{W},b::Chain{V,G,<:Chain}) where {W,G,V} = Chain{V,G}(
contraction(a::Single{W},b::Chain{V,G,<:Chain}) where {W,G,V} = Chain{V,G}(column(Ref(a).⋅value(b)))
contraction(x::Chain{V,G,<:Chain},y::Single{V,G}) where {V,G} = value(y)*x[bladeindex(mdims(V),UInt(basis(y)))]
contraction(x::Chain{V,G,<:Chain},y::Submanifold{V,G}) where {V,G} = x[bladeindex(mdims(V),UInt(y))]
contraction(a::Chain{V,L,<:Chain{V,G}},b::Chain{V,G,<:Chain{V}}) where {V,G,L} = Chain{V,G}(matmul(value(a),value(b)))
contraction(a::Chain{V,L,<:Chain{V,G},N},b::Chain{V,G,<:Chain{V},M}) where {V,G,L,N,M} = Chain{V,G}(contraction.(Ref(a),value(b)))
contraction(x::Chain{V,L,<:Chain{V,G},N},y::Chain{V,G,<:Chain{V,L},N}) where {L,N,V,G} = Chain{V,G}(contraction.(Ref(x),value(y)))
contraction(x::Chain{W,L,<:Chain{V,G},N},y::Chain{V,G,T,N}) where {W,L,N,V,G,T} = Chain{V,G}(matmul(value(x),value(y)))
contraction(x::Chain{W,L,<:Multivector{V},N},y::Chain{V,G,T,N}) where {W,L,N,V,G,T} = Multivector{V}(matmul(value(x),value(y)))
contraction(x::Multivector{W,<:Chain{V,G},N},y::Multivector{V,T,N}) where {W,N,V,G,T} = Chain{V,G}(matmul(value(x),value(y)))
Expand Down

2 comments on commit 90ba554

@chakravala
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register()

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/105641

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.8.18 -m "<description of version>" 90ba554674c92686c463ee6328296a68277be9e5
git push origin v0.8.18

Please sign in to comment.