Drawing error for ChromaticConvolution with pixel response

[FedericoBerlfein]

Noticed that Galsim throws an error when trying to draw a convolution between a chromatic PSF and the pixel response function. Here is a code snippet showing the issue:

import galsim
import galsim.roman as roman

filt = 'F184'
bandpass = roman.getBandpasses()[filt]
vega_sed = galsim.SED('vega.txt', 'nm', 'flambda')
scale = roman.pixel_scale/4
psfInt = roman.getPSF(8, filt, n_waves=10,  pupil_bin=8)
pixel_response = galsim.Pixel(roman.pixel_scale)
star = galsim.DeltaFunction()*vega_sed
galaxy = galsim.Gaussian(sigma=0.2)*vega_sed
eff_psf = galsim.Convolve(psfInt, pixel_response) 
star_convolved = galsim.Convolve(star, eff_psf)
galaxy_convolved = galsim.Convolve(galaxy, eff_psf)
star_image = star_convolved.drawImage(bandpass, nx=128, ny=128, scale=scale, method = 'no_pixel')
galaxy_image = galaxy_convolved.drawImage(bandpass, nx=128, ny=128, scale=scale, method = 'no_pixel')

The vega SED here is simply used as an example SED, but running this will throw the same error for the bottom two lines:
AttributeError: 'Pixel' object has no attribute '_fiducial_profile'.

Curiously there is a practical workaround for to avoid this error if you replace the PSF as a sum of the PSF at discrete wavelengths:

waves = np.linspace(bandpass.blue_limit, bandpass.red_limit, 10)
psfInt = roman.getPSF(8, filt, wavelength=waves[0],  pupil_bin=8)
for wave in waves[1:]:
    psf = roman.getPSF(8, filt, wavelength=wave,  pupil_bin=8)
    psfInt += psf

which does not result in an error. However, I don't think this implementation is correct.

Any comments or suggestions are greatly appreciated!

[rmandelb]

Just to elaborate slightly on the goal here:
We want to draw oversampled images for Roman (with smaller pixel scale than the native one) but with the correct pixel response for the native Roman pixel scale. The way this code is structured is meant to mirror demo13.py in terms of how it associates SEDs with objects (i.e., with the star and galaxy profiles, then convolve with a PSF from getPSF()) but because of the interest in drawing oversampled version, the pixel response is handled differently: here we convolve directly and then use drawImage in no_pixel mode, whereas there we use the implicit convolution that we get when using drawImage in auto mode. I was surprised that the above code snippet doesn't work.

Suggestions are welcome: is this code misuse or a bug in GalSim?

[rmandelb]

After some discussion I went back to confirm that the achromatic version of this process works - below I've made the PSF achromatic and did not multiply the star or galaxy by an SED, so they'd be achromatic as well, and the code runs without complaint:

filt = 'F184'
bandpass = roman.getBandpasses()[filt]
scale = roman.pixel_scale/4
psfInt = roman.getPSF(8, filt, wavelength=bandpass, pupil_bin=8) # achromatic, defined at effective wavelength
pixel_response = galsim.Pixel(roman.pixel_scale)
star = galsim.DeltaFunction()
galaxy = galsim.Gaussian(sigma=0.2)
eff_psf = galsim.Convolve(psfInt, pixel_response) 
star_convolved = galsim.Convolve(star, eff_psf)
galaxy_convolved = galsim.Convolve(galaxy, eff_psf)
star_image = star_convolved.drawImage(nx=128, ny=128, scale=scale, method = 'no_pixel')
galaxy_image = galaxy_convolved.drawImage(nx=128, ny=128, scale=scale, method = 'no_pixel')

[ztq1996]

I was getting around this issue with the following code:

filters = roman.getBandpasses()
flat_sed = galsim.SED(lambda x:1, 'nm', 'flambda').withFlux(1.,filters["J129"]) 
tophat = galsim.Pixel(wfirst.pixel_scale,flux = 1.0)
PSF = roman.getPSF(1, "J129", n_waves=10)

final = galsim.Convolve(tophat*flat_sed,PSF)

oversample_image = final.drawImage(bandpass=filters["J129"], scale=wfirst.pixel_scale/8, nx = 64, ny = 64, method = 'no_pixel')

