fix bug with metacal psf method "gauss"

[esheldon]

The psf determination was being done on the pre pixel rather than post pixel psf. This required a Convolve with the pixel afterward, which does not work as expected when the wcs is complex

Added a metacal accuracy test as well

[esheldon]

@beckermr is that failure obvious to you?

[beckermr]

I'll check the test thing in a bit. I still don't get why this fails or what is going on.

[esheldon]

Looks like the column det_bands was not expected

[beckermr]

So I am looking at this in detail and this code is wrong IMHO. IDK why it works, but it is not logically consistent.

The function _get_gauss_target_psf returns a galsim.Gaussian object. Then do_dilate applies a dilation based on the shear. Finally, the output of _get_dilated_psf is used in this block of code

        if key not in self._psf_cache:
            psf_grown = self._get_dilated_psf(shear, doshear=doshear)

            # this should carry over the wcs
            psf_grown_image = self.psf_image.copy()

            try:
                psf_grown.drawImage(
                    image=psf_grown_image,
                    method='no_pixel',  # pixel is already in psf
                )
            except RuntimeError as err:
                # argh, galsim uses generic exceptions
                raise GMixRangeError("galsim error: '%s'" % str(err))

            self._psf_cache[key] = (psf_grown_image, psf_grown)

Notice here the image is rendered with no_pixel. However, with the changes above, the pixel was never added back in. It doesn't matter that you kept the pixel when calling _get_gauss_target_psf since this function simply returns a raw galsim object that knows nothing of the pixel.

What I think is happening is that indeed the pixel handling has gone wrong and the change here happens to make things just right for the test case.

[esheldon]

With the new code the pixel is never taken out

[esheldon]

There was never a reason to take out the pixel before running that algorithm. In fact I think it was a bad idea to do so because the determination of the gaussian was affected by noise amplification.

[beckermr]

Sure but you do have to add the pixel in order for the rendering to be correct with no_pixel.

[beckermr]

Not taking it out doesn't matter. The returned reconv PSF is not an interpolated image. It is a Gaussian from galsim.

[esheldon]

This is the same thing we do with fitgauss. In that case we fit the pixelized PSF, dilate it and render with no pixel.

It is true that there could be problems with an undersampled PSF, but that's not unique to this method.

[beckermr]

Ahhhhhh. I think I understand now. The algorithm you have is correct, but I am not sure if it is correct for the reason stated (though I could not have understood you).

So the pixel is not "in the image" in the sense we mean normally.

The normal sense of this phrase is as follows. We mean specifically that if we fit a model to a pixelized image, then that model's values represent the value of the pixelized image at a given point. If we do this to a PSF image, then we can use this fit model as the PSF with the pixel. A similar process happens for interpolated images as PSFs in galsim. Here we could in principle decovolve the pixel and then render with galsim normally if we wanted to. In practice that procedure is not numerically stable.

What is happening in metacal for fitgauss here is different. Here in the case of fitgauss, we fit a Gaussian, dilate the model fit, and then use it to render images in Galsim using no_pixel. That dilated model fit is not the same as deconvolving the pixel from the PSF, dilating the per-pixel PSF, and then reconvolving the pixel. There is ofc formally a PSF that the dilated model represents, but that is just some other PSF.

Instead, when we use no_pixel to render the final images in fitgauss, I think what is happening is not that the actual WCS pixel is already in the images, but that instead we are simply not broadening PSF more than we need to when rendering.

So I think what you are saying is that we do not need to broaden the gauss PSF more since we now get the reconv PSF from the original PSF with a pixel.

In this case, I think what has happened is that the algorithm for making the gauss reconv PSF is not stable when run against a PSF image with the pixel removed. It must be erroneously returning a PSF that is too small by a large factor or something.

We should check that this last thing is indeed what happened before we merge.

[esheldon]

Right, when I said the pixel is in it I meant that it is broader than the pixelized PSF we started with, so should be valid.

But its not correct that the algorithm is just returning a too small PSF. Any time I reconvolve by the pixel I'm getting a bias.

For example I can just pick some good PSF by hand, say the one returned by fitgauss. If I convolve that by the pixel and use it I get a bias. This is what is going wrong.

[beckermr]

OK then there is another bug or assumption in the code we need to find.

Why would a reconv PSF that is a pixel convolved with some large Gaussian not work, but a reconv PSF that is simply that large gaussian be fine?

[beckermr]

Ahhhhhh. With a WCS that is not diagonal, the pixel is asymmetric. It also not true that the pixel is circularly symmetric even with a diagonal WCS. I am betting we have an additive effect that looks like a multiplicative one.

[esheldon]

It could be. I'd like to make that a separate issue for research, since there is no reason to think this doesn't work fine.

[beckermr]

Sounds good. Thanks for your patience with my questions!