sampling, extrema (#242) and preliminary support orthogonalising quasimatrices

Project.toml

@@ -47,9 +47,9 @@ FastTransforms = "0.16.6, 0.17"
 FillArrays = "1"
 HypergeometricFunctions = "0.3.4"
 InfiniteArrays = " 0.15"
-InfiniteLinearAlgebra = "0.10"
+InfiniteLinearAlgebra = "0.10.1"
 IntervalSets = "0.7"
-LazyArrays = "2.5.1"
+LazyArrays = "2.7"
 LazyBandedMatrices = "0.11"
 MutableArithmetics = "1"
 QuasiArrays = "0.12"

src/ClassicalOrthogonalPolynomials.jl

@@ -11,7 +11,7 @@ import Base: @_inline_meta, axes, getindex, unsafe_getindex, convert, prod, *, /
                 IndexStyle, IndexLinear, ==, OneTo, tail, similar, copyto!, copy, setindex,
                 first, last, Slice, size, length, axes, IdentityUnitRange, sum, _sum, cumsum,
                 to_indices, tail, getproperty, inv, show, isapprox, summary,
-                findall, searchsortedfirst, diff
+                findall, searchsortedfirst, diff, minimum, maximum, extrema
 import Base.Broadcast: materialize, BroadcastStyle, broadcasted, Broadcasted
 import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, adjointlayout,
                 sub_materialize, arguments, sub_paddeddata, paddeddata, AbstractPaddedLayout, PaddedColumns, resizedata!, LazyVector, ApplyLayout, call,

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), 1), 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,27 @@ 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
+
+####
+# min/max/extrema
+####
+function minimum(f::AbstractQuasiVector)
+    r = findall(iszero, diff(f))
+    min(first(f), minimum(f[r]), last(f))
+end
+
+function maximum(f::AbstractQuasiVector)
+    r = findall(iszero, diff(f))
+    max(first(f), maximum(f[r]), last(f))
+end
+
+function extrema(f::AbstractQuasiVector)
+    r = findall(iszero, diff(f))
+    extrema([first(f); f[r]; last(f)])
+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,24 @@ 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 =  cos.((0:5)' .* x)
+    @test (Legendre() \ A)[1:5,1] ≈ [1; zeros(4)] # test transform bug
+    @test (P \ A)[1:5,1] ≈ [sqrt(2); zeros(4)]
+    Q,R = qr(P \ A)
+end
+
+@testset "minimum/maximum/extrema (#242)" begin
+    f = expand(ChebyshevT(), x -> exp(x) * cos(100x.^2))
+    @test minimum(f) ≈ -2.682833127491678
+    @test maximum(f) ≈ 2.6401248792053362
+    @test extrema(f) == (minimum(f), maximum(f))
 end