PSF Response experiment

[andrevitorelli]

I'm opening this issue to document the investigation of using AutoMetaCal to calculate PSF responses.

The notebook can be found at PSF_experiment_notebook.ipynb.

The data is 100k galaxies with random shape noise (~.3), and a small amount of pixel noise. A constant shear is applied to all galaxies before convolution with a constant (slightly elliptical) psf and noise addition.

The PSF response code does:

1. Deconvolve the observed image by the input PSF model
2. Apply some shear g to the deconvolved image, and g_psf to the PSF model
3. reconvolve the sheared galaxy with the sheared PSF
4. get the measured ellipticity in the metacal image e_measured gradients of the measured ellipticity respective to the applied shears R, & R_psf
5. measure the ellipticity of the PSF e_psf

From the individual ellipticity and response estimates we calculate the shear as:

e = <R^-1> - (1)

The idea of the test is to compare this to not using the PSF response and reconvolving with a symmetrised, enlarged PSF.

For now, step 3. does not enlargen the PSF as ngmix does, but this can be easily implemented as there is code for this.

As a first test, I'm comparing this to ngmix. The first issue I encounter is that when using round input PSFs, I get unbiased final results ( residual m1 is compatible with zero). But when shearing the PSF, the procedure above (1) doesn't work.

[EiffL]

Thanks for opening this issue :-)

Let me start with a couple of questions:

• Have you checked whether in Huff&Mandelbaum they were correcting for the PSF in their shape measurement? (independently of metacalibrating the PSF ellipticity).

As a first test, I'm comparing this to ngmix

• Hummm how are you comparing with ngmix? are you using a version with PSF metacalibration?

• hummm is the multiplicative bias really the one we want to look at here?

[aguinot]

Apply some shear g to the deconvolved image, and g_psf to the PSF model

Juste to make sure, you the shear to the galaxy and PSF separately ? Like with a classic method you end up with 9 images (e.g. 1 noshear + 4 gal_shear (classic metacal) + 4 psf_shear (for psf response)).

[aguinot]

@EiffL a beginning of response.

• Have you checked whether in Huff&Mandelbaum they were correcting for the PSF in their shape measurement? (independently of metacalibrating the PSF ellipticity).

In Huff&Mandelbaum they used different shape measurement method including Re-Gauss (so with PSF correction). But I don't remember if they compute the PSF response.

As a first test, I'm comparing this to ngmix

• Hummm how are you comparing with ngmix? are you using a version with PSF metacalibration?

👍

• hummm is the multiplicative bias really the one we want to look at here?

That is a very good question. I didn't thought to much about it. Maybe we are more interested in the additive bias. But the PSF also contribute to the multiplicative bias but it is extremely small. Another thing is that the PSF leakage come mainly from a bad modeling of the PSF (see DES Y1 or DES SV) and the intern handling by metacal is "perfect" so I would expect to have a PSF response = 0. If we want to test I think we need to alter the input PSF. Meaning that we simulate the galaxy with a PSF, P, and we give a PSF, P+dp, to autometacal. But I have no idea how to relate "dp" to the PSF response.

[andrevitorelli]

Thanks for opening this issue :-)

Let me start with a couple of questions:

* Have you checked whether in Huff&Mandelbaum they were correcting for the PSF in their shape measurement? (independently of metacalibrating the PSF ellipticity).

Yes, I have and it is still not completely clear to me. It's different, though, because they reconvolve with a symmetric, enlarged PSF, while we're trying to reconvolve with the same PSF.

As a first test, I'm comparing this to ngmix

* Hummm how are you comparing with ngmix? are you using a version with PSF metacalibration?

Comparison with ngmix is a placeholder at the moment - I just do the standard ngmix test without any additional PSF correction. Additionally, it provides an initial validation by starting from zero PSF ellipticities.

* hummm is the multiplicative bias really the one we want to look at here?

Not only but it is a problem that we get biased m when using anisotropic psfs and trying to correct for it.

[andrevitorelli]

Apply some shear g to the deconvolved image, and g_psf to the PSF model

Juste to make sure, you the shear to the galaxy and PSF separately ? Like with a classic method you end up with 9 images (e.g. 1 noshear + 4 gal_shear (classic metacal) + 4 psf_shear (for psf response)).

