Skip to content

Initial implementation of Roman PSF#195

Merged
rmjarvis merged 21 commits intomainfrom
roman1
Mar 25, 2026
Merged

Initial implementation of Roman PSF#195
rmjarvis merged 21 commits intomainfrom
roman1

Conversation

@rmjarvis
Copy link
Copy Markdown
Owner

I think this set of commits is pretty much the minimum set of changes to implement the Roman PSF type. It's not quite everything that I have in the roman branch, but the main swaths of the code are here.

Things still to come in later PRs (to help ease the code review burden):

  1. aberration_interp="linear" option
  2. Running the SCAs in multiple processes
  3. chromatic=True (This technically has it, but it is too inefficient to be useful.)
  4. Misc changes I made to diagnostics and other parts of the code along the way.

@arunkannawadi arunkannawadi self-requested a review March 19, 2026 23:04
@nihardalal nihardalal self-assigned this Mar 20, 2026
@nihardalal nihardalal self-requested a review March 20, 2026 00:31
Copy link
Copy Markdown
Collaborator

@arunkannawadi arunkannawadi left a comment

Choose a reason for hiding this comment

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

I've taken only a quick pass for now, but wanted to leave the comments I have so far.

Comment thread piff/__init__.py Outdated

# Optics
from .optical_model import Optical, optical_templates
from .roman_psf import Roman, RomanOptics
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.

Unlike other items, this seems very Roman-specific. Do we want this in our namespace by default? Could we do something like you have for des below?

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.

Sure. That's a good idea. We can have a roman directory with this file and any others we decide to add relevant to the Roman PSF.

Comment thread piff/roman_psf.py Outdated
)


class Roman(Model):
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.

I find the naming unintuitive. If I saw a bare Roman class in code, I wouldn't know what that's actually referring to. RomanPSF would be better given the docstring.

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.

Agreed. These aren't great names. How about RomanOpticalModel and RomanOpticalPSF.

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.

I think I prefer RomanOptics over RomanOptical for the written type in the config file, so I think RomanOpticsPSF for the class name works better. But RomanOpticalModel is still good for the model class name.

Comment thread piff/roman_psf.py Outdated
)


class RomanOptics(PSF):
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.

Given that this is subclassing PSF, I wonder if this should be RomanOpticalPSF or something like that. I see other examples of the parent class' name appearing, like RomanSCAInterp when it inherits from Interp.

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.

A few small changes here and there, feel free to resolve as necessary. Also, in running the pytest for this branch, I also got an error in test_input.py in test_stars() where assert len(stars) == 90 fails on line 1683 under the "gratuitous coverage test"

I was running the tests with the wrong environment - all tests passed on my end with the correct environment, sorry for the confusion

Comment thread devel/roman/piff.yaml Outdated
select:

min_snr: 50 # reject very faint stars
#max_snr: 500 # don't over-weight very bright stars
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.

Could change max_snr to max_snr_weight - I think this gets the main point of the parameter across better

Suggested change
#max_snr: 500 # don't over-weight very bright stars
#max_snr_weight: 500 # don't over-weight very bright stars

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.

Good idea. I like this name. I'll keep the other one as deprecated, but I agree this is a better name.

