Changeset 2047

Timestamp:

03/15/11 18:31:04 (13 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:

• python/asapfitter.py (modified) (10 diffs)
• python/scantable.py (modified) (1 diff)
• src/STFitter.cpp (modified) (5 diffs)
• src/Scantable.cpp (modified) (22 diffs)
• src/Scantable.h (modified) (3 diffs)
• src/ScantableWrapper.h (modified) (2 diffs)
• src/python_Scantable.cpp (modified) (1 diff)

trunk/python/asapfitter.py

@@ -73 +73 @@
         Set the function to be fit.
         Parameters:
-            poly:    use a polynomial of the order given with nonlinear least squares fit
-            lpoly:   use polynomial of the order given with linear least squares fit
-            gauss:   fit the number of gaussian specified
-            lorentz: fit the number of lorentzian specified
+            poly:     use a polynomial of the order given with nonlinear least squares fit
+            lpoly:    use polynomial of the order given with linear least squares fit
+            gauss:    fit the number of gaussian specified
+            lorentz:  fit the number of lorentzian specified
+            sinusoid: fit the number of sinusoid specified
         Example:
             fitter.set_function(poly=3)  # will fit a 3rd order polynomial via nonlinear method
@@ -82 +83 @@
             fitter.set_function(gauss=2) # will fit two gaussians
             fitter.set_function(lorentz=2) # will fit two lorentzians
+            fitter.set_function(sinusoid=3) # will fit three sinusoids
         """
         #default poly order 0
@@ -109 +111 @@
             self.components = [ 3 for i in range(n) ]
             self.uselinear = False
+        elif kwargs.has_key('sinusoid'):
+            n = kwargs.get('sinusoid')
+            self.fitfunc = 'sinusoid'
+            self.fitfuncs = [ 'sinusoid' for i in range(n) ]
+            self.components = [ 3 for i in range(n) ]
+            self.uselinear = False
         else:
             msg = "Invalid function type."
@@ -204 +212 @@
               fixed:     a vector of which parameters are to be held fixed
                          (default is none)
-              component: in case of multiple gaussians/lorentzians,
-                         the index of the component
+              component: in case of multiple gaussians/lorentzians/sinusoidals,
+                         the index of the target component
         """
         component = None
@@ -220 +228 @@
             msg = "Please specify a fitting function first."
             raise RuntimeError(msg)
-        if (self.fitfunc == "gauss" or self.fitfunc == 'lorentz') and component is not None:
+        if (self.fitfunc == "gauss" or self.fitfunc == "lorentz" or self.fitfunc == "sinusoid") and component is not None:
             if not self.fitted and sum(self.fitter.getparameters()) == 0:
                 pars = _n_bools(len(self.components)*3, False)
-                fxd = _n_bools(len(pars), False)
+                fxd  = _n_bools(len(pars), False)
             else:
                 pars = list(self.fitter.getparameters())
-                fxd = list(self.fitter.getfixedparameters())
+                fxd  = list(self.fitter.getfixedparameters())
             i = 3*component
             pars[i:i+3] = params
-            fxd[i:i+3] = fixed
+            fxd[i:i+3]  = fixed
             params = pars
-            fixed = fxd
+            fixed  = fxd
         self.fitter.setparameters(params)
         if fixed is not None:
@@ -290 +298 @@
             d = {'params':[peak, centre, fwhm],
                  'fixed':[peakfixed, centrefixed, fwhmfixed]}
+            self.set_parameters(d, component)
+        else:
+            msg = "Please select a valid  component."
+            raise ValueError(msg)
+
+    @asaplog_post_dec
+    def set_sinusoid_parameters(self, ampl, period, x0,
+                             amplfixed=0, periodfixed=0,
+                             x0fixed=0,
+                             component=0):
+        """
+        Set the Parameters of a 'Sinusoidal' component, set with set_function.
+        Parameters:
+            ampl, period, x0:  The sinusoidal parameters
+            amplfixed,
+            periodfixed,
+            x0fixed:             Optional parameters to indicate if
+                                 the paramters should be held fixed during
+                                 the fitting process. The default is to keep
+                                 all parameters flexible.
+            component:           The number of the component (Default is the
+                                 component 0)
+        """
+        if self.fitfunc != "sinusoid":
+            msg = "Function only operates on Sinusoidal components."
+            raise ValueError(msg)
+        if 0 <= component < len(self.components):
+            d = {'params':[ampl, period, x0],
+                 'fixed': [amplfixed, periodfixed, x0fixed]}
             self.set_parameters(d, component)
         else:
@@ -338 +375 @@
         cerrs = errs
         if component is not None:
-            if self.fitfunc == "gauss" or self.fitfunc == "lorentz":
+            if self.fitfunc == "gauss" or self.fitfunc == "lorentz" or self.fitfunc == "sinusoid":
                 i = 3*component
                 if i < len(errs):
@@ -361 +398 @@
         area = []
         if component is not None:
-            if self.fitfunc == "gauss" or self.fitfunc == "lorentz":
+            if self.fitfunc == "poly" or self.fitfunc == "lpoly":
+                cpars = pars
+                cfixed = fixed
+                cerrs = errs
+            else:
                 i = 3*component
                 cpars = pars[i:i+3]
                 cfixed = fixed[i:i+3]
                 cerrs = errs[i:i+3]
-                a = self.get_area(component)
-                area = [a for i in range(3)]
-            else:
-                cpars = pars
-                cfixed = fixed
-                cerrs = errs
+                if self.fitfunc == "gauss" or self.fitfunc == "lorentz":
+                    a = self.get_area(component)
+                    area = [a for i in range(3)]
         else:
             cpars = pars
@@ -378 +416 @@
             if self.fitfunc == "gauss" or self.fitfunc == "lorentz":
                 for c in range(len(self.components)):
-                  a = self.get_area(c)
-                  area += [a for i in range(3)]
+                    a = self.get_area(c)
+                    area += [a for i in range(3)]
         fpars = self._format_pars(cpars, cfixed, errors and cerrs, area)
         asaplog.push(fpars)
@@ -387 +425 @@
     def _format_pars(self, pars, fixed, errors, area):
         out = ''
-        if self.fitfunc == 'poly':
+        if self.fitfunc == "poly" or self.fitfunc == "lpoly":
             c = 0
             for i in range(len(pars)):
                 fix = ""
                 if len(fixed) and fixed[i]: fix = "(fixed)"
-                if errors :
-                    out += '  p%d%s= %3.6f (%1.6f),' % (c,fix,pars[i], errors[i])
-                else:
-                    out += '  p%d%s= %3.6f,' % (c,fix,pars[i])
+                out += "  p%d%s= %3.6f" % (c, fix, pars[i])
+                if errors : out += " (%1.6f)" % errors[i]
+                out += ","
                 c+=1
             out = out[:-1]  # remove trailing ','
-        elif self.fitfunc == 'gauss' or self.fitfunc == 'lorentz':
+        else:
             i = 0
             c = 0
-            aunit = ''
-            ounit = ''
+            if self.fitfunc == "gauss" or self.fitfunc == "lorentz":
+                pnam = ["peak", "centre", "FWHM"]
+            elif self.fitfunc == "sinusoid":
+                pnam = ["amplitude", "period", "x0"]
+            aunit = ""
+            ounit = ""
             if self.data:
                 aunit = self.data.get_unit()
                 ounit = self.data.get_fluxunit()
             while i < len(pars):
-                if len(area):
-                    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)
-                else:
-                    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)
+                fix0 = fix1 = fix2 = ""
+                if i < len(fixed)-2:
+                    if fixed[i]:   fix0 = "(fixed)"
+                    if fixed[i+1]: fix1 = "(fixed)"
+                    if fixed[i+2]: fix2 = "(fixed)"
+                out += "  %2d: " % c
+                out += "%s%s = %3.3f %s, " % (pnam[0], fix0, pars[i],   ounit)
+                out += "%s%s = %3.3f %s, " % (pnam[1], fix1, pars[i+1], aunit)
+                out += "%s%s = %3.3f %s\n" % (pnam[2], fix2, pars[i+2], aunit)
+                if len(area): out += "      area = %3.3f %s %s\n" % (area[i], ounit, aunit)
                 c+=1
                 i+=3

