Changeset 2767


Ignore:
Timestamp:
02/08/13 22:23:22 (12 years ago)
Author:
WataruKawasaki
Message:

New Development: Yes

JIRA Issue: Yes CAS-4794

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): sd

Description: functions to apply/write STBaselineTable in which baseline parameters stored.


Location:
trunk
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • trunk/python/scantable.py

    r2761 r2767  
    173173def pack_progress_params(showprogress, minnrow):
    174174    return str(showprogress).lower() + ',' + str(minnrow)
     175
     176def pack_blinfo(blinfo, maxirow):
     177    """\
     178    convert a dictionary or a list of dictionaries of baseline info
     179    into a list of comma-separated strings.
     180    """
     181    if isinstance(blinfo, dict):
     182        res = do_pack_blinfo(blinfo, maxirow)
     183        return [res] if res != '' else []
     184    elif isinstance(blinfo, list) or isinstance(blinfo, tuple):
     185        res = []
     186        for i in xrange(len(blinfo)):
     187            resi = do_pack_blinfo(blinfo[i], maxirow)
     188            if resi != '':
     189                res.append(resi)
     190    return res
     191
     192def do_pack_blinfo(blinfo, maxirow):
     193    """\
     194    convert a dictionary of baseline info for a spectrum into
     195    a comma-separated string.
     196    """
     197    dinfo = {}
     198    for key in ['row', 'blfunc', 'masklist']:
     199        if blinfo.has_key(key):
     200            val = blinfo[key]
     201            if key == 'row':
     202                irow = val
     203            if isinstance(val, list) or isinstance(val, tuple):
     204                slval = []
     205                for i in xrange(len(val)):
     206                    if isinstance(val[i], list) or isinstance(val[i], tuple):
     207                        for j in xrange(len(val[i])):
     208                            slval.append(str(val[i][j]))
     209                    else:
     210                        slval.append(str(val[i]))
     211                sval = ",".join(slval)
     212            else:
     213                sval = str(val)
     214
     215            dinfo[key] = sval
     216        else:
     217            raise ValueError("'"+key+"' is missing in blinfo.")
     218
     219    if irow >= maxirow: return ''
     220   
     221    for key in ['order', 'npiece', 'nwave']:
     222        if blinfo.has_key(key):
     223            val = blinfo[key]
     224            if isinstance(val, list) or isinstance(val, tuple):
     225                slval = []
     226                for i in xrange(len(val)):
     227                    slval.append(str(val[i]))
     228                sval = ",".join(slval)
     229            else:
     230                sval = str(val)
     231            dinfo[key] = sval
     232
     233    blfunc = dinfo['blfunc']
     234    fspec_keys = {'poly': 'order', 'chebyshev': 'order', 'cspline': 'npiece', 'sinusoid': 'nwave'}
     235
     236    fspec_key = fspec_keys[blfunc]
     237    if not blinfo.has_key(fspec_key):
     238        raise ValueError("'"+fspec_key+"' is missing in blinfo.")
     239
     240    clip_params_n = 0
     241    for key in ['clipthresh', 'clipniter']:
     242        if blinfo.has_key(key):
     243            clip_params_n += 1
     244            dinfo[key] = str(blinfo[key])
     245
     246    if clip_params_n == 0:
     247        dinfo['clipthresh'] = '0.0'
     248        dinfo['clipniter']  = '0'
     249    elif clip_params_n != 2:
     250        raise ValueError("both 'clipthresh' and 'clipniter' must be given for n-sigma clipping.")
     251
     252    lf_params_n = 0
     253    for key in ['thresh', 'edge', 'chan_avg_limit']:
     254        if blinfo.has_key(key):
     255            lf_params_n += 1
     256            val = blinfo[key]
     257            if isinstance(val, list) or isinstance(val, tuple):
     258                slval = []
     259                for i in xrange(len(val)):
     260                    slval.append(str(val[i]))
     261                sval = ",".join(slval)
     262            else:
     263                sval = str(val)
     264            dinfo[key] = sval
     265
     266    if lf_params_n == 3:
     267        dinfo['use_linefinder'] = 'true'
     268    elif lf_params_n == 1:
     269        dinfo['use_linefinder'] = 'false'
     270        dinfo['thresh']         = ''
     271        dinfo['edge']           = ''
     272        dinfo['chan_avg_limit'] = ''
     273    else:
     274        raise ValueError("all of 'thresh', 'edge' and 'chan_avg_limit' must be given to use linefinder.")
     275   
     276    slblinfo = [dinfo['row'], blfunc, dinfo[fspec_key], dinfo['masklist'], \
     277                dinfo['clipthresh'], dinfo['clipniter'], \
     278                dinfo['use_linefinder'], dinfo['thresh'], dinfo['edge'], dinfo['chan_avg_limit']]
     279   
     280    return ":".join(slblinfo)
     281
     282def parse_fitresult(sres):
     283    """\
     284    Parse the returned value of apply_bltable() or sub_baseline() and
     285    extract row number, the best-fit coefficients and rms, then pack
     286    them into a dictionary.
     287    The input value is generated by Scantable::packFittingResults() and
     288    formatted as 'row:coeff[0],coeff[1],..,coeff[n-1]:rms'.
     289    """
     290    res = []
     291    for i in xrange(len(sres)):
     292        (srow, scoeff, srms) = sres[i].split(":")
     293        row = int(srow)
     294        rms = float(srms)
     295        lscoeff = scoeff.split(",")
     296        coeff = []
     297        for j in xrange(len(lscoeff)):
     298            coeff.append(float(lscoeff[j]))
     299        res.append({'row': row, 'coeff': coeff, 'rms': rms})
     300
     301    return res
    175302
    176303class scantable(Scantable):
     
    25562683
    25572684    @asaplog_post_dec
     2685    def apply_bltable(self, insitu=None, inbltable=None, outbltable=None, overwrite=None):
     2686        """\
     2687        Subtract baseline based on parameters written in Baseline Table.
     2688
     2689        Parameters:
     2690            insitu:        if False a new scantable is returned.
     2691                           Otherwise, the scaling is done in-situ
     2692                           The default is taken from .asaprc (False)
     2693            inbltable:     name of input baseline table. The row number of
     2694                           scantable and that of inbltable must be
     2695                           identical.
     2696            outbltable:    name of output baseline table where baseline
     2697                           parameters and fitting results recorded.
     2698                           default is ''(no output).
     2699            overwrite:     if True, overwrites the existing baseline table
     2700                           specified in outbltable.
     2701                           default is False.
     2702        """
     2703
     2704        try:
     2705            varlist = vars()
     2706            if inbltable      is None: raise ValueError("bltable missing.")
     2707            if outbltable     is None: outbltable     = ''
     2708            if overwrite      is None: overwrite      = False
     2709
     2710            if insitu is None: insitu = rcParams['insitu']
     2711            if insitu:
     2712                workscan = self
     2713            else:
     2714                workscan = self.copy()
     2715
     2716            sres = workscan._apply_bltable(inbltable,
     2717                                           outbltable,
     2718                                           os.path.exists(outbltable),
     2719                                           overwrite)
     2720            res = parse_fitresult(sres)
     2721
     2722            workscan._add_history('apply_bltable', varlist)
     2723
     2724            if insitu:
     2725                self._assign(workscan)
     2726                return res
     2727            else:
     2728                return {'scantable': workscan, 'fitresults': res}
     2729       
     2730        except RuntimeError, e:
     2731            raise_fitting_failure_exception(e)
     2732
     2733    @asaplog_post_dec
     2734    def sub_baseline(self, insitu=None, blinfo=None, bltable=None, overwrite=None):
     2735        """\
     2736        Subtract baseline based on parameters written in the input list.
     2737
     2738        Parameters:
     2739            insitu:        if False a new scantable is returned.
     2740                           Otherwise, the scaling is done in-situ
     2741                           The default is taken from .asaprc (False)
     2742            blinfo:        baseline parameter set stored in a dictionary
     2743                           or a list of dictionary. Each dictionary
     2744                           corresponds to each spectrum and must contain
     2745                           the following keys and values:
     2746                             'row': row number,
     2747                             'blfunc': function name. available ones include
     2748                                       'poly', 'chebyshev', 'cspline' and
     2749                                       'sinusoid',
     2750                             'order': maximum order of polynomial. needed
     2751                                      if blfunc='poly' or 'chebyshev',
     2752                             'npiece': number or piecewise polynomial.
     2753                                       needed if blfunc='cspline',
     2754                             'nwave': a list of sinusoidal wave numbers.
     2755                                      needed if blfunc='sinusoid', and
     2756                             'masklist': min-max windows for channel mask.
     2757                                         the specified ranges will be used
     2758                                         for fitting.
     2759            bltable:       name of output baseline table where baseline
     2760                           parameters and fitting results recorded.
     2761                           default is ''(no output).
     2762            overwrite:     if True, overwrites the existing baseline table
     2763                           specified in bltable.
     2764                           default is False.
     2765                           
     2766        Example:
     2767            sub_baseline(blinfo=[{'row':0, 'blfunc':'poly', 'order':5,
     2768                                  'masklist':[[10,350],[352,510]]},
     2769                                 {'row':1, 'blfunc':'cspline', 'npiece':3,
     2770                                  'masklist':[[3,16],[19,404],[407,511]]}
     2771                                  ])
     2772
     2773                the first spectrum (row=0) will be fitted with polynomial
     2774                of order=5 and the next one (row=1) will be fitted with cubic
     2775                spline consisting of 3 pieces.
     2776        """
     2777
     2778        try:
     2779            varlist = vars()
     2780            if blinfo         is None: blinfo         = []
     2781            if bltable        is None: bltable        = ''
     2782            if overwrite      is None: overwrite      = False
     2783
     2784            if insitu is None: insitu = rcParams['insitu']
     2785            if insitu:
     2786                workscan = self
     2787            else:
     2788                workscan = self.copy()
     2789
     2790            nrow = workscan.nrow()
     2791
     2792            in_blinfo = pack_blinfo(blinfo=blinfo, maxirow=nrow)
     2793
     2794            print "in_blinfo=< "+ str(in_blinfo)+" >"
     2795
     2796            sres = workscan._sub_baseline(in_blinfo,
     2797                                          bltable,
     2798                                          os.path.exists(bltable),
     2799                                          overwrite)
     2800            res = parse_fitresult(sres)
     2801           
     2802            workscan._add_history('sub_baseline', varlist)
     2803
     2804            if insitu:
     2805                self._assign(workscan)
     2806                return res
     2807            else:
     2808                return {'scantable': workscan, 'fitresults': res}
     2809
     2810        except RuntimeError, e:
     2811            raise_fitting_failure_exception(e)
     2812
     2813    @asaplog_post_dec
    25582814    def calc_aic(self, value=None, blfunc=None, order=None, mask=None,
    25592815                 whichrow=None, uselinefinder=None, edge=None,
     
    26382894                          clipniter=None, plot=None,
    26392895                          getresidual=None, showprogress=None,
    2640                           minnrow=None, outlog=None, blfile=None, csvformat=None):
     2896                          minnrow=None, outlog=None,
     2897                          blfile=None, csvformat=None,
     2898                          bltable=None):
    26412899        """\
    26422900        Return a scan which has been baselined (all rows) with sinusoidal
     
    26972955                           (default is '': no file/logger output)
    26982956            csvformat:     if True blfile is csv-formatted, default is False.
     2957            bltable:       name of a baseline table where fitting results
     2958                           (coefficients, rms, etc.) are to be written.
     2959                           if given, fitting results will NOT be output to
     2960                           scantable (insitu=True) or None will be
     2961                           returned (insitu=False).
     2962                           (default is "": no table output)
    26992963
    27002964        Example:
     
    27342998            if blfile        is None: blfile        = ''
    27352999            if csvformat     is None: csvformat     = False
    2736 
    2737             if csvformat:
    2738                 scsvformat = "T"
    2739             else:
    2740                 scsvformat = "F"
     3000            if bltable       is None: bltable       = ''
     3001
     3002            sapplyfft = 'true' if applyfft else 'false'
     3003            fftinfo = ','.join([sapplyfft, fftmethod.lower(), str(fftthresh).lower()])
     3004
     3005            scsvformat = 'T' if csvformat else 'F'
    27413006
    27423007            #CURRENTLY, PLOT=true is UNAVAILABLE UNTIL sinusoidal fitting is implemented as a fitter method.
    2743             workscan._sinusoid_baseline(mask, applyfft, fftmethod.lower(),
    2744                                         str(fftthresh).lower(),
     3008            workscan._sinusoid_baseline(mask,
     3009                                        fftinfo,
     3010                                        #applyfft, fftmethod.lower(),
     3011                                        #str(fftthresh).lower(),
    27453012                                        workscan._parse_wn(addwn),
    27463013                                        workscan._parse_wn(rejwn),
     
    27493016                                        pack_progress_params(showprogress,
    27503017                                                             minnrow),
    2751                                         outlog, scsvformat+blfile)
     3018                                        outlog, scsvformat+blfile,
     3019                                        bltable)
    27523020            workscan._add_history('sinusoid_baseline', varlist)
    2753            
    2754             if insitu:
    2755                 self._assign(workscan)
     3021
     3022            if bltable == '':
     3023                if insitu:
     3024                    self._assign(workscan)
     3025                else:
     3026                    return workscan
    27563027            else:
    2757                 return workscan
     3028                if not insitu:
     3029                    return None
    27583030           
    27593031        except RuntimeError, e:
     
    27683040                               chan_avg_limit=None, plot=None,
    27693041                               getresidual=None, showprogress=None,
    2770                                minnrow=None, outlog=None, blfile=None, csvformat=None):
     3042                               minnrow=None, outlog=None,
     3043                               blfile=None, csvformat=None,
     3044                               bltable=None):
    27713045        """\
    27723046        Return a scan which has been baselined (all rows) with sinusoidal
     
    28493123                            (default is "": no file/logger output)
    28503124            csvformat:      if True blfile is csv-formatted, default is False.
     3125            bltable:        name of a baseline table where fitting results
     3126                            (coefficients, rms, etc.) are to be written.
     3127                            if given, fitting results will NOT be output to
     3128                            scantable (insitu=True) or None will be
     3129                            returned (insitu=False).
     3130                            (default is "": no table output)
    28513131
    28523132        Example:
     
    28863166            if blfile         is None: blfile         = ''
    28873167            if csvformat      is None: csvformat      = False
    2888 
    2889             if csvformat:
    2890                 scsvformat = "T"
    2891             else:
    2892                 scsvformat = "F"
     3168            if bltable        is None: bltable        = ''
     3169
     3170            sapplyfft = 'true' if applyfft else 'false'
     3171            fftinfo = ','.join([sapplyfft, fftmethod.lower(), str(fftthresh).lower()])
     3172
     3173            scsvformat = 'T' if csvformat else 'F'
    28933174
    28943175            #CURRENTLY, PLOT=true is UNAVAILABLE UNTIL sinusoidal fitting is implemented as a fitter method.
    2895             workscan._auto_sinusoid_baseline(mask, applyfft,
    2896                                              fftmethod.lower(),
    2897                                              str(fftthresh).lower(),
     3176            workscan._auto_sinusoid_baseline(mask,
     3177                                             fftinfo,
    28983178                                             workscan._parse_wn(addwn),
    28993179                                             workscan._parse_wn(rejwn),
     
    29043184                                             pack_progress_params(showprogress,
    29053185                                                                  minnrow),
    2906                                              outlog, scsvformat+blfile)
     3186                                             outlog, scsvformat+blfile, bltable)
    29073187            workscan._add_history("auto_sinusoid_baseline", varlist)
    2908            
    2909             if insitu:
    2910                 self._assign(workscan)
     3188
     3189            if bltable == '':
     3190                if insitu:
     3191                    self._assign(workscan)
     3192                else:
     3193                    return workscan
    29113194            else:
    2912                 return workscan
     3195                if not insitu:
     3196                    return None
    29133197           
    29143198        except RuntimeError, e:
     
    29193203                         clipthresh=None, clipniter=None, plot=None,
    29203204                         getresidual=None, showprogress=None, minnrow=None,
    2921                          outlog=None, blfile=None, csvformat=None):
     3205                         outlog=None, blfile=None, csvformat=None,
     3206                         bltable=None):
    29223207        """\
    29233208        Return a scan which has been baselined (all rows) by cubic spline
     
    29493234                          (default is "": no file/logger output)
    29503235            csvformat:    if True blfile is csv-formatted, default is False.
     3236            bltable:      name of a baseline table where fitting results
     3237                          (coefficients, rms, etc.) are to be written.
     3238                          if given, fitting results will NOT be output to
     3239                          scantable (insitu=True) or None will be
     3240                          returned (insitu=False).
     3241                          (default is "": no table output)
    29513242
    29523243        Example:
     
    29813272            if outlog       is None: outlog       = False
    29823273            if blfile       is None: blfile       = ''
    2983             if csvformat     is None: csvformat     = False
    2984 
    2985             if csvformat:
    2986                 scsvformat = "T"
    2987             else:
    2988                 scsvformat = "F"
     3274            if csvformat    is None: csvformat    = False
     3275            if bltable      is None: bltable      = ''
     3276
     3277            scsvformat = 'T' if csvformat else 'F'
    29893278
    29903279            #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method.
    2991             workscan._cspline_baseline(mask, npiece, clipthresh, clipniter,
     3280            workscan._cspline_baseline(mask, npiece,
     3281                                       clipthresh, clipniter,
    29923282                                       getresidual,
    29933283                                       pack_progress_params(showprogress,
    29943284                                                            minnrow),
    2995                                        outlog, scsvformat+blfile)
     3285                                       outlog, scsvformat+blfile,
     3286                                       bltable)
    29963287            workscan._add_history("cspline_baseline", varlist)
    2997            
    2998             if insitu:
    2999                 self._assign(workscan)
     3288
     3289            if bltable == '':
     3290                if insitu:
     3291                    self._assign(workscan)
     3292                else:
     3293                    return workscan
    30003294            else:
    3001                 return workscan
     3295                if not insitu:
     3296                    return None
    30023297           
    30033298        except RuntimeError, e:
     
    30103305                              getresidual=None, plot=None,
    30113306                              showprogress=None, minnrow=None, outlog=None,
    3012                               blfile=None, csvformat=None):
     3307                              blfile=None, csvformat=None, bltable=None):
    30133308        """\
    30143309        Return a scan which has been baselined (all rows) by cubic spline
     
    30623357                            (default is "": no file/logger output)
    30633358            csvformat:      if True blfile is csv-formatted, default is False.
     3359            bltable:        name of a baseline table where fitting results
     3360                            (coefficients, rms, etc.) are to be written.
     3361                            if given, fitting results will NOT be output to
     3362                            scantable (insitu=True) or None will be
     3363                            returned (insitu=False).
     3364                            (default is "": no table output)
    30643365
    30653366        Example:
     
    30953396            if blfile         is None: blfile         = ''
    30963397            if csvformat      is None: csvformat      = False
    3097 
    3098             if csvformat:
    3099                 scsvformat = "T"
    3100             else:
    3101                 scsvformat = "F"
     3398            if bltable        is None: bltable        = ''
     3399
     3400            scsvformat = 'T' if csvformat else 'F'
    31023401
    31033402            #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method.
    3104             workscan._auto_cspline_baseline(mask, npiece, clipthresh,
    3105                                             clipniter,
     3403            workscan._auto_cspline_baseline(mask, npiece,
     3404                                            clipthresh, clipniter,
    31063405                                            normalise_edge_param(edge),
    31073406                                            threshold,
     
    31093408                                            pack_progress_params(showprogress,
    31103409                                                                 minnrow),
    3111                                             outlog, scsvformat+blfile)
     3410                                            outlog,
     3411                                            scsvformat+blfile,
     3412                                            bltable)
    31123413            workscan._add_history("auto_cspline_baseline", varlist)
    3113            
    3114             if insitu:
    3115                 self._assign(workscan)
     3414
     3415            if bltable == '':
     3416                if insitu:
     3417                    self._assign(workscan)
     3418                else:
     3419                    return workscan
    31163420            else:
    3117                 return workscan
     3421                if not insitu:
     3422                    return None
    31183423           
    31193424        except RuntimeError, e:
     
    31243429                           clipthresh=None, clipniter=None, plot=None,
    31253430                           getresidual=None, showprogress=None, minnrow=None,
    3126                            outlog=None, blfile=None, csvformat=None):
     3431                           outlog=None, blfile=None, csvformat=None,
     3432                           bltable=None):
    31273433        """\
    31283434        Return a scan which has been baselined (all rows) by Chebyshev polynomials.
     
    31533459                          (default is "": no file/logger output)
    31543460            csvformat:    if True blfile is csv-formatted, default is False.
     3461            bltable:      name of a baseline table where fitting results
     3462                          (coefficients, rms, etc.) are to be written.
     3463                          if given, fitting results will NOT be output to
     3464                          scantable (insitu=True) or None will be
     3465                          returned (insitu=False).
     3466                          (default is "": no table output)
    31553467
    31563468        Example:
     
    31743486                workscan = self.copy()
    31753487
    3176             #if mask         is None: mask         = [True for i in xrange(workscan.nchan())]
    31773488            if mask         is None: mask         = []
    31783489            if order        is None: order        = 5
     
    31853496            if outlog       is None: outlog       = False
    31863497            if blfile       is None: blfile       = ''
    3187             if csvformat     is None: csvformat     = False
    3188 
    3189             if csvformat:
    3190                 scsvformat = "T"
    3191             else:
    3192                 scsvformat = "F"
     3498            if csvformat    is None: csvformat    = False
     3499            if bltable      is None: bltable      = ''
     3500
     3501            scsvformat = 'T' if csvformat else 'F'
    31933502
    31943503            #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method.
    3195             workscan._chebyshev_baseline(mask, order, clipthresh, clipniter,
     3504            workscan._chebyshev_baseline(mask, order,
     3505                                         clipthresh, clipniter,
    31963506                                         getresidual,
    31973507                                         pack_progress_params(showprogress,
    31983508                                                              minnrow),
    3199                                          outlog, scsvformat+blfile)
     3509                                         outlog, scsvformat+blfile,
     3510                                         bltable)
    32003511            workscan._add_history("chebyshev_baseline", varlist)
    3201            
    3202             if insitu:
    3203                 self._assign(workscan)
     3512
     3513            if bltable == '':
     3514                if insitu:
     3515                    self._assign(workscan)
     3516                else:
     3517                    return workscan
    32043518            else:
    3205                 return workscan
     3519                if not insitu:
     3520                    return None
    32063521           
    32073522        except RuntimeError, e:
     
    32143529                              getresidual=None, plot=None,
    32153530                              showprogress=None, minnrow=None, outlog=None,
    3216                               blfile=None, csvformat=None):
     3531                              blfile=None, csvformat=None, bltable=None):
    32173532        """\
    32183533        Return a scan which has been baselined (all rows) by Chebyshev polynomials.
     
    32653580                            (default is "": no file/logger output)
    32663581            csvformat:      if True blfile is csv-formatted, default is False.
     3582            bltable:        name of a baseline table where fitting results
     3583                            (coefficients, rms, etc.) are to be written.
     3584                            if given, fitting results will NOT be output to
     3585                            scantable (insitu=True) or None will be
     3586                            returned (insitu=False).
     3587                            (default is "": no table output)
    32673588
    32683589        Example:
     
    32833604                workscan = self.copy()
    32843605           
    3285             #if mask           is None: mask           = [True for i in xrange(workscan.nchan())]
    32863606            if mask           is None: mask           = []
    32873607            if order          is None: order          = 5
     
    32983618            if blfile         is None: blfile         = ''
    32993619            if csvformat      is None: csvformat      = False
    3300 
    3301             if csvformat:
    3302                 scsvformat = "T"
    3303             else:
    3304                 scsvformat = "F"
     3620            if bltable        is None: bltable        = ''
     3621
     3622            scsvformat = 'T' if csvformat else 'F'
    33053623
    33063624            #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method.
    3307             workscan._auto_chebyshev_baseline(mask, order, clipthresh,
    3308                                               clipniter,
     3625            workscan._auto_chebyshev_baseline(mask, order,
     3626                                              clipthresh, clipniter,
    33093627                                              normalise_edge_param(edge),
    33103628                                              threshold,
     
    33123630                                              pack_progress_params(showprogress,
    33133631                                                                   minnrow),
    3314                                               outlog, scsvformat+blfile)
     3632                                              outlog, scsvformat+blfile,
     3633                                              bltable)
    33153634            workscan._add_history("auto_chebyshev_baseline", varlist)
    3316            
    3317             if insitu:
    3318                 self._assign(workscan)
     3635
     3636            if bltable == '':
     3637                if insitu:
     3638                    self._assign(workscan)
     3639                else:
     3640                    return workscan
    33193641            else:
    3320                 return workscan
     3642                if not insitu:
     3643                    return None
    33213644           
    33223645        except RuntimeError, e:
     
    33243647
    33253648    @asaplog_post_dec
    3326     def poly_baseline(self, mask=None, order=None, insitu=None, plot=None,
     3649    def poly_baseline(self, insitu=None, mask=None, order=None,
     3650                      clipthresh=None, clipniter=None, plot=None,
    33273651                      getresidual=None, showprogress=None, minnrow=None,
    3328                       outlog=None, blfile=None, csvformat=None):
     3652                      outlog=None, blfile=None, csvformat=None,
     3653                      bltable=None):
    33293654        """\
    33303655        Return a scan which has been baselined (all rows) by a polynomial.
     
    33353660            mask:         an optional mask
    33363661            order:        the order of the polynomial (default is 0)
     3662            clipthresh:   Clipping threshold. (default is 3.0, unit: sigma)
     3663            clipniter:    maximum number of iteration of 'clipthresh'-sigma
     3664                          clipping (default is 0)
    33373665            plot:         plot the fit and the residual. In this each
    33383666                          indivual fit has to be approved, by typing 'y'
     
    33503678                          (default is "": no file/logger output)
    33513679            csvformat:    if True blfile is csv-formatted, default is False.
     3680            bltable:      name of a baseline table where fitting results
     3681                          (coefficients, rms, etc.) are to be written.
     3682                          if given, fitting results will NOT be output to
     3683                          scantable (insitu=True) or None will be
     3684                          returned (insitu=False).
     3685                          (default is "": no table output)
    33523686
    33533687        Example:
     
    33673701                workscan = self.copy()
    33683702
    3369             #if mask         is None: mask         = [True for i in \
    3370             #                                           xrange(workscan.nchan())]
    33713703            if mask         is None: mask         = []
    33723704            if order        is None: order        = 0
     3705            if clipthresh   is None: clipthresh   = 3.0
     3706            if clipniter    is None: clipniter    = 0
    33733707            if plot         is None: plot         = False
    33743708            if getresidual  is None: getresidual  = True
     
    33763710            if minnrow      is None: minnrow      = 1000
    33773711            if outlog       is None: outlog       = False
    3378             if blfile       is None: blfile       = ""
     3712            if blfile       is None: blfile       = ''
    33793713            if csvformat    is None: csvformat    = False
    3380 
    3381             if csvformat:
    3382                 scsvformat = "T"
    3383             else:
    3384                 scsvformat = "F"
     3714            if bltable      is None: bltable      = ''
     3715
     3716            scsvformat = 'T' if csvformat else 'F'
    33853717
    33863718            if plot:
     
    34463778                    blf.close()
    34473779            else:
    3448                 workscan._poly_baseline(mask, order, getresidual,
     3780                workscan._poly_baseline(mask, order,
     3781                                        clipthresh, clipniter, #
     3782                                        getresidual,
    34493783                                        pack_progress_params(showprogress,
    34503784                                                             minnrow),
    3451                                         outlog, scsvformat+blfile)
     3785                                        outlog, scsvformat+blfile,
     3786                                        bltable)  #
    34523787           
    34533788            workscan._add_history("poly_baseline", varlist)
     
    34623797
    34633798    @asaplog_post_dec
    3464     def auto_poly_baseline(self, mask=None, order=None, edge=None,
    3465                            threshold=None, chan_avg_limit=None,
    3466                            plot=None, insitu=None,
    3467                            getresidual=None, showprogress=None,
    3468                            minnrow=None, outlog=None, blfile=None, csvformat=None):
     3799    def auto_poly_baseline(self, insitu=None, mask=None, order=None,
     3800                           clipthresh=None, clipniter=None,
     3801                           edge=None, threshold=None, chan_avg_limit=None,
     3802                           getresidual=None, plot=None,
     3803                           showprogress=None, minnrow=None, outlog=None,
     3804                           blfile=None, csvformat=None, bltable=None):
    34693805        """\
    34703806        Return a scan which has been baselined (all rows) by a polynomial.
     
    34783814            mask:           an optional mask retreived from scantable
    34793815            order:          the order of the polynomial (default is 0)
     3816            clipthresh:     Clipping threshold. (default is 3.0, unit: sigma)
     3817            clipniter:      maximum number of iteration of 'clipthresh'-sigma
     3818                            clipping (default is 0)
    34803819            edge:           an optional number of channel to drop at
    34813820                            the edge of spectrum. If only one value is
     
    35133852                            (default is "": no file/logger output)
    35143853            csvformat:      if True blfile is csv-formatted, default is False.
     3854            bltable:        name of a baseline table where fitting results
     3855                            (coefficients, rms, etc.) are to be written.
     3856                            if given, fitting results will NOT be output to
     3857                            scantable (insitu=True) or None will be
     3858                            returned (insitu=False).
     3859                            (default is "": no table output)
    35153860
    35163861        Example:
     
    35283873                workscan = self.copy()
    35293874
    3530             #if mask           is None: mask           = [True for i in xrange(workscan.nchan())]
    35313875            if mask           is None: mask           = []
    35323876            if order          is None: order          = 0
     3877            if clipthresh     is None: clipthresh     = 3.0
     3878            if clipniter      is None: clipniter      = 0
    35333879            if edge           is None: edge           = (0, 0)
    35343880            if threshold      is None: threshold      = 3
     
    35413887            if blfile         is None: blfile         = ''
    35423888            if csvformat      is None: csvformat      = False
    3543 
    3544             if csvformat:
    3545                 scsvformat = "T"
    3546             else:
    3547                 scsvformat = "F"
     3889            if bltable        is None: bltable        = ''
     3890
     3891            scsvformat = 'T' if csvformat else 'F'
    35483892
    35493893            edge = normalise_edge_param(edge)
     
    36153959                if outblfile: blf.close()
    36163960            else:
    3617                 workscan._auto_poly_baseline(mask, order, edge, threshold,
     3961                workscan._auto_poly_baseline(mask, order,
     3962                                             clipthresh, clipniter,
     3963                                             edge, threshold,
    36183964                                             chan_avg_limit, getresidual,
    36193965                                             pack_progress_params(showprogress,
    36203966                                                                  minnrow),
    3621                                              outlog, scsvformat+blfile)
    3622 
     3967                                             outlog, scsvformat+blfile,
     3968                                             bltable)
    36233969            workscan._add_history("auto_poly_baseline", varlist)
    3624            
    3625             if insitu:
    3626                 self._assign(workscan)
     3970
     3971            if bltable == '':
     3972                if insitu:
     3973                    self._assign(workscan)
     3974                else:
     3975                    return workscan
    36273976            else:
    3628                 return workscan
     3977                if not insitu:
     3978                    return None
    36293979           
    36303980        except RuntimeError, e:
  • trunk/src/STBaselineEnum.h

    r2737 r2767  
    2828class STBaselineFunc  {
    2929public:
    30   enum FuncName {Polynomial = 0,
     30  enum FuncName {Polynomial = 1,
    3131                 CSpline,
    3232                 Sinusoid,
  • trunk/src/STBaselineTable.cpp

    r2740 r2767  
    5151void STBaselineTable::setup()
    5252{
    53   table_.addColumn(ScalarColumnDesc<uInt>("NCHAN"));
     53  table_.addColumn(ScalarColumnDesc<Bool>("APPLY"));
    5454  table_.addColumn(ScalarColumnDesc<uInt>("FUNC_TYPE"));
    55   table_.addColumn(ArrayColumnDesc<uInt>("FUNC_PARAM"));
     55  table_.addColumn(ArrayColumnDesc<Int>("FUNC_PARAM"));
    5656  table_.addColumn(ArrayColumnDesc<Float>("FUNC_FPARAM"));
    57   table_.addColumn(ScalarColumnDesc<uInt>("CLIP_ITERATION"));
    58   table_.addColumn(ScalarColumnDesc<Float>("CLIP_THRESHOLD"));
    5957  table_.addColumn(ArrayColumnDesc<uInt>("MASKLIST"));
    6058  table_.addColumn(ArrayColumnDesc<Float>("RESULT"));
    6159  table_.addColumn(ScalarColumnDesc<Float>("RMS"));
     60  table_.addColumn(ScalarColumnDesc<uInt>("NCHAN"));
     61  table_.addColumn(ScalarColumnDesc<Float>("CLIP_THRESHOLD"));
     62  table_.addColumn(ScalarColumnDesc<uInt>("CLIP_ITERATION"));
     63  table_.addColumn(ScalarColumnDesc<Float>("LF_THRESHOLD"));
     64  table_.addColumn(ScalarColumnDesc<uInt>("LF_AVERAGE"));
     65  table_.addColumn(ArrayColumnDesc<uInt>("LF_EDGE"));
    6266
    6367  table_.rwKeywordSet().define("ApplyType", "BASELINE");
     
    6872void STBaselineTable::attachOptionalColumns()
    6973{
    70   nchanCol_.attach(table_, "NCHAN");
     74  applyCol_.attach(table_, "APPLY");
    7175  ftypeCol_.attach(table_, "FUNC_TYPE");
    7276  fparCol_.attach(table_, "FUNC_PARAM");
    7377  ffparCol_.attach(table_, "FUNC_FPARAM");
    74   citerCol_.attach(table_, "CLIP_ITERATION");
    75   cthresCol_.attach(table_, "CLIP_THRESHOLD");
    7678  maskCol_.attach(table_, "MASKLIST");
    7779  resCol_.attach(table_, "RESULT");
    7880  rmsCol_.attach(table_, "RMS");
     81  nchanCol_.attach(table_, "NCHAN");
     82  cthresCol_.attach(table_, "CLIP_THRESHOLD");
     83  citerCol_.attach(table_, "CLIP_ITERATION");
     84  lfthresCol_.attach(table_, "LF_THRESHOLD");
     85  lfavgCol_.attach(table_, "LF_AVERAGE");
     86  lfedgeCol_.attach(table_, "LF_EDGE");
    7987}
    8088
     
    9098                              uInt beamno, uInt ifno, uInt polno,
    9199                              uInt freqid, Double time,
    92                               uInt nchan,
     100                              Bool apply,
    93101                              STBaselineFunc::FuncName ftype,
    94                               Vector<uInt> fpar,
     102                              Vector<Int> fpar,
    95103                              Vector<Float> ffpar,
    96                               uInt citer,
    97                               Float cthres,
    98104                              Vector<uInt> mask,
    99105                              Vector<Float> res,
    100                               Float rms)
     106                              Float rms,
     107                              uInt nchan,
     108                              Float cthres,
     109                              uInt citer,
     110                              Float lfthres,
     111                              uInt lfavg,
     112                              Vector<uInt> lfedge)
    101113{
    102114  if (irow >= (uInt)nrow()) {
     
    110122
    111123  setbasedata(irow, scanno, cycleno, beamno, ifno, polno, freqid, time);
    112   nchanCol_.put(irow, nchan);
     124  applyCol_.put(irow, apply);
    113125  ftypeCol_.put(irow, uInt(ftype));
    114126  fparCol_.put(irow, fpar);
    115127  ffparCol_.put(irow, ffpar);
    116   citerCol_.put(irow, citer);
    117   cthresCol_.put(irow, cthres);
    118128  maskCol_.put(irow, mask);
    119129  resCol_.put(irow, res);
    120130  rmsCol_.put(irow, rms);
     131  nchanCol_.put(irow, nchan);
     132  cthresCol_.put(irow, cthres);
     133  citerCol_.put(irow, citer);
     134  lfthresCol_.put(irow, lfthres);
     135  lfavgCol_.put(irow, lfavg);
     136  lfedgeCol_.put(irow, lfedge);
    121137}
    122138
     
    124140                                 uInt beamno, uInt ifno, uInt polno,
    125141                                 uInt freqid, Double time,
    126                                  uInt nchan,
     142                                 Bool apply,
    127143                                 STBaselineFunc::FuncName ftype,
    128                                  Vector<uInt> fpar,
     144                                 Vector<Int> fpar,
    129145                                 Vector<Float> ffpar,
    130                                  uInt citer,
    131                                  Float cthres,
    132146                                 Vector<uInt> mask,
    133147                                 Vector<Float> res,
    134                                  Float rms)
     148                                 Float rms,
     149                                 uInt nchan,
     150                                 Float cthres,
     151                                 uInt citer,
     152                                 Float lfthres,
     153                                 uInt lfavg,
     154                                 Vector<uInt> lfedge)
    135155{
    136156  uInt irow = nrow();
    137157  table_.addRow(1, True);
    138158  setdata(irow, scanno, cycleno, beamno, ifno, polno, freqid, time,
    139           nchan, ftype, fpar, ffpar, citer, cthres, mask, res, rms);
    140 }
    141 
    142 Vector<STBaselineFunc::FuncName> STBaselineTable::getFunctionAsString()
     159          apply, ftype, fpar, ffpar, mask, res, rms,
     160          nchan, cthres, citer, lfthres, lfavg, lfedge);
     161}
     162
     163void STBaselineTable::appendbasedata(int scanno, int cycleno,
     164                                     int beamno, int ifno, int polno,
     165                                     int freqid, Double time)
     166{
     167  uInt irow = nrow();
     168  table_.addRow(1, True);
     169  setbasedata(irow, uInt(scanno), uInt(cycleno), uInt(beamno), uInt(ifno), uInt(polno), uInt(freqid), time);
     170}
     171
     172void STBaselineTable::setresult(casa::uInt irow,
     173                                casa::Vector<casa::Float> res,
     174                                casa::Float rms)
     175{
     176  resCol_.put(irow, res);
     177  rmsCol_.put(irow, rms);
     178}
     179
     180bool STBaselineTable::getApply(int irow)
     181{
     182  return (bool)applyCol_.get(irow);
     183}
     184
     185void STBaselineTable::setApply(int irow, bool apply)
     186{
     187  applyCol_.put(uInt(irow), Bool(apply));
     188}
     189
     190Vector<STBaselineFunc::FuncName> STBaselineTable::getFunctionNames()
    143191{
    144192  Vector<uInt> rawBlfuncColumn = ftypeCol_.getColumn();
     
    149197  }
    150198  return blfuncColumn;
     199}
     200
     201STBaselineFunc::FuncName STBaselineTable::getFunctionName(int irow)
     202{
     203  return STBaselineFunc::FuncName(ftypeCol_.get(irow));
     204}
     205
     206std::vector<int> STBaselineTable::getFuncParam(int irow)
     207{
     208  Vector<Int> uiparam = fparCol_.get(irow);
     209  std::vector<int> res(uiparam.size());
     210  for (uInt i = 0; i < res.size(); ++i) {
     211    res[i] = (int)uiparam[i];
     212  }
     213  return res;
     214}
     215
     216std::vector<bool> STBaselineTable::getMask(int irow)
     217{
     218  uInt nchan = getNChan(irow);
     219  Vector<uInt> masklist = maskCol_.get(irow);
     220  std::vector<int> masklist1(masklist.size());
     221  for (uInt i = 0; i < masklist1.size(); ++i) {
     222    masklist1[i] = (int)masklist[i];
     223  }
     224  return Scantable::getMaskFromMaskList(nchan, masklist1);
     225}
     226
     227uInt STBaselineTable::getNChan(int irow)
     228{
     229  return nchanCol_.get(irow);
    151230}
    152231
  • trunk/src/STBaselineTable.h

    r2740 r2767  
    4646               casa::uInt beamno, casa::uInt ifno, casa::uInt polno,
    4747               casa::uInt freqid, casa::Double time,
    48                casa::uInt nchan,
     48               casa::Bool apply,
    4949               STBaselineFunc::FuncName ftype,
    50                casa::Vector<casa::uInt> fpar,
     50               casa::Vector<casa::Int> fpar,
    5151               casa::Vector<casa::Float> ffpar,
    52                casa::uInt citer,
    53                casa::Float cthres,
    5452               casa::Vector<casa::uInt> mask,
    5553               casa::Vector<casa::Float> res,
    56                casa::Float rms);
     54               casa::Float rms,
     55               casa::uInt nchan,
     56               casa::Float cthres,
     57               casa::uInt citer,
     58               casa::Float lfthres,
     59               casa::uInt lfavg,
     60               casa::Vector<casa::uInt> lfedge);
    5761  void appenddata(casa::uInt scanno, casa::uInt cycleno,
    5862                  casa::uInt beamno, casa::uInt ifno, casa::uInt polno,
    5963                  casa::uInt freqid, casa::Double time,
    60                   casa::uInt nchan,
     64                  casa::Bool apply,
    6165                  STBaselineFunc::FuncName ftype,
    62                   casa::Vector<casa::uInt> fpar,
     66                  casa::Vector<casa::Int> fpar,
    6367                  casa::Vector<casa::Float> ffpar,
    64                   casa::uInt citer,
    65                   casa::Float cthres,
    6668                  casa::Vector<casa::uInt> mask,
    6769                  casa::Vector<casa::Float> res,
    68                   casa::Float rms);
     70                  casa::Float rms,
     71                  casa::uInt nchan,
     72                  casa::Float cthres,
     73                  casa::uInt citer,
     74                  casa::Float lfthres,
     75                  casa::uInt lfavg,
     76                  casa::Vector<casa::uInt> lfedge);
     77  void appendbasedata(int scanno, int cycleno,
     78                      int beamno, int ifno, int polno,
     79                      int freqid, casa::Double time);
     80  void setresult(casa::uInt irow,
     81                 casa::Vector<casa::Float> res,
     82                 casa::Float rms);
    6983  casa::uInt nchan(casa::uInt ifno);
     84  casa::Vector<casa::Bool> getApply() {return applyCol_.getColumn();}
     85  bool getApply(int irow);
    7086  casa::Vector<casa::uInt> getFunction() {return ftypeCol_.getColumn();}
    71   casa::Vector<STBaselineFunc::FuncName> getFunctionAsString();
    72   casa::Matrix<casa::uInt> getFuncParam() {return fparCol_.getColumn();}
     87  casa::Vector<STBaselineFunc::FuncName> getFunctionNames();
     88  STBaselineFunc::FuncName getFunctionName(int irow);
     89  casa::Matrix<casa::Int> getFuncParam() {return fparCol_.getColumn();}
     90  std::vector<int> getFuncParam(int irow);
    7391  casa::Matrix<casa::Float> getFuncFParam() {return ffparCol_.getColumn();}
    74   casa::Vector<casa::uInt> getClipIteration() {return citerCol_.getColumn();}
    75   casa::Vector<casa::Float> getClipThreshold() {return cthresCol_.getColumn();}
    7692  casa::Matrix<casa::uInt> getMaskList() {return maskCol_.getColumn();}
     93  std::vector<bool> getMask(int irow);
    7794  casa::Matrix<casa::Float> getResult() {return resCol_.getColumn();}
    7895  casa::Vector<casa::Float> getRms() {return rmsCol_.getColumn();}
     96  casa::Vector<casa::uInt> getNChan() {return nchanCol_.getColumn();}
     97  casa::uInt getNChan(int irow);
     98  casa::Vector<casa::Float> getClipThreshold() {return cthresCol_.getColumn();}
     99  casa::Vector<casa::uInt> getClipIteration() {return citerCol_.getColumn();}
     100  casa::Vector<casa::Float> getLineFinderThreshold() {return lfthresCol_.getColumn();}
     101  casa::Vector<casa::uInt> getLineFinderChanAvg() {return lfavgCol_.getColumn();}
     102  casa::Matrix<casa::uInt> getLineFinderEdge() {return lfedgeCol_.getColumn();}
     103  void setApply(int irow, bool apply);
    79104
    80105private:
    81106  static const casa::String name_ ;
    82   casa::ScalarColumn<casa::uInt> nchanCol_;
     107  casa::ScalarColumn<casa::Bool> applyCol_;
    83108  casa::ScalarColumn<casa::uInt> ftypeCol_;
    84   casa::ArrayColumn<casa::uInt> fparCol_;
     109  casa::ArrayColumn<casa::Int> fparCol_;
    85110  casa::ArrayColumn<casa::Float> ffparCol_;
    86   casa::ScalarColumn<casa::uInt> citerCol_;
    87   casa::ScalarColumn<casa::Float> cthresCol_;
    88111  casa::ArrayColumn<casa::uInt> maskCol_;
    89112  casa::ArrayColumn<casa::Float> resCol_;
    90113  casa::ScalarColumn<casa::Float> rmsCol_;
     114  casa::ScalarColumn<casa::uInt> nchanCol_;
     115  casa::ScalarColumn<casa::Float> cthresCol_;
     116  casa::ScalarColumn<casa::uInt> citerCol_;
     117  casa::ScalarColumn<casa::Float> lfthresCol_;
     118  casa::ScalarColumn<casa::uInt> lfavgCol_;
     119  casa::ArrayColumn<casa::uInt> lfedgeCol_;
    91120};
    92121
  • trunk/src/Scantable.cpp

    r2740 r2767  
    24902490}
    24912491
    2492 void Scantable::polyBaseline(const std::vector<bool>& mask, int order, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
    2493 {
     2492std::vector<std::string> Scantable::applyBaselineTable(const std::string& bltable, const std::string& outbltable, const bool outbltableexists, const bool overwrite)
     2493{
     2494  STBaselineTable* btin = new STBaselineTable(bltable);
     2495
     2496  Vector<Bool> applyCol = btin->getApply();
     2497  int nRowBl = applyCol.size();
     2498  if (nRowBl != nrow()) {
     2499    throw(AipsError("Scantable and bltable have different number of rows."));
     2500  }
     2501
     2502  std::vector<std::string> res;
     2503  res.clear();
     2504
     2505  bool outBaselineTable = ((outbltable != "") && (!outbltableexists || overwrite));
     2506  bool bltableidentical = (bltable == outbltable);
     2507  STBaselineTable* btout = new STBaselineTable(*this);
     2508  ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(table_, "TIME");
     2509  Vector<Double> timeSecCol = tcol->getColumn();
     2510
     2511  for (int whichrow = 0; whichrow < nRowBl; ++whichrow) {
     2512    if (applyCol[whichrow]) {
     2513      std::vector<float> spec = getSpectrum(whichrow);
     2514
     2515      //--channel mask--  maybe the following choices should be available
     2516      std::vector<bool> mask = btin->getMask(whichrow);  //use mask_bltable only (default)
     2517      //std::vector<bool> mask = getMask(whichrow);    // or mask_scantable only
     2518      //std::vector<bool> mask = getCompositeChanMask(whichrow, btin->getMask(whichrow)); // or mask_scantable && mask_bltable
     2519
     2520      STBaselineFunc::FuncName ftype = btin->getFunctionName(whichrow);
     2521      std::vector<int> fpar = btin->getFuncParam(whichrow);
     2522      std::vector<float> params;
     2523      float rms;
     2524      std::vector<float> resfit = doApplyBaselineTable(spec, mask, ftype, fpar, params, rms);
     2525
     2526      setSpectrum(resfit, whichrow);
     2527
     2528      res.push_back(packFittingResults(whichrow, params, rms));
     2529
     2530      if (outBaselineTable) {
     2531        Vector<Int> fparam(fpar.size());
     2532        for (uInt j = 0; j < fparam.size(); ++j) {
     2533          fparam[j] = (Int)fpar[j];
     2534        }
     2535        std::vector<int> maskList0;
     2536        maskList0.clear();
     2537        for (uInt j = 0; j < mask.size(); ++j) {
     2538          if (mask[j]) {
     2539            if ((j == 0)||(j == mask.size()-1)) {
     2540              maskList0.push_back(j);
     2541            } else {
     2542              if ((mask[j])&&(!mask[j-1])) {
     2543                maskList0.push_back(j);
     2544              }
     2545              if ((mask[j])&&(!mask[j+1])) {
     2546                maskList0.push_back(j);
     2547              }
     2548            }
     2549          }
     2550        }
     2551        Vector<uInt> maskList(maskList0.size());
     2552        for (uInt j = 0; j < maskList0.size(); ++j) {
     2553          maskList[j] = (uInt)maskList0[j];
     2554        }
     2555
     2556        if (outbltableexists) {
     2557          if (overwrite) {
     2558            if (bltableidentical) {
     2559              btin->setresult(uInt(whichrow), Vector<Float>(params), Float(rms));
     2560            } else {
     2561              btout->setresult(uInt(whichrow), Vector<Float>(params), Float(rms));
     2562            }
     2563          }
     2564        } else {
     2565          btout->appenddata(uInt(getScan(whichrow)), uInt(getCycle(whichrow)),
     2566                            uInt(getBeam(whichrow)), uInt(getIF(whichrow)), uInt(getPol(whichrow)),
     2567                            uInt(0), timeSecCol[whichrow], Bool(true), ftype, fparam,
     2568                            Vector<Float>(), maskList, Vector<Float>(params),
     2569                            Float(rms), uInt(spec.size()), Float(3.0), uInt(0),
     2570                            Float(0.0), uInt(0), Vector<uInt>());
     2571        }
     2572      }
     2573    }
     2574  }
     2575
     2576  if (outBaselineTable) {
     2577    if (bltableidentical) {
     2578      btin->save(outbltable);
     2579    } else {
     2580      btout->save(outbltable);
     2581    }
     2582  }
     2583  delete btin;
     2584  delete btout;
     2585
     2586  return res;
     2587}
     2588
     2589std::vector<std::string> Scantable::subBaseline(const std::vector<std::string>& blInfoList, const std::string& outbltable, const bool outbltableexists, const bool overwrite)
     2590{
     2591  int nRowBl = blInfoList.size();
     2592  int nRowSt = nrow();
     2593
     2594  std::vector<std::string> res;
     2595  res.clear();
     2596
     2597  bool outBaselineTable = ((outbltable != "") && (!outbltableexists || overwrite));
     2598  if ((outbltable != "") && outbltableexists && !overwrite) {
     2599    throw(AipsError("Cannot overwrite bltable. Set overwrite=True."));
     2600  }
     2601
     2602  STBaselineTable* bt = new STBaselineTable(*this);
     2603  ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(table_, "TIME");
     2604  Vector<Double> timeSecCol = tcol->getColumn();
     2605
     2606  if (outBaselineTable && !outbltableexists) {
     2607    for (int i = 0; i < nRowSt; ++i) {
     2608      bt->appendbasedata(getScan(i), getCycle(i), getBeam(i), getIF(i), getPol(i),
     2609                         0, timeSecCol[i]);
     2610      bt->setApply(i, false);
     2611    }
     2612  }
     2613
     2614  for (int i = 0; i < nRowBl; ++i) {
     2615    int irow;
     2616    STBaselineFunc::FuncName ftype;
     2617    std::vector<bool> mask;
     2618    std::vector<int> fpar;
     2619    float clipth;
     2620    int clipn;
     2621    bool uself;
     2622    float lfth;
     2623    std::vector<int> lfedge;
     2624    int lfavg;
     2625    parseBlInfo(blInfoList[i], irow, ftype, fpar, mask, clipth, clipn, uself, lfth, lfedge, lfavg);
     2626
     2627    if (irow < nRowSt) {
     2628      std::vector<float> spec = getSpectrum(irow);
     2629      std::vector<float> params;
     2630      float rms;
     2631      std::vector<bool> finalmask;
     2632
     2633      std::vector<float> resfit = doSubtractBaseline(spec, mask, ftype, fpar, params, rms, finalmask, clipth, clipn, uself, irow, lfth, lfedge, lfavg);
     2634      setSpectrum(resfit, irow);
     2635      res.push_back(packFittingResults(irow, params, rms));
     2636
     2637      if (outBaselineTable) {
     2638        Vector<Int> fparam(fpar.size());
     2639        for (uInt j = 0; j < fparam.size(); ++j) {
     2640          fparam[j] = (Int)fpar[j];
     2641        }
     2642        std::vector<int> maskList0;
     2643        maskList0.clear();
     2644        for (uInt j = 0; j < finalmask.size(); ++j) {
     2645          if (finalmask[j]) {
     2646            if ((j == 0)||(j == finalmask.size()-1)) {
     2647              maskList0.push_back(j);
     2648            } else {
     2649              if ((finalmask[j])&&(!finalmask[j-1])) {
     2650                maskList0.push_back(j);
     2651              }
     2652              if ((finalmask[j])&&(!finalmask[j+1])) {
     2653                maskList0.push_back(j);
     2654              }
     2655            }
     2656          }
     2657        }
     2658        Vector<uInt> maskList(maskList0.size());
     2659        for (uInt j = 0; j < maskList0.size(); ++j) {
     2660          maskList[j] = (uInt)maskList0[j];
     2661        }
     2662
     2663        bt->setdata(uInt(irow),
     2664                    uInt(getScan(irow)), uInt(getCycle(irow)),
     2665                    uInt(getBeam(irow)), uInt(getIF(irow)), uInt(getPol(irow)),
     2666                    uInt(0), timeSecCol[irow], Bool(true), ftype, fparam,
     2667                    Vector<Float>(), maskList, Vector<Float>(params),
     2668                    Float(rms), uInt(spec.size()), Float(clipth), uInt(clipn),
     2669                    Float(0.0), uInt(0), Vector<uInt>());
     2670      }
     2671
     2672    }
     2673  }
     2674
     2675  if (outBaselineTable) {
     2676    bt->save(outbltable);
     2677  }
     2678  delete bt;
     2679
     2680  return res;
     2681}
     2682
     2683std::vector<float> Scantable::doApplyBaselineTable(std::vector<float>& spec, std::vector<bool>& mask, const STBaselineFunc::FuncName ftype, std::vector<int>& fpar, std::vector<float>& params, float&rms)
     2684{
     2685  std::vector<bool> finalmask;
     2686  std::vector<int> lfedge;
     2687  return doSubtractBaseline(spec, mask, ftype, fpar, params, rms, finalmask, 0.0, 0, false, 0, 0.0, lfedge, 0);
     2688}
     2689
     2690  std::vector<float> Scantable::doSubtractBaseline(std::vector<float>& spec, std::vector<bool>& mask, const STBaselineFunc::FuncName ftype, std::vector<int>& fpar, std::vector<float>& params, float&rms, std::vector<bool>& finalmask, float clipth, int clipn, bool uself, int irow, float lfth, std::vector<int>& lfedge, int lfavg)
     2691{
     2692  //---replace mask with (mask AND mask_linefinder).
     2693  if (uself) {
     2694    STLineFinder lineFinder = STLineFinder();
     2695    lineFinder.setOptions(lfth, 3, lfavg);
     2696
     2697    int minEdgeSize = getIFNos().size()*2;
     2698    int edgeSize = lfedge.size();
     2699    std::vector<int> currentEdge;
     2700    currentEdge.clear();
     2701    if (edgeSize >= 2) {
     2702      int idx = 0;
     2703      if (edgeSize > 2) {
     2704        if (edgeSize < minEdgeSize) {
     2705          throw(AipsError("Length of edge element info is less than that of IFs"));
     2706        }
     2707        idx = 2 * getIF(irow);
     2708      }
     2709      currentEdge.push_back(lfedge[idx]);
     2710      currentEdge.push_back(lfedge[idx+1]);
     2711    } else {
     2712      throw(AipsError("Wrong length of edge element"));
     2713    }
     2714    lineFinder.setData(spec);
     2715    lineFinder.findLines(getCompositeChanMask(irow, mask), currentEdge, irow);
     2716
     2717    std::vector<bool> lfcompmask = lineFinder.getMask();
     2718    mask.clear();
     2719    mask.resize(lfcompmask.size());
     2720    for (uInt i = 0; i < mask.size(); ++i) {
     2721      mask[i] = lfcompmask[i];
     2722    }
     2723  }
     2724  //---
     2725
     2726  std::vector<float> res;
     2727  if (ftype == STBaselineFunc::Polynomial) {
     2728    res = doPolynomialFitting(spec, mask, fpar[0], params, rms, finalmask, clipth, clipn);
     2729  } else if (ftype == STBaselineFunc::Chebyshev) {
     2730    res = doChebyshevFitting(spec, mask, fpar[0], params, rms, finalmask, clipth, clipn);
     2731  } else if (ftype == STBaselineFunc::CSpline) {
     2732    if (fpar.size() > 1) { // reading from baseline table in which pieceEdges are already calculated and stored.
     2733      res = doCubicSplineFitting(spec, mask, fpar, params, rms, finalmask, clipth, clipn);
     2734    } else {               // usual cspline fitting by giving nPiece only. fpar will be replaced with pieceEdges.
     2735      res = doCubicSplineFitting(spec, mask, fpar[0], fpar, params, rms, finalmask, clipth, clipn);
     2736    }
     2737  } else if (ftype == STBaselineFunc::Sinusoid) {
     2738    res = doSinusoidFitting(spec, mask, fpar, params, rms, finalmask, clipth, clipn);
     2739  }
     2740
     2741  return res;
     2742}
     2743
     2744std::string Scantable::packFittingResults(const int irow, const std::vector<float>& params, const float rms)
     2745{
     2746  // returned value: "irow:params[0],params[1],..,params[n-1]:rms"
     2747  ostringstream os;
     2748  os << irow << ':';
     2749  for (uInt i = 0; i < params.size(); ++i) {
     2750    if (i > 0) {
     2751      os << ',';
     2752    }
     2753    os << params[i];
     2754  }
     2755  os << ':' << rms;
     2756
     2757  return os.str();
     2758}
     2759
     2760void Scantable::parseBlInfo(const std::string& blInfo, int& irow, STBaselineFunc::FuncName& ftype, std::vector<int>& fpar, std::vector<bool>& mask, float& thresClip, int& nIterClip, bool& useLineFinder, float& thresLF, std::vector<int>& edgeLF, int& avgLF)
     2761{
     2762  // The baseline info to be parsed must be column-delimited string like
     2763  // "0:chebyshev:5:3,5,169,174,485,487" where the elements are
     2764  // row number, funcType, funcOrder, maskList, clipThreshold, clipNIter,
     2765  // useLineFinder, lfThreshold, lfEdge and lfChanAvgLimit.
     2766
     2767  std::vector<string> res = splitToStringList(blInfo, ':');
     2768  if (res.size() < 4) {
     2769    throw(AipsError("baseline info has bad format")) ;
     2770  }
     2771
     2772  string ftype0, fpar0, masklist0, uself0, edge0;
     2773  std::vector<int> masklist;
     2774
     2775  stringstream ss;
     2776  ss << res[0];
     2777  ss >> irow;
     2778  ss.clear(); ss.str("");
     2779
     2780  ss << res[1];
     2781  ss >> ftype0;
     2782  if (ftype0 == "poly") {
     2783    ftype = STBaselineFunc::Polynomial;
     2784  } else if (ftype0 == "cspline") {
     2785    ftype = STBaselineFunc::CSpline;
     2786  } else if (ftype0 == "sinusoid") {
     2787    ftype = STBaselineFunc::Sinusoid;
     2788  } else if (ftype0 == "chebyshev") {
     2789    ftype = STBaselineFunc::Chebyshev;
     2790  } else {
     2791    throw(AipsError("invalid function type."));
     2792  }
     2793  ss.clear(); ss.str("");
     2794
     2795  ss << res[2];
     2796  ss >> fpar0;
     2797  fpar = splitToIntList(fpar0, ',');
     2798  ss.clear(); ss.str("");
     2799
     2800  ss << res[3];
     2801  ss >> masklist0;
     2802  mask = getMaskFromMaskList(nchan(getIF(irow)), splitToIntList(masklist0, ','));
     2803
     2804  ss << res[4];
     2805  ss >> thresClip;
     2806  ss.clear(); ss.str("");
     2807
     2808  ss << res[5];
     2809  ss >> nIterClip;
     2810  ss.clear(); ss.str("");
     2811
     2812  ss << res[6];
     2813  ss >> uself0;
     2814  if (uself0 == "true") {
     2815    useLineFinder = true;
     2816  } else {
     2817    useLineFinder = false;
     2818  }
     2819  ss.clear(); ss.str("");
     2820
     2821  if (useLineFinder) {
     2822    ss << res[7];
     2823    ss >> thresLF;
     2824    ss.clear(); ss.str("");
     2825
     2826    ss << res[8];
     2827    ss >> edge0;
     2828    edgeLF = splitToIntList(edge0, ',');
     2829    ss.clear(); ss.str("");
     2830
     2831    ss << res[9];
     2832    ss >> avgLF;
     2833    ss.clear(); ss.str("");
     2834  }
     2835
     2836}
     2837
     2838std::vector<int> Scantable::splitToIntList(const std::string& s, const char delim)
     2839{
     2840  istringstream iss(s);
     2841  string tmp;
     2842  int tmpi;
     2843  std::vector<int> res;
     2844  stringstream ss;
     2845  while (getline(iss, tmp, delim)) {
     2846    ss << tmp;
     2847    ss >> tmpi;
     2848    res.push_back(tmpi);
     2849    ss.clear(); ss.str("");
     2850  }
     2851
     2852  return res;
     2853}
     2854
     2855std::vector<string> Scantable::splitToStringList(const std::string& s, const char delim)
     2856{
     2857  istringstream iss(s);
     2858  std::string tmp;
     2859  std::vector<string> res;
     2860  while (getline(iss, tmp, delim)) {
     2861    res.push_back(tmp);
     2862  }
     2863
     2864  return res;
     2865}
     2866
     2867std::vector<bool> Scantable::getMaskFromMaskList(const int nchan, const std::vector<int>& masklist)
     2868{
     2869  if (masklist.size() % 2 != 0) {
     2870    throw(AipsError("masklist must have even number of elements."));
     2871  }
     2872
     2873  std::vector<bool> res(nchan);
     2874
     2875  for (int i = 0; i < nchan; ++i) {
     2876    res[i] = false;
     2877  }
     2878  for (uInt j = 0; j < masklist.size(); j += 2) {
     2879    for (int i = masklist[j]; i <= masklist[j+1]; ++i) {
     2880      res[i] = true;
     2881    }
     2882  }
     2883
     2884  return res;
     2885}
     2886
     2887void Scantable::polyBaseline(const std::vector<bool>& mask, int order, float thresClip, int nIterClip, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile, const std::string& bltable)
     2888{
     2889  /****
     2890  double TimeStart = mathutil::gettimeofday_sec();
     2891  ****/
     2892
    24942893  try {
    24952894    ofstream ofs;
     
    25112910    }
    25122911
    2513     Fitter fitter = Fitter();
    2514     fitter.setExpression("poly", order);
    2515     //fitter.setIterClipping(thresClip, nIterClip);
    2516 
    2517     int nRow = nrow();
    2518     std::vector<bool> chanMask;
    25192912    bool showProgress;
    25202913    int minNRow;
    25212914    parseProgressInfo(progressInfo, showProgress, minNRow);
    25222915
     2916    int nRow = nrow();
     2917    std::vector<bool> chanMask;
     2918    std::vector<bool> finalChanMask;
     2919    float rms;
     2920
     2921    //---
     2922    bool outBaselintTable = (bltable != "");
     2923    STBaselineTable* bt = new STBaselineTable(*this);
     2924    ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(table_, "TIME");
     2925    Vector<Double> timeSecCol = tcol->getColumn();
     2926    //---
     2927
    25232928    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
     2929      std::vector<float> sp = getSpectrum(whichrow);
    25242930      chanMask = getCompositeChanMask(whichrow, mask);
    2525       fitBaseline(chanMask, whichrow, fitter);
    2526       setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
    2527       outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "polyBaseline()", fitter);
     2931      std::vector<float> params(order+1);
     2932      int nClipped = 0;
     2933      std::vector<float> res = doPolynomialFitting(sp, chanMask, order, params, rms, finalChanMask, nClipped, thresClip, nIterClip, getResidual);
     2934
     2935      //---
     2936      if (outBaselintTable) {
     2937        Vector<Int> fparam(1);
     2938        fparam[0] = order;
     2939        std::vector<int> finalMaskList0;
     2940        finalMaskList0.clear();
     2941        for (uInt i = 0; i < finalChanMask.size(); ++i) {
     2942          if (finalChanMask[i]) {
     2943            if ((i == 0)||(i == finalChanMask.size()-1)) {
     2944              finalMaskList0.push_back(i);
     2945            } else {
     2946              if ((finalChanMask[i])&&(!finalChanMask[i-1])) {
     2947                finalMaskList0.push_back(i);
     2948              }
     2949              if ((finalChanMask[i])&&(!finalChanMask[i+1])) {
     2950                finalMaskList0.push_back(i);
     2951              }
     2952            }
     2953          }
     2954        }
     2955        Vector<uInt> finalMaskList(finalMaskList0.size());
     2956        for (uInt i = 0; i < finalMaskList0.size(); ++i) {
     2957          finalMaskList[i] = (uInt)finalMaskList0[i];
     2958        }
     2959
     2960        bt->appenddata(uInt(getScan(whichrow)), uInt(getCycle(whichrow)), uInt(getBeam(whichrow)), uInt(getIF(whichrow)), uInt(getPol(whichrow)),
     2961                       uInt(0),
     2962                       timeSecCol[whichrow],
     2963                       Bool(true),
     2964                       STBaselineFunc::Polynomial,
     2965                       fparam,
     2966                       Vector<Float>(),
     2967                       finalMaskList,
     2968                       Vector<Float>(params),
     2969                       Float(rms),
     2970                       uInt(sp.size()),
     2971                       Float(thresClip),
     2972                       uInt(nIterClip),
     2973                       Float(0.0),
     2974                       uInt(0),
     2975                       Vector<uInt>());
     2976      } else {
     2977        setSpectrum(res, whichrow);
     2978      }
     2979      //---
     2980
     2981      outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "chebyshevBaseline()", params, nClipped);
    25282982      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
    25292983    }
     2984   
     2985    //---
     2986    if (outBaselintTable) {
     2987      bt->save(bltable);
     2988    }
     2989    //---
     2990    delete bt;
    25302991
    25312992    if (outTextFile) ofs.close();
     
    25342995    throw;
    25352996  }
    2536 }
    2537 
    2538 void Scantable::autoPolyBaseline(const std::vector<bool>& mask, int order, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
     2997
     2998  /****
     2999  double TimeEnd = mathutil::gettimeofday_sec();
     3000  double elapse1 = TimeEnd - TimeStart;
     3001  std::cout << "poly-new   : " << elapse1 << " (sec.)" << endl;
     3002  ****/
     3003}
     3004
     3005void Scantable::autoPolyBaseline(const std::vector<bool>& mask, int order, float thresClip, int nIterClip, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile, const std::string& bltable)
    25393006{
    25403007  try {
     
    25573024    }
    25583025
    2559     Fitter fitter = Fitter();
    2560     fitter.setExpression("poly", order);
    2561     //fitter.setIterClipping(thresClip, nIterClip);
    2562 
    2563     int nRow = nrow();
    2564     std::vector<bool> chanMask;
     3026    bool showProgress;
     3027    int minNRow;
     3028    parseProgressInfo(progressInfo, showProgress, minNRow);
     3029
    25653030    int minEdgeSize = getIFNos().size()*2;
    25663031    STLineFinder lineFinder = STLineFinder();
    25673032    lineFinder.setOptions(threshold, 3, chanAvgLimit);
    25683033
    2569     bool showProgress;
    2570     int minNRow;
    2571     parseProgressInfo(progressInfo, showProgress, minNRow);
     3034    int nRow = nrow();
     3035    std::vector<bool> chanMask;
     3036    std::vector<bool> finalChanMask;
     3037    float rms;
     3038
     3039    //---
     3040    bool outBaselintTable = (bltable != "");
     3041    STBaselineTable* bt = new STBaselineTable(*this);
     3042    ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(table_, "TIME");
     3043    Vector<Double> timeSecCol = tcol->getColumn();
     3044    //---
    25723045
    25733046    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    2574 
    2575       //-------------------------------------------------------
    2576       //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
    2577       //-------------------------------------------------------
     3047      std::vector<float> sp = getSpectrum(whichrow);
     3048      //-- composote given channel mask and flagtra --------
    25783049      int edgeSize = edge.size();
    25793050      std::vector<int> currentEdge;
     3051      currentEdge.clear();
    25803052      if (edgeSize >= 2) {
    25813053        int idx = 0;
     
    25913063        throw(AipsError("Wrong length of edge element"));
    25923064      }
    2593       lineFinder.setData(getSpectrum(whichrow));
     3065      lineFinder.setData(sp);
    25943066      lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
    25953067      chanMask = lineFinder.getMask();
    2596       //-------------------------------------------------------
    2597 
    2598       fitBaseline(chanMask, whichrow, fitter);
    2599       setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
    2600 
    2601       outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoPolyBaseline()", fitter);
     3068      //-- composote given channel mask and flagtra END --------
     3069
     3070      std::vector<float> params(order+1);
     3071      int nClipped = 0;
     3072      std::vector<float> res = doPolynomialFitting(sp, chanMask, order, params, rms, finalChanMask, nClipped, thresClip, nIterClip, getResidual);
     3073
     3074      if (outBaselintTable) {
     3075        Vector<Int> fparam(1);
     3076        fparam[0] = order;
     3077        std::vector<int> finalMaskList0;
     3078        finalMaskList0.clear();
     3079        for (uInt i = 0; i < finalChanMask.size(); ++i) {
     3080          if (finalChanMask[i]) {
     3081            if ((i == 0)||(i == finalChanMask.size()-1)) {
     3082              finalMaskList0.push_back(i);
     3083            } else {
     3084              if ((finalChanMask[i])&&(!finalChanMask[i-1])) {
     3085                finalMaskList0.push_back(i);
     3086              }
     3087              if ((finalChanMask[i])&&(!finalChanMask[i+1])) {
     3088                finalMaskList0.push_back(i);
     3089              }
     3090            }
     3091          }
     3092        }
     3093        Vector<uInt> finalMaskList(finalMaskList0.size());
     3094        for (uInt i = 0; i < finalMaskList0.size(); ++i) {
     3095          finalMaskList[i] = (uInt)finalMaskList0[i];
     3096        }
     3097
     3098        Vector<uInt> cEdge(currentEdge.size());
     3099        for (uInt i = 0; i < currentEdge.size(); ++i) {
     3100          cEdge[i] = currentEdge[i];
     3101        }
     3102
     3103        bt->appenddata(uInt(getScan(whichrow)), uInt(getCycle(whichrow)), uInt(getBeam(whichrow)), uInt(getIF(whichrow)), uInt(getPol(whichrow)),
     3104                       uInt(0),
     3105                       timeSecCol[whichrow],
     3106                       Bool(true),
     3107                       STBaselineFunc::Polynomial,
     3108                       fparam,
     3109                       Vector<Float>(),
     3110                       finalMaskList,
     3111                       Vector<Float>(params),
     3112                       Float(rms),
     3113                       uInt(sp.size()),
     3114                       Float(thresClip),
     3115                       uInt(nIterClip),
     3116                       Float(threshold),
     3117                       uInt(chanAvgLimit),
     3118                       cEdge);
     3119      } else {
     3120        setSpectrum(res, whichrow);
     3121      }
     3122
     3123      outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoChebyshevBaseline()", params, nClipped);
    26023124      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
    26033125    }
     3126
     3127    //---
     3128    if (outBaselintTable) {
     3129      bt->save(bltable);
     3130    }
     3131    //---
     3132    delete bt;
    26043133
    26053134    if (outTextFile) ofs.close();
     
    26103139}
    26113140
    2612 void Scantable::chebyshevBaseline(const std::vector<bool>& mask, int order, float thresClip, int nIterClip, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
    2613 {
     3141void Scantable::chebyshevBaseline(const std::vector<bool>& mask, int order, float thresClip, int nIterClip, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile, const std::string& bltable)
     3142{
     3143  //
     3144  double TimeStart = mathutil::gettimeofday_sec();
     3145  //
     3146
    26143147  try {
    26153148    ofstream ofs;
     
    26413174
    26423175    //---
    2643     bool outBaselineParamTable = false;
     3176    bool outBaselintTable = (bltable != "");
    26443177    STBaselineTable* bt = new STBaselineTable(*this);
     3178    ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(table_, "TIME");
     3179    Vector<Double> timeSecCol = tcol->getColumn();
    26453180    //---
    26463181
     
    26503185      std::vector<float> params(order+1);
    26513186      int nClipped = 0;
    2652       std::vector<float> res = doChebyshevFitting(sp, chanMask, order, params, finalChanMask, rms, nClipped, thresClip, nIterClip, getResidual);
     3187      std::vector<float> res = doChebyshevFitting(sp, chanMask, order, params, rms, finalChanMask, nClipped, thresClip, nIterClip, getResidual);
    26533188
    26543189      //---
    2655       if (outBaselineParamTable) {
     3190      if (outBaselintTable) {
     3191        Vector<Int> fparam(1);
     3192        fparam[0] = order;
     3193        std::vector<int> finalMaskList0;
     3194        finalMaskList0.clear();
     3195        for (uInt i = 0; i < finalChanMask.size(); ++i) {
     3196          if (finalChanMask[i]) {
     3197            if ((i == 0)||(i == finalChanMask.size()-1)) {
     3198              finalMaskList0.push_back(i);
     3199            } else {
     3200              if ((finalChanMask[i])&&(!finalChanMask[i-1])) {
     3201                finalMaskList0.push_back(i);
     3202              }
     3203              if ((finalChanMask[i])&&(!finalChanMask[i+1])) {
     3204                finalMaskList0.push_back(i);
     3205              }
     3206            }
     3207          }
     3208        }
     3209        Vector<uInt> finalMaskList(finalMaskList0.size());
     3210        for (uInt i = 0; i < finalMaskList0.size(); ++i) {
     3211          finalMaskList[i] = (uInt)finalMaskList0[i];
     3212        }
     3213
    26563214        bt->appenddata(uInt(getScan(whichrow)), uInt(getCycle(whichrow)), uInt(getBeam(whichrow)), uInt(getIF(whichrow)), uInt(getPol(whichrow)),
    26573215                       uInt(0),
    2658                        Double(0.0), // <-- Double(getTime(whichrow, false)),
    2659                        uInt(nchan(whichrow)),
     3216                       timeSecCol[whichrow],
     3217                       Bool(true),
    26603218                       STBaselineFunc::Chebyshev,
    2661                        Vector<uInt>(1), // ==> MUST BE Vector<uInt> containing 'order'.
    2662                        Vector<Float>(), // ==> for ffpar. ** dummy **
     3219                       fparam,
     3220                       Vector<Float>(),
     3221                       finalMaskList,
     3222                       Vector<Float>(params),
     3223                       Float(rms),
     3224                       uInt(sp.size()),
     3225                       Float(thresClip),
    26633226                       uInt(nIterClip),
    2664                        Float(thresClip),
    2665                        Vector<uInt>(5), // <-- Vector<uInt>(finalChanMask),
    2666                        Vector<Float>(params),
    2667                        Float(rms));
     3227                       Float(0.0),
     3228                       uInt(0),
     3229                       Vector<uInt>());
    26683230      } else {
    26693231        setSpectrum(res, whichrow);
     
    26763238   
    26773239    //---
    2678     if (outBaselineParamTable) {
    2679       bt->save("chebyparamtable");
     3240    if (outBaselintTable) {
     3241      bt->save(bltable);
    26803242    }
    26813243    //---
     
    26873249    throw;
    26883250  }
     3251
     3252  double TimeEnd = mathutil::gettimeofday_sec();
     3253  double elapse1 = TimeEnd - TimeStart;
     3254  std::cout << "cheby   : " << elapse1 << " (sec.)" << endl;
     3255}
     3256
     3257void Scantable::autoChebyshevBaseline(const std::vector<bool>& mask, int order, float thresClip, int nIterClip, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile, const std::string& bltable)
     3258{
     3259  try {
     3260    ofstream ofs;
     3261    String coordInfo = "";
     3262    bool hasSameNchan = true;
     3263    bool outTextFile = false;
     3264    bool csvFormat = false;
     3265
     3266    if (blfile != "") {
     3267      csvFormat = (blfile.substr(0, 1) == "T");
     3268      ofs.open(blfile.substr(1).c_str(), ios::out | ios::app);
     3269      if (ofs) outTextFile = true;
     3270    }
     3271
     3272    if (outLogger || outTextFile) {
     3273      coordInfo = getCoordInfo()[0];
     3274      if (coordInfo == "") coordInfo = "channel";
     3275      hasSameNchan = hasSameNchanOverIFs();
     3276    }
     3277
     3278    bool showProgress;
     3279    int minNRow;
     3280    parseProgressInfo(progressInfo, showProgress, minNRow);
     3281
     3282    int minEdgeSize = getIFNos().size()*2;
     3283    STLineFinder lineFinder = STLineFinder();
     3284    lineFinder.setOptions(threshold, 3, chanAvgLimit);
     3285
     3286    int nRow = nrow();
     3287    std::vector<bool> chanMask;
     3288    std::vector<bool> finalChanMask;
     3289    float rms;
     3290
     3291    //---
     3292    bool outBaselintTable = (bltable != "");
     3293    STBaselineTable* bt = new STBaselineTable(*this);
     3294    ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(table_, "TIME");
     3295    Vector<Double> timeSecCol = tcol->getColumn();
     3296    //---
     3297
     3298    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
     3299      std::vector<float> sp = getSpectrum(whichrow);
     3300      //-- composote given channel mask and flagtra --------
     3301      int edgeSize = edge.size();
     3302      std::vector<int> currentEdge;
     3303      currentEdge.clear();
     3304      if (edgeSize >= 2) {
     3305        int idx = 0;
     3306        if (edgeSize > 2) {
     3307          if (edgeSize < minEdgeSize) {
     3308            throw(AipsError("Length of edge element info is less than that of IFs"));
     3309          }
     3310          idx = 2 * getIF(whichrow);
     3311        }
     3312        currentEdge.push_back(edge[idx]);
     3313        currentEdge.push_back(edge[idx+1]);
     3314      } else {
     3315        throw(AipsError("Wrong length of edge element"));
     3316      }
     3317      lineFinder.setData(sp);
     3318      lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
     3319      chanMask = lineFinder.getMask();
     3320      //-- composote given channel mask and flagtra END --------
     3321
     3322      std::vector<float> params(order+1);
     3323      int nClipped = 0;
     3324      std::vector<float> res = doChebyshevFitting(sp, chanMask, order, params, rms, finalChanMask, nClipped, thresClip, nIterClip, getResidual);
     3325
     3326      if (outBaselintTable) {
     3327        Vector<Int> fparam(1);
     3328        fparam[0] = order;
     3329        std::vector<int> finalMaskList0;
     3330        finalMaskList0.clear();
     3331        for (uInt i = 0; i < finalChanMask.size(); ++i) {
     3332          if (finalChanMask[i]) {
     3333            if ((i == 0)||(i == finalChanMask.size()-1)) {
     3334              finalMaskList0.push_back(i);
     3335            } else {
     3336              if ((finalChanMask[i])&&(!finalChanMask[i-1])) {
     3337                finalMaskList0.push_back(i);
     3338              }
     3339              if ((finalChanMask[i])&&(!finalChanMask[i+1])) {
     3340                finalMaskList0.push_back(i);
     3341              }
     3342            }
     3343          }
     3344        }
     3345        Vector<uInt> finalMaskList(finalMaskList0.size());
     3346        for (uInt i = 0; i < finalMaskList0.size(); ++i) {
     3347          finalMaskList[i] = (uInt)finalMaskList0[i];
     3348        }
     3349
     3350        Vector<uInt> cEdge(currentEdge.size());
     3351        for (uInt i = 0; i < currentEdge.size(); ++i) {
     3352          cEdge[i] = currentEdge[i];
     3353        }
     3354
     3355        bt->appenddata(uInt(getScan(whichrow)), uInt(getCycle(whichrow)), uInt(getBeam(whichrow)), uInt(getIF(whichrow)), uInt(getPol(whichrow)),
     3356                       uInt(0),
     3357                       timeSecCol[whichrow],
     3358                       Bool(true),
     3359                       STBaselineFunc::Chebyshev,
     3360                       fparam,
     3361                       Vector<Float>(),
     3362                       finalMaskList,
     3363                       Vector<Float>(params),
     3364                       Float(rms),
     3365                       uInt(sp.size()),
     3366                       Float(thresClip),
     3367                       uInt(nIterClip),
     3368                       Float(threshold),
     3369                       uInt(chanAvgLimit),
     3370                       cEdge);
     3371      } else {
     3372        setSpectrum(res, whichrow);
     3373      }
     3374
     3375      outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoChebyshevBaseline()", params, nClipped);
     3376      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
     3377    }
     3378
     3379    //---
     3380    if (outBaselintTable) {
     3381      bt->save(bltable);
     3382    }
     3383    //---
     3384    delete bt;
     3385
     3386    if (outTextFile) ofs.close();
     3387
     3388  } catch (...) {
     3389    throw;
     3390  }
    26893391}
    26903392
    26913393double Scantable::calculateModelSelectionCriteria(const std::string& valname, const std::string& blfunc, int order, const std::vector<bool>& inMask, int whichrow, bool useLineFinder, const std::vector<int>& edge, float threshold, int chanAvgLimit)
    26923394{
     3395  std::vector<float> sp = getSpectrum(whichrow);
     3396  std::vector<bool> chanMask;
     3397  chanMask.clear();
     3398
    26933399  if (useLineFinder) {
    26943400    int minEdgeSize = getIFNos().size()*2;
    26953401
    2696 
    26973402    int edgeSize = edge.size();
    26983403    std::vector<int> currentEdge;
     3404    currentEdge.clear();
    26993405    if (edgeSize >= 2) {
    27003406      int idx = 0;
     
    27123418
    27133419    STLineFinder lineFinder = STLineFinder();
    2714     std::vector<float> sp = getSpectrum(whichrow);
    27153420    lineFinder.setData(sp);
    27163421    lineFinder.setOptions(threshold, 3, chanAvgLimit);
    27173422    lineFinder.findLines(getCompositeChanMask(whichrow, inMask), currentEdge, whichrow);
    2718     std::vector<bool> chanMask = lineFinder.getMask();
    2719    
    2720     return doCalculateModelSelectionCriteria(valname, sp, chanMask, blfunc, order);
     3423    chanMask = lineFinder.getMask();
    27213424
    27223425  } else {
    2723     return doCalculateModelSelectionCriteria(valname, getSpectrum(whichrow), getCompositeChanMask(whichrow, inMask), blfunc, order);
    2724   }
    2725 
     3426
     3427    chanMask = getCompositeChanMask(whichrow, inMask);
     3428  }
     3429
     3430  return doCalculateModelSelectionCriteria(valname, sp, chanMask, blfunc, order);
    27263431}
    27273432
     
    27333438  float rms;
    27343439  int nClipped = 0;
    2735   float thresClip = 0.0;
    2736   int nIterClip = 0;
    27373440  std::vector<float> res;
    2738   if (blfunc == "chebyshev") {
     3441  if (blfunc == "poly") {
    27393442    nparam = order + 1;
    2740     res = doChebyshevFitting(spec, mask, order, params, finalChanMask, rms, nClipped, thresClip, nIterClip, true);
     3443    res = doPolynomialFitting(spec, mask, order, params, rms, finalChanMask, nClipped);
     3444  } else if (blfunc == "chebyshev") {
     3445    nparam = order + 1;
     3446    res = doChebyshevFitting(spec, mask, order, params, rms, finalChanMask, nClipped);
     3447  } else if (blfunc == "cspline") {
     3448    std::vector<int> pieceEdges;//(order+1);  //order = npiece
     3449    nparam = order + 3;
     3450    //params.resize(4*order);
     3451    res = doCubicSplineFitting(spec, mask, order, false, pieceEdges, params, rms, finalChanMask, nClipped);
    27413452  } else if (blfunc == "sinusoid") {
    27423453    std::vector<int> nWaves;
     
    27463457    }
    27473458    nparam = 2*order + 1;  // order = nwave
    2748     res = doSinusoidFitting(spec, mask, nWaves, params, nClipped, thresClip, nIterClip, true);
    2749   } else if (blfunc == "cspline") {
    2750     std::vector<int> pieceEdges(order+1);  //order = npiece
    2751     nparam = order + 3;
    2752     params.resize(4*order);
    2753     res = doCubicSplineFitting(spec, mask, order, pieceEdges, params, nClipped, thresClip, nIterClip, true);
     3459    res = doSinusoidFitting(spec, mask, nWaves, params, rms, finalChanMask, nClipped);
    27543460  } else {
    2755     throw(AipsError("blfunc must be chebyshev, cspline or sinusoid."));
     3461    throw(AipsError("blfunc must be poly, chebyshev, cspline or sinusoid."));
    27563462  }
    27573463
     
    27773483    double aic = nusedchan * (log(2.0 * PI * msq) + 1.0) + 2.0 * nparam;
    27783484
    2779     // Corrected AIC by Sugiura(1978)
     3485    // Corrected AIC by Sugiura(1978) (AICc)
    27803486    if (valname == "aicc") {
    27813487      if (nusedchan - nparam - 1 <= 0) {
     
    28033509}
    28043510
    2805 void Scantable::autoChebyshevBaseline(const std::vector<bool>& mask, int order, float thresClip, int nIterClip, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
    2806 {
    2807   try {
    2808     ofstream ofs;
    2809     String coordInfo = "";
    2810     bool hasSameNchan = true;
    2811     bool outTextFile = false;
    2812     bool csvFormat = false;
    2813 
    2814     if (blfile != "") {
    2815       csvFormat = (blfile.substr(0, 1) == "T");
    2816       ofs.open(blfile.substr(1).c_str(), ios::out | ios::app);
    2817       if (ofs) outTextFile = true;
    2818     }
    2819 
    2820     if (outLogger || outTextFile) {
    2821       coordInfo = getCoordInfo()[0];
    2822       if (coordInfo == "") coordInfo = "channel";
    2823       hasSameNchan = hasSameNchanOverIFs();
    2824     }
    2825 
    2826     int nRow = nrow();
    2827     std::vector<bool> chanMask;
    2828     int minEdgeSize = getIFNos().size()*2;
    2829     STLineFinder lineFinder = STLineFinder();
    2830     lineFinder.setOptions(threshold, 3, chanAvgLimit);
    2831 
    2832     bool showProgress;
    2833     int minNRow;
    2834     parseProgressInfo(progressInfo, showProgress, minNRow);
    2835 
    2836     std::vector<bool> finalChanMask;
    2837     float rms;
    2838 
    2839     for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    2840       std::vector<float> sp = getSpectrum(whichrow);
    2841 
    2842       //-------------------------------------------------------
    2843       //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
    2844       //-------------------------------------------------------
    2845       int edgeSize = edge.size();
    2846       std::vector<int> currentEdge;
    2847       if (edgeSize >= 2) {
    2848         int idx = 0;
    2849         if (edgeSize > 2) {
    2850           if (edgeSize < minEdgeSize) {
    2851             throw(AipsError("Length of edge element info is less than that of IFs"));
    2852           }
    2853           idx = 2 * getIF(whichrow);
    2854         }
    2855         currentEdge.push_back(edge[idx]);
    2856         currentEdge.push_back(edge[idx+1]);
    2857       } else {
    2858         throw(AipsError("Wrong length of edge element"));
    2859       }
    2860       //lineFinder.setData(getSpectrum(whichrow));
    2861       lineFinder.setData(sp);
    2862       lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
    2863       chanMask = lineFinder.getMask();
    2864       //-------------------------------------------------------
    2865 
    2866 
    2867       //fitBaseline(chanMask, whichrow, fitter);
    2868       //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
    2869       std::vector<float> params(order+1);
    2870       int nClipped = 0;
    2871       std::vector<float> res = doChebyshevFitting(sp, chanMask, order, params, finalChanMask, rms, nClipped, thresClip, nIterClip, getResidual);
    2872       setSpectrum(res, whichrow);
    2873 
    2874       outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoChebyshevBaseline()", params, nClipped);
    2875       showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
    2876     }
    2877 
    2878     if (outTextFile) ofs.close();
    2879 
    2880   } catch (...) {
    2881     throw;
    2882   }
    2883 }
    2884 
    2885   /*
    2886 double Scantable::getChebyshevPolynomial(int n, double x) {
    2887   if ((x < -1.0)||(x > 1.0)) {
    2888     throw(AipsError("out of definition range (-1 <= x <= 1)."));
    2889   } else if (n < 0) {
    2890     throw(AipsError("the order must be zero or positive."));
    2891   } else if (n == 0) {
     3511double Scantable::getNormalPolynomial(int n, double x) {
     3512  if (n == 0) {
    28923513    return 1.0;
    2893   } else if (n == 1) {
    2894     return x;
     3514  } else if (n > 0) {
     3515    double res = 1.0;
     3516    for (int i = 0; i < n; ++i) {
     3517      res *= x;
     3518    }
     3519    return res;
    28953520  } else {
    2896     return 2.0*x*getChebyshevPolynomial(n-1, x) - getChebyshevPolynomial(n-2, x);
    2897   }
    2898 }
    2899   */
     3521    if (x == 0.0) {
     3522      throw(AipsError("infinity result: x=0 given for negative power."));
     3523    } else {
     3524      return pow(x, (double)n);
     3525    }
     3526  }
     3527}
     3528
    29003529double Scantable::getChebyshevPolynomial(int n, double x) {
    29013530  if ((x < -1.0)||(x > 1.0)) {
     
    29393568}
    29403569
    2941 std::vector<float> Scantable::doChebyshevFitting(const std::vector<float>& data, const std::vector<bool>& mask, int order, std::vector<float>& params, std::vector<bool>& finalMask, float& rms, int& nClipped, float thresClip, int nIterClip, bool getResidual)
     3570std::vector<float> Scantable::doPolynomialFitting(const std::vector<float>& data, const std::vector<bool>& mask, int order, std::vector<float>& params, float& rms, std::vector<bool>& finalmask, float clipth, int clipn)
     3571{
     3572  int nClipped = 0;
     3573  return doPolynomialFitting(data, mask, order, params, rms, finalmask, nClipped, clipth, clipn);
     3574}
     3575
     3576std::vector<float> Scantable::doPolynomialFitting(const std::vector<float>& data, const std::vector<bool>& mask, int order, std::vector<float>& params, float& rms, std::vector<bool>& finalMask, int& nClipped, float thresClip, int nIterClip, bool getResidual)
     3577{
     3578  return doLeastSquareFitting(data, mask, &Scantable::getNormalPolynomial, order, params, rms, finalMask, nClipped, thresClip, nIterClip, getResidual);
     3579}
     3580
     3581std::vector<float> Scantable::doChebyshevFitting(const std::vector<float>& data, const std::vector<bool>& mask, int order, std::vector<float>& params, float& rms, std::vector<bool>& finalmask, float clipth, int clipn)
     3582{
     3583  int nClipped = 0;
     3584  return doChebyshevFitting(data, mask, order, params, rms, finalmask, nClipped, clipth, clipn);
     3585}
     3586
     3587std::vector<float> Scantable::doChebyshevFitting(const std::vector<float>& data, const std::vector<bool>& mask, int order, std::vector<float>& params, float& rms, std::vector<bool>& finalMask, int& nClipped, float thresClip, int nIterClip, bool getResidual)
     3588{
     3589  return doLeastSquareFitting(data, mask, &Scantable::getChebyshevPolynomial, order, params, rms, finalMask, nClipped, thresClip, nIterClip, getResidual);
     3590}
     3591
     3592std::vector<float> Scantable::doLeastSquareFitting(const std::vector<float>& data, const std::vector<bool>& mask, double (Scantable::*pfunc)(int, double), int order, std::vector<float>& params, float& rms, std::vector<bool>& finalMask, int& nClipped, float thresClip, int nIterClip, bool getResidual)
    29423593{
    29433594  if (data.size() != mask.size()) {
     
    29453596  }
    29463597  if (order < 0) {
    2947     throw(AipsError("maximum order of Chebyshev polynomial must not be negative."));
     3598    throw(AipsError("maximum order of polynomial must not be negative."));
    29483599  }
    29493600
     
    29713622  //          xArray.size() is nDOF and xArray[*].size() is nChan.
    29723623  //          Each xArray element are as follows:
     3624  //
     3625  //          (for normal polynomials)
     3626  //          xArray[0]   = {1.0,   1.0,   1.0,   ..., 1.0},
     3627  //          xArray[1]   = {0.0,   1.0,   2.0,   ..., (nChan-1)}
     3628  //          xArray[n-1] = ...,
     3629  //          xArray[n]   = {0.0^n, 1.0^n, 2.0^n, ..., (nChan-1)^n}
     3630  //          where (0 <= n <= order)
     3631  //
     3632  //          (for Chebyshev polynomials)
    29733633  //          xArray[0]   = {T0(-1), T0(2/(nChan-1)-1), T0(4/(nChan-1)-1), ..., T0(1)},
    29743634  //          xArray[n-1] = ...,
     
    29763636  //          where (0 <= n <= order),
    29773637  std::vector<std::vector<double> > xArray;
     3638  double xFactor, xShift;
     3639  if (pfunc == &Scantable::getChebyshevPolynomial) {
     3640    xFactor = 2.0/(double)(nChan - 1);
     3641    xShift  = -1.0;
     3642  } else {
     3643    xFactor = 1.0;
     3644    xShift  = 0.0;
     3645  }
    29783646  for (int i = 0; i < nDOF; ++i) {
    2979     double xFactor = 2.0/(double)(nChan - 1);
    29803647    std::vector<double> xs;
    29813648    xs.clear();
    29823649    for (int j = 0; j < nChan; ++j) {
    2983       xs.push_back(getChebyshevPolynomial(i, xFactor*(double)j-1.0));
     3650      xs.push_back((this->*pfunc)(i, xFactor*(double)j + xShift));
    29843651    }
    29853652    xArray.push_back(xs);
     
    30633730      }
    30643731    }
    3065     //compute a vector y which consists of the coefficients of the sinusoids forming the
    3066     //best-fit curves (a0,s1,c1,s2,c2,...), where a0 is constant and s* and c* are of sine
    3067     //and cosine functions, respectively.
     3732    //compute a vector y which consists of the coefficients of Chebyshev polynomials.
    30683733    std::vector<double> y;
    30693734    params.clear();
     
    31183783
    31193784  std::vector<float> result;
     3785  result.clear();
    31203786  if (getResidual) {
    31213787    for (int i = 0; i < nChan; ++i) {
     
    31313797}
    31323798
    3133 void Scantable::cubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
    3134 {
    3135   /*************
    3136   double totTimeStart, totTimeEnd, blTimeStart, blTimeEnd, ioTimeStart, ioTimeEnd, msTimeStart, msTimeEnd, seTimeStart, seTimeEnd, otTimeStart, otTimeEnd, prTimeStart, prTimeEnd;
    3137   double elapseMs = 0.0;
    3138   double elapseSe = 0.0;
    3139   double elapseOt = 0.0;
    3140   double elapsePr = 0.0;
    3141   double elapseBl = 0.0;
    3142   double elapseIo = 0.0;
    3143   totTimeStart = mathutil::gettimeofday_sec();
    3144   *************/
    3145 
     3799void Scantable::cubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile, const std::string& bltable)
     3800{
    31463801  try {
    31473802    ofstream ofs;
     
    31633818    }
    31643819
    3165     //Fitter fitter = Fitter();
    3166     //fitter.setExpression("cspline", nPiece);
    3167     //fitter.setIterClipping(thresClip, nIterClip);
    3168 
    31693820    bool showProgress;
    31703821    int minNRow;
     
    31733824    int nRow = nrow();
    31743825    std::vector<bool> chanMask;
    3175 
    3176     //--------------------------------
     3826    std::vector<bool> finalChanMask;
     3827    float rms;
     3828
     3829    //---
     3830    bool outBaselintTable = (bltable != "");
     3831    STBaselineTable* bt = new STBaselineTable(*this);
     3832    ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(table_, "TIME");
     3833    Vector<Double> timeSecCol = tcol->getColumn();
     3834    //---
     3835
     3836
    31773837    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    3178       /******************
    3179       ioTimeStart = mathutil::gettimeofday_sec();
    3180       **/
    31813838      std::vector<float> sp = getSpectrum(whichrow);
    3182       /**
    3183       ioTimeEnd = mathutil::gettimeofday_sec();
    3184       elapseIo += (double)(ioTimeEnd - ioTimeStart);
    3185       msTimeStart = mathutil::gettimeofday_sec();
    3186       ******************/
    3187 
    31883839      chanMask = getCompositeChanMask(whichrow, mask);
    3189 
    3190       /**
    3191       msTimeEnd = mathutil::gettimeofday_sec();
    3192       elapseMs += (double)(msTimeEnd - msTimeStart);
    3193       blTimeStart = mathutil::gettimeofday_sec();
    3194       **/
    3195 
    3196       //fitBaseline(chanMask, whichrow, fitter);
    3197       //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
    3198       std::vector<int> pieceEdges(nPiece+1);
    3199       std::vector<float> params(nPiece*4);
     3840      std::vector<int> pieceEdges;//(nPiece+1);
     3841      std::vector<float> params;//(nPiece*4);
    32003842      int nClipped = 0;
    3201       std::vector<float> res = doCubicSplineFitting(sp, chanMask, nPiece, pieceEdges, params, nClipped, thresClip, nIterClip, getResidual);
    3202 
    3203       /**
    3204       blTimeEnd = mathutil::gettimeofday_sec();
    3205       elapseBl += (double)(blTimeEnd - blTimeStart);
    3206       seTimeStart = mathutil::gettimeofday_sec();
    3207       **/
    3208 
    3209 
    3210       setSpectrum(res, whichrow);
    3211       //
    3212 
    3213       /**
    3214       seTimeEnd = mathutil::gettimeofday_sec();
    3215       elapseSe += (double)(seTimeEnd - seTimeStart);
    3216       otTimeStart = mathutil::gettimeofday_sec();
    3217       **/
     3843      std::vector<float> res = doCubicSplineFitting(sp, chanMask, nPiece, false, pieceEdges, params, rms, finalChanMask, nClipped, thresClip, nIterClip, getResidual);
     3844
     3845      //---
     3846      if (outBaselintTable) {
     3847        Vector<Int> fparam(nPiece);
     3848        for (uInt i = 0; i < fparam.size(); ++i) {
     3849          fparam[i] = pieceEdges[i];
     3850        }
     3851        std::vector<int> finalMaskList0;
     3852        finalMaskList0.clear();
     3853        for (uInt i = 0; i < finalChanMask.size(); ++i) {
     3854          if (finalChanMask[i]) {
     3855            if ((i == 0)||(i == finalChanMask.size()-1)) {
     3856              finalMaskList0.push_back(i);
     3857            } else {
     3858              if ((finalChanMask[i])&&(!finalChanMask[i-1])) {
     3859                finalMaskList0.push_back(i);
     3860              }
     3861              if ((finalChanMask[i])&&(!finalChanMask[i+1])) {
     3862                finalMaskList0.push_back(i);
     3863              }
     3864            }
     3865          }
     3866        }
     3867        Vector<uInt> finalMaskList(finalMaskList0.size());
     3868        for (uInt i = 0; i < finalMaskList0.size(); ++i) {
     3869          finalMaskList[i] = (uInt)finalMaskList0[i];
     3870        }
     3871
     3872        bt->appenddata(uInt(getScan(whichrow)), uInt(getCycle(whichrow)), uInt(getBeam(whichrow)), uInt(getIF(whichrow)), uInt(getPol(whichrow)),
     3873                       uInt(0),
     3874                       timeSecCol[whichrow],
     3875                       Bool(true),
     3876                       STBaselineFunc::CSpline,
     3877                       fparam,
     3878                       Vector<Float>(),
     3879                       finalMaskList,
     3880                       Vector<Float>(params),
     3881                       Float(rms),
     3882                       uInt(sp.size()),
     3883                       Float(thresClip),
     3884                       uInt(nIterClip),
     3885                       Float(0.0),
     3886                       uInt(0),
     3887                       Vector<uInt>());
     3888      } else {
     3889        setSpectrum(res, whichrow);
     3890      }
     3891      //---
    32183892
    32193893      outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()", pieceEdges, params, nClipped);
    3220 
    3221       /**
    3222       otTimeEnd = mathutil::gettimeofday_sec();
    3223       elapseOt += (double)(otTimeEnd - otTimeStart);
    3224       prTimeStart = mathutil::gettimeofday_sec();
    3225       **/
    3226 
    32273894      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
    3228 
    3229       /******************
    3230       prTimeEnd = mathutil::gettimeofday_sec();
    3231       elapsePr += (double)(prTimeEnd - prTimeStart);
    3232       ******************/
    3233     }
    3234     //--------------------------------
     3895    }
    32353896   
     3897    //---
     3898    if (outBaselintTable) {
     3899      bt->save(bltable);
     3900    }
     3901    //---
     3902    delete bt;
     3903
    32363904    if (outTextFile) ofs.close();
    32373905
     
    32393907    throw;
    32403908  }
    3241   /***************
    3242   totTimeEnd = mathutil::gettimeofday_sec();
    3243   std::cout << "io    : " << elapseIo << " (sec.)" << endl;
    3244   std::cout << "ms    : " << elapseMs << " (sec.)" << endl;
    3245   std::cout << "bl    : " << elapseBl << " (sec.)" << endl;
    3246   std::cout << "se    : " << elapseSe << " (sec.)" << endl;
    3247   std::cout << "ot    : " << elapseOt << " (sec.)" << endl;
    3248   std::cout << "pr    : " << elapsePr << " (sec.)" << endl;
    3249   std::cout << "total : " << (double)(totTimeEnd - totTimeStart) << " (sec.)" << endl;
    3250   ***************/
    3251 }
    3252 
    3253 void Scantable::autoCubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
     3909}
     3910
     3911void Scantable::autoCubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile, const std::string& bltable)
    32543912{
    32553913  try {
     
    32723930    }
    32733931
    3274     //Fitter fitter = Fitter();
    3275     //fitter.setExpression("cspline", nPiece);
    3276     //fitter.setIterClipping(thresClip, nIterClip);
    3277 
    3278     int nRow = nrow();
    3279     std::vector<bool> chanMask;
     3932    bool showProgress;
     3933    int minNRow;
     3934    parseProgressInfo(progressInfo, showProgress, minNRow);
     3935
    32803936    int minEdgeSize = getIFNos().size()*2;
    32813937    STLineFinder lineFinder = STLineFinder();
    32823938    lineFinder.setOptions(threshold, 3, chanAvgLimit);
    32833939
    3284     bool showProgress;
    3285     int minNRow;
    3286     parseProgressInfo(progressInfo, showProgress, minNRow);
     3940    int nRow = nrow();
     3941    std::vector<bool> chanMask;
     3942    std::vector<bool> finalChanMask;
     3943    float rms;
     3944
     3945    //---
     3946    bool outBaselintTable = (bltable != "");
     3947    STBaselineTable* bt = new STBaselineTable(*this);
     3948    ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(table_, "TIME");
     3949    Vector<Double> timeSecCol = tcol->getColumn();
     3950    //---
    32873951
    32883952    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    32893953      std::vector<float> sp = getSpectrum(whichrow);
    3290 
    3291       //-------------------------------------------------------
    3292       //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
    3293       //-------------------------------------------------------
     3954      //-- composote given channel mask and flagtra --------
    32943955      int edgeSize = edge.size();
    32953956      std::vector<int> currentEdge;
     3957      currentEdge.clear();
    32963958      if (edgeSize >= 2) {
    32973959        int idx = 0;
     
    33073969        throw(AipsError("Wrong length of edge element"));
    33083970      }
    3309       //lineFinder.setData(getSpectrum(whichrow));
    33103971      lineFinder.setData(sp);
    33113972      lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
    33123973      chanMask = lineFinder.getMask();
    3313       //-------------------------------------------------------
    3314 
    3315 
    3316       //fitBaseline(chanMask, whichrow, fitter);
    3317       //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
    3318       std::vector<int> pieceEdges(nPiece+1);
    3319       std::vector<float> params(nPiece*4);
     3974      //-- composote given channel mask and flagtra END --------
     3975
     3976      std::vector<int> pieceEdges;//(nPiece+1);
     3977      std::vector<float> params;//(nPiece*4);
    33203978      int nClipped = 0;
    3321       std::vector<float> res = doCubicSplineFitting(sp, chanMask, nPiece, pieceEdges, params, nClipped, thresClip, nIterClip, getResidual);
    3322       setSpectrum(res, whichrow);
    3323       //
     3979      std::vector<float> res = doCubicSplineFitting(sp, chanMask, nPiece, false, pieceEdges, params, rms, finalChanMask, nClipped, thresClip, nIterClip, getResidual);
     3980
     3981      if (outBaselintTable) {
     3982        Vector<Int> fparam(nPiece);
     3983        for (uInt i = 0; i < fparam.size(); ++i) {
     3984          fparam[i] = pieceEdges[i];
     3985        }
     3986        std::vector<int> finalMaskList0;
     3987        finalMaskList0.clear();
     3988        for (uInt i = 0; i < finalChanMask.size(); ++i) {
     3989          if (finalChanMask[i]) {
     3990            if ((i == 0)||(i == finalChanMask.size()-1)) {
     3991              finalMaskList0.push_back(i);
     3992            } else {
     3993              if ((finalChanMask[i])&&(!finalChanMask[i-1])) {
     3994                finalMaskList0.push_back(i);
     3995              }
     3996              if ((finalChanMask[i])&&(!finalChanMask[i+1])) {
     3997                finalMaskList0.push_back(i);
     3998              }
     3999            }
     4000          }
     4001        }
     4002        Vector<uInt> finalMaskList(finalMaskList0.size());
     4003        for (uInt i = 0; i < finalMaskList0.size(); ++i) {
     4004          finalMaskList[i] = (uInt)finalMaskList0[i];
     4005        }
     4006
     4007        Vector<uInt> cEdge(currentEdge.size());
     4008        for (uInt i = 0; i < currentEdge.size(); ++i) {
     4009          cEdge[i] = currentEdge[i];
     4010        }
     4011
     4012        bt->appenddata(uInt(getScan(whichrow)), uInt(getCycle(whichrow)), uInt(getBeam(whichrow)), uInt(getIF(whichrow)), uInt(getPol(whichrow)),
     4013                       uInt(0),
     4014                       timeSecCol[whichrow],
     4015                       Bool(true),
     4016                       STBaselineFunc::CSpline,
     4017                       fparam,
     4018                       Vector<Float>(),
     4019                       finalMaskList,
     4020                       Vector<Float>(params),
     4021                       Float(rms),
     4022                       uInt(sp.size()),
     4023                       Float(thresClip),
     4024                       uInt(nIterClip),
     4025                       Float(threshold),
     4026                       uInt(chanAvgLimit),
     4027                       cEdge);
     4028      } else {
     4029        setSpectrum(res, whichrow);
     4030      }
    33244031
    33254032      outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()", pieceEdges, params, nClipped);
     
    33274034    }
    33284035
     4036    //---
     4037    if (outBaselintTable) {
     4038      bt->save(bltable);
     4039    }
     4040    //---
     4041    delete bt;
     4042
    33294043    if (outTextFile) ofs.close();
    33304044
     
    33344048}
    33354049
    3336 std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data, const std::vector<bool>& mask, int nPiece, std::vector<int>& idxEdge, std::vector<float>& params, int& nClipped, float thresClip, int nIterClip, bool getResidual)
     4050std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data, const std::vector<bool>& mask, std::vector<int>& idxEdge, std::vector<float>& params, float& rms, std::vector<bool>& finalmask, float clipth, int clipn)
     4051{
     4052  int nClipped = 0;
     4053  return doCubicSplineFitting(data, mask, idxEdge.size()-1, true, idxEdge, params, rms, finalmask, nClipped, clipth, clipn);
     4054}
     4055
     4056std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data, const std::vector<bool>& mask, int nPiece, std::vector<int>& idxEdge, std::vector<float>& params, float& rms, std::vector<bool>& finalmask, float clipth, int clipn)
     4057{
     4058  int nClipped = 0;
     4059  return doCubicSplineFitting(data, mask, nPiece, false, idxEdge, params, rms, finalmask, nClipped, clipth, clipn);
     4060}
     4061
     4062std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data, const std::vector<bool>& mask, int nPiece, bool useGivenPieceBoundary, std::vector<int>& idxEdge, std::vector<float>& params, float& rms, std::vector<bool>& finalMask, int& nClipped, float thresClip, int nIterClip, bool getResidual)
    33374063{
    33384064  if (data.size() != mask.size()) {
     
    33444070
    33454071  int nChan = data.size();
     4072
     4073  finalMask.clear();
     4074  finalMask.resize(nChan);
     4075
    33464076  std::vector<int> maskArray(nChan);
    33474077  std::vector<int> x(nChan);
     
    33534083      j++;
    33544084    }
     4085    finalMask[i] = mask[i];
    33554086  }
    33564087  int initNData = j;
     
    33624093  int nElement = (int)(floor(floor((double)(initNData/nPiece))+0.5));
    33634094  std::vector<double> invEdge(nPiece-1);
    3364   idxEdge[0] = x[0];
     4095
     4096  if (useGivenPieceBoundary) {
     4097    if ((int)idxEdge.size() != nPiece+1) {
     4098      throw(AipsError("pieceEdge.size() must be equal to nPiece+1."));
     4099    }
     4100  } else {
     4101    idxEdge.clear();
     4102    idxEdge.resize(nPiece+1);
     4103    idxEdge[0] = x[0];
     4104  }
    33654105  for (int i = 1; i < nPiece; ++i) {
    33664106    int valX = x[nElement*i];
    3367     idxEdge[i] = valX;
     4107    if (!useGivenPieceBoundary) {
     4108      idxEdge[i] = valX;
     4109    }
    33684110    invEdge[i-1] = 1.0/(double)valX;
    33694111  }
    3370   idxEdge[nPiece] = x[initNData-1]+1;
     4112  if (!useGivenPieceBoundary) {
     4113    idxEdge[nPiece] = x[initNData-1]+1;
     4114  }
     4115
     4116  params.clear();
     4117  params.resize(nPiece*4);
    33714118
    33724119  int nData = initNData;
     
    33954142    //           identity matrix (right).
    33964143    // the right part is used to calculate the inverse matrix of the left part.
     4144
    33974145    double xMatrix[nDOF][2*nDOF];
    33984146    double zMatrix[nDOF];
     
    34964244      }
    34974245    }
     4246
    34984247    //compute a vector y which consists of the coefficients of the best-fit spline curves
    34994248    //(a0,a1,a2,a3(,b3,c3,...)), namely, the ones for the leftmost piece and the ones of
     
    35454294      }
    35464295    }
     4296
    35474297    if (idxEdge[nPiece] < nChan) {
    35484298      int n = idxEdge[nPiece]-1;
     
    35624312    }
    35634313
     4314    double stdDev = 0.0;
     4315    for (int i = 0; i < nChan; ++i) {
     4316      stdDev += residual[i]*residual[i]*(double)maskArray[i];
     4317    }
     4318    stdDev = sqrt(stdDev/(double)nData);
     4319    rms = (float)stdDev;
     4320
    35644321    if ((nClip == nIterClip) || (thresClip <= 0.0)) {
    35654322      break;
    35664323    } else {
    3567       double stdDev = 0.0;
    3568       for (int i = 0; i < nChan; ++i) {
    3569         stdDev += residual[i]*residual[i]*(double)maskArray[i];
    3570       }
    3571       stdDev = sqrt(stdDev/(double)nData);
    35724324     
    35734325      double thres = stdDev * thresClip;
     
    35764328        if (abs(residual[i]) >= thres) {
    35774329          maskArray[i] = 0;
     4330          finalMask[i] = false;
    35784331        }
    35794332        if (maskArray[i] > 0) {
     
    35864339        nData = newNData;
    35874340      }
     4341
    35884342    }
    35894343  }
     
    36054359}
    36064360
    3607 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)
    3608 {
     4361void Scantable::selectWaveNumbers(const int whichrow, const std::vector<bool>& chanMask, const std::string& fftInfo, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, std::vector<int>& nWaves)
     4362//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)
     4363{
     4364  bool applyFFT;
     4365  std::string fftMethod;
     4366  std::string fftThresh;
     4367
    36094368  nWaves.clear();
    36104369
     4370  parseFFTInfo(fftInfo, applyFFT, fftMethod, fftThresh);
    36114371  if (applyFFT) {
    36124372    string fftThAttr;
     
    36184378
    36194379  addAuxWaveNumbers(whichrow, addNWaves, rejectNWaves, nWaves);
     4380}
     4381
     4382void Scantable::parseFFTInfo(const std::string& fftInfo, bool& applyFFT, std::string& fftMethod, std::string& fftThresh)
     4383{
     4384  istringstream iss(fftInfo);
     4385  std::string tmp;
     4386  std::vector<string> res;
     4387  while (getline(iss, tmp, ',')) {
     4388    res.push_back(tmp);
     4389  }
     4390  if (res.size() < 3) {
     4391    throw(AipsError("wrong value in 'fftinfo' parameter")) ;
     4392  }
     4393  applyFFT = (res[0] == "true");
     4394  fftMethod = res[1];
     4395  fftThresh = res[2];
    36204396}
    36214397
     
    37024478{
    37034479  std::vector<int> tempAddNWaves, tempRejectNWaves;
     4480  tempAddNWaves.clear();
     4481  tempRejectNWaves.clear();
     4482
    37044483  for (uInt i = 0; i < addNWaves.size(); ++i) {
    37054484    tempAddNWaves.push_back(addNWaves[i]);
     
    37454524void Scantable::setWaveNumberListUptoNyquistFreq(const int whichrow, std::vector<int>& nWaves)
    37464525{
    3747   if ((nWaves.size() == 2)&&(nWaves[1] == -999)) {
    3748     int val = nWaves[0];
    3749     int nyquistFreq = nchan(getIF(whichrow))/2+1;
    3750     nWaves.clear();
    3751     if (val > nyquistFreq) {  // for safety, at least nWaves contains a constant; CAS-3759
    3752       nWaves.push_back(0);
    3753     }
    3754     while (val <= nyquistFreq) {
    3755       nWaves.push_back(val);
    3756       val++;
    3757     }
    3758   }
    3759 }
    3760 
    3761 void 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, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
     4526  int val = nWaves[0];
     4527  int nyquistFreq = nchan(getIF(whichrow))/2+1;
     4528  nWaves.clear();
     4529  if (val > nyquistFreq) {  // for safety, at least nWaves contains a constant; CAS-3759
     4530    nWaves.push_back(0);
     4531  }
     4532  while (val <= nyquistFreq) {
     4533    nWaves.push_back(val);
     4534    val++;
     4535  }
     4536}
     4537
     4538void Scantable::sinusoidBaseline(const std::vector<bool>& mask, const std::string& fftInfo, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, float thresClip, int nIterClip, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile, const std::string& bltable)
    37624539{
    37634540  try {
     
    37804557    }
    37814558
    3782     //Fitter fitter = Fitter();
    3783     //fitter.setExpression("sinusoid", nWaves);
    3784     //fitter.setIterClipping(thresClip, nIterClip);
     4559    bool showProgress;
     4560    int minNRow;
     4561    parseProgressInfo(progressInfo, showProgress, minNRow);
    37854562
    37864563    int nRow = nrow();
    37874564    std::vector<bool> chanMask;
    37884565    std::vector<int> nWaves;
    3789 
    3790     bool showProgress;
    3791     int minNRow;
    3792     parseProgressInfo(progressInfo, showProgress, minNRow);
     4566    std::vector<bool> finalChanMask;
     4567    float rms;
     4568
     4569    //---
     4570    bool outBaselintTable = (bltable != "");
     4571    STBaselineTable* bt = new STBaselineTable(*this);
     4572    ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(table_, "TIME");
     4573    Vector<Double> timeSecCol = tcol->getColumn();
     4574    //---
    37934575
    37944576    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
     4577      std::vector<float> sp = getSpectrum(whichrow);
    37954578      chanMask = getCompositeChanMask(whichrow, mask);
    3796       selectWaveNumbers(whichrow, chanMask, applyFFT, fftMethod, fftThresh, addNWaves, rejectNWaves, nWaves);
    3797 
    3798       //FOR DEBUGGING------------
    3799       /*
    3800       if (whichrow < 0) {// == nRow -1) {
    3801         cout << "+++ i=" << setw(3) << whichrow << ", IF=" << setw(2) << getIF(whichrow);
    3802         if (applyFFT) {
    3803           cout << "[ ";
    3804           for (uInt j = 0; j < nWaves.size(); ++j) {
    3805             cout << nWaves[j] << ", ";
    3806           }
    3807           cout << " ]    " << endl;
    3808         }
    3809         cout << flush;
    3810       }
    3811       */
    3812       //-------------------------
    3813 
    3814       //fitBaseline(chanMask, whichrow, fitter);
    3815       //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
     4579      selectWaveNumbers(whichrow, chanMask, fftInfo, addNWaves, rejectNWaves, nWaves);
     4580
    38164581      std::vector<float> params;
    38174582      int nClipped = 0;
    3818       std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, params, nClipped, thresClip, nIterClip, getResidual);
    3819       setSpectrum(res, whichrow);
    3820       //
     4583      std::vector<float> res = doSinusoidFitting(sp, chanMask, nWaves, params, rms, finalChanMask, nClipped, thresClip, nIterClip, getResidual);
     4584
     4585      //---
     4586      if (outBaselintTable) {
     4587        uInt nparam = nWaves.size();
     4588        Vector<Int> fparam(nparam);
     4589        for (uInt i = 0; i < nparam; ++i) {
     4590          fparam[i] = nWaves[i];
     4591        }
     4592        std::vector<int> finalMaskList0;
     4593        finalMaskList0.clear();
     4594        for (uInt i = 0; i < finalChanMask.size(); ++i) {
     4595          if (finalChanMask[i]) {
     4596            if ((i == 0)||(i == finalChanMask.size()-1)) {
     4597              finalMaskList0.push_back(i);
     4598            } else {
     4599              if ((finalChanMask[i])&&(!finalChanMask[i-1])) {
     4600                finalMaskList0.push_back(i);
     4601              }
     4602              if ((finalChanMask[i])&&(!finalChanMask[i+1])) {
     4603                finalMaskList0.push_back(i);
     4604              }
     4605            }
     4606          }
     4607        }
     4608        Vector<uInt> finalMaskList(finalMaskList0.size());
     4609        for (uInt i = 0; i < finalMaskList0.size(); ++i) {
     4610          finalMaskList[i] = (uInt)finalMaskList0[i];
     4611        }
     4612
     4613        bt->appenddata(uInt(getScan(whichrow)), uInt(getCycle(whichrow)), uInt(getBeam(whichrow)), uInt(getIF(whichrow)), uInt(getPol(whichrow)),
     4614                       uInt(0),
     4615                       timeSecCol[whichrow],
     4616                       Bool(true),
     4617                       STBaselineFunc::Sinusoid,
     4618                       fparam,
     4619                       Vector<Float>(),
     4620                       finalMaskList,
     4621                       Vector<Float>(params),
     4622                       Float(rms),
     4623                       uInt(sp.size()),
     4624                       Float(thresClip),
     4625                       uInt(nIterClip),
     4626                       Float(0.0),
     4627                       uInt(0),
     4628                       Vector<uInt>());
     4629      } else {
     4630        setSpectrum(res, whichrow);
     4631      }
     4632      //---
    38214633
    38224634      outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "sinusoidBaseline()", params, nClipped);
     
    38244636    }
    38254637
     4638    //---
     4639    if (outBaselintTable) {
     4640      bt->save(bltable);
     4641    }
     4642    //---
     4643    delete bt;
     4644
    38264645    if (outTextFile) ofs.close();
    38274646
     
    38314650}
    38324651
    3833 void 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, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
     4652void Scantable::autoSinusoidBaseline(const std::vector<bool>& mask, const std::string& fftInfo, 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, const std::string& progressInfo, const bool outLogger, const std::string& blfile, const std::string& bltable)
     4653//void 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, const std::string& progressInfo, const bool outLogger, const std::string& blfile, const std::string& bltable)
    38344654{
    38354655  try {
     
    38524672    }
    38534673
    3854     //Fitter fitter = Fitter();
    3855     //fitter.setExpression("sinusoid", nWaves);
    3856     //fitter.setIterClipping(thresClip, nIterClip);
     4674    bool showProgress;
     4675    int minNRow;
     4676    parseProgressInfo(progressInfo, showProgress, minNRow);
     4677
     4678    int minEdgeSize = getIFNos().size()*2;
     4679    STLineFinder lineFinder = STLineFinder();
     4680    lineFinder.setOptions(threshold, 3, chanAvgLimit);
    38574681
    38584682    int nRow = nrow();
    38594683    std::vector<bool> chanMask;
    38604684    std::vector<int> nWaves;
    3861 
    3862     int minEdgeSize = getIFNos().size()*2;
    3863     STLineFinder lineFinder = STLineFinder();
    3864     lineFinder.setOptions(threshold, 3, chanAvgLimit);
    3865 
    3866     bool showProgress;
    3867     int minNRow;
    3868     parseProgressInfo(progressInfo, showProgress, minNRow);
     4685    std::vector<bool> finalChanMask;
     4686    float rms;
     4687
     4688    //---
     4689    bool outBaselintTable = (bltable != "");
     4690    STBaselineTable* bt = new STBaselineTable(*this);
     4691    ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(table_, "TIME");
     4692    Vector<Double> timeSecCol = tcol->getColumn();
     4693    //---
    38694694
    38704695    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    3871 
    3872       //-------------------------------------------------------
    3873       //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
    3874       //-------------------------------------------------------
     4696      std::vector<float> sp = getSpectrum(whichrow);
     4697      //-- composote given channel mask and flagtra --------
    38754698      int edgeSize = edge.size();
    38764699      std::vector<int> currentEdge;
     4700      currentEdge.clear();
    38774701      if (edgeSize >= 2) {
    38784702        int idx = 0;
     
    38884712        throw(AipsError("Wrong length of edge element"));
    38894713      }
    3890       lineFinder.setData(getSpectrum(whichrow));
     4714      lineFinder.setData(sp);
    38914715      lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
    38924716      chanMask = lineFinder.getMask();
    3893       //-------------------------------------------------------
    3894 
    3895       selectWaveNumbers(whichrow, chanMask, applyFFT, fftMethod, fftThresh, addNWaves, rejectNWaves, nWaves);
    3896 
    3897       //fitBaseline(chanMask, whichrow, fitter);
    3898       //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
     4717      //-- composote given channel mask and flagtra END --------
     4718      selectWaveNumbers(whichrow, chanMask, fftInfo, addNWaves, rejectNWaves, nWaves);
     4719
    38994720      std::vector<float> params;
    39004721      int nClipped = 0;
    3901       std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, params, nClipped, thresClip, nIterClip, getResidual);
    3902       setSpectrum(res, whichrow);
    3903       //
     4722      std::vector<float> res = doSinusoidFitting(sp, chanMask, nWaves, params, rms, finalChanMask, nClipped, thresClip, nIterClip, getResidual);
     4723
     4724      //---
     4725      if (outBaselintTable) {
     4726        uInt nparam = nWaves.size();
     4727        Vector<Int> fparam(nparam);
     4728        for (uInt i = 0; i < nparam; ++i) {
     4729          fparam[i] = nWaves[i];
     4730        }
     4731        std::vector<int> finalMaskList0;
     4732        finalMaskList0.clear();
     4733        for (uInt i = 0; i < finalChanMask.size(); ++i) {
     4734          if (finalChanMask[i]) {
     4735            if ((i == 0)||(i == finalChanMask.size()-1)) {
     4736              finalMaskList0.push_back(i);
     4737            } else {
     4738              if ((finalChanMask[i])&&(!finalChanMask[i-1])) {
     4739                finalMaskList0.push_back(i);
     4740              }
     4741              if ((finalChanMask[i])&&(!finalChanMask[i+1])) {
     4742                finalMaskList0.push_back(i);
     4743              }
     4744            }
     4745          }
     4746        }
     4747        Vector<uInt> finalMaskList(finalMaskList0.size());
     4748        for (uInt i = 0; i < finalMaskList0.size(); ++i) {
     4749          finalMaskList[i] = (uInt)finalMaskList0[i];
     4750        }
     4751
     4752        Vector<uInt> cEdge(currentEdge.size());
     4753        for (uInt i = 0; i < currentEdge.size(); ++i) {
     4754          cEdge[i] = currentEdge[i];
     4755        }
     4756
     4757        bt->appenddata(uInt(getScan(whichrow)), uInt(getCycle(whichrow)), uInt(getBeam(whichrow)), uInt(getIF(whichrow)), uInt(getPol(whichrow)),
     4758                       uInt(0),
     4759                       timeSecCol[whichrow],
     4760                       Bool(true),
     4761                       STBaselineFunc::Sinusoid,
     4762                       fparam,
     4763                       Vector<Float>(),
     4764                       finalMaskList,
     4765                       Vector<Float>(params),
     4766                       Float(rms),
     4767                       uInt(sp.size()),
     4768                       Float(thresClip),
     4769                       uInt(nIterClip),
     4770                       Float(threshold),
     4771                       uInt(chanAvgLimit),
     4772                       cEdge);
     4773      } else {
     4774        setSpectrum(res, whichrow);
     4775      }
     4776      //---
    39044777
    39054778      outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoSinusoidBaseline()", params, nClipped);
     
    39074780    }
    39084781
     4782    //---
     4783    if (outBaselintTable) {
     4784      bt->save(bltable);
     4785    }
     4786    //---
     4787    delete bt;
     4788
    39094789    if (outTextFile) ofs.close();
    39104790
     
    39144794}
    39154795
    3916 std::vector<float> Scantable::doSinusoidFitting(const std::vector<float>& data, const std::vector<bool>& mask, const std::vector<int>& waveNumbers, std::vector<float>& params, int& nClipped, float thresClip, int nIterClip, bool getResidual)
     4796std::vector<float> Scantable::doSinusoidFitting(const std::vector<float>& data, const std::vector<bool>& mask, const std::vector<int>& waveNumbers, std::vector<float>& params, float& rms, std::vector<bool>& finalmask, float clipth, int clipn)
     4797{
     4798  int nClipped = 0;
     4799  return doSinusoidFitting(data, mask, waveNumbers, params, rms, finalmask, nClipped, clipth, clipn);
     4800}
     4801
     4802std::vector<float> Scantable::doSinusoidFitting(const std::vector<float>& data, const std::vector<bool>& mask, const std::vector<int>& waveNumbers, std::vector<float>& params, float& rms, std::vector<bool>& finalMask, int& nClipped, float thresClip, int nIterClip, bool getResidual)
    39174803{
    39184804  if (data.size() != mask.size()) {
     
    39394825
    39404826  int nChan = data.size();
     4827
     4828  finalMask.clear();
     4829  finalMask.resize(nChan);
     4830
    39414831  std::vector<int> maskArray;
    39424832  std::vector<int> x;
     
    39464836      x.push_back(i);
    39474837    }
     4838    finalMask[i] = mask[i];
    39484839  }
    39494840
     
    40864977    }
    40874978
     4979    double stdDev = 0.0;
     4980    for (int i = 0; i < nChan; ++i) {
     4981      stdDev += residual[i]*residual[i]*(double)maskArray[i];
     4982    }
     4983    stdDev = sqrt(stdDev/(double)nData);
     4984    rms = (float)stdDev;
     4985
    40884986    if ((nClip == nIterClip) || (thresClip <= 0.0)) {
    40894987      break;
    40904988    } else {
    4091       double stdDev = 0.0;
    4092       for (int i = 0; i < nChan; ++i) {
    4093         stdDev += residual[i]*residual[i]*(double)maskArray[i];
    4094       }
    4095       stdDev = sqrt(stdDev/(double)nData);
    4096      
     4989
    40974990      double thres = stdDev * thresClip;
    40984991      int newNData = 0;
     
    41004993        if (abs(residual[i]) >= thres) {
    41014994          maskArray[i] = 0;
     4995          finalMask[i] = false;
    41024996        }
    41034997        if (maskArray[i] > 0) {
     
    41105004        nData = newNData;
    41115005      }
     5006
    41125007    }
    41135008  }
  • trunk/src/Scantable.h

    r2737 r2767  
    2828#include <casa/Quanta/Quantum.h>
    2929#include <casa/Utilities/CountedPtr.h>
    30 
    3130#include <coordinates/Coordinates/SpectralCoordinate.h>
    32 
    3331#include <measures/TableMeasures/ScalarMeasColumn.h>
    34 
    3532#include <scimath/Mathematics/FFTServer.h>
    36 
    3733#include <tables/Tables/ArrayColumn.h>
    3834#include <tables/Tables/ScalarColumn.h>
    3935#include <tables/Tables/Table.h>
    4036
    41 
    4237#include "MathUtils.h"
     38#include "STBaselineEnum.h"
    4339#include "STFit.h"
    4440#include "STFitEntry.h"
     
    507503  void regridChannel( int nchan, double dnu, int irow ) ;
    508504  void regridChannel( int nchan, double dnu, double fmin, int irow ) ;
    509 
    510505  void regridSpecChannel( double dnu, int nchan=-1 ) ;
    511506
    512507  bool getFlagtraFast(casa::uInt whichrow);
    513508
    514   void polyBaseline(const std::vector<bool>& mask,
    515                     int order,
     509  std::vector<std::string> applyBaselineTable(const std::string& bltable, const std::string& outbltable, const bool outbltableexists, const bool overwrite);
     510  std::vector<std::string> subBaseline(const std::vector<std::string>& blInfoList, const std::string& outbltable, const bool outbltableexists, const bool overwrite);
     511  void polyBaseline(const std::vector<bool>& mask,
     512                    int order,
     513                    float thresClip,
     514                    int nIterClip,
    516515                    bool getResidual=true,
    517516                    const std::string& progressInfo="true,1000",
    518                     const bool outLogger=false,
    519                     const std::string& blfile="");
     517                    const bool outLogger=false,
     518                    const std::string& blfile="",
     519                    const std::string& bltable="");
    520520  void autoPolyBaseline(const std::vector<bool>& mask,
    521                         int order,
    522                         const std::vector<int>& edge,
    523                         float threshold=3.0,
    524                         int chanAvgLimit=1,
    525                         bool getResidual=true,
    526                         const std::string& progressInfo="true,1000",
    527                         const bool outLogger=false,
    528                         const std::string& blfile="");
     521                        int order,
     522                        float thresClip,
     523                        int nIterClip,
     524                        const std::vector<int>& edge,
     525                        float threshold=3.0,
     526                        int chanAvgLimit=1,
     527                        bool getResidual=true,
     528                        const std::string& progressInfo="true,1000",
     529                        const bool outLogger=false,
     530                        const std::string& blfile="",
     531                        const std::string& bltable="");
    529532  void chebyshevBaseline(const std::vector<bool>& mask,
    530533                         int order,
     
    534537                         const std::string& progressInfo="true,1000",
    535538                         const bool outLogger=false,
    536                          const std::string& blfile="");
     539                         const std::string& blfile="",
     540                         const std::string& bltable="");
    537541  void autoChebyshevBaseline(const std::vector<bool>& mask,
    538542                             int order,
     
    545549                             const std::string& progressInfo="true,1000",
    546550                             const bool outLogger=false,
    547                              const std::string& blfile="");
     551                             const std::string& blfile="",
     552                             const std::string& bltable="");
    548553  void cubicSplineBaseline(const std::vector<bool>& mask,
    549554                           int nPiece,
     
    553558                           const std::string& progressInfo="true,1000",
    554559                           const bool outLogger=false,
    555                            const std::string& blfile="");
     560                           const std::string& blfile="",
     561                           const std::string& bltable="");
    556562  void autoCubicSplineBaseline(const std::vector<bool>& mask,
    557563                               int nPiece,
     
    564570                               const std::string& progressInfo="true,1000",
    565571                               const bool outLogger=false,
    566                                const std::string& blfile="");
     572                               const std::string& blfile="",
     573                               const std::string& bltable="");
    567574  void sinusoidBaseline(const std::vector<bool>& mask,
    568                         const bool applyFFT,
    569                         const std::string& fftMethod,
    570                         const std::string& fftThresh,
     575                        const std::string& fftInfo,
    571576                        const std::vector<int>& addNWaves,
    572577                        const std::vector<int>& rejectNWaves,
     
    576581                        const std::string& progressInfo="true,1000",
    577582                        const bool outLogger=false,
    578                         const std::string& blfile="");
     583                        const std::string& blfile="",
     584                        const std::string& bltable="");
    579585  void autoSinusoidBaseline(const std::vector<bool>& mask,
    580                             const bool applyFFT,
    581                             const std::string& fftMethod,
    582                             const std::string& fftThresh,
     586                            const std::string& fftInfo,
    583587                            const std::vector<int>& addNWaves,
    584588                            const std::vector<int>& rejectNWaves,
     
    591595                            const std::string& progressInfo="true,1000",
    592596                            const bool outLogger=false,
    593                             const std::string& blfile="");
     597                            const std::string& blfile="",
     598                            const std::string& bltable="");
    594599  std::vector<float> execFFT(const int whichrow,
    595600                             const std::vector<bool>& inMask,
     
    628633                                         float threshold,
    629634                                         int chanAvgLimit);
     635  static std::vector<bool> getMaskFromMaskList(const int nchan,
     636                                               const std::vector<int>& masklist);
     637  static std::vector<int> splitToIntList(const std::string& str, const char delim);
     638  static std::vector<string> splitToStringList(const std::string& str, const char delim);
     639
    630640
    631641
     
    743753
    744754  void fitBaseline(const std::vector<bool>& mask, int whichrow, Fitter& fitter);
     755  double getNormalPolynomial(int n, double x);
    745756  double getChebyshevPolynomial(int n, double x);
     757  std::vector<float> doPolynomialFitting(const std::vector<float>& data,
     758                                         const std::vector<bool>& mask,
     759                                         int order,
     760                                         std::vector<float>& params,
     761                                         float& rms,
     762                                         std::vector<bool>& finalMask,
     763                                         float clipth,
     764                                         int clipn);
     765  std::vector<float> doPolynomialFitting(const std::vector<float>& data,
     766                                         const std::vector<bool>& mask,
     767                                         int order,
     768                                         std::vector<float>& params,
     769                                         float& rms,
     770                                         std::vector<bool>& finalMask,
     771                                         int& nClipped,
     772                                         float thresClip=3.0,
     773                                         int nIterClip=0,
     774                                         bool getResidual=true);
    746775  std::vector<float> doChebyshevFitting(const std::vector<float>& data,
     776                                        const std::vector<bool>& mask,
     777                                        int order,
     778                                        std::vector<float>& params,
     779                                        float& rms,
     780                                        std::vector<bool>& finalMask,
     781                                        float clipth,
     782                                        int clipn);
     783  std::vector<float> doChebyshevFitting(const std::vector<float>& data,
     784                                        const std::vector<bool>& mask,
     785                                        int order,
     786                                        std::vector<float>& params,
     787                                        float& rms,
     788                                        std::vector<bool>& finalMask,
     789                                        int& nClipped,
     790                                        float thresClip=3.0,
     791                                        int nIterClip=0,
     792                                        bool getResidual=true);
     793  std::vector<float> doLeastSquareFitting(const std::vector<float>& data,
     794                                          const std::vector<bool>& mask,
     795                                          double (Scantable::*pfunc)(int, double),
     796                                          int order,
     797                                          std::vector<float>& params,
     798                                          float& rms,
     799                                          std::vector<bool>& finalMask,
     800                                          int& nClipped,
     801                                          float thresClip=3.0,
     802                                          int nIterClip=0,
     803                                          bool getResidual=true);
     804  std::vector<float> doCubicSplineFitting(const std::vector<float>& data,
    747805                                          const std::vector<bool>& mask,
    748                                           int order,
     806                                          std::vector<int>& idxEdge,
    749807                                          std::vector<float>& params,
     808                                          float& rms,
    750809                                          std::vector<bool>& finalMask,
    751                                           float& rms,
    752                                           int& nClipped,
    753                                           float thresClip=3.0,
    754                                           int nIterClip=1,
    755                                           bool getResidual=true);
     810                                          float clipth,
     811                                          int clipn);
    756812  std::vector<float> doCubicSplineFitting(const std::vector<float>& data,
    757813                                          const std::vector<bool>& mask,
     
    759815                                          std::vector<int>& idxEdge,
    760816                                          std::vector<float>& params,
     817                                          float& rms,
     818                                          std::vector<bool>& finalMask,
     819                                          float clipth,
     820                                          int clipn);
     821  std::vector<float> doCubicSplineFitting(const std::vector<float>& data,
     822                                          const std::vector<bool>& mask,
     823                                          int nPiece,
     824                                          bool useGivenPieceBoundary,
     825                                          std::vector<int>& idxEdge,
     826                                          std::vector<float>& params,
     827                                          float& rms,
     828                                          std::vector<bool>& finalMask,
    761829                                          int& nClipped,
    762830                                          float thresClip=3.0,
    763                                           int nIterClip=1,
     831                                          int nIterClip=0,
    764832                                          bool getResidual=true);
    765833  std::vector<float> doSinusoidFitting(const std::vector<float>& data,
     
    767835                                       const std::vector<int>& waveNumbers,
    768836                                       std::vector<float>& params,
     837                                       float& rms,
     838                                       std::vector<bool>& finalMask,
     839                                       float clipth,
     840                                       int clipn);
     841  std::vector<float> doSinusoidFitting(const std::vector<float>& data,
     842                                       const std::vector<bool>& mask,
     843                                       const std::vector<int>& waveNumbers,
     844                                       std::vector<float>& params,
     845                                       float& rms,
     846                                       std::vector<bool>& finalMask,
    769847                                       int& nClipped,
    770848                                       float thresClip=3.0,
    771                                        int nIterClip=1,
     849                                       int nIterClip=0,
    772850                                       bool getResidual=true);
    773851  void selectWaveNumbers(const int whichrow,
    774852                         const std::vector<bool>& chanMask,
    775                          const bool applyFFT,
    776                          const std::string& fftMethod,
    777                          const std::string& fftThresh,
     853                         const std::string& fftInfo,
     854                         //const bool applyFFT,
     855                         //const std::string& fftMethod,
     856                         //const std::string& fftThresh,
    778857                         const std::vector<int>& addNWaves,
    779858                         const std::vector<int>& rejectNWaves,
    780859                         std::vector<int>& nWaves);
     860  void parseFFTInfo(const std::string& fftInfo,
     861                    bool& applyFFT,
     862                    std::string& fftMethod,
     863                    std::string& fftThresh);
    781864  void parseThresholdExpression(const std::string& fftThresh,
    782865                                std::string& fftThAttr,
     
    821904                                           int order);
    822905  double doGetRms(const std::vector<bool>& mask, const casa::Vector<casa::Float>& spec);
     906  std::string packFittingResults(const int irow, const std::vector<float>& params, const float rms);
     907  void parseBlInfo(const std::string& blInfo, int& whichrow, STBaselineFunc::FuncName& ftype, std::vector<int>& fpar, std::vector<bool>& mask, float& thresClip, int& nIterClip, bool& useLineFinder, float& thresLF, std::vector<int>& edgeLF, int& avgLF);
     908 std::vector<float> doApplyBaselineTable(std::vector<float>& spec, std::vector<bool>& mask, const STBaselineFunc::FuncName ftype, std::vector<int>& fpar, std::vector<float>& params, float&rms);
     909 std::vector<float> doSubtractBaseline(std::vector<float>& spec, std::vector<bool>& mask, const STBaselineFunc::FuncName ftype, std::vector<int>& fpar, std::vector<float>& params, float&rms, std::vector<bool>& finalmask, float clipth, int clipn, bool uself, int irow, float lfth, std::vector<int>& lfedge, int lfavg);
    823910
    824911};
    825 
    826912} // namespace
    827913
  • trunk/src/ScantableWrapper.h

    r2713 r2767  
    278278  { table_->regridSpecChannel( dnu, nchan ); }
    279279
    280   void polyBaseline(const std::vector<bool>& mask, int order, bool getresidual=true, const std::string& showprogress="true,1000", const bool outlog=false, const std::string& blfile="")
    281   { table_->polyBaseline(mask, order, getresidual, showprogress, outlog, blfile); }
    282 
    283   void autoPolyBaseline(const std::vector<bool>& mask, int order, const std::vector<int>& edge, float threshold=5.0, int chan_avg_limit=1, bool getresidual=true, const std::string& showprogress="true,1000", const bool outlog=false, const std::string& blfile="")
    284   { table_->autoPolyBaseline(mask, order, edge, threshold, chan_avg_limit, getresidual, showprogress, outlog, blfile); }
    285 
    286   void chebyshevBaseline(const std::vector<bool>& mask, int order, float clipthresh, int clipniter, bool getresidual=true, const std::string& showprogress="true,1000", const bool outlog=false, const std::string& blfile="")
    287   { table_->chebyshevBaseline(mask, order, clipthresh, clipniter, getresidual, showprogress, outlog, blfile); }
    288 
    289   void autoChebyshevBaseline(const std::vector<bool>& mask, int order, float clipthresh, int clipniter, const std::vector<int>& edge, float threshold=5.0, int chan_avg_limit=1, bool getresidual=true, const std::string& showprogress="true,1000", const bool outlog=false, const std::string& blfile="")
    290   { table_->autoChebyshevBaseline(mask, order, clipthresh, clipniter, edge, threshold, chan_avg_limit, getresidual, showprogress, outlog, blfile); }
    291 
    292   void cubicSplineBaseline(const std::vector<bool>& mask, int npiece, float clipthresh, int clipniter, bool getresidual=true, const std::string& showprogress="true,1000", const bool outlog=false, const std::string& blfile="")
    293   { table_->cubicSplineBaseline(mask, npiece, clipthresh, clipniter, getresidual, showprogress, outlog, blfile); }
    294 
    295   void autoCubicSplineBaseline(const std::vector<bool>& mask, int npiece, float clipthresh, int clipniter, const std::vector<int>& edge, float threshold=5.0, int chan_avg_limit=1, bool getresidual=true, const std::string& showprogress="true,1000", const bool outlog=false, const std::string& blfile="")
    296   { table_->autoCubicSplineBaseline(mask, npiece, clipthresh, clipniter, edge, threshold, chan_avg_limit, getresidual, showprogress, outlog, blfile); }
    297 
    298   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, const std::string& showprogress="true,1000", const bool outlog=false, const std::string& blfile="")
    299   { table_->sinusoidBaseline(mask, applyfft, fftmethod, fftthresh, addwn, rejwn, clipthresh, clipniter, getresidual, showprogress, outlog, blfile); }
    300 
    301   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, const std::string& showprogress="true,1000", const bool outlog=false, const std::string& blfile="")
    302   { table_->autoSinusoidBaseline(mask, applyfft, fftmethod, fftthresh, addwn, rejwn, clipthresh, clipniter, edge, threshold, chan_avg_limit, getresidual, showprogress, outlog, blfile); }
     280  std::vector<std::string> applyBaselineTable(const std::string& bltable, const std::string& outbltable, const bool outbltableexists, const bool overwrite)
     281  { return table_->applyBaselineTable(bltable, outbltable, outbltableexists, overwrite); }
     282  std::vector<std::string> subBaseline(const std::vector<std::string>& blinfo, const std::string& outbltable, const bool outbltableexists, const bool overwrite)
     283  { return table_->subBaseline(blinfo, outbltable, outbltableexists, overwrite); }
     284  void polyBaseline(const std::vector<bool>& mask, int order, float clipthresh, int clipniter, bool getresidual=true, const std::string& showprogress="true,1000", const bool outlog=false, const std::string& blfile="", const std::string& bltable="")
     285  { table_->polyBaseline(mask, order, clipthresh, clipniter, getresidual, showprogress, outlog, blfile, bltable); }
     286
     287  void autoPolyBaseline(const std::vector<bool>& mask, int order, float clipthresh, int clipniter, const std::vector<int>& edge, float threshold=5.0, int chan_avg_limit=1, bool getresidual=true, const std::string& showprogress="true,1000", const bool outlog=false, const std::string& blfile="", const std::string& bltable="")
     288  { table_->autoPolyBaseline(mask, order, clipthresh, clipniter, edge, threshold, chan_avg_limit, getresidual, showprogress, outlog, blfile, bltable); }
     289
     290  void chebyshevBaseline(const std::vector<bool>& mask, int order, float clipthresh, int clipniter, bool getresidual=true, const std::string& showprogress="true,1000", const bool outlog=false, const std::string& blfile="", const std::string& bltable="")
     291  { table_->chebyshevBaseline(mask, order, clipthresh, clipniter, getresidual, showprogress, outlog, blfile, bltable); }
     292
     293  void autoChebyshevBaseline(const std::vector<bool>& mask, int order, float clipthresh, int clipniter, const std::vector<int>& edge, float threshold=5.0, int chan_avg_limit=1, bool getresidual=true, const std::string& showprogress="true,1000", const bool outlog=false, const std::string& blfile="", const std::string& bltable="")
     294  { table_->autoChebyshevBaseline(mask, order, clipthresh, clipniter, edge, threshold, chan_avg_limit, getresidual, showprogress, outlog, blfile, bltable); }
     295
     296  void cubicSplineBaseline(const std::vector<bool>& mask, int npiece, float clipthresh, int clipniter, bool getresidual=true, const std::string& showprogress="true,1000", const bool outlog=false, const std::string& blfile="", const std::string& bltable="")
     297  { table_->cubicSplineBaseline(mask, npiece, clipthresh, clipniter, getresidual, showprogress, outlog, blfile, bltable); }
     298
     299  void autoCubicSplineBaseline(const std::vector<bool>& mask, int npiece, float clipthresh, int clipniter, const std::vector<int>& edge, float threshold=5.0, int chan_avg_limit=1, bool getresidual=true, const std::string& showprogress="true,1000", const bool outlog=false, const std::string& blfile="", const std::string& bltable="")
     300  { table_->autoCubicSplineBaseline(mask, npiece, clipthresh, clipniter, edge, threshold, chan_avg_limit, getresidual, showprogress, outlog, blfile, bltable); }
     301
     302  void sinusoidBaseline(const std::vector<bool>& mask, const std::string& fftinfo, const std::vector<int>& addwn, const std::vector<int>& rejwn, float clipthresh, int clipniter, bool getresidual=true, const std::string& showprogress="true,1000", const bool outlog=false, const std::string& blfile="", const std::string& bltable="")
     303  { table_->sinusoidBaseline(mask, fftinfo, addwn, rejwn, clipthresh, clipniter, getresidual, showprogress, outlog, blfile, bltable); }
     304
     305  void autoSinusoidBaseline(const std::vector<bool>& mask, const std::string& fftinfo, 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, const std::string& showprogress="true,1000", const bool outlog=false, const std::string& blfile="", const std::string& bltable="")
     306  { table_->autoSinusoidBaseline(mask, fftinfo, addwn, rejwn, clipthresh, clipniter, edge, threshold, chan_avg_limit, getresidual, showprogress, outlog, blfile, bltable); }
    303307
    304308  float getRms(const std::vector<bool>& mask, int whichrow)
  • trunk/src/python_Scantable.cpp

    r2713 r2767  
    145145    .def("_regrid_specchan", &ScantableWrapper::regridSpecChannel,
    146146         (boost::python::arg("nchan")=-1) )
     147    .def("_apply_bltable", &ScantableWrapper::applyBaselineTable)
     148    .def("_sub_baseline", &ScantableWrapper::subBaseline)
    147149    .def("_poly_baseline", &ScantableWrapper::polyBaseline)
    148150    .def("_auto_poly_baseline", &ScantableWrapper::autoPolyBaseline)
Note: See TracChangeset for help on using the changeset viewer.