Compute stellar asymmetry

README.md

@@ -18,11 +18,20 @@ The code is written in python and uses mpi4py for parallelism. Whichever
 install method you use, you will need MPI available so that mpi4py can be
 installed.
 
+We recommend cloning the repository, rather than installing SOAP directly from
+GitHub with pip. The example parameter files, test scripts, and the scripts used
+to generate the documentation are only available in the repository, and the
+commands in [Running SOAP](#running-soap) are run from the root of the
+repository. If you only want to import SOAP as a library, then it can be
+installed with `pip install git+https://github.com/SWIFTSIM/SOAP.git`.
+
 ### Quick install (serial HDF5)
 
-SOAP and its dependencies can be installed directly using the command
+SOAP and its dependencies can be installed using the commands
 ```
-pip install git+https://github.com/SWIFTSIM/SOAP.git
+git clone https://github.com/SWIFTSIM/SOAP.git
+cd SOAP
+pip install .
 ```
 This will usually install a serial version of h5py. SOAP will run correctly,
 but for large simulations the I/O will be slower.
@@ -36,11 +45,22 @@ source against that library before installing SOAP:
 ```
 pip install mpi4py
 export HDF5_MPI="ON"; export CC=mpicc; pip install --no-binary=h5py h5py
-pip install git+https://github.com/SWIFTSIM/SOAP.git
+git clone https://github.com/SWIFTSIM/SOAP.git
+cd SOAP
+pip install .
 ```
 If SOAP (and therefore serial h5py) is already installed, then add the flags
 `--no-cache-dir` and `--force-reinstall` when reinstalling h5py.
 
+### Optional dependencies
+
+Some properties require additional python packages. These properties are
+not computed unless they are explicitly enabled in the parameter file. To
+install the packages needed for all of these properties, run
+```
+pip install ".[extra_properties]"
+```
+
 ### Installation on COSMA
 
 If you are using the [COSMA system](https://cosma.readthedocs.io/en/latest/),
@@ -49,6 +69,9 @@ you can install an SOAP virtual environment by running
 
 ## Running SOAP
 
+The commands in this section should be run from the root of the SOAP
+repository.
+
 The command `./tests/run_small_volume.sh` will download a small example
 simulation, run the group membership and halo properties scripts on it.
 This uses the parameter file at `./tests/small_volume.yml`, and the
@@ -238,7 +261,8 @@ If your property is expensive to compute, or requires additional dependencies,
 set `opt_in_reason` in its `SOAP/property_table.py` entry to a short (20
 characters or less) reason. It will then only be calculated if it is explicitly
 enabled in the parameter file. Any additional dependencies must be imported
-within the `@lazy_property`, not at the top of the file.
+within the `@lazy_property`, not at the top of the file, and added to the
+`extra_properties` section of `pyproject.toml`.
 
 If SOAP does crash while evaluating your new property it will try to
 output the ID of the halo it was processing when it crashed. Then you

SOAP/particle_selection/aperture_properties.py

@@ -3649,6 +3649,166 @@ def StellarInertiaTensorReducedNoniterativeLuminosityWeighted(
             reduced=True, max_iterations=1
         )
 
+    @lazy_property
+    def shrinking_sphere_centre(
+        self, min_particles=200, shrink_factor=0.83, max_iter=256
+    ) -> unyt.unyt_array:
+        """
+        Estimate the galaxy center (center of mass) with iterative shrinking aperture.
+
+        Procedure:
+        1) Initialize center as the COM of all input particles.
+        2) Initialize radius as max(aperture_radius, farthest particle distance).
+        3) Recompute COM inside current aperture, recenter, then shrink radius by shrink_factor.
+        4) Stop when enclosed particle count is < min_particles, and return the last valid center.
+
+        Parameters
+        ----------
+        min_particles : int, default=200
+            Stop iteration when enclosed particle count is below this threshold, or 10% of particle count.
+        shrink_factor : float, default=0.83
+            Radius scaling factor applied at each iteration (0 < shrink_factor < 1).
+        max_iter : int, default=256
+            Maximum number of iterations
+        """
+
+        min_particles = min(min_particles, 0.1 * self.Nstar)
+
+        if self.Mstar == 0:
+            return None
+
+        centre = (self.star_mass_fraction[:, None] * self.pos_star).sum(axis=0)
+        dist = np.linalg.norm(self.pos_star - centre, axis=1)
+        radius = np.max(dist) * 1.01
+
+        for i in range(max_iter):
+            mask = dist <= radius
+            n_in = np.sum(mask)
+            if n_in < min_particles:
+                break
+
+            centre = (self.mass_star[mask, None] * self.pos_star[mask]).sum(axis=0)
+            centre /= np.sum(self.mass_star[mask])
+            dist = np.linalg.norm(self.pos_star - centre, axis=1)
+            radius *= shrink_factor
+
+        return centre
+
+    @lazy_property
+    def ShrinkingSphereCentre(self) -> unyt.unyt_array:
+        """
+        Centre computed by applying shrinking spheres method to stars
+        """
+        if self.Mstar == 0:
+            return None
+
+        return (self.shrinking_sphere_centre + self.centre) % self.boxsize
+
+    @lazy_property
+    def StellarAsymmetry(self):
+        return self.stellar_asymmetry()
+
+    @lazy_property
+    def StellarAsymmetryShrink(self):
+        return self.stellar_asymmetry(shrink_centre=True)
+
+    @lazy_property
+    def StellarAsymmetrySubsample(self):
+        return self.stellar_asymmetry(N=8)
+
+    @lazy_property
+    def StellarAsymmetry48(self):
+        return self.stellar_asymmetry(npix=48)
+
+    @lazy_property
+    def StellarAsymmetry192(self):
+        return self.stellar_asymmetry(npix=192)
+
+    def stellar_asymmetry(self, npix=12, N=None, shrink_centre=False):
+        """
+        Compute stellar asymmetry following https://arxiv.org/abs/1805.03210
+
+        Parameters
+        ----------
+        npix : int, default=12
+            Number of equal-area angular regions (HEALPix pixels) for directional
+            partitioning. Must satisfy npix = 12 * nside^2 where nside is a
+            power of 2 (e.g. 12, 48, 192, 768, ...).
+        N : int, optional
+            Randomly select 1/N of the original stellar particle set before
+            computing the asymmetry. If the total number of particles is smaller
+            than or equal to N, all particles are used.
+
+        """
+
+        import healpy as hp
+
+        if self.Mstar == 0:
+            return None
+
+        # Check we are using a valid value for npix
+        nside = int(round(np.sqrt(npix // 12)))
+        assert npix == 12 * nside**2
+        assert nside.bit_count() == 1
+
+        if N is None:
+            indices = slice(None)
+        else:
+            N = int(N)
+            if N <= 0:
+                raise ValueError("N must be a positive integer")
+            if self.Nstar <= N:
+                indices = slice(None)
+            else:
+                # Seed with the halo index so results are reproducible
+                rng = np.random.default_rng(int(self.index))
+                indices = rng.choice(self.Nstar, size=self.Nstar // N, replace=False)
+
+        mass_star = self.mass_star[indices]
+
+        # Centre using the shrinking sphere
+        if shrink_centre:
+            pos = self.pos_star[indices] - self.shrinking_sphere_centre
+        else:
+            pos = self.pos_star[indices]
+        r = np.linalg.norm(pos, axis=1)
+
+        # Remove particles close to the centre, they are symmetric
+        mask = r.to_value("kpc") < 0.1
+        if np.sum(mask):
+            pos = pos[np.logical_not(mask)]
+            r = r[np.logical_not(mask)]
+            mass_star = mass_star[np.logical_not(mask)]
+            if r.shape[0] == 0:
+                return np.float32(0)
+
+        if mass_star.shape[0] < 3:
+            return None
+
+        # Compute the mass in each pixel
+        vecs = (pos / r[:, None]).value
+        idx = hp.vec2pix(nside, vecs[:, 0], vecs[:, 1], vecs[:, 2])
+        # np.bincount will not return a unyt array
+        region_mass_msun = np.bincount(
+            idx,
+            weights=mass_star.to_value("Msun"),
+            minlength=npix,
+        )
+
+        # Create antipodal mapping
+        # Find the center vector of every pixel
+        vecs = hp.pix2vec(nside, np.arange(npix))
+        # Negate vectors to find antipodal points
+        anti_vecs = -np.array(vecs)
+        # Map those points back to pixel IDs
+        anti_indices = hp.vec2pix(nside, anti_vecs[0], anti_vecs[1], anti_vecs[2])
+
+        # Calculate asymmetry
+        mass_diff = np.abs(region_mass_msun - region_mass_msun[anti_indices])
+        asymmetry = np.sum(mass_diff) / (2.0 * mass_star.sum().to_value("Msun"))
+
+        return asymmetry
+
 
 class ApertureProperties(HaloProperty):
     """
@@ -3821,6 +3981,12 @@ class ApertureProperties(HaloProperty):
         "StellarInertiaTensorReducedLuminosityWeighted": True,
         "StellarInertiaTensorNoniterativeLuminosityWeighted": False,
         "StellarInertiaTensorReducedNoniterativeLuminosityWeighted": False,
+        "ShrinkingSphereCentre": False,
+        "StellarAsymmetry": False,
+        "StellarAsymmetryShrink": False,
+        "StellarAsymmetrySubsample": False,
+        "StellarAsymmetry48": False,
+        "StellarAsymmetry192": False,
     }
 
     property_list = {

SOAP/property_table.py

@@ -222,6 +222,13 @@ class PropertyTable:
         ],
         "footnote_compY.tex": ["compY", "compY_no_agn"],
         "footnote_dopplerB.tex": ["DopplerB"],
+        "footnote_asymmetry.tex": [
+            "StellarAsymmetry",
+            "StellarAsymmetryShrink",
+            "StellarAsymmetrySubsample",
+            "StellarAsymmetry48",
+            "StellarAsymmetry192",
+        ],
         "footnote_coreexcision.tex": [
             "Tgas_cy_weighted_core_excision",
             "Tgas_cy_weighted_core_excision_no_agn",
@@ -3754,6 +3761,84 @@ class PropertyTable:
             output_physical=False,
             a_scale_exponent=1,
         ),
+        "ShrinkingSphereCentre": Property(
+            name="ShrinkingSphereCentre",
+            shape=3,
+            dtype=np.float64,
+            unit="snap_length",
+            description="Shrinking sphere centre computed using stars.",
+            lossy_compression_filter="DScale6",
+            dmo_property=False,
+            particle_properties=["PartType4/Coordinates", "PartType4/Masses"],
+            output_physical=False,
+            a_scale_exponent=1,
+            opt_in_reason="Expensive to compute",
+        ),
+        "StellarAsymmetry": Property(
+            name="StellarAsymmetry",
+            shape=1,
+            dtype=np.float32,
+            unit="dimensionless",
+            description="Asymmetry of the stellar mass distribution around the halo centre, computed using 12 HEALPix pixels.",
+            lossy_compression_filter="FMantissa9",
+            dmo_property=False,
+            particle_properties=["PartType4/Coordinates", "PartType4/Masses"],
+            output_physical=True,
+            a_scale_exponent=0,
+            opt_in_reason="Requires healpy",
+        ),
+        "StellarAsymmetryShrink": Property(
+            name="StellarAsymmetryShrink",
+            shape=1,
+            dtype=np.float32,
+            unit="dimensionless",
+            description="As StellarAsymmetry, but centred on ShrinkingSphereCentre rather than the halo centre.",
+            lossy_compression_filter="FMantissa9",
+            dmo_property=False,
+            particle_properties=["PartType4/Coordinates", "PartType4/Masses"],
+            output_physical=True,
+            a_scale_exponent=0,
+            opt_in_reason="Requires healpy",
+        ),
+        "StellarAsymmetrySubsample": Property(
+            name="StellarAsymmetrySubsample",
+            shape=1,
+            dtype=np.float32,
+            unit="dimensionless",
+            description="As StellarAsymmetry, but computed using a random subsample of 1/8 of the star particles. All star particles are used if there are 8 or fewer.",
+            lossy_compression_filter="FMantissa9",
+            dmo_property=False,
+            particle_properties=["PartType4/Coordinates", "PartType4/Masses"],
+            output_physical=True,
+            a_scale_exponent=0,
+            opt_in_reason="Requires healpy",
+        ),
+        "StellarAsymmetry48": Property(
+            name="StellarAsymmetry48",
+            shape=1,
+            dtype=np.float32,
+            unit="dimensionless",
+            description="As StellarAsymmetry, but using 48 HEALPix pixels.",
+            lossy_compression_filter="FMantissa9",
+            dmo_property=False,
+            particle_properties=["PartType4/Coordinates", "PartType4/Masses"],
+            output_physical=True,
+            a_scale_exponent=0,
+            opt_in_reason="Requires healpy",
+        ),
+        "StellarAsymmetry192": Property(
+            name="StellarAsymmetry192",
+            shape=1,
+            dtype=np.float32,
+            unit="dimensionless",
+            description="As StellarAsymmetry, but using 192 HEALPix pixels.",
+            lossy_compression_filter="FMantissa9",
+            dmo_property=False,
+            particle_properties=["PartType4/Coordinates", "PartType4/Masses"],
+            output_physical=True,
+            a_scale_exponent=0,
+            opt_in_reason="Requires healpy",
+        ),
         "compY": Property(
             name="ComptonY",
             shape=1,

documentation/footnote_asymmetry.tex

@@ -0,0 +1,18 @@
+\paragraph{$^{$FOOTNOTE_NUMBER$}$The stellar asymmetry}\label{footnote:$FOOTNOTE_NUMBER$} is computed following
+\href{https://arxiv.org/abs/1805.03210}{arXiv:1805.03210}. The directions around the centre are divided into
+$N_{\rm{}pix}$ equal-area HEALPix pixels, and the asymmetry is
+
+\begin{equation}
+    A = \frac{1}{2 M_*} \sum_{i=1}^{N_{\rm{}pix}} \left| M_i - M_{\bar{i}} \right|,
+\end{equation}
+
+where $M_i$ is the stellar mass in pixel $i$, $M_{\bar{i}}$ is the stellar mass in the pixel in the opposite
+direction, and $M_*$ is the total stellar mass. Each pair of opposite pixels appears twice in the sum, so $A$
+ranges from 0 for a perfectly symmetric distribution to 1 if all the stellar mass is on one side of the centre.
+Star particles within 0.1 kpc of the centre are excluded, since their direction from the centre is poorly
+defined. The asymmetry is not calculated if there are fewer than 3 star particles.
+
+\verb+StellarAsymmetry+ uses $N_{\rm{}pix}=12$ and is centred on the halo centre. The other variants differ
+as follows: \verb+StellarAsymmetryShrink+ is centred on \verb+ShrinkingSphereCentre+,
+\verb+StellarAsymmetrySubsample+ uses a random subsample of 1/8 of the star particles, and
+\verb+StellarAsymmetry48+ and \verb+StellarAsymmetry192+ use $N_{\rm{}pix}=48$ and $N_{\rm{}pix}=192$.

pyproject.toml

@@ -32,9 +32,14 @@ dependencies = [
 ]
 
 [project.optional-dependencies]
+# Dependencies which are only required by some opt-in properties
+extra_properties = [
+    "healpy",
+]
 test = [
     "pytest-mpi",
     "numba",
+    "SOAP[extra_properties]",
 ]
 
 [project.urls]