trunk/python/scantable.py

@@ -2070 +2070 @@
 
     @asaplog_post_dec
+    def sinusoid_baseline(self, insitu=None, mask=None, minnwave=None, maxnwave=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
+            minnwave:   Minimum wave number in spectral window (default is 0)
+            maxnwave:   Maximum wave number in spectral window (default is 3)
+            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 combination of sinusoidal curves having
+            # wave numbers in spectral window from 1 to 10, 
+            # also with 3-sigma clipping, iteration up to 4 times
+            bscan = scan.sinusoid_baseline(maxnwave=10,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 minnwave is None: minnwave = 0
+        if maxnwave is None: maxnwave = 3
+        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._sinusoid_baseline(mask, minnwave, maxnwave, clipthresh, clipniter, outlog, blfile)
+            
+            workscan._add_history("sinusoid_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_sinusoid_baseline(self, insitu=None, mask=None, minnwave=None, maxnwave=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
+            minnwave:   Minimum wave number in spectral window (default is 0)
+            maxnwave:   Maximum wave number in spectral window (default is 3)
+            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_sinusoid_baseline(maxnwave=10, 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 minnwave is None: minnwave = 0
+        if maxnwave is None: maxnwave = 3
+        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_sinusoid_baseline(mask, minnwave, maxnwave, clipthresh, clipniter, curedge, threshold, chan_avg_limit, outlog, blfile)
+
+            workscan._add_history("auto_sinusoid_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 cspline_baseline(self, insitu=None, mask=None, npiece=None, clipthresh=None, clipniter=None, plot=None, outlog=None, blfile=None):
         """\

trunk/src/STFitter.cpp

@@ -39 +39 @@
 #include <scimath/Functionals/Gaussian1D.h>
 #include "Lorentzian1D.h"
+#include <scimath/Functionals/Sinusoid1D.h>
 #include <scimath/Functionals/Polynomial.h>
 #include <scimath/Mathematics/AutoDiff.h>
@@ -149 +150 @@
       funccomponents_.push_back(3);
     }
-  } else if (expr == "poly") {
-    funcs_.resize(1);
-    funcnames_.clear();
-    funccomponents_.clear();
-    funcs_[0] = new Polynomial<Float>(ncomp);
-      funcnames_.push_back(expr);
-      funccomponents_.push_back(ncomp);
   } else if (expr == "lorentz") {
     if (ncomp < 1) throw (AipsError("Need at least one lorentzian to fit."));
@@ -166 +160 @@
       funccomponents_.push_back(3);
     }
+  } else if (expr == "sinusoid") {
+    if (ncomp < 1) throw (AipsError("Need at least one sinusoid to fit."));
+    funcs_.resize(ncomp);
+    funcnames_.clear();
+    funccomponents_.clear();
+    for (Int k=0; k<ncomp; ++k) {
+      funcs_[k] = new Sinusoid1D<Float>();
+      funcnames_.push_back(expr);
+      funccomponents_.push_back(3);
+    }
+  } else if (expr == "poly") {
+    funcs_.resize(1);
+    funcnames_.clear();
+    funccomponents_.clear();
+    funcs_[0] = new Polynomial<Float>(ncomp);
+      funcnames_.push_back(expr);
+      funccomponents_.push_back(ncomp);
   } else {
-    //cerr << " compiled functions not yet implemented" << endl;
     LogIO os( LogOrigin( "Fitter", "setExpression()", WHERE ) ) ;
     os << LogIO::WARN << " compiled functions not yet implemented" << LogIO::POST;
@@ -244 +254 @@
             }
         }
+    } else if (dynamic_cast<Lorentzian1D<Float>* >(funcs_[0]) != 0) {
+        uInt count = 0;
+        for (uInt j=0; j < funcs_.nelements(); ++j) {
+            for (uInt i=0; i < funcs_[j]->nparameters(); ++i) {
+                (funcs_[j]->parameters())[i] = tmppar[count];
+                parameters_[count] = tmppar[count];
+                ++count;
+            }
+        }
+    } else if (dynamic_cast<Sinusoid1D<Float>* >(funcs_[0]) != 0) {
+        uInt count = 0;
+        for (uInt j=0; j < funcs_.nelements(); ++j) {
+            for (uInt i=0; i < funcs_[j]->nparameters(); ++i) {
+                (funcs_[j]->parameters())[i] = tmppar[count];
+                parameters_[count] = tmppar[count];
+                ++count;
+            }
+        }
     } else if (dynamic_cast<Polynomial<Float>* >(funcs_[0]) != 0) {
         for (uInt i=0; i < funcs_[0]->nparameters(); ++i) {
             parameters_[i] = tmppar[i];
             (funcs_[0]->parameters())[i] =  tmppar[i];
-        }
-    } else if (dynamic_cast<Lorentzian1D<Float>* >(funcs_[0]) != 0) {
-        uInt count = 0;
-        for (uInt j=0; j < funcs_.nelements(); ++j) {
-            for (uInt i=0; i < funcs_[j]->nparameters(); ++i) {
-                (funcs_[j]->parameters())[i] = tmppar[count];
-                parameters_[count] = tmppar[count];
-                ++count;
-            }
         }
     }
@@ -286 +305 @@
             }
         }
+    } else if (dynamic_cast<Lorentzian1D<Float>* >(funcs_[0]) != 0) {
+      uInt count = 0;
+        for (uInt j=0; j < funcs_.nelements(); ++j) {
+            for (uInt i=0; i < funcs_[j]->nparameters(); ++i) {
+                funcs_[j]->mask(i) = !fixed[count];
+                fixedpar_[count] = fixed[count];
+                ++count;
+            }
+        }
+    } else if (dynamic_cast<Sinusoid1D<Float>* >(funcs_[0]) != 0) {
+      uInt count = 0;
+        for (uInt j=0; j < funcs_.nelements(); ++j) {
+            for (uInt i=0; i < funcs_[j]->nparameters(); ++i) {
+                funcs_[j]->mask(i) = !fixed[count];
+                fixedpar_[count] = fixed[count];
+                ++count;
+            }
+        }
     } else if (dynamic_cast<Polynomial<Float>* >(funcs_[0]) != 0) {
         for (uInt i=0; i < funcs_[0]->nparameters(); ++i) {
             fixedpar_[i] = fixed[i];
             funcs_[0]->mask(i) =  !fixed[i];
-        }
-    } else if (dynamic_cast<Lorentzian1D<Float>* >(funcs_[0]) != 0) {
-      uInt count = 0;
-        for (uInt j=0; j < funcs_.nelements(); ++j) {
-            for (uInt i=0; i < funcs_[j]->nparameters(); ++i) {
-                funcs_[j]->mask(i) = !fixed[count];
-                fixedpar_[count] = fixed[count];
-                ++count;
-            }
         }
     }

trunk/src/Scantable.cpp

@@ -882 +882 @@
     uInt row = uInt(whichrow);
     stpol->setSpectra(getPolMatrix(row));
-    Float fang,fhand,parang;
+    Float fang,fhand;
     fang = focusTable_.getTotalAngle(mfocusidCol_(row));
     fhand = focusTable_.getFeedHand(mfocusidCol_(row));
@@ -1749 +1749 @@
 }
 
-bool Scantable::getFlagtraFast(int whichrow)
+bool Scantable::getFlagtraFast(uInt whichrow)
 {
   uChar flag;
   Vector<uChar> flags;
-  flagsCol_.get(uInt(whichrow), flags);
+  flagsCol_.get(whichrow, flags);
   flag = flags[0];
-  for (int i = 1; i < flags.size(); ++i) {
+  for (uInt i = 1; i < flags.size(); ++i) {
     flag &= flags[i];
   }
@@ -1761 +1761 @@
 }
 
-void Scantable::cubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, bool outLogger, const std::string& blfile) {
+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 hasSameNchan = true;
   bool outTextFile = false;
 
@@ -1777 +1777 @@
     if (coordInfo == "") coordInfo = "channel";
     hasSameNchan = hasSameNchanOverIFs();
-    firstIF = getIF(0);
-  }
-
-  //Fitter fitter = Fitter();
-  //fitter.setExpression("cspline", nPiece, thresClip, nIterClip);
+  }
+
+  Fitter fitter = Fitter();
+  fitter.setExpression("poly", order);
 
   int nRow = nrow();
@@ -1787 +1786 @@
 
   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;
-      }
-    }
-
+    fitBaseline(chanMask, whichrow, fitter);
+    setSpectrum(fitter.getResidual(), whichrow);
+    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "polyBaseline()", fitter);
   }
 
@@ -1816 +1795 @@
 }
 
-
-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)
+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 hasSameNchan = true;
   bool outTextFile = false;
 
@@ -1834 +1811 @@
     if (coordInfo == "") coordInfo = "channel";
     hasSameNchan = hasSameNchanOverIFs();
-    firstIF = getIF(0);
-  }
-
-  //Fitter fitter = Fitter();
-  //fitter.setExpression("cspline", nPiece, thresClip, nIterClip);
+  }
+
+  Fitter fitter = Fitter();
+  fitter.setExpression("poly", order);
 
   int nRow = nrow();
@@ -1871 +1847 @@
     //-------------------------------------------------------
 
-
-    //fitBaseline(chanMask, whichrow, fitter);
+    fitBaseline(chanMask, whichrow, fitter);
+    setSpectrum(fitter.getResidual(), whichrow);
+
+    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoPolyBaseline()", fitter);
+  }
+
+  if (outTextFile) ofs.close();
+}
+
+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 = true;
+  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();
+  }
+
+  //Fitter fitter = Fitter();
+  //fitter.setExpression("cspline", nPiece);
+
+  int nRow = nrow();
+  std::vector<bool> chanMask;
+
+  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
+    chanMask = getCompositeChanMask(whichrow, mask);
+    //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
     //setSpectrum(fitter.getResidual(), whichrow);
     std::vector<int> pieceRanges;
@@ -1878 +1887 @@
     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;
-      }
-    }
-
+    //
+
+    std::vector<bool>  fixed;
+    fixed.clear();
+    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()", pieceRanges, params, fixed);
   }
 
@@ -1899 +1897 @@
 }
 
+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 = true;
+  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();
+  }
+
+  //Fitter fitter = Fitter();
+  //fitter.setExpression("cspline", nPiece);
+
+  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, thresClip, nIterClip);
+    //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);
+    //
+
+    std::vector<bool> fixed;
+    fixed.clear();
+    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()", pieceRanges, params, fixed);
+  }
+
+  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) {
@@ -1921 +1987 @@
   int nElement = (int)(floor(floor((double)(nData/nPiece))+0.5));
 
-  std::vector<double> sectionList0, sectionList1, sectionListR;
-  sectionList0.push_back((double)x[0]);
+  std::vector<int> sectionList0, sectionList1;
+  std::vector<double> sectionListR;
+  sectionList0.push_back(x[0]);
   sectionRanges.clear();
   sectionRanges.push_back(x[0]);
   for (int i = 1; i < nPiece; ++i) {
-    double valX = (double)x[nElement*i];
+    int valX = 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));
+    sectionListR.push_back(1.0/(double)valX);
+
+    sectionRanges.push_back(valX-1);
+    sectionRanges.push_back(valX);
+  }
+  sectionList1.push_back(x[x.size()-1]+1);
   sectionRanges.push_back(x[x.size()-1]);
 
@@ -2111 +2178 @@
 }
 
-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)
-{
+  void Scantable::sinusoidBaseline(const std::vector<bool>& mask, int nMinWavesInSW, int nMaxWavesInSW, float thresClip, int nIterClip, bool outLogger, const std::string& blfile) {
   ofstream ofs;
   String coordInfo;
-  bool hasSameNchan;
-  int firstIF;
+  bool hasSameNchan = true;
   bool outTextFile = false;
 
@@ -2181 +2193 @@
     if (coordInfo == "") coordInfo = "channel";
     hasSameNchan = hasSameNchanOverIFs();
-    firstIF = getIF(0);
-  }
-
-  Fitter fitter = Fitter();
-  fitter.setExpression("poly", order);
+  }
+
+  //Fitter fitter = Fitter();
+  //fitter.setExpression("sinusoid", nMaxWavesInSW);
 
   int nRow = nrow();
@@ -2191 +2202 @@
 
   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;
-      }
-    }
-
+    //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
+    //setSpectrum(fitter.getResidual(), whichrow);
+    std::vector<int> pieceRanges;
+    std::vector<float> params;
+    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nMinWavesInSW, nMaxWavesInSW, params, thresClip, nIterClip);
+    setSpectrum(res, whichrow);
+    //
+
+    std::vector<bool>  fixed;
+    fixed.clear();
+    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "sinusoidBaseline()", pieceRanges, params, fixed);
   }
 
@@ -2216 +2219 @@
 }
 
-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)
+void 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)
 {
   ofstream ofs;
   String coordInfo;
-  bool hasSameNchan;
-  int firstIF;
+  bool hasSameNchan = true;
   bool outTextFile = false;
 
@@ -2233 +2235 @@
     if (coordInfo == "") coordInfo = "channel";
     hasSameNchan = hasSameNchanOverIFs();
-    firstIF = getIF(0);
-  }
-
-  Fitter fitter = Fitter();
-  fitter.setExpression("poly", order);
+  }
+
+  //Fitter fitter = Fitter();
+  //fitter.setExpression("sinusoid", nMaxWavesInSW);
 
   int nRow = nrow();
@@ -2271 +2272 @@
 
 
-    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;
-      }
-    }
-
+    //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
+    //setSpectrum(fitter.getResidual(), whichrow);
+    std::vector<int> pieceRanges;
+    std::vector<float> params;
+    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nMinWavesInSW, nMaxWavesInSW, params, thresClip, nIterClip);
+    setSpectrum(res, whichrow);
+    //
+
+    std::vector<bool> fixed;
+    fixed.clear();
+    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoSinusoidBaseline()", pieceRanges, params, fixed);
   }
 
@@ -2294 +2288 @@
 }
 
+std::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) {
+  if (nMaxWavesInSW < 1) {
+    throw(AipsError("maximum wave number must be a positive value"));
+  }
+  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();
+
+  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 = nMaxWavesInSW + 1;  //number of parameters to solve.
+
+  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 i = 0; i < currentNData; ++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 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 sinusoids forming the 
+    //best-fit curves (a0,b0,a1,b1,...), namely, a* and b* are of sine and cosine functions, 
+    //respectively. 
+    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 i = 0; i < nChan; ++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 ((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 (uInt 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);
+  uInt chanMaskSize = chanMask.size();
+  if (chanMaskSize != inMask.size()) {
+    throw(AipsError("different mask sizes"));
+  }
+  for (uInt 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();
+}
+*/
+
+/* for poly. the variations of outputFittingResult() should be merged into one eventually (2011/3/10 WK)  */
+void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, Fitter& fitter) {
+  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);
+
+    if (outLogger) {
+      LogIO ols(LogOrigin("Scantable", funcName, WHERE));
+      ols << formatBaselineParams(params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
+    }
+    if (outTextFile) {
+      ofs << formatBaselineParams(params, fixed, rms, masklist, whichrow, true) << flush;
+    }
+  }
+}
+
+/* for cspline. will be merged once cspline is available in fitter (2011/3/10 WK) */
+void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, const std::vector<int>& pieceEdges, const std::vector<float>& params, const std::vector<bool>& fixed) {
+  if (outLogger || outTextFile) {
+    float rms = getRms(chanMask, whichrow);
+    String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
+
+    if (outLogger) {
+      LogIO ols(LogOrigin("Scantable", funcName, WHERE));
+      ols << formatPiecewiseBaselineParams(pieceEdges, params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
+    }
+    if (outTextFile) {
+      ofs << formatPiecewiseBaselineParams(pieceEdges, params, fixed, rms, masklist, whichrow, true) << flush;
+    }
+  }
+}
+
+/* for sinusoid. will be merged once sinusoid is available in fitter (2011/3/10 WK) */
+void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, const std::vector<float>& params, const std::vector<bool>& fixed) {
+  if (outLogger || outTextFile) {
+    float rms = getRms(chanMask, whichrow);
+    String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
+
+    if (outLogger) {
+      LogIO ols(LogOrigin("Scantable", funcName, WHERE));
+      ols << formatBaselineParams(params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
+    }
+    if (outTextFile) {
+      ofs << formatBaselineParams(params, fixed, rms, masklist, whichrow, true) << flush;
+    }
+  }
+}
 
 float Scantable::getRms(const std::vector<bool>& mask, int whichrow) {
@@ -2302 +2566 @@
   float smean = 0.0;
   int n = 0;
-  for (int i = 0; i < spec.nelements(); ++i) {
+  for (uInt i = 0; i < spec.nelements(); ++i) {
     if (mask[i]) {
       mean += spec[i];
@@ -2344 +2608 @@
 
   if (verbose) {
-    oss << endl;
     oss << "Results of baseline fit" << endl;
     oss << "  rms = " << setprecision(6) << rms << endl;
@@ -2352 +2615 @@
 }
 
-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
+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 << ",";
+      }
+      std::string fix = ((fixed.size() > 0) && (fixed[i]) && verbose) ? "(fixed)" : "";
+      float vParam = params[i];
+      std::string equals = "= " + (String)((vParam < 0.0) ? "" : " ");
+      oss << "  p" << i << fix << equals << setprecision(6) << vParam;
+    }
+  } else {
+    oss << "  Not fitted";
+  }
+  oss << endl;
+
+  oss << formatBaselineParamsFooter(rms, verbose);
+  oss << flush;
+
+  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);
+      uInt nParam = params.size() / (ranges.size() / 2);
+      stringstream ss;
+      ss << ranges[ranges.size()-1] << flush;
+      int wRange = ss.str().size()*2+5;
+      ss.str("");
       for (uInt j = 0; j < ranges.size(); j+=2) {
-        oss << "[" << ranges[j] << "," << ranges[j+1] << "]";
+        ss << "  [" << ranges[j] << "," << ranges[j+1] << "]" << flush;
+        oss << setw(wRange) << left << ss.str();
+        ss.str("");
         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];
+          std::string fix = ((fixed.size() > 0) && (fixed[i]) && verbose) ? "(fixed)" : "";
+          float vParam = params[j/2*nParam+i];
+          std::string equals = "= " + (String)((vParam < 0.0) ? "" : " ");
+          oss << "  p" << i << fix << equals << setprecision(6) << vParam;
         }
+        oss << endl;
       }
     } else {
-      oss << "  ";
+      oss << "  " << endl;
     }
   } else {
-    oss << "  Not fitted";
+    oss << "  Not fitted" << endl;
   }
 
   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);
@@ -2471 +2699 @@
 }
 
-std::vector<int> Scantable::getMaskEdgeIndices(const std::vector<bool>& mask, bool getStartIndices)
+std::string Scantable::getMaskRangeList(const std::vector<bool>& mask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, bool verbose)
 {
   if (mask.size() < 2) {
     throw(AipsError("The mask elements should be > 1"));
   }
-  if (mask.size() != nchan()) {
+  int IF = getIF(whichrow);
+  if (mask.size() != (uInt)nchan(IF)) {
     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);
-    }
+  if (verbose) {
+    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=" << IF;
+    }
+    logOs << LogIO::POST;
+  }
+
+  std::vector<double> abcissa = getAbcissa(whichrow);
+  std::vector<int> edge = getMaskEdgeIndices(mask);
+
+  ostringstream oss;
+  oss.setf(ios::fixed);
+  oss << setprecision(1) << "[";
+  for (uInt i = 0; i < edge.size(); i+=2) {
+    if (i > 0) oss << ",";
+    oss << "[" << (float)abcissa[edge[i]] << "," << (float)abcissa[edge[i+1]] << "]";
+  }
+  oss << "]" << flush;
+
+  return String(oss);
+}
+
+std::vector<int> Scantable::getMaskEdgeIndices(const std::vector<bool>& mask)
+{
+  if (mask.size() < 2) {
+    throw(AipsError("The mask elements should be > 1"));
+  }
+
+  std::vector<int> out, startIndices, endIndices;
+  int maskSize = mask.size();
+
+  startIndices.clear();
+  endIndices.clear();
+
+  if (mask[0]) {
+    startIndices.push_back(0);
+  }
+  for (int i = 1; i < maskSize; ++i) {
+    if ((!mask[i-1]) && mask[i]) {
+      startIndices.push_back(i);
+    } else if (mask[i-1] && (!mask[i])) {
+      endIndices.push_back(i-1);
+    }
+  }
+  if (mask[maskSize-1]) {
+    endIndices.push_back(maskSize-1);
+  }
+
+  if (startIndices.size() != endIndices.size()) {
+    throw(AipsError("Inconsistent Mask Size: bad data?"));
+  }
+  for (uInt i = 0; i < startIndices.size(); ++i) {
+    if (startIndices[i] > endIndices[i]) {
+      throw(AipsError("Mask start index > mask end index"));
+    }
+  }
+
+  out.clear();
+  for (uInt i = 0; i < startIndices.size(); ++i) {
+    out.push_back(startIndices[i]);
+    out.push_back(endIndices[i]);
   }
 

trunk/src/Scantable.h

@@ -493 +493 @@
   void regridChannel( int nchan, double dnu, int irow ) ;
 
-  bool getFlagtraFast(int whichrow);
+  bool getFlagtraFast(casa::uInt whichrow);
 
   void polyBaseline(const std::vector<bool>& mask, 
@@ -521 +521 @@
                                bool outLogger=false, 
                                const std::string& blfile="");
+  void sinusoidBaseline(const std::vector<bool>& mask,
+                        int nMinWavesInSW,
+                        int nMaxWavesInSW,
+                        float thresClip, 
+                        int nIterClip,
+                        bool outLogger=false,
+                        const std::string& blfile="");
+  void autoSinusoidBaseline(const std::vector<bool>& mask,
+                            int nMinWavesInSW,
+                            int nMaxWavesInSW,
+                            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, 
@@ -657 +674 @@
                                           float thresClip=3.0, 
                                           int nIterClip=1);
+  std::vector<float> doSinusoidFitting(const std::vector<float>& data, 
+                                       const std::vector<bool>& mask,
+                                       int nMinWavesInSW,
+                                       int nMaxWavesInSW,
+                                       std::vector<float>& params,
+                                       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);
+                                bool hasSameNchan,
+                                bool verbose=false);
+  std::vector<int> getMaskEdgeIndices(const std::vector<bool>& mask);
   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 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);
+  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);
+  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);
 
   void applyChanFlag( casa::uInt whichrow, const std::vector<bool>& msk, casa::uChar flagval);

trunk/src/ScantableWrapper.h

@@ -269 +269 @@
   { table_->autoCubicSplineBaseline(mask, npiece, clipthresh, clipniter, edge, threshold, chan_avg_limit, outlog, blfile); }
 
+  void sinusoidBaseline(const std::vector<bool>& mask, int minnwave, int maxnwave, float clipthresh, int clipniter, bool outlog=false, const std::string& blfile="")
+  { table_->sinusoidBaseline(mask, minnwave, maxnwave, clipthresh, clipniter, outlog, blfile); }
+
+  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="")
+  { table_->autoSinusoidBaseline(mask, minnwave, maxnwave, clipthresh, clipniter, edge, threshold, chan_avg_limit, outlog, blfile); }
+
   float getRms(const std::vector<bool>& mask, int whichrow)
   { return table_->getRms(mask, whichrow); }
@@ -279 +285 @@
 
   bool getFlagtraFast(int whichrow=0) const
-  { return table_->getFlagtraFast(whichrow); }
+  { return table_->getFlagtraFast(casa::uInt(whichrow)); }
 
 

trunk/src/python_Scantable.cpp

@@ -146 +146 @@
     .def("_cspline_baseline", &ScantableWrapper::cubicSplineBaseline)
     .def("_auto_cspline_baseline", &ScantableWrapper::autoCubicSplineBaseline)
+    .def("_sinusoid_baseline", &ScantableWrapper::sinusoidBaseline)
+    .def("_auto_sinusoid_baseline", &ScantableWrapper::autoSinusoidBaseline)
     .def("get_rms", &ScantableWrapper::getRms)
     .def("format_blparams_row", &ScantableWrapper::formatBaselineParams)