Changeset 2082 for branches


Ignore:
Timestamp:
03/30/11 21:43:43 (14 years ago)
Author:
WataruKawasaki
Message:

merged bug fixes from trunk (r2081)

Location:
branches/casa-prerelease/pre-asap
Files:
15 edited

Legend:

Unmodified
Added
Removed

  • branches/casa-prerelease/pre-asap

    • Property svn:mergeinfo changed
      /trunkmerged: 2081

  • branches/casa-prerelease/pre-asap/Makefile

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • branches/casa-prerelease/pre-asap/SConstruct

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • branches/casa-prerelease/pre-asap/apps

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • branches/casa-prerelease/pre-asap/external-alma

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • branches/casa-prerelease/pre-asap/external-alma/atnf/pks/pks_maths.cc

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • branches/casa-prerelease/pre-asap/getsvnrev.sh

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • branches/casa-prerelease/pre-asap/python

  • branches/casa-prerelease/pre-asap/python/scantable.py

    r2047 r2082  
    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

  • branches/casa-prerelease/pre-asap/src

  • branches/casa-prerelease/pre-asap/src/SConscript

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • branches/casa-prerelease/pre-asap/src/STLineFinder.cpp

    r2012 r2082  
    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  }

  • branches/casa-prerelease/pre-asap/src/Scantable.cpp

    r2065 r2082  
    17641764{
    17651765  ofstream ofs;
    1766   String coordInfo;
     1766  String coordInfo = "";
    17671767  bool hasSameNchan = true;
    17681768  bool outTextFile = false;
    … …  
    17811781  Fitter fitter = Fitter();
    17821782  fitter.setExpression("poly", order);
     1783  //fitter.setIterClipping(thresClip, nIterClip);
    17831784
    17841785  int nRow = nrow();
    … …  
    17981799{
    17991800  ofstream ofs;
    1800   String coordInfo;
     1801  String coordInfo = "";
    18011802  bool hasSameNchan = true;
    18021803  bool outTextFile = false;
    … …  
    18151816  Fitter fitter = Fitter();
    18161817  fitter.setExpression("poly", order);
     1818  //fitter.setIterClipping(thresClip, nIterClip);
    18171819
    18181820  int nRow = nrow();
    … …  
    18561858}
    18571859
    1858 void Scantable::cubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, bool outLogger, const std::string& blfile) {
     1860void Scantable::cubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, bool outLogger, const std::string& blfile)
     1861{
    18591862  ofstream ofs;
    1860   String coordInfo;
     1863  String coordInfo = "";
    18611864  bool hasSameNchan = true;
    18621865  bool outTextFile = false;
    … …  
    18751878  //Fitter fitter = Fitter();
    18761879  //fitter.setExpression("cspline", nPiece);
     1880  //fitter.setIterClipping(thresClip, nIterClip);
    18771881
    18781882  int nRow = nrow();
    … …  
    18811885  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    18821886    chanMask = getCompositeChanMask(whichrow, mask);
    1883     //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
     1887    //fitBaseline(chanMask, whichrow, fitter);
    18841888    //setSpectrum(fitter.getResidual(), whichrow);
    1885     std::vector<int> pieceRanges;
     1889    std::vector<int> pieceEdges;
    18861890    std::vector<float> params;
    1887     std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, pieceRanges, params, nPiece, thresClip, nIterClip, true);
     1891    std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, nPiece, pieceEdges, params, thresClip, nIterClip, true);
    18881892    setSpectrum(res, whichrow);
    18891893    //
    18901894
    1891     std::vector<bool>  fixed;
    1892     fixed.clear();
    1893     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()", pieceRanges, params, fixed);
     1895    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()", pieceEdges, params);
    18941896  }
    18951897
    … …  
    19001902{
    19011903  ofstream ofs;
    1902   String coordInfo;
     1904  String coordInfo = "";
    19031905  bool hasSameNchan = true;
    19041906  bool outTextFile = false;
    … …  
    19171919  //Fitter fitter = Fitter();
    19181920  //fitter.setExpression("cspline", nPiece);
     1921  //fitter.setIterClipping(thresClip, nIterClip);
    19191922
    19201923  int nRow = nrow();
    … …  
    19501953
    19511954
    1952     //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
     1955    //fitBaseline(chanMask, whichrow, fitter);
    19531956    //setSpectrum(fitter.getResidual(), whichrow);
    1954     std::vector<int> pieceRanges;
     1957    std::vector<int> pieceEdges;
    19551958    std::vector<float> params;
    1956     std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, pieceRanges, params, nPiece, thresClip, nIterClip, true);
     1959    std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, nPiece, pieceEdges, params, thresClip, nIterClip, true);
    19571960    setSpectrum(res, whichrow);
    19581961    //
    19591962
    1960     std::vector<bool> fixed;
    1961     fixed.clear();
    1962     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()", pieceRanges, params, fixed);
     1963    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()", pieceEdges, params);
    19631964  }
    19641965
    … …  
    19661967}
    19671968
    1968 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) {
     1969std::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)
     1970{
     1971  if (data.size() != mask.size()) {
     1972    throw(AipsError("data and mask sizes are not identical"));
     1973  }
    19691974  if (nPiece < 1) {
    19701975    throw(AipsError("wrong number of the sections for fitting"));
    1971   }
    1972   if (data.size() != mask.size()) {
    1973     throw(AipsError("data and mask have different size"));
    19741976  }
    19751977
    … …  
    19841986  }
    19851987
    1986   int nData = x.size();
    1987   int nElement = (int)(floor(floor((double)(nData/nPiece))+0.5));
     1988  int initNData = x.size();
     1989
     1990  int nElement = (int)(floor(floor((double)(initNData/nPiece))+0.5));
    19881991  std::vector<double> invEdge;
    1989 
    19901992  idxEdge.clear();
    19911993  idxEdge.push_back(x[0]);
    1992 
    19931994  for (int i = 1; i < nPiece; ++i) {
    19941995    int valX = x[nElement*i];
    … …  
    19961997    invEdge.push_back(1.0/(double)valX);
    19971998  }
    1998 
    19991999  idxEdge.push_back(x[x.size()-1]+1);
    20002000
    2001   std::vector<double> x1, x2, x3, z1, x1z1, x2z1, x3z1, r1;
     2001  int nData = initNData;
     2002  int nDOF = nPiece + 3;  //number of parameters to solve, namely, 4+(nPiece-1).
     2003
     2004  std::vector<double> x1, x2, x3, z1, x1z1, x2z1, x3z1, r1, residual;
    20022005  for (int i = 0; i < nChan; ++i) {
    20032006    double di = (double)i;
    … …  
    20112014    x3z1.push_back(dD*di*di*di);
    20122015    r1.push_back(0.0);
    2013   }
    2014 
    2015   int currentNData = nData;
    2016   int nDOF = nPiece + 3;  //number of parameters to solve, namely, 4+(nPiece-1).
     2016    residual.push_back(0.0);
     2017  }
    20172018
    20182019  for (int nClip = 0; nClip < nIterClip+1; ++nClip) {
    … …  
    21282129      for (int i = idxEdge[n]; i < idxEdge[n+1]; ++i) {
    21292130        r1[i] = a0 + a1*x1[i] + a2*x2[i] + a3*x3[i];
     2131        residual[i] = z1[i] - r1[i];
    21302132      }
    21312133      params.push_back(a0);
    … …  
    21472149      break;
    21482150    } else {
    2149       std::vector<double> rz;
    21502151      double stdDev = 0.0;
    21512152      for (int i = 0; i < nChan; ++i) {
    2152         double val = abs(r1[i] - z1[i]);
    2153         rz.push_back(val);
    2154         stdDev += val*val*(double)maskArray[i];
     2153        stdDev += residual[i]*residual[i]*(double)maskArray[i];
    21552154      }
    21562155      stdDev = sqrt(stdDev/(double)nData);
    … …  
    21592158      int newNData = 0;
    21602159      for (int i = 0; i < nChan; ++i) {
    2161         if (rz[i] >= thres) {
     2160        if (abs(residual[i]) >= thres) {
    21622161          maskArray[i] = 0;
    21632162        }
    … …  
    21662165        }
    21672166      }
    2168       if (newNData == currentNData) {
     2167      if (newNData == nData) {
    21692168        break; //no more flag to add. iteration stops.
    21702169      } else {
    2171         currentNData = newNData;
     2170        nData = newNData;
    21722171      }
    21732172    }
    … …  
    21772176  if (getResidual) {
    21782177    for (int i = 0; i < nChan; ++i) {
    2179       result.push_back((float)(z1[i] - r1[i]));
     2178      result.push_back((float)residual[i]);
    21802179    }
    21812180  } else {
    … …  
    21882187}
    21892188
    2190 void Scantable::sinusoidBaseline(const std::vector<bool>& mask, int nMinWavesInSW, int nMaxWavesInSW, float thresClip, int nIterClip, bool outLogger, const std::string& blfile) {
     2189void 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)
     2190{
    21912191  ofstream ofs;
    2192   String coordInfo;
     2192  String coordInfo = "";
    21932193  bool hasSameNchan = true;
    21942194  bool outTextFile = false;
    … …  
    22062206
    22072207  //Fitter fitter = Fitter();
    2208   //fitter.setExpression("sinusoid", nMaxWavesInSW);
     2208  //fitter.setExpression("sinusoid", nWaves);
     2209  //fitter.setIterClipping(thresClip, nIterClip);
    22092210
    22102211  int nRow = nrow();
    … …  
    22132214  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    22142215    chanMask = getCompositeChanMask(whichrow, mask);
    2215     //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
     2216    //fitBaseline(chanMask, whichrow, fitter);
    22162217    //setSpectrum(fitter.getResidual(), whichrow);
    2217     std::vector<int> pieceRanges;
    22182218    std::vector<float> params;
    2219     std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nMinWavesInSW, nMaxWavesInSW, params, thresClip, nIterClip, true);
     2219    //std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, params, thresClip, nIterClip, true);
     2220    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, maxWaveLength, params, thresClip, nIterClip, getResidual);
    22202221    setSpectrum(res, whichrow);
    22212222    //
    22222223
    2223     std::vector<bool>  fixed;
    2224     fixed.clear();
    2225     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "sinusoidBaseline()", pieceRanges, params, fixed);
     2224    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "sinusoidBaseline()", params);
    22262225  }
    22272226
    … …  
    22292228}
    22302229
    2231 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)
     2230void 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)
    22322231{
    22332232  ofstream ofs;
    2234   String coordInfo;
     2233  String coordInfo = "";
    22352234  bool hasSameNchan = true;
    22362235  bool outTextFile = false;
    … …  
    22482247
    22492248  //Fitter fitter = Fitter();
    2250   //fitter.setExpression("sinusoid", nMaxWavesInSW);
     2249  //fitter.setExpression("sinusoid", nWaves);
     2250  //fitter.setIterClipping(thresClip, nIterClip);
    22512251
    22522252  int nRow = nrow();
    … …  
    22822282
    22832283
    2284     //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
     2284    //fitBaseline(chanMask, whichrow, fitter);
    22852285    //setSpectrum(fitter.getResidual(), whichrow);
    2286     std::vector<int> pieceRanges;
    22872286    std::vector<float> params;
    2288     std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nMinWavesInSW, nMaxWavesInSW, params, thresClip, nIterClip, true);
     2287    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, maxWaveLength, params, thresClip, nIterClip, getResidual);
    22892288    setSpectrum(res, whichrow);
    22902289    //
    22912290
    2292     std::vector<bool> fixed;
    2293     fixed.clear();
    2294     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoSinusoidBaseline()", pieceRanges, params, fixed);
     2291    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoSinusoidBaseline()", params);
    22952292  }
    22962293
    … …  
    22982295}
    22992296
    2300 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) {
     2297std::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)
     2298{
    23012299  if (data.size() != mask.size()) {
    2302     throw(AipsError("data and mask have different size"));
    2303   }
    2304   if ((nMinWavesInSW < 0) || (nMaxWavesInSW < 0)) {
     2300    throw(AipsError("data and mask sizes are not identical"));
     2301  }
     2302  if (data.size() < 2) {
     2303    throw(AipsError("data size is too short"));
     2304  }
     2305  if (waveNumbers.size() == 0) {
     2306    throw(AipsError("missing wave number info"));
     2307  }
     2308  std::vector<int> nWaves;  // sorted and uniqued array of wave numbers
     2309  nWaves.reserve(waveNumbers.size());
     2310  copy(waveNumbers.begin(), waveNumbers.end(), back_inserter(nWaves));
     2311  sort(nWaves.begin(), nWaves.end());
     2312  std::vector<int>::iterator end_it = unique(nWaves.begin(), nWaves.end());
     2313  nWaves.erase(end_it, nWaves.end());
     2314
     2315  int minNWaves = nWaves[0];
     2316  if (minNWaves < 0) {
    23052317    throw(AipsError("wave number must be positive or zero (i.e. constant)"));
    2306   } else {
    2307     if (nMaxWavesInSW < nMinWavesInSW) {
    2308       int temp = nMaxWavesInSW;
    2309       nMaxWavesInSW = nMinWavesInSW;
    2310       nMinWavesInSW = temp;
    2311     }
    2312   }
     2318  }
     2319  bool hasConstantTerm = (minNWaves == 0);
    23132320
    23142321  int nChan = data.size();
    … …  
    23222329  }
    23232330
    2324   int nData = x.size();
    2325 
    2326   std::vector<double> x1, x2, x3, z1, x1z1, x2z1, x3z1, r1;
     2331  int initNData = x.size();
     2332
     2333  int nData = initNData;
     2334  int nDOF = nWaves.size() * 2 - (hasConstantTerm ? 1 : 0);  //number of parameters to solve.
     2335
     2336  const double PI = 6.0 * asin(0.5); // PI (= 3.141592653...)
     2337  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.
     2338
     2339  // xArray : contains elemental values for computing the least-square matrix.
     2340  //          xArray.size() is nDOF and xArray[*].size() is nChan.
     2341  //          Each xArray element are as follows:
     2342  //          xArray[0]    = {1.0, 1.0, 1.0, ..., 1.0},
     2343  //          xArray[2n-1] = {sin(nPI/L*x[0]), sin(nPI/L*x[1]), ..., sin(nPI/L*x[nChan])},
     2344  //          xArray[2n]   = {cos(nPI/L*x[0]), cos(nPI/L*x[1]), ..., cos(nPI/L*x[nChan])},
     2345  //          where (1 <= n <= nMaxWavesInSW),
     2346  //          or,
     2347  //          xArray[2n-1] = {sin(wn[n]PI/L*x[0]), sin(wn[n]PI/L*x[1]), ..., sin(wn[n]PI/L*x[nChan])},
     2348  //          xArray[2n]   = {cos(wn[n]PI/L*x[0]), cos(wn[n]PI/L*x[1]), ..., cos(wn[n]PI/L*x[nChan])},
     2349  //          where wn[n] denotes waveNumbers[n] (1 <= n <= waveNumbers.size()).
     2350  std::vector<std::vector<double> > xArray;
     2351  if (hasConstantTerm) {
     2352    std::vector<double> xu;
     2353    for (int j = 0; j < nChan; ++j) {
     2354      xu.push_back(1.0);
     2355    }
     2356    xArray.push_back(xu);
     2357  }
     2358  for (uInt i = (hasConstantTerm ? 1 : 0); i < nWaves.size(); ++i) {
     2359    double xFactor = baseXFactor*(double)nWaves[i];
     2360    std::vector<double> xs, xc;
     2361    xs.clear();
     2362    xc.clear();
     2363    for (int j = 0; j < nChan; ++j) {
     2364      xs.push_back(sin(xFactor*(double)j));
     2365      xc.push_back(cos(xFactor*(double)j));
     2366    }
     2367    xArray.push_back(xs);
     2368    xArray.push_back(xc);
     2369  }
     2370
     2371  std::vector<double> z1, r1, residual;
    23272372  for (int i = 0; i < nChan; ++i) {
    2328     double di = (double)i;
    2329     double dD = (double)data[i];
    2330     x1.push_back(di);
    2331     x2.push_back(di*di);
    2332     x3.push_back(di*di*di);
    2333     z1.push_back(dD);
    2334     x1z1.push_back(di*dD);
    2335     x2z1.push_back(di*di*dD);
    2336     x3z1.push_back(di*di*di*dD);
     2373    z1.push_back((double)data[i]);
    23372374    r1.push_back(0.0);
    2338   }
    2339 
    2340   int currentNData = nData;
    2341   int nDOF = nMaxWavesInSW + 1;  //number of parameters to solve.
     2375    residual.push_back(0.0);
     2376  }
    23422377
    23432378  for (int nClip = 0; nClip < nIterClip+1; ++nClip) {
    2344     //xMatrix : horizontal concatenation of the least-sq. matrix (left) and an
    2345     //identity matrix (right).
    2346     //the right part is used to calculate the inverse matrix of the left part.
     2379    // xMatrix : horizontal concatenation of
     2380    //           the least-sq. matrix (left) and an
     2381    //           identity matrix (right).
     2382    // the right part is used to calculate the inverse matrix of the left part.
    23472383    double xMatrix[nDOF][2*nDOF];
    23482384    double zMatrix[nDOF];
    … …  
    23552391    }
    23562392
    2357     for (int i = 0; i < currentNData; ++i) {
    2358       if (maskArray[i] == 0) continue;
    2359       xMatrix[0][0] += 1.0;
    2360       xMatrix[0][1] += x1[i];
    2361       xMatrix[0][2] += x2[i];
    2362       xMatrix[0][3] += x3[i];
    2363       xMatrix[1][1] += x2[i];
    2364       xMatrix[1][2] += x3[i];
    2365       xMatrix[1][3] += x2[i]*x2[i];
    2366       xMatrix[2][2] += x2[i]*x2[i];
    2367       xMatrix[2][3] += x3[i]*x2[i];
    2368       xMatrix[3][3] += x3[i]*x3[i];
    2369       zMatrix[0] += z1[i];
    2370       zMatrix[1] += x1z1[i];
    2371       zMatrix[2] += x2z1[i];
    2372       zMatrix[3] += x3z1[i];
     2393    for (int k = 0; k < nChan; ++k) {
     2394      if (maskArray[k] == 0) continue;
     2395
     2396      for (int i = 0; i < nDOF; ++i) {
     2397        for (int j = i; j < nDOF; ++j) {
     2398          xMatrix[i][j] += xArray[i][k] * xArray[j][k];
     2399        }
     2400        zMatrix[i] += z1[k] * xArray[i][k];
     2401      }
    23732402    }
    23742403
    … …  
    24112440    }
    24122441    //compute a vector y which consists of the coefficients of the sinusoids forming the
    2413     //best-fit curves (a0,b0,a1,b1,...), namely, a* and b* are of sine and cosine functions,
    2414     //respectively.
     2442    //best-fit curves (a0,s1,c1,s2,c2,...), where a0 is constant and s* and c* are of sine
     2443    //and cosine functions, respectively.
    24152444    std::vector<double> y;
     2445    params.clear();
    24162446    for (int i = 0; i < nDOF; ++i) {
    24172447      y.push_back(0.0);
    … …  
    24192449        y[i] += xMatrix[i][nDOF+j]*zMatrix[j];
    24202450      }
    2421     }
    2422 
    2423     double a0 = y[0];
    2424     double a1 = y[1];
    2425     double a2 = y[2];
    2426     double a3 = y[3];
    2427     params.clear();
     2451      params.push_back(y[i]);
     2452    }
    24282453
    24292454    for (int i = 0; i < nChan; ++i) {
    2430       r1[i] = a0 + a1*x1[i] + a2*x2[i] + a3*x3[i];
    2431     }
    2432     params.push_back(a0);
    2433     params.push_back(a1);
    2434     params.push_back(a2);
    2435     params.push_back(a3);
    2436 
     2455      r1[i] = y[0];
     2456      for (int j = 1; j < nDOF; ++j) {
     2457        r1[i] += y[j]*xArray[j][i];
     2458      }
     2459      residual[i] = z1[i] - r1[i];
     2460    }
    24372461
    24382462    if ((nClip == nIterClip) || (thresClip <= 0.0)) {
    24392463      break;
    24402464    } else {
    2441       std::vector<double> rz;
    24422465      double stdDev = 0.0;
    24432466      for (int i = 0; i < nChan; ++i) {
    2444         double val = abs(r1[i] - z1[i]);
    2445         rz.push_back(val);
    2446         stdDev += val*val*(double)maskArray[i];
     2467        stdDev += residual[i]*residual[i]*(double)maskArray[i];
    24472468      }
    24482469      stdDev = sqrt(stdDev/(double)nData);
    … …  
    24512472      int newNData = 0;
    24522473      for (int i = 0; i < nChan; ++i) {
    2453         if (rz[i] >= thres) {
     2474        if (abs(residual[i]) >= thres) {
    24542475          maskArray[i] = 0;
    24552476        }
    … …  
    24582479        }
    24592480      }
    2460       if (newNData == currentNData) {
    2461         break;                   //no additional flag. finish iteration.
     2481      if (newNData == nData) {
     2482        break; //no more flag to add. iteration stops.
    24622483      } else {
    2463         currentNData = newNData;
     2484        nData = newNData;
    24642485      }
    24652486    }
    … …  
    24692490  if (getResidual) {
    24702491    for (int i = 0; i < nChan; ++i) {
    2471       result.push_back((float)(z1[i] - r1[i]));
     2492      result.push_back((float)residual[i]);
    24722493    }
    24732494  } else {
    … …  
    24822503void Scantable::fitBaseline(const std::vector<bool>& mask, int whichrow, Fitter& fitter)
    24832504{
    2484   std::vector<double> abcsd = getAbcissa(whichrow);
    2485   std::vector<float> abcs;
    2486   for (uInt i = 0; i < abcsd.size(); ++i) {
    2487     abcs.push_back((float)abcsd[i]);
     2505  std::vector<double> dAbcissa = getAbcissa(whichrow);
     2506  std::vector<float> abcissa;
     2507  for (uInt i = 0; i < dAbcissa.size(); ++i) {
     2508    abcissa.push_back((float)dAbcissa[i]);
    24882509  }
    24892510  std::vector<float> spec = getSpectrum(whichrow);
    24902511
    2491   fitter.setData(abcs, spec, mask);
     2512  fitter.setData(abcissa, spec, mask);
    24922513  fitter.lfit();
    24932514}
    … …  
    25522573
    25532574/* for cspline. will be merged once cspline is available in fitter (2011/3/10 WK) */
    2554 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) {
     2575void 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) {
    25552576  if (outLogger || outTextFile) {
    25562577    float rms = getRms(chanMask, whichrow);
    25572578    String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
     2579    std::vector<bool> fixed;
     2580    fixed.clear();
    25582581
    25592582    if (outLogger) {
    … …  
    25682591
    25692592/* for sinusoid. will be merged once sinusoid is available in fitter (2011/3/10 WK) */
    2570 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) {
     2593void 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) {
    25712594  if (outLogger || outTextFile) {
    25722595    float rms = getRms(chanMask, whichrow);
    25732596    String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
     2597    std::vector<bool> fixed;
     2598    fixed.clear();
    25742599
    25752600    if (outLogger) {

  • branches/casa-prerelease/pre-asap/src/Scantable.h

    r2065 r2082  
    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);

  • branches/casa-prerelease/pre-asap/src/ScantableWrapper.h

    r2047 r2082  
    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.