Changeset 2081

Timestamp:

03/30/11 21:30:28 (13 years ago)

Author:

WataruKawasaki

Message:

New Development: No

JIRA Issue: Yes CAS-2847

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): Scantable

Description: Scantable::sinusoidBaseline(), Scantable::autoSinusoidBaseline()

Location:

trunk

Files:

• python/scantable.py (modified) (10 diffs)
• src/STLineFinder.cpp (modified) (1 diff)
• src/Scantable.cpp (modified) (36 diffs)
• src/Scantable.h (modified) (5 diffs)
• src/ScantableWrapper.h (modified) (1 diff)

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)
+    def sinusoid_baseline(self, insitu=None, mask=None, nwave=None, maxwavelength=None,
+                          clipthresh=None, clipniter=None, plot=None, getresidual=None, outlog=None, blfile=None):
+        """\
+        Return a scan which has been baselined (all rows) by sinusoidal functions.
+        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
+            nwave:         the maximum wave number of sinusoids within
+                           maxwavelength * (spectral range).
+                           The default is 3 (i.e., sinusoids with wave 
+                           number of 0(=constant), 1, 2, and 3 are 
+                           used for fitting). Also it is possible to
+                           explicitly specify all the wave numbers to
+                           be used, by giving a list including them 
+                           (e.g. [0,1,2,15,16]).
+            maxwavelength: the longest sinusoidal wavelength. The
+                           default is 1.0 (unit: spectral range). 
+            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'
+            getresidual:   if False, returns best-fit values instead of
+                           residual. (default is True)
+            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, 
+            # wave numbers in spectral window up to 10, 
             # also with 3-sigma clipping, iteration up to 4 times
-            bscan = scan.sinusoid_baseline(maxnwave=10,clipthresh=3.0,clipniter=4)
+            bscan = scan.sinusoid_baseline(nwave=10,clipthresh=3.0,clipniter=4)
+
+        Note:
+            The best-fit parameter values output in logger and/or blfile are now
+            based on specunit of 'channel'. 
         """
         
@@ -2110 +2124 @@
         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 = ""
-
+        if mask          is None: mask          = [True for i in xrange(nchan)]
+        if nwave         is None: nwave         = 3
+        if maxwavelength is None: maxwavelength = 1.0
+        if clipthresh    is None: clipthresh    = 3.0
+        if clipniter     is None: clipniter     = 1
+        if plot          is None: plot          = False
+        if getresidual   is None: getresidual   = True
+        if outlog        is None: outlog        = False
+        if blfile        is None: blfile        = ""
+
+        if isinstance(nwave, int):
+            in_nwave = nwave
+            nwave = []
+            for i in xrange(in_nwave+1): nwave.append(i)
+        
         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)
+            #CURRENTLY, PLOT=true is UNAVAILABLE UNTIL sinusoidal fitting is implemented as a fitter method. 
+            workscan._sinusoid_baseline(mask, nwave, maxwavelength, clipthresh, clipniter, getresidual, outlog, blfile)
             
             workscan._add_history("sinusoid_baseline", varlist)
@@ -2142 +2162 @@
 
 
-    def auto_sinusoid_baseline(self, insitu=None, mask=None, minnwave=None, maxnwave=None,
+    def auto_sinusoid_baseline(self, insitu=None, mask=None, nwave=None, maxwavelength=None,
                                clipthresh=None, clipniter=None, edge=None, threshold=None,
-                               chan_avg_limit=None, plot=None, outlog=None, blfile=None):
+                               chan_avg_limit=None, plot=None, getresidual=None, outlog=None, blfile=None):
         """\
         Return a scan which has been baselined (all rows) by cubic spline
@@ -2152 +2172 @@
 
         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)
+            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
+            nwave:         the maximum wave number of sinusoids within
+                           maxwavelength * (spectral range).
+                           The default is 3 (i.e., sinusoids with wave 
+                           number of 0(=constant), 1, 2, and 3 are 
+                           used for fitting). Also it is possible to
+                           explicitly specify all the wave numbers to
+                           be used, by giving a list including them 
+                           (e.g. [0,1,2,15,16]).
+            maxwavelength: the longest sinusoidal wavelength. The
+                           default is 1.0 (unit: spectral range). 
+            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'
+            getresidual:   if False, returns best-fit values instead of
+                           residual. (default is True)
+            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)
+            bscan = scan.auto_sinusoid_baseline(nwave=10, insitu=False)
+        
+        Note:
+            The best-fit parameter values output in logger and/or blfile are now
+            based on specunit of 'channel'. 
         """
 
