demo.py 3.43 KB
 Martin Reinecke committed Mar 11, 2020 1 ``````import pyinterpol_ng `````` Martin Reinecke committed Feb 29, 2020 2 3 4 ``````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 10, 2020 7 ``````np.random.seed(48) `````` 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 Apr 11, 2020 13 14 15 ``````def random_alm(lmax, mmax, ncomp): res = np.random.uniform(-1., 1., (nalm(lmax, mmax), ncomp)) \ + 1j*np.random.uniform(-1., 1., (nalm(lmax, mmax), ncomp)) `````` Martin Reinecke committed Mar 02, 2020 16 `````` # make a_lm with m==0 real-valued `````` Martin Reinecke committed Apr 11, 2020 17 `````` res[0:lmax+1,:].imag = 0. `````` Martin Reinecke committed Mar 02, 2020 18 19 `````` return res `````` Martin Reinecke committed Mar 10, 2020 20 `````` `````` Martin Reinecke committed Mar 09, 2020 21 22 23 24 25 26 27 ``````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 `````` Martin Reinecke committed Mar 10, 2020 28 `````` `````` Martin Reinecke committed Mar 09, 2020 29 ``````def myalmdot(a1,a2,lmax,mmax,spin): `````` Martin Reinecke committed Mar 10, 2020 30 `````` return np.vdot(compress_alm(a1,lmax),compress_alm(np.conj(a2),lmax)) `````` Martin Reinecke committed Mar 02, 2020 31 `````` `````` Martin Reinecke committed Mar 10, 2020 32 33 34 35 36 37 38 39 40 41 `````` 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 Apr 24, 2020 42 43 ``````lmax=2048 kmax=8 `````` Martin Reinecke committed Apr 11, 2020 44 ``````ncomp=1 `````` Martin Reinecke committed Apr 24, 2020 45 46 ``````separate=False nptg = 10000000 `````` Martin Reinecke committed Apr 11, 2020 47 ``````ncomp2 = ncomp if separate else 1 `````` Martin Reinecke committed Apr 22, 2020 48 ``````epsilon = 1e-4 `````` Martin Reinecke committed Apr 24, 2020 49 50 51 52 ``````ofactor = 1.5 nthreads = 0 # use as many threads as available ncomp2 = ncomp if separate else 1 `````` Martin Reinecke committed Mar 02, 2020 53 54 55 `````` # get random sky a_lm # the a_lm arrays follow the same conventions as those in healpy `````` Martin Reinecke committed Apr 11, 2020 56 ``````slm = random_alm(lmax, lmax, ncomp) `````` Martin Reinecke committed Feb 29, 2020 57 `````` `````` Martin Reinecke committed Mar 03, 2020 58 ``````# build beam a_lm `````` Martin Reinecke committed Apr 11, 2020 59 ``````blm = random_alm(lmax, kmax, ncomp) `````` Martin Reinecke committed Mar 10, 2020 60 61 62 `````` t0=time.time() `````` Martin Reinecke committed Apr 14, 2020 63 ``````# build interpolator object for slm and blm `````` Martin Reinecke committed Apr 22, 2020 64 ``````foo = pyinterpol_ng.PyInterpolator(slm,blm,separate,lmax, kmax, epsilon=epsilon, ofactor=ofactor, nthreads=nthreads) `````` Martin Reinecke committed Mar 10, 2020 65 ``````print("setup time: ",time.time()-t0) `````` Martin Reinecke committed Apr 22, 2020 66 ``````print("support:",foo.support()) `````` Martin Reinecke committed Mar 10, 2020 67 68 69 70 71 72 73 ``````nth = lmax+1 nph = 2*lmax+1 # compute a convolved map at a fixed psi and compare it to a map convolved # "by hand" `````` Martin Reinecke committed Apr 24, 2020 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 ``````# 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=pyinterpol_ng.rotate_alm(blmfull[:,0], lmax, ptg[ith,0,2],ptg[ith,0,0],0) # for iph in range(nph): # rbeam=pyinterpol_ng.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() ptg=np.random.uniform(0.,1.,3*nptg).reshape(nptg,3) `````` Martin Reinecke committed Mar 09, 2020 100 101 102 ``````ptg[:,0]*=np.pi ptg[:,1]*=2*np.pi ptg[:,2]*=2*np.pi `````` Martin Reinecke committed Apr 11, 2020 103 104 ``````#foo = pyinterpol_ng.PyInterpolator(slm,blm,separate,lmax, kmax, epsilon=1e-6, nthreads=2) t0=time.time() `````` Martin Reinecke committed Mar 09, 2020 105 ``````bar=foo.interpol(ptg) `````` Martin Reinecke committed Apr 24, 2020 106 ``````del foo `````` Martin Reinecke committed Apr 11, 2020 107 108 ``````print("interpolation time: ", time.time()-t0) fake = np.random.uniform(0.,1., (ptg.shape[0],ncomp2)) `````` Martin Reinecke committed Apr 22, 2020 109 ``````foo2 = pyinterpol_ng.PyInterpolator(lmax, kmax, ncomp2, epsilon=epsilon, ofactor=ofactor, nthreads=nthreads) `````` Martin Reinecke committed Apr 11, 2020 110 111 112 113 114 115 116 117 118 ``````t0=time.time() foo2.deinterpol(ptg.reshape((-1,3)), fake) print("deinterpolation time: ", time.time()-t0) t0=time.time() bla=foo2.getSlm(blm) print("getSlm time: ", time.time()-t0) v1 = np.sum([myalmdot(slm[:,i], bla[:,i] , lmax, lmax, 0) for i in range(ncomp)]) v2 = np.sum([np.vdot(fake[:,i],bar[:,i]) for i in range(ncomp2)]) print(v1/v2-1.)``````