Ignore:
Timestamp:
02/25/11 16:51:50 (13 years ago)
Author:
WataruKawasaki
Message:

New Development: Yes

JIRA Issue: Yes (CAS-2373, CAS-2620)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: For Scantable::polyBaseline(), parameters and return value have been changed.

Test Programs:

Put in Release Notes: Yes

Module(s): sdbaseline, sd.linefinder

Description: (1) CAS-2373-related:

(1.1) Modified Scantable::polyBaseline() to have the row-based loop inside.

Now it fits and subtracts baseline for all rows and also output info
about the fitting result to logger and text file, while in the
previous version this method just did baseline fitting/subtraction
for one row only and had to be put inside a row-based loop at the
python side ("poly_baseline()" in scantable.py) and result output had
also to be controlled at the python side. Using a test data from NRO
45m telescope (348,000 rows, 512 channels), the processing time of
scantable.poly_baseline() has reduced from 130 minutes to 5-6 minutes.

(1.2) For accelerating the another polynomial fitting function, namely

scantable.auto_poly_baseline(), added a method
Scantable::autoPolyBaseline(). This basically does the same thing
with Scantable::polyBaseline(), but uses linefinfer also to
automatically flag the line regions.

(1.3) To make linefinder usable in Scantable class, added a method

linefinder.setdata(). This makes it possible to apply linefinder
for a float-list data given as a parameter, without setting scantable,
while it was indispensable to set scantable to use linefinger previously.

(2) CAS-2620-related:

Added Scantable::cubicSplineBaseline() and autoCubicSplineBaseline() for
fit baseline using the cubic spline function. Parameters include npiece
(number of polynomial pieces), clipthresh (clipping threshold), and
clipniter (maximum iteration number).



