Piracy breaks `expand(Fourier(), cos)`

[DanielVandH]

Loading HarmonicOrthogonalPolynomials breaks expand(Fourier(), f) (and probably other methods, but this is the only issue I keep hitting).

julia> using ClassicalOrthogonalPolynomials

julia> expand(Fourier(), cos)
Fourier{Float64}() * [-5.730183352904034e-17, 1.2877173848170242e-16, 0.9999999999999999, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  ]

julia> using HarmonicOrthogonalPolynomials

julia> expand(Fourier(), cos)
ERROR: MethodError: no method matching grid_axis(::BlockArrays.BlockedOneTo{Int64, InfiniteArrays.InfStepRange{…}}, ::Fourier{Float64}, ::Block{1, Int64})
The function `grid_axis` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  grid_axis(::Base.OneTo, ::Any, ::Block)
   @ ContinuumArrays C:\Users\danjv\.julia\packages\ContinuumArrays\Py5Q6\src\bases\bases.jl:185

Stacktrace:
  [1] grid_layout(::ContinuumArrays.BasisLayout, P::Fourier{Float64}, n::Block{1, Int64})
    @ ContinuumArrays C:\Users\danjv\.julia\packages\ContinuumArrays\Py5Q6\src\bases\bases.jl:183
  [2] grid(P::Fourier{Float64}, n::Block{1, Int64})
    @ ContinuumArrays C:\Users\danjv\.julia\packages\ContinuumArrays\Py5Q6\src\bases\bases.jl:177
  [3] plan_grid_transform(P::Fourier{Float64}, szs::Tuple{Block{1, Int64}}, dims::Int64)
    @ ContinuumArrays C:\Users\danjv\.julia\packages\ContinuumArrays\Py5Q6\src\bases\bases.jl:213
  [4] _sub_factorize(::Tuple{…}, ::Tuple{…}, ::QuasiArrays.SubQuasiArray{…}; kws::@Kwargs{})
    @ HarmonicOrthogonalPolynomials C:\Users\danjv\.julia\packages\HarmonicOrthogonalPolynomials\yHaVS\src\multivariateops.jl:94
  [5] _factorize(::ContinuumArrays.SubBasisLayout, ::QuasiArrays.SubQuasiArray{…}; kws::@Kwargs{})
    @ ContinuumArrays C:\Users\danjv\.julia\packages\ContinuumArrays\Py5Q6\src\bases\bases.jl:248
  [6] factorize(::QuasiArrays.SubQuasiArray{Float64, 2, Fourier{…}, Tuple{…}, false}; kws::@Kwargs{})
    @ QuasiArrays C:\Users\danjv\.julia\packages\QuasiArrays\31vF4\src\inv.jl:57
  [7] plan_ldiv(A::QuasiArrays.SubQuasiArray{…}, B::QuasiArrays.BroadcastQuasiVector{…})
    @ ContinuumArrays C:\Users\danjv\.julia\packages\ContinuumArrays\Py5Q6\src\bases\bases.jl:279
  [8] transform_ldiv_size(::Tuple{…}, A::QuasiArrays.SubQuasiArray{…}, B::QuasiArrays.BroadcastQuasiVector{…})
    @ ContinuumArrays C:\Users\danjv\.julia\packages\ContinuumArrays\Py5Q6\src\bases\bases.jl:282
  [9] transform_ldiv(A::QuasiArrays.SubQuasiArray{…}, B::QuasiArrays.BroadcastQuasiVector{…})
    @ ContinuumArrays C:\Users\danjv\.julia\packages\ContinuumArrays\Py5Q6\src\bases\bases.jl:283
 [10] adaptivetransform_ldiv(A::Fourier{Float64}, f::QuasiArrays.BroadcastQuasiVector{Float64, typeof(cos), Tuple{…}})
    @ ClassicalOrthogonalPolynomials C:\Users\danjv\.julia\packages\ClassicalOrthogonalPolynomials\CRJvN\src\adaptivetransform.jl:34
 [11] transform_ldiv_size(::Tuple{…}, A::Fourier{…}, f::QuasiArrays.BroadcastQuasiVector{…})
    @ ClassicalOrthogonalPolynomials C:\Users\danjv\.julia\packages\ClassicalOrthogonalPolynomials\CRJvN\src\adaptivetransform.jl:2
 [12] transform_ldiv(A::Fourier{Float64}, B::QuasiArrays.BroadcastQuasiVector{Float64, typeof(cos), Tuple{Inclusion{…}}})
    @ ContinuumArrays C:\Users\danjv\.julia\packages\ContinuumArrays\Py5Q6\src\bases\bases.jl:283
 [13] basis_ldiv_size
    @ C:\Users\danjv\.julia\packages\ContinuumArrays\Py5Q6\src\bases\bases.jl:332 [inlined]
 [14] copy
    @ C:\Users\danjv\.julia\packages\ContinuumArrays\Py5Q6\src\bases\bases.jl:324 [inlined]
 [15] materialize(M::ArrayLayouts.Ldiv{…}; kwds::@Kwargs{})
    @ ArrayLayouts C:\Users\danjv\.julia\packages\ArrayLayouts\7omvS\src\ldiv.jl:22
 [16] materialize
    @ C:\Users\danjv\.julia\packages\ArrayLayouts\7omvS\src\ldiv.jl:22 [inlined]
 [17] ldiv
    @ C:\Users\danjv\.julia\packages\ArrayLayouts\7omvS\src\ldiv.jl:98 [inlined]
 [18] \
    @ C:\Users\danjv\.julia\packages\QuasiArrays\31vF4\src\matmul.jl:34 [inlined]
 [19] transform(A::Fourier{Float64}, f::Function)
    @ ContinuumArrays C:\Users\danjv\.julia\packages\ContinuumArrays\Py5Q6\src\bases\bases.jl:295
 [20] expand(A::Fourier{Float64}, f::Function)
    @ ContinuumArrays C:\Users\danjv\.julia\packages\ContinuumArrays\Py5Q6\src\bases\bases.jl:306
 [21] top-level scope
    @ REPL[6]:1
Some type information was truncated. Use `show(err)` to see complete types.

I think the issue is the piracy in this line

https://github.com/JuliaApproximation/HarmonicOrthogonalPolynomials.jl/blob/00e3a5f91754f5a1a8ba2a976e7bbcbede5f6a0f/src/multivariateops.jl#L94

which would make sense why I'm only hitting it with Fourier() currently since it uses blocks (unlike e.g. Legendre()).

[dlfivefifty]

We should probably restructure support for block arrays as an extension in ContinuumArrays.jl

[DanielVandH]

Another temporary fix in the meantime on the ClassicalOrthogonalPolynomials.jl side could be to have a method like

function ContinuumArrays.plan_grid_transform(F::Fourier, szs::Tuple{Block{1,Int}}, dims=ntuple(identity, Val(1)))
    return ContinuumArrays.plan_grid_transform(F, last(axes(F, 2)[only(szs)]), dims)
end

to force the plan_grid_transform necessary with the pirated method to work (with the last(...) needed to convert the szs to an Int since the transform seems to work in that case).