cleanups and tests

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;
 
@@ -87,7 +87,6 @@ template<typename T> class Interpolator
           {
           if (j2>=nphi) j2-=nphi;
           tmp.v(i2,j) = sfct*tmp(i,j2);
-//tmp.v(i2,j)=0.;
           }
       r2r_fftpack(ftmp,ftmp,{0,1},true,true,1.,nthreads);
       auto fct = kernel.correction_factors(nphi, nphi0/2+1, nthreads);
@@ -98,14 +97,9 @@ template<typename T> class Interpolator
           tmp0.v(i,j) *= fct[(i+1)/2] * fct[(j+1)/2];
       // FFT to (theta, phi) domain on minimal grid
       r2r_fftpack(ftmp0,ftmp0,{1,0},false, false,1./(nphi0*nphi0),nthreads);
-arr.apply([](T &v){v=0.;});
       for (size_t i=0; i<ntheta0; ++i)
         for (size_t j=0; j<nphi0; ++j)
           arr.v(i,j) = tmp0(i,j);
-// adjoint of the extension to 2pi in theta
-//       for (size_t i=1; i+1<ntheta0; ++i)
-//         for (size_t j=0; j<nphi0; ++j)
-//           arr.v(i,j)*=2;
       }
 
   public:
@@ -142,9 +136,6 @@ arr.apply([](T &v){v=0.;});
       auto m1 = cube.template subarray<2>({supp,supp,0},{ntheta,nphi,0});
       sharp_alm2map(a1.Alms().data(), m1.vdata(), *ginfo, *ainfo, 0, nthreads);
       correct(m1,0);
-//  for (size_t i=1; i+1<ntheta; ++i)
-//    for (size_t j=0; j<nphi; ++j)
-//      m1.v(i,j)*=2;
       }
       for (size_t k=1; k<=kmax; ++k)
         {

interpol_ng/test/test_interpol_ng.py

@@ -39,7 +39,19 @@ def convolve(alm1, alm2, lmax):
     return job.map2alm(map)[0]*np.sqrt(4*np.pi)
 
 
-@pmp("lkmax", [(43,43),(2,1),(30,15),(512,2)])
+def compress_alm(alm,lmax):
+    res = np.empty(2*len(alm)-lmax-1, dtype=np.float64)
+    res[0:lmax+1] = alm[0:lmax+1].real
+    res[lmax+1::2] = np.sqrt(2)*alm[lmax+1:].real
+    res[lmax+2::2] = np.sqrt(2)*alm[lmax+1:].imag
+    return res
+
+
+def myalmdot(a1,a2,lmax,mmax,spin):
+    return np.vdot(compress_alm(a1,lmax),compress_alm(a2,lmax))
+
+
+@pmp("lkmax", [(43,43),(2,1),(30,15),(125,2)])
 def test_against_convolution(lkmax):
     lmax, kmax = lkmax
     slmT = random_alm(lmax, lmax)
@@ -62,3 +74,21 @@ def test_against_convolution(lkmax):
         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)
+
+@pmp("lkmax", [(43,43),(2,1),(30,15),(125,2)])
+def test_adjointness(lkmax):
+    lmax, kmax = lkmax
+    slmT = random_alm(lmax, lmax)
+    blmT = random_alm(lmax, kmax)
+    nptg=100000
+    ptg=np.random.uniform(0.,1.,nptg*3).reshape(nptg,3)
+    ptg[:,0]*=np.pi
+    ptg[:,1]*=2*np.pi
+    ptg[:,2]*=2*np.pi
+    foo = interpol_ng.PyInterpolator(slmT,blmT,lmax, kmax, epsilon=1e-6, nthreads=2)
+    inter1=foo.interpol(ptg)
+    fake = np.random.uniform(0.,1., ptg.shape[0])
+    foo2 = interpol_ng.PyInterpolator(lmax, kmax, epsilon=1e-6, nthreads=2)
+    foo2.deinterpol(ptg.reshape((-1,3)), fake)
+    bla=foo2.getSlm(blmT)
+    _assert_close(myalmdot(slmT, bla, lmax, lmax, 0), np.vdot(fake,inter1), 1e-12)