From 843b5c7813a0cdddcecbb0940450adb67af98771 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Thu, 11 Sep 2025 11:27:16 -0700 Subject: [PATCH 1/4] Add check for is_sigma positive. Add unitest --- simpeg_drivers/components/models.py | 32 +++++++++++++++++++++-------- 1 file changed, 24 insertions(+), 8 deletions(-) diff --git a/simpeg_drivers/components/models.py b/simpeg_drivers/components/models.py index eb3530202..7b9e8d98f 100644 --- a/simpeg_drivers/components/models.py +++ b/simpeg_drivers/components/models.py @@ -14,14 +14,12 @@ from typing import TYPE_CHECKING import numpy as np -from geoapps_utils.base import Driver +from geoapps_utils.utils.importing import GeoAppsError from geoapps_utils.utils.numerical import weighted_average -from geoapps_utils.utils.transformations import rotate_xyz from geoh5py.data import Data, FloatData, NumericData from geoh5py.data.data_type import GeometricDataValueMapType from geoh5py.objects import ObjectBase from simpeg.utils.mat_utils import ( - cartesian2amplitude_dip_azimuth, dip_azimuth2cartesian, mkvc, ) @@ -76,15 +74,25 @@ def __init__(self, driver: InversionDriver): self.is_sigma = self.driver.params.physical_property == "conductivity" self.is_vector = self.driver.params.inversion_type == "magnetic vector" - self._starting_model = InversionModel(driver, "starting_model") + self._starting_model = InversionModel( + driver, "starting_model", is_sigma=self.is_sigma + ) self._starting_inclination = InversionModel(driver, "starting_inclination") self._starting_declination = InversionModel(driver, "starting_declination") - self._reference_model = InversionModel(driver, "reference_model") + self._reference_model = InversionModel( + driver, "reference_model", is_sigma=self.is_sigma + ) self._reference_inclination = InversionModel(driver, "reference_inclination") self._reference_declination = InversionModel(driver, "reference_declination") - self._lower_bound = InversionModel(driver, "lower_bound") - self._upper_bound = InversionModel(driver, "upper_bound") - self._conductivity_model = InversionModel(driver, "conductivity_model") + self._lower_bound = InversionModel( + driver, "lower_bound", is_sigma=self.is_sigma + ) + self._upper_bound = InversionModel( + driver, "upper_bound", is_sigma=self.is_sigma + ) + self._conductivity_model = InversionModel( + driver, "conductivity_model", is_sigma=self.is_sigma + ) self._alpha_s = InversionModel(driver, "alpha_s") self._length_scale_x = InversionModel(driver, "length_scale_x") self._length_scale_y = InversionModel(driver, "length_scale_y") @@ -498,6 +506,7 @@ def __init__( driver: InversionDriver, model_type: str, trim_active_cells: bool = True, + is_sigma: bool = False, ): """ :param driver: InversionDriver object. @@ -507,6 +516,7 @@ def __init__( """ self.driver = driver self.model_type = model_type + self.is_sigma = is_sigma self.model: np.ndarray | None = None self.trim_active_cells = trim_active_cells self._initialize() @@ -522,6 +532,12 @@ def _initialize(self): model = self._get(self.model_type) if model is not None: + if self.is_sigma and np.any(model <= 0): + raise GeoAppsError( + f"All values in {self.model_type} must be positive when " + "inversion is in log-conductivity space." + ) + self.model = mkvc(model) if isinstance(self._fetch_reference(self.model_type), Data): From ab07b315a5545fa25e6fd0c7ba96592c5cacf14e Mon Sep 17 00:00:00 2001 From: dominiquef Date: Thu, 11 Sep 2025 11:33:47 -0700 Subject: [PATCH 2/4] Use GeoAppsError where needed --- .../components/factories/entity_factory.py | 3 ++- .../components/factories/survey_factory.py | 6 ++--- simpeg_drivers/driver.py | 22 +++++++++++++------ .../frequency_domain/options.py | 3 ++- simpeg_drivers/joint/joint_surveys/options.py | 3 ++- simpeg_drivers/line_sweep/driver.py | 3 ++- simpeg_drivers/options.py | 12 +++++----- 7 files changed, 32 insertions(+), 20 deletions(-) diff --git a/simpeg_drivers/components/factories/entity_factory.py b/simpeg_drivers/components/factories/entity_factory.py index 1222272a9..10b250b70 100644 --- a/simpeg_drivers/components/factories/entity_factory.py +++ b/simpeg_drivers/components/factories/entity_factory.py @@ -17,6 +17,7 @@ from typing import TYPE_CHECKING import numpy as np +from geoapps_utils.utils.importing import GeoAppsError from geoh5py.objects import ( CurrentElectrode, Curve, @@ -128,7 +129,7 @@ def _validate_large_loop_cells( Validate that the transmitter loops are counter-clockwise sorted and closed. """ if transmitter.receivers.tx_id_property is None: - raise ValueError( + raise GeoAppsError( "Transmitter ID property required for LargeLoopGroundTEMReceivers" ) diff --git a/simpeg_drivers/components/factories/survey_factory.py b/simpeg_drivers/components/factories/survey_factory.py index b98c47cd9..2f7c2c2e3 100644 --- a/simpeg_drivers/components/factories/survey_factory.py +++ b/simpeg_drivers/components/factories/survey_factory.py @@ -24,7 +24,7 @@ import numpy as np import simpeg.electromagnetics.time_domain as tdem -from geoh5py.objects.surveys.electromagnetics.airborne_fem import AirborneFEMReceivers +from geoapps_utils.utils.importing import GeoAppsError from geoh5py.objects.surveys.electromagnetics.ground_tem import ( LargeLoopGroundTEMTransmitters, ) @@ -207,14 +207,14 @@ def _tdem_arguments(self, data=None): if receivers.channels[-1] > ( receivers.waveform[:, 0].max() - receivers.timing_mark ): - raise ValueError( + raise GeoAppsError( f"The latest time channel {receivers.channels[-1]} exceeds " f"the waveform discretization. Revise waveform." ) if isinstance(transmitters, LargeLoopGroundTEMTransmitters): if receivers.tx_id_property is None: - raise ValueError( + raise GeoAppsError( "Transmitter ID property required for LargeLoopGroundTEMReceivers" ) diff --git a/simpeg_drivers/driver.py b/simpeg_drivers/driver.py index 41f450876..35f5bbb47 100644 --- a/simpeg_drivers/driver.py +++ b/simpeg_drivers/driver.py @@ -84,6 +84,9 @@ mlogger.setLevel(logging.WARNING) +logger = logging.getLogger("simpeg-drivers") + + class InversionDriver(Driver): _options_class = BaseForwardOptions | BaseInversionOptions _inversion_type: str | None = None @@ -646,14 +649,19 @@ def start( else: ifile = InputFile.read_ui_json(filepath, **kwargs) - if driver_class is None: - driver = cls.from_input_file(ifile) - else: - with ifile.data["geoh5"].open(mode="r+"): - params = driver_class._options_class.build(ifile) - driver = driver_class(params) + try: + if driver_class is None: + driver = cls.from_input_file(ifile) + else: + with ifile.data["geoh5"].open(mode="r+"): + params = driver_class._options_class.build(ifile) + driver = driver_class(params) + + driver.run() - driver.run() + except GeoAppsError as error: + logger.warning("\n\nApplicationError: %s\n\n", error) + sys.exit(1) return driver diff --git a/simpeg_drivers/electromagnetics/frequency_domain/options.py b/simpeg_drivers/electromagnetics/frequency_domain/options.py index ec7928b73..c92dc43a8 100644 --- a/simpeg_drivers/electromagnetics/frequency_domain/options.py +++ b/simpeg_drivers/electromagnetics/frequency_domain/options.py @@ -15,6 +15,7 @@ from pathlib import Path from typing import ClassVar, TypeAlias +from geoapps_utils.utils.importing import GeoAppsError from geoh5py.groups import PropertyGroup from geoh5py.objects import ( AirborneFEMReceivers, @@ -56,7 +57,7 @@ def tx_offsets(self): except KeyError as exception: msg = "Metadata must contain 'Frequency configurations' dictionary with 'Offset' data." - raise KeyError(msg) from exception + raise GeoAppsError(msg) from exception return tx_offsets diff --git a/simpeg_drivers/joint/joint_surveys/options.py b/simpeg_drivers/joint/joint_surveys/options.py index 296bdf77d..ac06aa3b7 100644 --- a/simpeg_drivers/joint/joint_surveys/options.py +++ b/simpeg_drivers/joint/joint_surveys/options.py @@ -14,6 +14,7 @@ from pathlib import Path from typing import ClassVar +from geoapps_utils.utils.importing import GeoAppsError from pydantic import model_validator from simpeg_drivers import assets_path @@ -38,7 +39,7 @@ class JointSurveysOptions(BaseJointOptions): def all_groups_same_physical_property(self): physical_properties = [k.options["physical_property"] for k in self.groups] if len(list(set(physical_properties))) > 1: - raise ValueError( + raise GeoAppsError( "All physical properties must be the same. " f"Provided SimPEG groups for {physical_properties}." ) diff --git a/simpeg_drivers/line_sweep/driver.py b/simpeg_drivers/line_sweep/driver.py index 20b25837f..c21435609 100644 --- a/simpeg_drivers/line_sweep/driver.py +++ b/simpeg_drivers/line_sweep/driver.py @@ -18,6 +18,7 @@ import numpy as np from geoapps_utils.param_sweeps.driver import SweepDriver, SweepParams from geoapps_utils.param_sweeps.generate import generate +from geoapps_utils.utils.importing import GeoAppsError from geoh5py.data import FilenameData from geoh5py.groups import SimPEGGroup from geoh5py.objects import DrapeModel, PotentialElectrode @@ -147,7 +148,7 @@ def collect_results(self): ) if not line_data: - raise ValueError(f"Line {line} not found in {survey.name}") + raise GeoAppsError(f"Line {line} not found in {survey.name}") line_indices = line_ids == line data = self.collect_line_data(survey, line_indices, data) diff --git a/simpeg_drivers/options.py b/simpeg_drivers/options.py index 057221536..e7e656e66 100644 --- a/simpeg_drivers/options.py +++ b/simpeg_drivers/options.py @@ -17,6 +17,7 @@ import numpy as np from geoapps_utils.base import Options +from geoapps_utils.utils.importing import GeoAppsError from geoh5py.data import ( BooleanData, DataAssociationEnum, @@ -28,9 +29,8 @@ from geoh5py.groups import PropertyGroup, SimPEGGroup, UIJsonGroup from geoh5py.objects import DrapeModel, Grid2D, Octree, Points from geoh5py.objects.surveys.electromagnetics.base import BaseEMSurvey -from geoh5py.shared.utils import fetch_active_workspace from geoh5py.ui_json import InputFile -from geoh5py.ui_json.templates import data_parameter +from geoh5py.ui_json.utils import fetch_active_workspace from pydantic import ( AliasChoices, BaseModel, @@ -89,7 +89,7 @@ class ActiveCellsOptions(BaseModel): @classmethod def at_least_one(cls, data): if all(v is None for v in data.values()): - raise ValueError("Must provide either topography or active model.") + raise GeoAppsError("Must provide either topography or active model.") return data @model_serializer(mode="wrap") @@ -208,7 +208,7 @@ def _component_name(self, component: str) -> str: @classmethod def mesh_cannot_be_rotated(cls, value: Octree): if isinstance(value, Octree) and value.rotation not in [0.0, None]: - raise ValueError( + raise GeoAppsError( "Rotated meshes are not supported. Please use a mesh with an angle of 0.0." ) return value @@ -520,13 +520,13 @@ class LineSelectionOptions(BaseModel): @classmethod def validate_cell_association(cls, value): if value.association is not DataAssociationEnum.CELL: - raise ValueError("Line identifier must be associated with cells.") + raise GeoAppsError("Line identifier must be associated with cells.") return value @model_validator(mode="after") def line_id_referenced(self): if self.line_id not in self.line_object.values: - raise ValueError("Line id isn't referenced in the line object.") + raise GeoAppsError("Line id isn't referenced in the line object.") return self From c6f200c714118cc283bd50b5cb96a911960cd9f2 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Thu, 11 Sep 2025 11:35:25 -0700 Subject: [PATCH 3/4] Add unitest --- tests/models_test.py | 42 ++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 40 insertions(+), 2 deletions(-) diff --git a/tests/models_test.py b/tests/models_test.py index c2f5421b5..acc6d4b7c 100644 --- a/tests/models_test.py +++ b/tests/models_test.py @@ -13,7 +13,7 @@ from pathlib import Path import numpy as np -from geoh5py import Workspace +import pytest from geoh5py.objects import Points from simpeg_drivers.components import ( @@ -21,6 +21,12 @@ InversionModel, InversionModelCollection, ) +from simpeg_drivers.electricals.direct_current.three_dimensions.driver import ( + DC3DForwardDriver, +) +from simpeg_drivers.electricals.direct_current.three_dimensions.options import ( + DC3DForwardOptions, +) from simpeg_drivers.options import ActiveCellsOptions from simpeg_drivers.potential_fields import MVIInversionOptions from simpeg_drivers.potential_fields.magnetic_vector.driver import ( @@ -43,7 +49,7 @@ def get_mvi_params(tmp_path: Path) -> MVIInversionOptions: mesh=MeshOptions(refinement=(2,)), model=ModelOptions(anomaly=0.05), ) - with get_workspace(tmp_path / "inversion_test.ui.geoh5") as geoh5: + with get_workspace(tmp_path / f"{__name__}.ui.geoh5") as geoh5: components = SyntheticsComponents(geoh5, options=opts) mesh = components.model.parent ref_inducing = mesh.add_data( @@ -81,6 +87,38 @@ def get_mvi_params(tmp_path: Path) -> MVIInversionOptions: return params +def get_dc_params(tmp_path: Path) -> MVIInversionOptions: + opts = SyntheticsComponentsOptions( + method="direct_current", + survey=SurveyOptions(n_stations=4, n_lines=2), + mesh=MeshOptions(refinement=(2,)), + model=ModelOptions(anomaly=0.05), + ) + with get_workspace(tmp_path / f"{__name__}.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + mesh = components.model.parent + params = DC3DForwardOptions.build( + geoh5=geoh5, + data_object=components.survey, + tmi_channel_bool=True, + mesh=mesh, + active_cells=ActiveCellsOptions(topography_object=components.topography), + starting_model=-1e-04, + ) + + return params + + +def test_negative_reference_model(tmp_path: Path): + params = get_dc_params(tmp_path) + geoh5 = params.geoh5 + with geoh5.open(): + driver = DC3DForwardDriver(params) + + with pytest.raises(ValueError, match="must be positive when"): + _ = driver.models.starting_model + + def test_zero_reference_model(tmp_path: Path): params = get_mvi_params(tmp_path) geoh5 = params.geoh5 From 22b6ded4023b27d36aba3cc16d3cf6bda55a9576 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Thu, 11 Sep 2025 14:44:53 -0700 Subject: [PATCH 4/4] Fix tests --- tests/models_test.py | 3 ++- tests/run_tests/driver_airborne_tem_test.py | 4 ++-- tests/run_tests/driver_joint_surveys_test.py | 2 +- 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/tests/models_test.py b/tests/models_test.py index acc6d4b7c..e080469b0 100644 --- a/tests/models_test.py +++ b/tests/models_test.py @@ -14,6 +14,7 @@ import numpy as np import pytest +from geoapps_utils.utils.importing import GeoAppsError from geoh5py.objects import Points from simpeg_drivers.components import ( @@ -115,7 +116,7 @@ def test_negative_reference_model(tmp_path: Path): with geoh5.open(): driver = DC3DForwardDriver(params) - with pytest.raises(ValueError, match="must be positive when"): + with pytest.raises(GeoAppsError, match="must be positive when"): _ = driver.models.starting_model diff --git a/tests/run_tests/driver_airborne_tem_test.py b/tests/run_tests/driver_airborne_tem_test.py index 86114119c..fe5b426b7 100644 --- a/tests/run_tests/driver_airborne_tem_test.py +++ b/tests/run_tests/driver_airborne_tem_test.py @@ -13,7 +13,7 @@ from pathlib import Path import numpy as np -from geoh5py.groups import SimPEGGroup +from geoapps_utils.utils.importing import GeoAppsError from geoh5py.workspace import Workspace from pytest import raises @@ -71,7 +71,7 @@ def test_bad_waveform(tmp_path: Path): params.data_object.channels[-1] = 1000.0 - with raises(ValueError, match="The latest time"): + with raises(GeoAppsError, match="The latest time"): _ = fwr_driver.inversion_data.survey diff --git a/tests/run_tests/driver_joint_surveys_test.py b/tests/run_tests/driver_joint_surveys_test.py index 3fe5c2804..052879063 100644 --- a/tests/run_tests/driver_joint_surveys_test.py +++ b/tests/run_tests/driver_joint_surveys_test.py @@ -199,7 +199,7 @@ def test_joint_surveys_conductivity_run( method="direct-current", survey=SurveyOptions(n_stations=4, n_lines=4, name="survey A"), mesh=MeshOptions(refinement=(2, 2, 2), name="mesh A"), - model=ModelOptions(anomaly=0.1, name="model A"), + model=ModelOptions(anomaly=0.1, background=0.01, name="model A"), active=SyntheticsActiveCellsOptions(name="active A"), )