From 163c12f4181bfff33e010c7312e837b3c8492107 Mon Sep 17 00:00:00 2001 From: Kyle Oman Date: Thu, 9 Oct 2025 10:02:10 +0100 Subject: [PATCH 01/19] Update README.md (#194) --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 048c839c..469fee13 100644 --- a/README.md +++ b/README.md @@ -19,7 +19,7 @@ steps are recommended for install ``` 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@soap_runtime +pip install git+https://github.com/SWIFTSIM/SOAP.git ``` ### Installation on COSMA From 6dd580f2b8d7dd7357c21b6079c740c66fd51648 Mon Sep 17 00:00:00 2001 From: robjmcgibbon <36136863+robjmcgibbon@users.noreply.github.com> Date: Fri, 10 Oct 2025 17:18:53 +0100 Subject: [PATCH 02/19] Update EAGLE sbatch script (#195) --- scripts/EAGLE.sh | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/scripts/EAGLE.sh b/scripts/EAGLE.sh index 9d5f5130..8f87ef72 100755 --- a/scripts/EAGLE.sh +++ b/scripts/EAGLE.sh @@ -9,6 +9,14 @@ #SBATCH --exclusive #SBATCH -t 02:00:00 # +# For L0025N0752 set ntasks=16 +# For L0100N1504 set ntasks=256 +# +# Install Hdecompose with: +# pip install git+ssh://git@github.com/kyleaoman/Hdecompose.git +# +# Download virtual snapshot script with: +# wget https://gitlab.cosma.dur.ac.uk/swift/swiftsim/-/raw/master/tools/create_virtual_snapshot.py sim_name='L0100N1504' snap_nr="028" @@ -17,7 +25,6 @@ z_suffix="z000p000" module purge module load python/3.12.4 gnu_comp/14.1.0 openmpi/5.0.3 parallel_hdf5/1.12.3 source openmpi-5.0.3-hdf5-1.12.3-env/bin/activate -pip install git+ssh://git@github.com/kyleaoman/Hdecompose.git ######## Link files to snap (to remove awful z suffix) @@ -65,12 +72,16 @@ mpirun -- python -u misc/hdecompose_hydrogen_fractions.py \ ######### Create virtual snapshot # Must be run from the snapshot directory itself or there will be issues with paths +soap_dir=$(pwd) cd "${output_dir}/swift_snapshots/swift_${snap_nr}" -wget https://gitlab.cosma.dur.ac.uk/swift/swiftsim/-/raw/master/tools/create_virtual_snapshot.py -python create_virtual_snapshot.py "snap_${snap_nr}.0.hdf5" -rm create_virtual_snapshot.py* +python "${soap_dir}/create_virtual_snapshot.py" "snap_${snap_nr}.0.hdf5" cd - +python compression/make_virtual_snapshot.py \ + "${output_dir}/swift_snapshots/swift_${snap_nr}/snap_${snap_nr}.hdf5" \ + "${output_dir}/SOAP_uncompressed/membership_${snap_nr}/membership_${snap_nr}.{file_nr}.hdf5" \ + "${output_dir}/SOAP_uncompressed/snap_${snap_nr}.hdf5" \ + ######### Run SOAP chunks=10 From 1a5c5f0de090d8045392d3624d2805ad7adb6dcb Mon Sep 17 00:00:00 2001 From: robjmcgibbon <36136863+robjmcgibbon@users.noreply.github.com> Date: Thu, 16 Oct 2025 17:47:20 +0100 Subject: [PATCH 03/19] Import os (#196) --- SOAP/group_membership.py | 1 + 1 file changed, 1 insertion(+) diff --git a/SOAP/group_membership.py b/SOAP/group_membership.py index 1b0db8ff..35cc7e0d 100644 --- a/SOAP/group_membership.py +++ b/SOAP/group_membership.py @@ -1,5 +1,6 @@ #!/bin/env python +import os import time import socket From 37bbadd9b8d70b521751340c3276ad1d2888f677 Mon Sep 17 00:00:00 2001 From: robjmcgibbon <36136863+robjmcgibbon@users.noreply.github.com> Date: Thu, 23 Oct 2025 15:50:09 +0100 Subject: [PATCH 04/19] Signed snapshot isolation (#197) --- SOAP/property_table.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/SOAP/property_table.py b/SOAP/property_table.py index 9666b4ac..0a8f541c 100644 --- a/SOAP/property_table.py +++ b/SOAP/property_table.py @@ -4635,7 +4635,7 @@ class PropertyTable: "HBTplus/SnapshotOfLastMaxMass": Property( name="HBTplus/SnapshotOfLastMaxMass", shape=1, - dtype=np.uint64, + dtype=np.int64, unit="dimensionless", description="Latest snapshot when this subhalo had its maximum mass.", lossy_compression_filter="None", @@ -4659,7 +4659,7 @@ class PropertyTable: "HBTplus/SnapshotOfLastMaxVmax": Property( name="HBTplus/SnapshotOfLastMaxVmax", shape=1, - dtype=np.uint64, + dtype=np.int64, unit="dimensionless", description="Latest snapshot when this subhalo had its largest maximum circular velocity.", lossy_compression_filter="None", @@ -4671,7 +4671,7 @@ class PropertyTable: "HBTplus/SnapshotOfLastIsolation": Property( name="HBTplus/SnapshotOfLastIsolation", shape=1, - dtype=np.uint64, + dtype=np.int64, unit="dimensionless", description="Latest snapshot when this subhalo was a central. -1 if the subhalo has always been a central.", lossy_compression_filter="None", From 10b59c3c8f11f616325d38fcc36599adf6649300 Mon Sep 17 00:00:00 2001 From: "Victor J. Forouhar Moreno" <60343271+VictorForouhar@users.noreply.github.com> Date: Mon, 27 Oct 2025 15:49:25 +0100 Subject: [PATCH 05/19] Rotation velocity and velocity dispersion additions (#192) * Lower required precision between cumsum and sum quantities Running on L200m6 crashed for half light radii (difference was 0.13% instead of less than 0.1%) * Implement helper function for mass-weighted rotational velocity * Define azimuthal velocity for subhalo apertures. * Add helper function for cylindrical velocity dispersion vector * Define cylindrical dispersion quantities for aperture properties * Add missing property declaration * Define luminosity-weighted kinematic properties Still need to think of how to handle different angular momentum vectors. * Change calculate_cylindrical_velocities to take reference pos and vel * Define cylindrical velocity calculation on a per-luminosity-band basis * Implementation of luminosity-weighted rotational velocity Testing pending. * Implement luminosity-weighted cylindrical dispersion routines. Tests to be done * Define luminosity weighted kinematics for aperture properties * Add new properties to parameter file * Add new properties to property table * Incorrect function call * Declare properties for bound subhalo * Fix indexing of arrays * Do not recentre star coordinates for cylindrical velocities We do not do so when computing the angular momentum vector, so this change makes it consistent with our choice to not re-centre. * Use stellar centre of mass as reference velocity Including the case for when we do luminosity-weighting. This reflects the same choice we made when computing luminosity-weighted angular momenta. * Use STELLAR CoM velocity as reference, not ALL CoM velocity. * Fix bug: accidental change of velocities within function * Run formatter * Update property table. --- .../particle_selection/aperture_properties.py | 171 ++++++++++++++++++ SOAP/particle_selection/subhalo_properties.py | 150 +++++++++++---- .../cylindrical_coordinates.py | 22 ++- SOAP/property_calculation/half_mass_radius.py | 2 +- .../kinematic_properties.py | 166 +++++++++++++++++ SOAP/property_table.py | 72 +++++++- parameter_files/COLIBRE_HYBRID.yml | 12 ++ parameter_files/COLIBRE_THERMAL.yml | 12 ++ 8 files changed, 567 insertions(+), 40 deletions(-) diff --git a/SOAP/particle_selection/aperture_properties.py b/SOAP/particle_selection/aperture_properties.py index 6a07f6cb..5d0cbe4b 100644 --- a/SOAP/particle_selection/aperture_properties.py +++ b/SOAP/particle_selection/aperture_properties.py @@ -147,6 +147,10 @@ def MetalFracStar(self): ) from SOAP.property_calculation.kinematic_properties import ( get_velocity_dispersion_matrix, + get_cylindrical_velocity_dispersion_vector_mass_weighted, + get_cylindrical_velocity_dispersion_vector_luminosity_weighted, + get_rotation_velocity_mass_weighted, + get_rotation_velocity_luminosity_weighted, get_angular_momentum, get_angular_momentum_and_kappa_corot_mass_weighted, get_angular_momentum_and_kappa_corot_luminosity_weighted, @@ -156,6 +160,9 @@ def MetalFracStar(self): get_inertia_tensor_mass_weighted, get_inertia_tensor_luminosity_weighted, ) +from SOAP.property_calculation.cylindrical_coordinates import ( + calculate_cylindrical_velocities, +) from SOAP.core.swift_cells import SWIFTCellGrid from SOAP.property_calculation.stellar_age_calculator import StellarAgeCalculator from SOAP.particle_filter.cold_dense_gas_filter import ColdDenseGasFilter @@ -1448,6 +1455,162 @@ def veldisp_matrix_star(self) -> unyt.unyt_array: self.star_mass_fraction, self.vel_star, self.vcom_star ) + @lazy_property + def star_cylindrical_velocities(self) -> unyt.unyt_array: + """ + Calculate the velocities of the star particles in cyclindrical + coordinates, where the axes are centred on the stellar CoM, + and the z axis is aligned with the stellar angular momentum. + """ + + # We need at least 2 particles to have an angular momentum vector + if self.Nstar < 2: + return None + + # This can happen if we have particles on top of each other + # or with the same velocity + if np.sum(self.Lstar) == 0: + return None + + # Get velocities in cylindrical coordinates + return calculate_cylindrical_velocities( + self.pos_star, + self.vel_star, + self.Lstar, + reference_velocity=self.vcom_star, + ) + + @lazy_property + def StellarRotationalVelocity(self) -> unyt.unyt_array: + """ + Mass-weighted average azimuthal velocity of stars. + """ + if (self.Nstar < 2) or (np.sum(self.Lstar) == 0): + return None + return get_rotation_velocity_mass_weighted( + self.mass_star, self.star_cylindrical_velocities[:, 1] + ) + + @lazy_property + def StellarCylindricalVelocityDispersionVector(self) -> unyt.unyt_array: + if (self.Nstar < 2) or (np.sum(self.Lstar) == 0): + return None + return get_cylindrical_velocity_dispersion_vector_mass_weighted( + self.mass_star, self.star_cylindrical_velocities + ) + + @lazy_property + def StellarCylindricalVelocityDispersion(self) -> unyt.unyt_array: + if self.StellarCylindricalVelocityDispersionVector is None: + return None + return np.sqrt((self.StellarCylindricalVelocityDispersionVector**2).sum() / 3) + + @lazy_property + def StellarCylindricalVelocityDispersionVertical(self) -> unyt.unyt_array: + if self.StellarCylindricalVelocityDispersionVector is None: + return None + return self.StellarCylindricalVelocityDispersionVector[2] + + @lazy_property + def StellarCylindricalVelocityDispersionDiscPlane(self) -> unyt.unyt_array: + if self.StellarCylindricalVelocityDispersionVector is None: + return None + return np.sqrt((self.StellarCylindricalVelocityDispersionVector[:2] ** 2).sum()) + + @lazy_property + def star_cylindrical_velocities_luminosity_weighted(self) -> unyt.unyt_array: + """ + Calculate the velocities of the star particles in cylindrical + coordinates, where the origin and reference frame are centred on the + stellar centre of light in each band. The z axis is aligned with the + stellar angular momentum obtained from each band. + """ + + # We need at least 2 particles to have an angular momentum vector + if self.Nstar < 2: + return None + + # This can happen if we have particles on top of each other + # or with the same velocity + if np.sum(self.Lstar_luminosity_weighted) == 0: + return None + + # We iterate over bands to use their own reference vector and luminosity- + # weighted centre of mass phase space coordinates. + cylindrical_velocities = ( + np.zeros( + ( + self.stellar_luminosities.shape[1], + self.stellar_luminosities.shape[0], + 3, + ) + ) + * self.vel_star.units + ) + for i_band, particle_luminosities_i_band in enumerate( + self.stellar_luminosities.T + ): + cylindrical_velocities[i_band] = calculate_cylindrical_velocities( + self.pos_star, + self.vel_star, + self.Lstar_luminosity_weighted[i_band * 3 : (1 + i_band) * 3], + reference_velocity=self.vcom_star, + ) + + return cylindrical_velocities + + @lazy_property + def StellarRotationalVelocityLuminosityWeighted(self) -> unyt.unyt_array: + if (self.Nstar < 2) or (np.sum(self.Lstar) == 0): + return None + return get_rotation_velocity_luminosity_weighted( + self.stellar_luminosities, + self.star_cylindrical_velocities_luminosity_weighted[:, :, 1], + ) + + @lazy_property + def StellarCylindricalVelocityDispersionVectorLuminosityWeighted( + self, + ) -> unyt.unyt_array: + if (self.Nstar < 2) or (np.sum(self.Lstar) == 0): + return None + return get_cylindrical_velocity_dispersion_vector_luminosity_weighted( + self.stellar_luminosities, + self.star_cylindrical_velocities_luminosity_weighted, + ) + + @lazy_property + def StellarCylindricalVelocityDispersionLuminosityWeighted(self) -> unyt.unyt_array: + if self.StellarCylindricalVelocityDispersionVectorLuminosityWeighted is None: + return None + return np.sqrt( + ( + self.StellarCylindricalVelocityDispersionVectorLuminosityWeighted**2 + ).sum(axis=1) + / 3 + ) + + @lazy_property + def StellarCylindricalVelocityDispersionVerticalLuminosityWeighted( + self, + ) -> unyt.unyt_array: + if self.StellarCylindricalVelocityDispersionVectorLuminosityWeighted is None: + return None + return self.StellarCylindricalVelocityDispersionVectorLuminosityWeighted[:, 2] + + @lazy_property + def StellarCylindricalVelocityDispersionDiscPlaneLuminosityWeighted( + self, + ) -> unyt.unyt_array: + if self.StellarCylindricalVelocityDispersionVectorLuminosityWeighted is None: + return None + return np.sqrt( + ( + self.StellarCylindricalVelocityDispersionVectorLuminosityWeighted[:, :2] + ** 2 + ).sum(axis=1) + ) + @lazy_property def KineticEnergyStars(self) -> unyt.unyt_quantity: """ @@ -3583,6 +3746,14 @@ class ApertureProperties(HaloProperty): "veldisp_matrix_gas": False, "veldisp_matrix_dm": False, "veldisp_matrix_star": False, + "StellarRotationalVelocity": False, + "StellarCylindricalVelocityDispersion": False, + "StellarCylindricalVelocityDispersionVertical": False, + "StellarCylindricalVelocityDispersionDiscPlane": False, + "StellarRotationalVelocityLuminosityWeighted": False, + "StellarCylindricalVelocityDispersionLuminosityWeighted": False, + "StellarCylindricalVelocityDispersionVerticalLuminosityWeighted": False, + "StellarCylindricalVelocityDispersionDiscPlaneLuminosityWeighted": False, "KineticEnergyGas": False, "KineticEnergyStars": False, "Mgas_SF": False, diff --git a/SOAP/particle_selection/subhalo_properties.py b/SOAP/particle_selection/subhalo_properties.py index a11231f5..bd2c94cd 100644 --- a/SOAP/particle_selection/subhalo_properties.py +++ b/SOAP/particle_selection/subhalo_properties.py @@ -33,6 +33,10 @@ get_angular_momentum_and_kappa_corot_luminosity_weighted, get_vmax, get_velocity_dispersion_matrix, + get_cylindrical_velocity_dispersion_vector_mass_weighted, + get_cylindrical_velocity_dispersion_vector_luminosity_weighted, + get_rotation_velocity_mass_weighted, + get_rotation_velocity_luminosity_weighted, ) from SOAP.property_calculation.inertia_tensors import ( get_inertia_tensor_mass_weighted, @@ -1389,7 +1393,7 @@ def star_mass_fraction(self) -> unyt.unyt_array: @lazy_property def star_cylindrical_velocities(self) -> unyt.unyt_array: """ - Calculate the velocities of the star particles in cyclindrical + Calculate the velocities of the star particles in cylindrical coordinates, where the axes are centred on the stellar CoM, and the z axis is aligned with the stellar angular momentum. """ @@ -1403,63 +1407,141 @@ def star_cylindrical_velocities(self) -> unyt.unyt_array: if np.sum(self.Lstar) == 0: return None - # Calculate the position of the stars relative to their CoM - pos = self.pos_star - (self.star_mass_fraction[:, None] * self.pos_star).sum( - axis=0 - ) - - # Calculate relative velocity of stars - vrel = self.vel_star - self.vcom[None, :] - # Get velocities in cylindrical coordinates return calculate_cylindrical_velocities( - pos, - vrel, + self.pos_star, + self.vel_star, self.Lstar, + reference_velocity=self.vcom_star, ) @lazy_property def StellarRotationalVelocity(self) -> unyt.unyt_array: if (self.Nstar < 2) or (np.sum(self.Lstar) == 0): return None - v_cylindrical = self.star_cylindrical_velocities - v_phi = v_cylindrical[:, 1] - return (self.star_mass_fraction * v_phi).sum() + return get_rotation_velocity_mass_weighted( + self.mass_star, self.star_cylindrical_velocities[:, 1] + ) @lazy_property - def stellar_cylindrical_squared_velocity_dispersion_vector(self) -> unyt.unyt_array: + def StellarCylindricalVelocityDispersionVector(self) -> unyt.unyt_array: if (self.Nstar < 2) or (np.sum(self.Lstar) == 0): return None - v_cylindrical = self.star_cylindrical_velocities + return get_cylindrical_velocity_dispersion_vector_mass_weighted( + self.mass_star, self.star_cylindrical_velocities + ) - # This implementation of standard deviation is more numerically stable than using - ^2 - mean_velocity = (self.star_mass_fraction[:, None] * v_cylindrical).sum(axis=0) - squared_velocity_dispersion = ( - self.star_mass_fraction[:, None] * (v_cylindrical - mean_velocity) ** 2 - ).sum(axis=0) + @lazy_property + def StellarCylindricalVelocityDispersion(self) -> unyt.unyt_array: + if self.StellarCylindricalVelocityDispersionVector is None: + return None + return np.sqrt((self.StellarCylindricalVelocityDispersionVector**2).sum() / 3) - return squared_velocity_dispersion + @lazy_property + def StellarCylindricalVelocityDispersionVertical(self) -> unyt.unyt_array: + if self.StellarCylindricalVelocityDispersionVector is None: + return None + return self.StellarCylindricalVelocityDispersionVector[2] @lazy_property - def StellarCylindricalVelocityDispersion(self) -> unyt.unyt_array: - if self.stellar_cylindrical_squared_velocity_dispersion_vector is None: + def StellarCylindricalVelocityDispersionDiscPlane(self) -> unyt.unyt_array: + if self.StellarCylindricalVelocityDispersionVector is None: + return None + return np.sqrt((self.StellarCylindricalVelocityDispersionVector[:2] ** 2).sum()) + + @lazy_property + def star_cylindrical_velocities_luminosity_weighted(self) -> unyt.unyt_array: + """ + Calculate the velocities of the star particles in cylindrical + coordinates, where the origin and reference frame are centred on the + stellar centre of light in each band. The z axis is aligned with the + stellar angular momentum obtained from each band. + """ + + # We need at least 2 particles to have an angular momentum vector + if self.Nstar < 2: + return None + + # This can happen if we have particles on top of each other + # or with the same velocity + if np.sum(self.Lstar_luminosity_weighted) == 0: + return None + + # We iterate over bands to use their own reference vector and luminosity- + # weighted centre of mass phase space coordinates. + cylindrical_velocities = ( + np.zeros( + ( + self.stellar_luminosities.shape[1], + self.stellar_luminosities.shape[0], + 3, + ) + ) + * self.vel_star.units + ) + for i_band, particle_luminosities_i_band in enumerate( + self.stellar_luminosities.T + ): + cylindrical_velocities[i_band] = calculate_cylindrical_velocities( + self.pos_star, + self.vel_star, + self.Lstar_luminosity_weighted[i_band * 3 : (1 + i_band) * 3], + reference_velocity=self.vcom_star, + ) + + return cylindrical_velocities + + @lazy_property + def StellarRotationalVelocityLuminosityWeighted(self) -> unyt.unyt_array: + if (self.Nstar < 2) or (np.sum(self.Lstar) == 0): + return None + + return get_rotation_velocity_luminosity_weighted( + self.stellar_luminosities, + self.star_cylindrical_velocities_luminosity_weighted[:, :, 1], + ) + + @lazy_property + def StellarCylindricalVelocityDispersionVectorLuminosityWeighted( + self, + ) -> unyt.unyt_array: + if (self.Nstar < 2) or (np.sum(self.Lstar) == 0): + return None + return get_cylindrical_velocity_dispersion_vector_luminosity_weighted( + self.stellar_luminosities, + self.star_cylindrical_velocities_luminosity_weighted, + ) + + @lazy_property + def StellarCylindricalVelocityDispersionLuminosityWeighted(self) -> unyt.unyt_array: + if self.StellarCylindricalVelocityDispersionVectorLuminosityWeighted is None: return None return np.sqrt( - self.stellar_cylindrical_squared_velocity_dispersion_vector.sum() / 3 + ( + self.StellarCylindricalVelocityDispersionVectorLuminosityWeighted**2 + ).sum(axis=1) + / 3 ) @lazy_property - def StellarCylindricalVelocityDispersionVertical(self) -> unyt.unyt_array: - if self.stellar_cylindrical_squared_velocity_dispersion_vector is None: + def StellarCylindricalVelocityDispersionVerticalLuminosityWeighted( + self, + ) -> unyt.unyt_array: + if self.StellarCylindricalVelocityDispersionVectorLuminosityWeighted is None: return None - return np.sqrt(self.stellar_cylindrical_squared_velocity_dispersion_vector[2]) + return self.StellarCylindricalVelocityDispersionVectorLuminosityWeighted[:, 2] @lazy_property - def StellarCylindricalVelocityDispersionDiscPlane(self) -> unyt.unyt_array: - if self.stellar_cylindrical_squared_velocity_dispersion_vector is None: + def StellarCylindricalVelocityDispersionDiscPlaneLuminosityWeighted( + self, + ) -> unyt.unyt_array: + if self.StellarCylindricalVelocityDispersionVectorLuminosityWeighted is None: return None return np.sqrt( - self.stellar_cylindrical_squared_velocity_dispersion_vector[:2].sum() + ( + self.StellarCylindricalVelocityDispersionVectorLuminosityWeighted[:, :2] + ** 2 + ).sum(axis=1) ) @lazy_property @@ -2307,10 +2389,14 @@ class SubhaloProperties(HaloProperty): "Lgas", "Ldm", "Lstar", + "StellarRotationalVelocity", "StellarCylindricalVelocityDispersion", "StellarCylindricalVelocityDispersionVertical", "StellarCylindricalVelocityDispersionDiscPlane", - "StellarRotationalVelocity", + "StellarRotationalVelocityLuminosityWeighted", + "StellarCylindricalVelocityDispersionLuminosityWeighted", + "StellarCylindricalVelocityDispersionVerticalLuminosityWeighted", + "StellarCylindricalVelocityDispersionDiscPlaneLuminosityWeighted", "kappa_corot_gas", "kappa_corot_star", "Lbaryons", diff --git a/SOAP/property_calculation/cylindrical_coordinates.py b/SOAP/property_calculation/cylindrical_coordinates.py index 223ac266..410b01c4 100644 --- a/SOAP/property_calculation/cylindrical_coordinates.py +++ b/SOAP/property_calculation/cylindrical_coordinates.py @@ -42,7 +42,9 @@ def build_rotation_matrix(z_target): return R -def calculate_cylindrical_velocities(positions, velocities, z_target): +def calculate_cylindrical_velocities( + positions, velocities, z_target, reference_position=None, reference_velocity=None +): """ Convert 3D Cartesian velocities to cylindrical coordinates (v_r, v_phi, v_z), after rotating the system such that the z-axis aligns with `z_target`. @@ -51,16 +53,30 @@ def calculate_cylindrical_velocities(positions, velocities, z_target): positions: (N, 3) array of particle positions in the original Cartesian frame. velocities: (N, 3) array of particle velocities in the original Cartesian frame. z_target: A 3-element vector indicating the new z-axis direction. + reference_position: (3,) array with a reference position on which to centre the Cartesian coordinate system. + reference_velocity: (3,) array with a reference velocity on which to centre the Cartesian reference frame. Returns: cyl_velocities: (N, 3) array of velocities in cylindrical coordinates: [v_r, v_phi, v_z] for each particle. """ + + # We need to declare a relative pos/vel array to not overwrite the original + # values. + if reference_position is None: + prel = positions + else: + prel = positions - reference_position + if reference_velocity is None: + vrel = velocities + else: + vrel = velocities - reference_velocity + R = build_rotation_matrix(z_target) # Rotate positions and velocities into new frame - positions_rot = positions @ R.T - velocities_rot = velocities @ R.T + positions_rot = prel @ R.T + velocities_rot = vrel @ R.T x = positions_rot[:, 0] y = positions_rot[:, 1] diff --git a/SOAP/property_calculation/half_mass_radius.py b/SOAP/property_calculation/half_mass_radius.py index 1daecf1d..13896b16 100644 --- a/SOAP/property_calculation/half_mass_radius.py +++ b/SOAP/property_calculation/half_mass_radius.py @@ -52,7 +52,7 @@ def get_half_weight_radius( # Consistency check. # np.sum() and np.cumsum() use different orders, so we have to allow for # some small difference. - if cumulative_weights[-1] < 0.999 * total_weight: + if cumulative_weights[-1] < 0.998 * total_weight: raise RuntimeError( "Weights sum up to less than the given total weight value:" f" cumulative_weights[-1] = {cumulative_weights[-1]}," diff --git a/SOAP/property_calculation/kinematic_properties.py b/SOAP/property_calculation/kinematic_properties.py index 14aeab68..1e5c900f 100644 --- a/SOAP/property_calculation/kinematic_properties.py +++ b/SOAP/property_calculation/kinematic_properties.py @@ -14,6 +14,80 @@ import unyt +def get_weighted_rotation_velocity( + particle_weights: unyt.unyt_array, particle_azimuthal_velocities: unyt.unyt_array +) -> unyt.unyt_quantity: + """ + Get the weighted average azimuthal velocity of a particle distribution. + + Parameters: + - particle_weights: unyt.unyt_array + Weight assigned to each particle. + - particle_azimuthal_velocities: unyt.unyt_array + Azimuthal velocity of each particle. + + Returns: + - Weighted average of the azimuthal velocity of particles. + """ + return (particle_weights * particle_azimuthal_velocities).sum() + + +def get_rotation_velocity_mass_weighted( + particle_masses, particle_azimuthal_velocities +) -> unyt.unyt_quantity: + """ + Return the mass-weighted average azimuthal velocity of a particle distribution. + + Parameters: + - particle_masses: unyt.unyt_array + Mass of particle. + - particle_azimuthal_velocities: unyt.unyt_array + Azimuthal velocity of each particle. + + Returns: + - Mass-weighted average of the azimuthal velocity of particles. + """ + mass_weights = particle_masses / particle_masses.sum() + return get_weighted_rotation_velocity(mass_weights, particle_azimuthal_velocities) + + +def get_rotation_velocity_luminosity_weighted( + particle_luminosities: unyt.unyt_array, + particle_azimuthal_velocities: unyt.unyt_array, +) -> unyt.unyt_array: + """ + Return the luminosity-weighted average azimuthal velocity of a particle distribution, for each + provided luminosity band. + + Parameters: + - particle_luminosities: unyt.unyt_array + Luminosity of each particle in the provided bands. + - particle_azimuthal_velocities: unyt.unyt_array + Azimuthal velocity of each particle, pre-computed for each luminosity band. + + Returns: + - Luminosity-weighted average of the azimuthal velocity of particles, with a value for each band. + """ + + number_luminosity_bands = particle_luminosities.shape[1] + rotation_velocities = ( + np.zeros(number_luminosity_bands) * particle_azimuthal_velocities.units + ) + + for i_band, ( + particle_luminosities_i_band, + particle_azimuthal_velocities_i_band, + ) in enumerate(zip(particle_luminosities.T, particle_azimuthal_velocities)): + luminosity_weights = ( + particle_luminosities_i_band / particle_luminosities_i_band.sum() + ) + rotation_velocities[i_band] = get_weighted_rotation_velocity( + luminosity_weights, particle_azimuthal_velocities_i_band + ) + + return rotation_velocities + + def get_velocity_dispersion_matrix( mass_fraction: unyt.unyt_array, velocity: unyt.unyt_array, @@ -53,6 +127,98 @@ def get_velocity_dispersion_matrix( return result +def get_weighted_cylindrical_velocity_dispersion_vector( + particle_weights: unyt.unyt_array, + particle_cylindrical_velocities: unyt.unyt_array, +) -> unyt.unyt_array: + """ + Compute the velocity dispersion along the radial, azimuthal and vertical + directions for the input particles and their specified weights. + + Parameters: + - particle_weights: unyt.unyt_array + Weight assigned to each particle. + - particle_cylindrical_velocities: unyt.unyt_array + Velocity of the particles in a cylindrical coordinate system. + + Returns a 3 element vector containing [sigma_r, sigma_phi, sigma_z]. + """ + + # This implementation of standard deviation is more numerically stable than using - ^2 + mean_velocity = (particle_weights[:, None] * particle_cylindrical_velocities).sum( + axis=0 + ) + squared_velocity_dispersion = ( + particle_weights[:, None] + * (particle_cylindrical_velocities - mean_velocity) ** 2 + ).sum(axis=0) + + return np.sqrt(squared_velocity_dispersion) + + +def get_cylindrical_velocity_dispersion_vector_mass_weighted( + particle_masses: unyt.unyt_array, + particle_cylindrical_velocities: unyt.unyt_array, +) -> unyt.unyt_array: + """ + Compute the mass-weighted velocity dispersion along the radial, azimuthal and vertical + directions for the input particles. + + Parameters: + - particle_masses: unyt.unyt_array + Mass of each particle. + - particle_cylindrical_velocities: unyt.unyt_array + Velocity of each particle in a cylindrical coordinate system. + + Returns a 3 element vector containing [sigma_r, sigma_phi, sigma_z]. + """ + mass_weights = particle_masses / particle_masses.sum() + return get_weighted_cylindrical_velocity_dispersion_vector( + mass_weights, particle_cylindrical_velocities + ) + + +def get_cylindrical_velocity_dispersion_vector_luminosity_weighted( + particle_luminosities: unyt.unyt_array, + particle_cylindrical_velocities: unyt.unyt_array, +) -> unyt.unyt_array: + """ + Compute the luminosity-weighted velocity dispersion along the radial, azimuthal and vertical + directions for the input particles. We return a vector for each luminosity band. + + Parameters: + - particle_luminosities: unyt.unyt_array + Luminosity of each particle in different luminosity bands. + - particle_cylindrical_velocities: unyt.unyt_array + Velocity of each particle in a cylindrical coordinate system, which varies between different + luminosity bands. + + Returns a 3 element vector for each luminosity band, which contains [sigma_r, sigma_phi, sigma_z]. + The velocity dispersion vectors for each band are appended to the same vector, hence the shape is + (number_luminosity_bands, 3). + """ + + number_luminosity_bands = particle_luminosities.shape[1] + velocity_dispersion_vectors = ( + np.zeros((number_luminosity_bands, 3)) * particle_cylindrical_velocities.units + ) + + for i_band, ( + particle_luminosities_i_band, + particle_cylindrical_velocities_i_band, + ) in enumerate(zip(particle_luminosities.T, particle_cylindrical_velocities)): + luminosity_weights = ( + particle_luminosities_i_band / particle_luminosities_i_band.sum() + ) + velocity_dispersion_vectors[i_band] = ( + get_weighted_cylindrical_velocity_dispersion_vector( + luminosity_weights, particle_cylindrical_velocities_i_band + ) + ) + + return velocity_dispersion_vectors + + def get_angular_momentum( mass: unyt.unyt_array, position: unyt.unyt_array, diff --git a/SOAP/property_table.py b/SOAP/property_table.py index 0a8f541c..583f291f 100644 --- a/SOAP/property_table.py +++ b/SOAP/property_table.py @@ -1657,7 +1657,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="snap_length/snap_time", - description="Mass weighted mean rotational velocity of the stars, in a cylindrical coordinate system where the axes are centred on the stellar CoM, and the z axis is aligned with the stellar angular momentum.", + description="Mass-weighted mean rotational velocity of the stars, in a cylindrical coordinate system where the axes are centred on the halo centre, the z axis is aligned with the mass-weighted stellar angular momentum, and stellar velocities are relative to the stellar centre of mass velocity.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=[ @@ -1673,7 +1673,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="snap_length/snap_time", - description="One-dimensional velocity dispersion of the stars computed in a cylindrical coordinate system where the axes are centred on the stellar CoM, and the z axis is aligned with the stellar angular momentum.", + description="One-dimensional mass-weighted velocity dispersion of the star, in a cylindrical coordinate system where the axes are centred on the halo centre, the z axis is aligned with the mass-weighted stellar angular momentum, and stellar velocities are relative to the stellar centre of mass velocity.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=[ @@ -1689,7 +1689,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="snap_length/snap_time", - description="Velocity dispersion perpendicular to the orbital plane of the stars, computed in a cylindrical coordinate system where the axes are centred on the stellar CoM, and the z axis is aligned with the stellar angular momentum.", + description="Mass-weighted velocity dispersion perpendicular to the orbital plane of the stars, in a cylindrical coordinate system where the axes are centred on the halo centre, the z axis is aligned with the mass-weighted stellar angular momentum, and stellar velocities are relative to the stellar centre of mass velocity.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=[ @@ -1705,7 +1705,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="snap_length/snap_time", - description="Total velocity dispersion in the orbital plane of the stars, computed in a cylindrical coordinate system where the axes are centred on the stellar CoM, and the z axis is aligned with the stellar angular momentum.", + description="Mass-weighted total velocity dispersion in the orbital plane of the stars, in a cylindrical coordinate system where the axes are centred on the halo centre, the z axis is aligned with the mass-weighted stellar angular momentum, and stellar velocities are relative to the stellar centre of mass velocity.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=[ @@ -1716,6 +1716,70 @@ class PropertyTable: output_physical=True, a_scale_exponent=0, ), + "StellarRotationalVelocityLuminosityWeighted": Property( + name="StellarRotationalVelocityLuminosityWeighted", + shape=9, # GAMA bands + dtype=np.float32, + unit="snap_length/snap_time", + description="Luminosity-weighted mean rotational velocity of the stars, in a cylindrical coordinate system where the axes are centred on the halo centre, the z axis is aligned with the luminosity-weighted stellar angular momentum, and stellar velocities are relative to the stellar centre of mass velocity.", + lossy_compression_filter="FMantissa9", + dmo_property=False, + particle_properties=[ + "PartType4/Coordinates", + "PartType4/Luminosities", + "PartType4/Velocities", + ], + output_physical=True, + a_scale_exponent=0, + ), + "StellarCylindricalVelocityDispersionLuminosityWeighted": Property( + name="StellarCylindricalVelocityDispersionLuminosityWeighted", + shape=9, # GAMA bands + dtype=np.float32, + unit="snap_length/snap_time", + description="One-dimensional luminosity-weighted velocity dispersion of the stars, in a cylindrical coordinate system where the axes are centred on the halo centre, the z axis is aligned with the luminosity-weighted stellar angular momentum, and stellar velocities are relative to the stellar centre of mass velocity.", + lossy_compression_filter="FMantissa9", + dmo_property=False, + particle_properties=[ + "PartType4/Coordinates", + "PartType4/Luminosities", + "PartType4/Velocities", + ], + output_physical=True, + a_scale_exponent=0, + ), + "StellarCylindricalVelocityDispersionVerticalLuminosityWeighted": Property( + name="StellarCylindricalVelocityDispersionVerticalLuminosityWeighted", + shape=9, # GAMA bands + dtype=np.float32, + unit="snap_length/snap_time", + description="Luminosity-weighted velocity dispersion perpendicular to the orbital plane of the stars, in a cylindrical coordinate system where the axes are centred on the halo centre, the z axis is aligned with the luminosity-weighted stellar angular momentum, and stellar velocities are relative to the stellar centre of mass velocity.", + lossy_compression_filter="FMantissa9", + dmo_property=False, + particle_properties=[ + "PartType4/Coordinates", + "PartType4/Luminosities", + "PartType4/Velocities", + ], + output_physical=True, + a_scale_exponent=0, + ), + "StellarCylindricalVelocityDispersionDiscPlaneLuminosityWeighted": Property( + name="StellarCylindricalVelocityDispersionDiscPlaneLuminosityWeighted", + shape=9, # GAMA bands + dtype=np.float32, + unit="snap_length/snap_time", + description="Luminosity-weighted total velocity dispersion in the orbital plane of the stars, in a cylindrical coordinate system where the axes are centred on the halo centre, the z axis is aligned with the luminosity-weighted stellar angular momentum, and stellar velocities are relative to the stellar centre of mass velocity.", + lossy_compression_filter="FMantissa9", + dmo_property=False, + particle_properties=[ + "PartType4/Coordinates", + "PartType4/Luminosities", + "PartType4/Velocities", + ], + output_physical=True, + a_scale_exponent=0, + ), "MaximumStellarBirthDensity": Property( name="MaximumStellarBirthDensity", shape=1, diff --git a/parameter_files/COLIBRE_HYBRID.yml b/parameter_files/COLIBRE_HYBRID.yml index a8de3712..0f9b3f6e 100644 --- a/parameter_files/COLIBRE_HYBRID.yml +++ b/parameter_files/COLIBRE_HYBRID.yml @@ -236,6 +236,14 @@ ApertureProperties: StellarMassFractionInMetals: true StellarMassFractionInOxygen: true StellarVelocityDispersionMatrix: general + StellarRotationalVelocity: false + StellarCylindricalVelocityDispersion: false + StellarCylindricalVelocityDispersionVertical: false + StellarCylindricalVelocityDispersionDiscPlane: false + StellarRotationalVelocityLuminosityWeighted: false + StellarCylindricalVelocityDispersionLuminosityWeighted: false + StellarCylindricalVelocityDispersionVerticalLuminosityWeighted: false + StellarCylindricalVelocityDispersionDiscPlaneLuminosityWeighted: false TotalMass: true TotalSNIaRate: true GasMassInColdDenseDiffuseMetals: @@ -730,6 +738,10 @@ SubhaloProperties: StellarCylindricalVelocityDispersion: general StellarCylindricalVelocityDispersionVertical: false StellarCylindricalVelocityDispersionDiscPlane: false + StellarRotationalVelocityLuminosityWeighted: false + StellarCylindricalVelocityDispersionLuminosityWeighted: false + StellarCylindricalVelocityDispersionVerticalLuminosityWeighted: false + StellarCylindricalVelocityDispersionDiscPlaneLuminosityWeighted: false aliases: PartType0/LastSNIIKineticFeedbackDensities: PartType0/DensitiesAtLastSupernovaEvent PartType0/LastSNIIThermalFeedbackDensities: PartType0/DensitiesAtLastSupernovaEvent diff --git a/parameter_files/COLIBRE_THERMAL.yml b/parameter_files/COLIBRE_THERMAL.yml index 8b2ddb4e..403267a5 100644 --- a/parameter_files/COLIBRE_THERMAL.yml +++ b/parameter_files/COLIBRE_THERMAL.yml @@ -236,6 +236,14 @@ ApertureProperties: StellarMassFractionInMetals: true StellarMassFractionInOxygen: true StellarVelocityDispersionMatrix: general + StellarCylindricalVelocityDispersion: false + StellarCylindricalVelocityDispersionVertical: false + StellarCylindricalVelocityDispersionDiscPlane: false + StellarCylindricalVelocityDispersionLuminosityWeighted: false + StellarCylindricalVelocityDispersionVerticalLuminosityWeighted: false + StellarCylindricalVelocityDispersionDiscPlaneLuminosityWeighted: false + StellarRotationalVelocity: false + StellarRotationalVelocityLuminosityWeighted: false TotalMass: true TotalSNIaRate: true GasMassInColdDenseDiffuseMetals: @@ -730,6 +738,10 @@ SubhaloProperties: StellarCylindricalVelocityDispersion: general StellarCylindricalVelocityDispersionVertical: false StellarCylindricalVelocityDispersionDiscPlane: false + StellarRotationalVelocityLuminosityWeighted: false + StellarCylindricalVelocityDispersionLuminosityWeighted: false + StellarCylindricalVelocityDispersionVerticalLuminosityWeighted: false + StellarCylindricalVelocityDispersionDiscPlaneLuminosityWeighted: false aliases: PartType0/LastSNIIKineticFeedbackDensities: PartType0/DensitiesAtLastSupernovaEvent PartType0/LastSNIIThermalFeedbackDensities: PartType0/DensitiesAtLastSupernovaEvent From e103b84d02bbe68220f90ac3bc34fefcf719631e Mon Sep 17 00:00:00 2001 From: robjmcgibbon <36136863+robjmcgibbon@users.noreply.github.com> Date: Mon, 8 Dec 2025 12:08:00 +0000 Subject: [PATCH 06/19] Add timing to matching script (#202) --- misc/match_group_membership.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/misc/match_group_membership.py b/misc/match_group_membership.py index 7424a9ba..081078b1 100644 --- a/misc/match_group_membership.py +++ b/misc/match_group_membership.py @@ -23,6 +23,7 @@ """ import argparse +import datetime import os from mpi4py import MPI @@ -319,6 +320,8 @@ def mpi_print(string, comm_rank): if __name__ == "__main__": + start_time = datetime.datetime.now() + parser = argparse.ArgumentParser( description=( "Script to match halos across runs by comparing particles" @@ -503,4 +506,6 @@ def mpi_print(string, comm_rank): phdf5.collective_write(file, "MatchCount2to1", match_count_21, comm=comm) phdf5.collective_write(file, "Consistent2to1", consistent_21, comm=comm) + comm.barrier() + mpi_print(f"Runtime: {datetime.datetime.now() - start_time}", comm_rank) mpi_print("Done!", comm_rank) From a7b3c4576055bdb2dfec1abb15303c89b212ad6c Mon Sep 17 00:00:00 2001 From: "Victor J. Forouhar Moreno" <60343271+VictorForouhar@users.noreply.github.com> Date: Wed, 14 Jan 2026 12:18:15 +0100 Subject: [PATCH 07/19] Spin parameter dates back to 2001 (#203) --- SOAP/particle_selection/SO_properties.py | 2 +- SOAP/particle_selection/subhalo_properties.py | 2 +- documentation/footnote_spin.tex | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/SOAP/particle_selection/SO_properties.py b/SOAP/particle_selection/SO_properties.py index 7faaf700..55746315 100644 --- a/SOAP/particle_selection/SO_properties.py +++ b/SOAP/particle_selection/SO_properties.py @@ -604,7 +604,7 @@ def spin_parameter(self) -> unyt.unyt_quantity: """ Spin parameter of all particles in the spherical overdensity. - Computed as in Bullock et al. (2021): + Computed as in Bullock et al. (2001): lambda = |Ltot| / (sqrt(2) * M * v_max * R) """ if self.Mtotpart == 0: diff --git a/SOAP/particle_selection/subhalo_properties.py b/SOAP/particle_selection/subhalo_properties.py index bd2c94cd..cbb68649 100644 --- a/SOAP/particle_selection/subhalo_properties.py +++ b/SOAP/particle_selection/subhalo_properties.py @@ -1032,7 +1032,7 @@ def spin_parameter(self) -> unyt.unyt_quantity: """ Spin parameter of all particles in the subhalo. - Computed as in Bullock et al. (2021): + Computed as in Bullock et al. (2001): lambda = |Ltot| / (sqrt(2) * M * v_max * R) Since a subhalo does not have a characteristic radius, R, we instead use diff --git a/documentation/footnote_spin.tex b/documentation/footnote_spin.tex index cf2ddc3c..f089562d 100644 --- a/documentation/footnote_spin.tex +++ b/documentation/footnote_spin.tex @@ -1,4 +1,4 @@ -\paragraph{$^{$FOOTNOTE_NUMBER$}$The spin parameter}\label{footnote:$FOOTNOTE_NUMBER$} is computed following Bullock et al. (2021): +\paragraph{$^{$FOOTNOTE_NUMBER$}$The spin parameter}\label{footnote:$FOOTNOTE_NUMBER$} is computed following Bullock et al. (2001): \begin{equation} \lambda{} = \frac{|\vec{L}_{\rm{}tot}|}{\sqrt{2}M v_{\rm{}max} R}, From e5aa7fa3226b5d7174b35757626b71c6dcc72bfc Mon Sep 17 00:00:00 2001 From: robjmcgibbon <36136863+robjmcgibbon@users.noreply.github.com> Date: Thu, 26 Feb 2026 10:53:51 +0000 Subject: [PATCH 08/19] Set version in requirements (#206) --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index 461cebf0..6cf2915d 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,6 +1,6 @@ mpi4py h5py -numpy>=2 +numpy>=2,<2.4 unyt>=3 astropy>=6 scipy From 86c3bfdb2c62b11f0a8cade67e4853575fe89354 Mon Sep 17 00:00:00 2001 From: robjmcgibbon <36136863+robjmcgibbon@users.noreply.github.com> Date: Mon, 16 Mar 2026 16:21:36 +0000 Subject: [PATCH 09/19] Calculate birth halo catalogue index (#204) * Virtual v1 * Combine update_vds_path with make_virtual_snapshot * Virtual snapshot creation using multiple auxilary files * BirthHaloCatalogueIndex v1 * Identify gas progenitor * Add ExSitu fraction property to SOAP * Add sbatch script * Remove extra parameter file * Format * Learn to spell * Add to BoundSubhalo * Handle named columns * Format * Move compression dir * Format --- README.md | 2 +- SOAP/catalogue_readers/read_hbtplus.py | 2 +- {compression => SOAP/compression}/README.md | 0 .../compression}/compress_soap_catalogue.py | 0 .../create_empty_SOAP_catalogue.py | 0 .../compression}/extract_filters.py | 0 {compression => SOAP/compression}/filters.yml | 0 SOAP/compression/make_virtual_snapshot.py | 362 ++++++++++++++++++ .../compression}/wrong_compression.yml | 0 SOAP/compute_halo_properties.py | 1 - SOAP/core/combine_chunks.py | 2 +- SOAP/core/halo_tasks.py | 1 - .../recently_heated_gas_filter.py | 1 - .../particle_selection/aperture_properties.py | 19 + SOAP/particle_selection/subhalo_properties.py | 17 + SOAP/property_table.py | 18 +- compression/make_virtual_snapshot.py | 192 ---------- compression/update_vds_paths.py | 118 ------ format.sh | 2 +- misc/check_subhalo_ranking.py | 1 - misc/compute_BirthHaloCatalogueIndex.py | 285 ++++++++++++++ parameter_files/COLIBRE_HYBRID.yml | 2 + parameter_files/COLIBRE_THERMAL.yml | 2 + scripts/COLIBRE/compress_group_membership.sh | 5 +- scripts/COLIBRE/compress_halo_properties.sh | 4 +- scripts/COLIBRE/compute_birth_index.sh | 47 +++ scripts/EAGLE.sh | 8 +- .../compress_halo_properties_L1000N1800.sh | 2 +- tests/COLIBRE/run_L0025N0188_Thermal.sh | 2 +- tests/test_SO_properties.py | 2 +- 30 files changed, 768 insertions(+), 329 deletions(-) rename {compression => SOAP/compression}/README.md (100%) rename {compression => SOAP/compression}/compress_soap_catalogue.py (100%) rename {compression => SOAP/compression}/create_empty_SOAP_catalogue.py (100%) rename {compression => SOAP/compression}/extract_filters.py (100%) rename {compression => SOAP/compression}/filters.yml (100%) create mode 100644 SOAP/compression/make_virtual_snapshot.py rename {compression => SOAP/compression}/wrong_compression.yml (100%) delete mode 100644 compression/make_virtual_snapshot.py delete mode 100644 compression/update_vds_paths.py create mode 100644 misc/compute_BirthHaloCatalogueIndex.py create mode 100755 scripts/COLIBRE/compute_birth_index.sh diff --git a/README.md b/README.md index 469fee13..7b72e845 100644 --- a/README.md +++ b/README.md @@ -116,7 +116,7 @@ The first is lossless compression via GZIP, the second is lossy compression. For the group membership files we only apply lossless compression. However, each property in the final SOAP catalogue has a lossy compression filter associated with it, which are set in `SOAP/property_table.py`. The script -`compression/compress_soap_catalogue.py` will apply both lossy and +`SOAP/compression/compress_soap_catalogue.py` will apply both lossy and lossless compression to SOAP catalogues. ### Documentation diff --git a/SOAP/catalogue_readers/read_hbtplus.py b/SOAP/catalogue_readers/read_hbtplus.py index d91c5849..4dff8d93 100644 --- a/SOAP/catalogue_readers/read_hbtplus.py +++ b/SOAP/catalogue_readers/read_hbtplus.py @@ -264,7 +264,7 @@ def read_hbtplus_catalogue( MassInMsunh = None VelInKmS = None sorted_file = None - (LengthInMpch, MassInMsunh, VelInKmS) = comm.bcast( + LengthInMpch, MassInMsunh, VelInKmS = comm.bcast( (LengthInMpch, MassInMsunh, VelInKmS) ) sorted_file = comm.bcast(sorted_file) diff --git a/compression/README.md b/SOAP/compression/README.md similarity index 100% rename from compression/README.md rename to SOAP/compression/README.md diff --git a/compression/compress_soap_catalogue.py b/SOAP/compression/compress_soap_catalogue.py similarity index 100% rename from compression/compress_soap_catalogue.py rename to SOAP/compression/compress_soap_catalogue.py diff --git a/compression/create_empty_SOAP_catalogue.py b/SOAP/compression/create_empty_SOAP_catalogue.py similarity index 100% rename from compression/create_empty_SOAP_catalogue.py rename to SOAP/compression/create_empty_SOAP_catalogue.py diff --git a/compression/extract_filters.py b/SOAP/compression/extract_filters.py similarity index 100% rename from compression/extract_filters.py rename to SOAP/compression/extract_filters.py diff --git a/compression/filters.yml b/SOAP/compression/filters.yml similarity index 100% rename from compression/filters.yml rename to SOAP/compression/filters.yml diff --git a/SOAP/compression/make_virtual_snapshot.py b/SOAP/compression/make_virtual_snapshot.py new file mode 100644 index 00000000..2311840e --- /dev/null +++ b/SOAP/compression/make_virtual_snapshot.py @@ -0,0 +1,362 @@ +#!/bin/env python + +import os.path +import h5py +import shutil +import numpy as np + + +class SafeDict(dict): + def __missing__(self, key): + # Return the key back in braces so it remains in the string + return "{" + key + "}" + + +def update_vds_paths(dset, modify_function): + """ + Modify the virtual paths of the specified dataset + + Note that querying the source dataspace and selection does not appear + to work (invalid pointer error from h5py) so here we assume that we're + referencing all of the source dataspace, which is correct for SWIFT + snapshots. + + dset: a h5py.Dataset object + modify_function: a function which takes the old path as its argument and + returns the new path + """ + + # Choose a temporary path for the new virtual dataset + path = dset.name + tmp_path = dset.name + ".__tmp__" + + # Build the creation property list for the new dataset + plist = h5py.h5p.create(h5py.h5p.DATASET_CREATE) + for vs in dset.virtual_sources(): + bounds = vs.vspace.get_select_bounds() + if bounds is not None: + lower, upper = bounds + size = np.asarray(upper, dtype=int) - np.asarray(lower, dtype=int) + 1 + src_space = h5py.h5s.create_simple(tuple(size)) + new_name = modify_function(vs.file_name) + plist.set_virtual( + vs.vspace, new_name.encode(), vs.dset_name.encode(), src_space + ) + + # Create the new dataset + tmp_dset = h5py.h5d.create( + dset.file["/"].id, + tmp_path.encode(), + dset.id.get_type(), + dset.id.get_space(), + dcpl=plist, + ) + tmp_dset = h5py.Dataset(tmp_dset) + for attr_name in dset.attrs: + tmp_dset.attrs[attr_name] = dset.attrs[attr_name] + + # Rename the new dataset + f = dset.file + del f[path] + f[path] = f[tmp_path] + del f[tmp_path] + + +def make_virtual_snapshot( + snapshot, + auxiliary_snapshots, + output_file, + absolute_paths=False, + discard_duplicate_datasets=False, +): + """ + Given a snapshot and auxiliary files, create + a new virtual snapshot with all datasets combine. + + snapshot: Path to the snapshot file + auxiliary_snapshots: List of auxiliary file patterns + output_file: Path to the output virtual snapshot + absolute_paths: If True, use absolute paths; if False, use relative paths + """ + + # Copy the input virtual snapshot to the output + shutil.copyfile(snapshot, output_file) + + # Open the output file + outfile = h5py.File(output_file, "r+") + + # Calculate directories for path updates + abs_snapshot_dir = os.path.abspath(os.path.dirname(snapshot)) + abs_auxiliary_dirs = [ + os.path.abspath(os.path.dirname(aux.format(file_nr=0))) + for aux in auxiliary_snapshots + ] + abs_output_dir = os.path.abspath(os.path.dirname(output_file)) + + if absolute_paths: + snapshot_dir = abs_snapshot_dir + auxiliary_dirs = abs_auxiliary_dirs + else: + snapshot_dir = os.path.relpath(abs_snapshot_dir, abs_output_dir) + auxiliary_dirs = [ + os.path.relpath(aux_dir, abs_output_dir) for aux_dir in abs_auxiliary_dirs + ] + + # Create path replacement functions + def make_replace_path(target_dir): + def replace_path(old_path): + basename = os.path.basename(old_path) + return os.path.join(target_dir, basename) + + return replace_path + + replace_snapshot_path = make_replace_path(snapshot_dir) + auxiliary_path_replacers = [make_replace_path(d) for d in auxiliary_dirs] + + all_auxiliary_datasets = {} + + for aux_index, auxiliary in enumerate(auxiliary_snapshots): + + # Check which datasets exist in the auxiliary files + # and store their attributes and datatype + filename = auxiliary.format(file_nr=0) + dset_attrs = {} + dset_dtype = {} + with h5py.File(filename, "r") as infile: + for ptype in range(7): + if not f"PartType{ptype}" in infile: + continue + dset_attrs[f"PartType{ptype}"] = {} + dset_dtype[f"PartType{ptype}"] = {} + for dset in infile[f"PartType{ptype}"].keys(): + attrs = dict(infile[f"PartType{ptype}/{dset}"].attrs) + dtype = infile[f"PartType{ptype}/{dset}"].dtype + + # Some auxiliary files are missing these attributes + if not "Value stored as physical" in attrs: + print(f"Setting comoving attrs for PartType{ptype}/{dset}") + attrs["Value stored as physical"] = [1] + attrs["Property can be converted to comoving"] = [0] + + # Add a flag that these datasets are stored in the auxiliary files + attrs["auxiliary file"] = [1] + + # Store the values we need for later + dset_attrs[f"PartType{ptype}"][dset] = attrs + dset_dtype[f"PartType{ptype}"][dset] = dtype + + # Check we don't have this dataset in any of the other auxiliary files + dset_path = f"PartType{ptype}/{dset}" + if dset_path in all_auxiliary_datasets: + other_file = all_auxiliary_datasets[f"PartType{ptype}/{dset}"] + raise ValueError( + f"{dset_path} is in {auxiliary} and {other_file}" + ) + all_auxiliary_datasets[dset_path] = auxiliary + + # Copy over the named column values, handling the case where we have + # dataset names that already exist in the original snapshot + for dset in infile.get("SubgridScheme/NamedColumns", []): + outfile_named_cols = outfile["SubgridScheme/NamedColumns"] + if dset in outfile_named_cols: + if discard_duplicate_datasets: + del outfile_named_cols[dset] + else: + outfile.move( + f"SubgridScheme/NamedColumns/{dset}", + f"SubgridScheme/NamedColumns/{dset}_snap", + ) + outfile_named_cols.create_dataset( + dset, + data=infile[f"SubgridScheme/NamedColumns/{dset}"], + ) + + # Loop over input auxiliary files to get dataset shapes + file_nr = 0 + filenames = [] + shapes = [] + counts = [] + while True: + filename = auxiliary.format(file_nr=file_nr) + if os.path.exists(filename): + filenames.append(filename) + with h5py.File(filename, "r") as infile: + shape = {} + count = {} + for ptype in range(7): + if f"PartType{ptype}" not in dset_attrs: + continue + shape[f"PartType{ptype}"] = {} + # Get the shape for each dataset + for dset in dset_attrs[f"PartType{ptype}"]: + s = infile[f"PartType{ptype}/{dset}"].shape + shape[f"PartType{ptype}"][dset] = s + # Get the number of particles in this chunk file + count[f"PartType{ptype}"] = s[0] + shapes.append(shape) + counts.append(count) + else: + break + file_nr += 1 + if file_nr == 0: + raise IOError(f"Failed to find files matching: {auxiliary}") + + # Loop over particle types in the output + for ptype in range(7): + if f"PartType{ptype}" not in dset_attrs: + continue + + # Create virtual layout for new datasets + layouts = {} + nr_parts = sum([count[f"PartType{ptype}"] for count in counts]) + for dset in dset_attrs[f"PartType{ptype}"]: + full_shape = list(shapes[0][f"PartType{ptype}"][dset]) + full_shape[0] = nr_parts + full_shape = tuple(full_shape) + dtype = dset_dtype[f"PartType{ptype}"][dset] + layouts[dset] = h5py.VirtualLayout(shape=full_shape, dtype=dtype) + + # Loop over input files + offset = 0 + for filename, count, shape in zip(filenames, counts, shapes): + n_part = count[f"PartType{ptype}"] + for dset in dset_attrs[f"PartType{ptype}"]: + layouts[dset][offset : offset + n_part] = h5py.VirtualSource( + filename, + f"PartType{ptype}/{dset}", + shape=shape[f"PartType{ptype}"][dset], + ) + offset += n_part + + # Create the virtual datasets, renaming datasets if they + # already exist in the snapshot + for dset, attrs in dset_attrs[f"PartType{ptype}"].items(): + if f"PartType{ptype}/{dset}" in outfile: + if discard_duplicate_datasets: + del outfile[f"PartType{ptype}/{dset}"] + else: + outfile.move( + f"PartType{ptype}/{dset}", f"PartType{ptype}/{dset}_snap" + ) + outfile.create_virtual_dataset( + f"PartType{ptype}/{dset}", layouts[dset], fillvalue=-999 + ) + for k, v in attrs.items(): + outfile[f"PartType{ptype}/{dset}"].attrs[k] = v + + # Update paths for this newly created auxiliary dataset + update_vds_paths( + outfile[f"PartType{ptype}/{dset}"], + auxiliary_path_replacers[aux_index], + ) + + # Copy GroupNr_bound to HaloCatalogueIndex, since + # that is the name in SOAP + if dset == "GroupNr_bound": + outfile.create_virtual_dataset( + f"PartType{ptype}/HaloCatalogueIndex", + layouts["GroupNr_bound"], + fillvalue=-999, + ) + for k, v in outfile[f"PartType{ptype}/GroupNr_bound"].attrs.items(): + outfile[f"PartType{ptype}/HaloCatalogueIndex"].attrs[k] = v + + # Update paths for HaloCatalogueIndex too + update_vds_paths( + outfile[f"PartType{ptype}/HaloCatalogueIndex"], + auxiliary_path_replacers[aux_index], + ) + + # Update paths for all original snapshot datasets + for ptype in range(7): + ptype_name = f"PartType{ptype}" + if ptype_name in outfile: + for dset_name in list(outfile[ptype_name].keys()): + dset = outfile[f"{ptype_name}/{dset_name}"] + if dset.is_virtual: + # Check if this is an auxiliary dataset (skip those, already handled) + if dset.attrs.get("auxiliary file", [0])[0] != 1: + # This is an original snapshot dataset + update_vds_paths(dset, replace_snapshot_path) + + # Done + outfile.close() + + +if __name__ == "__main__": + + import argparse + + # For description of parameters run the following: $ python make_virtual_snapshot.py --help + parser = argparse.ArgumentParser( + description=( + "Link SWIFT snapshots with SWIFT auxiliary snapshots (snapshot-like" + "files with the same number of particles in the same order as the" + "snapshot, but with less metadata), such as the SOAP memberships." + ) + ) + parser.add_argument( + "--virtual-snapshot", + type=str, + required=True, + help="Name of the SWIFT virtual snapshot file, e.g. snapshot_{snap_nr:04}.hdf5", + ) + parser.add_argument( + "--auxiliary-snapshots", + type=str, + nargs="+", + required=True, + help="One of more format strings for auxiliary files, e.g. membership_{snap_nr:04}.{file_nr}.hdf5", + ) + parser.add_argument( + "--output-file", + type=str, + required=True, + help="Name of the virtual snapshot to create, e.g. membership_{snap_nr:04}.hdf5", + ) + parser.add_argument( + "--snap-nr", + type=int, + required=False, + default=-1, + help="Snapshot number (default: -1). Not required if snap_nr is present in filenames passed.", + ) + parser.add_argument( + "--absolute-paths", + action="store_true", + help="Use absolute paths in the virtual dataset", + ) + parser.add_argument( + "--discard-duplicate-datasets", + action="store_true", + help=( + "This flag determines the behaviour when a dataset exists in both the original snapshot" + "and the auxilary file. By default the virtual file will rename the original snapshot" + "dataset as {dataset_name}_snap. If this flag is passed then the dataset from the original" + "snapshot will not be linked to" + ), + ) + args = parser.parse_args() + + print(f"Creating virtual snapshot") + for k, v in vars(args).items(): + print(f" {k}: {v}") + + # Substitute snap number + virtual_snapshot = args.virtual_snapshot.format(snap_nr=args.snap_nr) + output_file = args.output_file.format(snap_nr=args.snap_nr) + + # We don't want to replace {file_nr} for auxiliary snapshots + auxiliary_snapshots = [ + filename.format_map(SafeDict({"snap_nr": args.snap_nr})) + for filename in args.auxiliary_snapshots + ] + + # Make a new virtual snapshot with group info + make_virtual_snapshot( + virtual_snapshot, + auxiliary_snapshots, + output_file, + absolute_paths=args.absolute_paths, + discard_duplicate_datasets=args.discard_duplicate_datasets, + ) diff --git a/compression/wrong_compression.yml b/SOAP/compression/wrong_compression.yml similarity index 100% rename from compression/wrong_compression.yml rename to SOAP/compression/wrong_compression.yml diff --git a/SOAP/compute_halo_properties.py b/SOAP/compute_halo_properties.py index d9a9be12..46e72258 100644 --- a/SOAP/compute_halo_properties.py +++ b/SOAP/compute_halo_properties.py @@ -40,7 +40,6 @@ from SOAP.particle_filter.cold_dense_gas_filter import ColdDenseGasFilter from SOAP.particle_filter.recently_heated_gas_filter import RecentlyHeatedGasFilter - # Set numpy to raise divide by zero, overflow and invalid operation errors as exceptions np.seterr(divide="raise", over="raise", invalid="raise") diff --git a/SOAP/core/combine_chunks.py b/SOAP/core/combine_chunks.py index 1f969570..99c12ab8 100644 --- a/SOAP/core/combine_chunks.py +++ b/SOAP/core/combine_chunks.py @@ -420,7 +420,7 @@ def combine_chunks( fof_reg = None fof_com_unit = None fof_mass_unit = None - (fof_reg, fof_com_unit, fof_mass_unit) = comm_world.bcast( + fof_reg, fof_com_unit, fof_mass_unit = comm_world.bcast( (fof_reg, fof_com_unit, fof_mass_unit) ) diff --git a/SOAP/core/halo_tasks.py b/SOAP/core/halo_tasks.py index 53f91598..0261771c 100644 --- a/SOAP/core/halo_tasks.py +++ b/SOAP/core/halo_tasks.py @@ -10,7 +10,6 @@ from SOAP.particle_selection.halo_properties import SearchRadiusTooSmallError from SOAP.property_table import PropertyTable - # Factor by which to increase search radius when looking for density threshold SEARCH_RADIUS_FACTOR = 1.2 diff --git a/SOAP/particle_filter/recently_heated_gas_filter.py b/SOAP/particle_filter/recently_heated_gas_filter.py index 8cdb7721..63d4467d 100644 --- a/SOAP/particle_filter/recently_heated_gas_filter.py +++ b/SOAP/particle_filter/recently_heated_gas_filter.py @@ -16,7 +16,6 @@ requires knowledge of the cosmology. """ - from astropy.cosmology import w0waCDM, z_at_value import astropy.constants as const import astropy.units as astropy_units diff --git a/SOAP/particle_selection/aperture_properties.py b/SOAP/particle_selection/aperture_properties.py index 5d0cbe4b..a14c9be1 100644 --- a/SOAP/particle_selection/aperture_properties.py +++ b/SOAP/particle_selection/aperture_properties.py @@ -760,6 +760,24 @@ def TotalSNIaRate(self) -> unyt.unyt_quantity: self.star_mask_ap ].sum() + @lazy_property + def ExSituFraction(self) -> unyt.unyt_quantity: + """ + Mass fraction of bound stars that formed in a different subhalo. + """ + if self.Nstar == 0: + return None + + group_nr = self.get_dataset("PartType4/GroupNr_bound")[self.star_mask_all][ + self.star_mask_ap + ] + birth_group_nr = self.get_dataset("PartType4/BirthHaloCatalogueIndex")[ + self.star_mask_all + ][self.star_mask_ap] + ex_situ = group_nr != birth_group_nr + + return self.star_mass_fraction[ex_situ].sum() + @lazy_property def bh_mask_all(self) -> NDArray[bool]: """ @@ -3788,6 +3806,7 @@ class ApertureProperties(HaloProperty): "stellar_age_mw": False, "stellar_age_lw": False, "TotalSNIaRate": False, + "ExSituFraction": False, "HydrogenMass": False, "HeliumMass": False, "MolecularHydrogenMass": False, diff --git a/SOAP/particle_selection/subhalo_properties.py b/SOAP/particle_selection/subhalo_properties.py index cbb68649..c2294d38 100644 --- a/SOAP/particle_selection/subhalo_properties.py +++ b/SOAP/particle_selection/subhalo_properties.py @@ -481,6 +481,22 @@ def stellar_age_lw(self) -> unyt.unyt_array: Lrtot = Lr.sum() return ((Lr / Lrtot) * self.stellar_ages).sum() + @lazy_property + def ExSituFraction(self) -> unyt.unyt_quantity: + """ + Mass fraction of bound stars that formed in a different subhalo. + """ + if self.Nstar == 0: + return None + + group_nr = self.get_dataset("PartType4/GroupNr_bound")[self.star_mask_all] + birth_group_nr = self.get_dataset("PartType4/BirthHaloCatalogueIndex")[ + self.star_mask_all + ] + ex_situ = group_nr != birth_group_nr + + return self.star_mass_fraction[ex_situ].sum() + @lazy_property def bh_mask_all(self) -> NDArray[bool]: """ @@ -2455,6 +2471,7 @@ class SubhaloProperties(HaloProperty): "Lstar_luminosity_weighted", "stellar_age_mw", "stellar_age_lw", + "ExSituFraction", "Mgas_SF", "gasmetalfrac_SF", "MedianStellarBirthDensity", diff --git a/SOAP/property_table.py b/SOAP/property_table.py index 583f291f..1867517a 100644 --- a/SOAP/property_table.py +++ b/SOAP/property_table.py @@ -1905,6 +1905,22 @@ class PropertyTable: output_physical=True, a_scale_exponent=0, ), + "ExSituFraction": Property( + name="ExSituFraction", + shape=1, + dtype=np.float32, + unit="dimensionless", + description="Mass fraction of bound stars that formed in a different subhalo", + lossy_compression_filter="FMantissa9", + dmo_property=True, + particle_properties=[ + "PartType4/Masses", + "PartType4/BirthHaloCatalogueIndex", + "PartType4/GroupNr_bound", + ], + output_physical=True, + a_scale_exponent=0, + ), "Mgas": Property( name="GasMass", shape=1, @@ -5134,7 +5150,7 @@ def generate_tex_files(self, output_dir: str): # standalone table file footer tailstr = "\\end{document}" - # generate the auxilary documentation files + # generate the auxiliary documentation files with open(f"{output_dir}/timestamp.tex", "w") as ofile: ofile.write(get_version_string()) with open(f"{output_dir}/table.tex", "w") as ofile: diff --git a/compression/make_virtual_snapshot.py b/compression/make_virtual_snapshot.py deleted file mode 100644 index 5f6e37f8..00000000 --- a/compression/make_virtual_snapshot.py +++ /dev/null @@ -1,192 +0,0 @@ -#!/bin/env python - -import os.path -import h5py -import shutil - - -def make_virtual_snapshot(snapshot, membership, output_file, snap_nr): - """ - Given a FLAMINGO snapshot and group membership files, - create a new virtual snapshot with group info. - """ - - # Check which datasets exist in the membership files - # and store their attributes and datatype - filename = membership.format(file_nr=0, snap_nr=snap_nr) - dset_attrs = {} - dset_dtype = {} - with h5py.File(filename, "r") as infile: - for ptype in range(7): - if not f"PartType{ptype}" in infile: - continue - dset_attrs[f"PartType{ptype}"] = {} - dset_dtype[f"PartType{ptype}"] = {} - for dset in infile[f"PartType{ptype}"].keys(): - attrs = dict(infile[f"PartType{ptype}/{dset}"].attrs) - dtype = infile[f"PartType{ptype}/{dset}"].dtype - - # Some membership files are missing these attributes - if not "Value stored as physical" in attrs: - print(f"Setting comoving attrs for PartType{ptype}/{dset}") - attrs["Value stored as physical"] = [1] - attrs["Property can be converted to comoving"] = [0] - - # Add a flag that these are stored in the membership files - attrs["Auxilary file"] = [1] - - # Store the values we need for later - dset_attrs[f"PartType{ptype}"][dset] = attrs - dset_dtype[f"PartType{ptype}"][dset] = dtype - - # Copy the input virtual snapshot to the output - shutil.copyfile(snapshot, output_file) - - # Open the output file - outfile = h5py.File(output_file, "r+") - - # Loop over input membership files to get dataset shapes - file_nr = 0 - filenames = [] - shapes = [] - counts = [] - while True: - filename = membership.format(file_nr=file_nr, snap_nr=snap_nr) - if os.path.exists(filename): - filenames.append(filename) - with h5py.File(filename, "r") as infile: - shape = {} - count = {} - for ptype in range(7): - if f"PartType{ptype}" not in dset_attrs: - continue - shape[f"PartType{ptype}"] = {} - # Get the shape for each dataset - for dset in dset_attrs[f"PartType{ptype}"]: - s = infile[f"PartType{ptype}/{dset}"].shape - shape[f"PartType{ptype}"][dset] = s - # Get the number of particles in this chunk file - count[f"PartType{ptype}"] = s[0] - shapes.append(shape) - counts.append(count) - else: - break - file_nr += 1 - if file_nr == 0: - raise IOError(f"Failed to find files matching: {membership}") - - # Loop over particle types in the output - for ptype in range(7): - if f"PartType{ptype}" not in dset_attrs: - continue - - # Create virtual layout for new datasets - layouts = {} - nr_parts = sum([count[f"PartType{ptype}"] for count in counts]) - for dset in dset_attrs[f"PartType{ptype}"]: - full_shape = list(shapes[0][f"PartType{ptype}"][dset]) - full_shape[0] = nr_parts - full_shape = tuple(full_shape) - dtype = dset_dtype[f"PartType{ptype}"][dset] - layouts[dset] = h5py.VirtualLayout(shape=full_shape, dtype=dtype) - - # Loop over input files - offset = 0 - for filename, count, shape in zip(filenames, counts, shapes): - n_part = count[f"PartType{ptype}"] - for dset in dset_attrs[f"PartType{ptype}"]: - layouts[dset][offset : offset + n_part] = h5py.VirtualSource( - filename, - f"PartType{ptype}/{dset}", - shape=shape[f"PartType{ptype}"][dset], - ) - offset += n_part - - # Create the virtual datasets, renaming datasets if they - # already exist in the snapshot - for dset, attrs in dset_attrs[f"PartType{ptype}"].items(): - if f"PartType{ptype}/{dset}" in outfile: - outfile.move(f"PartType{ptype}/{dset}", f"PartType{ptype}/{dset}_snap") - outfile.create_virtual_dataset( - f"PartType{ptype}/{dset}", layouts[dset], fillvalue=-999 - ) - for k, v in attrs.items(): - outfile[f"PartType{ptype}/{dset}"].attrs[k] = v - - # Copy GroupNr_bound to HaloCatalogueIndex, since that is the name in SOAP - if dset == "GroupNr_bound": - outfile.create_virtual_dataset( - f"PartType{ptype}/HaloCatalogueIndex", - layouts["GroupNr_bound"], - fillvalue=-999, - ) - for k, v in outfile[f"PartType{ptype}/GroupNr_bound"].attrs.items(): - outfile[f"PartType{ptype}/HaloCatalogueIndex"].attrs[k] = v - - # Done - outfile.close() - - -if __name__ == "__main__": - - import argparse - from update_vds_paths import update_virtual_snapshot_paths - - # For description of parameters run the following: $ python make_virtual_snapshot.py --help - parser = argparse.ArgumentParser( - description=( - "Link SWIFT snapshots with SWIFT auxilary snapshots (snapshot-like" - "files with the same number of particles in the same order as the" - "snapshot, but with less metadata), such as the SOAP memberships" - ) - ) - parser.add_argument( - "virtual_snapshot", - type=str, - help="Name of the SWIFT virtual snapshot file, e.g. snapshot_{snap_nr:04}.hdf5", - ) - parser.add_argument( - "membership", - type=str, - help="Format string for membership files, e.g. membership_{snap_nr:04}.{file_nr}.hdf5", - ) - parser.add_argument( - "output_file", - type=str, - help="Name of the virtual snapshot to create, e.g. membership_{snap_nr:04}.hdf5", - ) - parser.add_argument( - "snap_nr", - type=int, - nargs="?", - default=-1, - help="Snapshot number (default: -1). Not required if snap_nr is present in filenames passed.", - ) - parser.add_argument( - "--absolute-paths", - action="store_true", - help="Use absolute paths in the virtual dataset", - ) - args = parser.parse_args() - - # Substitute snap number - virtual_snapshot = args.virtual_snapshot.format(snap_nr=args.snap_nr) - output_file = args.output_file.format(snap_nr=args.snap_nr) - - # Make a new virtual snapshot with group info - make_virtual_snapshot(virtual_snapshot, args.membership, output_file, args.snap_nr) - - # Set file paths for datasets - abs_snapshot_dir = os.path.abspath(os.path.dirname(virtual_snapshot)) - abs_membership_dir = os.path.abspath( - os.path.dirname(args.membership.format(snap_nr=args.snap_nr, file_nr=0)) - ) - if args.absolute_paths: - # Ensure all paths in the virtual file are absolute to avoid VDS prefix issues - # (we probably need to pick up datasets from two different directories) - update_virtual_snapshot_paths(output_file, abs_snapshot_dir, abs_membership_dir) - else: - abs_output_dir = os.path.abspath(os.path.dirname(output_file)) - rel_snapshot_dir = os.path.relpath(abs_snapshot_dir, abs_output_dir) - rel_membership_dir = os.path.relpath(abs_membership_dir, abs_output_dir) - update_virtual_snapshot_paths(output_file, rel_snapshot_dir, rel_membership_dir) diff --git a/compression/update_vds_paths.py b/compression/update_vds_paths.py deleted file mode 100644 index da24d57d..00000000 --- a/compression/update_vds_paths.py +++ /dev/null @@ -1,118 +0,0 @@ -#!/bin/env python - -import sys -import h5py -import numpy as np -import os.path - - -def update_vds_paths(dset, modify_function): - """ - Modify the virtual paths of the specified dataset - - Note that querying the source dataspace and selection does not appear - to work (invalid pointer error from h5py) so here we assume that we're - referencing all of the source dataspace, which is correct for SWIFT - snapshots. - - dset: a h5py.Dataset object - modify_function: a function which takes the old path as its argument and - returns the new path - """ - - # Choose a temporary path for the new virtual dataset - path = dset.name - tmp_path = dset.name + ".__tmp__" - - # Build the creation property list for the new dataset - plist = h5py.h5p.create(h5py.h5p.DATASET_CREATE) - for vs in dset.virtual_sources(): - bounds = vs.vspace.get_select_bounds() - if bounds is not None: - lower, upper = bounds - size = np.asarray(upper, dtype=int) - np.asarray(lower, dtype=int) + 1 - src_space = h5py.h5s.create_simple(tuple(size)) - new_name = modify_function(vs.file_name) - plist.set_virtual( - vs.vspace, new_name.encode(), vs.dset_name.encode(), src_space - ) - - # Create the new dataset - tmp_dset = h5py.h5d.create( - dset.file["/"].id, - tmp_path.encode(), - dset.id.get_type(), - dset.id.get_space(), - dcpl=plist, - ) - tmp_dset = h5py.Dataset(tmp_dset) - for attr_name in dset.attrs: - tmp_dset.attrs[attr_name] = dset.attrs[attr_name] - - # Rename the new dataset - f = dset.file - del f[path] - f[path] = f[tmp_path] - del f[tmp_path] - - -def update_virtual_snapshot_paths(filename, snapshot_dir=None, membership_dir=None): - """ - Add full paths to virtual datasets in the specified file - """ - f = h5py.File(filename, "r+") - - # Find all datasets in the file - all_datasets = [] - - def visit_datasets(name, obj): - if isinstance(obj, h5py.Dataset): - all_datasets.append(obj) - - f.visititems(visit_datasets) - - def replace_snapshot_path(old_path): - basename = os.path.basename(old_path) - return os.path.join(snapshot_dir, basename) - - def replace_membership_path(old_path): - basename = os.path.basename(old_path) - return os.path.join(membership_dir, basename) - - # Loop over datasets and update paths if necessary - for dset in all_datasets: - if dset.is_virtual: - name = dset.name.split("/")[-1] - # Check if the dataset comes from a membership file - if dset.attrs.get("Auxilary file", [0])[0] == 1: - if membership_dir is not None: - update_vds_paths(dset, replace_membership_path) - # Catch old datasets which didn't have the "Auxilary file" set - elif name in ( - "GroupNr_all", - "GroupNr_bound", - "Rank_bound", - "HaloCatalogueIndex", - "SpecificPotentialEnergies", - ): - if membership_dir is not None: - update_vds_paths(dset, replace_membership_path) - # Catch old case of FOF IDs from membership files - elif (name == "FOFGroupIDs") and ("PartType1/FOFGroupIDs_old" in f): - if membership_dir is not None: - update_vds_paths(dset, replace_membership_path) - # Data comes from the snapshot files - else: - if snapshot_dir is not None: - update_vds_paths(dset, replace_snapshot_path) - - f.close() - - -if __name__ == "__main__": - - filename = sys.argv[1] # Virtual snapshot file to update - snapshot_dir = sys.argv[2] # Directory with the real snapshot files - membership_dir = sys.argv[3] # Directory with the real membership files - - update_virtual_snapshot_paths(filename, snapshot_dir, membership_dir) diff --git a/format.sh b/format.sh index f532a46d..1c142b87 100755 --- a/format.sh +++ b/format.sh @@ -20,7 +20,7 @@ fi black="./black_formatting_env/bin/python3 -m black" # Make sure we don't try and format any virtual environments -files=$(echo {compression/*.py,misc/*.py,SOAP/*.py,SOAP/*/*.py,tests/*.py}) +files=$(echo {misc/*.py,SOAP/*.py,SOAP/*/*.py,tests/*.py}) # Run formatting if [[ "$1" == "--check" ]]; then diff --git a/misc/check_subhalo_ranking.py b/misc/check_subhalo_ranking.py index e10d5ad9..9b94af6f 100644 --- a/misc/check_subhalo_ranking.py +++ b/misc/check_subhalo_ranking.py @@ -3,7 +3,6 @@ import numpy as np import h5py - # Read VR IDs and positions filename = "/cosma8/data/dp004/flamingo/Runs/L1000N0900/HYDRO_FIDUCIAL_DATA/HYDRO_FIDUCIAL/VR/catalogue_0077/vr_catalogue_0077.properties.0" with h5py.File(filename, "r") as infile: diff --git a/misc/compute_BirthHaloCatalogueIndex.py b/misc/compute_BirthHaloCatalogueIndex.py new file mode 100644 index 00000000..8926db16 --- /dev/null +++ b/misc/compute_BirthHaloCatalogueIndex.py @@ -0,0 +1,285 @@ +#!/bin/env python + +""" +compute_BirthHaloCatalogueIndex.py + +This script produces an auxiliary snapshot which contains the subhalo +id each star was part of when it first formed. + +Usage: + + mpirun -- python -u misc/compute_BirthHaloCatalogueIndex.py \ + --snap-basename SNAP_BASENAME \ + --membership-basename MEMBERSHIP_BASENAME \ + --output-basename OUTPUT_FILENAME \ + --final-snap-nr FINAL_SNAP_NR + +Run "python misc/compute_BirthHaloCatalogueIndex.py -h" for a full description +of the arguments, and a list of optional arguments. + +""" + +import argparse +import datetime +import os + +from mpi4py import MPI + +comm = MPI.COMM_WORLD +comm_rank = comm.Get_rank() +comm_size = comm.Get_size() + +import h5py +import numpy as np + +import virgo.mpi.parallel_sort as psort +import virgo.mpi.parallel_hdf5 as phdf5 +from virgo.mpi.gather_array import gather_array + + +def load_particle_data(snap_basename, membership_basename, load_gas, comm): + """ + Load the particle IDs and halo membership for the particle types + we will use to match. Removes unbound particles. + """ + + particle_data = {} + + # Load particle IDs + snap_filename = snap_basename + ".{file_nr}.hdf5" + file = phdf5.MultiFile( + snap_filename, file_nr_attr=("Header", "NumFilesPerSnapshot"), comm=comm + ) + particle_data["PartType4/ParticleIDs"] = file.read("PartType4/ParticleIDs") + if load_gas: + particle_data["PartType0/ParticleIDs"] = file.read("PartType0/ParticleIDs") + + # Membership files don't have a header, so create a list of filenames + n_file = len(file.filenames) + membership_filenames = [f"{membership_basename}.{i}.hdf5" for i in range(n_file)] + # Load membership information + file = phdf5.MultiFile( + membership_filenames, file_nr_attr=("Header", "NumFilesPerSnapshot"), comm=comm + ) + particle_data["PartType4/GroupNr_bound"] = file.read("PartType4/GroupNr_bound") + if load_gas: + particle_data["PartType0/GroupNr_bound"] = file.read("PartType0/GroupNr_bound") + + # Check the two files are partitioned the same way + assert ( + particle_data["PartType4/GroupNr_bound"].shape + == particle_data["PartType4/ParticleIDs"].shape + ) + if load_gas: + assert ( + particle_data["PartType0/GroupNr_bound"].shape + == particle_data["PartType0/ParticleIDs"].shape + ) + + return particle_data + + +# Units for the dimensionless fields we will be saving +unit_attrs = { + "Conversion factor to CGS (not including cosmological corrections)": [1.0], + "Conversion factor to physical CGS (including cosmological corrections)": [1.0], + "U_I exponent": [0.0], + "U_L exponent": [0.0], + "U_M exponent": [0.0], + "U_t exponent": [0.0], + "U_T exponent": [0.0], + "a-scale exponent": [0.0], + "h-scale exponent": [0.0], + "Property can be converted to comoving": [0], + "Value stored as physical": [1], +} + + +def mpi_print(string, comm_rank): + if comm_rank == 0: + print(string) + + +if __name__ == "__main__": + + start_time = datetime.datetime.now() + + parser = argparse.ArgumentParser( + description=("Script to calculate BirthHaloCatalogueIndex of star particles"), + ) + parser.add_argument( + "--snap-basename", + type=str, + required=True, + help=( + "The basename of the snapshot files (the snapshot " + "name without the .{file_nr}.hdf5 suffix. Use " + "{snap_nr:04d} instead of the snapshot number)" + ), + ) + parser.add_argument( + "--membership-basename", + type=str, + required=True, + help="The basename of the membership files", + ) + parser.add_argument( + "--output-basename", + type=str, + required=True, + help="The basename of the output files", + ) + parser.add_argument( + "--final-snap-nr", + type=int, + required=True, + help=("Snapshot at which to load the particles"), + ) + parser.add_argument( + "--calculate-PreBirthHaloCatalogueIndex", + action="store_true", + help=( + "Whether to calculate and output the subhalo halo catalogue " + "index of the gas particle that formed each star" + ), + ) + + args = parser.parse_args() + + # Log the arguments + for k, v in vars(args).items(): + mpi_print(f" {k}: {v}", comm_rank) + + final_snap_basename = args.snap_basename.format(snap_nr=args.final_snap_nr) + final_membership_basename = args.membership_basename.format( + snap_nr=args.final_snap_nr + ) + mpi_print("Loading stars from final snapshot", comm_rank) + particle_data = load_particle_data( + final_snap_basename, + final_membership_basename, + False, + comm, + ) + star_particle_ids = particle_data["PartType4/ParticleIDs"] + star_birth_ids = particle_data["PartType4/GroupNr_bound"] + star_birth_ids[:] = -99 + star_first_snapshot = np.copy(star_birth_ids) + + if args.calculate_PreBirthHaloCatalogueIndex: + particle_data["PartType0/ParticleIDs"] = np.ones(0) + particle_data["PartType0/GroupNr_bound"] = np.ones(0) + star_prebirth_ids = np.copy(star_birth_ids) + + for snap_nr in range(0, args.final_snap_nr + 1): + + mpi_print(f"Loading data from snapshot {snap_nr}", comm_rank) + if args.calculate_PreBirthHaloCatalogueIndex: + # We need to keep the gas IDs from snapshot N-1 + gas_particle_ids = particle_data["PartType0/ParticleIDs"] + gas_group_nr = particle_data["PartType0/GroupNr_bound"] + snap_basename = args.snap_basename.format(snap_nr=snap_nr) + membership_basename = args.membership_basename.format(snap_nr=snap_nr) + particle_data = load_particle_data( + snap_basename, + membership_basename, + args.calculate_PreBirthHaloCatalogueIndex, + comm, + ) + + mpi_print(f"Matching stars", comm_rank) + # It would be quicker to make use of the BirthScaleFactors + # instead of checking all stars + idx = psort.parallel_match( + star_particle_ids[star_birth_ids == -99], + particle_data["PartType4/ParticleIDs"], + comm=comm, + ) + + new_birth_ids = psort.fetch_elements( + particle_data["PartType4/GroupNr_bound"], + idx[idx != -1], + comm=comm, + ) + + has_new_birth_id = star_birth_ids == -99 + has_new_birth_id[has_new_birth_id] = idx != -1 + star_birth_ids[has_new_birth_id] = new_birth_ids + star_first_snapshot[has_new_birth_id] = snap_nr + + if args.calculate_PreBirthHaloCatalogueIndex: + mpi_print(f"Matching gas", comm_rank) + # Identify the gas progenitor of the newly formed stars + gas_idx = psort.parallel_match( + star_particle_ids[has_new_birth_id], + gas_particle_ids, + comm=comm, + ) + # The gas progenitor may not exist for all stars due + # to particle splitting. Note this information is + # recoverable if required by using the SplitTrees. + new_prebirth_ids = -99 * np.ones_like(new_birth_ids) + new_prebirth_ids[gas_idx != -1] = psort.fetch_elements( + gas_group_nr, + gas_idx[gas_idx != -1], + comm=comm, + ) + star_prebirth_ids[has_new_birth_id] = new_prebirth_ids + + # Check we found a value for every star + assert np.sum(star_birth_ids == -99) == 0 + + # Set up what we want to output + output = { + "BirthHaloCatalogueIndex": star_birth_ids, + "FirstSnapshot": star_first_snapshot, + } + attrs = { + "BirthHaloCatalogueIndex": { + "Description": "The HaloCatalogueIndex of this particle at the first snapshot it appeared." + }, + "FirstSnapshot": { + "Description": "Index of the first simulation snapshot in which the star particle is present." + }, + } + attrs["BirthHaloCatalogueIndex"].update(unit_attrs) + attrs["FirstSnapshot"].update(unit_attrs) + if args.calculate_PreBirthHaloCatalogueIndex: + output["PreBirthHaloCatalogueIndex"] = star_prebirth_ids + attrs["PreBirthHaloCatalogueIndex"] = { + "Description": "The HaloCatalogueIndex of gas prognitor at the snapshot before the star formed. -99 if no gas progenitor is found." + } + attrs["PreBirthHaloCatalogueIndex"].update(unit_attrs) + + # Check the output directory exists + output_filename = ( + args.output_basename.format(snap_nr=args.final_snap_nr) + ".{file_nr}.hdf5" + ) + if comm_rank == 0: + output_dir = os.path.dirname(output_filename) + os.makedirs(output_dir, exist_ok=True) + comm.barrier() + + # Write the output + mpi_print("Writing output", comm_rank) + snap_file = phdf5.MultiFile( + final_snap_basename + ".{file_nr}.hdf5", + file_nr_attr=("Header", "NumFilesPerSnapshot"), + comm=comm, + ) + elements_per_file = snap_file.get_elements_per_file( + "ParticleIDs", group="PartType4" + ) + snap_file.write( + output, + elements_per_file, + filenames=output_filename, + mode="w", + group="PartType4", + attrs=attrs, + ) + + # Finished + comm.barrier() + mpi_print(f"Runtime: {datetime.datetime.now() - start_time}", comm_rank) + mpi_print("Done!", comm_rank) diff --git a/parameter_files/COLIBRE_HYBRID.yml b/parameter_files/COLIBRE_HYBRID.yml index 0f9b3f6e..4f823849 100644 --- a/parameter_files/COLIBRE_HYBRID.yml +++ b/parameter_files/COLIBRE_HYBRID.yml @@ -244,6 +244,7 @@ ApertureProperties: StellarCylindricalVelocityDispersionLuminosityWeighted: false StellarCylindricalVelocityDispersionVerticalLuminosityWeighted: false StellarCylindricalVelocityDispersionDiscPlaneLuminosityWeighted: false + ExSituFraction: false TotalMass: true TotalSNIaRate: true GasMassInColdDenseDiffuseMetals: @@ -742,6 +743,7 @@ SubhaloProperties: StellarCylindricalVelocityDispersionLuminosityWeighted: false StellarCylindricalVelocityDispersionVerticalLuminosityWeighted: false StellarCylindricalVelocityDispersionDiscPlaneLuminosityWeighted: false + ExSituFraction: false aliases: PartType0/LastSNIIKineticFeedbackDensities: PartType0/DensitiesAtLastSupernovaEvent PartType0/LastSNIIThermalFeedbackDensities: PartType0/DensitiesAtLastSupernovaEvent diff --git a/parameter_files/COLIBRE_THERMAL.yml b/parameter_files/COLIBRE_THERMAL.yml index 403267a5..7e2eb418 100644 --- a/parameter_files/COLIBRE_THERMAL.yml +++ b/parameter_files/COLIBRE_THERMAL.yml @@ -244,6 +244,7 @@ ApertureProperties: StellarCylindricalVelocityDispersionDiscPlaneLuminosityWeighted: false StellarRotationalVelocity: false StellarRotationalVelocityLuminosityWeighted: false + ExSituFraction: false TotalMass: true TotalSNIaRate: true GasMassInColdDenseDiffuseMetals: @@ -742,6 +743,7 @@ SubhaloProperties: StellarCylindricalVelocityDispersionLuminosityWeighted: false StellarCylindricalVelocityDispersionVerticalLuminosityWeighted: false StellarCylindricalVelocityDispersionDiscPlaneLuminosityWeighted: false + ExSituFraction: false aliases: PartType0/LastSNIIKineticFeedbackDensities: PartType0/DensitiesAtLastSupernovaEvent PartType0/LastSNIIThermalFeedbackDensities: PartType0/DensitiesAtLastSupernovaEvent diff --git a/scripts/COLIBRE/compress_group_membership.sh b/scripts/COLIBRE/compress_group_membership.sh index ee34965f..f04ec4c8 100644 --- a/scripts/COLIBRE/compress_group_membership.sh +++ b/scripts/COLIBRE/compress_group_membership.sh @@ -90,7 +90,10 @@ echo "Creating virtual snapshot" snapshot="${output_dir}/${sim}/snapshots/colibre_${snapnum}/colibre_${snapnum}.hdf5" membership="${output_filename}.{file_nr}.hdf5" virtual="${outbase}/colibre_with_SOAP_membership_${snapnum}.hdf5" -python compression/make_virtual_snapshot.py $snapshot $membership $virtual +python SOAP/compression/make_virtual_snapshot.py \ + --virtual-snapshot $snapshot \ + --auxiliary-snapshots $membership \ + --output-file $virtual echo "Setting virtual file to be read-only" chmod a=r "${virtual}" diff --git a/scripts/COLIBRE/compress_halo_properties.sh b/scripts/COLIBRE/compress_halo_properties.sh index b75a7632..5acbca36 100755 --- a/scripts/COLIBRE/compress_halo_properties.sh +++ b/scripts/COLIBRE/compress_halo_properties.sh @@ -29,7 +29,7 @@ output_dir="/cosma8/data/dp004/dc-mcgi1/COLIBRE/Runs" scratch_dir="/snap8/scratch/dp004/dc-mcgi1/COLIBRE/Runs" # compression script -script="./compression/compress_soap_catalogue.py" +script="./SOAP/compression/compress_soap_catalogue.py" # Which snapshot to do snapnum=`printf '%04d' ${SLURM_ARRAY_TASK_ID}` @@ -41,7 +41,7 @@ sim="${SLURM_JOB_NAME}" input_filename="${input_dir}/${sim}/SOAP_uncompressed/halo_properties_${snapnum}.hdf5" # Location and name of the output SOAP catalogue -outbase="${output_dir}/${sim}/SOAP" +outbase="${output_dir}/${sim}/SOAP-ExSitu" mkdir -p $outbase output_filename="${outbase}/halo_properties_${snapnum}.hdf5" diff --git a/scripts/COLIBRE/compute_birth_index.sh b/scripts/COLIBRE/compute_birth_index.sh new file mode 100755 index 00000000..35181cbb --- /dev/null +++ b/scripts/COLIBRE/compute_birth_index.sh @@ -0,0 +1,47 @@ +#!/bin/bash -l + +#SBATCH --cpus-per-task=1 +#SBATCH -o ./logs/birth_track_id_%j.out +#SBATCH -p cosma8 +#SBATCH -A dp004 +#SBATCH -J BirthHaloCatalogueIndex +#SBATCH --nodes=4 +#SBATCH -t 24:00:00 +# N0752: 1 node, 2 hours +# N1504: 1 node, 12 hours +# N3008: 4 nodes, 24 hours + +set -e + +# TODO: Set these values +base_dir="/cosma8/data/dp004/dc-mcgi1/COLIBRE/BirthHaloCatalogueIndex" +output_dir="/cosma8/data/dp004/dc-mcgi1/COLIBRE/BirthHaloCatalogueIndex" +sim="L0400N3008/Thermal" +snapnum="0127" + +snap_basename="${base_dir}/${sim}/snapshots/colibre_{snap_nr:04d}/colibre_{snap_nr:04d}" +membership_basename="${base_dir}/${sim}/SOAP-HBT/membership_{snap_nr:04d}/membership_{snap_nr:04d}" +output_basename="${output_dir}/${sim}/SOAP-ExSitu/birth_${snapnum}/birth_${snapnum}" + +mpirun -- python misc/compute_BirthHaloCatalogueIndex.py \ + --snap-basename ${snap_basename} \ + --membership-basename ${membership_basename} \ + --output-basename ${output_basename} \ + --final-snap-nr ${snapnum} \ + --calculate-PreBirthHaloCatalogueIndex + +chmod a=r "${output_basename}"* + +snapshot="${snap_basename}.hdf5" +membership="${membership_basename}.{file_nr}.hdf5" +output="${output_basename}.{file_nr}.hdf5" +virtual="${output_dir}/${sim}/SOAP-ExSitu/birth_${snapnum}.hdf5" +python SOAP/compression/make_virtual_snapshot.py \ + --virtual-snapshot "$snapshot" \ + --auxiliary-snapshots "$membership" "$output" \ + --output-file "$virtual" \ + --snap-nr "$snapnum" \ + +chmod a=r "${virtual}" + +echo "Job complete!" diff --git a/scripts/EAGLE.sh b/scripts/EAGLE.sh index 8f87ef72..9c83e6f4 100755 --- a/scripts/EAGLE.sh +++ b/scripts/EAGLE.sh @@ -77,10 +77,10 @@ cd "${output_dir}/swift_snapshots/swift_${snap_nr}" python "${soap_dir}/create_virtual_snapshot.py" "snap_${snap_nr}.0.hdf5" cd - -python compression/make_virtual_snapshot.py \ - "${output_dir}/swift_snapshots/swift_${snap_nr}/snap_${snap_nr}.hdf5" \ - "${output_dir}/SOAP_uncompressed/membership_${snap_nr}/membership_${snap_nr}.{file_nr}.hdf5" \ - "${output_dir}/SOAP_uncompressed/snap_${snap_nr}.hdf5" \ +python SOAP/compression/make_virtual_snapshot.py \ + --virtual-snapshot "${output_dir}/swift_snapshots/swift_${snap_nr}/snap_${snap_nr}.hdf5" \ + --auxiliary-snapshots "${output_dir}/SOAP_uncompressed/membership_${snap_nr}/membership_${snap_nr}.{file_nr}.hdf5" \ + --output-file "${output_dir}/SOAP_uncompressed/snap_${snap_nr}.hdf5" ######### Run SOAP diff --git a/scripts/FLAMINGO/L1000N1800/compress_halo_properties_L1000N1800.sh b/scripts/FLAMINGO/L1000N1800/compress_halo_properties_L1000N1800.sh index 60d202d0..ea032de7 100644 --- a/scripts/FLAMINGO/L1000N1800/compress_halo_properties_L1000N1800.sh +++ b/scripts/FLAMINGO/L1000N1800/compress_halo_properties_L1000N1800.sh @@ -29,7 +29,7 @@ output_dir="/cosma8/data/dp004/dc-mcgi1/FLAMINGO/Runs" scratch_dir="/snap8/scratch/dp004/dc-mcgi1/FLAMINGO/Runs" # compression script -script="./compression/compress_soap_catalogue.py" +script="./SOAP/compression/compress_soap_catalogue.py" # Which snapshot to do snapnum=`printf '%04d' ${SLURM_ARRAY_TASK_ID}` diff --git a/tests/COLIBRE/run_L0025N0188_Thermal.sh b/tests/COLIBRE/run_L0025N0188_Thermal.sh index 27a43e1e..58969b78 100755 --- a/tests/COLIBRE/run_L0025N0188_Thermal.sh +++ b/tests/COLIBRE/run_L0025N0188_Thermal.sh @@ -30,7 +30,7 @@ python tests/COLIBRE/create_parameters_file.py rm -r output/SOAP-tmp # Run SOAP on eight cores processing the selected halos. Use 'python3 -m pdb' to start in the debugger. -mpirun -np 8 python3 -u -m mpi4py SOAP/compute_halo_properties.py \ +mpirun -np 8 python SOAP/compute_halo_properties.py \ ./tests/COLIBRE/test_parameters.yml \ --halo-indices ${halo_indices} \ --sim-name=${sim} --snap-nr=${snapnum} --chunks=1 diff --git a/tests/test_SO_properties.py b/tests/test_SO_properties.py index 1ba8c03c..bd4459a4 100644 --- a/tests/test_SO_properties.py +++ b/tests/test_SO_properties.py @@ -419,7 +419,7 @@ def calculate_SO_properties_nfw_halo(seed, num_part, c): "crit", ) - (input_halo, data, rmax, Mtot, Npart, particle_numbers) = dummy_halos.gen_nfw_halo( + input_halo, data, rmax, Mtot, Npart, particle_numbers = dummy_halos.gen_nfw_halo( 100, c, num_part ) From d50455a6f35834490a7f87c7042365e8f74ea597 Mon Sep 17 00:00:00 2001 From: robjmcgibbon <36136863+robjmcgibbon@users.noreply.github.com> Date: Wed, 18 Mar 2026 12:01:35 +0000 Subject: [PATCH 10/19] Bump minor version --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 50f5b9c1..3704b2a2 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -6,7 +6,7 @@ build-backend = "setuptools.build_meta" [project] name = "SOAP" -version = "0.1" +version = "0.1.1" description = "MPI parallel Python code to compute properties of halos in SWIFT n-body simulations" readme = "README.md" requires-python = ">=3.10" From b57beeff44be6d4bbca8869031fedc81d5688054 Mon Sep 17 00:00:00 2001 From: robjmcgibbon <36136863+robjmcgibbon@users.noreply.github.com> Date: Mon, 30 Mar 2026 12:07:55 +0100 Subject: [PATCH 11/19] Add citation files (#210) * Add citation files * Update readme * Update readme * Add ads link --- CITATION.cff | 37 +++++++++++++++++++++++++++++++++++++ README.md | 4 +++- 2 files changed, 40 insertions(+), 1 deletion(-) create mode 100644 CITATION.cff diff --git a/CITATION.cff b/CITATION.cff new file mode 100644 index 00000000..b1e9dfc4 --- /dev/null +++ b/CITATION.cff @@ -0,0 +1,37 @@ +cff-version: 1.2.0 +message: "If you use this software, please cite it as below." +authors: + - family-names: "McGibbon" + given-names: "Robert" + - family-names: "Helly" + given-names: "John" + - family-names: "Schaye" + given-names: "Joop" + - family-names: "Schaller" + given-names: "Matthieu" + - family-names: "Vandenbroucke" + given-names: "Bert" +title: "SOAP: A Python Package for Calculating the Properties of Galaxies and Halos Formed in Cosmological Simulations" +doi: "10.21105/joss.08252" +date-released: "2025-07-01" +preferred-citation: + type: article + authors: + - family-names: "McGibbon" + given-names: "Robert" + - family-names: "Helly" + given-names: "John" + - family-names: "Schaye" + given-names: "Joop" + - family-names: "Schaller" + given-names: "Matthieu" + - family-names: "Vandenbroucke" + given-names: "Bert" + doi: "10.21105/joss.08252" + journal: "The Journal of Open Source Software" + month: 7 + year: 2025 + title: "SOAP: A Python Package for Calculating the Properties of Galaxies and Halos Formed in Cosmological Simulations" + volume: 10 + issue: 111 + start: 8252 diff --git a/README.md b/README.md index 7b72e845..73148e35 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # SOAP: Spherical Overdensity and Aperture Processor -[![DOI](https://joss.theoj.org/papers/10.21105/joss.08252/status.svg)](https://doi.org/10.21105/joss.08252) +[![DOI](https://joss.theoj.org/papers/10.21105/joss.08252/status.svg)](https://doi.org/10.21105/joss.08252) This repository contains programs which can be used to compute properties of halos in spherical apertures in [SWIFT](https://swift.strw.leidenuniv.nl/) snapshots. @@ -8,6 +8,8 @@ The resulting output halo catalogues can be read using the [swiftsimio](https://swiftsimio.readthedocs.io/en/latest/) python package. +Please cite SOAP using the [JOSS paper](https://ui.adsabs.harvard.edu/abs/2025JOSS...10.8252M) + ## Installation The code is written in python and uses mpi4py for parallelism. From b9b09afd62053f4f83d4e58b4c95bc31f0bb2d8e Mon Sep 17 00:00:00 2001 From: robjmcgibbon <36136863+robjmcgibbon@users.noreply.github.com> Date: Wed, 1 Apr 2026 13:26:10 +0100 Subject: [PATCH 12/19] Update metallicity descriptions (#211) --- SOAP/property_table.py | 48 +++++++++++++++++++++--------------------- 1 file changed, 24 insertions(+), 24 deletions(-) diff --git a/SOAP/property_table.py b/SOAP/property_table.py index 1867517a..d86fdfe0 100644 --- a/SOAP/property_table.py +++ b/SOAP/property_table.py @@ -4095,7 +4095,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="dimensionless", - description="Linear sum of the oxygen over hydrogen ratio of gas, multiplied with the gas mass.", + description="Linear sum of the oxygen over hydrogen ratio of gas, weighted by with the cold dense gas mass.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=[ @@ -4112,7 +4112,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="dimensionless", - description="Linear sum of the total nitrogen over oxygen ratio of gas, multiplied with the gas mass.", + description="Linear sum of the total nitrogen over oxygen ratio of gas, weighted by with the cold dense gas mass.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=[ @@ -4129,7 +4129,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="dimensionless", - description="Linear sum of the total carbon over oxygen ratio of gas, multiplied with the gas mass.", + description="Linear sum of the total carbon over oxygen ratio of gas, weighted by with the cold dense gas mass.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=[ @@ -4146,7 +4146,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="dimensionless", - description="Linear sum of the diffuse nitrogen over oxygen ratio of gas, multiplied with the gas mass.", + description="Linear sum of the diffuse nitrogen over oxygen ratio of gas, weighted by with the cold dense gas mass.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=[ @@ -4163,7 +4163,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="dimensionless", - description="Linear sum of the diffuse carbon over oxygen ratio of gas, multiplied with the gas mass.", + description="Linear sum of the diffuse carbon over oxygen ratio of gas, weighted by with the cold dense gas mass.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=[ @@ -4180,7 +4180,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="dimensionless", - description="Linear sum of the diffuse oxygen over hydrogen ratio of gas, multiplied with the gas mass.", + description="Linear sum of the diffuse oxygen over hydrogen ratio of gas, weighted by with the cold dense gas mass.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=[ @@ -4197,7 +4197,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="dimensionless", - description="Logarithmic sum of the diffuse nitrogen over oxygen ratio of gas, multiplied with the gas mass. Imposes a lower limit of 1.e-4 times solar N/O.", + description="Logarithmic sum of the diffuse nitrogen over oxygen ratio of gas, weighted by with the cold dense gas mass. Imposes a lower limit of 1.e-4 times solar N/O.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=[ @@ -4214,7 +4214,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="dimensionless", - description="Logarithmic sum of the diffuse nitrogen over oxygen ratio of gas, multiplied with the gas mass. Imposes a lower limit of 1.e-3 times solar N/O.", + description="Logarithmic sum of the diffuse nitrogen over oxygen ratio of gas, weighted by with the cold dense gas mass. Imposes a lower limit of 1.e-3 times solar N/O.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=[ @@ -4231,7 +4231,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="dimensionless", - description="Logarithmic sum of the diffuse carbon over oxygen ratio of gas, multiplied with the gas mass. Imposes a lower limit of 1.e-4 times solar C/O.", + description="Logarithmic sum of the diffuse carbon over oxygen ratio of gas, weighted by with the cold dense gas mass. Imposes a lower limit of 1.e-4 times solar C/O.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=[ @@ -4248,7 +4248,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="dimensionless", - description="Logarithmic sum of the diffuse carbon over oxygen ratio of gas, multiplied with the gas mass. Imposes a lower limit of 1.e-3 times solar C/O.", + description="Logarithmic sum of the diffuse carbon over oxygen ratio of gas, weighted by with the cold dense gas mass. Imposes a lower limit of 1.e-3 times solar C/O.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=[ @@ -4265,7 +4265,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="dimensionless", - description="Logarithmic sum of the diffuse oxygen over hydrogen ratio of gas, multiplied with the gas mass. Imposes a lower limit of 1.e-4 times solar O/H.", + description="Logarithmic sum of the diffuse oxygen over hydrogen ratio of gas, weighted by with the cold dense gas mass. Imposes a lower limit of 1.e-4 times solar O/H.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=[ @@ -4282,7 +4282,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="dimensionless", - description="Logarithmic sum of the diffuse oxygen over hydrogen ratio of gas, multiplied with the gas mass. Imposes a lower limit of 1.e-3 times solar O/H.", + description="Logarithmic sum of the diffuse oxygen over hydrogen ratio of gas, weighted by with the cold dense gas mass. Imposes a lower limit of 1.e-3 times solar O/H.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=[ @@ -4299,7 +4299,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="dimensionless", - description="Logarithmic sum of the diffuse oxygen over hydrogen ratio of atomic gas, multiplied with the gas mass. Imposes a lower limit of 1.e-4 times solar O/H.", + description="Logarithmic sum of the diffuse oxygen over hydrogen ratio of atomic gas, weighted by the atomic gas mass. Imposes a lower limit of 1.e-4 times solar O/H.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=[ @@ -4318,7 +4318,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="dimensionless", - description="Logarithmic sum of the diffuse oxygen over hydrogen ratio of atomic gas, multiplied with the gas mass. Imposes a lower limit of 1.e-3 times solar O/H.", + description="Logarithmic sum of the diffuse oxygen over hydrogen ratio of atomic gas, weighted by the atomic gas mass. Imposes a lower limit of 1.e-3 times solar O/H.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=[ @@ -4337,7 +4337,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="dimensionless", - description="Logarithmic sum of the diffuse oxygen over hydrogen ratio of molecular gas, multiplied with the gas mass. Imposes a lower limit of 1.e-4 times solar O/H.", + description="Logarithmic sum of the diffuse oxygen over hydrogen ratio of molecular gas, weighted by the molecular gas mass. Imposes a lower limit of 1.e-4 times solar O/H.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=[ @@ -4356,7 +4356,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="dimensionless", - description="Logarithmic sum of the diffuse oxygen over hydrogen ratio of molecular gas, multiplied with the gas mass. Imposes a lower limit of 1.e-3 times solar O/H.", + description="Logarithmic sum of the diffuse oxygen over hydrogen ratio of molecular gas, weighted by the molecular gas mass. Imposes a lower limit of 1.e-3 times solar O/H.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=[ @@ -4375,7 +4375,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="dimensionless", - description="Linear sum of the iron over hydrogen ratio of stars, multiplied with the stellar mass.", + description="Linear sum of the iron over hydrogen ratio of stars, weighted by the stellar mass.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=["PartType4/Masses", "PartType4/ElementMassFractions"], @@ -4387,7 +4387,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="dimensionless", - description="Logarithmic sum of the iron over hydrogen ratio of stars, multiplied with the stellar mass. Imposes a lower limit of 1.e-4 times solar Fe/H.", + description="Logarithmic sum of the iron over hydrogen ratio of stars, weighted by the stellar mass. Imposes a lower limit of 1.e-4 times solar Fe/H.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=["PartType4/Masses", "PartType4/ElementMassFractions"], @@ -4399,7 +4399,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="dimensionless", - description="Logarithmic sum of the iron over hydrogen ratio of stars, multiplied with the stellar mass. Imposes a lower limit of 1.e-3 times solar Fe/H.", + description="Logarithmic sum of the iron over hydrogen ratio of stars, weighted by the stellar mass. Imposes a lower limit of 1.e-3 times solar Fe/H.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=["PartType4/Masses", "PartType4/ElementMassFractions"], @@ -4411,7 +4411,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="dimensionless", - description="Linear sum of the magnesium over hydrogen ratio of stars, multiplied with the stellar mass.", + description="Linear sum of the magnesium over hydrogen ratio of stars, weighted by the stellar mass.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=["PartType4/Masses", "PartType4/ElementMassFractions"], @@ -4423,7 +4423,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="dimensionless", - description="Logarithmic sum of the magnesium over hydrogen ratio of stars, multiplied with the stellar mass. Imposes a lower limit of 1.e-4 times solar Fe/H.", + description="Logarithmic sum of the magnesium over hydrogen ratio of stars, weighted by the stellar mass. Imposes a lower limit of 1.e-4 times solar Fe/H.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=["PartType4/Masses", "PartType4/ElementMassFractions"], @@ -4435,7 +4435,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="dimensionless", - description="Logarithmic sum of the magnesium over hydrogen ratio of stars, multiplied with the stellar mass. Imposes a lower limit of 1.e-3 times solar Fe/H.", + description="Logarithmic sum of the magnesium over hydrogen ratio of stars, weighted by the stellar mass. Imposes a lower limit of 1.e-3 times solar Fe/H.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=["PartType4/Masses", "PartType4/ElementMassFractions"], @@ -4465,7 +4465,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="dimensionless", - description="Logarithmic sum of the iron over hydrogen ratio of stars, multiplied with the stellar mass, where only iron from SNIa is included. Imposes a lower limit of 1.e-4 times solar Fe/H.", + description="Logarithmic sum of the iron over hydrogen ratio of stars, weighted by the stellar mass, where only iron from SNIa is included. Imposes a lower limit of 1.e-4 times solar Fe/H.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=[ @@ -4481,7 +4481,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="dimensionless", - description="Sum of the iron over hydrogen ratio of stars, multiplied with the stellar mass, where only iron from SNIa is included.", + description="Sum of the iron over hydrogen ratio of stars, weighted by the stellar mass, where only iron from SNIa is included.", lossy_compression_filter="FMantissa9", dmo_property=False, particle_properties=[ From 7006a49e15ac3dc1118e754f2d9917662f52c077 Mon Sep 17 00:00:00 2001 From: robjmcgibbon <36136863+robjmcgibbon@users.noreply.github.com> Date: Thu, 9 Apr 2026 16:54:35 +0100 Subject: [PATCH 13/19] Update EAGLE script for WDM run (#212) * Update EAGLE script for WDM run * Format --- misc/convert_eagle.py | 23 ++++++++++++++++++----- scripts/EAGLE.sh | 6 +++++- 2 files changed, 23 insertions(+), 6 deletions(-) diff --git a/misc/convert_eagle.py b/misc/convert_eagle.py index 6066e2bc..3f818361 100644 --- a/misc/convert_eagle.py +++ b/misc/convert_eagle.py @@ -212,6 +212,7 @@ G = const_internal_header["newton_G"][0] critical_density = 3 * (H**2) / (8 * np.pi * G) cosmology_header = { + "Cosmological run": 1, "Omega_b": header.attrs["OmegaBaryon"], "Omega_m": header.attrs["Omega0"], "Omega_k": 0, @@ -322,6 +323,13 @@ "description": None, "conversion_factor": None, }, + "SmoothingLength": { + "swift_name": "SmoothingLengths", + "exponents": {"L": 1, "M": 0, "T": 0, "t": 0}, + "a_exponent": None, + "description": None, + "conversion_factor": None, + }, "StarFormationRate": { "swift_name": "StarFormationRates", "exponents": {"L": 0, "M": 1, "T": 0, "t": -1}, @@ -560,10 +568,15 @@ "conversion_factor" ] = conversion_factor - # Load DM mass from mass table + # DM mass can be a special case if "Mass" in properties.get(f"PartType1", {}): - dm_mass = infile["Header"].attrs["MassTable"][1] / h - properties["PartType1"]["Mass"]["conversion_factor"] = dm_mass + if "Mass" not in infile["PartType1"]: + # Load DM mass from mass table + dm_mass = infile["Header"].attrs["MassTable"][1] / h + properties["PartType1"]["Mass"]["conversion_factor"] = dm_mass + else: + # Treat DM mass as any other property + dm_mass = 0 # Get list of elements for ElementMassFractions if "ElementMassFractions" in properties.get(f"PartType0", {}): @@ -661,8 +674,8 @@ if comm_rank == 0: print(f"Converting PartType{ptype}/{prop}") - # DM particles all have the same mass, so are not saved in the snapshots - if (ptype == 1) and (prop == "Mass"): + if (ptype == 1) and (prop == "Mass") and (dm_mass == 0): + # DM particles all have the same mass, so are not saved in the snapshots arr = np.ones(pos.shape[0]) else: # Load data from file and sort according to cell structure diff --git a/scripts/EAGLE.sh b/scripts/EAGLE.sh index 9c83e6f4..72bc0511 100755 --- a/scripts/EAGLE.sh +++ b/scripts/EAGLE.sh @@ -15,8 +15,10 @@ # Install Hdecompose with: # pip install git+ssh://git@github.com/kyleaoman/Hdecompose.git # -# Download virtual snapshot script with: +# Download virtual snapshot script (to the main SOAP directory) with: # wget https://gitlab.cosma.dur.ac.uk/swift/swiftsim/-/raw/master/tools/create_virtual_snapshot.py +# +# Submit this script from the main SOAP directory ($ sbatch scripts/EAGLE.sh) sim_name='L0100N1504' snap_nr="028" @@ -29,6 +31,8 @@ source openmpi-5.0.3-hdf5-1.12.3-env/bin/activate ######## Link files to snap (to remove awful z suffix) sim_dir="/cosma7/data/Eagle/ScienceRuns/Planck1/${sim_name}/PE/REFERENCE/data" +# Note that if you update the output directory, you will also need to +# update the SOAP parameter file output_dir="/snap7/scratch/dp004/dc-mcgi1/SOAP_EAGLE/${sim_name}" sim_snap_dir="${sim_dir}/particledata_${snap_nr}_${z_suffix}" From 657286944e426e291f873b794c199734cd28f00f Mon Sep 17 00:00:00 2001 From: robjmcgibbon <36136863+robjmcgibbon@users.noreply.github.com> Date: Mon, 20 Apr 2026 12:40:45 +0100 Subject: [PATCH 14/19] Place large subhalos on their own chunk (#190) * Working? version * Add more options * Delete unused param file * Missing import * Set separate chunk thresholds * Update subfind radius --- SOAP/catalogue_readers/read_subfind.py | 4 +- SOAP/core/chunk_tasks.py | 2 + SOAP/core/combine_chunks.py | 1 + SOAP/core/domain_decomposition.py | 91 ++++++++++++++- SOAP/core/halo_centres.py | 39 ++++--- SOAP/core/soap_args.py | 20 +++- SOAP/core/swift_cells.py | 14 ++- documentation/SOAP.tex | 8 +- parameter_files/COLIBRE_HYBRID.yml | 7 ++ parameter_files/COLIBRE_THERMAL.yml | 7 ++ parameter_files/README.md | 9 ++ parameter_files/XRAY_ONLY.yml | 147 ------------------------- 12 files changed, 168 insertions(+), 181 deletions(-) delete mode 100644 parameter_files/XRAY_ONLY.yml diff --git a/SOAP/catalogue_readers/read_subfind.py b/SOAP/catalogue_readers/read_subfind.py index b6074dd2..14b5d10f 100644 --- a/SOAP/catalogue_readers/read_subfind.py +++ b/SOAP/catalogue_readers/read_subfind.py @@ -201,7 +201,7 @@ def read_gadget4_catalogue(comm, basename, a_unit, registry, boxsize): length_conversion = (gadget_length_unit / swift_pmpc).to(unyt.dimensionless) # Get position in comoving Mpc, assuming input position from Gadget is comoving - cofp = data["Subhalo/SubhaloPos"] * length_conversion * swift_cmpc + cofp = (data["Subhalo/SubhaloPos"] * length_conversion * swift_cmpc) % boxsize # Store central halo flag is_central = np.where(data["Subhalo/SubhaloRankInGr"] == 0, 1, 0) @@ -227,7 +227,7 @@ def read_gadget4_catalogue(comm, basename, a_unit, registry, boxsize): # Store initial search radius search_radius = ( - data["Subhalo/SubhaloHalfmassRad"] * length_conversion * swift_pmpc + (4 * data["Subhalo/SubhaloHalfmassRad"]) * length_conversion * swift_pmpc ) # different units from cofm, not a typo! local_halo = { diff --git a/SOAP/core/chunk_tasks.py b/SOAP/core/chunk_tasks.py index 413c4088..ae6e7fb3 100644 --- a/SOAP/core/chunk_tasks.py +++ b/SOAP/core/chunk_tasks.py @@ -381,6 +381,7 @@ def message(m): ) as outfile: results.collective_write(outfile, comm_have_results) comm_have_results.Free() + comm.barrier() # Store time taken for this task timings.append(task_time_all_iterations) @@ -394,6 +395,7 @@ def message(m): calc_names = sorted([hp.name for hp in self.halo_prop_list]) outfile.attrs["calc_names"] = calc_names outfile.attrs["Write complete"] = True + comm.barrier() # Return the names, dimensions and units of the quantities we computed # so that we can check they're consistent between chunks diff --git a/SOAP/core/combine_chunks.py b/SOAP/core/combine_chunks.py index 99c12ab8..63079d69 100644 --- a/SOAP/core/combine_chunks.py +++ b/SOAP/core/combine_chunks.py @@ -259,6 +259,7 @@ def combine_chunks( "Redshift", "RunName", "Scale-factor", + "Time", ]: header.attrs[attr] = cellgrid.swift_header_group[attr] header.attrs["Code"] = "SOAP" diff --git a/SOAP/core/domain_decomposition.py b/SOAP/core/domain_decomposition.py index d20f64d0..47d59d95 100644 --- a/SOAP/core/domain_decomposition.py +++ b/SOAP/core/domain_decomposition.py @@ -3,9 +3,10 @@ import numpy as np import virgo.util.peano as peano import virgo.mpi.parallel_sort as psort +from virgo.mpi.gather_array import gather_array -def peano_decomposition(boxsize, local_halo, nr_chunks, comm): +def peano_decomposition(boxsize, local_halo, nr_chunks, comm, separate_chunks): """ Gadget style domain decomposition using Peano-Hilbert curve. Allows an arbitrary number of chunks and tries to put equal @@ -23,6 +24,42 @@ def peano_decomposition(boxsize, local_halo, nr_chunks, comm): """ comm_rank = comm.Get_rank() + comm_size = comm.Get_size() + + # Handle halos which should be on their own chunk + # for memory reasons + nr_large_halo = 0 + total_nr_halos = comm.allreduce(local_halo["index"].shape[0]) + if separate_chunks: + nr_local_halos = local_halo["index"].shape[0] + local_halo_offset = comm.scan(nr_local_halos) - nr_local_halos + + # Identify the large halos + # Note that separate_chunks is sorted by 'n_bound_threshold' + n_bound_threshold = separate_chunks[-1]["n_bound_threshold"] + large_halo_mask = local_halo["nr_bound_part"] > n_bound_threshold + idx_large_halo = np.where(large_halo_mask)[0] + local_halo_offset + idx_large_halo = gather_array(idx_large_halo, root=comm_size - 1) + if comm_rank == comm_size - 1: + nr_large_halo = idx_large_halo.shape[0] + else: + idx_large_halo = np.zeros(0, dtype=np.int32) + nr_large_halo = comm.bcast(nr_large_halo, root=comm_size - 1) + + if nr_large_halo > 0: + if comm_rank == 0: + print(f"Placing {nr_large_halo} halos in custom chunks") + + # Move the data for these halos to the final rank + large_halo = {} + for name in local_halo: + large_halo[name] = psort.fetch_elements( + local_halo[name], idx_large_halo, comm=comm + ) + + # Remove the data for these halos from local_halo + for name in local_halo: + local_halo[name] = local_halo[name][~large_halo_mask] # Find size of grid to use to calculate PH keys centres = local_halo["cofp"] @@ -31,7 +68,11 @@ def peano_decomposition(boxsize, local_halo, nr_chunks, comm): grid_size = boxsize / cells_per_dimension nr_cells = cells_per_dimension**3 nr_halos = centres.shape[0] # number of halos on this rank - total_nr_halos = comm.allreduce(nr_halos) # number on all ranks + total_nr_small_halos = comm.allreduce(nr_halos) # number on all ranks + + # Reduce the number of chunks if necessary so that all chunks + # have at least one halo + nr_chunks = min(nr_chunks, total_nr_small_halos) if comm_rank == 0: print(f"Using Peano domain decomposition with bits={bits_per_dimension}") @@ -52,10 +93,50 @@ def peano_decomposition(boxsize, local_halo, nr_chunks, comm): for name in local_halo: local_halo[name] = psort.fetch_elements(local_halo[name], order, comm=comm) + # Handle the halos which should be on separate chunks + if (nr_large_halo > 0) and (comm_rank == comm_size - 1): + # Sort the halos based on the number of bound particles + argsort = np.argsort(large_halo["nr_bound_part"])[::-1] + for name in large_halo: + large_halo[name] = large_halo[name][argsort] + + # Loop through the halos and assign them to chunks + large_halo_chunk_size = [] + i_threshold = 0 + i_halo = 0 + while i_halo < nr_large_halo: + # The separate_chunks object is already sorted by n_bound_threshold + while ( + separate_chunks[i_threshold]["n_bound_threshold"] + > large_halo["nr_bound_part"][i_halo] + ): + i_threshold += 1 + c = min( + separate_chunks[i_threshold]["n_halo_per_chunk"], + nr_large_halo - i_halo, + ) + large_halo_chunk_size.append(c) + i_halo += c + large_halo_chunk_size = np.array(large_halo_chunk_size, dtype=int) + + # Add the large halos back in + for name in local_halo: + local_halo[name] = np.concatenate( + [ + local_halo[name], + large_halo[name], + ], + axis=0, + ) + else: + large_halo_chunk_size = np.zeros(0, dtype=int) + large_halo_chunk_size = comm.bcast(large_halo_chunk_size, root=comm_size - 1) + # Decide how many halos to put in each chunk - chunk_size = np.zeros(nr_chunks, dtype=int) - chunk_size[:] = total_nr_halos // nr_chunks - chunk_size[: total_nr_halos % nr_chunks] += 1 + chunk_size = np.zeros(nr_chunks + large_halo_chunk_size.shape[0], dtype=int) + chunk_size[:nr_chunks] = total_nr_small_halos // nr_chunks + chunk_size[: total_nr_small_halos % nr_chunks] += 1 + chunk_size[nr_chunks:] = large_halo_chunk_size assert np.sum(chunk_size) == total_nr_halos return chunk_size diff --git a/SOAP/core/halo_centres.py b/SOAP/core/halo_centres.py index 08e606f6..a7e49b2c 100644 --- a/SOAP/core/halo_centres.py +++ b/SOAP/core/halo_centres.py @@ -9,7 +9,6 @@ import unyt import virgo.util.match -import virgo.mpi.gather_array as g import virgo.mpi.parallel_sort as psort from SOAP.catalogue_readers import ( @@ -143,23 +142,25 @@ def __init__( for name in local_halo: local_halo[name] = psort.repartition(local_halo[name], ndesired, comm=comm) - # Store total number of halos - self.nr_local_halos = len(local_halo["index"]) - self.nr_halos = comm.allreduce(self.nr_local_halos, op=MPI.SUM) - - if (self.nr_halos == 0) and (comm_rank == 0): + # Exit if we don't have any halos + if (total_nr_halos == 0) and (comm_rank == 0): print("No halos found, aborting run") comm.Abort(1) - # Reduce the number of chunks if necessary so that all chunks have at least one halo - nr_chunks = min(args.chunks, self.nr_halos) - self.nr_chunks = nr_chunks - # Assign halos to chunk tasks: # This sorts the halos by chunk across all MPI ranks and returns the size of each chunk. chunk_size = domain_decomposition.peano_decomposition( - boxsize, local_halo, nr_chunks, comm + boxsize, + local_halo, + args.chunks, + comm, + args.separate_chunks, ) + self.nr_chunks = chunk_size.shape[0] + + # Store total number of halos + self.nr_local_halos = local_halo["index"].shape[0] + self.nr_halos = comm.allreduce(self.nr_local_halos, op=MPI.SUM) # Compute initial radius to read in about each halo local_halo["read_radius"] = local_halo["search_radius"].copy() @@ -235,9 +236,9 @@ def __init__( # Determine local offset to the first halo in each chunk. # This will be different on each MPI rank. - self.local_chunk_size = np.zeros(nr_chunks, dtype=int) - self.local_chunk_offset = np.zeros(nr_chunks, dtype=int) - for chunk_nr in range(nr_chunks): + self.local_chunk_size = np.zeros(self.nr_chunks, dtype=int) + self.local_chunk_offset = np.zeros(self.nr_chunks, dtype=int) + for chunk_nr in range(self.nr_chunks): # Find the range of local halos which are in this chunk (may be none) i1 = self.chunk_offset[chunk_nr] - self.local_halo_offset if i1 < 0: @@ -263,10 +264,12 @@ def __init__( # index in every chunk for which it has >0 halos. We then find the min and max of # each array element over all MPI ranks. chunk_min_rank = ( - np.ones(nr_chunks, dtype=int) * comm_size + np.ones(self.nr_chunks, dtype=int) * comm_size ) # One more than maximum rank - chunk_max_rank = np.ones(nr_chunks, dtype=int) - 1 # One less than minimum rank - for chunk_nr in range(nr_chunks): + chunk_max_rank = ( + np.ones(self.nr_chunks, dtype=int) - 1 + ) # One less than minimum rank + for chunk_nr in range(self.nr_chunks): if self.local_chunk_size[chunk_nr] > 0: chunk_min_rank[chunk_nr] = comm_rank chunk_max_rank[chunk_nr] = comm_rank @@ -278,7 +281,7 @@ def __init__( assert np.all(chunk_max_rank >= 0) # Check that chunk_[min|max]_rank is consistent with local_chunk_size - for chunk_nr in range(nr_chunks): + for chunk_nr in range(self.nr_chunks): assert ( comm_rank >= chunk_min_rank[chunk_nr] and comm_rank <= chunk_max_rank[chunk_nr] diff --git a/SOAP/core/soap_args.py b/SOAP/core/soap_args.py index 4e974e27..a8274831 100644 --- a/SOAP/core/soap_args.py +++ b/SOAP/core/soap_args.py @@ -140,8 +140,8 @@ def get_soap_args(comm): args.max_ranks_reading = all_args["Parameters"]["max_ranks_reading"] args.output_parameters = all_args["Parameters"]["output_parameters"] args.git_hash = all_args["git_hash"] - args.min_read_radius_cmpc = all_args["calculations"].get("min_read_radius_cmpc", 0) - args.calculations = all_args["calculations"] + args.calculations = all_args.get("calculations", {}) + args.min_read_radius_cmpc = args.calculations.get("min_read_radius_cmpc", 0) # Extra-input files which are optionally passed in the parameter file are # processed the same way as the membership files @@ -192,4 +192,20 @@ def get_soap_args(comm): print(f"Could not find FOF radius catalogue: {fof_filename}") comm.Abort(1) + # This really should be done in parameter_file.py + args.separate_chunks = args.calculations.get("separate_chunks", []) + if not isinstance(args.separate_chunks, list): + print("Invalid form for separate_chunks") + comm.Abort(1) + for threshold in args.separate_chunks: + if ("n_bound_threshold" not in threshold) or ( + "n_halo_per_chunk" not in threshold + ): + print("Invalid form for separate_chunks") + comm.Abort(1) + args.separate_chunks = sorted( + args.separate_chunks, + key=lambda x: -x["n_bound_threshold"], + ) + return args diff --git a/SOAP/core/swift_cells.py b/SOAP/core/swift_cells.py index 30061fe8..db42a7ae 100644 --- a/SOAP/core/swift_cells.py +++ b/SOAP/core/swift_cells.py @@ -93,11 +93,15 @@ def __call__(self, data, cache): file_start = self.file_offset file_end = self.file_offset + self.count - dataset.read_direct( - data[self.ptype][self.dataset].full, - np.s_[file_start:file_end, ...], - np.s_[mem_start:mem_end, ...], - ) + try: + dataset.read_direct( + data[self.ptype][self.dataset].full, + np.s_[file_start:file_end, ...], + np.s_[mem_start:mem_end, ...], + ) + except OSError as e: + print(f'Error reading {dataset_name}') + raise e def identify_datasets(filename, nr_files, ptypes, registry): diff --git a/documentation/SOAP.tex b/documentation/SOAP.tex index a625d5db..20e04eec 100644 --- a/documentation/SOAP.tex +++ b/documentation/SOAP.tex @@ -74,13 +74,17 @@ \section{Particle selection for different halo types} of the property calculated for the bound subhalo (e.g. the aperture cut for \verb+ExclusiveSphere/HalfMassRadiusTotal+ is given by the value of \verb+BoundSubhalo/HalfMassRadiusTotal+), meaning a different radial cut is used for each subhalo. -\paragraph{Inclusive sphere quantities (IS)} apply a aperture radii cut, the same as the exclusive sphere -quantities. However for the inclusive sphere we include all particles within the radius, regardless of their +\paragraph{Inclusive sphere quantities (IS)} apply a aperture radii cut, the same as the exclusive sphere +quantities. However for the inclusive sphere we include all particles within the radius, regardless of their membership status. If the aperture radius of an inclusive sphere variation is greater than the EncloseRadius (the maximum distance between a bound particle and the halo centre) of a subhalo, then for that subhalo no properties are computed for the variation. The quantities are stored within the group \verb+InclusiveSphere+. +For each subhalo InclusiveSphere variations are calculated for all apertures up to and including the first +aperture that exceeds the EncloseRadius(the maximum distance between any bound particle and the halo centre). +Apertures larger than this are not computed. + \paragraph{Exclusive projected quantities (EP)} are similar to exclusive sphere quantities, except that their aperture cut is applied in projection. For each radii there are three independent projections: along the x-, y- and z-axis. Along the projection axis, we do not apply any radial cut, meaning the depth corresponds to all particles diff --git a/parameter_files/COLIBRE_HYBRID.yml b/parameter_files/COLIBRE_HYBRID.yml index 4f823849..30670830 100644 --- a/parameter_files/COLIBRE_HYBRID.yml +++ b/parameter_files/COLIBRE_HYBRID.yml @@ -794,3 +794,10 @@ calculations: cold_dense_gas_filter: maximum_temperature_K: 3.16e4 minimum_hydrogen_number_density_cm3: 0.1 + separate_chunks: + - n_bound_threshold: 300000000 + n_halo_per_chunk: 1 + - n_bound_threshold: 100000000 + n_halo_per_chunk: 2 + - n_bound_threshold: 50000000 + n_halo_per_chunk: 5 diff --git a/parameter_files/COLIBRE_THERMAL.yml b/parameter_files/COLIBRE_THERMAL.yml index 7e2eb418..7d5ae39b 100644 --- a/parameter_files/COLIBRE_THERMAL.yml +++ b/parameter_files/COLIBRE_THERMAL.yml @@ -794,3 +794,10 @@ calculations: cold_dense_gas_filter: maximum_temperature_K: 3.16e4 minimum_hydrogen_number_density_cm3: 0.1 + separate_chunks: + - n_bound_threshold: 300000000 + n_halo_per_chunk: 1 + - n_bound_threshold: 100000000 + n_halo_per_chunk: 2 + - n_bound_threshold: 50000000 + n_halo_per_chunk: 5 diff --git a/parameter_files/README.md b/parameter_files/README.md index 6e058b95..e27a72ea 100644 --- a/parameter_files/README.md +++ b/parameter_files/README.md @@ -261,3 +261,12 @@ Contains information about how to run SOAP - **maximum_temperature_K**: Value above which gas is not considered to be cold - **minimum_hydrogen_number_density_cm3**: Value below which gas gas is not considered to be dense - **strict_halo_copy**: Optional, default False. When a halo has multiple ExclusiveSphere/ProjectedAperture halo types which encompass all the bound particles then we just copy across the values rather than recomputing them. There are a small number of properties for which this is not correct. If this flag is set then these properties are set to zero for the larger apertures instead of being copied across. +- **separate_chunks**: Optional, default []. SOAP processes subhalos in parallel, but this can cause memory issues if there are subhalos which take up a significant fraction of a node's memory. This parameter allows a list of dictionaries to be passed. Each dictionary must contain two keys: `n_bound_threshold` (which specifies the number of bound particles above which a subhalo should be treated differently) and `n_halo_per_chunk` (which gives the maximum number of subhalos of this size which can be placed on a single chunk). An example is +``` + separate_chunks: + - n_bound_threshold: 1000 + n_halo_per_chunk: 10 + - n_bound_threshold: 10000 + n_halo_per_chunk: 1 +``` +In this case any subhalo with more than 10000 bound particles would be placed on its own chunk, subhalos with more than 1000 (but less than 10000) bound particles would be grouped into sets of 10 and there will be a chunk for each group of 10, and all other subhalos would be chunked as normal. Note that grouping subhalos based on number of bound particles is not efficient since in general they will not share particles, so this parameter should only be used if required. The threshold values will depend on the system memory, and also on the number of properties being computed by SOAP. diff --git a/parameter_files/XRAY_ONLY.yml b/parameter_files/XRAY_ONLY.yml deleted file mode 100644 index f1dbd640..00000000 --- a/parameter_files/XRAY_ONLY.yml +++ /dev/null @@ -1,147 +0,0 @@ -# Values in this section are substituted into the other sections -Parameters: - sim_dir: /cosma8/data/dp004/colibre/Runs - output_dir: /snap8/scratch/dp004/dc-mcgi1/soap_xray_props - scratch_dir: /snap8/scratch/dp004/dc-mcgi1/soap_xray_props - -# Location of the Swift snapshots: -Snapshots: - # Use {snap_nr:04d} for the snapshot number and {file_nr} for the file number. - filename: "{sim_dir}/{sim_name}/snapshots/colibre_{snap_nr:04d}/colibre_{snap_nr:04d}.{file_nr}.hdf5" - -# Which halo finder we're using, and base name for halo finder output files -HaloFinder: - type: HBTplus - filename: "{sim_dir}/{sim_name}/HBTplus/{snap_nr:03d}/SubSnap_{snap_nr:03d}" - # fof_filename: "{sim_dir}/{sim_name}/fof/fof_output_{snap_nr:04d}/fof_output_{snap_nr:04d}.{file_nr}.hdf5" - #type: VR - #filename: "{sim_dir}/halo_{snap_nr:04d}" - #type: Subfind - #filename: "{sim_dir}/snapdir_{snap_nr:03d}/snapshot_{snap_nr:03d}" - -GroupMembership: - # Where to write the group membership files - filename: "{sim_dir}/{sim_name}/SOAP/membership_{snap_nr:04d}/membership_{snap_nr:04d}.{file_nr}.hdf5" - -HaloProperties: - # Where to write the halo properties file - filename: "{output_dir}/{sim_name}/SOAP_uncompressed/halo_properties_{snap_nr:04d}.hdf5" - # Where to write temporary chunk output - chunk_dir: "{scratch_dir}/{sim_name}/SOAP-tmp/" - -ApertureProperties: - properties: - StellarMass: true - variations: - exclusive_50_kpc: - inclusive: false - radius_in_kpc: 50.0 -ProjectedApertureProperties: - properties: - {} - variations: - {} -SOProperties: - properties: - XRayLuminosityDensityCut: true - XRayLuminosity: - snapshot: true - snipshot: false - XRayLuminosityWithoutRecentAGNHeating: - snapshot: true - snipshot: false - XRayLuminosityCoreExcision: - snapshot: true - snipshot: false - XRayLuminosityWithoutRecentAGNHeatingCoreExcision: - snapshot: true - snipshot: false - XRayLuminosityInRestframe: false - XRayLuminosityInRestframeWithoutRecentAGNHeating: false - XRayLuminosityInRestframeCoreExcision: false - XRayLuminosityInRestframeWithoutRecentAGNHeatingCoreExcision: false - XRayPhotonLuminosity: - snapshot: true - snipshot: false - XRayPhotonLuminosityWithoutRecentAGNHeating: - snapshot: true - snipshot: false - XRayPhotonLuminosityCoreExcision: - snapshot: true - snipshot: false - XRayPhotonLuminosityWithoutRecentAGNHeatingCoreExcision: - snapshot: true - snipshot: false - XRayPhotonLuminosityInRestframe: false - XRayPhotonLuminosityInRestframeWithoutRecentAGNHeating: false - XRayPhotonLuminosityInRestframeCoreExcision: false - XRayPhotonLuminosityInRestframeWithoutRecentAGNHeatingCoreExcision: false - TotalMass: true - XRayLuminosityNoSat: true - XRayLuminosityNoSatCoreExcision: true - variations: - 200_crit: - type: crit - value: 200.0 - 500_crit: - type: crit - value: 500.0 - core_excision_fraction: 0.15 -SubhaloProperties: - properties: - EncloseRadius: true - NumberOfBlackHoleParticles: true - NumberOfDarkMatterParticles: true - NumberOfGasParticles: true - NumberOfStarParticles: true - TotalMass: true -aliases: - PartType0/LastSNIIKineticFeedbackDensities: PartType0/DensitiesAtLastSupernovaEvent - PartType0/LastSNIIThermalFeedbackDensities: PartType0/DensitiesAtLastSupernovaEvent - snipshot: - PartType0/SpeciesFractions: PartType0/ReducedSpeciesFractions - PartType0/ElementMassFractions: PartType0/ReducedElementMassFractions - PartType0/LastSNIIKineticFeedbackDensities: PartType0/DensitiesAtLastSupernovaEvent - PartType0/LastSNIIThermalFeedbackDensities: PartType0/DensitiesAtLastSupernovaEvent -filters: - general: - limit: 100 - properties: - - BoundSubhalo/NumberOfGasParticles - - BoundSubhalo/NumberOfDarkMatterParticles - - BoundSubhalo/NumberOfStarParticles - - BoundSubhalo/NumberOfBlackHoleParticles - combine_properties: sum - baryon: - limit: 0 - properties: - - BoundSubhalo/NumberOfGasParticles - - BoundSubhalo/NumberOfStarParticles - combine_properties: sum - dm: - limit: 0 - properties: - - BoundSubhalo/NumberOfDarkMatterParticles - gas: - limit: 0 - properties: - - BoundSubhalo/NumberOfGasParticles - star: - limit: 0 - properties: - - BoundSubhalo/NumberOfStarParticles -defined_constants: - O_H_sun: 4.9e-4 - Fe_H_sun: 3.16e-5 - N_O_sun: 0.138 - C_O_sun: 0.549 - Mg_H_sun: 3.98e-5 -calculations: - calculate_missing_properties: false - stict_halo_copy: false - recently_heated_gas_filter: - delta_time_myr: 15 - use_AGN_delta_T: false - cold_dense_gas_filter: - maximum_temperature_K: 3.16e4 - minimum_hydrogen_number_density_cm3: 0.1 From 6d9bc69c1321ab1be89cf21a68e23f92420678df Mon Sep 17 00:00:00 2001 From: robjmcgibbon <36136863+robjmcgibbon@users.noreply.github.com> Date: Fri, 1 May 2026 10:41:22 +0100 Subject: [PATCH 15/19] Property table in rst format (#205) * Add rst generation * Move dropdown * Updates * Joop comments * Add links for citations * Output human units * Remove old footnote --- .gitignore | 1 + README.md | 2 + SOAP/property_table.py | 522 +++++++++++++++++- documentation/footnote_AngMom.tex | 6 +- documentation/footnote_MBH.tex | 3 +- documentation/footnote_Mnu.tex | 2 +- documentation/footnote_Tgas.tex | 7 +- documentation/footnote_Xray.tex | 4 +- documentation/footnote_circvel.tex | 5 +- documentation/footnote_com.tex | 2 +- documentation/footnote_compY.tex | 8 +- documentation/footnote_concentration.tex | 11 +- documentation/footnote_coreexcision.tex | 6 +- documentation/footnote_disc_fraction.tex | 3 +- documentation/footnote_dopplerB.tex | 6 +- documentation/footnote_elements.tex | 3 - documentation/footnote_halfmass.tex | 4 +- documentation/footnote_kappa.tex | 4 +- documentation/footnote_lum.tex | 16 +- .../footnote_progenitor_descendant.tex | 2 +- documentation/footnote_proj_veldisp.tex | 4 +- documentation/footnote_satfrac.tex | 6 +- documentation/footnote_spin.tex | 2 +- documentation/footnote_tensor.tex | 4 +- documentation/footnote_veldisp_matrix.tex | 4 +- documentation/property_table_intro.rst | 20 + 26 files changed, 584 insertions(+), 73 deletions(-) delete mode 100644 documentation/footnote_elements.tex create mode 100644 documentation/property_table_intro.rst diff --git a/.gitignore b/.gitignore index 3547d588..914afcd6 100644 --- a/.gitignore +++ b/.gitignore @@ -15,6 +15,7 @@ documentation/table.tex documentation/timestamp.tex documentation/variations_table.tex documentation/units.tex +documentation/property_table.rst tests/FLAMINGO/test_parameters.yml tests/COLIBRE/test_parameters.yml diff --git a/README.md b/README.md index 73148e35..48e1fc96 100644 --- a/README.md +++ b/README.md @@ -125,6 +125,8 @@ lossless compression to SOAP catalogues. A pdf describing the SOAP output can be generated. First run `SOAP/property_table.py` passing the parameter file used to run SOAP (to get the properties and halo types to include) and a snapshot (to get the units), e.g. `python SOAP/property_table.py parameter_files/COLIBRE_THERMAL.yml /cosma8/data/dp004/colibre/Runs/L0100N0752/Thermal/snapshots/colibre_0127/colibre_0127.hdf5`. This will generate a table containing all the properties which are enabled in the parameter file. To create the pdf run `cd documentation; pdflatex SOAP.tex; pdflatex SOAP.tex`. If you wish to see all possible properties then first run `python SOAP/property_table.py`, and then generate the pdf. +The `property_table.py` script also generates the file `documentation/property_table.rst`, which is a version of the property table that can be useful for sphnix websites. + ### Slurm scripts for running on COSMA The files in the `scripts` directory are made for running on cosma. diff --git a/SOAP/property_table.py b/SOAP/property_table.py index d86fdfe0..b2bb0446 100644 --- a/SOAP/property_table.py +++ b/SOAP/property_table.py @@ -143,14 +143,6 @@ class PropertyTable: "proj_veldisp_dm", "proj_veldisp_star", ], - "footnote_elements.tex": [ - "gasOfrac", - "gasOfrac_SF", - "gasFefrac", - "gasFefrac_SF", - "gasmetalfrac", - "gasmetalfrac_SF", - ], "footnote_halfmass.tex": [ "HalfMassRadiusTot", "HalfMassRadiusGas", @@ -1860,7 +1852,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="dimensionless", - description="Fraction of mass that is bound to a satellite in the same FOF group.", + description="Fraction of mass that is bound to a satellite in the same FoF group.", lossy_compression_filter="FMantissa9", dmo_property=True, particle_properties=[ @@ -1885,7 +1877,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="dimensionless", - description="Fraction of mass that is bound to a satellite outside this FOF group.", + description="Fraction of mass that is bound to a satellite outside this FoF group.", lossy_compression_filter="FMantissa9", dmo_property=True, particle_properties=[ @@ -4645,7 +4637,7 @@ class PropertyTable: shape=1, dtype=np.int64, unit="dimensionless", - description="ID of the host FOF halo of this subhalo. Hostless halos have HostFOFId == -1", + description="ID of the host FoF halo of this subhalo. Hostless halos have HostFOFId == -1", lossy_compression_filter="None", dmo_property=True, particle_properties=[], @@ -4766,7 +4758,7 @@ class PropertyTable: shape=3, dtype=np.float64, unit="snap_length", - description="Centre of mass of the host FOF halo of this subhalo. Zero for satellite and hostless subhalos.", + description="Centre of mass of the host FoF halo of this subhalo. Zero for satellite and hostless subhalos.", lossy_compression_filter="DScale6", dmo_property=True, particle_properties=[], @@ -4778,7 +4770,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="snap_mass", - description="Mass of the host FOF halo of this subhalo. Zero for satellite and hostless subhalos.", + description="Mass of the host FoF halo of this subhalo. Zero for satellite and hostless subhalos.", lossy_compression_filter="FMantissa9", dmo_property=True, particle_properties=[], @@ -4790,7 +4782,7 @@ class PropertyTable: shape=1, dtype=np.uint64, unit="dimensionless", - description="Number of particles in the host FOF halo of this subhalo. Zero for satellite and hostless subhalos.", + description="Number of particles in the host FoF halo of this subhalo. Zero for satellite and hostless subhalos.", lossy_compression_filter="None", dmo_property=True, particle_properties=[], @@ -4802,7 +4794,7 @@ class PropertyTable: shape=1, dtype=np.float32, unit="snap_length", - description="Radius of the particle furthest from the FOF centre of mass. Zero for satellite and hostless subhalos. Missing for older runs.", + description="Radius of the particle furthest from the FoF centre of mass. Zero for satellite and hostless subhalos. Missing for older runs.", lossy_compression_filter="FMantissa9", dmo_property=True, particle_properties=[], @@ -4948,7 +4940,9 @@ def add_properties(self, halo_property: HaloProperty, halo_type: str): units = units * unyt.Unit("a") ** prop.a_scale_exponent prop_unit = units.units.latex_repr.replace( "\\rm{km} \\cdot \\rm{kpc}", "\\rm{kpc} \\cdot \\rm{km}" - ).replace("\\frac{\\rm{km}^{2}}{\\rm{s}^{2}}", "\\rm{km}^{2} / \\rm{s}^{2}") + ).replace( + "\\frac{\\rm{km}^{2}}{\\rm{s}^{2}}", "\\rm{km}^{2} / \\rm{s}^{2}" + ).rstrip() prop_dtype = prop.dtype.__name__ @@ -5290,6 +5284,500 @@ def generate_tex_files(self, output_dir: str): vel_kms = (1 * unyt.snap_length / unyt.snap_time).to("km/s").value ofile.write(f"\\newcommand{{\\velbaseunit}}{{{vel_kms:.4g}}}\n") + ######### RST generation + + @staticmethod + def _camel_to_snake(name: str) -> str: + """ + Convert a CamelCase property name to snake_case, matching the + convention used by swiftsimio. + """ + + return re.sub("([a-z0-9])([A-Z])", r"\1_\2", name).lower() + + @staticmethod + def _latex_units_to_rst(latex_units: str) -> str: + """ + Convert a LaTeX units string (as produced by unyt) into RST math format + suitable for use in Sphinx documentation. + + Wraps the unit expression in :math:`...` and keeps the LaTeX content + as-is, since unyt already produces valid LaTeX math. + """ + if not latex_units or latex_units == "dimensionless": + return "dimensionless" + return f":math:`{latex_units}`" + + @staticmethod + def _compression_to_rst(compression_description: str) -> str: + """ + Convert a compression description string (which may contain LaTeX math) + into RST format. Dollar-sign-delimited math is replaced with :math:`...`. + """ + + def replace_math(match): + return f":math:`{match.group(1)}`" + + return re.sub(r"\$(.+?)\$", replace_math, compression_description) + + @staticmethod + def _latex_to_rst_math(text: str) -> str: + """ + Convert inline LaTeX math delimited by ``$...$`` in *text* to RST + ``:math:`...``` roles. Dollar signs that are not part of a math + expression are left unchanged. + """ + + def replace_math(match): + return f":math:`{match.group(1)}`" + + return re.sub(r"\$(.+?)\$", replace_math, text) + + # Map from internal snap_* unit names to human-readable physical equivalents. + # Values are validated against the live unyt registry by _validate_snap_unit_map. + _snap_unit_map = { + "snap_mass": "(1e10*Msun)", + "snap_length": "Mpc", + "snap_time": "(s*Mpc/km)", + "snap_temperature": "K", + } + + def _resolve_snap_units(self, prop_name: str, prop: dict) -> str: + """ + Return an RST-formatted unit string for *prop_name*, replacing the + internal ``snap_*`` base units with their physical equivalents. + + The substitution map is validated once against the live unyt registry + in :meth:`_validate_snap_unit_map` (called from ``generate_rst_files``). + """ + # Retrieve the raw unit expression from the canonical property definition + prop_obj = self.full_property_list[prop_name] + prop_unit = prop_obj.unit + + for snap_name, replacement in self._snap_unit_map.items(): + prop_unit = prop_unit.replace(snap_name, replacement) + + units = unyt.unyt_quantity(1, units=prop_unit) + latex = ( + units.units.latex_repr + .replace("\\rm{km} \\cdot \\rm{kpc}", "\\rm{kpc} \\cdot \\rm{km}") + .replace("\\frac{\\rm{km}^{2}}{\\rm{s}^{2}}", "\\rm{km}^{2} / \\rm{s}^{2}") + .replace(r"1.0 \times ", "") + .rstrip() + ) + + if (not prop_obj.output_physical) and prop_obj.a_scale_exponent: + a = prop_obj.a_scale_exponent + a_str = "a" if a == 1 else f"a^{{{a}}}" + latex = f"{a_str} \\cdot {latex}" if latex else a_str + + return self._latex_units_to_rst(latex) + + def _validate_snap_unit_map(self) -> None: + """ + Assert that every entry in ``_snap_unit_map`` matches the live unyt + registry. Called once at the start of ``generate_rst_files`` so the + assertions run exactly once rather than once per property. + """ + for snap_name, replacement in self._snap_unit_map.items(): + assert np.isclose( + unyt.unyt_quantity(1, snap_name), + unyt.unyt_quantity(1, replacement), + ), f"Unit mismatch: {snap_name} != {replacement}" + + def _get_output_types_rst(self, prop: dict) -> str: + """ + Return a highlighted string showing all halo type abbreviations. + + All five types are always shown, colour-coded by availability: + - ``:avail:`` (green) -- computed for all snapshots + - ``:snaponly:`` (blue) -- computed for snapshots only, not snipshots + - ``:unavail:`` (red) -- not computed + + All three roles must be registered in ``conf.py`` + + The abbreviations correspond to: + BS - BoundSubhalo (SubhaloProperties) + ES - ExclusiveSphere (ExclusiveSphereProperties) + IS - InclusiveSphere (InclusiveSphereProperties) + P - ProjectedAperture (ProjectedApertureProperties) + SO - SphericalOverdensity (SOProperties) + """ + type_map = [ + ("SubhaloProperties", "BS"), + ("ExclusiveSphereProperties", "ES"), + ("InclusiveSphereProperties", "IS"), + ("ProjectedApertureProperties", "EP"), + ("SOProperties", "SO"), + ] + tokens = [] + for halo_type, abbrev in type_map: + # A type is either absent, present for all outputs, or snapshot-only. + # The types list contains the exact string "SubhaloProperties" for + # full availability, or "SnapshotOnlySubhaloProperties" for + # snapshot-only — never both for the same halo type. + full_present = halo_type in prop["types"] + snap_only = f"SnapshotOnly{halo_type}" in prop["types"] + if full_present: + role = "avail" + elif snap_only: + role = "snaponly" + else: + role = "unavail" + tokens.append(f":{role}:`{abbrev}`") + return " ".join(tokens) + + def get_footnotes_rst(self, name: str) -> tuple[list[int], str]: + """ + Return RST footnote references for property *name*. + + Returns a tuple of: + - list of footnote numbers (ints) associated with this property + - string of RST hyperlinks, e.g. `` `¹ `_ ``, or ``""`` + + We use unicode superscript digits as link text, pointed at named RST + targets (``.. _footnote-N:``) in the Footnotes section. This is a + plain inline hyperlink with no role nesting, so it works correctly + inside ``list-table`` cells. + """ + footnote_nums = [] + for fnote in self.explanation.keys(): + if name in self.explanation[fnote]: + try: + i = self.footnotes.index(fnote) + except ValueError: + i = len(self.footnotes) + self.footnotes.append(fnote) + footnote_nums.append(i + 1) + if footnote_nums: + refs = " ".join( + f"`[{n}] `_" + for n in sorted(footnote_nums) + ) + return footnote_nums, refs + return [], "" + + def _read_footnote_rst(self, fnote_filename: str, footnote_number: int) -> str: + r""" + Read a LaTeX footnote file from the ``documentation/`` directory and + return an RST version suitable for embedding in the output RST file. + """ + import re + + filepath = f"documentation/{fnote_filename}" + try: + with open(filepath, "r") as fh: + text = fh.read() + except FileNotFoundError: + return f"*(Footnote file ``{fnote_filename}`` not found.)*" + + # Step 1: substitute the plain-text footnote number placeholder FIRST, + # before any dollar-sign processing, so that the literal string + # "$FOOTNOTE_NUMBER$" is replaced with e.g. "5" and can no longer be + # mistaken for a LaTeX math delimiter. + text = text.replace("$FOOTNOTE_NUMBER$", str(footnote_number)) + + # Step 2: substitute parameter-file-derived values (also plain text). + if "$LOG_COLD_GAS_TEMP$" in text: + params = self.parameters.get_cold_dense_params() + if params.get("initialised"): + T = f'{np.log10(params["maximum_temperature_K"]):.2g}' + text = text.replace("$LOG_COLD_GAS_TEMP$", T) + if "$LOG_COLD_GAS_DENSITY$" in text: + params = self.parameters.get_cold_dense_params() + if params.get("initialised"): + rho = f'{np.log10(params["minimum_hydrogen_number_density_cm3"]):.2g}' + text = text.replace("$LOG_COLD_GAS_DENSITY$", rho) + + # Step 3: convert \begin{equation}...\end{equation} blocks to + # RST ``.. math::`` directives. + # + # We stash each converted block behind a NUL-delimited placeholder so + # that subsequent LaTeX-stripping steps cannot corrupt the math content. + math_blocks: list[str] = [] + + def equation_to_rst(match): + math_content = match.group(1).strip() + # Indent each line of the math block by 3 spaces (directive body) + indented = "\n".join(f" {line}" for line in math_content.splitlines()) + block = f"\n\n.. math::\n\n{indented}\n\n" + placeholder = f"\x00MATHBLOCK{len(math_blocks)}\x00" + math_blocks.append(block) + return placeholder + + text = re.sub( + r"\\begin\{equation\}(.*?)\\end\{equation\}", + equation_to_rst, + text, + flags=re.DOTALL, + ) + + # Step 3b: convert \begin{enumerate}...\end{enumerate} blocks to + # RST auto-numbered lists. Each \item becomes a ``#.`` list entry. + # The converted block is stashed behind a placeholder so later + # stripping steps leave its content alone. + def enumerate_to_rst(match): + body = match.group(1) + # Split on \item, discard the empty string before the first \item + items = re.split(r"\\item\s*", body) + rst_items = [] + for item in items: + item = item.strip() + if item: + # Collapse internal newlines to spaces within each item + item = " ".join(item.split()) + # Convert inline math and non-breaking spaces within the + # item text now, before the block is stashed as a + # placeholder (the main conversion steps run before restore) + item = re.sub(r"\$([^\n\$]+?)\$", lambda m: f":math:`{m.group(1)}`", item) + item = item.replace("~", " ") + rst_items.append(f"#. {item}") + block = "\n\n" + "\n".join(rst_items) + "\n\n" + placeholder = f"\x00MATHBLOCK{len(math_blocks)}\x00" + math_blocks.append(block) + return placeholder + + text = re.sub( + r"\\begin\{enumerate\}(.*?)\\end\{enumerate\}", + enumerate_to_rst, + text, + flags=re.DOTALL, + ) + + # Step 4: process \paragraph{$^{N}$Title}\label{footnote:N} + # + # The paragraph command has the form: + # \paragraph{$^{N}$Title text}\label{footnote:N} rest of sentence... + # + # We want to: + # - Drop the superscript prefix "$^{N}$" (the number is already + # shown by the caller as **[N]**). + # - Render the title in bold. + # - Keep the rest of the sentence on the same line. + # - Strip the \label{...} command entirely. + # + # We use a lazy ``.*?`` match for the paragraph argument because the + # content may contain braces (e.g. the ``$^{5}$`` superscript prefix), + # which would trip up a ``[^}]*`` character class. + def paragraph_to_rst(match): + # match.group(1): everything inside \paragraph{...} + # match.group(2): the \label{...} target (discarded) + # match.group(3): the remainder of the line after \label{...} + inner = match.group(1) + rest = match.group(3).strip() + # Strip leading superscript: $^{N}$ or $^N$ + inner = re.sub(r"^\$\^{?[^}$]*}?\$\s*", "", inner, flags=re.DOTALL) + # Collapse any internal newlines in the title to a single space + inner = " ".join(inner.split()) + return f"**{inner}** {rest}" + + # re.DOTALL lets .*? in group 1 match newlines (title may wrap across + # lines). The final group uses [^\n]* to stay single-line so it only + # captures the rest of the \label{...} line, not the whole document. + text = re.sub( + r"\\paragraph\{(.*?)\}\\label\{[^}]*\}([ \t]*)([^\n]*)", + paragraph_to_rst, + text, + flags=re.DOTALL, + ) + + # Step 5: convert $...$ inline math → :math:`...` + # Use a non-greedy match that does not cross newlines, to avoid + # accidentally swallowing large sections of text. + text = re.sub(r"\$([^\n$]+?)\$", lambda m: f":math:`{m.group(1)}`", text) + + # Step 6: common LaTeX text commands → RST equivalents. + text = re.sub(r"\\verb\+([^+]+)\+", lambda m: f"``{m.group(1)}``", text) + text = re.sub( + r"\\(?:textit|emph)\{([^}]+)\}", lambda m: f"*{m.group(1)}*", text + ) + text = re.sub(r"\\textbf\{([^}]+)\}", lambda m: f"**{m.group(1)}**", text) + + # Step 7: strip remaining unrecognised LaTeX commands. + # \href{url}{text} is handled explicitly before the generic stripper, + # which would otherwise keep only the first argument (the URL) and + # discard the link text + # Commands with a brace argument: keep the argument text. + # Bare commands (no argument): remove entirely. + text = re.sub( + r"\\href\{([^}]+)\}\{([^}]+)\}", + lambda m: f"`{m.group(2)} <{m.group(1)}>`_", + text, + ) + text = re.sub(r"\\[a-zA-Z]+\{([^}]*)\}", r"\1", text) + text = re.sub(r"\\[a-zA-Z]+", "", text) + text = text.replace("~", " ") # LaTeX non-breaking space + + # Step 8: restore the protected .. math:: blocks now that all LaTeX + # stripping is complete. + for i, block in enumerate(math_blocks): + text = text.replace(f"\x00MATHBLOCK{i}\x00", block) + + # Step 9: tidy up whitespace. + text = re.sub(r"\n{3,}", "\n\n", text).strip() + + return text + + def _build_rst_table(self, prop_names: list, lines: list): + """ + Append RST list-table rows for the given properties to *lines*. + + Note: ``self.footnotes`` is populated as a side-effect (via + ``get_footnotes_rst``), so this must be called before the footnote + section is rendered. + """ + lines.append(".. list-table::") + lines.append(" :widths: 25 10 15 50") + lines.append(" :header-rows: 1") + lines.append("") + lines.append(" * - Name") + lines.append(" - Filter") + lines.append(" - Variations") + lines.append(" - Description") + + for prop_name in prop_names: + prop = self.properties[prop_name] + + # swiftsimio name + snake_name = ".".join( + self._camel_to_snake(part) for part in prop['name'].split("/") + ) + + # HDF5 output name (prepend InputHalos/ where applicable) + output_name = prop["name"] + if output_name.split("/")[0] in ("HBTplus", "VR", "FOF"): + output_name = "InputHalos/" + output_name + snake_name = 'input_halos_' + snake_name + + prop_units_rst = self._resolve_snap_units(prop_name, prop) + prop_comp_rst = self._compression_to_rst( + self.compression_description[prop["compression"]] + ) + prop_filter = prop['category'] + output_types = self._get_output_types_rst(prop) + if 'DummyProperties' in prop['types']: + prop_filter = 'basic' + output_types = r'\-' + description = prop["description"].format( + label="satisfying a spherical overdensity criterion.", + core_excision="excised core", + ) + + _, footnote_rst = self.get_footnotes_rst(prop_name) + display_name = f"``{snake_name}``" + if footnote_rst: + description = f"{description} {footnote_rst}" + + lines.append(f" * - .. dropdown:: {display_name}") + lines.append("") + lines.append(f" * **HDF5 name:** ``{output_name}``") + lines.append(f" * **Shape:** {prop['shape']}") + lines.append(f" * **Type:** {prop['dtype']}") + lines.append(f" * **Units:** {prop_units_rst}") + lines.append(f" * **Compression:** {prop_comp_rst}") + lines.append(f" - {prop_filter}") + lines.append(f" - {output_types}") + lines.append(f" - {description}") + + lines.append("") + + def generate_rst_files(self, output_dir: str): + """ + Generate a single RST documentation file (``property_table.rst``) for the + SPHYNIX documentation system. + + A preface to the tables is loaded from (``documentation/property_table_intro.rst``). + + The file contains four ``list-table`` directives separated by titles. + + RST footnotes referenced in the tables are appended at the bottom of + the file so that the superscript links in the Name column resolve + correctly. + """ + + self._validate_snap_unit_map() + + category_order = [ + "basic", + "general", + "gas", + "dm", + "star", + "baryon", + "InputHalos", + "VR", + "HBTplus", + "FOF", + "SOAP", + ] + prop_names_sorted = sorted( + self.properties.keys(), + key=lambda key: ( + category_order.index(self.properties[key]["category"]) + if self.properties[key]["category"] in category_order + else len(category_order), + self.properties[key]["name"].lower(), + ), + ) + + # Sort properties into the different tables + input_basic_props = [] + input_copied_props = [] + dmo_props = [] + hydro_props = [] + for n in prop_names_sorted: + if self.properties[n]['category'] == 'InputHalos': + input_basic_props.append(n) + elif 'DummyProperties' in self.properties[n]['types']: + input_copied_props.append(n) + elif self.properties[n]['dmo']: + dmo_props.append(n) + else: + hydro_props.append(n) + + # Reset footnote list; it is populated as a side-effect of + # _build_rst_table calling get_footnotes_rst for each property. + self.footnotes = [] + + # Introduction + with open(f'documentation/property_table_intro.rst') as file: + lines = [line.rstrip() for line in file.readlines()] + + # Tables + for props, table_title in [ + (input_basic_props, "Input halo properties"), + (dmo_props, "Dark matter only properties"), + (hydro_props, "Hydrodynamical properties"), + (input_copied_props, "Copied properties"), + ]: + + lines.append(table_title) + lines.append("-" * len(table_title)) + lines.append("") + self._build_rst_table(props, lines) + + # Footnotes + if self.footnotes: + fn_title = "Footnotes" + lines.append(fn_title) + lines.append("-" * len(fn_title)) + lines.append("") + # Each footnote is preceded by a named RST target ``.. _footnote-N:`` + # which is what the ``footnote-N_`` internal references in the table + # Name cells resolve to. + for i, fnote_filename in enumerate(self.footnotes): + fn_number = i + 1 + fn_text = self._read_footnote_rst(fnote_filename, fn_number) + lines.append(f".. _footnote-{fn_number}:") + lines.append("") + lines.append(f"**[{fn_number}]** {fn_text}") + lines.append("") + + os.makedirs(output_dir, exist_ok=True) + with open(f"{output_dir}/property_table.rst", "w") as ofile: + ofile.write("\n".join(lines)) + + class DummyProperties: """ @@ -5388,6 +5876,7 @@ def get_parameter_file_all_properties(): You must pass a parameter file and a snapshot to run this script """ + import re import sys import h5py import yaml @@ -5467,3 +5956,4 @@ def get_parameter_file_all_properties(): ) table.generate_tex_files("documentation") + table.generate_rst_files("documentation") diff --git a/documentation/footnote_AngMom.tex b/documentation/footnote_AngMom.tex index a0e6bfaf..a1e4a481 100644 --- a/documentation/footnote_AngMom.tex +++ b/documentation/footnote_AngMom.tex @@ -1,5 +1,5 @@ -\paragraph{$^{$FOOTNOTE_NUMBER$}$The angular momentum}\label{footnote:$FOOTNOTE_NUMBER$} of gas, dark matter and stars is computed relative to -the halo centre (cop) and the centre of mass velocity of that particular component, and not to the +\paragraph{$^{$FOOTNOTE_NUMBER$}$The angular momentum}\label{footnote:$FOOTNOTE_NUMBER$} of gas, dark matter, or stars is computed relative to +the halo centre and the centre of mass velocity of that particular component, and not to the total centre of mass velocity. The full expression is \begin{equation} @@ -22,5 +22,5 @@ \vec{v}_{\rm{}com,comp} = \frac{\sum_{i={\rm{}comp}} m_i \vec{v}_i}{\sum_{i={\rm{}comp}} m_i}. \end{equation} -For FLAMINGO, we also compute the angular momentum for baryons, where the sum is then over both gas and star +We also compute the angular momentum for baryons, where the sum is then over both gas and star particles. diff --git a/documentation/footnote_MBH.tex b/documentation/footnote_MBH.tex index 2aa71f9d..3132aa35 100644 --- a/documentation/footnote_MBH.tex +++ b/documentation/footnote_MBH.tex @@ -1,2 +1 @@ -\paragraph{$^{$FOOTNOTE_NUMBER$}$The most massive black hole}\label{footnote:$FOOTNOTE_NUMBER$} is identified based on the BH subgrid mass (i.e. -the same mass that goes into \verb+BlackHolesSubgridMass+). +\paragraph{$^{$FOOTNOTE_NUMBER$}$The most massive black hole}\label{footnote:$FOOTNOTE_NUMBER$} is identified based using subgrid masses of the black holes. diff --git a/documentation/footnote_Mnu.tex b/documentation/footnote_Mnu.tex index a94eef08..19efb023 100644 --- a/documentation/footnote_Mnu.tex +++ b/documentation/footnote_Mnu.tex @@ -6,5 +6,5 @@ M_{\nu{},{\rm{}NS}} = \sum_i m_i w_i + \frac{4\pi{}}{3} \rho{}_{\nu{}} R_{\rm{}SO}^3, \end{equation} -where $w_i$ are the neutrino weights (which can be negative), and $\rho{}_{\nu{}}$ is the background density +where $w_i$ are the neutrino weights (which can be negative), and $\rho_{\nu}$ is the background density of neutrinos that is also used in the SO radius calculation. The latter is obtained from the snapshot header. diff --git a/documentation/footnote_Tgas.tex b/documentation/footnote_Tgas.tex index fa74ed0c..0dad274f 100644 --- a/documentation/footnote_Tgas.tex +++ b/documentation/footnote_Tgas.tex @@ -8,13 +8,14 @@ that satisfy \begin{equation} - \verb+LastAGNFeedbackScaleFactors+_i \geq{} a - 15{\rm{}Myr} + 0.1 \Delta{}T_{\rm{}AGN} \leq{} T_i \leq{} 10^{0.3} \Delta{}T_{\rm{}AGN}, \end{equation} and \begin{equation} - 0.1 \Delta{}T_{\rm{}AGN} \leq{} T_i \leq{} 10^{0.3} \Delta{}T_{\rm{}AGN}, + \verb+LastAGNFeedbackScaleFactors+_i \geq{} a(t - 15{\rm{}Myr}) \end{equation} -using the same parameters as used internally by SWIFT and with $a$ the current scale factor. +where $\Delta T_{\rm{}AGN}$ is the same value as used internally by SWIFT, $t$ is the cosmic time, and $a()$ +gives the scale factor as a function of time. diff --git a/documentation/footnote_Xray.tex b/documentation/footnote_Xray.tex index 5d520053..692d2ada 100644 --- a/documentation/footnote_Xray.tex +++ b/documentation/footnote_Xray.tex @@ -4,7 +4,7 @@ order as in the snapshot: \begin{enumerate} - \item eRosita low/soft ($0.2-2.3$~keV) - \item eRosita high/hard ($2.3-8$~keV) + \item eROSITA low/soft ($0.2-2.3$~keV) + \item eROSITA high/hard ($2.3-8$~keV) \item ROSAT ($0.5-2$~keV) \end{enumerate} diff --git a/documentation/footnote_circvel.tex b/documentation/footnote_circvel.tex index 14b13bd9..d16f195e 100644 --- a/documentation/footnote_circvel.tex +++ b/documentation/footnote_circvel.tex @@ -7,9 +7,10 @@ where the cumulative mass $M(\leq{}r)$ includes all particles within the radius $r$, and includes the contribution of the particle(s) at $r=0$. The radius is computed relative to the halo centre. -The softened $v_{\rm{}max}$ value is calculated using the same method, except the particle -radius has a floor of the softening length. An alternative way to calculate $v_{\rm{}max}$ +The softened $v_{\rm{}max}$ value is calculated using the same method, except that the particle +radius has a floor equal to the softening length. An alternative way to calculate $v_{\rm{}max}$ is to estimate it from the halo concentration by assuming an NFW profile. We store the radius of the unsoftened maximum circular velocity. If the softened and unsoftened maximum circular velocities are equal, then their radii will also be equal. If the values are not equal, then the radius of the softened maximum circular velocity will be the simulation softening length. +Note that the softened vmax is $\textit{not}$ the $v_{\rm{}max}$ computed using a softened potential. diff --git a/documentation/footnote_com.tex b/documentation/footnote_com.tex index 746d0894..ba057533 100644 --- a/documentation/footnote_com.tex +++ b/documentation/footnote_com.tex @@ -1,2 +1,2 @@ \paragraph{$^{$FOOTNOTE_NUMBER$}$The centre of mass and centre of mass velocity}\label{footnote:$FOOTNOTE_NUMBER$} are computed using all -particle types except neutrinos (since neutrinos can never be bound to a halo). +particle types except neutrinos. diff --git a/documentation/footnote_compY.tex b/documentation/footnote_compY.tex index 61c42877..89777e0e 100644 --- a/documentation/footnote_compY.tex +++ b/documentation/footnote_compY.tex @@ -1,11 +1,11 @@ -\paragraph{$^{$FOOTNOTE_NUMBER$}$The Compton y parameter}\label{footnote:$FOOTNOTE_NUMBER$} is computed as in McCarthy et al. (2017): +\paragraph{$^{$FOOTNOTE_NUMBER$}$The Compton y parameter}\label{footnote:$FOOTNOTE_NUMBER$} is computed as in \href{https://ui.adsabs.harvard.edu/abs/2017MNRAS.465.2936M}{McCarthy et al. (2017)}: \begin{equation} y \, {d_A}^2(z) = \sum_i \frac{\sigma{}_T}{m_e c^2} n_{e,i} k_B T_{e,i} V_i, \end{equation} -where $d_A(z)$ is the angular diameter distance, $\sigma{}_T$ is the Thomson cross section, $m_e$ the electron mass, $c$ the speed of light and $k_B$ the +where $d_A(z)$ is the angular diameter distance, $\sigma_T$ is the Thomson cross section, $m_e$ the electron mass, $c$ the speed of light and $k_B$ the Boltzmann constant. $n_{e,i}$ and $T_{e,i}$ are the electron number density and electron temperature for gas -particle $i$, while $V_i=m_i/\rho{}_i$ is the SPH volume element that turns the sum over all particles $i$ +particle $i$, while $V_i=m_i/\rho_i$ is the SPH volume element that turns the sum over all particles $i$ within the inclusive sphere into a volume integral. Note that the snapshot already contains the individual -$y_i$ values for the SPH particles, computed from the cooling tables during the simulation. +$y_i$ values for the SPH particles. diff --git a/documentation/footnote_concentration.tex b/documentation/footnote_concentration.tex index c1feee4e..7dcc9284 100644 --- a/documentation/footnote_concentration.tex +++ b/documentation/footnote_concentration.tex @@ -1,12 +1,11 @@ \paragraph{$^{$FOOTNOTE_NUMBER$}$The concentration}\label{footnote:$FOOTNOTE_NUMBER$} is computed using the -method described in Wang et al. (2023), but using a fifth order polynomial fit to -the R1-concentration relation for $1`_. Clicking on each property name will open a dropdown box, which contains information about the dataset within the HDF5 file. The second column gives the filter applied to that property, as descibred in :doc:`property_filters`. The third column indicates the halo variations for which this property is available (green if the property is computed for a certain variation, red if not). The variations are as follows: + +* ``BS`` - :ref:`bound_subhalo_description` +* ``ES`` - :ref:`exclusive_sphere_description` +* ``IS`` - :ref:`inclusive_sphere_description` +* ``EP`` - :ref:`projected_aperture_description` +* ``SO`` - :ref:`spherical_overdensity_description` + +The final column gives a description of the property. Certain properties also contain a link to a footnote at the bottom of this page which gives a full description of how they were calculated. + + From 2561da13e8e0a8df7d356a2e5123be433b7bdba1 Mon Sep 17 00:00:00 2001 From: robjmcgibbon <36136863+robjmcgibbon@users.noreply.github.com> Date: Wed, 27 May 2026 09:24:29 +0100 Subject: [PATCH 16/19] Change attribute reading for numpy 2.4 (#215) * Change attribute reading * Update virgodc version --- .gitignore | 1 + SOAP/catalogue_readers/read_hbtplus.py | 8 ++++---- SOAP/core/halo_tasks.py | 2 +- pyproject.toml | 2 +- requirements.txt | 4 ++-- tests/dummy_halo_generator.py | 5 +++++ 6 files changed, 14 insertions(+), 8 deletions(-) diff --git a/.gitignore b/.gitignore index 914afcd6..707095ce 100644 --- a/.gitignore +++ b/.gitignore @@ -21,3 +21,4 @@ tests/FLAMINGO/test_parameters.yml tests/COLIBRE/test_parameters.yml tests/test_SO_radius_*.png test_data/* +test_SO_radius*png diff --git a/SOAP/catalogue_readers/read_hbtplus.py b/SOAP/catalogue_readers/read_hbtplus.py index 4dff8d93..b54c2047 100644 --- a/SOAP/catalogue_readers/read_hbtplus.py +++ b/SOAP/catalogue_readers/read_hbtplus.py @@ -33,7 +33,7 @@ def read_hbtplus_groupnr(basename, read_potential_energies=False, registry=None) if comm_rank == 0: if os.path.exists(hbt_filename(basename, 0)): with h5py.File(hbt_filename(basename, 0), "r") as infile: - nr_files = int(infile["NumberOfFiles"][...]) + nr_files = infile["NumberOfFiles"][0] sorted_file = False elif os.path.exists(basename): with h5py.File(basename, "r") as infile: @@ -241,9 +241,9 @@ def read_hbtplus_catalogue( have_units = False with h5py.File(filename, "r") as infile: if "Units" in infile: - LengthInMpch = float(infile["Units/LengthInMpch"][...]) - MassInMsunh = float(infile["Units/MassInMsunh"][...]) - VelInKmS = float(infile["Units/VelInKmS"][...]) + LengthInMpch = infile["Units/LengthInMpch"][0] + MassInMsunh = infile["Units/MassInMsunh"][0] + VelInKmS = infile["Units/VelInKmS"][0] have_units = True # Otherwise, will have to read the Parameters.log file if not (have_units): diff --git a/SOAP/core/halo_tasks.py b/SOAP/core/halo_tasks.py index 0261771c..6cb76e52 100644 --- a/SOAP/core/halo_tasks.py +++ b/SOAP/core/halo_tasks.py @@ -354,7 +354,7 @@ def process_halos( next_task.win.Lock(0) next_task.win.Fetch_and_op(one, task_to_do, 0) next_task.win.Unlock(0) - task_to_do = int(task_to_do) + task_to_do = int(task_to_do.flat[0]) # Execute the task, if there's one left if task_to_do < nr_halos: diff --git a/pyproject.toml b/pyproject.toml index 3704b2a2..9d22b7d8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -6,7 +6,7 @@ build-backend = "setuptools.build_meta" [project] name = "SOAP" -version = "0.1.1" +version = "0.1.2" description = "MPI parallel Python code to compute properties of halos in SWIFT n-body simulations" readme = "README.md" requires-python = ">=3.10" diff --git a/requirements.txt b/requirements.txt index 6cf2915d..1ce339de 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,11 +1,11 @@ mpi4py h5py -numpy>=2,<2.4 +numpy>=2 unyt>=3 astropy>=6 scipy matplotlib psutil -virgodc +virgodc>=1.0.3 numba pytest-mpi diff --git a/tests/dummy_halo_generator.py b/tests/dummy_halo_generator.py index 73e36aa8..a726da1f 100644 --- a/tests/dummy_halo_generator.py +++ b/tests/dummy_halo_generator.py @@ -213,6 +213,7 @@ def __init__(self): "BirthTemperatures", "SmoothedElementMassFractions", "IronMassFractionsFromSNIa", + "BirthHaloCatalogueIndex", ], "PartType5": [ "Coordinates", @@ -1170,6 +1171,10 @@ def get_random_halo( data["PartType4"]["GroupNr_all"] = groupnr_all[star_mask] data["PartType4"]["GroupNr_bound"] = groupnr_bound[star_mask] data["PartType4"]["FOFGroupIDs"] = fof_group_ids[star_mask] + # Some in-situ stars, some ex-situ + data["PartType4"]["BirthHaloCatalogueIndex"] = groupnr_bound[star_mask] + mask = np.random.random(Nstar) < 0.1 + data["PartType4"]["BirthHaloCatalogueIndex"][mask] = -1 # initial masses are always larger than the actual mass data["PartType4"]["InitialMasses"] = unyt.unyt_array( mass[star_mask].value * (1.0 + np.random.random(Nstar)), From 23494556d093ccfd58ea875c669e8fb64f2c946f Mon Sep 17 00:00:00 2001 From: robjmcgibbon <36136863+robjmcgibbon@users.noreply.github.com> Date: Wed, 3 Jun 2026 10:02:02 +0100 Subject: [PATCH 17/19] Update location of files for flamingo test (#216) --- SOAP/core/swift_cells.py | 2 +- SOAP/property_table.py | 73 +++++++++++++++-------------- tests/FLAMINGO/parameters_HYDRO.yml | 2 +- 3 files changed, 41 insertions(+), 36 deletions(-) diff --git a/SOAP/core/swift_cells.py b/SOAP/core/swift_cells.py index db42a7ae..62794afa 100644 --- a/SOAP/core/swift_cells.py +++ b/SOAP/core/swift_cells.py @@ -100,7 +100,7 @@ def __call__(self, data, cache): np.s_[mem_start:mem_end, ...], ) except OSError as e: - print(f'Error reading {dataset_name}') + print(f"Error reading {dataset_name}") raise e diff --git a/SOAP/property_table.py b/SOAP/property_table.py index b2bb0446..51a68ceb 100644 --- a/SOAP/property_table.py +++ b/SOAP/property_table.py @@ -4938,11 +4938,15 @@ def add_properties(self, halo_property: HaloProperty, halo_type: str): units = unyt.unyt_quantity(1, units=prop.unit) if not prop.output_physical: units = units * unyt.Unit("a") ** prop.a_scale_exponent - prop_unit = units.units.latex_repr.replace( - "\\rm{km} \\cdot \\rm{kpc}", "\\rm{kpc} \\cdot \\rm{km}" - ).replace( - "\\frac{\\rm{km}^{2}}{\\rm{s}^{2}}", "\\rm{km}^{2} / \\rm{s}^{2}" - ).rstrip() + prop_unit = ( + units.units.latex_repr.replace( + "\\rm{km} \\cdot \\rm{kpc}", "\\rm{kpc} \\cdot \\rm{km}" + ) + .replace( + "\\frac{\\rm{km}^{2}}{\\rm{s}^{2}}", "\\rm{km}^{2} / \\rm{s}^{2}" + ) + .rstrip() + ) prop_dtype = prop.dtype.__name__ @@ -5336,9 +5340,9 @@ def replace_math(match): # Map from internal snap_* unit names to human-readable physical equivalents. # Values are validated against the live unyt registry by _validate_snap_unit_map. _snap_unit_map = { - "snap_mass": "(1e10*Msun)", - "snap_length": "Mpc", - "snap_time": "(s*Mpc/km)", + "snap_mass": "(1e10*Msun)", + "snap_length": "Mpc", + "snap_time": "(s*Mpc/km)", "snap_temperature": "K", } @@ -5359,8 +5363,9 @@ def _resolve_snap_units(self, prop_name: str, prop: dict) -> str: units = unyt.unyt_quantity(1, units=prop_unit) latex = ( - units.units.latex_repr - .replace("\\rm{km} \\cdot \\rm{kpc}", "\\rm{kpc} \\cdot \\rm{km}") + units.units.latex_repr.replace( + "\\rm{km} \\cdot \\rm{kpc}", "\\rm{kpc} \\cdot \\rm{km}" + ) .replace("\\frac{\\rm{km}^{2}}{\\rm{s}^{2}}", "\\rm{km}^{2} / \\rm{s}^{2}") .replace(r"1.0 \times ", "") .rstrip() @@ -5450,10 +5455,7 @@ def get_footnotes_rst(self, name: str) -> tuple[list[int], str]: self.footnotes.append(fnote) footnote_nums.append(i + 1) if footnote_nums: - refs = " ".join( - f"`[{n}] `_" - for n in sorted(footnote_nums) - ) + refs = " ".join(f"`[{n}] `_" for n in sorted(footnote_nums)) return footnote_nums, refs return [], "" @@ -5529,7 +5531,9 @@ def enumerate_to_rst(match): # Convert inline math and non-breaking spaces within the # item text now, before the block is stashed as a # placeholder (the main conversion steps run before restore) - item = re.sub(r"\$([^\n\$]+?)\$", lambda m: f":math:`{m.group(1)}`", item) + item = re.sub( + r"\$([^\n\$]+?)\$", lambda m: f":math:`{m.group(1)}`", item + ) item = item.replace("~", " ") rst_items.append(f"#. {item}") block = "\n\n" + "\n".join(rst_items) + "\n\n" @@ -5640,24 +5644,24 @@ def _build_rst_table(self, prop_names: list, lines: list): # swiftsimio name snake_name = ".".join( - self._camel_to_snake(part) for part in prop['name'].split("/") + self._camel_to_snake(part) for part in prop["name"].split("/") ) # HDF5 output name (prepend InputHalos/ where applicable) output_name = prop["name"] if output_name.split("/")[0] in ("HBTplus", "VR", "FOF"): output_name = "InputHalos/" + output_name - snake_name = 'input_halos_' + snake_name + snake_name = "input_halos_" + snake_name prop_units_rst = self._resolve_snap_units(prop_name, prop) prop_comp_rst = self._compression_to_rst( self.compression_description[prop["compression"]] ) - prop_filter = prop['category'] + prop_filter = prop["category"] output_types = self._get_output_types_rst(prop) - if 'DummyProperties' in prop['types']: - prop_filter = 'basic' - output_types = r'\-' + if "DummyProperties" in prop["types"]: + prop_filter = "basic" + output_types = r"\-" description = prop["description"].format( label="satisfying a spherical overdensity criterion.", core_excision="excised core", @@ -5713,9 +5717,11 @@ def generate_rst_files(self, output_dir: str): prop_names_sorted = sorted( self.properties.keys(), key=lambda key: ( - category_order.index(self.properties[key]["category"]) - if self.properties[key]["category"] in category_order - else len(category_order), + ( + category_order.index(self.properties[key]["category"]) + if self.properties[key]["category"] in category_order + else len(category_order) + ), self.properties[key]["name"].lower(), ), ) @@ -5726,11 +5732,11 @@ def generate_rst_files(self, output_dir: str): dmo_props = [] hydro_props = [] for n in prop_names_sorted: - if self.properties[n]['category'] == 'InputHalos': + if self.properties[n]["category"] == "InputHalos": input_basic_props.append(n) - elif 'DummyProperties' in self.properties[n]['types']: + elif "DummyProperties" in self.properties[n]["types"]: input_copied_props.append(n) - elif self.properties[n]['dmo']: + elif self.properties[n]["dmo"]: dmo_props.append(n) else: hydro_props.append(n) @@ -5740,16 +5746,16 @@ def generate_rst_files(self, output_dir: str): self.footnotes = [] # Introduction - with open(f'documentation/property_table_intro.rst') as file: + with open(f"documentation/property_table_intro.rst") as file: lines = [line.rstrip() for line in file.readlines()] # Tables for props, table_title in [ - (input_basic_props, "Input halo properties"), - (dmo_props, "Dark matter only properties"), - (hydro_props, "Hydrodynamical properties"), - (input_copied_props, "Copied properties"), - ]: + (input_basic_props, "Input halo properties"), + (dmo_props, "Dark matter only properties"), + (hydro_props, "Hydrodynamical properties"), + (input_copied_props, "Copied properties"), + ]: lines.append(table_title) lines.append("-" * len(table_title)) @@ -5778,7 +5784,6 @@ def generate_rst_files(self, output_dir: str): ofile.write("\n".join(lines)) - class DummyProperties: """ Dummy HaloProperty object used to include properties which are not computed diff --git a/tests/FLAMINGO/parameters_HYDRO.yml b/tests/FLAMINGO/parameters_HYDRO.yml index cf7fa252..c57e712d 100644 --- a/tests/FLAMINGO/parameters_HYDRO.yml +++ b/tests/FLAMINGO/parameters_HYDRO.yml @@ -16,7 +16,7 @@ GroupMembership: filename: "{sim_dir}/{sim_name}/SOAP-HBT/membership_{snap_nr:04d}/membership_{snap_nr:04d}.{file_nr}.hdf5" ExtraInput: - xrays: "/cosma8/data/dp004/dc-mcgi1/FLAMINGO/Xray/{sim_name}/xray/flamingo_{snap_nr:04}/xray_{snap_nr:04}.{file_nr}.hdf5" + xrays: "/cosma8/data/dp004/dc-mcgi1/SOAP/TEST_DATA/Xray/{sim_name}/xray/flamingo_{snap_nr:04}/xray_{snap_nr:04}.{file_nr}.hdf5" HaloProperties: # Where to write the halo properties file From 9ff3902ccb7fed00819c1e5bb4d95123f6e42208 Mon Sep 17 00:00:00 2001 From: robjmcgibbon <36136863+robjmcgibbon@users.noreply.github.com> Date: Tue, 16 Jun 2026 10:29:37 +0100 Subject: [PATCH 18/19] Clean up EAGLE script (#217) * Handle DM in mass table * Various fixes * Format --- SOAP/catalogue_readers/read_subfind_eagle.py | 3 +- SOAP/compression/make_virtual_snapshot.py | 2 + misc/convert_eagle.py | 132 ++++++++++++------- parameter_files/EAGLE.yml | 2 +- scripts/EAGLE.sh | 68 ++++++++-- 5 files changed, 142 insertions(+), 65 deletions(-) diff --git a/SOAP/catalogue_readers/read_subfind_eagle.py b/SOAP/catalogue_readers/read_subfind_eagle.py index 31ff7fb0..4bc86ec0 100644 --- a/SOAP/catalogue_readers/read_subfind_eagle.py +++ b/SOAP/catalogue_readers/read_subfind_eagle.py @@ -117,7 +117,8 @@ def read_subfind_catalogue(comm, basename, a_unit, registry, boxsize): ) # Store initial search radius - search_radius = (5 * data["Subhalo/VmaxRadius"] / h) * swift_cmpc + search_radius_cmpc = np.minimum((5 * data["Subhalo/VmaxRadius"] / h), 5) + search_radius = search_radius_cmpc * swift_cmpc local_halo = { "cofp": cofp, diff --git a/SOAP/compression/make_virtual_snapshot.py b/SOAP/compression/make_virtual_snapshot.py index 2311840e..ea30ea1c 100644 --- a/SOAP/compression/make_virtual_snapshot.py +++ b/SOAP/compression/make_virtual_snapshot.py @@ -360,3 +360,5 @@ def replace_path(old_path): absolute_paths=args.absolute_paths, discard_duplicate_datasets=args.discard_duplicate_datasets, ) + + print(f"Done!") diff --git a/misc/convert_eagle.py b/misc/convert_eagle.py index 3f818361..cd468169 100644 --- a/misc/convert_eagle.py +++ b/misc/convert_eagle.py @@ -8,11 +8,15 @@ mpirun -- python convert_eagle.py \ --snapshot-basename=SNAPSHOT \ + --particledata-basename=PARTICLE_DATA \ + --subfind-basename=SUBFIND \ --output-basename=OUTPUT \ --membership-basename=MEMBERSHIP -where SNAPSHOT is the EAGLE snapshot (use the particledata_*** files, -since the normal snapshots don't store SubGroupNumber), and OUTPUT & +where SNAPSHOT is the EAGLE snapshot (the snapshots_*** files), +PARTICLE_DATA are the EAGLE membership files (the particledata_*** files, +since the normal snapshots don't store SubGroupNumber), SUBFIND are +the SubFind catalogues (the subfind_tab_*** files), and OUTPUT & MEMBERSHIP are the names of the output files. You must run with the same number of ranks as input files. @@ -82,11 +86,21 @@ "name without the .{file_nr}.hdf5 suffix)" ), ) +parser.add_argument( + "--particledata-basename", + type=str, + required=True, + help=( + "The basename for the particle data files (the files which " + "contain the GroupNumber and SubGroupNumber values for bound " + "particles. For EAGLE these are separate to the snapshots." + ), +) parser.add_argument( "--subfind-basename", type=str, required=True, - help=("The basename for the subfind files"), + help="The basename for the subfind files", ) parser.add_argument( "--output-basename", @@ -102,6 +116,7 @@ ) args = parser.parse_args() snap_filename = args.snap_basename + ".{file_nr}.hdf5" +particledata_filename = args.particledata_basename + ".{file_nr}.hdf5" subfind_filename = args.subfind_basename + ".{file_nr}.hdf5" output_filename = args.output_basename + ".{file_nr}.hdf5" membership_filename = args.membership_basename + ".{file_nr}.hdf5" @@ -132,7 +147,7 @@ h = comm.bcast(h) box_size_cmpc = comm.bcast(box_size_cmpc) -assert comm_size == n_file +assert comm_size <= n_file # Specify the unit system of the output SWIFT snapshot if comm_rank == 0: @@ -281,13 +296,6 @@ "description": None, "conversion_factor": None, }, - "GroupNumber": { - "swift_name": "FOFGroupIDs", - "exponents": {"L": 0, "M": 0, "T": 0, "t": 0}, - "a_exponent": None, - "description": None, - "conversion_factor": None, - }, "ParticleIDs": { "swift_name": "ParticleIDs", "exponents": {"L": 0, "M": 0, "T": 0, "t": 0}, @@ -380,13 +388,6 @@ "description": "Particle mass", "conversion_factor": None, }, - "GroupNumber": { - "swift_name": "FOFGroupIDs", - "exponents": {"L": 0, "M": 0, "T": 0, "t": 0}, - "a_exponent": None, - "description": None, - "conversion_factor": None, - }, "ParticleIDs": { "swift_name": "ParticleIDs", "exponents": {"L": 0, "M": 0, "T": 0, "t": 0}, @@ -417,13 +418,6 @@ "description": None, "conversion_factor": None, }, - "GroupNumber": { - "swift_name": "FOFGroupIDs", - "exponents": {"L": 0, "M": 0, "T": 0, "t": 0}, - "a_exponent": None, - "description": None, - "conversion_factor": None, - }, "ParticleIDs": { "swift_name": "ParticleIDs", "exponents": {"L": 0, "M": 0, "T": 0, "t": 0}, @@ -506,13 +500,6 @@ "description": None, "conversion_factor": None, }, - "GroupNumber": { - "swift_name": "FOFGroupIDs", - "exponents": {"L": 0, "M": 0, "T": 0, "t": 0}, - "a_exponent": None, - "description": None, - "conversion_factor": None, - }, "ParticleIDs": { "swift_name": "ParticleIDs", "exponents": {"L": 0, "M": 0, "T": 0, "t": 0}, @@ -569,14 +556,13 @@ ] = conversion_factor # DM mass can be a special case + properties["dm_mass_in_table"] = False if "Mass" in properties.get(f"PartType1", {}): if "Mass" not in infile["PartType1"]: # Load DM mass from mass table dm_mass = infile["Header"].attrs["MassTable"][1] / h properties["PartType1"]["Mass"]["conversion_factor"] = dm_mass - else: - # Treat DM mass as any other property - dm_mass = 0 + properties["dm_mass_in_table"] = True # Get list of elements for ElementMassFractions if "ElementMassFractions" in properties.get(f"PartType0", {}): @@ -608,6 +594,9 @@ snap_file = phdf5.MultiFile( snap_filename, file_nr_attr=("Header", "NumFilesPerSnapshot"), comm=comm ) +particledata_file = phdf5.MultiFile( + particledata_filename, file_nr_attr=("Header", "NumFilesPerSnapshot"), comm=comm +) # Load the SubFind catalogue and create an array that links the GroupNumber # and SubGroupNumber of a subhalo to its index within the subhalo catalogue # (for creating the membership files) @@ -657,6 +646,15 @@ elements_per_file[1:] -= elements_per_file[:-1] assert np.sum(elements_per_file) == np.sum(cell_counts[ptype]) + # Each rank writes the files assigned to it by MultiFile, so it + # must hold the concatenation of the data of those files. + elements_per_rank = np.zeros(comm_size, dtype=elements_per_file.dtype) + for rank in range(comm_size): + first = snap_file.first_file_on_rank[rank] + num = snap_file.num_files_on_rank[rank] + elements_per_rank[rank] = np.sum(elements_per_file[first : first + num]) + assert np.sum(elements_per_rank) == np.sum(elements_per_file) + # Calculate offsets of the first particle in each cell cell_files[ptype] = np.repeat(np.arange(n_file), cells_per_file) absolute_offset = np.cumsum(cell_counts[ptype]) - cell_counts[ptype] @@ -674,7 +672,7 @@ if comm_rank == 0: print(f"Converting PartType{ptype}/{prop}") - if (ptype == 1) and (prop == "Mass") and (dm_mass == 0): + if (ptype == 1) and (prop == "Mass") and properties["dm_mass_in_table"]: # DM particles all have the same mass, so are not saved in the snapshots arr = np.ones(pos.shape[0]) else: @@ -704,7 +702,7 @@ attrs.update(unit_attrs) # Write to the output file - arr = psort.repartition(arr, elements_per_file, comm=comm) + arr = psort.repartition(arr, elements_per_rank, comm=comm) if create_output_file: mode = "w" create_output_file = False @@ -751,7 +749,7 @@ attrs.update(unit_attrs) # Write to the output file - arr = psort.repartition(arr, elements_per_file, comm=comm) + arr = psort.repartition(arr, elements_per_rank, comm=comm) if create_output_file: mode = "w" create_output_file = False @@ -766,44 +764,78 @@ attrs={"ElementMassFractions": attrs}, ) + # Load the GroupNumber and SubGroupNumber of particles by matching + # the snapshot_* files with the particledata_* files + snap_ids = snap_file.read(f"PartType{ptype}/ParticleIDs") + particledata_ids = particledata_file.read(f"PartType{ptype}/ParticleIDs") + idx = psort.parallel_match(snap_ids, particledata_ids, comm=comm) + + # EAGLE uses a value of 2^30 to indicate unbound particles + # Particles missing from the particledata_* files are always unbound + particledata_sub_group_nr = particledata_file.read( + f"PartType{ptype}/SubGroupNumber" + ) + sub_group_nr = 1073741824 * np.ones(snap_ids.shape[0], dtype=np.int32) + sub_group_nr[idx != -1] = psort.fetch_elements( + particledata_sub_group_nr, + idx[idx != -1], + comm=comm, + ) + + # Negative values indicate that the particle is not part of a FoF, + # but is within the SO group of a FoF (the FoF it is part of is the positive + # value, e.g. -10 means it is within the SO of FoF 10) + particledata_group_nr = particledata_file.read(f"PartType{ptype}/GroupNumber") + particledata_group_nr[particledata_group_nr < 0] = 1073741824 + group_nr = 1073741824 * np.ones(snap_ids.shape[0], dtype=np.int32) + group_nr[idx != -1] = psort.fetch_elements( + particledata_group_nr, + idx[idx != -1], + comm=comm, + ) + # Create a subhalo id for each particle by combining the # group number and subgroup number - sub_group = snap_file.read(f"PartType{ptype}/SubGroupNumber") - subhalo = snap_file.read(f"PartType{ptype}/GroupNumber").astype(np.int64) + subhalo = group_nr.astype(np.int64) subhalo <<= 32 - subhalo += sub_group.astype(np.int64) - # Indicate unbound particles with -1 - bound = sub_group != 1073741824 + subhalo += sub_group_nr.astype(np.int64) + # For SOAP we want unbound particles to be indicated with a value of -1 + bound = sub_group_nr != 1073741824 subhalo[np.logical_not(bound)] = -1 # Get SubFind index of bound particles subhalo[bound] = psort.parallel_match(subhalo[bound], subfind_id, comm=comm) assert np.all(subhalo[bound] != -1) - # Sort, add units, and write to file (same as for other properties) + # Sort values spatially (same as for other properties) subhalo = psort.fetch_elements(subhalo, order, comm=comm) + group_nr = psort.fetch_elements(group_nr, order, comm=comm) + # Add units, and write to file units = unyt.Unit("dimensionless", registry=reg) unit_attrs = swift_units.attributes_from_units(units, False, 0) - attrs = { + subhalo_attrs = { "Description": ( "Unique identifier of the subhalo this particle is " "bound to. This is a combination of the GroupNumber and" "the SubGroupNumber. -1 if the particle is not bound" ) } - attrs.update(unit_attrs) - subhalo = psort.repartition(subhalo, elements_per_file, comm=comm) + subhalo_attrs.update(unit_attrs) + fof_attrs = {"Description": "FoF group number particle is in"} + fof_attrs.update(unit_attrs) + subhalo = psort.repartition(subhalo, elements_per_rank, comm=comm) + group_nr = psort.repartition(group_nr, elements_per_rank, comm=comm) if create_membership_file: mode = "w" create_membership_file = False else: mode = "r+" snap_file.write( - {"GroupNr_bound": subhalo}, + {"GroupNr_bound": subhalo, "FOFGroupIDs": group_nr}, elements_per_file, filenames=membership_filename, mode=mode, group=f"PartType{ptype}", - attrs={"GroupNr_bound": attrs}, + attrs={"GroupNr_bound": subhalo_attrs, "FOFGroupIDs": fof_attrs}, ) # Add headers to the snapshots @@ -814,7 +846,7 @@ header = outfile.create_group("Header") for name, value in swift_header.items(): header.attrs[name] = value - n_part = np.zeros(max(ptypes) + 1) + n_part = np.zeros(max(ptypes) + 1, dtype=np.int64) for ptype in ptypes: n_part[ptype] = outfile[f"PartType{ptype}/Coordinates"].shape[0] header.attrs["NumPart_ThisFile"] = n_part diff --git a/parameter_files/EAGLE.yml b/parameter_files/EAGLE.yml index 5a1d05d8..2b825167 100644 --- a/parameter_files/EAGLE.yml +++ b/parameter_files/EAGLE.yml @@ -9,7 +9,7 @@ Parameters: # Location of the Swift snapshots: Snapshots: # Use {snap_nr:04d} for the snapshot number and {file_nr} for the file number. - filename: "{sim_dir}/{sim_name}/swift_snapshots/swift_{snap_nr:03d}/snap_{snap_nr:03d}.{file_nr}.hdf5" + filename: "{sim_dir}/{sim_name}/snapshots/snap_{snap_nr:03d}/snap_{snap_nr:03d}.{file_nr}.hdf5" # Which halo finder we're using, and base name for halo finder output files HaloFinder: diff --git a/scripts/EAGLE.sh b/scripts/EAGLE.sh index 72bc0511..3ae44945 100755 --- a/scripts/EAGLE.sh +++ b/scripts/EAGLE.sh @@ -9,7 +9,9 @@ #SBATCH --exclusive #SBATCH -t 02:00:00 # -# For L0025N0752 set ntasks=16 +# For L0025N0752/REFERENCE set ntasks=16 +# For L0025N0752/RECALIBRATED set ntasks=32 +# For L0025N0752/WDM set ntasks=32 # For L0100N1504 set ntasks=256 # # Install Hdecompose with: @@ -20,7 +22,18 @@ # # Submit this script from the main SOAP directory ($ sbatch scripts/EAGLE.sh) -sim_name='L0100N1504' +# sim_name='L0025N0752/PE/REFERENCE' +# sim_name='L0025N0752/PE/RECALIBRATED' +sim_name='L0100N1504/PE/REFERENCE' +# sim_name='L0025N0752/EAGLE_WDM' + +sim_dir="/cosma7/data/Eagle/ScienceRuns/Planck1/${sim_name}/data" +# sim_dir="/snap7/scratch/dp004/dc-mcgi1/EAGLE_WDM/data" + +# Note that if you update the output directory, you will +# also need to update the SOAP parameter file +output_dir="/snap7/scratch/dp004/dc-mcgi1/SOAP_EAGLE/${sim_name}" + snap_nr="028" z_suffix="z000p000" @@ -30,13 +43,19 @@ source openmpi-5.0.3-hdf5-1.12.3-env/bin/activate ######## Link files to snap (to remove awful z suffix) -sim_dir="/cosma7/data/Eagle/ScienceRuns/Planck1/${sim_name}/PE/REFERENCE/data" -# Note that if you update the output directory, you will also need to -# update the SOAP parameter file -output_dir="/snap7/scratch/dp004/dc-mcgi1/SOAP_EAGLE/${sim_name}" +sim_snap_dir="${sim_dir}/snapshot_${snap_nr}_${z_suffix}" +output_snap_dir="${output_dir}/gadget_snapshots/snapshot_${snap_nr}" +mkdir -p $output_snap_dir +i=0 +while [[ -e "${sim_snap_dir}/snap_${snap_nr}_${z_suffix}.${i}.hdf5" ]]; do + old_name="${sim_snap_dir}/snap_${snap_nr}_${z_suffix}.${i}.hdf5" + new_name="${output_snap_dir}/snap_${snap_nr}.${i}.hdf5" + ln -s $old_name $new_name + ((i++)) +done sim_snap_dir="${sim_dir}/particledata_${snap_nr}_${z_suffix}" -output_snap_dir="${output_dir}/gadget_snapshots/snapshot_${snap_nr}" +output_snap_dir="${output_dir}/gadget_membership/snapshot_${snap_nr}" mkdir -p $output_snap_dir i=0 while [[ -e "${sim_snap_dir}/eagle_subfind_particles_${snap_nr}_${z_suffix}.${i}.hdf5" ]]; do @@ -63,28 +82,41 @@ set -e mpirun -- python -u misc/convert_eagle.py \ --snap-basename "${output_dir}/gadget_snapshots/snapshot_${snap_nr}/snap_${snap_nr}" \ + --particledata-basename "${output_dir}/gadget_membership/snapshot_${snap_nr}/snap_${snap_nr}" \ --subfind-basename "${output_group_dir}/subfind_tab_${snap_nr}" \ - --output-basename "${output_dir}/swift_snapshots/swift_${snap_nr}/snap_${snap_nr}" \ + --output-basename "${output_dir}/snapshots/snap_${snap_nr}/snap_${snap_nr}" \ --membership-basename "${output_dir}/SOAP_uncompressed/membership_${snap_nr}/membership_${snap_nr}" ######### Estimate SpeciesFraction of hydrogen mpirun -- python -u misc/hdecompose_hydrogen_fractions.py \ - --snap-basename "${output_dir}/swift_snapshots/swift_${snap_nr}/snap_${snap_nr}" \ + --snap-basename "${output_dir}/snapshots/snap_${snap_nr}/snap_${snap_nr}" \ --output-basename "${output_dir}/species_fractions/swift_${snap_nr}/snap_${snap_nr}" ######### Create virtual snapshot # Must be run from the snapshot directory itself or there will be issues with paths soap_dir=$(pwd) -cd "${output_dir}/swift_snapshots/swift_${snap_nr}" +cd "${output_dir}/snapshots/snap_${snap_nr}" python "${soap_dir}/create_virtual_snapshot.py" "snap_${snap_nr}.0.hdf5" cd - +######### Repack memebership files and create virtual snapshot + +input_filename="${output_dir}/SOAP_uncompressed/membership_${snap_nr}/membership_${snap_nr}" +output_filename="${output_dir}/SOAP/membership_${snap_nr}/membership_${snap_nr}" +mkdir -p "${output_dir}/SOAP/membership_${snap_nr}" + +nr_files=`ls -1 ${input_filename}.*.hdf5 | wc -l` +nr_files_minus_one=$(( ${nr_files} - 1 )) + +seq 0 ${nr_files_minus_one} | xargs -I {} -P 16 bash -c \ + "h5repack -i ${input_filename}.{}.hdf5 -o ${output_filename}.{}.hdf5 -l CHUNK=10000 -f GZIP=4" + python SOAP/compression/make_virtual_snapshot.py \ - --virtual-snapshot "${output_dir}/swift_snapshots/swift_${snap_nr}/snap_${snap_nr}.hdf5" \ - --auxiliary-snapshots "${output_dir}/SOAP_uncompressed/membership_${snap_nr}/membership_${snap_nr}.{file_nr}.hdf5" \ - --output-file "${output_dir}/SOAP_uncompressed/snap_${snap_nr}.hdf5" + --virtual-snapshot "${output_dir}/snapshots/snap_${snap_nr}/snap_${snap_nr}.hdf5" \ + --auxiliary-snapshots "${output_dir}/SOAP/membership_${snap_nr}/membership_${snap_nr}.{file_nr}.hdf5" "${output_dir}/species_fractions/swift_${snap_nr}/snap_${snap_nr}.{file_nr}.hdf5" \ + --output-file "${output_dir}/SOAP/snap_with_SOAP_membership_${snap_nr}.hdf5" ######### Run SOAP @@ -94,6 +126,16 @@ mpirun -- python3 -u -m mpi4py SOAP/compute_halo_properties.py \ parameter_files/EAGLE.yml \ --sim-name=${sim_name} --snap-nr=${snap_nr} --chunks=${chunks} +######### Compress SOAP (run on single node) + +module purge +module load python/3.12.4 + +mpirun -np 28 python -u SOAP/compression/compress_soap_catalogue.py \ + "${output_dir}/SOAP_uncompressed/halo_properties_${snap_nr}.hdf5" \ + "${output_dir}/SOAP/halo_properties_${snap_nr}.hdf5" \ + "${output_dir}/SOAP_compression_tmp" + ############## echo "Job complete!" From 15e5299ae48cf095c35653f3f450e71b22328089 Mon Sep 17 00:00:00 2001 From: robjmcgibbon <36136863+robjmcgibbon@users.noreply.github.com> Date: Fri, 14 Aug 2026 13:02:55 +0100 Subject: [PATCH 19/19] Pass halo indices from a file (#218) * Flush errors * Pass halo indices as a file --- README.md | 61 +++- SOAP/catalogue_readers/read_hbtplus.py | 2 +- SOAP/compression/compress_soap_catalogue.py | 2 +- .../create_empty_SOAP_catalogue.py | 1 - SOAP/compute_halo_properties.py | 8 +- SOAP/core/chunk_tasks.py | 2 +- SOAP/core/combine_args.py | 34 ++ SOAP/core/combine_chunks.py | 5 - SOAP/core/halo_centres.py | 38 ++- SOAP/core/soap_args.py | 120 +++++++- SOAP/core/swift_cells.py | 5 +- SOAP/group_membership.py | 2 +- misc/recalculate_xrays.py | 2 +- tests/COLIBRE/find_halo_ids.py | 12 +- tests/COLIBRE/halo_indices_0092.txt | 290 ++++++++++++++++++ tests/COLIBRE/run_L0025N0188_Thermal.sh | 7 +- tests/test_halo_indices_file.py | 64 ++++ tests/test_read_subfind.py | 2 +- tests/test_shared_mesh.py | 5 +- 19 files changed, 608 insertions(+), 54 deletions(-) create mode 100644 tests/COLIBRE/halo_indices_0092.txt create mode 100644 tests/test_halo_indices_file.py diff --git a/README.md b/README.md index 48e1fc96..add6a57f 100644 --- a/README.md +++ b/README.md @@ -102,14 +102,66 @@ The optional `--max-ranks-reading` flag determines how many MPI ranks per node read the snapshot. This can be used to avoid overloading the file system. The default value is 32. +### Selecting which subhalos to process + +By default SOAP calculates properties for every subhalo in the input catalogue. +There are a number of flags which can be used to process only a subset of them. + +The `--centrals-only` flag discards satellites, so that only central subhalos are +processed. + +Individual subhalos can be selected with the `--halo-indices` flag. +This specifies the index of the required subhalos in the halo catalogue, which is +the quantity written to `InputHalos/HaloCatalogueIndex`. + +For larger numbers of subhalos the indices can be listed in a text file, which is +passed with the `--halo-indices-file` flag. +The file must contain one index per line. Blank lines, and anything following a +`#`, are ignored. Duplicate indices are discarded. + +### Command line arguments + +The arguments listed here are passed on the command line, and cannot be set in +the parameter file. + +`SOAP/group_membership.py`: + +| Argument | Description | +| --- | --- | +| `config_file` | Name of the yaml parameter file. Required | +| `--sim-name` | Name of the simulation to process | +| `--snap-nr` | Snapshot number to process | + +`SOAP/compute_halo_properties.py`: + +| Argument | Default | Description | +| --- | --- | --- | +| `config_file` | | Name of the yaml parameter file. Required | +| `--sim-name` | | Name of the simulation to process | +| `--snap-nr` | | Snapshot number to process | +| `--chunks` | 1 | Number of chunks to split the volume into. Should be at least the number of compute nodes | +| `--dmo` | off | Run in dark matter only mode, skipping any hydro-only properties | +| `--centrals-only` | off | Only process central halos, discarding satellites. See [Selecting which subhalos to process](#selecting-which-subhalos-to-process) | +| `--halo-indices` | | Only process the listed halo indices. See [Selecting which subhalos to process](#selecting-which-subhalos-to-process) | +| `--halo-indices-file` | | Only process the halo indices listed in the given file. See [Selecting which subhalos to process](#selecting-which-subhalos-to-process) | +| `--max-halos` | 0 (all) | Only process the first N halos in the catalogue. See [Debugging](#debugging) | +| `--record-halo-timings` | off | Record the time taken to process each halo. See [Timing](#timing) | +| `--record-property-timings` | off | Record the time taken to calculate each property. This doubles the size of the output catalogue. See [Timing](#timing) | +| `--reference-snapshot` | | Number of a snapshot which contains all particle types. Used to determine the datasets and units of any particle types which are missing from the snapshot being processed, e.g. stars or black holes at high redshift | +| `--snipshot` / `--snapshot` | auto | Force snipshot or snapshot mode. By default this is determined from the value of `SelectOutput` in the snapshot header | +| `--profile` | 0 | Run with profiling. See [Profiling](#profiling) | +| `--max-ranks-reading` | 32 | Number of MPI ranks per node which read snapshot data. Can be reduced to avoid overloading the file system | +| `--output-parameters` | | Where to write the parameters used by this run, in yaml format | + ### Parameter files To run either of the programs a parameters file must be passed. This contains information including the input and output directories, the halo finder to use, which halo definitions to use, and which properties to calculate for each halo definition. A description -of all possible fields, and a number of example parameter files -can be found in the `parameters_files` directory. +of all possible fields can be found in +[`parameter_files/README.md`](parameter_files/README.md), alongside a number +of example parameter files. ### Compression @@ -194,8 +246,9 @@ mpirun. The `-Werror` flag is useful for making pdb stop on warnings. E.g. division by zero in the halo property calculations will be caught. -It is also possible to select individual halos to process with the `--halo-indices` -flag. This specifies the index of the required halos in the halo catalogue. E.g. +If SOAP crashes while processing a particular halo it will try to report the +index of that halo, which can then be re-run on its own with the +`--halo-indices` flag, e.g. ``` python3 -Werror -m pdb ./compute_halo_properties.py --halo-indices 1 2 3 ... ``` diff --git a/SOAP/catalogue_readers/read_hbtplus.py b/SOAP/catalogue_readers/read_hbtplus.py index b54c2047..6e0e8b65 100644 --- a/SOAP/catalogue_readers/read_hbtplus.py +++ b/SOAP/catalogue_readers/read_hbtplus.py @@ -41,7 +41,7 @@ def read_hbtplus_groupnr(basename, read_potential_energies=False, registry=None) nr_files = 1 sorted_file = True else: - print(f"No HBT files found for basename {basename}") + print(f"No HBT files found for basename {basename}", flush=True) comm.Abort() else: nr_files = None diff --git a/SOAP/compression/compress_soap_catalogue.py b/SOAP/compression/compress_soap_catalogue.py index c19dadf1..6114f419 100644 --- a/SOAP/compression/compress_soap_catalogue.py +++ b/SOAP/compression/compress_soap_catalogue.py @@ -206,7 +206,7 @@ def assign_datasets(nr_files, nr_ranks, comm_rank): datasets = h5copy.dsets.copy() except Exception as e: - print(f"Error: {e}") + print(f"Error: {e}", flush=True) comm.Abort(1) else: tmp_dir = None diff --git a/SOAP/compression/create_empty_SOAP_catalogue.py b/SOAP/compression/create_empty_SOAP_catalogue.py index 17079f8e..f3367713 100644 --- a/SOAP/compression/create_empty_SOAP_catalogue.py +++ b/SOAP/compression/create_empty_SOAP_catalogue.py @@ -115,7 +115,6 @@ def __call__(self, name, h5obj): elif name == "Parameters": for attr in self.ifile[name].attrs: self.ofile[name].attrs[attr] = self.ifile[name].attrs[attr] - self.ofile[name].attrs["halo_indices"] = np.array([], dtype="int64") self.ofile[name].attrs["snapshot_nr"] = self.snapnum else: for attr in self.ifile[name].attrs: diff --git a/SOAP/compute_halo_properties.py b/SOAP/compute_halo_properties.py index 46e72258..61f19aab 100644 --- a/SOAP/compute_halo_properties.py +++ b/SOAP/compute_halo_properties.py @@ -126,7 +126,7 @@ def compute_halo_properties(): swift_filename, extra_input, swift_filename_ref, extra_input_ref ) except Exception as err_msg: - print(err_msg) + print(err_msg, flush=True) # Thrown if there are issues with the input files comm_world.Abort(1) parsec_cgs = cellgrid.constants["parsec"] @@ -234,7 +234,7 @@ def compute_halo_properties(): # We require BoundSubhalo since it's used for filters if comm_world_rank == 0: if "SubhaloProperties" not in parameter_file.parameters: - print("SubhaloProperties must be in the parameter file") + print("SubhaloProperties must be in the parameter file", flush=True) comm_world.Abort(1) halo_prop_list.append( subhalo_properties.SubhaloProperties( @@ -544,7 +544,7 @@ def compute_halo_properties(): try: os.makedirs(os.path.dirname(args.output_file), exist_ok=True) except OSError as e: - print(f"Error creating output directory: {e}") + print(f"Error creating output directory: {e}", flush=True) comm_world.Abort(1) comm_world.barrier() @@ -586,7 +586,7 @@ def compute_halo_properties(): try: os.makedirs(scratch_file_dir, exist_ok=True) except OSError as e: - print(f"Error creating scratch directory: {e}") + print(f"Error creating scratch directory: {e}", flush=True) comm_world.Abort(1) comm_world.barrier() diff --git a/SOAP/core/chunk_tasks.py b/SOAP/core/chunk_tasks.py index ae6e7fb3..986deae4 100644 --- a/SOAP/core/chunk_tasks.py +++ b/SOAP/core/chunk_tasks.py @@ -244,7 +244,7 @@ def message(m): try: cellgrid.check_datasets_exist(properties, self.halo_prop_list) except KeyError as err_msg: - print(err_msg) + print(err_msg, flush=True) comm.Abort(1) else: properties = None diff --git a/SOAP/core/combine_args.py b/SOAP/core/combine_args.py index 8f06f184..89407131 100644 --- a/SOAP/core/combine_args.py +++ b/SOAP/core/combine_args.py @@ -4,6 +4,30 @@ from virgo.util.partial_formatter import PartialFormatter +# Arguments which must be passed on the command line. The Parameters section +# of the config file is intended for values which are substituted into the +# other sections, so we don't allow these to be set there. Note that most of +# them would be silently ignored if they were, since command line arguments +# which have a default value always take precedence over the config file. +COMMAND_LINE_ONLY_PARAMETERS = frozenset( + ( + "config_file", + "chunks", + "dmo", + "centrals_only", + "record_halo_timings", + "record_property_timings", + "max_halos", + "halo_indices", + "halo_indices_file", + "profile", + "max_ranks_reading", + "output_parameters", + "snipshot", + "snapshot", + ) +) + def combine_arguments(command_line_args, config_file): """ @@ -21,6 +45,16 @@ def combine_arguments(command_line_args, config_file): with open(config_file, "r") as infile: config_file_args = yaml.safe_load(infile) + # Check the config file doesn't set arguments which must be passed + # on the command line + invalid = COMMAND_LINE_ONLY_PARAMETERS.intersection(config_file_args["Parameters"]) + if invalid: + raise ValueError( + "The following cannot be set in the Parameters section of the " + f"config file, they must be passed on the command line: " + f"{', '.join(sorted(invalid))}" + ) + # Combine the two all_args = {"Parameters": {}} for name in config_file_args["Parameters"]: diff --git a/SOAP/core/combine_chunks.py b/SOAP/core/combine_chunks.py index 63079d69..90d5478f 100644 --- a/SOAP/core/combine_chunks.py +++ b/SOAP/core/combine_chunks.py @@ -224,11 +224,6 @@ def combine_chunks( params.attrs["centrals_only"] = 0 if args.centrals_only == False else 1 calc_names = sorted([hp.name for hp in halo_prop_list]) params.attrs["calculations"] = calc_names - params.attrs["halo_indices"] = ( - args.halo_indices - if args.halo_indices is not None - else np.ndarray(0, dtype=int) - ) if recently_heated_gas_filter.initialised: recently_heated_gas_metadata = recently_heated_gas_filter.get_metadata() recently_heated_gas_params = params.create_group( diff --git a/SOAP/core/halo_centres.py b/SOAP/core/halo_centres.py index a7e49b2c..a5538c75 100644 --- a/SOAP/core/halo_centres.py +++ b/SOAP/core/halo_centres.py @@ -106,14 +106,36 @@ def __init__( del halo_data # Only keep halos in the supplied list of halo IDs. - if (args.halo_indices is not None) and (local_halo["index"].shape[0]): + if args.halo_indices is not None: halo_indices = np.asarray(args.halo_indices, dtype=np.int64) - keep = np.zeros_like(local_halo["index"], dtype=bool) - matching_index = virgo.util.match.match(halo_indices, local_halo["index"]) - have_match = matching_index >= 0 - keep[matching_index[have_match]] = True - for name in local_halo: - local_halo[name] = local_halo[name][keep, ...] + nr_requested = halo_indices.shape[0] + have_match = np.zeros(nr_requested, dtype=np.int8) + if local_halo["index"].shape[0]: + keep = np.zeros_like(local_halo["index"], dtype=bool) + matching_index = virgo.util.match.match( + halo_indices, local_halo["index"] + ) + have_match[:] = matching_index >= 0 + keep[matching_index[have_match.astype(bool)]] = True + for name in local_halo: + local_halo[name] = local_halo[name][keep, ...] + + # Report any requested halos which are not in the catalogue. This + # is collective, so it must be done on all ranks. + comm.Allreduce(MPI.IN_PLACE, have_match, op=MPI.MAX) + have_match = have_match.astype(bool) + nr_matched = np.sum(have_match) + if comm_rank == 0: + print(f"Matched {nr_matched} of {nr_requested} requested halo indices") + if nr_matched < nr_requested: + # Written to the directory SOAP is being run from, and + # named after the output catalogue + output_file = sub_snapnum(args.output_file, args.snapshot_nr) + basename = os.path.basename(output_file) + basename = os.path.splitext(basename)[0] + filename = f"{basename}_unmatched_halo_indices.txt" + np.savetxt(filename, halo_indices[~have_match], fmt="%d") + print(f"WARNING: wrote unmatched halo indices to {filename}") # Discard satellites, if necessary if args.centrals_only: @@ -144,7 +166,7 @@ def __init__( # Exit if we don't have any halos if (total_nr_halos == 0) and (comm_rank == 0): - print("No halos found, aborting run") + print("No halos found, aborting run", flush=True) comm.Abort(1) # Assign halos to chunk tasks: diff --git a/SOAP/core/soap_args.py b/SOAP/core/soap_args.py index a8274831..0b50ed04 100644 --- a/SOAP/core/soap_args.py +++ b/SOAP/core/soap_args.py @@ -4,13 +4,60 @@ import os import subprocess import sys +import warnings from mpi4py import MPI +import numpy as np from virgo.mpi.util import MPIArgumentParser +from virgo.util.partial_formatter import PartialFormatter from . import combine_args +def get_halo_indices(parameters): + """ + Return the indices of the halos to process, or None if all halos should + be processed. The indices are either passed on the command line, or are + read from a file which contains one index per line. Blank lines and lines + starting with "#" are ignored. + + Returns a sorted array of the unique indices which were requested. + """ + + filename = parameters["halo_indices_file"] + if filename is not None: + # Substitute the snapshot number into the filename + pf = PartialFormatter() + filename = pf.format(filename, snap_nr=parameters["snap_nr"], file_nr=None) + if not os.path.exists(filename): + raise ValueError(f"Unable to find halo indices file: {filename}") + try: + with warnings.catch_warnings(): + # Empty files generate a warning, but we handle them below + warnings.filterwarnings( + "ignore", message="loadtxt: input contained no data" + ) + # converters=int prevents non-integer values from being silently + # truncated, which np.loadtxt would otherwise do + halo_indices = np.loadtxt( + filename, dtype=np.int64, ndmin=1, comments="#", converters=int + ) + except Exception as e: + raise ValueError(f"Unable to read halo indices file {filename}: {e}") + if halo_indices.ndim != 1: + raise ValueError( + f"Halo indices file must contain one index per line: {filename}" + ) + if halo_indices.shape[0] == 0: + raise ValueError(f"Halo indices file is empty: {filename}") + elif parameters["halo_indices"] is not None: + halo_indices = np.asarray(parameters["halo_indices"], dtype=np.int64) + else: + return None + + return np.unique(halo_indices) + + def get_git_hash() -> str: try: return ( @@ -46,13 +93,18 @@ def get_soap_args(comm): metavar="N", type=int, default=1, - help="Splits volume into N chunks and each compute node processes one chunk at a time", + help="Splits volume into N chunks and each compute node processes one chunk " + "at a time. Should be at least the number of nodes (default: 1)", ) parser.add_argument( - "--dmo", action="store_true", help="Run in dark matter only mode" + "--dmo", + action="store_true", + help="Run in dark matter only mode, skipping any hydro-only properties", ) parser.add_argument( - "--centrals-only", action="store_true", help="Only process central halos" + "--centrals-only", + action="store_true", + help="Only process central halos, discarding satellites", ) parser.add_argument( "--record-halo-timings", @@ -62,7 +114,8 @@ def get_soap_args(comm): parser.add_argument( "--record-property-timings", action="store_true", - help="Record time taken to process each property", + help="Record time taken to process each property. This doubles the size of " + "the output catalogue", ) parser.add_argument( "--max-halos", @@ -71,15 +124,25 @@ def get_soap_args(comm): default=0, help="(For debugging) only process the first N halos in the catalogue", ) - parser.add_argument( + halo_index_group = parser.add_mutually_exclusive_group() + halo_index_group.add_argument( "--halo-indices", nargs="*", type=int, help="Only process the specified halo indices", ) + halo_index_group.add_argument( + "--halo-indices-file", + type=str, + help="Only process the halo indices listed in the specified file, which " + "must contain one index per line. The snapshot number is substituted " + "into the filename, e.g. halo_indices_{snap_nr:04d}.txt", + ) parser.add_argument( "--reference-snapshot", - help="Specify reference snapshot number containing all particle types", + help="Specify reference snapshot number containing all particle types. " + "Used to determine the datasets and units of any particle types which " + "are missing from the snapshot being processed", metavar="N", type=int, ) @@ -88,27 +151,48 @@ def get_soap_args(comm): metavar="LEVEL", type=int, default=0, - help="Run with profiling (0=off, 1=first MPI rank only, 2=all ranks)", + help="Run with profiling (0=off, 1=first MPI rank only, 2=all ranks) " + "(default: 0)", ) parser.add_argument( "--max-ranks-reading", type=int, default=32, - help="Number of ranks per node reading snapshot data", + help="Number of ranks per node reading snapshot data. Can be reduced to " + "avoid overloading the file system (default: 32)", ) parser.add_argument( "--output-parameters", type=str, default="", - help="Where to write the used parameters", + help="Where to write the parameters used by this run, in yaml format", + ) + parser.add_argument( + "--snipshot", + action="store_true", + help="Run in snipshot mode, overriding the value of SelectOutput in the " + "snapshot header", + ) + parser.add_argument( + "--snapshot", + action="store_true", + help="Run in snapshot mode, overriding the value of SelectOutput in the " + "snapshot header", ) - parser.add_argument("--snipshot", action="store_true", help="Run in snipshot mode") - parser.add_argument("--snapshot", action="store_true", help="Run in snapshot mode") all_args = parser.parse_args() # Combine with parameters from configuration file if comm.Get_rank() == 0: - all_args = combine_args.combine_arguments(all_args, all_args.config_file) + try: + all_args = combine_args.combine_arguments(all_args, all_args.config_file) + # Halo indices are read on this rank and broadcast as an array, + # since there can be a large number of them + all_args["Parameters"]["halo_indices"] = get_halo_indices( + all_args["Parameters"] + ) + except ValueError as e: + print(e, flush=True) + comm.Abort(1) all_args["git_hash"] = get_git_hash() else: all_args = None @@ -173,7 +257,7 @@ def get_soap_args(comm): while not os.path.exists(dirname): dirname = os.path.dirname(dirname) if not os.access(dirname, os.W_OK): - print("Can't write to output directory") + print("Can't write to output directory", flush=True) comm.Abort(1) # Check if the FOF files exist if args.fof_group_filename != "": @@ -181,7 +265,7 @@ def get_soap_args(comm): snap_nr=args.snapshot_nr, file_nr=0 ) if not os.path.exists(fof_filename): - print(f"Could not find FOF group catalogue: {fof_filename}") + print(f"Could not find FOF group catalogue: {fof_filename}", flush=True) comm.Abort(1) if args.fof_radius_filename != "": assert args.fof_group_filename != "" @@ -189,19 +273,21 @@ def get_soap_args(comm): snap_nr=args.snapshot_nr, file_nr=0 ) if not os.path.exists(fof_filename): - print(f"Could not find FOF radius catalogue: {fof_filename}") + print( + f"Could not find FOF radius catalogue: {fof_filename}", flush=True + ) comm.Abort(1) # This really should be done in parameter_file.py args.separate_chunks = args.calculations.get("separate_chunks", []) if not isinstance(args.separate_chunks, list): - print("Invalid form for separate_chunks") + print("Invalid form for separate_chunks", flush=True) comm.Abort(1) for threshold in args.separate_chunks: if ("n_bound_threshold" not in threshold) or ( "n_halo_per_chunk" not in threshold ): - print("Invalid form for separate_chunks") + print("Invalid form for separate_chunks", flush=True) comm.Abort(1) args.separate_chunks = sorted( args.separate_chunks, diff --git a/SOAP/core/swift_cells.py b/SOAP/core/swift_cells.py index 62794afa..53e25ee0 100644 --- a/SOAP/core/swift_cells.py +++ b/SOAP/core/swift_cells.py @@ -442,7 +442,10 @@ def verify_extra_input(self, comm): dset = list(extra_metadata[parttype].keys())[0] npart_extra = extra_file[f"{parttype}/{dset}"].shape[0] if npart_snapshot[parttype] != npart_extra: - print(f"Incorrect number of {parttype} in {extra_filename}") + print( + f"Incorrect number of {parttype} in {extra_filename}", + flush=True, + ) comm.Abort(1) def check_datasets_exist(self, required_datasets, halo_prop_list): diff --git a/SOAP/group_membership.py b/SOAP/group_membership.py index 35cc7e0d..0f03b0ff 100644 --- a/SOAP/group_membership.py +++ b/SOAP/group_membership.py @@ -210,7 +210,7 @@ def main(): try: os.makedirs(os.path.dirname(output_filename), exist_ok=True) except OSError as e: - print(f"Error creating output directory: {e}") + print(f"Error creating output directory: {e}", flush=True) comm.Abort(1) comm.barrier() diff --git a/misc/recalculate_xrays.py b/misc/recalculate_xrays.py index 0e3661e4..293cb9e5 100644 --- a/misc/recalculate_xrays.py +++ b/misc/recalculate_xrays.py @@ -218,7 +218,7 @@ def recalculate_xrays(snap_file, output_filename, units, xray_calculator): try: os.makedirs(os.path.dirname(output_filename), exist_ok=True) except OSError as e: - print(f"Error creating output directory: {e}") + print(f"Error creating output directory: {e}", flush=True) comm.Abort(1) comm.barrier() diff --git a/tests/COLIBRE/find_halo_ids.py b/tests/COLIBRE/find_halo_ids.py index 7d36bf3f..9af61980 100755 --- a/tests/COLIBRE/find_halo_ids.py +++ b/tests/COLIBRE/find_halo_ids.py @@ -17,9 +17,11 @@ def find_halo_indices(sim, snap_nr, boxsize): index = f["InputHalos/HaloCatalogueIndex"][()] is_central = f["InputHalos/IsCentral"][()] nstar = f['BoundSubhalo/NumberOfStarParticles'][:] + # Diagnostics are written to stderr, so that stdout can be redirected + # to a file containing only the halo indices if np.sum(is_central[mask]) == 0: - print('No centrals loaded') - print(f'Max number of stars: {np.max(nstar[mask])}') + print('No centrals loaded', file=sys.stderr) + print(f'Max number of stars: {np.max(nstar[mask])}', file=sys.stderr) return index[mask] @@ -29,5 +31,7 @@ def find_halo_indices(sim, snap_nr, boxsize): boxsize = float(sys.argv[3]) indices = find_halo_indices(sim, snap_nr, boxsize) - indices_list = " ".join([str(i) for i in indices]) - print(indices_list) + # Print one index per line, so that the output can be redirected to a + # file which can be passed to SOAP with --halo-indices-file + for i in indices: + print(i) diff --git a/tests/COLIBRE/halo_indices_0092.txt b/tests/COLIBRE/halo_indices_0092.txt new file mode 100644 index 00000000..f30a6b3e --- /dev/null +++ b/tests/COLIBRE/halo_indices_0092.txt @@ -0,0 +1,290 @@ +# Halo indices to do: all halos with x<5, y<5, and z<5 cMpc in snap 92 +# Generated with `python tests/COLIBRE/find_halo_ids.py L0025N0188/Thermal 92 5` +17079 +20065 +22326 +27035 +27037 +34951 +36275 +39305 +40463 +40495 +44619 +45938 +45939 +48451 +49639 +49646 +49657 +51919 +51938 +53031 +56226 +57389 +60474 +61533 +62532 +64282 +67777 +67801 +68350 +69437 +69897 +70019 +70400 +71432 +71975 +72461 +72943 +73459 +73932 +73939 +73962 +74036 +74440 +75916 +7 +819 +1684 +1689 +2295 +3231 +3232 +5123 +5928 +5954 +6828 +6853 +6863 +7859 +8930 +8932 +10191 +11482 +11496 +11499 +11507 +12942 +12951 +14444 +16066 +16096 +16097 +19060 +19077 +20095 +21208 +21212 +21218 +23496 +25868 +25870 +25872 +26584 +27042 +28291 +29610 +30883 +30942 +32234 +32268 +33601 +34967 +34986 +36282 +37688 +37689 +37699 +37710 +37711 +39117 +41890 +41892 +43234 +44574 +44576 +44589 +44598 +44599 +45947 +45948 +47198 +47200 +47208 +48449 +48461 +49654 +50787 +50788 +50798 +51921 +51967 +51974 +54152 +54699 +55224 +55227 +56283 +56284 +57383 +57385 +57386 +57392 +59473 +59476 +62536 +62546 +62547 +63591 +64274 +64917 +66673 +66675 +66697 +67779 +67781 +68342 +68911 +69426 +69430 +69432 +70953 +70955 +72942 +72948 +72950 +73452 +73931 +73933 +74460 +75423 +75424 +39078 +47222 +61547 +73961 +821 +6842 +12933 +14454 +14459 +16071 +27065 +28285 +29595 +29612 +33595 +36311 +39103 +43263 +45950 +45954 +53065 +54100 +56277 +60517 +61497 +62524 +62554 +63595 +63987 +64286 +64920 +64921 +65527 +66694 +67786 +70404 +71417 +71420 +71968 +73941 +74448 +74945 +75903 +3380 +3381 +5116 +5128 +5133 +9138 +14494 +14732 +16055 +16241 +18094 +19219 +19222 +24819 +24846 +28564 +28566 +28569 +31201 +32240 +33597 +34944 +37694 +43240 +45927 +47203 +48440 +48446 +49637 +50805 +53029 +53033 +54146 +55183 +55186 +55189 +55195 +55405 +56258 +60703 +61536 +61759 +63763 +64294 +64307 +65631 +66120 +66122 +67241 +68344 +68905 +69436 +69438 +70394 +70954 +70957 +71430 +72453 +72454 +72962 +72963 +73935 +73945 +73963 +74438 +74441 +75415 +75918 +5939 +18100 +22345 +32291 +55169 +55194 +71418 +72951 +30929 +57361 +62589 +67803 +36297 +50822 +67788 +68347 +69899 diff --git a/tests/COLIBRE/run_L0025N0188_Thermal.sh b/tests/COLIBRE/run_L0025N0188_Thermal.sh index 58969b78..f182df03 100755 --- a/tests/COLIBRE/run_L0025N0188_Thermal.sh +++ b/tests/COLIBRE/run_L0025N0188_Thermal.sh @@ -20,8 +20,9 @@ sim="L0025N0188/Thermal" # Snapshot number to do snapnum=0092 -# Halo indices to do: all halos with x<5, y<5, and z<5 cMpc in snap 92 -halo_indices="17079 20065 22326 27035 27037 34951 36275 39305 40463 40495 44619 45938 45939 48451 49639 49646 49657 51919 51938 53031 56226 57389 60474 61533 62532 64282 67777 67801 68350 69437 69897 70019 70400 71432 71975 72461 72943 73459 73932 73939 73962 74036 74440 75916 7 819 1684 1689 2295 3231 3232 5123 5928 5954 6828 6853 6863 7859 8930 8932 10191 11482 11496 11499 11507 12942 12951 14444 16066 16096 16097 19060 19077 20095 21208 21212 21218 23496 25868 25870 25872 26584 27042 28291 29610 30883 30942 32234 32268 33601 34967 34986 36282 37688 37689 37699 37710 37711 39117 41890 41892 43234 44574 44576 44589 44598 44599 45947 45948 47198 47200 47208 48449 48461 49654 50787 50788 50798 51921 51967 51974 54152 54699 55224 55227 56283 56284 57383 57385 57386 57392 59473 59476 62536 62546 62547 63591 64274 64917 66673 66675 66697 67779 67781 68342 68911 69426 69430 69432 70953 70955 72942 72948 72950 73452 73931 73933 74460 75423 75424 39078 47222 61547 73961 821 6842 12933 14454 14459 16071 27065 28285 29595 29612 33595 36311 39103 43263 45950 45954 53065 54100 56277 60517 61497 62524 62554 63595 63987 64286 64920 64921 65527 66694 67786 70404 71417 71420 71968 73941 74448 74945 75903 3380 3381 5116 5128 5133 9138 14494 14732 16055 16241 18094 19219 19222 24819 24846 28564 28566 28569 31201 32240 33597 34944 37694 43240 45927 47203 48440 48446 49637 50805 53029 53033 54146 55183 55186 55189 55195 55405 56258 60703 61536 61759 63763 64294 64307 65631 66120 66122 67241 68344 68905 69436 69438 70394 70954 70957 71430 72453 72454 72962 72963 73935 73945 73963 74438 74441 75415 75918 5939 18100 22345 32291 55169 55194 71418 72951 30929 57361 62589 67803 36297 50822 67788 68347 69899" +# File containing the halo indices to do, one index per line. The snapshot +# number is substituted into the filename. +halo_indices_file="tests/COLIBRE/halo_indices_{snap_nr:04d}.txt" # Create parameters files python tests/COLIBRE/create_parameters_file.py @@ -32,6 +33,6 @@ rm -r output/SOAP-tmp # Run SOAP on eight cores processing the selected halos. Use 'python3 -m pdb' to start in the debugger. mpirun -np 8 python SOAP/compute_halo_properties.py \ ./tests/COLIBRE/test_parameters.yml \ - --halo-indices ${halo_indices} \ + --halo-indices-file ${halo_indices_file} \ --sim-name=${sim} --snap-nr=${snapnum} --chunks=1 diff --git a/tests/test_halo_indices_file.py b/tests/test_halo_indices_file.py new file mode 100644 index 00000000..9f580a2b --- /dev/null +++ b/tests/test_halo_indices_file.py @@ -0,0 +1,64 @@ +#!/bin/env python + +import numpy as np +import pytest + +from SOAP.core.soap_args import get_halo_indices + + +def make_parameters(tmp_path, contents, snap_nr=10): + """ + Write a halo indices file and return the parameters required by + get_halo_indices + """ + filename = tmp_path / f"halo_indices_{snap_nr:04d}.txt" + filename.write_text(contents) + return { + "halo_indices": None, + "halo_indices_file": str(tmp_path / "halo_indices_{snap_nr:04d}.txt"), + "snap_nr": snap_nr, + } + + +def test_no_halo_indices(): + parameters = {"halo_indices": None, "halo_indices_file": None, "snap_nr": 10} + assert get_halo_indices(parameters) is None + + +def test_halo_indices_command_line(): + parameters = { + "halo_indices": [3, 1, 2, 1], + "halo_indices_file": None, + "snap_nr": 10, + } + assert np.array_equal(get_halo_indices(parameters), [1, 2, 3]) + + +def test_halo_indices_file_single_index(tmp_path): + parameters = make_parameters(tmp_path, "7\n") + halo_indices = get_halo_indices(parameters) + assert halo_indices.shape == (1,) + assert halo_indices[0] == 7 + + +def test_halo_indices_file_comments_and_blank_lines(tmp_path): + parameters = make_parameters(tmp_path, "# a comment\n\n1\n2 # another comment\n\n") + assert np.array_equal(get_halo_indices(parameters), [1, 2]) + + +def test_halo_indices_file_empty(tmp_path): + parameters = make_parameters(tmp_path, "# no indices here\n") + with pytest.raises(ValueError, match="empty"): + get_halo_indices(parameters) + + +def test_halo_indices_file_multiple_columns(tmp_path): + parameters = make_parameters(tmp_path, "1 2\n3 4\n") + with pytest.raises(ValueError, match="one index per line"): + get_halo_indices(parameters) + + +def test_halo_indices_file_not_integer(tmp_path): + parameters = make_parameters(tmp_path, "1\n2.5\n") + with pytest.raises(ValueError, match="Unable to read"): + get_halo_indices(parameters) diff --git a/tests/test_read_subfind.py b/tests/test_read_subfind.py index 91c07de4..b876d0f7 100644 --- a/tests/test_read_subfind.py +++ b/tests/test_read_subfind.py @@ -75,7 +75,7 @@ def test_read_gadget4_groupnr(): if comm_rank == 0: print(f"Number of groups from fof_subhalo_tab = {nr_groups_from_subtab}") if nr_groups_from_subtab != nr_groups_from_grnr: - print("Number of groups does not agree!") + print("Number of groups does not agree!", flush=True) comm.Abort(1) # Ensure nbound arrays are partitioned the same way diff --git a/tests/test_shared_mesh.py b/tests/test_shared_mesh.py index 6a425385..ec3e25af 100644 --- a/tests/test_shared_mesh.py +++ b/tests/test_shared_mesh.py @@ -134,7 +134,10 @@ def periodic_distance_squared(pos, centre): if nr_failures == 0: print(f" OK") else: - print(f" {nr_failures} of {nr_queries*comm_size} queries FAILED") + print( + f" {nr_failures} of {nr_queries*comm_size} queries FAILED", + flush=True, + ) comm.Abort(1)