Implement Spergel (2010) galaxy profile model

[mdschneider]

Motivation: We are interested in speeding up the rendering of galaxy images for use in MCMC algorithms to infer galaxy model parameters given a postage stamp image. An alternative galaxy model parameterization may be helpful.

The slowest part of the prediction of a galaxy image is the convolution with the PSF. Here's a plan for moving the PSF convolution into a pre-computation step, leaving MCMC steps to a simple linear operation for galaxy image predictions,

• Spergel (2010) demonstrates a Sersic-like profile, but with the important feature that it can be Taylor expanded in ellipticity and relative size in such a way that all model parameters are in coefficients of the expansion. Then, the PSF convolution can be pre-computed with each term in the expansion and subsequent model updates are a simple linear operation.
• If the PSF is represented by a Zernike expansion of the exit pupil wavefront, then it’s possible to derive an expansion of the PSF in terms of the Zernike coefficients. Then we can treat PSF models in pre-computation analogous to the Spergel galaxy profile model.
• A dependency graph for these features is here.
• A first step is to implement the functional form from Spergel (2010) as the new SBSpergel class, mimicking SBSersic.

While GalSim predictions can be sped up by changing GSParams settings, we would like to avoid identifying the combination of accuracy and speed tradeoffs for each application if a re-parameterization can be made to work.

Credit for benchmarks and planning: @jmeyers314

[rmandelb]

What kind of interplay between C++ and python layers do you anticipate needing for the second part of this (the PSF Zernike expansion)? Right now, all Zernike-related calculations happen in the python layer, so it seems like what you're describing could require substantial refactoring of the code in order to be feasible?

[mdschneider]

The chief requirement I am assuming above is that the PSF can be expanded in a series in such a way that any PSF model parameters are confined to the linear coefficients of the series. Then the terms with spatial dependence in the series can be individually convolved with the terms in the series expansion of the galaxy model.

By assuming a Zernike expansion of the exit pupil wavefront, it is possible to explicitly derive the series expansion of the PSF that has as coefficients the wavefront Zernike coefficients. See, for example, equations (7) and (8) here.

This approach will require new GSObjects (and corresponding SBProfiles) representing each term in the series expansion of the PSF. These objects are intended to be passed to the Convolve method along with the terms in the series expansion of the galaxy profile. So, I don't think we will need significant code refactoring. However, we may not directly utilize the existing Zernike wavefront methods (but these should be very important for validation).

This approach might work for wavefronts that are not well described by a low-order Zernike expansion. But I have only understood how a finite series might be useful for the case of a low-order Zernike expansion of the wavefront. So, crucially, this is not an approach for pre-computing atmosphere PSF convolutions.

[mdschneider]

Branch started for this issue: c2b9afa4a30c6ed382e2aec9e9eb9265adf820f7

[rmandelb]

Michael - I didn't get to fully read through and grok your message, but since you're already coding, I wanted to mention something you might find useful. If you follow the instructions in devel/commit-msg, then git will automatically append the branch name to each commit message for you, which means that the commits will show up on this issue page automatically (after you push them). That can be really useful for later discussions of the code that might happen on this issue.

[jmeyers314]

@mdschneider, I added a unit test for generating a few circularly symmetric Spergel profiles of varying indices. They await your completion of a "half_light_radius" input keyword.

[jmeyers314]

I'm guessing we can ignore the info cacheing here since the Spergel profile is analytic both in real space and Fourier space: no need to cache the results of Hankel transforms and such.

[rmjarvis]

You might still want the Info class even without the need for caching the Hankel transform. cf. SBExponential, which mostly uses its Info class to cache the calculations needed for photon shooting.

But if that is analytic too, then you can skip all the Info stuff here.

[jmeyers314]

So I think I did more harm than good in b0eadfe. Looking for some git-fu: anyone know how to excise just the one commit from this branch?

-Josh

[rmjarvis]

The simple revert function: git revert b0eadfe should be fine in your case. BTW, when things get more complicated, this is a useful site to check out:

http://sethrobertson.github.io/GitFixUm/fixup.html

But I believe your "adventure" ends with git revert.

[jmeyers314]

I think I got Fourier mode drawing basically working for this now. GalSim output matches my straight python test script. A couple of interesting features popped up:

