From 2f4fa485b207b7b8018728a698cb904a75517b94 Mon Sep 17 00:00:00 2001 From: benjamink Date: Tue, 2 Dec 2025 10:14:39 -0800 Subject: [PATCH 1/5] Update size of tilemap for mvi in joint inversion --- simpeg_drivers/joint/joint_surveys/driver.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/simpeg_drivers/joint/joint_surveys/driver.py b/simpeg_drivers/joint/joint_surveys/driver.py index c0906560..f807c4d2 100644 --- a/simpeg_drivers/joint/joint_surveys/driver.py +++ b/simpeg_drivers/joint/joint_surveys/driver.py @@ -74,7 +74,10 @@ def validate_create_models(self): def wires(self): """Model projections""" if self._wires is None: - wires = [maps.IdentityMap(nP=self.models.n_active) for _ in self.drivers] + wires = [ + maps.IdentityMap(nP=driver.n_values * driver.n_blocks) + for driver in self.drivers + ] self._wires = wires return self._wires From 253f2fd2f27e98689502c105fc2c7de78995890b Mon Sep 17 00:00:00 2001 From: benjamink Date: Tue, 2 Dec 2025 10:33:28 -0800 Subject: [PATCH 2/5] use n_active in identitymap, not n_values --- simpeg_drivers/joint/joint_surveys/driver.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/simpeg_drivers/joint/joint_surveys/driver.py b/simpeg_drivers/joint/joint_surveys/driver.py index f807c4d2..fdb52abf 100644 --- a/simpeg_drivers/joint/joint_surveys/driver.py +++ b/simpeg_drivers/joint/joint_surveys/driver.py @@ -75,7 +75,7 @@ def wires(self): """Model projections""" if self._wires is None: wires = [ - maps.IdentityMap(nP=driver.n_values * driver.n_blocks) + maps.IdentityMap(nP=self.models.n_active * driver.n_blocks) for driver in self.drivers ] self._wires = wires From aa26839653dacfab7a615cf24ed11b99b6d96555 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Wed, 7 Jan 2026 10:03:57 -0500 Subject: [PATCH 3/5] Change mechanics for mapping models in joint surveys. Add unit test --- simpeg_drivers/joint/driver.py | 3 + simpeg_drivers/joint/joint_surveys/driver.py | 23 ++++-- simpeg_drivers/utils/synthetics/driver.py | 4 +- tests/run_tests/driver_joint_surveys_test.py | 87 +++++++++++++++++++- 4 files changed, 107 insertions(+), 10 deletions(-) diff --git a/simpeg_drivers/joint/driver.py b/simpeg_drivers/joint/driver.py index bc93f59e..c3f9ac7e 100644 --- a/simpeg_drivers/joint/driver.py +++ b/simpeg_drivers/joint/driver.py @@ -129,6 +129,7 @@ def initialize(self): self.models.active_cells = global_actives for driver, wire in zip(self.drivers, self.wires, strict=True): logger.info("Initializing driver %s", driver.params.name) + # Create a projection from global mesh to driver specific mesh projection = TileMap( self.inversion_mesh.mesh, global_actives, @@ -139,6 +140,8 @@ def initialize(self): tile_map = projection * wire driver.params.active_model = None driver.models.active_cells = projection.local_active + + # Keep a copy on the top combo/future for saving directives and model creation driver.data_misfit.model_map = tile_map multipliers = [] diff --git a/simpeg_drivers/joint/joint_surveys/driver.py b/simpeg_drivers/joint/joint_surveys/driver.py index fdb52abf..02b45979 100644 --- a/simpeg_drivers/joint/joint_surveys/driver.py +++ b/simpeg_drivers/joint/joint_surveys/driver.py @@ -40,19 +40,25 @@ def __init__(self, params: JointSurveysOptions): def validate_create_models(self): """Check if all models were provided, otherwise use the first driver models.""" + # Create projection for first driver to global mesh + mapping = maps.TileMap( + self.inversion_mesh.mesh, + self.models.active_cells, + self.drivers[0].inversion_mesh.mesh, + enforce_active=False, + ) + projection = mapping.deriv(np.ones(self.models.n_active)).T + norm = np.array(np.sum(projection, axis=1)).flatten() + for model_type in self.models.model_types: model = getattr(self.models, model_type) if model is not None or getattr(self.drivers[0].models, model_type) is None: continue model_local_values = getattr(self.drivers[0].models, model_type) - projection = ( - self.drivers[0] - .data_misfit.model_map.deriv(np.ones(self.models.n_active)) - .T - ) - norm = np.array(np.sum(projection, axis=1)).flatten() - model = (projection * model_local_values) / (norm + 1e-8) + model = ( + projection * model_local_values[: self.drivers[0].models.n_active] + ) / (norm + 1e-8) if self.drivers[0].models.is_sigma and model_type in [ "starting_model", @@ -70,6 +76,9 @@ def validate_create_models(self): getattr(self.models, f"_{model_type}").model = model + # For MVI, set is_vector from first driver + self.models.is_vector = self.drivers[0].models.is_vector + @property def wires(self): """Model projections""" diff --git a/simpeg_drivers/utils/synthetics/driver.py b/simpeg_drivers/utils/synthetics/driver.py index abd7f714..7dc4c889 100644 --- a/simpeg_drivers/utils/synthetics/driver.py +++ b/simpeg_drivers/utils/synthetics/driver.py @@ -86,7 +86,7 @@ def mesh(self): @property def active(self): - entity = self.geoh5.get_entity(self.options.active.name)[0] + entity = self.mesh.get_entity(self.options.active.name)[0] assert isinstance(entity, FloatData | type(None)) if entity is None: entity = get_active(self.mesh, self.topography) @@ -95,7 +95,7 @@ def active(self): @property def model(self): - entity = self.geoh5.get_entity(self.options.model.name)[0] + entity = self.mesh.get_entity(self.options.model.name)[0] assert isinstance(entity, FloatData | type(None)) if entity is None: assert self.options is not None diff --git a/tests/run_tests/driver_joint_surveys_test.py b/tests/run_tests/driver_joint_surveys_test.py index 29537085..6facb2b6 100644 --- a/tests/run_tests/driver_joint_surveys_test.py +++ b/tests/run_tests/driver_joint_surveys_test.py @@ -13,7 +13,7 @@ import numpy as np from geoh5py.objects import Octree from geoh5py.workspace import Workspace -from simpeg.directives import SavePropertyGroup +from simpeg.directives import SaveModelGeoH5, SavePropertyGroup from simpeg_drivers.electricals.direct_current.three_dimensions.driver import ( DC3DInversionDriver, @@ -31,6 +31,8 @@ GravityInversionOptions, ) from simpeg_drivers.potential_fields.gravity.driver import GravityInversionDriver +from simpeg_drivers.potential_fields.magnetic_vector.driver import MVIInversionDriver +from simpeg_drivers.potential_fields.magnetic_vector.options import MVIInversionOptions from simpeg_drivers.utils.synthetics.driver import ( SyntheticsComponents, ) @@ -195,6 +197,89 @@ def test_joint_surveys_inv_run( check_target(output, target_run) +def test_joint_surveys_mvi_run(tmp_path, anomaly=0.05): + drivers = [] + + with Workspace.create(tmp_path / f"{__name__}.geoh5") as geoh5: + for ii in range(1, 3): + opts = SyntheticsComponentsOptions( + method="magnetic_vector", + survey=SurveyOptions( + n_stations=3**ii, + n_lines=3**ii, + drape=5.0, + name=f"Survey Driver[{ii}]", + ), + mesh=MeshOptions(refinement=(2**ii, 2, 2), name=f"Mesh Driver[{ii}]"), + model=ModelOptions(anomaly=anomaly), + ) + components = SyntheticsComponents(geoh5, options=opts) + survey = components.survey + obs, uncrt = survey.add_data( + { + "TMI": {"values": np.random.randn(survey.n_vertices)}, + "Uncertainty": {"values": np.ones(survey.n_vertices) * 1e-3}, + } + ) + + # Add an inclination model on the first driver only to test handling of + # models from the main driver + if ii == 1: + model = components.model.values + model[model > 0] = 45.0 + model[model <= 0] = 90.0 + inc_mod = components.mesh.add_data( + {"Inclination Model": {"values": model}} + ) + else: + inc_mod = None + + params = MVIInversionOptions.build( + geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, + tmi_channel=obs, + tmi_uncertainty=uncrt, + inducing_field_strength=45000, + inducing_field_inclination=90.0, + inducing_field_declination=0.0, + data_object=survey, + starting_model=components.model, + starting_inclination=inc_mod, + reference_model=0.0, + ) + drivers.append(MVIInversionDriver(params)) + + # Run the inverse + joint_params = JointSurveysOptions.build( + geoh5=geoh5, + active_cells=ActiveCellsOptions(topography_object=components.topography), + group_a=drivers[0].out_group, + group_b=drivers[1].out_group, + starting_model=0.01, + # Default to Conductivity (S/m) + ) + + driver = JointSurveyDriver(joint_params) + assert np.isclose(driver.models.reference_model[0], 0) # Took it from driver_A + assert driver.models.starting_model.shape == (driver.models.n_active * 3,) + assert np.isclose( + driver.models.starting_model.max(), 0.01 * np.cos(np.deg2rad(45.0)) + ) + + # Test saving the starting models on each mesh (open file to validate) + assert ( + len( + [ + directive.write(0, driver.models.starting_model) + for directive in driver.directives.directive_list + if isinstance(directive, SaveModelGeoH5) + ] + ) + == 3 + ) + + def test_joint_surveys_conductivity_run( tmp_path, ): From 2c7d5e15c1d259c19b7d3e82934a0b12d8918cba Mon Sep 17 00:00:00 2001 From: dominiquef Date: Wed, 7 Jan 2026 11:22:35 -0500 Subject: [PATCH 4/5] Remove checker for mvi group --- simpeg_drivers/joint/joint_surveys/options.py | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/simpeg_drivers/joint/joint_surveys/options.py b/simpeg_drivers/joint/joint_surveys/options.py index e0bf5a09..edc0bbe9 100644 --- a/simpeg_drivers/joint/joint_surveys/options.py +++ b/simpeg_drivers/joint/joint_surveys/options.py @@ -49,19 +49,6 @@ class JointSurveysOptions(BaseJointOptions): models: JointSurveysModelOptions - @field_validator("group_a", "group_b", "group_c") - @classmethod - def no_mvi_groups(cls, val): - if val is None: - return val - - if "magnetic vector" in val.options.get("inversion_type", ""): - raise ValueError( - f"Joint inversion doesn't currently support MVI data as passed in " - f"the group: {val.name}." - ) - return val - @model_validator(mode="after") def all_groups_same_physical_property(self): physical_properties = [k.options["physical_property"] for k in self.groups] From 7b2a06a756c2168fa36a371b2fc528da3ca03280 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Wed, 7 Jan 2026 11:32:08 -0500 Subject: [PATCH 5/5] Fix typing --- simpeg_drivers/utils/synthetics/driver.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/simpeg_drivers/utils/synthetics/driver.py b/simpeg_drivers/utils/synthetics/driver.py index 9a148bff..c4f6adfd 100644 --- a/simpeg_drivers/utils/synthetics/driver.py +++ b/simpeg_drivers/utils/synthetics/driver.py @@ -10,7 +10,7 @@ import numpy as np from geoh5py import Workspace -from geoh5py.data import FloatData +from geoh5py.data import BooleanData, FloatData from geoh5py.objects import DrapeModel, ObjectBase, Octree, Surface from simpeg_drivers.utils.synthetics.meshes.factory import get_mesh @@ -87,7 +87,7 @@ def mesh(self): @property def active(self): entity = self.mesh.get_entity(self.options.active.name)[0] - assert isinstance(entity, FloatData | type(None)) + assert isinstance(entity, BooleanData | type(None)) if entity is None: entity = get_active(self.mesh, self.topography) self._active = entity