reintroduce visual comparison in demo

interpol_ng/demo.py

@@ -4,7 +4,7 @@ import pysharp
 import time
 import matplotlib.pyplot as plt
 
-np.random.seed(2099)
+np.random.seed(48)
 
 def nalm(lmax, mmax):
     return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)
@@ -17,6 +17,7 @@ def random_alm(lmax, mmax):
     res[0:lmax+1].imag = 0.
     return res
 
+
 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
@@ -24,10 +25,21 @@ def compress_alm(alm,lmax):
     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(np.conj(a2),lmax))
 
-lmax=30
+
+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)
+
+
+lmax=60
 kmax=13
 
 
@@ -37,6 +49,44 @@ slmT = random_alm(lmax, lmax)
 
 # build beam a_lm
 blmT = random_alm(lmax, kmax)
+
+
+t0=time.time()
+# build interpolator object for slmT and blmT
+foo = interpol_ng.PyInterpolator(slmT,blmT,lmax, kmax, epsilon=1e-6, nthreads=2)
+print("setup time: ",time.time()-t0)
+nth = lmax+1
+nph = 2*lmax+1
+
+
+# compute a convolved map at a fixed psi and compare it to a map convolved
+# "by hand"
+
+ptg = np.zeros((nth,nph,3))
+ptg[:,:,0] = (np.pi*(0.5+np.arange(nth))/nth).reshape((-1,1))
+ptg[:,:,1] = (2*np.pi*(0.5+np.arange(nph))/nph).reshape((1,-1))
+ptg[:,:,2] = np.pi*0.2
+t0=time.time()
+# do the actual interpolation
+bar=foo.interpol(ptg.reshape((-1,3))).reshape((nth,nph))
+print("interpolation time: ", time.time()-t0)
+plt.subplot(2,2,1)
+plt.imshow(bar.reshape((nth,nph)))
+bar2 = np.zeros((nth,nph))
+blmTfull = np.zeros(slmT.size)+0j
+blmTfull[0:blmT.size] = blmT
+for ith in range(nth):
+    rbeamth=interpol_ng.rotate_alm(blmTfull, lmax, ptg[ith,0,2],ptg[ith,0,0],0)
+    for iph in range(nph):
+        rbeam=interpol_ng.rotate_alm(rbeamth, lmax, 0, 0, ptg[ith,iph,1])
+        bar2[ith,iph] = convolve(slmT, rbeam, lmax).real
+plt.subplot(2,2,2)
+plt.imshow(bar2)
+plt.subplot(2,2,3)
+plt.imshow((bar2-bar.reshape((nth,nph))))
+plt.show()
+
+
 ptg=np.random.uniform(0.,1.,3*1000000).reshape(1000000,3)
 ptg[:,0]*=np.pi
 ptg[:,1]*=2*np.pi

interpol_ng/test_adjoint.py

-import interpol_ng
-import numpy as np
-import pypocketfft as ppf
-import matplotlib.pyplot as plt
-import pysharp
-
-np.random.seed(42)
-lmax=1000
-mmax=lmax
-kmax=10
-
-
-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 theta_extend(arr, spin):
-    nth, nph = arr.shape
-    arr2 = np.zeros(((nth-1)*2,nph))
-    arr2[0:nth,:] = arr
-    arr2[nth:,:] = np.roll(arr[nth-2:0:-1,:],nph//2,axis=1)
-    if spin&1:
-        arr2 = -arr2
-    return arr2
-
-def mydot(a1,a2,spin):
-    return np.vdot(theta_extend(a1,spin),theta_extend(a2,spin))
-
-
-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))
-spin=0
-inter = interpol_ng.PyInterpolator(lmax, kmax, epsilon=1e-6, nthreads=1)
-nphi0, ntheta0, nphi, ntheta = inter.Nphi0(), inter.Ntheta0(), inter.Nphi(), inter.Ntheta()
-
-alm0 = random_alm(lmax, mmax)
-job = pysharp.sharpjob_d()
-job.set_triangular_alm_info(lmax, lmax)
-job.set_cc_geometry(ntheta0, nphi0)
-in0 = job.alm2map(alm0).reshape((ntheta0,nphi0))
-out0=inter.test_correct(in0,spin)
-print(out0.shape,out0.strides)
-job.set_cc_geometry(ntheta, nphi)
-out0x=np.copy(out0)
-out0x[1:-1,:]*=2
-aout0=job.alm2map_adjoint(out0x.reshape((-1,)))
-
-alm1 = random_alm(lmax, mmax)
-in1 = job.alm2map(alm1).reshape((ntheta,nphi))
-out1=inter.test_decorrect(in1,spin)
-print(out1.shape,out1.strides)
-job.set_cc_geometry(ntheta0, nphi0)
-out1x=np.copy(out1)
-out1x[1:-1,:]*=2
-aout1=job.alm2map_adjoint(out1x.reshape((-1,)))
-
-print(myalmdot(alm0,aout1,lmax,lmax,spin)/myalmdot(alm1,aout0,lmax,lmax,spin))
-print(nphi0, ntheta0, nphi, ntheta)