transform by factorization (#35)

* transform by factorization

* update quasiarrays

* Update ContinuumArrays.jl

* Update Project.toml

* add WeightedBasis

* mass matrix C.O.V.

* Update bases.jl

* Update Project.toml

* Increase coverage

* v0.2

Author: dlfivefifty

.gitignore

@@ -3,3 +3,4 @@
 *.jl.mem
 deps/deps.jl
 .DS_Store
+Manifest.toml

.travis.yml

@@ -5,8 +5,6 @@ os:
   - osx
   - windows
 julia:
-  - 1.1
-  - 1.2
   - 1.3
   - nightly
 matrix:

Project.toml

@@ -1,24 +1,27 @@
 name = "ContinuumArrays"
 uuid = "7ae1f121-cc2c-504b-ac30-9b923412ae5c"
-version = "0.1"
+version = "0.2"
 
 [deps]
+ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
+FastTransforms = "057dd010-8810-581a-b7be-e3fc3b93f78c"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
 
 [compat]
-BandedMatrices = "0.14.1"
+BandedMatrices = "0.14.2"
 FillArrays = "0.8.2"
 IntervalSets = "0.3.2"
-LazyArrays = "0.14.7, 0.15"
-QuasiArrays = "0.0.6, 0.1"
-julia = "1.1"
+LazyArrays = "0.14.11"
+QuasiArrays = "0.1"
+julia = "1.3"
 
 [extras]
+FastTransforms = "057dd010-8810-581a-b7be-e3fc3b93f78c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]

src/ContinuumArrays.jl

@@ -15,7 +15,7 @@ import FillArrays: AbstractFill, getindex_value, SquareEye
 import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclusion, SubQuasiArray,
                     QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat, quasimulapplystyle,
                     ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle,
-                    LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle, quasildivapplystyle
+                    LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle, quasildivapplystyle, _factorize
 
 export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, fullmaterialize, ℵ₁, Inclusion, Basis, WeightedBasis, grid, transform, affine
 
@@ -181,9 +181,10 @@ affine(a, b) = affine(Inclusion(a), Inclusion(b))
 
 const QInfAxes = Union{Inclusion,AbstractAffineQuasiVector}
 
+
 sub_materialize(_, V::AbstractQuasiArray, ::Tuple{QInfAxes}) = V
 sub_materialize(_, V::AbstractQuasiArray, ::Tuple{QInfAxes,QInfAxes}) = V
-sub_materialize(_, V::AbstractQuasiArray, ::Tuple{<:Any,QInfAxes}) = V
+sub_materialize(_, V::AbstractQuasiArray, ::Tuple{Any,QInfAxes}) = V
 sub_materialize(_, V::AbstractQuasiArray, ::Tuple{QInfAxes,Any}) = V
 
 

src/bases/bases.jl

@@ -9,6 +9,7 @@ abstract type AbstractBasisLayout <: MemoryLayout end
 struct BasisLayout <: AbstractBasisLayout end
 struct SubBasisLayout <: AbstractBasisLayout end
 struct MappedBasisLayout <: AbstractBasisLayout end
+struct WeightedBasisLayout <: AbstractBasisLayout end
 
 abstract type AbstractAdjointBasisLayout <: MemoryLayout end
 struct AdjointBasisLayout <: AbstractAdjointBasisLayout end
@@ -21,8 +22,8 @@ MemoryLayout(::Type{<:Weight}) = WeightLayout()
 adjointlayout(::Type, ::BasisLayout) = AdjointBasisLayout()
 adjointlayout(::Type, ::SubBasisLayout) = AdjointSubBasisLayout()
 adjointlayout(::Type, ::MappedBasisLayout) = AdjointMappedBasisLayout()
-broadcastlayout(::Type{typeof(*)}, ::WeightLayout, ::BasisLayout) = BasisLayout()
-broadcastlayout(::Type{typeof(*)}, ::WeightLayout, ::SubBasisLayout) = SubBasisLayout()
+broadcastlayout(::Type{typeof(*)}, ::WeightLayout, ::BasisLayout) = WeightedBasisLayout()
+broadcastlayout(::Type{typeof(*)}, ::WeightLayout, ::SubBasisLayout) = WeightedBasisLayout()
 
 combine_mul_styles(::AbstractBasisLayout) = LazyQuasiArrayApplyStyle()
 combine_mul_styles(::AbstractAdjointBasisLayout) = LazyQuasiArrayApplyStyle()
@@ -88,18 +89,46 @@ end
 _grid(_, P) = error("Overload Grid")
 _grid(::MappedBasisLayout, P) = igetindex.(Ref(parentindices(P)[1]), grid(demap(P)))
 _grid(::SubBasisLayout, P) = grid(parent(P))
+_grid(::WeightedBasisLayout, P) = grid(last(P.args))
 grid(P) = _grid(MemoryLayout(typeof(P)), P)
 
-function transform(L)
+struct TransformFactorization{T,Grid,Plan,IPlan} <: Factorization{T}
+    grid::Grid
+    plan::Plan
+    iplan::IPlan
+end
+
+TransformFactorization(grid, plan) = 
+    TransformFactorization{promote_type(eltype(grid),eltype(plan)),typeof(grid),typeof(plan),Nothing}(grid, plan, nothing)
+
+
+TransformFactorization(grid, ::Nothing, iplan) = 
+    TransformFactorization{promote_type(eltype(grid),eltype(iplan)),typeof(grid),Nothing,typeof(iplan)}(grid, nothing, iplan)
+
+grid(T::TransformFactorization) = T.grid    
+
+\(a::TransformFactorization{<:Any,<:Any,Nothing}, b::AbstractQuasiVector) = a.iplan \  convert(Array, b[a.grid])
+\(a::TransformFactorization, b::AbstractQuasiVector) = a.plan * convert(Array, b[a.grid])
+
+\(a::TransformFactorization{<:Any,<:Any,Nothing}, b::AbstractVector) = a.iplan \  b
+\(a::TransformFactorization, b::AbstractVector) = a.plan * b
+
+function _factorize(::AbstractBasisLayout, L)
     p = grid(L)
-    p,L[p,:]
+    TransformFactorization(p, nothing, factorize(L[p,:]))
 end
 
-function transform_ldiv(A, B, _)
-    p,T = transform(A)
-    T \ convert(Array, B[p])
+struct ProjectionFactorization{T, FAC<:Factorization{T}, INDS} <: Factorization{T}
+    F::FAC
+    inds::INDS
 end
 
+\(a::ProjectionFactorization, b::AbstractQuasiVector) = (a.F \ b)[a.inds]
+\(a::ProjectionFactorization, b::AbstractVector) = (a.F \ b)[a.inds]
+
+_factorize(::SubBasisLayout, L) = ProjectionFactorization(factorize(parent(L)), parentindices(L)[2])
+
+transform_ldiv(A, B, _) = factorize(A) \ B
 transform_ldiv(A, B) = transform_ldiv(A, B, axes(A))
 
 copy(L::Ldiv{<:AbstractBasisLayout,<:Any,<:Any,<:AbstractQuasiVector}) =
@@ -116,7 +145,7 @@ end
 
 
 function copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(*)},<:AbstractQuasiMatrix,<:AbstractQuasiVector})