Yes, exactly. This is one of the main advantages we want to get from AutoMetaCal. We shear the galaxy and the PSF separately, reconvolve them. But we don't have to make any of the 8 additional images beyond the "noshear", because we get the derivatives of perturbations at noshear.

[EiffL]

@andrevitorelli regarding Huff and Mandelbaum, are you sure they reconvolve by an isotropized PSF?

https://arxiv.org/abs/1702.02600

They state for instance in that paper that:

If the shape measurement algorithm does not perfectly remove PSF ellipticity, then the shapes measured from shearing the PSF (e+,PSF and e−,PSF) permit calculation of at least the linear-order residual PSF ellipticity biases

Implying that 1. They use a method that removes most of the PSF contribution, with something like regauss, and are metacalibrating the residuals. 2. They don't reconvolve by an isotropized PSF because there wouldn't need to do such a correction in that case

[aguinot]

I might missing something but in your code I see that you are tracking only 1 shear tensor. If you apply the shear independently on the galaxy and the PSF shouldn't you have 2 ? or doing it in 2 passes ?

[EiffL]

there is one shear tensor but it has 4 entries, 2 for the image shear, 2 for the PSF shear

[andrevitorelli]

I still can't figure out if the correct is to do like what I'm doing it:

e = <R^-1> -

or:

e = <R_{psf_removed}^-1><e_measured - R_psf e_psf>

Note that each observed shape is obs = intrinsic + shear + psf, from which an average should give shear + psf. So, in the first eq, we're getting an debiased () estimator for shear + psf, from which we subtract the psf contribution by doing - < Rpsf > .

In the second option, which was what I thought the paper implied (e_obs - Rpsf e_psf is just Kaiser95), it would have to been nested. We'd do:

e_g = method(metacal_img(img,psf,g)) #shears only the deconvolved galaxy image
e_gs = method(metacal_psf_img(img,psf,g_psf)) #shears only the psf
e_psf = method(psf) #this seems still biased
R_psf = d e_gs / d g_psf
e = e_g - R_psf e_psf
R = de/dg

However, it's now clear to me the way I'm calculating the residual bias is wrong (after a very helpful discussion with @b-remy ). I have only one shear point which is observed averaging all 100k galaxies, so I can't distinguish m from c. I'll be fixing this today.

[EiffL]

hummm I don't understand why your second method looks like this. I do agree that you want to compute the response of the PSF-corrected ellipticity, and subtract a residual response to PSF ellipticity. So something like your second option.

So I think the pseudo code should be something like this:

# where method is a function of both image and psf (after extra shearing)
e = autometacal(image, psf, method, de_shear, de_psf) 

R = grad(e, de_shear)
R_psf = grad(e, de_psf)

e_corrected = R^{-1} (e  - R_psf e_psf)

no ?

[andrevitorelli]

I don't understand how R in your case would be a response of e - R_psf e_psf, since it's being calculated on e....

Or well, maybe I do, since d (e - Rpsf e_psf)/dg = de/dg - d(Rpsf e_psf)/dg = de/dg, since the second term does not depend on g...

[EiffL]

Yeah R_shear corrects the multiplicative bias, we can apply on an "additive bias-subtracted" e

[andrevitorelli]

So, just to report, I decided to build up the full experiment first testing on images without shape noise. I create 20 "fields" of 10k galaxies. Each of these fields have a single true shear applied to them, and I add some pixel noise. We can then 1) look at the residuals of the estimated shear both calibrated and uncalibrated, and 2) vary the PSF and repeat 1). The nb I'm using right now is here: PSF_experiment_notebook-shapenoiseless.ipynb When I get shape noise back, I'll do so by changing to the "galgen" dataset and then to the "CFIS"-like dataset.

[aguinot]

It looks very nice!
Also, it is hard to comment on the results without error bars. Is it possible to run this notebook on collab for example? Maybe @EiffL knows? (I just wonder if all the library are available, like the modification done for the interpolation)

[andrevitorelli]

Is it possible to run this notebook on collab for example?

Not at the moment, unfortunately, but I'm trying to address this.

