Kronecker operators on BivariateFuns (#229)

* kronecker operators on bivariate fun

* remove unused methods

* bugfix in getindex

* version bump to v0.7.19

Author: jishnub

Project.toml

@@ -1,6 +1,6 @@
 name = "ApproxFunBase"
 uuid = "fbd15aa5-315a-5a7d-a8a4-24992e37be05"
-version = "0.7.18"
+version = "0.7.19"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"

src/Multivariate/ProductFun.jl

@@ -5,6 +5,29 @@
 
 export ProductFun
 
+"""
+    ProductFun(f, space::TensorSpace; [tol=eps()])
+
+Represent a bivariate function `f(x,y)` as a univariate expansion over the second space,
+with the coefficients being functions in the first space.
+```math
+f\\left(x,y\\right)=\\sum_{i}f_{i}\\left(x\\right)b_{i}\\left(y\\right),
+```
+where ``b_{i}\\left(y\\right)`` represents the ``i``-th basis function in the space over ``y``.
+
+# Examples
+```jldoctest
+julia> P = ProductFun((x,y)->x*y, Chebyshev() ⊗ Chebyshev());
+
+julia> P(0.1, 0.2) ≈ 0.1 * 0.2
+true
+
+julia> coefficients(P) # power only at the (1,1) Chebyshev mode
+2×2 Matrix{Float64}:
+ 0.0  0.0
+ 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
@@ -22,6 +45,22 @@ Base.size(f::ProductFun) = (size(f,1),size(f,2))
 
 ## Construction in an AbstractProductSpace via a Matrix of coefficients
 
+"""
+    ProductFun(coeffs::AbstractMatrix{T}, sp::AbstractProductSpace; [tol=100eps(T)], [chopping=false]) where {T<:Number}
+
+Represent a bivariate function `f` in terms of the coefficient matrix `coeffs`,
+where the coefficients are obtained using a bivariate
+transform of the function `f` in the basis `sp`.
+
+# Examples
+```jldoctest
+julia> P = ProductFun([0 0; 0 1], Chebyshev() ⊗ Chebyshev()) # corresponds to (x,y) -> x*y
+ProductFun on Chebyshev() ⊗ Chebyshev()
+
+julia> P(0.1, 0.2) ≈ 0.1 * 0.2
+true
+```
+"""
 function ProductFun(cfs::AbstractMatrix{T},sp::AbstractProductSpace{Tuple{S,V},DD};
   tol::Real=100eps(T),chopping::Bool=false) where {S<:UnivariateSpace,V<:UnivariateSpace,T<:Number,DD}
     if chopping
@@ -36,8 +75,25 @@ function ProductFun(cfs::AbstractMatrix{T},sp::AbstractProductSpace{Tuple{S,V},D
 end
 
 ## Construction in a ProductSpace via a Vector of Funs
-
-function ProductFun(M::AbstractVector{VFun{S,T}},dy::V) where {S<:UnivariateSpace,V<:UnivariateSpace,T<:Number}
+"""
+    ProductFun(M::AbstractVector{<:Fun{<:UnivariateSpace}}, sp::UnivariateSpace)
+
+Represent a bivariate function `f(x,y)` in terms of the univariate coefficient functions from `M`.
+The function `f` may be reconstructed as
+```math
+f\\left(x,y\\right)=\\sum_{i}M_{i}\\left(x\\right)b_{i}\\left(y\\right),
+```
+where ``b_{i}\\left(y\\right)`` represents the ``i``-th basis function for the space `sp`.
+
+# Examples
+```jldoctest
+julia> P = ProductFun([zeros(Chebyshev()), Fun(Chebyshev())], Chebyshev()); # corresponds to (x,y)->x*y
+
+julia> P(0.1, 0.2) ≈ 0.1 * 0.2
+true
+```
+"""
+function ProductFun(M::AbstractVector{VFun{S,T}}, dy::V) where {S<:UnivariateSpace,V<:UnivariateSpace,T<:Number}
     prodsp = ProductSpace(map(space, M), dy)
     ProductFun{S,V,typeof(prodsp),T}(copy(M), prodsp)
 end

src/Operators/Operator.jl

@@ -472,12 +472,13 @@ true
 ```
 """
 getindex(B::Operator,f::Fun) = B*Multiplication(domainspace(B),f)
-getindex(B::Operator,f::LowRankFun{S,M,SS,T}) where {S,M,SS,T} = mapreduce(i->f.A[i]*B[f.B[i]],+,1:rank(f))
-getindex(B::Operator{BT},f::ProductFun{S,V,SS,T}) where {BT,S,V,SS,T} =
-    mapreduce(i->f.coefficients[i]*B[Fun(f.space[2],[zeros(promote_type(BT,T),i-1);
-                                            one(promote_type(BT,T))])],
-                +,1:length(f.coefficients))
-
+getindex(B::Operator,f::LowRankFun) = mapreduce(((fAi,fBi),) -> fAi * B[fBi], +, zip(f.A, f.B))
+function getindex(B::Operator{BT}, f::ProductFun{S,V,SS,T}) where {BT,S,V,SS,T}
+    TBF = promote_type(BT,T)
+    sp2 = factors(f.space)[2]
+    mapreduce(((ind, fi),)-> fi * B[Fun(sp2, [zeros(TBF,i-1); one(TBF)])], +,
+                enumerate(f.coefficients))
+end
 
 
 # Convenience for wrapper ops

src/PDE/KroneckerOperator.jl

@@ -353,9 +353,9 @@ function BandedBlockBandedMatrix(S::SubOperator{T,KroneckerOperator{SS,V,DS,RS,
     A,B = KO.ops
 
 
-    AA = strictconvert(BandedMatrix, view(A, Block(1):last(KR),Block(1):last(JR)))::BandedMatrix{eltype(S)}
+    AA = strictconvert(BandedMatrix, view(A, Block(1):last(KR),Block(1):last(JR)))
     Al,Au = bandwidths(AA)
-    BB = strictconvert(BandedMatrix, view(B, Block(1):last(KR),Block(1):last(JR)))::BandedMatrix{eltype(S)}
+    BB = strictconvert(BandedMatrix, view(B, Block(1):last(KR),Block(1):last(JR)))
     Bl,Bu = bandwidths(BB)
     λ,μ = subblockbandwidths(ret)
 
@@ -364,7 +364,7 @@ function BandedBlockBandedMatrix(S::SubOperator{T,KroneckerOperator{SS,V,DS,RS,
         Bs = view(ret, K, J)
         l = min(Al,Bu+n-m,λ)
         u = min(Au,Bl+m-n,μ)
-        @inbounds for j=1:m, k=max(1,j-u):min(n,j+l)
+        for j=1:m, k=max(1,j-u):min(n,j+l)
             a = AA[k,j]
             b = BB[n-k+1,m-j+1]
             c = a*b
@@ -422,3 +422,25 @@ function mul_coefficients(A::SubOperator{T,KKO,Tuple{UnitRange{Int},UnitRange{In
     M = P[KR,JR]
     M*pad(b, size(M,2))
 end
+
+Base.getindex(A::KroneckerOperator, B::MultivariateFun) = A[Fun(B)]
+function Base.getindex(K::KroneckerOperator, f::LowRankFun)
+    op1, op2 = K.ops
+    mapreduce(((A,B),)-> op1[A] ⊗ op2[B], +, zip(f.A, f.B))
+end
+function Base.getindex(K::KroneckerOperator, B::ProductFun)
+    op1, op2 = K.ops
+    S2 = factors(B.space)[2]
+    T = cfstype(B)
+    mapreduce(((ind, fi),)-> op1[fi] ⊗ op2[Fun(S2, [zeros(T, ind-1); one(T)])], +,
+        enumerate(B.coefficients))
+end
+for F in [:LowRankFun, :ProductFun, :MultivariateFun]
+    for O in [:DerivativeWrapper, :DefiniteIntegralWrapper]
+        @eval Base.getindex(K::$O{<:KroneckerOperator}, f::$F) = K.op[f]
+        @eval (*)(A::KroneckerOperator, B::$F) = A * Fun(B)
+        @eval (*)(A::$O{<:KroneckerOperator}, B::$F) = A.op * B
+        @eval (*)(A::$F, B::KroneckerOperator) = Fun(A) * B
+        @eval (*)(A::$F, B::$O{<:KroneckerOperator}) = Fun(A) * B.op
+    end
+end

src/hacks.jl

@@ -18,15 +18,22 @@ Fun(f::Fun) = f # Fun of Fun should be like a conversion
 """
     Fun(f)
 
-Return `Fun(f, Chebyshev())`
+Return `Fun(f, space)` by choosing an appropriate `space` for the function.
+For univariate functions, `space` is chosen to be `Chebyshev()`, whereas for
+multivariate functions, it is a tensor product of `Chebyshev()` spaces.
 
 # Examples
 ```jldoctest
-julia> f = Fun(x->x^2)
+julia> f = Fun(x -> x^2)
 Fun(Chebyshev(), [0.5, 0.0, 0.5])
 
 julia> f(0.1) == (0.1)^2
 true
+
+julia> f = Fun((x,y) -> x + y);
+
+julia> f(0.1, 0.2) ≈ 0.3
+true
 ```
 """
 function Fun(fin::Function)