Changeset 1907

Timestamp:

08/26/10 19:46:26 (14 years ago)

Author:

WataruKawasaki

Message:

New Development: No

JIRA Issue: Yes CAS-1937,CAS-2373

Ready to Release: Yes

Interface Changes: Yes

What Interface Changed: A new parameter 'batch' was added to

sd.scantable.poly_baseline(), while
'uselin' was removed.

Test Programs:

Put in Release Notes: Yes

Module(s): sdbaseline

Description: A faster version of sd.scantable.poly_baseline().

Location:

trunk

Files:

• python/scantable.py (modified) (20 diffs)
• src/STMath.cpp (modified) (1 diff)
• src/STMath.h (modified) (1 diff)
• src/STMathWrapper.h (modified) (1 diff)
• src/Scantable.cpp (modified) (1 diff)
• src/Scantable.h (modified) (4 diffs)
• src/ScantableWrapper.h (modified) (1 diff)
• src/python_STMath.cpp (modified) (1 diff)
• src/python_Scantable.cpp (modified) (1 diff)

trunk/python/scantable.py

@@ -16 +16 @@
 from asap.coordinate import coordinate
 from asap.utils import _n_bools, mask_not, mask_and, mask_or, page
+from asap.asapfitter import fitter
 
 
@@ -486 +487 @@
 
     @asaplog_post_dec
-    def stats(self, stat='stddev', mask=None, form='3.3f'):
+    def stats(self, stat='stddev', mask=None, form='3.3f', row=None):
         """\
         Determine the specified statistic of the current beam/if/pol
@@ -502 +503 @@
             form:    format string to print statistic values
 
-        Example::
-
+            row:     row number of spectrum to process.
+                     (default is None: for all rows)
+
+        Example:
             scan.set_unit('channel')
             msk = scan.create_mask([100, 200], [500, 600])
@@ -521 +524 @@
             getchan = True
             statvals = []
-        if not rtnabc: statvals = self._math._stats(self, mask, stat)
+        if not rtnabc:
+            if row == None:
+                statvals = self._math._stats(self, mask, stat)
+            else:
+                statvals = self._math._statsrow(self, mask, stat, int(row))
 
         #def cb(i):
@@ -533 +540 @@
         #outvec = []
         sep = '-'*50
-        for i in range(self.nrow()):
+
+        if row == None:
+            rows = xrange(self.nrow())
+        elif isinstance(row, int):
+            rows = [ row ]
+
+        for i in rows:
             refstr = ''
             statunit= ''
@@ -549 +562 @@
             out += 'Scan[%d] (%s) ' % (self.getscan(i), src)
             out += 'Time[%s]:\n' % (tm)
-            if self.nbeam(-1) > 1:
-                out +=  ' Beam[%d] ' % (self.getbeam(i))
-            if self.nif(-1) > 1: out +=  ' IF[%d] ' % (self.getif(i))
-            if self.npol(-1) > 1: out +=  ' Pol[%d] ' % (self.getpol(i))
+            if self.nbeam(-1) > 1: out +=  ' Beam[%d] ' % (self.getbeam(i))
+            if self.nif(-1) > 1:   out +=  ' IF[%d] ' % (self.getif(i))
+            if self.npol(-1) > 1:  out +=  ' Pol[%d] ' % (self.getpol(i))
             #outvec.append(callback(i))
-            #out += ('= %'+form) % (outvec[i]) +'   '+refstr+'\n'
-            out += ('= %'+form) % (statvals[i]) +'   '+refstr+'\n'
+            if len(rows) > 1:
+                # out += ('= %'+form) % (outvec[i]) +'   '+refstr+'\n'
+                out += ('= %'+form) % (statvals[i]) +'   '+refstr+'\n'
+            else:
+                # out += ('= %'+form) % (outvec[0]) +'   '+refstr+'\n'
+                out += ('= %'+form) % (statvals[0]) +'   '+refstr+'\n'
             out +=  sep+"\n"
-
 
         import os
@@ -1837 +1852 @@
 
     @asaplog_post_dec
-    def poly_baseline(self, mask=None, order=0, plot=False, uselin=False,
-                      insitu=None):
+    def old_poly_baseline(self, mask=None, order=0, plot=False, uselin=False, insitu=None, rows=None):
         """\
         Return a scan which has been baselined (all rows) by a polynomial.
-
+        
         Parameters:
 
@@ -1858 +1872 @@
                         The default is taken from .asaprc (False)
 
