Changeset 2047 for trunk


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

New Development: Yes

JIRA Issue: Yes CAS-2847

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: Yes

Module(s): scantable

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


Location:
trunk
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/python/asapfitter.py

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

  • trunk/python/scantable.py

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

  • trunk/src/STFitter.cpp

    r1932 r2047  
    3939#include <scimath/Functionals/Gaussian1D.h>
    4040#include "Lorentzian1D.h"
     41#include <scimath/Functionals/Sinusoid1D.h>
    4142#include <scimath/Functionals/Polynomial.h>
    4243#include <scimath/Mathematics/AutoDiff.h>
    … …  
    149150      funccomponents_.push_back(3);
    150151    }
    151   } else if (expr == "poly") {
    152     funcs_.resize(1);
    153     funcnames_.clear();
    154     funccomponents_.clear();
    155     funcs_[0] = new Polynomial<Float>(ncomp);
    156       funcnames_.push_back(expr);
    157       funccomponents_.push_back(ncomp);
    158152  } else if (expr == "lorentz") {
    159153    if (ncomp < 1) throw (AipsError("Need at least one lorentzian to fit."));
    … …  
    166160      funccomponents_.push_back(3);
    167161    }
     162  } else if (expr == "sinusoid") {
     163    if (ncomp < 1) throw (AipsError("Need at least one sinusoid to fit."));
     164    funcs_.resize(ncomp);
     165    funcnames_.clear();
     166    funccomponents_.clear();
     167    for (Int k=0; k<ncomp; ++k) {
     168      funcs_[k] = new Sinusoid1D<Float>();
     169      funcnames_.push_back(expr);
     170      funccomponents_.push_back(3);
     171    }
     172  } else if (expr == "poly") {
     173    funcs_.resize(1);
     174    funcnames_.clear();
     175    funccomponents_.clear();
     176    funcs_[0] = new Polynomial<Float>(ncomp);
     177      funcnames_.push_back(expr);
     178      funccomponents_.push_back(ncomp);
    168179  } else {
    169     //cerr << " compiled functions not yet implemented" << endl;
    170180    LogIO os( LogOrigin( "Fitter", "setExpression()", WHERE ) ) ;
    171181    os << LogIO::WARN << " compiled functions not yet implemented" << LogIO::POST;
    … …  
    244254            }
    245255        }
     256    } else if (dynamic_cast<Lorentzian1D<Float>* >(funcs_[0]) != 0) {
     257        uInt count = 0;
     258        for (uInt j=0; j < funcs_.nelements(); ++j) {
     259            for (uInt i=0; i < funcs_[j]->nparameters(); ++i) {
     260                (funcs_[j]->parameters())[i] = tmppar[count];
     261                parameters_[count] = tmppar[count];
     262                ++count;
     263            }
     264        }
     265    } else if (dynamic_cast<Sinusoid1D<Float>* >(funcs_[0]) != 0) {
     266        uInt count = 0;
     267        for (uInt j=0; j < funcs_.nelements(); ++j) {
     268            for (uInt i=0; i < funcs_[j]->nparameters(); ++i) {
     269                (funcs_[j]->parameters())[i] = tmppar[count];
     270                parameters_[count] = tmppar[count];
     271                ++count;
     272            }
     273        }
    246274    } else if (dynamic_cast<Polynomial<Float>* >(funcs_[0]) != 0) {
    247275        for (uInt i=0; i < funcs_[0]->nparameters(); ++i) {
    248276            parameters_[i] = tmppar[i];
    249277            (funcs_[0]->parameters())[i] =  tmppar[i];
    250         }
    251     } else if (dynamic_cast<Lorentzian1D<Float>* >(funcs_[0]) != 0) {
    252         uInt count = 0;
    253         for (uInt j=0; j < funcs_.nelements(); ++j) {
    254             for (uInt i=0; i < funcs_[j]->nparameters(); ++i) {
    255                 (funcs_[j]->parameters())[i] = tmppar[count];
    256                 parameters_[count] = tmppar[count];
    257                 ++count;
    258             }
    259278        }
    260279    }
    … …  
    286305            }
    287306        }
     307    } else if (dynamic_cast<Lorentzian1D<Float>* >(funcs_[0]) != 0) {
     308      uInt count = 0;
     309        for (uInt j=0; j < funcs_.nelements(); ++j) {
     310            for (uInt i=0; i < funcs_[j]->nparameters(); ++i) {
     311                funcs_[j]->mask(i) = !fixed[count];
     312                fixedpar_[count] = fixed[count];
     313                ++count;
     314            }
     315        }
     316    } else if (dynamic_cast<Sinusoid1D<Float>* >(funcs_[0]) != 0) {
     317      uInt count = 0;
     318        for (uInt j=0; j < funcs_.nelements(); ++j) {
     319            for (uInt i=0; i < funcs_[j]->nparameters(); ++i) {
     320                funcs_[j]->mask(i) = !fixed[count];
     321                fixedpar_[count] = fixed[count];
     322                ++count;
     323            }
     324        }
    288325    } else if (dynamic_cast<Polynomial<Float>* >(funcs_[0]) != 0) {
    289326        for (uInt i=0; i < funcs_[0]->nparameters(); ++i) {
    290327            fixedpar_[i] = fixed[i];
    291328            funcs_[0]->mask(i) =  !fixed[i];
    292         }
    293     } else if (dynamic_cast<Lorentzian1D<Float>* >(funcs_[0]) != 0) {
    294       uInt count = 0;
    295         for (uInt j=0; j < funcs_.nelements(); ++j) {
    296             for (uInt i=0; i < funcs_[j]->nparameters(); ++i) {
    297                 funcs_[j]->mask(i) = !fixed[count];
    298                 fixedpar_[count] = fixed[count];
    299                 ++count;
    300             }
    301329        }
    302330    }

  • trunk/src/Scantable.cpp

    r2036 r2047  
    882882    uInt row = uInt(whichrow);
    883883    stpol->setSpectra(getPolMatrix(row));
    884     Float fang,fhand,parang;
     884    Float fang,fhand;
    885885    fang = focusTable_.getTotalAngle(mfocusidCol_(row));
    886886    fhand = focusTable_.getFeedHand(mfocusidCol_(row));
    … …  
    17491749}
    17501750
    1751 bool Scantable::getFlagtraFast(int whichrow)
     1751bool Scantable::getFlagtraFast(uInt whichrow)
    17521752{
    17531753  uChar flag;
    17541754  Vector<uChar> flags;
    1755   flagsCol_.get(uInt(whichrow), flags);
     1755  flagsCol_.get(whichrow, flags);
    17561756  flag = flags[0];
    1757   for (int i = 1; i < flags.size(); ++i) {
     1757  for (uInt i = 1; i < flags.size(); ++i) {
    17581758    flag &= flags[i];
    17591759  }
    … …  
    17611761}
    17621762
    1763 void Scantable::cubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, bool outLogger, const std::string& blfile) {
     1763void Scantable::polyBaseline(const std::vector<bool>& mask, int order, bool outLogger, const std::string& blfile)
     1764{
    17641765  ofstream ofs;
    17651766  String coordInfo;
    1766   bool hasSameNchan;
    1767   int firstIF;
     1767  bool hasSameNchan = true;
    17681768  bool outTextFile = false;
    17691769
    … …  
    17771777    if (coordInfo == "") coordInfo = "channel";
    17781778    hasSameNchan = hasSameNchanOverIFs();
    1779     firstIF = getIF(0);
    1780   }
    1781 
    1782   //Fitter fitter = Fitter();
    1783   //fitter.setExpression("cspline", nPiece, thresClip, nIterClip);
     1779  }
     1780
     1781  Fitter fitter = Fitter();
     1782  fitter.setExpression("poly", order);
    17841783
    17851784  int nRow = nrow();
    … …  
    17871786
    17881787  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    1789 
    17901788    chanMask = getCompositeChanMask(whichrow, mask);
    1791     //fitBaseline(chanMask, whichrow, fitter);
    1792     //setSpectrum(fitter.getResidual(), whichrow);
    1793     std::vector<int> pieceRanges;
    1794     std::vector<float> params;
    1795     std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, pieceRanges, params, nPiece, thresClip, nIterClip);
    1796     setSpectrum(res, whichrow);
    1797 
    1798     if (outLogger || outTextFile) {
    1799       std::vector<bool>  fixed;
    1800       fixed.clear();
    1801       float rms = getRms(chanMask, whichrow);
    1802       String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan, firstIF, true);
    1803 
    1804       if (outLogger) {
    1805         LogIO ols(LogOrigin("Scantable", "cubicSplineBaseline()", WHERE));
    1806         ols << formatPiecewiseBaselineParams(pieceRanges, params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
    1807       }
    1808       if (outTextFile) {
    1809         ofs << formatPiecewiseBaselineParams(pieceRanges, params, fixed, rms, masklist, whichrow, true) << flush;
    1810       }
    1811     }
    1812 
     1789    fitBaseline(chanMask, whichrow, fitter);
     1790    setSpectrum(fitter.getResidual(), whichrow);
     1791    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "polyBaseline()", fitter);
    18131792  }
    18141793
    … …  
    18161795}
    18171796
    1818 
    1819 void Scantable::autoCubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool outLogger, const std::string& blfile)
     1797void Scantable::autoPolyBaseline(const std::vector<bool>& mask, int order, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool outLogger, const std::string& blfile)
    18201798{
    18211799  ofstream ofs;
    18221800  String coordInfo;
    1823   bool hasSameNchan;
    1824   int firstIF;
     1801  bool hasSameNchan = true;
    18251802  bool outTextFile = false;
    18261803
    … …  
    18341811    if (coordInfo == "") coordInfo = "channel";
    18351812    hasSameNchan = hasSameNchanOverIFs();
    1836     firstIF = getIF(0);
    1837   }
    1838 
    1839   //Fitter fitter = Fitter();
    1840   //fitter.setExpression("cspline", nPiece, thresClip, nIterClip);
     1813  }
     1814
     1815  Fitter fitter = Fitter();
     1816  fitter.setExpression("poly", order);
    18411817
    18421818  int nRow = nrow();
    … …  
    18711847    //-------------------------------------------------------
    18721848
    1873 
    1874     //fitBaseline(chanMask, whichrow, fitter);
     1849    fitBaseline(chanMask, whichrow, fitter);
     1850    setSpectrum(fitter.getResidual(), whichrow);
     1851
     1852    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoPolyBaseline()", fitter);
     1853  }
     1854
     1855  if (outTextFile) ofs.close();
     1856}
     1857
     1858void Scantable::cubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, bool outLogger, const std::string& blfile) {
     1859  ofstream ofs;
     1860  String coordInfo;
     1861  bool hasSameNchan = true;
     1862  bool outTextFile = false;
     1863
     1864  if (blfile != "") {
     1865    ofs.open(blfile.c_str(), ios::out | ios::app);
     1866    if (ofs) outTextFile = true;
     1867  }
     1868
     1869  if (outLogger || outTextFile) {
     1870    coordInfo = getCoordInfo()[0];
     1871    if (coordInfo == "") coordInfo = "channel";
     1872    hasSameNchan = hasSameNchanOverIFs();
     1873  }
     1874
     1875  //Fitter fitter = Fitter();
     1876  //fitter.setExpression("cspline", nPiece);
     1877
     1878  int nRow = nrow();
     1879  std::vector<bool> chanMask;
     1880
     1881  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
     1882    chanMask = getCompositeChanMask(whichrow, mask);
     1883    //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
    18751884    //setSpectrum(fitter.getResidual(), whichrow);
    18761885    std::vector<int> pieceRanges;
    … …  
    18781887    std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, pieceRanges, params, nPiece, thresClip, nIterClip);
    18791888    setSpectrum(res, whichrow);
    1880 
    1881     if (outLogger || outTextFile) {
    1882       std::vector<bool> fixed;
    1883       fixed.clear();
    1884       float rms = getRms(chanMask, whichrow);
    1885       String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan, firstIF, true);
    1886 
    1887       if (outLogger) {
    1888         LogIO ols(LogOrigin("Scantable", "autoCubicSplineBaseline()", WHERE));
    1889         ols << formatPiecewiseBaselineParams(pieceRanges, params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
    1890       }
    1891       if (outTextFile) {
    1892         ofs << formatPiecewiseBaselineParams(pieceRanges, params, fixed, rms, masklist, whichrow, true) << flush;
    1893       }
    1894     }
    1895 
     1889    //
     1890
     1891    std::vector<bool>  fixed;
     1892    fixed.clear();
     1893    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()", pieceRanges, params, fixed);
    18961894  }
    18971895
    … …  
    18991897}
    19001898
     1899void Scantable::autoCubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool outLogger, const std::string& blfile)
     1900{
     1901  ofstream ofs;
     1902  String coordInfo;
     1903  bool hasSameNchan = true;
     1904  bool outTextFile = false;
     1905
     1906  if (blfile != "") {
     1907    ofs.open(blfile.c_str(), ios::out | ios::app);
     1908    if (ofs) outTextFile = true;
     1909  }
     1910
     1911  if (outLogger || outTextFile) {
     1912    coordInfo = getCoordInfo()[0];
     1913    if (coordInfo == "") coordInfo = "channel";
     1914    hasSameNchan = hasSameNchanOverIFs();
     1915  }
     1916
     1917  //Fitter fitter = Fitter();
     1918  //fitter.setExpression("cspline", nPiece);
     1919
     1920  int nRow = nrow();
     1921  std::vector<bool> chanMask;
     1922  int minEdgeSize = getIFNos().size()*2;
     1923  STLineFinder lineFinder = STLineFinder();
     1924  lineFinder.setOptions(threshold, 3, chanAvgLimit);
     1925
     1926  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
     1927
     1928    //-------------------------------------------------------
     1929    //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
     1930    //-------------------------------------------------------
     1931    int edgeSize = edge.size();
     1932    std::vector<int> currentEdge;
     1933    if (edgeSize >= 2) {
     1934      int idx = 0;
     1935      if (edgeSize > 2) {
     1936        if (edgeSize < minEdgeSize) {
     1937          throw(AipsError("Length of edge element info is less than that of IFs"));
     1938        }
     1939        idx = 2 * getIF(whichrow);
     1940      }
     1941      currentEdge.push_back(edge[idx]);
     1942      currentEdge.push_back(edge[idx+1]);
     1943    } else {
     1944      throw(AipsError("Wrong length of edge element"));
     1945    }
     1946    lineFinder.setData(getSpectrum(whichrow));
     1947    lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
     1948    chanMask = lineFinder.getMask();
     1949    //-------------------------------------------------------
     1950
     1951
     1952    //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
     1953    //setSpectrum(fitter.getResidual(), whichrow);
     1954    std::vector<int> pieceRanges;
     1955    std::vector<float> params;
     1956    std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, pieceRanges, params, nPiece, thresClip, nIterClip);
     1957    setSpectrum(res, whichrow);
     1958    //
     1959
     1960    std::vector<bool> fixed;
     1961    fixed.clear();
     1962    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()", pieceRanges, params, fixed);
     1963  }
     1964
     1965  if (outTextFile) ofs.close();
     1966}
    19011967
    19021968std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data, const std::vector<bool>& mask, std::vector<int>& sectionRanges, std::vector<float>& params, int nPiece, float thresClip, int nIterClip) {
    … …  
    19211987  int nElement = (int)(floor(floor((double)(nData/nPiece))+0.5));
    19221988
    1923   std::vector<double> sectionList0, sectionList1, sectionListR;
    1924   sectionList0.push_back((double)x[0]);
     1989  std::vector<int> sectionList0, sectionList1;
     1990  std::vector<double> sectionListR;
     1991  sectionList0.push_back(x[0]);
    19251992  sectionRanges.clear();
    19261993  sectionRanges.push_back(x[0]);
    19271994  for (int i = 1; i < nPiece; ++i) {
    1928     double valX = (double)x[nElement*i];
     1995    int valX = x[nElement*i];
    19291996    sectionList0.push_back(valX);
    19301997    sectionList1.push_back(valX);
    1931     sectionListR.push_back(1.0/valX);
    1932 
    1933     sectionRanges.push_back(x[nElement*i]-1);
    1934     sectionRanges.push_back(x[nElement*i]);
    1935   }
    1936   sectionList1.push_back((double)(x[x.size()-1]+1));
     1998    sectionListR.push_back(1.0/(double)valX);
     1999
     2000    sectionRanges.push_back(valX-1);
     2001    sectionRanges.push_back(valX);
     2002  }
     2003  sectionList1.push_back(x[x.size()-1]+1);
    19372004  sectionRanges.push_back(x[x.size()-1]);
    19382005
    … …  
    21112178}
    21122179
    2113 void Scantable::fitBaseline(const std::vector<bool>& mask, int whichrow, Fitter& fitter)
    2114 {
    2115   std::vector<double> abcsd = getAbcissa(whichrow);
    2116   std::vector<float> abcs;
    2117   for (int i = 0; i < abcsd.size(); ++i) {
    2118     abcs.push_back((float)abcsd[i]);
    2119   }
    2120   std::vector<float> spec = getSpectrum(whichrow);
    2121 
    2122   fitter.setData(abcs, spec, mask);
    2123   fitter.lfit();
    2124 }
    2125 
    2126 std::vector<bool> Scantable::getCompositeChanMask(int whichrow, const std::vector<bool>& inMask)
    2127 {
    2128   std::vector<bool> chanMask = getMask(whichrow);
    2129   int chanMaskSize = chanMask.size();
    2130   if (chanMaskSize != inMask.size()) {
    2131     throw(AipsError("different mask sizes"));
    2132   }
    2133   for (int i = 0; i < chanMaskSize; ++i) {
    2134     chanMask[i] = chanMask[i] && inMask[i];
    2135   }
    2136 
    2137   return chanMask;
    2138 }
    2139 
    2140 /*
    2141 std::vector<bool> Scantable::getCompositeChanMask(int whichrow, const std::vector<bool>& inMask, const std::vector<int>& edge, const int minEdgeSize, STLineFinder& lineFinder)
    2142 {
    2143   int edgeSize = edge.size();
    2144   std::vector<int> currentEdge;
    2145   if (edgeSize >= 2) {
    2146       int idx = 0;
    2147       if (edgeSize > 2) {
    2148         if (edgeSize < minEdgeSize) {
    2149           throw(AipsError("Length of edge element info is less than that of IFs"));
    2150         }
    2151         idx = 2 * getIF(whichrow);
    2152       }
    2153       currentEdge.push_back(edge[idx]);
    2154       currentEdge.push_back(edge[idx+1]);
    2155   } else {
    2156     throw(AipsError("Wrong length of edge element"));
    2157   }
    2158 
    2159   lineFinder.setData(getSpectrum(whichrow));
    2160   lineFinder.findLines(getCompositeChanMask(whichrow, inMask), currentEdge, whichrow);
    2161 
    2162   return lineFinder.getMask();
    2163 }
    2164 */
    2165 
    2166 void Scantable::polyBaseline(const std::vector<bool>& mask, int order, bool outLogger, const std::string& blfile)
    2167 {
     2180  void Scantable::sinusoidBaseline(const std::vector<bool>& mask, int nMinWavesInSW, int nMaxWavesInSW, float thresClip, int nIterClip, bool outLogger, const std::string& blfile) {
    21682181  ofstream ofs;
    21692182  String coordInfo;
    2170   bool hasSameNchan;
    2171   int firstIF;
     2183  bool hasSameNchan = true;
    21722184  bool outTextFile = false;
    21732185
    … …  
    21812193    if (coordInfo == "") coordInfo = "channel";
    21822194    hasSameNchan = hasSameNchanOverIFs();
    2183     firstIF = getIF(0);
    2184   }
    2185 
    2186   Fitter fitter = Fitter();
    2187   fitter.setExpression("poly", order);
     2195  }
     2196
     2197  //Fitter fitter = Fitter();
     2198  //fitter.setExpression("sinusoid", nMaxWavesInSW);
    21882199
    21892200  int nRow = nrow();
    … …  
    21912202
    21922203  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    2193 
    21942204    chanMask = getCompositeChanMask(whichrow, mask);
    2195     fitBaseline(chanMask, whichrow, fitter);
    2196     setSpectrum(fitter.getResidual(), whichrow);
    2197    
    2198     if (outLogger || outTextFile) {
    2199       std::vector<float> params = fitter.getParameters();
    2200       std::vector<bool>  fixed  = fitter.getFixedParameters();
    2201       float rms = getRms(chanMask, whichrow);
    2202       String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan, firstIF, true);
    2203 
    2204       if (outLogger) {
    2205         LogIO ols(LogOrigin("Scantable", "polyBaseline()", WHERE));
    2206         ols << formatBaselineParams(params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
    2207       }
    2208       if (outTextFile) {
    2209         ofs << formatBaselineParams(params, fixed, rms, masklist, whichrow, true) << flush;
    2210       }
    2211     }
    2212 
     2205    //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
     2206    //setSpectrum(fitter.getResidual(), whichrow);
     2207    std::vector<int> pieceRanges;
     2208    std::vector<float> params;
     2209    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nMinWavesInSW, nMaxWavesInSW, params, thresClip, nIterClip);
     2210    setSpectrum(res, whichrow);
     2211    //
     2212
     2213    std::vector<bool>  fixed;
     2214    fixed.clear();
     2215    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "sinusoidBaseline()", pieceRanges, params, fixed);
    22132216  }
    22142217
    … …  
    22162219}
    22172220
    2218 void Scantable::autoPolyBaseline(const std::vector<bool>& mask, int order, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool outLogger, const std::string& blfile)
     2221void Scantable::autoSinusoidBaseline(const std::vector<bool>& mask, int nMinWavesInSW, int nMaxWavesInSW, float thresClip, int nIterClip, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool outLogger, const std::string& blfile)
    22192222{
    22202223  ofstream ofs;
    22212224  String coordInfo;
    2222   bool hasSameNchan;
    2223   int firstIF;
     2225  bool hasSameNchan = true;
    22242226  bool outTextFile = false;
    22252227
    … …  
    22332235    if (coordInfo == "") coordInfo = "channel";
    22342236    hasSameNchan = hasSameNchanOverIFs();
    2235     firstIF = getIF(0);
    2236   }
    2237 
    2238   Fitter fitter = Fitter();
    2239   fitter.setExpression("poly", order);
     2237  }
     2238
     2239  //Fitter fitter = Fitter();
     2240  //fitter.setExpression("sinusoid", nMaxWavesInSW);
    22402241
    22412242  int nRow = nrow();
    … …  
    22712272
    22722273
    2273     fitBaseline(chanMask, whichrow, fitter);
    2274     setSpectrum(fitter.getResidual(), whichrow);
    2275 
    2276     if (outLogger || outTextFile) {
    2277       std::vector<float> params = fitter.getParameters();
    2278       std::vector<bool>  fixed  = fitter.getFixedParameters();
    2279       float rms = getRms(chanMask, whichrow);
    2280       String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan, firstIF, true);
    2281 
    2282       if (outLogger) {
    2283         LogIO ols(LogOrigin("Scantable", "autoPolyBaseline()", WHERE));
    2284         ols << formatBaselineParams(params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
    2285       }
    2286       if (outTextFile) {
    2287         ofs << formatBaselineParams(params, fixed, rms, masklist, whichrow, true) << flush;
    2288       }
    2289     }
    2290 
     2274    //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
     2275    //setSpectrum(fitter.getResidual(), whichrow);
     2276    std::vector<int> pieceRanges;
     2277    std::vector<float> params;
     2278    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nMinWavesInSW, nMaxWavesInSW, params, thresClip, nIterClip);
     2279    setSpectrum(res, whichrow);
     2280    //
     2281
     2282    std::vector<bool> fixed;
     2283    fixed.clear();
     2284    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoSinusoidBaseline()", pieceRanges, params, fixed);
    22912285  }
    22922286
    … …  
    22942288}
    22952289
     2290std::vector<float> Scantable::doSinusoidFitting(const std::vector<float>& data, const std::vector<bool>& mask, int nMinWavesInSW, int nMaxWavesInSW, std::vector<float>& params, float thresClip, int nIterClip) {
     2291  if (nMaxWavesInSW < 1) {
     2292    throw(AipsError("maximum wave number must be a positive value"));
     2293  }
     2294  if (data.size() != mask.size()) {
     2295    throw(AipsError("data and mask have different size"));
     2296  }
     2297
     2298  int nChan = data.size();
     2299  std::vector<int> maskArray;
     2300  std::vector<int> x;
     2301  for (int i = 0; i < nChan; ++i) {
     2302    maskArray.push_back(mask[i] ? 1 : 0);
     2303    if (mask[i]) {
     2304      x.push_back(i);
     2305    }
     2306  }
     2307
     2308  int nData = x.size();
     2309
     2310  std::vector<double> x1, x2, x3, z1, x1z1, x2z1, x3z1, r1;
     2311  for (int i = 0; i < nChan; ++i) {
     2312    x1.push_back((double)i);
     2313    x2.push_back((double)(i*i));
     2314    x3.push_back((double)(i*i*i));
     2315    z1.push_back((double)data[i]);
     2316    x1z1.push_back(((double)i)*(double)data[i]);
     2317    x2z1.push_back(((double)(i*i))*(double)data[i]);
     2318    x3z1.push_back(((double)(i*i*i))*(double)data[i]);
     2319    r1.push_back(0.0);
     2320  }
     2321
     2322  int currentNData = nData;
     2323  int nDOF = nMaxWavesInSW + 1;  //number of parameters to solve.
     2324
     2325  for (int nClip = 0; nClip < nIterClip+1; ++nClip) {
     2326    //xMatrix : horizontal concatenation of the least-sq. matrix (left) and an
     2327    //identity matrix (right).
     2328    //the right part is used to calculate the inverse matrix of the left part.
     2329    double xMatrix[nDOF][2*nDOF];
     2330    double zMatrix[nDOF];
     2331    for (int i = 0; i < nDOF; ++i) {
     2332      for (int j = 0; j < 2*nDOF; ++j) {
     2333        xMatrix[i][j] = 0.0;
     2334      }
     2335      xMatrix[i][nDOF+i] = 1.0;
     2336      zMatrix[i] = 0.0;
     2337    }
     2338
     2339    for (int i = 0; i < currentNData; ++i) {
     2340      if (maskArray[i] == 0) continue;
     2341      xMatrix[0][0] += 1.0;
     2342      xMatrix[0][1] += x1[i];
     2343      xMatrix[0][2] += x2[i];
     2344      xMatrix[0][3] += x3[i];
     2345      xMatrix[1][1] += x2[i];
     2346      xMatrix[1][2] += x3[i];
     2347      xMatrix[1][3] += x2[i]*x2[i];
     2348      xMatrix[2][2] += x2[i]*x2[i];
     2349      xMatrix[2][3] += x3[i]*x2[i];
     2350      xMatrix[3][3] += x3[i]*x3[i];
     2351      zMatrix[0] += z1[i];
     2352      zMatrix[1] += x1z1[i];
     2353      zMatrix[2] += x2z1[i];
     2354      zMatrix[3] += x3z1[i];
     2355    }
     2356
     2357    for (int i = 0; i < nDOF; ++i) {
     2358      for (int j = 0; j < i; ++j) {
     2359        xMatrix[i][j] = xMatrix[j][i];
     2360      }
     2361    }
     2362
     2363    std::vector<double> invDiag;
     2364    for (int i = 0; i < nDOF; ++i) {
     2365      invDiag.push_back(1.0/xMatrix[i][i]);
     2366      for (int j = 0; j < nDOF; ++j) {
     2367        xMatrix[i][j] *= invDiag[i];
     2368      }
     2369    }
     2370
     2371    for (int k = 0; k < nDOF; ++k) {
     2372      for (int i = 0; i < nDOF; ++i) {
     2373        if (i != k) {
     2374          double factor1 = xMatrix[k][k];
     2375          double factor2 = xMatrix[i][k];
     2376          for (int j = k; j < 2*nDOF; ++j) {
     2377            xMatrix[i][j] *= factor1;
     2378            xMatrix[i][j] -= xMatrix[k][j]*factor2;
     2379            xMatrix[i][j] /= factor1;
     2380          }
     2381        }
     2382      }
     2383      double xDiag = xMatrix[k][k];
     2384      for (int j = k; j < 2*nDOF; ++j) {
     2385        xMatrix[k][j] /= xDiag;
     2386      }
     2387    }
     2388   
     2389    for (int i = 0; i < nDOF; ++i) {
     2390      for (int j = 0; j < nDOF; ++j) {
     2391        xMatrix[i][nDOF+j] *= invDiag[j];
     2392      }
     2393    }
     2394    //compute a vector y which consists of the coefficients of the sinusoids forming the
     2395    //best-fit curves (a0,b0,a1,b1,...), namely, a* and b* are of sine and cosine functions,
     2396    //respectively.
     2397    std::vector<double> y;
     2398    for (int i = 0; i < nDOF; ++i) {
     2399      y.push_back(0.0);
     2400      for (int j = 0; j < nDOF; ++j) {
     2401        y[i] += xMatrix[i][nDOF+j]*zMatrix[j];
     2402      }
     2403    }
     2404
     2405    double a0 = y[0];
     2406    double a1 = y[1];
     2407    double a2 = y[2];
     2408    double a3 = y[3];
     2409    params.clear();
     2410
     2411    for (int i = 0; i < nChan; ++i) {
     2412      r1[i] = a0 + a1*x1[i] + a2*x2[i] + a3*x3[i];
     2413    }
     2414    params.push_back(a0);
     2415    params.push_back(a1);
     2416    params.push_back(a2);
     2417    params.push_back(a3);
     2418
     2419
     2420    if ((nClip == nIterClip) || (thresClip <= 0.0)) {
     2421      break;
     2422    } else {
     2423      std::vector<double> rz;
     2424      double stdDev = 0.0;
     2425      for (int i = 0; i < nChan; ++i) {
     2426        double val = abs(r1[i] - z1[i]);
     2427        rz.push_back(val);
     2428        stdDev += val*val*(double)maskArray[i];
     2429      }
     2430      stdDev = sqrt(stdDev/(double)nData);
     2431     
     2432      double thres = stdDev * thresClip;
     2433      int newNData = 0;
     2434      for (int i = 0; i < nChan; ++i) {
     2435        if (rz[i] >= thres) {
     2436          maskArray[i] = 0;
     2437        }
     2438        if (maskArray[i] > 0) {
     2439          newNData++;
     2440        }
     2441      }
     2442      if (newNData == currentNData) {
     2443        break;                   //no additional flag. finish iteration.
     2444      } else {
     2445        currentNData = newNData;
     2446      }
     2447    }
     2448  }
     2449
     2450  std::vector<float> residual;
     2451  for (int i = 0; i < nChan; ++i) {
     2452    residual.push_back((float)(z1[i] - r1[i]));
     2453  }
     2454  return residual;
     2455
     2456}
     2457
     2458void Scantable::fitBaseline(const std::vector<bool>& mask, int whichrow, Fitter& fitter)
     2459{
     2460  std::vector<double> abcsd = getAbcissa(whichrow);
     2461  std::vector<float> abcs;
     2462  for (uInt i = 0; i < abcsd.size(); ++i) {
     2463    abcs.push_back((float)abcsd[i]);
     2464  }
     2465  std::vector<float> spec = getSpectrum(whichrow);
     2466
     2467  fitter.setData(abcs, spec, mask);
     2468  fitter.lfit();
     2469}
     2470
     2471std::vector<bool> Scantable::getCompositeChanMask(int whichrow, const std::vector<bool>& inMask)
     2472{
     2473  std::vector<bool> chanMask = getMask(whichrow);
     2474  uInt chanMaskSize = chanMask.size();
     2475  if (chanMaskSize != inMask.size()) {
     2476    throw(AipsError("different mask sizes"));
     2477  }
     2478  for (uInt i = 0; i < chanMaskSize; ++i) {
     2479    chanMask[i] = chanMask[i] && inMask[i];
     2480  }
     2481
     2482  return chanMask;
     2483}
     2484
     2485/*
     2486std::vector<bool> Scantable::getCompositeChanMask(int whichrow, const std::vector<bool>& inMask, const std::vector<int>& edge, const int minEdgeSize, STLineFinder& lineFinder)
     2487{
     2488  int edgeSize = edge.size();
     2489  std::vector<int> currentEdge;
     2490  if (edgeSize >= 2) {
     2491      int idx = 0;
     2492      if (edgeSize > 2) {
     2493        if (edgeSize < minEdgeSize) {
     2494          throw(AipsError("Length of edge element info is less than that of IFs"));
     2495        }
     2496        idx = 2 * getIF(whichrow);
     2497      }
     2498      currentEdge.push_back(edge[idx]);
     2499      currentEdge.push_back(edge[idx+1]);
     2500  } else {
     2501    throw(AipsError("Wrong length of edge element"));
     2502  }
     2503
     2504  lineFinder.setData(getSpectrum(whichrow));
     2505  lineFinder.findLines(getCompositeChanMask(whichrow, inMask), currentEdge, whichrow);
     2506
     2507  return lineFinder.getMask();
     2508}
     2509*/
     2510
     2511/* for poly. the variations of outputFittingResult() should be merged into one eventually (2011/3/10 WK)  */
     2512void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, Fitter& fitter) {
     2513  if (outLogger || outTextFile) {
     2514    std::vector<float> params = fitter.getParameters();
     2515    std::vector<bool>  fixed  = fitter.getFixedParameters();
     2516    float rms = getRms(chanMask, whichrow);
     2517    String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
     2518
     2519    if (outLogger) {
     2520      LogIO ols(LogOrigin("Scantable", funcName, WHERE));
     2521      ols << formatBaselineParams(params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
     2522    }
     2523    if (outTextFile) {
     2524      ofs << formatBaselineParams(params, fixed, rms, masklist, whichrow, true) << flush;
     2525    }
     2526  }
     2527}
     2528
     2529/* for cspline. will be merged once cspline is available in fitter (2011/3/10 WK) */
     2530void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, const std::vector<int>& pieceEdges, const std::vector<float>& params, const std::vector<bool>& fixed) {
     2531  if (outLogger || outTextFile) {
     2532    float rms = getRms(chanMask, whichrow);
     2533    String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
     2534
     2535    if (outLogger) {
     2536      LogIO ols(LogOrigin("Scantable", funcName, WHERE));
     2537      ols << formatPiecewiseBaselineParams(pieceEdges, params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
     2538    }
     2539    if (outTextFile) {
     2540      ofs << formatPiecewiseBaselineParams(pieceEdges, params, fixed, rms, masklist, whichrow, true) << flush;
     2541    }
     2542  }
     2543}
     2544
     2545/* for sinusoid. will be merged once sinusoid is available in fitter (2011/3/10 WK) */
     2546void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, const std::vector<float>& params, const std::vector<bool>& fixed) {
     2547  if (outLogger || outTextFile) {
     2548    float rms = getRms(chanMask, whichrow);
     2549    String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
     2550
     2551    if (outLogger) {
     2552      LogIO ols(LogOrigin("Scantable", funcName, WHERE));
     2553      ols << formatBaselineParams(params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
     2554    }
     2555    if (outTextFile) {
     2556      ofs << formatBaselineParams(params, fixed, rms, masklist, whichrow, true) << flush;
     2557    }
     2558  }
     2559}
    22962560
    22972561float Scantable::getRms(const std::vector<bool>& mask, int whichrow) {
    … …  
    23022566  float smean = 0.0;
    23032567  int n = 0;
    2304   for (int i = 0; i < spec.nelements(); ++i) {
     2568  for (uInt i = 0; i < spec.nelements(); ++i) {
    23052569    if (mask[i]) {
    23062570      mean += spec[i];
    … …  
    23442608
    23452609  if (verbose) {
    2346     oss << endl;
    23472610    oss << "Results of baseline fit" << endl;
    23482611    oss << "  rms = " << setprecision(6) << rms << endl;
    … …  
    23522615}
    23532616
    2354 std::string Scantable::formatPiecewiseBaselineParams(const std::vector<int>& ranges, const std::vector<float>& params, const std::vector<bool>& fixed, float rms, const std::string& masklist, int whichrow, bool verbose) const
     2617std::string Scantable::formatBaselineParams(const std::vector<float>& params, const std::vector<bool>& fixed, float rms, const std::string& masklist, int whichrow, bool verbose) const
    23552618{
    23562619  ostringstream oss;
    23572620  oss << formatBaselineParamsHeader(whichrow, masklist, verbose);
    23582621
     2622  if (params.size() > 0) {
     2623    for (uInt i = 0; i < params.size(); ++i) {
     2624      if (i > 0) {
     2625        oss << ",";
     2626      }
     2627      std::string fix = ((fixed.size() > 0) && (fixed[i]) && verbose) ? "(fixed)" : "";
     2628      float vParam = params[i];
     2629      std::string equals = "= " + (String)((vParam < 0.0) ? "" : " ");
     2630      oss << "  p" << i << fix << equals << setprecision(6) << vParam;
     2631    }
     2632  } else {
     2633    oss << "  Not fitted";
     2634  }
     2635  oss << endl;
     2636
     2637  oss << formatBaselineParamsFooter(rms, verbose);
     2638  oss << flush;
     2639
     2640  return String(oss);
     2641}
     2642
     2643std::string Scantable::formatPiecewiseBaselineParams(const std::vector<int>& ranges, const std::vector<float>& params, const std::vector<bool>& fixed, float rms, const std::string& masklist, int whichrow, bool verbose) const
     2644{
     2645  ostringstream oss;
     2646  oss << formatBaselineParamsHeader(whichrow, masklist, verbose);
     2647
    23592648  if ((params.size() > 0) && (ranges.size() > 0)) {
    23602649    if (((ranges.size() % 2) == 0) && ((params.size() % (ranges.size() / 2)) == 0)) {
    2361       int nParam = params.size() / (ranges.size() / 2);
     2650      uInt nParam = params.size() / (ranges.size() / 2);
     2651      stringstream ss;
     2652      ss << ranges[ranges.size()-1] << flush;
     2653      int wRange = ss.str().size()*2+5;
     2654      ss.str("");
    23622655      for (uInt j = 0; j < ranges.size(); j+=2) {
    2363         oss << "[" << ranges[j] << "," << ranges[j+1] << "]";
     2656        ss << "  [" << ranges[j] << "," << ranges[j+1] << "]" << flush;
     2657        oss << setw(wRange) << left << ss.str();
     2658        ss.str("");
    23642659        for (uInt i = 0; i < nParam; ++i) {
    23652660          if (i > 0) {
    23662661            oss << ",";
    23672662          }
    2368           String fix = ((fixed.size() > 0) && (fixed[i]) && verbose) ? "(fixed)" : "";
    2369           oss << "  p" << i << fix << "= " << setprecision(6) << params[j/2*nParam+i];
     2663          std::string fix = ((fixed.size() > 0) && (fixed[i]) && verbose) ? "(fixed)" : "";
     2664          float vParam = params[j/2*nParam+i];
     2665          std::string equals = "= " + (String)((vParam < 0.0) ? "" : " ");
     2666          oss << "  p" << i << fix << equals << setprecision(6) << vParam;
    23702667        }
     2668        oss << endl;
    23712669      }
    23722670    } else {
    2373       oss << "  ";
     2671      oss << "  " << endl;
    23742672    }
    23752673  } else {
    2376     oss << "  Not fitted";
     2674    oss << "  Not fitted" << endl;
    23772675  }
    23782676
    23792677  oss << formatBaselineParamsFooter(rms, verbose);
    23802678  oss << flush;
    2381 
    2382   return String(oss);
    2383 }
    2384 
    2385 std::string Scantable::formatBaselineParams(const std::vector<float>& params, const std::vector<bool>& fixed, float rms, const std::string& masklist, int whichrow, bool verbose) const
    2386 {
    2387   ostringstream oss;
    2388   oss << formatBaselineParamsHeader(whichrow, masklist, verbose);
    2389 
    2390   if (params.size() > 0) {
    2391     for (uInt i = 0; i < params.size(); ++i) {
    2392       if (i > 0) {
    2393         oss << ",";
    2394       }
    2395       String fix = ((fixed.size() > 0) && (fixed[i]) && verbose) ? "(fixed)" : "";
    2396       oss << "  p" << i << fix << "= " << setprecision(6) << params[i];
    2397     }
    2398   } else {
    2399     oss << "  Not fitted";
    2400   }
    2401 
    2402   oss << formatBaselineParamsFooter(rms, verbose);
    2403   oss << flush;
    2404 
    2405   return String(oss);
    2406 }
    2407 
    2408 std::string Scantable::getMaskRangeList(const std::vector<bool>& mask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, int firstIF, bool silent)
    2409 {
    2410   if (mask.size() < 2) {
    2411     throw(AipsError("The mask elements should be > 1"));
    2412   }
    2413   if (mask.size() != nchan()) {
    2414     throw(AipsError("Number of channels in scantable != number of mask elements"));
    2415   }
    2416 
    2417   std::vector<int> startIdx = getMaskEdgeIndices(mask, true);
    2418   std::vector<int> endIdx   = getMaskEdgeIndices(mask, false);
    2419 
    2420   if (startIdx.size() != endIdx.size()) {
    2421     throw(AipsError("Numbers of mask start and mask end are not identical"));
    2422   }
    2423   for (int i = 0; i < startIdx.size(); ++i) {
    2424     if (startIdx[i] > endIdx[i]) {
    2425       throw(AipsError("Mask start index > mask end index"));
    2426     }
    2427   }
    2428 
    2429   if (!silent) {
    2430     LogIO logOs(LogOrigin("Scantable", "getMaskRangeList()", WHERE));
    2431     logOs << LogIO::WARN << "The current mask window unit is " << coordInfo;
    2432     if (!hasSameNchan) {
    2433       logOs << endl << "This mask is only valid for IF=" << firstIF;
    2434     }
    2435     logOs << LogIO::POST;
    2436   }
    2437 
    2438   std::vector<double> abcissa = getAbcissa(whichrow);
    2439   ostringstream oss;
    2440   oss.setf(ios::fixed);
    2441   oss << setprecision(1) << "[";
    2442   for (int i = 0; i < startIdx.size(); ++i) {
    2443     std::vector<float> aMaskRange;
    2444     aMaskRange.push_back((float)abcissa[startIdx[i]]);
    2445     aMaskRange.push_back((float)abcissa[endIdx[i]]);
    2446     sort(aMaskRange.begin(), aMaskRange.end());
    2447     if (i > 0) oss << ",";
    2448     oss << "[" << aMaskRange[0] << "," << aMaskRange[1] << "]";
    2449   }
    2450   oss << "]" << flush;
    24512679
    24522680  return String(oss);
    … …  
    24712699}
    24722700
    2473 std::vector<int> Scantable::getMaskEdgeIndices(const std::vector<bool>& mask, bool getStartIndices)
     2701std::string Scantable::getMaskRangeList(const std::vector<bool>& mask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, bool verbose)
    24742702{
    24752703  if (mask.size() < 2) {
    24762704    throw(AipsError("The mask elements should be > 1"));
    24772705  }
    2478   if (mask.size() != nchan()) {
     2706  int IF = getIF(whichrow);
     2707  if (mask.size() != (uInt)nchan(IF)) {
    24792708    throw(AipsError("Number of channels in scantable != number of mask elements"));
    24802709  }
    24812710
    2482   std::vector<int> out;
    2483   int endIdx = mask.size() - 1;
    2484 
    2485   if (getStartIndices) {
    2486     if (mask[0]) {
    2487       out.push_back(0);
    2488     }
    2489     for (int i = 0; i < endIdx; ++i) {
    2490       if ((!mask[i]) && mask[i+1]) {
    2491         out.push_back(i+1);
    2492       }
    2493     }
    2494   } else {
    2495     for (int i = 0; i < endIdx; ++i) {
    2496       if (mask[i] && (!mask[i+1])) {
    2497         out.push_back(i);
    2498       }
    2499     }
    2500     if (mask[endIdx]) {
    2501       out.push_back(endIdx);
    2502     }
     2711  if (verbose) {
     2712    LogIO logOs(LogOrigin("Scantable", "getMaskRangeList()", WHERE));
     2713    logOs << LogIO::WARN << "The current mask window unit is " << coordInfo;
     2714    if (!hasSameNchan) {
     2715      logOs << endl << "This mask is only valid for IF=" << IF;
     2716    }
     2717    logOs << LogIO::POST;
     2718  }
     2719
     2720  std::vector<double> abcissa = getAbcissa(whichrow);
     2721  std::vector<int> edge = getMaskEdgeIndices(mask);
     2722
     2723  ostringstream oss;
     2724  oss.setf(ios::fixed);
     2725  oss << setprecision(1) << "[";
     2726  for (uInt i = 0; i < edge.size(); i+=2) {
     2727    if (i > 0) oss << ",";
     2728    oss << "[" << (float)abcissa[edge[i]] << "," << (float)abcissa[edge[i+1]] << "]";
     2729  }
     2730  oss << "]" << flush;
     2731
     2732  return String(oss);
     2733}
     2734
     2735std::vector<int> Scantable::getMaskEdgeIndices(const std::vector<bool>& mask)
     2736{
     2737  if (mask.size() < 2) {
     2738    throw(AipsError("The mask elements should be > 1"));
     2739  }
     2740
     2741  std::vector<int> out, startIndices, endIndices;
     2742  int maskSize = mask.size();
     2743
     2744  startIndices.clear();
     2745  endIndices.clear();
     2746
     2747  if (mask[0]) {
     2748    startIndices.push_back(0);
     2749  }
     2750  for (int i = 1; i < maskSize; ++i) {
     2751    if ((!mask[i-1]) && mask[i]) {
     2752      startIndices.push_back(i);
     2753    } else if (mask[i-1] && (!mask[i])) {
     2754      endIndices.push_back(i-1);
     2755    }
     2756  }
     2757  if (mask[maskSize-1]) {
     2758    endIndices.push_back(maskSize-1);
     2759  }
     2760
     2761  if (startIndices.size() != endIndices.size()) {
     2762    throw(AipsError("Inconsistent Mask Size: bad data?"));
     2763  }
     2764  for (uInt i = 0; i < startIndices.size(); ++i) {
     2765    if (startIndices[i] > endIndices[i]) {
     2766      throw(AipsError("Mask start index > mask end index"));
     2767    }
     2768  }
     2769
     2770  out.clear();
     2771  for (uInt i = 0; i < startIndices.size(); ++i) {
     2772    out.push_back(startIndices[i]);
     2773    out.push_back(endIndices[i]);
    25032774  }
    25042775

  • trunk/src/Scantable.h

    r2012 r2047  
    493493  void regridChannel( int nchan, double dnu, int irow ) ;
    494494
    495   bool getFlagtraFast(int whichrow);
     495  bool getFlagtraFast(casa::uInt whichrow);
    496496
    497497  void polyBaseline(const std::vector<bool>& mask,
    … …  
    521521                               bool outLogger=false,
    522522                               const std::string& blfile="");
     523  void sinusoidBaseline(const std::vector<bool>& mask,
     524                        int nMinWavesInSW,
     525                        int nMaxWavesInSW,
     526                        float thresClip,
     527                        int nIterClip,
     528                        bool outLogger=false,
     529                        const std::string& blfile="");
     530  void autoSinusoidBaseline(const std::vector<bool>& mask,
     531                            int nMinWavesInSW,
     532                            int nMaxWavesInSW,
     533                            float thresClip,
     534                            int nIterClip,
     535                            const std::vector<int>& edge,
     536                            float threshold=3.0,
     537                            int chanAvgLimit=1,
     538                            bool outLogger=false,
     539                            const std::string& blfile="");
    523540  float getRms(const std::vector<bool>& mask, int whichrow);
    524541  std::string formatBaselineParams(const std::vector<float>& params,
    … …  
    657674                                          float thresClip=3.0,
    658675                                          int nIterClip=1);
     676  std::vector<float> doSinusoidFitting(const std::vector<float>& data,
     677                                       const std::vector<bool>& mask,
     678                                       int nMinWavesInSW,
     679                                       int nMaxWavesInSW,
     680                                       std::vector<float>& params,
     681                                       float thresClip=3.0,
     682                                       int nIterClip=1);
    659683  bool hasSameNchanOverIFs();
    660684  std::string getMaskRangeList(const std::vector<bool>& mask,
    661685                                int whichrow,
    662686                                const casa::String& coordInfo,
    663                                 bool hasSameNchan,
    664                                 int firstIF,
    665                                 bool silent=false);
    666   std::vector<int> getMaskEdgeIndices(const std::vector<bool>& mask, bool getStartIndices=true);
     687                                bool hasSameNchan,
     688                                bool verbose=false);
     689  std::vector<int> getMaskEdgeIndices(const std::vector<bool>& mask);
    667690  std::string formatBaselineParamsHeader(int whichrow, const std::string& masklist, bool verbose) const;
    668691  std::string formatBaselineParamsFooter(float rms, bool verbose) const;
    669692  std::vector<bool> getCompositeChanMask(int whichrow, const std::vector<bool>& inMask);
    670693  //std::vector<bool> getCompositeChanMask(int whichrow, const std::vector<bool>& inMask, const std::vector<int>& edge, const int minEdgeSize, STLineFinder& lineFinder);
     694  void outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, std::ofstream& ofs, const casa::String& funcName, Fitter& fitter);
     695  void outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, std::ofstream& ofs, const casa::String& funcName, const std::vector<int>& pieceEdges, const std::vector<float>& params, const std::vector<bool>& fixed);
     696  void outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, std::ofstream& ofs, const casa::String& funcName, const std::vector<float>& params, const std::vector<bool>& fixed);
    671697
    672698  void applyChanFlag( casa::uInt whichrow, const std::vector<bool>& msk, casa::uChar flagval);

  • trunk/src/ScantableWrapper.h

    r2012 r2047  
    269269  { table_->autoCubicSplineBaseline(mask, npiece, clipthresh, clipniter, edge, threshold, chan_avg_limit, outlog, blfile); }
    270270
     271  void sinusoidBaseline(const std::vector<bool>& mask, int minnwave, int maxnwave, float clipthresh, int clipniter, bool outlog=false, const std::string& blfile="")
     272  { table_->sinusoidBaseline(mask, minnwave, maxnwave, clipthresh, clipniter, outlog, blfile); }
     273
     274  void autoSinusoidBaseline(const std::vector<bool>& mask, int minnwave, int maxnwave, float clipthresh, int clipniter, const std::vector<int>& edge, float threshold=5.0, int chan_avg_limit=1, bool outlog=false, const std::string& blfile="")
     275  { table_->autoSinusoidBaseline(mask, minnwave, maxnwave, clipthresh, clipniter, edge, threshold, chan_avg_limit, outlog, blfile); }
     276
    271277  float getRms(const std::vector<bool>& mask, int whichrow)
    272278  { return table_->getRms(mask, whichrow); }
    … …  
    279285
    280286  bool getFlagtraFast(int whichrow=0) const
    281   { return table_->getFlagtraFast(whichrow); }
     287  { return table_->getFlagtraFast(casa::uInt(whichrow)); }
    282288
    283289

  • trunk/src/python_Scantable.cpp

    r2012 r2047  
    146146    .def("_cspline_baseline", &ScantableWrapper::cubicSplineBaseline)
    147147    .def("_auto_cspline_baseline", &ScantableWrapper::autoCubicSplineBaseline)
     148    .def("_sinusoid_baseline", &ScantableWrapper::sinusoidBaseline)
     149    .def("_auto_sinusoid_baseline", &ScantableWrapper::autoSinusoidBaseline)
    148150    .def("get_rms", &ScantableWrapper::getRms)
    149151    .def("format_blparams_row", &ScantableWrapper::formatBaselineParams)
Note: See TracChangeset for help on using the changeset viewer.