DESC icon indicating copy to clipboard operation
DESC copied to clipboard

Make Dommaschk factorials more accurate

Open dpanici opened this issue 1 year ago • 5 comments

n! / m!

n! = gamma(n+1) gammaln(n+1) = log(n!) exp(gammaln(n+1) - gammaln(m+1)) = n!/m!

dpanici avatar Sep 27 '24 15:09 dpanici

might fix some numerical issues I suspect are happening at higher resolutions in the dommaschk

dpanici avatar Sep 27 '24 15:09 dpanici

I don't know about your numerical issue but I remember that using factorials is much more accurate than gammaln, since latter uses some approximation instead of actually calculating the factorials. For non-integral numbers, that is the only choice but if n and m are integers, I would recommend using factorials. This was used in one of the Zernike PRs and improved the accuracy substantially, probably #1037

YigitElma avatar Sep 30 '24 00:09 YigitElma

yea but they need to be differentiable unfortunately. However... since they are a basis, I think I could make some sort of transform matrix for them and pre-compute the basis... I should look into that instead of doing this (as I just checked and I already use gammaln actually)

dpanici avatar Sep 30 '24 04:09 dpanici

Actually I cannot do that since R,Z,phi are the coordiantes and so we never have a fixed grid in them. I can, however, pre-compute all the factorials and use math.factorial for that, since those are purely functions of the mode numbers of the basis... yea I can do that at initialiazaiton of the Dommaschk field or something.

make a self._constants or something that will get passed into the dommaschk_potential call that contains the various non-coordinate-dependent and non-coefficient-dependent parts of the potential (might though be a bit more complex, as they will need to likely be np arrays of the same length as the number of modes in the basis, since each different mode has a different arg to these constant parts)

dpanici avatar Sep 30 '24 04:09 dpanici

possibly rising factorial

dpanici avatar Sep 30 '24 18:09 dpanici