-        Example::
-
+            rows:       row numbers of spectra to be processed.
+                        (default is None: for all rows)
+        
+        Example:
             # return a scan baselined by a third order polynomial,
             # not using a mask
@@ -1872 +1888 @@
         varlist = vars()
         if mask is None:
-            mask = [True for i in xrange(self.nchan(-1))]
-
-        from asap.asapfitter import fitter
+            mask = [True for i in xrange(self.nchan())]
+
         try:
             f = fitter()
@@ -1882 +1897 @@
                 f.set_function(poly=order)
 
-            rows = range(workscan.nrow())
+            if rows == None:
+                rows = xrange(workscan.nrow())
+            elif isinstance(rows, int):
+                rows = [ rows ]
+            
             if len(rows) > 0:
                 self.blpars = []
-
+                self.masklists = []
+                self.actualmask = []
+            
             for r in rows:
-                # take into account flagtra info (CAS-1434)
-                flagtra = workscan._getmask(r)
-                actualmask = mask[:]
-                if len(actualmask) == 0:
-                    actualmask = list(flagtra[:])
-                else:
-                    if len(actualmask) != len(flagtra):
-                        raise RuntimeError, "Mask and flagtra have different length"
-                    else:
-                        for i in range(0, len(actualmask)):
-                            actualmask[i] = actualmask[i] and flagtra[i]
-                f.set_scan(workscan, actualmask)
                 f.x = workscan._getabcissa(r)
                 f.y = workscan._getspectrum(r)
+                f.mask = mask_and(mask, workscan._getmask(r))    # (CAS-1434)
                 f.data = None
                 f.fit()
@@ -1908 +1918 @@
                     if x.upper() == 'N':
                         self.blpars.append(None)
+                        self.masklists.append(None)
+                        self.actualmask.append(None)
                         continue
                 workscan._setspectrum(f.fitter.getresidual(), r)
                 self.blpars.append(f.get_parameters())
+                self.masklists.append(workscan.get_masklist(f.mask, row=r))
+                self.actualmask.append(f.mask)
 
             if plot:
@@ -1925 +1939 @@
 
 
-    def auto_poly_baseline(self, mask=[], edge=(0, 0), order=0,
+    def poly_baseline(self, mask=None, order=0, plot=False, batch=False, insitu=None, rows=None):
+        """\
+        Return a scan which has been baselined (all rows) by a polynomial.
+        Parameters:
+            mask:       an optional mask
+            order:      the order of the polynomial (default is 0)
+            plot:       plot the fit and the residual. In this each
+                        indivual fit has to be approved, by typing 'y'
+                        or 'n'. Ignored if batch = True. 
+            batch:      if True a faster algorithm is used and logs
+                        including the fit results are not output
+                        (default is False) 
+            insitu:     if False a new scantable is returned.
+                        Otherwise, the scaling is done in-situ
+                        The default is taken from .asaprc (False)
+            rows:       row numbers of spectra to be processed.
+                        (default is None: for all rows)
+        Example:
+            # return a scan baselined by a third order polynomial,
+            # not using a mask
+            bscan = scan.poly_baseline(order=3)
+        """
+        if insitu is None: insitu = rcParams["insitu"]
+        if insitu:
+            workscan = self
+        else:
+            workscan = self.copy()
+
+        varlist = vars()
+        nchan = workscan.nchan()
+        
+        if mask is None:
+            mask = [True for i in xrange(nchan)]
+
+        try:
+            if rows == None:
+                rows = xrange(workscan.nrow())
+            elif isinstance(rows, int):
+                rows = [ rows ]
+            
+            if len(rows) > 0:
+                self.blpars = []
+                self.masklists = []
+                self.actualmask = []
+
+            if batch:
+                for r in rows:
+                    workscan._poly_baseline(mask, order, r)
+            elif plot:
+                f = fitter()
+                f.set_function(lpoly=order)
+                for r in rows:
+                    f.x = workscan._getabcissa(r)
+                    f.y = workscan._getspectrum(r)
+                    f.mask = mask_and(mask, workscan._getmask(r))    # (CAS-1434)
+                    f.data = None
+                    f.fit()
+                    
+                    f.plot(residual=True)
+                    accept_fit = raw_input("Accept fit ( [y]/n ): ")
+                    if accept_fit.upper() == "N":
+                        self.blpars.append(None)
+                        self.masklists.append(None)
+                        self.actualmask.append(None)
+                        continue
+                    workscan._setspectrum(f.fitter.getresidual(), r)
+                    self.blpars.append(f.get_parameters())
+                    self.masklists.append(workscan.get_masklist(f.mask, row=r))
+                    self.actualmask.append(f.mask)
+                    
+                f._p.unmap()
+                f._p = None
+            else:
+                import array
+                for r in rows:
+                    pars = array.array("f", [0.0 for i in range(order+1)])
+                    pars_adr = pars.buffer_info()[0]
+                    pars_len = pars.buffer_info()[1]
+                    
+                    errs = array.array("f", [0.0 for i in range(order+1)])
+                    errs_adr = errs.buffer_info()[0]
+                    errs_len = errs.buffer_info()[1]
+                    
+                    fmsk = array.array("i", [1 for i in range(nchan)])
+                    fmsk_adr = fmsk.buffer_info()[0]
+                    fmsk_len = fmsk.buffer_info()[1]
+                    
+                    workscan._poly_baseline(mask, order, r, pars_adr, pars_len, errs_adr, errs_len, fmsk_adr, fmsk_len)
+                    
+                    params = pars.tolist()
+                    fmtd = ""
+                    for i in xrange(len(params)): fmtd += "  p%d= %3.6f," % (i, params[i])
+                    fmtd = fmtd[:-1]  # remove trailing ","
+                    errors = errs.tolist()
+                    fmask = fmsk.tolist()
+                    for i in xrange(len(fmask)): fmask[i] = (fmask[i] > 0)    # transform (1/0) -> (True/False)
+                    
+                    self.blpars.append({"params":params, "fixed":[], "formatted":fmtd, "errors":errors})
+                    self.masklists.append(workscan.get_masklist(fmask, r))
+                    self.actualmask.append(fmask)
+                    
+                    asaplog.push(str(fmtd))
+            
+            workscan._add_history("poly_baseline", varlist)
+            print_log()
+            
+            if insitu:
+                self._assign(workscan)
+            else:
+                return workscan
+            
+        except RuntimeError:
+            msg = "The fit failed, possibly because it didn't converge."
+            if rcParams["verbose"]:
+                print_log()
+                asaplog.push(str(msg))
+                print_log("ERROR")
+                return
+            else:
+                raise RuntimeError(msg)
+
+
+    def auto_poly_baseline(self, mask=None, edge=(0, 0), order=0,
                            threshold=3, chan_avg_limit=1, plot=False,
-                           insitu=None):
+                           insitu=None, rows=None):
         """\
         Return a scan which has been baselined (all rows) by a polynomial.
@@ -1968 +2104 @@
                         Otherwise, the scaling is done in-situ
                         The default is taken from .asaprc (False)
+            rows:       row numbers of spectra to be processed.
+                        (default is None: for all rows)
 
 
@@ -1977 +2115 @@
         if insitu is None: insitu = rcParams['insitu']
         varlist = vars()
-        from asap.asapfitter import fitter
         from asap.asaplinefind import linefinder
         from asap import _is_sequence_or_number as _is_valid
@@ -2001 +2138 @@
             curedge = edge;
 
+        if not insitu:
+            workscan = self.copy()
+        else:
+            workscan = self
+
         # setup fitter
         f = fitter()
-        f.set_function(poly=order)
+        f.set_function(lpoly=order)
 
         # setup line finder
@@ -2009 +2151 @@
         fl.set_options(threshold=threshold,avg_limit=chan_avg_limit)
 
-        if not insitu:
-            workscan = self.copy()
-        else:
-            workscan = self
-
         fl.set_scan(workscan)
 
-        rows = range(workscan.nrow())
+        if mask is None:
+            mask = _n_bools(workscan.nchan(), True)
+        
+        if rows is None:
+            rows = xrange(workscan.nrow())
+        elif isinstance(rows, int):
+            rows = [ rows ]
+        
         # Save parameters of baseline fits & masklists as a class attribute.
         # NOTICE: It does not reflect changes in scantable!
@@ -2022 +2166 @@
             self.blpars=[]
             self.masklists=[]
+            self.actualmask=[]
         asaplog.push("Processing:")
         for r in rows:
@@ -2036 +2181 @@
                     curedge = edge[workscan.getif(r)]
 
-            # take into account flagtra info (CAS-1434)
-            flagtra = workscan._getmask(r)
-            actualmask = mask[:]
-            if len(actualmask) == 0:
-                actualmask = list(flagtra[:])
-            else:
-                if len(actualmask) != len(flagtra):
-                    raise RuntimeError, "Mask and flagtra have different length"
-                else:
-                    for i in range(0, len(actualmask)):
-                        actualmask[i] = actualmask[i] and flagtra[i]
+            actualmask = mask_and(mask, workscan._getmask(r))    # (CAS-1434)
 
             # setup line finder
             fl.find_lines(r, actualmask, curedge)
-            outmask=fl.get_mask()
-            f.set_scan(workscan, fl.get_mask())
+            
             f.x = workscan._getabcissa(r)
             f.y = workscan._getspectrum(r)
+            f.mask = fl.get_mask()
             f.data = None
             f.fit()
 
             # Show mask list
-            masklist=workscan.get_masklist(fl.get_mask(),row=r)
+            masklist=workscan.get_masklist(f.mask, row=r)
             msg = "mask range: "+str(masklist)
             asaplog.push(msg, False)
@@ -2068 +2203 @@
                     self.blpars.append(None)
                     self.masklists.append(None)
+                    self.actualmask.append(None)
                     continue
 
@@ -2073 +2209 @@
             self.blpars.append(f.get_parameters())
             self.masklists.append(masklist)
+            self.actualmask.append(f.mask)
         if plot:
             f._p.unmap()

trunk/src/STMath.cpp

@@ -1796 +1796 @@
     out.push_back(outstat);
   }
+  return out;
+}
+
+std::vector< float > STMath::statisticRow( const CountedPtr< Scantable > & in,
+                                        const std::vector< bool > & mask,
+                                        const std::string& which,
+                                        int row )
+{
+
+  Vector<Bool> m(mask);
+  const Table& tab = in->table();
+  ROArrayColumn<Float> specCol(tab, "SPECTRA");
+  ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
+  std::vector<float> out;
+
+  Vector<Float> spec; specCol.get(row, spec);
+  Vector<uChar> flag; flagCol.get(row, flag);
+  MaskedArray<Float> ma  = maskedArray(spec, flag);
+  float outstat = 0.0;
+  if ( spec.nelements() == m.nelements() ) {
+    outstat = mathutil::statistics(which, ma(m));
+  } else {
+    outstat = mathutil::statistics(which, ma);
+  }
+  out.push_back(outstat);
+
   return out;
 }

trunk/src/STMath.h

@@ -277 +277 @@
                                const std::vector<bool>& mask,
                                const std::string& which);
+
+  std::vector<float> statisticRow(const casa::CountedPtr<Scantable>& in,
+                               const std::vector<bool>& mask,
+                               const std::string& which,
+                               int row);
 
   std::vector< int > minMaxChan(const casa::CountedPtr<Scantable>& in,

trunk/src/STMathWrapper.h

@@ -132 +132 @@
   { return STMath::statistic(in.getCP(), mask, which); }
 
+  std::vector<float> statisticRow(const ScantableWrapper& in,
+                               const std::vector<bool>& mask,
+                               const std::string& which,
+                               int row)
+  { return STMath::statisticRow(in.getCP(), mask, which, row); }
+
   std::vector<int> minMaxChan(const ScantableWrapper& in,
                                const std::vector<bool>& mask,

trunk/src/Scantable.cpp

@@ -1712 +1712 @@
 }
 
+bool Scantable::getFlagtraFast(int whichrow)
+{
+  uChar flag;
+  Vector<uChar> flags;
+  flagsCol_.get(uInt(whichrow), flags);
+  for (int i = 0; i < flags.size(); i++) {
+    if (i==0) {
+      flag = flags[i];
+    }
+    else {
+      flag &= flags[i];
+    }
+    return ((flag >> 7) == 1);
+   }
+}
+
+void Scantable::doPolyBaseline(const std::vector<bool>& mask, int order, int rowno, Fitter& fitter)
+{
+  fitter.setExpression("poly", order);
+
+  std::vector<double> abcsd = getAbcissa(rowno);
+  std::vector<float> abcs;
+  for (int i = 0; i < abcsd.size(); i++) {
+    abcs.push_back((float)abcsd[i]);
+  }
+  std::vector<float> spec = getSpectrum(rowno);
+  std::vector<bool> fmask = getMask(rowno);
+  for (int i = 0; i < fmask.size(); i++) {
+    fmask[i] = fmask[i] && mask[i];
+  }
+  fitter.setData(abcs, spec, fmask);
+
+  fitter.lfit();
+}
+
+void Scantable::polyBaseline(const std::vector<bool>& mask, int order, int rowno)
+{
+  Fitter fitter = Fitter();
+  doPolyBaseline(mask, order, rowno, fitter);
+  setSpectrum(fitter.getResidual(), rowno);
+}
+
+void Scantable::polyBaseline(const std::vector<bool>& mask, int order, int rowno, int pars_ptr, int pars_size, int errs_ptr, int errs_size, int fmask_ptr, int fmask_size)
+{
+  Fitter fitter = Fitter();
+  doPolyBaseline(mask, order, rowno, fitter);
+  setSpectrum(fitter.getResidual(), rowno);
+
+  std::vector<float> pars = fitter.getParameters();
+  if (pars_size != pars.size()) {
+    throw(AipsError("wrong pars size"));
+  }
+  float *ppars = reinterpret_cast<float*>(pars_ptr);
+  for (int i = 0; i < pars_size; i++) {
+    ppars[i] = pars[i];
+  }
+
+  std::vector<float> errs = fitter.getErrors();
+  if (errs_size != errs.size()) {
+    throw(AipsError("wrong errors size"));
+  }
+  float *perrs = reinterpret_cast<float*>(errs_ptr);
+  for (int i = 0; i < errs_size; i++) {
+    perrs[i] = errs[i];
+  }
+
+  if (fmask_size != fmask.size()) {
+    throw(AipsError("wrong fmask size"));
+  }
+  int *pfmask = reinterpret_cast<int*>(fmask_ptr);
+  for (int i = 0; i < fmask_size; i++) {
+    pfmask[i] = (fmask[i] ? 1 : 0);
+  }
+}
+
 }
 //namespace asap

trunk/src/Scantable.h

@@ -47 +47 @@
 #include "STFit.h"
 #include "STFitEntry.h"
+#include "STFitter.h"
 
 namespace asap {
@@ -486 +487 @@
   void regridChannel( int nchan, double dnu, int irow ) ;
 
+  bool getFlagtraFast(int whichrow);
+
+  void polyBaseline(const std::vector<bool>& mask, int order, int rowno);
+  void polyBaseline(const std::vector<bool>& mask, int order, int rowno, int pars_ptr, int pars_size, int errs_ptr, int errs_size, int fmask_ptr, int fmask_size);
+
 
 private:
@@ -506 +512 @@
    */
   std::string formatDirection(const casa::MDirection& md) const;
-
 
   /**
@@ -598 +603 @@
                                                       const casa::String&,
                                                       const casa::Array<T2>&);
+
+  void doPolyBaseline(const std::vector<bool>& mask, int order, int rowno, Fitter& fitter);
 };
 

trunk/src/ScantableWrapper.h

@@ -250 +250 @@
   { table_->reshapeSpectrum( nmin, nmax ); }
 
+  void polyBaseline(const std::vector<bool>& mask, int order, int rowno, int pars_ptr, int pars_size, int errs_ptr, int errs_size, int fmask_ptr, int fmask_size)
+  { table_->polyBaseline(mask, order, rowno, pars_ptr, pars_size, errs_ptr, errs_size, fmask_ptr, fmask_size); }
+
+  void polyBaseline(const std::vector<bool>& mask, int order, int rowno)
+  { table_->polyBaseline(mask, order, rowno); }
+
+  bool getFlagtraFast(int whichrow=0) const
+    { return table_->getFlagtraFast(whichrow); }
+
+
 private:
   casa::CountedPtr<Scantable> table_;

trunk/src/python_STMath.cpp

@@ -59 +59 @@
         .def("_dofs", &STMathWrapper::dofs)
         .def("_stats", &STMathWrapper::statistic)
+        .def("_statsrow", &STMathWrapper::statisticRow)
         .def("_minmaxchan", &STMathWrapper::minMaxChan)
         .def("_freqswitch", &STMathWrapper::freqSwitch)

trunk/src/python_Scantable.cpp

@@ -139 +139 @@
          (boost::python::arg("nmin")=-1, 
           boost::python::arg("nmax")=-1) )
+    .def("_poly_baseline", &ScantableWrapper::polyBaseline)
+    .def("_getflagtrafast", &ScantableWrapper::getFlagtraFast,
+         (boost::python::arg("whichrow")=0) )
   ;
 };