Comment thread devel/roman/piff.yaml Outdated
Comment thread piff/roman/roman_psf.py
mean = np.array([self.sca_mean[s] for s in sca])
data = np.array(
list(zip(sca, mean)),
dtype=[('sca', int), ('mean', float, (len(self.global_mean),))],
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 line confused me at first - might be good to document to say "writes an n_param length array for each sca"

Comment thread piff/roman/roman_psf.py
each SCA. For each ``(sca, params)`` state we build corner PSFs at the four SCA corners and
interpolate between them at each star position. This keeps the optical model much faster when
many stars on a chip share similar parameter vectors.
"""
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.

Seems like a reasonable strategy. I do know that Cycle9 Zernikes are actually measured at 5 points across the chip - 4 corners and the center. Does it make sense to use all 5 points and do some kind of grid interpolation? (I guess the improvement is marginal, but its nice to have a one to one between the cycle9 data and the drawing in this case)

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.

OK. I could think about how to interpolate that. 5 is a bit awkward for interpolation. I guess you get a single quadratic term? 4 gives you (1,x,y,xy). 5 gives either (1,x,y,x^2,y^2) or (1,x,y,xy,x^2+y^2) I guess. Or we could do linear in each triangle?

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.

Here's a function that we use in psfsim to interpolate the Zernikes. I think scipy.interpolate.griddata can also work

def altgriddata(points, values, xi):
    """
    Quadratic fitting function.

    This is a substitute for ``scipy.interpolate.griddata`` when we have 5 points
    and want to fit a quadratic, minus the x^2-y^2 term.

    Paramters
    ---------
    points : np.ndarray of float
        The coordinates of the 5 points where we have values; shape (5, 2).
    values : np.ndarray of flat
        The values to interpolate; shape (5,).
    xi : (float, float)
        The location to interpolate to.

    Returns
    -------
    float
        The interpolated value.

    """

    # Set up the linear system.
    # Order of coefficients is 1, x, y, x*y, x^2+y^2.
    bval = np.zeros((5, 5))  # bval[i,j] is the ith basis function at jth point

    bval[0, :] = 1.0
    bval[1, :] = points[:, 0]
    bval[2, :] = points[:, 1]
    bval[3, :] = points[:, 0] * points[:, 1]
    bval[4, :] = points[:, 0] ** 2 + points[:, 1] ** 2

    mat = bval @ bval.T
    vec = bval @ values

    coefs = np.linalg.solve(mat, vec)

    return (
        coefs[0]
        + coefs[1] * xi[0]
        + coefs[2] * xi[1]
        + coefs[3] * xi[0] * xi[1]
        + coefs[4] * (xi[0] ** 2 + xi[1] ** 2)
    )

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.

Right. Thanks. So (1,x,y,xy,x^2+y^2) then. I can implement that. I'll do it in a different PR though.

Comment thread piff/roman_psf.py Outdated
Comment thread tests/test_roman.py
Comment thread tests/test_roman.py
Comment thread tests/test_roman.py
Comment thread tests/test_roman.py
Comment thread tests/test_roman.py
assert isinstance(psf2.interp, RomanSCAInterp)
assert psf2.interp.global_mean is None
assert psf2.interp.sca_mean == {}

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.

I think these tests look okay. I haven't been able to check for comprehensive coverage given the size of the file, but I think something like codecov would do a better job than me.

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.

I do use codecov. Looks like I missed one line of coverage when I pulled out this subset. I think the roman branch has 100% coverage. But I guess I missed one test line when I was assembling this subset.

cf. https://github.com/rmjarvis/Piff/pull/195/checks?check_run_id=67957433502

rmjarvis and others added 7 commits March 20, 2026 15:11
Co-authored-by: Nihar Dalal <140463192+nihardalal@users.noreply.github.com>
Co-authored-by: Nihar Dalal <140463192+nihardalal@users.noreply.github.com>
Co-authored-by: Nihar Dalal <140463192+nihardalal@users.noreply.github.com>
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.

I think everything looks okay now, so good to merge whenever

@rmjarvis
Copy link
Copy Markdown
Owner Author

I think I'm up to date with responding to suggestions.
@arunkannawadi Did you want to look at this some more before I merge?

@arunkannawadi
Copy link
Copy Markdown
Collaborator

I can take another look, but if I don't get back to you by EOD tomorrow, don't wait up for me.

@rmjarvis rmjarvis merged commit 6190aa3 into main Mar 25, 2026
9 checks passed
@rmjarvis rmjarvis deleted the roman1 branch March 25, 2026 21:03
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.

3 participants