Skip to content

Commit

Permalink
implemented Betti numbers and Euler characteristic
Browse files Browse the repository at this point in the history
  • Loading branch information
@chakravala
chakravala committed Oct 15, 2019
1 parent d784b4d commit 315d3b4
Show file tree
Hide file tree
Showing 4 changed files with 121 additions and 49 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "Grassmann"
uuid = "4df31cd9-4c27-5bea-88d0-e6a7146666d8"
authors = ["Michael Reed"]
version = "0.3.2"
version = "0.3.3"

[deps]
AbstractLattices = "398f06c4-4d28-53ec-89ca-5b2656b7603d"
Expand Down
68 changes: 56 additions & 12 deletions src/Grassmann.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -410,20 +410,21 @@ end
# ParaAlgebra

using Leibniz
import Leibniz: ∂, d, ∇
export ∇, ∂, d, ↑, ↓
import Leibniz: ∂, d, ∇, Δ
export ∇, Δ, ∂, d, ↑, ↓

generate_products(:(Leibniz.Operator),:svec)

@pure function (W::Signature{N})(d::Leibniz.Derivation{T,O}) where {N,T,O}
O < 1 && (return SChain{Int,W,1}(ones(Int,grade(W))))
C = mixedmode(W)<0
V = diffvars(W)≠0 ? W : tangent(W,O<2 ? 2 : O,C ? Int(ndims(W)/2) : ndims(W))
V = diffvars(W)≠0 ? W : tangent(W,O,C ? Int(ndims(W)/2) : ndims(W))
G,D = grade(V),diffvars(V)==1
G2 = (C ? Int(G/2) : G)-1
∇ = sum([getbasis(V,1<<(D ? G : k+G))*getbasis(V,1<<k) for k ∈ 0:G2])
isone(O) && (return ∇)
x = (∇⋅∇)^div(isodd(O) ? O-1 : O,2)
isodd(O) ? sum([x*(getbasis(V,1<<(k+G))*getbasis(V,1<<k)) for k ∈ 0:G2]) : x
isodd(O) ? sum([(x*getbasis(V,1<<(k+G)))*getbasis(V,1<<k) for k ∈ 0:G2]) : x
end

∂(ω::T) where T<:TensorAlgebra{V} where V = ω⋅V(∇)
Expand Down Expand Up @@ -454,6 +455,47 @@ function ↓(ω::T) where T<:TensorAlgebra{V} where V
end
end

## skeleton

export skeleton

absym(t) = abs(t)
absym(t::Basis) = t
absym(t::T) where T<:TensorTerm{V,G} where {V,G} = SBlade{V,G}(absym(value(t)),basis(t))
absym(t::SChain{T,V,G}) where {T,V,G} = SChain{T,V,G}(absym.(value(t)))
absym(t::MChain{T,V,G}) where {T,V,G} = MChain{T,V,G}(absym.(value(t)))
absym(t::MultiVector{T,V}) where {T,V} = MultiVector{T,V}(absym.(value(t)))

function skeleton(x::T) where T<:TensorTerm{V} where V
X,B = absym(x),bits(basis(x))
count_ones(symmetricmask(V,B,B)[1])>1 ? X + subcomplex(absym(∂(X))) : X
end
skeleton(x::S) where {S<:TensorMixed{T,V} where T} where V = absym(x) + subcomplex(absym(∂(x)))
function subcomplex(x::S) where {S<:TensorMixed{T,V} where T} where V
N,G,g = ndims(V),grade(x),0
ib = indexbasis(N,G)
for k ∈ 1:binomial(N,G)
if !iszero(x.v[k]) && count_ones(symmetricmask(V,ib[k],ib[k])[1]) > 0
g += skeleton(SBlade{V,G}(x.v[k],getbasis(V,ib[k])))
end
end
return g
end
function subcomplex(x::MultiVector{T,V} where T) where V
N,g = ndims(V),0
for i ∈ 2:N
R = binomsum(N,i)
ib = indexbasis(N,i)
for k ∈ 1:binomial(N,i)
if !iszero(x.v[k+R])
G = count_ones(symmetricmask(V,ib[k],ib[k])[1])
G>0 && (g += skeleton(SBlade{V,i}(x.v[k+R],getbasis(V,ib[k]))))
end
end
end
return g
end

