From 06328e7f7532918be6fa378b651f2538801b3b36 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Thu, 11 Sep 2025 12:15:28 -0700 Subject: [PATCH 1/2] Fix options.model reference --- simpeg_drivers/components/factories/directives_factory.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/simpeg_drivers/components/factories/directives_factory.py b/simpeg_drivers/components/factories/directives_factory.py index 6301d2ff3..40e369965 100644 --- a/simpeg_drivers/components/factories/directives_factory.py +++ b/simpeg_drivers/components/factories/directives_factory.py @@ -332,9 +332,11 @@ def vector_inversion_directive(self): """Directive to update vector model.""" if self._vector_inversion_directive is None and "vector" in self.factory_type: reference_angles = ( - getattr(self.driver.params, "reference_model", None) is not None, - getattr(self.driver.params, "reference_inclination", None) is not None, - getattr(self.driver.params, "reference_declination", None) is not None, + getattr(self.driver.params.models, "reference_model", None) is not None, + getattr(self.driver.params.models, "reference_inclination", None) + is not None, + getattr(self.driver.params.models, "reference_declination", None) + is not None, ) self._vector_inversion_directive = directives.VectorInversion( From 10ce18d6340329e87ed99871310268e79814b2a2 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Thu, 11 Sep 2025 14:14:00 -0700 Subject: [PATCH 2/2] Fix reference model output. Add unitest --- simpeg_drivers/components/models.py | 14 ++++---- tests/run_tests/driver_mvi_test.py | 53 +++++++++++++++++++++++++++++ 2 files changed, 59 insertions(+), 8 deletions(-) diff --git a/simpeg_drivers/components/models.py b/simpeg_drivers/components/models.py index eb3530202..7c5049841 100644 --- a/simpeg_drivers/components/models.py +++ b/simpeg_drivers/components/models.py @@ -212,18 +212,16 @@ def reference_model(self) -> np.ndarray | None: ref_model = mref.copy() - if ( - self.is_sigma - and self.driver.params.models.model_type == "Resistivity (Ohm-m)" - ): - ref_model = 1 / ref_model + if self.is_sigma: + if self.driver.params.models.model_type == "Resistivity (Ohm-m)": + ref_model = 1 / ref_model - ref_model = np.log(ref_model) if self.is_sigma else ref_model + ref_model = np.log(ref_model) if self.is_vector: field_vecs = dip_azimuth2cartesian( - self.starting_inclination, - self.starting_declination, + self.reference_inclination, + self.reference_declination, ) ref_model = (field_vecs.T * ref_model).flatten() diff --git a/tests/run_tests/driver_mvi_test.py b/tests/run_tests/driver_mvi_test.py index 5308d2635..5db8981b9 100644 --- a/tests/run_tests/driver_mvi_test.py +++ b/tests/run_tests/driver_mvi_test.py @@ -20,6 +20,7 @@ from geoh5py.objects import Curve from geoh5py.workspace import Workspace +from simpeg_drivers.components.factories import DirectivesFactory from simpeg_drivers.potential_fields import ( MVIForwardOptions, MVIInversionOptions, @@ -180,6 +181,58 @@ def test_magnetic_vector_run( check_target(output, target_mvi_run) +def test_magnetic_vector_reference( + tmp_path: Path, + n_grid_points=3, + refinement=(2,), +): + # Run the forward + opts = SyntheticsComponentsOptions( + method="magnetic_vector", + survey=SurveyOptions( + n_stations=n_grid_points, n_lines=n_grid_points, drape=5.0 + ), + mesh=MeshOptions(refinement=refinement), + model=ModelOptions(anomaly=0.05), + ) + with get_workspace(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + + tmi = components.survey.add_data( + {"tmi": {"values": np.random.randn(components.survey.n_vertices)}} + ) + inducing_field = (50000.0, 90.0, 0.0) + params = MVIInversionOptions.build( + geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, + inducing_field_strength=inducing_field[0], + inducing_field_inclination=inducing_field[1], + inducing_field_declination=inducing_field[2], + tmi_channel=tmi, + tmi_uncertainty=5.0, + data_object=components.survey, + starting_model=components.model, + reference_model=1.0, + reference_inclination=30, + reference_declination=0, + ) + driver = MVIInversionDriver(params) + + directives = DirectivesFactory(driver) + assert np.all(directives.vector_inversion_directive.reference_angles) + assert np.all(driver.models.reference_inclination == 30) + assert np.all(driver.models.reference_declination == 0) + + np.allclose( + np.kron( + np.r_[0, np.cos(-np.deg2rad(30)), np.sin(-np.deg2rad(30))], + np.ones(driver.models.n_active), + ), + driver.models.reference_model, + ) + + if __name__ == "__main__": # Full run test_magnetic_vector_fwr_run(Path("./"), n_grid_points=20, refinement=(4, 8))