Commit ef3bd6c6 authored by Martin Reinecke's avatar Martin Reinecke

move uncategorized functions to

parent 2fdb003b
...@@ -157,3 +157,11 @@ as the `wgridder` component. ...@@ -157,3 +157,11 @@ as the `wgridder` component.
- in combination these two aspects allow extremely accurate gridding/degridding - in combination these two aspects allow extremely accurate gridding/degridding
operations (L2 error compared to explicit DFTs can go below 1e-12) with operations (L2 error compared to explicit DFTs can go below 1e-12) with
reasonable resource consumption reasonable resource consumption
ducc.misc
---------
Various unsorted functionality which will hopefully be categorized in the
future.
import ducc_0_1.sht as sht import ducc_0_1.sht as sht
import ducc_0_1.misc as misc
import numpy as np import numpy as np
from time import time from time import time
...@@ -33,7 +34,7 @@ map = job.alm2map(alm) ...@@ -33,7 +34,7 @@ map = job.alm2map(alm)
print("time for map synthesis: {}s".format(time()-t0)) print("time for map synthesis: {}s".format(time()-t0))
nlat2 = 2*lmax+3 nlat2 = 2*lmax+3
t0=time() t0=time()
map2 = sht.upsample_to_cc(map.reshape((nlat,nlon)), nlat2, False, False) map2 = misc.upsample_to_cc(map.reshape((nlat,nlon)), nlat2, False, False)
print("time for upsampling: {}s".format(time()-t0)) print("time for upsampling: {}s".format(time()-t0))
job.set_cc_geometry(nlat2, nlon) job.set_cc_geometry(nlat2, nlon)
t0=time() t0=time()
...@@ -55,7 +56,7 @@ map = job.alm2map(alm) ...@@ -55,7 +56,7 @@ map = job.alm2map(alm)
print("time for map synthesis: {}s".format(time()-t0)) print("time for map synthesis: {}s".format(time()-t0))
nlat2 = 2*lmax+3 nlat2 = 2*lmax+3
t0=time() t0=time()
map2 = sht.upsample_to_cc(map.reshape((nlat,nlon)), nlat2, True, True) map2 = misc.upsample_to_cc(map.reshape((nlat,nlon)), nlat2, True, True)
print("time for upsampling: {}s".format(time()-t0)) print("time for upsampling: {}s".format(time()-t0))
job.set_cc_geometry(nlat2, nlon) job.set_cc_geometry(nlat2, nlon)
t0=time() t0=time()
...@@ -78,7 +79,7 @@ map = job.alm2map(alm) ...@@ -78,7 +79,7 @@ map = job.alm2map(alm)
print("time for map synthesis: {}s".format(time()-t0)) print("time for map synthesis: {}s".format(time()-t0))
nlat2 = 2*lmax+3 nlat2 = 2*lmax+3
t0=time() t0=time()
map2 = sht.upsample_to_cc(map.reshape((nlat,nlon)), nlat2, False, True) map2 = misc.upsample_to_cc(map.reshape((nlat,nlon)), nlat2, False, True)
print("time for upsampling: {}s".format(time()-t0)) print("time for upsampling: {}s".format(time()-t0))
job.set_cc_geometry(nlat2, nlon) job.set_cc_geometry(nlat2, nlon)
t0=time() t0=time()
......
...@@ -31,8 +31,12 @@ ...@@ -31,8 +31,12 @@
#include <pybind11/numpy.h> #include <pybind11/numpy.h>
#include <vector> #include <vector>
#include "mr_util/infra/mav.h"
#include "mr_util/math/fft.h"
#include "mr_util/math/constants.h" #include "mr_util/math/constants.h"
#include "mr_util/math/gl_integrator.h" #include "mr_util/math/gl_integrator.h"
#include "mr_util/bindings/pybind_utils.h"
#include "python/alm.h"
namespace mr { namespace mr {
...@@ -65,6 +69,93 @@ a_d_c GL_thetas(int64_t nlat) ...@@ -65,6 +69,93 @@ a_d_c GL_thetas(int64_t nlat)
return res; return res;
} }
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 move(alm);
}
void upsample_to_cc(const mav<double,2> &in, bool has_np, bool has_sp,
mav<double,2> &out)
{
size_t ntheta_in = in.shape(0),
ntheta_out = out.shape(0),
nphi = in.shape(1);
MR_assert(out.shape(1)==nphi, "phi dimensions must be equal");
MR_assert((nphi&1)==0, "nphi must be even");
size_t nrings_in = 2*ntheta_in-has_np-has_sp;
size_t nrings_out = 2*ntheta_out-2;
MR_assert(nrings_out>=nrings_in, "number of rings must increase");
constexpr size_t delta=128;
for (size_t js=0; js<nphi; js+=delta)
{
size_t je = min(js+delta, nphi);
mav<double,2> tmp({nrings_out,je-js});
fmav<double> ftmp(tmp);
mav<double,2> tmp2(tmp.vdata(),{nrings_in, je-js}, true);
fmav<double> ftmp2(tmp2);
// enhance to "double sphere"
if (has_np)
for (size_t j=js; j<je; ++j)
tmp2.v(0,j-js) = in(0,j);
if (has_sp)
for (size_t j=js; j<je; ++j)
tmp2.v(ntheta_in-1,j-js) = in(ntheta_in-1,j);
for (size_t i=has_np, i2=nrings_in-1; i+has_sp<ntheta_in; ++i,--i2)
for (size_t j=js,j2=js+nphi/2; j<je; ++j,++j2)
{
if (j2>=nphi) j2-=nphi;
tmp2.v(i,j-js) = in(i,j);
tmp2.v(i2,j-js) = in(i,j2);
}
// FFT in theta direction
r2r_fftpack(ftmp2,ftmp2,{0},true,true,1./nrings_in,0);
if (!has_np) // shift
{
double ang = -pi/nrings_in;
for (size_t i=1; i<ntheta_in; ++i)
{
complex<double> rot(cos(i*ang),sin(i*ang));
for (size_t j=js; j<je; ++j)
{
complex<double> ctmp(tmp2(2*i-1,j-js),tmp2(2*i,j-js));
ctmp *= rot;
tmp2.v(2*i-1,j-js) = ctmp.real();
tmp2.v(2*i ,j-js) = ctmp.imag();
}
}
}
// zero-padding
for (size_t i=nrings_in; i<nrings_out; ++i)
for (size_t j=js; j<je; ++j)
tmp.v(i,j-js) = 0;
// FFT back
r2r_fftpack(ftmp,ftmp,{0},false,false,1.,0);
// copy to output map
for (size_t i=0; i<ntheta_out; ++i)
for (size_t j=js; j<je; ++j)
out.v(i,j) = tmp(i,j-js);
}
}
py::array py_upsample_to_cc(const py::array &in, size_t nrings_out, bool has_np,
bool has_sp, py::object &out_)
{
auto in2 = to_mav<double,2>(in);
auto out = get_optional_Pyarr<double>(out_, {nrings_out,size_t(in.shape(1))});
auto out2 = to_mav<double,2>(out,true);
MR_assert(out2.writable(),"x1");
upsample_to_cc(in2, has_np, has_sp, out2);
return move(out);
}
const char *misc_DS = R"""( const char *misc_DS = R"""(
Various unsorted utilities Various unsorted utilities
)"""; )""";
...@@ -77,6 +168,12 @@ void add_misc(py::module &msup) ...@@ -77,6 +168,12 @@ void add_misc(py::module &msup)
m.def("GL_weights",&GL_weights, "nlat"_a, "nlon"_a); m.def("GL_weights",&GL_weights, "nlat"_a, "nlon"_a);
m.def("GL_thetas",&GL_thetas, "nlat"_a); m.def("GL_thetas",&GL_thetas, "nlat"_a);
m.def("rotate_alm", &pyrotate_alm<double>, "alm"_a, "lmax"_a, "psi"_a, "theta"_a,
"phi"_a);
m.def("upsample_to_cc",&py_upsample_to_cc, "in"_a, "nrings_out"_a,
"has_np"_a, "has_sp"_a, "out"_a=py::none());
} }
} }
......
...@@ -37,8 +37,6 @@ ...@@ -37,8 +37,6 @@
#include "mr_util/sharp/sharp_almhelpers.h" #include "mr_util/sharp/sharp_almhelpers.h"
#include "mr_util/infra/string_utils.h" #include "mr_util/infra/string_utils.h"
#include "mr_util/infra/error_handling.h" #include "mr_util/infra/error_handling.h"
#include "mr_util/infra/mav.h"
#include "mr_util/math/fft.h"
#include "mr_util/math/constants.h" #include "mr_util/math/constants.h"
#include "mr_util/bindings/pybind_utils.h" #include "mr_util/bindings/pybind_utils.h"
...@@ -193,80 +191,6 @@ Python interface for some of the libsharp functionality ...@@ -193,80 +191,6 @@ Python interface for some of the libsharp functionality
Error conditions are reported by raising exceptions. Error conditions are reported by raising exceptions.
)"""; )""";
void upsample_to_cc(const mav<double,2> &in, bool has_np, bool has_sp,
mav<double,2> &out)
{
size_t ntheta_in = in.shape(0),
ntheta_out = out.shape(0),
nphi = in.shape(1);
MR_assert(out.shape(1)==nphi, "phi dimensions must be equal");
MR_assert((nphi&1)==0, "nphi must be even");
size_t nrings_in = 2*ntheta_in-has_np-has_sp;
size_t nrings_out = 2*ntheta_out-2;
MR_assert(nrings_out>=nrings_in, "number of rings must increase");
constexpr size_t delta=128;
for (size_t js=0; js<nphi; js+=delta)
{
size_t je = min(js+delta, nphi);
mav<double,2> tmp({nrings_out,je-js});
fmav<double> ftmp(tmp);
mav<double,2> tmp2(tmp.vdata(),{nrings_in, je-js}, true);
fmav<double> ftmp2(tmp2);
// enhance to "double sphere"
if (has_np)
for (size_t j=js; j<je; ++j)
tmp2.v(0,j-js) = in(0,j);
if (has_sp)
for (size_t j=js; j<je; ++j)
tmp2.v(ntheta_in-1,j-js) = in(ntheta_in-1,j);
for (size_t i=has_np, i2=nrings_in-1; i+has_sp<ntheta_in; ++i,--i2)
for (size_t j=js,j2=js+nphi/2; j<je; ++j,++j2)
{
if (j2>=nphi) j2-=nphi;
tmp2.v(i,j-js) = in(i,j);
tmp2.v(i2,j-js) = in(i,j2);
}
// FFT in theta direction
r2r_fftpack(ftmp2,ftmp2,{0},true,true,1./nrings_in,0);
if (!has_np) // shift
{
double ang = -pi/nrings_in;
for (size_t i=1; i<ntheta_in; ++i)
{
complex<double> rot(cos(i*ang),sin(i*ang));
for (size_t j=js; j<je; ++j)
{
complex<double> ctmp(tmp2(2*i-1,j-js),tmp2(2*i,j-js));
ctmp *= rot;
tmp2.v(2*i-1,j-js) = ctmp.real();
tmp2.v(2*i ,j-js) = ctmp.imag();
}
}
}
// zero-padding
for (size_t i=nrings_in; i<nrings_out; ++i)
for (size_t j=js; j<je; ++j)
tmp.v(i,j-js) = 0;
// FFT back
r2r_fftpack(ftmp,ftmp,{0},false,false,1.,0);
// copy to output map
for (size_t i=0; i<ntheta_out; ++i)
for (size_t j=js; j<je; ++j)
out.v(i,j) = tmp(i,j-js);
}
}
py::array py_upsample_to_cc(const py::array &in, size_t nrings_out, bool has_np,
bool has_sp, py::object &out_)
{
auto in2 = to_mav<double,2>(in);
auto out = get_optional_Pyarr<double>(out_, {nrings_out,size_t(in.shape(1))});
auto out2 = to_mav<double,2>(out,true);
MR_assert(out2.writable(),"x1");
upsample_to_cc(in2, has_np, has_sp, out2);
return move(out);
}
void add_sht(py::module &msup) void add_sht(py::module &msup)
{ {
using namespace pybind11::literals; using namespace pybind11::literals;
...@@ -299,8 +223,6 @@ void add_sht(py::module &msup) ...@@ -299,8 +223,6 @@ void add_sht(py::module &msup)
.def("alm2map_spin", &py_sharpjob<double>::alm2map_spin,"alm"_a,"spin"_a) .def("alm2map_spin", &py_sharpjob<double>::alm2map_spin,"alm"_a,"spin"_a)
.def("map2alm_spin", &py_sharpjob<double>::map2alm_spin,"map"_a,"spin"_a) .def("map2alm_spin", &py_sharpjob<double>::map2alm_spin,"map"_a,"spin"_a)
.def("__repr__", &py_sharpjob<double>::repr); .def("__repr__", &py_sharpjob<double>::repr);
m.def("upsample_to_cc",&py_upsample_to_cc, "in"_a, "nrings_out"_a,
"has_np"_a, "has_sp"_a, "out"_a=py::none());
} }
} }
......
...@@ -3,6 +3,7 @@ import pytest ...@@ -3,6 +3,7 @@ import pytest
from numpy.testing import assert_ from numpy.testing import assert_
import ducc_0_1.totalconvolve as totalconvolve import ducc_0_1.totalconvolve as totalconvolve
import ducc_0_1.sht as sht import ducc_0_1.sht as sht
import ducc_0_1.misc as misc
pmp = pytest.mark.parametrize pmp = pytest.mark.parametrize
...@@ -74,7 +75,7 @@ def test_against_convolution(lkmax, ncomp, separate): ...@@ -74,7 +75,7 @@ def test_against_convolution(lkmax, ncomp, separate):
res2 = np.zeros((nptg, ncomp)) res2 = np.zeros((nptg, ncomp))
for c in range(ncomp): for c in range(ncomp):
for i in range(nptg): for i in range(nptg):
rbeam=totalconvolve.rotate_alm(blm2[:,c], lmax, ptg[i,2],ptg[i,0],ptg[i,1]) rbeam=misc.rotate_alm(blm2[:,c], lmax, ptg[i,2],ptg[i,0],ptg[i,1])
res2[i,c] = convolve(slm[:,c], rbeam, lmax).real res2[i,c] = convolve(slm[:,c], rbeam, lmax).real
if separate: if separate:
_assert_close(res1, res2, 1e-7) _assert_close(res1, res2, 1e-7)
......
...@@ -5,7 +5,7 @@ ...@@ -5,7 +5,7 @@
#include <pybind11/pybind11.h> #include <pybind11/pybind11.h>
#include <pybind11/numpy.h> #include <pybind11/numpy.h>
#include "totalconvolve.h" #include "python/totalconvolve.h"
namespace mr { namespace mr {
...@@ -79,20 +79,6 @@ void makevec_v(py::array &inp, int64_t lmax, int64_t kmax, vector<Alm<complex<T> ...@@ -79,20 +79,6 @@ void makevec_v(py::array &inp, int64_t lmax, int64_t kmax, vector<Alm<complex<T>
} }
}; };
#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 move(alm);
}
#endif
constexpr const char *totalconvolve_DS = R"""( constexpr const char *totalconvolve_DS = R"""(
Python interface for total convolution/interpolation library Python interface for total convolution/interpolation library
...@@ -262,10 +248,6 @@ void add_totalconvolve(py::module &msup) ...@@ -262,10 +248,6 @@ void add_totalconvolve(py::module &msup)
.def ("deinterpol", &inter_f::pydeinterpol, deinterpol_DS, "ptg"_a, "data"_a) .def ("deinterpol", &inter_f::pydeinterpol, deinterpol_DS, "ptg"_a, "data"_a)
.def ("getSlm", &inter_f::pygetSlm, getSlm_DS, "beam"_a) .def ("getSlm", &inter_f::pygetSlm, getSlm_DS, "beam"_a)
.def ("support", &inter_f::support); .def ("support", &inter_f::support);
#if 1
m.def("rotate_alm", &pyrotate_alm<double>, "alm"_a, "lmax"_a, "psi"_a, "theta"_a,
"phi"_a);
#endif
m.def("epsilon_guess", &epsilon_guess, "support"_a, "ofactor"_a); m.def("epsilon_guess", &epsilon_guess, "support"_a, "ofactor"_a);
} }
......
...@@ -20,9 +20,8 @@ ...@@ -20,9 +20,8 @@
#include "mr_util/sharp/sharp.h" #include "mr_util/sharp/sharp.h"
#include "mr_util/sharp/sharp_almhelpers.h" #include "mr_util/sharp/sharp_almhelpers.h"
#include "mr_util/sharp/sharp_geomhelpers.h" #include "mr_util/sharp/sharp_geomhelpers.h"
#include "alm.h" #include "python/alm.h"
#include "mr_util/math/fft.h" #include "mr_util/math/fft.h"
#include "mr_util/bindings/pybind_utils.h"
namespace mr { namespace mr {
...@@ -116,7 +115,7 @@ template<typename T> void convolve_1d(const fmav<T> &in, ...@@ -116,7 +115,7 @@ template<typename T> void convolve_1d(const fmav<T> &in,
using detail_fft::convolve_1d; using detail_fft::convolve_1d;
namespace detail_interpol_ng { namespace detail_totalconvolve {
using namespace std; using namespace std;
...@@ -904,8 +903,8 @@ double epsilon_guess(size_t support, double ofactor) ...@@ -904,8 +903,8 @@ double epsilon_guess(size_t support, double ofactor)
} }
using detail_interpol_ng::Interpolator; using detail_totalconvolve::Interpolator;
using detail_interpol_ng::epsilon_guess; using detail_totalconvolve::epsilon_guess;
} }
......
...@@ -22,7 +22,7 @@ ...@@ -22,7 +22,7 @@
#include <pybind11/pybind11.h> #include <pybind11/pybind11.h>
#include <pybind11/numpy.h> #include <pybind11/numpy.h>
#include "mr_util/bindings/pybind_utils.h" #include "mr_util/bindings/pybind_utils.h"
#include "gridder_cxx.h" #include "python/gridder_cxx.h"
namespace mr { namespace mr {
......
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