Matching of observations (#127)

* Rename file

* Add indents

* Update imports

* Add counting

* Docs

* Add nb

* Rename nb

* Update nb

* Add PV processing

* Update nb

* Add Pantheon+groups

* Update submission scripts

* Add Pantheon+zSN

* Update nb

* Edit param

* Matchin SFI

* Update nb

* Fix path bug

* Add list of clusters

* Update imports

* Update imports

* Add cartesian & mass of clusters

* Add observation to halo matching

* Add nb

* Add inverse CDF

* Add import

* Update nb

* Add comments

csiborgtools/__init__.py

@@ -14,6 +14,7 @@
 # 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 from csiborgtools import clustering, field, flow, halo, match, read, summary    # noqa
 
+from .clusters import clusters                                                  # noqa
 from .utils import (center_of_mass, delta2ncells, number_counts,                # noqa
                     periodic_distance, periodic_distance_two_points,            # noqa
                     binned_statistic, cosine_similarity, fprint,                # noqa
@@ -58,18 +59,3 @@ class SDSSxALFALFA:
         survey = read.SDSS(fpath, h=1, sel_steps=sel_steps)
         survey.name = "SDSSxALFALFA"
         return survey
-
-
-###############################################################################
-#                              Clusters                                       #
-###############################################################################
-
-clusters = {"Virgo": read.ObservedCluster(RA=hms_to_degrees(12, 27),
-                                          dec=dms_to_degrees(12, 43),
-                                          dist=16.5 * 0.7,
-                                          name="Virgo"),
-            "Fornax": read.ObservedCluster(RA=hms_to_degrees(3, 38),
-                                           dec=dms_to_degrees(-35, 27),
-                                           dist=19 * 0.7,
-                                           name="Fornax"),
-            }

csiborgtools/clusters.py

@@ -0,0 +1,105 @@
+# Copyright (C) 2024 Richard Stiskalek
+# This program is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the
+# Free Software Foundation; either version 3 of the License, or (at your
+# option) any later version.
+#
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
+# Public License for more details.
+#
+# 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.
+"""
+Database of a few nearby observed clusters. Can be augmented with the list
+compiled in https://arxiv.org/abs/2402.01834 or some eROSITA clusters?
+"""
+from csiborgtools.read import ObservedCluster
+from .utils import hms_to_degrees, dms_to_degrees
+
+
+# https://arxiv.org/abs/astro-ph/0702510
+# https://arxiv.org/abs/2002.12820
+# https://en.wikipedia.org/wiki/Virgo_Cluster
+VIRGO = ObservedCluster(
+    RA=hms_to_degrees(12, 27), dec=dms_to_degrees(12, 43), dist=16.5 * 0.73,
+    mass=6.3e14 * 0.73, name="Virgo")
+
+# https://arxiv.org/abs/astro-ph/0702320
+# https://en.wikipedia.org/wiki/Fornax_Cluster
+FORNAX = ObservedCluster(
+    RA=hms_to_degrees(3, 35), dec=-35.7, dist=19.3 * 0.7,
+    mass=7e13 * 0.73, name="Fornax")
+
+# https://en.wikipedia.org/wiki/Coma_Cluster
+# https://arxiv.org/abs/2311.08603
+COMA = ObservedCluster(
+    RA=hms_to_degrees(12, 59, 48.7), dec=dms_to_degrees(27, 58, 50),
+    dist=102.975 * 0.705, mass=1.2e15 * 0.73, name="Coma")
+
+# https://en.wikipedia.org/wiki/Perseus_Cluster
+# https://ui.adsabs.harvard.edu/abs/2020MNRAS.494.1681A/abstract
+PERSEUS = ObservedCluster(
+    RA=hms_to_degrees(3, 18), dec=dms_to_degrees(41, 30),
+    dist=73.6 * 0.705, mass=1.2e15 * 0.7, name="Perseus")
+
+# https://en.wikipedia.org/wiki/Centaurus_Cluster
+# Not sure about the mass, couldn't find a good estimate. Some paper claimed
+# 3e13 Msun, but that seems a little low?
+CENTAURUS = ObservedCluster(
+    RA=hms_to_degrees(12, 48, 51.8), dec=dms_to_degrees(-41, 18, 21),
+    dist=52.4 * 0.705, mass=2e14 * 0.7, name="Centaurus")
+
+# https://en.wikipedia.org/wiki/Shapley_Supercluster
+# https://arxiv.org/abs/0805.0596
+SHAPLEY = ObservedCluster(
+    RA=hms_to_degrees(13, 25), dec=dms_to_degrees(-30),
+    dist=136, mass=1e16 * 0.7, name="Shapley")
+
+# https://en.wikipedia.org/wiki/Norma_Cluster
+# https://arxiv.org/abs/0706.2227
+NORMA = ObservedCluster(
+    RA=hms_to_degrees(16, 15, 32.8), dec=dms_to_degrees(-60, 53, 30),
+    dist=67.8 * 0.705, mass=1e15 * 0.7, name="Norma")
+
+# Wikipedia seems to give the wrong distance.
+# https://en.wikipedia.org/wiki/Leo_Cluster
+# https://arxiv.org/abs/astro-ph/0406367
+LEO = ObservedCluster(
+    RA=hms_to_degrees(11, 44, 36.5), dec=dms_to_degrees(19, 43, 32),
+    dist=91.3 * 0.705, mass=7e14 * 0.7, name="Leo")
+
+# https://en.wikipedia.org/wiki/Hydra_Cluster
+HYDRA = ObservedCluster(
+    RA=hms_to_degrees(9, 18), dec=dms_to_degrees(-12, 5),
+    dist=58.3 * 0.705, mass=4e15 * 0.7, name="Hydra")
+
+# I think this is Pisces? Not very sure about its mass.
+# https://en.wikipedia.org/wiki/Abell_262
+# https://arxiv.org/abs/0911.1774
+PISCES = ObservedCluster(
+    RA=hms_to_degrees(1, 52, 50.4), dec=dms_to_degrees(36, 8, 46),
+    dist=68.8 * 0.705, mass=2e14 * 0.7, name="Pisces")
+
+# This one is in the ZOA
+# https://en.wikipedia.org/wiki/Ophiuchus_Supercluster
+# https://arxiv.org/abs/1509.00986
+OPICHIUS = ObservedCluster(
+    RA=hms_to_degrees(17, 10, 0), dec=dms_to_degrees(-22),
+    dist=83.4, mass=1e15 * 0.7, name="Ophiuchus")
+
+
+clusters = {"Virgo": VIRGO,
+            "Fornax": FORNAX,
+            "Coma": COMA,
+            "Perseus": PERSEUS,
+            "Centaurus": CENTAURUS,
+            "Shapley": SHAPLEY,
+            "Norma": NORMA,
+            "Leo": LEO,
+            "Hydra": HYDRA,
+            "Pisces": PISCES,
+            "Opichius": OPICHIUS,
+            }

csiborgtools/flow/flow_model.py

@@ -212,7 +212,12 @@ class DataLoader:
             raise ValueError("Invalid simulation index.")
         nsim = nsims[ksim]
 
-        with File(paths.field_los(simname, catalogue), 'r') as f:
+        if "Pantheon+" in catalogue:
+            fpath = paths.field_los(simname, "Pantheon+")
+        else:
+            fpath = paths.field_los(simname, catalogue)
+
+        with File(fpath, 'r') as f:
             has_smoothed = True if f[f"density_{nsim}"].ndim > 2 else False
             if has_smoothed and (ksmooth is None or not isinstance(ksmooth, int)):  # noqa
                 raise ValueError("The output contains smoothed field but "
@@ -236,8 +241,13 @@ class DataLoader:
                 for key in f.keys():
                     arr[key] = f[key][:]
         elif catalogue in ["LOSS", "Foundation", "SFI_gals", "2MTF",
-                           "Pantheon+", "SFI_gals_masked", "SFI_groups"]:
+                           "Pantheon+", "SFI_gals_masked", "SFI_groups",
+                           "Pantheon+_groups", "Pantheon+_groups_zSN",
+                           "Pantheon+_zSN"]:
             with File(catalogue_fpath, 'r') as f:
+                if "Pantheon+" in catalogue:
+                    grp = f["Pantheon+"]
+                else:
                     grp = f[catalogue]
 
                 dtype = [(key, np.float32) for key in grp.keys()]
@@ -1381,15 +1391,26 @@ def get_model(loader, zcmb_max=None, verbose=True):
             los_overdensity[mask], los_velocity[mask], RA[mask], dec[mask],
             zCMB[mask], mB[mask], x1[mask], c[mask], e_mB[mask], e_x1[mask],
             e_c[mask], loader.rdist, loader._Omega_m)
-    elif kind == "Pantheon+":
+    elif "Pantheon+" in kind:
         keys = ["RA", "DEC", "zCMB", "mB", "x1", "c", "biasCor_m_b", "mBERR",
-                "x1ERR", "cERR", "biasCorErr_m_b"]
+                "x1ERR", "cERR", "biasCorErr_m_b", "zCMB_SN", "zCMB_Group"]
 
