diff --git a/README.md b/README.md index dbc66807..494f37c7 100644 --- a/README.md +++ b/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 diff --git a/SOAP/particle_selection/aperture_properties.py b/SOAP/particle_selection/aperture_properties.py index 625579c4..54029599 100644 --- a/SOAP/particle_selection/aperture_properties.py +++ b/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 = { diff --git a/SOAP/property_table.py b/SOAP/property_table.py index 5c1ce372..ad9aad7b 100644 --- a/SOAP/property_table.py +++ b/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, diff --git a/documentation/footnote_asymmetry.tex b/documentation/footnote_asymmetry.tex new file mode 100644 index 00000000..bec8af86 --- /dev/null +++ b/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$. diff --git a/pyproject.toml b/pyproject.toml index 44f13618..3847a191 100644 --- a/pyproject.toml +++ b/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]