From 0f2d96fa98bde1b0e5fb03ca15f555ec715ad480 Mon Sep 17 00:00:00 2001 From: Kevin Wolz Date: Tue, 4 Jun 2024 13:32:43 +0100 Subject: [PATCH] add unit tests --- pymaster/field.py | 2 +- pymaster/tests/test_field_catalog.py | 230 +++++++++++++++++++++++++++ pymaster/tests/test_master.py | 26 +++ 3 files changed, 257 insertions(+), 1 deletion(-) create mode 100644 pymaster/tests/test_field_catalog.py diff --git a/pymaster/field.py b/pymaster/field.py index 76d4a186..9091242d 100644 --- a/pymaster/field.py +++ b/pymaster/field.py @@ -757,7 +757,7 @@ def __init__(self, positions, weights, field, lmax, lmax_mask=-1, self.field = np.array(field, dtype=np.float64) # Sanity checks if spin is None: - spin = 0 if field.ndim == 1 else 2 + spin = 0 if self.field.ndim == 1 else 2 if spin and np.shape(self.field) != (2, len(weights)): raise ValueError("Field has wrong shape.") if spin == 0 and (self.field.ndim != 1 diff --git a/pymaster/tests/test_field_catalog.py b/pymaster/tests/test_field_catalog.py new file mode 100644 index 00000000..b39b9091 --- /dev/null +++ b/pymaster/tests/test_field_catalog.py @@ -0,0 +1,230 @@ +import numpy as np +import healpy as hp +import pymaster as nmt +import warnings +import sys + + +class FieldTesterCatalog(object): + def __init__(self): + # This is to avoid showing an ugly warning that + # has nothing to do with pymaster + if (sys.version_info > (3, 1)): + warnings.simplefilter("ignore", ResourceWarning) + nside = 64 + self.ncat = hp.nside2npix(nside) + self.f = np.zeros([3, self.ncat]) + self.w = np.ones(self.ncat) + self.lmax = 2*nside - 1 + + th, ph = hp.pix2ang(nside, np.arange(self.ncat)) + self.p = np.array([th, ph]) + sth = np.sin(th) + cth = np.cos(th) + # Re(Y_22) + self.f[0] = np.sqrt(15./2./np.pi)*sth**2*np.cos(2*ph) + # _2Y^E_20 + _2Y^B_30 + self.f[1] = -np.sqrt(15./2./np.pi)*sth**2/4. + self.f[2] = -np.sqrt(105./2./np.pi)*cth*sth**2/2. + + +FT = FieldTesterCatalog() + + +def test_field_catalog_compatibility(): + # Different field values + f_rand = np.random.rand(FT.ncat)*0.1 + 1 + f0 = nmt.NmtFieldCatalog(FT.p, FT.w, FT.f[0], FT.lmax) + f1 = nmt.NmtFieldCatalog(FT.p, FT.w, f_rand, FT.lmax) + assert f0.is_compatible(f1) + + # Different positions + p1 = np.array([np.random.permutation(FT.p[0]), FT.p[1]]) + f1 = nmt.NmtFieldCatalog(p1, FT.w, f_rand, FT.lmax) + assert f0.is_compatible(f1) + + # Different weights + w_rand = np.random.rand(FT.ncat)*0.1 + 1 + f1 = nmt.NmtFieldCatalog(FT.p, w_rand, f_rand, FT.lmax) + assert f0.is_compatible(f1) + + # Different lmax + print(f0.ainfo_mask.lmax) + f1 = nmt.NmtFieldCatalog(FT.p, w_rand, f_rand, lmax=111) + assert not f0.is_compatible(f1) + + # Different lmax_mask + f1 = nmt.NmtFieldCatalog(FT.p, w_rand, f_rand, FT.lmax, lmax_mask=300) + assert not f0.is_compatible(f1) + + +def test_field_catalog_init(): + # Checks correct initialization of positions (lon/lat and theta/phi), + # weights, fields, different spins. + Ncat = 100 + lmax = 10 + np.random.seed(5675) + w0 = np.zeros(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) + lon_rad = 2*np.pi*np.random.rand(Ncat) + lon_deg = 180./np.pi*lon_rad + lat_deg = 90. - 180./np.pi*col_rad + + for ndim, vals in zip([1, 2], [val_s0, val_s2]): + f = nmt.NmtFieldCatalog([col_rad, lon_rad], w0, vals, + lmax, field_is_weighted=True) + assert np.array_equal(f.field, vals) + assert ndim == f.field.ndim + f = nmt.NmtFieldCatalog([col_rad, lon_rad], w0, vals, + lmax, field_is_weighted=False) + assert not np.any(f.field) + assert ndim == f.field.ndim + + for vals in [val_s0, val_s2]: + f1 = nmt.NmtFieldCatalog([col_rad, lon_rad], w0, vals, + lmax, lonlat=False) + f2 = nmt.NmtFieldCatalog([lon_deg, lat_deg], w0, vals, + lmax, lonlat=True) + assert np.array_equal(f1.field, f2.field) + + +def test_field_catalog_Nw(): + # Check if Nw matches mask pcl in noise-dominated regime. + Ncat = 100 + lmax = 100 + lmin = 50 + Nsims = 10 + w1 = np.ones(Ncat) + + Nw = [] + Nw_true = Ncat/(4.*np.pi) + + pcl = [] + 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) + 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()) + ell = np.arange(len(pcl_mask)) + msk = ell > lmin + pcl.append(pcl_mask[msk]) + assert f.Nw == Nw_true + Nw.append(f.Nw) + + y = np.mean(np.array(pcl) - np.array(Nw)[:, None], axis=0) + yerr = np.std(np.array(pcl) - np.array(Nw)[:, None], axis=0)/np.sqrt(Nsims) + + assert all(np.logical_and(v < 5, v > -5) for v in y/yerr) + + +def test_field_catalog_Nf(): + # Check if Nf is mathematically correct and only subtracted + # for equal fields. + Ncat = 100 + lmax = 100 + + 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) + lon_rad = 2*np.pi*np.random.rand(Ncat) + + for val in [val_s0, val_s2]: + nd = val.ndim + idx = -1 if nd == 1 else (0, -1) + + # Perturb val to extract Nf correction + val2 = np.copy(val) + val2[idx] *= (1. + 1.e-6) + f1 = nmt.NmtFieldCatalog([col_rad, lon_rad], np.ones(Ncat), val, lmax) + f2 = nmt.NmtFieldCatalog([col_rad, lon_rad], np.ones(Ncat), val2, lmax) + + cl1 = nmt.compute_coupled_cell(f1, f1) + cl2 = nmt.compute_coupled_cell(f1, f2) + nl = np.shape(cl1)[-1] + diff_true = np.sum(f1.field**2)/(4*np.pi*nd)*np.eye(nd).reshape(nd**2) + diff_true = diff_true[:, None]*np.ones((nd**2, nl)) + + assert not np.any(np.fabs(cl2 - cl1 - diff_true) > 1.e-6) + + +def test_field_catalog_alm(): + # Compare catalog-based alms with analytical input and make sure they + # are equal up to numerical accuracy. + nside = 32 + npix = int(hp.nside2npix(nside)) + pixel_area = np.pi*4./npix + lmax = 20 + + msk = np.ones(npix) + mps = np.zeros([3, npix]) + th, ph = hp.pix2ang(nside, np.arange(npix)) + sth = np.sin(th) + cth = np.cos(th) + # Re(Y_22) + mps[0] = np.sqrt(15./2./np.pi)*sth**2*np.cos(2*ph) + # _2Y^E_20 + _2Y^B_30 + mps[1] = -np.sqrt(15./2./np.pi)*sth**2/4. + mps[2] = -np.sqrt(105./2./np.pi)*cth*sth**2/2. + + # spin 0 + f0_cat = nmt.NmtFieldCatalog([th, ph], np.ones(npix), mps[0], lmax) + f0_map = nmt.NmtField(msk, [mps[0]], lmax=lmax) + + # spin 2 + f2_cat = nmt.NmtFieldCatalog([th, ph], np.ones(npix), + np.array([mps[1], mps[2]]), lmax) + f2_map = nmt.NmtField(msk, [mps[1], mps[2]], lmax=lmax) + + alms_map = np.array([f0_map.get_alms()[0], + f2_map.get_alms()[0], + f2_map.get_alms()[1]]) + 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_in = np.zeros_like(alms_map) + 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.fabs(np.real(alms_cat - alms_in)[f_idx]) < 1.e-3) + assert np.all(np.fabs(np.imag(alms_cat - alms_in)[f_idx]) < 1.e-3) + + +def test_field_catalog_errors(): + import pytest + + with pytest.raises(ValueError): # Incorrect field size + nmt.NmtFieldCatalog(np.zeros((2, 100)), np.zeros(99), np.zeros(99), 10) + with pytest.raises(ValueError): # Passing 3 fields + nmt.NmtFieldCatalog([[0., 0.], [1., 1.]], [1., 1.], + [[1., 1.], [1., 1.], [1., 1.]], 10) + with pytest.raises(ValueError): # Passing single angle + nmt.NmtFieldCatalog([[0., 0.]], [1., 1.], [1., 1.], 10) + with pytest.raises(ValueError): # Trash angles (th, phi) + nmt.NmtFieldCatalog([[-1., 0.], [1., 1.]], [1., 1.], [1., 1.], 10) + with pytest.raises(ValueError): # Trash angles (lon, lat) + nmt.NmtFieldCatalog([[-45., 0.], [30., 120.]], [1., 1.], [1., 1.], 10, + lonlat=True) + with pytest.raises(ValueError): # Passing crap beam + nmt.NmtFieldCatalog([[0., 0.], [1., 1.]], [1., 1.], [1., 1.], 10, + beam=1) + # Automatically assign spin = 0 for a single field + f = nmt.NmtFieldCatalog([[0., 0.], [1., 1.]], [1., 1.], [1., 1.], 10) + assert (f.spin == 0) + # Automatically assign spin = 2 for 2 fields + f = nmt.NmtFieldCatalog([[0., 0.], [1., 1.]], [1., 1.], + [[1., 1.], [1., 1.]], 10) + assert (f.spin == 2) + with pytest.raises(ValueError): # Spin = 0 but 2 maps + f = nmt.NmtFieldCatalog([[0., 0.], [1., 1.]], [1., 1.], + [[1., 1.], [1., 1.]], 10, spin=0) + with pytest.raises(ValueError): # Spin = 2 but single map + f = nmt.NmtFieldCatalog([[0., 0.], [1., 1.]], [1., 1.], + [1., 1.], 10, spin=2) diff --git a/pymaster/tests/test_master.py b/pymaster/tests/test_master.py index 13710510..abca54e8 100644 --- a/pymaster/tests/test_master.py +++ b/pymaster/tests/test_master.py @@ -700,3 +700,29 @@ def test_unbinned_mcm_io(): w2.get_bandpower_windows() os.system("rm test/wspc.fits") + + +def test_ducc_catalog2alm(): + from utils import _catalog2alm_ducc0 as cat2alm + nside = 128 + npix = int(hp.nside2npix(nside)) + pixel_area = np.pi*4./npix + lmax = 128 + + # input alms + alm_in = np.zeros((3, hp.Alm.getidx(lmax, lmax, lmax)+1), dtype="complex") + alm_in[0, hp.Alm.getidx(lmax, 2, 2)] = 2. + alm_in[1, hp.Alm.getidx(lmax, 2, 0)] = 1. + alm_in[2, hp.Alm.getidx(lmax, 3, 1)] = 2. + + map = hp.alm2map(alm_in, nside, lmax=lmax, mmax=lmax) + th, ph = hp.pix2ang(nside, np.arange(npix)) + + # spin 0 + alm_cat_0 = pixel_area*cat2alm(map[0], np.array([th, ph]), 0, lmax) + + # spin 2 + alm_cat_2 = pixel_area*cat2alm(map[1:3], np.array([th, ph]), 2, lmax) + + assert np.all(np.absolute(alm_cat_0[0] - alm_in[0]) < 1E-4) + assert np.all(np.absolute(alm_cat_2[:2] - alm_in[1:3]) < 1E-4)