Add Quijote (#61)

* Rename paths object

* Remove redshift calculation

* Explicit keywrod arg

* Rename box units

* Basic renaming

* Little docs

* Rename paths

* add imports

* Sort imports

* Add Quijote cat

* Split boxes

* add Quijote path

* Add origin argument

* Update nbs

csiborgtools/__init__.py

@@ -14,7 +14,8 @@
 # 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 from csiborgtools import clustering, field, fits, match, read  # noqa
 
-# Arguments to csiborgtools.read.CSiBORGPaths.
+# Arguments to csiborgtools.read.Paths.
 paths_glamdring = {"srcdir": "/mnt/extraspace/hdesmond/",
                    "postdir": "/mnt/extraspace/rstiskalek/CSiBORG/",
+                   "quijote_dir": "/mnt/extraspace/rstiskalek/Quijote",
                    }

csiborgtools/field/density.py

@@ -52,7 +52,7 @@ class BaseField(ABC):
 
         Returns
         -------
-        box : :py:class:`csiborgtools.units.BoxUnits`
+        box : :py:class:`csiborgtools.units.CSiBORGBox`
         """
         return self._box
 
@@ -94,7 +94,7 @@ class DensityField(BaseField):
 
     Parameters
     ----------
-    box : :py:class:`csiborgtools.read.BoxUnits`
+    box : :py:class:`csiborgtools.read.CSiBORGBox`
         The simulation box information and transformations.
     MAS : str
         Mass assignment scheme. Options are Options are: 'NGP' (nearest grid
@@ -204,7 +204,7 @@ class VelocityField(BaseField):
 
     Parameters
     ----------
-    box : :py:class:`csiborgtools.read.BoxUnits`
+    box : :py:class:`csiborgtools.read.CSiBORGBox`
         The simulation box information and transformations.
     MAS : str
         Mass assignment scheme. Options are Options are: 'NGP' (nearest grid
@@ -295,7 +295,7 @@ class PotentialField(BaseField):
 
     Parameters
     ----------
-    box : :py:class:`csiborgtools.read.BoxUnits`
+    box : :py:class:`csiborgtools.read.CSiBORGBox`
         The simulation box information and transformations.
     MAS : str
         Mass assignment scheme. Options are Options are: 'NGP' (nearest grid
@@ -334,7 +334,7 @@ class TidalTensorField(BaseField):
 
     Parameters
     ----------
-    box : :py:class:`csiborgtools.read.BoxUnits`
+    box : :py:class:`csiborgtools.read.CSiBORGBox`
         The simulation box information and transformations.
     MAS : str
         Mass assignment scheme. Options are Options are: 'NGP' (nearest grid

csiborgtools/field/interp.py

@@ -67,7 +67,7 @@ def evaluate_sky(*fields, pos, box, isdeg=True):
     pos : 2-dimensional array of shape `(n_samples, 3)`
         Spherical coordinates to evaluate the field. Columns are distance,
         right ascension, declination, respectively.
-    box : :py:class:`csiborgtools.read.BoxUnits`
+    box : :py:class:`csiborgtools.read.CSiBORGBox`
         The simulation box information and transformations.
     isdeg : bool, optional
         Whether `ra` and `dec` are in degres. By default `True`.

csiborgtools/fits/halo.py

@@ -67,7 +67,7 @@ class BaseStructure(ABC):
 
         Returns
         -------
-        box : :py:class:`csiborgtools.units.BoxUnits`
+        box : :py:class:`csiborgtools.units.CSiBORGBox`
         """
         return self._box
 
@@ -280,7 +280,7 @@ class Clump(BaseStructure):
         Particle array. Must contain `['x', 'y', 'z', 'vx', 'vy', 'vz', 'M']`.
     info : structured array
         Array containing information from the clump finder.
-    box : :py:class:`csiborgtools.read.BoxUnits`
+    box : :py:class:`csiborgtools.read.CSiBORGBox`
         Box units object.
     """
 
@@ -301,7 +301,7 @@ class Halo(BaseStructure):
         Particle array. Must contain `['x', 'y', 'z', 'vx', 'vy', 'vz', 'M']`.
     info : structured array
         Array containing information from the clump finder.
-    box : :py:class:`csiborgtools.read.BoxUnits`
+    box : :py:class:`csiborgtools.read.CSiBORGBox`
         Box units object.
     """
 

csiborgtools/read/__init__.py

@@ -12,14 +12,14 @@
 # You should have received a copy of the GNU General Public License along
 # with this program; if not, write to the Free Software Foundation, Inc.,
 # 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
-from .box_units import BoxUnits  # noqa
+from .box_units import CSiBORGBox, QuijoteBox  # noqa
-from .halo_cat import ClumpsCatalogue, HaloCatalogue  # noqa
+from .halo_cat import ClumpsCatalogue, HaloCatalogue, QuijoteHaloCatalogue  # noqa
 from .knn_summary import kNNCDFReader  # noqa
 from .obs import (SDSS, MCXCClusters, PlanckClusters, TwoMPPGalaxies,  # noqa
                   TwoMPPGroups)
 from .overlap_summary import (NPairsOverlap, PairOverlap,  # noqa
                               binned_resample_mean)
-from .paths import CSiBORGPaths  # noqa
+from .paths import Paths  # noqa
 from .pk_summary import PKReader  # noqa
 from .readsim import (MmainReader, ParticleReader, halfwidth_mask,  # noqa
                       load_clump_particles, load_parent_particles, read_initcm)

csiborgtools/read/box_units.py

@@ -15,14 +15,15 @@
 """
 Simulation box unit transformations.
 """
+from abc import ABC
+
 import numpy
 from astropy import constants, units
 from astropy.cosmology import LambdaCDM
-from scipy.interpolate import interp1d
 
 from .readsim import ParticleReader
 
-# Map of unit conversions
+# Map of CSiBORG unit conversions
 CONV_NAME = {
     "length": ["x", "y", "z", "peak_x", "peak_y", "peak_z", "Rs", "rmin",
                "rmax", "r200c", "r500c", "r200m", "x0", "y0", "z0",
@@ -33,53 +34,18 @@ CONV_NAME = {
     "density": ["rho0"]}
 
 
-class BoxUnits:
+###############################################################################
-    r"""
+#                              Base box                                       #
-    Box units class for converting between box and physical units.
+###############################################################################
 
-    Paramaters
+
-    ----------
+class BaseBox(ABC):
-    nsnap : int
+    """
-        Snapshot index.
+    Base class for box units.
-    nsim : int
-        IC realisation index.
-    paths : py:class`csiborgtools.read.CSiBORGPaths`
-        CSiBORG paths object.
     """
     _name = "box_units"
     _cosmo = None
 
-    def __init__(self, nsnap, nsim, paths):
-        """
-        Read in the snapshot info file and set the units from it.
-        """
-        partreader = ParticleReader(paths)
-        info = partreader.read_info(nsnap, nsim)
-        pars = [
-            "boxlen",
-            "time",
-            "aexp",
-            "H0",
-            "omega_m",
-            "omega_l",
-            "omega_k",
-            "omega_b",
-            "unit_l",
-            "unit_d",
-            "unit_t",
-        ]
-        for par in pars:
-            setattr(self, "_" + par, float(info[par]))
-
-        self._cosmo = LambdaCDM(
-            H0=self._H0,
-            Om0=self._omega_m,
-            Ode0=self._omega_l,
-            Tcmb0=2.725 * units.K,
-            Ob0=self._omega_b,
-        )
-        self._Msuncgs = constants.M_sun.cgs.value  # Solar mass in grams
-
     @property
     def cosmo(self):
         """
@@ -89,6 +55,8 @@ class BoxUnits:
         -------
         cosmo : `astropy.cosmology.LambdaCDM`
         """
+        if self._cosmo is None:
+            raise ValueError("Cosmology not set.")
         return self._cosmo
 
     @property
@@ -101,7 +69,7 @@ class BoxUnits:
         -------
         H0 : float
         """
-        return self._H0
+        return self.cosmo.H0.value
 
     @property
     def h(self):
@@ -112,7 +80,54 @@ class BoxUnits:
         -------
         h : float
         """
-        return self._H0 / 100
+        return self.H0 / 100
+
+    @property
+    def Om0(self):
+        r"""
+        The matter density parameter.
+
+        Returns
+        -------
+        Om0 : float
+        """
+        return self.cosmo.Om0
+
+
+###############################################################################
+#                              CSiBORG box                                    #
+###############################################################################
+
+
+class CSiBORGBox(BaseBox):
+    r"""
+    CSiBORG box units class for converting between box and physical units.
+
+    Paramaters
+    ----------
+    nsnap : int
+        Snapshot index.
+    nsim : int
+        IC realisation index.
+    paths : py:class`csiborgtools.read.Paths`
+        CSiBORG paths object.
+    """
+
+    def __init__(self, nsnap, nsim, paths):
+        """
+        Read in the snapshot info file and set the units from it.
+        """
+        partreader = ParticleReader(paths)
+        info = partreader.read_info(nsnap, nsim)
+        pars = ["boxlen", "time", "aexp", "H0", "omega_m", "omega_l",
+                "omega_k", "omega_b", "unit_l", "unit_d", "unit_t"]
+        for par in pars:
+            setattr(self, "_" + par, float(info[par]))
+
+        self._cosmo = LambdaCDM(H0=self._H0, Om0=self._omega_m,
+                                Ode0=self._omega_l, Tcmb0=2.725 * units.K,
+                                Ob0=self._omega_b)
+        self._Msuncgs = constants.M_sun.cgs.value  # Solar mass in grams
 
     @property
     def box_G(self):
@@ -157,7 +172,6 @@ class BoxUnits:
         -------
         rhoc : float
         """
-
         return 3 * self.box_H0**2 / (8 * numpy.pi * self.box_G)
 
     def box2kpc(self, length):
@@ -244,83 +258,6 @@ class BoxUnits:
         """
         return vel * (1e-2 * self._unit_l / self._unit_t / self._aexp) * 1e-3
 
-    def box2cosmoredshift(self, dist):
-        r"""
-        Convert the box comoving distance to cosmological redshift.
-
-        NOTE: this likely is already the observed redshift.
-
-        Parameters
-        ----------
-        dist : float
-            Distance in box units.
-
-        Returns
-        -------
-        cosmo_redshift : foat
-            Cosmological redshift.
-        """
-        x = numpy.linspace(0.0, 1.0, 5001)
-        y = self.cosmo.comoving_distance(x)
-        return interp1d(y, x)(self.box2mpc(dist))
-
-    def box2pecredshift(self, vx, vy, vz, px, py, pz, p0x=0, p0y=0, p0z=0):
-        """
-        Convert the box phase-space information to a peculiar redshift.
-
-        NOTE: there is some confusion about this.
-
-        Parameters
-        ----------
-        vx, vy, vz : 1-dimensional arrays
-            Cartesian velocity components.
-        px, py, pz : 1-dimensional arrays
-            Cartesian position vectors components.
-        p0x, p0y, p0z : floats
-            Centre of the box coordinates. By default 0, in which it is assumed
-            that the coordinates are already centred.
-
-        Returns
-        -------
-        pec_redshift : 1-dimensional array
-            Peculiar redshift.
-        """
-        # Peculiar velocity along the radial distance
-        r = numpy.vstack([px - p0x, py - p0y, pz - p0z]).T
-        norm = numpy.sum(r**2, axis=1) ** 0.5
-        v = numpy.vstack([vx, vy, vz]).T
-        vpec = numpy.sum(r * v, axis=1) / norm
-        # Ratio between the peculiar velocity and speed of light
-        x = self.box2vel(vpec) / constants.c.value
-        # Doppler shift
-        return numpy.sqrt((1 + x) / (1 - x)) - 1
-
-    def box2obsredshift(self, vx, vy, vz, px, py, pz, p0x=0, p0y=0, p0z=0):
-        """
-        Convert the box phase-space information to an 'observed' redshift.
-
-        NOTE: there is some confusion about this.
-
-        Parameters
-        ----------
-        vx, vy, vz : 1-dimensional arrays
-            Cartesian velocity components.
-        px, py, pz : 1-dimensional arrays
-            Cartesian position vectors components.
-        p0x, p0y, p0z : floats
-            Centre of the box coordinates. By default 0, in which it is assumed
-            that the coordinates are already centred.
-
-        Returns
-        -------
-        obs_redshift : 1-dimensional array
-            Observed redshift.
-        """
-        r = numpy.vstack([px - p0x, py - p0y, pz - p0z]).T
-        zcosmo = self.box2cosmoredshift(numpy.sum(r**2, axis=1) ** 0.5)
-        zpec = self.box2pecredshift(vx, vy, vz, px, py, pz, p0x, p0y, p0z)
-        return (1 + zpec) * (1 + zcosmo) - 1
-
     def solarmass2box(self, mass):
         r"""
         Convert mass from :math:`M_\odot` (with :math:`h=0.705`) to box units.
@@ -391,7 +328,7 @@ class BoxUnits:
         return (density / self._unit_d * self._Msuncgs
                 / (units.Mpc.to(units.cm)) ** 3)
 
-    def convert_from_boxunits(self, data, names):
+    def convert_from_box(self, data, names):
         r"""
         Convert columns named `names` in array `data` from box units to
         physical units, such that
@@ -421,13 +358,12 @@ class BoxUnits:
         transforms = {"length": self.box2mpc,
                       "mass": self.box2solarmass,
                       "velocity": self.box2vel,
-                      "density": self.box2dens,
+                      "density": self.box2dens}
-                      }
 
         for name in names:
             # Check that the name is even in the array
             if name not in data.dtype.names:
-                raise ValueError("Name `{}` not in `data` array.".format(name))
+                raise ValueError(f"Name `{name}` not in `data` array.")
 
             # Convert
             found = False
@@ -439,11 +375,36 @@ class BoxUnits:
             # If nothing found
             if not found:
                 raise NotImplementedError(
-                    "Conversion of `{}` is not defined.".format(name)
+                    f"Conversion of `{name}` is not defined.")
-                )
 
             # Center at the observer
             if name in ["peak_x", "peak_y", "peak_z", "x0", "y0", "z0"]:
                 data[name] -= transforms["length"](0.5)
 
         return data
+
+
+###############################################################################
+#                      Quijote fiducial cosmology box                         #
+###############################################################################
+
+
+class QuijoteBox(BaseBox):
+    """
+    Quijote fiducial cosmology box.
+
+    Parameters
+    ----------
+    nsnap : int
+        Snapshot number.
+    **kwargs : dict
+        Empty keyword arguments. For backwards compatibility.
+    """
+
+    def __init__(self, nsnap, **kwargs):
+        zdict = {4: 0.0, 3: 0.5, 2: 1.0, 1: 2.0, 0: 3.0}
+        assert nsnap in zdict.keys(), f"`nsnap` must be in {zdict.keys()}."
+        self._aexp = 1 / (1 + zdict[nsnap])
+
+        self._cosmo = LambdaCDM(H0=67.11, Om0=0.3175, Ode0=0.6825,
+                                Tcmb0=2.725 * units.K, Ob0=0.049)

csiborgtools/read/halo_cat.py

@@ -12,24 +12,29 @@
 # You should have received a copy of the GNU General Public License along
 # with this program; if not, write to the Free Software Foundation, Inc.,
 # 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
-"""CSiBORG halo and clumps catalogues."""
+"""
-from abc import ABC
+Simulation catalogues:
+    - CSiBORG: halo and clump catalogue.
+    - Quijote: halo catalogue.
+"""
+from abc import ABC, abstractproperty
+from os.path import join
 
 import numpy
+from readfof import FoF_catalog
 from sklearn.neighbors import NearestNeighbors
 
-from .box_units import BoxUnits
+from .box_units import CSiBORGBox, QuijoteBox
-from .paths import CSiBORGPaths
+from .paths import Paths
 from .readsim import ParticleReader
-from .utils import (add_columns, cartesian_to_radec, flip_cols,
+from .utils import (add_columns, cartesian_to_radec, cols_to_structured,
-                    radec_to_cartesian, real2redshift)
+                    flip_cols, radec_to_cartesian, real2redshift)
 
 
 class BaseCatalogue(ABC):
     """
     Base (sub)halo catalogue.
     """
-
     _data = None
     _paths = None
     _nsim = None
@@ -51,14 +56,25 @@ class BaseCatalogue(ABC):
     def nsim(self, nsim):
         self._nsim = nsim
 
-    @property
+    @abstractproperty
-    def paths(self):
+    def nsnap(self):
         """
-        CSiBORG paths manager.
+        Catalogue's snapshot index.
 
         Returns
         -------
-        paths : :py:class:`csiborgtools.read.CSiBORGPaths`
+        nsnap : int
+        """
+        pass
+
+    @property
+    def paths(self):
+        """
+        Paths manager.
+
+        Returns
+        -------
+        paths : :py:class:`csiborgtools.read.Paths`
         """
         if self._paths is None:
             raise RuntimeError("`paths` is not set!")
@@ -66,7 +82,7 @@ class BaseCatalogue(ABC):
 
     @paths.setter
     def paths(self, paths):
-        assert isinstance(paths, CSiBORGPaths)
+        assert isinstance(paths, Paths)
         self._paths = paths
 
     @property
@@ -82,36 +98,16 @@ class BaseCatalogue(ABC):
             raise RuntimeError("Catalogue data not loaded!")
         return self._data
 
-    @property
+    @abstractproperty
-    def nsnap(self):
-        """
-        Catalogue's snapshot index corresponding to the maximum simulation
-        snapshot index.
-
-        Returns
-        -------
-        nsnap : int
-        """
-        return max(self.paths.get_snapshots(self.nsim))
-
-    @property
     def box(self):
         """
-        CSiBORG box object handling unit conversion.
+        Box object.
 
         Returns
         -------
-        box : :py:class:`csiborgtools.units.BoxUnits`
+        box : instance of :py:class:`csiborgtools.units.BoxUnits`
         """
-        return BoxUnits(self.nsnap, self.nsim, self.paths)
+        pass
-
-    @box.setter
-    def box(self, box):
-        try:
-            assert box._name == "box_units"
-            self._box = box
-        except AttributeError as err:
-            raise TypeError from err
 
     def position(self, in_initial=False, cartesian=True):
         r"""
@@ -143,8 +139,8 @@ class BaseCatalogue(ABC):
         return pos
 
     def velocity(self):
-        """
+        r"""
-        Cartesian velocity components in box units.
+        Cartesian velocity components in :math:`\mathrm{km} / \mathrm{s}`.
 
         Returns
         -------
@@ -204,7 +200,7 @@ class BaseCatalogue(ABC):
         knn : :py:class:`sklearn.neighbors.NearestNeighbors`
         """
         knn = NearestNeighbors()
-        return knn.fit(self.position(in_initial))
+        return knn.fit(self.position(in_initial=in_initial))
 
     def radius_neigbours(self, X, radius, in_initial):
         r"""
@@ -299,20 +295,56 @@ class BaseCatalogue(ABC):
 
     @property
     def keys(self):
-        """Catalogue keys."""
+        """
+        Catalogue keys.
+
+        Returns
+        -------
+        keys : list of strings
+        """
         return self.data.dtype.names
 
     def __getitem__(self, key):
-        initpars = ["x0", "y0", "z0"]
+        if key not in self.keys:
-        if key in initpars and key not in self.keys:
+            raise KeyError(f"Key '{key}' not in catalogue.")
-            raise RuntimeError("Initial positions are not set!")
         return self.data[key]
 
     def __len__(self):
         return self.data.size
 
 
-class ClumpsCatalogue(BaseCatalogue):
+###############################################################################
+#                       CSiBORG  base catalogue                               #
+###############################################################################
+
+
+class BaseCSiBORG(BaseCatalogue):
+    """
+    Base CSiBORG catalogue class.
+    """
+
+    @property
+    def nsnap(self):
+        return max(self.paths.get_snapshots(self.nsim))
+
+    @property
+    def box(self):
+        """
+        CSiBORG box object handling unit conversion.
+
+        Returns
+        -------
+        box : instance of :py:class:`csiborgtools.units.BaseBox`
+        """
+        return CSiBORGBox(self.nsnap, self.nsim, self.paths)
+
+
+###############################################################################
+#                        CSiBORG clumps catalogue                             #
+###############################################################################
+
+
+class ClumpsCatalogue(BaseCSiBORG):
     r"""
     Clumps catalogue, defined in the final snapshot.
 
@@ -320,8 +352,8 @@ class ClumpsCatalogue(BaseCatalogue):
     ----------
     sim : int
         IC realisation index.
-    paths : py:class`csiborgtools.read.CSiBORGPaths`
+    paths : py:class`csiborgtools.read.Paths`
-        CSiBORG paths object.
+        Paths object.
     maxdist : float, optional
         The maximum comoving distance of a halo. By default
         :math:`155.5 / 0.705 ~ \mathrm{Mpc}` with assumed :math:`h = 0.705`,
@@ -363,7 +395,7 @@ class ClumpsCatalogue(BaseCatalogue):
             names = ["x", "y", "z", "mass_cl", "totpartmass", "rho0", "r200c",
                      "r500c", "m200c", "m500c", "r200m", "m200m",
                      "vx", "vy", "vz"]
-            self._data = self.box.convert_from_boxunits(self._data, names)
+            self._data = self.box.convert_from_box(self._data, names)
             if maxdist is not None:
                 dist = numpy.sqrt(self._data["x"]**2 + self._data["y"]**2
                                   + self._data["z"]**2)
@@ -383,7 +415,12 @@ class ClumpsCatalogue(BaseCatalogue):
         return self["index"] == self["parent"]
 
 
-class HaloCatalogue(BaseCatalogue):
+###############################################################################
+#                        CSiBORG halo catalogue                               #
+###############################################################################
+
+
+class HaloCatalogue(BaseCSiBORG):
     r"""
     Halo catalogue, i.e. parent halos with summed substructure, defined in the
     final snapshot.
@@ -392,8 +429,8 @@ class HaloCatalogue(BaseCatalogue):
     ----------
     nsim : int
         IC realisation index.
-    paths : py:class`csiborgtools.read.CSiBORGPaths`
+    paths : py:class`csiborgtools.read.Paths`
-        CSiBORG paths object.
+        Paths object.
     maxdist : float, optional
         The maximum comoving distance of a halo. By default
         :math:`155.5 / 0.705 ~ \mathrm{Mpc}` with assumed :math:`h = 0.705`,
@@ -455,11 +492,11 @@ class HaloCatalogue(BaseCatalogue):
             names = ["x", "y", "z", "M", "totpartmass", "rho0", "r200c",
                      "r500c", "m200c", "m500c", "r200m", "m200m",
                      "vx", "vy", "vz"]
-            self._data = self.box.convert_from_boxunits(self._data, names)
+            self._data = self.box.convert_from_box(self._data, names)
 
             if load_initial:
                 names = ["x0", "y0", "z0", "lagpatch"]
-                self._data = self.box.convert_from_boxunits(self._data, names)
+                self._data = self.box.convert_from_box(self._data, names)
 
             if maxdist is not None:
                 dist = numpy.sqrt(self._data["x"]**2 + self._data["y"]**2
@@ -480,3 +517,112 @@ class HaloCatalogue(BaseCatalogue):
         if self._clumps_cat is None:
             raise ValueError("`clumps_cat` is not loaded.")
         return self._clumps_cat
+
+
+###############################################################################
+#                         Quijote halo catalogue                              #
+###############################################################################
+
+
+class QuijoteHaloCatalogue(BaseCatalogue):
+    """
+    Quijote halo catalogue.
+
+    Parameters
+    ----------
+    nsim : int
+        IC realisation index.
+    paths : py:class`csiborgtools.read.Paths`
+        Paths object.
+    nsnap : int
+        Snapshot index.
+    origin : len-3 tuple, optional
+        Where to place the origin of the box. By default the centre of the box.
+        In units of :math:`cMpc`.
+    maxdist : float, optional
+        The maximum comoving distance of a halo in the new reference frame, in
+        units of :math:`cMpc`.
+    minmass : len-2 tuple
+        Minimum mass. The first element is the catalogue key and the second is
+        the value.
+    rawdata : bool, optional
+        Whether to return the raw data. In this case applies no cuts and
+        transformations.
+    **kwargs : dict
+        Keyword arguments for backward compatibility.
+    """
+    _nsnap = None
+
+    def __init__(self, nsim, paths, nsnap,
+                 origin=[500 / 0.6711, 500 / 0.6711, 500 / 0.6711],
+                 maxdist=None, minmass=("group_mass", 1e12), rawdata=False,
+                 **kwargs):
+        self.paths = paths
+        self.nsnap = nsnap
+        fpath = join(self.paths.quijote_dir, "halos", str(nsim))
+        fof = FoF_catalog(fpath, self.nsnap, long_ids=False, swap=False,
+                          SFR=False, read_IDs=False)
+
+        cols = [("x", numpy.float32), ("y", numpy.float32),
+                ("z", numpy.float32), ("vx", numpy.float32),
+                ("vy", numpy.float32), ("vz", numpy.float32),
+                ("group_mass", numpy.float32), ("npart", numpy.int32)]
+        data = cols_to_structured(fof.GroupLen.size, cols)
+
+        pos = fof.GroupPos / 1e3 / self.box.h
+        for i in range(3):
+            pos -= origin[i]
+        vel = fof.GroupVel * (1 + self.redshift)
+        for i, p in enumerate(["x", "y", "z"]):
+            data[p] = pos[:, i]
+            data["v" + p] = vel[:, i]
+        data["group_mass"] = fof.GroupMass * 1e10 / self.box.h
+        data["npart"] = fof.GroupLen
+
+        if not rawdata:
+            if maxdist is not None:
+                pos = numpy.vstack([data["x"], data["y"], data["z"]]).T
+                data = data[numpy.linalg.norm(pos, axis=1) < maxdist]
+            if minmass is not None:
+                data = data[data[minmass[0]] > minmass[1]]
+
+        self._data = data
+
+    @property
+    def nsnap(self):
+        """
+        Snapshot number.
+
+        Returns
+        -------
+        nsnap : int
+        """
+        return self._nsnap
+
+    @nsnap.setter
+    def nsnap(self, nsnap):
+        assert nsnap in [0, 1, 2, 3, 4]
+        self._nsnap = nsnap
+
+    @property
+    def redshift(self):
+        """
+        Redshift of the snapshot.
+
+        Returns
+        -------
+        redshift : float
+        """
+        z_dict = {4: 0.0, 3: 0.5, 2: 1.0, 1: 2.0, 0: 3.0}
+        return z_dict[self.nsnap]
+
+    @property
+    def box(self):
+        """
+        Quijote box object.
+
+        Returns
+        -------
+        box : instance of :py:class:`csiborgtools.units.BaseBox`
+        """
+        return QuijoteBox(self.nsnap)

csiborgtools/read/knn_summary.py

@@ -24,7 +24,7 @@ class kNNCDFReader:
 
     Parameters
     ----------
-    paths : py:class`csiborgtools.read.CSiBORGPaths`
+    paths : py:class`csiborgtools.read.Paths`
     """
     _paths = None
 
@@ -38,13 +38,12 @@ class kNNCDFReader:
 
         Parameters
         ----------
-        paths : py:class`csiborgtools.read.CSiBORGPaths`
+        paths : py:class`csiborgtools.read.Paths`
         """
         return self._paths
 
     @paths.setter
     def paths(self, paths):
-        # assert isinstance(paths, CSiBORGPaths)  # REMOVE
         self._paths = paths
 
     def read(self, run, kind, rmin=None, rmax=None, to_clip=True):

csiborgtools/read/overlap_summary.py

@@ -36,7 +36,7 @@ class PairOverlap:
         Halo catalogue corresponding to the reference simulation.
     catx : :py:class:`csiborgtools.read.HaloCatalogue`
         Halo catalogue corresponding to the cross simulation.
-    paths : py:class`csiborgtools.read.CSiBORGPaths`
+    paths : py:class`csiborgtools.read.Paths`
         CSiBORG paths object.
     """
     _cat0 = None
@@ -59,7 +59,7 @@ class PairOverlap:
             Halo catalogue corresponding to the reference simulation.
         catx : :py:class:`csiborgtools.read.HaloCatalogue`
             Halo catalogue corresponding to the cross simulation.
-        paths : py:class`csiborgtools.read.CSiBORGPaths`
+        paths : py:class`csiborgtools.read.Paths`
             CSiBORG paths object.
 
         Returns
@@ -268,8 +268,8 @@ class PairOverlap:
         assert (norm_kind is None
                 or norm_kind in ("r200c", "ref_patch", "sum_patch"))
         # Get positions either in the initial or final snapshot
-        pos0 = self.cat0().position(in_initial)
+        pos0 = self.cat0().position(in_initial=in_initial)
-        posx = self.catx().position(in_initial)
+        posx = self.catx().position(in_initial=in_initial)
 
         # Get the normalisation array if applicable
         if norm_kind == "r200c":
@@ -474,7 +474,7 @@ class NPairsOverlap:
         Single reference simulation halo catalogue.
     catxs : list of :py:class:`csiborgtools.read.HaloCatalogue`
         List of cross simulation halo catalogues.
-    paths : py:class`csiborgtools.read.CSiBORGPaths`
+    paths : py:class`csiborgtools.read.Paths`
         CSiBORG paths object.
     verbose : bool, optional
         Verbosity flag for loading the overlap objects.

csiborgtools/read/paths.py

@@ -21,29 +21,32 @@ from warnings import warn
 import numpy
 
 
-class CSiBORGPaths:
+class Paths:
     """
-    Paths manager for CSiBORG IC realisations.
+    Paths manager for CSiBORG and Quijote simulations.
 
-     Parameters
+    Parameters
-     ----------
+    ----------
-     srcdir : str
+    srcdir : str, optional
-         Path to the folder where the RAMSES outputs are stored.
+        Path to the folder where the RAMSES outputs are stored.
-     postdir: str
+    postdir: str, optional
-         Path to the folder where post-processed files are stored.
+        Path to the folder where post-processed files are stored.
+    quiote_dir : str, optional
+        Path to the folder where Quijote simulations are stored.
     """
-
     _srcdir = None
     _postdir = None
+    _quijote_dir = None
 
-    def __init__(self, srcdir=None, postdir=None):
+    def __init__(self, srcdir=None, postdir=None, quijote_dir=None):
         self.srcdir = srcdir
         self.postdir = postdir
+        self.quijote_dir = quijote_dir
 
     @staticmethod
     def _check_directory(path):
         if not isdir(path):
-            raise IOError("Invalid directory `{}`!".format(path))
+            raise IOError(f"Invalid directory `{path}`!")
 
     @property
     def srcdir(self):
@@ -65,6 +68,26 @@ class CSiBORGPaths:
         self._check_directory(path)
         self._srcdir = path
 
+    @property
+    def quijote_dir(self):
+        """
+        Path to the folder where Quijote simulations are stored.
+
+        Returns
+        -------
+        path : str
+        """
+        if self._quijote_dir is None:
+            raise ValueError("`quijote_dir` is not set!")
+        return self._quijote_dir
+
+    @quijote_dir.setter
+    def quijote_dir(self, path):
+        if path is None:
+            return
+        self._check_directory(path)
+        self._quijote_dir = path
+
     @property
     def postdir(self):
         """

csiborgtools/read/readsim.py

@@ -22,7 +22,7 @@ import numpy
 from scipy.io import FortranFile
 from tqdm import tqdm, trange
 
-from .paths import CSiBORGPaths
+from .paths import Paths
 from .utils import cols_to_structured
 
 ###############################################################################
@@ -36,7 +36,7 @@ class ParticleReader:
 
     Parameters
     ----------
-    paths : py:class`csiborgtools.read.CSiBORGPaths`
+    paths : py:class`csiborgtools.read.Paths`
     """
     _paths = None
 
@@ -50,13 +50,13 @@ class ParticleReader:
 
         Parameters
         ----------
-        paths : py:class`csiborgtools.read.CSiBORGPaths`
+        paths : py:class`csiborgtools.read.Paths`
         """
         return self._paths
 
     @paths.setter
     def paths(self, paths):
-        assert isinstance(paths, CSiBORGPaths)
+        assert isinstance(paths, Paths)
         self._paths = paths
 
     def read_info(self, nsnap, nsim):
@@ -396,17 +396,17 @@ class ParticleReader:
 
 class MmainReader:
     """
-    Object to generate the summed substructure catalogue.
+    Object to generate the summed substructure CSiBORG PHEW catalogue.
 
     Parameters
     ----------
-    paths : :py:class:`csiborgtools.read.CSiBORGPaths`
+    paths : :py:class:`csiborgtools.read.Paths`
         Paths objects.
     """
     _paths = None
 
     def __init__(self, paths):
-        assert isinstance(paths, CSiBORGPaths)
+        assert isinstance(paths, Paths)
         self._paths = paths
 
     @property

csiborgtools/read/tpcf_summary.py

@@ -16,7 +16,7 @@
 import joblib
 import numpy
 
-from .paths import CSiBORGPaths
+from .paths import Paths
 
 
 class TPCFReader:
@@ -25,7 +25,7 @@ class TPCFReader:
 
     Parameters
     ----------
-    paths : py:class`csiborgtools.read.CSiBORGPaths`
+    paths : py:class`csiborgtools.read.Paths`
     """
     _paths = None
 
@@ -39,13 +39,13 @@ class TPCFReader:
 
         Parameters
         ----------
-        paths : py:class`csiborgtools.read.CSiBORGPaths`
+        paths : py:class`csiborgtools.read.Paths`
         """
         return self._paths
 
     @paths.setter
     def paths(self, paths):
-        assert isinstance(paths, CSiBORGPaths)
+        assert isinstance(paths, Paths)
         self._paths = paths
 
     def read(self, run):

csiborgtools/read/utils.py

@@ -93,7 +93,7 @@ def real2redshift(pos, vel, origin, box, in_box_units, periodic_wrap=True,
         Cartesian velocity components.
     origin : 1-dimensional array `(3,)`
         Origin of the coordinate system in the `pos` reference frame.
-    box : py:class:`csiborg.read.BoxUnits`
+    box : py:class:`csiborg.read.CSiBORGBox`
         Box units.
     in_box_units: bool
         Whether `pos` and `vel` are in box units. If not, position is assumed

notebooks/perseus.ipynb

@@ -50,7 +50,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "paths = csiborgtools.read.CSiBORGPaths(**csiborgtools.paths_glamdring)\n",
+    "paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)\n",
     "nsims = paths.get_ics(False)"
    ]
   },

notebooks/playground_field.ipynb

@@ -97,9 +97,9 @@
    "source": [
     "grid = 256\n",
     "length = 677.05\n",
-    "paths = csiborgtools.read.CSiBORGPaths()\n",
+    "paths = csiborgtools.read.Paths()\n",
     "paths.set_info(7444, paths.get_maximum_snapshot(7444))\n",
-    "box = csiborgtools.units.BoxUnits(paths)\n",
+    "box = csiborgtools.units.CSiBORGBox(paths)\n",
     "MAS = \"CIC\"\n",
     "field = csiborgtools.field.DensityField(particles, length, box, MAS)"
    ]

notebooks/plot_correlation.ipynb

@@ -52,7 +52,7 @@
    },
    "outputs": [],
    "source": [
-    "paths = csiborgtools.read.CSiBORGPaths()\n",
+    "paths = csiborgtools.read.Paths()\n",
     "hw = 0.05\n",
     "darkskyps = np.genfromtxt(\"../plots/darksky.csv\", delimiter=\", \")\n",
     "\n",
@@ -3699,7 +3699,7 @@
    },
    "outputs": [],
    "source": [
-    "paths = csiborgtools.read.CSiBORGPaths()\n",
+    "paths = csiborgtools.read.Paths()\n",
     "cat = csiborgtools.read.CombinedCatalogue(paths, min_m500=1e13, max_dist=210)"
    ]
   },

notebooks/plot_mass_function.ipynb

@@ -55,7 +55,7 @@
     "\n",
     "data = data[(data[\"npart\"] > 100) & np.isfinite(data[\"m200\"])]\n",
     "\n",
-    "boxunits = csiborgtools.units.BoxUnits(Nsnap, simpath)"
+    "CSiBORGBox = csiborgtools.units.CSiBORGBox(Nsnap, simpath)"
    ]
   },
   {
@@ -96,7 +96,7 @@
     }
    ],
    "source": [
-    "boxunits.box2kpc(0.21) * 1e-3"
+    "CSiBORGBox.box2kpc(0.21) * 1e-3"
    ]
   },
   {
@@ -1120,7 +1120,7 @@
     }
    ],
    "source": [
-    "mass = boxunits.box2solarmass(data[\"m200\"])\n",
+    "mass = CSiBORGBox.box2solarmass(data[\"m200\"])\n",
     "\n"
    ]
   },
@@ -1220,7 +1220,7 @@
    "source": [
     "# n = 584\n",
     "# xs = csiborgtools.fits.pick_single_clump(n, parts, part_clumps, clumps)\n",
-    "# halo = csiborgtools.fits.Clump.from_arrays(*xs, rhoc=boxunits.box_rhoc)\n",
+    "# halo = csiborgtools.fits.Clump.from_arrays(*xs, rhoc=CSiBORGBox.box_rhoc)\n",
     "# print(halo.Npart)"
    ]
   },

notebooks/test_mmain.ipynb

@@ -69,7 +69,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "paths = csiborgtools.read.CSiBORGPaths(**csiborgtools.paths_glamdring)\n",
+    "paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)\n",
     "# cat = csiborgtools.read.ClumpsCatalogue(7444, paths)"
    ]
   },

scripts/cluster_crosspk.py

@@ -48,8 +48,8 @@ rank = comm.Get_rank()
 nproc = comm.Get_size()
 MAS = "CIC"  # mass asignment scheme
 
-paths = csiborgtools.read.CSiBORGPaths(**csiborgtools.paths_glamdring)
+paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
-box = csiborgtools.read.BoxUnits(paths)
+box = csiborgtools.read.CSiBORGBox(paths)
 reader = csiborgtools.read.ParticleReader(paths)
 ics = paths.get_ics()
 nsims = len(ics)

