Commit 1da6d79c authored by Martin Reinecke's avatar Martin Reinecke
Browse files

pimp the demo script

parent 7a81a8d0
import pyinterpol_ng
import numpy as np
import pysharp
import time
import matplotlib.pyplot as plt
np.random.seed(48)
......@@ -39,11 +37,11 @@ def convolve(alm1, alm2, lmax):
return job.map2alm(map)[0]*np.sqrt(4*np.pi)
lmax=60
lmax=1024
kmax=13
ncomp=1
separate=False
nptg = 1000000
separate=True
nptg = 50000000
epsilon = 1e-4
ofactor = 1.5
nthreads = 0 # use as many threads as available
......@@ -61,57 +59,70 @@ blm = random_alm(lmax, kmax, ncomp)
t0=time.time()
# build interpolator object for slm and blm
foo = pyinterpol_ng.PyInterpolator(slm,blm,separate,lmax, kmax, epsilon=epsilon, ofactor=ofactor, nthreads=nthreads)
print("setup time: ",time.time()-t0)
print("support:",foo.support())
t1 = time.time()-t0
print("Convolving sky and beam with lmax=mmax={}, kmax={}".format(lmax,kmax))
print("Interpolation taking place with a maximum error of {}\n"
"and an oversampling factor of {}".format(epsilon, ofactor))
supp = foo.support()
print("(resulting in a kernel support size of {}x{})".format(supp,supp))
if ncomp == 1:
print("One component")
else:
print("{} components, which are {}coadded".format(ncomp, "not " if separate else ""))
print("\nDouble precision convolution/interpolation:")
print("preparation of interpolation grid: {}s".format(t1))
t0=time.time()
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=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)
ptg[:,0]*=np.pi
ptg[:,1]*=2*np.pi
ptg[:,2]*=2*np.pi
#foo = pyinterpol_ng.PyInterpolator(slm,blm,separate,lmax, kmax, epsilon=1e-6, nthreads=2)
t0=time.time()
bar=foo.interpol(ptg)
del foo
print("interpolation time: ", time.time()-t0)
print("Interpolating {} random angle triplets: {}s".format(nptg, time.time() -t0))
t0=time.time()
fake = np.random.uniform(0.,1., (ptg.shape[0],ncomp2))
foo2 = pyinterpol_ng.PyInterpolator(lmax, kmax, ncomp2, epsilon=epsilon, ofactor=ofactor, nthreads=nthreads)
t0=time.time()
foo2.deinterpol(ptg.reshape((-1,3)), fake)
print("deinterpolation time: ", time.time()-t0)
print("Adjoint interpolation: {}s".format(time.time() -t0))
t0=time.time()
bla=foo2.getSlm(blm)
print("getSlm time: ", time.time()-t0)
del foo2
print("Computing s_lm: {}s".format(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.)
print("Adjointness error: {}".format(v1/v2-1.))
# build interpolator object for slm and blm
t0=time.time()
foo_f = pyinterpol_ng.PyInterpolator_f(slm.astype(np.complex64),blm.astype(np.complex64),separate,lmax, kmax, epsilon=epsilon, ofactor=ofactor, nthreads=nthreads)
print("\nSingle precision convolution/interpolation:")
print("preparation of interpolation grid: {}s".format(time.time()-t0))
ptgf = ptg.astype(np.float32)
del ptg
fake_f = fake.astype(np.float32)
del fake
t0=time.time()
bar_f=foo_f.interpol(ptgf)
del foo_f
print("Interpolating {} random angle triplets: {}s".format(nptg, time.time() -t0))
foo2_f = pyinterpol_ng.PyInterpolator_f(lmax, kmax, ncomp2, epsilon=epsilon, ofactor=ofactor, nthreads=nthreads)
t0=time.time()
foo2_f.deinterpol(ptgf.reshape((-1,3)), fake_f)
print("Adjoint interpolation: {}s".format(time.time() -t0))
t0=time.time()
bla_f=foo2_f.getSlm(blm.astype(np.complex64))
del foo2_f
print("Computing s_lm: {}s".format(time.time() -t0))
v1 = np.sum([myalmdot(slm[:,i], bla_f[:,i] , lmax, lmax, 0) for i in range(ncomp)])
v2 = np.sum([np.vdot(fake_f[:,i],bar_f[:,i]) for i in range(ncomp2)])
print("Adjointness error: {}".format(v1/v2-1.))
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment