Extend enclosed density to all CB2 (#126)

* Add borg2 flag

* add borg2_all

* Undo a comment?

* ADd borg2

* Update notebook

* Update nb

* Add external halo catalogue

* Add MDPL2

* Add BORG2 density profile plot

* Add more params

* Add comments

csiborgtools/params.py

@@ -92,6 +92,7 @@ def simname2boxsize(simname):
          "csiborg2_random": 676.6,
          "borg1": 677.7,
          "borg2": 676.6,
+         "borg2_all": 676.6,
          "quijote": 1000.,
          "TNG300-1": 205.,
          "Carrick2015": 400.,
@@ -123,6 +124,8 @@ def simname2Omega_m(simname):
          "csiborg2_random": 0.3111,
          "csiborg2_varysmall": 0.3111,
          "borg1": 0.307,
+         "borg2": 0.3111,
+         "borg2_all": 0.3111,
          "Carrick2015": 0.3,
          }
 
@@ -144,6 +147,7 @@ paths_glamdring = {
     "borg1_dir": "/mnt/users/hdesmond/BORG_final",
     "borg2_dir": "/mnt/extraspace/rstiskalek/BORG_STOPYRA_2023",
     "tng300_1_dir": "/mnt/extraspace/rstiskalek/TNG300-1/",
+    "aux_cat_dir": "/mnt/extraspace/rstiskalek/catalogs",
     }
 
 

csiborgtools/read/__init__.py

@@ -14,7 +14,8 @@
 # 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 from .catalogue import (CSiBORG1Catalogue, CSiBORG2Catalogue,                   # noqa
                         CSiBORG2SUBFINDCatalogue,                               # noqa
-                        CSiBORG2MergerTreeReader, QuijoteCatalogue)             # noqa
+                        CSiBORG2MergerTreeReader, QuijoteCatalogue,             # noqa
+                        MDPL2Catalogue)                                         # noqa
 from .snapshot import (CSiBORG1Snapshot, CSiBORG2Snapshot, QuijoteSnapshot,     # noqa
                        CSiBORG1Field, CSiBORG2Field, QuijoteField, BORG2Field,  # noqa
                        BORG1Field, TNG300_1Field)                               # noqa

csiborgtools/read/catalogue.py

@@ -582,6 +582,29 @@ class BaseCatalogue(ABC):
         """
         return self._properties + self._custom_keys
 
+    def pick_fiducial_observer(self, n, rmax):
+        r"""
+        Select a new fiducial observer in the box.
+
+        Parameters
+        ----------
+        n : int
+            Fiducial observer index.
+        rmax : float
+            Max. distance from the fiducial obs. in :math:`\mathrm{cMpc} / h`.
+        """
+        self.clear_cache()
+        print(fiducial_observers(self.boxsize, rmax))
+        self.observer_location = fiducial_observers(self.boxsize, rmax)[n]
+        self.observer_velocity = None
+
+        if self._bounds is None:
+            bounds = {"dist": (0, rmax)}
+        else:
+            bounds = {**self._bounds, "dist": (0, rmax)}
+
+        self._make_mask(bounds)
+
     def __getitem__(self, key):
         # For internal calls we don't want to load the filtered data and use
         # the __ prefixed keys. The internal calls are not being cached.
@@ -1366,27 +1389,83 @@ class QuijoteCatalogue(BaseCatalogue):
         fpath = self.paths.initial_lagpatch(self.nsim, self.simname)
         return numpy.load(fpath)["lagpatch_size"]
 
-    def pick_fiducial_observer(self, n, rmax):
+
+###############################################################################
+#                     External halo catalogues                                #
+###############################################################################
+
+class MDPL2Catalogue(BaseCatalogue):
     r"""
-        Select a new fiducial observer in the box.
+    XXX
 
     Parameters
     ----------
