Merge branch 'main' into convolution-tests

healpix_convolution/kernels/gaussian.py

@@ -6,6 +6,16 @@
 from healpix_convolution.neighbours import neighbours
 
 
+def compute_ring(resolution, sigma, truncate, kernel_size):
+    if kernel_size is not None:
+        ring = int(kernel_size / 2)
+    else:
+        cell_distance = hp.nside2resol(2**resolution, arcmin=False)
+        ring = int((truncate * sigma / cell_distance) // 2)
+
+    return ring if ring >= 1 else 1
+
+
 def gaussian_function(distances, sigma, *, mask=None):
     sigma2 = sigma * sigma
     phi_x = np.exp(-0.5 / sigma2 * distances**2)
@@ -60,12 +70,7 @@ def gaussian_kernel(
 
     cell_ids = np.reshape(cell_ids, (-1,))
 
-    # TODO: figure out whether there is a better way of defining the units of `sigma`
-    if kernel_size is not None:
-        ring = int(kernel_size / 2)
-    else:
-        cell_distance = hp.nside2resol(2**resolution, arcmin=False)
-        ring = int((truncate * sigma / cell_distance) // 2)
+    ring = compute_ring(resolution, sigma, truncate, kernel_size)
 
     nb = neighbours(
         cell_ids, resolution=resolution, indexing_scheme=indexing_scheme, ring=ring

healpix_convolution/tests/test_kernels.py

@@ -48,6 +48,47 @@ def test_create_sparse(cell_ids, neighbours, weights):
     assert actual.shape == expected_shape
 
 
+def fit_polynomial(x, y, deg):
+    mask = y > 0
+    x_ = x[mask]
+    y_ = y[mask]
+
+    p = np.polynomial.Polynomial.fit(x_, np.log(y_), deg=deg)
+    return p.convert()
+
+
+def reconstruct_sigma(
+    cell_ids,
+    kernel,
+    *,
+    resolution,
+    indexing_scheme,
+    sigma,
+    truncate=4.0,
+    kernel_size=None,
+):
+    from healpix_convolution import angular_distances, neighbours
+
+    ring = np_kernels.gaussian.compute_ring(resolution, sigma, truncate, kernel_size)
+
+    nb = neighbours(
+        cell_ids, resolution=resolution, indexing_scheme=indexing_scheme, ring=ring
+    )
+    distances = angular_distances(
+        nb, resolution=resolution, indexing_scheme=indexing_scheme
+    )
+
+    _, distances_ = np_kernels.common.create_sparse(cell_ids, nb, distances)
+
+    x = distances_.todense()
+    y = kernel.todense()
+
+    polynomials = [
+        fit_polynomial(x[n, :], y[n, :], deg=2) for n in range(cell_ids.size)
+    ]
+    return np.array([np.sqrt(-1 / 2 / p.coef[2]) for p in polynomials])
+
+
 class TestGaussian:
     @pytest.mark.parametrize(
         ["cell_ids", "kwargs"],
@@ -93,6 +134,8 @@ def test_gaussian_kernel(self, cell_ids, kwargs):
         np.testing.assert_allclose(kernel_sum.todense(), 1)
 
         # try determining the sigma from the values for better tests
+        reconstructed_sigma = reconstruct_sigma(cell_ids, actual, **kwargs)
+        np.testing.assert_allclose(reconstructed_sigma, kwargs["sigma"])
 
     @pytest.mark.parametrize(
         ["cell_ids", "kwargs", "error", "pattern"],
@@ -248,3 +291,13 @@ def test_gaussian_kernel(self, obj, kwargs):
 
         assert actual.count() == actual.size
         np.testing.assert_allclose(kernel_sum.data.todense(), 1)
+
+        grid_info = obj_.dggs.grid_info
+        reconstructed_sigma = reconstruct_sigma(
+            obj.cell_ids.data,
+            actual.data,
+            resolution=grid_info.resolution,
+            indexing_scheme=grid_info.indexing_scheme,
+            **kwargs,
+        )
+        np.testing.assert_allclose(reconstructed_sigma, kwargs["sigma"])