Merge pull request #45 from BirendraDh/master

New tests were added.

pyHalo/realization_extensions.py

@@ -503,8 +503,8 @@ def _get_fluctuation_halos(realization, fluctuation_amplitude, fluctuation_size,
 
     return fluctuations
 
-def corr_kappa_with_mask(kappa_map, map_size, r, mu, apply_mask=True, r_min=0,
-                         r_max=None, normalization=True):
+def corr_kappa_with_mask(kappa_map, map_size, r, mu, apply_mask=True, r_min=0.5, r_max=1.5, normalization=False): 
+
     """
     This function computes the two-point correlation function from a convergence map.
 
@@ -515,7 +515,8 @@ def corr_kappa_with_mask(kappa_map, map_size, r, mu, apply_mask=True, r_min=0,
             of the angles between 0 and 180 degrees.
     :param apply_mask: if True, apply the mask on the convergence map.
     :param r_min: inner radius of mask in units of grid coordinates
-    :param r_max: outer radius of mask in units of grid coordinates
+    :param r_max: outer radius of mask in units of grid coordinates. If r_max = None, the size of the convergence map's outside 
+            boundary becomes the mask's outer boundary.
     :param normalization: if True, apply normalization to the correlation function.
     :return: the two-point correlation function on the (mu, r) coordinate grid.
     """
@@ -550,11 +551,9 @@ def corr_kappa_with_mask(kappa_map, map_size, r, mu, apply_mask=True, r_min=0,
     if apply_mask == True:
         mask = mask_annular(center_x, center_y, XX_, YY_, r_min, r_max)
         mask_interp = RectBivariateSpline(X_, Y_, mask, kx=1, ky=1, s=0)
-
     else:
         mask = np.ones(XX_.shape)
-
-
+
     kappa_interp = RectBivariateSpline(X_, Y_, kappa_map, kx=5, ky=5, s=0)
 
     corr = np.zeros((r.shape[0], mu.shape[0]))

tests/test_realization_extensions.py

@@ -1,6 +1,6 @@
 import pytest
 from pyHalo.single_realization import SingleHalo
-from pyHalo.realization_extensions import RealizationExtensions, corr_kappa_with_mask, _xi_l, _xi_l_to_Pk_l, fit_correlation_multipole
+from pyHalo.realization_extensions import RealizationExtensions, corr_kappa_with_mask, xi_l, xi_l_to_Pk_l, fit_correlation_multipole
 from pyHalo.Cosmology.cosmology import Cosmology
 from scipy.interpolate import interp1d
 import numpy.testing as npt
@@ -417,6 +417,11 @@ def kappa_GRF(delta_pix, num_pix, alpha):
         kappa = kappa_GRF(delta_pix, npix, alpha)
 
         corr = corr_kappa_with_mask(kappa, window_size, r, mu, apply_mask = False, r_min = 0, r_max = None, normalization = False)
+        corr_mask = corr_kappa_with_mask(kappa, window_size, r, mu, apply_mask = True, r_min = 0.5, r_max = None, normalization = False)
+        corr_mask_ann = corr_kappa_with_mask(kappa, window_size, r, mu, apply_mask = True, r_min = 0.5, r_max = 1.5, normalization = False)
+        corr_norm = corr_kappa_with_mask(kappa, window_size, r, mu, apply_mask = False, r_min = 0, r_max = None, normalization = True)
+        corr_mask_norm = corr_kappa_with_mask(kappa, window_size, r, mu, apply_mask = True, r_min = 0.5, r_max = None, normalization = True)
+        corr_mask_ann_norm = corr_kappa_with_mask(kappa, window_size, r, mu, apply_mask = True, r_min = 0.5, r_max = 1.5, normalization = True)
 
         xi_0_real = delta_pix**(2-alpha)/(2*np.pi*r)
 
@@ -425,8 +430,14 @@ def kappa_GRF(delta_pix, num_pix, alpha):
         xi_l_grid = np.transpose([xi_0_real] *mu.shape[0])
 
         corr_real = xi_l_grid*T_l_grid
-
-        npt.assert_array_almost_equal(corr_real,corr, decimal=2)
+        corr_real_norm = np.linalg.norm(corr_real, 1)*corr_real
+
+        npt.assert_array_almost_equal(corr_real, corr, decimal=2)
+        npt.assert_array_almost_equal(corr_real, corr_mask, decimal=2)
+        npt.assert_array_almost_equal(corr_real, corr_mask_ann, decimal=2)
+        npt.assert_array_almost_equal(corr_real_norm,corr_norm, 1)
+        npt.assert_array_almost_equal(corr_real_norm,corr_mask_norm, 1)
+        npt.assert_array_almost_equal(corr_real_norm,corr_mask_ann_norm, 1)
 
     def test_xi_l(self):
         mu = np.linspace(-1, 1, 100)

tests/test_utilities.py

@@ -1,10 +1,12 @@
 from lenstronomy.LensModel.lens_model import LensModel
 import numpy.testing as npt
+from scipy.integrate import simps
 import numpy as np
 from pyHalo.utilities import interpolate_ray_paths, de_broglie_wavelength, delta_sigma, ITSampling, \
     inverse_transform_sampling, delta_kappa, nfw_velocity_dispersion
 from pyHalo.Cosmology.cosmology import Cosmology
 from pyHalo.Halos.lens_cosmo import LensCosmo
+from pyHalo.utilities import mask_annular
 import pytest
 
 class TestUtilities(object):
@@ -87,6 +89,29 @@ def test_ITSampling(self):
         npt.assert_almost_equal(np.mean(x_samples)/mu, 1.0, 2)
         npt.assert_almost_equal(np.std(x_samples)/sigma, 1.0, 2)
 
+    def test_mask_annular(self):
+        npix = 1000
+        window_size = 4
+        delta_pix = window_size/npix
+        _R = np.linspace(-window_size/2, window_size/2, npix)
+        XX, YY = np.meshgrid(_R, _R)
+        r_min, r_max = 0.5, 1.5
+
+        no_mask = mask_annular(0, 0, XX, YY, r_min = 0, r_max = None)
+        mask_1 = mask_annular(0, 0, XX, YY, r_min, r_max = None)
+        mask_2 = mask_annular(0, 0, XX, YY, r_min, r_max)
+
+        area_no_mask = simps(simps(no_mask, _R, axis=0), _R, axis=-1)
+        area_mask_1 = simps(simps(mask_1, _R, axis=0), _R, axis=-1)
+        area_mask_2 = simps(simps(mask_2, _R, axis=0), _R, axis=-1)
+
+        area_no_mask_real = window_size**2
+        area_mask_1_real = window_size**2 - np.pi*r_min**2
+        area_mask_2_real = np.pi*(r_max**2 - r_min**2)
+
+        npt.assert_array_almost_equal(area_no_mask_real,area_no_mask, decimal=3)
+        npt.assert_array_almost_equal(area_mask_1_real,area_mask_1, decimal=3)
+        npt.assert_array_almost_equal(area_mask_2_real,area_mask_2, decimal=3)
 
 if __name__ == '__main__':