Merge pull request #373 from MiraGeoscience/GEOPY-2801

GEOPY-2801: PGI: Add checks for reference model present on inversion groups. Return clean GeoAppsError if not

Author: domfournier

simpeg_drivers/joint/joint_cross_gradient/driver.py

@@ -13,7 +13,9 @@
 
 from __future__ import annotations
 
+import sys
 from itertools import combinations
+from pathlib import Path
 
 from geoh5py.shared.utils import fetch_active_workspace
 from simpeg import maps
@@ -71,3 +73,8 @@ def get_regularization(self):
                     )
 
         return ComboObjectiveFunction(objfcts=reg_list, multipliers=multipliers)
+
+
+if __name__ == "__main__":
+    file = Path(sys.argv[1]).resolve()
+    JointCrossGradientDriver.start_dask_run(file)

simpeg_drivers/joint/joint_petrophysics/driver.py

@@ -10,7 +10,11 @@
 
 from __future__ import annotations
 
+import sys
+from pathlib import Path
+
 import numpy as np
+from geoapps_utils.utils.importing import GeoAppsError
 from geoh5py.shared.utils import fetch_active_workspace
 from simpeg import directives, maps, utils
 from simpeg.objective_function import ComboObjectiveFunction
@@ -158,6 +162,12 @@ def means(self) -> np.ndarray:
         """
         means = []
         for mapping in self.mapping:
+            if self.models.reference_model is None:
+                raise GeoAppsError(
+                    "A reference model must be set and active on each inversion driver "
+                    "to determine the means of the Gaussian mixture model.\n"
+                    "Please revise the input options of individual drivers."
+                )
             model_vec = mapping @ self.models.reference_model
             unit_mean = []
             for uid in self.geo_units:
@@ -220,3 +230,8 @@ def _overload_regularization(self, regularization: ComboObjectiveFunction):
             reg.alpha_s = 0.0
 
         return reg_list, multipliers
+
+
+if __name__ == "__main__":
+    file = Path(sys.argv[1]).resolve()
+    JointPetrophysicsDriver.start_dask_run(file)

simpeg_drivers/joint/joint_surveys/driver.py

@@ -11,7 +11,9 @@
 
 from __future__ import annotations
 
+import sys
 from logging import getLogger
+from pathlib import Path
 
 import numpy as np
 from geoh5py.shared.utils import fetch_active_workspace
@@ -144,3 +146,7 @@ def directives(self):
 
 JointSurveysDriver.n_values = InversionDriver.n_values
 JointSurveysDriver.mapping = InversionDriver.mapping
+
+if __name__ == "__main__":
+    file = Path(sys.argv[1]).resolve()
+    JointSurveysDriver.start_dask_run(file)

tests/run_tests/driver_joint_pgi_homogeneous_test.py

@@ -13,6 +13,8 @@
 from pathlib import Path
 
 import numpy as np
+import pytest
+from geoapps_utils import GeoAppsError
 from geoh5py.groups.property_group import GroupTypeEnum, PropertyGroup
 from geoh5py.objects import Octree, Points
 from geoh5py.workspace import Workspace
@@ -145,10 +147,10 @@ def test_homogeneous_fwr_run(
 def test_homogeneous_run(
     tmp_path: Path,
     max_iterations=1,
-    pytest=True,
+    use_pytest=True,
 ):
     workpath = tmp_path / "inversion_test.ui.geoh5"
-    if pytest:
+    if use_pytest:
         workpath = (
             tmp_path.parent / "test_homogeneous_fwr_run0" / "inversion_test.ui.geoh5"
         )
@@ -228,27 +230,28 @@ def test_homogeneous_run(
                     inducing_field_declination=INDUCING_FIELD[2],
                     data_object=survey,
                     starting_model=ref_model,
-                    reference_model=ref_model,
+                    reference_model=None,
                     tile_spatial=1,
                     tmi_channel=data,
                     tmi_uncertainty=5e0,
                 )
                 drivers.append(MagneticInversionDriver(params))
 
-        # Test if single group is valid
-        params = JointPetrophysicsOptions.build(
-            topography_object=topography,
-            geoh5=geoh5,
-            group_a=drivers[0].out_group,
-            mesh=global_mesh,
-            petrophysical_model=petrophysics,
-        )
-        driver = JointPetrophysicsDriver(params)
-        assert len(driver.data_misfit.objfcts) == 1
-        assert driver.data_misfit.multipliers == [1.0]
+            if len(drivers) == 1:
+                # Test if single group is valid
+                params = JointPetrophysicsOptions.build(
+                    topography_object=topography,
+                    geoh5=geoh5,
+                    group_a=drivers[0].out_group,
+                    mesh=global_mesh,
+                    petrophysical_model=petrophysics,
+                )
+                driver = JointPetrophysicsDriver(params)
+                assert len(driver.data_misfit.objfcts) == 1
+                assert driver.data_misfit.multipliers == [1.0]
 
         # Re-build full
-        params = JointPetrophysicsOptions.build(
+        joint_params = JointPetrophysicsOptions.build(
             topography_object=topography,
             geoh5=geoh5,
             group_a=drivers[0].out_group,
@@ -264,10 +267,22 @@ def test_homogeneous_run(
             initial_beta_ratio=1e2,
             max_global_iterations=max_iterations,
         )
-        driver = JointPetrophysicsDriver(params)
+        driver = JointPetrophysicsDriver(joint_params)
+
+        with pytest.raises(
+            GeoAppsError, match="A reference model must be set and active on each"
+        ):
+            _ = driver.means
+
+        # Re-instate
+        params.models.reference_model = ref_model
+        params.out_group = None
+        new_driver = MagneticInversionDriver(params)
+        joint_params.group_b = new_driver.out_group
+        driver = JointPetrophysicsDriver(joint_params)
         driver.run()
 
-    if pytest:
+    if use_pytest:
         with Workspace(driver.params.geoh5.h5file) as run_ws:
             output = get_inversion_output(
                 driver.params.geoh5.h5file, driver.out_group.uid
@@ -293,5 +308,5 @@ def test_homogeneous_run(
     test_homogeneous_run(
         Path("./"),
         max_iterations=20,
-        pytest=False,
+        use_pytest=False,
     )