Keep the floating dtype of the input in the host round trips

CHANGES.rst

@@ -80,9 +80,10 @@ New Features
   cleaned by ``cosmicray_lacosmic`` with a NumPy ``float64`` ``pssl`` (such
   as ``np.median`` of ``float64`` data) now stays ``float32``, as it already
   did for every other type of ``pssl``. Likewise, ``cosmicray_lacosmic`` with
-  ``gain_apply=True`` now returns astroscrappy's ``float32`` data for any
-  gain, as it already did for a gain of 1; previously any other gain
-  promoted the data to ``float64``. [#930, #933, #935, #1011]
+  ``gain_apply=True`` now returns ``float32`` data for integer and
+  ``float32`` input whatever the gain, as it already did for a gain of 1;
+  previously any other gain promoted the data to ``float64``.
+  [#930, #933, #935, #1011, #1023]
 
 Other Changes and Additions
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -131,6 +132,12 @@ Bug Fixes
 - ``wcs_project`` no longer returns a NumPy-backed ``CCDData`` for non-NumPy
   input: the reprojected data and mask are converted back to the array
   namespace and device of the input image. [#930, #1011]
+- ``wcs_project``, ``subtract_overscan`` with a ``model`` and
+  ``cosmicray_lacosmic`` now return floating data in its input dtype rather
+  than in the dtype of the NumPy-only library that does the work. [#1023]
+- ``cosmicray_lacosmic`` now accepts a ``float64`` ``invar``, which it
+  passes to astroscrappy as ``float32`` instead of letting astroscrappy
+  raise ``ValueError``. [#1032]
 - ``flat_correct`` now checks the flat's mask with ``xp.any`` instead of the
   ``mask.any()`` method, which arrays of spec-only array-API namespaces do
   not have. [#1005]

ccdproc/core.py

@@ -215,22 +215,40 @@ def _from_numpy(arr, like=None, *, xp=None):
     array or None
         ``arr`` as an array of ``xp`` on the device of ``like`` (on the
         default device of ``xp`` if ``like`` is `None`), or `None` if
-        ``arr`` is `None`. The result has the dtype of ``arr``, not of
-        ``like``, unless ``xp`` cannot represent it: JAX without 64-bit
-        mode silently converts float64 to float32.
+        ``arr`` is `None`. If ``arr`` and ``like`` both have a real
+        floating dtype, the result has the dtype of ``like``; otherwise,
+        including when there is no ``like``, it has the dtype of ``arr``,
+        unless ``xp`` cannot represent it: JAX without 64-bit mode
+        silently converts float64 to float32.
 
     Notes
     -----
     This is the inverse of `_to_numpy`: it returns the result of a
     NumPy-only operation to the caller's array namespace and device, so
     that the caller never receives a NumPy array in place of what it
     passed in.
+
+    The NumPy-only libraries decide the dtype of their own output:
+    ``reproject`` and ``astropy.modeling`` always return float64 and
+    astroscrappy always returns float32. Floating input keeps its dtype
+    instead (#1023), so a floating result is cast to the floating dtype of
+    ``like``. A boolean result, such as a mask, and the result for integer
+    input keep their own dtype: casting a floating result back to an
+    integer dtype would wrap negative values and could not hold the NaN
+    that ``wcs_project`` needs.
     """
     if arr is None:
         return None
     xp = xp or array_api_compat.array_namespace(like)
-    device = None if like is None else array_api_compat.device(like)
-    return xp.asarray(arr, device=device)
+    device = None
+    dtype = None
+    if like is not None:
+        device = array_api_compat.device(like)
+        if np.isdtype(arr.dtype, "real floating") and xp.isdtype(
+            like.dtype, "real floating"
+        ):
+            dtype = like.dtype
+    return xp.asarray(arr, dtype=dtype, device=device)
 
 
 def _ccddata_from_numpy(ccd, like=None, *, xp=None):
@@ -258,16 +276,24 @@ def _ccddata_from_numpy(ccd, like=None, *, xp=None):
     -----
     Arrays read from a FITS file may be big-endian, which only NumPy
     accepts, so each is put in native byte order with `_native_numpy`
-    before `_from_numpy` converts it.
+    before it is converted.
+
+    The conversion does not go through `_from_numpy`, whose floating dtype
+    rule is for results computed from ``like`` (#1023). A file is not such
+    a result, so its arrays keep the dtype they were read with: a
+    ``float64`` file combined with a ``float32`` image in memory must not
+    lose precision before ``combine`` casts it to its ``dtype``.
 
     The mask is left in NumPy: astropy's `~astropy.nddata.CCDData` mask
     setter converts any mask to NumPy, so a caller that needs the mask in
     ``xp`` moves it there itself.
     """
-    ccd.data = _from_numpy(_native_numpy(ccd.data), like=like, xp=xp)
+    xp = xp or array_api_compat.array_namespace(like)
+    device = None if like is None else array_api_compat.device(like)
+    ccd.data = xp.asarray(_native_numpy(ccd.data), device=device)
     if ccd.uncertainty is not None:
-        ccd.uncertainty.array = _from_numpy(
-            _native_numpy(ccd.uncertainty.array), like=like, xp=xp
+        ccd.uncertainty.array = xp.asarray(
+            _native_numpy(ccd.uncertainty.array), device=device
         )
     return ccd
 
@@ -1291,7 +1317,12 @@ def subtract_overscan(
         of = fitting.LinearLSQFitter()
         yarr = np.arange(oscan_np.shape[0])
         fitted = of(model, yarr, oscan_np)
-        oscan = _from_numpy(fitted(yarr), like=ccd.data, xp=xp)
+        # Convert like the reduced overscan rather than like ccd.data: it is
+        # always floating, so the fit comes back in the dtype the median and
+        # mean paths subtract, including for integer data, where that dtype
+        # depends on the array library (for example float64 for uint16 on
+        # NumPy, float32 on JAX).
+        oscan = _from_numpy(fitted(yarr), like=oscan, xp=xp)
         if overscan_axis == 1:
             oscan = xp.reshape(oscan, (oscan.size, 1))
         else:
@@ -1873,6 +1904,9 @@ def wcs_project(ccd, target_wcs, target_shape=None, order="bilinear"):
     data and mask are copied to the host, and the reprojected data and mask
     are copied back to the array namespace and device of the input; the
     copy is announced with a `HostCopyWarning`.
+
+    Floating data keeps its dtype, even though ``reproject`` works in
+    float64.
     """
     from astropy.nddata.ccddata import _generate_wcs_and_update_header
     from reproject import reproject_interp
@@ -2863,8 +2897,12 @@ def cosmicray_lacosmic(
 
         if gain_apply:
             if nccd.uncertainty is not None:
+                # Give the gain the uncertainty's dtype so that multiplying
+                # by it does not change that dtype.
                 gain_value = xp.asarray(
-                    float(gain.value), device=array_api_compat.device(_ccd.data)
+                    float(gain.value),
+                    dtype=xp.result_type(nccd.uncertainty.array.dtype, xp.float32),
+                    device=array_api_compat.device(_ccd.data),
                 )
                 gain_corrected = _ccd.multiply(
                     gain_value, xp=xp, handle_mask=xp.logical_or
@@ -2965,6 +3003,12 @@ def _lacosmic_on_host(
         for name, value in (("inbkg", inbkg), ("invar", invar))
         if value is not None
     }
+    # astroscrappy copies inbkg to float32 but not invar, and its cleaning
+    # routines only accept float32, so a float64 invar would make it raise.
+    if _is_array(invar):
+        background_kwargs["invar"] = background_kwargs["invar"].astype(
+            np.float32, copy=False
+        )
 
     # pssl is added on the host: an integer array plus a float pssl is
     # refused by some namespaces (array-api-strict), and adding it on the

ccdproc/tests/pytest_fixtures.py

@@ -35,6 +35,7 @@ def ccd_data(
     data_scale=DEFAULT_DATA_SCALE,
     data_mean=DEFAULT_DATA_MEAN,
     rng_seed=DEFAULT_SEED,
+    dtype=None,
 ):
     """
     Return a CCDData object with units of ADU.
@@ -51,6 +52,10 @@ def ccd_data(
 
     The mean can be changed with the marker @pytest.marker.scale(m) on the
     test function, where m is the desired mean.
+
+    If ``dtype``, the name of a dtype of the array library under test such
+    as ``"float32"`` or ``"uint16"``, is given, the random values are cast
+    to it.
     """
     # Need the import here to avoid circular import issues
     from ..conftest import testing_array_device as xp_device
@@ -65,8 +70,12 @@ def ccd_data(
 
     data = rng.normal(loc=mean, size=[size, size], scale=scale)
 
+    data = xp.asarray(data, device=xp_device)
+    if dtype is not None:
+        data = xp.astype(data, getattr(xp, dtype))
+
     fake_meta = {"my_key": 42, "your_key": "not 42"}
-    ccd = CCDData(xp.asarray(data, device=xp_device), unit=u.adu)
+    ccd = CCDData(data, unit=u.adu)
     ccd.header = fake_meta
     return ccd
 

ccdproc/tests/test_combiner.py

@@ -1803,6 +1803,32 @@ def test_combine_mixed_list_with_masked_file(tmp_path, file_first):
     assert not np.any(_to_numpy(result.mask))
 
 
+@pytest.mark.parametrize("file_first", [True, False], ids=["file-first", "ccd-first"])
+def test_combine_mixed_list_keeps_float64_file_precision(tmp_path, file_first):
+    """
+    A ``float64`` file combined with a ``float32`` image in memory keeps
+    its ``float64`` values until ``combine`` casts it to its ``dtype``.
+
+    Notes
+    -----
+    Guards the interaction of #1023 with #1025. ``combine`` reads files onto
+    the device of the image in memory with ``_ccddata_from_numpy``. If that
+    applied ``_from_numpy``'s floating dtype rule, meant for results
+    computed from an input, the file would first be cast to the ``float32``
+    of the image in memory, rounding ``1 + 2**-30`` to ``1``.
+    """
+    value = 1 + 2**-30
+    path = tmp_path / "float64.fits"
+    CCDData(np.full((3, 3), value), unit=u.adu).write(path)
+    in_memory = CCDData(to_xp(np.zeros((3, 3), dtype=np.float32)), unit=u.adu)
+    img_list = [str(path), in_memory] if file_first else [in_memory, str(path)]
+
+    result = combine(img_list)
+
+    assert result.data.dtype == xp.float64
+    np.testing.assert_array_equal(_to_numpy(result.data), np.full((3, 3), value / 2))
+
+
 def test_combine_files_keeps_uncertainty_dtype(tmp_path):
     """
     ``combine`` of files with a float32 uncertainty into ``array_package``

ccdproc/tests/test_cosmicray.py

@@ -404,6 +404,28 @@ def test_cosmicray_lacosmic_invar_inbkg(new_args):
         cosmicray_lacosmic(ccd_data, sigclip=5.9, **new_args)
 
 
+def test_cosmicray_lacosmic_float64_invar():
+    """
+    A float64 ``invar`` finds the cosmic rays instead of making astroscrappy
+    raise.
+
+    Notes
+    -----
+    A variance image built with NumPy defaults is float64. astroscrappy
+    works in float32 and does not convert ``invar`` itself, so it raised a
+    ``ValueError`` about the item size of the buffer.
+    """
+    ccd_data = ccd_data_func(data_scale=DATA_SCALE)
+    add_cosmicrays(ccd_data, DATA_SCALE, threshold=10, ncrays=NCRAYS)
+    # The data are near zero, so their variance is the read noise squared.
+    readnoise = 6.5
+    invar = xp.full(ccd_data.shape, readnoise**2, dtype=xp.float64, device=xp_device)
+
+    result = cosmicray_lacosmic(ccd_data, sigclip=5.9, readnoise=readnoise, invar=invar)
+
+    assert count_true(result.mask) == NCRAYS
+
+
 def test_cosmicray_median_check_data():
     with pytest.raises(TypeError):
         ndata, crarr = cosmicray_median(10, thresh=5, mbox=11, error_image=DATA_SCALE)
@@ -929,6 +951,95 @@ def no_cosmics(data, **_kwargs):
     assert_allclose(_to_numpy(cleaned), frame, rtol=1e-6)
 
 
+@pytest.mark.parametrize("array_input", [True, False])
+@pytest.mark.parametrize("dtype", ["float32", "float64"])
+def test_cosmicray_lacosmic_keeps_floating_dtype(dtype, array_input):
+    """
+    Floating data comes back from ``cosmicray_lacosmic`` in its own dtype,
+    for both a ``CCDData`` and a bare array (#1023).
+
+    Notes
+    -----
+    astroscrappy always returns float32, and its result used to be returned
+    as it was, so float64 data came back as float32.
+    """
+    ccd_data = ccd_data_func(data_scale=DATA_SCALE, data_mean=1000.0, dtype=dtype)
+    add_cosmicrays(ccd_data, DATA_SCALE, threshold=10, ncrays=NCRAYS)
+    ccd = ccd_data.data if array_input else ccd_data
+
+    result = cosmicray_lacosmic(ccd, sigclip=5.9, gain=2.0)
+
+    cleaned = result[0] if array_input else result.data
+    assert cleaned.dtype == getattr(xp, dtype)
+
+
+@pytest.mark.parametrize("gain", [1.0, 2.0])
+@pytest.mark.parametrize("dtype", ["uint16", "int32"])
+def test_cosmicray_lacosmic_integer_input_gives_float32(dtype, gain):
+    """
+    Integer data comes back from ``cosmicray_lacosmic`` as float32, whatever
+    the integer type and whatever the gain.
+
+    Notes
+    -----
+    float32 is the only dtype astroscrappy provides. This was settled in
+    #1023: floating input keeps its dtype, but integer input is not
+    promoted to float64 just to match other ccdproc functions. The integer
+    types are the ones NumPy promotes differently with float32 (uint16
+    gives float32, int32 gives float64), so a result left to type promotion
+    would differ between them. In 2.5.1 a gain other than 1 promoted the
+    result to float64 by accident, through a NumPy float64 scalar; the
+    ``gain`` parametrization pins that it no longer does.
+    """
+    ccd_data = ccd_data_func(data_scale=DATA_SCALE, data_mean=1000.0, dtype=dtype)
+    add_cosmicrays(ccd_data, DATA_SCALE, threshold=10, ncrays=NCRAYS)
+
+    result = cosmicray_lacosmic(ccd_data, sigclip=5.9, gain=gain)
+
+    assert result.data.dtype == xp.float32
+
+
+@pytest.mark.parametrize(
+    "dtype",
+    [
+        "float32",
+        "float64",
+        pytest.param(
+            "uint16",
+            marks=pytest.mark.backend_xfail(
+                "array-api-strict",
+                reason="array-api-strict does not promote integer data with a "
+                "floating scalar, so the gain cannot be applied to it "
+                "(https://github.com/astropy/ccdproc/issues/1034)",
+            ),
+        ),
+    ],
+)
+def test_cosmicray_lacosmic_gain_apply_keeps_uncertainty_dtype(dtype):
+    """
+    With ``gain_apply=True`` the uncertainty keeps its own dtype, whatever
+    the dtype of the data (#1023).
+
+    Notes
+    -----
+    The uncertainty is multiplied by the gain, which used to be made a
+    float64 array, so a float32 uncertainty came back as float64 on every
+    backend except JAX. The gain now has the uncertainty's dtype.
+    """
+    ccd_data = ccd_data_func(
+        data_size=40, data_scale=DATA_SCALE, data_mean=1000.0, dtype=dtype
+    )
+    ccd_data.uncertainty = StdDevUncertainty(
+        xp.full(ccd_data.shape, DATA_SCALE, dtype=xp.float32, device=xp_device)
+    )
+
+    result = cosmicray_lacosmic(
+        ccd_data, gain=2.0 * u.electron / u.adu, gain_apply=True
+    )
+
+    assert result.uncertainty.array.dtype == xp.float32
+
+
 def test_cosmicray_median_mask_shape_mismatch():
     # CCDData and MaskedArray validate the mask shape themselves, so exercise
     # the shared helper directly.