Preserve non-float dtypes in Quantity columns from serialized mixins - #12505

👋 Thank you for your draft pull request! Do you know that you can use [ci skip] or [skip ci] in your commit messages to skip running continuous integration tests until you are ready?

There's a bit of a problem with this:

In [5]: t = Table([Time([1,2], format='cxcsec')])
In [6]: t.write('junk.ecsv', overwrite=True)
In [7]: t2 = Table.read('junk.ecsv')
In [8]: t2['col0'].info.dtype
Out[8]: dtype('float64')   # <== dtype('O') in main

For mixin classes like Time, the dtype is left unset which then defaults to object.

I think a better approach is to change how Quantity serializes itself to add a dtype object attribute:

Out[12]: 
<QTable length=2>
 col0
  m  
int64
-----
    1
    2

In [13]: ascii.write(qt, format='ecsv')
# %ECSV 1.0
# ---
# datatype:
# - {name: col0, unit: m, datatype: int64}
# meta: !!omap
# - __serialized_columns__:
#     col0:
#       __class__: astropy.units.quantity.Quantity
#       unit: !astropy.units.Unit {unit: m}
#       value: !astropy.table.SerializedColumn {name: col0}
#       dtype: int64                                         <=== ADD THIS
# schema: astropy-2.0
col0
1
2

This seems to work:

diff --git a/astropy/table/serialize.py b/astropy/table/serialize.py
index 5dfa402818..62d69ccf44 100644
--- a/astropy/table/serialize.py
+++ b/astropy/table/serialize.py
@@ -342,7 +342,6 @@ def _construct_mixin_from_columns(new_name, obj_attrs, out):
 
         # Now copy the relevant attributes
         for attr, nontrivial in (('unit', lambda x: x not in (None, '')),
-                                 ('dtype', lambda x: x is not None),
                                  ('format', lambda x: x is not None),
                                  ('description', lambda x: x is not None),
                                  ('meta', lambda x: x)):
diff --git a/astropy/units/quantity.py b/astropy/units/quantity.py
index 2c60005a84..35b01159e5 100644
--- a/astropy/units/quantity.py
+++ b/astropy/units/quantity.py
@@ -171,10 +171,16 @@ class QuantityInfo(QuantityInfoBase):
     required when the object is used as a mixin column within a table, but can
     be used as a general way to store meta information.
     """
-    _represent_as_dict_attrs = ('value', 'unit')
     _construct_from_dict_args = ['value']
     _represent_as_dict_primary_data = 'value'
 
+    def _represent_as_dict(self, attrs=None):
+        q = self._parent
+        out = {'value': q.value,
+               'unit': q.unit,
+               'dtype': q.dtype.name}
+        return out
+
     def new_like(self, cols, length, metadata_conflicts='warn', name=None):
         """
         Return a new Quantity instance which is consistent with the

In [3]: qt = QTable([u.Quantity([1,2], unit=u.m, dtype=int)])

In [4]: qt.write('junk.ecsv', overwrite=True)

In [5]: qt2 = Table.read('junk.ecsv')

In [6]: qt2
Out[6]: 
<Table length=2>
 col0
  m  
int64
-----
    1
    2

(astropy) ➜  astropy git:(decode-dtype-mixins) ✗ cat junk.ecsv
# %ECSV 1.0
# ---
# datatype:
# - {name: col0, unit: m, datatype: int64}
# meta: !!omap
# - __serialized_columns__:
#     col0:
#       __class__: astropy.units.quantity.Quantity
#       dtype: int64
#       unit: !astropy.units.Unit {unit: m}
#       value: !astropy.table.SerializedColumn {name: col0}
# schema: astropy-2.0
col0
1
2

Certainly a bug that the dtype does not round-trip. I like @taldcroft's solution! I wondered whether it is worth not storing the dtype if it float64 (i.e., the default), but checking quickly, I saw that Column also stores its dtype in the header. So, I think this is better generally - really seems more like an oversight that we did not do this already!

This seems to work:

+    def _represent_as_dict(self, attrs=None):
+        q = self._parent
+        out = {'value': q.value,
+               'unit': q.unit,
+               'dtype': q.dtype.name}
+        return out
+