scripts/cluster_knn_auto.py

@@ -49,7 +49,7 @@ with open("../scripts/knn_auto.yml", "r") as file:
 
 Rmax = 155 / 0.705  # Mpc (h = 0.705) high resolution region radius
 totvol = 4 * numpy.pi * Rmax**3 / 3
-paths = csiborgtools.read.CSiBORGPaths(**csiborgtools.paths_glamdring)
+paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
 ics = paths.get_ics()
 knncdf = csiborgtools.clustering.kNN_CDF()
 

scripts/cluster_knn_cross.py

@@ -48,7 +48,7 @@ with open("../scripts/knn_cross.yml", "r") as file:
     config = yaml.safe_load(file)
 
 Rmax = 155 / 0.705  # Mpc (h = 0.705) high resolution region radius
-paths = csiborgtools.read.CSiBORGPaths(**csiborgtools.paths_glamdring)
+paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
 ics = paths.get_ics()
 knncdf = csiborgtools.clustering.kNN_CDF()
 

scripts/cluster_tcpf_auto.py

@@ -47,7 +47,7 @@ with open("../scripts/tpcf_auto.yml", "r") as file:
     config = yaml.safe_load(file)
 
 Rmax = 155 / 0.705  # Mpc (h = 0.705) high resolution region radius
