diff --git a/SOAP/compute_halo_properties.py b/SOAP/compute_halo_properties.py index 1cf0279f..3649c702 100644 --- a/SOAP/compute_halo_properties.py +++ b/SOAP/compute_halo_properties.py @@ -164,6 +164,9 @@ def compute_halo_properties(): cellgrid.snapshot_datasets.setup_defined_constants( parameter_file.get_defined_constants() ) + # Tell the parameter file which datasets are in the input files, so that + # properties which cannot be computed can be skipped or reported + parameter_file.set_available_datasets(cellgrid.snapshot_datasets.datasets_in_file) parameter_file.record_property_timings = args.record_property_timings # Try to load parameters for RecentlyHeatedGasFilter. If a property that uses the @@ -473,6 +476,7 @@ def compute_halo_properties(): if args.record_property_timings: print("Storing processing time for each property") parameter_file.print_unregistered_properties(halo_prop_list, dmo=args.dmo) + parameter_file.print_skipped_properties(halo_prop_list, dmo=args.dmo) parameter_file.print_invalid_properties(halo_prop_list) parameter_file.print_variation_warnings() if not parameter_file.renclose_enabled(): @@ -481,6 +485,12 @@ def compute_halo_properties(): ) category_filter.print_filters() + # Properties enabled in the parameter file must be computed, so abort + # if the input files do not contain the datasets they require + parameter_file.print_uncomputable_properties() + if len(parameter_file.uncomputable_properties): + comm_world.Abort(1) + # Ensure output dir exists if comm_world_rank == 0: try: diff --git a/SOAP/core/parameter_file.py b/SOAP/core/parameter_file.py index 59d62b5a..4e8d392a 100644 --- a/SOAP/core/parameter_file.py +++ b/SOAP/core/parameter_file.py @@ -15,6 +15,28 @@ from SOAP import property_table +# Lazily built map from the output name of a property to its entry in the +# property table. The property table itself is keyed on an internal name. +_PROPERTY_BY_NAME = None + + +def _property_by_name(name: str): + """ + Look up a property in the property table by its output name. + + Returns None if there is no property with this name, which happens for + properties that are not defined by the property table (e.g. the dummy + properties used when generating the documentation). + """ + global _PROPERTY_BY_NAME + if _PROPERTY_BY_NAME is None: + _PROPERTY_BY_NAME = { + prop.name: prop + for prop in property_table.PropertyTable.full_property_list.values() + } + return _PROPERTY_BY_NAME.get(name) + + # Known parameter file structure, used by check_schema to flag typos. A value # of None means the keys directly under that section are user-named or free-form # and are not checked; a set lists the only keys allowed directly under that @@ -104,6 +126,24 @@ def __init__( # on a single rank via print_variation_warnings() self.variation_warnings = [] + # Datasets present in the input files, as {particle type: set of dataset + # names}. While this is None no property is considered uncomputable, + # which is the case when building a parameter file from a dictionary. + self.available_datasets = None + + # Properties which are not in the parameter file, and which are not + # calculated because the input files lack the datasets they need. + # Only used when calculate_missing_properties is True. + self.skipped_properties = set() + + # Properties which are enabled in the parameter file, but which cannot + # be calculated because the input files lack the datasets they need. + self.uncomputable_properties = {} + + # Names of the properties used by the filters defined in the parameter + # file, generated on demand by _filter_property_names() + self.filter_property_names = None + def get_parameters(self) -> Dict: """ Get a copy of the parameter dictionary. @@ -148,6 +188,63 @@ def _validate_filter_name(self, filter_name, context: str) -> None: f'the "filters" section of the parameter file' ) + def set_available_datasets(self, datasets_in_file: Dict) -> None: + """ + Record which datasets are present in the input files + + Parameters: + - datasets_in_file: Dict + Dictionary of the datasets present in the snapshot and extra-input + files, as {particle type: set of dataset names}. + """ + self.available_datasets = datasets_in_file + + def missing_datasets(self, property_name: str) -> List[str]: + """ + Get the datasets which are required to compute the given property, but + which are not present in the input files. + + Particle types which are absent from the input files entirely are not + considered. This matches the check done in SWIFTCellGrid.check_datasets_exist. + + Returns the (aliased) names of the missing datasets. An empty list is + returned if the property can be computed. + """ + if self.available_datasets is None: + return [] + prop = _property_by_name(property_name) + if prop is None: + return [] + missing = [] + for dataset in prop.particle_properties: + ptype, name = self.get_particle_property(dataset) + # Skip particle types which are not in the input files at all + if ptype not in self.available_datasets: + continue + dataset_name = f"{ptype}/{name}" + if (name not in self.available_datasets[ptype]) and ( + dataset_name not in missing + ): + missing.append(dataset_name) + return missing + + def _filter_property_names(self) -> set: + """ + Get the names of the properties used by the filters defined in the + parameter file. + + These properties are never skipped, since the filters that use them + would otherwise be undefined. + """ + if self.filter_property_names is None: + self.filter_property_names = set() + for filter_info in self.get_filters().values(): + for prop in filter_info.get("properties", []): + # Filters name properties by their full path in the output, + # e.g. BoundSubhalo/NumberOfGasParticles + self.filter_property_names.add(prop.split("/")[-1]) + return self.filter_property_names + def get_property_filters(self, base_halo_type: str, full_list: List[str]) -> Dict: """ Get a dictionary with the filter that should be applied to each @@ -176,18 +273,18 @@ def get_property_filters(self, base_halo_type: str, full_list: List[str]) -> Dic if not base_halo_type in self.parameters: self.parameters[base_halo_type] = {} - # Handle the case where no properties are listed for the halo type + # Handle the case where no properties are listed for the halo type. Each + # property is then treated as if it were missing from the parameter file. if not "properties" in self.parameters[base_halo_type]: self.parameters[base_halo_type]["properties"] = {} - for property in full_list: - self.parameters[base_halo_type]["properties"][ - property - ] = self.calculate_missing_properties() + listed = self.parameters[base_halo_type]["properties"] filters = {} for property in full_list: + # Datasets this property needs which are not in the input files + missing = self.missing_datasets(property) # Check if property is listed in the parameter file for this base_halo_type - if property in self.parameters[base_halo_type]["properties"]: - filter_name = self.parameters[base_halo_type]["properties"][property] + if property in listed: + filter_name = listed[property] # filter_name will a dict if we want different behaviour # for snapshots/snipshots if isinstance(filter_name, dict): @@ -200,15 +297,27 @@ def get_property_filters(self, base_halo_type: str, full_list: List[str]) -> Dic if filter_name == True: filter_name = "basic" filters[property] = filter_name + # An uncomputable property enabled in the parameter file was asked for + # explicitly, so we abort rather than quietly omitting it + if filter_name and missing: + self.uncomputable_properties[property] = missing # Property is not listed in the parameter file for this base_halo_type + elif not self.calculate_missing_properties(): + filters[property] = False + elif missing and property not in self._filter_property_names(): + # The property was not asked for explicitly and cannot be + # computed, so it is skipped. Properties used by a filter are + # never skipped, since the filter would then be undefined. + filters[property] = False + listed[property] = False + self.skipped_properties.add(property) else: - if self.calculate_missing_properties(): - filters[property] = "basic" - self.parameters[base_halo_type]["properties"][property] = "basic" - if self.unregistered_parameters is not None: - self.unregistered_parameters.add((base_halo_type, property)) - else: - filters[property] = False + filters[property] = "basic" + listed[property] = "basic" + if self.unregistered_parameters is not None: + self.unregistered_parameters.add((base_halo_type, property)) + if missing: + self.uncomputable_properties[property] = missing if isinstance(filters[property], str): self._validate_filter_name( filters[property], f"{base_halo_type}/{property}" @@ -241,12 +350,9 @@ def print_unregistered_properties( # In a DMO run, drop properties that will be skipped because they are # not DMO properties, so the printed list matches the output if dmo and halo_prop_list is not None: - dmo_flag = {} - for halo_type in halo_prop_list: - for prop in halo_type.property_list.values(): - dmo_flag[(halo_type.base_halo_type, prop.name)] = prop.dmo_property + non_dmo_names = self._non_dmo_property_names(halo_prop_list) unregistered = { - entry for entry in unregistered if dmo_flag.get(entry, True) + entry for entry in unregistered if entry[1] not in non_dmo_names } if len(unregistered): @@ -256,6 +362,72 @@ def print_unregistered_properties( for base_halo_type, property in sorted(unregistered): print(f" {base_halo_type.ljust(30)}{property}") + def _non_dmo_property_names(self, halo_prop_list) -> set: + """ + Get the names of the properties which are not calculated in a DMO run. + + Properties which are not found in halo_prop_list are not included, since + we cannot tell whether they would be calculated. + + Parameters: + - halo_prop_list: List + List of the halo property calculations that are enabled. + """ + names = set() + for halo_type in halo_prop_list: + for prop in halo_type.property_list.values(): + if not prop.dmo_property: + names.add(prop.name) + return names + + def print_skipped_properties(self, halo_prop_list=None, dmo: bool = False) -> None: + """ + Print a list of the properties which are not in the parameter file, and + which are not calculated because the input files lack the datasets they + need. + + Each property is listed once, rather than once per halo type, since + whether a dataset is present does not depend on the halo type. + + In a DMO run the property calculators skip any property that is not + flagged as a DMO property. + """ + + # A property that is explicitly enabled in the parameter file + # is reported as uncomputable instead + skipped = set(self.skipped_properties) - set(self.uncomputable_properties) + + # In a DMO run, drop properties that would be skipped anyway because + # they are not DMO properties + if dmo and halo_prop_list is not None: + skipped -= self._non_dmo_property_names(halo_prop_list) + + if len(skipped): + print( + "Not computing the following properties as required datasets are missing:" + ) + for property in sorted(skipped): + print(f" {property}") + + def print_uncomputable_properties(self) -> None: + """ + Print a list of the properties which are enabled in the parameter file, + but which cannot be calculated because the input files lack the datasets + they need. Only the missing datasets are listed for each property, not + all of the datasets it requires. + """ + if not len(self.uncomputable_properties): + return + print( + f"Cannot compute {len(self.uncomputable_properties)} properties enabled in " + f"the parameter file (only the missing datasets are listed for each one):", + flush=True, + ) + for property in sorted(self.uncomputable_properties): + print(f" {property}", flush=True) + for dataset in self.uncomputable_properties[property]: + print(f" {dataset}", flush=True) + def print_invalid_properties(self, halo_prop_list) -> None: """ Print a list of any properties in the parameter file that are not present in diff --git a/SOAP/core/swift_cells.py b/SOAP/core/swift_cells.py index 53e25ee0..bb3fbcbc 100644 --- a/SOAP/core/swift_cells.py +++ b/SOAP/core/swift_cells.py @@ -451,9 +451,15 @@ def verify_extra_input(self, comm): def check_datasets_exist(self, required_datasets, halo_prop_list): # Check we have all the fields needed for each property # Doing it at this point rather than in masked cells since we want - # to output a list of properties that require the missing fields - for ptype in set(self.ptypes).intersection(set(required_datasets.keys())): - for name in required_datasets[ptype]: + # to output a list of properties that require the missing fields. + # Properties which cannot be computed are normally skipped or reported + # when the parameter file is resolved, so reaching this point indicates + # a problem with an alias or with a declared dependency. + missing_datasets = [] + for ptype in sorted( + set(self.ptypes).intersection(set(required_datasets.keys())) + ): + for name in sorted(required_datasets[ptype]): # Note that the field names in required_datasets have already had # any aliases applied, so we can check the raw files themselves in_extra = (self.extra_filenames is not None) and ( @@ -461,21 +467,34 @@ def check_datasets_exist(self, required_datasets, halo_prop_list): ) in_snap = name in self.snap_metadata[ptype] if not (in_extra or in_snap): - dataset = f"{ptype}/{name}" - print(f"The following properties require {dataset}:") - full_property_list = property_table.PropertyTable.full_property_list - for k, v in full_property_list.items(): - # Skip property if it doesn't require this dataset - if dataset not in v.particle_properties: - continue - # Only print if the property is being calculated for some halo type - for halo_prop in halo_prop_list: - if halo_prop.property_filters.get(v.name, False): - print(f" {v.name}") - break - raise KeyError( - f"Can't find required dataset {dataset} in input file(s)!" - ) + missing_datasets.append(f"{ptype}/{name}") + + if not missing_datasets: + return + + # Report the missing datasets + full_property_list = property_table.PropertyTable.full_property_list + for dataset in missing_datasets: + print(f"The following properties require {dataset}:") + n_properties = 0 + for k, v in full_property_list.items(): + # Skip property if it doesn't require this dataset + if dataset not in v.particle_properties: + continue + # Only print if the property is being calculated for some halo type + for halo_prop in halo_prop_list: + if halo_prop.property_filters.get(v.name, False): + print(f" {v.name}") + n_properties += 1 + break + if n_properties == 0: + # No enabled property lists this dataset, so it is one that SOAP + # always reads. Disabling properties will not help here. + print(f" (none, {dataset} is required for every calculation)") + raise KeyError( + "Can't find required dataset(s) " + f"{', '.join(missing_datasets)} in input file(s)!" + ) def prepare_read(self, ptype, mask): """ diff --git a/parameter_files/README.md b/parameter_files/README.md index d351abf6..d13a4363 100644 --- a/parameter_files/README.md +++ b/parameter_files/README.md @@ -245,8 +245,8 @@ defined_constants: Contains information about how to run SOAP -- **min_read_radius_cmpc**: Optional. Using the input halo catalogues SOAP makes an initial guess of the radius around each halo to read in. This value can be set so SOAP will read a minimum radius by default, which can be useful if large SOs are being calculated. -- **calculate_missing_properties**: Optional, default True. If set to true then SOAP will calculate any properties which are not listed in the parameter file. If set to false then SOAP will ignore these properties +- **calculate_missing_properties**: Optional, default True. If set to true then SOAP will calculate any properties which are not listed in the parameter file, provided the input files contain the datasets those properties require. Properties which cannot be calculated are skipped, and are listed at the start of the run. If set to false then SOAP will ignore any property which is not listed in the parameter file. + - **reduced_snapshots**: Optional. We create reduced snapshots where we keep the particles within the virial radius of certain objects. The values here determine which halos to keep. - **min_halo_mass**: The minimumum M200 halo mass to keep - **halo_bin_size_dex**: The size of the halo mass bins @@ -258,6 +258,7 @@ 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. +- **min_read_radius_cmpc**: Optional. Using the input halo catalogues SOAP makes an initial guess of the radius around each halo to read in. This value can be set so SOAP will read a minimum radius by default, which can be useful if large SOs are being calculated. - **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: diff --git a/tests/test_parameter_file.py b/tests/test_parameter_file.py index f6c2bdef..ef7a15cf 100644 --- a/tests/test_parameter_file.py +++ b/tests/test_parameter_file.py @@ -16,7 +16,11 @@ def make_parameter_file( - section=None, calculate_missing_properties=True, snipshot=False, filters=None + section=None, + calculate_missing_properties=True, + snipshot=False, + filters=None, + available_datasets=None, ): """ Build a ParameterFile with a single ApertureProperties section. @@ -28,7 +32,10 @@ def make_parameter_file( parameters["ApertureProperties"] = section if filters is not None: parameters["filters"] = filters - return ParameterFile(parameter_dictionary=parameters, snipshot=snipshot) + pf = ParameterFile(parameter_dictionary=parameters, snipshot=snipshot) + if available_datasets is not None: + pf.set_available_datasets(available_datasets) + return pf VARIATIONS = {"exclusive_50_kpc": {"radius_in_kpc": 50.0, "inclusive": False}} @@ -323,3 +330,51 @@ def test_check_schema_reports_all_errors_together(): message = str(excinfo.value) assert "Nonsense" in message assert "oops" in message + + +def test_auto_enabled_property_is_skipped_when_datasets_missing(): + # GasMass requires PartType0/Masses, DustMass additionally requires + # PartType0/TotalDustMassFractions, so only DustMass cannot be computed + pf = make_parameter_file( + section={"properties": {}, "variations": VARIATIONS}, + available_datasets={"PartType0": {"Masses"}}, + ) + filters = pf.get_property_filters("ApertureProperties", ["GasMass", "DustMass"]) + + assert filters["GasMass"] == "basic" + assert filters["DustMass"] == False + assert pf.skipped_properties == {"DustMass"} + assert pf.uncomputable_properties == {} + # The property must be recorded as disabled, so that the used parameters + # file does not claim SOAP computed something it did not + assert pf.parameters["ApertureProperties"]["properties"]["DustMass"] == False + + +def test_enabled_property_with_missing_datasets_is_uncomputable(): + # A property asked for by name is never quietly skipped + pf = make_parameter_file( + section={"properties": {"DustMass": True}, "variations": VARIATIONS}, + available_datasets={"PartType0": {"Masses"}}, + ) + filters = pf.get_property_filters("ApertureProperties", ["GasMass", "DustMass"]) + + assert filters["DustMass"] == "basic" + assert pf.skipped_properties == set() + # Only the missing datasets are reported, not every dataset it requires + assert pf.uncomputable_properties == { + "DustMass": ["PartType0/TotalDustMassFractions"] + } + + +def test_absent_particle_type_is_not_treated_as_missing(): + # A particle type which is absent entirely (e.g. gas in a DMO run) is + # handled elsewhere, so it must not make properties look uncomputable + pf = make_parameter_file( + section={"properties": {}, "variations": VARIATIONS}, + available_datasets={"PartType1": {"Masses"}}, + ) + filters = pf.get_property_filters("ApertureProperties", ["GasMass", "DustMass"]) + + assert filters["DustMass"] == "basic" + assert pf.skipped_properties == set() + assert pf.uncomputable_properties == {}