Added error messages for fitting failure

Author: bobmyhill

burnman/optimize/nonlinear_fitting.py

@@ -339,13 +339,29 @@ def _update_beta(lmbda):
         return (current_params - new_params) / mod_params
 
     for n_it in range(max_lm_iterations):
-        # update the parameters with a LM iteration
-        f_delta_beta = _update_beta(lm_damping)
-        max_f = np.max(np.abs(f_delta_beta))
-        if verbose:
-            print(f"Iteration {n_it}: max param change = {max_f:.2e}")
-        if max_f < param_tolerance:
-            break
+        try:
+            # update the parameters with a LM iteration
+            f_delta_beta = _update_beta(lm_damping)
+            max_f = np.max(np.abs(f_delta_beta))
+
+            if np.isnan(max_f):
+                raise ValueError(
+                    "The Levenberg-Marquardt update for "
+                    f"Iteration {n_it} was non-numerical."
+                )
+
+            if verbose:
+                print(f"Iteration {n_it}: max param change = {max_f:.2e}")
+            if max_f < param_tolerance:
+                break
+        except Exception:
+            raise Exception(
+                f"During non-linear fitting, Iteration {n_it} produced an "
+                "exception. This is probably due to numerical failure of "
+                "the input model. "
+                "Consider imposing bounds on fitting or priors on your "
+                "parameter values to prevent this behaviour."
+            )
 
     J = model.jacobian
     r = model.weighted_residuals

tests/test_fitting.py

@@ -214,6 +214,46 @@ def test_fit_bounded_PVT_data(self):
         self.assertEqual(len(s[1]), 3)
         self.assertEqual(len(s[2]), 3)
 
+    def test_PTV_data_eos_numerical_failure(self):
+        fo = burnman.minerals.HP_2011_ds62.fo()
+
+        pressures = np.linspace(1.0e9, 2.0e9, 8)
+        temperatures = np.ones_like(pressures) * fo.params["T_0"]
+
+        fo.set_state(1.0e9, fo.params["T_0"])
+        V0 = fo.V
+        fo.set_state(2.0e9, fo.params["T_0"])
+        V1 = fo.V
+        volumes = V0 + (V1 - V0) * (pressures - 1.0e9) / 1.0e9
+        PTV = np.array([pressures, temperatures, volumes]).T
+        params = ["V_0", "K_0", "Kprime_0"]
+
+        try:
+            _ = burnman.eos_fitting.fit_PTV_data(fo, params, PTV, verbose=False)
+            self.assertTrue(False)
+        except Exception:
+            self.assertTrue(True)
+
+    def test_PTV_data_eos_exception(self):
+        fo = burnman.minerals.SLB_2011.forsterite()
+
+        pressures = np.linspace(1.0e9, 100.0e9, 8)
+        temperatures = np.ones_like(pressures) * fo.params["T_0"]
+
+        fo.set_state(1.0e9, fo.params["T_0"])
+        V0 = fo.V
+        fo.set_state(2.0e9, fo.params["T_0"])
+        V1 = fo.V
+        volumes = V0 + (V1 - V0) * (pressures - 1.0e9) / 1.0e9
+        PTV = np.array([pressures, temperatures, volumes]).T
+        params = ["V_0", "K_0", "Kprime_0"]
+
+        try:
+            _ = burnman.eos_fitting.fit_PTV_data(fo, params, PTV, verbose=False)
+            self.assertTrue(False)
+        except Exception:
+            self.assertTrue(True)
+
     def test_bounded_solution_fitting(self):
         solution = burnman.minerals.SLB_2011.mg_fe_olivine()
         solution.set_state(1.0e5, 300.0)