Changeset 2193
- Timestamp:
- 06/14/11 20:47:29 (14 years ago)
- Location:
- trunk/src
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/Scantable.cpp
r2189 r2193 1819 1819 void Scantable::polyBaseline(const std::vector<bool>& mask, int order, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile) 1820 1820 { 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 } 1856 1861 } 1857 1862 1858 1863 void 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) 1859 1864 { 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); 1902 1910 } 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 } 1923 1933 } 1924 1934 1925 1935 void 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) 1926 1936 { 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 } 1968 1984 } 1969 1985 1970 1986 void 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) 1971 1987 { 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); 2014 2033 } 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 2065 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, int& nClipped, float thresClip, int nIterClip, bool getResidual) 2044 2066 { 2045 2067 if (data.size() != mask.size()) { … … 2061 2083 2062 2084 int initNData = x.size(); 2085 if (initNData < nPiece) { 2086 throw(AipsError("too few non-flagged channels")); 2087 } 2063 2088 2064 2089 int nElement = (int)(floor(floor((double)(initNData/nPiece))+0.5)); … … 2107 2132 2108 2133 for (int n = 0; n < nPiece; ++n) { 2134 int nUseDataInPiece = 0; 2109 2135 for (int i = idxEdge[n]; i < idxEdge[n+1]; ++i) { 2110 2136 … … 2142 2168 } 2143 2169 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))); 2144 2184 } 2145 2185 } … … 2246 2286 } 2247 2287 } 2288 2289 nClipped = initNData - nData; 2248 2290 2249 2291 std::vector<float> result; … … 2386 2428 void 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) 2387 2429 { 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) { 2424 2467 cout << "[ "; 2425 2468 for (uInt j = 0; j < nWaves.size(); ++j) { … … 2427 2470 } 2428 2471 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 } 2446 2494 } 2447 2495 2448 2496 void 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) 2449 2497 { 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); 2494 2545 } 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 2577 std::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) 2524 2578 { 2525 2579 if (data.size() != mask.size()) { … … 2617 2671 } 2618 2672 2673 int nUseData = 0; 2619 2674 for (int k = 0; k < nChan; ++k) { 2620 2675 if (maskArray[k] == 0) continue; … … 2626 2681 zMatrix[i] += z1[k] * xArray[i][k]; 2627 2682 } 2683 2684 nUseData++; 2685 } 2686 2687 if (nUseData < 1) { 2688 throw(AipsError("all channels clipped or masked. can't execute fitting anymore.")); 2628 2689 } 2629 2690 … … 2713 2774 } 2714 2775 2776 nClipped = initNData - nData; 2777 2715 2778 std::vector<float> result; 2716 2779 if (getResidual) { … … 2791 2854 if (outLogger) { 2792 2855 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 ; 2794 2857 } 2795 2858 if (outTextFile) { 2796 ofs << formatBaselineParams(params, fixed, rms, masklist, whichrow, true) << flush;2859 ofs << formatBaselineParams(params, fixed, rms, -1, masklist, whichrow, true) << flush; 2797 2860 } 2798 2861 } … … 2800 2863 2801 2864 /* 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 )2865 void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, const std::vector<int>& edge, const std::vector<float>& params, const int nClipped) 2803 2866 { 2804 2867 if (outLogger || outTextFile) { … … 2810 2873 if (outLogger) { 2811 2874 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 ; 2813 2876 } 2814 2877 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; 2816 2879 } 2817 2880 } … … 2819 2882 2820 2883 /* 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 )2884 void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, const std::vector<float>& params, const int nClipped) 2822 2885 { 2823 2886 if (outLogger || outTextFile) { … … 2829 2892 if (outLogger) { 2830 2893 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 ; 2832 2895 } 2833 2896 if (outTextFile) { 2834 ofs << formatBaselineParams(params, fixed, rms, masklist, whichrow, true) << flush;2897 ofs << formatBaselineParams(params, fixed, rms, nClipped, masklist, whichrow, true) << flush; 2835 2898 } 2836 2899 } … … 2854 2917 if (nInterval == 0) nInterval++; 2855 2918 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) { 2862 2920 printf("\r"); //go to the head of line 2863 2921 printf("\x1b[31m\x1b[1m"); //set red color, highlighted … … 2866 2924 fflush(NULL); 2867 2925 } 2926 2868 2927 if (nProcessed == nTotal - 1) { 2869 2928 printf("\r\x1b[K"); //clear 2870 2929 fflush(NULL); 2871 2930 } 2931 2872 2932 } 2873 2933 } … … 2954 3014 } 2955 3015 2956 std::string Scantable::formatBaselineParamsFooter(float rms, bool verbose) const3016 std::string Scantable::formatBaselineParamsFooter(float rms, int nClipped, bool verbose) const 2957 3017 { 2958 3018 ostringstream oss; … … 2961 3021 oss << "Results of baseline fit" << endl; 2962 3022 oss << " rms = " << setprecision(6) << rms << endl; 3023 if (nClipped >= 0) { 3024 oss << " Number of clipped channels = " << nClipped << endl; 3025 } 2963 3026 for (int i = 0; i < 60; ++i) { 2964 3027 oss << "-"; … … 2974 3037 const std::vector<bool>& fixed, 2975 3038 float rms, 3039 int nClipped, 2976 3040 const std::string& masklist, 2977 3041 int whichrow, … … 3004 3068 3005 3069 oss << endl; 3006 oss << formatBaselineParamsFooter(rms, verbose);3070 oss << formatBaselineParamsFooter(rms, nClipped, verbose); 3007 3071 3008 3072 return String(oss); … … 3011 3075 } 3012 3076 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) const3077 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 3014 3078 { 3015 3079 int nOutParam = (int)(params.size()); … … 3019 3083 return(" Not fitted"); 3020 3084 } else if (nPiece < 0) { 3021 return formatBaselineParams(params, fixed, rms, masklist, whichrow, verbose);3085 return formatBaselineParams(params, fixed, rms, nClipped, masklist, whichrow, verbose); 3022 3086 } else if (nPiece < 1) { 3023 3087 return(" Bad count of the piece edge info"); … … 3039 3103 ss << " [" << ranges[i] << "," << (ranges[i+1]-1) << "]"; 3040 3104 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); 3045 3109 3046 3110 return String(oss); -
trunk/src/Scantable.h
r2189 r2193 568 568 const std::vector<bool>& fixed, 569 569 float rms, 570 int nClipped, 570 571 const std::string& masklist, 571 572 int whichrow, … … 578 579 const std::vector<bool>& fixed, 579 580 float rms, 581 int nClipped, 580 582 const std::string& masklist, 581 583 int whichrow, … … 701 703 std::vector<int>& idxEdge, 702 704 std::vector<float>& params, 705 int& nClipped, 703 706 float thresClip=3.0, 704 707 int nIterClip=1, … … 708 711 const std::vector<int>& waveNumbers, 709 712 std::vector<float>& params, 713 int& nClipped, 710 714 float thresClip=3.0, 711 715 int nIterClip=1, … … 741 745 std::vector<int> getMaskEdgeIndices(const std::vector<bool>& mask); 742 746 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; 744 748 std::vector<bool> getCompositeChanMask(int whichrow, const std::vector<bool>& inMask); 745 749 //std::vector<bool> getCompositeChanMask(int whichrow, const std::vector<bool>& inMask, const std::vector<int>& edge, const int minEdgeSize, STLineFinder& lineFinder); 746 750 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); 749 753 void parseProgressInfo(const std::string& progressInfo, bool& showProgress, int& minNRow); 750 754 void showProgressOnTerminal(const int nProcessed, const int nTotal, const bool showProgress=true, const int nTotalThreshold=1000); -
trunk/src/ScantableWrapper.h
r2189 r2193 286 286 287 287 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); } 289 289 290 290 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); } 292 292 293 293 bool getFlagtraFast(int whichrow=0) const
Note:
See TracChangeset
for help on using the changeset viewer.