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)