Changeset 2012

Timestamp:

02/25/11 16:51:50 (13 years ago)

Author:

WataruKawasaki

Message:

New Development: Yes

JIRA Issue: Yes (CAS-2373, CAS-2620)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: For Scantable::polyBaseline(), parameters and return value have been changed.

Test Programs:

Put in Release Notes: Yes

Module(s): sdbaseline, sd.linefinder

Description: (1) CAS-2373-related:

(1.1) Modified Scantable::polyBaseline() to have the row-based loop inside.

Now it fits and subtracts baseline for all rows and also output info
about the fitting result to logger and text file, while in the
previous version this method just did baseline fitting/subtraction
for one row only and had to be put inside a row-based loop at the
python side ("poly_baseline()" in scantable.py) and result output had
also to be controlled at the python side. Using a test data from NRO
45m telescope (348,000 rows, 512 channels), the processing time of
scantable.poly_baseline() has reduced from 130 minutes to 5-6 minutes.

(1.2) For accelerating the another polynomial fitting function, namely

scantable.auto_poly_baseline(), added a method
Scantable::autoPolyBaseline(). This basically does the same thing
with Scantable::polyBaseline(), but uses linefinfer also to
automatically flag the line regions.

(1.3) To make linefinder usable in Scantable class, added a method

linefinder.setdata(). This makes it possible to apply linefinder
for a float-list data given as a parameter, without setting scantable,
while it was indispensable to set scantable to use linefinger previously.

(2) CAS-2620-related:

Added Scantable::cubicSplineBaseline() and autoCubicSplineBaseline() for
fit baseline using the cubic spline function. Parameters include npiece
(number of polynomial pieces), clipthresh (clipping threshold), and
clipniter (maximum iteration number).

Location:

trunk

Files:

• python/asaplinefind.py (modified) (1 diff)
• python/scantable.py (modified) (3 diffs)
• src/STLineFinder.cpp (modified) (8 diffs)
• src/STLineFinder.h (modified) (2 diffs)
• src/Scantable.cpp (modified) (3 diffs)
• src/Scantable.h (modified) (3 diffs)
• src/ScantableWrapper.h (modified) (2 diffs)
• src/python_STLineFinder.cpp (modified) (1 diff)
• src/python_Scantable.cpp (modified) (1 diff)

trunk/python/asaplinefind.py

@@ -101 +101 @@
         self.finder.setscan(scan)
 
