Giter VIP home page Giter VIP logo

orthogonalpolynomialsquasi.jl's Introduction

OrthogonalPolynomialsQuasi.jl

A package for representing orthogonal polynomials as quasi arrays

Build Status codecov

This package has been superseded by ClassicalOrthogonalPolynomials.jl

This package implements classical orthogonal polynomials as quasi-arrays where one one axes is continuous and the other axis is discrete (countably infinite), as implemented in QuasiArrays.jl and ContinuumArrays.jl.

julia> using OrthogonalPolynomialsQuasi, ContinuumArrays

julia> P = Legendre(); # Legendre polynomials

julia> size(P) # uncountable ∞ x countable ∞
(ℵ₁, ∞)

julia> axes(P) # essentially (-1..1, 1:∞), Inclusion plays the same role as Slice
(Inclusion(-1.0..1.0 (Chebyshev)), OneToInf())

julia> P[0.1,1:10] # [P_0(0.1), …, P_9(0.1)]
10-element Array{Float64,1}:
  1.0                
  0.1                
 -0.485              
 -0.14750000000000002
  0.3379375          
  0.17882875         
 -0.2488293125       
 -0.19949294375000004
  0.180320721484375  
  0.21138764183593753

julia> @time P[range(-1,1; length=10_000), 1:10_000]; # construct 10_000^2 Vandermonde matrix
  1.624796 seconds (10.02 k allocations: 1.491 GiB, 6.81% gc time)

This also works for associated Legendre polynomials as weighted Ultraspherical polynomials:

julia> associatedlegendre(m) = ((-1)^m*prod(1:2:(2m-1)))*(UltrasphericalWeight((m+1)/2).*Ultraspherical(m+1/2))
associatedlegendre (generic function with 1 method)

julia> associatedlegendre(2)[0.1,1:10]
10-element Array{Float64,1}:
   2.9699999999999998
   1.4849999999999999
  -6.9052500000000006
  -5.041575          
  10.697754375       
  10.8479361375      
 -13.334647528125    
 -18.735466024687497 
  13.885467170308594 
  28.220563705988674 

p-Finite Element Method

The language of quasi-arrays gives a natural framework for constructing p-finite element methods. The convention is that adjoint-products are understood as inner products over the axes with uniform weight. Thus to solve Poisson's equation using its weak formulation with Dirichlet conditions we can expand in a weighted Jacobi basis:

julia> P¹¹ = Jacobi(1.0,1.0); # Quasi-matrix of Jacobi polynomials

julia> w = JacobiWeight(1.0,1.0); # quasi-vector correspoinding to (1-x^2)

julia> w[0.1]  (1-0.1^2)
true

julia> S = w .* P¹¹; # Quasi-matrix of weighted Jacobi polynomials

julia> D = Derivative(axes(S,1)); # quasi-matrix corresponding to derivative

julia> Δ = (D*S)'*(D*S) # weak laplacian corresponding to inner products of weighted Jacobi polynomials×∞ LazyArrays.ApplyArray{Float64,2,typeof(*),Tuple{Adjoint{Int64,BandedMatrices.BandedMatrix{Int64,Adjoint{Int64,InfiniteArrays.InfStepRange{Int64,Int64}},InfiniteArrays.OneToInf{Int64}}},LazyArrays.BroadcastArray{Float64,2,typeof(*),Tuple{LazyArrays.BroadcastArray{Float64,1,typeof(/),Tuple{Int64,InfiniteArrays.InfStepRange{Int64,Int64}}},BandedMatrices.BandedMatrix{Int64,Adjoint{Int64,InfiniteArrays.InfStepRange{Int64,Int64}},InfiniteArrays.OneToInf{Int64}}}}}} with indices OneToInf()×OneToInf():
 2.66667                                                        
         6.4                                                     
             10.2857                                             
                     14.2222                                     
                             18.1818                             
                                     22.1538                    
                                             26.1333             
                                                     30.1176     
                                                                
                                                                
                                                               
                                                                
                                                                     

orthogonalpolynomialsquasi.jl's People

Contributors

dlfivefifty avatar github-actions[bot] avatar juliatagbot avatar mikaelslevinsky avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

orthogonalpolynomialsquasi.jl's Issues

operator==operator throws an error

julia> M=(Inclusion(-1..1) .* Jacobi(1,1)).args[2] # multiplication operator
∞×∞ BandedMatrices.BandedMatrix{Float64,LazyArrays.ApplyArray{Float64,2,typeof(vcat),Tuple{Adjoint{Float64,LazyArrays.BroadcastArray{Float64,1,typeof(/),Tuple{LazyArrays.BroadcastArray{Float64,1,typeof(*),Tuple{InfiniteArrays.InfStepRange{Float64,Float64},InfiniteArrays.InfStepRange{Float64,Float64}}},LazyArrays.BroadcastArray{Float64,1,typeof(*),Tuple{InfiniteArrays.InfStepRange{Float64,Float64},InfiniteArrays.InfStepRange{Float64,Float64}}}}}},LazyArrays.ApplyArray{Float64,2,typeof(hcat),Tuple{Float64,Adjoint{Float64,LazyArrays.BroadcastArray{Float64,1,typeof(/),Tuple{Float64,LazyArrays.BroadcastArray{Float64,1,typeof(*),Tuple{InfiniteArrays.InfStepRange{Float64,Float64},InfiniteArrays.InfStepRange{Float64,Float64}}}}}}}},LazyArrays.ApplyArray{Float64,2,typeof(hcat),Tuple{Float64,Adjoint{Float64,LazyArrays.BroadcastArray{Float64,1,typeof(/),Tuple{LazyArrays.BroadcastArray{Float64,1,typeof(*),Tuple{InfiniteArrays.InfStepRange{Int64,Int64},InfiniteArrays.InfStepRange{Float64,Float64}}},LazyArrays.BroadcastArray{Float64,1,typeof(*),Tuple{InfiniteArrays.InfStepRange{Float64,Float64},InfiniteArrays.InfStepRange{Float64,Float64}}}}}}}}}},InfiniteArrays.OneToInf{Int64}} with indices OneToInf()×OneToInf():
 0.0  0.4                                                                            
 0.5  0.0       0.428571                                                                 
     0.533333  0.0       0.444444                                                  
              0.535714  0.0       0.454545                                         
                       0.533333  0.0       0.461538                                
                                0.530303  0.0       0.466667                          
                                                                                             

