Multivariate functions (#185)

* first version for multivariate functions of type trivialtensorizer

* performance improvements

* added simple test cases

* fix test and extension

* comment broken test out for old Julia versions

* fixed ordering of tensorization

* simd change and vector type

* small changes

* enabled broken test

* fixed unbound parameters

* type fix

* type fixes

* add Polynomials for testing in Project toml

* invoked 2D case

* move multivariate tests to orthogonal polynomial repository

* delete multivariate

* fixed block type and rewritten trivialtensor iterator

* remove Orthogonal polynomials from Project toml

* adding TensorIteratorFun

* changed project toml

* rename TrivialTensorFun and add check if all spaces in Tensorspace are trivial

* performance optimization

* added ProductFun which is not working as a comment

Co-authored-by: Benjamin Zanger <[email protected]>
Co-authored-by: Jishnu Bhattacharya <[email protected]>

Author: benjione

Project.toml

@@ -8,6 +8,7 @@ BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
 Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
+Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
 DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
@@ -32,6 +33,7 @@ BandedMatrices = "0.16, 0.17"
 BlockArrays = "0.14, 0.15, 0.16"
 BlockBandedMatrices = "0.10, 0.11"
 Calculus = "0.5"
+Combinatorics = "1.0.2"
 DSP = "0.6, 0.7"
 DomainSets = "0.5"
 DualNumbers = "0.6.2"

src/ApproxFunBase.jl

@@ -35,6 +35,8 @@ import Base.Broadcast: BroadcastStyle, Broadcasted, AbstractArrayStyle,
 
 import Statistics: mean
 
+import Combinatorics: multiexponents
+
 import LinearAlgebra: BlasInt, BlasFloat, norm, ldiv!, mul!, det, cross,
               qr, qr!, rank, isdiag, istril, istriu, issymmetric,
               Tridiagonal, diagm, diagm_container, factorize,

src/Multivariate/Multivariate.jl

@@ -30,6 +30,7 @@ include("VectorFun.jl")
 include("TensorSpace.jl")
 include("LowRankFun.jl")
 include("ProductFun.jl")
+include("TrivialTensorFun.jl")
 
 
 arglength(f)=length(Base.uncompressed_ast(f.code.def).args[1])

src/Multivariate/ProductFun.jl

@@ -5,6 +5,7 @@
 
 export ProductFun
 
+
 """
     ProductFun(f, space::TensorSpace; [tol=eps()])
 
@@ -28,6 +29,7 @@ julia> coefficients(P) # power only at the (1,1) Chebyshev mode
  0.0  1.0
 ```
 """
+
 struct ProductFun{S<:UnivariateSpace,V<:UnivariateSpace,SS<:AbstractProductSpace,T} <: BivariateFun{T}
     coefficients::Vector{VFun{S,T}}     # coefficients are in x
     space::SS
@@ -268,7 +270,8 @@ function coefficients(f::ProductFun,ox::Space,oy::Space)
 end
 
 (f::ProductFun)(x,y) = evaluate(f,x,y)
-(f::ProductFun)(x,y,z) = evaluate(f,x,y,z)
+# ProductFun does only support BivariateFunctions, this function below just does not work
+# (f::ProductFun)(x,y,z) = evaluate(f,x,y,z)
 
 coefficients(f::ProductFun,ox::TensorSpace) = coefficients(f,ox[1],ox[2])
 
@@ -303,7 +306,6 @@ canonicalevaluate(f::ProductFun,xx::AbstractVector,yy::AbstractVector) =
 
 
 evaluate(f::ProductFun,x,y) = canonicalevaluate(f,tocanonical(f,x,y)...)
-evaluate(f::ProductFun,x,y,z) = canonicalevaluate(f,tocanonical(f,x,y,z)...)
 
 # TensorSpace does not use map
 evaluate(f::ProductFun{S,V,SS,T},x::Number,::Colon) where {S<:UnivariateSpace,V<:UnivariateSpace,SS<:TensorSpace,T} =
@@ -313,6 +315,7 @@ evaluate(f::ProductFun{S,V,SS,T},x::Number,y::Number) where {S<:UnivariateSpace,
     evaluate(f,x,:)(y)
 
 
+
 evaluate(f::ProductFun,x) = evaluate(f,x...)
 
 *(c::Number,f::F) where {F<:ProductFun} = F(c*f.coefficients,f.space)

src/Multivariate/TensorSpace.jl

@@ -24,17 +24,76 @@ struct Tensorizer{DMS<:Tuple}
     blocks::DMS
 end
 
+const Tensorizer2D{AA, BB} = Tensorizer{Tuple{AA, BB}}
 const TrivialTensorizer{d} = Tensorizer{NTuple{d,Ones{Int,1,Tuple{OneToInf{Int}}}}}
 
 Base.eltype(a::Tensorizer) = NTuple{length(a.blocks),Int}
 Base.eltype(::Tensorizer{<:NTuple{d}}) where {d} = NTuple{d,Int}
 dimensions(a::Tensorizer) = map(sum,a.blocks)
 Base.length(a::Tensorizer) = mapreduce(sum,*,a.blocks)
 
+
+function start(a::TrivialTensorizer{d}) where {d}
+    if d==2
+        return invoke(start, Tuple{Tensorizer2D}, a)
+    else
+        # ((block_dim_1, block_dim_2,...), (itaration_number, iterator, iterator_state)), (itemssofar, length)
+        return (ones(Int, d),(0, nothing, nothing)), (0,length(a))
+    end
+end
+
+function next(a::TrivialTensorizer{d}, iterator_tuple) where {d}
+
+    if d==2
+        return invoke(next, Tuple{Tensorizer2D, Tuple}, a, iterator_tuple)
+    end
+
+    (block, (j, iterator, iter_state)), (i,tot) = iterator_tuple
+
+
+    @inline function check_block_finished()
+        if iterator === nothing
+            return true
+        end
+        # there are N-1 over d-1 combinations in a block
+        amount_combinations_block = binomial(sum(block)-1, d-1)
+        # check if all combinations have been iterated over
+        amount_combinations_block <= j
+    end
+
+    ret = reverse(block)
+
+    if check_block_finished()   # end of new block
+
+        # set up iterator for new block
+        current_sum = sum(block)
+        iterator = multiexponents(d, current_sum+1-d)
+        iter_state = nothing
+        j = 0
+    end
+
+    # increase block, or initialize new block
+    res, iter_state = iterate(iterator, iter_state)
+    block .= res.+1
+    j = j+1
+
+    ret, ((block, (j, iterator, iter_state)), (i,tot))
+end
+
+
+function done(a::TrivialTensorizer{d}, iterator_tuple) where {d}
+    if d==2
+        return invoke(done, Tuple{Tensorizer2D, Tuple}, a, iterator_tuple)
+    end
+    (_, (i,tot)) = iterator_tuple
+    return i ≥ tot
+end
+
+
 # (blockrow,blockcol), (subrow,subcol), (rowshift,colshift), (numblockrows,numblockcols), (itemssofar, length)
-start(a::Tensorizer{Tuple{AA,BB}}) where {AA,BB} = (1,1), (1,1), (0,0), (a.blocks[1][1],a.blocks[2][1]), (0,length(a))
+start(a::Tensorizer2D{AA, BB}) where {AA,BB} = (1,1), (1,1), (0,0), (a.blocks[1][1],a.blocks[2][1]), (0,length(a))
 
-function next(a::Tensorizer{Tuple{AA,BB}}, ((K,J), (k,j), (rsh,csh), (n,m), (i,tot))) where {AA,BB}
+function next(a::Tensorizer2D{AA, BB}, ((K,J), (k,j), (rsh,csh), (n,m), (i,tot))) where {AA,BB}
     ret = k+rsh,j+csh
     if k==n && j==m  # end of block
         if J == 1 || K == length(a.blocks[1])   # end of new block
@@ -59,7 +118,7 @@ function next(a::Tensorizer{Tuple{AA,BB}}, ((K,J), (k,j), (rsh,csh), (n,m), (i,t
 end
 
 
-done(a::Tensorizer, ((K,J), (k,j), (rsh,csh), (n,m), (i,tot))) = i ≥ tot
+done(a::Tensorizer2D, ((K,J), (k,j), (rsh,csh), (n,m), (i,tot))) = i ≥ tot
 
 iterate(a::Tensorizer) = next(a, start(a))
 function iterate(a::Tensorizer, st)
@@ -104,6 +163,14 @@ block(ci::CachedIterator{T,TrivialTensorizer{2}},k::Int) where {T} =
 block(::TrivialTensorizer{2},n::Int) =
     Block(floor(Integer,sqrt(2n) + 1/2))
 
+function block(::TrivialTensorizer{d},n::Int) where {d}
+    order::Int = 0
+    while binomial(order+d, d) < n
+        order = order + 1
+    end
+    return Block(order+1)
+end
+
 block(sp::Tensorizer{<:Tuple{<:AbstractFill{S},<:AbstractFill{T}}},n::Int) where {S,T} =
     Block(floor(Integer,sqrt(2floor(Integer,(n-1)/(getindex_value(sp.blocks[1])*getindex_value(sp.blocks[2])))+1) + 1/2))
 _cumsum(x) = cumsum(x)
@@ -211,6 +278,10 @@ struct TensorSpace{SV,D,R} <:AbstractProductSpace{SV,D,R}
     spaces::SV
 end
 
+# Tensorspace of 2 univariate spaces
+const TensorSpace2D{AA, BB, D,R} = TensorSpace{<:Tuple{AA, BB}, D, R} where {AA<:UnivariateSpace, BB<:UnivariateSpace}
+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)...)
 
@@ -473,6 +544,8 @@ end
 
 fromtensor(S::Space,M::AbstractMatrix) = fromtensor(tensorizer(S),M)
 totensor(S::Space,M::AbstractVector) = totensor(tensorizer(S),M)
+totensor(SS::TensorSpace{<:NTuple{d}},M::AbstractVector) where {d} = 
+        if d>2; totensoriterator(tensorizer(SS),M) else totensor(tensorizer(SS),M) end
 
 function fromtensor(it::Tensorizer,M::AbstractMatrix)
     n,m=size(M)
@@ -496,19 +569,27 @@ function totensor(it::Tensorizer,M::AbstractVector)
     B=block(it,n)
     ds = dimensions(it)
 
+    #ret=zeros(eltype(M),[sum(it.blocks[i][1:min(B.n[1],length(it.blocks[i]))]) for i=1:length(it.blocks)]...)
+    
     ret=zeros(eltype(M),sum(it.blocks[1][1:min(B.n[1],length(it.blocks[1]))]),
                         sum(it.blocks[2][1:min(B.n[1],length(it.blocks[2]))]))
+
     k=1
-    for (K,J) in it
+    for index in it
         if k > n
             break
         end
-        ret[K,J] = M[k]
+        ret[index...] = M[k]
         k += 1
     end
     ret
 end
 
+@inline function totensoriterator(it::TrivialTensorizer{d},M::AbstractVector) where {d}
+    B=block(it,length(M))
+    return it, M, B
+end
+
 for OP in (:block,:blockstart,:blockstop)
     @eval begin
         $OP(s::TensorSpace, ::PosInfinity) = ℵ₀
@@ -542,10 +623,12 @@ end
 
 itransform(sp::TensorSpace,cfs::AbstractVector) = vec(itransform!(sp,coefficientmatrix(Fun(sp,cfs))))
 
-evaluate(f::AbstractVector,S::AbstractProductSpace,x) = ProductFun(totensor(S,f),S)(x...)
-evaluate(f::AbstractVector,S::AbstractProductSpace,x,y) = ProductFun(totensor(S,f),S)(x,y)
-
+# 2D evaluation functions
+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...)
 
 coefficientmatrix(f::Fun{<:AbstractProductSpace}) = totensor(space(f),f.coefficients)
 

src/Multivariate/TrivialTensorFun.jl

@@ -0,0 +1,45 @@
+
+
+
+struct TrivialTensorFun{d, SS<:TensorSpaceND{d}, T<:Number} <: MultivariateFun{T, d}
+    space::SS
+    coefficients::Vector{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
+    TrivialTensorFun(sp, cfs, iter, blk)
+end
+
+(f::TrivialTensorFun)(x...) = evaluate(f, x...)
+
+# TensorSpace evaluation
+function evaluate(f::TrivialTensorFun{d, SS, T},x...) where {d, SS, T}
+    highest_order = f.orders.n[1]
+    n = length(f.coefficients)
+
+    # this could be lazy evaluated for the sparse case
+    A = T[Fun(f.space.spaces[i], [zeros(T, k);1])(x[i]) for k=0:highest_order, i=1:d]
+    result::T = 0
+    coef_counter::Int = 1
+    for i in f.iterator
+        tmp = f.coefficients[coef_counter]
+        if tmp != 0
+            tmp_res = 1
+            for k=1:d
+                tmp_res *= A[i[k], k]
+            end
+            result += tmp * tmp_res
+        end
+        coef_counter += 1
+        if coef_counter > n
+            break
+        end
+    end
+    return result
+end