RCAL-511 Implement Uneven ramp fitting #175

stscieisenhamer · 2023-06-22T12:23:57Z

This PR implements the Casertino, et.al., 2021 algorithm for fitting unevenly spaced ramps. Initial implementation based primarily on code from romanisim.

Checklist

added entry in CHANGES.rst (either in Bug Fixes or Changes to API)
updated relevant tests
updated relevant documentation
updated relevant milestone(s)
added relevant label(s)

codecov · 2023-06-22T12:26:10Z

Codecov Report

Patch coverage: 62.82% and project coverage change: -0.30% ⚠️

Comparison is base (b5bd209) 74.41% compared to head (2004ac5) 74.11%.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #175      +/-   ##
==========================================
- Coverage   74.41%   74.11%   -0.30%     
==========================================
  Files          29       33       +4     
  Lines        5944     6100     +156     
==========================================
+ Hits         4423     4521      +98     
- Misses       1521     1579      +58

Files Changed	Coverage Δ
setup.py	`0.00% <0.00%> (ø)`
src/stcal/ramp_fitting/ols_cas22_fit.py	`40.47% <40.47%> (ø)`
src/stcal/ramp_fitting/ols_cas22_util.py	`100.00% <100.00%> (ø)`
tests/test_ramp_fitting_cas22.py	`100.00% <100.00%> (ø)`

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

stscieisenhamer · 2023-06-26T19:21:20Z

Developer note: None of the model and reference file handling has yet to be added. Currently planned on handling this in romancal. Part of the reason is to get the basic code out for further development of the jump algorithm, which is closely related

stscieisenhamer · 2023-06-27T12:07:14Z

Just realized there was a commit that was neglected. Will be re-applying it shortly.

kmacdonald-stsci

Many functions are too long and complicated.

If reasonable, attempts should be made to keep variable names across the code base the same. For example, nreads and variations were often used interchangeably with group and frame time, making the code confusing to follow. Also, not sure why use the variable name resultant, rather than group, since they seem to be used for the same thing.

The code appears to properly implement the algorithm from Casertano's paper.

kmacdonald-stsci · 2023-06-27T17:35:48Z

src/stcal/ramp_fitting/ols_cas21.pyx

+        for j in range(resstart[i], resend[i] + 1):
+            # Casertano+22 Eq. 37
+            # Note that we are replacing ww with kk to save memory; we don't
+            # need ww again.


Should this say "we don't need kk again"?

kmacdonald-stsci · 2023-06-27T17:39:01Z

src/stcal/ramp_fitting/ols_cas21.pyx

+        The first resultant in this ramp
+    resend : np.ndarray[nramp]
+        The last resultant in this ramp.
+    """


This function is 100 lines long and does many things. In general functions should be able to completely fit on a screen, so roughly 35-45 lines, and shorter if possible. It would be easier to read and maintain if broken down into smaller functions.

kmacdonald-stsci · 2023-06-27T18:42:40Z

src/stcal/ramp_fitting/ols_cas21_util.py

+    """
+    firstreads = np.array([x[0] for x in ma_table])
+    nreads = np.array([x[1] for x in ma_table])
+    meantimes = read_time * firstreads + read_time * (nreads - 1) / 2


Is read_time frame time or group time? In other ramp fitting software, read_time or something similar was used interchangeably with group and frame time, causing confusion in how the code works.

kmacdonald-stsci · 2023-06-27T18:43:43Z

src/stcal/ramp_fitting/ols_cas21.pyx

+
+    Parameters
+    ----------
+    resultants : np.ndarry[nresultants, npixel]


Is there a reason to use the term "resultant", rather than "group"? It would be good to try to keep nomenclature across files the same, when it makes sense.

Roman and JWST are using different terminology for this. Roman terms this "resultant", so its not clear how to name this.

I understand that, but this is STCAL. It's supposed to be telescope agnostic. Developing parallel code bases with different nomenclature for the same thing undermines being agnostic.

kmacdonald-stsci · 2023-06-27T18:47:53Z

tests/test_ramp_fitting_cas21.py

+    assert result == expected
+
+
+def test_ramp(test_table=None):


There is a lot of things being tested in this one test. It should be broken into smaller tests, with each test testing fewer things.

kmacdonald-stsci · 2023-06-27T18:55:23Z

src/stcal/ramp_fitting/ols_cas21_fit.py

+
+    Returns
+    -------
+    par : np.ndarray[..., 2] (float)


Why use the variable name par?

kmacdonald-stsci · 2023-06-27T18:58:31Z

