Chromatic Real Galaxy

[jmeyers314]

GalSim has both RealGalaxy objects, and ChromaticObjects, but doesn't have ChromaticRealGalaxy (yet!).

@rmjarvis, @rmandelb, and I talked a bit about how one might take high-res HST images of a galaxy through 2+ filters, deconvolve the known (chromatic) HST PSF, and then reconvolve by a different chromatic PSF (to presumably less resolution than HST, c.f. RealGalaxy). This would be especially helpful for studies of galaxies with color-gradients for LSST/Euclid/etc.

Semboloni++13 section 5.2 equation 25 briefly lays this out for the case where the SED at each point in the galaxy is linear in wavelength. We think we can generalize this to the case of generic SEDs, for instance, typical bulge and disk SEDs.

This issue is to develop the framework for a ChromaticRealGalaxy object. For the moment, I think it's reasonable to limit the scope to just the framework. In particular, this issue will not attempt to actually assemble a catalog of HST images and PSFs; we'll do that later, and in the meantime we can test with simulated GalSim images.

[jmeyers314]

Hey @rmjarvis, @rmandelb, @gbernstein, @barnabytprowe.

I'm in the middle of figuring out how to propagate correlated noise through RealChromaticGalaxy, which is under construction in this branch. See some notes on this here. I believe the right way to handle the noise propagation is to work Fourier mode by Fourier mode. However, I'm finding that I'm unsure about the right way to consistently compute the required Fourier coefficients for 3 different objects: the input noise power spectrum (given a CorrelatedNoise instance), the input HST images, and the (SED*bandpass)-weighted input PSFs.

In particular, I see that the Fourier transforms in CorrelatedNoise are done via np.fft, rather than using the stored profile's .drawKImage() method. What's the reason for this? I created the script pasted below to investigate the difference between drawK and np.fft, and realized quickly that for the amplitudes at least, the difference is precisely the left panel of Bernstein&Gruen2014 Figure 1. (So why not use drawKImage with x_interpolant='sinc'?) I don't understand the phase differences, though I suppose these don't matter for the noise generation in CorrelatedNoise since new noise gets generated with uniformly random phase anyway.

Any thoughts about np.fft vs drawKImage(x_interpolant='sinc') (or other x_interpolant?) for RealChromaticGalaxy?

-Josh

Script:

import matplotlib.pyplot as plt
import numpy as np
import galsim

rng = galsim.BaseDeviate(1)
xi = galsim.getCOSMOSNoise(rng=rng)
scale = 0.03

Nx = 256
Nx2 = Nx/2

img = galsim.Image(Nx, Nx, scale=scale)
img.addNoise(xi)

fig = plt.figure()
ax = fig.add_subplot(111)

fig2 = plt.figure()
ax2 = fig2.add_subplot(111)

for interpolant in ['linear', 'cubic', 'quintic',
                    'lanczos3', 'lanczos4', 'lanczos5',
                    'sinc']:

    # Fourier transform via drawKImage
    stepk = 2.0*np.pi/(Nx*scale)
    maxk = np.pi/scale
    ii = galsim.InterpolatedImage(img, calculate_stepk=stepk, calculate_maxk=maxk,
                                  x_interpolant=interpolant)
    re, im = ii.drawKImage(nx=Nx, ny=Nx, scale=stepk)
    drawKout = re.array + 1j * im.array
    # The alternating sign flips on the following line seems to help match the phase (but not
    # perfectly), though I don't know why.
    drawKout[::2, ::2] *= -1
    drawKk = (np.arange(Nx)-Nx2)*stepk

    # Fourier transform via np.fft
    fftout = np.fft.fft2(img.array)
    fftk = np.fft.fftfreq(Nx, scale) * 2.0 * np.pi
    fftout = np.fft.fftshift(fftout)  # Center k=(0,0)
    fftk = np.fft.fftshift(fftk)
    assert np.allclose(drawKk, fftk)  # Make sure we've queried the same modes.
    ax.plot(fftk, 1.0-np.abs(drawKout[Nx2, :])/np.abs(fftout[Nx2, :]), label=interpolant)

    phases = np.angle(drawKout[Nx2, :]) - np.angle(fftout[Nx2, :])
    phases = phases % (2 * np.pi)
    ax2.plot(fftk, phases, label=interpolant)
