From 473b3488cb76119de30333e83cac9fb7392d4b4c Mon Sep 17 00:00:00 2001 From: Michael Reed <18372368+chakravala@users.noreply.github.com> Date: Thu, 17 Sep 2020 19:01:38 -0400 Subject: [PATCH] reconfigured tensor dispatch #50 #81 --- Project.toml | 2 +- src/algebra.jl | 47 +++++++++++++++++++++++++--- src/composite.jl | 31 +++++++----------- src/forms.jl | 50 ++++++++++++++++------------- src/multivectors.jl | 76 +++++++++++++++++++++++++++++++++++++++++---- 5 files changed, 155 insertions(+), 51 deletions(-) diff --git a/Project.toml b/Project.toml index 11e252f..d4f12e7 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Grassmann" uuid = "4df31cd9-4c27-5bea-88d0-e6a7146666d8" authors = ["Michael Reed"] -version = "0.7.1" +version = "0.7.2" [deps] AbstractTensors = "a8e43f4a-99b7-5565-8bf1-0165161caaea" diff --git a/src/algebra.jl b/src/algebra.jl index 64d437d..792bcbf 100644 --- a/src/algebra.jl +++ b/src/algebra.jl @@ -105,7 +105,10 @@ end export ∧, ∨, ⊗ -⊗(a::A,b::B) where {A<:TensorAlgebra,B<:TensorAlgebra} = a∧b +#⊗(a::A,b::B) where {A<:TensorAlgebra,B<:TensorAlgebra} = a∧b +⊗(a::A,b::B) where {A<:TensorGraded,B<:TensorGraded} = Dyadic(a,b) +⊗(a::A,b::B) where {A<:TensorGraded,B<:TensorGraded{V,0} where V} = a*b +⊗(a::A,b::B) where {A<:TensorGraded{V,0} where V,B<:TensorGraded} = a*b ## regressive product: (L = grade(a) + grade(b); (-1)^(L*(L-mdims(V)))*⋆(⋆(a)∧⋆(b))) @@ -202,20 +205,56 @@ 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::Proj,b::TensorGraded) = a.v⊗(a.v⋅b) +contraction(a::Dyadic,b::TensorGraded) = a.x⊗(a.y⋅b) +contraction(a::TensorGraded,b::Dyadic) = (a⋅b.x)⊗b.y +contraction(a::TensorGraded,b::Proj) = (a⋅b.v)⊗b.v +contraction(a::Dyadic,b::Dyadic) = (a.x*(a.y⋅b.x))⊗b.y +contraction(a::Dyadic,b::Proj) = (a.x*(a.y⋅b.v))⊗b.v +contraction(a::Proj,b::Dyadic) = (a.v*(a.v⋅b.x))⊗b.y +contraction(a::Proj,b::Proj) = (a.v*(a.v⋅b.v))⊗b.v +contraction(a::Dyadic{V},b::TensorGraded{V,0}) where V = Dyadic{V}(a.x*b,a.y) +contraction(a::Proj{V},b::TensorGraded{V,0}) where V = valuetype(b)<:Complex ? Proj{V}(a.v*sqrt(b)) : Dyadic{V}(a.v*b,a.v) +contraction(a::Proj{V,<:Chain{V,1,<:TensorNested}},b::TensorGraded{V,0}) where V = Proj(Chain{V,1}(contraction.(value(a.v),b))) +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,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)) +Base.:(:)(a::Chain{W,1,<:Dyadic{V}},b::Chain{V,1}) where {W,V} = sum(value(a).⋅Ref(b)) +Base.:(:)(a::Chain{W,1,<:Proj{V}},b::Chain{V,1}) where {W,V} = sum(broadcast(⋅,value(a),Ref(b))) + ++(a::Proj{V}...) where V = Proj(Chain(a...)) ++(a::Dyadic{V}...) where V = Proj(Chain(a...)) ++(a::TensorNested{V}...) where V = Proj(Chain(Dyadic.(a)...)) ++(a::Proj{W,<:Chain{W,1,<:TensorNested{V}}} where W,b::TensorNested{V}) where V = +(value(a.v)...,b) ++(a::TensorNested{V},b::Proj{W,<:Chain{W,1,<:TensorNested{V}}} where W) where V = +(a,value(b.v)...) ++(a::Proj{M,<:Chain{M,1,<:TensorNested{V}}} where M,b::Proj{W,<:Chain{W,1,<:TensorNested{V}}} where W) where V = +(value(a.v)...,value(b.v)...) + +-(a::TensorNested) where V = -1a +-(a::TensorNested,b::TensorNested) where V = a+(-b) +*(a::Number,b::TensorNested{V}) where V = (a*one(V))*b +*(a::TensorNested{V},b::Number) where V = a*(b*one(V)) +@inline *(a::TensorGraded{V,0},b::TensorNested{V}) where V = b⋅a +@inline *(a::TensorNested{V},b::TensorGraded{V,0}) where V = a⋅b +@inline *(a::TensorGraded{V,0},b::Proj{V,<:Chain{V,1,<:TensorNested}}) where V = Proj{V}(a*b.v) +@inline *(a::Proj{V,<:Chain{V,1,<:TensorNested}},b::TensorGraded{V,0}) where V = Proj{V}(a.v*b) +Base.:∘(a::A,b::B) where {A<:TensorAlgebra,B<:TensorAlgebra} = a⋅b # dyadic identity element -Base.:+(g::Chain{V,1,<:Chain{V,1}},t::LinearAlgebra.UniformScaling{Bool}) where V = t+g +Base.:+(t::LinearAlgebra.UniformScaling,g::TensorNested) = t+DyadicChain(g) +Base.:+(g::TensorNested,t::LinearAlgebra.UniformScaling) = DyadicChain(g)+t Base.:+(g::Chain{V,1,<:Chain{V,1}},t::LinearAlgebra.UniformScaling) where V = t+g -Base.:-(g::Chain{V,1,<:Chain{V,1}},t::LinearAlgebra.UniformScaling{Bool}) where V = t+g -Base.:-(g::Chain{V,1,<:Chain{V,1}},t::LinearAlgebra.UniformScaling) where V = t+g +Base.:-(t::LinearAlgebra.UniformScaling,g::TensorNested) = t-DyadicChain(g) +Base.:-(g::TensorNested,t::LinearAlgebra.UniformScaling) = DyadicChain(g)-t @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.:-(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)]))) ## cross product diff --git a/src/composite.jl b/src/composite.jl index 4e7fe38..a097980 100644 --- a/src/composite.jl +++ b/src/composite.jl @@ -90,7 +90,7 @@ end function Base.exp(t::T) where T<:TensorGraded S,B = T<:SubManifold,T<:TensorTerm if B && isnull(t) - return one(V) + return one(Manifold(t)) elseif isR301(Manifold(t)) && grade(t)==2 # && abs(t[0])<1e-9 && !options.over u = sqrt(abs(abs2(t)[1])) u<1e-5 && (return 1+t) @@ -294,12 +294,12 @@ for (logfast,expf) ∈ ((:log_fast,:exp),(:logh_fast,:exph)) @eval function $logfast(t::T) where T<:TensorAlgebra V = Manifold(t) term = zero(V) - norm = FixedVector{2}(0.,0.) + nrm = FixedVector{2}(0.,0.) while true en = $expf(term) term -= 2(en-t)/(en+t) - @inbounds norm .= (norm[2],norm(term)) - @inbounds norm[1] ≈ norm[2] && break + @inbounds nrm .= (nrm[2],norm(term)) + @inbounds nrm[1] ≈ nrm[2] && break end return term end @@ -678,39 +678,32 @@ Base.rand(::AbstractRNG,::SamplerType{MultiVector}) = rand(MultiVector{rand(Mani Base.rand(::AbstractRNG,::SamplerType{MultiVector{V}}) where V = MultiVector{V}(DirectSum.orand(svec(mdims(V),Float64))) Base.rand(::AbstractRNG,::SamplerType{MultiVector{V,T}}) where {V,T} = MultiVector{V}(rand(svec(mdims(V),T))) -export Orthotope, Orthogrid +export Orthogrid -struct Orthotope{V,T} - min::Chain{V,1,T} - max::Chain{V,1,T} -end - -(::Base.Colon)(min::Chain{V,1,T},max::Chain{V,1,T}) where {V,T} = Orthotope{V,T}(min,max) - -struct Orthogrid{V,T} # <: TensorGraded{V,1} mess up collect? - x::Orthotope{V,T} +@computed struct Orthogrid{V,T} # <: TensorGraded{V,1} mess up collect? + v::Dyadic{V,Chain{V,1,T,mdims(V)},Chain{V,1,T,mdims(V)}} n::Chain{V,1,Int} s::Chain{V,1,Float64} end -Orthogrid{V,T}(x,n) where {V,T} = Orthogrid{V,T}(x,n,Chain{V,1}(value(x.max-x.min)./(value(n)-1))) +Orthogrid{V,T}(v,n) where {V,T} = Orthogrid{V,T}(v,n,Chain{V,1}(value(v.x-v.y)./(value(n)-1))) -Base.show(io::IO,t::Orthogrid) = println('(',t.x.min,"):(",t.s,"):(",t.x.max,')') +Base.show(io::IO,t::Orthogrid) = println('(',t.v.x,"):(",t.s,"):(",t.v.y,')') zeroinf(f) = iszero(f) ? Inf : f -(::Base.Colon)(min::Chain{V,1,T},step::Chain{V,1,T},max::Chain{V,1,T}) where {V,T} = Orthogrid{V,T}(min:max,Chain{V,1}(Int.(round.(value(max-min)./zeroinf.(value(step))))+1),step) +(::Base.Colon)(min::Chain{V,1,T},step::Chain{V,1,T},max::Chain{V,1,T}) where {V,T} = Orthogrid{V,T}(min⊗max,Chain{V,1}(Int.(round.(value(max-min)./zeroinf.(value(step)))).+1),step) Base.iterate(t::Orthogrid) = (getindex(t,1),1) Base.iterate(t::Orthogrid,state) = (s=state+1; s≤length(t) ? (getindex(t,s),s) : nothing) @pure Base.eltype(::Type{Orthogrid{V,T}}) where {V,T} = Chain{V,1,T,mdims(V)} @pure Base.step(t::Orthogrid) = value(t.s) -@pure Base.size(t::Orthogrid) = value(t.n).data +@pure Base.size(t::Orthogrid) = value(t.n).v @pure Base.length(t::Orthogrid) = prod(size(t)) @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)+(Values(i)-1).*step(t)) +@pure Base.getindex(t::Orthogrid{V},i::Vararg{Int}) where V = Chain{V,1}(value(t.v.x)+(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 3cfc8b7..0ce514b 100644 --- a/src/forms.jl +++ b/src/forms.jl @@ -20,24 +20,27 @@ W = Manifold(w) if isbasis(W) if Q == V - if G == M == 1 - x = bits(W) - X = isdyadic(V) ? x>>Int(mdims(V)/2) : x - Y = 0≠X ? X : x - out = :(@inbounds b.v[bladeindex($(mdims(V)),Y)]) - return :(Simplex{V}(V[intlog(Y)+1] ? -($out) : $out,SubManifold{V}())) - elseif G == 1 && M == 2 - (!isdyadic(V)) && :(throw(error("wrong basis"))) - ib,(m1,m2) = indexbasis(N,1),DirectSum.eval_shift(W) - :(@inbounds $(V[m2] ? :(-(b.v[m2])) : :(b.v[m2]))*getbasis(V,ib[m1])) + if isdyadic(V) + if G == M == 1 + x = bits(W) + X = isdyadic(V) ? x>>Int(mdims(V)/2) : x + Y = 0≠X ? X : x + out = :(@inbounds b.v[bladeindex($(mdims(V)),Y)]) + return :(Simplex{V}(V[intlog(Y)+1] ? -($out) : $out,SubManifold{V}())) + elseif G == 1 && M == 2 + ib,(m1,m2) = indexbasis(N,1),DirectSum.eval_shift(W) + :(@inbounds $(V[m2] ? :(-(b.v[m2])) : :(b.v[m2]))*getbasis(V,ib[m1])) + else + :(throw(error("not yet possible"))) + end else - :(throw(error("not yet possible"))) + return :(contraction(w,b)) end else - :(interform(w,b)) + return :(interform(w,b)) end elseif V==W - return :b + return V===W ? :b : :(Chain{w,G,T}(value(b))) elseif W⊆V if G == 1 ind = Values{mdims(W),Int}(indices(bits(W),mdims(V))) @@ -85,9 +88,10 @@ end (W::Signature)(b::MultiVector{V,T}) where {V,T} = SubManifold(W)(b) function (W::SubManifold{Q,M,S})(m::MultiVector{V,T}) where {Q,M,V,S,T} if isbasis(W) - throw(error("MultiVector forms not yet supported")) + isdyadic(V) && throw(error("MultiVector forms not yet supported")) + return V==W ? contraction(W,m) : interform(W,m) elseif V==W - return m + return V===W ? m : MultiVector{W,T}(value(m)) elseif W⊆V out,N = zeros(choicevec(M,valuetype(m))),mdims(V) bs = binomsum_set(N) @@ -185,9 +189,10 @@ end ## Chain forms (a::Chain)(b::T) where {T<:TensorAlgebra} = interform(a,b) +(a::Chain{V,1,<:Manifold} where V)(b::T) where {T<:TensorAlgebra} = contraction(a,b) @eval begin function (a::Chain{V,2})(b::Chain{V,1}) where V - (!isdyadic(V)) && throw(error("wrong basis")) + (!isdyadic(V)) && (return contraction(a,b)) $(insert_expr((:N,:t,:df,:di))...) out = zero(mvec(N,1,t)) for Q ∈ 1:Int(N/2) @@ -198,7 +203,7 @@ end end return Chain{V,1}(out) end - function Chain{V,T}(b::Matrix{T}) where {V,T} + function Chain{V,T}(b::AbstractMatrix{T}) where {V,T} (!isdyadic(V)) && throw(error("$V does not support this conversion")) $(insert_expr((:N,:M))...) size(b) ≠ (M,M) && throw(error("dimension mismatch")) @@ -216,6 +221,7 @@ end # more forms function (a::Chain{V,1})(b::SubManifold{V,1}) where V + (!isdyadic(V)) && (return contraction(a,b)) x = bits(b) X = isdyadic(V) ? x<2^mdims(V) ? x : X @@ -224,7 +230,7 @@ function (a::Chain{V,1})(b::SubManifold{V,1}) where V end @eval begin function (a::Chain{V,2,T})(b::SubManifold{V,1}) where {V,T} - (!isdyadic(V)) && throw(error("wrong basis")) + (!isdyadic(V)) && (return contraction(a,b)) $(insert_expr((:N,:df,:di))...) Q = bladeindex(N,bits(b)) @inbounds m,val = df[Q][1],df[Q][2]*value(b) @@ -235,6 +241,7 @@ end return Chain{V,1}(out) end function (a::Chain{V,1})(b::Simplex{V,1}) where V + (!isdyadic(V)) && (return contraction(a,b)) $(insert_expr((:t,))...) x = bits(b) X = isdyadic(V) ? x<>Int(mdims(V)/2) : x @@ -251,13 +259,13 @@ end Simplex{V}((a.v*V[intlog(Y)+1]*out)::t,SubManifold{V}()) end function (a::Simplex{V,2})(b::Chain{V,1}) where V - (!isdyadic(V)) && throw(error("wrong basis")) + (!isdyadic(V)) && (return contraction(a,b)) $(insert_expr((:N,:t))...) ib,(m1,m2) = indexbasis(N,1),DirectSum.eval_shift(a) @inbounds ((V[m2]*a.v*b.v[m2])::t)*getbasis(V,ib[m1]) end function (a::Chain{V,2})(b::Simplex{V,1}) where V - (!isdyadic(V)) && throw(error("wrong basis")) + (!isdyadic(V)) && (return contraction(a,b)) $(insert_expr((:N,:t,:df,:di))...) Q = bladeindex(N,bits(b)) out = zero(mvec(N,1,T)) @@ -268,6 +276,7 @@ end return Chain{V,1}(out) end function (a::Chain{V,1})(b::Chain{V,1}) where V + (!isdyadic(V)) && (return contraction(a,b)) $(insert_expr((:N,:M,:t,:df))...) out = zero(t) for Q ∈ 1:M @@ -276,4 +285,3 @@ end return Simplex{V}(out::t,SubManifold{V}()) end end - diff --git a/src/multivectors.jl b/src/multivectors.jl index 3b18478..7711a59 100644 --- a/src/multivectors.jl +++ b/src/multivectors.jl @@ -7,6 +7,16 @@ export TensorTerm, TensorGraded, TensorMixed, SubManifold, Simplex, MultiVector, import AbstractTensors: TensorTerm, TensorGraded, TensorMixed import Leibniz: grade +export TensorNested +abstract type TensorNested{V} <: Manifold{V} end + +for op ∈ (:(Base.:+),:(Base.:-)) + @eval begin + $op(a::A,b::B) where {A<:TensorNested,B<:TensorAlgebra} = $op(DyadicChain(a),b) + $op(a::A,b::B) where {A<:TensorAlgebra,B<:TensorNested} = $op(a,DyadicChain(b)) + end +end + # symbolic print types import Leibniz: Fields, parval, mixed @@ -54,7 +64,10 @@ end 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 +DyadicProduct{V,W,T,N} = Chain{V,1,Chain{W,1,T,N},N} +DyadicChain{V,T,N} = DyadicProduct{V,V,T,N} + +export Chain, DyadicProduct, DyadicChain 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] @@ -69,16 +82,18 @@ Base.zero(::Chain{V,G,T}) where {V,G,T} = Chain{V,G}(zeros(svec(mdims(V),G,T))) 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))) -@generated _transpose(t::FixedVector{N,<:Chain{V,1}},W=V) where {N,V} = :(Chain{V,1}(transpose_row.(Ref(t),$(Values{mdims(V)}(1:mdims(V))),W))) +@generated _transpose(t::Values{N,<:Chain{V,1}},W=V) where {N,V} = :(Chain{V,1}(transpose_row.(Ref(t),$(list(1,mdims(V))),W))) +@generated _transpose(t::FixedVector{N,<:Chain{V,1}},W=V) where {N,V} = :(Chain{V,1}(transpose_row.(Ref(t),$(list(1,mdims(V))),W))) Base.transpose(t::Chain{V,1,<:Chain{V,1}}) where V = _transpose(value(t)) Base.transpose(t::Chain{V,1,<:Chain{W,1}}) where {V,W} = _transpose(value(t),V) +showparens(tmv) = (!(tmv<:TensorTerm||tmv<:Projector)) && |(broadcast(<:,tmv,parval)...) + function show(io::IO, m::Chain{V,G,T}) where {V,G,T} ib = indexbasis(mdims(V),G) @inbounds tmv = typeof(m.v[1]) if |(broadcast(<:,tmv,parsym)...) - par = (!(tmv<:TensorTerm)) && |(broadcast(<:,tmv,parval)...) + par = showparens(tmv) @inbounds par ? print(io,"(",m.v[1],")") : print(io,m.v[1]) else @inbounds print(io,m.v[1]) @@ -88,7 +103,7 @@ function show(io::IO, m::Chain{V,G,T}) where {V,G,T} @inbounds mvs = m.v[k] tmv = typeof(mvs) if |(broadcast(<:,tmv,parsym)...) - par = (!(tmv<:TensorTerm)) && |(broadcast(<:,tmv,parval)...) + par = showparens(tmv) par ? print(io," + (",mvs,")") : print(io," + ",mvs) else sbm = signbit(mvs) @@ -263,7 +278,7 @@ function show(io::IO, m::MultiVector{V,T}) where {V,T} @inbounds if mvs ≠ 0 tmv = typeof(mvs) if |(broadcast(<:,tmv,parsym)...) - par = (!(tmv<:TensorTerm)) && |(broadcast(<:,tmv,parval)...) + par = showparens(tmv) par ? print(io," + (",mvs,")") : print(io," + ",mvs) else sba = signbit(mvs) @@ -317,6 +332,55 @@ end getindex(m::MultiVector,i::T) where T<:AbstractVector{<:SubManifold} = getindex.(m,i) getindex(m::MultiVector{V},i::SubManifold{V}) where V = m[basisindex(mdims(V),UInt(i))] +# Dyadic + +export Projector, Dyadic, Proj + +struct Projector{V,T} <: TensorNested{V} + v::T + Projector{V,T}(v::T) where {T<:Manifold{V}} where V = new{V,T}(v) + Projector{V}(v::T) where {T<:Manifold{V}} where V = new{V,T}(v) +end + +const Proj = Projector + +Proj(v::T) where T<:TensorGraded{V} where V = Proj{V}(v/abs(v)) +Proj(v::Chain{V,1,<:TensorNested}) where V = Proj{V}(v) + +(P::Projector)(x) = contraction(P,x) + +getindex(P::Proj,i::Int,j::Int) = P.v[i]*P.v[j] +getindex(P::Proj{V,<:Chain{V,1,<:TensorNested}} where V,i::Int,j::Int) = sum(getindex.(value(P.v),i,j)) + +show(io::IO,P::Projector) = print(io,"Proj(",P.v,")") + +DyadicChain{V,T}(P::Proj{V,T}) where {V,T} = outer(P.v,P.v) +DyadicChain{V,T}(P::Proj{V,T}) where {V,T<:Chain{V,1,<:TensorNested}} = sum(DyadicChain.(value(P.v))) +DyadicChain{V}(P::Proj{V,T}) where {V,T} = DyadicChain{V,T}(P) +DyadicChain(P::Proj{V,T}) where {V,T} = DyadicChain{V,T}(P) + +struct Dyadic{V,X,Y} <: TensorNested{V} + x::X + y::Y + Dyadic{V,X,Y}(x::X,y::Y) where {X<:TensorGraded,Y<:TensorGraded{V}} where V = new{V,X,Y}(x,y) + Dyadic{V}(x::X,y::Y) where {X<:TensorGraded,Y<:TensorGraded{V}} where V = new{V,X,Y}(x,y) +end + +Dyadic(x::X,y::Y) where {X<:TensorGraded,Y<:TensorGraded{V}} where V = Dyadic{V}(x,y) +Dyadic(P::Projector) = Dyadic(P.v,P.v) +Dyadic(D::Dyadic) = D + +(P::Dyadic)(x) = contraction(P,x) + +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 = DyadicProduct{V}(p) +DyadicProduct(P::Dyadic{V}) where V = DyadicProduct{V}(P) +DyadicProduct{V}(P::Dyadic{V}) where V = outer(P.x,P.y) + ## Generic import Base: isinf, isapprox