From 34d985b7eecbdd71ef28cc2a6c92f4b45ebedeb4 Mon Sep 17 00:00:00 2001 From: damonge Date: Tue, 11 Jun 2024 18:41:26 +0100 Subject: [PATCH] test for clustering --- pymaster/tests/test_field_catalog.py | 78 +++++++++++++++++++++++----- 1 file changed, 66 insertions(+), 12 deletions(-) diff --git a/pymaster/tests/test_field_catalog.py b/pymaster/tests/test_field_catalog.py index e3cc6433..456e8ef6 100644 --- a/pymaster/tests/test_field_catalog.py +++ b/pymaster/tests/test_field_catalog.py @@ -28,7 +28,6 @@ def test_field_catalog_compatibility(): assert f0.is_compatible(f1) # Different lmax - print(f0.ainfo_mask.lmax) f1 = nmt.NmtFieldCatalog(p, w_rand, f_rand, lmax=111) assert not f0.is_compatible(f1) @@ -46,7 +45,7 @@ def test_field_catalog_init(): w1 = np.random.rand(Ncat) val_s0 = np.random.rand(Ncat) - 0.5 val_s2 = np.random.rand(2, Ncat) - 0.5 - col_rad = np.pi*np.random.rand(Ncat) + col_rad = np.arccos(-1+2*np.random.rand(Ncat)) lon_rad = 2*np.pi*np.random.rand(Ncat) lon_deg = np.rad2deg(lon_rad) lat_deg = np.rad2deg(np.pi/2. - col_rad) @@ -85,7 +84,7 @@ def test_field_catalog_Nw(): for i in range(Nsims): np.random.seed(5675 + i) val_s0 = np.random.rand(Ncat) - 0.5 - col_rad = np.pi*np.random.rand(Ncat) + col_rad = np.arccos(-1+2*np.random.rand(Ncat)) lon_rad = 2*np.pi*np.random.rand(Ncat) f = nmt.NmtFieldCatalog([col_rad, lon_rad], w1, val_s0, lmax) pcl_mask = hp.alm2cl(f.get_mask_alms()) @@ -110,7 +109,7 @@ def test_field_catalog_Nf(): np.random.seed(5675) val_s0 = np.random.rand(Ncat) - 0.5 val_s2 = np.random.rand(2, Ncat) - 0.5 - col_rad = np.pi*np.random.rand(Ncat) + col_rad = np.arccos(-1+2*np.random.rand(Ncat)) lon_rad = 2*np.pi*np.random.rand(Ncat) for val in [val_s0, val_s2]: @@ -135,7 +134,7 @@ def test_field_catalog_alm(): # are equal up to numerical accuracy. nside = 32 npix = int(hp.nside2npix(nside)) - pixel_area = np.pi*4./npix + pixel_area = hp.nside2pixarea(nside) lmax = 20 mps = np.zeros([3, npix]) @@ -152,20 +151,17 @@ def test_field_catalog_alm(): f0_cat = nmt.NmtFieldCatalog([th, ph], np.ones(npix), mps[0], lmax) # spin 2 - f2_cat = nmt.NmtFieldCatalog([th, ph], np.ones(npix), - np.array([mps[1], mps[2]]), lmax) + f2_cat = nmt.NmtFieldCatalog([th, ph], np.ones(npix), mps[1:], lmax) - alms_cat = np.array([pixel_area*f0_cat.get_alms()[0], - pixel_area*f2_cat.get_alms()[0], - pixel_area*f2_cat.get_alms()[1]]) + alms_cat = pixel_area * np.array(list(f0_cat.get_alms()) + + list(f2_cat.get_alms())) alms_in = np.zeros_like(alms_cat) alms_in[0, hp.Alm.getidx(lmax, 2, 2)] = 2. alms_in[1, hp.Alm.getidx(lmax, 2, 0)] = 1. alms_in[2, hp.Alm.getidx(lmax, 3, 0)] = 2. - for f_idx in range(3): - assert np.all(np.absolute(alms_cat - alms_in)[f_idx] < 1.e-3) + assert np.all(np.absolute(alms_cat - alms_in) < 1.e-3) def test_field_catalog_errors(): @@ -230,3 +226,61 @@ def test_field_catalog_errors(): nmt.NmtFieldCatalog([[0., 0.], [1., 1.]], [1., 1.], [[1., 1.], [1., 1.]], 10) nmt.utils.HAVE_DUCC = True + + +def test_field_catalog_clustering_poisson(): + # Checks that a purely Poisson catalog has a power spectrum + # that is close to zero. + + # Create random catalog + nside = 64 + npix = hp.nside2npix(nside) + pos_ran = np.array([np.arccos(-1+2*np.random.rand(4*npix)), + 2*np.pi*np.random.rand(4*npix)]) + w_ran = np.ones(4*npix) + # We'll do a full-sky version and a half-sky version + goodhalfran = pos_ran[0] <= np.pi/2 + + # Create dummy field for workspaces + b = nmt.NmtBin.from_nside_linear(nside, 4) + f = nmt.NmtFieldCatalogClustering(pos_ran, w_ran, pos_ran, w_ran, + lmax=3*nside-1) + w = nmt.NmtWorkspace() + w.compute_coupling_matrix(f, f, b) + # Half-sky version + f = nmt.NmtFieldCatalogClustering(pos_ran[:, goodhalfran], + w_ran[goodhalfran], + pos_ran[:, goodhalfran], + w_ran[goodhalfran], + lmax=3*nside-1) + wh = nmt.NmtWorkspace() + wh.compute_coupling_matrix(f, f, b) + + nsims = 10 + ncat = npix//4 + cls_full = [] + cls_half = [] + w_dat = np.ones(ncat) + for i in range(nsims): + # Generate sim and get C_ell + pos_dat = np.array([np.arccos(-1+2*np.random.rand(ncat)), + 2*np.pi*np.random.rand(ncat)]) + f = nmt.NmtFieldCatalogClustering(pos_dat, w_dat, pos_ran, w_ran, + lmax=3*nside-1) + cls_full.append(w.decouple_cell(nmt.compute_coupled_cell(f, f))) + + # Half-sky version + goodhalfdat = pos_dat[0] <= np.pi/2 + f = nmt.NmtFieldCatalogClustering(pos_dat[:, goodhalfdat], + w_dat[goodhalfdat], + pos_ran[:, goodhalfran], + w_ran[goodhalfran], + lmax=3*nside-1) + cls_half.append(w.decouple_cell(nmt.compute_coupled_cell(f, f))) + cls_full = np.array(cls_full).squeeze() + cls_half = np.array(cls_half).squeeze() + + for c in [cls_full, cls_half]: + x = np.mean(c, axis=0) + s = np.std(c, axis=0)/np.sqrt(nsims) + assert np.all(np.fabs(x) < 5*s)