Fix a few bugs and investigate issue with C compilers in mac (#238)

* fixed two easy bugs

* catalog2alm parallel

* better CAR tests

Author: damonge

pymaster/tests/test_master_car.py

@@ -3,91 +3,54 @@
 import pymaster as nmt
 
 
-class WorkspaceTesterCAR(object):
+class CARTester(object):
     def __init__(self):
         from astropy.io import fits
         from astropy.wcs import WCS
 
-        # Read mask
         hdul = fits.open("test/benchmarks/msk_car.fits")
         self.msk = hdul[0].data
-        # Set up coordinates
         self.wcs = WCS(hdul[0].header)
-        self.nx, self.ny = self.wcs.pixel_shape
         hdul.close()
-        # Read maps
+
         hdul = fits.open("test/benchmarks/mps_car.fits")
-        self.mps = hdul[0].data
+        self.mp = hdul[0].data
         hdul.close()
 
-        self.minfo = nmt.NmtMapInfo(self.wcs, (self.ny, self.nx))
-        self.lmax = self.minfo.get_lmax()
-        self.nlb = 50
-        self.npix = self.minfo.npix
-        self.b = nmt.NmtBin.from_lmax_linear(self.lmax, self.nlb)
-        (self.l, self.cltt, self.clte,
-         self.clee, self.clbb, self.cltb,
-         self.cleb) = np.loadtxt("test/benchmarks/pspy_cls.txt", unpack=True)
-        self.f0 = nmt.NmtField(self.msk, [self.mps[0]], wcs=self.wcs, n_iter=0)
-        self.f0_half = nmt.NmtField(self.msk[:self.ny//2, :self.nx//2],
-                                    [self.mps[0][:self.ny//2, :self.nx//2]],
-                                    wcs=self.wcs, n_iter=0)
-        self.b_half = nmt.NmtBin.from_lmax_linear(self.lmax//2, self.nlb)
-        self.b_doub = nmt.NmtBin.from_lmax_linear(self.lmax*2, self.nlb)
-        self.n_good = np.zeros([1, (self.lmax+1)])
-        self.n_bad = np.zeros([2, (self.lmax+1)])
-        self.n_half = np.zeros([1, (self.lmax//2+1)])
-
-
-WT = WorkspaceTesterCAR()
-
-
-@pytest.mark.skipif(True, reason='slow')
-def test_workspace_car_master():
-    f0 = nmt.NmtField(WT.msk, [WT.mps[0]],
-                      wcs=WT.wcs, n_iter=0)
-    f2 = nmt.NmtField(WT.msk, [WT.mps[1], WT.mps[2]],
-                      wcs=WT.wcs, n_iter=0)
-    f = [f0, f2]
-
-    for ip1 in range(2):
-        for ip2 in range(ip1, 2):
-            if ip1 == ip2 == 0:
-                cl_bm = np.array([WT.cltt])
-            elif ip1 == ip2 == 1:
-                cl_bm = np.array([WT.clee, WT.cleb,
-                                  WT.cleb, WT.clbb])
-            else:
-                cl_bm = np.array([WT.clte, WT.cltb])
-            w = nmt.NmtWorkspace.from_fields(f[ip1], f[ip2], WT.b)
-            cl = w.decouple_cell(nmt.compute_coupled_cell(f[ip1], f[ip2]))
-            assert np.amax(np.fabs(cl-cl_bm)) <= 1E-10
-
-    # TEB
-    w = nmt.NmtWorkspace.from_fields(f[0], f[1], WT.b,
-                                     is_teb=True)
-    c00 = nmt.compute_coupled_cell(f[0], f[0])
-    c02 = nmt.compute_coupled_cell(f[0], f[1])
-    c22 = nmt.compute_coupled_cell(f[1], f[1])
-    cl = w.decouple_cell(np.array([c00[0], c02[0], c02[1],
-                                   c22[0], c22[1], c22[2],
-                                   c22[3]]))
-    cl_bm = np.array([WT.cltt, WT.clte, WT.cltb,
-                      WT.clee, WT.cleb, WT.cleb,
-                      WT.clbb])
-    assert np.amax(np.fabs(cl-cl_bm)) <= 1E-10
-
-
-@pytest.mark.skipif(False, reason='slow')
+        self.leff, self.cl_bm = np.loadtxt(
+            "test/benchmarks/bm_car.txt", unpack=True)
+
+        self.lmax = 100
+        self.f0 = nmt.NmtField(self.msk, [self.mp],
+                               wcs=self.wcs, lmax=self.lmax)
+        self.b = nmt.NmtBin.from_lmax_linear(self.lmax, nlb=5)
+
+
+WT = CARTester()
+
+
+def test_cl_car():
+    w = nmt.NmtWorkspace.from_fields(WT.f0, WT.f0, WT.b)
+    cl = w.decouple_cell(nmt.compute_coupled_cell(WT.f0,
+                                                  WT.f0))[0]
+
+    assert np.amax(np.fabs(cl-WT.cl_bm)) <= 1E-10
+
+
 def test_workspace_car_methods():
     w = nmt.NmtWorkspace.from_fields(WT.f0, WT.f0,
                                      WT.b)  # OK init
     # Incompatible bandpowers
     with pytest.raises(ValueError):
-        w.compute_coupling_matrix(WT.f0, WT.f0,
-                                  WT.b_doub)
+        bd = nmt.NmtBin.from_lmax_linear(2*WT.lmax, nlb=5)
+        w.compute_coupling_matrix(WT.f0, WT.f0, bd)
+
     # Incompatible resolutions
-    assert not WT.f0.is_compatible(WT.f0_half)
+    nx, ny = WT.wcs.pixel_shape
+    fhalf = nmt.NmtField(WT.msk[:ny//2, :nx//2],
+                         [WT.mp[:ny//2, :nx//2]],
+                         wcs=WT.wcs, lmax=WT.lmax)
+    assert not WT.f0.is_compatible(fhalf)
 
     # Wrong fields for TEB
     with pytest.raises(RuntimeError):
@@ -97,32 +60,36 @@ def test_workspace_car_methods():
     w.compute_coupling_matrix(WT.f0, WT.f0, WT.b)
 
     # Test couple_cell
-    c = w.couple_cell(WT.n_good)
+    n_good = np.ones([1, WT.lmax+1])
+    n_bad = np.ones([2, WT.lmax+1])
+    n_half = np.ones([1, WT.lmax//2+1])
+    c = w.couple_cell(n_good)
     assert c.shape == (1, w.wsp.lmax+1)
     with pytest.raises(ValueError):
-        w.couple_cell(WT.n_bad)
+        w.couple_cell(n_bad)
     with pytest.raises(ValueError):
-        w.couple_cell(WT.n_half)
+        w.couple_cell(n_half)
 
     # Test decouple_cell
-    c = w.decouple_cell(WT.n_good)
+    c = w.decouple_cell(n_good)
     assert c.shape == (1, WT.b.bin.n_bands)
     with pytest.raises(ValueError):
-        w.decouple_cell(WT.n_bad)
+        w.couple_cell(n_bad)
     with pytest.raises(ValueError):
-        w.decouple_cell(WT.n_half)
-    c = w.decouple_cell(WT.n_good, cl_bias=WT.n_good,
-                        cl_noise=WT.n_good)
+        w.couple_cell(n_half)
+    c = w.decouple_cell(n_good, cl_bias=n_good,
+                        cl_noise=n_good)
     assert c.shape == (1, WT.b.bin.n_bands)
     with pytest.raises(ValueError):
-        w.decouple_cell(WT.n_good, cl_bias=WT.n_good,
-                        cl_noise=WT.n_bad)
+        w.decouple_cell(n_good, cl_bias=n_good,
+                        cl_noise=n_bad)
     with pytest.raises(ValueError):
-        w.decouple_cell(WT.n_good, cl_bias=WT.n_good,
-                        cl_noise=WT.n_half)
+        w.decouple_cell(n_good, cl_bias=n_good,
+                        cl_noise=n_half)
     with pytest.raises(ValueError):
-        w.decouple_cell(WT.n_good, cl_bias=WT.n_bad,
-                        cl_noise=WT.n_good)
+        w.decouple_cell(n_good, cl_bias=n_bad,
+                        cl_noise=n_good)
     with pytest.raises(ValueError):
-        w.decouple_cell(WT.n_good, cl_bias=WT.n_half,
-                        cl_noise=WT.n_good)
+        w.decouple_cell(n_good, cl_bias=n_half,
+                        cl_noise=n_good)
+    return

pymaster/utils.py

@@ -849,11 +849,12 @@ def _map2alm_healpy(map, spin, sht_info, alm_info):
           'healpy': _alm2map_healpy}
 
 
-def _alm2catalog_ducc0(alms, positions, spin, lmax):
+def _alm2catalog_ducc0(alms, positions, spin, lmax, nthreads=0):
     alms = np.atleast_2d(alms)
     values = ducc0.sht.synthesis_general(alm=alms, spin=spin,
                                          lmax=lmax, loc=positions.T,
-                                         epsilon=1E-5)
+                                         epsilon=1E-5,
+                                         nthreads=nthreads)
     return values
 
 

pyproject.toml

@@ -7,24 +7,21 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "pymaster"
-version = "2.5"
+version = "2.5.1"
 authors = [
   {name="David Alonso", email="[email protected]"}
 ]
 description = "Library for pseudo-Cl computation"
 readme = {file = "README.md", content-type = "text/markdown"}
 license = {file = "LICENSE"}
-requires-python = ">=3.8"
+requires-python = ">=3.10"
 dependencies = [
 	"healpy",
 	"ducc0"
 ]
 classifiers=[
 	'License :: OSI Approved :: BSD License',
-	'Programming Language :: Python :: 2',
-	'Programming Language :: Python :: 2.7',
 	'Programming Language :: Python :: 3',
-	'Programming Language :: Python :: 3.6',
 	'Operating System :: Unix',
 	'Operating System :: MacOS'
 ]

test/benchmarks/bm_car.txt

@@ -0,0 +1,19 @@
+4.000000000000000000e+00 5.510475822947054347e-02
+9.000000000000000000e+00 4.207912493588251979e-02
+1.400000000000000000e+01 4.075652624298726351e-02
+1.900000000000000000e+01 3.532924063446543639e-02
+2.400000000000000000e+01 2.892971272739440838e-02
+2.900000000000000000e+01 2.816949340835114488e-02
+3.400000000000000000e+01 2.007970877841726634e-02
+3.900000000000000000e+01 1.992746104802381457e-02
+4.400000000000000000e+01 1.846897729604191768e-02
+4.900000000000000000e+01 1.750155971048967737e-02
+5.400000000000000000e+01 1.551095715555601134e-02
+5.900000000000000000e+01 1.640164463469220499e-02
+6.400000000000000000e+01 1.279488746149882658e-02
+6.900000000000000000e+01 1.284001677258483388e-02
+7.400000000000000000e+01 1.268205528084791721e-02
+7.900000000000000000e+01 1.205514296421096414e-02
+8.400000000000000000e+01 1.143024725294552839e-02
+8.900000000000000000e+01 1.033587911480781697e-02
+9.400000000000000000e+01 9.231566612771795188e-03

test/benchmarks/generate_benchmarks.py

@@ -92,6 +92,30 @@ def getmaskapoana_car(w,aps,fsk=0.1,dec0=-50,ra0=0.) :
 #write_flat_map("msk_car_small.fits",mask,wcs,"mask")
 
 
+
+from pixell import enmap, utils, enplot, curvedsky
+
+shape, wcs = enmap.fullsky_geometry(res=np.deg2rad(0.5))
+lmax = 100
+ls = np.arange(lmax+1)
+cl_in = 1/(ls+10)
+
+mp = curvedsky.rand_map(shape, wcs, ps=cl_in)
+b = nmt.NmtBin.from_lmax_linear(lmax, nlb=5)
+cl1 = b.bin_cell(hp.alm2cl(curvedsky.map2alm(mp, lmax=lmax)))
+
+mp.write("mps_car_b.fits")
+
+msk = enmap.ones(shape, wcs)
+msk.write("msk_car_b.fits")
+
+np.savetxt("bm_car.txt",
+           np.transpose([b.get_effective_ells(), cl1]))
+
+
+
+
+
 ##########################
 # Flat-sky stuff
 wcs,msk=read_flat_map("msk_flat.fits")