def zeros(*args,**kwargs):
def ToShared(Name,A):
def DelArray(Name):
def ListNames():
def DelAll(key=None):
def GiveArray(Name):
def DicoToShared(Prefix,Dico,DelInput=False):
def SharedToDico(Prefix):
    def __init__(self,prefixName,Dico):
def SharedObjectToDico(SObject):
def PackListArray(Name,LArray):
def UnPackListArray(Name):
    T=ClassTimeIt.ClassTimeIt("   SharedToDico")
    T=ClassTimeIt.ClassTimeIt("   SharedToDico")
    T=ClassTimeIt.ClassTimeIt("   SharedToDico")
    S=SharedArray.create(Name,(1+NArray+NDimTot+TotSize,),dtype=dS)
    S[0]=NArray
        A=LArray[i]
        A=LArray[i]
    S=GiveArray(Name)
        A=S[idx:idx+size].reshape(shape)
class SharedDicoDescriptor(): class ClassVisServer():
    def __init__(self,MSName,
                 ColName="DATA",
                 TChunkSize=1,
                 TVisSizeMin=1,
                 DicoSelectOptions={},
                 LofarBeam=None,
                 AddNoiseJy=None,IdSharedMem="",
                 SM=None,NCPU=None,
                 Robust=2,Weighting="Natural",
                 WeightUVMinMax=None, WTUV=1.0,
                 GD=None,GDImag=None):

        self.GD=GD
        self.GDImag=GDImag
        self.CalcGridBasedFlags=False
        self.IdSharedMem=IdSharedMem
        PrefixShared="%sSharedVis"%self.IdSharedMem
        self.AddNoiseJy=AddNoiseJy
        self.ReInitChunkCount()
        self.TMemChunkSize=TChunkSize
        self.TVisSizeMin=TVisSizeMin
        self.MSName=MSName
        self.SM=SM
        self.NCPU=NCPU
        self.VisWeights=None
        self.CountPickle=0
        self.ColName=ColName
        self.DicoSelectOptions=DicoSelectOptions
        self.SharedNames=[]
        self.PrefixShared=PrefixShared
        self.VisInSharedMem = (PrefixShared!=None)
        self.LofarBeam=LofarBeam
        self.ApplyBeam=False
        self.DicoClusterDirs_Descriptor=None
        self.Robust=Robust
        self.Weighting=Weighting
        self.DomainsMachine=ClassJonesDomains.ClassJonesDomains()
        self.BeamTimes=np.array([],np.float64)
        self.Init()
        self.dTimesVisMin=self.TVisSizeMin
        self.CurrentVisTimes_SinceStart_Sec=0.,0.
        self.iCurrentVisTime=0
        self.WeightUVMinMax=WeightUVMinMax
        self.WTUV=WTUV
        # self.LoadNextVisChunk()

        # self.TEST_TLIST=[]

    def setSM(self,SM):
        self.SM=SM
        rac,decc=self.MS.radec
        if self.SM.Type=="Catalog":
            self.SM.Calc_LM(rac,decc)

            if self.GD!=None:
                if self.GD["PreApply"]["PreApplySols"][0]!="":
                    CJ=ClassReCluster.ClassReCluster(self.GD)
                    CJ.ReClusterSkyModel(self.SM,self.MS.MSName)
            

    # def SetBeam(self,LofarBeam):
    #     self.BeamMode,self.DtBeamMin,self.BeamRAs,self.BeamDECs = LofarBeam
    #     useArrayFactor=("A" in self.BeamMode)
    #     useElementBeam=("E" in self.BeamMode)
    #     self.MS.LoadSR(useElementBeam=useElementBeam,useArrayFactor=useArrayFactor)
    #     self.ApplyBeam=True
    #     stop

    def Init(self,PointingID=0,NChanJones=1):
        #MSName=self.MDC.giveMS(PointingID).MSName
        kwargs={}
        DecorrMode=""
        ReadUVWDT=False
        if self.GD!=None:
            kwargs["Field"]=self.GD["DataSelection"]["FieldID"]
            kwargs["ChanSlice"]=self.GD["DataSelection"]["ChanSlice"]
            kwargs["DDID"]=self.GD["DataSelection"]["DDID"]
            DecorrMode=self.GD["SkyModel"]["Decorrelation"]
            ReadUVWDT=(("T" in DecorrMode) or ("F" in DecorrMode))
            
        self.ReadUVWDT=ReadUVWDT
        
        ToRaDec=None
        if self.GD is not None and "GDImage" in list(self.GD.keys()):
            ToRaDec=self.GD["GDImage"]["Image"]["PhaseCenterRADEC"]

            if not ToRaDec: # catch empty string
                ToRaDec=None
            
        if ToRaDec=="align":
            log.print(ModColor.Str("kMS does not understand align mode for PhaseCenterRADEC, setting to None..."))
            ToRaDec=None
            #raise RuntimeError("incorrect BeamAt setting: use Facet or Tessel")
        
        MS=ClassMS.ClassMS(self.MSName,Col=self.ColName,DoReadData=False,ReadUVWDT=ReadUVWDT,GD=self.GD,ToRADEC=ToRaDec,**kwargs)
        
        TimesInt=np.arange(0,MS.DTh,self.TMemChunkSize).tolist()
        if not(MS.DTh+1./3600 in TimesInt): TimesInt.append(MS.DTh+1./3600)
        self.TimesInt=TimesInt
        self.NTChunk=len(self.TimesInt)-1
        self.MS=MS
        
        self.DicoMergeStations={}
        
        if self.GD and self.GD["Compression"]["MergeStations"] is not None:
            MergeStations=self.GD["Compression"]["MergeStations"]
            ListMergeNames=[]
            ListMergeStations=[]
            for Name in MergeStations:
                for iAnt in range(MS.na):
                    if Name in MS.StationNames[iAnt]:
                        ListMergeStations.append(iAnt)
                        ListMergeNames.append(MS.StationNames[iAnt])
                        
            log.print("Merging into a single station %s"%str(ListMergeNames))
            self.DicoMergeStations["ListMergeNames"]=ListMergeNames
            self.DicoMergeStations["ListBLMerge"]=[(a0,a1) for a0 in ListMergeStations for a1 in ListMergeStations if a0!=a1]
            self.DicoMergeStations["ListMergeStations"]=ListMergeStations
                        


        ######################################################
        ## Taken from ClassLOFARBeam in DDFacet
        
        ChanWidth=abs(self.MS.ChanWidth.ravel()[0])
        ChanFreqs=self.MS.ChanFreq.flatten()
        if self.GD!=None:
            NChanJones=self.GD["Solvers"]["NChanSols"]
        if NChanJones==0:
            NChanJones=self.MS.NSPWChan
        ChanEdges=np.linspace(ChanFreqs.min()-ChanWidth/2.,ChanFreqs.max()+ChanWidth/2.,NChanJones+1)

        FreqDomains=[[ChanEdges[iF],ChanEdges[iF+1]] for iF in range(NChanJones)]
        FreqDomains=np.array(FreqDomains)
        self.SolsFreqDomains=FreqDomains
        self.NChanJones=NChanJones
        


        MeanFreqJonesChan=(FreqDomains[:,0]+FreqDomains[:,1])/2.
        log.print("Center of frequency domains [MHz]: %s"%str((MeanFreqJonesChan/1e6).tolist()))
        DFreq=np.abs(self.MS.ChanFreq.reshape((self.MS.NSPWChan,1))-MeanFreqJonesChan.reshape((1,NChanJones)))
        self.VisToSolsChanMapping=np.argmin(DFreq,axis=1)
        #log.print(("VisToSolsChanMapping %s"%ListToStr(self.VisToSolsChanMapping)))

        
        self.SolsToVisChanMapping=[]
        for iChanSol in range(NChanJones):
            ind=np.where(self.VisToSolsChanMapping==iChanSol)[0] 
            self.SolsToVisChanMapping.append((ind[0],ind[-1]+1))
        #log.print(("SolsToVisChanMapping %s"%ListToStr(self.SolsToVisChanMapping)))
        

        # ChanDegrid
        FreqBands=ChanEdges
        self.FreqBandsMean=(FreqBands[0:-1]+FreqBands[1::])/2.
        self.FreqBandsMin=FreqBands[0:-1].copy()
        self.FreqBandsMax=FreqBands[1::].copy()
        
        NChanDegrid = NChanJones
        MS=self.MS
        ChanDegridding=np.linspace(FreqBands.min(),FreqBands.max(),NChanDegrid+1)
        FreqChanDegridding=(ChanDegridding[1::]+ChanDegridding[0:-1])/2.
        self.FreqChanDegridding=FreqChanDegridding
        NChanDegrid=FreqChanDegridding.size
        NChanMS=MS.ChanFreq.size
        DChan=np.abs(MS.ChanFreq.reshape((NChanMS,1))-FreqChanDegridding.reshape((1,NChanDegrid)))
        ThisMappingDegrid=np.argmin(DChan,axis=1)
        self.MappingDegrid=ThisMappingDegrid
        #log.print("Mapping degrid: %s"%ListToStr(self.MappingDegrid))
        #NpShared.ToShared("%sMappingDegrid"%self.IdSharedMem,self.MappingDegrid)

        ######################################################
        
        
        # self.CalcWeigths()

        #TimesVisMin=np.arange(0,MS.DTh*60.,self.TVisSizeMin).tolist()
        #if not(MS.DTh*60. in TimesVisMin): TimesVisMin.append(MS.DTh*60.)
        #self.TimesVisMin=np.array(TimesVisMin)




    def CalcWeigths(self,FOV=5.):

        if self.VisWeights!=None: return
        
        uvw,WEIGHT,flags=self.GiveAllUVW()
        u,v,w=uvw.T

        freq=np.mean(self.MS.ChanFreq)
        uvdist=np.sqrt(u**2+v**2)
        uvmax=np.max(uvdist)
        CellSizeRad=res=1./(uvmax*freq/3.e8)
        npix=(FOV*np.pi/180)/res
        
        npix=np.min([npix,30000])

        ImShape=(1,1,npix,npix)
        #VisWeights=WEIGHT[:,0]#np.ones((uvw.shape[0],),dtype=np.float32)
        VisWeights=np.ones((uvw.shape[0],),dtype=np.float32)
        Robust=self.Robust

        # WeightMachine=ClassWeighting.ClassWeighting(ImShape,res)
        # self.VisWeights=WeightMachine.CalcWeights(uvw,VisWeights,Robust=Robust,
        #                                           Weighting=self.Weighting)

        ######################
        #uvw,WEIGHT,flags=self.GiveAllUVW()

        if self.GD is not None:
            if self.GD["Weighting"]["WeightInCol"] is not None and self.GD["Weighting"]["WeightInCol"]!="":
                log.print("Using column %s to compute the weights"%self.GD["Weighting"]["WeightInCol"])
                VisWeights=WEIGHT
            else:
                VisWeights=np.ones((uvw.shape[0],self.MS.ChanFreq.size),dtype=np.float32)
        else:
            VisWeights=np.ones((uvw.shape[0],self.MS.ChanFreq.size),dtype=np.float32)
        
        if np.max(VisWeights)==0.:
            log.print("All imaging weights are 0, setting them to ones")
            VisWeights.fill(1)

        if self.SM.Type=="Image":
            ImShape=self.PaddedFacetShape
            CellSizeRad=self.CellSizeRad

        WeightMachine=ClassWeighting.ClassWeighting(ImShape,CellSizeRad)#res)
        VisWeights=WeightMachine.CalcWeights(uvw,VisWeights,flags,self.MS.ChanFreq,
                                             Robust=Robust,
                                             Weighting=self.Weighting)

        if self.WeightUVMinMax is not None:
            uvmin,uvmax=self.WeightUVMinMax
            log.print('Giving full weight to data in range %f - %f km' % (uvmin, uvmax))
            uvmin*=1000
            uvmax*=1000
            filter=(uvdist<uvmin) | (uvdist>uvmax)
            log.print('Downweighting %i out of %i visibilities' % (np.sum(filter),len(uvdist)))
            VisWeights[filter]*=self.WTUV

        MeanW=np.mean(VisWeights[VisWeights!=0.])
        VisWeights/=MeanW
        log.print('Min weight is %f max is %f' % (np.min(VisWeights),np.max(VisWeights)))
        #VisWeight[VisWeight==0.]=1.
        self.VisWeights=VisWeights

 
    def ReInitChunkCount(self):
        self.CurrentMemTimeChunk=0

    def GiveNextVis(self):

        #log.print( "GiveNextVis")

        t0_bef,t1_bef=self.CurrentVisTimes_SinceStart_Sec
        t0_sec,t1_sec=t1_bef,t1_bef+60.*self.dTimesVisMin

        its_t0,its_t1=self.MS.CurrentChunkTimeRange_SinceT0_sec
        t1_sec=np.min([its_t1,t1_sec])

        self.iCurrentVisTime+=1
        self.CurrentVisTimes_SinceStart_Minutes = t0_sec/60.,t1_sec/60.
        self.CurrentVisTimes_SinceStart_Sec     = t0_sec,t1_sec

        #log.print(("(t0_sec,t1_sec,t1_bef,t1_bef+60.*self.dTimesVisMin)",t0_sec,t1_sec,t1_bef,t1_bef+60.*self.dTimesVisMin))
        #log.print(("(its_t0,its_t1)",its_t0,its_t1))
        #log.print(("self.CurrentVisTimes_SinceStart_Minutes",self.CurrentVisTimes_SinceStart_Minutes))


        #print(self.CurrentVisTimes_SinceStart_Minutes)
        
        if (t0_sec>=its_t1):
            return "EndChunk"

        if not self.have_data:
            return "AllFlaggedThisTime"

        t0_MS=self.MS.F_tstart
        t0_sec+=t0_MS
        t1_sec+=t0_MS
        self.CurrentVisTimes_MS_Sec=t0_sec,t1_sec
        
        D=self.ThisDataChunk

        # Calculate uvw speed for time spearing

        Tmax=self.TimeMemChunkRange_sec_Since70[1]#self.ThisDataChunk["times"][-1]
        # time selection
        indRowsThisChunk=np.where((self.ThisDataChunk["times"]>=t0_sec)&(self.ThisDataChunk["times"]<t1_sec))[0]
        # np.save("indRowsThisChunk.npy",indRowsThisChunk)
        # indRowsThisChunk=np.load("indRowsThisChunk.npy")
        
        if indRowsThisChunk.shape[0]==0:
            if t0_sec>=Tmax:
                return "EndChunk"
            else:
                return "AllFlaggedThisTime"
            
        DATA={}
        DATA["indRowsThisChunk"]=indRowsThisChunk
        for key in D.keys():
            if type(D[key])!=np.ndarray: continue
            if not(key in ['times', 'A1', 'A0', 'flags', 'uvw', 'data', 'Map_VisToJones_Time', "UVW_dt",#"IndexTimesThisChunk", 
                           "W"]):             
                DATA[key]=D[key]
            else:
                DATA[key]=D[key][indRowsThisChunk]

        #############################
        ### data selection
        #############################
        flags=DATA["flags"]
        uvw=DATA["uvw"]
        data=DATA["data"]
        A0=DATA["A0"]
        A1=DATA["A1"]
        times=DATA["times"]
        W=DATA["W"]
        Map_VisToJones_Time=DATA["Map_VisToJones_Time"]
        indRowsThisChunk=DATA["indRowsThisChunk"]
        if self.ReadUVWDT: duvw_dt=DATA["UVW_dt"]

        # IndexTimesThisChunk=DATA["IndexTimesThisChunk"]

        for Field in self.DicoSelectOptions.keys():
            if Field=="UVRangeKm":
                if self.DicoSelectOptions[Field]==None: break
                d0,d1=self.DicoSelectOptions[Field]

                d0*=1e3
                d1*=1e3
                u,v,w=uvw.T
                duv=np.sqrt(u**2+v**2)
                #ind=np.where((duv<d0)|(duv>d1))[0]
                ind=np.where((duv>d0)&(duv<d1))[0]
                
                flags=flags[ind]
                data=data[ind]
                A0=A0[ind]
                A1=A1[ind]
                uvw=uvw[ind]
                times=times[ind]

                #IndexTimesThisChunk=IndexTimesThisChunk[ind]
                W=W[ind]
                Map_VisToJones_Time=Map_VisToJones_Time[ind]
                indRowsThisChunk=indRowsThisChunk[ind]
                if self.ReadUVWDT: duvw_dt=duvw_dt[ind]

        for A in self.FlagAntNumber:
            ind=np.where((A0!=A)&(A1!=A))[0]
            flags=flags[ind]
            data=data[ind]
            A0=A0[ind]
            A1=A1[ind]
            uvw=uvw[ind]
            times=times[ind]
            # IndexTimesThisChunk=IndexTimesThisChunk[ind]
            W=W[ind]
            Map_VisToJones_Time=Map_VisToJones_Time[ind]
            indRowsThisChunk=indRowsThisChunk[ind]
            if self.ReadUVWDT: duvw_dt=duvw_dt[ind]
        
        if self.GD["DataSelection"]["FillFactor"]!=1.:
            Mask=np.random.rand(flags.shape[0])<self.GD["DataSelection"]["FillFactor"]
            ind=np.where(Mask)[0]
            flags=flags[ind]
            data=data[ind]
            A0=A0[ind]
            A1=A1[ind]
            uvw=uvw[ind]
            times=times[ind]
            # IndexTimesThisChunk=IndexTimesThisChunk[ind]
            W=W[ind]
            Map_VisToJones_Time=Map_VisToJones_Time[ind]
            indRowsThisChunk=indRowsThisChunk[ind]
            if self.ReadUVWDT: duvw_dt=duvw_dt[ind]
            
            

        ind=np.where(A0!=A1)[0]
        flags=flags[ind,:,:]
        data=data[ind,:,:]
        A0=A0[ind]
        A1=A1[ind]
        uvw=uvw[ind,:]
        times=times[ind]
        #IndexTimesThisChunk=IndexTimesThisChunk[ind]
        W=W[ind]
        Map_VisToJones_Time=Map_VisToJones_Time[ind]
        indRowsThisChunk=indRowsThisChunk[ind]
        if self.ReadUVWDT: duvw_dt=duvw_dt[ind]
                
        DATA["flags"]=flags
        DATA["rac_decc"]=np.array([self.MS.rac,self.MS.decc])
        DATA["uvw"]=uvw
        DATA["data"]=data
        DATA["A0"]=A0
        DATA["A1"]=A1
        DATA["times"]=times
        #DATA["IndexTimesThisChunk"]=IndexTimesThisChunk
        DATA["W"]=W
        DATA["Map_VisToJones_Time"]=Map_VisToJones_Time
        DATA["indRowsThisChunk"]=indRowsThisChunk
        
        if self.ReadUVWDT: DATA["UVW_dt"]=duvw_dt
                                
        #DATA["UVW_dt"]=self.MS.Give_dUVW_dt(times,A0,A1)
        
        if DATA["flags"].size==0:
            return "AllFlaggedThisTime"
        fFlagged=np.count_nonzero(DATA["flags"])/float(DATA["flags"].size)
        #print fFlagged
        if fFlagged>0.9:
            # log.print( "AllFlaggedThisTime [%f%%]"%(fFlagged*100))
            return "AllFlaggedThisTime"
        
        #if fFlagged==0.:
        #    stop
        # it0=np.min(DATA["IndexTimesThisChunk"])
        # it1=np.max(DATA["IndexTimesThisChunk"])+1
        # DATA["UVW_RefAnt"]=self.ThisDataChunk["UVW_RefAnt"][it0:it1,:,:]
        #
        # # PM=ClassPredict(NCPU=self.NCPU,IdMemShared=self.IdSharedMem)
        # # DATA["Kp"]=PM.GiveKp(DATA,self.SM)

        self.ClearSharedMemory()
        DATA=self.PutInShared(DATA)

        self.SharedVis_Descriptor=NpShared.SharedDicoDescriptor(self.PrefixShared,DATA)

        DATA["A0A1"]=(DATA["A0"],DATA["A1"])

        self.PreApplyJones_Descriptor=None
        if "PreApplyJones" in D.keys():
            NpShared.DicoToShared("%sPreApplyJones"%self.IdSharedMem,D["PreApplyJones"])
            self.PreApplyJones_Descriptor=NpShared.SharedDicoDescriptor("%sPreApplyJones"%self.IdSharedMem,D["PreApplyJones"])


        #it0=np.min(DATA["IndexTimesThisChunk"])
        #it1=np.max(DATA["IndexTimesThisChunk"])+1
        #DATA["UVW_RefAnt"]=self.ThisDataChunk["UVW_RefAnt"][it0:it1,:,:]
        

        #print
        #print self.MS.ROW0,self.MS.ROW1
        #t0=np.min(DATA["times"])-self.MS.F_tstart
        #t1=np.max(DATA["times"])-self.MS.F_tstart
        #self.TEST_TLIST+=sorted(list(set(DATA["times"].tolist())))

        
        return DATA

    def setGridProps(self,Cell,nx):
        self.Cell=Cell
        self.nx=nx
        self.CalcGridBasedFlags=True



    def giveDataSizeAntenna(self):
        t=table(self.MS.MSName,ack=False)
        uvw=t.getcol("UVW")
        flags=t.getcol("FLAG")
        A0,A1=t.getcol("ANTENNA1"),t.getcol("ANTENNA2")
        t.close()
        NVisPerAnt=np.zeros(self.MS.na,np.float64)
        Field="UVRangeKm"
        self.fracNVisPerAnt=np.ones_like(NVisPerAnt)
        NVis=flags[flags==0].size
        if NVis==0:
            log.print( ModColor.Str("Hummm - All the data is flagged!!!"))
            return
        
        if self.DicoSelectOptions[Field] is not None:
            d0,d1=self.DicoSelectOptions[Field]
            
            d0*=1e3
            d1*=1e3
            u,v,w=uvw.T
            duv=np.sqrt(u**2+v**2)
            #ind=np.where((duv<d0)|(duv>d1))[0]
            ind=np.where((duv>d0)&(duv<d1))[0]
            
            flags=flags[ind]
            A0=A0[ind]
            A1=A1[ind]
            uvw=uvw[ind]

            for iAnt in range(self.MS.na):
                NVisPerAnt[iAnt]=np.where((A0==iAnt)|(A1==iAnt))[0].size

            self.fracNVisPerAnt=NVisPerAnt/np.max(NVisPerAnt)
            log.print("Fraction of data per antenna for covariance estimate: %s"%str(self.fracNVisPerAnt.tolist()))


            u,v,w=uvw.T
            d=np.sqrt(u**2+v**2)
            Compactness=np.zeros((self.MS.na,),np.float32)
            for iAnt in range(self.MS.na):
                ind=np.where((A0==iAnt)|(A1==iAnt))[0]
                if ind.size==0:
                    Compactness[iAnt]=1.e-3
                    continue
                Compactness[iAnt]=np.mean(d[ind])
            self.Compactness=Compactness/np.max(Compactness)
            log.print("Compactness: %s"%str(self.Compactness.tolist()))



            
            NVisSel=flags[flags==0].size
            log.print("Total fraction of remaining data after uv-cut: %5.2f %%"%(100*NVisSel/float(NVis)))


    def LoadNextVisChunk(self):
        MS=self.MS

        # bug out when we hit the buffers
        if self.CurrentMemTimeChunk >= self.NTChunk:
            log.print( ModColor.Str("Reached end of observations"))
            self.ReInitChunkCount()
            return "EndOfObservation"

        # get current chunk boundaries
        iT0,iT1=self.CurrentMemTimeChunk,self.CurrentMemTimeChunk+1
        self.CurrentMemTimeChunk+=1

        log.print( "Reading next data chunk in [%5.2f, %5.2f] hours (column %s)"%(self.TimesInt[iT0],self.TimesInt[iT1],MS.ColName))
        self.have_data = MS.ReadData(t0=self.TimesInt[iT0],t1=self.TimesInt[iT1],ReadWeight=True)
        
        if not self.have_data:
            self.CurrentVisTimes_SinceStart_Sec=self.TimesInt[iT0]*3600.,self.TimesInt[iT1]*3600.
            self.CurrentVisTimes_MS_Sec=self.TimesInt[iT0]*3600.+self.MS.F_tstart,self.TimesInt[iT1]*3600.+self.MS.F_tstart
            log.print( ModColor.Str("this data chunk is empty"))
            return "Empty"

        #log.print( "    Rows= [%i, %i]"%(MS.ROW0,MS.ROW1))
        #print float(MS.ROW0)/MS.nbl,float(MS.ROW1)/MS.nbl

        ###############################
        MS=self.MS

        #self.TimeMemChunkRange_sec=MS.times_all[0],MS.times_all[-1]
        self.TimeMemChunkRange_sec=self.TimesInt[iT0]*3600.,self.TimesInt[iT1]*3600.
        self.TimeMemChunkRange_sec_Since70=self.TimesInt[iT0]*3600.+self.MS.F_tstart,self.TimesInt[iT1]*3600.+self.MS.F_tstart
        
        #log.print(("!!!!!!!",self.TimeMemChunkRange_sec))
        times=MS.times_all
        data=MS.data
        A0=MS.A0
        A1=MS.A1
        uvw=MS.uvw
        flags=MS.flag_all
        freqs=MS.ChanFreq.flatten()
        nbl=MS.nbl
        dfreqs=MS.dFreq
        duvw_dt=MS.uvw_dt

        # if Nchan>1:
        #     DoRevertChans=(freqs.flatten()[0]>freqs.flatten()[-1])
        # if self.DoRevertChans:
        #     print ModColor.Str("  ====================== >> Revert Channel order!")
        #     wavelength_chan=wavelength_chan[0,::-1]
        #     freqs=freqs[0,::-1]
        #     self.dFreq=np.abs(self.dFreq)
        
        
        



        #flags.fill(0)

        # f=(np.random.rand(*flags.shape)>0.5)
        # flags[f]=1
        # data[flags]=1e6

        # iAFlag=12
        # ind=np.where((A0==iAFlag)|(A1==iAFlag))[0]
        # flags[ind,:,:]=1
        
        Equidistant=IsChanEquidistant(freqs)
        if freqs.size>1:
            if Equidistant:
                log.print( "Channels are equidistant, can go fast")
            else:
                log.print( ModColor.Str("Channels are not equidistant, cannot go fast"))
        
        MS=self.MS


        if self.SM.Type=="Image":
            u,v,w=uvw.T
            wmax=self.GD["GDImage"]["CF"]["wmax"]
            wmaxkm=wmax/1000.
            log.print( "Flagging baselines with w > %f km"%(wmaxkm))
            C=299792458.
            fmax=self.MS.ChanFreq.ravel()[-1]

            ind=np.where(np.abs(w)>wmax)[0]
            flags[ind,:,:]=1

            f=ind.size/float(flags.shape[0])
            log.print( "  w-Flagged %5.1f%% of the data"%(100*f))



            # data=data[ind]
            # A0=A0[ind]
            # A1=A1[ind]
            # uvw=uvw[ind]
            # times=times[ind]
            # W=W[ind]
            # MapJones=MapJones[ind]
            # indRowsThisChunk=indRowsThisChunk[ind]


        #log.print("::::!!!!!!!!!!!!!!!")
        self.ThresholdFlag=1.#0.9
        self.FlagAntNumber=[]


        ########################################
        
        for A in range(MS.na):
            ind=np.where((MS.A0==A)|(MS.A1==A))[0]
            fA=MS.flag_all[ind].ravel()
            if ind.size==0:
                log.print( "Antenna #%2.2i[%s] is not in the MS"%(A,MS.StationNames[A]))
                self.FlagAntNumber.append(A)
                continue
                
            nf=np.count_nonzero(fA)
            
            Frac=nf/float(fA.size)
            if Frac>self.ThresholdFlag:
                log.print( "Taking antenna #%2.2i[%s] out of the solve (~%4.1f%% of flagged data, more than %4.1f%%)"%\
                    (A,MS.StationNames[A],Frac*100,self.ThresholdFlag*100))
                self.FlagAntNumber.append(A)
                
        if self.CalcGridBasedFlags:
            Cell=self.Cell
            nx=self.nx
            MS=self.MS
            u,v,w=MS.uvw.T
            d=np.sqrt(u**2+v**2)
            CellRad=(Cell/3600.)*np.pi/180
            #_,_,nx,ny=GridShape

            # ###
            # S=CellRad*nx
            # C=3e8
            # freqs=MS.ChanFreq
            # x=d.reshape((d.size,1))*(freqs.reshape((1,freqs.size))/C)*S
            # fA_all=(x>(nx/2))
            # ###
            
            C=3e8
            freqs=MS.ChanFreq.flatten()
            x=d.reshape((d.size,1))*(freqs.reshape((1,freqs.size))/C)*CellRad
            fA_all=(x>(1./2))
            
            for A in range(MS.na):
                ind=np.where((MS.A0==A)|(MS.A1==A))[0]
                fA=fA_all[ind].ravel()
                nf=np.count_nonzero(fA)
                if fA.size==0:
                    Frac=1.0
                else:
                    Frac=nf/float(fA.size)
                if Frac>self.ThresholdFlag:
                    log.print( "Taking antenna #%2.2i[%s] out of the solve (~%4.1f%% of out-grid data, more than %4.1f%%)"%\
                               (A,MS.StationNames[A],Frac*100,self.ThresholdFlag*100))
                    self.FlagAntNumber.append(A)




        if "FlagAnts" in self.DicoSelectOptions.keys():
            FlagAnts=self.DicoSelectOptions["FlagAnts"]
            for Name in FlagAnts:
                for iAnt in range(MS.na):
                    if Name in MS.StationNames[iAnt]:
                        log.print( "Taking antenna #%2.2i[%s] out of the solve"%(iAnt,MS.StationNames[iAnt]))
                        self.FlagAntNumber.append(iAnt)



                        
        if "DistMaxToCore" in self.DicoSelectOptions.keys():
            DMax=self.DicoSelectOptions["DistMaxToCore"]*1e3
            X,Y,Z=MS.StationPos.T
            Xm,Ym,Zm=np.median(MS.StationPos,axis=0).flatten().tolist()
            D=np.sqrt((X-Xm)**2+(Y-Ym)**2+(Z-Zm)**2)
            ind=np.where(D>DMax)[0]

            for iAnt in ind.tolist():
                log.print("Taking antenna #%2.2i[%s] out of the solve (distance to core: %.1f km)"%(iAnt,MS.StationNames[iAnt],D[iAnt]/1e3))
                self.FlagAntNumber.append(iAnt)
            



        #############################
        #############################

        ind=np.where(np.isnan(data))
        flags[ind]=1
        
        # ## debug
        # ind=np.where((A0==0)&(A1==1))[0]
        # flags=flags[ind]
        # data=data[ind]
        # A0=A0[ind]
        # A1=A1[ind]
        # uvw=uvw[ind]
        # times=times[ind]
        # ##


        if self.AddNoiseJy!=None:
            data+=(self.AddNoiseJy/np.sqrt(2.))*(np.random.randn(*data.shape)+1j*np.random.randn(*data.shape))
            stop
        # Building uvw infos
        #################################################
        #################################################
        # log.print( "Building uvw infos .... ")
        # Luvw=np.zeros((MS.times.size,MS.na,3),uvw.dtype)
        # AntRef=0
        # indexTimes=np.zeros((times.size,),np.int64)
        # iTime=0

        # # UVW per antenna
        # #Times_all_32=np.float32(times-MS.times[0])
        # #Times32=np.float32(MS.times-MS.times[0])
        # irow=0
        # for ThisTime in MS.times:#Times32:

        #     T= ClassTimeIt.ClassTimeIt("VS")
        #     T.disable()
        #     #ind=np.where(Times_all_32[irow::]==ThisTime)[0]
        #     ind=np.where(times[irow::]==ThisTime)[0]+irow
            
        #     irow+=ind.size
        #     T.timeit("0b")

        #     indAnt=np.where(A0[ind]==AntRef)[0]
        #     ThisUVW0=uvw[ind][indAnt].copy()
        #     Ant0=A1[ind][indAnt].copy()
        #     T.timeit("1")

        #     indAnt=np.where(A1[ind]==AntRef)[0]
        #     ThisUVW1=-uvw[ind][indAnt].copy()
        #     Ant1=A0[ind][indAnt].copy()
        #     ThisUVW=np.concatenate((ThisUVW1,ThisUVW0[1::]))

        #     T.timeit("2")

        #     #AA=np.concatenate((Ant1,Ant0[1::]))
        #     Luvw[iTime,:,:]=ThisUVW[:,:]
        
        #     T.timeit("3")

        #     #Luvw.append(ThisUVW)
        #     indexTimes[ind]=iTime
        #     iTime+=1
        # log.print( "     .... Done ")
        #################################################
        #################################################


        # Dt_UVW_dt=1.*3600
        # t0=times[0]
        # t1=times[-1]
        # All_UVW_dt=[]
        # Times=np.arange(t0,t1,Dt_UVW_dt).tolist()
        # if not(t1 in Times): Times.append(t1+1)

        # All_UVW_dt=np.array([],np.float32).reshape((0,3))
        # for it in range(len(Times)-1):
        #     t0=Times[it]
        #     t1=Times[it+1]
        #     tt=(t0+t1)/2.
        #     indRows=np.where((times>=t0)&(times<t1))[0]
        #     All_UVW_dt=np.concatenate((All_UVW_dt,self.MS.Give_dUVW_dt(tt,A0[indRows],A1[indRows])))

            

        # UVW_dt=All_UVW_dt



        PredictedData=np.zeros_like(data)
        Indices=np.arange(PredictedData.size).reshape(PredictedData.shape)
        NpShared.ToShared("%sPredictedData"%self.IdSharedMem,PredictedData)
        NpShared.ToShared("%sIndicesData"%self.IdSharedMem,Indices)
        
        PredictedDataGains=np.zeros_like(data)
        IndicesGains=np.arange(PredictedDataGains.size).reshape(PredictedDataGains.shape)
        NpShared.ToShared("%sPredictedDataGains"%self.IdSharedMem,PredictedDataGains)
        NpShared.ToShared("%sIndicesDataGains"%self.IdSharedMem,IndicesGains)
        
        
        
        #NpShared.PackListArray("%sUVW_Ants"%self.IdSharedMem,Luvw)
        #self.UVW_RefAnt=NpShared.ToShared("%sUVW_RefAnt"%self.IdSharedMem,Luvw)
        #self.IndexTimes=NpShared.ToShared("%sIndexTimes"%self.IdSharedMem,indexTimes)
        ThisDataChunk={"times":times,
                       "freqs":freqs,
                       "dfreqs":dfreqs,
                       "freqs_full":freqs,
                       "dfreqs_full":dfreqs,
                       #"A0A1":(A0[ind],A1[ind]),
                       #"A0A1":(A0,A1),
                       "A0":A0,
                       "A1":A1,
                       "uvw":uvw,
                       "flags":flags,
                       "nbl":nbl,
                       "na":MS.na,
                       "data":data,
                       "ROW0":MS.ROW0,
                       "ROW1":MS.ROW1,
                       "infos":np.array([MS.na,MS.TimeInterVal[0]]),
                       #"IndexTimesThisChunk":indexTimes,
                       #"UVW_RefAnt": Luvw,
                       "W":self.VisWeights[MS.ROW0:MS.ROW1],
                       #"IndRows_All_UVW_dt":IndRows_All_UVW_dt,
                       "UVW_dt":duvw_dt
                     }


        self.ThisDataChunk=ThisDataChunk#NpShared.DicoToShared("%sThisDataChunk"%self.IdSharedMem,ThisDataChunk)
        #self.UpdateCompression()
        #self.ThisDataChunk["Map_VisToJones_Time"]=np.zeros(([],),np.int32)
        self.ThisDataChunk["Map_VisToJones_Time"]=np.zeros((times.size,),np.int32)


        ListDicoPreApply=[]
        DoPreApplyJones=False
        if self.GD is not None:
            if self.GD["Beam"]["BeamModel"] is not None:
                
                if self.GD["Beam"]["BeamAt"].lower() == "tessel":
                    log.print("Estimating Beam directions at the center of the tesselated areas")
                    RA,DEC=self.SM.ClusterCat.ra,self.SM.ClusterCat.dec
                elif self.GD["Beam"]["BeamAt"].lower() == "facet":
                    log.print("Estimating Beam directions at the center of the individual facets areas")
                    RA=np.array([self.SM.DicoImager[iFacet]["RaDec"][0] for iFacet in range(len(self.SM.DicoImager))])
                    DEC=np.array([self.SM.DicoImager[iFacet]["RaDec"][1] for iFacet in range(len(self.SM.DicoImager))])
                else:
                    raise RuntimeError("incorrect BeamAt setting: use Facet or Tessel")
                
                if self.GD["Beam"]["BeamModel"]=="LOFAR":
                    NDir=RA.size
                    self.DtBeamMin=self.GD["Beam"]["DtBeamMin"]
                    useArrayFactor=("A" in self.GD["Beam"]["LOFARBeamMode"])
                    useElementBeam=("E" in self.GD["Beam"]["LOFARBeamMode"])
                    self.MS.LoadSR(useElementBeam=useElementBeam,useArrayFactor=useArrayFactor)
                    log.print( "Update LOFAR beam in %i directions [Dt = %3.1f min] ... "%(NDir,self.DtBeamMin))
                    DtBeamSec=self.DtBeamMin*60
                    tmin,tmax=np.min(times)-MS.dt/2.,np.max(times)+MS.dt/2.
                    # TimesBeam=np.arange(np.min(times),np.max(times),DtBeamSec).tolist()
                    # if not(tmax in TimesBeam): TimesBeam.append(tmax)
                    NTimesBeam=round((tmax-tmin)/DtBeamSec)
                    NTimesBeam=int(np.max([2,NTimesBeam]))


                    TimesBeam=np.linspace(np.min(times)-1,np.max(times)+1,NTimesBeam).tolist()
                    TimesBeam=np.array(TimesBeam)

                    T0s=TimesBeam[:-1]
                    T1s=TimesBeam[1:]
                    Tm=(T0s+T1s)/2.
                    
                    self.BeamTimes=TimesBeam
                    # print "!!!!!!!!!!!!!!!!!!!!"
                    # T0s=MS.F_times-MS.dt/2.
                    # T1s=MS.F_times+MS.dt/2.
                    # Tm=MS.F_times

                    # from killMS.Other.rad2hmsdms import rad2hmsdms
                    # for i in range(RA.size): 
                    #     ra,dec=RA[i],DEC[i]
                    #     print rad2hmsdms(ra,Type="ra").replace(" ",":"),rad2hmsdms(dec,Type="dec").replace(" ",".")

                    Beam=np.zeros((Tm.size,NDir,self.MS.na,self.MS.NSPWChan,2,2),np.complex64)
                    for itime in range(Tm.size):
                        ThisTime=Tm[itime]
                        Beam[itime]=self.MS.GiveBeam(ThisTime,RA,DEC)
    
                    # # Beam[:,76,:,:,0,0]=20.
                    # # Beam[:,76,:,:,0,1]=0.
                    # # Beam[:,76,:,:,1,0]=0.
                    # # Beam[:,76,:,:,1,1]=20.
                    # Beam[:,:,:,:,0,0]=20.
                    # Beam[:,:,:,:,0,1]=0.
                    # Beam[:,:,:,:,1,0]=0.
                    # Beam[:,:,:,:,1,1]=20.


                    DicoBeam={}
                    DicoBeam["t0"]=T0s
                    DicoBeam["t1"]=T1s
                    DicoBeam["tm"]=Tm
                    DicoBeam["Jones"]=Beam

                    ChanWidth=self.MS.ChanWidth.ravel()[0]
                    ChanFreqs=self.MS.ChanFreq.flatten()

                    self.DomainsMachine.AddFreqDomains(DicoBeam,ChanFreqs,ChanWidth)
                    
                    NChanBeam=self.GD["Beam"]["NChanBeamPerMS"]
                    if NChanBeam==0:
                        NChanBeam=self.MS.NSPWChan
                    FreqDomainsOut=self.DomainsMachine.GiveFreqDomains(ChanFreqs,ChanWidth,NChanJones=NChanBeam)
                    self.DomainsMachine.AverageInFreq(DicoBeam,FreqDomainsOut)

                    ###### Normalise
                    rac,decc=self.MS.OriginalRadec
                    if self.GD["Beam"]["CenterNorm"]==1:

                        Beam=DicoBeam["Jones"]
                        Beam0=np.zeros((Tm.size,1,self.MS.na,self.MS.NSPWChan,2,2),np.complex64)
                        for itime in range(Tm.size):
                            ThisTime=Tm[itime]
                            Beam0[itime]=self.MS.GiveBeam(ThisTime,np.array([rac]),np.array([decc]))

                            
                        DicoBeamCenter={}
                        DicoBeamCenter["t0"]=T0s
                        DicoBeamCenter["t1"]=T1s
                        DicoBeamCenter["tm"]=Tm
                        DicoBeamCenter["Jones"]=Beam0
                        self.DomainsMachine.AddFreqDomains(DicoBeamCenter,ChanFreqs,ChanWidth)
                        self.DomainsMachine.AverageInFreq(DicoBeamCenter,FreqDomainsOut)
                        Beam0=DicoBeamCenter["Jones"]
                        Beam0inv=ModLinAlg.BatchInverse(Beam0)
                        nt,nd,_,_,_,_=Beam.shape
                        Ones=np.ones((nt,nd, 1, 1, 1, 1),np.float32)
                        Beam0inv=Beam0inv*Ones
                        DicoBeam["Jones"]=ModLinAlg.BatchDot(Beam0inv,Beam)

                    ###### 




                    #nt,nd,na,nch,_,_= Beam.shape
                    #Beam=np.mean(Beam,axis=3).reshape((nt,nd,na,1,2,2))
                    
                    #DicoBeam["ChanMap"]=np.zeros((nch))
                    ListDicoPreApply.append(DicoBeam)

                    DoPreApplyJones=True
                    log.print( "       .... done Update LOFAR beam ")
                elif self.GD["Beam"]["BeamModel"] == "FITS" or self.GD["Beam"]["BeamModel"] == "ATCA":
                    NDir = RA.size
                    self.DtBeamMin = self.GD["Beam"]["DtBeamMin"]

                    if self.GD["Beam"]["BeamModel"] == "FITS":
                        from DDFacet.Data.ClassFITSBeam import ClassFITSBeam as ClassDDFBeam
                    elif self.GD["Beam"]["BeamModel"] == "ATCA":
                        from DDFacet.Data.ClassATCABeam import ClassATCABeam as ClassDDFBeam
                        
                    # make fake opts dict (DDFacet clss expects slightly different option names)
                    opts = self.GD["Beam"]
                    opts["NBand"] = opts["NChanBeamPerMS"]
                    ddfbeam = ClassDDFBeam(self.MS, opts)

                    TimesBeam = np.array(ddfbeam.getBeamSampleTimes(times))
                    FreqDomains = ddfbeam.getFreqDomains()
                    nfreq_dom = FreqDomains.shape[0]

                    log.print( "Update %s beam in %i dirs, %i times, %i freqs ... " % (self.GD["Beam"]["BeamModel"],NDir, len(TimesBeam), nfreq_dom))

                    T0s = TimesBeam[:-1]
                    T1s = TimesBeam[1:]
                    Tm = (T0s + T1s) / 2.

                    self.BeamTimes = TimesBeam

                    Beam = np.zeros((Tm.size, NDir, self.MS.na, FreqDomains.shape[0], 2, 2), np.complex64)
                    for itime, tm in enumerate(Tm):
                        Beam[itime] = ddfbeam.evaluateBeam(tm, RA, DEC)

                    DicoBeam = {}
                    DicoBeam["t0"] = T0s
                    DicoBeam["t1"] = T1s
                    DicoBeam["tm"] = Tm
                    DicoBeam["Jones"] = Beam
                    DicoBeam["FreqDomain"] = FreqDomains

                    ###### Normalise
                    #rac, decc = self.MS.radec
                    rac,decc=self.MS.OriginalRadec
                    if self.GD["Beam"]["CenterNorm"] == 1:

                        Beam = DicoBeam["Jones"]
                        Beam0 = np.zeros((Tm.size, 1, self.MS.na, nfreq_dom, 2, 2), np.complex64)
                        for itime, tm in enumerate(Tm):
                            Beam0[itime] = ddfbeam.evaluateBeam(tm, np.array([rac]), np.array([decc]))

                        DicoBeamCenter = {}
                        DicoBeamCenter["t0"] = T0s
                        DicoBeamCenter["t1"] = T1s
                        DicoBeamCenter["tm"] = Tm
                        DicoBeamCenter["Jones"] = Beam0
                        DicoBeamCenter["FreqDomain"] = FreqDomains
                        Beam0inv = ModLinAlg.BatchInverse(Beam0)
                        nt, nd, _, _, _, _ = Beam.shape
                        Ones = np.ones((nt, nd, 1, 1, 1, 1), np.float32)
                        Beam0inv = Beam0inv * Ones
                        DicoBeam["Jones"] = ModLinAlg.BatchDot(Beam0inv, Beam)

                    ######




                    # nt,nd,na,nch,_,_= Beam.shape
                    # Beam=np.mean(Beam,axis=3).reshape((nt,nd,na,1,2,2))

                    # DicoBeam["ChanMap"]=np.zeros((nch))
                    ListDicoPreApply.append(DicoBeam)

                    DoPreApplyJones = True
                    log.print( "       .... done Update beam ")
                elif self.GD["Beam"]["BeamModel"] == "GMRT":
                    NDir = RA.size
                    self.DtBeamMin = self.GD["Beam"]["DtBeamMin"]

                    from DDFacet.Data.ClassGMRTBeam import ClassGMRTBeam
                    # make fake opts dict (DDFacet clss expects slightly different option names)
                    opts = self.GD["Beam"]
                    opts["NBand"] = opts["NChanBeamPerMS"]
                    gmrtbeam = ClassGMRTBeam(self.MS, opts)

                    TimesBeam = np.array(gmrtbeam.getBeamSampleTimes(times))
                    FreqDomains = gmrtbeam.getFreqDomains()
                    nfreq_dom = FreqDomains.shape[0]

                    log.print( "Update GMRT beam in %i dirs, %i times, %i freqs ... " % (NDir, len(TimesBeam), nfreq_dom))

                    T0s = TimesBeam[:-1]
                    T1s = TimesBeam[1:]
                    Tm = (T0s + T1s) / 2.

                    self.BeamTimes = TimesBeam

                    Beam = np.zeros((Tm.size, NDir, self.MS.na, FreqDomains.shape[0], 2, 2), np.complex64)
                    for itime, tm in enumerate(Tm):
                        Beam[itime] = gmrtbeam.GiveInstrumentBeam(tm, RA, DEC)

                    DicoBeam = {}
                    DicoBeam["t0"] = T0s
                    DicoBeam["t1"] = T1s
                    DicoBeam["tm"] = Tm
                    DicoBeam["Jones"] = Beam
                    DicoBeam["FreqDomain"] = FreqDomains

                    ###### Normalise
                    #rac, decc = self.MS.radec
                    rac,decc=self.MS.OriginalRadec
                    if self.GD["Beam"]["CenterNorm"] == 1:

                        Beam = DicoBeam["Jones"]
                        Beam0 = np.zeros((Tm.size, 1, self.MS.na, nfreq_dom, 2, 2), np.complex64)
                        for itime, tm in enumerate(Tm):
                            Beam0[itime] = gmrtbeam.evaluateBeam(tm, np.array([rac]), np.array([decc]))

                        DicoBeamCenter = {}
                        DicoBeamCenter["t0"] = T0s
                        DicoBeamCenter["t1"] = T1s
                        DicoBeamCenter["tm"] = Tm
                        DicoBeamCenter["Jones"] = Beam0
                        DicoBeamCenter["FreqDomain"] = FreqDomains
                        Beam0inv = ModLinAlg.BatchInverse(Beam0)
                        nt, nd, _, _, _, _ = Beam.shape
                        Ones = np.ones((nt, nd, 1, 1, 1, 1), np.float32)
                        Beam0inv = Beam0inv * Ones
                        DicoBeam["Jones"] = ModLinAlg.BatchDot(Beam0inv, Beam)

                    ######




                    # nt,nd,na,nch,_,_= Beam.shape
                    # Beam=np.mean(Beam,axis=3).reshape((nt,nd,na,1,2,2))

                    # DicoBeam["ChanMap"]=np.zeros((nch))
                    ListDicoPreApply.append(DicoBeam)

                    DoPreApplyJones = True
                    log.print( "       .... done Update GMRT beam ")

            if self.GD["PreApply"]["PreApplySols"][0]!="":
                ModeList=self.GD["PreApply"]["PreApplyMode"]
                if ModeList==[""]: ModeList=["AP"]*len(self.GD["PreApply"]["PreApplySols"])
                for SolFile,Mode in zip(self.GD["PreApply"]["PreApplySols"],ModeList):
                    log.print( "Loading solution file %s in %s mode"%(SolFile,Mode))

                    if (SolFile!="")&(not(".npz" in SolFile)):
                        Method=SolFile
                        ThisMSName=reformat.reformat(os.path.abspath(self.MS.MSName),LastSlash=False)
                        SolFileLoad="%s/killMS.%s.sols.npz"%(ThisMSName,Method)
                        if self.GD["Solutions"]["SolsDir"]:
                            _MSName=reformat.reformat(self.MSName).split("/")[-2]
                            DirName="%s%s"%(reformat.reformat(self.GD["Solutions"]["SolsDir"]),_MSName)
                            SolFileLoad="%s/killMS.%s.sols.npz"%(DirName,SolFile)
                    else:
                        SolFileLoad=SolFile

                    S=np.load(SolFileLoad)
                    Sols=S["Sols"]
                    nt,nch,na,nd,_,_=Sols["G"].shape
                    
                    DicoSols={}
                    DicoSols["t0"]=Sols["t0"]
                    DicoSols["t1"]=Sols["t1"]
                    DicoSols["tm"]=(Sols["t0"]+Sols["t1"])/2.
                    DicoSols["Jones"]=np.swapaxes(Sols["G"],1,3).reshape((nt,nd,na,nch,2,2))
                    DicoSols["FreqDomain"]=S["FreqDomains"]
                    if not("A" in Mode):
                        ind=(DicoSols["Jones"]!=0.)
                        DicoSols["Jones"][ind]/=np.abs(DicoSols["Jones"][ind])
                    if not("P" in Mode):
                        dtype=DicoSols["Jones"].dtype
                        DicoSols["Jones"]=(np.abs(DicoSols["Jones"]).astype(dtype)).copy()

                    #DicoSols["Jones"]=Sols["G"].reshape((nt,nd,na,1,2,2))
                    ListDicoPreApply.append(DicoSols)
                    DoPreApplyJones=True


                    
            
        if DoPreApplyJones:
            DicoJones=ListDicoPreApply[0]
            DomainsMachine=self.DomainsMachine

            if self.SM.Type=="Image":
                DomainsMachine.setFacetToDirMapping([self.SM.DicoImager[iFacet]["iDirJones"] for iFacet in range(len(self.SM.DicoImager))])
            for DicoJones1 in ListDicoPreApply[1::]:
                DicoJones=DomainsMachine.MergeJones(DicoJones1,DicoJones)
            
            DomainsMachine.AddVisToJonesMapping(DicoJones,self.ThisDataChunk["times"],self.ThisDataChunk["freqs"])
            
            
            # ind=np.zeros((times.size,),np.int32)
            # #nt,na,nd,_,_,_=Beam.shape
            # ii=0
            # for it in range(nt):
            #     t0=DicoJones["t0"][it]
            #     t1=DicoJones["t1"][it]
            #     indMStime=np.where((times>=t0)&(times<t1))[0]
            #     indMStime=np.ones((indMStime.size,),np.int32)*it
            #     ind[ii:ii+indMStime.size]=indMStime[:]
            #     ii+=indMStime.size
            # TimeMapping=ind

            #DicoJones["ChanMap"]=self.VisToJonesChanMapping
            #self.ThisDataChunk["MapJones"]=TimeMapping



            self.ThisDataChunk["PreApplyJones"]=DicoJones

            DicoClusterDirs={}
            DicoClusterDirs["l"]=self.SM.ClusterCat.l
            DicoClusterDirs["m"]=self.SM.ClusterCat.m
            DicoClusterDirs["ra"]=self.SM.ClusterCat.ra
            DicoClusterDirs["dec"]=self.SM.ClusterCat.dec
            DicoClusterDirs["I"]=self.SM.ClusterCat.SumI
            DicoClusterDirs["Cluster"]=self.SM.ClusterCat.Cluster
            
            NpShared.DicoToShared("%sDicoClusterDirs"%self.IdSharedMem,DicoClusterDirs)
            self.DicoClusterDirs_Descriptor=NpShared.SharedDicoDescriptor("%sDicoClusterDirs"%self.IdSharedMem,DicoClusterDirs)

            self.ThisDataChunk["Map_VisToJones_Time"]=self.ThisDataChunk["PreApplyJones"]["Map_VisToJones_Time"]

        return "LoadOK"


    def setFOV(self,FullImShape,PaddedFacetShape,FacetShape,CellSizeRad):
        self.FullImShape=FullImShape
        self.PaddedFacetShape=PaddedFacetShape
        self.FacetShape=FacetShape
        self.CellSizeRad=CellSizeRad

    def UpdateCompression(self):
        ThisMSName=reformat.reformat(os.path.abspath(self.MS.MSName),LastSlash=False)

        D=self.ThisDataChunk

        DATA={}
        DATA["flags"]=D["flags"]
        DATA["data"]=D["data"]
        DATA["uvw"]=D["uvw"]
        DATA["A0"]=D["A0"]
        DATA["A1"]=D["A1"]
        DATA["times"]=D["times"]

        GD=self.GD["GDImage"]
        if GD["Compression"]["CompDeGridMode"]:
            MapName="%s/Mapping.CompDeGrid.npy"%ThisMSName
            try:
                FinalMapping=np.load(MapName)
            except:
                if GD["Compression"]["CompDeGridFOV"]=="Facet":
                    _,_,nx,ny=self.FacetShape
                elif GD["Compression"]["CompDeGridFOV"]=="Full":
                    _,_,nx,ny=self.FullImShape
                FOV=self.CellSizeRad*nx*(np.sqrt(2.)/2.)*180./np.pi
                SmearMapMachine=ClassSmearMapping.ClassSmearMapping(self.MS,radiusDeg=FOV,Decorr=(1.-GD["Compression"]["CompDeGridDecorr"]),IdSharedMem=self.IdSharedMem,NCPU=self.NCPU)
                #SmearMapMachine.BuildSmearMapping(DATA)
                FinalMapping,fact=SmearMapMachine.BuildSmearMappingParallel(DATA)
                np.save(MapName,FinalMapping)
                log.print( ModColor.Str("  Effective compression [DeGrid]:   %.2f%%"%fact,col="green"))

            Map=NpShared.ToShared("%sMappingSmearing.DeGrid"%(self.IdSharedMem),FinalMapping)


    def GiveAllUVW(self):
        t=self.MS.GiveMainTable()
        uvw=t.getcol("UVW")
        if self.GD is not None:
            if self.GD["Weighting"]["WeightInCol"] is not None and self.GD["Weighting"]["WeightInCol"]!="":
                WEIGHT=t.getcol(self.GD["Weighting"]["WeightInCol"])
            else:
                WEIGHT=t.getcol("WEIGHT")
        else:
            WEIGHT=t.getcol("WEIGHT")

        F=t.getcol("FLAG")
        t.close()

        return uvw,WEIGHT,F


    def ClearSharedMemory(self):
        NpShared.DelAll(self.PrefixShared)
        NpShared.DelAll("%sDicoData"%self.IdSharedMem)
        NpShared.DelAll("%sKernelMat"%self.IdSharedMem)

        # for Name in self.SharedNames:
        #     NpShared.DelArray(Name)
        self.SharedNames=[]

    def PutInShared(self,Dico):
        DicoOut={}
        for key in Dico.keys():
            if type(Dico[key])!=np.ndarray: continue
            #print "%s.%s"%(self.PrefixShared,key)
            Shared=NpShared.ToShared("%s.%s"%(self.PrefixShared,key),Dico[key])
            DicoOut[key]=Shared
            self.SharedNames.append("%s.%s"%(self.PrefixShared,key))
            
        return DicoOut def invertChol(A):