@@ -2205 +2238 @@
         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 mask           is None: mask           = [True for i in xrange(nchan)]
+        if nwave          is None: nwave          = 3
+        if maxwavelength  is None: maxwavelength  = 1.0
+        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 = ""
+        if plot           is None: plot           = False
+        if getresidual    is None: getresidual    = True
+        if outlog         is None: outlog         = False
+        if blfile         is None: blfile         = ""
 
         outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
@@ -2221 +2255 @@
         from asap.asaplinefind import linefinder
         from asap import _is_sequence_or_number as _is_valid
+
+        if isinstance(nwave, int):
+            in_nwave = nwave
+            nwave = []
+            for i in xrange(in_nwave+1): nwave.append(i)
 
         if not (isinstance(edge, list) or isinstance(edge, tuple)): edge = [ edge ]
@@ -2244 +2283 @@
 
         try:
-            #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method. 
+            #CURRENTLY, PLOT=true is UNAVAILABLE UNTIL sinusoidal fitting is implemented as a fitter method. 
             if individualedge:
                 curedge = []
@@ -2250 +2289 @@
                     curedge += edge[i]
                 
-            workscan._auto_sinusoid_baseline(mask, minnwave, maxnwave, clipthresh, clipniter, curedge, threshold, chan_avg_limit, outlog, blfile)
+            workscan._auto_sinusoid_baseline(mask, nwave, maxwavelength, clipthresh, clipniter, curedge, threshold, chan_avg_limit, getresidual, outlog, blfile)
 
             workscan._add_history("auto_sinusoid_baseline", varlist)
@@ -2295 +2334 @@
             # also with 3-sigma clipping, iteration up to 4 times
             bscan = scan.cspline_baseline(npiece=2,clipthresh=3.0,clipniter=4)
+        
+        Note:
+            The best-fit parameter values output in logger and/or blfile are now
+            based on specunit of 'channel'. 
         """
         
@@ -2388 +2431 @@
         Example:
             bscan = scan.auto_cspline_baseline(npiece=3, insitu=False)
+        
+        Note:
+            The best-fit parameter values output in logger and/or blfile are now
+            based on specunit of 'channel'. 
         """
 

trunk/src/STLineFinder.cpp

@@ -1203 +1203 @@
     vel = scan->getAbcissa(last_row_used);
   } else {
-    for (int i = 0; i < spectrum.nelements(); ++i)
+    for (uInt i = 0; i < spectrum.nelements(); ++i)
       vel.push_back((double)i);
   }

trunk/src/Scantable.cpp

