Ignore:
Timestamp:
03/15/11 18:31:04 (13 years ago)
Author:
WataruKawasaki
Message:

New Development: Yes

JIRA Issue: Yes CAS-2847

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: Yes

Module(s): scantable

Description: added {auto_}sinusoid_baseline() for sinusoidal baseline fitting. also minor bug fixes for asapfitter.


File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/Scantable.cpp

    r2036 r2047  
    882882    uInt row = uInt(whichrow);
    883883    stpol->setSpectra(getPolMatrix(row));
    884     Float fang,fhand,parang;
     884    Float fang,fhand;
    885885    fang = focusTable_.getTotalAngle(mfocusidCol_(row));
    886886    fhand = focusTable_.getFeedHand(mfocusidCol_(row));
    … …  
    17491749}
    17501750
    1751 bool Scantable::getFlagtraFast(int whichrow)
     1751bool Scantable::getFlagtraFast(uInt whichrow)
    17521752{
    17531753  uChar flag;
    17541754  Vector<uChar> flags;
    1755   flagsCol_.get(uInt(whichrow), flags);
     1755  flagsCol_.get(whichrow, flags);
    17561756  flag = flags[0];
    1757   for (int i = 1; i < flags.size(); ++i) {
     1757  for (uInt i = 1; i < flags.size(); ++i) {
    17581758    flag &= flags[i];
    17591759  }
    … …  
    17611761}
    17621762
    1763 void Scantable::cubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, bool outLogger, const std::string& blfile) {
     1763void Scantable::polyBaseline(const std::vector<bool>& mask, int order, bool outLogger, const std::string& blfile)
     1764{
    17641765  ofstream ofs;
    17651766  String coordInfo;
    1766   bool hasSameNchan;
    1767   int firstIF;
     1767  bool hasSameNchan = true;
    17681768  bool outTextFile = false;
    17691769
    … …  
    17771777    if (coordInfo == "") coordInfo = "channel";
    17781778    hasSameNchan = hasSameNchanOverIFs();
    1779     firstIF = getIF(0);
    1780   }
    1781 
    1782   //Fitter fitter = Fitter();
    1783   //fitter.setExpression("cspline", nPiece, thresClip, nIterClip);
     1779  }
     1780
     1781  Fitter fitter = Fitter();
     1782  fitter.setExpression("poly", order);
    17841783
    17851784  int nRow = nrow();
    … …  
    17871786
    17881787  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    1789 
    17901788    chanMask = getCompositeChanMask(whichrow, mask);
    1791     //fitBaseline(chanMask, whichrow, fitter);
    1792     //setSpectrum(fitter.getResidual(), whichrow);
    1793     std::vector<int> pieceRanges;
    1794     std::vector<float> params;
    1795     std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, pieceRanges, params, nPiece, thresClip, nIterClip);
    1796     setSpectrum(res, whichrow);
    1797 
    1798     if (outLogger || outTextFile) {
    1799       std::vector<bool>  fixed;
    1800       fixed.clear();
    1801       float rms = getRms(chanMask, whichrow);
    1802       String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan, firstIF, true);
    1803 
    1804       if (outLogger) {
    1805         LogIO ols(LogOrigin("Scantable", "cubicSplineBaseline()", WHERE));
    1806         ols << formatPiecewiseBaselineParams(pieceRanges, params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
    1807       }
    1808       if (outTextFile) {
    1809         ofs << formatPiecewiseBaselineParams(pieceRanges, params, fixed, rms, masklist, whichrow, true) << flush;
    1810       }
    1811     }
    1812 
     1789    fitBaseline(chanMask, whichrow, fitter);
     1790    setSpectrum(fitter.getResidual(), whichrow);
     1791    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "polyBaseline()", fitter);
    18131792  }
    18141793
    … …  
    18161795}
    18171796
    1818 
    1819 void 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)
     1797void Scantable::autoPolyBaseline(const std::vector<bool>& mask, int order, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool outLogger, const std::string& blfile)
    18201798{
    18211799  ofstream ofs;
    18221800  String coordInfo;
    1823   bool hasSameNchan;
    1824   int firstIF;
     1801  bool hasSameNchan = true;
    18251802  bool outTextFile = false;
    18261803
    … …  
    18341811    if (coordInfo == "") coordInfo = "channel";
    18351812    hasSameNchan = hasSameNchanOverIFs();
    1836     firstIF = getIF(0);
    1837   }
    1838 
    1839   //Fitter fitter = Fitter();
    1840   //fitter.setExpression("cspline", nPiece, thresClip, nIterClip);
     1813  }
     1814
     1815  Fitter fitter = Fitter();
     1816  fitter.setExpression("poly", order);
    18411817
    18421818  int nRow = nrow();
    … …  
    18711847    //-------------------------------------------------------
    18721848
    1873 
    1874     //fitBaseline(chanMask, whichrow, fitter);
     1849    fitBaseline(chanMask, whichrow, fitter);
     1850    setSpectrum(fitter.getResidual(), whichrow);
     1851
     1852    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoPolyBaseline()", fitter);
     1853  }
     1854
     1855  if (outTextFile) ofs.close();
     1856}
     1857
     1858void Scantable::cubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, bool outLogger, const std::string& blfile) {
     1859  ofstream ofs;
     1860  String coordInfo;
     1861  bool hasSameNchan = true;
     1862  bool outTextFile = false;
     1863
     1864  if (blfile != "") {
     1865    ofs.open(blfile.c_str(), ios::out | ios::app);
     1866    if (ofs) outTextFile = true;
     1867  }
     1868
     1869  if (outLogger || outTextFile) {
     1870    coordInfo = getCoordInfo()[0];
     1871    if (coordInfo == "") coordInfo = "channel";
     1872    hasSameNchan = hasSameNchanOverIFs();
     1873  }
     1874
     1875  //Fitter fitter = Fitter();
     1876  //fitter.setExpression("cspline", nPiece);
     1877
     1878  int nRow = nrow();
     1879  std::vector<bool> chanMask;
     1880
     1881  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
     1882    chanMask = getCompositeChanMask(whichrow, mask);
     1883    //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
    18751884    //setSpectrum(fitter.getResidual(), whichrow);
    18761885    std::vector<int> pieceRanges;
    … …  
    18781887    std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, pieceRanges, params, nPiece, thresClip, nIterClip);
    18791888    setSpectrum(res, whichrow);
    1880 
    1881     if (outLogger || outTextFile) {
    1882       std::vector<bool> fixed;
    1883       fixed.clear();
    1884       float rms = getRms(chanMask, whichrow);
    1885       String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan, firstIF, true);
    1886 
    1887       if (outLogger) {
    1888         LogIO ols(LogOrigin("Scantable", "autoCubicSplineBaseline()", WHERE));
    1889         ols << formatPiecewiseBaselineParams(pieceRanges, params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
    1890       }
    1891       if (outTextFile) {
    1892         ofs << formatPiecewiseBaselineParams(pieceRanges, params, fixed, rms, masklist, whichrow, true) << flush;
    1893       }
    1894     }
    1895 
     1889    //
     1890
     1891    std::vector<bool>  fixed;
     1892    fixed.clear();
     1893    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()", pieceRanges, params, fixed);
    18961894  }
    18971895
    … …  
    18991897}
    19001898
     1899void 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)
     1900{
     1901  ofstream ofs;
     1902  String coordInfo;
     1903  bool hasSameNchan = true;
     1904  bool outTextFile = false;
     1905
     1906  if (blfile != "") {
     1907    ofs.open(blfile.c_str(), ios::out | ios::app);
     1908    if (ofs) outTextFile = true;
     1909  }
     1910
     1911  if (outLogger || outTextFile) {
     1912    coordInfo = getCoordInfo()[0];
     1913    if (coordInfo == "") coordInfo = "channel";
     1914    hasSameNchan = hasSameNchanOverIFs();
     1915  }
     1916
     1917  //Fitter fitter = Fitter();
     1918  //fitter.setExpression("cspline", nPiece);
     1919
     1920  int nRow = nrow();
     1921  std::vector<bool> chanMask;
     1922  int minEdgeSize = getIFNos().size()*2;
     1923  STLineFinder lineFinder = STLineFinder();
     1924  lineFinder.setOptions(threshold, 3, chanAvgLimit);
     1925
     1926  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
     1927
     1928    //-------------------------------------------------------
     1929    //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
     1930    //-------------------------------------------------------
     1931    int edgeSize = edge.size();
     1932    std::vector<int> currentEdge;
     1933    if (edgeSize >= 2) {
     1934      int idx = 0;
     1935      if (edgeSize > 2) {
     1936        if (edgeSize < minEdgeSize) {
     1937          throw(AipsError("Length of edge element info is less than that of IFs"));
     1938        }
     1939        idx = 2 * getIF(whichrow);
     1940      }
     1941      currentEdge.push_back(edge[idx]);
     1942      currentEdge.push_back(edge[idx+1]);
     1943    } else {
     1944      throw(AipsError("Wrong length of edge element"));
     1945    }
     1946    lineFinder.setData(getSpectrum(whichrow));
     1947    lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
     1948    chanMask = lineFinder.getMask();
     1949    //-------------------------------------------------------
     1950
     1951
     1952    //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
     1953    //setSpectrum(fitter.getResidual(), whichrow);
     1954    std::vector<int> pieceRanges;
     1955    std::vector<float> params;
     1956    std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, pieceRanges, params, nPiece, thresClip, nIterClip);
     1957    setSpectrum(res, whichrow);
     1958    //
     1959
     1960    std::vector<bool> fixed;
     1961    fixed.clear();
     1962    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()", pieceRanges, params, fixed);
     1963  }
     1964
     1965  if (outTextFile) ofs.close();
     1966}
    19011967
    19021968std::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) {
    … …  
    19211987  int nElement = (int)(floor(floor((double)(nData/nPiece))+0.5));
    19221988
    1923   std::vector<double> sectionList0, sectionList1, sectionListR;
    1924   sectionList0.push_back((double)x[0]);
     1989  std::vector<int> sectionList0, sectionList1;
     1990  std::vector<double> sectionListR;
     1991  sectionList0.push_back(x[0]);
    19251992  sectionRanges.clear();
    19261993  sectionRanges.push_back(x[0]);
    19271994  for (int i = 1; i < nPiece; ++i) {
    1928     double valX = (double)x[nElement*i];
     1995    int valX = x[nElement*i];
    19291996    sectionList0.push_back(valX);
    19301997    sectionList1.push_back(valX);
    1931     sectionListR.push_back(1.0/valX);
    1932 
    1933     sectionRanges.push_back(x[nElement*i]-1);
    1934     sectionRanges.push_back(x[nElement*i]);
    1935   }
    1936   sectionList1.push_back((double)(x[x.size()-1]+1));
     1998    sectionListR.push_back(1.0/(double)valX);
     1999
     2000    sectionRanges.push_back(valX-1);
     2001    sectionRanges.push_back(valX);
     2002  }
     2003  sectionList1.push_back(x[x.size()-1]+1);
    19372004  sectionRanges.push_back(x[x.size()-1]);
    19382005
    … …  
    21112178}
    21122179
    2113 void Scantable::fitBaseline(const std::vector<bool>& mask, int whichrow, Fitter& fitter)
    2114 {
    2115   std::vector<double> abcsd = getAbcissa(whichrow);
    2116   std::vector<float> abcs;
    2117   for (int i = 0; i < abcsd.size(); ++i) {
    2118     abcs.push_back((float)abcsd[i]);
    2119   }
    2120   std::vector<float> spec = getSpectrum(whichrow);
    2121 
    2122   fitter.setData(abcs, spec, mask);
    2123   fitter.lfit();
    2124 }
    2125 
    2126 std::vector<bool> Scantable::getCompositeChanMask(int whichrow, const std::vector<bool>& inMask)
    2127 {
    2128   std::vector<bool> chanMask = getMask(whichrow);
    2129   int chanMaskSize = chanMask.size();
    2130   if (chanMaskSize != inMask.size()) {
    2131     throw(AipsError("different mask sizes"));
    2132   }
    2133   for (int i = 0; i < chanMaskSize; ++i) {
    2134     chanMask[i] = chanMask[i] && inMask[i];
    2135   }
    2136 
    2137   return chanMask;
    2138 }
    2139 
    2140 /*
    2141 std::vector<bool> Scantable::getCompositeChanMask(int whichrow, const std::vector<bool>& inMask, const std::vector<int>& edge, const int minEdgeSize, STLineFinder& lineFinder)
    2142 {
    2143   int edgeSize = edge.size();
    2144   std::vector<int> currentEdge;
    2145   if (edgeSize >= 2) {
    2146       int idx = 0;
    2147       if (edgeSize > 2) {
    2148         if (edgeSize < minEdgeSize) {
    2149           throw(AipsError("Length of edge element info is less than that of IFs"));
    2150         }
    2151         idx = 2 * getIF(whichrow);
    2152       }
    2153       currentEdge.push_back(edge[idx]);
    2154       currentEdge.push_back(edge[idx+1]);
    2155   } else {
    2156     throw(AipsError("Wrong length of edge element"));
    2157   }
    2158 
    2159   lineFinder.setData(getSpectrum(whichrow));
    2160   lineFinder.findLines(getCompositeChanMask(whichrow, inMask), currentEdge, whichrow);
    2161 
    2162   return lineFinder.getMask();
    2163 }
    2164 */
    2165 
    2166 void Scantable::polyBaseline(const std::vector<bool>& mask, int order, bool outLogger, const std::string& blfile)
    2167 {
     2180  void Scantable::sinusoidBaseline(const std::vector<bool>& mask, int nMinWavesInSW, int nMaxWavesInSW, float thresClip, int nIterClip, bool outLogger, const std::string& blfile) {
    21682181  ofstream ofs;
    21692182  String coordInfo;
    2170   bool hasSameNchan;
    2171   int firstIF;
     2183  bool hasSameNchan = true;
    21722184  bool outTextFile = false;
    21732185
    … …  
    21812193    if (coordInfo == "") coordInfo = "channel";
    21822194    hasSameNchan = hasSameNchanOverIFs();
    2183     firstIF = getIF(0);
    2184   }
    2185 
    2186   Fitter fitter = Fitter();
    2187   fitter.setExpression("poly", order);
     2195  }
     2196
     2197  //Fitter fitter = Fitter();
     2198  //fitter.setExpression("sinusoid", nMaxWavesInSW);
    21882199
    21892200  int nRow = nrow();
    … …  
    21912202
    21922203  for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    2193 
    21942204    chanMask = getCompositeChanMask(whichrow, mask);
    2195     fitBaseline(chanMask, whichrow, fitter);
    2196     setSpectrum(fitter.getResidual(), whichrow);
    2197    
    2198     if (outLogger || outTextFile) {
    2199       std::vector<float> params = fitter.getParameters();
    2200       std::vector<bool>  fixed  = fitter.getFixedParameters();
    2201       float rms = getRms(chanMask, whichrow);
    2202       String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan, firstIF, true);
    2203 
    2204       if (outLogger) {
    2205         LogIO ols(LogOrigin("Scantable", "polyBaseline()", WHERE));
    2206         ols << formatBaselineParams(params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
    2207       }
    2208       if (outTextFile) {
    2209         ofs << formatBaselineParams(params, fixed, rms, masklist, whichrow, true) << flush;
    2210       }
    2211     }
    2212 
     2205    //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
     2206    //setSpectrum(fitter.getResidual(), whichrow);
     2207    std::vector<int> pieceRanges;
     2208    std::vector<float> params;
     2209    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nMinWavesInSW, nMaxWavesInSW, params, thresClip, nIterClip);
     2210    setSpectrum(res, whichrow);
     2211    //
     2212
     2213    std::vector<bool>  fixed;
     2214    fixed.clear();
     2215    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "sinusoidBaseline()", pieceRanges, params, fixed);
    22132216  }
    22142217
    … …  
    22162219}
    22172220
    2218 void Scantable::autoPolyBaseline(const std::vector<bool>& mask, int order, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool outLogger, const std::string& blfile)
     2221void Scantable::autoSinusoidBaseline(const std::vector<bool>& mask, int nMinWavesInSW, int nMaxWavesInSW, float thresClip, int nIterClip, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool outLogger, const std::string& blfile)
    22192222{
    22202223  ofstream ofs;
    22212224  String coordInfo;
    2222   bool hasSameNchan;
    2223   int firstIF;
     2225  bool hasSameNchan = true;
    22242226  bool outTextFile = false;
    22252227
    … …  
    22332235    if (coordInfo == "") coordInfo = "channel";
    22342236    hasSameNchan = hasSameNchanOverIFs();
    2235     firstIF = getIF(0);
    2236   }
    2237 
    2238   Fitter fitter = Fitter();
    2239   fitter.setExpression("poly", order);
     2237  }
     2238
     2239  //Fitter fitter = Fitter();
     2240  //fitter.setExpression("sinusoid", nMaxWavesInSW);
    22402241
    22412242  int nRow = nrow();
    … …  
    22712272
    22722273
    2273     fitBaseline(chanMask, whichrow, fitter);
    2274     setSpectrum(fitter.getResidual(), whichrow);
    2275 
    2276     if (outLogger || outTextFile) {
    2277       std::vector<float> params = fitter.getParameters();
    2278       std::vector<bool>  fixed  = fitter.getFixedParameters();
    2279       float rms = getRms(chanMask, whichrow);
    2280       String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan, firstIF, true);
    2281 
    2282       if (outLogger) {
    2283         LogIO ols(LogOrigin("Scantable", "autoPolyBaseline()", WHERE));
    2284         ols << formatBaselineParams(params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
    2285       }
    2286       if (outTextFile) {
    2287         ofs << formatBaselineParams(params, fixed, rms, masklist, whichrow, true) << flush;
    2288       }
    2289     }
    2290 
     2274    //fitBaseline(chanMask, whichrow, fitter, thresClip, nIterClip);
     2275    //setSpectrum(fitter.getResidual(), whichrow);
     2276    std::vector<int> pieceRanges;
     2277    std::vector<float> params;
     2278    std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nMinWavesInSW, nMaxWavesInSW, params, thresClip, nIterClip);
     2279    setSpectrum(res, whichrow);
     2280    //
     2281
     2282    std::vector<bool> fixed;
     2283    fixed.clear();
     2284    outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoSinusoidBaseline()", pieceRanges, params, fixed);
    22912285  }
    22922286
    … …  
    22942288}
    22952289
     2290std::vector<float> Scantable::doSinusoidFitting(const std::vector<float>& data, const std::vector<bool>& mask, int nMinWavesInSW, int nMaxWavesInSW, std::vector<float>& params, float thresClip, int nIterClip) {
     2291  if (nMaxWavesInSW < 1) {
     2292    throw(AipsError("maximum wave number must be a positive value"));
     2293  }
     2294  if (data.size() != mask.size()) {
     2295    throw(AipsError("data and mask have different size"));
     2296  }
     2297
     2298  int nChan = data.size();
     2299  std::vector<int> maskArray;
     2300  std::vector<int> x;
     2301  for (int i = 0; i < nChan; ++i) {
     2302    maskArray.push_back(mask[i] ? 1 : 0);
     2303    if (mask[i]) {
     2304      x.push_back(i);
     2305    }
     2306  }
     2307
     2308  int nData = x.size();
     2309
     2310  std::vector<double> x1, x2, x3, z1, x1z1, x2z1, x3z1, r1;
     2311  for (int i = 0; i < nChan; ++i) {
     2312    x1.push_back((double)i);
     2313    x2.push_back((double)(i*i));
     2314    x3.push_back((double)(i*i*i));
     2315    z1.push_back((double)data[i]);
     2316    x1z1.push_back(((double)i)*(double)data[i]);
     2317    x2z1.push_back(((double)(i*i))*(double)data[i]);
     2318    x3z1.push_back(((double)(i*i*i))*(double)data[i]);
     2319    r1.push_back(0.0);
     2320  }
     2321
     2322  int currentNData = nData;
     2323  int nDOF = nMaxWavesInSW + 1;  //number of parameters to solve.
     2324
     2325  for (int nClip = 0; nClip < nIterClip+1; ++nClip) {
     2326    //xMatrix : horizontal concatenation of the least-sq. matrix (left) and an
     2327    //identity matrix (right).
     2328    //the right part is used to calculate the inverse matrix of the left part.
     2329    double xMatrix[nDOF][2*nDOF];
     2330    double zMatrix[nDOF];
     2331    for (int i = 0; i < nDOF; ++i) {
     2332      for (int j = 0; j < 2*nDOF; ++j) {
     2333        xMatrix[i][j] = 0.0;
     2334      }
     2335      xMatrix[i][nDOF+i] = 1.0;
     2336      zMatrix[i] = 0.0;
     2337    }
     2338
     2339    for (int i = 0; i < currentNData; ++i) {
     2340      if (maskArray[i] == 0) continue;
     2341      xMatrix[0][0] += 1.0;
     2342      xMatrix[0][1] += x1[i];
     2343      xMatrix[0][2] += x2[i];
     2344      xMatrix[0][3] += x3[i];
     2345      xMatrix[1][1] += x2[i];
     2346      xMatrix[1][2] += x3[i];
     2347      xMatrix[1][3] += x2[i]*x2[i];
     2348      xMatrix[2][2] += x2[i]*x2[i];
     2349      xMatrix[2][3] += x3[i]*x2[i];
     2350      xMatrix[3][3] += x3[i]*x3[i];
     2351      zMatrix[0] += z1[i];
     2352      zMatrix[1] += x1z1[i];
     2353      zMatrix[2] += x2z1[i];
     2354      zMatrix[3] += x3z1[i];
     2355    }
     2356
     2357    for (int i = 0; i < nDOF; ++i) {
     2358      for (int j = 0; j < i; ++j) {
     2359        xMatrix[i][j] = xMatrix[j][i];
     2360      }
     2361    }
     2362
     2363    std::vector<double> invDiag;
     2364    for (int i = 0; i < nDOF; ++i) {
     2365      invDiag.push_back(1.0/xMatrix[i][i]);
     2366      for (int j = 0; j < nDOF; ++j) {
     2367        xMatrix[i][j] *= invDiag[i];
     2368      }
     2369    }
     2370
     2371    for (int k = 0; k < nDOF; ++k) {
     2372      for (int i = 0; i < nDOF; ++i) {
     2373        if (i != k) {
     2374          double factor1 = xMatrix[k][k];
     2375          double factor2 = xMatrix[i][k];
     2376          for (int j = k; j < 2*nDOF; ++j) {
     2377            xMatrix[i][j] *= factor1;
     2378            xMatrix[i][j] -= xMatrix[k][j]*factor2;
     2379            xMatrix[i][j] /= factor1;
     2380          }
     2381        }
     2382      }
     2383      double xDiag = xMatrix[k][k];
     2384      for (int j = k; j < 2*nDOF; ++j) {
     2385        xMatrix[k][j] /= xDiag;
     2386      }
     2387    }
     2388   
     2389    for (int i = 0; i < nDOF; ++i) {
     2390      for (int j = 0; j < nDOF; ++j) {
     2391        xMatrix[i][nDOF+j] *= invDiag[j];
     2392      }
     2393    }
     2394    //compute a vector y which consists of the coefficients of the sinusoids forming the
     2395    //best-fit curves (a0,b0,a1,b1,...), namely, a* and b* are of sine and cosine functions,
     2396    //respectively.
     2397    std::vector<double> y;
     2398    for (int i = 0; i < nDOF; ++i) {
     2399      y.push_back(0.0);
     2400      for (int j = 0; j < nDOF; ++j) {
     2401        y[i] += xMatrix[i][nDOF+j]*zMatrix[j];
     2402      }
     2403    }
     2404
     2405    double a0 = y[0];
     2406    double a1 = y[1];
     2407    double a2 = y[2];
     2408    double a3 = y[3];
     2409    params.clear();
     2410
     2411    for (int i = 0; i < nChan; ++i) {
     2412      r1[i] = a0 + a1*x1[i] + a2*x2[i] + a3*x3[i];
     2413    }
     2414    params.push_back(a0);
     2415    params.push_back(a1);
     2416    params.push_back(a2);
     2417    params.push_back(a3);
     2418
     2419
     2420    if ((nClip == nIterClip) || (thresClip <= 0.0)) {
     2421      break;
     2422    } else {
     2423      std::vector<double> rz;
     2424      double stdDev = 0.0;
     2425      for (int i = 0; i < nChan; ++i) {
     2426        double val = abs(r1[i] - z1[i]);
     2427        rz.push_back(val);
     2428        stdDev += val*val*(double)maskArray[i];
     2429      }
     2430      stdDev = sqrt(stdDev/(double)nData);
     2431     
     2432      double thres = stdDev * thresClip;
     2433      int newNData = 0;
     2434      for (int i = 0; i < nChan; ++i) {
     2435        if (rz[i] >= thres) {
     2436          maskArray[i] = 0;
     2437        }
     2438        if (maskArray[i] > 0) {
     2439          newNData++;
     2440        }
     2441      }
     2442      if (newNData == currentNData) {
     2443        break;                   //no additional flag. finish iteration.
     2444      } else {
     2445        currentNData = newNData;
     2446      }
     2447    }
     2448  }
     2449
     2450  std::vector<float> residual;
     2451  for (int i = 0; i < nChan; ++i) {
     2452    residual.push_back((float)(z1[i] - r1[i]));
     2453  }
     2454  return residual;
     2455
     2456}
     2457
     2458void Scantable::fitBaseline(const std::vector<bool>& mask, int whichrow, Fitter& fitter)
     2459{
     2460  std::vector<double> abcsd = getAbcissa(whichrow);
     2461  std::vector<float> abcs;
     2462  for (uInt i = 0; i < abcsd.size(); ++i) {
     2463    abcs.push_back((float)abcsd[i]);
     2464  }
     2465  std::vector<float> spec = getSpectrum(whichrow);
     2466
     2467  fitter.setData(abcs, spec, mask);
     2468  fitter.lfit();
     2469}
     2470
     2471std::vector<bool> Scantable::getCompositeChanMask(int whichrow, const std::vector<bool>& inMask)
     2472{
     2473  std::vector<bool> chanMask = getMask(whichrow);
     2474  uInt chanMaskSize = chanMask.size();
     2475  if (chanMaskSize != inMask.size()) {
     2476    throw(AipsError("different mask sizes"));
     2477  }
     2478  for (uInt i = 0; i < chanMaskSize; ++i) {
     2479    chanMask[i] = chanMask[i] && inMask[i];
     2480  }
     2481
     2482  return chanMask;
     2483}
     2484
     2485/*
     2486std::vector<bool> Scantable::getCompositeChanMask(int whichrow, const std::vector<bool>& inMask, const std::vector<int>& edge, const int minEdgeSize, STLineFinder& lineFinder)
     2487{
     2488  int edgeSize = edge.size();
     2489  std::vector<int> currentEdge;
     2490  if (edgeSize >= 2) {
     2491      int idx = 0;
     2492      if (edgeSize > 2) {
     2493        if (edgeSize < minEdgeSize) {
     2494          throw(AipsError("Length of edge element info is less than that of IFs"));
     2495        }
     2496        idx = 2 * getIF(whichrow);
     2497      }
     2498      currentEdge.push_back(edge[idx]);
     2499      currentEdge.push_back(edge[idx+1]);
     2500  } else {
     2501    throw(AipsError("Wrong length of edge element"));
     2502  }
     2503
     2504  lineFinder.setData(getSpectrum(whichrow));
     2505  lineFinder.findLines(getCompositeChanMask(whichrow, inMask), currentEdge, whichrow);
     2506
     2507  return lineFinder.getMask();
     2508}
     2509*/
     2510
     2511/* for poly. the variations of outputFittingResult() should be merged into one eventually (2011/3/10 WK)  */
     2512void 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, Fitter& fitter) {
     2513  if (outLogger || outTextFile) {
     2514    std::vector<float> params = fitter.getParameters();
     2515    std::vector<bool>  fixed  = fitter.getFixedParameters();
     2516    float rms = getRms(chanMask, whichrow);
     2517    String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
     2518
     2519    if (outLogger) {
     2520      LogIO ols(LogOrigin("Scantable", funcName, WHERE));
     2521      ols << formatBaselineParams(params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
     2522    }
     2523    if (outTextFile) {
     2524      ofs << formatBaselineParams(params, fixed, rms, masklist, whichrow, true) << flush;
     2525    }
     2526  }
     2527}
     2528
     2529/* for cspline. will be merged once cspline is available in fitter (2011/3/10 WK) */
     2530void 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>& pieceEdges, const std::vector<float>& params, const std::vector<bool>& fixed) {
     2531  if (outLogger || outTextFile) {
     2532    float rms = getRms(chanMask, whichrow);
     2533    String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
     2534
     2535    if (outLogger) {
     2536      LogIO ols(LogOrigin("Scantable", funcName, WHERE));
     2537      ols << formatPiecewiseBaselineParams(pieceEdges, params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
     2538    }
     2539    if (outTextFile) {
     2540      ofs << formatPiecewiseBaselineParams(pieceEdges, params, fixed, rms, masklist, whichrow, true) << flush;
     2541    }
     2542  }
     2543}
     2544
     2545/* for sinusoid. will be merged once sinusoid is available in fitter (2011/3/10 WK) */
     2546void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, const std::vector<float>& params, const std::vector<bool>& fixed) {
     2547  if (outLogger || outTextFile) {
     2548    float rms = getRms(chanMask, whichrow);
     2549    String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
     2550
     2551    if (outLogger) {
     2552      LogIO ols(LogOrigin("Scantable", funcName, WHERE));
     2553      ols << formatBaselineParams(params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
     2554    }
     2555    if (outTextFile) {
     2556      ofs << formatBaselineParams(params, fixed, rms, masklist, whichrow, true) << flush;
     2557    }
     2558  }
     2559}
    22962560
    22972561float Scantable::getRms(const std::vector<bool>& mask, int whichrow) {
    … …  
    23022566  float smean = 0.0;
    23032567  int n = 0;
    2304   for (int i = 0; i < spec.nelements(); ++i) {
     2568  for (uInt i = 0; i < spec.nelements(); ++i) {
    23052569    if (mask[i]) {
    23062570      mean += spec[i];
    … …  
    23442608
    23452609  if (verbose) {
    2346     oss << endl;
    23472610    oss << "Results of baseline fit" << endl;
    23482611    oss << "  rms = " << setprecision(6) << rms << endl;
    … …  
    23522615}
    23532616
    2354 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
     2617std::string Scantable::formatBaselineParams(const std::vector<float>& params, const std::vector<bool>& fixed, float rms, const std::string& masklist, int whichrow, bool verbose) const
    23552618{
    23562619  ostringstream oss;
    23572620  oss << formatBaselineParamsHeader(whichrow, masklist, verbose);
    23582621
     2622  if (params.size() > 0) {
     2623    for (uInt i = 0; i < params.size(); ++i) {
     2624      if (i > 0) {
     2625        oss << ",";
     2626      }
     2627      std::string fix = ((fixed.size() > 0) && (fixed[i]) && verbose) ? "(fixed)" : "";
     2628      float vParam = params[i];
     2629      std::string equals = "= " + (String)((vParam < 0.0) ? "" : " ");
     2630      oss << "  p" << i << fix << equals << setprecision(6) << vParam;
     2631    }
     2632  } else {
     2633    oss << "  Not fitted";
     2634  }
     2635  oss << endl;
     2636
     2637  oss << formatBaselineParamsFooter(rms, verbose);
     2638  oss << flush;
     2639
     2640  return String(oss);
     2641}
     2642
     2643std::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
     2644{
     2645  ostringstream oss;
     2646  oss << formatBaselineParamsHeader(whichrow, masklist, verbose);
     2647
    23592648  if ((params.size() > 0) && (ranges.size() > 0)) {
    23602649    if (((ranges.size() % 2) == 0) && ((params.size() % (ranges.size() / 2)) == 0)) {
    2361       int nParam = params.size() / (ranges.size() / 2);
     2650      uInt nParam = params.size() / (ranges.size() / 2);
     2651      stringstream ss;
     2652      ss << ranges[ranges.size()-1] << flush;
     2653      int wRange = ss.str().size()*2+5;
     2654      ss.str("");
    23622655      for (uInt j = 0; j < ranges.size(); j+=2) {
    2363         oss << "[" << ranges[j] << "," << ranges[j+1] << "]";
     2656        ss << "  [" << ranges[j] << "," << ranges[j+1] << "]" << flush;
     2657        oss << setw(wRange) << left << ss.str();
     2658        ss.str("");
    23642659        for (uInt i = 0; i < nParam; ++i) {
    23652660          if (i > 0) {
    23662661            oss << ",";
    23672662          }
    2368           String fix = ((fixed.size() > 0) && (fixed[i]) && verbose) ? "(fixed)" : "";
    2369           oss << "  p" << i << fix << "= " << setprecision(6) << params[j/2*nParam+i];
     2663          std::string fix = ((fixed.size() > 0) && (fixed[i]) && verbose) ? "(fixed)" : "";
     2664          float vParam = params[j/2*nParam+i];
     2665          std::string equals = "= " + (String)((vParam < 0.0) ? "" : " ");
     2666          oss << "  p" << i << fix << equals << setprecision(6) << vParam;
    23702667        }
     2668        oss << endl;
    23712669      }
    23722670    } else {
    2373       oss << "  ";
     2671      oss << "  " << endl;
    23742672    }
    23752673  } else {
    2376     oss << "  Not fitted";
     2674    oss << "  Not fitted" << endl;
    23772675  }
    23782676
    23792677  oss << formatBaselineParamsFooter(rms, verbose);
    23802678  oss << flush;
    2381 
    2382   return String(oss);
    2383 }
    2384 
    2385 std::string Scantable::formatBaselineParams(const std::vector<float>& params, const std::vector<bool>& fixed, float rms, const std::string& masklist, int whichrow, bool verbose) const
    2386 {
    2387   ostringstream oss;
    2388   oss << formatBaselineParamsHeader(whichrow, masklist, verbose);
    2389 
    2390   if (params.size() > 0) {
    2391     for (uInt i = 0; i < params.size(); ++i) {
    2392       if (i > 0) {
    2393         oss << ",";
    2394       }
    2395       String fix = ((fixed.size() > 0) && (fixed[i]) && verbose) ? "(fixed)" : "";
    2396       oss << "  p" << i << fix << "= " << setprecision(6) << params[i];
    2397     }
    2398   } else {
    2399     oss << "  Not fitted";
    2400   }
    2401 
    2402   oss << formatBaselineParamsFooter(rms, verbose);
    2403   oss << flush;
    2404 
    2405   return String(oss);
    2406 }
    2407 
    2408 std::string Scantable::getMaskRangeList(const std::vector<bool>& mask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, int firstIF, bool silent)
    2409 {
    2410   if (mask.size() < 2) {
    2411     throw(AipsError("The mask elements should be > 1"));
    2412   }
    2413   if (mask.size() != nchan()) {
    2414     throw(AipsError("Number of channels in scantable != number of mask elements"));
    2415   }
    2416 
    2417   std::vector<int> startIdx = getMaskEdgeIndices(mask, true);
    2418   std::vector<int> endIdx   = getMaskEdgeIndices(mask, false);
    2419 
    2420   if (startIdx.size() != endIdx.size()) {
    2421     throw(AipsError("Numbers of mask start and mask end are not identical"));
    2422   }
    2423   for (int i = 0; i < startIdx.size(); ++i) {
    2424     if (startIdx[i] > endIdx[i]) {
    2425       throw(AipsError("Mask start index > mask end index"));
    2426     }
    2427   }
    2428 
    2429   if (!silent) {
    2430     LogIO logOs(LogOrigin("Scantable", "getMaskRangeList()", WHERE));
    2431     logOs << LogIO::WARN << "The current mask window unit is " << coordInfo;
    2432     if (!hasSameNchan) {
    2433       logOs << endl << "This mask is only valid for IF=" << firstIF;
    2434     }
    2435     logOs << LogIO::POST;
    2436   }
    2437 
    2438   std::vector<double> abcissa = getAbcissa(whichrow);
    2439   ostringstream oss;
    2440   oss.setf(ios::fixed);
    2441   oss << setprecision(1) << "[";
    2442   for (int i = 0; i < startIdx.size(); ++i) {
    2443     std::vector<float> aMaskRange;
    2444     aMaskRange.push_back((float)abcissa[startIdx[i]]);
    2445     aMaskRange.push_back((float)abcissa[endIdx[i]]);
    2446     sort(aMaskRange.begin(), aMaskRange.end());
    2447     if (i > 0) oss << ",";
    2448     oss << "[" << aMaskRange[0] << "," << aMaskRange[1] << "]";
    2449   }
    2450   oss << "]" << flush;
    24512679
    24522680  return String(oss);
    … …  
    24712699}
    24722700
    2473 std::vector<int> Scantable::getMaskEdgeIndices(const std::vector<bool>& mask, bool getStartIndices)
     2701std::string Scantable::getMaskRangeList(const std::vector<bool>& mask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, bool verbose)
    24742702{
    24752703  if (mask.size() < 2) {
    24762704    throw(AipsError("The mask elements should be > 1"));
    24772705  }
    2478   if (mask.size() != nchan()) {
     2706  int IF = getIF(whichrow);
     2707  if (mask.size() != (uInt)nchan(IF)) {
    24792708    throw(AipsError("Number of channels in scantable != number of mask elements"));
    24802709  }
    24812710
    2482   std::vector<int> out;
    2483   int endIdx = mask.size() - 1;
    2484 
    2485   if (getStartIndices) {
    2486     if (mask[0]) {
    2487       out.push_back(0);
    2488     }
    2489     for (int i = 0; i < endIdx; ++i) {
    2490       if ((!mask[i]) && mask[i+1]) {
    2491         out.push_back(i+1);
    2492       }
    2493     }
    2494   } else {
    2495     for (int i = 0; i < endIdx; ++i) {
    2496       if (mask[i] && (!mask[i+1])) {
    2497         out.push_back(i);
    2498       }
    2499     }
    2500     if (mask[endIdx]) {
    2501       out.push_back(endIdx);
    2502     }
     2711  if (verbose) {
     2712    LogIO logOs(LogOrigin("Scantable", "getMaskRangeList()", WHERE));
     2713    logOs << LogIO::WARN << "The current mask window unit is " << coordInfo;
     2714    if (!hasSameNchan) {
     2715      logOs << endl << "This mask is only valid for IF=" << IF;
     2716    }
     2717    logOs << LogIO::POST;
     2718  }
     2719
     2720  std::vector<double> abcissa = getAbcissa(whichrow);
     2721  std::vector<int> edge = getMaskEdgeIndices(mask);
     2722
     2723  ostringstream oss;
     2724  oss.setf(ios::fixed);
     2725  oss << setprecision(1) << "[";
     2726  for (uInt i = 0; i < edge.size(); i+=2) {
     2727    if (i > 0) oss << ",";
     2728    oss << "[" << (float)abcissa[edge[i]] << "," << (float)abcissa[edge[i+1]] << "]";
     2729  }
     2730  oss << "]" << flush;
     2731
     2732  return String(oss);
     2733}
     2734
     2735std::vector<int> Scantable::getMaskEdgeIndices(const std::vector<bool>& mask)
     2736{
     2737  if (mask.size() < 2) {
     2738    throw(AipsError("The mask elements should be > 1"));
     2739  }
     2740
     2741  std::vector<int> out, startIndices, endIndices;
     2742  int maskSize = mask.size();
     2743
     2744  startIndices.clear();
     2745  endIndices.clear();
     2746
     2747  if (mask[0]) {
     2748    startIndices.push_back(0);
     2749  }
     2750  for (int i = 1; i < maskSize; ++i) {
     2751    if ((!mask[i-1]) && mask[i]) {
     2752      startIndices.push_back(i);
     2753    } else if (mask[i-1] && (!mask[i])) {
     2754      endIndices.push_back(i-1);
     2755    }
     2756  }
     2757  if (mask[maskSize-1]) {
     2758    endIndices.push_back(maskSize-1);
     2759  }
     2760
     2761  if (startIndices.size() != endIndices.size()) {
     2762    throw(AipsError("Inconsistent Mask Size: bad data?"));
     2763  }
     2764  for (uInt i = 0; i < startIndices.size(); ++i) {
     2765    if (startIndices[i] > endIndices[i]) {
     2766      throw(AipsError("Mask start index > mask end index"));
     2767    }
     2768  }
     2769
     2770  out.clear();
     2771  for (uInt i = 0; i < startIndices.size(); ++i) {
     2772    out.push_back(startIndices[i]);
     2773    out.push_back(endIndices[i]);
    25032774  }
    25042775
Note: See TracChangeset for help on using the changeset viewer.