GEOPY-2544: Add mechanics to re-start inversion if results are present in the SimPEGGroups

simpeg/dask/inverse_problem.py

@@ -57,7 +57,7 @@ def get_dpred(self, m, f=None, return_residuals=False):
 def dask_evalFunction(self, m, return_g=True, return_H=True):
     """evalFunction(m, return_g=True, return_H=True)"""
 
-    if not np.allclose(self.model, m):
+    if not np.allclose(self.model, m) or getattr(self, "residuals", None) is None:
         self.model = m
         self.dpred, self.residuals = self.get_dpred(m, return_residuals=True)
 

simpeg/directives/_directives.py

@@ -2691,11 +2691,13 @@ def __init__(
         path: pathlib.Path | None = None,
         nesting: list[list] | None = None,
         target_chi: float = 1.0,
+        cooling_factor: float = 2.0,
         headers: list[str] | None = None,
         **kwargs,
     ):
-        self.last_beta = None
+
         self.chi_factors = None
+        self.cooling_factor = cooling_factor
         self.target_chi = target_chi
         self.nesting = nesting
         self.headers = headers
@@ -2706,7 +2708,9 @@ def __init__(
         self.filepath = path
 
         super().__init__(**kwargs)
-
+        self._last_beta = None
+        self._multipliers: np.ndarray | None = None
+        self._scalings: np.ndarray | None = None
         self._log_array: np.ndarray | None = None
 
     @property
@@ -2737,11 +2741,7 @@ def append_sub_indices(elements, header):
         return self._log_array
 
     def initialize(self):
-        self.last_beta = self.invProb.beta
-        self.multipliers = self.invProb.dmisfit.multipliers
-        self.scalings = np.ones_like(self.multipliers)  # Everyone gets a fair chance
         self.misfit_tree_indices = self.parse_by_nested_levels(self.nesting)
-
         self.write_log()
 
     def scale_by_level(
@@ -2806,7 +2806,11 @@ def scale_by_level(
         return scaling_vector
 
     def endIter(self):
-        ratio = self.invProb.beta / self.last_beta
+        if self._last_beta is None:
+            ratio = 1.0 / self.cooling_factor
+        else:
+            ratio = self.invProb.beta / self._last_beta
+
         nested_residuals = self.parse_by_nested_levels(
             self.nesting, self.invProb.residuals
         )
@@ -2820,7 +2824,7 @@ def endIter(self):
 
         # Normalize total phi_d with scalings
         self.invProb.dmisfit.multipliers = self.multipliers * self.scalings
-        self.last_beta = self.invProb.beta
+        self._last_beta = self.invProb.beta
 
         # Log the scaling factors
         self.write_log()
@@ -2890,9 +2894,39 @@ def write_log(self):
                 f,
                 self.log_array,
                 header="Iterations - Scaling per misfit",
-                fmt=["%d"] + ["%0.2e"] * (len(self._log_array.dtype) - 1),
+                fmt=["%d"] + ["%0.3e"] * (len(self._log_array.dtype) - 1),
             )
 
+    @property
+    def multipliers(self) -> np.ndarray:
+        """Static misfit multipliers set at the start of the inversion."""
+        if self._multipliers is None:
+            self._multipliers = self.invProb.dmisfit.multipliers
+
+        return self._multipliers
+
+    @multipliers.setter
+    def multipliers(self, value):
+        if not isinstance(value, np.ndarray):
+            raise TypeError("Multipliers must be a numpy array.")
+        self._multipliers = value
+
+    @property
+    def scalings(self) -> np.ndarray:
+        """Get the current scaling factors."""
+        if self._scalings is None:
+            self._scalings = np.ones_like(
+                self.invProb.dmisfit.multipliers
+            )  # Everyone gets a fair chance
+
+        return self._scalings
+
+    @scalings.setter
+    def scalings(self, value):
+        if not isinstance(value, np.ndarray):
+            raise TypeError("Scaling factors must be a numpy array.")
+        self._scalings = value
+
 
 def compute_JtJdiags(data_misfit, m):
     if hasattr(data_misfit, "getJtJdiag"):

simpeg/directives/_save_geoh5.py

@@ -55,7 +55,8 @@ def initialize(self):
         if self.open_geoh5 and not getattr(self._workspace, "_geoh5", None):
             self._workspace.open(mode="r+")
 
-        self.write(0)
+        if getattr(self.opt, "iter", 0) == 0:
+            self.write(0)
 
         if self.close_geoh5:
             self._workspace.close()
@@ -279,25 +280,27 @@ def write(self, iteration: int, values: list[np.ndarray] = None):  # flake8: noq
                     channel_name, base_name = self.get_names(
                         component, label, iteration
                     )
-
+                    data_type = self.data_type[component].get(ii, None)
                     data = h5_object.add_data(
                         {
                             channel_name: {
                                 "association": self.association,
                                 "values": values,
+                                "entity_type": data_type,
                             }
                         }
                     )
-                    # Re-assign the data type
-                    if ii not in self.data_type[component].keys():
+                    if data_type is None:
+                        type_name = f"{self._attribute_type}"
+
+                        if component:
+                            type_name += f"_{component}"
+
+                        if label:
+                            type_name += f"_{label}"
+
                         self.data_type[component][ii] = data.entity_type
-                        type_name = f"{self._attribute_type}_{component}" + f"_{label}"
                         data.entity_type.name = type_name
-                    else:
-                        data.entity_type = w_s.find_type(
-                            self.data_type[component][ii].uid,
-                            type(self.data_type[component][ii]),
-                        )
 
 
 class SaveModelGeoH5(SaveArrayGeoH5):
@@ -421,9 +424,6 @@ def write(self, iteration: int, **_):
 
             with open(filepath, "a", encoding="utf-8") as file:
                 date_time = datetime.now().strftime("%b-%d-%Y:%H:%M:%S")
-
-                if len(log) == 2:  # First iteration with 0th iter
-                    file.write(f"{0} " + " ".join(log[0]) + f" {date_time}\n")
                 file.write(f"{iteration-1} " + " ".join(log[-1]) + f" {date_time}\n")
 
         self.save_log()

simpeg/directives/_sim_directives.py

@@ -397,9 +397,15 @@ def endIter(self):
         """
         End of iteration update.
         """
+        self.scale_on_current_model(self.opt.xc)
+
+    def scale_on_current_model(self, model):
+        """
+        Scale the cross-gradient term based on the current iteration.
+        """
         max_deriv = []
         for _, wire in self.cross_gradient.wire_map.maps:
-            component = wire * self.opt.xc
+            component = wire * model
             max_deriv.append(
                 (component.max() - component.min())
                 / self.cross_gradient.regularization_mesh.base_length**2.0

simpeg/optimization.py

@@ -220,15 +220,15 @@ class IterationPrinters(object):
         "title": "f",
         "value": lambda M: M.f,
         "width": 10,
-        "format": lambda v: f"{v:1.2e}",
+        "format": lambda v: f"{v:1.3e}",
     }
     norm_g = {
         "title": "|proj(x-g)-x|",
         "value": lambda M: (
             None if M.iter == 0 else norm(M.projection(M.xc - M.g) - M.xc)
         ),
         "width": 15,
-        "format": lambda v: f"{v:1.2e}",
+        "format": lambda v: f"{v:1.3e}",
     }
     totalLS = {
         "title": "LS",
@@ -247,7 +247,7 @@ class IterationPrinters(object):
         "title": "ft",
         "value": lambda M: M._LS_ft,
         "width": 10,
-        "format": lambda v: f"{v:1.2e}",
+        "format": lambda v: f"{v:1.3e}",
     }
     LS_t = {
         "title": "t",
@@ -259,14 +259,14 @@ class IterationPrinters(object):
         "title": "f + alp*g.T*p",
         "value": lambda M: M.f + M.LSreduction * M._LS_descent,
         "width": 16,
-        "format": lambda v: f"{v:1.2e}",
+        "format": lambda v: f"{v:1.3e}",
     }
     LS_WolfeCurvature = {
         "title": "alp*g.T*p",
         "str": "%d :    ft     = %1.4e >= alp*descent     = %1.4e",
         "value": lambda M: M.LScurvature * M._LS_descent,
         "width": 16,
-        "format": lambda v: f"{v:1.2e}",
+        "format": lambda v: f"{v:1.3e}",
     }
 
     itType = {
@@ -298,44 +298,44 @@ class IterationPrinters(object):
         "title": "beta",
         "value": lambda M: M.parent.beta,
         "width": 10,
-        "format": lambda v: f"{v:1.2e}",
+        "format": lambda v: f"{v:1.3e}",
     }
     phi_d = {
         "title": "phi_d",
         "value": lambda M: M.parent.phi_d,
         "width": 10,
-        "format": lambda v: f"{v:1.2e}",
+        "format": lambda v: f"{v:1.3e}",
     }
     phi_m = {
         "title": "phi_m",
         "value": lambda M: M.parent.phi_m,
         "width": 10,
-        "format": lambda v: f"{v:1.2e}",
+        "format": lambda v: f"{v:1.3e}",
     }
 
     phi_s = {
         "title": "phi_s",
         "value": lambda M: M.parent.phi_s,
         "width": 10,
-        "format": lambda v: f"{v:1.2e}",
+        "format": lambda v: f"{v:1.3e}",
     }
     phi_x = {
         "title": "phi_x",
         "value": lambda M: M.parent.phi_x,
         "width": 10,
-        "format": lambda v: f"{v:1.2e}",
+        "format": lambda v: f"{v:1.3e}",
     }
     phi_y = {
         "title": "phi_y",
         "value": lambda M: M.parent.phi_y,
         "width": 10,
-        "format": lambda v: f"{v:1.2e}",
+        "format": lambda v: f"{v:1.3e}",
     }
     phi_z = {
         "title": "phi_z",
         "value": lambda M: M.parent.phi_z,
         "width": 10,
-        "format": lambda v: f"{v:1.2e}",
+        "format": lambda v: f"{v:1.3e}",
     }
 
     iterationCG = {
@@ -349,14 +349,14 @@ class IterationPrinters(object):
         "title": "CG |Ax-b|/|b|",
         "value": lambda M: getattr(M, "cg_rel_resid", None),
         "width": 15,
-        "format": lambda v: f"{v:1.2e}",
+        "format": lambda v: f"{v:1.3e}",
     }
 
     iteration_CG_abs_residual = {
         "title": "CG |Ax-b|",
         "value": lambda M: getattr(M, "cg_abs_resid", None),
         "width": 11,
-        "format": lambda v: f"{v:1.2e}",
+        "format": lambda v: f"{v:1.3e}",
     }
 
 
@@ -397,7 +397,8 @@ class Minimize(object):
 
     def __init__(self, **kwargs):
         set_kwargs(self, **kwargs)
-
+        self.iter = 0
+        self.f0 = None
         self.stoppersLS = [
             StoppingCriteria.armijoGoldstein,
             StoppingCriteria.iterationLS,
@@ -546,8 +547,6 @@ def startup(self, x0: np.ndarray) -> None:
         x0 : numpy.ndarray
             initial x
         """
-
-        self.iter = 0
         self.iterLS = 0
         self.stopNextIteration = False
 
@@ -644,7 +643,7 @@ def finish(self) -> None:
         pass
 
     def stoppingCriteria(self, inLS: bool = False) -> bool:
-        if self.iter == 0:
+        if self.f0 is None:
             self.f0 = self.f
             self.g0 = self.g
         return check_stoppers(self, self.stoppers if not inLS else self.stoppersLS)