-        n : int
+    nsim : int
-            Fiducial observer index.
+        IC realisation index.
-        rmax : float
+    paths : py:class`csiborgtools.read.Paths`, optional
-            Max. distance from the fiducial obs. in :math:`\mathrm{cMpc} / h`.
+        Paths object.
+    snapshot : subclass of py:class:`BaseSnapshot`, optional
+        Snapshot object corresponding to the catalogue.
+    bounds : dict
+        Parameter bounds; keys as parameter names, values as (min, max)
+        tuples. Use `dist` for radial distance, `None` for no bound.
+    observer_velocity : array, optional
+        Observer's velocity in :math:`\mathrm{km} / \mathrm{s}`.
+    cache_maxsize : int, optional
+        Maximum number of cached arrays.
     """
-        self.clear_cache()
+    def __init__(self, paths=None, bounds=None, cache_maxsize=64):
-        self.observer_location = fiducial_observers(self.box.boxsize, rmax)[n]
+        boxsize = 1000.
-        self.observer_velocity = None
+        super().__init__()
+        x0 = boxsize / 2
+        super().init_with_snapshot(
+            "MDPL2", 0, 125, paths, None, bounds, boxsize, [x0, x0, x0], None,
+            False, cache_maxsize)
 
-        if self._bounds is None:
+        self._custom_keys = []
-            bounds = {"dist": (0, rmax)}
+        self._bounds = bounds
-        else:
-            bounds = {**self._bounds, "dist": (0, rmax)}
 
-        self._make_mask(bounds)
+    def _read_fof_catalogue(self, kind):
+        fpath = self.paths.external_halo_catalogue(self.simname)
+
+        with File(fpath, 'r') as f:
+            if kind == "index":
+                return numpy.arange(len(f["x"]))
+
+            if kind not in f.keys():
+                raise ValueError(f"FoF catalogue key '{kind}' not available. Available keys are: {list(f.keys())}")  # noqa
+            out = f[kind][...]
+        return out
+
+    @property
+    def coordinates(self):
+        return numpy.vstack(
+            [self._read_fof_catalogue(key) for key in ["x", "y", "z"]]).T
+
+    @property
+    def velocities(self):
+        return numpy.vstack(
+            [self._read_fof_catalogue(key) for key in ["vx", "vy", "vz"]]).T
+
+    @property
+    def totmass(self):
+        return self._read_fof_catalogue("mass")
+
+    @property
+    def npart(self):
+        raise RuntimeError("Number of particles is not available.")
+
+    @property
+    def index(self):
+        return self._read_fof_catalogue("index")
+
+    @property
+    def lagpatch_coordinates(self):
+        raise RuntimeError("Lagrangian patch information is not available")
+
+    @property
+    def lagpatch_radius(self):
+        raise RuntimeError("Lagrangian patch information is not available")
 
 
 ###############################################################################

csiborgtools/read/paths.py

@@ -59,6 +59,8 @@ class Paths:
         Path to the BORG2 simulation directory.
     tng300_1_dir : str
         Path to the TNG300-1 simulation directory.
+    aux_cat_dir : str
+        Path to the directory containing auxiliary catalogues.
     """
     def __init__(self,
                  csiborg1_srcdir,
@@ -69,7 +71,8 @@ class Paths:
                  quijote_dir,
                  borg1_dir,
                  borg2_dir,
-                 tng300_1_dir
+                 tng300_1_dir,
+                 aux_cat_dir
                  ):
         self.csiborg1_srcdir = csiborg1_srcdir
         self.csiborg2_main_srcdir = csiborg2_main_srcdir
@@ -80,6 +83,7 @@ class Paths:
         self.borg2_dir = borg2_dir
         self.tng300_1_dir = tng300_1_dir
         self.postdir = postdir
+        self.aux_cat_dir = aux_cat_dir
 
     def get_ics(self, simname):
         """
@@ -100,6 +104,9 @@ class Paths:
         elif simname == "csiborg2_main" or simname == "borg2":
             files = glob(join(self.csiborg2_main_srcdir, "chain_*"))
             files = [int(search(r'chain_(\d+)', f).group(1)) for f in files]
+        elif simname == "borg2_all":
+            files = glob(join(self.borg2_dir, "mcmc_*"))
+            files = [int(search(r'mcmc_(\d+)', f).group(1)) for f in files]
         elif simname == "csiborg2_random":
             files = glob(join(self.csiborg2_random_srcdir, "chain_*"))
             files = [int(search(r'chain_(\d+)', f).group(1)) for f in files]
@@ -246,6 +253,24 @@ class Paths:
         else:
             raise ValueError(f"Unknown simulation name `{simname}`.")
 
+    def external_halo_catalogue(self, name):
+        """
+        Path to an external halo catalogue.
+
+        Parameters
+        ----------
+        name : str
+            Catalogue name.
+
+        Returns
+        -------
+        str
+        """
+        if name == "MDPL2":
+            return join(self.aux_cat_dir, "MDPL2_FOF_125.hdf5")
+        else:
+            raise ValueError(f"Unknown external FOF catalogue `{name}`.")
+
     def initial_lagpatch(self, nsim, simname):
         """
         Path to the Lagrangain patch information of a simulation for halos

csiborgtools/read/snapshot.py

@@ -969,7 +969,7 @@ class BORG2Field(BaseField):
         rho_mean = omega0 * 277.53662724583074  # h^2 Msun / kpc^3
         field += 1
         field *= rho_mean
-        # return field
+        return field
 
     def velocity_field(self, MAS, grid):
         raise RuntimeError("The velocity field is not available.")

csiborgtools/utils.py