Thanks, this looks much neater and more concise than my approach. And I guess keeping datatype to record the format the columns were written in, while meta stores the format and dtype they are supposed to be returned in, has some logic to it. The only drawback would be that this does not decode files written with earlier versions.
And I found one test failure on the Longitude mixin, due to the wrap_angle not being included with the definition above.
A naïve attempt to just add it to the dict as

             out['wrap_angle'] = {'value': q.wrap_angle.value, 'unit': q.wrap_angle.unit}

fails on decoding with a ValueError: Invalid character at col 0 in angle 'unit', since this is not fully reproducing the previous description

COMMENT     lon:
COMMENT       __class__: astropy.coordinates.angles.Longitude
COMMENT       unit: &id001 !astropy.units.Unit {unit: deg}
COMMENT       value: !astropy.table.SerializedColumn {name: lon}
COMMENT       wrap_angle: !astropy.coordinates.Angle
COMMENT         unit: *id001
COMMENT         value: 360.0

but instead outputs

COMMENT     lon:
COMMENT       __class__: astropy.coordinates.angles.Longitude
COMMENT       dtype: float32
COMMENT       unit: &id001 !astropy.units.Unit {unit: deg}
COMMENT       value: !astropy.table.SerializedColumn {name: lon}
COMMENT       wrap_angle:
COMMENT         unit: *id001
COMMENT         value: 360.0

i.e. the !astropy.coordinates.Angle still needs to be reconstructed, but how? Also, need to make sure there are no other attributes missed. I did not find any failures pointing to others in the io, coordinates or units tests though.

But I do see several more failures in table I cannot really make sense of:

___________________________________________________________________________________ test_vstack ____________________________________________________________________________________

    def test_vstack():
        """
        Vstack tables with mixin cols.
        """
        t1 = QTable(MIXIN_COLS)
        t2 = QTable(MIXIN_COLS)
        with pytest.raises(NotImplementedError):
>           vstack([t1, t2])

opt/lib/python3.10/site-packages/astropy/table/tests/test_mixin.py:504: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
opt/lib/python3.10/site-packages/astropy/table/operations.py:651: in vstack
    out = _vstack(tables, join_type, col_name_map, metadata_conflicts)
opt/lib/python3.10/site-packages/astropy/table/operations.py:1396: in _vstack
    col = col_cls.info.new_like(cols, n_rows, metadata_conflicts, out_name)
opt/lib/python3.10/site-packages/astropy/units/quantity.py:224: in new_like
    out = self._construct_from_dict(map)
opt/lib/python3.10/site-packages/astropy/utils/data_info.py:387: in _construct_from_dict
    args = [map.pop(attr) for attr in self._construct_from_dict_args]
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

.0 = <list_iterator object at 0x12adeda20>

>   args = [map.pop(attr) for attr in self._construct_from_dict_args]
E   KeyError: 'value'

opt/lib/python3.10/site-packages/astropy/utils/data_info.py:387: KeyError
...
============================================================================= short test summary info ==============================================================================
FAILED opt/lib/python3.10/site-packages/astropy/table/tests/test_mixin.py::test_vstack - KeyError: 'value'
FAILED opt/lib/python3.10/site-packages/astropy/table/tests/test_operations.py::TestJoin::test_col_meta_merge[QTable] - KeyError: 'value'
FAILED opt/lib/python3.10/site-packages/astropy/table/tests/test_operations.py::TestVStack::test_col_meta_merge_inner[QTable] - KeyError: 'value'
FAILED opt/lib/python3.10/site-packages/astropy/table/tests/test_operations.py::TestVStack::test_mixin_functionality[latitude] - KeyError: 'value'
FAILED opt/lib/python3.10/site-packages/astropy/table/tests/test_operations.py::TestVStack::test_mixin_functionality[longitude] - KeyError: 'value'
FAILED opt/lib/python3.10/site-packages/astropy/table/tests/test_operations.py::TestVStack::test_mixin_functionality[quantity] - KeyError: 'value'
FAILED opt/lib/python3.10/site-packages/astropy/table/tests/test_operations.py::test_stack_columns - KeyError: 'value'
FAILED opt/lib/python3.10/site-packages/astropy/table/tests/test_operations.py::test_mixin_join_regression - KeyError: 'value'

(all with the same KeyError in map.pop(attr) for attr in self._construct_from_dict_args, but why is 'value' not present in those cases???)

Stealing the _represent_as_dict implementation from representations seems to do.

