Complete cosmo_array numpy support and add cosmo_quantity

.github/workflows/pytest.yml

@@ -16,7 +16,7 @@ jobs:
       with:
         python-version: ${{ matrix.python-version }}
     - name: Cache pip
-      uses: actions/cache@v2
+      uses: actions/cache@v4
       with:
         # This path is specific to Ubuntu
         path: ~/.cache/pip

docs/source/conf.py

@@ -40,6 +40,7 @@
     "sphinx.ext.napoleon",
     "sphinx.ext.mathjax",
     "sphinx.ext.autosummary",
+    "sphinx.ext.intersphinx",
 ]
 
 # Add any paths that contain templates here, relative to this directory.
@@ -98,3 +99,13 @@ def run_apidoc(_):
 
 def setup(app):
     app.connect("builder-inited", run_apidoc)
+
+
+intersphinx_mapping = dict(
+    numpy=("https://numpy.org/doc/stable/", None),
+    numba=("https://numba.readthedocs.io/en/stable/", None),
+    unyt=("https://unyt.readthedocs.io/en/stable/", None),
+    scipy=("https://docs.scipy.org/doc/scipy/", None),
+    swiftgalaxy=("https://swiftsimio.readthedocs.io/en/latest/", None),
+    velociraptor=("https://velociraptor-python.readthedocs.io/en/latest/", None),
+)

docs/source/cosmo_array/index.rst

@@ -0,0 +1,91 @@
+The ``cosmo_array``
+===================
+
+:mod:`swiftsimio` uses a customized class based on the :class:`~unyt.array.unyt_array`
+to store data arrays. The :class:`~swiftsimio.objects.cosmo_array` has all of the same
+functionality as the :class:`~unyt.array.unyt_array`, but also adds information about
+data transformation between physical and comoving coordinates, and descriptive metadata.
+
+For instance, should you ever need to know what a dataset represents, you can
+ask for a description by accessing the ``name`` attribute:
+
+.. code-block:: python
+
+   print(rho_gas.name)
+
+which will output ``Co-moving mass densities of the particles``.
+
+The cosmology information is stored in three attributes:
+
+ + ``comoving``
+ + ``cosmo_factor``
+ + ``valid_transform``
+
+The ``comoving`` attribute specifies whether the array is a physical (``True``) or
+comoving (``False``) quantity, while the ``cosmo_factor`` stores the expression needed
+to convert back and forth between comoving and physical quantities and the value of
+the scale factor. The conversion factors can be accessed like this:
+
+.. code-block:: python
+
+   # Conversion factor to make the densities a physical quantity
+   print(rho_gas.cosmo_factor.a_factor)
+   physical_rho_gas = rho_gas.cosmo_factor.a_factor * rho_gas
+
+   # Symbolic scale-factor expression
+   print(rho_gas.cosmo_factor.expr)
+
+which will output ``132651.002785671`` and ``a**(-3.0)``. Converting an array to/from physical/comoving
+is done with the :meth:`~swiftsimio.objects.cosmo_array.to_physical`, :meth:`~swiftsimio.objects.cosmo_array.to_comoving`, :meth:`~swiftsimio.objects.cosmo_array.convert_to_physical` and :meth:`~swiftsimio.objects.cosmo_array.to_comoving` methods, for instance:
+
+.. code-block:: python
+
+   physical_rho_gas = rho_gas.to_physical()
+
+   # Convert in-place
+   rho_gas.convert_to_physical()
+
+The ``valid_transform`` is a boolean flag that is set to ``False`` for some arrays that don't make sense to convert to comoving.
+
+:class:`~swiftsimio.objects.cosmo_array` supports array arithmetic and the entire :mod:`numpy` range of functions. Attempting to combine arrays (e.g. by addition) will validate the cosmology information first. The implementation is designed to be permissive: it will only raise exceptions when a genuinely invalid combination is encountered, but is tolerant of missing cosmology information. When one argument in a relevant operation (like addition, for example) is not a :class:`~swiftsimio.objects.cosmo_array` the attributes of the :class:`~swiftsimio.objects.cosmo_array` will be assumed for both arguments. In such cases a warning is produced stating that this assumption has been made.
+
+.. note::
+
+   :class:`~swiftsimio.objects.cosmo_array` and the related :class:`~swiftsimio.objects.cosmo_quantity` are now intended to support all :mod:`numpy` functions, propagating units and cosmology information correctly through mathematical operations. Try making a histogram with weights and ``density=True`` with :func:`numpy.histogram`! There are a large number of functions and a very large number of possible parameter combinations, so some corner cases may have been missed in development. Please report any errors or unexpected results using github issues or other channels so that they can be fixed. Currently :mod:`scipy` functions are not supported (although some might "just work"). Requests to support specific functions can be accommodated.
+
+To make the most of the utility offered by the :class:`~swiftsimio.objects.cosmo_array` class, it is helpful to know how to create your own. A good template for this looks like:
+
+.. code-block:: python
+
+   import unyt as u
+   from swiftsimio.objects import cosmo_array, cosmo_factor
+
+   # suppose the scale factor is 0.5 and it scales as a**1, then: 
+   my_cosmo_array = cosmo_array(
+       [1, 2, 3],
+       u.Mpc,
+       comoving=True,
+       scale_factor=0.5,  # a=0.5, i.e. z=1
+       scale_exponent=1,  # distances scale as a**1, so the scale exponent is 1
+   )
+   # consider getting the scale factor from metadata when applicable, i.e. replace:
+   # scale_factor=0.5
+   # with:
+   # scale_factor=data.metadata.a
+
+There is also a very similar :class:`~swiftsimio.objects.cosmo_quantity` class designed for scalar values,
+analogous to the :class:`~unyt.array.unyt_quantity`. You may encounter this being returned by :mod:`numpy` functions. Cosmology-aware scalar values can be initialized similarly:
+
+.. code-block:: python
+
+   import unyt as u
+   from swiftsimio.objects import cosmo_quantity, cosmo_factor
+
+   my_cosmo_quantity = cosmo_quantity(
+       2,
+       u.Mpc,
+       comoving=False,
+       scale_factor=0.5,
+       cosmo_factor=1,
+   )
+

