Refine tensor transform interface (#170)

* Refine tensor transform interface

* work on tests

* Make plan_transform and grid only overloads

* tests pass

* coverage

* add tests

* add docs

* Update index.md

* Update bases.jl

* Update bases.jl

* Update test_splines.jl

* fix tests

* increase cov

* plan_transform(P, Block(K,J))

* plan_grid_transform with Block(K,J)

* Update bases.jl

Author: dlfivefifty

.github/workflows/documentation.yml

@@ -0,0 +1,26 @@
+name: Documentation
+
+on:
+  push:
+    branches:
+      - master # update to match your development branch (master, main, dev, trunk, ...)
+    tags: '*'
+  pull_request:
+
+jobs:
+  build:
+    permissions:
+      contents: write
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/setup-julia@v1
+        with:
+          version: '1.9'
+      - name: Install dependencies
+        run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
+      - name: Build and deploy
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # If authenticating with GitHub Actions token
+          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # If authenticating with SSH deploy key
+        run: julia --project=docs/ docs/make.jl

Project.toml

@@ -1,6 +1,6 @@
 name = "ContinuumArrays"
 uuid = "7ae1f121-cc2c-504b-ac30-9b923412ae5c"
-version = "0.16.4"
+version = "0.17"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"

README.md

@@ -3,6 +3,7 @@ A package for representing quasi arrays with continuous dimensions
 
 [![Build Status](https://github.com/JuliaApproximation/ContinuumArrays.jl/workflows/CI/badge.svg)](https://github.com/JuliaApproximation/ContinuumArrays.jl/actions)
 [![codecov](https://codecov.io/gh/JuliaApproximation/ContinuumArrays.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/JuliaApproximation/ContinuumArrays.jl)
+[![docs](https://img.shields.io/badge/docs-dev-blue.svg)](https://JuliaApproximation.github.io/ContinuumArrays.jl/dev)
 
 
 A quasi array as implemented in [QuasiArrays.jl](https://github.com/JuliaApproximation/QuasiArrays.jl) is a 

docs/Project.toml

@@ -0,0 +1,5 @@
+[deps]
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+
+[compat]
+Documenter = "1"

docs/make.jl

@@ -0,0 +1,12 @@
+using Documenter, ContinuumArrays
+
+makedocs(
+    modules = [ContinuumArrays],
+    sitename="ContinuumArrays.jl",
+    pages = Any[
+        "Home" => "index.md"])
+
+deploydocs(
+    repo = "github.com/JuliaApproximation/ContinuumArrays.jl.git",
+    push_preview = true
+)

docs/src/index.md

@@ -0,0 +1,75 @@
+# ContinuumArrays.jl
+A Julia package for working with bases as quasi-arrays
+
+## The Basis interface
+
+To add your own bases, subtype `Basis` and overload the following routines:
+
+1. `axes(::MyBasis) = (Inclusion(mydomain), OneTo(mylength))`
+2. `grid(::MyBasis, ::Integer)`
+3. `getindex(::MyBasis, x, ::Integer)`
+4. `==(::MyBasis, ::MyBasis)`
+
+Optional overloads are:
+
+1. `plan_transform(::MyBasis, sizes::NTuple{N,Int}, region)` to plan a transform, for a tensor
+product of the specified sizes, similar to `plan_fft`.
+2. `diff(::MyBasis, dims=1)` to support differentiation and `Derivative`. 
+3. `grammatrix(::MyBasis)` to support `Q'Q`. 
+4. `ContinuumArrays._sum(::MyBasis, dims)` and `ContinuumArrays._cumsum(::MyBasis, dims)` to support definite and indefinite integeration.
+
+
+## Routines
+
+
+```@docs
+Derivative
+```
+```@docs
+transform
+```
+```@docs
+expand
+```
+```@docs
+grid
+```
+
+
+## Interal Routines
+
+```@docs
+ContinuumArrays.TransformFactorization
+```
+
+```@docs
+ContinuumArrays.AbstractConcatBasis
+```
+```@docs
+ContinuumArrays.MulPlan
+```
+```@docs
+ContinuumArrays.PiecewiseBasis
+```
+```@docs
+ContinuumArrays.MappedFactorization
+```
+```@docs
+ContinuumArrays.basis
+```
+```@docs
+ContinuumArrays.InvPlan
+```
+```@docs
+ContinuumArrays.VcatBasis
+```
+```@docs
+ContinuumArrays.WeightedFactorization
+```
+```@docs
+ContinuumArrays.HvcatBasis
+```
+```@docs
+ContinuumArrays.ProjectionFactorization
+```
+```

src/bases/bases.jl

@@ -153,50 +153,68 @@ copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(*)}}) = _broadcast_mul
 copy(L::Ldiv{<:MappedBasisLayouts,BroadcastLayout{typeof(*)}}) = _broadcast_mul_ldiv(map(MemoryLayout,arguments(L.B)), L.A, L.B)
 
 
-# expansion
-grid_layout(_, P, n...) = error("Overload Grid")
 
-grid_layout(::MappedBasisLayout, P, n...) = invmap(parentindices(P)[1])[grid(demap(P), n...)]
-grid_layout(::SubBasisLayout, P::AbstractQuasiMatrix, n) = grid(parent(P), maximum(parentindices(P)[2][n]))
-grid_layout(::SubBasisLayout, P::AbstractQuasiMatrix) = grid(parent(P), maximum(parentindices(P)[2]))
-grid_layout(::WeightedBasisLayouts, P, n...) = grid(unweighted(P), n...)
 
 
 """
     grid(P, n...)
 
 Creates a grid of points. if `n` is unspecified it will
 be sufficient number of points to determine `size(P,2)`
-coefficients. Otherwise its enough points to determine `n`
-coefficients.
+coefficients. If `n` is an integer or `Block` its enough points to determine `n`
+coefficients. If `n` is a tuple then it returns a tuple of grids corresponding to a
+tensor-product. That is, a 5⨱6 2D transform would be
+```julia
+(x,y) = grid(P, (5,6))
+plan_transform(P, (5,6)) * f.(x, y')
+```
+and a 5×6×7 3D transform would be
+```julia
+(x,y,z) = grid(P, (5,6,7))
+plan_transform(P, (5,6,7)) * f.(x, y', reshape(z,1,1,:))
+```
 """
-grid(P, n...) = grid_layout(MemoryLayout(P), P, n...)
+grid(P, n::Block{1}) = grid_layout(MemoryLayout(P), P, n...)
+grid(P, n::Integer) = grid_layout(MemoryLayout(P), P, n...)
+grid(L, B::Block) = grid(L, Block.(B.n)) # grid(L, Block(2,3)) == grid(L, (Block(2), Block(3))
+grid(L, ns::Tuple) = grid.(Ref(L), ns)
+grid(L) = grid(L, size(L,2))
 
+grid_layout(_, P, n) = grid_axis(axes(P,2), P, n)
 
-# values(f) = 
+grid_axis(::OneTo, P, n::Block) = grid(P, size(P,2))
 
+grid_layout(::MappedBasisLayout, P, n) = invmap(parentindices(P)[1])[grid(demap(P), n)]
+grid_layout(::SubBasisLayout, P::AbstractQuasiMatrix, n) = grid(parent(P), parentindices(P)[2][n])
+grid_layout(::WeightedBasisLayouts, P, n) = grid(unweighted(P), n)
 
 
-function plan_grid_transform_layout(lay, L, szs::NTuple{N,Int}, dims=1:N) where N
-    p = grid(L)
-    p, InvPlan(factorize(L[p,:]), dims)
-end
+# Default transform is just solve least squares on a grid
+# note this computes the grid an extra time.
+mapfactorize(L, n::Integer) = factorize(L[grid(L,n),OneTo(n)])
+mapfactorize(L, n::Block{1}) = factorize(L[grid(L,n),Block.(OneTo(Int(n)))])
 
-function plan_grid_transform_layout(::MappedBasisLayout, L, szs::NTuple{N,Int}, dims=1:N) where N
-    x,F = plan_grid_transform(demap(L), szs, dims)
-    invmap(parentindices(L)[1])[x], F
+mapfactorize(L, ns) = map(n -> mapfactorize(L,n), ns)
+function plan_transform_layout(lay, L, szs::NTuple{N,Union{Int,Block{1}}}, dims=ntuple(identity,Val(N))) where N
+    dimsz = getindex.(Ref(szs), dims) # get the sizes of transformed dimensions
+    InvPlan(mapfactorize(L, dimsz), dims)
 end
+plan_transform_layout(::MappedBasisLayout, L, szs::NTuple{N,Union{Int,Block{1}}}, dims=ntuple(identity,Val(N))) where N = plan_transform(demap(L), szs, dims)
+plan_transform(L, szs::NTuple{N,Union{Int,Block{1}}}, dims=ntuple(identity,Val(N))) where N = plan_transform_layout(MemoryLayout(L), L, szs, dims)
 
-plan_grid_transform(L, szs::NTuple{N,Int}, dims=1:N) where N = plan_grid_transform_layout(MemoryLayout(L), L, szs, dims)
+plan_transform(L, arr::AbstractArray, dims...) = plan_transform(L, size(arr), dims...)
+plan_transform(L, lng::Union{Integer,Block{1}}, dims...) = plan_transform(L, (lng,), dims...)
+plan_transform(L) = plan_transform(L, size(L,2))
 
-plan_grid_transform(L, arr::AbstractArray, dims...) = plan_grid_transform(L, size(arr), dims...)
-plan_grid_transform(L, lng::Union{Integer,Block{1}}, dims...) = plan_grid_transform(L, (lng,), dims...)
+plan_transform(L, B::Block, dims...) = plan_transform(L, Block.(B.n), dims...) # grid(L, Block(2,3)) == grid(L, (Block(2), Block(3))
+    
 
-plan_transform(P, szs, dims...) = plan_grid_transform(P, szs, dims...)[2]
 
-_factorize(::AbstractBasisLayout, L, dims...; kws...) =
-    TransformFactorization(plan_grid_transform(L, (size(L,2), dims...), 1)...)
+plan_grid_transform(P, szs::NTuple{N,Union{Integer,Block{1}}}, dims=ntuple(identity,Val(N))) where N = grid(P, getindex.(Ref(szs), dims)), plan_transform(P, szs, dims)
+plan_grid_transform(P, lng::Union{Integer,Block{1}}, dims=1) = plan_grid_transform(P, (lng,), dims)
+plan_grid_transform(P, B::Block{N}, dims=ntuple(identity,Val(N))) where N = plan_grid_transform(P, Block.(B.n), dims)
 
+_factorize(::AbstractBasisLayout, L, dims...; kws...) = TransformFactorization(plan_grid_transform(L, (size(L,2), dims...), 1)...)
 
 
 """

src/bases/splines.jl

@@ -47,8 +47,9 @@ function getindex(B::HeavisideSpline{T}, x::Number, k::Int) where T
     return zero(T)
 end
 
-grid(L::HeavisideSpline, n...) = L.points[1:end-1] .+ diff(L.points)/2
-grid(L::LinearSpline, n...) = L.points
+# Splines sample same number of points regardless of length.
+grid(L::HeavisideSpline, ::Integer) = L.points[1:end-1] .+ diff(L.points)/2
+grid(L::LinearSpline, ::Integer) = L.points
 
 ## Sub-bases
 

src/plans.jl

@@ -31,52 +31,53 @@ end
 
 Takes a factorization and supports it applied to different dimensions.
 """
-struct InvPlan{T, Fact, Dims} <: Plan{T}
-    factorization::Fact
+struct InvPlan{T, Facts<:Tuple, Dims} <: Plan{T}
+    factorizations::Facts
     dims::Dims
 end
 
-InvPlan(fact, dims) = InvPlan{eltype(fact), typeof(fact), typeof(dims)}(fact, dims)
+InvPlan(fact::Tuple, dims) = InvPlan{eltype(fact), typeof(fact), typeof(dims)}(fact, dims)
+InvPlan(fact, dims) = InvPlan((fact,), dims)
 
-size(F::InvPlan, k...) = size(F.factorization, k...)
+size(F::InvPlan) = size.(F.factorizations, 1)
 
 
-function *(P::InvPlan{<:Any,<:Any,Int}, x::AbstractVector)
+function *(P::InvPlan{<:Any,<:Tuple,Int}, x::AbstractVector)
     @assert P.dims == 1
-    P.factorization \ x
+    only(P.factorizations) \ x # Only a single factorization when dims isa Int
 end
 
-function *(P::InvPlan{<:Any,<:Any,Int}, X::AbstractMatrix)
+function *(P::InvPlan{<:Any,<:Tuple,Int}, X::AbstractMatrix)
     if P.dims == 1
-        P.factorization \ X
+        only(P.factorizations) \ X  # Only a single factorization when dims isa Int
     else
         @assert P.dims == 2
-        permutedims(P.factorization \ permutedims(X))
+        permutedims(only(P.factorizations) \ permutedims(X))
     end
 end
 
-function *(P::InvPlan{<:Any,<:Any,Int}, X::AbstractArray{<:Any,3})
+function *(P::InvPlan{<:Any,<:Tuple,Int}, X::AbstractArray{<:Any,3})
     Y = similar(X)
     if P.dims == 1
         for j in axes(X,3)
-            Y[:,:,j] = P.factorization \ X[:,:,j]
+            Y[:,:,j] = only(P.factorizations) \ X[:,:,j]
         end
     elseif P.dims == 2
         for k in axes(X,1)
-            Y[k,:,:] = P.factorization \ X[k,:,:]
+            Y[k,:,:] = only(P.factorizations) \ X[k,:,:]
         end
     else
         @assert P.dims == 3
         for k in axes(X,1), j in axes(X,2)
-            Y[k,j,:] = P.factorization \ X[k,j,:]
+            Y[k,j,:] = only(P.factorizations) \ X[k,j,:]
         end
     end
     Y
 end
 
 function *(P::InvPlan, X::AbstractArray)
     for d in P.dims
-        X = InvPlan(P.factorization, d) * X
+        X = InvPlan(P.factorizations[d], d) * X
     end
     X
 end
@@ -87,54 +88,55 @@ end
 
 Takes a matrix and supports it applied to different dimensions.
 """
-struct MulPlan{T, Fact, Dims} <: Plan{T}
-    matrix::Fact
+struct MulPlan{T, Fact<:Tuple, Dims} <: Plan{T}
+    matrices::Fact
     dims::Dims
 end
 
-MulPlan(fact, dims) = MulPlan{eltype(fact), typeof(fact), typeof(dims)}(fact, dims)
+MulPlan(mats::Tuple, dims) = MulPlan{eltype(mats), typeof(mats), typeof(dims)}(mats, dims)
+MulPlan(mats::AbstractMatrix, dims) = MulPlan((mats,), dims)
 
-function *(P::MulPlan{<:Any,<:Any,Int}, x::AbstractVector)
+function *(P::MulPlan{<:Any,<:Tuple,Int}, x::AbstractVector)
     @assert P.dims == 1
-    P.matrix * x
+    only(P.matrices) * x
 end
 
-function *(P::MulPlan{<:Any,<:Any,Int}, X::AbstractMatrix)
+function *(P::MulPlan{<:Any,<:Tuple,Int}, X::AbstractMatrix)
     if P.dims == 1
-        P.matrix * X
+        only(P.matrices) * X
     else
         @assert P.dims == 2
-        permutedims(P.matrix * permutedims(X))
+        permutedims(only(P.matrices) * permutedims(X))
     end
 end
 
-function *(P::MulPlan{<:Any,<:Any,Int}, X::AbstractArray{<:Any,3})
+function *(P::MulPlan{<:Any,<:Tuple,Int}, X::AbstractArray{<:Any,3})
     Y = similar(X)
     if P.dims == 1
         for j in axes(X,3)
-            Y[:,:,j] = P.matrix * X[:,:,j]
+            Y[:,:,j] = only(P.matrices) * X[:,:,j]
         end
     elseif P.dims == 2
         for k in axes(X,1)
-            Y[k,:,:] = P.matrix * X[k,:,:]
+            Y[k,:,:] = only(P.matrices) * X[k,:,:]
         end
     else
         @assert P.dims == 3
         for k in axes(X,1), j in axes(X,2)
-            Y[k,j,:] = P.matrix * X[k,j,:]
+            Y[k,j,:] = only(P.matrices) * X[k,j,:]
         end
     end
     Y
 end
 
 function *(P::MulPlan, X::AbstractArray)
     for d in P.dims
-        X = MulPlan(P.matrix, d) * X
+        X = MulPlan(P.matrices[d], d) * X
     end
     X
 end
 
-*(A::AbstractMatrix, P::MulPlan) = MulPlan(A*P.matrix, P.dims)
+*(A::AbstractMatrix, P::MulPlan) = MulPlan(Ref(A) .* P.matrices, P.dims)
 
-inv(P::MulPlan) = InvPlan(factorize(P.matrix), P.dims)
-inv(P::InvPlan) = MulPlan(P.factorization, P.dims)
+inv(P::MulPlan) = InvPlan(map(factorize,P.matrices), P.dims)
+inv(P::InvPlan) = MulPlan(convert.(Matrix,P.factorizations), P.dims)