Merge branch 'main' into saturation-dilution

CHANGES.rst

@@ -10,6 +10,10 @@ general
 
 - Fix regression tests for PSF fitting methods [#872]
 
+- Fix regression test ``compare_asdf`` function replacing use of
+  ``asdf.commands.diff`` with ``deepdiff`` and add ``deepdiff`` as
+  a test dependency [#868]
+
 ramp_fitting
 ------------
 
@@ -57,6 +61,8 @@ ramp_fitting
 ------------
 - Update unit tests for stcal 1.4.0 [#725]
 
+- Adjust ramp slopes and associated unceratinties for gain. [#804]
+
 refpix
 ------
 

pyproject.toml

@@ -53,6 +53,7 @@ test = [
     'pytest >=4.6.0',
     'pytest-astropy',
     'metrics_logger >= 0.1.0',
+    'deepdiff',
 ]
 dev = [
     'romancal[docs,test]',

romancal/ramp_fitting/ramp_fit_step.py

@@ -103,14 +103,41 @@ def ols(self, input_model, readnoise_model, gain_model):
             max_cores,
             dqflags.pixel,
         )
+
+        if image_info is not None:
+            # JWST has a separate GainScale step.  This isn't in Roman.
+            # We can adjust the gains here instead.
+            # data, dq, var_poisson, var_rnoise, err = image_info
+            # This step isn't needed for ols_cas22, which works directly in
+            # electrons.
+            image_info[0][...] *= gain_model.data.value  # data
+            # no dq multiplication!
+            image_info[2][...] *= gain_model.data.value**2  # var_poisson
+            image_info[3][...] *= gain_model.data.value**2  # var_rnoise
+            image_info[4][...] *= gain_model.data.value  # err
+
         readnoise_model.close()
-        gain_model.close()
 
         # Save the OLS optional fit product, if it exists
         if opt_info is not None:
+            # We should scale these by the gain, too.
+            # opt_info = (slope, sigslope, var_poisson, var_readnoise, yint,
+            # sig_yint,
+            # ped_int, weights, cr_mag_seg
+            opt_info[0][...] *= gain_model.data.value  # slope
+            opt_info[1][...] *= gain_model.data.value  # sigslope
+            opt_info[2][...] *= gain_model.data.value**2  # var_poisson
+            opt_info[3][...] *= gain_model.data.value**2  # var_readnoise
+            opt_info[4][...] *= gain_model.data.value  # yint
+            opt_info[5][...] *= gain_model.data.value  # sigyint
+            opt_info[6][...] *= gain_model.data.value  # pedestal
+            # no change to weights due to gain
+            opt_info[8][...] *= gain_model.data.value  # crmag
             opt_model = create_optional_results_model(input_model, opt_info)
             self.save_model(opt_model, "fitopt", output_file=self.opt_name)
 
+        gain_model.close()
+
         # All pixels saturated, therefore returning an image file with zero data
         if image_info is None:
             log.info("All pixels are saturated. Returning a zeroed-out image.")

romancal/ramp_fitting/tests/test_ramp_fit_ols.py

@@ -41,22 +41,22 @@ def test_ols_multicore_ramp_fit_match(make_data):
     """Test various core amount calculation"""
     model, override_gain, override_readnoise = make_data
 
-    # gain or read noise are also modified in place in an important way (!)
+    # read noise is also modified in place in an important way (!)
     # so we make copies here so that we can get agreement.
     out_model = RampFitStep.call(
         model.copy(),  # model1 is modified in place now.
         algorithm="ols",
         maximum_cores="none",
         override_gain=override_gain,
-        override_readnoise=override_readnoise,
+        override_readnoise=override_readnoise.copy(),
     )
 
     all_out_model = RampFitStep.call(
         model.copy(),  # model1 is modified in place now.
         algorithm="ols",
         maximum_cores="all",
         override_gain=override_gain,
-        override_readnoise=override_readnoise,
+        override_readnoise=override_readnoise.copy(),
     )
 
     # Original ramp parameters
@@ -112,6 +112,7 @@ def test_ols_one_group_small_buffer_fit(max_cores, make_data):
     model, override_gain, override_readnoise = make_data
 
     model.data[0, 15, 10] = 10.0 * model.data.unit  # add single CR
+    override_gain.data[15, 10] = 2 * override_gain.data.unit
 
     out_model = RampFitStep.call(
         model,
@@ -124,7 +125,7 @@ def test_ols_one_group_small_buffer_fit(max_cores, make_data):
     data = out_model.data.value
 
     # Index changes due to trimming of reference pixels
-    np.testing.assert_allclose(data[11, 6], -1.0e-5, 1e-5)
+    np.testing.assert_allclose(data[11, 6], 10 * 2, 1e-5)
 
 
 @pytest.mark.parametrize("max_cores", MAXIMUM_CORES)
@@ -241,6 +242,7 @@ def generate_ramp_model(shape, deltatime=1):
     dm_ramp.groupdq = gdq
     dm_ramp.err = u.Quantity(err, u.DN, dtype=np.float32)
 
+    dm_ramp.meta.exposure.group_time = deltatime
     dm_ramp.meta.exposure.frame_time = deltatime
     dm_ramp.meta.exposure.ngroups = shape[0]
     dm_ramp.meta.exposure.nframes = 1

romancal/regtest/conftest.py

@@ -203,10 +203,16 @@ def rtdata_module(artifactory_repos, envopt, request, jail):
 @pytest.fixture
 def ignore_asdf_paths():
     ignore_attr = [
-        "meta.[date, filename]",
+        # generic asdf stuff that will contain program version numbers
+        # and other things that will almost certainly change in every case
         "asdf_library",
         "history",
-        "cal_logs",
+        # roman-specific stuff to ignore
+        "roman.cal_logs",
+        "roman.meta.date",
+        # roman.meta.filename is used by the ExposurePipeline so should likely
+        # not be ignored here
+        # "roman.meta.filename",
     ]
 
     return {"ignore": ignore_attr}

romancal/regtest/regtestdata.py

@@ -5,18 +5,23 @@
 import sys
 from difflib import unified_diff
 from glob import glob as _sys_glob
-from io import StringIO
 from pathlib import Path
 
 import asdf
+import astropy.time
+import deepdiff
+import gwcs
+import numpy as np
 import requests
-from asdf.commands.diff import diff as asdf_diff
+from astropy.units import Quantity
 from ci_watson.artifactory_helpers import (
     BigdataError,
     check_url,
     get_bigdata,
     get_bigdata_root,
 )
+from deepdiff.operator import BaseOperator
+from gwcs.wcstools import grid_from_bounding_box
 
 # from romancal.lib.suffix import replace_suffix
 from romancal.stpipe import RomanStep
@@ -525,8 +530,176 @@ def _data_glob_url(*url_parts, root=None):
     return url_paths
 
 
-def compare_asdf(result, truth, **kwargs):
-    f = StringIO()
-    asdf_diff([result, truth], minimal=False, iostream=f, **kwargs)
-    if f.getvalue():
-        f.getvalue()
+class NDArrayTypeOperator(BaseOperator):
+    def __init__(self, rtol=1e-05, atol=1e-08, equal_nan=True, **kwargs):
+        super().__init__(**kwargs)
+        self.rtol = rtol
+        self.atol = atol
+        self.equal_nan = equal_nan
+
+    def give_up_diffing(self, level, diff_instance):
+        a, b = level.t1, level.t2
+        meta = {}
+        if a.shape != b.shape:
+            meta["shapes"] = [a.shape, b.shape]
+        if a.dtype != b.dtype:
+            meta["dtypes"] = [a.dtype, b.dtype]
+        if isinstance(a, Quantity) and isinstance(b, Quantity):
+            if a.unit != b.unit:
+                meta["units"] = [a.unit, b.unit]
+        if not meta:  # only compare if shapes and dtypes match
+            if not np.allclose(
+                a, b, rtol=self.rtol, atol=self.atol, equal_nan=self.equal_nan
+            ):
+                abs_diff = np.abs(a - b)
+                index = np.unravel_index(np.nanargmax(abs_diff), a.shape)
+                meta["worst_abs_diff"] = {
+                    "index": index,
+                    "value": abs_diff[index],
+                }
+                with np.errstate(invalid="ignore", divide="ignore"):
+                    # 0 / 0 == nan and produces an 'invalid' error
+                    # 1 / 0 == inf and produces a 'divide' error
+                    # ignore these here for computing the fractional diff
+                    fractional_diff = np.abs(a / b)
+                index = np.unravel_index(np.nanargmax(fractional_diff), a.shape)
+                meta["worst_fractional_diff"] = {
+                    "index": index,
+                    "value": fractional_diff[index],
+                }
+                meta["abs_diff"] = np.nansum(np.abs(a - b))
+                meta["n_diffs"] = np.count_nonzero(
+                    np.isclose(
+                        a, b, rtol=self.rtol, atol=self.atol, equal_nan=self.equal_nan
+                    )
+                )
+        if meta:
+            diff_instance.custom_report_result("arrays_differ", level, meta)
+        return True
+
+
+class TimeOperator(BaseOperator):
+    def give_up_diffing(self, level, diff_instance):
+        if level.t1 != level.t2:
+            diff_instance.custom_report_result(
+                "times_differ",
+                level,
+                {
+                    "difference": level.t1 - level.t2,
+                },
+            )
+        return True
+
+
+def _wcs_to_ra_dec(wcs):
+    x, y = grid_from_bounding_box(wcs.bounding_box)
+    return wcs(x, y)
+
+
+class WCSOperator(BaseOperator):
+    def give_up_diffing(self, level, diff_instance):
+        # for comparing wcs instances this function evaluates
+        # each wcs and compares the resulting ra and dec outputs
+        # TODO should we compare the bounding boxes?
+        ra_a, dec_a = _wcs_to_ra_dec(level.t1)
+        ra_b, dec_b = _wcs_to_ra_dec(level.t2)
+        meta = {}
+        for name, a, b in [("ra", ra_a, ra_b), ("dec", dec_a, dec_b)]:
+            # TODO do we want to do something fancier than allclose?
+            if not np.allclose(a, b):
+                meta[name] = {
+                    "abs_diff": np.abs(a - b),
+                }
+        if meta:
+            diff_instance.custom_report_result(
+                "wcs_differ",
+                level,
+                meta,
+            )
+        return True
+
+
+class DiffResult:
+    def __init__(self, diff):
+        self.diff = diff
+
+    @property
+    def identical(self):
+        return not self.diff
+
+    def report(self, **kwargs):
+        return pprint.pformat(self.diff)
+
+
+def compare_asdf(result, truth, ignore=None, rtol=1e-05, atol=1e-08, equal_nan=True):
+    """
+    Compare 2 asdf files: result and truth. Note that this comparison is
+    asymmetric (swapping result and truth will give a different result).
+
+    Parameters
+    ----------
+
+    result : str
+        Filename of result asdf file
+
+    truth : str
+        Filename of truth asdf file
+
+    ignore : list
+        List of tree node paths to ignore during the comparison
+
+    rtol : float
+        rtol argument passed to `numpyp.allclose`
+
+    atol : float
+        atol argument passed to `numpyp.allclose`
+
+    equal_nan : bool
+        Ignore nan inequality
+
+    Returns
+    -------
+
+    diff_result : DiffResult
+        result of the comparison
+    """
+    exclude_paths = []
+    ignore = ignore or []
+    for path in ignore:
+        key_path = "".join([f"['{k}']" for k in path.split(".")])
+        exclude_paths.append(f"root{key_path}")
+    operators = [
+        NDArrayTypeOperator(
+            rtol, atol, equal_nan, types=[asdf.tags.core.NDArrayType, np.ndarray]
+        ),
+        TimeOperator(types=[astropy.time.Time]),
+        WCSOperator(types=[gwcs.WCS]),
+    ]
+    # warnings can be seen in regtest runs which indicate
+    # that ddtrace logs are evaluated at times after the below
+    # with statement exits resulting in access attempts on the
+    # closed asdf file. To try and avoid that we disable
+    # lazy loading and memmory mapping
+    open_kwargs = {
+        "lazy_load": False,
+        "copy_arrays": True,
+    }
+    with asdf.open(result, **open_kwargs) as af0, asdf.open(
+        truth, **open_kwargs
+    ) as af1:
+        # swap the inputs here so DeepDiff(truth, result)
+        # this will create output with 'new_value' referring to
+        # the value in the result and 'old_value' referring to the truth
+        diff = deepdiff.DeepDiff(
+            af1.tree,
+            af0.tree,
+            ignore_nan_inequality=equal_nan,
+            custom_operators=operators,
+            exclude_paths=exclude_paths,
+        )
+        # the conversion between NDArrayType and ndarray adds a bunch
+        # of type changes, ignore these for now.
+        # TODO Ideally we could find a way to remove just the NDArrayType ones
+        if "type_changes" in diff:
+            del diff["type_changes"]
+        return DiffResult(diff)

romancal/regtest/test_dark_current.py

@@ -20,7 +20,8 @@ def test_dark_current_subtraction_step(rtdata, ignore_asdf_paths):
     output = "r0000101001001001001_01101_0001_WFI01_darkcurrent.asdf"
     rtdata.output = output
     rtdata.get_truth(f"truth/WFI/image/{output}")
-    assert compare_asdf(rtdata.output, rtdata.truth, **ignore_asdf_paths) is None
+    diff = compare_asdf(rtdata.output, rtdata.truth, **ignore_asdf_paths)
+    assert diff.identical, diff.report()
 
 
 @pytest.mark.bigdata
@@ -40,7 +41,8 @@ def test_dark_current_outfile_step(rtdata, ignore_asdf_paths):
     output = "Test_darkcurrent.asdf"
     rtdata.output = output
     rtdata.get_truth(f"truth/WFI/image/{output}")
-    assert compare_asdf(rtdata.output, rtdata.truth, **ignore_asdf_paths) is None
+    diff = compare_asdf(rtdata.output, rtdata.truth, **ignore_asdf_paths)
+    assert diff.identical, diff.report()
 
 
 @pytest.mark.bigdata
@@ -61,7 +63,8 @@ def test_dark_current_outfile_suffix(rtdata, ignore_asdf_paths):
     output = "Test_darkcurrent.asdf"
     rtdata.output = output
     rtdata.get_truth(f"truth/WFI/image/{output}")
-    assert compare_asdf(rtdata.output, rtdata.truth, **ignore_asdf_paths) is None
+    diff = compare_asdf(rtdata.output, rtdata.truth, **ignore_asdf_paths)
+    assert diff.identical, diff.report()
 
 
 @pytest.mark.bigdata
@@ -83,4 +86,5 @@ def test_dark_current_output(rtdata, ignore_asdf_paths):
     output = "r0000101001001001001_01101_0001_WFI01_darkcurrent.asdf"
     rtdata.output = output
     rtdata.get_truth(f"truth/WFI/image/{output}")
-    assert compare_asdf(rtdata.output, rtdata.truth, **ignore_asdf_paths) is None
+    diff = compare_asdf(rtdata.output, rtdata.truth, **ignore_asdf_paths)
+    assert diff.identical, diff.report()

romancal/regtest/test_jump_det.py

@@ -28,4 +28,5 @@ def test_jump_detection_step(rtdata, ignore_asdf_paths):
     output = "r0000101001001001001_01101_0001_WFI01_jump.asdf"
     rtdata.output = output
     rtdata.get_truth(f"truth/WFI/image/{output}")
-    assert compare_asdf(rtdata.output, rtdata.truth, **ignore_asdf_paths) is None
+    diff = compare_asdf(rtdata.output, rtdata.truth, **ignore_asdf_paths)
+    assert diff.identical, diff.report()