src/stcal/ramp_fitting/ols_cas21_fit.py

+        the read noise, Poisson source noise, and total noise.
+    """
+
+    # Get the Multi-accum table, either as given or from the read pattern


This is a very long function, over 150 lines and should be broken down into smaller functions.

schlafly · 2023-06-27T19:45:31Z

FWIW, some of the nomenclature answers:

read_time: I think this is frame time in JWST. It's one read of the array and is ~3.04 s for Roman. The reads/frames get averaged into resultants/groups.
resultant vs. group: yes, these are the same things. No, I don't know why Roman adopted the "resultant" language and JWST the "group" language, but "Roman" is pretty insistent on using resultants as names in place of groups. I don't feel strongly about global search and replace resultant for group in this codebase, but more broadly we will just have to live with the fact that the two missions went with different words for the same quantities.
par and var were short for parameters and variances, but these could be given other names.

ddavis-stsci

Modulo the given comments, LGTM.
Since I don't see JWST using this for their ramp analysis so the Roman specific nomenclature should not be too confusing.

stscieisenhamer · 2023-07-13T14:50:56Z

Finally addressing comments. @kmacdonald-stsci : Agree with most everything.

Driver of not doing too much refactoring/naming conventions are primarily that this algorithm is still undergoing algorithmic/science discussion. Kept the implementation close to the original so that discussion would not be confused by refactor/rename issues. Also, the JWST use is still down-the-road so the terminology issue, though definitely needing to be addressed, will be so later on.

Of course, this is all modulo what RCAL would like to see in the current PR.

schlafly · 2023-07-13T21:01:05Z

Hi Jonathan,

This looks good. Thinking about next steps:

I would go ahead and get rid of the utility code supporting the ramp fit interpolator that has been removed. I'll keep it in the simulator if we ever need it. I think removing it would mean dropping
- construct_covar
- construct_ramp_fitting_matrices
- construct_ki_and_variances
- ki_and_variance_grid
- resultants_to_differences
- associated tests
Let's either pull in simulate_many_ramps for test_simulated_ramps or drop that test. I think I'm in favor of pulling it in; that is a pretty good test.
I'm happy to sign on to s/resultants/r/groups to improve terminology consistency, etc..

stscieisenhamer · 2023-07-20T18:14:38Z

Some external feedback on the terminology issue: At the 20230720 TIPS, the terminology issue was brought up, of which Casertano defended the current uneven/Roman centric form. Unless there is strong reason to change in this PR, I propose to deal with it in a later PR.

WilliamJamieson · 2023-07-20T18:29:07Z

Some external feedback on the terminology issue: At the 20230720 TIPS, the terminology issue was brought up, of which Casertano defended the current uneven/Roman centric form. Unless there is strong reason to change in this PR, I propose to deal with it in a later PR.

IMOP I think terminology should follow whatever the reference paper uses if possible as we may have to refer back to the paper in the future. It is unfortunate that the two telescopes could not agree on terminology.

stscieisenhamer · 2023-07-20T18:47:19Z

@schlafly The simulate_many_ramps has further romanisim and galsim dependencies. The current solution allows full testing if the relevant packages are present. If there is strong desire to have this test execute on CI, romanisim can be a test dependency. Otherwise, galsim.GaussianDeviate(seed) will need to be dealt with, plus copy of at least two other functions.

schlafly · 2023-07-20T20:21:35Z

Oh no, you're right, I don't want to bring in l1.apportion_counts_to_resultants. That's more detailed than we need for this. Let me rig up a simulate_many_ramps replacement that would be appropriate here; it's not so bad.

schlafly · 2023-07-20T20:31:22Z

src/stcal/ramp_fitting/ols_cas22_util.py

+#   For Roman, the read time of the detectors is a fixed value and is currently
+#   backed into code. Will need to refactor to consider the more general case.
+#   Used to deconstruct the MultiAccum tables into integration times.
+READ_TIME = 3.04


Oops, this is Roman specific. We need to pull this out as an argument.

Seems there was some recent discussion of this value, potentially coming from PPS. I am not seeing anything like this in the schema. What is the process to retrieve from PPS and add as a meta value to the Level 1 file?

also note: this is already an argument. the above is the default. However, I will move this into romancal and make the argument required.

This is the value to use: https://github.com/spacetelescope/rad/blob/main/src/rad/resources/schemas/exposure-1.0.0.yaml#L268
I don't know if we're populating it yet usually.

Sorry I missed that you were just using FRAME_TIME as a default. I agree with you that it makes more sense to have that on the romancal side.

schlafly · 2023-07-20T20:37:55Z

With the caveat that I haven't tested anything, so I'm sure it's full of bugs, here's a simpler version of simulate_many_ramps...

def simulate_many_ramps(ntrial=100, flux=100, readnoise=5, ma_table=None):
    if ma_table is None:
        ma_table = [[1, 4], [5, 1], [6, 3], [9, 10], [19, 3], [22, 15]]
    nread = np.array([x[1] for x in ma_table])
    tij = ma_table_to_tij(ma_table)
    resultants = np.zeros((len(ma_table), ntrial), dtype='f4')
    buf = np.zeros(ntrial, dtype='i4')
    for i, ti in enumerate(tij):
        for t0 in ti:
            buf += np.random.poisson(READ_TIME * flux, ntrial)
        resultants[i] = (buf / len(ti)).astype('f4')
    resultants += np.random.randn(len(ma_table), ntrial) * (
        readnoise / np.sqrt(nread)).reshape(len(ma_table), 1)
    return (ma_table, flux, readnoise, resultants)

This is a good simulation of a ramp in the context of the assumptions inherent in the algorithm. That's better than trying to bring over stuff from l1.apportion_counts_to_resultants, sorry for the bad pointer.

Changes: - fix circular import between `matable_fit` and `matable_fit_cas2022`. - Update tests from `romanisim` to work completely within `stcal`. - Allow test to use `romanisim` if present.

Initial implementation done. However, untested and tests need to be implemented.

This reverts commit a34fc59.

This reverts commit 1d310a2.

This reverts commit 63512ec.

stscieisenhamer · 2023-08-03T02:38:12Z

Got the romanisim test dependency completed. However, I want to point out a C-ism that was causing issues. In
ols_cas22.pyx:85, which is

fabs((tbar[start + i] - tbarmid) / tscale) ** weight_power)

generates an error concerning not being able to coerce a complex number to double. This error appeared after taking in the memory enhancement changes Eddy made in the romanisim version. The condition when this occurs seems to be when tbar[start + i] == tbarmid. I can only guess that the floating math is creating a -0.0 situation. Note that this is not reproducible with pure Python.

The error does suggest setting of a cython parameter, cython.cpow(True), to mitigate the situation. Doing so, as best as I can tell, does remediate the exception, and produces the expected results.

schlafly · 2023-08-03T02:46:18Z

Oops, I hit this one too. This happened when cython 3 came out a week or two ago. I tried pinging you here but should have followed up. spacetelescope/romanisim#67 (comment)

stscieisenhamer requested a review from WilliamJamieson June 22, 2023 12:23

github-actions bot added installation ramp_fitting testing labels Jun 22, 2023

stscieisenhamer force-pushed the rcal-511-rampfit branch from f38cdbc to f27f0b4 Compare June 26, 2023 18:51

stscieisenhamer marked this pull request as ready for review June 26, 2023 19:16

stscieisenhamer requested a review from a team as a code owner June 26, 2023 19:16

stscieisenhamer requested review from schlafly and ddavis-stsci June 26, 2023 19:17

stscieisenhamer changed the title ~~RCAL-511 Implement Uneven ramp fitting~~ WIP: RCAL-511 Implement Uneven ramp fitting Jun 27, 2023

stscieisenhamer marked this pull request as draft June 27, 2023 12:06

kmacdonald-stsci reviewed Jun 27, 2023

View reviewed changes

ddavis-stsci approved these changes Jun 30, 2023

View reviewed changes

stscieisenhamer force-pushed the rcal-511-rampfit branch from 014cbc7 to 44a7109 Compare July 10, 2023 15:31

stscieisenhamer mentioned this pull request Jul 13, 2023

RCAL-511 Inititial implementation of the Uneven Ramp fitting spacetelescope/romancal#779

Merged

4 tasks

stscieisenhamer force-pushed the rcal-511-rampfit branch from 0832df9 to 553dd45 Compare July 13, 2023 14:41

stscieisenhamer marked this pull request as ready for review July 13, 2023 14:44

kmacdonald-stsci approved these changes Jul 17, 2023

View reviewed changes

schlafly reviewed Jul 20, 2023

View reviewed changes

Jonathan Eisenhamer and others added 22 commits August 2, 2023 22:24

initial ma ramp fitting tests implemented and passing

63d3a33

fix circular import

31ecb2d

Changes: - fix circular import between `matable_fit` and `matable_fit_cas2022`. - Update tests from `romanisim` to work completely within `stcal`. - Allow test to use `romanisim` if present.

rename algorithm from matable to ols_cas21

0bbbfc2

implement conversion from read pattern to multi-accum table

6a0d546

add conversion from multi-accum tables to read patterns

c6fd767

Update fitting routines to take both multi-accum and read pattern

60ea1ef

fixed doc test format

4e48a0e

add read pattern to RampData

deb7251

refactor to use match statements

033958a

WIP: Integrating CAS21 algorithm into the ramp_fit api

c90125c

Initial implementation done. However, untested and tests need to be implemented.

Revert "WIP: Integrating CAS21 algorithm into the ramp_fit api"

8c459d3

This reverts commit a34fc59.

Revert "refactor to use match statements"

6629a6e

This reverts commit 1d310a2.

Revert "add read pattern to RampData"

287866e

This reverts commit 63512ec.

fix (most) ruff issues and tox setup

3e3b5f2

remove stray import

4b0ab20

update changelog

71019bc

Rename references to the Casertano to year 2022

3f89365

update ramp fitting from romanisim memory efficiency work

5aacab1

Raise exception when read_pattern does not match resultant data

c585b2c

remove code unused by the ramp fitting code

13cd212

Parameterize the read time

88585af

Remove romanisim test dependency and protect against complex results

36a1136

stscieisenhamer force-pushed the rcal-511-rampfit branch from ef8b39f to 36a1136 Compare August 3, 2023 02:25

stscieisenhamer requested a review from schlafly August 3, 2023 02:38

simplify the variances return

10b8398

stscieisenhamer changed the title ~~WIP: RCAL-511 Implement Uneven ramp fitting~~ RCAL-511 Implement Uneven ramp fitting Aug 3, 2023

handle special case where there are no or only one resultant

2004ac5

stscieisenhamer merged commit 0c7cb32 into spacetelescope:main Aug 11, 2023
16 of 18 checks passed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

RCAL-511 Implement Uneven ramp fitting #175

RCAL-511 Implement Uneven ramp fitting #175

stscieisenhamer commented Jun 22, 2023 •

edited

Loading

codecov bot commented Jun 22, 2023 •

edited

Loading

stscieisenhamer commented Jun 26, 2023

stscieisenhamer commented Jun 27, 2023

kmacdonald-stsci left a comment

kmacdonald-stsci Jun 27, 2023

kmacdonald-stsci Jun 27, 2023

kmacdonald-stsci Jun 27, 2023

kmacdonald-stsci Jun 27, 2023

WilliamJamieson Jul 6, 2023

kmacdonald-stsci Jul 6, 2023

kmacdonald-stsci Jun 27, 2023

kmacdonald-stsci Jun 27, 2023

kmacdonald-stsci Jun 27, 2023

schlafly commented Jun 27, 2023

ddavis-stsci left a comment

stscieisenhamer commented Jul 13, 2023

schlafly commented Jul 13, 2023

stscieisenhamer commented Jul 20, 2023

WilliamJamieson commented Jul 20, 2023

stscieisenhamer commented Jul 20, 2023

schlafly commented Jul 20, 2023

schlafly Jul 20, 2023

stscieisenhamer Jul 21, 2023

stscieisenhamer Jul 21, 2023

schlafly Jul 21, 2023

schlafly Jul 21, 2023

schlafly commented Jul 20, 2023

stscieisenhamer commented Aug 3, 2023

schlafly commented Aug 3, 2023

RCAL-511 Implement Uneven ramp fitting #175

RCAL-511 Implement Uneven ramp fitting #175

Conversation

stscieisenhamer commented Jun 22, 2023 • edited Loading

codecov bot commented Jun 22, 2023 • edited Loading

Codecov Report

stscieisenhamer commented Jun 26, 2023

stscieisenhamer commented Jun 27, 2023

kmacdonald-stsci left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

schlafly commented Jun 27, 2023

ddavis-stsci left a comment

Choose a reason for hiding this comment

stscieisenhamer commented Jul 13, 2023

schlafly commented Jul 13, 2023

stscieisenhamer commented Jul 20, 2023

WilliamJamieson commented Jul 20, 2023

stscieisenhamer commented Jul 20, 2023

schlafly commented Jul 20, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

schlafly commented Jul 20, 2023

stscieisenhamer commented Aug 3, 2023

schlafly commented Aug 3, 2023

stscieisenhamer commented Jun 22, 2023 •

edited

Loading

codecov bot commented Jun 22, 2023 •

edited

Loading