Ignore:
Timestamp:
06/07/11 23:49:53 (13 years ago)
Author:
WataruKawasaki
Message:

New Development: Yes

JIRA Issue: Yes CAS-3149

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: scantable.*sinusoid_baseline() params

Test Programs:

Put in Release Notes: Yes

Module(s):

Description: (1) Implemented an automated sinusoidal fitting functionality

(2) FFT available with scantable.fft()
(3) fixed a bug of parsing 'edge' param used by linefinder.
(4) a function to show progress status for row-based iterations.


File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/python/scantable.py

    r2178 r2186  
    7676    else:
    7777        return False
     78
     79def normalise_edge_param(edge):
     80    """\
     81    Convert a given edge value to a one-dimensional array that can be
     82    given to baseline-fitting/subtraction functions.
     83    The length of the output value will be an even because values for
     84    the both sides of spectra are to be contained for each IF. When
     85    the length is 2, the values will be applied to all IFs. If the length
     86    is larger than 2, it will be 2*ifnos().
     87    Accepted format of edge include:
     88            * an integer - will be used for both sides of spectra of all IFs.
     89                  e.g. 10 is converted to [10,10]
     90            * an empty list/tuple [] - converted to [0, 0] and used for all IFs.
     91            * a list/tuple containing an integer - same as the above case.
     92                  e.g. [10] is converted to [10,10]
     93            * a list/tuple containing two integers - will be used for all IFs.
     94                  e.g. [5,10] is output as it is. no need to convert.
     95            * a list/tuple of lists/tuples containing TWO integers -
     96                  each element of edge will be used for each IF.
     97                  e.g. [[5,10],[15,20]] - [5,10] for IF[0] and [15,20] for IF[1].
     98                 
     99                  If an element contains the same integer values, the input 'edge'
     100                  parameter can be given in a simpler shape in the following cases:
     101                      ** when len(edge)!=2
     102                          any elements containing the same values can be replaced
     103                          to single integers.
     104                          e.g. [[15,15]] can be simplified to [15] (or [15,15] or 15 also).
     105                          e.g. [[1,1],[2,2],[3,3]] can be simplified to [1,2,3].
     106                      ** when len(edge)=2
     107                          care is needed for this case: ONLY ONE of the
     108                          elements can be a single integer,
     109                          e.g. [[5,5],[10,10]] can be simplified to [5,[10,10]]
     110                               or [[5,5],10], but can NOT be simplified to [5,10].
     111                               when [5,10] given, it is interpreted as
     112                               [[5,10],[5,10],[5,10],...] instead, as shown before.
     113    """
     114    from asap import _is_sequence_or_number as _is_valid
     115    if isinstance(edge, list) or isinstance(edge, tuple):
     116        for edgepar in edge:
     117            if not _is_valid(edgepar, int):
     118                raise ValueError, "Each element of the 'edge' tuple has \
     119                                   to be a pair of integers or an integer."
     120            if isinstance(edgepar, list) or isinstance(edgepar, tuple):
     121                if len(edgepar) != 2:
     122                    raise ValueError, "Each element of the 'edge' tuple has \
     123                                       to be a pair of integers or an integer."
     124    else:
     125        if not _is_valid(edge, int):
     126            raise ValueError, "Parameter 'edge' has to be an integer or a \
     127                               pair of integers specified as a tuple. \
     128                               Nested tuples are allowed \
     129                               to make individual selection for different IFs."
     130       
     131
     132    if isinstance(edge, int):
     133        edge = [ edge, edge ]                 # e.g. 3   => [3,3]
     134    elif isinstance(edge, list) or isinstance(edge, tuple):
     135        if len(edge) == 0:
     136            edge = [0, 0]                     # e.g. []  => [0,0]
     137        elif len(edge) == 1:
     138            if isinstance(edge[0], int):
     139                edge = [ edge[0], edge[0] ]   # e.g. [1] => [1,1]
     140
     141    commonedge = True
     142    if len(edge) > 2: commonedge = False
     143    else:
     144        for edgepar in edge:
     145            if isinstance(edgepar, list) or isinstance(edgepar, tuple):
     146                commonedge = False
     147                break
     148
     149    if commonedge:
     150        if len(edge) > 1:
     151            norm_edge = edge
     152        else:
     153            norm_edge = edge + edge
     154    else:
     155        norm_edge = []
     156        for edgepar in edge:
     157            if isinstance(edgepar, int):
     158                norm_edge += [edgepar, edgepar]
     159            else:
     160                norm_edge += edgepar
     161
     162    return norm_edge
     163
     164def raise_fitting_failure_exception(e):
     165    msg = "The fit failed, possibly because it didn't converge."
     166    if rcParams["verbose"]:
     167        asaplog.push(str(e))
     168        asaplog.push(str(msg))
     169    else:
     170        raise RuntimeError(str(e)+'\n'+msg)
    78171   
     172
    79173class scantable(Scantable):
    80174    """\
    … …  
    11201214
    11211215    @asaplog_post_dec
    1122     def fft(self, getrealimag=False, rowno=[]):
     1216    def fft(self, rowno=[], mask=[], getrealimag=False):
    11231217        """\
    11241218        Apply FFT to the spectra.
    … …  
    11261220
    11271221        Parameters:
     1222
     1223            rowno:          The row number(s) to be processed. int, list
     1224                            and tuple are accepted. By default ([]), FFT
     1225                            is applied to the whole data.
     1226
     1227            mask:           Auxiliary channel mask(s). Given as a boolean
     1228                            list, it is applied to all specified rows.
     1229                            A list of boolean lists can also be used to
     1230                            apply different masks. In the latter case, the
     1231                            length of 'mask' must be the same as that of
     1232                            'rowno'. The default is [].
    11281233       
    11291234            getrealimag:    If True, returns the real and imaginary part
    … …  
    11331238                            phase (argument, in unit of radian).
    11341239
    1135             rowno:          The row number(s) to be processed. int, list
    1136                             and tuple are accepted. By default ([]), FFT
    1137                             is applied to the whole data.
    1138 
    11391240        Returns:
    11401241
    1141             A dictionary containing two keys, i.e., 'real' and 'imag' for
    1142             getrealimag = True, or 'ampl' and 'phase' for getrealimag = False,
     1242            A list of dictionaries containing the results for each spectrum.
     1243            Each dictionary contains two values, the real and the imaginary
     1244            parts when getrealimag = True, or the amplitude(absolute value)
     1245            and the phase(argument) when getrealimag = False. The key for
     1246            these values are 'real' and 'imag', or 'ampl' and 'phase',
    11431247            respectively.
    1144             The value for each key is a list of lists containing the FFT
    1145             results from each spectrum.
    1146            
    1147         """
    1148         if getrealimag:
    1149             res = {"real":[], "imag":[]}  # return real and imaginary values
    1150         else:
    1151             res = {"ampl":[], "phase":[]} # return amplitude and phase(argument)
    1152        
     1248        """
    11531249        if isinstance(rowno, int):
    11541250            rowno = [rowno]
    11551251        elif not (isinstance(rowno, list) or isinstance(rowno, tuple)):
    1156             raise TypeError("The row number(s) should be int, list or tuple.")
     1252            raise TypeError("The row number(s) must be int, list or tuple.")
     1253
     1254        if len(rowno) == 0: rowno = [i for i in xrange(self.nrow())]
     1255
     1256        if not (isinstance(mask, list) or isinstance(mask, tuple)):
     1257            raise TypeError("The mask must be a boolean list or a list of boolean list.")
     1258        if len(mask) == 0: mask = [True for i in xrange(self.nchan())]
     1259        if isinstance(mask[0], bool): mask = [mask]
     1260        elif not (isinstance(mask[0], list) or isinstance(mask[0], tuple)):
     1261            raise TypeError("The mask must be a boolean list or a list of boolean list.")
     1262
     1263        usecommonmask = (len(mask) == 1)
     1264        if not usecommonmask:
     1265            if len(mask) != len(rowno):
     1266                raise ValueError("When specifying masks for each spectrum, the numbers of them must be identical.")
     1267        for amask in mask:
     1268            if len(amask) != self.nchan():
     1269                raise ValueError("The spectra and the mask have different length.")
    11571270       
    1158         nrow = len(rowno)
    1159         if nrow == 0: nrow = self.nrow()  # by default, calculate for all rows.
    1160        
    1161         fspec = self._math._fft(self, rowno, getrealimag)
    1162         nspec = len(fspec)/nrow
    1163        
    1164         i = 0
    1165         while (i < nrow):
    1166             v1 = []
    1167             v2 = []
     1271        res = []
     1272
     1273        imask = 0
     1274        for whichrow in rowno:
     1275            fspec = self._fft(whichrow, mask[imask], getrealimag)
     1276            nspec = len(fspec)
    11681277           
    1169             j = 0
    1170             while (j < nspec):
    1171                 k = i*nspec + j
    1172                 v1.append(fspec[k])
    1173                 v2.append(fspec[k+1])
    1174                 j += 2
     1278            i=0
     1279            v1=[]
     1280            v2=[]
     1281            reselem = {"real":[],"imag":[]} if getrealimag else {"ampl":[],"phase":[]}
     1282           
     1283            while (i < nspec):
     1284                v1.append(fspec[i])
     1285                v2.append(fspec[i+1])
     1286                i+=2
    11751287           
    11761288            if getrealimag:
    1177                 res["real"].append(v1)
    1178                 res["imag"].append(v2)
     1289                reselem["real"]  += v1
     1290                reselem["imag"]  += v2
    11791291            else:
    1180                 res["ampl"].append(v1)
    1181                 res["phase"].append(v2)
     1292                reselem["ampl"]  += v1
     1293                reselem["phase"] += v2
    11821294           
    1183             i += 1
    1184 
     1295            res.append(reselem)
     1296           
     1297            if not usecommonmask: imask += 1
     1298       
    11851299        return res
    11861300
    … …  
    21312245        else: return s
    21322246
    2133 
    2134     @asaplog_post_dec
    2135     def sinusoid_baseline(self, insitu=None, mask=None, nwave=None, maxwavelength=None,
    2136                           clipthresh=None, clipniter=None, plot=None, getresidual=None, outlog=None, blfile=None):
     2247    @asaplog_post_dec
     2248    def _parse_wn(self, wn):
     2249        if isinstance(wn, list) or isinstance(wn, tuple):
     2250            return wn
     2251        elif isinstance(wn, int):
     2252            return [ wn ]
     2253        elif isinstance(wn, str):
     2254            if '-' in wn:                            # case 'a-b' : return [a,a+1,...,b-1,b]
     2255                val = wn.split('-')
     2256                val = [int(val[0]), int(val[1])]
     2257                val.sort()
     2258                res = [i for i in xrange(val[0], val[1]+1)]
     2259            elif wn[:2] == '<=' or wn[:2] == '=<':   # cases '<=a','=<a' : return [0,1,...,a-1,a]
     2260                val = int(wn[2:])+1
     2261                res = [i for i in xrange(val)]
     2262            elif wn[-2:] == '>=' or wn[-2:] == '=>': # cases 'a>=','a=>' : return [0,1,...,a-1,a]
     2263                val = int(wn[:-2])+1
     2264                res = [i for i in xrange(val)]
     2265            elif wn[0] == '<':                       # case '<a' :         return [0,1,...,a-2,a-1]
     2266                val = int(wn[1:])
     2267                res = [i for i in xrange(val)]
     2268            elif wn[-1] == '>':                      # case 'a>' :         return [0,1,...,a-2,a-1]
     2269                val = int(wn[:-1])
     2270                res = [i for i in xrange(val)]
     2271            elif wn[:2] == '>=' or wn[:2] == '=>':   # cases '>=a','=>a' : return [a,a+1,...,a_nyq]
     2272                val = int(wn[2:])
     2273                res = [i for i in xrange(val, self.nchan()/2+1)]
     2274            elif wn[-2:] == '<=' or wn[-2:] == '=<': # cases 'a<=','a=<' : return [a,a+1,...,a_nyq]
     2275                val = int(wn[:-2])
     2276                res = [i for i in xrange(val, self.nchan()/2+1)]
     2277            elif wn[0] == '>':                       # case '>a' :         return [a+1,a+2,...,a_nyq]
     2278                val = int(wn[1:])+1
     2279                res = [i for i in xrange(val, self.nchan()/2+1)]
     2280            elif wn[-1] == '<':                      # case 'a<' :         return [a+1,a+2,...,a_nyq]
     2281                val = int(wn[:-1])+1
     2282                res = [i for i in xrange(val, self.nchan()/2+1)]
     2283
     2284            return res
     2285        else:
     2286            msg = 'wrong value given for addwn/rejwn'
     2287            raise RuntimeError(msg)
     2288
     2289
     2290    @asaplog_post_dec
     2291    def sinusoid_baseline(self, insitu=None, mask=None, applyfft=None, fftmethod=None, fftthresh=None,
     2292                          addwn=None, rejwn=None, clipthresh=None, clipniter=None, plot=None,
     2293                          getresidual=None, outlog=None, blfile=None):
    21372294        """\
    21382295        Return a scan which has been baselined (all rows) with sinusoidal functions.
    21392296        Parameters:
    2140             insitu:        If False a new scantable is returned.
     2297            insitu:        if False a new scantable is returned.
    21412298                           Otherwise, the scaling is done in-situ
    21422299                           The default is taken from .asaprc (False)
    2143             mask:          An optional mask
    2144             nwave:         the maximum wave number of sinusoids within
    2145                            maxwavelength * (spectral range).
    2146                            The default is 3 (i.e., sinusoids with wave
    2147                            number of 0(=constant), 1, 2, and 3 are
    2148                            used for fitting). Also it is possible to
    2149                            explicitly specify all the wave numbers to
    2150                            be used, by giving a list including them
    2151                            (e.g. [0,1,2,15,16]).
    2152             maxwavelength: the longest sinusoidal wavelength. The
    2153                            default is 1.0 (unit: spectral range).
     2300            mask:          an optional mask
     2301            applyfft:      if True use some method, such as FFT, to find
     2302                           strongest sinusoidal components in the wavenumber
     2303                           domain to be used for baseline fitting.
     2304                           default is True.
     2305            fftmethod:     method to find the strong sinusoidal components.
     2306                           now only 'fft' is available and it is the default.
     2307            fftthresh:     the threshold to select wave numbers to be used for
     2308                           fitting from the distribution of amplitudes in the
     2309                           wavenumber domain.
     2310                           both float and string values accepted.
     2311                           given a float value, the unit is set to sigma.
     2312                           for string values, allowed formats include:
     2313                               'xsigma' or 'x' (= x-sigma level. e.g., '3sigma'), or
     2314                               'topx' (= the x strongest ones, e.g. 'top5').
     2315                           default is 3.0 (unit: sigma).
     2316            addwn:         the additional wave numbers to be used for fitting.
     2317                           list or integer value is accepted to specify every
     2318                           wave numbers. also string value can be used in case
     2319                           you need to specify wave numbers in a certain range,
     2320                           e.g., 'a-b' (= a, a+1, a+2, ..., b-1, b),
     2321                                 '<a'  (= 0,1,...,a-2,a-1),
     2322                                 '>=a' (= a, a+1, ... up to the maximum wave
     2323                                        number corresponding to the Nyquist
     2324                                        frequency for the case of FFT).
     2325                           default is [].
     2326            rejwn:         the wave numbers NOT to be used for fitting.
     2327                           can be set just as addwn but has higher priority:
     2328                           wave numbers which are specified both in addwn
     2329                           and rejwn will NOT be used. default is [].
    21542330            clipthresh:    Clipping threshold. (default is 3.0, unit: sigma)
    21552331            clipniter:     maximum number of iteration of 'clipthresh'-sigma clipping (default is 0)
    … …  
    21642340            blfile:        Name of a text file in which the best-fit
    21652341                           parameter values to be written
    2166                            (default is "": no file/logger output)
     2342                           (default is '': no file/logger output)
    21672343
    21682344        Example:
    … …  
    21702346            # wave numbers in spectral window up to 10,
    21712347            # also with 3-sigma clipping, iteration up to 4 times
    2172             bscan = scan.sinusoid_baseline(nwave=10,clipthresh=3.0,clipniter=4)
     2348            bscan = scan.sinusoid_baseline(addwn='<=10',clipthresh=3.0,clipniter=4)
    21732349
    21742350        Note:
    … …  
    21772353        """
    21782354       
    2179         varlist = vars()
     2355        try:
     2356            varlist = vars()
    21802357       
    2181         if insitu is None: insitu = rcParams["insitu"]
    2182         if insitu:
    2183             workscan = self
    2184         else:
    2185             workscan = self.copy()
    2186 
    2187         nchan = workscan.nchan()
    2188        
    2189         if mask          is None: mask          = [True for i in xrange(nchan)]
    2190         if nwave         is None: nwave         = 3
    2191         if maxwavelength is None: maxwavelength = 1.0
    2192         if clipthresh    is None: clipthresh    = 3.0
    2193         if clipniter     is None: clipniter     = 0
    2194         if plot          is None: plot          = False
    2195         if getresidual   is None: getresidual   = True
    2196         if outlog        is None: outlog        = False
    2197         if blfile        is None: blfile        = ""
    2198 
    2199         if isinstance(nwave, int):
    2200             in_nwave = nwave
    2201             nwave = []
    2202             for i in xrange(in_nwave+1): nwave.append(i)
    2203        
    2204         outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
    2205        
    2206         try:
     2358            if insitu is None: insitu = rcParams['insitu']
     2359            if insitu:
     2360                workscan = self
     2361            else:
     2362                workscan = self.copy()
     2363           
     2364            if mask          is None: mask          = [True for i in xrange(workscan.nchan())]
     2365            if applyfft      is None: applyfft      = True
     2366            if fftmethod     is None: fftmethod     = 'fft'
     2367            if fftthresh     is None: fftthresh     = 3.0
     2368            if addwn         is None: addwn         = []
     2369            if rejwn         is None: rejwn         = []
     2370            if clipthresh    is None: clipthresh    = 3.0
     2371            if clipniter     is None: clipniter     = 0
     2372            if plot          is None: plot          = False
     2373            if getresidual   is None: getresidual   = True
     2374            if outlog        is None: outlog        = False
     2375            if blfile        is None: blfile        = ''
     2376
    22072377            #CURRENTLY, PLOT=true is UNAVAILABLE UNTIL sinusoidal fitting is implemented as a fitter method.
    2208             workscan._sinusoid_baseline(mask, nwave, maxwavelength, clipthresh, clipniter, getresidual, outlog, blfile)
    2209            
    2210             workscan._add_history("sinusoid_baseline", varlist)
     2378            workscan._sinusoid_baseline(mask, applyfft, fftmethod.lower(), str(fftthresh).lower(), workscan._parse_wn(addwn), workscan._parse_wn(rejwn), clipthresh, clipniter, getresidual, outlog, blfile)
     2379            workscan._add_history('sinusoid_baseline', varlist)
    22112380           
    22122381            if insitu:
    … …  
    22162385           
    22172386        except RuntimeError, e:
    2218             msg = "The fit failed, possibly because it didn't converge."
    2219             if rcParams["verbose"]:
    2220                 asaplog.push(str(e))
    2221                 asaplog.push(str(msg))
    2222                 return
    2223             else:
    2224                 raise RuntimeError(str(e)+'\n'+msg)
    2225 
    2226 
    2227     def auto_sinusoid_baseline(self, insitu=None, mask=None, nwave=None, maxwavelength=None,
    2228                                clipthresh=None, clipniter=None, edge=None, threshold=None,
     2387            raise_fitting_failure_exception(e)
     2388
     2389
     2390    @asaplog_post_dec
     2391    def auto_sinusoid_baseline(self, insitu=None, mask=None, applyfft=None, fftmethod=None, fftthresh=None,
     2392                               addwn=None, rejwn=None, clipthresh=None, clipniter=None, edge=None, threshold=None,
    22292393                               chan_avg_limit=None, plot=None, getresidual=None, outlog=None, blfile=None):
    22302394        """\
    … …  
    22382402                           The default is taken from .asaprc (False)
    22392403            mask:          an optional mask retreived from scantable
    2240             nwave:         the maximum wave number of sinusoids within
    2241                            maxwavelength * (spectral range).
    2242                            The default is 3 (i.e., sinusoids with wave
    2243                            number of 0(=constant), 1, 2, and 3 are
    2244                            used for fitting). Also it is possible to
    2245                            explicitly specify all the wave numbers to
    2246                            be used, by giving a list including them
    2247                            (e.g. [0,1,2,15,16]).
    2248             maxwavelength: the longest sinusoidal wavelength. The
    2249                            default is 1.0 (unit: spectral range).
     2404            applyfft:      if True use some method, such as FFT, to find
     2405                           strongest sinusoidal components in the wavenumber
     2406                           domain to be used for baseline fitting.
     2407                           default is True.
     2408            fftmethod:     method to find the strong sinusoidal components.
     2409                           now only 'fft' is available and it is the default.
     2410            fftthresh:     the threshold to select wave numbers to be used for
     2411                           fitting from the distribution of amplitudes in the
     2412                           wavenumber domain.
     2413                           both float and string values accepted.
     2414                           given a float value, the unit is set to sigma.
     2415                           for string values, allowed formats include:
     2416                               'xsigma' or 'x' (= x-sigma level. e.g., '3sigma'), or
     2417                               'topx' (= the x strongest ones, e.g. 'top5').
     2418                           default is 3.0 (unit: sigma).
     2419            addwn:         the additional wave numbers to be used for fitting.
     2420                           list or integer value is accepted to specify every
     2421                           wave numbers. also string value can be used in case
     2422                           you need to specify wave numbers in a certain range,
     2423                           e.g., 'a-b' (= a, a+1, a+2, ..., b-1, b),
     2424                                 '<a'  (= 0,1,...,a-2,a-1),
     2425                                 '>=a' (= a, a+1, ... up to the maximum wave
     2426                                        number corresponding to the Nyquist
     2427                                        frequency for the case of FFT).
     2428                           default is [].
     2429            rejwn:         the wave numbers NOT to be used for fitting.
     2430                           can be set just as addwn but has higher priority:
     2431                           wave numbers which are specified both in addwn
     2432                           and rejwn will NOT be used. default is [].
    22502433            clipthresh:    Clipping threshold. (default is 3.0, unit: sigma)
    22512434            clipniter:     maximum number of iteration of 'clipthresh'-sigma clipping (default is 0)
    … …  
    22832466
    22842467        Example:
    2285             bscan = scan.auto_sinusoid_baseline(nwave=10, insitu=False)
     2468            bscan = scan.auto_sinusoid_baseline(addwn='<=10', insitu=False)
    22862469       
    22872470        Note:
    … …  
    22902473        """
    22912474
    2292         varlist = vars()
    2293 
    2294         if insitu is None: insitu = rcParams['insitu']
    2295         if insitu:
    2296             workscan = self
    2297         else:
    2298             workscan = self.copy()
    2299 
    2300         nchan = workscan.nchan()
    2301        
    2302         if mask           is None: mask           = [True for i in xrange(nchan)]
    2303         if nwave          is None: nwave          = 3
    2304         if maxwavelength  is None: maxwavelength  = 1.0
    2305         if clipthresh     is None: clipthresh     = 3.0
    2306         if clipniter      is None: clipniter      = 0
    2307         if edge           is None: edge           = (0,0)
    2308         if threshold      is None: threshold      = 3
    2309         if chan_avg_limit is None: chan_avg_limit = 1
    2310         if plot           is None: plot           = False
    2311         if getresidual    is None: getresidual    = True
    2312         if outlog         is None: outlog         = False
    2313         if blfile         is None: blfile         = ""
    2314 
    2315         outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
    2316        
    2317         from asap.asaplinefind import linefinder
    2318         from asap import _is_sequence_or_number as _is_valid
    2319 
    2320         if isinstance(nwave, int):
    2321             in_nwave = nwave
    2322             nwave = []
    2323             for i in xrange(in_nwave+1): nwave.append(i)
    2324 
    2325         if not (isinstance(edge, list) or isinstance(edge, tuple)): edge = [ edge ]
    2326         individualedge = False;
    2327         if len(edge) > 1: individualedge = isinstance(edge[0], list) or isinstance(edge[0], tuple)
    2328 
    2329         if individualedge:
    2330             for edgepar in edge:
    2331                 if not _is_valid(edgepar, int):
    2332                     raise ValueError, "Each element of the 'edge' tuple has \
    2333                                        to be a pair of integers or an integer."
    2334         else:
    2335             if not _is_valid(edge, int):
    2336                 raise ValueError, "Parameter 'edge' has to be an integer or a \
    2337                                    pair of integers specified as a tuple. \
    2338                                    Nested tuples are allowed \
    2339                                    to make individual selection for different IFs."
    2340 
    2341             if len(edge) > 1:
    2342                 curedge = edge
     2475        try:
     2476            varlist = vars()
     2477
     2478            if insitu is None: insitu = rcParams['insitu']
     2479            if insitu:
     2480                workscan = self
    23432481            else:
    2344                 curedge = edge + edge
    2345 
    2346         try:
     2482                workscan = self.copy()
     2483           
     2484            if mask           is None: mask           = [True for i in xrange(workscan.nchan())]
     2485            if applyfft       is None: applyfft       = True
     2486            if fftmethod      is None: fftmethod      = 'fft'
     2487            if fftthresh      is None: fftthresh      = 3.0
     2488            if addwn          is None: addwn          = []
     2489            if rejwn          is None: rejwn          = []
     2490            if clipthresh     is None: clipthresh     = 3.0
     2491            if clipniter      is None: clipniter      = 0
     2492            if edge           is None: edge           = (0,0)
     2493            if threshold      is None: threshold      = 3
     2494            if chan_avg_limit is None: chan_avg_limit = 1
     2495            if plot           is None: plot           = False
     2496            if getresidual    is None: getresidual    = True
     2497            if outlog         is None: outlog         = False
     2498            if blfile         is None: blfile         = ''
     2499
    23472500            #CURRENTLY, PLOT=true is UNAVAILABLE UNTIL sinusoidal fitting is implemented as a fitter method.
    2348             if individualedge:
    2349                 curedge = []
    2350                 for i in xrange(len(edge)):
    2351                     curedge += edge[i]
    2352                
    2353             workscan._auto_sinusoid_baseline(mask, nwave, maxwavelength, clipthresh, clipniter, curedge, threshold, chan_avg_limit, getresidual, outlog, blfile)
    2354 
     2501            workscan._auto_sinusoid_baseline(mask, applyfft, fftmethod.lower(), str(fftthresh).lower(), workscan._parse_wn(addwn), workscan._parse_wn(rejwn), clipthresh, clipniter, normalise_edge_param(edge), threshold, chan_avg_limit, getresidual, outlog, blfile)
    23552502            workscan._add_history("auto_sinusoid_baseline", varlist)
    23562503           
    … …  
    23612508           
    23622509        except RuntimeError, e:
    2363             msg = "The fit failed, possibly because it didn't converge."
    2364             if rcParams["verbose"]:
    2365                 asaplog.push(str(e))
    2366                 asaplog.push(str(msg))
    2367                 return
    2368             else:
    2369                 raise RuntimeError(str(e)+'\n'+msg)
    2370 
     2510            raise_fitting_failure_exception(e)
    23712511
    23722512    @asaplog_post_dec
    … …  
    24052545        """
    24062546       
    2407         varlist = vars()
    2408        
    2409         if insitu is None: insitu = rcParams["insitu"]
    2410         if insitu:
    2411             workscan = self
    2412         else:
    2413             workscan = self.copy()
    2414 
    2415         nchan = workscan.nchan()
    2416        
    2417         if mask        is None: mask        = [True for i in xrange(nchan)]
    2418         if npiece      is None: npiece      = 2
    2419         if clipthresh  is None: clipthresh  = 3.0
    2420         if clipniter   is None: clipniter   = 0
    2421         if plot        is None: plot        = False
    2422         if getresidual is None: getresidual = True
    2423         if outlog      is None: outlog      = False
    2424         if blfile      is None: blfile      = ""
    2425 
    2426         outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
    2427        
    24282547        try:
     2548            varlist = vars()
     2549           
     2550            if insitu is None: insitu = rcParams['insitu']
     2551            if insitu:
     2552                workscan = self
     2553            else:
     2554                workscan = self.copy()
     2555
     2556            if mask        is None: mask        = [True for i in xrange(workscan.nchan())]
     2557            if npiece      is None: npiece      = 2
     2558            if clipthresh  is None: clipthresh  = 3.0
     2559            if clipniter   is None: clipniter   = 0
     2560            if plot        is None: plot        = False
     2561            if getresidual is None: getresidual = True
     2562            if outlog      is None: outlog      = False
     2563            if blfile      is None: blfile      = ''
     2564
    24292565            #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method.
    24302566            workscan._cspline_baseline(mask, npiece, clipthresh, clipniter, getresidual, outlog, blfile)
    2431            
    24322567            workscan._add_history("cspline_baseline", varlist)
    24332568           
    … …  
    24382573           
    24392574        except RuntimeError, e:
    2440             msg = "The fit failed, possibly because it didn't converge."
    2441             if rcParams["verbose"]:
    2442                 asaplog.push(str(e))
    2443                 asaplog.push(str(msg))
    2444                 return
    2445             else:
    2446                 raise RuntimeError(str(e)+'\n'+msg)
    2447 
    2448 
     2575            raise_fitting_failure_exception(e)
     2576
     2577    @asaplog_post_dec
    24492578    def auto_cspline_baseline(self, insitu=None, mask=None, npiece=None, clipthresh=None,
    24502579                              clipniter=None, edge=None, threshold=None,
    … …  
    25052634        """
    25062635
    2507         varlist = vars()
    2508 
    2509         if insitu is None: insitu = rcParams['insitu']
    2510         if insitu:
    2511             workscan = self
    2512         else:
    2513             workscan = self.copy()
    2514 
    2515         nchan = workscan.nchan()
    2516        
    2517         if mask           is None: mask           = [True for i in xrange(nchan)]
    2518         if npiece         is None: npiece         = 2
    2519         if clipthresh     is None: clipthresh     = 3.0
    2520         if clipniter      is None: clipniter      = 0
    2521         if edge           is None: edge           = (0, 0)
    2522         if threshold      is None: threshold      = 3
    2523         if chan_avg_limit is None: chan_avg_limit = 1
    2524         if plot           is None: plot           = False
    2525         if getresidual    is None: getresidual    = True
    2526         if outlog         is None: outlog         = False
    2527         if blfile         is None: blfile         = ""
    2528 
    2529         outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
    2530        
    2531         from asap.asaplinefind import linefinder
    2532         from asap import _is_sequence_or_number as _is_valid
    2533 
    2534         if not (isinstance(edge, list) or isinstance(edge, tuple)): edge = [ edge ]
    2535         individualedge = False;
    2536         if len(edge) > 1: individualedge = isinstance(edge[0], list) or isinstance(edge[0], tuple)
    2537 
    2538         if individualedge:
    2539             for edgepar in edge:
    2540                 if not _is_valid(edgepar, int):
    2541                     raise ValueError, "Each element of the 'edge' tuple has \
    2542                                        to be a pair of integers or an integer."
    2543         else:
    2544             if not _is_valid(edge, int):
    2545                 raise ValueError, "Parameter 'edge' has to be an integer or a \
    2546                                    pair of integers specified as a tuple. \
    2547                                    Nested tuples are allowed \
    2548                                    to make individual selection for different IFs."
    2549 
    2550             if len(edge) > 1:
    2551                 curedge = edge
     2636        try:
     2637            varlist = vars()
     2638
     2639            if insitu is None: insitu = rcParams['insitu']
     2640            if insitu:
     2641                workscan = self
    25522642            else:
    2553                 curedge = edge + edge
    2554 
    2555         try:
    2556             #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method.
    2557             if individualedge:
    2558                 curedge = []
    2559                 for i in xrange(len(edge)):
    2560                     curedge += edge[i]
    2561                
    2562             workscan._auto_cspline_baseline(mask, npiece, clipthresh, clipniter, curedge, threshold, chan_avg_limit, getresidual, outlog, blfile)
    2563 
     2643                workscan = self.copy()
     2644           
     2645            if mask           is None: mask           = [True for i in xrange(workscan.nchan())]
     2646            if npiece         is None: npiece         = 2
     2647            if clipthresh     is None: clipthresh     = 3.0
     2648            if clipniter      is None: clipniter      = 0
     2649            if edge           is None: edge           = (0, 0)
     2650            if threshold      is None: threshold      = 3
     2651            if chan_avg_limit is None: chan_avg_limit = 1
     2652            if plot           is None: plot           = False
     2653            if getresidual    is None: getresidual    = True
     2654            if outlog         is None: outlog         = False
     2655            if blfile         is None: blfile         = ''
     2656
     2657            #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method.
     2658            workscan._auto_cspline_baseline(mask, npiece, clipthresh, clipniter, normalise_edge_param(edge), threshold, chan_avg_limit, getresidual, outlog, blfile)
    25642659            workscan._add_history("auto_cspline_baseline", varlist)
    25652660           
    … …  
    25702665           
    25712666        except RuntimeError, e:
    2572             msg = "The fit failed, possibly because it didn't converge."
    2573             if rcParams["verbose"]:
    2574                 asaplog.push(str(e))
    2575                 asaplog.push(str(msg))
    2576                 return
    2577             else:
    2578                 raise RuntimeError(str(e)+'\n'+msg)
    2579 
     2667            raise_fitting_failure_exception(e)
    25802668
    25812669    @asaplog_post_dec
    … …  
    26062694        """
    26072695       
    2608         varlist = vars()
     2696        try:
     2697            varlist = vars()
    26092698       
    2610         if insitu is None: insitu = rcParams["insitu"]
    2611         if insitu:
    2612             workscan = self
    2613         else:
    2614             workscan = self.copy()
    2615 
    2616         nchan = workscan.nchan()
    2617        
    2618         if mask        is None: mask        = [True for i in xrange(nchan)]
    2619         if order       is None: order       = 0
    2620         if plot        is None: plot        = False
    2621         if getresidual is None: getresidual = True
    2622         if outlog      is None: outlog      = False
    2623         if blfile      is None: blfile      = ""
    2624 
    2625         outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
    2626        
    2627         try:
    2628             rows = xrange(workscan.nrow())
    2629            
    2630             #if len(rows) > 0: workscan._init_blinfo()
     2699            if insitu is None: insitu = rcParams["insitu"]
     2700            if insitu:
     2701                workscan = self
     2702            else:
     2703                workscan = self.copy()
     2704
     2705            if mask        is None: mask        = [True for i in xrange(workscan.nchan())]
     2706            if order       is None: order       = 0
     2707            if plot        is None: plot        = False
     2708            if getresidual is None: getresidual = True
     2709            if outlog      is None: outlog      = False
     2710            if blfile      is None: blfile      = ""
    26312711
    26322712            if plot:
     2713                outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
    26332714                if outblfile: blf = open(blfile, "a")
    26342715               
    26352716                f = fitter()
    26362717                f.set_function(lpoly=order)
     2718               
     2719                rows = xrange(workscan.nrow())
     2720                #if len(rows) > 0: workscan._init_blinfo()
     2721               
    26372722                for r in rows:
    26382723                    f.x = workscan._getabcissa(r)
    … …  
    26732758           
    26742759        except RuntimeError, e:
    2675             msg = "The fit failed, possibly because it didn't converge."
    2676             if rcParams["verbose"]:
    2677                 asaplog.push(str(e))
    2678                 asaplog.push(str(msg))
    2679                 return
    2680             else:
    2681                 raise RuntimeError(str(e)+'\n'+msg)
    2682 
    2683 
     2760            raise_fitting_failure_exception(e)
     2761
     2762    @asaplog_post_dec
    26842763    def auto_poly_baseline(self, insitu=None, mask=None, order=None, edge=None, threshold=None,
    26852764                           chan_avg_limit=None, plot=None, getresidual=None, outlog=None, blfile=None):
    … …  
    27312810        """
    27322811
    2733         varlist = vars()
    2734 
    2735         if insitu is None: insitu = rcParams['insitu']
    2736         if insitu:
    2737             workscan = self
    2738         else:
    2739             workscan = self.copy()
    2740 
    2741         nchan = workscan.nchan()
    2742        
    2743         if mask           is None: mask           = [True for i in xrange(nchan)]
    2744         if order          is None: order          = 0
    2745         if edge           is None: edge           = (0, 0)
    2746         if threshold      is None: threshold      = 3
    2747         if chan_avg_limit is None: chan_avg_limit = 1
    2748         if plot           is None: plot           = False
    2749         if getresidual    is None: getresidual    = True
    2750         if outlog         is None: outlog         = False
    2751         if blfile         is None: blfile         = ""
    2752 
    2753         outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
    2754        
    2755         from asap.asaplinefind import linefinder
    2756         from asap import _is_sequence_or_number as _is_valid
    2757 
    2758         if not (isinstance(edge, list) or isinstance(edge, tuple)): edge = [ edge ]
    2759         individualedge = False;
    2760         if len(edge) > 1: individualedge = isinstance(edge[0], list) or isinstance(edge[0], tuple)
    2761 
    2762         if individualedge:
    2763             for edgepar in edge:
    2764                 if not _is_valid(edgepar, int):
    2765                     raise ValueError, "Each element of the 'edge' tuple has \
    2766                                        to be a pair of integers or an integer."
    2767         else:
    2768             if not _is_valid(edge, int):
    2769                 raise ValueError, "Parameter 'edge' has to be an integer or a \
    2770                                    pair of integers specified as a tuple. \
    2771                                    Nested tuples are allowed \
    2772                                    to make individual selection for different IFs."
    2773 
    2774             if len(edge) > 1:
    2775                 curedge = edge
     2812        try:
     2813            varlist = vars()
     2814
     2815            if insitu is None: insitu = rcParams['insitu']
     2816            if insitu:
     2817                workscan = self
    27762818            else:
    2777                 curedge = edge + edge
    2778 
    2779         try:
    2780             rows = xrange(workscan.nrow())
    2781            
    2782             #if len(rows) > 0: workscan._init_blinfo()
     2819                workscan = self.copy()
     2820
     2821            if mask           is None: mask           = [True for i in xrange(workscan.nchan())]
     2822            if order          is None: order          = 0
     2823            if edge           is None: edge           = (0, 0)
     2824            if threshold      is None: threshold      = 3
     2825            if chan_avg_limit is None: chan_avg_limit = 1
     2826            if plot           is None: plot           = False
     2827            if getresidual    is None: getresidual    = True
     2828            if outlog         is None: outlog         = False
     2829            if blfile         is None: blfile         = ''
     2830
     2831            edge = normalise_edge_param(edge)
    27832832
    27842833            if plot:
     2834                outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
    27852835                if outblfile: blf = open(blfile, "a")
    27862836               
     2837                from asap.asaplinefind import linefinder
    27872838                fl = linefinder()
    27882839                fl.set_options(threshold=threshold,avg_limit=chan_avg_limit)
    27892840                fl.set_scan(workscan)
     2841               
    27902842                f = fitter()
    27912843                f.set_function(lpoly=order)
    27922844
     2845                rows = xrange(workscan.nrow())
     2846                #if len(rows) > 0: workscan._init_blinfo()
     2847               
    27932848                for r in rows:
    2794                     if individualedge:
    2795                         if len(edge) <= workscan.getif(r):
    2796                             raise RuntimeError, "Number of edge elements appear to " \
    2797                                   "be less than the number of IFs"
    2798                         else:
    2799                             curedge = edge[workscan.getif(r)]
    2800 
    2801                     fl.find_lines(r, mask_and(mask, workscan._getmask(r)), curedge)  # (CAS-1434)
     2849                    idx = 2*workscan.getif(r)
     2850                    fl.find_lines(r, mask_and(mask, workscan._getmask(r)), edge[idx:idx+2])  # (CAS-1434)
    28022851
    28032852                    f.x = workscan._getabcissa(r)
    … …  
    28272876
    28282877                if outblfile: blf.close()
    2829                
    28302878            else:
    2831                 if individualedge:
    2832                     curedge = []
    2833                     for i in xrange(len(edge)):
    2834                         curedge += edge[i]
    2835                
    2836                 workscan._auto_poly_baseline(mask, order, curedge, threshold, chan_avg_limit, getresidual, outlog, blfile)
     2879                workscan._auto_poly_baseline(mask, order, edge, threshold, chan_avg_limit, getresidual, outlog, blfile)
    28372880
    28382881            workscan._add_history("auto_poly_baseline", varlist)
    … …  
    28442887           
    28452888        except RuntimeError, e:
    2846             msg = "The fit failed, possibly because it didn't converge."
    2847             if rcParams["verbose"]:
    2848                 asaplog.push(str(e))
    2849                 asaplog.push(str(msg))
    2850                 return
    2851             else:
    2852                 raise RuntimeError(str(e)+'\n'+msg)
    2853 
    2854 
    2855     ### OBSOLETE ##################################################################
    2856     @asaplog_post_dec
    2857     def old_poly_baseline(self, mask=None, order=0, plot=False, uselin=False, insitu=None, rows=None):
    2858         """
    2859         Return a scan which has been baselined (all rows) by a polynomial.
    2860        
    2861         Parameters:
    2862 
    2863             mask:       an optional mask
    2864 
    2865             order:      the order of the polynomial (default is 0)
    2866 
    2867             plot:       plot the fit and the residual. In this each
    2868                         indivual fit has to be approved, by typing 'y'
    2869                         or 'n'
    2870 
    2871             uselin:     use linear polynomial fit
    2872 
    2873             insitu:     if False a new scantable is returned.
    2874                         Otherwise, the scaling is done in-situ
    2875                         The default is taken from .asaprc (False)
    2876 
    2877             rows:       row numbers of spectra to be processed.
    2878                         (default is None: for all rows)
    2879        
    2880         Example:
    2881             # return a scan baselined by a third order polynomial,
    2882             # not using a mask
    2883             bscan = scan.poly_baseline(order=3)
    2884 
    2885         """
    2886         if insitu is None: insitu = rcParams['insitu']
    2887         if not insitu:
    2888             workscan = self.copy()
    2889         else:
    2890             workscan = self
    2891         varlist = vars()
    2892         if mask is None:
    2893             mask = [True for i in xrange(self.nchan())]
    2894 
    2895         try:
    2896             f = fitter()
    2897             if uselin:
    2898                 f.set_function(lpoly=order)
    2899             else:
    2900                 f.set_function(poly=order)
    2901 
    2902             if rows == None:
    2903                 rows = xrange(workscan.nrow())
    2904             elif isinstance(rows, int):
    2905                 rows = [ rows ]
    2906            
    2907             if len(rows) > 0:
    2908                 self.blpars = []
    2909                 self.masklists = []
    2910                 self.actualmask = []
    2911            
    2912             for r in rows:
    2913                 f.x = workscan._getabcissa(r)
    2914                 f.y = workscan._getspectrum(r)
    2915                 f.mask = mask_and(mask, workscan._getmask(r))    # (CAS-1434)
    2916                 f.data = None
    2917                 f.fit()
    2918                 if plot:
    2919                     f.plot(residual=True)
    2920                     x = raw_input("Accept fit ( [y]/n ): ")
    2921                     if x.upper() == 'N':
    2922                         self.blpars.append(None)
    2923                         self.masklists.append(None)
    2924                         self.actualmask.append(None)
    2925                         continue
    2926                 workscan._setspectrum(f.fitter.getresidual(), r)
    2927                 self.blpars.append(f.get_parameters())
    2928                 self.masklists.append(workscan.get_masklist(f.mask, row=r, silent=True))
    2929                 self.actualmask.append(f.mask)
    2930 
    2931             if plot:
    2932                 f._p.unmap()
    2933                 f._p = None
    2934             workscan._add_history("poly_baseline", varlist)
    2935             if insitu:
    2936                 self._assign(workscan)
    2937             else:
    2938                 return workscan
    2939         except RuntimeError:
    2940             msg = "The fit failed, possibly because it didn't converge."
    2941             raise RuntimeError(msg)
     2889            raise_fitting_failure_exception(e)
    29422890
    29432891    def _init_blinfo(self):
Note: See TracChangeset for help on using the changeset viewer.