@@ -1778 +1778 @@
 {
   ofstream ofs;
-  String coordInfo;
+  String coordInfo = "";
   bool hasSameNchan = true;
   bool outTextFile = false;
@@ -1795 +1795 @@
   Fitter fitter = Fitter();
   fitter.setExpression("poly", order);
+  //fitter.setIterClipping(thresClip, nIterClip);
 
   int nRow = nrow();
@@ -1812 +1813 @@
 {
   ofstream ofs;
-  String coordInfo;
+  String coordInfo = "";
   bool hasSameNchan = true;
   bool outTextFile = false;
@@ -1829 +1830 @@
   Fitter fitter = Fitter();
   fitter.setExpression("poly", order);
+  //fitter.setIterClipping(thresClip, nIterClip);
 
   int nRow = nrow();
@@ -1870 +1872 @@
 }
 
-void Scantable::cubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, bool outLogger, const std::string& blfile) {
+void Scantable::cubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, bool outLogger, const std::string& blfile)
+{
   ofstream ofs;
-  String coordInfo;
+  String coordInfo = "";
   bool hasSameNchan = true;
   bool outTextFile = false;
@@ -1889 +1892 @@
   //Fitter fitter = Fitter();
   //fitter.setExpression("cspline", nPiece);
+  //fitter.setIterClipping(thresClip, nIterClip);
 
   int nRow = nrow();
@@ -1895 +1899 @@
   for (int whichrow = 0; whichrow < nRow; ++whichrow) {
     chanMask = getCompositeChanMask(whichrow, mask);
-    //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
+    //fitBaseline(chanMask, whichrow, fitter);
     //setSpectrum(fitter.getResidual(), whichrow);
-    std::vector<int> pieceRanges;
+    std::vector<int> pieceEdges;
     std::vector<float> params;
-    std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, pieceRanges, params, nPiece, thresClip, nIterClip, true);
+    std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, nPiece, pieceEdges, params, thresClip, nIterClip, true);
     setSpectrum(res, whichrow);
     //
 
-    std::vector<bool>  fixed;
-    fixed.clear();
-    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()", pieceRanges, params, fixed);
+    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()", pieceEdges, params);
   }
 
@@ -1914 +1916 @@
 {
   ofstream ofs;
-  String coordInfo;
+  String coordInfo = "";
   bool hasSameNchan = true;
   bool outTextFile = false;
@@ -1931 +1933 @@
   //Fitter fitter = Fitter();
   //fitter.setExpression("cspline", nPiece);
+  //fitter.setIterClipping(thresClip, nIterClip);
 
   int nRow = nrow();
@@ -1964 +1967 @@
 
 
-    //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
+    //fitBaseline(chanMask, whichrow, fitter);
     //setSpectrum(fitter.getResidual(), whichrow);
-    std::vector<int> pieceRanges;
+    std::vector<int> pieceEdges;
     std::vector<float> params;
-    std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, pieceRanges, params, nPiece, thresClip, nIterClip, true);
+    std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, nPiece, pieceEdges, params, thresClip, nIterClip, true);
     setSpectrum(res, whichrow);
     //
 
-    std::vector<bool> fixed;
-    fixed.clear();
-    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()", pieceRanges, params, fixed);
+    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()", pieceEdges, params);
   }
 
@@ -1980 +1981 @@
 }
 
-std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data, const std::vector<bool>& mask, std::vector<int>& idxEdge, std::vector<float>& params, int nPiece, float thresClip, int nIterClip, bool getResidual) {
+std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data, const std::vector<bool>& mask, int nPiece, std::vector<int>& idxEdge, std::vector<float>& params, float thresClip, int nIterClip, bool getResidual)
+{
+  if (data.size() != mask.size()) {
+    throw(AipsError("data and mask sizes are not identical"));
+  }
   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"));
   }
 
@@ -1998 +2000 @@
   }
 
-  int nData = x.size();
-  int nElement = (int)(floor(floor((double)(nData/nPiece))+0.5));
+  int initNData = x.size();
+
+  int nElement = (int)(floor(floor((double)(initNData/nPiece))+0.5));
   std::vector<double> invEdge;
-
   idxEdge.clear();
   idxEdge.push_back(x[0]);
-
   for (int i = 1; i < nPiece; ++i) {
     int valX = x[nElement*i];
@@ -2010 +2011 @@
     invEdge.push_back(1.0/(double)valX);
   }
-
   idxEdge.push_back(x[x.size()-1]+1);
 
