improved dyadic product operations

Project.toml

@@ -10,7 +10,6 @@ DirectSum = "22fd7b30-a8c0-5bf2-aabe-97783860d07c"
 Leibniz = "edad4870-8a01-11e9-2d75-8f02e448fc59"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
-StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
@@ -20,7 +19,6 @@ Leibniz = "0.1"
 DirectSum = "0.7"
 AbstractTensors = "0.5.2"
 ComputedFieldTypes = "0.1"
-StaticArrays = "0"
 Requires = "1"
 
 [extras]

src/Grassmann.jl

@@ -5,7 +5,7 @@ module Grassmann
 
 using SparseArrays, ComputedFieldTypes
 using AbstractTensors, Leibniz, DirectSum, Requires
-import AbstractTensors: SVector, MVector, SizedVector, clifford
+import AbstractTensors: Values, Variables, FixedVector, clifford
 
 export ⊕, ℝ, @V_str, @S_str, @D_str, Manifold, SubManifold, Signature, DiagonalForm, value
 export @basis, @basis_str, @dualbasis, @dualbasis_str, @mixedbasis, @mixedbasis_str, Λ
@@ -18,7 +18,7 @@ import Leibniz: Bits, bit2int, indexbits, indices, diffvars, diffmask
 import Leibniz: symmetricmask, indexstring, indexsymbol, combo, digits_fast
 
 import Leibniz: hasconformal, hasinf2origin, hasorigin2inf, g_zero, g_one, bits
-import AbstractTensors: valuetype, scalar, isscalar
+import AbstractTensors: valuetype, scalar, isscalar, ⊗
 import AbstractTensors: vector, isvector, bivector, isbivector, volume, isvolume
 
 ## cache
@@ -62,7 +62,7 @@ end
 ⊘(x::Derivation,y::T) where T<:TensorAlgebra{V} where V = V(x)⊘y
 
 @pure function (V::Signature{N})(d::Leibniz.Derivation{T,O}) where {N,T,O}
-    (O<1||diffvars(V)==0) && (return Chain{V,1,Int}(ones(Int,N)))
+    (O<1||diffvars(V)==0) && (return Chain{V,1,Int}(d.v.λ*ones(Values{N,Int})))
     G,D,C = grade(V),diffvars(V)==1,isdyadic(V)
     G2 = (C ? Int(G/2) : G)-1
     ∇ = sum([getbasis(V,1<<(D ? G : k+G))*getbasis(V,1<<k) for k ∈ 0:G2])
@@ -73,7 +73,7 @@ end
 
 @pure function (M::SubManifold{W,N})(d::Leibniz.Derivation{T,O}) where {W,N,T,O}
     V = isbasis(M) ? W : M
-    (O<1||diffvars(V)==0) && (return Chain{V,1,Int}(ones(Int,N)))
+    (O<1||diffvars(V)==0) && (return Chain{V,1,Int}(d.v.λ*ones(Values{N,Int})))
     G,D,C = grade(V),diffvars(V)==1,isdyadic(V)
     G2 = (C ? Int(G/2) : G)-1
     ∇ = sum([getbasis(V,1<<(D ? G : k+G))*getbasis(V,1<<k) for k ∈ 0:G2])
@@ -93,18 +93,18 @@ function boundary_rank(t,d=gdims(t))
     for k ∈ 2:length(out)-1
         @inbounds out[k] = min(out[k],d[k+1])
     end
-    return SVector(out)
+    return Values(out)
 end
 
 function boundary_null(t)
     d = gdims(t)
     r = boundary_rank(t,d)
     l = length(d)
-    out = zeros(MVector{l,Int})
+    out = zeros(Variables{l,Int})
     for k ∈ 1:l-1
         @inbounds out[k] = d[k+1] - r[k]
     end
