ExpansionLayout-Lazy ambiguity (#150)

* ExpansionLayout-Lazy ambiguity

* Use expansion for Derivative

* Update operators.jl

* Update bases.jl

* Support _sum

* __sum -> sum_layout

* Update Project.toml

* Drop old Julia versions

* test sum and dot and basis

* increase coverage

Author: dlfivefifty

.github/workflows/ci.yml

@@ -10,9 +10,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.6'
-          - '1'
-          - '^1.9.0-0'
+          - '1.9'
         os:
           - ubuntu-latest
           - macOS-latest

Project.toml

@@ -1,6 +1,6 @@
 name = "ContinuumArrays"
 uuid = "7ae1f121-cc2c-504b-ac30-9b923412ae5c"
-version = "0.12.6"
+version = "0.13"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
@@ -19,21 +19,21 @@ RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
-AbstractFFTs = "1"
-ArrayLayouts = "0.7.7, 0.8, 1"
-BandedMatrices = "0.16, 0.17"
+AbstractFFTs = "1.0"
+ArrayLayouts = "1.0"
+BandedMatrices = "0.17"
 BlockArrays = "0.16"
-DomainSets = "0.5, 0.6"
-FastTransforms = "0.13, 0.14"
-FillArrays = "0.12, 0.13, 1"
+DomainSets = "0.6"
+FastTransforms = "0.15"
+FillArrays = "1.0"
 InfiniteArrays = "0.12"
 Infinities = "0.1"
-IntervalSets = "0.5, 0.6, 0.7"
-LazyArrays = "0.22, 1"
-QuasiArrays = "0.9.6"
+IntervalSets = "0.7"
+LazyArrays = "1.0"
+QuasiArrays = "0.10"
 RecipesBase = "1.0"
 StaticArrays = "1.0"
-julia = "1.6"
+julia = "1.9"
 
 [extras]
 Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

src/ContinuumArrays.jl

@@ -8,7 +8,7 @@ import Base.Broadcast: materialize, BroadcastStyle, broadcasted, Broadcasted
 import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, colsupport, combine_mul_styles, AbstractArrayApplyStyle,
                         adjointlayout, arguments, _mul_arguments, call, broadcastlayout, layout_getindex, UnknownLayout,
                         sublayout, sub_materialize, ApplyLayout, BroadcastLayout, combine_mul_styles, applylayout,
-                        simplifiable, _simplify, AbstractLazyLayout, PaddedLayout
+                        simplifiable, _simplify, AbstractLazyLayout, PaddedLayout, simplify, Dot
 import LinearAlgebra: pinv, inv, dot, norm2, ldiv!, mul!
 import BandedMatrices: AbstractBandedLayout, _BandedMatrix
 import BlockArrays: block, blockindex, unblock, blockedrange, _BlockedUnitRange, _BlockArray
@@ -18,7 +18,7 @@ import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclu
                     QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat, QuasiArrayLayout,
                     ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle, AbstractQuasiLazyLayout,
                     LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle, _factorize, _cutdim,
-                    AbstractQuasiFill, UnionDomain, __sum, _cumsum, __cumsum, applylayout, _equals, layout_broadcasted, PolynomialLayout
+                    AbstractQuasiFill, UnionDomain, sum_size, sum_layout, _cumsum, cumsum_layout, applylayout, _equals, layout_broadcasted, PolynomialLayout, _dot
 import InfiniteArrays: Infinity, InfAxes
 import AbstractFFTs: Plan
 
@@ -47,6 +47,8 @@ function dot(x::Inclusion{T,<:AbstractInterval}, y::Inclusion{V,<:AbstractInterv
     convert(TV, b^3 - a^3)/3
 end
 
+sum(x::Inclusion{T,<:AbstractInterval}) where T = convert(T, width(x.domain))
+
 
 include("maps.jl")
 
@@ -97,4 +99,17 @@ include("bases/bases.jl")
 
 include("plotting.jl")
 
+###
+# sum/dot
+###
+
+sum_size(::Tuple{InfiniteCardinal{1}}, a, dims) = _sum(expand(a), dims)
+_dot(::InfiniteCardinal{1}, a, b) = dot(expand(a), expand(b))
+function copy(d::Dot{<:ExpansionLayout,<:ExpansionLayout,<:AbstractQuasiArray,<:AbstractQuasiArray})
+    a,b = d.A,d.B
+    P,c = basis(a),coefficients(a)
+    Q,d = basis(b),coefficients(b)
+    c' * (P'Q) * d
+end
+
 end

src/bases/bases.jl

@@ -280,6 +280,19 @@ A / A \\ f.(axes(A,1))
 """
 expand(A, f) = A * transform(A, f)
 
+"""
+    expand(v)
+
+finds a natural basis for a quasi-vector and expands
+in that basis.
+"""
+function expand(v)
+    P = basis(v)
+    ApplyQuasiArray(*, P, P \ v)
+end
+
+
+
 copy(L::Ldiv{<:AbstractBasisLayout}) = transform_ldiv(L.A, L.B)
 # TODO: redesign to use simplifiable(\, A, B)
 copy(L::Ldiv{<:AbstractBasisLayout,ApplyLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) = transform_ldiv(L.A, L.B)
@@ -318,9 +331,19 @@ _factorize(::WeightedBasisLayouts, wS, dims...; kws...) = WeightedFactorization(
 ##
 
 struct ExpansionLayout{Lay} <: AbstractLazyLayout end
-applylayout(::Type{typeof(*)}, ::Lay, ::Union{PaddedLayout,AbstractStridedLayout}) where Lay <: AbstractBasisLayout = ExpansionLayout{Lay}()
+applylayout(::Type{typeof(*)}, ::Lay, ::Union{PaddedLayout,AbstractStridedLayout,ZerosLayout}) where Lay <: AbstractBasisLayout = ExpansionLayout{Lay}()
+
+"""
+    basis(v)
+
+gives a basis for expanding given quasi-vector.
+"""
+basis(v) = basis_layout(MemoryLayout(v), v)
+
+basis_layout(::ExpansionLayout, v::ApplyQuasiArray{<:Any,N,typeof(*)}) where N = v.args[1]
+basis_layout(lay, v) = basis_axes(axes(v,1), v) # allow choosing a basis based on axes
+basis_axes(ax, v) = error("Overload for $ax")
 
-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)...)
 
@@ -333,6 +356,7 @@ end
 LazyArrays._mul_arguments(::ExpansionLayout, A) = LazyArrays._mul_arguments(ApplyLayout{typeof(*)}(), A)
 copy(L::Ldiv{Bas,<:ExpansionLayout}) where Bas<:AbstractBasisLayout = copy(Ldiv{Bas,ApplyLayout{typeof(*)}}(L.A, L.B))
 copy(L::Mul{<:ExpansionLayout,Lay}) where Lay = copy(Mul{ApplyLayout{typeof(*)},Lay}(L.A, L.B))
+copy(L::Mul{<:ExpansionLayout,Lay}) where Lay<:AbstractLazyLayout = copy(Mul{ApplyLayout{typeof(*)},Lay}(L.A, L.B))
 
 function _broadcastbasis(::typeof(+), _, _, a, b)
     try
@@ -374,8 +398,8 @@ function layout_broadcasted(::Tuple{Vararg{ExpansionLayout}}, ::typeof(+), fs...
     P * +(_plus_P_ldiv_Ps_cs(P, Ps, cs)...)  # +((Ref(P) .\ Ps .* cs)...)
 end
 
-function layout_broadcasted(::NTuple{2,ExpansionLayout}, ::typeof(*), a, f)
-    axes(a,1) == axes(f,1) || throw(DimensionMismatch())
+function layout_broadcasted(::Tuple{Any,ExpansionLayout}, ::typeof(*), a, f)
+    axes(a)[1] == axes(f)[1] || throw(DimensionMismatch())
     P,c = arguments(f)
     (a .* P) * c
 end
@@ -575,24 +599,22 @@ end
 ####
 # sum
 ####
-
-
-function __sum(::SubBasisLayout, Vm, dims)
+function sum_layout(::SubBasisLayout, Vm, dims)
     @assert dims == 1
     sum(parent(Vm); dims=dims)[:,parentindices(Vm)[2]]
 end
 
-__sum(::AdjointBasisLayout, Vm::AbstractQuasiMatrix, dims) = permutedims(sum(Vm'; dims=(isone(dims) ? 2 : 1)))
+sum_layout(::AdjointBasisLayout, Vm::AbstractQuasiMatrix, dims) = permutedims(sum(Vm'; dims=(isone(dims) ? 2 : 1)))
 
 
-function __sum(::MappedBasisLayouts, V, dims)
+function sum_layout(::MappedBasisLayouts, V, dims)
     kr = basismap(V)
     @assert kr isa AbstractAffineQuasiVector
     sum(demap(V); dims=dims)/kr.A
 end
 
-__sum(::ExpansionLayout, A, dims) = __sum(ApplyLayout{typeof(*)}(), A, dims)
-__cumsum(::ExpansionLayout, A, dims) = __cumsum(ApplyLayout{typeof(*)}(), A, dims)
+sum_layout(::ExpansionLayout, A, dims) = sum_layout(ApplyLayout{typeof(*)}(), A, dims)
+cumsum_layout(::ExpansionLayout, A, dims) = cumsum_layout(ApplyLayout{typeof(*)}(), A, dims)
 
 include("basisconcat.jl")
 include("basiskron.jl")

src/bases/splines.jl

@@ -11,7 +11,7 @@ Spline{o}(S::Spline) where {o} = Spline{o}(S.points)
 
 for Typ in (:LinearSpline, :HeavisideSpline)
     STyp = string(Typ)
-    @eval function summary(io::IO, L::$Typ)
+    @eval function show(io::IO, L::$Typ)
         print(io, "$($STyp)(")
         print(IOContext(io, :limit=>true), L.points)
         print(io,")")

src/maps.jl

@@ -7,7 +7,7 @@ A subtype of `Map` is used as a one-to-one map between two domains
 via `view`. The domain of the map `m` is `axes(m,1)` and the range
 is `union(m)`.
 
-Maps must also overload `invmap` to give the inverse of the map, which 
+Maps must also overload `invmap` to give the inverse of the map, which
 is equivalent to `invmap(m)[x] == findfirst(isequal(x), m)`.
 """
 
@@ -40,6 +40,7 @@ end
 # Affine map represents A*x .+ b
 abstract type AbstractAffineQuasiVector{T,AA,X,B} <: Map{T} end
 
+show(io::IO, a::AbstractAffineQuasiVector) = summary(io, a)
 summary(io::IO, a::AbstractAffineQuasiVector) = print(io, "$(a.A) * $(a.x) .+ ($(a.b))")
 
 MemoryLayout(::Type{<:AbstractAffineQuasiVector}) = PolynomialLayout()
@@ -152,5 +153,11 @@ const AffineMappedQuasiMatrix = SubQuasiArray{<:Any, 2, <:Any, <:Tuple{AbstractA
 _sum(V::AffineMappedQuasiVector, ::Colon) = parentindices(V)[1].A \ sum(parent(V))
 
 # pretty print for bases
-summary(io::IO, P::AffineMappedQuasiMatrix) = print(io, "$(parent(P)) affine mapped to $(parentindices(P)[1].x.domain)")
-summary(io::IO, P::AffineMappedQuasiVector) = print(io, "$(parent(P)) affine mapped to $(parentindices(P)[1].x.domain)")
+show(io::IO, a::AffineMap) = summary(io, a)
+summary(io::IO, a::AffineMap) = print(io, "Affine map from $(a.domain) to $(a.range)")
+
+show(io::IO, P::Union{AffineMappedQuasiVector,AffineMappedQuasiMatrix}) = print(io, "$(parent(P)) affine mapped to $(parentindices(P)[1].x.domain)")
+function summary(io::IO, P::Union{AffineMappedQuasiVector,AffineMappedQuasiMatrix})
+    summary(io, parent(P))
+    print(io, " affine mapped to $(parentindices(P)[1].x.domain)")
+end

src/operators.jl

@@ -111,6 +111,7 @@ Derivative{T}(domain) where T = Derivative{T}(Inclusion(domain))
 
 Derivative(L::AbstractQuasiMatrix) = Derivative(axes(L,1))
 
+show(io::IO, a::Derivative) = summary(io, a)
 function summary(io::IO, D::Derivative)
     print(io, "Derivative(")
     summary(io,D.axis)
@@ -154,3 +155,9 @@ Identity(d::Inclusion) = QuasiDiagonal(d)
 
 @simplify *(D::Derivative, x::Inclusion) = ones(promote_type(eltype(D),eltype(x)), x)
 @simplify *(D::Derivative, c::AbstractQuasiFill) = zeros(promote_type(eltype(D),eltype(c)), axes(c,1))
+# @simplify *(D::Derivative, x::AbstractQuasiVector) = D * expand(x)
+
+struct OperatorLayout <: AbstractLazyLayout end
+MemoryLayout(::Type{<:Derivative}) = OperatorLayout()
+# copy(M::Mul{OperatorLayout, <:ExpansionLayout}) = simplify(M)
+# copy(M::Mul{OperatorLayout, <:AbstractLazyLayout}) = M.A * expand(M.B)

test/test_chebyshev.jl

@@ -1,8 +1,9 @@
-using ContinuumArrays, LinearAlgebra, FastTransforms, QuasiArrays, ArrayLayouts, Test
+using ContinuumArrays, LinearAlgebra, FastTransforms, QuasiArrays, ArrayLayouts, Base64, LazyArrays, Test
 import ContinuumArrays: Basis, Weight, Map, LazyQuasiArrayStyle, TransformFactorization,
                         ExpansionLayout, checkpoints, MappedBasisLayout, MappedWeightedBasisLayout,
-                        SubWeightedBasisLayout, WeightedBasisLayout, WeightLayout
+                        SubWeightedBasisLayout, WeightedBasisLayout, WeightLayout, basis
 
+using IntervalSets: AbstractInterval
 """
 This is a simple implementation of Chebyshev for testing. Use ClassicalOrthogonalPolynomials
     for the real implementation.
@@ -24,10 +25,23 @@ Base.getindex(::ChebyshevWeight, x::Float64) = 1/sqrt(1-x^2)
 Base.getindex(w::ChebyshevWeight, ::Inclusion) = w # TODO: make automatic
 
 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)
+ContinuumArrays.basis_axes(::Inclusion{<:Any,<:AbstractInterval}, v) = Chebyshev(100)
+function ContinuumArrays._sum(T::Chebyshev, dims)
+    n = 2:size(T,2)-1
+    [2; 0; @. ((1/(n+1) - 1/(n-1)) - ((-1)^(n+1)/(n+1) - (-1)^(n-1)/(n-1)))/2]'
+end
 
 # 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)
 
+ContinuumArrays.@simplify function *(A::QuasiAdjoint{<:Any,<:Chebyshev}, B::Chebyshev)
+    m,n = size(A,1),size(B,2)
+    T = promote_type(eltype(A), eltype(B))
+    f = (k,j) -> isodd(j-k) ? zero(T) : -(((1 + (-1)^(j + k))*(-1 + j^2 + k^2))/(j^4 + (-1 + k^2)^2 - 2j^2*(1 + k^2)))
+    f.(0:m-1, (0:n-1)')
+end
+
+
 struct QuadraticMap{T} <: Map{T} end
 struct InvQuadraticMap{T} <: Map{T} end
 
@@ -89,6 +103,7 @@ ContinuumArrays.invmap(::InvQuadraticMap{T}) where T = QuadraticMap{T}()
         y = affine(0..1, x)
 
         @test summary(T[y,:]) == "Chebyshev affine mapped to 0..1"
+        @test stringmime("text/plain", T[y,:]) == "Chebyshev(5) affine mapped to 0..1"
         @test MemoryLayout(wT[y,:]) isa MappedWeightedBasisLayout
         @test MemoryLayout(w[y] .* T[y,:]) isa MappedWeightedBasisLayout
         @test wT[y,:][[0.1,0.2],1:5] == (w[y] .* T[y,:])[[0.1,0.2],1:5] == (w .* T[:,1:5])[y,:][[0.1,0.2],:]
@@ -135,4 +150,17 @@ ContinuumArrays.invmap(::InvQuadraticMap{T}) where T = QuadraticMap{T}()
         ã = T * (T \ a)
         @test T \ (ã .* ã) ≈ [1.5,1,0.5,0,0]
     end
+
+    @testset "sum/dot" begin
+        @test sum(x) ≡ 2.0
+        @test sum(exp.(x)) ≈ ℯ - 1/ℯ
+        @test dot(x, x) ≈ 2/3
+        @test dot(exp.(x), x) ≈ 2/ℯ
+    end
+
+    @testset "Expansion * Lazy" begin
+        f = T * collect(1.0:5)
+        @test (f * ones(1,4))[0.1,:] == fill(f[0.1],4)
+        @test (f * BroadcastArray(exp, (1:4)'))[0.1,:] ≈ f[0.1] * exp.(1:4)
+    end
 end

test/test_maps.jl

@@ -74,6 +74,6 @@ import ContinuumArrays: AffineQuasiVector
 
     @testset "show" begin
         @test stringmime("text/plain", y) == "2.0 * Inclusion(0..1) .+ (-1.0)"
-        @test stringmime("text/plain", a) == "0.4 * Inclusion(-2..3) .+ (-0.2)"
+        @test stringmime("text/plain", a) == "Affine map from Inclusion(-2..3) to Inclusion(-1..1)"
     end
 end