sampling and orthogonalising quasimatrices

Author: dlfivefifty

src/adaptivetransform.jl

@@ -7,6 +7,9 @@ ContinuumArrays.transform_ldiv_size(::Tuple{<:Any,InfiniteCardinal{0}}, A::Abstr
 padchop!(cfs, tol, ax...) = pad(chop!(cfs, tol), ax...)
 padchop(cfs, tol, ax...) = pad(chop(cfs, tol), ax...)
 
+padrows(a::AbstractVector, ax) = pad(a, ax)
+padrows(a::AbstractMatrix, ax) = pad(a, ax, :)
+
 # ax will impose block structure for us
 padchop!(cfs::BlockedVector, tol, ax...) = padchop!(cfs.blocks, tol, ax...)
 
@@ -56,7 +59,6 @@ end
 function adaptivetransform_ldiv(A::AbstractQuasiArray{U}, f::AbstractQuasiMatrix{V}) where {U,V}
     T = promote_type(eltype(U),eltype(V))
 
-    bx = axes(f,2)
     r = checkpoints(A)
     fr = f[r,:]
     maxabsfr = norm(fr,Inf)
@@ -69,12 +71,12 @@ function adaptivetransform_ldiv(A::AbstractQuasiArray{U}, f::AbstractQuasiMatrix
         cfs = An \ f
         maxabsc = maximum(abs, cfs)
         if maxabsc == 0 && maxabsfr == 0
-            return pad(similar(cfs,0,size(cfs,2)), ax, bx)
+            return pad(similar(cfs,0,size(cfs,2)), ax, :)
         end
 
         if maximum(abs,@views(cfs[end-2:end,:])) < 10tol*maxabsc
             n = size(cfs,1)
-            c = padchop(cfs, tol, ax, bx)
+            c = padchop(cfs, tol, ax, :)
             un = A * c # expansion
         # we allow for transformed coefficients being a different size
         ##TODO: how to do scaling for unnormalized bases like Jacobi?

src/classical/jacobi.jl

@@ -288,7 +288,7 @@ ldiv(P::Jacobi{V}, f::Inclusion{T}) where {T,V} = _op_ldiv(P, f)
 ldiv(P::Jacobi{V}, f::AbstractQuasiFill{T,1}) where {T,V} = _op_ldiv(P, f)
 function transform_ldiv(P::Jacobi{V}, f::AbstractQuasiArray) where V
     T = ChebyshevT{V}()
-    pad(FastTransforms.th_cheb2jac(paddeddata(T \ f), P.a, P.b, 1), axes(P,2), tail(axes(f))...)
+    padrows(FastTransforms.th_cheb2jac(paddeddata(T \ f), P.a, P.b, 1), axes(P,2))
 end
 
 

src/classical/legendre.jl

@@ -189,7 +189,7 @@ ldiv(P::Legendre{V}, f::AbstractQuasiFill{T,1}) where {T,V} = _op_ldiv(P, f)
 function transform_ldiv(::Legendre{V}, f::Union{AbstractQuasiVector,AbstractQuasiMatrix}) where V
     T = ChebyshevT{V}()
     dat = transform_ldiv(T, f)
-    pad(th_cheb2leg(paddeddata(dat)), axes(dat)...)
+    padrows(th_cheb2leg(paddeddata(dat)), axes(dat,1))
 end
 
 plan_transform(::Legendre{T}, szs::NTuple{N,Int}, dims...) where {T,N} = JacobiTransformPlan(FastTransforms.plan_th_cheb2leg!(T, szs, dims...), plan_chebyshevtransform(T, szs, dims...))

src/roots.jl

@@ -49,4 +49,9 @@ function _searchsortedfirst(::ExpansionLayout{<:AbstractOPLayout}, f, x; iterati
     end
     (a+b)/2
 end
-searchsortedfirst(f::AbstractQuasiVector, x; kwds...) = _searchsortedfirst(MemoryLayout(f), f, x; kwds...)
+searchsortedfirst(f::AbstractQuasiVector, x; kwds...) = _searchsortedfirst(MemoryLayout(f), f, x; kwds...)
+
+function sample(f::AbstractQuasiVector, n...)
+    g = cumsum(f)
+    searchsortedfirst.(Ref(g/last(g)), rand(n...))
+end

test/test_roots.jl

@@ -1,4 +1,7 @@
-using ClassicalOrthogonalPolynomials, Test
+using ClassicalOrthogonalPolynomials, Random, Test
+using ClassicalOrthogonalPolynomials: sample
+
+Random.seed!(5)
 
 @testset "roots" begin
     T = Chebyshev()
@@ -8,4 +11,15 @@ using ClassicalOrthogonalPolynomials, Test
 
     g = x -> x + 0.001cos(x)
     @test searchsortedfirst(expand(T, g), 0.1) ≈ searchsortedfirst(expand(P, g), 0.1) ≈ findall(iszero, expand(T, x -> g(x)-0.1))[1]
+end
+
+@testset "sample" begin
+    f = expand(Chebyshev(), exp)
+    @test sum(sample(f, 100_000))/100_000 ≈ 0.31 atol=1E-2
+end
+
+@testset "det point sampling" begin
+    P = Normalized(Legendre()); x = axes(P,1)
+    A =  [(1 .+ x) (1 .+ x.^3)]
+    Q,R = qr(P \ A)
 end