Merge pull request #494 from o-smirnov/oms-deploy

Fixes #463 #466 #472 #360 #475 #477 #479 #480 #484

DDFacet/Array/shared_dict.py

@@ -106,6 +106,9 @@ def __del__(self):
         if self._delete_items:
             self.delete()
 
+    def is_writeable(self):
+        return self._readwrite
+
     def readwrite(self):
         if not self._load:
             raise RuntimeError("SharedDict %s attached without load permissions" % self.path)

DDFacet/Data/ClassFITSBeam.py

@@ -49,7 +49,10 @@
 class ClassFITSBeam (object):
     def __init__ (self, ms, opts):
         self.ms = ms
-        self.filename = opts["FITSFile"]
+        # filename is potentially a list (frequencies will be matched)
+        self.beamsets = opts["FITSFile"]
+        if type(self.beamsets) is not list:
+            self.beamsets = self.beamsets.split(',')
         self.pa_inc = opts["FITSParAngleIncDeg"]
         self.time_inc = opts["DtBeamMin"]
         self.nchan = opts["NBand"]
@@ -79,18 +82,17 @@ def __init__ (self, ms, opts):
             if len(feed) != 2:
                 raise ValueError,"FITSFeed parameter must be two characters (e.g. 'xy')"
             feed = feed.lower()
-            CORRS = [ a+b for a in feed for b in feed ]
-            print>>log,"polarization basis specified by FITSFeed parameter: %s"%" ".join(CORRS)
+            self.corrs = [ a+b for a in feed for b in feed ]
+            print>>log,"polarization basis specified by FITSFeed parameter: %s"%" ".join(self.corrs)
         else:
             # NB: need to check correlation names better. This assumes four correlations in that order!
             if "x" in self.ms.CorrelationNames[0].lower():
-                CORRS = "xx","xy","yx","yy"
+                self.corrs = "xx","xy","yx","yy"
                 print>>log,"polarization basis is linear (MS corrs: %s)"%" ".join(self.ms.CorrelationNames)
             else:
-                CORRS = "rr","rl","lr","ll"
+                self.corrs = "rr","rl","lr","ll"
                 print>>log,"polarization basis is circular (MS corrs: %s)"%" ".join(self.ms.CorrelationNames)
         # Following code is nicked from Cattery/Siamese/OMS/pybeams_fits.py
-        REIM = "re","im";
         REALIMAG = dict(re="real",im="imag");
 
         # get the Cattery: if an explicit path to Cattery set, use this and import Siamese directly
@@ -119,28 +121,43 @@ def make_beam_filename (filename_pattern,corr,reim):
                       realimag=REALIMAG[reim].lower(),REALIMAG=REALIMAG[reim].upper(),
                       RealImag=REALIMAG[reim].title());
 
-        filename_real = [];
-        filename_imag = [];
-        for corr in CORRS:
-            # make FITS images or nulls for real and imaginary part
-            filename_real.append(make_beam_filename(self.filename,corr,'re'))
-            filename_imag.append(make_beam_filename(self.filename,corr,'im'))
-
-        # load beam interpolator
-        self.vbs = []
-        for filename_pair in zip(filename_real,filename_imag):
-            if filename_pair in ClassFITSBeam._vb_cache:
-                vb = ClassFITSBeam._vb_cache[filename_pair]
-                print>>log,"beam patterns %s %s already in memory"%filename_pair
-            else:
-                print>>log,"loading beam patterns %s %s"%filename_pair
-                vb = InterpolatedBeams.LMVoltageBeam(
-                    verbose=opts["FITSVerbosity"],
-                    l_axis=opts["FITSLAxis"], m_axis=opts["FITSMAxis"]
-                )  # verbose, XY must come from options
-                vb.read(*filename_pair)
-                ClassFITSBeam._vb_cache[filename_pair] = vb
-            self.vbs.append(vb)
+        self.vbs = {}
+
+        # now, self.beamsets specifies a list of filename patterns. We need to find the one with the closest
+        # frequency coverage
+
+        for corr in self.corrs:
+            beamlist = []
+            for beamset in self.beamsets:
+                filenames = make_beam_filename(beamset, corr, 're'), make_beam_filename(beamset, corr, 'im')
+                # get interpolator from cache, or create object
+                vb = ClassFITSBeam._vb_cache.get(filenames)
+                if vb is None:
+                    print>> log, "loading beam patterns %s %s" % filenames
+                    ClassFITSBeam._vb_cache[filenames] = vb = InterpolatedBeams.LMVoltageBeam(
+                        verbose=opts["FITSVerbosity"],
+                        l_axis=opts["FITSLAxis"], m_axis=opts["FITSMAxis"]
+                    )  # verbose, XY must come from options
+                    vb.read(*filenames)
+                else:
+                    print>> log, "beam patterns %s %s already in memory" % filenames
+                # find frequency "distance". If beam frequency range completely overlaps MS frequency range,
+                # this is 0, otherwise a positive number
+                distance = max(vb._freqgrid[0] - self.freqs[0], 0) + \
+                           max(self.freqs[-1] - vb._freqgrid[-1], 0)
+                beamlist.append((distance, vb, filenames))
+            # select beams with smallest distance
+            dist0, vb, filenames = sorted(beamlist)[0]
+            if len(beamlist) > 1:
+                if dist0 == 0:
+                    print>> log, "beam patterns %s %s overlap the frequency coverage" % filenames
+                else:
+                    print>> log, "beam patterns %s %s are closest to the frequency coverage (%.1f MHz max separation)" % (
+                                    filenames[0], filenames[1], dist0*1e-6)
+                print>>log,"  MS coverage is %.1f to %.1f GHz, beams are %.1f to %.1f MHz"%(
+                    self.freqs[0]*1e-6, self.freqs[-1]*1e-6, vb._freqgrid[0]*1e-6, vb._freqgrid[-1]*1e-6)
+            self.vbs[corr] = vb
+
 
     _vb_cache = {}
 
@@ -216,7 +233,7 @@ def evaluateBeam (self, t0, ra, dec):
         # print>>log,m
 
         # get interpolated values. Output shape will be [ndir,nfreq]
-        beamjones = [ self.vbs[i].interpolate(l,m,freq=self.freqs,freqaxis=1) for i in range(4) ]
+        beamjones = [ self.vbs[corr].interpolate(l,m,freq=self.freqs,freqaxis=1) for corr in self.corrs ]
 
         # now make output matrix
         jones = numpy.zeros((ndir,self.ms.na,len(self.freqs),2,2),dtype=numpy.complex64)

DDFacet/Data/ClassJones.py

@@ -179,8 +179,8 @@ def MakeSols(self, StrType, DATA, quiet=False):
         if StrType == "Beam":
 
             if self.FacetMachine is not None:
-                if not(self.HasKillMSSols):
-                    print>>log, "  Getting Jones directions from Nones"
+                if not(self.HasKillMSSols) or self.GD["Beam"]["At"] == "facet":
+                    print>>log, "  Getting beam Jones directions from facets"
                     DicoImager = self.FacetMachine.DicoImager
                     NFacets = len(DicoImager)
                     self.ClusterCatBeam = self.FacetMachine.JonesDirCat
@@ -192,7 +192,7 @@ def MakeSols(self, StrType, DATA, quiet=False):
                     DicoClusterDirs["I"] = self.ClusterCatBeam.I
                     DicoClusterDirs["Cluster"] = self.ClusterCatBeam.Cluster
                 else:
-                    print>>log, "  Getting Jones directions from DDE-solutions"
+                    print>>log, "  Getting beam Jones directions from DDE solution tessels"
                     DicoClusterDirs = self.DicoClusterDirs_kMS
                     NDir = DicoClusterDirs["l"].size
                     self.ClusterCatBeam = np.zeros(

DDFacet/Data/ClassMS.py

@@ -51,62 +51,6 @@
 # except:
 #     print>>log, ModColor.Str("Could not import lofar.stationresponse")
 
-def obs_detail(filename,field=0):
-
-    results={}
-    object=""
-
-    try:
-        to = table(filename+ '/OBSERVATION', readonly=True, ack=False)
-    except RuntimeError:
-        to = None
-    try:
-        tf = table(filename+ '/FIELD', readonly=True, ack=False)
-    except RuntimeError:
-        tf = None
-    if tf is not None and to is not None:
-        print >>log, 'Read observing details successfully'
-    else:
-        print >>log, 'Some observing details missing'
-
-    # Stuff relying on an OBSERVATION table:
-    if to is not None:
-        # Time
-        tm = Time(to[0]['TIME_RANGE']/86400.0,
-                  scale="utc",
-                  format='mjd')
-        results['DATE-OBS'] = tm[0].iso.split()[0]
-
-        # Object
-        try:
-            object = to[0]['LOFAR_TARGET'][0]
-        except:
-            pass
-
-        # Telescope
-        telescope=to[0]['TELESCOPE_NAME']
-        results['TELESCOP'] = telescope
-
-        # observer
-        observer = to[0]['OBSERVER']
-        results['OBSERVER'] = observer
-
-    if not object and tf is not None:
-        object = tf[field]['NAME'] 
-
-    if object:
-        results['OBJECT'] = object
-
-    # Time now
-    tn = Time(time.time(),format='unix')
-    results['DATE-MAP'] = tn.iso.split()[0]
-
-    if to is not None:
-        to.close()
-    if tf is not None:
-        tf.close()
-
-    return results
 
 class ClassMS():
     def __init__(self,MSname,Col="DATA",zero_flag=True,ReOrder=False,EqualizeFlag=False,DoPrint=True,DoReadData=True,
@@ -145,8 +89,7 @@ def __init__(self,MSname,Col="DATA",zero_flag=True,ReOrder=False,EqualizeFlag=Fa
         if self.ToRADEC is "": self.ToRADEC=None
 
         self.AverageSteps=AverageTimeFreq
-        MSname= reformat.reformat(os.path.abspath(MSname), LastSlash=False)
-        self.MSName=MSname
+        self.MSName = MSName = reformat.reformat(os.path.abspath(MSname), LastSlash=False)
         self.ColName=Col
         self.ChanSlice = ChanSlice or slice(None)
         self.zero_flag=zero_flag
@@ -193,8 +136,63 @@ def __init__(self,MSname,Col="DATA",zero_flag=True,ReOrder=False,EqualizeFlag=Fa
         if GetBeam:
             self.LoadSR()
 
-        if get_obs_detail:
-            self.obs_detail=obs_detail(self.MSName, field=self.Field)
+    def get_obs_details(self):
+        """Gets observer details from MS, for FITS header mainly"""
+
+        results = {}
+        object = ""
+
+        try:
+            to = table(self.MSName + '/OBSERVATION', readonly=True, ack=False)
+        except RuntimeError:
+            to = None
+        try:
+            tf = table(self.MSName + '/FIELD', readonly=True, ack=False)
+        except RuntimeError:
+            tf = None
+        if tf is not None and to is not None:
+            print >> log, 'Read observing details from %s'%self.MSName
+        else:
+            print >> log, 'Some observing details in %s missing'%self.MSName
+
+        # Stuff relying on an OBSERVATION table:
+        if to is not None:
+            # Time
+            tm = Time(to[0]['TIME_RANGE'] / 86400.0,
+                      scale="utc",
+                      format='mjd')
+            results['DATE-OBS'] = tm[0].iso.split()[0]
+
+            # Object
+            try:
+                object = to[0]['LOFAR_TARGET'][0]
+            except:
+                pass
+
+            # Telescope
+            telescope = to[0]['TELESCOPE_NAME']
+            results['TELESCOP'] = telescope
+
+            # observer
+            observer = to[0]['OBSERVER']
+            results['OBSERVER'] = observer
+
+        if not object and tf is not None:
+            object = tf[self.Field]['NAME']
+
+        if object:
+            results['OBJECT'] = object
+
+        # Time now
+        tn = Time(time.time(), format='unix')
+        results['DATE-MAP'] = tn.iso.split()[0]
+
+        if to is not None:
+            to.close()
+        if tf is not None:
+            tf.close()
+
+        return results
 
     def GiveMainTable (self,**kw):
         """Returns main MS table, applying TaQL selection if any"""
@@ -902,11 +900,28 @@ def RemoveStation(self):
         self.nbl=(na*(na-1))/2+na
 
 
+    # static member caching DDID/FIELD_ID lookups
+    _ddid_field_cache = {}
+
     def ReadMSInfo(self,DoPrint=True):
         T= ClassTimeIt.ClassTimeIt()
         T.enableIncr()
         T.disable()
 
+        # quick check if DDID and FieldId is present at all. This is much faster than running a full query
+        # (which GiveMainTable() does), helps when many MSs are specified, many of them missing DDIDs
+        if self.MSName in ClassMS._ddid_field_cache:
+            ddid_fields = ClassMS._ddid_field_cache.get(self.MSName)
+        else:
+            maintab = table(self.MSName, ack=False)
+            ddid_fields = set(zip(maintab.getcol("FIELD_ID"), maintab.getcol("DATA_DESC_ID")))
+            ClassMS._ddid_field_cache[self.MSName] = ddid_fields
+
+        self.empty = (self.Field,self.DDID) not in ddid_fields
+        if self.empty:
+            print>>log, ModColor.Str("MS %s (field %d, ddid %d): no rows, skipping"%(self.MSName, self.Field, self.DDID))
+            return
+
         # open main table
         table_all = self.GiveMainTable()
         self.empty = not table_all.nrows()

DDFacet/Data/ClassSmearMapping.py

@@ -584,6 +584,7 @@ def GiveBlocksRowsListBL_old(a0, a1, DATA, dPhi, l, channel_mapping):
                 blocklist.append(ThiDesc)
                 BlocksRowsListBL += (ThiDesc)
                 NBlocksTotBL += 1
+            # import pdb; pdb.set_trace()
             NChanBlockMax = 1e3
             CurrentRows = []
             duvtot = 0