julia> M==M
ERROR: TypeError: in new, expected Tuple{Int64,Int64}, got a value of type Tuple{InfiniteArrays.Infinity,InfiniteArrays.Infinity}
Stacktrace:
 [1] CartesianIndex{2}(::Tuple{InfiniteArrays.Infinity,InfiniteArrays.Infinity}) at .\multidimensional.jl:65
 [2] CartesianIndex(::Tuple{InfiniteArrays.Infinity,InfiniteArrays.Infinity}) at .\multidimensional.jl:68
 [3] last(::CartesianIndices{2,Tuple{InfiniteArrays.OneToInf{Int64},InfiniteArrays.OneToInf{Int64}}}) at .\multidimensional.jl:389
 [4] iterate at .\multidimensional.jl:347 [inlined]
 [5] iterate(::BandedMatrices.BandedMatrix{Float64,LazyArrays.ApplyArray{Float64,2,typeof(vcat),Tuple{Adjoint{Float64,LazyArrays.BroadcastArray{Float64,1,typeof(/),Tuple{LazyArrays.BroadcastArray{Float64,1,typeof(*),Tuple{InfiniteArrays.InfStepRange{Float64,Float64},InfiniteArrays.InfStepRange{Float64,Float64}}},LazyArrays.BroadcastArray{Float64,1,typeof(*),Tuple{InfiniteArrays.InfStepRange{Float64,Float64},InfiniteArrays.InfStepRange{Float64,Float64}}}}}},LazyArrays.ApplyArray{Float64,2,typeof(hcat),Tuple{Float64,Adjoint{Float64,LazyArrays.BroadcastArray{Float64,1,typeof(/),Tuple{Float64,LazyArrays.BroadcastArray{Float64,1,typeof(*),Tuple{InfiniteArrays.InfStepRange{Float64,Float64},InfiniteArrays.InfStepRange{Float64,Float64}}}}}}}},LazyArrays.ApplyArray{Float64,2,typeof(hcat),Tuple{Float64,Adjoint{Float64,LazyArrays.BroadcastArray{Float64,1,typeof(/),Tuple{LazyArrays.BroadcastArray{Float64,1,typeof(*),Tuple{InfiniteArrays.InfStepRange{Int64,Int64},InfiniteArrays.InfStepRange{Float64,Float64}}},LazyArrays.BroadcastArray{Float64,1,typeof(*),Tuple{InfiniteArrays.InfStepRange{Float64,Float64},InfiniteArrays.InfStepRange{Float64,Float64}}}}}}}}}},InfiniteArrays.OneToInf{Int64}}, ::Tuple{CartesianIndices{2,Tuple{InfiniteArrays.OneToInf{Int64},InfiniteArrays.OneToInf{Int64}}}}) at .\abstractarray.jl:984 (repeats 2 times)
 [6] _zip_iterate_some at .\iterators.jl:352 [inlined]
 [7] _zip_iterate_all at .\iterators.jl:344 [inlined]
 [8] iterate at .\iterators.jl:334 [inlined]
 [9] ==(::BandedMatrices.BandedMatrix{Float64,LazyArrays.ApplyArray{Float64,2,typeof(vcat),Tuple{Adjoint{Float64,LazyArrays.BroadcastArray{Float64,1,typeof(/),Tuple{LazyArrays.BroadcastArray{Float64,1,typeof(*),Tuple{InfiniteArrays.InfStepRange{Float64,Float64},InfiniteArrays.InfStepRange{Float64,Float64}}},LazyArrays.BroadcastArray{Float64,1,typeof(*),Tuple{InfiniteArrays.InfStepRange{Float64,Float64},InfiniteArrays.InfStepRange{Float64,Float64}}}}}},LazyArrays.ApplyArray{Float64,2,typeof(hcat),Tuple{Float64,Adjoint{Float64,LazyArrays.BroadcastArray{Float64,1,typeof(/),Tuple{Float64,LazyArrays.BroadcastArray{Float64,1,typeof(*),Tuple{InfiniteArrays.InfStepRange{Float64,Float64},InfiniteArrays.InfStepRange{Float64,Float64}}}}}}}},LazyArrays.ApplyArray{Float64,2,typeof(hcat),Tuple{Float64,Adjoint{Float64,LazyArrays.BroadcastArray{Float64,1,typeof(/),Tuple{LazyArrays.BroadcastArray{Float64,1,typeof(*),Tuple{InfiniteArrays.InfStepRange{Int64,Int64},InfiniteArrays.InfStepRange{Float64,Float64}}},LazyArrays.BroadcastArray{Float64,1,typeof(*),Tuple{InfiniteArrays.InfStepRange{Float64,Float64},InfiniteArrays.InfStepRange{Float64,Float64}}}}}}}}}},InfiniteArrays.OneToInf{Int64}}, ::BandedMatrices.BandedMatrix{Float64,LazyArrays.ApplyArray{Float64,2,typeof(vcat),Tuple{Adjoint{Float64,LazyArrays.BroadcastArray{Float64,1,typeof(/),Tuple{LazyArrays.BroadcastArray{Float64,1,typeof(*),Tuple{InfiniteArrays.InfStepRange{Float64,Float64},InfiniteArrays.InfStepRange{Float64,Float64}}},LazyArrays.BroadcastArray{Float64,1,typeof(*),Tuple{InfiniteArrays.InfStepRange{Float64,Float64},InfiniteArrays.InfStepRange{Float64,Float64}}}}}},LazyArrays.ApplyArray{Float64,2,typeof(hcat),Tuple{Float64,Adjoint{Float64,LazyArrays.BroadcastArray{Float64,1,typeof(/),Tuple{Float64,LazyArrays.BroadcastArray{Float64,1,typeof(*),Tuple{InfiniteArrays.InfStepRange{Float64,Float64},InfiniteArrays.InfStepRange{Float64,Float64}}}}}}}},LazyArrays.ApplyArray{Float64,2,typeof(hcat),Tuple{Float64,Adjoint{Float64,LazyArrays.BroadcastArray{Float64,1,typeof(/),Tuple{LazyArrays.BroadcastArray{Float64,1,typeof(*),Tuple{InfiniteArrays.InfStepRange{Int64,Int64},InfiniteArrays.InfStepRange{Float64,Float64}}},LazyArrays.BroadcastArray{Float64,1,typeof(*),Tuple{InfiniteArrays.InfStepRange{Float64,Float64},InfiniteArrays.InfStepRange{Float64,Float64}}}}}}}}}},InfiniteArrays.OneToInf{Int64}}) at .\abstractarray.jl:1864
 [10] top-level scope at none:1

