diff --git a/DDFacet/Gridder/GridderSmearPols.cc b/DDFacet/Gridder/GridderSmearPols.cc index fc65aba6..672d3405 100644 --- a/DDFacet/Gridder/GridderSmearPols.cc +++ b/DDFacet/Gridder/GridderSmearPols.cc @@ -68,8 +68,9 @@ namespace DDF { } } } - - void pyGridderWPol(py::array_t, py::array::c_style>& np_grid, + + template + void pyGridderWPol(py::array_t, py::array::c_style>& np_grid, const py::array_t, py::array::c_style>& vis, const py::array_t& uvw, const py::array_t& flags, @@ -89,8 +90,7 @@ namespace DDF { const py::list& LSmearing, const py::array_t& np_ChanMapping, const py::array_t& LDataCorrFormat, - const py::array_t& LExpectedOutStokes - ) + const py::array_t& LExpectedOutStokes) { using svec = vector; const svec stokeslookup = {"undef","I","Q","U","V","RR","RL","LR","LL","XX","XY","YX","YY"}; @@ -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(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(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; @@ -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 } @@ -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 } @@ -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 } @@ -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 } @@ -179,7 +182,7 @@ 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 } @@ -187,6 +190,7 @@ namespace DDF { throw std::invalid_argument("Cannot convert input correlations to desired output Stokes parameters."); } + py::array_t, py::array::c_style> & pyDeGridderWPol( const py::array_t, py::array::c_style>& np_grid, py::array_t, py::array::c_style>& np_vis, @@ -263,8 +267,10 @@ namespace DDF { &pyAccumulateWeightsOntoGrid); m.def("pyAccumulateWeightsOntoGridNoSem", &pyAccumulateWeightsOntoGridNoSem); - m.def("pyGridderWPol", - &pyGridderWPol); + m.def("pyGridderWPol32", + &pyGridderWPol); + m.def("pyGridderWPol64", + &pyGridderWPol); m.def("pyDeGridderWPol", &pyDeGridderWPol, py::return_value_policy::take_ownership); diff --git a/DDFacet/Gridder/gridder.h b/DDFacet/Gridder/gridder.h index 5abd3a15..e85ae170 100644 --- a/DDFacet/Gridder/gridder.h +++ b/DDFacet/Gridder/gridder.h @@ -75,8 +75,9 @@ namespace DDF { Vis[0] += VisMeas[0]*Weight; } } - template - void gridder(py::array_t, py::array::c_style>& grid, + template + void gridder(py::array_t, py::array::c_style>& grid, const py::array_t, py::array::c_style>& vis, const py::array_t& uvw, const py::array_t& flags, @@ -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 *griddata = grid.mutable_data(0); const fcmplx *visdata = vis.data(0); @@ -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* __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++); diff --git a/DDFacet/Imager/ClassDDEGridMachine.py b/DDFacet/Imager/ClassDDEGridMachine.py index 035e1491..d46e4db9 100644 --- a/DDFacet/Imager/ClassDDEGridMachine.py +++ b/DDFacet/Imager/ClassDDEGridMachine.py @@ -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) @@ -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) @@ -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, @@ -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( diff --git a/DDFacet/Imager/ClassFacetMachine.py b/DDFacet/Imager/ClassFacetMachine.py index 8e74fc8e..5399e1ae 100644 --- a/DDFacet/Imager/ClassFacetMachine.py +++ b/DDFacet/Imager/ClassFacetMachine.py @@ -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) @@ -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, diff --git a/DDFacet/Imager/ModCF.py b/DDFacet/Imager/ModCF.py index 2576d5cd..1b82eb95 100644 --- a/DDFacet/Imager/ModCF.py +++ b/DDFacet/Imager/ModCF.py @@ -292,7 +292,6 @@ def Give_dn(l0, m0, rad=1., order=4): return Cl, Cm, C.flatten() - class ClassWTermModified(): def __init__(self, Cell=10, Sup=15, Nw=11, wmax=30000, Npix=101, Freqs=np.array([100.e6]), OverS=11, lmShift=None, mode="compute", @@ -315,8 +314,8 @@ def __init__(self, Cell=10, Sup=15, Nw=11, wmax=30000, Npix=101, Freqs=np.array( "dict" to load CFs from store_dict IDFacet: """ - self.Nw = int(Nw) + self.Nw = max(1, Nw) self.Cell = Cell self.Sup = Sup self.Nw = Nw @@ -334,6 +333,8 @@ def __init__(self, Cell=10, Sup=15, Nw=11, wmax=30000, Npix=101, Freqs=np.array( if compute_cf: cf_dict["wmax"] = self.wmax = wmax self.InitSphe() + # enable only AIPS faceting + # we are using linspace so linspace(...,...,0) will result int a real-valued AA filter self.InitW() dS = np.float32 cf_dict["Sphe"] = dS(self.ifzfCF.real) @@ -444,115 +445,19 @@ def InitW(self): # print self.IDFacet,l0,m0,self.Cv,self.Cu # print "done FIT" - - for i in range(Nw): - #print>>log, "%i/%i"%(i,Nw) - if not(Sups[i] % 2): - Sups[i] += 1 - dummy, dymmy, ThisSphe = self.SpheM.MakeSphe(Sups[i]) - wl = w[i]/waveMin - - # ############## - # l,m=np.mgrid[-lrad:lrad:Npix*1j,-lrad:lrad:Npix*1j] - # n_1=np.sqrt(1.-l**2-m**2)-1 - # WTrue=np.exp(-2.*1j*np.pi*wl*(n_1))*self.ifzfCF - # ############## - - DX = 2*lrad/Sups[i] - l, m = np.mgrid[-lrad+DX/2:lrad-DX/2:Sups[i] - * 1j, -lrad+DX/2:lrad-DX/2:Sups[i]*1j] - # l,m=np.mgrid[-lrad:lrad:Sups[i]*1j,-lrad:lrad:Sups[i]*1j] - # n_1=np.sqrt(1.-l**2-m**2)-1 - n_1 = ModFitPoly2D.polyval2d(l, m, CoefPoly) - # n_1=n_1.T[::-1,:] - T.timeit("3a") - - # stop - # n_1=np.sqrt(1.-(l-l0)**2-(m-m0)**2)-1 - # n_1=(1./np.sqrt(1.-l0**2-m0**2))*(l0*l+m0*m) - W = np.exp(-2.*1j*np.pi*wl*(n_1)) - #import pylab - - # pylab.clf() - # pylab.imshow(np.angle(W),interpolation="nearest",extent=(l.min(),l.max(),m.min(),m.max()),vmin=-np.pi,vmax=np.pi) - # pylab.draw() - # pylab.show(False) - # pylab.pause(0.1) - # T.timeit("3b") - - # #### - # W.fill(1.) - # #### - W *= np.abs(ThisSphe) - - # # #################### - # # fW=fft2(W) - # # zfW=ZeroPad(fW,outshape=Npix) - # # WFull=ifft2(zfW) - - # # # #fW=ifft2(W) - # # # print W.shape - # fWTrue=fft2(WTrue) - # cfWTrue=fft2(np.conj(WTrue)) - # nc,_=fWTrue.shape - # xc=(nc-1)/2 - # dx=(Sups[i]-1)/2 - # x0,x1=xc-dx,xc+dx+1 - # #ZfWTrue=np.zeros_like(fWTrue) - # #ZfWTrue[x0:x1,x0:x1]=fWTrue[x0:x1,x0:x1] - # fzW=fWTrue[x0:x1,x0:x1] - # fzWconj=cfWTrue[x0:x1,x0:x1] - - # if_fzW=ifft2(fzW) - # cif_fzW=ifft2(fzWconj) - # z_if_fzW=ZeroPad(if_fzW,outshape=if_fzW.shape[0]*self.OverS) - # z_cif_fzW=ZeroPad(cif_fzW,outshape=if_fzW.shape[0]*self.OverS) - # f_z_if_fzW=fft2(z_if_fzW) - # f_z_cif_fzW=fft2(z_cif_fzW) - # # ifZfWTrue=ifft2(ZfWTrue) - - # # iffWTrue=ifft2(fWTrue[x0:x1,x0:x1]) - # # iffW=ifft2(fW) - - # # pylab.clf() - # # pylab.subplot(1,2,1) - # # pylab.imshow(np.real(fzW),interpolation="nearest")#,extent=(l.min(),l.max(),m.min(),m.max()),vmin=-np.pi,vmax=np.pi) - # # pylab.subplot(1,2,2) - # # pylab.imshow(np.real(f_z_if_fzW),interpolation="nearest")#,extent=(l.min(),l.max(),m.min(),m.max()),vmin=-np.pi,vmax=np.pi) - # # pylab.draw() - # # pylab.show(False) - # # pylab.pause(0.1) - # # fzW=f_z_if_fzW - # # fzWconj=f_z_cif_fzW - # # # #################### - - # T.timeit("3c") - - # # W=ThisSphe - + if Nw <= 1: + if not(Sups[0] % 2): + Sups[0] += 1 + dummy, dymmy, ThisSphe = self.SpheM.MakeSphe(Sups[0]) + W = np.abs(ThisSphe) zW = ZeroPad(W, outshape=W.shape[0]*self.OverS) - - # T.timeit("3d") - - # #### - # # W=np.abs(W) - # #### - zWconj = np.conj(zW) - - # #FFTWMachine=ModFFTW.FFTW_2Donly(W.shape,W.dtype, ncores = 1) - # #W=FFTWMachine.fft(W) - # #Wconj=FFTWMachine.fft(Wconj) fzW = fft2(zW) fzWconj = fft2(zWconj) - # T.timeit("3e") fzW = np.complex64(fzW).copy() fzWconj = np.complex64(fzWconj).copy() - # fzW.fill(2+3*1j) - # fzWconj.fill(2+3*1j) - fzW = self.GiveReorgCF(fzW) fzWconj = self.GiveReorgCF(fzWconj) @@ -560,7 +465,123 @@ def InitW(self): fzWconj = np.require(fzWconj.copy(), requirements=["A", "C"]) Wplanes.append(fzW) WplanesConj.append(fzWconj) - # T.timeit("3f") + else: + for i in range(Nw): + #print>>log, "%i/%i"%(i,Nw) + if not(Sups[i] % 2): + Sups[i] += 1 + dummy, dymmy, ThisSphe = self.SpheM.MakeSphe(Sups[i]) + wl = w[i]/waveMin + + # ############## + # l,m=np.mgrid[-lrad:lrad:Npix*1j,-lrad:lrad:Npix*1j] + # n_1=np.sqrt(1.-l**2-m**2)-1 + # WTrue=np.exp(-2.*1j*np.pi*wl*(n_1))*self.ifzfCF + # ############## + + DX = 2*lrad/Sups[i] + l, m = np.mgrid[-lrad+DX/2:lrad-DX/2:Sups[i] + * 1j, -lrad+DX/2:lrad-DX/2:Sups[i]*1j] + # l,m=np.mgrid[-lrad:lrad:Sups[i]*1j,-lrad:lrad:Sups[i]*1j] + # n_1=np.sqrt(1.-l**2-m**2)-1 + n_1 = ModFitPoly2D.polyval2d(l, m, CoefPoly) + # n_1=n_1.T[::-1,:] + T.timeit("3a") + + # stop + # n_1=np.sqrt(1.-(l-l0)**2-(m-m0)**2)-1 + # n_1=(1./np.sqrt(1.-l0**2-m0**2))*(l0*l+m0*m) + W = np.exp(-2.*1j*np.pi*wl*(n_1)) + #import pylab + + # pylab.clf() + # pylab.imshow(np.angle(W),interpolation="nearest",extent=(l.min(),l.max(),m.min(),m.max()),vmin=-np.pi,vmax=np.pi) + # pylab.draw() + # pylab.show(False) + # pylab.pause(0.1) + # T.timeit("3b") + + # #### + # W.fill(1.) + # #### + W *= np.abs(ThisSphe) + + # # #################### + # # fW=fft2(W) + # # zfW=ZeroPad(fW,outshape=Npix) + # # WFull=ifft2(zfW) + + # # # #fW=ifft2(W) + # # # print W.shape + # fWTrue=fft2(WTrue) + # cfWTrue=fft2(np.conj(WTrue)) + # nc,_=fWTrue.shape + # xc=(nc-1)/2 + # dx=(Sups[i]-1)/2 + # x0,x1=xc-dx,xc+dx+1 + # #ZfWTrue=np.zeros_like(fWTrue) + # #ZfWTrue[x0:x1,x0:x1]=fWTrue[x0:x1,x0:x1] + # fzW=fWTrue[x0:x1,x0:x1] + # fzWconj=cfWTrue[x0:x1,x0:x1] + + # if_fzW=ifft2(fzW) + # cif_fzW=ifft2(fzWconj) + # z_if_fzW=ZeroPad(if_fzW,outshape=if_fzW.shape[0]*self.OverS) + # z_cif_fzW=ZeroPad(cif_fzW,outshape=if_fzW.shape[0]*self.OverS) + # f_z_if_fzW=fft2(z_if_fzW) + # f_z_cif_fzW=fft2(z_cif_fzW) + # # ifZfWTrue=ifft2(ZfWTrue) + + # # iffWTrue=ifft2(fWTrue[x0:x1,x0:x1]) + # # iffW=ifft2(fW) + + # # pylab.clf() + # # pylab.subplot(1,2,1) + # # pylab.imshow(np.real(fzW),interpolation="nearest")#,extent=(l.min(),l.max(),m.min(),m.max()),vmin=-np.pi,vmax=np.pi) + # # pylab.subplot(1,2,2) + # # pylab.imshow(np.real(f_z_if_fzW),interpolation="nearest")#,extent=(l.min(),l.max(),m.min(),m.max()),vmin=-np.pi,vmax=np.pi) + # # pylab.draw() + # # pylab.show(False) + # # pylab.pause(0.1) + # # fzW=f_z_if_fzW + # # fzWconj=f_z_cif_fzW + # # # #################### + + # T.timeit("3c") + + # # W=ThisSphe + + zW = ZeroPad(W, outshape=W.shape[0]*self.OverS) + + # T.timeit("3d") + + # #### + # # W=np.abs(W) + # #### + + zWconj = np.conj(zW) + + # #FFTWMachine=ModFFTW.FFTW_2Donly(W.shape,W.dtype, ncores = 1) + # #W=FFTWMachine.fft(W) + # #Wconj=FFTWMachine.fft(Wconj) + fzW = fft2(zW) + fzWconj = fft2(zWconj) + + # T.timeit("3e") + fzW = np.complex64(fzW).copy() + fzWconj = np.complex64(fzWconj).copy() + + # fzW.fill(2+3*1j) + # fzWconj.fill(2+3*1j) + + fzW = self.GiveReorgCF(fzW) + fzWconj = self.GiveReorgCF(fzWconj) + + fzW = np.require(fzW.copy(), requirements=["A", "C"]) + fzWconj = np.require(fzWconj.copy(), requirements=["A", "C"]) + Wplanes.append(fzW) + WplanesConj.append(fzWconj) + # T.timeit("3f") self.Wplanes = Wplanes self.WplanesConj = WplanesConj diff --git a/DDFacet/Parset/DefaultParset.cfg b/DDFacet/Parset/DefaultParset.cfg index ba896cf4..e0c63d06 100755 --- a/DDFacet/Parset/DefaultParset.cfg +++ b/DDFacet/Parset/DefaultParset.cfg @@ -121,7 +121,7 @@ DecorrLocation = Edge # where decorrelation is estimated #options:Ce _Help = Imager convolution function settings OverS = 11 # Oversampling factor. #type:int #metavar:N Support = 7 # CF support size. #type:int #metavar:N -Nw = 100 # Number of w-planes. #type:int #metavar:PLANES +Nw = 100 # Number of w-planes. Setting this to 1 enables AIPS style faceting. #type:int #metavar:PLANES wmax = 0 # Maximum w coordinate. Visibilities with larger w will not be gridded. If 0, no maximum is imposed. #type:float #metavar:METERS