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

remove Peano; we don't actually need it

parent a45f5bb3
......@@ -89,133 +89,6 @@ template<typename T> inline T fmodulo (T v1, T v2)
return (tmp==v2) ? T(0) : tmp;
}
//
// Utilities for converting between Cartesian coordinates and Peano index
//
static const uint16_t utab[] = {
#define Z(a) 0x##a##0, 0x##a##1, 0x##a##4, 0x##a##5
#define Y(a) Z(a##0), Z(a##1), Z(a##4), Z(a##5)
#define X(a) Y(a##0), Y(a##1), Y(a##4), Y(a##5)
X(0),X(1),X(4),X(5)
#undef X
#undef Y
#undef Z
};
uint32_t coord2morton2D_32 (uint32_t x, uint32_t y)
{
typedef uint32_t I;
return (I)(utab[x&0xff]) | ((I)(utab[(x>>8)&0xff])<<16)
| ((I)(utab[y&0xff])<<1) | ((I)(utab[(y>>8)&0xff])<<17);
}
static const uint8_t m2p2D_1[4][4] = {
{ 4, 1, 11, 2},{0,15, 5, 6},{10,9,3,12},{14,7,13,8}};
static uint8_t m2p2D_3[4][64];
static const uint8_t p2m2D_1[4][4] = {
{ 4, 1, 3, 10},{0,6,7,13},{15,9,8,2},{11,14,12,5}};
static uint8_t p2m2D_3[4][64];
static int peano2d_done=0;
static void init_peano2d (void)
{
peano2d_done=1;
for (int d=0; d<4; ++d)
for (uint32_t p=0; p<64; ++p)
{
unsigned rot = d;
uint32_t v = p<<26;
uint32_t res = 0;
for (int i=0; i<3; ++i)
{
unsigned tab=m2p2D_1[rot][v>>30];
v<<=2;
res = (res<<2) | (tab&0x3);
rot = tab>>2;
}
m2p2D_3[d][p]=res|(rot<<6);
}
for (int d=0; d<4; ++d)
for (uint32_t p=0; p<64; ++p)
{
unsigned rot = d;
uint32_t v = p<<26;
uint32_t res = 0;
for (int i=0; i<3; ++i)
{
unsigned tab=p2m2D_1[rot][v>>30];
v<<=2;
res = (res<<2) | (tab&0x3);
rot = tab>>2;
}
p2m2D_3[d][p]=res|(rot<<6);
}
}
uint32_t morton2peano2D_32(uint32_t v, int bits)
{
if (!peano2d_done) init_peano2d();
unsigned rot = 0;
uint32_t res = 0;
v<<=32-(bits<<1);
int i=0;
for (; i<bits-2; i+=3)
{
unsigned tab=m2p2D_3[rot][v>>26];
v<<=6;
res = (res<<6) | (tab&0x3fu);
rot = tab>>6;
}
for (; i<bits; ++i)
{
unsigned tab=m2p2D_1[rot][v>>30];
v<<=2;
res = (res<<2) | (tab&0x3);
rot = tab>>2;
}
return res;
}
//
// Utilities for indirect sorting (argsort)
//
template<typename It, typename Comp> class IdxComp__
{
private:
It begin;
Comp comp;
public:
IdxComp__ (It begin_, Comp comp_): begin(begin_), comp(comp_) {}
bool operator() (std::size_t a, std::size_t b) const
{ return comp(*(begin+a),*(begin+b)); }
};
/*! Performs an indirect sort on the supplied iterator range and returns in
\a idx a \a vector containing the indices of the smallest, second smallest,
third smallest, etc. element, according to \a comp. */
template<typename It, typename T2, typename Comp>
inline void buildIndex (It begin, It end, std::vector<T2> &idx, Comp comp)
{
using namespace std;
T2 num=end-begin;
idx.resize(num);
for (T2 i=0; i<num; ++i) idx[i] = i;
stable_sort (idx.begin(),idx.end(),IdxComp__<It,Comp>(begin,comp));
}
/*! Performs an indirect sort on the supplied iterator range and returns in
\a idx a \a vector containing the indices of the smallest, second smallest,
third smallest, etc. element. */
template<typename It, typename T2> inline void buildIndex (It begin, It end,
std::vector<T2> &idx)
{
using namespace std;
typedef typename iterator_traits<It>::value_type T;
buildIndex(begin,end,idx,less<T>());
}
//
// Utilities for Gauss-Legendre quadrature
//
......@@ -424,22 +297,17 @@ template<typename T> class Baselines
private:
vector<UVW<T>> coord;
vector<T> scaling;
vector<bool> flags;
size_t nrows, nchan;
size_t channelbits, channelmask;
public:
Baselines(const pyarr_c<T> &coord_, const pyarr_c<T> &scaling_,
const pyarr_c<bool> &flags_)
Baselines(const pyarr_c<T> &coord_, const pyarr_c<T> &scaling_)
{
myassert(coord_.ndim()==2, "coord array must be 2D");
myassert(coord_.shape(1)==3, "coord.shape[1] must be 3");
myassert(scaling_.ndim()==1, "scaling array must be 1D");
nrows = coord_.shape(0);
nchan = scaling_.shape(0);
myassert(flags_.ndim()==2, "flags array must be 2D");
myassert(size_t(flags_.shape(0))==nrows, "flags.shape[0] must be nrows");
myassert(size_t(flags_.shape(1))==nchan, "flags.shape[1] must be chan");
scaling.resize(nchan);
for (size_t i=0; i<nchan; ++i)
scaling[i] = scaling_.data()[i];
......@@ -447,31 +315,11 @@ template<typename T> class Baselines
auto cood = coord_.data();
for (size_t i=0; i<coord.size(); ++i)
coord[i] = UVW<T>(cood[3*i], cood[3*i+1], cood[3*i+2]);
flags.resize(nrows*nchan);
auto pflags=flags_.data();
for (size_t i=0; i<nrows*nchan; ++i)
flags[i] = (pflags[i]!=0);
channelbits = bits_needed(nchan);
channelmask = (size_t(1)<<channelbits)-1;
auto rowbits = bits_needed(nrows);
myassert(rowbits+channelbits<=8*sizeof(uint32_t), "Ti too small");
}
vector<uint32_t> getIndices(int chbegin, int chend, T wmin,
T wmax) const
{
if (chbegin<0) chbegin=0;
if (chend<0) chend=nchan;
vector<uint32_t> odata;
for (size_t i=0; i<nrows; ++i)
for (int j=chbegin; j<chend; ++j)
if (!flags[i*nchan + j])
{
auto w = coord[i].w*scaling[j];
if ((w>=wmin) && (w<wmax))
odata.push_back((i<<channelbits)+j);
}
return odata;
}
UVW<T> effectiveCoord(uint32_t index) const
{ return coord[index>>channelbits]*scaling[index&channelmask]; }
......@@ -483,6 +331,8 @@ template<typename T> class Baselines
{ return index&channelmask; }
size_t offset(uint32_t index) const
{ return (index>>channelbits)*nchan + (index&channelmask); }
uint32_t Index(size_t irow, size_t ichan) const
{ return (irow<<channelbits) + ichan; }
template<typename T2> pyarr_c<T2> ms2vis(const pyarr_c<T2> &ms_,
const pyarr_c<uint32_t> &idx_) const
......@@ -552,7 +402,6 @@ template<typename T> class GridderConfig
size_t nx_dirty, ny_dirty;
T ucorr, vcorr;
size_t w, nsafe, nu, nv;
size_t peano_level;
vector<T> cfu, cfv;
public:
......@@ -562,7 +411,6 @@ template<typename T> class GridderConfig
ucorr(1./urange), vcorr(1./vrange),
w(get_w(epsilon)), nsafe((w+1)/2),
nu(max(2*nsafe,2*nx_dirty)), nv(max(2*nsafe,2*ny_dirty)),
peano_level(max(1,bits_needed(max(nu, nv)+2*nsafe)-logsquare)),
cfu(nx_dirty), cfv(ny_dirty)
{
myassert((nx_dirty&1)==0, "nx_dirty must be even");
......@@ -585,22 +433,9 @@ template<typename T> class GridderConfig
size_t Nu() const { return nu; }
size_t Nv() const { return nv; }
size_t W() const { return w; }
size_t Nsafe() const { return nsafe; }
T Ucorr() const { return ucorr; }
T Vcorr() const { return vcorr; }
size_t coord2peano(const UVW<T> &coord) const
{
double u=fmodulo(coord.u*ucorr, T(1))*nu,
v=fmodulo(coord.v*vcorr, T(1))*nv;
int iu0 = int(u-w*0.5 + 1 + nu) - nu;
if (iu0+w>nu+nsafe) iu0 = nu+nsafe-w;
iu0+=nsafe;
int iv0 = int(v-w*0.5 + 1 + nv) - nv;
if (iv0+w>nv+nsafe) iv0 = nv+nsafe-w;
iv0+=nsafe;
iu0>>=logsquare;
iv0>>=logsquare;
return morton2peano2D_32(coord2morton2D_32(iu0,iv0),peano_level);
}
pyarr_c<T> grid2dirty(const pyarr_c<T> &grid) const
{
myassert(grid.ndim()==2, "grid must be a 2D array");
......@@ -918,19 +753,76 @@ template<typename T> pyarr_c<complex<T>> grid2vis(const Baselines<T> &baselines,
}
template<typename T> pyarr_c<uint32_t> getIndices(const Baselines<T> &baselines,
const GridderConfig<T> &gconf, int chbegin, int chend, T wmin, T wmax)
const GridderConfig<T> &gconf, const pyarr_c<bool> &flags_, int chbegin,
int chend, T wmin, T wmax)
{
auto idx = baselines.getIndices(chbegin, chend, wmin, wmax);
vector<size_t> peano(idx.size());
for (size_t i=0; i<peano.size(); ++i)
peano[i] = gconf.coord2peano(baselines.effectiveCoord(idx[i]));
vector<size_t> newind;
buildIndex(peano.begin(), peano.end(), newind);
peano=vector<size_t>(); // deallocate
auto res = makearray<uint32_t>({idx.size()});
size_t nrow=baselines.Nrows(),
nchan=baselines.Nchannels(),
nu=gconf.Nu(),
nv=gconf.Nv(),
nsafe=gconf.Nsafe(),
w=gconf.W();
T ucorr=gconf.Ucorr(),
vcorr=gconf.Vcorr();
if (chbegin<0) chbegin=0;
if (chend<0) chend=nchan;
myassert(flags_.ndim()==2, "flags must be 2D");
myassert(size_t(flags_.shape(0))==nrow, "bad flags dimension");
myassert(size_t(flags_.shape(1))==nchan, "bad flags dimension");
auto flags = flags_.data();
constexpr int side=1<<logsquare;
size_t nbu = (gconf.Nu()+1+side-1) >> logsquare,
nbv = (gconf.Nv()+1+side-1) >> logsquare;
vector<uint32_t> bincnt(nbu*nbv, 0);
for (size_t irow=0; irow<nrow; ++irow)
for (int ichan=chbegin; ichan<chend; ++ichan)
if (!flags[irow*nchan + ichan])
{
auto idx = baselines.Index(irow, ichan);
auto uvw = baselines.effectiveCoord(idx);
if ((uvw.w>=wmin) && (uvw.w<wmax))
{
auto u=fmodulo(uvw.u*ucorr, T(1))*nu,
v=fmodulo(uvw.v*vcorr, T(1))*nv;
int iu0 = int(u-w*0.5 + 1 + nu) - nu;
if (iu0+w>nu+nsafe) iu0 = nu+nsafe-w;
iu0+=nsafe;
int iv0 = int(v-w*0.5 + 1 + nv) - nv;
if (iv0+w>nv+nsafe) iv0 = nv+nsafe-w;
iv0+=nsafe;
iu0>>=logsquare;
iv0>>=logsquare;
++bincnt[nbv*iu0 + iv0];
}
}
vector<uint32_t> acc(bincnt.size()+1);
acc[0] = 0;
for (size_t i=0; i<bincnt.size(); ++i)
acc[i+1] = acc[i] + bincnt[i];
auto res = makearray<uint32_t>({acc.back()});
auto iout = res.mutable_data();
for (size_t i=0; i<idx.size(); ++i)
iout[i] = idx[newind[i]];
for (size_t irow=0; irow<nrow; ++irow)
for (int ichan=chbegin; ichan<chend; ++ichan)
if (!flags[irow*nchan + ichan])
{
auto idx = baselines.Index(irow, ichan);
auto uvw = baselines.effectiveCoord(idx);
if ((uvw.w>=wmin) && (uvw.w<wmax))
{
auto u=fmodulo(uvw.u*ucorr, T(1))*nu,
v=fmodulo(uvw.v*vcorr, T(1))*nv;
int iu0 = int(u-w*0.5 + 1 + nu) - nu;
if (iu0+w>nu+nsafe) iu0 = nu+nsafe-w;
iu0+=nsafe;
int iv0 = int(v-w*0.5 + 1 + nv) - nv;
if (iv0+w>nv+nsafe) iv0 = nv+nsafe-w;
iv0+=nsafe;
iu0>>=logsquare;
iv0>>=logsquare;
iout[acc[nbv*iu0 + iv0]] = idx;
++acc[nbv*iu0 + iv0];
}
}
return res;
}
......@@ -943,8 +835,6 @@ coord: np.array((nrows, 3), dtype=np.float)
u, v and w coordinates for each row
scaling: np.array((nchannels,), dtype=np.float)
scaling factor for u, v, w for each individual channel
flags: np.array((nrows, nchannels), dtype=np.bool)
"True" indicates that the value should not be used
)""";
const char *BL_ms2vis_DS = R"""(
......@@ -1051,6 +941,8 @@ baselines: Baselines
gconf: GridderConf
the GridderConf object to be used with the returned indices.
(used to optimize the ordering of the indices)
flags: np.array((nrows, nchannels), dtype=np.bool)
"True" indicates that the value should not be used
chbegin: int
first channel to use (-1: start with the first available channel)
chend: int
......@@ -1117,8 +1009,8 @@ PYBIND11_MODULE(nifty_gridder, m)
using namespace pybind11::literals;
py::class_<Baselines<double>> (m, "Baselines", Baselines_DS)
.def(py::init<const pyarr_c<double> &, const pyarr_c<double> &,
const pyarr_c<bool> &>(), "coord"_a, "scaling"_a, "flags"_a)
.def(py::init<const pyarr_c<double> &, const pyarr_c<double> &>(),
"coord"_a, "scaling"_a)
.def ("Nrows",&Baselines<double>::Nrows)
.def ("Nchannels",&Baselines<double>::Nchannels)
.def ("ms2vis",&Baselines<double>::ms2vis<complex<double>>, BL_ms2vis_DS, "ms"_a, "idx"_a)
......@@ -1133,14 +1025,14 @@ PYBIND11_MODULE(nifty_gridder, m)
.def("Nv", &GridderConfig<double>::Nv)
.def("grid2dirty", &GridderConfig<double>::grid2dirty, grid2dirty_DS, "grid"_a)
.def("dirty2grid", &GridderConfig<double>::dirty2grid, dirty2grid_DS, "dirty"_a);
m.def("getIndices", getIndices<double>, getIndices_DS, "baselines"_a,
"gconf"_a, "chbegin"_a=-1, "chend"_a=-1, "wmin"_a=-1e30, "wmax"_a=1e30);
m.def("getIndices", getIndices<double>, getIndices_DS, "baselines"_a, "gconf"_a,
"flags"_a, "chbegin"_a=-1, "chend"_a=-1, "wmin"_a=-1e30, "wmax"_a=1e30);
m.def("vis2grid",&vis2grid<double>, vis2grid_DS, "baselines"_a, "gconf"_a, "idx"_a, "vis"_a);
m.def("grid2vis",&grid2vis<double>, grid2vis_DS, "baselines"_a, "gconf"_a, "idx"_a, "grid"_a);
py::class_<Baselines<float>> (m, "Baselines_f", Baselines_DS)
.def(py::init<const pyarr_c<float> &, const pyarr_c<float> &,
const pyarr_c<bool> &>(), "coord"_a, "scaling"_a, "flags"_a)
.def(py::init<const pyarr_c<float> &, const pyarr_c<float> &>(),
"coord"_a, "scaling"_a)
.def ("Nrows",&Baselines<float>::Nrows)
.def ("Nchannels",&Baselines<float>::Nchannels)
.def ("ms2vis",&Baselines<float>::ms2vis<complex<float>>, BL_ms2vis_DS, "ms"_a, "idx"_a)
......@@ -1154,8 +1046,8 @@ PYBIND11_MODULE(nifty_gridder, m)
.def("Nv", &GridderConfig<float>::Nv)
.def("grid2dirty", &GridderConfig<float>::grid2dirty, "grid"_a)
.def("dirty2grid", &GridderConfig<float>::dirty2grid, "dirty"_a);
m.def("getIndices_f", getIndices<float>, "baselines"_a, "gconf"_a, "chbegin"_a=-1,
"chend"_a=-1, "wmin"_a=-1e30, "wmax"_a=1e30);
m.def("getIndices_f", getIndices<float>, "baselines"_a, "gconf"_a,
"flags"_a, "chbegin"_a=-1, "chend"_a=-1, "wmin"_a=-1e30, "wmax"_a=1e30);
m.def("vis2grid_f",&vis2grid<float>, "baselines"_a, "gconf"_a, "idx"_a, "vis"_a);
m.def("grid2vis_f",&grid2vis<float>, "baselines"_a, "gconf"_a, "idx"_a, "grid"_a);
}
Supports Markdown
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