Redesign grid to take in an n corresponding to the number of coefficients needed (#138)

* have grid work without subarray

* Update test_chebyshev.jl

* simplify grid overloading

* v0.12

* arr -> szs

* tests pass

* minor changes

* increase cov

* mul vectors

Author: dlfivefifty

Project.toml

@@ -1,6 +1,7 @@
 name = "ContinuumArrays"
 uuid = "7ae1f121-cc2c-504b-ac30-9b923412ae5c"
-version = "0.11.2"
+version = "0.12"
+
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

src/ContinuumArrays.jl

@@ -21,7 +21,7 @@ import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclu
                     AbstractQuasiFill, UnionDomain, __sum, _cumsum, __cumsum, applylayout, _equals, layout_broadcasted, PolynomialLayout
 import InfiniteArrays: Infinity, InfAxes
 
-export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, ℵ₁, Inclusion, Basis, grid, plotgrid, affine, .., transform, expand, plan_transform
+export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, ℵ₁, Inclusion, Basis, grid, plotgrid, affine, .., transform, expand, plan_transform, basis, coefficients 
 
 
 

src/bases/bases.jl

@@ -30,6 +30,9 @@ broadcastlayout(::Type{typeof(*)}, ::WeightLayout, ::Basis) where Basis<:Abstrac
 sublayout(::WeightLayout, _) = WeightLayout()
 sublayout(::AbstractBasisLayout, ::Type{<:Tuple{Map,AbstractVector}}) = MappedBasisLayout()
 
+# copy with an Inclusion can not be materialized
+copy(V::SubQuasiArray{<:Any,N,<:Basis,<:Tuple{Inclusion,Vararg{Any}}, trfl}) where {N,trfl}  = V
+
 
 ## Weighted basis interface
 unweighted(P::BroadcastQuasiMatrix{<:Any,typeof(*),<:Tuple{AbstractQuasiVector,AbstractQuasiMatrix}}) = last(P.args)
@@ -148,12 +151,26 @@ copy(L::Ldiv{<:MappedBasisLayouts,BroadcastLayout{typeof(*)},<:Any,<:AbstractQua
 
 
 # expansion
-_grid(_, P) = error("Overload Grid")
+_grid(_, P, n...) = error("Overload Grid")
+
+_grid(::MappedBasisLayout, P, n...) = invmap(parentindices(P)[1])[grid(demap(P), n...)]
+_grid(::SubBasisLayout, P::AbstractQuasiMatrix, n) = grid(parent(P), maximum(parentindices(P)[2][n]))
+_grid(::SubBasisLayout, P::AbstractQuasiMatrix) = grid(parent(P), maximum(parentindices(P)[2]))
+_grid(::WeightedBasisLayouts, P, n...) = grid(unweighted(P), n...)
 
-_grid(::MappedBasisLayout, P) = invmap(parentindices(P)[1])[grid(demap(P))]
-_grid(::SubBasisLayout, P) = grid(parent(P))
-_grid(::WeightedBasisLayouts, P) = grid(unweighted(P))
-grid(P) = _grid(MemoryLayout(P), P)
+
+"""
+    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.
+"""
+grid(P, n...) = _grid(MemoryLayout(P), P, n...)
+
+
+# values(f) = 
 
 
 struct TransformFactorization{T,Grid,Plan} <: Factorization{T}
@@ -239,15 +256,18 @@ function *(P::InvPlan, X::AbstractArray)
 end
 
 
-function plan_grid_transform(L, arr, dims=1:ndims(arr))
+function plan_grid_transform(L, szs::NTuple{N,Int}, dims=1:N) where N
     p = grid(L)
     p, InvPlan(factorize(L[p,:]), dims)
 end
 
-plan_transform(P, arr, dims...) = plan_grid_transform(P, arr, dims...)[2]
+plan_grid_transform(L, arr::AbstractArray{<:Any,N}, dims=1:N) where N = 
+    plan_grid_transform(L, size(arr), dims)
+
+plan_transform(P, szs, dims...) = plan_grid_transform(P, szs, dims...)[2]
 
 _factorize(::AbstractBasisLayout, L, dims...; kws...) =
-    TransformFactorization(plan_grid_transform(L, Array{eltype(L)}(undef, size(L,2), dims...), 1)...)
+    TransformFactorization(plan_grid_transform(L, (size(L,2), dims...), 1)...)
 
 
 
@@ -273,7 +293,7 @@ _sub_factorize(::Tuple{Any,Int}, (kr,jr)::Tuple{Any,OneTo}, L, dims...; kws...)
 
 # ∞-dimensional parents need to use transforms. For now we assume the size of the transform is equal to the size of the truncation
 _sub_factorize(::Tuple{Any,Any}, (kr,jr)::Tuple{Any,OneTo}, L, dims...; kws...) =
-    TransformFactorization(plan_grid_transform(parent(L), Array{eltype(L)}(undef, last(jr), dims...), 1)...)
+    TransformFactorization(plan_grid_transform(parent(L), (last(jr), dims...), 1)...)
 
 # If jr is not OneTo we project
 _sub_factorize(::Tuple{Any,Any}, (kr,jr), L, dims...; kws...) =
@@ -377,6 +397,7 @@ applylayout(::Type{typeof(*)}, ::Lay, ::Union{PaddedLayout,AbstractStridedLayout
 
 basis(v::ApplyQuasiArray{<:Any,N,typeof(*)}) where N = v.args[1]
 coefficients(v::ApplyQuasiArray{<:Any,N,typeof(*),<:Tuple{Any,Any}}) where N = v.args[2]
+coefficients(v::ApplyQuasiArray{<:Any,N,typeof(*),<:Tuple{Any,Any,Vararg{Any}}}) where N = ApplyArray(*, tail(v.args)...)
 
 
 function unweighted(lay::ExpansionLayout, a)

src/bases/splines.jl

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

src/plotting.jl

@@ -1,27 +1,45 @@
 
+const MAX_PLOT_POINTS = 10_000 # above this rendering is too slow
 
-_mul_plotgrid(_, args) = grid(first(args))
-_mul_plotgrid(::Tuple{Any,PaddedLayout}, (P,c)) = plotgrid(P[:,colsupport(c)])
+
+"""
+    plotgrid(P, n...)
+
+returns a grid of points suitable for plotting. This may include
+endpoints or singular points not included in `grid`. `n` specifies
+the number of coefficients.
+"""
+
+plotgrid(P, n...) = _plotgrid(MemoryLayout(P), P, n...)
+_plotgrid(lay, P, n=size(P,2)) = grid(P, min(n,MAX_PLOT_POINTS))
+
+_plotgrid(::WeightedBasisLayouts, wP, n...) = plotgrid(unweighted(wP), n...)
+_plotgrid(::MappedBasisLayout, P, n...) = invmap(parentindices(P)[1])[plotgrid(demap(P), n...)]
+_plotgrid(::SubBasisLayout, P::AbstractQuasiMatrix, n) = plotgrid(parent(P), maximum(parentindices(P)[2][n]))
+_plotgrid(::SubBasisLayout, P::AbstractQuasiMatrix) = plotgrid(parent(P), maximum(parentindices(P)[2]))
+
+
+_mul_plotgrid(_, args) = _plotgrid(UnknownLayout(), first(args))
+_mul_plotgrid(::Tuple{Any,PaddedLayout}, (P,c)) = plotgrid(P, maximum(colsupport(c)))
 
 function _plotgrid(lay::ExpansionLayout, P)
     args = arguments(lay,P)
     _mul_plotgrid(map(MemoryLayout,args), args)
 end
 
-_plotgrid(_, P) = grid(P)
-
-_plotgrid(::WeightedBasisLayouts, wP) = plotgrid(unweighted(wP))
-_plotgrid(::MappedBasisLayout, P) = invmap(parentindices(P)[1])[plotgrid(demap(P))]
-
-plotgrid(g) = _plotgrid(MemoryLayout(g), g)
-
 _split_svec(x) = (x,)
 _split_svec(x::AbstractArray{<:StaticVector{2}}) = (map(first,x), map(last,x))
 
 plotvalues(g::AbstractQuasiVector, x) = g[x]
 plotvalues(g::AbstractQuasiMatrix, x) = g[x,:]
+plotvalues(g::AbstractQuasiArray) = plotvalues(g, plotgrid(g))
 
-@recipe function f(g::AbstractQuasiArray)
+function plotgridvalues(g)
     x = plotgrid(g)
-    tuple(_split_svec(x)..., plotvalues(g,x))
+    x, plotvalues(g,x)
+end
+
+@recipe function f(g::AbstractQuasiArray)
+    x,v = plotgridvalues(g)
+    tuple(_split_svec(x)..., v)
 end

test/runtests.jl

@@ -1,7 +1,7 @@
 using ContinuumArrays, QuasiArrays, IntervalSets, DomainSets, FillArrays, LinearAlgebra, BandedMatrices, InfiniteArrays, Test, Base64, RecipesBase
 import ContinuumArrays: ℵ₁, materialize, AffineQuasiVector, BasisLayout, AdjointBasisLayout, SubBasisLayout, ℵ₁,
                         MappedBasisLayout, AdjointMappedBasisLayout, MappedWeightedBasisLayout, TransformFactorization, Weight, WeightedBasisLayout, SubWeightedBasisLayout, WeightLayout,
-                        basis, invmap, Map, checkpoints, _plotgrid, mul
+                        basis, invmap, Map, checkpoints, _plotgrid, mul, plotvalues
 import QuasiArrays: SubQuasiArray, MulQuasiMatrix, Vec, Inclusion, QuasiDiagonal, LazyQuasiArrayApplyStyle, LazyQuasiArrayStyle
 import LazyArrays: MemoryLayout, ApplyStyle, Applied, colsupport, arguments, ApplyLayout, LdivStyle, MulStyle
 
@@ -89,9 +89,16 @@ include("test_basisconcat.jl")
     rep = RecipesBase.apply_recipe(Dict{Symbol, Any}(), L)
     @test rep[1].args == (L.points,L[L.points,:])
 
+    rep = RecipesBase.apply_recipe(Dict{Symbol, Any}(), L[:,1:3])
+    @test rep[1].args == (L.points,L[L.points,1:3])
+
+    @test plotgrid(L[:,1:3],3) == grid(L[:,1:3]) == grid(L[:,1:3],3) == L.points
+    
+
     u = L*randn(6)
     rep = RecipesBase.apply_recipe(Dict{Symbol, Any}(), u)
     @test rep[1].args == (L.points,u[L.points])
+    @test plotvalues(u) == u[plotgrid(u)]
 
     @testset "padded" begin
         u = L * Vcat(rand(3), Zeros(3))

test/test_chebyshev.jl

@@ -1,9 +1,11 @@
-using ContinuumArrays, LinearAlgebra, FastTransforms, QuasiArrays, Test
-import ContinuumArrays: Basis, Weight, Map, LazyQuasiArrayStyle, TransformFactorization, ExpansionLayout
+using ContinuumArrays, LinearAlgebra, FastTransforms, QuasiArrays, ArrayLayouts, Test
+import ContinuumArrays: Basis, Weight, Map, LazyQuasiArrayStyle, TransformFactorization,
+                        ExpansionLayout, checkpoints, MappedBasisLayout, MappedWeightedBasisLayout,
+                        SubWeightedBasisLayout, WeightedBasisLayout, WeightLayout
 
 """
 This is a simple implementation of Chebyshev for testing. Use ClassicalOrthogonalPolynomials
-for the real implementation.
+    for the real implementation.
 """
 struct Chebyshev <: Basis{Float64}
     n::Int
@@ -14,17 +16,14 @@ struct ChebyshevWeight <: Weight{Float64} end
 Base.:(==)(::Chebyshev, ::Chebyshev) = true
 Base.:(==)(::ChebyshevWeight, ::ChebyshevWeight) = true
 Base.axes(T::Chebyshev) = (Inclusion(-1..1), Base.OneTo(T.n))
-ContinuumArrays.grid(T::Chebyshev) = chebyshevpoints(Float64, T.n, Val(1))
+ContinuumArrays.grid(T::Chebyshev, n...) = chebyshevpoints(Float64, T.n, Val(1))
 Base.axes(T::ChebyshevWeight) = (Inclusion(-1..1),)
 
 Base.getindex(::Chebyshev, x::Float64, n::Int) = cos((n-1)*acos(x))
 Base.getindex(::ChebyshevWeight, x::Float64) = 1/sqrt(1-x^2)
 Base.getindex(w::ChebyshevWeight, ::Inclusion) = w # TODO: make automatic
 
-LinearAlgebra.factorize(L::Chebyshev) =
-    TransformFactorization(grid(L), plan_chebyshevtransform(Array{Float64}(undef, size(L,2))))
-LinearAlgebra.factorize(L::Chebyshev, n) =
-    TransformFactorization(grid(L), plan_chebyshevtransform(Array{Float64}(undef, size(L,2),n),1))
+ContinuumArrays.plan_grid_transform(L::Chebyshev, szs::NTuple{N,Int}, dims=1:N) where N = grid(L), plan_chebyshevtransform(Array{eltype(L)}(undef, szs...), dims)
 
 # This is wrong but just for tests
 QuasiArrays.layout_broadcasted(::Tuple{ExpansionLayout,Any}, ::typeof(*), a::ApplyQuasiVector{<:Any,typeof(*),<:Tuple{Chebyshev,Any}}, b::Chebyshev) = b * Matrix(I, 5, 5)

test/test_splines.jl

@@ -1,7 +1,7 @@
 using ContinuumArrays, LinearAlgebra, Base64, FillArrays, QuasiArrays, BandedMatrices, Test
 using QuasiArrays: ApplyQuasiArray, ApplyStyle, MemoryLayout, mul, MulQuasiMatrix, Vec
 import LazyArrays: MulStyle, LdivStyle, arguments, applied, apply
-import ContinuumArrays: basis, AdjointBasisLayout, ExpansionLayout, BasisLayout, SubBasisLayout, AdjointMappedBasisLayout, MappedBasisLayout, coefficients
+import ContinuumArrays: basis, AdjointBasisLayout, ExpansionLayout, BasisLayout, SubBasisLayout, AdjointMappedBasisLayout, MappedBasisLayout
 
 @testset "Splines" begin
     @testset "HeavisideSpline" begin
@@ -28,6 +28,8 @@ import ContinuumArrays: basis, AdjointBasisLayout, ExpansionLayout, BasisLayout,
         @test MemoryLayout(typeof(H)) == BasisLayout()
         @test ApplyStyle(*, typeof(H), typeof([1,2])) isa MulStyle
 
+        @test copy(H[:,1:2]) == H[:,1:2]
+
         f = H*[1,2]
         @test f isa ApplyQuasiArray
         @test axes(f) == (Inclusion(1.0..3.0),)
@@ -495,4 +497,10 @@ import ContinuumArrays: basis, AdjointBasisLayout, ExpansionLayout, BasisLayout,
         end
         @test PX ≈ X
     end
+
+    @testset "Mul coefficients" begin
+        L = LinearSpline(0:5)
+        u = ApplyQuasiArray(*, L, randn(6,5), randn(5))
+        @test coefficients(u) ≈ L \ u
+    end
 end