Skip to content

Commit

Permalink
Double precision and AIPS style-gridding (#809)
Browse files Browse the repository at this point in the history 
* Add double precision gridding mode switching

* Ensure Wplanes==1 gives real valued filters
  • Loading branch information
@bennahugo
bennahugo authored Feb 11, 2022
1 parent 88d5713 commit 1e9186e
Show file tree
Hide file tree
Showing 6 changed files with 198 additions and 163 deletions.
54 changes: 30 additions & 24 deletions DDFacet/Gridder/GridderSmearPols.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,8 +68,9 @@ namespace DDF {
}
}
}

void pyGridderWPol(py::array_t<std::complex<float>, py::array::c_style>& np_grid,

template<typename gridtype>
void pyGridderWPol(py::array_t<std::complex<gridtype>, py::array::c_style>& np_grid,
const py::array_t<std::complex<float>, py::array::c_style>& vis,
const py::array_t<double, py::array::c_style>& uvw,
const py::array_t<bool, py::array::c_style>& flags,
Expand All @@ -89,8 +90,7 @@ namespace DDF {
const py::list& LSmearing,
const py::array_t<int32_t, py::array::c_style>& np_ChanMapping,
const py::array_t<uint16_t, py::array::c_style>& LDataCorrFormat,
const py::array_t<uint16_t, py::array::c_style>& LExpectedOutStokes
)
const py::array_t<uint16_t, py::array::c_style>& LExpectedOutStokes)
{
using svec = vector<string>;
const svec stokeslookup = {"undef","I","Q","U","V","RR","RL","LR","LL","XX","XY","YX","YY"};
Expand All @@ -111,10 +111,13 @@ namespace DDF {
throw std::invalid_argument("Only accepts I,Q,U,V as polarization output type");
expstokes[i] = stokeslookup[polid];
}
#define callgridder(stokesgrid, nVisPol) \
#define callgridder32(stokesgrid, nVisPol) \
{\
gridder::gridder<readcorr, mulaccum, stokesgrid>(np_grid, vis, uvw, flags, weights, sumwt, bool(dopsf), Lcfs, LcfsConj, WInfos, increment, freqs, Lmaps, LJones, SmearMapping, Sparsification, LOptimisation,LSmearing,np_ChanMapping, expstokes); \
done=true;\
gridder::gridder<readcorr, mulaccum, stokesgrid, gridtype>(np_grid, vis, uvw, flags, weights, sumwt,\
bool(dopsf), Lcfs, LcfsConj, WInfos, increment,\
freqs, Lmaps, LJones, SmearMapping, Sparsification,\
LOptimisation,LSmearing,np_ChanMapping, expstokes); \
done=true;\
}
using namespace DDF::gridder::policies;
bool done=false;
Expand All @@ -123,15 +126,15 @@ namespace DDF {
#define readcorr gridder::policies::Read_4Corr
#define mulaccum gridder::policies::Mulaccum_4Corr
if (expstokes==svec{"I"})
callgridder(I_from_XXXYYXYY, 1)
callgridder32(I_from_XXXYYXYY, 1)
else if (expstokes==svec{"I", "Q"})
callgridder(IQ_from_XXXYYXYY, 2)
callgridder32(IQ_from_XXXYYXYY, 2)
else if (expstokes==svec{"I", "V"})
callgridder(IV_from_XXXYYXYY, 2)
callgridder32(IV_from_XXXYYXYY, 2)
else if (expstokes==svec{"Q", "U"})
callgridder(QU_from_XXXYYXYY, 2)
callgridder32(QU_from_XXXYYXYY, 2)
else if (expstokes==svec{"I", "Q", "U", "V"})
callgridder(IQUV_from_XXXYYXYY, 4)
callgridder32(IQUV_from_XXXYYXYY, 4)
#undef readcorr
#undef mulaccum
}
Expand All @@ -140,15 +143,15 @@ namespace DDF {
#define readcorr gridder::policies::Read_4Corr
#define mulaccum gridder::policies::Mulaccum_4Corr
if (expstokes==svec{"I"})
callgridder(I_from_RRRLLRLL, 1)
callgridder32(I_from_RRRLLRLL, 1)
else if (expstokes==svec{"I", "Q"})
callgridder(IQ_from_RRRLLRLL, 2)
callgridder32(IQ_from_RRRLLRLL, 2)
else if (expstokes==svec{"I", "V"})
callgridder(IV_from_RRRLLRLL, 2)
callgridder32(IV_from_RRRLLRLL, 2)
else if (expstokes==svec{"Q", "U"})
callgridder(QU_from_RRRLLRLL, 2)
callgridder32(QU_from_RRRLLRLL, 2)
else if (expstokes==svec{"I", "Q", "U", "V"})
callgridder(IQUV_from_RRRLLRLL, 4)
callgridder32(IQUV_from_RRRLLRLL, 4)
#undef readcorr
#undef mulaccum
}
Expand All @@ -157,9 +160,9 @@ namespace DDF {
#define readcorr gridder::policies::Read_2Corr_Pad
#define mulaccum gridder::policies::Mulaccum_2Corr_Unpad
if (expstokes==svec{"I"})
callgridder(I_from_XXYY, 1)
callgridder32(I_from_XXYY, 1)
else if (expstokes==svec{"I", "Q"})
callgridder(IQ_from_XXYY, 2)
callgridder32(IQ_from_XXYY, 2)
#undef readcorr
#undef mulaccum
}
Expand All @@ -168,9 +171,9 @@ namespace DDF {
#define readcorr gridder::policies::Read_2Corr_Pad
#define mulaccum gridder::policies::Mulaccum_2Corr_Unpad
if (expstokes==svec{"I"})
callgridder(I_from_RRLL, 1)
callgridder32(I_from_RRLL, 1)
else if (expstokes==svec{"I", "V"})
callgridder(IV_from_RRLL, 2)
callgridder32(IV_from_RRLL, 2)
#undef readcorr
#undef mulaccum
}
Expand All @@ -179,14 +182,15 @@ namespace DDF {
#define readcorr gridder::policies::Read_1Corr_Pad
#define mulaccum gridder::policies::Mulaccum_1Corr_Unpad
if (expstokes==svec{"I"})
callgridder(I_from_I, 1)
callgridder32(I_from_I, 1)
#undef readcorr
#undef mulaccum
}
if (!done)
throw std::invalid_argument("Cannot convert input correlations to desired output Stokes parameters.");
}


py::array_t<std::complex<float>, py::array::c_style> & pyDeGridderWPol(
const py::array_t<std::complex<float>, py::array::c_style>& np_grid,
py::array_t<std::complex<float>, py::array::c_style>& np_vis,
Expand Down Expand Up @@ -263,8 +267,10 @@ namespace DDF {
&pyAccumulateWeightsOntoGrid);
m.def("pyAccumulateWeightsOntoGridNoSem",
&pyAccumulateWeightsOntoGridNoSem);
m.def("pyGridderWPol",
&pyGridderWPol);
m.def("pyGridderWPol32",
&pyGridderWPol<float>);
m.def("pyGridderWPol64",
&pyGridderWPol<double>);
m.def("pyDeGridderWPol",
&pyDeGridderWPol,
py::return_value_policy::take_ownership);
Expand Down
9 changes: 5 additions & 4 deletions DDFacet/Gridder/gridder.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,8 +75,9 @@ namespace DDF {
Vis[0] += VisMeas[0]*Weight;
}
}
template<policies::ReadCorrType readcorr, policies::MulaccumType mulaccum, policies::StokesGridType stokesgrid>
void gridder(py::array_t<std::complex<float>, py::array::c_style>& grid,
template<policies::ReadCorrType readcorr, policies::MulaccumType mulaccum, policies::StokesGridType stokesgrid,
typename accum_grid_type>
void gridder(py::array_t<std::complex<accum_grid_type>, py::array::c_style>& grid,
const py::array_t<std::complex<float>, py::array::c_style>& vis,
const py::array_t<double, py::array::c_style>& uvw,
const py::array_t<bool, py::array::c_style>& flags,
Expand Down Expand Up @@ -118,7 +119,7 @@ namespace DDF {
const int nGridY = int(grid.shape(2));
const int nGridPol = int(grid.shape(1));
const int nGridChan = int(grid.shape(0));
fcmplx *griddata = grid.mutable_data(0);
std::complex<accum_grid_type> *griddata = grid.mutable_data(0);

const fcmplx *visdata = vis.data(0);

Expand Down Expand Up @@ -355,7 +356,7 @@ namespace DDF {
const size_t goff = size_t((gridChan*nGridPol + ipol) * nGridX*nGridY);
const dcmplx VisVal =stokes_vis[ipol];
const fcmplx* __restrict__ cf0 = cfsdata + cfoff;
fcmplx* __restrict__ gridPtr = griddata + goff + (locy-supy)*nGridX + locx;
std::complex<accum_grid_type>* __restrict__ gridPtr = griddata + goff + (locy-supy)*nGridX + locx;
for (int sy=-supy; sy<=supy; ++sy, gridPtr+=nGridX)
for (int sx=-supx; sx<=supx; ++sx)
gridPtr[sx] += VisVal * dcmplx(*cf0++);
Expand Down
61 changes: 33 additions & 28 deletions DDFacet/Imager/ClassDDEGridMachine.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -366,7 +366,6 @@ def __init__(self,
elif Precision == "D":
self.dtype = np.complex128

self.dtype = np.complex64
T.timeit("0")
Padding = GD["Facets"]["Padding"]
self.NonPaddedNpix, Npix = EstimateNpix(Npix, Padding)
Expand Down Expand Up @@ -784,32 +783,36 @@ def put(self, times, uvw, visIn, flag, A0A1, W=None,
ChanEquidistant,
self.SkyType,
self.PolModeID]
_pyGridderSmear.pyGridderWPol(Grid,
vis,
uvw,
flag,
W,
SumWeigths,
DoPSF,
self.WTerm.Wplanes,
self.WTerm.WplanesConj,
np.array([self.WTerm.RefWave,
self.WTerm.wmax,
len(self.WTerm.Wplanes),
self.WTerm.OverS],
dtype=np.float64),
self.incr.astype(np.float64),
freqs,
[self.PolMap,
FacetInfos],
ParamJonesList,
self._bda_grid,
sparsification if sparsification is not None else np.array([], dtype=np.bool),
OptimisationInfos,
self.LSmear,
np.int32(ChanMapping),
np.array(self.DataCorrelationFormat).astype(np.uint16),
np.array(self.ExpectedOutputStokes).astype(np.uint16))
if self.dtype == np.complex64:
gridfun = _pyGridderSmear.pyGridderWPol32
else:
gridfun = _pyGridderSmear.pyGridderWPol64
gridfun(Grid,
vis,
uvw,
flag,
W,
SumWeigths,
DoPSF,
self.WTerm.Wplanes,
self.WTerm.WplanesConj,
np.array([self.WTerm.RefWave,
self.WTerm.wmax,
len(self.WTerm.Wplanes),
self.WTerm.OverS],
dtype=np.float64),
self.incr.astype(np.float64),
freqs,
[self.PolMap,
FacetInfos],
ParamJonesList,
self._bda_grid,
sparsification if sparsification is not None else np.array([], dtype=np.bool),
OptimisationInfos,
self.LSmear,
np.int32(ChanMapping),
np.array(self.DataCorrelationFormat).astype(np.uint16),
np.array(self.ExpectedOutputStokes).astype(np.uint16))

T.timeit("gridder")
T.timeit("grid %d" % self.IDFacet)
Expand All @@ -819,6 +822,8 @@ def put(self, times, uvw, visIn, flag, A0A1, W=None,
ChanEquidistant,
self.SkyType,
self.PolModeID]
if self.dtype != np.complex64:
raise RuntimeError("Classic BDA does not support double accumulation gridding")
_pyGridderSmearClassic.pyGridderWPol(Grid,
vis,
uvw,
Expand Down Expand Up @@ -863,7 +868,7 @@ def CheckTypes(
A1=None,
Jones=None):
if not isinstance(Grid, type(None)):
if not(Grid.dtype == np.complex64):
if not(Grid.dtype == self.dtype):
raise NameError(
'Grid.dtype %s %s' %
(str(
Expand Down
6 changes: 4 additions & 2 deletions DDFacet/Imager/ClassFacetMachine.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,10 +94,13 @@ def __init__(self,
self.stitchedType = np.float32 # cleaning requires float32
elif Precision == "D":
self.dtype = np.complex128
self.CType = np.complex64
self.CType = np.complex128
self.FType = np.float64
self.stitchedType = np.float32 # cleaning requires float32

if int(GD["CF"]["Nw"]) <= 1:
print("Disabling W-projection. Enabling AIPS-style faceting", file=log)

self.DoDDE = False
if Sols is not None:
self.setSols(Sols)
Expand Down Expand Up @@ -332,7 +335,6 @@ def AppendFacet(self, iFacet, l0, m0, diam):
self.DicoImager[iFacet]["lmSol"] = lSol[iSol], mSol[iSol]
self.DicoImager[iFacet]["radecSol"] = raSol[iSol], decSol[iSol]
self.DicoImager[iFacet]["iSol"] = iSol

DicoConfigGM = {"NPix": NpixFacet,
"Cell": self.GD["Image"]["Cell"],
"ChanFreq": self.ChanFreq,
Expand Down
Loading

0 comments on commit 1e9186e

Please sign in to comment.