Evaluation of SH/Zernike/Triangle expansions/polynomials?

[dlfivefifty]

@MikaelSlevinsky Is there a routine for pointwise evaluation for these bases?

CC @tsgut

[MikaelSlevinsky]

There's 1D clenshaw! for evaluating an expansion in a function set that satisfies a three-term recurrence. Can supply one initial condition as phi0 (the other is phi_{-1} == 0). There's also horner!.

The API is subject to change because I think they need so-called transposed versions (transposed Clenshaw algorithm is for evaluating coefficients from function samples).

Example code is https://github.com/JuliaApproximation/SphericalHarmonics.jl/issues/3#issuecomment-620726510

[dlfivefifty]

I guess the answer is no? That is, we'd have to implement for each case the example code?

[MikaelSlevinsky]

What kind of evaluation do you need? A single point? A Cartesian-product grid? A list of points?

[dlfivefifty]

A single point, and a list of points, with support for in-place. (Would the best way to do a list of points to loop over single points? This seems memory ideal and good for parallelisation.)

[MikaelSlevinsky]

Note that Cartesian-product is already supported with specific grids.

I guess a list of points can benefit from multithreading (and planning recurrence coefficients) but probably not SIMD?

Aliasing for in-place methods only partially fills the point sets in more than one dimension unless the in-placing is a vector field.

[MikaelSlevinsky]

I thought you had a 2D clenshaw for triangle polynomials?

[dlfivefifty]

I do, was hoping you had a better one 🤣

[MikaelSlevinsky]

What we really need is a feature-complete C library for nonuniform FFTs (and DCTs and DSTs). Two options include https://github.com/flatironinstitute/finufft and https://github.com/danfortunato/nufftw. It seems that neither one completely replicates the FFTW interface at the moment.

My 1D clenshaw code was really a case study in code generation for SIMD kernels so that I don't write them all by hand. If your points aren't Cartesian-product after some transformation, then I don't think I can do what you're asking better in C than in Julia, perhaps slapping on a Threads.@threads on a loop over the points.

[dlfivefifty]

NFFT seems massive overkill when the grid is just a point

[MikaelSlevinsky]

Exactly, but how important is performance for only one point?

[dlfivefifty]

Touche, but a very fast multithreaded 1-pt method would likely outperform NFFT for moderately large use cases

[dlfivefifty]

I'm going to start cleaning up clenshaw from ApproxFun and put it in FastTransforms.jl.

Do you have any idea why I pre-allocated the temporary vectors in ClenshawPlan?

https://github.com/JuliaApproximation/ApproxFunBase.jl/blob/master/src/LinearAlgebra/clenshaw.jl

This seems completely nuts to me...why would we ever want these stored in memory at all?

[dlfivefifty]

You should document your code...how do I know what clenshaw!(c::Vector{Float32}, x::Vector{Float32}, f::Vector{Float32}) does??

[MikaelSlevinsky]

At least there's tests :evergreen_tree:.

That's first kind Chebyshev summation. You can alias x with f if the points were created just to be overwritten.

There's also the OP-like clenshaw which has the method clenshaw!(c, A, B, C, x, phi0, f). A, B, and C are the DLMF recurrences and phi0 is the initial condition which it could be like phi0(x_j) = sin(x_j). Here, f can be aliased with either x or phi0 but not both (excepting the special case where phi0(x_j) == x_j).

[dlfivefifty]

Ok will add docs in a PR