docs/source/getting_started/index.rst

@@ -109,8 +109,8 @@ In the above it's important to note the following:
 + Only the density and temperatures (corresponding to the ``PartType0/Densities`` and
   ``PartType0/Temperatures``) datasets are read in.
 + That data is only read in once the
-  :meth:`swiftsimio.objects.cosmo_array.convert_to_cgs` method is called.
-+ :meth:`swiftsimio.objects.cosmo_array.convert_to_cgs` converts data in-place;
+  :meth:`~swiftsimio.objects.cosmo_array.convert_to_cgs` method is called.
++ :meth:`~swiftsimio.objects.cosmo_array.convert_to_cgs` converts data in-place;
   i.e. it returns `None`.
 + The data is cached and not re-read in when ``plt.scatter`` is called.
 

docs/source/index.rst

@@ -17,8 +17,10 @@ snapshots to enable partial reading.
 
    getting_started/index
    loading_data/index
+   cosmo_array/index
    masking/index
    visualisation/index
+   soap/index
    velociraptor/index
    creating_initial_conditions/index
    statistics/index

docs/source/loading_data/index.rst

@@ -247,47 +247,6 @@ Then, to access individual columns (in this case element abundances):
    data.gas.element_mass_fractions.silicon
 
 
-Non-unyt properties
--------------------
-
-Each data array has some custom properties that are not present within the base
-:obj:`unyt.unyt_array` class. We create our own version of this in
-:obj:`swiftsimio.objects.cosmo_array`, which allows each dataset to contain
-its own cosmology and name properties.
-
-For instance, should you ever need to know what a dataset represents, you can
-ask for a description:
-
-.. code-block:: python
-
-   print(rho_gas.name)
-
-which will output ``Co-moving mass densities of the particles``. They include
-scale-factor information, too, through the ``cosmo_factor`` object,
-
-.. code-block:: python
-
-   # Conversion factor to make the densities a physical quantity
-   print(rho_gas.cosmo_factor.a_factor)
-   physical_rho_gas = rho_gas.cosmo_factor.a_factor * rho_gas
-
-   # Symbolic scale-factor expression
-   print(rho_gas.cosmo_factor.expr)
-
-which will output ``132651.002785671`` and ``a**(-3.0)``. This is an easy way
-to convert your co-moving values to physical ones.
-
-An even easier way to convert your properties to physical is to use the
-built-in ``to_physical`` and ``convert_to_physical`` methods, as follows:
-
-.. code-block:: python
-
-   physical_rho_gas = rho_gas.to_physical()
-
-   # Convert in-place
-   rho_gas.convert_to_physical()
-
-
 User-defined particle types
 ---------------------------
 
@@ -310,38 +269,3 @@ in SWIFT will be automatically read.
        "extra_test.hdf5",
    )
 