-  std::vector<double> x1, x2, x3, z1, x1z1, x2z1, x3z1, r1;
+  int nData = initNData;
+  int nDOF = nPiece + 3;  //number of parameters to solve, namely, 4+(nPiece-1).
+
+  std::vector<double> x1, x2, x3, z1, x1z1, x2z1, x3z1, r1, residual;
   for (int i = 0; i < nChan; ++i) {
     double di = (double)i;
@@ -2025 +2028 @@
     x3z1.push_back(dD*di*di*di);
     r1.push_back(0.0);
-  }
-
-  int currentNData = nData;
-  int nDOF = nPiece + 3;  //number of parameters to solve, namely, 4+(nPiece-1).
+    residual.push_back(0.0);
+  }
 
   for (int nClip = 0; nClip < nIterClip+1; ++nClip) {
@@ -2142 +2143 @@
       for (int i = idxEdge[n]; i < idxEdge[n+1]; ++i) {
         r1[i] = a0 + a1*x1[i] + a2*x2[i] + a3*x3[i];
+        residual[i] = z1[i] - r1[i];
       }
       params.push_back(a0);
@@ -2161 +2163 @@
       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 += residual[i]*residual[i]*(double)maskArray[i];
       }
       stdDev = sqrt(stdDev/(double)nData);
@@ -2173 +2172 @@
       int newNData = 0;
       for (int i = 0; i < nChan; ++i) {
-        if (rz[i] >= thres) {
+        if (abs(residual[i]) >= thres) {
           maskArray[i] = 0;
         }
@@ -2180 +2179 @@
         }
       }
-      if (newNData == currentNData) {
+      if (newNData == nData) {
         break; //no more flag to add. iteration stops.
       } else {
-        currentNData = newNData;
+        nData = newNData;
       }
     }
@@ -2191 +2190 @@
   if (getResidual) {
     for (int i = 0; i < nChan; ++i) {
-      result.push_back((float)(z1[i] - r1[i]));
+      result.push_back((float)residual[i]);
     }
   } else {
@@ -2202 +2201 @@
 }
 
-void Scantable::sinusoidBaseline(const std::vector<bool>& mask, int nMinWavesInSW, int nMaxWavesInSW, float thresClip, int nIterClip, bool outLogger, const std::string& blfile) {
+void Scantable::sinusoidBaseline(const std::vector<bool>& mask, const std::vector<int>& nWaves, float maxWaveLength, float thresClip, int nIterClip, bool getResidual, bool outLogger, const std::string& blfile)
+{
   ofstream ofs;
-  String coordInfo;
+  String coordInfo = "";
   bool hasSameNchan = true;
   bool outTextFile = false;
@@ -2220 +2220 @@
 
   //Fitter fitter = Fitter();
-  //fitter.setExpression("sinusoid", nMaxWavesInSW);
+  //fitter.setExpression("sinusoid", nWaves);
+  //fitter.setIterClipping(thresClip, nIterClip);
 
   int nRow = nrow();
@@ -2227 +2228 @@
   for (int whichrow = 0; whichrow < nRow; ++whichrow) {
     chanMask = getCompositeChanMask(whichrow, mask);
-    //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
+    //fitBaseline(chanMask, whichrow, fitter);
     //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, true);
+    //std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, params, thresClip, nIterClip, true);
+    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, maxWaveLength, params, thresClip, nIterClip, getResidual);
     setSpectrum(res, whichrow);
     //
 
-    std::vector<bool>  fixed;
-    fixed.clear();
-    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "sinusoidBaseline()", pieceRanges, params, fixed);
+    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "sinusoidBaseline()", params);
   }
 
@@ -2243 +2242 @@
 }
 
-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)
+void Scantable::autoSinusoidBaseline(const std::vector<bool>& mask, const std::vector<int>& nWaves, float maxWaveLength, float thresClip, int nIterClip, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool getResidual, bool outLogger, const std::string& blfile)
 {
   ofstream ofs;
-  String coordInfo;
+  String coordInfo = "";
   bool hasSameNchan = true;
   bool outTextFile = false;
@@ -2262 +2261 @@
 
   //Fitter fitter = Fitter();
-  //fitter.setExpression("sinusoid", nMaxWavesInSW);
+  //fitter.setExpression("sinusoid", nWaves);
+  //fitter.setIterClipping(thresClip, nIterClip);
 
   int nRow = nrow();
@@ -2296 +2296 @@
 
 
-    //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
+    //fitBaseline(chanMask, whichrow, fitter);
     //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, true);
+    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, maxWaveLength, params, thresClip, nIterClip, getResidual);
     setSpectrum(res, whichrow);
     //
 