Also, it is hard to comment on the results without error bars.

The cells that are running right now (the ones with numpy save commands), do get error bars on estimated shears, so when it finishes I intend to come with some plots to show what we're getting.

[andrevitorelli]

So, it seems that the PSF correction is wrong for some reason. I also updated the notebook in the link above.

[andrevitorelli]

No improvements since this last plot, but one nice thing. While I was doing some noiseless tests, I decided to see where the assumption of linear response broke down for our test. And this is it:

This is the same test with pixel (but not shape) noise:

[andrevitorelli]

Ok, so... a few tests seem to show me that at least for PSF ellipticities (g convention) <0.05, the R_psf is higher by around 14%, and this is stable for some orders of magnitude... I don't understand what it means, but wanted to log it here.

[aguinot]

Maybe you could try to make a plot of the measured PSF ellipticity vs the true one. If your PSF is gaussian and you are using moments the match should be perfect. We might learn something from it.

[andrevitorelli]

I have been trying to nail down how to use model fitting with ngmix, but I am doing something wrong. Can you take a quick look, @aguinot ? PSF_response/notebooks/psf_metacal_ngmix.ipynb

I think it might be some stupid bug in my code to get the R matrix* out of the result dictionary. (I know the code isn't the best here, but I did it so to do minimal changes from a code that worked before, with the gaussmom result dic, which is different...)

def get_metacal_response_ngmix(resdict):
  step=0.01

  #shear
  g0s = np.array([resdict['noshear']['g'][0], resdict['noshear']['g'][1]])
  g1p = np.array([resdict['1p']['g'][0], resdict['1p']['g'][1]])
  g1m = np.array([resdict['1m']['g'][0], resdict['1m']['g'][1]])
  g2p = np.array([resdict['2p']['g'][0], resdict['2p']['g'][1]])
  g2m = np.array([resdict['2m']['g'][0], resdict['2m']['g'][1]])    

  R11 = (g1p[0]-g1m[0])/(2*step)
  R21 = (g1p[1]-g1m[1])/(2*step) 
  R12 = (g2p[0]-g2m[0])/(2*step)
  R22 = (g2p[1]-g2m[1])/(2*step)

  #PSF
  g1p_psf = np.array([resdict['1p_psf']['g'][0], resdict['1p_psf']['g'][1]])
  g1m_psf = np.array([resdict['1m_psf']['g'][0], resdict['1m_psf']['g'][1]])
  g2p_psf = np.array([resdict['2p_psf']['g'][0], resdict['2p_psf']['g'][1]])
  g2m_psf = np.array([resdict['2m_psf']['g'][0], resdict['2m_psf']['g'][1]])    

  R11_psf = (g1p_psf[0]-g1m_psf[0])/(2*step)
  R21_psf = (g1p_psf[1]-g1m_psf[1])/(2*step) 
  R12_psf = (g2p_psf[0]-g2m_psf[0])/(2*step)
  R22_psf = (g2p_psf[1]-g2m_psf[1])/(2*step)  

  ellip_dict = {
    'noshear':g0s,
    '1p':g1p,
    '1m':g1m,
    '2p':g2p,
    '2m':g2m,
    '1p_psf':g1p_psf,
    '1m_psf':g1m_psf,
    '2p_psf':g2p_psf,
    '2m_psf':g2m_psf,    
  } 

  R = np.array(
    [[R11,R12],
     [R21,R22]])

  Rpsf = np.array(
    [[R11_psf,R12_psf],
     [R21_psf,R22_psf]])

  return ellip_dict, R, Rpsf

• or not, since the R matrix results and the measured ellipticities of the sheared images seem strange also...

[aguinot]

I don't really see something wrong in your code, appart from how you define the jacobian. The center of the jacobian doesn't need to be integers. In ngmix it is defined as (size_x + 1)/2. I don't know if that will change a lot your results though.. (it might since there is 1 pixel offset and that's the boundary if your prior) Could you post the row output of ngmix for a noiseless galaxy with the expected value to see if we can spot something?

[andrevitorelli]

Most recent results are here, after a gist by @EiffL : PSF Experiment 2.0 notebook

[andrevitorelli]

