MultivariateOrthogonalPolynomials.jl
MultivariateOrthogonalPolynomials.jl copied to clipboard
JacobiTriangle() expansions failing in strange way
There is some strange behavior when expanding functions on triangles, apparently when it comes to globally defined variables? Here is a minimal example:
julia> using MultivariateOrthogonalPolynomials
julia> T = JacobiTriangle()
JacobiTriangle(0, 0, 0)
julia> xy = axes(T,1)
Inclusion(UnitSimplex(Val(2)))
julia> x, y = first.(xy), last.(xy)
(first.(Inclusion(UnitSimplex(Val(2)))), last.(Inclusion(UnitSimplex(Val(2)))))
julia> α = 2.
2.0
julia> K(x,y) = α*x*y
K (generic function with 1 method)
julia> K(α,x,y) = α*x*y
K (generic function with 2 methods)
julia> T \ (@. α*x*y) # this works
ℵ₀-blocked ℵ₀-element LazyBandedMatrices.DiagTrav{Float64, 2, LazyArrays.ApplyArray{Float64, 2, typeof(Base.setindex), Tuple{FillArrays.Zeros{Float64, 2, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Matrix{Float64}, Base.OneTo{Int64}, Base.OneTo{Int64}}}}:
0.1666666666666667
───────────────────────
0.06666666666666667
0.19999999999999998
───────────────────────
-0.09999999999999999
0.19999999999999993
-6.699811370772778e-18
───────────────────────
0.0
-6.148984834775253e-18
-1.6452652161292338e-18
⋮
julia> T \ (@. K(α,x,y)) # this works
ℵ₀-blocked ℵ₀-element LazyBandedMatrices.DiagTrav{Float64, 2, LazyArrays.ApplyArray{Float64, 2, typeof(Base.setindex), Tuple{FillArrays.Zeros{Float64, 2, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Matrix{Float64}, Base.OneTo{Int64}, Base.OneTo{Int64}}}}:
0.1666666666666667
───────────────────────
0.06666666666666667
0.19999999999999998
───────────────────────
-0.09999999999999999
0.19999999999999993
-6.699811370772778e-18
───────────────────────
0.0
-6.148984834775253e-18
-1.6452652161292338e-18
⋮
julia> T \ (@. K(x,y)) # this fails
ERROR: MethodError: no method matching *(::FastTransforms.FTPlan{Float64, 2, FastTransforms.TRIANGLEANALYSIS}, ::Matrix{Any})
Closest candidates are:
*(::Any, ::Any, ::Any, ::Any...) at operators.jl:591
*(::Union{SparseArrays.AbstractSparseMatrixCSC{TA, Ti}, SubArray{TA, 2, <:SparseArrays.AbstractSparseMatrixCSC{TA, Ti}, Tuple{Base.Slice{Base.OneTo{Int64}}, I}} where I<:AbstractUnitRange} where Ti, ::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, StridedMatrix, BitMatrix}) where TA at ~/Documents/Backends/julia-1.8.5/share/julia/stdlib/v1.8/SparseArrays/src/linalg.jl:52
*(::LinearAlgebra.Hermitian{<:Any, <:ArrayLayouts.LayoutMatrix}, ::AbstractMatrix) at ~/.julia/packages/ArrayLayouts/4IG3b/src/mul.jl:263
...
Stacktrace:
[1] *(T::MultivariateOrthogonalPolynomials.TriPlan{Float64}, F::Matrix{Any})
@ MultivariateOrthogonalPolynomials ~/.julia/packages/MultivariateOrthogonalPolynomials/tNdGo/src/triangle.jl:661
[2] \(a::ContinuumArrays.TransformFactorization{Any, Matrix{StaticArraysCore.SVector{2, Float64}}, MultivariateOrthogonalPolynomials.TriPlan{Float64}}, b::QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}})
@ ContinuumArrays ~/.julia/packages/ContinuumArrays/pWRm2/src/plans.jl:26
[3] transform_ldiv(V::QuasiArrays.SubQuasiArray{Float64, 2, JacobiTriangle{Float64, Int64}, Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, BlockArrays.BlockSlice{BlockArrays.BlockRange{1, Tuple{Base.OneTo{Int64}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, Base.OneTo{Int64}}}}}, false}, B::QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}}, #unused#::Tuple{Infinities.InfiniteCardinal{1}, Int64})
@ HarmonicOrthogonalPolynomials ~/.julia/packages/HarmonicOrthogonalPolynomials/51vXS/src/multivariateops.jl:91
[4] transform_ldiv(A::QuasiArrays.SubQuasiArray{Float64, 2, JacobiTriangle{Float64, Int64}, Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, BlockArrays.BlockSlice{BlockArrays.BlockRange{1, Tuple{Base.OneTo{Int64}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, Base.OneTo{Int64}}}}}, false}, B::QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}})
@ ContinuumArrays ~/.julia/packages/ContinuumArrays/pWRm2/src/bases/bases.jl:258
[5] transform_ldiv_if_columns(L::ArrayLayouts.Ldiv{ContinuumArrays.SubBasisLayout, LazyArrays.BroadcastLayout{typeof(K)}, QuasiArrays.SubQuasiArray{Float64, 2, JacobiTriangle{Float64, Int64}, Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, BlockArrays.BlockSlice{BlockArrays.BlockRange{1, Tuple{Base.OneTo{Int64}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, Base.OneTo{Int64}}}}}, false}, QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}}}, #unused#::Base.OneTo{Int64})
@ ContinuumArrays ~/.julia/packages/ContinuumArrays/pWRm2/src/bases/bases.jl:289
[6] transform_ldiv_if_columns(L::ArrayLayouts.Ldiv{ContinuumArrays.SubBasisLayout, LazyArrays.BroadcastLayout{typeof(K)}, QuasiArrays.SubQuasiArray{Float64, 2, JacobiTriangle{Float64, Int64}, Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, BlockArrays.BlockSlice{BlockArrays.BlockRange{1, Tuple{Base.OneTo{Int64}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, Base.OneTo{Int64}}}}}, false}, QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}}})
@ ContinuumArrays ~/.julia/packages/ContinuumArrays/pWRm2/src/bases/bases.jl:290
[7] copy(L::ArrayLayouts.Ldiv{ContinuumArrays.SubBasisLayout, LazyArrays.BroadcastLayout{typeof(K)}, QuasiArrays.SubQuasiArray{Float64, 2, JacobiTriangle{Float64, Int64}, Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, BlockArrays.BlockSlice{BlockArrays.BlockRange{1, Tuple{Base.OneTo{Int64}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, Base.OneTo{Int64}}}}}, false}, QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}}})
@ ContinuumArrays ~/.julia/packages/ContinuumArrays/pWRm2/src/bases/bases.jl:291
[8] materialize(M::ArrayLayouts.Ldiv{ContinuumArrays.SubBasisLayout, LazyArrays.BroadcastLayout{typeof(K)}, QuasiArrays.SubQuasiArray{Float64, 2, JacobiTriangle{Float64, Int64}, Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, BlockArrays.BlockSlice{BlockArrays.BlockRange{1, Tuple{Base.OneTo{Int64}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, Base.OneTo{Int64}}}}}, false}, QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}}}; kwds::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ ArrayLayouts ~/.julia/packages/ArrayLayouts/4IG3b/src/ldiv.jl:22
[9] materialize
@ ~/.julia/packages/ArrayLayouts/4IG3b/src/ldiv.jl:22 [inlined]
[10] #ldiv#18
@ ~/.julia/packages/ArrayLayouts/4IG3b/src/ldiv.jl:86 [inlined]
[11] ldiv
@ ~/.julia/packages/ArrayLayouts/4IG3b/src/ldiv.jl:86 [inlined]
[12] \
@ ~/.julia/packages/QuasiArrays/tF3vb/src/matmul.jl:34 [inlined]
[13] adaptivetransform_ldiv(A::JacobiTriangle{Float64, Int64}, f::QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}})
@ ClassicalOrthogonalPolynomials ~/.julia/packages/ClassicalOrthogonalPolynomials/ix48P/src/adaptivetransform.jl:35
[14] transform_ldiv(A::JacobiTriangle{Float64, Int64}, f::QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}}, #unused#::Tuple{Infinities.InfiniteCardinal{1}, Infinities.InfiniteCardinal{0}})
@ ClassicalOrthogonalPolynomials ~/.julia/packages/ClassicalOrthogonalPolynomials/ix48P/src/adaptivetransform.jl:2
[15] transform_ldiv(A::JacobiTriangle{Float64, Int64}, B::QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}})
@ ContinuumArrays ~/.julia/packages/ContinuumArrays/pWRm2/src/bases/bases.jl:258
[16] transform_ldiv_if_columns(L::ArrayLayouts.Ldiv{HarmonicOrthogonalPolynomials.MultivariateOPLayout{2}, LazyArrays.BroadcastLayout{typeof(K)}, JacobiTriangle{Float64, Int64}, QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}}}, #unused#::Base.OneTo{Int64})
@ ContinuumArrays ~/.julia/packages/ContinuumArrays/pWRm2/src/bases/bases.jl:289
[17] transform_ldiv_if_columns(L::ArrayLayouts.Ldiv{HarmonicOrthogonalPolynomials.MultivariateOPLayout{2}, LazyArrays.BroadcastLayout{typeof(K)}, JacobiTriangle{Float64, Int64}, QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}}})
@ ContinuumArrays ~/.julia/packages/ContinuumArrays/pWRm2/src/bases/bases.jl:290
[18] copy(L::ArrayLayouts.Ldiv{HarmonicOrthogonalPolynomials.MultivariateOPLayout{2}, LazyArrays.BroadcastLayout{typeof(K)}, JacobiTriangle{Float64, Int64}, QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}}})
@ ContinuumArrays ~/.julia/packages/ContinuumArrays/pWRm2/src/bases/bases.jl:291
[19] #materialize#11
@ ~/.julia/packages/ArrayLayouts/4IG3b/src/ldiv.jl:22 [inlined]
[20] materialize(M::ArrayLayouts.Ldiv{HarmonicOrthogonalPolynomials.MultivariateOPLayout{2}, LazyArrays.BroadcastLayout{typeof(K)}, JacobiTriangle{Float64, Int64}, QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}}})
@ ArrayLayouts ~/.julia/packages/ArrayLayouts/4IG3b/src/ldiv.jl:22
[21] ldiv(A::JacobiTriangle{Float64, Int64}, B::QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}}; kwds::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ ArrayLayouts ~/.julia/packages/ArrayLayouts/4IG3b/src/ldiv.jl:86
[22] ldiv
@ ~/.julia/packages/ArrayLayouts/4IG3b/src/ldiv.jl:86 [inlined]
[23] \(A::JacobiTriangle{Float64, Int64}, B::QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}})
@ QuasiArrays ~/.julia/packages/QuasiArrays/tF3vb/src/matmul.jl:34
[24] top-level scope
@ ~/Documents/Projects/SparseVolterraExamples.jl/examples/Example 0 - Some linear Volterra integral equation examples.jl:82
Any idea what is happening here?
This follow from the type instability of global variable. From your example
julia> typeof(K.(α,x,y))
BroadcastQuasiVector{Float64, typeof(K), Tuple{Float64, BroadcastQuasiVector{Float64, typeof(first), Tuple{Inclusion{SVector{2, Float64}, EuclideanUnitSimplex{2, Float64, :closed}}}}, BroadcastQuasiVector{Float64, typeof(last), Tuple{Inclusion{SVector{2, Float64}, EuclideanUnitSimplex{2, Float64, :closed}}}}}}
whereas
julia> typeof(K.(x,y))
BroadcastQuasiVector{Any, typeof(K), Tuple{BroadcastQuasiVector{Float64, typeof(first), Tuple{Inclusion{SVector{2, Float64}, EuclideanUnitSimplex{2, Float64, :closed}}}}, BroadcastQuasiVector{Float64, typeof(last), Tuple{Inclusion{SVector{2, Float64}, EuclideanUnitSimplex{2, Float64, :closed}}}}}}
Note the Any
in the second line of code. However if α is declare as a constant
const α =2.0
it works
julia> T \ (@. K(x,y))
ℵ₀-blocked ℵ₀-element LazyBandedMatrices.DiagTrav{Float64, 2, LazyArrays.ApplyArray{Float64, 2, typeof(Base.setindex), Tuple{FillArrays.Zeros{Float64, 2, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Matrix{Float64}, Base.OneTo{Int64}, Base.OneTo{Int64}}}}:
0.1666666666666667
───────────────────────
0.06666666666666667
0.19999999999999998
───────────────────────
-0.09999999999999999
0.19999999999999993
-6.699811370772778e-18
───────────────────────
0.0
-6.148984834775253e-18
-1.6452652161292338e-18
0.0
───────────────────────
⋮
Note I fixed essentially the same bug back in August: https://github.com/JuliaApproximation/ClassicalOrthogonalPolynomials.jl/commit/2e33f50b7326706775597c75f284eb57e3dc0597
So I suspect something similar needs to be done here, just need to step through the transforms and see where the promotion to Any
happens.