The same issue happens with conversion operator, etc. I'm not sure which package causes this exactly so I just put this issue here.

Problem with JacobiWeight?

α = -0.5
β = -0.5
w = JacobiWeight(α,β)
p = LanczosPolynomial(w)
x = axes(w,1)
b = 2
ϕw = JacobiWeight(α+1,β+1).*(b .- x)
pϕ = LanczosPolynomial(ϕw)

The following doesn't give a diagonal operator:
pϕ' * (ϕw .* pϕ)

p' * (w .* p)
gives
not implemented for OrthogonalPolynomialsQuasi.Jacobi{Float64}(0.0, 0.0) and OrthogonalPolynomialsQuasi.Jacobi{Float64}(-0.5, -0.5)

Add normalized polynomials

I was going to propose trying to stick to linear algebra language as much as possible, in which case normalisation is really just a QR. That is:

P = Legendre()
Q,D = qr(P)
qr(P) isa OrthogonalPolynomialQR{Float64}
@test Q isa Normalized{Float64,Legendre{Float64}}
@test D isa Diagonal
x,n = 0.1,5
@test P[x,n] == (Q*D)[x,n]

For other OPs we would need to weight them first, to avoid dealing with inner products:

T = Chebyshev()
wT = SqrtWeighted(T)
@test wT isa SqrtWeighted{Float64,Chebyshev{Float64}}
wQ,D = qr(wT)
@test wQ isa SqrtWeighted{Float64,Normalized{Float64,Chebyshev{Float64}}}
@test wQ[x,n] == (T*D)[x,n]/sqrt(1-x^2)

But this seems to suggest that Normalized{T,<:OrthogonalPolynomial} <: OrthogonalPolynomial is the more fundamental type.

StackOverflow with quadratically shifted Legendre and Jacobi polynomials

I am encountering a StackOverflowError using the latest versions of this package but for some reason this worked in some previous versions (some of them quite old, like it definitely works for v0.1.3)

Consider the following code which I think should allow us to evaluate and expand in a basis like P(2x^2-1). Let's stick to Legendre polynomials for now, but the same happens for Jacobi polynomials.

using ContinuumArrays, OrthogonalPolynomialsQuasi, QuasiArrays, DomainSets, ArrayLayouts, Test, Plots
import ContinuumArrays: AbstractBasisLayout, MappedBasisLayout

struct QuadraticQuasiVector{T} <: AbstractQuasiVector{T} end

function QuadraticQuasiVector()
    QuadraticQuasiVector{Float64}()
end

function Base.getindex(A::QuadraticQuasiVector, k::Inclusion)
    @assert k.domain == UnitInterval()
    A
end

function Base.getindex(::QuadraticQuasiVector, r::Number)
    2r^2-1
end

function Base.axes(::QuadraticQuasiVector{T}) where T
    (Inclusion(UnitInterval()),)
end

function ContinuumArrays.igetindex(d::QuadraticQuasiVector, x::Number)
    sqrt((x+1)/2)
end

function Base.checkindex(::Type{Bool}, inds::Inclusion{<:Any,<:ChebyshevInterval}, r::QuadraticQuasiVector)
    true
end

function ContinuumArrays.sublayout(::AbstractBasisLayout, ::Type{<:Tuple{<:QuadraticQuasiVector,<:AbstractUnitRange}})
    MappedBasisLayout()
end

and we can then run:

julia> P = Legendre()[QuadraticQuasiVector(),:]
view(Legendre{Float64}, QuadraticQuasiVector{Float64}, :) with eltype Float64

julia> x = axes(P,1)
Inclusion(0.0..1.0 (Unit))

Unfortunately in recent versions we then get the following:

julia> P \ (2*x.^2 .-1)
ERROR: LoadError: StackOverflowError:
Stacktrace:
 [1] objectid at ./reflection.jl:312 [inlined]
 [2] dataids at /home/timon/.julia/packages/QuasiArrays/iY9M0/src/abstractquasiarray.jl:492 [inlined]
 [3] mightalias at /home/timon/.julia/packages/QuasiArrays/iY9M0/src/abstractquasiarray.jl:480 [inlined]
 [4] unalias at /home/timon/.julia/packages/QuasiArrays/iY9M0/src/abstractquasiarray.jl:469 [inlined]
 [5] #18 at /home/timon/.julia/packages/QuasiArrays/iY9M0/src/subquasiarray.jl:67 [inlined]
 [6] map at ./tuple.jl:158 [inlined]
 [7] view at /home/timon/.julia/packages/QuasiArrays/iY9M0/src/subquasiarray.jl:67 [inlined]
 [8] layout_getindex at /home/timon/.julia/packages/ArrayLayouts/9PQtL/src/ArrayLayouts.jl:108 [inlined]
 [9] _getindex at /home/timon/.julia/packages/QuasiArrays/iY9M0/src/abstractquasiarray.jl:386 [inlined]
 [10] _getindex at /home/timon/.julia/packages/QuasiArrays/iY9M0/src/abstractquasiarray.jl:377 [inlined]
 [11] getindex at /home/timon/.julia/packages/QuasiArrays/iY9M0/src/abstractquasiarray.jl:372 [inlined]
 [12] demap(::SubQuasiArray{Float64,2,Legendre{Float64},Tuple{QuadraticQuasiVector{Float64},Base.Slice{InfiniteArrays.OneToInf{Int64}}},false}) at /home/timon/.julia/packages/ContinuumArrays/jbgDD/src/bases/bases.jl:295 (repeats 79983 times)