More general comment: there have been quite a few issues about Quantity automatically converting to float, and whether or not that is desired. My sense remains that generally it is desired, as an integer quantity makes little physical sense (apart from spin!). I think there is a lot to be said for avoiding silent truncation errors from, e.g., in-place operations.

Given this, I wonder if it isn't better to use some kind of option for the cases where one really wants to preserve the data type with which the data are stored internally. After all, conversion to float is not that different from applying BZERO and BSCALE to integer image data.

More general comment: there have been quite a few issues about Quantity automatically converting to float, and whether or not that is desired. My sense remains that generally it is desired, as an integer quantity makes little physical sense (apart from spin!). I think there is a lot to be said for avoiding silent truncation errors from, e.g., in-place operations.

Yes, though to me the biggest concern seemed that many in-place operations would simply fail on integer quantities.

Put it another way, operations that may automatically cast to float, in order of decreasing desirability, to me would be

1. u.Quantity(arr)
2. QTable(tab), QTable([columns])
3. QTable.read(file)
4. QTable(tab, dtype=tab.dtype)
5. QTable.read(file, do_not_scale_table_data=True)

where the last option would offer a yet to be implemented option analogous to do_not_scale_image_data for ensuring that the data are read in as stored. That was actually a kind of minimal solution I should have pointed out a bit more in the initial comment. But the original specifically issue arose from the large memory footprint of data read in from file, so there should be some form of control at that step.

@dhomeier - good to see your summary - I agree with the order, and that (4) and (5) for sure should respect the dtype. I think I'm actually also on-board with (3) , though maybe prefer (5) with the keyword argument to control it (as in the difference between (2) and (4)).

Also, stopping the auto-conversion to float (i.e., (3)) may well break code that relies on QTable doing that...

This version implements (3) - (5) now.

There are pros and cons for adding 'f' (and 'c'?) to the dtype kinds to be serialised; the wrap_angle case pointed out above is probably not very relevant, and both keeping read-in float columns in big-endian format and automatically converting them to native has its own points in favour.

Also, stopping the auto-conversion to float (i.e., (3)) may well break code that relies on QTable doing that...

That's probably still worth discussing; I'd think that since creating a QTable with integer quantities already requires some very explicit steps, people would generally be aware of the format. But it might be safer to keep the old behaviour as default.
On adding an option to QTable.read() for controlling this: both outcomes can already be effected relatively easily by forcing conversion with QTable(Table.read(file)) in the current version, or
t = Table.read(file); qt = QTable(t, dtype=t.dtype) if conversion is made the default.

I've checked locally this works with #12509, so this should be safe to rebase once either PR is merged. Need to revert the yaml docstring changes anyway.
Pinging @saimn and @nstarman for opinions.

I've checked locally this works with #12509, so this should be safe to rebase once either PR is merged. Need to revert the yaml docstring changes anyway.

Merged! So good to rebase.

More general comment: there have been quite a few issues about Quantity automatically converting to float, and whether or not that is desired. My sense remains that generally it is desired, as an integer quantity makes little physical sense (apart from spin!). I think there is a lot to be said for avoiding silent truncation errors from, e.g., in-place operations.

I have more mixed feelings about this. When wearing my science hat I agree completely: rarely is the int dtype desirable and if it's necessary should be explicitly stated.

• stays int: Quantity([1, 2, 3], dtype=int)
• upcast to float: Quantity([1, 2, 3])

When wearing my comp-sci hat, the thought of auto-upcasting is less attractive. Quantity is an extension of ndarray and this is a departure from their treatment of dtype. At the very least we should more explicitly document this. I don't feel I was sufficiently aware of these issues until reading the discussion in this PR.

Given this, I wonder if it isn't better to use some kind of option for the cases where one really wants to preserve the data type with which the data are stored internally. After all, conversion to float is not that different from applying BZERO and BSCALE to integer image data.

I do think there might be a comfortable middle ground:

For a list x = [1, 2, 3], using the following constructor

The Ellipsis option would essentially mean "pass it on". Numpy would be the end-point here, so a Longitude(x, dtype=...) -> Angle -> SpecificTypeQuantity -> Quantity -> ndarray, and only all the Quantity subclasses understand to delegate dtype to super, until it hits numpy. When passed to numpy the Ellipsis would be changed to None to use numpy's dtype detection machinery.

