diff --git a/Project.toml b/Project.toml index 0f6ebfe..f437e67 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" diff --git a/src/ClassicalOrthogonalPolynomials.jl b/src/ClassicalOrthogonalPolynomials.jl index 6f7f0db..6a7db81 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/adaptivetransform.jl b/src/adaptivetransform.jl index a623ad5..68d7836 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 8eaa23c..022e3a2 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 745daac..ec4154f 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), 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/src/roots.jl b/src/roots.jl index 5f99bc2..7d0210c 100644 --- a/src/roots.jl +++ b/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...) \ 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 + +#### +# 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 479baf9..58c9a5a 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,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 \ No newline at end of file