Merge pull request #102 from MTgeophysics/ak_occam2d_fixes_June2026

Update data.py and use of IRLs

Author: kujaku11

mtpy/modeling/occam2d/data.py

@@ -504,6 +504,9 @@ def read_response_file(self, response_fn=None, data_fn=None):
 
         entries = []
 
+        # initialise at zero
+        self.dataframe["t_zy"] = 0.0 + 0 * 1j
+
         for line in data_list:
             res_log = False
             # try:
@@ -556,35 +559,70 @@ def read_response_file(self, response_fn=None, data_fn=None):
             self.dataframe[key + "_model_error"] = 0.0
 
     def _group_df(self):
+
         df = self.dataframe
-        for station in np.unique(df["station"]):
-            for period in np.unique(df["period"]):
-                filt = np.all(
-                    [
-                        df["station"] == station,
-                        df["period"] == period,
-                    ],
-                    axis=0,
-                )
-                if np.any(filt):
-                    for key in ["res_xy", "res_yx", "phase_xy", "phase_yx"]:
-                        df.loc[filt, key] = np.unique(df.loc[filt, key])[0]
-                        df.loc[filt, key + "_model_error"] = np.unique(
-                            df.loc[filt, key + "_model_error"]
-                        )[0]
-                    tx_real_vals = np.unique(np.real(df.loc[filt, "t_zy"]))
-                    tx_imag_vals = np.unique(np.imag(df.loc[filt, "t_zy"]))
-
-                    if np.any(tx_real_vals != 0) or np.any(tx_imag_vals != 0):
-                        tx_real_val, tx_imag_val = 0.0, 0.0
-                        if len(tx_real_vals[tx_real_vals != 0]) > 0:
-                            tx_real_val = tx_real_vals[tx_real_vals != 0][0]
-                        if len(tx_imag_vals[tx_imag_vals != 0]) > 0:
-                            tx_imag_val = tx_imag_vals[tx_imag_vals != 0][0]
-
-                        df.loc[filt, "t_zy"] = tx_real_val + 1j * tx_imag_val
-
-                    df.drop(df.index[filt][1:], inplace=True)
+
+        # 1. Helper function for the complex number logic
+        def process_t_zy(series):
+            # Drop NaNs and convert to a standard Python list
+            vals = series.dropna().tolist()
+
+            reals = []
+            imags = []
+
+            for v in vals:
+                try:
+                    # Force the value into a standard complex type safely
+                    c = complex(v)
+                    if c.real != 0.0:
+                        reals.append(c.real)
+                    if c.imag != 0.0:
+                        imags.append(c.imag)
+                except (TypeError, ValueError):
+                    # Safely ignore any weird strings or corrupted data
+                    continue
+
+            # sort the values, so [0] grabs the minimum.
+            r = sorted(set(reals))[0] if reals else 0.0
+            i = sorted(set(imags))[0] if imags else 0.0
+
+            return complex(r, i)
+
+        # 2. Identify the specific columns
+        keys = ["res_xy", "res_yx", "phase_xy", "phase_yx"]
+        error_keys = [f"{k}_model_error" for k in keys]
+        special_cols = keys + error_keys
+
+        # 3. Build the aggregation dictionary with Python Lambdas
+        agg_funcs = {}
+        for col in df.columns:
+            if col in ["station", "period"]:
+                continue
+            elif col == "t_zy":
+                agg_funcs[col] = process_t_zy
+            elif col in special_cols:
+                # Bypass Cython by using a lambda for minimum
+                agg_funcs[col] = lambda x: x.min()
+            else:
+                # Bypass Cython by explicitly asking for the first index
+                agg_funcs[col] = lambda x: x.iloc[0]
+
+        # 4. Perform the GroupBy operation
+        self.dataframe = df.groupby(["station", "period"], as_index=False).agg(
+            agg_funcs
+        )
+
+    # def _group_df(self):
+    #     print("grouping df")
+
+    #     # A custom Python function that does exactly what .first() does:
+    #     # It drops NaNs, and grabs the first available value.
+    #     def get_first_valid(series):
+    #         valid_data = series.dropna()
+    #         return valid_data.iloc[0] if not valid_data.empty else np.nan
+
+    #     # .agg() with a custom lambda forces pandas to stay in Python space
+    #     self.dataframe = self.dataframe.groupby(['station', 'period'], as_index=False).agg(get_first_valid)
 
     def _get_model_mode_from_data(self, res_log):
         """Get inversion mode from the data.

mtpy/modeling/simpeg/recipes/inversion_2d.py

@@ -120,6 +120,9 @@ def __init__(
         self.beta_cooling_rate = 1
 
         self.target_misfit_chi_factor = 1
+        self.irls_chifact_start = 2.0
+        self.irls_chifact_target = 1.0
+        self.irls_cooling_factor = 1.5
 
         self.saved_model_outputs = directives.SaveOutputDictEveryIteration()
         self.saved_model_outputs.outDict = {}
@@ -314,7 +317,7 @@ def regularization(self):
         )
 
         if self.use_irls:
-            reg.norms = np.c_[self.p_s, self.p_y, self.p_z]
+            reg.norms = [self.p_s, self.p_y, self.p_z]
 
         return reg
 
@@ -395,10 +398,14 @@ def directives(self):
         """
 
         if self.use_irls:
-            IRLS = directives.Update_IRLS(
-                max_irls_iterations=self.max_iteration_irls,
-                minGNiter=self.minimum_gauss_newton_iterations,
+            IRLS = directives.UpdateIRLS(
+                max_irls_iterations=self.max_iterations_irls,
+                irls_cooling_factor=self.irls_cooling_factor,
                 f_min_change=self.f_min_change,
+                cooling_rate=self.beta_cooling_rate,
+                cooling_factor=self.beta_cooling_factor,
+                chifact_start=self.irls_chifact_start,
+                chifact_target=self.irls_chifact_target,
             )
             return [
                 IRLS,
@@ -410,7 +417,7 @@ def directives(self):
                 self.starting_beta,
                 self.beta_schedule,
                 self.saved_model_outputs,
-                # self.target_misfit,
+                self.target_misfit,
             ]
 
     def run_inversion(self):