diff --git a/docs/changes/2188.feature.md b/docs/changes/2188.feature.md new file mode 100644 index 0000000000..c33cd13fcf --- /dev/null +++ b/docs/changes/2188.feature.md @@ -0,0 +1 @@ +Allow NSHOW to follow a pre-defined power law for simulation production grid definition. diff --git a/docs/source/api-reference/production_configuration.md b/docs/source/api-reference/production_configuration.md index a1f7aef886..acb6b17b7f 100644 --- a/docs/source/api-reference/production_configuration.md +++ b/docs/source/api-reference/production_configuration.md @@ -52,6 +52,17 @@ the calculation of the number of events to be simulated given a pre-determined m :members: ``` + +(job-spec-builder)= + +## job_spec_builder + +```{eval-rst} +.. automodule:: production_configuration.job_spec_builder + :members: +``` + + (interpolation-handler)= ## interpolation_handler diff --git a/src/simtools/applications/simulate_prod_htcondor_generator.py b/src/simtools/applications/simulate_prod_htcondor_generator.py index e4dc97f76b..8933b79a97 100644 --- a/src/simtools/applications/simulate_prod_htcondor_generator.py +++ b/src/simtools/applications/simulate_prod_htcondor_generator.py @@ -51,6 +51,11 @@ Base directory for simulation output packages passed through as ``pack_for_grid_register``. priority (int, optional) Job priority (default: 1). +nshow_power_index (float, optional) + Power-law index used to scale the baseline ``nshow`` for each ``energy_range`` entry. +nshow_reference_energy (quantity, optional) + Reference energy used for ``nshow`` power-law scaling (e.g. ``"100 GeV"``). + Required when ``nshow_power_index`` is configured. (all other command line arguments are identical to those of :ref:`simulate_prod`). @@ -97,6 +102,26 @@ def _add_arguments(parser): required=False, default=None, ) + parser.add_argument( + "--nshow_power_index", + help=( + "Power-law index used to scale the baseline nshow with the geometric-mean energy " + "of each energy_range entry." + ), + type=float, + required=False, + default=None, + ) + parser.add_argument( + "--nshow_reference_energy", + help=( + "Reference energy for nshow power-law scaling (for example: '100 GeV'). " + "Required together with --nshow_power_index." + ), + type=str, + required=False, + default=None, + ) def main(): diff --git a/src/simtools/job_execution/htcondor_script_generator.py b/src/simtools/job_execution/htcondor_script_generator.py index 705ee7868d..164a0ea74c 100644 --- a/src/simtools/job_execution/htcondor_script_generator.py +++ b/src/simtools/job_execution/htcondor_script_generator.py @@ -8,33 +8,19 @@ """ -import ast -import itertools import logging import re from pathlib import Path import astropy.units as u -import simtools.version as simtools_version -from simtools.configuration import defaults +from simtools.production_configuration.job_spec_builder import ( + build_job_specs, + resolve_array_layout_name, +) _logger = logging.getLogger(__name__) -_GRID_AXES = [ - "primary", - "azimuth_angle", - "zenith_angle", - "model_version", - "corsika_le_interaction", - "corsika_he_interaction", -] - -_GRID_AXIS_DEFAULTS = { - "corsika_le_interaction": defaults.CORSIKA_LE_INTERACTION, - "corsika_he_interaction": defaults.CORSIKA_HE_INTERACTION, -} - _PARAMS_FIELDS = [ "apptainer_label", "primary", @@ -56,45 +42,6 @@ ] -def _normalize_to_list(value): - """Normalize scalar values to lists of length one.""" - if isinstance(value, list): - return value - if isinstance(value, tuple): - return list(value) - return [value] - - -def _normalize_grid_axes(args_dict): - """Return normalized grid axes for cartesian product expansion.""" - return { - axis: ( - _normalize_to_list(args_dict[axis]) - if axis in args_dict and args_dict[axis] is not None - else [_GRID_AXIS_DEFAULTS[axis]] - if axis in _GRID_AXIS_DEFAULTS - else [None] - ) - for axis in _GRID_AXES - } - - -def _normalize_energy_ranges(energy_range): - """Normalize energy range argument to a list of (e_min, e_max) pairs.""" - if isinstance(energy_range, tuple) and len(energy_range) == 2: - return [energy_range] - - if isinstance(energy_range, list): - if len(energy_range) == 2 and all(hasattr(item, "to") for item in energy_range): - return [(energy_range[0], energy_range[1])] - if all(isinstance(item, (list, tuple)) and len(item) == 2 for item in energy_range): - return [tuple(item) for item in energy_range] - - raise ValueError( - "energy_range must be one pair (e_min, e_max) or a list of (e_min, e_max) pairs." - ) - - def _resolve_apptainer_images(apptainer_image_arg): """ Resolve and validate apptainer image configuration. @@ -179,141 +126,11 @@ def _sanitize_label_for_params(label): return re.sub(r"\s+", "_", str(label).strip()) -def _resolve_array_layout_name(array_layout_name, model_version): - """Resolve array layout configuration for a specific model version.""" - if isinstance(array_layout_name, list) and len(array_layout_name) == 1: - array_layout_name = array_layout_name[0] - - # Configurator/_fill_config stringifies dict values when rebuilding argparse arguments. - if isinstance(array_layout_name, str) and array_layout_name.strip().startswith("{"): - try: - parsed_layout = ast.literal_eval(array_layout_name) - if isinstance(parsed_layout, dict): - array_layout_name = parsed_layout - except (SyntaxError, ValueError): - return array_layout_name - - if not isinstance(array_layout_name, dict) or list(array_layout_name) != ["by_version"]: - return array_layout_name - - resolved = simtools_version.resolve_by_version( - {"array_layout_name": array_layout_name}, model_version - ) - return resolved["array_layout_name"] - - -def _get_energy_range_for_zenith_angle(zenith_angle, energy_range_pair, corsika_limits): - """ - Return a zenith-dependent energy range pair or None to skip the simulation step. - - Dummy function - to be implemented. - """ - _ = (zenith_angle, corsika_limits) - return energy_range_pair - - -def _get_core_scatter_max_for_zenith_angle(zenith_angle, core_scatter, corsika_limits): - """Return zenith-dependent max core-scatter value. - - Dummy function - to be implemented. - """ - _ = (zenith_angle, corsika_limits) - return core_scatter[1] - - -def _get_nshow_for_energy_range_and_zenith_angle( - energy_range_pair, zenith_angle, nshow, corsika_limits -): - """Return nshow that may depend on energy range and zenith angle. - - Dummy function - to be implemented. - """ - _ = (energy_range_pair, zenith_angle, corsika_limits) - return nshow - - -def _build_job_specs(args_dict, image_labels): - """Build backend-agnostic job specs from comparison and production grids.""" - grid_axes = _normalize_grid_axes(args_dict) - energy_ranges = _normalize_energy_ranges(args_dict["energy_range"]) - base_pack_dir = args_dict.get("simulation_output") or "simtools-output" - corsika_limits = args_dict.get("corsika_limits") - core_scatter = args_dict["core_scatter"] - nshow = args_dict["nshow"] - - combinations = list( - itertools.product( - grid_axes["primary"], - grid_axes["azimuth_angle"], - grid_axes["zenith_angle"], - grid_axes["model_version"], - grid_axes["corsika_le_interaction"], - grid_axes["corsika_he_interaction"], - energy_ranges, - ) - ) - - number_of_runs = args_dict.get("number_of_runs", 1) - run_number = int(args_dict.get("run_number") or 1) - - job_specs = [] - for label in image_labels: - params_label = _sanitize_label_for_params(label) - row_index = 0 - for ( - primary, - azimuth, - zenith, - model_version, - corsika_le, - corsika_he, - energy_range_pair, - ) in combinations: - selected_energy_range_pair = energy_range_pair - selected_core_scatter_max = core_scatter[1] - selected_nshow = nshow - - if corsika_limits is not None: - selected_energy_range_pair = _get_energy_range_for_zenith_angle( - zenith, energy_range_pair, corsika_limits - ) - if selected_energy_range_pair is None: - continue - selected_core_scatter_max = _get_core_scatter_max_for_zenith_angle( - zenith, core_scatter, corsika_limits - ) - selected_nshow = _get_nshow_for_energy_range_and_zenith_angle( - selected_energy_range_pair, zenith, nshow, corsika_limits - ) - - for _ in range(number_of_runs): - job_specs.append( - { - "apptainer_label": str(label), - "primary": primary, - "azimuth_angle": azimuth, - "zenith_angle": zenith, - "model_version": model_version, - "array_layout_name": args_dict.get("array_layout_name"), - "corsika_le_interaction": corsika_le, - "corsika_he_interaction": corsika_he, - "energy_min": selected_energy_range_pair[0], - "energy_max": selected_energy_range_pair[1], - "core_scatter_max": selected_core_scatter_max, - "nshow": selected_nshow, - "pack_for_grid_register": f"{base_pack_dir}/{params_label}", - "run_number": run_number + row_index, - } - ) - row_index += 1 - return job_specs - - def _group_job_specs_by_label(job_specs): """Group job specs by apptainer image label.""" grouped = {} for job_spec in job_specs: - label = job_spec["apptainer_label"] + label = job_spec["image_label"] grouped.setdefault(label, []).append(job_spec) return grouped @@ -322,7 +139,7 @@ def _write_params_file(params_file_path, label_job_specs): """Write parameter file consumed by HTCondor queue-from syntax.""" with open(params_file_path, "w", encoding="utf-8") as params_file_handle: for job_spec in label_job_specs: - array_layout_name = _resolve_array_layout_name( + array_layout_name = resolve_array_layout_name( job_spec["array_layout_name"], job_spec["model_version"] ) @@ -337,7 +154,7 @@ def _write_params_file(params_file_path, label_job_specs): ) row = [ - _format_param_value(job_spec["apptainer_label"], "apptainer_label"), + _format_param_value(job_spec["image_label"], "apptainer_label"), _format_param_value(job_spec["primary"], "primary"), _format_param_value(job_spec["azimuth_angle"], "azimuth_angle"), _format_param_value(job_spec["zenith_angle"], "zenith_angle"), @@ -368,7 +185,7 @@ def generate_submission_script(args_dict): Arguments dictionary. """ apptainer_images = _resolve_apptainer_images(args_dict["apptainer_image"]) - job_specs = _build_job_specs(args_dict, list(apptainer_images.keys())) + job_specs = build_job_specs(args_dict, list(apptainer_images.keys())) grouped_job_specs = _group_job_specs_by_label(job_specs) work_dir = Path(args_dict["output_path"]) diff --git a/src/simtools/production_configuration/job_spec_builder.py b/src/simtools/production_configuration/job_spec_builder.py new file mode 100644 index 0000000000..341ee913f1 --- /dev/null +++ b/src/simtools/production_configuration/job_spec_builder.py @@ -0,0 +1,245 @@ +"""Build backend-agnostic job specifications for production submissions.""" + +import ast +import itertools + +import numpy as np +from astropy import units as u + +import simtools.version as simtools_version +from simtools.configuration import defaults + +_GRID_AXES = [ + "primary", + "azimuth_angle", + "zenith_angle", + "model_version", + "corsika_le_interaction", + "corsika_he_interaction", +] + +_GRID_AXIS_DEFAULTS = { + "corsika_le_interaction": defaults.CORSIKA_LE_INTERACTION, + "corsika_he_interaction": defaults.CORSIKA_HE_INTERACTION, +} + + +def normalize_to_list(value): + """Normalize scalar values to lists of length one.""" + if isinstance(value, list): + return value + if isinstance(value, tuple): + return list(value) + return [value] + + +def normalize_grid_axes(args_dict): + """Return normalized grid axes for cartesian product expansion.""" + return { + axis: ( + normalize_to_list(args_dict[axis]) + if axis in args_dict and args_dict[axis] is not None + else [_GRID_AXIS_DEFAULTS[axis]] + if axis in _GRID_AXIS_DEFAULTS + else [None] + ) + for axis in _GRID_AXES + } + + +def normalize_energy_ranges(energy_range): + """Normalize energy range argument to a list of ``(e_min, e_max)`` pairs.""" + if isinstance(energy_range, tuple) and len(energy_range) == 2: + return [energy_range] + + if isinstance(energy_range, list): + if len(energy_range) == 2 and all(hasattr(item, "to") for item in energy_range): + return [(energy_range[0], energy_range[1])] + if all(isinstance(item, (list, tuple)) and len(item) == 2 for item in energy_range): + return [tuple(item) for item in energy_range] + + raise ValueError( + "energy_range must be one pair (e_min, e_max) or a list of (e_min, e_max) pairs." + ) + + +def resolve_array_layout_name(array_layout_name, model_version): + """Resolve array layout configuration for a specific model version.""" + if isinstance(array_layout_name, list) and len(array_layout_name) == 1: + array_layout_name = array_layout_name[0] + + if isinstance(array_layout_name, str) and array_layout_name.strip().startswith("{"): + try: + parsed_layout = ast.literal_eval(array_layout_name) + if isinstance(parsed_layout, dict): + array_layout_name = parsed_layout + except (SyntaxError, ValueError): + return array_layout_name + + if not isinstance(array_layout_name, dict) or list(array_layout_name) != ["by_version"]: + return array_layout_name + + resolved = simtools_version.resolve_by_version( + {"array_layout_name": array_layout_name}, model_version + ) + return resolved["array_layout_name"] + + +def get_energy_range_for_zenith_angle(zenith_angle, energy_range_pair, corsika_limits): + """Return a zenith-dependent energy range pair or None to skip the simulation step.""" + _ = (zenith_angle, corsika_limits) + return energy_range_pair + + +def get_core_scatter_max_for_zenith_angle(zenith_angle, core_scatter, corsika_limits): + """Return zenith-dependent max core-scatter value.""" + _ = (zenith_angle, corsika_limits) + return core_scatter[1] + + +def calculate_log_energy_midpoint(energy_range_pair): + """Return the geometric-mean energy for an energy range pair.""" + energy_min, energy_max = energy_range_pair + + if not isinstance(energy_min, u.Quantity) or not isinstance(energy_max, u.Quantity): + raise TypeError("energy_range_pair must contain astropy Quantity values.") + + energy_min_tev = energy_min.to(u.TeV) + energy_max_tev = energy_max.to(u.TeV) + + if energy_min_tev <= 0 * u.TeV or energy_max_tev <= 0 * u.TeV: + raise ValueError("Energy range values must be strictly positive.") + + mean_log_energy = np.mean( + [ + np.log10(energy_min_tev.value), + np.log10(energy_max_tev.value), + ] + ) + return 10**mean_log_energy * u.TeV + + +def calculate_scaled_nshow( + energy_range_pair, + baseline_nshow, + nshow_power_index=None, + reference_energy=None, +): + """Return an energy-dependent nshow value.""" + if baseline_nshow < 1: + raise ValueError("baseline_nshow must be a positive integer.") + + if nshow_power_index is None: + return baseline_nshow + + if reference_energy is None: + raise ValueError("reference_energy is required when nshow_power_index is configured.") + + midpoint_energy = calculate_log_energy_midpoint(energy_range_pair) + scaling_factor = (midpoint_energy / reference_energy.to(midpoint_energy.unit)).to_value( + u.dimensionless_unscaled + ) ** nshow_power_index + scaled_nshow = int(np.ceil(baseline_nshow * scaling_factor)) + + if scaled_nshow < 1: + raise ValueError("Scaled nshow must be at least 1.") + + return scaled_nshow + + +def _select_energy_and_core_scatter_for_job( + zenith, energy_range_pair, core_scatter, corsika_limits +): + """Return selected energy range and core scatter maximum for a job spec row.""" + selected_energy_range_pair = energy_range_pair + selected_core_scatter_max = core_scatter[1] + + if corsika_limits is None: + return selected_energy_range_pair, selected_core_scatter_max + + selected_energy_range_pair = get_energy_range_for_zenith_angle( + zenith, energy_range_pair, corsika_limits + ) + if selected_energy_range_pair is None: + return None, None + + selected_core_scatter_max = get_core_scatter_max_for_zenith_angle( + zenith, core_scatter, corsika_limits + ) + return selected_energy_range_pair, selected_core_scatter_max + + +def build_job_specs(args_dict, image_labels): + """Build backend-agnostic job specs from comparison and production grids.""" + grid_axes = normalize_grid_axes(args_dict) + energy_ranges = normalize_energy_ranges(args_dict["energy_range"]) + base_pack_dir = args_dict.get("simulation_output") or "simtools-output" + corsika_limits = args_dict.get("corsika_limits") + core_scatter = args_dict["core_scatter"] + nshow = args_dict["nshow"] + nshow_power_index = args_dict.get("nshow_power_index") + reference_energy = args_dict.get("nshow_reference_energy") + if nshow_power_index is not None and reference_energy is not None: + reference_energy = u.Quantity(reference_energy) + + combinations = list( + itertools.product( + grid_axes["primary"], + grid_axes["azimuth_angle"], + grid_axes["zenith_angle"], + grid_axes["model_version"], + grid_axes["corsika_le_interaction"], + grid_axes["corsika_he_interaction"], + energy_ranges, + ) + ) + + number_of_runs = args_dict.get("number_of_runs", 1) + run_number = int(args_dict.get("run_number") or 1) + + job_specs = [] + for label in image_labels: + row_index = 0 + for ( + primary, + azimuth, + zenith, + model_version, + corsika_le, + corsika_he, + energy_range_pair, + ) in combinations: + ( + selected_energy_range_pair, + selected_core_scatter_max, + ) = _select_energy_and_core_scatter_for_job( + zenith, energy_range_pair, core_scatter, corsika_limits + ) + if selected_energy_range_pair is None: + continue + + selected_nshow = calculate_scaled_nshow( + selected_energy_range_pair, nshow, nshow_power_index, reference_energy + ) + + for _ in range(number_of_runs): + job_specs.append( + { + "image_label": str(label), + "primary": primary, + "azimuth_angle": azimuth, + "zenith_angle": zenith, + "model_version": model_version, + "array_layout_name": args_dict.get("array_layout_name"), + "corsika_le_interaction": corsika_le, + "corsika_he_interaction": corsika_he, + "energy_min": selected_energy_range_pair[0], + "energy_max": selected_energy_range_pair[1], + "core_scatter_max": selected_core_scatter_max, + "nshow": selected_nshow, + "pack_for_grid_register": f"{base_pack_dir}/{label!s}", + "run_number": run_number + row_index, + } + ) + row_index += 1 + return job_specs diff --git a/tests/unit_tests/applications/test_simulate_prod_htcondor_generator.py b/tests/unit_tests/applications/test_simulate_prod_htcondor_generator.py index bfa76c9bdb..38288feef8 100644 --- a/tests/unit_tests/applications/test_simulate_prod_htcondor_generator.py +++ b/tests/unit_tests/applications/test_simulate_prod_htcondor_generator.py @@ -1,8 +1,11 @@ #!/usr/bin/python3 +import argparse from types import SimpleNamespace from unittest.mock import patch +import pytest + import simtools.applications.simulate_prod_htcondor_generator as app @@ -24,3 +27,22 @@ def test_main_uses_standard_build_application( "simulation_configuration": {"software": None, "corsika_configuration": ["all"]}, } mock_generate_submission_script.assert_called_once_with({"output_path": "htcondor_submit"}) + + +def test_add_arguments_registers_nshow_scaling_arguments(): + parser = argparse.ArgumentParser() + + app._add_arguments(parser) + args = parser.parse_args( + [ + "--number_of_runs", + "1", + "--nshow_power_index", + "-0.5", + "--nshow_reference_energy", + "100 GeV", + ] + ) + + assert args.nshow_power_index == pytest.approx(-0.5) + assert args.nshow_reference_energy == "100 GeV" diff --git a/tests/unit_tests/job_execution/test_htcondor_script_generator.py b/tests/unit_tests/job_execution/test_htcondor_script_generator.py index ffef52b577..3b94547ddb 100644 --- a/tests/unit_tests/job_execution/test_htcondor_script_generator.py +++ b/tests/unit_tests/job_execution/test_htcondor_script_generator.py @@ -5,7 +5,7 @@ import pytest from simtools.job_execution.htcondor_script_generator import ( - _build_job_specs, + _format_param_value, _format_quantity, _get_submit_file, _get_submit_script, @@ -14,6 +14,7 @@ _write_params_file, generate_submission_script, ) +from simtools.production_configuration.job_spec_builder import build_job_specs @pytest.fixture @@ -66,6 +67,44 @@ def test_generate_submission_script(mock_is_file, mock_chmod, mock_open, mock_mk mock_chmod.assert_called_once_with(0o755) +@mock.patch("simtools.job_execution.htcondor_script_generator.Path.mkdir") +@mock.patch("simtools.job_execution.htcondor_script_generator.open", new_callable=mock.mock_open) +@mock.patch("simtools.job_execution.htcondor_script_generator.Path.chmod") +@mock.patch("simtools.job_execution.htcondor_script_generator.Path.is_file", return_value=True) +def test_generate_submission_script_writes_label_specific_files( + mock_is_file, mock_chmod, mock_open, mock_mkdir, args_dict +): + args_dict["output_path"] = "/test_output" + args_dict["htcondor_log_path"] = "/custom_logs" + args_dict["number_of_runs"] = 1 + args_dict["apptainer_image"] = { + "prod 7.0.0": "/path/to/prod.sif", + "beta": "/path/to/beta.sif", + } + + generate_submission_script(args_dict) + + work_dir = Path(args_dict["output_path"]) + + mock_is_file.assert_has_calls( + [mock.call(), mock.call()], + any_order=False, + ) + mock_open.assert_any_call( + work_dir / "simulate_prod.submit.prod_7.0.0.condor", "w", encoding="utf-8" + ) + mock_open.assert_any_call( + work_dir / "simulate_prod.submit.prod_7.0.0.params.txt", "w", encoding="utf-8" + ) + mock_open.assert_any_call(work_dir / "simulate_prod.submit.beta.condor", "w", encoding="utf-8") + mock_open.assert_any_call( + work_dir / "simulate_prod.submit.beta.params.txt", "w", encoding="utf-8" + ) + mock_open.assert_any_call(work_dir / "simulate_prod.submit.sh", "w", encoding="utf-8") + mock_mkdir.assert_any_call(parents=True, exist_ok=True) + mock_chmod.assert_called_once_with(0o755) + + @mock.patch("simtools.job_execution.htcondor_script_generator.Path.is_file", return_value=False) @mock.patch("simtools.job_execution.htcondor_script_generator.open", new_callable=mock.mock_open) def test_generate_submission_script_raises_for_missing_apptainer_image( @@ -191,6 +230,17 @@ def test_resolve_apptainer_images_empty_dict(): _resolve_apptainer_images({}) +def test_resolve_apptainer_images_raises_for_truthy_empty_mapping(): + class TruthyEmptyDict(dict): + def __bool__(self): + return True + + with pytest.raises( + ValueError, match="At least one apptainer image label/path must be configured" + ): + _resolve_apptainer_images(TruthyEmptyDict()) + + def test_resolve_apptainer_images_invalid_type(tmp_test_directory): with pytest.raises( TypeError, match="apptainer_image must be a string path or a label-to-path dictionary" @@ -223,60 +273,17 @@ def test_format_quantity_full_coverage(): assert unit is None +def test_format_param_value_raises_for_missing_required_value(): + with pytest.raises(ValueError, match="Missing required value for field 'primary'"): + _format_param_value(None, "primary") + + def test_sanitize_label_for_filename(): assert _sanitize_label_for_filename(" my label:7/0*0? ") == "my_label_7_0_0_" assert _sanitize_label_for_filename("v7.0.0-beta_1") == "v7.0.0-beta_1" assert _sanitize_label_for_filename(42) == "42" -def test_build_job_specs_expands_model_version_list(args_dict): - args_dict["model_version"] = ["6.3.0", "7.0.0"] - - job_specs = _build_job_specs(args_dict, ["7.0.0"]) - model_versions = {job_spec["model_version"] for job_spec in job_specs} - - assert model_versions == {"6.3.0", "7.0.0"} - assert len(job_specs) == 2 * args_dict["number_of_runs"] - - -def test_build_job_specs_expands_energy_range_list_of_pairs(args_dict): - args_dict["number_of_runs"] = 1 - args_dict["energy_range"] = [ - (30 * u.GeV, 30 * u.GeV), - (300 * u.GeV, 300 * u.GeV), - ] - - job_specs = _build_job_specs(args_dict, ["7.0.0"]) - energy_pairs = {(job_spec["energy_min"], job_spec["energy_max"]) for job_spec in job_specs} - - assert len(job_specs) == 2 - assert energy_pairs == { - (30 * u.GeV, 30 * u.GeV), - (300 * u.GeV, 300 * u.GeV), - } - - -def test_build_job_specs_uses_default_interaction_models_when_missing(args_dict): - args_dict.pop("corsika_le_interaction") - args_dict.pop("corsika_he_interaction") - - job_specs = _build_job_specs(args_dict, ["7.0.0"]) - - assert {job_spec["corsika_le_interaction"] for job_spec in job_specs} == {"urqmd"} - assert {job_spec["corsika_he_interaction"] for job_spec in job_specs} == {"epos"} - - -def test_build_job_specs_increments_run_number(args_dict): - args_dict["number_of_runs"] = 2 - args_dict["run_number"] = 10 - args_dict["model_version"] = ["6.3.0", "7.0.0"] - - job_specs = _build_job_specs(args_dict, ["7.0.0"]) - run_numbers = [job_spec["run_number"] for job_spec in job_specs] - - assert run_numbers == [10, 11, 12, 13] - - def test_write_params_file_resolves_array_layout_name_by_model_version( args_dict, tmp_test_directory ): @@ -290,7 +297,7 @@ def test_write_params_file_resolves_array_layout_name_by_model_version( } params_file_path = Path(tmp_test_directory) / "params.txt" - job_specs = _build_job_specs(args_dict, ["7.0.0"]) + job_specs = build_job_specs(args_dict, ["7.0.0"]) _write_params_file(params_file_path, job_specs) @@ -313,7 +320,7 @@ def test_write_params_file_resolves_stringified_by_version_layout(args_dict, tmp ) params_file_path = Path(tmp_test_directory) / "params.txt" - job_specs = _build_job_specs(args_dict, ["7.0.0"]) + job_specs = build_job_specs(args_dict, ["7.0.0"]) _write_params_file(params_file_path, job_specs) @@ -327,7 +334,7 @@ def test_write_params_file_keeps_energy_units(tmp_test_directory): params_file_path = Path(tmp_test_directory) / "params.txt" label_job_specs = [ { - "apptainer_label": "7.0.0", + "image_label": "7.0.0", "primary": "gamma", "azimuth_angle": 0 * u.deg, "zenith_angle": 20 * u.deg, @@ -361,7 +368,7 @@ def test_write_params_file_replaces_whitespace_in_apptainer_label(tmp_test_direc params_file_path = Path(tmp_test_directory) / "params.txt" label_job_specs = [ { - "apptainer_label": "grid label 7.0.0", + "image_label": "grid label 7.0.0", "primary": "gamma", "azimuth_angle": 0 * u.deg, "zenith_angle": 20 * u.deg, @@ -384,32 +391,3 @@ def test_write_params_file_replaces_whitespace_in_apptainer_label(tmp_test_direc "grid_label_7.0.0 gamma 0.0 20.0 30.0 GeV 10.0 TeV 200.0 m " "1000 7.0.0 CTAO-North-Alpha urqmd epos 10 simtools-output/grid_label_7.0.0\n" ) - - -@mock.patch( - "simtools.job_execution.htcondor_script_generator._get_energy_range_for_zenith_angle", - return_value=None, -) -def test_build_job_specs_skips_entries_when_energy_range_is_none(mock_energy_range, args_dict): - args_dict["corsika_limits"] = "limits.ecsv" - args_dict["number_of_runs"] = 1 - - job_specs = _build_job_specs(args_dict, ["7.0.0"]) - - assert job_specs == [] - mock_energy_range.assert_called_once() - - -@mock.patch( - "simtools.job_execution.htcondor_script_generator._get_nshow_for_energy_range_and_zenith_angle", - return_value=777, -) -def test_build_job_specs_uses_dummy_nshow_when_corsika_limits_set(mock_nshow, args_dict): - args_dict["corsika_limits"] = "limits.ecsv" - args_dict["number_of_runs"] = 1 - - job_specs = _build_job_specs(args_dict, ["7.0.0"]) - - assert len(job_specs) == 1 - assert job_specs[0]["nshow"] == 777 - mock_nshow.assert_called_once() diff --git a/tests/unit_tests/production_configuration/test_job_spec_builder.py b/tests/unit_tests/production_configuration/test_job_spec_builder.py new file mode 100644 index 0000000000..1e5321df13 --- /dev/null +++ b/tests/unit_tests/production_configuration/test_job_spec_builder.py @@ -0,0 +1,252 @@ +from unittest import mock + +import astropy.units as u +import pytest +from astropy.tests.helper import assert_quantity_allclose + +from simtools.production_configuration.job_spec_builder import ( + build_job_specs, + calculate_log_energy_midpoint, + calculate_scaled_nshow, + normalize_energy_ranges, + normalize_grid_axes, + normalize_to_list, + resolve_array_layout_name, +) + + +@pytest.fixture +def args_dict(): + return { + "azimuth_angle": 45 * u.deg, + "zenith_angle": 20 * u.deg, + "energy_range": [1 * u.GeV, 10 * u.GeV], + "core_scatter": [0, 100 * u.m], + "model_version": "v1.0", + "array_layout_name": "test_layout", + "primary": "gamma", + "nshow": 1000, + "run_number": 1, + "number_of_runs": 10, + "corsika_le_interaction": "urqmd", + "corsika_he_interaction": "epos", + } + + +def test_normalize_energy_ranges_expands_list_of_pairs(): + energy_ranges = normalize_energy_ranges( + [ + (30 * u.GeV, 30 * u.GeV), + (300 * u.GeV, 300 * u.GeV), + ] + ) + + assert len(energy_ranges) == 2 + for actual, expected in zip( + energy_ranges, + [ + (30 * u.GeV, 30 * u.GeV), + (300 * u.GeV, 300 * u.GeV), + ], + ): + assert_quantity_allclose(actual[0], expected[0]) + assert_quantity_allclose(actual[1], expected[1]) + + +def test_normalize_to_list_converts_tuple_values(): + assert normalize_to_list((1, 2)) == [1, 2] + + +def test_normalize_grid_axes_uses_defaults_and_none_for_missing_axes(): + grid_axes = normalize_grid_axes({"primary": "gamma"}) + + assert grid_axes["primary"] == ["gamma"] + assert grid_axes["azimuth_angle"] == [None] + assert grid_axes["zenith_angle"] == [None] + assert grid_axes["model_version"] == [None] + assert grid_axes["corsika_le_interaction"] == ["urqmd"] + assert grid_axes["corsika_he_interaction"] == ["epos"] + + +def test_normalize_energy_ranges_accepts_single_tuple_pair(): + energy_ranges = normalize_energy_ranges((30 * u.GeV, 300 * u.GeV)) + assert len(energy_ranges) == 1 + assert_quantity_allclose(energy_ranges[0][0], 30 * u.GeV) + assert_quantity_allclose(energy_ranges[0][1], 300 * u.GeV) + + +def test_normalize_energy_ranges_raises_for_invalid_shape(): + with pytest.raises(ValueError, match="energy_range must be one pair"): + normalize_energy_ranges([30 * u.GeV, 300]) + + +def test_calculate_log_energy_midpoint(): + midpoint_energy = calculate_log_energy_midpoint((1 * u.GeV, 100 * u.GeV)) + + assert midpoint_energy.to_value(u.GeV) == pytest.approx(10.0) + + +def test_calculate_log_energy_midpoint_raises_for_non_quantity_values(): + with pytest.raises(TypeError, match="energy_range_pair must contain astropy Quantity values"): + calculate_log_energy_midpoint((1, 100 * u.GeV)) + + +def test_calculate_log_energy_midpoint_raises_for_non_positive_values(): + with pytest.raises(ValueError, match="Energy range values must be strictly positive"): + calculate_log_energy_midpoint((0 * u.GeV, 100 * u.GeV)) + + +def test_calculate_scaled_nshow_returns_baseline_without_power_index(): + scaled_nshow = calculate_scaled_nshow((10 * u.GeV, 100 * u.GeV), 50) + + assert scaled_nshow == 50 + + +def test_calculate_scaled_nshow_scales_against_reference_energy(): + scaled_nshow = calculate_scaled_nshow( + (100 * u.GeV, 100 * u.GeV), + 100, + nshow_power_index=-1.0, + reference_energy=10 * u.GeV, + ) + + assert scaled_nshow == 10 + + +def test_calculate_scaled_nshow_uses_ceil_for_fractional_result(): + scaled_nshow = calculate_scaled_nshow( + (100 * u.GeV, 100 * u.GeV), + 10, + nshow_power_index=-2.0, + reference_energy=10 * u.GeV, + ) + + assert scaled_nshow == 1 + + +def test_calculate_scaled_nshow_raises_for_invalid_baseline(): + with pytest.raises(ValueError, match="baseline_nshow must be a positive integer"): + calculate_scaled_nshow((10 * u.GeV, 100 * u.GeV), 0) + + +def test_calculate_scaled_nshow_requires_reference_energy_when_scaled(): + with pytest.raises(ValueError, match="reference_energy is required"): + calculate_scaled_nshow((10 * u.GeV, 100 * u.GeV), 10, nshow_power_index=-1.0) + + +def test_calculate_scaled_nshow_raises_when_scaled_result_drops_below_one(): + with pytest.raises(ValueError, match="Scaled nshow must be at least 1"): + calculate_scaled_nshow( + (10 * u.GeV, 10 * u.GeV), + 1, + nshow_power_index=1.0, + reference_energy=-100 * u.GeV, + ) + + +def test_build_job_specs_expands_model_version_list(args_dict): + args_dict["model_version"] = ["6.3.0", "7.0.0"] + + job_specs = build_job_specs(args_dict, ["7.0.0"]) + model_versions = {job_spec["model_version"] for job_spec in job_specs} + + assert model_versions == {"6.3.0", "7.0.0"} + assert len(job_specs) == 2 * args_dict["number_of_runs"] + + +def test_build_job_specs_scales_nshow_by_energy_range(args_dict): + args_dict["number_of_runs"] = 1 + args_dict["nshow"] = 100 + args_dict["nshow_power_index"] = -1.0 + args_dict["nshow_reference_energy"] = 10 * u.GeV + args_dict["energy_range"] = [ + (10 * u.GeV, 10 * u.GeV), + (100 * u.GeV, 100 * u.GeV), + ] + + job_specs = build_job_specs(args_dict, ["7.0.0"]) + nshow_values = [job_spec["nshow"] for job_spec in job_specs] + assert nshow_values[0] == pytest.approx(100) + assert nshow_values[1] == pytest.approx(10) + + +def test_build_job_specs_requires_reference_energy_for_nshow_scaling(args_dict): + args_dict["number_of_runs"] = 1 + args_dict["nshow_power_index"] = -1.0 + args_dict["nshow_reference_energy"] = None + + with pytest.raises(ValueError, match="reference_energy is required"): + build_job_specs(args_dict, ["7.0.0"]) + + +def test_build_job_specs_uses_default_interaction_models_when_missing(args_dict): + args_dict.pop("corsika_le_interaction") + args_dict.pop("corsika_he_interaction") + + job_specs = build_job_specs(args_dict, ["7.0.0"]) + + assert {job_spec["corsika_le_interaction"] for job_spec in job_specs} == {"urqmd"} + assert {job_spec["corsika_he_interaction"] for job_spec in job_specs} == {"epos"} + + +def test_build_job_specs_increments_run_number(args_dict): + args_dict["number_of_runs"] = 2 + args_dict["run_number"] = 10 + args_dict["model_version"] = ["6.3.0", "7.0.0"] + + job_specs = build_job_specs(args_dict, ["7.0.0"]) + run_numbers = [job_spec["run_number"] for job_spec in job_specs] + + assert run_numbers == [10, 11, 12, 13] + + +@mock.patch( + "simtools.production_configuration.job_spec_builder.get_energy_range_for_zenith_angle", + return_value=None, +) +def test_build_job_specs_skips_entries_when_energy_range_is_none(mock_energy_range, args_dict): + args_dict["corsika_limits"] = "limits.ecsv" + args_dict["number_of_runs"] = 1 + + job_specs = build_job_specs(args_dict, ["7.0.0"]) + + assert job_specs == [] + mock_energy_range.assert_called_once() + + +@mock.patch( + "simtools.production_configuration.job_spec_builder.calculate_scaled_nshow", + return_value=777, +) +def test_build_job_specs_uses_dummy_nshow_when_corsika_limits_set(mock_nshow, args_dict): + args_dict["corsika_limits"] = "limits.ecsv" + args_dict["number_of_runs"] = 1 + + job_specs = build_job_specs(args_dict, ["7.0.0"]) + + assert len(job_specs) == 1 + assert job_specs[0]["nshow"] == 777 + mock_nshow.assert_called_once() + + +def test_resolve_array_layout_name_resolves_stringified_by_version_layout(): + array_layout_name = str( + { + "by_version": { + "<7.0.0": "alpha", + ">=7.0.0": "CTAO-North-Alpha", + } + } + ) + + assert resolve_array_layout_name(array_layout_name, "7.0.0") == "CTAO-North-Alpha" + + +def test_resolve_array_layout_name_unwraps_single_item_list(): + assert resolve_array_layout_name(["CTAO-North-Alpha"], "7.0.0") == "CTAO-North-Alpha" + + +def test_resolve_array_layout_name_keeps_invalid_stringified_dict(): + invalid_layout = "{not valid" + + assert resolve_array_layout_name(invalid_layout, "7.0.0") == invalid_layout