improved vectorfield and finalized JuliaCon paper, fixed #37

Project.toml

@@ -23,6 +23,9 @@ AbstractTensors = "0.3"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+GeometryTypes = "4d00f742-c7ba-57c2-abde-4428a4b178cb"
+LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
+GraphPlot = "a2cc645c-3eea-5389-862e-a155d0052231"
 
 [targets]
 test = ["Test"]

paper/ref.bib

@@ -15,6 +15,19 @@ @article{chappell-iqbal-gunn-abbott
 	year = {2011},
 }
 
+@article{taylor-algebra,
+	author = {Aaron D. Schutte},
+	title = {A nilpotent algebra approach to Lagrangian mechanics and constrained motion},
+	year = {2016},
+}
+
+@article{doran-hestenes-sommen-acker,
+	author = {Chris Doran, David Hestenes, F. Sommen, and N. van Acker},
+	title = {Lie groups as spin groups},
+	journal = {J. Math. Phys.},
+	year = {1993},
+}
+
 @article{ivancevic,
 	author = {Vladimir and Tijana Ivancevic},
 	title = {Undergraduate Lecture Notes in De Rahm-Hodge Theory},
@@ -37,43 +50,49 @@ @article{hestenes
 }
 
 @book{garling,
-	author = {D.J.H. Garling},
+	author = {Garling},
 	title = {Clifford Algebras: An Introduction},
 	year = {2011},
 }
 
 @book{sobczyk,
 	author = {Garret Sobczyk},
 	title = {New Foundations in Mathematics: The Geometric Concept of Number},
+	publisher = {Springer},
 	year = {2013},
 }
 
 @book{browne,
 	author = {John Browne},
 	title = {Grassmann Algebra, Volume 1: Foundations},
+	publisher = {Barnard Publishing},
 	year = {2011},
 }
 
 @book{grassmann-2,
 	author = {Hermann Grassmann},
 	title = {Extension Theory (Ausdehnungslehre 1862)},
+	publisher = {AMS},
 	year = {2000},
 }
 
 @book{artin,
 	author = {Emil Artin},
 	title = {Geometric Algebra},
+	publisher = {Interscience},
 	year = {1957},
 }
 
 @book{shahshahani,
 	author = {Siavash Shahshahani},
 	title = {An Introductory Course on Differentiable Manifolds},
+	publisher = {Dover},
 	year = {2016},
 }
 
 @book{bishop-goldberg,
 	author = {Richard L. Bishop and Samuel I. Goldberg},
 	title = {Tensor Analysis on Manifolds},
+	publisher = {Macmillan},
 	year = {1968},
 }

src/Grassmann.jl

@@ -428,17 +428,15 @@ end
 ∂(ω::T) where T<:TensorAlgebra{V} where V = ω⋅V(∇)
 d(ω::T) where T<:TensorAlgebra{V} where V = V(∇)∧ω
 
-@pure ℙ(V) = ((i,o)=(hasinf(V),hasorigin(V));i+o==2 ? V : (i+o==0 ? S"∞∅"⊕V : V))
-
 function ↑(ω::T) where T<:TensorAlgebra{V} where V
-    PV = ℙ(V)
-    G = Λ(PV)
-    return if hasinf(PV) && hasorigin(PV)
+    !(hasinf(V)||hasorigin(V)) && (return ω)
+    G = Λ(V)
+    return if hasinf(V) && hasorigin(V)
         ((G.v∞/2)*ω^2+G.v∅)+ω
     else
         ω2 = ω^2
         iω2 = inv(ω2+1)
-        (hasinf(PV) ? G.v∞ : G.v∅)*(ω2-1)*iω2 + 2*iω2*ω
+        (hasinf(V) ? G.v∞ : G.v∅)*(ω2-1)*iω2 + 2*iω2*ω
     end
 end
 function ↑(ω,b)
@@ -453,12 +451,12 @@ function ↑(ω,p,m)
 end
 
 function ↓(ω::T) where T<:TensorAlgebra{V} where V
-    PV = ℙ(V)
-    G = Λ(PV)
-    return if hasinf(PV) && hasorigin(PV)
+    !(hasinf(V)||hasorigin(V)) && (return ω)
+    G = Λ(V)
+    return if hasinf(V) && hasorigin(V)
         inv(G.v∞∅)*(G.v∞∅∧ω)/(-ω⋅G.v∞)
     else
-        b = hasinf(PV) ? G.v∞ : G.v∅
+        b = hasinf(V) ? G.v∞ : G.v∅
         ((ω∧b)*b)/(1-b⋅ω)
     end
 end
@@ -494,8 +492,9 @@ function chain(t::S,::Val{T}=Val{true}()) where S<:TensorTerm{V} where {V,T}
 end
 path(t) = chain(t,Val{false}())
 
-𝒫(t::T) where T<:TensorAlgebra = skeleton(t,Val{false}())
-subcomplex(x::S,v=Val{true}()) where S<:TensorAlgebra = skeleton(absym(∂(x)),v)
+@inline (::Leibniz.Derivation)(x::T,v=Val{true}()) where T<:TensorAlgebra = skeleton(x,v)
+𝒫(t::T) where T<:TensorAlgebra = Δ(t,Val{false}())
+subcomplex(x::S,v=Val{true}()) where S<:TensorAlgebra = Δ(absym(∂(x)),v)
 function skeleton(x::S,v::Val{T}=Val{true}()) where S<:TensorTerm{V} where {V,T}
     B = bits(basis(x))
     count_ones(symmetricmask(V,B,B)[1])>0 ? absym(x)+skeleton(absym(∂(x)),v) : (T ? g_zero(V) : absym(x))
@@ -601,10 +600,10 @@ 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)
-        @pure ptype(::Point{N,T} where N) where T = T
+        @pure ptype(::GeometryTypes.Point{N,T} where N) where T = T
         export points, vectorfield
         points(f,V=identity;r=-2π:0.0001:2π) = [GeometryTypes.Point(V(vector(f(t)))) for t ∈ r]
-        vectorfield(t,V=vectorspace(t)) = p->GeometryTypes.Point(V(vector(SChain{ptype(p),V,1}(p.data)⊘t)))
+        vectorfield(t,V=vectorspace(t),W=V) = p->GeometryTypes.Point(V(vector(↓(↑((V∪vectorspace(t))(SChain{ptype(p),W,1}(p.data)))⊘t))))
     end
     #@require AbstractPlotting="537997a7-5e4e-5d89-9595-2241ea00577e" nothing
     #@require Makie="ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" nothing

src/algebra.jl

@@ -526,7 +526,7 @@ export ⊘
 
 Sandwich product: ω⊘η = (~ω)⊖η⊖ω
 """
-⊘(x::TensorAlgebra{V},y::TensorAlgebra{V}) where V = inv(y) * x * involute(y)
+⊘(x::TensorAlgebra{V},y::TensorAlgebra{V}) where V = diffvars(V)≠0 ? (~y)*x*involute(y) : inv(y)*x*involute(y)
 ⊘(x::TensorAlgebra{V},y::TensorAlgebra{W}) where {V,W} = interop(⊘,x,y)
 
 """

src/multivectors.jl

@@ -671,6 +671,8 @@ end
 #@pure supblade(N,S,B) = bladeindex(N,expandbits(N,S,B))
 #@pure supmulti(N,S,B) = basisindex(N,expandbits(N,S,B))
 
+@pure subvert(::SubManifold{M,V,S} where {M,V}) where S = S
+
 @pure function mixed(V::M,ibk::UInt) where M<:Manifold
     N,D,VC = ndims(V),diffvars(V),mixedmode(V)
     return if D≠0
@@ -687,10 +689,11 @@ end
     WC,VC = mixedmode(W),mixedmode(V)
     #if ((C1≠C2)&&(C1≥0)&&(C2≥0))
     #    return V0
+    B = typeof(V)<:SubManifold ? expandbits(ndims(W),subvert(V),bits(b)) : bits(b)
     if WC<0 && VC≥0
-        getbasis(W,mixed(V,bits(b)))
+        getbasis(W,mixed(V,B))
     elseif WC≥0 && VC≥0
-        getbasis(W,bits(b))
+        getbasis(W,B)
     else
         throw(error("arbitrary Manifold intersection not yet implemented."))
     end
@@ -717,10 +720,11 @@ for Chain ∈ MSC
             ib = indexbasis(N,G)
             for k ∈ 1:length(ib)
                 @inbounds if b[k] ≠ 0
+                    @inbounds B = typeof(V)<:SubManifold ? expandbits(M,subvert(V),ib[k]) : ib[k]
                     if WC<0 && VC≥0
