GEOPY-1531: r2

simpeg/dask/data_misfit.py

@@ -14,7 +14,7 @@ def dask_call(self, m, f=None):
     Distributed :obj:`simpeg.data_misfit.L2DataMisfit.__call__`
     """
     R = self.W * self.residual(m, f=f)
-    phi_d = 0.5 * da.dot(R, R)
+    phi_d = da.dot(R, R)
     if not isinstance(phi_d, np.ndarray):
         return compute(self, phi_d)
     return phi_d

simpeg/dask/inverse_problem.py

@@ -139,7 +139,7 @@ def dask_evalFunction(self, m, return_g=True, return_H=True):
     phi_d = 0
     for (mult, objfct), pred in zip(self.dmisfit, self.dpred):
         residual = objfct.W * (objfct.data.dobs - pred)
-        phi_d += 0.5 * mult * np.vdot(residual, residual)
+        phi_d += mult * np.vdot(residual, residual)
 
     phi_d = np.asarray(phi_d)
     # print(self.dpred[0])

simpeg/dask/simulation.py

@@ -172,6 +172,7 @@ def Jmatrix(self):
     if getattr(self, "_Jmatrix", None) is None:
         if self.workers is None:
             self._Jmatrix = self.compute_J()
+            self._G = self._Jmatrix
         else:
             client = get_client()  # Assumes a Client already exists
 

simpeg/directives/__init__.py

@@ -120,6 +120,7 @@
     VectorInversion,
     SaveIterationsGeoH5,
     ProjectSphericalBounds,
+    ScaleMisfitMultipliers,
 )
 
 from .pgi_directives import (

simpeg/directives/directives.py

@@ -3426,3 +3426,70 @@ def endIter(self):
             for directive in directiveList:
                 if not isinstance(directive, SaveIterationsGeoH5):
                     directive.endIter()
+
+
+class ScaleMisfitMultipliers(InversionDirective):
+    """
+    Scale the misfits by a factor.
+    """
+
+    def __init__(self, chifact_target, out_group: SimPEGGroup, **kwargs):
+        self.last_beta = None
+        self.chifact_target = chifact_target
+        self.out_group = out_group
+
+        dirpath = Path(self.out_group.workspace.h5file).parent
+        self.filepath = dirpath / "ChiFactors.log"
+
+        super().__init__(**kwargs)
+
+    def initialize(self):
+        self.last_beta = self.invProb.beta
+        self.multipliers = self.invProb.dmisfit.multipliers
+
+        with open(self.filepath, "w", encoding="utf-8") as f:
+            f.write(
+                "Iterations\t"
+                + "\t".join(
+                    f"[{objfct.name}]" for objfct in self.invProb.dmisfit.objfcts
+                )
+            )
+            f.write("\n")
+
+    def endIter(self):
+        ratio = self.invProb.beta / self.last_beta
+        chi_factors = []
+        phi_ds = []
+        for objfct, pred in zip(self.invProb.dmisfit.objfcts, self.invProb.dpred):
+            residual = objfct.W * (objfct.data.dobs - pred)
+            phi_d = np.vdot(residual, residual)
+            chi_factors.append(phi_d / objfct.nD)
+            phi_ds.append(phi_d)
+
+        phi_ds = np.asarray(phi_ds)
+        chi_factors = np.asarray(chi_factors)
+        scalings = chi_factors / chi_factors.max()
+
+        # Force beta ratio scaling if below target
+        scalings[chi_factors < 1] *= ratio
+
+        # Normalize total phi_d with scalings
+        multipliers = (
+            self.multipliers
+            * scalings
+            * phi_ds.sum()
+            / (self.multipliers * phi_ds * scalings).sum()
+        )
+
+        with open(self.filepath, "a", encoding="utf-8") as f:
+            f.write(
+                f"{self.opt.iter}\t"
+                + "\t".join(
+                    f"{multi:.2e}*{chi:.2e}"
+                    for multi, chi in zip(multipliers, chi_factors)
+                )
+                + "\n"
+            )
+
+        self.invProb.dmisfit.multipliers = multipliers.tolist()
+        self.last_beta = self.invProb.beta