• For nu < 0, the surface brightness diverges at the origin, so I had to hack the 'sb' normalization unit test for this position.
• related to this, maxK gets really crazy huge for nu close to -1. For half_light_radius=1, and nu=-0.9, I compute maxK=3e13! Even Sersic(n=6.2) only gets up to ~600. I sanity checked nu=0.5, which is an exponential and I get the same maxK as galsim.Exponential, so I think I'm doing this right. Hopefully not important once convolutions are included.

Next step, photon-shooting.

[rmandelb]

It is certainly the case that we've run into issues when people try to draw deVauc profiles without a convolution, and our standard solution is to suggest convolving with at least a small PSF. So I'm guessing that fix should work for you. It sounds like something to mention in the docstring.

[jmeyers314]

I think the first part of this issue (drawing the Spergel profile without any Taylor expansion trickery) is almost ready for a PR. The last remaining question I have is what to do for photon-shooting when nu<0 (in which case the s.b. profile diverges at the origin). I think this should be possible in principle, but in practice the rejection sampling near the origin ends up in an infinite loop waiting for a deviate smaller than 1/inf. Any suggestions?

Some other items of note:

• Drawing a nu=-0.6 Spergel profile is about 3x faster than drawing a n=4 Sersic profile, which has about the same shape near the half_light_radius. A nu=0.5 Spergel takes the same amount of time as an n=1 Sersic (both exponentials).
• I think I found an error in the Spergel(2010) Taylor expansion which when corrected increases the number of terms in the expansion rather significantly, so I'm a bit less optimistic now about precomputing all the convolutions.
• I noticed a couple weeks ago that the Spergel profile is the F.T. of the Moffat profile. Not sure how to take advantage of that but I thought it was interesting. :)

[rmjarvis]

in practice the rejection sampling near the origin ends up in an infinite loop waiting for a deviate smaller than 1/inf. Any suggestions?

Where exactly is the problem happening? Making the OneDimensionalDeviate? Or when drawing from the interval (Interval::drawWithin)?

Anyway, some things you could try:

• Don't go all the way to 0 for the range you provide to the OneDimensionalDeviate constructor. i.e. range[0] = epsilon. You can figure out what epsilon is appropriate based on shoot_accuracy, which dictates how much flux you are allowed to miss when shooting.
• Similarly, you can have the radial function that you give to OneDimensionalDeviate flatten off at the origin in a way that misses less than shoot_accuracy worth of flux. In this case, you would want to have range be [ 0, epsilon, shoot_maxr ], since you'll want to let it know that the profile changes non-smoothly at whatever epsilon where you flatten it out.
• Or rather than miss any flux, you can have the central region at r < epsilon have a higher surface brightness such that the total flux there matches what the real profile would have. Just the s.b. would be discontinuous at r=epsilon (so definitely need to include that value in the range). The appropriate value for epsilon then would be where 2pi int_0^epsilon( r abs(f_true(r) - f_approx(r)) ) < shoot_accuracy. i.e. the integrated absolute error in the flux.
• You could try to get even more fancy by trying to come up with a function to use at r < epsilon that has the right flux and the right value at r=epsilon, but has f(0) < infinity (and presumably would have a discontinuous first derivative at r=epsilon). I'm not sure that this is worth it, but you could think about it.

I noticed a couple weeks ago that the Spergel profile is the F.T. of the Moffat profile. Not sure how to take advantage of that but I thought it was interesting. :)

Aside from using the analytic form in kValue, probably not much else to do. The Moffat class already uses this in reverse. We used to use a tabulated FT, before I discovered the analytic version -- it was a significant speed up to switch to using cyl_bessel_k.

[jmeyers314]

Where exactly is the problem happening? Making the OneDimensionalDeviate? Or when drawing from the interval (Interval::drawWithin)?

Problem was in Interval::drawWithin. I went ahead and implemented your fancy suggestion above, which turned out to not be too difficult. Seems to work fine now.

[rmandelb]

@jmeyers314 or @mdschneider , should this issue be closed? Or did you have a further set of work to be done on this issue after the initial PR?

[jmeyers314]

@jmeyers314 or @mdschneider , should this issue be closed? Or did you have a further set of work to be done on this issue after the initial PR?

I had originally envisioned keeping this open until #628 is merged, as the original description by @mdschneider included the series implementation of the Spergel and PSF profiles. Maybe it's sufficient to just keep 628 open though.

[rmandelb]

It's up to you, really, but it seems to me that having two open issues for basically the same body of work isn't necessary.

[jmeyers314]

I guess that's fair. I'm certainly not working in branch 616 anymore. Closing.