juliaapproximation / classicalorthogonalpolynomials.jl Goto Github PK
View Code? Open in Web Editor NEWA Julia package for classical orthogonal polynomials and expansions
License: MIT License
classicalorthogonalpolynomials.jl's People
Contributors
![github-actions[bot] avatar](https://avatars.githubusercontent.com/in/15368?v=4 "github-actions[bot]")
Stargazers
Watchers
classicalorthogonalpolynomials.jl's Issues
What's the structure of "OffHilbert"?
@MikaelSlevinsky Any thoughts on how we o]should treat the structure of "OffHilbert"?
At the moment it's just treated as dense, which is not so nice for adaptive linear algebra. I could introduce a notion of "compact".
Make AbstractOPLayout <: AbstractPolynomialLayout
QuasiArrays.jl now has PolynomialLayout, which currently only supports Inclusion but is designed to represent general polynomials. This package has AbstractOPLayout which should be a special class of polynomial layout.
This helps make sense of how WeightedBasisLayout should work. E.g. equality can be made precise (if w1 and w2 are not polynomial and P and Q are non-zero polynomial quasimatrices then w1 .* P == w2 .* Q iff w1 == w2 and P == Q)
Fourier{BigFloat}
Is there a way to work with Fourier and BigFloat? It goes via FFTW which is causing some issues.
using ClassicalOrthogonalPolynomials
F = Fourier{BigFloat}()
F[:, 1:10] \ sin.(axes(F,1))
ERROR: MethodError: no method matching mul!(::Vector{BigFloat}, ::FFTW.r2rFFTWPlan{Float64, Vector{Int32}, false, 1, Int64}, ::Vector{BigFloat}, ::Bool, ::Bool)
Closest candidates are:
mul!(::StridedVecOrMat, ::SparseArrays.AbstractSparseMatrixCSC, ::Union{LinearAlgebra.Adjoint{<:Any, <:Union{LinearAlgebra.LowerTriangular, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedMatrix, BitMatrix}}, LinearAlgebra.LowerTriangular, LinearAlgebra.Transpose{<:Any, <:Union{LinearAlgebra.LowerTriangular, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedMatrix, BitMatrix}}, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedVector, StridedMatrix, BitMatrix, BitVector}, ::Number, ::Number)
@ SparseArrays C:\Users\john.papad\AppData\Local\Programs\Julia-1.9.2\share\julia\stdlib\v1.9\SparseArrays\src\linalg.jl:30
mul!(::StridedVecOrMat, ::LinearAlgebra.Adjoint{<:Any, <:SparseArrays.AbstractSparseMatrixCSC}, ::Union{LinearAlgebra.Adjoint{<:Any, <:Union{LinearAlgebra.LowerTriangular, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedMatrix, BitMatrix}}, LinearAlgebra.LowerTriangular, LinearAlgebra.Transpose{<:Any, <:Union{LinearAlgebra.LowerTriangular, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedMatrix, BitMatrix}}, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedVector, StridedMatrix, BitMatrix, BitVector}, ::Number, ::Number)
@ SparseArrays C:\Users\john.papad\AppData\Local\Programs\Julia-1.9.2\share\julia\stdlib\v1.9\SparseArrays\src\linalg.jl:56
mul!(::StridedVecOrMat, ::LinearAlgebra.Transpose{<:Any, <:SparseArrays.AbstractSparseMatrixCSC}, ::Union{LinearAlgebra.Adjoint{<:Any, <:Union{LinearAlgebra.LowerTriangular, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedMatrix, BitMatrix}}, LinearAlgebra.LowerTriangular, LinearAlgebra.Transpose{<:Any, <:Union{LinearAlgebra.LowerTriangular, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedMatrix, BitMatrix}}, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedVector, StridedMatrix, BitMatrix, BitVector}, ::Number, ::Number)
@ SparseArrays C:\Users\john.papad\AppData\Local\Programs\Julia-1.9.2\share\julia\stdlib\v1.9\SparseArrays\src\linalg.jl:56
...
Stacktrace:
[1] mul!
@ C:\Users\john.papad\AppData\Local\Programs\Julia-1.9.2\share\julia\stdlib\v1.9\LinearAlgebra\src\matmul.jl:276 [inlined]
[2] mul!(ret::Vector{BigFloat}, F::ClassicalOrthogonalPolynomials.ShuffledR2HC{BigFloat, FFTW.r2rFFTWPlan{Float64, Vector{Int32}, false, 1, Int64}}, bin::Vector{Float64})
@ ClassicalOrthogonalPolynomials C:\Users\john.papad\.julia\packages\ClassicalOrthogonalPolynomials\KZ9ds\src\classical\fourier.jl:128
[3] *(F::ClassicalOrthogonalPolynomials.ShuffledR2HC{BigFloat, FFTW.r2rFFTWPlan{Float64, Vector{Int32}, false, 1, Int64}}, b::Vector{Float64})
@ ClassicalOrthogonalPolynomials C:\Users\john.papad\.julia\packages\ClassicalOrthogonalPolynomials\KZ9ds\src\classical\fourier.jl:133
[4] \(a::ContinuumArrays.TransformFactorization{BigFloat, Vector{Float64}, ClassicalOrthogonalPolynomials.ShuffledR2HC{BigFloat, FFTW.r2rFFTWPlan{Float64, Vector{Int32}, false, 1, Int64}}}, b::QuasiArrays.BroadcastQuasiVector{Float64, typeof(f), Tuple{Inclusion{Float64, DomainSets.RealNumbers}}})
@ ContinuumArrays C:\Users\john.papad\.julia\packages\ContinuumArrays\6Nbuj\src\plans.jl:26
[5] \(a::ContinuumArrays.ProjectionFactorization{BigFloat, ContinuumArrays.TransformFactorization{BigFloat, Vector{Float64}, ClassicalOrthogonalPolynomials.ShuffledR2HC{BigFloat, FFTW.r2rFFTWPlan{Float64, Vector{Int32}, false, 1, Int64}}}, UnitRange{Int64}}, b::QuasiArrays.BroadcastQuasiVector{Float64, typeof(f), Tuple{Inclusion{Float64, DomainSets.RealNumbers}}})
@ ContinuumArrays C:\Users\john.papad\.julia\packages\ContinuumArrays\6Nbuj\src\bases\bases.jl:213
[6] transform_ldiv_size(#unused#::Tuple{Infinities.InfiniteCardinal{1}, Int64}, A::QuasiArrays.SubQuasiArray{BigFloat, 2, Fourier{BigFloat}, Tuple{Inclusion{Float64, DomainSets.RealNumbers}, UnitRange{Int64}}, false}, B::QuasiArrays.BroadcastQuasiVector{Float64, typeof(f), Tuple{Inclusion{Float64, DomainSets.RealNumbers}}})
@ ContinuumArrays C:\Users\john.papad\.julia\packages\ContinuumArrays\6Nbuj\src\bases\bases.jl:260
[7] transform_ldiv(A::QuasiArrays.SubQuasiArray{BigFloat, 2, Fourier{BigFloat}, Tuple{Inclusion{Float64, DomainSets.RealNumbers}, UnitRange{Int64}}, false}, B::QuasiArrays.BroadcastQuasiVector{Float64, typeof(f), Tuple{Inclusion{Float64, DomainSets.RealNumbers}}})
@ ContinuumArrays C:\Users\john.papad\.julia\packages\ContinuumArrays\6Nbuj\src\bases\bases.jl:261
[8] basis_ldiv_size
@ C:\Users\john.papad\.julia\packages\ContinuumArrays\6Nbuj\src\bases\bases.jl:310 [inlined]
[9] copy
@ C:\Users\john.papad\.julia\packages\ContinuumArrays\6Nbuj\src\bases\bases.jl:302 [inlined]
[10] #materialize#13
@ C:\Users\john.papad\.julia\packages\ArrayLayouts\tD7jV\src\ldiv.jl:22 [inlined]
[11] materialize
@ C:\Users\john.papad\.julia\packages\ArrayLayouts\tD7jV\src\ldiv.jl:22 [inlined]
[12] #ldiv#20
@ C:\Users\john.papad\.julia\packages\ArrayLayouts\tD7jV\src\ldiv.jl:98 [inlined]
[13] ldiv
@ C:\Users\john.papad\.julia\packages\ArrayLayouts\tD7jV\src\ldiv.jl:98 [inlined]
[14] \(A::QuasiArrays.SubQuasiArray{BigFloat, 2, Fourier{BigFloat}, Tuple{Inclusion{Float64, DomainSets.RealNumbers}, UnitRange{Int64}}, false}, B::QuasiArrays.BroadcastQuasiVector{Float64, typeof(f), Tuple{Inclusion{Float64, DomainSets.RealNumbers}}})
@ QuasiArrays C:\Users\john.papad\.julia\packages\QuasiArrays\GYEdf\src\matmul.jl:34
[15] top-level scope
@ Untitled-1:6
Make adaptive transform a Vcat instead of a cached vector?
Sometimes we use the transform to build operators, and the bandwidth is determined by the number of non-zeros. In this case, cached vectors are not so useful as anytime we access them they grow, and the bandwidth of the operators grows.
But is there any reason we would need coefficient vectors to be mutable (in their zeros) by default? A user could call Base.copymutable(P \ f) if they really wanted to mutate.
Add Associated(::OrthogonalPolynomial)
I think we want to support e.g. Associated(Legendre()) for the associated OPs. It should be straight forward;
struct Associated{T, PP} <: OrthogonalPolynomial{T}
P::PP
end
function jacobimatrix(Q::Associated)
X = jacobimatrix(Q.P)
c,a,b = subdiagonaldata(X), diagonaldata(X), supdiagonaldata(X)
Tridiagonal(c[2:end], a[2:end], b[2:end])
end@MikaelSlevinsky Are there any other features I should add? Do you have that the transforms implemented yet?
We could also support Associated(Legendre(), m) for more general case.
Make P[0.1,10] allocation free
At the moment we unnecessarily build P[0.1,1:10] but the use of forwardrecurrence would avoid this
`points(Ultraspherical(1, -1..1), n)` returns points on Chebyshev grid of the first kind and not the second kind
using ApproxFunOrthogonalPolynomials
# For reference
chebyshevpoints(Float64, 3, Val(1)) # [0.866, 0.0, -0.866]
chebyshevpoints(Float64, 3, Val(2)) # [1.0, 0.0, -1.0]
# These should be points corresponding to the Chebyshev polynomials of the second kind, I thought?
points(Ultraspherical(1, -1..1), 3) # [0.866, 0.0, -0.866]Add roots
We can port the code
https://github.com/JuliaApproximation/ApproxFunOrthogonalPolynomials.jl/blob/master/src/roots.jl
I think call it findall(iszero,f).
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!
How to get BigFloat accurate conversion operators?
For OPs on y⁴ = 1 - x⁴ I need conversion operators between ChebyshevT() and certain semiclassical OPs, e.g.,
T = ChebyshevT{BigFloat}()
x = axes(T,1)
P₁ = LanczosPolynomial( (1 .+ x.^2).^(0.5) )
R₁ = P₁\TThis is accurate to only Float64 precision:
xv = big"0.1"
T[xv,3] - ( P₁[xv,1]*R₁[1,3] + P₁[xv,3]*R₁[3,3] )
-6.05859383536196727705447569248701020544252033186250536144991618812281214296214e-17Another example:
P = Legendre{BigFloat}()
R = P\TThe entry R[3,3] is 4/3 but it's accurate to only Float64 precision:
R[3,3]
1.333333333333333250903485233170513237354072211061085640237818614115573313606597does not work with `Symbols.jl` due to type casting (it seems)
using Symbolics
using ClassicalOrthogonalPolynomials
P = Legendre()
@variables z
P[z,1:10]breaks with
MethodError: no method matching Float64(::Symbolics.Num)
Closest candidates are:
(::Type{T})(::Real, !Matched::RoundingMode) where T<:AbstractFloat at rounding.jl:200
(::Type{T})(::T) where T<:Number at boot.jl:760
(::Type{T})(!Matched::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50
...
convert(::Type{Float64}, ::Symbolics.Num)@number.jl:7
to_quasi_index(::Type{Float64}, ::Symbolics.Num)@indices.jl:103
to_quasi_index(::ClassicalOrthogonalPolynomials.Legendre{Float64}, ::Type, ::Symbolics.Num)@indices.jl:111
to_indices@indices.jl:118[inlined]
to_indices@multidimensional.jl:413[inlined]
view@subquasiarray.jl:67[inlined]
layout_getindex@ArrayLayouts.jl:108[inlined]
getindex(::ClassicalOrthogonalPolynomials.Legendre{Float64}, ::Symbolics.Num, ::UnitRange{Int64})@clenshaw.jl:67
top-level scope@Local: 1[inlined]long compilation time for Jacobi(m,n) \ Jacobi(0,0)
julia> using ClassicalOrthogonalPolynomials
julia> @time Jacobi(2,2) \ Jacobi(0,0);
0.206673 seconds (326.39 k allocations: 19.069 MiB, 98.56% compilation time)
julia> @time Jacobi(5,5) \ Jacobi(0,0);
0.435400 seconds (788.60 k allocations: 43.161 MiB, 98.64% compilation time)
julia> @time Jacobi(10,10) \ Jacobi(0,0);
16.587771 seconds (26.72 M allocations: 1.142 GiB, 1.86% gc time, 99.94% compilation time)
julia> @time Jacobi(11,11) \ Jacobi(0,0);
46.804348 seconds (75.86 M allocations: 3.206 GiB, 1.57% gc time, 100.00% compilation time)
julia> @time Jacobi(11,11) \ Jacobi(0,0);
0.000060 seconds (43 allocations: 20.266 KiB)This could possibly be caused by LazyArrays.ApplyArray using NTuple.
NaN results from Jacobi expansion
julia> @time Chebyshev()\JacobiWeight(0,0.5)
0.027141 seconds (3.08 k allocations: 3.701 MiB)
ℵ₀-element LazyArrays.CachedArray{Float64, 1, Vector{Float64}, FillArrays.Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}} with indices OneToInf():
0.9003163161786569
0.600210877394969
-0.12004217544451247
0.05144664659444724
-0.028581470311092174
⋮
julia> @time Jacobi(-0.5,-0.5)\JacobiWeight(0,0.5)
14.685838 seconds (3.11 k allocations: 4.702 MiB)
ℵ₀-element LazyArrays.CachedArray{Float64, 1, Vector{Float64}, FillArrays.Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}} with indices OneToInf():
NaN
NaN
NaN
NaN
NaN
⋮
julia> @time Legendre()\JacobiWeight(0,0.5)
15.032195 seconds (3.09 k allocations: 4.702 MiB)
ℵ₀-element LazyArrays.CachedArray{Float64, 1, Vector{Float64}, FillArrays.Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}} with indices OneToInf():
0.9428090415820628
0.5656854249492391
-0.13468700594029615
0.0628539361054728
-0.036732819801901045
⋮
julia> @time Jacobi(0,0)\JacobiWeight(0,0.5)
14.450656 seconds (3.11 k allocations: 4.702 MiB)
ℵ₀-element LazyArrays.CachedArray{Float64, 1, Vector{Float64}, FillArrays.Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}} with indices OneToInf():
NaN
NaN
NaN
NaN
NaN
⋮Shifted Legendre{BigFloat}() error
I'm seeing some strange behavior for Legendre polynomials with BigFloat precision shifted to 0..1. Quite possible that something about my syntax isn't how it is intended to be used but it might also be bug.
Here's a simple example that reproduces the error:
julia> # float64
julia> P = Legendre()[affine(Inclusion(0..1), axes(Legendre(),1)), :]
Legendre{Float64} affine mapped to 0..1
julia> x = axes(P,1)
Inclusion(0..1)
julia> P \ sin.(x)
ℵ₀-element LazyArrays.CachedArray{Float64,1,Array{Float64,1},Zeros{Float64,1,Tuple{InfiniteArrays.OneToInf{Int64}}}} with indices OneToInf():
0.4596976941318605
0.42791899124296007
-0.03924363301542014
-0.007212191228055754
0.0002821450243386184
2.8742703093836133e-5
-7.146539529763204e-7
-5.0363463994812964e-8
9.178337836736593e-10
4.944523575181289e-11
-7.110265998998354e-13
-3.106481507574827e-14
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
⋮
julia> # bigfloat
julia> Pbig = Legendre{BigFloat}()[affine(Inclusion(BigInt(0)..BigInt(1)), axes(Legendre{BigFloat}(),1)), :]
Legendre{BigFloat} affine mapped to 0..1
julia> xbig = axes(Pbig,1)
Inclusion(0..1)
julia> Pbig \ sin.(xbig)
ERROR: cannot resize array with shared data
Stacktrace:
[1] _deleteend! at ./array.jl:901 [inlined]
[2] resize! at ./array.jl:1090 [inlined]
[3] chop!(::Array{BigFloat,1}, ::BigFloat) at /home/timon/.julia/packages/ClassicalOrthogonalPolynomials/YpJSc/src/ClassicalOrthogonalPolynomials.jl:71
[4] adaptivetransform_ldiv(::SubQuasiArray{BigFloat,2,Legendre{BigFloat},Tuple{ContinuumArrays.AffineMap{BigFloat,Inclusion{BigFloat,IntervalSets.Interval{:closed,:closed,BigInt}},Inclusion{BigFloat,ChebyshevInterval{BigFloat}}},Base.Slice{InfiniteArrays.OneToInf{Int64}}},false}, ::BroadcastQuasiArray{BigFloat,1,typeof(sin),Tuple{Inclusion{BigFloat,IntervalSets.Interval{:closed,:closed,BigInt}}}}) at /home/timon/.julia/packages/ClassicalOrthogonalPolynomials/YpJSc/src/ClassicalOrthogonalPolynomials.jl:104
[5] transform_ldiv at /home/timon/.julia/packages/ClassicalOrthogonalPolynomials/YpJSc/src/ClassicalOrthogonalPolynomials.jl:64 [inlined]
[6] transform_ldiv at /home/timon/.julia/packages/ContinuumArrays/D2vwQ/src/bases/bases.jl:178 [inlined]
[7] copy at /home/timon/.julia/packages/ContinuumArrays/D2vwQ/src/bases/bases.jl:180 [inlined]
[8] materialize(::Ldiv{ContinuumArrays.MappedBasisLayout,LazyArrays.BroadcastLayout{typeof(sin)},SubQuasiArray{BigFloat,2,Legendre{BigFloat},Tuple{ContinuumArrays.AffineMap{BigFloat,Inclusion{BigFloat,IntervalSets.Interval{:closed,:closed,BigInt}},Inclusion{BigFloat,ChebyshevInterval{BigFloat}}},Base.Slice{InfiniteArrays.OneToInf{Int64}}},false},BroadcastQuasiArray{BigFloat,1,typeof(sin),Tuple{Inclusion{BigFloat,IntervalSets.Interval{:closed,:closed,BigInt}}}}}) at /home/timon/.julia/packages/ArrayLayouts/7LZ91/src/ldiv.jl:22
[9] ldiv at /home/timon/.julia/packages/ArrayLayouts/7LZ91/src/ldiv.jl:86 [inlined]
[10] \(::SubQuasiArray{BigFloat,2,Legendre{BigFloat},Tuple{ContinuumArrays.AffineMap{BigFloat,Inclusion{BigFloat,IntervalSets.Interval{:closed,:closed,BigInt}},Inclusion{BigFloat,ChebyshevInterval{BigFloat}}},Base.Slice{InfiniteArrays.OneToInf{Int64}}},false}, ::BroadcastQuasiArray{BigFloat,1,typeof(sin),Tuple{Inclusion{BigFloat,IntervalSets.Interval{:closed,:closed,BigInt}}}}) at /home/timon/.julia/packages/QuasiArrays/Kytgm/src/matmul.jl:34
[11] top-level scope at REPL[77]:1
julia> # but this appears to work
julia> Pbig[:,1:20] \ sin.(xbig)
20-element Array{BigFloat,1}:
0.459697694131860282599063392557024291715022853422108373904832662310088291651591
0.4279189912429598877122041074528643054476364481970978558557517258602675648095019
-0.03924363301542034337341694172731222575488871592208301467640058895035480424366479
-0.007212191228055742704936278556706084716636823072871192246094339888763058783659867
0.000282145024338535280236421222451671910237190823901368171354926761601708489205559
2.874270309358435249957361731260692515947881520105291315102420408358148420958176e-05
-7.146539530749001084413538229440980152040709571024099798175596892344455384798109e-07
-5.036346369111653274248451118754811843127054009898515258792902259125311582439218e-08
9.178336969399690650953052093298209802341083885853023198941221702769575165208465e-10
4.944521184304420861236925934788833915496435733298095933603120844076503786764353e-11
-7.112115015695451912101596871391910849417223830934967825666746546919715622639409e-13
-3.101024774116237937183881067754829492183211029272491093718995995982405455400059e-14
3.684185721261002718333096225757674226208193004456196754596387429774500100182845e-16
1.349202245083837855737893815880041795386546710643032360063629130209568728057782e-17
-1.365328931605316752701716780091452534277698239440733900315583084467897549345079e-19
-4.310049800103575884451734511728530188560316941533188009045818391107703584309708e-21
3.798549690418122116191713913664570774313410824956465599814127761088843678732203e-23
1.053718404267767678750397050992597701028600476518485629180358929672362729284939e-24
-8.224287969317306291875699136774388048268971744012539584301297942139600288084304e-27
-2.035023861243751730600511331372270370404905435888710306620612749522021466923814e-28Ambiguities in ClassicalOrthogonalPolynomials and its dependencies
I've submitted a PR to QuantumOptics.jl which replaces it's own implementation of Hermite polynomials with this package (qojulia/QuantumOpticsBase.jl#137). However adding ClassicalOrthogonalPolynomials as a dependency caused the Aqua tests for ambiguities to break in QuantumOpticsBase due 150 ambiguities total brought in by ClassicalOrthogonalPolynomials and it's dependencies. This can be checked by running
julia> using Test
julia> using ClassicalOrthogonalPolynomials
julia> Test.detect_ambiguities(Base, Core)
The set of packages which produce ambiguities are:
- ArrayLayouts
- BandedMatrices
- BlockArrays
- BlockBandedMatrices
- ClassicalOrthogonalPolynomials
- ContinuumArrays
- DomainSets
- InfiniteArrays
- Infinities
- LazyArrays
- QuasiArrays
- SemiseparableMatrices
Unfortunately there does not appear to be a way to disable ambiguities on a per package basis. It appears all these packages are all under the JuliaApproximation, JuliaArrays, and JuliaMath organizations with overlapping developers so hopefully there's already some awareness of these ambiguities and maybe some plan to eventually fix them?
Quadratically shifted bases
I want to use
It seems a Chebyhev quadratic shift is implemented in the tests, see
I think this should be low effort to get working for at least general Jacobi polynomials. Any objections to me adding this functionality to the package itself rather than just have it defined in tests?
Once we do this I would also be happy to type up some documentation on shifted bases. 😁
Identity mapping wT -> wU
More of a question than an issue. The following works just fine:
Weighted(ChebyshevT())\Weighted(ChebyshevU())
but
Weighted(ChebyshevU())\Weighted(ChebyshevT())
throws an ERROR: OutOfMemoryError() message. Is this expected behaviour? Or do we want something similar to T\U?
extrapolation syntax?
It's sometimes helpful to use an OP on a domain bigger than the support of its weight. What syntax to use?
We are already using Base.unsafe_getindex for extrapolation at a point, so we could support, e.g.,
T_ext = Base.unsafe_getindex(Chebyshev(), Inclusion(-2..2), :)though perhaps that's not explicit enough.
There's also the question of what type T_ext should be. SuperQuasiArray (as opposed to SubQuasiArray)?
Stackoverflow for P'P
Seems to happen specifically for normalized polynomials.
This works:
julia> P = Jacobi(0,1)
Jacobi(0.0, 1.0)
julia> P'P
((inv(ℵ₀×ℵ₀ LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}} with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf())' with indices OneToInf()×OneToInf()) * (ℵ₀×ℵ₀ Diagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, InfiniteArrays.InfStepRange{Int64, Int64}}}} with indices OneToInf()×OneToInf()) * (inv(ℵ₀×ℵ₀ LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}} with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf():
2.0 -1.0 0.666667 -0.5 0.4 -0.333333 0.285714 -0.25 0.222222 …
-1.0 2.0 -1.33333 1.0 -0.8 0.666667 -0.571429 0.5 -0.444444
0.666667 -1.33333 2.0 -1.5 1.2 -1.0 0.857143 -0.75 0.666667
-0.5 1.0 -1.5 2.0 -1.6 1.33333 -1.14286 1.0 -0.888889
0.4 -0.8 1.2 -1.6 2.0 -1.66667 1.42857 -1.25 1.11111
-0.333333 0.666667 -1.0 1.33333 -1.66667 2.0 -1.71429 1.5 -1.33333 …
0.285714 -0.571429 0.857143 -1.14286 1.42857 -1.71429 2.0 -1.75 1.55556
-0.25 0.5 -0.75 1.0 -1.25 1.5 -1.75 2.0 -1.77778
⋮ ⋮ ⋱ This fails:
julia> P = Normalized(Jacobi(0,1))
Normalized(Jacobi(0.0, 1.0))
julia> P'P
ERROR: StackOverflowError:
Stacktrace:
[1] axes
@ ./abstractarray.jl:74 [inlined]
[2] axes(D::Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}})
@ ArrayLayouts ~/.julia/packages/ArrayLayouts/CNneg/src/memorylayout.jl:716
[3] axes
@ ./abstractarray.jl:74 [inlined]
[4] _check_mul_axes(A::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, B::Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}})
@ ArrayLayouts ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:91
[5] check_mul_axes
@ ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:93 [inlined]
[6] instantiate
@ ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:105 [inlined]
[7] materialize(M::ArrayLayouts.Mul{LazyArrays.InvLayout{ArrayLayouts.BidiagonalLayout{LazyArrays.LazyLayout, LazyArrays.LazyLayout}}, ArrayLayouts.DiagonalLayout{LazyArrays.LazyLayout}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}})
@ ArrayLayouts ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:109
[8] mul(A::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, B::Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}})
@ ArrayLayouts ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:110
[9] _twoarg_simplify(#unused#::typeof(*), #unused#::Val{true}, a::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, b::Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}})
@ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:295
[10] _simplify(#unused#::typeof(*), a::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, b::Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}})
@ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:293
[11] simplify(::typeof(*), ::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, ::Vararg{Any})
@ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:292
[12] simplify(M::ArrayLayouts.Mul{LazyArrays.InvLayout{ArrayLayouts.BidiagonalLayout{LazyArrays.LazyLayout, LazyArrays.LazyLayout}}, ArrayLayouts.DiagonalLayout{LazyArrays.LazyLayout}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}})
@ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:302
[13] copy(M::ArrayLayouts.Mul{LazyArrays.InvLayout{ArrayLayouts.BidiagonalLayout{LazyArrays.LazyLayout, LazyArrays.LazyLayout}}, ArrayLayouts.DiagonalLayout{LazyArrays.LazyLayout}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}})
@ LazyBandedMatrices ~/.julia/packages/LazyBandedMatrices/ZaJ9N/src/LazyBandedMatrices.jl:903
--- the last 7 lines are repeated 4041 more times ---
[28301] materialize(M::ArrayLayouts.Mul{LazyArrays.InvLayout{ArrayLayouts.BidiagonalLayout{LazyArrays.LazyLayout, LazyArrays.LazyLayout}}, ArrayLayouts.DiagonalLayout{LazyArrays.LazyLayout}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}})
@ ArrayLayouts ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:109
[28302] mul(A::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, B::Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}})
@ ArrayLayouts ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:110
[28303] _twoarg_simplify(#unused#::typeof(*), #unused#::Val{true}, a::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, b::Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}})
@ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:295
[28304] _simplify
@ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:293 [inlined]
[28305] _tail_simplify
@ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:299 [inlined]
[28306] _most_simplify
@ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:298 [inlined]
[28307] _simplify
@ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:296 [inlined]
[28308] _tail_simplify
@ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:299 [inlined]
[28309] _most_simplify(#unused#::Val{false}, args::Tuple{Adjoint{Float64, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, Diagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, InfiniteArrays.InfStepRange{Int64, Int64}}}}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}})
@ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:298
[28310] _simplify(::typeof(*), ::Adjoint{Float64, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, ::Diagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, InfiniteArrays.InfStepRange{Int64, Int64}}}}, ::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, ::Vararg{Any})
@ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:296
--- the last 3 lines are repeated 1 more time ---
[28314] simplify
@ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:292 [inlined]
[28315] simplify
@ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:302 [inlined]
[28316] copy
@ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:307 [inlined]
[28317] materialize
@ ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:109 [inlined]
[28318] mul
@ ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:110 [inlined]
[28319] *(A::LazyArrays.ApplyArray{Float64, 2, typeof(*), Tuple{Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}, Adjoint{Float64, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, Diagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, InfiniteArrays.InfStepRange{Int64, Int64}}}}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}}, B::Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}})
@ ArrayLayouts ~/.julia/packages/ArrayLayouts/CNneg/src/ArrayLayouts.jl:194
[28320] afoldl(::typeof(*), ::LazyArrays.ApplyArray{Float64, 2, typeof(*), Tuple{Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}, Adjoint{Float64, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, Diagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, InfiniteArrays.InfStepRange{Int64, Int64}}}}}}, ::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, ::Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}})
@ Base ./operators.jl:549
[28321] *(::Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}, ::Adjoint{Float64, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, ::Diagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, InfiniteArrays.InfStepRange{Int64, Int64}}}}, ::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, ::Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}})
@ Base ./operators.jl:591
[28322] _most_simplify(#unused#::Val{true}, args::Tuple{Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}, QuasiArrays.QuasiAdjoint{Float64, Jacobi{Float64}}, Jacobi{Float64}, Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}})
@ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:297
[28323] _simplify(::typeof(*), ::Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}, ::QuasiArrays.QuasiAdjoint{Float64, Jacobi{Float64}}, ::Jacobi{Float64}, ::Vararg{Any})
@ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:296
[28324] simplify
@ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:292 [inlined]
[28325] simplify(M::ArrayLayouts.Mul{ContinuumArrays.AdjointBasisLayout{ClassicalOrthogonalPolynomials.NormalizedOPLayout{ClassicalOrthogonalPolynomials.OPLayout}}, ClassicalOrthogonalPolynomials.NormalizedOPLayout{ClassicalOrthogonalPolynomials.OPLayout}, QuasiArrays.QuasiAdjoint{Float64, Normalized{Float64, Jacobi{Float64}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}}, Normalized{Float64, Jacobi{Float64}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}})
@ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:302
[28326] copy
@ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:306 [inlined]
[28327] materialize
@ ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:109 [inlined]
[28328] mul
@ ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:110 [inlined]
[28329] *(A::QuasiArrays.QuasiAdjoint{Float64, Normalized{Float64, Jacobi{Float64}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}}, B::Normalized{Float64, Jacobi{Float64}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}})
@ QuasiArrays ~/.julia/packages/QuasiArrays/tF3vb/src/matmul.jl:23A similar but quite possibly unrelated issue also exists for finite sets of polynomials but here it happens even when not normalized:
julia> P = Jacobi(0,1)[:,1:10]
view(Jacobi(0.0, 1.0), Inclusion(-1.0..1.0 (Chebyshev)), 1:10) with eltype Float64
julia> P'P
ERROR: StackOverflowError:
Stacktrace:
[1] simplify(M::ArrayLayouts.Mul{LazyArrays.InvLayout{ArrayLayouts.BidiagonalLayout{LazyArrays.LazyLayout, LazyArrays.LazyLayout}}, BandedMatrices.BandedColumns{ArrayLayouts.OnesLayout}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}}})
@ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:302
[2] copy(M::ArrayLayouts.Mul{LazyArrays.InvLayout{ArrayLayouts.BidiagonalLayout{LazyArrays.LazyLayout, LazyArrays.LazyLayout}}, BandedMatrices.BandedColumns{ArrayLayouts.OnesLayout}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}}})
@ LazyBandedMatrices ~/.julia/packages/LazyBandedMatrices/ZaJ9N/src/LazyBandedMatrices.jl:899
[3] materialize(M::ArrayLayouts.Mul{LazyArrays.InvLayout{ArrayLayouts.BidiagonalLayout{LazyArrays.LazyLayout, LazyArrays.LazyLayout}}, BandedMatrices.BandedColumns{ArrayLayouts.OnesLayout}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}}})
@ ArrayLayouts ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:109
[4] mul(A::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, B::BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}})
@ ArrayLayouts ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:110
[5] _twoarg_simplify(#unused#::typeof(*), #unused#::Val{true}, a::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, b::BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}})
@ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:295
[6] _simplify(#unused#::typeof(*), a::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, b::BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}})
@ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:293
[7] simplify(::typeof(*), ::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, ::Vararg{Any})
@ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:292
--- the last 7 lines are repeated 4794 more times ---
[33566] simplify(M::ArrayLayouts.Mul{LazyArrays.InvLayout{ArrayLayouts.BidiagonalLayout{LazyArrays.LazyLayout, LazyArrays.LazyLayout}}, BandedMatrices.BandedColumns{ArrayLayouts.OnesLayout}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}}})
@ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:302
[33567] copy(M::ArrayLayouts.Mul{LazyArrays.InvLayout{ArrayLayouts.BidiagonalLayout{LazyArrays.LazyLayout, LazyArrays.LazyLayout}}, BandedMatrices.BandedColumns{ArrayLayouts.OnesLayout}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}}})
@ LazyBandedMatrices ~/.julia/packages/LazyBandedMatrices/ZaJ9N/src/LazyBandedMatrices.jl:899
[33568] materialize(M::ArrayLayouts.Mul{LazyArrays.InvLayout{ArrayLayouts.BidiagonalLayout{LazyArrays.LazyLayout, LazyArrays.LazyLayout}}, BandedMatrices.BandedColumns{ArrayLayouts.OnesLayout}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}}})
@ ArrayLayouts ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:109
[33569] mul(A::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, B::BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}})
@ ArrayLayouts ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:110
[33570] _twoarg_simplify(#unused#::typeof(*), #unused#::Val{true}, a::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, b::BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}})
@ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:295
[33571] _simplify
@ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:293 [inlined]
[33572] _tail_simplify
@ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:299 [inlined]
[33573] _most_simplify
@ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:298 [inlined]
[33574] _simplify
@ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:296 [inlined]
[33575] _tail_simplify
@ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:299 [inlined]
[33576] _most_simplify(#unused#::Val{false}, args::Tuple{Adjoint{Float64, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, Diagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, InfiniteArrays.InfStepRange{Int64, Int64}}}}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}}})
@ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:298
[33577] _simplify(::typeof(*), ::Adjoint{Float64, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, ::Diagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, InfiniteArrays.InfStepRange{Int64, Int64}}}}, ::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, ::Vararg{Any})
@ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:296
--- the last 3 lines are repeated 1 more time ---
[33581] simplify
@ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:292 [inlined]
[33582] simplify
@ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:302 [inlined]
[33583] copy
@ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:307 [inlined]
[33584] materialize
@ ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:109 [inlined]
[33585] mul
@ ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:110 [inlined]
[33586] *
@ ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:198 [inlined]
[33587] afoldl(op::typeof(*), a::LazyArrays.ApplyArray{Float64, 2, typeof(*), Tuple{Adjoint{Int64, BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}}}, Adjoint{Float64, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, Diagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, InfiniteArrays.InfStepRange{Int64, Int64}}}}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}}, bs::BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}})
@ Base ./operators.jl:548
[33588] *
@ ./operators.jl:591 [inlined]
[33589] _default_materialize(A::LazyArrays.Applied{LazyArrays.DefaultArrayApplyStyle, typeof(*), Tuple{Adjoint{Int64, BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}}}, QuasiArrays.QuasiAdjoint{Float64, Jacobi{Float64}}, Jacobi{Float64}, BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}}}}, #unused#::LazyArrays.Applied{LazyArrays.DefaultArrayApplyStyle, typeof(*), Tuple{Adjoint{Int64, BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}}}, QuasiArrays.QuasiAdjoint{Float64, Jacobi{Float64}}, Jacobi{Float64}, BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}}}})
@ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/lazyapplying.jl:87
[33590] copy
@ ~/.julia/packages/LazyArrays/NYra8/src/lazyapplying.jl:90 [inlined]
[33591] materialize
@ ~/.julia/packages/LazyArrays/NYra8/src/lazyapplying.jl:82 [inlined]
[33592] apply
@ ~/.julia/packages/LazyArrays/NYra8/src/lazyapplying.jl:79 [inlined]
[33593] copy(M::ArrayLayouts.Mul{ContinuumArrays.AdjointBasisLayout{ContinuumArrays.SubBasisLayout}, ContinuumArrays.SubBasisLayout, QuasiArrays.QuasiAdjoint{Float64, QuasiArrays.SubQuasiArray{Float64, 2, Jacobi{Float64}, Tuple{Inclusion{Float64, ChebyshevInterval{Float64}}, UnitRange{Int64}}, false}}, QuasiArrays.SubQuasiArray{Float64, 2, Jacobi{Float64}, Tuple{Inclusion{Float64, ChebyshevInterval{Float64}}, UnitRange{Int64}}, false}})
@ ContinuumArrays ~/.julia/packages/ContinuumArrays/sZfkh/src/bases/bases.jl:632
[33594] materialize
@ ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:109 [inlined]
[33595] mul
@ ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:110 [inlined]
[33596] *(A::QuasiArrays.QuasiAdjoint{Float64, QuasiArrays.SubQuasiArray{Float64, 2, Jacobi{Float64}, Tuple{Inclusion{Float64, ChebyshevInterval{Float64}}, UnitRange{Int64}}, false}}, B::QuasiArrays.SubQuasiArray{Float64, 2, Jacobi{Float64}, Tuple{Inclusion{Float64, ChebyshevInterval{Float64}}, UnitRange{Int64}}, false})
@ QuasiArrays ~/.julia/packages/QuasiArrays/tF3vb/src/matmul.jl:23
Make p-FEM more versatile
This should work:
julia> P = JacobiWeight(1,1) .* Jacobi(1,1)[:,1:n]
(1-x)^1 * (1+x)^1 on -1..1 .* view(Jacobi(1.0, 1.0), Inclusion(-1.0..1.0 (Chebyshev)), 1:10) with eltype Float64
julia> P'P
(QuasiArrays.QuasiAdjoint{Float64, QuasiArrays.BroadcastQuasiMatrix{Float64, typeof(*), Tuple{JacobiWeight{Int64}, QuasiArrays.SubQuasiArray{Float64, 2, Jacobi{Float64}, Tuple{Inclusion{Float64, ChebyshevInterval{Float64}}, UnitRange{Int64}}, false}}}}) * ((1-x)^1 * (1+x)^1 on -1..1 .* Jacobi(1.0, 1.0)) * (∞×10 BandedMatrix{Int64} with bandwidths (0, 0) with data 1×10 Ones{Int64} with indices OneToInf()×Base.OneTo(10))Instead we have to do:
julia> P = JacobiWeight(1,1) .* Jacobi(1,1)
(1-x)^1 * (1+x)^1 on -1..1 .* Jacobi(1.0, 1.0)
julia> (P'P)[1:n,1:n]
10×10 BandedMatrix{Float64} with bandwidths (2, 2):
1.06667 0.0 -0.228571 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
0.0 0.609524 5.55112e-17 -0.203175 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
-0.228571 0.0 0.457143 0.0 -0.17316 ⋅ ⋅ ⋅ ⋅ ⋅
⋅ -0.203175 0.0 0.369408 0.0 -0.149184 ⋅ ⋅ ⋅ ⋅
⋅ ⋅ -0.17316 -1.38778e-17 0.3108 2.77556e-17 -0.130536 ⋅ ⋅ ⋅
⋅ ⋅ ⋅ -0.149184 0.0 0.268531 -2.77556e-17 -0.115837 ⋅ ⋅
⋅ ⋅ ⋅ ⋅ -0.130536 -4.16334e-17 0.236501 5.55112e-17 -0.104025 ⋅
⋅ ⋅ ⋅ ⋅ ⋅ -0.115837 1.38778e-17 0.211352 0.0 -0.0943535
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -0.104025 0.0 0.191066 1.38778e-17
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -0.0943535 -1.38778e-17 0.174349
```Chebyshev() \ exp.(im*x) errors
julia> Chebyshev() \ exp.(im*x)
ERROR: InexactError: Float64(0.5443479390547952 + 0.8388595360647676im)
Stacktrace:
[1] Real
@ ./complex.jl:44 [inlined]
[2] convert
@ ./number.jl:7 [inlined]
[3] setindex!
@ ./array.jl:969 [inlined]
[4] _unsafe_copyto!(dest::Vector{Float64}, doffs::Int64, src::Vector{ComplexF64}, soffs::Int64, n::Int64)
@ Base ./array.jl:250
[5] unsafe_copyto!
@ ./array.jl:304 [inlined]
[6] _copyto_impl!
@ ./array.jl:327 [inlined]
[7] copyto!
@ ./array.jl:314 [inlined]
[8] copyto!
@ ./array.jl:339 [inlined]
[9] circcopy!(dest::Vector{Float64}, src::Vector{ComplexF64})
@ Base ./multidimensional.jl:1244
[10] copy1
@ ~/.julia/packages/AbstractFFTs/0uOAT/src/definitions.jl:50 [inlined]
[11] *
@ ~/.julia/packages/AbstractFFTs/0uOAT/src/definitions.jl:220 [inlined]
[12] \(a::ContinuumArrays.TransformFactorization{Float64, ChebyshevGrid{1, Float64}, FastTransforms.ChebyshevTransformPlan{Float64, 1, Vector{Int32}, false, 1, Int64}}, b::QuasiArrays.BroadcastQuasiVector{ComplexF64, typeof(exp), Tuple{ContinuumArrays.AffineQuasiVector{ComplexF64, ComplexF64, Inclusion{Float64, ChebyshevInterval{Float64}}, ComplexF64}}})
@ ContinuumArrays ~/.julia/packages/ContinuumArrays/2jOj3/src/plans.jl:26
[13] transform_ldiv
@ ~/.julia/packages/ContinuumArrays/2jOj3/src/bases/bases.jl:257 [inlined]
[14] transform_ldiv
@ ~/.julia/packages/ContinuumArrays/2jOj3/src/bases/bases.jl:258 [inlined]
[15] transform_ldiv_if_columns
@ ~/.julia/packages/ContinuumArrays/2jOj3/src/bases/bases.jl:289 [inlined]
[16] transform_ldiv_if_columns
@ ~/.julia/packages/ContinuumArrays/2jOj3/src/bases/bases.jl:290 [inlined]
[17] copy
@ ~/.julia/packages/ContinuumArrays/2jOj3/src/bases/bases.jl:291 [inlined]
[18] #materialize#13
@ ~/.julia/packages/ArrayLayouts/ph6IB/src/ldiv.jl:22 [inlined]
[19] materialize
@ ~/.julia/packages/ArrayLayouts/ph6IB/src/ldiv.jl:22 [inlined]
[20] #ldiv#20
@ ~/.julia/packages/ArrayLayouts/ph6IB/src/ldiv.jl:86 [inlined]
[21] ldiv
@ ~/.julia/packages/ArrayLayouts/ph6IB/src/ldiv.jl:86 [inlined]
[22] \
@ ~/.julia/packages/QuasiArrays/EK26C/src/matmul.jl:34 [inlined]
[23] adaptivetransform_ldiv(A::ChebyshevT{Float64}, f::QuasiArrays.BroadcastQuasiVector{ComplexF64, typeof(exp), Tuple{ContinuumArrays.AffineQuasiVector{ComplexF64, ComplexF64, Inclusion{Float64, ChebyshevInterval{Float64}}, ComplexF64}}})
@ ClassicalOrthogonalPolynomials ~/.julia/packages/ClassicalOrthogonalPolynomials/sveBE/src/adaptivetransform.jl:34
[24] transform_ldiv
@ ~/.julia/packages/ClassicalOrthogonalPolynomials/sveBE/src/adaptivetransform.jl:2 [inlined]
[25] transform_ldiv
@ ~/.julia/packages/ContinuumArrays/2jOj3/src/bases/bases.jl:258 [inlined]
[26] transform_ldiv_if_columns
@ ~/.julia/packages/ContinuumArrays/2jOj3/src/bases/bases.jl:289 [inlined]
[27] transform_ldiv_if_columns
@ ~/.julia/packages/ContinuumArrays/2jOj3/src/bases/bases.jl:290 [inlined]
[28] copy
@ ~/.julia/packages/ContinuumArrays/2jOj3/src/bases/bases.jl:291 [inlined]
[29] materialize(M::ArrayLayouts.Ldiv{ClassicalOrthogonalPolynomials.OPLayout, LazyArrays.BroadcastLayout{typeof(exp)}, ChebyshevT{Float64}, QuasiArrays.BroadcastQuasiVector{ComplexF64, typeof(exp), Tuple{ContinuumArrays.AffineQuasiVector{ComplexF64, ComplexF64, Inclusion{Float64, ChebyshevInterval{Float64}}, ComplexF64}}}}; kwds::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ ArrayLayouts ~/.julia/packages/ArrayLayouts/ph6IB/src/ldiv.jl:22
[30] materialize(M::ArrayLayouts.Ldiv{ClassicalOrthogonalPolynomials.OPLayout, LazyArrays.BroadcastLayout{typeof(exp)}, ChebyshevT{Float64}, QuasiArrays.BroadcastQuasiVector{ComplexF64, typeof(exp), Tuple{ContinuumArrays.AffineQuasiVector{ComplexF64, ComplexF64, Inclusion{Float64, ChebyshevInterval{Float64}}, ComplexF64}}}})
@ ArrayLayouts ~/.julia/packages/ArrayLayouts/ph6IB/src/ldiv.jl:22
[31] ldiv(A::ChebyshevT{Float64}, B::QuasiArrays.BroadcastQuasiVector{ComplexF64, typeof(exp), Tuple{ContinuumArrays.AffineQuasiVector{ComplexF64, ComplexF64, Inclusion{Float64, ChebyshevInterval{Float64}}, ComplexF64}}}; kwds::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ ArrayLayouts ~/.julia/packages/ArrayLayouts/ph6IB/src/ldiv.jl:86
[32] ldiv
@ ~/.julia/packages/ArrayLayouts/ph6IB/src/ldiv.jl:86 [inlined]
[33] \(A::ChebyshevT{Float64}, B::QuasiArrays.BroadcastQuasiVector{ComplexF64, typeof(exp), Tuple{ContinuumArrays.AffineQuasiVector{ComplexF64, ComplexF64, Inclusion{Float64, ChebyshevInterval{Float64}}, ComplexF64}}})
@ QuasiArrays ~/.julia/packages/QuasiArrays/EK26C/src/matmul.jl:34
[34] top-level scope
@ REPL[39]:1It's because the transform needs to take into account the type of RHS.
How to cite the package?
@dlfivefifty What's the preferred way to cite this package as well as others in its ecosystem such as https://github.com/JuliaApproximation/ContinuumArrays.jl and https://github.com/JuliaApproximation/QuasiArrays.jl ? I want to specifically refer to these packages in my thesis, hence the question.
Some considerations:
- Any papers associated with them to cite (either along with the repo or instead of)?
- Perhaps it would be good to consider either setting up a Zenodo DOI or leaving some notes in the readme file on which papers to cite instead?
What to do with axes of associated?
Sometimes the Associated OPs contain discrete points not present in the original OPs...
So what do we do for Hilbert transform? It's defined in terms of Associated OPs but we wouldn't want the support to increase.
Perhaps we should properly support restrictions? This would involve introducing QuasiSlice instead of assuming Inclusion is playing that role
Lanczos Jacobi matrices and `^`
Another bug related to JuliaApproximation/SemiclassicalOrthogonalPolynomials.jl#70
Jacobi matrices obtained from Lanczos don't appear to play nice with ^:
julia> using ClassicalOrthogonalPolynomials
julia> P = Legendre()
Legendre()
julia> x = axes(P,1)
Inclusion(-1.0..1.0 (Chebyshev))
julia> T = LanczosPolynomial(@.(2-x),P)
LanczosPolynomial{Float64} with weight with singularities LegendreWeight{Float64}
julia> X = jacobimatrix(T)
ℵ₀×ℵ₀ SymTridiagonal{Float64} Jacobi operator from Lanczos:
-0.166667 0.552771 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ …
0.552771 0.0212121 0.519377 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ 0.519377 0.00679908 0.508021 ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ 0.508021 0.00221202 0.504252 ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ 0.504252 0.000936283 0.502645 ⋅ ⋅ ⋅
⋮ ⋮ ⋱
julia> X^1
ERROR: MethodError: Cannot `convert` an object of type LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}} to an object of type ApplyArray{Float64, 2, typeof(*), Tuple{LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}, LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}}}
Closest candidates are:
convert(::Type{T}, ::Factorization) where T<:AbstractArray at ~/Documents/Backends/julia-1.8.2/share/julia/stdlib/v1.8/LinearAlgebra/src/factorization.jl:58
convert(::Type{T}, ::QuasiArrays.QuasiIteratorsMD.QuasiCartesianIndex{1}) where T at ~/.julia/packages/QuasiArrays/tF3vb/src/multidimensional.jl:131
convert(::Type{T}, ::T) where T<:AbstractArray at abstractarray.jl:16
...
Stacktrace:
[1] to_power_type(x::LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}})
@ Base ./intfuncs.jl:250
[2] power_by_squaring(x_::LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}, p::Int64)
@ Base ./intfuncs.jl:265
[3] _power_by_squaring
@ ~/.julia/packages/ArrayLayouts/4IG3b/src/mul.jl:337 [inlined]
[4] power_by_squaring
@ ~/.julia/packages/ArrayLayouts/4IG3b/src/mul.jl:345 [inlined]
[5] ^(A::LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}, p::Int64)
@ LinearAlgebra ~/Documents/Backends/julia-1.8.2/share/julia/stdlib/v1.8/LinearAlgebra/src/dense.jl:454
[6] literal_pow(f::typeof(^), x::LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}, #unused#::Val{1})
@ Base ./intfuncs.jl:340
[7] top-level scope
@ ~/Documents/Projects/ClassicalOrthogonalPolynomials.jl/test/exp.jl:267
julia> X^2
ERROR: MethodError: Cannot `convert` an object of type LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}} to an object of type ApplyArray{Float64, 2, typeof(*), Tuple{LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}, LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}}}
Closest candidates are:
convert(::Type{T}, ::Factorization) where T<:AbstractArray at ~/Documents/Backends/julia-1.8.2/share/julia/stdlib/v1.8/LinearAlgebra/src/factorization.jl:58
convert(::Type{T}, ::QuasiArrays.QuasiIteratorsMD.QuasiCartesianIndex{1}) where T at ~/.julia/packages/QuasiArrays/tF3vb/src/multidimensional.jl:131
convert(::Type{T}, ::T) where T<:AbstractArray at abstractarray.jl:16
...
Stacktrace:
[1] to_power_type(x::LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}})
@ Base ./intfuncs.jl:250
[2] power_by_squaring(x_::LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}, p::Int64)
@ Base ./intfuncs.jl:265
[3] _power_by_squaring
@ ~/.julia/packages/ArrayLayouts/4IG3b/src/mul.jl:337 [inlined]
[4] power_by_squaring
@ ~/.julia/packages/ArrayLayouts/4IG3b/src/mul.jl:345 [inlined]
[5] ^(A::LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}, p::Int64)
@ LinearAlgebra ~/Documents/Backends/julia-1.8.2/share/julia/stdlib/v1.8/LinearAlgebra/src/dense.jl:454
[6] literal_pow(f::typeof(^), x::LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}, #unused#::Val{2})
@ Base ./intfuncs.jl:340
[7] top-level scope
@ ~/Documents/Projects/ClassicalOrthogonalPolynomials.jl/test/exp.jl:268Multiplication and addition work fine:
julia> X*X
(ℵ₀×ℵ₀ SymTridiagonal{Float64} Jacobi operator from Lanczos) * (ℵ₀×ℵ₀ SymTridiagonal{Float64} Jacobi operator from Lanczos) with indices OneToInf()×OneToInf():
0.333333 -0.080403 0.287096 ⋅ ⋅ ⋅ ⋅ ⋅ …
-0.080403 0.575758 0.0145484 0.263854 ⋅ ⋅ ⋅ ⋅
0.287096 0.0145484 0.527884 0.00457783 0.256171 ⋅ ⋅ ⋅
⋅ 0.263854 0.00457783 0.512361 0.00158754 0.25346 ⋅ ⋅
⋅ ⋅ 0.256171 0.00158754 0.506923 0.000714292 0.252233 ⋅
⋮ ⋮ ⋱
julia> X+X
ℵ₀×ℵ₀ LazyBandedMatrices.SymTridiagonal{Float64, BroadcastVector{Float64, typeof(+), Tuple{ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}}, BroadcastVector{Float64, typeof(+), Tuple{ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}}} with indices OneToInf()×OneToInf():
-0.333333 1.10554 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ …
1.10554 0.0424242 1.03875 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ 1.03875 0.0135982 1.01604 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ 1.01604 0.00442404 1.0085 ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ 1.0085 0.00187257 1.00529 ⋅ ⋅ ⋅ ⋅
⋮ ⋮ ⋱ In the meantime, there are many workarounds. E.g.:
julia> X = Symmetric(BandedMatrix(0 => X.dv, 1 => X.ev))
ℵ₀×ℵ₀ Symmetric{Float64, BandedMatrix{Float64, ApplyArray{Float64, 2, typeof(*), Tuple{Matrix{Float64}, ApplyArray{Float64, 2, typeof(vcat), Tuple{ApplyArray{Float64, 2, typeof(hcat), Tuple{Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.ReshapedArray{Float64, 2, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, Tuple{Int64, Infinities.InfiniteCardinal{0}}, Tuple{}}}}, ApplyArray{Float64, 2, typeof(hcat), Tuple{Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.ReshapedArray{Float64, 2, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, Tuple{Int64, Infinities.InfiniteCardinal{0}}, Tuple{}}}}}}}}, InfiniteArrays.OneToInf{Int64}}} with indices OneToInf()×OneToInf():
-0.166667 0.552771 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ …
0.552771 0.0212121 0.519377 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ 0.519377 0.00679908 0.508021 ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ 0.508021 0.00221202 0.504252 ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ 0.504252 0.000936283 0.502645 ⋅ ⋅ ⋅
⋮ ⋮ ⋱
julia> X^1
ℵ₀×ℵ₀ Symmetric{Float64, BandedMatrix{Float64, ApplyArray{Float64, 2, typeof(*), Tuple{Matrix{Float64}, ApplyArray{Float64, 2, typeof(vcat), Tuple{ApplyArray{Float64, 2, typeof(hcat), Tuple{Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.ReshapedArray{Float64, 2, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, Tuple{Int64, Infinities.InfiniteCardinal{0}}, Tuple{}}}}, ApplyArray{Float64, 2, typeof(hcat), Tuple{Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.ReshapedArray{Float64, 2, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, Tuple{Int64, Infinities.InfiniteCardinal{0}}, Tuple{}}}}}}}}, InfiniteArrays.OneToInf{Int64}}} with indices OneToInf()×OneToInf():
-0.166667 0.552771 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ …
0.552771 0.0212121 0.519377 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ 0.519377 0.00679908 0.508021 ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ 0.508021 0.00221202 0.504252 ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ 0.504252 0.000936283 0.502645 ⋅ ⋅ ⋅
⋮ ⋮ ⋱
julia> X^2
ℵ₀×ℵ₀ Symmetric{Float64, ApplyArray{Float64, 2, typeof(*), Tuple{Symmetric{Float64, BandedMatrix{Float64, ApplyArray{Float64, 2, typeof(*), Tuple{Matrix{Float64}, ApplyArray{Float64, 2, typeof(vcat), Tuple{ApplyArray{Float64, 2, typeof(hcat), Tuple{Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.ReshapedArray{Float64, 2, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, Tuple{Int64, Infinities.InfiniteCardinal{0}}, Tuple{}}}}, ApplyArray{Float64, 2, typeof(hcat), Tuple{Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.ReshapedArray{Float64, 2, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, Tuple{Int64, Infinities.InfiniteCardinal{0}}, Tuple{}}}}}}}}, InfiniteArrays.OneToInf{Int64}}}, Symmetric{Float64, BandedMatrix{Float64, ApplyArray{Float64, 2, typeof(*), Tuple{Matrix{Float64}, ApplyArray{Float64, 2, typeof(vcat), Tuple{ApplyArray{Float64, 2, typeof(hcat), Tuple{Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.ReshapedArray{Float64, 2, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, Tuple{Int64, Infinities.InfiniteCardinal{0}}, Tuple{}}}}, ApplyArray{Float64, 2, typeof(hcat), Tuple{Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.ReshapedArray{Float64, 2, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, Tuple{Int64, Infinities.InfiniteCardinal{0}}, Tuple{}}}}}}}}, InfiniteArrays.OneToInf{Int64}}}}}} with indices OneToInf()×OneToInf():
0.333333 -0.080403 0.287096 ⋅ ⋅ ⋅ ⋅ ⋅ …
-0.080403 0.575758 0.0145484 0.263854 ⋅ ⋅ ⋅ ⋅
0.287096 0.0145484 0.527884 0.00457783 0.256171 ⋅ ⋅ ⋅
⋅ 0.263854 0.00457783 0.512361 0.00158754 0.25346 ⋅ ⋅
⋅ ⋅ 0.256171 0.00158754 0.506923 0.000714292 0.252233 ⋅
⋮ ⋮ ⋱ Ultraspherical{BigFloat} expansion broken
Solutions:
- Use GenericLinearAlgebra.jl
- Use Chebyshev->Ultraspherical transform
julia> C = Ultraspherical{BigFloat}(2)
Ultraspherical{BigFloat, Int64}
julia> x = axes(C,2)
OneToInf()
julia> C \ exp.(x)^C
julia> x = axes(C,1)
Inclusion(-1.0..1.0 (Chebyshev))
julia> C \ exp.(x)
ERROR: MethodError: no method matching eigen!(::SymTridiagonal{BigFloat, Vector{BigFloat}})
Closest candidates are:
eigen!(::StridedMatrix{T}; permute, scale, sortby) where T<:Union{Float32, Float64} at /Users/sheehanolver/Projects/julia-1.6/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/eigen.jl:148
eigen!(::StridedMatrix{T}; permute, scale, sortby) where T<:Union{ComplexF32, ComplexF64} at /Users/sheehanolver/Projects/julia-1.6/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/eigen.jl:171
eigen!(::StridedMatrix{T}, ::StridedMatrix{T}; sortby) where T<:Union{Float32, Float64} at /Users/sheehanolver/Projects/julia-1.6/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/eigen.jl:427
...
Stacktrace:
[1] eigen(A::SymTridiagonal{BigFloat, Vector{BigFloat}})
@ LinearAlgebra ~/Projects/julia-1.6/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/tridiag.jl:281
[2] eigen(A::LazyBandedMatrices.SymTridiagonal{BigFloat, FillArrays.Zeros{BigFloat, 1, Tuple{Base.OneTo{Int64}}}, Vector{BigFloat}})
@ LazyBandedMatrices ~/.julia/packages/LazyBandedMatrices/IX6Zy/src/tridiag.jl:723
[3] golubwelsch(X::LazyBandedMatrices.Tridiagonal{BigFloat, Vector{BigFloat}, FillArrays.Zeros{BigFloat, 1, Tuple{Base.OneTo{Int64}}}, Vector{BigFloat}})
@ ClassicalOrthogonalPolynomials ~/Projects/ClassicalOrthogonalPolynomials.jl/src/ClassicalOrthogonalPolynomials.jl:287
[4] golubwelsch(V::QuasiArrays.SubQuasiArray{BigFloat, 2, Ultraspherical{BigFloat, Int64}, Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}, Base.OneTo{Int64}}, false})
@ ClassicalOrthogonalPolynomials ~/Projects/ClassicalOrthogonalPolynomials.jl/src/ClassicalOrthogonalPolynomials.jl:292
[5] factorize(L::QuasiArrays.SubQuasiArray{BigFloat, 2, Ultraspherical{BigFloat, Int64}, Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}, Base.OneTo{Int64}}, false})
@ ClassicalOrthogonalPolynomials ~/Projects/ClassicalOrthogonalPolynomials.jl/src/ClassicalOrthogonalPolynomials.jl:304
[6] transform_ldiv(A::QuasiArrays.SubQuasiArray{BigFloat, 2, Ultraspherical{BigFloat, Int64}, Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}, Base.OneTo{Int64}}, false}, B::QuasiArrays.BroadcastQuasiVector{BigFloat, typeof(exp), Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}}}, #unused#::Tuple{Infinities.InfiniteCardinal{1}, Int64})
@ ContinuumArrays ~/Projects/ContinuumArrays.jl/src/bases/bases.jl:255
[7] transform_ldiv(A::QuasiArrays.SubQuasiArray{BigFloat, 2, Ultraspherical{BigFloat, Int64}, Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}, Base.OneTo{Int64}}, false}, B::QuasiArrays.BroadcastQuasiVector{BigFloat, typeof(exp), Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}}})
@ ContinuumArrays ~/Projects/ContinuumArrays.jl/src/bases/bases.jl:256
[8] copy(L::ArrayLayouts.Ldiv{ContinuumArrays.SubBasisLayout, LazyArrays.BroadcastLayout{typeof(exp)}, QuasiArrays.SubQuasiArray{BigFloat, 2, Ultraspherical{BigFloat, Int64}, Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}, Base.OneTo{Int64}}, false}, QuasiArrays.BroadcastQuasiVector{BigFloat, typeof(exp), Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}}}})
@ ContinuumArrays ~/Projects/ContinuumArrays.jl/src/bases/bases.jl:262
[9] materialize
@ ~/.julia/packages/ArrayLayouts/3lnt1/src/ldiv.jl:22 [inlined]
[10] ldiv
@ ~/.julia/packages/ArrayLayouts/3lnt1/src/ldiv.jl:86 [inlined]
[11] \
@ ~/Projects/QuasiArrays.jl/src/matmul.jl:34 [inlined]
[12] adaptivetransform_ldiv(A::Ultraspherical{BigFloat, Int64}, f::QuasiArrays.BroadcastQuasiVector{BigFloat, typeof(exp), Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}}})
@ ClassicalOrthogonalPolynomials ~/Projects/ClassicalOrthogonalPolynomials.jl/src/ClassicalOrthogonalPolynomials.jl:92
[13] transform_ldiv(A::Ultraspherical{BigFloat, Int64}, f::QuasiArrays.BroadcastQuasiVector{BigFloat, typeof(exp), Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}}}, #unused#::Tuple{Infinities.InfiniteCardinal{1}, Infinities.InfiniteCardinal{0}})
@ ClassicalOrthogonalPolynomials ~/Projects/ClassicalOrthogonalPolynomials.jl/src/ClassicalOrthogonalPolynomials.jl:62
[14] transform_ldiv(A::Ultraspherical{BigFloat, Int64}, B::QuasiArrays.BroadcastQuasiVector{BigFloat, typeof(exp), Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}}})
@ ContinuumArrays ~/Projects/ContinuumArrays.jl/src/bases/bases.jl:256
[15] copy(L::ArrayLayouts.Ldiv{ContinuumArrays.BasisLayout, LazyArrays.BroadcastLayout{typeof(exp)}, Ultraspherical{BigFloat, Int64}, QuasiArrays.BroadcastQuasiVector{BigFloat, typeof(exp), Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}}}})
@ ContinuumArrays ~/Projects/ContinuumArrays.jl/src/bases/bases.jl:262
[16] materialize(M::ArrayLayouts.Ldiv{ContinuumArrays.BasisLayout, LazyArrays.BroadcastLayout{typeof(exp)}, Ultraspherical{BigFloat, Int64}, QuasiArrays.BroadcastQuasiVector{BigFloat, typeof(exp), Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}}}})
@ ArrayLayouts ~/.julia/packages/ArrayLayouts/3lnt1/src/ldiv.jl:22
[17] ldiv
@ ~/.julia/packages/ArrayLayouts/3lnt1/src/ldiv.jl:86 [inlined]
[18] \(A::Ultraspherical{BigFloat, Int64}, B::QuasiArrays.BroadcastQuasiVector{BigFloat, typeof(exp), Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}}})
@ QuasiArrays ~/Projects/QuasiArrays.jl/src/matmul.jl:34
[19] top-level scope
@ REPL[130]:1Conversion operator from LanczosPolynomial
x = Inclusion(-1.0..1)
a = 1.5
ϕ = x.^4 - (a^2 + 1)*x.^2 .+ a^2
Pϕ = Normalized(LanczosPolynomial(ϕ))
P = Normalized(Legendre())
Cϕ = Pϕ\P
MethodError: copy(::ArrayLayouts.Ldiv{ContinuumArrays.MappedBasisLayout, ContinuumArrays.BasisLayout, QuasiArrays.SubQuasiArray{Float64, 2, Legendre{Float64}, Tuple{ContinuumArrays.AffineMap{Float64, Inclusion{Float64, IntervalSets.ClosedInterval{Float64}}, Inclusion{Float64, ChebyshevInterval{Float64}}}, Base.Slice{InfiniteArrays.OneToInf{Int64}}}, false}, Legendre{Float64}}) is ambiguous. Candidates:
copy(P::ArrayLayouts.Ldiv{var"#s13", var"#s12", AType, BType} where {var"#s13"<:Union{ContinuumArrays.MappedBasisLayout, ContinuumArrays.MappedWeightedBasisLayout}, var"#s12"<:LazyArrays.AbstractLazyLayout, AType, BType}) in ContinuumArrays at C:\Users\mfaso\.julia\packages\ContinuumArrays\gK7Yu\src\bases\bases.jl:88
copy(P::ArrayLayouts.Ldiv{var"#s13", var"#s12", AType, BType} where {var"#s13"<:ContinuumArrays.AbstractBasisLayout, var"#s12"<:ContinuumArrays.AbstractBasisLayout, AType, BType}) in ContinuumArrays at C:\Users\mfaso\.julia\packages\ContinuumArrays\gK7Yu\src\bases\bases.jl:71
Possible fix, define
copy(::ArrayLayouts.Ldiv{var"#s13", var"#s12", AType, BType} where {var"#s13"<:Union{ContinuumArrays.MappedBasisLayout, ContinuumArrays.MappedWeightedBasisLayout}, var"#s12"<:ContinuumArrays.AbstractBasisLayout, AType, BType})
Stacktrace:
[1] materialize
@ ~\.julia\packages\ArrayLayouts\pPGEV\src\ldiv.jl:22 [inlined]
[2] ldiv
@ ~\.julia\packages\ArrayLayouts\pPGEV\src\ldiv.jl:86 [inlined]
[3] \
@ ~\.julia\packages\QuasiArrays\z4zwW\src\matmul.jl:34 [inlined]
[4] _copy_ldiv_ldiv
@ ~\.julia\packages\LazyArrays\Wldxo\src\linalg\inv.jl:105 [inlined]
[5] copy
@ ~\.julia\packages\LazyArrays\Wldxo\src\linalg\inv.jl:107 [inlined]
[6] copy
@ ~\.julia\packages\ClassicalOrthogonalPolynomials\63Ert\src\normalized.jl:115 [inlined]
[7] materialize(M::ArrayLayouts.Ldiv{ClassicalOrthogonalPolynomials.NormalizedBasisLayout{ContinuumArrays.MappedBasisLayout}, ContinuumArrays.BasisLayout, Normalized{Float64, QuasiArrays.SubQuasiArray{Float64, 2, Legendre{Float64}, Tuple{ContinuumArrays.AffineMap{Float64, Inclusion{Float64, IntervalSets.ClosedInterval{Float64}}, Inclusion{Float64, ChebyshevInterval{Float64}}}, Base.Slice{InfiniteArrays.OneToInf{Int64}}}, false}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}, Legendre{Float64}})
@ ArrayLayouts ~\.julia\packages\ArrayLayouts\pPGEV\src\ldiv.jl:22
[8] ldiv
@ ~\.julia\packages\ArrayLayouts\pPGEV\src\ldiv.jl:86 [inlined]
[9] \
@ ~\.julia\packages\QuasiArrays\z4zwW\src\matmul.jl:34 [inlined]
[10] \(Q::LanczosPolynomial{Float64, QuasiArrays.ApplyQuasiVector{Float64, typeof(*), Tuple{Normalized{Float64, QuasiArrays.SubQuasiArray{Float64, 2, Legendre{Float64}, Tuple{ContinuumArrays.AffineMap{Float64, Inclusion{Float64, IntervalSets.ClosedInterval{Float64}}, Inclusion{Float64, ChebyshevInterval{Float64}}}, Base.Slice{InfiniteArrays.OneToInf{Int64}}}, false}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}, LazyArrays.CachedArray{Float64, 1, Vector{Float64}, FillArrays.Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, Normalized{Float64, QuasiArrays.SubQuasiArray{Float64, 2, Legendre{Float64}, Tuple{ContinuumArrays.AffineMap{Float64, Inclusion{Float64, IntervalSets.ClosedInterval{Float64}}, Inclusion{Float64, ChebyshevInterval{Float64}}}, Base.Slice{InfiniteArrays.OneToInf{Int64}}}, false}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}}, A::Legendre{Float64})
@ ClassicalOrthogonalPolynomials ~\.julia\packages\ClassicalOrthogonalPolynomials\63Ert\src\lanczos.jl:218
[11] copy
@ ~\.julia\packages\ClassicalOrthogonalPolynomials\63Ert\src\normalized.jl:110 [inlined]
[12] materialize(M::ArrayLayouts.Ldiv{ClassicalOrthogonalPolynomials.NormalizedBasisLayout{ClassicalOrthogonalPolynomials.LanczosLayout}, ClassicalOrthogonalPolynomials.NormalizedBasisLayout{ContinuumArrays.BasisLayout}, Normalized{Float64, LanczosPolynomial{Float64, QuasiArrays.ApplyQuasiVector{Float64, typeof(*), Tuple{Normalized{Float64, QuasiArrays.SubQuasiArray{Float64, 2, Legendre{Float64}, Tuple{ContinuumArrays.AffineMap{Float64, Inclusion{Float64, IntervalSets.ClosedInterval{Float64}}, Inclusion{Float64, ChebyshevInterval{Float64}}}, Base.Slice{InfiniteArrays.OneToInf{Int64}}}, false}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}, LazyArrays.CachedArray{Float64, 1, Vector{Float64}, FillArrays.Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, Normalized{Float64, QuasiArrays.SubQuasiArray{Float64, 2, Legendre{Float64}, Tuple{ContinuumArrays.AffineMap{Float64, Inclusion{Float64, IntervalSets.ClosedInterval{Float64}}, Inclusion{Float64, ChebyshevInterval{Float64}}}, Base.Slice{InfiniteArrays.OneToInf{Int64}}}, false}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}, Normalized{Float64, Legendre{Float64}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}})
@ ArrayLayouts ~\.julia\packages\ArrayLayouts\pPGEV\src\ldiv.jl:22
[13] ldiv
@ ~\.julia\packages\ArrayLayouts\pPGEV\src\ldiv.jl:86 [inlined]
[14] \(A::Normalized{Float64, LanczosPolynomial{Float64, QuasiArrays.ApplyQuasiVector{Float64, typeof(*), Tuple{Normalized{Float64, QuasiArrays.SubQuasiArray{Float64, 2, Legendre{Float64}, Tuple{ContinuumArrays.AffineMap{Float64, Inclusion{Float64, IntervalSets.ClosedInterval{Float64}}, Inclusion{Float64, ChebyshevInterval{Float64}}}, Base.Slice{InfiniteArrays.OneToInf{Int64}}}, false}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}, LazyArrays.CachedArray{Float64, 1, Vector{Float64}, FillArrays.Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, Normalized{Float64, QuasiArrays.SubQuasiArray{Float64, 2, Legendre{Float64}, Tuple{ContinuumArrays.AffineMap{Float64, Inclusion{Float64, IntervalSets.ClosedInterval{Float64}}, Inclusion{Float64, ChebyshevInterval{Float64}}}, Base.Slice{InfiniteArrays.OneToInf{Int64}}}, false}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}, B::Normalized{Float64, Legendre{Float64}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}})
@ QuasiArrays ~\.julia\packages\QuasiArrays\z4zwW\src\matmul.jl:34
[15] top-level scope
@ In[30]:6
[16] eval
@ .\boot.jl:360 [inlined]
[17] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base .\loading.jl:1094Unify LegendreWeight and QuasiOnes
They are the same thing though the Weight interface might be an issue
Cholesky Jacobi matrices are insanely slow
The culprit is that UX is a lazy multiplication so the following line doesn't do anything:
Here is the slow code:
julia> r = 0.5
0.5
julia> lmin, lmax = (1-sqrt(r))^2, (1+sqrt(r))^2
(0.08578643762690492, 2.914213562373095)
julia> U = chebyshevu(lmin..lmax);
julia> Q = OrthogonalPolynomial(x -> 1/(2π) * sqrt((lmax-x)*(x-lmin))/(x*r), U);
julia> J = jacobimatrix(Q);
julia> @time J[1:1000,1:1000]
36.732166 seconds (178.69 k allocations: 1.624 GiB, 0.29% gc time)
1000×1000 BandedMatrix{Float64} with bandwidths (1, 1):
1.0 0.707107 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
0.707107 1.5 0.707107 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ 0.707107 1.5 0.707107 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ 0.707107 1.5 0.707107 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ 0.707107 1.5 0.707107 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ 0.707107 1.5 0.707107 ⋅ ⋅ ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋮ ⋮ ⋮ ⋱ ⋮ ⋮
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.707107 1.5 0.707107 ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.707107 1.5 0.707107 ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.707107 1.5 0.707107 ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.707107 1.5 0.707107
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.707107 1.5UndefVarError in docs
What's going on in the docs that almost all lines of code evaluate to UndefVarError?
Laguerre
Is this package expected to, at some point, include the Laguerre polynomials (for evaluation), or are they already in some other package?
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
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.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.