Changeset 2081 for trunk


Ignore:
Timestamp:
03/30/11 21:30:28 (14 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:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/python/scantable.py

    r2047 r2081  
    20702070
    20712071    @asaplog_post_dec
    2072     def sinusoid_baseline(self, insitu=None, mask=None, minnwave=None, maxnwave=None,
    2073                           clipthresh=None, clipniter=None, plot=None, outlog=None, blfile=None):
    2074         """\
    2075         Return a scan which has been baselined (all rows) by cubic spline function (piecewise cubic polynomial).
    2076         Parameters:
    2077             insitu:     If False a new scantable is returned.
    2078                         Otherwise, the scaling is done in-situ
    2079                         The default is taken from .asaprc (False)
    2080             mask:       An optional mask
    2081             minnwave:   Minimum wave number in spectral window (default is 0)
    2082             maxnwave:   Maximum wave number in spectral window (default is 3)
    2083             clipthresh: Clipping threshold. (default is 3.0, unit: sigma)
    2084             clipniter:  maximum number of iteration of 'clipthresh'-sigma clipping (default is 1)
    2085             plot:   *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
    2086                         plot the fit and the residual. In this each
    2087                         indivual fit has to be approved, by typing 'y'
    2088                         or 'n'
    2089             outlog:     Output the coefficients of the best-fit
    2090                         function to logger (default is False)
    2091             blfile:     Name of a text file in which the best-fit
    2092                         parameter values to be written
    2093                         (default is "": no file/logger output)
     2072    def sinusoid_baseline(self, insitu=None, mask=None, nwave=None, maxwavelength=None,
     2073                          clipthresh=None, clipniter=None, plot=None, getresidual=None, outlog=None, blfile=None):
     2074        """\
     2075        Return a scan which has been baselined (all rows) by sinusoidal functions.
     2076        Parameters:
     2077            insitu:        If False a new scantable is returned.
     2078                           Otherwise, the scaling is done in-situ
     2079                           The default is taken from .asaprc (False)
     2080            mask:          An optional mask
     2081            nwave:         the maximum wave number of sinusoids within
     2082                           maxwavelength * (spectral range).
     2083                           The default is 3 (i.e., sinusoids with wave
     2084                           number of 0(=constant), 1, 2, and 3 are
     2085                           used for fitting). Also it is possible to
     2086                           explicitly specify all the wave numbers to
     2087                           be used, by giving a list including them
     2088                           (e.g. [0,1,2,15,16]).
     2089            maxwavelength: the longest sinusoidal wavelength. The
     2090                           default is 1.0 (unit: spectral range).
     2091            clipthresh:    Clipping threshold. (default is 3.0, unit: sigma)
     2092            clipniter:     maximum number of iteration of 'clipthresh'-sigma clipping (default is 1)
     2093            plot:      *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
     2094                           plot the fit and the residual. In this each
     2095                           indivual fit has to be approved, by typing 'y'
     2096                           or 'n'
     2097            getresidual:   if False, returns best-fit values instead of
     2098                           residual. (default is True)
     2099            outlog:        Output the coefficients of the best-fit
     2100                           function to logger (default is False)
     2101            blfile:        Name of a text file in which the best-fit
     2102                           parameter values to be written
     2103                           (default is "": no file/logger output)
    20942104
    20952105        Example:
    20962106            # return a scan baselined by a combination of sinusoidal curves having
    2097             # wave numbers in spectral window from 1 to 10,
     2107            # wave numbers in spectral window up to 10,
    20982108            # also with 3-sigma clipping, iteration up to 4 times
    2099             bscan = scan.sinusoid_baseline(maxnwave=10,clipthresh=3.0,clipniter=4)
     2109            bscan = scan.sinusoid_baseline(nwave=10,clipthresh=3.0,clipniter=4)
     2110
     2111        Note:
     2112            The best-fit parameter values output in logger and/or blfile are now
     2113            based on specunit of 'channel'.
    21002114        """
    21012115       
    … …  
    21102124        nchan = workscan.nchan()
    21112125       
    2112         if mask is None: mask = [True for i in xrange(nchan)]
    2113         if minnwave is None: minnwave = 0
    2114         if maxnwave is None: maxnwave = 3
    2115         if clipthresh is None: clipthresh = 3.0
    2116         if clipniter is None: clipniter = 1
    2117         if plot is None: plot = False
    2118         if outlog is None: outlog = False
    2119         if blfile is None: blfile = ""
    2120 
     2126        if mask          is None: mask          = [True for i in xrange(nchan)]
     2127        if nwave         is None: nwave         = 3
     2128        if maxwavelength is None: maxwavelength = 1.0
     2129        if clipthresh    is None: clipthresh    = 3.0
     2130        if clipniter     is None: clipniter     = 1
     2131        if plot          is None: plot          = False
     2132        if getresidual   is None: getresidual   = True
     2133        if outlog        is None: outlog        = False
     2134        if blfile        is None: blfile        = ""
     2135
     2136        if isinstance(nwave, int):
     2137            in_nwave = nwave
     2138            nwave = []
     2139            for i in xrange(in_nwave+1): nwave.append(i)
     2140       
    21212141        outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
    21222142       
    21232143        try:
    2124             #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method.
    2125             workscan._sinusoid_baseline(mask, minnwave, maxnwave, clipthresh, clipniter, outlog, blfile)
     2144            #CURRENTLY, PLOT=true is UNAVAILABLE UNTIL sinusoidal fitting is implemented as a fitter method.
     2145            workscan._sinusoid_baseline(mask, nwave, maxwavelength, clipthresh, clipniter, getresidual, outlog, blfile)
    21262146           
    21272147            workscan._add_history("sinusoid_baseline", varlist)
    … …  
    21422162
    21432163
    2144     def auto_sinusoid_baseline(self, insitu=None, mask=None, minnwave=None, maxnwave=None,
     2164    def auto_sinusoid_baseline(self, insitu=None, mask=None, nwave=None, maxwavelength=None,
    21452165                               clipthresh=None, clipniter=None, edge=None, threshold=None,
    2146                                chan_avg_limit=None, plot=None, outlog=None, blfile=None):
     2166                               chan_avg_limit=None, plot=None, getresidual=None, outlog=None, blfile=None):
    21472167        """\
    21482168        Return a scan which has been baselined (all rows) by cubic spline
    … …  
    21522172
    21532173        Parameters:
    2154             insitu:     if False a new scantable is returned.
    2155                         Otherwise, the scaling is done in-situ
    2156                         The default is taken from .asaprc (False)
    2157             mask:       an optional mask retreived from scantable
    2158             minnwave:   Minimum wave number in spectral window (default is 0)
    2159             maxnwave:   Maximum wave number in spectral window (default is 3)
    2160             clipthresh: Clipping threshold. (default is 3.0, unit: sigma)
    2161             clipniter:  maximum number of iteration of 'clipthresh'-sigma clipping (default is 1)
    2162             edge:       an optional number of channel to drop at
    2163                         the edge of spectrum. If only one value is
    2164                         specified, the same number will be dropped
    2165                         from both sides of the spectrum. Default
    2166                         is to keep all channels. Nested tuples
    2167                         represent individual edge selection for
    2168                         different IFs (a number of spectral channels
    2169                         can be different)
    2170             threshold:  the threshold used by line finder. It is
    2171                         better to keep it large as only strong lines
    2172                         affect the baseline solution.
    2173             chan_avg_limit:
    2174                         a maximum number of consequtive spectral
    2175                         channels to average during the search of
    2176                         weak and broad lines. The default is no
    2177                         averaging (and no search for weak lines).
    2178                         If such lines can affect the fitted baseline
    2179                         (e.g. a high order polynomial is fitted),
    2180                         increase this parameter (usually values up
    2181                         to 8 are reasonable). Most users of this
    2182                         method should find the default value sufficient.
    2183             plot:   *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
    2184                         plot the fit and the residual. In this each
    2185                         indivual fit has to be approved, by typing 'y'
    2186                         or 'n'
    2187             outlog:     Output the coefficients of the best-fit
    2188                         function to logger (default is False)
    2189             blfile:     Name of a text file in which the best-fit
    2190                         parameter values to be written
    2191                         (default is "": no file/logger output)
     2174            insitu:        if False a new scantable is returned.
     2175                           Otherwise, the scaling is done in-situ
     2176                           The default is taken from .asaprc (False)
     2177            mask:          an optional mask retreived from scantable
     2178            nwave:         the maximum wave number of sinusoids within
     2179                           maxwavelength * (spectral range).
     2180                           The default is 3 (i.e., sinusoids with wave
     2181                           number of 0(=constant), 1, 2, and 3 are
     2182                           used for fitting). Also it is possible to
     2183                           explicitly specify all the wave numbers to
     2184                           be used, by giving a list including them
     2185                           (e.g. [0,1,2,15,16]).
     2186            maxwavelength: the longest sinusoidal wavelength. The
     2187                           default is 1.0 (unit: spectral range).
     2188            clipthresh:    Clipping threshold. (default is 3.0, unit: sigma)
     2189            clipniter:     maximum number of iteration of 'clipthresh'-sigma clipping (default is 1)
     2190            edge:          an optional number of channel to drop at
     2191                           the edge of spectrum. If only one value is
     2192                           specified, the same number will be dropped
     2193                           from both sides of the spectrum. Default
     2194                           is to keep all channels. Nested tuples
     2195                           represent individual edge selection for
     2196                           different IFs (a number of spectral channels
     2197                           can be different)
     2198            threshold:     the threshold used by line finder. It is
     2199                           better to keep it large as only strong lines
     2200                           affect the baseline solution.
     2201            chan_avg_limit:a maximum number of consequtive spectral
     2202                           channels to average during the search of
     2203                           weak and broad lines. The default is no
     2204                           averaging (and no search for weak lines).
     2205                           If such lines can affect the fitted baseline
     2206                           (e.g. a high order polynomial is fitted),
     2207                           increase this parameter (usually values up
     2208                           to 8 are reasonable). Most users of this
     2209                           method should find the default value sufficient.
     2210            plot:      *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
     2211                           plot the fit and the residual. In this each
     2212                           indivual fit has to be approved, by typing 'y'
     2213                           or 'n'
     2214            getresidual:   if False, returns best-fit values instead of
     2215                           residual. (default is True)
     2216            outlog:        Output the coefficients of the best-fit
     2217                           function to logger (default is False)
     2218            blfile:        Name of a text file in which the best-fit
     2219                           parameter values to be written
     2220                           (default is "": no file/logger output)
    21922221
    21932222        Example:
    2194             bscan = scan.auto_sinusoid_baseline(maxnwave=10, insitu=False)
     2223            bscan = scan.auto_sinusoid_baseline(nwave=10, insitu=False)
     2224       
     2225        Note:
     2226            The best-fit parameter values output in logger and/or blfile are now
     2227            based on specunit of 'channel'.
    21952228        """
    21962229
    … …  
    22052238        nchan = workscan.nchan()
    22062239       
    2207         if mask is None: mask = [True for i in xrange(nchan)]
    2208         if minnwave is None: minnwave = 0
    2209         if maxnwave is None: maxnwave = 3
    2210         if clipthresh is None: clipthresh = 3.0
    2211         if clipniter is None: clipniter = 1
    2212         if edge is None: edge = (0, 0)
    2213         if threshold is None: threshold = 3
     2240        if mask           is None: mask          = [True for i in xrange(nchan)]
     2241        if nwave          is None: nwave          = 3
     2242        if maxwavelength  is None: maxwavelength  = 1.0
     2243        if clipthresh     is None: clipthresh    = 3.0
     2244        if clipniter      is None: clipniter      = 1
     2245        if edge           is None: edge           = (0,0)
     2246        if threshold      is None: threshold      = 3
    22142247        if chan_avg_limit is None: chan_avg_limit = 1
    2215         if plot is None: plot = False
    2216         if outlog is None: outlog = False
    2217         if blfile is None: blfile = ""
     2248        if plot           is None: plot           = False
     2249        if getresidual    is None: getresidual    = True
     2250        if outlog         is None: outlog         = False
     2251        if blfile         is None: blfile         = ""
    22182252
    22192253        outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
    … …  
    22212255        from asap.asaplinefind import linefinder
    22222256        from asap import _is_sequence_or_number as _is_valid
     2257
     2258        if isinstance(nwave, int):
     2259            in_nwave = nwave
     2260            nwave = []
     2261            for i in xrange(in_nwave+1): nwave.append(i)
    22232262
    22242263        if not (isinstance(edge, list) or isinstance(edge, tuple)): edge = [ edge ]
    … …  
    22442283
    22452284        try:
    2246             #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method.
     2285            #CURRENTLY, PLOT=true is UNAVAILABLE UNTIL sinusoidal fitting is implemented as a fitter method.
    22472286            if individualedge:
    22482287                curedge = []
    … …  
    22502289                    curedge += edge[i]
    22512290               
    2252             workscan._auto_sinusoid_baseline(mask, minnwave, maxnwave, clipthresh, clipniter, curedge, threshold, chan_avg_limit, outlog, blfile)
     2291            workscan._auto_sinusoid_baseline(mask, nwave, maxwavelength, clipthresh, clipniter, curedge, threshold, chan_avg_limit, getresidual, outlog, blfile)
    22532292
    22542293            workscan._add_history("auto_sinusoid_baseline", varlist)
    … …  
    22952334            # also with 3-sigma clipping, iteration up to 4 times
    22962335            bscan = scan.cspline_baseline(npiece=2,clipthresh=3.0,clipniter=4)
     2336       
     2337        Note:
     2338            The best-fit parameter values output in logger and/or blfile are now
     2339            based on specunit of 'channel'.
    22972340        """
    22982341       
    … …  
    23882431        Example:
    23892432            bscan = scan.auto_cspline_baseline(npiece=3, insitu=False)
     2433       
     2434        Note:
     2435            The best-fit parameter values output in logger and/or blfile are now
     2436            based on specunit of 'channel'.
    23902437        """
    23912438

  • trunk/src/STLineFinder.cpp

    r2012 r2081  
    12031203    vel = scan->getAbcissa(last_row_used);
    12041204  } else {
    1205     for (int i = 0; i < spectrum.nelements(); ++i)
     1205    for (uInt i = 0; i < spectrum.nelements(); ++i)
    12061206      vel.push_back((double)i);
    12071207  }

  • trunk/src/Scantable.cpp

    r2080 r2081  
    17781778{
    17791779  ofstream ofs;
    1780   String coordInfo;
     1780  String coordInfo = "";
    17811781  bool hasSameNchan = true;
    17821782  bool outTextFile = false;
    … …  
    17951795  Fitter fitter = Fitter();
    17961796  fitter.setExpression("poly", order);
     1797  //fitter.setIterClipping(thresClip, nIterClip);
    17971798
    17981799  int nRow = nrow();
    … …  
    18121813{
    18131814  ofstream ofs;
    1814   String coordInfo;
     1815  String coordInfo = "";
    18151816  bool hasSameNchan = true;
    18161817  bool outTextFile = false;
    … …  
    18291830  Fitter fitter = Fitter();
    18301831  fitter.setExpression("poly", order);
     1832  //fitter.setIterClipping(thresClip, nIterClip);
    18311833
    18321834  int nRow = nrow();
    … …  
    18701872}
    18711873
    1872 void Scantable::cubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, bool outLogger, const std::string& blfile) {
     1874void Scantable::cubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, bool outLogger, const std::string& blfile)
     1875{
    18731876  ofstream ofs;
    1874   String coordInfo;
     1877  String coordInfo = "";
    18751878  bool hasSameNchan = true;
    18761879  bool outTextFile = false;
    … …  
    18891892  //Fitter fitter = Fitter();
    18901893  //fitter.setExpression("cspline", nPiece);
     1894  //fitter.setIterClipping(thresClip, nIterClip);
    18911895
    18921896  int nRow = nrow();
    … …  
    18951899  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    18961900    chanMask = getCompositeChanMask(whichrow, mask);
    1897     //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
     1901    //fitBaseline(chanMask, whichrow, fitter);
    18981902    //setSpectrum(fitter.getResidual(), whichrow);
    1899     std::vector<int> pieceRanges;
     1903    std::vector<int> pieceEdges;
    19001904    std::vector<float> params;
    1901     std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, pieceRanges, params, nPiece, thresClip, nIterClip, true);
     1905    std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, nPiece, pieceEdges, params, thresClip, nIterClip, true);
    19021906    setSpectrum(res, whichrow);
    19031907    //
    19041908
    1905     std::vector<bool>  fixed;
    1906     fixed.clear();
    1907     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()", pieceRanges, params, fixed);
     1909    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()", pieceEdges, params);
    19081910  }
    19091911
    … …  
    19141916{
    19151917  ofstream ofs;
    1916   String coordInfo;
     1918  String coordInfo = "";
    19171919  bool hasSameNchan = true;
    19181920  bool outTextFile = false;
    … …  
    19311933  //Fitter fitter = Fitter();
    19321934  //fitter.setExpression("cspline", nPiece);
     1935  //fitter.setIterClipping(thresClip, nIterClip);
    19331936
    19341937  int nRow = nrow();
    … …  
    19641967
    19651968
    1966     //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
     1969    //fitBaseline(chanMask, whichrow, fitter);
    19671970    //setSpectrum(fitter.getResidual(), whichrow);
    1968     std::vector<int> pieceRanges;
     1971    std::vector<int> pieceEdges;
    19691972    std::vector<float> params;
    1970     std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, pieceRanges, params, nPiece, thresClip, nIterClip, true);
     1973    std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, nPiece, pieceEdges, params, thresClip, nIterClip, true);
    19711974    setSpectrum(res, whichrow);
    19721975    //
    19731976
    1974     std::vector<bool> fixed;
    1975     fixed.clear();
    1976     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()", pieceRanges, params, fixed);
     1977    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()", pieceEdges, params);
    19771978  }
    19781979
    … …  
    19801981}
    19811982
    1982 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) {
     1983std::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)
     1984{
     1985  if (data.size() != mask.size()) {
     1986    throw(AipsError("data and mask sizes are not identical"));
     1987  }
    19831988  if (nPiece < 1) {
    19841989    throw(AipsError("wrong number of the sections for fitting"));
    1985   }
    1986   if (data.size() != mask.size()) {
    1987     throw(AipsError("data and mask have different size"));
    19881990  }
    19891991
    … …  
    19982000  }
    19992001
    2000   int nData = x.size();
    2001   int nElement = (int)(floor(floor((double)(nData/nPiece))+0.5));
     2002  int initNData = x.size();
     2003
     2004  int nElement = (int)(floor(floor((double)(initNData/nPiece))+0.5));
    20022005  std::vector<double> invEdge;
    2003 
    20042006  idxEdge.clear();
    20052007  idxEdge.push_back(x[0]);
    2006 
    20072008  for (int i = 1; i < nPiece; ++i) {
    20082009    int valX = x[nElement*i];
    … …  
    20102011    invEdge.push_back(1.0/(double)valX);
    20112012  }
    2012 
    20132013  idxEdge.push_back(x[x.size()-1]+1);
    20142014
    2015   std::vector<double> x1, x2, x3, z1, x1z1, x2z1, x3z1, r1;
     2015  int nData = initNData;
     2016  int nDOF = nPiece + 3;  //number of parameters to solve, namely, 4+(nPiece-1).
     2017
     2018  std::vector<double> x1, x2, x3, z1, x1z1, x2z1, x3z1, r1, residual;
    20162019  for (int i = 0; i < nChan; ++i) {
    20172020    double di = (double)i;
    … …  
    20252028    x3z1.push_back(dD*di*di*di);
    20262029    r1.push_back(0.0);
    2027   }
    2028 
    2029   int currentNData = nData;
    2030   int nDOF = nPiece + 3;  //number of parameters to solve, namely, 4+(nPiece-1).
     2030    residual.push_back(0.0);
     2031  }
    20312032
    20322033  for (int nClip = 0; nClip < nIterClip+1; ++nClip) {
    … …  
    21422143      for (int i = idxEdge[n]; i < idxEdge[n+1]; ++i) {
    21432144        r1[i] = a0 + a1*x1[i] + a2*x2[i] + a3*x3[i];
     2145        residual[i] = z1[i] - r1[i];
    21442146      }
    21452147      params.push_back(a0);
    … …  
    21612163      break;
    21622164    } else {
    2163       std::vector<double> rz;
    21642165      double stdDev = 0.0;
    21652166      for (int i = 0; i < nChan; ++i) {
    2166         double val = abs(r1[i] - z1[i]);
    2167         rz.push_back(val);
    2168         stdDev += val*val*(double)maskArray[i];
     2167        stdDev += residual[i]*residual[i]*(double)maskArray[i];
    21692168      }
    21702169      stdDev = sqrt(stdDev/(double)nData);
    … …  
    21732172      int newNData = 0;
    21742173      for (int i = 0; i < nChan; ++i) {
    2175         if (rz[i] >= thres) {
     2174        if (abs(residual[i]) >= thres) {
    21762175          maskArray[i] = 0;
    21772176        }
    … …  
    21802179        }
    21812180      }
    2182       if (newNData == currentNData) {
     2181      if (newNData == nData) {
    21832182        break; //no more flag to add. iteration stops.
    21842183      } else {
    2185         currentNData = newNData;
     2184        nData = newNData;
    21862185      }
    21872186    }
    … …  
    21912190  if (getResidual) {
    21922191    for (int i = 0; i < nChan; ++i) {
    2193       result.push_back((float)(z1[i] - r1[i]));
     2192      result.push_back((float)residual[i]);
    21942193    }
    21952194  } else {
    … …  
    22022201}
    22032202
    2204 void Scantable::sinusoidBaseline(const std::vector<bool>& mask, int nMinWavesInSW, int nMaxWavesInSW, float thresClip, int nIterClip, bool outLogger, const std::string& blfile) {
     2203void 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)
     2204{
    22052205  ofstream ofs;
    2206   String coordInfo;
     2206  String coordInfo = "";
    22072207  bool hasSameNchan = true;
    22082208  bool outTextFile = false;
    … …  
    22202220
    22212221  //Fitter fitter = Fitter();
    2222   //fitter.setExpression("sinusoid", nMaxWavesInSW);
     2222  //fitter.setExpression("sinusoid", nWaves);
     2223  //fitter.setIterClipping(thresClip, nIterClip);
    22232224
    22242225  int nRow = nrow();
    … …  
    22272228  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    22282229    chanMask = getCompositeChanMask(whichrow, mask);
    2229     //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
     2230    //fitBaseline(chanMask, whichrow, fitter);
    22302231    //setSpectrum(fitter.getResidual(), whichrow);
    2231     std::vector<int> pieceRanges;
    22322232    std::vector<float> params;
    2233     std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nMinWavesInSW, nMaxWavesInSW, params, thresClip, nIterClip, true);
     2233    //std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, params, thresClip, nIterClip, true);
     2234    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, maxWaveLength, params, thresClip, nIterClip, getResidual);
    22342235    setSpectrum(res, whichrow);
    22352236    //
    22362237
    2237     std::vector<bool>  fixed;
    2238     fixed.clear();
    2239     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "sinusoidBaseline()", pieceRanges, params, fixed);
     2238    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "sinusoidBaseline()", params);
    22402239  }
    22412240
    … …  
    22432242}
    22442243
    2245 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)
     2244void 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)
    22462245{
    22472246  ofstream ofs;
    2248   String coordInfo;
     2247  String coordInfo = "";
    22492248  bool hasSameNchan = true;
    22502249  bool outTextFile = false;
    … …  
    22622261
    22632262  //Fitter fitter = Fitter();
    2264   //fitter.setExpression("sinusoid", nMaxWavesInSW);
     2263  //fitter.setExpression("sinusoid", nWaves);
     2264  //fitter.setIterClipping(thresClip, nIterClip);
    22652265
    22662266  int nRow = nrow();
    … …  
    22962296
    22972297
    2298     //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
     2298    //fitBaseline(chanMask, whichrow, fitter);
    22992299    //setSpectrum(fitter.getResidual(), whichrow);
    2300     std::vector<int> pieceRanges;
    23012300    std::vector<float> params;
    2302     std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nMinWavesInSW, nMaxWavesInSW, params, thresClip, nIterClip, true);
     2301    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, maxWaveLength, params, thresClip, nIterClip, getResidual);
    23032302    setSpectrum(res, whichrow);
    23042303    //
    23052304
    2306     std::vector<bool> fixed;
    2307     fixed.clear();
    2308     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoSinusoidBaseline()", pieceRanges, params, fixed);
     2305    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoSinusoidBaseline()", params);
    23092306  }
    23102307
    … …  
    23122309}
    23132310
    2314 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) {
     2311std::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)
     2312{
    23152313  if (data.size() != mask.size()) {
    2316     throw(AipsError("data and mask have different size"));
    2317   }
    2318   if ((nMinWavesInSW < 0) || (nMaxWavesInSW < 0)) {
     2314    throw(AipsError("data and mask sizes are not identical"));
     2315  }
     2316  if (data.size() < 2) {
     2317    throw(AipsError("data size is too short"));
     2318  }
     2319  if (waveNumbers.size() == 0) {
     2320    throw(AipsError("missing wave number info"));
     2321  }
     2322  std::vector<int> nWaves;  // sorted and uniqued array of wave numbers
     2323  nWaves.reserve(waveNumbers.size());
     2324  copy(waveNumbers.begin(), waveNumbers.end(), back_inserter(nWaves));
     2325  sort(nWaves.begin(), nWaves.end());
     2326  std::vector<int>::iterator end_it = unique(nWaves.begin(), nWaves.end());
     2327  nWaves.erase(end_it, nWaves.end());
     2328
     2329  int minNWaves = nWaves[0];
     2330  if (minNWaves < 0) {
    23192331    throw(AipsError("wave number must be positive or zero (i.e. constant)"));
    2320   } else {
    2321     if (nMaxWavesInSW < nMinWavesInSW) {
    2322       int temp = nMaxWavesInSW;
    2323       nMaxWavesInSW = nMinWavesInSW;
    2324       nMinWavesInSW = temp;
    2325     }
    2326   }
     2332  }
     2333  bool hasConstantTerm = (minNWaves == 0);
    23272334
    23282335  int nChan = data.size();
    … …  
    23362343  }
    23372344
    2338   int nData = x.size();
    2339 
    2340   std::vector<double> x1, x2, x3, z1, x1z1, x2z1, x3z1, r1;
     2345  int initNData = x.size();
     2346
     2347  int nData = initNData;
     2348  int nDOF = nWaves.size() * 2 - (hasConstantTerm ? 1 : 0);  //number of parameters to solve.
     2349
     2350  const double PI = 6.0 * asin(0.5); // PI (= 3.141592653...)
     2351  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.
     2352
     2353  // xArray : contains elemental values for computing the least-square matrix.
     2354  //          xArray.size() is nDOF and xArray[*].size() is nChan.
     2355  //          Each xArray element are as follows:
     2356  //          xArray[0]    = {1.0, 1.0, 1.0, ..., 1.0},
     2357  //          xArray[2n-1] = {sin(nPI/L*x[0]), sin(nPI/L*x[1]), ..., sin(nPI/L*x[nChan])},
     2358  //          xArray[2n]   = {cos(nPI/L*x[0]), cos(nPI/L*x[1]), ..., cos(nPI/L*x[nChan])},
     2359  //          where (1 <= n <= nMaxWavesInSW),
     2360  //          or,
     2361  //          xArray[2n-1] = {sin(wn[n]PI/L*x[0]), sin(wn[n]PI/L*x[1]), ..., sin(wn[n]PI/L*x[nChan])},
     2362  //          xArray[2n]   = {cos(wn[n]PI/L*x[0]), cos(wn[n]PI/L*x[1]), ..., cos(wn[n]PI/L*x[nChan])},
     2363  //          where wn[n] denotes waveNumbers[n] (1 <= n <= waveNumbers.size()).
     2364  std::vector<std::vector<double> > xArray;
     2365  if (hasConstantTerm) {
     2366    std::vector<double> xu;
     2367    for (int j = 0; j < nChan; ++j) {
     2368      xu.push_back(1.0);
     2369    }
     2370    xArray.push_back(xu);
     2371  }
     2372  for (uInt i = (hasConstantTerm ? 1 : 0); i < nWaves.size(); ++i) {
     2373    double xFactor = baseXFactor*(double)nWaves[i];
     2374    std::vector<double> xs, xc;
     2375    xs.clear();
     2376    xc.clear();
     2377    for (int j = 0; j < nChan; ++j) {
     2378      xs.push_back(sin(xFactor*(double)j));
     2379      xc.push_back(cos(xFactor*(double)j));
     2380    }
     2381    xArray.push_back(xs);
     2382    xArray.push_back(xc);
     2383  }
     2384
     2385  std::vector<double> z1, r1, residual;
    23412386  for (int i = 0; i < nChan; ++i) {
    2342     double di = (double)i;
    2343     double dD = (double)data[i];
    2344     x1.push_back(di);
    2345     x2.push_back(di*di);
    2346     x3.push_back(di*di*di);
    2347     z1.push_back(dD);
    2348     x1z1.push_back(di*dD);
    2349     x2z1.push_back(di*di*dD);
    2350     x3z1.push_back(di*di*di*dD);
     2387    z1.push_back((double)data[i]);
    23512388    r1.push_back(0.0);
    2352   }
    2353 
    2354   int currentNData = nData;
    2355   int nDOF = nMaxWavesInSW + 1;  //number of parameters to solve.
     2389    residual.push_back(0.0);
     2390  }
    23562391
    23572392  for (int nClip = 0; nClip < nIterClip+1; ++nClip) {
    2358     //xMatrix : horizontal concatenation of the least-sq. matrix (left) and an
    2359     //identity matrix (right).
    2360     //the right part is used to calculate the inverse matrix of the left part.
     2393    // xMatrix : horizontal concatenation of
     2394    //           the least-sq. matrix (left) and an
     2395    //           identity matrix (right).
     2396    // the right part is used to calculate the inverse matrix of the left part.
    23612397    double xMatrix[nDOF][2*nDOF];
    23622398    double zMatrix[nDOF];
    … …  
    23692405    }
    23702406
    2371     for (int i = 0; i < currentNData; ++i) {
    2372       if (maskArray[i] == 0) continue;
    2373       xMatrix[0][0] += 1.0;
    2374       xMatrix[0][1] += x1[i];
    2375       xMatrix[0][2] += x2[i];
    2376       xMatrix[0][3] += x3[i];
    2377       xMatrix[1][1] += x2[i];
    2378       xMatrix[1][2] += x3[i];
    2379       xMatrix[1][3] += x2[i]*x2[i];
    2380       xMatrix[2][2] += x2[i]*x2[i];
    2381       xMatrix[2][3] += x3[i]*x2[i];
    2382       xMatrix[3][3] += x3[i]*x3[i];
    2383       zMatrix[0] += z1[i];
    2384       zMatrix[1] += x1z1[i];
    2385       zMatrix[2] += x2z1[i];
    2386       zMatrix[3] += x3z1[i];
     2407    for (int k = 0; k < nChan; ++k) {
     2408      if (maskArray[k] == 0) continue;
     2409
     2410      for (int i = 0; i < nDOF; ++i) {
     2411        for (int j = i; j < nDOF; ++j) {
     2412          xMatrix[i][j] += xArray[i][k] * xArray[j][k];
     2413        }
     2414        zMatrix[i] += z1[k] * xArray[i][k];
     2415      }
    23872416    }
    23882417
    … …  
    24252454    }
    24262455    //compute a vector y which consists of the coefficients of the sinusoids forming the
    2427     //best-fit curves (a0,b0,a1,b1,...), namely, a* and b* are of sine and cosine functions,
    2428     //respectively.
     2456    //best-fit curves (a0,s1,c1,s2,c2,...), where a0 is constant and s* and c* are of sine
     2457    //and cosine functions, respectively.
    24292458    std::vector<double> y;
     2459    params.clear();
    24302460    for (int i = 0; i < nDOF; ++i) {
    24312461      y.push_back(0.0);
    … …  
    24332463        y[i] += xMatrix[i][nDOF+j]*zMatrix[j];
    24342464      }
    2435     }
    2436 
    2437     double a0 = y[0];
    2438     double a1 = y[1];
    2439     double a2 = y[2];
    2440     double a3 = y[3];
    2441     params.clear();
     2465      params.push_back(y[i]);
     2466    }
    24422467
    24432468    for (int i = 0; i < nChan; ++i) {
    2444       r1[i] = a0 + a1*x1[i] + a2*x2[i] + a3*x3[i];
    2445     }
    2446     params.push_back(a0);
    2447     params.push_back(a1);
    2448     params.push_back(a2);
    2449     params.push_back(a3);
    2450 
     2469      r1[i] = y[0];
     2470      for (int j = 1; j < nDOF; ++j) {
     2471        r1[i] += y[j]*xArray[j][i];
     2472      }
     2473      residual[i] = z1[i] - r1[i];
     2474    }
    24512475
    24522476    if ((nClip == nIterClip) || (thresClip <= 0.0)) {
    24532477      break;
    24542478    } else {
    2455       std::vector<double> rz;
    24562479      double stdDev = 0.0;
    24572480      for (int i = 0; i < nChan; ++i) {
    2458         double val = abs(r1[i] - z1[i]);
    2459         rz.push_back(val);
    2460         stdDev += val*val*(double)maskArray[i];
     2481        stdDev += residual[i]*residual[i]*(double)maskArray[i];
    24612482      }
    24622483      stdDev = sqrt(stdDev/(double)nData);
    … …  
    24652486      int newNData = 0;
    24662487      for (int i = 0; i < nChan; ++i) {
    2467         if (rz[i] >= thres) {
     2488        if (abs(residual[i]) >= thres) {
    24682489          maskArray[i] = 0;
    24692490        }
    … …  
    24722493        }
    24732494      }
    2474       if (newNData == currentNData) {
    2475         break;                   //no additional flag. finish iteration.
     2495      if (newNData == nData) {
     2496        break; //no more flag to add. iteration stops.
    24762497      } else {
    2477         currentNData = newNData;
     2498        nData = newNData;
    24782499      }
    24792500    }
    … …  
    24832504  if (getResidual) {
    24842505    for (int i = 0; i < nChan; ++i) {
    2485       result.push_back((float)(z1[i] - r1[i]));
     2506      result.push_back((float)residual[i]);
    24862507    }
    24872508  } else {
    … …  
    24962517void Scantable::fitBaseline(const std::vector<bool>& mask, int whichrow, Fitter& fitter)
    24972518{
    2498   std::vector<double> abcsd = getAbcissa(whichrow);
    2499   std::vector<float> abcs;
    2500   for (uInt i = 0; i < abcsd.size(); ++i) {
    2501     abcs.push_back((float)abcsd[i]);
     2519  std::vector<double> dAbcissa = getAbcissa(whichrow);
     2520  std::vector<float> abcissa;
     2521  for (uInt i = 0; i < dAbcissa.size(); ++i) {
     2522    abcissa.push_back((float)dAbcissa[i]);
    25022523  }
    25032524  std::vector<float> spec = getSpectrum(whichrow);
    25042525
    2505   fitter.setData(abcs, spec, mask);
     2526  fitter.setData(abcissa, spec, mask);
    25062527  fitter.lfit();
    25072528}
    … …  
    25662587
    25672588/* for cspline. will be merged once cspline is available in fitter (2011/3/10 WK) */
    2568 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) {
     2589void 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) {
    25692590  if (outLogger || outTextFile) {
    25702591    float rms = getRms(chanMask, whichrow);
    25712592    String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
     2593    std::vector<bool> fixed;
     2594    fixed.clear();
    25722595
    25732596    if (outLogger) {
    … …  
    25822605
    25832606/* for sinusoid. will be merged once sinusoid is available in fitter (2011/3/10 WK) */
    2584 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) {
     2607void 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) {
    25852608  if (outLogger || outTextFile) {
    25862609    float rms = getRms(chanMask, whichrow);
    25872610    String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
     2611    std::vector<bool> fixed;
     2612    fixed.clear();
    25882613
    25892614    if (outLogger) {

  • trunk/src/Scantable.h

    r2064 r2081  
    522522                               const std::string& blfile="");
    523523  void sinusoidBaseline(const std::vector<bool>& mask,
    524                         int nMinWavesInSW,
    525                         int nMaxWavesInSW,
     524                        const std::vector<int>& nWaves,
     525                        float maxWaveLength,
    526526                        float thresClip,
    527527                        int nIterClip,
     528                        bool getResidual=true,
    528529                        bool outLogger=false,
    529530                        const std::string& blfile="");
    530531  void autoSinusoidBaseline(const std::vector<bool>& mask,
    531                             int nMinWavesInSW,
    532                             int nMaxWavesInSW,
     532                            const std::vector<int>& nWaves,
     533                            float maxWaveLength,
    533534                            float thresClip,
    534535                            int nIterClip,
    … …  
    536537                            float threshold=3.0,
    537538                            int chanAvgLimit=1,
     539                            bool getResidual=true,
    538540                            bool outLogger=false,
    539541                            const std::string& blfile="");
    … …  
    672674  std::vector<float> doCubicSplineFitting(const std::vector<float>& data,
    673675                                          const std::vector<bool>& mask,
     676                                          int nPiece,
    674677                                          std::vector<int>& idxEdge,
    675678                                          std::vector<float>& params,
    676                                           int nPiece=2,
    677679                                          float thresClip=3.0,
    678680                                          int nIterClip=1,
    … …  
    680682  std::vector<float> doSinusoidFitting(const std::vector<float>& data,
    681683                                       const std::vector<bool>& mask,
    682                                        int nMinWavesInSW,
    683                                        int nMaxWavesInSW,
     684                                       const std::vector<int>& waveNumbers,
     685                                       float maxWaveLength,
    684686                                       std::vector<float>& params,
    685687                                       float thresClip=3.0,
    … …  
    698700  //std::vector<bool> getCompositeChanMask(int whichrow, const std::vector<bool>& inMask, const std::vector<int>& edge, const int minEdgeSize, STLineFinder& lineFinder);
    699701  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);
    700   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);
    701   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);
     702  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);
     703  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);
    702704
    703705  void applyChanFlag( casa::uInt whichrow, const std::vector<bool>& msk, casa::uChar flagval);

  • trunk/src/ScantableWrapper.h

    r2047 r2081  
    269269  { table_->autoCubicSplineBaseline(mask, npiece, clipthresh, clipniter, edge, threshold, chan_avg_limit, outlog, blfile); }
    270270
    271   void sinusoidBaseline(const std::vector<bool>& mask, int minnwave, int maxnwave, float clipthresh, int clipniter, bool outlog=false, const std::string& blfile="")
    272   { table_->sinusoidBaseline(mask, minnwave, maxnwave, clipthresh, clipniter, outlog, blfile); }
    273 
    274   void autoSinusoidBaseline(const std::vector<bool>& mask, int minnwave, int maxnwave, float clipthresh, int clipniter, const std::vector<int>& edge, float threshold=5.0, int chan_avg_limit=1, bool outlog=false, const std::string& blfile="")
    275   { table_->autoSinusoidBaseline(mask, minnwave, maxnwave, clipthresh, clipniter, edge, threshold, chan_avg_limit, outlog, blfile); }
     271  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="")
     272  { table_->sinusoidBaseline(mask, nwave, maxwavelength, clipthresh, clipniter, getresidual, outlog, blfile); }
     273
     274  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="")
     275  { table_->autoSinusoidBaseline(mask, nwave, maxwavelength, clipthresh, clipniter, edge, threshold, chan_avg_limit, getresidual, outlog, blfile); }
    276276
    277277  float getRms(const std::vector<bool>& mask, int whichrow)
Note: See TracChangeset for help on using the changeset viewer.