From 525e4128cb6fac22314fec1310deb72a3ef7515c Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 29 Sep 2025 03:50:42 +0100 Subject: [PATCH 1/4] sampling and orthogonalising quasimatrices --- src/adaptivetransform.jl | 8 +++++--- src/classical/jacobi.jl | 2 +- src/classical/legendre.jl | 2 +- src/roots.jl | 7 ++++++- test/test_roots.jl | 16 +++++++++++++++- 5 files changed, 28 insertions(+), 7 deletions(-) diff --git a/src/adaptivetransform.jl b/src/adaptivetransform.jl index a623ad5d..68d7836f 100644 --- a/src/adaptivetransform.jl +++ b/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? diff --git a/src/classical/jacobi.jl b/src/classical/jacobi.jl index 8eaa23c0..022e3a28 100644 --- a/src/classical/jacobi.jl +++ b/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 diff --git a/src/classical/legendre.jl b/src/classical/legendre.jl index 745daace..396b5caa 100644 --- a/src/classical/legendre.jl +++ b/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...)) diff --git a/src/roots.jl b/src/roots.jl index 5f99bc28..d9df2717 100644 --- a/src/roots.jl +++ b/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...) \ No newline at end of file +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 \ No newline at end of file diff --git a/test/test_roots.jl b/test/test_roots.jl index 479baf94..1bf33a25 100644 --- a/test/test_roots.jl +++ b/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 \ No newline at end of file From cf832d094c91a70148799a3fbf9dea5982fd7016 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Fri, 3 Oct 2025 22:38:59 +0100 Subject: [PATCH 2/4] Fix #242 --- src/ClassicalOrthogonalPolynomials.jl | 2 +- src/roots.jl | 18 ++++++++++++++++++ test/test_roots.jl | 7 +++++++ 3 files changed, 26 insertions(+), 1 deletion(-) diff --git a/src/ClassicalOrthogonalPolynomials.jl b/src/ClassicalOrthogonalPolynomials.jl index 6f7f0db2..6a7db811 100644 --- a/src/ClassicalOrthogonalPolynomials.jl +++ b/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, diff --git a/src/roots.jl b/src/roots.jl index d9df2717..7d0210c9 100644 --- a/src/roots.jl +++ b/src/roots.jl @@ -54,4 +54,22 @@ searchsortedfirst(f::AbstractQuasiVector, x; kwds...) = _searchsortedfirst(Memor 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 \ No newline at end of file diff --git a/test/test_roots.jl b/test/test_roots.jl index 1bf33a25..aec7b1de 100644 --- a/test/test_roots.jl +++ b/test/test_roots.jl @@ -22,4 +22,11 @@ end P = Normalized(Legendre()); x = axes(P,1) A = [(1 .+ x) (1 .+ x.^3)] 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 \ No newline at end of file From d885d091460c35aca79ff87501860fa84b8fa16e Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sat, 4 Oct 2025 14:29:50 +0100 Subject: [PATCH 3/4] Update Project.toml --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 0f6ebfe2..f437e676 100644 --- a/Project.toml +++ b/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" From c8cc5e7372a4cc2ad5c5b22e2818e22693fa5d80 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sat, 4 Oct 2025 21:26:35 +0100 Subject: [PATCH 4/4] Fixlegendre transform on matrices --- src/classical/legendre.jl | 2 +- test/test_roots.jl | 4 +++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/classical/legendre.jl b/src/classical/legendre.jl index 396b5caa..ec4154f3 100644 --- a/src/classical/legendre.jl +++ b/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) - padrows(th_cheb2leg(paddeddata(dat)), axes(dat,1)) + 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...)) diff --git a/test/test_roots.jl b/test/test_roots.jl index aec7b1de..58c9a5ac 100644 --- a/test/test_roots.jl +++ b/test/test_roots.jl @@ -20,7 +20,9 @@ end @testset "det point sampling" begin P = Normalized(Legendre()); x = axes(P,1) - A = [(1 .+ x) (1 .+ x.^3)] + 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