@@ -35,6 +35,19 @@ def center_of_mass(particle_positions, particles_mass, boxsize):
     """
     Calculate the center of mass of a halo while assuming periodic boundary
     conditions of a cubical box.
+
+    Parameters
+    ----------
+    particle_positions : 2-dimensional array of shape `(nparticles, 3)`
+        Particle positions in the box.
+    particles_mass : 1-dimensional array of shape `(nparticles,)`
+        Particle masses.
+    boxsize : float
+        Box size.
+
+    Returns
+    -------
+    1-dimensional array of shape `(3,)`
     """
     cm = np.zeros(3, dtype=particle_positions.dtype)
     totmass = sum(particles_mass)
@@ -60,28 +73,46 @@ def periodic_distance(points, reference_point, boxsize):
     """
     Compute the 3D distance between multiple points and a reference point using
     periodic boundary conditions.
+
+    Parameters
+    ----------
+    points : 2-dimensional array of shape `(npoints, 3)`
+        Points to calculate the distance from.
+    reference_point : 1-dimensional array of shape `(3,)`
+        Reference point.
+    boxsize : float
+        Box size.
+
+    Returns
+    -------
+    1-dimensional array of shape `(npoints,)`
     """
     npoints = len(points)
-    half_box = boxsize / 2
 
     dist = np.zeros(npoints, dtype=points.dtype)
     for i in range(npoints):
-        for j in range(3):
+        dist[i] = periodic_distance_two_points(
-            dist_1d = abs(points[i, j] - reference_point[j])
+            points[i], reference_point, boxsize)
-
-            if dist_1d > (half_box):
-                dist_1d = boxsize - dist_1d
-
-            dist[i] += dist_1d**2
-
-        dist[i] = dist[i]**0.5
 
     return dist
 
 
 @jit(nopython=True, fastmath=True, boundscheck=False)
 def periodic_distance_two_points(p1, p2, boxsize):
-    """Compute the 3D distance between two points in a periodic box."""
+    """
+    Compute the 3D distance between two points in a periodic box.
+
+    Parameters
+    ----------
+    p1, p2 : 1-dimensional array of shape `(3,)`
+        Points to calculate the distance between.
+    boxsize : float
+        Box size.
+
+    Returns
+    -------
+    float
+    """
     half_box = boxsize / 2
 
     dist = 0
@@ -98,7 +129,20 @@ def periodic_distance_two_points(p1, p2, boxsize):
 
 @jit(nopython=True, boundscheck=False)
 def periodic_wrap_grid(pos, boxsize=1):
-    """Wrap positions in a periodic box."""
+    """
+    Wrap positions in a periodic box. Overwrites the input array.
+
+    Parameters
+    ----------
+    pos : 2-dimensional array of shape `(npoints, 3)`
+        Positions to wrap.
+    boxsize : float, optional
+        Box size.
+
+    Returns
+    -------
+    2-dimensional array of shape `(npoints, 3)`
+    """
     for n in range(pos.shape[0]):
         for i in range(3):
             if pos[n, i] > boxsize:
@@ -113,6 +157,15 @@ def periodic_wrap_grid(pos, boxsize=1):
 def delta2ncells(field):
     """
     Calculate the number of cells in `field` that are non-zero.
+
+    Parameters
+    ----------
+    field : 3-dimensional array of shape `(nx, ny, nz)`
+        Field to calculate the number of non-zero cells.
+
+    Returns
+    -------
+    int
     """
     tot = 0
     imax, jmax, kmax = field.shape
@@ -124,20 +177,44 @@ def delta2ncells(field):
     return tot
 
 
-def cartesian_to_radec(X):
+def cartesian_to_radec(X, return_degrees=True, origin=[0., 0., 0.]):
     """
     Calculate the radial distance, RA [0, 360) deg and dec [-90, 90] deg.
+
+    Parameters
+    ----------
+    X : 2-dimensional array of shape `(npoints, 3)`
+        Cartesian coordinates.
+    return_degrees : bool, optional
+        Whether to return the angles in degrees.
+    origin : 1-dimensional array of shape `(3,)`, optional
+        Origin of the coordinate system.
+
+    Returns
+    -------
+    out : 2-dimensional array of shape `(npoints, 3)`
+        Spherical coordinates: distance, RA and dec.
     """
     x, y, z = X[:, 0], X[:, 1], X[:, 2]
 
+    x -= origin[0]
+    y -= origin[1]
+    z -= origin[2]
+
     dist = np.linalg.norm(X, axis=1)
     dec = np.arcsin(z / dist)
     ra = np.arctan2(y, x)
     ra[ra < 0] += 2 * np.pi
 
+    if return_degrees:
         ra *= 180 / np.pi
         dec *= 180 / np.pi
 
+    # Place the origin back
+    x += origin[0]
+    y += origin[1]
+    z += origin[2]
+
     return np.vstack([dist, ra, dec]).T
 
 