in expression starting at /home/timon/Documents/Projects/TSG_PhD_2/project - equilibrium measures/2 - arbitrary dimensional balls/test/test_quadJacobi.jl:54

julia> P[:,Base.OneTo(10)] \ (2*x.^2 .-1)
ERROR: LoadError: StackOverflowError:
Stacktrace:
 [1] objectid at ./reflection.jl:312 [inlined]
 [2] dataids at /home/timon/.julia/packages/QuasiArrays/iY9M0/src/abstractquasiarray.jl:492 [inlined]
 [3] mightalias at /home/timon/.julia/packages/QuasiArrays/iY9M0/src/abstractquasiarray.jl:480 [inlined]
 [4] unalias at /home/timon/.julia/packages/QuasiArrays/iY9M0/src/abstractquasiarray.jl:469 [inlined]
 [5] #18 at /home/timon/.julia/packages/QuasiArrays/iY9M0/src/subquasiarray.jl:67 [inlined]
 [6] map at ./tuple.jl:158 [inlined]
 [7] view at /home/timon/.julia/packages/QuasiArrays/iY9M0/src/subquasiarray.jl:67 [inlined]
 [8] layout_getindex at /home/timon/.julia/packages/ArrayLayouts/9PQtL/src/ArrayLayouts.jl:108 [inlined]
 [9] _getindex at /home/timon/.julia/packages/QuasiArrays/iY9M0/src/abstractquasiarray.jl:386 [inlined]
 [10] _getindex at /home/timon/.julia/packages/QuasiArrays/iY9M0/src/abstractquasiarray.jl:377 [inlined]
 [11] getindex at /home/timon/.julia/packages/QuasiArrays/iY9M0/src/abstractquasiarray.jl:372 [inlined]
 [12] demap(::SubQuasiArray{Float64,2,Legendre{Float64},Tuple{QuadraticQuasiVector{Float64},Base.Slice{InfiniteArrays.OneToInf{Int64}}},false}) at /home/timon/.julia/packages/ContinuumArrays/jbgDD/src/bases/bases.jl:295 (repeats 79983 times)
in expression starting at /home/timon/Documents/Projects/TSG_PhD_2/project - equilibrium measures/2 - arbitrary dimensional balls/test/test_quadJacobi.jl:52

whereas on the aforementioned version we can do it e.g.:

julia> P[:,Base.OneTo(10)]\(2 *x.^2 .-1)
10-element Array{Float64,1}:
 -1.1102230246251565e-16
  0.9999999999999999
  7.998262167527769e-17
 -5.61254717377705e-17
 -8.822163343216031e-17
  9.47501003660476e-17
  7.407798083866165e-17
 -1.4987585997444753e-16
 -3.4940086295796035e-17
  1.182456355770199e-16

Maybe this has a quick fix somewhere or I'm doing something wrong? Fixes or pointers appreciated.

weight = product of polynomials

x = Inclusion(ChebyshevInterval())
ϕ = x.^3 - 2*x.^2 - x .+ 2
LanczosPolynomial(ϕ)

works, but writing ϕ as ().*() gives

ϕ = (1 .- x.^2).*(2 .- x)
LanczosPolynomial(ϕ)

MethodError: no method matching singularitiesbroadcast(::typeof(*), ::LegendreWeight{Float64}, ::LegendreWeight{Float64})
Closest candidates are:
  singularitiesbroadcast(::typeof(*), ::LegendreWeight, !Matched::OrthogonalPolynomialsQuasi.JacobiWeight) at C:\Users\mfaso\.julia\packages\OrthogonalPolynomialsQuasi\bIORM\src\jacobi.jl:79
  singularitiesbroadcast(::typeof(*), !Matched::OrthogonalPolynomialsQuasi.JacobiWeight, ::LegendreWeight) at C:\Users\mfaso\.julia\packages\OrthogonalPolynomialsQuasi\bIORM\src\jacobi.jl:80
  singularitiesbroadcast(::typeof(*), !Matched::OrthogonalPolynomialsQuasi.NoSingularities, ::LegendreWeight) at In[3]:8
  ...