-paths = csiborgtools.read.CSiBORGPaths()
+paths = csiborgtools.read.Paths()
 ics = paths.get_ics()
 tpcf = csiborgtools.clustering.Mock2PCF()
 

scripts/field_derived.py

@@ -45,7 +45,7 @@ parser.add_argument("--grid", type=int, help="Grid resolution.")
 parser.add_argument("--in_rsp", type=lambda x: bool(strtobool(x)),
                     help="Calculate from the RSP density field?")
 args = parser.parse_args()
-paths = csiborgtools.read.CSiBORGPaths(**csiborgtools.paths_glamdring)
+paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
 
 if args.ics is None or args.ics[0] == -1:
     ics = paths.get_ics()
@@ -59,7 +59,7 @@ for i in csiborgtools.fits.split_jobs(len(ics), nproc)[rank]:
         print(f"{datetime.now()}: rank {rank} working on simulation {nsim}.",
               flush=True)
     nsnap = max(paths.get_snapshots(nsim))
-    box = csiborgtools.read.BoxUnits(nsnap, nsim, paths)
+    box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths)
     density_gen = csiborgtools.field.DensityField(box, args.MAS)
 
     rho = numpy.load(paths.field_path("density", args.MAS, args.grid, nsim,

scripts/field_fromparts.py

@@ -46,7 +46,7 @@ parser.add_argument("--grid", type=int, help="Grid resolution.")
 parser.add_argument("--in_rsp", type=lambda x: bool(strtobool(x)),
                     help="Calculate the density field in redshift space?")
 args = parser.parse_args()
-paths = csiborgtools.read.CSiBORGPaths(**csiborgtools.paths_glamdring)
+paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
 mpart = 1.1641532e-10  # Particle mass in CSiBORG simulations.
 
 if args.ics is None or args.ics[0] == -1:
@@ -61,7 +61,7 @@ for i in csiborgtools.fits.split_jobs(len(ics), nproc)[rank]:
           flush=True)
 
     nsnap = max(paths.get_snapshots(nsim))
-    box = csiborgtools.read.BoxUnits(nsnap, nsim, paths)
+    box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths)
     parts = csiborgtools.read.read_h5(paths.particles_path(nsim))["particles"]
 
     if args.kind == "density":

