Skip to content
GitLab
Commit 4ae3295b authored by Martin Reinecke's avatar Martin Reinecke
Browse files

fixes, tweaks and documentation

parent 051b454d
Pipeline #72710 passed with stages
in 11 minutes and 53 seconds
Hide whitespace changes
Inline Side-by-side
pyinterpol_ng/demo.py
View file @ 4ae3295b
......@@ -10,11 +10,11 @@ def nalm(lmax, mmax):
return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)
def random_alm(lmax, mmax):
res = np.random.uniform(-1., 1., nalm(lmax, mmax)) \
+ 1j*np.random.uniform(-1., 1., nalm(lmax, mmax))
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))
# make a_lm with m==0 real-valued
res[0:lmax+1].imag = 0.
res[0:lmax+1,:].imag = 0.
return res
......@@ -41,19 +41,21 @@ def convolve(alm1, alm2, lmax):
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
slmT = random_alm(lmax, lmax)
slm = random_alm(lmax, lmax, ncomp)
# build beam a_lm
blmT = random_alm(lmax, kmax)
blm = random_alm(lmax, kmax, ncomp)
t0=time.time()
# build interpolator object for slmT and blmT
foo = pyinterpol_ng.PyInterpolator(slmT.reshape((-1,1)),blmT.reshape((-1,1)),True,lmax, kmax, epsilon=1e-6, nthreads=2)
foo = pyinterpol_ng.PyInterpolator(slm,blm,separate,lmax, kmax, epsilon=1e-6, nthreads=2)
print("setup time: ",time.time()-t0)
nth = lmax+1
nph = 2*lmax+1
......@@ -68,22 +70,22 @@ 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))
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.reshape((nth,nph)))
plt.imshow(bar[:,:,0])
bar2 = np.zeros((nth,nph))
blmTfull = np.zeros(slmT.size)+0j
blmTfull[0:blmT.size] = blmT
blmfull = np.zeros(slm.shape)+0j
blmfull[0:blm.shape[0],:] = blm
for ith in range(nth):
rbeamth=pyinterpol_ng.rotate_alm(blmTfull, lmax, ptg[ith,0,2],ptg[ith,0,0],0)
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(slmT, rbeam, lmax).real
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.reshape((nth,nph))))
plt.imshow(bar2-bar[:,:,0])
plt.show()
......@@ -91,12 +93,18 @@ ptg=np.random.uniform(0.,1.,3*1000000).reshape(1000000,3)
ptg[:,0]*=np.pi
ptg[:,1]*=2*np.pi
ptg[:,2]*=2*np.pi
foo = pyinterpol_ng.PyInterpolator(slmT.reshape((-1,1)),blmT.reshape((-1,1)),True,lmax, kmax, epsilon=1e-6, nthreads=2)
#foo = pyinterpol_ng.PyInterpolator(slm,blm,separate,lmax, kmax, epsilon=1e-6, nthreads=2)
t0=time.time()
bar=foo.interpol(ptg)
fake = np.random.uniform(0.,1., ptg.shape[0])
foo2 = pyinterpol_ng.PyInterpolator(lmax, kmax, 1, epsilon=1e-6, nthreads=2)
foo2.deinterpol(ptg.reshape((-1,3)), fake.reshape((-1,1)))
bla=foo2.getSlm(blmT.reshape((-1,1))).reshape(-1)
print(myalmdot(slmT, bla, lmax, lmax, 0))
print(np.vdot(fake,bar))
print(myalmdot(slmT, bla, lmax, lmax, 0)/np.vdot(fake,bar))
print("interpolation time: ", time.time()-t0)
fake = np.random.uniform(0.,1., (ptg.shape[0],ncomp2))
foo2 = pyinterpol_ng.PyInterpolator(lmax, kmax, ncomp2, epsilon=1e-6, nthreads=2)
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.)
pyinterpol_ng/interpol_ng.h
View file @ 4ae3295b
......@@ -155,6 +155,7 @@ template<typename T> class Interpolator
ncomp(separate ? slm.size() : 1),
cube({ntheta+2*supp, nphi+2*supp, 2*kmax+1, ncomp})
{
MR_assert((ncomp==1)||(ncomp==3), "currently only 1 or 3 components allowed");
MR_assert(slm.size()==blm.size(), "inconsistent slm and blm vectors");
for (size_t i=0; i<slm.size(); ++i)
{
......@@ -265,6 +266,7 @@ template<typename T> class Interpolator
ncomp(ncomp_),
cube({ntheta+2*supp, nphi+2*supp, 2*kmax+1, ncomp_})
{
MR_assert((ncomp==1)||(ncomp==3), "currently only 1 or 3 components allowed");
MR_assert((supp<=ntheta) && (supp<=nphi), "support too large!");
cube.apply([](T &v){v=0.;});
}
......@@ -283,7 +285,6 @@ template<typename T> class Interpolator
{
vector<T> wt(supp), wp(supp);
vector<T> psiarr(2*kmax+1);
vector<T> val(ncomp);
while (auto rng=sched.getNext()) for(auto ind=rng.lo; ind<rng.hi; ++ind)
{
size_t i=idx[ind];
......@@ -307,15 +308,37 @@ template<typename T> class Interpolator
cnpsi=cnpsi*cpsi - snpsi*spsi;
snpsi=tmp;
}
for (size_t m=0; m<ncomp; ++m)
val[m]=0;
for (size_t j=0; j<supp; ++j)
for (size_t k=0; k<supp; ++k)
for (size_t l=0; l<2*kmax+1; ++l)
for (size_t m=0; m<ncomp; ++m)
val[m] += cube(i0+j,i1+k,l,m)*wt[j]*wp[k]*psiarr[l];
for (size_t m=0; m<ncomp; ++m)
res.v(i,m) = val[m];
if (ncomp==1)
{
double vv=0.;
for (size_t j=0; j<supp; ++j)
for (size_t k=0; k<supp; ++k)
{
double t0=wt[j]*wp[k];
for (size_t l=0; l<2*kmax+1; ++l)
vv += cube(i0+j,i1+k,l,0)*t0*psiarr[l];
}
res.v(i,0) = vv;
}
else // ncomp==3
{
double v0=0., v1=0., v2=0.;
for (size_t j=0; j<supp; ++j)
for (size_t k=0; k<supp; ++k)
{
double t0 = wt[j]*wp[k];
for (size_t l=0; l<2*kmax+1; ++l)
{
auto tmp = t0*psiarr[l];
v0 += cube(i0+j,i1+k,l,0)*tmp;
v1 += cube(i0+j,i1+k,l,1)*tmp;
v2 += cube(i0+j,i1+k,l,2)*tmp;
}
}
res.v(i,0) = v0;
res.v(i,1) = v1;
res.v(i,2) = v2;
}
}
});
}
......@@ -341,7 +364,6 @@ template<typename T> class Interpolator
size_t b_theta=99999999999999, b_phi=9999999999999999;
vector<T> wt(supp), wp(supp);
vector<T> psiarr(2*kmax+1);
vector<T> val(ncomp);
while (auto rng=sched.getNext()) for(auto ind=rng.lo; ind<rng.hi; ++ind)
{
size_t i=idx[ind];
......@@ -353,8 +375,6 @@ template<typename T> class Interpolator
size_t i1 = size_t(f1+1.);
for (size_t t=0; t<supp; ++t)
wp[t] = kernel((t+i1-f1)*delta - 1);
for (size_t m=0; m<ncomp; ++m)
val[m] = data(i,m);
psiarr[0]=1.;
double psi=ptg(i,2);
double cpsi=cos(psi), spsi=sin(psi);
......@@ -385,11 +405,30 @@ template<typename T> class Interpolator
locks.v(b_theta+1,b_phi).lock();
locks.v(b_theta+1,b_phi+1).lock();
}
for (size_t j=0; j<supp; ++j)
for (size_t k=0; k<supp; ++k)
for (size_t l=0; l<2*kmax+1; ++l)
for (size_t m=0; m<ncomp; ++m)
cube.v(i0+j,i1+k,l,m) += val[m]*wt[j]*wp[k]*psiarr[l];
if (ncomp==1)
{
double val = data(i,0);
for (size_t j=0; j<supp; ++j)
for (size_t k=0; k<supp; ++k)
for (size_t l=0; l<2*kmax+1; ++l)
cube.v(i0+j,i1+k,l,0) += val*wt[j]*wp[k]*psiarr[l];
}
else // ncomp==3
{
double v0=data(i,0), v1=data(i,1), v2=data(i,2);
for (size_t j=0; j<supp; ++j)
for (size_t k=0; k<supp; ++k)
{
double t0 = wt[j]*wp[k];
for (size_t l=0; l<2*kmax+1; ++l)
{
double tmp = t0*psiarr[l];
cube.v(i0+j,i1+k,l,0) += v0*tmp;
cube.v(i0+j,i1+k,l,1) += v1*tmp;
cube.v(i0+j,i1+k,l,2) += v2*tmp;
}
}
}
}
if (b_theta<locks.shape(0)) // unlock
{
......
pyinterpol_ng/pyinterpol_ng.cc
View file @ 4ae3295b
......@@ -117,14 +117,14 @@ Constructor for interpolation mode
Parameters
----------
sky : numpy.ndarray(numpy.complex)
spherical harmonic coefficients of the sky
beam : numpy.ndarray(numpy.complex)
spherical harmonic coefficients of the sky
sky : numpy.ndarray((nalm_sky, ncomp), dtype=numpy.complex)
spherical harmonic coefficients of the sky. ncomp can be 1 or 3.
beam : numpy.ndarray((nalm_beam, ncomp), dtype=numpy.complex)
spherical harmonic coefficients of the beam. ncomp can be 1 or 3
separate : bool
whether contributions of individual components should be added together.
lmax : int
maximum l in the coefficient arays
mmax : int
maximum m in the sky coefficients
kmax : int
maximum azimuthal moment in the beam coefficients
epsilon : float
......@@ -139,10 +139,11 @@ Parameters
----------
lmax : int
maximum l in the coefficient arays
mmax : int
maximum m in the sky coefficients
kmax : int
maximum azimuthal moment in the beam coefficients
ncomp : int
the number of components which are going to input to `deinterpol`.
Can be 1 or 3.
epsilon : float
desired accuracy for the interpolation; a typical value is 1e-5
nthreads : the number of threads to use for computation
......@@ -153,7 +154,7 @@ Computes the interpolated values for a given set of angle triplets
Parameters
----------
ptg : numpy.ndarray(numpy.float64) of shape(N,3)
ptg : numpy.ndarray((N, 3), dtype=numpy.float64)
theta, phi and psi angles (in radian) for N pointings
theta must be in the range [0; pi]
phi must be in the range [0; 2pi]
......@@ -161,8 +162,10 @@ ptg : numpy.ndarray(numpy.float64) of shape(N,3)
Returns
-------
numpy.array(numpy.float64) of shape (N,)
numpy.array((N, n2), dtype=numpy.float64)
the interpolated values
n2 is either 1 (if separate=True was used in the constructor) or the
second dimension of the input slm and blm arrays (otherwise)
Notes
-----
......@@ -177,13 +180,14 @@ data cube.
Parameters
----------
ptg : numpy.ndarray(numpy.float64) of shape(N,3)
ptg : numpy.ndarray((N,3), dtype=numpy.float64)
theta, phi and psi angles (in radian) for N pointings
theta must be in the range [0; pi]
phi must be in the range [0; 2pi]
psi should be in the range [-2pi; 2pi]
data : numpy.ndarray(numpy.float64) of shape (N,)
data : numpy.ndarray((N, n2), dtype=numpy.float64)
the interpolated values
n2 must match the `ncomp` value specified in the constructor.
Notes
-----
......@@ -194,18 +198,19 @@ Notes
constexpr const char *getSlm_DS = R"""(
Returns a set of sky spherical hamonic coefficients resulting from adjoint
interplation
interpolation
Parameters
----------
blmT : numpy.array(numpy.complex)
beam : numpy.array(nalm_beam, nbeam), dtype=numpy.complex)
spherical harmonic coefficients of the beam with lmax and kmax defined
in the constructor call
nbeam must match the ncomp specified in the constructor, unless ncomp was 1.
Returns
-------
numpy.array(numpy.complex):
spherical harmonic coefficients of the sky with lmax and mmax defined
numpy.array(nalm_sky, nbeam), dtype=numpy.complex):
spherical harmonic coefficients of the sky with lmax defined
in the constructor call
Notes
......@@ -230,7 +235,7 @@ PYBIND11_MODULE(pyinterpol_ng, m)
"lmax"_a, "kmax"_a, "ncomp"_a, "epsilon"_a, "nthreads"_a)
.def ("interpol", &PyInterpolator<double>::pyinterpol, interpol_DS, "ptg"_a)
.def ("deinterpol", &PyInterpolator<double>::pydeinterpol, deinterpol_DS, "ptg"_a, "data"_a)
.def ("getSlm", &PyInterpolator<double>::pygetSlm, getSlm_DS, "blmT"_a);
.def ("getSlm", &PyInterpolator<double>::pygetSlm, getSlm_DS, "beam"_a);
#if 1
m.def("rotate_alm", &pyrotate_alm<double>, "alm"_a, "lmax"_a, "psi"_a, "theta"_a,
"phi"_a);
......
Supports Markdown
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