convolution with rotated PSF gives different results

[esheldon]

If I convolve with a round gaussian PSF I get a different answer from when I convolve with the same round PSF rotated by 90 degrees

import galsim
import numpy as np


def compare_images(
    im1, im2,
    label1=None,
    label2=None,
    title=None,
    cmap='inferno',
    fname=None,
    show=False,
):
    import numpy as np
    import matplotlib.pyplot as mplt
    from mpl_toolkits.axes_grid1 import make_axes_locatable

    with mplt.style.context('dark_background'):
        fig, axs = mplt.subplots(ncols=2, nrows=2)
        if title is not None:
            fig.suptitle(title)

        if label1 is not None:
            axs[0, 0].set_title(label1)
        if label2 is not None:
            axs[0, 1].set_title(label2)

        if label1 is not None and label2 is not None:
            difftit = f'{label2} - {label1}'
        else:
            difftit = 'diff'

        axs[1, 0].set_title(difftit)
        axs[1, 1].set(xlabel=difftit)

        diff = im2 - im1
        maxdiff = np.abs(diff).max()
        diffrng = (-maxdiff, maxdiff)

        for ax, im, rng in zip(
            [axs[0, 0], axs[0, 1], axs[1, 0]],
            [im1, im2, diff],
            [None, None, diffrng],
        ):
            kw = {}
            if rng is not None:
                kw['vmin'] = rng[0]
                kw['vmax'] = rng[1]

            cim = ax.imshow(im, cmap=cmap, **kw)
            divider = make_axes_locatable(ax)
            cax = divider.append_axes("right", size="5%", pad=0.05)
            fig.colorbar(cim, cax=cax)

        bins = np.linspace(-maxdiff, maxdiff, 30)
        difftit = f'{label2} - {label1}'
        axs[1, 1].hist(
            diff.ravel(), bins=bins, color='sandybrown',
        )

        axs[1, 1].set(xlabel=difftit)

        fig.tight_layout(h_pad=1)
        if show:
            mplt.show()

        if fname is not None:
            fig.savefig(fname)

        mplt.close(fig)


def main():

    scale = 0.2
    flux = 1
    stamp_size = 35

    dims = [stamp_size] * 2

    angle90 = galsim.Angle(90 * galsim.degrees)

    psf = galsim.Gaussian(fwhm=0.8)
    psf_rot = psf.rotate(angle90)

    obj0 = galsim.Gaussian(fwhm=1.2, flux=flux)

    obj_convolved = galsim.Convolve(psf, obj0)
    obj_convolved_rot = galsim.Convolve(psf_rot, obj0)

    dtype = np.float64
    image1 = obj_convolved.drawImage(
        ny=dims[0],
        nx=dims[1],
        scale=scale,
        dtype=dtype,
    ).array

    image2 = obj_convolved_rot.drawImage(
        ny=dims[0],
        nx=dims[1],
        scale=scale,
        dtype=dtype,
    ).array
    print('image1 tot:', image1.sum())
    print('image2 tot:', image2.sum())
    print('tot diff:', image2.sum() - image1.sum())

    compare_images(image1, image2, fname='compare.png')


main()

[jmeyers314]

I bet it's within the guarantees of kvalue_accuracy. Swapping in

gsp = galsim.GSParams(kvalue_accuracy=1e-8)

obj_convolved = galsim.Convolve(psf, obj0, gsparams=gsp)

obj_convolved_rot = galsim.Convolve(psf_rot, obj0, gsparams=gsp)

yields

[rmjarvis]

For some value of "inconsistent" I suppose. But yeah, they are consistent with the accuracy guarantees of GalSim, such as the one Josh suggested.

[esheldon]

What's happening behind the scenes that results in a difference here?

[rmjarvis]

Probably our use of one of these two values. Trying to avoid exp when possible.

https://github.com/GalSim-developers/GalSim/blob/releases/2.6/src/SBGaussian.cpp#L62