[WIP] Add quasi-discrete spherical Bessel transform

[sethaxen]

This PR is built on #2 and adds the spherical Hankel transform discussed in #1. I ended up implementing the hyperspherical case actually, since it was only marginally more difficult than the cylindrical and spherical cases.

I still need to write tests and document the transform. I also need to verify that the approximation at the bottom of https://chrisbrahms.github.io/Hankel.jl/dev/derivations/#Integration-of-functions-1 follows in the general case.

[sethaxen]

Okay, all of the QDSHT functionality should be implemented and tested, so this PR is ready for review.

I still need to document the QDSHT type and a few methods (but these will be nearly identical to those for QDHT), and something about QDSHT should go in the docs; how would you like that to be structured?

[chrisbrahms]

Nice work!

Since you've generalised this over dimensions as well, and just using n=1 results in the QDHT, it seems to me that having QDHT and QDSHT be separate entities makes little sense. So QDHT could just become the overarching type for everything, defaulting to p=0 and n=1, and perhaps an alias QDSHT could default to the spherical case. Or am I missing something? The only thing to check would be that performance is identical (I use QDHTs in hot loops) but I can't see why it wouldn't be.

re: docs, is there a reference that shows how to derive the SHT as a re-parameterisation of the HT? If not, ideally you'd take the time to add that to the derivations page. Since your new code encompasses pretty much everything from before, you could just edit the main docs page with the more general definitions and leave the QDHT as an important special case.

[sethaxen]

Nice work!

Thanks! Your code, tests, and derivation were very nice, so it made it relatively easy to extend.

Since you've generalised this over dimensions as well, and just using n=1 results in the QDHT, it seems to me that having QDHT and QDSHT be separate entities makes little sense. So QDHT could just become the overarching type for everything, defaulting to p=0 and n=1, and perhaps an alias QDSHT could default to the spherical case. Or am I missing something?

No, you're quite right. The two should be entirely equivalent, and this is also currently tested.

The only thing to check would be that performance is identical (I use QDHTs in hot loops) but I can't see why it wouldn't be.

As for performance, there might be a slight decrease in performance moving from QDHT to QDSHT with n=1 from control flow checking n and a few powers, but it will be small, and it will be entirely handled in set-up. The only other performance decrease could come if you call onaxis in a loop, since we don't have the J00 constant anymore and have to compute it ourselves. But I'll test to make sure. Adding n to the type signature for QDHT would handle this and make aliasing easier.

re: docs, is there a reference that shows how to derive the SHT as a re-parameterisation of the HT? If not, ideally you'd take the time to add that to the derivations page.

No, there's not. There is a reference that generalizes the QDHT to the half-integer order (it actually predates the Yu paper but IIRC didn't do as extensive error analysis): https://doi.org/10.1016/0010-4655(87)90204-9. And no papers that I've seen handled the hyperspherical bessel function case, so I can add those derivations.

Since your new code encompasses pretty much everything from before, you could just edit the main docs page with the more general definitions and leave the QDHT as an important special case.

Sounds good! I'll get started harmonizing the two. I'll probably along the way remove the AbstractQDHT introduced in #2, since I don't anticipate another QDHT plan.

[sethaxen]

Actually, since the proposed changes would interact with a lot of the changes in #2, I'll wait until that PR is finalized and merged.

[chrisbrahms]

re: #2 vs #3, doesn't the proposed unification of QDHT and QDHST make most of #2 unnecessary? Essentially, you'd just edit QDHT to handle arbitrary dimensions and orders instead of adding abstraction and splitting into 2 types. The reorganisation into qdht.jl and utils.jl still makes sense to me. The only thing we'd lose with this approach is the specificity of the docstrings. I'm open to counterarguments though!

As for performance, there might be a slight decrease in performance moving from QDHT to QDSHT with n=1 from control flow checking n and a few powers, but it will be small, and it will be entirely handled in set-up. The only other performance decrease could come if you call onaxis in a loop, since we don't have the J00 constant anymore and have to compute it ourselves. But I'll test to make sure. Adding n to the type signature for QDHT would handle this and make aliasing easier.

Slowing down the constructor a little bit isn't a problem, at least not for me. It's the actual transform that needs to be fast (hence the convoluted definitions of dot!). Having both p and n part of the type signature absolutely makes sense, and then we could even keep the J00 constant and use that in the special case of n=1, p=0.

No, there's not. There is a reference that generalizes the QDHT to the half-integer order (it actually predates the Yu paper but IIRC didn't do as extensive error analysis): https://doi.org/10.1016/0010-4655(87)90204-9. And no papers that I've seen handled the hyperspherical bessel function case, so I can add those derivations.

Yes, that would be great. Makes me think that somebody should publish that...?

[sethaxen]

re: #2 vs #3, doesn't the proposed unification of QDHT and QDHST make most of #2 unnecessary? Essentially, you'd just edit QDHT to handle arbitrary dimensions and orders instead of adding abstraction and splitting into 2 types. The reorganisation into qdht.jl and utils.jl still makes sense to me. The only thing we'd lose with this approach is the specificity of the docstrings. I'm open to counterarguments though!

Good point, yeah, I'll open a new PR that copies over the minimal necessary changes from each PR.

Perhaps you already did this, but I also stumbled upon a more formal way to explain the approximation in https://chrisbrahms.github.io/Hankel.jl/dev/derivations/#Integration-of-functions-1 and (probably) also quantify its error. I'll include it in the derivation.

No, there's not. There is a reference that generalizes the QDHT to the half-integer order (it actually predates the Yu paper but IIRC didn't do as extensive error analysis): https://doi.org/10.1016/0010-4655(87)90204-9. And no papers that I've seen handled the hyperspherical bessel function case, so I can add those derivations.

Yes, that would be great. Makes me think that somebody should publish that...?

It'd be a pretty simple little paper and (probably?) not too hard to write, but I'll start with updating derivation.md and then see if I (we) want to polish it into a paper.

[chrisbrahms]

Good point, yeah, I'll open a new PR that copies over the minimal necessary changes from each PR.

Good plan!

Perhaps you already did this, but I also stumbled upon a more formal way to explain the approximation in https://chrisbrahms.github.io/Hankel.jl/dev/derivations/#Integration-of-functions-1 and (probably) also quantify its error. I'll include it in the derivation.

I basically found out that this works by accident, verified it numerically and did a little digging but couldn't find a more formal explanation. Adding that would be very useful!

[sethaxen]

I basically found out that this works by accident, verified it numerically and did a little digging but couldn't find a more formal explanation. Adding that would be very useful!

Yeah, I stumbled upon it by accident too. It turns out that the inner sum is a truncated expansion of the Fourier-Bessel series of the function f(r) = 1/2. Hence the absolute error of using 1/2 is just the truncation error.