Changeset 1907 for trunk/python

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().

File:

• trunk/python/scantable.py (modified) (20 diffs)

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()