-
-Halo Catalogues
----------------
-
-SWIFT-compatible halo catalogues, such as those written with SOAP, can be
-loaded entirely transparently with ``swiftsimio``. It is generally possible
-to use all of the functionality (masking, visualisation, etc.) that is used
-with snapshots with these files, assuming the files conform to the
-correct metadata standard.
-
-An example SOAP file is available at
-``http://virgodb.cosma.dur.ac.uk/swift-webstorage/IOExamples/soap_example.hdf5``
-
-You can load SOAP files as follows:
-
-.. code-block:: python
-
-   from swiftsimio import load
-
-   catalogue = load("soap_example.hdf5")
-
-   print(catalogue.spherical_overdensity_200_mean.total_mass)
-
-   # >>> [  591.      328.5     361.      553.      530.      507.      795.
-   #        574.      489.5     233.75      0.     1406.      367.5    2308.
-   #        ...
-   #        0.      534.        0.      191.75   1450.      600.      290.   ] 10000000000.0*Msun (Physical)
-
-What's going on here? Under the hood, ``swiftsimio`` has a discrimination function
-between different metadata types, based upon a property stored in the HDF5 file,
-``Header/OutputType``. If this is set to ``FullVolume``, we have a snapshot,
-and use the :obj:`swiftsimio.metadata.objects.SWIFTSnapshotMetadata`
-class. If it is ``SOAP``, we use
-:obj:`swiftsimio.metadata.objects.SWIFTSOAPMetadata`, which instructs
-``swiftsimio`` to read slightly different properties from the HDF5 file.

docs/source/soap/index.rst

@@ -0,0 +1,41 @@
+Halo Catalogues & SOAP integration
+==================================
+
+SWIFT-compatible halo catalogues, such as those written with SOAP, can be
+loaded entirely transparently with ``swiftsimio``. It is generally possible
+to use all of the functionality (masking, visualisation, etc.) that is used
+with snapshots with these files, assuming the files conform to the
+correct metadata standard.
+
+An example SOAP file is available at
+``http://virgodb.cosma.dur.ac.uk/swift-webstorage/IOExamples/soap_example.hdf5``
+
+You can load SOAP files as follows:
+
+.. code-block:: python
+
+   from swiftsimio import load
+
+   catalogue = load("soap_example.hdf5")
+
+   print(catalogue.spherical_overdensity_200_mean.total_mass)
+
+   # >>> [  591.      328.5     361.      553.      530.      507.      795.
+   #        574.      489.5     233.75      0.     1406.      367.5    2308.
+   #        ...
+   #        0.      534.        0.      191.75   1450.      600.      290.   ] 10000000000.0*Msun (Physical)
+
+What's going on here? Under the hood, ``swiftsimio`` has a discrimination function
+between different metadata types, based upon a property stored in the HDF5 file,
+``Header/OutputType``. If this is set to ``FullVolume``, we have a snapshot,
+and use the :obj:`swiftsimio.metadata.objects.SWIFTSnapshotMetadata`
+class. If it is ``SOAP``, we use
+:obj:`swiftsimio.metadata.objects.SWIFTSOAPMetadata`, which instructs
+``swiftsimio`` to read slightly different properties from the HDF5 file.
+
+swiftgalaxy
+-----------
+
+The :mod:`swiftgalaxy` companion package to :mod:`swiftsimio` offers further integration with halo catalogues in SOAP, Caesar and Velociraptor formats (so far). It greatly simplifies efficient loading of particles belonging to an object from a catalogue, and additional tools that are useful when working with a galaxy or other localized collection of particles. Refer to the `swiftgalaxy documentation`_ for details.
+
+.. _swiftgalaxy documentation: https://swiftgalaxy.readthedocs.io/en/latest/

docs/source/visualisation/index.rst

@@ -1,12 +1,12 @@
 Visualisation
 =============
 
-:mod:`swiftsimio` provides visualisation routines accelerated with the
-:mod:`numba` module. They work without this module, but we strongly recommend
-installing it for the best performance (1000x+ speedups). These are provided
-in the :mod:`swiftismio.visualisation` sub-modules.
+:mod:`swiftsimio` provides visualisation routines in the
+:mod:`swiftsimio.visualisation` sub-module.  They are accelerated with the
+:mod:`numba` module. They can work without :mod:`numba`, but we strongly recommend
+installing it for the best performance (1000x+ speedups).
 
-The three built-in rendering types (described below) have the following
+The four built-in rendering types (described below) have the following
 common interface:
 
 .. code-block:: python
@@ -30,6 +30,7 @@ additional functionality.
    projection
    slice
    volume_render
+   ray_trace
    power_spectra
    tools