Stacktrace:
 [1] singularities(::LazyArrays.BroadcastLayout{typeof(*)}, ::BroadcastQuasiArray{Float64,1,typeof(*),Tuple{BroadcastQuasiArray{Float64,1,typeof(-),Tuple{Int64,BroadcastQuasiArray{Float64,1,typeof(Base.literal_pow),Tuple{Base.RefValue{typeof(^)},Inclusion{Float64,ChebyshevInterval{Float64}},Base.RefValue{Val{2}}}}}},ContinuumArrays.AffineQuasiVector{Float64,Float64,Inclusion{Float64,ChebyshevInterval{Float64}},Float64}}}) at C:\Users\mfaso\.julia\packages\OrthogonalPolynomialsQuasi\bIORM\src\OrthogonalPolynomialsQuasi.jl:165
 [2] singularities(::BroadcastQuasiArray{Float64,1,typeof(*),Tuple{BroadcastQuasiArray{Float64,1,typeof(-),Tuple{Int64,BroadcastQuasiArray{Float64,1,typeof(Base.literal_pow),Tuple{Base.RefValue{typeof(^)},Inclusion{Float64,ChebyshevInterval{Float64}},Base.RefValue{Val{2}}}}}},ContinuumArrays.AffineQuasiVector{Float64,Float64,Inclusion{Float64,ChebyshevInterval{Float64}},Float64}}}) at C:\Users\mfaso\.julia\packages\OrthogonalPolynomialsQuasi\bIORM\src\OrthogonalPolynomialsQuasi.jl:167
 [3] LanczosPolynomial(::BroadcastQuasiArray{Float64,1,typeof(*),Tuple{BroadcastQuasiArray{Float64,1,typeof(-),Tuple{Int64,BroadcastQuasiArray{Float64,1,typeof(Base.literal_pow),Tuple{Base.RefValue{typeof(^)},Inclusion{Float64,ChebyshevInterval{Float64}},Base.RefValue{Val{2}}}}}},ContinuumArrays.AffineQuasiVector{Float64,Float64,Inclusion{Float64,ChebyshevInterval{Float64}},Float64}}}) at C:\Users\mfaso\.julia\packages\OrthogonalPolynomialsQuasi\bIORM\src\lanczos.jl:207
 [4] top-level scope at In[27]:5
 [5] include_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1091
 [6] execute_code(::String, ::String) at C:\Users\mfaso\.julia\packages\IJulia\a1SNk\src\execute_request.jl:27
 [7] execute_request(::ZMQ.Socket, ::IJulia.Msg) at C:\Users\mfaso\.julia\packages\IJulia\a1SNk\src\execute_request.jl:86
 [8] #invokelatest#1 at .\essentials.jl:710 [inlined]
 [9] invokelatest at .\essentials.jl:709 [inlined]
 [10] eventloop(::ZMQ.Socket) at C:\Users\mfaso\.julia\packages\IJulia\a1SNk\src\eventloop.jl:8
 [11] (::IJulia.var"#15#18")() at .\task.jl:356

Jacobi(1.0,1.0)[-1,1:∞] throws an error

I tested changing -1 to 0 and the problems still exist.

julia> J=Jacobi(1.0,1.0)
Jacobi{Float64}(1.0, 1.0)

julia> J[-1,:]
∞-element view(::Jacobi{Float64}, -1.0, :) with eltype Float64 with indices OneToInf():
  1.0
 -2.0
  3.0
 -4.0
  5.000000000000001
  

julia> J[-1,1:∞]
ERROR: StackOverflowError:
Stacktrace:
 [1] getindex(::Jacobi{Float64}, ::Int64, ::InfiniteArrays.OneToInf{Int64}) at C:\Users\DELL\.julia\packages\OrthogonalPolynomialsQuasi\llrOw\src\OrthogonalPolynomialsQuasi.jl:234 (repeats 79984 times)