-    p,T = transform(L.A)
+    p,T = factorize(L.A)
     T \ L.B[p]
 end
 
@@ -126,6 +155,25 @@ end
 #     *(arguments(S)...)
 
 
+
+# mass matrix
+# y = p(x), dy = p'(x) * dx
+# \int_a^b f(y) g(y) dy = \int_{-1}^1 f(p(x))*g(p(x)) * p'(x) dx
+
+
+_sub_getindex(A, kr, jr) = A[kr, jr]
+_sub_getindex(A, ::Slice, ::Slice) = A
+
+function copy(M::QMul2{<:QuasiAdjoint{<:Any,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}},
+                        <:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}})
+    Ac, B = M.args
+    A = Ac'
+    PA,PB = parent(A),parent(B)
+    kr,jr = parentindices(B)
+    _sub_getindex((PA'PB)/kr.A,parentindices(A)[2],jr)
+end
+
+
 # Differentiation of sub-arrays
 function copy(M::QMul2{<:Derivative,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:Inclusion,<:Any}}})
     A, B = M.args
@@ -155,7 +203,7 @@ end
 
 # we represent as a Mul with a banded matrix
 sublayout(::AbstractBasisLayout, ::Type{<:Tuple{<:Inclusion,<:AbstractUnitRange}}) = SubBasisLayout()