scripts/fit_halos.py

@@ -42,7 +42,7 @@ parser.add_argument("--kind", type=str, choices=["halos", "clumps"])
 parser.add_argument("--ics", type=int, nargs="+", default=None,
                     help="IC realisations. If `-1` processes all simulations.")
 args = parser.parse_args()
-paths = csiborgtools.read.CSiBORGPaths(**csiborgtools.paths_glamdring)
+paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
 partreader = csiborgtools.read.ParticleReader(paths)
 nfwpost = csiborgtools.fits.NFWPosterior()
 
@@ -105,7 +105,7 @@ for nsim in [ics[i] for i in jobs]:
     print(f"{datetime.now()}: rank {rank} calculating simulation `{nsim}`.",
           flush=True)
     nsnap = max(paths.get_snapshots(nsim))
-    box = csiborgtools.read.BoxUnits(nsnap, nsim, paths)
+    box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths)
 
     # Particle archive
     f = csiborgtools.read.read_h5(paths.particles_path(nsim))

scripts/fit_init.py

@@ -22,7 +22,6 @@ from datetime import datetime
 
 import numpy
 from mpi4py import MPI
-
 from tqdm import tqdm
 
 try:
@@ -45,7 +44,7 @@ parser = ArgumentParser()
 parser.add_argument("--ics", type=int, nargs="+", default=None,
                     help="IC realisations. If `-1` processes all simulations.")
 args = parser.parse_args()
