Cosmo Unified I/O - #11559

Hello @nstarman 👋! It looks like you've made some changes in your pull request, so I've checked the code again for style.

There are no PEP8 style issues with this pull request - thanks! 🎉

Comment last updated at 2021-07-01 20:27:31 UTC

👋 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?

@nstarman - if you could point me to where the MaskedConstant managed to come in, I can try to help figure out where it came from!

Thanks!!

In this branch the following will result in a TypeError: unhashable type: 'MaskedConstant'

from astropy.cosmology import WMAP5
WMAP5.write("temp.csv", format="ascii.csv")
WMAP5.__class__.read("temp.csv", format="ascii.csv")

Actually, the following formats all have this problem:
"rdb","ascii.fast_rdb","ascii.basic","ascii.fast_basic","ascii","ascii.rdb","ascii.fast_tab","ascii.csv","csv","ascii.fast_commented_header","pandas.csv","ascii.commented_header","ascii.fast_csv","ascii.tab",

OK, so at least QTable.read("temp.csv", format="ascii.csv") works fine, so the problem at least is local to cosmology...

And name is a MaskedConstant, which is from numpy (though I use it too). I think the problem is that the name column is empty, which gets seen as a masked entry.

@nstarman - I see this PR is using the generalized TableRead/Write classes to allow writing the cosmology parameters to any Table-supported format. Clever! My question is whether nominally supporting all possible formats is really a good idea. I had thought you were just going to do ECSV, with the understanding that PyYAML will be a required dependency of astropy.

