Changeset 2012 for trunk/src


Ignore:
Timestamp:
02/25/11 16:51:50 (14 years ago)
Author:
WataruKawasaki
Message:

New Development: Yes

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

Ready for Test: Yes

Interface Changes: Yes

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

Test Programs:

Put in Release Notes: Yes

Module(s): sdbaseline, sd.linefinder

Description: (1) CAS-2373-related:

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

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

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

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

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

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

(2) CAS-2620-related:

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



Location:
trunk/src
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/STLineFinder.cpp

    r1670 r2012  
    925925  itsNoiseBox = in_noise_box;
    926926  itsUseMedian = in_median;
     927
     928  useScantable = true;
    927929}
    928930
    … …  
    934936  scan=in_scan.getCP();
    935937  AlwaysAssert(!scan.null(),AipsError);
    936 
     938}
     939
     940// set spectrum data to work with. this is a method to allow linefinder work
     941// without setting scantable for the purpose of using linefinder inside some
     942// method in scantable class. (Dec 22, 2010 by W.Kawasaki)
     943void STLineFinder::setData(const std::vector<float> &in_spectrum)
     944{
     945  spectrum = Vector<Float>(in_spectrum);
     946  useScantable = false;
    937947}
    938948
    … …  
    948958                const casa::uInt &whichRow) throw(casa::AipsError)
    949959{
    950   if (scan.null())
     960  if (useScantable && scan.null())
    951961      throw AipsError("STLineFinder::findLines - a scan should be set first,"
    952962                      " use set_scan");
    953963
    954   uInt nchan = scan->nchan(scan->getIF(whichRow));
     964  uInt nchan = useScantable ? scan->nchan(scan->getIF(whichRow)) : spectrum.nelements();
    955965  // set up mask and edge rejection
    956966  // no mask given...
    … …  
    962972  }
    963973  if (mask.nelements()!=nchan)
    964       throw AipsError("STLineFinder::findLines - in_scan and in_mask have different"
    965             "number of spectral channels.");
     974      throw AipsError("STLineFinder::findLines - in_scan and in_mask, or in_spectrum "
     975                      "and in_mask have different number of spectral channels.");
    966976
    967977  // taking flagged channels into account
    968   vector<bool> flaggedChannels = scan->getMask(whichRow);
    969   if (flaggedChannels.size()) {
     978  if (useScantable) {
     979    vector<bool> flaggedChannels = scan->getMask(whichRow);
     980    if (flaggedChannels.size()) {
    970981      // there is a mask set for this row
    971982      if (flaggedChannels.size() != mask.nelements()) {
    972           throw AipsError("STLineFinder::findLines - internal inconsistency: number of mask elements do not match the number of channels");
     983          throw AipsError("STLineFinder::findLines - internal inconsistency: number of "
     984                          "mask elements do not match the number of channels");
    973985      }
    974986      for (size_t ch = 0; ch<mask.nelements(); ++ch) {
    975987           mask[ch] &= flaggedChannels[ch];
    976988      }
     989    }
    977990  }
    978991
    … …  
    10151028      throw AipsError("STLineFinder::findLines - noise_box is supposed to be at least 2 elements");
    10161029
    1017   spectrum.resize();
    1018   spectrum = Vector<Float>(scan->getSpectrum(whichRow));
     1030  if (useScantable) {
     1031    spectrum.resize();
     1032    spectrum = Vector<Float>(scan->getSpectrum(whichRow));
     1033  }
    10191034
    10201035  lines.resize(0); // search from the scratch
    … …  
    11411156{
    11421157  try {
    1143        if (scan.null())
    1144            throw AipsError("STLineFinder::getMask - a scan should be set first,"
    1145                       " use set_scan followed by find_lines");
    1146        DebugAssert(mask.nelements()==scan->getChannels(last_row_used), AipsError);
    1147        /*
    1148        if (!lines.size())
    1149            throw AipsError("STLineFinder::getMask - one have to search for "
     1158    if (useScantable) {
     1159      if (scan.null())
     1160        throw AipsError("STLineFinder::getMask - a scan should be set first,"
     1161                        " use set_scan followed by find_lines");
     1162      DebugAssert(mask.nelements()==scan->getChannels(last_row_used), AipsError);
     1163    }
     1164    /*
     1165    if (!lines.size())
     1166       throw AipsError("STLineFinder::getMask - one have to search for "
    11501167                           "lines first, use find_lines");
    1151        */
    1152        std::vector<bool> res_mask(mask.nelements());
    1153        // iterator through lines
    1154        std::list<std::pair<int,int> >::const_iterator cli=lines.begin();
    1155        for (int ch=0;ch<int(res_mask.size());++ch) {
    1156             if (ch<edge.first || ch>=edge.second) res_mask[ch]=false;
    1157             else if (!mask[ch]) res_mask[ch]=false;
    1158             else {
    1159                     res_mask[ch]=!invert; // no line by default
    1160                     if (cli!=lines.end())
    1161                         if (ch>=cli->first && ch<cli->second)
    1162                              res_mask[ch]=invert; // this is a line
    1163             }
    1164             if (cli!=lines.end())
    1165                 if (ch>=cli->second) {
    1166                     ++cli; // next line in the list
    1167                 }
    1168        }
    1169        return res_mask;
     1168    */
     1169    std::vector<bool> res_mask(mask.nelements());
     1170    // iterator through lines
     1171    std::list<std::pair<int,int> >::const_iterator cli=lines.begin();
     1172    for (int ch=0;ch<int(res_mask.size());++ch) {
     1173      if (ch<edge.first || ch>=edge.second) res_mask[ch]=false;
     1174      else if (!mask[ch]) res_mask[ch]=false;
     1175      else {
     1176        res_mask[ch]=!invert; // no line by default
     1177        if (cli!=lines.end())
     1178          if (ch>=cli->first && ch<cli->second)
     1179            res_mask[ch]=invert; // this is a line
     1180      }
     1181      if (cli!=lines.end())
     1182        if (ch>=cli->second)
     1183          ++cli; // next line in the list
     1184    }
     1185    return res_mask;
    11701186  }
    11711187  catch (const AipsError &ae) {
    1172       throw;
     1188    throw;
    11731189  }
    11741190  catch (const exception &ex) {
    1175       throw AipsError(String("STLineFinder::getMask - STL error: ")+ex.what());
     1191    throw AipsError(String("STLineFinder::getMask - STL error: ")+ex.what());
    11761192  }
    11771193}
    … …  
    11821198                             const throw(casa::AipsError)
    11831199{
    1184   // convert to required abscissa units
    1185   std::vector<double> vel=scan->getAbcissa(last_row_used);
     1200  std::vector<double> vel;
     1201  if (useScantable) {
     1202    // convert to required abscissa units
     1203    vel = scan->getAbcissa(last_row_used);
     1204  } else {
     1205    for (int i = 0; i < spectrum.nelements(); ++i)
     1206      vel.push_back((double)i);
     1207  }
    11861208  std::vector<int> ranges=getLineRangesInChannels();
    11871209  std::vector<double> res(ranges.size());
    … …  
    12021224{
    12031225  try {
    1204        if (scan.null())
    1205            throw AipsError("STLineFinder::getLineRangesInChannels - a scan should be set first,"
    1206                       " use set_scan followed by find_lines");
    1207        DebugAssert(mask.nelements()==scan->getChannels(last_row_used), AipsError);
    1208 
    1209        if (!lines.size())
    1210            throw AipsError("STLineFinder::getLineRangesInChannels - one have to search for "
    1211                            "lines first, use find_lines");
    1212 
    1213        std::vector<int> res(2*lines.size());
    1214        // iterator through lines & result
    1215        std::list<std::pair<int,int> >::const_iterator cli=lines.begin();
    1216        std::vector<int>::iterator ri=res.begin();
    1217        for (;cli!=lines.end() && ri!=res.end();++cli,++ri) {
    1218             *ri=cli->first;
    1219             if (++ri!=res.end())
    1220                 *ri=cli->second-1;
    1221        }
    1222        return res;
    1223   }
    1224   catch (const AipsError &ae) {
    1225       throw;
    1226   }
    1227   catch (const exception &ex) {
    1228       throw AipsError(String("STLineFinder::getLineRanges - STL error: ")+ex.what());
     1226    if (useScantable) {
     1227      if (scan.null())
     1228        throw AipsError("STLineFinder::getLineRangesInChannels - a scan should be set first,"
     1229                        " use set_scan followed by find_lines");
     1230      DebugAssert(mask.nelements()==scan->getChannels(last_row_used), AipsError);
     1231    }
     1232
     1233    if (!lines.size())
     1234      throw AipsError("STLineFinder::getLineRangesInChannels - one have to search for "
     1235                      "lines first, use find_lines");
     1236
     1237    std::vector<int> res(2*lines.size());
     1238    // iterator through lines & result
     1239    std::list<std::pair<int,int> >::const_iterator cli = lines.begin();
     1240    std::vector<int>::iterator ri = res.begin();
     1241    for (; cli != lines.end() && ri != res.end(); ++cli,++ri) {
     1242      *ri = cli->first;
     1243      if (++ri != res.end())
     1244        *ri = cli->second - 1;
     1245    }
     1246    return res;
     1247  } catch (const AipsError &ae) {
     1248    throw;
     1249  } catch (const exception &ex) {
     1250    throw AipsError(String("STLineFinder::getLineRanges - STL error: ") + ex.what());
    12291251  }
    12301252}

  • trunk/src/STLineFinder.h

    r1644 r2012  
    164164   // set the scan to work with (in_scan parameter)
    165165   void setScan(const ScantableWrapper &in_scan) throw(casa::AipsError);
     166
     167   // set spectrum data to work with. this is a method to allow linefinder work
     168   // without setting scantable for the purpose of using linefinder inside some
     169   // method in scantable class. (Dec. 22, 2010 by W.Kawasaki)
     170   void setData(const std::vector<float> &in_spectrum);
    166171
    167172   // search for spectral lines in a row specified by whichRow
    … …  
    257262   // the lowest 80% of deviations (default)
    258263   casa::Bool itsUseMedian;
     264
     265   // true if spectra and mask data should be provided from
     266   // scantable (default = true)
     267   bool useScantable;
    259268};
    260269

  • trunk/src/Scantable.cpp

    r2005 r2012  
    6767#include "STPolStokes.h"
    6868#include "STAttr.h"
     69#include "STLineFinder.h"
    6970#include "MathUtils.h"
    7071
    … …  
    12171218}
    12181219**/
     1220
    12191221void Scantable::setRestFrequencies( const vector<std::string>& name )
    12201222{
    … …  
    17491751  Vector<uChar> flags;
    17501752  flagsCol_.get(uInt(whichrow), flags);
    1751   for (int i = 0; i < flags.size(); i++) {
    1752     if (i==0) {
    1753       flag = flags[i];
    1754     }
    1755     else {
    1756       flag &= flags[i];
    1757     }
    1758     return ((flag >> 7) == 1);
    1759    }
    1760 }
    1761 
    1762 void Scantable::doPolyBaseline(const std::vector<bool>& mask, int order, int rowno, Fitter& fitter)
    1763 {
     1753  flag = flags[0];
     1754  for (int i = 1; i < flags.size(); ++i) {
     1755    flag &= flags[i];
     1756  }
     1757  return ((flag >> 7) == 1);
     1758}
     1759
     1760void Scantable::cubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, bool outLogger, const std::string& blfile) {
     1761  ofstream ofs;
     1762  String coordInfo;
     1763  bool hasSameNchan;
     1764  int firstIF;
     1765  bool outTextFile = false;
     1766
     1767  if (blfile != "") {
     1768    ofs.open(blfile.c_str(), ios::out | ios::app);
     1769    if (ofs) outTextFile = true;
     1770  }
     1771
     1772  if (outLogger || outTextFile) {
     1773    coordInfo = getCoordInfo()[0];
     1774    if (coordInfo == "") coordInfo = "channel";
     1775    hasSameNchan = hasSameNchanOverIFs();
     1776    firstIF = getIF(0);
     1777  }
     1778
     1779  //Fitter fitter = Fitter();
     1780  //fitter.setExpression("cspline", nPiece, thresClip, nIterClip);
     1781
     1782  int nRow = nrow();
     1783  std::vector<bool> chanMask;
     1784
     1785  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
     1786
     1787    chanMask = getCompositeChanMask(whichrow, mask);
     1788    //fitBaseline(chanMask, whichrow, fitter);
     1789    //setSpectrum(fitter.getResidual(), whichrow);
     1790    std::vector<int> pieceRanges;
     1791    std::vector<float> params;
     1792    std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, pieceRanges, params, nPiece, thresClip, nIterClip);
     1793    setSpectrum(res, whichrow);
     1794
     1795    if (outLogger || outTextFile) {
     1796      std::vector<bool>  fixed;
     1797      fixed.clear();
     1798      float rms = getRms(chanMask, whichrow);
     1799      String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan, firstIF, true);
     1800
     1801      if (outLogger) {
     1802        LogIO ols(LogOrigin("Scantable", "cubicSplineBaseline()", WHERE));
     1803        ols << formatPiecewiseBaselineParams(pieceRanges, params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
     1804      }
     1805      if (outTextFile) {
     1806        ofs << formatPiecewiseBaselineParams(pieceRanges, params, fixed, rms, masklist, whichrow, true) << flush;
     1807      }
     1808    }
     1809
     1810  }
     1811
     1812  if (outTextFile) ofs.close();
     1813}
     1814
     1815
     1816void Scantable::autoCubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool outLogger, const std::string& blfile)
     1817{
     1818  ofstream ofs;
     1819  String coordInfo;
     1820  bool hasSameNchan;
     1821  int firstIF;
     1822  bool outTextFile = false;
     1823
     1824  if (blfile != "") {
     1825    ofs.open(blfile.c_str(), ios::out | ios::app);
     1826    if (ofs) outTextFile = true;
     1827  }
     1828
     1829  if (outLogger || outTextFile) {
     1830    coordInfo = getCoordInfo()[0];
     1831    if (coordInfo == "") coordInfo = "channel";
     1832    hasSameNchan = hasSameNchanOverIFs();
     1833    firstIF = getIF(0);
     1834  }
     1835
     1836  //Fitter fitter = Fitter();
     1837  //fitter.setExpression("cspline", nPiece, thresClip, nIterClip);
     1838
     1839  int nRow = nrow();
     1840  std::vector<bool> chanMask;
     1841  int minEdgeSize = getIFNos().size()*2;
     1842  STLineFinder lineFinder = STLineFinder();
     1843  lineFinder.setOptions(threshold, 3, chanAvgLimit);
     1844
     1845  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
     1846
     1847    //-------------------------------------------------------
     1848    //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
     1849    //-------------------------------------------------------
     1850    int edgeSize = edge.size();
     1851    std::vector<int> currentEdge;
     1852    if (edgeSize >= 2) {
     1853      int idx = 0;
     1854      if (edgeSize > 2) {
     1855        if (edgeSize < minEdgeSize) {
     1856          throw(AipsError("Length of edge element info is less than that of IFs"));
     1857        }
     1858        idx = 2 * getIF(whichrow);
     1859      }
     1860      currentEdge.push_back(edge[idx]);
     1861      currentEdge.push_back(edge[idx+1]);
     1862    } else {
     1863      throw(AipsError("Wrong length of edge element"));
     1864    }
     1865    lineFinder.setData(getSpectrum(whichrow));
     1866    lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
     1867    chanMask = lineFinder.getMask();
     1868    //-------------------------------------------------------
     1869
     1870
     1871    //fitBaseline(chanMask, whichrow, fitter);
     1872    //setSpectrum(fitter.getResidual(), whichrow);
     1873    std::vector<int> pieceRanges;
     1874    std::vector<float> params;
     1875    std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, pieceRanges, params, nPiece, thresClip, nIterClip);
     1876    setSpectrum(res, whichrow);
     1877
     1878    if (outLogger || outTextFile) {
     1879      std::vector<bool> fixed;
     1880      fixed.clear();
     1881      float rms = getRms(chanMask, whichrow);
     1882      String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan, firstIF, true);
     1883
     1884      if (outLogger) {
     1885        LogIO ols(LogOrigin("Scantable", "autoCubicSplineBaseline()", WHERE));
     1886        ols << formatPiecewiseBaselineParams(pieceRanges, params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
     1887      }
     1888      if (outTextFile) {
     1889        ofs << formatPiecewiseBaselineParams(pieceRanges, params, fixed, rms, masklist, whichrow, true) << flush;
     1890      }
     1891    }
     1892
     1893  }
     1894
     1895  if (outTextFile) ofs.close();
     1896}
     1897
     1898
     1899std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data, const std::vector<bool>& mask, std::vector<int>& sectionRanges, std::vector<float>& params, int nPiece, float thresClip, int nIterClip) {
     1900  if (nPiece < 1) {
     1901    throw(AipsError("wrong number of the sections for fitting"));
     1902  }
     1903  if (data.size() != mask.size()) {
     1904    throw(AipsError("data and mask have different size"));
     1905  }
     1906
     1907  int nChan = data.size();
     1908  std::vector<int> maskArray;
     1909  std::vector<int> x;
     1910  for (int i = 0; i < nChan; ++i) {
     1911    maskArray.push_back(mask[i] ? 1 : 0);
     1912    if (mask[i]) {
     1913      x.push_back(i);
     1914    }
     1915  }
     1916
     1917  int nData = x.size();
     1918  int nElement = (int)(floor(floor((double)(nData/nPiece))+0.5));
     1919
     1920  std::vector<double> sectionList0, sectionList1, sectionListR;
     1921  sectionList0.push_back((double)x[0]);
     1922  sectionRanges.clear();
     1923  sectionRanges.push_back(x[0]);
     1924  for (int i = 1; i < nPiece; ++i) {
     1925    double valX = (double)x[nElement*i];
     1926    sectionList0.push_back(valX);
     1927    sectionList1.push_back(valX);
     1928    sectionListR.push_back(1.0/valX);
     1929
     1930    sectionRanges.push_back(x[nElement*i]-1);
     1931    sectionRanges.push_back(x[nElement*i]);
     1932  }
     1933  sectionList1.push_back((double)(x[x.size()-1]+1));
     1934  sectionRanges.push_back(x[x.size()-1]);
     1935
     1936  std::vector<double> x1, x2, x3, z1, x1z1, x2z1, x3z1, r1;
     1937  for (int i = 0; i < nChan; ++i) {
     1938    x1.push_back((double)i);
     1939    x2.push_back((double)(i*i));
     1940    x3.push_back((double)(i*i*i));
     1941    z1.push_back((double)data[i]);
     1942    x1z1.push_back(((double)i)*(double)data[i]);
     1943    x2z1.push_back(((double)(i*i))*(double)data[i]);
     1944    x3z1.push_back(((double)(i*i*i))*(double)data[i]);
     1945    r1.push_back(0.0);
     1946  }
     1947
     1948  int currentNData = nData;
     1949  int nDOF = nPiece + 3;  //number of parameters to solve, namely, 4+(nPiece-1).
     1950
     1951  for (int nClip = 0; nClip < nIterClip+1; ++nClip) {
     1952    //xMatrix : horizontal concatenation of the least-sq. matrix (left) and an
     1953    //identity matrix (right).
     1954    //the right part is used to calculate the inverse matrix of the left part.
     1955    double xMatrix[nDOF][2*nDOF];
     1956    double zMatrix[nDOF];
     1957    for (int i = 0; i < nDOF; ++i) {
     1958      for (int j = 0; j < 2*nDOF; ++j) {
     1959        xMatrix[i][j] = 0.0;
     1960      }
     1961      xMatrix[i][nDOF+i] = 1.0;
     1962      zMatrix[i] = 0.0;
     1963    }
     1964
     1965    for (int n = 0; n < nPiece; ++n) {
     1966      for (int i = sectionList0[n]; i < sectionList1[n]; ++i) {
     1967        if (maskArray[i] == 0) continue;
     1968        xMatrix[0][0] += 1.0;
     1969        xMatrix[0][1] += x1[i];
     1970        xMatrix[0][2] += x2[i];
     1971        xMatrix[0][3] += x3[i];
     1972        xMatrix[1][1] += x2[i];
     1973        xMatrix[1][2] += x3[i];
     1974        xMatrix[1][3] += x2[i]*x2[i];
     1975        xMatrix[2][2] += x2[i]*x2[i];
     1976        xMatrix[2][3] += x3[i]*x2[i];
     1977        xMatrix[3][3] += x3[i]*x3[i];
     1978        zMatrix[0] += z1[i];
     1979        zMatrix[1] += x1z1[i];
     1980        zMatrix[2] += x2z1[i];
     1981        zMatrix[3] += x3z1[i];
     1982        for (int j = 0; j < n; ++j) {
     1983          double q = 1.0 - x1[i]*sectionListR[j];
     1984          q = q*q*q;
     1985          xMatrix[0][j+4] += q;
     1986          xMatrix[1][j+4] += q*x1[i];
     1987          xMatrix[2][j+4] += q*x2[i];
     1988          xMatrix[3][j+4] += q*x3[i];
     1989          for (int k = 0; k < j; ++k) {
     1990            double r = 1.0 - x1[i]*sectionListR[k];
     1991            r = r*r*r;
     1992            xMatrix[k+4][j+4] += r*q;
     1993          }
     1994          xMatrix[j+4][j+4] += q*q;
     1995          zMatrix[j+4] += q*z1[i];
     1996        }
     1997      }
     1998    }
     1999
     2000    for (int i = 0; i < nDOF; ++i) {
     2001      for (int j = 0; j < i; ++j) {
     2002        xMatrix[i][j] = xMatrix[j][i];
     2003      }
     2004    }
     2005
     2006    std::vector<double> invDiag;
     2007    for (int i = 0; i < nDOF; ++i) {
     2008      invDiag.push_back(1.0/xMatrix[i][i]);
     2009      for (int j = 0; j < nDOF; ++j) {
     2010        xMatrix[i][j] *= invDiag[i];
     2011      }
     2012    }
     2013
     2014    for (int k = 0; k < nDOF; ++k) {
     2015      for (int i = 0; i < nDOF; ++i) {
     2016        if (i != k) {
     2017          double factor1 = xMatrix[k][k];
     2018          double factor2 = xMatrix[i][k];
     2019          for (int j = k; j < 2*nDOF; ++j) {
     2020            xMatrix[i][j] *= factor1;
     2021            xMatrix[i][j] -= xMatrix[k][j]*factor2;
     2022            xMatrix[i][j] /= factor1;
     2023          }
     2024        }
     2025      }
     2026      double xDiag = xMatrix[k][k];
     2027      for (int j = k; j < 2*nDOF; ++j) {
     2028        xMatrix[k][j] /= xDiag;
     2029      }
     2030    }
     2031   
     2032    for (int i = 0; i < nDOF; ++i) {
     2033      for (int j = 0; j < nDOF; ++j) {
     2034        xMatrix[i][nDOF+j] *= invDiag[j];
     2035      }
     2036    }
     2037    //compute a vector y which consists of the coefficients of the best-fit spline curves
     2038    //(a0,a1,a2,a3(,b3,c3,...)), namely, the ones for the leftmost piece and the ones of
     2039    //cubic terms for the other pieces (in case nPiece>1).
     2040    std::vector<double> y;
     2041    for (int i = 0; i < nDOF; ++i) {
     2042      y.push_back(0.0);
     2043      for (int j = 0; j < nDOF; ++j) {
     2044        y[i] += xMatrix[i][nDOF+j]*zMatrix[j];
     2045      }
     2046    }
     2047
     2048    double a0 = y[0];
     2049    double a1 = y[1];
     2050    double a2 = y[2];
     2051    double a3 = y[3];
     2052    params.clear();
     2053
     2054    for (int n = 0; n < nPiece; ++n) {
     2055      for (int i = sectionList0[n]; i < sectionList1[n]; ++i) {
     2056        r1[i] = a0 + a1*x1[i] + a2*x2[i] + a3*x3[i];
     2057      }
     2058      params.push_back(a0);
     2059      params.push_back(a1);
     2060      params.push_back(a2);
     2061      params.push_back(a3);
     2062
     2063      if (n == nPiece-1) break;
     2064
     2065      double d = y[4+n];
     2066      a0 += d;
     2067      a1 -= 3.0*d*sectionListR[n];
     2068      a2 += 3.0*d*sectionListR[n]*sectionListR[n];
     2069      a3 -= d*sectionListR[n]*sectionListR[n]*sectionListR[n];
     2070    }
     2071
     2072    if ((nClip == nIterClip) || (thresClip <= 0.0)) {
     2073      break;
     2074    } else {
     2075      std::vector<double> rz;
     2076      double stdDev = 0.0;
     2077      for (int i = 0; i < nChan; ++i) {
     2078        double val = abs(r1[i] - z1[i]);
     2079        rz.push_back(val);
     2080        stdDev += val*val*(double)maskArray[i];
     2081      }
     2082      stdDev = sqrt(stdDev/(double)nData);
     2083     
     2084      double thres = stdDev * thresClip;
     2085      int newNData = 0;
     2086      for (int i = 0; i < nChan; ++i) {
     2087        if (rz[i] >= thres) {
     2088          maskArray[i] = 0;
     2089        }
     2090        if (maskArray[i] > 0) {
     2091          newNData++;
     2092        }
     2093      }
     2094      if (newNData == currentNData) {
     2095        break;                   //no additional flag. finish iteration.
     2096      } else {
     2097        currentNData = newNData;
     2098      }
     2099    }
     2100  }
     2101
     2102  std::vector<float> residual;
     2103  for (int i = 0; i < nChan; ++i) {
     2104    residual.push_back((float)(z1[i] - r1[i]));
     2105  }
     2106  return residual;
     2107
     2108}
     2109
     2110void Scantable::fitBaseline(const std::vector<bool>& mask, int whichrow, Fitter& fitter)
     2111{
     2112  std::vector<double> abcsd = getAbcissa(whichrow);
     2113  std::vector<float> abcs;
     2114  for (int i = 0; i < abcsd.size(); ++i) {
     2115    abcs.push_back((float)abcsd[i]);
     2116  }
     2117  std::vector<float> spec = getSpectrum(whichrow);
     2118
     2119  fitter.setData(abcs, spec, mask);
     2120  fitter.lfit();
     2121}
     2122
     2123std::vector<bool> Scantable::getCompositeChanMask(int whichrow, const std::vector<bool>& inMask)
     2124{
     2125  std::vector<bool> chanMask = getMask(whichrow);
     2126  int chanMaskSize = chanMask.size();
     2127  if (chanMaskSize != inMask.size()) {
     2128    throw(AipsError("different mask sizes"));
     2129  }
     2130  for (int i = 0; i < chanMaskSize; ++i) {
     2131    chanMask[i] = chanMask[i] && inMask[i];
     2132  }
     2133
     2134  return chanMask;
     2135}
     2136
     2137/*
     2138std::vector<bool> Scantable::getCompositeChanMask(int whichrow, const std::vector<bool>& inMask, const std::vector<int>& edge, const int minEdgeSize, STLineFinder& lineFinder)
     2139{
     2140  int edgeSize = edge.size();
     2141  std::vector<int> currentEdge;
     2142  if (edgeSize >= 2) {
     2143      int idx = 0;
     2144      if (edgeSize > 2) {
     2145        if (edgeSize < minEdgeSize) {
     2146          throw(AipsError("Length of edge element info is less than that of IFs"));
     2147        }
     2148        idx = 2 * getIF(whichrow);
     2149      }
     2150      currentEdge.push_back(edge[idx]);
     2151      currentEdge.push_back(edge[idx+1]);
     2152  } else {
     2153    throw(AipsError("Wrong length of edge element"));
     2154  }
     2155
     2156  lineFinder.setData(getSpectrum(whichrow));
     2157  lineFinder.findLines(getCompositeChanMask(whichrow, inMask), currentEdge, whichrow);
     2158
     2159  return lineFinder.getMask();
     2160}
     2161*/
     2162
     2163void Scantable::polyBaseline(const std::vector<bool>& mask, int order, bool outLogger, const std::string& blfile)
     2164{
     2165  ofstream ofs;
     2166  String coordInfo;
     2167  bool hasSameNchan;
     2168  int firstIF;
     2169  bool outTextFile = false;
     2170
     2171  if (blfile != "") {
     2172    ofs.open(blfile.c_str(), ios::out | ios::app);
     2173    if (ofs) outTextFile = true;
     2174  }
     2175
     2176  if (outLogger || outTextFile) {
     2177    coordInfo = getCoordInfo()[0];
     2178    if (coordInfo == "") coordInfo = "channel";
     2179    hasSameNchan = hasSameNchanOverIFs();
     2180    firstIF = getIF(0);
     2181  }
     2182
     2183  Fitter fitter = Fitter();
    17642184  fitter.setExpression("poly", order);
    17652185
    1766   std::vector<double> abcsd = getAbcissa(rowno);
    1767   std::vector<float> abcs;
    1768   for (int i = 0; i < abcsd.size(); i++) {
    1769     abcs.push_back((float)abcsd[i]);
    1770   }
    1771   std::vector<float> spec = getSpectrum(rowno);
     2186  int nRow = nrow();
     2187  std::vector<bool> chanMask;
     2188
     2189  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
     2190
     2191    chanMask = getCompositeChanMask(whichrow, mask);
     2192    fitBaseline(chanMask, whichrow, fitter);
     2193    setSpectrum(fitter.getResidual(), whichrow);
     2194   
     2195    if (outLogger || outTextFile) {
     2196      std::vector<float> params = fitter.getParameters();
     2197      std::vector<bool>  fixed  = fitter.getFixedParameters();
     2198      float rms = getRms(chanMask, whichrow);
     2199      String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan, firstIF, true);
     2200
     2201      if (outLogger) {
     2202        LogIO ols(LogOrigin("Scantable", "polyBaseline()", WHERE));
     2203        ols << formatBaselineParams(params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
     2204      }
     2205      if (outTextFile) {
     2206        ofs << formatBaselineParams(params, fixed, rms, masklist, whichrow, true) << flush;
     2207      }
     2208    }
     2209
     2210  }
     2211
     2212  if (outTextFile) ofs.close();
     2213}
     2214
     2215void Scantable::autoPolyBaseline(const std::vector<bool>& mask, int order, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool outLogger, const std::string& blfile)
     2216{
     2217  ofstream ofs;
     2218  String coordInfo;
     2219  bool hasSameNchan;
     2220  int firstIF;
     2221  bool outTextFile = false;
     2222
     2223  if (blfile != "") {
     2224    ofs.open(blfile.c_str(), ios::out | ios::app);
     2225    if (ofs) outTextFile = true;
     2226  }
     2227
     2228  if (outLogger || outTextFile) {
     2229    coordInfo = getCoordInfo()[0];
     2230    if (coordInfo == "") coordInfo = "channel";
     2231    hasSameNchan = hasSameNchanOverIFs();
     2232    firstIF = getIF(0);
     2233  }
     2234
     2235  Fitter fitter = Fitter();
     2236  fitter.setExpression("poly", order);
     2237
     2238  int nRow = nrow();
     2239  std::vector<bool> chanMask;
     2240  int minEdgeSize = getIFNos().size()*2;
     2241  STLineFinder lineFinder = STLineFinder();
     2242  lineFinder.setOptions(threshold, 3, chanAvgLimit);
     2243
     2244  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
     2245
     2246    //-------------------------------------------------------
     2247    //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
     2248    //-------------------------------------------------------
     2249    int edgeSize = edge.size();
     2250    std::vector<int> currentEdge;
     2251    if (edgeSize >= 2) {
     2252      int idx = 0;
     2253      if (edgeSize > 2) {
     2254        if (edgeSize < minEdgeSize) {
     2255          throw(AipsError("Length of edge element info is less than that of IFs"));
     2256        }
     2257        idx = 2 * getIF(whichrow);
     2258      }
     2259      currentEdge.push_back(edge[idx]);
     2260      currentEdge.push_back(edge[idx+1]);
     2261    } else {
     2262      throw(AipsError("Wrong length of edge element"));
     2263    }
     2264    lineFinder.setData(getSpectrum(whichrow));
     2265    lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
     2266    chanMask = lineFinder.getMask();
     2267    //-------------------------------------------------------
     2268
     2269
     2270    fitBaseline(chanMask, whichrow, fitter);
     2271    setSpectrum(fitter.getResidual(), whichrow);
     2272
     2273    if (outLogger || outTextFile) {
     2274      std::vector<float> params = fitter.getParameters();
     2275      std::vector<bool>  fixed  = fitter.getFixedParameters();
     2276      float rms = getRms(chanMask, whichrow);
     2277      String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan, firstIF, true);
     2278
     2279      if (outLogger) {
     2280        LogIO ols(LogOrigin("Scantable", "autoPolyBaseline()", WHERE));
     2281        ols << formatBaselineParams(params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
     2282      }
     2283      if (outTextFile) {
     2284        ofs << formatBaselineParams(params, fixed, rms, masklist, whichrow, true) << flush;
     2285      }
     2286    }
     2287
     2288  }
     2289
     2290  if (outTextFile) ofs.close();
     2291}
     2292
     2293
     2294float Scantable::getRms(const std::vector<bool>& mask, int whichrow) {
     2295  Vector<Float> spec;
     2296  specCol_.get(whichrow, spec);
     2297
     2298  float mean = 0.0;
     2299  float smean = 0.0;
     2300  int n = 0;
     2301  for (int i = 0; i < spec.nelements(); ++i) {
     2302    if (mask[i]) {
     2303      mean += spec[i];
     2304      smean += spec[i]*spec[i];
     2305      n++;
     2306    }
     2307  }
     2308
     2309  mean /= (float)n;
     2310  smean /= (float)n;
     2311
     2312  return sqrt(smean - mean*mean);
     2313}
     2314
     2315
     2316std::string Scantable::formatBaselineParamsHeader(int whichrow, const std::string& masklist, bool verbose) const
     2317{
     2318  ostringstream oss;
     2319
     2320  if (verbose) {
     2321    for (int i = 0; i < 60; ++i) {
     2322      oss << "-";
     2323    }
     2324    oss << endl;
     2325    oss <<  " Scan[" << getScan(whichrow)  << "]";
     2326    oss <<  " Beam[" << getBeam(whichrow)  << "]";
     2327    oss <<    " IF[" << getIF(whichrow)    << "]";
     2328    oss <<   " Pol[" << getPol(whichrow)   << "]";
     2329    oss << " Cycle[" << getCycle(whichrow) << "]: " << endl;
     2330    oss << "Fitter range = " << masklist << endl;
     2331    oss << "Baseline parameters" << endl;
     2332    oss << flush;
     2333  }
     2334
     2335  return String(oss);
     2336}
     2337
     2338std::string Scantable::formatBaselineParamsFooter(float rms, bool verbose) const
     2339{
     2340  ostringstream oss;
     2341
     2342  if (verbose) {
     2343    oss << endl;
     2344    oss << "Results of baseline fit" << endl;
     2345    oss << "  rms = " << setprecision(6) << rms << endl;
     2346  }
     2347
     2348  return String(oss);
     2349}
     2350
     2351std::string Scantable::formatPiecewiseBaselineParams(const std::vector<int>& ranges, const std::vector<float>& params, const std::vector<bool>& fixed, float rms, const std::string& masklist, int whichrow, bool verbose) const
     2352{
     2353  ostringstream oss;
     2354  oss << formatBaselineParamsHeader(whichrow, masklist, verbose);
     2355
     2356  if ((params.size() > 0) && (ranges.size() > 0)) {
     2357    if (((ranges.size() % 2) == 0) && ((params.size() % (ranges.size() / 2)) == 0)) {
     2358      int nParam = params.size() / (ranges.size() / 2);
     2359      for (uInt j = 0; j < ranges.size(); j+=2) {
     2360        oss << "[" << ranges[j] << "," << ranges[j+1] << "]";
     2361        for (uInt i = 0; i < nParam; ++i) {
     2362          if (i > 0) {
     2363            oss << ",";
     2364          }
     2365          String fix = ((fixed.size() > 0) && (fixed[i]) && verbose) ? "(fixed)" : "";
     2366          oss << "  p" << i << fix << "= " << setprecision(6) << params[j/2*nParam+i];
     2367        }
     2368      }
     2369    } else {
     2370      oss << "  ";
     2371    }
     2372  } else {
     2373    oss << "  Not fitted";
     2374  }
     2375
     2376  oss << formatBaselineParamsFooter(rms, verbose);
     2377  oss << flush;
     2378
     2379  return String(oss);
     2380}
     2381
     2382std::string Scantable::formatBaselineParams(const std::vector<float>& params, const std::vector<bool>& fixed, float rms, const std::string& masklist, int whichrow, bool verbose) const
     2383{
     2384  ostringstream oss;
     2385  oss << formatBaselineParamsHeader(whichrow, masklist, verbose);
     2386
     2387  if (params.size() > 0) {
     2388    for (uInt i = 0; i < params.size(); ++i) {
     2389      if (i > 0) {
     2390        oss << ",";
     2391      }
     2392      String fix = ((fixed.size() > 0) && (fixed[i]) && verbose) ? "(fixed)" : "";
     2393      oss << "  p" << i << fix << "= " << setprecision(6) << params[i];
     2394    }
     2395  } else {
     2396    oss << "  Not fitted";
     2397  }
     2398
     2399  oss << formatBaselineParamsFooter(rms, verbose);
     2400  oss << flush;
     2401
     2402  return String(oss);
     2403}
     2404
     2405std::string Scantable::getMaskRangeList(const std::vector<bool>& mask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, int firstIF, bool silent)
     2406{
     2407  if (mask.size() < 2) {
     2408    throw(AipsError("The mask elements should be > 1"));
     2409  }
     2410  if (mask.size() != nchan()) {
     2411    throw(AipsError("Number of channels in scantable != number of mask elements"));
     2412  }
     2413
     2414  std::vector<int> startIdx = getMaskEdgeIndices(mask, true);
     2415  std::vector<int> endIdx   = getMaskEdgeIndices(mask, false);
     2416
     2417  if (startIdx.size() != endIdx.size()) {
     2418    throw(AipsError("Numbers of mask start and mask end are not identical"));
     2419  }
     2420  for (int i = 0; i < startIdx.size(); ++i) {
     2421    if (startIdx[i] > endIdx[i]) {
     2422      throw(AipsError("Mask start index > mask end index"));
     2423    }
     2424  }
     2425
     2426  if (!silent) {
     2427    LogIO logOs(LogOrigin("Scantable", "getMaskRangeList()", WHERE));
     2428    logOs << LogIO::WARN << "The current mask window unit is " << coordInfo;
     2429    if (!hasSameNchan) {
     2430      logOs << endl << "This mask is only valid for IF=" << firstIF;
     2431    }
     2432    logOs << LogIO::POST;
     2433  }
     2434
     2435  std::vector<double> abcissa = getAbcissa(whichrow);
     2436  ostringstream oss;
     2437  oss.setf(ios::fixed);
     2438  oss << setprecision(1) << "[";
     2439  for (int i = 0; i < startIdx.size(); ++i) {
     2440    std::vector<float> aMaskRange;
     2441    aMaskRange.push_back((float)abcissa[startIdx[i]]);
     2442    aMaskRange.push_back((float)abcissa[endIdx[i]]);
     2443    sort(aMaskRange.begin(), aMaskRange.end());
     2444    if (i > 0) oss << ",";
     2445    oss << "[" << aMaskRange[0] << "," << aMaskRange[1] << "]";
     2446  }
     2447  oss << "]" << flush;
     2448
     2449  return String(oss);
     2450}
     2451
     2452bool Scantable::hasSameNchanOverIFs()
     2453{
     2454  int nIF = nif(-1);
     2455  int nCh;
     2456  int totalPositiveNChan = 0;
     2457  int nPositiveNChan = 0;
     2458
     2459  for (int i = 0; i < nIF; ++i) {
     2460    nCh = nchan(i);
     2461    if (nCh > 0) {
     2462      totalPositiveNChan += nCh;
     2463      nPositiveNChan++;
     2464    }
     2465  }
     2466
     2467  return (totalPositiveNChan == (nPositiveNChan * nchan(0)));
     2468}
     2469
     2470std::vector<int> Scantable::getMaskEdgeIndices(const std::vector<bool>& mask, bool getStartIndices)
     2471{
     2472  if (mask.size() < 2) {
     2473    throw(AipsError("The mask elements should be > 1"));
     2474  }
     2475  if (mask.size() != nchan()) {
     2476    throw(AipsError("Number of channels in scantable != number of mask elements"));
     2477  }
     2478
     2479  std::vector<int> out;
     2480  int endIdx = mask.size() - 1;
     2481
     2482  if (getStartIndices) {
     2483    if (mask[0]) {
     2484      out.push_back(0);
     2485    }
     2486    for (int i = 0; i < endIdx; ++i) {
     2487      if ((!mask[i]) && mask[i+1]) {
     2488        out.push_back(i+1);
     2489      }
     2490    }
     2491  } else {
     2492    for (int i = 0; i < endIdx; ++i) {
     2493      if (mask[i] && (!mask[i+1])) {
     2494        out.push_back(i);
     2495      }
     2496    }
     2497    if (mask[endIdx]) {
     2498      out.push_back(endIdx);
     2499    }
     2500  }
     2501
     2502  return out;
     2503}
     2504
     2505
     2506/*
     2507STFitEntry Scantable::polyBaseline(const std::vector<bool>& mask, int order, int rowno)
     2508{
     2509  Fitter fitter = Fitter();
     2510  fitter.setExpression("poly", order);
     2511
    17722512  std::vector<bool> fmask = getMask(rowno);
    17732513  if (fmask.size() != mask.size()) {
    17742514    throw(AipsError("different mask sizes"));
    17752515  }
    1776   for (int i = 0; i < fmask.size(); i++) {
     2516  for (int i = 0; i < fmask.size(); ++i) {
    17772517    fmask[i] = fmask[i] && mask[i];
    17782518  }
    1779   fitter.setData(abcs, spec, fmask);
    1780 
    1781   fitter.lfit();
    1782 }
    1783 
    1784 void Scantable::polyBaselineBatch(const std::vector<bool>& mask, int order)
    1785 {
    1786   Fitter fitter = Fitter();
    1787   for (uInt rowno=0; rowno < nrow(); ++rowno) {
    1788     doPolyBaseline(mask, order, rowno, fitter);
    1789     setSpectrum(fitter.getResidual(), rowno);
    1790   }
    1791 }
    1792 
    1793 STFitEntry Scantable::polyBaseline(const std::vector<bool>& mask, int order, int rowno)
    1794 {
    1795   Fitter fitter = Fitter();
    1796   doPolyBaseline(mask, order, rowno, fitter);
     2519
     2520  fitBaseline(fmask, rowno, fitter);
    17972521  setSpectrum(fitter.getResidual(), rowno);
    17982522  return fitter.getFitEntry();
    17992523}
     2524*/
    18002525
    18012526}

  • trunk/src/Scantable.h

    r2005 r2012  
    1616#include <string>
    1717#include <vector>
     18#include <iostream>
     19#include <fstream>
    1820// AIPS++
    1921#include <casa/aips.h>
    … …  
    493495  bool getFlagtraFast(int whichrow);
    494496
    495   void polyBaselineBatch(const std::vector<bool>& mask, int order);
    496   STFitEntry polyBaseline(const std::vector<bool>& mask, int order, int rowno);
     497  void polyBaseline(const std::vector<bool>& mask,
     498                    int order,
     499                    bool outLogger=false,
     500                    const std::string& blfile="");
     501  void autoPolyBaseline(const std::vector<bool>& mask,
     502                        int order,
     503                        const std::vector<int>& edge,
     504                        float threshold=3.0,
     505                        int chanAvgLimit=1,
     506                        bool outLogger=false,
     507                        const std::string& blfile="");
     508  void cubicSplineBaseline(const std::vector<bool>& mask,
     509                           int nPiece,
     510                           float thresClip,
     511                           int nIterClip,
     512                           bool outLogger=false,
     513                           const std::string& blfile="");
     514  void autoCubicSplineBaseline(const std::vector<bool>& mask,
     515                               int nPiece,
     516                               float thresClip,
     517                               int nIterClip,
     518                               const std::vector<int>& edge,
     519                               float threshold=3.0,
     520                               int chanAvgLimit=1,
     521                               bool outLogger=false,
     522                               const std::string& blfile="");
     523  float getRms(const std::vector<bool>& mask, int whichrow);
     524  std::string formatBaselineParams(const std::vector<float>& params,
     525                                   const std::vector<bool>& fixed,
     526                                   float rms,
     527                                   const std::string& masklist,
     528                                   int whichrow,
     529                                   bool verbose=false) const;
     530  std::string formatPiecewiseBaselineParams(const std::vector<int>& ranges,
     531                                            const std::vector<float>& params,
     532                                            const std::vector<bool>& fixed,
     533                                            float rms,
     534                                            const std::string& masklist,
     535                                            int whichrow,
     536                                            bool verbose=false) const;
     537
    497538
    498539private:
    … …  
    608649                                                      const casa::Array<T2>&);
    609650
    610   void doPolyBaseline(const std::vector<bool>& mask, int order, int rowno, Fitter& fitter);
     651  void fitBaseline(const std::vector<bool>& mask, int whichrow, Fitter& fitter);
     652  std::vector<float> doCubicSplineFitting(const std::vector<float>& data,
     653                                          const std::vector<bool>& mask,
     654                                          std::vector<int>& sectionRanges,
     655                                          std::vector<float>& params,
     656                                          int nPiece=2,
     657                                          float thresClip=3.0,
     658                                          int nIterClip=1);
     659  bool hasSameNchanOverIFs();
     660  std::string getMaskRangeList(const std::vector<bool>& mask,
     661                                int whichrow,
     662                                const casa::String& coordInfo,
     663                                bool hasSameNchan,
     664                                int firstIF,
     665                                bool silent=false);
     666  std::vector<int> getMaskEdgeIndices(const std::vector<bool>& mask, bool getStartIndices=true);
     667  std::string formatBaselineParamsHeader(int whichrow, const std::string& masklist, bool verbose) const;
     668  std::string formatBaselineParamsFooter(float rms, bool verbose) const;
     669  std::vector<bool> getCompositeChanMask(int whichrow, const std::vector<bool>& inMask);
     670  //std::vector<bool> getCompositeChanMask(int whichrow, const std::vector<bool>& inMask, const std::vector<int>& edge, const int minEdgeSize, STLineFinder& lineFinder);
    611671
    612672  void applyChanFlag( casa::uInt whichrow, const std::vector<bool>& msk, casa::uChar flagval);

  • trunk/src/ScantableWrapper.h

    r1994 r2012  
    230230  { return table_->getFit(whichrow); }
    231231
    232   void calculateAZEL() { table_->calculateAZEL(); };
     232  void calculateAZEL() { table_->calculateAZEL(); }
    233233
    234234  std::vector<std::string> columnNames() const
    … …  
    257257  { table_->reshapeSpectrum( nmin, nmax ); }
    258258
    259   STFitEntry polyBaseline(const std::vector<bool>& mask, int order, int rowno)
    260   { return table_->polyBaseline(mask, order, rowno); }
    261 
    262   void polyBaselineBatch(const std::vector<bool>& mask, int order)
    263   { table_->polyBaselineBatch(mask, order); }
     259  void polyBaseline(const std::vector<bool>& mask, int order, bool outlog=false, const std::string& blfile="")
     260  { table_->polyBaseline(mask, order, outlog, blfile); }
     261
     262  void autoPolyBaseline(const std::vector<bool>& mask, int order, const std::vector<int>& edge, float threshold=5.0, int chan_avg_limit=1, bool outlog=false, const std::string& blfile="")
     263  { table_->autoPolyBaseline(mask, order, edge, threshold, chan_avg_limit, outlog, blfile); }
     264
     265  void cubicSplineBaseline(const std::vector<bool>& mask, int npiece, float clipthresh, int clipniter, bool outlog=false, const std::string& blfile="")
     266  { table_->cubicSplineBaseline(mask, npiece, clipthresh, clipniter, outlog, blfile); }
     267
     268  void autoCubicSplineBaseline(const std::vector<bool>& mask, int npiece, float clipthresh, int clipniter, const std::vector<int>& edge, float threshold=5.0, int chan_avg_limit=1, bool outlog=false, const std::string& blfile="")
     269  { table_->autoCubicSplineBaseline(mask, npiece, clipthresh, clipniter, edge, threshold, chan_avg_limit, outlog, blfile); }
     270
     271  float getRms(const std::vector<bool>& mask, int whichrow)
     272  { return table_->getRms(mask, whichrow); }
     273
     274  std::string formatBaselineParams(const std::vector<float>& params, const std::vector<bool>& fixed, float rms, const std::string& masklist, int whichrow, bool verbose=false)
     275  { return table_->formatBaselineParams(params, fixed, rms, masklist, whichrow, verbose); }
     276
     277  std::string formatPiecewiseBaselineParams(const std::vector<int>& ranges, const std::vector<float>& params, const std::vector<bool>& fixed, float rms, const std::string& masklist, int whichrow, bool verbose=false)
     278  { return table_->formatPiecewiseBaselineParams(ranges, params, fixed, rms, masklist, whichrow, verbose); }
    264279
    265280  bool getFlagtraFast(int whichrow=0) const
    266     { return table_->getFlagtraFast(whichrow); }
     281  { return table_->getFlagtraFast(whichrow); }
     282
    267283
    268284

  • trunk/src/python_STLineFinder.cpp

    r894 r2012  
    4242         .def("setoptions",&STLineFinder::setOptions)
    4343         .def("setscan",&STLineFinder::setScan)
     44         .def("setdata",&STLineFinder::setData)
    4445         .def("findlines",&STLineFinder::findLines)
    4546         .def("getmask",&STLineFinder::getMask)

  • trunk/src/python_Scantable.cpp

    r1947 r2012  
    143143          boost::python::arg("nmax")=-1) )
    144144    .def("_poly_baseline", &ScantableWrapper::polyBaseline)
    145     .def("_poly_baseline_batch", &ScantableWrapper::polyBaselineBatch)
     145    .def("_auto_poly_baseline", &ScantableWrapper::autoPolyBaseline)
     146    .def("_cspline_baseline", &ScantableWrapper::cubicSplineBaseline)
     147    .def("_auto_cspline_baseline", &ScantableWrapper::autoCubicSplineBaseline)
     148    .def("get_rms", &ScantableWrapper::getRms)
     149    .def("format_blparams_row", &ScantableWrapper::formatBaselineParams)
     150    .def("format_piecewise_blparams_row", &ScantableWrapper::formatPiecewiseBaselineParams)
    146151    .def("_getflagtrafast", &ScantableWrapper::getFlagtraFast,
    147152         (boost::python::arg("whichrow")=0) )
     153    //.def("_sspline_baseline", &ScantableWrapper::smoothingSplineBaseline,
     154    // (boost::python::arg("thres")=3.0,
     155    //  boost::python::arg("niter")=1) )
     156    //.def("_test_cin", &ScantableWrapper::testCin)
    148157  ;
    149158};
Note: See TracChangeset for help on using the changeset viewer.