Commit 3eae85ce authored by Martin Reinecke's avatar Martin Reinecke

fully generic synthesis

parent 6a00d0f0
Pipeline #97802 passed with stages
in 22 minutes and 13 seconds
......@@ -170,6 +170,37 @@ void Pyresample_theta(const py::array &legi_, bool npi, bool spi,
resample_theta(legi, npi, spi, lego, npo, spo, spin, nthreads);
}
template<typename T> py::array full_synthesis2(const py::array &alm_,
size_t lmax, const py::array &mval_, const py::array &mstart_,
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_, false);
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_, true);
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);
synthesis(alm, lmax, mval, mstart, map, theta, phi0, nphi, ringstart, spin, nthreads);
return map_;
}
py::array Pyfull_synthesis(const py::array &alm_, size_t lmax,
const py::array &mval_, const py::array &mstart_,
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_synthesis2<float>(alm_, lmax, mval_, mstart_, map_, theta_,
phi0_, nphi_, ringstart_, spin, nthreads);
else if (isPyarr<complex<double>>(alm_))
return full_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>;
......@@ -324,6 +355,8 @@ void add_sht(py::module_ &msup)
#if 1
// 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("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);
m.def("leg2alm", &Pyleg2alm, "leg"_a, "theta"_a, "lmax"_a, "spin"_a, "mval"_a=None, "mstart"_a=None, "nthreads"_a=1, "out"_a=None);
......
......@@ -1389,7 +1389,7 @@ template<typename T> void alm2leg( // associated Legendre transform
for (size_t l=m; l<lmin; ++l)
almtmp.v(l,ialm) = 0;
for (size_t l=lmin; l<=lmax; ++l)
almtmp.v(l,ialm) = alm(ialm,mstart(mi)+l)*norm_l[l];
almtmp.v(l,ialm) = alm(ialm,mstart(mi)+l)*T(norm_l[l]);
almtmp.v(lmax+1,ialm) = 0;
}
gen.prepare(m);
......@@ -1442,14 +1442,14 @@ 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)*norm_l[l];
alm.v(ialm,mstart(mi)+l) = almtmp(l,ialm)*T(norm_l[l]);
}
}); /* end of parallel region */
}
template<typename T> void leg2map( // FFT
const mav<complex<T>,3> &leg, // (ncomp, nrings, mmax+1)
mav<complex<T>,2> &map, // (ncomp, pix)
mav<T,2> &map, // (ncomp, pix)
const mav<size_t,1> &nphi, // (nrings)
const mav<size_t,1> &offset, // (nrings)
const mav<double,1> &phi0, // (nrings)
......@@ -1477,14 +1477,14 @@ template<typename T> void leg2map( // FFT
auto ltmp = subarray<1>(leg, {icomp, ith, 0}, {0, 0, MAXIDX});
helper.phase2ring (nphi(ith),phi0(ith),ringtmp,mmax,ltmp);
for (size_t i=0; i<nphi(ith); ++i)
map.v(icomp,offset(ith)+i) = ringtmp(i);
map.v(icomp,offset(ith)+i) = T(ringtmp(i+1));
}
}
}); /* end of parallel region */
}
template<typename T> void map2leg( // FFT
const mav<complex<T>,2> &map, // (ncomp, pix)
const mav<T,2> &map, // (ncomp, pix)
mav<complex<T>,3> &leg, // (ncomp, nrings, mmax+1)
const mav<size_t,1> &nphi, // (nrings)
const mav<size_t,1> &offset, // (nrings)
......@@ -1741,6 +1741,72 @@ 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, *)
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)
{
size_t nm = mval.shape(0);
MR_assert(nm>0, "need at least one m value");
MR_assert(nm==mstart.shape(0), "nm mismatch");
// check that mval are a gapless sequence starting at 0
size_t mmax = nm-1;
{
vector<bool> m_present(mmax+1, false);
for (size_t im=0; im<=mmax; ++im)
{
MR_assert(mval(im)<=mmax, "m value too large");
MR_assert(!m_present[mval(im)], "m value supplied more than once");
m_present[mval(im)] = true;
}
}
MR_assert(lmax>=mmax, "lmax must be >= mmax");
size_t nrings = theta.shape(0);
MR_assert(nrings>0, "need at least one ring");
MR_assert((phi0.shape(0)==nrings) &&
(nphi.shape(0)==nrings) &&
(ringstart.shape(0)==nrings),
"inconsistency in the number of rings");
size_t ncomp = 1+(spin>0);
MR_assert((alm.shape(0)==ncomp) && (map.shape(0)==ncomp),
"inconsistent number of components");
// 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});
alm2leg(alm, leg, theta, mval, mstart, lmax, spin, nthreads, ALM2MAP);
leg2map(leg, map, nphi, ringstart, phi0, nthreads);
}
template void synthesis(
const mav<complex<double>,2> &alm, size_t lmax,
const mav<size_t,1> &mval, // (nm)
const mav<size_t,1> &mstart, // (nm)
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 synthesis(
const mav<complex<float>,2> &alm, size_t lmax,
const mav<size_t,1> &mval, // (nm)
const mav<size_t,1> &mstart, // (nm)
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
}}
......@@ -297,8 +297,8 @@ struct ringhelper
length=nph;
}
}
DUCC0_NOINLINE void phase2ring (size_t nph, double phi0,
mav<double,1> &data, size_t mmax, const mav<dcmplx,1> &phase)
template<typename T> DUCC0_NOINLINE void phase2ring (size_t nph,
double phi0, mav<double,1> &data, size_t mmax, const mav<complex<T>,1> &phase)
{
update (nph, mmax, phi0);
......@@ -313,7 +313,7 @@ struct ringhelper
else
for (size_t m=0; m<=mmax; ++m)
{
dcmplx tmp = phase(m)*shiftarr[m];
dcmplx tmp = dcmplx(phase(m))*shiftarr[m];
data.v(2*m)=tmp.real();
data.v(2*m+1)=tmp.imag();
}
......@@ -345,8 +345,8 @@ struct ringhelper
data.v(1)=data(0);
plan->exec(&(data.v(1)), 1., false);
}
DUCC0_NOINLINE void ring2phase (size_t nph, double phi0,
mav<double,1> &data, size_t mmax, mav<dcmplx,1> &phase)
template<typename T> DUCC0_NOINLINE void ring2phase (size_t nph, double phi0,
mav<double,1> &data, size_t mmax, mav<complex<T>,1> &phase)
{
update (nph, mmax, -phi0);
......@@ -358,10 +358,10 @@ struct ringhelper
{
if (norot)
for (size_t m=0; m<=mmax; ++m)
phase.v(m) = dcmplx(data(2*m), data(2*m+1));
phase.v(m) = complex<T>(T(data(2*m)), T(data(2*m+1)));
else
for (size_t m=0; m<=mmax; ++m)
phase.v(m) = dcmplx(data(2*m), data(2*m+1)) * shiftarr[m];
phase.v(m) = complex<T>(dcmplx(data(2*m), data(2*m+1)) * shiftarr[m]);
}
else
{
......@@ -374,7 +374,7 @@ struct ringhelper
val = dcmplx(data(2*(nph-idx)), -data(2*(nph-idx)+1));
if (!norot)
val *= shiftarr[m];
phase.v(m)=val;
phase.v(m)=complex<T>(val);
}
}
}
......@@ -404,14 +404,14 @@ template<typename T> void alm2leg( // associated Legendre transform
template<typename T> void leg2map( // FFT
const mav<complex<T>,3> &leg, // (ncomp, nrings, mmax+1)
mav<complex<T>,2> &map, // (ncomp, pix)
mav<T,2> &map, // (ncomp, pix)
const mav<size_t,1> &nphi, // (nrings)
const mav<size_t,1> &offset, // (nrings)
const mav<double,1> &phi0, // (nrings)
size_t nthreads);
template<typename T> void map2leg( // FFT
const mav<complex<T>,2> &map, // (ncomp, pix)
const mav<T,2> &map, // (ncomp, pix)
mav<complex<T>,3> &leg, // (ncomp, nrings, mmax+1)
const mav<size_t,1> &nphi, // (nrings)
const mav<size_t,1> &offset, // (nrings)
......@@ -435,6 +435,34 @@ void resample_theta(const mav<complex<double>,2> &legi, bool npi, bool spi,
#endif
// fully general map synthesis
// conditions:
// - ncomp = 1+(spin>0)
// - all mval together must form a gapless sequence starting from 0
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);
// - check that mval are a gapless sequence starting at 0
// - check that lmax>=mmax
// - (check that mstart has consistent values?)
// - check that theta are in [0;pi]
// - (check that ringstart and nphi are consistent)
// - check that ncomp and spin are consistent
// - is theta grid equidistant and (CC, F!, MW)?
// - if no, run standard 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 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