I've added a module that can approximate a Sersic profile (either de Vaucouleurs or exponential) as a mixture of Gaussians (MoG) using Hogg and Lang 2013 and analytically compute the second moments. This will be used heavily in the SprinkledRealizer (lenses in DC2) and CosmoDC2Realizer (non-lenses in DC2) subclasses. Because I'll most likely be using the module in a vectorized way, e.g. across all the lenses at a time, I tried to be flexible in the type of input arguments the module can take so it can compute second moments on both arrays and scalars.
I confused myself for a whole day trying to figure out the right normalization! So I've also added a tutorial with some sanity checks on the Sersic-as-MoG approximation to supplement the documentation.
It turns out the approximations aren't very close when the profiles are rendered on an image grid. I wasn't able to do a very rigorous "goodness of fit" test between the analytical Sersic profile and the MoG. Maybe I can do a chi-square test on the 1D radial slice of the MoG. The "badness" evaluation in Hogg and Lang 2013 seems to do something like that, but I couldn't understand the way they scaled the 1D chi-squared value.
Vectorization will be good for speed as long as we have enough memory! We'll have up to ~1B non-lenses and ~100K(?) lenses so, at some point, parallelizing using Spark might become necessary. Whereas I've only implemented the vectorized version here, I plan to write a version with two explicit for-loops (over the Sersic profiles to approximate and over the Gaussian components of each Sersic profile).
TL;DR This PR contains:
A new utility module that computes the analytical second moments of a Sersic profile by approximating it as a MoG, using vectorization
A draft of a tutorial that explains the approximation and does some simple checks to ensure that the original Sersic and the MoG are close
I've added a module that can approximate a Sersic profile (either de Vaucouleurs or exponential) as a mixture of Gaussians (MoG) using Hogg and Lang 2013 and analytically compute the second moments. This will be used heavily in the
SprinkledRealizer(lenses in DC2) andCosmoDC2Realizer(non-lenses in DC2) subclasses. Because I'll most likely be using the module in a vectorized way, e.g. across all the lenses at a time, I tried to be flexible in the type of input arguments the module can take so it can compute second moments on both arrays and scalars. I confused myself for a whole day trying to figure out the right normalization! So I've also added a tutorial with some sanity checks on the Sersic-as-MoG approximation to supplement the documentation. It turns out the approximations aren't very close when the profiles are rendered on an image grid. I wasn't able to do a very rigorous "goodness of fit" test between the analytical Sersic profile and the MoG. Maybe I can do a chi-square test on the 1D radial slice of the MoG. The "badness" evaluation in Hogg and Lang 2013 seems to do something like that, but I couldn't understand the way they scaled the 1D chi-squared value. Vectorization will be good for speed as long as we have enough memory! We'll have up to ~1B non-lenses and ~100K(?) lenses so, at some point, parallelizing using Spark might become necessary. Whereas I've only implemented the vectorized version here, I plan to write a version with two explicit for-loops (over the Sersic profiles to approximate and over the Gaussian components of each Sersic profile). TL;DR This PR contains: