From 615ade0a5463e8cae79172179ea7868027ab8399 Mon Sep 17 00:00:00 2001 From: benjamink Date: Wed, 8 Apr 2026 09:26:53 -0700 Subject: [PATCH 01/43] test deployment strategy for leroi air (.exe in .conda-env/Scripts folder with conftest.py added to ensure the PATH variable is accessible in pytest runs). --- tests/plate_simulation/runtest/conftest.py | 33 ++++++++++++++++++++ tests/plate_simulation/runtest/leroi_test.py | 15 +++++++++ 2 files changed, 48 insertions(+) create mode 100644 tests/plate_simulation/runtest/conftest.py create mode 100644 tests/plate_simulation/runtest/leroi_test.py diff --git a/tests/plate_simulation/runtest/conftest.py b/tests/plate_simulation/runtest/conftest.py new file mode 100644 index 00000000..dddedfa3 --- /dev/null +++ b/tests/plate_simulation/runtest/conftest.py @@ -0,0 +1,33 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2026 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +from __future__ import annotations + +import os +import sys +from pathlib import Path + +import pytest + + +@pytest.fixture(autouse=True, scope="session") +def conda_scripts_on_path(): + """Ensure the conda env Scripts directory is on PATH. + PyCharm (and other IDEs) launch pytest with the configured Python interpreter + but do not fully activate the conda environment, so the interpreter's sibling + Scripts/ directory is absent from PATH. This fixture adds it once per session + so that subprocess calls to executables installed there (e.g. LeroiAir550_JR) + resolve correctly without hard-coding any paths in the tests themselves. + """ + scripts_dir = str(Path(sys.executable).parent / "Scripts") + original_path = os.environ.get("PATH", "") + if scripts_dir not in original_path.split(os.pathsep): + os.environ["PATH"] = scripts_dir + os.pathsep + original_path + yield + os.environ["PATH"] = original_path diff --git a/tests/plate_simulation/runtest/leroi_test.py b/tests/plate_simulation/runtest/leroi_test.py new file mode 100644 index 00000000..f15e82f1 --- /dev/null +++ b/tests/plate_simulation/runtest/leroi_test.py @@ -0,0 +1,15 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2026 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +import subprocess + + +def test_leroi_executable(): + subprocess.run("LeroiAir550_JR", check=False) From d763e455ce1def824c991c17631159a178efd219 Mon Sep 17 00:00:00 2001 From: benjamink Date: Wed, 8 Apr 2026 10:00:23 -0700 Subject: [PATCH 02/43] update to run from .bat file with relative pathing --- tests/plate_simulation/runtest/leroi_test.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/tests/plate_simulation/runtest/leroi_test.py b/tests/plate_simulation/runtest/leroi_test.py index f15e82f1..d45d3d9d 100644 --- a/tests/plate_simulation/runtest/leroi_test.py +++ b/tests/plate_simulation/runtest/leroi_test.py @@ -11,5 +11,7 @@ import subprocess -def test_leroi_executable(): - subprocess.run("LeroiAir550_JR", check=False) +def test_leroi_executable(tmp_path): + control_file = tmp_path / "test.cfl" + control_file.touch() + subprocess.run("F2.bat test output", cwd=tmp_path, shell=True, check=False) From 26ea236d7764956f8d77a14e1dd76518d52c932e Mon Sep 17 00:00:00 2001 From: benjamink Date: Wed, 15 Apr 2026 08:49:11 -0700 Subject: [PATCH 03/43] initial structure --- .../uijson/plate_simulation.ui.json | 6 + simpeg_drivers/plate_simulation/driver.py | 59 +++-- .../plate_simulation/leroi_air/__init__.py | 9 + .../plate_simulation/leroi_air/driver.py | 50 ++++ .../plate_simulation/leroi_air/interface.py | 244 ++++++++++++++++++ .../plate_simulation/leroi_air/options.py | 94 +++++++ simpeg_drivers/plate_simulation/options.py | 1 + tests/plate_simulation/leroi_air/__init__.py | 9 + .../plate_simulation/leroi_air/driver_test.py | 12 + tests/plate_simulation/runtest/leroi_test.py | 4 + 10 files changed, 466 insertions(+), 22 deletions(-) create mode 100644 simpeg_drivers/plate_simulation/leroi_air/__init__.py create mode 100644 simpeg_drivers/plate_simulation/leroi_air/driver.py create mode 100644 simpeg_drivers/plate_simulation/leroi_air/interface.py create mode 100644 simpeg_drivers/plate_simulation/leroi_air/options.py create mode 100644 tests/plate_simulation/leroi_air/__init__.py create mode 100644 tests/plate_simulation/leroi_air/driver_test.py diff --git a/simpeg_drivers-assets/uijson/plate_simulation.ui.json b/simpeg_drivers-assets/uijson/plate_simulation.ui.json index a7c5037c..f2b6fcd2 100644 --- a/simpeg_drivers-assets/uijson/plate_simulation.ui.json +++ b/simpeg_drivers-assets/uijson/plate_simulation.ui.json @@ -18,6 +18,12 @@ "enabled": true, "tooltip": "Forward modelling SimPEG Group with at least the topography and survey set" }, + "use_leroi": { + "main": true, + "label": "Use leroi", + "value": false, + "tooltip": "If checked, plate simulation will use LeroiAir to simulate the plate response." + }, "name": { "main": true, "label": "Label", diff --git a/simpeg_drivers/plate_simulation/driver.py b/simpeg_drivers/plate_simulation/driver.py index 575f1348..960edfb0 100644 --- a/simpeg_drivers/plate_simulation/driver.py +++ b/simpeg_drivers/plate_simulation/driver.py @@ -84,28 +84,12 @@ def run(self) -> InversionDriver: @property def simulation_driver(self) -> InversionDriver: if self._simulation_driver is None: - with fetch_active_workspace(self.params.geoh5, mode="r+"): - self.simulation_parameters.mesh = self.mesh - self.simulation_parameters.models.starting_model = self.model - - if not isinstance( - self.simulation_parameters.active_cells.topography_object, - Surface | Points, - ): - raise ValueError( - "The topography object of the forward simulation must be a 'Surface'." - ) - - self.simulation_parameters.out_group = None - driver_class = driver_class_from_name( - self.simulation_parameters.inversion_type, forward_only=True - ) - self._simulation_driver = driver_class( - self.simulation_parameters, - client=self._client, - workers=self._workers, - ) - self._simulation_driver.out_group.parent = self._out_group + run_command = self.simulation_parameters.run_command + is_tem = "time_domain.forward" in run_command + if self.params.use_leroi & is_tem: + self._simulation_driver = self._get_leroi_driver() + else: + self._simulation_driver = self._get_simpeg_driver() return self._simulation_driver @@ -312,6 +296,37 @@ def start_dask_run( """ start_dask_run(cls, json_path, n_workers=n_workers, n_threads=n_threads) + def _get_simpeg_driver(self): + + with fetch_active_workspace(self.params.geoh5, mode="r+"): + self.simulation_parameters.mesh = self.mesh + self.simulation_parameters.models.starting_model = self.model + + if not isinstance( + self.simulation_parameters.active_cells.topography_object, + Surface | Points, + ): + raise ValueError( + "The topography object of the forward simulation must be a 'Surface'." + ) + + self.simulation_parameters.out_group = None + driver_class = driver_class_from_name( + self.simulation_parameters.inversion_type, forward_only=True + ) + self._simulation_driver = driver_class( + self.simulation_parameters, + client=self._client, + workers=self._workers, + ) + self._simulation_driver.out_group.parent = self._out_group + + + def _get_leroi_driver(self): + leroi_opts = LeroiAirOptions.from_plate_simulation(self.params) + driver = LeroiAirDriver(leroi_opts) + return driver + if __name__ == "__main__": file = Path(sys.argv[1]) diff --git a/simpeg_drivers/plate_simulation/leroi_air/__init__.py b/simpeg_drivers/plate_simulation/leroi_air/__init__.py new file mode 100644 index 00000000..bd9d6d2b --- /dev/null +++ b/simpeg_drivers/plate_simulation/leroi_air/__init__.py @@ -0,0 +1,9 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2026 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' diff --git a/simpeg_drivers/plate_simulation/leroi_air/driver.py b/simpeg_drivers/plate_simulation/leroi_air/driver.py new file mode 100644 index 00000000..51c6ea08 --- /dev/null +++ b/simpeg_drivers/plate_simulation/leroi_air/driver.py @@ -0,0 +1,50 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2026 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +from pathlib import Path + +from geoh5py import Workspace + +from simpeg_drivers.plate_simulation.models.options import PlateOptions + +from .interface import LeroiAirInterface +from .options import BackgroundOptions, LeroiAirOptions, ModellingOptions, OutputOptions + + +class LeroiAirDriver: + def __init__(self, options: LeroiAirOptions): + self.options = options + + def run(self): + opts = LeroiAirOptions( + title="test", + background=BackgroundOptions( + basement_thickness=5000, + basement_resistivity=1000, + ), + modelling=ModellingOptions(offtime=3.1, cell_size=10), + output=OutputOptions(channel="all") + ) + + with Workspace("dom_waveform_600Ohmm_bkgr_and_plate.geoh5", mode="r") as geoh5: + survey = geoh5.get_entity("survey")[1] + plate = PlateOptions( + reference=[0.0, 0.0, -20.0], + strike_length=80., + dip_length=100., + thickness=5., + dip_direction=90., + dip=90., + resistivity=1., + ) + interface = LeroiAirInterface(geoh5, survey, plate, opts) + + interface.format_cfl_file() + interface.write_cfl_file(Path('LeroiAir.cfl')) \ No newline at end of file diff --git a/simpeg_drivers/plate_simulation/leroi_air/interface.py b/simpeg_drivers/plate_simulation/leroi_air/interface.py new file mode 100644 index 00000000..8fba7cc2 --- /dev/null +++ b/simpeg_drivers/plate_simulation/leroi_air/interface.py @@ -0,0 +1,244 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2026 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +from __future__ import annotations + +from pathlib import Path +from typing import Any + +import numpy as np +from geoh5py import Workspace +from geoh5py.objects import AirborneFEMReceivers, AirborneTEMReceivers +from numpy import array_split + +from simpeg_drivers.plate_simulation.models.options import PlateOptions + +from .options import LeroiAirOptions + + +class LeroiAirInterface: + """Interface for running LeroiAir from geoh5py objects.""" + + version: str = "8.0" + + def __init__( + self, + geoh5: Workspace, + opts: LeroiAirOptions + ): + self.geoh5 = geoh5 + self.opts = opts + + + + @property + def conductivity_thicknesses(self) -> np.ndarray: + """All conductivity thicknesses.""" + layer_sigma = self.opts.layer_thicknesses * (1 / np.array(self.opts.layer_resistivities)) + plate_sigma = ( + [g.width for g in self.opts.plate_geometries] + * (1 / np.array(self.opts.plate_resistivities)) + ) + return np.hstack([layer_sigma, plate_sigma]) + + @property + def resistivities(self) -> np.ndarray: + """All resistivities.""" + return np.hstack([self.layer_resistivities, self.plate_resistivities]) + + @property + def aliased_values(self) -> dict[str, Any]: + """Serves .cfl input file aliases and corresponding data to line formatter.""" + return { + "TDFD": 1 if self.opts.domain == "time" else 2, + "DO3D": 0 if self.opts.layered_earth_only else 1, + "PRFL": 1, + "ISTOP": 0, + "ISW": 1, + "NSX": len(self.opts.ontime_waveform), + "STEP": 0 if self.opts.modelling.magnetic_field == "dBdt" else 1, + "UNITS": 1, + "NCHNL": len(self.opts.channels), + "KRXW": 2, + "REFTYM": self.opts.timing_mark, + "OFFTIME": self.opts.offtime, + "TXON": self.opts.ontime_waveform[:, 0], + "TXAMP": self.opts.ontime_waveform[:, 1], + "TOPN": self.opts.timing_mark + np.array([0.0] + self.opts.channels[: -1]), + "TCLS": self.opts.timing_mark + np.array(self.opts.channels), + "TMS": self.opts.timing_mark + np.array(self.opts.channels), + "WIDTH": np.array([0.3, 0.3, 0.9]), + "TXCLN": 0., + "CMP": 3, + "KPPM": 0, + "NPPF": 3, + "TXAREA": 1., + "NTRN": 1, + "ZRX0": 0., + "XRX0": 0., + "YRX0": 0., + "NSTAT": len(self.opts.locations), + "SURVEY": 2, + "BAROMTRC": 1, + "LINE_TAG": 0, + "EAST": self.opts.locations[:, 0], + "NORTH": self.opts.locations[:, 1], + "ALT": 13 * np.ones(len(self.opts.locations)), + "NLAYER": 2, + "NPLATE": 1, + "NLITH": 3, + "GND_LVL": 0.0, + "RES": self.resistivities, + "SIG_T": self.conductivity_thicknesses, + "RMU": np.ones_like(self.resistivities), + "REPS": np.ones_like(self.resistivities), + "CHRG": np.zeros_like(self.resistivities), + "CTAU": np.zeros_like(self.resistivities), + "CFREQ": np.ones_like(self.resistivities), + "LITH": np.array([1, 2]) if self.opts.background.has_overburden else np.array([1]), + "LITHP": 3 if self.opts.background.has_overburden else 2, + "THICK": self.layer_properties[:, 1], + "CELLW": self.opts.modelling.cell_size, + "IPLATE": 1, + "CNTR_East": self.plate.reference[0], + "CNTR_North": self.plate.reference[1], + "PLTOP": self.plate.reference[2], + "PLNGTH": self.plate.strike_length, + "DPWDTH": self.plate.dip_length, + "DZM": self.plate.dip_direction, + "DIP": self.plate.dip, + } + + def format_line(self, params: list[str]) -> str: + """format a string from a list of params and the retrieved values.""" + values = [str(self.aliased_values[k]) for k in params] + return f"{" ".join(values)} \t ! {', '.join(params)}" + + def format_line_from_array(self, param: str): + values = [str(k) for k in self.aliased_values[param]] + return f"{" ".join(values)} \t ! {param}" + + def format_multi_line(self, params: str | list[str]) -> str: + """Format a multi-line string a column, or row oriented array.""" + if isinstance(params, str): + values = array_split(self.aliased_values[param], 10) + else: + values = np.column_stack([self.aliased_values[k] for k in params]) + return self._format_multi_line(values) + "\t ! " + ", ".join(params) + + @property + def record_2(self): + return self.format_line(["TDFD", "DO3D", "PRFL", "ISTOP"]) + + @property + def record_3(self): + return self.format_line(["ISW", "NSX", "STEP", "UNITS", "NCHNL", "KRXW", "OFFTIME"]) + + @property + def record_4(self): + return self.format_multi_line(["TXON", "TXAMP"]) + + @property + def record_5(self): + return self.format_line_from_array("TMS") + # return self.format_multi_line(["TOPN", "TCLS"]) + + @property + def record_6(self): + return self.format_line_from_array("WIDTH") + + @property + def record_7(self): + return self.format_line(["TXCLN", "CMP", "KPPM"]) + + @property + def record_7p1(self): + return self.format_line(["NPPF"]) + + @property + def record_7p2(self): + return self.format_line(["TXAREA", "NTRN"]) + + @property + def record_8(self): + return self.format_line(["ZRX0", "XRX0", "YRX0"]) + + @property + def record_9(self): + return self.format_line(["NSTAT", "SURVEY", "BAROMTRC", "LINE_TAG"]) + + @property + def record_9p1(self): + return self.format_multi_line(["EAST", "NORTH", "ALT"]) + + @property + def record_10(self): + return self.format_line(["NLAYER", "NPLATE", "NLITH", "GND_LVL"]) + + @property + def record_11(self): + return self.format_multi_line(["RES", "SIG_T", "RMU", "REPS", "CHRG", "CTAU", "CFREQ"]) + + @property + def record_12(self): + return self.format_multi_line(["LITH", "THICK"]) + + @property + def record_13(self): + return self.format_line(["CELLW"]) + + @property + def record_14(self): + return self.format_line(["LITHP", "CNTR_East", "CNTR_North", "PLTOP"]) + + @property + def record_15(self): + return self.format_line(["PLNGTH", "DPWDTH", "DZM", "DIP"]) + + def format_cfl_file(self): + """ + Generates lines of text for an .cfl input file to run LeroiAir. + + Collects approriate 'Records' and adds lines to the input file one by one. + """ + lines = [] + lines.append(self.opts.title) + lines.append(self.record_2) + lines.append(self.record_3) + lines.append(self.record_4) + lines.append(self.record_5) + lines.append(self.record_6) + lines.append(self.record_7) + if self.aliased_values["KPPM"] > 0: + lines.append(self.record_7p1) + lines.append(self.record_7p2) + lines.append(self.record_8) + lines.append(self.record_9) + lines.append(self.record_9p1) + lines.append(self.record_10) + lines.append(self.record_11) + lines.append(self.record_12) + lines.append(self.record_13) + lines.append(self.record_14) + lines.append(self.record_15) + + return "\n".join(lines) + "\n" + + def write_cfl_file(self, filepath: Path): + with open(filepath, mode='w', encoding="utf-8") as f: + return f.write(self.format_cfl_file()) + + + + def _format_multi_line(self, values: np.ndarray) -> str: + """Format a multi-line string from array.""" + lines = [] + for row in values: + lines.append(f"{' '.join([str(k) for k in row])}") + return "\n".join(lines) \ No newline at end of file diff --git a/simpeg_drivers/plate_simulation/leroi_air/options.py b/simpeg_drivers/plate_simulation/leroi_air/options.py new file mode 100644 index 00000000..e16d42b2 --- /dev/null +++ b/simpeg_drivers/plate_simulation/leroi_air/options.py @@ -0,0 +1,94 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2026 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +from __future__ import annotations + +import numpy as np +from dataclasses import dataclass +from typing import Literal, Self + + +from geoh5py.objects import Points +from geoapps_utils.modelling.plates import PlateModel + +from simpeg_drivers.plate_simulation.options import PlateSimulationOptions + + +@dataclass +class LeroiAirOptions: + survey: Points + layer_resistivities: list[float] + layer_thicknesses: list[float] + plate_resistivities: list[float] + plate_geometries: list[PlateModel] + cell_size: float = 10.0 + magnetic_field: Literal["dBdt", "B"] = "dBdt" + domain: Literal["time", "frequency"] = "time" + layered_earth_only: bool = False + + @classmethod + def from_plate_simulation_options( + cls, + options: PlateSimulationOptions + ) -> Self: + + return cls( + survey = options.simulations_parameters().survey, + layer_resistivities = [ + options.model.overburden_options.overburden_property, + options.background + ], + layer_thicknesses = [ + options.model.overburden_options.thickness, + 9999 + ], + plate_resistivities = [options.model.plate.plate_property], + plate_geometries = [options.model.plate.geometry], + magnetic_field = "dBdt" if "dBdt" in options.data_units else "B" + ) + + @property + def locations(self) -> np.ndarray: + return self.survey.vertices + + @property + def waveform(self) -> np.ndarray: + return self.survey.waveform + + @property + def channels(self) -> np.ndarray: + return self.survey.channels + + @property + def timing_mark(self) -> float: + return self.survey.timing_mark + + @property + def units(self): + return self.survey.units + + @property + def offtime(self) -> float: + half_cycle = 1 / (2 * self.survey.metadata["frequency"]) + ontime = self.waveform[np.where(self.waveform[:, 1] == 0)[0][1], 0] + return half_cycle - ontime + + @property + def ontime_waveform(self) -> np.ndarray: + """slice the waveform for the on times.""" + ontime_waveform = self.opts.waveform[~self._offtime_mask(), :] + endpoint = self.opts.waveform[self._offtime_mask()][0, :] + return np.vstack([ontime_waveform, endpoint]) + + def _offtime_mask(self): + """Returns a mask to slice the offtimes from the waveform array.""" + ind = [bool(np.isclose(k, 0)) for k in self.waveform[:, 1]] + ind[0] = False + return np.array(ind) \ No newline at end of file diff --git a/simpeg_drivers/plate_simulation/options.py b/simpeg_drivers/plate_simulation/options.py index ddbe3d32..8564468a 100644 --- a/simpeg_drivers/plate_simulation/options.py +++ b/simpeg_drivers/plate_simulation/options.py @@ -78,6 +78,7 @@ class PlateSimulationOptions(Options): mesh: MeshOptions model: ModelOptions simulation: SimPEGGroup | UIJsonGroup + use_leroi: bool = False def simulation_parameters(self) -> BaseForwardOptions: """ diff --git a/tests/plate_simulation/leroi_air/__init__.py b/tests/plate_simulation/leroi_air/__init__.py new file mode 100644 index 00000000..bd9d6d2b --- /dev/null +++ b/tests/plate_simulation/leroi_air/__init__.py @@ -0,0 +1,9 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2026 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' diff --git a/tests/plate_simulation/leroi_air/driver_test.py b/tests/plate_simulation/leroi_air/driver_test.py new file mode 100644 index 00000000..029ecfdd --- /dev/null +++ b/tests/plate_simulation/leroi_air/driver_test.py @@ -0,0 +1,12 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2026 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + + +def test_leroi_options(tmp_path): \ No newline at end of file diff --git a/tests/plate_simulation/runtest/leroi_test.py b/tests/plate_simulation/runtest/leroi_test.py index d45d3d9d..3ebc0a4c 100644 --- a/tests/plate_simulation/runtest/leroi_test.py +++ b/tests/plate_simulation/runtest/leroi_test.py @@ -15,3 +15,7 @@ def test_leroi_executable(tmp_path): control_file = tmp_path / "test.cfl" control_file.touch() subprocess.run("F2.bat test output", cwd=tmp_path, shell=True, check=False) + + +def test_leroi_run(tmp_path): + From 43f410e951ca33d087cfd5223ccd0375435facac Mon Sep 17 00:00:00 2001 From: benjamink Date: Wed, 15 Apr 2026 11:36:41 -0700 Subject: [PATCH 04/43] options tested --- .../plate_simulation/leroi_air/interface.py | 97 ++++++++----------- .../plate_simulation/leroi_air/options.py | 84 +++++++++++----- tests/plate_simulation/leroi_air/__init__.py | 39 ++++++++ .../leroi_air/interface_test.py | 9 ++ .../leroi_air/options_test.py | 28 ++++++ 5 files changed, 176 insertions(+), 81 deletions(-) create mode 100644 tests/plate_simulation/leroi_air/interface_test.py create mode 100644 tests/plate_simulation/leroi_air/options_test.py diff --git a/simpeg_drivers/plate_simulation/leroi_air/interface.py b/simpeg_drivers/plate_simulation/leroi_air/interface.py index 8fba7cc2..fecc4777 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/interface.py +++ b/simpeg_drivers/plate_simulation/leroi_air/interface.py @@ -27,31 +27,10 @@ class LeroiAirInterface: version: str = "8.0" - def __init__( - self, - geoh5: Workspace, - opts: LeroiAirOptions - ): + def __init__(self, geoh5: Workspace, opts: LeroiAirOptions): self.geoh5 = geoh5 self.opts = opts - - - @property - def conductivity_thicknesses(self) -> np.ndarray: - """All conductivity thicknesses.""" - layer_sigma = self.opts.layer_thicknesses * (1 / np.array(self.opts.layer_resistivities)) - plate_sigma = ( - [g.width for g in self.opts.plate_geometries] - * (1 / np.array(self.opts.plate_resistivities)) - ) - return np.hstack([layer_sigma, plate_sigma]) - - @property - def resistivities(self) -> np.ndarray: - """All resistivities.""" - return np.hstack([self.layer_resistivities, self.plate_resistivities]) - @property def aliased_values(self) -> dict[str, Any]: """Serves .cfl input file aliases and corresponding data to line formatter.""" @@ -70,19 +49,19 @@ def aliased_values(self) -> dict[str, Any]: "OFFTIME": self.opts.offtime, "TXON": self.opts.ontime_waveform[:, 0], "TXAMP": self.opts.ontime_waveform[:, 1], - "TOPN": self.opts.timing_mark + np.array([0.0] + self.opts.channels[: -1]), + "TOPN": self.opts.timing_mark + np.array([0.0] + self.opts.channels[:-1]), "TCLS": self.opts.timing_mark + np.array(self.opts.channels), "TMS": self.opts.timing_mark + np.array(self.opts.channels), "WIDTH": np.array([0.3, 0.3, 0.9]), - "TXCLN": 0., + "TXCLN": 0.0, "CMP": 3, "KPPM": 0, "NPPF": 3, - "TXAREA": 1., + "TXAREA": 1.0, "NTRN": 1, - "ZRX0": 0., - "XRX0": 0., - "YRX0": 0., + "ZRX0": 0.0, + "XRX0": 0.0, + "YRX0": 0.0, "NSTAT": len(self.opts.locations), "SURVEY": 2, "BAROMTRC": 1, @@ -94,35 +73,35 @@ def aliased_values(self) -> dict[str, Any]: "NPLATE": 1, "NLITH": 3, "GND_LVL": 0.0, - "RES": self.resistivities, - "SIG_T": self.conductivity_thicknesses, - "RMU": np.ones_like(self.resistivities), - "REPS": np.ones_like(self.resistivities), - "CHRG": np.zeros_like(self.resistivities), - "CTAU": np.zeros_like(self.resistivities), - "CFREQ": np.ones_like(self.resistivities), - "LITH": np.array([1, 2]) if self.opts.background.has_overburden else np.array([1]), - "LITHP": 3 if self.opts.background.has_overburden else 2, - "THICK": self.layer_properties[:, 1], - "CELLW": self.opts.modelling.cell_size, + "RES": self.opts.resistivities, + "SIG_T": self.opts.conductivity_thicknesses, + "RMU": np.ones_like(self.opts.resistivities), + "REPS": np.ones_like(self.opts.resistivities), + "CHRG": np.zeros_like(self.opts.resistivities), + "CTAU": np.zeros_like(self.opts.resistivities), + "CFREQ": np.ones_like(self.opts.resistivities), + "LITH": np.array([1, 2]), + "LITHP": 3, + "THICK": self.opts.layer_thicknesses, + "CELLW": self.opts.cell_size, "IPLATE": 1, - "CNTR_East": self.plate.reference[0], - "CNTR_North": self.plate.reference[1], - "PLTOP": self.plate.reference[2], - "PLNGTH": self.plate.strike_length, - "DPWDTH": self.plate.dip_length, - "DZM": self.plate.dip_direction, - "DIP": self.plate.dip, + "CNTR_East": self.opts.plate.reference[0], + "CNTR_North": self.opts.plate.reference[1], + "PLTOP": self.opts.plate.reference[2], + "PLNGTH": self.opts.plate.strike_length, + "DPWDTH": self.opts.plate.dip_length, + "DZM": self.opts.plate.dip_direction, + "DIP": self.opts.plate.dip, } def format_line(self, params: list[str]) -> str: """format a string from a list of params and the retrieved values.""" values = [str(self.aliased_values[k]) for k in params] - return f"{" ".join(values)} \t ! {', '.join(params)}" + return f"{' '.join(values)} \t ! {', '.join(params)}" def format_line_from_array(self, param: str): values = [str(k) for k in self.aliased_values[param]] - return f"{" ".join(values)} \t ! {param}" + return f"{' '.join(values)} \t ! {param}" def format_multi_line(self, params: str | list[str]) -> str: """Format a multi-line string a column, or row oriented array.""" @@ -138,7 +117,9 @@ def record_2(self): @property def record_3(self): - return self.format_line(["ISW", "NSX", "STEP", "UNITS", "NCHNL", "KRXW", "OFFTIME"]) + return self.format_line( + ["ISW", "NSX", "STEP", "UNITS", "NCHNL", "KRXW", "OFFTIME"] + ) @property def record_4(self): @@ -183,7 +164,9 @@ def record_10(self): @property def record_11(self): - return self.format_multi_line(["RES", "SIG_T", "RMU", "REPS", "CHRG", "CTAU", "CFREQ"]) + return self.format_multi_line( + ["RES", "SIG_T", "RMU", "REPS", "CHRG", "CTAU", "CFREQ"] + ) @property def record_12(self): @@ -205,7 +188,7 @@ def format_cfl_file(self): """ Generates lines of text for an .cfl input file to run LeroiAir. - Collects approriate 'Records' and adds lines to the input file one by one. + Collects approriate 'Records' and adds lines to the input file one by one. """ lines = [] lines.append(self.opts.title) @@ -230,15 +213,13 @@ def format_cfl_file(self): return "\n".join(lines) + "\n" - def write_cfl_file(self, filepath: Path): - with open(filepath, mode='w', encoding="utf-8") as f: - return f.write(self.format_cfl_file()) - - - def _format_multi_line(self, values: np.ndarray) -> str: """Format a multi-line string from array.""" lines = [] for row in values: lines.append(f"{' '.join([str(k) for k in row])}") - return "\n".join(lines) \ No newline at end of file + return "\n".join(lines) + + def write_cfl_file(self, filepath: Path): + with open(filepath, mode="w", encoding="utf-8") as f: + return f.write(self.format_cfl_file()) diff --git a/simpeg_drivers/plate_simulation/leroi_air/options.py b/simpeg_drivers/plate_simulation/leroi_air/options.py index e16d42b2..fa4e4a13 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/options.py +++ b/simpeg_drivers/plate_simulation/leroi_air/options.py @@ -10,20 +10,19 @@ from __future__ import annotations -import numpy as np from dataclasses import dataclass from typing import Literal, Self - -from geoh5py.objects import Points +import numpy as np from geoapps_utils.modelling.plates import PlateModel +from geoh5py.objects.surveys.electromagnetics.airborne_tem import AirborneTEMReceivers from simpeg_drivers.plate_simulation.options import PlateSimulationOptions @dataclass class LeroiAirOptions: - survey: Points + survey: AirborneTEMReceivers layer_resistivities: list[float] layer_thicknesses: list[float] plate_resistivities: list[float] @@ -34,61 +33,100 @@ class LeroiAirOptions: layered_earth_only: bool = False @classmethod - def from_plate_simulation_options( - cls, - options: PlateSimulationOptions - ) -> Self: + def from_plate_simulation_options(cls, options: PlateSimulationOptions) -> Self: return cls( - survey = options.simulations_parameters().survey, - layer_resistivities = [ + survey=options.simulations_parameters().survey, + layer_resistivities=[ options.model.overburden_options.overburden_property, - options.background - ], - layer_thicknesses = [ - options.model.overburden_options.thickness, - 9999 + options.background, ], - plate_resistivities = [options.model.plate.plate_property], - plate_geometries = [options.model.plate.geometry], - magnetic_field = "dBdt" if "dBdt" in options.data_units else "B" + layer_thicknesses=[options.model.overburden_options.thickness, 9999], + plate_resistivities=[options.model.plate.plate_property], + plate_geometries=[options.model.plate.geometry], + magnetic_field="dBdt" if "dBdt" in options.data_units else "B", ) @property def locations(self) -> np.ndarray: + """Survey receiver locations.""" return self.survey.vertices @property def waveform(self) -> np.ndarray: + """Survey transmitter waveform.""" return self.survey.waveform + @property + def frequency(self) -> float: + """ + Transmitter frequency. + + Can be set by the user in the metadata. If not set, it will be assumed + that the waveform property contains a full halfcycle and the frequency + will be calculated as the reciprocal of the provided waveform time span. + """ + frequency = self.survey.metadata.get("frequency", None) + if frequency is None: + half_cycle_time = self.waveform[-1, 0] - self.waveform[0, 0] + frequency = 1 / (2 * half_cycle_time) + + return frequency + + # TODO use units to convert time to milliseconds used by leroi. + @property def channels(self) -> np.ndarray: + """Time channel midpoints referenced from the timing_mark.""" return self.survey.channels + @property + def channel_widths(self) -> np.ndarray: + """Time channel widths.""" + channel_widths = self.survey.metadata.get("channel_widths", None) + if channel_widths is None: + channel_widths = np.diff([0] + self.channels) + return channel_widths + @property def timing_mark(self) -> float: + """Reference point for timing of the channels.""" return self.survey.timing_mark @property def units(self): + """Units of the time channels.""" return self.survey.units @property def offtime(self) -> float: - half_cycle = 1 / (2 * self.survey.metadata["frequency"]) + half_cycle = 1 / (2 * self.frequency) ontime = self.waveform[np.where(self.waveform[:, 1] == 0)[0][1], 0] return half_cycle - ontime @property def ontime_waveform(self) -> np.ndarray: - """slice the waveform for the on times.""" - ontime_waveform = self.opts.waveform[~self._offtime_mask(), :] - endpoint = self.opts.waveform[self._offtime_mask()][0, :] + """On-time waveform including leading and trailing 0 current times.""" + ontime_waveform = self.waveform[~self._offtime_mask(), :] + endpoint = self.waveform[self._offtime_mask()][0, :] return np.vstack([ontime_waveform, endpoint]) + @property + def resistivities(self) -> np.ndarray: + """All resistivities.""" + return np.hstack([self.layer_resistivities, self.plate_resistivities]) + + @property + def conductivity_thicknesses(self) -> np.ndarray: + """All conductivity thicknesses.""" + layer_sigma = self.layer_thicknesses * (1 / np.array(self.layer_resistivities)) + plate_sigma = [g.width for g in self.plate_geometries] * ( + 1 / np.array(self.plate_resistivities) + ) + return np.hstack([layer_sigma, plate_sigma]) + def _offtime_mask(self): """Returns a mask to slice the offtimes from the waveform array.""" ind = [bool(np.isclose(k, 0)) for k in self.waveform[:, 1]] ind[0] = False - return np.array(ind) \ No newline at end of file + return np.array(ind) diff --git a/tests/plate_simulation/leroi_air/__init__.py b/tests/plate_simulation/leroi_air/__init__.py index bd9d6d2b..a34ccbe9 100644 --- a/tests/plate_simulation/leroi_air/__init__.py +++ b/tests/plate_simulation/leroi_air/__init__.py @@ -7,3 +7,42 @@ # (see LICENSE file at the root of this source code package). ' # ' # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +import numpy as np +from geoapps_utils.modelling.plates import PlateModel +from geoh5py.shared.utils import fetch_active_workspace + +from simpeg_drivers.plate_simulation.leroi_air.options import LeroiAirOptions +from simpeg_drivers.utils.synthetics.surveys.time_domain.airborne_tdem import ( + generate_airborne_tdem_survey, +) + + +def generate_plate_options(workspace): + x = np.linspace(-1000, 1000, 81) + y = np.linspace(-1000, 1000, 81) + X, Y = np.meshgrid(x, y) + Z = np.zeros_like(X) + + with fetch_active_workspace(workspace) as geoh5: + survey = generate_airborne_tdem_survey(geoh5, X=X, Y=Y, Z=Z) + layer_resistivities = [1500.0, 2000.0, 5000.0] + layer_thicknesses = [50.0, 1000.0, 2000.0] + plate_resistivities = [100.0] + plate_geometries = [ + PlateModel( + strike_length=200, + dip_length=100, + width=10, + direction=90, + dip=45, + ) + ] + opts = LeroiAirOptions( + survey=survey, + layer_resistivities=layer_resistivities, + layer_thicknesses=layer_thicknesses, + plate_resistivities=plate_resistivities, + plate_geometries=plate_geometries, + ) + return opts diff --git a/tests/plate_simulation/leroi_air/interface_test.py b/tests/plate_simulation/leroi_air/interface_test.py new file mode 100644 index 00000000..bd9d6d2b --- /dev/null +++ b/tests/plate_simulation/leroi_air/interface_test.py @@ -0,0 +1,9 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2026 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' diff --git a/tests/plate_simulation/leroi_air/options_test.py b/tests/plate_simulation/leroi_air/options_test.py new file mode 100644 index 00000000..ee593b87 --- /dev/null +++ b/tests/plate_simulation/leroi_air/options_test.py @@ -0,0 +1,28 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2026 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +import numpy as np +from geoapps_utils.modelling.plates import PlateModel +from geoh5py import Workspace + +from . import generate_plate_options + + +def test_options(tmp_path): + + with Workspace(tmp_path / "test.geoh5") as geoh5: + opts = generate_plate_options(workspace=geoh5) + + assert np.allclose(opts.channel_widths, [0.3, 0.3, 0.6]) + assert np.isclose(opts.frequency, 0.142857143) + assert np.isclose(opts.offtime, 1.6) + assert np.allclose(opts.ontime_waveform, opts.waveform[0:8, :]) + assert np.allclose( + opts.conductivity_thicknesses, np.array([0.0333333, 0.5, 0.4, 0.1]) + ) From 1176e1d0b494f5e7513c8884f2a1b589ee3a694e Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 15 Apr 2026 22:10:28 +0000 Subject: [PATCH 05/43] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- simpeg_drivers/plate_simulation/driver.py | 5 ++--- .../plate_simulation/leroi_air/driver.py | 16 ++++++++-------- tests/plate_simulation/leroi_air/driver_test.py | 2 +- tests/plate_simulation/runtest/leroi_test.py | 1 - 4 files changed, 11 insertions(+), 13 deletions(-) diff --git a/simpeg_drivers/plate_simulation/driver.py b/simpeg_drivers/plate_simulation/driver.py index 960edfb0..1ff12169 100644 --- a/simpeg_drivers/plate_simulation/driver.py +++ b/simpeg_drivers/plate_simulation/driver.py @@ -303,8 +303,8 @@ def _get_simpeg_driver(self): self.simulation_parameters.models.starting_model = self.model if not isinstance( - self.simulation_parameters.active_cells.topography_object, - Surface | Points, + self.simulation_parameters.active_cells.topography_object, + Surface | Points, ): raise ValueError( "The topography object of the forward simulation must be a 'Surface'." @@ -321,7 +321,6 @@ def _get_simpeg_driver(self): ) self._simulation_driver.out_group.parent = self._out_group - def _get_leroi_driver(self): leroi_opts = LeroiAirOptions.from_plate_simulation(self.params) driver = LeroiAirDriver(leroi_opts) diff --git a/simpeg_drivers/plate_simulation/leroi_air/driver.py b/simpeg_drivers/plate_simulation/leroi_air/driver.py index 51c6ea08..e5757784 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/driver.py +++ b/simpeg_drivers/plate_simulation/leroi_air/driver.py @@ -30,21 +30,21 @@ def run(self): basement_resistivity=1000, ), modelling=ModellingOptions(offtime=3.1, cell_size=10), - output=OutputOptions(channel="all") + output=OutputOptions(channel="all"), ) with Workspace("dom_waveform_600Ohmm_bkgr_and_plate.geoh5", mode="r") as geoh5: survey = geoh5.get_entity("survey")[1] plate = PlateOptions( reference=[0.0, 0.0, -20.0], - strike_length=80., - dip_length=100., - thickness=5., - dip_direction=90., - dip=90., - resistivity=1., + strike_length=80.0, + dip_length=100.0, + thickness=5.0, + dip_direction=90.0, + dip=90.0, + resistivity=1.0, ) interface = LeroiAirInterface(geoh5, survey, plate, opts) interface.format_cfl_file() - interface.write_cfl_file(Path('LeroiAir.cfl')) \ No newline at end of file + interface.write_cfl_file(Path("LeroiAir.cfl")) diff --git a/tests/plate_simulation/leroi_air/driver_test.py b/tests/plate_simulation/leroi_air/driver_test.py index 029ecfdd..c4fbffea 100644 --- a/tests/plate_simulation/leroi_air/driver_test.py +++ b/tests/plate_simulation/leroi_air/driver_test.py @@ -9,4 +9,4 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -def test_leroi_options(tmp_path): \ No newline at end of file +def test_leroi_options(tmp_path): diff --git a/tests/plate_simulation/runtest/leroi_test.py b/tests/plate_simulation/runtest/leroi_test.py index 3ebc0a4c..ef1041ea 100644 --- a/tests/plate_simulation/runtest/leroi_test.py +++ b/tests/plate_simulation/runtest/leroi_test.py @@ -18,4 +18,3 @@ def test_leroi_executable(tmp_path): def test_leroi_run(tmp_path): - From 77da1d54ee923be777b45088c772668e15b913cd Mon Sep 17 00:00:00 2001 From: benjamink Date: Thu, 16 Apr 2026 14:24:05 -0700 Subject: [PATCH 06/43] Filled in unit tests, working on runtest --- .pre-commit-config.yaml | 10 +- simpeg_drivers/plate_simulation/driver.py | 10 +- .../plate_simulation/leroi_air/driver.py | 58 ++--- .../plate_simulation/leroi_air/interface.py | 208 +++++++++++++----- .../plate_simulation/leroi_air/options.py | 34 ++- tests/plate_simulation/leroi_air/conftest.py | 74 +++++++ .../plate_simulation/leroi_air/driver_test.py | 31 ++- .../leroi_air/interface_test.py | 177 +++++++++++++++ tests/plate_simulation/runtest/leroi_test.py | 96 ++++++++ 9 files changed, 600 insertions(+), 98 deletions(-) create mode 100644 tests/plate_simulation/leroi_air/conftest.py diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 7aecb968..834c406d 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -65,7 +65,15 @@ repos: language: system require_serial: true # pylint does its own parallelism types: [python] - exclude: ^(devtools|docs)/ + exclude: ^(devtools|docs|tests)/ + - id: pylint-tests + name: pylint (tests) + entry: .\\devtools\\conda_env_pylint.bat + language: system + require_serial: true + args: [--rcfile=pylintrc, --disable=protected-access] + types: [python] + files: ^tests/ - repo: https://github.com/codespell-project/codespell rev: v2.4.1 hooks: diff --git a/simpeg_drivers/plate_simulation/driver.py b/simpeg_drivers/plate_simulation/driver.py index 960edfb0..1fb15868 100644 --- a/simpeg_drivers/plate_simulation/driver.py +++ b/simpeg_drivers/plate_simulation/driver.py @@ -29,6 +29,8 @@ validate_workers, ) from simpeg_drivers.options import BaseForwardOptions, ModelTypeEnum +from simpeg_drivers.plate_simulation.leroi_air.driver import LeroiAirDriver +from simpeg_drivers.plate_simulation.leroi_air.options import LeroiAirOptions from simpeg_drivers.plate_simulation.models.events import Anomaly, Erosion, Overburden from simpeg_drivers.plate_simulation.models.series import DikeSwarm, Geology from simpeg_drivers.plate_simulation.options import PlateSimulationOptions @@ -303,8 +305,8 @@ def _get_simpeg_driver(self): self.simulation_parameters.models.starting_model = self.model if not isinstance( - self.simulation_parameters.active_cells.topography_object, - Surface | Points, + self.simulation_parameters.active_cells.topography_object, + Surface | Points, ): raise ValueError( "The topography object of the forward simulation must be a 'Surface'." @@ -321,10 +323,12 @@ def _get_simpeg_driver(self): ) self._simulation_driver.out_group.parent = self._out_group + return self.simulation_driver def _get_leroi_driver(self): - leroi_opts = LeroiAirOptions.from_plate_simulation(self.params) + leroi_opts = LeroiAirOptions.from_plate_simulation_options(self.params) driver = LeroiAirDriver(leroi_opts) + driver.out_group = self._out_group return driver diff --git a/simpeg_drivers/plate_simulation/leroi_air/driver.py b/simpeg_drivers/plate_simulation/leroi_air/driver.py index 51c6ea08..9ee051f5 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/driver.py +++ b/simpeg_drivers/plate_simulation/leroi_air/driver.py @@ -8,43 +8,51 @@ # ' # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +import subprocess from pathlib import Path -from geoh5py import Workspace - -from simpeg_drivers.plate_simulation.models.options import PlateOptions +from geoh5py.groups import UIJsonGroup from .interface import LeroiAirInterface -from .options import BackgroundOptions, LeroiAirOptions, ModellingOptions, OutputOptions +from .options import LeroiAirOptions class LeroiAirDriver: def __init__(self, options: LeroiAirOptions): self.options = options + self._interface: LeroiAirInterface | None = None + self.out_group: UIJsonGroup | None = None + + @property + def interface(self) -> LeroiAirInterface: + if self._interface is None: + self._interface = LeroiAirInterface(self.options) + return self._interface + + @property + def project_path(self) -> Path: + return self.options.survey.workspace.h5file.parent def run(self): - opts = LeroiAirOptions( - title="test", - background=BackgroundOptions( - basement_thickness=5000, - basement_resistivity=1000, - ), - modelling=ModellingOptions(offtime=3.1, cell_size=10), - output=OutputOptions(channel="all") + self.interface.write_cfl_file(self.project_path / "LeroiAir.cfl") + + result = subprocess.run( + ["LeroiAir550_JR", "LeroiAir"], + cwd=self.project_path, + capture_output=True, + text=True, + check=False, ) - with Workspace("dom_waveform_600Ohmm_bkgr_and_plate.geoh5", mode="r") as geoh5: - survey = geoh5.get_entity("survey")[1] - plate = PlateOptions( - reference=[0.0, 0.0, -20.0], - strike_length=80., - dip_length=100., - thickness=5., - dip_direction=90., - dip=90., - resistivity=1., + if result.returncode != 0: + raise RuntimeError( + f"LeroiAir failed with return code {result.returncode}.\n" + f"stderr:\n{result.stderr}\n" + f"stdout:\n{result.stdout}" ) - interface = LeroiAirInterface(geoh5, survey, plate, opts) - interface.format_cfl_file() - interface.write_cfl_file(Path('LeroiAir.cfl')) \ No newline at end of file + outfile = self.project_path / "LeroiAir.out" + self.interface.save_to_geoh5( + outfile=outfile, + out_group=self.out_group, + ) diff --git a/simpeg_drivers/plate_simulation/leroi_air/interface.py b/simpeg_drivers/plate_simulation/leroi_air/interface.py index fecc4777..c9dc5e2d 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/interface.py +++ b/simpeg_drivers/plate_simulation/leroi_air/interface.py @@ -10,15 +10,15 @@ from __future__ import annotations from pathlib import Path -from typing import Any +from typing import Any, Literal import numpy as np -from geoh5py import Workspace -from geoh5py.objects import AirborneFEMReceivers, AirborneTEMReceivers +from geoh5py.groups import UIJsonGroup +from geoh5py.objects import MaxwellPlate +from geoh5py.objects.maxwell_plate import PlateGeometry, PlatePosition +from geoh5py.shared.utils import fetch_active_workspace from numpy import array_split -from simpeg_drivers.plate_simulation.models.options import PlateOptions - from .options import LeroiAirOptions @@ -27,9 +27,9 @@ class LeroiAirInterface: version: str = "8.0" - def __init__(self, geoh5: Workspace, opts: LeroiAirOptions): - self.geoh5 = geoh5 + def __init__(self, opts: LeroiAirOptions): self.opts = opts + self.float_precision: int = 4 @property def aliased_values(self) -> dict[str, Any]: @@ -41,7 +41,7 @@ def aliased_values(self) -> dict[str, Any]: "ISTOP": 0, "ISW": 1, "NSX": len(self.opts.ontime_waveform), - "STEP": 0 if self.opts.modelling.magnetic_field == "dBdt" else 1, + "STEP": 0 if self.opts.magnetic_field == "dBdt" else 1, "UNITS": 1, "NCHNL": len(self.opts.channels), "KRXW": 2, @@ -49,10 +49,8 @@ def aliased_values(self) -> dict[str, Any]: "OFFTIME": self.opts.offtime, "TXON": self.opts.ontime_waveform[:, 0], "TXAMP": self.opts.ontime_waveform[:, 1], - "TOPN": self.opts.timing_mark + np.array([0.0] + self.opts.channels[:-1]), - "TCLS": self.opts.timing_mark + np.array(self.opts.channels), "TMS": self.opts.timing_mark + np.array(self.opts.channels), - "WIDTH": np.array([0.3, 0.3, 0.9]), + "WIDTH": self.opts.channel_widths, "TXCLN": 0.0, "CMP": 3, "KPPM": 0, @@ -62,16 +60,16 @@ def aliased_values(self) -> dict[str, Any]: "ZRX0": 0.0, "XRX0": 0.0, "YRX0": 0.0, - "NSTAT": len(self.opts.locations), + "NSTAT": self.opts.n_stations, "SURVEY": 2, "BAROMTRC": 1, "LINE_TAG": 0, "EAST": self.opts.locations[:, 0], "NORTH": self.opts.locations[:, 1], - "ALT": 13 * np.ones(len(self.opts.locations)), - "NLAYER": 2, - "NPLATE": 1, - "NLITH": 3, + "ALT": 13 * np.ones(self.opts.n_stations), # TODO - take from survey + "NLAYER": self.opts.n_layers, + "NPLATE": self.opts.n_plates, + "NLITH": self.opts.n_layers + self.opts.n_plates, "GND_LVL": 0.0, "RES": self.opts.resistivities, "SIG_T": self.opts.conductivity_thicknesses, @@ -80,115 +78,137 @@ def aliased_values(self) -> dict[str, Any]: "CHRG": np.zeros_like(self.opts.resistivities), "CTAU": np.zeros_like(self.opts.resistivities), "CFREQ": np.ones_like(self.opts.resistivities), - "LITH": np.array([1, 2]), - "LITHP": 3, + "LITH": 1 + np.arange(self.opts.n_layers, dtype=int), + "LITHP": 1 + np.arange(self.opts.n_plates, dtype=int) + self.opts.n_layers, "THICK": self.opts.layer_thicknesses, "CELLW": self.opts.cell_size, "IPLATE": 1, - "CNTR_East": self.opts.plate.reference[0], - "CNTR_North": self.opts.plate.reference[1], - "PLTOP": self.opts.plate.reference[2], - "PLNGTH": self.opts.plate.strike_length, - "DPWDTH": self.opts.plate.dip_length, - "DZM": self.opts.plate.dip_direction, - "DIP": self.opts.plate.dip, + "CNTR_East": [g.easting for g in self.opts.plate_geometries], + "CNTR_North": [g.northing for g in self.opts.plate_geometries], + "PLTOP": [g.elevation for g in self.opts.plate_geometries], + "PLNGTH": [g.strike_length for g in self.opts.plate_geometries], + "DPWDTH": [g.dip_length for g in self.opts.plate_geometries], + "DZM": [g.direction for g in self.opts.plate_geometries], + "DIP": [g.dip for g in self.opts.plate_geometries], } + def _format_value(self, value: int | float) -> str: + """Format a scalar value, preserving integer formatting and applying float precision only when needed.""" + match value: + case int() | np.integer(): + return str(int(value)) + case float() | np.floating(): + return self._format_float(value) + case _: + return str(value) + + def _format_float(self, value: float) -> str: + """Format a float, truncating to float_precision only when needed.""" + _, _, decimals = str(value).partition(".") + if len(decimals) > self.float_precision: + return f"{value:.{self.float_precision}f}" + return str(value) + def format_line(self, params: list[str]) -> str: """format a string from a list of params and the retrieved values.""" - values = [str(self.aliased_values[k]) for k in params] + values = [self._format_value(self.aliased_values[k]) for k in params] return f"{' '.join(values)} \t ! {', '.join(params)}" def format_line_from_array(self, param: str): - values = [str(k) for k in self.aliased_values[param]] + """Format a line string from an array.""" + values = [self._format_value(k) for k in self.aliased_values[param]] return f"{' '.join(values)} \t ! {param}" def format_multi_line(self, params: str | list[str]) -> str: - """Format a multi-line string a column, or row oriented array.""" + """Format a multi-line string from a column, or row oriented array.""" if isinstance(params, str): - values = array_split(self.aliased_values[param], 10) + rows = [ + [v] + for chunk in array_split(self.aliased_values[params], 10) + for v in chunk + ] else: - values = np.column_stack([self.aliased_values[k] for k in params]) - return self._format_multi_line(values) + "\t ! " + ", ".join(params) + columns = [self.aliased_values[k] for k in params] + rows = [list(row) for row in zip(*columns, strict=True)] + return self._format_rows(rows) + "\t ! " + ", ".join(params) @property - def record_2(self): + def record_2(self) -> str: return self.format_line(["TDFD", "DO3D", "PRFL", "ISTOP"]) @property - def record_3(self): + def record_3(self) -> str: return self.format_line( ["ISW", "NSX", "STEP", "UNITS", "NCHNL", "KRXW", "OFFTIME"] ) @property - def record_4(self): + def record_4(self) -> str: return self.format_multi_line(["TXON", "TXAMP"]) @property - def record_5(self): + def record_5(self) -> str: return self.format_line_from_array("TMS") - # return self.format_multi_line(["TOPN", "TCLS"]) @property - def record_6(self): + def record_6(self) -> str: return self.format_line_from_array("WIDTH") @property - def record_7(self): + def record_7(self) -> str: return self.format_line(["TXCLN", "CMP", "KPPM"]) @property - def record_7p1(self): + def record_7p1(self) -> str: return self.format_line(["NPPF"]) @property - def record_7p2(self): + def record_7p2(self) -> str: return self.format_line(["TXAREA", "NTRN"]) @property - def record_8(self): + def record_8(self) -> str: return self.format_line(["ZRX0", "XRX0", "YRX0"]) @property - def record_9(self): + def record_9(self) -> str: return self.format_line(["NSTAT", "SURVEY", "BAROMTRC", "LINE_TAG"]) @property - def record_9p1(self): + def record_9p1(self) -> str: return self.format_multi_line(["EAST", "NORTH", "ALT"]) @property - def record_10(self): + def record_10(self) -> str: return self.format_line(["NLAYER", "NPLATE", "NLITH", "GND_LVL"]) @property - def record_11(self): + def record_11(self) -> str: return self.format_multi_line( ["RES", "SIG_T", "RMU", "REPS", "CHRG", "CTAU", "CFREQ"] ) @property - def record_12(self): + def record_12(self) -> str: return self.format_multi_line(["LITH", "THICK"]) @property - def record_13(self): + def record_13(self) -> str: return self.format_line(["CELLW"]) @property - def record_14(self): - return self.format_line(["LITHP", "CNTR_East", "CNTR_North", "PLTOP"]) + def record_14(self) -> str: + return self.format_multi_line(["LITHP", "CNTR_East", "CNTR_North", "PLTOP"]) @property - def record_15(self): - return self.format_line(["PLNGTH", "DPWDTH", "DZM", "DIP"]) + def record_15(self) -> str: + return self.format_multi_line(["PLNGTH", "DPWDTH", "DZM", "DIP"]) - def format_cfl_file(self): + def format_cfl_file(self) -> str: """ Generates lines of text for an .cfl input file to run LeroiAir. - Collects approriate 'Records' and adds lines to the input file one by one. + Collects appropriate 'Records' and adds lines to the input file one by one. """ lines = [] lines.append(self.opts.title) @@ -213,13 +233,79 @@ def format_cfl_file(self): return "\n".join(lines) + "\n" - def _format_multi_line(self, values: np.ndarray) -> str: - """Format a multi-line string from array.""" - lines = [] - for row in values: - lines.append(f"{' '.join([str(k) for k in row])}") - return "\n".join(lines) + def _format_rows(self, rows: list[list]) -> str: + """Format a multi-line string from a list of rows.""" + return "\n".join(" ".join(self._format_value(v) for v in row) for row in rows) - def write_cfl_file(self, filepath: Path): + def write_cfl_file(self, filepath: Path) -> None: + """Write the formatted .cfl input file to disk.""" with open(filepath, mode="w", encoding="utf-8") as f: - return f.write(self.format_cfl_file()) + f.write(self.format_cfl_file()) + + _COMPONENT_ANCHORS: dict[str, str] = { + "x": "TRANSVERSE COMPONENT", + "y": "IN-LINE COMPONENT", + "z": "VERTICAL COMPONENT", + } + + def _slice_data_lines(self, lines: list[str], anchor: str) -> list[str]: + """Slice the station data rows that follow the given section header.""" + anchor_idx = next(i for i, line in enumerate(lines) if anchor in line) + chunk = lines[anchor_idx:] + data_start = ( + next( + i + for i, line in enumerate(chunk) + if all(k in line for k in ["EAST", "NORTH", "ALT"]) + ) + + 2 + ) + return chunk[data_start : data_start + self.opts.n_stations] + + def _extract_data( + self, outfile: str | Path, component: Literal["x", "y", "z"] + ) -> np.ndarray: + """Extract channel data for a single component from a LeroiAir .out file.""" + lines = Path(outfile).read_text(encoding="utf-8", errors="replace").splitlines() + data_lines = self._slice_data_lines(lines, self._COMPONENT_ANCHORS[component]) + return np.array([line.split() for line in data_lines], dtype=float)[:, 4:] + + def save_to_geoh5(self, outfile: str | Path, out_group: UIJsonGroup): + """Save LeroiAir simulated data on the provided survey to geoh5.""" + + crossline_data = self._extract_data(outfile=outfile, component="x") + inline_data = self._extract_data(outfile=outfile, component="y") + vertical_data = self._extract_data(outfile=outfile, component="z") + + with fetch_active_workspace(self.opts.survey.workspace, mode="r+") as geoh5: + survey = self.opts.survey.copy(parent=out_group, copy_children=False) + data = survey.add_data( + { + "fwd inline [0]": {"values": inline_data[:, 0]}, + "fwd inline [1]": {"values": inline_data[:, 1]}, + "fwd inline [2]": {"values": inline_data[:, 2]}, + "fwd crossline [0]": {"values": crossline_data[:, 0]}, + "fwd crossline [1]": {"values": crossline_data[:, 1]}, + "fwd crossline [2]": {"values": crossline_data[:, 2]}, + "fwd vertical [0]": {"values": vertical_data[:, 0]}, + "fwd vertical [1]": {"values": vertical_data[:, 1]}, + "fwd vertical [2]": {"values": vertical_data[:, 2]}, + } + ) + survey.create_property_group(name="inline", properties=data[:3]) + survey.create_property_group(name="crossline", properties=data[3:6]) + survey.create_property_group(name="vertical", properties=data[6:]) + + for plate in self.opts.plate_geometries: + position = PlatePosition( + x=plate.easting, y=plate.northing, z=plate.elevation + ) + geometry = PlateGeometry( + position=position, + length=plate.strike_length, + width=plate.dip_length, + thickness=plate.width, + dip_direction=plate.direction, + dip=plate.dip, + ) + MaxwellPlate.create(geoh5, geometry=geometry) diff --git a/simpeg_drivers/plate_simulation/leroi_air/options.py b/simpeg_drivers/plate_simulation/leroi_air/options.py index fa4e4a13..529cdc3b 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/options.py +++ b/simpeg_drivers/plate_simulation/leroi_air/options.py @@ -34,24 +34,44 @@ class LeroiAirOptions: @classmethod def from_plate_simulation_options(cls, options: PlateSimulationOptions) -> Self: - + simulation_options = options.simulation_parameters() return cls( - survey=options.simulations_parameters().survey, + survey=simulation_options.data_object, layer_resistivities=[ options.model.overburden_options.overburden_property, - options.background, + options.model.background, ], layer_thicknesses=[options.model.overburden_options.thickness, 9999], - plate_resistivities=[options.model.plate.plate_property], - plate_geometries=[options.model.plate.geometry], - magnetic_field="dBdt" if "dBdt" in options.data_units else "B", + plate_resistivities=[options.model.plate_options.plate_property], + plate_geometries=[options.model.plate_options.geometry], + magnetic_field="dBdt" if "dBdt" in simulation_options.data_units else "B", ) + @property + def title(self) -> str: + """Provides a generic title for all LeroiAir simulations.""" + return "LeroiAir modelling for plate-simulation package." + @property def locations(self) -> np.ndarray: """Survey receiver locations.""" return self.survey.vertices + @property + def n_stations(self) -> int: + """Number of survey stations at which time channel data will be simulated""" + return len(self.locations) + + @property + def n_layers(self) -> int: + """Number of background layers.""" + return len(self.layer_resistivities) + + @property + def n_plates(self) -> int: + """Number of plates.""" + return len(self.plate_geometries) + @property def waveform(self) -> np.ndarray: """Survey transmitter waveform.""" @@ -96,7 +116,7 @@ def timing_mark(self) -> float: @property def units(self): """Units of the time channels.""" - return self.survey.units + return self.survey.unit @property def offtime(self) -> float: diff --git a/tests/plate_simulation/leroi_air/conftest.py b/tests/plate_simulation/leroi_air/conftest.py new file mode 100644 index 00000000..0fe8d76e --- /dev/null +++ b/tests/plate_simulation/leroi_air/conftest.py @@ -0,0 +1,74 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2026 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +from unittest.mock import MagicMock + +import pytest +from geoh5py import Workspace + +from simpeg_drivers.plate_simulation.leroi_air.interface import LeroiAirInterface + +from . import generate_plate_options + + +N_CHANNELS = 3 +N_STATIONS = 2 + +FAKE_OUT = """\ + VERTICAL COMPONENT - nT/s + + TRANSMITTER POSITION CHNL 1 CHNL 2 CHNL 3 + EAST NORTH ALT 2.300 2.600 3.200 + + 1 -100 0 13 13.3 14.4 15.5 + 2 -87 0 13 16.6 17.7 18.8 +------------------------------------------------------------------------------------- + + + IN-LINE COMPONENT - nT/s + + TRANSMITTER POSITION CHNL 1 CHNL 2 CHNL 3 + EAST NORTH ALT 2.300 2.600 3.200 + + 1 -100 0 13 7.7 8.8 9.9 + 2 -87 0 13 10.0 11.1 12.2 +------------------------------------------------------------------------------------- + + + TRANSVERSE COMPONENT - nT/s + + TRANSMITTER POSITION CHNL 1 CHNL 2 CHNL 3 + EAST NORTH ALT 2.300 2.600 3.200 + + 1 -100 0 13 1.1 2.2 3.3 + 2 -87 0 13 4.4 5.5 6.6 +------------------------------------------------------------------------------------- +""" + + +@pytest.fixture +def mock_interface(): + opts = MagicMock() + opts.n_stations = N_STATIONS + return LeroiAirInterface(opts=opts) + + +@pytest.fixture +def real_interface(tmp_path): + with Workspace(tmp_path / "test.geoh5") as geoh5: + opts = generate_plate_options(geoh5) + return LeroiAirInterface(opts=opts) + + +@pytest.fixture +def fake_outfile(tmp_path): + path = tmp_path / "leroi.out" + path.write_text(FAKE_OUT, encoding="utf-8") + return path diff --git a/tests/plate_simulation/leroi_air/driver_test.py b/tests/plate_simulation/leroi_air/driver_test.py index 029ecfdd..bb7efa94 100644 --- a/tests/plate_simulation/leroi_air/driver_test.py +++ b/tests/plate_simulation/leroi_air/driver_test.py @@ -8,5 +8,34 @@ # ' # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +import subprocess +from unittest.mock import MagicMock, patch -def test_leroi_options(tmp_path): \ No newline at end of file +import pytest +from geoh5py import Workspace + +from simpeg_drivers.plate_simulation.leroi_air.driver import LeroiAirDriver + +from . import generate_plate_options + + +def test_leroi_air_driver_run_raises_on_subprocess_failure(tmp_path): + with Workspace(tmp_path / "test.geoh5") as geoh5: + opts = generate_plate_options(workspace=geoh5) + + driver = LeroiAirDriver(opts) + driver._interface = MagicMock() + failed_result = MagicMock( + spec=subprocess.CompletedProcess, returncode=1, stderr="err", stdout="out" + ) + + with patch( + "simpeg_drivers.plate_simulation.leroi_air.driver.subprocess.run", + return_value=failed_result, + ): + with pytest.raises(RuntimeError) as exc_info: + driver.run() + + assert "return code 1" in str(exc_info.value) + assert "err" in str(exc_info.value) + assert "out" in str(exc_info.value) diff --git a/tests/plate_simulation/leroi_air/interface_test.py b/tests/plate_simulation/leroi_air/interface_test.py index bd9d6d2b..203bec3f 100644 --- a/tests/plate_simulation/leroi_air/interface_test.py +++ b/tests/plate_simulation/leroi_air/interface_test.py @@ -7,3 +7,180 @@ # (see LICENSE file at the root of this source code package). ' # ' # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +import numpy as np +import pytest +from geoh5py import Workspace + +from simpeg_drivers.plate_simulation.leroi_air.interface import LeroiAirInterface + +from . import generate_plate_options +from .conftest import FAKE_OUT, N_CHANNELS, N_STATIONS + + +def test_line_formatting(tmp_path): + with Workspace(tmp_path / "test.geoh5") as geoh5: + opts = generate_plate_options(geoh5) + + interface = LeroiAirInterface(opts=opts) + + line = interface.format_line(["TXCLN", "CMP", "KPPM"]) + assert line == "0.0 3 0 \t ! TXCLN, CMP, KPPM" + + line = interface.format_line_from_array("TMS") + assert line == "2.3 2.6 3.2 \t ! TMS" + + multi_line = interface.format_multi_line(["TMS", "WIDTH"]) + assert multi_line == "2.3 0.3\n2.6 0.3\n3.2 0.6\t ! TMS, WIDTH" + + +def test_format_cfl_file(tmp_path): + with Workspace(tmp_path / "test.geoh5") as geoh5: + opts = generate_plate_options(geoh5) + + interface = LeroiAirInterface(opts=opts) + cfl_str = interface.format_cfl_file() + + # title + assert opts.title in cfl_str + + # record_2: TDFD, DO3D, PRFL, ISTOP + assert "1 1 1 0 \t ! TDFD, DO3D, PRFL, ISTOP" in cfl_str + + # record_3: ISW, NSX, STEP, UNITS, NCHNL, KRXW, OFFTIME + assert "1 8 0 1 3 2 1.6 \t ! ISW, NSX, STEP, UNITS, NCHNL, KRXW, OFFTIME" in cfl_str + + # record_4: TXON, TXAMP (multi-line waveform table) + assert ( + "0.0 0.0\n" + "0.5 0.3333\n" + "1.0 0.6667\n" + "1.5 1.0\n" + "1.6 0.9\n" + "1.7 0.6000\n" + "1.8 0.3000\n" + "1.9 0.0\t ! TXON, TXAMP" + ) in cfl_str + + # record_5: TMS (time channel midpoints) + assert "2.3 2.6 3.2 \t ! TMS" in cfl_str + + # record_6: WIDTH (time channel widths) + assert "0.3 0.3 0.6 \t ! WIDTH" in cfl_str + + # record_7: TXCLN, CMP, KPPM + assert "0.0 3 0 \t ! TXCLN, CMP, KPPM" in cfl_str + + # record_7p1 (NPPF) is only written when KPPM > 0; default KPPM=0 so absent + assert "3 \t ! NPPF" not in cfl_str + + # record_7p2: TXAREA, NTRN + assert "1.0 1 \t ! TXAREA, NTRN" in cfl_str + + # record_8: ZRX0, XRX0, YRX0 + assert "0.0 0.0 0.0 \t ! ZRX0, XRX0, YRX0" in cfl_str + + # record_9: NSTAT, SURVEY, BAROMTRC, LINE_TAG + assert "6561 2 1 0 \t ! NSTAT, SURVEY, BAROMTRC, LINE_TAG" in cfl_str + + # record_9p1: EAST, NORTH, ALT — too long to assert explicitly + + # record_10: NLAYER, NPLATE, NLITH, GND_LVL + assert "3 1 4 0.0 \t ! NLAYER, NPLATE, NLITH, GND_LVL" in cfl_str + + # record_11: RES, SIG_T, RMU, REPS, CHRG, CTAU, CFREQ (multi-line lithology table) + assert ( + "1500.0 0.0333 1.0 1.0 0.0 0.0 1.0\n" + "2000.0 0.5 1.0 1.0 0.0 0.0 1.0\n" + "5000.0 0.4 1.0 1.0 0.0 0.0 1.0\n" + "100.0 0.1 1.0 1.0 0.0 0.0 1.0\t ! RES, SIG_T, RMU, REPS, CHRG, CTAU, CFREQ" + ) in cfl_str + + # record_12: LITH, THICK (multi-line layer table) + assert ("1 50.0\n2 1000.0\n3 2000.0\t ! LITH, THICK") in cfl_str + + # record_13: CELLW + assert "10.0 \t ! CELLW" in cfl_str + + # record_14: LITHP, CNTR_East, CNTR_North, PLTOP + assert "4 0.0 0.0 0.0\t ! LITHP, CNTR_East, CNTR_North, PLTOP" in cfl_str + + # record_15: PLNGTH, DPWDTH, DZM, DIP + assert "200.0 100.0 90.0 45.0\t ! PLNGTH, DPWDTH, DZM, DIP" in cfl_str + + +def test_slice_data_lines_returns_correct_station_rows(mock_interface): + lines = FAKE_OUT.splitlines() + anchor = LeroiAirInterface._COMPONENT_ANCHORS["x"] + rows = mock_interface._slice_data_lines(lines, anchor) + assert len(rows) == N_STATIONS + assert rows[0].split()[4] == "1.1" + assert rows[1].split()[4] == "4.4" + + +def test_extract_data_returns_channel_columns_only(mock_interface, fake_outfile): + data = mock_interface._extract_data(fake_outfile, component="x") + assert data.shape == (N_STATIONS, N_CHANNELS) + + +def test_extract_data_parses_transverse_component_values(mock_interface, fake_outfile): + data = mock_interface._extract_data(fake_outfile, component="x") + np.testing.assert_array_almost_equal(data[0], [1.1, 2.2, 3.3]) + np.testing.assert_array_almost_equal(data[1], [4.4, 5.5, 6.6]) + + +def test_extract_data_parses_inline_component_values(mock_interface, fake_outfile): + data = mock_interface._extract_data(fake_outfile, component="y") + np.testing.assert_array_almost_equal(data[0], [7.7, 8.8, 9.9]) + np.testing.assert_array_almost_equal(data[1], [10.0, 11.1, 12.2]) + + +def test_extract_data_parses_vertical_component_values(mock_interface, fake_outfile): + data = mock_interface._extract_data(fake_outfile, component="z") + np.testing.assert_array_almost_equal(data[0], [13.3, 14.4, 15.5]) + np.testing.assert_array_almost_equal(data[1], [16.6, 17.7, 18.8]) + + +def test_format_value_python_int(mock_interface): + assert mock_interface._format_value(3) == "3" + + +def test_format_value_numpy_int(mock_interface): + assert mock_interface._format_value(np.int64(7)) == "7" + + +def test_format_value_short_float(mock_interface): + assert mock_interface._format_value(1.5) == "1.5" + + +def test_format_value_long_float(mock_interface): + assert mock_interface._format_value(1.123456789) == "1.1235" + + +def test_format_value_numpy_float(mock_interface): + assert mock_interface._format_value(np.float64(1.123456789)) == "1.1235" + + +def test_format_value_non_numeric(mock_interface): + assert mock_interface._format_value("hello") == "hello" + + +def test_format_float_within_precision(mock_interface): + assert mock_interface._format_float(3.14) == "3.14" + + +def test_format_float_exceeds_precision(mock_interface): + assert mock_interface._format_float(3.141592653) == "3.1416" + + +def test_format_float_custom_precision(mock_interface): + mock_interface.float_precision = 2 + assert mock_interface._format_float(3.141592653) == "3.14" + + +def test_format_rows_single_row(mock_interface): + assert mock_interface._format_rows([[1, 2.5, 3]]) == "1 2.5 3" + + +def test_format_rows_multiple_rows(mock_interface): + assert mock_interface._format_rows([[1, 2], [3, 4]]) == "1 2\n3 4" diff --git a/tests/plate_simulation/runtest/leroi_test.py b/tests/plate_simulation/runtest/leroi_test.py index 3ebc0a4c..39ad853e 100644 --- a/tests/plate_simulation/runtest/leroi_test.py +++ b/tests/plate_simulation/runtest/leroi_test.py @@ -10,6 +10,33 @@ import subprocess +import numpy as np +from geoapps_utils.modelling.plates import PlateModel +from geoh5py.groups import SimPEGGroup + +from simpeg_drivers.electromagnetics.time_domain import ( + TDEMForwardOptions, +) +from simpeg_drivers.plate_simulation.driver import PlateSimulationDriver +from simpeg_drivers.plate_simulation.models.options import ( + ModelOptions, + OverburdenOptions, + PlateOptions, +) +from simpeg_drivers.plate_simulation.options import MeshOptions, PlateSimulationOptions +from simpeg_drivers.utils.synthetics.driver import SyntheticsComponents +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions as SyntheticsMeshOptions, +) +from simpeg_drivers.utils.synthetics.options import ( + ModelOptions as SyntheticsModelOptions, +) +from simpeg_drivers.utils.synthetics.options import ( + SurveyOptions, + SyntheticsComponentsOptions, +) +from tests.utils.targets import get_workspace + def test_leroi_executable(tmp_path): control_file = tmp_path / "test.cfl" @@ -18,4 +45,73 @@ def test_leroi_executable(tmp_path): def test_leroi_run(tmp_path): + opts = SyntheticsComponentsOptions( + method="airborne tdem", + survey=SurveyOptions(n_stations=8, n_lines=8, drape=40.0), + mesh=SyntheticsMeshOptions(), + model=SyntheticsModelOptions(background=1.0 / 2000.0), + ) + with get_workspace(tmp_path / "leroi_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + + mesh_params = MeshOptions( + u_cell_size=25.0, + v_cell_size=25.0, + w_cell_size=25.0, + padding_distance=1000.0, + depth_core=600.0, + max_distance=200.0, + survey_refinement=[2, 4], + topography_refinement=[0, 1], + plate_refinement=[2, 2], + ) + + overburden_params = OverburdenOptions( + thickness=50.0, + overburden_property=1500.0, # overburden resistivity (ohm-m) + ) + + plate_params = PlateOptions( + name="plate", + geometry=PlateModel( + easting=0.0, + northing=0.0, + elevation=-50.0, + width=10.0, + strike_length=200.0, + dip_length=100.0, + dip=80.0, + direction=90.0, + ), + plate_property=1.0, # plate resistivity (ohm-m) + ) + + model_params = ModelOptions( + background=2000.0, # background resistivity (ohm-m) + overburden_options=overburden_params, + plate_options=plate_params, + ) + + options = TDEMForwardOptions.build( + topography_object=components.topography, + data_object=components.survey, + geoh5=geoh5, + starting_model=1.0 / 2000.0, # background conductivity (S/m) + z_channel_bool=True, + ) + + tdem_group = SimPEGGroup.create(geoh5) + tdem_group.options = options.serialize() + + params = PlateSimulationOptions( + title="test_leroi", + run_command="run", + geoh5=geoh5, + mesh=mesh_params, + model=model_params, + simulation=tdem_group, + use_leroi=True, + ) + driver = PlateSimulationDriver(params) + driver.run() From 9f851e1ebb45817f21b358e54fa75a376f39b2eb Mon Sep 17 00:00:00 2001 From: benjamink Date: Fri, 17 Apr 2026 09:24:27 -0700 Subject: [PATCH 07/43] fix plate tests for new origin convention --- tests/plate_simulation/models/plates_test.py | 22 +++++++------------- 1 file changed, 8 insertions(+), 14 deletions(-) diff --git a/tests/plate_simulation/models/plates_test.py b/tests/plate_simulation/models/plates_test.py index aaecf84b..3cdac382 100644 --- a/tests/plate_simulation/models/plates_test.py +++ b/tests/plate_simulation/models/plates_test.py @@ -58,15 +58,9 @@ def test_vertical_east_striking_plate(tmp_path): vertical_east_striking.extent[1, 2] - vertical_east_striking.extent[0, 2], 500.0, ) - assert ( - vertical_east_striking.vertices[:, 0].mean() == 0.0 # pylint: disable=no-member - ) - assert ( - vertical_east_striking.vertices[:, 1].mean() == 0.0 # pylint: disable=no-member - ) - assert ( - vertical_east_striking.vertices[:, 2].mean() == 0.0 # pylint: disable=no-member - ) + assert np.isclose(vertical_east_striking.vertices[:, 0].mean(), 0.0) + assert np.isclose(vertical_east_striking.vertices[:, 1].mean(), 0.0) + assert np.isclose(vertical_east_striking.vertices[:, 2].mean(), -250.0) def test_dipping_plates_all_quadrants(tmp_path): @@ -112,11 +106,11 @@ def test_replicate_even(tmp_path): assert plates[1].surface(workspace).vertices is not None assert np.allclose( plates[0].surface(workspace).vertices.mean(axis=0), - np.array([-5.0, 0.0, 0.0]), + np.array([-5.0, 0.5, 0.0]), ) assert np.allclose( plates[1].surface(workspace).vertices.mean(axis=0), - np.array([5.0, 0.0, 0.0]), + np.array([5.0, 0.5, 0.0]), ) @@ -139,13 +133,13 @@ def test_replicate_odd(tmp_path): assert plates[2].surface(workspace).vertices is not None assert np.allclose( plates[0].surface(workspace).vertices.mean(axis=0), - np.array([0.0, -5.0, 0.0]), + np.array([0.0, -4.5, 0.0]), ) assert np.allclose( plates[1].surface(workspace).vertices.mean(axis=0), - np.array([0.0, 0.0, 0.0]), + np.array([0.0, 0.5, 0.0]), ) assert np.allclose( plates[2].surface(workspace).vertices.mean(axis=0), - np.array([0.0, 5.0, 0.0]), + np.array([0.0, 5.5, 0.0]), ) From 58c795688263a12150efea0906859ce845c581db Mon Sep 17 00:00:00 2001 From: benjamink Date: Fri, 17 Apr 2026 11:17:54 -0700 Subject: [PATCH 08/43] bulk up unit tests for options class --- .../plate_simulation/leroi_air/interface.py | 2 +- .../plate_simulation/leroi_air/options.py | 33 +++++++- tests/plate_simulation/leroi_air/__init__.py | 4 +- tests/plate_simulation/leroi_air/conftest.py | 8 ++ .../leroi_air/options_test.py | 78 +++++++++++++++++-- 5 files changed, 114 insertions(+), 11 deletions(-) diff --git a/simpeg_drivers/plate_simulation/leroi_air/interface.py b/simpeg_drivers/plate_simulation/leroi_air/interface.py index a6ff4803..c53c8e46 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/interface.py +++ b/simpeg_drivers/plate_simulation/leroi_air/interface.py @@ -65,7 +65,7 @@ def aliased_values(self) -> dict[str, Any]: "LINE_TAG": 0, "EAST": self.opts.locations[:, 0], "NORTH": self.opts.locations[:, 1], - "ALT": 13 * np.ones(self.opts.n_stations), # TODO - take from survey + "ALT": self.opts.drape_height, "NLAYER": self.opts.n_layers, "NPLATE": self.opts.n_plates, "NLITH": self.opts.n_layers + self.opts.n_plates, diff --git a/simpeg_drivers/plate_simulation/leroi_air/options.py b/simpeg_drivers/plate_simulation/leroi_air/options.py index 45727b42..9d47010e 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/options.py +++ b/simpeg_drivers/plate_simulation/leroi_air/options.py @@ -16,7 +16,9 @@ from geoapps_utils.modelling.plates import PlateModel from geoh5py.objects.surveys.electromagnetics.airborne_tem import AirborneTEMReceivers from pydantic import BaseModel, ConfigDict, Field, field_validator, model_validator +from scipy.interpolate import LinearNDInterpolator +from simpeg_drivers.components.topography import InversionTopography from simpeg_drivers.plate_simulation.options import PlateSimulationOptions @@ -26,6 +28,7 @@ class LeroiAirOptions(BaseModel): model_config = ConfigDict(arbitrary_types_allowed=True) survey: AirborneTEMReceivers + topo: np.ndarray layer_resistivities: list[float] layer_thicknesses: list[float] plate_resistivities: list[float] @@ -57,9 +60,22 @@ def validate_plate_lengths_match(self) -> Self: @classmethod def from_plate_simulation_options(cls, options: PlateSimulationOptions) -> Self: """Construct from a :class:`PlateSimulationOptions` instance.""" + simulation_options = options.simulation_parameters() + survey = simulation_options.data_object.copy() + + if simulation_options.active_cells.topography_object is None: + raise NotImplementedError( + "Passing active cells directly does not currently work for simulating " + "plates using the LeroiAir option. Simulation options must contain a " + "topography object." + ) + + topo_xyz = InversionTopography(survey.workspace, simulation_options).locations + return cls( survey=simulation_options.data_object, + topo=topo_xyz, layer_resistivities=[ options.model.overburden_options.overburden_property, options.model.background, @@ -78,7 +94,17 @@ def title(self) -> str: @property def locations(self) -> np.ndarray: """Survey receiver locations.""" - return self.survey.vertices + return self.survey.locations + + @property + def drape_height(self) -> np.ndarray: + """Survey height over topography.""" + + survey_locs = self.locations.copy() + topo_interp = LinearNDInterpolator(self.topo[:, :2], self.topo[:, 2]) + topo_at_survey_locations = topo_interp(survey_locs[:, 0], survey_locs[:, 1]) + + return survey_locs[:, 2] - topo_at_survey_locations @property def n_stations(self) -> int: @@ -155,7 +181,10 @@ def offtime(self) -> float: are not accounted for by the waveform. """ half_cycle = 1 / (2 * self.frequency) - ontime = self.waveform[np.where(self.waveform[:, 1] == 0)[0][1], 0] + zero_current_ind = np.where(self.waveform[:, 1] == 0)[0] + first_zero = 1 if self.waveform[0, 1] == 0.0 else 0 + ontime = self.waveform[zero_current_ind][first_zero, 1] + return half_cycle - ontime @property diff --git a/tests/plate_simulation/leroi_air/__init__.py b/tests/plate_simulation/leroi_air/__init__.py index a34ccbe9..bebd5700 100644 --- a/tests/plate_simulation/leroi_air/__init__.py +++ b/tests/plate_simulation/leroi_air/__init__.py @@ -22,7 +22,7 @@ def generate_plate_options(workspace): x = np.linspace(-1000, 1000, 81) y = np.linspace(-1000, 1000, 81) X, Y = np.meshgrid(x, y) - Z = np.zeros_like(X) + Z = np.full_like(X, 20.0) with fetch_active_workspace(workspace) as geoh5: survey = generate_airborne_tdem_survey(geoh5, X=X, Y=Y, Z=Z) @@ -38,8 +38,10 @@ def generate_plate_options(workspace): dip=45, ) ] + topo = np.column_stack([X.flatten(), Y.flatten(), np.zeros(X.size)]) opts = LeroiAirOptions( survey=survey, + topo=topo, layer_resistivities=layer_resistivities, layer_thicknesses=layer_thicknesses, plate_resistivities=plate_resistivities, diff --git a/tests/plate_simulation/leroi_air/conftest.py b/tests/plate_simulation/leroi_air/conftest.py index 41efd082..9101c456 100644 --- a/tests/plate_simulation/leroi_air/conftest.py +++ b/tests/plate_simulation/leroi_air/conftest.py @@ -68,6 +68,14 @@ def real_interface(tmp_path): return LeroiAirInterface(opts=opts) +@pytest.fixture(scope="module") +def plate_options(tmp_path_factory): + tmp_path = tmp_path_factory.mktemp("leroi_air_options") + with Workspace(tmp_path / "test.geoh5") as geoh5: + opts = generate_plate_options(workspace=geoh5) + return opts + + @pytest.fixture def fake_outfile(tmp_path): path = tmp_path / "leroi.out" diff --git a/tests/plate_simulation/leroi_air/options_test.py b/tests/plate_simulation/leroi_air/options_test.py index ee593b87..c2ee42cb 100644 --- a/tests/plate_simulation/leroi_air/options_test.py +++ b/tests/plate_simulation/leroi_air/options_test.py @@ -8,21 +8,85 @@ # ' # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' import numpy as np -from geoapps_utils.modelling.plates import PlateModel +import pytest from geoh5py import Workspace +from pydantic import ValidationError + +from simpeg_drivers.plate_simulation.leroi_air.options import LeroiAirOptions from . import generate_plate_options -def test_options(tmp_path): +def test_options_channel_widths(plate_options): + assert np.allclose(plate_options.channel_widths, [0.3, 0.3, 0.6]) + + +def test_options_channel_widths_from_metadata(tmp_path): + explicit_widths = [0.1, 0.2, 0.3] + with Workspace(tmp_path / "test.geoh5") as geoh5: + opts = generate_plate_options(workspace=geoh5) + meta = opts.survey.metadata + meta["channel_widths"] = explicit_widths + opts.survey.metadata = meta + + assert np.allclose(opts.channel_widths, explicit_widths) + + +def test_options_frequency(plate_options): + assert np.isclose(plate_options.frequency, 0.142857143) + +def test_options_frequency_from_metadata(tmp_path): with Workspace(tmp_path / "test.geoh5") as geoh5: opts = generate_plate_options(workspace=geoh5) + meta = opts.survey.metadata + meta["frequency"] = 0.15 + opts.survey.metadata = meta + + assert np.isclose(opts.frequency, 0.15) + + +def test_options_offtime(plate_options): + assert np.isclose(plate_options.offtime, 1.6) - assert np.allclose(opts.channel_widths, [0.3, 0.3, 0.6]) - assert np.isclose(opts.frequency, 0.142857143) - assert np.isclose(opts.offtime, 1.6) - assert np.allclose(opts.ontime_waveform, opts.waveform[0:8, :]) + +def test_options_ontime_waveform(plate_options): + assert np.allclose(plate_options.ontime_waveform, plate_options.waveform[0:8, :]) + + +def test_options_conductivity_thicknesses(plate_options): assert np.allclose( - opts.conductivity_thicknesses, np.array([0.0333333, 0.5, 0.4, 0.1]) + plate_options.conductivity_thicknesses, np.array([0.0333333, 0.5, 0.4, 0.1]) ) + + +def test_options_resistivities(plate_options): + assert np.allclose(plate_options.resistivities, [1500.0, 2000.0, 5000.0, 100.0]) + + +def test_options_drape_height(plate_options): + assert np.allclose(plate_options.drape_height, 20.0) + + +def test_options_validate_layer_lengths_mismatch(plate_options): + with pytest.raises(ValidationError, match="layer_resistivities"): + LeroiAirOptions( + survey=plate_options.survey, + topo=plate_options.topo, + layer_resistivities=[1500.0], + layer_thicknesses=[50.0, 1000.0], + plate_resistivities=plate_options.plate_resistivities, + plate_geometries=plate_options.plate_geometries, + ) + + +def test_options_validate_plate_lengths_mismatch(plate_options): + with pytest.raises(ValidationError, match="plate_resistivities"): + LeroiAirOptions( + survey=plate_options.survey, + topo=plate_options.topo, + layer_resistivities=plate_options.layer_resistivities, + layer_thicknesses=plate_options.layer_thicknesses, + plate_resistivities=[100.0, 200.0], + plate_geometries=plate_options.plate_geometries, + ) From 551e8af19d15b213269e59d9d181353a895b38bf Mon Sep 17 00:00:00 2001 From: benjamink Date: Fri, 17 Apr 2026 16:03:34 -0700 Subject: [PATCH 09/43] Update the runtests --- simpeg_drivers/plate_simulation/driver.py | 64 +++++--- .../plate_simulation/leroi_air/interface.py | 2 +- tests/plate_simulation/runtest/leroi_test.py | 137 ++++++++++-------- tests/utils/runtests.py | 24 +++ 4 files changed, 143 insertions(+), 84 deletions(-) create mode 100644 tests/utils/runtests.py diff --git a/simpeg_drivers/plate_simulation/driver.py b/simpeg_drivers/plate_simulation/driver.py index 1fb15868..31b2b99d 100644 --- a/simpeg_drivers/plate_simulation/driver.py +++ b/simpeg_drivers/plate_simulation/driver.py @@ -75,7 +75,11 @@ def run(self) -> InversionDriver: """Create octree mesh, fill model, and simulate.""" with fetch_active_workspace(self.params.geoh5, mode="r+"): + self._organize_out_group() + self.simulation_parameters.mesh = self.mesh + self.simulation_parameters.models.starting_model = self.model self.simulation_driver.run() + self._update_simulation_options() self.update_monitoring_directory(self._out_group) logger.info("done.") @@ -181,7 +185,7 @@ def make_mesh(self) -> Octree: name=self.params.mesh.name, ) - mesh.parent = self._out_group + mesh.parent = self.params.simulation return mesh @@ -298,30 +302,48 @@ def start_dask_run( """ start_dask_run(cls, json_path, n_workers=n_workers, n_threads=n_threads) - def _get_simpeg_driver(self): + def _organize_out_group(self) -> None: + """ + Place the simulation group inside out_group and copy topography there. - with fetch_active_workspace(self.params.geoh5, mode="r+"): - self.simulation_parameters.mesh = self.mesh - self.simulation_parameters.models.starting_model = self.model + Also updates the active-cells topography reference to the new copy so + that subsequent option serialization points to the object inside the + plate-simulation output group. + """ + self.params.simulation.parent = self._out_group + topo_copy = self.topography.copy(parent=self._out_group, copy_children=True) + self.simulation_parameters.active_cells.topography_object = topo_copy - if not isinstance( - self.simulation_parameters.active_cells.topography_object, - Surface | Points, - ): - raise ValueError( - "The topography object of the forward simulation must be a 'Surface'." - ) + def _update_simulation_options(self) -> None: + """ + Serialize current mesh, model, and topography into the simulation group. - self.simulation_parameters.out_group = None - driver_class = driver_class_from_name( - self.simulation_parameters.inversion_type, forward_only=True - ) - self._simulation_driver = driver_class( - self.simulation_parameters, - client=self._client, - workers=self._workers, + This keeps the TDEMForward group runnable as a standalone SimPEG + forward simulation for comparison with the LeroiAir result. + """ + self.simulation_parameters.out_group = self.params.simulation + self.simulation_parameters.update_out_group_options() + + def _get_simpeg_driver(self): + + if not isinstance( + self.simulation_parameters.active_cells.topography_object, + Surface | Points, + ): + raise ValueError( + "The topography object of the forward simulation must be a 'Surface'." ) - self._simulation_driver.out_group.parent = self._out_group + + self.simulation_parameters.out_group = None + driver_class = driver_class_from_name( + self.simulation_parameters.inversion_type, forward_only=True + ) + self._simulation_driver = driver_class( + self.simulation_parameters, + client=self._client, + workers=self._workers, + ) + self._simulation_driver.out_group.parent = self.params.simulation return self.simulation_driver diff --git a/simpeg_drivers/plate_simulation/leroi_air/interface.py b/simpeg_drivers/plate_simulation/leroi_air/interface.py index c53c8e46..3f4393ea 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/interface.py +++ b/simpeg_drivers/plate_simulation/leroi_air/interface.py @@ -307,4 +307,4 @@ def save_to_geoh5(self, outfile: str | Path, out_group: UIJsonGroup): dip_direction=plate.direction, dip=plate.dip, ) - MaxwellPlate.create(geoh5, geometry=geometry) + MaxwellPlate.create(geoh5, geometry=geometry, parent=out_group) diff --git a/tests/plate_simulation/runtest/leroi_test.py b/tests/plate_simulation/runtest/leroi_test.py index 39ad853e..e424b5c2 100644 --- a/tests/plate_simulation/runtest/leroi_test.py +++ b/tests/plate_simulation/runtest/leroi_test.py @@ -8,16 +8,15 @@ # ' # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -import subprocess - import numpy as np from geoapps_utils.modelling.plates import PlateModel +from geoh5py import Workspace from geoh5py.groups import SimPEGGroup +from geoh5py.objects import AirborneTEMReceivers, MaxwellPlate, Octree from simpeg_drivers.electromagnetics.time_domain import ( TDEMForwardOptions, ) -from simpeg_drivers.plate_simulation.driver import PlateSimulationDriver from simpeg_drivers.plate_simulation.models.options import ( ModelOptions, OverburdenOptions, @@ -35,83 +34,97 @@ SurveyOptions, SyntheticsComponentsOptions, ) +from tests.utils.runtests import run_driver_from_ui_json from tests.utils.targets import get_workspace -def test_leroi_executable(tmp_path): - control_file = tmp_path / "test.cfl" - control_file.touch() - subprocess.run("F2.bat test output", cwd=tmp_path, shell=True, check=False) - - def test_leroi_run(tmp_path): - opts = SyntheticsComponentsOptions( - method="airborne tdem", - survey=SurveyOptions(n_stations=8, n_lines=8, drape=40.0), - mesh=SyntheticsMeshOptions(), - model=SyntheticsModelOptions(background=1.0 / 2000.0), - ) - with get_workspace(tmp_path / "leroi_test.ui.geoh5") as geoh5: - components = SyntheticsComponents(geoh5, options=opts) - mesh_params = MeshOptions( - u_cell_size=25.0, - v_cell_size=25.0, - w_cell_size=25.0, - padding_distance=1000.0, - depth_core=600.0, - max_distance=200.0, - survey_refinement=[2, 4], - topography_refinement=[0, 1], - plate_refinement=[2, 2], + with Workspace(tmp_path / "leroi_test.geoh5") as geoh5: + geometry = PlateModel( + easting=0.0, + northing=0.0, + elevation=-50.0, + width=10.0, + strike_length=200.0, + dip_length=100.0, + dip=80.0, + direction=90.0, ) - overburden_params = OverburdenOptions( - thickness=50.0, - overburden_property=1500.0, # overburden resistivity (ohm-m) - ) - - plate_params = PlateOptions( - name="plate", - geometry=PlateModel( - easting=0.0, - northing=0.0, - elevation=-50.0, - width=10.0, - strike_length=200.0, - dip_length=100.0, - dip=80.0, - direction=90.0, + opts = SyntheticsComponentsOptions( + method="airborne tdem", + survey=SurveyOptions(width=800.0, n_stations=32, n_lines=4, drape=40.0), + mesh=SyntheticsMeshOptions(), + model=SyntheticsModelOptions( + anomaly=1 / 10, + background=1 / 5000, + plate=geometry, ), - plate_property=1.0, # plate resistivity (ohm-m) - ) - - model_params = ModelOptions( - background=2000.0, # background resistivity (ohm-m) - overburden_options=overburden_params, - plate_options=plate_params, ) - + components = SyntheticsComponents(geoh5, options=opts) + out_group = SimPEGGroup.create(geoh5, name="TEM forward") options = TDEMForwardOptions.build( topography_object=components.topography, data_object=components.survey, geoh5=geoh5, - starting_model=1.0 / 2000.0, # background conductivity (S/m) + starting_model=components.model, z_channel_bool=True, + out_group=out_group, ) - - tdem_group = SimPEGGroup.create(geoh5) - tdem_group.options = options.serialize() + options.update_out_group_options() params = PlateSimulationOptions( - title="test_leroi", - run_command="run", geoh5=geoh5, - mesh=mesh_params, - model=model_params, - simulation=tdem_group, + mesh=MeshOptions( + u_cell_size=25.0, + v_cell_size=25.0, + w_cell_size=25.0, + padding_distance=1000.0, + depth_core=600.0, + max_distance=200.0, + survey_refinement=[2, 4], + topography_refinement=[0, 1], + plate_refinement=[2, 2], + ), + model=ModelOptions( + background=2000.0, + overburden_options=OverburdenOptions( + thickness=50.0, + overburden_property=1500.0, + ), + plate_options=PlateOptions( + name="plate", + geometry=geometry, + plate_property=1.0, + ), + ), + simulation=options.out_group, use_leroi=True, ) - driver = PlateSimulationDriver(params) - driver.run() + run_driver_from_ui_json(params) + + with Workspace(tmp_path / "leroi_test.geoh5") as geoh5: + out_group = geoh5.get_entity("TEM forward")[0] + assert out_group is not None + + survey = out_group.get_entity("survey")[0] + assert isinstance(survey, AirborneTEMReceivers) + + maxwell_plate = out_group.get_entity("Maxwell Plate")[0] + assert isinstance(maxwell_plate, MaxwellPlate) + + expected_channels = [ + "fwd inline [0]", + "fwd inline [1]", + "fwd inline [2]", + "fwd crossline [0]", + "fwd crossline [1]", + "fwd crossline [2]", + "fwd vertical [0]", + "fwd vertical [1]", + "fwd vertical [2]", + ] + for channel in expected_channels: + assert survey.get_entity(channel)[0] is not None diff --git a/tests/utils/runtests.py b/tests/utils/runtests.py new file mode 100644 index 00000000..cab10791 --- /dev/null +++ b/tests/utils/runtests.py @@ -0,0 +1,24 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2026 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +from geoapps_utils.base import Driver, Options +from geoapps_utils.run import run_from_uijson + + +def run_driver_from_ui_json(params: Options, name: str = "runtest.ui.json") -> Driver: + """ + Serialize params to a ui.json file and run. + + :param params: Options instance to serialize. + :param name: Name of the ui.json file to write. + """ + path = params.geoh5.h5file.parent / name + params.write_ui_json(path) + return run_from_uijson(path) From 3a07436306247ad89ab07170b2ff4f6f5777188e Mon Sep 17 00:00:00 2001 From: benjamink Date: Tue, 21 Apr 2026 08:43:01 -0700 Subject: [PATCH 10/43] don't create mesh or model if use_leroi --- .../factories/directives_factory.py | 2 +- simpeg_drivers/plate_simulation/driver.py | 14 +++--- simpeg_drivers/plate_simulation/options.py | 46 +++++++++++++------ 3 files changed, 41 insertions(+), 21 deletions(-) diff --git a/simpeg_drivers/components/factories/directives_factory.py b/simpeg_drivers/components/factories/directives_factory.py index 94834769..42f73b2e 100644 --- a/simpeg_drivers/components/factories/directives_factory.py +++ b/simpeg_drivers/components/factories/directives_factory.py @@ -246,7 +246,7 @@ def save_iteration_log_files(self): """""" if self._save_iteration_log_files is None and self.driver.logger: self._save_iteration_log_files = directives.SaveLogFilesGeoH5( - self.driver.out_group, self.driver.logger.start_date_time + self.driver.out_group ) return self._save_iteration_log_files diff --git a/simpeg_drivers/plate_simulation/driver.py b/simpeg_drivers/plate_simulation/driver.py index 31b2b99d..b88a2dc3 100644 --- a/simpeg_drivers/plate_simulation/driver.py +++ b/simpeg_drivers/plate_simulation/driver.py @@ -76,8 +76,6 @@ def run(self) -> InversionDriver: with fetch_active_workspace(self.params.geoh5, mode="r+"): self._organize_out_group() - self.simulation_parameters.mesh = self.mesh - self.simulation_parameters.models.starting_model = self.model self.simulation_driver.run() self._update_simulation_options() self.update_monitoring_directory(self._out_group) @@ -90,9 +88,7 @@ def run(self) -> InversionDriver: @property def simulation_driver(self) -> InversionDriver: if self._simulation_driver is None: - run_command = self.simulation_parameters.run_command - is_tem = "time_domain.forward" in run_command - if self.params.use_leroi & is_tem: + if self.params.use_leroi: self._simulation_driver = self._get_leroi_driver() else: self._simulation_driver = self._get_simpeg_driver() @@ -102,11 +98,17 @@ def simulation_driver(self) -> InversionDriver: @property def simulation_parameters(self) -> BaseForwardOptions: if self._simulation_parameters is None: - self._simulation_parameters = self.params.simulation_parameters() + self._simulation_parameters = self.params.simulation_parameters.copy() + if self._simulation_parameters.physical_property == "conductivity": self._simulation_parameters.models.model_type = ( ModelTypeEnum.resistivity ) + + if not self.params.use_leroi: + self._simulation_parameters.mesh = self.mesh + self._simulation_parameters.models.starting_model = self.model + return self._simulation_parameters @property diff --git a/simpeg_drivers/plate_simulation/options.py b/simpeg_drivers/plate_simulation/options.py index 8564468a..1d731014 100644 --- a/simpeg_drivers/plate_simulation/options.py +++ b/simpeg_drivers/plate_simulation/options.py @@ -15,6 +15,7 @@ from geoapps_utils.base import Options from geoh5py.groups import SimPEGGroup, UIJsonGroup from geoh5py.ui_json import InputFile +from pydantic import model_validator from simpeg_drivers import assets_path from simpeg_drivers.electricals.direct_current.three_dimensions.options import ( @@ -79,7 +80,18 @@ class PlateSimulationOptions(Options): model: ModelOptions simulation: SimPEGGroup | UIJsonGroup use_leroi: bool = False - + _simulation_parameters: BaseForwardOptions | None = None + + @model_validator(mode="before") + @classmethod + def use_leroi_em_only(cls, data) -> bool: + run_command = data["simulation"].options["run_command"] + is_tem = "time_domain.forward" in run_command + use_leroi = data.get("use_leroi", False) & is_tem + data["use_leroi"] = use_leroi + return data + + @property def simulation_parameters(self) -> BaseForwardOptions: """ Create SimPEG parameters from the simulation options. @@ -87,21 +99,27 @@ def simulation_parameters(self) -> BaseForwardOptions: A new SimPEGGroup is created inside the out_group to store the result of the forward simulation. """ - simulation_options = deepcopy(self.simulation.options) - simulation_options["geoh5"] = self.geoh5 + if self._simulation_parameters is None: + simulation_options = deepcopy(self.simulation.options) + simulation_options["geoh5"] = self.geoh5 + + # TODO replace InputFile.data with UIJson.to_params + input_file = InputFile(ui_json=simulation_options, validate=False) + if input_file.ui_json is None: + raise ValueError("Input file must have ui_json set.") - input_file = InputFile(ui_json=simulation_options, validate=False) - if input_file.ui_json is None: - raise ValueError("Input file must have ui_json set.") + input_file.ui_json["mesh"]["value"] = None - input_file.ui_json["mesh"]["value"] = None + if input_file.data is None: + raise ValueError("Input file data must be set.") - if input_file.data is None: - raise ValueError("Input file data must be set.") + inversion_type = input_file.data["inversion_type"] - if input_file.data["inversion_type"] in PARAM_MAP: - return PARAM_MAP[input_file.data["inversion_type"]].build(input_file.data) + if inversion_type in PARAM_MAP: + self._simulation_parameters = PARAM_MAP[inversion_type].build( + input_file.data + ) + else: + raise NotImplementedError(f"Unknown inversion type: {inversion_type}") - raise NotImplementedError( - f"Unknown inversion type: {input_file.data['inversion_type']}" - ) + return self._simulation_parameters From 999e98569f2e002042249897262c40467c079856 Mon Sep 17 00:00:00 2001 From: benjamink Date: Tue, 21 Apr 2026 09:12:38 -0700 Subject: [PATCH 11/43] simplify and remove _organize_out_group --- simpeg_drivers/plate_simulation/driver.py | 14 +------------- 1 file changed, 1 insertion(+), 13 deletions(-) diff --git a/simpeg_drivers/plate_simulation/driver.py b/simpeg_drivers/plate_simulation/driver.py index f59a7e48..c9c901b4 100644 --- a/simpeg_drivers/plate_simulation/driver.py +++ b/simpeg_drivers/plate_simulation/driver.py @@ -75,7 +75,7 @@ def run(self) -> InversionDriver: """Create octree mesh, fill model, and simulate.""" with fetch_active_workspace(self.params.geoh5, mode="r+"): - self._organize_out_group() + self.params.simulation.parent = self._out_group self.simulation_driver.run() self._update_simulation_options() self.update_monitoring_directory(self._out_group) @@ -298,18 +298,6 @@ def start_dask_run( """ start_dask_run(cls, json_path, n_workers=n_workers, n_threads=n_threads) - def _organize_out_group(self) -> None: - """ - Place the simulation group inside out_group and copy topography there. - - Also updates the active-cells topography reference to the new copy so - that subsequent option serialization points to the object inside the - plate-simulation output group. - """ - self.params.simulation.parent = self._out_group - topo_copy = self.topography.copy(parent=self._out_group, copy_children=True) - self.simulation_parameters.active_cells.topography_object = topo_copy - def _update_simulation_options(self) -> None: """ Serialize current mesh, model, and topography into the simulation group. From 45e420203991cc3883dfdb7aad9530260f5df56b Mon Sep 17 00:00:00 2001 From: benjamink Date: Tue, 21 Apr 2026 09:55:21 -0700 Subject: [PATCH 12/43] Add Maxwell Plate creation for non-leroi runs --- simpeg_drivers/plate_simulation/driver.py | 4 ++++ .../plate_simulation/leroi_air/interface.py | 18 ++---------------- .../plate_simulation/leroi_air/options.py | 2 +- tests/plate_simulation/runtest/leroi_test.py | 2 +- 4 files changed, 8 insertions(+), 18 deletions(-) diff --git a/simpeg_drivers/plate_simulation/driver.py b/simpeg_drivers/plate_simulation/driver.py index c9c901b4..8748d6a3 100644 --- a/simpeg_drivers/plate_simulation/driver.py +++ b/simpeg_drivers/plate_simulation/driver.py @@ -145,6 +145,10 @@ def plates(self) -> list[Plate]: self.params.model.plate_options.spacing, self.params.model.plate_options.geometry.direction, ) + with fetch_active_workspace(self.params.geoh5, mode="r+") as geoh5: + for plate in self._plates: + plate.to_maxwell_plate(geoh5, parent=self._out_group) + return self._plates @property diff --git a/simpeg_drivers/plate_simulation/leroi_air/interface.py b/simpeg_drivers/plate_simulation/leroi_air/interface.py index 3f4393ea..41056e8f 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/interface.py +++ b/simpeg_drivers/plate_simulation/leroi_air/interface.py @@ -84,7 +84,7 @@ def aliased_values(self) -> dict[str, Any]: "IPLATE": 1, "CNTR_East": [g.easting for g in self.opts.plate_geometries], "CNTR_North": [g.northing for g in self.opts.plate_geometries], - "PLTOP": [g.elevation for g in self.opts.plate_geometries], + "PLTOP": [-1 * g.elevation for g in self.opts.plate_geometries], "PLNGTH": [g.strike_length for g in self.opts.plate_geometries], "DPWDTH": [g.dip_length for g in self.opts.plate_geometries], "DZM": [g.direction for g in self.opts.plate_geometries], @@ -276,7 +276,7 @@ def save_to_geoh5(self, outfile: str | Path, out_group: UIJsonGroup): inline_data = self._extract_data(outfile=outfile, component="y") vertical_data = self._extract_data(outfile=outfile, component="z") - with fetch_active_workspace(self.opts.survey.workspace, mode="r+") as geoh5: + with fetch_active_workspace(self.opts.survey.workspace, mode="r+"): survey = self.opts.survey.copy(parent=out_group, copy_children=False) data = survey.add_data( { @@ -294,17 +294,3 @@ def save_to_geoh5(self, outfile: str | Path, out_group: UIJsonGroup): survey.create_property_group(name="inline", properties=data[:3]) survey.create_property_group(name="crossline", properties=data[3:6]) survey.create_property_group(name="vertical", properties=data[6:]) - - for plate in self.opts.plate_geometries: - position = PlatePosition( - x=plate.easting, y=plate.northing, z=plate.elevation - ) - geometry = PlateGeometry( - position=position, - length=plate.strike_length, - width=plate.dip_length, - thickness=plate.width, - dip_direction=plate.direction, - dip=plate.dip, - ) - MaxwellPlate.create(geoh5, geometry=geometry, parent=out_group) diff --git a/simpeg_drivers/plate_simulation/leroi_air/options.py b/simpeg_drivers/plate_simulation/leroi_air/options.py index 9d47010e..13a415fe 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/options.py +++ b/simpeg_drivers/plate_simulation/leroi_air/options.py @@ -61,7 +61,7 @@ def validate_plate_lengths_match(self) -> Self: def from_plate_simulation_options(cls, options: PlateSimulationOptions) -> Self: """Construct from a :class:`PlateSimulationOptions` instance.""" - simulation_options = options.simulation_parameters() + simulation_options = options.simulation_parameters survey = simulation_options.data_object.copy() if simulation_options.active_cells.topography_object is None: diff --git a/tests/plate_simulation/runtest/leroi_test.py b/tests/plate_simulation/runtest/leroi_test.py index e424b5c2..850d35b4 100644 --- a/tests/plate_simulation/runtest/leroi_test.py +++ b/tests/plate_simulation/runtest/leroi_test.py @@ -44,7 +44,7 @@ def test_leroi_run(tmp_path): geometry = PlateModel( easting=0.0, northing=0.0, - elevation=-50.0, + elevation=100.0, width=10.0, strike_length=200.0, dip_length=100.0, From 16e9eeccbf6e77274b79a65d0a9b1925471a0fb8 Mon Sep 17 00:00:00 2001 From: benjamink Date: Tue, 21 Apr 2026 12:38:56 -0700 Subject: [PATCH 13/43] Results in the right place: foward group under plate_simulation group --- simpeg_drivers/components/meshes.py | 2 ++ simpeg_drivers/plate_simulation/driver.py | 20 +++++-------------- simpeg_drivers/plate_simulation/options.py | 1 + .../plate_simulation/runtest/gravity_test.py | 8 +++----- 4 files changed, 11 insertions(+), 20 deletions(-) diff --git a/simpeg_drivers/components/meshes.py b/simpeg_drivers/components/meshes.py index aa0a4436..b9dff812 100644 --- a/simpeg_drivers/components/meshes.py +++ b/simpeg_drivers/components/meshes.py @@ -103,6 +103,8 @@ def get_entity(self) -> Octree | DrapeModel: "No mesh provided. Creating optimized mesh from data and topography." ) mesh_entity = self._auto_mesh() + elif self.params.mesh.parent == self.params.out_group: # plate-simulation + mesh_entity = self.params.mesh else: mesh_entity = self.params.mesh.copy( parent=self.params.out_group, copy_children=False diff --git a/simpeg_drivers/plate_simulation/driver.py b/simpeg_drivers/plate_simulation/driver.py index 8748d6a3..7dfa1d62 100644 --- a/simpeg_drivers/plate_simulation/driver.py +++ b/simpeg_drivers/plate_simulation/driver.py @@ -75,9 +75,9 @@ def run(self) -> InversionDriver: """Create octree mesh, fill model, and simulate.""" with fetch_active_workspace(self.params.geoh5, mode="r+"): - self.params.simulation.parent = self._out_group self.simulation_driver.run() - self._update_simulation_options() + self.simulation_parameters.update_out_group_options() + self.params.simulation.parent = self._out_group self.update_monitoring_directory(self._out_group) logger.info("done.") @@ -89,6 +89,7 @@ def run(self) -> InversionDriver: def simulation_driver(self) -> InversionDriver: if self._simulation_driver is None: if self.params.use_leroi: + _ = self.plates # Saves MaxwellPlate(s) when no octree/model self._simulation_driver = self._get_leroi_driver() else: self._simulation_driver = self._get_simpeg_driver() @@ -302,16 +303,6 @@ def start_dask_run( """ start_dask_run(cls, json_path, n_workers=n_workers, n_threads=n_threads) - def _update_simulation_options(self) -> None: - """ - Serialize current mesh, model, and topography into the simulation group. - - This keeps the TDEMForward group runnable as a standalone SimPEG - forward simulation for comparison with the LeroiAir result. - """ - self.simulation_parameters.out_group = self.params.simulation - self.simulation_parameters.update_out_group_options() - def _get_simpeg_driver(self): if not isinstance( @@ -322,7 +313,6 @@ def _get_simpeg_driver(self): "The topography object of the forward simulation must be a 'Surface'." ) - self.simulation_parameters.out_group = None driver_class = driver_class_from_name( self.simulation_parameters.inversion_type, forward_only=True ) @@ -331,14 +321,14 @@ def _get_simpeg_driver(self): client=self._client, workers=self._workers, ) - self._simulation_driver.out_group.parent = self.params.simulation return self.simulation_driver def _get_leroi_driver(self): + # assert self.params.geoh5.get_entity("TEM forward")[0].parent.name == "Plate Simulation" leroi_opts = LeroiAirOptions.from_plate_simulation_options(self.params) driver = LeroiAirDriver(leroi_opts) - driver.out_group = self._out_group + driver.out_group = self.params.simulation return driver diff --git a/simpeg_drivers/plate_simulation/options.py b/simpeg_drivers/plate_simulation/options.py index 1d731014..2f0afcdc 100644 --- a/simpeg_drivers/plate_simulation/options.py +++ b/simpeg_drivers/plate_simulation/options.py @@ -102,6 +102,7 @@ def simulation_parameters(self) -> BaseForwardOptions: if self._simulation_parameters is None: simulation_options = deepcopy(self.simulation.options) simulation_options["geoh5"] = self.geoh5 + simulation_options["out_group"] = self.simulation # TODO replace InputFile.data with UIJson.to_params input_file = InputFile(ui_json=simulation_options, validate=False) diff --git a/tests/plate_simulation/runtest/gravity_test.py b/tests/plate_simulation/runtest/gravity_test.py index 4f7d618d..a098e940 100644 --- a/tests/plate_simulation/runtest/gravity_test.py +++ b/tests/plate_simulation/runtest/gravity_test.py @@ -84,16 +84,14 @@ def test_gravity_plate_simulation(tmp_path): starting_model=0.1, ) - gravity_inversion = SimPEGGroup.create(geoh5) - gravity_inversion.options = options.serialize() + gravity_forward = SimPEGGroup.create(geoh5, name="gravity forward") + gravity_forward.options = options.serialize() params = PlateSimulationOptions( - title="test", - run_command="run", geoh5=geoh5, mesh=mesh_params, model=model_params, - simulation=gravity_inversion, + simulation=gravity_forward, ) driver = PlateSimulationDriver(params) driver.run() From 8c792fca09d7978f5996414dc373111c3bdb3ac0 Mon Sep 17 00:00:00 2001 From: benk-mira <81254271+benk-mira@users.noreply.github.com> Date: Tue, 21 Apr 2026 13:02:32 -0700 Subject: [PATCH 14/43] Update tests/plate_simulation/runtest/leroi_test.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- tests/plate_simulation/runtest/leroi_test.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/plate_simulation/runtest/leroi_test.py b/tests/plate_simulation/runtest/leroi_test.py index 850d35b4..0bd4a39c 100644 --- a/tests/plate_simulation/runtest/leroi_test.py +++ b/tests/plate_simulation/runtest/leroi_test.py @@ -35,7 +35,6 @@ SyntheticsComponentsOptions, ) from tests.utils.runtests import run_driver_from_ui_json -from tests.utils.targets import get_workspace def test_leroi_run(tmp_path): From 49ad96781752fdfe533e8d096ee6b322225fba9f Mon Sep 17 00:00:00 2001 From: benjamink Date: Tue, 21 Apr 2026 13:32:02 -0700 Subject: [PATCH 15/43] fix save_to_geoh5 to handle n channels --- simpeg_drivers/plate_simulation/driver.py | 3 +- .../plate_simulation/leroi_air/interface.py | 45 ++++++++----------- .../plate_simulation/leroi_air/options.py | 2 +- .../leroi_air/interface_test.py | 2 +- tests/plate_simulation/runtest/driver_test.py | 2 +- tests/plate_simulation/runtest/leroi_test.py | 9 ++-- 6 files changed, 27 insertions(+), 36 deletions(-) diff --git a/simpeg_drivers/plate_simulation/driver.py b/simpeg_drivers/plate_simulation/driver.py index 7dfa1d62..2aa4ff30 100644 --- a/simpeg_drivers/plate_simulation/driver.py +++ b/simpeg_drivers/plate_simulation/driver.py @@ -322,10 +322,9 @@ def _get_simpeg_driver(self): workers=self._workers, ) - return self.simulation_driver + return self._simulation_driver def _get_leroi_driver(self): - # assert self.params.geoh5.get_entity("TEM forward")[0].parent.name == "Plate Simulation" leroi_opts = LeroiAirOptions.from_plate_simulation_options(self.params) driver = LeroiAirDriver(leroi_opts) driver.out_group = self.params.simulation diff --git a/simpeg_drivers/plate_simulation/leroi_air/interface.py b/simpeg_drivers/plate_simulation/leroi_air/interface.py index 41056e8f..4d9ecf91 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/interface.py +++ b/simpeg_drivers/plate_simulation/leroi_air/interface.py @@ -240,9 +240,9 @@ def write_cfl_file(self, filepath: Path) -> None: f.write(self.format_cfl_file()) _COMPONENT_ANCHORS: dict[str, str] = { - "x": "TRANSVERSE COMPONENT", - "y": "IN-LINE COMPONENT", - "z": "VERTICAL COMPONENT", + "crossline": "TRANSVERSE COMPONENT", + "inline": "IN-LINE COMPONENT", + "vertical": "VERTICAL COMPONENT", } def _find_data_start(self, chunk: list[str]) -> int: @@ -262,35 +262,26 @@ def _slice_data_lines(self, lines: list[str], anchor: str) -> list[str]: return chunk[data_start : data_start + self.opts.n_stations] def _extract_data( - self, outfile: str | Path, component: Literal["x", "y", "z"] + self, outfile: str | Path, component: Literal["inline", "crossline", "vertical"] ) -> np.ndarray: """Extract channel data for a single component from a LeroiAir .out file.""" lines = Path(outfile).read_text(encoding="utf-8", errors="replace").splitlines() data_lines = self._slice_data_lines(lines, self._COMPONENT_ANCHORS[component]) return np.array([line.split() for line in data_lines], dtype=float)[:, 4:] - def save_to_geoh5(self, outfile: str | Path, out_group: UIJsonGroup): + def save_to_geoh5(self, outfile: str | Path, out_group): """Save LeroiAir simulated data on the provided survey to geoh5.""" - crossline_data = self._extract_data(outfile=outfile, component="x") - inline_data = self._extract_data(outfile=outfile, component="y") - vertical_data = self._extract_data(outfile=outfile, component="z") - - with fetch_active_workspace(self.opts.survey.workspace, mode="r+"): - survey = self.opts.survey.copy(parent=out_group, copy_children=False) - data = survey.add_data( - { - "fwd inline [0]": {"values": inline_data[:, 0]}, - "fwd inline [1]": {"values": inline_data[:, 1]}, - "fwd inline [2]": {"values": inline_data[:, 2]}, - "fwd crossline [0]": {"values": crossline_data[:, 0]}, - "fwd crossline [1]": {"values": crossline_data[:, 1]}, - "fwd crossline [2]": {"values": crossline_data[:, 2]}, - "fwd vertical [0]": {"values": vertical_data[:, 0]}, - "fwd vertical [1]": {"values": vertical_data[:, 1]}, - "fwd vertical [2]": {"values": vertical_data[:, 2]}, - } - ) - survey.create_property_group(name="inline", properties=data[:3]) - survey.create_property_group(name="crossline", properties=data[3:6]) - survey.create_property_group(name="vertical", properties=data[6:]) + survey = self.opts.survey.copy(parent=out_group, copy_children=False) + for component in "inline", "crossline", "vertical": + data = self._extract_data(outfile=outfile, component=component) + + with fetch_active_workspace(self.opts.survey.workspace, mode="r+"): + entities = survey.add_data( + { + f"fwd {component} [{i}]": {"values": data[:, i]} + for i in range(len(self.opts.channels)) + } + ) + + survey.create_property_group(name=component, properties=entities) diff --git a/simpeg_drivers/plate_simulation/leroi_air/options.py b/simpeg_drivers/plate_simulation/leroi_air/options.py index 13a415fe..c5792f66 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/options.py +++ b/simpeg_drivers/plate_simulation/leroi_air/options.py @@ -183,7 +183,7 @@ def offtime(self) -> float: half_cycle = 1 / (2 * self.frequency) zero_current_ind = np.where(self.waveform[:, 1] == 0)[0] first_zero = 1 if self.waveform[0, 1] == 0.0 else 0 - ontime = self.waveform[zero_current_ind][first_zero, 1] + ontime = self.waveform[zero_current_ind][first_zero, 0] return half_cycle - ontime diff --git a/tests/plate_simulation/leroi_air/interface_test.py b/tests/plate_simulation/leroi_air/interface_test.py index c718d0ce..e876f425 100644 --- a/tests/plate_simulation/leroi_air/interface_test.py +++ b/tests/plate_simulation/leroi_air/interface_test.py @@ -103,7 +103,7 @@ def test_format_cfl_file(tmp_path): assert "10.0 \t ! CELLW" in cfl_str # record_14: LITHP, CNTR_East, CNTR_North, PLTOP - assert "4 0.0 0.0 0.0\t ! LITHP, CNTR_East, CNTR_North, PLTOP" in cfl_str + assert "4 0.0 0.0 -0.0\t ! LITHP, CNTR_East, CNTR_North, PLTOP" in cfl_str # record_15: PLNGTH, DPWDTH, DZM, DIP assert "200.0 100.0 90.0 45.0\t ! PLNGTH, DPWDTH, DZM, DIP" in cfl_str diff --git a/tests/plate_simulation/runtest/driver_test.py b/tests/plate_simulation/runtest/driver_test.py index 297237c4..5bb7575d 100644 --- a/tests/plate_simulation/runtest/driver_test.py +++ b/tests/plate_simulation/runtest/driver_test.py @@ -94,7 +94,7 @@ def test_plate_simulation_params_from_input_file(tmp_path, caplog): assert "Overburden thickness exceeds the plate depth" in caplog.text assert isinstance(params.simulation, SimPEGGroup) - simulation_parameters = params.simulation_parameters() + simulation_parameters = params.simulation_parameters assert simulation_parameters.inversion_type == "gravity" assert simulation_parameters.forward_only diff --git a/tests/plate_simulation/runtest/leroi_test.py b/tests/plate_simulation/runtest/leroi_test.py index 850d35b4..5fc8d54b 100644 --- a/tests/plate_simulation/runtest/leroi_test.py +++ b/tests/plate_simulation/runtest/leroi_test.py @@ -106,13 +106,14 @@ def test_leroi_run(tmp_path): run_driver_from_ui_json(params) with Workspace(tmp_path / "leroi_test.geoh5") as geoh5: - out_group = geoh5.get_entity("TEM forward")[0] - assert out_group is not None + plate_simulation_group = geoh5.get_entity("Plate Simulation")[0] + forward_group = geoh5.get_entity("TEM forward")[0] + assert forward_group.parent == plate_simulation_group - survey = out_group.get_entity("survey")[0] + survey = forward_group.get_entity("survey")[0] assert isinstance(survey, AirborneTEMReceivers) - maxwell_plate = out_group.get_entity("Maxwell Plate")[0] + maxwell_plate = plate_simulation_group.get_entity("Maxwell Plate")[0] assert isinstance(maxwell_plate, MaxwellPlate) expected_channels = [ From f7739bc7c2a11a807999d083857fb09af547a9f2 Mon Sep 17 00:00:00 2001 From: benk-mira <81254271+benk-mira@users.noreply.github.com> Date: Tue, 21 Apr 2026 13:34:08 -0700 Subject: [PATCH 16/43] Update simpeg_drivers/plate_simulation/options.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- simpeg_drivers/plate_simulation/options.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/simpeg_drivers/plate_simulation/options.py b/simpeg_drivers/plate_simulation/options.py index 2f0afcdc..0e9cbefe 100644 --- a/simpeg_drivers/plate_simulation/options.py +++ b/simpeg_drivers/plate_simulation/options.py @@ -84,10 +84,10 @@ class PlateSimulationOptions(Options): @model_validator(mode="before") @classmethod - def use_leroi_em_only(cls, data) -> bool: + def use_leroi_em_only(cls, data) -> dict: run_command = data["simulation"].options["run_command"] is_tem = "time_domain.forward" in run_command - use_leroi = data.get("use_leroi", False) & is_tem + use_leroi = data.get("use_leroi", False) and is_tem data["use_leroi"] = use_leroi return data From c69414c7b61533cea418f73bb482b8e496b38488 Mon Sep 17 00:00:00 2001 From: benk-mira <81254271+benk-mira@users.noreply.github.com> Date: Tue, 21 Apr 2026 13:34:47 -0700 Subject: [PATCH 17/43] Update simpeg_drivers/plate_simulation/leroi_air/interface.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- simpeg_drivers/plate_simulation/leroi_air/interface.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/simpeg_drivers/plate_simulation/leroi_air/interface.py b/simpeg_drivers/plate_simulation/leroi_air/interface.py index 4d9ecf91..ffa56122 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/interface.py +++ b/simpeg_drivers/plate_simulation/leroi_air/interface.py @@ -14,8 +14,6 @@ import numpy as np from geoh5py.groups import UIJsonGroup -from geoh5py.objects import MaxwellPlate -from geoh5py.objects.maxwell_plate import PlateGeometry, PlatePosition from geoh5py.shared.utils import fetch_active_workspace from .options import LeroiAirOptions From 602eca4b9e6dc7e09d2a371f12c71b0376889161 Mon Sep 17 00:00:00 2001 From: benjamink Date: Wed, 22 Apr 2026 13:43:41 -0700 Subject: [PATCH 18/43] initialize simulation_parameters on construction with lazy load mesh/model as required by simpeg runs --- simpeg_drivers/plate_simulation/driver.py | 109 ++++++++++-------- .../plate_simulation/leroi_air/driver.py | 10 +- .../plate_simulation/leroi_air/interface.py | 27 ++--- .../plate_simulation/leroi_air/options.py | 32 ++--- simpeg_drivers/plate_simulation/options.py | 62 +++++----- simpeg_drivers/utils/utils.py | 58 ++++++++-- .../leroi_air/interface_test.py | 10 +- tests/plate_simulation/runtest/driver_test.py | 11 -- .../plate_simulation/runtest/gravity_test.py | 94 +++++++-------- tests/plate_simulation/runtest/leroi_test.py | 7 +- 10 files changed, 229 insertions(+), 191 deletions(-) diff --git a/simpeg_drivers/plate_simulation/driver.py b/simpeg_drivers/plate_simulation/driver.py index 2aa4ff30..d3c32c76 100644 --- a/simpeg_drivers/plate_simulation/driver.py +++ b/simpeg_drivers/plate_simulation/driver.py @@ -11,6 +11,7 @@ from __future__ import annotations import sys +from copy import deepcopy from pathlib import Path import numpy as np @@ -21,6 +22,7 @@ from geoh5py.data import FloatData, ReferencedData from geoh5py.objects import Octree, Points, Surface from geoh5py.shared.utils import fetch_active_workspace +from geoh5py.ui_json.input_file import InputFile from simpeg_drivers.driver import ( InversionDriver, @@ -35,7 +37,11 @@ from simpeg_drivers.plate_simulation.models.series import DikeSwarm, Geology from simpeg_drivers.plate_simulation.options import PlateSimulationOptions from simpeg_drivers.utils.synthetics.meshes import get_octree_mesh -from simpeg_drivers.utils.utils import start_dask_run, validate_out_group +from simpeg_drivers.utils.utils import ( + driver_class_from_dict, + start_dask_run, + validate_out_group, +) logger = get_logger(__name__, propagate=False) @@ -66,19 +72,19 @@ def __init__( self._survey: Points | None = None self._mesh: Octree | None = None self._model: FloatData | None = None - self._simulation_parameters: BaseForwardOptions | None = None - self._simulation_driver: InversionDriver | None = None self._client: Client | bool = validate_client(client) self._workers: list[tuple[str]] = validate_workers(self._client, workers) + self._simulation_driver: InversionDriver | None = None + self.simulation_parameters: BaseForwardOptions = self._initialize_forward_opts() def run(self) -> InversionDriver: """Create octree mesh, fill model, and simulate.""" - with fetch_active_workspace(self.params.geoh5, mode="r+"): + with fetch_active_workspace(self.simulation_parameters.geoh5, mode="r+"): self.simulation_driver.run() self.simulation_parameters.update_out_group_options() - self.params.simulation.parent = self._out_group - self.update_monitoring_directory(self._out_group) + # self.params.simulation.parent = self._out_group + # self.update_monitoring_directory(self._out_group) logger.info("done.") logger.handlers.clear() @@ -96,28 +102,9 @@ def simulation_driver(self) -> InversionDriver: return self._simulation_driver - @property - def simulation_parameters(self) -> BaseForwardOptions: - if self._simulation_parameters is None: - self._simulation_parameters = self.params.simulation_parameters.copy() - - if self._simulation_parameters.physical_property == "conductivity": - self._simulation_parameters.models.model_type = ( - ModelTypeEnum.resistivity - ) - - if not self.params.use_leroi: - self._simulation_parameters.mesh = self.mesh - self._simulation_parameters.models.starting_model = self.model - - return self._simulation_parameters - @property def survey(self): - if self._survey is None: - self._survey = self.simulation_parameters.data_object - - return self._survey + return self.simulation_parameters.data_object @property def topography(self) -> Surface | Points: @@ -152,22 +139,6 @@ def plates(self) -> list[Plate]: return self._plates - @property - def mesh(self) -> Octree: - """Returns an octree mesh built from mesh parameters.""" - if self._mesh is None: - self._mesh = self.make_mesh() - - return self._mesh - - @property - def model(self) -> FloatData: - """Returns the model built from model parameters.""" - if self._model is None: - self._model = self.make_model() - - return self._model - def make_mesh(self) -> Octree: """ Build specialized mesh for plate simulation from parameters. @@ -186,7 +157,7 @@ def make_mesh(self) -> Octree: name=self.params.mesh.name, ) - mesh.parent = self.params.simulation + mesh.parent = self.simulation_parameters.out_group return mesh @@ -215,7 +186,7 @@ def make_model(self) -> FloatData: scenario = Geology( workspace=self.params.geoh5, - mesh=self.mesh, + mesh=self.simulation_parameters.mesh, background=self.params.model.background, history=[dikes, overburden, erosion], ) @@ -228,7 +199,7 @@ def make_model(self) -> FloatData: if physical_property == "conductivity": physical_property = "resistivity" - model = self.mesh.add_data( + model = self.simulation_parameters.mesh.add_data( { "geology": { "type": "referenced", @@ -244,7 +215,7 @@ def make_model(self) -> FloatData: for k, v in physical_property_map.items(): starting_model_values[geology == k] = v - starting_model = self.mesh.add_data( + starting_model = self.simulation_parameters.mesh.add_data( {"starting_model": {"values": starting_model_values}} ) @@ -316,6 +287,8 @@ def _get_simpeg_driver(self): driver_class = driver_class_from_name( self.simulation_parameters.inversion_type, forward_only=True ) + self.simulation_parameters.mesh = self.make_mesh() + self.simulation_parameters.models.starting_model = self.make_model() self._simulation_driver = driver_class( self.simulation_parameters, client=self._client, @@ -325,11 +298,51 @@ def _get_simpeg_driver(self): return self._simulation_driver def _get_leroi_driver(self): - leroi_opts = LeroiAirOptions.from_plate_simulation_options(self.params) + leroi_opts = LeroiAirOptions.from_plate_simulation_options( + self.params.model, self.simulation_parameters + ) driver = LeroiAirDriver(leroi_opts) - driver.out_group = self.params.simulation return driver + def _collect_simulation_opts(self) -> BaseForwardOptions: + """Collect template simulation options.""" + + simulation_options = deepcopy(self.params.simulation.options) + simulation_options["geoh5"] = self.params.geoh5 + + # TODO replace InputFile.data with UIJson.to_params + input_file = InputFile(ui_json=simulation_options, validate=False) + driver = driver_class_from_dict(input_file.data) + + return driver._params_class.build(input_file.data) # pylint: disable=protected-access + + def _initialize_forward_opts(self) -> BaseForwardOptions: + """Initialize the forward simulation options with mesh and model.""" + + opts = self._collect_simulation_opts() + + update = {} + models_update = {} + if opts.physical_property == "conductivity": + # TODO: validate this logic + models_update["model_type"] = ModelTypeEnum.resistivity + if not self.params.use_leroi: + update["mesh"] = None + models_update["starting_model"] = None + update["models"] = opts.models.model_copy(update=models_update) + + with fetch_active_workspace(self.params.geoh5, mode="r+"): + out_group = opts.out_group.copy( + parent=self.out_group, + copy_children=False, + copy_relatives=False, + ) + update["out_group"] = out_group + forward_opts = opts.model_copy(update=update) + forward_opts.update_out_group_options() + + return forward_opts + if __name__ == "__main__": file = Path(sys.argv[1]) diff --git a/simpeg_drivers/plate_simulation/leroi_air/driver.py b/simpeg_drivers/plate_simulation/leroi_air/driver.py index f9d89a65..4e2d4838 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/driver.py +++ b/simpeg_drivers/plate_simulation/leroi_air/driver.py @@ -13,6 +13,8 @@ from geoh5py.groups import UIJsonGroup +from simpeg_drivers.utils.utils import validate_out_group + from .interface import LeroiAirInterface from .options import LeroiAirOptions @@ -20,11 +22,13 @@ class LeroiAirDriver: """Orchestrates a LeroiAir forward simulation from input preparation to geoh5 output.""" - def __init__(self, options: LeroiAirOptions): + def __init__( + self, + options: LeroiAirOptions, + ) -> None: """Initialize with simulation options.""" self.options = options self._interface: LeroiAirInterface | None = None - self.out_group: UIJsonGroup | None = None @property def interface(self) -> LeroiAirInterface: @@ -61,5 +65,5 @@ def run(self) -> None: self.run_leroi() self.interface.save_to_geoh5( outfile=self.project_path / "LeroiAir.out", - out_group=self.out_group, + out_group=self.options.out_group, ) diff --git a/simpeg_drivers/plate_simulation/leroi_air/interface.py b/simpeg_drivers/plate_simulation/leroi_air/interface.py index 4d9ecf91..755336f0 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/interface.py +++ b/simpeg_drivers/plate_simulation/leroi_air/interface.py @@ -272,16 +272,17 @@ def _extract_data( def save_to_geoh5(self, outfile: str | Path, out_group): """Save LeroiAir simulated data on the provided survey to geoh5.""" - survey = self.opts.survey.copy(parent=out_group, copy_children=False) - for component in "inline", "crossline", "vertical": - data = self._extract_data(outfile=outfile, component=component) - - with fetch_active_workspace(self.opts.survey.workspace, mode="r+"): - entities = survey.add_data( - { - f"fwd {component} [{i}]": {"values": data[:, i]} - for i in range(len(self.opts.channels)) - } - ) - - survey.create_property_group(name=component, properties=entities) + with out_group.workspace.open(mode="r+"): + survey = self.opts.survey.copy(parent=out_group, copy_children=False) + for component in "inline", "crossline", "vertical": + data = self._extract_data(outfile=outfile, component=component) + + with fetch_active_workspace(self.opts.survey.workspace, mode="r+"): + entities = survey.add_data( + { + f"fwd {component} [{i}]": {"values": data[:, i]} + for i in range(len(self.opts.channels)) + } + ) + + survey.create_property_group(name=component, properties=entities) diff --git a/simpeg_drivers/plate_simulation/leroi_air/options.py b/simpeg_drivers/plate_simulation/leroi_air/options.py index c5792f66..678d840f 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/options.py +++ b/simpeg_drivers/plate_simulation/leroi_air/options.py @@ -14,12 +14,13 @@ import numpy as np from geoapps_utils.modelling.plates import PlateModel +from geoh5py.groups import SimPEGGroup from geoh5py.objects.surveys.electromagnetics.airborne_tem import AirborneTEMReceivers -from pydantic import BaseModel, ConfigDict, Field, field_validator, model_validator +from pydantic import BaseModel, ConfigDict, model_validator from scipy.interpolate import LinearNDInterpolator from simpeg_drivers.components.topography import InversionTopography -from simpeg_drivers.plate_simulation.options import PlateSimulationOptions +from simpeg_drivers.plate_simulation.options import ModelOptions class LeroiAirOptions(BaseModel): @@ -38,6 +39,7 @@ class LeroiAirOptions(BaseModel): domain: Literal["time", "frequency"] = "time" layered_earth_only: bool = False float_precision: int = 4 + out_group: SimPEGGroup | None = None @model_validator(mode="after") def validate_layer_lengths_match(self) -> Self: @@ -58,32 +60,34 @@ def validate_plate_lengths_match(self) -> Self: return self @classmethod - def from_plate_simulation_options(cls, options: PlateSimulationOptions) -> Self: + def from_plate_simulation_options( + cls, model: ModelOptions, simulation: SimPEGGroup + ) -> Self: """Construct from a :class:`PlateSimulationOptions` instance.""" - simulation_options = options.simulation_parameters - survey = simulation_options.data_object.copy() + survey = simulation.data_object - if simulation_options.active_cells.topography_object is None: + if simulation.active_cells.topography_object is None: raise NotImplementedError( "Passing active cells directly does not currently work for simulating " "plates using the LeroiAir option. Simulation options must contain a " "topography object." ) - topo_xyz = InversionTopography(survey.workspace, simulation_options).locations + topo_xyz = InversionTopography(survey.workspace, simulation).locations return cls( - survey=simulation_options.data_object, + survey=survey, topo=topo_xyz, layer_resistivities=[ - options.model.overburden_options.overburden_property, - options.model.background, + model.overburden_options.overburden_property, + model.background, ], - layer_thicknesses=[options.model.overburden_options.thickness, 9999], - plate_resistivities=[options.model.plate_options.plate_property], - plate_geometries=[options.model.plate_options.geometry], - magnetic_field="dBdt" if "dBdt" in simulation_options.data_units else "B", + layer_thicknesses=[model.overburden_options.thickness, 9999], + plate_resistivities=[model.plate_options.plate_property], + plate_geometries=[model.plate_options.geometry], + magnetic_field="dBdt" if "dBdt" in simulation.data_units else "B", + out_group=simulation.out_group, ) @property diff --git a/simpeg_drivers/plate_simulation/options.py b/simpeg_drivers/plate_simulation/options.py index 2f0afcdc..60d75d03 100644 --- a/simpeg_drivers/plate_simulation/options.py +++ b/simpeg_drivers/plate_simulation/options.py @@ -8,7 +8,6 @@ # ' # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -from copy import deepcopy from pathlib import Path from typing import ClassVar @@ -34,7 +33,6 @@ MTForwardOptions, ) from simpeg_drivers.natural_sources.tipper.options import TipperForwardOptions -from simpeg_drivers.options import BaseForwardOptions from simpeg_drivers.potential_fields.gravity.options import GravityForwardOptions from simpeg_drivers.potential_fields.magnetic_vector import ( MagneticVectorForwardOptions, @@ -80,7 +78,7 @@ class PlateSimulationOptions(Options): model: ModelOptions simulation: SimPEGGroup | UIJsonGroup use_leroi: bool = False - _simulation_parameters: BaseForwardOptions | None = None + # _simulation_parameters: BaseForwardOptions | None = None @model_validator(mode="before") @classmethod @@ -91,36 +89,28 @@ def use_leroi_em_only(cls, data) -> bool: data["use_leroi"] = use_leroi return data - @property - def simulation_parameters(self) -> BaseForwardOptions: - """ - Create SimPEG parameters from the simulation options. - - A new SimPEGGroup is created inside the out_group to store the - result of the forward simulation. - """ - if self._simulation_parameters is None: - simulation_options = deepcopy(self.simulation.options) - simulation_options["geoh5"] = self.geoh5 - simulation_options["out_group"] = self.simulation - - # TODO replace InputFile.data with UIJson.to_params - input_file = InputFile(ui_json=simulation_options, validate=False) - if input_file.ui_json is None: - raise ValueError("Input file must have ui_json set.") - - input_file.ui_json["mesh"]["value"] = None - - if input_file.data is None: - raise ValueError("Input file data must be set.") - - inversion_type = input_file.data["inversion_type"] - - if inversion_type in PARAM_MAP: - self._simulation_parameters = PARAM_MAP[inversion_type].build( - input_file.data - ) - else: - raise NotImplementedError(f"Unknown inversion type: {inversion_type}") - - return self._simulation_parameters + # @property + # def simulation_parameters(self) -> BaseForwardOptions: + # """ + # Create SimPEG parameters from the simulation options. + # + # A new SimPEGGroup is created inside the out_group to store the + # result of the forward simulation. + # """ + # if self.simulation_parameters is None: + # simulation_options = deepcopy(self.simulation.options) + # simulation_options["geoh5"] = self.geoh5 + # simulation_options["out_group"] = self.simulation + # + # # TODO replace InputFile.data with UIJson.to_params + # input_file = InputFile(ui_json=simulation_options, validate=False) + # if input_file.ui_json is None: + # raise ValueError("Input file must have ui_json set.") + # + # input_file.ui_json["mesh"]["value"] = None + # + # if input_file.data is None: + # raise ValueError("Input file data must be set.") + # + # driver = driver_class_from_dict(input_file.data) + # return driver._params_class.build(input_file.data) # pylint: disable=protected-access diff --git a/simpeg_drivers/utils/utils.py b/simpeg_drivers/utils/utils.py index fc4c3016..0cf07f08 100644 --- a/simpeg_drivers/utils/utils.py +++ b/simpeg_drivers/utils/utils.py @@ -24,8 +24,9 @@ from dask.distributed import Client, LocalCluster, performance_report from discretize import TensorMesh, TreeMesh from discretize.utils import mesh_utils -from geoapps_utils.base import Options -from geoapps_utils.run import load_ui_json_as_dict +from geoapps_utils import GeoAppsError +from geoapps_utils.base import Driver, Options +from geoapps_utils.run import fetch_driver_class_from_string, load_ui_json_as_dict from geoapps_utils.utils.locations import mask_under_horizon from geoapps_utils.utils.numerical import running_mean, traveling_salesman from geoh5py import Workspace @@ -613,14 +614,7 @@ def simpeg_group_to_driver(group: SimPEGGroup, workspace: Workspace) -> Inversio ui_json["geoh5"] = workspace ifile = InputFile(ui_json=ui_json) - forward_only = ui_json.get("forward_only", False) - mod_name, classes = DRIVER_MAP.get(ui_json["inversion_type"]) - if forward_only: - class_name = classes.get("forward", classes["inversion"]) - else: - class_name = classes.get("inversion") - module = __import__(mod_name, fromlist=[class_name]) - inversion_driver = getattr(module, class_name) + inversion_driver = driver_class_from_dict(ifile.ui_json) ifile.set_data_value("out_group", group) params = inversion_driver._params_class.build(ifile) # pylint: disable=protected-access @@ -758,3 +752,47 @@ def start_dask_run( ps = pstats.Stats(profiler, stream=s) ps.sort_stats("cumulative") ps.print_stats() + + +def driver_class_from_name( + name: str, forward_only: bool = False +) -> type[InversionDriver]: + """ + Get the driver class from the inversion type name. + + TODO: Only for backward compatibility. To be deprecated in future versions. + + :param name: The inversion type name. + :param forward_only: Whether to forward only the forward operators. + + :return: InversionDriver or ForwardDriver class. + """ + if name not in DRIVER_MAP: + msg = f"Inversion type '{name}' is not supported." + msg += f" Valid inversions are: {(*list(DRIVER_MAP),)}." + raise NotImplementedError(msg) + + mod_name, classes = DRIVER_MAP.get(name) + class_name = classes.get("inversion") + if forward_only: + class_name = classes.get("forward", class_name) + + module = __import__(mod_name, fromlist=[class_name]) + return getattr(module, class_name) + + +def driver_class_from_dict(data: dict) -> type[InversionDriver | Driver]: + + inversion_type = data.get("inversion_type", None) + + if inversion_type: + forward_only = data.get("forward_only", False) + return driver_class_from_name(inversion_type, forward_only=forward_only) + + if data.get("run_command", None): + return fetch_driver_class_from_string(data["run_command"]) + + raise GeoAppsError( + "Could not find a driver for the given 'inversion_type' or 'run_command' parameter. " + "Please review the 'run_command' in the UI JSON." + ) diff --git a/tests/plate_simulation/leroi_air/interface_test.py b/tests/plate_simulation/leroi_air/interface_test.py index e876f425..02a9de09 100644 --- a/tests/plate_simulation/leroi_air/interface_test.py +++ b/tests/plate_simulation/leroi_air/interface_test.py @@ -111,7 +111,7 @@ def test_format_cfl_file(tmp_path): def test_slice_data_lines_returns_correct_station_rows(mock_interface): lines = FAKE_OUT.splitlines() - anchor = LeroiAirInterface._COMPONENT_ANCHORS["x"] + anchor = LeroiAirInterface._COMPONENT_ANCHORS["crossline"] rows = mock_interface._slice_data_lines(lines, anchor) assert len(rows) == N_STATIONS assert rows[0].split()[4] == "1.1" @@ -119,24 +119,24 @@ def test_slice_data_lines_returns_correct_station_rows(mock_interface): def test_extract_data_returns_channel_columns_only(mock_interface, fake_outfile): - data = mock_interface._extract_data(fake_outfile, component="x") + data = mock_interface._extract_data(fake_outfile, component="crossline") assert data.shape == (N_STATIONS, N_CHANNELS) def test_extract_data_parses_transverse_component_values(mock_interface, fake_outfile): - data = mock_interface._extract_data(fake_outfile, component="x") + data = mock_interface._extract_data(fake_outfile, component="crossline") np.testing.assert_array_almost_equal(data[0], [1.1, 2.2, 3.3]) np.testing.assert_array_almost_equal(data[1], [4.4, 5.5, 6.6]) def test_extract_data_parses_inline_component_values(mock_interface, fake_outfile): - data = mock_interface._extract_data(fake_outfile, component="y") + data = mock_interface._extract_data(fake_outfile, component="inline") np.testing.assert_array_almost_equal(data[0], [7.7, 8.8, 9.9]) np.testing.assert_array_almost_equal(data[1], [10.0, 11.1, 12.2]) def test_extract_data_parses_vertical_component_values(mock_interface, fake_outfile): - data = mock_interface._extract_data(fake_outfile, component="z") + data = mock_interface._extract_data(fake_outfile, component="vertical") np.testing.assert_array_almost_equal(data[0], [13.3, 14.4, 15.5]) np.testing.assert_array_almost_equal(data[1], [16.6, 17.7, 18.8]) diff --git a/tests/plate_simulation/runtest/driver_test.py b/tests/plate_simulation/runtest/driver_test.py index 5bb7575d..026c1746 100644 --- a/tests/plate_simulation/runtest/driver_test.py +++ b/tests/plate_simulation/runtest/driver_test.py @@ -94,17 +94,6 @@ def test_plate_simulation_params_from_input_file(tmp_path, caplog): assert "Overburden thickness exceeds the plate depth" in caplog.text assert isinstance(params.simulation, SimPEGGroup) - simulation_parameters = params.simulation_parameters - - assert simulation_parameters.inversion_type == "gravity" - assert simulation_parameters.forward_only - assert simulation_parameters.geoh5.h5file == geoh5.h5file - assert ( - simulation_parameters.active_cells.topography_object.uid - == components.topography.uid - ) - assert simulation_parameters.data_object.uid == components.survey.uid - assert isinstance(params.mesh, MeshOptions) assert params.mesh.u_cell_size == 10.0 assert params.mesh.v_cell_size == 10.0 diff --git a/tests/plate_simulation/runtest/gravity_test.py b/tests/plate_simulation/runtest/gravity_test.py index a098e940..bf0b8f0a 100644 --- a/tests/plate_simulation/runtest/gravity_test.py +++ b/tests/plate_simulation/runtest/gravity_test.py @@ -12,7 +12,6 @@ from geoapps_utils.modelling.plates import PlateModel from geoh5py.groups import SimPEGGroup -from simpeg_drivers.plate_simulation.driver import PlateSimulationDriver from simpeg_drivers.plate_simulation.models.options import ( ModelOptions, OverburdenOptions, @@ -31,69 +30,70 @@ SurveyOptions, SyntheticsComponentsOptions, ) +from tests.utils.runtests import run_driver_from_ui_json from tests.utils.targets import get_workspace def test_gravity_plate_simulation(tmp_path): - opts = SyntheticsComponentsOptions( - method="gravity", - survey=SurveyOptions(n_stations=8, n_lines=8, drape=5.0), - mesh=SyntheticsMeshOptions(), - model=SyntheticsModelOptions(anomaly=0.0), - ) - with get_workspace(tmp_path / "inversion_test.ui.geoh5") as geoh5: - components = SyntheticsComponents(geoh5, options=opts) - mesh_params = MeshOptions( - u_cell_size=10.0, - v_cell_size=10.0, - w_cell_size=10.0, - padding_distance=1500.0, - depth_core=600.0, - max_distance=200.0, - survey_refinement=[4, 6], - topography_refinement=[0, 1], - plate_refinement=[4, 2], - ) - overburden_params = OverburdenOptions(thickness=50.0, overburden_property=0.2) - - plate_params = PlateOptions( - name="plate", - geometry=PlateModel( - elevation=100.0, - width=100.0, - strike_length=100.0, - dip_length=100.0, - dip=0.0, - direction=0.0, - ), - plate_property=0.5, + with get_workspace(tmp_path / "inversion_test.ui.geoh5") as geoh5: + geometry = PlateModel( + elevation=100.0, + width=100.0, + strike_length=100.0, + dip_length=100.0, + dip=0.0, + direction=0.0, ) - model_params = ModelOptions( - name="density", - background=0.0, - overburden_options=overburden_params, - plate_options=plate_params, + opts = SyntheticsComponentsOptions( + method="gravity", + survey=SurveyOptions(n_stations=8, n_lines=8, drape=5.0), + mesh=SyntheticsMeshOptions(), + model=SyntheticsModelOptions(anomaly=0.0), ) + components = SyntheticsComponents(geoh5, options=opts) + out_group = SimPEGGroup.create(geoh5, name="Gravity forward") options = GravityForwardOptions.build( topography_object=components.topography, data_object=components.survey, geoh5=geoh5, starting_model=0.1, + out_group=out_group, ) - - gravity_forward = SimPEGGroup.create(geoh5, name="gravity forward") - gravity_forward.options = options.serialize() + options.update_out_group_options() params = PlateSimulationOptions( geoh5=geoh5, - mesh=mesh_params, - model=model_params, - simulation=gravity_forward, + mesh=MeshOptions( + u_cell_size=10.0, + v_cell_size=10.0, + w_cell_size=10.0, + padding_distance=1500.0, + depth_core=600.0, + max_distance=200.0, + survey_refinement=[4, 6], + topography_refinement=[0, 1], + plate_refinement=[4, 2], + ), + model=ModelOptions( + name="density", + background=0.0, + overburden_options=OverburdenOptions( + thickness=50.0, overburden_property=0.2 + ), + plate_options=PlateOptions( + name="plate", + geometry=geometry, + plate_property=0.5, + ), + ), + simulation=options.out_group, ) - driver = PlateSimulationDriver(params) - driver.run() - assert np.nanmax(driver.model.values) == 0.5 + driver = run_driver_from_ui_json(params) + + assert ( + np.nanmax(driver.simulation_parameters.models.starting_model.values) == 0.5 + ) diff --git a/tests/plate_simulation/runtest/leroi_test.py b/tests/plate_simulation/runtest/leroi_test.py index 1be9e8a3..c65374fe 100644 --- a/tests/plate_simulation/runtest/leroi_test.py +++ b/tests/plate_simulation/runtest/leroi_test.py @@ -8,7 +8,6 @@ # ' # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -import numpy as np from geoapps_utils.modelling.plates import PlateModel from geoh5py import Workspace from geoh5py.groups import SimPEGGroup @@ -35,11 +34,12 @@ SyntheticsComponentsOptions, ) from tests.utils.runtests import run_driver_from_ui_json +from tests.utils.targets import get_workspace def test_leroi_run(tmp_path): - with Workspace(tmp_path / "leroi_test.geoh5") as geoh5: + with get_workspace(tmp_path / "leroi_test.geoh5") as geoh5: geometry = PlateModel( easting=0.0, northing=0.0, @@ -106,8 +106,7 @@ def test_leroi_run(tmp_path): with Workspace(tmp_path / "leroi_test.geoh5") as geoh5: plate_simulation_group = geoh5.get_entity("Plate Simulation")[0] - forward_group = geoh5.get_entity("TEM forward")[0] - assert forward_group.parent == plate_simulation_group + forward_group = plate_simulation_group.get_entity("TEM forward")[0] survey = forward_group.get_entity("survey")[0] assert isinstance(survey, AirborneTEMReceivers) From cb275a4c208ee2c2df479d87392930ca741ce31c Mon Sep 17 00:00:00 2001 From: benjamink Date: Thu, 23 Apr 2026 10:43:18 -0700 Subject: [PATCH 19/43] fix sweep test --- simpeg_drivers/plate_simulation/driver.py | 3 ++- simpeg_drivers/plate_simulation/options.py | 27 ---------------------- 2 files changed, 2 insertions(+), 28 deletions(-) diff --git a/simpeg_drivers/plate_simulation/driver.py b/simpeg_drivers/plate_simulation/driver.py index d3c32c76..df607d25 100644 --- a/simpeg_drivers/plate_simulation/driver.py +++ b/simpeg_drivers/plate_simulation/driver.py @@ -332,7 +332,8 @@ def _initialize_forward_opts(self) -> BaseForwardOptions: update["models"] = opts.models.model_copy(update=models_update) with fetch_active_workspace(self.params.geoh5, mode="r+"): - out_group = opts.out_group.copy( + out_group = validate_out_group(opts) + out_group = out_group.copy( parent=self.out_group, copy_children=False, copy_relatives=False, diff --git a/simpeg_drivers/plate_simulation/options.py b/simpeg_drivers/plate_simulation/options.py index 6d3c846b..5ea80a6f 100644 --- a/simpeg_drivers/plate_simulation/options.py +++ b/simpeg_drivers/plate_simulation/options.py @@ -78,7 +78,6 @@ class PlateSimulationOptions(Options): model: ModelOptions simulation: SimPEGGroup | UIJsonGroup use_leroi: bool = False - # _simulation_parameters: BaseForwardOptions | None = None @model_validator(mode="before") @classmethod @@ -88,29 +87,3 @@ def use_leroi_em_only(cls, data) -> dict: use_leroi = data.get("use_leroi", False) and is_tem data["use_leroi"] = use_leroi return data - - # @property - # def simulation_parameters(self) -> BaseForwardOptions: - # """ - # Create SimPEG parameters from the simulation options. - # - # A new SimPEGGroup is created inside the out_group to store the - # result of the forward simulation. - # """ - # if self.simulation_parameters is None: - # simulation_options = deepcopy(self.simulation.options) - # simulation_options["geoh5"] = self.geoh5 - # simulation_options["out_group"] = self.simulation - # - # # TODO replace InputFile.data with UIJson.to_params - # input_file = InputFile(ui_json=simulation_options, validate=False) - # if input_file.ui_json is None: - # raise ValueError("Input file must have ui_json set.") - # - # input_file.ui_json["mesh"]["value"] = None - # - # if input_file.data is None: - # raise ValueError("Input file data must be set.") - # - # driver = driver_class_from_dict(input_file.data) - # return driver._params_class.build(input_file.data) # pylint: disable=protected-access From d12133a6b5d751626747fd454ce6dba717e236d8 Mon Sep 17 00:00:00 2001 From: benjamink Date: Thu, 23 Apr 2026 10:52:40 -0700 Subject: [PATCH 20/43] cleanup --- simpeg_drivers/plate_simulation/driver.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/simpeg_drivers/plate_simulation/driver.py b/simpeg_drivers/plate_simulation/driver.py index df607d25..9bc16e19 100644 --- a/simpeg_drivers/plate_simulation/driver.py +++ b/simpeg_drivers/plate_simulation/driver.py @@ -83,8 +83,7 @@ def run(self) -> InversionDriver: with fetch_active_workspace(self.simulation_parameters.geoh5, mode="r+"): self.simulation_driver.run() self.simulation_parameters.update_out_group_options() - # self.params.simulation.parent = self._out_group - # self.update_monitoring_directory(self._out_group) + self.update_monitoring_directory(self._out_group) logger.info("done.") logger.handlers.clear() From bba8f2246bee7a54cc6a49e391165543ba6db5ca Mon Sep 17 00:00:00 2001 From: benjamink Date: Thu, 23 Apr 2026 11:06:47 -0700 Subject: [PATCH 21/43] update monitoring directory without crash on read-only workspace --- simpeg_drivers/plate_simulation/driver.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/simpeg_drivers/plate_simulation/driver.py b/simpeg_drivers/plate_simulation/driver.py index 9bc16e19..960909ec 100644 --- a/simpeg_drivers/plate_simulation/driver.py +++ b/simpeg_drivers/plate_simulation/driver.py @@ -83,7 +83,8 @@ def run(self) -> InversionDriver: with fetch_active_workspace(self.simulation_parameters.geoh5, mode="r+"): self.simulation_driver.run() self.simulation_parameters.update_out_group_options() - self.update_monitoring_directory(self._out_group) + with fetch_active_workspace(self.params.geoh5, mode="r+"): + self.update_monitoring_directory(self._out_group) logger.info("done.") logger.handlers.clear() From 84e3e34b7a30491d57799ecb9a61b0b538db3ad8 Mon Sep 17 00:00:00 2001 From: benk-mira <81254271+benk-mira@users.noreply.github.com> Date: Thu, 23 Apr 2026 13:29:55 -0700 Subject: [PATCH 22/43] Update simpeg_drivers/plate_simulation/leroi_air/options.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- simpeg_drivers/plate_simulation/leroi_air/options.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/simpeg_drivers/plate_simulation/leroi_air/options.py b/simpeg_drivers/plate_simulation/leroi_air/options.py index 678d840f..b9b55d1e 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/options.py +++ b/simpeg_drivers/plate_simulation/leroi_air/options.py @@ -158,7 +158,7 @@ def channel_widths(self) -> np.ndarray: """Time channel widths.""" channel_widths = self.survey.metadata.get("channel_widths", None) if channel_widths is None: - channel_widths = np.diff([0] + self.channels) + channel_widths = np.diff(np.r_[0, self.channels]) return channel_widths @property From b9c8ac62aed392103bcf1e6b03fe405cb436cdcb Mon Sep 17 00:00:00 2001 From: benk-mira <81254271+benk-mira@users.noreply.github.com> Date: Fri, 24 Apr 2026 08:13:06 -0700 Subject: [PATCH 23/43] Update simpeg_drivers/plate_simulation/leroi_air/__init__.py Co-authored-by: domfournier --- simpeg_drivers/plate_simulation/leroi_air/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/simpeg_drivers/plate_simulation/leroi_air/__init__.py b/simpeg_drivers/plate_simulation/leroi_air/__init__.py index bd9d6d2b..1b788624 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/__init__.py +++ b/simpeg_drivers/plate_simulation/leroi_air/__init__.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2026 Mira Geoscience Ltd. ' +# Copyright (c) 2023-2026 Mira Geoscience Ltd. ' # ' # This file is part of simpeg-drivers package. ' # ' From 2c1f4e1b960cceeac680cdb080cd4aa8d1ce6856 Mon Sep 17 00:00:00 2001 From: benjamink Date: Fri, 24 Apr 2026 09:49:40 -0700 Subject: [PATCH 24/43] forward and plate-simulation each get a copy of octree with geology saved on plate-simulation mesh and starting_model saved on forward mesh. Fix mesh naming --- simpeg_drivers/plate_simulation/driver.py | 11 +++++------ tests/plate_simulation/runtest/gravity_test.py | 1 + 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/simpeg_drivers/plate_simulation/driver.py b/simpeg_drivers/plate_simulation/driver.py index 68de9411..db9c909f 100644 --- a/simpeg_drivers/plate_simulation/driver.py +++ b/simpeg_drivers/plate_simulation/driver.py @@ -148,17 +148,16 @@ def make_mesh(self) -> Octree: logger.info("making the mesh...") with fetch_active_workspace(self.params.geoh5, mode="r+") as geoh5: surfaces = [p.surface(geoh5) for p in self.plates] - mesh = get_octree_mesh( + self._mesh = get_octree_mesh( opts=self.params.mesh, survey=self.survey, topography=self.simulation_parameters.active_cells.topography_object, plates=surfaces, - name=self.params.mesh.name, + name="Octree", ) + self._mesh.parent = self._out_group - mesh.parent = self.simulation_parameters.out_group - - return mesh + return self._mesh.copy(parent=self.simulation_parameters.out_group) def make_model(self) -> FloatData: """Create background + plate and overburden model from parameters.""" @@ -198,7 +197,7 @@ def make_model(self) -> FloatData: if physical_property == "conductivity": physical_property = "resistivity" - model = self.simulation_parameters.mesh.add_data( + model = self._mesh.add_data( { "geology": { "type": "referenced", diff --git a/tests/plate_simulation/runtest/gravity_test.py b/tests/plate_simulation/runtest/gravity_test.py index bf0b8f0a..25401f6c 100644 --- a/tests/plate_simulation/runtest/gravity_test.py +++ b/tests/plate_simulation/runtest/gravity_test.py @@ -76,6 +76,7 @@ def test_gravity_plate_simulation(tmp_path): survey_refinement=[4, 6], topography_refinement=[0, 1], plate_refinement=[4, 2], + name="mesh", ), model=ModelOptions( name="density", From 00153a3161234ab2194339a291be6f517f8ade8c Mon Sep 17 00:00:00 2001 From: benjamink Date: Fri, 24 Apr 2026 10:18:45 -0700 Subject: [PATCH 25/43] even better organization for simpeg runs --- simpeg_drivers/components/meshes.py | 4 ++-- simpeg_drivers/plate_simulation/driver.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/simpeg_drivers/components/meshes.py b/simpeg_drivers/components/meshes.py index b9dff812..8e2d63f2 100644 --- a/simpeg_drivers/components/meshes.py +++ b/simpeg_drivers/components/meshes.py @@ -103,8 +103,8 @@ def get_entity(self) -> Octree | DrapeModel: "No mesh provided. Creating optimized mesh from data and topography." ) mesh_entity = self._auto_mesh() - elif self.params.mesh.parent == self.params.out_group: # plate-simulation - mesh_entity = self.params.mesh + # elif self.params.mesh.parent == self.params.out_group: # plate-simulation + # mesh_entity = self.params.mesh else: mesh_entity = self.params.mesh.copy( parent=self.params.out_group, copy_children=False diff --git a/simpeg_drivers/plate_simulation/driver.py b/simpeg_drivers/plate_simulation/driver.py index db9c909f..6d2f3cad 100644 --- a/simpeg_drivers/plate_simulation/driver.py +++ b/simpeg_drivers/plate_simulation/driver.py @@ -157,7 +157,7 @@ def make_mesh(self) -> Octree: ) self._mesh.parent = self._out_group - return self._mesh.copy(parent=self.simulation_parameters.out_group) + return self._mesh def make_model(self) -> FloatData: """Create background + plate and overburden model from parameters.""" From 517ce5f258213dbcdecc951b3c2522a78a7ce6cc Mon Sep 17 00:00:00 2001 From: benjamink Date: Fri, 24 Apr 2026 10:21:12 -0700 Subject: [PATCH 26/43] cleanup --- simpeg_drivers/components/meshes.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/simpeg_drivers/components/meshes.py b/simpeg_drivers/components/meshes.py index 8e2d63f2..aa0a4436 100644 --- a/simpeg_drivers/components/meshes.py +++ b/simpeg_drivers/components/meshes.py @@ -103,8 +103,6 @@ def get_entity(self) -> Octree | DrapeModel: "No mesh provided. Creating optimized mesh from data and topography." ) mesh_entity = self._auto_mesh() - # elif self.params.mesh.parent == self.params.out_group: # plate-simulation - # mesh_entity = self.params.mesh else: mesh_entity = self.params.mesh.copy( parent=self.params.out_group, copy_children=False From 207a532cb8ad101741b1230312df146374d54b70 Mon Sep 17 00:00:00 2001 From: benjamink Date: Fri, 24 Apr 2026 12:19:59 -0700 Subject: [PATCH 27/43] split the interface class into input and output classes --- .../plate_simulation/leroi_air/driver.py | 20 ++----- .../plate_simulation/leroi_air/interface.py | 12 +++- tests/plate_simulation/leroi_air/conftest.py | 19 ++++--- .../leroi_air/interface_test.py | 56 ++++++++++--------- 4 files changed, 54 insertions(+), 53 deletions(-) diff --git a/simpeg_drivers/plate_simulation/leroi_air/driver.py b/simpeg_drivers/plate_simulation/leroi_air/driver.py index 4e2d4838..dc73d710 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/driver.py +++ b/simpeg_drivers/plate_simulation/leroi_air/driver.py @@ -11,11 +11,7 @@ import subprocess from pathlib import Path -from geoh5py.groups import UIJsonGroup - -from simpeg_drivers.utils.utils import validate_out_group - -from .interface import LeroiAirInterface +from .interface import LeroiAirInput, LeroiAirOutput from .options import LeroiAirOptions @@ -28,14 +24,8 @@ def __init__( ) -> None: """Initialize with simulation options.""" self.options = options - self._interface: LeroiAirInterface | None = None - - @property - def interface(self) -> LeroiAirInterface: - """Lazy-initialized interface for formatting and parsing LeroiAir files.""" - if self._interface is None: - self._interface = LeroiAirInterface(self.options) - return self._interface + self.input: LeroiAirInput = LeroiAirInput(options) + self.output: LeroiAirOutput = LeroiAirOutput(options) @property def project_path(self) -> Path: @@ -61,9 +51,9 @@ def run_leroi(self) -> subprocess.CompletedProcess: def run(self) -> None: """Write input, run LeroiAir, and save simulated data to geoh5.""" - self.interface.write_cfl_file(self.project_path / "LeroiAir.cfl") + self.input.write_cfl_file(self.project_path / "LeroiAir.cfl") self.run_leroi() - self.interface.save_to_geoh5( + self.output.save_to_geoh5( outfile=self.project_path / "LeroiAir.out", out_group=self.options.out_group, ) diff --git a/simpeg_drivers/plate_simulation/leroi_air/interface.py b/simpeg_drivers/plate_simulation/leroi_air/interface.py index 73b8c11c..17472296 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/interface.py +++ b/simpeg_drivers/plate_simulation/leroi_air/interface.py @@ -19,13 +19,12 @@ from .options import LeroiAirOptions -class LeroiAirInterface: - """Interface for running LeroiAir from geoh5py objects.""" +class LeroiAirInput: + """LeroiAir control file formatting from geoh5py objects and options.""" version: str = "8.0" def __init__(self, opts: LeroiAirOptions): - """Initialize with simulation options.""" self.opts = opts @property @@ -237,12 +236,19 @@ def write_cfl_file(self, filepath: Path) -> None: with open(filepath, mode="w", encoding="utf-8") as f: f.write(self.format_cfl_file()) + +class LeroiAirOutput: + """LeroiAir output file parsing and saving to geoh5py objects/data.""" + _COMPONENT_ANCHORS: dict[str, str] = { "crossline": "TRANSVERSE COMPONENT", "inline": "IN-LINE COMPONENT", "vertical": "VERTICAL COMPONENT", } + def __init__(self, opts: LeroiAirOptions): + self.opts = opts + def _find_data_start(self, chunk: list[str]) -> int: """Return the index of the first station data row within a section chunk.""" header_idx = next( diff --git a/tests/plate_simulation/leroi_air/conftest.py b/tests/plate_simulation/leroi_air/conftest.py index 9101c456..8ed2318f 100644 --- a/tests/plate_simulation/leroi_air/conftest.py +++ b/tests/plate_simulation/leroi_air/conftest.py @@ -13,7 +13,10 @@ import pytest from geoh5py import Workspace -from simpeg_drivers.plate_simulation.leroi_air.interface import LeroiAirInterface +from simpeg_drivers.plate_simulation.leroi_air.interface import ( + LeroiAirInput, + LeroiAirOutput, +) from . import generate_plate_options @@ -54,18 +57,18 @@ @pytest.fixture -def mock_interface(): +def mock_input(): opts = MagicMock() - opts.n_stations = N_STATIONS opts.float_precision = 4 - return LeroiAirInterface(opts=opts) + return LeroiAirInput(opts) @pytest.fixture -def real_interface(tmp_path): - with Workspace(tmp_path / "test.geoh5") as geoh5: - opts = generate_plate_options(geoh5) - return LeroiAirInterface(opts=opts) +def mock_output(): + opts = MagicMock() + opts.n_stations = N_STATIONS + opts.float_precision = 4 + return LeroiAirOutput(opts=opts) @pytest.fixture(scope="module") diff --git a/tests/plate_simulation/leroi_air/interface_test.py b/tests/plate_simulation/leroi_air/interface_test.py index 02a9de09..7ce933c6 100644 --- a/tests/plate_simulation/leroi_air/interface_test.py +++ b/tests/plate_simulation/leroi_air/interface_test.py @@ -9,10 +9,12 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' import numpy as np -import pytest from geoh5py import Workspace -from simpeg_drivers.plate_simulation.leroi_air.interface import LeroiAirInterface +from simpeg_drivers.plate_simulation.leroi_air.interface import ( + LeroiAirInput, + LeroiAirOutput, +) from . import generate_plate_options from .conftest import FAKE_OUT, N_CHANNELS, N_STATIONS @@ -22,15 +24,15 @@ def test_line_formatting(tmp_path): with Workspace(tmp_path / "test.geoh5") as geoh5: opts = generate_plate_options(geoh5) - interface = LeroiAirInterface(opts=opts) + leroi_input = LeroiAirInput(opts=opts) - line = interface.format_line(["TXCLN", "CMP", "KPPM"]) + line = leroi_input.format_line(["TXCLN", "CMP", "KPPM"]) assert line == "0.0 3 0 \t ! TXCLN, CMP, KPPM" - line = interface.format_line("TMS") + line = leroi_input.format_line("TMS") assert line == "2.3 2.6 3.2 \t ! TMS" - multi_line = interface.format_multi_line(["TMS", "WIDTH"]) + multi_line = leroi_input.format_multi_line(["TMS", "WIDTH"]) assert multi_line == "2.3 0.3\n2.6 0.3\n3.2 0.6\t ! TMS, WIDTH" @@ -38,8 +40,8 @@ def test_format_cfl_file(tmp_path): with Workspace(tmp_path / "test.geoh5") as geoh5: opts = generate_plate_options(geoh5) - interface = LeroiAirInterface(opts=opts) - cfl_str = interface.format_cfl_file() + leroi_input = LeroiAirInput(opts=opts) + cfl_str = leroi_input.format_cfl_file() # title assert opts.title in cfl_str @@ -109,45 +111,45 @@ def test_format_cfl_file(tmp_path): assert "200.0 100.0 90.0 45.0\t ! PLNGTH, DPWDTH, DZM, DIP" in cfl_str -def test_slice_data_lines_returns_correct_station_rows(mock_interface): +def test_slice_data_lines_returns_correct_station_rows(mock_output): lines = FAKE_OUT.splitlines() - anchor = LeroiAirInterface._COMPONENT_ANCHORS["crossline"] - rows = mock_interface._slice_data_lines(lines, anchor) + anchor = LeroiAirOutput._COMPONENT_ANCHORS["crossline"] + rows = mock_output._slice_data_lines(lines, anchor) assert len(rows) == N_STATIONS assert rows[0].split()[4] == "1.1" assert rows[1].split()[4] == "4.4" -def test_extract_data_returns_channel_columns_only(mock_interface, fake_outfile): - data = mock_interface._extract_data(fake_outfile, component="crossline") +def test_extract_data_returns_channel_columns_only(mock_output, fake_outfile): + data = mock_output._extract_data(fake_outfile, component="crossline") assert data.shape == (N_STATIONS, N_CHANNELS) -def test_extract_data_parses_transverse_component_values(mock_interface, fake_outfile): - data = mock_interface._extract_data(fake_outfile, component="crossline") +def test_extract_data_parses_transverse_component_values(mock_output, fake_outfile): + data = mock_output._extract_data(fake_outfile, component="crossline") np.testing.assert_array_almost_equal(data[0], [1.1, 2.2, 3.3]) np.testing.assert_array_almost_equal(data[1], [4.4, 5.5, 6.6]) -def test_extract_data_parses_inline_component_values(mock_interface, fake_outfile): - data = mock_interface._extract_data(fake_outfile, component="inline") +def test_extract_data_parses_inline_component_values(mock_output, fake_outfile): + data = mock_output._extract_data(fake_outfile, component="inline") np.testing.assert_array_almost_equal(data[0], [7.7, 8.8, 9.9]) np.testing.assert_array_almost_equal(data[1], [10.0, 11.1, 12.2]) -def test_extract_data_parses_vertical_component_values(mock_interface, fake_outfile): - data = mock_interface._extract_data(fake_outfile, component="vertical") +def test_extract_data_parses_vertical_component_values(mock_output, fake_outfile): + data = mock_output._extract_data(fake_outfile, component="vertical") np.testing.assert_array_almost_equal(data[0], [13.3, 14.4, 15.5]) np.testing.assert_array_almost_equal(data[1], [16.6, 17.7, 18.8]) -def test_format_value_int(mock_interface): - assert mock_interface._format_value(3) == "3" - assert mock_interface._format_value(np.int64(7)) == "7" +def test_format_value_int(mock_input): + assert mock_input._format_value(3) == "3" + assert mock_input._format_value(np.int64(7)) == "7" -def test_format_value_float(mock_interface): - assert mock_interface._format_value(1.5) == "1.5" - assert mock_interface._format_value(1.123456789) == "1.1235" - mock_interface.opts.float_precision = 2 - assert mock_interface._format_value(3.141592653) == "3.14" +def test_format_value_float(mock_input): + assert mock_input._format_value(1.5) == "1.5" + assert mock_input._format_value(1.123456789) == "1.1235" + mock_input.opts.float_precision = 2 + assert mock_input._format_value(3.141592653) == "3.14" From b66248e499a2b74e835d60c911ccb94c1b3f66db Mon Sep 17 00:00:00 2001 From: benjamink Date: Fri, 24 Apr 2026 13:30:43 -0700 Subject: [PATCH 28/43] refactor leroi options to limit scope of options passed to output half of interface class --- .../plate_simulation/leroi_air/driver.py | 11 +- .../plate_simulation/leroi_air/interface.py | 42 ++-- .../plate_simulation/leroi_air/options.py | 209 +++++++++--------- tests/plate_simulation/leroi_air/__init__.py | 7 +- .../leroi_air/options_test.py | 24 +- 5 files changed, 157 insertions(+), 136 deletions(-) diff --git a/simpeg_drivers/plate_simulation/leroi_air/driver.py b/simpeg_drivers/plate_simulation/leroi_air/driver.py index dc73d710..056d9c92 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/driver.py +++ b/simpeg_drivers/plate_simulation/leroi_air/driver.py @@ -11,7 +11,7 @@ import subprocess from pathlib import Path -from .interface import LeroiAirInput, LeroiAirOutput +from .interface import LeroiAirInterface from .options import LeroiAirOptions @@ -24,13 +24,12 @@ def __init__( ) -> None: """Initialize with simulation options.""" self.options = options - self.input: LeroiAirInput = LeroiAirInput(options) - self.output: LeroiAirOutput = LeroiAirOutput(options) + self.interface = LeroiAirInterface(options) @property def project_path(self) -> Path: """Directory containing the geoh5 workspace file.""" - return self.options.survey.workspace.h5file.parent + return self.options.survey.entity.workspace.h5file.parent def run_leroi(self) -> subprocess.CompletedProcess: """Run the LeroiAir executable and raise on non-zero exit.""" @@ -51,9 +50,9 @@ def run_leroi(self) -> subprocess.CompletedProcess: def run(self) -> None: """Write input, run LeroiAir, and save simulated data to geoh5.""" - self.input.write_cfl_file(self.project_path / "LeroiAir.cfl") + self.interface.input.write_cfl_file(self.project_path / "LeroiAir.cfl") self.run_leroi() - self.output.save_to_geoh5( + self.interface.output.save_to_geoh5( outfile=self.project_path / "LeroiAir.out", out_group=self.options.out_group, ) diff --git a/simpeg_drivers/plate_simulation/leroi_air/interface.py b/simpeg_drivers/plate_simulation/leroi_air/interface.py index 17472296..c3ad9e42 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/interface.py +++ b/simpeg_drivers/plate_simulation/leroi_air/interface.py @@ -9,14 +9,20 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' from __future__ import annotations +from functools import cached_property from pathlib import Path from typing import Any, Literal import numpy as np -from geoh5py.groups import UIJsonGroup from geoh5py.shared.utils import fetch_active_workspace -from .options import LeroiAirOptions +from .options import LeroiAirOptions, SurveyOptions + + +class LeroiAirInterface: + def __init__(self, opts: LeroiAirOptions): + self.input = LeroiAirInput(opts) + self.output = LeroiAirOutput(opts.survey) class LeroiAirInput: @@ -27,7 +33,7 @@ class LeroiAirInput: def __init__(self, opts: LeroiAirOptions): self.opts = opts - @property + @cached_property def aliased_values(self) -> dict[str, Any]: """Serves .cfl input file aliases and corresponding data to line formatter.""" return { @@ -36,17 +42,17 @@ def aliased_values(self) -> dict[str, Any]: "PRFL": 1, "ISTOP": 0, "ISW": 1, - "NSX": len(self.opts.ontime_waveform), + "NSX": len(self.opts.survey.ontime_waveform), "STEP": 0 if self.opts.magnetic_field == "dBdt" else 1, "UNITS": 1, - "NCHNL": len(self.opts.channels), + "NCHNL": len(self.opts.survey.channels), "KRXW": 2, - "REFTYM": self.opts.timing_mark, - "OFFTIME": self.opts.offtime, - "TXON": self.opts.ontime_waveform[:, 0], - "TXAMP": self.opts.ontime_waveform[:, 1], - "TMS": self.opts.timing_mark + np.array(self.opts.channels), - "WIDTH": self.opts.channel_widths, + "REFTYM": self.opts.survey.timing_mark, + "OFFTIME": self.opts.survey.offtime, + "TXON": self.opts.survey.ontime_waveform[:, 0], + "TXAMP": self.opts.survey.ontime_waveform[:, 1], + "TMS": self.opts.survey.timing_mark + np.array(self.opts.survey.channels), + "WIDTH": self.opts.survey.channel_widths, "TXCLN": 0.0, "CMP": 3, "KPPM": 0, @@ -56,13 +62,13 @@ def aliased_values(self) -> dict[str, Any]: "ZRX0": 0.0, "XRX0": 0.0, "YRX0": 0.0, - "NSTAT": self.opts.n_stations, + "NSTAT": self.opts.survey.n_stations, "SURVEY": 2, "BAROMTRC": 1, "LINE_TAG": 0, - "EAST": self.opts.locations[:, 0], - "NORTH": self.opts.locations[:, 1], - "ALT": self.opts.drape_height, + "EAST": self.opts.survey.locations[:, 0], + "NORTH": self.opts.survey.locations[:, 1], + "ALT": self.opts.survey.drape_height(self.opts.topo), "NLAYER": self.opts.n_layers, "NPLATE": self.opts.n_plates, "NLITH": self.opts.n_layers + self.opts.n_plates, @@ -246,7 +252,7 @@ class LeroiAirOutput: "vertical": "VERTICAL COMPONENT", } - def __init__(self, opts: LeroiAirOptions): + def __init__(self, opts: SurveyOptions): self.opts = opts def _find_data_start(self, chunk: list[str]) -> int: @@ -277,11 +283,11 @@ def save_to_geoh5(self, outfile: str | Path, out_group): """Save LeroiAir simulated data on the provided survey to geoh5.""" with out_group.workspace.open(mode="r+"): - survey = self.opts.survey.copy(parent=out_group, copy_children=False) + survey = self.opts.entity.copy(parent=out_group, copy_children=False) for component in "inline", "crossline", "vertical": data = self._extract_data(outfile=outfile, component=component) - with fetch_active_workspace(self.opts.survey.workspace, mode="r+"): + with fetch_active_workspace(self.opts.entity.workspace, mode="r+"): entities = survey.add_data( { f"fwd {component} [{i}]": {"values": data[:, i]} diff --git a/simpeg_drivers/plate_simulation/leroi_air/options.py b/simpeg_drivers/plate_simulation/leroi_air/options.py index b9b55d1e..a3029270 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/options.py +++ b/simpeg_drivers/plate_simulation/leroi_air/options.py @@ -23,12 +23,120 @@ from simpeg_drivers.plate_simulation.options import ModelOptions +class SurveyOptions(BaseModel): + model_config = ConfigDict(arbitrary_types_allowed=True) + + entity: AirborneTEMReceivers + + @property + def waveform(self) -> np.ndarray: + """Survey transmitter waveform.""" + return self.entity.waveform + + @property + def frequency(self) -> float: + """ + Transmitter frequency. + + Can be set by the user in the metadata. If not set, it will be assumed + that the waveform property contains a full halfcycle and the frequency + will be calculated as the reciprocal of the provided waveform time span. + """ + frequency = self.entity.metadata.get("frequency", None) + if frequency is None: + half_cycle_time = self.waveform[-1, 0] - self.waveform[0, 0] + frequency = 1 / (2 * half_cycle_time) + + return frequency + + @property + def channels(self) -> np.ndarray: + """Time channel midpoints referenced from the timing_mark.""" + return self.entity.channels + + @property + def channel_widths(self) -> np.ndarray: + """Time channel widths.""" + channel_widths = self.entity.metadata.get("channel_widths", None) + if channel_widths is None: + channel_widths = np.diff(np.r_[0, self.channels]) + return channel_widths + + @property + def timing_mark(self) -> float: + """Reference point for timing of the channels.""" + return self.entity.timing_mark + + @property + def units(self) -> str: + """Units of the time channels.""" + return self.entity.unit + + @property + def _ontime(self) -> float: + """Time at which the transmitter current turns off.""" + return float(self.waveform[self._offtime_mask(), 0][0]) + + @property + def offtime(self) -> float: + """ + Time at which the transmitter current is zero. + + This offtime is based on system frequency and may include times that + are not accounted for by the waveform. + """ + half_cycle = 1 / (2 * self.frequency) + zero_current_ind = np.where(self.waveform[:, 1] == 0)[0] + first_zero = 1 if self.waveform[0, 1] == 0.0 else 0 + ontime = self.waveform[zero_current_ind][first_zero, 0] + + return half_cycle - ontime + + @property + def ontime_waveform(self) -> np.ndarray: + """On-time waveform including leading and trailing 0 current times.""" + ontime_waveform = self.waveform[~self._offtime_mask(), :] + endpoint = self.waveform[self._offtime_mask()][0, :] + return np.vstack([ontime_waveform, endpoint]) + + @property + def locations(self) -> np.ndarray: + """Survey receiver locations.""" + return self.entity.locations + + @property + def n_stations(self) -> int: + """Number of survey stations at which time channel data will be simulated.""" + return len(self.locations) + + def drape_height(self, topo: np.ndarray) -> np.ndarray: + """ + Survey height over topography. + + :param topo: Topography array of x, y, elevation. + + :returns Array of survey height over topography. + """ + + survey_locs = self.entity.locations.copy() + topo_interp = LinearNDInterpolator(topo[:, :2], topo[:, 2]) + topo_at_survey_locations = topo_interp(survey_locs[:, 0], survey_locs[:, 1]) + + return survey_locs[:, 2] - topo_at_survey_locations + + def _offtime_mask(self) -> np.ndarray: + """Mask selecting off-time rows from the waveform array.""" + mask = np.isclose(self.waveform[:, 1], 0) + mask[0] = False + return mask + + class LeroiAirOptions(BaseModel): """Configuration for a LeroiAir airborne TEM forward simulation.""" model_config = ConfigDict(arbitrary_types_allowed=True) - survey: AirborneTEMReceivers + survey: SurveyOptions topo: np.ndarray layer_resistivities: list[float] layer_thicknesses: list[float] @@ -77,7 +185,7 @@ def from_plate_simulation_options( topo_xyz = InversionTopography(survey.workspace, simulation).locations return cls( - survey=survey, + survey=SurveyOptions(entity=survey), topo=topo_xyz, layer_resistivities=[ model.overburden_options.overburden_property, @@ -95,26 +203,6 @@ def title(self) -> str: """Provides a generic title for all LeroiAir simulations.""" return "LeroiAir modelling for plate-simulation package." - @property - def locations(self) -> np.ndarray: - """Survey receiver locations.""" - return self.survey.locations - - @property - def drape_height(self) -> np.ndarray: - """Survey height over topography.""" - - survey_locs = self.locations.copy() - topo_interp = LinearNDInterpolator(self.topo[:, :2], self.topo[:, 2]) - topo_at_survey_locations = topo_interp(survey_locs[:, 0], survey_locs[:, 1]) - - return survey_locs[:, 2] - topo_at_survey_locations - - @property - def n_stations(self) -> int: - """Number of survey stations at which time channel data will be simulated.""" - return len(self.locations) - @property def n_layers(self) -> int: """Number of background layers.""" @@ -125,79 +213,8 @@ def n_plates(self) -> int: """Number of plates.""" return len(self.plate_geometries) - @property - def waveform(self) -> np.ndarray: - """Survey transmitter waveform.""" - return self.survey.waveform - - @property - def frequency(self) -> float: - """ - Transmitter frequency. - - Can be set by the user in the metadata. If not set, it will be assumed - that the waveform property contains a full halfcycle and the frequency - will be calculated as the reciprocal of the provided waveform time span. - """ - frequency = self.survey.metadata.get("frequency", None) - if frequency is None: - half_cycle_time = self.waveform[-1, 0] - self.waveform[0, 0] - frequency = 1 / (2 * half_cycle_time) - - return frequency - # TODO use units to convert time to milliseconds used by leroi. - @property - def channels(self) -> np.ndarray: - """Time channel midpoints referenced from the timing_mark.""" - return self.survey.channels - - @property - def channel_widths(self) -> np.ndarray: - """Time channel widths.""" - channel_widths = self.survey.metadata.get("channel_widths", None) - if channel_widths is None: - channel_widths = np.diff(np.r_[0, self.channels]) - return channel_widths - - @property - def timing_mark(self) -> float: - """Reference point for timing of the channels.""" - return self.survey.timing_mark - - @property - def units(self) -> str: - """Units of the time channels.""" - return self.survey.unit - - @property - def _ontime(self) -> float: - """Time at which the transmitter current turns off.""" - return float(self.waveform[self._offtime_mask(), 0][0]) - - @property - def offtime(self) -> float: - """ - Time at which the transmitter current is zero. - - This offtime is based on system frequency and may include times that - are not accounted for by the waveform. - """ - half_cycle = 1 / (2 * self.frequency) - zero_current_ind = np.where(self.waveform[:, 1] == 0)[0] - first_zero = 1 if self.waveform[0, 1] == 0.0 else 0 - ontime = self.waveform[zero_current_ind][first_zero, 0] - - return half_cycle - ontime - - @property - def ontime_waveform(self) -> np.ndarray: - """On-time waveform including leading and trailing 0 current times.""" - ontime_waveform = self.waveform[~self._offtime_mask(), :] - endpoint = self.waveform[self._offtime_mask()][0, :] - return np.vstack([ontime_waveform, endpoint]) - @property def resistivities(self) -> np.ndarray: """All resistivities.""" @@ -213,9 +230,3 @@ def conductivity_thicknesses(self) -> np.ndarray: sigma.append([g.width for g in self.plate_geometries] * plate_conductivities) return np.hstack(sigma) - - def _offtime_mask(self) -> np.ndarray: - """Mask selecting off-time rows from the waveform array.""" - mask = np.isclose(self.waveform[:, 1], 0) - mask[0] = False - return mask diff --git a/tests/plate_simulation/leroi_air/__init__.py b/tests/plate_simulation/leroi_air/__init__.py index bebd5700..b938fdce 100644 --- a/tests/plate_simulation/leroi_air/__init__.py +++ b/tests/plate_simulation/leroi_air/__init__.py @@ -12,7 +12,10 @@ from geoapps_utils.modelling.plates import PlateModel from geoh5py.shared.utils import fetch_active_workspace -from simpeg_drivers.plate_simulation.leroi_air.options import LeroiAirOptions +from simpeg_drivers.plate_simulation.leroi_air.options import ( + LeroiAirOptions, + SurveyOptions, +) from simpeg_drivers.utils.synthetics.surveys.time_domain.airborne_tdem import ( generate_airborne_tdem_survey, ) @@ -40,7 +43,7 @@ def generate_plate_options(workspace): ] topo = np.column_stack([X.flatten(), Y.flatten(), np.zeros(X.size)]) opts = LeroiAirOptions( - survey=survey, + survey=SurveyOptions(entity=survey), topo=topo, layer_resistivities=layer_resistivities, layer_thicknesses=layer_thicknesses, diff --git a/tests/plate_simulation/leroi_air/options_test.py b/tests/plate_simulation/leroi_air/options_test.py index c2ee42cb..baa3461e 100644 --- a/tests/plate_simulation/leroi_air/options_test.py +++ b/tests/plate_simulation/leroi_air/options_test.py @@ -18,40 +18,42 @@ def test_options_channel_widths(plate_options): - assert np.allclose(plate_options.channel_widths, [0.3, 0.3, 0.6]) + assert np.allclose(plate_options.survey.channel_widths, [0.3, 0.3, 0.6]) def test_options_channel_widths_from_metadata(tmp_path): explicit_widths = [0.1, 0.2, 0.3] with Workspace(tmp_path / "test.geoh5") as geoh5: opts = generate_plate_options(workspace=geoh5) - meta = opts.survey.metadata + meta = opts.survey.entity.metadata meta["channel_widths"] = explicit_widths - opts.survey.metadata = meta + opts.survey.entity.metadata = meta - assert np.allclose(opts.channel_widths, explicit_widths) + assert np.allclose(opts.survey.channel_widths, explicit_widths) def test_options_frequency(plate_options): - assert np.isclose(plate_options.frequency, 0.142857143) + assert np.isclose(plate_options.survey.frequency, 0.142857143) def test_options_frequency_from_metadata(tmp_path): with Workspace(tmp_path / "test.geoh5") as geoh5: opts = generate_plate_options(workspace=geoh5) - meta = opts.survey.metadata + meta = opts.survey.entity.metadata meta["frequency"] = 0.15 - opts.survey.metadata = meta + opts.survey.entity.metadata = meta - assert np.isclose(opts.frequency, 0.15) + assert np.isclose(opts.survey.frequency, 0.15) def test_options_offtime(plate_options): - assert np.isclose(plate_options.offtime, 1.6) + assert np.isclose(plate_options.survey.offtime, 1.6) def test_options_ontime_waveform(plate_options): - assert np.allclose(plate_options.ontime_waveform, plate_options.waveform[0:8, :]) + assert np.allclose( + plate_options.survey.ontime_waveform, plate_options.survey.waveform[0:8, :] + ) def test_options_conductivity_thicknesses(plate_options): @@ -65,7 +67,7 @@ def test_options_resistivities(plate_options): def test_options_drape_height(plate_options): - assert np.allclose(plate_options.drape_height, 20.0) + assert np.allclose(plate_options.survey.drape_height(plate_options.topo), 20.0) def test_options_validate_layer_lengths_mismatch(plate_options): From e8603493265493864aa65359b86dfbc626805a3c Mon Sep 17 00:00:00 2001 From: benjamink Date: Fri, 24 Apr 2026 13:40:49 -0700 Subject: [PATCH 29/43] title as field --- simpeg_drivers/plate_simulation/leroi_air/options.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/simpeg_drivers/plate_simulation/leroi_air/options.py b/simpeg_drivers/plate_simulation/leroi_air/options.py index a3029270..98c24df1 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/options.py +++ b/simpeg_drivers/plate_simulation/leroi_air/options.py @@ -136,6 +136,7 @@ class LeroiAirOptions(BaseModel): model_config = ConfigDict(arbitrary_types_allowed=True) + title: str = "LeroiAir modelling for plate-simulation package." survey: SurveyOptions topo: np.ndarray layer_resistivities: list[float] @@ -198,11 +199,6 @@ def from_plate_simulation_options( out_group=simulation.out_group, ) - @property - def title(self) -> str: - """Provides a generic title for all LeroiAir simulations.""" - return "LeroiAir modelling for plate-simulation package." - @property def n_layers(self) -> int: """Number of background layers.""" From 337b00be1ef44ad2e44c20af07a22cd49ea23b4f Mon Sep 17 00:00:00 2001 From: benjamink Date: Fri, 24 Apr 2026 14:13:28 -0700 Subject: [PATCH 30/43] clean up all the fetch_active_workspace calls --- simpeg_drivers/plate_simulation/driver.py | 43 ++++++++----------- .../plate_simulation/leroi_air/interface.py | 25 +++++------ 2 files changed, 30 insertions(+), 38 deletions(-) diff --git a/simpeg_drivers/plate_simulation/driver.py b/simpeg_drivers/plate_simulation/driver.py index 6d2f3cad..a773aff2 100644 --- a/simpeg_drivers/plate_simulation/driver.py +++ b/simpeg_drivers/plate_simulation/driver.py @@ -79,11 +79,9 @@ def __init__( def run(self) -> InversionDriver: """Create octree mesh, fill model, and simulate.""" - with fetch_active_workspace(self.simulation_parameters.geoh5, mode="r+"): - self.simulation_driver.run() - self.simulation_parameters.update_out_group_options() - with fetch_active_workspace(self.params.geoh5, mode="r+"): - self.update_monitoring_directory(self._out_group) + self.simulation_driver.run() + self.simulation_parameters.update_out_group_options() + self.update_monitoring_directory(self._out_group) logger.info("done.") logger.handlers.clear() @@ -132,9 +130,8 @@ def plates(self) -> list[Plate]: self.params.model.plate_options.spacing, self.params.model.plate_options.geometry.direction, ) - with fetch_active_workspace(self.params.geoh5, mode="r+") as geoh5: - for plate in self._plates: - plate.to_maxwell_plate(geoh5, parent=self._out_group) + for plate in self._plates: + plate.to_maxwell_plate(self.params.geoh5, parent=self._out_group) return self._plates @@ -146,15 +143,14 @@ def make_mesh(self) -> Octree: """ logger.info("making the mesh...") - with fetch_active_workspace(self.params.geoh5, mode="r+") as geoh5: - surfaces = [p.surface(geoh5) for p in self.plates] - self._mesh = get_octree_mesh( - opts=self.params.mesh, - survey=self.survey, - topography=self.simulation_parameters.active_cells.topography_object, - plates=surfaces, - name="Octree", - ) + surfaces = [p.surface(self.params.geoh5) for p in self.plates] + self._mesh = get_octree_mesh( + opts=self.params.mesh, + survey=self.survey, + topography=self.simulation_parameters.active_cells.topography_object, + plates=surfaces, + name="Octree", + ) self._mesh.parent = self._out_group return self._mesh @@ -327,13 +323,12 @@ def _initialize_forward_opts(self) -> BaseForwardOptions: models_update["starting_model"] = None update["models"] = opts.models.model_copy(update=models_update) - with fetch_active_workspace(self.params.geoh5, mode="r+"): - out_group = validate_out_group(opts) - out_group = out_group.copy( - parent=self.out_group, - copy_children=False, - copy_relatives=False, - ) + out_group = validate_out_group(opts) + out_group = out_group.copy( + parent=self.out_group, + copy_children=False, + copy_relatives=False, + ) update["out_group"] = out_group forward_opts = opts.model_copy(update=update) forward_opts.update_out_group_options() diff --git a/simpeg_drivers/plate_simulation/leroi_air/interface.py b/simpeg_drivers/plate_simulation/leroi_air/interface.py index c3ad9e42..550f99ce 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/interface.py +++ b/simpeg_drivers/plate_simulation/leroi_air/interface.py @@ -282,17 +282,14 @@ def _extract_data( def save_to_geoh5(self, outfile: str | Path, out_group): """Save LeroiAir simulated data on the provided survey to geoh5.""" - with out_group.workspace.open(mode="r+"): - survey = self.opts.entity.copy(parent=out_group, copy_children=False) - for component in "inline", "crossline", "vertical": - data = self._extract_data(outfile=outfile, component=component) - - with fetch_active_workspace(self.opts.entity.workspace, mode="r+"): - entities = survey.add_data( - { - f"fwd {component} [{i}]": {"values": data[:, i]} - for i in range(len(self.opts.channels)) - } - ) - - survey.create_property_group(name=component, properties=entities) + survey = self.opts.entity.copy(parent=out_group, copy_children=False) + for component in "inline", "crossline", "vertical": + data = self._extract_data(outfile=outfile, component=component) + entities = survey.add_data( + { + f"fwd {component} [{i}]": {"values": data[:, i]} + for i in range(len(self.opts.channels)) + } + ) + + survey.create_property_group(name=component, properties=entities) From 5d3dbea127017c73944c7f5f584751d0c6ca6b6c Mon Sep 17 00:00:00 2001 From: benjamink Date: Mon, 27 Apr 2026 13:49:08 -0700 Subject: [PATCH 31/43] use leroi -> UseLeroiAir --- simpeg_drivers-assets/uijson/plate_simulation.ui.json | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/simpeg_drivers-assets/uijson/plate_simulation.ui.json b/simpeg_drivers-assets/uijson/plate_simulation.ui.json index 45a2fc84..dbbf357b 100644 --- a/simpeg_drivers-assets/uijson/plate_simulation.ui.json +++ b/simpeg_drivers-assets/uijson/plate_simulation.ui.json @@ -20,7 +20,7 @@ }, "use_leroi": { "main": true, - "label": "Use leroi", + "label": "Use LeroiAir", "value": false, "tooltip": "If checked, plate simulation will use LeroiAir to simulate the plate response." }, From 2587551201e1d5c1d71c024142ebeb05ff07f4c3 Mon Sep 17 00:00:00 2001 From: benjamink Date: Mon, 27 Apr 2026 13:51:03 -0700 Subject: [PATCH 32/43] fix copyrights --- simpeg_drivers/plate_simulation/leroi_air/driver.py | 2 +- simpeg_drivers/plate_simulation/leroi_air/interface.py | 2 +- simpeg_drivers/plate_simulation/leroi_air/options.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/simpeg_drivers/plate_simulation/leroi_air/driver.py b/simpeg_drivers/plate_simulation/leroi_air/driver.py index 056d9c92..1539e7a7 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/driver.py +++ b/simpeg_drivers/plate_simulation/leroi_air/driver.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2026 Mira Geoscience Ltd. ' +# Copyright (c) 2023-2026 Mira Geoscience Ltd. ' # ' # This file is part of simpeg-drivers package. ' # ' diff --git a/simpeg_drivers/plate_simulation/leroi_air/interface.py b/simpeg_drivers/plate_simulation/leroi_air/interface.py index 550f99ce..bd66e4ba 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/interface.py +++ b/simpeg_drivers/plate_simulation/leroi_air/interface.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2026 Mira Geoscience Ltd. ' +# Copyright (c) 2023-2026 Mira Geoscience Ltd. ' # ' # This file is part of simpeg-drivers package. ' # ' diff --git a/simpeg_drivers/plate_simulation/leroi_air/options.py b/simpeg_drivers/plate_simulation/leroi_air/options.py index 98c24df1..751e9667 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/options.py +++ b/simpeg_drivers/plate_simulation/leroi_air/options.py @@ -1,5 +1,5 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2026 Mira Geoscience Ltd. ' +# Copyright (c) 2023-2026 Mira Geoscience Ltd. ' # ' # This file is part of simpeg-drivers package. ' # ' From 3bb7798096c95a25162c39c978884faacbb88807 Mon Sep 17 00:00:00 2001 From: benjamink Date: Tue, 28 Apr 2026 14:49:07 -0700 Subject: [PATCH 33/43] try to bypass github run of leroi test only --- tests/plate_simulation/runtest/leroi_test.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/tests/plate_simulation/runtest/leroi_test.py b/tests/plate_simulation/runtest/leroi_test.py index c65374fe..3da33648 100644 --- a/tests/plate_simulation/runtest/leroi_test.py +++ b/tests/plate_simulation/runtest/leroi_test.py @@ -8,6 +8,9 @@ # ' # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +import os + +import pytest from geoapps_utils.modelling.plates import PlateModel from geoh5py import Workspace from geoh5py.groups import SimPEGGroup @@ -37,6 +40,13 @@ from tests.utils.targets import get_workspace +IN_GITHUB_ACTIONS = os.getenv("GITHUB_ACTIONS") == "true" + + +@pytest.mark.skipif( + IN_GITHUB_ACTIONS, + reason="Test requires a LeroiAir550_JR.exe in the environments Scripts directory that is not stored in GitHub.", +) def test_leroi_run(tmp_path): with get_workspace(tmp_path / "leroi_test.geoh5") as geoh5: From 3241db350ea8d8cb2a522dd5799dcaf2f2dfba8c Mon Sep 17 00:00:00 2001 From: benjamink Date: Tue, 28 Apr 2026 15:11:23 -0700 Subject: [PATCH 34/43] remove survey derived properties in favor of direct access --- .../plate_simulation/leroi_air/interface.py | 13 +++--- .../plate_simulation/leroi_air/options.py | 45 +++++-------------- .../leroi_air/options_test.py | 3 +- 3 files changed, 19 insertions(+), 42 deletions(-) diff --git a/simpeg_drivers/plate_simulation/leroi_air/interface.py b/simpeg_drivers/plate_simulation/leroi_air/interface.py index bd66e4ba..3ca11d5c 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/interface.py +++ b/simpeg_drivers/plate_simulation/leroi_air/interface.py @@ -45,13 +45,14 @@ def aliased_values(self) -> dict[str, Any]: "NSX": len(self.opts.survey.ontime_waveform), "STEP": 0 if self.opts.magnetic_field == "dBdt" else 1, "UNITS": 1, - "NCHNL": len(self.opts.survey.channels), + "NCHNL": len(self.opts.survey.entity.channels), "KRXW": 2, - "REFTYM": self.opts.survey.timing_mark, + "REFTYM": self.opts.survey.entity.timing_mark, "OFFTIME": self.opts.survey.offtime, "TXON": self.opts.survey.ontime_waveform[:, 0], "TXAMP": self.opts.survey.ontime_waveform[:, 1], - "TMS": self.opts.survey.timing_mark + np.array(self.opts.survey.channels), + "TMS": self.opts.survey.entity.timing_mark + + np.array(self.opts.survey.entity.channels), "WIDTH": self.opts.survey.channel_widths, "TXCLN": 0.0, "CMP": 3, @@ -66,8 +67,8 @@ def aliased_values(self) -> dict[str, Any]: "SURVEY": 2, "BAROMTRC": 1, "LINE_TAG": 0, - "EAST": self.opts.survey.locations[:, 0], - "NORTH": self.opts.survey.locations[:, 1], + "EAST": self.opts.survey.entity.locations[:, 0], + "NORTH": self.opts.survey.entity.locations[:, 1], "ALT": self.opts.survey.drape_height(self.opts.topo), "NLAYER": self.opts.n_layers, "NPLATE": self.opts.n_plates, @@ -288,7 +289,7 @@ def save_to_geoh5(self, outfile: str | Path, out_group): entities = survey.add_data( { f"fwd {component} [{i}]": {"values": data[:, i]} - for i in range(len(self.opts.channels)) + for i in range(len(self.opts.entity.channels)) } ) diff --git a/simpeg_drivers/plate_simulation/leroi_air/options.py b/simpeg_drivers/plate_simulation/leroi_air/options.py index 751e9667..33f6e3c5 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/options.py +++ b/simpeg_drivers/plate_simulation/leroi_air/options.py @@ -28,11 +28,6 @@ class SurveyOptions(BaseModel): entity: AirborneTEMReceivers - @property - def waveform(self) -> np.ndarray: - """Survey transmitter waveform.""" - return self.entity.waveform - @property def frequency(self) -> float: """ @@ -44,38 +39,23 @@ def frequency(self) -> float: """ frequency = self.entity.metadata.get("frequency", None) if frequency is None: - half_cycle_time = self.waveform[-1, 0] - self.waveform[0, 0] + half_cycle_time = self.entity.waveform[-1, 0] - self.entity.waveform[0, 0] frequency = 1 / (2 * half_cycle_time) return frequency - @property - def channels(self) -> np.ndarray: - """Time channel midpoints referenced from the timing_mark.""" - return self.entity.channels - @property def channel_widths(self) -> np.ndarray: """Time channel widths.""" channel_widths = self.entity.metadata.get("channel_widths", None) if channel_widths is None: - channel_widths = np.diff(np.r_[0, self.channels]) + channel_widths = np.diff(np.r_[0, self.entity.channels]) return channel_widths - @property - def timing_mark(self) -> float: - """Reference point for timing of the channels.""" - return self.entity.timing_mark - - @property - def units(self) -> str: - """Units of the time channels.""" - return self.entity.unit - @property def _ontime(self) -> float: """Time at which the transmitter current turns off.""" - return float(self.waveform[self._offtime_mask(), 0][0]) + return float(self.entity.waveform[self._offtime_mask(), 0][0]) @property def offtime(self) -> float: @@ -86,28 +66,23 @@ def offtime(self) -> float: are not accounted for by the waveform. """ half_cycle = 1 / (2 * self.frequency) - zero_current_ind = np.where(self.waveform[:, 1] == 0)[0] - first_zero = 1 if self.waveform[0, 1] == 0.0 else 0 - ontime = self.waveform[zero_current_ind][first_zero, 0] + zero_current_ind = np.where(self.entity.waveform[:, 1] == 0)[0] + first_zero = 1 if self.entity.waveform[0, 1] == 0.0 else 0 + ontime = self.entity.waveform[zero_current_ind][first_zero, 0] return half_cycle - ontime @property def ontime_waveform(self) -> np.ndarray: """On-time waveform including leading and trailing 0 current times.""" - ontime_waveform = self.waveform[~self._offtime_mask(), :] - endpoint = self.waveform[self._offtime_mask()][0, :] + ontime_waveform = self.entity.waveform[~self._offtime_mask(), :] + endpoint = self.entity.waveform[self._offtime_mask()][0, :] return np.vstack([ontime_waveform, endpoint]) - @property - def locations(self) -> np.ndarray: - """Survey receiver locations.""" - return self.entity.locations - @property def n_stations(self) -> int: """Number of survey stations at which time channel data will be simulated.""" - return len(self.locations) + return len(self.entity.locations) def drape_height(self, topo: np.ndarray) -> np.ndarray: """ @@ -126,7 +101,7 @@ def drape_height(self, topo: np.ndarray) -> np.ndarray: def _offtime_mask(self) -> np.ndarray: """Mask selecting off-time rows from the waveform array.""" - mask = np.isclose(self.waveform[:, 1], 0) + mask = np.isclose(self.entity.waveform[:, 1], 0) mask[0] = False return mask diff --git a/tests/plate_simulation/leroi_air/options_test.py b/tests/plate_simulation/leroi_air/options_test.py index baa3461e..c400fced 100644 --- a/tests/plate_simulation/leroi_air/options_test.py +++ b/tests/plate_simulation/leroi_air/options_test.py @@ -52,7 +52,8 @@ def test_options_offtime(plate_options): def test_options_ontime_waveform(plate_options): assert np.allclose( - plate_options.survey.ontime_waveform, plate_options.survey.waveform[0:8, :] + plate_options.survey.ontime_waveform, + plate_options.survey.entity.waveform[0:8, :], ) From f6f8a8fc089ca8d6b0728cdc12faa015d44d6d3b Mon Sep 17 00:00:00 2001 From: benjamink Date: Wed, 29 Apr 2026 09:17:24 -0700 Subject: [PATCH 35/43] improve docstrings --- .../plate_simulation/leroi_air/interface.py | 70 ++++++++++++++++--- 1 file changed, 59 insertions(+), 11 deletions(-) diff --git a/simpeg_drivers/plate_simulation/leroi_air/interface.py b/simpeg_drivers/plate_simulation/leroi_air/interface.py index 3ca11d5c..db175e0c 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/interface.py +++ b/simpeg_drivers/plate_simulation/leroi_air/interface.py @@ -96,7 +96,11 @@ def aliased_values(self) -> dict[str, Any]: } def _format_value(self, value: int | float) -> str: - """Format a scalar value, preserving integer formatting and applying float precision only when needed.""" + """ + Format a scalar value as an integer of bounded precision float. + + :param value: Value to format. + """ match value: case int() | np.integer(): return str(int(value)) @@ -106,30 +110,50 @@ def _format_value(self, value: int | float) -> str: return str(value) def _format_float(self, value: float) -> str: - """Format a float, truncating to float_precision only when needed.""" + """ + Format a float, truncating to float_precision only when needed. + + :param value: Value to format as a bounded precision float. + """ _, _, decimals = str(value).partition(".") if len(decimals) > self.opts.float_precision: return f"{value:.{self.opts.float_precision}f}" return str(value) def _format_scalar_params(self, params: list[str]) -> str: - """Format one scalar value per param onto a single line.""" + """ + Format one scalar value per param onto a single line. + + :param params: Parameter names to format a single line of scalar data. + """ values = [self._format_value(self.aliased_values[k]) for k in params] return f"{' '.join(values)} \t ! {', '.join(params)}" def _format_vector_param(self, param: str) -> str: - """Format all elements of a single vector param onto a single line.""" + """ + Format all elements of a single vector param onto a single line. + + :param param: Parameter name to format as a single line of vector data. + """ values = [self._format_value(v) for v in self.aliased_values[param]] return f"{' '.join(values)} \t ! {param}" def format_line(self, params: str | list[str]) -> str: - """Format one or more param values on a single line.""" + """ + Format one or more param values on a single line. + + :param params: Parameter names to format as a single line. + """ if isinstance(params, str): return self._format_vector_param(params) return self._format_scalar_params(params) def format_multi_line(self, params: list[str]) -> str: - """Format one or more vector param values as a row-per-entry table.""" + """ + Format one or more vector param values as a row-per-entry table. + + :param params: Parameter names to format as a multi-line. + """ columns = [self.aliased_values[k] for k in params] rows = [ " ".join(self._format_value(v) for v in row) @@ -239,7 +263,11 @@ def format_cfl_file(self) -> str: return "\n".join(lines) + "\n" def write_cfl_file(self, filepath: Path) -> None: - """Write the formatted .cfl input file to disk.""" + """ + Write the formatted .cfl input file to disk. + + :param filepath: Path where the .cfl file will be written. + """ with open(filepath, mode="w", encoding="utf-8") as f: f.write(self.format_cfl_file()) @@ -257,7 +285,11 @@ def __init__(self, opts: SurveyOptions): self.opts = opts def _find_data_start(self, chunk: list[str]) -> int: - """Return the index of the first station data row within a section chunk.""" + """ + Return the index of the first station data row within a section chunk. + + :param chunk: Chunk of data containing a header to index lines from. + """ header_idx = next( i for i, line in enumerate(chunk) @@ -266,7 +298,12 @@ def _find_data_start(self, chunk: list[str]) -> int: return header_idx + 2 def _slice_data_lines(self, lines: list[str], anchor: str) -> list[str]: - """Slice the station data rows that follow the given section header.""" + """ + Slice the station data rows that follow the given section header. + + :param lines: Lines to slice. + :param anchor: String marking the start of a chunk of output data. + """ anchor_idx = next(i for i, line in enumerate(lines) if anchor in line) chunk = lines[anchor_idx:] data_start = self._find_data_start(chunk) @@ -275,13 +312,24 @@ def _slice_data_lines(self, lines: list[str], anchor: str) -> list[str]: def _extract_data( self, outfile: str | Path, component: Literal["inline", "crossline", "vertical"] ) -> np.ndarray: - """Extract channel data for a single component from a LeroiAir .out file.""" + """ + Extract channel data for a single component from a LeroiAir .out file. + + :param outfile: Path to the output file from a LeroiAir run. + :param component: Component to extract. + """ lines = Path(outfile).read_text(encoding="utf-8", errors="replace").splitlines() data_lines = self._slice_data_lines(lines, self._COMPONENT_ANCHORS[component]) return np.array([line.split() for line in data_lines], dtype=float)[:, 4:] def save_to_geoh5(self, outfile: str | Path, out_group): - """Save LeroiAir simulated data on the provided survey to geoh5.""" + """ + Save LeroiAir simulated data on the provided survey to geoh5. + + :param outfile: Path to output file from a LeroiAir run. + :param out_group: Group where a copy of the survey will be saved + along with all the data computed by LeroiAir. + """ survey = self.opts.entity.copy(parent=out_group, copy_children=False) for component in "inline", "crossline", "vertical": From bea6cf49255010047c5593c2bc7274286e858bdb Mon Sep 17 00:00:00 2001 From: benjamink Date: Wed, 29 Apr 2026 09:22:23 -0700 Subject: [PATCH 36/43] undo disable protected-access for all tests --- .pre-commit-config.yaml | 10 +--------- 1 file changed, 1 insertion(+), 9 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index e1016162..bdd3023b 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -65,15 +65,7 @@ repos: language: system require_serial: true # pylint does its own parallelism types: [python] - exclude: ^(devtools|docs|tests)/ - - id: pylint-tests - name: pylint (tests) - entry: .\\devtools\\conda_env_pylint.bat - language: system - require_serial: true - args: [--rcfile=pylintrc, --disable=protected-access] - types: [python] - files: ^tests/ + exclude: ^(devtools|docs)/ - repo: https://github.com/codespell-project/codespell rev: v2.4.2 hooks: From f194b2a83f37757577baecd0bca4632805ff17a8 Mon Sep 17 00:00:00 2001 From: benjamink Date: Wed, 29 Apr 2026 09:33:09 -0700 Subject: [PATCH 37/43] disable protected-access on pylint for interface_test --- tests/plate_simulation/leroi_air/interface_test.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/tests/plate_simulation/leroi_air/interface_test.py b/tests/plate_simulation/leroi_air/interface_test.py index 7ce933c6..676a0d8b 100644 --- a/tests/plate_simulation/leroi_air/interface_test.py +++ b/tests/plate_simulation/leroi_air/interface_test.py @@ -20,6 +20,9 @@ from .conftest import FAKE_OUT, N_CHANNELS, N_STATIONS +# pylint: disable=protected-access + + def test_line_formatting(tmp_path): with Workspace(tmp_path / "test.geoh5") as geoh5: opts = generate_plate_options(geoh5) From 3af27d6e7e11295859d87a4b63f5e2513cabcfe1 Mon Sep 17 00:00:00 2001 From: benjamink Date: Wed, 29 Apr 2026 09:34:19 -0700 Subject: [PATCH 38/43] also driver_test --- tests/plate_simulation/leroi_air/driver_test.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/tests/plate_simulation/leroi_air/driver_test.py b/tests/plate_simulation/leroi_air/driver_test.py index bb7efa94..106b4b92 100644 --- a/tests/plate_simulation/leroi_air/driver_test.py +++ b/tests/plate_simulation/leroi_air/driver_test.py @@ -19,6 +19,9 @@ from . import generate_plate_options +# pylint: disable=protected-access + + def test_leroi_air_driver_run_raises_on_subprocess_failure(tmp_path): with Workspace(tmp_path / "test.geoh5") as geoh5: opts = generate_plate_options(workspace=geoh5) From 33e39bcd7926667dcf8a6cde3a7fdb7b1ed3baeb Mon Sep 17 00:00:00 2001 From: benjamink Date: Wed, 29 Apr 2026 10:57:01 -0700 Subject: [PATCH 39/43] Convert any timing parameters (timing_mark, channeld, waveform[:, 0]) to milliseconds if survey.unit is not ms --- .../plate_simulation/leroi_air/interface.py | 9 ++- .../plate_simulation/leroi_air/options.py | 56 +++++++++++++++---- .../leroi_air/options_test.py | 2 +- 3 files changed, 49 insertions(+), 18 deletions(-) diff --git a/simpeg_drivers/plate_simulation/leroi_air/interface.py b/simpeg_drivers/plate_simulation/leroi_air/interface.py index db175e0c..c89505f6 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/interface.py +++ b/simpeg_drivers/plate_simulation/leroi_air/interface.py @@ -45,14 +45,13 @@ def aliased_values(self) -> dict[str, Any]: "NSX": len(self.opts.survey.ontime_waveform), "STEP": 0 if self.opts.magnetic_field == "dBdt" else 1, "UNITS": 1, - "NCHNL": len(self.opts.survey.entity.channels), + "NCHNL": len(self.opts.survey.channels), "KRXW": 2, - "REFTYM": self.opts.survey.entity.timing_mark, + "REFTYM": self.opts.survey.timing_mark, "OFFTIME": self.opts.survey.offtime, "TXON": self.opts.survey.ontime_waveform[:, 0], "TXAMP": self.opts.survey.ontime_waveform[:, 1], - "TMS": self.opts.survey.entity.timing_mark - + np.array(self.opts.survey.entity.channels), + "TMS": self.opts.survey.timing_mark + np.array(self.opts.survey.channels), "WIDTH": self.opts.survey.channel_widths, "TXCLN": 0.0, "CMP": 3, @@ -337,7 +336,7 @@ def save_to_geoh5(self, outfile: str | Path, out_group): entities = survey.add_data( { f"fwd {component} [{i}]": {"values": data[:, i]} - for i in range(len(self.opts.entity.channels)) + for i in range(len(self.opts.channels)) } ) diff --git a/simpeg_drivers/plate_simulation/leroi_air/options.py b/simpeg_drivers/plate_simulation/leroi_air/options.py index 33f6e3c5..c3695650 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/options.py +++ b/simpeg_drivers/plate_simulation/leroi_air/options.py @@ -20,14 +20,44 @@ from scipy.interpolate import LinearNDInterpolator from simpeg_drivers.components.topography import InversionTopography +from simpeg_drivers.electromagnetics.time_domain.options import TDEMForwardOptions from simpeg_drivers.plate_simulation.options import ModelOptions +TIME_UNIT_CONVERSION = { + "Seconds (s)": 1e3, + "Milliseconds (ms)": 1, + "Microseconds (us)": 1e-3, + "Nanoseconds (ns)": 1e-6, +} + + class SurveyOptions(BaseModel): model_config = ConfigDict(arbitrary_types_allowed=True) entity: AirborneTEMReceivers + @property + def waveform(self) -> np.ndarray: + """Waveform in expected units of milliseconds.""" + waveform = self.entity.waveform.copy() + waveform[:, 0] *= TIME_UNIT_CONVERSION[self.entity.unit] + return waveform + + @property + def timing_mark(self) -> float: + """Timing mark in expected units of milliseconds.""" + timing_mark = self.entity.timing_mark + timing_mark *= TIME_UNIT_CONVERSION[self.entity.unit] + return timing_mark + + @property + def channels(self) -> np.ndarray: + """Channels in expected units of milliseconds.""" + channels = self.entity.channels.copy() + channels *= TIME_UNIT_CONVERSION[self.entity.unit] + return channels + @property def frequency(self) -> float: """ @@ -39,8 +69,8 @@ def frequency(self) -> float: """ frequency = self.entity.metadata.get("frequency", None) if frequency is None: - half_cycle_time = self.entity.waveform[-1, 0] - self.entity.waveform[0, 0] - frequency = 1 / (2 * half_cycle_time) + half_cycle_seconds = (self.waveform[-1, 0] - self.waveform[0, 0]) / 1e3 + frequency = 1 / (2 * half_cycle_seconds) return frequency @@ -49,13 +79,15 @@ def channel_widths(self) -> np.ndarray: """Time channel widths.""" channel_widths = self.entity.metadata.get("channel_widths", None) if channel_widths is None: - channel_widths = np.diff(np.r_[0, self.entity.channels]) + channel_widths = np.diff(np.r_[0, self.channels]) + else: + channel_widths *= TIME_UNIT_CONVERSION[self.entity.unit] return channel_widths @property def _ontime(self) -> float: """Time at which the transmitter current turns off.""" - return float(self.entity.waveform[self._offtime_mask(), 0][0]) + return float(self.waveform[self._offtime_mask(), 0][0]) @property def offtime(self) -> float: @@ -65,18 +97,18 @@ def offtime(self) -> float: This offtime is based on system frequency and may include times that are not accounted for by the waveform. """ - half_cycle = 1 / (2 * self.frequency) - zero_current_ind = np.where(self.entity.waveform[:, 1] == 0)[0] - first_zero = 1 if self.entity.waveform[0, 1] == 0.0 else 0 - ontime = self.entity.waveform[zero_current_ind][first_zero, 0] + half_cycle = 1000 / (2 * self.frequency) + zero_current_ind = np.where(self.waveform[:, 1] == 0)[0] + first_zero = 1 if self.waveform[0, 1] == 0.0 else 0 + ontime = self.waveform[zero_current_ind][first_zero, 0] return half_cycle - ontime @property def ontime_waveform(self) -> np.ndarray: """On-time waveform including leading and trailing 0 current times.""" - ontime_waveform = self.entity.waveform[~self._offtime_mask(), :] - endpoint = self.entity.waveform[self._offtime_mask()][0, :] + ontime_waveform = self.waveform[~self._offtime_mask(), :] + endpoint = self.waveform[self._offtime_mask()][0, :] return np.vstack([ontime_waveform, endpoint]) @property @@ -101,7 +133,7 @@ def drape_height(self, topo: np.ndarray) -> np.ndarray: def _offtime_mask(self) -> np.ndarray: """Mask selecting off-time rows from the waveform array.""" - mask = np.isclose(self.entity.waveform[:, 1], 0) + mask = np.isclose(self.waveform[:, 1], 0) mask[0] = False return mask @@ -145,7 +177,7 @@ def validate_plate_lengths_match(self) -> Self: @classmethod def from_plate_simulation_options( - cls, model: ModelOptions, simulation: SimPEGGroup + cls, model: ModelOptions, simulation: TDEMForwardOptions ) -> Self: """Construct from a :class:`PlateSimulationOptions` instance.""" diff --git a/tests/plate_simulation/leroi_air/options_test.py b/tests/plate_simulation/leroi_air/options_test.py index c400fced..93c3e380 100644 --- a/tests/plate_simulation/leroi_air/options_test.py +++ b/tests/plate_simulation/leroi_air/options_test.py @@ -33,7 +33,7 @@ def test_options_channel_widths_from_metadata(tmp_path): def test_options_frequency(plate_options): - assert np.isclose(plate_options.survey.frequency, 0.142857143) + assert np.isclose(plate_options.survey.frequency, 142.857143) def test_options_frequency_from_metadata(tmp_path): From 960544209b42dfb01f5dc1726cef3c21adb8c784 Mon Sep 17 00:00:00 2001 From: benjamink Date: Wed, 29 Apr 2026 16:43:13 -0700 Subject: [PATCH 40/43] add normalization to convert nT -> T, and simplify dbdt/b field selection --- simpeg_drivers/plate_simulation/leroi_air/driver.py | 1 + .../plate_simulation/leroi_air/interface.py | 11 ++++++++--- simpeg_drivers/plate_simulation/leroi_air/options.py | 4 ++-- 3 files changed, 11 insertions(+), 5 deletions(-) diff --git a/simpeg_drivers/plate_simulation/leroi_air/driver.py b/simpeg_drivers/plate_simulation/leroi_air/driver.py index 1539e7a7..143dd2ef 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/driver.py +++ b/simpeg_drivers/plate_simulation/leroi_air/driver.py @@ -55,4 +55,5 @@ def run(self) -> None: self.interface.output.save_to_geoh5( outfile=self.project_path / "LeroiAir.out", out_group=self.options.out_group, + normalization=1e-9, ) diff --git a/simpeg_drivers/plate_simulation/leroi_air/interface.py b/simpeg_drivers/plate_simulation/leroi_air/interface.py index c89505f6..4357164f 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/interface.py +++ b/simpeg_drivers/plate_simulation/leroi_air/interface.py @@ -14,6 +14,7 @@ from typing import Any, Literal import numpy as np +from geoh5py.groups import UIJsonGroup from geoh5py.shared.utils import fetch_active_workspace from .options import LeroiAirOptions, SurveyOptions @@ -43,7 +44,7 @@ def aliased_values(self) -> dict[str, Any]: "ISTOP": 0, "ISW": 1, "NSX": len(self.opts.survey.ontime_waveform), - "STEP": 0 if self.opts.magnetic_field == "dBdt" else 1, + "STEP": int(self.opts.step), "UNITS": 1, "NCHNL": len(self.opts.survey.channels), "KRXW": 2, @@ -321,13 +322,17 @@ def _extract_data( data_lines = self._slice_data_lines(lines, self._COMPONENT_ANCHORS[component]) return np.array([line.split() for line in data_lines], dtype=float)[:, 4:] - def save_to_geoh5(self, outfile: str | Path, out_group): + def save_to_geoh5( + self, outfile: str | Path, out_group: UIJsonGroup, normalization: float = 1 + ): """ Save LeroiAir simulated data on the provided survey to geoh5. :param outfile: Path to output file from a LeroiAir run. :param out_group: Group where a copy of the survey will be saved along with all the data computed by LeroiAir. + :param normalization: Normalization multiplied against the data + before saving. """ survey = self.opts.entity.copy(parent=out_group, copy_children=False) @@ -335,7 +340,7 @@ def save_to_geoh5(self, outfile: str | Path, out_group): data = self._extract_data(outfile=outfile, component=component) entities = survey.add_data( { - f"fwd {component} [{i}]": {"values": data[:, i]} + f"fwd {component} [{i}]": {"values": data[:, i] * normalization} for i in range(len(self.opts.channels)) } ) diff --git a/simpeg_drivers/plate_simulation/leroi_air/options.py b/simpeg_drivers/plate_simulation/leroi_air/options.py index c3695650..5ae02c3d 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/options.py +++ b/simpeg_drivers/plate_simulation/leroi_air/options.py @@ -151,7 +151,7 @@ class LeroiAirOptions(BaseModel): plate_resistivities: list[float] plate_geometries: list[PlateModel] cell_size: float = 10.0 - magnetic_field: Literal["dBdt", "B"] = "dBdt" + step: bool = True domain: Literal["time", "frequency"] = "time" layered_earth_only: bool = False float_precision: int = 4 @@ -202,7 +202,7 @@ def from_plate_simulation_options( layer_thicknesses=[model.overburden_options.thickness, 9999], plate_resistivities=[model.plate_options.plate_property], plate_geometries=[model.plate_options.geometry], - magnetic_field="dBdt" if "dBdt" in simulation.data_units else "B", + step="dBdt" in simulation.data_units, out_group=simulation.out_group, ) From 8b1818dc4b3c74af630ccd3c3b2aa44d8ce02648 Mon Sep 17 00:00:00 2001 From: benjamink Date: Wed, 29 Apr 2026 16:49:45 -0700 Subject: [PATCH 41/43] remove todo for time unit conversion --- simpeg_drivers/plate_simulation/leroi_air/options.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/simpeg_drivers/plate_simulation/leroi_air/options.py b/simpeg_drivers/plate_simulation/leroi_air/options.py index 5ae02c3d..c2f78f68 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/options.py +++ b/simpeg_drivers/plate_simulation/leroi_air/options.py @@ -216,8 +216,6 @@ def n_plates(self) -> int: """Number of plates.""" return len(self.plate_geometries) - # TODO use units to convert time to milliseconds used by leroi. - @property def resistivities(self) -> np.ndarray: """All resistivities.""" From 5191acb85beefda05d812ff32b07013cc30f09be Mon Sep 17 00:00:00 2001 From: benjamink Date: Wed, 29 Apr 2026 17:00:14 -0700 Subject: [PATCH 42/43] fix test --- simpeg_drivers/plate_simulation/leroi_air/options.py | 2 +- tests/plate_simulation/leroi_air/interface_test.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/simpeg_drivers/plate_simulation/leroi_air/options.py b/simpeg_drivers/plate_simulation/leroi_air/options.py index c2f78f68..834018ef 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/options.py +++ b/simpeg_drivers/plate_simulation/leroi_air/options.py @@ -202,7 +202,7 @@ def from_plate_simulation_options( layer_thicknesses=[model.overburden_options.thickness, 9999], plate_resistivities=[model.plate_options.plate_property], plate_geometries=[model.plate_options.geometry], - step="dBdt" in simulation.data_units, + step="dBdt" in simulation.get("data_units", "dBdt"), out_group=simulation.out_group, ) diff --git a/tests/plate_simulation/leroi_air/interface_test.py b/tests/plate_simulation/leroi_air/interface_test.py index 676a0d8b..124f81b5 100644 --- a/tests/plate_simulation/leroi_air/interface_test.py +++ b/tests/plate_simulation/leroi_air/interface_test.py @@ -53,7 +53,7 @@ def test_format_cfl_file(tmp_path): assert "1 1 1 0 \t ! TDFD, DO3D, PRFL, ISTOP" in cfl_str # record_3: ISW, NSX, STEP, UNITS, NCHNL, KRXW, OFFTIME - assert "1 8 0 1 3 2 1.6 \t ! ISW, NSX, STEP, UNITS, NCHNL, KRXW, OFFTIME" in cfl_str + assert "1 8 1 1 3 2 1.6 \t ! ISW, NSX, STEP, UNITS, NCHNL, KRXW, OFFTIME" in cfl_str # record_4: TXON, TXAMP (multi-line waveform table) assert ( From 33bcc43b6b4c16a93002fa1cc86847e6e6d784e4 Mon Sep 17 00:00:00 2001 From: benjamink Date: Thu, 30 Apr 2026 11:47:50 -0700 Subject: [PATCH 43/43] out_group is not optional for options class --- simpeg_drivers/plate_simulation/leroi_air/options.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/simpeg_drivers/plate_simulation/leroi_air/options.py b/simpeg_drivers/plate_simulation/leroi_air/options.py index 834018ef..b8486bc8 100644 --- a/simpeg_drivers/plate_simulation/leroi_air/options.py +++ b/simpeg_drivers/plate_simulation/leroi_air/options.py @@ -144,6 +144,7 @@ class LeroiAirOptions(BaseModel): model_config = ConfigDict(arbitrary_types_allowed=True) title: str = "LeroiAir modelling for plate-simulation package." + out_group: SimPEGGroup survey: SurveyOptions topo: np.ndarray layer_resistivities: list[float] @@ -155,7 +156,6 @@ class LeroiAirOptions(BaseModel): domain: Literal["time", "frequency"] = "time" layered_earth_only: bool = False float_precision: int = 4 - out_group: SimPEGGroup | None = None @model_validator(mode="after") def validate_layer_lengths_match(self) -> Self: