Add WeightedFactorization (#81)

* Add WeightedFactorization

* Update runtests.jl

* to_indices for BlockIndex arrays

* increase coverage

* Update runtests.jl

Project.toml

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

src/ContinuumArrays.jl

@@ -74,6 +74,9 @@ end
 @inline to_indices(A::AbstractQuasiArray, inds, I::Tuple{BlockIndex{1}, Vararg{Any}}) =
     (inds[1][I[1]], to_indices(A, _maybetail(inds), tail(I))...)
 
+@inline to_indices(A::AbstractQuasiArray, inds, I::Tuple{AbstractArray{<:BlockIndex{1}}, Vararg{Any}}) =
+    (inds[1][I[1]], to_indices(A, _maybetail(inds), tail(I))...)    
+
 checkpoints(d::AbstractInterval{T}) where T = width(d) .* SVector{3,float(T)}(0.823972,0.01,0.3273484) .+ leftendpoint(d)
 checkpoints(x::Inclusion) = checkpoints(x.domain)
 checkpoints(A::AbstractQuasiMatrix) = checkpoints(axes(A,1))

src/bases/bases.jl

@@ -40,6 +40,8 @@ sublayout(::AbstractBasisLayout, ::Type{<:Tuple{Map,AbstractUnitRange}}) = Mappe
 ## Weighted basis interface
 unweightedbasis(P::BroadcastQuasiMatrix{<:Any,typeof(*),<:Tuple{AbstractQuasiVector,AbstractQuasiMatrix}}) = last(P.args)
 unweightedbasis(V::SubQuasiArray) = view(unweightedbasis(parent(V)), parentindices(V)...)
+weight(P::BroadcastQuasiMatrix{<:Any,typeof(*),<:Tuple{AbstractQuasiVector,AbstractQuasiMatrix}}) = first(P.args)
+weight(V::SubQuasiArray) = weight(parent(V))[parentindices(V)[1]]
 
 
 
@@ -181,6 +183,15 @@ copy(L::Ldiv{<:AbstractBasisLayout,<:Any,<:Any,<:AbstractQuasiVector}) =
 copy(L::Ldiv{<:AbstractBasisLayout,ApplyLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) =
     transform_ldiv(L.A, L.B)
 
+struct WeightedFactorization{T, WW, FAC<:Factorization{T}} <: Factorization{T}
+    w::WW
+    F::FAC
+end
+
+_factorize(::WeightedBasisLayouts, wS) = WeightedFactorization(weight(wS), factorize(unweightedbasis(wS)))
+
+
+\(F::WeightedFactorization, b::AbstractQuasiVector) = F.F \ (b ./ F.w)
 
 ##
 # Algebra

test/runtests.jl

@@ -1,7 +1,7 @@
 using ContinuumArrays, QuasiArrays, LazyArrays, IntervalSets, FillArrays, LinearAlgebra, BandedMatrices, FastTransforms, InfiniteArrays, Test, Base64
 import ContinuumArrays: ℵ₁, materialize, AffineQuasiVector, BasisLayout, AdjointBasisLayout, SubBasisLayout, ℵ₁,
                         MappedBasisLayout, AdjointMappedBasisLayout, MappedWeightedBasisLayout, TransformFactorization, Weight, WeightedBasisLayout, SubWeightedBasisLayout, WeightLayout,
-                        Expansion, basis, invmap, Map
+                        Expansion, basis, invmap, Map, checkpoints
 import QuasiArrays: SubQuasiArray, MulQuasiMatrix, Vec, Inclusion, QuasiDiagonal, LazyQuasiArrayApplyStyle, LazyQuasiArrayStyle
 import LazyArrays: MemoryLayout, ApplyStyle, Applied, colsupport, arguments, ApplyLayout, LdivStyle, MulStyle
 
@@ -425,6 +425,7 @@ 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))))
@@ -451,9 +452,13 @@ ContinuumArrays.invmap(::InvQuadraticMap{T}) where T = QuadraticMap{T}()
     w = ChebyshevWeight()
     wT = w .* T
     x = axes(T,1)
-    F = factorize(T)
-    g = grid(F)
-    @test T \ exp.(x) == F \ exp.(x) == F \ exp.(g) == chebyshevtransform(exp.(g), Val(1))
+
+    @testset "basics" begin
+        F = factorize(T)
+        g = grid(F)
+        @test T \ exp.(x) == F \ exp.(x) == F \ exp.(g) == chebyshevtransform(exp.(g), Val(1))
+        @test all(checkpoints(T) .∈ Ref(axes(T,1)))
+    end
 
     @testset "Weighted" begin
         @test MemoryLayout(w) isa WeightLayout
@@ -469,6 +474,10 @@ ContinuumArrays.invmap(::InvQuadraticMap{T}) where T = QuadraticMap{T}()
         @test ContinuumArrays.unweightedbasis(wT) ≡ T
         @test ContinuumArrays.unweightedbasis(wT2) ≡ T[:,2:4]
         @test ContinuumArrays.unweightedbasis(wT3) ≡ T[:,2:4]
+
+        @test ContinuumArrays.weight(wT) ≡ ContinuumArrays.weight(wT2) ≡ ContinuumArrays.weight(wT3) ≡ w
+
+        @test wT \ @.(exp(x) / sqrt(1-x^2)) ≈ T \ exp.(x)
     end
     @testset "Mapped" begin
         y = affine(0..1, x)
@@ -508,6 +517,7 @@ ContinuumArrays.invmap(::InvQuadraticMap{T}) where T = QuadraticMap{T}()
 
     @testset "Broadcasted" begin
         T = Chebyshev(5)
+        F = factorize(T)
         x = axes(T,1)
         a = 1 .+ x .+ x.^2
         # The following are wrong, just testing dispatch

test/test_basisconcat.jl

@@ -61,12 +61,15 @@ import ContinuumArrays: PiecewiseBasis, VcatBasis, HvcatBasis
         S1 = LinearSpline(0:1)
         S2 = LinearSpline(0:0.5:1)
         S = HvcatBasis(2, S1, S2, S2, S1)
+        D = Derivative(axes(S,1))
 
         @test S == S
 
         @test S[0.1, 1] == [S1[0.1,1] 0; 0 0]
         @test S[0.1,Block(1)] == [[S1[0.1,1] 0; 0 0], [S1[0.1,2] 0; 0 0]]
+        @test S[0.1,Block(1)[1]] == [S1[0.1,1] 0; 0 0]
+        @test S[0.1,getindex.(Block(1),1:2)] == [[S1[0.1,1] 0; 0 0], [S1[0.1,2] 0; 0 0]]
         D = Derivative(axes(S,1))
-        @test_broken (D*S)[0.1,1]
+        @test_broken (D*S)[0.1,1] # throws error
     end
 end