Commit 382a00f2 authored by Martin Reinecke's avatar Martin Reinecke

add adjoint synthesis

parent 3eae85ce
Pipeline #97889 passed with stages
in 22 minutes and 26 seconds
......@@ -201,6 +201,37 @@ py::array Pyfull_synthesis(const py::array &alm_, size_t lmax,
phi0_, nphi_, ringstart_, spin, nthreads);
MR_fail("type matching failed: 'alm' has neither type 'c8' nor 'c16'");
}
template<typename T> py::array full_adjoint_synthesis2(py::array &alm_,
size_t lmax, const py::array &mval_, const py::array &mstart_,
const py::array &map_, const py::array &theta_, const py::array &phi0_,
const py::array &nphi_, const py::array &ringstart_, size_t spin,
size_t nthreads)
{
auto alm = to_mav<complex<T>,2>(alm_, true);
auto mval = to_mav<size_t,1>(mval_, false);
auto mstart = to_mav<size_t,1>(mstart_, false);
auto map = to_mav<T,2>(map_, false);
auto theta = to_mav<double,1>(theta_, false);
auto phi0 = to_mav<double,1>(phi0_, false);
auto nphi = to_mav<size_t,1>(nphi_, false);
auto ringstart = to_mav<size_t,1>(ringstart_, false);
adjoint_synthesis(alm, lmax, mval, mstart, map, theta, phi0, nphi, ringstart, spin, nthreads);
return alm_;
}
py::array Pyfull_adjoint_synthesis(py::array &alm_, size_t lmax,
const py::array &mval_, const py::array &mstart_,
const py::array &map_, const py::array &theta_, const py::array &phi0_,
const py::array &nphi_, const py::array &ringstart_, size_t spin,
size_t nthreads)
{
if (isPyarr<complex<float>>(alm_))
return full_adjoint_synthesis2<float>(alm_, lmax, mval_, mstart_, map_, theta_,
phi0_, nphi_, ringstart_, spin, nthreads);
else if (isPyarr<complex<double>>(alm_))
return full_adjoint_synthesis2<double>(alm_, lmax, mval_, mstart_, map_, theta_,
phi0_, nphi_, ringstart_, spin, nthreads);
MR_fail("type matching failed: 'alm' has neither type 'c8' nor 'c16'");
}
#endif
using a_d = py::array_t<double>;
......@@ -356,6 +387,7 @@ void add_sht(py::module_ &msup)
// m.def("synthesis", &Pysynthesis, "type"_a, "alm"_a, "map"_a, "lmax"_a, "mmax"_a, "spin"_a);
// m.def("synthesis", &Pysynthesis, "alm"_a, "map"_a, "lmax"_a, "mmax"_a, "spin"_a, "theta"_a, "nphi"_a, "phi0"_a, "offset"_a);
m.def("synthesis", &Pyfull_synthesis, "alm"_a, "lmax"_a, "mval"_a, "mstart"_a, "map"_a, "theta"_a, "phi0"_a, "nphi"_a, "ringstart"_a, "spin"_a, "nthreads"_a);
m.def("adjoint_synthesis", &Pyfull_adjoint_synthesis, "alm"_a, "lmax"_a, "mval"_a, "mstart"_a, "map"_a, "theta"_a, "phi0"_a, "nphi"_a, "ringstart"_a, "spin"_a, "nthreads"_a);
m.def("get_gridweights", &Pyget_gridweights, "type"_a, "nrings"_a);
m.def("alm2leg", &Pyalm2leg, "alm"_a, "theta"_a, "lmax"_a, "spin"_a, "mval"_a=None, "mstart"_a=None, "nthreads"_a=1, "out"_a=None);
......
......@@ -1442,7 +1442,7 @@ template<typename T> void leg2alm( // associated Legendre transform
alm.v(ialm,mstart(mi)+l) = 0;
for (size_t l=lmin; l<=lmax; ++l)
for (size_t ialm=0; ialm<ncomp; ++ialm)
alm.v(ialm,mstart(mi)+l) = almtmp(l,ialm)*T(norm_l[l]);
alm.v(ialm,mstart(mi)+l) = complex<T>(almtmp(l,ialm)*norm_l[l]);
}
}); /* end of parallel region */
}
......@@ -1742,18 +1742,17 @@ void prep_for_analysis2(mav<complex<double>,3> &leg, size_t spin, size_t nthread
}
}
template<typename T> void synthesis(
const mav<complex<T>,2> &alm, // (ncomp, *)
void sanity_checks(
const mav_info<2> &alm, // (ncomp, *)
size_t lmax,
const mav<size_t,1> &mval, // (nm)
const mav<size_t,1> &mstart, // (nm)
mav<T,2> &map, // (ncomp, *)
const mav_info<2> &map, // (ncomp, *)
const mav<double,1> &theta, // (nrings)
const mav<double,1> &phi0, // (nrings)
const mav_info<1> &phi0, // (nrings)
const mav<size_t,1> &nphi, // (nrings)
const mav<size_t,1> &ringstart, // (nrings)
size_t spin,
size_t nthreads)
size_t spin)
{
size_t nm = mval.shape(0);
MR_assert(nm>0, "need at least one m value");
......@@ -1779,9 +1778,25 @@ template<typename T> void synthesis(
size_t ncomp = 1+(spin>0);
MR_assert((alm.shape(0)==ncomp) && (map.shape(0)==ncomp),
"inconsistent number of components");
}
template<typename T> void synthesis(
const mav<complex<T>,2> &alm, // (ncomp, *)
size_t lmax,
const mav<size_t,1> &mval, // (nm)
const mav<size_t,1> &mstart, // (nm)
mav<T,2> &map, // (ncomp, *)
const mav<double,1> &theta, // (nrings)
const mav<double,1> &phi0, // (nrings)
const mav<size_t,1> &nphi, // (nrings)
const mav<size_t,1> &ringstart, // (nrings)
size_t spin,
size_t nthreads)
{
sanity_checks(alm, lmax, mval, mstart, map, theta, phi0, nphi, ringstart, spin);
// just doing standard synthesis now, in the future we can use faster methods
// for some of the theta-equidistant grids here
mav<complex<T>,3> leg({ncomp,nrings,nm});
mav<complex<T>,3> leg({alm.shape(0),theta.shape(0),mval.shape(0)});
alm2leg(alm, leg, theta, mval, mstart, lmax, spin, nthreads, ALM2MAP);
leg2map(leg, map, nphi, ringstart, phi0, nthreads);
}
......@@ -1807,6 +1822,48 @@ template void synthesis(
const mav<size_t,1> &ringstart, // (nrings)
size_t spin,
size_t nthreads);
template<typename T> void adjoint_synthesis(
mav<complex<T>,2> &alm, // (ncomp, *)
size_t lmax,
const mav<size_t,1> &mval, // (nm)
const mav<size_t,1> &mstart, // (nm)
const mav<T,2> &map, // (ncomp, *)
const mav<double,1> &theta, // (nrings)
const mav<double,1> &phi0, // (nrings)
const mav<size_t,1> &nphi, // (nrings)
const mav<size_t,1> &ringstart, // (nrings)
size_t spin,
size_t nthreads)
{
sanity_checks(alm, lmax, mval, mstart, map, theta, phi0, nphi, ringstart, spin);
// just doing standard synthesis now, in the future we can use faster methods
// for some of the theta-equidistant grids here
mav<complex<T>,3> leg({alm.shape(0),theta.shape(0),mval.shape(0)});
map2leg(map, leg, nphi, ringstart, phi0, nthreads);
leg2alm(leg, alm, theta, mval, mstart, lmax, spin, nthreads);
}
template void adjoint_synthesis(
mav<complex<double>,2> &alm, size_t lmax,
const mav<size_t,1> &mval, // (nm)
const mav<size_t,1> &mstart, // (nm)
const mav<double,2> &map, // (ncomp, *)
const mav<double,1> &theta, // (nrings)
const mav<double,1> &phi0, // (nrings)
const mav<size_t,1> &nphi, // (nrings)
const mav<size_t,1> &ringstart, // (nrings)
size_t spin,
size_t nthreads);
template void adjoint_synthesis(
mav<complex<float>,2> &alm, size_t lmax,
const mav<size_t,1> &mval, // (nm)
const mav<size_t,1> &mstart, // (nm)
const mav<float,2> &map, // (ncomp, *)
const mav<double,1> &theta, // (nrings)
const mav<double,1> &phi0, // (nrings)
const mav<size_t,1> &nphi, // (nrings)
const mav<size_t,1> &ringstart, // (nrings)
size_t spin,
size_t nthreads);
#endif
}}
......@@ -462,6 +462,18 @@ template<typename T> void synthesis(
// - check lmax and nrings whether accelerated synthesis makes sense
// - if yes, run synthesis on grid with minimal nrings and upsample
// - otherwise run standard synthesis
template<typename T> void adjoint_synthesis(
mav<complex<T>,2> &alm, // (ncomp, *)
size_t lmax,
const mav<size_t,1> &mval, // (nm)
const mav<size_t,1> &mstart, // (nm)
const mav<T,2> &map, // (ncomp, *)
const mav<double,1> &theta, // (nrings)
const mav<double,1> &phi0, // (nrings)
const mav<size_t,1> &nphi, // (nrings)
const mav<size_t,1> &ringstart, // (nrings)
size_t spin,
size_t nthreads);
template<typename T> void synthesis(const mav<complex<T>,2> &alm, size_t lmax,
mav<T,3> &map, size_t spin, const string &geometry, size_t nthreads)
......
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