1294 from images

[FedericoBerlfein]

I have added a classmethod to the InterpolatedChromaticObject class called from_images() that takes in as input a list of Galsim Image objects and a list of wavelengths. This will allow the user to created aInterpolatedChromaticObject directly from images. The method also takes as arguments a list of stepk and maxk values if the user wishes to pass these directly to the InterpolatedImage object that is created within the method. Any other keyword that the InterpolatedImage class accepts can in addition be passed as kwargs to the method. In addition, unit tests have been added to the file test_chromatic.py that compare this initialization method to the "default" method. Any feedback/suggestions are highly appreciated!

[rmandelb]

Hi Federico - thanks, this looks quite good to me. I have a few suggestions that you can commit directly to fix minor issues with the docstrings. The only substantive comments are about places where the code makes assumptions that don't seem to be documented or tested. Other than that, the set of unit tests looks good to me.

[rmandelb]

I think we should check whether they are actually sorted and raise an exception if they aren't. Should also document this requirement in the docstring.

[rmjarvis]

A nice API we use in GalSim a lot is to have the main function start with doing any sanity checks like this and then call a "leading underscore" version of the function. So InterpolatedChromaticObject._from_images in this case, which does the rest of the work. Then if a user wants to avoid the sanity checks for efficiency reasons when they know they are all guaranteed to pass, then the can just directly call _from_images instead.

[FedericoBerlfein]

Thanks for the suggestion! I I have added this main function do the sanity checks and then have it call the "leading underscore". I haven't pushed these changes yet due to a question about the other comment about the stepk, maxk.

[rmandelb]

Are we assuming all images have the same pixel scale? If so, that needs to be documented and potentially tested.

[FedericoBerlfein]

Clarification added to documentation and tested in unit tests.

[rmjarvis]

Looks good. A few comments.

[rmjarvis]

A nice API we use in GalSim a lot is to have the main function start with doing any sanity checks like this and then call a "leading underscore" version of the function. So InterpolatedChromaticObject._from_images in this case, which does the rest of the work. Then if a user wants to avoid the sanity checks for efficiency reasons when they know they are all guaranteed to pass, then the can just directly call _from_images instead.

[rmjarvis]

These aren't right I think. Or at least, they are much simpler (so both less conservative in some respects but also less efficient in others) than what InterpolatedImage does to compute these. See the functions InterpolatedImage._getStepK and _getMaxK for the more sophisticated calculation of these. Most importantly, this stepk will usually be much smaller than required, and probably maxk will be larger than required. This will make subsequent draw operations much slower than they need to be.

Some of the calculation of this can be a little slow, since it needs to do some calculations based on the actual profile in the image to get the right values. But probably the resulting values for stepk and maxk are going to be nearly the same for all the images in something like this, so we can try to just do it once and use the same values for all the images here.

So I think you could make a regular InterpolatedImage from either the first image in the list, or maybe the sum of all of them (to get an average profile), let it compute stepk and maxk in its way with the sophisticated calculation, and then use those values of stepk and maxk for all of the objects here using the _force parameters. That would probably work well and save a bunch of nearly equivalent calculations.

[FedericoBerlfein]

Thanks for the explanation Mike. We initially chose to make this the default stepk, maxk because we had noticed some potential issues when using the stepk, maxk the InterpolatedImage calculates. The problem we were seeing was that when drawing images at a different pixel scale from the input images given to the from_images function, the image was being truncated after a certain pixel length, causing to truncate some features of the PSF (such as the spikes). As an example using the roman.getPSF, I did quick test to show this potential problem and would love to hear your opinions on this. I essentially test the from_images method using the stepk and maxk obtained from creating an average profile of the input images and letting InterpolatedImage calculate the stepk and maxk from this image to use for all other input images. I then test the "conservative" stepk, maxk from the code snippet of this comment, and compare both from using the roman.getPSF directly. The image below shows the differences and the cut I had mentioned before. In specific, I noitced that this truncation is related to the difference in stepk, not maxk. Let me know what you think of this, if needed I can also provide the exact code I used to generate these images, I haven't pushed the changes made to include the stepk, maxk from the average profile yet as I am still testing it here.

[FedericoBerlfein]

Maybe to explain in a bit more detail what went into making the image above. I started from getting the PSF for the Y band in Roman and then generating individual oversampled images at specfic wavelengths. Then I create the InterpolatedChromaticObject using the two different ways of calculating stepk, maxk I described above. Note that in this version of the code I am using locally, I changed the default to be the average profile method proposed by Mike. Finally, I draw the all images at a different pixel scale from the one provided to the images:

filt = 'Y106'
bandpass = roman_filters[filt]
# get roman PSF
psf_chrom = roman.getPSF(1, filt, n_waves = 10,  pupil_bin=8)
waves = psf_chrom.waves
psf_images = []
# get oversampled images at set wavelengths
for wave in waves:
    psf_wave = roman.getPSF(1, filt, wavelength=wave,  pupil_bin=8)
    #psf_chrom = roman.getPSF(1, filt, n_waves = 5,  pupil_bin=8)
    psf_img = psf_wave.drawImage(nx=480, ny=480, scale=scale, method = 'no_pixel')
    psf_images.append(psf_img)

N = 480
# default stepk and maxk from pixel scale
stepk  = 2*np.pi/(scale*N)
maxk = np.pi/scale
    