-    return SVector(out)
+    return Values(out)
 end
 
 """
@@ -116,11 +116,11 @@ function betti(t::T) where T<:TensorAlgebra
     d = gdims(t)
     r = boundary_rank(t,d)
     l = length(d)-1
-    out = zeros(MVector{l,Int})
+    out = zeros(Variables{l,Int})
     for k ∈ 1:l
         @inbounds out[k] = d[k+1] - r[k] - r[k+1]
     end
-    return SVector(out)
+    return Values(out)
 end
 
 @generated function ↑(ω::T) where T<:TensorAlgebra
@@ -277,12 +277,12 @@ function SparseArrays.sparse(t,cols=columns(t))
     return A
 end
 
-edges(t,cols::SVector) = edges(t,adjacency(t,cols))
+edges(t,cols::Values) = edges(t,adjacency(t,cols))
 function edges(t,adj=adjacency(t))
     ndims(t) == 2 && (return t)
     N = ndims(Manifold(t)); M = points(t)(list(N-1,N)...)
     f = findall(x->!iszero(x),LinearAlgebra.triu(adj))
-    [Chain{M,1}(SVector{2,Int}(f[n].I)) for n ∈ 1:length(f)]
+    [Chain{M,1}(Values{2,Int}(f[n].I)) for n ∈ 1:length(f)]
 end
 
 function facetsinterior(t::Vector{<:Chain{V}}) where V
@@ -323,8 +323,8 @@ function faces(t::Vector{<:Chain{V}},h,::Val{N},g=identity) where {V,N}
     N == 0 && (return [Chain{W,1}(list(1,N))],Int[sum(h)])
     out = Chain{W,1,Int,N}[]
     bnd = Int[]
-    vec = zeros(MVector{ndims(V),Int})
-    val = N+1==ndims(V) ? ∂(Manifold(points(t))(list(1,N+1))(I)) : ones(SVector{binomial(ndims(V),N)})
+    vec = zeros(Variables{ndims(V),Int})
+    val = N+1==ndims(V) ? ∂(Manifold(points(t))(list(1,N+1))(I)) : ones(Values{binomial(ndims(V),N)})
     for i ∈ 1:length(t)
         vec[:] = value(t[i])
         par = DirectSum.indexparity!(vec)
@@ -344,8 +344,8 @@ function faces(t::Vector{<:Chain{V}},h,::Val{N},g=identity) where {V,N}
 end
 
 ∂(t::ChainBundle) = ∂(value(t))
-∂(t::SVector{N,<:Tuple}) where N = ∂.(t)
-∂(t::SVector{N,<:Vector}) where N = ∂.(t)
+∂(t::Values{N,<:Tuple}) where N = ∂.(t)
+∂(t::Values{N,<:Vector}) where N = ∂.(t)
 ∂(t::Tuple{Vector{<:Chain},Vector{Int}}) = ∂(t[1],t[2])
 ∂(t::Vector{<:Chain},u::Vector{Int}) = (f=facets(t,u); f[1][findall(x->!iszero(x),f[2])])
 ∂(t::Vector{<:Chain}) = ndims(t)≠3 ? (f=facetsinterior(t); f[1][setdiff(1:length(f[1]),f[2])]) : edges(t,adjacency(t).%2)
@@ -527,9 +527,9 @@ function __init__()
             c,f = GeometryBasics.coordinates(m),GeometryBasics.faces(m)
             s = size(eltype(c))[1]+1; V = SubManifold(ℝ^s) # s
             n = size(eltype(f))[1]
-            p = ChainBundle([Chain{V,1}(SVector{s,Float64}(1.0,k...)) for k ∈ c])
+            p = ChainBundle([Chain{V,1}(Values{s,Float64}(1.0,k...)) for k ∈ c])
             M = s ≠ n ? p(list(s-n+1,s)) : p
-            t = ChainBundle([Chain{M,1}(SVector{n,Int}(k)) for k ∈ f])
+            t = ChainBundle([Chain{M,1}(Values{n,Int}(k)) for k ∈ f])
             return (p,ChainBundle(∂(t)),t)
         end
         @pure ptype(::GeometryBasics.Point{N,T} where N) where T = T
@@ -692,10 +692,10 @@ function __init__()
         end
         function initmesh(tio::TetGen.TetgenIO, command = "Qp")
             r = TetGen.tetrahedralize(tio, command); V = SubManifold(ℝ^4)
-            p = ChainBundle([Chain{V,1}(SVector{4,Float64}(1.0,k...)) for k ∈ r.points])
-            t = Chain{p,1}.(SVector{4,Int}.(r.tetrahedra))
-            e = Chain{p(2,3,4),1}.(SVector{3,Int}.(r.trifaces))
-            # Chain{p(2,3),1}.(SVector{2,Int}.(r.edges)
+            p = ChainBundle([Chain{V,1}(Values{4,Float64}(1.0,k...)) for k ∈ r.points])
+            t = Chain{p,1}.(Values{4,Int}.(r.tetrahedra))
+            e = Chain{p(2,3,4),1}.(Values{3,Int}.(r.trifaces))
+            # Chain{p(2,3),1}.(Values{2,Int}.(r.edges)
             return p,ChainBundle(e),ChainBundle(t)
         end
         function TetGen.tetrahedralize(mesh::ChainBundle, command = "Qp";

src/algebra.jl

@@ -9,8 +9,8 @@ import Leibniz: loworder, isnull
 
 ## mutating operations
 
-@pure tvec(N,G,t::Symbol=:Any) = :(SVector{$(binomial(N,G)),$t})
-@pure tvec(N,t::Symbol=:Any) = :(SVector{$(1<<N),$t})
+@pure tvec(N,G,t::Symbol=:Any) = :(Values{$(binomial(N,G)),$t})
+@pure tvec(N,t::Symbol=:Any) = :(Values{$(1<<N),$t})
 @pure tvec(N,μ::Bool) = tvec(N,μ ? :Any : :t)
 
 import Leibniz: g_one, g_zero, Field, ExprField
@@ -353,10 +353,10 @@ Exterior product as defined by the anti-symmetric quotient Λ≡⊗/~
 @inline ∧(a::X,b::Y) where {X<:TensorAlgebra,Y<:TensorAlgebra} = interop(∧,a,b)
 @inline ∧(a::TensorAlgebra{V},b::UniformScaling{T}) where {V,T<:Field} = a∧V(b)
 @inline ∧(a::UniformScaling{T},b::TensorAlgebra{V}) where {V,T<:Field} = V(a)∧b
-@generated ∧(t::T) where T<:SVector{N} where N = wedges([:(t[$i]) for i ∈ 1:N])
-@generated ∧(t::T) where T<:SizedVector{N} where N = wedges([:(t[$i]) for i ∈ 1:N])
-∧(::SVector{0,<:Chain{V}}) where V = one(V) # ∧() = 1
-∧(::SizedVector{0,<:Chain{V}}) where V = one(V)
+@generated ∧(t::T) where T<:Values{N} where N = wedges([:(t[$i]) for i ∈ 1:N])
+@generated ∧(t::T) where T<:FixedVector{N} where N = wedges([:(t[$i]) for i ∈ 1:N])
+∧(::Values{0,<:Chain{V}}) where V = one(V) # ∧() = 1
+∧(::FixedVector{0,<:Chain{V}}) where V = one(V)
 ∧(t::Chain{V,1,<:Chain} where V) = ∧(value(t))
 ∧(a::X,b::Y,c::Z...) where {X<:TensorAlgebra,Y<:TensorAlgebra,Z<:TensorAlgebra} = ∧(a∧b,c...)
 
@@ -423,10 +423,10 @@ Regressive product as defined by the DeMorgan's law: ∨(ω...) = ⋆⁻¹(∧(
 @inline ∨(a::X,b::Y) where {X<:TensorAlgebra,Y<:TensorAlgebra} = interop(∨,a,b)
 @inline ∨(a::TensorAlgebra{V},b::UniformScaling{T}) where {V,T<:Field} = a∨V(b)
 @inline ∨(a::UniformScaling{T},b::TensorAlgebra{V}) where {V,T<:Field} = V(a)∨b
-@generated ∨(t::T) where T<:SVector = Expr(:call,:∨,[:(t[$k]) for k ∈ 1:length(t)]...)
-@generated ∨(t::T) where T<:SizedVector = Expr(:call,:∨,[:(t[$k]) for k ∈ 1:length(t)]...)
-∨(::SVector{0,<:Chain{V}}) where V = SubManifold(V) # ∨() = I
-∨(::SizedVector{0,<:Chain{V}}) where V = SubManifold(V)
+@generated ∨(t::T) where T<:Values = Expr(:call,:∨,[:(t[$k]) for k ∈ 1:length(t)]...)
+@generated ∨(t::T) where T<:FixedVector = Expr(:call,:∨,[:(t[$k]) for k ∈ 1:length(t)]...)
+∨(::Values{0,<:Chain{V}}) where V = SubManifold(V) # ∨() = I
+∨(::FixedVector{0,<:Chain{V}}) where V = SubManifold(V)
 ∨(t::Chain{V,1,<:Chain} where V) = ∧(value(t))
 ∨(a::X,b::Y,c::Z...) where {X<:TensorAlgebra,Y<:TensorAlgebra,Z<:TensorAlgebra} = ∨(a∨b,c...)
 
@@ -481,9 +481,27 @@ Interior (right) contraction product: ω⋅η = ω∨⋆η
 
 # dyadic products
 
-contraction(a::Chain{W,G,<:Chain},b::Chain{V,1,<:Chain}) where {W,G,V} = Chain{V,1}(a.⋅value(b))
+export outer
+
+outer(a::Leibniz.Derivation,b::Chain{V,1}) where V= outer(V(a),b)
+outer(a::Chain{W},b::Leibniz.Derivation{T,1}) where {W,T} = outer(a,W(b))
+outer(a::Chain{W},b::Chain{V,1}) where {W,V} = Chain{V,1}(a.*value(b))
+
+contraction(a::Chain{W,G},b::Chain{V,1,<:Chain}) where {W,G,V} = Chain{V,1}(column(Ref(a).⋅value(b)))
+contraction(a::Chain{W,G,<:Chain},b::Chain{V,1,<:Chain}) where {W,G,V} = Chain{V,1}(Ref(a).⋅value(b))
 Base.:(:)(a::Chain{V,1,<:Chain},b::Chain{V,1,<:Chain}) where V = sum(value(a).⋅value(b))
 
+# dyadic identity element
+
+@generated Base.:+(g::Chain{V,1,<:Chain{V,1}},t::LinearAlgebra.UniformScaling{Bool}) where V = :(Chain{V,1}(value(g).+$(getalgebra(V).b[Grassmann.list(2,mdims(V)+1)])))
+@generated Base.:+(g::Chain{V,1,<:Chain{V,1}},t::LinearAlgebra.UniformScaling) where V = :(Chain{V,1}(value(g).+t.λ*$(getalgebra(V).b[Grassmann.list(2,mdims(V)+1)])))
+@generated Base.:+(t::LinearAlgebra.UniformScaling{Bool},g::Chain{V,1,<:Chain{V,1}}) where V = :(Chain{V,1}($(getalgebra(V).b[Grassmann.list(2,mdims(V)+1)]).+value(g)))
+@generated Base.:+(t::LinearAlgebra.UniformScaling,g::Chain{V,1,<:Chain{V,1}}) where V = :(Chain{V,1}(t.λ*$(getalgebra(V).b[Grassmann.list(2,mdims(V)+1)]).+value(g)))
+@generated Base.:-(g::Chain{V,1,<:Chain{V,1}},t::LinearAlgebra.UniformScaling{Bool}) where V = :(Chain{V,1}(value(g).-$(getalgebra(V).b[Grassmann.list(2,mdims(V)+1)])))
+@generated Base.:-(g::Chain{V,1,<:Chain{V,1}},t::LinearAlgebra.UniformScaling) where V = :(Chain{V,1}(value(g).-t.λ*$(getalgebra(V).b[Grassmann.list(2,mdims(V)+1)])))
+@generated Base.:-(t::LinearAlgebra.UniformScaling{Bool},g::Chain{V,1,<:Chain{V,1}}) where V = :(Chain{V,1}($(getalgebra(V).b[Grassmann.list(2,mdims(V)+1)]).-value(g)))
+@generated Base.:-(t::LinearAlgebra.UniformScaling,g::Chain{V,1,<:Chain{V,1}}) where V = :(Chain{V,1}(t.λ*$(getalgebra(V).b[Grassmann.list(2,mdims(V)+1)]).-value(g)))
+
 ## cross product
 
 import LinearAlgebra: cross, ×
@@ -794,7 +812,7 @@ for (op,eop) ∈ ((:+,:(+=)),(:-,:(-=)))
                 end
             end
             G = A|B
-            MultiGrade{V,G}(SVector{count_ones(G),TensorGraded{V}}(out))
+            MultiGrade{V,G}(Values{count_ones(G),TensorGraded{V}}(out))
         end
         function $op(a::MultiGrade{V,A},b::T) where T<:TensorGraded{V,B} where {V,A,B}
             N = mdims(V)
@@ -817,9 +835,9 @@ for (op,eop) ∈ ((:+,:(+=)),(:-,:(-=)))
                 end
             end
             G = A|(UInt(1)<<B)
-            MultiGrade{V,G}(SVector{count_ones(G),TensorGraded{V}}(out))
+            MultiGrade{V,G}(Values{count_ones(G),TensorGraded{V}}(out))
         end
-        $op(a::SparseChain{V,A},b::T) where T<:TensorGraded{V,B} where {V,A,B} = MultiGrade{V,(UInt(1)<<A)|(UInt(1)<<B)}(A<B ? SVector(a,b) : SVector(b,a))
+        $op(a::SparseChain{V,A},b::T) where T<:TensorGraded{V,B} where {V,A,B} = MultiGrade{V,(UInt(1)<<A)|(UInt(1)<<B)}(A<B ? Values(a,b) : Values(b,a))
     end
 end