Refactor generate velocity model parameters

pyproject.toml

@@ -17,7 +17,6 @@ dependencies = [
   "oq_wrapper>=2025.12.3",
   "qcore-utils>=2025.12.1",
   "source_modelling>=2025.12.1",
-
   # Data Formats
   "geopandas",
   "pandas[parquet, hdf5]",
@@ -29,20 +28,22 @@ dependencies = [
   "numpy",
   "scipy",
   "shapely",
-
   # CLI
   "tqdm",
   "typer",
-
   # Misc.
   "requests",  # For gcmt-to-realisation
   "schema",    # For loading realisations
   "structlog", # Logging.
   "psutil",    # To get the CPU affinity for jobs
+
 ]
 
 [project.optional-dependencies]
-test = ["pytest"]
+test = [
+  "pytest>=6.0.0",            # required for the tool.pytest section
+  "hypothesis[numpy]>=6.0.0",
+]
 types = ["pandas-stubs", "types-geopandas", "types-requests", "scipy-stubs"]
 dev = ["ruff", "deptry", "ty", "numpydoc"]
 
@@ -53,7 +54,7 @@ pep621_dev_dependency_groups = ["test", "dev"]
 nshm2022-to-realisation = "workflow.scripts.nshm2022_to_realisation:app"
 gcmt-to-realisation = "workflow.scripts.gcmt_to_realisation:app"
 realisation-to-srf = "workflow.scripts.realisation_to_srf:app"
-generate-velocity-model-parameters = "workflow.scripts.generate_velocity_model_parameters:app"
+generate-domain = "workflow.scripts.generate_domain:app"
 generate-velocity-model = "workflow.scripts.generate_velocity_model:app"
 generate-station-coordinates = "workflow.scripts.generate_station_coordinates:app"
 generate-model-coordinates = "workflow.scripts.generate_model_coordinates:app"
@@ -76,6 +77,8 @@ workflow = "workflow"
 
 [tool.setuptools_scm]
 
+[tool.pytest]
+markers = ["slow: mark test as slow."]
 
 [tool.ruff.lint]
 extend-select = [

tests/test_cli.py

@@ -1,60 +1,59 @@
-from collections.abc import Callable
+import importlib
+import pkgutil
+from types import ModuleType
 
-import pytest
+from pytest import Metafunc
 from typer import Typer
 from typer.testing import CliRunner
 
-from workflow.scripts import (
-    bb_sim,
-    check_domain,
-    check_srf,
-    copy_velocity_model_parameters,
-    create_e3d_par,
-    gcmt_auto_simulate,
-    gcmt_to_realisation,
-    generate_rupture_propagation,
-    generate_station_coordinates,
-    generate_stoch,
-    generate_velocity_model,
-    generate_velocity_model_parameters,
-    hf_sim,
-    im_calc,
-    import_realisation,
-    lf_to_xarray,
-    nshm2022_to_realisation,
-    realisation_to_srf,
-)
-
-
-@pytest.mark.parametrize(
-    "script",
-    [
-        bb_sim,
-        check_domain,
-        check_srf,
-        copy_velocity_model_parameters,
-        create_e3d_par,
-        gcmt_auto_simulate,
-        gcmt_to_realisation,
-        generate_rupture_propagation,
-        generate_station_coordinates,
-        generate_stoch,
-        generate_velocity_model,
-        generate_velocity_model_parameters,
-        lf_to_xarray,
-        hf_sim,
-        im_calc,
-        import_realisation,
-        nshm2022_to_realisation,
-        realisation_to_srf,
-    ],
-)
-def test_invocation_of_script(script: Callable) -> None:
-    """Basic check that the scripts can be invoked."""
+import workflow.scripts as scripts_package
+
+EXCLUDE_MODULES = set()
+
+
+def collect_script_modules() -> list[ModuleType]:
+    """
+    Dynamically discovers all modules within the workflow.scripts package.
+    """
+    modules = []
+    # Iterates through all modules in the package directory
+    for loader, module_name, is_pkg in pkgutil.iter_modules(scripts_package.__path__):
+        if module_name in EXCLUDE_MODULES:
+            continue
+
+        # Construct full module path (e.g., workflow.scripts.bb_sim)
+        full_module_name = f"{scripts_package.__name__}.{module_name}"
+        module = importlib.import_module(full_module_name)
+        modules.append(module)
+    return modules
+
+
+def pytest_generate_tests(metafunc: Metafunc) -> None:
+    """
+    Generate tests dynamically based on discovered modules.
+    """
+    if "script_module" in metafunc.fixturenames:
+        found_modules = collect_script_modules()
+        # Create readable IDs from the module names (e.g., 'bb_sim')
+        ids = [m.__name__.split(".")[-1] for m in found_modules]
+        metafunc.parametrize("script_module", found_modules, ids=ids)
+
+
+def test_invocation_of_script(script_module: ModuleType) -> None:
+    """
+    Test that each discovered module has a Typer app and responds to --help.
+    """
     runner = CliRunner()
-    # The following satisifies the type checker.
-    app = getattr(script, "app", None)
-    assert isinstance(app, Typer)
+
+    assert hasattr(script_module, "app"), (
+        f"Module {script_module.__name__} is missing the 'app' attribute."
+    )
+
+    app = getattr(script_module, "app")
+
+    assert isinstance(app, Typer), (
+        f"'app' in {script_module.__name__} should be a Typer instance, but got {type(app)}."
+    )
 
     result = runner.invoke(app, ["--help"])
     assert result.exit_code == 0

tests/test_generate_domain.py

@@ -0,0 +1,112 @@
+from types import SimpleNamespace
+
+import numpy as np
+import pytest
+import shapely
+from hypothesis import given
+from hypothesis import strategies as st
+
+from source_modelling import sources
+from workflow.realisations import (
+    Magnitudes,
+    Rakes,
+    SourceConfig,
+    VelocityModelParameters,
+)
+from workflow.scripts import generate_domain
+
+
+# Slow because of openquake import
+@pytest.mark.slow
+def test_significant_duration_calculation() -> None:
+    """Basic integration test for OQW, tests that the interface works as we expect still."""
+    ds595 = generate_domain.get_significant_duration(
+        magnitude=6.5, distance=100.0, vs30=500.0, rake=180.0, z1pt0=5.0
+    )
+    assert isinstance(ds595, float)  # Should not be lying about the type
+    assert np.isfinite(ds595), "Ds595 is invalid"
+    # Check that ds595 is not in log-space or similar
+    assert ds595 > 0, "Ds595 should be positive"
+    assert ds595 < 1000, "Ds595 unrealistically high"
+    # Combined with type checking this should be enough to catch most of these problems.
+
+
+def test_simulation_max_depth_increases_with_mw() -> None:
+    depth_mw5p0 = generate_domain.simulation_max_depth(magnitude=5.0, bottom_depth=10.0)
+    depth_mw7p0 = generate_domain.simulation_max_depth(magnitude=6.0, bottom_depth=10.0)
+    depth_mw9p0 = generate_domain.simulation_max_depth(magnitude=9.0, bottom_depth=10.0)
+    assert depth_mw5p0 < depth_mw7p0 < depth_mw9p0, (
+        "Depths do not increase with magnitude"
+    )
+
+
+@given(
+    magnitude=st.floats(min_value=3.5, max_value=9.0),
+    depth=st.floats(min_value=1.0, max_value=250),
+)
+def test_simulation_max_depth_more_than_bottom_depth(
+    magnitude: float, depth: float
+) -> None:
+    simulation_depth = generate_domain.simulation_max_depth(magnitude, depth)
+    assert np.isfinite(simulation_depth), "Invalid simulation depth"
+    assert depth <= simulation_depth <= 350, "Sensible simulation depths are applied"
+
+
+def test_estimate_domain_contains_fault_geometry() -> None:
+    fault_coords = [(100000, 100000), (110000, 100000)]
+    fault_geom = shapely.LineString(fault_coords)
+
+    mock_fault = SimpleNamespace(geometry=fault_geom)
+    source_config = SimpleNamespace(source_geometries={"fault_a": mock_fault})
+
+    rrups = {"fault_a": 5000.0}
+
+    nz_outline = shapely.box(0, 0, 500000, 500000)
+
+    result_domain = generate_domain.estimate_domain(
+        source_config=source_config,  # type: ignore[invalid-argument-type]
+        rrups=rrups,
+        nz_outline=nz_outline,
+        fault_buffer=2000.0,
+    )
+
+    assert result_domain.polygon.contains(fault_geom), (
+        f"Domain polygon should contain the original fault geometry.\n"
+        f"Fault bounds: {fault_geom.bounds}\n"
+        f"Domain bounds: {result_domain.polygon.bounds}"
+    )
+
+
+# Slow because of openquake import
+@pytest.mark.slow
+def test_generate_domain() -> None:
+    """Basic E2E test to check that domain generation works without crashing or producing a silly domain."""
+    source = sources.Point(
+        np.array([-43.0, -172.0, 10000.0]),
+        length_m=1000,
+        width_m=1000,
+        strike=90.0,
+        dip=45.0,
+        dip_dir=180.0,
+    )
+    source_config = SourceConfig(dict(source=source))
+    magnitudes = Magnitudes(dict(source=6.0))
+    rakes = Rakes(dict(source=180.0))
+
+    velocity_model_parameters = VelocityModelParameters(
+        min_vs=500.0,
+        version="2.09",
+        topo_type="BULLDOZED",
+        ds_multiplier=1.2,
+        vs30=500.0,
+        fault_buffer=2000,
+        s_wave_velocity=3500,
+        rrup_interpolants=np.array([[5.0, 8.0], [50.0, 50.0]]),
+    )
+    domain_parameters = generate_domain.generate_domain(
+        source_config,
+        magnitudes,
+        rakes,
+        velocity_model_parameters,
+    )
+    assert shapely.contains(domain_parameters.domain.polygon, source.geometry)

tests/test_realisation.py

@@ -245,8 +245,9 @@ def test_velocity_model(tmp_path: Path) -> None:
         topo_type="SQUASHED_TAPERED",
         ds_multiplier=1.2,
         vs30=300.0,
+        fault_buffer=2000.0,
         s_wave_velocity=3500.0,
-        pgv_interpolants=np.ones(shape=(2, 2), dtype=np.float32),
+        rrup_interpolants=np.ones(shape=(2, 2), dtype=np.float32),
     )
     realisation_ffp = tmp_path / "realisation.json"
     velocity_model.write_to_realisation(realisation_ffp)
@@ -258,8 +259,9 @@ def test_velocity_model(tmp_path: Path) -> None:
                 "topo_type": "SQUASHED_TAPERED",
                 "ds_multiplier": 1.2,
                 "vs30": 300.0,
+                "fault_buffer": 2000.0,
                 "s_wave_velocity": 3500.0,
-                "pgv_interpolants": [[1, 1], [1, 1]],
+                "rrup_interpolants": [[1, 1], [1, 1]],
             }
         }
 

tests/test_utils.py

@@ -45,6 +45,7 @@ def test_raises_when_no_cpu_info() -> None:
         utils.get_available_cores()
 
 
+@pytest.mark.slow
 def test_read_nz_coastline() -> None:
     gdf = utils.read_nz_coastline()
     assert isinstance(gdf, gpd.GeoDataFrame)
@@ -80,3 +81,64 @@ def test_get_nz_outline_polygon_selects_two_largest() -> None:
     expected_union = shapely.transform(expected_union, lambda x: x[:, ::-1])
     # Check the result is equal to the union of the two largest polygons
     assert shapely.area(shapely.symmetric_difference(result, expected_union)) < 1e-4
+
+
+def test_dict_zip_basic_two_dicts() -> None:
+    """Test standard zipping of two dictionaries with matching keys."""
+    d1 = {"a": 1, "b": 2}
+    d2 = {"a": "apple", "b": "banana"}
+
+    expected = {"a": (1, "apple"), "b": (2, "banana")}
+    assert utils.dict_zip(d1, d2) == expected
+
+
+def test_dict_zip_three_dicts() -> None:
+    """Test standard zipping of three dictionaries."""
+    d1 = {"a": 1}
+    d2 = {"a": 2}
+    d3 = {"a": 3}
+
+    assert utils.dict_zip(d1, d2, d3) == {"a": (1, 2, 3)}
+
+
+def test_dict_zip_strict_mismatch_raises_error() -> None:
+    """Test that strict=True raises ValueError when keys don't match exactly."""
+    d1 = {"a": 1, "b": 2}
+    d2 = {"a": 1}  # Missing 'b'
+
+    with pytest.raises(ValueError, match="Keys in dictionaries are not all the same"):
+        utils.dict_zip(d1, d2, strict=True)
+
+
+def test_dict_zip_non_strict_intersection() -> None:
+    """Test that strict=False returns the intersection of keys."""
+    d1 = {"a": 1, "b": 2, "c": 3}
+    d2 = {"a": 10, "b": 20, "d": 40}
+
+    # Only 'a' and 'b' are in both
+    result = utils.dict_zip(d1, d2, strict=False)
+    assert set(result.keys()) == {"a", "b"}
+    assert result["a"] == (1, 10)
+    assert result["b"] == (2, 20)
+
+
+def test_dict_zip_empty_input() -> None:
+    """Test behaviour with no dictionaries provided."""
+    assert utils.dict_zip() == {}
+
+
+def test_dict_zip_single_dict() -> None:
+    """Test behaviour with a single dictionary."""
+    d1 = {"a": 1, "b": 2}
+    assert utils.dict_zip(d1) == {"a": (1,), "b": (2,)}
+
+
+def test_dict_zip_identical_keys_different_order() -> None:
+    """Test that key order in input doesn't cause strict mode to fail."""
+    d1 = {"a": 1, "b": 2}
+    d2 = {"b": 20, "a": 10}
+
+    # Should not raise ValueError even though insertion order differs
+    result = utils.dict_zip(d1, d2, strict=True)
+    assert result["a"] == (1, 10)
+    assert result["b"] == (2, 20)