inverse transforms (#59)

* inverse transforms

* 2D inverse plans for PiecewisePolynomial

* ContinuousPolynomial{1}

* 2D ContinuousPolynomial{1} transform

* C^(1) inverse transform

* differentiation works

* Dirichlet plan_transform

* Update Project.toml

* tests pass

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "PiecewiseOrthogonalPolynomials"
 uuid = "4461d12d-4663-4550-8580-cb764c85e20f"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.2.5"
+version = "0.3"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -21,11 +21,11 @@ QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
 
 [compat]
 ArrayLayouts = "1.0.12"
-BandedMatrices = "1"
+BandedMatrices = "1.3"
 BlockArrays = "0.16.25"
-BlockBandedMatrices = "0.12.2"
-ClassicalOrthogonalPolynomials = "0.12"
-ContinuumArrays = "0.17"
+BlockBandedMatrices = "0.12.9"
+ClassicalOrthogonalPolynomials = "0.12.1"
+ContinuumArrays = "0.17.1"
 FillArrays = "1.0"
 InfiniteArrays = "0.13"
 InfiniteLinearAlgebra = "0.7"

examples/2d.jl

@@ -0,0 +1,75 @@
+###
+# This example shows how we can expand functions on a 2d rectangle using plan_transform. 
+# This is a "manual mode" interface that might in the future be cleaned up to a user-facing
+# interface.
+###
+
+
+using PiecewiseOrthogonalPolynomials, Plots, Test
+using PiecewiseOrthogonalPolynomials: plan_grid_transform, plan_transform
+
+r = range(-1,1; length=3)
+P = ContinuousPolynomial{0}(r)
+
+# degrees in each dimension
+M,N = Block(80),Block(81)
+((x,y),pl) = plan_grid_transform(P, (M,N))
+
+
+
+# x and y are matrices of points in each dimension, with the columns corresponding
+# to different elements . This is needed to leverage
+# fast 1D transforms acting on arrays.  Thus to transform we make a tensor grid as
+# a 4-tensor. The way to do this is as follows:
+
+f = (x,y) -> exp(x*cos(10y))
+F = f.(x, reshape(y, 1, 1, size(y)...))
+
+# F contains our function samples. to plot we need to reduce back to vector/matrices.
+# But an issue is that the sample points are in the reverse order. Thus we need to do
+# some reordering:
+surface(vec(x[end:-1:1,:]), vec(y[end:-1:1,:]), reshape(F[end:-1:1,:,end:-1:1,:], length(x), length(y))')
+
+
+# We can now transform the samples to coefficient space:
+X = pl * F
+
+# We approximate the original function to high-order
+@test P[0.1,Block(1):M]' * X * P[0.2,Block(1):N] ≈ f(0.1,0.2)
+
+# Further we can transform back:
+@test pl \ X ≈ F
+
+
+# If we use ContinuousPolynomial{1} we can differentiate. This can be done as follows:
+
+C = ContinuousPolynomial{1}(r)
+
+# Compute coefficients of f in tensor product of C^(1)
+# This has one less block so we decrease M and N by one:
+
+pl_C = plan_transform(C, (M-1,N-1))
+F_C = f.(x, reshape(y, 1, 1, size(y,1),:))
+X_C = pl_C * F_C
+
+@test C[0.1,Block(1):M-1]' * X_C * C[0.2,Block(1):N-1] ≈ f(0.1,0.2)
+
+# Make the 1D differentiation matrices from C to P:
+
+D_x = (P \ diff(C))[Block(1):M, Block(1):M-1]
+D_y = (P \ diff(C))[Block(1):N, Block(1):N-1]
+
+
+# We now compare this to an analytical derivative at a poiunt
+
+f_xy = (x,y) -> -10sin(10y)*exp(x*cos(10y)) - 10cos(10y)*x*sin(10y)*exp(x*cos(10y))
+@test P[0.1,Block(1):M]'*(D_x*X_C*D_y')*P[0.2,Block(1):N] ≈ f_xy(0.1,0.2)
+
+# We can transform back to get the values on a large grid:
+
+F_xy = pl \ (D_x*X_C*D_y')
+
+@test F_xy ≈ f_xy.(x, reshape(y,1,1,size(y)...))
+
+surface(vec(x[end:-1:1,:]), vec(y[end:-1:1,:]), reshape(F_xy[end:-1:1,:,end:-1:1,:], length(x), length(y))')
+

src/PiecewiseOrthogonalPolynomials.jl

@@ -6,10 +6,10 @@ import BandedMatrices: _BandedMatrix
 import BlockArrays: BlockSlice, block, blockindex, blockvec
 import BlockBandedMatrices: _BandedBlockBandedMatrix, AbstractBandedBlockBandedMatrix, subblockbandwidths, blockbandwidths, AbstractBandedBlockBandedLayout, layout_replace_in_print_matrix
 import ClassicalOrthogonalPolynomials: grid, plotgrid, ldiv, pad, adaptivetransform_ldiv, Plan, singularities, basis_singularities, singularitiesbroadcast
-import ContinuumArrays: @simplify, factorize, TransformFactorization, AbstractBasisLayout, MemoryLayout, layout_broadcasted, ExpansionLayout, basis, plan_transform, plan_grid_transform, grammatrix, weaklaplacian
+import ContinuumArrays: @simplify, factorize, TransformFactorization, AbstractBasisLayout, MemoryLayout, layout_broadcasted, ExpansionLayout, basis, plan_transform, plan_grid_transform, grammatrix, weaklaplacian, MulPlan, InvPlan
 import LazyArrays: paddeddata
 import LazyBandedMatrices: BlockBroadcastMatrix, BlockVec, BandedLazyLayouts, AbstractLazyBandedBlockBandedLayout, UpperOrLowerTriangular
-import Base: axes, getindex, +, -, *, /, ==, \, OneTo, oneto, replace_in_print_matrix, copy, diff, getproperty, adjoint, transpose, tail, _sum
+import Base: size, axes, getindex, +, -, *, /, ==, \, OneTo, oneto, replace_in_print_matrix, copy, diff, getproperty, adjoint, transpose, tail, _sum, inv
 import LinearAlgebra: BlasInt
 import MatrixFactorizations: reversecholcopy
 

src/continuouspolynomial.jl

@@ -51,10 +51,113 @@ plan_transform(P::ContinuousPolynomial{0}, szs::NTuple{N,Union{Int,Block{1}}}, d
 for grd in (:grid, :plotgrid)
     @eval begin
         $grd(C::ContinuousPolynomial, n::Block{1}) = $grd(PiecewisePolynomial(C), n) 
-        $grd(C::ContinuousPolynomial, n::Int) = $grd(PiecewisePolynomial(C), n) 
+        $grd(C::ContinuousPolynomial, n::Int) = $grd(C, findblock(axes(C,2), n))
     end
 end
 
+struct ContinuousPolynomialTransform{T, Pl<:Plan{T}, RRs, Dims} <: Plan{T}
+    legendretransform::Pl
+    R::RRs
+    dims::Dims
+end
+
+
+
+
+function plan_transform(C::ContinuousPolynomial{1,T}, Ns::NTuple{N,Block{1}}, dims=ntuple(identity,Val(N))) where {N,T}
+    P = Legendre{T}()
+    W = Weighted(Jacobi{T}(1,1))
+    # TODO: this is unnecessarily not triangular which makes it much slower than necessary and prevents in-place.
+    # However, the speed of the Legendre transform will far exceed this so its not high priority.
+    # Probably using Ultraspherical(-1/2) would be better.
+    R̃ = [[T[1 1; -1 1]/2; Zeros{T}(∞,2)] (P \ W)]
+    ContinuousPolynomialTransform(plan_transform(ContinuousPolynomial{0}(C), Ns .+ 1, dims).F, 
+                                  InvPlan(map(N -> R̃[1:Int(N),1:Int(N)], Ns .+ 1), _doubledims(dims...)),
+                                  dims)
+end
+
+_tensorsize2contblocks() = ()
+_tensorsize2contblocks(M,N, Ns...) = (Vcat(N+1, Fill(N, M-2)), _tensorsize2contblocks(Ns...)...)
+
+
+function *(Pl::ContinuousPolynomialTransform{T,<:Any,<:Any,Int}, X::AbstractMatrix{T}) where T
+    dat = Pl.R * (Pl.legendretransform*X)
+    cfs = PseudoBlockArray{T}(undef,  _tensorsize2contblocks(size(X)...)...)
+    dims = Pl.dims
+    @assert dims == 1
+
+    M,N = size(X,1), size(X,2)
+    if size(dat,1) ≥ 1
+        cfs[Block(1)[1]] = dat[1,1]
+        for j = 1:N-1
+            isapprox(dat[2,j], dat[1,j+1]; atol=1000*M*eps()) || throw(ArgumentError("Discontinuity in data on order of $(abs(dat[2,j]- dat[1,j+1]))."))
+        end
+        for j = 1:N
+            cfs[Block(1)[j+1]] = dat[2,j]
+        end
+    end
+    cfs[Block.(2:M-1)] .= vec(dat[3:end,:]')
+    cfs
+end
+
+function \(Pl::ContinuousPolynomialTransform{T,<:Any,<:Any,Int}, cfs::AbstractVector{T}) where T
+    dims = Pl.dims
+    @assert dims == 1
+    
+    M,N = blocksize(cfs,1)+1, size(axes(cfs,1)[Block(1)],1)-1
+    dat = Matrix{T}(undef, M, N)
+    
+    if M ≥ 1
+        dat[1,1] = cfs[Block(1)[1]]
+        for j = 1:N-1
+            dat[2,j] = dat[1,j+1] = cfs[Block(1)[j+1]]
+        end
+        dat[2,end] = cfs[Block(1)[N+1]]
+    end
+
+    for j = 1:N, k = 3:M
+        dat[k,j] = cfs[Block(k-1)[j]]
+    end
+
+    Pl.legendretransform \ (Pl.R \ dat)
+end
+
+
+
+function _contpolyinds2blocks(k, j)
+    k == 1 && return Block(1)[j]
+    k == 2 && return Block(1)[j+1]
+    Block(k-1)[j]
+end
+function *(Pl::ContinuousPolynomialTransform{T}, X::AbstractArray{T,4}) where T
+    dat = Pl.R * (Pl.legendretransform*X)
+    cfs = PseudoBlockArray{T}(undef,  _tensorsize2contblocks(size(X)...)...)
+    dims = Pl.dims
+    @assert dims == (1,2)
+
+    M,N,O,P = size(X)
+    for k = 1:M, j = 1:N, l = 1:O, m = 1:P
+        cfs[_contpolyinds2blocks(k,j), _contpolyinds2blocks(l,m)] = dat[k,j,l,m]
+    end
+    cfs
+end
+
+function \(Pl::ContinuousPolynomialTransform{T}, cfs::AbstractMatrix{T}) where T
+    M,N = blocksize(cfs,1)+1, size(axes(cfs,1)[Block(1)],1)-1
+    O,P = blocksize(cfs,2)+1, size(axes(cfs,2)[Block(1)],1)-1
+
+    dat = Array{T}(undef,  M, N, O, P)
+    dims = Pl.dims
+    @assert dims == (1,2)
+
+    for k = 1:M, j = 1:N, l = 1:O, m = 1:P
+        dat[k,j,l,m] = cfs[_contpolyinds2blocks(k,j), _contpolyinds2blocks(l,m)]
+    end
+
+    Pl.legendretransform \ (Pl.R \ dat)
+end
+
+
 function adaptivetransform_ldiv(Q::ContinuousPolynomial{1,V}, f::AbstractQuasiVector) where V
     T = promote_type(V, eltype(f))
     C₀ = ContinuousPolynomial{0,V}(Q)
@@ -88,12 +191,7 @@ end
 adaptivetransform_ldiv(Q::ContinuousPolynomial{1,V}, f::AbstractQuasiMatrix) where V =
     BlockBroadcastArray(hcat, (Q \ f[:,j] for j = axes(f,2))...)
 
-function grid(V::SubQuasiArray{T,2,<:ContinuousPolynomial{1},<:Tuple{Inclusion,BlockSlice}}) where {T}
-    P = parent(V)
-    _, JR = parentindices(V)
-    pts = P.points
-    grid(view(PiecewisePolynomial(Weighted(Jacobi{T}(1, 1)), pts), :, JR))
-end
+grid(P::ContinuousPolynomial{1}, M::Block{1}) = grid(ContinuousPolynomial{0}(P), M + 1)
 
 #######
 # Conversion

src/dirichlet.jl

@@ -48,12 +48,120 @@ function \(P::ContinuousPolynomial, Q::DirichletPolynomial)
     (P \ C) * (C \ Q)
 end
 
+###
+# transforms
+###
+
 
 function adaptivetransform_ldiv(Q::DirichletPolynomial{V}, f::AbstractQuasiVector) where V
     C = ContinuousPolynomial{1}(Q)
     (Q \ C) * (C \ f)
 end
 
+struct DirichletPolynomialTransform{T, Pl<:Plan{T}, RRs, Dims} <: Plan{T}
+    legendretransform::Pl
+    R::RRs
+    dims::Dims
+end
+
+DirichletPolynomialTransform(pl::ContinuousPolynomialTransform) =
+    DirichletPolynomialTransform(pl.legendretransform, pl.R, pl.dims)
+
+plan_transform(C::DirichletPolynomial{T}, Ns::NTuple{N,Block{1}}, dims=ntuple(identity,Val(N))) where {N,T} =
+    DirichletPolynomialTransform(plan_transform(ContinuousPolynomial{1}(C), Ns, dims))
+
+_tensorsize2dirichletblocks() = ()
+_tensorsize2dirichletblocks(M,N, Ns...) = (Vcat(N-1, Fill(N, M-2)), _tensorsize2dirichletblocks(Ns...)...)
+    
+    
+function *(Pl::DirichletPolynomialTransform{T,<:Any,<:Any,Int}, X::AbstractMatrix{T}) where T
+    dat = Pl.R * (Pl.legendretransform*X)
+    cfs = PseudoBlockArray{T}(undef,  _tensorsize2dirichletblocks(size(X)...)...)
+    dims = Pl.dims
+    @assert dims == 1
+
+    M,N = size(X,1), size(X,2)
+    if size(dat,1) ≥ 1
+        for j = 1:N-1
+            cfs[Block(1)[j]] = dat[2,j]
+        end
+    end
+    cfs[Block.(2:M-1)] .= vec(dat[3:end,:]')
+    cfs
+end
+
+function \(Pl::DirichletPolynomialTransform{T,<:Any,<:Any,Int}, cfs::AbstractVector{T}) where T
+    dims = Pl.dims
+    @assert dims == 1
+    
+    M,N = blocksize(cfs,1)+1, size(axes(cfs,1)[Block(1)],1)+1
+    dat = Matrix{T}(undef, M, N)
+    
+    if M ≥ 1
+        dat[1,1] = zero(T)
+        for j = 1:N-1
+            dat[2,j] = dat[1,j+1] = cfs[Block(1)[j]]
+        end
+        dat[2,end] = zero(T)
+    end
+
+    for j = 1:N, k = 3:M
+        dat[k,j] = cfs[Block(k-1)[j]]
+    end
+
+    Pl.legendretransform \ (Pl.R \ dat)
+end
+    
+    
+    
+function _dirichletpolyinds2blocks(k, j)
+    k == 1 && return Block(1)[j-1]
+    k == 2 && return Block(1)[j]
+    Block(k-1)[j]
+end
+function *(Pl::DirichletPolynomialTransform{T}, X::AbstractArray{T,4}) where T
+    dat = Pl.R * (Pl.legendretransform*X)
+    cfs = PseudoBlockArray{T}(undef,  _tensorsize2dirichletblocks(size(X)...)...)
+    dims = Pl.dims
+    @assert dims == (1,2)
+
+    M,N,O,P = size(X)
+    for k = 1:M, j = 1:N, l = 1:O, m = 1:P
+        k == j == 1 && continue
+        k == 2 && j == N && continue
+        l == m == 1 && continue
+        l == 2 && m == P && continue
+        cfs[_dirichletpolyinds2blocks(k,j), _dirichletpolyinds2blocks(l,m)] = dat[k,j,l,m]
+    end
+    cfs
+end
+    
+function \(Pl::DirichletPolynomialTransform{T}, cfs::AbstractMatrix{T}) where T
+    M,N = blocksize(cfs,1)+1, size(axes(cfs,1)[Block(1)],1)+1
+    O,P = blocksize(cfs,2)+1, size(axes(cfs,2)[Block(1)],1)+1
+
+    dat = Array{T}(undef,  M, N, O, P)
+    dims = Pl.dims
+    @assert dims == (1,2)
+
+    for k = 1:M, j = 1:N, l = 1:O, m = 1:P
+        if k == j == 1 ||
+           (k == 2 && j == N) ||
+            l == m == 1 ||
+           (l == 2 && m == P)
+           dat[k,j,l,m] = zero(T)
+        else
+            dat[k,j,l,m] = cfs[_dirichletpolyinds2blocks(k,j), _dirichletpolyinds2blocks(l,m)]
+        end
+    end
+
+    Pl.legendretransform \ (Pl.R \ dat)
+end
+    
+###
+# weaklaplacian
+###    
+
 function weaklaplacian(C::DirichletPolynomial{T,<:AbstractRange}) where T
     r = C.points
     N = length(r)
@@ -102,8 +210,8 @@ end
 
 for grd in (:grid, :plotgrid)
     @eval begin
-        $grd(C::DirichletPolynomial, n::Integer) = $grd(PiecewisePolynomial(C), n)
-        $grd(C::DirichletPolynomial, n::Block{1}) = $grd(PiecewisePolynomial(C), n)
+        $grd(C::DirichletPolynomial, n::Integer) = $grd(ContinuousPolynomial{1}(C.points), n)
+        $grd(C::DirichletPolynomial, n::Block{1}) = $grd(ContinuousPolynomial{1}(C.points), n)
     end
 end