SpectralKit.jl is not really intended to work with coefficients of the polynomials directly (it uses Chebyshev bases). It is possible but you have to iterate. Its intended uses case is getting a linear combination with a given coefficients, allocation-free. Eg

```
using SpectralKit, StaticArrays
dom = BoundedLinear(0, 200)
basis = smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(4, 2), Val(2)) ∘
coordinate_transformations(dom, dom)
θ = randn(dimension(basis)); # random coefficients
linear_combination(basis, θ, (1, 1))
```

The last line should not allocate when inside a function. You can iterate the last argument on `CartesianIndices`

, as in (continuing the example)

```
julia> using BenchmarkTools
julia> cids = CartesianIndices((200, 200))
CartesianIndices((200, 200))
julia> buffer = Vector{Float64}(undef, length(cids));
julia> function eval_at_θ!(buffer, basis, θ, cids)
for (i, x) in enumerate(cids)
buffer[i] = linear_combination(basis, θ, Tuple(x))
end
buffer
end
eval_at_θ! (generic function with 2 methods)
julia> @ballocated eval_at_θ!($buffer, $basis, $θ, $cids)
0
```

Note that I am never collecting the CartesianIndices, since they iterate just fine online.

One caveat: using Smolyak basis as tensor basis (as above) is not ideal and is not what this package has been optimized for. If there is a demand, I could introduce a `TensorBasis`

, it is much simple as it does not have to jump through all the hoops of Smolyak indexing.