Changeset 1907 for trunk/python


Ignore:
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:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/python/scantable.py

    r1904 r1907  
    1616from asap.coordinate import coordinate
    1717from asap.utils import _n_bools, mask_not, mask_and, mask_or, page
     18from asap.asapfitter import fitter
    1819
    1920
     
    486487
    487488    @asaplog_post_dec
    488     def stats(self, stat='stddev', mask=None, form='3.3f'):
     489    def stats(self, stat='stddev', mask=None, form='3.3f', row=None):
    489490        """\
    490491        Determine the specified statistic of the current beam/if/pol
     
    502503            form:    format string to print statistic values
    503504
    504         Example::
    505 
     505            row:     row number of spectrum to process.
     506                     (default is None: for all rows)
     507
     508        Example:
    506509            scan.set_unit('channel')
    507510            msk = scan.create_mask([100, 200], [500, 600])
     
    521524            getchan = True
    522525            statvals = []
    523         if not rtnabc: statvals = self._math._stats(self, mask, stat)
     526        if not rtnabc:
     527            if row == None:
     528                statvals = self._math._stats(self, mask, stat)
     529            else:
     530                statvals = self._math._statsrow(self, mask, stat, int(row))
    524531
    525532        #def cb(i):
     
    533540        #outvec = []
    534541        sep = '-'*50
    535         for i in range(self.nrow()):
     542
     543        if row == None:
     544            rows = xrange(self.nrow())
     545        elif isinstance(row, int):
     546            rows = [ row ]
     547
     548        for i in rows:
    536549            refstr = ''
    537550            statunit= ''
     
    549562            out += 'Scan[%d] (%s) ' % (self.getscan(i), src)
    550563            out += 'Time[%s]:\n' % (tm)
    551             if self.nbeam(-1) > 1:
    552                 out +=  ' Beam[%d] ' % (self.getbeam(i))
    553             if self.nif(-1) > 1: out +=  ' IF[%d] ' % (self.getif(i))
    554             if self.npol(-1) > 1: out +=  ' Pol[%d] ' % (self.getpol(i))
     564            if self.nbeam(-1) > 1: out +=  ' Beam[%d] ' % (self.getbeam(i))
     565            if self.nif(-1) > 1:   out +=  ' IF[%d] ' % (self.getif(i))
     566            if self.npol(-1) > 1:  out +=  ' Pol[%d] ' % (self.getpol(i))
    555567            #outvec.append(callback(i))
    556             #out += ('= %'+form) % (outvec[i]) +'   '+refstr+'\n'
    557             out += ('= %'+form) % (statvals[i]) +'   '+refstr+'\n'
     568            if len(rows) > 1:
     569                # out += ('= %'+form) % (outvec[i]) +'   '+refstr+'\n'
     570                out += ('= %'+form) % (statvals[i]) +'   '+refstr+'\n'
     571            else:
     572                # out += ('= %'+form) % (outvec[0]) +'   '+refstr+'\n'
     573                out += ('= %'+form) % (statvals[0]) +'   '+refstr+'\n'
    558574            out +=  sep+"\n"
    559 
    560575
    561576        import os
     
    18371852
    18381853    @asaplog_post_dec
    1839     def poly_baseline(self, mask=None, order=0, plot=False, uselin=False,
    1840                       insitu=None):
     1854    def old_poly_baseline(self, mask=None, order=0, plot=False, uselin=False, insitu=None, rows=None):
    18411855        """\
    18421856        Return a scan which has been baselined (all rows) by a polynomial.
    1843 
     1857       
    18441858        Parameters:
    18451859
     
    18581872                        The default is taken from .asaprc (False)
    18591873
    1860         Example::
    1861 
     1874            rows:       row numbers of spectra to be processed.
     1875                        (default is None: for all rows)
     1876       
     1877        Example:
    18621878            # return a scan baselined by a third order polynomial,
    18631879            # not using a mask
     
    18721888        varlist = vars()
    18731889        if mask is None:
    1874             mask = [True for i in xrange(self.nchan(-1))]
    1875 
    1876         from asap.asapfitter import fitter
     1890            mask = [True for i in xrange(self.nchan())]
     1891
    18771892        try:
    18781893            f = fitter()
     
    18821897                f.set_function(poly=order)
    18831898
    1884             rows = range(workscan.nrow())
     1899            if rows == None:
     1900                rows = xrange(workscan.nrow())
     1901            elif isinstance(rows, int):
     1902                rows = [ rows ]
     1903           
    18851904            if len(rows) > 0:
    18861905                self.blpars = []
    1887 
     1906                self.masklists = []
     1907                self.actualmask = []
     1908           
    18881909            for r in rows:
    1889                 # take into account flagtra info (CAS-1434)
    1890                 flagtra = workscan._getmask(r)
    1891                 actualmask = mask[:]
    1892                 if len(actualmask) == 0:
    1893                     actualmask = list(flagtra[:])
    1894                 else:
    1895                     if len(actualmask) != len(flagtra):
    1896                         raise RuntimeError, "Mask and flagtra have different length"
    1897                     else:
    1898                         for i in range(0, len(actualmask)):
    1899                             actualmask[i] = actualmask[i] and flagtra[i]
    1900                 f.set_scan(workscan, actualmask)
    19011910                f.x = workscan._getabcissa(r)
    19021911                f.y = workscan._getspectrum(r)
     1912                f.mask = mask_and(mask, workscan._getmask(r))    # (CAS-1434)
    19031913                f.data = None
    19041914                f.fit()
     
    19081918                    if x.upper() == 'N':
    19091919                        self.blpars.append(None)
     1920                        self.masklists.append(None)
     1921                        self.actualmask.append(None)
    19101922                        continue
    19111923                workscan._setspectrum(f.fitter.getresidual(), r)
    19121924                self.blpars.append(f.get_parameters())
     1925                self.masklists.append(workscan.get_masklist(f.mask, row=r))
     1926                self.actualmask.append(f.mask)
    19131927
    19141928            if plot:
     
    19251939
    19261940
    1927     def auto_poly_baseline(self, mask=[], edge=(0, 0), order=0,
     1941    def poly_baseline(self, mask=None, order=0, plot=False, batch=False, insitu=None, rows=None):
     1942        """\
     1943        Return a scan which has been baselined (all rows) by a polynomial.
     1944        Parameters:
     1945            mask:       an optional mask
     1946            order:      the order of the polynomial (default is 0)
     1947            plot:       plot the fit and the residual. In this each
     1948                        indivual fit has to be approved, by typing 'y'
     1949                        or 'n'. Ignored if batch = True.
     1950            batch:      if True a faster algorithm is used and logs
     1951                        including the fit results are not output
     1952                        (default is False)
     1953            insitu:     if False a new scantable is returned.
     1954                        Otherwise, the scaling is done in-situ
     1955                        The default is taken from .asaprc (False)
     1956            rows:       row numbers of spectra to be processed.
     1957                        (default is None: for all rows)
     1958        Example:
     1959            # return a scan baselined by a third order polynomial,
     1960            # not using a mask
     1961            bscan = scan.poly_baseline(order=3)
     1962        """
     1963        if insitu is None: insitu = rcParams["insitu"]
     1964        if insitu:
     1965            workscan = self
     1966        else:
     1967            workscan = self.copy()
     1968
     1969        varlist = vars()
     1970        nchan = workscan.nchan()
     1971       
     1972        if mask is None:
     1973            mask = [True for i in xrange(nchan)]
     1974
     1975        try:
     1976            if rows == None:
     1977                rows = xrange(workscan.nrow())
     1978            elif isinstance(rows, int):
     1979                rows = [ rows ]
     1980           
     1981            if len(rows) > 0:
     1982                self.blpars = []
     1983                self.masklists = []
     1984                self.actualmask = []
     1985
     1986            if batch:
     1987                for r in rows:
     1988                    workscan._poly_baseline(mask, order, r)
     1989            elif plot:
     1990                f = fitter()
     1991                f.set_function(lpoly=order)
     1992                for r in rows:
     1993                    f.x = workscan._getabcissa(r)
     1994                    f.y = workscan._getspectrum(r)
     1995                    f.mask = mask_and(mask, workscan._getmask(r))    # (CAS-1434)
     1996                    f.data = None
     1997                    f.fit()
     1998                   
     1999                    f.plot(residual=True)
     2000                    accept_fit = raw_input("Accept fit ( [y]/n ): ")
     2001                    if accept_fit.upper() == "N":
     2002                        self.blpars.append(None)
     2003                        self.masklists.append(None)
     2004                        self.actualmask.append(None)
     2005                        continue
     2006                    workscan._setspectrum(f.fitter.getresidual(), r)
     2007                    self.blpars.append(f.get_parameters())
     2008                    self.masklists.append(workscan.get_masklist(f.mask, row=r))
     2009                    self.actualmask.append(f.mask)
     2010                   
     2011                f._p.unmap()
     2012                f._p = None
     2013            else:
     2014                import array
     2015                for r in rows:
     2016                    pars = array.array("f", [0.0 for i in range(order+1)])
     2017                    pars_adr = pars.buffer_info()[0]
     2018                    pars_len = pars.buffer_info()[1]
     2019                   
     2020                    errs = array.array("f", [0.0 for i in range(order+1)])
     2021                    errs_adr = errs.buffer_info()[0]
     2022                    errs_len = errs.buffer_info()[1]
     2023                   
     2024                    fmsk = array.array("i", [1 for i in range(nchan)])
     2025                    fmsk_adr = fmsk.buffer_info()[0]
     2026                    fmsk_len = fmsk.buffer_info()[1]
     2027                   
     2028                    workscan._poly_baseline(mask, order, r, pars_adr, pars_len, errs_adr, errs_len, fmsk_adr, fmsk_len)
     2029                   
     2030                    params = pars.tolist()
     2031                    fmtd = ""
     2032                    for i in xrange(len(params)): fmtd += "  p%d= %3.6f," % (i, params[i])
     2033                    fmtd = fmtd[:-1]  # remove trailing ","
     2034                    errors = errs.tolist()
     2035                    fmask = fmsk.tolist()
     2036                    for i in xrange(len(fmask)): fmask[i] = (fmask[i] > 0)    # transform (1/0) -> (True/False)
     2037                   
     2038                    self.blpars.append({"params":params, "fixed":[], "formatted":fmtd, "errors":errors})
     2039                    self.masklists.append(workscan.get_masklist(fmask, r))
     2040                    self.actualmask.append(fmask)
     2041                   
     2042                    asaplog.push(str(fmtd))
     2043           
     2044            workscan._add_history("poly_baseline", varlist)
     2045            print_log()
     2046           
     2047            if insitu:
     2048                self._assign(workscan)
     2049            else:
     2050                return workscan
     2051           
     2052        except RuntimeError:
     2053            msg = "The fit failed, possibly because it didn't converge."
     2054            if rcParams["verbose"]:
     2055                print_log()
     2056                asaplog.push(str(msg))
     2057                print_log("ERROR")
     2058                return
     2059            else:
     2060                raise RuntimeError(msg)
     2061
     2062
     2063    def auto_poly_baseline(self, mask=None, edge=(0, 0), order=0,
    19282064                           threshold=3, chan_avg_limit=1, plot=False,
    1929                            insitu=None):
     2065                           insitu=None, rows=None):
    19302066        """\
    19312067        Return a scan which has been baselined (all rows) by a polynomial.
     
    19682104                        Otherwise, the scaling is done in-situ
    19692105                        The default is taken from .asaprc (False)
     2106            rows:       row numbers of spectra to be processed.
     2107                        (default is None: for all rows)
    19702108
    19712109
     
    19772115        if insitu is None: insitu = rcParams['insitu']
    19782116        varlist = vars()
    1979         from asap.asapfitter import fitter
    19802117        from asap.asaplinefind import linefinder
    19812118        from asap import _is_sequence_or_number as _is_valid
     
    20012138            curedge = edge;
    20022139
     2140        if not insitu:
     2141            workscan = self.copy()
     2142        else:
     2143            workscan = self
     2144
    20032145        # setup fitter
    20042146        f = fitter()
    2005         f.set_function(poly=order)
     2147        f.set_function(lpoly=order)
    20062148
    20072149        # setup line finder
     
    20092151        fl.set_options(threshold=threshold,avg_limit=chan_avg_limit)
    20102152
    2011         if not insitu:
    2012             workscan = self.copy()
    2013         else:
    2014             workscan = self
    2015 
    20162153        fl.set_scan(workscan)
    20172154
    2018         rows = range(workscan.nrow())
     2155        if mask is None:
     2156            mask = _n_bools(workscan.nchan(), True)
     2157       
     2158        if rows is None:
     2159            rows = xrange(workscan.nrow())
     2160        elif isinstance(rows, int):
     2161            rows = [ rows ]
     2162       
    20192163        # Save parameters of baseline fits & masklists as a class attribute.
    20202164        # NOTICE: It does not reflect changes in scantable!
     
    20222166            self.blpars=[]
    20232167            self.masklists=[]
     2168            self.actualmask=[]
    20242169        asaplog.push("Processing:")
    20252170        for r in rows:
     
    20362181                    curedge = edge[workscan.getif(r)]
    20372182
    2038             # take into account flagtra info (CAS-1434)
    2039             flagtra = workscan._getmask(r)
    2040             actualmask = mask[:]
    2041             if len(actualmask) == 0:
    2042                 actualmask = list(flagtra[:])
    2043             else:
    2044                 if len(actualmask) != len(flagtra):
    2045                     raise RuntimeError, "Mask and flagtra have different length"
    2046                 else:
    2047                     for i in range(0, len(actualmask)):
    2048                         actualmask[i] = actualmask[i] and flagtra[i]
     2183            actualmask = mask_and(mask, workscan._getmask(r))    # (CAS-1434)
    20492184
    20502185            # setup line finder
    20512186            fl.find_lines(r, actualmask, curedge)
    2052             outmask=fl.get_mask()
    2053             f.set_scan(workscan, fl.get_mask())
     2187           
    20542188            f.x = workscan._getabcissa(r)
    20552189            f.y = workscan._getspectrum(r)
     2190            f.mask = fl.get_mask()
    20562191            f.data = None
    20572192            f.fit()
    20582193
    20592194            # Show mask list
    2060             masklist=workscan.get_masklist(fl.get_mask(),row=r)
     2195            masklist=workscan.get_masklist(f.mask, row=r)
    20612196            msg = "mask range: "+str(masklist)
    20622197            asaplog.push(msg, False)
     
    20682203                    self.blpars.append(None)
    20692204                    self.masklists.append(None)
     2205                    self.actualmask.append(None)
    20702206                    continue
    20712207
     
    20732209            self.blpars.append(f.get_parameters())
    20742210            self.masklists.append(masklist)
     2211            self.actualmask.append(f.mask)
    20752212        if plot:
    20762213            f._p.unmap()
Note: See TracChangeset for help on using the changeset viewer.