cleanup, tests

interpol_ng/demo.py

@@ -41,8 +41,7 @@ kmax=lmax
 # the a_lm arrays follow the same conventions as those in healpy
 slmT = random_alm(lmax, mmax)
 
-# build beam a_lm (pencil beam for now)
-blmT = deltabeam(lmax,kmax)
+# build beam a_lm
 blmT = random_alm(lmax, mmax)
 
 t0=time.time()
@@ -66,7 +65,7 @@ bar2 = np.zeros((nth,nph))
 for ith in range(nth):
     for iph in range(nph):
         rbeam=interpol_ng.rotate_alm(blmT, lmax, ptg[ith,iph,2],ptg[ith,iph,0],ptg[ith,iph,1])
-        bar2[ith,iph] = convolve(slmT, rbeam, lmax)
+        bar2[ith,iph] = convolve(slmT, rbeam, lmax).real
 plt.subplot(2,2,2)
 plt.imshow(bar2)
 plt.subplot(2,2,3)

interpol_ng/interpol_ng.cc

@@ -17,7 +17,7 @@
 #include "alm.h"
 #include "mr_util/math/fft.h"
 #include "mr_util/bindings/pybind_utils.h"
-#include <iostream>
+
 using namespace std;
 using namespace mr;
 
@@ -107,12 +107,10 @@ template<typename T> class Interpolator
               a1(l,m)=a2(l,m)=0.;
             else
               {
-              a1(l,m) = slmT(l,m)*blmT(l,k).real()*T(spinsign*lnorm[l]);
+              auto tmp = blmT(l,k)*T(spinsign*lnorm[l]);
+              a1(l,m) = slmT(l,m)*tmp.real();
               if (k>0)
-                {
-                complex<T> tmp = slmT(l,m)*complex<T>(0,blmT(l,k).imag())*T(-spinsign*lnorm[l]);
-                a2(l,m) = complex<T>(-tmp.imag(), tmp.real());
-                }
+                a2(l,m) = slmT(l,m)*tmp.imag();
               }
             }
         size_t kidx1 = (k==0) ? 0 : 2*k-1,

interpol_ng/test/test_interpol_ng.py

+import numpy as np
+import pytest
+from numpy.testing import assert_
+import interpol_ng
+import pysharp
+
+pmp = pytest.mark.parametrize
+
+
+def _l2error(a, b):
+    return np.sqrt(np.sum(np.abs(a-b)**2)/np.sum(np.abs(a)**2))
+
+
+def _assert_close(a, b, epsilon):
+    err = _l2error(a, b)
+    if (err >= epsilon):
+        print("Error: {} > {}".format(err, epsilon))
+    assert_(err<epsilon)
+
+
+def nalm(lmax, mmax):
+    return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)
+
+
+def random_alm(lmax, mmax):
+    res = np.random.uniform(-1., 1., nalm(lmax, mmax)) \
+     + 1j*np.random.uniform(-1., 1., nalm(lmax, mmax))
+    # make a_lm with m==0 real-valued
+    res[0:lmax+1].imag = 0.
+    return res
+
+
+def convolve(alm1, alm2, lmax):
+    job = pysharp.sharpjob_d()
+    job.set_triangular_alm_info(lmax, lmax)
+    job.set_gauss_geometry(lmax+1, 2*lmax+1)
+    map = job.alm2map(alm1)*job.alm2map(alm2)
+    job.set_triangular_alm_info(0,0)
+    return job.map2alm(map)[0]*np.sqrt(4*np.pi)
+
+
+@pmp("lkmax", [(43,43),(2,1),(30,15),(512,2)])
+def test_against_convolution(lkmax):
+    lmax, kmax = lkmax
+    slmT = random_alm(lmax, lmax)
+    blmT = random_alm(lmax, kmax)
+
+    inter = interpol_ng.PyInterpolator(slmT, blmT, lmax, kmax, epsilon=1e-8,
+                                       nthreads=2)
+    nptg = 50
+    ptg = np.zeros((nptg,3))
+    ptg[:,0] = np.random.uniform(0, np.pi, nptg)
+    ptg[:,1] = np.random.uniform(0, 2*np.pi, nptg)
+    ptg[:,2] = np.random.uniform(-np.pi, np.pi, nptg)
+
+    res1 = inter.interpol(ptg)
+
+    res2 = np.zeros(nptg)
+    blmT2 = np.zeros(nalm(lmax,lmax))+0j
+    blmT2[0:blmT.shape[0]] = blmT
+    for i in range(nptg):
+        rbeam=interpol_ng.rotate_alm(blmT2, lmax, ptg[i,2],ptg[i,0],ptg[i,1])
+        res2[i] = convolve(slmT, rbeam, lmax).real
+    _assert_close(res1, res2, 1e-7)