demo.py 2.94 KB
 Martin Reinecke committed Feb 29, 2020 1 2 3 4 ``````import interpol_ng import numpy as np import pysharp import time `````` Martin Reinecke committed Mar 02, 2020 5 ``````import matplotlib.pyplot as plt `````` Martin Reinecke committed Feb 29, 2020 6 `````` `````` Martin Reinecke committed Mar 08, 2020 7 ``````np.random.seed(20) `````` Martin Reinecke committed Feb 29, 2020 8 9 10 11 `````` def nalm(lmax, mmax): return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax) `````` Martin Reinecke committed Feb 29, 2020 12 `````` `````` Martin Reinecke committed Mar 02, 2020 13 14 15 16 17 18 19 ``````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 `````` Martin Reinecke committed Mar 09, 2020 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 ``````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)) 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)) `````` Martin Reinecke committed Mar 02, 2020 41 `````` `````` Martin Reinecke committed Feb 29, 2020 42 ``````def deltabeam(lmax,kmax): `````` Martin Reinecke committed Feb 29, 2020 43 `````` beam=np.zeros(nalm(lmax, kmax))+0j `````` Martin Reinecke committed Feb 29, 2020 44 45 46 47 `````` for l in range(lmax+1): beam[l] = np.sqrt((2*l+1.)/(4*np.pi)) return beam `````` Martin Reinecke committed Mar 03, 2020 48 49 50 51 52 53 54 ``````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) `````` Martin Reinecke committed Feb 29, 2020 55 `````` `````` Martin Reinecke committed Mar 08, 2020 56 ``````lmax=40 `````` Martin Reinecke committed Mar 02, 2020 57 ``````mmax=lmax `````` Martin Reinecke committed Mar 08, 2020 58 ``````kmax=10 `````` Martin Reinecke committed Mar 02, 2020 59 60 61 62 `````` # get random sky a_lm # the a_lm arrays follow the same conventions as those in healpy `````` Martin Reinecke committed Mar 02, 2020 63 ``````slmT = random_alm(lmax, mmax) `````` Martin Reinecke committed Feb 29, 2020 64 `````` `````` Martin Reinecke committed Mar 03, 2020 65 ``````# build beam a_lm `````` Martin Reinecke committed Mar 08, 2020 66 ``````blmT = random_alm(lmax, kmax) `````` Martin Reinecke committed Feb 29, 2020 67 68 `````` t0=time.time() `````` Martin Reinecke committed Mar 02, 2020 69 ``````# build interpolator object for slmT and blmT `````` Martin Reinecke committed Mar 08, 2020 70 ``````foo = interpol_ng.PyInterpolator(slmT,blmT,lmax, kmax, epsilon=1e-6, nthreads=1) `````` Martin Reinecke committed Mar 02, 2020 71 ``````print("setup time: ",time.time()-t0) `````` Martin Reinecke committed Mar 03, 2020 72 ``````nth = 2*lmax+1 `````` Martin Reinecke committed Feb 29, 2020 73 ``````nph = 2*mmax+1 `````` Martin Reinecke committed Mar 08, 2020 74 ``````#exit() `````` Martin Reinecke committed Feb 29, 2020 75 ``````ptg = np.zeros((nth,nph,3)) `````` Martin Reinecke committed Mar 03, 2020 76 ``````ptg[:,:,0] = (np.pi*np.arange(nth)/(nth-1)).reshape((-1,1)) `````` Martin Reinecke committed Mar 02, 2020 77 ``````ptg[:,:,1] = (2*np.pi*np.arange(nph)/nph).reshape((1,-1)) `````` Martin Reinecke committed Mar 03, 2020 78 ``````ptg[:,:,2] = np.pi*0.2 `````` Martin Reinecke committed Feb 29, 2020 79 ``````t0=time.time() `````` Martin Reinecke committed Mar 02, 2020 80 ``````# do the actual interpolation `````` Martin Reinecke committed Feb 29, 2020 81 ``````bar=foo.interpol(ptg.reshape((-1,3))).reshape((nth,nph)) `````` Martin Reinecke committed Mar 02, 2020 82 ``````print("interpolation time: ", time.time()-t0) `````` Martin Reinecke committed Mar 03, 2020 83 84 85 ``````plt.subplot(2,2,1) plt.imshow(bar.reshape((nth,nph))) bar2 = np.zeros((nth,nph)) `````` Martin Reinecke committed Mar 08, 2020 86 87 ``````blmTfull = np.zeros(slmT.size)+0j blmTfull[0:blmT.size] = blmT `````` Martin Reinecke committed Mar 03, 2020 88 ``````for ith in range(nth): `````` Martin Reinecke committed Mar 09, 2020 89 `````` rbeamth=interpol_ng.rotate_alm(blmTfull, lmax, ptg[ith,0,2],ptg[ith,0,0],0) `````` Martin Reinecke committed Mar 03, 2020 90 `````` for iph in range(nph): `````` Martin Reinecke committed Mar 09, 2020 91 `````` rbeam=interpol_ng.rotate_alm(rbeamth, lmax, 0, 0, ptg[ith,iph,1]) `````` Martin Reinecke committed Mar 03, 2020 92 `````` bar2[ith,iph] = convolve(slmT, rbeam, lmax).real `````` Martin Reinecke committed Mar 03, 2020 93 94 95 96 ``````plt.subplot(2,2,2) plt.imshow(bar2) plt.subplot(2,2,3) plt.imshow((bar2-bar.reshape((nth,nph)))) `````` Martin Reinecke committed Feb 29, 2020 97 ``````plt.show() `````` Martin Reinecke committed Mar 08, 2020 98 99 100 101 102 `````` fake = np.random.uniform(-1.,1., bar.size) foo2 = interpol_ng.PyInterpolator(lmax, kmax, epsilon=1e-6, nthreads=2) foo2.deinterpol(ptg.reshape((-1,3)), fake) bla=foo2.getSlm(blmT) `````` Martin Reinecke committed Mar 09, 2020 103 ``````print(myalmdot(slmT, bla, lmax, lmax, 0)) `````` Martin Reinecke committed Mar 08, 2020 104 ``print(np.vdot(fake,bar))``