From e1fb80d93ab0b2d7a2e6d5956e1327237d088214 Mon Sep 17 00:00:00 2001 From: Michael Reed <18372368+chakravala@users.noreply.github.com> Date: Fri, 28 Aug 2020 16:05:48 -0400 Subject: [PATCH] improved dyadic product operations --- Project.toml | 2 - src/Grassmann.jl | 42 ++++++++++----------- src/algebra.jl | 46 ++++++++++++++++------- src/composite.jl | 92 ++++++++++++++++++++++----------------------- src/forms.jl | 2 +- src/multivectors.jl | 60 ++++++++++++++++------------- src/parity.jl | 4 +- 7 files changed, 134 insertions(+), 114 deletions(-) diff --git a/Project.toml b/Project.toml index fe6ea8f..9d70563 100644 --- a/Project.toml +++ b/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] diff --git a/src/Grassmann.jl b/src/Grassmann.jl index 370db42..f78f4ea 100644 --- a/src/Grassmann.jl +++ b/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<!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"; diff --git a/src/algebra.jl b/src/algebra.jl index 22f5359..10b5104 100644 --- a/src/algebra.jl +++ b/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<1 S += term @@ -36,7 +36,7 @@ end term = zeros(mvec(N,t)) S .= B out .= value(b^2)/2 - norms = SizedVector{3}(nb,norm(out),norm(term)) + norms = FixedVector{3}(nb,norm(out),norm(term)) k::Int = 3 @inbounds while (norms[2]1) && k ≤ 10000 S += out @@ -80,7 +80,7 @@ function Base.exp(t::MultiVector,::Val{hint}) where hint end end -@pure isR301(V::DiagonalForm) = DirectSum.diagonalform(V) == SVector(1,1,1,0) +@pure isR301(V::DiagonalForm) = DirectSum.diagonalform(V) == Values(1,1,1,0) @pure isR301(::SubManifold{V}) where V = isR301(V) @pure isR301(V) = false @@ -129,7 +129,7 @@ function qlog(w::T,x::Int=10000) where T<:TensorAlgebra w2,f = w^2,norm(w) prod = w*w2 S,term = w,prod/3 - norms = SizedVector{3}(f,norm(term),f) + norms = FixedVector{3}(f,norm(term),f) k::Int = 5 @inbounds while (norms[2]1) && k ≤ x S += term @@ -156,7 +156,7 @@ end # http://www.netlib.org/cephes/qlibdoc.html#qlog S .= value(b) out .= value(b*w2) term .= out/3 - norms = SizedVector{3}(f,norm(term),f) + norms = FixedVector{3}(f,norm(term),f) k::Int = 5 @inbounds while (norms[2]1) && k ≤ x S += term @@ -197,7 +197,7 @@ function Base.cosh(t::T) where T<:TensorAlgebra τ = t^2 S,term = τ/2,(τ^2)/24 f = norm(S) - norms = SizedVector{3}(f,norm(term),f) + norms = FixedVector{3}(f,norm(term),f) k::Int = 6 @inbounds while norms[2]1 S += term @@ -222,7 +222,7 @@ end term = zeros(mvec(N,t)) S .= value(τ)/2 out .= value((τ^2))/24 - norms = SizedVector{3}(norm(S),norm(out),norm(term)) + norms = FixedVector{3}(norm(S),norm(out),norm(term)) k::Int = 6 @inbounds while (norms[2]1) && k ≤ 10000 S += out @@ -246,7 +246,7 @@ function Base.sinh(t::T) where T<:TensorAlgebra V = Manifold(t) τ,f = t^2,norm(t) S,term = t,(t*τ)/6 - norms = SizedVector{3}(f,norm(term),f) + norms = FixedVector{3}(f,norm(term),f) k::Int = 5 @inbounds while norms[2]1 S += term @@ -271,7 +271,7 @@ end term = zeros(mvec(N,t)) S .= value(b) out .= value(b*τ)/6 - norms = SizedVector{3}(norm(S),norm(out),norm(term)) + norms = FixedVector{3}(norm(S),norm(out),norm(term)) k::Int = 5 @inbounds while (norms[2]1) && k ≤ 10000 S += out @@ -294,7 +294,7 @@ for (logfast,expf) ∈ ((:log_fast,:exp),(:logh_fast,:exph)) @eval function $logfast(t::T) where T<:TensorAlgebra V = Manifold(t) term = zero(V) - norm = SizedVector{2}(0.,0.) + norm = FixedVector{2}(0.,0.) while true en = $expf(term) term -= 2(en-t)/(en+t) @@ -337,12 +337,12 @@ end=# function Cramer(N::Int,j=0) t = j ≠ 0 ? :T : :t - x,y = SVector{N}([Symbol(:x,i) for i ∈ 1:N]),SVector{N}([Symbol(:y,i) for i ∈ 1:N]) + x,y = Values{N}([Symbol(:x,i) for i ∈ 1:N]),Values{N}([Symbol(:y,i) for i ∈ 1:N]) xy = [:(($(x[1+i]),$(y[1+i])) = ($(x[i])∧$t[$(1+i-j)],$t[end-$i]∧$(y[i]))) for i ∈ 1:N-1-j] return x,y,xy end -@generated function Base.:\(t::SVector{M,<:Chain{V,1}},v::Chain{V,1}) where {M,V} +@generated function Base.:\(t::Values{M,<:Chain{V,1}},v::Chain{V,1}) where {M,V} W = M≠mdims(V) ? SubManifold(M) : V; N = M-1 if M == 1 && (V === ℝ1 || V == 1) return :(Chain{V,1}(Values(v[1]/t[1][1]))) @@ -377,16 +377,16 @@ end :(Chain{$W,1}(column($(Expr(:call,:.⋅,out,:(Ref(detx))))./abs2(detx))))) end -@generated function Base.in(v::Chain{V,1},t::SVector{N,<:Chain{V,1}}) where {V,N} +@generated function Base.in(v::Chain{V,1},t::Values{N,<:Chain{V,1}}) where {V,N} if N == mdims(V) x,y,xy = Grassmann.Cramer(N-1) mid = [:(s==signbit(($(x[i])∧v∧$(y[end-i]))[1])) for i ∈ 1:N-2] - out = SVector(:(s==signbit((v∧$(y[end]))[1])),mid...,:(s==signbit(($(x[end])∧v)[1]))) + out = Values(:(s==signbit((v∧$(y[end]))[1])),mid...,:(s==signbit(($(x[end])∧v)[1]))) return Expr(:block,:((x1,y1)=(t[1],t[end])),xy...,:(s=signbit((t[1]∧$(y[end]))[1])),ands(out)) else x,y,xy = Grassmann.Cramer(N-1,1) mid = [:(signscalar(($(x[i])∧(v-x1)∧$(y[end-i]))/d)) for i ∈ 1:N-2] - out = SVector(:(signscalar((v∧∧(vectors(t,v)))/d)),mid...,:(signscalar(($(x[end])∧(v-x1))/d))) + out = Values(:(signscalar((v∧∧(vectors(t,v)))/d)),mid...,:(signscalar(($(x[end])∧(v-x1))/d))) return Expr(:block,:(T=vectors(t)),:((x1,y1)=(t[1],T[end])),xy..., :($(x[end])=$(x[end-1])∧T[end-1];d=$(x[end])∧T[end]),ands(out)) end @@ -398,11 +398,11 @@ end M > mdims(V) && (return :(tt = _transpose(t,$W); tt⋅inv(Chain{$W,1}(t)⋅tt))) x,y,xy = Grassmann.Cramer(N) val = if iseven(N) - Expr(:call,:SVector,y[end],[:($(y[end-i])∧$(x[i])) for i ∈ 1:N-1]...,x[end]) + Expr(:call,:Values,y[end],[:($(y[end-i])∧$(x[i])) for i ∈ 1:N-1]...,x[end]) elseif M≠mdims(V) - Expr(:call,:SVector,y[end],[:($(iseven(i) ? :+ : :-)($(y[end-i])∧$(x[i]))) for i ∈ 1:N-1]...,:(-$(x[end]))) + Expr(:call,:Values,y[end],[:($(iseven(i) ? :+ : :-)($(y[end-i])∧$(x[i]))) for i ∈ 1:N-1]...,:(-$(x[end]))) else - Expr(:call,:SVector,:(-$(y[end])),[:($(isodd(i) ? :+ : :-)($(y[end-i])∧$(x[i]))) for i ∈ 1:N-1]...,x[end]) + Expr(:call,:Values,:(-$(y[end])),[:($(isodd(i) ? :+ : :-)($(y[end-i])∧$(x[i]))) for i ∈ 1:N-1]...,x[end]) end out = if M≠mdims(V) :(vector.($(Expr(:call,:./,val,:((t[1]∧$(y[end]))))))) @@ -412,16 +412,16 @@ end return Expr(:block,:((x1,y1)=(t[1],t[end])),xy...,:(_transpose($out,$W))) end -@generated function grad(T::SVector{M,<:Chain{V,1}}) where {M,V} +@generated function grad(T::Values{M,<:Chain{V,1}}) where {M,V} W = M≠mdims(V) ? SubManifold(M) : V; N = mdims(V)-1 M < mdims(V) && (return :(ct = Chain{$W,1}(T); map(↓(V),ct⋅inv(_transpose(T,$W)⋅ct)))) x,y,xy = Grassmann.Cramer(N) val = if iseven(N) - Expr(:call,:SVector,[:($(y[end-i])∧$(x[i])) for i ∈ 1:N-1]...,x[end]) + Expr(:call,:Values,[:($(y[end-i])∧$(x[i])) for i ∈ 1:N-1]...,x[end]) elseif M≠mdims(V) - Expr(:call,:SVector,y[end],[:($(iseven(i) ? :+ : :-)($(y[end-i])∧$(x[i]))) for i ∈ 1:N-1]...,:(-$(x[end]))) + Expr(:call,:Values,y[end],[:($(iseven(i) ? :+ : :-)($(y[end-i])∧$(x[i]))) for i ∈ 1:N-1]...,:(-$(x[end]))) else - Expr(:call,:SVector,[:($(isodd(i) ? :+ : :-)($(y[end-i])∧$(x[i]))) for i ∈ 1:N-1]...,x[end]) + Expr(:call,:Values,[:($(isodd(i) ? :+ : :-)($(y[end-i])∧$(x[i]))) for i ∈ 1:N-1]...,x[end]) end out = if M≠mdims(V) :(vector.($(Expr(:call,:./,val,:((t[1]∧$(y[end]))))))) @@ -431,7 +431,7 @@ end return Expr(:block,:(t=_transpose(T,$W)),:((x1,y1)=(t[1],t[end])),xy...,:(_transpose($out,↓(V)))) end -@generated function Base.:\(t::SVector{N,<:Chain{M,1}},v::Chain{V,1}) where {N,M,V} +@generated function Base.:\(t::Values{N,<:Chain{M,1}},v::Chain{V,1}) where {N,M,V} W = M≠mdims(V) ? SubManifold(N) : V if mdims(M) > mdims(V) :(ct=Chain{$W,1}(t); ct⋅(inv(_transpose(t,$W)⋅ct)⋅v)) @@ -453,7 +453,7 @@ INV(m::Chain{V,1,<:Chain{V,1}}) where V = Chain{V,1,Chain{V,1}}(inv(SMatrix(m))) export vandermonde @generated approx(x,y::Chain{V}) where V = :(polynom(x,$(Val(mdims(V))))⋅y) -approx(x,y::SVector{N}) where N = value(polynom(x,Val(N)))⋅y +approx(x,y::Values{N}) where N = value(polynom(x,Val(N)))⋅y approx(x,y::AbstractVector) = [x^i for i ∈ 0:length(y)-1]⋅y vandermonde(x::Array,y::Array,N::Int) = vandermonde(x,N)\y[:] # compute ((inv(X'*X))*X')*y @@ -468,8 +468,8 @@ end vandermonde(x,y,V) = (length(x)≠mdims(V) ? _vandermonde(x,V) : vandermonde(x,V))\y vandermonde(x,V) = transpose(_vandermonde(x,V)) _vandermonde(x::Chain,V) = _vandermonde(value(x),V) -@generated _vandermonde(x::SVector{N},V) where N = :(Chain{$(SubManifold(N)),1}(polynom.(x,$(Val(mdims(V)))))) -@generated polynom(x,::Val{N}) where N = Expr(:call,:(Chain{$(SubManifold(N)),1}),Expr(:call,:SVector,[:(x^$i) for i ∈ 0:N-1]...)) +@generated _vandermonde(x::Values{N},V) where N = :(Chain{$(SubManifold(N)),1}(polynom.(x,$(Val(mdims(V)))))) +@generated polynom(x,::Val{N}) where N = Expr(:call,:(Chain{$(SubManifold(N)),1}),Expr(:call,:Values,[:(x^$i) for i ∈ 0:N-1]...)) function vandermondeinterp(x,y,V,grid) # grid=384 coef = vandermonde(x,y,V) # Vandermonde ((inv(X'*X))*X')*y @@ -486,11 +486,11 @@ end quote p = points(t) M,A = ↓(Manifold(p)),p[c[1]] - Chain{M,1}.($(Expr(:.,:SVector,v))) + Chain{M,1}.($(Expr(:.,:Values,v))) end end -@pure list(a::Int,b::Int) = SVector{max(0,b-a+1),Int}(a:b...) -vectors(x::SVector{N,<:Chain{V}},y=x[1]) where {N,V} = ↓(V).(x[list(2,N)].-y) +@pure list(a::Int,b::Int) = Values{max(0,b-a+1),Int}(a:b...) +vectors(x::Values{N,<:Chain{V}},y=x[1]) where {N,V} = ↓(V).(x[list(2,N)].-y) vectors(x::Chain{V,1},y=x[1]) where V = vectors(value(x),y) #point(x,y=x[1]) = y∧∧(vectors(x)) @@ -538,12 +538,13 @@ volumes(m) = mdims(Manifold(m))≠2 ? volumes(m,detsimplex(m)) : edgelength.(val detsimplex(m::Vector{<:Chain{V}}) where V = ∧(m)/factorial(mdims(V)-1) detsimplex(m::ChainBundle) = detsimplex(value(m)) mean(m::T) where T<:AbstractVector{<:Chain} = sum(m)/length(m) -mean(m::T) where T<:SVector = sum(m)/length(m) +mean(m::T) where T<:Values = sum(m)/length(m) mean(m::Chain{V,1,<:Chain} where V) = mean(value(m)) -barycenter(m::SVector{N,<:Chain}) where N = (s=sum(m);s/s[1]) +barycenter(m::Values{N,<:Chain}) where N = (s=sum(m);s/s[1]) barycenter(m::Vector{<:Chain}) = (s=sum(m);s/s[1]) barycenter(m::Chain{V,1,<:Chain} where V) = barycenter(value(m)) -curl(m::SVector{N,<:Chain{V}} where N) where V = curl(Chain{V,1}(m)) +curl(m::FixedVector{N,<:Chain{V}} where N) where V = curl(Chain{V,1}(m)) +curl(m::Values{N,<:Chain{V}} where N) where V = curl(Chain{V,1}(m)) curl(m::T) where T<:TensorAlgebra = Manifold(m)(∇)×m LinearAlgebra.det(t::Chain{V,1,<:Chain} where V) = ∧(t) LinearAlgebra.det(m::Vector{<:Chain{V}}) where V = ∧(m) @@ -581,12 +582,12 @@ function refinemesh!(::R,p::ChainBundle{W},e,t,η,_=nothing) where {W,R<:Abstrac p = points(t) x,T,V = value(p),value(t),Manifold(p) for i ∈ η - push!(x,Chain{V,1}(SVector(1,(x[i+1][2]+x[i][2])/2))) + push!(x,Chain{V,1}(Values(1,(x[i+1][2]+x[i][2])/2))) end sort!(x,by=x->x[2]); submesh!(p) - e[end] = Chain{p(2),1}(SVector(length(x))) + e[end] = Chain{p(2),1}(Values(length(x))) for i ∈ length(t)+2:length(x) - push!(T,Chain{p,1}(SVector{2,Int}(i-1,i))) + push!(T,Chain{p,1}(Values{2,Int}(i-1,i))) end end @@ -649,25 +650,20 @@ if VERSION >= v"1.4" (T = promote_type(TA, Chain{V,G,𝕂,X}); mul!(similar(B, T, (size(transA, 1), size(B, 2))), transA, B, 1, 0)) end -#=@generated function StaticArrays._diff(::Size{S}, a::SVector{Q,<:Chain}, ::Val{D}) where {S,D,Q} - N = length(S) - Snew = ([n==D ? S[n]-1 : S[n] for n = 1:N]...,) - +@generated function AbstractTensors._diff(::Val{N}, a::Values{Q,<:Chain}, ::Val{1}) where {N,Q} + Snew = N-1 exprs = Array{Expr}(undef, Snew) - itr = [1:n for n = Snew] - - for i1 = Base.product(itr...) + for i1 = Base.product(1:Snew) i2 = copy([i1...]) - i2[D] = i1[D] + 1 + i2[1] = i1[1] + 1 exprs[i1...] = :(a[$(i2...)] - a[$(i1...)]) end - return quote Base.@_inline_meta - T = eltype(a) - @inbounds return similar_type(a, T, Size($Snew))(tuple($(exprs...))) + elements = tuple($(exprs...)) + @inbounds return AbstractTensors.similar_type(a, eltype(a), Val($Snew))(elements) end -end=# +end Base.map(fn, x::MultiVector{V}) where V = MultiVector{V}(map(fn, value(x))) Base.map(fn, x::Chain{V,G}) where {V,G} = Chain{V,G}(map(fn,value(x))) @@ -715,7 +711,7 @@ Base.iterate(t::Orthogrid,state) = (s=state+1; s≤length(t) ? (getindex(t,s),s) @pure Base.lastindex(t::Orthogrid) = length(t) @pure Base.lastindex(t::Orthogrid,i::Int) = size(t)[i] @pure Base.getindex(t::Orthogrid,i::CartesianIndex) = getindex(t,i.I...) -@pure Base.getindex(t::Orthogrid{V},i::Vararg{Int}) where V = Chain{V,1}(value(t.x.min)+(SVector(i)-1).*step(t)) +@pure Base.getindex(t::Orthogrid{V},i::Vararg{Int}) where V = Chain{V,1}(value(t.x.min)+(Values(i)-1).*step(t)) Base.IndexStyle(::Orthogrid) = IndexCartesian() function Base.getindex(A::Orthogrid, I::Int) diff --git a/src/forms.jl b/src/forms.jl index 2ac14cc..3e0b9ff 100644 --- a/src/forms.jl +++ b/src/forms.jl @@ -40,7 +40,7 @@ return :b elseif W⊆V if G == 1 - ind = SVector{mdims(W),Int}(indices(bits(W),mdims(V))) + ind = Values{mdims(W),Int}(indices(bits(W),mdims(V))) :(@inbounds Chain{w,1,T}(b.v[$ind])) else quote out,N = zeros(choicevec(M,G,valuetype(b))),mdims(V) diff --git a/src/multivectors.jl b/src/multivectors.jl index e885dbd..daf9598 100644 --- a/src/multivectors.jl +++ b/src/multivectors.jl @@ -21,7 +21,7 @@ export UniformScaling, I, points ## Chain{V,G,𝕂} @computed struct Chain{V,G,𝕂} <: TensorGraded{V,G} - v::SVector{binomial(mdims(V),G),𝕂} + v::Values{binomial(mdims(V),G),𝕂} end """ @@ -30,27 +30,29 @@ end Chain type with pseudoscalar `V::Manifold`, grade/rank `G::Int`, scalar field `𝕂::Type`. """ Chain{V,G}(val::S) where {V,G,S<:AbstractVector{𝕂}} where 𝕂 = Chain{V,G,𝕂}(val) -#Chain{V,G}(args::𝕂...) where {V,G,𝕂} = Chain{V,G}(SVector{binomial(mdims(V),G)}(args...)) +Chain{V}(val::S) where {V,S<:TupleVector{N,𝕂}} where {N,𝕂} = Chain{V,1,𝕂}(val) +Chain(val::S) where S<:TupleVector{N,𝕂} where {N,𝕂} = Chain{SubManifold(N),1,𝕂}(val) +#Chain{V,G}(args::𝕂...) where {V,G,𝕂} = Chain{V,G}(Values{binomial(mdims(V),G)}(args...)) @generated function Chain{V,G}(args::𝕂...) where {V,G,𝕂} bg = binomial(mdims(V),G) - ref = SVector{bg}([:(args[$i]) for i ∈ 1:bg]) - :(Chain{V,G}($(Expr(:call,:(SVector{$bg,𝕂}),ref...)))) + ref = Values{bg}([:(args[$i]) for i ∈ 1:bg]) + :(Chain{V,G}($(Expr(:call,:(Values{$bg,𝕂}),ref...)))) end @generated function Chain{V}(args::𝕂...) where {V,𝕂} - bg = mdims(V); ref = SVector{bg}([:(args[$i]) for i ∈ 1:bg]) - :(Chain{V,1}($(Expr(:call,:(SVector{$bg,𝕂}),ref...)))) + bg = mdims(V); ref = Values{bg}([:(args[$i]) for i ∈ 1:bg]) + :(Chain{V,1}($(Expr(:call,:(Values{$bg,𝕂}),ref...)))) end @generated function Chain(args::𝕂...) where 𝕂 N = length(args) V = SubManifold(N) - ref = SVector{N}([:(args[$i]) for i ∈ 1:N]) - :(Chain{$V,1}($(Expr(:call,:(SVector{$N,𝕂}),ref...)))) + ref = Values{N}([:(args[$i]) for i ∈ 1:N]) + :(Chain{$V,1}($(Expr(:call,:(Values{$N,𝕂}),ref...)))) end -Chain(v::Chain{V,G,𝕂}) where {V,G,𝕂} = Chain{V,G}(SVector{binomial(mdims(V),G),𝕂}(v.v)) -Chain{𝕂}(v::Chain{V,G}) where {V,G,𝕂} = Chain{V,G}(SVector{binomial(mdims(V),G),𝕂}(v.v)) +Chain(v::Chain{V,G,𝕂}) where {V,G,𝕂} = Chain{V,G}(Values{binomial(mdims(V),G),𝕂}(v.v)) +Chain{𝕂}(v::Chain{V,G}) where {V,G,𝕂} = Chain{V,G}(Values{binomial(mdims(V),G),𝕂}(v.v)) export Chain getindex(m::Chain,i::Int) = m.v[i] @@ -64,7 +66,7 @@ Base.firstindex(m::Chain) = 1 Base.zero(::Type{<:Chain{V,G,T}}) where {V,G,T} = Chain{V,G}(zeros(svec(mdims(V),G,T))) Base.zero(::Chain{V,G,T}) where {V,G,T} = Chain{V,G}(zeros(svec(mdims(V),G,T))) -transpose_row(t::SVector{N,<:Chain{V}},i,W=V) where {N,V} = Chain{W,1}(getindex.(t,i)) +transpose_row(t::Values{N,<:Chain{V}},i,W=V) where {N,V} = Chain{W,1}(getindex.(t,i)) transpose_row(t::FixedVector{N,<:Chain{V}},i,W=V) where {N,V} = Chain{W,1}(getindex.(t,i)) transpose_row(t::Chain{V,1,<:Chain},i) where V = transpose_row(value(t),i,V) @generated _transpose(t::Values{N,<:Chain{V,1}},W=V) where {N,V} = :(Chain{V,1}(transpose_row.(Ref(t),$(Values{mdims(V)}(1:mdims(V))),W))) @@ -134,6 +136,12 @@ end ==(a::T,b::Chain{V}) where T<:TensorTerm{V} where V = b==a ==(a::Chain{V},b::T) where T<:TensorTerm{V} where V = prod(0==value(b).==value(a)) +Base.ones(::Type{Chain{V,G,T,X}}) where {V,G,T,X} = Chain{V,G,T}(ones(Values{X,T})) +Base.ones(::Type{Chain{V,G,T,X}}) where {V,G,T<:Chain,X} = Chain{V,G,T}(ones.(ntuple(n->T,mdims(V)))) +⊗(a::Type{<:Chain{V}},b::Type{<:Chain{W}}) where {V,W} = Chain{V,1,Chain{W,1,Float64,mdims(W)},mdims(V)} +⊗(a::Type{<:Chain{V,1}},b::Type{<:Chain{W,1}}) where {V,W} = Chain{V,1,Chain{W,1,Float64,mdims(W)},mdims(V)} +⊗(a::Type{<:Chain{V,1}},b::Type{<:Chain{W,1,T}}) where {V,W,T} = Chain{V,1,Chain{W,1,T,mdims(W)},mdims(V)} + """ ChainBundle{V,G,P} <: Manifold{V} <: TensorAlgebra{V} @@ -156,7 +164,7 @@ end @pure bundle(::ChainBundle{V,G,T,P} where {V,G,T}) where P = P @pure deletebundle!(V) = deletebundle!(bundle(V)) @pure function deletebundle!(P::Int) - bundle_cache[P] = [Chain{ℝ^0,0,Int}(SVector(0))] + bundle_cache[P] = [Chain{ℝ^0,0,Int}(Values(0))] end @pure isbundle(::ChainBundle) = true @pure isbundle(t) = false @@ -204,7 +212,7 @@ Base.show(io::IO,m::ChainBundle) = print(io,showbundle(m),length(m)) ## MultiVector{V,𝕂} @computed struct MultiVector{V,T} <: TensorMixed{V} - v::SVector{1<!iszero(x),t[2])]) function Leibniz.gdims(t::Vector{<:Chain}) - out = zeros(MVector{mdims(Manifold(points(t)))+1,Int}) + out = zeros(Variables{mdims(Manifold(points(t)))+1,Int}) out[mdims(Manifold(t))+1] = length(t) return out end -function Leibniz.gdims(t::SVector{N,<:Vector}) where N - out = zeros(MVector{mdims(points(t[1]))+1,Int}) +function Leibniz.gdims(t::Values{N,<:Vector}) where N + out = zeros(Variables{mdims(points(t[1]))+1,Int}) for i ∈ 1:N out[mdims(Manifold(t[i]))+1] = length(t[i]) end return out end -function Leibniz.gdims(t::SVector{N,<:Tuple}) where N - out = zeros(MVector{mdims(points(t[1][1]))+1,Int}) +function Leibniz.gdims(t::Values{N,<:Tuple}) where N + out = zeros(Variables{mdims(points(t[1][1]))+1,Int}) for i ∈ 1:N out[mdims(Manifold(t[i][1]))+1] = length(t[i][1]) end @@ -494,7 +502,7 @@ function Leibniz.gdims(t::SVector{N,<:Tuple}) where N end function Leibniz.gdims(t::MultiVector{V}) where V N = mdims(V) - out = zeros(MVector{N+1,Int}) + out = zeros(Variables{N+1,Int}) bs = binomsum_set(N) for G ∈ 0:N ib = indexbasis(N,G) @@ -505,8 +513,8 @@ function Leibniz.gdims(t::MultiVector{V}) where V return out end -Leibniz.χ(t::SVector{N,<:Vector}) where N = (B=gdims(t);sum([B[t]*(-1)^t for t ∈ 1:length(B)])) -Leibniz.χ(t::SVector{N,<:Tuple}) where N = (B=gdims(t);sum([B[t]*(-1)^t for t ∈ 1:length(B)])) +Leibniz.χ(t::Values{N,<:Vector}) where N = (B=gdims(t);sum([B[t]*(-1)^t for t ∈ 1:length(B)])) +Leibniz.χ(t::Values{N,<:Tuple}) where N = (B=gdims(t);sum([B[t]*(-1)^t for t ∈ 1:length(B)])) ## Adjoint diff --git a/src/parity.jl b/src/parity.jl index 1a2b341..4d0ef28 100644 --- a/src/parity.jl +++ b/src/parity.jl @@ -205,8 +205,8 @@ export odd, even, angular, radial, ₊, ₋, ǂ @pure signbit(V::T) where T<:Manifold{N} where N = (ib=indexbasis(N); parity.(Ref(V),ib,ib)) @pure signbit(V::T,G) where T<:Manifold{N} where N = (ib=indexbasis(N,G); parity.(Ref(V),ib,ib)) -@pure angular(V::T) where T<:Manifold = SVector(findall(signbit(V))...) -@pure radial(V::T) where T<:Manifold = SVector(findall(.!signbit(V))...) +@pure angular(V::T) where T<:Manifold = Values(findall(signbit(V))...) +@pure radial(V::T) where T<:Manifold = Values(findall(.!signbit(V))...) @pure angular(V::T,G) where T<:Manifold = findall(signbit(V,G)) @pure radial(V::T,G) where T<:Manifold = findall(.!signbit(V,G))