Skip to content

Add two options for how aberrations vary over SCAs#197

Merged
rmjarvis merged 17 commits intomainfrom
roman3
Apr 3, 2026
Merged

Add two options for how aberrations vary over SCAs#197
rmjarvis merged 17 commits intomainfrom
roman3

Conversation

@rmjarvis
Copy link
Copy Markdown
Owner

@rmjarvis rmjarvis commented Apr 1, 2026

The PR adds two features to the RomanOptics PSF class.

  1. nominal_interp = "five_points". The project supplied nominal aberrations on each SCA are supplied at 5 points for each SCA: the four corners plus the center of the SCA. The default pattern for approximating the PSF in Piff is to just use the 4 corner values and take the PSF at any other position to be a weighted sum of these. @nihardalal pointed out that we could use all 5 locations where the project gives us information, including the central point instead. So setting this parameter to five_points rather than bilinear (still the default) does that.

Piff does this a little bit different from what GalSim does when calculating the PSF at an arbitrary location. GalSim does a bilinear interpolation of the aberrations themselves, and then computes the PSF from that. Piff computes the 4 or 5 PSFs with the project aberration values and then interpolates the whole PSF. I don't really have a good feel for which is more accurate. But the way Piff does it is computationally much more tractable, since constructing the PSF is relatively expensive, so using the interpolated (either bilinear or five_points) PSF is much faster.

In my test exposure, bilinear performed slightly better than five_points, but there wasn't much difference. We'll see how this goes on the real data eventually. I think this will at least be something we'll want to be able to test as to which one gives better results.

  1. aberration_interp = linear. In addition to the nominal aberrations, Piff fits for an offset to these. An "extra_aberrations" array. The default is to have a single extra_aberrations array for each SCA. This is aberration_interp="constant", and that is the default. This PR adds the option to have them be a linear function of x,y across the SCA. This option gave a worse fit on the test exposure I've been using, so constant is still the default, but we will definitely want to have this option available to us when testing the models on real data. The issue is that there are 3x as many degrees of freedom with this fit, and the data aren't strongly constraining. (The prior needs to be quite tight to avoid it going off the rails.) But with real data, we may be able to have very tight priors on linear patterns based on an ensemble of exposures, in which case this option may become more viable.

@nihardalal nihardalal self-requested a review April 2, 2026 01:02
@rmjarvis
Copy link
Copy Markdown
Owner Author

rmjarvis commented Apr 2, 2026

By the way, no need to do a careful code review of the two diagnostic scripts in the devel/roman directory. These are scripts that Codex wrote for evaluating the quality of the model solutions for the piff.yaml runs there.

Copy link
Copy Markdown
Collaborator

@nihardalal nihardalal left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good to go. I think you caught all the variable renames, and the unit tests seem pretty good to me as well. As I mentioned in the in-line reviews, we should definitely do more testing of five-points vs bilinear (either on other exposures or real data)

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Did a brief scan, looks fine!

Comment thread devel/roman/piff.yaml

min_snr: 50 # reject very faint stars
max_snr_weight: 500 # don't over-weight very bright stars
#max_mask_pixels: 10 # how many pixels in stamp are allowed to be masked.
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Did you mean for this comment to be removed?

Copy link
Copy Markdown
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, I had tried it a long time ago, and it wasn't helpful. So I commented it out.
I finally decided it wasn't worth keeping here as an advertisement that it was possible. My guess is we won't want it for Roman. I haven't been using it for this test run.

Comment thread piff/roman/roman_psf.py
def clear_cache(self):
self._corner_cache = {}

def _get_roman_five_point_data(self, sca):
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Makes sense to me!

Comment thread piff/roman/roman_psf.py
for i, (pt, p) in enumerate(zip(points, sample_params)):
extra = self._make_extra_aberrations(p)
if self.nominal_interp == 'five_point' and i == 4:
extra = extra + self._get_roman_five_point_data(sca)['center_delta']
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks correct to me - indexing seems fine/taken proper care of

Comment thread piff/roman/roman_psf.py
fx * fy,
)

def _five_point_weights(self, star):
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks correct!

Comment thread tests/test_roman.py


@timer
def test_five_point_weights():
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This test looks fine to me!

Comment thread tests/test_roman.py
),
)
fitted = model.fit(fit_star)
prof = model.getProfile(truth_params, star=star)
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a lot more compact and readable, good change

Comment thread tests/test_roman.py
aberration_prior_sigma=[0.1, 0.2, 0.3],
)
np.testing.assert_allclose(linear_vec.prior_sigma, np.tile([0.1, 0.2, 0.3], 3))

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks like a comprehensive test set

Comment thread tests/test_roman.py

@timer
def test_roman_fit_many():
def test_linear_prior_io():
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good to me!

Comment thread tests/test_roman.py

# With direct GalSim truth generation, both nominal interpolation options should converge
# to reasonable fits, but neither is expected to recover truth_params exactly
# for the reason mentioned above. (Profile interpolation != aberration interpolation.)
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should definitely do some more investigation into this (worse fit but better chi^2) in the near future

Copy link
Copy Markdown
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This wasn't the case in the big devel/roman test. There the overall chisq was worse with five_point. Not outrageously so, but enough that I left the canonical value there as bilinear.

@rmjarvis
Copy link
Copy Markdown
Owner Author

rmjarvis commented Apr 3, 2026

Thanks Nihar! Merging.

@rmjarvis rmjarvis merged commit 2ae505e into main Apr 3, 2026
9 checks passed
@rmjarvis rmjarvis deleted the roman3 branch April 3, 2026 16:58
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants