Issue a warning if exiting without convergence

Author: leouieda

harmonica/_euler_methods.py

@@ -8,6 +8,8 @@
 Classes for Euler-based source location estimation.
 """
 
+import warnings
+
 import numpy as np
 import scipy as sp
 import verde.base as vdb
@@ -391,7 +393,7 @@ def _fit(self, coordinates, data, weights):
         data_misfit = np.linalg.norm(residuals * data_weights)
         merit = data_misfit + self.euler_misfit_balance * euler_misfit
         for _ in range(self.max_iterations):
-            parameter_step, data_step, cofactor_new = self._newton_step(
+            parameter_step, data_step, cofactor = self._newton_step(
                 coordinates,
                 data_observed,
                 data_predicted,
@@ -404,25 +406,21 @@ def _fit(self, coordinates, data, weights):
             # Update metrics
             euler = self._eulers_equation(coordinates, data_predicted, parameters)
             residuals = data_observed - data_predicted
-            new_euler_misfit = np.linalg.norm(euler)
-            new_data_misfit = np.linalg.norm(residuals * data_weights)
-            new_merit = new_data_misfit + self.euler_misfit_balance * new_euler_misfit
-            # Check if was increasing the merit function, which means this solution is
-            # worse than the last.
-            if new_merit > merit:
-                # If it's bad, walk back the last step
-                data_predicted -= data_step
-                parameters -= parameter_step
-                break
-            cofactor = cofactor_new
-            # Update tracked metrics
-            euler_misfit = new_euler_misfit
-            data_misfit = new_data_misfit
+            euler_misfit = np.linalg.norm(euler)
+            data_misfit = np.linalg.norm(residuals * data_weights)
+            new_merit = data_misfit + self.euler_misfit_balance * euler_misfit
             merit_change = abs((merit - new_merit) / merit)
             merit = new_merit
             # Check for convergence
             if merit_change < self.tol:
                 break
+        else:
+            message = (
+                "Euler Inversion exited because maximum number of iterations was reached"
+                " and not because the algorithm converged. Consider increasing the"
+                " maximum number of iterations."
+            )
+            warnings.warn(message, stacklevel=2)
         # Save output attributes
         self.location_ = parameters[:3]
         self.base_level_ = parameters[3]

harmonica/tests/test_euler_methods.py

@@ -59,8 +59,7 @@ def test_euler_inversion_with_numeric_derivatives(euler):
     npt.assert_allclose(euler.location_[:2], dipole_coordinates[:2], rtol=0.01)
     npt.assert_allclose(euler.location_[2], dipole_coordinates[2], rtol=0.05)
     npt.assert_allclose(euler.base_level_, true_base_level, rtol=0.01)
-    if hasattr(euler, "structural_index_"):
-        assert euler.structural_index_ == 3
+    assert euler.structural_index_ == 3
 
 
 @pytest.mark.parametrize(
@@ -107,8 +106,45 @@ def test_euler_inversion_with_analytic_derivatives(euler):
     npt.assert_allclose(euler.location_[:2], masses_coordinates[:2], rtol=0.01)
     npt.assert_allclose(euler.location_[2], masses_coordinates[2], rtol=0.05)
     npt.assert_allclose(euler.base_level_, base_level, rtol=0.01)
-    if hasattr(euler, "structural_index_"):
-        assert euler.structural_index_ == 2
+    assert euler.structural_index_ == 2
+
+
+
+@pytest.mark.parametrize(
+    "euler",
+    [
+        EulerInversion(structural_index=3, max_iterations=1),
+        EulerInversion(max_iterations=1),
+    ],
+    ids=("fixed_SI", "estimate_SI"),
+)
+def test_euler_inversion_convergence_warning(euler):
+    "Check that a warning is raised when the method exits without convergence"
+    # Add dipole source
+    dipole_coordinates = (10e3, 15e3, -10e3)
+    inc, dec = -40, 15
+    dipole_moments = magnetic_angles_to_vec(1.0e14, inc, dec)
+    region = [-100e3, 100e3, -80e3, 80e3]
+    coordinates = vd.grid_coordinates(region, spacing=500, extra_coords=500)
+    b = dipole_magnetic(coordinates, dipole_coordinates, dipole_moments, field="b")
+    # Add a fixed base level
+    true_base_level = 200  # nT
+    anomaly = total_field_anomaly(b, inc, dec) + true_base_level
+    anomaly += np.random.default_rng(42).normal(0, 20, anomaly.shape)
+    grid = vd.make_xarray_grid(
+        coordinates, anomaly, data_names="tfa", extra_coords_names="upward"
+    )
+    grid["d_east"] = derivative_easting(grid.tfa)
+    grid["d_north"] = derivative_northing(grid.tfa)
+    grid["d_up"] = derivative_upward(grid.tfa)
+    table = vd.grid_to_table(grid)
+    coordinates = (table.easting, table.northing, table.upward)
+    with pytest.warns(UserWarning, match="Euler Inversion exited"):
+        euler.fit(
+            (table.easting, table.northing, table.upward),
+            (table.tfa, table.d_east, table.d_north, table.d_up),
+        )
+
 
 
 def test_euler_deconvolution_with_numeric_derivatives():
@@ -139,8 +175,6 @@ def test_euler_deconvolution_with_numeric_derivatives():
     )
     npt.assert_allclose(euler.location_, dipole_coordinates, rtol=0.01)
     npt.assert_allclose(euler.base_level_, true_base_level, rtol=0.01)
-    if hasattr(euler, "structural_index_"):
-        npt.assert_allclose(euler.structural_index_, 3, rtol=0.01)
 
 
 def test_euler_deconvolution_with_analytic_derivatives():
@@ -172,5 +206,3 @@ def test_euler_deconvolution_with_analytic_derivatives():
     )
     npt.assert_allclose(euler.location_, masses_coordinates, atol=1.0e-3, rtol=1.0e-3)
     npt.assert_allclose(euler.base_level_, 0.0, atol=1.0e-3, rtol=1.0e-3)
-    if hasattr(euler, "structural_index_"):
-        npt.assert_allclose(euler.structural_index_, 2, rtol=0.01)