Skip to content

GitLab

Commit b9190aa6 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

rename etc.

parent 9e9294d4
Pipeline #70641 passed with stages
in 9 minutes and 5 seconds
Hide whitespace changes
Inline Side-by-side
interpol_ng/MANIFEST.in
View file @ b9190aa6
prune mr_util
include alm.h
include interpol_ng.h
include mr_util/math/fft1d.h
include mr_util/math/fft.h
include mr_util/infra/threading.h
......
interpol_ng/README.md
View file @ b9190aa6
interpol_ng
===========
pyinterpol_ng
=============
Library for high-accuracy 4pi convolution on the sphere, which generates a
total convolution data cube from a set of sky and beam `a_lm` and computes
......
interpol_ng/demo.py
View file @ b9190aa6
import interpol_ng
import pyinterpol_ng
import numpy as np
import pysharp
import time
......@@ -53,7 +53,7 @@ blmT = random_alm(lmax, kmax)
t0=time.time()
# build interpolator object for slmT and blmT
foo = interpol_ng.PyInterpolator(slmT,blmT,lmax, kmax, epsilon=1e-6, nthreads=2)
foo = pyinterpol_ng.PyInterpolator(slmT,blmT,lmax, kmax, epsilon=1e-6, nthreads=2)
print("setup time: ",time.time()-t0)
nth = lmax+1
nph = 2*lmax+1
......@@ -76,9 +76,9 @@ 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)
rbeamth=pyinterpol_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])
rbeam=pyinterpol_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)
......@@ -91,10 +91,10 @@ 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 = interpol_ng.PyInterpolator(slmT,blmT,lmax, kmax, epsilon=1e-6, nthreads=2)
foo = pyinterpol_ng.PyInterpolator(slmT,blmT,lmax, kmax, epsilon=1e-6, nthreads=2)
bar=foo.interpol(ptg)
fake = np.random.uniform(0.,1., ptg.shape[0])
foo2 = interpol_ng.PyInterpolator(lmax, kmax, epsilon=1e-6, nthreads=2)
foo2 = pyinterpol_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))
......
interpol_ng/interpol_ng.cc → interpol_ng/interpol_ng.h
View file @ b9190aa6
......@@ -3,10 +3,12 @@
* Author: Martin Reinecke
*/
#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>
#ifndef MRUTIL_INTERPOL_NG_H
#define MRUTIL_INTERPOL_NG_H
#include <vector>
#include <complex>
#include <cmath>
#include "mr_util/math/constants.h"
#include "mr_util/math/gl_integrator.h"
#include "mr_util/math/es_kernel.h"
......@@ -18,12 +20,11 @@
#include "mr_util/math/fft.h"
#include "mr_util/bindings/pybind_utils.h"
using namespace std;
using namespace mr;
namespace mr {
namespace py = pybind11;
namespace detail_interpol_ng {
namespace {
using namespace std;
constexpr double ofmin=1.5;
......@@ -143,7 +144,7 @@ template<typename T> class Interpolator
kmax(blmT.Mmax()),
nphi0(2*good_size_real(lmax+1)),
ntheta0(nphi0/2+1),
nphi(max<size_t>(20,2*good_size_real(size_t((2*lmax+1)*ofmin/2.)))),
nphi(std::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)),
......@@ -160,7 +161,7 @@ template<typename T> class Interpolator
vector<double>lnorm(lmax+1);
for (size_t i=0; i<=lmax; ++i)
lnorm[i]=sqrt(4*pi/(2*i+1.));
lnorm[i]=std::sqrt(4*pi/(2*i+1.));
{
for (size_t m=0; m<=lmax; ++m)
......@@ -216,7 +217,7 @@ template<typename T> class Interpolator
kmax(kmax_),
nphi0(2*good_size_real(lmax+1)),
ntheta0(nphi0/2+1),
nphi(max<size_t>(20,2*good_size_real(size_t((2*lmax+1)*ofmin/2.)))),
nphi(std::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)),
......@@ -391,7 +392,7 @@ template<typename T> class Interpolator
vector<double>lnorm(lmax+1);
for (size_t i=0; i<=lmax; ++i)
lnorm[i]=sqrt(4*pi/(2*i+1.));
lnorm[i]=std::sqrt(4*pi/(2*i+1.));
{
auto m1 = cube.template subarray<2>({supp,supp,0},{ntheta,nphi,0});
......@@ -425,79 +426,10 @@ template<typename T> class Interpolator
}
};
template<typename T> class PyInterpolator: public Interpolator<T>
{
protected:
using Interpolator<T>::lmax;
using Interpolator<T>::kmax;
using Interpolator<T>::interpolx;
using Interpolator<T>::deinterpolx;
using Interpolator<T>::getSlmx;
}
public:
PyInterpolator(const py::array &slmT, const py::array &blmT,
int64_t lmax, int64_t kmax, double epsilon, int nthreads=0)
: Interpolator<T>(Alm<complex<T>>(to_mav<complex<T>,1>(slmT), lmax, lmax),
Alm<complex<T>>(to_mav<complex<T>,1>(blmT), lmax, kmax),
epsilon, nthreads) {}
PyInterpolator(int64_t lmax, int64_t kmax, double epsilon, int nthreads=0)
: Interpolator<T>(lmax, kmax, epsilon, nthreads) {}
py::array interpol(const py::array &ptg) const
{
auto ptg2 = to_mav<T,2>(ptg);
auto res = make_Pyarr<double>({ptg2.shape(0)});
auto res2 = to_mav<double,1>(res,true);
interpolx(ptg2, res2);
return res;
}
using detail_interpol_ng::Interpolator;
void deinterpol(const py::array &ptg, const py::array &data)
{
auto ptg2 = to_mav<T,2>(ptg);
auto data2 = to_mav<T,1>(data);
deinterpolx(ptg2, data2);
}
py::array getSlm(const py::array &blmT_)
{
auto res = make_Pyarr<complex<T>>({Alm_Base::Num_Alms(lmax, lmax)});
Alm<complex<T>> blmT(to_mav<complex<T>,1>(blmT_, false), lmax, kmax);
auto slmT_=to_mav<complex<T>,1>(res, true);
Alm<complex<T>> slmT(slmT_, lmax, lmax);
getSlmx(blmT, slmT);
return res;
}
};
#if 1
template<typename T> py::array pyrotate_alm(const py::array &alm_, int64_t lmax,
double psi, double theta, double phi)
{
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
} // unnamed namespace
PYBIND11_MODULE(interpol_ng, m)
{
using namespace pybind11::literals;
}
py::class_<PyInterpolator<double>> (m, "PyInterpolator")
.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(py::init<int64_t, int64_t, double, int>(),
"lmax"_a, "kmax"_a, "epsilon"_a, "nthreads"_a)
.def ("interpol", &PyInterpolator<double>::interpol, "ptg"_a)
.def ("deinterpol", &PyInterpolator<double>::deinterpol, "ptg"_a, "data"_a)
.def ("getSlm", &PyInterpolator<double>::getSlm, "blmT"_a);
#if 1
m.def("rotate_alm", &pyrotate_alm<double>, "alm"_a, "lmax"_a, "psi"_a, "theta"_a,
"phi"_a);
#endif
}
interpol_ng/pyinterpol_ng.cc 0 → 100644
View file @ b9190aa6
/*
* Copyright (C) 2020 Max-Planck-Society
* Author: Martin Reinecke
*/
#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>
#include "interpol_ng.h"
using namespace std;
using namespace mr;
namespace py = pybind11;
namespace {
template<typename T> class PyInterpolator: public Interpolator<T>
{
protected:
using Interpolator<T>::lmax;
using Interpolator<T>::kmax;
using Interpolator<T>::interpolx;
using Interpolator<T>::deinterpolx;
using Interpolator<T>::getSlmx;
public:
PyInterpolator(const py::array &slmT, const py::array &blmT,
int64_t lmax, int64_t kmax, double epsilon, int nthreads=0)
: Interpolator<T>(Alm<complex<T>>(to_mav<complex<T>,1>(slmT), lmax, lmax),
Alm<complex<T>>(to_mav<complex<T>,1>(blmT), lmax, kmax),
epsilon, nthreads) {}
PyInterpolator(int64_t lmax, int64_t kmax, double epsilon, int nthreads=0)
: Interpolator<T>(lmax, kmax, epsilon, nthreads) {}
py::array interpol(const py::array &ptg) const
{
auto ptg2 = to_mav<T,2>(ptg);
auto res = make_Pyarr<double>({ptg2.shape(0)});
auto res2 = to_mav<double,1>(res,true);
interpolx(ptg2, res2);
return res;
}
void deinterpol(const py::array &ptg, const py::array &data)
{
auto ptg2 = to_mav<T,2>(ptg);
auto data2 = to_mav<T,1>(data);
deinterpolx(ptg2, data2);
}
py::array getSlm(const py::array &blmT_)
{
auto res = make_Pyarr<complex<T>>({Alm_Base::Num_Alms(lmax, lmax)});
Alm<complex<T>> blmT(to_mav<complex<T>,1>(blmT_, false), lmax, kmax);
auto slmT_=to_mav<complex<T>,1>(res, true);
Alm<complex<T>> slmT(slmT_, lmax, lmax);
getSlmx(blmT, slmT);
return res;
}
};
#if 1
template<typename T> py::array pyrotate_alm(const py::array &alm_, int64_t lmax,
double psi, double theta, double phi)
{
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 tmp = Alm<complex<T>>(a2,lmax,lmax);
rotate_alm(tmp, psi, theta, phi);
return alm;
}
#endif
} // unnamed namespace
PYBIND11_MODULE(pyinterpol_ng, m)
{
using namespace pybind11::literals;
py::class_<PyInterpolator<double>> (m, "PyInterpolator")
.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(py::init<int64_t, int64_t, double, int>(),
"lmax"_a, "kmax"_a, "epsilon"_a, "nthreads"_a)
.def ("interpol", &PyInterpolator<double>::interpol, "ptg"_a)
.def ("deinterpol", &PyInterpolator<double>::deinterpol, "ptg"_a, "data"_a)
.def ("getSlm", &PyInterpolator<double>::getSlm, "blmT"_a);
#if 1
m.def("rotate_alm", &pyrotate_alm<double>, "alm"_a, "lmax"_a, "psi"_a, "theta"_a,
"phi"_a);
#endif
}
interpol_ng/setup.py
View file @ b9190aa6
......@@ -32,9 +32,9 @@ else:
# if you don't want debugging info, add "-s" to python_module_link_args
def get_extension_modules():
return [Extension('interpol_ng',
return [Extension('pyinterpol_ng',
language='c++',
sources=['interpol_ng.cc',
sources=['pyinterpol_ng.cc',
'mr_util/infra/threading.cc',
'mr_util/sharp/sharp.cc',
'mr_util/sharp/sharp_core.cc',
......@@ -64,6 +64,7 @@ def get_extension_modules():
'mr_util/sharp/sharp_geomhelpers.h',
'mr_util/sharp/sharp_almhelpers.h',
'setup.py',
'interpol_ng.h',
'alm.h'],
include_dirs=include_dirs,
define_macros=define_macros,
......@@ -71,7 +72,7 @@ def get_extension_modules():
extra_link_args=python_module_link_args)]
setup(name='interpol_ng',
setup(name='pyinterpol_ng',
version='0.0.1',
description='Python bindings for total convolution and interpolation library',
include_package_data=True,
......
interpol_ng/test/test_interpol_ng.py
View file @ b9190aa6
import numpy as np
import pytest
from numpy.testing import assert_
import interpol_ng
import pyinterpol_ng
import pysharp
pmp = pytest.mark.parametrize
......@@ -57,7 +57,7 @@ def test_against_convolution(lkmax):
slmT = random_alm(lmax, lmax)
blmT = random_alm(lmax, kmax)
inter = interpol_ng.PyInterpolator(slmT, blmT, lmax, kmax, epsilon=1e-8,
inter = pyinterpol_ng.PyInterpolator(slmT, blmT, lmax, kmax, epsilon=1e-8,
nthreads=2)
nptg = 50
ptg = np.zeros((nptg,3))
......@@ -71,7 +71,7 @@ def test_against_convolution(lkmax):
blmT2 = np.zeros(nalm(lmax,lmax))+0j
blmT2[0:blmT.shape[0]] = blmT
for i in range(nptg):
rbeam=interpol_ng.rotate_alm(blmT2, lmax, ptg[i,2],ptg[i,0],ptg[i,1])
rbeam=pyinterpol_ng.rotate_alm(blmT2, lmax, ptg[i,2],ptg[i,0],ptg[i,1])
res2[i] = convolve(slmT, rbeam, lmax).real
_assert_close(res1, res2, 1e-7)
......@@ -85,10 +85,10 @@ def test_adjointness(lkmax):
ptg[:,0]*=np.pi
ptg[:,1]*=2*np.pi
ptg[:,2]*=2*np.pi
foo = interpol_ng.PyInterpolator(slmT,blmT,lmax, kmax, epsilon=1e-6, nthreads=2)
foo = pyinterpol_ng.PyInterpolator(slmT,blmT,lmax, kmax, epsilon=1e-6, nthreads=2)
inter1=foo.interpol(ptg)
fake = np.random.uniform(0.,1., ptg.shape[0])
foo2 = interpol_ng.PyInterpolator(lmax, kmax, epsilon=1e-6, nthreads=2)
foo2 = pyinterpol_ng.PyInterpolator(lmax, kmax, epsilon=1e-6, nthreads=2)
foo2.deinterpol(ptg.reshape((-1,3)), fake)
bla=foo2.getSlm(blmT)
_assert_close(myalmdot(slmT, bla, lmax, lmax, 0), np.vdot(fake,inter1), 1e-12)
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