From my perspective if you are looking to invent a new way of serializing the cosmology, then less is more and you should pick exactly one format (noting that most of the Table formats don't support serializing units or metadata).

Disclaimer: I've never actually used the cosmology subpackage, so I may be missing something.

My question is whether nominally supporting all possible formats is really a good idea. I had thought you were just going to do ECSV, with the understanding that PyYAML will be a required dependency of astropy.
From my perspective if you are looking to invent a new way of serializing the cosmology, then less is more and you should pick exactly one format (noting that most of the Table formats don't support serializing units or metadata).

My original intent was to only implement ECSV and JSON.
ECSV will be the preferred method of Cosmology I/O because it supports units, metadata, and arrays of cosmologies.
I feel JSON is also necessary for eventual to_, from_ with CLASS, CAMB, and other cosmology packages that don't use Astropy. I haven't examined this yet in depth, but it appears that the pandas.json table I/O might be sufficient for this use case. 🤞.

When I switched from a home-brewed solution to using the generalized TableRead/Write my thought was that it would be easier maintenance-wise as the methods will always be up-to-date with the Table I/O and this would also ensure a consistent function signature with Table.read/write for writing ECSV and JSON.

The cosmology package has been careful to allow for float and Quantity input for all arguments. If a user is not interested in roundtrip-ability (on the metadata) nor automatic unit handling, then basically any Table format is technically valid input, even if not advised. I agree that less is more, but by subclassing `TableRead/Write``, which is kind of necessary for the ECSV and JSON anyways, all the other table formats are taken along for the ride.

The ECSV is most definitely the preferred format as it allows for saving citations.
Astropy realizations will be saved as ECSV and auto-import of realizations (#11517 ) will only support ECSV and JSON.
But I want very human-readable formats like fixed-width tab-separated tables or HTML tables or even latex to also be valid input since those are more commonly embedded online / in a paper.

Edit:

Perhaps a good solution would be to have to/from_table methods. read and write would then be limited to ECSV, JSON (pandas.json), and also for read: table instances, and dictionaries (the latter two are already supported). Now anyone brave enough to initialize a cosmology realization from an HTML table would go file -> Table.read -> Cosmology.from_table. And for writing a Cosmology to HTML would do cosmo_realization.to_table -> table.write -> file.

I think this is a good solution to allow for flexible I/O and also keep the recommended formats prominent.

@nstarman - some thoughts...

First, the 🤞 plan for supporting JSON won't easily work in general since I see that m_nu is an array parameter in Planck15. Pandas does not support multidimensional columns. In fact ECSV does not support multidim columns before 4.3.

Is there really a strong use-case for tables of cosmologies such that it needs to be supported as a code one-liner instead of potentially writing a few lines of code? In particular I think current behavior for Cosmology.read, where it returns a list of Cosmology objects if the table has more than 1 row, is undesirable. The unified I/O API is that cls.read() should return an instance of cls, so returning a list is unexpected.

I'm not sure that using this unified I/O read infrastructure with formats like latex and HTML will be useful in the wild. Basically it seems like the column names are likely to differ enough (esp. in the latex case) to end up requiring some by-hand munging. Then you just have e.g.

pars = Table.read('cosmo.tex')[0]
c_pars['H0'] = pars['H_0']
...
c = Cosmology(**c_pars)

I recommend against including functionality that is "not advised" or partially working, since somebody will start doing that and then come to report problems. Ideally everything that goes into the core should be completely robust and bulletproof.

About human readability, I think the ECSV is sufficiently readable:

H0 Om0 Tcmb0 Neff m_nu Ob0 name
67.74 0.3075 2.7255 3.046 [0.0,0.0,0.06] 0.0486 null

With astropy installed a quick showtable plank15.ecsv will give a nice result. (Remember that the fixed_width format does not support multidimensional columns). Non-astropy users can pull the table into TOPCAT or their favorite spreadsheet.

I've really come around to a philosophy of favoring minimal solutions that are simple. Then wait to see if users to submit issues like "I really need a fixed-format output for cosmology parameters".

When I switched from a home-brewed solution to using the generalized TableRead/Write my thought was that it would be easier maintenance-wise as the methods will always be up-to-date with the Table I/O and this would also ensure a consistent function signature with Table.read/write for writing ECSV and JSON.

The recent case of QDP is a counter-example, where a new ASCII format got added that doesn't work with Cosmology. But more to the the point, I don't think of Cosmology as being Table-like. Certainly having read/write methods in the unified I/O scheme is a very good thing, but inheriting from Table is not the right way IMO.

At the end my vote would be to support exactly two formats, each of which can serialize a single cosmology:

• ECSV (1-row table). Might be useful to enforce canonical units for readers that are not able to parse the header.
• JSON, where this is output as a simple dict of parameters. All values are in canonical units.

I've not thought through much of this, but @taldcroft's comment that cosmology is not a table did make me think that perhaps one should instead think of cosmology as an object that can be serialized and thus can be meta data on anything else. Conceptually, that seems to make more sense; e.g., it could be meta data on a table of supernova observations and inferred luminosity distances, etc.

I'm not quite sure how to extend that to also serialize just the cosmology. I huess it would need an .info attribute and a fake shape=(1,), which seems a bit odd.

Anyway, overall very much 👍 to not reinventing the wheel on a yet another serialization method that is astropy.cosmology-specific, but let's think a bit more what would be the best way to approach this.

@mhvk - yes, putting .info seems like an interesting option to consider. And certainly serializing Cosmology to YAML would be low-hanging fruit and allow Cosmology objects to be included in any meta. The main issue there is that this YAML output is quite astropy-specific. (In particular a numpy array gets serialized to an encoded bytestring).

As far as a universally-accessible serialization format for a cosmology... First is there any prior art? Beyond that, honestly a simple JSON dict seems like the best choice, mostly driven by array-valued parameters. Simple metadata like references can be easily stored.

With ECSV there is the problem of m_nu that shows up as [0.0,0.0,0.06]. Any vanilla CSV reader will happily read that as the string "[0.0,0.0,0.06]", but this doesn't help most people.

First, the 🤞 plan for supporting JSON won't easily work in general since I see that m_nu is an array parameter in Planck15. Pandas does not support multidimensional columns.

Yeah in further testing this is proving troublesome. Before I switched to the TableRead/Write I had a simple to/from_json using the built-in json library. It was essentially a 1-liner to pass to the from_mapping method that is also necessary for ECSV (table <-> dict <-> cosmology).

Is there really a strong use-case for tables of cosmologies such that it needs to be supported as a code one-liner instead of potentially writing a few lines of code? In particular I think current behavior for Cosmology.read, where it returns a list of Cosmology objects if the table has more than 1 row, is undesirable. The unified I/O API is that cls.read() should return an instance of cls, so returning a list is unexpected.

I think this is the fine balance between over-engineering and not wanting to back myself into a corner later on. In prior research I've had to run cosmological MCMC chains. Many programs then allow one to move along the chain, performing calculations. While these chains are usually specific to the program running them because the parameters aren't labelled, if re-represented as a Table my thought is astropy.cosmology could do much the same thing. The list implementation was a temporary measure before I implement a CosmologySet or some such. But the more I think about it, I think this is the wrong approach entirely. Currently Cosmology only takes scalar inputs on the parameters (except m_nu). What if it were vectorized and subscriptable, so the following was possible FlatLambdaCDM(H0=[67, 68, 69, 70], ...).
now an MCMC chain of cosmological parameters can be read in by Cosmology and have calculations performed on all / any subset of the chain.

I'm not sure that using this unified I/O read infrastructure with formats like latex and HTML will be useful in the wild. Basically it seems like the column names are likely to differ enough (esp. in the latex case) to end up requiring some by-hand munging. Then you just have e.g.
I recommend against including functionality that is "not advised" or partially working, since somebody will start doing that and then come to report problems. Ideally everything that goes into the core should be completely robust and bulletproof.

Good points both. I think that this really applies to the read method. Having something in the write method work for outputting latex, HTML, etc is still very useful for papers and such. So long as it's documented that Cosmology -> not ECSV/JSON is a one way trip this obviates the issue. Below I propose having the write class have a to_table method.

About human readability, I think the ECSV is sufficiently readable:

It's very readable for single cosmologies. Less so for many rows with different numbers of decimal places.
But agreed this is probably an edge case.

With astropy installed a quick showtable plank15.ecsv will give a nice result. (Remember that the fixed_width format does not support multidimensional columns). Non-astropy users can pull the table into TOPCAT or their favorite spreadsheet.
I've really come around to a philosophy of favoring minimal solutions that are simple. Then wait to see if users to submit issues like "I really need a fixed-format output for cosmology parameters".

👍

Certainly having read/write methods in the unified I/O scheme is a very good thing, but inheriting from Table is not the right way IMO.
At the end my vote would be to support exactly two formats, each of which can serialize a single cosmology:

• ECSV (1-row table). Might be useful to enforce canonical units for readers that are not able to parse the header.
• JSON, where this is output as a simple dict of parameters. All values are in canonical units.

I'm convinced. My proposal is this:
For now:

• Implement the read/write using registry.UnifiedReadWrite, not Table's subclasses. The read method will be restricted to ECSV and JSON. The write method will be similarly restricted.
• For ECSV it will be a 1-row table, but read allows for the row to be specified as an argument so that a file with many cosmologies can still be read in.
• In addition to the standard __call__, the write class will also expose the two utility methods to_mapping and to_table so that cosmo.write.to_mapping and cosmo.write.to_table can convert a cosmology to a dict/table. What users do with a Table is at their discretion and any possible write failures are at the Table level, not buried in the cosmology code.
• In addition to the standard __call__ (for files), the read class will also expose the two utility methods from_mappingand from_table so that cosmo.read.from_mapping and cosmo.write.from_table can instantiate a cosmology.

For later, when Cosmology is vectorized:

• Allow I/O of N-row tables

[@mhvk] Perhaps one should instead think of cosmology as an object that can be serialized and thus can be meta data on anything else. Conceptually, that seems to make more sense; e.g., it could be meta data on a table of supernova observations and inferred luminosity distances, etc.
I'm not quite sure how to extend that to also serialize just the cosmology. I [g]uess it would need an .info attribute and a fake shape=(1,), which seems a bit odd.

I really like this! I'm definitely interested in implementing this. Correct me if I'm misinterpreting, but this would be in addition to the above serialization implementation (esp for a "vectorized" Cosmology)? I think the shape meshes nicely with vectorizing Cosmology.

[@taldcroft ] As far as a universally-accessible serialization format for a cosmology... First is there any prior art? Beyond that, honestly a simple JSON dict seems like the best choice, mostly driven by array-valued parameters. Simple metadata like references can be easily stored.

I have not encountered a universally-accessible serialization format. I think they are all bespoke. CLASS uses a .ini YAML file. I agree JSON is the easiest, especially since cosmologies now store their parameters in a dictionary (._init_parameters)
and this (+meta) is trivially exported to JSON.

@mhvk @taldcroft. Thanks for the detailed reviews and the tough questions. I appreciate all the insight into the astropy.io package and how to integrate the read/write registries into the astropy.cosmology package.

I’m going to try to tackle all points (#11559 (review) and similar) in one long-form response

Goals

1. A lossless reading/writing of cosmology objects from/to files in human readable formats. Metadata must be preserved.
2. A minimal-loss format using only built-in packages for future interoperability with 3rd party packages.
3. Easy extensibility, for both users and developers, to other file formats. In particular, readers / writers to 3rd party package formats.
4. Not muddy the namespace — e.g. having a to_X, from_X for each format.
5. Transformation methods of objects (e.g. Cosmology <-> Table), that likewise satisfy the above goals of being minimally lossy (hopefully lossless), easily extensible, and not being burdensome on the namespace. This will enable Cosmology <-> other package object.

Solutions

• ECSV is the clear file format choice to accomplish goal 1.
• JSON appears to be the best choice to accomplish goal 2.
• The I/O registry is excellent for goal 3 and 5. The method pairread/write manage everything.
• The I/O registry is also great for goal 4! To the I/O registry, there’s essentially no difference between a file and an object. (the only change needed is a single added return to not break the pass-through chain for write methods.)

Everything will use the Unified I/O methods, including transforming to/from objects like dictionaries and Tables. The descriptors on Cosmology will be the only way to read/write, with the io module made private and only used for registering in the readers/writers.

For example (the property names from_, to_ are definitely NOT final):

class Cosmology:
    # Unified I/O object transformation methods from ``astropy.cosmology._io``
    # The 'from_' method is NotImplemented in subclasses.
    from_ = UnifiedReadWriteMethod(CosmologyFrom)  # can’t use ‘from =’
    to_ = UnifiedReadWriteMethod(CosmologyTo)

    # Unified I/O read and write methods from ``astropy.cosmology._io``
    # The 'read' method is NotImplemented in subclasses.
    read = UnifiedReadWriteMethod(CosmologyRead)
    write = UnifiedReadWriteMethod(CosmologyWrite)

Now converting to/from a mapping/table is as easy as

from astropy.cosmology import Cosmology, Planck18

mapping = Planck18.to_("mapping")
cosmo = Cosmology.from_(mapping)
assert cosmo == Planck18  # True

table = Planck18.to_("table")
cosmo = Cosmology.from_(table)
assert cosmo == Planck18  # True

from astropy.table import vstack
newtable = vstack([table, table])
cosmo = Cosmology.from_(newtable, index=1)
assert cosmo == Planck18  # True

And new formats can be easily added, without adding anything to the Cosmology object itself! Also it’s quite literally a registry so 3rd party packages can contribute formats — e.g. to_("cosmo_jax").

This PR now contains a reference implementation of the above.

Discussion

[Just using Table I/O]

Cosmology file format requirements are different than that of Tables. For one, other popular packages employ formats like .ini that are incompatible with Table, but should be readable by Cosmology. Second, Cosmology needs to be able to serialize None, which cannot be done in some Table formats. Last, column names must correspond to Cosmology (sub)class arguments, but formats like latex and html should be math-formatted, so there must be some custom mapping between the two.

Partially it is because this PR does (too) many things in one go:
I think it would be much easier if we could do things in turn, separating table/ECSV from JSON...

Agreed, this PR has become too large. I’m going to split it up into pieces, if we agree that the above is the right approach. This PR branch is called cosmo_io_JSON, so I’m going to make this PR for JSON and precede it with a PRs for mapping instances and then do PRs for Table instances and ECSV files.

extend the custom JSON encoder we already have in utils.misc

I would like to support this in addition to base JSON. For a followup.

make pyyaml a required depen[den]cy for 5.0 and stick with ECSV…

Done on the first part and ECSV will be the preferred format here. JSON is need for interoperability with non-astropy-supporting packages

I wonder if eventually that isn't better done through a new .info property (this really depends on whether Cosmology could perhaps gain a shape).

I will be adding an .info property, but this will use the Cosmology.to_/from_. Also when Cosmology is eventually vectorized it will also gain a shape.

I think the bottom line is a question of what Cosmology really is; arguably it is closer to a Model than a SkyCoord…

Agreed! I want to eventually make stuff like H0 into something analogous to Parameter from astropy.modeling.

As a new maintainer for cosmology I think it is reasonable to give @nstarman some leeway to implement his vision for the interface. This version of the PR is already the 2nd entire rewrite based on my feedback, so I do think it's a bit late to ask for a 3rd rewrite as this has been in work for quite some time.

Thanks! But I think third time's the charm :).

I’m still doubtful about the methods on read and write - the functions are defined and exposed in io, why not just make that the primary method?

I wanted them for 2 reasons:

1. Convenience.
2. Having a method that doesn’t require importing anything from astropy so that 3rd party packages can have from_astropy without astropy as a dependency.

I agree that having 2 ways of doing things was an implementation mistake. Now it’s all in the class method.

@nstarman - thanks, that really helps! You definitely completely sold me on the read/write aspect, with gaining the ability to read and write data in formats used by other packages particularly compelling. I would actually suggest to start with that, since it also helps for your other very nice idea to store the parameters for the standard cosmologies for astropy itself.

For JSON, I can see the case too (and it fits with read/write anyway), but mostly if there is already some use outside astropy - if there isn't, I think there should first be some discussion with other packages, including a careful look at whether there is an existing format that is "good enough" (obligatory link: https://xkcd.com/927/).

The use case for the methods is fairly clear but the implementation I'm much less sure about; name pollution depends on the number of to/from pairs that you expect.

Anyway, for what it is worth, my suggested ordering would be to start with a slightly narrowed down PR from this one with just read/write. Next may depend on whether there is an existing human-readable standard that could encapsulate all our cosmological parameters. If so, incorporate it and use it to store ours. If not, use ECSV.

p.s. I'll add that I really appreciate your willingness to discuss and re-do - @taldcroft and I often similarly take our time to converge to a solution (like .info for mixin columns and their serialization), but I've similarly enjoyed our collaboration on representations, and you stand out for not getting irritated (or at least not showing it!). It's part of what makes open-source projects enjoyable!

but I've similarly enjoyed our collaboration on representations, and you stand out for not getting irritated (or at least not showing it!).

Same! And the former.

The use case for the methods is fairly clear but the implementation I'm much less sure about; name pollution depends on the number of to/from pairs that you expect.

Could you elaborate about the drawbacks?

Re name pollution, I expect to have to/from_mapping, to/from_table, to/from_class, to/from_camb, to/from_colossus, to/from_cosmosis, to/from_cosmo_jax, and probably more (maybe to/from_numpy for MCMC chains and to/from_arviz for probability stuff). The to_format('mapping/table/...'), from_format(obj) compacts everything to these two methods. In particular, if I can work with the other packages to get these registrations implemented on their side, that would dramatically reduce the maintenance burden for offering interoperability. My hope is that other packages would want to register themselves for easy use with, e.g. an MCMC routine.

It is partially names - to_format('camb') would be nicer (to me) than to_('camb'), and I can see 'mapping' and 'table' could fit in there. There's a question how best to do it. If you want to_format, that means you'll have to import those packages to create the class, and you're stuck with ensuring you initialize them correctly. For from_format instead, you can just have knowledge about their attributes, and in principle, you could even do it in Cosmology.__new__. We've had somewhat similar discussions about compatibility between different quantity/unit implementations, where my sense was that each package should be responsible primarily for being able to import from other packages, not also for exporting to them. But obviously, in the end it depends on the relation to the other package: for pandas, table does have to and from.

Ideally, of course, you have some form of interchange format - perhaps something like __array__ or __array_interface__ for numerical arrays - though my guess is that an ASCII on-disk format may be the easiest to start with, especially if there is already a format supported by multiple packages.

to_format('camb') would be nicer (to me) than to_('camb')

Agreed. The to_ was temporary while I tried different names. to/from_format is the best thus far.

If you want to_format, that means you'll have to import those packages to create the class, and you're stuck with ensuring you initialize them correctly

My hope is to actually do none/few/only the 'from_' of these in Astropy itself, but work with package developers to register in a reader/writer/identifier. The to/from_mapping makes this all relatively straightforward. One package (https://bdiemer.bitbucket.io/colossus/cosmology_cosmology.html) already has a to/from Astropy, so that would be a good place to start for the inverse methods.

you could even do it in Cosmology.new

I like that as an initializer. I think the constructor might look something like this (pseudo) code.

def __new__(cls, *args, **kwargs):
    if args[0].__module__ is not astropy:
        return cls.from_format(args[0])
    else:
        ...

an ASCII on-disk format may be the easiest to start with, especially if there is already a format supported by multiple packages.

Unfortunately, though many support .ini files, that is such a flexible format that all of them do it quite incompatibly.

Ah, yes, now I finally see the connection with readers... A registration mechanism would definitely seem the way to go!

This PR will be implemented in a number of smaller PRs. Thanks @mhvk and @taldcroft for the great discussion and feedback that will be motivating these PRs.