Commit 4607010f authored by Martin Reinecke's avatar Martin Reinecke
Browse files

adjointness looks really good except for theta borders

parent bb0aff75
Pipeline #70469 passed with stages
in 8 minutes and 58 seconds
......@@ -4,7 +4,7 @@ import pysharp
import time
import matplotlib.pyplot as plt
np.random.seed(20)
np.random.seed(20909)
def nalm(lmax, mmax):
return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)
......@@ -27,35 +27,12 @@ def compress_alm(alm,lmax):
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))
def deltabeam(lmax,kmax):
beam=np.zeros(nalm(lmax, kmax))+0j
for l in range(lmax+1):
beam[l] = np.sqrt((2*l+1.)/(4*np.pi))
return beam
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=40
lmax=512
mmax=lmax
kmax=10
kmax=12
# get random sky a_lm
......@@ -64,41 +41,14 @@ slmT = random_alm(lmax, mmax)
# build beam a_lm
blmT = random_alm(lmax, kmax)
t0=time.time()
# build interpolator object for slmT and blmT
ptg = np.array([[1.2,0.01,1.],[1.4,0.7,2.]])
foo = interpol_ng.PyInterpolator(slmT,blmT,lmax, kmax, epsilon=1e-6, nthreads=1)
print("setup time: ",time.time()-t0)
nth = 2*lmax+1
nph = 2*mmax+1
#exit()
ptg = np.zeros((nth,nph,3))
ptg[:,:,0] = (np.pi*np.arange(nth)/(nth-1)).reshape((-1,1))
ptg[:,:,1] = (2*np.pi*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()
fake = np.random.uniform(-1.,1., bar.size)
bar=foo.interpol(ptg)
print(foo.Nphi(),foo.Nphi0())
fake = np.random.uniform(0.,1., ptg.shape[0])
foo2 = interpol_ng.PyInterpolator(lmax, kmax, epsilon=1e-6, nthreads=2)
foo2.deinterpol(ptg.reshape((-1,3)), fake)
bla=foo2.getSlm(blmT)
print(myalmdot(slmT, bla, lmax, lmax, 0))
print(myalmdot(slmT, np.conj(bla), lmax, lmax, 0))
print(np.vdot(fake,bar))
print(myalmdot(slmT, np.conj(bla), lmax, lmax, 0)/np.vdot(fake,bar))
......@@ -142,6 +142,9 @@ arr.apply([](T &v){v=0.;});
auto m1 = cube.template subarray<2>({supp,supp,0},{ntheta,nphi,0});
sharp_alm2map(a1.Alms().data(), m1.vdata(), *ginfo, *ainfo, 0, nthreads);
correct(m1,0);
// for (size_t i=1; i+1<ntheta; ++i)
// for (size_t j=0; j<nphi; ++j)
// m1.v(i,j)*=2;
}
for (size_t k=1; k<=kmax; ++k)
{
......@@ -357,6 +360,11 @@ arr.apply([](T &v){v=0.;});
{
auto m1 = cube.template subarray<2>({supp,supp,0},{ntheta,nphi,0});
decorrect(m1,0);
for (size_t j=0; j<nphi0; ++j)
{
m1.v(0,j)*=0.5;;
m1.v(ntheta0-1,j)*=0.5;;
}
sharp_alm2map_adjoint(a1.Alms().vdata(), m1.data(), *ginfo, *ainfo, 0, nthreads);
for (size_t m=0; m<=lmax; ++m)
for (size_t l=m; l<=lmax; ++l)
......@@ -369,6 +377,16 @@ arr.apply([](T &v){v=0.;});
auto m2 = cube.template subarray<2>({supp,supp,2*k },{ntheta,nphi,0});
decorrect(m1,k);
decorrect(m2,k);
for (size_t j=0; j<nphi0; ++j)
{
m1.v(0,j)*=0.5;;
m1.v(ntheta0-1,j)*=0.5;;
}
for (size_t j=0; j<nphi0; ++j)
{
m2.v(0,j)*=0.5;;
m2.v(ntheta0-1,j)*=0.5;;
}
sharp_alm2map_spin_adjoint(k, a1.Alms().vdata(), a2.Alms().vdata(), m1.data(),
m2.data(), *ginfo, *ainfo, 0, nthreads);
......@@ -379,7 +397,7 @@ arr.apply([](T &v){v=0.;});
{
auto tmp = -2.*conj(blmT(l,k))*T(lnorm[l]);
slmT(l,m) += conj(a1(l,m))*tmp.real();
slmT(l,m) += conj(a2(l,m))*tmp.imag();
slmT(l,m) -= conj(a2(l,m))*tmp.imag();
}
}
}
......
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