Make the RomanOptics chromatic=True option more effective and more efficient. by rmjarvis · Pull Request #198 · rmjarvis/Piff

rmjarvis · 2026-04-05T11:05:16Z

This PR is focused on the implementation of the chromatic=True option for RomanOptics. It includes the following pieces:

Added the bandpass specification in input, rather than as a detail of the RomanOptics PSF.
If the bandpass is a galsim.roman version, then RomanOptics class can get the filter name from that. If not (e.g. if we want to specify a different version that is more accurate perhaps), then filter_name can be specified manually for use in the getPSF call.
The input SED is immediately multiplied by the bandpass, and that combination is thinned directly, rather than only thinning the SED and still using the unthinned Bandpass, which was still quite slow.
Since we've already multiplied by the filter throughput, the draw command needs to use a flat version of the bandpass with just the end points, and throughput = 1 in between. This is called flat_bandpass in the code.
I added a more aggressive kind of thinning, which specifies a maximum number of samples to use in the SED. E.g. sed_max_samples=3 keeps only 3 points, so piecewise linear with 2 segments. I think given how uncertain our SED knowledge will be in practice, this will be acceptable on real data. Either sed_tol or sed_max_samples or both are available for us to try out on the data.
Fixed a few bugs that turned up related to using ChromaticObjects in various places.
Fixed an error in both RomanOpticalModel and OpticalModel. They should both have _centered=True, not False. With _centered=False, they don't re-centroid properly with uncertain star positions. This hadn't turned up, since both of them are usually used in combination with other components (e.g. in a SumPSF), which allows for recentroiding. It turned up in one of my tests here, but the error is unrelated to the chromatic option.

The config file in devel/roman using chromatic=True is about 5-6 times slower than the achromatic version. This is with just a linear approximation to the SED (sed_max_samples=2) -- it's almost double that with 3 samples. This seems reasonable enough given the extra work intrinsic to drawing chromatically, but there might be further optimizations available to try. E.g. I haven't tried using photon shooting (for either one), which might be worth trying despite the fact that it results in noisy images (which makes the empirical derivatives tricky, but there might be ways around that). I also don't have a good feel for how much the different SED approximations matter in practice, since we haven't tested yet with realistic SEDs.

I'll also reiterate that even with chromatic=False, the final PSF model is chromatic. Just the fitting is done achromatically with the PSF being evaluated at just a single wavelength (the bandpass effective wavelength). I'm open to adjustments to the naming of this parameter, since it's potentially confusing as it currently stands.

… #1355

…andpass

This reverts commit 01eace8.

HyeongHan · 2026-04-05T19:16:09Z

+            bp_wave_native = _bandpass_native_waves(newx[in_band])
+            tp_native[in_band] = bandpass._tp(bp_wave_native) / bandpass.wave_factor
+            nz = tp_native != 0.  # Don't divide by 0.
+            assert np.all(newf[~nz] == 0.)


It's an temporary compatibility patch, but for clarity an error message can be added something like "thinned sed has non-zero values where throughput is zero"

Comments on this file should probably be migrated to GalSim-developers/GalSim#1355.
(I no longer have this patch in the PR.) But also, I don't think it's possible for this assert to trigger. That's why it's an assert, not a normal error.

HyeongHan · 2026-04-05T19:20:28Z

