From ae101e9b68710ceae3dc4a71c5fca8454e01de44 Mon Sep 17 00:00:00 2001 From: Benjamin Zanger Date: Tue, 7 Mar 2023 17:13:13 +0100 Subject: [PATCH 1/9] first steps toward operators for higher dimension TensorSpaces --- src/Multivariate/TensorSpace.jl | 8 +++--- src/PDE/KroneckerOperator.jl | 46 ++++++++++++++++++++------------- 2 files changed, 32 insertions(+), 22 deletions(-) diff --git a/src/Multivariate/TensorSpace.jl b/src/Multivariate/TensorSpace.jl index f7efab66..f4f2b9c5 100644 --- a/src/Multivariate/TensorSpace.jl +++ b/src/Multivariate/TensorSpace.jl @@ -327,11 +327,11 @@ dimension(sp::TensorSpace) = mapreduce(dimension,*,sp.spaces) ==(A::TensorSpace{<:NTuple{N,Space}}, B::TensorSpace{<:NTuple{N,Space}}) where {N} = factors(A) == factors(B) -conversion_rule(a::TensorSpace{<:NTuple{2,Space}}, b::TensorSpace{<:NTuple{2,Space}}) = - conversion_type(a.spaces[1],b.spaces[1]) ⊗ conversion_type(a.spaces[2],b.spaces[2]) +conversion_rule(a::TensorSpace{<:NTuple{N,Space}}, b::TensorSpace{<:NTuple{N,Space}}) where {N} = + mapreduce((a,b)->conversion_type(a,b),⊗,a.spaces,b.spaces) -maxspace_rule(a::TensorSpace{<:NTuple{2,Space}}, b::TensorSpace{<:NTuple{2,Space}}) = - maxspace(a.spaces[1],b.spaces[1]) ⊗ maxspace(a.spaces[2],b.spaces[2]) +maxspace_rule(a::TensorSpace{<:NTuple{N,Space}}, b::TensorSpace{<:NTuple{N,Space}}) where {N} = + mapreduce((a,b)->maxspace(a,b),⊗,a.spaces,b.spaces) function spacescompatible(A::TensorSpace{<:NTuple{N,Space}}, B::TensorSpace{<:NTuple{N,Space}}) where {N} _spacescompatible(factors(A), factors(B)) diff --git a/src/PDE/KroneckerOperator.jl b/src/PDE/KroneckerOperator.jl index 635fc996..75dbd75f 100644 --- a/src/PDE/KroneckerOperator.jl +++ b/src/PDE/KroneckerOperator.jl @@ -6,21 +6,31 @@ export KroneckerOperator # KroneckerOperator gives the kronecker product of two 1D operators ######### -struct KroneckerOperator{S,V,DS,RS,DI,RI,T} <: Operator{T} - ops::Tuple{S,V} +struct KroneckerOperator{S, V, DS,RS,DI,RI,T} <: Operator{T} + ops::Tuple domainspace::DS rangespace::RS domaintensorizer::DI rangetensorizer::RI end +# function KroneckerOperator{S,V,DS,RS,DI,RI,T}(ops::NTuple{<:Any, <:Conversion}, ds, rs, di, ri) +# KroneckerOperator{DS,RS,DI,RI,T}(ops,ds,rs,di,ri) +# end + KroneckerOperator(A,B,ds::Space,rs::Space,di,ri) = KroneckerOperator{typeof(A),typeof(B),typeof(ds),typeof(rs),typeof(di),typeof(ri), promote_type(eltype(A),eltype(B))}((A,B),ds,rs,di,ri) +KroneckerOperator(A::Tuple,ds::Space,rs::Space,di,ri) = + KroneckerOperator{typeof(A[1]), typeof(A[2]), typeof(ds),typeof(rs),typeof(di),typeof(ri), + promote_type(eltype(A[1]),eltype(A[2]))}(A,ds,rs,di,ri) KroneckerOperator(A,B,ds::Space,rs::Space) = KroneckerOperator(A,B,ds,rs, CachedIterator(tensorizer(ds)),CachedIterator(tensorizer(rs))) +KroneckerOperator(A::Tuple, ds::Space, rs::Space) = KroneckerOperator(A,ds,rs, + CachedIterator(tensorizer(ds)),CachedIterator(tensorizer(rs))) + function KroneckerOperator(A,B) ds=domainspace(A)⊗domainspace(B) rs=rangespace(A)⊗rangespace(B) @@ -262,25 +272,18 @@ Base.transpose(S::ConstantTimesOperator) = sp.c*transpose(S.op) ### Calculus -#TODO: general dimension -function Derivative(S::TensorSpace{<:Any,<:EuclideanDomain{2}}, order) - @assert length(order)==2 - if order[1]==0 - Dy=Derivative(S.spaces[2],order[2]) - K=Operator(I,S.spaces[1])⊗Dy - elseif order[2]==0 - Dx=Derivative(S.spaces[1],order[1]) - K=Dx⊗Operator(I,S.spaces[2]) - else - Dx=Derivative(S.spaces[1],order[1]) - Dy=Derivative(S.spaces[2],order[2]) - K=Dx⊗Dy - end - DerivativeWrapper(K,order,S) +function Derivative(S::TensorSpace{<:Any,<:EuclideanDomain}, order) + @assert length(order)==length(S.spaces) + @inline Derivative_or_I(i) = order[i]>0 ? Derivative(S.spaces[i], order[i]) : Operator(I,S.spaces[i]) + DerivativeWrapper(mapreduce(i->Derivative_or_I(i),⊗,1:length(order)), order, S) end - DefiniteIntegral(S::TensorSpace) = DefiniteIntegralWrapper(mapreduce(DefiniteIntegral,⊗,S.spaces)) +function DefiniteIntegral(S::TensorSpace, dim::Vector{Int}) + @assert length(dim)==length(S.spaces) + @inline Int_or_I(i) = 1==dim[i] ? DefiniteIntegral(S.spaces[i]) : Operator(I,S.spaces[i]) + DefiniteIntegralWrapper(mapreduce(i->Int_or_I(i),⊗,1:length(dim))) +end @@ -401,6 +404,13 @@ function Conversion(a::TensorSpace2D,b::TensorSpace2D) ConversionWrapper(strictconvert(Operator{T}, K)) end +function Conversion(a::TensorSpace,b::TensorSpace) + C = map(Conversion,a.spaces,b.spaces) + K = KroneckerOperator(C, a, b) + T = promote_type(prectype(a),prectype(b)) + ConversionWrapper(strictconvert(Operator{T}, K)) +end + function Multiplication(f::Fun{<:TensorSpace}, S::TensorSpace) lr = LowRankFun(f) MAs = map(a->Multiplication(a,S.spaces[1]), lr.A) From 50bd00aa91a68de05aa341d29723e644f29a9ca7 Mon Sep 17 00:00:00 2001 From: Benjamin Zanger Date: Tue, 7 Mar 2023 23:07:57 +0100 Subject: [PATCH 2/9] updated blocklengths and new signature for Kroneckeroperator --- src/Multivariate/TensorSpace.jl | 12 ++++++++++-- src/PDE/KroneckerOperator.jl | 24 ++++++++++-------------- test/runtests.jl | 6 +++--- 3 files changed, 23 insertions(+), 19 deletions(-) diff --git a/src/Multivariate/TensorSpace.jl b/src/Multivariate/TensorSpace.jl index f4f2b9c5..868bbd24 100644 --- a/src/Multivariate/TensorSpace.jl +++ b/src/Multivariate/TensorSpace.jl @@ -215,7 +215,7 @@ blocklength(it,k::Block) = blocklength(it,k.n[1]) blocklength(it,k::BlockRange) = blocklength(it,Int.(k)) blocklengths(::TrivialTensorizer{2}) = 1:∞ - +blocklengths(::TrivialTensorizer{d}) where {d} = binomial.((1:∞).+(d-2), d-1) blocklengths(it::Tensorizer) = tensorblocklengths(it.blocks...) @@ -304,7 +304,15 @@ const TensorSpace2D{AA, BB, D,R} = TensorSpace{<:Tuple{AA, BB}, D, R} where {AA< const TensorSpaceND{d, D, R} = TensorSpace{<:NTuple{d, <:UnivariateSpace}, D, R} tensorizer(sp::TensorSpace) = Tensorizer(map(blocklengths,sp.spaces)) -blocklengths(S::TensorSpace) = tensorblocklengths(map(blocklengths,S.spaces)...) +function blocklengths(S::TensorSpace) + list_blocks = map(blocklengths,S.spaces) + if all(x->x == Ones{Int}(ℵ₀), list_blocks) + d = length(S.spaces) + return binomial.((1:∞).+(d-2), d-1) + else + return tensorblocklengths(list_blocks) + end +end # the evaluation is *, so the type will be the same as * diff --git a/src/PDE/KroneckerOperator.jl b/src/PDE/KroneckerOperator.jl index 75dbd75f..356894ee 100644 --- a/src/PDE/KroneckerOperator.jl +++ b/src/PDE/KroneckerOperator.jl @@ -6,24 +6,20 @@ export KroneckerOperator # KroneckerOperator gives the kronecker product of two 1D operators ######### -struct KroneckerOperator{S, V, DS,RS,DI,RI,T} <: Operator{T} - ops::Tuple +struct KroneckerOperator{TT, DS,RS,DI,RI,T} <: Operator{T} + ops::TT domainspace::DS rangespace::RS domaintensorizer::DI rangetensorizer::RI end -# function KroneckerOperator{S,V,DS,RS,DI,RI,T}(ops::NTuple{<:Any, <:Conversion}, ds, rs, di, ri) -# KroneckerOperator{DS,RS,DI,RI,T}(ops,ds,rs,di,ri) -# end - KroneckerOperator(A,B,ds::Space,rs::Space,di,ri) = - KroneckerOperator{typeof(A),typeof(B),typeof(ds),typeof(rs),typeof(di),typeof(ri), + KroneckerOperator{typeof((A, B)),typeof(ds),typeof(rs),typeof(di),typeof(ri), promote_type(eltype(A),eltype(B))}((A,B),ds,rs,di,ri) KroneckerOperator(A::Tuple,ds::Space,rs::Space,di,ri) = - KroneckerOperator{typeof(A[1]), typeof(A[2]), typeof(ds),typeof(rs),typeof(di),typeof(ri), + KroneckerOperator{typeof(A), typeof(ds),typeof(rs),typeof(di),typeof(ri), promote_type(eltype(A[1]),eltype(A[2]))}(A,ds,rs,di,ri) KroneckerOperator(A,B,ds::Space,rs::Space) = KroneckerOperator(A,B,ds,rs, @@ -67,7 +63,7 @@ function convert(::Type{Operator{T}},K::KroneckerOperator) where T<:Number K else ops = map(Operator{T}, K.ops) - KroneckerOperator{map(typeof, ops)..., + KroneckerOperator{typeof(ops), typeof(K.domainspace),typeof(K.rangespace), typeof(K.domaintensorizer),typeof(K.rangetensorizer),T}(ops, K.domainspace,K.rangespace, @@ -348,9 +344,9 @@ const Trivial2DTensorizer = CachedIterator{Tuple{Int,Int}, # This routine is an efficient version of KroneckerOperator for the case of # tensor product of trivial blocks -function BandedBlockBandedMatrix(S::SubOperator{T,KroneckerOperator{SS,V,DS,RS, +function BandedBlockBandedMatrix(S::SubOperator{T,KroneckerOperator{TT,DS,RS, Trivial2DTensorizer,Trivial2DTensorizer,T}, - Tuple{BlockRange1,BlockRange1}}) where {SS,V,DS,RS,T} + Tuple{BlockRange1,BlockRange1}}) where {TT,DS,RS,T} KR,JR = parentindices(S) KR_i, JR_i = Int.(KR), Int.(JR) @@ -382,9 +378,9 @@ function BandedBlockBandedMatrix(S::SubOperator{T,KroneckerOperator{SS,V,DS,RS, ret end -convert(::Type{BandedBlockBandedMatrix}, S::SubOperator{T,KroneckerOperator{SS,V,DS,RS, +convert(::Type{BandedBlockBandedMatrix}, S::SubOperator{T,KroneckerOperator{TT,DS,RS, Trivial2DTensorizer,Trivial2DTensorizer,T}, - Tuple{BlockRange1,BlockRange1}}) where {SS,V,DS,RS,T} = + Tuple{BlockRange1,BlockRange1}}) where {TT,DS,RS,T} = BandedBlockBandedMatrix(S) ## TensorSpace operators @@ -459,7 +455,7 @@ end # if the second operator is a constant, we may scale the first operator, # and apply it on the coefficients -function (*)(ko::KroneckerOperator{<:Operator, <:ConstantOperator}, pf::ProductFun) +function (*)(ko::KroneckerOperator{<:Tuple{<:Operator, <:ConstantOperator}}, pf::ProductFun) O1, O2 = ko.ops O12 = O2.λ * O1 ProductFun(map(x -> O12*x, pf.coefficients), factor(pf.space, 2)) diff --git a/test/runtests.jl b/test/runtests.jl index 56d6b999..f31947d9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -491,11 +491,11 @@ end @testset "kronecker type stability" begin T = ApproxFunBase.ToeplitzOperator(Float64[1,2], Float64[1,2]) TT = T ⊗ T - TypeExp = Union{ApproxFunBase.SubOperator{Float64, KroneckerOperator{ToeplitzOperator{Float64}, ToeplitzOperator{Float64}, TensorSpace{Tuple{SequenceSpace, SequenceSpace}, DomainSets.VcatDomain{2, Int64, (1, 1), Tuple{ApproxFunBase.PositiveIntegers, ApproxFunBase.PositiveIntegers}}, Union{}}, TensorSpace{Tuple{SequenceSpace, SequenceSpace}, DomainSets.VcatDomain{2, Int64, (1, 1), Tuple{ApproxFunBase.PositiveIntegers, ApproxFunBase.PositiveIntegers}}, Union{}}, ApproxFunBase.CachedIterator{Tuple{Int64, Int64}, ApproxFunBase.Tensorizer{Tuple{FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, ApproxFunBase.CachedIterator{Tuple{Int64, Int64}, ApproxFunBase.Tensorizer{Tuple{FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, Float64}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, Tuple{Int64, Int64}, Tuple{Infinities.InfiniteCardinal{0}, Infinities.InfiniteCardinal{0}}}, ApproxFunBase.SubOperator{Float64, KroneckerOperator{ToeplitzOperator{Float64}, ToeplitzOperator{Float64}, TensorSpace{Tuple{SequenceSpace, SequenceSpace}, DomainSets.VcatDomain{2, Int64, (1, 1), Tuple{ApproxFunBase.PositiveIntegers, ApproxFunBase.PositiveIntegers}}, Union{}}, TensorSpace{Tuple{SequenceSpace, SequenceSpace}, DomainSets.VcatDomain{2, Int64, (1, 1), Tuple{ApproxFunBase.PositiveIntegers, ApproxFunBase.PositiveIntegers}}, Union{}}, ApproxFunBase.CachedIterator{Tuple{Int64, Int64}, ApproxFunBase.Tensorizer{Tuple{FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, ApproxFunBase.CachedIterator{Tuple{Int64, Int64}, ApproxFunBase.Tensorizer{Tuple{FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, Float64}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, Tuple{Int64, Int64}, Tuple{Int64, Int64}}} + TypeExp = Union{ApproxFunBase.SubOperator{Float64, KroneckerOperator{Tuple{ToeplitzOperator{Float64}, ToeplitzOperator{Float64}}, TensorSpace{Tuple{SequenceSpace, SequenceSpace}, DomainSets.VcatDomain{2, Int64, (1, 1), Tuple{ApproxFunBase.PositiveIntegers, ApproxFunBase.PositiveIntegers}}, Union{}}, TensorSpace{Tuple{SequenceSpace, SequenceSpace}, DomainSets.VcatDomain{2, Int64, (1, 1), Tuple{ApproxFunBase.PositiveIntegers, ApproxFunBase.PositiveIntegers}}, Union{}}, ApproxFunBase.CachedIterator{Tuple{Int64, Int64}, ApproxFunBase.Tensorizer{Tuple{FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, ApproxFunBase.CachedIterator{Tuple{Int64, Int64}, ApproxFunBase.Tensorizer{Tuple{FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, Float64}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, Tuple{Int64, Int64}, Tuple{Infinities.InfiniteCardinal{0}, Infinities.InfiniteCardinal{0}}}, ApproxFunBase.SubOperator{Float64, KroneckerOperator{Tuple{ToeplitzOperator{Float64}, ToeplitzOperator{Float64}}, TensorSpace{Tuple{SequenceSpace, SequenceSpace}, DomainSets.VcatDomain{2, Int64, (1, 1), Tuple{ApproxFunBase.PositiveIntegers, ApproxFunBase.PositiveIntegers}}, Union{}}, TensorSpace{Tuple{SequenceSpace, SequenceSpace}, DomainSets.VcatDomain{2, Int64, (1, 1), Tuple{ApproxFunBase.PositiveIntegers, ApproxFunBase.PositiveIntegers}}, Union{}}, ApproxFunBase.CachedIterator{Tuple{Int64, Int64}, ApproxFunBase.Tensorizer{Tuple{FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, ApproxFunBase.CachedIterator{Tuple{Int64, Int64}, ApproxFunBase.Tensorizer{Tuple{FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, Float64}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, Tuple{Int64, Int64}, Tuple{Int64, Int64}}} @inferred TypeExp view(T, 1:1, 1:1) - TypeExp2 = Union{ApproxFunBase.SubOperator{Float64, KroneckerOperator{ToeplitzOperator{Float64}, ToeplitzOperator{Float64}, TensorSpace{Tuple{SequenceSpace, SequenceSpace}, DomainSets.VcatDomain{2, Int64, (1, 1), Tuple{ApproxFunBase.PositiveIntegers, ApproxFunBase.PositiveIntegers}}, Union{}}, TensorSpace{Tuple{SequenceSpace, SequenceSpace}, DomainSets.VcatDomain{2, Int64, (1, 1), Tuple{ApproxFunBase.PositiveIntegers, ApproxFunBase.PositiveIntegers}}, Union{}}, ApproxFunBase.CachedIterator{Tuple{Int64, Int64}, ApproxFunBase.Tensorizer{Tuple{FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, ApproxFunBase.CachedIterator{Tuple{Int64, Int64}, ApproxFunBase.Tensorizer{Tuple{FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, Float64}, Tuple{StepRange{Int64, Int64}, StepRange{Int64, Int64}}, Tuple{Int64, Int64}, Tuple{Infinities.InfiniteCardinal{0}, Infinities.InfiniteCardinal{0}}}, ApproxFunBase.SubOperator{Float64, KroneckerOperator{ToeplitzOperator{Float64}, ToeplitzOperator{Float64}, TensorSpace{Tuple{SequenceSpace, SequenceSpace}, DomainSets.VcatDomain{2, Int64, (1, 1), Tuple{ApproxFunBase.PositiveIntegers, ApproxFunBase.PositiveIntegers}}, Union{}}, TensorSpace{Tuple{SequenceSpace, SequenceSpace}, DomainSets.VcatDomain{2, Int64, (1, 1), Tuple{ApproxFunBase.PositiveIntegers, ApproxFunBase.PositiveIntegers}}, Union{}}, ApproxFunBase.CachedIterator{Tuple{Int64, Int64}, ApproxFunBase.Tensorizer{Tuple{FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, ApproxFunBase.CachedIterator{Tuple{Int64, Int64}, ApproxFunBase.Tensorizer{Tuple{FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, Float64}, Tuple{StepRange{Int64, Int64}, StepRange{Int64, Int64}}, Tuple{Int64, Int64}, Tuple{Int64, Int64}}} + TypeExp2 = Union{ApproxFunBase.SubOperator{Float64, KroneckerOperator{Tuple{ToeplitzOperator{Float64}, ToeplitzOperator{Float64}}, TensorSpace{Tuple{SequenceSpace, SequenceSpace}, DomainSets.VcatDomain{2, Int64, (1, 1), Tuple{ApproxFunBase.PositiveIntegers, ApproxFunBase.PositiveIntegers}}, Union{}}, TensorSpace{Tuple{SequenceSpace, SequenceSpace}, DomainSets.VcatDomain{2, Int64, (1, 1), Tuple{ApproxFunBase.PositiveIntegers, ApproxFunBase.PositiveIntegers}}, Union{}}, ApproxFunBase.CachedIterator{Tuple{Int64, Int64}, ApproxFunBase.Tensorizer{Tuple{FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, ApproxFunBase.CachedIterator{Tuple{Int64, Int64}, ApproxFunBase.Tensorizer{Tuple{FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, Float64}, Tuple{StepRange{Int64, Int64}, StepRange{Int64, Int64}}, Tuple{Int64, Int64}, Tuple{Infinities.InfiniteCardinal{0}, Infinities.InfiniteCardinal{0}}}, ApproxFunBase.SubOperator{Float64, KroneckerOperator{Tuple{ToeplitzOperator{Float64}, ToeplitzOperator{Float64}}, TensorSpace{Tuple{SequenceSpace, SequenceSpace}, DomainSets.VcatDomain{2, Int64, (1, 1), Tuple{ApproxFunBase.PositiveIntegers, ApproxFunBase.PositiveIntegers}}, Union{}}, TensorSpace{Tuple{SequenceSpace, SequenceSpace}, DomainSets.VcatDomain{2, Int64, (1, 1), Tuple{ApproxFunBase.PositiveIntegers, ApproxFunBase.PositiveIntegers}}, Union{}}, ApproxFunBase.CachedIterator{Tuple{Int64, Int64}, ApproxFunBase.Tensorizer{Tuple{FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, ApproxFunBase.CachedIterator{Tuple{Int64, Int64}, ApproxFunBase.Tensorizer{Tuple{FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, Float64}, Tuple{StepRange{Int64, Int64}, StepRange{Int64, Int64}}, Tuple{Int64, Int64}, Tuple{Int64, Int64}}} @inferred TypeExp2 view(T, 1:1:1, 1:1:1) - TypeExp3 = Union{ApproxFunBase.SubOperator{Float64, KroneckerOperator{ToeplitzOperator{Float64}, ToeplitzOperator{Float64}, TensorSpace{Tuple{SequenceSpace, SequenceSpace}, DomainSets.VcatDomain{2, Int64, (1, 1), Tuple{ApproxFunBase.PositiveIntegers, ApproxFunBase.PositiveIntegers}}, Union{}}, TensorSpace{Tuple{SequenceSpace, SequenceSpace}, DomainSets.VcatDomain{2, Int64, (1, 1), Tuple{ApproxFunBase.PositiveIntegers, ApproxFunBase.PositiveIntegers}}, Union{}}, ApproxFunBase.CachedIterator{Tuple{Int64, Int64}, ApproxFunBase.Tensorizer{Tuple{Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, ApproxFunBase.CachedIterator{Tuple{Int64, Int64}, ApproxFunBase.Tensorizer{Tuple{Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, Float64}, Tuple{InfiniteArrays.InfUnitRange{Int64}, InfiniteArrays.InfUnitRange{Int64}}, Tuple{Infinities.Infinity, Infinities.Infinity}, Tuple{Infinities.InfiniteCardinal{0}, Infinities.InfiniteCardinal{0}}}, ApproxFunBase.SubOperator{Float64, KroneckerOperator{ToeplitzOperator{Float64}, ToeplitzOperator{Float64}, TensorSpace{Tuple{SequenceSpace, SequenceSpace}, DomainSets.VcatDomain{2, Int64, (1, 1), Tuple{ApproxFunBase.PositiveIntegers, ApproxFunBase.PositiveIntegers}}, Union{}}, TensorSpace{Tuple{SequenceSpace, SequenceSpace}, DomainSets.VcatDomain{2, Int64, (1, 1), Tuple{ApproxFunBase.PositiveIntegers, ApproxFunBase.PositiveIntegers}}, Union{}}, ApproxFunBase.CachedIterator{Tuple{Int64, Int64}, ApproxFunBase.Tensorizer{Tuple{Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, ApproxFunBase.CachedIterator{Tuple{Int64, Int64}, ApproxFunBase.Tensorizer{Tuple{Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, Float64}, Tuple{InfiniteArrays.InfUnitRange{Int64}, InfiniteArrays.InfUnitRange{Int64}}, Tuple{Infinities.Infinity, Infinities.Infinity}, Tuple{Int64, Int64}}} + TypeExp3 = Union{ApproxFunBase.SubOperator{Float64, KroneckerOperator{Tuple{ToeplitzOperator{Float64}, ToeplitzOperator{Float64}}, TensorSpace{Tuple{SequenceSpace, SequenceSpace}, DomainSets.VcatDomain{2, Int64, (1, 1), Tuple{ApproxFunBase.PositiveIntegers, ApproxFunBase.PositiveIntegers}}, Union{}}, TensorSpace{Tuple{SequenceSpace, SequenceSpace}, DomainSets.VcatDomain{2, Int64, (1, 1), Tuple{ApproxFunBase.PositiveIntegers, ApproxFunBase.PositiveIntegers}}, Union{}}, ApproxFunBase.CachedIterator{Tuple{Int64, Int64}, ApproxFunBase.Tensorizer{Tuple{Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, ApproxFunBase.CachedIterator{Tuple{Int64, Int64}, ApproxFunBase.Tensorizer{Tuple{Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, Float64}, Tuple{InfiniteArrays.InfUnitRange{Int64}, InfiniteArrays.InfUnitRange{Int64}}, Tuple{Infinities.Infinity, Infinities.Infinity}, Tuple{Infinities.InfiniteCardinal{0}, Infinities.InfiniteCardinal{0}}}, ApproxFunBase.SubOperator{Float64, KroneckerOperator{Tuple{ToeplitzOperator{Float64}, ToeplitzOperator{Float64}}, TensorSpace{Tuple{SequenceSpace, SequenceSpace}, DomainSets.VcatDomain{2, Int64, (1, 1), Tuple{ApproxFunBase.PositiveIntegers, ApproxFunBase.PositiveIntegers}}, Union{}}, TensorSpace{Tuple{SequenceSpace, SequenceSpace}, DomainSets.VcatDomain{2, Int64, (1, 1), Tuple{ApproxFunBase.PositiveIntegers, ApproxFunBase.PositiveIntegers}}, Union{}}, ApproxFunBase.CachedIterator{Tuple{Int64, Int64}, ApproxFunBase.Tensorizer{Tuple{Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, ApproxFunBase.CachedIterator{Tuple{Int64, Int64}, ApproxFunBase.Tensorizer{Tuple{Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, Float64}, Tuple{InfiniteArrays.InfUnitRange{Int64}, InfiniteArrays.InfUnitRange{Int64}}, Tuple{Infinities.Infinity, Infinities.Infinity}, Tuple{Int64, Int64}}} @inferred TypeExp3 view(T, 1:∞, 1:∞) end end From 1398dd0f07889ae37174d642c60c0ad4908c6eb4 Mon Sep 17 00:00:00 2001 From: Benjamin Zanger Date: Wed, 8 Mar 2023 10:45:02 +0100 Subject: [PATCH 3/9] rework blocklengths for trivialtensor spaces --- src/Caching/blockbanded.jl | 6 ++++++ src/Multivariate/TensorSpace.jl | 9 +++++++-- 2 files changed, 13 insertions(+), 2 deletions(-) diff --git a/src/Caching/blockbanded.jl b/src/Caching/blockbanded.jl index dcbb8120..d923bf07 100644 --- a/src/Caching/blockbanded.jl +++ b/src/Caching/blockbanded.jl @@ -22,6 +22,12 @@ end diagblockshift(a,b) = error("Developer: Not implemented for blocklengths $a, $b") +function diagblockshift(a::BroadcastArray{<:Int, 1}, b::BroadcastArray{<:Int, 1}) + if a.f == b.f + return 0 + end + error("Broadcastvectors of blocklengths $a, $b are not implemented.") +end function diagblockshift(a::AbstractRange, b::AbstractRange) @assert step(a) == step(b) first(b)-first(a) diff --git a/src/Multivariate/TensorSpace.jl b/src/Multivariate/TensorSpace.jl index 868bbd24..901998b8 100644 --- a/src/Multivariate/TensorSpace.jl +++ b/src/Multivariate/TensorSpace.jl @@ -215,7 +215,12 @@ blocklength(it,k::Block) = blocklength(it,k.n[1]) blocklength(it,k::BlockRange) = blocklength(it,Int.(k)) blocklengths(::TrivialTensorizer{2}) = 1:∞ -blocklengths(::TrivialTensorizer{d}) where {d} = binomial.((1:∞).+(d-2), d-1) + +## anonymous function needed in order to compare if two blocklenghts are equal +_blocklengths_trivialTensorizer(d) = let d=d + x->binomial(x+(d-2), d-1) +end +blocklengths(::TrivialTensorizer{d}) where {d} = _blocklengths_trivialTensorizer(d).(1:∞) blocklengths(it::Tensorizer) = tensorblocklengths(it.blocks...) @@ -308,7 +313,7 @@ function blocklengths(S::TensorSpace) list_blocks = map(blocklengths,S.spaces) if all(x->x == Ones{Int}(ℵ₀), list_blocks) d = length(S.spaces) - return binomial.((1:∞).+(d-2), d-1) + return _blocklengths_trivialTensorizer(d).(1:∞) else return tensorblocklengths(list_blocks) end From 194e1996017a83cc301efce792e75ef05c9934e5 Mon Sep 17 00:00:00 2001 From: Benjamin Zanger Date: Wed, 8 Mar 2023 18:22:34 +0100 Subject: [PATCH 4/9] the right bandwidths in Kroneckeroperator are still missing --- src/Multivariate/TensorSpace.jl | 47 +++++++++++++++++++++++-- src/Multivariate/TrivialTensorFun.jl | 10 +++--- src/PDE/KroneckerOperator.jl | 52 ++++++++++++++++++---------- 3 files changed, 82 insertions(+), 27 deletions(-) diff --git a/src/Multivariate/TensorSpace.jl b/src/Multivariate/TensorSpace.jl index 901998b8..a2f434c8 100644 --- a/src/Multivariate/TensorSpace.jl +++ b/src/Multivariate/TensorSpace.jl @@ -31,6 +31,7 @@ end const InfOnes = Ones{Int,1,Tuple{OneToInf{Int}}} const Tensorizer2D{AA, BB} = Tensorizer{Tuple{AA, BB}} const TrivialTensorizer{d} = Tensorizer{NTuple{d,InfOnes}} +const TrivialConstantTensorizer{d} = Tensorizer{<:NTuple{d,SVector{1, <:Int}}} # for all dimensions constant eltype(::Type{<:Tensorizer{<:Tuple{Vararg{Any,N}}}}) where {N} = NTuple{N,Int} dimensions(a::Tensorizer) = map(sum,a.blocks) @@ -40,6 +41,24 @@ Base.IteratorSize(::Type{Tensorizer{T}}) where {T<:Tuple} = _IteratorSize(T) Base.keys(a::Tensorizer) = oneto(length(a)) + +function start(a::TrivialConstantTensorizer{d}) where {d} + @assert length(a) == 1 + block = ntuple(one, d) + return (block, (0,1)) +end + +function next(a::TrivialConstantTensorizer{d}, iterator_tuple) where {d} + (block, (i,tot)) = iterator_tuple + ret = block + ret, (block, (i+1,tot)) +end + +function done(a::TrivialConstantTensorizer, iterator_tuple)::Bool + i, tot = last(iterator_tuple) + return i ≥ tot +end + function start(a::TrivialTensorizer{d}) where {d} # ((block_dim_1, block_dim_2,...), (itaration_number, iterator, iterator_state)), (itemssofar, length) block = ntuple(one, d) @@ -315,7 +334,7 @@ function blocklengths(S::TensorSpace) d = length(S.spaces) return _blocklengths_trivialTensorizer(d).(1:∞) else - return tensorblocklengths(list_blocks) + return tensorblocklengths(list_blocks...) end end @@ -634,7 +653,7 @@ function totensor(it::Tensorizer,M::AbstractVector) ret end -@inline function totensoriterator(it::TrivialTensorizer{d},M::AbstractVector) where {d} +@inline function totensoriterator(it::Union{TrivialTensorizer{d}, TrivialConstantTensorizer{d}} ,M::AbstractVector) where {d} B=block(it,length(M)) return it, M, B end @@ -677,7 +696,29 @@ evaluate(f::AbstractVector,S::TensorSpace2D,x) = ProductFun(totensor(S,f), S)(x. evaluate(f::AbstractVector,S::TensorSpace2D,x,y) = ProductFun(totensor(S,f),S)(x,y) # ND evaluation functions of Trivial Spaces -evaluate(f::AbstractVector,S::TensorSpaceND,x) = TrivialTensorFun(totensor(S, f)..., S)(x...) +not_const_spaces_indices(S) = filter!(i->i≠0, map(i->S.spaces[i] isa ConstantSpace ? 0 : i,1:length(S.spaces))) +function evaluate(f::AbstractVector,S::TensorSpaceND,x) + if !any(s->s isa ConstantSpace, S.spaces) + return TrivialTensorFun(totensor(S, f)..., S)(x...) + end + not_cons_indices = not_const_spaces_indices(S) + xmod = if length(x) == length(not_cons_indices) + x + else + x[not_cons_indices] + end + (length(not_cons_indices) == 0) && return f[1] + S_new = reduce(⊗, S.spaces[not_cons_indices]) + if length(not_cons_indices) > 2 + return TrivialTensorFun(totensor(S_new, f)..., S_new)(x...) + elseif length(S_new) == 2 + return ProductFun(totensor(S_new, f), S_new)(x...) + elseif length(S_new) == 1 + return Fun(S_new[1], f)(x) + else + error("This should not happen") + end +end coefficientmatrix(f::Fun{<:AbstractProductSpace}) = totensor(space(f),f.coefficients) diff --git a/src/Multivariate/TrivialTensorFun.jl b/src/Multivariate/TrivialTensorFun.jl index 61363a8e..6790e5c7 100644 --- a/src/Multivariate/TrivialTensorFun.jl +++ b/src/Multivariate/TrivialTensorFun.jl @@ -3,16 +3,14 @@ struct TrivialTensorFun{d, SS<:TensorSpaceND{d}, T<:Number} <: MultivariateFun{T, d} space::SS - coefficients::Vector{T} + coefficients::AbstractVector{T} iterator::TrivialTensorizer{d} orders::Block{1, Int} end -function TrivialTensorFun(iter::TrivialTensorizer{d},cfs::Vector{T},blk::Block, sp::TensorSpaceND{d}) where {T<:Number,d} - if any(map(dimension, sp.spaces).!=ℵ₀) - error("This Space is not a Trivial Tensor space!") - end +function TrivialTensorFun(iter::TrivialTensorizer{d}, + cfs::AbstractVector{T}, blk::Block, sp::TensorSpaceND{d}) where {T<:Number,d} TrivialTensorFun(sp, cfs, iter, blk) end @@ -20,7 +18,7 @@ end # TensorSpace evaluation function evaluate(f::TrivialTensorFun{d, SS, T},x...) where {d, SS, T} - highest_order = f.orders.n[1] + highest_order = f.orders.n[1]-1 n = length(f.coefficients) # this could be lazy evaluated for the sparse case diff --git a/src/PDE/KroneckerOperator.jl b/src/PDE/KroneckerOperator.jl index 356894ee..e6a6c6c7 100644 --- a/src/PDE/KroneckerOperator.jl +++ b/src/PDE/KroneckerOperator.jl @@ -27,6 +27,22 @@ KroneckerOperator(A,B,ds::Space,rs::Space) = KroneckerOperator(A,B,ds,rs, KroneckerOperator(A::Tuple, ds::Space, rs::Space) = KroneckerOperator(A,ds,rs, CachedIterator(tensorizer(ds)),CachedIterator(tensorizer(rs))) +# DON'T nest KroneckerOperators +function KroneckerOperator(A::KroneckerOperator,B::KroneckerOperator) + ds=domainspace(A)⊗domainspace(B) + rs=rangespace(A)⊗rangespace(B) + KroneckerOperator((A.ops..., B.ops...), ds, rs) +end +function KroneckerOperator(A::KroneckerOperator,B) + ds=domainspace(A)⊗domainspace(B) + rs=rangespace(A)⊗rangespace(B) + KroneckerOperator((A.ops..., B), ds, rs) +end +function KroneckerOperator(A,B::KroneckerOperator) + ds=domainspace(A)⊗domainspace(B) + rs=rangespace(A)⊗rangespace(B) + KroneckerOperator((A, B.ops...), ds, rs) +end function KroneckerOperator(A,B) ds=domainspace(A)⊗domainspace(B) rs=rangespace(A)⊗rangespace(B) @@ -46,15 +62,13 @@ KroneckerOperator(K::KroneckerOperator) = K KroneckerOperator(T::TimesOperator) = mapfoldr(op -> KroneckerOperator(op), *, T.ops) function promotedomainspace(K::KroneckerOperator,ds::TensorSpace) - A=promotedomainspace(K.ops[1],ds.spaces[1]) - B=promotedomainspace(K.ops[2],ds.spaces[2]) - KroneckerOperator(A,B,ds,rangespace(A)⊗rangespace(B)) + A = (i->promoterangespace(K.ops[i], ds.spaces[i]) for i=1:length(K.ops)) + KroneckerOperator(A,ds,rangespace(K)) ## TODO: maybe rangespace(K) has to be replaces by rangespace(A[1])⊗... end function promoterangespace(K::KroneckerOperator,rs::TensorSpace) - A=promoterangespace(K.ops[1],rs.spaces[1]) - B=promoterangespace(K.ops[2],rs.spaces[2]) - KroneckerOperator(A,B,domainspace(K),rs) + A = (i->promoterangespace(K.ops[i], rs.spaces[i]) for i=1:length(K.ops)) + KroneckerOperator(A,domainspace(K),rs) end @@ -130,11 +144,12 @@ end _isbandedblockbanded(::Tuple{}) = true blockbandwidths(K::KroneckerOperator) = - (blockbandwidth(K.ops[1],1)+blockbandwidth(K.ops[2],1), - blockbandwidth(K.ops[1],2)+blockbandwidth(K.ops[2],2)) + (mapreduce(k->blockbandwidth(k,1), +, K.ops), + mapreduce(k->blockbandwidth(k,2), +, K.ops)) # If each block were in turn BlockBandedMatrix, these would # be the bandwidths +# TODO: How does this work for multiple Ops? subblock_blockbandwidths(K::KroneckerOperator) = (max(blockbandwidth(K.ops[1],1),blockbandwidth(K.ops[2],2)) , max(blockbandwidth(K.ops[1],2),blockbandwidth(K.ops[2],1))) @@ -188,20 +203,24 @@ domaintensorizer(K::KroneckerOperator) = K.domaintensorizer rangetensorizer(K::KroneckerOperator) = K.rangetensorizer +# For 2 ops, we can do this # we suport 4-indexing with KroneckerOperator # If A is K x J and B is N x M, then w # index to match KO=reshape(kron(B,A),N,K,M,J) # that is # KO[k,n,j,m] = A[k,j]*B[n,m] + # TODO: arbitrary number of ops +# We get tuples of arbitraty length from the tensorizers -getindex(KO::KroneckerOperator,k::Integer,n::Integer,j::Integer,m::Integer) = - KO.ops[1][k,j]*KO.ops[2][n,m] +## this should not be used anymore, can be deleted +# getindex(KO::KroneckerOperator,k::Integer,n::Integer,j::Integer,m::Integer) = +# KO.ops[1][k,j]*KO.ops[2][n,m] function getindex(KO::KroneckerOperator,kin::Integer,jin::Integer) - j,m=KO.domaintensorizer[jin] - k,n=KO.rangetensorizer[kin] - KO[k,n,j,m] + domain_tuple=KO.domaintensorizer[jin] + range_tuple=KO.rangetensorizer[kin] + mapreduce((k,i_in,j_out)->k[i_in,j_out], *, KO.ops, range_tuple, domain_tuple) end function getindex(KO::KroneckerOperator,k::Integer) @@ -218,11 +237,8 @@ end function *(A::KroneckerOperator, B::KroneckerOperator) dspB = domainspace(B) rspA = rangespace(A) - A1, A2 = A.ops - B1, B2 = B.ops - AB1 = A1 * B1 - AB2 = A2 * B2 - KroneckerOperator(AB1, AB2, dspB, rspA) + AB = Tuple([a*b for (a,b) in zip(A.ops, B.ops)]) + KroneckerOperator(AB, dspB, rspA) end From 97a2f8249b5aa20c05a9b32cfb3583d2d830d92f Mon Sep 17 00:00:00 2001 From: Benjamin Zanger Date: Wed, 8 Mar 2023 18:44:11 +0100 Subject: [PATCH 5/9] drop blockbanded for Kroneckeroperators for more than 2D2 --- src/PDE/KroneckerOperator.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/PDE/KroneckerOperator.jl b/src/PDE/KroneckerOperator.jl index e6a6c6c7..517e9fce 100644 --- a/src/PDE/KroneckerOperator.jl +++ b/src/PDE/KroneckerOperator.jl @@ -136,7 +136,7 @@ for f in [:isblockbanded, :israggedbelow, :isdiag] $(_f)(::Tuple{}) = true end end -isbandedblockbanded(K::KroneckerOperator) = _isbandedblockbanded(K.ops) +isbandedblockbanded(K::KroneckerOperator) = length(K.ops)>2 ? false : _isbandedblockbanded(K.ops) isbandedblockbandedcheck(op) = isbanded(op) && isinf(size(op,1)) && isinf(size(op,2)) function _isbandedblockbanded(ops::Tuple) isbandedblockbandedcheck(first(ops)) && _isbandedblockbanded(Base.tail(ops)) From 0ec33a10b6982b5b559e984332860aeaa630b856 Mon Sep 17 00:00:00 2001 From: Benjamin Zanger Date: Thu, 9 Mar 2023 11:27:49 +0100 Subject: [PATCH 6/9] iterator for mixed constant trivial spaces --- src/Multivariate/TensorSpace.jl | 73 +++++++++++++++++++++++++++++---- 1 file changed, 66 insertions(+), 7 deletions(-) diff --git a/src/Multivariate/TensorSpace.jl b/src/Multivariate/TensorSpace.jl index a2f434c8..ad4a6860 100644 --- a/src/Multivariate/TensorSpace.jl +++ b/src/Multivariate/TensorSpace.jl @@ -30,8 +30,10 @@ end const InfOnes = Ones{Int,1,Tuple{OneToInf{Int}}} const Tensorizer2D{AA, BB} = Tensorizer{Tuple{AA, BB}} +const MixedTrivConstTensorizer{d} = Tensorizer{<:Tuple{Vararg{Union{InfOnes, SVector{1, <:Int}},d}}} # const or trivial +# TrivialTensorizer and ConstantTensorizer are special cases of MixedTrivConstTensorizer const TrivialTensorizer{d} = Tensorizer{NTuple{d,InfOnes}} -const TrivialConstantTensorizer{d} = Tensorizer{<:NTuple{d,SVector{1, <:Int}}} # for all dimensions constant +const ConstantTensorizer{d} = Tensorizer{<:NTuple{d,SVector{1, <:Int}}} # for all dimensions constant eltype(::Type{<:Tensorizer{<:Tuple{Vararg{Any,N}}}}) where {N} = NTuple{N,Int} dimensions(a::Tensorizer) = map(sum,a.blocks) @@ -42,19 +44,19 @@ Base.IteratorSize(::Type{Tensorizer{T}}) where {T<:Tuple} = _IteratorSize(T) Base.keys(a::Tensorizer) = oneto(length(a)) -function start(a::TrivialConstantTensorizer{d}) where {d} +function start(a::ConstantTensorizer{d}) where {d} @assert length(a) == 1 block = ntuple(one, d) return (block, (0,1)) end -function next(a::TrivialConstantTensorizer{d}, iterator_tuple) where {d} +function next(a::ConstantTensorizer{d}, iterator_tuple) where {d} (block, (i,tot)) = iterator_tuple ret = block ret, (block, (i+1,tot)) end -function done(a::TrivialConstantTensorizer, iterator_tuple)::Bool +function done(a::ConstantTensorizer, iterator_tuple)::Bool i, tot = last(iterator_tuple) return i ≥ tot end @@ -97,21 +99,71 @@ function next(a::TrivialTensorizer{d}, iterator_tuple) where {d} ret, ((block, (j, iterator, iter_state)), (i,tot)) end - function done(a::TrivialTensorizer, iterator_tuple)::Bool i, tot = last(iterator_tuple) return i ≥ tot end +function start(a::MixedTrivConstTensorizer{d}) where {d} + # const indices are always left to be one + relevant_ind = filter!(i->i≠0, map(i->a.blocks[i] isa SVector{1, <:Int} ? 0 : i,1:length(a.blocks))) + real_d = length(relevant_ind) + # ((block_dim_1, block_dim_2,...), (itaration_number, iterator, iterator_state)), (itemssofar, length) + block = ones(Int, real_d) + return (block, (relevant_ind, real_d, 0, nothing, nothing)), (0,length(a)) +end + +function next(a::MixedTrivConstTensorizer{d}, iterator_tuple) where {d} + (block, (relevant_ind, real_d, j, iterator, iter_state)), (i,tot) = iterator_tuple + + @inline function check_block_finished(j, iterator, block) + if iterator === nothing + return true + end + # there are N-1 over d-1 combinations in a block + amount_combinations_block = binomial(sum(block)-1, real_d-1) + # check if all combinations have been iterated over + amount_combinations_block <= j + end + + ret_vec = ones(Int, d) + ret_vec[relevant_ind] = reverse(block) + ret = Tuple(SVector{d}(ret_vec)) + + if check_block_finished(j, iterator, block) # end of new block + # set up iterator for new block + current_sum = sum(block) + iterator = multiexponents(real_d, current_sum+1-real_d) + iter_state = nothing + j = 0 + end + + # increase block, or initialize new block + _res, iter_state = iterate(iterator, iter_state) + # res = Tuple(SVector{real_d}(_res)) + block = _res.+1 + j = j+1 + + ret, ((block, (relevant_ind, real_d, j, iterator, iter_state)), (i,tot)) +end + +function done(a::MixedTrivConstTensorizer, iterator_tuple)::Bool + i, tot = last(iterator_tuple) + return i ≥ tot +end + + # (blockrow,blockcol), (subrow,subcol), (rowshift,colshift), (numblockrows,numblockcols), (itemssofar, length) start(a::Tensorizer2D) = _start(a) start(a::TrivialTensorizer{2}) = _start(a) +start(a::MixedTrivConstTensorizer{2}) = _start(a) _start(a) = (1,1, 1,1, 0,0, a.blocks[1][1],a.blocks[2][1]), (0,length(a)) next(a::Tensorizer2D, state) = _next(a, state::typeof(_start(a))) next(a::TrivialTensorizer{2}, state) = _next(a, state::typeof(_start(a))) +next(a::MixedTrivConstTensorizer{2}, state) = _next(a, state::typeof(_start(a))) function _next(a, st) (K,J, k,j, rsh,csh, n,m), (i,tot) = st @@ -140,6 +192,7 @@ end done(a::Tensorizer2D, state) = _done(a, state::typeof(_start(a))) done(a::TrivialTensorizer{2}, state) = _done(a, state::typeof(_start(a))) +done(a::MixedTrivConstTensorizer{2}, state) = _done(a, state::typeof(_start(a))) function _done(a, st)::Bool i, tot = last(st) @@ -240,7 +293,13 @@ _blocklengths_trivialTensorizer(d) = let d=d x->binomial(x+(d-2), d-1) end blocklengths(::TrivialTensorizer{d}) where {d} = _blocklengths_trivialTensorizer(d).(1:∞) - +blocklengths(::ConstantTensorizer) = SVector(1) +function blocklengths(a::MixedTrivConstTensorizer{d}) where {d} + real_d = mapreduce(bl->bl isa SVector{1, <:Int} ? 0 : 1, +, a.blocks) + real_d == 0 && return SVector(1) + real_d == 1 && return 1:∞ + return _blocklengths_trivialTensorizer(real_d).(1:∞) +end blocklengths(it::Tensorizer) = tensorblocklengths(it.blocks...) blocklengths(it::CachedIterator) = blocklengths(it.iterator) @@ -653,7 +712,7 @@ function totensor(it::Tensorizer,M::AbstractVector) ret end -@inline function totensoriterator(it::Union{TrivialTensorizer{d}, TrivialConstantTensorizer{d}} ,M::AbstractVector) where {d} +@inline function totensoriterator(it::MixedTrivConstTensorizer{d} ,M::AbstractVector) where {d} B=block(it,length(M)) return it, M, B end From 8869fec7421bbbfc3df00dabfea5a7637c034cad Mon Sep 17 00:00:00 2001 From: Benjamin Zanger Date: Thu, 9 Mar 2023 15:04:56 +0100 Subject: [PATCH 7/9] add Integral for multivariate case --- src/PDE/KroneckerOperator.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/PDE/KroneckerOperator.jl b/src/PDE/KroneckerOperator.jl index 517e9fce..b0237c33 100644 --- a/src/PDE/KroneckerOperator.jl +++ b/src/PDE/KroneckerOperator.jl @@ -290,6 +290,13 @@ function Derivative(S::TensorSpace{<:Any,<:EuclideanDomain}, order) DerivativeWrapper(mapreduce(i->Derivative_or_I(i),⊗,1:length(order)), order, S) end +function Integral(S::TensorSpace{<:Any,<:EuclideanDomain}, order) + @assert length(order)==length(S.spaces) + @assert max(order...)<=1 + @inline Integral_or_I(i) = order[i]>0 ? Integral(S.spaces[i]) : Operator(I,S.spaces[i]) + IntegralWrapper(mapreduce(i->Integral_or_I(i),⊗,1:length(order))) +end + DefiniteIntegral(S::TensorSpace) = DefiniteIntegralWrapper(mapreduce(DefiniteIntegral,⊗,S.spaces)) function DefiniteIntegral(S::TensorSpace, dim::Vector{Int}) @assert length(dim)==length(S.spaces) From 97c5c2808be1f4a62dddd97930f5c30557ff019f Mon Sep 17 00:00:00 2001 From: Benjamin Zanger Date: Mon, 20 Mar 2023 13:16:30 +0100 Subject: [PATCH 8/9] fix in promotedomainspace --- src/PDE/KroneckerOperator.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/PDE/KroneckerOperator.jl b/src/PDE/KroneckerOperator.jl index b0237c33..f0c3d853 100644 --- a/src/PDE/KroneckerOperator.jl +++ b/src/PDE/KroneckerOperator.jl @@ -63,7 +63,8 @@ KroneckerOperator(T::TimesOperator) = mapfoldr(op -> KroneckerOperator(op), *, T function promotedomainspace(K::KroneckerOperator,ds::TensorSpace) A = (i->promoterangespace(K.ops[i], ds.spaces[i]) for i=1:length(K.ops)) - KroneckerOperator(A,ds,rangespace(K)) ## TODO: maybe rangespace(K) has to be replaces by rangespace(A[1])⊗... + range = reduce(⊗, [rangespace(a) for a∈A]) + KroneckerOperator(A,ds,range) end function promoterangespace(K::KroneckerOperator,rs::TensorSpace) From 40fbf7ae77907f0fccfbe678d5857fea3ca798ec Mon Sep 17 00:00:00 2001 From: Benjamin Zanger Date: Mon, 20 Mar 2023 14:16:21 +0100 Subject: [PATCH 9/9] fix tuple of promote Kronercker spaces --- src/PDE/KroneckerOperator.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/PDE/KroneckerOperator.jl b/src/PDE/KroneckerOperator.jl index f0c3d853..4adc5ff5 100644 --- a/src/PDE/KroneckerOperator.jl +++ b/src/PDE/KroneckerOperator.jl @@ -62,13 +62,13 @@ KroneckerOperator(K::KroneckerOperator) = K KroneckerOperator(T::TimesOperator) = mapfoldr(op -> KroneckerOperator(op), *, T.ops) function promotedomainspace(K::KroneckerOperator,ds::TensorSpace) - A = (i->promoterangespace(K.ops[i], ds.spaces[i]) for i=1:length(K.ops)) - range = reduce(⊗, [rangespace(a) for a∈A]) - KroneckerOperator(A,ds,range) + A = ntuple(i->promotedomainspace(K.ops[i], ds.spaces[i]), length(K.ops)) + rs = reduce(⊗, [rangespace(a) for a∈A]) + KroneckerOperator(A,ds,rs) end function promoterangespace(K::KroneckerOperator,rs::TensorSpace) - A = (i->promoterangespace(K.ops[i], rs.spaces[i]) for i=1:length(K.ops)) + A = ntuple(i->promoterangespace(K.ops[i], rs.spaces[i]), length(K.ops)) KroneckerOperator(A,domainspace(K),rs) end