-                        @inbounds setblade!(out,b[k],mixed(V,ib[k]),Dimension{M}())
+                        @inbounds setblade!(out,b[k],mixed(V,B),Dimension{M}())
                     elseif WC≥0 && VC≥0
-                        @inbounds setblade!(out,b[k],ib[k],Dimension{M}())
+                        @inbounds setblade!(out,b[k],B,Dimension{M}())
                     else
                         throw(error("arbitrary Manifold intersection not yet implemented."))
                     end
@@ -737,7 +741,7 @@ for Chain ∈ MSC
             for k ∈ 1:length(ib)
                 @inbounds if b[k] ≠ 0
                     @inbounds if count_ones(ib[k]&S) == G
-                        @inbounds setblade!(out,b[k],lowerbits(N,S,ib[k]),Dimension{M}())
+                        @inbounds setblade!(out,b[k],lowerbits(M,S,ib[k]),Dimension{M}())
                     end
                 end
             end
@@ -760,10 +764,11 @@ function (W::Signature)(m::MultiVector{T,V}) where {T,V}
         for k ∈ 1:length(ib)
             @inbounds s = k+bs[i]
             @inbounds if m.v[s] ≠ 0
+                @inbounds B = typeof(V)<:SubManifold ? expandbits(M,subvert(V),ib[k]) : ib[k]
                 if WC<0 && VC≥0
-                    @inbounds setmulti!(out,m.v[s],mixed(V,ib[k]),Dimension{M}())
+                    @inbounds setmulti!(out,m.v[s],mixed(V,B),Dimension{M}())
                 elseif WC≥0 && VC≥0
-                    @inbounds setmulti!(out,m.v[s],ib[k],Dimension{M}())
+                    @inbounds setmulti!(out,m.v[s],B,Dimension{M}())
                 else
                     throw(error("arbitrary Manifold intersection not yet implemented."))
                 end
@@ -822,3 +827,7 @@ Base.copysign(x::Simplex{V,G,B,T},y::Simplex{V,G,B,T}) where {V,G,B,T} = Simplex
     end
     return A
 end
+
+# genfun
+
+

src/utilities.jl

@@ -244,7 +244,7 @@ const lowerbits_extra = Dict{UInt,Dict{UInt,UInt}}[]
     end
 end
 
-#=const expandbits_cache = Dict{UInt,Dict{UInt,UInt}}[]
+const expandbits_cache = Dict{UInt,Dict{UInt,UInt}}[]
 @pure expandbits_calc(N,S,B) = bit2int(indexbits(N,indices(S,N)[indices(B,N)]))
 @pure function expandbits(N,S,B)
     for k ∈ length(expandbits_cache)+1:N
@@ -253,6 +253,30 @@ end
     @inbounds !haskey(expandbits_cache[N],S) && push!(expandbits_cache[N],S=>Dict{UInt,UInt}())
     @inbounds !haskey(expandbits_cache[N][S],B) && push!(expandbits_cache[N][S],B=>expandbits_calc(N,S,B))
     @inbounds expandbits_cache[N][S][B]
+end
+
+#=const expandbits_cache = Vector{Vector{UInt}}[]
+const expandbits_extra = Dict{UInt,Dict{UInt,UInt}}[]
+@pure expandbits_calc(N,S,B,k=indices(S,N)) = bit2int(indexbits(N,k[indices(B,N)]))
+@pure function expandbits(N,S,B)
+    if N>cache_limit
+        n = N-cache_limit
+        for k ∈ length(expandbits_extra)+1:n
+            push!(expandbits_extra,Dict{UInt,Dict{UInt,UInt}}())
+        end
+        @inbounds !haskey(expandbits_extra[n],S) && push!(expandbits_extra[n],S=>Dict{UInt,UInt}())
+        @inbounds !haskey(expandbits_extra[n][S],B) && push!(expandbits_extra[n][S],B=>expandbits_calc(N,S,B))
+        @inbounds expandbits_extra[n][S][B]
+    else
+        for k ∈ length(expandbits_cache)+1:min(N,cache_limit)
+            push!(expandbits_cache,Vector{Int}[])
+        end
+        for s ∈ length(expandbits_cache[N])+1:S
+            k = indices(S,N)
+            push!(expandbits_cache[N],[expandbits_calc(N,s,d,k) for d ∈ UInt(0):UInt(1)<<(N+1)-1])
+        end
+        @inbounds expandbits_cache[N][S][B+1]
+    end
 end=#