Changeset 2186 for trunk


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.


Location:
trunk
Files:
9 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):

  • trunk/src/MathUtils.cpp

    r2163 r2186  
    248248  }
    249249}
     250
     251void mathutil::doZeroOrderInterpolation(casa::Vector<casa::Float>& data,
     252                                        std::vector<bool>& mask) {
     253  int fstart = -1;
     254  int fend = -1;
     255  for (uInt i = 0; i < mask.size(); ++i) {
     256    if (!mask[i]) {
     257      fstart = i;
     258      while (!mask[i] && i < mask.size()) {
     259        fend = i;
     260        i++;
     261      }
     262    }
     263
     264    // execute interpolation as the following criteria:
     265    // (1) for a masked region inside the spectrum, replace the spectral
     266    //     values with the mean of those at the two channels just outside
     267    //     the both edges of the masked region.
     268    // (2) for a masked region at the spectral edge, replace the values
     269    //     with the one at the nearest non-masked channel.
     270    //     (ZOH, but bilateral)
     271    Float interp = 0.0;
     272    if (fstart-1 > 0) {
     273      interp = data[fstart-1];
     274      if (fend+1 < Int(data.nelements())) {
     275        interp = (interp + data[fend+1]) / 2.0;
     276      }
     277    } else {
     278      interp = data[fend+1];
     279    }
     280    if (fstart > -1 && fend > -1) {
     281      for (int j = fstart; j <= fend; ++j) {
     282        data[j] = interp;
     283      }
     284    }
     285
     286    fstart = -1;
     287    fend = -1;
     288  }
     289}

  • trunk/src/MathUtils.h

    r1819 r2186  
    6868 * @param hwidth half-width of the smoothing window
    6969 */
    70  void runningMedian(casa::Vector<casa::Float>& out,
     70void runningMedian(casa::Vector<casa::Float>& out,
    7171                   casa::Vector<casa::Bool>& outflag,
    7272                   const casa::Vector<casa::Float>& in,
    … …  
    7474                   float hwidth);
    7575
    76  void polyfit(casa::Vector<casa::Float>& out,
    77                    casa::Vector<casa::Bool>& outmask,
    78                    const casa::Vector<casa::Float>& in,
    79                    const casa::Vector<casa::Bool>& mask,
    80                    float hwidth, int order);
     76void polyfit(casa::Vector<casa::Float>& out,
     77             casa::Vector<casa::Bool>& outmask,
     78             const casa::Vector<casa::Float>& in,
     79             const casa::Vector<casa::Bool>& mask,
     80             float hwidth, int order);
    8181
    8282// Generate specified statistic
    … …  
    8585
    8686// Return a position of min or max value
    87  casa::IPosition minMaxPos(const casa::String& which,
     87casa::IPosition minMaxPos(const casa::String& which,
    8888                 const casa::MaskedArray<casa::Float>& data);
    8989
    … …  
    106106casa::Vector<casa::String> toVectorString(const std::vector<std::string>& in);
    107107
     108void doZeroOrderInterpolation(casa::Vector<casa::Float>& data,
     109                              std::vector<bool>& mask);
     110
    108111}
    109112

  • trunk/src/STMath.cpp

    r2177 r2186  
    4949#include <atnf/PKSIO/SrcType.h>
    5050
    51 #include "MathUtils.h"
    5251#include "RowAccumulator.h"
    5352#include "STAttr.h"
    … …  
    28272826  std::vector<float> res;
    28282827  Table tab = in->table();
    2829   ArrayColumn<Float> specCol(tab, "SPECTRA");
    2830   ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
    2831   FFTServer<Float,Complex> ffts;
     2828  std::vector<bool> mask;
    28322829
    28332830  if (whichrow.size() < 1) {          // for all rows (by default)
    28342831    int nrow = int(tab.nrow());
    28352832    for (int i = 0; i < nrow; ++i) {
    2836       Vector<Float> spec = specCol(i);
    2837       Vector<uChar> flag = flagCol(i);
    2838       doFFT(res, ffts, getRealImag, spec, flag);
     2833      res = in->execFFT(i, mask, getRealImag);
    28392834    }
    28402835  } else {                           // for specified rows
    28412836    for (uInt i = 0; i < whichrow.size(); ++i) {
    2842       Vector<Float> spec = specCol(whichrow[i]);
    2843       Vector<uChar> flag = flagCol(whichrow[i]);
    2844       doFFT(res, ffts, getRealImag, spec, flag);
     2837      res = in->execFFT(i, mask, getRealImag);
    28452838    }
    28462839  }
    … …  
    28492842}
    28502843
    2851 void asap::STMath::doFFT( std::vector<float>& out,
    2852                           FFTServer<Float,Complex>& ffts,
    2853                           bool getRealImag,
    2854                           Vector<Float>& spec,
    2855                           Vector<uChar>& flag )
    2856 {
    2857   doZeroOrderInterpolation(spec, flag);
    2858 
    2859   Vector<Complex> fftout;
    2860   ffts.fft0(fftout, spec);
    2861 
    2862   float norm = float(2.0/double(spec.size()));
    2863   if (getRealImag) {
    2864     for (uInt j = 0; j < fftout.size(); ++j) {
    2865       out.push_back(real(fftout[j])*norm);
    2866       out.push_back(imag(fftout[j])*norm);
    2867     }
    2868   } else {
    2869     for (uInt j = 0; j < fftout.size(); ++j) {
    2870       out.push_back(abs(fftout[j])*norm);
    2871       out.push_back(arg(fftout[j]));
    2872     }
    2873   }
    2874 
    2875 }
    2876 
    2877 void asap::STMath::doZeroOrderInterpolation( Vector<Float>& spec,
    2878                                              Vector<uChar>& flag )
    2879 {
    2880   int fstart = -1;
    2881   int fend = -1;
    2882   for (uInt i = 0; i < flag.nelements(); ++i) {
    2883     if (flag[i] > 0) {
    2884       fstart = i;
    2885       while (flag[i] > 0 && i < flag.nelements()) {
    2886         fend = i;
    2887         i++;
    2888       }
    2889     }
    2890 
    2891     // execute interpolation as the following criteria:
    2892     // (1) for a masked region inside the spectrum, replace the spectral
    2893     //     values with the mean of those at the two channels just outside
    2894     //     the both edges of the masked region.
    2895     // (2) for a masked region at the spectral edge, replace the values
    2896     //     with the one at the nearest non-masked channel.
    2897     //     (ZOH, but bilateral)
    2898     Float interp = 0.0;
    2899     if (fstart-1 > 0) {
    2900       interp = spec[fstart-1];
    2901       if (fend+1 < Int(spec.nelements())) {
    2902         interp = (interp + spec[fend+1]) / 2.0;
    2903       }
    2904     } else {
    2905       interp = spec[fend+1];
    2906     }
    2907     if (fstart > -1 && fend > -1) {
    2908       for (int j = fstart; j <= fend; ++j) {
    2909         spec[j] = interp;
    2910       }
    2911     }
    2912 
    2913     fstart = -1;
    2914     fend = -1;
    2915   }
    2916 }
    29172844
    29182845CountedPtr<Scantable>
    … …  
    29392866      Vector<Float> spec = specCol(i);
    29402867      Vector<uChar> flag = flagCol(i);
    2941       doZeroOrderInterpolation(spec, flag);
     2868      std::vector<bool> mask;
     2869      for (uInt j = 0; j < flag.nelements(); ++j) {
     2870        mask.push_back(!(flag[j]>0));
     2871      }
     2872      mathutil::doZeroOrderInterpolation(spec, mask);
    29422873
    29432874      Vector<Complex> lags;

  • trunk/src/STMath.h

    r2177 r2186  
    2121#include <casa/Utilities/CountedPtr.h>
    2222
    23 #include <scimath/Mathematics/FFTServer.h>
    2423#include <scimath/Mathematics/InterpolateArray1D.h>
    2524
    … …  
    433432                                        int index ) ;
    434433  double getMJD( string strtime ) ;
    435   void doFFT( std::vector<float>& out,
    436               casa::FFTServer< casa::Float, casa::Complex >& ffts,
    437               bool getRealImag,
    438               casa::Vector<casa::Float>& spec,
    439               casa::Vector<casa::uChar>& flag ) ;
    440   void doZeroOrderInterpolation( casa::Vector<casa::Float>& spec,
    441                                  casa::Vector<casa::uChar>& flag) ;
    442434
    443435  bool insitu_;

  • trunk/src/Scantable.cpp

    r2163 r2186  
    1111//
    1212#include <map>
    13 #include <fstream>
     13
     14#include <atnf/PKSIO/SrcType.h>
    1415
    1516#include <casa/aips.h>
     17#include <casa/iomanip.h>
    1618#include <casa/iostream.h>
    17 #include <casa/iomanip.h>
     19#include <casa/OS/File.h>
    1820#include <casa/OS/Path.h>
    19 #include <casa/OS/File.h>
    2021#include <casa/Arrays/Array.h>
     22#include <casa/Arrays/ArrayAccessor.h>
     23#include <casa/Arrays/ArrayLogical.h>
    2124#include <casa/Arrays/ArrayMath.h>
    2225#include <casa/Arrays/MaskArrMath.h>
    23 #include <casa/Arrays/ArrayLogical.h>
    24 #include <casa/Arrays/ArrayAccessor.h>
     26#include <casa/Arrays/Slice.h>
    2527#include <casa/Arrays/Vector.h>
    2628#include <casa/Arrays/VectorSTLIterator.h>
    27 #include <casa/Arrays/Slice.h>
    2829#include <casa/BasicMath/Math.h>
    2930#include <casa/BasicSL/Constants.h>
     31#include <casa/Containers/RecordField.h>
     32#include <casa/Logging/LogIO.h>
    3033#include <casa/Quanta/MVAngle.h>
    31 #include <casa/Containers/RecordField.h>
     34#include <casa/Quanta/MVTime.h>
    3235#include <casa/Utilities/GenSort.h>
    33 #include <casa/Logging/LogIO.h>
    34 
    35 #include <tables/Tables/TableParse.h>
    36 #include <tables/Tables/TableDesc.h>
    37 #include <tables/Tables/TableCopy.h>
    38 #include <tables/Tables/SetupNewTab.h>
    39 #include <tables/Tables/ScaColDesc.h>
    40 #include <tables/Tables/ArrColDesc.h>
    41 #include <tables/Tables/TableRow.h>
    42 #include <tables/Tables/TableVector.h>
    43 #include <tables/Tables/TableIter.h>
    44 
    45 #include <tables/Tables/ExprNode.h>
    46 #include <tables/Tables/TableRecord.h>
    47 #include <casa/Quanta/MVTime.h>
    48 #include <casa/Quanta/MVAngle.h>
    49 #include <measures/Measures/MeasRef.h>
    50 #include <measures/Measures/MeasTable.h>
     36
     37#include <coordinates/Coordinates/CoordinateUtil.h>
     38
    5139// needed to avoid error in .tcc
    5240#include <measures/Measures/MCDirection.h>
    5341//
    5442#include <measures/Measures/MDirection.h>
     43#include <measures/Measures/MEpoch.h>
    5544#include <measures/Measures/MFrequency.h>
    56 #include <measures/Measures/MEpoch.h>
     45#include <measures/Measures/MeasRef.h>
     46#include <measures/Measures/MeasTable.h>
     47#include <measures/TableMeasures/ScalarMeasColumn.h>
     48#include <measures/TableMeasures/TableMeasDesc.h>
    5749#include <measures/TableMeasures/TableMeasRefDesc.h>
    5850#include <measures/TableMeasures/TableMeasValueDesc.h>
    59 #include <measures/TableMeasures/TableMeasDesc.h>
    60 #include <measures/TableMeasures/ScalarMeasColumn.h>
    61 #include <coordinates/Coordinates/CoordinateUtil.h>
    62 
    63 #include <atnf/PKSIO/SrcType.h>
    64 #include "Scantable.h"
    65 #include "STPolLinear.h"
    66 #include "STPolCircular.h"
    67 #include "STPolStokes.h"
     51
     52#include <tables/Tables/ArrColDesc.h>
     53#include <tables/Tables/ExprNode.h>
     54#include <tables/Tables/ScaColDesc.h>
     55#include <tables/Tables/SetupNewTab.h>
     56#include <tables/Tables/TableCopy.h>
     57#include <tables/Tables/TableDesc.h>
     58#include <tables/Tables/TableIter.h>
     59#include <tables/Tables/TableParse.h>
     60#include <tables/Tables/TableRecord.h>
     61#include <tables/Tables/TableRow.h>
     62#include <tables/Tables/TableVector.h>
     63
     64#include "MathUtils.h"
    6865#include "STAttr.h"
    6966#include "STLineFinder.h"
    70 #include "MathUtils.h"
     67#include "STPolCircular.h"
     68#include "STPolLinear.h"
     69#include "STPolStokes.h"
     70#include "Scantable.h"
    7171
    7272using namespace casa;
    … …  
    18471847    setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
    18481848    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "polyBaseline()", fitter);
     1849    showProgressOnTerminal(whichrow, nRow);
    18491850  }
    18501851
    … …  
    19091910
    19101911    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoPolyBaseline()", fitter);
     1912    showProgressOnTerminal(whichrow, nRow);
    19111913  }
    19121914
    … …  
    19501952
    19511953    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()", pieceEdges, params);
     1954    showProgressOnTerminal(whichrow, nRow);
    19521955  }
    19531956
    … …  
    20182021
    20192022    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()", pieceEdges, params);
     2023    showProgressOnTerminal(whichrow, nRow);
    20202024  }
    20212025
    … …  
    22432247}
    22442248
    2245 void Scantable::sinusoidBaseline(const std::vector<bool>& mask, const std::vector<int>& nWaves, float maxWaveLength, float thresClip, int nIterClip, bool getResidual, bool outLogger, const std::string& blfile)
     2249  void Scantable::selectWaveNumbers(const int whichrow, const std::vector<bool>& chanMask, const bool applyFFT, const std::string& fftMethod, const std::string& fftThresh, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, std::vector<int>& nWaves)
     2250{
     2251  nWaves.clear();
     2252
     2253  if (applyFFT) {
     2254    string fftThAttr;
     2255    float fftThSigma;
     2256    int fftThTop;
     2257    parseThresholdExpression(fftThresh, fftThAttr, fftThSigma, fftThTop);
     2258    doSelectWaveNumbers(whichrow, chanMask, fftMethod, fftThSigma, fftThTop, fftThAttr, nWaves);
     2259  }
     2260
     2261  addAuxWaveNumbers(addNWaves, rejectNWaves, nWaves);
     2262}
     2263
     2264void Scantable::parseThresholdExpression(const std::string& fftThresh, std::string& fftThAttr, float& fftThSigma, int& fftThTop)
     2265{
     2266  uInt idxSigma = fftThresh.find("sigma");
     2267  uInt idxTop   = fftThresh.find("top");
     2268
     2269  if (idxSigma == fftThresh.size() - 5) {
     2270    std::istringstream is(fftThresh.substr(0, fftThresh.size() - 5));
     2271    is >> fftThSigma;
     2272    fftThAttr = "sigma";
     2273  } else if (idxTop == 0) {
     2274    std::istringstream is(fftThresh.substr(3));
     2275    is >> fftThTop;
     2276    fftThAttr = "top";
     2277  } else {
     2278    bool isNumber = true;
     2279    for (uInt i = 0; i < fftThresh.size()-1; ++i) {
     2280      char ch = (fftThresh.substr(i, 1).c_str())[0];
     2281      if (!(isdigit(ch) || (fftThresh.substr(i, 1) == "."))) {
     2282        isNumber = false;
     2283        break;
     2284      }
     2285    }
     2286    if (isNumber) {
     2287      std::istringstream is(fftThresh);
     2288      is >> fftThSigma;
     2289      fftThAttr = "sigma";
     2290    } else {
     2291      throw(AipsError("fftthresh has a wrong value"));
     2292    }
     2293  }
     2294}
     2295
     2296void Scantable::doSelectWaveNumbers(const int whichrow, const std::vector<bool>& chanMask, const std::string& fftMethod, const float fftThSigma, const int fftThTop, const std::string& fftThAttr, std::vector<int>& nWaves)
     2297{
     2298  std::vector<float> fspec;
     2299  if (fftMethod == "fft") {
     2300    fspec = execFFT(whichrow, chanMask, false, true);
     2301  //} else if (fftMethod == "lsp") {
     2302  //  fspec = lombScarglePeriodogram(whichrow);
     2303  }
     2304
     2305  if (fftThAttr == "sigma") {
     2306    float mean  = 0.0;
     2307    float mean2 = 0.0;
     2308    for (uInt i = 0; i < fspec.size(); ++i) {
     2309      mean  += fspec[i];
     2310      mean2 += fspec[i]*fspec[i];
     2311    }
     2312    mean  /= float(fspec.size());
     2313    mean2 /= float(fspec.size());
     2314    float thres = mean + fftThSigma * float(sqrt(mean2 - mean*mean));
     2315
     2316    for (uInt i = 0; i < fspec.size(); ++i) {
     2317      if (fspec[i] >= thres) {
     2318        nWaves.push_back(i);
     2319      }
     2320    }
     2321
     2322  } else if (fftThAttr == "top") {
     2323    for (int i = 0; i < fftThTop; ++i) {
     2324      float max = 0.0;
     2325      int maxIdx = 0;
     2326      for (uInt j = 0; j < fspec.size(); ++j) {
     2327        if (fspec[j] > max) {
     2328          max = fspec[j];
     2329          maxIdx = j;
     2330        }
     2331      }
     2332      nWaves.push_back(maxIdx);
     2333      fspec[maxIdx] = 0.0;
     2334    }
     2335
     2336  }
     2337
     2338  if (nWaves.size() > 1) {
     2339    sort(nWaves.begin(), nWaves.end());
     2340  }
     2341}
     2342
     2343void Scantable::addAuxWaveNumbers(const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, std::vector<int>& nWaves)
     2344{
     2345  for (uInt i = 0; i < addNWaves.size(); ++i) {
     2346    bool found = false;
     2347    for (uInt j = 0; j < nWaves.size(); ++j) {
     2348      if (nWaves[j] == addNWaves[i]) {
     2349        found = true;
     2350        break;
     2351      }
     2352    }
     2353    if (!found) nWaves.push_back(addNWaves[i]);
     2354  }
     2355
     2356  for (uInt i = 0; i < rejectNWaves.size(); ++i) {
     2357    for (std::vector<int>::iterator j = nWaves.begin(); j != nWaves.end(); ) {
     2358      if (*j == rejectNWaves[i]) {
     2359        j = nWaves.erase(j);
     2360      } else {
     2361        ++j;
     2362      }
     2363    }
     2364  }
     2365
     2366  if (nWaves.size() > 1) {
     2367    sort(nWaves.begin(), nWaves.end());
     2368    unique(nWaves.begin(), nWaves.end());
     2369  }
     2370}
     2371
     2372void Scantable::sinusoidBaseline(const std::vector<bool>& mask, const bool applyFFT, const std::string& fftMethod, const std::string& fftThresh, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, float thresClip, int nIterClip, bool getResidual, bool outLogger, const std::string& blfile)
    22462373{
    22472374  ofstream ofs;
    … …  
    22672394  int nRow = nrow();
    22682395  std::vector<bool> chanMask;
     2396  std::vector<int> nWaves;
    22692397
    22702398  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    22712399    chanMask = getCompositeChanMask(whichrow, mask);
     2400    selectWaveNumbers(whichrow, chanMask, applyFFT, fftMethod, fftThresh, addNWaves, rejectNWaves, nWaves);
     2401
     2402    //FOR DEBUGGING------------
     2403    if (whichrow < 0) {// == nRow -1) {
     2404      cout << "+++ i=" << setw(3) << whichrow << ", IF=" << setw(2) << getIF(whichrow);
     2405      if (applyFFT) {
     2406          cout << "[ ";
     2407          for (uInt j = 0; j < nWaves.size(); ++j) {
     2408            cout << nWaves[j] << ", ";
     2409          }
     2410          cout << " ]    " << endl;
     2411      }
     2412      cout << flush;
     2413    }
     2414    //-------------------------
     2415
    22722416    //fitBaseline(chanMask, whichrow, fitter);
    22732417    //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
    22742418    std::vector<float> params;
    2275     std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, maxWaveLength, params, thresClip, nIterClip, getResidual);
     2419    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, params, thresClip, nIterClip, getResidual);
    22762420    setSpectrum(res, whichrow);
    22772421    //
    22782422
    22792423    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "sinusoidBaseline()", params);
     2424    showProgressOnTerminal(whichrow, nRow);
    22802425  }
    22812426
    … …  
    22832428}
    22842429
    2285 void Scantable::autoSinusoidBaseline(const std::vector<bool>& mask, const std::vector<int>& nWaves, float maxWaveLength, float thresClip, int nIterClip, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool getResidual, bool outLogger, const std::string& blfile)
     2430void Scantable::autoSinusoidBaseline(const std::vector<bool>& mask, const bool applyFFT, const std::string& fftMethod, const std::string& fftThresh, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, float thresClip, int nIterClip, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool getResidual, bool outLogger, const std::string& blfile)
    22862431{
    22872432  ofstream ofs;
    … …  
    23072452  int nRow = nrow();
    23082453  std::vector<bool> chanMask;
     2454  std::vector<int> nWaves;
     2455
    23092456  int minEdgeSize = getIFNos().size()*2;
    23102457  STLineFinder lineFinder = STLineFinder();
    … …  
    23362483    //-------------------------------------------------------
    23372484
     2485    selectWaveNumbers(whichrow, chanMask, applyFFT, fftMethod, fftThresh, addNWaves, rejectNWaves, nWaves);
    23382486
    23392487    //fitBaseline(chanMask, whichrow, fitter);
    23402488    //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
    23412489    std::vector<float> params;
    2342     std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, maxWaveLength, params, thresClip, nIterClip, getResidual);
     2490    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, params, thresClip, nIterClip, getResidual);
    23432491    setSpectrum(res, whichrow);
    23442492    //
    23452493
    23462494    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoSinusoidBaseline()", params);
     2495    showProgressOnTerminal(whichrow, nRow);
    23472496  }
    23482497
    … …  
    23502499}
    23512500
    2352 std::vector<float> Scantable::doSinusoidFitting(const std::vector<float>& data, const std::vector<bool>& mask, const std::vector<int>& waveNumbers, float maxWaveLength, std::vector<float>& params, float thresClip, int nIterClip, bool getResidual)
     2501std::vector<float> Scantable::doSinusoidFitting(const std::vector<float>& data, const std::vector<bool>& mask, const std::vector<int>& waveNumbers, std::vector<float>& params, float thresClip, int nIterClip, bool getResidual)
    23532502{
    23542503  if (data.size() != mask.size()) {
    … …  
    23592508  }
    23602509  if (waveNumbers.size() == 0) {
    2361     throw(AipsError("missing wave number info"));
     2510    throw(AipsError("no wave numbers given"));
    23622511  }
    23632512  std::vector<int> nWaves;  // sorted and uniqued array of wave numbers
    … …  
    23902539
    23912540  const double PI = 6.0 * asin(0.5); // PI (= 3.141592653...)
    2392   double baseXFactor = 2.0*PI/(double)maxWaveLength/(double)(nChan-1);  //the denominator (nChan-1) should be changed to (xdata[nChan-1]-xdata[0]) for accepting x-values given in velocity or frequency when this function is moved to fitter. (2011/03/30 WK)
     2541  double baseXFactor = 2.0*PI/(double)(nChan-1);  //the denominator (nChan-1) should be changed to (xdata[nChan-1]-xdata[0]) for accepting x-values given in velocity or frequency when this function is moved to fitter. (2011/03/30 WK)
    23932542
    23942543  // xArray : contains elemental values for computing the least-square matrix.
    … …  
    25712720std::vector<bool> Scantable::getCompositeChanMask(int whichrow, const std::vector<bool>& inMask)
    25722721{
    2573   std::vector<bool> chanMask = getMask(whichrow);
    2574   uInt chanMaskSize = chanMask.size();
    2575   if (chanMaskSize != inMask.size()) {
    2576     throw(AipsError("different mask sizes"));
    2577   }
    2578   for (uInt i = 0; i < chanMaskSize; ++i) {
    2579     chanMask[i] = chanMask[i] && inMask[i];
    2580   }
    2581 
    2582   return chanMask;
     2722  std::vector<bool> mask = getMask(whichrow);
     2723  uInt maskSize = mask.size();
     2724  if (maskSize != inMask.size()) {
     2725    throw(AipsError("mask sizes are not the same."));
     2726  }
     2727  for (uInt i = 0; i < maskSize; ++i) {
     2728    mask[i] = mask[i] && inMask[i];
     2729  }
     2730
     2731  return mask;
    25832732}
    25842733
    … …  
    26102759
    26112760/* for poly. the variations of outputFittingResult() should be merged into one eventually (2011/3/10 WK)  */
    2612 void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, Fitter& fitter) {
     2761void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, Fitter& fitter)
     2762{
    26132763  if (outLogger || outTextFile) {
    26142764    std::vector<float> params = fitter.getParameters();
    … …  
    26282778
    26292779/* for cspline. will be merged once cspline is available in fitter (2011/3/10 WK) */
    2630 void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, const std::vector<int>& edge, const std::vector<float>& params) {
     2780void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, const std::vector<int>& edge, const std::vector<float>& params)
     2781{
    26312782  if (outLogger || outTextFile) {
    26322783    float rms = getRms(chanMask, whichrow);
    … …  
    26462797
    26472798/* for sinusoid. will be merged once sinusoid is available in fitter (2011/3/10 WK) */
    2648 void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, const std::vector<float>& params) {
     2799void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, const std::vector<float>& params)
     2800{
    26492801  if (outLogger || outTextFile) {
    26502802    float rms = getRms(chanMask, whichrow);
    … …  
    26632815}
    26642816
    2665 float Scantable::getRms(const std::vector<bool>& mask, int whichrow) {
     2817void Scantable::showProgressOnTerminal(const int nProcessed, const int nTotal, const int nTotalThreshold)
     2818{
     2819  if (nTotal >= nTotalThreshold) {
     2820    int nInterval = int(floor(double(nTotal)/100.0));
     2821    if (nInterval == 0) nInterval++;
     2822
     2823    if (nProcessed == 0) {
     2824      printf("\x1b[31m\x1b[1m");             //set red color, highlighted
     2825      printf("[  0%%]");
     2826      printf("\x1b[39m\x1b[0m");             //default attributes
     2827      fflush(NULL);
     2828    } else if (nProcessed % nInterval == 0) {
     2829      printf("\r\x1b[1C");                   //go to the 2nd column
     2830      printf("\x1b[31m\x1b[1m");             //set red color, highlighted
     2831      printf("%3d", (int)(100.0*(double(nProcessed+1))/(double(nTotal))) );
     2832      printf("\x1b[39m\x1b[0m");             //default attributes
     2833      printf("\x1b[2C");                     //go to the end of line
     2834      fflush(NULL);
     2835    }
     2836    if (nProcessed == nTotal - 1) {
     2837      printf("\r\x1b[K");                    //clear
     2838      fflush(NULL);
     2839    }
     2840  }
     2841}
     2842
     2843std::vector<float> Scantable::execFFT(const int whichrow, const std::vector<bool>& inMask, bool getRealImag, bool getAmplitudeOnly)
     2844{
     2845  std::vector<bool>  mask = getMask(whichrow);
     2846
     2847  if (inMask.size() > 0) {
     2848    uInt maskSize = mask.size();
     2849    if (maskSize != inMask.size()) {
     2850      throw(AipsError("mask sizes are not the same."));
     2851    }
     2852    for (uInt i = 0; i < maskSize; ++i) {
     2853      mask[i] = mask[i] && inMask[i];
     2854    }
     2855  }
     2856
     2857  Vector<Float> spec = getSpectrum(whichrow);
     2858  mathutil::doZeroOrderInterpolation(spec, mask);
     2859
     2860  FFTServer<Float,Complex> ffts;
     2861  Vector<Complex> fftres;
     2862  ffts.fft0(fftres, spec);
     2863
     2864  std::vector<float> res;
     2865  float norm = float(2.0/double(spec.size()));
     2866
     2867  if (getRealImag) {
     2868    for (uInt i = 0; i < fftres.size(); ++i) {
     2869      res.push_back(real(fftres[i])*norm);
     2870      res.push_back(imag(fftres[i])*norm);
     2871    }
     2872  } else {
     2873    for (uInt i = 0; i < fftres.size(); ++i) {
     2874      res.push_back(abs(fftres[i])*norm);
     2875      if (!getAmplitudeOnly) res.push_back(arg(fftres[i]));
     2876    }
     2877  }
     2878
     2879  return res;
     2880}
     2881
     2882
     2883float Scantable::getRms(const std::vector<bool>& mask, int whichrow)
     2884{
    26662885  Vector<Float> spec;
    26672886  specCol_.get(whichrow, spec);
    … …  
    26852904
    26862905
    2687 std::string Scantable::formatBaselineParamsHeader(int whichrow,
    2688                                                   const std::string& masklist,
    2689                                                   bool verbose) const
     2906std::string Scantable::formatBaselineParamsHeader(int whichrow, const std::string& masklist, bool verbose) const
    26902907{
    26912908  ostringstream oss;
    … …  
    27222939}
    27232940
    2724   std::string Scantable::formatBaselineParams(const std::vector<float>& params,
    2725                                               const std::vector<bool>& fixed,
    2726                                               float rms,
    2727                                               const std::string& masklist,
    2728                                               int whichrow,
    2729                                               bool verbose,
    2730                                               int start, int count,
    2731                                               bool resetparamid) const
     2941std::string Scantable::formatBaselineParams(const std::vector<float>& params,
     2942                                            const std::vector<bool>& fixed,
     2943                                            float rms,
     2944                                            const std::string& masklist,
     2945                                            int whichrow,
     2946                                            bool verbose,
     2947                                            int start, int count,
     2948                                            bool resetparamid) const
    27322949{
    27332950  int nParam = (int)(params.size());
    … …  
    29043121}
    29053122
    2906 /*
    2907 STFitEntry Scantable::polyBaseline(const std::vector<bool>& mask, int order, int rowno)
    2908 {
    2909   Fitter fitter = Fitter();
    2910   fitter.setExpression("poly", order);
    2911 
    2912   std::vector<bool> fmask = getMask(rowno);
    2913   if (fmask.size() != mask.size()) {
    2914     throw(AipsError("different mask sizes"));
    2915   }
    2916   for (int i = 0; i < fmask.size(); ++i) {
    2917     fmask[i] = fmask[i] && mask[i];
    2918   }
    2919 
    2920   fitBaseline(fmask, rowno, fitter);
    2921   setSpectrum(fitter.getResidual(), rowno);
    2922   return fitter.getFitEntry();
    2923 }
    2924 */
    29253123
    29263124}

  • trunk/src/Scantable.h

    r2163 r2186  
    1414
    1515// STL
     16#include <fstream>
     17#include <iostream>
     18#include <sstream>
    1619#include <string>
    1720#include <vector>
    18 #include <iostream>
    19 #include <fstream>
    2021// AIPS++
    2122#include <casa/aips.h>
     23#include <casa/Arrays/MaskedArray.h>
     24#include <casa/Arrays/Vector.h>
     25#include <casa/BasicSL/String.h>
    2226#include <casa/Containers/Record.h>
    23 #include <casa/Arrays/MaskedArray.h>
    24 #include <casa/BasicSL/String.h>
     27#include <casa/Exceptions/Error.h>
     28#include <casa/Quanta/Quantum.h>
    2529#include <casa/Utilities/CountedPtr.h>
    2630
    27 #include <tables/Tables/Table.h>
     31#include <coordinates/Coordinates/SpectralCoordinate.h>
     32
     33#include <measures/TableMeasures/ScalarMeasColumn.h>
     34
     35#include <scimath/Mathematics/FFTServer.h>
     36
    2837#include <tables/Tables/ArrayColumn.h>
    2938#include <tables/Tables/ScalarColumn.h>
    30 
    31 #include <measures/TableMeasures/ScalarMeasColumn.h>
    32 
    33 #include <coordinates/Coordinates/SpectralCoordinate.h>
    34 
    35 #include <casa/Arrays/Vector.h>
    36 #include <casa/Quanta/Quantum.h>
    37 
    38 #include <casa/Exceptions/Error.h>
     39#include <tables/Tables/Table.h>
    3940
    4041#include "Logger.h"
    41 #include "STHeader.h"
    42 #include "STFrequencies.h"
    43 #include "STWeather.h"
    44 #include "STFocus.h"
    45 #include "STTcal.h"
    46 #include "STMolecules.h"
    47 #include "STSelector.h"
    48 #include "STHistory.h"
    49 #include "STPol.h"
     42#include "MathUtils.h"
    5043#include "STFit.h"
    5144#include "STFitEntry.h"
    5245#include "STFitter.h"
     46#include "STFocus.h"
     47#include "STFrequencies.h"
     48#include "STHeader.h"
     49#include "STHistory.h"
     50#include "STMolecules.h"
     51#include "STPol.h"
     52#include "STSelector.h"
     53#include "STTcal.h"
     54#include "STWeather.h"
    5355
    5456namespace asap {
    … …  
    528530                               const std::string& blfile="");
    529531  void sinusoidBaseline(const std::vector<bool>& mask,
    530                         const std::vector<int>& nWaves,
    531                         float maxWaveLength,
     532                        const bool applyFFT,
     533                        const std::string& fftMethod,
     534                        const std::string& fftThresh,
     535                        const std::vector<int>& addNWaves,
     536                        const std::vector<int>& rejectNWaves,
    532537                        float thresClip,
    533538                        int nIterClip,
    … …  
    536541                        const std::string& blfile="");
    537542  void autoSinusoidBaseline(const std::vector<bool>& mask,
    538                             const std::vector<int>& nWaves,
    539                             float maxWaveLength,
     543                            const bool applyFFT,
     544                            const std::string& fftMethod,
     545                            const std::string& fftThresh,
     546                            const std::vector<int>& addNWaves,
     547                            const std::vector<int>& rejectNWaves,
    540548                            float thresClip,
    541549                            int nIterClip,
    … …  
    546554                            bool outLogger=false,
    547555                            const std::string& blfile="");
     556  std::vector<float> execFFT(const int whichrow,
     557                             const std::vector<bool>& inMask,
     558                             bool getRealImag=false,
     559                             bool getAmplitudeOnly=false);
    548560  float getRms(const std::vector<bool>& mask, int whichrow);
    549561  std::string formatBaselineParams(const std::vector<float>& params,
    … …  
    689701                                       const std::vector<bool>& mask,
    690702                                       const std::vector<int>& waveNumbers,
    691                                        float maxWaveLength,
    692703                                       std::vector<float>& params,
    693704                                       float thresClip=3.0,
    694705                                       int nIterClip=1,
    695706                                       bool getResidual=true);
     707  void selectWaveNumbers(const int whichrow,
     708                         const std::vector<bool>& chanMask,
     709                         const bool applyFFT,
     710                         const std::string& fftMethod,
     711                         const std::string& fftThresh,
     712                         const std::vector<int>& addNWaves,
     713                         const std::vector<int>& rejectNWaves,
     714                         std::vector<int>& nWaves);
     715  void parseThresholdExpression(const std::string& fftThresh,
     716                                std::string& fftThAttr,
     717                                float& fftThSigma,
     718                                int& fftThTop);
     719  void doSelectWaveNumbers(const int whichrow,
     720                           const std::vector<bool>& chanMask,
     721                           const std::string& fftMethod,
     722                           const float fftThSigma,
     723                           const int fftThTop,
     724                           const std::string& fftThAttr,
     725                           std::vector<int>& nWaves);
     726  void addAuxWaveNumbers(const std::vector<int>& addNWaves,
     727                         const std::vector<int>& rejectNWaves,
     728                         std::vector<int>& nWaves);
    696729  bool hasSameNchanOverIFs();
    697730  std::string getMaskRangeList(const std::vector<bool>& mask,
    … …  
    708741  void outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, std::ofstream& ofs, const casa::String& funcName, const std::vector<int>& edge, const std::vector<float>& params);
    709742  void outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, std::ofstream& ofs, const casa::String& funcName, const std::vector<float>& params);
     743  void showProgressOnTerminal(const int nProcessed, const int nTotal, const int nTotalThreshold=1000);
    710744
    711745  void applyChanFlag( casa::uInt whichrow, const std::vector<bool>& msk, casa::uChar flagval);

  • trunk/src/ScantableWrapper.h

    r2178 r2186  
    276276  { table_->autoCubicSplineBaseline(mask, npiece, clipthresh, clipniter, edge, threshold, chan_avg_limit, getresidual, outlog, blfile); }
    277277
    278   void sinusoidBaseline(const std::vector<bool>& mask, const std::vector<int>& nwave, float maxwavelength, float clipthresh, int clipniter, bool getresidual=true, bool outlog=false, const std::string& blfile="")
    279   { table_->sinusoidBaseline(mask, nwave, maxwavelength, clipthresh, clipniter, getresidual, outlog, blfile); }
    280 
    281   void autoSinusoidBaseline(const std::vector<bool>& mask, const std::vector<int>& nwave, float maxwavelength, float clipthresh, int clipniter, const std::vector<int>& edge, float threshold=5.0, int chan_avg_limit=1, bool getresidual=true, bool outlog=false, const std::string& blfile="")
    282   { table_->autoSinusoidBaseline(mask, nwave, maxwavelength, clipthresh, clipniter, edge, threshold, chan_avg_limit, getresidual, outlog, blfile); }
     278  void sinusoidBaseline(const std::vector<bool>& mask, const bool applyfft, const std::string& fftmethod, const std::string& fftthresh, const std::vector<int>& addwn, const std::vector<int>& rejwn, float clipthresh, int clipniter, bool getresidual=true, bool outlog=false, const std::string& blfile="")
     279  { table_->sinusoidBaseline(mask, applyfft, fftmethod, fftthresh, addwn, rejwn, clipthresh, clipniter, getresidual, outlog, blfile); }
     280
     281  void autoSinusoidBaseline(const std::vector<bool>& mask, const bool applyfft, const std::string& fftmethod, const std::string& fftthresh, const std::vector<int>& addwn, const std::vector<int>& rejwn, float clipthresh, int clipniter, const std::vector<int>& edge, float threshold=5.0, int chan_avg_limit=1, bool getresidual=true, bool outlog=false, const std::string& blfile="")
     282  { table_->autoSinusoidBaseline(mask, applyfft, fftmethod, fftthresh, addwn, rejwn, clipthresh, clipniter, edge, threshold, chan_avg_limit, getresidual, outlog, blfile); }
    283283
    284284  float getRms(const std::vector<bool>& mask, int whichrow)
    … …  
    293293  bool getFlagtraFast(int whichrow=0) const
    294294  { return table_->getFlagtraFast(casa::uInt(whichrow)); }
     295
     296  std::vector<float> execFFT(int whichrow, const std::vector<bool>& mask, bool getRealImag=false, bool getAmplitudeOnly=false)
     297  { return table_->execFFT(whichrow, mask, getRealImag, getAmplitudeOnly); }
    295298
    296299

  • trunk/src/python_Scantable.cpp

    r2178 r2186  
    154154    .def("_getflagtrafast", &ScantableWrapper::getFlagtraFast,
    155155         (boost::python::arg("whichrow")=0) )
     156    .def("_fft", &ScantableWrapper::execFFT)
    156157    //.def("_sspline_baseline", &ScantableWrapper::smoothingSplineBaseline)
    157158    //.def("_test_cin", &ScantableWrapper::testCin)
Note: See TracChangeset for help on using the changeset viewer.