# create flat SED and make interpolated object using default or conservative stepk, maxk
flat_sed = galsim.SED('1',wave_type = 'nm', flux_type = 'flambda' )
int_obj = galsim.InterpolatedChromaticObject.from_images(psf_images,  waves)
int_obj_conservative = galsim.InterpolatedChromaticObject.from_images(psf_images,  waves,
                                                                      _force_stepk =stepk,_force_maxk =maxk )
conv_obj = galsim.Convolve(int_obj, galsim.DeltaFunction()*flat_sed)
conv_obj_conservative = galsim.Convolve(int_obj_conservative, galsim.DeltaFunction()*flat_sed)

# draw images using 
drawing_scale = roman.pixel_scale
test_img = conv_obj.drawImage(bandpass, nx=100, ny=100, scale=drawing_scale, method = 'auto').array
test_img /= np.sum(test_img)

test_img_conservative = conv_obj_conservative.drawImage(bandpass, nx=100, ny=100, scale=drawing_scale, method = 'auto').array
test_img_conservative /= np.sum(test_img_conservative)

# true PSF using getPSF() object directly
true_obj = galsim.Convolve(psf_chrom, galsim.DeltaFunction()*flat_sed)
true_img = true_obj.drawImage(bandpass, nx=100, ny=100, scale=drawing_scale, method = 'auto').array
true_img /= np.sum(true_img)

[rmjarvis]

The knob we have in GalSim to control the level of approximation used for stepk is folding_threshold. Smaller values mean large area in real space.

I think what you are effectively doing in your "conservative" version is the equivalent of lowering folding_threshold significantly. If you try making the InterpolatedImage with gsparams=galsim.GSParams(folding_threshold=1.e-3), I suspect you will get similar results.

The advantage of conforming to this paradigm is that users can then be consistent in terms of the level of folding in k-space for the PSF and galaxy profiles when you convolve them together, rather than have one of them much more precise than the other, which would be wasted precision. If you set folding_threshold for a Convolve, then it gets applied to all constituent profiles.

See also the docs about folding_threshold for the Roman PSF here, and also this PR where I show examples of the PSF with different folding_threshold and pupil_bin values.

[FedericoBerlfein]

Yes that solved the issue, thanks Mike! In that case, I will keep the default to be the stepk and maxk from the average profile from your suggestion. Do you think its worth adding some documentation about setting folding_threshold when using this method of initialization? If so I can add it.

[rmjarvis]

I don't think the reasons for using it here are different from other places in GalSim. It's only going to be a cosmetic difference for very bright things. When used with typical galaxy profiles, the default folding_threshold will be just fine for the vast majority of use cases. E.g. in the Rubin-Roman sims, Troxel had a special setup for the bright stars so they looked pretty. We have something similar in imsim. But when used for galaxies, the regular gsparams are fine.

[FedericoBerlfein]

Got it, I will try to push the changes discussed today.

[FedericoBerlfein]

I just pushed the changes discussed here. Mainly, I changed the main function from_images to do some sanity checks first (check pixel scale is the same for all images, check dimensions for all images match) and do sorting of the wavelengths and images. It then calls the _from_images function which has the actual functionality. It now calculates the default stepk and maxk from the average profile of the input images.

[FedericoBerlfein]

@rmjarvis I noticed that the code does not hit 100% coverage on codecov, however, the line that its telling me is not covered is a commented line. Not sure why this is happening or how to fix it. Let me know if you have any thoughts on this. I have also made changes based on previous conversations that are ready for review whenever you have time, let me know of any further comments or changes!

[rmjarvis]

I'm not sure why codecov is reporting the wrong line, but when I run coverage locally, this is the line that isn't covered.

[FedericoBerlfein]

I see, thats odd, I test those errors in lines 1970-1980 of test_chromatic.py. Not sure why it says they are not covered.

[rmjarvis]

The trick I use when building up coverage is to add a line before the one I want to check is covered that looks like:

raise '1'

Then if I don't get a TypeError: exceptions must derive from BaseException, then I know it never got to that point.
Then you can keep adjusting the test until you get that error, at which point remove the added line, and you'll know it's getting to the line you want covered.
If there doesn't seem to be any way to hit this line, then it may be that the condition you are testing is already implicitly tested elsewhere, in which case this test is unnecessary and can be removed.

[FedericoBerlfein]

Thanks that was helpful in finding the issue. It should be solved now.

[rmjarvis]

FWIW, our assert_raises is identical to np.testing.assert_raises, so you could have left it as the context version if you prefer that for clarity. But this is fine too.

[rmjarvis]

If someone uses images with a non-PixelScale wcs, this will give a possibly obtuse error message. Maybe worth adding:

if images[0].wcs is None or not images[0].wcs._isPixelScale:
    raise GalSimValueError("images must have a pixel scale wcs.")

[FedericoBerlfein]

Added.

[rmjarvis]

I think both of these should be GalSimValueError, not GalSimNotImplementedError.
cf. our description of when to use each error type at the top of errors.py. I think both of these feel more like something that will never work, rather than something that could be implemented down the road.

[FedericoBerlfein]

Corrected.

[rmjarvis]

This should have a doc string, as it is not just an implementation detail, but rather something that users can be expected to be stable API. Just more dangerous than the non-underscore version in terms of the inputs. Our usual docstring would be something like

Equivalent to InterpolatedChromaticObject.from_images, but without the sanity checks.  
Also, the waves must be given in sorted order.

[FedericoBerlfein]

Added.

[rmjarvis]

Suggested change

	wave: list of wavelength values in nanometers.
	waves: list of wavelength values in nanometers.

[FedericoBerlfein]

Resolved.

[rmjarvis]

Looks good. Thanks Federico.