JuliaApproximation/ClassicalOrthogonalPolynomials.jl

Can't compute scalar times squared jacobmatrix

Closed this issue · 3 comments

DanielVandH commented

Expressions like 1jacobimatrix(P)^2 seem to never complete. Not sure if it's because of something in this package or not - tried getting a specific matrix using just LazyArrays or LazyBandedMatrices but can't seem to get a better MWE of this.

julia> using ClassicalOrthogonalPolynomials

julia> P = Jacobi(0.2, 0.5);

julia> X = jacobimatrix(P);

julia> 1X; # works fine

julia> 2X # works fine
ℵ₀×ℵ₀ LazyBandedMatrices.Tridiagonal{Float64, LazyArrays.ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Int64, Int64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}, Int64}}}}, LazyArrays.ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{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}}}}}, Int64}}} with indices OneToInf()×OneToInf():
 0.222222  0.720721    ⋅          ⋅           ⋅          …
 1.48148   0.0330969  0.821202    ⋅           ⋅
  ⋅        1.24209    0.0133376  0.868385     ⋅
  ⋅         ⋅         1.16261    0.00720535  0.895841
  ⋅         ⋅          ⋅         1.12256     0.00451176
  ⋅         ⋅          ⋅          ⋅          1.09837     …
  ⋅         ⋅          ⋅          ⋅           ⋅
  ⋅         ⋅          ⋅          ⋅           ⋅
  ⋅         ⋅          ⋅          ⋅           ⋅
  ⋅         ⋅          ⋅          ⋅           ⋅
  ⋅         ⋅          ⋅          ⋅           ⋅          …
  ⋅         ⋅          ⋅          ⋅           ⋅
  ⋅         ⋅          ⋅          ⋅           ⋅
 ⋮                                                       ⋱

julia> 1X^2 # never ends
ERROR: InterruptException:
Stacktrace:
  [1] cache_filldata!(::LazyArrays.CachedArray{…}, ::UnitRange{…}, ::Base.OneTo{…})
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:222
  [2] resizedata!(::ArrayLayouts.DenseColumnMajor, ::ArrayLayouts.ZerosLayout, ::LazyArrays.CachedArray{…}, ::Int64, ::Int64)
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:263
  [3] resizedata!
    @ C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:218 [inlined]
  [4] setindex!
    @ C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:82 [inlined]
  [5] _setindex!
    @ .\abstractarray.jl:1426 [inlined]
  [6] setindex!
    @ .\abstractarray.jl:1396 [inlined]
  [7] macro expansion
    @ .\broadcast.jl:1004 [inlined]
  [8] macro expansion
    @ .\simdloop.jl:77 [inlined]
  [9] copyto!
    @ .\broadcast.jl:1003 [inlined]
 [10] copyto!
    @ .\broadcast.jl:956 [inlined]
 [11] copy
    @ .\broadcast.jl:928 [inlined]
 [12] materialize
    @ .\broadcast.jl:903 [inlined]
 [13] broadcast(::typeof(*), ::Int64, ::LazyBandedMatrices.Tridiagonal{…})
    @ Base.Broadcast .\broadcast.jl:841
 [14] broadcasted(::LazyArrays.LazyArrayStyle{…}, ::typeof(*), a::Int64, b::LazyArrays.ApplyArray{…})
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\linalg\mul.jl:284
 [15] broadcasted
    @ .\broadcast.jl:1347 [inlined]
 [16] broadcast_preserving_zero_d
    @ .\broadcast.jl:891 [inlined]
 [17] *(A::Int64, B::LazyArrays.ApplyArray{Float64, 2, typeof(*), Tuple{…}})
    @ Base .\arraymath.jl:21
 [18] top-level scope
    @ REPL[8]:1
Some type information was truncated. Use `show(err)` to see complete types.

julia> typeof(X^2)
LazyArrays.ApplyArray{Float64, 2, typeof(*), Tuple{LazyBandedMatrices.Tridiagonal{Float64, LazyArrays.ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Int64, Int64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, LazyArrays.ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{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}}}}}}, LazyBandedMatrices.Tridiagonal{Float64, LazyArrays.ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Int64, Int64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, LazyArrays.ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{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}}}}}}}}

julia> typeof(X)
LazyBandedMatrices.Tridiagonal{Float64, LazyArrays.ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Int64, Int64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, LazyArrays.ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{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}}}}}}
DanielVandH commented

Might be an issue with broadcasting somehow, since e.g.

julia> 1 .+ X

also never completes.

dlfivefifty commented

Yes I started debugging and looks like 1 .* X is the MWE. It should return:

julia> BroadcastArray(*, 1, X)
(Int64) .* (ℵ₀×ℵ₀ LazyBandedMatrices.Tridiagonal{Float64, ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, BroadcastVector{Float64, typeof(/), Tuple{BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Int64, Int64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, BroadcastVector{Float64, typeof(/), Tuple{Float64, BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, BroadcastVector{Float64, typeof(/), Tuple{BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}} with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf():
 0.111111  0.36036     ⋅          …  
 0.740741  0.0165485  0.410601       
  ⋅        0.621047   0.00666878     
  ⋅         ⋅         0.581304       
  ⋅         ⋅          ⋅             
  ⋅         ⋅          ⋅          …  
  ⋅         ⋅          ⋅             
  ⋅         ⋅          ⋅             
  ⋅         ⋅          ⋅             
  ⋅         ⋅          ⋅             
 ⋮                                ⋱
dlfivefifty commented

It looks like the issue is:

julia> Base.BroadcastStyle(typeof(X))
LinearAlgebra.StructuredMatrixStyle{Tridiagonal{Float64, ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, BroadcastVector{Float64, typeof(/), Tuple{Float64, BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}}}()

This should be returning some sort of LazyArrayStyle.

Unfortunately it's currently overloaded in LazyBandedMatrices.jl so will have to wait until LazyArrays v2.0 to avoid having to fix it twice.

  • 👍1