Repository: astropy/astropy

Issue:
Impossible to input np.longdouble into custom time formats
I am trying to write a custom time format "mjd_long" which accepts and returns np.longdouble values upon conversion, to make nanosecond accuracy easier. Unfortunately during the preparation stage, before values are handed off to the custom code, they are passed through _make_array, which forcibly converts any non-approved types, including np.longdouble, into np.float64, losing accuracy in the process. As far as I can tell it is a simple matter of allowing higher-precision numerical dtypes (np.float96, np.float128, np.longdouble) through untouched. There is some delicacy about platforms like Windows/MSVC that do not have extended-precision types, where I think some of the type names are not defined. But as far as I can tell this is the only obstacle to writing a custom mjd_long type. (I think with some care plain mjd objects could be made to respect the extra precision of np.longdouble as well, but that's another issue.) 

```python
def _make_array(val, copy=False):
    """
    Take ``val`` and convert/reshape to an array.  If ``copy`` is `True`
    then copy input values.

    Returns
    -------
    val : ndarray
        Array version of ``val``.
    """
    val = np.array(val, copy=copy, subok=True)

    # Allow only float64, string or object arrays as input
    # (object is for datetime, maybe add more specific test later?)
    # This also ensures the right byteorder for float64 (closes #2942).
    if not (val.dtype == np.float64 or val.dtype.kind in 'OSUMaV'):
        val = np.asanyarray(val, dtype=np.float64)

    return val
```


Discussion:
aarchiba: > Thanks! Some smaller comments, but also a general one: ideally, one could just put in higher precision floats and it would "just work". I _think_ it may be possible to do so by changing `day_frac` to special-case higher-precision - it may be possible to do so by adjusting `two_sum` with some calls that would be no-op for regular float (NOTE: NOT TESTED!)
> 
> ```
> def two_sum(a, b):
>     x = np.asanyarray(a + b, float)
>     eb = x - a  # bvirtual in Shewchuk
>     ea = x - eb  # avirtual in Shewchuk
>     eb = b - eb  # broundoff in Shewchuk
>     ea = a - ea  # aroundoff in Shewchuk
>     return x, (ea + eb).astype(float, copy=False)
> ```
> 
> It may be safer to just pre-condition the arguments, a bit like done in `quantity_day_frac`, doing a careful `divmod(val[12], 1)` before adding up the values and residuals.

This would absolutely call for testing with hypothesis, once we figured out exactly what the criterion was for successful computation. 

Well, okay, I think two_sum should do exactly what it does regardless of type it is passed: compute the sum to maximum precision available with the types. But making day_frac operate appropriately would take some doing, especially if we want it to behave well when also passed Quantities. And definitely testing.
mhvk: Yes, makes sense to have `two_sum` be type-agnostic. There is already a `quantity_day_frac` so there is precedent... But, yes, substantial testing needed...
aarchiba: > Yes, makes sense to have `two_sum` be type-agnostic. There is already a `quantity_day_frac` so there is precedent... But, yes, substantial testing needed...

Sorry for the repetition, but my experience with hypothesis in getting the kinks out of PINT's numerical conversion code was extremely positive. It nailed down some difficult-to-spot bugs - including finding leap-second-day problems before I had even told it about leap second days - and programming it was easy (and had that delightful pieces-clicking-together feeling I get with Haskell).
mhvk: I'm all for trying hypothesis on `Time` - I think it is a much easier target than units.
aarchiba: Maybe they're off-topic for this PR but I also found myself poking around to confirm that there is no data validation on the jd1 and jd2 values set in Time objects, and there is no validation of types returned by .value - but there are strict and hard-to-explain requirements on them. This means that implementing custom time formats involves a lot of impenetrable error messages. I'd like to add some validation on both ends; ensure that a Time object's jd1 and jd2 are always arrays of float64s of the same shape no matter what a custom format tries to set them to, and ensure that .value always returns an array of that same shape of whatevers, or a singleton whatever, as appropriate; also that implementors of custom time formats get helpful error messages if their .value methods don't do this correctly. (It'd sort of be nice to be able to provide element-by-element conversion funcions that would automatically be broadcast; these would be too slow for actual numeric arrays but for object arrays they're not any slower and handling arrays of arbitrary shape is awkward; see https://github.com/nanograv/PINT/blob/a4b33eae5fce21a6010e8665bfeecf9caec8cf96/src/pint/pulsar_mjd.py#L501 for the mess.)
aarchiba: So uh.

Right now we quietly ignore it if someone passes a val2 when creating a Time object for a string format. That is, the user provides some information in val2, and astropy throws it away and doesn't tell them. I do not approve.

I added some ValueErrors in this case, only to find that astropy has a test for explicitly this behaviour! Astropy's tests *demand* that we silently throw away information explicitly provided by the user. What? Can I change that? Please?
mhvk: @aarchiba - I agree ignoring `val2` rather than checking it is `None` as appropriate is bad, and have meant to address that myself. So, happy to change, but let's do that in a separate PR!

@taldcroft  and I also discussed that somewhat of a design mistake is to have `jd1` and `jd2` on the `TimeFormat` class - it would make more sense to just let it do the transformations to and from input. 

p.s. Your `mjds_to_str` would seem rather like `Time._shaped_like_input`...
aarchiba: > @aarchiba - I agree ignoring `val2` rather than checking it is `None` as appropriate is bad, and have meant to address that myself. So, happy to change, but let's do that in a separate PR!

I'll pull it out. I have a bad habit of fixing minor problems as they fly past, which gets them fixed but makes my PRs messy.
 
> @taldcroft and I also discussed that somewhat of a design mistake is to have `jd1` and `jd2` on the `TimeFormat` class - it would make more sense to just let it do the transformations to and from input.

This does seem awkward. We are sort of stuck with it, though properties allow for some deviousness.

> p.s. Your `mjds_to_str` would seem rather like `Time._shaped_like_input`...

Ah, yes, probably. I've had to write similar code too many times now, to get functions to behave like ufuncs in terms of correctly maintaining array-ness (scalar->scalar, array->array including zero-dimensional) and shape. It's worse needing to have multiple input arrays that are broadcast together. I thought numpy.vectorize might get me this for free but at least for strings it didn't quite. Modern numpy.vectorize can let a vector-capable function act on an entire one-dimensional array, at least, which broadens its applicability.
mhvk: On ufuncs, see https://github.com/numpy/numpy/issues/13105 - would definitely be good to have a solution...
aarchiba: > On ufuncs, see [numpy/numpy#13105](https://github.com/numpy/numpy/issues/13105) - would definitely be good to have a solution...

I'm not too hopeful, I remember frustration and argument about numpy zero-dimensional arrays versus scalars on the numpy mailing lists more than a decade ago. But it sure would be nice.
aarchiba: Validation moved to #9373 
aarchiba: > > On ufuncs, see [numpy/numpy#13105](https://github.com/numpy/numpy/issues/13105) - would definitely be good to have a solution...
> 
> I'm not too hopeful, I remember frustration and argument about numpy zero-dimensional arrays versus scalars on the numpy mailing lists more than a decade ago. But it sure would be nice.

One of my ancient-times fixes was to record the input array shape (or that it wasn't an array), reshape it into a 1d array, do the vectorized calculation, then reshape it into the same shape (and unpack if it wasn't an array) on output. We could do that, though as we allow Time objects to be indexed and reshaped it would be more complicated. We could at least record whether the input was an array or not so as to know whether to strip off the array boxing on output.
mhvk: I think the shaping, etc., and possibly returning to `np.float64` instead of python float, could be done as part of `_shaped_like_input`. Though I'm not sure in the `Time` case it matters all that much whether it is a `np.float64` or a `float` - I can only see arguments for never returning a python type, just to ensure output always has a `shape`, etc.
mhvk: @aarchiba - on `day_frac`: like `two_sum`, it currently works fine on any type of float - it may make sense to do the conversion at the end of it - in particular, just before where the second rounding is done:
```
day = day.astype(np.float64, copy=False)
frac = frac.astype(np.float64, copy=False)
excess = np.round(frac)
...
```
aarchiba: > @aarchiba - on `day_frac`: like `two_sum`, it currently works fine on any type of float - it may make sense to do the conversion at the end of it - in particular, just before where the second rounding is done:

The day_frac with multiplication/division is not type-agnostic, because two_prod is not - it has the "split" function with a magic constant that would need to be adapted for each type (which I think is not too hard).

More locally, I discovered in writing #9375 that you can't run jd1,jd2 through day_frac upon storage because some tests check for (and presumably some users rely upon?) their not being copied.

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index 706f15e95d..baaed90cce 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -369,6 +369,9 @@ astropy.time
 - Time formats that do not use ``val2`` now raise ValueError instead of
   silently ignoring a provided value. [#9373]
 
+- Custom time formats can now accept floating-point types with extended
+  precision. Existing time formats raise exceptions rather than discarding
+  extended precision through conversion to ordinary floating-point. [#9368]
 
 astropy.timeseries
 ^^^^^^^^^^^^^^^^^^
diff --git a/astropy/time/core.py b/astropy/time/core.py
index a32c42088f..a8c0846def 100644
--- a/astropy/time/core.py
+++ b/astropy/time/core.py
@@ -2285,7 +2285,11 @@ def _make_array(val, copy=False):
     # Allow only float64, string or object arrays as input
     # (object is for datetime, maybe add more specific test later?)
     # This also ensures the right byteorder for float64 (closes #2942).
-    if not (val.dtype == np.float64 or val.dtype.kind in 'OSUMaV'):
+    if val.dtype.kind == "f" and val.dtype.itemsize >= np.dtype(np.float64).itemsize:
+        pass
+    elif val.dtype.kind in 'OSUMaV':
+        pass
+    else:
         val = np.asanyarray(val, dtype=np.float64)
 
     return val
diff --git a/astropy/time/formats.py b/astropy/time/formats.py
index 14a1e7c7f3..edecd83416 100644
--- a/astropy/time/formats.py
+++ b/astropy/time/formats.py
@@ -184,9 +184,9 @@ class TimeFormat(metaclass=TimeFormatMeta):
     def _check_val_type(self, val1, val2):
         """Input value validation, typically overridden by derived classes"""
         # val1 cannot contain nan, but val2 can contain nan
-        ok1 = (val1.dtype == np.double and np.all(np.isfinite(val1)) or
+        ok1 = (val1.dtype.kind == 'f' and val1.dtype.itemsize == 8 and np.all(np.isfinite(val1)) or
                val1.size == 0)
-        ok2 = val2 is None or (val2.dtype == np.double and
+        ok2 = val2 is None or (val2.dtype.kind == 'f' and val2.dtype.itemsize == 8 and
                                not np.any(np.isinf(val2))) or val2.size == 0
         if not (ok1 and ok2):
             raise TypeError('Input values for {} class must be finite doubles'
diff --git a/docs/modeling/add-units.rst b/docs/modeling/add-units.rst
index 79cd735935..a02143839e 100644
--- a/docs/modeling/add-units.rst
+++ b/docs/modeling/add-units.rst
@@ -30,8 +30,6 @@ defaults in the class definition, e.g.::
 
 Note that these are all optional.
 
-.. _models_input_units:
-
 ``input_units``
 ^^^^^^^^^^^^^^^
 
@@ -63,7 +61,7 @@ also be specified as an attribute rather than a property in simple cases::
 ``return_units``
 ^^^^^^^^^^^^^^^^
 
-Similarly to :ref:`models_input_units`, this should be dictionary that maps the return
+Similarly to ``return_units``, this should be dictionary that maps the return
 values of a model to units. If :meth:`~astropy.modeling.Model.evaluate` was called
 with quantities but returns unitless values, the units are added to the output.
 If the return values are quantities in different units, they are converted to
@@ -87,8 +85,8 @@ equivalencies, for example::
 
     input_units_equivalencies = {'nu': u.spectral()}
 
-``_input_units_allow_dimensionless``
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+``input_units_allow_dimensionless``
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 If set to `True`, values that are plain scalars or Numpy arrays can be passed to
 evaluate even if ``input_units`` specifies that the input should have units. It
@@ -101,7 +99,7 @@ maps input names to a Boolean to enable passing dimensionless values to
 Fitting
 -------
 
-To allow models with parameters that have units to be fitted to data with units,
+To allow models with parameters that have units to be fit to data with units,
 you will need to add a method called ``_parameter_units_for_data_units`` to your
 model class. This should take two arguments ``input_units`` and
 ``output_units`` - ``input_units`` will be set to a dictionary with
diff --git a/docs/modeling/changes_for_4.rst b/docs/modeling/changes_for_4.rst
index 18d0ad4d0f..4c3c296044 100644
--- a/docs/modeling/changes_for_4.rst
+++ b/docs/modeling/changes_for_4.rst
@@ -100,14 +100,13 @@ unless you want to understand why some changes were necessary.
 - Many of the non-public methods have changed, particularly for
   compound models.
 
-- The ``_CompoundModelMeta`` metaclass no longer exists.
+- The `_CompoundModelMeta` metaclass no longer exists.
 
+** Technical Details **
 
-Technical Details
-=================
 
-Parameter-related changes
--------------------------
+** Parameter-related changes **
+
 
 Previously Python descriptors were used to define parameters. The
 drawback of descriptors is that all instances of the models that
@@ -154,8 +153,8 @@ Since the parameter no longer is required to link to a specific
 model, it holds any dimension corresponding to the model_set
 size; when it previously did not.
 
-Compound Model Implementation Changes
--------------------------------------
+** Compound Model Implementation Changes **
+
 
 Compound models previously were implemented using a metaclass
 for the compound model while also inheriting from the Model


Test patch:
diff --git a/astropy/time/tests/test_custom_formats.py b/astropy/time/tests/test_custom_formats.py
index 4b8d09b953..1a41cde10b 100644
--- a/astropy/time/tests/test_custom_formats.py
+++ b/astropy/time/tests/test_custom_formats.py
@@ -1,7 +1,13 @@
-import pytest
 from itertools import count
 
+import pytest
+
+import numpy as np
+
+import astropy.units as u
+from astropy._erfa import DJM0
 from astropy.time import Time, TimeFormat
+from astropy.time.utils import day_frac
 
 
 class SpecificException(ValueError):
@@ -76,6 +82,8 @@ def test_custom_time_format_fine(custom_format_name):
 
     t = Time.now()
     getattr(t, custom_format_name)
+    t2 = Time(7, 9, format=custom_format_name)
+    getattr(t2, custom_format_name)
 
 
 def test_custom_time_format_forgot_property(custom_format_name):
@@ -127,3 +135,56 @@ def test_custom_time_format_problematic_name():
 
     finally:
         Time.FORMATS.pop("sort", None)
+
+
+def test_mjd_longdouble_preserves_precision(custom_format_name):
+    class CustomMJD(TimeFormat):
+        name = custom_format_name
+
+        def _check_val_type(self, val, val2):
+            val = np.longdouble(val)
+            if val2 is not None:
+                raise ValueError("Only one value permitted")
+            return val, 0
+
+        def set_jds(self, val, val2):
+            mjd1 = np.float64(np.floor(val))
+            mjd2 = np.float64(val - mjd1)
+            self.jd1, self.jd2 = day_frac(mjd1 + DJM0, mjd2)
+
+        @property
+        def value(self):
+            mjd1, mjd2 = day_frac(self.jd1 - DJM0, self.jd2)
+            return np.longdouble(mjd1) + np.longdouble(mjd2)
+
+    m = 58000.0
+    t = Time(m, format=custom_format_name)
+    t2 = Time(m + 2 * m * np.finfo(np.longdouble).eps, format=custom_format_name)
+    assert t != t2
+    assert isinstance(getattr(t, custom_format_name), np.longdouble)
+    assert getattr(t, custom_format_name) != getattr(t2, custom_format_name)
+
+
+def test_mjd_unit_validation():
+    with pytest.raises(u.UnitConversionError):
+        Time(58000 * u.m, format="mjd")
+
+
+def test_mjd_unit_conversion():
+    assert Time(58000 * u.day, format="mjd") == Time(58000 * u.day, format="mjd")
+    assert Time(58000 * u.day, format="mjd") != Time(58000 * u.s, format="mjd")
+    assert Time(58000 * u.day, format="mjd") == Time(58000 * 86400 * u.s, format="mjd")
+
+
+@pytest.mark.parametrize("f", ["mjd", "unix", "cxcsec"])
+def test_existing_types_refuse_longdoubles(f):
+    t = np.longdouble(getattr(Time(58000, format="mjd"), f))
+    t2 = t + np.finfo(np.longdouble).eps * 2 * t
+    try:
+        tm = Time(np.longdouble(t), format=f)
+    except ValueError:
+        # Time processing makes it always ValueError not TypeError
+        return
+    else:
+        # accepts long doubles, better preserve accuracy!
+        assert Time(np.longdouble(t2), format=f) != tm