added AIC and BIC

Author: ramav87

sidpy/proc/fitter_refactor.py

@@ -181,7 +181,7 @@ def residuals(p, x, y):
                     raise ValueError(f"Unknown cov_mode={cov_mode!r}.")
 
             if return_metrics:
-                out_parts.append(np.array([np.nan, np.nan], dtype=np.float64))
+                out_parts.append(np.array([np.nan, np.nan, np.nan, np.nan], dtype=np.float64))
 
             return np.hstack(out_parts) if len(out_parts) > 1 else params
 
@@ -199,6 +199,8 @@ def residuals(p, x, y):
 
         # --- Metrics (default) appending ---
         if return_metrics:
+            n = int(y_input.size)
+            k = int(params.size)
             sse = float(np.sum(res.fun ** 2))
             rmse = float(np.sqrt(sse / y_input.size)) if y_input.size > 0 else np.nan
             y_mean = float(np.mean(y_input))
@@ -208,7 +210,17 @@ def residuals(p, x, y):
             else:
                 r2 = 1.0 - (sse / sst)
 
-            out_parts.append(np.array([r2, rmse], dtype=np.float64))
+            # AIC/BIC (up to an additive constant)
+            if n > 0 and sse > 0.0:
+                ll_term = n * np.log(sse / n)
+                aic = float(ll_term + 2.0 * k)
+                bic = float(ll_term + k * np.log(n))
+            else:
+                aic = np.nan
+                bic = np.nan
+                
+            # Order: [r2, rmse, aic, bic]
+            out_parts.append(np.array([r2, rmse, aic, bic], dtype=np.float64))
 
         # --- Covariance (optional) appending ---
         if cov_mode is not None:
@@ -504,7 +516,7 @@ def do_fit(self, guesses=None, use_kmeans=False, n_clusters=10,
         else:
             raise ValueError(f"Unknown cov_mode={cov_mode!r}. Use None, 'full', 'diag', or 'stderr'.")
 
-        metrics_size = 2 if return_metrics else 0
+        metrics_size = 4 if return_metrics else 0
         out_dim = self.num_params + cov_size + metrics_size
 
 
@@ -627,19 +639,21 @@ def set_spatial_dims(dset):
 
         # --- 2) Metrics dataset ---
         if return_metrics:
-            metrics_payload = fit_dask_array[..., cursor:cursor + 2]
+            metrics_payload = fit_dask_array[..., cursor:cursor + 4]
+            cursor += 4
 
             sid_metrics = sid.Dataset.from_array(metrics_payload, title=f"{self.dataset.title}_fit_metrics")
             set_spatial_dims(sid_metrics)
             sid_metrics.set_dimension(
                 len(self.spat_dims),
-                sid.Dimension(np.arange(2), name='metrics', quantity='index', dimension_type='spectral')
+                sid.Dimension(np.arange(4), name='metrics', quantity='index', dimension_type='spectral')
             )
             sid_metrics.metadata = dict(self.metadata).copy()
-            sid_metrics.metadata.setdefault("fit_parameters", {}).update({"metrics": ["r2", "rmse"]})
+            sid_metrics.metadata.setdefault("fit_parameters", {}).update({"metrics": ["r2", "rmse", "aic", "bic"]})
             sid_metrics.provenance = {'sidpy': {'generated_from': self.dataset.title, 'parent_fit': sid_params.title}}
             out.append(sid_metrics)
 
+
         # --- 3) Covariance dataset (optional) ---
         if cov_size > 0:
             cov_payload = fit_dask_array[..., cursor:cursor + cov_size]

tests/proc/test_fitter_refactor.py

@@ -357,15 +357,52 @@ def test_metrics_perfect_fit_r2_one_rmse_zero(self):
         res_params, res_metrics = fitter.do_fit(use_kmeans=False, loss='linear', return_metrics=True)
 
         assert res_params.shape == (nx, ny, 6)
-        assert res_metrics.shape == (nx, ny, 2)
+        assert res_metrics.shape == (nx, ny, 4)
 
         m = np.asarray(res_metrics)
         r2 = m[..., 0]
         rmse = m[..., 1]
-
+        assert np.isfinite(m[..., 2]).any()  # AIC
+        assert np.isfinite(m[..., 3]).any()  # BIC
         assert np.nanmin(r2) > 0.999
         assert np.nanmax(rmse) < 1e-6
 
+    def test_aic_bic_increase_with_noise_on_same_model():
+        import sidpy as sid
+        rng = np.random.default_rng(0)
+
+        nx, ny = 2, 2
+        nkx, nky = 16, 16
+        x_vec = np.linspace(-10, 10, nkx)
+        y_vec = np.linspace(-10, 10, nky)
+        true_p = np.array([1.8, 1.0, -1.0, 2.4, 1.7, 0.05], dtype=float)
+
+        clean2d = gaussian_2d((x_vec, y_vec), *true_p)
+        clean = np.tile(clean2d, (nx, ny, 1, 1))
+        sigma = 0.02 * np.max(np.abs(clean2d))
+        noisy = clean + rng.normal(0.0, sigma, size=clean.shape)
+
+        def mk(arr, title):
+            d = sid.Dataset.from_array(arr, title=title)
+            d.set_dimension(0, sid.Dimension(np.arange(nx), name='X', dimension_type='spatial'))
+            d.set_dimension(1, sid.Dimension(np.arange(ny), name='Y', dimension_type='spatial'))
+            d.set_dimension(2, sid.Dimension(x_vec, name='kx', dimension_type='spectral'))
+            d.set_dimension(3, sid.Dimension(y_vec, name='ky', dimension_type='spectral'))
+            return d
+
+        def run(arr):
+            fitter = SidpyFitterRefactor(mk(arr, "tmp"), gaussian_2d, gaussian_2d_guess, ind_dims=(2, 3), num_params=6)
+            fitter.setup_calc()
+            _, met = fitter.do_fit(use_kmeans=False, loss='linear')
+            m = np.asarray(met)
+            return np.nanmedian(m[..., 2]), np.nanmedian(m[..., 3])  # AIC, BIC
+
+        aic_clean, bic_clean = run(clean)
+        aic_noisy, bic_noisy = run(noisy)
+
+        assert aic_noisy > aic_clean
+        assert bic_noisy > bic_clean
+
 
     def test_metrics_noisy_fit_lower_r2_higher_rmse_than_clean(self):