File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/python/scantable.py

    r2008 r2012  
    18981898        else: return s
    18991899
     1900
     1901    @asaplog_post_dec
     1902    def cspline_baseline(self, insitu=None, mask=None, npiece=None, clipthresh=None, clipniter=None, plot=None, outlog=None, blfile=None):
     1903        """\
     1904        Return a scan which has been baselined (all rows) by cubic spline function (piecewise cubic polynomial).
     1905        Parameters:
     1906            insitu:     If False a new scantable is returned.
     1907                        Otherwise, the scaling is done in-situ
     1908                        The default is taken from .asaprc (False)
     1909            mask:       An optional mask
     1910            npiece:     Number of pieces. (default is 2)
     1911            clipthresh: Clipping threshold. (default is 3.0, unit: sigma)
     1912            clipniter:  maximum number of iteration of 'clipthresh'-sigma clipping (default is 1)
     1913            plot:   *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
     1914                        plot the fit and the residual. In this each
     1915                        indivual fit has to be approved, by typing 'y'
     1916                        or 'n'
     1917            outlog:     Output the coefficients of the best-fit
     1918                        function to logger (default is False)
     1919            blfile:     Name of a text file in which the best-fit
     1920                        parameter values to be written
     1921                        (default is "": no file/logger output)
     1922
     1923        Example:
     1924            # return a scan baselined by a cubic spline consisting of 2 pieces (i.e., 1 internal knot),
     1925            # also with 3-sigma clipping, iteration up to 4 times
     1926            bscan = scan.cspline_baseline(npiece=2,clipthresh=3.0,clipniter=4)
     1927        """
     1928       
     1929        varlist = vars()
     1930       
     1931        if insitu is None: insitu = rcParams["insitu"]
     1932        if insitu:
     1933            workscan = self
     1934        else:
     1935            workscan = self.copy()
     1936
     1937        nchan = workscan.nchan()
     1938       
     1939        if mask is None: mask = [True for i in xrange(nchan)]
     1940        if npiece is None: npiece = 2
     1941        if clipthresh is None: clipthresh = 3.0
     1942        if clipniter is None: clipniter = 1
     1943        if plot is None: plot = False
     1944        if outlog is None: outlog = False
     1945        if blfile is None: blfile = ""
     1946
     1947        outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
     1948       
     1949        try:
     1950            #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method.
     1951            workscan._cspline_baseline(mask, npiece, clipthresh, clipniter, outlog, blfile)
     1952           
     1953            workscan._add_history("cspline_baseline", varlist)
     1954           
     1955            if insitu:
     1956                self._assign(workscan)
     1957            else:
     1958                return workscan
     1959           
     1960        except RuntimeError, e:
     1961            msg = "The fit failed, possibly because it didn't converge."
     1962            if rcParams["verbose"]:
     1963                asaplog.push(str(e))
     1964                asaplog.push(str(msg))
     1965                return
     1966            else:
     1967                raise RuntimeError(str(e)+'\n'+msg)
     1968
     1969
     1970    def auto_cspline_baseline(self, insitu=None, mask=None, npiece=None, clipthresh=None,
     1971                              clipniter=None, edge=None, threshold=None,
     1972                              chan_avg_limit=None, plot=None, outlog=None, blfile=None):
     1973        """\
     1974        Return a scan which has been baselined (all rows) by cubic spline
     1975        function (piecewise cubic polynomial).
     1976        Spectral lines are detected first using linefinder and masked out
     1977        to avoid them affecting the baseline solution.
     1978
     1979        Parameters:
     1980            insitu:     if False a new scantable is returned.
     1981                        Otherwise, the scaling is done in-situ
     1982                        The default is taken from .asaprc (False)
     1983            mask:       an optional mask retreived from scantable
     1984            npiece:     Number of pieces. (default is 2)
     1985            clipthresh: Clipping threshold. (default is 3.0, unit: sigma)
     1986            clipniter:  maximum number of iteration of 'clipthresh'-sigma clipping (default is 1)
     1987            edge:       an optional number of channel to drop at
     1988                        the edge of spectrum. If only one value is
     1989                        specified, the same number will be dropped
     1990                        from both sides of the spectrum. Default
     1991                        is to keep all channels. Nested tuples
     1992                        represent individual edge selection for
     1993                        different IFs (a number of spectral channels
     1994                        can be different)
     1995            threshold:  the threshold used by line finder. It is
     1996                        better to keep it large as only strong lines
     1997                        affect the baseline solution.
     1998            chan_avg_limit:
     1999                        a maximum number of consequtive spectral
     2000                        channels to average during the search of
     2001                        weak and broad lines. The default is no
     2002                        averaging (and no search for weak lines).
     2003                        If such lines can affect the fitted baseline
     2004                        (e.g. a high order polynomial is fitted),
     2005                        increase this parameter (usually values up
     2006                        to 8 are reasonable). Most users of this
     2007                        method should find the default value sufficient.
     2008            plot:   *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
     2009                        plot the fit and the residual. In this each
     2010                        indivual fit has to be approved, by typing 'y'
     2011                        or 'n'
     2012            outlog:     Output the coefficients of the best-fit
     2013                        function to logger (default is False)
     2014            blfile:     Name of a text file in which the best-fit
     2015                        parameter values to be written
     2016                        (default is "": no file/logger output)
     2017
     2018        Example:
     2019            bscan = scan.auto_cspline_baseline(npiece=3, insitu=False)
     2020        """
     2021
     2022        varlist = vars()
     2023
     2024        if insitu is None: insitu = rcParams['insitu']
     2025        if insitu:
     2026            workscan = self
     2027        else:
     2028            workscan = self.copy()
     2029
     2030        nchan = workscan.nchan()
     2031       
     2032        if mask is None: mask = [True for i in xrange(nchan)]
     2033        if npiece is None: npiece = 2
     2034        if clipthresh is None: clipthresh = 3.0
     2035        if clipniter is None: clipniter = 1
     2036        if edge is None: edge = (0, 0)
     2037        if threshold is None: threshold = 3
     2038        if chan_avg_limit is None: chan_avg_limit = 1
     2039        if plot is None: plot = False
     2040        if outlog is None: outlog = False
     2041        if blfile is None: blfile = ""
     2042
     2043        outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
     2044       
     2045        from asap.asaplinefind import linefinder
     2046        from asap import _is_sequence_or_number as _is_valid
     2047
     2048        if not (isinstance(edge, list) or isinstance(edge, tuple)): edge = [ edge ]
     2049        individualedge = False;
     2050        if len(edge) > 1: individualedge = isinstance(edge[0], list) or isinstance(edge[0], tuple)
     2051
     2052        if individualedge:
     2053            for edgepar in edge:
     2054                if not _is_valid(edgepar, int):
     2055                    raise ValueError, "Each element of the 'edge' tuple has \
     2056                                       to be a pair of integers or an integer."
     2057        else:
     2058            if not _is_valid(edge, int):
     2059                raise ValueError, "Parameter 'edge' has to be an integer or a \
     2060                                   pair of integers specified as a tuple. \
     2061                                   Nested tuples are allowed \
     2062                                   to make individual selection for different IFs."
     2063
     2064            if len(edge) > 1:
     2065                curedge = edge
     2066            else:
     2067                curedge = edge + edge
     2068
     2069        try:
     2070            #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method.
     2071            if individualedge:
     2072                curedge = []
     2073                for i in xrange(len(edge)):
     2074                    curedge += edge[i]
     2075               
     2076            workscan._auto_cspline_baseline(mask, npiece, clipthresh, clipniter, curedge, threshold, chan_avg_limit, outlog, blfile)
     2077
     2078            workscan._add_history("auto_cspline_baseline", varlist)
     2079           
     2080            if insitu:
     2081                self._assign(workscan)
     2082            else:
     2083                return workscan
     2084           
     2085        except RuntimeError, e:
     2086            msg = "The fit failed, possibly because it didn't converge."
     2087            if rcParams["verbose"]:
     2088                asaplog.push(str(e))
     2089                asaplog.push(str(msg))
     2090                return
     2091            else:
     2092                raise RuntimeError(str(e)+'\n'+msg)
     2093
     2094
     2095    @asaplog_post_dec
     2096    def poly_baseline(self, insitu=None, mask=None, order=None, plot=None, outlog=None, blfile=None):
     2097        """\
     2098        Return a scan which has been baselined (all rows) by a polynomial.
     2099        Parameters:
     2100            insitu:     if False a new scantable is returned.
     2101                        Otherwise, the scaling is done in-situ
     2102                        The default is taken from .asaprc (False)
     2103            mask:       an optional mask
     2104            order:      the order of the polynomial (default is 0)
     2105            plot:       plot the fit and the residual. In this each
     2106                        indivual fit has to be approved, by typing 'y'
     2107                        or 'n'
     2108            outlog:     Output the coefficients of the best-fit
     2109                        function to logger (default is False)
     2110            blfile:     Name of a text file in which the best-fit
     2111                        parameter values to be written
     2112                        (default is "": no file/logger output)
     2113
     2114        Example:
     2115            # return a scan baselined by a third order polynomial,
     2116            # not using a mask
     2117            bscan = scan.poly_baseline(order=3)
     2118        """
     2119       
     2120        varlist = vars()
     2121       
     2122        if insitu is None: insitu = rcParams["insitu"]
     2123        if insitu:
     2124            workscan = self
     2125        else:
     2126            workscan = self.copy()
     2127
     2128        nchan = workscan.nchan()
     2129       
     2130        if mask is None: mask = [True for i in xrange(nchan)]
     2131        if order is None: order = 0
     2132        if plot is None: plot = False
     2133        if outlog is None: outlog = False
     2134        if blfile is None: blfile = ""
     2135
     2136        outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
     2137       
     2138        try:
     2139            rows = xrange(workscan.nrow())
     2140           
     2141            #if len(rows) > 0: workscan._init_blinfo()
     2142
     2143            if plot:
     2144                if outblfile: blf = open(blfile, "a")
     2145               
     2146                f = fitter()
     2147                f.set_function(lpoly=order)
     2148                for r in rows:
     2149                    f.x = workscan._getabcissa(r)
     2150                    f.y = workscan._getspectrum(r)
     2151                    f.mask = mask_and(mask, workscan._getmask(r))    # (CAS-1434)
     2152                    f.data = None
     2153                    f.fit()
     2154                   
     2155                    f.plot(residual=True)
     2156                    accept_fit = raw_input("Accept fit ( [y]/n ): ")
     2157                    if accept_fit.upper() == "N":
     2158                        #workscan._append_blinfo(None, None, None)
     2159                        continue
     2160                   
     2161                    blpars = f.get_parameters()
     2162                    masklist = workscan.get_masklist(f.mask, row=r, silent=True)
     2163                    #workscan._append_blinfo(blpars, masklist, f.mask)
     2164                    workscan._setspectrum(f.fitter.getresidual(), r)
     2165                   
     2166                    if outblfile:
     2167                        rms = workscan.get_rms(f.mask, r)
     2168                        dataout = workscan.format_blparams_row(blpars["params"], blpars["fixed"], rms, str(masklist), r, True)
     2169                        blf.write(dataout)
     2170
     2171                f._p.unmap()
     2172                f._p = None
     2173
     2174                if outblfile: blf.close()
     2175            else:
     2176                workscan._poly_baseline(mask, order, outlog, blfile)
     2177           
     2178            workscan._add_history("poly_baseline", varlist)
     2179           
     2180            if insitu:
     2181                self._assign(workscan)
     2182            else:
     2183                return workscan
     2184           
     2185        except RuntimeError, e:
     2186            msg = "The fit failed, possibly because it didn't converge."
     2187            if rcParams["verbose"]:
     2188                asaplog.push(str(e))
     2189                asaplog.push(str(msg))
     2190                return
     2191            else:
     2192                raise RuntimeError(str(e)+'\n'+msg)
     2193
     2194
     2195    def auto_poly_baseline(self, insitu=None, mask=None, order=None, edge=None, threshold=None,
     2196                           chan_avg_limit=None, plot=None, outlog=None, blfile=None):
     2197        """\
     2198        Return a scan which has been baselined (all rows) by a polynomial.
     2199        Spectral lines are detected first using linefinder and masked out
     2200        to avoid them affecting the baseline solution.
     2201
     2202        Parameters:
     2203            insitu:     if False a new scantable is returned.
     2204                        Otherwise, the scaling is done in-situ
     2205                        The default is taken from .asaprc (False)
     2206            mask:       an optional mask retreived from scantable
     2207            order:      the order of the polynomial (default is 0)
     2208            edge:       an optional number of channel to drop at
     2209                        the edge of spectrum. If only one value is
     2210                        specified, the same number will be dropped
     2211                        from both sides of the spectrum. Default
     2212                        is to keep all channels. Nested tuples
     2213                        represent individual edge selection for
     2214                        different IFs (a number of spectral channels
     2215                        can be different)
     2216            threshold:  the threshold used by line finder. It is
     2217                        better to keep it large as only strong lines
     2218                        affect the baseline solution.
     2219            chan_avg_limit:
     2220                        a maximum number of consequtive spectral
     2221                        channels to average during the search of
     2222                        weak and broad lines. The default is no
     2223                        averaging (and no search for weak lines).
     2224                        If such lines can affect the fitted baseline
     2225                        (e.g. a high order polynomial is fitted),
     2226                        increase this parameter (usually values up
     2227                        to 8 are reasonable). Most users of this
     2228                        method should find the default value sufficient.
     2229            plot:       plot the fit and the residual. In this each
     2230                        indivual fit has to be approved, by typing 'y'
     2231                        or 'n'
     2232            outlog:     Output the coefficients of the best-fit
     2233                        function to logger (default is False)
     2234            blfile:     Name of a text file in which the best-fit
     2235                        parameter values to be written
     2236                        (default is "": no file/logger output)
     2237
     2238        Example:
     2239            bscan = scan.auto_poly_baseline(order=7, insitu=False)
     2240        """
     2241
     2242        varlist = vars()
     2243
     2244        if insitu is None: insitu = rcParams['insitu']
     2245        if insitu:
     2246            workscan = self
     2247        else:
     2248            workscan = self.copy()
     2249
     2250        nchan = workscan.nchan()
     2251       
     2252        if mask is None: mask = [True for i in xrange(nchan)]
     2253        if order is None: order = 0
     2254        if edge is None: edge = (0, 0)
     2255        if threshold is None: threshold = 3
     2256        if chan_avg_limit is None: chan_avg_limit = 1
     2257        if plot is None: plot = False
     2258        if outlog is None: outlog = False
     2259        if blfile is None: blfile = ""
     2260
     2261        outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
     2262       
     2263        from asap.asaplinefind import linefinder
     2264        from asap import _is_sequence_or_number as _is_valid
     2265
     2266        if not (isinstance(edge, list) or isinstance(edge, tuple)): edge = [ edge ]
     2267        individualedge = False;
     2268        if len(edge) > 1: individualedge = isinstance(edge[0], list) or isinstance(edge[0], tuple)
     2269
     2270        if individualedge:
     2271            for edgepar in edge:
     2272                if not _is_valid(edgepar, int):
     2273                    raise ValueError, "Each element of the 'edge' tuple has \
     2274                                       to be a pair of integers or an integer."
     2275        else:
     2276            if not _is_valid(edge, int):
     2277                raise ValueError, "Parameter 'edge' has to be an integer or a \
     2278                                   pair of integers specified as a tuple. \
     2279                                   Nested tuples are allowed \
     2280                                   to make individual selection for different IFs."
     2281
     2282            if len(edge) > 1:
     2283                curedge = edge
     2284            else:
     2285                curedge = edge + edge
     2286
     2287        try:
     2288            rows = xrange(workscan.nrow())
     2289           
     2290            #if len(rows) > 0: workscan._init_blinfo()
     2291
     2292            if plot:
     2293                if outblfile: blf = open(blfile, "a")
     2294               
     2295                fl = linefinder()
     2296                fl.set_options(threshold=threshold,avg_limit=chan_avg_limit)
     2297                fl.set_scan(workscan)
     2298                f = fitter()
     2299                f.set_function(lpoly=order)
     2300
     2301                for r in rows:
     2302                    if individualedge:
     2303                        if len(edge) <= workscan.getif(r):
     2304                            raise RuntimeError, "Number of edge elements appear to " \
     2305                                  "be less than the number of IFs"
     2306                        else:
     2307                            curedge = edge[workscan.getif(r)]
     2308
     2309                    fl.find_lines(r, mask_and(mask, workscan._getmask(r)), curedge)  # (CAS-1434)
     2310
     2311                    f.x = workscan._getabcissa(r)
     2312                    f.y = workscan._getspectrum(r)
     2313                    f.mask = fl.get_mask()
     2314                    f.data = None
     2315                    f.fit()
     2316
     2317                    f.plot(residual=True)
     2318                    accept_fit = raw_input("Accept fit ( [y]/n ): ")
     2319                    if accept_fit.upper() == "N":
     2320                        #workscan._append_blinfo(None, None, None)
     2321                        continue
     2322
     2323                    blpars = f.get_parameters()
     2324                    masklist = workscan.get_masklist(f.mask, row=r, silent=True)
     2325                    #workscan._append_blinfo(blpars, masklist, f.mask)
     2326                    workscan._setspectrum(f.fitter.getresidual(), r)
     2327
     2328                    if outblfile:
     2329                        rms = workscan.get_rms(f.mask, r)
     2330                        dataout = workscan.format_blparams_row(blpars["params"], blpars["fixed"], rms, str(masklist), r, True)
     2331                        blf.write(dataout)
     2332                   
     2333                f._p.unmap()
     2334                f._p = None
     2335
     2336                if outblfile: blf.close()
     2337               
     2338            else:
     2339                if individualedge:
     2340                    curedge = []
     2341                    for i in xrange(len(edge)):
     2342                        curedge += edge[i]
     2343               
     2344                workscan._auto_poly_baseline(mask, order, curedge, threshold, chan_avg_limit, outlog, blfile)
     2345
     2346            workscan._add_history("auto_poly_baseline", varlist)
     2347           
     2348            if insitu:
     2349                self._assign(workscan)
     2350            else:
     2351                return workscan
     2352           
     2353        except RuntimeError, e:
     2354            msg = "The fit failed, possibly because it didn't converge."
     2355            if rcParams["verbose"]:
     2356                asaplog.push(str(e))
     2357                asaplog.push(str(msg))
     2358                return
     2359            else:
     2360                raise RuntimeError(str(e)+'\n'+msg)
     2361
     2362
     2363    ### OBSOLETE ##################################################################
    19002364    @asaplog_post_dec
    19012365    def old_poly_baseline(self, mask=None, order=0, plot=False, uselin=False, insitu=None, rows=None):
    1902         """\
     2366        """
    19032367        Return a scan which has been baselined (all rows) by a polynomial.
    19042368       
     
    19852449            raise RuntimeError(msg)
    19862450
    1987     @asaplog_post_dec
    1988     def poly_baseline(self, mask=None, order=0, plot=False, batch=False, insitu=None, rows=None):
    1989         """\
    1990         Return a scan which has been baselined (all rows) by a polynomial.
    1991         Parameters:
    1992             mask:       an optional mask
    1993             order:      the order of the polynomial (default is 0)
    1994             plot:       plot the fit and the residual. In this each
    1995                         indivual fit has to be approved, by typing 'y'
    1996                         or 'n'. Ignored if batch = True.
    1997             batch:      if True a faster algorithm is used and logs
    1998                         including the fit results are not output
    1999                         (default is False)
    2000             insitu:     if False a new scantable is returned.
    2001                         Otherwise, the scaling is done in-situ
    2002                         The default is taken from .asaprc (False)
    2003             rows:       row numbers of spectra to be baselined.
    2004                         (default is None: for all rows)
    2005         Example:
    2006             # return a scan baselined by a third order polynomial,
    2007             # not using a mask
    2008             bscan = scan.poly_baseline(order=3)
    2009         """
     2451    def _init_blinfo(self):
     2452        """\
     2453        Initialise the following three auxiliary members:
     2454           blpars : parameters of the best-fit baseline,
     2455           masklists : mask data (edge positions of masked channels) and
     2456           actualmask : mask data (in boolean list),
     2457        to keep for use later (including output to logger/text files).
     2458        Used by poly_baseline() and auto_poly_baseline() in case of
     2459        'plot=True'.
     2460        """
     2461        self.blpars = []
     2462        self.masklists = []
     2463        self.actualmask = []
     2464        return
     2465
     2466    def _append_blinfo(self, data_blpars, data_masklists, data_actualmask):
     2467        """\
     2468        Append baseline-fitting related info to blpars, masklist and
     2469        actualmask.
     2470        """
     2471        self.blpars.append(data_blpars)
     2472        self.masklists.append(data_masklists)
     2473        self.actualmask.append(data_actualmask)
     2474        return
    20102475       
    2011         varlist = vars()
    2012        
    2013         if insitu is None: insitu = rcParams["insitu"]
    2014         if insitu:
    2015             workscan = self
    2016         else:
    2017             workscan = self.copy()
    2018 
    2019         nchan = workscan.nchan()
    2020        
    2021         if mask is None:
    2022             mask = [True for i in xrange(nchan)]
    2023 
    2024         try:
    2025             if rows == None:
    2026                 rows = xrange(workscan.nrow())
    2027             elif isinstance(rows, int):
    2028                 rows = [ rows ]
    2029            
    2030             if len(rows) > 0:
    2031                 workscan.blpars = []
    2032                 workscan.masklists = []
    2033                 workscan.actualmask = []
    2034 
    2035             if batch:
    2036                 workscan._poly_baseline_batch(mask, order)
    2037             elif plot:
    2038                 f = fitter()
    2039                 f.set_function(lpoly=order)
    2040                 for r in rows:
    2041                     f.x = workscan._getabcissa(r)
    2042                     f.y = workscan._getspectrum(r)
    2043                     f.mask = mask_and(mask, workscan._getmask(r))    # (CAS-1434)
    2044                     f.data = None
    2045                     f.fit()
    2046                    
    2047                     f.plot(residual=True)
    2048                     accept_fit = raw_input("Accept fit ( [y]/n ): ")
    2049                     if accept_fit.upper() == "N":
    2050                         self.blpars.append(None)
    2051                         self.masklists.append(None)
    2052                         self.actualmask.append(None)
    2053                         continue
    2054                     workscan._setspectrum(f.fitter.getresidual(), r)
    2055                     workscan.blpars.append(f.get_parameters())
    2056                     workscan.masklists.append(workscan.get_masklist(f.mask, row=r))
    2057                     workscan.actualmask.append(f.mask)
    2058                    
    2059                 f._p.unmap()
    2060                 f._p = None
    2061             else:
    2062                 for r in rows:
    2063                     fitparams = workscan._poly_baseline(mask, order, r)
    2064                     params = fitparams.getparameters()
    2065                     fmtd = ", ".join(["p%d = %3.6f" % (i, v) for i, v in enumerate(params)])
    2066                     errors = fitparams.geterrors()
    2067                     fmask = mask_and(mask, workscan._getmask(r))
    2068 
    2069                     workscan.blpars.append({"params":params,
    2070                                             "fixed": fitparams.getfixedparameters(),
    2071                                             "formatted":fmtd, "errors":errors})
    2072                     workscan.masklists.append(workscan.get_masklist(fmask, r, silent=True))
    2073                     workscan.actualmask.append(fmask)
    2074                    
    2075                     asaplog.push(fmtd)
    2076            
    2077             workscan._add_history("poly_baseline", varlist)
    2078            
    2079             if insitu:
    2080                 self._assign(workscan)
    2081             else:
    2082                 return workscan
    2083            
    2084         except RuntimeError, e:
    2085             msg = "The fit failed, possibly because it didn't converge."
    2086             if rcParams["verbose"]:
    2087                 asaplog.push(str(e))
    2088                 asaplog.push(str(msg))
    2089                 return
    2090             else:
    2091                 raise RuntimeError(str(e)+'\n'+msg)
    2092 
    2093 
    2094     def auto_poly_baseline(self, mask=None, edge=(0, 0), order=0,
    2095                            threshold=3, chan_avg_limit=1, plot=False,
    2096                            insitu=None, rows=None):
    2097         """\
    2098         Return a scan which has been baselined (all rows) by a polynomial.
    2099         Spectral lines are detected first using linefinder and masked out
    2100         to avoid them affecting the baseline solution.
    2101 
    2102         Parameters:
    2103 
    2104             mask:       an optional mask retreived from scantable
    2105 
    2106             edge:       an optional number of channel to drop at the edge of
    2107                         spectrum. If only one value is
    2108                         specified, the same number will be dropped from
    2109                         both sides of the spectrum. Default is to keep
    2110                         all channels. Nested tuples represent individual
    2111                         edge selection for different IFs (a number of spectral
    2112                         channels can be different)
    2113 
    2114             order:      the order of the polynomial (default is 0)
    2115 
    2116             threshold:  the threshold used by line finder. It is better to
    2117                         keep it large as only strong lines affect the
    2118                         baseline solution.
    2119 
    2120             chan_avg_limit:
    2121                         a maximum number of consequtive spectral channels to
    2122                         average during the search of weak and broad lines.
    2123                         The default is no averaging (and no search for weak
    2124                         lines). If such lines can affect the fitted baseline
    2125                         (e.g. a high order polynomial is fitted), increase this
    2126                         parameter (usually values up to 8 are reasonable). Most
    2127                         users of this method should find the default value
    2128                         sufficient.
    2129 
    2130             plot:       plot the fit and the residual. In this each
    2131                         indivual fit has to be approved, by typing 'y'
    2132                         or 'n'
    2133 
    2134             insitu:     if False a new scantable is returned.
    2135                         Otherwise, the scaling is done in-situ
    2136                         The default is taken from .asaprc (False)
    2137             rows:       row numbers of spectra to be processed.
    2138                         (default is None: for all rows)
    2139 
    2140 
    2141         Example::
    2142 
    2143             scan2 = scan.auto_poly_baseline(order=7, insitu=False)
    2144 
    2145         """
    2146         if insitu is None: insitu = rcParams['insitu']
    2147         varlist = vars()
    2148         from asap.asaplinefind import linefinder
    2149         from asap import _is_sequence_or_number as _is_valid
    2150 
    2151         # check whether edge is set up for each IF individually
    2152         individualedge = False;
    2153         if len(edge) > 1:
    2154             if isinstance(edge[0], list) or isinstance(edge[0], tuple):
    2155                 individualedge = True;
    2156 
    2157         if not _is_valid(edge, int) and not individualedge:
    2158             raise ValueError, "Parameter 'edge' has to be an integer or a \
    2159             pair of integers specified as a tuple. Nested tuples are allowed \
    2160             to make individual selection for different IFs."
    2161 
    2162         curedge = (0, 0)
    2163         if individualedge:
    2164             for edgepar in edge:
    2165                 if not _is_valid(edgepar, int):
    2166                     raise ValueError, "Each element of the 'edge' tuple has \
    2167                                        to be a pair of integers or an integer."
    2168         else:
    2169             curedge = edge;
    2170 
    2171         if not insitu:
    2172             workscan = self.copy()
    2173         else:
    2174             workscan = self
    2175 
    2176         # setup fitter
    2177         f = fitter()
    2178         f.set_function(lpoly=order)
    2179 
    2180         # setup line finder
    2181         fl = linefinder()
    2182         fl.set_options(threshold=threshold,avg_limit=chan_avg_limit)
    2183 
    2184         fl.set_scan(workscan)
    2185 
    2186         if mask is None:
    2187             mask = _n_bools(workscan.nchan(), True)
    2188        
    2189         if rows is None:
    2190             rows = xrange(workscan.nrow())
    2191         elif isinstance(rows, int):
    2192             rows = [ rows ]
    2193        
    2194         # Save parameters of baseline fits & masklists as a class attribute.
    2195         # NOTICE: It does not reflect changes in scantable!
    2196         if len(rows) > 0:
    2197             self.blpars=[]
    2198             self.masklists=[]
    2199             self.actualmask=[]
    2200         asaplog.push("Processing:")
    2201         for r in rows:
    2202             msg = " Scan[%d] Beam[%d] IF[%d] Pol[%d] Cycle[%d]" % \
    2203                 (workscan.getscan(r), workscan.getbeam(r), workscan.getif(r), \
    2204                  workscan.getpol(r), workscan.getcycle(r))
    2205             asaplog.push(msg, False)
    2206 
    2207             # figure out edge parameter
    2208             if individualedge:
    2209                 if len(edge) >= workscan.getif(r):
    2210                     raise RuntimeError, "Number of edge elements appear to " \
    2211                                         "be less than the number of IFs"
    2212                     curedge = edge[workscan.getif(r)]
    2213 
    2214             actualmask = mask_and(mask, workscan._getmask(r))    # (CAS-1434)
    2215 
    2216             # setup line finder
    2217             fl.find_lines(r, actualmask, curedge)
    2218            
    2219             f.x = workscan._getabcissa(r)
    2220             f.y = workscan._getspectrum(r)
    2221             f.mask = fl.get_mask()
    2222             f.data = None
    2223             f.fit()
    2224 
    2225             # Show mask list
    2226             masklist=workscan.get_masklist(f.mask, row=r, silent=True)
    2227             msg = "mask range: "+str(masklist)
    2228             asaplog.push(msg, False)
    2229 
    2230             if plot:
    2231                 f.plot(residual=True)
    2232                 x = raw_input("Accept fit ( [y]/n ): ")
    2233                 if x.upper() == 'N':
    2234                     self.blpars.append(None)
    2235                     self.masklists.append(None)
    2236                     self.actualmask.append(None)
    2237                     continue
    2238 
    2239             workscan._setspectrum(f.fitter.getresidual(), r)
    2240             self.blpars.append(f.get_parameters())
    2241             self.masklists.append(masklist)
    2242             self.actualmask.append(f.mask)
    2243         if plot:
    2244             f._p.unmap()
    2245             f._p = None
    2246         workscan._add_history("auto_poly_baseline", varlist)
    2247         if insitu:
    2248             self._assign(workscan)
    2249         else:
    2250             return workscan
    2251 
    22522476    @asaplog_post_dec
    22532477    def rotate_linpolphase(self, angle):
     
    27442968            self.set_freqframe(rcParams['scantable.freqframe'])
    27452969
     2970
    27462971    def __getitem__(self, key):
    27472972        if key < 0:
Note: See TracChangeset for help on using the changeset viewer.