Changeset 2193


Ignore:
Timestamp:
06/14/11 20:47:29 (13 years ago)
Author:
WataruKawasaki
Message:

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s):

Description: modified Scantable::*Baseline() so that the number of clipped channels be output.


Location:
trunk/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/Scantable.cpp

    r2189 r2193  
    18191819void Scantable::polyBaseline(const std::vector<bool>& mask, int order, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
    18201820{
    1821   ofstream ofs;
    1822   String coordInfo = "";
    1823   bool hasSameNchan = true;
    1824   bool outTextFile = false;
    1825 
    1826   if (blfile != "") {
    1827     ofs.open(blfile.c_str(), ios::out | ios::app);
    1828     if (ofs) outTextFile = true;
    1829   }
    1830 
    1831   if (outLogger || outTextFile) {
    1832     coordInfo = getCoordInfo()[0];
    1833     if (coordInfo == "") coordInfo = "channel";
    1834     hasSameNchan = hasSameNchanOverIFs();
    1835   }
    1836 
    1837   Fitter fitter = Fitter();
    1838   fitter.setExpression("poly", order);
    1839   //fitter.setIterClipping(thresClip, nIterClip);
    1840 
    1841   int nRow = nrow();
    1842   std::vector<bool> chanMask;
    1843   bool showProgress;
    1844   int minNRow;
    1845   parseProgressInfo(progressInfo, showProgress, minNRow);
    1846 
    1847   for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    1848     chanMask = getCompositeChanMask(whichrow, mask);
    1849     fitBaseline(chanMask, whichrow, fitter);
    1850     setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
    1851     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "polyBaseline()", fitter);
    1852     showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
    1853   }
    1854 
    1855   if (outTextFile) ofs.close();
     1821  try {
     1822    ofstream ofs;
     1823    String coordInfo = "";
     1824    bool hasSameNchan = true;
     1825    bool outTextFile = false;
     1826
     1827    if (blfile != "") {
     1828      ofs.open(blfile.c_str(), ios::out | ios::app);
     1829      if (ofs) outTextFile = true;
     1830    }
     1831
     1832    if (outLogger || outTextFile) {
     1833      coordInfo = getCoordInfo()[0];
     1834      if (coordInfo == "") coordInfo = "channel";
     1835      hasSameNchan = hasSameNchanOverIFs();
     1836    }
     1837
     1838    Fitter fitter = Fitter();
     1839    fitter.setExpression("poly", order);
     1840    //fitter.setIterClipping(thresClip, nIterClip);
     1841
     1842    int nRow = nrow();
     1843    std::vector<bool> chanMask;
     1844    bool showProgress;
     1845    int minNRow;
     1846    parseProgressInfo(progressInfo, showProgress, minNRow);
     1847
     1848    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
     1849      chanMask = getCompositeChanMask(whichrow, mask);
     1850      fitBaseline(chanMask, whichrow, fitter);
     1851      setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
     1852      outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "polyBaseline()", fitter);
     1853      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
     1854    }
     1855
     1856    if (outTextFile) ofs.close();
     1857
     1858  } catch (...) {
     1859    throw;
     1860  }
    18561861}
    18571862
    18581863void Scantable::autoPolyBaseline(const std::vector<bool>& mask, int order, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
    18591864{
    1860   ofstream ofs;
    1861   String coordInfo = "";
    1862   bool hasSameNchan = true;
    1863   bool outTextFile = false;
    1864 
    1865   if (blfile != "") {
    1866     ofs.open(blfile.c_str(), ios::out | ios::app);
    1867     if (ofs) outTextFile = true;
    1868   }
    1869 
    1870   if (outLogger || outTextFile) {
    1871     coordInfo = getCoordInfo()[0];
    1872     if (coordInfo == "") coordInfo = "channel";
    1873     hasSameNchan = hasSameNchanOverIFs();
    1874   }
    1875 
    1876   Fitter fitter = Fitter();
    1877   fitter.setExpression("poly", order);
    1878   //fitter.setIterClipping(thresClip, nIterClip);
    1879 
    1880   int nRow = nrow();
    1881   std::vector<bool> chanMask;
    1882   int minEdgeSize = getIFNos().size()*2;
    1883   STLineFinder lineFinder = STLineFinder();
    1884   lineFinder.setOptions(threshold, 3, chanAvgLimit);
    1885 
    1886   bool showProgress;
    1887   int minNRow;
    1888   parseProgressInfo(progressInfo, showProgress, minNRow);
    1889 
    1890   for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    1891 
    1892     //-------------------------------------------------------
    1893     //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
    1894     //-------------------------------------------------------
    1895     int edgeSize = edge.size();
    1896     std::vector<int> currentEdge;
    1897     if (edgeSize >= 2) {
    1898       int idx = 0;
    1899       if (edgeSize > 2) {
    1900         if (edgeSize < minEdgeSize) {
    1901           throw(AipsError("Length of edge element info is less than that of IFs"));
     1865  try {
     1866    ofstream ofs;
     1867    String coordInfo = "";
     1868    bool hasSameNchan = true;
     1869    bool outTextFile = false;
     1870
     1871    if (blfile != "") {
     1872      ofs.open(blfile.c_str(), ios::out | ios::app);
     1873      if (ofs) outTextFile = true;
     1874    }
     1875
     1876    if (outLogger || outTextFile) {
     1877      coordInfo = getCoordInfo()[0];
     1878      if (coordInfo == "") coordInfo = "channel";
     1879      hasSameNchan = hasSameNchanOverIFs();
     1880    }
     1881
     1882    Fitter fitter = Fitter();
     1883    fitter.setExpression("poly", order);
     1884    //fitter.setIterClipping(thresClip, nIterClip);
     1885
     1886    int nRow = nrow();
     1887    std::vector<bool> chanMask;
     1888    int minEdgeSize = getIFNos().size()*2;
     1889    STLineFinder lineFinder = STLineFinder();
     1890    lineFinder.setOptions(threshold, 3, chanAvgLimit);
     1891
     1892    bool showProgress;
     1893    int minNRow;
     1894    parseProgressInfo(progressInfo, showProgress, minNRow);
     1895
     1896    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
     1897
     1898      //-------------------------------------------------------
     1899      //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
     1900      //-------------------------------------------------------
     1901      int edgeSize = edge.size();
     1902      std::vector<int> currentEdge;
     1903      if (edgeSize >= 2) {
     1904        int idx = 0;
     1905        if (edgeSize > 2) {
     1906          if (edgeSize < minEdgeSize) {
     1907            throw(AipsError("Length of edge element info is less than that of IFs"));
     1908          }
     1909          idx = 2 * getIF(whichrow);
    19021910        }
    1903         idx = 2 * getIF(whichrow);
    1904       }
    1905       currentEdge.push_back(edge[idx]);
    1906       currentEdge.push_back(edge[idx+1]);
    1907     } else {
    1908       throw(AipsError("Wrong length of edge element"));
    1909     }
    1910     lineFinder.setData(getSpectrum(whichrow));
    1911     lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
    1912     chanMask = lineFinder.getMask();
    1913     //-------------------------------------------------------
    1914 
    1915     fitBaseline(chanMask, whichrow, fitter);
    1916     setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
    1917 
    1918     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoPolyBaseline()", fitter);
    1919     showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
    1920   }
    1921 
    1922   if (outTextFile) ofs.close();
     1911        currentEdge.push_back(edge[idx]);
     1912        currentEdge.push_back(edge[idx+1]);
     1913      } else {
     1914        throw(AipsError("Wrong length of edge element"));
     1915      }
     1916      lineFinder.setData(getSpectrum(whichrow));
     1917      lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
     1918      chanMask = lineFinder.getMask();
     1919      //-------------------------------------------------------
     1920
     1921      fitBaseline(chanMask, whichrow, fitter);
     1922      setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
     1923
     1924      outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoPolyBaseline()", fitter);
     1925      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
     1926    }
     1927
     1928    if (outTextFile) ofs.close();
     1929
     1930  } catch (...) {
     1931    throw;
     1932  }
    19231933}
    19241934
    19251935void Scantable::cubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
    19261936{
    1927   ofstream ofs;
    1928   String coordInfo = "";
    1929   bool hasSameNchan = true;
    1930   bool outTextFile = false;
    1931 
    1932   if (blfile != "") {
    1933     ofs.open(blfile.c_str(), ios::out | ios::app);
    1934     if (ofs) outTextFile = true;
    1935   }
    1936 
    1937   if (outLogger || outTextFile) {
    1938     coordInfo = getCoordInfo()[0];
    1939     if (coordInfo == "") coordInfo = "channel";
    1940     hasSameNchan = hasSameNchanOverIFs();
    1941   }
    1942 
    1943   //Fitter fitter = Fitter();
    1944   //fitter.setExpression("cspline", nPiece);
    1945   //fitter.setIterClipping(thresClip, nIterClip);
    1946 
    1947   int nRow = nrow();
    1948   std::vector<bool> chanMask;
    1949   bool showProgress;
    1950   int minNRow;
    1951   parseProgressInfo(progressInfo, showProgress, minNRow);
    1952 
    1953   for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    1954     chanMask = getCompositeChanMask(whichrow, mask);
    1955     //fitBaseline(chanMask, whichrow, fitter);
    1956     //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
    1957     std::vector<int> pieceEdges;
    1958     std::vector<float> params;
    1959     std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, nPiece, pieceEdges, params, thresClip, nIterClip, getResidual);
    1960     setSpectrum(res, whichrow);
    1961     //
    1962 
    1963     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()", pieceEdges, params);
    1964     showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
    1965   }
    1966 
    1967   if (outTextFile) ofs.close();
     1937  try {
     1938    ofstream ofs;
     1939    String coordInfo = "";
     1940    bool hasSameNchan = true;
     1941    bool outTextFile = false;
     1942
     1943    if (blfile != "") {
     1944      ofs.open(blfile.c_str(), ios::out | ios::app);
     1945      if (ofs) outTextFile = true;
     1946    }
     1947
     1948    if (outLogger || outTextFile) {
     1949      coordInfo = getCoordInfo()[0];
     1950      if (coordInfo == "") coordInfo = "channel";
     1951      hasSameNchan = hasSameNchanOverIFs();
     1952    }
     1953
     1954    //Fitter fitter = Fitter();
     1955    //fitter.setExpression("cspline", nPiece);
     1956    //fitter.setIterClipping(thresClip, nIterClip);
     1957
     1958    int nRow = nrow();
     1959    std::vector<bool> chanMask;
     1960    bool showProgress;
     1961    int minNRow;
     1962    parseProgressInfo(progressInfo, showProgress, minNRow);
     1963
     1964    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
     1965      chanMask = getCompositeChanMask(whichrow, mask);
     1966      //fitBaseline(chanMask, whichrow, fitter);
     1967      //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
     1968      std::vector<int> pieceEdges;
     1969      std::vector<float> params;
     1970      int nClipped = 0;
     1971      std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, nPiece, pieceEdges, params, nClipped, thresClip, nIterClip, getResidual);
     1972      setSpectrum(res, whichrow);
     1973      //
     1974
     1975      outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()", pieceEdges, params, nClipped);
     1976      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
     1977    }
     1978
     1979    if (outTextFile) ofs.close();
     1980
     1981  } catch (...) {
     1982    throw;
     1983  }
    19681984}
    19691985
    19701986void Scantable::autoCubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
    19711987{
    1972   ofstream ofs;
    1973   String coordInfo = "";
    1974   bool hasSameNchan = true;
    1975   bool outTextFile = false;
    1976 
    1977   if (blfile != "") {
    1978     ofs.open(blfile.c_str(), ios::out | ios::app);
    1979     if (ofs) outTextFile = true;
    1980   }
    1981 
    1982   if (outLogger || outTextFile) {
    1983     coordInfo = getCoordInfo()[0];
    1984     if (coordInfo == "") coordInfo = "channel";
    1985     hasSameNchan = hasSameNchanOverIFs();
    1986   }
    1987 
    1988   //Fitter fitter = Fitter();
    1989   //fitter.setExpression("cspline", nPiece);
    1990   //fitter.setIterClipping(thresClip, nIterClip);
    1991 
    1992   int nRow = nrow();
    1993   std::vector<bool> chanMask;
    1994   int minEdgeSize = getIFNos().size()*2;
    1995   STLineFinder lineFinder = STLineFinder();
    1996   lineFinder.setOptions(threshold, 3, chanAvgLimit);
    1997 
    1998   bool showProgress;
    1999   int minNRow;
    2000   parseProgressInfo(progressInfo, showProgress, minNRow);
    2001 
    2002   for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    2003 
    2004     //-------------------------------------------------------
    2005     //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
    2006     //-------------------------------------------------------
    2007     int edgeSize = edge.size();
    2008     std::vector<int> currentEdge;
    2009     if (edgeSize >= 2) {
    2010       int idx = 0;
    2011       if (edgeSize > 2) {
    2012         if (edgeSize < minEdgeSize) {
    2013           throw(AipsError("Length of edge element info is less than that of IFs"));
     1988  try {
     1989    ofstream ofs;
     1990    String coordInfo = "";
     1991    bool hasSameNchan = true;
     1992    bool outTextFile = false;
     1993
     1994    if (blfile != "") {
     1995      ofs.open(blfile.c_str(), ios::out | ios::app);
     1996      if (ofs) outTextFile = true;
     1997    }
     1998
     1999    if (outLogger || outTextFile) {
     2000      coordInfo = getCoordInfo()[0];
     2001      if (coordInfo == "") coordInfo = "channel";
     2002      hasSameNchan = hasSameNchanOverIFs();
     2003    }
     2004
     2005    //Fitter fitter = Fitter();
     2006    //fitter.setExpression("cspline", nPiece);
     2007    //fitter.setIterClipping(thresClip, nIterClip);
     2008
     2009    int nRow = nrow();
     2010    std::vector<bool> chanMask;
     2011    int minEdgeSize = getIFNos().size()*2;
     2012    STLineFinder lineFinder = STLineFinder();
     2013    lineFinder.setOptions(threshold, 3, chanAvgLimit);
     2014
     2015    bool showProgress;
     2016    int minNRow;
     2017    parseProgressInfo(progressInfo, showProgress, minNRow);
     2018
     2019    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
     2020
     2021      //-------------------------------------------------------
     2022      //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
     2023      //-------------------------------------------------------
     2024      int edgeSize = edge.size();
     2025      std::vector<int> currentEdge;
     2026      if (edgeSize >= 2) {
     2027        int idx = 0;
     2028        if (edgeSize > 2) {
     2029          if (edgeSize < minEdgeSize) {
     2030            throw(AipsError("Length of edge element info is less than that of IFs"));
     2031          }
     2032          idx = 2 * getIF(whichrow);
    20142033        }
    2015         idx = 2 * getIF(whichrow);
    2016       }
    2017       currentEdge.push_back(edge[idx]);
    2018       currentEdge.push_back(edge[idx+1]);
    2019     } else {
    2020       throw(AipsError("Wrong length of edge element"));
    2021     }
    2022     lineFinder.setData(getSpectrum(whichrow));
    2023     lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
    2024     chanMask = lineFinder.getMask();
    2025     //-------------------------------------------------------
    2026 
    2027 
    2028     //fitBaseline(chanMask, whichrow, fitter);
    2029     //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
    2030     std::vector<int> pieceEdges;
    2031     std::vector<float> params;
    2032     std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, nPiece, pieceEdges, params, thresClip, nIterClip, getResidual);
    2033     setSpectrum(res, whichrow);
    2034     //
    2035 
    2036     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()", pieceEdges, params);
    2037     showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
    2038   }
    2039 
    2040   if (outTextFile) ofs.close();
    2041 }
    2042 
    2043 std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data, const std::vector<bool>& mask, int nPiece, std::vector<int>& idxEdge, std::vector<float>& params, float thresClip, int nIterClip, bool getResidual)
     2034        currentEdge.push_back(edge[idx]);
     2035        currentEdge.push_back(edge[idx+1]);
     2036      } else {
     2037        throw(AipsError("Wrong length of edge element"));
     2038      }
     2039      lineFinder.setData(getSpectrum(whichrow));
     2040      lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
     2041      chanMask = lineFinder.getMask();
     2042      //-------------------------------------------------------
     2043
     2044
     2045      //fitBaseline(chanMask, whichrow, fitter);
     2046      //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
     2047      std::vector<int> pieceEdges;
     2048      std::vector<float> params;
     2049      int nClipped = 0;
     2050      std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, nPiece, pieceEdges, params, nClipped, thresClip, nIterClip, getResidual);
     2051      setSpectrum(res, whichrow);
     2052      //
     2053
     2054      outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()", pieceEdges, params, nClipped);
     2055      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
     2056    }
     2057
     2058    if (outTextFile) ofs.close();
     2059
     2060  } catch (...) {
     2061    throw;
     2062  }
     2063}
     2064
     2065std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data, const std::vector<bool>& mask, int nPiece, std::vector<int>& idxEdge, std::vector<float>& params, int& nClipped, float thresClip, int nIterClip, bool getResidual)
    20442066{
    20452067  if (data.size() != mask.size()) {
    … …  
    20612083
    20622084  int initNData = x.size();
     2085  if (initNData < nPiece) {
     2086    throw(AipsError("too few non-flagged channels"));
     2087  }
    20632088
    20642089  int nElement = (int)(floor(floor((double)(initNData/nPiece))+0.5));
    … …  
    21072132
    21082133    for (int n = 0; n < nPiece; ++n) {
     2134      int nUseDataInPiece = 0;
    21092135      for (int i = idxEdge[n]; i < idxEdge[n+1]; ++i) {
    21102136
    … …  
    21422168        }
    21432169
     2170        nUseDataInPiece++;
     2171      }
     2172
     2173      if (nUseDataInPiece < 1) {
     2174        std::vector<string> suffixOfPieceNumber(4);
     2175        suffixOfPieceNumber[0] = "th";
     2176        suffixOfPieceNumber[1] = "st";
     2177        suffixOfPieceNumber[2] = "nd";
     2178        suffixOfPieceNumber[3] = "rd";
     2179        int idxNoDataPiece = (n % 10 <= 3) ? n : 0;
     2180        ostringstream oss;
     2181        oss << "all channels clipped or masked in " << n << suffixOfPieceNumber[idxNoDataPiece];
     2182        oss << " piece of the spectrum. can't execute fitting anymore.";
     2183        throw(AipsError(String(oss)));
    21442184      }
    21452185    }
    … …  
    22462286    }
    22472287  }
     2288
     2289  nClipped = initNData - nData;
    22482290
    22492291  std::vector<float> result;
    … …  
    23862428void Scantable::sinusoidBaseline(const std::vector<bool>& mask, const bool applyFFT, const std::string& fftMethod, const std::string& fftThresh, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, float thresClip, int nIterClip, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
    23872429{
    2388   ofstream ofs;
    2389   String coordInfo = "";
    2390   bool hasSameNchan = true;
    2391   bool outTextFile = false;
    2392 
    2393   if (blfile != "") {
    2394     ofs.open(blfile.c_str(), ios::out | ios::app);
    2395     if (ofs) outTextFile = true;
    2396   }
    2397 
    2398   if (outLogger || outTextFile) {
    2399     coordInfo = getCoordInfo()[0];
    2400     if (coordInfo == "") coordInfo = "channel";
    2401     hasSameNchan = hasSameNchanOverIFs();
    2402   }
    2403 
    2404   //Fitter fitter = Fitter();
    2405   //fitter.setExpression("sinusoid", nWaves);
    2406   //fitter.setIterClipping(thresClip, nIterClip);
    2407 
    2408   int nRow = nrow();
    2409   std::vector<bool> chanMask;
    2410   std::vector<int> nWaves;
    2411 
    2412   bool showProgress;
    2413   int minNRow;
    2414   parseProgressInfo(progressInfo, showProgress, minNRow);
    2415 
    2416   for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    2417     chanMask = getCompositeChanMask(whichrow, mask);
    2418     selectWaveNumbers(whichrow, chanMask, applyFFT, fftMethod, fftThresh, addNWaves, rejectNWaves, nWaves);
    2419 
    2420     //FOR DEBUGGING------------
    2421     if (whichrow < 0) {// == nRow -1) {
    2422       cout << "+++ i=" << setw(3) << whichrow << ", IF=" << setw(2) << getIF(whichrow);
    2423       if (applyFFT) {
     2430  try {
     2431    ofstream ofs;
     2432    String coordInfo = "";
     2433    bool hasSameNchan = true;
     2434    bool outTextFile = false;
     2435
     2436    if (blfile != "") {
     2437      ofs.open(blfile.c_str(), ios::out | ios::app);
     2438      if (ofs) outTextFile = true;
     2439    }
     2440
     2441    if (outLogger || outTextFile) {
     2442      coordInfo = getCoordInfo()[0];
     2443      if (coordInfo == "") coordInfo = "channel";
     2444      hasSameNchan = hasSameNchanOverIFs();
     2445    }
     2446
     2447    //Fitter fitter = Fitter();
     2448    //fitter.setExpression("sinusoid", nWaves);
     2449    //fitter.setIterClipping(thresClip, nIterClip);
     2450
     2451    int nRow = nrow();
     2452    std::vector<bool> chanMask;
     2453    std::vector<int> nWaves;
     2454
     2455    bool showProgress;
     2456    int minNRow;
     2457    parseProgressInfo(progressInfo, showProgress, minNRow);
     2458
     2459    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
     2460      chanMask = getCompositeChanMask(whichrow, mask);
     2461      selectWaveNumbers(whichrow, chanMask, applyFFT, fftMethod, fftThresh, addNWaves, rejectNWaves, nWaves);
     2462
     2463      //FOR DEBUGGING------------
     2464      if (whichrow < 0) {// == nRow -1) {
     2465        cout << "+++ i=" << setw(3) << whichrow << ", IF=" << setw(2) << getIF(whichrow);
     2466        if (applyFFT) {
    24242467          cout << "[ ";
    24252468          for (uInt j = 0; j < nWaves.size(); ++j) {
    … …  
    24272470          }
    24282471          cout << " ]    " << endl;
    2429       }
    2430       cout << flush;
    2431     }
    2432     //-------------------------
    2433 
    2434     //fitBaseline(chanMask, whichrow, fitter);
    2435     //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
    2436     std::vector<float> params;
    2437     std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, params, thresClip, nIterClip, getResidual);
    2438     setSpectrum(res, whichrow);
    2439     //
    2440 
    2441     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "sinusoidBaseline()", params);
    2442     showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
    2443   }
    2444 
    2445   if (outTextFile) ofs.close();
     2472        }
     2473        cout << flush;
     2474      }
     2475      //-------------------------
     2476
     2477      //fitBaseline(chanMask, whichrow, fitter);
     2478      //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
     2479      std::vector<float> params;
     2480      int nClipped = 0;
     2481      std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, params, nClipped, thresClip, nIterClip, getResidual);
     2482      setSpectrum(res, whichrow);
     2483      //
     2484
     2485      outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "sinusoidBaseline()", params, nClipped);
     2486      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
     2487    }
     2488
     2489    if (outTextFile) ofs.close();
     2490
     2491  } catch (...) {
     2492    throw;
     2493  }
    24462494}
    24472495
    24482496void Scantable::autoSinusoidBaseline(const std::vector<bool>& mask, const bool applyFFT, const std::string& fftMethod, const std::string& fftThresh, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, float thresClip, int nIterClip, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
    24492497{
    2450   ofstream ofs;
    2451   String coordInfo = "";
    2452   bool hasSameNchan = true;
    2453   bool outTextFile = false;
    2454 
    2455   if (blfile != "") {
    2456     ofs.open(blfile.c_str(), ios::out | ios::app);
    2457     if (ofs) outTextFile = true;
    2458   }
    2459 
    2460   if (outLogger || outTextFile) {
    2461     coordInfo = getCoordInfo()[0];
    2462     if (coordInfo == "") coordInfo = "channel";
    2463     hasSameNchan = hasSameNchanOverIFs();
    2464   }
    2465 
    2466   //Fitter fitter = Fitter();
    2467   //fitter.setExpression("sinusoid", nWaves);
    2468   //fitter.setIterClipping(thresClip, nIterClip);
    2469 
    2470   int nRow = nrow();
    2471   std::vector<bool> chanMask;
    2472   std::vector<int> nWaves;
    2473 
    2474   int minEdgeSize = getIFNos().size()*2;
    2475   STLineFinder lineFinder = STLineFinder();
    2476   lineFinder.setOptions(threshold, 3, chanAvgLimit);
    2477 
    2478   bool showProgress;
    2479   int minNRow;
    2480   parseProgressInfo(progressInfo, showProgress, minNRow);
    2481 
    2482   for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    2483 
    2484     //-------------------------------------------------------
    2485     //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
    2486     //-------------------------------------------------------
    2487     int edgeSize = edge.size();
    2488     std::vector<int> currentEdge;
    2489     if (edgeSize >= 2) {
    2490       int idx = 0;
    2491       if (edgeSize > 2) {
    2492         if (edgeSize < minEdgeSize) {
    2493           throw(AipsError("Length of edge element info is less than that of IFs"));
     2498  try {
     2499    ofstream ofs;
     2500    String coordInfo = "";
     2501    bool hasSameNchan = true;
     2502    bool outTextFile = false;
     2503
     2504    if (blfile != "") {
     2505      ofs.open(blfile.c_str(), ios::out | ios::app);
     2506      if (ofs) outTextFile = true;
     2507    }
     2508
     2509    if (outLogger || outTextFile) {
     2510      coordInfo = getCoordInfo()[0];
     2511      if (coordInfo == "") coordInfo = "channel";
     2512      hasSameNchan = hasSameNchanOverIFs();
     2513    }
     2514
     2515    //Fitter fitter = Fitter();
     2516    //fitter.setExpression("sinusoid", nWaves);
     2517    //fitter.setIterClipping(thresClip, nIterClip);
     2518
     2519    int nRow = nrow();
     2520    std::vector<bool> chanMask;
     2521    std::vector<int> nWaves;
     2522
     2523    int minEdgeSize = getIFNos().size()*2;
     2524    STLineFinder lineFinder = STLineFinder();
     2525    lineFinder.setOptions(threshold, 3, chanAvgLimit);
     2526
     2527    bool showProgress;
     2528    int minNRow;
     2529    parseProgressInfo(progressInfo, showProgress, minNRow);
     2530
     2531    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
     2532
     2533      //-------------------------------------------------------
     2534      //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
     2535      //-------------------------------------------------------
     2536      int edgeSize = edge.size();
     2537      std::vector<int> currentEdge;
     2538      if (edgeSize >= 2) {
     2539        int idx = 0;
     2540        if (edgeSize > 2) {
     2541          if (edgeSize < minEdgeSize) {
     2542            throw(AipsError("Length of edge element info is less than that of IFs"));
     2543          }
     2544          idx = 2 * getIF(whichrow);
    24942545        }
    2495         idx = 2 * getIF(whichrow);
    2496       }
    2497       currentEdge.push_back(edge[idx]);
    2498       currentEdge.push_back(edge[idx+1]);
    2499     } else {
    2500       throw(AipsError("Wrong length of edge element"));
    2501     }
    2502     lineFinder.setData(getSpectrum(whichrow));
    2503     lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
    2504     chanMask = lineFinder.getMask();
    2505     //-------------------------------------------------------
    2506 
    2507     selectWaveNumbers(whichrow, chanMask, applyFFT, fftMethod, fftThresh, addNWaves, rejectNWaves, nWaves);
    2508 
    2509     //fitBaseline(chanMask, whichrow, fitter);
    2510     //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
    2511     std::vector<float> params;
    2512     std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, params, thresClip, nIterClip, getResidual);
    2513     setSpectrum(res, whichrow);
    2514     //
    2515 
    2516     outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoSinusoidBaseline()", params);
    2517     showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
    2518   }
    2519 
    2520   if (outTextFile) ofs.close();
    2521 }
    2522 
    2523 std::vector<float> Scantable::doSinusoidFitting(const std::vector<float>& data, const std::vector<bool>& mask, const std::vector<int>& waveNumbers, std::vector<float>& params, float thresClip, int nIterClip, bool getResidual)
     2546        currentEdge.push_back(edge[idx]);
     2547        currentEdge.push_back(edge[idx+1]);
     2548      } else {
     2549        throw(AipsError("Wrong length of edge element"));
     2550      }
     2551      lineFinder.setData(getSpectrum(whichrow));
     2552      lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
     2553      chanMask = lineFinder.getMask();
     2554      //-------------------------------------------------------
     2555
     2556      selectWaveNumbers(whichrow, chanMask, applyFFT, fftMethod, fftThresh, addNWaves, rejectNWaves, nWaves);
     2557
     2558      //fitBaseline(chanMask, whichrow, fitter);
     2559      //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
     2560      std::vector<float> params;
     2561      int nClipped = 0;
     2562      std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, params, nClipped, thresClip, nIterClip, getResidual);
     2563      setSpectrum(res, whichrow);
     2564      //
     2565
     2566      outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoSinusoidBaseline()", params, nClipped);
     2567      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
     2568    }
     2569
     2570    if (outTextFile) ofs.close();
     2571
     2572  } catch (...) {
     2573    throw;
     2574  }
     2575}
     2576
     2577std::vector<float> Scantable::doSinusoidFitting(const std::vector<float>& data, const std::vector<bool>& mask, const std::vector<int>& waveNumbers, std::vector<float>& params, int& nClipped, float thresClip, int nIterClip, bool getResidual)
    25242578{
    25252579  if (data.size() != mask.size()) {
    … …  
    26172671    }
    26182672
     2673    int nUseData = 0;
    26192674    for (int k = 0; k < nChan; ++k) {
    26202675      if (maskArray[k] == 0) continue;
    … …  
    26262681        zMatrix[i] += z1[k] * xArray[i][k];
    26272682      }
     2683
     2684      nUseData++;
     2685    }
     2686
     2687    if (nUseData < 1) {
     2688        throw(AipsError("all channels clipped or masked. can't execute fitting anymore."));     
    26282689    }
    26292690
    … …  
    27132774  }
    27142775
     2776  nClipped = initNData - nData;
     2777
    27152778  std::vector<float> result;
    27162779  if (getResidual) {
    … …  
    27912854    if (outLogger) {
    27922855      LogIO ols(LogOrigin("Scantable", funcName, WHERE));
    2793       ols << formatBaselineParams(params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
     2856      ols << formatBaselineParams(params, fixed, rms, -1, masklist, whichrow, false) << LogIO::POST ;
    27942857    }
    27952858    if (outTextFile) {
    2796       ofs << formatBaselineParams(params, fixed, rms, masklist, whichrow, true) << flush;
     2859      ofs << formatBaselineParams(params, fixed, rms, -1, masklist, whichrow, true) << flush;
    27972860    }
    27982861  }
    … …  
    28002863
    28012864/* for cspline. will be merged once cspline is available in fitter (2011/3/10 WK) */
    2802 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)
     2865void 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 int nClipped)
    28032866{
    28042867  if (outLogger || outTextFile) {
    … …  
    28102873    if (outLogger) {
    28112874      LogIO ols(LogOrigin("Scantable", funcName, WHERE));
    2812       ols << formatPiecewiseBaselineParams(edge, params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
     2875      ols << formatPiecewiseBaselineParams(edge, params, fixed, rms, nClipped, masklist, whichrow, false) << LogIO::POST ;
    28132876    }
    28142877    if (outTextFile) {
    2815       ofs << formatPiecewiseBaselineParams(edge, params, fixed, rms, masklist, whichrow, true) << flush;
     2878      ofs << formatPiecewiseBaselineParams(edge, params, fixed, rms, nClipped, masklist, whichrow, true) << flush;
    28162879    }
    28172880  }
    … …  
    28192882
    28202883/* for sinusoid. will be merged once sinusoid is available in fitter (2011/3/10 WK) */
    2821 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)
     2884void 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 int nClipped)
    28222885{
    28232886  if (outLogger || outTextFile) {
    … …  
    28292892    if (outLogger) {
    28302893      LogIO ols(LogOrigin("Scantable", funcName, WHERE));
    2831       ols << formatBaselineParams(params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
     2894      ols << formatBaselineParams(params, fixed, rms, nClipped, masklist, whichrow, false) << LogIO::POST ;
    28322895    }
    28332896    if (outTextFile) {
    2834       ofs << formatBaselineParams(params, fixed, rms, masklist, whichrow, true) << flush;
     2897      ofs << formatBaselineParams(params, fixed, rms, nClipped, masklist, whichrow, true) << flush;
    28352898    }
    28362899  }
    … …  
    28542917    if (nInterval == 0) nInterval++;
    28552918
    2856     if (nProcessed == 0) {
    2857       printf("\x1b[31m\x1b[1m");             //set red color, highlighted
    2858       printf("[  0%%]");
    2859       printf("\x1b[39m\x1b[0m");             //set default attributes
    2860       fflush(NULL);
    2861     } else if (nProcessed % nInterval == 0) {
     2919    if (nProcessed % nInterval == 0) {
    28622920      printf("\r");                          //go to the head of line
    28632921      printf("\x1b[31m\x1b[1m");             //set red color, highlighted
    … …  
    28662924      fflush(NULL);
    28672925    }
     2926
    28682927    if (nProcessed == nTotal - 1) {
    28692928      printf("\r\x1b[K");                    //clear
    28702929      fflush(NULL);
    28712930    }
     2931
    28722932  }
    28732933}
    … …  
    29543014}
    29553015
    2956 std::string Scantable::formatBaselineParamsFooter(float rms, bool verbose) const
     3016std::string Scantable::formatBaselineParamsFooter(float rms, int nClipped, bool verbose) const
    29573017{
    29583018  ostringstream oss;
    … …  
    29613021    oss << "Results of baseline fit" << endl;
    29623022    oss << "  rms = " << setprecision(6) << rms << endl;
     3023    if (nClipped >= 0) {
     3024      oss << "  Number of clipped channels = " << nClipped << endl;
     3025    }
    29633026    for (int i = 0; i < 60; ++i) {
    29643027      oss << "-";
    … …  
    29743037                                            const std::vector<bool>& fixed,
    29753038                                            float rms,
     3039                                            int nClipped,
    29763040                                            const std::string& masklist,
    29773041                                            int whichrow,
    … …  
    30043068
    30053069    oss << endl;
    3006     oss << formatBaselineParamsFooter(rms, verbose);
     3070    oss << formatBaselineParamsFooter(rms, nClipped, verbose);
    30073071
    30083072    return String(oss);
    … …  
    30113075}
    30123076
    3013 std::string Scantable::formatPiecewiseBaselineParams(const std::vector<int>& ranges, const std::vector<float>& params, const std::vector<bool>& fixed, float rms, const std::string& masklist, int whichrow, bool verbose) const
     3077  std::string Scantable::formatPiecewiseBaselineParams(const std::vector<int>& ranges, const std::vector<float>& params, const std::vector<bool>& fixed, float rms, int nClipped, const std::string& masklist, int whichrow, bool verbose) const
    30143078{
    30153079  int nOutParam = (int)(params.size());
    … …  
    30193083    return("  Not fitted");
    30203084  } else if (nPiece < 0) {
    3021     return formatBaselineParams(params, fixed, rms, masklist, whichrow, verbose);
     3085    return formatBaselineParams(params, fixed, rms, nClipped, masklist, whichrow, verbose);
    30223086  } else if (nPiece < 1) {
    30233087    return("  Bad count of the piece edge info");
    … …  
    30393103      ss << "  [" << ranges[i] << "," << (ranges[i+1]-1) << "]";
    30403104      oss << left << setw(wRange) << ss.str();
    3041       oss << formatBaselineParams(params, fixed, rms, masklist, whichrow, false, i*nParam, nParam, true);
    3042     }
    3043 
    3044     oss << formatBaselineParamsFooter(rms, verbose);
     3105      oss << formatBaselineParams(params, fixed, rms, 0, masklist, whichrow, false, i*nParam, nParam, true);
     3106    }
     3107
     3108    oss << formatBaselineParamsFooter(rms, nClipped, verbose);
    30453109
    30463110    return String(oss);

  • trunk/src/Scantable.h

    r2189 r2193  
    568568                                   const std::vector<bool>& fixed,
    569569                                   float rms,
     570                                   int nClipped,
    570571                                   const std::string& masklist,
    571572                                   int whichrow,
    … …  
    578579                                            const std::vector<bool>& fixed,
    579580                                            float rms,
     581                                            int nClipped,
    580582                                            const std::string& masklist,
    581583                                            int whichrow,
    … …  
    701703                                          std::vector<int>& idxEdge,
    702704                                          std::vector<float>& params,
     705                                          int& nClipped,
    703706                                          float thresClip=3.0,
    704707                                          int nIterClip=1,
    … …  
    708711                                       const std::vector<int>& waveNumbers,
    709712                                       std::vector<float>& params,
     713                                       int& nClipped,
    710714                                       float thresClip=3.0,
    711715                                       int nIterClip=1,
    … …  
    741745  std::vector<int> getMaskEdgeIndices(const std::vector<bool>& mask);
    742746  std::string formatBaselineParamsHeader(int whichrow, const std::string& masklist, bool verbose) const;
    743   std::string formatBaselineParamsFooter(float rms, bool verbose) const;
     747  std::string formatBaselineParamsFooter(float rms, int nClipped, bool verbose) const;
    744748  std::vector<bool> getCompositeChanMask(int whichrow, const std::vector<bool>& inMask);
    745749  //std::vector<bool> getCompositeChanMask(int whichrow, const std::vector<bool>& inMask, const std::vector<int>& edge, const int minEdgeSize, STLineFinder& lineFinder);
    746750  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);
    747   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);
    748   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);
     751  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 int nClipped);
     752  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 int nClipped);
    749753  void parseProgressInfo(const std::string& progressInfo, bool& showProgress, int& minNRow);
    750754  void showProgressOnTerminal(const int nProcessed, const int nTotal, const bool showProgress=true, const int nTotalThreshold=1000);

  • trunk/src/ScantableWrapper.h

    r2189 r2193  
    286286
    287287  std::string formatBaselineParams(const std::vector<float>& params, const std::vector<bool>& fixed, float rms, const std::string& masklist, int whichrow, bool verbose=false)
    288   { return table_->formatBaselineParams(params, fixed, rms, masklist, whichrow, verbose); }
     288  { return table_->formatBaselineParams(params, fixed, rms, -1, masklist, whichrow, verbose); }
    289289
    290290  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)
    291   { return table_->formatPiecewiseBaselineParams(ranges, params, fixed, rms, masklist, whichrow, verbose); }
     291  { return table_->formatPiecewiseBaselineParams(ranges, params, fixed, rms, -1, masklist, whichrow, verbose); }
    292292
    293293  bool getFlagtraFast(int whichrow=0) const
Note: See TracChangeset for help on using the changeset viewer.