-        RA, dec, zCMB, mB, x1, c, bias_corr_mB, e_mB, e_x1, e_c, e_bias_corr_mB = (loader.cat[k] for k in keys)  # noqa
+        RA, dec, zCMB, mB, x1, c, bias_corr_mB, e_mB, e_x1, e_c, e_bias_corr_mB, zCMB_SN, zCMB_Group = (loader.cat[k] for k in keys)  # noqa
         mB -= bias_corr_mB
         e_mB = np.sqrt(e_mB**2 + e_bias_corr_mB**2)
 
         mask = (zCMB < zcmb_max)
+
+        if kind == "Pantheon+_groups":
+            mask &= np.isfinite(zCMB_Group)
+
+        if kind == "Pantheon+_groups_zSN":
+            mask &= np.isfinite(zCMB_Group)
+            zCMB = zCMB_SN
+
+        if kind == "Pantheon+_zSN":
+            zCMB = zCMB_SN
+
         model = SN_PV_validation_model(
             los_overdensity[mask], los_velocity[mask], RA[mask], dec[mask],
             zCMB[mask], mB[mask], x1[mask], c[mask], e_mB[mask], e_x1[mask],

csiborgtools/match/__init__.py

@@ -12,5 +12,6 @@
 # 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 .match import (ParticleOverlap, RealisationsMatcher, calculate_overlap,    # noqa
+from .overlap import (ParticleOverlap, RealisationsMatcher, calculate_overlap,  # noqa
                       find_neighbour, matching_max)                             # noqa
+from .obs_to_box import (MatchingProbability, MatchCatalogues)                  # noqa

csiborgtools/match/obs_to_box.py

@@ -0,0 +1,396 @@
+# Copyright (C) 2024 Richard Stiskalek
+# This program is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the
+# Free Software Foundation; either version 3 of the License, or (at your
+# option) any later version.
+#
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
+# Public License for more details.
+#
+# 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.
+"""
+Code to match observations to a constrained simulation.
+"""
+from abc import ABC
+
+import numpy as np
+from colossus.cosmology import cosmology
+from colossus.lss import mass_function
+from scipy.interpolate import interp1d
+from sklearn.neighbors import NearestNeighbors
+from tqdm import trange
+
+###############################################################################
+#                     Matching probability class                              #
+###############################################################################
+
+
+class BaseMatchingProbability(ABC):
+    """Base class for `MatchingProbability`."""
+
+    @property
+    def halo_pos(self):
+        """
+        Halo positions in the constrained simulation.
+
+        Returns
+        -------
+        2-dimensional array of shape `(n, 3)`
+        """
+        return self._halo_pos
+
+    @halo_pos.setter
+    def halo_pos(self, x):
+        if not isinstance(x, np.ndarray) and x.ndim == 2 and x.shape[1] == 3:
+            raise ValueError("Invalid halo positions.")
+        self._halo_pos = x
+
+    @property
+    def halo_log_mass(self):
+        """
+        Halo log mass in the constrained simulation.
+
+        Returns
+        -------
+        1-dimensional array of shape `(n,)`
+        """
+        return self._halo_log_mass
+
+    @halo_log_mass.setter
+    def halo_log_mass(self, x):
+        if not isinstance(x, np.ndarray) and x.ndim == 1 and len(x) != len(self.halo_pos):  # noqa
+            raise ValueError("Invalid halo log mass.")
+        self._halo_log_mass = x
+
+    @property
+    def nhalo(self):
+        """"
+        Number of haloes in the constrained simulation that are used for
+        matching.
+
+        Returns
+        -------
+        int
+        """
+        return self.halo_log_mass.size
+
+    def HMF(self, log_mass):
+        """
+        Evaluate the halo mass function at a given mass.
+
+        Parameters
+        ----------
+        log_mass : float
+            Logarithmic mass of the halo in `Msun / h`.
+
+        Returns
+        -------
+        HMF : float
+            The HMF in `h^3 Mpc^-3 dex^-1`.
+        """
+        return self._hmf(log_mass)
+
+
+class MatchingProbability(BaseMatchingProbability):
+    """"
+    Matching probability by calculating the CDF of finding a halo of a certain
+    mass at a given distance from a reference point. Calibrated against a HMF,
+    by assuming that the haloes are uniformly distributed. This is only
+    approximate treatment, as the haloes are not uniformly distributed, however
+    it is sufficient for the present purposes.
+
+    NOTE: The method currently does not account for uncertainty in distance.
+
+    Parameters
+    ----------
+    halo_pos : 2-dimensional array of shape `(n, 3)`
+        Halo positions in the constrained simulation in `Mpc / h`.
+    halo_log_mass : 1-dimensional array of shape `(n,)`
+        Halo log mass in the constrained simulation in `Msun / h`.
+    mdef : str, optional
+        Definition of the halo mass. Default is 'fof'.
+    cosmo_params : dict, optional
+        Cosmological parameters of the constrained simulation.
+    """
+    def __init__(self, halo_pos, halo_log_mass, mdef="fof",
+                 cosmo_params={'flat': True, 'H0': 67.66, 'Om0': 0.3111,
+                               'Ob0': 0.0489, 'sigma8': 0.8101, 'ns': 0.9665}):
+        self.halo_pos = halo_pos
+        self.halo_log_mass = halo_log_mass
+
+        # Define the kNN object and fit it to the halo positions, so that we
+        # can quickly query distances to an arbitrary point.
+        self._knn = NearestNeighbors()
+        self._knn.fit(halo_pos)
+
+        # Next, get the HMF from colossus and create its interpolant.
+        cosmology.addCosmology("myCosmo", **cosmo_params)
+        cosmology.setCosmology("myCosmo")
+
+        x = np.logspace(10, 16, 10000)
+        y = mass_function.massFunction(
+            x, 0.0, mdef=mdef, model="angulo12", q_out="dndlnM") * np.log(10)
+        self._hmf = interp1d(np.log10(x), y, kind="cubic")
+
+    def pdf(self, r, log_mass):
+        """
+        Calculate the PDF of finding a halo of a given mass at a given distance
+        from a random point.
+
+        Parameters
+        ----------
+        r : float
+            Distance from the random point in `Mpc / h`.
+        log_mass : float
+            Logarithmic mass of the halo in `Msun / h`.
+
+        Returns
+        -------
+        float
+        """
+        nd = self.HMF(log_mass)
+        return 4 * np.pi * r**2 * nd * np.exp(-4 / 3 * np.pi * r**3 * nd)
+
+    def cdf(self, r, log_mass):
+        """
+        Calculate the CDF of finding a halo of a given mass at a given distance
+        from a random point.
+
+        Parameters
+        ----------
+        r : float
+            Distance from the random point in `Mpc / h`.
+        log_mass : float
+            Logarithmic mass of the halo in `Msun / h`.
+
+        Returns
+        -------
+        float
+        """
+        nd = self.HMF(log_mass)
+        return 1 - np.exp(-4 / 3 * np.pi * r**3 * nd)
+
+    def inverse_cdf(self, cdf, log_mass):
+        """
+        Calculate the inverse CDF of finding a halo of a given mass at a given
+        distance from a random point.
+
+        Parameters
+        ----------
+        cdf : float
+            CDF of finding a halo of a given mass at a given distance.
+        log_mass : float
+            Logarithmic mass of the halo in `Msun / h`.
+
+        Returns
+        -------
+        float
+        """
+        nd = self.HMF(log_mass)
+        return (np.log(1 - cdf) / (-4 / 3 * np.pi * nd))**(1 / 3)
+
+    def cdf_per_halo(self, refpos, ref_log_mass=None, rmax=50,
+                     return_full=True):
+        """
+        Calculate the CDF per each halo in the constrained simulation.
+
+        Parameters
+        ----------
+        refpos : 1-dimensional array of shape `(3,)`
+            Reference position in `Mpc / h`.
+        ref_log_mass : float, optional
+            Reference log mass, used to calculate the difference in log mass
+            between the reference and each halo.
+        rmax : float, optional
+            Maximum distance from the reference point to consider. Below this,
+            the CDF is simply set to 1.
+        return_full : bool, optional
+            If `True`, return the CDF, dlogmass and indxs for all haloes,
+            otherwise return only the haloes within `rmax`.
+
+        Returns
+        -------
+        cdf : 1-dimensional array of shape `(nhalo,)`
+            CDF per halo.
+        dlogmass : 1-dimensional array of shape `(nhalo,)`
+            Difference in log mass between the reference and each halo.
+        indxs : 1-dimensional array of shape `(nhalo,)`
+            Indices of the haloes.
+        """
+        if not (isinstance(refpos, np.ndarray) and refpos.ndim == 1):
+            raise ValueError("Invalid reference position.")
+        if ref_log_mass is not None and not isinstance(ref_log_mass, (float, int, np.float32, np.float64)):  # noqa
+            raise ValueError("Invalid reference log mass.")
+
+        # Use the kNN  to pick out the haloes within `rmax` of the reference
+        # point.
+        dist, indxs = self._knn.radius_neighbors(
+            refpos.reshape(-1, 3), rmax, return_distance=True)
+        dist, indxs = dist[0], indxs[0]
+
+        cdf_ = self.cdf(dist, self.halo_log_mass[indxs])
+        if ref_log_mass is not None:
+            dlogmass_ = self.halo_log_mass[indxs] - ref_log_mass
+        else:
+            dlogmass_ = None
+
+        if return_full:
+            cdf = np.ones(self.nhalo)
+            cdf[indxs] = cdf_
+            if ref_log_mass is not None:
+                dlogmass = np.full(self.nhalo, np.infty)
+                dlogmass[indxs] = dlogmass_
+            else:
+                dlogmass = dlogmass_
+
+            indxs = np.arange(self.nhalo)
+        else:
+            cdf, dlogmass = cdf_, dlogmass_
+
+        return cdf, dlogmass, indxs
+
+    def match_halo(self, refpos, ref_log_mass, pvalue_threshold=0.005,
+                   max_absdlogmass=1., rmax=50, verbose=True,
+                   catalogue_index=0):
+        """
+        Match a halo in the constrained simulation to a reference halo.
+        Considers match the highest significance halo within `rmax` and
+        within `max_absdlogmass` of the reference halo mass. In case of no
+        match, returns `None`.
+
+        Parameters
+        ----------
+        refpos : 1-dimensional array of shape `(3,)`
+            Reference position.
+        ref_log_mass : float
+            Reference log mass.
+        pvalue_threshold : float, optional
+            Threshold for the CDF to be considered a match.
+        max_absdlogmass : float, optional
+            Maximum difference in log mass between the reference and the
+            matched halo.
+        rmax : float, optional
+            Maximum distance from the reference point to consider.
+        verbose : bool, optional
+            If `True`, print information about the match.
+        catalogue_index : int, optional
+            Optional catalogue index for more informative printing.
+
+        Returns
+        -------
+        cdf : float, or None
+            CDF of the matched halo (significance), if any.
+        index : int, or None
+            Index of the matched halo, if any.
+        """
+        cdf, dlogmass, indxs = self.cdf_per_halo(
+            refpos, ref_log_mass, rmax, return_full=False)
+
+        dlogmass = np.abs(dlogmass)
+        ks = np.argsort(cdf)
+        cdf, dlogmass, indxs = cdf[ks], dlogmass[ks], indxs[ks]
+
+        matches = np.where(
+            (cdf < pvalue_threshold) & (dlogmass < max_absdlogmass))[0]
+
+        if len(matches) == 0:
+            return None, None
+
+        if verbose and len(matches) > 1:
+            print(f"Found {len(matches)} plausible matches in catalogue {catalogue_index}.")  # noqa
+            for i, k in enumerate(matches):
+                j = indxs[k]
+                logM = self.halo_log_mass[j]
+                dx = np.linalg.norm(self.halo_pos[j] - refpos)
+                print(f"    {i + 1}: CDF = {cdf[k]:.3e}, index = {j}, logM = {logM:.3f} Msun / h, dx = {dx:.3f} Mpc / h.")  # noqa
+
+            print(flush=True)
+
+        k = matches[0]
+        return cdf[k], indxs[k]
+
+
+class MatchCatalogues:
+    """
+    A wrapper for `MatchingProbability` that allows to match observed clusters
+    to haloes in multiple catalogues.
+
+    Parameters
+    ----------
+    catalogues : list
+        List of halo catalogues of constrained simulations.
+    cosmo_params : dict, optional
+        Cosmological parameters of the constrained simulation to calculate
+        the corresponding FOF mass function.
+    """
+    def __init__(self, catalogues,
+                 cosmo_params={'flat': True, 'H0': 67.66, 'Om0': 0.3111,
+                               'Ob0': 0.0489, 'sigma8': 0.8101, 'ns': 0.9665}):
+        mdef = "fof"
+        self._catalogues = catalogues
+        self._prob_models = [None] * len(catalogues)
+
+        for i in trange(len(catalogues)):
+            pos = catalogues[i]["cartesian_pos"]
+            log_mass = np.log10(catalogues[i]["totmass"])
+
+            self._prob_models[i] = MatchingProbability(
+                pos, log_mass, mdef, cosmo_params)
+
+    def __getitem__(self, index):
+        return self._prob_models[index]
+
+    def __len__(self):
+        return len(self._catalogues)
+
+    def __call__(self, refpos, ref_log_mass, pvalue_threshold=0.05,
+                 max_absdlogmass=1., rmax=50, verbose=True):
+        """
+        Calculate the CDFs of finding a halo of a certain mass at a given
+        distance from a reference point for all catalogues.
+
+        Parameters
+        ----------
+        refpos : 1-dimensional array of shape `(3,)`
+            Reference position.
+        ref_log_mass : float
+            Reference log mass.
+        pvalue_threshold : float, optional
+            Threshold for the CDF to be considered a match.
+        max_absdlogmass : float, optional
+            Maximum difference in log mass between the reference and the
+            matched halo.
+        rmax : float, optional
+            Maximum distance from the reference point to consider.
+        verbose : bool, optional
+            Verbosity flag.
+
+        Returns
+        -------
+        cdfs : dict
+            Dictionary of CDFs per halo, with keys being the simulation
+            indices.
+        indxs : dict
+            Dictionary of indices of the matched haloes, with keys being the
+            simulation indices.
+        """
+        cdfs, indxs = {}, {}
+        for i in trange(len(self), desc="Matching catalogues",
+                        disable=not verbose):
+            cdf, indx = self._prob_models[i].match_halo(
+                refpos, ref_log_mass, pvalue_threshold, max_absdlogmass, rmax,
+                verbose, i)
+
+            if cdf is not None:
+                cdfs[i] = cdf
+                indxs[i] = indx
+
+        n = len(self) - len(cdfs)
+        if n > 0 and verbose:
+            print(f"Failed to assign haloes in {n} catalogues.")
+
+        return cdfs, indxs

csiborgtools/read/obs.py

@@ -26,6 +26,8 @@ from astropy.io import fits
 from astropy.cosmology import FlatLambdaCDM
 from scipy import constants
 
+from ..utils import radec_to_cartesian
+
 
 ###############################################################################
 #                           Text survey base class                            #
@@ -756,6 +758,42 @@ class BaseSingleObservation(ABC):
 
         self._spherical_pos = pos
 
+    @property
+    def mass(self):
+        """
+        Total mass estimate in Msun / h.
+
+        Returns
+        -------
+        float
+        """
+        if self._mass is None:
+            raise ValueError("`mass` is not set!")
+        return self._mass
+
+    @mass.setter
+    def mass(self, mass):
+        if not isinstance(mass, (int, float)):
+            raise ValueError("`mass` must be a float.")
+        self._mass = mass
+
+    def cartesian_pos(self, boxsize):
+        """
+        Cartesian position of the observation in Mpc / h, assuming the observer
+        is in the centre of the box.
+
+        Parameters
+        ----------
+        boxsize : float
+            Box size in Mpc / h.
+
+        Returns
+        -------
+        1-dimensional array of shape (3,)
+        """
+        return radec_to_cartesian(
+            self.spherical_pos.reshape(1, 3)).reshape(-1,) + boxsize / 2
+
     @property
     def name(self):
         """
@@ -788,13 +826,16 @@ class ObservedCluster(BaseSingleObservation):
         Declination in degrees.
     dist : float
         Distance in Mpc / h.
+    mass : float
+        Total mass estimate in Msun / h.
     name : str
         Cluster name.
     """
-    def __init__(self, RA, dec, dist, name):
+    def __init__(self, RA, dec, dist, mass, name):
         super().__init__()
         self.name = name
         self.spherical_pos = [dist, RA, dec]
+        self.mass = mass
 
 
 ###############################################################################

notebooks/flow/flow_calibration.py

@@ -38,7 +38,7 @@ def read_samples(catalogue, simname, ksmooth, include_calibration=False,
     Vx, Vy, Vz, beta, sigma_v, alpha = [], [], [], [], [], []
     BIC, AIC, logZ, chi2 = [], [], [], []
 
-    if catalogue in ["LOSS", "Foundation", "Pantheon+"]:
+    if catalogue in ["LOSS", "Foundation"] or "Pantheon+" in catalogue:
         alpha_cal, beta_cal, mag_cal, e_mu_intrinsic = [], [], [], []
     elif catalogue in ["2MTF", "SFI_gals", "SFI_gals_masked"]:
         a, b, e_mu_intrinsic = [], [], []
@@ -86,7 +86,7 @@ def read_samples(catalogue, simname, ksmooth, include_calibration=False,
             except KeyError:
                 chi2.append([0.])
 
-            if catalogue in ["LOSS", "Foundation", "Pantheon+"]:
+            if catalogue in ["LOSS", "Foundation"]  or "Pantheon+" in catalogue:  # noqa
                 alpha_cal.append(f[f"sim_{nsim}/alpha_cal"][:])
                 beta_cal.append(f[f"sim_{nsim}/beta_cal"][:])
                 mag_cal.append(f[f"sim_{nsim}/mag_cal"][:])
@@ -106,7 +106,7 @@ def read_samples(catalogue, simname, ksmooth, include_calibration=False,
 
     gof = np.hstack(BIC), np.hstack(AIC), np.hstack(logZ), np.hstack(chi2)
 
-    if catalogue in ["LOSS", "Foundation", "Pantheon+"]:
+    if catalogue in ["LOSS", "Foundation"] or "Pantheon+" in catalogue:
         alpha_cal, beta_cal, mag_cal, e_mu_intrinsic = np.hstack(alpha_cal), np.hstack(beta_cal), np.hstack(mag_cal), np.hstack(e_mu_intrinsic)  # noqa
     elif catalogue in ["2MTF", "SFI_gals", "SFI_gals_masked"]:
         a, b, e_mu_intrinsic = np.hstack(a), np.hstack(b), np.hstack(e_mu_intrinsic)  # noqa
@@ -118,6 +118,7 @@ def read_samples(catalogue, simname, ksmooth, include_calibration=False,
         raise ValueError(f"Catalogue {catalogue} not recognized.")
 
     # Calculate magnitude of V_ext
+
     Vmag = np.sqrt(Vx**2 + Vy**2 + Vz**2)
     # Calculate direction in galactic coordinates of V_ext
     V = np.vstack([Vx, Vy, Vz]).T
@@ -128,7 +129,7 @@ def read_samples(catalogue, simname, ksmooth, include_calibration=False,
     names = ["alpha", "beta", "Vmag", "l", "b", "sigma_v"]
 
     if include_calibration:
-        if catalogue in ["LOSS", "Foundation", "Pantheon+"]:
+        if catalogue in ["LOSS", "Foundation"] or "Pantheon+" in catalogue:
             data += [alpha_cal, beta_cal, mag_cal, e_mu_intrinsic]
             names += ["alpha_cal", "beta_cal", "mag_cal", "e_mu_intrinsic"]
         elif catalogue in ["2MTF", "SFI_gals", "SFI_gals_masked"]:

notebooks/match_observation/match_clusters.py

@@ -0,0 +1,25 @@
+import csiborgtools
+from tqdm import tqdm
+
+
+def open_cat(nsim, simname, bounds):
+    if "csiborg1" in simname:
+        cat = csiborgtools.read.CSiBORG1Catalogue(nsim, bounds=bounds)
+    elif "csiborg2" in simname:
+        cat = csiborgtools.read.CSiBORG2Catalogue(
+            nsim, 99, simname.split("_")[-1], bounds=bounds)
+    else:
+        raise ValueError(f"Unknown simulation name: {simname}.")
+
+    return cat
+
+
+def open_cats(simname, bounds):
+    paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
+    nsims = paths.get_ics(simname)
+    catalogues = [None] * len(nsims)
+
+    for i, nsim in enumerate(tqdm(nsims, desc="Opening catalogues")):
+        catalogues[i] = open_cat(nsim, simname, bounds)
+
+    return catalogues

scripts/field_los.py

@@ -61,7 +61,7 @@ def get_los(catalogue_name, simname, comm):
         if catalogue_name in ["LOSS", "Foundation", "SFI_gals",
                               "SFI_gals_masked", "SFI_groups", "2MTF",
                               "Pantheon+"]:
-            fpath = join(folder, "PV_compilation_Supranta2019.hdf5")
+            fpath = join(folder, "PV_compilation.hdf5")
             with File(fpath, 'r') as f:
                 grp = f[catalogue_name]
                 RA = grp["RA"][:]
@@ -286,7 +286,7 @@ if __name__ == "__main__":
 
     rmax = 200
     dr = 0.5
-    smooth_scales = [0, 2]
+    smooth_scales = [0, 2, 4]
 
     comm = MPI.COMM_WORLD
     paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)

scripts/flow_validation.py

@@ -52,8 +52,10 @@ def get_model(args, nsim_iterator, get_model_kwargs):
     if args.catalogue == "A2":
         fpath = join(folder, "A2.h5")
     elif args.catalogue in ["LOSS", "Foundation", "Pantheon+", "SFI_gals",
-                            "2MTF", "SFI_groups", "SFI_gals_masked"]:
+                            "2MTF", "SFI_groups", "SFI_gals_masked",
-        fpath = join(folder, "PV_compilation_Supranta2019.hdf5")
+                            "Pantheon+_groups", "Pantheon+_groups_zSN",
+                            "Pantheon+_zSN"]:
+        fpath = join(folder, "PV_compilation.hdf5")
     elif "CB2_" in args.catalogue:
         kind = args.catalogue.split("_")[-1]
         fpath = join(folder, f"PV_mock_CB2_17417_{kind}.hdf5")

scripts/flow_validation.sh

@@ -7,8 +7,8 @@ queue="berg"
 env="/mnt/users/rstiskalek/csiborgtools/venv_csiborg/bin/python"
 file="flow_validation.py"
 
-catalogue="CB2_small"
+catalogue="Pantheon+_groups"
-simname="csiborg2_main"
+simname="csiborg2_varysmall"
 
 
 pythoncm="$env $file --catalogue $catalogue --simname $simname --ksmooth $ksmooth"