-sublayout(::BasisLayout, ::Type{<:Tuple{<:AbstractAffineQuasiVector,<:AbstractUnitRange}}) = MappedBasisLayout()
+sublayout(::AbstractBasisLayout, ::Type{<:Tuple{<:AbstractAffineQuasiVector,<:AbstractUnitRange}}) = MappedBasisLayout()
 
 @inline sub_materialize(::AbstractBasisLayout, V::AbstractQuasiArray) = V
 @inline sub_materialize(::AbstractBasisLayout, V::AbstractArray) = V

test/runtests.jl

@@ -1,5 +1,6 @@
-using ContinuumArrays, QuasiArrays, LazyArrays, IntervalSets, FillArrays, LinearAlgebra, BandedMatrices, Test
-import ContinuumArrays: ℵ₁, materialize, SimplifyStyle, AffineQuasiVector, BasisLayout, AdjointBasisLayout, SubBasisLayout, MappedBasisLayout, igetindex
+using ContinuumArrays, QuasiArrays, LazyArrays, IntervalSets, FillArrays, LinearAlgebra, BandedMatrices, FastTransforms, Test
+import ContinuumArrays: ℵ₁, materialize, SimplifyStyle, AffineQuasiVector, BasisLayout, AdjointBasisLayout, SubBasisLayout, 
+                        MappedBasisLayout, igetindex, TransformFactorization, Weight, WeightedBasisLayout
 import QuasiArrays: SubQuasiArray, MulQuasiMatrix, Vec, Inclusion, QuasiDiagonal, LazyQuasiArrayApplyStyle, LazyQuasiArrayStyle
 import LazyArrays: MemoryLayout, ApplyStyle, Applied, colsupport, arguments, ApplyLayout, LdivApplyStyle
 
@@ -228,6 +229,7 @@ end
         L = LinearSpline([1,2,3])
         x = axes(L,1)
         @test (L \ x) == [1,2,3]
+        @test factorize(L[:,2:end-1]) isa ContinuumArrays.ProjectionFactorization
 
         L = LinearSpline(range(0,1; length=10_000))
         x = axes(L,1)
@@ -347,6 +349,8 @@ end
     y = 2x .- 1
     L = LinearSpline(range(-1,stop=1,length=10))
     @test L[y,:][0.1,:] == L[2*0.1-1,:]
+    @test L[y,:]'L[y,:] isa SymTridiagonal
+    @test L[y,:]'L[y,:] == 1/2*(L'L)
 
     D = Derivative(axes(L,1))
     H = HeavisideSpline(L.points)
@@ -368,6 +372,7 @@ end
     @test B[0.1,1] == L[2*0.1-1,2]
     @test B\B == Eye(8)
     @test L[y,:] \ B == Eye(10)[:,2:end-1]
+    @test B'B == (1/2)*(L'L)[2:end-1,2:end-1]
 end
 
 @testset "diff" begin
@@ -382,4 +387,33 @@ end
     K = x .- x'
     @test K[0.1,0.2] == K[Inclusion(0..0.5), Inclusion(0..0.5)][0.1,0.2] == 0.1 - 0.2
     @test_throws BoundsError K[Inclusion(0..0.5), Inclusion(0..0.5)][1,1]
+end
+
+struct Chebyshev <: Basis{Float64}
+    n::Int
+end
+
+struct ChebyshevWeight <: Weight{Float64} end
+
+
+Base.axes(T::Chebyshev) = (Inclusion(-1..1), Base.OneTo(T.n))
+ContinuumArrays.grid(T::Chebyshev) = chebyshevpoints(Float64, T.n; kind=1)
+Base.axes(T::ChebyshevWeight) = (Inclusion(-1..1),)
+
+LinearAlgebra.factorize(L::Chebyshev) =
+    TransformFactorization(grid(L), plan_chebyshevtransform(Array{Float64}(undef, size(L,2))))
+
+@testset "Chebyshev" begin
+    T = Chebyshev(5)
+    x = axes(T,1)
+    F = factorize(T)
+    g = grid(F)
+    @test T \ exp.(x) == F \ exp.(x) == F \ exp.(g) == chebyshevtransform(exp.(g); kind=1)
+
+    w = ChebyshevWeight()
+    wT = w .* T
+    wT2 = w .* T[:,2:4]
+    @test MemoryLayout(typeof(wT)) == WeightedBasisLayout()
+    @test MemoryLayout(typeof(wT2)) == WeightedBasisLayout()
+    @test grid(wT) == grid(wT2) == grid(T)
 end