ax.set_xlim(0, maxk)
ax.set_ylim(-0.01, 0.01)
ax.legend(loc='best')
ax.set_ylabel("1 - abs(drawK) / abs(fft)")
ax.set_xlabel("k")
fig.show()

ax2.legend(loc='best')
ax2.set_ylabel("arg(drawK) - arg(fft)")
ax2.set_xlabel("k")
fig2.show()
plt.show()

[rmjarvis]

I see that the Fourier transforms in CorrelatedNoise are done via np.fft, rather than using the stored profile's .drawKImage() method. What's the reason for this?

I think the answer to this is simply that this is what Barney was familiar with, and he wrote the code for CorrelatedNoise. I don't think anyone has investigated any of these differences, so it is certainly quite possible that drawKImage would be more accurate and/or more efficient.

See some notes on this here.

I didn't try to follow all the math here yet, but I think the broad brush strokes seem right to me. Also, an appropriate place for this would be in devel/modules. We have some other latex writeups of some of the math we use in the code there.

the input noise power spectrum (given a CorrelatedNoise instance),

The one thing that I was wondering though is whether we can assume that the noise in each of the input HST images is the same. Is that necessarily valid? I think that's an assumption of your formalism, but it seems like we might at some point have observations of a galaxy in different bands from different telescopes. Trying to think about how to compute the net noise profile of the final object with an arbitrary SED seems complicated though, so maybe we can ignore this for now...

[jmeyers314]

I don't think anyone has investigated any of these differences, so it is certainly quite possible that drawKImage would be more accurate and/or more efficient.

Good to know. I'll try some timing comparisons too in that case. I think the salient point is that the FT of the correlation function should correctly predict the variances of the image FTs. This can be checked via Monte Carlo.

an appropriate place for this would be in devel/modules

Good idea. Will migrate.

The one thing that I was wondering though is whether we can assume that the noise in each of the input HST images is the same.

I think you're right that this is a bad assumption, but I don't think I actually assumed it :) Note the subscript 'i's on the xi's and Pks.

The part that currently worries me the most is estimating the covariances of the weighted least squares solutions. I haven't convinced myself yet whether or not the formula given here is a complete description of the uncertainties in the fit or just a Gaussian approximation. If it's complete, then I think it will be possible to properly propagate the uncertainties.

[rmjarvis]

I think you're right that this is a bad assumption, but I don't think I actually assumed it :) Note the subscript 'i's on the xi's and Pks.

Ah. Good. I had missed that. Or misunderstood, really. Thanks.

[rmandelb]

In particular, I see that the Fourier transforms in CorrelatedNoise are done via np.fft, rather than using the stored profile's .drawKImage() method. What's the reason for this?

I'm almost positive that the code in CorrelatedNoise predates the existence of drawKImage.

[rmjarvis]

I'm almost positive that the code in CorrelatedNoise predates the existence of drawKImage.

Ha. Good point. Yes, you're right of course. Another good reason for Barney not to have used that.

[rmandelb]

And to address Josh's main point: I see no reason not to update CorrelatedNoise to use this newer functionality that wasn't available when the code was originally written. We do have a pretty good set of unit tests that should be able to pick up on anything that gets broken by doing this (Josh - if you try this I would recommend running test_correlatednoise.py from the command line, to get some more stringent tests).

I do wonder about speed issues though; maybe I'm off base but I usually think of sinc interpolation as being slow if you want to do it right.

[barnabytprowe]

Hi guys! Funnily enough, I couldn't myself remember if something called drawKImage did exist or not at the time I wrote the correlated noise stuff and was quite happy to go along with Mike's original explanation, it was highly plausible to me anyway :)

Just wanted to let you know that as this is material which I was pretty intimately involved in writing and testing (the latter involving tackling a few tricky, conceptual fourier space issues with the help of others on the team) I'm happy to help as much as I can in explaining anything in the old code that needs it - I'll try keep a closer eye on my gmail to help with this.

@rmandelb I think you're right - the sinc interpolation will be slow. Actually, at the heart of how correlated noise stuff is used in galsim is an approximate rendering of the noise correlation function using the blinear interpolant (see, e.g., lines around ~990 in correlated_noise.py... there's also a Github discussion of this very point somewhere, I'll try to find it...). This is a fudge, one that works well enough in practice so far, as evidenced by our tests, but that would fail if you wanted to make an image of HST correlated noise at super HST resolution, for example.

(The problem is that an image noise correlation function as a mathematical surface brightness profile is really a load of delta functions convolved with galsim.Pixel objects... The power spectrum is infinitely repeating. The DFT handles this all naturally, once you've got a reasonably accurate image of the correlation function to start from. It may be that drawKImage does too, but I'm not really knowledgeable about what it does!)

[barnabytprowe]

Phew, took a bit of finding, but here it is:

https://github.com/GalSim-developers/GalSim/pull/452#discussion-diff-5701561

Please excuse my being very slow to enlighten myself as to the issue in the discussion above!! It's almost painful to read :)

But if you work through the discussion I think you might get the point. Using drawKImage might provide a really elegant solution to this issue, or it might run aground against the same issue: the numerical difficulty of seeking to accurately describe a discrete, sharp-edged object (the correlation function of noise in an image) using our continuous, surface brightness profile machinery.

[jmeyers314]

Hey @barnabytprowe! Thanks for the link. I'll definitely take a long look.

In the meantime, I can report that I've found the right combination of np.fft.fftshift and np.fft.ifftshift to make the numpy and drawK phases match, so I'm a bit less concerned about mixing the two methods now. Hopefully I'll be able to report on some end-to-end noise propagation tests in the coming few days.

[jmeyers314]

There's some discussion of using drawK inside CorrelatedNoise here: https://github.com/GalSim-developers/GalSim/issues/430#issuecomment-22287517. Looks like at least one stumbling block at the time was what to do with non-square images, which imply the need for Fourier grids with non-square "pixels". I think this would potentially be fixable by handling more generic wcs's in drawKImage; something like JacobianWCS([dkx, 0, 0, dky]).

However, I think I'm now close to getting the noise right for ChromaticRealGalaxy without switching to drawKImage inside CorrelatedNoise. Here are some results for an end-to-end test, where I generate 1000 "Euclid" images given 2 * 1000 "HST" images of pure noise in 2 bands. I then measure the noise correlations in the output images and compare these to a prediction that only depends on the input noise models, input PSFs and output PSF.

First in 2D, the observed correlation function:

Next, the predicted correlation function (this one's deterministic, no pseudo-random numbers involved):

And a slice across the middle row of each for comparison:

This is a bit slow to run right now, ~90 min, so I'll have to think up something simpler for a unit test.

[jmeyers314]

So transferring the code from my nice pretty ipython notebook into GalSim has turned out to be a bit trickier than I thought it would be... Hopefully I'll figure that soonish though.

In the meantime, I've placed my notes into devel/modules/CGNotes.tex as suggested by @rmjarvis , and also added a section attempting to push the Bernstein&Gruen interpolation equations through the Wiener-Khinchin theorem to see if I could robustly derive the (co-)variance of interpolated Fourier amplitudes given a discretely sampled auto-covariance function and asserted real/Fourier interpolation kernels. The results of that (assuming I did the math right), is that even if the discrete version of image noise is stationary, and can hence be described by a discrete covariance function / power spectrum, a continuously interpolated noise image will possess non-stationary noise, or equivalently, that the interpolated Fourier amplitudes are correlated. In hindsight, I think this makes sense, since, for instance, I would expect the variance of interpolated samples to always be less than the variance of the samples themselves. I'm not sure if this helps or not with determining the "right" way to propagate correlated noise, however.

[jmeyers314]

Issue defeated. At last. ( #687 ).

[plazas]

Nice!