function __init__()
@require Reduce="93e0c654-6965-5f22-aba9-9c1ae6b3c259" begin
*(a::Reduce.RExpr,b::Basis{V}) where V = SBlade{V}(a,b)
Expand All @@ -479,7 +521,7 @@ function __init__()
@require GaloisFields="8d0d7f98-d412-5cd4-8397-071c807280aa" generate_algebra(:GaloisFields,:AbstractGaloisField)
@require LightGraphs="093fc24a-ae57-5d10-9952-331d41423f4d" begin
function LightGraphs.SimpleDiGraph(x::T,g=LightGraphs.SimpleDiGraph(grade(V))) where T<:TensorTerm{V} where V
ind = (signbit(value(x)) ? reverse : identity)(Grassmann.indices(basis(x)))
ind = (signbit(value(x)) ? reverse : identity)(indices(basis(x)))
grade(x) == 2 ? LightGraphs.add_edge!(g,ind...) : graph(∂(x),g)
return g
end
Expand All @@ -488,10 +530,10 @@ function __init__()
end
function graph(x::S,g=LightGraphs.SimpleDiGraph(grade(V))) where {S<:TensorMixed{T,V} where T} where V
N,G = ndims(V),grade(x)
ib = Grassmann.indexbasis(N,G)
for k ∈ 1:Grassmann.binomial(N,G)
ib = indexbasis(N,G)
for k ∈ 1:binomial(N,G)
if !iszero(x.v[k])
B = Grassmann.symmetricmask(V,ib[k],ib[k])[1]
B = symmetricmask(V,ib[k],ib[k])[1]
count_ones(B) ≠1 && LightGraphs.SimpleDiGraph(x.v[k]*getbasis(V,B),g)
end
end
Expand All @@ -500,11 +542,11 @@ function __init__()
function graph(x::MultiVector{T,V} where T,g=LightGraphs.SimpleDiGraph(grade(V))) where V
N = ndims(V)
for i ∈ 2:N
R = Grassmann.binomsum(N,i)
ib = Grassmann.indexbasis(N,i)
for k ∈ 1:Grassmann.binomial(N,i)
R = binomsum(N,i)
ib = indexbasis(N,i)
for k ∈ 1:binomial(N,i)
if !iszero(x.v[k+R])
B = Grassmann.symmetricmask(V,ib[k],ib[k])[1]
B = symmetricmask(V,ib[k],ib[k])[1]
count_ones(B) ≠ 1 && LightGraphs.SimpleDiGraph(x.v[k+R]*getbasis(V,B),g)
end
end
Expand All @@ -519,6 +561,8 @@ function __init__()
Base.convert(::Type{GeometryTypes.Point},t::MChain{T,V,G}) where {T,V,G} = G == 1 ? GeometryTypes.Point(value(vector(t))) : GeometryTypes.Point(zeros(T,ndims(V))...)
Base.convert(::Type{GeometryTypes.Point},t::SChain{T,V,G}) where {T,V,G} = G == 1 ? GeometryTypes.Point(value(vector(t))) : GeometryTypes.Point(zeros(T,ndims(V))...)
GeometryTypes.Point(t::T) where T<:TensorAlgebra = convert(GeometryTypes.Point,t)
export points
points(f;V=identity,r=-2π:0.0001:2π) = [Point(V(Grassmann.vector(f(t)))) for t ∈ r]
end
#@require AbstractPlotting="537997a7-5e4e-5d89-9595-2241ea00577e" nothing
#@require Makie="ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" nothing
Expand Down
10 changes: 8 additions & 2 deletions src/algebra.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -535,6 +535,12 @@ Sandwich product: ω>>>η = ω⊖η⊖(~ω)
"""
>>>(x::TensorAlgebra{V},y::TensorAlgebra{V}) where V = x * y * ~x

## linear algebra

export ⟂, ∥

∥(a,b) = iszero(a∧b)

### Product Algebra Constructor

function generate_products(Field=Field,VEC=:mvec,MUL=:*,ADD=:+,SUB=:-,CONJ=:conj,PAR=false)
Expand Down Expand Up @@ -1313,13 +1319,13 @@ function generate_products(Field=Field,VEC=:mvec,MUL=:*,ADD=:+,SUB=:-,CONJ=:conj
end
return MChain{t,V,G}(out)
end
function $op(a::$Chain{T,V,G},b::MultiVector{S,V}) where {T<:$Field,V,G,S}
function $op(a::$Chain{T,V,G},b::MultiVector{S,V}) where {T<:$Field,V,G,S<:$Field}
$(insert_expr((:N,:t,:r,:bng),VEC)...)
out = convert($VEC(N,t),$(bcast(bop,:(copy(value(b,$VEC(N,t))),))))
@inbounds $(add_val(eop,:(out[r+1:r+bng]),:(value(a,$VEC(N,G,t))),bop))
return MultiVector{t,V}(out)
end
function $op(a::MultiVector{T,V},b::$Chain{S,V,G}) where {T<:$Field,V,G,S}
function $op(a::MultiVector{T,V},b::$Chain{S,V,G}) where {T<:$Field,V,G,S<:$Field}
$(insert_expr((:N,:t,:r,:bng),VEC)...)
out = copy(value(a,$VEC(N,t)))
@inbounds $(add_val(eop,:(out[r+1:r+bng]),:(value(b,$VEC(N,G,t))),bop))
Expand Down
90 changes: 56 additions & 34 deletions src/multivectors.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,15 +109,19 @@ for Blade ∈ MSB
@pure $Blade(b::Basis{V,G}) where {V,G} = $Blade{V}(b)
@pure $Blade{V}(b::Basis{V,G}) where {V,G} = $Blade{V,G,b,Int}(1)
$Blade{V}(v::T) where {V,T} = $Blade{V,0,Basis{V}(),T}(v)
$Blade{V}(v::TensorTerm{V}) where V = v
$Blade{V}(v::S) where S<:TensorTerm where V = v
$Blade{V,G,B}(v::T) where {V,G,B,T} = $Blade{V,G,B,T}(v)
$Blade(v,b::TensorTerm{V}) where V = $Blade{V}(v,b)
$Blade{V}(v,b::S) where S<:TensorMixed where V = v*b
$Blade{V}(v::T,b::Basis{V,G}) where {V,G,T} = $Blade{V,G}(v,b)
$Blade(v,b::S) where S<:TensorTerm{V} where V = $Blade{V}(v,b)
$Blade{V}(v,b::S) where S<:TensorAlgebra where V = v*b
$Blade{V}(v,b::Basis{V,G}) where {V,G} = $Blade{V,G}(v,b)
$Blade{V}(v,b::Basis{W,G}) where {V,W,G} = $Blade{V,G}(v,b)
function $Blade{V,G}(v::T,b::Basis{V,G}) where {V,G,T}
order(v)+order(b)>diffmode(V) ? zero(V) : $Blade{V,G,b,T}(v)
end
function $Blade{V,G}(v::T,b::Basis{V,G}) where T<:TensorTerm{V} where {V,G}
function $Blade{V,G}(v::T,b::Basis{W,G}) where {V,W,G,T}
order(v)+order(b)>diffmode(V) ? zero(V) : $Blade{V,G,V(b),T}(v)
end
function $Blade{V,G}(v::T,b::Basis{V,G}) where T<:TensorTerm where {V,G}
order(v)+order(b)>diffmode(V) ? zero(V) : $Blade{V,G,b,Any}(v)
end
function $Blade{V,G,B}(b::T) where T<:TensorTerm{V} where {V,G,B}
Expand Down Expand Up @@ -485,7 +489,7 @@ MultiVector(v::MultiGrade{V}) where V = MultiVector{promote_type(typeval.(v.v)..
import Base: isinf, isapprox
import DirectSum: grade
import AbstractTensors: scalar, involute, unit, even, odd
export basis, grade, hasinf, hasorigin, isorigin, scalar, norm
export basis, grade, hasinf, hasorigin, isorigin, scalar, norm, betti, χ

const VBV = Union{MBlade,SBlade,MChain,SChain,MultiVector}

Expand Down Expand Up @@ -516,6 +520,37 @@ valuetype(t::MultiGrade) = promote_type(valuetype.(terms(t))...)
@pure hasorigin(t::Union{MBlade,SBlade}) = hasorigin(basis(t))
@pure hasorigin(m::TensorAlgebra) = hasorigin(vectorspace(m))

@inline χ(V,b::UInt,t) = iszero(t) ? 0 : isodd(count_ones(symmetricmask(V,b,b)[1])) ? 1 : -1
χ(t::T) where T<:TensorTerm{V,G} where {V,G} = χ(V,bits(basis(t)),t)
χ(t::T) where T<:TensorAlgebra{V} where V = (B=betti(t);sum([B[t]*(-1)^t for t ∈ 1:length(B)]))

function betti(t::T) where T<:TensorTerm{V} where V
B,N = bits(basis(t)),ndims(V)
g = count_ones(symmetricmask(V,B,B)[1])
MVector{N+1,Int}([g==G ? abs(χ(t)) : 0 for G ∈ 0:N])
end
function betti(t::T) where T<:TensorAlgebra{V} where V
N = ndims(V)
out = zeros(MVector{N+1,Int})
ib = indexbasis(N,grade(t))
for k ∈ 1:length(ib)
@inbounds t.v[k] ≠ 0 && (out[count_ones(symmetricmask(V,ib[k],ib[k])[1])+1] += 1)
end
return out
end
function betti(t::MultiVector{T,V} where T) where V
N = ndims(V)
out = zeros(MVector{N+1,Int})
bs = binomsum_set(N)
for G ∈ 0:N
ib = indexbasis(N,G)
for k ∈ 1:length(ib)
@inbounds t.v[k+bs[G+1]] ≠ 0 && (out[count_ones(symmetricmask(V,ib[k],ib[k])[1])+1] += 1)
end
end
return out
end

for A ∈ (:TensorTerm,MSC...), B ∈ (:TensorTerm,MSC...)
@eval isapprox(a::S,b::T) where {S<:$A,T<:$B} = vectorspace(a)==vectorspace(b) && (grade(a)==grade(b) ? norm(a)≈norm(b) : (iszero(a) && iszero(b)))
end
Expand Down Expand Up @@ -606,23 +641,24 @@ end
#@pure supblade(N,S,B) = bladeindex(N,expandbits(N,S,B))
#@pure supmulti(N,S,B) = basisindex(N,expandbits(N,S,B))

@pure function mixed(V::M,ibk::UInt) where M<:Manifold
N,D,VC = ndims(V),diffvars(V),mixedmode(V)
return if D≠0
A,B = ibk&(UInt(1)<<(N-D)-1),ibk&DirectSum.diffmask(V)
VC>0 ? (A<<(N-D))|(B<<N) : A|(B<<(N-D))
else
VC>0 ? ibk<<N : ibk
end
end

@pure function (W::Signature)(b::Basis{V}) where V
V==W && (return b)
!(V⊆W) && throw(error("cannot convert from $(V) to $(W)"))
WC,VC = mixedmode(W),mixedmode(V)
#if ((C1≠C2)&&(C1≥0)&&(C2≥0))
# return V0
if WC<0 && VC≥0
N = ndims(V)
dm = diffvars(V)
if dm≠0
D = DirectSum.diffmask(V)
m = (~D)&bits(b)
d = (D&bits(b))<<(N-dm+(VC>0 ? dm : 0))
return getbasis(W,(VC>0 ? m<<(N-dm) : m)⊻d)
else
return getbasis(W,VC>0 ? bits(b)<<(N-dm) : bits(b))
end
getbasis(W,mixed(V,bits(b)))
elseif WC≥0 && VC≥0
getbasis(W,bits(b))
else
Expand All @@ -646,20 +682,13 @@ for Chain ∈ MSC
WC,VC = mixedmode(W),mixedmode(V)
#if ((C1≠C2)&&(C1≥0)&&(C2≥0))
# return V0
N,M,D = ndims(V),ndims(W),diffvars(V)
N,M = ndims(V),ndims(W)
out = zeros(choicevec(M,G,valuetype(b)))
ib = indexbasis(N,G)
for k ∈ 1:length(ib)
@inbounds if b[k] ≠ 0
if WC<0 && VC≥0
@inbounds ibk = ib[k]
if D≠0
A,B = ibk&(UInt(1)<<(N-D)-1),ibk&((UInt(1)<<D-1)<<(N-D))
ibk = VC>0 ? (A<<D)|(B<<N) : A|(B<<(N-D))
else
ibk = VC>0 ? ibk<<N : ibk
end
@inbounds setblade!(out,b[k],ibk,Dimension{M}())
@inbounds setblade!(out,b[k],mixed(V,ib[k]),Dimension{M}())
elseif WC≥0 && VC≥0
@inbounds setblade!(out,b[k],ib[k],Dimension{M}())
else
Expand Down Expand Up @@ -693,7 +722,7 @@ function (W::Signature)(m::MultiVector{T,V}) where {T,V}
WC,VC = mixedmode(W),mixedmode(V)
#if ((C1≠C2)&&(C1≥0)&&(C2≥0))
# return V0
N,M,D = ndims(V),ndims(W),diffvars(V)
N,M = ndims(V),ndims(W)
out = zeros(choicevec(M,valuetype(m)))
bs = binomsum_set(N)
for i ∈ 1:N+1
Expand All @@ -702,14 +731,7 @@ function (W::Signature)(m::MultiVector{T,V}) where {T,V}
@inbounds s = k+bs[i]
@inbounds if m.v[s] ≠ 0
if WC<0 && VC≥0
@inbounds ibk = ib[k]
if D≠0
A,B = ibk&(UInt(1)<<(N-D)-1),ibk&((UInt(1)<<D-1)<<(N-D))
ibk = VC>0 ? (A<<D)|(B<<N) : A|(B<<(N-D))
else
ibk = VC>0 ? ibk<<N : ibk
end
@inbounds setmulti!(out,m.v[s],ibk,Dimension{M}())
@inbounds setmulti!(out,m.v[s],mixed(V,ib[k]),Dimension{M}())
elseif WC≥0 && VC≥0
@inbounds setmulti!(out,m.v[s],ib[k],Dimension{M}())
else
Expand Down

0 comments on commit 315d3b4

Please sign in to comment.