-    std::vector<bool> fixed;
-    fixed.clear();
-    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoSinusoidBaseline()", pieceRanges, params, fixed);
+    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoSinusoidBaseline()", params);
   }
 
@@ -2312 +2309 @@
 }
 
-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, bool getResidual) {
+std::vector<float> Scantable::doSinusoidFitting(const std::vector<float>& data, const std::vector<bool>& mask, const std::vector<int>& waveNumbers, float maxWaveLength, std::vector<float>& params, float thresClip, int nIterClip, bool getResidual)
+{
   if (data.size() != mask.size()) {
-    throw(AipsError("data and mask have different size"));
-  }
-  if ((nMinWavesInSW < 0) || (nMaxWavesInSW < 0)) {
+    throw(AipsError("data and mask sizes are not identical"));
+  }
+  if (data.size() < 2) {
+    throw(AipsError("data size is too short"));
+  }
+  if (waveNumbers.size() == 0) {
+    throw(AipsError("missing wave number info"));
+  }
+  std::vector<int> nWaves;  // sorted and uniqued array of wave numbers
+  nWaves.reserve(waveNumbers.size());
+  copy(waveNumbers.begin(), waveNumbers.end(), back_inserter(nWaves));
+  sort(nWaves.begin(), nWaves.end());
+  std::vector<int>::iterator end_it = unique(nWaves.begin(), nWaves.end());
+  nWaves.erase(end_it, nWaves.end());
+
+  int minNWaves = nWaves[0];
+  if (minNWaves < 0) {
     throw(AipsError("wave number must be positive or zero (i.e. constant)"));
-  } else {
-    if (nMaxWavesInSW < nMinWavesInSW) {
-      int temp = nMaxWavesInSW;
-      nMaxWavesInSW = nMinWavesInSW;
-      nMinWavesInSW = temp;
-    }
-  }
+  }
+  bool hasConstantTerm = (minNWaves == 0);
 
   int nChan = data.size();
@@ -2336 +2343 @@
   }
 
