Merge pull request #13 from SWIFTSIM/support_new_soap_format

Update for new SOAP catalogue format

README.md

@@ -35,7 +35,8 @@ Usage:
 ./morphology-pipeline \
   -C/--config      <configuration directory> \
   -i/--input       <input folder that contains the snapshots and catalogues> \
-  -s/--snapshots   <input snapshot> \
+  -s/--snapshots   <input swift snapshot containing particle properties> \
+  -g/--groups      <input swift snapshot containing particle membership> \
   -c/--catalogues  <input catalogue> \
   -o/--output      <directory where output images and web pages are stored> \
  [-n/--run-names   <name for the run>] \
@@ -44,7 +45,9 @@ Usage:
  [-d/--debug] \
  [-l/--lazy]
 ```
-(the last four arguments are optional). `--debug` outputs additional debugging information, while
+(the last four arguments are optional). If you have a virtual snapshot which also contains
+membership information than `-s` and `-g` can be pointed at the same file.
+`--debug` outputs additional debugging information, while
 `--lazy` simply reads an existing meta-data file and creates the plots and web pages, without
 recomputing any quantities. This is the same behaviour as the comparison mode, but for a single run.
 

morphology-pipeline

@@ -11,8 +11,8 @@ Usage:
   ./morphology-pipeline \
     -C/--config      <configuration directory> \
     -i/--input       <input folder that contains the snapshots and catalogues> \
-    -s/--snapshots   <input snapshot> \
-    -g/--groups      <input particle membership> \
+    -s/--snapshots   <input swift snapshot containing particle properties> \
+    -g/--groups      <input swift snapshot containing particle membership> \
     -c/--catalogues  <input catalogue> \
     -o/--output      <directory where output images and web pages are stored> \
    [-n/--run-names   <name for the run>] \
@@ -38,17 +38,17 @@ parser = ap.ArgumentParser(
     ),
     epilog=(
         "Example usage:\n"
-        "  ./morphology-pipeline \ \n"
-        "    -C/--config      <configuration directory> \ \n"
-        "    -i/--input       <input folder that contains the snapshots and catalogues> \ \n"
-        "    -s/--snapshots   <input snapshot> \ \n"
-        "    -g/--groups      <input particle membership> \ \n"
-        "    -c/--catalogues  <input catalogue> \ \n"
-        "    -o/--output      <directory where output images and web pages are stored> \ \n"
-        "   [-n/--run-names   <name for the run>] \ \n"
-        "   [-j/--num-of-cpus <number of parallel processes to use>] \ \n"
-        "   [-m/--metadata    <prefix for the meta-data output file] \ \n"
-        "   [-d/--debug] \ \n"
+        "  ./morphology-pipeline \\ \n"
+        "    -C/--config      <configuration directory> \\ \n"
+        "    -i/--input       <input folder that contains the snapshots and catalogues> \\ \n"
+        "    -s/--snapshots   <input swift snapshot containing particle properties> \\ \n"
+        "    -g/--groups      <input swift snapshot containing particle membership> \\ \n"
+        "    -c/--catalogues  <input catalogue> \\ \n"
+        "    -o/--output      <directory where output images and web pages are stored> \\ \n"
+        "   [-n/--run-names   <name for the run>] \\ \n"
+        "   [-j/--num-of-cpus <number of parallel processes to use>] \\ \n"
+        "   [-m/--metadata    <prefix for the meta-data output file] \\ \n"
+        "   [-d/--debug] \\ \n"
         "   [-l/--lazy]"
     ),
     # we need this to preserve the line breaks in the epilog
@@ -376,7 +376,7 @@ if __name__ == "__main__":
 
             # make sure galaxies with plots are processed first, since they take
             # longer than average
-            arglist = sorted(arglist, key=lambda x: x[7], reverse=True)
+            arglist = sorted(arglist, key=lambda x: x[11], reverse=True)
 
             # determine the appropriate number of parallel processes to use:
             #  - the number requested by the user

morpholopy/filtered_catalogue.py

@@ -42,9 +42,8 @@ def __init__(
         """
         Constructor.
 
-        Reads the catalogue and only keeps galaxies with
-        structure type 10 (centrals) that have a stellar and gas
-        mass within the provided limits.
+        Reads the catalogue and only keeps central galaxies with
+        a stellar and gas mass within the provided limits.
 
         Parameters:
          - catalogue: VelociraptorCatalogue
@@ -79,7 +78,7 @@ def __init__(
         # compute the mask
         mask = (
             (Mstar >= minimum_mass_stars)
-            & (catalogue.get_quantity("structure_type.structuretype") == 10)
+            & (catalogue.get_quantity("structure_type.structuretype") == 1)
             & (Mgas > minimum_mass_gas)
         )
         # turn the mask into a list of indices

morpholopy/galaxy_data.py

@@ -435,7 +435,6 @@ def process_galaxy(
     # (re)load the catalogue to read the properties of this particular galaxy
     # we need to disregard the units to avoid problems with the SFR unit
     catalogue = load_catalogue(catalogue_filename, disregard_units=True)
-    membership_data = h5.File(halo_membership_filename, "r")
 
     swift_mask = swiftsimio.mask(snapshot_filename, spatial_only=True)
     # SWIFT data is stored in comoving units, so we need to un-correct
@@ -450,27 +449,33 @@ def process_galaxy(
         / length_factor
     )
 
-    spatial_mask = [
+    load_region = [
         [halo_x - r_size, halo_x + r_size],
         [halo_y - r_size, halo_y + r_size],
         [halo_z - r_size, halo_z + r_size],
     ]
+    swift_mask.constrain_spatial(load_region)
 
-    swift_mask.constrain_spatial(spatial_mask)
-    data = swiftsimio.load(snapshot_filename)
+    membership_data = swiftsimio.load(halo_membership_filename, mask=swift_mask)
+    data = swiftsimio.load(snapshot_filename, mask=swift_mask)
+
+    # The value group_nr_bound in the membership files corresponds to the index
+    # of the object in the original subhalo catalogue. In general this is not
+    # equal to the index of the object in the final SOAP catalogue.
+    halo_catalogue_index = catalogue.get_quantity("inputhalos.halocatalogueindex")[
+        galaxy_index
+    ]
 
     particle_name_masks = {}
-    for particle_name, particle_type in zip(
-        data.metadata.present_particle_names, data.metadata.present_particle_types
-    ):
+    for particle_name in ["gas", "dark_matter", "stars", "black_holes"]:
         mask_per_type = (
-            membership_data[f"PartType{particle_type}"]["GroupNr_bound"][:]
-            == galaxy_index
+            getattr(getattr(membership_data, particle_name), "group_nr_bound").value
+            == halo_catalogue_index
         )
         particle_name_masks[particle_name] = mask_per_type
-    MaskTuple = namedtuple("MaskCollection", data.metadata.present_particle_names)
+    MaskTuple = namedtuple("MaskCollection", data.metadata.present_group_names)
     mask = MaskTuple(**particle_name_masks)
-    membership_data.close()
+    del membership_data
 
     # mask out all particles that are actually bound to the galaxy
     # while at it: convert everything to physical coordinates

morpholopy/galaxy_plots.py

@@ -16,7 +16,7 @@
 from swiftsimio.visualisation.projection import project_gas
 from swiftsimio.visualisation.projection import project_pixel_grid
 from swiftsimio.visualisation.slice import kernel_gamma
-from swiftsimio.visualisation.smoothing_length_generation import (
+from swiftsimio.visualisation.smoothing_length.generate import (
     generate_smoothing_lengths,
 )
 from astropy.visualization import make_lupton_rgb