/* * Copyright (C) 2020 Max-Planck-Society * Author: Martin Reinecke */ #include #include #include #include #include "mr_util/math/constants.h" #include "mr_util/math/gl_integrator.h" #include "mr_util/math/es_kernel.h" #include "mr_util/infra/mav.h" #include "mr_util/sharp/sharp.h" #include "mr_util/sharp/sharp_almhelpers.h" #include "mr_util/sharp/sharp_geomhelpers.h" #include "alm.h" #include "mr_util/math/fft.h" #include "mr_util/bindings/pybind_utils.h" using namespace std; using namespace mr; namespace py = pybind11; namespace { template class Interpolator { protected: size_t lmax, kmax, ppring, ntheta, nphi; double ofactor; size_t supp; ES_Kernel kernel; mav cube; // the data cube (theta, phi, 2*mbeam+1[, IQU]) void correct(mav &arr) { mav tmp({nphi,nphi}); fmav ftmp(tmp); for (size_t i=0; i=nphi) j2-=nphi; tmp.v(i2,j) = arr(i,j2); } // FFT r2r_fftpack(ftmp,ftmp,{0,1},true,true,1.,0); auto fct = kernel.correction_factors(nphi, nphi/2+2, 0); for (size_t i=0; i> &slmT, const Alm> &blmT, double epsilon) : lmax(slmT.Lmax()), kmax(blmT.Mmax()), ppring(2*good_size_real(2*lmax+1)), ntheta(ppring/2+1), nphi(ppring), ofactor(double(ppring)/(2*lmax+1)), supp(ES_Kernel::get_supp(epsilon, ofactor)), kernel(supp, ofactor, 1), cube({ntheta+2*supp, nphi+2*supp, 2*kmax+1}) { MR_assert((supp<=ntheta) && (supp<=nphi), "support too large!"); MR_assert(slmT.Mmax()==lmax, "Sky lmax must be equal to Sky mmax"); MR_assert(blmT.Lmax()==lmax, "Sky and beam lmax must be equal"); Alm> a1(lmax, lmax), a2(lmax,lmax); auto ginfo = sharp_make_cc_geom_info(ntheta,nphi,0,cube.stride(1),cube.stride(0)); auto ainfo = sharp_make_triangular_alm_info(lmax,lmax,1); vectorlnorm(lmax+1); for (size_t i=0; i<=lmax; ++i) lnorm[i]=sqrt(4*pi/(2*i+1.)); for (size_t k=0; k<=kmax; ++k) { double spinsign = (k==0) ? 1. : -1.; for (size_t m=0; m<=lmax; ++m) { T mfac=T((m&1) ? -1.:1.); for (size_t l=m; l<=lmax; ++l) { if (l v1=slmT(l,m)*blmT(l,k), v2=conj(slmT(l,m))*blmT(l,k)*mfac; a1(l,m) = (v1+conj(v2)*mfac)*T(0.5*spinsign*lnorm[l]); if (k>0) { complex tmp = (v1-conj(v2)*mfac)*T(-spinsign*0.5*lnorm[l]); a2(l,m) = complex(-tmp.imag(), tmp.real()); } } } } 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); mav m1(&cube.v(supp,supp,kidx1),{ntheta,nphi},{cube.stride(0),cube.stride(1)},true); mav m2(&cube.v(supp,supp,kidx2),{ntheta,nphi},{cube.stride(0),cube.stride(1)},true); if (k==0) 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, 0, 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;}); } // fill border regions for (size_t i=0; i=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); } for (size_t i=0; i &ptg, mav &res) const { MR_assert(ptg.shape(0)==res.shape(0), "dimension mismatch"); MR_assert(ptg.shape(1)==3, "second dimension must have length 3"); vector wt(supp), wp(supp); vector psiarr(2*kmax+1); for (size_t i=0; i class PyInterpolator: public Interpolator { public: PyInterpolator(const py::array &slmT, const py::array &blmT, int64_t lmax, int64_t kmax, double epsilon) : Interpolator(Alm>(to_mav,1>(slmT), lmax, lmax), Alm>(to_mav,1>(blmT), lmax, kmax), epsilon) {} using Interpolator::interpolx; py::array interpol(const py::array &ptg) { auto ptg2 = to_mav(ptg); auto res = make_Pyarr({ptg2.shape(0)}); auto res2 = to_mav(res,true); interpolx(ptg2, res2); return res; } }; } // unnamed namespace PYBIND11_MODULE(interpol_ng, m) { using namespace pybind11::literals; py::class_> (m, "PyInterpolator") .def(py::init(), "sky"_a, "beam"_a, "lmax"_a, "kmax"_a, "epsilon"_a) .def ("interpol", &PyInterpolator::interpol, "ptg"_a); }