Continuing the investigations, I tested some settings between almost no to low noise, and some different levels of anisotropy of the psf.

TL;DR: We have a fairly constant bias (~0.014) in the diagonal elements of the Rpsf that account for our inability to correct for the additive bias. The shear response difference between ngmix and amc is not very meaningful, as we do get correct multiplicative biases residuals after correction.

Some points:

• No noise = 1e-6 noise on the observed images
• Low noise = 1e-4 "
• PSF images have a 1e3 lower level of noise.
• The second histograms for each test have a typo in their labels, it should be R, not Rpsf.

1. No noise, anisotropic PSF

2. No noise, isotropic PSF

3. Low noise, anisotropic PSF

4. Low noise, isotropic PSF

[andrevitorelli]

The code is here: u/andrezvitorelli/PSF_response/notebooks/psf_metacal_ngmix_moments.ipynb

[andrevitorelli]

Some additional info:

• Removing fixnoise from ngmix doesn't help much. ngmix does get worse results, but not as bad as ours and not closer to ours.
• Using a dilated PSF in our code didn't solve the isse (dilate code below)

def dilate(img,factor,interpolator="bernsteinquintic"):
  """ Dilate images by some factor, preserving the center. 

  Args:
    img: tf tensor containing [batch_size, nx, ny, channels] images
    factor: tf tensor containing dilation factor (factor >= 1), either 1 for all, or batch_size

  Returns:
    dilated: tf tensor containing [batch_size, nx, ny, channels] images dilated by factor around the centre
  """
  img = tf.convert_to_tensor(img,dtype=tf.float32)
  batch_size, nx, ny, _ = img.get_shape()

  #x
  sampling_x = tf.linspace(0.,tf.cast(nx,tf.float32)-1.,nx)[tf.newaxis]
  centred_sampling_x = sampling_x - nx//2
  batched_sampling_x = tf.repeat(centred_sampling_x,batch_size,axis=0)
  rescale_sampling_x = tf.transpose(batched_sampling_x) / factor
  reshift_sampling_x = tf.transpose(rescale_sampling_x)+nx//2
  #y
  sampling_y = tf.linspace(0.,tf.cast(ny,tf.float32)-1.,ny)[tf.newaxis]
  centred_sampling_y = sampling_y - ny//2
  batched_sampling_y = tf.repeat(centred_sampling_y,batch_size,axis=0)
  rescale_sampling_y = tf.transpose(batched_sampling_y) / factor
  reshift_sampling_y = tf.transpose(rescale_sampling_y)+ny//2

  meshx = tf.transpose(tf.repeat([reshift_sampling_x],nx,axis=0),[1,0,2])
  meshy = tf.transpose(tf.transpose(tf.repeat([reshift_sampling_y],ny,axis=0)),[1,0,2])
  warp = tf.transpose(tf.stack([meshx,meshy]),[1,2,3,0])

  dilated= resampler(img,warp,interpolator)

  return tf.transpose(tf.transpose(dilated) /factor**2)

[andrevitorelli]

I have been investigating the issue with the pixel response and the psf, as suggested by @aguinot

Notes to self, what ngmix does is this:

1. Deconvolve galaxy by psf gal/psf
2. Shear deconvolved gal s(gal/psf,g_gal)
3. Deconvolve psf from pixel response psf/pix
4. Dilate no pix psf D(psf/pix)
5. Shear dilated no pix psf s(D(psf/pix),g_psf)
6. Reconvolve pixel response to the sheared and dilated psf s(D(psf/pix))*pix
7. Reconvolve with galaxy s(D(psf/pix),g_psf)*pix*s(gal/psf,g_gal)

[andrevitorelli]

Possible good news. By removing the line that deconvolves the pixel from the PSF, ngmix and amc get similar results after correction (not before, though - it might be because of other things I needed to edit on ngmix code to fully remove all mentions of pixel response):

I just altered the line from

self.psf_int_nopix = galsim.Convolve([psf_int, self.pixel_inv])

to

self.psf_int_nopix = psf_int

[andrevitorelli]

Notes to self 2:

• galsim.Pixel is a galsim.Box with height=width=scale and flux=1
• galsim.Box is defined in SBBox.cpp