Changeset 907

Timestamp:

03/20/06 15:17:06 (19 years ago)

Author:

vor010

Message:

LineFinder & auto_poly_baseline: a support of
new scantable format has been added. Now findLines accept a mask and edge
parameters, which can therefore be different for different rows. Constructor
need now just a scan table. Boost types removed from STLineFinder. auto_poly_baseline can accept a nested tuple of edges, which would be interpreted as
different edge parameters for different IFs

Location:

trunk

Files:

• python/asaplinefind.py (modified) (3 diffs)
• python/scantable.py (modified) (4 diffs)
• src/Makefile (modified) (2 diffs)
• src/STLineFinder.cpp (modified) (10 diffs)
• src/STLineFinder.h (modified) (1 diff)

trunk/python/asaplinefind.py

@@ -7 +7 @@
     Example:
        fl=linefinder()
-       fl.set_scan(sc,edge=(50,))
+       fl.set_scan(sc)
        fl.set_options(threshold=3)
-       nlines=fl.find_lines()
+       nlines=fl.find_lines(edge=(50,0))
        if nlines!=0:
           print "Found ",nlines," spectral lines"
@@ -65 +65 @@
         return
 
-    def set_scan(self, scan, mask=None, edge=(0,0)):
+    def set_scan(self, scan):
         """
         Set the 'data' (scantable) to work with.
         Parameters:
              scan:    a scantable
-             mask:       an optional mask retreived from scantable
-             edge:       an optional number of channel to drop at
-                         the edge of spectrum. If only one value is
+        """
+        if not scan:
+           raise RuntimeError, 'Please give a correct scan'
+        return
+    def find_lines(self,nRow=0,mask=None,edge=(0,0)):
+        """
+        Search for spectral lines in the scan assigned in set_scan.
+        Parameters:
+             nRow:       a row in the scantable to work with
+             mask:       an optional mask (e.g. retreived from scantable)
+             edge:       an optional number of channels to drop at
+                         the edge of the spectrum. If only one value is
                          specified, the same number will be dropped from
                          both sides of the spectrum. Default is to keep
                          all channels
+        A number of lines found will be returned
         """
-        if not scan:
-           raise RuntimeError, 'Please give a correct scan'
-        if not scan._check_ifs():
-           raise RuntimeError, 'IFs with different numbers of channels are not yet supported'
-
         if isinstance(edge,int):
            edge=(edge,)
@@ -96 +101 @@
         if mask is None:
             from numarray import ones
-            self.finder.setscan(scan,ones(scan.nchan(-1)),tuple(edge))
+            return self.finder.findlines(ones(scan.nchan(nRow)),list(edge),nRow)
         else:
-            self.finder.setscan(scan,mask,tuple(edge))
-        return
-    def find_lines(self,nRow=0):
-        """
-        Search for spectral lines in the scan assigned in set_scan.
-        Current Beam/IF/Pol is used, Row is specified by parameter
-        A number of lines found will be returned
-        """
-        return self.finder.findlines(nRow)
+            return self.finder.setscan(mask,list(edge),nRow)
     def get_mask(self,invert=False):
         """

trunk/python/scantable.py

@@ -1033 +1033 @@
                         specified, the same number will be dropped from
                         both sides of the spectrum. Default is to keep
-                        all channels
+                        all channels. Nested tuples represent individual
+                        edge selection for different IFs (a number of spectral
+                        channels can be different)
             order:      the order of the polynomial (default is 0)
             threshold:  the threshold used by line finder. It is better to
@@ -1051 +1053 @@
         from asap import _is_sequence_or_number as _is_valid
 
-        if not _is_valid(edge, int):
+        # check whether edge is set up for each IF individually
+        individualEdge = False;
+        if len(edge)>1:
+           if isinstance(edge[0],list) or isinstance(edge[0],tuple):
+           individualEdge = True;
+
+        if not _is_valid(edge, int) and not individualEdge:
             raise RuntimeError, "Parameter 'edge' has to be an integer or a \
-            pair of integers specified as a tuple"
+            pair of integers specified as a tuple. Nested tuples are allowed \
+            to make individual selection for different IFs."
+        
+        curedge = (0,0)
+        if individualEdge:
+           for edge_par in edge:
+               if not _is_valid(edge,int):
+                  raise RuntimeError, "Each element of the 'edge' tuple has \
+                  to be a pair of integers or an integer."
+        else:
+           curedge = edge;
 
         # setup fitter
@@ -1067 +1085 @@
         else:
             workscan=self
+
+        fl.set_scan(workscan)
 
         rows=range(workscan.nrow())
@@ -1074 +1094 @@
             msg = " Scan[%d] Beam[%d] IF[%d] Pol[%d] Cycle[%d]" %        (workscan.getscan(r),workscan.getbeam(r),workscan.getif(r),workscan.getpol(r), workscan.getcycle(r))
             asaplog.push(msg, False)
-            fl.set_scan(workscan, mask, edge)
-            fl.find_lines(r)
+
+            # figure out edge parameter
+            if individualEdge:
+               if len(edge)>=workscan.getif(r):
+                  raise RuntimeError, "Number of edge elements appear to be less than the number of IFs"
+                  curedge = edge[workscan.getif(r)]
+            
+            # setup line finder
+            fl.find_lines(r,mask,curedge)
             f.set_scan(workscan, fl.get_mask())
             f.x = workscan._getabcissa(r)

trunk/src/Makefile

@@ -1 @@
-TARGET   := /tmp/asap/_asap.so
+TARGET   := ../tmp/asap/_asap.so
 
 # the casa environment AIPSPATH has to be defined
@@ -64 +64 @@
 # python only 2.3 has been tested
 PYVERSION := 2.3
-PYTHONROOT := /usr
+#PYTHONROOT := /usr
+PYTHONROOT := /export/DEBIANlocal
 PYTHONINC := -I$(PYTHONROOT)/include/python$(PYVERSION)
-PYTHONLIB := -L$(PYTHONROOT)/lib -lpython$(PYVERSION)
+PYTHONLIB := -L$(PYTHONROOT)/lib -lpython$(PYVERSION) -L/export/DEBIANlocal/gnu/lib/
 
 # has to be build with same g++ version as casa

trunk/src/STLineFinder.cpp

@@ -150 +150 @@
                                            // the value - current mean
                                            // (used to search wings)
+   casa::Int last_sign;                    // a sign (+1, -1 or 0) of the
+                                           // last point of the detected line
+                                           //
 public:
 
@@ -189 +192 @@
    void processCurLine(const casa::Vector<casa::Bool> &mask)
                                                  throw(casa::AipsError);
+   
+   // get the sign of runningBox->aboveMean(). The RunningBox pointer
+   // should be defined
+   casa::Int getAboveMeanSign() const throw();
 };
 
@@ -370 +377 @@
 }
 
+// get the sign of runningBox->aboveMean(). The RunningBox pointer
+// should be defined
+casa::Int LFAboveThreshold::getAboveMeanSign() const throw()
+{
+  const Float buf=running_box->aboveMean();
+  if (buf>0) return 1;
+  if (buf<0) return -1;
+  return 0;
+}
+
 
 // process a channel: update cur_line and is_detected before and
@@ -377 +394 @@
 {
   try {
-       if (detect) {
-           if (is_detected_before)
-               cur_line.second=running_box->getChannel()+1;
-           else {
-               is_detected_before=True;
-               cur_line.first=running_box->getChannel();
-               cur_line.second=running_box->getChannel()+1;
-           }
-       } else processCurLine(mask);
+       if (is_detected_before) {
+          // we have to check that the current detection has the
+          // same sign of running_box->aboveMean
+          // otherwise it could be a spurious detection
+          if (last_sign && last_sign!=getAboveMeanSign())
+              detect=False;
+          }
+          if (detect) {
+              last_sign=getAboveMeanSign();
+              if (is_detected_before)
+                  cur_line.second=running_box->getChannel()+1;
+              else {
+                  is_detected_before=True;
+                  cur_line.first=running_box->getChannel();
+                  cur_line.second=running_box->getChannel()+1;
+              }
+          } else processCurLine(mask);
   }
   catch (const AipsError &ae) {
@@ -490 +515 @@
                   threshold*offline_variance), mask);
            else processCurLine(mask); // just finish what was accumulated before
+
+           signs[ch]=getAboveMeanSign();
+           // os<<ch<<" "<<spectrum[ch]<<" "<<fabs(running_box->aboveMean())<<" "<<
+           //   threshold*offline_variance<<endl;
+           
            const Float buf=running_box->aboveMean();
            if (buf>0) signs[ch]=1;
@@ -635 +665 @@
 STLineFinder::~STLineFinder() throw(AipsError) {}
 
-// set scan to work with (in_scan parameter), associated mask (in_mask
-// parameter) and the edge channel rejection (in_edge parameter)
+// set scan to work with (in_scan parameter)
+void STLineFinder::setScan(const ScantableWrapper &in_scan) throw(AipsError)
+{
+  scan=in_scan.getCP();
+  AlwaysAssert(!scan.null(),AipsError);
+ 
+}
+
+// search for spectral lines. Number of lines found is returned
+// in_edge and in_mask control channel rejection for a given row
 //   if in_edge has zero length, all channels chosen by mask will be used
 //   if in_edge has one element only, it represents the number of
@@ -642 +680 @@
 //   in_edge is introduced for convinience, although all functionality
 //   can be achieved using a spectrum mask only
-void STLineFinder::setScan(const ScantableWrapper &in_scan,
-               const std::vector<bool> &in_mask,
-               const boost::python::tuple &in_edge) throw(AipsError)
-{
-  try {
-       scan=in_scan.getCP();
-       AlwaysAssert(!scan.null(),AipsError);
-
-       mask=in_mask;
-       if (mask.nelements()!=scan->nchan())
-           throw AipsError("STLineFinder::setScan - in_scan and in_mask have different"
-                           "number of spectral channels.");
-
-       // number of elements in the in_edge tuple
-       int n=extract<int>(in_edge.attr("__len__")());
-       if (n>2 || n<0)
-           throw AipsError("STLineFinder::setScan - the length of the in_edge parameter"
-                           "should not exceed 2");
-       if (!n) {
-           // all spectra, no rejection
-           edge.first=0;
-           edge.second=scan->nchan();
-       } else {
-           edge.first=extract<int>(in_edge.attr("__getitem__")(0));
-           if (edge.first<0)
-               throw AipsError("STLineFinder::setScan - the in_edge parameter has a negative"
-                                "number of channels to drop");
-           if (edge.first>=scan->nchan())
-               throw AipsError("STLineFinder::setScan - all channels are rejected by the in_edge parameter");
-           if (n==2) {
-               edge.second=extract<int>(in_edge.attr("__getitem__")(1));
-               if (edge.second<0)
-                   throw AipsError("STLineFinder::setScan - the in_edge parameter has a negative"
-                                   "number of channels to drop");
-               edge.second=scan->nchan()-edge.second;
-           } else edge.second=scan->nchan()-edge.first;
-           if (edge.second<0 || (edge.first>=edge.second))
-               throw AipsError("STLineFinder::setScan - all channels are rejected by the in_edge parameter");
-       }
-  }
-  catch(const AipsError &ae) {
-       // setScan is unsuccessfull, reset scan/mask/edge
-       scan=CountedConstPtr<Scantable>(); // null pointer
-       mask.resize(0);
-       edge=pair<int,int>(0,0);
-       throw;
-  }
-}
-
-// search for spectral lines. Number of lines found is returned
-int STLineFinder::findLines(const casa::uInt &whichRow) throw(casa::AipsError)
+int STLineFinder::findLines(const std::vector<bool> &in_mask,
+                const std::vector<int> &in_edge,
+                const casa::uInt &whichRow) throw(casa::AipsError)
 {
   const int minboxnchan=4;
@@ -698 +688 @@
       throw AipsError("STLineFinder::findLines - a scan should be set first,"
                       " use set_scan");
-  DebugAssert(mask.nelements()==scan->nchan(), AipsError);
-  int max_box_nchan=int(scan->nchan()*box_size); // number of channels in running
+  
+  // set up mask and edge rejection
+  mask=in_mask;
+  if (mask.nelements()!=scan->nchan(whichRow))
+      throw AipsError("STLineFinder::findLines - in_scan and in_mask have different"
+            "number of spectral channels.");
+
+  // number of elements in in_edge
+  if (in_edge.size()>2)
+      throw AipsError("STLineFinder::findLines - the length of the in_edge parameter"
+                      "should not exceed 2");
+      if (!in_edge.size()) {
+           // all spectra, no rejection
+           edge.first=0;
+           edge.second=scan->nchan(whichRow);
+      } else {
+           edge.first=in_edge[0];
+           if (edge.first<0)
+               throw AipsError("STLineFinder::findLines - the in_edge parameter has a negative"
+                                "number of channels to drop");
+           if (edge.first>=scan->nchan(whichRow))
+               throw AipsError("STLineFinder::findLines - all channels are rejected by the in_edge parameter");
+           if (in_edge.size()==2) {
+               edge.second=in_edge[1];
+               if (edge.second<0)
+                   throw AipsError("STLineFinder::findLines - the in_edge parameter has a negative"
+                                   "number of channels to drop");
+               edge.second=scan->nchan(whichRow)-edge.second;
+           } else edge.second=scan->nchan(whichRow)-edge.first;
+           if (edge.second<0 || (edge.first>=edge.second))
+               throw AipsError("STLineFinder::findLines - all channels are rejected by the in_edge parameter");
+       }
+     
+  //
+  int max_box_nchan=int(scan->nchan(whichRow)*box_size); // number of channels in running
                                                  // box
   if (max_box_nchan<2)
@@ -833 +856 @@
            throw AipsError("STLineFinder::getMask - a scan should be set first,"
                       " use set_scan followed by find_lines");
-       DebugAssert(mask.nelements()==scan->nchan(), AipsError);
+       DebugAssert(mask.nelements()==scan->nchan(last_row_used), AipsError);
        /*
        if (!lines.size())
@@ -892 +915 @@
            throw AipsError("STLineFinder::getLineRangesInChannels - a scan should be set first,"
                       " use set_scan followed by find_lines");
-       DebugAssert(mask.nelements()==scan->nchan(), AipsError);
+       DebugAssert(mask.nelements()==scan->nchan(last_row_used), AipsError);
 
        if (!lines.size())

trunk/src/STLineFinder.h

@@ -158 +158 @@
                    const casa::Float &in_box_size=0.2) throw();
 
-   // set the scan to work with (in_scan parameter), associated mask (in_mask
-   // parameter) and the edge channel rejection (in_edge parameter)
-   //   if in_edge has zero length, all channels chosen by mask will be used
+   // set the scan to work with (in_scan parameter)
+   void setScan(const ScantableWrapper &in_scan) throw(casa::AipsError);
+
+   // search for spectral lines in a row specified by whichRow 
+   // in_mask and in_edge parameters control channel rejection
+   //  if in_edge has zero length, all channels chosen by mask will be used
    //   if in_edge has one element only, it represents the number of
    //      channels to drop from both sides of the spectrum
    //   in_edge is introduced for convinience, although all functionality
    //   can be achieved using a spectrum mask only
-   void setScan(const ScantableWrapper &in_scan,
-                const std::vector<bool> &in_mask,
-                const boost::python::tuple &in_edge) throw(casa::AipsError);
-
-   // search for spectral lines for a row specified by whichRow and
-   // Beam/IF/Pol specified by current cursor set for the scantable
    // Number of lines found is returned
-   int findLines(const casa::uInt &whichRow = 0) throw(casa::AipsError);
+   int findLines(const std::vector<bool> &in_mask,
+                 const std::vector<int> &in_edge = std::vector<int>(),
+                 const casa::uInt &whichRow = 0) throw(casa::AipsError);
 
    // get the mask to mask out all lines that have been found (default)