diff --git a/DDFacet/Array/PrintRecArray.py b/DDFacet/Array/PrintRecArray.py index 657e9c62..474ad713 100644 --- a/DDFacet/Array/PrintRecArray.py +++ b/DDFacet/Array/PrintRecArray.py @@ -26,12 +26,17 @@ def remove_field_name(a, name): return b from prettytable import PrettyTable -def Print(CatIn,RemoveFieldName='ChanFreq'): +def Print(CatIn,RemoveFieldName='ChanFreq',HideList=None): if RemoveFieldName in CatIn.dtype.names: Cat=remove_field_name(CatIn, RemoveFieldName) else: Cat=CatIn + Cat=Cat.copy() + if HideList is not None: + for field in HideList: + Cat=remove_field_name(Cat, field) + x = PrettyTable(Cat.dtype.names) for row in Cat: x.add_row(row) print x diff --git a/DDFacet/Array/shared_dict.py b/DDFacet/Array/shared_dict.py index 588da015..fdb542bf 100644 --- a/DDFacet/Array/shared_dict.py +++ b/DDFacet/Array/shared_dict.py @@ -59,7 +59,7 @@ def load(self): except: print "Error loading item %s" % self.path traceback.print_exc() - return SharedDict.ItemLoadError(path, sys.exc_info()) + return SharedDict.ItemLoadError(self.path, sys.exc_info()) class SharedArrayProxy (ItemProxy): def load_impl(self): diff --git a/DDFacet/Data/ClassBeamMean.py b/DDFacet/Data/ClassBeamMean.py index 31c837a6..875d124d 100644 --- a/DDFacet/Data/ClassBeamMean.py +++ b/DDFacet/Data/ClassBeamMean.py @@ -42,7 +42,9 @@ def __init__(self,VS): self.DoCentralNorm=self.GD["Beam"]["CenterNorm"] self.SmoothBeam=None self.CheckCache() - if self.CacheValid: return + if self.CacheValid: + MyLogger.setLoud(["ClassJones","ClassLOFARBeam"]) + return #self.GD["Beam"]["CenterNorm"]=0 @@ -173,7 +175,6 @@ def StackBeam(self,ThisMSData,iDir): WW=Ws**2 T.timeit("4") - print>>log, 'WW.shape is', WW.shape WW=WW.reshape((1,ind.size,MSnchan)) T.timeit("5") JJsq=WW*JJ**2 diff --git a/DDFacet/Data/ClassFITSBeam.py b/DDFacet/Data/ClassFITSBeam.py index 875dfbf3..0b9f1b24 100644 --- a/DDFacet/Data/ClassFITSBeam.py +++ b/DDFacet/Data/ClassFITSBeam.py @@ -70,8 +70,9 @@ def __init__ (self, ms, opts): if not self.nchan: self.nchan = len(self.freqs) else: - chanstep = max(1, len(self.freqs) / self.nchan) - self.freqs = self.freqs[chanstep/2::chanstep] + cw = self.ms.ChanWidth.ravel() + fq = np.linspace(self.freqs[0]-cw[0]/2, self.freqs[-1]+cw[-1]/2, self.nchan+1) + self.freqs = (fq[:-1] + fq[1:])/2 feed = opts["FITSFeed"] if feed: @@ -179,6 +180,7 @@ def getFreqDomains (self): df = (self.freqs[1]-self.freqs[0])/2 if len(self.freqs)>1 else self.freqs[0] domains[:,0] = self.freqs-df domains[:,1] = self.freqs+df +# import pdb; pdb.set_trace() return domains def evaluateBeam (self, t0, ra, dec): diff --git a/DDFacet/Data/ClassJones.py b/DDFacet/Data/ClassJones.py index b1084f58..f2d6aea1 100644 --- a/DDFacet/Data/ClassJones.py +++ b/DDFacet/Data/ClassJones.py @@ -29,6 +29,7 @@ import ClassLOFARBeam import ClassFITSBeam # import ClassSmoothJones is not used anywhere, should be able to remove it +import tables class ClassJones(): @@ -50,7 +51,7 @@ def InitDDESols(self, DATA, quiet=False): if SolsFile != "": self.ApplyCal = True self.JonesNormSolsFile_killMS, valid = self.MS.cache.checkCache( - "JonesNorm_killMS.npz", + "JonesNorm_killMS", dict(VisData=GD["Data"], DDESolutions=GD["DDESolutions"], DataSelection=self.GD["Selection"], @@ -62,7 +63,7 @@ def InitDDESols(self, DATA, quiet=False): DicoSols, TimeMapping, DicoClusterDirs = self.DiskToSols(self.JonesNormSolsFile_killMS) else: DicoSols, TimeMapping, DicoClusterDirs = self.MakeSols("killMS", DATA, quiet=quiet) - self.MS.cache.saveCache("JonesNorm_killMS.npz") + self.MS.cache.saveCache("JonesNorm_killMS") DATA["killMS"] = dict(Jones=DicoSols, TimeMapping=TimeMapping, Dirs=DicoClusterDirs) self.DicoClusterDirs_kMS=DicoClusterDirs @@ -121,20 +122,26 @@ def SolsToDisk(self, OutName, DicoSols, DicoClusterDirs, TimeMapping): # np.savez(self.JonesNorm_killMS,l=l,m=m,I=I,Cluster=Cluster,t0=t0,t1=t1,tm=tm,Jones=Jones,TimeMapping=TimeMapping) - np.savez(file(OutName, "w"), + os.system("touch %s"%OutName) + np.savez(file("%s.npz"%OutName, "w"), l=l, m=m, I=I, Cluster=Cluster, t0=t0, t1=t1, tm=tm, ra=ra,dec=dec, - Jones=Jones, TimeMapping=TimeMapping, VisToJonesChanMapping=VisToJonesChanMapping) + np.save(file("%s.npy"%OutName, "w"), + Jones) + def DiskToSols(self, InName): # SolsFile_killMS=np.load(self.JonesNorm_killMS) # print>>log, " Loading %s"%InName # print "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!",InName - SolsFile = np.load(InName) - print>>log, " %s loaded" % InName + SolsFile = np.load("%s.npz"%InName) + print>>log, " %s.npz loaded" % InName + Jones=np.load("%s.npy"%InName) + print>>log, " %s.npy loaded" % InName + DicoClusterDirs = {} DicoClusterDirs["l"] = SolsFile["l"] @@ -147,13 +154,14 @@ def DiskToSols(self, InName): DicoSols["t0"] = SolsFile["t0"] DicoSols["t1"] = SolsFile["t1"] DicoSols["tm"] = SolsFile["tm"] - DicoSols["Jones"] = SolsFile["Jones"] + + DicoSols["Jones"] = Jones + DicoSols["VisToJonesChanMapping"] = SolsFile["VisToJonesChanMapping"] TimeMapping = SolsFile["TimeMapping"] return DicoSols, TimeMapping, DicoClusterDirs def MakeSols(self, StrType, DATA, quiet=False): - print>>log, "Build solution Dico for %s" % StrType if StrType == "killMS": @@ -317,6 +325,7 @@ def GiveKillMSSols(self): for File, ThisGlobalMode, ThisJonesMode in zip( SolsFileList, GlobalNormList, JonesNormList): + DicoClusterDirs, DicoSols, VisToJonesChanMapping = self.GiveKillMSSols_SingleFile( File, GlobalMode=ThisGlobalMode, JonesMode=ThisJonesMode) print>>log, " VisToJonesChanMapping: %s" % str(VisToJonesChanMapping) @@ -330,19 +339,8 @@ def GiveKillMSSols(self): return DicoClusterDirs, DicoJones - def GiveKillMSSols_SingleFile( - self, - SolsFile, - JonesMode="AP", - GlobalMode=""): - + def ReadNPZ(self,SolsFile): print>>log, " Loading solution file %s" % (SolsFile) - if not(".npz" in SolsFile): - Method = SolsFile - ThisMSName = reformat.reformat( - os.path.abspath(self.MS.MSName), - LastSlash=False) - SolsFile = "%s/killMS.%s.sols.npz" % (ThisMSName, Method) self.ApplyCal = True DicoSolsFile = np.load(SolsFile) @@ -360,6 +358,15 @@ def GiveKillMSSols_SingleFile( DicoClusterDirs["I"] = ClusterCat.SumI DicoClusterDirs["Cluster"] = ClusterCat.Cluster + Sols = DicoSolsFile["Sols"] + Sols = Sols.view(np.recarray) + DicoSols = {} + DicoSols["t0"] = Sols.t0 + DicoSols["t1"] = Sols.t1 + DicoSols["tm"] = (Sols.t1+Sols.t0)/2. + nt, nf, na, nd, _, _ = Sols.G.shape + G = np.swapaxes(Sols.G, 1, 3).reshape((nt, nd, na, nf, 2, 2)) + if "FreqDomains" in DicoSolsFile.keys(): FreqDomains = DicoSolsFile["FreqDomains"] VisToJonesChanMapping = self.GiveVisToJonesChanMapping(FreqDomains) @@ -368,24 +375,113 @@ def GiveKillMSSols_SingleFile( self.BeamTimes_kMS = DicoSolsFile["BeamTimes"] - Sols = DicoSolsFile["Sols"] - Sols = Sols.view(np.recarray) + return VisToJonesChanMapping,DicoClusterDirs,DicoSols,G + + + def ReadH5(self,SolsFile): + print>>log, " Loading H5 solution file %s" % (SolsFile) + + self.ApplyCal = True + H=tables.open_file(SolsFile) + raNode,decNode=H.root.sol000.source[:]["dir"].T + times=H.root.sol000.tec000.time[:] + lFacet, mFacet = self.FacetMachine.CoordMachine.radec2lm(raNode, decNode) + # nt, na, nd, 1 + tec=H.root.sol000.tec000.val[:] + scphase=H.root.sol000.scalarphase000.val[:] + H.close() + del(H) + + DicoClusterDirs = {} + DicoClusterDirs["l"] = lFacet + DicoClusterDirs["m"] = mFacet + DicoClusterDirs["ra"] = raNode + DicoClusterDirs["dec"] = decNode + DicoClusterDirs["I"] = np.ones((lFacet.size,),np.float32) + DicoClusterDirs["Cluster"] = np.arange(lFacet.size) + + ClusterCat=np.zeros((lFacet.size,),dtype=[('Name','|S200'), + ('ra',np.float),('dec',np.float), + ('l',np.float),('m',np.float), + ('SumI',np.float),("Cluster",int)]) + ClusterCat=ClusterCat.view(np.recarray) + ClusterCat.l=lFacet + ClusterCat.m=mFacet + ClusterCat.ra=raNode + ClusterCat.dec=decNode + ClusterCat.I=DicoClusterDirs["I"] + ClusterCat.Cluster=DicoClusterDirs["Cluster"] + self.ClusterCat = ClusterCat + + + + dts=times[1::]-times[0:-1] + if not np.max(np.abs(dts-dts[0]))<0.1: + raise ValueError("The solutions dt should be the same") + dt=dts[0] + t0=times-dt/2. + t1=times+dt/2. DicoSols = {} - DicoSols["t0"] = Sols.t0 - DicoSols["t1"] = Sols.t1 - DicoSols["tm"] = (Sols.t1+Sols.t0)/2. - nt, nf, na, nd, _, _ = Sols.G.shape - G = np.swapaxes(Sols.G, 1, 3).reshape((nt, nd, na, nf, 2, 2)) + DicoSols["t0"] = t0 + DicoSols["t1"] = t1 + DicoSols["tm"] = (t1+t0)/2. + + + nt, na, nd, _=tec.shape + tecvals=tec.reshape((nt,na,nd,1)) + #freqs=self.FacetMachine.VS.GlobalFreqs.reshape((1,1,1,-1)) + freqs=self.MS.ChanFreq.ravel() + + scphase=scphase.reshape((nt,na,nd,1)) + freqs=freqs.reshape((1,1,1,-1)) + phase = (-8.4479745e9 * tecvals/freqs) + scphase + # nt,na,nd,nf,1 + phase=np.swapaxes(phase,1,2) + nf=freqs.size + G=np.zeros((nt, nd, na, nf, 2, 2),np.complex64) + z=np.exp(1j*phase) + G[:,:,:,:,0,0]=z + G[:,:,:,:,1,1]=z + + VisToJonesChanMapping = np.int32(np.arange(self.MS.NSPWChan,)) + + #self.BeamTimes_kMS = DicoSolsFile["BeamTimes"] + + return VisToJonesChanMapping,DicoClusterDirs,DicoSols,G + + + + + def GiveKillMSSols_SingleFile( + self, + SolsFile, + JonesMode="AP", + GlobalMode=""): + + + if not ".h5" in SolsFile: + if not(".npz" in SolsFile): + Method = SolsFile + ThisMSName = reformat.reformat( + os.path.abspath(self.MS.MSName), + LastSlash=False) + SolsFile = "%s/killMS.%s.sols.npz" % (ThisMSName, Method) + VisToJonesChanMapping,DicoClusterDirs,DicoSols,G=self.ReadNPZ(SolsFile) + else: + VisToJonesChanMapping,DicoClusterDirs,DicoSols,G=self.ReadH5(SolsFile) + + nt, nd, na, nf, _, _ = G.shape # G[:,:,:,:,0,0]=0. # G[:,:,:,:,1,1]=0. # G[:,0,:,:,0,0]=1. # G[:,0,:,:,1,1]=1. - # G.fill(0) # print>>log, "!!!!!!!!!!!!!!" # #G[:,:,:,:,1,1]=G[:,:,:,:,0,0] - # G[:,:,:,:,0,0]=G[:,:,:,:,1,1] + # G.fill(0) + # G[:,:,:,:,0,0]=1 + # G[:,:,:,:,1,1]=1 # print>>log, "!!!!!!!!!!!!!!" if GlobalMode == "MeanAbsAnt": @@ -409,16 +505,16 @@ def GiveKillMSSols_SingleFile( G[:, :, :, :, 1, 1] /= gmean_abs if GlobalMode == "BLBased": - print>>log, " Normalising by the mean of the amplitude (against time, freq, antenna)" - gmean_abs = np.mean(np.mean( - np.mean( - np.abs(G[:, :, :, :, 0, 0]), - axis=0), - axis=1), - axis=1) - gmean_abs = gmean_abs.reshape((1, nd, 1, 1)) - G[:, :, :, :, 0, 0] /= gmean_abs - G[:, :, :, :, 1, 1] /= gmean_abs + # print>>log, " Normalising by the mean of the amplitude (against time, freq, antenna)" + # gmean_abs = np.mean(np.mean( + # np.mean( + # np.abs(G[:, :, :, :, 0, 0]), + # axis=0), + # axis=1), + # axis=1) + # gmean_abs = gmean_abs.reshape((1, nd, 1, 1)) + # G[:, :, :, :, 0, 0] /= gmean_abs + # G[:, :, :, :, 1, 1] /= gmean_abs print>>log, " Extracting correction factor per-baseline" #(nt, nd, na, nf, 2, 2) @@ -438,6 +534,7 @@ def GiveKillMSSols_SingleFile( G[:,iDir,:,:,1,1]=G[:,iDir,:,:,1,1]/gu + if GlobalMode == "SumBLBased": print>>log, " Normalising by the mean of the amplitude (against time, freq, antenna)" gmean_abs = np.mean(np.mean( diff --git a/DDFacet/Data/ClassMS.py b/DDFacet/Data/ClassMS.py index 035d29f6..9e15eec4 100755 --- a/DDFacet/Data/ClassMS.py +++ b/DDFacet/Data/ClassMS.py @@ -35,6 +35,7 @@ from DDFacet.Other.CacheManager import CacheManager from DDFacet.Array import NpShared import sidereal +from DDFacet.Array import PrintRecArray import datetime import DDFacet.ToolsDir.ModRotate @@ -70,7 +71,9 @@ def obs_detail(filename,field=0): # Stuff relying on an OBSERVATION table: if to is not None: # Time - tm = Time(to[0]['TIME_RANGE']/86400.0,format='mjd') + tm = Time(to[0]['TIME_RANGE']/86400.0, + scale="utc", + format='mjd') results['DATE-OBS'] = tm[0].iso.split()[0] # Object @@ -593,10 +596,13 @@ def ReadData(self,DATA,row0,row1, # cached from previous run) # In auto cache mode, cache key is the start time of the process. The cache is thus reset when first # touched, so we read the MS on the first major cycle, and cache subsequently. - cache_key = dict(time=self._start_time) + # cache_key = dict(time=self._start_time) # @o-smirnov: why not that? # cache_key = dict(data=self.GD["Data"]) + cache_key = dict(data=self.GD["Data"], + selection=self.GD["Selection"], + Comp=self.GD["Comp"]) metadata_path, metadata_valid = self.cache.checkCache("A0A1UVWT.npz", cache_key, ignore_key=(use_cache=="force")) else: metadata_valid = False @@ -904,8 +910,8 @@ def ReadMSInfo(self,DoPrint=True): table_all = self.GiveMainTable() self.empty = not table_all.nrows() if self.empty: - return print>>log, ModColor.Str("MS %s (field %d, ddid %d): no rows, skipping"%(self.MSName, self.Field, self.DDID)) + return # raise RuntimeError,"no rows in MS %s, check your Field/DDID/TaQL settings"%(self.MSName) #print MSname+'/ANTENNA' @@ -1190,6 +1196,8 @@ def UpdateFlags(self, flags, uvw, data, A0, A1, times): antenna_rows = [(A0 == A) | (A1 == A) for A in xrange(self.na)] # print>>log," row index formed" + + antenna_flagfrac = [flags1[rows].sum() / float(flags1[rows].size or 1) for rows in antenna_rows] print>> log, " flagged fractions per antenna: %s" % " ".join(["%.2f" % frac for frac in antenna_flagfrac]) @@ -1223,8 +1231,10 @@ def __str__(self): ll.append(" - Selection: %s, channels: %s" % (ModColor.Str(str(self.TaQL), col="green"), self.ChanSlice)) ll.append(" - Phase centre (field %d): (ra, dec)=(%s, %s) "%(self.Field, rad2hmsdms(self.rarad,Type="ra").replace(" ",":")\ ,rad2hmsdms(self.decrad,Type="dec").replace(" ","."))) - ll.append(" - Frequency = %s MHz"%str(self.reffreq/1e6)) + ll.append(" - Frequency = %s MHz"%str(np.mean(self.ChanFreq)/1e6)) ll.append(" - Wavelength = %5.2f meters"%(np.mean(self.wavelength_chan))) + Freqs=3.e8/self.wavelength_chan.ravel()/1e6 + ll.append(" - Bandwidth = %5.2f MHz"%(np.max(Freqs)-np.min(Freqs))) ll.append(" - Time bin = %4.1f seconds"%(self.dt)) ll.append(" - Total Integration time = %6.2f hours"%self.DTh) ll.append(" - Number of antenna = %i"%self.na) @@ -1249,10 +1259,11 @@ def PutVisColumn(self, colname, vis, row0, row1, likecol="DATA", sort_index=None vis = vis[:,::-1,:] print>>log, "writing column %s rows %d:%d"%(colname,row0,row1) t = self.GiveMainTable(readonly=False, ack=False) + # if sorting rows, rearrange vis array back into MS order # if not sorting, then using slice(None) for row has no effect if sort_index is not None: - reverse_index = np.empty(self.nRowRead,dtype=int) + reverse_index = np.empty(nrow,dtype=int) reverse_index[sort_index] = np.arange(0,nrow,dtype=int) else: reverse_index = slice(None) @@ -1262,9 +1273,9 @@ def PutVisColumn(self, colname, vis, row0, row1, likecol="DATA", sort_index=None try: vis0 = t.getcol(colname, row0, nrow) except RuntimeError: - vis0 = t.getcol("DATA", row0, snrow) - vis0[reverse_index, self.ChanSlice, :] = vis - t.putcol(colname, vis0, row0, now) + vis0 = t.getcol("DATA", row0, nrow) + vis0[:, self.ChanSlice, :] = vis[reverse_index, :, :] + t.putcol(colname, vis0, row0, nrow) else: t.putcol(colname, vis[reverse_index,:,:], row0, nrow) t.close() @@ -1398,13 +1409,14 @@ def CopyCol(self,Colin,Colout): def AddCol(self,ColName,LikeCol="DATA",quiet=False): t=table(self.MSName,readonly=False,ack=False) - if (ColName in t.colnames() and not self.GD["Predict"]["Overwrite"]): + if (ColName in t.colnames()):# and not self.GD["Predict"]["Overwrite"]): if not quiet: print>>log, " Column %s already in %s"%(ColName,self.MSName) t.close() return - elif (ColName in t.colnames() and self.GD["Predict"]["Overwrite"]): - t.removecols(ColName) + # elif (ColName in t.colnames() and self.GD["Predict"]["Overwrite"]): + # t.removecols(ColName) + print>>log, " Putting column %s in %s"%(ColName,self.MSName) desc=t.getcoldesc(LikeCol) desc["name"]=ColName @@ -1451,17 +1463,24 @@ def PutNewCol(self,Name,LikeCol="CORRECTED_DATA"): t.close() - def Rotate(self,DATA,RotateType=["uvw","vis"]): + def Rotate(self,DATA,RotateType=["uvw","vis"],Sense="ToTarget",DataFieldName="data"): #DDFacet.ToolsDir.ModRotate.Rotate(self,radec) - StrRAOld = rad2hmsdms(self.OldRadec[0],Type="ra").replace(" ",":") - StrDECOld = rad2hmsdms(self.OldRadec[1],Type="dec").replace(" ",".") - StrRA = rad2hmsdms(self.NewRadec[0],Type="ra").replace(" ",":") - StrDEC = rad2hmsdms(self.NewRadec[1],Type="dec").replace(" ",".") - print>>log, "Rotate %s"%(",".join(RotateType)) + if Sense=="ToTarget": + ra0,dec0=self.OldRadec + ra1,dec1=self.NewRadec + elif Sense=="ToPhaseCenter": + ra0,dec0=self.NewRadec + ra1,dec1=self.OldRadec + + StrRAOld = rad2hmsdms(ra0,Type="ra").replace(" ",":") + StrDECOld = rad2hmsdms(dec0,Type="dec").replace(" ",".") + StrRA = rad2hmsdms(ra1,Type="ra").replace(" ",":") + StrDEC = rad2hmsdms(dec1,Type="dec").replace(" ",".") + print>>log, "Rotate %s [Mode = %s]"%(",".join(RotateType),Sense) print>>log, " from [%s, %s]"%(StrRAOld,StrDECOld) print>>log, " to [%s, %s]"%(StrRA,StrDEC) - DDFacet.ToolsDir.ModRotate.Rotate2(self.OldRadec,self.NewRadec,DATA["uvw"],DATA["data"],self.wavelength_chan, + DDFacet.ToolsDir.ModRotate.Rotate2((ra0,dec0),(ra1,dec1),DATA["uvw"],DATA[DataFieldName],self.wavelength_chan, RotateType=RotateType) diff --git a/DDFacet/Data/ClassVisServer.py b/DDFacet/Data/ClassVisServer.py index af42c40b..c092bdff 100644 --- a/DDFacet/Data/ClassVisServer.py +++ b/DDFacet/Data/ClassVisServer.py @@ -139,6 +139,7 @@ def Init(self, PointingID=0): DicoSelectOptions = self.DicoSelectOptions, get_obs_detail=get_detail) if MS.empty: + print>>log,"" continue self.ListMS.append(MS) # accumulate global set of frequencies, and min/max frequency @@ -387,6 +388,10 @@ def visPutColumnHandler (self, DATA, field, column, likecol): iMS, iChunk = DATA["iMS"], DATA["iChunk"] ms = self.ListMS[iMS] row0, row1 = ms.getChunkRow0Row1()[iChunk] + if ms.ToRADEC is not None: + ms.Rotate(DATA,RotateType=["vis"],Sense="ToPhaseCenter",DataFieldName=field) + + ms.PutVisColumn(column, DATA[field], row0, row1, likecol=likecol, sort_index=DATA["sort_index"]) def collectPutColumnResults(self): @@ -430,7 +435,7 @@ def startChunkLoadInBackground(self): # tell the IO thread to start loading the chunk APP.runJob(self._next_chunk_name, self._handler_LoadVisChunk, args=(self._next_chunk_name, self.iCurrentMS, self.iCurrentChunk), - io=0) + io=0)#,serial=True) return self._next_chunk_label def collectLoadedChunk(self, start_next=True): @@ -753,7 +758,7 @@ def _CalcWeights_handler(self): for ichunk in xrange(len(ms.getChunkRow0Row1())): APP.runJob("FinalizeWeights:%d:%d" % (ims, ichunk), self._finalizeWeights_handler, args=(self._weight_grid.readonly(), - self._weight_dict[ims][ichunk].readwrite()), + self._weight_dict.readwrite(),ims,ichunk,self._uvmax), counter=self._weightjob_counter, collect_result=False) APP.awaitJobCounter(self._weightjob_counter, progress="Finalize weights") # delete stuff @@ -865,8 +870,8 @@ def _accumulateWeights_handler (self, wg, msw, ims, ichunk, freqs, cell, npix, n y = uv[:, 1, :] x += xymax # offset, since X grid starts at -xymax # convert to index array -- this gives the number of the uv-bin on the grid - index = msw.addSharedArray("index", (uv.shape[0], len(freqs)), np.int64) - #index = np.zeros((uv.shape[0], len(freqs)), np.int64) + #index = msw.addSharedArray("index", (uv.shape[0], len(freqs)), np.int64) + index = np.zeros((uv.shape[0], len(freqs)), np.int64) index[...] = y * npixx + x # if we're in per-band weighting mode, then adjust the index to refer to each band's grid if nbands > 1: @@ -881,12 +886,54 @@ def _accumulateWeights_handler (self, wg, msw, ims, ichunk, freqs, cell, npix, n # print>>log,weights,index _pyGridderSmearPols.pyAccumulateWeightsOntoGrid(wg["grid"], weights.ravel(), index.ravel()) - def _finalizeWeights_handler(self, wg, msw): + def _finalizeWeights_handler(self, wg, mswAll,ims,ichunk,uvmax): + msw=mswAll[ims][ichunk] if "weight" in msw: + ms = self.ListMS[ims] + row0, row1 = ms.getChunkRow0Row1()[ichunk] + + msfreqs = ms.ChanFreq + freqs=msfreqs + nrows = row1 - row0 + chanslice = ms.ChanSlice + if not nrows: + return + tab = ms.GiveMainTable() + uvw = tab.getcol("UVW", row0, nrows) + tab.close() + nch, npol, npixIm, _ = self.FullImShape + FOV = self.CellSizeRad * npixIm + nbands = self.NFreqBands + cell = 1. / (self.Super * FOV) + xymax = int(math.floor(uvmax / cell)) + 1 + # grid will be from [-xymax,xymax] in U and [0,xymax] in V + npixx = xymax * 2 + 1 + npixy = xymax + 1 + npix = npixx * npixy + + uv = uvw[:,:2] + # flip sign of negative v values -- we'll only grid the top half of the plane + uv[uv[:, 1] < 0] *= -1 + # convert u/v to lambda, and then to pixel offset + uv = uv[..., np.newaxis] * freqs[np.newaxis, np.newaxis, :] / _cc + uv = np.floor(uv / cell).astype(int) + # u is offset, v isn't since it's the top half + x = uv[:, 0, :] + y = uv[:, 1, :] + x += xymax # offset, since X grid starts at -xymax + # convert to index array -- this gives the number of the uv-bin on the grid + index = np.zeros((uv.shape[0], len(freqs)), np.int64) + index[...] = y * npixx + x + + weight = msw["weight"] + #index[weight == 0] = 0 if self.Weighting != "natural": grid = wg["grid"].reshape((wg["grid"].size,)) - weight /= grid[msw["index"]] + #weight /= grid[msw["index"]] + index[index>=len(grid)]=0 + weight /= grid[index] + np.save(msw["cachepath"], weight) msw.delete_item("weight") if "index" in msw: diff --git a/DDFacet/Imager/ClassCasaImage.py b/DDFacet/Imager/ClassCasaImage.py index f9528f7f..801deb37 100644 --- a/DDFacet/Imager/ClassCasaImage.py +++ b/DDFacet/Imager/ClassCasaImage.py @@ -35,7 +35,16 @@ def FileToArray(FileName,CorrT): """ Read a FITS FileName file to an array """ hdu=fits.open(FileName) NormImage=hdu[0].data - nch,npol,_,_=NormImage.shape + + + if len(NormImage.shape)==4: + nch,npol,_,_=NormImage.shape + else: + nch,nx,ny=NormImage.shape + npol=1 + Sh=(nch,1,nx,ny) + NormImage=NormImage.reshape(Sh) + if CorrT: for ch in range(nch): for pol in range(npol): diff --git a/DDFacet/Imager/ClassDeconvMachine.py b/DDFacet/Imager/ClassDeconvMachine.py index 2bec6cd3..eeedee1e 100644 --- a/DDFacet/Imager/ClassDeconvMachine.py +++ b/DDFacet/Imager/ClassDeconvMachine.py @@ -36,6 +36,7 @@ from DDFacet.Other import MyLogger import traceback from DDFacet.ToolsDir.ModToolBox import EstimateNpix +from DDFacet.ToolsDir.ClassAdaptShape import ClassAdaptShape import copy from DDFacet.Other import AsyncProcessPool from DDFacet.Other.AsyncProcessPool import APP @@ -240,6 +241,18 @@ def Init(self): from DDFacet.Imager.HOGBOM import ClassImageDeconvMachineHogbom self.DeconvMachine=ClassImageDeconvMachineHogbom.ClassImageDeconvMachine(**MinorCycleConfig) print>>log,"Using Hogbom algorithm" + elif self.GD["Deconv"]["Mode"]=="MORESANE": + if MinorCycleConfig["ImagePolDescriptor"] != ["I"]: + raise NotImplementedError("Multi-polarization is not supported in MORESANE") + from DDFacet.Imager.MORESANE import ClassImageDeconvMachineMoresane + self.DeconvMachine=ClassImageDeconvMachineMoresane.ClassImageDeconvMachine(MainCache=self.VS.maincache, **MinorCycleConfig) + print>>log,"Using MORESANE algorithm" + elif self.GD["Deconv"]["Mode"]=="MUFFIN": + if MinorCycleConfig["ImagePolDescriptor"] != ["I"]: + raise NotImplementedError("Multi-polarization is not supported in MORESANE") + from DDFacet.Imager.MUFFIN import ClassImageDeconvMachineMUFFIN + self.DeconvMachine=ClassImageDeconvMachineMUFFIN.ClassImageDeconvMachine(MainCache=self.VS.maincache, **MinorCycleConfig) + print>>log,"Using MUFFIN algorithm" else: raise NotImplementedError("Unknown --Deconvolution-Mode setting '%s'" % self.GD["Deconv"]["Mode"]) self.ImageNoiseMachine.setMainCache(self.VS.maincache) @@ -366,7 +379,7 @@ def _checkForCachedDirty (self, sparsify, key=None): cachepath = self.VS.maincache.getElementPath("Dirty") valid = os.path.exists(cachepath) if not valid: - print>> log, ModColor.Str("Can't force-read cached dirty %s: does not exist", col="red") + print>> log, ModColor.Str("Can't force-read cached dirty %s: does not exist" % cachepath, col="red") raise RuntimeError("--Cache-Dirty forcedirty in effect, but no cached dirty image found") print>> log, ModColor.Str("Forcing reading the cached dirty image", col="red") writecache = False @@ -375,7 +388,7 @@ def _checkForCachedDirty (self, sparsify, key=None): valid = os.path.exists(cachepath) if not valid: - print>> log, ModColor.Str("Can't force-read cached last residual %s: does not exist", col="red") + print>> log, ModColor.Str("Can't force-read cached last residual %s: does not exist" % cachepath, col="red") raise RuntimeError("--Cache-Dirty forceresidual in effect, but no cached residual image found") print>> log, ModColor.Str("Forcing reading the cached last residual image", col="red") @@ -398,6 +411,23 @@ def _loadCachedPSF (self, cachepath): self.PSFGaussPars=self.DicoImagesPSF["PSFGaussPars"] self.PSFSidelobes=self.DicoImagesPSF["PSFSidelobes"] (self.FWHMBeamAvg, self.PSFGaussParsAvg, self.PSFSidelobesAvg)=self.DicoImagesPSF["EstimatesAvgPSF"] + + # #########################" + # Needed if cached PSF is there but --Output-RestoringBeam set differently + forced_beam=self.GD["Output"]["RestoringBeam"] + if forced_beam is not None: + if isinstance(forced_beam,float) or isinstance(forced_beam,int): + forced_beam=[float(forced_beam),float(forced_beam),0] + elif len(forced_beam)==1: + forced_beam=[forced_beam[0],forced_beam[0],0] + f_beam=(forced_beam[0]/3600.0,forced_beam[1]/3600.0,forced_beam[2]) + FWHMFact = 2. * np.sqrt(2. * np.log(2.)) + f_gau=(np.deg2rad(f_beam[0])/FWHMFact,np.deg2rad(f_beam[1])/FWHMFact,np.deg2rad(f_beam[2])) + print>>log, 'Will use user-specified beam: bmaj=%f, bmin=%f, bpa=%f degrees' % f_beam + beam, gausspars = f_beam, f_gau + self.FWHMBeamAvg, self.PSFGaussParsAvg = beam, gausspars + # #########################" + self.HasFittedPSFBeam=True @@ -536,9 +566,10 @@ def GiveDirty(self, psf=False, sparsify=0, last_cycle=False): if self.DicoDirty["JonesNorm"] is not None: self.FacetMachine.setNormImages(self.DicoDirty) - self.FacetMachinePSF.setNormImages(self.DicoDirty) - self.MeanJonesNorm = self.FacetMachinePSF.MeanJonesNorm - self.JonesNorm = self.FacetMachinePSF.JonesNorm + self.MeanJonesNorm = self.FacetMachine.MeanJonesNorm + self.JonesNorm = self.FacetMachine.JonesNorm + if self.FacetMachinePSF is not None: + self.FacetMachinePSF.setNormImages(self.DicoDirty) elif self.DicoImagesPSF["JonesNorm"] is not None: self.FacetMachine.setNormImages(self.DicoImagesPSF) self.FacetMachinePSF.setNormImages(self.DicoImagesPSF) @@ -584,9 +615,21 @@ def GiveDirty(self, psf=False, sparsify=0, last_cycle=False): # get loaded chunk from I/O thread, schedule next chunk # self.VS.startChunkLoadInBackground() DATA = self.VS.collectLoadedChunk(start_next=True) + if type(DATA) is str: print>>log,ModColor.Str("no more data: %s"%DATA, col="red") break + + # Allow for predict mode when a residual only is computed + predict_colname = None + if self.GD["Output"]["Mode"]=="Dirty": + predict_colname = self.GD["Predict"]["ColName"] + if self.DoDirtySub and predict_colname: + predict = DATA.addSharedArray("predict", DATA["datashape"], DATA["datatype"]) + visdata = DATA["data"] + np.copyto(predict, visdata) + + # None weights indicates an all-flagged chunk: go on to the next chunk if DATA["Weights"] is None: continue @@ -598,8 +641,8 @@ def GiveDirty(self, psf=False, sparsify=0, last_cycle=False): model_freqs = DATA["FreqMappingDegrid"] if not np.array_equal(model_freqs, current_model_freqs): ModelImage = self.FacetMachine.setModelImage(self.ModelMachine.GiveModelImage(model_freqs)) - # self.FacetMachine.ToCasaImage(ModelImage,ImageName="%s.model"%(self.BaseName), - # Fits=True,Stokes=self.VS.StokesConverter.RequiredStokesProducts()) + self.FacetMachine.ToCasaImage(ModelImage,ImageName="%s.model"%(self.BaseName), + Fits=True,Stokes=self.VS.StokesConverter.RequiredStokesProducts()) current_model_freqs = model_freqs print>> log, "model image @%s MHz (min,max) = (%f, %f)" % ( str(model_freqs / 1e6), ModelImage.min(), ModelImage.max()) @@ -609,6 +652,13 @@ def GiveDirty(self, psf=False, sparsify=0, last_cycle=False): self.FacetMachine.getChunkInBackground(DATA) self.FacetMachine.collectDegriddingResults() + if predict_colname: + predict -= visdata + # schedule jobs for saving visibilities, then start reading next chunk (both are on io queue) + self.VS.startVisPutColumnInBackground(DATA, "predict", predict_colname, likecol=self.GD["Data"]["ColName"]) + + + # crude but we need it here, since FacetMachine computes/loads CFs, which FacetMachinePSF uses. # so even if we're not using FM to make a dirty, we still need this call to make sure the CFs come in. self.FacetMachine.awaitInitCompletion() @@ -662,6 +712,8 @@ def GiveDirty(self, psf=False, sparsify=0, last_cycle=False): if psf and not psf_valid: self._finalizeComputedPSF(self.FacetMachinePSF, psf_writecache and psf_cachepath) + # self.SaveDirtyProducts() + # This call needs to be here to attach the cached smooth beam to FacetMachine if it exists # and if dirty has been initialised from cache self.FacetMachine.finaliseSmoothBeam() @@ -761,9 +813,9 @@ def GivePredict(self, subtract=False, from_fits=True): CleanMaskImage=None CleanMaskImageName=self.GD["Mask"]["External"] - if CleanMaskImageName is not None and CleanMaskImageName is not "": - print>>log,ModColor.Str("Will use mask image %s for the predict" % CleanMaskImageName) - CleanMaskImage = np.bool8(ClassCasaImage.FileToArray(CleanMaskImageName,True)) + # if CleanMaskImageName is not None and CleanMaskImageName is not "": + # print>>log,ModColor.Str("Will use mask image %s for the predict"%CleanMaskImageName) + # CleanMaskImage = np.bool8(ClassCasaImage.FileToArray(CleanMaskImageName,True)) modelfile = self.GD["Predict"]["FromImage"] @@ -771,6 +823,13 @@ def GivePredict(self, subtract=False, from_fits=True): if modelfile is not None and modelfile is not "": print>>log,ModColor.Str("Reading image file for the predict: %s" % modelfile) FixedModelImage = ClassCasaImage.FileToArray(modelfile,True) + nch,npol,_,NPix=self.FacetMachine.OutImShape + nchModel,npolModel,_,NPixModel=FixedModelImage.shape + if NPixModel!=NPix: + print>>log,ModColor.Str("Model image spatial shape does not match DDFacet settings [%i vs %i]"%(FixedModelImage.shape[-1],NPix)) + CA=ClassAdaptShape(FixedModelImage) + FixedModelImage=CA.giveOutIm(NPix) + if len(FixedModelImage.shape) != 4: raise RuntimeError("Expect FITS file with 4 axis: NX, NY, NPOL, NCH. Cannot continue.") nch, npol, ny, nx = FixedModelImage.shape @@ -824,7 +883,24 @@ def GivePredict(self, subtract=False, from_fits=True): print>> log, "reusing model image from previous chunk" else: if ModelImage is None: - ModelImage = self.FacetMachine.setModelImage(FixedModelImage) + nch=model_freqs.size + nchModel=FixedModelImage.shape[0] + ThisChFixedModelImage=FixedModelImage # so it's initialized + if nch!=nchModel: + print>>log,ModColor.Str("Model image spectral shape does not match DDFacet settings [%i vs %i]"%(nchModel,nch)) + if nchModel>nch: + print>>log,ModColor.Str(" taking the model's %i first channels only"%(nch)) + ThisChFixedModelImage=FixedModelImage[0:nch].copy() + else: + print>>log,ModColor.Str(" Replicating %i-times the 1st channel"%(nch)) + ThisChFixedModelImage=FixedModelImage[0].reshape((1,npol,NPix,NPix))*np.ones((DATA["ChanMappingDegrid"].size,1,1,1)) + self.FacetMachine.ToCasaImage(ThisChFixedModelImage, + ImageName="%s.cube.model"%(self.BaseName), + Fits=True, + Freqs=model_freqs, + Stokes=self.VS.StokesConverter.RequiredStokesProducts()) + ModelImage = self.FacetMachine.setModelImage(ThisChFixedModelImage) + if self.GD["Predict"]["MaskSquare"]: # MaskInside: choose mask inside (0) or outside (1) @@ -868,7 +944,7 @@ def GivePredict(self, subtract=False, from_fits=True): # Stokes=self.VS.StokesConverter.RequiredStokesProducts()) - if self.PredictMode == "BDA-degrid" or self.PredictMode == "DeGridder": # latter for backwards compatibility + if self.PredictMode == "BDA-degrid" or self.PredictMode == "Classic": # latter for backwards compatibility self.FacetMachine.getChunkInBackground(DATA) elif self.PredictMode == "Montblanc": from ClassMontblancMachine import ClassMontblancMachine @@ -993,7 +1069,21 @@ def main(self,NMajor=None): except: pass + + ### self.ModelMachine.ToFile(self.DicoModelName) + # ### + model_freqs=np.array([self.RefFreq],np.float64) + ModelImage = self.FacetMachine.setModelImage(self.DeconvMachine.GiveModelImage(model_freqs)) + # write out model image, if asked to + current_model_freqs = model_freqs + print>>log,"model image @%s MHz (min,max) = (%f, %f)"%(str(model_freqs/1e6),ModelImage.min(),ModelImage.max()) + if "o" in self._saveims: + self.FacetMachine.ToCasaImage(ModelImage, ImageName="%s.model%2.2i" % (self.BaseName, iMajor), + Fits=True, Freqs=current_model_freqs, + Stokes=self.VS.StokesConverter.RequiredStokesProducts()) + # stop + # ### ## returned with nothing done in minor cycle? Break out if not update_model or iMajor == NMajor: @@ -1347,7 +1437,8 @@ def fitSinglePSF(self, PSF, off, label="mean"): off = min(off, x[0], nx-x[0], y[0], ny-y[0]) print>> log, "Fitting %s PSF in a [%i,%i] box ..." % (label, off * 2, off * 2) P = PSF[0, x[0] - off:x[0] + off, y[0] - off:y[0] + off].copy() - (sidelobes), (bmaj, bmin, theta) = ModFitPSF.FindSidelobe(P) + bmaj, bmin, theta = ModFitPSF.FitCleanBeam(P) + sidelobes = ModFitPSF.FindSidelobe(P) print>>log, "PSF max is %f"%P.max() FWHMFact = 2. * np.sqrt(2. * np.log(2.)) @@ -1385,8 +1476,8 @@ def FitPSF(self): if forced_beam is not None: FWHMFact = 2. * np.sqrt(2. * np.log(2.)) - if isinstance(forced_beam,float): - forced_beam=[forced_beam,forced_beam,0] + if isinstance(forced_beam,float) or isinstance(forced_beam,int): + forced_beam=[float(forced_beam),float(forced_beam),0] elif len(forced_beam)==1: forced_beam=[forced_beam[0],forced_beam[0],0] f_beam=(forced_beam[0]/3600.0,forced_beam[1]/3600.0,forced_beam[2]) @@ -1539,7 +1630,7 @@ def RestoreAndShift(self): valid = os.path.exists(dirty_cachepath) if not valid: - print>> log, ModColor.Str("Can't force-read cached last residual %s: does not exist", col="red") + print>> log, ModColor.Str("Can't force-read cached last residual %s: does not exist" % dirty_cachepath, col="red") raise RuntimeError("--Cache-Dirty forceresidual in effect, but no cached residual image found") print>> log, ModColor.Str("Forcing reading the cached last residual image", col="red") @@ -1550,8 +1641,8 @@ def RestoreAndShift(self): cachepath = self.VS.maincache.getElementPath("PSF") valid = os.path.exists(cachepath) if not valid: - print>> log, ModColor.Str("Can't force-read cached last residual %s: does not exist", col="red") - raise RuntimeError("--Cache-Dirty forceresidual in effect, but no cached residual image found") + print>> log, ModColor.Str("Can't force-read cached PSF %s: does not exist" % cachepath, col="red") + raise RuntimeError("--Cache-PSF force in effect, but no cached PSF image found") print>> log, ModColor.Str("Forcing to read the cached PSF", col="red") self.DicoImagesPSF = shared_dict.create("FMPSF_AllImages") self.DicoImagesPSF.restore(cachepath) @@ -1850,7 +1941,13 @@ def alphamap(): label = 'alphamap' if label not in _images: _images.addSharedArray(label, intmodel().shape, np.float32) - _images[label] = ModelMachine.FreqMachine.alpha_map.reshape(intmodel().shape) + + # ############################## + # # Reverting for issue458 + #_images[label] = ModelMachine.FreqMachine.alpha_map.reshape(intmodel().shape) + _images[label] = ModelMachine.GiveSpectralIndexMap() + # ############################## + return _images[label] def alphaconvmap(): label = 'alphaconvmap' @@ -1989,10 +2086,13 @@ def alphaconvmap(): # Alpha image if "A" in self._saveims and self.VS.MultiFreqMode: - _images['alphaconvmap'] = alphaconvmap() - APP.runJob("save:alphaconv", self._saveImage_worker, io=0, args=(_images.readwrite(), 'alphaconvmap',), kwargs=dict( - ImageName="%s.alphaconv" % self.BaseName, Fits=True, delete=True, beam=self.FWHMBeamAvg, - Stokes=self.VS.StokesConverter.RequiredStokesProducts())) + # ############################## + # # Reverting for issue458 + # _images['alphaconvmap'] = alphaconvmap() + # APP.runJob("save:alphaconv", self._saveImage_worker, io=0, args=(_images.readwrite(), 'alphaconvmap',), kwargs=dict( + # ImageName="%s.alphaconv" % self.BaseName, Fits=True, delete=True, beam=self.FWHMBeamAvg, + # Stokes=self.VS.StokesConverter.RequiredStokesProducts())) + # ############################## _images['alphamap'] = alphamap() APP.runJob("save:alpha", self._saveImage_worker, io=0, args=(_images.readwrite(), 'alphamap',), kwargs=dict( ImageName="%s.alpha" % self.BaseName, Fits=True, delete=True, beam=self.FWHMBeamAvg, diff --git a/DDFacet/Imager/ClassFacetMachine.py b/DDFacet/Imager/ClassFacetMachine.py index 197d440f..e5a940b0 100644 --- a/DDFacet/Imager/ClassFacetMachine.py +++ b/DDFacet/Imager/ClassFacetMachine.py @@ -620,17 +620,24 @@ def initCFInBackground (self, other_fm=None): # subprocesses will place W-terms etc. here. Reset this first. self._CF = shared_dict.create("CFPSF" if self.DoPSF else "CF") # check if w-kernels, spacial weights, etc. are cached - cachekey = dict(ImagerCF=self.GD["CF"], - ImagerMainFacet=self.GD["Image"], - Facets=self.GD["Facets"], - RIME=self.GD["RIME"]) - cachename = self._cf_cachename = "CF" - # in oversize-PSF mode, make separate cache for PSFs - if self.DoPSF and self.Oversize != 1: - cachename = self._cf_cachename = "CFPSF" - cachekey["Oversize"] = self.Oversize - # check cache - cachepath, cachevalid = self.VS.maincache.checkCache(cachename, cachekey, directory=True) + + if self.GD["Cache"]["CacheCF"]: + cachekey = dict(ImagerCF=self.GD["CF"], + ImagerMainFacet=self.GD["Image"], + Facets=self.GD["Facets"], + RIME=self.GD["RIME"], + DDESolutions={"DDSols":self.GD["DDESolutions"]["DDSols"]}) + cachename = self._cf_cachename = "CF" + # in oversize-PSF mode, make separate cache for PSFs + if self.DoPSF and self.Oversize != 1: + cachename = self._cf_cachename = "CFPSF" + cachekey["Oversize"] = self.Oversize + # check cache + cachepath, cachevalid = self.VS.maincache.checkCache(cachename, cachekey, directory=True) + else: + print>>log,ModColor.Str("Explicitly not caching nor using cache for the Convolution Function") + cachepath, cachevalid="",False + # up to workers to load/save cache for iFacet in self.DicoImager.iterkeys(): facet_dict = self._CF.addSubdict(iFacet) @@ -720,14 +727,14 @@ def awaitInitCompletion (self): # mark cache as safe for res in workers_res: Type,path,iFacet=res - if Type=="compute": + if Type=="compute" and self.GD["Cache"]["CacheCF"]: #print iFacet facet_dict=self._CF[iFacet] d={} for key in facet_dict.keys(): d[key]=facet_dict[key] np.savez(file(path, "w"), **d) - self.VS.maincache.saveCache(self._cf_cachename) + if self.GD["Cache"]["CacheCF"]: self.VS.maincache.saveCache(self._cf_cachename) self.IsDDEGridMachineInit = True def setCasaImage(self, ImageName=None, Shape=None, Freqs=None, Stokes=["I"]): @@ -1206,7 +1213,7 @@ def FacetsToIm_Channel(self, kind="Dirty",ChanSel=None): NFacets=len(self.DicoImager.keys()) pBAR.render(0, NFacets) - # numexpr.set_num_threads(self.GD["Parallel"]["NCPU"]) # done in DDF.py + numexpr.set_num_threads(self.GD["Parallel"]["NCPU"]) # done in DDF.py for iFacet in self.DicoImager.keys(): diff --git a/DDFacet/Imager/ClassFacetMachineTessel.py b/DDFacet/Imager/ClassFacetMachineTessel.py index 8c37f048..c814e3ab 100755 --- a/DDFacet/Imager/ClassFacetMachineTessel.py +++ b/DDFacet/Imager/ClassFacetMachineTessel.py @@ -37,6 +37,7 @@ import os from DDFacet.ToolsDir.ModToolBox import EstimateNpix from DDFacet.Other import ModColor +import tables from DDFacet.Imager.ClassImToGrid import ClassImToGrid from matplotlib.path import Path @@ -97,7 +98,7 @@ def setFacetsLocs(self): if isinstance(SolsFile, list): SolsFile = self.GD["DDESolutions"]["DDSols"][0] - if SolsFile and (not (".npz" in SolsFile)): + if SolsFile and (not (".npz" in SolsFile)) and (not (".h5" in SolsFile)): Method = SolsFile ThisMSName = reformat.reformat( os.path.abspath(MSName), LastSlash=False) @@ -112,13 +113,20 @@ def setFacetsLocs(self): raNode = ClusterNodes.ra decNode = ClusterNodes.dec lFacet, mFacet = self.CoordMachine.radec2lm(raNode, decNode) - elif SolsFile != "": + elif ".npz" in SolsFile: print>> log, "Taking facet directions from solutions file: %s" % SolsFile ClusterNodes = np.load(SolsFile)["ClusterCat"] ClusterNodes = ClusterNodes.view(np.recarray) raNode = ClusterNodes.ra decNode = ClusterNodes.dec lFacet, mFacet = self.CoordMachine.radec2lm(raNode, decNode) + elif ".h5" in SolsFile: + print>> log, "Taking facet directions from HDF5 solutions file: %s" % SolsFile + H=tables.open_file(SolsFile) + raNode,decNode=H.root.sol000.source[:]["dir"].T + lFacet, mFacet = self.CoordMachine.radec2lm(raNode, decNode) + H.close() + del(H) else: print>> log, "Taking facet directions from regular grid" regular_grid = True diff --git a/DDFacet/Imager/ClassFrequencyMachine.py b/DDFacet/Imager/ClassFrequencyMachine.py index b2f2bef1..6b7e24a1 100644 --- a/DDFacet/Imager/ClassFrequencyMachine.py +++ b/DDFacet/Imager/ClassFrequencyMachine.py @@ -20,7 +20,7 @@ import numpy as np from DDFacet.ToolsDir import ClassRRGP -import matplotlib.pyplot as plt +# import matplotlib.pyplot as plt from DDFacet.Other import MyPickle class ClassFrequencyMachine(object): diff --git a/DDFacet/Imager/MORESANE/ClassImageDeconvMachineMoresane.py b/DDFacet/Imager/MORESANE/ClassImageDeconvMachineMoresane.py new file mode 100644 index 00000000..c55af3ae --- /dev/null +++ b/DDFacet/Imager/MORESANE/ClassImageDeconvMachineMoresane.py @@ -0,0 +1,305 @@ +''' +DDFacet, a facet-based radio imaging package +Copyright (C) 2013-2016 Cyril Tasse, l'Observatoire de Paris, +SKA South Africa, Rhodes University + +This program is free software; you can redistribute it and/or +modify it under the terms of the GNU General Public License +as published by the Free Software Foundation; either version 2 +of the License, or (at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program; if not, write to the Free Software +Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +''' + +import os +import numpy as np +from DDFacet.Other import MyLogger +from DDFacet.Other import ModColor +log=MyLogger.getLogger("ClassImageDeconvMachine") +from DDFacet.Array import NpParallel +from DDFacet.Array import NpShared +from DDFacet.ToolsDir import ModFFTW +from DDFacet.ToolsDir import ModToolBox +from DDFacet.Other import ClassTimeIt +from pyrap.images import image +from DDFacet.Imager.ClassPSFServer import ClassPSFServer +from DDFacet.Other.progressbar import ProgressBar +from DDFacet.Imager import ClassGainMachine +from DDFacet.Other import MyPickle +import multiprocessing +import time +from DDFacet.Imager.MORESANE.ClassMoresaneSingleSlice import ClassMoresaneSingleSlice +from DDFacet.Array import shared_dict +from DDFacet.ToolsDir import ClassSpectralFunctions +from scipy.optimize import least_squares +from DDFacet.ToolsDir.GiveEdges import GiveEdges + + +class ClassImageDeconvMachine(): + def __init__(self,GD=None,ModelMachine=None,RefFreq=None,*args,**kw): + self.GD=GD + self.ModelMachine = ModelMachine + self.RefFreq=RefFreq + if self.ModelMachine.DicoModel["Type"]!="MORESANE": + raise ValueError("ModelMachine Type should be MORESANE") + self.MultiFreqMode=(self.GD["Freq"]["NBand"]>1) + + def SetPSF(self,DicoVariablePSF): + self.PSFServer=ClassPSFServer(self.GD) + DicoVariablePSF=shared_dict.attach(DicoVariablePSF.path)#["CubeVariablePSF"] + self.PSFServer.setDicoVariablePSF(DicoVariablePSF) + self.PSFServer.setRefFreq(self.ModelMachine.RefFreq) + self.DicoVariablePSF=DicoVariablePSF + self.setFreqs(self.PSFServer.DicoMappingDesc) + + def setMaskMachine(self,MaskMachine): + self.MaskMachine=MaskMachine + + def setFreqs(self,DicoMappingDesc): + self.DicoMappingDesc=DicoMappingDesc + if self.DicoMappingDesc is None: return + self.SpectralFunctionsMachine=ClassSpectralFunctions.ClassSpectralFunctions(self.DicoMappingDesc,RefFreq=self.DicoMappingDesc["RefFreq"])#,BeamEnable=False) + self.SpectralFunctionsMachine.CalcFluxBands() + + def GiveModelImage(self,*args): return self.ModelMachine.GiveModelImage(*args) + + def Update(self,DicoDirty,**kwargs): + """ + Method to update attributes from ClassDeconvMachine + """ + #Update image dict + self.SetDirty(DicoDirty) + + def ToFile(self, fname): + """ + Write model dict to file + """ + self.ModelMachine.ToFile(fname) + + def FromFile(self, fname): + """ + Read model dict from file SubtractModel + """ + self.ModelMachine.FromFile(fname) + + def FromDico(self, DicoName): + """ + Read in model dict + """ + self.ModelMachine.FromDico(DicoName) + + def setSideLobeLevel(self,SideLobeLevel,OffsetSideLobe): + self.SideLobeLevel=SideLobeLevel + self.OffsetSideLobe=OffsetSideLobe + + def Init(self,**kwargs): + self.SetPSF(kwargs["PSFVar"]) + if "PSFSideLobes" not in self.DicoVariablePSF.keys(): + self.DicoVariablePSF["PSFSideLobes"]=kwargs["PSFAve"] + self.setSideLobeLevel(kwargs["PSFAve"][0], kwargs["PSFAve"][1]) + self.ModelMachine.setRefFreq(kwargs["RefFreq"]) + # store grid and degrid freqs for ease of passing to MSMF + #print kwargs["GridFreqs"],kwargs["DegridFreqs"] + self.GridFreqs=kwargs["GridFreqs"] + self.DegridFreqs=kwargs["DegridFreqs"] + self.ModelMachine.setFreqMachine(kwargs["GridFreqs"], kwargs["DegridFreqs"]) + + def SetDirty(self,DicoDirty): + self.DicoDirty=DicoDirty + self._Dirty=self.DicoDirty["ImageCube"] + self._MeanDirty=self.DicoDirty["MeanImage"] + NPSF=self.PSFServer.NPSF + _,_,NDirty,_=self._Dirty.shape + off=(NPSF-NDirty)/2 + self.DirtyExtent=(off,off+NDirty,off,off+NDirty) + self.ModelMachine.setModelShape(self._Dirty.shape) + + def AdaptArrayShape(self,A,Nout): + nch,npol,Nin,_=A.shape + if Nin==Nout: + return A + elif Nin>Nout: + # dx=Nout/2 + # B=np.zeros((nch,npol,Nout,Nout),A.dtype) + # print>>log," Adapt shapes: %s -> %s"%(str(A.shape),str(B.shape)) + # B[:]=A[...,Nin/2-dx:Nin/2+dx+1,Nin/2-dx:Nin/2+dx+1] + + N0=A.shape[-1] + xc0=yc0=N0/2 + N1=Nout + xc1=yc1=N1/2 + Aedge,Bedge=GiveEdges((xc0,yc0),N0,(xc1,yc1),N1) + x0d,x1d,y0d,y1d=Aedge + x0p,x1p,y0p,y1p=Bedge + B=A[...,x0d:x1d,y0d:y1d] + + return B + else: + return A + + def giveSliceCut(self,A,Nout): + nch,npol,Nin,_=A.shape + if Nin==Nout: + slice(None) + elif Nin>Nout: + N0=A.shape[-1] + xc0=yc0=N0/2 + if Nout%2==0: + x0d,x1d=xc0-Nout/2,xc0+Nout/2 + else: + x0d,x1d=xc0-Nout/2,xc0+Nout/2+1 + return slice(x0d,x1d) + else: + return None + + + + def updateModelMachine(self,ModelMachine): + self.ModelMachine=ModelMachine + if self.ModelMachine.RefFreq!=self.RefFreq: + raise ValueError("freqs should be equal") + + def updateMask(self,Mask): + nx,ny=Mask.shape + self._MaskArray = np.zeros((1,1,nx,ny),np.bool8) + self._MaskArray[0,0,:,:]=Mask[:,:] + + + def Deconvolve(self): + + + if self._Dirty.shape[-1]!=self._Dirty.shape[-2]: + # print "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" + # print self._Dirty.shape + # print "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" + return "MaxIter", True, True + + + dirty=self._Dirty + nch,npol,_,_=dirty.shape + Model=np.zeros_like(dirty) + + _,_,xp,yp=np.where(self._MeanDirty==np.max(self._MeanDirty)) + self.PSFServer.setLocation(xp,yp) + self.iFacet=self.PSFServer.iFacet + + psf,_=self.PSFServer.GivePSF() + + Nout=np.min([dirty.shape[-1],psf.shape[-1]]) + + if Nout%2!=0: Nout-=1 + + s_dirty_cut=self.giveSliceCut(dirty,Nout) + s_psf_cut=self.giveSliceCut(psf,2*Nout) + + if s_psf_cut is None: + print>>log, ModColor.Str("Could not adapt psf shape to 2*dirty shape!") + print>>log, ModColor.Str(" shapes are (dirty, psf) = [%s, %s]"%(str(dirty.shape),str(psf.shape))) + s_psf_cut=self.giveSliceCut(psf,Nout) + + + for ch in range(nch): + #print dirty[ch,0,s_dirty_cut,s_dirty_cut].shape + #print psf[ch,0,s_psf_cut,s_psf_cut].shape + CM=ClassMoresaneSingleSlice(dirty[ch,0,s_dirty_cut,s_dirty_cut], + psf[ch,0,s_psf_cut,s_psf_cut], + mask=None,GD=None) + model,resid=CM.giveModelResid(major_loop_miter=self.GD["MORESANE"]["NMajorIter"], + minor_loop_miter=self.GD["MORESANE"]["NMinorIter"], + loop_gain=self.GD["MORESANE"]["Gain"], + sigma_level=self.GD["MORESANE"]["SigmaCutLevel"],# tolerance=1., + enforce_positivity=self.GD["MORESANE"]["ForcePositive"]) + Model[ch,0,s_dirty_cut,s_dirty_cut]=model[:,:] + + # import pylab + # pylab.clf() + # pylab.subplot(2,2,1) + # pylab.imshow(dirty[ch,0,SliceDirty,SliceDirty],interpolation="nearest") + # pylab.colorbar() + + # pylab.subplot(2,2,2) + # pylab.imshow(psf[ch,0,SlicePSF,SlicePSF],interpolation="nearest") + # pylab.colorbar() + + # pylab.subplot(2,2,3) + # pylab.imshow(model,interpolation="nearest") + # pylab.colorbar() + + # pylab.subplot(2,2,4) + # pylab.imshow(resid,interpolation="nearest") + # pylab.colorbar() + + # pylab.draw() + # pylab.show() + + + # print + # print np.max(np.max(Model,axis=-1),axis=-1) + # print + # print + + + + + #_,_,nx,ny=Model.shape + #Model=np.mean(Model,axis=0).reshape((1,1,nx,ny)) + + #Model.fill(0) + #Model[:,:,xp,yp]=self._Dirty[:,:,xp,yp] + + if self.MultiFreqMode: + S,Al=self.DoSpectralFit(Model) + self.ModelMachine.setModel(S,0) + self.ModelMachine.setModel(Al,1) + else: + self.ModelMachine.setModel(Model,0) + + + return "MaxIter", True, True # stop deconvolution but do update model + + def DoSpectralFit(self,Model): + + def GiveResid(X,F,iFacet): + R=np.zeros_like(F) + S0,Alpha=X + for iBand in range(R.size): + #print iBand,self.SpectralFunctionsMachine.IntExpFunc(Alpha=np.array([0.]),iChannel=iBand,iFacet=iFacet) + R[iBand]=F[iBand]-S0*self.SpectralFunctionsMachine.IntExpFunc(Alpha=np.array([Alpha]).ravel(),iChannel=iBand,iFacet=iFacet) + + #stop + return R + + + nx,ny=Model.shape[-2],Model.shape[-1] + S=np.zeros((1,1,nx,ny),np.float32) + Al=np.zeros((1,1,nx,ny),np.float32) + + for iPix in range(Model.shape[-2]): + for jPix in range(Model.shape[-1]): + F=Model[:,0,iPix,jPix] + + JonesNorm=(self.DicoDirty["JonesNorm"][:,:,iPix,jPix]).reshape((-1,1,1,1)) + #W=self.DicoDirty["WeightChansImages"] + #JonesNorm=np.sum(JonesNorm*W.reshape((-1,1,1,1)),axis=0).reshape((1,1,1,1)) + + #F=F/np.sqrt(JonesNorm).ravel() + F0=np.mean(F) + if F0==0: + continue + + x0=(F0,-0.8) + + #print self.iFacet,iPix,jPix,F,F0 + X=least_squares(GiveResid, x0, args=(F,self.iFacet),ftol=1e-3,gtol=1e-3,xtol=1e-3) + x=X['x'] + S[0,0,iPix,jPix]=x[0] + Al[0,0,iPix,jPix]=x[1] + + return S,Al diff --git a/DDFacet/Imager/MORESANE/ClassImageDeconvMachineSSD.py b/DDFacet/Imager/MORESANE/ClassImageDeconvMachineSSD.py deleted file mode 100644 index 5d7257cd..00000000 --- a/DDFacet/Imager/MORESANE/ClassImageDeconvMachineSSD.py +++ /dev/null @@ -1,1053 +0,0 @@ -''' -DDFacet, a facet-based radio imaging package -Copyright (C) 2013-2016 Cyril Tasse, l'Observatoire de Paris, -SKA South Africa, Rhodes University - -This program is free software; you can redistribute it and/or -modify it under the terms of the GNU General Public License -as published by the Free Software Foundation; either version 2 -of the License, or (at your option) any later version. - -This program is distributed in the hope that it will be useful, -but WITHOUT ANY WARRANTY; without even the implied warranty of -MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -GNU General Public License for more details. - -You should have received a copy of the GNU General Public License -along with this program; if not, write to the Free Software -Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. -''' - -import numpy as np -from DDFacet.Other import MyLogger -from DDFacet.Other import ModColor -log=MyLogger.getLogger("ClassImageDeconvMachine") -from DDFacet.Array import NpParallel -from DDFacet.Array import NpShared -from DDFacet.ToolsDir import ModFFTW -from DDFacet.ToolsDir import ModToolBox -from DDFacet.Other import ClassTimeIt -from DDFacet.Imager import ClassMultiScaleMachine -from pyrap.images import image -from DDFacet.Imager.ClassPSFServer import ClassPSFServer -import sys -from DDFacet.Other.progressbar import ProgressBar -from DDFacet.Imager import ClassGainMachine -from SkyModel.PSourceExtract import ClassIslands -from SkyModel.PSourceExtract import ClassIncreaseIsland - -from DDFacet.Other import MyPickle - -import multiprocessing -import time - -MyLogger.setSilent("ClassArrayMethodGA") -MyLogger.setSilent("ClassIsland") - - - -class ClassImageDeconvMachine(): - def __init__(self,Gain=0.3, - MaxMinorIter=100,NCPU=6, - CycleFactor=2.5,FluxThreshold=None,RMSFactor=3,PeakFactor=0, - GD=None,SearchMaxAbs=1,CleanMaskImage=None,IdSharedMem="",ModelMachine=None, - **kw # absorb any unknown keywords arguments into this - ): - #self.im=CasaImage - - self.SearchMaxAbs=SearchMaxAbs - self.ModelImage=None - self.MaxMinorIter=MaxMinorIter - self.NCPU=NCPU - self.Chi2Thr=10000 - self.MaskArray=None - self.GD=GD - self.IdSharedMem=IdSharedMem - self.SubPSF=None - self.MultiFreqMode=(self.GD["Freq"]["NBand"]>1) - self.FluxThreshold = FluxThreshold - self.CycleFactor = CycleFactor - self.RMSFactor = RMSFactor - self.PeakFactor = PeakFactor - self.GainMachine=ClassGainMachine.ClassGainMachine(GainMin=Gain) - # reset overall iteration counter - self._niter = 0 - self.PSFCross=None - - self.IslandDeconvMode=self.GD["SSD"]["IslandDeconvMode"] # "GA" or "Moresane" or "Sasir" - if ModelMachine is None: - if self.IslandDeconvMode == "GA": - from DDFacet.Imager.GA import ClassModelMachineGA - self.ModelMachine = ClassModelMachineGA.ClassModelMachine(self.GD, GainMachine=self.GainMachine) - elif self.IslandDeconvMode == "Moresane": - from DDFacet.Imager.MORESANE import ClassModelMachineMORESANE - self.ModelMachine = ClassModelMachineMORESANE.ClassModelMachine(self.GD, GainMachine=self.GainMachine) - elif self.IslandDeconvMode == "Sasir": - raise NotImplementedError("ClassModelMachineSASIR is not implemented") - else: - # Trusting the user to pass correct ModelMachine for deconv algo - self.ModelMachine = ModelMachine - - if CleanMaskImage is not None: - print>>log, "Reading mask image: %s"%CleanMaskImage - MaskArray=image(CleanMaskImage).getdata() - nch,npol,_,_=MaskArray.shape - self._MaskArray=np.zeros(MaskArray.shape,np.bool8) - for ch in range(nch): - for pol in range(npol): - self._MaskArray[ch,pol,:,:]=np.bool8(1-MaskArray[ch,pol].T[::-1].copy())[:,:] - self.MaskArray=self._MaskArray[0] - self.IslandArray=np.zeros_like(self._MaskArray) - self.IslandHasBeenDone=np.zeros_like(self._MaskArray) - else: - print>>log, "You have to provide a mask image for SSD deconvolution" - - - - - def GiveModelImage(self,*args): return self.ModelMachine.GiveModelImage(*args) - - def setSideLobeLevel(self,SideLobeLevel,OffsetSideLobe): - self.SideLobeLevel=SideLobeLevel - self.OffsetSideLobe=OffsetSideLobe - - - def SetPSF(self,DicoVariablePSF): - self.PSFServer=ClassPSFServer(self.GD) - DicoVariablePSF["CubeVariablePSF"]=NpShared.ToShared("%s.CubeVariablePSF"%self.IdSharedMem,DicoVariablePSF["CubeVariablePSF"]) - self.PSFServer.setDicoVariablePSF(DicoVariablePSF) - #self.DicoPSF=DicoPSF - self.DicoVariablePSF=DicoVariablePSF - #self.NChannels=self.DicoDirty["NChannels"] - self.ModelMachine.setRefFreq(self.PSFServer.RefFreq) #,self.PSFServer.AllFreqs) - - def Init(self,**kwargs): - self.SetPSF(kwargs["PSFVar"]) - self.setSideLobeLevel(kwargs["PSFAve"][0], kwargs["PSFAve"][1]) - self.InitSSD() - - def InitSSD(self): - pass - - def AdaptArrayShape(self,A,Nout): - nch,npol,Nin,_=A.shape - if Nin==Nout: - return A - elif Nin>log,"Adapt mask shape" - self._MaskArray=self.AdaptArrayShape(self._MaskArray,NDirty) - self.MaskArray=self._MaskArray[0] - self.IslandArray=np.zeros_like(self._MaskArray) - self.IslandHasBeenDone=np.zeros_like(self._MaskArray) - - self.DirtyExtent=(off,off+NDirty,off,off+NDirty) - - if self.ModelImage is None: - self._ModelImage=np.zeros_like(self._Dirty) - self.ModelMachine.setModelShape(self._Dirty.shape) - if self.MaskArray is None: - self._MaskArray=np.zeros(self._Dirty.shape,dtype=np.bool8) - self.IslandArray=np.zeros_like(self._MaskArray) - self.IslandHasBeenDone=np.zeros_like(self._MaskArray) - - def CalcCrossIslandPSF(self,ListIslands): - print>>log," calculating global islands cross-contamination" - PSF=np.mean(self.PSFServer.DicoVariablePSF["MeanFacetPSF"][:,0],axis=0)#self.PSFServer.DicoVariablePSF["MeanFacetPSF"][0,0] - - nPSF,_=PSF.shape - xcPSF,ycPSF=nPSF/2,nPSF/2 - - IN=lambda x: ((x>=0)&(x=0)&(dx=0)&(dyTh)[0] - - - #Th=0.3 - #Flux=self.CrossFluxContrib[iIsland,iIsland] - #C0=(self.CrossFluxContrib[iIsland] > Flux*Th) - #indNearbyIsland=np.where(C0)[0] - - ii=0 - #print DicoIsland.keys() - #print>>log,"Looking around island #%i"%(iIsland) - for jIsland in indNearbyIsland: - #if jIsland in DicoIsland.keys(): - try: - Island=DicoIsland[jIsland] - #print>>log," merging island #%i -> #%i"%(jIsland,iIsland) - del(DicoIsland[jIsland]) - SubIslands=self.GiveNearbyIsland(DicoIsland,jIsland) - if SubIslands is not None: - Island+=SubIslands - return Island - except: - continue - - - #print>>log," could not find island #%i"%(iIsland) - - return None - -# JG Mod - def ToSquareIslands(self,ListIslands): - NIslands=len(ListIslands) - DicoSquareIslands={} - Dirty=self.DicoDirty["MeanImage"] - - for iIsland in range(NIslands): - x, y = np.array(ListIslands[iIsland]).T - minx,maxx=x.min(),x.max() - miny,maxy=y.min(),y.max() - nxisl,nyisl=maxx-minx+1,maxy-miny+1 # size of the smallest square around the island - npix=np.max([nxisl,nyisl]) - - if npix %2 == 0: - npix+=1 - - xc,yc=np.round((maxx+minx)/2),np.round((maxy+miny)/2) # compute island centre from minx,maxx,miny,maxy in the - - locxc,locyc=npix/2,npix/2 # center of the square around island - - SquareIsland=np.zeros((npix,npix),np.float32) - SquareIsland[0:npix,0:npix]=Dirty[0,0,xc-(npix-1)/2:xc+(npix-1)/2+1,yc-(npix-1)/2:yc+(npix-1)/2+1] - - SquareIsland_Mask=np.zeros((npix,npix),np.float32) - SquareIsland_Mask[locxc+(x-xc),locyc+(y-yc)]=1 # put 1 on the Island - #SquareIsland[locxc+(minx-xc):locxc+(maxx-xc)+1,locyc+(miny-yc):locyc+(maxy-yc)+1]=Dirty[0,0,minx:maxx+1,miny:maxy+1] - - DicoSquareIslands[iIsland]={"IslandCenter":(xc,yc), "IslandSquareData":SquareIsland, "IslandSquareMask":SquareIsland_Mask} - - return DicoSquareIslands - - def SquareIslandtoIsland(self, Model,ThisSquarePixList,ThisPixList): - ### Build ThisPixList from Model, in the reference frame of the Dirty - - xc,yc = ThisSquarePixList['IslandCenter'] # island center in original dirty - ListSquarePix_Data = ThisSquarePixList['IslandSquareData'] # square image of the dirty around Island center - ListSquarePix_Mask = ThisSquarePixList['IslandSquareMask'] # Corresponding square mask image - - NIslandPix=len(ThisPixList) - - Mod_x,Mod_y=Model.shape - SquarePix_x,SquarePix_y=ListSquarePix_Data.shape - - if Mod_x != SquarePix_x or Mod_y != SquarePix_y: - raise NameError('Mismatch between output Model image dims and original Square image dims. Please check if the even to uneven correction worked.') - - FluxV = [] - NewThisPixList = [] - for tmpcoor in ThisPixList: - currentx = tmpcoor[0] - currenty = tmpcoor[1] - x_loc_coor = (currentx - xc) + SquarePix_x / 2 # coordinates in the small Model image - y_loc_coor = (currenty - yc) + SquarePix_y / 2 # coordinates in the small Model image - if ListSquarePix_Mask[x_loc_coor, y_loc_coor] == 1: # if it is not masked (e.g. part of the island) - FluxV.append(ListSquarePix_Data[x_loc_coor, y_loc_coor]) - NewThisPixList.append([currentx, currenty]) - - return np.array(FluxV), np.array(NewThisPixList) - - def CalcCrossIslandFlux(self,ListIslands): - if self.PSFCross is None: - self.CalcCrossIslandPSF(ListIslands) - NIslands=len(ListIslands) - print>>log," grouping cross contaminating islands..." - - MaxIslandFlux=np.zeros((NIslands,),np.float32) - DicoIsland={} - - Dirty=self.DicoDirty["MeanImage"] - - - for iIsland in range(NIslands): - - x0,y0=np.array(ListIslands[iIsland]).T - PixVals0=Dirty[0,0,x0,y0] - MaxIslandFlux[iIsland]=np.max(PixVals0) - DicoIsland[iIsland]=ListIslands[iIsland] - - self.CrossFluxContrib=self.PSFCross*MaxIslandFlux.reshape((1,NIslands)) - - - NDone=0 - NJobs=NIslands - pBAR= ProgressBar(Title=" Group islands") - pBAR.disable() - pBAR.render(0, '%4i/%i' % (0,NJobs)) - - Th=0.05 - ListIslandMerged=[] - for iIsland in range(NIslands): - NDone+=1 - intPercent=int(100* NDone / float(NJobs)) - pBAR.render(intPercent, '%4i/%i' % (NDone,NJobs)) - - ThisIsland=self.GiveNearbyIsland(DicoIsland,iIsland) - - # indiIsland=np.where((self.PSFCross[iIsland])>Th)[0] - # ThisIsland=[] - # #print "Island #%i: %s"%(iIsland,str(np.abs(self.PSFCross[iIsland]))) - # for jIsland in indiIsland: - # if not(jIsland in DicoIsland.keys()): - # #print>>log," island #%i not there "%(jIsland) - # continue - # #print>>log," Putting island #%i in #%i"%(jIsland,iIsland) - # for iPix in range(len(DicoIsland[jIsland])): - # ThisIsland.append(DicoIsland[jIsland][iPix]) - # del(DicoIsland[jIsland]) - - - if ThisIsland is not None: - ListIslandMerged.append(ThisIsland) - - print>>log," have grouped %i --> %i islands"%(NIslands, len(ListIslandMerged)) - - return ListIslandMerged - - - - def SearchIslands(self,Threshold): - print>>log,"Searching Islands" - Dirty=self.DicoDirty["MeanImage"] - self.IslandArray[0,0]=(Dirty[0,0]>Threshold)|(self.IslandArray[0,0]) - #MaskImage=(self.IslandArray[0,0])&(np.logical_not(self._MaskArray[0,0])) - #MaskImage=(np.logical_not(self._MaskArray[0,0])) - MaskImage=(np.logical_not(self._MaskArray[0,0])) - Islands=ClassIslands.ClassIslands(Dirty[0,0],MaskImage=MaskImage, - MinPerIsland=0,DeltaXYMin=0) - Islands.FindAllIslands() - - ListIslands=Islands.LIslands - - print>>log," found %i islands"%len(ListIslands) - dx=self.GD["GAClean"]["NEnlargePars"] - if dx>0: - print>>log," increase their sizes by %i pixels"%dx - IncreaseIslandMachine=ClassIncreaseIsland.ClassIncreaseIsland() - for iIsland in range(len(ListIslands)):#self.NIslands): - ListIslands[iIsland]=IncreaseIslandMachine.IncreaseIsland(ListIslands[iIsland],dx=dx) - - - ListIslands=self.CalcCrossIslandFlux(ListIslands) - - - # FluxIslands=[] - # for iIsland in range(len(ListIslands)): - # x,y=np.array(ListIslands[iIsland]).T - # FluxIslands.append(np.sum(Dirty[0,0,x,y])) - # ind=np.argsort(np.array(FluxIslands))[::-1] - - # ListIslandsSort=[ListIslands[i] for i in ind] - - - # ListIslands=self.CalcCrossIslandFlux(ListIslandsSort) - self.ListIslands=[] - - for iIsland in range(len(ListIslands)): - x,y=np.array(ListIslands[iIsland]).T - PixVals=Dirty[0,0,x,y] - DoThisOne=False - - MaxIsland=np.max(np.abs(PixVals)) - - if (MaxIsland>(3.*self.RMS))|(MaxIsland>Threshold): - self.ListIslands.append(ListIslands[iIsland]) - - # ############################### - # if np.max(np.abs(PixVals))>Threshold: - # DoThisOne=True - # self.IslandHasBeenDone[0,0,x,y]=1 - # if ((DoThisOne)|self.IslandHasBeenDone[0,0,x[0],y[0]]): - # self.ListIslands.append(ListIslands[iIsland]) - # ############################### - - - self.NIslands=len(self.ListIslands) - print>>log," selected %i islands [out of %i] with peak flux > %.3g Jy"%(self.NIslands,len(ListIslands),Threshold) - - - self.ListIslands_Keep=self.ListIslands - - Sz=np.array([len(self.ListIslands[iIsland]) for iIsland in range(self.NIslands)]) - ind=np.argsort(Sz)[::-1] - - ListIslandsOut=[self.ListIslands[ind[i]] for i in ind] - self.ListIslands=ListIslandsOut - - # MORESANE MOD TO MAKE SQUARE IMAGE AROUND EACH ISLAND - if self.IslandDeconvMode == "Moresane" or self.IslandDeconvMode == "Sasir": # convert island to square image to pass to MORESANE & SASIR - self.ListSquareIslands = self.ToSquareIslands(self.ListIslands) - - - def setChannel(self,ch=0): - self.Dirty=self._MeanDirty[ch] - self.ModelImage=self._ModelImage[ch] - self.MaskArray=self._MaskArray[ch] - - - def GiveThreshold(self,Max): - return ((self.CycleFactor-1.)/4.*(1.-self.SideLobeLevel)+self.SideLobeLevel)*Max if self.CycleFactor else 0 - - def Deconvolve(self, *args, **kwargs): - - if self.GD['SSD']['Parallel'] == True: - print>>log,"Using Deconvolve Parallel for GA" - return self.DeconvolveParallel(*args,**kwargs) - else: - print>>log,"Using Deconvolve Serial for GA" - return self.DeconvolveSerial(*args, **kwargs) - - def DeconvolveSerial(self, ch=0): - """ - Runs minor cycle over image channel 'ch'. - initMinor is number of minor iteration (keeps continuous count through major iterations) - Nminor is max number of minor iteration - - Returns tuple of: return_code,continue,updated - where return_code is a status string; - continue is True if another cycle should be executed; - update is True if model has been updated (note that update=False implies continue=False) - """ - from DDFacet.Imager.GA.ClassEvolveGA import ClassEvolveGA - from DDFacet.Imager.MORESANE.ClassMoresane import ClassMoresane - from DDFacet.Imager.SASIR.ClassSasir import ClassSasir - if self._niter >= self.MaxMinorIter: - return "MaxIter", False, False - - self.setChannel(ch) - - _,npix,_=self.Dirty.shape - xc=(npix)/2 - - npol,_,_=self.Dirty.shape - - m0,m1=self.Dirty[0].min(),self.Dirty[0].max() - # pylab.clf() - # pylab.subplot(1,2,1) - # pylab.imshow(self.Dirty[0],interpolation="nearest",vmin=m0,vmax=m1) - # pylab.draw() - # pylab.show(False) - # pylab.pause(0.1) - - DoAbs=int(self.GD["Deconv"]["AllowNegative"]) - print>>log, " Running minor cycle [MinorIter = %i/%i, SearchMaxAbs = %i]"%(self._niter,self.MaxMinorIter,DoAbs) - - NPixStats=1000 - RandomInd=np.int64(np.random.rand(NPixStats)*npix**2) - RMS=np.std(np.real(self.Dirty.ravel()[RandomInd])) - self.RMS=RMS - - self.GainMachine.SetRMS(RMS) - - Fluxlimit_RMS = self.RMSFactor*RMS - - x,y,MaxDirty=NpParallel.A_whereMax(self.Dirty,NCPU=self.NCPU,DoAbs=DoAbs,Mask=self.MaskArray) - #MaxDirty=np.max(np.abs(self.Dirty)) - #Fluxlimit_SideLobe=MaxDirty*(1.-self.SideLobeLevel) - #Fluxlimit_Sidelobe=self.CycleFactor*MaxDirty*(self.SideLobeLevel) - Fluxlimit_Peak = MaxDirty*self.PeakFactor - Fluxlimit_Sidelobe = self.GiveThreshold(MaxDirty) - - mm0,mm1=self.Dirty.min(),self.Dirty.max() - - # work out uper threshold - StopFlux = max(Fluxlimit_Peak, Fluxlimit_RMS, Fluxlimit_Sidelobe, Fluxlimit_Peak, self.FluxThreshold) - - print>>log, " Dirty image peak flux = %10.6f Jy [(min, max) = (%.3g, %.3g) Jy]"%(MaxDirty,mm0,mm1) - print>>log, " RMS-based threshold = %10.6f Jy [rms = %.3g Jy; RMS factor %.1f]"%(Fluxlimit_RMS, RMS, self.RMSFactor) - print>>log, " Sidelobe-based threshold = %10.6f Jy [sidelobe = %.3f of peak; cycle factor %.1f]"%(Fluxlimit_Sidelobe,self.SideLobeLevel,self.CycleFactor) - print>>log, " Peak-based threshold = %10.6f Jy [%.3f of peak]"%(Fluxlimit_Peak,self.PeakFactor) - print>>log, " Absolute threshold = %10.6f Jy"%(self.FluxThreshold) - print>>log, " Stopping flux = %10.6f Jy [%.3f of peak ]"%(StopFlux,StopFlux/MaxDirty) - - - MaxModelInit=np.max(np.abs(self.ModelImage)) - - - # Fact=4 - # self.BookKeepShape=(npix/Fact,npix/Fact) - # BookKeep=np.zeros(self.BookKeepShape,np.float32) - # NPixBook,_=self.BookKeepShape - # FactorBook=float(NPixBook)/npix - - T=ClassTimeIt.ClassTimeIt() - T.disable() - - x,y,ThisFlux=NpParallel.A_whereMax(self.Dirty,NCPU=self.NCPU,DoAbs=DoAbs,Mask=self.MaskArray) - - if ThisFlux < StopFlux: - print>>log, ModColor.Str(" Initial maximum peak %g Jy below threshold, we're done here" % (ThisFlux),col="green" ) - return "FluxThreshold", False, False - - self.SearchIslands(StopFlux) - - for iIsland in range(self.NIslands): - ThisPixList=self.ListIslands[iIsland] - if self.IslandDeconvMode == "Moresane" or self.IslandDeconvMode == "Sasir": # convert island to square image to pass to MORESANE & SASIR - ThisSquarePixList=self.ListSquareIslands[iIsland] - - print>>log," Fitting island #%4.4i with %i pixels"%(iIsland,len(ThisPixList)) - - XY=np.array(ThisPixList,dtype=np.float32) - xm,ym=np.int64(np.mean(np.float32(XY),axis=0)) - - FacetID=self.PSFServer.giveFacetID2(xm,ym) - PSF=self.DicoVariablePSF["CubeVariablePSF"][FacetID] - - # self.DicoVariablePSF["CubeMeanVariablePSF"][FacetID] - - # FreqsInfo={"freqs":self.DicoVariablePSF["freqs"], - # "WeightChansImages":self.DicoVariablePSF["WeightChansImages"]} - - FreqsInfo=self.PSFServer.DicoMappingDesc - - - nchan,npol,_,_=self._Dirty.shape - JonesNorm=(self.DicoDirty["NormData"][:,:,xm,ym]).reshape((nchan,npol,1,1)) - W=self.DicoDirty["WeightChansImages"] - JonesNorm=np.sum(JonesNorm*W.reshape((nchan,1,1,1)),axis=0).reshape((1,npol,1,1)) - - - IslandBestIndiv=self.ModelMachine.GiveIndividual(ThisPixList) - - # COMMENT TO RUN IN DDF, UNCOMMENT TO TEST WITH TESTXXDeconv.py under XX folder - ################################ - #DicoSave={"Dirty":self._Dirty, - # "PSF":PSF, - # "FreqsInfo":FreqsInfo, - #"DicoMappingDesc":self.PSFServer.DicoMappingDesc, - # "ListPixData":ThisPixList, - # "ListPixParms":ThisPixList, - ## "ListSquarePix": ThisSquarePixList, ## When Moresane or SASIR are used - # "IslandBestIndiv":IslandBestIndiv, - # "GD":self.GD, - # "FacetID":FacetID} - #print "saving" - #MyPickle.Save(DicoSave, "SaveTest") - #print "saving ok" - #ipdb.set_trace() - - ################################ - - nch=nchan - self.FreqsInfo=FreqsInfo - WeightMeanJonesBand=self.FreqsInfo["MeanJonesBand"][FacetID].reshape((nch,1,1,1)) - WeightMueller=WeightMeanJonesBand.ravel() - WeightMuellerSignal=WeightMueller*self.FreqsInfo["WeightChansImages"].ravel() - - - if self.IslandDeconvMode=="GA": - CEv=ClassEvolveGA(self._Dirty,PSF,FreqsInfo,ListPixParms=ThisPixList, - ListPixData=ThisPixList,IslandBestIndiv=IslandBestIndiv, - GD=self.GD,WeightFreqBands=WeightMuellerSignal) - Model=CEv.main(NGen=100,DoPlot=False) - - if self.IslandDeconvMode=="Moresane" or self.IslandDeconvMode=="Sasir": - - # 0) Load Island info (center and square data) - ListSquarePix_Center = ThisSquarePixList['IslandCenter'] # island center in original dirty - ListSquarePix_Data = ThisSquarePixList['IslandSquareData'] # square image of the dirty around Island center - ListSquarePix_Mask= ThisSquarePixList['IslandSquareMask'] # Corresponding square mask image - xisland, yisland = ListSquarePix_Data.shape # size of the square postage stamp around island - - if self.IslandDeconvMode == "Moresane": - # 1) Shape PSF and Dirty to have even number of pixels (required by Moresane) - # DEAL WITH SQUARE DATA OF ISLAND IF UNEVEN - PSF_monochan=np.squeeze(PSF[0,0,:,:]) - # Single Channel for the moment - PSF_monochan_nx,PSF_monochan_ny=PSF_monochan.shape - - cropped_square_to_even = False - if xisland % 2 != 0: - # PSFCrop_even = np.zeros((xisland+1, xisland+1)) - # PSFCrop_even[:-1, :-1] = np.squeeze(PSFCrop) - Dirty_even = np.zeros((xisland - 1, xisland - 1)) - Dirty_even[:, :] = ListSquarePix_Data[:-1, :-1] - - PSF2_monochan=self.PSFServer.CropPSF(PSF_monochan,2*(xisland-1)+1) # PSF uneven cropped to double uneven dirty island - cropped_square_to_even = True - else: - Dirty_even = ListSquarePix_Data - PSF2_monochan=self.PSFServer.CropPSF(PSF_monochan,2*(xisland)+1) # PSF uneven cropped to double even dirty island - - PSF2_monochan_nx, PSF2_monochan_ny = PSF2_monochan.shape - - # DEAL WITH PSF IF UNEVEN (WILL ALWAYS BE UNEVEN EXCEPT IN PYMORESANE) - if PSF2_monochan_nx % 2 != 0: - PSF2_monochan_even = np.zeros((PSF2_monochan_nx-1, PSF2_monochan_nx-1)) - PSF2_monochan_even = PSF2_monochan[:-1,:-1] - else: - PSF2_monochan_even = PSF2_monochan - - - # 2) Run the actual MinorCycle algo - #Moresane = ClassMoresane(Dirty_even, PSF_monochan_even, DictMoresaneParms, GD=self.GD) - Moresane= ClassMoresane(Dirty_even, PSF2_monochan_even, self.GD['MORESANE'], GD=self.GD) - Model_Square,Residuals=Moresane.main() - - # 3) Apply Island mask to model to get rid of regions outside the island. - if cropped_square_to_even: # then restore the model to its original uneven dimension - Model_Square_uneven = np.zeros((xisland,xisland)) - Model_Square_uneven[:-1, :-1] = Model_Square - Model_Square = Model_Square_uneven - - Model_Square *= ListSquarePix_Mask # masking outside the island - - # 4) Convert back to Island format ( "S" and ThisPixList ) - NewModel, NewThisPixList=self.SquareIslandtoIsland(Model_Square, ThisSquarePixList, ThisPixList) - - Model = NewModel - ThisPixList = NewThisPixList - - if self.IslandDeconvMode == "Sasir": - # DO SASIR STUFF - # INCOMPLETE - - Sasir = ClassSasir(Dirty, PSF, self.GD['SASIR'], GD=self.GD) - Model_Square, Residuals = Sasir.main() - - # INCOMPLETE - - # Common command for every MinorCycle deconvolution algo - self.ModelMachine.AppendIsland(ThisPixList, Model) - - return "MaxIter", True, True # stop deconvolution but do update model - - - - - - def DeconvolveParallel(self, ch=0): - if self._niter >= self.MaxMinorIter: - return "MaxIter", False, False - - self.setChannel(ch) - - _,npix,_=self.Dirty.shape - xc=(npix)/2 - - npol,_,_=self.Dirty.shape - - m0,m1=self.Dirty[0].min(),self.Dirty[0].max() - - DoAbs=int(self.GD["Deconv"]["AllowNegative"]) - print>>log, " Running minor cycle [MinorIter = %i/%i, SearchMaxAbs = %i]"%(self._niter,self.MaxMinorIter,DoAbs) - - NPixStats=1000 - RandomInd=np.int64(np.random.rand(NPixStats)*npix**2) - RMS=np.std(np.real(self.Dirty.ravel()[RandomInd])) - self.RMS=RMS - - self.GainMachine.SetRMS(RMS) - - Fluxlimit_RMS = self.RMSFactor*RMS - - x,y,MaxDirty=NpParallel.A_whereMax(self.Dirty,NCPU=self.NCPU,DoAbs=DoAbs,Mask=self.MaskArray) - #MaxDirty=np.max(np.abs(self.Dirty)) - #Fluxlimit_SideLobe=MaxDirty*(1.-self.SideLobeLevel) - #Fluxlimit_Sidelobe=self.CycleFactor*MaxDirty*(self.SideLobeLevel) - Fluxlimit_Peak = MaxDirty*self.PeakFactor - Fluxlimit_Sidelobe = self.GiveThreshold(MaxDirty) - - mm0,mm1=self.Dirty.min(),self.Dirty.max() - - # work out uper threshold - StopFlux = max(Fluxlimit_Peak, Fluxlimit_RMS, Fluxlimit_Sidelobe, Fluxlimit_Peak, self.FluxThreshold) - - print>>log, " Dirty image peak flux = %10.6f Jy [(min, max) = (%.3g, %.3g) Jy]"%(MaxDirty,mm0,mm1) - print>>log, " RMS-based threshold = %10.6f Jy [rms = %.3g Jy; RMS factor %.1f]"%(Fluxlimit_RMS, RMS, self.RMSFactor) - print>>log, " Sidelobe-based threshold = %10.6f Jy [sidelobe = %.3f of peak; cycle factor %.1f]"%(Fluxlimit_Sidelobe,self.SideLobeLevel,self.CycleFactor) - print>>log, " Peak-based threshold = %10.6f Jy [%.3f of peak]"%(Fluxlimit_Peak,self.PeakFactor) - print>>log, " Absolute threshold = %10.6f Jy"%(self.FluxThreshold) - print>>log, " Stopping flux = %10.6f Jy [%.3f of peak ]"%(StopFlux,StopFlux/MaxDirty) - - - MaxModelInit=np.max(np.abs(self.ModelImage)) - - - # Fact=4 - # self.BookKeepShape=(npix/Fact,npix/Fact) - # BookKeep=np.zeros(self.BookKeepShape,np.float32) - # NPixBook,_=self.BookKeepShape - # FactorBook=float(NPixBook)/npix - - T=ClassTimeIt.ClassTimeIt() - T.disable() - - x,y,ThisFlux=NpParallel.A_whereMax(self.Dirty,NCPU=self.NCPU,DoAbs=DoAbs,Mask=self.MaskArray) - - if ThisFlux < StopFlux: - print>>log, ModColor.Str(" Initial maximum peak %g Jy below threshold, we're done here" % (ThisFlux),col="green" ) - return "FluxThreshold", False, False - - self.SearchIslands(StopFlux) - - - - # ================== Parallel part - NCPU=self.NCPU - work_queue = multiprocessing.Queue() - - - ListBestIndiv=[] - - NJobs=self.NIslands - T=ClassTimeIt.ClassTimeIt(" ") - T.disable() - for iIsland in range(self.NIslands): - # print "%i/%i"%(iIsland,self.NIslands) - ThisPixList=self.ListIslands[iIsland] - XY=np.array(ThisPixList,dtype=np.float32) - xm,ym=np.mean(np.float32(XY),axis=0) - T.timeit("xm,ym") - nchan,npol,_,_=self._Dirty.shape - JonesNorm=(self.DicoDirty["NormData"][:,:,xm,ym]).reshape((nchan,npol,1,1)) - W=self.DicoDirty["WeightChansImages"] - JonesNorm=np.sum(JonesNorm*W.reshape((nchan,1,1,1)),axis=0).reshape((1,npol,1,1)) - T.timeit("JonesNorm") - - IslandBestIndiv=self.ModelMachine.GiveIndividual(ThisPixList) - T.timeit("GiveIndividual") - ListBestIndiv.append(IslandBestIndiv) - FacetID=self.PSFServer.giveFacetID2(xm,ym) - T.timeit("FacetID") - - DicoOrder={"iIsland":iIsland, - "FacetID":FacetID, - "JonesNorm":JonesNorm} - - ListOrder=[iIsland,FacetID,JonesNorm.flat[0]] - - - work_queue.put(ListOrder) - T.timeit("Put") - - SharedListIsland="%s.ListIslands"%(self.IdSharedMem) - ListArrayIslands=[np.array(self.ListIslands[iIsland]) for iIsland in range(self.NIslands)] - NpShared.PackListArray(SharedListIsland,ListArrayIslands) - T.timeit("Pack0") - SharedBestIndiv="%s.ListBestIndiv"%(self.IdSharedMem) - NpShared.PackListArray(SharedBestIndiv,ListBestIndiv) - T.timeit("Pack1") - - - workerlist=[] - - # List_Result_queue=[] - # for ii in range(NCPU): - # List_Result_queue.append(multiprocessing.JoinableQueue()) - - - result_queue=multiprocessing.Queue() - - - for ii in range(NCPU): - W=WorkerDeconvIsland(work_queue, - result_queue, - # List_Result_queue[ii], - self.GD, - IdSharedMem=self.IdSharedMem, - FreqsInfo=self.PSFServer.DicoMappingDesc) - workerlist.append(W) - workerlist[ii].start() - - - print>>log, "Evolving %i generations of %i sourcekin"%(self.GD["GAClean"]["NMaxGen"],self.GD["GAClean"]["NSourceKin"]) - pBAR= ProgressBar(Title=" Evolve pop.") - #pBAR.disable() - pBAR.render(0, '%4i/%i' % (0,NJobs)) - - iResult=0 - while iResult < NJobs: - DicoResult=None - # for result_queue in List_Result_queue: - # if result_queue.qsize()!=0: - # try: - # DicoResult=result_queue.get_nowait() - - # break - # except: - - # pass - # #DicoResult=result_queue.get() - if result_queue.qsize()!=0: - try: - DicoResult=result_queue.get_nowait() - except: - pass - #DicoResult=result_queue.get() - - - if DicoResult is None: - time.sleep(0.05) - continue - - iResult+=1 - NDone=iResult - intPercent=int(100* NDone / float(NJobs)) - pBAR.render(intPercent, '%4i/%i' % (NDone,NJobs)) - - if DicoResult["Success"]: - iIsland=DicoResult["iIsland"] - ThisPixList=self.ListIslands[iIsland] - SharedIslandName="%s.FitIsland_%5.5i"%(self.IdSharedMem,iIsland) - Model=NpShared.GiveArray(SharedIslandName) - self.ModelMachine.AppendIsland(ThisPixList,Model) - NpShared.DelArray(SharedIslandName) - - - - for ii in range(NCPU): - workerlist[ii].shutdown() - workerlist[ii].terminate() - workerlist[ii].join() - - - return "MaxIter", True, True # stop deconvolution but do update model - - def Update(self, DicoDirty, **kwargs): - """ - Method to update attributes from ClassDeconvMachine - """ - # Update image dict - self.SetDirty(DicoDirty) - - def ToFile(self, fname): - """ - Write model dict to file - """ - self.ModelMachine.ToFile(fname) - - def FromFile(self, fname): - """ - Read model dict from file SubtractModel - """ - self.ModelMachine.FromFile(fname) - - def FromDico(self, DicoName): - """ - Read in model dict - """ - self.ModelMachine.FromDico(DicoName) - - ################################################################################### - ################################################################################### - - def GiveEdges(self,(xc0,yc0),N0,(xc1,yc1),N1): - M_xc=xc0 - M_yc=yc0 - NpixMain=N0 - F_xc=xc1 - F_yc=yc1 - NpixFacet=N1 - - ## X - M_x0=M_xc-NpixFacet/2 - x0main=np.max([0,M_x0]) - dx0=x0main-M_x0 - x0facet=dx0 - - M_x1=M_xc+NpixFacet/2 - x1main=np.min([NpixMain-1,M_x1]) - dx1=M_x1-x1main - x1facet=NpixFacet-dx1 - x1main+=1 - ## Y - M_y0=M_yc-NpixFacet/2 - y0main=np.max([0,M_y0]) - dy0=y0main-M_y0 - y0facet=dy0 - - M_y1=M_yc+NpixFacet/2 - y1main=np.min([NpixMain-1,M_y1]) - dy1=M_y1-y1main - y1facet=NpixFacet-dy1 - y1main+=1 - - Aedge=[x0main,x1main,y0main,y1main] - Bedge=[x0facet,x1facet,y0facet,y1facet] - return Aedge,Bedge - - - def SubStep(self,(dx,dy),LocalSM): - npol,_,_=self.Dirty.shape - x0,x1,y0,y1=self.DirtyExtent - xc,yc=dx,dy - N0=self.Dirty.shape[-1] - N1=LocalSM.shape[-1] - Aedge,Bedge=self.GiveEdges((xc,yc),N0,(N1/2,N1/2),N1) - factor=-1. - nch,npol,nx,ny=LocalSM.shape - x0d,x1d,y0d,y1d=Aedge - x0p,x1p,y0p,y1p=Bedge - self._Dirty[:,:,x0d:x1d,y0d:y1d]-=LocalSM[:,:,x0p:x1p,y0p:y1p] - W=np.float32(self.DicoDirty["WeightChansImages"]) - self._MeanDirty[0,:,x0d:x1d,y0d:y1d]-=np.sum(LocalSM[:,:,x0p:x1p,y0p:y1p]*W.reshape((W.size,1,1,1)),axis=0) - - -#=============================================== -#=============================================== -#=============================================== -#=============================================== - -class WorkerDeconvIsland(multiprocessing.Process): - def __init__(self, - work_queue, - result_queue, - GD, - IdSharedMem=None, - FreqsInfo=None, - MultiFreqMode=False): - multiprocessing.Process.__init__(self) - self.MultiFreqMode=MultiFreqMode - self.work_queue = work_queue - self.result_queue = result_queue - self.kill_received = False - self.exit = multiprocessing.Event() - self.GD=GD - self.IdSharedMem=IdSharedMem - self.FreqsInfo=FreqsInfo - self.CubeVariablePSF=NpShared.GiveArray("%s.CubeVariablePSF"%self.IdSharedMem) - self._Dirty=NpShared.GiveArray("%s.Dirty.ImagData"%self.IdSharedMem) - #self.WeightFreqBands=WeightFreqBands - - def shutdown(self): - self.exit.set() - - - def run(self): - from DDFacet.Imager.GA.ClassEvolveGA import ClassEvolveGA - while not self.kill_received: - #gc.enable() - try: - iIsland,FacetID,JonesNorm = self.work_queue.get() - except: - break - - - # iIsland=DicoOrder["iIsland"] - # FacetID=DicoOrder["FacetID"] - - # JonesNorm=DicoOrder["JonesNorm"] - - SharedListIsland="%s.ListIslands"%(self.IdSharedMem) - ThisPixList=NpShared.UnPackListArray(SharedListIsland)[iIsland].tolist() - - SharedBestIndiv="%s.ListBestIndiv"%(self.IdSharedMem) - IslandBestIndiv=NpShared.UnPackListArray(SharedBestIndiv)[iIsland] - - PSF=self.CubeVariablePSF[FacetID] - NGen=self.GD["GAClean"]["NMaxGen"] - NIndiv=self.GD["GAClean"]["NSourceKin"] - - ListPixParms=ThisPixList - ListPixData=ThisPixList - dx=self.GD["GAClean"]["NEnlargeData"] - if dx>0: - IncreaseIslandMachine=ClassIncreaseIsland.ClassIncreaseIsland() - ListPixData=IncreaseIslandMachine.IncreaseIsland(ListPixData,dx=dx) - - - # if island lies inside image - try: - nch=self.FreqsInfo["MeanJonesBand"][FacetID].size - WeightMeanJonesBand=self.FreqsInfo["MeanJonesBand"][FacetID].reshape((nch,1,1,1)) - WeightMueller=WeightMeanJonesBand.ravel() - WeightMuellerSignal=np.sqrt(WeightMueller*self.FreqsInfo["WeightChansImages"].ravel()) - - CEv=ClassEvolveGA(self._Dirty, - PSF, - self.FreqsInfo, - ListPixParms=ListPixParms, - ListPixData=ListPixData, - IslandBestIndiv=IslandBestIndiv,#*np.sqrt(JonesNorm), - GD=self.GD, - WeightFreqBands=WeightMuellerSignal) - Model=CEv.main(NGen=NGen,NIndiv=NIndiv,DoPlot=False) - - Model=np.array(Model).copy()#/np.sqrt(JonesNorm) - #Model*=CEv.ArrayMethodsMachine.Gain - - del(CEv) - - NpShared.ToShared("%s.FitIsland_%5.5i"%(self.IdSharedMem,iIsland),Model) - - #print "Current process: %s [%s left]"%(str(multiprocessing.current_process()),str(self.work_queue.qsize())) - - self.result_queue.put({"Success":True,"iIsland":iIsland}) - except Exception,e: - print "Exception : %s"%str(e) - - self.result_queue.put({"Success":False}) - diff --git a/DDFacet/Imager/MORESANE/ClassModelMachineMORESANE.py b/DDFacet/Imager/MORESANE/ClassModelMachineMORESANE.py index 2736246e..f5c1a3c0 100644 --- a/DDFacet/Imager/MORESANE/ClassModelMachineMORESANE.py +++ b/DDFacet/Imager/MORESANE/ClassModelMachineMORESANE.py @@ -17,12 +17,11 @@ along with this program; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. ''' - import numpy as np from DDFacet.Other import MyLogger from DDFacet.Other import ClassTimeIt from DDFacet.Other import ModColor -log=MyLogger.getLogger("ClassModelMachineMoresane") +log=MyLogger.getLogger("ClassModelMachineSSD") from DDFacet.Array import NpParallel from DDFacet.Array import ModLinAlg from DDFacet.ToolsDir import ModFFTW @@ -30,7 +29,7 @@ from DDFacet.Other import ClassTimeIt from DDFacet.Other import MyPickle from DDFacet.Other import reformat - +from DDFacet.Imager import ClassFrequencyMachine from DDFacet.ToolsDir.GiveEdges import GiveEdges from DDFacet.Imager import ClassModelMachine as ClassModelMachinebase from DDFacet.ToolsDir import ModFFTW @@ -39,67 +38,63 @@ from pyrap.images import image from SkyModel.Sky import ClassSM import os +import copy class ClassModelMachine(ClassModelMachinebase.ClassModelMachine): def __init__(self,*args,**kwargs): ClassModelMachinebase.ClassModelMachine.__init__(self, *args, **kwargs) - # self.GD=GD - # if Gain is None: - # self.Gain=self.GD["Deconv"]["Gain"] - # else: - # self.Gain=Gain - # self.GainMachine=GainMachine - # self.DicoSMStacked={} - # self.DicoSMStacked["Comp"]={} - if self.GD is not None: - self.SolveParam = self.GD["MORESANE"]["MOSolvePars"] - print>>log,"Solved parameters: %s"%(str(self.SolveParam)) - self.NParam=len(self.SolveParam) - - - def setRefFreq(self,RefFreq): # ,AllFreqs): + self.RefFreq=None + self.DicoModel={} + self.DicoModel["Type"]="MORESANE" + + def setRefFreq(self,RefFreq,Force=False):#,AllFreqs): + if self.RefFreq is not None and not Force: + print>>log,ModColor.Str("Reference frequency already set to %f MHz"%(self.RefFreq/1e6)) + return self.RefFreq=RefFreq - self.DicoSMStacked["RefFreq"]=RefFreq -# self.DicoSMStacked["AllFreqs"]=np.array(AllFreqs) - # print "ModelMachine:",self.RefFreq, self.DicoSMStacked["RefFreq"], self.DicoSMStacked["AllFreqs"] - - + self.DicoModel["RefFreq"]=RefFreq + + def setFreqMachine(self,GridFreqs, DegridFreqs): + # Initialise the Frequency Machine + self.FreqMachine = ClassFrequencyMachine.ClassFrequencyMachine(GridFreqs, DegridFreqs, self.DicoModel["RefFreq"], self.GD) def ToFile(self,FileName,DicoIn=None): print>>log, "Saving dico model to %s"%FileName if DicoIn is None: - D=self.DicoSMStacked + D=self.DicoModel else: D=DicoIn - #D["PM"]=self.PM + D["GD"]=self.GD D["ModelShape"]=self.ModelShape - D["Type"]="GA" # WILL CHANGE WHEN REFACTORING MODELMACHINE - D["SolveParam"]=self.SolveParam + D["Type"]="MORESANE" MyPickle.Save(D,FileName) + def giveDico(self): + D=self.DicoModel + D["GD"]=self.GD + D["ModelShape"]=self.ModelShape + D["Type"]="MORESANE" + return D def FromFile(self,FileName): - print>>log, "Reading dico model from %s"%FileName - self.DicoSMStacked=MyPickle.Load(FileName) - self.FromDico(self.DicoSMStacked) + print>>log, "Reading dico model from file %s"%FileName + self.DicoModel=MyPickle.Load(FileName) + self.FromDico(self.DicoModel) - def FromDico(self,DicoSMStacked): - #self.PM=self.DicoSMStacked["PM"] - self.DicoSMStacked=DicoSMStacked - self.RefFreq=self.DicoSMStacked["RefFreq"] - self.ModelShape=self.DicoSMStacked["ModelShape"] - try: - self.SolveParam=self.DicoSMStacked["SolveParam"] - except: - print>>log, "SolveParam is not in the keyword lists of DicoSMStacked" - print>>log, " setting SolveParam to [S]" - self.SolveParam=["S"] + def FromDico(self,DicoModel): + print>>log, "Reading dico model from dico with %i components"%len(DicoModel["Comp"]) + #self.PM=self.DicoModel["PM"] + self.DicoModel=DicoModel + self.RefFreq=self.DicoModel["RefFreq"] + self.ModelShape=self.DicoModel["ModelShape"] - self.NParam=len(self.SolveParam) + + + def setModelShape(self,ModelShape): self.ModelShape=ModelShape @@ -107,124 +102,49 @@ def setThreshold(self,Th): self.Th=Th - def GiveIndividual(self,ListPixParms): - NParms=self.NParam - OutArr=np.zeros((NParms,len(ListPixParms)),np.float32) - DicoComp=self.DicoSMStacked["Comp"] - - for iPix in range(len(ListPixParms)): - x,y=ListPixParms[iPix] - - xy=x,y - try: - Vals=DicoComp[xy]["Vals"][0] - OutArr[:,iPix]=Vals[:] - del(DicoComp[xy]) - except: - pass - - return OutArr.flatten() - - - def AppendIsland(self,ListPixParms,V,JonesNorm=None): - ListPix=ListPixParms - Vr=V.reshape((self.NParam,V.size/self.NParam)) - NPixListParms=len(ListPixParms) - - - #S=self.PM.ArrayToSubArray(V,Type="S") - - #S[np.abs(S)>log,((x,y),iSol,ThisSol) - - iS=np.where(SolveParam=="S")[0] - S=ThisSol[iS] - - ThisAlpha=0 - iAlpha=np.where(SolveParam=="Alpha")[0] - if iAlpha.size!=0: - ThisAlpha=ThisSol[iAlpha] - - for ch in range(nchan): - Flux=S*(FreqIn[ch]/RefFreq)**(ThisAlpha) - - for pol in range(npol): - ModelImage[ch,pol,x,y]+=Flux - - # vmin,vmax=np.min(self._MeanDirtyOrig[0,0]),np.max(self._MeanDirtyOrig[0,0]) - # vmin,vmax=-1,1 - # #vmin,vmax=np.min(ModelImage),np.max(ModelImage) - # pylab.clf() - # ax=pylab.subplot(1,3,1) - # pylab.imshow(self._MeanDirtyOrig[0,0],interpolation="nearest",vmin=vmin,vmax=vmax) - # pylab.subplot(1,3,2,sharex=ax,sharey=ax) - # pylab.imshow(self._MeanDirty[0,0],interpolation="nearest",vmin=vmin,vmax=vmax) - # pylab.colorbar() - # pylab.subplot(1,3,3,sharex=ax,sharey=ax) - # pylab.imshow( ModelImage[0,0],interpolation="nearest",vmin=vmin,vmax=vmax) - # pylab.colorbar() - # pylab.draw() - # pylab.show(False) - # print np.max(ModelImage[0,0]) - # # stop + if out is not None: + if out.shape != (nchan,npol,nx,ny) or out.dtype != np.float32: + raise RuntimeError("supplied image has incorrect type (%s) or shape (%s)" % (out.dtype, out.shape)) + ModelImage = out + else: + ModelImage = np.zeros((nchan,npol,nx,ny),dtype=np.float32) + + if 0 in self.DicoModel.keys(): + C0=self.DicoModel[0] + else: + C0=0 + if 1 in self.DicoModel.keys(): + C1=self.DicoModel[1] + else: + C1=0 + + ModelImage[:,:,:,:]=C0*(FreqIn.reshape((-1,1,1,1))/self.RefFreq)**C1 return ModelImage @@ -233,168 +153,20 @@ def GiveModelImage(self,FreqIn=None): def setListComponants(self,ListScales): self.ListScales=ListScales + def GiveSpectralIndexMap(self, threshold=0.1, save_dict=True): + # Get the model image + IM = self.GiveModelImage(self.FreqMachine.Freqsp) + nchan, npol, Nx, Ny = IM.shape + # Fit the alpha map + self.FreqMachine.FitAlphaMap(IM[:, 0, :, :], threshold=threshold) # should set threshold based on SNR of final residual + if save_dict: + FileName = self.GD['Output']['Name'] + ".Dicoalpha" + print>>log, "Saving componentwise SPI map to %s"%FileName - def GiveSpectralIndexMap(self,CellSizeRad=1.,GaussPars=[(1,1,0)],DoConv=True): - - - dFreq=1e6 - f0=self.DicoSMStacked["AllFreqs"].min() - f1=self.DicoSMStacked["AllFreqs"].max() - M0=self.GiveModelImage(f0) - M1=self.GiveModelImage(f1) - if DoConv: - M0=ModFFTW.ConvolveGaussian(M0,CellSizeRad=CellSizeRad,GaussPars=GaussPars) - M1=ModFFTW.ConvolveGaussian(M1,CellSizeRad=CellSizeRad,GaussPars=GaussPars) - - Np=1000 - indx,indy=np.int64(np.random.rand(Np)*M0.shape[0]),np.int64(np.random.rand(Np)*M0.shape[1]) - med=np.median(np.abs(M0[:,:,indx,indy])) - - Mask=((M1>100*med)&(M0>100*med)) - alpha=np.zeros_like(M0) - alpha[Mask]=(np.log(M0[Mask])-np.log(M1[Mask]))/(np.log(f0/f1)) - return alpha - - - def CleanNegComponants(self,box=20,sig=3,RemoveNeg=True): - print>>log, "Cleaning model dictionary from negative components with (box, sig) = (%i, %i)"%(box,sig) - ModelImage=self.GiveModelImage(self.DicoSMStacked["RefFreq"])[0,0] - - Min=scipy.ndimage.filters.minimum_filter(ModelImage,(box,box)) - Min[Min>0]=0 - Min=-Min - - if RemoveNeg==False: - Lx,Ly=np.where((ModelImage>log, " Removing neg components too" - Lx,Ly=np.where( ((ModelImage>log, " Component at (%i, %i) not in dict "%key - - def CleanMaskedComponants(self,MaskName): - print>>log, "Cleaning model dictionary from masked components using %s"%(MaskName) - im=image(MaskName) - MaskArray=im.getdata()[0,0].T[::-1] - for (x,y) in self.DicoSMStacked["Comp"].keys(): - if MaskArray[x,y]==0: - del(self.DicoSMStacked["Comp"][(x,y)]) - - def ToNPYModel(self,FitsFile,SkyModel,BeamImage=None): - #R=ModRegFile.RegToNp(PreCluster) - #R.Read() - #R.Cluster() - #PreClusterCat=R.CatSel - #ExcludeCat=R.CatExclude - - - AlphaMap=self.GiveSpectralIndexMap(DoConv=False) - ModelMap=self.GiveModelImage() - nch,npol,_,_=ModelMap.shape - - for ch in range(nch): - for pol in range(npol): - ModelMap[ch,pol]=ModelMap[ch,pol][::-1]#.T - AlphaMap[ch,pol]=AlphaMap[ch,pol][::-1]#.T - - if BeamImage is not None: - ModelMap*=(BeamImage) - + MyPickle.Save(self.FreqMachine.alpha_dict, FileName) - im=image(FitsFile) - pol,freq,decc,rac=im.toworld((0,0,0,0)) + return self.FreqMachine.weighted_alpha_map.reshape((1, 1, Nx, Ny)) - Lx,Ly=np.where(ModelMap[0,0]!=0) - - X=np.array(Lx) - Y=np.array(Ly) - - #pol,freq,decc1,rac1=im.toworld((0,0,1,0)) - dx=abs(im.coordinates().dict()["direction0"]["cdelt"][0]) - - SourceCat=np.zeros((X.size,),dtype=[('Name','|S200'),('ra',np.float),('dec',np.float),('Sref',np.float),('I',np.float),('Q',np.float),\ - ('U',np.float),('V',np.float),('RefFreq',np.float),('alpha',np.float),('ESref',np.float),\ - ('Ealpha',np.float),('kill',np.int),('Cluster',np.int),('Type',np.int),('Gmin',np.float),\ - ('Gmaj',np.float),('Gangle',np.float),("Select",np.int),('l',np.float),('m',np.float),("Exclude",bool), - ("X",np.int32),("Y",np.int32)]) - SourceCat=SourceCat.view(np.recarray) - - IndSource=0 - - SourceCat.RefFreq[:]=self.DicoSMStacked["RefFreq"] - _,_,nx,ny=ModelMap.shape - - for iSource in range(X.shape[0]): - x_iSource,y_iSource=X[iSource],Y[iSource] - _,_,dec_iSource,ra_iSource=im.toworld((0,0,y_iSource,x_iSource)) - SourceCat.ra[iSource]=ra_iSource - SourceCat.dec[iSource]=dec_iSource - SourceCat.X[iSource]=(nx-1)-X[iSource] - SourceCat.Y[iSource]=Y[iSource] - - #print self.DicoSMStacked["Comp"][(SourceCat.X[iSource],SourceCat.Y[iSource])] - # SourceCat.Cluster[IndSource]=iCluster - Flux=ModelMap[0,0,x_iSource,y_iSource] - Alpha=AlphaMap[0,0,x_iSource,y_iSource] - # print iSource,"/",X.shape[0],":",x_iSource,y_iSource,Flux,Alpha - SourceCat.I[iSource]=Flux - SourceCat.alpha[iSource]=Alpha - - - SourceCat=(SourceCat[SourceCat.ra!=0]).copy() - np.save(SkyModel,SourceCat) - self.AnalyticSourceCat=ClassSM.ClassSM(SkyModel) - - def DelAllComp(self): - for key in self.DicoSMStacked["Comp"].keys(): - del(self.DicoSMStacked["Comp"][key]) - - - def PutBackSubsComps(self): - #if self.GD["Data"]["RestoreDico"] is None: return - - SolsFile=self.GD["DDESolutions"]["DDSols"] - if not(".npz" in SolsFile): - Method=SolsFile - ThisMSName=reformat.reformat(os.path.abspath(self.GD["Data"]["MS"]),LastSlash=False) - SolsFile="%s/killMS.%s.sols.npz"%(ThisMSName,Method) - DicoSolsFile=np.load(SolsFile) - SourceCat=DicoSolsFile["SourceCatSub"] - SourceCat=SourceCat.view(np.recarray) - #RestoreDico=self.GD["Data"]["RestoreDico"] - RestoreDico=DicoSolsFile["ModelName"][()][0:-4]+".DicoModel" - - print>>log, "Adding previously subtracted components" - ModelMachine0=ClassModelMachine(self.GD) - - - ModelMachine0.FromFile(RestoreDico) - - - - _,_,nx0,ny0=ModelMachine0.DicoSMStacked["ModelShape"] - - _,_,nx1,ny1=self.ModelShape - dx=nx1-nx0 - - - for iSource in range(SourceCat.shape[0]): - x0=SourceCat.X[iSource] - y0=SourceCat.Y[iSource] - - x1=x0+dx - y1=y0+dx - - if not((x1,y1) in self.DicoSMStacked["Comp"].keys()): - self.DicoSMStacked["Comp"][(x1,y1)]=ModelMachine0.DicoSMStacked["Comp"][(x0,y0)] - else: - self.DicoSMStacked["Comp"][(x1,y1)]+=ModelMachine0.DicoSMStacked["Comp"][(x0,y0)] - diff --git a/DDFacet/Imager/MORESANE/ClassMoresane.py b/DDFacet/Imager/MORESANE/ClassMoresane.py index 4052c1bc..8670fd4b 100644 --- a/DDFacet/Imager/MORESANE/ClassMoresane.py +++ b/DDFacet/Imager/MORESANE/ClassMoresane.py @@ -38,95 +38,29 @@ class ClassMoresane(FI): # inherits from FitsImage but overriding __init__ to get rid of FITS file processing - def __init__(self,Dirty,PSF,DictMoresaneParms,GD=None): - - # manage Dirty - self.Dirty=Dirty - self.PSF=PSF - self.DictMoresaneParms=DictMoresaneParms - self.InitMoresane(self.DictMoresaneParms) - - def InitMoresane(self,DictMoresaneParms): - # load Moresane parms - - self.singlerun=DictMoresaneParms['SingleRun'] - self.startscale=DictMoresaneParms['StartScale'] - self.stopscale=DictMoresaneParms['StopScale'] - self.subregion=DictMoresaneParms['SubRegion'] - self.scalecount=DictMoresaneParms['ScaleCount'] - self.sigmalevel=DictMoresaneParms['SigmaLevel'] - self.loopgain=DictMoresaneParms['LoopGain'] - self.tolerance=DictMoresaneParms['Tolerance'] - self.accuracy=DictMoresaneParms['Accuracy'] - self.majorloopmiter=DictMoresaneParms['MajorLoopMiter'] - self.minorloopmiter=DictMoresaneParms['MinorLoopMiter'] - self.allongpu=DictMoresaneParms['AllOnGpu'] - self.decommode=DictMoresaneParms['DecomMode'] - self.corecount=DictMoresaneParms['CoreCount'] - self.convdevice=DictMoresaneParms['ConvDevice'] - self.convmode=DictMoresaneParms['ConvMode'] - self.extractionmode=DictMoresaneParms['ExtractionMode'] - self.enforcepositivity=DictMoresaneParms['EnforcePositivity'] - self.edgesuppression=DictMoresaneParms['EdgeSuppression'] - self.edgeoffset=DictMoresaneParms['EdgeOffset'] - self.fluxthreshold=DictMoresaneParms['FluxThreshold'] - self.negcomp=DictMoresaneParms['NegComp'] - self.edgeexcl=DictMoresaneParms['EdgeExcl'] - self.intexcl=DictMoresaneParms['IntExcl'] - - # Init - mask_name= None - self.mask_name = mask_name - - if self.mask_name is not None: - self.mask = fits.open("{}".format(mask_name))[0].data + def __init__(self,dirty,psf,mask=None,GD=None): + self.dirty_data = dirty + self.psf_data = psf + + self.mask_name=None + if mask is not None: + self.mask_name="NumpyMask" + self.mask = mask self.mask = self.mask.reshape(self.mask.shape[-2], self.mask.shape[-1]) - self.mask = self.mask / np.max(self.mask) - self.mask = fftconvolve(self.mask, np.ones([5, 5]), mode="same") - self.mask = self.mask / np.max(self.mask) + self.mask = self.mask/np.max(self.mask) + self.mask = fftconvolve(self.mask,np.ones([5,5]),mode="same") + self.mask = self.mask/np.max(self.mask) - self.dirty_data=self.Dirty - self.psf_data=self.PSF self.dirty_data_shape = self.dirty_data.shape self.psf_data_shape = self.psf_data.shape self.complete = False - self.model = np.zeros_like(self.Dirty) - self.residual = np.copy(self.Dirty) - self.restored = np.zeros_like(self.Dirty) - - - def main(self): - - # Proper Moresane run - if self.singlerun: - print "Single run" - self.moresane(self.subregion, self.scalecount, self.sigmalevel, self.loopgain, self.tolerance, self.accuracy, - self.majorloopmiter, self.minorloopmiter, self.allongpu, self.decommode, self.corecount, - self.convdevice, self.convmode, self.extractionmode, self.enforcepositivity, - self.edgesuppression, self.edgeoffset,self.fluxthreshold, self.negcomp, self.edgeexcl, self.intexcl) - else: - print "By scale" - self.moresane_by_scale(self.startscale, self.stopscale, self.subregion, self.sigmalevel, self.loopgain, self.tolerance, self.accuracy, - self.majorloopmiter, self.minorloopmiter, self.allongpu, self.decommode, self.corecount, - self.convdevice, self.convmode, self.extractionmode, self.enforcepositivity, - self.edgesuppression, self.edgeoffset,self.fluxthreshold, self.negcomp, self.edgeexcl, self.intexcl) - return self.model,self.residual #IslandModel - - def residuals(self): - return self.residual - - def GiveCLEANBeam(self, PSF, cellsize): - cellsizeindeg = cellsize * 1. / 3600 # to convert arcsec in degrees - fakePSFFitsHeader = {"CDELT1": cellsizeindeg, "CDELT2": cellsizeindeg} - - clean_beam, beam_params = beam_fit(PSF, fakePSFFitsHeader) - return clean_beam, beam_params - - def GiveBeamArea(self, beam_params, cellsize): - print beam_params - cellsizedeg = cellsize * 1. / 3600 - bx, by, _ = beam_params # bx and by in degrees - px, py = cellsizedeg, cellsizedeg # in degrees - BeamArea = np.pi * bx * by / (4 * np.log(2) * px * py) - return BeamArea + self.model = np.zeros_like(self.dirty_data) + self.residual = np.copy(self.dirty_data) + + + + def giveModelResid(self,*args,**kwargs): + self.moresane(*args,**kwargs) + return self.model,self.residual + diff --git a/DDFacet/Imager/MORESANE/ClassMoresaneSingleSlice.py b/DDFacet/Imager/MORESANE/ClassMoresaneSingleSlice.py new file mode 100644 index 00000000..8c369873 --- /dev/null +++ b/DDFacet/Imager/MORESANE/ClassMoresaneSingleSlice.py @@ -0,0 +1,56 @@ +''' +DDFacet, a facet-based radio imaging package +Copyright (C) 2013-2016 Cyril Tasse, l'Observatoire de Paris, +SKA South Africa, Rhodes University + +This program is free software; you can redistribute it and/or +modify it under the terms of the GNU General Public License +as published by the Free Software Foundation; either version 2 +of the License, or (at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program; if not, write to the Free Software +Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +''' + +import numpy as np +import logging +import time + +from DDFacet.Other import MyLogger +from pymoresane.main import FitsImage as FI + +class ClassMoresaneSingleSlice(FI): + def __init__(self,dirty,psf,mask=None,GD=None): + self.dirty_data = dirty + self.psf_data = psf + + #MyLogger.setSilent(["pymoresane.main"]) + self.mask_name=None + if mask is not None: + self.mask_name="NumpyMask" + self.mask = mask + self.mask = self.mask.reshape(self.mask.shape[-2], self.mask.shape[-1]) + self.mask = self.mask/np.max(self.mask) + self.mask = fftconvolve(self.mask,np.ones([5,5]),mode="same") + self.mask = self.mask/np.max(self.mask) + + self.dirty_data_shape = self.dirty_data.shape + self.psf_data_shape = self.psf_data.shape + + self.complete = False + self.model = np.zeros_like(self.dirty_data) + self.residual = np.copy(self.dirty_data) + + + + def giveModelResid(self,*args,**kwargs): + with np.errstate(divide='ignore'): + self.moresane(*args,**kwargs) + return self.model,self.residual + diff --git a/DDFacet/Imager/MSMF/ClassImageDeconvMachineMSMF.py b/DDFacet/Imager/MSMF/ClassImageDeconvMachineMSMF.py index 40ba393d..0ab61075 100644 --- a/DDFacet/Imager/MSMF/ClassImageDeconvMachineMSMF.py +++ b/DDFacet/Imager/MSMF/ClassImageDeconvMachineMSMF.py @@ -537,6 +537,19 @@ def SubStep(self,(dx,dy),LocalSM): # #unc print Bedge # # print self.Dirty[0,x0d:x1d,y0d:y1d] + def Plot(self): + import pylab + pylab.clf() + pylab.subplot(1,3,1) + pylab.imshow(self._CubeDirty[0,0]) + pylab.colorbar() + pylab.subplot(1,3,2) + pylab.imshow(self._CubeDirty[1,0]) + pylab.colorbar() + pylab.draw() + pylab.show() + + def updateRMS(self): _,npol,npix,_ = self._MeanDirty.shape NPixStats = self.GD["Deconv"]["NumRMSSamples"] diff --git a/DDFacet/Imager/MSMF/ClassModelMachineMSMF.py b/DDFacet/Imager/MSMF/ClassModelMachineMSMF.py index 6e0ecac6..d393abae 100644 --- a/DDFacet/Imager/MSMF/ClassModelMachineMSMF.py +++ b/DDFacet/Imager/MSMF/ClassModelMachineMSMF.py @@ -158,65 +158,70 @@ def AppendComponentToDictStacked(self,key,Fpol,Sols,pol_array_index=0): def setListComponants(self,ListScales): self.ListScales=ListScales - def GiveSpectralIndexMap(self, threshold=0.1, save_dict=True): - # Get the model image - IM = self.GiveModelImage(self.FreqMachine.Freqsp) - nchan, npol, Nx, Ny = IM.shape - - # Fit the alpha map - self.FreqMachine.FitAlphaMap(IM[:, 0, :, :], - threshold=threshold) # should set threshold based on SNR of final residual - - if save_dict: - FileName = self.GD['Output']['Name'] + ".Dicoalpha" - print>> log, "Saving componentwise SPI map to %s" % FileName - - MyPickle.Save(self.FreqMachine.alpha_dict, FileName) - - return self.FreqMachine.weighted_alpha_map.reshape((1, 1, Nx, Ny)) - - # def GiveSpectralIndexMap(self,CellSizeRad=1.,GaussPars=[(1,1,0)],DoConv=True,MaxSpi=100,MaxDR=1e+6): - # dFreq=1e6 - # #f0=self.DicoSMStacked["AllFreqs"].min() - # #f1=self.DicoSMStacked["AllFreqs"].max() - # RefFreq=self.DicoSMStacked["RefFreq"] - # f0=RefFreq/1.5 - # f1=RefFreq*1.5 - # - # M0=self.GiveModelImage(f0) - # M1=self.GiveModelImage(f1) - # if DoConv: - # M0=ModFFTW.ConvolveGaussian(M0,CellSizeRad=CellSizeRad,GaussPars=GaussPars) - # M1=ModFFTW.ConvolveGaussian(M1,CellSizeRad=CellSizeRad,GaussPars=GaussPars) - # - # # compute threshold for alpha computation by rounding DR threshold to .1 digits (i.e. 1.65e-6 rounds to 1.7e-6) - # if not np.all(M0==0): - # minmod = float("%.1e"%(np.max(np.abs(M0))/MaxDR)) - # else: - # minmod=1e-6 - # - # # mask out pixels above threshold - # mask=(M1>log,"computing alpha map for model pixels above %.1e Jy (based on max DR setting of %g)"%(minmod,MaxDR) - # M0[mask]=minmod - # M1[mask]=minmod - # #with np.errstate(invalid='ignore'): - # # alpha = (np.log(M0)-np.log(M1))/(np.log(f0/f1)) - # # print - # # print np.min(M0),np.min(M1),minmod - # # print - # alpha = (np.log(M0)-np.log(M1))/(np.log(f0/f1)) - # alpha[mask] = 0 - # - # # mask out |alpha|>MaxSpi. These are not physically meaningful anyway - # mask = alpha>MaxSpi - # alpha[mask] = MaxSpi - # masked = mask.any() - # mask = alpha<-MaxSpi - # alpha[mask] = -MaxSpi - # if masked or mask.any(): - # print>>log,ModColor.Str("WARNING: some alpha pixels outside +/-%g. Masking them."%MaxSpi,col="red") - # return alpha + # def GiveSpectralIndexMap(self, threshold=0.1, save_dict=True): + # # Get the model image + # IM = self.GiveModelImage(self.FreqMachine.Freqsp) + # nchan, npol, Nx, Ny = IM.shape + # # Fit the alpha map + # self.FreqMachine.FitAlphaMap(IM[:, 0, :, :], + # threshold=threshold) # should set threshold based on SNR of final residual + # if save_dict: + # FileName = self.GD['Output']['Name'] + ".Dicoalpha" + # print>> log, "Saving componentwise SPI map to %s" % FileName + # MyPickle.Save(self.FreqMachine.alpha_dict, FileName) + # #return self.FreqMachine.weighted_alpha_map.reshape((1, 1, Nx, Ny)) + # return self.FreqMachine.alpha_map.reshape((1, 1, Nx, Ny)) + + def GiveSpectralIndexMap(self,CellSizeRad=1.,GaussPars=[(1,1,0)],DoConv=True,MaxSpi=100,MaxDR=1e+6,threshold=None): + dFreq=1e6 + #f0=self.DicoSMStacked["AllFreqs"].min() + #f1=self.DicoSMStacked["AllFreqs"].max() + RefFreq=self.DicoSMStacked["RefFreq"] + f0=RefFreq/1.5 + f1=RefFreq*1.5 + + M0=self.GiveModelImage(f0) + M1=self.GiveModelImage(f1) + if DoConv: + #M0=ModFFTW.ConvolveGaussian(M0,CellSizeRad=CellSizeRad,GaussPars=GaussPars) + #M1=ModFFTW.ConvolveGaussian(M1,CellSizeRad=CellSizeRad,GaussPars=GaussPars) + #M0,_=ModFFTW.ConvolveGaussianWrapper(M0,Sig=GaussPars[0][0]/CellSizeRad) + #M1,_=ModFFTW.ConvolveGaussianWrapper(M1,Sig=GaussPars[0][0]/CellSizeRad) + M0,_=ModFFTW.ConvolveGaussianScipy(M0,Sig=GaussPars[0][0]/CellSizeRad) + M1,_=ModFFTW.ConvolveGaussianScipy(M1,Sig=GaussPars[0][0]/CellSizeRad) + + + #print M0.shape,M1.shape + # compute threshold for alpha computation by rounding DR threshold to .1 digits (i.e. 1.65e-6 rounds to 1.7e-6) + if threshold is not None: + minmod = threshold + elif not np.all(M0==0): + minmod = float("%.1e"%(np.max(np.abs(M0))/MaxDR)) + else: + minmod=1e-6 + + # mask out pixels above threshold + mask=(M1>log,"computing alpha map for model pixels above %.1e Jy (based on max DR setting of %g)"%(minmod,MaxDR) + M0[mask]=minmod + M1[mask]=minmod + #with np.errstate(invalid='ignore'): + # alpha = (np.log(M0)-np.log(M1))/(np.log(f0/f1)) + # print + # print np.min(M0),np.min(M1),minmod + # print + alpha = (np.log(M0)-np.log(M1))/(np.log(f0/f1)) + alpha[mask] = 0 + + # mask out |alpha|>MaxSpi. These are not physically meaningful anyway + mask = alpha>MaxSpi + alpha[mask] = MaxSpi + masked = mask.any() + mask = alpha<-MaxSpi + alpha[mask] = -MaxSpi + if masked or mask.any(): + print>>log,ModColor.Str("WARNING: some alpha pixels outside +/-%g. Masking them."%MaxSpi,col="red") + return alpha def GiveModelList(self): """ diff --git a/DDFacet/Imager/MSMF/ClassMultiScaleMachine.py b/DDFacet/Imager/MSMF/ClassMultiScaleMachine.py index cbabd659..616a0b11 100644 --- a/DDFacet/Imager/MSMF/ClassMultiScaleMachine.py +++ b/DDFacet/Imager/MSMF/ClassMultiScaleMachine.py @@ -223,7 +223,10 @@ def SetDirty(self,DicoDirty): self.DicoDirty=DicoDirty #self.NChannels=self.DicoDirty["NChannels"] - + PeakSearchImage=self.DicoDirty["MeanImage"] + NPixStats = 1000 + IndStats=np.int64(np.linspace(0,PeakSearchImage.size-1,NPixStats)) + self.RMS=np.std(np.real(PeakSearchImage.ravel()[IndStats])) self._Dirty=self.DicoDirty["ImageCube"] self._MeanDirty=self.DicoDirty["MeanImage"] @@ -1014,7 +1017,11 @@ def GiveLocalSM(self,(x,y),Fpol): PeakMeanOrigDirty=MeanOrigDirty[xc0[0],yc0[0]] dirtyVec=dirtyVec.copy() Mask=np.zeros(WVecPSF.shape,np.bool8) - for iIter in range(10): + + T=ClassTimeIt.ClassTimeIt() + T.disable() + NNLSStep=10 + for iIter in range(1000): A=W*BM y=W*dirtyVec d=dirtyVec.reshape((nchan,1,nxp,nyp))[:,0] @@ -1022,17 +1029,23 @@ def GiveLocalSM(self,(x,y),Fpol): FactNorm=np.abs(PeakMeanOrigDirty/PeakMeanOrigResid) if np.isnan(FactNorm) or np.isinf(FactNorm): + #print "Cond1 %i"%iIter Sol=np.zeros((A.shape[1],),dtype=np.float32) break - + T.timeit("0") if 1.2*sig)&(Resid==MaxResid)) + _,xc1,yc1=np.where((Resid>self.GD["HMP"]["OuterSpaceTh"]*sig)&(Resid==MaxResid)) dirtyVecSub=d Sol=x @@ -1082,9 +1097,12 @@ def GiveLocalSM(self,(x,y),Fpol): # If source is contaminating, substract it with the delta (with alpha=0) + T.timeit("3") if xc1.size>0 and MaxResid>Peak/100.: CentralPixel=(xc1[0]==xc0[0] and yc1[0]==yc0[0]) - if CentralPixel: break + if CentralPixel: + #print "CondCentralPix %i"%iIter + break F=Resid[:,xc1[0],yc1[0]] dx,dy=nxp/2-xc1[0],nyp/2-yc1[0] _,_,nxPSF,nyPSF=self.SubPSF.shape @@ -1099,25 +1117,17 @@ def GiveLocalSM(self,(x,y),Fpol): ThisPSF=self.SubPSF[:,0,x0p:x1p,y0p:y1p] _,nxThisPSF,nyThisPSF=ThisPSF.shape - # # find the optimal flux value for the two cross contaminating sources case - # al=np.abs(ThisPSF[:,nxThisPSF/2,nyThisPSF/2]) - # MeanAl=np.mean(al) - # if 0.010, MaxResid>Peak/100.: ",xc1.size>0, MaxResid>Peak/100. + + x,_=scipy.optimize.nnls(A, y.ravel()) DoBreak=True # ####### debug @@ -1140,13 +1150,14 @@ def GiveLocalSM(self,(x,y),Fpol): # pylab.colorbar() # pylab.subplot(2,3,5) # pylab.imshow(dirtyVecSub[0],interpolation="nearest") - # pylab.title("NewResid") + # pylab.title("NewDirty") # pylab.colorbar() # pylab.draw() # pylab.show(False) # pylab.pause(0.1) # ##################### + T.timeit("4") if DoBreak: break diff --git a/DDFacet/Imager/MUFFIN/ClassImageDeconvMachineMUFFIN.py b/DDFacet/Imager/MUFFIN/ClassImageDeconvMachineMUFFIN.py new file mode 100644 index 00000000..1e492860 --- /dev/null +++ b/DDFacet/Imager/MUFFIN/ClassImageDeconvMachineMUFFIN.py @@ -0,0 +1,303 @@ +''' +DDFacet, a facet-based radio imaging package +Copyright (C) 2013-2016 Cyril Tasse, l'Observatoire de Paris, +SKA South Africa, Rhodes University + +This program is free software; you can redistribute it and/or +modify it under the terms of the GNU General Public License +as published by the Free Software Foundation; either version 2 +of the License, or (at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program; if not, write to the Free Software +Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +''' + +import os +import numpy as np +from DDFacet.Other import MyLogger +from DDFacet.Other import ModColor +log=MyLogger.getLogger("ClassImageDeconvMachine") +from DDFacet.Array import NpParallel +from DDFacet.Array import NpShared +from DDFacet.ToolsDir import ModFFTW +from DDFacet.ToolsDir import ModToolBox +from DDFacet.Other import ClassTimeIt +from pyrap.images import image +from DDFacet.Imager.ClassPSFServer import ClassPSFServer +from DDFacet.Other.progressbar import ProgressBar +from DDFacet.Imager import ClassGainMachine +from DDFacet.Other import MyPickle +import multiprocessing +import time +#from DDFacet.Imager.MORESANE.ClassMoresaneSingleSlice import ClassMoresaneSingleSlice +##from DDFacet.Imager.MUFFIN.easy_muffin.easy_muffin_py.deconv3d import EasyMuffin +from easy_muffin_py.deconv3d import EasyMuffin +from DDFacet.Array import shared_dict +from DDFacet.ToolsDir import ClassSpectralFunctions +from scipy.optimize import least_squares +from DDFacet.ToolsDir.GiveEdges import GiveEdges + + +class ClassImageDeconvMachine(): + def __init__(self,GD=None,ModelMachine=None,RefFreq=None,*args,**kw): + self.GD=GD + self.ModelMachine = ModelMachine + self.RefFreq=RefFreq + if self.ModelMachine.DicoModel["Type"]!="MUFFIN": + raise ValueError("ModelMachine Type should be MUFFIN") + self.MultiFreqMode=(self.GD["Freq"]["NBand"]>1) + + def SetPSF(self,DicoVariablePSF): + self.PSFServer=ClassPSFServer(self.GD) + DicoVariablePSF=shared_dict.attach(DicoVariablePSF.path)#["CubeVariablePSF"] + self.PSFServer.setDicoVariablePSF(DicoVariablePSF) + self.PSFServer.setRefFreq(self.ModelMachine.RefFreq) + self.DicoVariablePSF=DicoVariablePSF + self.setFreqs(self.PSFServer.DicoMappingDesc) + + def setMaskMachine(self,MaskMachine): + self.MaskMachine=MaskMachine + + def setFreqs(self,DicoMappingDesc): + self.DicoMappingDesc=DicoMappingDesc + if self.DicoMappingDesc is None: return + self.SpectralFunctionsMachine=ClassSpectralFunctions.ClassSpectralFunctions(self.DicoMappingDesc,RefFreq=self.DicoMappingDesc["RefFreq"])#,BeamEnable=False) + self.SpectralFunctionsMachine.CalcFluxBands() + + def GiveModelImage(self,*args): return self.ModelMachine.GiveModelImage(*args) + + def Update(self,DicoDirty,**kwargs): + """ + Method to update attributes from ClassDeconvMachine + """ + #Update image dict + self.SetDirty(DicoDirty) + + def ToFile(self, fname): + """ + Write model dict to file + """ + self.ModelMachine.ToFile(fname) + + def FromFile(self, fname): + """ + Read model dict from file SubtractModel + """ + self.ModelMachine.FromFile(fname) + + def FromDico(self, DicoName): + """ + Read in model dict + """ + self.ModelMachine.FromDico(DicoName) + + def setSideLobeLevel(self,SideLobeLevel,OffsetSideLobe): + self.SideLobeLevel=SideLobeLevel + self.OffsetSideLobe=OffsetSideLobe + + def Init(self,**kwargs): + self.SetPSF(kwargs["PSFVar"]) + if "PSFSideLobes" not in self.DicoVariablePSF.keys(): + self.DicoVariablePSF["PSFSideLobes"]=kwargs["PSFAve"] + self.setSideLobeLevel(kwargs["PSFAve"][0], kwargs["PSFAve"][1]) + self.ModelMachine.setRefFreq(kwargs["RefFreq"]) + # store grid and degrid freqs for ease of passing to MSMF + #print kwargs["GridFreqs"],kwargs["DegridFreqs"] + self.GridFreqs=kwargs["GridFreqs"] + self.DegridFreqs=kwargs["DegridFreqs"] + self.ModelMachine.setFreqMachine(kwargs["GridFreqs"], kwargs["DegridFreqs"]) + + def SetDirty(self,DicoDirty): + self.DicoDirty=DicoDirty + self._Dirty=self.DicoDirty["ImageCube"] + self._MeanDirty=self.DicoDirty["MeanImage"] + NPSF=self.PSFServer.NPSF + _,_,NDirty,_=self._Dirty.shape + off=(NPSF-NDirty)/2 + self.DirtyExtent=(off,off+NDirty,off,off+NDirty) + self.ModelMachine.setModelShape(self._Dirty.shape) + + def AdaptArrayShape(self,A,Nout): + nch,npol,Nin,_=A.shape + if Nin==Nout: + return A + elif Nin>Nout: + # dx=Nout/2 + # B=np.zeros((nch,npol,Nout,Nout),A.dtype) + # print>>log," Adapt shapes: %s -> %s"%(str(A.shape),str(B.shape)) + # B[:]=A[...,Nin/2-dx:Nin/2+dx+1,Nin/2-dx:Nin/2+dx+1] + + N0=A.shape[-1] + xc0=yc0=N0/2 + N1=Nout + xc1=yc1=N1/2 + Aedge,Bedge=GiveEdges((xc0,yc0),N0,(xc1,yc1),N1) + x0d,x1d,y0d,y1d=Aedge + x0p,x1p,y0p,y1p=Bedge + B=A[...,x0d:x1d,y0d:y1d] + + return B + else: + stop + return None + + def updateModelMachine(self,ModelMachine): + self.ModelMachine=ModelMachine + if self.ModelMachine.RefFreq!=self.RefFreq: + raise ValueError("freqs should be equal") + + def updateMask(self,Mask): + nx,ny=Mask.shape + self._MaskArray = np.zeros((1,1,nx,ny),np.bool8) + self._MaskArray[0,0,:,:]=Mask[:,:] + + + def Deconvolve(self): + + if self._Dirty.shape[-1]!=self._Dirty.shape[-2]: + # print "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" + # print self._Dirty.shape + # print "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" + return "MaxIter", True, True + + dirty=self._Dirty + nch,npol,nx,ny=dirty.shape + Model=np.zeros_like(dirty) + + _,_,xp,yp=np.where(self._MeanDirty==np.max(self._MeanDirty)) + self.PSFServer.setLocation(xp,yp) + self.iFacet=self.PSFServer.iFacet + psf,_=self.PSFServer.GivePSF() + nxPSF=psf.shape[-1] + nxDirty=dirty.shape[-1] + + Nout=np.min([dirty.shape[-1],psf.shape[-1]]) + dirty=self.AdaptArrayShape(dirty,Nout) + SliceDirty=slice(0,None) + if dirty.shape[-1]%2!=0: + SliceDirty=slice(0,-1) + + d=dirty[:,:,SliceDirty,SliceDirty] + psf=self.AdaptArrayShape(psf,d.shape[-1]) + + SlicePSF=slice(0,None) + if psf.shape[-1]%2!=0: + SlicePSF=slice(0,-1) + + p=psf[:,:,SlicePSF,SlicePSF] + + dirty_MUFFIN = np.squeeze(d[:,0,:,:]) + dirty_MUFFIN = dirty_MUFFIN.transpose((2,1,0)) + + psf_MUFFIN = np.squeeze(p[:,0,:,:]) + psf_MUFFIN = psf_MUFFIN.transpose((2,1,0)) + + EM = EasyMuffin(mu_s=self.GD['MUFFIN']['mu_s'], + mu_l=self.GD['MUFFIN']['mu_l'], + nb=self.GD['MUFFIN']['nb'], + truesky=dirty_MUFFIN, + psf=psf_MUFFIN, + dirty=dirty_MUFFIN) + EM.loop(nitermax=self.GD['MUFFIN']['NMinorIter']) + + + nxModel=dirty_MUFFIN.shape[0] + Aedge,Bedge=GiveEdges((nxModel/2,nxModel/2),nxModel,(nxDirty/2,nxDirty/2),nxDirty) + x0,x1,y0,y1=Bedge + + Model=np.zeros((nxDirty,nxDirty,nch)) + Model[x0:x1,y0:y1,:]=EM.x + self.ModelMachine.setMUFFINModel(Model) + + # if self._Dirty.shape[-1]!=self._Dirty.shape[-2]: + # # print "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" + # # print self._Dirty.shape + # # print "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" + # return "MaxIter", True, True + + + # dirty=self._Dirty + # nch,npol,_,_=dirty.shape + # Model=np.zeros_like(dirty) + + # _,_,xp,yp=np.where(self._MeanDirty==np.max(self._MeanDirty)) + # self.PSFServer.setLocation(xp,yp) + # self.iFacet=self.PSFServer.iFacet + + # psf,_=self.PSFServer.GivePSF() + + # Nout=np.min([dirty.shape[-1],psf.shape[-1]]) + # dirty=self.AdaptArrayShape(dirty,Nout) + # SliceDirty=slice(0,None) + # if dirty.shape[-1]%2!=0: + # SliceDirty=slice(0,-1) + + # d=dirty[:,:,SliceDirty,SliceDirty] + # psf=self.AdaptArrayShape(psf,d.shape[-1]*2) + + # SlicePSF=slice(0,None) + # if psf.shape[-1]%2!=0: + # SlicePSF=slice(0,-1) + + # p=psf[:,:,SlicePSF,SlicePSF] + # if p.shape[-1]!=2*d.shape[-1]: + # print "!!!!!!!!!!!!!!!!!!!!!!!!!" + # print "Could not adapt psf shape to 2*dirty shape!!!!!!!!!!!!!!!!!!!!!!!!!" + # print p.shape[-1],d.shape[-1] + # print "!!!!!!!!!!!!!!!!!!!!!!!!!" + # psf=self.AdaptArrayShape(psf,d.shape[-1]) + # SlicePSF=SliceDirty + + # for ch in range(nch): + # CM=ClassMoresaneSingleSlice(dirty[ch,0,SliceDirty,SliceDirty],psf[ch,0,SlicePSF,SlicePSF],mask=None,GD=None) + # model,resid=CM.giveModelResid(major_loop_miter=self.GD["MORESANE"]["NMajorIter"], + # minor_loop_miter=self.GD["MORESANE"]["NMinorIter"], + # loop_gain=self.GD["MORESANE"]["Gain"], + # sigma_level=self.GD["MORESANE"]["SigmaCutLevel"],# tolerance=1., + # enforce_positivity=self.GD["MORESANE"]["ForcePositive"]) + # Model[ch,0,SliceDirty,SliceDirty]=model[:,:] + + # import pylab + # pylab.clf() + # pylab.subplot(2,2,1) + # pylab.imshow(dirty[ch,0,SliceDirty,SliceDirty],interpolation="nearest") + # pylab.colorbar() + + # pylab.subplot(2,2,2) + # pylab.imshow(psf[ch,0,SlicePSF,SlicePSF],interpolation="nearest") + # pylab.colorbar() + + # pylab.subplot(2,2,3) + # pylab.imshow(model,interpolation="nearest") + # pylab.colorbar() + + # pylab.subplot(2,2,4) + # pylab.imshow(resid,interpolation="nearest") + # pylab.colorbar() + + # pylab.draw() + # pylab.show() + + + # print + # print np.max(np.max(Model,axis=-1),axis=-1) + # print + # print + + + + + #_,_,nx,ny=Model.shape + #Model=np.mean(Model,axis=0).reshape((1,1,nx,ny)) + + #Model.fill(0) + #Model[:,:,xp,yp]=self._Dirty[:,:,xp,yp] + + + return "MaxIter", True, True # stop deconvolution but do update model diff --git a/DDFacet/Imager/MUFFIN/ClassModelMachineMUFFIN.py b/DDFacet/Imager/MUFFIN/ClassModelMachineMUFFIN.py new file mode 100644 index 00000000..c0063c23 --- /dev/null +++ b/DDFacet/Imager/MUFFIN/ClassModelMachineMUFFIN.py @@ -0,0 +1,145 @@ +''' +DDFacet, a facet-based radio imaging package +Copyright (C) 2013-2016 Cyril Tasse, l'Observatoire de Paris, +SKA South Africa, Rhodes University + +This program is free software; you can redistribute it and/or +modify it under the terms of the GNU General Public License +as published by the Free Software Foundation; either version 2 +of the License, or (at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program; if not, write to the Free Software +Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +''' +import numpy as np +from DDFacet.Other import MyLogger +from DDFacet.Other import ClassTimeIt +from DDFacet.Other import ModColor +log=MyLogger.getLogger("ClassModelMachineSSD") +from DDFacet.Array import NpParallel +from DDFacet.Array import ModLinAlg +from DDFacet.ToolsDir import ModFFTW +from DDFacet.ToolsDir import ModToolBox +from DDFacet.Other import ClassTimeIt +from DDFacet.Other import MyPickle +from DDFacet.Other import reformat +from DDFacet.Imager import ClassFrequencyMachine +from DDFacet.ToolsDir.GiveEdges import GiveEdges +from DDFacet.Imager import ClassModelMachine as ClassModelMachinebase +from DDFacet.ToolsDir import ModFFTW +import scipy.ndimage +from SkyModel.Sky import ModRegFile +from pyrap.images import image +from SkyModel.Sky import ClassSM +import os +import copy + +class ClassModelMachine(ClassModelMachinebase.ClassModelMachine): + def __init__(self,*args,**kwargs): + ClassModelMachinebase.ClassModelMachine.__init__(self, *args, **kwargs) + self.RefFreq=None + self.DicoModel={} + self.DicoModel["Type"]="MUFFIN" + self.ModelImageCube=0 + + def setRefFreq(self,RefFreq,Force=False):#,AllFreqs): + if self.RefFreq is not None and not Force: + print>>log,ModColor.Str("Reference frequency already set to %f MHz"%(self.RefFreq/1e6)) + return + self.RefFreq=RefFreq + self.DicoModel["RefFreq"]=RefFreq + + def setFreqMachine(self,GridFreqs, DegridFreqs): + # Initialise the Frequency Machine + self.FreqMachine = ClassFrequencyMachine.ClassFrequencyMachine(GridFreqs, DegridFreqs, self.DicoModel["RefFreq"], self.GD) + self.GridFreqs=GridFreqs + + def ToFile(self,FileName,DicoIn=None): + print>>log, "Saving dico model to %s"%FileName + if DicoIn is None: + D=self.DicoModel + else: + D=DicoIn + + D["GD"]=self.GD + D["ModelShape"]=self.ModelShape + D["FreqsCube"]=self.GridFreqs + D["Type"]="MUFFIN" + + MyPickle.Save(D,FileName) + + def giveDico(self): + D=self.DicoModel + D["GD"]=self.GD + D["ModelShape"]=self.ModelShape + D["Type"]="MUFFIN" + return D + + def FromFile(self,FileName): + print>>log, "Reading dico model from file %s"%FileName + self.DicoModel=MyPickle.Load(FileName) + self.FromDico(self.DicoModel) + + + def FromDico(self,DicoModel): + print>>log, "Reading dico model from dico with %i components"%len(DicoModel["Comp"]) + #self.PM=self.DicoModel["PM"] + self.DicoModel=DicoModel + self.RefFreq=self.DicoModel["RefFreq"] + self.ModelShape=self.DicoModel["ModelShape"] + + + def setModelShape(self,ModelShape): + self.ModelShape=ModelShape + + def setThreshold(self,Th): + self.Th=Th + + def setMUFFINModel(self,ModelImageCube): + self.ModelImageCube += ModelImageCube + + def GiveModelImage(self,FreqIn=None): + + FreqIn = np.array([FreqIn.ravel()]).flatten() + + iFreqSliceMuffinCube = np.argmin(np.abs(FreqIn.reshape(-1,1)-self.GridFreqs.reshape(1,-1)),axis=1) + + nchan = FreqIn.size + _,npol,nx,ny=self.ModelShape + ModelImage = np.zeros((nchan,npol,nx,ny),dtype=np.float32) + + ModelImage_ = self.ModelImageCube[:,:,iFreqSliceMuffinCube] + ModelImage_ = ModelImage_.transpose(2,1,0) + + ModelImage[:,0,:,:] = ModelImage_ + + return ModelImage + + + + def setListComponants(self,ListScales): + self.ListScales=ListScales + + def GiveSpectralIndexMap(self, threshold=0.1, save_dict=True): + # Get the model image + IM = self.GiveModelImage(self.FreqMachine.Freqsp) + nchan, npol, Nx, Ny = IM.shape + + # Fit the alpha map + self.FreqMachine.FitAlphaMap(IM[:, 0, :, :], threshold=threshold) # should set threshold based on SNR of final residual + + if save_dict: + FileName = self.GD['Output']['Name'] + ".Dicoalpha" + print>>log, "Saving componentwise SPI map to %s"%FileName + + MyPickle.Save(self.FreqMachine.alpha_dict, FileName) + + return self.FreqMachine.weighted_alpha_map.reshape((1, 1, Nx, Ny)) + + diff --git a/DDFacet/Imager/MUFFIN/__init__.py b/DDFacet/Imager/MUFFIN/__init__.py new file mode 100644 index 00000000..6c677ea9 --- /dev/null +++ b/DDFacet/Imager/MUFFIN/__init__.py @@ -0,0 +1,19 @@ +''' +DDFacet, a facet-based radio imaging package +Copyright (C) 2013-2016 Cyril Tasse, l'Observatoire de Paris, +SKA South Africa, Rhodes University + +This program is free software; you can redistribute it and/or +modify it under the terms of the GNU General Public License +as published by the Free Software Foundation; either version 2 +of the License, or (at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program; if not, write to the Free Software +Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +''' \ No newline at end of file diff --git a/DDFacet/Imager/ModModelMachine.py b/DDFacet/Imager/ModModelMachine.py index ae697a88..9885608d 100644 --- a/DDFacet/Imager/ModModelMachine.py +++ b/DDFacet/Imager/ModModelMachine.py @@ -40,6 +40,7 @@ def __init__(self,GD=None): self.SSDMM = None self.MSMFMM = None self.MORSANEMM = None + self.MUFFINMM = None self.HOGBOMMM = None def GiveInitialisedMMFromFile(self,FileName): @@ -114,12 +115,18 @@ def GiveMM(self,Mode=None): if self.MORSANEMM is None: print>> log, "Initialising MORESANE model machine" from DDFacet.Imager.MORESANE import ClassModelMachineMORESANE - self.MORESANEMM = ClassModelMachineMORESANE.ClassModelMachine( - self.GD, - GainMachine= ClassGainMachine.ClassGainMachine(GainMin=self.GD["MORESANE"]["loopgain"])) + self.MORESANEMM = ClassModelMachineMORESANE.ClassModelMachine(self.GD,GainMachine=ClassGainMachine.ClassGainMachine()) else: print>> log, "MORSANE model machine already initialised" return self.MORESANEMM + elif Mode == "MUFFIN": + if self.MUFFINMM is None: + print>> log, "Initialising MUFFIN model machine" + from DDFacet.Imager.MUFFIN import ClassModelMachineMUFFIN + self.MUFFINMM = ClassModelMachineMUFFIN.ClassModelMachine(self.GD,GainMachine=ClassGainMachine.ClassGainMachine()) + else: + print>> log, "MUFFIN model machine already initialised" + return self.MUFFINMM elif Mode == "Hogbom": if self.HOGBOMMM is None: print>> log, "Initialising HOGBOM model machine" diff --git a/DDFacet/Imager/SSD/ClassArrayMethodSSD.py b/DDFacet/Imager/SSD/ClassArrayMethodSSD.py index e9b4d161..3c81a79a 100644 --- a/DDFacet/Imager/SSD/ClassArrayMethodSSD.py +++ b/DDFacet/Imager/SSD/ClassArrayMethodSSD.py @@ -875,10 +875,10 @@ def runSingleMutation(self,DicoJob): Name="%sIsland_%5.5i_Individual_%4.4i"%(self.IdSharedMem,self.iIsland,iIndividual) individual=NpShared.GiveArray(Name) - Mut_pFlux, Mut_p0, Mut_pMove=DicoJob["mutConfig"] + Mut_pFlux, Mut_p0, Mut_pMove, Mut_pScale, Mut_pOffset=DicoJob["mutConfig"] individualOut,=self.MutMachine.mutGaussian(individual.copy(), - Mut_pFlux, Mut_p0, Mut_pMove) + Mut_pFlux, Mut_p0, Mut_pMove, Mut_pScale, Mut_pOffset) individual[:]=individualOut[:] self.result_queue.put({"Success": True, diff --git a/DDFacet/Imager/SSD/ClassImageDeconvMachineSSD.py b/DDFacet/Imager/SSD/ClassImageDeconvMachineSSD.py index 60a4635e..da2b2bff 100644 --- a/DDFacet/Imager/SSD/ClassImageDeconvMachineSSD.py +++ b/DDFacet/Imager/SSD/ClassImageDeconvMachineSSD.py @@ -35,7 +35,7 @@ from DDFacet.Other import MyPickle import multiprocessing import time -import ClassInitSSDModel +import ClassInitSSDModelHMP from DDFacet.Imager.SSD.GA.ClassEvolveGA import ClassEvolveGA from DDFacet.Imager.SSD.MCMC.ClassMetropolis import ClassMetropolis #try: # Genetic Algo @@ -227,8 +227,8 @@ def SearchIslands(self,Threshold): # ############################# print>>log," selected %i islands [out of %i] with peak flux > %.3g Jy"%(len(ListIslandsFiltered),len(ListIslands),Threshold) ListIslands=ListIslandsFiltered - #ListIslands=[np.load("errIsland_000000.keep.npy").tolist()] - + #ListIslands=[np.load("errIsland_000524.npy").tolist()] + ListIslands=IslandDistanceMachine.CalcCrossIslandFlux(ListIslands) ListIslands=IslandDistanceMachine.ConvexifyIsland(ListIslands) ListIslands=IslandDistanceMachine.MergeIslands(ListIslands) @@ -251,7 +251,7 @@ def SearchIslands(self,Threshold): def InitMSMF(self): self.DicoInitIndiv={} - if self.GD["SSDClean"]["MinSizeInitHMP"]==-1: return + if self.GD["GAClean"]["MinSizeInit"]==-1: return DoAbs=int(self.GD["Deconv"]["AllowNegative"]) print>>log, " Running minor cycle [MinorIter = %i/%i, SearchMaxAbs = %i]"%(self._niter,self.MaxMinorIter,DoAbs) @@ -286,30 +286,39 @@ def InitMSMF(self): self.ListSizeIslands.append(dd) - ListDoMSMFIslandsInit=[True if self.ListSizeIslands[iIsland]>=self.GD["SSDClean"]["MinSizeInitHMP"] else False for iIsland in range(len(self.ListIslands))] + ListDoIslandsInit=[True if self.ListSizeIslands[iIsland]>=self.GD["GAClean"]["MinSizeInit"] else False for iIsland in range(len(self.ListIslands))] #ListDoMSMFIslandsInit=[True if iIsland==16 else False for iIsland in range(len(self.ListIslands))] - print>>log," selected %i islands larger than %i pixels for HMP initialisation"%(np.count_nonzero(ListDoMSMFIslandsInit),self.GD["SSDClean"]["MinSizeInitHMP"]) - + print>>log," selected %i islands larger than %i pixels for initialisation"%(np.count_nonzero(ListDoIslandsInit),self.GD["GAClean"]["MinSizeInit"]) - - - - if np.count_nonzero(ListDoMSMFIslandsInit)>0: - InitMachine=ClassInitSSDModel.ClassInitSSDModelParallel(self.GD, - self.DicoVariablePSF, - self.DicoDirty, - self.ModelMachine.RefFreq, - self.GridFreqs, - self.DegridFreqs, - MainCache=self.maincache, - NCPU=self.NCPU, - IdSharedMem=self.IdSharedMem) + if np.count_nonzero(ListDoIslandsInit)>0: + if self.GD["GAClean"]["InitType"]=="HMP": + InitMachine=ClassInitSSDModelHMP.ClassInitSSDModelParallel(self.GD, + self.DicoVariablePSF, + self.DicoDirty, + self.ModelMachine.RefFreq, + self.GridFreqs, + self.DegridFreqs, + MainCache=self.maincache, + NCPU=self.NCPU, + IdSharedMem=self.IdSharedMem) + elif self.GD["GAClean"]["InitType"]=="MORESANE": + + import ClassInitSSDModelMoresane + InitMachine=ClassInitSSDModelMoresane.ClassInitSSDModelParallel(self.GD, + self.DicoVariablePSF, + self.DicoDirty, + self.ModelMachine.RefFreq, + self.GridFreqs, + self.DegridFreqs, + MainCache=self.maincache, + NCPU=self.NCPU, + IdSharedMem=self.IdSharedMem) InitMachine.setSSDModelImage(ModelImage) - self.DicoInitIndiv=InitMachine.giveDicoInitIndiv(self.ListIslands,ListDoIsland=ListDoMSMFIslandsInit) + self.DicoInitIndiv=InitMachine.giveDicoInitIndiv(self.ListIslands,ListDoIsland=ListDoIslandsInit) def setChannel(self,ch=0): @@ -458,12 +467,12 @@ def DeconvListIsland(self,ListIslands,ParallelMode="OverIsland",ListInitIslands= NCPU=np.min([NCPU,NIslands]) Parallel=True ParallelPerIsland=False + StopWhenQueueEmpty=True elif ParallelMode=="PerIsland": - NCPU=self.NCPU - Parallel=False + NCPU=1#self.NCPU + Parallel=True ParallelPerIsland=True - - StopWhenQueueEmpty=True + StopWhenQueueEmpty=True # ######### Debug # ParallelPerIsland=False diff --git a/DDFacet/Imager/SSD/ClassInitSSDModel.py b/DDFacet/Imager/SSD/ClassInitSSDModelHMP.py similarity index 93% rename from DDFacet/Imager/SSD/ClassInitSSDModel.py rename to DDFacet/Imager/SSD/ClassInitSSDModelHMP.py index 3b3aec17..4546c074 100644 --- a/DDFacet/Imager/SSD/ClassInitSSDModel.py +++ b/DDFacet/Imager/SSD/ClassInitSSDModelHMP.py @@ -146,21 +146,18 @@ def __init__(self,GD,DicoVariablePSF,DicoDirty,RefFreq,GridFreqs,DegridFreqs, self.GD=GD self.GD["Parallel"]["NCPU"]=1 #self.GD["HMP"]["Alpha"]=[0,0,1]#-1.,1.,5] - self.GD["HMP"]["Alpha"]=[-1.,1.,5] + self.GD["HMP"]["Alpha"]=[-3.,1.,5] self.GD["Deconv"]["Mode"]="HMP" self.GD["Deconv"]["CycleFactor"]=0 - self.GD["Deconv"]["PeakFactor"]=0.01 - self.GD["Deconv"]["RMSFactor"]=3. - self.GD["Deconv"]["Gain"]=.1 - self.GD["Deconv"]["AllowNegative"]=False + self.GD["Deconv"]["PeakFactor"]=0.0 + self.GD["Deconv"]["RMSFactor"]=self.GD["GAClean"]["RMSFactorInitHMP"] - self.GD["Deconv"]["MaxMinorIter"]=10000 + self.GD["Deconv"]["Gain"]=self.GD["GAClean"]["GainInitHMP"] + self.GD["Deconv"]["AllowNegative"]=self.GD["GAClean"]["AllowNegativeInitHMP"] + self.GD["Deconv"]["MaxMinorIter"]=self.GD["GAClean"]["MaxMinorIterInitHMP"] - if self.GD["SSDClean"]["ScalesInitHMP"] is not None: - self.GD["HMP"]["Scales"]=self.GD["SSDClean"]["ScalesInitHMP"] - else: - self.GD["HMP"]["Scales"]=[0,1,2,4,8,16,24,32] + self.GD["HMP"]["Scales"]=self.GD["GAClean"]["ScalesInitHMP"] self.GD["HMP"]["Ratios"]=[] #self.GD["MultiScale"]["Ratios"]=[] @@ -364,6 +361,7 @@ def giveModel(self,ListPixParms): #print "update" #time.sleep(30) self.DeconvMachine.Deconvolve(UpdateRMS=False) + #self.DeconvMachine.Plot() T.timeit("deconv %s"%str(self.DicoSubDirty["ImageCube"].shape)) #print "deconv" #time.sleep(30) @@ -391,6 +389,7 @@ def giveModel(self,ListPixParms): x,y=self.ArrayPixParms.T + # PSF,MeanPSF=self.DeconvMachine.PSFServer.GivePSF() # ConvModel=ClassConvMachineImages(PSF).giveConvModel(ModelImage*np.ones((self.NFreqBands,1,1,1))) # #T.timeit("Conv1") @@ -410,23 +409,43 @@ def giveModel(self,ListPixParms): # # pylab.draw() # # pylab.show(False) # # stop - + + # ModelOnes=np.zeros_like(ModelImage) + # ModelOnes[:,:,x,y]=1 + # ConvModelOnes=ClassConvMachineImages(PSF).giveConvModel(ModelOnes*np.ones((self.NFreqBands,1,1,1))) # SumConvModel=np.sum(ConvModel[:,:,x,y]) + # SumConvModelOnes=np.sum(ConvModelOnes[:,:,x,y]) # SumResid=np.sum(self.DeconvMachine._CubeDirty[:,:,x,y]) # SumConvModel=np.max([SumConvModel,1e-6]) - # factor=(SumResid+SumConvModel)/SumConvModel - # fMult=1. - # if 1.1: + AModel=self.ModelMachine.GiveSpectralIndexMap()[0,0,x,y] + else: + AModel=np.zeros_like(SModel) T.timeit("spec index") return SModel,AModel diff --git a/DDFacet/Imager/SSD/ClassInitSSDModelMoresane.py b/DDFacet/Imager/SSD/ClassInitSSDModelMoresane.py new file mode 100644 index 00000000..edcf313f --- /dev/null +++ b/DDFacet/Imager/SSD/ClassInitSSDModelMoresane.py @@ -0,0 +1,555 @@ +import numpy as np +from DDFacet.Imager.MORESANE import ClassImageDeconvMachineMoresane +import copy +from DDFacet.ToolsDir.GiveEdges import GiveEdges +from DDFacet.ToolsDir.GiveEdges import GiveEdgesDissymetric +from DDFacet.Imager.ClassPSFServer import ClassPSFServer +from DDFacet.Imager.ModModelMachine import ClassModModelMachine +import multiprocessing +from DDFacet.Other import ClassTimeIt +from DDFacet.Other.progressbar import ProgressBar +import time +from DDFacet.Array import NpShared +from DDFacet.Other import MyLogger +log=MyLogger.getLogger("ClassInitSSDModel") +import traceback +from DDFacet.Other import ModColor +from ClassConvMachine import ClassConvMachineImages +from DDFacet.Imager import ClassMaskMachine + +SilentModules=["ClassPSFServer","ClassImageDeconvMachine","GiveModelMachine","ClassModelMachineMoresane","ClassModelMachineSSD","pymoresane.main"] + +class ClassInitSSDModelParallel(): + def __init__(self,GD,DicoVariablePSF,DicoDirty,RefFreq,GridFreqs,DegridFreqs,MainCache=None,NCPU=1,IdSharedMem=""): + self.DicoVariablePSF=DicoVariablePSF + self.DicoDirty=DicoDirty + GD=copy.deepcopy(GD) + self.RefFreq=RefFreq + self.GridFreqs=GridFreqs + self.DegridFreqs=DegridFreqs + self.MainCache=MainCache + self.GD=GD + self.NCPU=NCPU + self.IdSharedMem=IdSharedMem + print>>log,"Initialise MORESANE machine" + self.InitMachine=ClassInitSSDModel(self.GD, + self.DicoVariablePSF, + self.DicoDirty, + self.RefFreq, + self.GridFreqs, + self.DegridFreqs, + MainCache=self.MainCache, + IdSharedMem=self.IdSharedMem) + + def setSSDModelImage(self,ModelImage): + self.ModelImage=ModelImage + + def giveDicoInitIndiv(self,ListIslands,ListDoIsland=None,Parallel=True): + NCPU=self.NCPU + work_queue = multiprocessing.JoinableQueue() + ListIslands=ListIslands#[300:308] + DoIsland=True + + + + for iIsland in range(len(ListIslands)): + if ListDoIsland is not None: + DoIsland=ListDoIsland[iIsland] + if DoIsland: work_queue.put({"iIsland":iIsland}) + + result_queue=multiprocessing.JoinableQueue() + NJobs=work_queue.qsize() + workerlist=[] + + MyLogger.setSilent(SilentModules) + #MyLogger.setLoud(SilentModules) + + #MyLogger.setLoud("ClassImageDeconvMachineMSMF") + + DicoHMPFunctions=None + + print>>log,"Launch MORESANE workers" + for ii in range(NCPU): + W = WorkerInitMSMF(work_queue, + result_queue, + self.GD, + self.DicoVariablePSF, + self.DicoDirty, + self.RefFreq, + self.GridFreqs, + self.DegridFreqs, + self.MainCache, + self.ModelImage, + ListIslands, + self.IdSharedMem, + DicoHMPFunctions) + workerlist.append(W) + if Parallel: + workerlist[ii].start() + + timer = ClassTimeIt.ClassTimeIt() + pBAR = ProgressBar(Title=" MORESANing islands ") + #pBAR.disable() + pBAR.render(0, NJobs) + iResult = 0 + if not Parallel: + for ii in range(NCPU): + workerlist[ii].run() # just run until all work is completed + + self.DicoInitIndiv={} + while iResult < NJobs: + DicoResult = None + if result_queue.qsize() != 0: + try: + DicoResult = result_queue.get() + except: + pass + + if DicoResult == None: + time.sleep(0.5) + continue + + if DicoResult["Success"]: + iResult+=1 + NDone=iResult + + pBAR.render(NDone,NJobs) + + iIsland=DicoResult["iIsland"] + NameDico="%sDicoInitIsland_%5.5i"%(self.IdSharedMem,iIsland) + Dico=NpShared.SharedToDico(NameDico) + self.DicoInitIndiv[iIsland]=copy.deepcopy(Dico) + NpShared.DelAll(NameDico) + + + + if Parallel: + for ii in range(NCPU): + workerlist[ii].shutdown() + workerlist[ii].terminate() + workerlist[ii].join() + + #MyLogger.setLoud(["pymoresane.main"]) + #MyLogger.setLoud(["ClassImageDeconvMachineMSMF","ClassPSFServer","ClassMultiScaleMachine","GiveModelMachine","ClassModelMachineMSMF"]) + return self.DicoInitIndiv + +###################################################################################################### + +class ClassInitSSDModel(): + def __init__(self,GD,DicoVariablePSF,DicoDirty,RefFreq,GridFreqs,DegridFreqs, + MainCache=None, + IdSharedMem="", + DoWait=False, + DicoHMPFunctions=None): + self.DicoVariablePSF=DicoVariablePSF + self.DicoDirty=DicoDirty + GD=copy.deepcopy(GD) + self.RefFreq=RefFreq + self.GridFreqs=GridFreqs + self.DegridFreqs=DegridFreqs + self.GD=GD + self.GD["Parallel"]["NCPU"]=1 + #self.GD["HMP"]["Alpha"]=[0,0,1]#-1.,1.,5] + self.GD["HMP"]["Alpha"]=[-1.,1.,5] + self.GD["Deconv"]["Mode"]="MORESANE" + + self.GD["Deconv"]["CycleFactor"]=0 + self.GD["Deconv"]["PeakFactor"]=0.0 + self.GD["Deconv"]["RMSFactor"]=self.GD["GAClean"]["RMSFactorInitHMP"] + + self.GD["Deconv"]["Gain"]=self.GD["GAClean"]["GainInitHMP"] + self.GD["Deconv"]["AllowNegative"]=self.GD["GAClean"]["AllowNegativeInitHMP"] + self.GD["Deconv"]["MaxMinorIter"]=self.GD["GAClean"]["MaxMinorIterInitHMP"] + + MyLogger.setSilent(SilentModules) + + self.GD["HMP"]["Scales"]=self.GD["GAClean"]["ScalesInitHMP"] + + self.GD["HMP"]["Ratios"]=[] + #self.GD["MultiScale"]["Ratios"]=[] + self.GD["HMP"]["NTheta"]=4 + + self.GD["HMP"]["SolverMode"]="NNLS" + #self.GD["MultiScale"]["SolverMode"]="PI" + + self.NFreqBands=len(DicoVariablePSF["freqs"]) + MinorCycleConfig=dict(self.GD["Deconv"]) + MinorCycleConfig["NCPU"]=self.GD["Parallel"]["NCPU"] + MinorCycleConfig["NFreqBands"]=self.NFreqBands + MinorCycleConfig["GD"] = self.GD + MinorCycleConfig["GridFreqs"] = self.GridFreqs + MinorCycleConfig["DegridFreqs"] = self.DegridFreqs + + #MinorCycleConfig["RefFreq"] = self.RefFreq + + ModConstructor = ClassModModelMachine(self.GD) + ModelMachine = ModConstructor.GiveMM(Mode="MORESANE") + ModelMachine.setRefFreq(self.RefFreq) + MinorCycleConfig["ModelMachine"]=ModelMachine + #MinorCycleConfig["CleanMaskImage"]=None + self.MinorCycleConfig=MinorCycleConfig + self.DeconvMachine=ClassImageDeconvMachineMoresane.ClassImageDeconvMachine(MainCache=MainCache, + CacheSharedMode=True, + ParallelMode=False, + RefFreq=self.RefFreq, + CacheFileName="HMP_Init", + IdSharedMem=IdSharedMem, + **self.MinorCycleConfig) + self.GD["Mask"]["Auto"]=False + self.GD["Mask"]["External"]=None + self.MaskMachine=ClassMaskMachine.ClassMaskMachine(self.GD) + self.DeconvMachine.setMaskMachine(self.MaskMachine) + + self.DicoHMPFunctions=None + if self.DicoHMPFunctions is not None: + self.DeconvMachine.set_DicoHMPFunctions(self.DicoHMPFunctions) + + self.Margin=50 + self.DicoDirty=DicoDirty + self.Dirty=DicoDirty["ImageCube"] + self.MeanDirty=DicoDirty["MeanImage"] + + #print "Start 3" + self.DeconvMachine.Init(PSFVar=self.DicoVariablePSF,PSFAve=self.DicoVariablePSF["PSFSideLobes"], + GridFreqs=self.GridFreqs,DegridFreqs=self.DegridFreqs,DoWait=DoWait,RefFreq=self.RefFreq) + + if DoWait: + print "IINit3" + time.sleep(10) + print "Start 4" + + self.DeconvMachine.Update(self.DicoDirty,DoSetMask=False) + if DoWait: + print "IINit4" + time.sleep(10) + + #self.DeconvMachine.updateRMS() + + #self.DicoBasicModelMachine=copy.deepcopy(self.DeconvMachine.ModelMachine.DicoSMStacked) + + def setSubDirty(self,ListPixParms): + T=ClassTimeIt.ClassTimeIt("InitSSD.setSubDirty") + T.disable() + + x,y=np.array(ListPixParms).T + x0,x1=x.min(),x.max()+1 + y0,y1=y.min(),y.max()+1 + dx=x1-x0+self.Margin + dy=y1-y0+self.Margin + Size=np.max([dx,dy]) + if Size%2==0: Size+=1 + _,_,N0,_=self.Dirty.shape + + xc0,yc0=int((x1+x0)/2.),int((y1+y0)/2.) + self.xy0=xc0,yc0 + self.DeconvMachine.PSFServer.setLocation(*self.xy0) + + N1=Size + xc1=yc1=N1/2 + Aedge,Bedge=GiveEdges((xc0,yc0),N0,(xc1,yc1),N1) + x0d,x1d,y0d,y1d=Aedge + x0p,x1p,y0p,y1p=Bedge + self.SubDirty=self.Dirty[:,:,x0d:x1d,y0d:y1d].copy() + T.timeit("0") + self.blc=(x0d,y0d) + self.DeconvMachine.PSFServer.setBLC(self.blc) + _,_,nx,ny=self.SubDirty.shape + ArrayPixParms=np.array(ListPixParms) + ArrayPixParms[:,0]-=x0d + ArrayPixParms[:,1]-=y0d + self.ArrayPixParms=ArrayPixParms + self.DicoSubDirty={} + for key in self.DicoDirty.keys(): + if key in ["ImageCube", "MeanImage",'FacetNorm',"JonesNorm"]: + self.DicoSubDirty[key]=self.DicoDirty[key][...,x0d:x1d,y0d:y1d].copy() + else: + self.DicoSubDirty[key]=self.DicoDirty[key] + + T.timeit("1") + # ModelImage=np.zeros_like(self.Dirty) + # ModelImage[:,:,N0/2,N0/2]=10 + # ModelImage[:,:,N0/2+3,N0/2]=10 + # ModelImage[:,:,N0/2-2,N0/2-1]=10 + # self.setSSDModelImage(ModelImage) + + # Mask=np.zeros((nx,ny),np.bool8) + # Mask[x,y]=1 + # self.SubMask=Mask + + + x,y=ArrayPixParms.T + Mask=np.zeros(self.DicoSubDirty["ImageCube"].shape[-2::],np.bool8) + Mask[x,y]=1 + self.SubMask=Mask + + + if self.SSDModelImage is not None: + self.SubSSDModelImage=self.SSDModelImage[:,:,x0d:x1d,y0d:y1d].copy() + for ch in range(self.NFreqBands): + self.SubSSDModelImage[ch,0][np.logical_not(self.SubMask)]=0 + self.addSubModelToSubDirty() + T.timeit("2") + + + def setSSDModelImage(self,ModelImage): + self.SSDModelImage=ModelImage + + def giveConvModel(self,SubModelImage): + + PSF,MeanPSF=self.DeconvMachine.PSFServer.GivePSF() + ConvModel=ClassConvMachineImages(PSF).giveConvModel(SubModelImage) + + # ConvModel=np.zeros_like(SubModelImage) + # nch,_,N0x,N0y=ConvModel.shape + # indx,indy=np.where(SubModelImage[0,0]!=0) + # xc,yc=N0x/2,N0y/2 + # N1=PSF.shape[-1] + # #T.timeit("0") + # for i,j in zip(indx.tolist(),indy.tolist()): + # ThisPSF=np.roll(np.roll(PSF,i-xc,axis=-2),j-yc,axis=-1) + # Aedge,Bedge=GiveEdgesDissymetric((xc,yc),(N0x,N0y),(N1/2,N1/2),(N1,N1)) + # x0d,x1d,y0d,y1d=Aedge + # x0p,x1p,y0p,y1p=Bedge + # ConvModel[...,x0d:x1d,y0d:y1d]+=ThisPSF[...,x0p:x1p,y0p:y1p]*SubModelImage[...,i,j].reshape((-1,1,1,1)) + # #T.timeit("1 %s"%(str(ConvModel.shape))) + + return ConvModel + + + + def addSubModelToSubDirty(self): + T=ClassTimeIt.ClassTimeIt("InitSSD.addSubModelToSubDirty") + T.disable() + ConvModel=self.giveConvModel(self.SubSSDModelImage) + _,_,N0x,N0y=ConvModel.shape + MeanConvModel=np.mean(ConvModel,axis=0).reshape((1,1,N0x,N0y)) + self.DicoSubDirty["ImageCube"]+=ConvModel + self.DicoSubDirty['MeanImage']+=MeanConvModel + #print "MAX=",np.max(self.DicoSubDirty['MeanImage']) + T.timeit("2") + + # import pylab + # pylab.clf() + # ax=pylab.subplot(1,3,1) + # pylab.imshow(self.SubSSDModelImage[0,0],interpolation="nearest") + # pylab.subplot(1,3,2,sharex=ax,sharey=ax) + # pylab.imshow(PSF[0,0],interpolation="nearest") + # pylab.subplot(1,3,3,sharex=ax,sharey=ax) + # pylab.imshow(ConvModel[0,0],interpolation="nearest") + # pylab.draw() + # pylab.show(False) + # pylab.pause(0.1) + + + def giveModel(self,ListPixParms): + T=ClassTimeIt.ClassTimeIt("giveModel") + T.disable() + self.setSubDirty(ListPixParms) + T.timeit("setsub") + ModConstructor = ClassModModelMachine(self.GD) + ModelMachine = ModConstructor.GiveMM(Mode=self.GD["Deconv"]["Mode"]) + #print "ModelMachine" + #time.sleep(30) + T.timeit("giveMM") + self.ModelMachine=ModelMachine + #self.ModelMachine.DicoSMStacked=self.DicoBasicModelMachine + self.ModelMachine.setRefFreq(self.RefFreq,Force=True) + self.ModelMachine.setFreqMachine(self.GridFreqs,self.DegridFreqs) + self.MinorCycleConfig["ModelMachine"] = ModelMachine + self.ModelMachine.setModelShape(self.SubDirty.shape) + #self.ModelMachine.setListComponants(self.DeconvMachine.ModelMachine.ListScales) + T.timeit("setlistcomp") + + self.DeconvMachine.Update(self.DicoSubDirty,DoSetMask=False) + #self.DeconvMachine.updateMask(np.logical_not(self.SubMask)) + self.DeconvMachine.updateModelMachine(ModelMachine) + #self.DeconvMachine.resetCounter() + T.timeit("update") + #print "update" + #time.sleep(30) + self.DeconvMachine.Deconvolve() + T.timeit("deconv %s"%str(self.DicoSubDirty["ImageCube"].shape)) + #print "deconv" + #time.sleep(30) + + ModelImage=self.ModelMachine.GiveModelImage() + T.timeit("getmodel") + + # import pylab + # pylab.clf() + # pylab.subplot(2,2,1) + # pylab.imshow(self.DicoDirty["MeanImage"][0,0,:,:],interpolation="nearest") + # pylab.colorbar() + # pylab.subplot(2,2,2) + # pylab.imshow(self.DicoSubDirty["MeanImage"][0,0,:,:],interpolation="nearest") + # pylab.colorbar() + # pylab.subplot(2,2,3) + # pylab.imshow(self.SubMask,interpolation="nearest") + # pylab.colorbar() + # pylab.subplot(2,2,4) + # pylab.imshow(ModelImage[0,0],interpolation="nearest") + # pylab.colorbar() + # pylab.draw() + # pylab.show(False) + # pylab.pause(0.1) + + + x,y=self.ArrayPixParms.T + + # PSF,MeanPSF=self.DeconvMachine.PSFServer.GivePSF() + # ConvModel=ClassConvMachineImages(PSF).giveConvModel(ModelImage*np.ones((self.NFreqBands,1,1,1))) + # #T.timeit("Conv1") + # #print "done1" + # #ConvModel=self.giveConvModel(ModelImage*np.ones((self.NFreqBands,1,1,1))) + # # print "done2" + # # T.timeit("Conv2") + # # import pylab + # # pylab.clf() + # # pylab.subplot(1,3,1) + # # pylab.imshow(ConvModel[0,0],interpolation="nearest") + # # pylab.subplot(1,3,2) + # # pylab.imshow(ConvModel1[0,0],interpolation="nearest") + # # pylab.subplot(1,3,3) + # # pylab.imshow((ConvModel-ConvModel1)[0,0],interpolation="nearest") + # # pylab.colorbar() + # # pylab.draw() + # # pylab.show(False) + # # stop + + # ModelOnes=np.zeros_like(ModelImage) + # ModelOnes[:,:,x,y]=1 + # ConvModelOnes=ClassConvMachineImages(PSF).giveConvModel(ModelOnes*np.ones((self.NFreqBands,1,1,1))) + + # SumConvModel=np.sum(ConvModel[:,:,x,y]) + # SumConvModelOnes=np.sum(ConvModelOnes[:,:,x,y]) + # SumResid=np.sum(self.DeconvMachine._CubeDirty[:,:,x,y]) + + # SumConvModel=np.max([SumConvModel,1e-6]) + + # factor=(SumResid+SumConvModel)/SumConvModel + + + # ############### + #fMult=1. + #if 1.>log," Convexify islands" ListConvexIslands=[] for Island in ListIslands: + points=np.array(Island) + x,y=points.T + Cx=(np.abs(x.min()-x.max())==0) + Cy=(np.abs(y.min()-y.max())==0) + if (x.size<=3) or Cx or Cy: + ListConvexIslands.append(Island) + continue try: - points=np.array(Island) hull = ConvexHull(points) Contour = np.array( [hull.points[hull.vertices, 0], hull.points[hull.vertices, 1]]) poly2 = Contour.T - x,y=points.T x0,x1=x.min(),x.max() y0,y1=y.min(),y.max() diff --git a/DDFacet/Imager/SSD/ClassModelMachineSSD.py b/DDFacet/Imager/SSD/ClassModelMachineSSD.py index bbd0ad30..ce8a6ffc 100644 --- a/DDFacet/Imager/SSD/ClassModelMachineSSD.py +++ b/DDFacet/Imager/SSD/ClassModelMachineSSD.py @@ -326,44 +326,59 @@ def GiveModelImage(self,FreqIn=None,out=None): def setListComponants(self,ListScales): self.ListScales=ListScales - def GiveSpectralIndexMap(self, threshold=0.1, save_dict=True): - # Get the model image - IM = self.GiveModelImage(self.FreqMachine.Freqsp) - nchan, npol, Nx, Ny = IM.shape - - # Fit the alpha map - self.FreqMachine.FitAlphaMap(IM[:, 0, :, :], threshold=threshold) # should set threshold based on SNR of final residual - - if save_dict: - FileName = self.GD['Output']['Name'] + ".Dicoalpha" - print>>log, "Saving componentwise SPI map to %s"%FileName - - MyPickle.Save(self.FreqMachine.alpha_dict, FileName) - - return self.FreqMachine.weighted_alpha_map.reshape((1, 1, Nx, Ny)) - - - # def GiveSpectralIndexMap(self,CellSizeRad=1.,GaussPars=[(1,1,0)],DoConv=True): - # - # - # dFreq=1e6 - # RefFreq=self.DicoSMStacked["RefFreq"] - # f0=RefFreq/1.5#self.DicoSMStacked["AllFreqs"].min() - # f1=RefFreq*1.5#self.DicoSMStacked["AllFreqs"].max() - # M0=self.GiveModelImage(f0) - # M1=self.GiveModelImage(f1) - # if DoConv: - # M0=ModFFTW.ConvolveGaussian(M0,CellSizeRad=CellSizeRad,GaussPars=GaussPars) - # M1=ModFFTW.ConvolveGaussian(M1,CellSizeRad=CellSizeRad,GaussPars=GaussPars) - # - # Np=1000 - # indx,indy=np.int64(np.random.rand(Np)*M0.shape[0]),np.int64(np.random.rand(Np)*M0.shape[1]) - # med=np.median(np.abs(M0[:,:,indx,indy])) - # - # Mask=((M1>100*med)&(M0>100*med)) - # alpha=np.zeros_like(M0) - # alpha[Mask]=(np.log(M0[Mask])-np.log(M1[Mask]))/(np.log(f0/f1)) - # return alpha + # def GiveSpectralIndexMap(self, threshold=0.1, save_dict=True): + # # Get the model image + # IM = self.GiveModelImage(self.FreqMachine.Freqsp) + # nchan, npol, Nx, Ny = IM.shape + # # Fit the alpha map + # self.FreqMachine.FitAlphaMap(IM[:, 0, :, :], threshold=threshold) # should set threshold based on SNR of final residual + # if save_dict: + # FileName = self.GD['Output']['Name'] + ".Dicoalpha" + # print>>log, "Saving componentwise SPI map to %s"%FileName + # MyPickle.Save(self.FreqMachine.alpha_dict, FileName) + # return self.FreqMachine.weighted_alpha_map.reshape((1, 1, Nx, Ny)) + + + def GiveSpectralIndexMap(self,CellSizeRad=1.,GaussPars=[(1,1,0)],DoConv=True,MaxSpi=100,MaxDR=1e+6,threshold=None): + + dFreq=1e6 + RefFreq=self.DicoSMStacked["RefFreq"] + f0=RefFreq/1.5#self.DicoSMStacked["AllFreqs"].min() + f1=RefFreq*1.5#self.DicoSMStacked["AllFreqs"].max() + M0=self.GiveModelImage(f0) + M1=self.GiveModelImage(f1) + if DoConv: + #M0=ModFFTW.ConvolveGaussian(M0,CellSizeRad=CellSizeRad,GaussPars=GaussPars) + #M1=ModFFTW.ConvolveGaussian(M1,CellSizeRad=CellSizeRad,GaussPars=GaussPars) + #M0,_=ModFFTW.ConvolveGaussianWrapper(M0,Sig=GaussPars[0][0]/CellSizeRad) + #M1,_=ModFFTW.ConvolveGaussianWrapper(M1,Sig=GaussPars[0][0]/CellSizeRad) + M0,_=ModFFTW.ConvolveGaussianScipy(M0,Sig=GaussPars[0][0]/CellSizeRad) + M1,_=ModFFTW.ConvolveGaussianScipy(M1,Sig=GaussPars[0][0]/CellSizeRad) + + + # compute threshold for alpha computation by rounding DR threshold to .1 digits (i.e. 1.65e-6 rounds to 1.7e-6) + if threshold is not None: + minmod = threshold + elif not np.all(M0==0): + minmod = float("%.1e"%(np.max(np.abs(M0))/MaxDR)) + else: + minmod=1e-6 + + # mask out pixels above threshold + mask=(M1>log,"computing alpha map for model pixels above %.1e Jy (based on max DR setting of %g)"%(minmod,MaxDR) + M0[mask]=minmod + M1[mask]=minmod + alpha = (np.log(M0)-np.log(M1))/(np.log(f0/f1)) + alpha[mask] = 0 + + # Np=1000 + # indx,indy=np.int64(np.random.rand(Np)*M0.shape[0]),np.int64(np.random.rand(Np)*M0.shape[1]) + # med=np.median(np.abs(M0[:,:,indx,indy])) + # Mask=((M1>100*med)&(M0>100*med)) + # alpha=np.zeros_like(M0) + # alpha[Mask]=(np.log(M0[Mask])-np.log(M1[Mask]))/(np.log(f0/f1)) + return alpha def RemoveNegComponants(self): @@ -381,6 +396,8 @@ def RemoveNegComponants(self): def FilterNegComponants(self,box=20,sig=3,RemoveNeg=True): print>>log, "Cleaning model dictionary from negative components with (box, sig) = (%i, %i)"%(box,sig) + + print>>log," Number of componants before filtering: %i"%len(self.DicoSMStacked["Comp"]) ModelImage=self.GiveModelImage(self.DicoSMStacked["RefFreq"])[0,0] Min=scipy.ndimage.filters.minimum_filter(ModelImage,(box,box)) @@ -399,6 +416,9 @@ def FilterNegComponants(self,box=20,sig=3,RemoveNeg=True): del(self.DicoSMStacked["Comp"][key]) except: print>>log, " Component at (%i, %i) not in dict "%key + print>>log," Number of componants after filtering: %i"%len(self.DicoSMStacked["Comp"]) + + def CleanMaskedComponants(self,MaskName): print>>log, "Cleaning model dictionary from masked components using %s"%(MaskName) diff --git a/DDFacet/Imager/SSD/ClassMutate.py b/DDFacet/Imager/SSD/ClassMutate.py index b0326569..0b23ba2f 100644 --- a/DDFacet/Imager/SSD/ClassMutate.py +++ b/DDFacet/Imager/SSD/ClassMutate.py @@ -94,7 +94,7 @@ def setData(self,DicoData): - def mutGaussian(self,individual, pFlux, p0, pMove,FactorAccelerate=1.): + def mutGaussian(self,individual, pFlux, p0, pMove, pScale, pOffset,FactorAccelerate=1.): #return individual, T= ClassTimeIt.ClassTimeIt() T.disable() @@ -124,7 +124,7 @@ def mutGaussian(self,individual, pFlux, p0, pMove,FactorAccelerate=1.): T.timeit("start2") - PMat=np.array([0.,pFlux, p0, pMove]) + PMat=np.array([0.,pFlux, p0, pMove, pScale, pOffset]) PMat/=np.sum(PMat) PMat=np.cumsum(PMat) @@ -140,26 +140,30 @@ def mutGaussian(self,individual, pFlux, p0, pMove,FactorAccelerate=1.): if Type==0: NMax=int(np.max([3.,10])) N=int(random.uniform(1, NMax)) + indR=sorted(list(set(np.int32(np.random.rand(N)*NNonZero).tolist()))) + indSel=ind[indR] # zero a pixel elif Type==1: N=np.max([(NNonZero/10),1]) + indR=sorted(list(set(np.int32(np.random.rand(N)*NNonZero).tolist()))) + indSel=ind[indR] # move a pixel - else: + elif Type==2: NMax=int(np.max([3.,self.PM.NPixListParms/10])) NMax=np.min([NMax,10]) N=int(random.uniform(1, NMax)) - #N=int(random.uniform(1, individual.shape[1]/100)) - - # InReg=random.uniform(-1,1) - # if InReg<0: - # InReg=-1 - - + indR=sorted(list(set(np.int32(np.random.rand(N)*NNonZero).tolist()))) + indSel=ind[indR] + elif Type==3: + FactorScale=1.+np.random.randn(1)[0]*0.01 + indSel=ind#np.arange(Af.size) + elif Type==4: + SMin=np.min(np.abs(Af[Af!=0.])) + Offset=np.random.randn(1)[0]*SMin + indSel=ind#np.arange(Af.size) #print pFlux, p0, pMove #print "PPPPPP",PMat,RType,Type - indR=sorted(list(set(np.int32(np.random.rand(N)*NNonZero).tolist()))) - indSel=ind[indR] #Type=0 #print "Type:",Type,RType @@ -212,6 +216,11 @@ def mutGaussian(self,individual, pFlux, p0, pMove,FactorAccelerate=1.): # for iReg in [1,3,5,7]: # individual=self.MovePix(individual,iPix,Flux,InReg=iReg) + if Type==3: + Af[iPix]*=FactorScale + if Type==4: + Af[iPix]+=Offset + if "GSig" in self.PM.SolveParam: GSig=self.PM.ArrayToSubArray(individual,"GSig") GSig[GSig<0]=0 diff --git a/DDFacet/Imager/SSD/GA/ClassEvolveGA.py b/DDFacet/Imager/SSD/GA/ClassEvolveGA.py index ca7bc533..16eac685 100644 --- a/DDFacet/Imager/SSD/GA/ClassEvolveGA.py +++ b/DDFacet/Imager/SSD/GA/ClassEvolveGA.py @@ -88,8 +88,8 @@ def InitEvolutionAlgo(self): # toolbox.register("mate", tools.cxOrdered) # toolbox.register("mutate", tools.mutGaussian, indpb=0.3, mu=0.0, sigma=.1) #toolbox.register("mutate", self.ArrayMethodsMachine.mutGaussian, pFlux=0.3, p0=0.3, pMove=0.3) - self.MutConfig=pFlux,p0,pMove=0.5,0.5,0.5 - toolbox.register("mutate", self.ArrayMethodsMachine.mutGaussian, pFlux=0.2, p0=0.5, pMove=0.2) + self.MutConfig=pFlux,p0,pMove,pScale,pOffset=0.5,0.5,0.5,0.5,0.5 + toolbox.register("mutate", self.ArrayMethodsMachine.mutGaussian, pFlux=0.2, p0=0.5, pMove=0.2, pScale=0.2, pOffset=0.2) toolbox.register("select", tools.selTournament, tournsize=3) #toolbox.register("select", Select.selTolTournament, tournsize=3, Tol=4) diff --git a/DDFacet/Parset/DefaultParset.cfg b/DDFacet/Parset/DefaultParset.cfg index 02cc5518..d6e0b127 100755 --- a/DDFacet/Parset/DefaultParset.cfg +++ b/DDFacet/Parset/DefaultParset.cfg @@ -136,7 +136,7 @@ Affinity = 1 # pin processes to cores. Alternatively specifies a list of length NCPU. Alternatively "disable" to disable affinity settings Alternatively "enable_ht" uses stepping of 1 (equivalent to Parallel.Affinity=1), will use all vthreads - the obvious exception is if HT is disabled at BIOS level - Alternatively "disable_ht" autodetects the NUMA layout of the chip for Debian-based systems and don't use both vthreads per core + Alternatively "disable_ht" autodetects the NUMA layout of the chip for Debian-based systems and dont use both vthreads per core Use 1 if unsure. MainProcessAffinity = 0 # this should be set to a core that is not used by forked processes, this option is ignored when using option "disable or disable_ht" for Parallel.Affinity @@ -150,8 +150,9 @@ VisData = auto # Cache visibility data and flags at run LastResidual = 1 # Cache last residual data (at end of last minor cycle) #type:bool Dir = # Directory to store caches in. Default is to keep cache next to the MS, but this can cause performance issues with e.g. NFS volumes. If you have fast local storage, point to it. %metavar:DIR -DirWisdomFFTW = ~/.fftw_wisdom # Directory in which to store the FFTW wisdom files -ResetWisdom = 0 # Reset Wisdom file #type:bool +DirWisdomFFTW = ~/.fftw_wisdom # Directory in which to store the FFTW wisdom files +ResetWisdom = 0 # Reset Wisdom file #type:bool +CacheCF = True [Beam] _Help = Apply E-Jones (beam) during imaging @@ -183,24 +184,24 @@ NDegridBand = 0 # Number of image bands for degridding. 0 means degrid each [DDESolutions] _Help = Apply DDE solutions during imaging (@cyriltasse please document this section) -DDSols = -GlobalNorm = None # MeanAbs -JonesNormList = AP -JonesMode = Full # #options:Scalar|Diag|Full -DDModeGrid = AP -DDModeDeGrid = AP -ScaleAmpGrid = 0 -ScaleAmpDeGrid = 0 -CalibErr = 10. -Type = Nearest # #options:Krigging|Nearest -Scale = 1. # #metavar:DEG -gamma = 4. -RestoreSub = False -ReWeightSNR = 0. +DDSols = # Name of the DDE solution file +GlobalNorm = None # Option to normalise the Jones matrices (options: MeanAbs, MeanAbsAnt, BLBased or SumBLBased). See code for more detail +JonesNormList = AP # Deprecated? +JonesMode = Full # #options:Scalar|Diag|Full +DDModeGrid = AP # In the gridding step, apply Jones matrices Amplitude (A) or Phase (P) or Amplitude&Phase (AP) +DDModeDeGrid = AP # In the degridding step, apply Jones matrices Amplitude (A) or Phase (P) or Amplitude&Phase (AP) +ScaleAmpGrid = 0 # Deprecated? +ScaleAmpDeGrid = 0 # Deprecated? +CalibErr = 10. # Deprecated? +Type = Nearest # Deprecated? #options:Krigging|Nearest +Scale = 1. # Deprecated? #metavar:DEG +gamma = 4. # Deprecated? +RestoreSub = False # Deprecated? +ReWeightSNR = 0. # Deprecated? [Deconv] _Help = Common deconvolution options. Not all of these apply to all deconvolution modes -Mode = HMP # Deconvolution algorithm. #options:HMP|Hogbom|SSD|GAClean +Mode = HMP # Deconvolution algorithm. #options:HMP|Hogbom|SSD MaxMajorIter = 20 # Max number of major cycles. #type:int #metavar:N MaxMinorIter = 20000 # Max number of (overall) minor cycle iterations (HMP, Hogbom). #type:int #metavar:N AllowNegative = 1 # Allow negative components (HMP, Hogbom). #type:bool @@ -249,6 +250,7 @@ Support = 0 # Basis function support size. If 0, determined PeakWeightImage = None # weigh the peak finding by given image Kappa = 0 # Regularization parameter. If stddev of per-alpha solutions exceeds the maximum solution amplitude divided by Kappa, forces a fully-regularized solution. Use 0 for no such regularization. #type:float +OuterSpaceTh = 2. [Hogbom] PolyFitOrder = 3 # polynomial order for frequency fitting @@ -273,42 +275,33 @@ ConvFFTSwitch = 1000 NEnlargePars = 0 NEnlargeData = 2 RestoreMetroSwitch = 0 -MinSizeInitHMP = -1 -ScalesInitHMP = None MinMaxGroupDistance = [10,50] [GAClean] -NSourceKin = 100 -NMaxGen = 100 +NSourceKin = 50 +NMaxGen = 50 +MinSizeInit = 10 +InitType = HMP +ScalesInitHMP = [0,1,2,4,8,16,24,32] +GainInitHMP = 0.1 +MaxMinorIterInitHMP = 10000 +AllowNegativeInitHMP = False +RMSFactorInitHMP = 3. [MORESANE] -_Help = MORESANE deconvolution mode settings. @JulienNGirard please document, and maybe rename using InitialCaps convention. -MOSolvePars = [S] # (default=[S]): Solving parameters for MORESANE (only Flux & monochannel for the moment) -Accuracy = 1e-6 # (default=1e-6): Threshold on the standard deviation of the residual noise. Exit main loop when this threshold is reached. -AllOnGpu = False # (default=False): Boolean specifier to toggle all gpu modes on. -ConvDevice = cpu # (default='cpu'): Specifier for device to be used - cpu or gpu. -ConvMode = linear # (default='linear'): Specifier for convolution mode - linear or circular. -CoreCount = 1 # (default=1): For multiprocessing, specifies the number of cores. -DecomMode = ser # (default='ser'): (ser) = serial, (mp) = multiprocessing, (gpu) implementation of the IUWT decomposition -EdgeExcl = 0 # (default=0) : Width (in pixel) of the area to discard from border -EdgeOffset = 0 # (default=0): Numeric value for an additional user-specified number of edge pixels to be ignored. This is added to the minimum suppression. -EdgeSuppression = False # (default=False): Boolean specifier for whether or not the edges are to be suppressed. -EnforcePositivity = False # (default=False): Boolean specifier for whether or not a model must be strictly positive. -ExtractionMode = cpu # (default='cpu'): Specifier for mode to be used - cpu or gpu. -FluxThreshold = 0 # (default=0): Float value, assumed to be in Jy, which specifies an approximate convolution depth. -IntExcl = 0 # (default=0) : Width of the box (in pixel) to discard from the center. -LogLevel = INFO # (default=INFO) : Logging level in the python logger. -LoopGain = 0.5 # (default=0.1): Loop gain for MORESANE deconvolution. -MajorLoopMiter = 100 # (default=100): Maximum number of iterations allowed in the major loop. Exit condition. -MinorLoopMiter = 30 # (default=30): Maximum number of iterations allowed in the minor loop. Serves as an exit condition when the SNR is does not reach a maximum. -NegComp = False # Enable negative components in model -ScaleCount = None # (default=None): Maximum scale to be considered - maximum scale considered during initialisation. -SigmaLevel = 1 # (default=4): Number of sigma at which thresholding is to be performed. -SingleRun = False # (default=False). If false, will run moresane sequentially scale by scale from startscale to stopscale, else, will use scalecount as a maximum scale. -StartScale = 1 # Starting wavelet scale in IUWT -StopScale = 20 # Stopping wavelet scale in IUWT -SubRegion = None # (default=None): Size, in pixels, of the central region to be analyzed and deconvolved. -Tolerance = 0.75 # (default=0.75): Tolerance level for object extraction. Significant objects contain wavelet coefficients greater than the tolerance multiplied by the maximum wavelet coefficient in the scale under consideration. +_Help = PyMoresane internal options +NMajorIter = 200 # Maximum number of iterations allowed in the major loop. Exit condition. +NMinorIter = 200 # Maximum number of iterations allowed in the minor loop. Serves as an exit condition when the SNR is does not reach a maximum. +Gain = 0.1 # Loop gain for the deconvolution. +ForcePositive = True # Boolean specifier for whether or not a model must be strictly positive. +SigmaCutLevel = 1 # Number of sigma at which thresholding is to be performed. + +[MUFFIN] +_Help = MUFFIN internal options +mu_s = 0.1 # Regularization parameter for spatial regularization +mu_l = 0.2 # Regularization parameter for spectral regularization +nb = (8,0) # nb must be a tuple of wavelets for dwt, or a list of 2 integer for IUWT where the first integer is the scale of decomposition (default:8) and the second integer is the scale or recomposition (default:0) +NMinorIter = 200 # Maximum number of iterations [Log] _Help = Options related to logging @@ -337,6 +330,7 @@ ParsetVersion = 0.2 # parset version number, for migration purposes. + ## NB: Changes introduced for issue #255: ## ## diff --git a/DDFacet/Restore.py b/DDFacet/Restore.py index dfefa16f..a7e11c20 100755 --- a/DDFacet/Restore.py +++ b/DDFacet/Restore.py @@ -28,13 +28,21 @@ from DDFacet.Other import MyLogger from DDFacet.ToolsDir import ModFFTW +from DDFacet.Other import AsyncProcessPool +from DDFacet.Other.AsyncProcessPool import APP, WorkerProcessError +from DDFacet.Other import Multiprocessing + from DDFacet.Other import MyLogger from DDFacet.Other import MyPickle +from DDFacet.ToolsDir.rad2hmsdms import rad2hmsdms log=MyLogger.getLogger("ClassRestoreMachine") - +import scipy.signal +import scipy.stats import multiprocessing NCPU_default=str(int(0.75*multiprocessing.cpu_count())) + + def read_options(): desc="""DDFacet """ @@ -43,16 +51,22 @@ def read_options(): group = optparse.OptionGroup(opt, "* Data selection options") group.add_option('--BaseImageName',help='') group.add_option('--ResidualImage',help='',type="str",default="") - group.add_option('--BeamPix',type='float',help='',default=5) + group.add_option('--BeamPix',type='float',help='',default=5.) group.add_option('--SmoothMode',type='int',help='0 = use default beam, 1 = use smooth beam, default 0',default=0) group.add_option('--MakeCorrected',type='int',help='0 = no normalization correction, 1 = make corrected image, default 1',default=1) group.add_option('--MaskName',type="str",help='',default=5) group.add_option('--NBands',type="int",help='',default=1) group.add_option('--CleanNegComp',type="int",help='',default=0) + group.add_option('--Mode',type="str",help='',default="App") + group.add_option('--RandomCat',type="int",help='',default=0) + group.add_option('--RandomCat_TotalToPeak',type=float,help='',default=1.) + group.add_option('--RandomCat_CountsFile',type=str,help='',default=None) + group.add_option('--ZeroNegComp',type="int",help='',default=0) group.add_option('--DoAlpha',type="int",help='',default=0) group.add_option('--OutName',type="str",help='',default="") group.add_option('--PSFCache',type="str",help='',default="") group.add_option('--NCPU',type="int",help='',default=NCPU_default) + group.add_option('--AddNoise',type=float,help='',default=0) opt.add_option_group(group) options, arguments = opt.parse_args() @@ -61,29 +75,33 @@ def read_options(): return options + + + + class ClassRestoreMachine(): def __init__(self,BaseImageName,BeamPix=5,ResidualImName="",DoAlpha=1, - MaskName="",CleanNegComp=False,NBands=1,OutName="", + MaskName="",CleanNegComp=False, + NBands=1, SmoothMode=0,MakeCorrected=1,options=None): self.DoAlpha=DoAlpha self.BaseImageName=BaseImageName self.BeamPix=BeamPix self.NBands=NBands - self.OutName=OutName + self.OutName=options.OutName self.options=options self.SmoothMode=SmoothMode self.MakeCorrected=MakeCorrected - + self.header_dict={} FileDicoModel="%s.DicoModel"%BaseImageName # ClassModelMachine,DicoModel=GiveModelMachine(FileDicoModel) # self.ModelMachine=ClassModelMachine(Gain=0.1) # self.ModelMachine.FromDico(DicoModel) + print>>log,"Building model machine" ModConstructor = ClassModModelMachine() self.ModelMachine=ModConstructor.GiveInitialisedMMFromFile(FileDicoModel) - - if MaskName!="": self.ModelMachine.CleanMaskedComponants(MaskName) if CleanNegComp: @@ -98,13 +116,15 @@ def __init__(self,BaseImageName,BeamPix=5,ResidualImName="",DoAlpha=1, ResidualImName=FitsFile="%s.app.residual.fits"%BaseImageName else: ResidualImName=FitsFile=ResidualImName + if self.MakeCorrected: if self.SmoothMode: - NormImageName="%s.SmoothNorm.fits"%BaseImageName + NormImageName="%s.MeanSmoothNorm.fits"%BaseImageName else: NormImageName="%s.Norm.fits"%BaseImageName + print>>log,"Reading residual image" self.FitsFile=FitsFile im=image(FitsFile) @@ -119,6 +139,7 @@ def __init__(self,BaseImageName,BeamPix=5,ResidualImName="",DoAlpha=1, nchan,npol,_,_=self.ResidualData.shape testImage=np.zeros_like(self.ResidualData) + print>>log,"Transposing residual..." if ResidualImName!="": for ch in range(nchan): for pol in range(npol): @@ -126,8 +147,10 @@ def __init__(self,BaseImageName,BeamPix=5,ResidualImName="",DoAlpha=1, if self.MakeCorrected: + print>>log,"Reading beam..." SqrtNormImage=np.zeros_like(self.ResidualData) imNorm=image(NormImageName).getdata() + print>>log,"Transposing beam..." for ch in range(nchan): for pol in range(npol): SqrtNormImage[ch,pol,:,:]=np.sqrt(imNorm[ch,pol,:,:].T[::-1,:]) @@ -141,24 +164,39 @@ def __init__(self,BaseImageName,BeamPix=5,ResidualImName="",DoAlpha=1, self.Residual=testImage self.SqrtNormImage=SqrtNormImage + def killWorkers(self): + print>>log, "Killing workers" + APP.terminate() + APP.shutdown() + Multiprocessing.cleanupShm() + def Restore(self): print>>log, "Create restored image" - ModelMachine=self.ModelMachine FWHMFact=2.*np.sqrt(2.*np.log(2.)) + FWHMFact=2.*np.sqrt(2.*np.log(2.)) + BeamPix=self.BeamPix/FWHMFact sigma_x, sigma_y=BeamPix,BeamPix theta=0. bmaj=np.max([sigma_x, sigma_y])*self.CellArcSec*FWHMFact bmin=np.min([sigma_x, sigma_y])*self.CellArcSec*FWHMFact - self.FWHMBeam=(bmaj/3600.,bmin/3600.,theta) + + #bmaj=bmin=0.001666666666666667*3600 + #sigma_x= + self.FWHMBeam=(bmaj/3600./np.sqrt(2.),bmin/3600./np.sqrt(2.),theta) self.PSFGaussPars = (sigma_x*self.CellSizeRad, sigma_y*self.CellSizeRad, theta) + #print "!!!!!!!!!!!!!!!!!!!!" + #self.PSFGaussPars = (BeamPix,BeamPix,0) + + + RefFreq=self.ModelMachine.RefFreq df=RefFreq*0.5 @@ -169,7 +207,7 @@ def Restore(self): import os IdSharedMem=str(int(os.getpid()))+"." - MeanModelImage=ModelMachine.GiveModelImage(RefFreq) + MeanModelImage=self.ModelMachine.GiveModelImage(RefFreq) # #imNorm=image("6SBc.KAFCA.restoredNew.fits.6SBc.KAFCA.restoredNew.fits.MaskLarge.fits").getdata() # imNorm=image("6SB.KAFCA.GA.BIC_00.AP.dirty.fits.mask.fits").getdata() @@ -239,16 +277,78 @@ def Restore(self): ListModelIm=[] #print C/np.array(Lambda) # restored image + for l in Lambda: freq=C/l - print>>log,"Get ModelImage... " - ModelImage=ModelMachine.GiveModelImage(freq) + + if self.options.RandomCat: + print>>log,"Create random catalog... " + ModelImage=self.GiveRandomModelIm() + else: + print>>log,"Get ModelImage... " + ModelImage=self.ModelMachine.GiveModelImage(freq) + + if self.options.ZeroNegComp: + print>>log,"Zeroing negative componants... " + ModelImage[ModelImage<0]=0 ListModelIm.append(ModelImage) - print>>log," ModelImage to apparent flux... " - ModelImage=ModelImage*self.SqrtNormImage + + + if self.options.Mode=="App": + print>>log," ModelImage to apparent flux... " + ModelImage=ModelImage*self.SqrtNormImage print>>log,"Convolve... " print>>log," MinMax = [%f , %f] @ freq = %f MHz"%(ModelImage.min(),ModelImage.max(),freq/1e6) - RestoredImage=ModFFTW.ConvolveGaussian(ModelImage,CellSizeRad=self.CellSizeRad,GaussPars=[self.PSFGaussPars]) + #RestoredImage=ModFFTW.ConvolveGaussianScipy(ModelImage,CellSizeRad=self.CellSizeRad,GaussPars=[self.PSFGaussPars]) + + + if self.options.AddNoise>0.: + print>>log,"Adding Noise... " + ModelImage+=np.random.randn(*ModelImage.shape)*self.options.AddNoise + + RestoredImage,_=ModFFTW.ConvolveGaussianWrapper(ModelImage,Sig=BeamPix) +# ======= +# #RestoredImage,_=ModFFTW.ConvolveGaussianWrapper(ModelImage,Sig=BeamPix) + +# def GiveGauss(Sig0,Sig1): +# npix=20*int(np.sqrt(Sig0**2+Sig1**2)) +# if not npix%2: npix+=1 +# dx=npix/2 +# x,y=np.mgrid[-dx:dx:npix*1j,-dx:dx:npix*1j] +# dsq=x**2+y**2 +# return Sig0**2/(Sig0**2+Sig1**2)*np.exp(-dsq/(2.*(Sig0**2+Sig1**2))) +# R2=np.zeros_like(ModelImage) + +# Sig0=BeamPix/np.sqrt(2.) +# if self.options.RandomCat: +# Sig1=(self.options.RandomCat_SigFactor-1.)*Sig0 +# else: +# Sig1=0. +# nch,npol,_,_=ModelImage.shape +# for ch in range(nch): +# in1=ModelImage[ch,0] +# R2[ch,0,:,:]=scipy.signal.fftconvolve(in1,GiveGauss(Sig0,Sig1), mode='same').real +# RestoredImage=R2 + +# self.header_dict["GSIGMA"]=Sig0 + + +# # print np.max(np.abs(R2-RestoredImage)) +# # import pylab +# # ax=pylab.subplot(1,3,1) +# # pylab.imshow(RestoredImage[0,0],interpolation="nearest") +# # pylab.colorbar() +# # pylab.subplot(1,3,2,sharex=ax,sharey=ax) +# # pylab.imshow(R2[0,0],interpolation="nearest") +# # pylab.colorbar() +# # pylab.subplot(1,3,3,sharex=ax,sharey=ax) +# # pylab.imshow((RestoredImage-R2)[0,0],interpolation="nearest") +# # pylab.colorbar() +# # pylab.show() +# # stop + +# >>>>>>> 0457182a873da89a2758f4be8a18f55cefd88e44 + RestoredImageRes=RestoredImage+self.Residual ListRestoredIm.append(RestoredImageRes) RestoredImageResCorr=RestoredImageRes/self.SqrtNormImage @@ -269,24 +369,33 @@ def Restore(self): ImageName="%s.restoredNew"%self.BaseImageName ImageNameCorr="%s.restoredNew.corr"%self.BaseImageName ImageNameModel="%s.model"%self.BaseImageName + ImageNameModelConv="%s.modelConv"%self.BaseImageName else: ImageName=self.OutName ImageNameCorr=self.OutName+".corr" + ImageNameModel="%s.model"%self.OutName + ImageNameModelConv="%s.modelConv"%self.OutName - CasaImage=ClassCasaImage.ClassCasaimage(ImageNameModel,RestoredImageRes.shape,self.Cell,self.radec)#Lambda=(Lambda0,dLambda,self.NBands)) + CasaImage=ClassCasaImage.ClassCasaimage(ImageNameModel,RestoredImageRes.shape,self.Cell,self.radec,header_dict=self.header_dict)#Lambda=(Lambda0,dLambda,self.NBands)) CasaImage.setdata(ModelImage,CorrT=True) CasaImage.setBeam(self.FWHMBeam) CasaImage.ToFits() CasaImage.close() - CasaImage=ClassCasaImage.ClassCasaimage(ImageName,RestoredImageRes.shape,self.Cell,self.radec)#,Lambda=(Lambda0,dLambda,self.NBands)) + CasaImage=ClassCasaImage.ClassCasaimage(ImageName,RestoredImageRes.shape,self.Cell,self.radec,Freqs=C/np.array(Lambda).ravel(),header_dict=self.header_dict)#,Lambda=(Lambda0,dLambda,self.NBands)) CasaImage.setdata(RestoredImageRes,CorrT=True) CasaImage.setBeam(self.FWHMBeam) CasaImage.ToFits() CasaImage.close() + + CasaImage=ClassCasaImage.ClassCasaimage(ImageNameModelConv,RestoredImage.shape,self.Cell,self.radec,Freqs=C/np.array(Lambda).ravel(),header_dict=self.header_dict)#,Lambda=(Lambda0,dLambda,self.NBands)) + CasaImage.setdata(RestoredImage,CorrT=True) + CasaImage.setBeam(self.FWHMBeam) + CasaImage.ToFits() + CasaImage.close() if self.MakeCorrected: - CasaImage=ClassCasaImage.ClassCasaimage(ImageNameCorr,RestoredImageResCorr.shape,self.Cell,self.radec)#,Lambda=(Lambda0,dLambda,self.NBands)) + CasaImage=ClassCasaImage.ClassCasaimage(ImageNameCorr,RestoredImageResCorr.shape,self.Cell,self.radec,Freqs=C/np.array(Lambda).ravel(),header_dict=self.header_dict)#,Lambda=(Lambda0,dLambda,self.NBands)) CasaImage.setdata(RestoredImageResCorr,CorrT=True) CasaImage.setBeam(self.FWHMBeam) CasaImage.ToFits() @@ -304,7 +413,7 @@ def Restore(self): # Alpha image if self.DoAlpha: print>>log,"Get Index Map... " - IndexMap=ModelMachine.GiveSpectralIndexMap(CellSizeRad=self.CellSizeRad,GaussPars=[self.PSFGaussPars]) + IndexMap=self.ModelMachine.GiveSpectralIndexMap(CellSizeRad=self.CellSizeRad,GaussPars=[self.PSFGaussPars]) ImageName="%s.alphaNew"%self.BaseImageName print>>log," Save... " CasaImage=ClassCasaImage.ClassCasaimage(ImageName,ModelImage.shape,self.Cell,self.radec) @@ -313,6 +422,111 @@ def Restore(self): CasaImage.close() print>>log," Done. " + def GiveRandomModelIm(self): +# np.random.seed(0) + SModel,NModel=np.load(self.options.RandomCat_CountsFile).T + ind=np.argsort(SModel) + SModel=SModel[ind] + NModel=NModel[ind] + #SModel/=1e3 + NModel/=SModel**(5/2.) + def GiveNPerOmega(s): + xp=np.interp(np.log10(s), np.log10(SModel), np.log10(NModel), left=None, right=None) + # import pylab + # pylab.clf() + # pylab.plot(np.log10(SModel), np.log10(NModel)) + # pylab.scatter(np.log10(s),xp) + # pylab.draw() + # pylab.show() + # pylab.pause(0.1) + return 10**xp + + std=np.std(self.Residual.flat[np.int64(np.random.rand(1000)*self.Residual.size)]) + nbin=10000 + smin=2.*std + smax=10. + LogS=np.linspace(np.log10(smin),np.log10(smax),nbin) + + Model=np.zeros_like(self.Residual) + Omega=self.Residual.size*self.CellSizeRad**2 + nx=Model.shape[-1] + im=image(self.FitsFile) + Lra=[] + Ldec=[] + LS=[] + SRA=[] + SDEC=[] + + f,p,_,_=im.toworld((0,0,0,0)) + for iBin in range(nbin-1): + ThisS=(10**LogS[iBin]+10**LogS[iBin+1])/2. + dx=10**LogS[iBin+1]-10**LogS[iBin] + n=int(scipy.stats.poisson.rvs(GiveNPerOmega(ThisS)*Omega*dx))#int(round(GiveNPerOmega(ThisS)*Omega*dx)) + indx=np.array([np.int64(np.random.rand(n)*nx)]).ravel() + indy=np.array([np.int64(np.random.rand(n)*nx)]).ravel() + s0,s1=10**LogS[iBin],10**LogS[iBin+1] + RandS=np.random.rand(n)*(s1-s0)+s0 + Model[0,0,indy,indx]=RandS + for iS in range(indx.size): + _,_,dec,ra=im.toworld((0,0,indy[iS],indx[iS])) + Lra.append(ra) + Ldec.append(dec) + LS.append(RandS[iS]) + #SRA.append(rad2hmsdms(ra,Type="ra").replace(" ",":")) + #SDEC.append(rad2hmsdms(dec,Type="dec").replace(" ",":")) + + #Cat=np.zeros((len(Lra),),dtype=[("ra",np.float64),("StrRA","S200"),("dec",np.float64),("StrDEC","S200"),("S",np.float64)]) + Cat=np.zeros((len(Lra),),dtype=[("ra",np.float64),("dec",np.float64),("S",np.float64)]) + Cat=Cat.view(np.recarray) + Cat.ra=np.array(Lra) + Cat.dec=np.array(Ldec) + #Cat.StrRA=np.array(SRA) + #Cat.StrDEC=np.array(SDEC) + Cat.S=np.array(LS) + CatName="%s.cat.npy"%self.OutName + print>>log,"Saving simulated catalog as %s"%CatName + np.save(CatName,Cat) + + ModelOut=np.zeros_like(Model) + ModelOut[0,0]=Model[0,0].T[::-1] + + p=self.options.RandomCat_TotalToPeak + if p>1.: + nx=101 + x,y=np.mgrid[-nx:nx+1,-nx:nx+1] + r2=x**2+y**2 + def G(sig): + C0=1./(2.*np.pi*sig**2) + C=C0*np.exp(-r2/(2.*sig**2)) + C/=np.sum(C) + return C + ListSig=np.linspace(0.001,10.,100) + TotToPeak=np.array([1./np.max(G(s)) for s in ListSig]) + sig=np.interp(self.options.RandomCat_TotalToPeak,TotToPeak,ListSig) + print>>log,"Found a sig of %f"%sig + Gaussian=G(sig) + ModelOut[0,0]=scipy.signal.fftconvolve(ModelOut[0,0], G(sig), mode='same') + + FWHMFact=2.*np.sqrt(2.*np.log(2.)) + BeamPix=self.BeamPix/FWHMFact + Model=G(sig).reshape((1,1,x.shape[0],x.shape[0])) + ConvModel,_=ModFFTW.ConvolveGaussianWrapper(Model,Sig=BeamPix) + self.SimulObsPeak=np.max(ConvModel) + print>>log," Gaussian Peak: %f"%np.max(Gaussian) + print>>log," Gaussian Int : %f"%np.sum(Gaussian) + print>>log," Obs peak : %f"%self.SimulObsPeak + self.header_dict["OPKRATIO"]=self.SimulObsPeak + self.header_dict["GSIGMA"]=sig + self.header_dict["RTOTPK"]=self.options.RandomCat_TotalToPeak + + + else: + self.header_dict["OPKRATIO"]=1. + self.header_dict["GSIGMA"]=0. + self.header_dict["RTOTPK"]=1. + + + return ModelOut def test(): @@ -332,11 +546,13 @@ def main(options=None): DoAlpha=options.DoAlpha, NBands=options.NBands, CleanNegComp=options.CleanNegComp, - OutName=options.OutName, SmoothMode=options.SmoothMode, MakeCorrected=options.MakeCorrected, options=options) - return CRM.Restore() + CRM.Restore() + CRM.killWorkers() + + diff --git a/DDFacet/Tests/FastUnitTests/TestFitter.py b/DDFacet/Tests/FastUnitTests/TestFitter.py index 3d50ddec..adf20936 100644 --- a/DDFacet/Tests/FastUnitTests/TestFitter.py +++ b/DDFacet/Tests/FastUnitTests/TestFitter.py @@ -105,7 +105,7 @@ def testFitSinc(): cleanBeam = gauss2d(inpars, circle=0, rotate=1, vheight=0)(xx, yy) wnd = 200 psfWnd = psf[(maxAtCrd[0] - wnd):(maxAtCrd[0] + wnd + 1),(maxAtCrd[1] - wnd):(maxAtCrd[1] + wnd + 1)] - (lev, fnull), (bmaj, bmin, th) = fitter.FindSidelobe(psfWnd) + lev, fnull = fitter.FindSidelobe(psfWnd) assert np.isclose(lev,np.sinc(5/2.0),rtol=1e-2,atol=1e-8) # just check if it found the first sidelobe level diff --git a/DDFacet/ToolsDir/ClassAdaptShape.py b/DDFacet/ToolsDir/ClassAdaptShape.py new file mode 100644 index 00000000..b0a25040 --- /dev/null +++ b/DDFacet/ToolsDir/ClassAdaptShape.py @@ -0,0 +1,41 @@ +import numpy as np +from DDFacet.Other import MyLogger +log=MyLogger.getLogger("ClassAdaptShape") +from DDFacet.Other import ModColor + +class ClassAdaptShape(): + def __init__(self,A): + self.A=A + + + def giveOutIm(self,Nout): + A=self.A + nch,npol,Nin,_=A.shape + + if Nin==Nout: + return A.copy() + elif Nin>Nout: + B=np.zeros((nch,npol,Nout,Nout),A.dtype) + print>>log,ModColor.Str(" Output image smaller than input image") + print>>log,ModColor.Str(" adapting %s --> %s"%(str(A.shape),str(B.shape))) + # Input image larger than requested + N0=A.shape[-1] + xc0=yc0=N0/2 + x0d,x1d=xc0-Nout/2,xc0-Nout/2+Nout + s=slice(x0d,x1d) + B[:,:,:,:]=A[:,:,s,s] + return B + else: + B=np.zeros((nch,npol,Nout,Nout),A.dtype) + # Input image smaller - has to zero pad + print>>log,ModColor.Str(" Output image larger than input image") + print>>log,ModColor.Str(" adapting %s --> %s"%(str(A.shape),str(B.shape))) + Na=A.shape[-1] + Nb=B.shape[-1] + xa=Na/2 + xb=Nb/2 + x0d,x1d=xb-xa,xb-xa+Na + s=slice(x0d,x1d) + B[:,:,s,s]=A[:,:,:,:] + return B + diff --git a/DDFacet/ToolsDir/ClassRRGP.py b/DDFacet/ToolsDir/ClassRRGP.py index 8412da5c..50127837 100644 --- a/DDFacet/ToolsDir/ClassRRGP.py +++ b/DDFacet/ToolsDir/ClassRRGP.py @@ -31,7 +31,7 @@ from scipy.linalg import solve_triangular as soltri from scipy import optimize as opt from DDFacet.ToolsDir import ClassGP -import matplotlib.pyplot as plt +#import matplotlib.pyplot as plt import time @@ -241,4 +241,4 @@ def mean_and_cov(self, theta): fbar = np.dot(self.Phip, np.dot(Linv.T, np.dot(Linv, np.dot(self.Phi.T, self.y)))) self.fcovcoeffs = theta[2] ** 2 * np.dot(Linv.T, Linv) covf = theta[2] ** 2 * np.dot(self.Phip, np.dot(Linv.T, np.dot(Linv, self.Phip.T))) - return fbar, covf \ No newline at end of file + return fbar, covf diff --git a/DDFacet/ToolsDir/ClassSpectralFunctions.py b/DDFacet/ToolsDir/ClassSpectralFunctions.py index e9b2be63..2921db6b 100644 --- a/DDFacet/ToolsDir/ClassSpectralFunctions.py +++ b/DDFacet/ToolsDir/ClassSpectralFunctions.py @@ -53,7 +53,6 @@ def setBeamFactors(self): - def GiveBeamFactorsFacet(self,iFacet): @@ -107,6 +106,13 @@ def GiveFreqBandsFluxRatio(self,iFacet,Alpha): return FreqBandsFluxRatio + + def CalcFluxBands(self,NAlpha=21): + Alphas=np.linspace(-1,1,NAlpha) + self.FluxBands=np.zeros((self.NFacets,NAlpha,self.NFreqBand),np.float32) + for iFacet in range(self.NFacets): + self.FluxBands[iFacet]=self.GiveFreqBandsFluxRatio(iFacet,Alphas) + def IntExpFunc(self,S0=1.,Alpha=0.,iChannel=0,iFacet=0): diff --git a/DDFacet/ToolsDir/GiveEdges.py b/DDFacet/ToolsDir/GiveEdges.py index 4b8f9c7e..a690aa39 100644 --- a/DDFacet/ToolsDir/GiveEdges.py +++ b/DDFacet/ToolsDir/GiveEdges.py @@ -20,7 +20,7 @@ import numpy as np -def GiveEdges((xc0,yc0),N0,(xc1,yc1),N1): +def GiveEdges((xc0,yc0),N0,(xc1,yc1),N1,Parity=None): M_xc=xc0 M_yc=yc0 NpixMain=N0 @@ -51,8 +51,28 @@ def GiveEdges((xc0,yc0),N0,(xc1,yc1),N1): y1facet=NpixFacet-dy1 y1main+=1 + IsEven=(lambda x: x%2==0) + + if Parity is not None: + if Parity=="Even": + F=lambda x: IsEven(x) + elif Parity=="Odd": + F=lambda x: not IsEven(x) + + dx=x1main-x0main + dy=y1main-y0main + if F(dx): + x1main-=1 + x1facet-=1 + + if F(dy): + y1main-=1 + y1facet-=1 + Aedge=[x0main,x1main,y0main,y1main] Bedge=[x0facet,x1facet,y0facet,y1facet] + + return Aedge,Bedge diff --git a/DDFacet/ToolsDir/ModFFTW.py b/DDFacet/ToolsDir/ModFFTW.py index e7c2372e..05683d00 100644 --- a/DDFacet/ToolsDir/ModFFTW.py +++ b/DDFacet/ToolsDir/ModFFTW.py @@ -281,24 +281,24 @@ def GiveGauss(Npix,CellSizeRad=None,GaussPars=(0.,0.,0.),dtype=np.float32,parall #Gauss/=np.sum(Gauss) return Gauss -#def ConvolveGaussianScipy(Ain0,Sig=1.,GaussPar=None): -# #warnings.warn("deprecated: this wont work for small ffts...", -# # DeprecationWarning) -# Npix=int(2*8*Sig) -# if Npix%2==0: Npix+=1 -# x0=Npix/2 -# x,y=np.mgrid[-x0:x0:Npix*1j,-x0:x0:Npix*1j] -# #in2=np.exp(-(x**2+y**2)/(2.*Sig**2)) -# if GaussPar is None: -# GaussPar=(Sig,Sig,0) -# in2=Gaussian.Gaussian2D(x,y,GaussPar=GaussPar) -# -# nch,npol,_,_=Ain0.shape -# Out=np.zeros_like(Ain0) -# for ch in range(nch): -# in1=Ain0[ch,0] -# Out[ch,0,:,:]=scipy.signal.fftconvolve(in1, in2, mode='same').real -# return Out,in2 +def ConvolveGaussianScipy(Ain0,Sig=1.,GaussPar=None): + #warnings.warn("deprecated: this wont work for small ffts...", + # DeprecationWarning) + Npix=int(2*8*Sig) + if Npix%2==0: Npix+=1 + x0=Npix/2 + x,y=np.mgrid[-x0:x0:Npix*1j,-x0:x0:Npix*1j] + #in2=np.exp(-(x**2+y**2)/(2.*Sig**2)) + if GaussPar is None: + GaussPar=(Sig,Sig,0) + in2=Gaussian.Gaussian2D(x,y,GaussPar=GaussPar) + + nch,npol,_,_=Ain0.shape + Out=np.zeros_like(Ain0) + for ch in range(nch): + in1=Ain0[ch,0] + Out[ch,0,:,:]=scipy.signal.fftconvolve(in1, in2, mode='same').real + return Out,in2 def ConvolveGaussianWrapper(Ain0,Sig=1.0,GaussPar=None): # a drop-in replacement for ConvolveGaussianScipy which uses diff --git a/DDFacet/ToolsDir/ModFitPSF.py b/DDFacet/ToolsDir/ModFitPSF.py index 74953229..3c7010d1 100644 --- a/DDFacet/ToolsDir/ModFitPSF.py +++ b/DDFacet/ToolsDir/ModFitPSF.py @@ -20,6 +20,8 @@ import numpy as np import gaussfitter2 +import scipy.ndimage.measurements + from DDFacet.Other import MyLogger log= MyLogger.getLogger("FitPSF") @@ -34,7 +36,12 @@ def FitCleanBeam(PSF): sigma_y of gaussian in pixels theta: The rotation angle (radians) of the gaussian from the y axis counter-clockwise (beware of the periodicity) """ - popt=gaussfitter2.gaussfit(PSF, vheight=0, return_all=0, rotate=1) + # zero everything except the main lobe + PSF1 = PSF.copy() + labels = scipy.ndimage.measurements.label(PSF1>=0)[0] + amax = np.unravel_index(np.argmax(PSF), PSF.shape) + PSF1[labels != labels[amax]] = 0 + popt=gaussfitter2.gaussfit(PSF1, vheight=0, return_all=0, rotate=1) amp, xo, yo, sigma_x, sigma_y, theta=popt theta = np.deg2rad(theta) gmaj, gmin = (0, 0) @@ -79,26 +86,18 @@ def FindSidelobe(PSF): """ # Peak of PSF (may be less than unity because of DD terms) # Assume there is only one peak - max_psf = np.max(PSF) - amax_psf = np.argwhere(PSF==max_psf)[0] + amax = np.unravel_index(np.argmax(PSF), PSF.shape) nx,ny=PSF.shape - - # Fit cleanbeam to position - popt = FitCleanBeam(PSF) - - # Create a clean beam: - x = np.arange(0, nx) - y = np.arange(0, ny) - x, y = np.meshgrid(x, y) - bestFit = [max_psf, amax_psf[0], amax_psf[1], popt[0], popt[1], np.rad2deg(popt[2]) - 90] - clean_beam = gaussfitter2.twodgaussian(bestFit, circle=0, rotate=1, vheight=0)(x,y) - - # Create a residual beam without the primary lobe - sidelobes = PSF - clean_beam - amax_sl = np.argwhere(sidelobes == np.max(sidelobes))[0] - sidelobe_dist = np.sqrt(np.sum((amax_sl - amax_psf)**2)) - sidelobe_level = sidelobes[amax_sl[0], amax_sl[1]] - return (sidelobe_level, sidelobe_dist), (popt[0], popt[1], popt[2]) + + PSF1 = PSF.copy() + labels = scipy.ndimage.measurements.label(PSF1>=0)[0] + PSF1[labels == labels[amax]] = 0 + + amax_sl = np.unravel_index(np.argmax(PSF1), PSF1.shape) + sidelobe_dist = np.sqrt((amax_sl[0] - amax[0])**2 + (amax_sl[1] - amax[1])**2) + sidelobe_level = PSF1[amax_sl] + + return sidelobe_level, sidelobe_dist if __name__ == "__main__": from matplotlib import pyplot as plt diff --git a/DDFacet/ToolsDir/gaussfitter2.py b/DDFacet/ToolsDir/gaussfitter2.py index ef1a82ba..5894d8da 100644 --- a/DDFacet/ToolsDir/gaussfitter2.py +++ b/DDFacet/ToolsDir/gaussfitter2.py @@ -35,15 +35,17 @@ def moments (data,circle,rotate,vheight): the gaussian parameters of a 2D distribution by calculating its moments. Depending on the input parameters, will only output a subset of the above""" + data = data.copy() + data[data<0] = 0 #clamp the negatives to 0 to make sure the moments computation doesn't fall over total = data.sum() X, Y = indices(data.shape) x = (X*data).sum()/total y = (Y*data).sum()/total - col = data[:, int(y)].copy() - col[where(col < 0)] = 0.0 #clamp the negatives to 0 to make sure the moments computation doesn't fall over + col = data[:, int(y)] +# col[where(col < 0)] = 0.0 #clamp the negatives to 0 to make sure the moments computation doesn't fall over width_x = sqrt(abs((arange(col.size)-y)**2*col).sum()/col.sum()) - row = data[int(x), :].copy() - row[where(row < 0)] = 0.0 # clamp the negatives to 0 to make sure the moments computation doesn't fall over + row = data[int(x), :] +# row[where(row < 0)] = 0.0 # clamp the negatives to 0 to make sure the moments computation doesn't fall over width_y = sqrt(abs((arange(row.size)-x)**2*row).sum()/row.sum()) width = ( width_x + width_y ) / 2. height = stats.mode(data.ravel())[0][0] if vheight else 0; diff --git a/Dockerfile b/Dockerfile index cc75deda..7646fff6 100644 --- a/Dockerfile +++ b/Dockerfile @@ -48,6 +48,7 @@ ENV DEB_DEPENCENDIES \ python2.7-dev \ libboost-all-dev \ libcfitsio3-dev \ + libhdf5-dev \ wcslib-dev \ libatlas-dev \ liblapack-dev \ diff --git a/setup.py b/setup.py index 755f1093..b4e6c20e 100755 --- a/setup.py +++ b/setup.py @@ -102,6 +102,8 @@ def define_scripts(): }, packages=[pkg, skymodel_pkg], install_requires=[ + "nose >= 1.3.7", + "Cython >= 0.25.2", "numpy >= 1.11.0", "SharedArray >= 2.0.2", "Polygon2 >= 2.0.8", @@ -120,8 +122,9 @@ def define_scripts(): "owlcat >= 1.4.2", "astLib >= 0.8.0", "psutil >= 5.2.2", - "astro-tigger >= 1.3.5", "py-cpuinfo >= 3.2.0", + "tables >= 3.3.0", + "PrettyTable >= 0.7.2" ], include_package_data=True, zip_safe=False,