Merge pull request #337 from MiraGeoscience/GEOPY-2659

GEOPY-2659: MVI: Add small value to reference amplitude if reference are used

Author: domfournier

simpeg_drivers/components/models.py

@@ -233,6 +233,7 @@ def reference_model(self) -> np.ndarray | None:
                 self.reference_inclination,
                 self.reference_declination,
             )
+            ref_model += 1e-8  # To avoid zeroing the angles
             ref_model = (field_vecs.T * ref_model).flatten()
 
         return ref_model
@@ -271,11 +272,11 @@ def lower_bound(self) -> np.ndarray | None:
         else:
             bound_model = self._lower_bound.model
 
-        if (
-            self.driver.params.inversion_type == "magnetic vector"
-            and self._upper_bound.model is not None
-        ):
-            bound_model = -self._upper_bound.model
+        if self.driver.params.inversion_type == "magnetic vector":
+            bound_model = None
+
+            if self._upper_bound.model is not None:
+                bound_model = -self._upper_bound.model
 
         if bound_model is None:
             lbound = np.full(self.n_active, -np.inf)

tests/run_tests/driver_mvi_test.py

@@ -20,6 +20,7 @@
 from geoh5py.groups.property_group import GroupTypeEnum
 from geoh5py.objects import Curve
 from geoh5py.workspace import Workspace
+from simpeg.utils.mat_utils import cartesian2amplitude_dip_azimuth
 
 from simpeg_drivers.components.factories import DirectivesFactory
 from simpeg_drivers.potential_fields import (
@@ -213,7 +214,7 @@ def test_magnetic_vector_reference(
             tmi_uncertainty=5.0,
             data_object=components.survey,
             starting_model=components.model,
-            reference_model=1.0,
+            reference_model=0.0,
             reference_inclination=30,
             reference_declination=0,
         )
@@ -224,13 +225,10 @@ def test_magnetic_vector_reference(
     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,
-    )
+    ref_model = driver.models.reference_model
+    ref_spherical = cartesian2amplitude_dip_azimuth(ref_model.reshape(-1, 3, order="F"))
+    np.allclose(ref_spherical[0, 1], 30)
+    np.allclose(ref_spherical[0, 2], 0)
 
 
 if __name__ == "__main__":