julia> J[-1,:][2:∞] # I waited for some minutes before interrupting Julia
ERROR: InterruptException:
Stacktrace:
 [1] similar(::Array{LazyArrays.LazyLayout,1}, ::Type) at .\array.jl:358
 [2] copyto_nonleaf!(::Array{LazyArrays.LazyLayout,1}, ::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Tuple{Base.OneTo{Int64}},Type{ArrayLayouts.MemoryLayout},Tuple{Base.Broadcast.Extruded{Array{Any,1},Tuple{Bool},Tuple{Int64}}}}, ::Base.OneTo{Int64}, ::Int64, ::Int64) at .\broadcast.jl:1010
 [3] copy at .\broadcast.jl:858 [inlined]
 [4] materialize at .\broadcast.jl:820 [inlined]
 [5] tuple_type_memorylayouts(::Type{Tuple{InfiniteArrays.InfStepRange{Int64,Int64},Float64,Float64,Int64}}) at C:\Users\DELL\.julia\packages\LazyArrays\MAEFm\src\lazybroadcasting.jl:109
 [6] MemoryLayout at C:\Users\DELL\.julia\packages\LazyArrays\MAEFm\src\lazybroadcasting.jl:119 [inlined]
 [7] call(::BroadcastArray{Float64,1,typeof(+),Tuple{InfiniteArrays.InfStepRange{Int64,Int64},Float64,Float64,Int64}}) at C:\Users\DELL\.julia\packages\LazyArrays\MAEFm\src\lazyapplying.jl:22
 [8] Base.Broadcast.Broadcasted(::BroadcastArray{Float64,1,typeof(+),Tuple{InfiniteArrays.InfStepRange{Int64,Int64},Float64,Float64,Int64}}) at C:\Users\DELL\.julia\packages\LazyArrays\MAEFm\src\lazybroadcasting.jl:43
 [9] axes(::BroadcastArray{Float64,1,typeof(+),Tuple{InfiniteArrays.InfStepRange{Int64,Int64},Float64,Float64,Int64}}) at C:\Users\DELL\.julia\packages\LazyArrays\MAEFm\src\lazybroadcasting.jl:48       
 [10] combine_axes at .\broadcast.jl:473 [inlined]
 [11] instantiate(::Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Nothing,typeof(*),Tuple{BroadcastArray{Float64,1,typeof(+),Tuple{InfiniteArrays.InfStepRange{Int64,Int64},Float64,Float64,Int64}},BroadcastArray{Float64,1,typeof(+),Tuple{InfiniteArrays.InfStepRange{Int64,Int64},Float64,Float64,Int64}}}}) at .\broadcast.jl:256
 [12] Base.Broadcast.Broadcasted(::BroadcastArray{Float64,1,typeof(*),Tuple{BroadcastArray{Float64,1,typeof(+),Tuple{InfiniteArrays.InfStepRange{Int64,Int64},Float64,Float64,Int64}},BroadcastArray{Float64,1,typeof(+),Tuple{InfiniteArrays.InfStepRange{Int64,Int64},Float64,Float64,Int64}}}}) at C:\Users\DELL\.julia\packages\LazyArrays\MAEFm\src\lazybroadcasting.jl:43
 [13] axes(::BroadcastArray{Float64,1,typeof(*),Tuple{BroadcastArray{Float64,1,typeof(+),Tuple{InfiniteArrays.InfStepRange{Int64,Int64},Float64,Float64,Int64}},BroadcastArray{Float64,1,typeof(+),Tuple{InfiniteArrays.InfStepRange{Int64,Int64},Float64,Float64,Int64}}}}) at C:\Users\DELL\.julia\packages\LazyArrays\MAEFm\src\lazybroadcasting.jl:48
 [14] newindex at .\broadcast.jl:537 [inlined]
 [15] _broadcast_getindex at .\broadcast.jl:584 [inlined]
 [16] _getindex at .\broadcast.jl:628 [inlined]
 [17] _getindex at .\broadcast.jl:627 [inlined]
 [18] _broadcast_getindex at .\broadcast.jl:603 [inlined]
 [19] getindex(::Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(/),Tuple{BroadcastArray{Float64,1,typeof(*),Tuple{Int64,InfiniteArrays.InfUnitRange{Int64},BroadcastArray{Float64,1,typeof(+),Tuple{InfiniteArrays.InfUnitRange{Int64},Float64,Float64,Int64}}}},BroadcastArray{Float64,1,typeof(*),Tuple{BroadcastArray{Float64,1,typeof(+),Tuple{InfiniteArrays.InfStepRange{Int64,Int64},Float64,Float64,Int64}},BroadcastArray{Float64,1,typeof(+),Tuple{InfiniteArrays.InfStepRange{Int64,Int64},Float64,Float64,Int64}}}}}}, ::Int64) at .\broadcast.jl:564
 [20] getindex(::BroadcastArray{Float64,1,typeof(/),Tuple{BroadcastArray{Float64,1,typeof(*),Tuple{Int64,InfiniteArrays.InfUnitRange{Int64},BroadcastArray{Float64,1,typeof(+),Tuple{InfiniteArrays.InfUnitRange{Int64},Float64,Float64,Int64}}}},BroadcastArray{Float64,1,typeof(*),Tuple{BroadcastArray{Float64,1,typeof(+),Tuple{InfiniteArrays.InfStepRange{Int64,Int64},Float64,Float64,Int64}},BroadcastArray{Float64,1,typeof(+),Tuple{InfiniteArrays.InfStepRange{Int64,Int64},Float64,Float64,Int64}}}}}}, ::Int64) at C:\Users\DELL\.julia\packages\LazyArrays\MAEFm\src\lazybroadcasting.jl:52
 [21] vcat_getindex at C:\Users\DELL\.julia\packages\LazyArrays\MAEFm\src\lazyconcat.jl:51 [inlined]
 [22] getindex at C:\Users\DELL\.julia\packages\LazyArrays\MAEFm\src\lazyconcat.jl:68 [inlined]
 [23] forwardrecurrence!(::Array{Float64,1}, ::BroadcastArray{Float64,1,typeof(/),Tuple{BroadcastArray{Float64,1,typeof(*),Tuple{Int64,InfiniteArrays.InfStepRange{Float64,Float64},InfiniteArrays.InfStepRange{Float64,Float64}}},BroadcastArray{Float64,1,typeof(*),Tuple{BroadcastArray{Float64,1,typeof(+),Tuple{InfiniteArrays.InfStepRange{Int64,Int64},Float64,Float64}},BroadcastArray{Float64,1,typeof(+),Tuple{InfiniteArrays.InfStepRange{Int64,Int64},Float64,Float64,Int64}}}}}}, ::LazyArrays.ApplyArray{Float64,1,typeof(vcat),Tuple{Float64,BroadcastArray{Float64,1,typeof(/),Tuple{Float64,BroadcastArray{Float64,1,typeof(*),Tuple{InfiniteArrays.InfStepRange{Float64,Float64},InfiniteArrays.InfStepRange{Float64,Float64}}}}}}}, ::LazyArrays.ApplyArray{Float64,1,typeof(vcat),Tuple{Float64,BroadcastArray{Float64,1,typeof(/),Tuple{BroadcastArray{Float64,1,typeof(*),Tuple{Int64,InfiniteArrays.InfUnitRange{Int64},BroadcastArray{Float64,1,typeof(+),Tuple{InfiniteArrays.InfUnitRange{Int64},Float64,Float64,Int64}}}},BroadcastArray{Float64,1,typeof(*),Tuple{BroadcastArray{Float64,1,typeof(+),Tuple{InfiniteArrays.InfStepRange{Int64,Int64},Float64,Float64,Int64}},BroadcastArray{Float64,1,typeof(+),Tuple{InfiniteArrays.InfStepRange{Int64,Int64},Float64,Float64,Int64}}}}}}}}, ::Float64, ::Int64) at C:\Users\DELL\.julia\packages\OrthogonalPolynomialsQuasi\llrOw\src\OrthogonalPolynomialsQuasi.jl:146
 [24] forwardrecurrence!(::Array{Float64,1}, ::BroadcastArray{Float64,1,typeof(/),Tuple{BroadcastArray{Float64,1,typeof(*),Tuple{Int64,InfiniteArrays.InfStepRange{Float64,Float64},InfiniteArrays.InfStepRange{Float64,Float64}}},BroadcastArray{Float64,1,typeof(*),Tuple{BroadcastArray{Float64,1,typeof(+),Tuple{InfiniteArrays.InfStepRange{Int64,Int64},Float64,Float64}},BroadcastArray{Float64,1,typeof(+),Tuple{InfiniteArrays.InfStepRange{Int64,Int64},Float64,Float64,Int64}}}}}}, ::LazyArrays.ApplyArray{Float64,1,typeof(vcat),Tuple{Float64,BroadcastArray{Float64,1,typeof(/),Tuple{Float64,BroadcastArray{Float64,1,typeof(*),Tuple{InfiniteArrays.InfStepRange{Float64,Float64},InfiniteArrays.InfStepRange{Float64,Float64}}}}}}}, ::LazyArrays.ApplyArray{Float64,1,typeof(vcat),Tuple{Float64,BroadcastArray{Float64,1,typeof(/),Tuple{BroadcastArray{Float64,1,typeof(*),Tuple{Int64,InfiniteArrays.InfUnitRange{Int64},BroadcastArray{Float64,1,typeof(+),Tuple{InfiniteArrays.InfUnitRange{Int64},Float64,Float64,Int64}}}},BroadcastArray{Float64,1,typeof(*),Tuple{BroadcastArray{Float64,1,typeof(+),Tuple{InfiniteArrays.InfStepRange{Int64,Int64},Float64,Float64,Int64}},BroadcastArray{Float64,1,typeof(+),Tuple{InfiniteArrays.InfStepRange{Int64,Int64},Float64,Float64,Int64}}}}}}}}, ::Float64) at C:\Users\DELL\.julia\packages\OrthogonalPolynomialsQuasi\llrOw\src\OrthogonalPolynomialsQuasi.jl:136
 [25] getindex(::Jacobi{Float64}, ::Float64, ::Base.OneTo{Int64}) at C:\Users\DELL\.julia\packages\OrthogonalPolynomialsQuasi\llrOw\src\OrthogonalPolynomialsQuasi.jl:200
 [26] getindex at C:\Users\DELL\.julia\packages\OrthogonalPolynomialsQuasi\llrOw\src\OrthogonalPolynomialsQuasi.jl:240 [inlined]
 [27] getindex at C:\Users\DELL\.julia\packages\QuasiArrays\4gS8F\src\subquasiarray.jl:215 [inlined]
 [28] macro expansion at .\multidimensional.jl:756 [inlined]
 [29] macro expansion at .\cartesian.jl:64 [inlined]
 [30] macro expansion at .\multidimensional.jl:751 [inlined]
 [31] _unsafe_getindex! at .\multidimensional.jl:747 [inlined]
 [32] _unsafe_getindex(::IndexCartesian, ::SubArray{Float64,1,Jacobi{Float64},Tuple{Float64,Base.Slice{InfiniteArrays.OneToInf{Int64}}},false}, ::InfiniteArrays.InfUnitRange{Int64}) at .\multidimensional.jl:741
 [33] _getindex at .\multidimensional.jl:727 [inlined]
 [34] getindex(::SubArray{Float64,1,Jacobi{Float64},Tuple{Float64,Base.Slice{InfiniteArrays.OneToInf{Int64}}},false}, ::InfiniteArrays.InfUnitRange{Int64}) at .\abstractarray.jl:980
 [35] top-level scope at none:0

Chebyshev weight times polynomial

w₁ = ChebyshevTWeight()
x = axes(w1,1)
w₂ = OrthogonalPolynomialsQuasi.JacobiWeight(-0.5,-0.5)
ϕ = x.^3 - 2*x.^2 - x .+ 2
LanczosPolynomial.*w₂)

works but

LanczosPolynomial.*w₁)
MethodError: no method matching singularitiesbroadcast(::typeof(*), ::LegendreWeight{Float64}, ::OrthogonalPolynomialsQuasi.ChebyshevWeight{1,Float64})
Closest candidates are:
  singularitiesbroadcast(::typeof(*), ::LegendreWeight, !Matched::OrthogonalPolynomialsQuasi.JacobiWeight) at C:\Users\mfaso\.julia\packages\OrthogonalPolynomialsQuasi\bIORM\src\jacobi.jl:79
  singularitiesbroadcast(::typeof(*), ::LegendreWeight, !Matched::OrthogonalPolynomialsQuasi.NoSingularities) at In[3]:9
  singularitiesbroadcast(::Any, ::LegendreWeight) at C:\Users\mfaso\.julia\packages\OrthogonalPolynomialsQuasi\bIORM\src\jacobi.jl:59
  ...

Stacktrace:
 [1] singularities(::LazyArrays.BroadcastLayout{typeof(*)}, ::BroadcastQuasiArray{Float64,1,typeof(*),Tuple{BroadcastQuasiArray{Float64,1,typeof(+),Tuple{BroadcastQuasiArray{Float64,1,typeof(-),Tuple{BroadcastQuasiArray{Float64,1,typeof(-),Tuple{BroadcastQuasiArray{Float64,1,typeof(Base.literal_pow),Tuple{Base.RefValue{typeof(^)},Inclusion{Float64,ChebyshevInterval{Float64}},Base.RefValue{Val{3}}}},BroadcastQuasiArray{Float64,1,typeof(*),Tuple{Int64,BroadcastQuasiArray{Float64,1,typeof(Base.literal_pow),Tuple{Base.RefValue{typeof(^)},Inclusion{Float64,ChebyshevInterval{Float64}},Base.RefValue{Val{2}}}}}}}},Inclusion{Float64,ChebyshevInterval{Float64}}}},Int64}},OrthogonalPolynomialsQuasi.ChebyshevWeight{1,Float64}}}) at C:\Users\mfaso\.julia\packages\OrthogonalPolynomialsQuasi\bIORM\src\OrthogonalPolynomialsQuasi.jl:165
 [2] singularities(::BroadcastQuasiArray{Float64,1,typeof(*),Tuple{BroadcastQuasiArray{Float64,1,typeof(+),Tuple{BroadcastQuasiArray{Float64,1,typeof(-),Tuple{BroadcastQuasiArray{Float64,1,typeof(-),Tuple{BroadcastQuasiArray{Float64,1,typeof(Base.literal_pow),Tuple{Base.RefValue{typeof(^)},Inclusion{Float64,ChebyshevInterval{Float64}},Base.RefValue{Val{3}}}},BroadcastQuasiArray{Float64,1,typeof(*),Tuple{Int64,BroadcastQuasiArray{Float64,1,typeof(Base.literal_pow),Tuple{Base.RefValue{typeof(^)},Inclusion{Float64,ChebyshevInterval{Float64}},Base.RefValue{Val{2}}}}}}}},Inclusion{Float64,ChebyshevInterval{Float64}}}},Int64}},OrthogonalPolynomialsQuasi.ChebyshevWeight{1,Float64}}}) at C:\Users\mfaso\.julia\packages\OrthogonalPolynomialsQuasi\bIORM\src\OrthogonalPolynomialsQuasi.jl:167
 [3] LanczosPolynomial(::BroadcastQuasiArray{Float64,1,typeof(*),Tuple{BroadcastQuasiArray{Float64,1,typeof(+),Tuple{BroadcastQuasiArray{Float64,1,typeof(-),Tuple{BroadcastQuasiArray{Float64,1,typeof(-),Tuple{BroadcastQuasiArray{Float64,1,typeof(Base.literal_pow),Tuple{Base.RefValue{typeof(^)},Inclusion{Float64,ChebyshevInterval{Float64}},Base.RefValue{Val{3}}}},BroadcastQuasiArray{Float64,1,typeof(*),Tuple{Int64,BroadcastQuasiArray{Float64,1,typeof(Base.literal_pow),Tuple{Base.RefValue{typeof(^)},Inclusion{Float64,ChebyshevInterval{Float64}},Base.RefValue{Val{2}}}}}}}},Inclusion{Float64,ChebyshevInterval{Float64}}}},Int64}},OrthogonalPolynomialsQuasi.ChebyshevWeight{1,Float64}}}) at C:\Users\mfaso\.julia\packages\OrthogonalPolynomialsQuasi\bIORM\src\lanczos.jl:207
 [4] top-level scope at In[29]:1
 [5] include_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1091
 [6] execute_code(::String, ::String) at C:\Users\mfaso\.julia\packages\IJulia\a1SNk\src\execute_request.jl:27
 [7] execute_request(::ZMQ.Socket, ::IJulia.Msg) at C:\Users\mfaso\.julia\packages\IJulia\a1SNk\src\execute_request.jl:86
 [8] #invokelatest#1 at .\essentials.jl:710 [inlined]
 [9] invokelatest at .\essentials.jl:709 [inlined]
 [10] eventloop(::ZMQ.Socket) at C:\Users\mfaso\.julia\packages\IJulia\a1SNk\src\eventloop.jl:8
 [11] (::IJulia.var"#15#18")() at .\task.jl:356