-  int nData = x.size();
-
-  std::vector<double> x1, x2, x3, z1, x1z1, x2z1, x3z1, r1;
+  int initNData = x.size();
+
+  int nData = initNData;
+  int nDOF = nWaves.size() * 2 - (hasConstantTerm ? 1 : 0);  //number of parameters to solve.
+
+  const double PI = 6.0 * asin(0.5); // PI (= 3.141592653...)
+  double baseXFactor = 2.0*PI/(double)maxWaveLength/(double)(nChan-1);  //the denominator (nChan-1) should be changed to (xdata[nChan-1]-xdata[0]) for accepting x-values given in velocity or frequency when this function is moved to fitter.
+
+  // xArray : contains elemental values for computing the least-square matrix. 
+  //          xArray.size() is nDOF and xArray[*].size() is nChan. 
+  //          Each xArray element are as follows:
+  //          xArray[0]    = {1.0, 1.0, 1.0, ..., 1.0}, 
+  //          xArray[2n-1] = {sin(nPI/L*x[0]), sin(nPI/L*x[1]), ..., sin(nPI/L*x[nChan])}, 
+  //          xArray[2n]   = {cos(nPI/L*x[0]), cos(nPI/L*x[1]), ..., cos(nPI/L*x[nChan])}, 
+  //          where (1 <= n <= nMaxWavesInSW),
+  //          or, 
+  //          xArray[2n-1] = {sin(wn[n]PI/L*x[0]), sin(wn[n]PI/L*x[1]), ..., sin(wn[n]PI/L*x[nChan])}, 
+  //          xArray[2n]   = {cos(wn[n]PI/L*x[0]), cos(wn[n]PI/L*x[1]), ..., cos(wn[n]PI/L*x[nChan])}, 
+  //          where wn[n] denotes waveNumbers[n] (1 <= n <= waveNumbers.size()).
+  std::vector<std::vector<double> > xArray;
+  if (hasConstantTerm) {
+    std::vector<double> xu;
+    for (int j = 0; j < nChan; ++j) {
+      xu.push_back(1.0);
+    }
+    xArray.push_back(xu);
+  }
+  for (uInt i = (hasConstantTerm ? 1 : 0); i < nWaves.size(); ++i) {
+    double xFactor = baseXFactor*(double)nWaves[i];
+    std::vector<double> xs, xc;
+    xs.clear();
+    xc.clear();
+    for (int j = 0; j < nChan; ++j) {
+      xs.push_back(sin(xFactor*(double)j));
+      xc.push_back(cos(xFactor*(double)j));
+    }
+    xArray.push_back(xs);
+    xArray.push_back(xc);
+  }
+
+  std::vector<double> z1, r1, residual;
   for (int i = 0; i < nChan; ++i) {
-    double di = (double)i;
-    double dD = (double)data[i];
-    x1.push_back(di);
-    x2.push_back(di*di);
-    x3.push_back(di*di*di);
-    z1.push_back(dD);
-    x1z1.push_back(di*dD);
-    x2z1.push_back(di*di*dD);
-    x3z1.push_back(di*di*di*dD);
+    z1.push_back((double)data[i]);
     r1.push_back(0.0);
-  }
-
-  int currentNData = nData;
-  int nDOF = nMaxWavesInSW + 1;  //number of parameters to solve.
+    residual.push_back(0.0);
+  }
 
   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. 
+    // 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];
@@ -2369 +2405 @@
     }
 
-    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 k = 0; k < nChan; ++k) {
+      if (maskArray[k] == 0) continue;
+
+      for (int i = 0; i < nDOF; ++i) {
+        for (int j = i; j < nDOF; ++j) {
+          xMatrix[i][j] += xArray[i][k] * xArray[j][k];
+        }
+        zMatrix[i] += z1[k] * xArray[i][k];
+      }
     }
 
@@ -2425 +2454 @@
     }
     //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. 
+    //best-fit curves (a0,s1,c1,s2,c2,...), where a0 is constant and s* and c* are of sine 
+    //and cosine functions, respectively. 
     std::vector<double> y;
+    params.clear();
     for (int i = 0; i < nDOF; ++i) {
       y.push_back(0.0);
@@ -2433 +2463 @@
         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();
+      params.push_back(y[i]);
+    }
 
     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);
-
+      r1[i] = y[0];
+      for (int j = 1; j < nDOF; ++j) {
+        r1[i] += y[j]*xArray[j][i];
+      }
+      residual[i] = z1[i] - r1[i];
+    }
 
     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 += residual[i]*residual[i]*(double)maskArray[i];
       }
       stdDev = sqrt(stdDev/(double)nData);
@@ -2465 +2486 @@
       int newNData = 0;
       for (int i = 0; i < nChan; ++i) {
-        if (rz[i] >= thres) {
+        if (abs(residual[i]) >= thres) {
           maskArray[i] = 0;
         }
@@ -2472 +2493 @@
         }
       }
-      if (newNData == currentNData) {
-        break;                   //no additional flag. finish iteration.
+      if (newNData == nData) {
+        break; //no more flag to add. iteration stops.
       } else {
-        currentNData = newNData;
+        nData = newNData;
       }
     }
