Spergelets

[jmeyers314]

The next step in implementing the series representation of the Spergel profile (see issue #616 and PR #625) is the implementation of the series basis functions (which for now I'm deeming "Spergelets") and some framework for picking the appropriate coefficients for each Spergelet. I've already more or less got the Spergelets drawable on a local branch, but haven't started the rest of the framework yet. Here's my proposal though:

I think it would be nice if the implementation of the Spergel profile via Spergelets shares most of the same API as the directly implemented Spergel profile. For example:

psf = galsim.Moffat(beta=3, fwhm=0.67)

# Direct draw
gal1 = Spergel(nu=0.0, half_light_radius=0.5).shear(g1=0.2, g2=-0.3)
conv1 = galsim.Convolve(gal1, psf)
img1 = conv1.drawImage(nx=32, ny=32, scale=0.2)

# Series draw initiated by `jmax` keyword indicating order of the series expansion
# and `r1` indicating the fiducial radius about which to expand.
gal2 = Spergel(nu=0.0, half_light_radius=0.5, jmax=3, r1=0.4).shear(g1=0.2, g2=-0.3)
conv2 = galsim.Convolve(gal2, psf)
img2 = conv2.drawImage(nx=32, ny=32, scale=0.2)

My current plan to accomplish this, which is pretty close to the plan @mdschneider laid out in #616, is to create a GSObject subclass Series from which SpergelSeries would derive. SpergelSeries would then have methods: getCoeff((radial_index, azimuthal_index)), and getBasisFunc((radial_index, azimuthal_index)), and also override transformation methods like shear, dilate, etc. to implement these by changing the coefficients of the expansion. Similar to ChromaticConvolution, I would also imagine a SeriesConvolution subclass of Series with a drawImage method that does something like:

def SeriesConvolution.drawImage():
    #1) figure out nx, ny, scale
    #2) if not already precomputed, compute outer product of convolutions (at right 
         size/scale) of all series' terms and store for potential reuse (this uses 
         getBasisFunc).  For Spergelets, it's probably also useful to index these 
         convolutions by the fiducial profile size `r1` mentioned above.
    #3) determine outer product of all series' coefficients (uses getCoeff())
    #4) return appropriate linear combination of coefficients and basis functions.

There are at least two potential drawbacks that I can see:

• The Spergel (2010) expansion does not expand in centroid position. This might not be a problem if the use case is to sample from the likelihood function since under reasonable assumptions the centroid conditional probability can be computed from the autocorrelation of the model image and the observed image as in Miller++13 (LensFit).
• I can't see a way to avoid specifying an image size and scale. One could create an InterpolatedImage and redraw, I suppose, but that almost certainly destroys the speedup acquired from precomputing convolutions and making everything linear.

[mdschneider]

This looks like a very good plan.

Here are couple of further questions:

• Storing the pre-computed basis functions and their convolutions with the PSF might get complicated as the order of the series, jmax, increases. Does this merit consideration of the dictionary or other data structure for storing all this information? Is it possible that looking up saved basis functions could contribute significantly to the computation time of SeriesConvolution.drawImage()?
• To clarify, I think @jmeyers314 is proposing to assert a fixed centroid in the Spergelets implementation with the understanding that sampling of galaxy model parameters will proceed in two distinct steps: 1) sample from the Spergelet model likelihood, 2) sample galaxy centroid with all other galaxy model parameters fixed using the conditional probability determined from the auto-correlation of the data and model as in Miller et al. (2007) section 3.1.

[rmandelb]

Just a quick comment on the user interface.

Perhaps I am being overly influenced by our recent discussion about the ChromaticObject interpolation stuff, but an alternative to the user interface that you proposed would be something like

# Direct draw
gal1 = Spergel(nu=0.0, half_light_radius=0.5).shear(g1=0.2, g2=-0.3)
conv1 = galsim.Convolve(gal1, psf)
img1 = conv1.drawImage(nx=32, ny=32, scale=0.2)

# Series draw 
gal1.setupSeries(jmax=3)
conv2 = galsim.Convolve(gal1, psf)
img2 = conv2.drawImage(nx=32, ny=32, scale=0.2)

So the basic idea is that it's the same object, but you're changing how you do the internal calculations and image rendering to use the series approximation.