(v0.4.3) Dependency issue

The current release v0.4.3 causes many other packages to downgrade

(@JuliaPro_v1.5.4-1) pkg> status
Status `C:\Users\pty\.julia\environments\JuliaPro_v1.5.4-1\Project.toml`
  [4c555306] ArrayLayouts v0.6.5
  [c52e3926] Atom v0.12.30 ⚲
  [aae01518] BandedMatrices v0.16.8
  [8e7c35d0] BlockArrays v0.15.2
  [ffab5731] BlockBandedMatrices v0.10.4
  [057dd010] FastTransforms v0.12.3
  [1a297f60] FillArrays v0.11.7
  [7073ff75] IJulia v1.23.2
  [cde9dba0] InfiniteLinearAlgebra v0.5.5
  [4138dd39] JLD v0.12.3
  [e5e0dc1b] Juno v0.8.4 ⚲
  [8e2b3108] KahanSummation v0.2.0
  [5078a376] LazyArrays v0.21.3
  [d7e5e226] LazyBandedMatrices v0.5.5
  [23992714] MAT v0.10.1
  [4722fa14] PkgAuthentication v1.1.1
  [91a5bcdd] Plots v1.12.0
  [90137ffa] StaticArrays v1.1.1
  [44d3d7a6] Weave v0.10.7

(@JuliaPro_v1.5.4-1) pkg> add OrthogonalPolynomialsQuasi
  Resolving package versions...
Updating `C:\Users\pty\.julia\environments\JuliaPro_v1.5.4-1\Project.toml`
  [4c555306]  ArrayLayouts v0.6.5  v0.5.4
  [8e7c35d0]  BlockArrays v0.15.2  v0.14.5
  [ffab5731]  BlockBandedMatrices v0.10.4  v0.10.3
  [057dd010]  FastTransforms v0.12.3  v0.11.3
  [cde9dba0]  InfiniteLinearAlgebra v0.5.5  v0.4.8
  [5078a376]  LazyArrays v0.21.3  v0.20.5
  [d7e5e226]  LazyBandedMatrices v0.5.5  v0.4.7
  [aa41a628] + OrthogonalPolynomialsQuasi v0.4.3
Updating `C:\Users\pty\.julia\environments\JuliaPro_v1.5.4-1\Manifest.toml`
  [4fba245c] + ArrayInterface v3.1.7
  [4c555306]  ArrayLayouts v0.6.5  v0.5.4
  [8e7c35d0]  BlockArrays v0.15.2  v0.14.5
  [ffab5731]  BlockBandedMatrices v0.10.4  v0.10.3
  [7ae1f121] + ContinuumArrays v0.4.2
  [5b8099bc] + DomainSets v0.4.2
  [da5c29d0] + EllipsisNotation v1.1.0
  [057dd010]  FastTransforms v0.12.3  v0.11.3
  [34b6f7d7]  FastTransforms_jll v0.5.1+0  v0.4.1+0
  [615f187c] + IfElse v0.1.0
  [4858937d]  InfiniteArrays v0.10.6  v0.9.5
  [cde9dba0]  InfiniteLinearAlgebra v0.5.5  v0.4.8
  [e1ba4f0e] - Infinities v0.1.0
  [8197267c] + IntervalSets v0.5.3
  [5078a376]  LazyArrays v0.21.3  v0.20.5
  [d7e5e226]  LazyBandedMatrices v0.5.5  v0.4.7
  [aa41a628] + OrthogonalPolynomialsQuasi v0.4.3
  [c4ea9172] + QuasiArrays v0.3.9
  [aedffcd0] + Static v0.2.4

Support P^(a,b) magic?

julia> Base.:^(P::Legendre, (a,b)) = Jacobi(a,b)

julia> P = Legendre()
Legendre{Float64}

julia> P^(0.1,0.2)
Jacobi{Float64}

Rename to ClassicalOrthogonalPolynomials.jl?

This seems more in line with MultivariateOrthogonalPolynomials.jl, SemiclassicalOrthogonalPolynomials.jl.

Eventually, this will displace/merge with ApproxFunOrthogonalPolynomials.jl?

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.