@@ -2483 +2504 @@
   if (getResidual) {
     for (int i = 0; i < nChan; ++i) {
-      result.push_back((float)(z1[i] - r1[i]));
+      result.push_back((float)residual[i]);
     }
   } else {
@@ -2496 +2517 @@
 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<double> dAbcissa = getAbcissa(whichrow);
+  std::vector<float> abcissa;
+  for (uInt i = 0; i < dAbcissa.size(); ++i) {
+    abcissa.push_back((float)dAbcissa[i]);
   }
   std::vector<float> spec = getSpectrum(whichrow);
 
-  fitter.setData(abcs, spec, mask);
+  fitter.setData(abcissa, spec, mask);
   fitter.lfit();
 }
@@ -2566 +2587 @@
 
 /* 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>& edge, const std::vector<float>& params, const std::vector<bool>& fixed) {
+void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, const std::vector<int>& edge, const std::vector<float>& params) {
   if (outLogger || outTextFile) {
     float rms = getRms(chanMask, whichrow);
     String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
+    std::vector<bool> fixed;
+    fixed.clear();
 
     if (outLogger) {
@@ -2582 +2605 @@
 
 /* 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) {
+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) {
   if (outLogger || outTextFile) {
     float rms = getRms(chanMask, whichrow);
     String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
+    std::vector<bool> fixed;
+    fixed.clear();
 
     if (outLogger) {

trunk/src/Scantable.h

@@ -522 +522 @@
                                const std::string& blfile="");
   void sinusoidBaseline(const std::vector<bool>& mask,
-                        int nMinWavesInSW,
-                        int nMaxWavesInSW,
+                        const std::vector<int>& nWaves,
+                        float maxWaveLength,
                         float thresClip, 
                         int nIterClip,
+                        bool getResidual=true,
                         bool outLogger=false,
                         const std::string& blfile="");
   void autoSinusoidBaseline(const std::vector<bool>& mask,
-                            int nMinWavesInSW,
-                            int nMaxWavesInSW,
+                            const std::vector<int>& nWaves,
+                            float maxWaveLength,
                             float thresClip, 
                             int nIterClip,
@@ -536 +537 @@
                             float threshold=3.0, 
                             int chanAvgLimit=1, 
+                            bool getResidual=true,
                             bool outLogger=false, 
                             const std::string& blfile="");
@@ -672 +674 @@
   std::vector<float> doCubicSplineFitting(const std::vector<float>& data, 
                                           const std::vector<bool>& mask,
+                                          int nPiece,
                                           std::vector<int>& idxEdge,
                                           std::vector<float>& params,
-                                          int nPiece=2,
                                           float thresClip=3.0, 
                                           int nIterClip=1,
@@ -680 +682 @@
   std::vector<float> doSinusoidFitting(const std::vector<float>& data, 
                                        const std::vector<bool>& mask,
-                                       int nMinWavesInSW,
-                                       int nMaxWavesInSW,
+                                       const std::vector<int>& waveNumbers,
+                                       float maxWaveLength, 
                                        std::vector<float>& params,
                                        float thresClip=3.0, 
@@ -698 +700 @@
   //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>& edge, 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 outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, std::ofstream& ofs, const casa::String& funcName, const std::vector<int>& edge, const std::vector<float>& params);
+  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);
 
   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); }
+  void sinusoidBaseline(const std::vector<bool>& mask, const std::vector<int>& nwave, float maxwavelength, float clipthresh, int clipniter, bool getresidual=true, bool outlog=false, const std::string& blfile="")
+  { table_->sinusoidBaseline(mask, nwave, maxwavelength, clipthresh, clipniter, getresidual, outlog, blfile); }
+
+  void autoSinusoidBaseline(const std::vector<bool>& mask, const std::vector<int>& nwave, float maxwavelength, float clipthresh, int clipniter, const std::vector<int>& edge, float threshold=5.0, int chan_avg_limit=1, bool getresidual=true, bool outlog=false, const std::string& blfile="")
+  { table_->autoSinusoidBaseline(mask, nwave, maxwavelength, clipthresh, clipniter, edge, threshold, chan_avg_limit, getresidual, outlog, blfile); }
 
   float getRms(const std::vector<bool>& mask, int whichrow)