+    def set_data(self, spectrum):
+        """
+        Set the 'data' (spectrum) to work with
+        Parameters: a method to allow linefinder work without setting scantable
+        for the purpose of using linefinder inside some method in scantable
+        class. (Dec 22, 2010 by W.Kawasaki)
+        """
+        if isinstance(spectrum, list) or isinstance(spectrum, tuple):
+            if not isinstance(spectrum[0], float):
+                raise RuntimeError, "Parameter 'spectrum' has to be a list or a tuple of float"
+        else:
+            raise RuntimeError, "Parameter 'spectrum' has to be a list or a tuple of float"
+        self.finder.setdata(spectrum)
+        
     def find_lines(self,nRow=0,mask=[],edge=(0,0)):
         """

trunk/python/scantable.py

@@ -1898 +1898 @@
         else: return s
 
+
+    @asaplog_post_dec
+    def cspline_baseline(self, insitu=None, mask=None, npiece=None, clipthresh=None, clipniter=None, plot=None, outlog=None, blfile=None):
+        """\
+        Return a scan which has been baselined (all rows) by cubic spline function (piecewise cubic polynomial).
+        Parameters:
+            insitu:     If False a new scantable is returned.
+                        Otherwise, the scaling is done in-situ
+                        The default is taken from .asaprc (False)
+            mask:       An optional mask
+            npiece:     Number of pieces. (default is 2)
+            clipthresh: Clipping threshold. (default is 3.0, unit: sigma)
+            clipniter:  maximum number of iteration of 'clipthresh'-sigma clipping (default is 1)
+            plot:   *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
+                        plot the fit and the residual. In this each
+                        indivual fit has to be approved, by typing 'y'
+                        or 'n'
+            outlog:     Output the coefficients of the best-fit
+                        function to logger (default is False)
+            blfile:     Name of a text file in which the best-fit
+                        parameter values to be written
+                        (default is "": no file/logger output)
+
+        Example:
+            # return a scan baselined by a cubic spline consisting of 2 pieces (i.e., 1 internal knot),
+            # also with 3-sigma clipping, iteration up to 4 times
+            bscan = scan.cspline_baseline(npiece=2,clipthresh=3.0,clipniter=4)
+        """
+        
+        varlist = vars()
+        
+        if insitu is None: insitu = rcParams["insitu"]
+        if insitu:
+            workscan = self
+        else:
+            workscan = self.copy()
+
+        nchan = workscan.nchan()
+        
+        if mask is None: mask = [True for i in xrange(nchan)]
+        if npiece is None: npiece = 2
+        if clipthresh is None: clipthresh = 3.0
+        if clipniter is None: clipniter = 1
+        if plot is None: plot = False
+        if outlog is None: outlog = False
+        if blfile is None: blfile = ""
+
+        outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
+        
+        try:
+            #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method. 
+            workscan._cspline_baseline(mask, npiece, clipthresh, clipniter, outlog, blfile)
+            
+            workscan._add_history("cspline_baseline", varlist)
+            
+            if insitu:
+                self._assign(workscan)
+            else:
+                return workscan
+            
+        except RuntimeError, e:
+            msg = "The fit failed, possibly because it didn't converge."
+            if rcParams["verbose"]:
+                asaplog.push(str(e))
+                asaplog.push(str(msg))
+                return
+            else:
+                raise RuntimeError(str(e)+'\n'+msg)
+
+
+    def auto_cspline_baseline(self, insitu=None, mask=None, npiece=None, clipthresh=None,
+                              clipniter=None, edge=None, threshold=None,
+                              chan_avg_limit=None, plot=None, outlog=None, blfile=None):
+        """\
+        Return a scan which has been baselined (all rows) by cubic spline
+        function (piecewise cubic polynomial).
+        Spectral lines are detected first using linefinder and masked out
+        to avoid them affecting the baseline solution.
+
+        Parameters:
+            insitu:     if False a new scantable is returned.
+                        Otherwise, the scaling is done in-situ
+                        The default is taken from .asaprc (False)
+            mask:       an optional mask retreived from scantable
+            npiece:     Number of pieces. (default is 2)
+            clipthresh: Clipping threshold. (default is 3.0, unit: sigma)
+            clipniter:  maximum number of iteration of 'clipthresh'-sigma clipping (default is 1)
+            edge:       an optional number of channel to drop at
+                        the edge of spectrum. If only one value is
+                        specified, the same number will be dropped
+                        from both sides of the spectrum. Default
+                        is to keep all channels. Nested tuples
+                        represent individual edge selection for
+                        different IFs (a number of spectral channels
+                        can be different)
+            threshold:  the threshold used by line finder. It is
+                        better to keep it large as only strong lines
+                        affect the baseline solution.
+            chan_avg_limit:
+                        a maximum number of consequtive spectral
+                        channels to average during the search of
+                        weak and broad lines. The default is no
+                        averaging (and no search for weak lines).
+                        If such lines can affect the fitted baseline
+                        (e.g. a high order polynomial is fitted),
+                        increase this parameter (usually values up
+                        to 8 are reasonable). Most users of this
+                        method should find the default value sufficient.
+            plot:   *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
+                        plot the fit and the residual. In this each
+                        indivual fit has to be approved, by typing 'y'
+                        or 'n'
+            outlog:     Output the coefficients of the best-fit
+                        function to logger (default is False)
+            blfile:     Name of a text file in which the best-fit
+                        parameter values to be written
+                        (default is "": no file/logger output)
+
+        Example:
+            bscan = scan.auto_cspline_baseline(npiece=3, insitu=False)
+        """
+
+        varlist = vars()
+
+        if insitu is None: insitu = rcParams['insitu']
+        if insitu:
+            workscan = self
+        else:
+            workscan = self.copy()
+
+        nchan = workscan.nchan()
+        
+        if mask is None: mask = [True for i in xrange(nchan)]
+        if npiece is None: npiece = 2
+        if clipthresh is None: clipthresh = 3.0
+        if clipniter is None: clipniter = 1
+        if edge is None: edge = (0, 0)
+        if threshold is None: threshold = 3
+        if chan_avg_limit is None: chan_avg_limit = 1
+        if plot is None: plot = False
+        if outlog is None: outlog = False
+        if blfile is None: blfile = ""
+
+        outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
+        
+        from asap.asaplinefind import linefinder
+        from asap import _is_sequence_or_number as _is_valid
+
+        if not (isinstance(edge, list) or isinstance(edge, tuple)): edge = [ edge ]
+        individualedge = False;
+        if len(edge) > 1: individualedge = isinstance(edge[0], list) or isinstance(edge[0], tuple)
+
+        if individualedge:
+            for edgepar in edge:
+                if not _is_valid(edgepar, int):
+                    raise ValueError, "Each element of the 'edge' tuple has \
+                                       to be a pair of integers or an integer."
+        else:
+            if not _is_valid(edge, int):
+                raise ValueError, "Parameter 'edge' has to be an integer or a \
+                                   pair of integers specified as a tuple. \
+                                   Nested tuples are allowed \
+                                   to make individual selection for different IFs."
+
+            if len(edge) > 1:
+                curedge = edge
+            else:
+                curedge = edge + edge
+
+        try:
+            #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method. 
+            if individualedge:
+                curedge = []
+                for i in xrange(len(edge)):
+                    curedge += edge[i]
+                
+            workscan._auto_cspline_baseline(mask, npiece, clipthresh, clipniter, curedge, threshold, chan_avg_limit, outlog, blfile)
+
+            workscan._add_history("auto_cspline_baseline", varlist)
+            
+            if insitu:
+                self._assign(workscan)
+            else:
+                return workscan
+            
+        except RuntimeError, e:
+            msg = "The fit failed, possibly because it didn't converge."
+            if rcParams["verbose"]:
+                asaplog.push(str(e))
+                asaplog.push(str(msg))
+                return
+            else:
+                raise RuntimeError(str(e)+'\n'+msg)
+
+
+    @asaplog_post_dec
+    def poly_baseline(self, insitu=None, mask=None, order=None, plot=None, outlog=None, blfile=None):
+        """\
+        Return a scan which has been baselined (all rows) by a polynomial.
+        Parameters:
+            insitu:     if False a new scantable is returned.
+                        Otherwise, the scaling is done in-situ
+                        The default is taken from .asaprc (False)
+            mask:       an optional mask
+            order:      the order of the polynomial (default is 0)
+            plot:       plot the fit and the residual. In this each
+                        indivual fit has to be approved, by typing 'y'
+                        or 'n'
+            outlog:     Output the coefficients of the best-fit
+                        function to logger (default is False)
+            blfile:     Name of a text file in which the best-fit
+                        parameter values to be written
+                        (default is "": no file/logger output)
+
+        Example:
+            # return a scan baselined by a third order polynomial,
+            # not using a mask
+            bscan = scan.poly_baseline(order=3)
+        """
+        
+        varlist = vars()
+        
+        if insitu is None: insitu = rcParams["insitu"]
+        if insitu:
+            workscan = self
+        else:
+            workscan = self.copy()
+
+        nchan = workscan.nchan()
+        
+        if mask is None: mask = [True for i in xrange(nchan)]
+        if order is None: order = 0
+        if plot is None: plot = False
+        if outlog is None: outlog = False
+        if blfile is None: blfile = ""
+
+        outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
+        
+        try:
+            rows = xrange(workscan.nrow())
+            
+            #if len(rows) > 0: workscan._init_blinfo()
+
+            if plot:
+                if outblfile: blf = open(blfile, "a")
+                
+                f = fitter()
+                f.set_function(lpoly=order)
+                for r in rows:
+                    f.x = workscan._getabcissa(r)
+                    f.y = workscan._getspectrum(r)
+                    f.mask = mask_and(mask, workscan._getmask(r))    # (CAS-1434)
+                    f.data = None
+                    f.fit()
+                    
+                    f.plot(residual=True)
+                    accept_fit = raw_input("Accept fit ( [y]/n ): ")
+                    if accept_fit.upper() == "N":
+                        #workscan._append_blinfo(None, None, None)
+                        continue
+                    
+                    blpars = f.get_parameters()
+                    masklist = workscan.get_masklist(f.mask, row=r, silent=True)
+                    #workscan._append_blinfo(blpars, masklist, f.mask)
+                    workscan._setspectrum(f.fitter.getresidual(), r)
+                    
+                    if outblfile:
+                        rms = workscan.get_rms(f.mask, r)
+                        dataout = workscan.format_blparams_row(blpars["params"], blpars["fixed"], rms, str(masklist), r, True)
+                        blf.write(dataout)
+
+                f._p.unmap()
+                f._p = None
+
+                if outblfile: blf.close()
+            else:
+                workscan._poly_baseline(mask, order, outlog, blfile)
+            
+            workscan._add_history("poly_baseline", varlist)
+            
+            if insitu:
+                self._assign(workscan)
+            else:
+                return workscan
+            
+        except RuntimeError, e:
+            msg = "The fit failed, possibly because it didn't converge."
+            if rcParams["verbose"]:
+                asaplog.push(str(e))
+                asaplog.push(str(msg))
+                return
+            else:
+                raise RuntimeError(str(e)+'\n'+msg)
+
+
+    def auto_poly_baseline(self, insitu=None, mask=None, order=None, edge=None, threshold=None,
+                           chan_avg_limit=None, plot=None, outlog=None, blfile=None):
+        """\
+        Return a scan which has been baselined (all rows) by a polynomial.
+        Spectral lines are detected first using linefinder and masked out
+        to avoid them affecting the baseline solution.
+
+        Parameters:
+            insitu:     if False a new scantable is returned.
+                        Otherwise, the scaling is done in-situ
+                        The default is taken from .asaprc (False)
+            mask:       an optional mask retreived from scantable
+            order:      the order of the polynomial (default is 0)
+            edge:       an optional number of channel to drop at
+                        the edge of spectrum. If only one value is
+                        specified, the same number will be dropped
+                        from both sides of the spectrum. Default
+                        is to keep all channels. Nested tuples
+                        represent individual edge selection for
+                        different IFs (a number of spectral channels
+                        can be different)
+            threshold:  the threshold used by line finder. It is
+                        better to keep it large as only strong lines
+                        affect the baseline solution.
+            chan_avg_limit:
+                        a maximum number of consequtive spectral
+                        channels to average during the search of
+                        weak and broad lines. The default is no
+                        averaging (and no search for weak lines).
+                        If such lines can affect the fitted baseline
+                        (e.g. a high order polynomial is fitted),
+                        increase this parameter (usually values up
+                        to 8 are reasonable). Most users of this
+                        method should find the default value sufficient.
+            plot:       plot the fit and the residual. In this each
+                        indivual fit has to be approved, by typing 'y'
+                        or 'n'
+            outlog:     Output the coefficients of the best-fit
+                        function to logger (default is False)
+            blfile:     Name of a text file in which the best-fit
+                        parameter values to be written
+                        (default is "": no file/logger output)
+
+        Example:
+            bscan = scan.auto_poly_baseline(order=7, insitu=False)
+        """
+
+        varlist = vars()
+
+        if insitu is None: insitu = rcParams['insitu']
+        if insitu:
+            workscan = self
+        else:
+            workscan = self.copy()
+
+        nchan = workscan.nchan()
+        
+        if mask is None: mask = [True for i in xrange(nchan)]
+        if order is None: order = 0
+        if edge is None: edge = (0, 0)
+        if threshold is None: threshold = 3
+        if chan_avg_limit is None: chan_avg_limit = 1
+        if plot is None: plot = False
+        if outlog is None: outlog = False
+        if blfile is None: blfile = ""
+
+        outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
+        
+        from asap.asaplinefind import linefinder
+        from asap import _is_sequence_or_number as _is_valid
+
+        if not (isinstance(edge, list) or isinstance(edge, tuple)): edge = [ edge ]
+        individualedge = False;
+        if len(edge) > 1: individualedge = isinstance(edge[0], list) or isinstance(edge[0], tuple)
+
+        if individualedge:
+            for edgepar in edge:
+                if not _is_valid(edgepar, int):
+                    raise ValueError, "Each element of the 'edge' tuple has \
+                                       to be a pair of integers or an integer."
+        else:
+            if not _is_valid(edge, int):
+                raise ValueError, "Parameter 'edge' has to be an integer or a \
+                                   pair of integers specified as a tuple. \
+                                   Nested tuples are allowed \
+                                   to make individual selection for different IFs."
+
+            if len(edge) > 1:
+                curedge = edge
+            else:
+                curedge = edge + edge
+
+        try:
+            rows = xrange(workscan.nrow())
+            
+            #if len(rows) > 0: workscan._init_blinfo()
+
+            if plot:
+                if outblfile: blf = open(blfile, "a")
+                
+                fl = linefinder()
+                fl.set_options(threshold=threshold,avg_limit=chan_avg_limit)
+                fl.set_scan(workscan)
+                f = fitter()
+                f.set_function(lpoly=order)
+
+                for r in rows:
+                    if individualedge:
+                        if len(edge) <= workscan.getif(r):
+                            raise RuntimeError, "Number of edge elements appear to " \
+                                  "be less than the number of IFs"
+                        else:
+                            curedge = edge[workscan.getif(r)]
+
+                    fl.find_lines(r, mask_and(mask, workscan._getmask(r)), curedge)  # (CAS-1434)
+
+                    f.x = workscan._getabcissa(r)
+                    f.y = workscan._getspectrum(r)
+                    f.mask = fl.get_mask()
+                    f.data = None
+                    f.fit()
+
+                    f.plot(residual=True)
+                    accept_fit = raw_input("Accept fit ( [y]/n ): ")
+                    if accept_fit.upper() == "N":
+                        #workscan._append_blinfo(None, None, None)
+                        continue
+
+                    blpars = f.get_parameters()
+                    masklist = workscan.get_masklist(f.mask, row=r, silent=True)
+                    #workscan._append_blinfo(blpars, masklist, f.mask)
+                    workscan._setspectrum(f.fitter.getresidual(), r)
+
+                    if outblfile:
+                        rms = workscan.get_rms(f.mask, r)
+                        dataout = workscan.format_blparams_row(blpars["params"], blpars["fixed"], rms, str(masklist), r, True)
+                        blf.write(dataout)
+                    
+                f._p.unmap()
+                f._p = None
+
+                if outblfile: blf.close()
+                
+            else:
+                if individualedge:
+                    curedge = []
+                    for i in xrange(len(edge)):
+                        curedge += edge[i]
+                
+                workscan._auto_poly_baseline(mask, order, curedge, threshold, chan_avg_limit, outlog, blfile)
+
+            workscan._add_history("auto_poly_baseline", varlist)
+            
+            if insitu:
+                self._assign(workscan)
+            else:
+                return workscan
+            
+        except RuntimeError, e:
+            msg = "The fit failed, possibly because it didn't converge."
+            if rcParams["verbose"]:
+                asaplog.push(str(e))
+                asaplog.push(str(msg))
+                return
+            else:
+                raise RuntimeError(str(e)+'\n'+msg)
+
+
+    ### OBSOLETE ##################################################################
     @asaplog_post_dec
     def old_poly_baseline(self, mask=None, order=0, plot=False, uselin=False, insitu=None, rows=None):
-        """\
+        """
         Return a scan which has been baselined (all rows) by a polynomial.
         
@@ -1985 +2449 @@
             raise RuntimeError(msg)
 
-    @asaplog_post_dec
-    def poly_baseline(self, mask=None, order=0, plot=False, batch=False, insitu=None, rows=None):
-        """\
-        Return a scan which has been baselined (all rows) by a polynomial.
-        Parameters:
-            mask:       an optional mask
-            order:      the order of the polynomial (default is 0)
-            plot:       plot the fit and the residual. In this each
-                        indivual fit has to be approved, by typing 'y'
-                        or 'n'. Ignored if batch = True. 
-            batch:      if True a faster algorithm is used and logs
-                        including the fit results are not output
-                        (default is False) 
-            insitu:     if False a new scantable is returned.
-                        Otherwise, the scaling is done in-situ
-                        The default is taken from .asaprc (False)
-            rows:       row numbers of spectra to be baselined.
-                        (default is None: for all rows)
-        Example:
-            # return a scan baselined by a third order polynomial,
-            # not using a mask
-            bscan = scan.poly_baseline(order=3)
-        """
+    def _init_blinfo(self):
+        """\
+        Initialise the following three auxiliary members:
+           blpars : parameters of the best-fit baseline, 
+           masklists : mask data (edge positions of masked channels) and 
+           actualmask : mask data (in boolean list), 
+        to keep for use later (including output to logger/text files). 
+        Used by poly_baseline() and auto_poly_baseline() in case of
+        'plot=True'. 
+        """
+        self.blpars = []
+        self.masklists = []
+        self.actualmask = []
+        return
+
+    def _append_blinfo(self, data_blpars, data_masklists, data_actualmask):
+        """\
+        Append baseline-fitting related info to blpars, masklist and
+        actualmask. 
+        """
+        self.blpars.append(data_blpars)
+        self.masklists.append(data_masklists)
+        self.actualmask.append(data_actualmask)
+        return
         
-        varlist = vars()
-        
-        if insitu is None: insitu = rcParams["insitu"]
-        if insitu:
-            workscan = self
-        else:
-            workscan = self.copy()
-
-        nchan = workscan.nchan()
-        
-        if mask is None:
-            mask = [True for i in xrange(nchan)]
-
-        try:
-            if rows == None:
-                rows = xrange(workscan.nrow())
-            elif isinstance(rows, int):
-                rows = [ rows ]
-            
-            if len(rows) > 0:
-                workscan.blpars = []
-                workscan.masklists = []
-                workscan.actualmask = []
-
-            if batch:
-                workscan._poly_baseline_batch(mask, order)
-            elif plot:
-                f = fitter()
-                f.set_function(lpoly=order)
-                for r in rows:
-                    f.x = workscan._getabcissa(r)
-                    f.y = workscan._getspectrum(r)
-                    f.mask = mask_and(mask, workscan._getmask(r))    # (CAS-1434)
-                    f.data = None
-                    f.fit()
-                    
-                    f.plot(residual=True)
-                    accept_fit = raw_input("Accept fit ( [y]/n ): ")
-                    if accept_fit.upper() == "N":
-                        self.blpars.append(None)
-                        self.masklists.append(None)
-                        self.actualmask.append(None)
-                        continue
-                    workscan._setspectrum(f.fitter.getresidual(), r)
-                    workscan.blpars.append(f.get_parameters())
-                    workscan.masklists.append(workscan.get_masklist(f.mask, row=r))
-                    workscan.actualmask.append(f.mask)
-                    
-                f._p.unmap()
-                f._p = None
-            else:
-                for r in rows:
-                    fitparams = workscan._poly_baseline(mask, order, r)
-                    params = fitparams.getparameters()
-                    fmtd = ", ".join(["p%d = %3.6f" % (i, v) for i, v in enumerate(params)])
-                    errors = fitparams.geterrors()
-                    fmask = mask_and(mask, workscan._getmask(r))
-
-                    workscan.blpars.append({"params":params,
-                                            "fixed": fitparams.getfixedparameters(),
-                                            "formatted":fmtd, "errors":errors})
-                    workscan.masklists.append(workscan.get_masklist(fmask, r, silent=True))
-                    workscan.actualmask.append(fmask)
-                    
-                    asaplog.push(fmtd)
-            
-            workscan._add_history("poly_baseline", varlist)
-            
-            if insitu:
-                self._assign(workscan)
-            else:
-                return workscan
-            
-        except RuntimeError, e:
-            msg = "The fit failed, possibly because it didn't converge."
-            if rcParams["verbose"]:
-                asaplog.push(str(e))
-                asaplog.push(str(msg))
-                return
-            else:
-                raise RuntimeError(str(e)+'\n'+msg)
-
-
-    def auto_poly_baseline(self, mask=None, edge=(0, 0), order=0,
-                           threshold=3, chan_avg_limit=1, plot=False,
-                           insitu=None, rows=None):
-        """\
-        Return a scan which has been baselined (all rows) by a polynomial.
-        Spectral lines are detected first using linefinder and masked out
-        to avoid them affecting the baseline solution.
-
-        Parameters:
-
-            mask:       an optional mask retreived from scantable
-
-            edge:       an optional number of channel to drop at the edge of
-                        spectrum. If only one value is
-                        specified, the same number will be dropped from
-                        both sides of the spectrum. Default is to keep
-                        all channels. Nested tuples represent individual
-                        edge selection for different IFs (a number of spectral
-                        channels can be different)
-
-            order:      the order of the polynomial (default is 0)
-
-            threshold:  the threshold used by line finder. It is better to
-                        keep it large as only strong lines affect the
-                        baseline solution.
-
-            chan_avg_limit:
-                        a maximum number of consequtive spectral channels to
-                        average during the search of weak and broad lines.
-                        The default is no averaging (and no search for weak
-                        lines). If such lines can affect the fitted baseline
-                        (e.g. a high order polynomial is fitted), increase this
-                        parameter (usually values up to 8 are reasonable). Most
-                        users of this method should find the default value
-                        sufficient.
-
-            plot:       plot the fit and the residual. In this each
-                        indivual fit has to be approved, by typing 'y'
-                        or 'n'
-
-            insitu:     if False a new scantable is returned.
-                        Otherwise, the scaling is done in-situ
-                        The default is taken from .asaprc (False)
-            rows:       row numbers of spectra to be processed.
-                        (default is None: for all rows)
-
-
-        Example::
-
-            scan2 = scan.auto_poly_baseline(order=7, insitu=False)
-
-        """
-        if insitu is None: insitu = rcParams['insitu']
-        varlist = vars()
-        from asap.asaplinefind import linefinder
-        from asap import _is_sequence_or_number as _is_valid
-
-        # check whether edge is set up for each IF individually
-        individualedge = False;
-        if len(edge) > 1:
-            if isinstance(edge[0], list) or isinstance(edge[0], tuple):
-                individualedge = True;
-
-        if not _is_valid(edge, int) and not individualedge:
-            raise ValueError, "Parameter 'edge' has to be an integer or a \
-            pair of integers specified as a tuple. Nested tuples are allowed \
-            to make individual selection for different IFs."
-
-        curedge = (0, 0)
-        if individualedge:
-            for edgepar in edge:
-                if not _is_valid(edgepar, int):
-                    raise ValueError, "Each element of the 'edge' tuple has \
-                                       to be a pair of integers or an integer."
-        else:
-            curedge = edge;
-
-        if not insitu:
-            workscan = self.copy()
-        else:
-            workscan = self
-
-        # setup fitter
-        f = fitter()
-        f.set_function(lpoly=order)
-
-        # setup line finder
-        fl = linefinder()
-        fl.set_options(threshold=threshold,avg_limit=chan_avg_limit)
-
-        fl.set_scan(workscan)
-
-        if mask is None:
-            mask = _n_bools(workscan.nchan(), True)
-        
-        if rows is None:
-            rows = xrange(workscan.nrow())
-        elif isinstance(rows, int):
-            rows = [ rows ]
-        
-        # Save parameters of baseline fits & masklists as a class attribute.
-        # NOTICE: It does not reflect changes in scantable!
-        if len(rows) > 0:
-            self.blpars=[]
-            self.masklists=[]
-            self.actualmask=[]
-        asaplog.push("Processing:")
-        for r in rows:
-            msg = " Scan[%d] Beam[%d] IF[%d] Pol[%d] Cycle[%d]" % \
-                (workscan.getscan(r), workscan.getbeam(r), workscan.getif(r), \
-                 workscan.getpol(r), workscan.getcycle(r))
-            asaplog.push(msg, False)
-
-            # figure out edge parameter
-            if individualedge:
-                if len(edge) >= workscan.getif(r):
-                    raise RuntimeError, "Number of edge elements appear to " \
-                                        "be less than the number of IFs"
-                    curedge = edge[workscan.getif(r)]
-
-            actualmask = mask_and(mask, workscan._getmask(r))    # (CAS-1434)
-
-            # setup line finder
-            fl.find_lines(r, actualmask, curedge)
-            
-            f.x = workscan._getabcissa(r)
-            f.y = workscan._getspectrum(r)
-            f.mask = fl.get_mask()
-            f.data = None
-            f.fit()
-
-            # Show mask list
-            masklist=workscan.get_masklist(f.mask, row=r, silent=True)
-            msg = "mask range: "+str(masklist)
-            asaplog.push(msg, False)
-
-            if plot:
-                f.plot(residual=True)
-                x = raw_input("Accept fit ( [y]/n ): ")
-                if x.upper() == 'N':
-                    self.blpars.append(None)
-                    self.masklists.append(None)
-                    self.actualmask.append(None)
-                    continue
-
-            workscan._setspectrum(f.fitter.getresidual(), r)
-            self.blpars.append(f.get_parameters())
-            self.masklists.append(masklist)
-            self.actualmask.append(f.mask)
-        if plot:
-            f._p.unmap()
-            f._p = None
-        workscan._add_history("auto_poly_baseline", varlist)
-        if insitu:
-            self._assign(workscan)
-        else:
-            return workscan
-
     @asaplog_post_dec
     def rotate_linpolphase(self, angle):
@@ -2744 +2968 @@
             self.set_freqframe(rcParams['scantable.freqframe'])
 
+
     def __getitem__(self, key):
         if key < 0:

trunk/src/STLineFinder.cpp

@@ -925 +925 @@
   itsNoiseBox = in_noise_box;
   itsUseMedian = in_median;
+
+  useScantable = true;
 }
 
@@ -934 +936 @@
   scan=in_scan.getCP();
   AlwaysAssert(!scan.null(),AipsError);
-
+}
+
+// set spectrum data to work with. this is a method to allow linefinder work 
+// without setting scantable for the purpose of using linefinder inside some 
+// method in scantable class. (Dec 22, 2010 by W.Kawasaki)
+void STLineFinder::setData(const std::vector<float> &in_spectrum)
+{
+  spectrum = Vector<Float>(in_spectrum);
+  useScantable = false;
 }
 
@@ -948 +958 @@
                 const casa::uInt &whichRow) throw(casa::AipsError)
 {
-  if (scan.null())
+  if (useScantable && scan.null())
       throw AipsError("STLineFinder::findLines - a scan should be set first,"
                       " use set_scan");
 
-  uInt nchan = scan->nchan(scan->getIF(whichRow));
+  uInt nchan = useScantable ? scan->nchan(scan->getIF(whichRow)) : spectrum.nelements();
   // set up mask and edge rejection
   // no mask given...
@@ -962 +972 @@
   }
   if (mask.nelements()!=nchan)
-      throw AipsError("STLineFinder::findLines - in_scan and in_mask have different"
-            "number of spectral channels.");
+      throw AipsError("STLineFinder::findLines - in_scan and in_mask, or in_spectrum "
+                      "and in_mask have different number of spectral channels.");
 
   // taking flagged channels into account
-  vector<bool> flaggedChannels = scan->getMask(whichRow);
-  if (flaggedChannels.size()) {
+  if (useScantable) {
+    vector<bool> flaggedChannels = scan->getMask(whichRow);
+    if (flaggedChannels.size()) {
       // there is a mask set for this row
       if (flaggedChannels.size() != mask.nelements()) {
-          throw AipsError("STLineFinder::findLines - internal inconsistency: number of mask elements do not match the number of channels");
+          throw AipsError("STLineFinder::findLines - internal inconsistency: number of "
+                          "mask elements do not match the number of channels");
       }
       for (size_t ch = 0; ch<mask.nelements(); ++ch) {
            mask[ch] &= flaggedChannels[ch];
       }
+    }
   }
 
@@ -1015 +1028 @@
       throw AipsError("STLineFinder::findLines - noise_box is supposed to be at least 2 elements");
 
-  spectrum.resize();
-  spectrum = Vector<Float>(scan->getSpectrum(whichRow));
+  if (useScantable) {
+    spectrum.resize();
+    spectrum = Vector<Float>(scan->getSpectrum(whichRow));
+  }
 
   lines.resize(0); // search from the scratch
@@ -1141 +1156 @@
 {
   try {
-       if (scan.null())
-           throw AipsError("STLineFinder::getMask - a scan should be set first,"
-                      " use set_scan followed by find_lines");
-       DebugAssert(mask.nelements()==scan->getChannels(last_row_used), AipsError);
-       /*
-       if (!lines.size())
-           throw AipsError("STLineFinder::getMask - one have to search for "
+    if (useScantable) {
+      if (scan.null())
+        throw AipsError("STLineFinder::getMask - a scan should be set first,"
+                        " use set_scan followed by find_lines");
+      DebugAssert(mask.nelements()==scan->getChannels(last_row_used), AipsError);
+    }
+    /*
+    if (!lines.size())
+       throw AipsError("STLineFinder::getMask - one have to search for "
                            "lines first, use find_lines");
-       */
-       std::vector<bool> res_mask(mask.nelements());
-       // iterator through lines
-       std::list<std::pair<int,int> >::const_iterator cli=lines.begin();
-       for (int ch=0;ch<int(res_mask.size());++ch) {
-            if (ch<edge.first || ch>=edge.second) res_mask[ch]=false;
-            else if (!mask[ch]) res_mask[ch]=false;
-            else {
-                    res_mask[ch]=!invert; // no line by default
-                    if (cli!=lines.end())
-                        if (ch>=cli->first && ch<cli->second)
-                             res_mask[ch]=invert; // this is a line
-            }
-            if (cli!=lines.end())
-                if (ch>=cli->second) {
-                    ++cli; // next line in the list
-                }
-       }
-       return res_mask;
+    */
+    std::vector<bool> res_mask(mask.nelements());
+    // iterator through lines
+    std::list<std::pair<int,int> >::const_iterator cli=lines.begin();
+    for (int ch=0;ch<int(res_mask.size());++ch) {
+      if (ch<edge.first || ch>=edge.second) res_mask[ch]=false;
+      else if (!mask[ch]) res_mask[ch]=false;
+      else {
+        res_mask[ch]=!invert; // no line by default
+        if (cli!=lines.end())
+          if (ch>=cli->first && ch<cli->second)
+            res_mask[ch]=invert; // this is a line
+      }
+      if (cli!=lines.end())
+        if (ch>=cli->second)
+          ++cli; // next line in the list
+    }
+    return res_mask;
   }
   catch (const AipsError &ae) {
-      throw;
+    throw;
   }
   catch (const exception &ex) {
-      throw AipsError(String("STLineFinder::getMask - STL error: ")+ex.what());
+    throw AipsError(String("STLineFinder::getMask - STL error: ")+ex.what());
   }
 }
@@ -1182 +1198 @@
                              const throw(casa::AipsError)
 {
-  // convert to required abscissa units
-  std::vector<double> vel=scan->getAbcissa(last_row_used);
+  std::vector<double> vel;
+  if (useScantable) {
+    // convert to required abscissa units
+    vel = scan->getAbcissa(last_row_used);
+  } else {
+    for (int i = 0; i < spectrum.nelements(); ++i)
+      vel.push_back((double)i);
+  }
   std::vector<int> ranges=getLineRangesInChannels();
   std::vector<double> res(ranges.size());
@@ -1202 +1224 @@
 {
   try {
-       if (scan.null())
-           throw AipsError("STLineFinder::getLineRangesInChannels - a scan should be set first,"
-                      " use set_scan followed by find_lines");
-       DebugAssert(mask.nelements()==scan->getChannels(last_row_used), AipsError);
-
-       if (!lines.size())
-           throw AipsError("STLineFinder::getLineRangesInChannels - one have to search for "
-                           "lines first, use find_lines");
-
-       std::vector<int> res(2*lines.size());
-       // iterator through lines & result
-       std::list<std::pair<int,int> >::const_iterator cli=lines.begin();
-       std::vector<int>::iterator ri=res.begin();
-       for (;cli!=lines.end() && ri!=res.end();++cli,++ri) {
-            *ri=cli->first;
-            if (++ri!=res.end())
-                *ri=cli->second-1;
-       }
-       return res;
-  }
-  catch (const AipsError &ae) {
-      throw;
-  }
-  catch (const exception &ex) {
-      throw AipsError(String("STLineFinder::getLineRanges - STL error: ")+ex.what());
+    if (useScantable) {
+      if (scan.null())
+        throw AipsError("STLineFinder::getLineRangesInChannels - a scan should be set first,"
+                        " use set_scan followed by find_lines");
+      DebugAssert(mask.nelements()==scan->getChannels(last_row_used), AipsError);
+    }
+
+    if (!lines.size())
+      throw AipsError("STLineFinder::getLineRangesInChannels - one have to search for "
+                      "lines first, use find_lines");
+
+    std::vector<int> res(2*lines.size());
+    // iterator through lines & result
+    std::list<std::pair<int,int> >::const_iterator cli = lines.begin();
+    std::vector<int>::iterator ri = res.begin();
+    for (; cli != lines.end() && ri != res.end(); ++cli,++ri) {
+      *ri = cli->first;
+      if (++ri != res.end())
+        *ri = cli->second - 1;
+    }
+    return res;
+  } catch (const AipsError &ae) {
+    throw;
+  } catch (const exception &ex) {
+    throw AipsError(String("STLineFinder::getLineRanges - STL error: ") + ex.what());
   }
 }

trunk/src/STLineFinder.h

@@ -164 +164 @@
    // set the scan to work with (in_scan parameter)
    void setScan(const ScantableWrapper &in_scan) throw(casa::AipsError);
+
+   // set spectrum data to work with. this is a method to allow linefinder work 
+   // without setting scantable for the purpose of using linefinder inside some 
+   // method in scantable class. (Dec. 22, 2010 by W.Kawasaki)
+   void setData(const std::vector<float> &in_spectrum);
 
    // search for spectral lines in a row specified by whichRow 
@@ -257 +262 @@
    // the lowest 80% of deviations (default)
    casa::Bool itsUseMedian;
+
+   // true if spectra and mask data should be provided from 
+   // scantable (default = true)
+   bool useScantable;
 };
 

trunk/src/Scantable.cpp

@@ -67 +67 @@
 #include "STPolStokes.h"
 #include "STAttr.h"
+#include "STLineFinder.h"
 #include "MathUtils.h"
 
@@ -1217 +1218 @@
 }
 **/
+
 void Scantable::setRestFrequencies( const vector<std::string>& name )
 {
@@ -1749 +1751 @@
   Vector<uChar> flags;
   flagsCol_.get(uInt(whichrow), flags);
-  for (int i = 0; i < flags.size(); i++) {
-    if (i==0) {
-      flag = flags[i];
-    }
-    else {
-      flag &= flags[i];
-    }
-    return ((flag >> 7) == 1);
-   }
-}
-
-void Scantable::doPolyBaseline(const std::vector<bool>& mask, int order, int rowno, Fitter& fitter)
-{
+  flag = flags[0];
+  for (int i = 1; i < flags.size(); ++i) {
+    flag &= flags[i];
+  }
+  return ((flag >> 7) == 1);
+}
+
+void Scantable::cubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, bool outLogger, const std::string& blfile) {
+  ofstream ofs;
+  String coordInfo;
+  bool hasSameNchan;
+  int firstIF;
+  bool outTextFile = false;
+
+  if (blfile != "") {
+    ofs.open(blfile.c_str(), ios::out | ios::app);
+    if (ofs) outTextFile = true;
+  }
+
+  if (outLogger || outTextFile) {
+    coordInfo = getCoordInfo()[0];
+    if (coordInfo == "") coordInfo = "channel";
+    hasSameNchan = hasSameNchanOverIFs();
+    firstIF = getIF(0);
+  }
+
+  //Fitter fitter = Fitter();
+  //fitter.setExpression("cspline", nPiece, thresClip, nIterClip);
+
+  int nRow = nrow();
+  std::vector<bool> chanMask;
+
+  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
+
+    chanMask = getCompositeChanMask(whichrow, mask);
+    //fitBaseline(chanMask, whichrow, fitter);
+    //setSpectrum(fitter.getResidual(), whichrow);
+    std::vector<int> pieceRanges;
+    std::vector<float> params;
+    std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, pieceRanges, params, nPiece, thresClip, nIterClip);
+    setSpectrum(res, whichrow);
+
+    if (outLogger || outTextFile) {
+      std::vector<bool>  fixed;
+      fixed.clear();
+      float rms = getRms(chanMask, whichrow);
+      String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan, firstIF, true);
+
+      if (outLogger) {
+        LogIO ols(LogOrigin("Scantable", "cubicSplineBaseline()", WHERE));
+        ols << formatPiecewiseBaselineParams(pieceRanges, params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
+      }
+      if (outTextFile) {
+        ofs << formatPiecewiseBaselineParams(pieceRanges, params, fixed, rms, masklist, whichrow, true) << flush;
+      }
+    }
+
+  }
+
+  if (outTextFile) ofs.close();
+}
+
+
+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)
+{
+  ofstream ofs;
+  String coordInfo;
+  bool hasSameNchan;
+  int firstIF;
+  bool outTextFile = false;
+
+  if (blfile != "") {
+    ofs.open(blfile.c_str(), ios::out | ios::app);
+    if (ofs) outTextFile = true;
+  }
+
+  if (outLogger || outTextFile) {
+    coordInfo = getCoordInfo()[0];
+    if (coordInfo == "") coordInfo = "channel";
+    hasSameNchan = hasSameNchanOverIFs();
+    firstIF = getIF(0);
+  }
+
+  //Fitter fitter = Fitter();
+  //fitter.setExpression("cspline", nPiece, thresClip, nIterClip);
+
+  int nRow = nrow();
+  std::vector<bool> chanMask;
+  int minEdgeSize = getIFNos().size()*2;
+  STLineFinder lineFinder = STLineFinder();
+  lineFinder.setOptions(threshold, 3, chanAvgLimit);
+
+  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
+
+    //-------------------------------------------------------
+    //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
+    //-------------------------------------------------------
+    int edgeSize = edge.size();
+    std::vector<int> currentEdge;
+    if (edgeSize >= 2) {
+      int idx = 0;
+      if (edgeSize > 2) {
+        if (edgeSize < minEdgeSize) {
+          throw(AipsError("Length of edge element info is less than that of IFs"));
+        }
+        idx = 2 * getIF(whichrow);
+      }
+      currentEdge.push_back(edge[idx]);
+      currentEdge.push_back(edge[idx+1]);
+    } else {
+      throw(AipsError("Wrong length of edge element"));
+    }
+    lineFinder.setData(getSpectrum(whichrow));
+    lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
+    chanMask = lineFinder.getMask();
+    //-------------------------------------------------------
+
+
+    //fitBaseline(chanMask, whichrow, fitter);
+    //setSpectrum(fitter.getResidual(), whichrow);
+    std::vector<int> pieceRanges;
+    std::vector<float> params;
+    std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, pieceRanges, params, nPiece, thresClip, nIterClip);
+    setSpectrum(res, whichrow);
+
+    if (outLogger || outTextFile) {
+      std::vector<bool> fixed;
+      fixed.clear();
+      float rms = getRms(chanMask, whichrow);
+      String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan, firstIF, true);
+
+      if (outLogger) {
+        LogIO ols(LogOrigin("Scantable", "autoCubicSplineBaseline()", WHERE));
+        ols << formatPiecewiseBaselineParams(pieceRanges, params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
+      }
+      if (outTextFile) {
+        ofs << formatPiecewiseBaselineParams(pieceRanges, params, fixed, rms, masklist, whichrow, true) << flush;
+      }
+    }
+
+  }
+
+  if (outTextFile) ofs.close();
+}
+
+
+std::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) {
+  if (nPiece < 1) {
+    throw(AipsError("wrong number of the sections for fitting"));
+  }
+  if (data.size() != mask.size()) {
+    throw(AipsError("data and mask have different size"));
+  }
+
+  int nChan = data.size();
+  std::vector<int> maskArray;
+  std::vector<int> x;
+  for (int i = 0; i < nChan; ++i) {
+    maskArray.push_back(mask[i] ? 1 : 0);
+    if (mask[i]) {
+      x.push_back(i);
+    }
+  }
+
+  int nData = x.size();
+  int nElement = (int)(floor(floor((double)(nData/nPiece))+0.5));
+
+  std::vector<double> sectionList0, sectionList1, sectionListR;
+  sectionList0.push_back((double)x[0]);
+  sectionRanges.clear();
+  sectionRanges.push_back(x[0]);
+  for (int i = 1; i < nPiece; ++i) {
+    double valX = (double)x[nElement*i];
+    sectionList0.push_back(valX);
+    sectionList1.push_back(valX);
+    sectionListR.push_back(1.0/valX);
+
+    sectionRanges.push_back(x[nElement*i]-1);
+    sectionRanges.push_back(x[nElement*i]);
+  }
+  sectionList1.push_back((double)(x[x.size()-1]+1));
+  sectionRanges.push_back(x[x.size()-1]);
+
+  std::vector<double> x1, x2, x3, z1, x1z1, x2z1, x3z1, r1;
+  for (int i = 0; i < nChan; ++i) {
+    x1.push_back((double)i);
+    x2.push_back((double)(i*i));
+    x3.push_back((double)(i*i*i));
+    z1.push_back((double)data[i]);
+    x1z1.push_back(((double)i)*(double)data[i]);
+    x2z1.push_back(((double)(i*i))*(double)data[i]);
+    x3z1.push_back(((double)(i*i*i))*(double)data[i]);
+    r1.push_back(0.0);
+  }
+
+  int currentNData = nData;
+  int nDOF = nPiece + 3;  //number of parameters to solve, namely, 4+(nPiece-1).
+
+  for (int nClip = 0; nClip < nIterClip+1; ++nClip) {
+    //xMatrix : horizontal concatenation of the least-sq. matrix (left) and an 
+    //identity matrix (right).
+    //the right part is used to calculate the inverse matrix of the left part. 
+    double xMatrix[nDOF][2*nDOF];
+    double zMatrix[nDOF];
+    for (int i = 0; i < nDOF; ++i) {
+      for (int j = 0; j < 2*nDOF; ++j) {
+        xMatrix[i][j] = 0.0;
+      }
+      xMatrix[i][nDOF+i] = 1.0;
+      zMatrix[i] = 0.0;
+    }
+
+    for (int n = 0; n < nPiece; ++n) {
+      for (int i = sectionList0[n]; i < sectionList1[n]; ++i) {
+        if (maskArray[i] == 0) continue;
+        xMatrix[0][0] += 1.0;
+        xMatrix[0][1] += x1[i];
+        xMatrix[0][2] += x2[i];
+        xMatrix[0][3] += x3[i];
+        xMatrix[1][1] += x2[i];
+        xMatrix[1][2] += x3[i];
+        xMatrix[1][3] += x2[i]*x2[i];
+        xMatrix[2][2] += x2[i]*x2[i];
+        xMatrix[2][3] += x3[i]*x2[i];
+        xMatrix[3][3] += x3[i]*x3[i];
+        zMatrix[0] += z1[i];
+        zMatrix[1] += x1z1[i];
+        zMatrix[2] += x2z1[i];
+        zMatrix[3] += x3z1[i];
+        for (int j = 0; j < n; ++j) {
+          double q = 1.0 - x1[i]*sectionListR[j];
+          q = q*q*q;
+          xMatrix[0][j+4] += q;
+          xMatrix[1][j+4] += q*x1[i];
+          xMatrix[2][j+4] += q*x2[i];
+          xMatrix[3][j+4] += q*x3[i];
+          for (int k = 0; k < j; ++k) {
+            double r = 1.0 - x1[i]*sectionListR[k];
+            r = r*r*r;
+            xMatrix[k+4][j+4] += r*q;
+          }
+          xMatrix[j+4][j+4] += q*q;
+          zMatrix[j+4] += q*z1[i];
+        }
+      }
+    }
+
+    for (int i = 0; i < nDOF; ++i) {
+      for (int j = 0; j < i; ++j) {
+        xMatrix[i][j] = xMatrix[j][i];
+      }
+    }
+
+    std::vector<double> invDiag;
+    for (int i = 0; i < nDOF; ++i) {
+      invDiag.push_back(1.0/xMatrix[i][i]);
+      for (int j = 0; j < nDOF; ++j) {
+        xMatrix[i][j] *= invDiag[i];
+      }
+    }
+
+    for (int k = 0; k < nDOF; ++k) {
+      for (int i = 0; i < nDOF; ++i) {
+        if (i != k) {
+          double factor1 = xMatrix[k][k];
+          double factor2 = xMatrix[i][k];
+          for (int j = k; j < 2*nDOF; ++j) {
+            xMatrix[i][j] *= factor1;
+            xMatrix[i][j] -= xMatrix[k][j]*factor2;
+            xMatrix[i][j] /= factor1;
+          }
+        }
+      }
+      double xDiag = xMatrix[k][k];
+      for (int j = k; j < 2*nDOF; ++j) {
+        xMatrix[k][j] /= xDiag;
+      }
+    }
+    
+    for (int i = 0; i < nDOF; ++i) {
+      for (int j = 0; j < nDOF; ++j) {
+        xMatrix[i][nDOF+j] *= invDiag[j];
+      }
+    }
+    //compute a vector y which consists of the coefficients of the best-fit spline curves 
+    //(a0,a1,a2,a3(,b3,c3,...)), namely, the ones for the leftmost piece and the ones of 
+    //cubic terms for the other pieces (in case nPiece>1).
+    std::vector<double> y;
+    for (int i = 0; i < nDOF; ++i) {
+      y.push_back(0.0);
+      for (int j = 0; j < nDOF; ++j) {
+        y[i] += xMatrix[i][nDOF+j]*zMatrix[j];
+      }
+    }
+
+    double a0 = y[0];
+    double a1 = y[1];
+    double a2 = y[2];
+    double a3 = y[3];
+    params.clear();
+
+    for (int n = 0; n < nPiece; ++n) {
+      for (int i = sectionList0[n]; i < sectionList1[n]; ++i) {
+        r1[i] = a0 + a1*x1[i] + a2*x2[i] + a3*x3[i];
+      }
+      params.push_back(a0);
+      params.push_back(a1);
+      params.push_back(a2);
+      params.push_back(a3);
+
+      if (n == nPiece-1) break;
+
+      double d = y[4+n];
+      a0 += d;
+      a1 -= 3.0*d*sectionListR[n];
+      a2 += 3.0*d*sectionListR[n]*sectionListR[n];
+      a3 -= d*sectionListR[n]*sectionListR[n]*sectionListR[n];
+    }
+
+    if ((nClip == nIterClip) || (thresClip <= 0.0)) {
+      break;
+    } else {
+      std::vector<double> rz;
+      double stdDev = 0.0;
+      for (int i = 0; i < nChan; ++i) {
+        double val = abs(r1[i] - z1[i]);
+        rz.push_back(val);
+        stdDev += val*val*(double)maskArray[i];
+      }
+      stdDev = sqrt(stdDev/(double)nData);
+      
+      double thres = stdDev * thresClip;
+      int newNData = 0;
+      for (int i = 0; i < nChan; ++i) {
+        if (rz[i] >= thres) {
+          maskArray[i] = 0;
+        }
+        if (maskArray[i] > 0) {
+          newNData++;
+        }
+      }
+      if (newNData == currentNData) {
+        break;                   //no additional flag. finish iteration.
+      } else {
+        currentNData = newNData;
+      }
+    }
+  }
+
+  std::vector<float> residual;
+  for (int i = 0; i < nChan; ++i) {
+    residual.push_back((float)(z1[i] - r1[i]));
+  }
+  return residual;
+
+}
+
+void Scantable::fitBaseline(const std::vector<bool>& mask, int whichrow, Fitter& fitter)
+{
+  std::vector<double> abcsd = getAbcissa(whichrow);
+  std::vector<float> abcs;
+  for (int i = 0; i < abcsd.size(); ++i) {
+    abcs.push_back((float)abcsd[i]);
+  }
+  std::vector<float> spec = getSpectrum(whichrow);
+
+  fitter.setData(abcs, spec, mask);
+  fitter.lfit();
+}
+
+std::vector<bool> Scantable::getCompositeChanMask(int whichrow, const std::vector<bool>& inMask)
+{
+  std::vector<bool> chanMask = getMask(whichrow);
+  int chanMaskSize = chanMask.size();
+  if (chanMaskSize != inMask.size()) {
+    throw(AipsError("different mask sizes"));
+  }
+  for (int i = 0; i < chanMaskSize; ++i) {
+    chanMask[i] = chanMask[i] && inMask[i];
+  }
+
+  return chanMask;
+}
+
+/*
+std::vector<bool> Scantable::getCompositeChanMask(int whichrow, const std::vector<bool>& inMask, const std::vector<int>& edge, const int minEdgeSize, STLineFinder& lineFinder)
+{
+  int edgeSize = edge.size();
+  std::vector<int> currentEdge;
+  if (edgeSize >= 2) {
+      int idx = 0;
+      if (edgeSize > 2) {
+        if (edgeSize < minEdgeSize) {
+          throw(AipsError("Length of edge element info is less than that of IFs"));
+        }
+        idx = 2 * getIF(whichrow);
+      }
+      currentEdge.push_back(edge[idx]);
+      currentEdge.push_back(edge[idx+1]);
+  } else {
+    throw(AipsError("Wrong length of edge element"));
+  }
+
+  lineFinder.setData(getSpectrum(whichrow));
+  lineFinder.findLines(getCompositeChanMask(whichrow, inMask), currentEdge, whichrow);
+
+  return lineFinder.getMask();
+}
+*/
+
+void Scantable::polyBaseline(const std::vector<bool>& mask, int order, bool outLogger, const std::string& blfile)
+{
+  ofstream ofs;
+  String coordInfo;
+  bool hasSameNchan;
+  int firstIF;
+  bool outTextFile = false;
+
+  if (blfile != "") {
+    ofs.open(blfile.c_str(), ios::out | ios::app);
+    if (ofs) outTextFile = true;
+  }
+
+  if (outLogger || outTextFile) {
+    coordInfo = getCoordInfo()[0];
+    if (coordInfo == "") coordInfo = "channel";
+    hasSameNchan = hasSameNchanOverIFs();
+    firstIF = getIF(0);
+  }
+
+  Fitter fitter = Fitter();
   fitter.setExpression("poly", order);
 
-  std::vector<double> abcsd = getAbcissa(rowno);
-  std::vector<float> abcs;
-  for (int i = 0; i < abcsd.size(); i++) {
-    abcs.push_back((float)abcsd[i]);
-  }
-  std::vector<float> spec = getSpectrum(rowno);
+  int nRow = nrow();
+  std::vector<bool> chanMask;
+
+  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
+
+    chanMask = getCompositeChanMask(whichrow, mask);
+    fitBaseline(chanMask, whichrow, fitter);
+    setSpectrum(fitter.getResidual(), whichrow);
+    
+    if (outLogger || outTextFile) {
+      std::vector<float> params = fitter.getParameters();
+      std::vector<bool>  fixed  = fitter.getFixedParameters();
+      float rms = getRms(chanMask, whichrow);
+      String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan, firstIF, true);
+
+      if (outLogger) {
+        LogIO ols(LogOrigin("Scantable", "polyBaseline()", WHERE));
+        ols << formatBaselineParams(params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
+      }
+      if (outTextFile) {
+        ofs << formatBaselineParams(params, fixed, rms, masklist, whichrow, true) << flush;
+      }
+    }
+
+  }
+
+  if (outTextFile) ofs.close();
+}
+
+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)
+{
+  ofstream ofs;
+  String coordInfo;
+  bool hasSameNchan;
+  int firstIF;
+  bool outTextFile = false;
+
+  if (blfile != "") {
+    ofs.open(blfile.c_str(), ios::out | ios::app);
+    if (ofs) outTextFile = true;
+  }
+
+  if (outLogger || outTextFile) {
+    coordInfo = getCoordInfo()[0];
+    if (coordInfo == "") coordInfo = "channel";
+    hasSameNchan = hasSameNchanOverIFs();
+    firstIF = getIF(0);
+  }
+
+  Fitter fitter = Fitter();
+  fitter.setExpression("poly", order);
+
+  int nRow = nrow();
+  std::vector<bool> chanMask;
+  int minEdgeSize = getIFNos().size()*2;
+  STLineFinder lineFinder = STLineFinder();
+  lineFinder.setOptions(threshold, 3, chanAvgLimit);
+
+  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
+
+    //-------------------------------------------------------
+    //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
+    //-------------------------------------------------------
+    int edgeSize = edge.size();
+    std::vector<int> currentEdge;
+    if (edgeSize >= 2) {
+      int idx = 0;
+      if (edgeSize > 2) {
+        if (edgeSize < minEdgeSize) {
+          throw(AipsError("Length of edge element info is less than that of IFs"));
+        }
+        idx = 2 * getIF(whichrow);
+      }
+      currentEdge.push_back(edge[idx]);
+      currentEdge.push_back(edge[idx+1]);
+    } else {
+      throw(AipsError("Wrong length of edge element"));
+    }
+    lineFinder.setData(getSpectrum(whichrow));
+    lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
+    chanMask = lineFinder.getMask();
+    //-------------------------------------------------------
+
+
+    fitBaseline(chanMask, whichrow, fitter);
+    setSpectrum(fitter.getResidual(), whichrow);
+
+    if (outLogger || outTextFile) {
+      std::vector<float> params = fitter.getParameters();
+      std::vector<bool>  fixed  = fitter.getFixedParameters();
+      float rms = getRms(chanMask, whichrow);
+      String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan, firstIF, true);
+
+      if (outLogger) {
+        LogIO ols(LogOrigin("Scantable", "autoPolyBaseline()", WHERE));
+        ols << formatBaselineParams(params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
+      }
+      if (outTextFile) {
+        ofs << formatBaselineParams(params, fixed, rms, masklist, whichrow, true) << flush;
+      }
+    }
+
+  }
+
+  if (outTextFile) ofs.close();
+}
+
+
+float Scantable::getRms(const std::vector<bool>& mask, int whichrow) {
+  Vector<Float> spec;
+  specCol_.get(whichrow, spec);
+
+  float mean = 0.0;
+  float smean = 0.0;
+  int n = 0;
+  for (int i = 0; i < spec.nelements(); ++i) {
+    if (mask[i]) {
+      mean += spec[i];
+      smean += spec[i]*spec[i];
+      n++;
+    }
+  }
+
+  mean /= (float)n;
+  smean /= (float)n;
+
+  return sqrt(smean - mean*mean);
+}
+
+
+std::string Scantable::formatBaselineParamsHeader(int whichrow, const std::string& masklist, bool verbose) const
+{
+  ostringstream oss;
+
+  if (verbose) {
+    for (int i = 0; i < 60; ++i) {
+      oss << "-";
+    }
+    oss << endl;
+    oss <<  " Scan[" << getScan(whichrow)  << "]";
+    oss <<  " Beam[" << getBeam(whichrow)  << "]";
+    oss <<    " IF[" << getIF(whichrow)    << "]";
+    oss <<   " Pol[" << getPol(whichrow)   << "]";
+    oss << " Cycle[" << getCycle(whichrow) << "]: " << endl;
+    oss << "Fitter range = " << masklist << endl;
+    oss << "Baseline parameters" << endl;
+    oss << flush;
+  }
+
+  return String(oss);
+}
+
+std::string Scantable::formatBaselineParamsFooter(float rms, bool verbose) const
+{
+  ostringstream oss;
+
+  if (verbose) {
+    oss << endl;
+    oss << "Results of baseline fit" << endl;
+    oss << "  rms = " << setprecision(6) << rms << endl;
+  }
+
+  return String(oss);
+}
+
+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
+{
+  ostringstream oss;
+  oss << formatBaselineParamsHeader(whichrow, masklist, verbose);
+
+  if ((params.size() > 0) && (ranges.size() > 0)) {
+    if (((ranges.size() % 2) == 0) && ((params.size() % (ranges.size() / 2)) == 0)) {
+      int nParam = params.size() / (ranges.size() / 2);
+      for (uInt j = 0; j < ranges.size(); j+=2) {
+        oss << "[" << ranges[j] << "," << ranges[j+1] << "]";
+        for (uInt i = 0; i < nParam; ++i) {
+          if (i > 0) {
+            oss << ",";
+          }
+          String fix = ((fixed.size() > 0) && (fixed[i]) && verbose) ? "(fixed)" : "";
+          oss << "  p" << i << fix << "= " << setprecision(6) << params[j/2*nParam+i];
+        }
+      }
+    } else {
+      oss << "  ";
+    }
+  } else {
+    oss << "  Not fitted";
+  }
+
+  oss << formatBaselineParamsFooter(rms, verbose);
+  oss << flush;
+
+  return String(oss);
+}
+
+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
+{
+  ostringstream oss;
+  oss << formatBaselineParamsHeader(whichrow, masklist, verbose);
+
+  if (params.size() > 0) {
+    for (uInt i = 0; i < params.size(); ++i) {
+      if (i > 0) {
+        oss << ",";
+      }
+      String fix = ((fixed.size() > 0) && (fixed[i]) && verbose) ? "(fixed)" : "";
+      oss << "  p" << i << fix << "= " << setprecision(6) << params[i];
+    }
+  } else {
+    oss << "  Not fitted";
+  }
+
+  oss << formatBaselineParamsFooter(rms, verbose);
+  oss << flush;
+
+  return String(oss);
+}
+
+std::string Scantable::getMaskRangeList(const std::vector<bool>& mask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, int firstIF, bool silent)
+{
+  if (mask.size() < 2) {
+    throw(AipsError("The mask elements should be > 1"));
+  }
+  if (mask.size() != nchan()) {
+    throw(AipsError("Number of channels in scantable != number of mask elements"));
+  }
+
+  std::vector<int> startIdx = getMaskEdgeIndices(mask, true);
+  std::vector<int> endIdx   = getMaskEdgeIndices(mask, false);
+
+  if (startIdx.size() != endIdx.size()) {
+    throw(AipsError("Numbers of mask start and mask end are not identical"));
+  }
+  for (int i = 0; i < startIdx.size(); ++i) {
+    if (startIdx[i] > endIdx[i]) {
+      throw(AipsError("Mask start index > mask end index"));
+    }
+  }
+
+  if (!silent) {
+    LogIO logOs(LogOrigin("Scantable", "getMaskRangeList()", WHERE));
+    logOs << LogIO::WARN << "The current mask window unit is " << coordInfo;
+    if (!hasSameNchan) {
+      logOs << endl << "This mask is only valid for IF=" << firstIF;
+    }
+    logOs << LogIO::POST;
+  }
+
+  std::vector<double> abcissa = getAbcissa(whichrow);
+  ostringstream oss;
+  oss.setf(ios::fixed);
+  oss << setprecision(1) << "[";
+  for (int i = 0; i < startIdx.size(); ++i) {
+    std::vector<float> aMaskRange;
+    aMaskRange.push_back((float)abcissa[startIdx[i]]);
+    aMaskRange.push_back((float)abcissa[endIdx[i]]);
+    sort(aMaskRange.begin(), aMaskRange.end());
+    if (i > 0) oss << ",";
+    oss << "[" << aMaskRange[0] << "," << aMaskRange[1] << "]";
+  }
+  oss << "]" << flush;
+
+  return String(oss);
+}
+
+bool Scantable::hasSameNchanOverIFs()
+{
+  int nIF = nif(-1);
+  int nCh;
+  int totalPositiveNChan = 0;
+  int nPositiveNChan = 0;
+
+  for (int i = 0; i < nIF; ++i) {
+    nCh = nchan(i);
+    if (nCh > 0) {
+      totalPositiveNChan += nCh;
+      nPositiveNChan++;
+    }
+  }
+
+  return (totalPositiveNChan == (nPositiveNChan * nchan(0)));
+}
+
+std::vector<int> Scantable::getMaskEdgeIndices(const std::vector<bool>& mask, bool getStartIndices)
+{
+  if (mask.size() < 2) {
+    throw(AipsError("The mask elements should be > 1"));
+  }
+  if (mask.size() != nchan()) {
+    throw(AipsError("Number of channels in scantable != number of mask elements"));
+  }
+
+  std::vector<int> out;
+  int endIdx = mask.size() - 1;
+
+  if (getStartIndices) {
+    if (mask[0]) {
+      out.push_back(0);
+    }
+    for (int i = 0; i < endIdx; ++i) {
+      if ((!mask[i]) && mask[i+1]) {
+        out.push_back(i+1);
+      }
+    }
+  } else {
+    for (int i = 0; i < endIdx; ++i) {
+      if (mask[i] && (!mask[i+1])) {
+        out.push_back(i);
+      }
+    }
+    if (mask[endIdx]) {
+      out.push_back(endIdx);
+    }
+  }
+
+  return out;
+}
+
+
+/*
+STFitEntry Scantable::polyBaseline(const std::vector<bool>& mask, int order, int rowno)
+{
+  Fitter fitter = Fitter();
+  fitter.setExpression("poly", order);
+
   std::vector<bool> fmask = getMask(rowno);
   if (fmask.size() != mask.size()) {
     throw(AipsError("different mask sizes"));
   }
-  for (int i = 0; i < fmask.size(); i++) {
+  for (int i = 0; i < fmask.size(); ++i) {
     fmask[i] = fmask[i] && mask[i];
   }
-  fitter.setData(abcs, spec, fmask);
-
-  fitter.lfit();
-}
-
-void Scantable::polyBaselineBatch(const std::vector<bool>& mask, int order)
-{
-  Fitter fitter = Fitter();
-  for (uInt rowno=0; rowno < nrow(); ++rowno) {
-    doPolyBaseline(mask, order, rowno, fitter);
-    setSpectrum(fitter.getResidual(), rowno);
-  }
-}
-
-STFitEntry Scantable::polyBaseline(const std::vector<bool>& mask, int order, int rowno)
-{
-  Fitter fitter = Fitter();
-  doPolyBaseline(mask, order, rowno, fitter);
+
+  fitBaseline(fmask, rowno, fitter);
   setSpectrum(fitter.getResidual(), rowno);
   return fitter.getFitEntry();
 }
+*/
 
 }

trunk/src/Scantable.h

@@ -16 +16 @@
 #include <string>
 #include <vector>
+#include <iostream>
+#include <fstream>
 // AIPS++
 #include <casa/aips.h>
@@ -493 +495 @@
   bool getFlagtraFast(int whichrow);
 
-  void polyBaselineBatch(const std::vector<bool>& mask, int order);
-  STFitEntry polyBaseline(const std::vector<bool>& mask, int order, int rowno);
+  void polyBaseline(const std::vector<bool>& mask, 
+                    int order, 
+                    bool outLogger=false, 
+                    const std::string& blfile="");
+  void autoPolyBaseline(const std::vector<bool>& mask,
+                        int order, 
+                        const std::vector<int>& edge, 
+                        float threshold=3.0, 
+                        int chanAvgLimit=1, 
+                        bool outLogger=false, 
+                        const std::string& blfile="");
+  void cubicSplineBaseline(const std::vector<bool>& mask,
+                           int nPiece,
+                           float thresClip, 
+                           int nIterClip,
+                           bool outLogger=false,
+                           const std::string& blfile="");
+  void autoCubicSplineBaseline(const std::vector<bool>& mask,
+                               int nPiece,
+                               float thresClip, 
+                               int nIterClip,
+                               const std::vector<int>& edge, 
+                               float threshold=3.0, 
+                               int chanAvgLimit=1, 
+                               bool outLogger=false, 
+                               const std::string& blfile="");
+  float getRms(const std::vector<bool>& mask, int whichrow);
+  std::string formatBaselineParams(const std::vector<float>& params, 
+                                   const std::vector<bool>& fixed, 
+                                   float rms, 
+                                   const std::string& masklist, 
+                                   int whichrow, 
+                                   bool verbose=false) const;
+  std::string 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=false) const;
+
 
 private:
@@ -608 +649 @@
                                                       const casa::Array<T2>&);
 
-  void doPolyBaseline(const std::vector<bool>& mask, int order, int rowno, Fitter& fitter);
+  void fitBaseline(const std::vector<bool>& mask, int whichrow, Fitter& fitter);
+  std::vector<float> doCubicSplineFitting(const std::vector<float>& data, 
+                                          const std::vector<bool>& mask,
+                                          std::vector<int>& sectionRanges,
+                                          std::vector<float>& params,
+                                          int nPiece=2,
+                                          float thresClip=3.0, 
+                                          int nIterClip=1);
+  bool hasSameNchanOverIFs();
+  std::string getMaskRangeList(const std::vector<bool>& mask, 
+                                int whichrow, 
+                                const casa::String& coordInfo, 
+                                bool hasSameNchan, 
+                                int firstIF, 
+                                bool silent=false);
+  std::vector<int> getMaskEdgeIndices(const std::vector<bool>& mask, bool getStartIndices=true);
+  std::string formatBaselineParamsHeader(int whichrow, const std::string& masklist, bool verbose) const;
+  std::string formatBaselineParamsFooter(float rms, bool verbose) const;
+  std::vector<bool> getCompositeChanMask(int whichrow, const std::vector<bool>& inMask);
+  //std::vector<bool> getCompositeChanMask(int whichrow, const std::vector<bool>& inMask, const std::vector<int>& edge, const int minEdgeSize, STLineFinder& lineFinder);
 
   void applyChanFlag( casa::uInt whichrow, const std::vector<bool>& msk, casa::uChar flagval);

trunk/src/ScantableWrapper.h

@@ -230 +230 @@
   { return table_->getFit(whichrow); }
 
-  void calculateAZEL() { table_->calculateAZEL(); };
+  void calculateAZEL() { table_->calculateAZEL(); }
 
   std::vector<std::string> columnNames() const
@@ -257 +257 @@
   { table_->reshapeSpectrum( nmin, nmax ); }
 
-  STFitEntry polyBaseline(const std::vector<bool>& mask, int order, int rowno)
-  { return table_->polyBaseline(mask, order, rowno); }
-
-  void polyBaselineBatch(const std::vector<bool>& mask, int order)
-  { table_->polyBaselineBatch(mask, order); }
+  void polyBaseline(const std::vector<bool>& mask, int order, bool outlog=false, const std::string& blfile="")
+  { table_->polyBaseline(mask, order, outlog, blfile); }
+
+  void autoPolyBaseline(const std::vector<bool>& mask, int order, const std::vector<int>& edge, float threshold=5.0, int chan_avg_limit=1, bool outlog=false, const std::string& blfile="")
+  { table_->autoPolyBaseline(mask, order, edge, threshold, chan_avg_limit, outlog, blfile); }
+
+  void cubicSplineBaseline(const std::vector<bool>& mask, int npiece, float clipthresh, int clipniter, bool outlog=false, const std::string& blfile="")
+  { table_->cubicSplineBaseline(mask, npiece, clipthresh, clipniter, outlog, blfile); }
+
+  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 outlog=false, const std::string& blfile="")
+  { table_->autoCubicSplineBaseline(mask, npiece, clipthresh, clipniter, edge, threshold, chan_avg_limit, outlog, blfile); }
+
+  float getRms(const std::vector<bool>& mask, int whichrow)
+  { return table_->getRms(mask, whichrow); }
+
+  std::string formatBaselineParams(const std::vector<float>& params, const std::vector<bool>& fixed, float rms, const std::string& masklist, int whichrow, bool verbose=false)
+  { return table_->formatBaselineParams(params, fixed, rms, masklist, whichrow, verbose); }
+
+  std::string 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=false)
+  { return table_->formatPiecewiseBaselineParams(ranges, params, fixed, rms, masklist, whichrow, verbose); }
 
   bool getFlagtraFast(int whichrow=0) const
-    { return table_->getFlagtraFast(whichrow); }
+  { return table_->getFlagtraFast(whichrow); }
+
 
 

trunk/src/python_STLineFinder.cpp

@@ -42 +42 @@
          .def("setoptions",&STLineFinder::setOptions)
          .def("setscan",&STLineFinder::setScan)
+         .def("setdata",&STLineFinder::setData)
          .def("findlines",&STLineFinder::findLines)
          .def("getmask",&STLineFinder::getMask)

trunk/src/python_Scantable.cpp

@@ -143 +143 @@
           boost::python::arg("nmax")=-1) )
     .def("_poly_baseline", &ScantableWrapper::polyBaseline)
-    .def("_poly_baseline_batch", &ScantableWrapper::polyBaselineBatch)
+    .def("_auto_poly_baseline", &ScantableWrapper::autoPolyBaseline)
+    .def("_cspline_baseline", &ScantableWrapper::cubicSplineBaseline)
+    .def("_auto_cspline_baseline", &ScantableWrapper::autoCubicSplineBaseline)
+    .def("get_rms", &ScantableWrapper::getRms)
+    .def("format_blparams_row", &ScantableWrapper::formatBaselineParams)
+    .def("format_piecewise_blparams_row", &ScantableWrapper::formatPiecewiseBaselineParams)
     .def("_getflagtrafast", &ScantableWrapper::getFlagtraFast,
          (boost::python::arg("whichrow")=0) )
+    //.def("_sspline_baseline", &ScantableWrapper::smoothingSplineBaseline,
+    // (boost::python::arg("thres")=3.0,
+    //  boost::python::arg("niter")=1) )
+    //.def("_test_cin", &ScantableWrapper::testCin)
   ;
 };