totalconvolve_accuracy.py 2.46 KB
 Martin Reinecke committed Jun 05, 2020 1 2 3 4 5 6 7 ``````import ducc_0_1.totalconvolve as totalconvolve import numpy as np import ducc_0_1.sht as sht import ducc_0_1.misc as misc import time import matplotlib.pyplot as plt `````` Martin Reinecke committed Jun 08, 2020 8 ``````rng = np.random.default_rng(42) `````` Martin Reinecke committed Jun 05, 2020 9 10 11 12 13 14 `````` def nalm(lmax, mmax): return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax) def random_alm(lmax, mmax, ncomp): `````` Martin Reinecke committed Jun 08, 2020 15 16 `````` res = rng.uniform(-1., 1., (nalm(lmax, mmax), ncomp)) \ + 1j*rng.uniform(-1., 1., (nalm(lmax, mmax), ncomp)) `````` Martin Reinecke committed Jun 05, 2020 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 `````` # make a_lm with m==0 real-valued 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 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(np.conj(a2),lmax)) def convolve(alm1, alm2, lmax): job = sht.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 ncomp=1 separate=True ncomp2 = ncomp if separate else 1 # get random sky a_lm # the a_lm arrays follow the same conventions as those in healpy slm = random_alm(lmax, lmax, ncomp) # build beam a_lm blm = random_alm(lmax, kmax, ncomp) t0=time.time() # build interpolator object for slm and blm foo = totalconvolve.PyInterpolator(slm,blm,separate,lmax, kmax, epsilon=1e-4, 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,ncomp2)) print("interpolation time: ", time.time()-t0) plt.subplot(2,2,1) plt.imshow(bar[:,:,0]) bar2 = np.zeros((nth,nph)) blmfull = np.zeros(slm.shape)+0j blmfull[0:blm.shape[0],:] = blm for ith in range(nth): rbeamth=misc.rotate_alm(blmfull[:,0], lmax, ptg[ith,0,2],ptg[ith,0,0],0) for iph in range(nph): rbeam=misc.rotate_alm(rbeamth, lmax, 0, 0, ptg[ith,iph,1]) bar2[ith,iph] = convolve(slm[:,0], rbeam, lmax).real plt.subplot(2,2,2) plt.imshow(bar2) plt.subplot(2,2,3) plt.imshow(bar2-bar[:,:,0]) plt.show()``````