def invertLU(A):
def sqrtSVD(A,Rank=None):
def BatchInverse(A,H=False):
def BatchH(A):
def BatchDot(A,B):
def BatchDot2(A,B):
def PlotMatSVD(A,s,Ainv):
def invSVD(A):
def SVDw(A):
def EigClean(A):
def Dot_ListBlockMat_Mat(ListBlocks,Mat):
def Dot_ListBlockMat_Mat_Iregular(ListBlocks,Mat):
def test_Dot_ListBlockMat_Mat():
def test_Dot_ListBlockMat_Mat_Big():
    L=np.linalg.cholesky(A)
    A=(A+A.T)/2.
    A=A.reshape((A.size//4,2,2))
    A=A.reshape((A.size//4,2,2))
    AH=np.zeros_like(A)
    AH[:,0,0]=a0
    AH[:,1,1]=d0
    AH[:,0,1]=b0
    AH[:,1,0]=c0
    AH=AH.reshape(shapeOut)
    A=A.reshape((A.size//4,2,2))
    B=B.reshape((B.size//4,2,2))
    C=np.zeros_like(A)
    C00=C[:,0,0]
    C01=C[:,0,1]
    C10=C[:,1,0]
    C11=C[:,1,1]
    C00[:]=a0*a1+c0*b1
    C01[:]=a0*c1+c0*d1
    C10[:]=b0*a1+d0*b1
    C11[:]=b0*c1+d0*d1
    C=C.reshape(shapeOut)
    A=A.reshape((NDir_a,nf,na,2,2))
    B=B.reshape((NDir_b,nf,na,2,2))
    C=np.zeros_like(A)
    C00=C[:,:,:,0,0]
    C01=C[:,:,:,0,1]
    C10=C[:,:,:,1,0]
    C11=C[:,:,:,1,1]
    C00[:,:,:]=a0*a1+c0*b1
    C01[:,:,:]=a0*c1+c0*d1
    C10[:,:,:]=b0*a1+d0*b1
    C11[:,:,:]=b0*c1+d0*d1
    C=C.reshape(shapeOut)
    I=np.dot(Ainv,A)
    B=np.random.randn(nblocks*n,m)
    T=ClassTimeIt.ClassTimeIt()
    R0=Dot_ListBlockMat_Mat(ListBlocks,B)
    R1=np.dot(BlocksMat,B)
    R2=Dot_ListBlockMat_Mat_Iregular(ListBlocks,B)
    B=np.random.randn(nblocks*n,m)
    T=ClassTimeIt.ClassTimeIt()
    R0=Dot_ListBlockMat_Mat(ListBlocks,B) HEADER = '\033[95m'
OKBLUE = '\033[94m'
OKGREEN = '\033[92m'
WARNING = '\033[93m'
FAIL = '\033[91m'
ENDC = '\033[0m'
    HEADER = ''
    OKBLUE = ''
    OKGREEN = ''
    WARNING = ''
    FAIL = ''
    ENDC = ''
def Str(strin0,col="red",Bold=True):
def Sep(strin=None,D=1):
def Title(strin,Big=False):
def disable(): class ClassTimeIt():
    def __init__(self,name="",f=1.):
        self.t0=timemod.time()
        self.IsEnable=True
        if name=="":
            self.name=name
        else:
            self.name=name+": "
        self.IsEnableIncr=False
        self.Counter=""
        self.f=f

    def reinit(self):
        self.t0=timemod.time()

    def timeit(self,stri=" Time",hms=False):
        if self.IsEnable==False: return
        t1=timemod.time()
        dt=self.f*(t1-self.t0)
        if not(hms):
            Sout= "  * %s%s %s : %7.5fs"%(self.name,stri,str(self.Counter),dt)
            if self.IsEnableIncr: self.Counter+=1
        else:
            ss=(dt)/60.
            m=int(ss)
            s=(ss-m)*60.
            Sout= "  * %s computation time: %i min. %4.1f sec."%(stri,m,s)
        self.t0=t1
        if DoLog:
            log.print( Sout)
        else:
            print(Sout)

        return dt

    def timeitHMS(self,stri=" Time"):
        t1=timemod.time()
        self.t0=t1

    def disable(self):
        self.IsEnable=False

    def enableIncr(self,incr=1):
        self.IsEnableIncr=True
        self.Counter=0

    def AddDt(self,Var):
        t1=timemod.time()
        dt=t1-self.t0
        Var=Var+dt
        self.t0=t1
        return Var class Counter():
    def __init__(self,N=1):
        if N==0: N=1
        self.N=N
        self.i=-1

    def __call__(self):
        i=self.i
        #print "%i -> %i"%(i,i+1)
        self.i+=1
        Cond=(self.i%self.N==0)
        return Cond class ClassModelEvolution():
    def __init__(self,iAnt,iChanSol,WeightType="exp",WeigthScale=1,order=1,StepStart=5,BufferNPoints=10,sigQ=0.01,DoEvolve=True,IdSharedMem=""):
        self.WeightType=WeightType 
        self.WeigthScale=WeigthScale*60. #in min
        self.order=order
        self.StepStart=StepStart
        self.iAnt=iAnt
        self.iChanSol=iChanSol
        self.BufferNPoints=BufferNPoints
        self.sigQ=sigQ
        self.DoEvolve=DoEvolve
        self.IdSharedMem=IdSharedMem

    def Evolve0(self,Gin,Pa,kapa=1.):
        done=NpShared.GiveArray("%sSolsArray_done"%self.IdSharedMem)
        
        
        indDone=np.where(done==1)[0]
        #print kapa
        #print type(NpShared.GiveArray("%sSharedCovariance_Q"%self.IdSharedMem))
        Q=kapa*NpShared.GiveArray("%sSharedCovariance_Q"%self.IdSharedMem)[self.iChanSol,self.iAnt]
        #print indDone.size
        #print "mean",np.mean(Q)


        ##########
        # Ptot=Pa+Q
        # #nt,_,_,_=Gin.shape
        # #print Gin.shape
        # g=Gin
        # gg=g.ravel()
        # #gg+=(np.random.randn(*gg.shape)+1j*np.random.randn(*gg.shape))*np.sqrt(np.diag(Ptot))/np.sqrt(2.)

        # # print Pa.shape,Q.shape
        # # print np.diag(Pa)
        # # print np.diag(Q)
        # # print np.diag(Ptot)
        # # print

        # return Ptot
        ##############

        # return Pa+Q
        if indDone.size<2: return Pa+Q
        
        # #########################
        # take scans where no solve has been done into account
        G=NpShared.GiveArray("%sSolsArray_G"%self.IdSharedMem)[indDone][:,self.iChanSol,self.iAnt,:,0,0]
        Gm=np.mean(G,axis=-1)
        dG=Gm[:-1]-Gm[1:]
        done0=NpShared.GiveArray("%sSolsArray_done"%self.IdSharedMem)
        done1=np.zeros((done0.size,),int)
        done1[1:1+dG.size]=(dG!=0)
        done1[0]=1
        done=(done0&done1)
        indDone=np.where(done==1)[0]
        if indDone.size<2: return Pa+Q
        # #########################
        
        t0=NpShared.GiveArray("%sSolsArray_t0"%self.IdSharedMem)[indDone]
        t1=NpShared.GiveArray("%sSolsArray_t1"%self.IdSharedMem)[indDone]
        tm=NpShared.GiveArray("%sSolsArray_tm"%self.IdSharedMem)[indDone]


        G=NpShared.GiveArray("%sSolsArray_G"%self.IdSharedMem)[indDone][:,self.iChanSol,self.iAnt,:,:,:]

        
        nt,nd,npolx,npoly=G.shape

        #if nt<=self.StepStart: return None

        if nt>self.BufferNPoints:
            G=G[-self.BufferNPoints::,:,:,:]
            tm=tm[-self.BufferNPoints::]

        G=G.copy()

        nt,_,_,_=G.shape
        NPars=nd*npolx*npoly
        G=G.reshape((nt,NPars))

        F=np.ones((NPars,),G.dtype)
        PaOut=np.zeros_like(Pa)

        tm0=tm.copy()
        tm0=np.abs(tm-tm[-1])
        w=np.exp(-tm0/self.WeigthScale)
        w/=np.sum(w)
        w=w[::-1]

        for iPar in range(NPars):
            #g_t=G[:,iPar][-1]
            #ThisG=Gin.ravel()[iPar]
            #ratio=1.+(ThisG-g_t)/g_t
            g_t=G[:,iPar]
            ThisG=Gin.ravel()[iPar]
            #ratio=1.+np.std(g_t)
            #norm=np.max([np.abs(np.mean(g_t))
            #ratio=np.cov(g_t)/Pa[iPar,iPar]
            #print np.cov(g_t),Pa[iPar,iPar],ratio
            ratio=np.abs(ThisG-np.mean(g_t))/np.sqrt(Pa[iPar,iPar]+Q[iPar,iPar])
            
            #ratio=np.abs(g_t[-1]-np.mean(g_t))/np.sqrt(Pa[iPar,iPar])#+Q[iPar,iPar]))
            diff=np.sum(w*(ThisG-g_t))/np.sum(w)
            ratio=np.abs(diff)/np.sqrt(Pa[iPar,iPar]+Q[iPar,iPar])
            F[iPar]=1.#ratio#/np.sqrt(2.)

            diff=ThisG-g_t[-1]#np.sum(w*(ThisG-g_t))/np.sum(w)
            #PaOut[iPar,iPar]=np.abs(diff)**2+Pa[iPar,iPar]+Q[iPar,iPar]
            PaOut[iPar,iPar]=np.abs(diff)**2+Pa[iPar,iPar]+Q[iPar,iPar]

        # Q=np.diag(np.ones((PaOut.shape[0],)))*(self.sigQ**2)

        PaOut=F.reshape((NPars,1))*Pa*F.reshape((1,NPars)).conj()+Q
        #print(F)
        #print(Q)
        # stop
        return PaOut
        


    def Evolve(self,xEst,Pa,CurrentTime):
        done=NpShared.GiveArray("%sSolsArray_done"%self.IdSharedMem)
        indDone=np.where(done==1)[0]
        Q=NpShared.GiveArray("%sSharedCovariance_Q"%self.IdSharedMem)[self.iAnt]

        t0=NpShared.GiveArray("%sSolsArray_t0"%self.IdSharedMem)[indDone]
        t1=NpShared.GiveArray("%sSolsArray_t1"%self.IdSharedMem)[indDone]
        tm=NpShared.GiveArray("%sSolsArray_tm"%self.IdSharedMem)[indDone]


        G=NpShared.GiveArray("%sSolsArray_G"%self.IdSharedMem)[indDone][:,self.iAnt,:,:,:]

        
        nt,nd,npol,_=G.shape

        if nt<=self.StepStart: return None,None

        if nt>self.BufferNPoints:
            G=G[-nt::,:,:,:]
            tm=tm[-nt::]

        G=G.copy()
        tm0=tm.copy()
        tm0=tm-tm[-1]
        ThisTime=CurrentTime-tm[-1]

        nt,_,_,_=G.shape
        NPars=nd*npol*npol
        G=G.reshape((nt,NPars))


        Gout=np.zeros((nd*npol*npol),dtype=G.dtype)
        Gout[:]=G[-1]


        F=np.ones((NPars,),G.dtype)
        if self.DoEvolve:
            if self.WeightType=="exp":
                w=np.exp(-tm0/self.WeigthScale)
                w/=np.sum(w)
                w=w[::-1]
            dx=1e-6
            for iPar in range(NPars):
                g_t=G[:,iPar]
                g_r=g_t.real.copy()
                g_i=g_t.imag.copy()
                
                ####
                z_r0 = np.polyfit(tm0, g_r, self.order, w=w)
                z_i0 = np.polyfit(tm0, g_i, self.order, w=w)
                poly_r = np.poly1d(z_r0)
                poly_i = np.poly1d(z_i0)
                x0_r=poly_r(ThisTime)
                x0_i=poly_i(ThisTime)
                Gout[iPar]=x0_r+1j*x0_i
                
                ####
                g_r[-1]+=dx
                g_i[-1]+=dx
                z_r1 = np.polyfit(tm0, g_r, self.order, w=w)
                z_i1 = np.polyfit(tm0, g_i, self.order, w=w)
                poly_r = np.poly1d(z_r1)
                poly_i = np.poly1d(z_i1)
                x1_r=poly_r(ThisTime)
                x1_i=poly_i(ThisTime)
                
                # dz=((x0_r-x1_r)+1j*(x0_i-x1_i))/dx

                xlast=G[-1][iPar]
                dz=((x0_r-xlast.real)+1j*(x0_i-xlast.imag))/np.sqrt((Pa[iPar,iPar]+Q[iPar,iPar]))
                F[iPar]=dz/np.sqrt(2.)

            # if self.iAnt==0:
            #     xx=np.linspace(tm0.min(),tm0.max(),100)
            #     pylab.clf()
            #     pylab.plot(tm0, g_r)
            #     pylab.plot(xx, poly_r(xx))
            #     pylab.scatter([ThisTime],[x1_r])
            #     pylab.draw()
            #     pylab.show(False)
            #     pylab.pause(0.1)
            #     print F
            
        # if self.iAnt==0:
        #     pylab.clf()
        #     pylab.imshow(np.diag(F).real,interpolation="nearest")
        #     pylab.draw()
        #     pylab.show(False)
        #     pylab.pause(0.1)


        
        #Pa=P[self.iAnt]
        PaOut=np.zeros_like(Pa)
        #Q=np.diag(np.ones((PaOut.shape[0],)))*(self.sigQ**2)
        PaOut=F.reshape((NPars,1))*Pa*F.reshape((1,NPars)).conj()+Q
        
        
        
        Gout=Gout.reshape((nd,npol,npol))
        print(np.diag(PaOut))

        return Gout,PaOut class ClassSolverLM(ClassJacobianAntenna):
    def __init__(self, *args, **kwargs):
        ClassJacobianAntenna.__init__(self, *args, **kwargs)

    def PrepareJHJ_LM(self):

        self.L_JHJinv=[]
        if self.DataAllFlagged:
            return

        for ipol in range(self.NJacobBlocks_X):

            M=self.L_JHJ[ipol]
            if self.DoTikhonov:
                self.LambdaTkNorm=self.LambdaTk*np.mean(np.abs(np.diag(M)))
                
                # Lin.shape= (self.NDir,self.NJacobBlocks_X,self.NJacobBlocks_Y)
                Linv=np.diag(self.Linv[:,ipol,:].ravel())
                Linv*=self.LambdaTkNorm/(1.+self.LambdaLM)
                M2=M+Linv

                # pylab.clf()
                # pylab.subplot(1,2,1)
                # pylab.imshow(np.abs(M),interpolation="nearest")
                # pylab.colorbar()
                # pylab.subplot(1,2,2)
                # pylab.imshow(np.abs(Linv),interpolation="nearest")
                # pylab.colorbar()
                # pylab.draw()
                # pylab.show(False)
                # pylab.pause(0.1)
                # stop

            else:
                M2=M


            JHJinv=ModLinAlg.invSVD(M2)
            #JHJinv=ModLinAlg.invSVD(self.JHJ)
            self.L_JHJinv.append(JHJinv)

    def doLMStep(self,Gains):
        # if self.iAnt==55:
        #     print(self.iAnt,Gains)
        T=ClassTimeIt.ClassTimeIt("doLMStep")
        T.disable()

        #Gains.fill(1.)
        
#         A=np.random.randn(10000,100)+1j*np.random.randn(10000,100)
#         B=np.random.randn(10000,100)+1j*np.random.randn(10000,100)
#         AT=A.T#.conj().copy()
# #        AT=A.T
#         # A=np.require(A,requirements='F_CONTIGUOUS')
#         # AT=np.require(AT,requirements='F_CONTIGUOUS')
#         # A=np.require(A,requirements='F')
#         # AT=np.require(AT,requirements='F')


#         T=ClassTimeIt.ClassTimeIt("doLMStep")
#         for i in range(20):
#             np.dot(AT,B)

#         T.timeit("%i"%i)

        
        if not(self.HasKernelMatrix):
            self.CalcKernelMatrix()
            self.SelectChannelKernelMat()
            T.timeit("CalcKernelMatrix")

        Ga=self.GiveSubVecGainAnt(Gains)


        if self.DoCompress:
            flags_key="flags_flat_avg"
            data_key="data_flat_avg"
            if self.DoMergeStations:
                flags_key="flags_flat_avg_merged"
                data_key="data_flat_avg_merged"
        else:
            flags_key="flags_flat"
            data_key="data_flat"
            
        f=(self.DicoData[flags_key]==0)
        
        # ind=np.where(f)[0]
        # if self.iAnt==56:
        #     print ind.size/float(f.size),np.abs(Gains[self.iAnt,0,0,0])

        
        if self.DataAllFlagged:
            return Ga.reshape((self.NDir,self.NJacobBlocks_X,self.NJacobBlocks_Y)),None,{"std":-1.,"max":-1.,"kapa":None}



        # if ind.size==0:
        #     return Ga.reshape((self.NDir,self.NJacobBlocks_X,self.NJacobBlocks_Y)),None,{"std":-1.,"max":-1.,"kapa":None}


        z=self.DicoData[data_key]#self.GiveDataVec()
        

        self.CalcJacobianAntenna(Gains)
        T.timeit("CalcJacobianAntenna")
        self.PrepareJHJ_LM()
        T.timeit("PrepareJHJ_L")



        T.timeit("GiveSubVecGainAnt")
        Jx=self.J_x(Ga)
        T.timeit("Jx")
        zr=z-Jx

        zr[self.DicoData[flags_key]]=0
        T.timeit("resid")

        # JH_z_0=np.load("LM.npz")["JH_z"]
        # x1_0=np.load("LM.npz")["x1"]
        # z_0=np.load("LM.npz")["z"]
        # Jx_0=np.load("LM.npz")["Jx"]




        InfoNoise={"std":np.std(zr[f]),"max":np.max(np.abs(zr[f])),"kapa":None}


        JH_z=self.JH_z(zr)
        T.timeit("JH_z")
        #self.JHJinv=ModLinAlg.invSVD(self.JHJ)
        #self.JHJinv=np.linalg.inv(self.JHJ)
        xi=Ga.flatten()
        T.timeit("self.JHJinv_x")
        

        if self.DoTikhonov:
            self.LambdaTkNorm
            Gi=xi.reshape((self.NDir,self.NJacobBlocks_X,self.NJacobBlocks_Y))
            JH_z=JH_z.reshape((self.NDir,self.NJacobBlocks_X,self.NJacobBlocks_Y))
            for polIndex in range(self.NJacobBlocks_X):
                gireg=Gi[:,polIndex,:]
                #gi=JH_z[:,polIndex,:]
                x0reg=self.X0[:,polIndex,:]
                Linv=(self.Linv[:,polIndex,:])
                JH_z[:,polIndex,:]-=self.LambdaTkNorm*Linv*(gireg-x0reg)
        dx = (1./(1.+self.LambdaLM)) * self.JHJinv_x(JH_z)

        
        
        # #print self.iAnt
        # if True:#self.iAnt==55:
        #     f=(self.DicoData[flags_key]==0)
        #     import pylab
        #     fig=pylab.figure(1)
        #     op0=np.abs
        #     op1=np.real
        #     pylab.clf()
        #     pylab.subplot(1,3,1)
        #     #pylab.plot(op0(z[f])[::1]**2)#[::11])
        #     pylab.plot( op1( z[f])[::1] )#[::11])
        #     #pylab.ylim(0,800)
        #     pylab.subplot(1,3,2)
        #     #pylab.plot(op0(Jx[f])[::1]**2)#[::11])
        #     pylab.plot( op1(Jx[f])[::1] )#[::11])
        #     #pylab.ylim(0,800)
        #     pylab.subplot(1,3,3)
        #     #pylab.plot(op0(zr[f])[::1])#[::11])
        #     pylab.plot( op1(zr[f])[::1] )#[::11])
        #     #pylab.ylim(-30,30)
        #     pylab.draw()
        #     pylab.show(block=False)
        #     pylab.pause(0.1)
            
        #     # iF=0
        #     # while True:
        #     #     fName="Graph_%i_%i.png"%(self.iAnt,iF)
        #     #     import os
        #     #     if not os.path.isfile(fName):
        #     #         break
        #     #     else:
        #     #         iF+=1
        #     # fig.savefig(fName)
            


        # # pylab.figure(2)
        # # pylab.clf()
        # # #pylab.plot((z)[::11])
        # # #pylab.plot((Jx-Jx_0)[::11])
        # # #pylab.plot(zr[::11])
        # # #pylab.plot(JH_z.flatten())
        # # #pylab.plot(JH_z_0.flatten())
        # # pylab.plot(x1.flatten())
        # # pylab.plot(x1_0.flatten())
        # # pylab.draw()
        # # pylab.show(False)
        # # pylab.pause(0.1)

        # stop
        # # np.savez("LM",JH_z=JH_z,x1=x1,z=z,Jx=Jx)
 
        # print JH_z.shape

        dx+=xi
        del(self.LJacob)
        T.timeit("rest")
        # print self.iAnt,np.mean(x1),x1.size,ind.size
        
        xout=dx.reshape((self.NDir,self.NJacobBlocks_X,self.NJacobBlocks_Y))
        # print self.iAnt,xout.ravel()
        return xout,None,InfoNoise class ClassSolverEKF(ClassJacobianAntenna):
    def __init__(self, *args, **kwargs):
        ClassJacobianAntenna.__init__(self, *args, **kwargs)

    def GivePaPol(self,Pa_in,ipol):
        PaPol=Pa_in.reshape((self.NDir,self.NJacobBlocks_X,self.NJacobBlocks_Y,self.NDir,self.NJacobBlocks_X,self.NJacobBlocks_Y))
        PaPol=PaPol[:,ipol,:,:,ipol,:].reshape((self.NDir*self.NJacobBlocks_Y,self.NDir*self.NJacobBlocks_Y))
        return PaPol

    def setQxInvPol(self):
        QxInv=np.ones((self.NDir,1,self.NJacobBlocks_Y,self.NDir,1,self.NJacobBlocks_Y),np.float32)
        QxInv*=1./(self.AmpQx**2)
        QxInv=QxInv.reshape((self.NDir*self.NJacobBlocks_Y,self.NDir*self.NJacobBlocks_Y))
        self.LQxInv=[QxInv,QxInv]

    def PrepareJHJ_EKF(self,Pa_in,rms):
        self.L_JHJinv=[]
        incr=1
        # pylab.figure(1)
        # pylab.clf()
        # pylab.imshow(np.abs(self.JHJ),interpolation="nearest")
        # pylab.draw()
        # pylab.show(False)
        # pylab.pause(0.1)

        for ipol in range(self.NJacobBlocks_X):
            PaPol=self.GivePaPol(Pa_in,ipol)
            Pinv=ModLinAlg.invSVD(PaPol)
            JHJ=self.L_JHJ[ipol]#*(1./rms**2)
            JHJ+=Pinv
            if self.DoReg:
                JHJ+=self.LQxInv[ipol]*(self.gamma**2)
            JHJinv=ModLinAlg.invSVD(JHJ)
            self.L_JHJinv.append(JHJinv)




    def CalcKapa_i(self,yr,Pa,rms):
        kapaout=0
        for ipol in range(self.NJacobBlocks_X):
            J=self.LJacob[ipol]
            PaPol=self.GivePaPol(Pa,ipol)
            pa=np.abs(np.diag(PaPol))
            pa=pa.reshape(1,pa.size)
            JP=J*pa
            trJPJH=np.sum(np.abs(JP*J.conj()))
            trYYH=np.sum(np.abs(yr)**2)
            Np=np.where(self.DicoData["flags_flat"]==0)[0].size
            Take=(self.DicoData["flags_flat"]==0)
            trR=np.sum(self.R_flat[Take])#Np*rms**2
            kapa=np.abs((trYYH-trR)/trJPJH)
            kapaout+=np.sqrt(kapa)
            #if self.iAnt==0:
            #    print("old",self.iAnt,rms,np.sqrt(kapa),trYYH,trR,trJPJH,pa)

        kapaout=np.max([1.,kapaout])
        return kapaout


    def CalcKapa_i_new(self,yr,Pa,rms):
        kapaout=0
        T=ClassTimeIt.ClassTimeIt("    Kapa")
        T.disable()
        iT=0
        for ipol in range(self.NJacobBlocks_X):
            J=self.LJacob[ipol]
            PaPol=self.GivePaPol(Pa,ipol)
            pa=np.abs(np.diag(PaPol))
            pa=pa.reshape(1,pa.size)
            T.timeit(iT); iT+=1
            nrow,_=J.shape
            flags=(self.DicoData["flags_flat"][ipol]==0)
            T.timeit(iT); iT+=1

            Weigths=self.Weights_flat[ipol].reshape((nrow,1))
            T.timeit(iT); iT+=1

            Jw=Weigths*J
            JP=Jw*pa
            T.timeit(iT); iT+=1

            trJPJH=np.sum(np.abs(JP[flags]*Jw[flags].conj()))
            T.timeit(iT); iT+=1

            YYH=np.abs(yr[ipol,flags])**2

            T.timeit(iT); iT+=1
            # Np=np.where(self.DicoData["flags_flat"]==0)[0].size
            # Take=(self.DicoData["flags_flat"]==0)
            
            T.timeit(iT); iT+=1
            R=(self.R_flat[ipol][flags])#Np*rms**2
            
            ww=(Weigths[flags]).ravel()
            
            trYYH_R=np.sum(ww**2*(YYH-R))
            T.timeit(iT); iT+=1
            kapa=np.abs(trYYH_R/trJPJH)
            #kapa=1
            kapaout+=np.sqrt(kapa)
            #if self.iAnt==0:
            #    print("new",self.iAnt,rms,np.sqrt(kapa),trYYH_R,trJPJH,pa)
        kapaout=np.max([1.,kapaout])

        return kapaout

    def ApplyK_vec(self,zr,rms,Pa,DoReg=True):

        

        # pylab.figure(1)
        # pylab.clf()
        # incr=1
        # pylab.subplot(2,1,1)
        # pylab.plot((z.flatten())[::incr].real)
        # pylab.plot((Jx.flatten())[::incr].real)
        # pylab.plot(zr.flatten()[::incr].real)
        # pylab.subplot(2,1,2)
        # pylab.plot((z.flatten())[::incr].imag)
        # pylab.plot((Jx.flatten())[::incr].imag)
        # pylab.plot(zr.flatten()[::incr].imag)
        # pylab.draw()
        # pylab.show(False)
        # pylab.pause(0.1)
        # stop


        # pylab.figure(1)
        # pylab.clf()
        # incr=1
        # pylab.subplot(2,1,1)
        # pylab.imshow(np.abs(self.JHJ))
        # pylab.subplot(2,1,2)
        # pylab.imshow(np.abs(self.JHJinv))
        # pylab.draw()
        # pylab.show(False)
        # pylab.pause(0.1)
        # stop


        # if self.iAnt==5:
        #     pylab.figure(1)
        #     pylab.clf()
        #     incr=1
        #     f=(self.DicoData["flags_flat"]==0)
        #     pylab.subplot(2,1,1)
        #     pylab.plot((z1[f].flatten())[::incr].real)
        #     pylab.plot((zr[f].flatten())[::incr].real)
        #     pylab.plot((z1-zr)[f].flatten()[::incr].real)
        #     pylab.subplot(2,1,2)
        #     pylab.plot((z1[f].flatten())[::incr].imag)
        #     pylab.plot((zr[f].flatten())[::incr].imag)
        #     pylab.plot((z1-zr)[f].flatten()[::incr].imag)
        #     pylab.draw()
        #     pylab.show(False)
        #     pylab.pause(0.1)

        Rinv_zr=self.Rinv_flat*zr
        JH_z=self.JH_z(Rinv_zr)

        x1 = self.JHJinv_x(JH_z)
        if (self.DoReg)&(DoReg):
            self.G0=np.ones_like(self.Ga)#*200.
            dx1a=self.Msq_x(self.LQxInv,(self.Ga-self.G0))
            dx1b = self.gamma*self.JHJinv_x(dx1a)
            print("x'_0:",dx1b)
            x1+=dx1b
            print("x':",x1)

        z1 = self.J_x(x1)

        zr-=z1
        zr*=self.Rinv_flat
        x2=self.JH_z(zr)

        if (self.DoReg)&(DoReg):
            xr=self.Ga.ravel()-self.G0.ravel()-self.gamma*x1.ravel()
            dx2=self.gamma*self.Msq_x(self.LQxInv,xr)
            x2+=dx2.reshape(x2.shape)

        x3=[]
        for ipol in range(self.NJacobBlocks_X):
            PaPol=self.GivePaPol(Pa,ipol)
            #print(PaPol,PaPol.shape)


            if self.TypeDot=="Numpy":
                Prod=np.dot(PaPol,x2[:,ipol,:].flatten())
            elif self.TypeDot=="SSE":
                X2=x2[:,ipol,:].flatten()
                X2=X2.reshape((1,X2.size))
                Prod=NpDotSSE.dot_A_BT(PaPol.copy(),X2)
            

            x3.append(Prod.reshape((self.NDir,1,self.NJacobBlocks_Y)))

            

        x3=np.concatenate(x3,axis=1)
        #x3=np.swapaxes(x3,1,2)

        return x3
    

        
    def doEKFStep(self,Gains,P,evP,rms,Gains0Iter=None):
        T=ClassTimeIt.ClassTimeIt("    EKF")
        T.disable()
        if not(self.HasKernelMatrix):
            self.CalcKernelMatrix(rms)
            self.SelectChannelKernelMat()
            T.timeit("CalcKernelMatrix")

        # print(self.iAnt,"MMMM",np.max(Gains),np.max(P),np.max(evP))
        # print(self.iAnt,"MMMM",np.max(Gains),np.max(P),np.max(evP))
        # print(self.iAnt,"MMMM",np.max(Gains),np.max(P),np.max(evP))

            
        z=self.DicoData["data_flat"]#self.GiveDataVec()

        f=(self.DicoData["flags_flat"]==0)
        ind=np.where(f)[0]
        Pa=P[self.iAnt]
        Ga=self.GiveSubVecGainAnt(Gains)
        self.Ga=Ga
        self.rms=rms

        #if self.iAnt==1:
        #    print(evP.ravel())
        self.rmsFromData=None
        if ind.size==0 or self.DataAllFlagged or self.ZeroSizedData:
            return Ga.reshape((self.NDir,self.NJacobBlocks_X,self.NJacobBlocks_Y)),Pa,{"std":-1.,"max":-1.,"kapa":-1.},0
        
        if self.DoReg:
            self.setQxInvPol()
        
        self.CalcJacobianAntenna(Gains)
        T.timeit("Jacob")
        
        # if Gains0Iter!=None:
        #     Ga=self.GiveSubVecGainAnt(Gains0Iter)

        Jx=self.J_x(Ga)
        T.timeit("J_x")

        

        self.PrepareJHJ_EKF(Pa,rms)
        T.timeit("PrepareJHJ")

        # estimate x
        zr=(z-Jx)
        zr[self.DicoData["flags_flat"]]=0
        
        T.timeit("Resid")
        
        kapa=self.CalcKapa_i(zr,Pa,rms)
        
        # try:
        #     kapa=self.CalcKapa_i_new(zr,Pa,rms)
        #     self.DicoData["Ga"]=Ga
        #     self.DicoData["zr"]=zr
        #     self.DicoData["rms"]=rms
        #     self.DicoData["KernelMat_AllChan"]=self.KernelMat_AllChan
        #     #killMS.Other.MyPickle.Save(self.DicoData,"DicoData_%i.pickle"%self.iAnt)
        #     #print(np.array([1])/0)
        # except:
        #     self.DicoData["Ga"]=Ga
        #     self.DicoData["zr"]=zr
        #     self.DicoData["rms"]=rms
        #     self.DicoData["Pa"]=Pa
        #     self.DicoData["KernelMat_AllChan"]=self.KernelMat_AllChan
        #     killMS.Other.MyPickle.Save(self.DicoData,"DicoData_crash_%i.pickle"%self.iAnt)
        #     stop

        # Weighted std estimate 
        zrs=zr[f]
        ws=np.absolute(self.DicoData["Rinv_flat"][f])
        std=np.sqrt(np.sum(ws*np.absolute(zrs)**2)/np.sum(ws))

        # # Original std estimate 
        # std=np.std(zr[f])


        InfoNoise={"std":std,"max":np.max(np.abs(zr[f])),"kapa":kapa}
        #print(self.iAnt,InfoNoise)
        #T.timeit("kapa")

        self.rmsFromData=np.std(zr[f])
        T.timeit("rmsFromData")

        # if np.isnan(self.rmsFromData):
        #     print(zr)
        #     print(zr[f])
        #     print(self.rmsFromData)
        #     stop

        # if self.iAnt==51:
        #     #self.DicoData["flags_flat"].fill(0)
        #     f=(self.DicoData["flags_flat"]==0)
        #     fig=pylab.figure(2)
        #     pylab.clf()
        #     pylab.plot((z[f]))#[::11])#[::11])
        #     pylab.plot((Jx[f]))#[::11])#[::11])

        #     #pylab.plot(zr[f])#[::11])#[::11])
        #     #pylab.draw()
        #     ifile=0
        #     while True:
        #         fname="png/png.%5.5i.png"%ifile
        #         if os.path.isfile(fname) :
        #             fig.savefig(fname)
        #             break
        #         ifile+=1
                
                
        #     #pylab.show(False)
        #     #pylab.pause(0.1)
        #     #stop


        x3=self.ApplyK_vec(zr,rms,Pa)
        
        T.timeit("ApplyK_vec")
        x0=Ga.flatten()
        x4=x0+self.LambdaKF*x3.flatten()

        # estimate P

        #Pa_new1=Pa-np.dot(evPa,Pa)
        evPa=evP[self.iAnt]
        Pa_new1=np.dot(evPa,Pa)
        #Pa_new1=Pa


        T.timeit("EstimateP")
        # ##################
        # for iPar in range(Pa.shape[0]):
        #     J_Px=self.J_x(Pa[iPar,:])
        #     xP=self.ApplyK_vec(J_Px,rms,Pa)
        #     evPa[iPar,:]=xP.flatten()
        # evPa= Pa-evPa
        # Pa_new1=evPa

        del(self.LJacob)
        T.timeit("Rest")
        
        # if self.iAnt==0:
        #     print(x4,Pa_new1,InfoNoise,evPa,Pa)

        # if np.max(np.abs(x4))>10:
        #     self.DicoData["Ga"]=Ga
        #     self.DicoData["Jx"]=Jx
        #     self.DicoData["zr"]=zr
        #     self.DicoData["rms"]=rms
        #     self.DicoData["Pa"]=Pa
        #     self.DicoData["x4"]=x4
        #     self.DicoData["ch0ch1"]=(self.ch0,self.ch1)
        #     self.DicoData["KernelMat_AllChan"]=self.KernelMat_AllChan
        #     killMS.Other.MyPickle.Save(self.DicoData,"DicoData_diverge_%i.pickle"%self.iAnt)
        #     stop

        return x4.reshape((self.NDir,self.NJacobBlocks_X,self.NJacobBlocks_Y)),Pa_new1,InfoNoise,1


    def CalcMatrixEvolveCov(self,Gains,P,rms):
        #print("EVOLVE!!!!!!!!!")
        if not(self.HasKernelMatrix):
            self.CalcKernelMatrix(rms)
            self.SelectChannelKernelMat()
        if self.LQxInv==None:
            self.setQxInvPol()
#            self.CalcKernelMatrix(rms)
        self.CalcJacobianAntenna(Gains)
        Pa=P[self.iAnt]
        NPars=Pa.shape[0]
        if self.DataAllFlagged:
            # print(self.iAnt,"OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO")
            # print(self.iAnt,"OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO")
            # print(self.iAnt,"OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO")
            PaOnes=np.diag(np.ones((NPars,),self.CType))
            return PaOnes
        
        self.PrepareJHJ_EKF(Pa,rms)
        PaOnes=np.diag(np.ones((NPars,),self.CType))

        evPa=np.zeros_like(Pa)

        for iPar in range(Pa.shape[0]):
            J_Px=self.J_x(PaOnes[iPar,:])
            xP=self.ApplyK_vec(J_Px,rms,Pa,DoReg=False)
            evPa[iPar,:]=xP.flatten()


        evPa= PaOnes-evPa#(np.diag(np.diag(Pa-Pa_new)))#Pa-Pa_new#np.abs(np.diag(np.diag(Pa-Pa_new)))
        evPa=np.diag(np.diag(evPa))
        # print("CalcMatrixEvolveCov",self.iAnt,np.diag(evPa).max())
        # print("CalcMatrixEvolveCov",self.iAnt,np.diag(evPa).max())
        # print("CalcMatrixEvolveCov",self.iAnt,np.diag(evPa).max())
        # print("CalcMatrixEvolveCov",self.iAnt,np.diag(evPa).max())


        #print(evPa.min(),evPa.real.min())
        #print("=========Ev",self.iAnt,evPa)
        #if self.iAnt==0:
        #    print(evPa,Gains.ravel(),P.ravel())
        return evPa class ClassSolPredictMachine():
    def __init__(self,GD):
        self.GD=GD
        # FileName=self.GD["KAFCA"]["EvolutionSolFile"]
        # if not ".npz" in FileName:
        #     FileName="%s/killMS.%s.sols.npz"%(self.GD["VisData"]["MSName"],FileName)


        SolsDir=GD["Solutions"]["SolsDir"]
        SolsName=self.GD["KAFCA"]["EvolutionSolFile"]
        MSName=os.path.abspath(self.GD["VisData"]["MSName"])
        if SolsDir is None:
            FileName="%skillMS.%s.sols.npz"%(reformat.reformat(options.MSName),SolsName)
        else:
            _MSName=reformat.reformat(MSName).split("/")[-2]
            DirName=os.path.abspath("%s%s"%(reformat.reformat(SolsDir),_MSName))
            if not os.path.isdir(DirName):
                os.makedirs(DirName)
            FileName="%s/killMS.%s.sols.npz"%(DirName,SolsName)




        log.print( "Reading solution file %s"%FileName)
        self.DicoSols=np.load(FileName)
        self.Sols=self.DicoSols["Sols"]
        self.Sols=self.Sols.view(np.recarray)
        self.FreqDomains=self.DicoSols["FreqDomains"]
        self.MeanFreqDomains=np.mean(self.FreqDomains,axis=1)
        self.t0=self.Sols.t0
        self.t1=self.Sols.t1
        self.tmean=(self.t0+self.t1)/2.
        self.G=self.Sols.G
        self.ClusterCat=self.DicoSols["ClusterCat"]
        self.ClusterCat=self.ClusterCat.view(np.recarray)
        self.raSol=self.ClusterCat.ra
        self.decSol=self.ClusterCat.dec

        # nt,nch,na,nd,_,_=self.G.shape
        

    def GiveClosestSol(self,t,freq_domains,ants,ras,decs):
        d=D(ras,self.raSol,decs,self.decSol)
        ind_dir=np.argmin(d,axis=1)
        
        mean_freqs=np.mean(freq_domains,axis=1)
        dfreq=D(mean_freqs,self.MeanFreqDomains)
        ind_freq=np.argmin(dfreq,axis=1)
        
        dtime=D(np.array([t]),self.tmean)
        ind_time=np.argmin(dtime,axis=1)

        # print ind_time
        # print self.G.shape
        # # nChan,na,nd,npolx,npoly

        nch=ind_freq.size
        ndir=ind_dir.size
        na=ants.size
        Gout=np.zeros((nch,na,ndir,2,2),np.complex64)
        ind_pol=np.arange(4)

        # Gout0=np.zeros((nch,na,ndir,2,2),np.complex64)
        # for ch in ind_freq.ravel():
        #     for ant in ants.ravel():
        #         for d in ind_dir.ravel():
        #             for polx in range(2):
        #                 for poly in range(2):
        #                     Gout0[ch,ant,d,polx,poly]=self.G[ind_time[0]][ch,ant,d,polx,poly]

        ind_time=ind_time.reshape((-1,1,1,1,1))
        ind_freq=ind_freq.reshape((1,-1,1,1,1))
        ind_ant=ants.reshape((1,1,-1,1,1))
        ind_dir=ind_dir.reshape((1,1,1,-1,1))
        ind_pol=ind_pol.reshape((1,1,1,1,-1))

        nt_G,nch_G,na_G,nd_G,_,_=self.G.shape
        index=ind_time*nch_G*na_G*nd_G*2*2 + ind_freq*na_G*nd_G*2*2 + ind_ant*nd_G*2*2 + ind_dir*2*2 + ind_pol
        Gout.flat[:]=self.G.flat[index.ravel()]

        

        # # ind_freq=ind_freq.reshape((-1,1,1,1))
        # # ind_ant=ants.reshape((1,-1,1,1))
        # # ind_dir=ind_dir.reshape((1,1,-1,1))
        # # ind_pol=ind_pol.reshape((1,1,1,-1))
        # # index= ind_freq*na*ndir*2*2 + ind_ant*ndir*2*2 + ind_dir*2*2 + ind_pol
        # # Gout.flat[:]=self.G[ind_time[0]].flat[index.ravel()]
        # print np.max(Gout0-Gout)
        # print
        # print
        # print
        # #Gout=Gout0
        
        return Gout class ClassPredict():
    def __init__(self,Precision="S",NCPU=6,IdMemShared=None,DoSmearing="",BeamAtFacet=False,LExp=None,LSinc=None):
        self.NCPU=NCPU
        ne.set_num_threads(self.NCPU)
        if Precision=="D":
            self.CType=np.complex128
            self.FType=np.float64
        if Precision=="S":
            self.CType=np.complex64
            self.FType=np.float32
        self.DoSmearing=DoSmearing
        self.IdSharedMem=IdMemShared
        self._BeamAtFacet = BeamAtFacet
        
        Np=2
        if self.DoSmearing!=0:
            if (("F" in DoSmearing)or("T" in DoSmearing)): Np=100000


            
        if LExp==None:
            x=np.linspace(0.,10,Np)
            Exp=np.float32(np.exp(-x))
            LExp=[Exp,x[1]-x[0]]
        self.LExp=LExp

        if LSinc==None:
            x=np.linspace(0.,10,Np)
            Sinc=np.zeros(x.shape,np.float32)
            Sinc[0]=1.
            Sinc[1::]=np.sin(x[1::])/(x[1::])
            LSinc=[Sinc,x[1]-x[0]]
            
            #phi=1.471034
            #d=int(phi/(x[1]-x[0]))
            #print Sinc[d],np.sin(phi)/(phi)
            #stop
        self.LSinc=LSinc



    def ApplyCal(self,DicoData,ApplyTimeJones,iCluster):
        D=ApplyTimeJones
        Jones=D["Jones"]
        JonesH=D["JonesH"]
        lt0,lt1=D["t0"],D["t1"]
        ColOutDir=DicoData["data"]
        A0=DicoData["A0"]
        A1=DicoData["A1"]
        times=DicoData["times"]
        na=int(DicoData["infos"][0])

        
        for it in range(lt0.size):
            
            t0,t1=lt0[it],lt1[it]
            ind=np.where((times>=t0)&(times<t1))[0]

            if ind.size==0: continue
            data=ColOutDir[ind]
            # flags=DicoData["flags"][ind]
            A0sel=A0[ind]
            A1sel=A1[ind]
            
            if "Map_VisToJones_Freq" in ApplyTimeJones.keys():
                ChanMap=ApplyTimeJones["Map_VisToJones_Freq"]
            else:
                ChanMap=range(nf)

            for ichan in range(len(ChanMap)):
                JChan=ChanMap[ichan]
                if iCluster!=-1:
                    J0=Jones[it,iCluster,:,JChan,:,:].reshape((na,4))
                    JH0=JonesH[it,iCluster,:,JChan,:,:].reshape((na,4))
                else:
                    J0=np.mean(Jones[it,:,:,JChan,:,:],axis=1).reshape((na,4))
                    JH0=np.mean(JonesH[it,:,:,JChan,:,:],axis=1).reshape((na,4))
                
                J=ModLinAlg.BatchInverse(J0)
                JH=ModLinAlg.BatchInverse(JH0)

                data[:,ichan,:]=ModLinAlg.BatchDot(J[A0sel,:],data[:,ichan,:])
                data[:,ichan,:]=ModLinAlg.BatchDot(data[:,ichan,:],JH[A1sel,:])

                # Abs_g0=(np.abs(J0[A0sel,0])<Threshold)
                # Abs_g1=(np.abs(JH0[A1sel,0])<Threshold)
                # flags[Abs_g0,ichan,:]=True
                # flags[Abs_g1,ichan,:]=True

            ColOutDir[ind]=data[:]

            # DicoData["flags"][ind]=flags[:]


    def GiveParamJonesList(self,DicoJonesMatricesIn,A0,A1):
        if "DicoApplyJones" in DicoJonesMatricesIn.keys():
            DicoJonesMatrices=DicoJonesMatricesIn["DicoApplyJones"]
        else:
            DicoJonesMatrices=DicoJonesMatricesIn
        JonesMatrices=np.complex64(DicoJonesMatrices["Jones"])
        Map_VisToJones_Time=np.int32(DicoJonesMatrices["Map_VisToJones_Time"])
        # print DicoJonesMatrices.keys()
        # print DicoJonesMatricesIn.keys()

        Map_VisToJones_Freq=np.int32(DicoJonesMatrices["Map_VisToJones_Freq"])
        #MapJones=np.int32(np.arange(A0.shape[0]))
        #print DicoJonesMatrices["MapJones"].shape
        #stop
        A0=np.int32(A0)
        A1=np.int32(A1)
        #print Map_VisToJones_Time
        #print Map_VisToJones_Freq
        ParamJonesList=[Map_VisToJones_Time,A0,A1,JonesMatrices,Map_VisToJones_Freq]
        # print sorted(list(set(ParamJonesList[0].tolist())))
        return ParamJonesList

    def GiveCovariance(self,DicoData,ApplyTimeJones,SM):
        D=ApplyTimeJones
        Beam=D["Beam"]
        BeamH=D["BeamH"]
        lt0,lt1=D["t0"],D["t1"]
        A0=DicoData["A0"]
        A1=DicoData["A1"]
        times=DicoData["times"]
        na=DicoData["infos"][0]


        #print "predict...",times[0],times[-1]
        #Predict=self.predictKernelPolCluster(DicoData,SM,ApplyTimeJones=ApplyTimeJones)
        #print "done..."


        Resid=DicoData["resid"]#-Predict
        flags=DicoData["flags"]

        ListDirection=SM.Dirs

        #import pylab
        #pylab.clf()
        #import ppgplot

        MaxMat=np.zeros(Resid.shape,dtype=np.float32)
        CVis=np.zeros_like(Resid)

        


        for iCluster in ListDirection:
            CVis.fill(0)
            ParamJonesList=self.GiveParamJonesList(ApplyTimeJones,A0,A1)
            ParamJonesList=ParamJonesList+[iCluster]
            predict.CorrVis(Resid,CVis,ParamJonesList)

            CVis[flags==1]=0.
            aCVis=np.abs(CVis)
            
            #print "In direction %i: (std, std_abs, med_abs, max_abs)=(%f, %f, %f, %f)"%(iCluster,np.std(CVis),np.std(aCVis),np.median(aCVis),np.max(aCVis))
            ind=(aCVis>MaxMat)
            MaxMat[ind]=aCVis[ind]

            del(aCVis)


        nrow,nch,_=Resid.shape
        W=DicoData["W"]#np.ones((nrow,nch),np.float32)
        rms=findrms.findrms(MaxMat)
        diff=np.abs(MaxMat-np.median(MaxMat))/rms
        cond=(diff>3.)
        ind=np.any(cond,axis=2)
        W[ind]=0.


    def GiveResidAntCovariance(self,DicoData,ApplyTimeJones,SM):
        D=ApplyTimeJones
        Jones=D["Jones"]
        JonesH=D["JonesH"]
        lt0,lt1=D["t0"],D["t1"]
        ColOutDir=DicoData["data"]
        A0=DicoData["A0"]
        A1=DicoData["A1"]
        times=DicoData["times"]
        na=DicoData["infos"][0]
        Sigma=ApplyTimeJones["Stats"]
        W=DicoData["W"]#np.ones((nrow,nch),np.float32)
        nrow,nch=W.shape

        # print Jones.shape
        rmsAllAnts=Sigma[:,:,:,0]


        rmsAllAnts=rmsAllAnts[rmsAllAnts>0.]
        
        if len(rmsAllAnts)==0:
            W.fill(1)
            return
        rms=np.min(rmsAllAnts)
        #print rmsAllAnts,rms
        S=Sigma[:,:,:,0]
        S[S==0]=1e6
        S[S==-1]=1e6


        for it in range(lt0.size):
            
            t0,t1=lt0[it],lt1[it]
            ind=np.where((times>=t0)&(times<t1))[0]

            if ind.size==0: continue
            data=ColOutDir[ind]
            # flags=DicoData["flags"][ind]

            # print ind,t0,t1

            A0sel=A0[ind]
            A1sel=A1[ind]
            
            ThisStats=np.sqrt(np.abs(Sigma[it][0,:,0]**2-rms**2))
            # try:
            #     ThisStats=np.sqrt(np.abs(Sigma[it][:,0]**2-rms**2))
            # except:
            #     print "S",Sigma[it][:,0]
            #     print "rms",rms
            #     print "Resid",(Sigma[it][:,0]**2-rms**2)
            #     #ApplyTimeJones["Stats"]

            # print "t:",it
            # print "0",np.min(Sigma[it][:,0]**2-rms**2)
            # print "1",np.min(Sigma[it][:,0]**2)
            # print "2",(rms**2)

            # nt,nd,na,1,2,2
            Jabs=np.abs(Jones[it,:,:,0,0,0])
            J=np.mean(Jabs,axis=0)
            # print "================="
            # print A0sel
            # print Jones.shape,Jabs.shape,J.shape
            J0=J[A0sel]
            J1=J[A1sel]
            
            sig0=ThisStats[A0sel]
            sig1=ThisStats[A1sel]

            fact=1.
            w=1./(rms**2+(sig0)**(2*fact)+(sig1)**(2*fact))#+sig0*sig1)
            for ich in range(nch):
                W[ind,ich]=w[:]




    # def GiveGM(self,iFacet,SM):
    #     GridMachine=ClassDDEGridMachine.ClassDDEGridMachine(SM.GD,
    #                                                         SM.DicoImager[iFacet]["DicoConfigGM"]["ChanFreq"],
    #                                                         SM.DicoImager[iFacet]["DicoConfigGM"]["Npix"],
    #                                                         lmShift=SM.DicoImager[iFacet]["lmShift"],
    #                                                         IdSharedMem=self.IdSharedMem,
    #                                                         IDFacet=SM.DicoImager[iFacet]["IDFacet"],
    #                                                         SpheNorm=False)

    def GiveGM(self,iFacet,SM):
        """
        Factory: Initializes a gridding machine for this facet
        Args:
            iFacet: index of facet

        Returns:
            grid machine instance
        """

        # GridMachine=ClassDDEGridMachine.ClassDDEGridMachine(SM.GD,#RaDec=self.DicoImager[iFacet]["RaDec"],
        #                                                     SM.DicoImager[iFacet]["DicoConfigGM"]["ChanFreq"],
        #                                                     SM.DicoImager[iFacet]["DicoConfigGM"]["NPix"],
        #                                                     lmShift=SM.DicoImager[iFacet]["lmShift"],
        #                                                     IdSharedMem=IdSharedMem,
        #                                                     IdSharedMemData=IdSharedMemData,
        #                                                     FacetDataCache=FacetDataCache,
        #                                                     ChunkDataCache=ChunkDataCache,
        #                                                     IDFacet=SM.DicoImager[iFacet]["IDFacet"],
        #                                                     SpheNorm=False)
        #                                                     #,
        #                                                     #NFreqBands=self.VS.NFreqBands,
        #                                                     #DataCorrelationFormat=self.VS.StokesConverter.AvailableCorrelationProductsIds(),
        #                                                     #ExpectedOutputStokes=self.VS.StokesConverter.RequiredStokesProductsIds(),
        #                                                     #ListSemaphores=self.ListSemaphores)        

        SpheNorm = False
        FacetInfo = SM.DicoImager[iFacet]
        IDFacet=FacetInfo["IDFacet"]
        cf_dict=shared_dict.attach(SM.Path["cf_dict_path"])[IDFacet]
        #print iFacet,IDFacet
        GridMachine= ClassDDEGridMachine.ClassDDEGridMachine(SM.GD,
                                                             FacetInfo["DicoConfigGM"]["ChanFreq"],
                                                             FacetInfo["DicoConfigGM"]["NPix"],
                                                             lmShift=FacetInfo["lmShift"],
                                                             IDFacet=IDFacet,
                                                             SpheNorm=SpheNorm, 
                                                             NFreqBands=SM.NFreqBands,
                                                             DataCorrelationFormat=SM.AvailableCorrelationProductsIds,
                                                             ExpectedOutputStokes=SM.RequiredStokesProductsIds,
                                                             cf_dict=cf_dict)

        return GridMachine

    def InitGM(self,SM):
        self.DicoGM={}
        LFacets=SM.DicoImager.keys()
        for iFacet in LFacets:
            #print "Initialising Facet %i"%iFacet
            self.DicoGM[iFacet]=self.GiveGM(iFacet,SM)


    def predictKernelPolCluster(self,DicoData,SM,**kwargs):
        if SM.Type=="Catalog":
            return self.predictKernelPolClusterCatalog(DicoData,SM,**kwargs)
        elif SM.Type=="Image":
            return self.predictKernelPolClusterImage(DicoData,SM,**kwargs)



    def predictKernelPolClusterCatalog(self,DicoData,SM,iDirection=None,ApplyJones=None,ApplyTimeJones=None,Noise=None):
        self.DicoData=DicoData
        self.SourceCat=SM.SourceCat
        self.SM=SM

        freq=DicoData["freqs_full"]
        times=DicoData["times"]
        nf=freq.size
        na=DicoData["infos"][0]
        
        nrows=DicoData["A0"].size
        DataOut=np.zeros((nrows,nf,4),self.CType)
        if nrows==0: return DataOut
        
        self.freqs=freq
        self.wave=299792458./self.freqs
        
        if iDirection!=None:
            ListDirection=[iDirection]
        else:
            ListDirection=SM.Dirs#range(SM.NDir)
        
        A0=DicoData["A0"]
        A1=DicoData["A1"]
        if ApplyJones is not None:
            #print "!!!!!",ApplyJones.shape
            #print "!!!!!",ApplyJones.shape
            #print "!!!!!",ApplyJones.shape
            na,NDir,_=ApplyJones.shape
            Jones=np.swapaxes(ApplyJones,0,1)
            Jones=Jones.reshape((NDir,na,4))
            JonesH=ModLinAlg.BatchH(Jones)

        TSmear=0.
        FSmear=0.
        DT=DicoData["infos"][1]
        UVW_dt=DicoData["uvw"]
        if self.DoSmearing:
            if "T" in self.DoSmearing:
                TSmear=1.
                UVW_dt=DicoData["UVW_dt"]
            if "F" in self.DoSmearing:
                FSmear=1.




        # self.SourceCat.m[:]=0
        # self.SourceCat.l[:]=0.1
        # self.SourceCat.I[:]=10
        # self.SourceCat.alpha[:]=0

        # DataOut=DataOut[1:2]
        # self.DicoData["uvw"]=self.DicoData["uvw"][1:2]
        # self.DicoData["A0"]=self.DicoData["A0"][1:2]
        # self.DicoData["A1"]=self.DicoData["A1"][1:2]
        # self.DicoData["IndexTimesThisChunk"]=self.DicoData["IndexTimesThisChunk"][1:2]
        # self.SourceCat=self.SourceCat[0:1]

        
        ColOutDir=np.zeros(DataOut.shape,np.complex64)

        for iCluster in ListDirection:
            ColOutDir.fill(0)

            indSources=np.where(self.SourceCat.Cluster==iCluster)[0]
            T=ClassTimeIt("predict")
            T.disable()
            ### new
            SourceCat=self.SourceCat[indSources].copy()
            #l=np.ones((1,),dtype=np.float64)#,float64(SourceCat.l).copy()
            l=np.require(SourceCat.l, dtype=np.float64, requirements=["A","C"])
            m=np.require(SourceCat.m, dtype=np.float64, requirements=["A","C"])
            
            #m=SourceCat.m#np.float64(SourceCat.m).copy()
            I=np.float32(SourceCat.I)
            Gmaj=np.float32(SourceCat.Gmaj)
            Gmin=np.float32(SourceCat.Gmin)
            GPA=np.float32(SourceCat.Gangle)
            alpha=np.float32(SourceCat.alpha)
            WaveL=np.float64(299792458./self.freqs)
            WaveL=np.require(WaveL, dtype=np.float64, requirements=["A","C"])

            flux=np.float32(SourceCat.I)
            alpha=SourceCat.alpha
            dnu=np.float32(self.DicoData["dfreqs_full"])
            f0=(self.freqs/SourceCat.RefFreq[0])
            fluxFreq=np.float32(flux.reshape((flux.size,1))*(f0.reshape((1,f0.size)))**(alpha.reshape((alpha.size,1))))

            LSM=[l,m,fluxFreq,Gmin,Gmaj,GPA]
            LFreqs=[WaveL,np.float32(self.freqs),dnu]
            LUVWSpeed=[UVW_dt,DT]

            LSmearMode=[FSmear,TSmear]
            T.timeit("init")

            AllowEqualiseChan=IsChanEquidistant(DicoData["freqs_full"])

            if ApplyTimeJones!=None:

                #predict.predictJones(ColOutDir,(DicoData["uvw"]),LFreqs,LSM,LUVWSpeed,LSmearMode,ParamJonesList)

                #ColOutDir.fill(0)
                #predict.predictJones(ColOutDir,(DicoData["uvw"]),LFreqs,LSM,LUVWSpeed,LSmearMode,ParamJonesList,0)
                #d0=ColOutDir.copy()
                #ColOutDir.fill(0)

                # predict.predictJones(ColOutDir,(DicoData["uvw"]),LFreqs,LSM,LUVWSpeed,LSmearMode,ParamJonesList,AllowEqualiseChan)
                # ColOutDir0=ColOutDir.copy()
                # ColOutDir.fill(0)

                #predict.predictJones2(ColOutDir,(DicoData["uvw"]),LFreqs,LSM,LUVWSpeed,LSmearMode,ParamJonesList,AllowEqualiseChan)
                #print LSmearMode
                predict.predictJones2_Gauss(ColOutDir,(DicoData["uvw"]),LFreqs,LSM,LUVWSpeed,LSmearMode,AllowEqualiseChan,
                                            self.LExp,
                                            self.LSinc)
                
                T.timeit("predict")

                ParamJonesList=self.GiveParamJonesList(ApplyTimeJones,A0,A1)
                ParamJonesList=ParamJonesList+[iCluster]

                
                predict.ApplyJones(ColOutDir,ParamJonesList)
                T.timeit("apply")


                # print ColOutDir

                #d1=ColOutDir.copy()
                #ind=np.where(d0!=0)
                #print np.max((d0-d1)[ind]/(d0[ind]))
                #stop
            else:
                #predict.predict(ColOutDir,(DicoData["uvw"]),LFreqs,LSM,LUVWSpeed,LSmearMode)
                #AllowEqualiseChan=0
                #predict.predict(ColOutDir,(DicoData["uvw"]),LFreqs,LSM,LUVWSpeed,LSmearMode,AllowEqualiseChan)
                #d0=ColOutDir.copy()
                #ColOutDir.fill(0)

                
                predict.predictJones2_Gauss(ColOutDir,(DicoData["uvw"]),LFreqs,LSM,LUVWSpeed,LSmearMode,AllowEqualiseChan,
                                            self.LExp,
                                            self.LSinc)
                #print ColOutDir
                #predict.predict(ColOutDir,(DicoData["uvw"]),LFreqs,LSM,LUVWSpeed,LSmearMode,AllowEqualiseChan)
                T.timeit("predict")
                # d0=ColOutDir.copy()
                # ColOutDir.fill(0)

                # predict_np19.predict(ColOutDir,(DicoData["uvw"]),LFreqs,LSM,LUVWSpeed,LSmearMode,AllowEqualiseChan)
                # print ColOutDir,d0
                # d1=ColOutDir.copy()
                # ind=np.where(d0!=0)
                # print np.max((d0-d1)[ind]/(d0[ind]))
                # stop


            del(l,m,I,SourceCat,alpha,WaveL,flux,dnu,f0,fluxFreq,LSM,LFreqs)
            T.timeit("del")


                # d1=ColOutDir
                # ind=np.where(d0!=0)
                # print np.max((d0-d1)[ind]/(d0[ind]))
                # stop
                



            if Noise!=None:
                ColOutDir+=Noise*(np.random.randn(*ColOutDir.shape)+1j*np.random.randn(*ColOutDir.shape))
                stop
            DataOut+=ColOutDir
            T.timeit("add")
        #del(LFreqs,LSM,LUVWSpeed,LSmearMode)
        #del(ColOutDir)

        return DataOut


    ######################################################
    ######################################################
    ######################################################


    def predictKernelPolClusterImage(self,DicoData,SM,iDirection=None,ApplyJones=None,ApplyTimeJones=None,Noise=None,ForceNoDecorr=False):

        T=ClassTimeIt("predictKernelPolClusterImage")
        T.disable()
        self.DicoData=DicoData
        self.SM=SM
        
        
        freq=DicoData["freqs_full"]
        times=DicoData["times"]
        nf=freq.size
        na=DicoData["infos"][0]
        
        nrows=DicoData["A0"].size
        T.timeit("0")
        DataOut=np.zeros((nrows,nf,4),self.CType)
        if nrows==0: return DataOut
        
        self.freqs=freq
        self.wave=299792458./self.freqs
        
        if iDirection!=None:
            ListDirection=[iDirection]
        else:
            ListDirection=SM.Dirs#range(SM.NDir)
        
        A0=DicoData["A0"]
        A1=DicoData["A1"]
        # if ApplyJones!=None:
        #     na,NDir,_=ApplyJones.shape
        #     Jones=np.swapaxes(ApplyJones,0,1)
        #     Jones=Jones.reshape((NDir,na,4))
        #     JonesH=ModLinAlg.BatchH(Jones)

        TSmear=0.
        FSmear=0.


        if self.DoSmearing!=0:
            if "T" in self.DoSmearing:
                TSmear=1.
            if "F" in self.DoSmearing:
                FSmear=1.

        # self.SourceCat.m[:]=0
        # self.SourceCat.l[:]=0.1
        # self.SourceCat.I[:]=10
        # self.SourceCat.alpha[:]=0

        # DataOut=DataOut[1:2]
        # self.DicoData["uvw"]=self.DicoData["uvw"][1:2]
        # self.DicoData["A0"]=self.DicoData["A0"][1:2]
        # self.DicoData["A1"]=self.DicoData["A1"][1:2]
        # self.DicoData["IndexTimesThisChunk"]=self.DicoData["IndexTimesThisChunk"][1:2]
        # self.SourceCat=self.SourceCat[0:1]

        DT=DicoData["infos"][1]
        #UVW_dt=DicoData["UVW_dt"]
        Dnu=DicoData["dfreqs_full"][0]

        ColOutDir=np.zeros(DataOut.shape,np.complex64)
        DATA=DicoData

        T.timeit("1")
        ListFacets=SM.DicoJonesDirToFacet[iDirection]["FacetsIDs"]

        for iFacet in ListFacets:
            if SM.DicoImager[iFacet]["SumFlux"]==0: continue
            GridMachine=self.DicoGM[iFacet]#self.GiveGM(iFacet,SM)
            T.timeit("2: GM")
            uvwThis=DATA["uvw"]
            flagsThis=DATA["flags_image"]
            times=DATA["times"]
            A0=DATA["A0"]
            A1=DATA["A1"]
            A0A1=A0,A1
            freqs=DATA["freqs_full"]


            DicoJonesMatrices=None
            #DicoJonesMatrices=ApplyTimeJones



            # ModelSharedMemName="%sModelImage.Facet_%3.3i"%(self.IdSharedMem,iFacet)
            # print "Facet %i: take model image %s"%(iFacet,ModelSharedMemName)
            # ModelIm = NpShared.GiveArray(ModelSharedMemName)

            #ModelIm = NpShared.UnPackListArray("%sGrids"%self.IdSharedMem)[iFacet]
            ModelIm = SM._model_dict[iFacet]["FacetGrid"]
            
            ChanMapping=np.int32(SM.ChanMappingDegrid)
            # print ChanMapping
            
            GridMachine.LSmear=[]
            DecorrMode = SM.GD["RIME"]["DecorrMode"]
            CondSmear=(not ForceNoDecorr) and (('F' in DecorrMode) | ("T" in DecorrMode))
            if CondSmear:
                #print "DOSMEAR",ForceNoDecorr, (('F' in DecorrMode) | ("T" in DecorrMode))
                uvw_dt = DicoData["UVW_dt"]#DATA["uvw_dt"]
                lm_min=None
                if SM.GD["RIME"]["DecorrLocation"]=="Edge":
                    lm_min=SM.DicoImager[iFacet]["lm_min"]
                GridMachine.setDecorr(uvw_dt, DT, Dnu, SmearMode=SM.GD["RIME"]["DecorrMode"], lm_min=lm_min)
                

            T.timeit("2: Stuff")

            # print """
            # print times.shape,times.dtype
            # print uvwThis.shape,uvwThis.dtype
            # print ColOutDir.shape,ColOutDir.dtype
            # print flagsThis.shape,flagsThis.dtype
            # print ModelIm.shape,ModelIm.dtype"""
            
            # print times.shape,times.dtype
            # print uvwThis.shape,uvwThis.dtype
            # print ColOutDir.shape,ColOutDir.dtype
            # print flagsThis.shape,flagsThis.dtype
            # print ModelIm.shape,ModelIm.dtype
            # stop
            vis=GridMachine.get(times,uvwThis,ColOutDir,flagsThis,A0A1,ModelIm,DicoJonesMatrices=DicoJonesMatrices,freqs=freqs,
                                ImToGrid=False,ChanMapping=ChanMapping)
            T.timeit("2: Predict")
            # get() is substracting
            if ApplyTimeJones is not None and self._BeamAtFacet:
                ParamJonesList=self.GiveParamJonesList(ApplyTimeJones,A0,A1)
                ParamJonesList=ParamJonesList+[iFacet]

                # print "facet"
                # import killMS.Other.rad2hmsdms
                # RA,DEC=SM.DicoImager[iFacet]["RaDec"]
                # sra=killMS.Other.rad2hmsdms.rad2hmsdms(RA,Type="ra")
                # sdec=killMS.Other.rad2hmsdms.rad2hmsdms(DEC,Type="dec")
                # print iFacet,sra,sdec
                
                predict.ApplyJones(ColOutDir,ParamJonesList)


            DataOut-=ColOutDir
            ColOutDir.fill(0)
            
        if ApplyTimeJones is not None and not self._BeamAtFacet:
            #print "tessel"
            #print "apply in direction %i"%iDirection
            ParamJonesList=self.GiveParamJonesList(ApplyTimeJones,A0,A1)
            ParamJonesList=ParamJonesList+[iDirection]
            predict.ApplyJones(DataOut,ParamJonesList)
            T.timeit("apply")
            


        T.timeit("2: End")


        return DataOut
import numpy as np
import time
        import pylab
        import pylab
        import time
                        import pylab
import multiprocessing
from killMS.Array import NpShared
from killMS.Data import ClassVisServer
from killMS.Array import ModLinAlg
from DDFacet.Other import logger
from killMS.Other import ModColor
from DDFacet.Other.progressbar import ProgressBar
from killMS.Other import ClassTimeIt
from killMS.Other import Counter
from .ClassEvolve import ClassModelEvolution
from itertools import product as ItP
from killMS.Wirtinger import ClassSolverLM
from killMS.Wirtinger import ClassSolverEKF
from killMS.Wirtinger import ClassSolPredictMachine
from SkyModel.Sky import ClassSM
            from killMS.Data import ClassBeam
        from Plot import Graph
from killMS.Predict.PredictGaussPoints_NumExpr5 import ClassPredict
        # print(NDone)
        # if NDone>=867:
        #     FileName="CurrentSols.%i.npy"%NDone
        #     log.print( "Save Solutions in file: %s"%FileName)
        #     log.print( "Save Solutions in file: %s"%FileName)
        #     log.print( "Save Solutions in file: %s"%FileName)
        #     Sols=self.GiveSols()
            
        #     np.save(FileName,Sols.G.copy())
        #     if self.SolverType=="KAFCA":
        #         # np.save("P.%s.npy"%self.GD["Solutions"]['OutSolsName'],np.array(self.PListKeep))
        #         FName="P.%s.%i.npz"%(self.GD["Solutions"]['OutSolsName'],NDone)
        #         log.print( "Save PQ in file: %s"%FName)
        #         log.print( "Save PQ in file: %s"%FName)
        #         log.print( "Save PQ in file: %s"%FName)
        #         np.savez(FName,
        #                  ListP=np.array(self.PListKeep),
        #                  ListQ=np.array(self.QListKeep))
        
        # if self.SolverType=="KAFCA":
        #     self.PListKeep.append(self.P.copy())
        #     self.QListKeep.append(self.Q.copy())
        #     np.save("P.%s.npy"%self.GD["Solutions"]['OutSolsName'],np.array(self.PListKeep))
        #     np.savez("P.%s.npz"%self.GD["Solutions"]['OutSolsName'],
        #              ListP=np.array(self.PListKeep),
        #              ListQ=np.array(self.QListKeep))




#======================================
class WorkerAntennaLM(multiprocessing.Process):
    def __init__(self,
                 work_queue,
                 result_queue,SM,PolMode,SolverType,IdSharedMem,ConfigJacobianAntenna=None,
                 GD=None,JonesToVisChanMapping=None,SM_Compress=None):
        multiprocessing.Process.__init__(self)
        self.work_queue = work_queue
        self.result_queue = result_queue
        self.kill_received = False
        self.exit = multiprocessing.Event()
        self.SM=SM

        self.SM_Compress=SM_Compress
        self.PolMode=PolMode
        self.SolverType=SolverType
        self.IdSharedMem=IdSharedMem
        self.ConfigJacobianAntenna=ConfigJacobianAntenna
        self.GD=GD
        self.JonesToVisChanMapping=JonesToVisChanMapping

        self.InitPM()

        #self.DoCalcEvP=DoCalcEvP
        #self.ThisTime=ThisTime
        #self.e,=kwargs["args"]
        

    def InitPM(self):

        x=np.linspace(0.,15,100000)
        Exp=np.float32(np.exp(-x))
        LExp=[Exp,x[1]-x[0]]

        Sinc=np.zeros(x.shape,np.float32)
        Sinc[0]=1.
        Sinc[1::]=np.sin(x[1::])/(x[1::])
        LSinc=[Sinc,x[1]-x[0]]

        
        self.PM=ClassPredict(Precision="S",DoSmearing=self.GD["SkyModel"]["Decorrelation"],
                             IdMemShared=self.IdSharedMem,
                             LExp=LExp,LSinc=LSinc,
                             BeamAtFacet=(self.GD["Beam"]["BeamAt"].lower() == "facet"))

        if self.GD["ImageSkyModel"]["BaseImageName"]!="" and self.GD["SkyModel"]["SkyModelCol"] is None:
            self.PM.InitGM(self.SM)

        self.PM_Compress=None
        if self.SM_Compress:
            self.PM_Compress=ClassPredict(Precision="S",
                                          DoSmearing=self.GD["SkyModel"]["Decorrelation"],
                                          IdMemShared=self.IdSharedMem,
                                          LExp=LExp,LSinc=LSinc)
            
    def shutdown(self):
        self.exit.set()
    def run(self):

        ####################
        # Parallel
        while not self.kill_received:# and not self.work_queue.qsize()==0:
            iAnt,iChanSol,DoCalcEvP,ThisTime,rms,DoEvP,DoFullPredict,SharedDicoDescriptors = self.work_queue.get()

            
        # ####################
        # # Serial
        # while not self.kill_received and not self.work_queue.qsize()==0:
        #     #iAnt,iChanSol,DoCalcEvP,ThisTime,rms,DoEvP,DoFullPredict,SharedDicoDescriptors = self.work_queue.get(True,2)
        #     #iAnt,iChanSol,DoCalcEvP,ThisTime,rms,DoEvP,DoFullPredict,SharedDicoDescriptors = self.work_queue.get()
        #     # print(iAnt,iChanSol,DoCalcEvP,ThisTime,rms,DoEvP,DoFullPredict,SharedDicoDescriptors)
        #     try:
        #         iAnt,iChanSol,DoCalcEvP,ThisTime,rms,DoEvP,DoFullPredict,SharedDicoDescriptors = self.work_queue.get()
        #     except:
        #         break
        #     # #self.e.wait()
        # ####################

            ch0,ch1=self.JonesToVisChanMapping[iChanSol]

            T0=time.time()
            T=ClassTimeIt.ClassTimeIt("  Worker Ant=%2.2i"%iAnt)
            T.disable()
            # if DoCalcEvP:
            #     T.disable()
            # print SharedDicoDescriptors

            if self.SolverType=="CohJones":
                SolverClass=ClassSolverLM.ClassSolverLM
            elif self.SolverType=="KAFCA":

                SolverClass=ClassSolverEKF.ClassSolverEKF