Changeset 2047 for trunk/python


Ignore:
Timestamp:
03/15/11 18:31:04 (14 years ago)
Author:
WataruKawasaki
Message:

New Development: Yes

JIRA Issue: Yes CAS-2847

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: Yes

Module(s): scantable

Description: added {auto_}sinusoid_baseline() for sinusoidal baseline fitting. also minor bug fixes for asapfitter.


Location:
trunk/python
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/python/asapfitter.py

    r1938 r2047  
    7373        Set the function to be fit.
    7474        Parameters:
    75             poly:    use a polynomial of the order given with nonlinear least squares fit
    76             lpoly:   use polynomial of the order given with linear least squares fit
    77             gauss:   fit the number of gaussian specified
    78             lorentz: fit the number of lorentzian specified
     75            poly:     use a polynomial of the order given with nonlinear least squares fit
     76            lpoly:    use polynomial of the order given with linear least squares fit
     77            gauss:    fit the number of gaussian specified
     78            lorentz:  fit the number of lorentzian specified
     79            sinusoid: fit the number of sinusoid specified
    7980        Example:
    8081            fitter.set_function(poly=3)  # will fit a 3rd order polynomial via nonlinear method
    … …  
    8283            fitter.set_function(gauss=2) # will fit two gaussians
    8384            fitter.set_function(lorentz=2) # will fit two lorentzians
     85            fitter.set_function(sinusoid=3) # will fit three sinusoids
    8486        """
    8587        #default poly order 0
    … …  
    109111            self.components = [ 3 for i in range(n) ]
    110112            self.uselinear = False
     113        elif kwargs.has_key('sinusoid'):
     114            n = kwargs.get('sinusoid')
     115            self.fitfunc = 'sinusoid'
     116            self.fitfuncs = [ 'sinusoid' for i in range(n) ]
     117            self.components = [ 3 for i in range(n) ]
     118            self.uselinear = False
    111119        else:
    112120            msg = "Invalid function type."
    … …  
    204212              fixed:     a vector of which parameters are to be held fixed
    205213                         (default is none)
    206               component: in case of multiple gaussians/lorentzians,
    207                          the index of the component
     214              component: in case of multiple gaussians/lorentzians/sinusoidals,
     215                         the index of the target component
    208216        """
    209217        component = None
    … …  
    220228            msg = "Please specify a fitting function first."
    221229            raise RuntimeError(msg)
    222         if (self.fitfunc == "gauss" or self.fitfunc == 'lorentz') and component is not None:
     230        if (self.fitfunc == "gauss" or self.fitfunc == "lorentz" or self.fitfunc == "sinusoid") and component is not None:
    223231            if not self.fitted and sum(self.fitter.getparameters()) == 0:
    224232                pars = _n_bools(len(self.components)*3, False)
    225                 fxd = _n_bools(len(pars), False)
     233                fxd  = _n_bools(len(pars), False)
    226234            else:
    227235                pars = list(self.fitter.getparameters())
    228                 fxd = list(self.fitter.getfixedparameters())
     236                fxd  = list(self.fitter.getfixedparameters())
    229237            i = 3*component
    230238            pars[i:i+3] = params
    231             fxd[i:i+3] = fixed
     239            fxd[i:i+3]  = fixed
    232240            params = pars
    233             fixed = fxd
     241            fixed  = fxd
    234242        self.fitter.setparameters(params)
    235243        if fixed is not None:
    … …  
    290298            d = {'params':[peak, centre, fwhm],
    291299                 'fixed':[peakfixed, centrefixed, fwhmfixed]}
     300            self.set_parameters(d, component)
     301        else:
     302            msg = "Please select a valid  component."
     303            raise ValueError(msg)
     304
     305    @asaplog_post_dec
     306    def set_sinusoid_parameters(self, ampl, period, x0,
     307                             amplfixed=0, periodfixed=0,
     308                             x0fixed=0,
     309                             component=0):
     310        """
     311        Set the Parameters of a 'Sinusoidal' component, set with set_function.
     312        Parameters:
     313            ampl, period, x0:  The sinusoidal parameters
     314            amplfixed,
     315            periodfixed,
     316            x0fixed:             Optional parameters to indicate if
     317                                 the paramters should be held fixed during
     318                                 the fitting process. The default is to keep
     319                                 all parameters flexible.
     320            component:           The number of the component (Default is the
     321                                 component 0)
     322        """
     323        if self.fitfunc != "sinusoid":
     324            msg = "Function only operates on Sinusoidal components."
     325            raise ValueError(msg)
     326        if 0 <= component < len(self.components):
     327            d = {'params':[ampl, period, x0],
     328                 'fixed': [amplfixed, periodfixed, x0fixed]}
    292329            self.set_parameters(d, component)
    293330        else:
    … …  
    338375        cerrs = errs
    339376        if component is not None:
    340             if self.fitfunc == "gauss" or self.fitfunc == "lorentz":
     377            if self.fitfunc == "gauss" or self.fitfunc == "lorentz" or self.fitfunc == "sinusoid":
    341378                i = 3*component
    342379                if i < len(errs):
    … …  
    361398        area = []
    362399        if component is not None:
    363             if self.fitfunc == "gauss" or self.fitfunc == "lorentz":
     400            if self.fitfunc == "poly" or self.fitfunc == "lpoly":
     401                cpars = pars
     402                cfixed = fixed
     403                cerrs = errs
     404            else:
    364405                i = 3*component
    365406                cpars = pars[i:i+3]
    366407                cfixed = fixed[i:i+3]
    367408                cerrs = errs[i:i+3]
    368                 a = self.get_area(component)
    369                 area = [a for i in range(3)]
    370             else:
    371                 cpars = pars
    372                 cfixed = fixed
    373                 cerrs = errs
     409                if self.fitfunc == "gauss" or self.fitfunc == "lorentz":
     410                    a = self.get_area(component)
     411                    area = [a for i in range(3)]
    374412        else:
    375413            cpars = pars
    … …  
    378416            if self.fitfunc == "gauss" or self.fitfunc == "lorentz":
    379417                for c in range(len(self.components)):
    380                   a = self.get_area(c)
    381                   area += [a for i in range(3)]
     418                    a = self.get_area(c)
     419                    area += [a for i in range(3)]
    382420        fpars = self._format_pars(cpars, cfixed, errors and cerrs, area)
    383421        asaplog.push(fpars)
    … …  
    387425    def _format_pars(self, pars, fixed, errors, area):
    388426        out = ''
    389         if self.fitfunc == 'poly':
     427        if self.fitfunc == "poly" or self.fitfunc == "lpoly":
    390428            c = 0
    391429            for i in range(len(pars)):
    392430                fix = ""
    393431                if len(fixed) and fixed[i]: fix = "(fixed)"
    394                 if errors :
    395                     out += '  p%d%s= %3.6f (%1.6f),' % (c,fix,pars[i], errors[i])
    396                 else:
    397                     out += '  p%d%s= %3.6f,' % (c,fix,pars[i])
     432                out += "  p%d%s= %3.6f" % (c, fix, pars[i])
     433                if errors : out += " (%1.6f)" % errors[i]
     434                out += ","
    398435                c+=1
    399436            out = out[:-1]  # remove trailing ','
    400         elif self.fitfunc == 'gauss' or self.fitfunc == 'lorentz':
     437        else:
    401438            i = 0
    402439            c = 0
    403             aunit = ''
    404             ounit = ''
     440            if self.fitfunc == "gauss" or self.fitfunc == "lorentz":
     441                pnam = ["peak", "centre", "FWHM"]
     442            elif self.fitfunc == "sinusoid":
     443                pnam = ["amplitude", "period", "x0"]
     444            aunit = ""
     445            ounit = ""
    405446            if self.data:
    406447                aunit = self.data.get_unit()
    407448                ounit = self.data.get_fluxunit()
    408449            while i < len(pars):
    409                 if len(area):
    410                     out += '  %2d: peak = %3.3f %s , centre = %3.3f %s, FWHM = %3.3f %s\n      area = %3.3f %s %s\n' % (c,pars[i],ounit,pars[i+1],aunit,pars[i+2],aunit, area[i],ounit,aunit)
    411                 else:
    412                     out += '  %2d: peak = %3.3f %s , centre = %3.3f %s, FWHM = %3.3f %s\n' % (c,pars[i],ounit,pars[i+1],aunit,pars[i+2],aunit,ounit,aunit)
     450                fix0 = fix1 = fix2 = ""
     451                if i < len(fixed)-2:
     452                    if fixed[i]:   fix0 = "(fixed)"
     453                    if fixed[i+1]: fix1 = "(fixed)"
     454                    if fixed[i+2]: fix2 = "(fixed)"
     455                out += "  %2d: " % c
     456                out += "%s%s = %3.3f %s, " % (pnam[0], fix0, pars[i],   ounit)
     457                out += "%s%s = %3.3f %s, " % (pnam[1], fix1, pars[i+1], aunit)
     458                out += "%s%s = %3.3f %s\n" % (pnam[2], fix2, pars[i+2], aunit)
     459                if len(area): out += "      area = %3.3f %s %s\n" % (area[i], ounit, aunit)
    413460                c+=1
    414461                i+=3

  • trunk/python/scantable.py

    r2029 r2047  
    20702070
    20712071    @asaplog_post_dec
     2072    def sinusoid_baseline(self, insitu=None, mask=None, minnwave=None, maxnwave=None,
     2073                          clipthresh=None, clipniter=None, plot=None, outlog=None, blfile=None):
     2074        """\
     2075        Return a scan which has been baselined (all rows) by cubic spline function (piecewise cubic polynomial).
     2076        Parameters:
     2077            insitu:     If False a new scantable is returned.
     2078                        Otherwise, the scaling is done in-situ
     2079                        The default is taken from .asaprc (False)
     2080            mask:       An optional mask
     2081            minnwave:   Minimum wave number in spectral window (default is 0)
     2082            maxnwave:   Maximum wave number in spectral window (default is 3)
     2083            clipthresh: Clipping threshold. (default is 3.0, unit: sigma)
     2084            clipniter:  maximum number of iteration of 'clipthresh'-sigma clipping (default is 1)
     2085            plot:   *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
     2086                        plot the fit and the residual. In this each
     2087                        indivual fit has to be approved, by typing 'y'
     2088                        or 'n'
     2089            outlog:     Output the coefficients of the best-fit
     2090                        function to logger (default is False)
     2091            blfile:     Name of a text file in which the best-fit
     2092                        parameter values to be written
     2093                        (default is "": no file/logger output)
     2094
     2095        Example:
     2096            # return a scan baselined by a combination of sinusoidal curves having
     2097            # wave numbers in spectral window from 1 to 10,
     2098            # also with 3-sigma clipping, iteration up to 4 times
     2099            bscan = scan.sinusoid_baseline(maxnwave=10,clipthresh=3.0,clipniter=4)
     2100        """
     2101       
     2102        varlist = vars()
     2103       
     2104        if insitu is None: insitu = rcParams["insitu"]
     2105        if insitu:
     2106            workscan = self
     2107        else:
     2108            workscan = self.copy()
     2109
     2110        nchan = workscan.nchan()
     2111       
     2112        if mask is None: mask = [True for i in xrange(nchan)]
     2113        if minnwave is None: minnwave = 0
     2114        if maxnwave is None: maxnwave = 3
     2115        if clipthresh is None: clipthresh = 3.0
     2116        if clipniter is None: clipniter = 1
     2117        if plot is None: plot = False
     2118        if outlog is None: outlog = False
     2119        if blfile is None: blfile = ""
     2120
     2121        outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
     2122       
     2123        try:
     2124            #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method.
     2125            workscan._sinusoid_baseline(mask, minnwave, maxnwave, clipthresh, clipniter, outlog, blfile)
     2126           
     2127            workscan._add_history("sinusoid_baseline", varlist)
     2128           
     2129            if insitu:
     2130                self._assign(workscan)
     2131            else:
     2132                return workscan
     2133           
     2134        except RuntimeError, e:
     2135            msg = "The fit failed, possibly because it didn't converge."
     2136            if rcParams["verbose"]:
     2137                asaplog.push(str(e))
     2138                asaplog.push(str(msg))
     2139                return
     2140            else:
     2141                raise RuntimeError(str(e)+'\n'+msg)
     2142
     2143
     2144    def auto_sinusoid_baseline(self, insitu=None, mask=None, minnwave=None, maxnwave=None,
     2145                               clipthresh=None, clipniter=None, edge=None, threshold=None,
     2146                               chan_avg_limit=None, plot=None, outlog=None, blfile=None):
     2147        """\
     2148        Return a scan which has been baselined (all rows) by cubic spline
     2149        function (piecewise cubic polynomial).
     2150        Spectral lines are detected first using linefinder and masked out
     2151        to avoid them affecting the baseline solution.
     2152
     2153        Parameters:
     2154            insitu:     if False a new scantable is returned.
     2155                        Otherwise, the scaling is done in-situ
     2156                        The default is taken from .asaprc (False)
     2157            mask:       an optional mask retreived from scantable
     2158            minnwave:   Minimum wave number in spectral window (default is 0)
     2159            maxnwave:   Maximum wave number in spectral window (default is 3)
     2160            clipthresh: Clipping threshold. (default is 3.0, unit: sigma)
     2161            clipniter:  maximum number of iteration of 'clipthresh'-sigma clipping (default is 1)
     2162            edge:       an optional number of channel to drop at
     2163                        the edge of spectrum. If only one value is
     2164                        specified, the same number will be dropped
     2165                        from both sides of the spectrum. Default
     2166                        is to keep all channels. Nested tuples
     2167                        represent individual edge selection for
     2168                        different IFs (a number of spectral channels
     2169                        can be different)
     2170            threshold:  the threshold used by line finder. It is
     2171                        better to keep it large as only strong lines
     2172                        affect the baseline solution.
     2173            chan_avg_limit:
     2174                        a maximum number of consequtive spectral
     2175                        channels to average during the search of
     2176                        weak and broad lines. The default is no
     2177                        averaging (and no search for weak lines).
     2178                        If such lines can affect the fitted baseline
     2179                        (e.g. a high order polynomial is fitted),
     2180                        increase this parameter (usually values up
     2181                        to 8 are reasonable). Most users of this
     2182                        method should find the default value sufficient.
     2183            plot:   *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
     2184                        plot the fit and the residual. In this each
     2185                        indivual fit has to be approved, by typing 'y'
     2186                        or 'n'
     2187            outlog:     Output the coefficients of the best-fit
     2188                        function to logger (default is False)
     2189            blfile:     Name of a text file in which the best-fit
     2190                        parameter values to be written
     2191                        (default is "": no file/logger output)
     2192
     2193        Example:
     2194            bscan = scan.auto_sinusoid_baseline(maxnwave=10, insitu=False)
     2195        """
     2196
     2197        varlist = vars()
     2198
     2199        if insitu is None: insitu = rcParams['insitu']
     2200        if insitu:
     2201            workscan = self
     2202        else:
     2203            workscan = self.copy()
     2204
     2205        nchan = workscan.nchan()
     2206       
     2207        if mask is None: mask = [True for i in xrange(nchan)]
     2208        if minnwave is None: minnwave = 0
     2209        if maxnwave is None: maxnwave = 3
     2210        if clipthresh is None: clipthresh = 3.0
     2211        if clipniter is None: clipniter = 1
     2212        if edge is None: edge = (0, 0)
     2213        if threshold is None: threshold = 3
     2214        if chan_avg_limit is None: chan_avg_limit = 1
     2215        if plot is None: plot = False
     2216        if outlog is None: outlog = False
     2217        if blfile is None: blfile = ""
     2218
     2219        outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
     2220       
     2221        from asap.asaplinefind import linefinder
     2222        from asap import _is_sequence_or_number as _is_valid
     2223
     2224        if not (isinstance(edge, list) or isinstance(edge, tuple)): edge = [ edge ]
     2225        individualedge = False;
     2226        if len(edge) > 1: individualedge = isinstance(edge[0], list) or isinstance(edge[0], tuple)
     2227
     2228        if individualedge:
     2229            for edgepar in edge:
     2230                if not _is_valid(edgepar, int):
     2231                    raise ValueError, "Each element of the 'edge' tuple has \
     2232                                       to be a pair of integers or an integer."
     2233        else:
     2234            if not _is_valid(edge, int):
     2235                raise ValueError, "Parameter 'edge' has to be an integer or a \
     2236                                   pair of integers specified as a tuple. \
     2237                                   Nested tuples are allowed \
     2238                                   to make individual selection for different IFs."
     2239
     2240            if len(edge) > 1:
     2241                curedge = edge
     2242            else:
     2243                curedge = edge + edge
     2244
     2245        try:
     2246            #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method.
     2247            if individualedge:
     2248                curedge = []
     2249                for i in xrange(len(edge)):
     2250                    curedge += edge[i]
     2251               
     2252            workscan._auto_sinusoid_baseline(mask, minnwave, maxnwave, clipthresh, clipniter, curedge, threshold, chan_avg_limit, outlog, blfile)
     2253
     2254            workscan._add_history("auto_sinusoid_baseline", varlist)
     2255           
     2256            if insitu:
     2257                self._assign(workscan)
     2258            else:
     2259                return workscan
     2260           
     2261        except RuntimeError, e:
     2262            msg = "The fit failed, possibly because it didn't converge."
     2263            if rcParams["verbose"]:
     2264                asaplog.push(str(e))
     2265                asaplog.push(str(msg))
     2266                return
     2267            else:
     2268                raise RuntimeError(str(e)+'\n'+msg)
     2269
     2270
     2271    @asaplog_post_dec
    20722272    def cspline_baseline(self, insitu=None, mask=None, npiece=None, clipthresh=None, clipniter=None, plot=None, outlog=None, blfile=None):
    20732273        """\
Note: See TracChangeset for help on using the changeset viewer.