@@ -145,6 +222,15 @@ def radec_to_cartesian(X):
     """
     Calculate Cartesian coordinates from radial distance, RA [0, 360) deg  and
     dec [-90, 90] deg.
+
+    Parameters
+    ----------
+    X : 2-dimensional array of shape `(npoints, 3)`
+        Spherical coordinates: distance, RA and dec.
+
+    Returns
+    -------
+    2-dimensional array of shape `(npoints, 3)`
     """
     dist, ra, dec = X[:, 0], X[:, 1], X[:, 2]
 
@@ -175,23 +261,43 @@ def radec_to_galactic(ra, dec):
 
 
 @jit(nopython=True, fastmath=True, boundscheck=False)
-def great_circle_distance(x1, x2):
+def great_circle_distance(x1, x2, in_degrees=True):
     """
     Great circle distance between two points on a sphere, defined by RA and
     dec, both in degrees.
+
+    Parameters
+    ----------
+    x1, x2 : 1-dimensional arrays of shape `(2,)`
+        RA and dec in degrees.
+    in_degrees : bool, optional
+        Whether the input is in degrees.
+
+    Returns
+    -------
+    float
     """
     ra1, dec1 = x1
     ra2, dec2 = x2
 
+    if in_degrees:
         ra1 *= np.pi / 180
         dec1 *= np.pi / 180
         ra2 *= np.pi / 180
         dec2 *= np.pi / 180
 
-    return 180 / np.pi * np.arccos(
+    dist = np.arccos(np.sin(dec1) * np.sin(dec2)
-        np.sin(dec1) * np.sin(dec2)
+                     + np.cos(dec1) * np.cos(dec2) * np.cos(ra1 - ra2))
-        + np.cos(dec1) * np.cos(dec2) * np.cos(ra1 - ra2)
+
-        )
+    # Convert to degrees and ensure the inputs are unchanged.
+    if in_degrees:
+        dist *= 180 / np.pi
+        ra1 *= 180 / np.pi
+        dec1 *= 180 / np.pi
+        ra2 *= 180 / np.pi
+        dec2 *= 180 / np.pi
+
+    return dist
 
 
 def cosine_similarity(x, y):
@@ -316,6 +422,17 @@ def real2redshift(pos, vel, observer_location, observer_velocity, boxsize,
 def number_counts(x, bin_edges):
     """
     Calculate counts of samples in bins.
+
+    Parameters
+    ----------
+    x : 1-dimensional array
+        Samples to bin.
+    bin_edges : 1-dimensional array
+        Bin edges.
+
+    Returns
+    -------
+    1-dimensional array
     """
     out = np.full(bin_edges.size - 1, np.nan, dtype=np.float32)
     for i in range(bin_edges.size - 1):
@@ -326,6 +443,21 @@ def number_counts(x, bin_edges):
 def binned_statistic(x, y, left_edges, bin_width, statistic):
     """
     Calculate a binned statistic.
+
+    Parameters
+    ----------
+    x, y : 1-dimensional arrays
+        Values by which to bin and calculate the statistic on, respectively.
+    left_edges : 1-dimensional array
+        Left edges of the bins.
+    bin_width : float
+        Width of the bins.
+    statistic : callable
+        Function to calculate the statistic, must be `f(x)`.
+
+    Returns
+    -------
+    1-dimensional array
     """
     out = np.full(left_edges.size, np.nan, dtype=x.dtype)
 

scripts/field_bulk.py

@@ -148,7 +148,7 @@ def main_borg(args, folder):
         if args.simname == "borg1":
             reader = csiborgtools.read.BORG1Field(nsim)
             field = reader.density_field()
-        elif args.simname == "borg2":
+        elif args.simname == "borg2" or args.simname == "borg2_all":
             reader = csiborgtools.read.BORG2Field(nsim)
             field = reader.density_field()
         else:
@@ -204,7 +204,7 @@ def main_csiborg(args, folder):
 if __name__ == "__main__":
     parser = ArgumentParser()
     parser.add_argument("--simname", type=str, help="Simulation name.",
-                        choices=["csiborg1", "csiborg2_main", "csiborg2_varysmall", "csiborg2_random", "borg1", "borg2"])  # noqa
+                        choices=["csiborg1", "csiborg2_main", "csiborg2_varysmall", "csiborg2_random", "borg1", "borg2", "borg2_all"])  # noqa
     args = parser.parse_args()
 
     folder = "/mnt/extraspace/rstiskalek/csiborg_postprocessing/field_shells"

scripts/field_bulk.sh

@@ -1,9 +1,9 @@
 nthreads=1
-memory=40
+memory=12
 on_login=0
-queue="cmb"
+queue="berg"
 env="/mnt/zfsusers/rstiskalek/csiborgtools/venv_csiborg/bin/python"
-file="mass_enclosed.py"
+file="field_bulk.py"
 
 simname=${1}