Basically, the trick is to apply the SED to the pixel response function, rather than the galaxy or star... Obviously the SED belongs to the pixel, not the galaxy, Duh

[ztq1996]

Similarly, if you want to make a chromatic galaxy image, you have to do the convolution twice, and apply the SED at the second convolution:

gal = galsim.Gaussian(sigma = 0.2).shear(e1 = 0.0, e2 = 0.3)
tophat = galsim.Pixel(wfirst.pixel_scale,flux = 1.0)
PSF = roman.getPSF(1, "J129", n_waves=10)
flat_sed = galsim.SED(lambda x:1, 'nm', 'flambda').withFlux(1.,filters["J129"]) 


gal_tophat = galsim.Convolve(tophat,gal)
gal_tophat_PSF = galsim.Convolve(gal_tophat*flat_sed,PSF)
gal_image = gal_tophat_PSF.drawImage(bandpass=filters["J129"], scale=wfirst.pixel_scale/8, nx = 64, ny = 64, method = 'no_pixel')

Note that gal can be any kind of galaxy profile. If you need to do bulge+disk with different SEDs, do it per component:

bulge = galsim.Sersic(half_light_radius = 0.2, n = 4)
disk = galsim.Sersic(half_light_radius = 0.5, n=1) 

bulge_tophat = galsim.Convolve(tophat,bulge)
disk_tophat = galsim.Convolve(tophat,disk)

gal_tophat_PSF = galsim.Convolve(disk_tophat*flat_sed,PSF) + galsim.Convolve(bulge_tophat*flat_sed,PSF)

[rmjarvis]

I believe I fixed this in #1306.

[FedericoBerlfein]

Thanks for looking into this Mike! I'm glad to see the solution was simple.

[rmandelb]

Great - thank you so much for digging into this, Mike!

@FedericoBerlfein have you been able to confirm that your code now works on this branch, and the results make sense?

[FedericoBerlfein]

Yes I just checked and everything was working well. One thing that is worth noting is that the previous workaround mentioned by @ztq1996 for drawing the PSF is not equivalent to using this change. Meaning that convolving the Pixel response times the SED with the PSF is not equivalent to convolving the PSF with the pixel response and then convolving with the DeltaFunction() times the SED. To summarize, if you were to run the code below, Approach 1 and Approach 3 give the same output image, but not Approach 2.

import galsim
import galsim.roman as roman
filt = 'F184'
bandpass = roman.getBandpasses()[filt]
vega_sed = galsim.SED('vega.txt', 'nm', 'flambda')
scale = roman.pixel_scale/4
psfInt = roman.getPSF(8, filt, n_waves=10,  pupil_bin=8)

## Approach 1: Corrected approach from #1302 branch
pixel_response = galsim.Pixel(roman.pixel_scale)
star = galsim.DeltaFunction()*vega_sed
eff_psf = galsim.Convolve(psfInt, pixel_response) 
star_convolved = galsim.Convolve(star, eff_psf)
star_image = star_convolved.drawImage(bandpass, nx=128, ny=128, scale=scale, method = 'no_pixel')

## Approach 2: Previous workaround for stars
pixel_response = galsim.Pixel(roman.pixel_scale)
star_convolved = galsim.Convolve(pixel_response*vega_sed, psfInt)
star_image_2 = star_convolved.drawImage(bandpass, nx=128, ny=128, scale=scale, method = 'no_pixel')

## Approach 3: Correct workaround for stars
pixel_response = galsim.Pixel(roman.pixel_scale)
star = galsim.DeltaFunction()
star_pixresp = galsim.Convolve(star, pixel_response)
star_convolved = galsim.Convolve(star_pixresp*vega_sed, psfInt)
star_image_3 = star_convolved.drawImage(bandpass, nx=128, ny=128, scale=scale, method = 'no_pixel')```

[rmjarvis]

I'm confused. I think all three of those are identical. (On branch '#1302' that is.)

print(np.mean(star_image_2.array - star_image.array))
print(np.mean(star_image_3.array - star_image.array))

prints

0.0
0.0

What am I missing?

[FedericoBerlfein]

Ah yes you are right! Ignore my previous comment. When I was comparing on my notebook I was using a variable called star_image2 rather than star_image_2 from a different experiment I did that isn't shown on the code snippet. Silly mistake. All three approaches are equivalent. Sorry for the confusion.