-paths = csiborgtools.read.CSiBORGPaths(**csiborgtools.paths_glamdring)
+paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
 partreader = csiborgtools.read.ParticleReader(paths)
 
 if args.ics is None or args.ics[0] == -1:

scripts/fit_profiles.py

@@ -45,7 +45,7 @@ nproc = comm.Get_size()
 if nproc > 1:
     raise NotImplementedError("MPI is not implemented implemented yet.")
 
-paths = csiborgtools.read.CSiBORGPaths(**csiborgtools.paths_glamdring)
+paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
 cols_collect = [("r", numpy.float32), ("M", numpy.float32)]
 if args.ics is None or args.ics == -1:
     nsims = paths.get_ics()
@@ -59,7 +59,7 @@ for i, nsim in enumerate(nsims):
         now = datetime.now()
         print(f"{now}: calculating {i}th simulation `{nsim}`.", flush=True)
     nsnap = max(paths.get_snapshots(nsim))
-    box = csiborgtools.read.BoxUnits(nsnap, nsim, paths)
+    box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths)
 
     f = csiborgtools.read.read_h5(paths.particles_path(nsim))
     particles = f["particles"]

scripts/match_all.py

@@ -54,7 +54,7 @@ def get_combs():
     Get the list of all pairs of simulations, then permute them with a known
     seed to minimise loading the same files simultaneously.
     """
-    paths = csiborgtools.read.CSiBORGPaths(**csiborgtools.paths_glamdring)
+    paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
     ics = paths.get_ics()
     combs = list(combinations(ics, 2))
     Random(42).shuffle(combs)

scripts/match_singlematch.py

@@ -32,7 +32,7 @@ except ModuleNotFoundError:
 def pair_match(nsim0, nsimx, sigma, smoothen, verbose):
     from csiborgtools.read import HaloCatalogue, read_h5
 
-    paths = csiborgtools.read.CSiBORGPaths(**csiborgtools.paths_glamdring)
+    paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
     smooth_kwargs = {"sigma": sigma, "mode": "constant", "cval": 0.0}
     overlapper = csiborgtools.match.ParticleOverlap()
     matcher = csiborgtools.match.RealisationsMatcher()

scripts/pre_dumppart.py

@@ -49,7 +49,7 @@ parser.add_argument("--ics", type=int, nargs="+", default=None,
 args = parser.parse_args()
 
 verbose = nproc == 1
-paths = csiborgtools.read.CSiBORGPaths(**csiborgtools.paths_glamdring)
+paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
 partreader = csiborgtools.read.ParticleReader(paths)
 # Keep "ID" as the last column!
 pars_extract = ['x', 'y', 'z', 'vx', 'vy', 'vz', 'M', "ID"]

scripts/pre_mmain.py

@@ -33,7 +33,7 @@ comm = MPI.COMM_WORLD
 rank = comm.Get_rank()
 nproc = comm.Get_size()
 
-paths = csiborgtools.read.CSiBORGPaths(**csiborgtools.paths_glamdring)
+paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
 mmain_reader = csiborgtools.read.MmainReader(paths)
 
 

scripts/pre_sortinit.py

@@ -18,9 +18,9 @@ snapshot ordering, which is sorted by the clump IDs.
 """
 from argparse import ArgumentParser
 from datetime import datetime
+from gc import collect
 
 import h5py
-from gc import collect
 import numpy
 from mpi4py import MPI
 
@@ -44,7 +44,7 @@ parser = ArgumentParser()
 parser.add_argument("--ics", type=int, nargs="+", default=None,
                     help="IC realisations. If `-1` processes all simulations.")
 args = parser.parse_args()
-paths = csiborgtools.read.CSiBORGPaths(**csiborgtools.paths_glamdring)
+paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
 partreader = csiborgtools.read.ParticleReader(paths)
 # NOTE: ID has to be the last column.
 pars_extract = ["x", "y", "z", "M", "ID"]