+        # Note that this is thinning in native units, not nm and photons/nm.
+        interpolant = (self.interpolant if not isinstance(self._spec, LookupTable)
+                       else self._spec.interpolant)
+        newx, newf = utilities.thin_tabulated_values(


it can be renamed as thinned_lambda thinned_flux but if this module is going to be removed anyway please ignore this comment

HyeongHan · 2026-04-05T19:35:49Z

-        self.wcs = wcs
-        self.pointing = pointing
+        if self.wcs is None:
+            logger.error("WARNING: wcs and pointing should now be given in process, not fit.")


If it's warning, shouldn't it be logger.warning?
If it's an error, the message should be ERROR.
self.pointing is None condition also be added?

All such messages in Piff are done through the logger mechanism, not python's warnings module. Given that, logger.error vs logger.warning dictates what verbosity level it shows up. logger.error's are always written, which I think makes sense for deprecation warnings. The user should fix it. Warnings that are merely about something going a little weird in the processing, but might be ok are logger.warning.

HyeongHan · 2026-04-05T20:17:51Z


+    def set_context(self, wcs, pointing, bandpass):
+        super().set_context(wcs, pointing, bandpass)
+        if isinstance(self.components, list):


if self.components is not None: can be more straight forward.

None wouldn't work. It's either a list or an int. The latter is when it is still in the process of being read in from a file, but _finish_read hasn't been called yet. The other check for whether the read is complete uses this, so I'm inclined to keep it this way for consistency.

nihardalal

Looks good from my end - I think I have a couple questions/comments here and there, feel free to answer as/if necessary

nihardalal · 2026-04-20T15:35:48Z

        return self._interpolate_samples(star, sample_profiles)

-    def _draw_model_image_from_samples(self, star, sample_profiles, convert_func=None):
+    def _draw_model_image_from_samples(self, star, sample_profiles, convert_func):


Why is convert_func a mandatory argument now?

I generally prefer my backend functions (I.e. not ones that users would use) to not use default arguments. That way when I change anything about them, I immediately catch any uses where I forgot to update to a new signature, rather than silently using a default value that isn't appropriate. So this should not have used a default =None in the first place.

nihardalal · 2026-04-20T15:39:04Z

+        self.bandpass = bandpass
+        self.flat_bandpass = make_flat(bandpass)
+
+    def _require_bandpass(self):


Maybe consider renaming this method to _get_bandpass? This is just a matter of preference, but require suggests a bool output to me, and we use this more like a get

The only thing this does is validate that the PSF object has a bandpass set (which happens upstream of the Roman-specific processing). Returning that value is secondary. Maybe the better change is just to remove the return value.

I think it is called in a place where it does actually return the value, so feel free to leave as is.

nihardalal · 2026-04-20T15:41:14Z

        # Only call set_num if they are actually built.
        if isinstance(self.model, Model):
            self.model.set_num(num)
-        if isinstance(self.interp, Interp):


Isn't the interp instance also potentially not real during the read?

Yes, but they are real or not real in tandem. So I simplified this to a single check. If Model is real, so is interp, and they can both do their set_num now. Otherwise they are both fake and will be dealt with by the read function.

nihardalal · 2026-04-20T15:55:17Z

+        g = np.linalg.lstsq(Aw, fw, rcond=None)[0]
+
+        # g is now the vector that gives the best least squares approximation f'.
+        return g


Makes sense to me!

nihardalal · 2026-04-20T15:58:08Z

+            resid[selected] = -1.0  # Don't reselect something we already have.
+            bisect.insort(selected, np.argmax(resid))
+
+        return candidate_wave[selected]


Also looks good to me!

nihardalal · 2026-04-20T16:01:32Z

+        wave = np.array(sed.wave_list, dtype=float)
+        if len(wave) == 0:
+            # If the sed doesn't have a wave_list, we can't thin it.  (This is very rare!)
+            sed = sed.withFlux(1.0, flat_bandpass)


I think this could potentially be a failure mode that merits some extra documentation? I can imagine us accidentally having SEDs saved in an incorrect format and then this being some weird behavior that shows up (but I might be overthinking)

It's pretty rare, but I think it's not something that would merit an error. It requires both the SED and the bandpass to be lambda expressions, rather than lookup tables read from a file. This is technically allowed, although it would require a user doing some very specific special case work. There is a test that covers this branch (https://github.com/rmjarvis/Piff/pull/198/changes#diff-4650d01edc49b02fbf29d980f24442b85e599827c64340d1ec43820d757e511eR2246), but I'm not really worried about this happening by accident in realistic scenarios.

nihardalal · 2026-04-20T18:18:25Z

Tests look good to me here!

nihardalal · 2026-04-20T18:23:48Z

+            assert 0 < np.abs(star.fit.center[0]) < 0.05  # centroid shift should be small.
+            assert 0 < np.abs(star.fit.center[1]) < 0.05
+
+


Tests looks good and comprehensive

rmjarvis added 18 commits April 3, 2026 15:05

Add mokey-patch to get the bandpass option to SED.thin from GalSim PR…

b30096f

… #1355

Add input.bandpass option and store in PSF class with wcs and pointing.

b237700

Add PSF.set_context as a cleaner pattern for setting wcs, pointing, b…

82267c8

…andpass

Store seds as unit-flux normalized, which is normally what we'll want.

c387637

Add sed_tol option

7fcf9eb

Have Roman PSF (normally) get filter_name from input bandpass.

7bbdcca

Construct effective sed = sed*bandpass to thin and store

fdb9f12

Add a real test of the accuracy of chromatic RomanPSF fitting

a1ec916

Need to use Add, not Sum when components can be chromatic

933e6ce

Make reflux work correctly for chromatic objects

b60a3d6

Use a picklable flat_bandpass

990b44a

Update piff.yaml for dev run with separate files for chrom/achrom

022545e

Remove galsim_patch.py, since no longer necessary after 3be3599

94f2299

This reverts commit 01eace8.

Both OpticalModel and RomanOpticalModel should have _centered=True

56e2ddf

Add test of chromatic psf using reflux

30815fb

Add sed_max_samples option

766f158

this is the best chromatic run so far

312cb3e

typos

921ba6f

HyeongHan reviewed Apr 5, 2026

View reviewed changes

HyeongHan approved these changes Apr 7, 2026

View reviewed changes

coverage

fb6e1d2

nihardalal approved these changes Apr 20, 2026

View reviewed changes

rmjarvis merged commit cc23aaf into main Apr 24, 2026
9 checks passed

rmjarvis deleted the roman4 branch April 24, 2026 21:29

		assert 0 < np.abs(star.fit.center[0]) < 0.05 # centroid shift should be small.
		assert 0 < np.abs(star.fit.center[1]) < 0.05

Conversation

rmjarvis commented Apr 5, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

HyeongHan Apr 5, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

nihardalal left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

HyeongHan Apr 5, 2026 •

edited

Loading