The advantage of this approach is that there is a simple means of specifying the dtype without knowing what it is. It's essentially equivalent to

x = [1, 2, 3]
arr = np.array(x)
q = Quantity(arr, dtype=arr.dtype)

but allows the user to skip the explicit construction of an intermediate array.

If we apply this consistently, I hope it can solve (or at least alleviate) @dhomeier's list of upcast situations (#12505 (comment)). 🤞
Each class can have appropriate (& well documented) default rules for upcasting, and allow for an easy override to force the dtypes to be kept unchanged.

Rebased (it works!) to get the discussion going again.

When wearing my comp-sci hat, the thought of auto-upcasting is less attractive. Quantity is an extension of ndarray and this is a departure from their treatment of dtype. At the very least we should more explicitly document this. I don't feel I was sufficiently aware of these issues until reading the discussion in this PR.

The unexpected upcasting was also the initial issue triggering this, so the use case is certainly relevant enough.
I do quite like that proposal @nstarman!

The Ellipsis option would essentially mean "pass it on". Numpy would be the end-point here, so a Longitude(x, dtype=...) -> Angle -> SpecificTypeQuantity -> Quantity -> ndarray, and only all the Quantity subclasses understand to delegate dtype to super, until it hits numpy. When passed to numpy the Ellipsis would be changed to None to use numpy's dtype detection machinery.

The advantage of this approach is that there is a simple means of specifying the dtype without knowing what it is. It's essentially equivalent to

x = [1, 2, 3]
arr = np.array(x)
q = Quantity(arr, dtype=arr.dtype)

but allows the user to skip the explicit construction of an intermediate array.

Just wondering how many in places code needs to be changed to handle the dtype='...' setting...

And for Quantity that would need a decision whether '...' should become the new default (still a change in behaviour), or keep that at None+upcast. But even being able to simply call Quantity(array_like, dtype='...') without having to bother what the actual dtype is would be a great advantage IMO.

One point from #15447 which I opened and then learned that this PR already exists: I'm looking at specific fits tables that are the standard for all of X-ray astronomy, a so-called ARF/RMF/PHA files. All three of those file types are tables where most columns have good, physical fits units, e.g. keV or counts or counts/s or cm**2; thus I'd like to read those files as a QTable. However, there are other columns that need to remain integers; in particular the "channel". A channel is by definition a discrete quantity and downstream software will fail with floating point numbers of channels; also channel numbers are frequently used as indexes into the other vector-valued columns (e.g. in an RMF, there a "response matrix". Since this matrix is very sparse, it's not stored in a rectangular form, but instead only the non-zero elements of a row are stored in the file and other columns give the integer indices for that, e.g. a row in the matrix might be non-zero for elements 567 to 681, and then there would be a column that says N_start=567 - obviously, that column should be an integer).

Without explaining the format in more detail here, there is no way to represent this as a QTable without the ability to have some of the columns as an int.

@hamogu - yes, this PR was a good one that we should try to revive.

One question: does the channel column have a unit?

Also, for a possible rebase: what @nstarman ended up implementing was that u.Quantity(..., dtype=None) is what is needed to avoid the upcast (not the ... mentioned in the comments above; the change that we ended up with is to change the default to dtype=np.inexact to indicate auto-convert to float).

p.s. Apologies, but I won't be able to spend time here before feature-freeze: I'm stuck writing my next 5-year research grant until Nov 1st.

It has a unit set in the fits header that does not parse with astropy, probably because it's not part of the main fits standard, but only defined by NASA's HEASARC fits working group (formerly "Office of guest investigator programs") (e.g. https://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/spectra/ogip_92_007/node7.html for PHA type I files)

WARNING: UnitsWarning: 'channel' did not parse as fits unit: At col 0, Unit 'channel' not supported by the FITS standard.

(Links given for future reference, not because it's needed to follow all of them to understand what I need Quantity to do to make this work.)

Just rebasing this with the black/merge conflicts resolved; I'll need more time to look into @mhvk's last round of comments, too.

Removing the milestone for now, as it seems unlikely this will make 7.0

I would love to have the VOTable example from #17963 boiled down into some test cases to be included here.

@dhomeier - will you have time to try to push this to the finish line? I think we were getting close...

@dhomeier @taldcroft - is this still relevant/needed now #20273 is merged?

Superseded by #20273.