Skip to content

GitLab

Commit 975f7603 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

many fixes, but still not correct

parent e3f2e9b3
Hide whitespace changes
Inline Side-by-side
interpol_ng/alm.h
View file @ 975f7603
......@@ -8,7 +8,7 @@
#ifndef MRUTIL_ALM_H
#define MRUTIL_ALM_H
#if 0
#if 1
#include <complex>
#include <cmath>
#include "mr_util/infra/threading.h"
......@@ -143,7 +143,7 @@ template<typename T> class Alm: public Alm_Base
}
};
#if 0
#if 1
/*! Class for calculation of the Wigner matrix at arbitrary angles, using Risbo
recursion in a way that can be OpenMP-parallelised. This approach uses more
memory and is slightly slower than wigner_d_risbo_scalar. */
......@@ -262,7 +262,7 @@ size_t nthreads=1;
using detail_alm::Alm_Base;
using detail_alm::Alm;
#if 0
#if 1
using detail_alm::rotate_alm;
#endif
}
......
interpol_ng/demo.py
View file @ 975f7603
......@@ -4,7 +4,7 @@ import pysharp
import time
import matplotlib.pyplot as plt
np.random.seed(42)
np.random.seed(48)
def nalm(lmax, mmax):
return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)
......@@ -24,10 +24,17 @@ def deltabeam(lmax,kmax):
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=2047
lmax=10
mmax=lmax
kmax=2 # doesn't make any sense for the beam we are using, but just for demonstration purposes ...
kmax=lmax
# get random sky a_lm
......@@ -36,36 +43,35 @@ slmT = random_alm(lmax, mmax)
# build beam a_lm (pencil beam for now)
blmT = deltabeam(lmax,kmax)
blmT = random_alm(lmax, mmax)
#blmT[:]=0
#blmT[lmax+1]=1j
blmT[:].imag=0
t0=time.time()
# build interpolator object for slmT and blmT
foo = interpol_ng.PyInterpolator(slmT,blmT,lmax, kmax, epsilon=1e-6, nthreads=2)
print("setup time: ",time.time()-t0)
# assemble a set of (theta, phi, psi) pointings, at which we want to evaluate
# the interpolated signal.
# For testig purposes, we choose theta and phi such that they form a
# Gauss-Legendre grid of sufficient resolution to recover the input a_lm
nth = lmax+1
nth = 2*lmax+1
nph = 2*mmax+1
ptg = np.zeros((nth,nph,3))
th, _ = np.polynomial.legendre.leggauss(nth)
th = np.arccos(-th)
ptg[:,:,0] = th.reshape((-1,1))
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] = 0
ptg[:,:,2] = 1.
t0=time.time()
# do the actual interpolation
bar=foo.interpol(ptg.reshape((-1,3))).reshape((nth,nph))
print("interpolation time: ", time.time()-t0)
# get a_lm back from interpolation results
job = pysharp.sharpjob_d()
job.set_triangular_alm_info(lmax, mmax)
job.set_gauss_geometry(nth, nph)
alm2 = job.map2alm(bar.reshape((-1,)))
#compare with original a_lm
plt.plot(np.abs(alm2-slmT))
plt.suptitle("Deviations between original and reconstructed a_lm")
plt.subplot(2,2,1)
plt.imshow(bar.reshape((nth,nph)))
bar2 = np.zeros((nth,nph))
for ith in range(nth):
for iph in range(nph):
rbeam=interpol_ng.rotate_alm(blmT, lmax, ptg[ith,iph,2],ptg[ith,iph,0],ptg[ith,iph,1])
bar2[ith,iph] = convolve(slmT, rbeam, lmax)
plt.subplot(2,2,2)
plt.imshow(bar2)
plt.subplot(2,2,3)
plt.imshow((bar2-bar.reshape((nth,nph))))
plt.show()
interpol_ng/interpol_ng.cc
View file @ 975f7603
......@@ -17,7 +17,7 @@
#include "alm.h"
#include "mr_util/math/fft.h"
#include "mr_util/bindings/pybind_utils.h"
#include <iostream>
using namespace std;
using namespace mr;
......@@ -37,8 +37,9 @@ template<typename T> class Interpolator
ES_Kernel kernel;
mav<T,3> cube; // the data cube (theta, phi, 2*mbeam+1[, IQU])
void correct(mav<T,2> &arr)
void correct(mav<T,2> &arr, int spin)
{
double sfct = (spin&1) ? -1 : 1;
mav<T,2> tmp({nphi,nphi});
auto tmp0=tmp.template subarray<2>({0,0},{nphi0, nphi0});
fmav<T> ftmp0(tmp0);
......@@ -50,10 +51,10 @@ template<typename T> class Interpolator
for (size_t j=0,j2=nphi0/2; j<nphi0; ++j,++j2)
{
if (j2>=nphi0) j2-=nphi0;
tmp0.v(i2,j) = arr(i,j2);
tmp0.v(i2,j) = sfct*arr(i,j2);
}
// FFT to frequency domain on minimal grid
r2r_fftpack(ftmp0,ftmp0,{0,1},true,true,1./(nphi0*nphi0),nthreads);
r2r_fftpack(ftmp0,ftmp0,{1,0},true,true,1./(nphi0*nphi0),nthreads);
auto fct = kernel.correction_factors(nphi, nphi0/2+1, nthreads);
for (size_t i=0; i<nphi0; ++i)
for (size_t j=0; j<nphi0; ++j)
......@@ -76,7 +77,7 @@ template<typename T> class Interpolator
kmax(blmT.Mmax()),
nphi0(2*good_size_real(lmax+1)),
ntheta0(nphi0/2+1),
nphi(2*good_size_real(size_t((2*lmax+1)*ofmin/2.))),
nphi(max<size_t>(20,2*good_size_real(size_t((2*lmax+1)*ofmin/2.)))),
ntheta(nphi/2+1),
nthreads(nthreads_),
ofactor(double(nphi)/(2*lmax+1)),
......@@ -115,29 +116,26 @@ template<typename T> class Interpolator
}
size_t kidx1 = (k==0) ? 0 : 2*k-1,
kidx2 = (k==0) ? 0 : 2*k;
auto quadrant=k%4;
if (quadrant&1)
swap(kidx1, kidx2);
auto m1 = cube.template subarray<2>({supp,supp,kidx1},{ntheta,nphi,0});
auto m2 = cube.template subarray<2>({supp,supp,kidx2},{ntheta,nphi,0});
if (k==0)
sharp_alm2map(a1.Alms().data(), m1.vdata(), *ginfo, *ainfo, nthreads, nullptr, nullptr);
sharp_alm2map(a1.Alms().data(), m1.vdata(), *ginfo, *ainfo, 0, nullptr, nullptr);
else
sharp_alm2map_spin(k, a1.Alms().data(), a2.Alms().data(), m1.vdata(), m2.vdata(), *ginfo, *ainfo, nthreads, nullptr, nullptr);
correct(m1);
if (k!=0) correct(m2);
if ((quadrant==1)||(quadrant==2)) m1.apply([](T &v){v=-v;});
if ((k>0) &&((quadrant==0)||(quadrant==1))) m2.apply([](T &v){v=-v;});
sharp_alm2map_spin(k, a1.Alms().data(), a2.Alms().data(), m1.vdata(), m2.vdata(), *ginfo, *ainfo, 0, nullptr, nullptr);
correct(m1,k);
if (k!=0) correct(m2,k);
if (k!=0)
{ m1.apply([](T &v){v*=2;}); m2.apply([](T &v){v*=2;}); }
}
// fill border regions
for (size_t i=0; i<supp; ++i)
for (size_t j=0, j2=nphi/2; j<nphi; ++j,++j2)
for (size_t k=0; k<cube.shape(2); ++k)
{
double fct = (((k+1)/2)&1) ? -1 : 1;
if (j2>=nphi) j2-=nphi;
cube.v(supp-1-i,j2+supp,k) = cube(supp+1+i,j+supp,k);
cube.v(supp+ntheta+i,j2+supp,k) = cube(supp+ntheta-2-i, j+supp,k);
cube.v(supp-1-i,j2+supp,k) = fct*cube(supp+1+i,j+supp,k);
cube.v(supp+ntheta+i,j2+supp,k) = fct*cube(supp+ntheta-2-i, j+supp,k);
}
for (size_t i=0; i<ntheta+2*supp; ++i)
for (size_t j=0; j<supp; ++j)
......@@ -164,8 +162,8 @@ template<typename T> class Interpolator
vector<vector<size_t>> mapper(nct*ncp);
for (size_t i=0; i<ptg.shape(0); ++i)
{
size_t itheta=size_t(ptg(i,0)/pi*nct),
iphi=size_t(ptg(i,1)/(2*pi)*ncp);
size_t itheta=min(nct-1,size_t(ptg(i,0)/pi*nct)),
iphi=min(ncp-1,size_t(ptg(i,1)/(2*pi)*ncp));
mapper[itheta*ncp+iphi].push_back(i);
}
size_t cnt=0;
......@@ -177,7 +175,6 @@ template<typename T> class Interpolator
for (size_t i=0; i<idx.size(); ++i)
idx[i]=i;
#endif
execStatic(idx.size(), nthreads, 0, [&](Scheduler &sched)
{
vector<T> wt(supp), wp(supp);
......@@ -195,12 +192,13 @@ template<typename T> class Interpolator
wp[t] = kernel((t+i1-f1)*delta - 1);
double val=0;
psiarr[0]=1.;
double cpsi=cos(ptg(i,2)), spsi=sin(ptg(i,2));
double psi=ptg(i,2);
double cpsi=cos(psi), spsi=sin(psi);
double cnpsi=cpsi, snpsi=spsi;
for (size_t l=1; l<=kmax; ++l)
{
psiarr[2*l-1]=cnpsi;
psiarr[2*l]=-snpsi;
psiarr[2*l]=snpsi;
const double tmp = snpsi*cpsi + cnpsi*spsi;
cnpsi=cnpsi*cpsi - snpsi*spsi;
snpsi=tmp;
......@@ -233,13 +231,17 @@ template<typename T> class PyInterpolator: public Interpolator<T>
}
};
#if 0
#if 1
template<typename T> py::array pyrotate_alm(const py::array &alm_, int64_t lmax,
double psi, double theta, double phi)
{
Alm<complex<T>> alm(to_mav<complex<T>,1>(alm_), lmax, lmax);
rotate_alm(alm, psi, theta, phi);
return alm_;
auto a1 = to_mav<complex<T>,1>(alm_);
auto alm = make_Pyarr<complex<T>>({a1.shape(0)});
auto a2 = to_mav<complex<T>,1>(alm,true);
for (size_t i=0; i<a1.shape(0); ++i) a2.v(i)=a1(i);
auto blah = Alm<complex<T>>(a2,lmax,lmax);
rotate_alm(blah, psi, theta, phi);
return alm;
}
#endif
......@@ -253,7 +255,7 @@ PYBIND11_MODULE(interpol_ng, m)
.def(py::init<const py::array &, const py::array &, int64_t, int64_t, double, int>(),
"sky"_a, "beam"_a, "lmax"_a, "kmax"_a, "epsilon"_a, "nthreads"_a)
.def ("interpol", &PyInterpolator<double>::interpol, "ptg"_a);
#if 0
#if 1
m.def("rotate_alm", &pyrotate_alm<double>, "alm"_a, "lmax"_a, "psi"_a, "theta"_a,
"phi"_a);
#endif
......
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