Error when expanding one triangle basis function in another

[MikaelSlevinsky]

I thought this would work, but there's a type issue

julia> n, k, a, b, c = 5, 3, 0.123, 0.456, 0.789
(5, 3, 0.123, 0.456, 0.789)

julia> P = JacobiTriangle()
JacobiTriangle(0, 0, 0)

julia> Pabc = JacobiTriangle(a, b, c)
JacobiTriangle(0.123, 0.456, 0.789)

julia> f = expand(P, splat((x,y) -> Pabc[SVector(x, y), Block(n+1)[k+1]]))
ERROR: MethodError: no method matching *(::FastTransforms.FTPlan{Float64, 2, FastTransforms.TRIANGLEANALYSIS}, ::Matrix{Any})

Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...)
   @ Base operators.jl:587
  *(::FastTransforms.IChebyshevTransformPlan{T, kind, K, false, N}, ::AbstractArray{T, N}) where {T, kind, K, N}
   @ FastTransforms ~/.julia/dev/FastTransforms/src/chebyshevtransform.jl:325
  *(::FillArrays.AbstractZeros{T, 2} where T, ::AbstractMatrix)
   @ FillArrays ~/.julia/packages/FillArrays/Ksvco/src/fillalgebra.jl:97
  ...

Stacktrace:
  [1] *(T::MultivariateOrthogonalPolynomials.TriPlan{Float64}, F::Matrix{Any})
    @ MultivariateOrthogonalPolynomials ~/.julia/dev/MultivariateOrthogonalPolynomials/src/triangle.jl:823
  [2] \(a::ContinuumArrays.TransformFactorization{…}, b::QuasiArrays.BroadcastQuasiVector{…})
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/yY1lG/src/plans.jl:26
  [3] transform_ldiv_size(::Tuple{…}, A::QuasiArrays.SubQuasiArray{…}, B::QuasiArrays.BroadcastQuasiVector{…})
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/yY1lG/src/bases/bases.jl:341
  [4] transform_ldiv(A::QuasiArrays.SubQuasiArray{…}, B::QuasiArrays.BroadcastQuasiVector{…})
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/yY1lG/src/bases/bases.jl:342
  [5] adaptivetransform_ldiv(A::JacobiTriangle{Float64, Int64}, f::QuasiArrays.BroadcastQuasiVector{Any, Base.Splat{var"#55#56"}, Tuple{Inclusion{…}}})
    @ ClassicalOrthogonalPolynomials ~/.julia/packages/ClassicalOrthogonalPolynomials/cd3TK/src/adaptivetransform.jl:37
  [6] transform_ldiv_size(::Tuple{…}, A::JacobiTriangle{…}, f::QuasiArrays.BroadcastQuasiVector{…})
    @ ClassicalOrthogonalPolynomials ~/.julia/packages/ClassicalOrthogonalPolynomials/cd3TK/src/adaptivetransform.jl:2
  [7] transform_ldiv(A::JacobiTriangle{Float64, Int64}, B::QuasiArrays.BroadcastQuasiVector{Any, Base.Splat{var"#55#56"}, Tuple{Inclusion{…}}})
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/yY1lG/src/bases/bases.jl:342
  [8] basis_ldiv_size
    @ ~/.julia/packages/ContinuumArrays/yY1lG/src/bases/bases.jl:393 [inlined]
  [9] copy
    @ ~/.julia/packages/ContinuumArrays/yY1lG/src/bases/bases.jl:385 [inlined]
 [10] materialize(M::ArrayLayouts.Ldiv{…}; kwds::@Kwargs{})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/RCHQ8/src/ldiv.jl:22
 [11] materialize
    @ ~/.julia/packages/ArrayLayouts/RCHQ8/src/ldiv.jl:22 [inlined]
 [12] ldiv
    @ ~/.julia/packages/ArrayLayouts/RCHQ8/src/ldiv.jl:98 [inlined]
 [13] \
    @ ~/.julia/packages/QuasiArrays/epj8s/src/matmul.jl:34 [inlined]
 [14] transform(A::JacobiTriangle{Float64, Int64}, f::Function)
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/yY1lG/src/bases/bases.jl:354
 [15] expand(A::JacobiTriangle{Float64, Int64}, f::Function)
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/yY1lG/src/bases/bases.jl:365
 [16] top-level scope
    @ REPL[126]:1
Some type information was truncated. Use `show(err)` to see complete types.

[MikaelSlevinsky]

Same error when replacing the quasi array with a function

julia> p = (n,k,a,b,c,x,y) -> jacobip(n-k,2k+b+c+1,a,2x-1) * (1-x)^k * jacobip(k,c,b,2y/(1-x)-1)
#59 (generic function with 1 method)

julia> f = expand(P, splat((x,y) -> p(n, k, a, b, c, x, y) ))
ERROR: MethodError: no method matching *(::FastTransforms.FTPlan{Float64, 2, FastTransforms.TRIANGLEANALYSIS}, ::Matrix{Any})

[MikaelSlevinsky]

Oh this is related to #140. In short, type instability of global variables causes the appearance of Any coefficients

[dlfivefifty]

The quick workaround is to declare the type:

julia> f = expand(P, splat((x,y) -> Pabc[SVector(x, y), Block(n+1)[k+1]] :: Float64))
JacobiTriangle(0, 0, 0) * [0.01884925224170753, 0.028935187711955863, 0.10672035718183336, -0.009142604629232859, 0.2958315403308107, 0.049882481929621275, -0.05954645634502389, 0.336771695932404, 0.1677915614316821, 0.14376137194923191  …  ]

Note in the new comprehension support in QuasiArrays.jl we choose the type based on evaluation on a random point:

https://github.com/JuliaApproximation/QuasiArrays.jl/blob/b92610c7262939f4f81f87d751572143b0fccb70/src/quasireducedim.jl#L284

This idea could be used more.