Changeset 2773


Ignore:
Timestamp:
02/25/13 15:49:06 (12 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: Yes

Module(s): sd

Description: optimisation/refactoring the baselining functions for scantable.


Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/STBaselineTable.cpp

    r2767 r2773  
    137137}
    138138
     139void STBaselineTable::appenddata(int scanno, int cycleno,
     140                                 int beamno, int ifno, int polno,
     141                                 int freqid, Double time,
     142                                 bool apply,
     143                                 STBaselineFunc::FuncName ftype,
     144                                 int fpar,
     145                                 vector<float> ffpar,
     146                                 Vector<uInt> mask,
     147                                 vector<float> res,
     148                                 float rms,
     149                                 int nchan,
     150                                 float cthres,
     151                                 int citer,
     152                                 float lfthres,
     153                                 int lfavg,
     154                                 vector<int> lfedge)
     155{
     156  vector<int> fparam(1);
     157  fparam[0] = fpar;
     158
     159  appenddata(scanno, cycleno, beamno, ifno, polno, freqid, time,
     160             apply, ftype, fparam, ffpar, mask, res, rms, nchan,
     161             cthres, citer, lfthres, lfavg, lfedge);
     162}
     163
     164void STBaselineTable::appenddata(int scanno, int cycleno,
     165                                 int beamno, int ifno, int polno,
     166                                 int freqid, Double time,
     167                                 bool apply,
     168                                 STBaselineFunc::FuncName ftype,
     169                                 vector<int> fpar,
     170                                 vector<float> ffpar,
     171                                 Vector<uInt> mask,
     172                                 vector<float> res,
     173                                 float rms,
     174                                 int nchan,
     175                                 float cthres,
     176                                 int citer,
     177                                 float lfthres,
     178                                 int lfavg,
     179                                 vector<int> lfedge)
     180{
     181  Vector<Int> fparam(fpar.size());
     182  for (uInt i = 0; i < fpar.size(); ++i) {
     183    fparam[i] = fpar[i];
     184  }
     185  Vector<uInt> edge(lfedge.size());
     186  for (uInt i = 0; i < lfedge.size(); ++i) {
     187    edge[i] = lfedge[i];
     188  }
     189  appenddata(uInt(scanno), uInt(cycleno), uInt(beamno),
     190             uInt(ifno), uInt(polno), uInt(freqid), time,
     191             Bool(apply), ftype, fparam, Vector<Float>(ffpar),
     192             mask, Vector<Float>(res), Float(rms), uInt(nchan),
     193             Float(cthres), uInt(citer),
     194             Float(lfthres), uInt(lfavg), edge);
     195}
     196
    139197void STBaselineTable::appenddata(uInt scanno, uInt cycleno,
    140198                                 uInt beamno, uInt ifno, uInt polno,

  • trunk/src/STBaselineTable.h

    r2767 r2773  
    5959               casa::uInt lfavg,
    6060               casa::Vector<casa::uInt> lfedge);
     61  void appenddata(int scanno, int cycleno,
     62                  int beamno, int ifno, int polno,
     63                  int freqid, casa::Double time,
     64                  bool apply,
     65                  STBaselineFunc::FuncName ftype,
     66                  vector<int> fpar,
     67                  vector<float> ffpar,
     68                  casa::Vector<casa::uInt> mask,
     69                  vector<float> res,
     70                  float rms,
     71                  int nchan,
     72                  float cthres,
     73                  int citer,
     74                  float lfthres,
     75                  int lfavg,
     76                  vector<int> lfedge);
     77  void appenddata(int scanno, int cycleno,
     78                  int beamno, int ifno, int polno,
     79                  int freqid, casa::Double time,
     80                  bool apply,
     81                  STBaselineFunc::FuncName ftype,
     82                  int fpar,
     83                  vector<float> ffpar,
     84                  casa::Vector<casa::uInt> mask,
     85                  vector<float> res,
     86                  float rms,
     87                  int nchan,
     88                  float cthres,
     89                  int citer,
     90                  float lfthres,
     91                  int lfavg,
     92                  vector<int> lfedge);
    6193  void appenddata(casa::uInt scanno, casa::uInt cycleno,
    6294                  casa::uInt beamno, casa::uInt ifno, casa::uInt polno,

  • trunk/src/Scantable.cpp

    r2767 r2773  
    24922492std::vector<std::string> Scantable::applyBaselineTable(const std::string& bltable, const std::string& outbltable, const bool outbltableexists, const bool overwrite)
    24932493{
    2494   STBaselineTable* btin = new STBaselineTable(bltable);
    2495 
    2496   Vector<Bool> applyCol = btin->getApply();
     2494  STBaselineTable btin = STBaselineTable(bltable);
     2495
     2496  Vector<Bool> applyCol = btin.getApply();
    24972497  int nRowBl = applyCol.size();
    24982498  if (nRowBl != nrow()) {
    … …  
    25052505  bool outBaselineTable = ((outbltable != "") && (!outbltableexists || overwrite));
    25062506  bool bltableidentical = (bltable == outbltable);
    2507   STBaselineTable* btout = new STBaselineTable(*this);
    2508   ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(table_, "TIME");
    2509   Vector<Double> timeSecCol = tcol->getColumn();
     2507  STBaselineTable btout = STBaselineTable(*this);
     2508  ROScalarColumn<Double> tcol = ROScalarColumn<Double>(table_, "TIME");
     2509  Vector<Double> timeSecCol = tcol.getColumn();
    25102510
    25112511  for (int whichrow = 0; whichrow < nRowBl; ++whichrow) {
    … …  
    25132513      std::vector<float> spec = getSpectrum(whichrow);
    25142514
    2515       //--channel mask--  maybe the following choices should be available
    2516       std::vector<bool> mask = btin->getMask(whichrow);  //use mask_bltable only (default)
    2517       //std::vector<bool> mask = getMask(whichrow);    // or mask_scantable only
    2518       //std::vector<bool> mask = getCompositeChanMask(whichrow, btin->getMask(whichrow)); // or mask_scantable && mask_bltable
    2519 
    2520       STBaselineFunc::FuncName ftype = btin->getFunctionName(whichrow);
    2521       std::vector<int> fpar = btin->getFuncParam(whichrow);
     2515      std::vector<bool> mask = btin.getMask(whichrow);  //use mask_bltable only
     2516
     2517      STBaselineFunc::FuncName ftype = btin.getFunctionName(whichrow);
     2518      std::vector<int> fpar = btin.getFuncParam(whichrow);
    25222519      std::vector<float> params;
    25232520      float rms;
    … …  
    25292526
    25302527      if (outBaselineTable) {
    2531         Vector<Int> fparam(fpar.size());
    2532         for (uInt j = 0; j < fparam.size(); ++j) {
    2533           fparam[j] = (Int)fpar[j];
    2534         }
    2535         std::vector<int> maskList0;
    2536         maskList0.clear();
    2537         for (uInt j = 0; j < mask.size(); ++j) {
    2538           if (mask[j]) {
    2539             if ((j == 0)||(j == mask.size()-1)) {
    2540               maskList0.push_back(j);
    2541             } else {
    2542               if ((mask[j])&&(!mask[j-1])) {
    2543                 maskList0.push_back(j);
    2544               }
    2545               if ((mask[j])&&(!mask[j+1])) {
    2546                 maskList0.push_back(j);
    2547               }
    2548             }
    2549           }
    2550         }
    2551         Vector<uInt> maskList(maskList0.size());
    2552         for (uInt j = 0; j < maskList0.size(); ++j) {
    2553           maskList[j] = (uInt)maskList0[j];
    2554         }
    2555 
    25562528        if (outbltableexists) {
    25572529          if (overwrite) {
    25582530            if (bltableidentical) {
    2559               btin->setresult(uInt(whichrow), Vector<Float>(params), Float(rms));
     2531              btin.setresult(uInt(whichrow), Vector<Float>(params), Float(rms));
    25602532            } else {
    2561               btout->setresult(uInt(whichrow), Vector<Float>(params), Float(rms));
     2533              btout.setresult(uInt(whichrow), Vector<Float>(params), Float(rms));
    25622534            }
    25632535          }
    25642536        } else {
    2565           btout->appenddata(uInt(getScan(whichrow)), uInt(getCycle(whichrow)),
    2566                             uInt(getBeam(whichrow)), uInt(getIF(whichrow)), uInt(getPol(whichrow)),
    2567                             uInt(0), timeSecCol[whichrow], Bool(true), ftype, fparam,
    2568                             Vector<Float>(), maskList, Vector<Float>(params),
    2569                             Float(rms), uInt(spec.size()), Float(3.0), uInt(0),
    2570                             Float(0.0), uInt(0), Vector<uInt>());
     2537          btout.appenddata(getScan(whichrow), getCycle(whichrow), getBeam(whichrow),
     2538                           getIF(whichrow), getPol(whichrow), 0, timeSecCol[whichrow],
     2539                           true, ftype, fpar, std::vector<float>(),
     2540                           getMaskListFromMask(mask), params, rms, spec.size(),
     2541                           3.0, 0, 0.0, 0, std::vector<int>());
    25712542        }
    25722543      }
    … …  
    25762547  if (outBaselineTable) {
    25772548    if (bltableidentical) {
    2578       btin->save(outbltable);
     2549      btin.save(outbltable);
    25792550    } else {
    2580       btout->save(outbltable);
    2581     }
    2582   }
    2583   delete btin;
    2584   delete btout;
     2551      btout.save(outbltable);
     2552    }
     2553  }
    25852554
    25862555  return res;
    … …  
    26002569  }
    26012570
    2602   STBaselineTable* bt = new STBaselineTable(*this);
    2603   ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(table_, "TIME");
    2604   Vector<Double> timeSecCol = tcol->getColumn();
     2571  STBaselineTable bt = STBaselineTable(*this);
     2572  ROScalarColumn<Double> tcol = ROScalarColumn<Double>(table_, "TIME");
     2573  Vector<Double> timeSecCol = tcol.getColumn();
    26052574
    26062575  if (outBaselineTable && !outbltableexists) {
    26072576    for (int i = 0; i < nRowSt; ++i) {
    2608       bt->appendbasedata(getScan(i), getCycle(i), getBeam(i), getIF(i), getPol(i),
     2577      bt.appendbasedata(getScan(i), getCycle(i), getBeam(i), getIF(i), getPol(i),
    26092578                         0, timeSecCol[i]);
    2610       bt->setApply(i, false);
     2579      bt.setApply(i, false);
    26112580    }
    26122581  }
    … …  
    26402609          fparam[j] = (Int)fpar[j];
    26412610        }
    2642         std::vector<int> maskList0;
    2643         maskList0.clear();
    2644         for (uInt j = 0; j < finalmask.size(); ++j) {
    2645           if (finalmask[j]) {
    2646             if ((j == 0)||(j == finalmask.size()-1)) {
    2647               maskList0.push_back(j);
    2648             } else {
    2649               if ((finalmask[j])&&(!finalmask[j-1])) {
    2650                 maskList0.push_back(j);
    2651               }
    2652               if ((finalmask[j])&&(!finalmask[j+1])) {
    2653                 maskList0.push_back(j);
    2654               }
    2655             }
    2656           }
    2657         }
    2658         Vector<uInt> maskList(maskList0.size());
    2659         for (uInt j = 0; j < maskList0.size(); ++j) {
    2660           maskList[j] = (uInt)maskList0[j];
    2661         }
    2662 
    2663         bt->setdata(uInt(irow),
     2611
     2612        bt.setdata(uInt(irow),
    26642613                    uInt(getScan(irow)), uInt(getCycle(irow)),
    26652614                    uInt(getBeam(irow)), uInt(getIF(irow)), uInt(getPol(irow)),
    26662615                    uInt(0), timeSecCol[irow], Bool(true), ftype, fparam,
    2667                     Vector<Float>(), maskList, Vector<Float>(params),
     2616                    Vector<Float>(), getMaskListFromMask(finalmask), Vector<Float>(params),
    26682617                    Float(rms), uInt(spec.size()), Float(clipth), uInt(clipn),
    26692618                    Float(0.0), uInt(0), Vector<uInt>());
    … …  
    26742623
    26752624  if (outBaselineTable) {
    2676     bt->save(outbltable);
    2677   }
    2678   delete bt;
     2625    bt.save(outbltable);
     2626  }
    26792627
    26802628  return res;
    26812629}
    26822630
    2683 std::vector<float> Scantable::doApplyBaselineTable(std::vector<float>& spec, std::vector<bool>& mask, const STBaselineFunc::FuncName ftype, std::vector<int>& fpar, std::vector<float>& params, float&rms)
     2631std::vector<float> Scantable::doApplyBaselineTable(std::vector<float>& spec,
     2632                                                   std::vector<bool>& mask,
     2633                                                   const STBaselineFunc::FuncName ftype,
     2634                                                   std::vector<int>& fpar,
     2635                                                   std::vector<float>& params,
     2636                                                   float&rms)
    26842637{
    26852638  std::vector<bool> finalmask;
    … …  
    26882641}
    26892642
    2690   std::vector<float> Scantable::doSubtractBaseline(std::vector<float>& spec, std::vector<bool>& mask, const STBaselineFunc::FuncName ftype, std::vector<int>& fpar, std::vector<float>& params, float&rms, std::vector<bool>& finalmask, float clipth, int clipn, bool uself, int irow, float lfth, std::vector<int>& lfedge, int lfavg)
    2691 {
    2692   //---replace mask with (mask AND mask_linefinder).
     2643  std::vector<float> Scantable::doSubtractBaseline(std::vector<float>& spec,
     2644                                                   std::vector<bool>& mask,
     2645                                                   const STBaselineFunc::FuncName ftype,
     2646                                                   std::vector<int>& fpar,
     2647                                                   std::vector<float>& params,
     2648                                                   float&rms,
     2649                                                   std::vector<bool>& finalmask,
     2650                                                   float clipth,
     2651                                                   int clipn,
     2652                                                   bool uself,
     2653                                                   int irow,
     2654                                                   float lfth,
     2655                                                   std::vector<int>& lfedge,
     2656                                                   int lfavg)
     2657{
     2658  //--- replacing mask with (mask AND mask_linefinder).
    26932659  if (uself) {
    26942660    STLineFinder lineFinder = STLineFinder();
    … …  
    28852851}
    28862852
    2887 void Scantable::polyBaseline(const std::vector<bool>& mask, int order, float thresClip, int nIterClip, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile, const std::string& bltable)
    2888 {
    2889   /****
     2853Vector<uInt> Scantable::getMaskListFromMask(const std::vector<bool>& mask)
     2854{
     2855  std::vector<int> masklist;
     2856  masklist.clear();
     2857
     2858  for (uInt i = 0; i < mask.size(); ++i) {
     2859    if (mask[i]) {
     2860      if ((i == 0)||(i == mask.size()-1)) {
     2861        masklist.push_back(i);
     2862      } else {
     2863        if ((mask[i])&&(!mask[i-1])) {
     2864          masklist.push_back(i);
     2865        }
     2866        if ((mask[i])&&(!mask[i+1])) {
     2867          masklist.push_back(i);
     2868        }
     2869      }
     2870    }
     2871  }
     2872
     2873  Vector<uInt> res(masklist.size());
     2874  for (uInt i = 0; i < masklist.size(); ++i) {
     2875    res[i] = (uInt)masklist[i];
     2876  }
     2877
     2878  return res;
     2879}
     2880
     2881void Scantable::initialiseBaselining(const std::string& blfile,
     2882                                     ofstream& ofs,
     2883                                     const bool outLogger,
     2884                                     bool& outTextFile,
     2885                                     bool& csvFormat,
     2886                                     String& coordInfo,
     2887                                     bool& hasSameNchan,
     2888                                     const std::string& progressInfo,
     2889                                     bool& showProgress,
     2890                                     int& minNRow,
     2891                                     Vector<Double>& timeSecCol)
     2892{
     2893  csvFormat = false;
     2894  outTextFile = false;
     2895
     2896  if (blfile != "") {
     2897    csvFormat = (blfile.substr(0, 1) == "T");
     2898    ofs.open(blfile.substr(1).c_str(), ios::out | ios::app);
     2899    if (ofs) outTextFile = true;
     2900  }
     2901
     2902  coordInfo = "";
     2903  hasSameNchan = true;
     2904
     2905  if (outLogger || outTextFile) {
     2906    coordInfo = getCoordInfo()[0];
     2907    if (coordInfo == "") coordInfo = "channel";
     2908    hasSameNchan = hasSameNchanOverIFs();
     2909  }
     2910
     2911  parseProgressInfo(progressInfo, showProgress, minNRow);
     2912
     2913  ROScalarColumn<Double> tcol = ROScalarColumn<Double>(table_, "TIME");
     2914  timeSecCol = tcol.getColumn();
     2915}
     2916
     2917void Scantable::finaliseBaselining(const bool outBaselineTable,
     2918                                   STBaselineTable* pbt,
     2919                                   const string& bltable,
     2920                                   const bool outTextFile,
     2921                                   ofstream& ofs)
     2922{
     2923  if (outBaselineTable) {
     2924    pbt->save(bltable);
     2925  }
     2926
     2927  if (outTextFile) ofs.close();
     2928}
     2929
     2930void Scantable::initLineFinder(const std::vector<int>& edge,
     2931                               const float threshold,
     2932                               const int chanAvgLimit,
     2933                               STLineFinder& lineFinder)
     2934{
     2935  if (edge.size() < getIFNos().size()*2) {
     2936    throw(AipsError("Length of edge element info is less than that of IFs"));
     2937  }
     2938
     2939  lineFinder.setOptions(threshold, 3, chanAvgLimit);
     2940}
     2941
     2942void Scantable::polyBaseline(const std::vector<bool>& mask, int order,
     2943                             float thresClip, int nIterClip,
     2944                             bool getResidual,
     2945                             const std::string& progressInfo,
     2946                             const bool outLogger, const std::string& blfile,
     2947                             const std::string& bltable)
     2948{
     2949  /****/
    28902950  double TimeStart = mathutil::gettimeofday_sec();
    2891   ****/
     2951  /****/
    28922952
    28932953  try {
    28942954    ofstream ofs;
    2895     String coordInfo = "";
    2896     bool hasSameNchan = true;
    2897     bool outTextFile = false;
    2898     bool csvFormat = false;
    2899 
    2900     if (blfile != "") {
    2901       csvFormat = (blfile.substr(0, 1) == "T");
    2902       ofs.open(blfile.substr(1).c_str(), ios::out | ios::app);
    2903       if (ofs) outTextFile = true;
    2904     }
    2905 
    2906     if (outLogger || outTextFile) {
    2907       coordInfo = getCoordInfo()[0];
    2908       if (coordInfo == "") coordInfo = "channel";
    2909       hasSameNchan = hasSameNchanOverIFs();
    2910     }
    2911 
    2912     bool showProgress;
     2955    String coordInfo;
     2956    bool hasSameNchan, outTextFile, csvFormat, showProgress;
    29132957    int minNRow;
    2914     parseProgressInfo(progressInfo, showProgress, minNRow);
    2915 
    29162958    int nRow = nrow();
    2917     std::vector<bool> chanMask;
    2918     std::vector<bool> finalChanMask;
     2959    std::vector<bool> chanMask, finalChanMask;
    29192960    float rms;
    2920 
    2921     //---
    2922     bool outBaselintTable = (bltable != "");
    2923     STBaselineTable* bt = new STBaselineTable(*this);
    2924     ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(table_, "TIME");
    2925     Vector<Double> timeSecCol = tcol->getColumn();
    2926     //---
     2961    bool outBaselineTable = (bltable != "");
     2962    STBaselineTable bt = STBaselineTable(*this);
     2963    Vector<Double> timeSecCol;
     2964
     2965    initialiseBaselining(blfile, ofs, outLogger, outTextFile, csvFormat,
     2966                         coordInfo, hasSameNchan,
     2967                         progressInfo, showProgress, minNRow,
     2968                         timeSecCol);
     2969
     2970    std::vector<int> nChanNos;
     2971    std::vector<std::vector<std::vector<double> > > modelReservoir;
     2972    modelReservoir = getPolynomialModelReservoir(order,
     2973                                                 &Scantable::getNormalPolynomial,
     2974                                                 nChanNos);
    29272975
    29282976    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    29292977      std::vector<float> sp = getSpectrum(whichrow);
    29302978      chanMask = getCompositeChanMask(whichrow, mask);
    2931       std::vector<float> params(order+1);
     2979
     2980      std::vector<float> params;
    29322981      int nClipped = 0;
    2933       std::vector<float> res = doPolynomialFitting(sp, chanMask, order, params, rms, finalChanMask, nClipped, thresClip, nIterClip, getResidual);
    2934 
    2935       //---
    2936       if (outBaselintTable) {
    2937         Vector<Int> fparam(1);
    2938         fparam[0] = order;
    2939         std::vector<int> finalMaskList0;
    2940         finalMaskList0.clear();
    2941         for (uInt i = 0; i < finalChanMask.size(); ++i) {
    2942           if (finalChanMask[i]) {
    2943             if ((i == 0)||(i == finalChanMask.size()-1)) {
    2944               finalMaskList0.push_back(i);
    2945             } else {
    2946               if ((finalChanMask[i])&&(!finalChanMask[i-1])) {
    2947                 finalMaskList0.push_back(i);
    2948               }
    2949               if ((finalChanMask[i])&&(!finalChanMask[i+1])) {
    2950                 finalMaskList0.push_back(i);
    2951               }
    2952             }
    2953           }
    2954         }
    2955         Vector<uInt> finalMaskList(finalMaskList0.size());
    2956         for (uInt i = 0; i < finalMaskList0.size(); ++i) {
    2957           finalMaskList[i] = (uInt)finalMaskList0[i];
    2958         }
    2959 
    2960         bt->appenddata(uInt(getScan(whichrow)), uInt(getCycle(whichrow)), uInt(getBeam(whichrow)), uInt(getIF(whichrow)), uInt(getPol(whichrow)),
    2961                        uInt(0),
    2962                        timeSecCol[whichrow],
    2963                        Bool(true),
    2964                        STBaselineFunc::Polynomial,
    2965                        fparam,
    2966                        Vector<Float>(),
    2967                        finalMaskList,
    2968                        Vector<Float>(params),
    2969                        Float(rms),
    2970                        uInt(sp.size()),
    2971                        Float(thresClip),
    2972                        uInt(nIterClip),
    2973                        Float(0.0),
    2974                        uInt(0),
    2975                        Vector<uInt>());
     2982      std::vector<float> res = doLeastSquareFitting(sp, chanMask,
     2983                                   modelReservoir[getIdxOfNchan(sp.size(), nChanNos)],
     2984                                   params, rms, finalChanMask,
     2985                                   nClipped, thresClip, nIterClip, getResidual);
     2986
     2987      if (outBaselineTable) {
     2988        bt.appenddata(getScan(whichrow), getCycle(whichrow), getBeam(whichrow),
     2989                      getIF(whichrow), getPol(whichrow), 0, timeSecCol[whichrow],
     2990                      true, STBaselineFunc::Polynomial, order, std::vector<float>(),
     2991                      getMaskListFromMask(finalChanMask), params, rms, sp.size(),
     2992                      thresClip, nIterClip, 0.0, 0, std::vector<int>());
    29762993      } else {
    29772994        setSpectrum(res, whichrow);
    29782995      }
    2979       //---
    2980 
    2981       outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "chebyshevBaseline()", params, nClipped);
     2996
     2997      outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow,
     2998                          coordInfo, hasSameNchan, ofs, "polyBaseline()",
     2999                          params, nClipped);
    29823000      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
    29833001    }
    2984    
    2985     //---
    2986     if (outBaselintTable) {
    2987       bt->save(bltable);
    2988     }
    2989     //---
    2990     delete bt;
    2991 
    2992     if (outTextFile) ofs.close();
     3002
     3003    finaliseBaselining(outBaselineTable, &bt, bltable, outTextFile, ofs);
    29933004
    29943005  } catch (...) {
    … …  
    29963007  }
    29973008
    2998   /****
     3009  /****/
    29993010  double TimeEnd = mathutil::gettimeofday_sec();
    30003011  double elapse1 = TimeEnd - TimeStart;
    30013012  std::cout << "poly-new   : " << elapse1 << " (sec.)" << endl;
    3002   ****/
    3003 }
    3004 
    3005 void Scantable::autoPolyBaseline(const std::vector<bool>& mask, int order, 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, const std::string& bltable)
     3013  /****/
     3014}
     3015
     3016void Scantable::autoPolyBaseline(const std::vector<bool>& mask, int order,
     3017                                 float thresClip, int nIterClip,
     3018                                 const std::vector<int>& edge,
     3019                                 float threshold, int chanAvgLimit,
     3020                                 bool getResidual,
     3021                                 const std::string& progressInfo,
     3022                                 const bool outLogger, const std::string& blfile,
     3023                                 const std::string& bltable)
    30063024{
    30073025  try {
    30083026    ofstream ofs;
    3009     String coordInfo = "";
    3010     bool hasSameNchan = true;
    3011     bool outTextFile = false;
    3012     bool csvFormat = false;
    3013 
    3014     if (blfile != "") {
    3015       csvFormat = (blfile.substr(0, 1) == "T");
    3016       ofs.open(blfile.substr(1).c_str(), ios::out | ios::app);
    3017       if (ofs) outTextFile = true;
    3018     }
    3019 
    3020     if (outLogger || outTextFile) {
    3021       coordInfo = getCoordInfo()[0];
    3022       if (coordInfo == "") coordInfo = "channel";
    3023       hasSameNchan = hasSameNchanOverIFs();
    3024     }
    3025 
    3026     bool showProgress;
     3027    String coordInfo;
     3028    bool hasSameNchan, outTextFile, csvFormat, showProgress;
    30273029    int minNRow;
    3028     parseProgressInfo(progressInfo, showProgress, minNRow);
    3029 
    3030     int minEdgeSize = getIFNos().size()*2;
     3030    int nRow = nrow();
     3031    std::vector<bool> chanMask, finalChanMask;
     3032    float rms;
     3033    bool outBaselineTable = (bltable != "");
     3034    STBaselineTable bt = STBaselineTable(*this);
     3035    Vector<Double> timeSecCol;
    30313036    STLineFinder lineFinder = STLineFinder();
    3032     lineFinder.setOptions(threshold, 3, chanAvgLimit);
    3033 
    3034     int nRow = nrow();
    3035     std::vector<bool> chanMask;
    3036     std::vector<bool> finalChanMask;
    3037     float rms;
    3038 
    3039     //---
    3040     bool outBaselintTable = (bltable != "");
    3041     STBaselineTable* bt = new STBaselineTable(*this);
    3042     ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(table_, "TIME");
    3043     Vector<Double> timeSecCol = tcol->getColumn();
    3044     //---
     3037
     3038    initialiseBaselining(blfile, ofs, outLogger, outTextFile, csvFormat,
     3039                         coordInfo, hasSameNchan,
     3040                         progressInfo, showProgress, minNRow,
     3041                         timeSecCol);
     3042
     3043    initLineFinder(edge, threshold, chanAvgLimit, lineFinder);
     3044
     3045    std::vector<int> nChanNos;
     3046    std::vector<std::vector<std::vector<double> > > modelReservoir;
     3047    modelReservoir = getPolynomialModelReservoir(order,
     3048                                                 &Scantable::getNormalPolynomial,
     3049                                                 nChanNos);
    30453050
    30463051    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    30473052      std::vector<float> sp = getSpectrum(whichrow);
    3048       //-- composote given channel mask and flagtra --------
    3049       int edgeSize = edge.size();
    30503053      std::vector<int> currentEdge;
    3051       currentEdge.clear();
    3052       if (edgeSize >= 2) {
    3053         int idx = 0;
    3054         if (edgeSize > 2) {
    3055           if (edgeSize < minEdgeSize) {
    3056             throw(AipsError("Length of edge element info is less than that of IFs"));
    3057           }
    3058           idx = 2 * getIF(whichrow);
    3059         }
    3060         currentEdge.push_back(edge[idx]);
    3061         currentEdge.push_back(edge[idx+1]);
    3062       } else {
    3063         throw(AipsError("Wrong length of edge element"));
    3064       }
    3065       lineFinder.setData(sp);
    3066       lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
    3067       chanMask = lineFinder.getMask();
    3068       //-- composote given channel mask and flagtra END --------
    3069 
    3070       std::vector<float> params(order+1);
     3054      chanMask = getCompositeChanMask(whichrow, mask, edge, currentEdge, lineFinder);
     3055
     3056      std::vector<float> params;
    30713057      int nClipped = 0;
    3072       std::vector<float> res = doPolynomialFitting(sp, chanMask, order, params, rms, finalChanMask, nClipped, thresClip, nIterClip, getResidual);
    3073 
    3074       if (outBaselintTable) {
    3075         Vector<Int> fparam(1);
    3076         fparam[0] = order;
    3077         std::vector<int> finalMaskList0;
    3078         finalMaskList0.clear();
    3079         for (uInt i = 0; i < finalChanMask.size(); ++i) {
    3080           if (finalChanMask[i]) {
    3081             if ((i == 0)||(i == finalChanMask.size()-1)) {
    3082               finalMaskList0.push_back(i);
    3083             } else {
    3084               if ((finalChanMask[i])&&(!finalChanMask[i-1])) {
    3085                 finalMaskList0.push_back(i);
    3086               }
    3087               if ((finalChanMask[i])&&(!finalChanMask[i+1])) {
    3088                 finalMaskList0.push_back(i);
    3089               }
    3090             }
    3091           }
    3092         }
    3093         Vector<uInt> finalMaskList(finalMaskList0.size());
    3094         for (uInt i = 0; i < finalMaskList0.size(); ++i) {
    3095           finalMaskList[i] = (uInt)finalMaskList0[i];
    3096         }
    3097 
    3098         Vector<uInt> cEdge(currentEdge.size());
    3099         for (uInt i = 0; i < currentEdge.size(); ++i) {
    3100           cEdge[i] = currentEdge[i];
    3101         }
    3102 
    3103         bt->appenddata(uInt(getScan(whichrow)), uInt(getCycle(whichrow)), uInt(getBeam(whichrow)), uInt(getIF(whichrow)), uInt(getPol(whichrow)),
    3104                        uInt(0),
    3105                        timeSecCol[whichrow],
    3106                        Bool(true),
    3107                        STBaselineFunc::Polynomial,
    3108                        fparam,
    3109                        Vector<Float>(),
    3110                        finalMaskList,
    3111                        Vector<Float>(params),
    3112                        Float(rms),
    3113                        uInt(sp.size()),
    3114                        Float(thresClip),
    3115                        uInt(nIterClip),
    3116                        Float(threshold),
    3117                        uInt(chanAvgLimit),
    3118                        cEdge);
     3058      std::vector<float> res = doLeastSquareFitting(sp, chanMask,
     3059                                   modelReservoir[getIdxOfNchan(sp.size(), nChanNos)],
     3060                                   params, rms, finalChanMask,
     3061                                   nClipped, thresClip, nIterClip, getResidual);
     3062
     3063      if (outBaselineTable) {
     3064        bt.appenddata(getScan(whichrow), getCycle(whichrow), getBeam(whichrow),
     3065                      getIF(whichrow), getPol(whichrow), 0, timeSecCol[whichrow],
     3066                      true, STBaselineFunc::Polynomial, order, std::vector<float>(),
     3067                      getMaskListFromMask(finalChanMask), params, rms, sp.size(),
     3068                      thresClip, nIterClip, threshold, chanAvgLimit, currentEdge);
    31193069      } else {
    31203070        setSpectrum(res, whichrow);
    31213071      }
    31223072
    3123       outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoChebyshevBaseline()", params, nClipped);
     3073      outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow,
     3074                          coordInfo, hasSameNchan, ofs, "autoPolyBaseline()",
     3075                          params, nClipped);
    31243076      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
    31253077    }
    31263078
    3127     //---
    3128     if (outBaselintTable) {
    3129       bt->save(bltable);
    3130     }
    3131     //---
    3132     delete bt;
    3133 
    3134     if (outTextFile) ofs.close();
     3079    finaliseBaselining(outBaselineTable, &bt, bltable, outTextFile, ofs);
    31353080
    31363081  } catch (...) {
    … …  
    31393084}
    31403085
    3141 void Scantable::chebyshevBaseline(const std::vector<bool>& mask, int order, float thresClip, int nIterClip, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile, const std::string& bltable)
     3086void Scantable::chebyshevBaseline(const std::vector<bool>& mask, int order,
     3087                                  float thresClip, int nIterClip,
     3088                                  bool getResidual,
     3089                                  const std::string& progressInfo,
     3090                                  const bool outLogger, const std::string& blfile,
     3091                                  const std::string& bltable)
    31423092{
    31433093  //
    … …  
    31473097  try {
    31483098    ofstream ofs;
    3149     String coordInfo = "";
    3150     bool hasSameNchan = true;
    3151     bool outTextFile = false;
    3152     bool csvFormat = false;
    3153 
    3154     if (blfile != "") {
    3155       csvFormat = (blfile.substr(0, 1) == "T");
    3156       ofs.open(blfile.substr(1).c_str(), ios::out | ios::app);
    3157       if (ofs) outTextFile = true;
    3158     }
    3159 
    3160     if (outLogger || outTextFile) {
    3161       coordInfo = getCoordInfo()[0];
    3162       if (coordInfo == "") coordInfo = "channel";
    3163       hasSameNchan = hasSameNchanOverIFs();
    3164     }
    3165 
    3166     bool showProgress;
     3099    String coordInfo;
     3100    bool hasSameNchan, outTextFile, csvFormat, showProgress;
    31673101    int minNRow;
    3168     parseProgressInfo(progressInfo, showProgress, minNRow);
    3169 
    31703102    int nRow = nrow();
    3171     std::vector<bool> chanMask;
    3172     std::vector<bool> finalChanMask;
     3103    std::vector<bool> chanMask, finalChanMask;
    31733104    float rms;
    3174 
    3175     //---
    3176     bool outBaselintTable = (bltable != "");
    3177     STBaselineTable* bt = new STBaselineTable(*this);
    3178     ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(table_, "TIME");
    3179     Vector<Double> timeSecCol = tcol->getColumn();
    3180     //---
     3105    bool outBaselineTable = (bltable != "");
     3106    STBaselineTable bt = STBaselineTable(*this);
     3107    Vector<Double> timeSecCol;
     3108
     3109    initialiseBaselining(blfile, ofs, outLogger, outTextFile, csvFormat,
     3110                         coordInfo, hasSameNchan,
     3111                         progressInfo, showProgress, minNRow,
     3112                         timeSecCol);
     3113
     3114    std::vector<int> nChanNos;
     3115    std::vector<std::vector<std::vector<double> > > modelReservoir;
     3116    modelReservoir = getPolynomialModelReservoir(order,
     3117                                                 &Scantable::getChebyshevPolynomial,
     3118                                                 nChanNos);
    31813119
    31823120    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    31833121      std::vector<float> sp = getSpectrum(whichrow);
    31843122      chanMask = getCompositeChanMask(whichrow, mask);
    3185       std::vector<float> params(order+1);
     3123
     3124      std::vector<float> params;
    31863125      int nClipped = 0;
    3187       std::vector<float> res = doChebyshevFitting(sp, chanMask, order, params, rms, finalChanMask, nClipped, thresClip, nIterClip, getResidual);
    3188 
    3189       //---
    3190       if (outBaselintTable) {
    3191         Vector<Int> fparam(1);
    3192         fparam[0] = order;
    3193         std::vector<int> finalMaskList0;
    3194         finalMaskList0.clear();
    3195         for (uInt i = 0; i < finalChanMask.size(); ++i) {
    3196           if (finalChanMask[i]) {
    3197             if ((i == 0)||(i == finalChanMask.size()-1)) {
    3198               finalMaskList0.push_back(i);
    3199             } else {
    3200               if ((finalChanMask[i])&&(!finalChanMask[i-1])) {
    3201                 finalMaskList0.push_back(i);
    3202               }
    3203               if ((finalChanMask[i])&&(!finalChanMask[i+1])) {
    3204                 finalMaskList0.push_back(i);
    3205               }
    3206             }
    3207           }
    3208         }
    3209         Vector<uInt> finalMaskList(finalMaskList0.size());
    3210         for (uInt i = 0; i < finalMaskList0.size(); ++i) {
    3211           finalMaskList[i] = (uInt)finalMaskList0[i];
    3212         }
    3213 
    3214         bt->appenddata(uInt(getScan(whichrow)), uInt(getCycle(whichrow)), uInt(getBeam(whichrow)), uInt(getIF(whichrow)), uInt(getPol(whichrow)),
    3215                        uInt(0),
    3216                        timeSecCol[whichrow],
    3217                        Bool(true),
    3218                        STBaselineFunc::Chebyshev,
    3219                        fparam,
    3220                        Vector<Float>(),
    3221                        finalMaskList,
    3222                        Vector<Float>(params),
    3223                        Float(rms),
    3224                        uInt(sp.size()),
    3225                        Float(thresClip),
    3226                        uInt(nIterClip),
    3227                        Float(0.0),
    3228                        uInt(0),
    3229                        Vector<uInt>());
     3126      std::vector<float> res = doLeastSquareFitting(sp, chanMask,
     3127                                   modelReservoir[getIdxOfNchan(sp.size(), nChanNos)],
     3128                                   params, rms, finalChanMask,
     3129                                   nClipped, thresClip, nIterClip, getResidual);
     3130
     3131      if (outBaselineTable) {
     3132        bt.appenddata(getScan(whichrow), getCycle(whichrow), getBeam(whichrow),
     3133                      getIF(whichrow), getPol(whichrow), 0, timeSecCol[whichrow],
     3134                      true, STBaselineFunc::Chebyshev, order, std::vector<float>(),
     3135                      getMaskListFromMask(finalChanMask), params, rms, sp.size(),
     3136                      thresClip, nIterClip, 0.0, 0, std::vector<int>());
    32303137      } else {
    32313138        setSpectrum(res, whichrow);
    32323139      }
    3233       //---
    3234 
    3235       outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "chebyshevBaseline()", params, nClipped);
     3140
     3141      outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow,
     3142                          coordInfo, hasSameNchan, ofs, "chebyshevBaseline()",
     3143                          params, nClipped);
    32363144      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
    32373145    }
    32383146   
    3239     //---
    3240     if (outBaselintTable) {
    3241       bt->save(bltable);
    3242     }
    3243     //---
    3244     delete bt;
    3245 
    3246     if (outTextFile) ofs.close();
     3147    finaliseBaselining(outBaselineTable, &bt, bltable, outTextFile, ofs);
    32473148
    32483149  } catch (...) {
    … …  
    32553156}
    32563157
    3257 void Scantable::autoChebyshevBaseline(const std::vector<bool>& mask, int order, 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, const std::string& bltable)
     3158void Scantable::autoChebyshevBaseline(const std::vector<bool>& mask, int order,
     3159                                      float thresClip, int nIterClip,
     3160                                      const std::vector<int>& edge,
     3161                                      float threshold, int chanAvgLimit,
     3162                                      bool getResidual,
     3163                                      const std::string& progressInfo,
     3164                                      const bool outLogger, const std::string& blfile,
     3165                                      const std::string& bltable)
    32583166{
    32593167  try {
    32603168    ofstream ofs;
    3261     String coordInfo = "";
    3262     bool hasSameNchan = true;
    3263     bool outTextFile = false;
    3264     bool csvFormat = false;
    3265 
    3266     if (blfile != "") {
    3267       csvFormat = (blfile.substr(0, 1) == "T");
    3268       ofs.open(blfile.substr(1).c_str(), ios::out | ios::app);
    3269       if (ofs) outTextFile = true;
    3270     }
    3271 
    3272     if (outLogger || outTextFile) {
    3273       coordInfo = getCoordInfo()[0];
    3274       if (coordInfo == "") coordInfo = "channel";
    3275       hasSameNchan = hasSameNchanOverIFs();
    3276     }
    3277 
    3278     bool showProgress;
     3169    String coordInfo;
     3170    bool hasSameNchan, outTextFile, csvFormat, showProgress;
    32793171    int minNRow;
    3280     parseProgressInfo(progressInfo, showProgress, minNRow);
    3281 
    3282     int minEdgeSize = getIFNos().size()*2;
     3172    int nRow = nrow();
     3173    std::vector<bool> chanMask, finalChanMask;
     3174    float rms;
     3175    bool outBaselineTable = (bltable != "");
     3176    STBaselineTable bt = STBaselineTable(*this);
     3177    Vector<Double> timeSecCol;
    32833178    STLineFinder lineFinder = STLineFinder();
    3284     lineFinder.setOptions(threshold, 3, chanAvgLimit);
    3285 
    3286     int nRow = nrow();
    3287     std::vector<bool> chanMask;
    3288     std::vector<bool> finalChanMask;
    3289     float rms;
    3290 
    3291     //---
    3292     bool outBaselintTable = (bltable != "");
    3293     STBaselineTable* bt = new STBaselineTable(*this);
    3294     ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(table_, "TIME");
    3295     Vector<Double> timeSecCol = tcol->getColumn();
    3296     //---
     3179
     3180    initialiseBaselining(blfile, ofs, outLogger, outTextFile, csvFormat,
     3181                         coordInfo, hasSameNchan,
     3182                         progressInfo, showProgress, minNRow,
     3183                         timeSecCol);
     3184
     3185    initLineFinder(edge, threshold, chanAvgLimit, lineFinder);
     3186
     3187    std::vector<int> nChanNos;
     3188    std::vector<std::vector<std::vector<double> > > modelReservoir;
     3189    modelReservoir = getPolynomialModelReservoir(order,
     3190                                                 &Scantable::getChebyshevPolynomial,
     3191                                                 nChanNos);
    32973192
    32983193    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    32993194      std::vector<float> sp = getSpectrum(whichrow);
    3300       //-- composote given channel mask and flagtra --------
    3301       int edgeSize = edge.size();
    33023195      std::vector<int> currentEdge;
    3303       currentEdge.clear();
    3304       if (edgeSize >= 2) {
    3305         int idx = 0;
    3306         if (edgeSize > 2) {
    3307           if (edgeSize < minEdgeSize) {
    3308             throw(AipsError("Length of edge element info is less than that of IFs"));
    3309           }
    3310           idx = 2 * getIF(whichrow);
    3311         }
    3312         currentEdge.push_back(edge[idx]);
    3313         currentEdge.push_back(edge[idx+1]);
    3314       } else {
    3315         throw(AipsError("Wrong length of edge element"));
    3316       }
    3317       lineFinder.setData(sp);
    3318       lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
    3319       chanMask = lineFinder.getMask();
    3320       //-- composote given channel mask and flagtra END --------
    3321 
    3322       std::vector<float> params(order+1);
     3196      chanMask = getCompositeChanMask(whichrow, mask, edge, currentEdge, lineFinder);
     3197
     3198      std::vector<float> params;
    33233199      int nClipped = 0;
    3324       std::vector<float> res = doChebyshevFitting(sp, chanMask, order, params, rms, finalChanMask, nClipped, thresClip, nIterClip, getResidual);
    3325 
    3326       if (outBaselintTable) {
    3327         Vector<Int> fparam(1);
    3328         fparam[0] = order;
    3329         std::vector<int> finalMaskList0;
    3330         finalMaskList0.clear();
    3331         for (uInt i = 0; i < finalChanMask.size(); ++i) {
    3332           if (finalChanMask[i]) {
    3333             if ((i == 0)||(i == finalChanMask.size()-1)) {
    3334               finalMaskList0.push_back(i);
    3335             } else {
    3336               if ((finalChanMask[i])&&(!finalChanMask[i-1])) {
    3337                 finalMaskList0.push_back(i);
    3338               }
    3339               if ((finalChanMask[i])&&(!finalChanMask[i+1])) {
    3340                 finalMaskList0.push_back(i);
    3341               }
    3342             }
    3343           }
    3344         }
    3345         Vector<uInt> finalMaskList(finalMaskList0.size());
    3346         for (uInt i = 0; i < finalMaskList0.size(); ++i) {
    3347           finalMaskList[i] = (uInt)finalMaskList0[i];
    3348         }
    3349 
    3350         Vector<uInt> cEdge(currentEdge.size());
    3351         for (uInt i = 0; i < currentEdge.size(); ++i) {
    3352           cEdge[i] = currentEdge[i];
    3353         }
    3354 
    3355         bt->appenddata(uInt(getScan(whichrow)), uInt(getCycle(whichrow)), uInt(getBeam(whichrow)), uInt(getIF(whichrow)), uInt(getPol(whichrow)),
    3356                        uInt(0),
    3357                        timeSecCol[whichrow],
    3358                        Bool(true),
    3359                        STBaselineFunc::Chebyshev,
    3360                        fparam,
    3361                        Vector<Float>(),
    3362                        finalMaskList,
    3363                        Vector<Float>(params),
    3364                        Float(rms),
    3365                        uInt(sp.size()),
    3366                        Float(thresClip),
    3367                        uInt(nIterClip),
    3368                        Float(threshold),
    3369                        uInt(chanAvgLimit),
    3370                        cEdge);
     3200      std::vector<float> res = doLeastSquareFitting(sp, chanMask,
     3201                                   modelReservoir[getIdxOfNchan(sp.size(), nChanNos)],
     3202                                   params, rms, finalChanMask,
     3203                                   nClipped, thresClip, nIterClip, getResidual);
     3204
     3205      if (outBaselineTable) {
     3206        bt.appenddata(getScan(whichrow), getCycle(whichrow), getBeam(whichrow),
     3207                      getIF(whichrow), getPol(whichrow), 0, timeSecCol[whichrow],
     3208                      true, STBaselineFunc::Chebyshev, order, std::vector<float>(),
     3209                      getMaskListFromMask(finalChanMask), params, rms, sp.size(),
     3210                      thresClip, nIterClip, threshold, chanAvgLimit, currentEdge);
    33713211      } else {
    33723212        setSpectrum(res, whichrow);
    33733213      }
    33743214
    3375       outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoChebyshevBaseline()", params, nClipped);
     3215      outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow,
     3216                          coordInfo, hasSameNchan, ofs, "autoChebyshevBaseline()",
     3217                          params, nClipped);
    33763218      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
    33773219    }
    33783220
    3379     //---
    3380     if (outBaselintTable) {
    3381       bt->save(bltable);
    3382     }
    3383     //---
    3384     delete bt;
    3385 
    3386     if (outTextFile) ofs.close();
     3221    finaliseBaselining(outBaselineTable, &bt, bltable, outTextFile, ofs);
    33873222
    33883223  } catch (...) {
    … …  
    34483283    std::vector<int> pieceEdges;//(order+1);  //order = npiece
    34493284    nparam = order + 3;
    3450     //params.resize(4*order);
    34513285    res = doCubicSplineFitting(spec, mask, order, false, pieceEdges, params, rms, finalChanMask, nClipped);
    34523286  } else if (blfunc == "sinusoid") {
    … …  
    35683402}
    35693403
    3570 std::vector<float> Scantable::doPolynomialFitting(const std::vector<float>& data, const std::vector<bool>& mask, int order, std::vector<float>& params, float& rms, std::vector<bool>& finalmask, float clipth, int clipn)
     3404std::vector<float> Scantable::doPolynomialFitting(const std::vector<float>& data,
     3405                                                  const std::vector<bool>& mask,
     3406                                                  int order,
     3407                                                  std::vector<float>& params,
     3408                                                  float& rms,
     3409                                                  std::vector<bool>& finalmask,
     3410                                                  float clipth,
     3411                                                  int clipn)
    35713412{
    35723413  int nClipped = 0;
    … …  
    35743415}
    35753416
    3576 std::vector<float> Scantable::doPolynomialFitting(const std::vector<float>& data, const std::vector<bool>& mask, int order, std::vector<float>& params, float& rms, std::vector<bool>& finalMask, int& nClipped, float thresClip, int nIterClip, bool getResidual)
    3577 {
    3578   return doLeastSquareFitting(data, mask, &Scantable::getNormalPolynomial, order, params, rms, finalMask, nClipped, thresClip, nIterClip, getResidual);
    3579 }
    3580 
    3581 std::vector<float> Scantable::doChebyshevFitting(const std::vector<float>& data, const std::vector<bool>& mask, int order, std::vector<float>& params, float& rms, std::vector<bool>& finalmask, float clipth, int clipn)
     3417std::vector<float> Scantable::doPolynomialFitting(const std::vector<float>& data,
     3418                                                  const std::vector<bool>& mask,
     3419                                                  int order,
     3420                                                  std::vector<float>& params,
     3421                                                  float& rms,
     3422                                                  std::vector<bool>& finalMask,
     3423                                                  int& nClipped,
     3424                                                  float thresClip,
     3425                                                  int nIterClip,
     3426                                                  bool getResidual)
     3427{
     3428  return doLeastSquareFitting(data, mask,
     3429                              getPolynomialModel(order, data.size(), &Scantable::getNormalPolynomial),
     3430                              params, rms, finalMask,
     3431                              nClipped, thresClip, nIterClip,
     3432                              getResidual);
     3433}
     3434
     3435std::vector<float> Scantable::doChebyshevFitting(const std::vector<float>& data,
     3436                                                 const std::vector<bool>& mask,
     3437                                                 int order,
     3438                                                 std::vector<float>& params,
     3439                                                 float& rms,
     3440                                                 std::vector<bool>& finalmask,
     3441                                                 float clipth,
     3442                                                 int clipn)
    35823443{
    35833444  int nClipped = 0;
    … …  
    35853446}
    35863447
    3587 std::vector<float> Scantable::doChebyshevFitting(const std::vector<float>& data, const std::vector<bool>& mask, int order, std::vector<float>& params, float& rms, std::vector<bool>& finalMask, int& nClipped, float thresClip, int nIterClip, bool getResidual)
    3588 {
    3589   return doLeastSquareFitting(data, mask, &Scantable::getChebyshevPolynomial, order, params, rms, finalMask, nClipped, thresClip, nIterClip, getResidual);
    3590 }
    3591 
    3592 std::vector<float> Scantable::doLeastSquareFitting(const std::vector<float>& data, const std::vector<bool>& mask, double (Scantable::*pfunc)(int, double), int order, std::vector<float>& params, float& rms, std::vector<bool>& finalMask, int& nClipped, float thresClip, int nIterClip, bool getResidual)
    3593 {
    3594   if (data.size() != mask.size()) {
    3595     throw(AipsError("data and mask sizes are not identical"));
    3596   }
    3597   if (order < 0) {
    3598     throw(AipsError("maximum order of polynomial must not be negative."));
    3599   }
    3600 
    3601   int nChan = data.size();
    3602 
    3603   finalMask.clear();
    3604   finalMask.resize(nChan);
    3605 
    3606   std::vector<int> maskArray;
    3607   std::vector<int> x;
    3608   for (int i = 0; i < nChan; ++i) {
    3609     maskArray.push_back(mask[i] ? 1 : 0);
    3610     if (mask[i]) {
    3611       x.push_back(i);
    3612     }
    3613     finalMask[i] = mask[i];
    3614   }
    3615 
    3616   int initNData = x.size();
    3617 
    3618   int nData = initNData;
    3619   int nDOF = order + 1;  //number of parameters to solve.
    3620 
    3621   // xArray : contains elemental values for computing the least-square matrix.
    3622   //          xArray.size() is nDOF and xArray[*].size() is nChan.
    3623   //          Each xArray element are as follows:
     3448std::vector<float> Scantable::doChebyshevFitting(const std::vector<float>& data,
     3449                                                 const std::vector<bool>& mask,
     3450                                                 int order,
     3451                                                 std::vector<float>& params,
     3452                                                 float& rms,
     3453                                                 std::vector<bool>& finalMask,
     3454                                                 int& nClipped,
     3455                                                 float thresClip,
     3456                                                 int nIterClip,
     3457                                                 bool getResidual)
     3458{
     3459  return doLeastSquareFitting(data, mask,
     3460                              getPolynomialModel(order, data.size(), &Scantable::getChebyshevPolynomial),
     3461                              params, rms, finalMask,
     3462                              nClipped, thresClip, nIterClip,
     3463                              getResidual);
     3464}
     3465
     3466std::vector<std::vector<double> > Scantable::getPolynomialModel(int order, int nchan, double (Scantable::*pfunc)(int, double))
     3467{
     3468  // model  : contains model values for computing the least-square matrix.
     3469  //          model.size() is nmodel and model[*].size() is nchan.
     3470  //          Each model element are as follows:
    36243471  //
    36253472  //          (for normal polynomials)
    3626   //          xArray[0]   = {1.0,   1.0,   1.0,   ..., 1.0},
    3627   //          xArray[1]   = {0.0,   1.0,   2.0,   ..., (nChan-1)}
    3628   //          xArray[n-1] = ...,
    3629   //          xArray[n]   = {0.0^n, 1.0^n, 2.0^n, ..., (nChan-1)^n}
     3473  //          model[0]   = {1.0,   1.0,   1.0,   ..., 1.0},
     3474  //          model[1]   = {0.0,   1.0,   2.0,   ..., (nchan-1)}
     3475  //          model[n-1] = ...,
     3476  //          model[n]   = {0.0^n, 1.0^n, 2.0^n, ..., (nchan-1)^n}
    36303477  //          where (0 <= n <= order)
    36313478  //
    36323479  //          (for Chebyshev polynomials)
    3633   //          xArray[0]   = {T0(-1), T0(2/(nChan-1)-1), T0(4/(nChan-1)-1), ..., T0(1)},
    3634   //          xArray[n-1] = ...,
    3635   //          xArray[n]   = {Tn(-1), Tn(2/(nChan-1)-1), Tn(4/(nChan-1)-1), ..., Tn(1)}
     3480  //          model[0]   = {T0(-1), T0(2/(nchan-1)-1), T0(4/(nchan-1)-1), ..., T0(1)},
     3481  //          model[n-1] = ...,
     3482  //          model[n]   = {Tn(-1), Tn(2/(nchan-1)-1), Tn(4/(nchan-1)-1), ..., Tn(1)}
    36363483  //          where (0 <= n <= order),
    3637   std::vector<std::vector<double> > xArray;
    3638   double xFactor, xShift;
     3484
     3485  int nmodel = order + 1;
     3486  std::vector<std::vector<double> > model(nmodel, std::vector<double>(nchan));
     3487
     3488  double stretch, shift;
    36393489  if (pfunc == &Scantable::getChebyshevPolynomial) {
    3640     xFactor = 2.0/(double)(nChan - 1);
    3641     xShift  = -1.0;
     3490    stretch = 2.0/(double)(nchan - 1);
     3491    shift   = -1.0;
    36423492  } else {
    3643     xFactor = 1.0;
    3644     xShift  = 0.0;
     3493    stretch = 1.0;
     3494    shift   = 0.0;
     3495  }
     3496
     3497  for (int i = 0; i < nmodel; ++i) {
     3498    for (int j = 0; j < nchan; ++j) {
     3499      model[i][j] = (this->*pfunc)(i, stretch*(double)j + shift);
     3500    }
     3501  }
     3502
     3503  return model;
     3504}
     3505
     3506std::vector<std::vector<std::vector<double> > > Scantable::getPolynomialModelReservoir(int order,
     3507                                                                                       double (Scantable::*pfunc)(int, double),
     3508                                                                                       std::vector<int>& nChanNos)
     3509{
     3510  std::vector<std::vector<std::vector<double> > > res;
     3511  res.clear();
     3512  nChanNos.clear();
     3513
     3514  std::vector<uint> ifNos = getIFNos();
     3515  for (uint i = 0; i < ifNos.size(); ++i) {
     3516    int currNchan = nchan(ifNos[i]);
     3517    bool hasDifferentNchan = (i == 0);
     3518    for (uint j = 0; j < i; ++j) {
     3519      if (currNchan != nchan(ifNos[j])) {
     3520        hasDifferentNchan = true;
     3521        break;
     3522      }
     3523    }
     3524    if (hasDifferentNchan) {
     3525      res.push_back(getPolynomialModel(order, currNchan, pfunc));
     3526      nChanNos.push_back(currNchan);
     3527    }
     3528  }
     3529
     3530  return res;
     3531}
     3532
     3533std::vector<float> Scantable::doLeastSquareFitting(const std::vector<float>& data,
     3534                                                   const std::vector<bool>& mask,
     3535                                                   const std::vector<std::vector<double> >& model,
     3536                                                   std::vector<float>& params,
     3537                                                   float& rms,
     3538                                                   std::vector<bool>& finalMask,
     3539                                                   int& nClipped,
     3540                                                   float thresClip,
     3541                                                   int nIterClip,
     3542                                                   bool getResidual)
     3543{
     3544  int nDOF = model.size();
     3545  int nChan = data.size();
     3546
     3547  if (nDOF == 0) {
     3548    throw(AipsError("no model data given"));
     3549  }
     3550  if (nChan < 2) {
     3551    throw(AipsError("data size is too few"));
     3552  }
     3553  if (nChan != (int)mask.size()) {
     3554    throw(AipsError("data and mask sizes are not identical"));
    36453555  }
    36463556  for (int i = 0; i < nDOF; ++i) {
    3647     std::vector<double> xs;
    3648     xs.clear();
    3649     for (int j = 0; j < nChan; ++j) {
    3650       xs.push_back((this->*pfunc)(i, xFactor*(double)j + xShift));
    3651     }
    3652     xArray.push_back(xs);
    3653   }
    3654 
    3655   std::vector<double> z1, r1, residual;
     3557    if (nChan != (int)model[i].size()) {
     3558      throw(AipsError("data and model sizes are not identical"));
     3559    }
     3560  }
     3561
     3562  params.clear();
     3563  params.resize(nDOF);
     3564
     3565  finalMask.clear();
     3566  finalMask.resize(nChan);
     3567
     3568  std::vector<int> maskArray(nChan);
     3569  int j = 0;
    36563570  for (int i = 0; i < nChan; ++i) {
    3657     z1.push_back((double)data[i]);
    3658     r1.push_back(0.0);
    3659     residual.push_back(0.0);
     3571    maskArray[i] = mask[i] ? 1 : 0;
     3572    if (mask[i]) {
     3573      j++;
     3574    }
     3575    finalMask[i] = mask[i];
     3576  }
     3577
     3578  int initNData = j;
     3579  int nData = initNData;
     3580
     3581  std::vector<double> z1(nChan), r1(nChan), residual(nChan);
     3582  for (int i = 0; i < nChan; ++i) {
     3583    z1[i] = (double)data[i];
     3584    r1[i] = 0.0;
     3585    residual[i] = 0.0;
    36603586  }
    36613587
    … …  
    36813607      for (int i = 0; i < nDOF; ++i) {
    36823608        for (int j = i; j < nDOF; ++j) {
    3683           xMatrix[i][j] += xArray[i][k] * xArray[j][k];
     3609          xMatrix[i][j] += model[i][k] * model[j][k];
    36843610        }
    3685         zMatrix[i] += z1[k] * xArray[i][k];
     3611        zMatrix[i] += z1[k] * model[i][k];
    36863612      }
    36873613
    … …  
    36993625    }
    37003626
    3701     std::vector<double> invDiag;
     3627    std::vector<double> invDiag(nDOF);
    37023628    for (int i = 0; i < nDOF; ++i) {
    3703       invDiag.push_back(1.0/xMatrix[i][i]);
     3629      invDiag[i] = 1.0 / xMatrix[i][i];
    37043630      for (int j = 0; j < nDOF; ++j) {
    37053631        xMatrix[i][j] *= invDiag[i];
    … …  
    37113637        if (i != k) {
    37123638          double factor1 = xMatrix[k][k];
     3639          double invfactor1 = 1.0 / factor1;
    37133640          double factor2 = xMatrix[i][k];
    37143641          for (int j = k; j < 2*nDOF; ++j) {
    37153642            xMatrix[i][j] *= factor1;
    37163643            xMatrix[i][j] -= xMatrix[k][j]*factor2;
    3717             xMatrix[i][j] /= factor1;
     3644            xMatrix[i][j] *= invfactor1;
    37183645          }
    37193646        }
    37203647      }
    3721       double xDiag = xMatrix[k][k];
     3648      double invXDiag = 1.0 / xMatrix[k][k];
    37223649      for (int j = k; j < 2*nDOF; ++j) {
    3723         xMatrix[k][j] /= xDiag;
     3650        xMatrix[k][j] *= invXDiag;
    37243651      }
    37253652    }
    … …  
    37303657      }
    37313658    }
    3732     //compute a vector y which consists of the coefficients of Chebyshev polynomials.
    3733     std::vector<double> y;
    3734     params.clear();
     3659    //compute a vector y in which coefficients of the best-fit
     3660    //model functions are stored.
     3661    //in case of polynomials, y consists of (a0,a1,a2,...)
     3662    //where ai is the coefficient of the term x^i.
     3663    //in case of sinusoids, y consists of (a0,s1,c1,s2,c2,...)
     3664    //where a0 is constant term and s* and c* are of sine
     3665    //and cosine functions, respectively.
     3666    std::vector<double> y(nDOF);
    37353667    for (int i = 0; i < nDOF; ++i) {
    3736       y.push_back(0.0);
     3668      y[i] = 0.0;
    37373669      for (int j = 0; j < nDOF; ++j) {
    37383670        y[i] += xMatrix[i][nDOF+j]*zMatrix[j];
    37393671      }
    3740       params.push_back(y[i]);
     3672      params[i] = (float)y[i];
    37413673    }
    37423674
    … …  
    37443676      r1[i] = y[0];
    37453677      for (int j = 1; j < nDOF; ++j) {
    3746         r1[i] += y[j]*xArray[j][i];
     3678        r1[i] += y[j]*model[j][i];
    37473679      }
    37483680      residual[i] = z1[i] - r1[i];
    … …  
    37823714  nClipped = initNData - nData;
    37833715
    3784   std::vector<float> result;
    3785   result.clear();
     3716  std::vector<float> result(nChan);
    37863717  if (getResidual) {
    37873718    for (int i = 0; i < nChan; ++i) {
    3788       result.push_back((float)residual[i]);
     3719      result[i] = (float)residual[i];
    37893720    }
    37903721  } else {
    37913722    for (int i = 0; i < nChan; ++i) {
    3792       result.push_back((float)r1[i]);
     3723      result[i] = (float)r1[i];
    37933724    }
    37943725  }
    … …  
    37973728}
    37983729
    3799 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, const std::string& bltable)
    3800 {
     3730void Scantable::cubicSplineBaseline(const std::vector<bool>& mask, int nPiece,
     3731                                    float thresClip, int nIterClip,
     3732                                    bool getResidual,
     3733                                    const std::string& progressInfo,
     3734                                    const bool outLogger, const std::string& blfile,
     3735                                    const std::string& bltable)
     3736{
     3737  /****/
     3738  double TimeStart = mathutil::gettimeofday_sec();
     3739  /****/
     3740
    38013741  try {
    38023742    ofstream ofs;
    3803     String coordInfo = "";
    3804     bool hasSameNchan = true;
    3805     bool outTextFile = false;
    3806     bool csvFormat = false;
    3807 
    3808     if (blfile != "") {
    3809       csvFormat = (blfile.substr(0, 1) == "T");
    3810       ofs.open(blfile.substr(1).c_str(), ios::out | ios::app);
    3811       if (ofs) outTextFile = true;
    3812     }
    3813 
    3814     if (outLogger || outTextFile) {
    3815       coordInfo = getCoordInfo()[0];
    3816       if (coordInfo == "") coordInfo = "channel";
    3817       hasSameNchan = hasSameNchanOverIFs();
    3818     }
    3819 
    3820     bool showProgress;
     3743    String coordInfo;
     3744    bool hasSameNchan, outTextFile, csvFormat, showProgress;
    38213745    int minNRow;
    3822     parseProgressInfo(progressInfo, showProgress, minNRow);
    3823 
    38243746    int nRow = nrow();
    3825     std::vector<bool> chanMask;
    3826     std::vector<bool> finalChanMask;
     3747    std::vector<bool> chanMask, finalChanMask;
    38273748    float rms;
    3828 
    3829     //---
    3830     bool outBaselintTable = (bltable != "");
    3831     STBaselineTable* bt = new STBaselineTable(*this);
    3832     ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(table_, "TIME");
    3833     Vector<Double> timeSecCol = tcol->getColumn();
    3834     //---
    3835 
     3749    bool outBaselineTable = (bltable != "");
     3750    STBaselineTable bt = STBaselineTable(*this);
     3751    Vector<Double> timeSecCol;
     3752
     3753    initialiseBaselining(blfile, ofs, outLogger, outTextFile, csvFormat,
     3754                         coordInfo, hasSameNchan,
     3755                         progressInfo, showProgress, minNRow,
     3756                         timeSecCol);
     3757
     3758    std::vector<int> nChanNos;
     3759    std::vector<std::vector<std::vector<double> > > modelReservoir;
     3760    modelReservoir = getPolynomialModelReservoir(3,
     3761                                                 &Scantable::getNormalPolynomial,
     3762                                                 nChanNos);
    38363763
    38373764    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    38383765      std::vector<float> sp = getSpectrum(whichrow);
    38393766      chanMask = getCompositeChanMask(whichrow, mask);
    3840       std::vector<int> pieceEdges;//(nPiece+1);
    3841       std::vector<float> params;//(nPiece*4);
     3767
     3768      std::vector<int> pieceEdges;
     3769      std::vector<float> params;
    38423770      int nClipped = 0;
    3843       std::vector<float> res = doCubicSplineFitting(sp, chanMask, nPiece, false, pieceEdges, params, rms, finalChanMask, nClipped, thresClip, nIterClip, getResidual);
    3844 
    3845       //---
    3846       if (outBaselintTable) {
    3847         Vector<Int> fparam(nPiece);
    3848         for (uInt i = 0; i < fparam.size(); ++i) {
    3849           fparam[i] = pieceEdges[i];
    3850         }
    3851         std::vector<int> finalMaskList0;
    3852         finalMaskList0.clear();
    3853         for (uInt i = 0; i < finalChanMask.size(); ++i) {
    3854           if (finalChanMask[i]) {
    3855             if ((i == 0)||(i == finalChanMask.size()-1)) {
    3856               finalMaskList0.push_back(i);
    3857             } else {
    3858               if ((finalChanMask[i])&&(!finalChanMask[i-1])) {
    3859                 finalMaskList0.push_back(i);
    3860               }
    3861               if ((finalChanMask[i])&&(!finalChanMask[i+1])) {
    3862                 finalMaskList0.push_back(i);
    3863               }
    3864             }
    3865           }
    3866         }
    3867         Vector<uInt> finalMaskList(finalMaskList0.size());
    3868         for (uInt i = 0; i < finalMaskList0.size(); ++i) {
    3869           finalMaskList[i] = (uInt)finalMaskList0[i];
    3870         }
    3871 
    3872         bt->appenddata(uInt(getScan(whichrow)), uInt(getCycle(whichrow)), uInt(getBeam(whichrow)), uInt(getIF(whichrow)), uInt(getPol(whichrow)),
    3873                        uInt(0),
    3874                        timeSecCol[whichrow],
    3875                        Bool(true),
    3876                        STBaselineFunc::CSpline,
    3877                        fparam,
    3878                        Vector<Float>(),
    3879                        finalMaskList,
    3880                        Vector<Float>(params),
    3881                        Float(rms),
    3882                        uInt(sp.size()),
    3883                        Float(thresClip),
    3884                        uInt(nIterClip),
    3885                        Float(0.0),
    3886                        uInt(0),
    3887                        Vector<uInt>());
     3771      std::vector<float> res = doCubicSplineLeastSquareFitting(sp, chanMask,
     3772                                   modelReservoir[getIdxOfNchan(sp.size(), nChanNos)],
     3773                                   nPiece, false, pieceEdges, params, rms, finalChanMask,
     3774                                   nClipped, thresClip, nIterClip, getResidual);
     3775
     3776      if (outBaselineTable) {
     3777        bt.appenddata(getScan(whichrow), getCycle(whichrow), getBeam(whichrow),
     3778                      getIF(whichrow), getPol(whichrow), 0, timeSecCol[whichrow],
     3779                      true, STBaselineFunc::CSpline, pieceEdges, std::vector<float>(),
     3780                      getMaskListFromMask(finalChanMask), params, rms, sp.size(),
     3781                      thresClip, nIterClip, 0.0, 0, std::vector<int>());
    38883782      } else {
    38893783        setSpectrum(res, whichrow);
    38903784      }
    3891       //---
    3892 
    3893       outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()", pieceEdges, params, nClipped);
     3785
     3786      outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow,
     3787                          coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()",
     3788                          pieceEdges, params, nClipped);
    38943789      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
    38953790    }
    38963791   
    3897     //---
    3898     if (outBaselintTable) {
    3899       bt->save(bltable);
    3900     }
    3901     //---
    3902     delete bt;
    3903 
    3904     if (outTextFile) ofs.close();
     3792    finaliseBaselining(outBaselineTable, &bt, bltable, outTextFile, ofs);
    39053793
    39063794  } catch (...) {
    39073795    throw;
    39083796  }
    3909 }
    3910 
    3911 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, const std::string& bltable)
     3797
     3798  /****/
     3799  double TimeEnd = mathutil::gettimeofday_sec();
     3800  double elapse1 = TimeEnd - TimeStart;
     3801  std::cout << "cspline-new   : " << elapse1 << " (sec.)" << endl;
     3802  /****/
     3803}
     3804
     3805void Scantable::autoCubicSplineBaseline(const std::vector<bool>& mask, int nPiece,
     3806                                        float thresClip, int nIterClip,
     3807                                        const std::vector<int>& edge,
     3808                                        float threshold, int chanAvgLimit,
     3809                                        bool getResidual,
     3810                                        const std::string& progressInfo,
     3811                                        const bool outLogger, const std::string& blfile,
     3812                                        const std::string& bltable)
    39123813{
    39133814  try {
    39143815    ofstream ofs;
    3915     String coordInfo = "";
    3916     bool hasSameNchan = true;
    3917     bool outTextFile = false;
    3918     bool csvFormat = false;
    3919 
    3920     if (blfile != "") {
    3921       csvFormat = (blfile.substr(0, 1) == "T");
    3922       ofs.open(blfile.substr(1).c_str(), ios::out | ios::app);
    3923       if (ofs) outTextFile = true;
    3924     }
    3925 
    3926     if (outLogger || outTextFile) {
    3927       coordInfo = getCoordInfo()[0];
    3928       if (coordInfo == "") coordInfo = "channel";
    3929       hasSameNchan = hasSameNchanOverIFs();
    3930     }
    3931 
    3932     bool showProgress;
     3816    String coordInfo;
     3817    bool hasSameNchan, outTextFile, csvFormat, showProgress;
    39333818    int minNRow;
    3934     parseProgressInfo(progressInfo, showProgress, minNRow);
    3935 
    3936     int minEdgeSize = getIFNos().size()*2;
     3819    int nRow = nrow();
     3820    std::vector<bool> chanMask, finalChanMask;
     3821    float rms;
     3822    bool outBaselineTable = (bltable != "");
     3823    STBaselineTable bt = STBaselineTable(*this);
     3824    Vector<Double> timeSecCol;
    39373825    STLineFinder lineFinder = STLineFinder();
    3938     lineFinder.setOptions(threshold, 3, chanAvgLimit);
    3939 
    3940     int nRow = nrow();
    3941     std::vector<bool> chanMask;
    3942     std::vector<bool> finalChanMask;
    3943     float rms;
    3944 
    3945     //---
    3946     bool outBaselintTable = (bltable != "");
    3947     STBaselineTable* bt = new STBaselineTable(*this);
    3948     ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(table_, "TIME");
    3949     Vector<Double> timeSecCol = tcol->getColumn();
    3950     //---
     3826
     3827    initialiseBaselining(blfile, ofs, outLogger, outTextFile, csvFormat,
     3828                         coordInfo, hasSameNchan,
     3829                         progressInfo, showProgress, minNRow,
     3830                         timeSecCol);
     3831
     3832    initLineFinder(edge, threshold, chanAvgLimit, lineFinder);
     3833
     3834    std::vector<int> nChanNos;
     3835    std::vector<std::vector<std::vector<double> > > modelReservoir;
     3836    modelReservoir = getPolynomialModelReservoir(3,
     3837                                                 &Scantable::getNormalPolynomial,
     3838                                                 nChanNos);
    39513839
    39523840    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    39533841      std::vector<float> sp = getSpectrum(whichrow);
    3954       //-- composote given channel mask and flagtra --------
    3955       int edgeSize = edge.size();
    39563842      std::vector<int> currentEdge;
    3957       currentEdge.clear();
    3958       if (edgeSize >= 2) {
    3959         int idx = 0;
    3960         if (edgeSize > 2) {
    3961           if (edgeSize < minEdgeSize) {
    3962             throw(AipsError("Length of edge element info is less than that of IFs"));
    3963           }
    3964           idx = 2 * getIF(whichrow);
    3965         }
    3966         currentEdge.push_back(edge[idx]);
    3967         currentEdge.push_back(edge[idx+1]);
    3968       } else {
    3969         throw(AipsError("Wrong length of edge element"));
    3970       }
    3971       lineFinder.setData(sp);
    3972       lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
    3973       chanMask = lineFinder.getMask();
    3974       //-- composote given channel mask and flagtra END --------
    3975 
    3976       std::vector<int> pieceEdges;//(nPiece+1);
    3977       std::vector<float> params;//(nPiece*4);
     3843      chanMask = getCompositeChanMask(whichrow, mask, edge, currentEdge, lineFinder);
     3844
     3845      std::vector<int> pieceEdges;
     3846      std::vector<float> params;
    39783847      int nClipped = 0;
    3979       std::vector<float> res = doCubicSplineFitting(sp, chanMask, nPiece, false, pieceEdges, params, rms, finalChanMask, nClipped, thresClip, nIterClip, getResidual);
    3980 
    3981       if (outBaselintTable) {
    3982         Vector<Int> fparam(nPiece);
    3983         for (uInt i = 0; i < fparam.size(); ++i) {
    3984           fparam[i] = pieceEdges[i];
    3985         }
    3986         std::vector<int> finalMaskList0;
    3987         finalMaskList0.clear();
    3988         for (uInt i = 0; i < finalChanMask.size(); ++i) {
    3989           if (finalChanMask[i]) {
    3990             if ((i == 0)||(i == finalChanMask.size()-1)) {
    3991               finalMaskList0.push_back(i);
    3992             } else {
    3993               if ((finalChanMask[i])&&(!finalChanMask[i-1])) {
    3994                 finalMaskList0.push_back(i);
    3995               }
    3996               if ((finalChanMask[i])&&(!finalChanMask[i+1])) {
    3997                 finalMaskList0.push_back(i);
    3998               }
    3999             }
    4000           }
    4001         }
    4002         Vector<uInt> finalMaskList(finalMaskList0.size());
    4003         for (uInt i = 0; i < finalMaskList0.size(); ++i) {
    4004           finalMaskList[i] = (uInt)finalMaskList0[i];
    4005         }
    4006 
    4007         Vector<uInt> cEdge(currentEdge.size());
    4008         for (uInt i = 0; i < currentEdge.size(); ++i) {
    4009           cEdge[i] = currentEdge[i];
    4010         }
    4011 
    4012         bt->appenddata(uInt(getScan(whichrow)), uInt(getCycle(whichrow)), uInt(getBeam(whichrow)), uInt(getIF(whichrow)), uInt(getPol(whichrow)),
    4013                        uInt(0),
    4014                        timeSecCol[whichrow],
    4015                        Bool(true),
    4016                        STBaselineFunc::CSpline,
    4017                        fparam,
    4018                        Vector<Float>(),
    4019                        finalMaskList,
    4020                        Vector<Float>(params),
    4021                        Float(rms),
    4022                        uInt(sp.size()),
    4023                        Float(thresClip),
    4024                        uInt(nIterClip),
    4025                        Float(threshold),
    4026                        uInt(chanAvgLimit),
    4027                        cEdge);
     3848      std::vector<float> res = doCubicSplineLeastSquareFitting(sp, chanMask,
     3849                                   modelReservoir[getIdxOfNchan(sp.size(), nChanNos)],
     3850                                   nPiece, false, pieceEdges, params, rms, finalChanMask,
     3851                                   nClipped, thresClip, nIterClip, getResidual);
     3852
     3853      if (outBaselineTable) {
     3854        bt.appenddata(getScan(whichrow), getCycle(whichrow), getBeam(whichrow),
     3855                      getIF(whichrow), getPol(whichrow), 0, timeSecCol[whichrow],
     3856                      true, STBaselineFunc::CSpline, pieceEdges, std::vector<float>(),
     3857                      getMaskListFromMask(finalChanMask), params, rms, sp.size(),
     3858                      thresClip, nIterClip, threshold, chanAvgLimit, currentEdge);
    40283859      } else {
    40293860        setSpectrum(res, whichrow);
    40303861      }
    40313862
    4032       outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()", pieceEdges, params, nClipped);
     3863      outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow,
     3864                          coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()",
     3865                          pieceEdges, params, nClipped);
    40333866      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
    40343867    }
    40353868
    4036     //---
    4037     if (outBaselintTable) {
    4038       bt->save(bltable);
    4039     }
    4040     //---
    4041     delete bt;
    4042 
    4043     if (outTextFile) ofs.close();
     3869    finaliseBaselining(outBaselineTable, &bt, bltable, outTextFile, ofs);
    40443870
    40453871  } catch (...) {
    … …  
    40483874}
    40493875
    4050 std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data, const std::vector<bool>& mask, std::vector<int>& idxEdge, std::vector<float>& params, float& rms, std::vector<bool>& finalmask, float clipth, int clipn)
     3876std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data,
     3877                                                   const std::vector<bool>& mask,
     3878                                                   std::vector<int>& idxEdge,
     3879                                                   std::vector<float>& params,
     3880                                                   float& rms,
     3881                                                   std::vector<bool>& finalmask,
     3882                                                   float clipth,
     3883                                                   int clipn)
    40513884{
    40523885  int nClipped = 0;
    … …  
    40543887}
    40553888
    4056 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& rms, std::vector<bool>& finalmask, float clipth, int clipn)
     3889std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data,
     3890                                                   const std::vector<bool>& mask,
     3891                                                   int nPiece,
     3892                                                   std::vector<int>& idxEdge,
     3893                                                   std::vector<float>& params,
     3894                                                   float& rms,
     3895                                                   std::vector<bool>& finalmask,
     3896                                                   float clipth,
     3897                                                   int clipn)
    40573898{
    40583899  int nClipped = 0;
    … …  
    40603901}
    40613902
    4062 std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data, const std::vector<bool>& mask, int nPiece, bool useGivenPieceBoundary, std::vector<int>& idxEdge, std::vector<float>& params, float& rms, std::vector<bool>& finalMask, int& nClipped, float thresClip, int nIterClip, bool getResidual)
    4063 {
    4064   if (data.size() != mask.size()) {
    4065     throw(AipsError("data and mask sizes are not identical"));
     3903std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data,
     3904                                                   const std::vector<bool>& mask,
     3905                                                   int nPiece,
     3906                                                   bool useGivenPieceBoundary,
     3907                                                   std::vector<int>& idxEdge,
     3908                                                   std::vector<float>& params,
     3909                                                   float& rms,
     3910                                                   std::vector<bool>& finalMask,
     3911                                                   int& nClipped,
     3912                                                   float thresClip,
     3913                                                   int nIterClip,
     3914                                                   bool getResidual)
     3915{
     3916  return doCubicSplineLeastSquareFitting(data, mask,
     3917                                         getPolynomialModel(3, data.size(), &Scantable::getNormalPolynomial),
     3918                                         nPiece, useGivenPieceBoundary, idxEdge,
     3919                                         params, rms, finalMask,
     3920                                         nClipped, thresClip, nIterClip,
     3921                                         getResidual);
     3922}
     3923
     3924std::vector<float> Scantable::doCubicSplineLeastSquareFitting(const std::vector<float>& data,
     3925                                                              const std::vector<bool>& mask,
     3926                                                              const std::vector<std::vector<double> >& model,
     3927                                                              int nPiece,
     3928                                                              bool useGivenPieceBoundary,
     3929                                                              std::vector<int>& idxEdge,
     3930                                                              std::vector<float>& params,
     3931                                                              float& rms,
     3932                                                              std::vector<bool>& finalMask,
     3933                                                              int& nClipped,
     3934                                                              float thresClip,
     3935                                                              int nIterClip,
     3936                                                              bool getResidual)
     3937{
     3938  int nDOF = nPiece + 3;  //number of independent parameters to solve, namely, 4+(nPiece-1).
     3939  int nModel = model.size();
     3940  int nChan = data.size();
     3941
     3942  if (nModel != 4) {
     3943    throw(AipsError("model size must be 4."));
    40663944  }
    40673945  if (nPiece < 1) {
    40683946    throw(AipsError("number of the sections must be one or more"));
    40693947  }
    4070 
    4071   int nChan = data.size();
     3948  if (nChan < 2*nPiece) {
     3949    throw(AipsError("data size is too few"));
     3950  }
     3951  if (nChan != (int)mask.size()) {
     3952    throw(AipsError("data and mask sizes are not identical"));
     3953  }
     3954  for (int i = 0; i < nModel; ++i) {
     3955    if (nChan != (int)model[i].size()) {
     3956      throw(AipsError("data and model sizes are not identical"));
     3957    }
     3958  }
     3959
     3960  params.clear();
     3961  params.resize(nPiece*nModel);
    40723962
    40733963  finalMask.clear();
    … …  
    40853975    finalMask[i] = mask[i];
    40863976  }
     3977
    40873978  int initNData = j;
     3979  int nData = initNData;
    40883980
    40893981  if (initNData < nPiece) {
    … …  
    41144006  }
    41154007
    4116   params.clear();
    4117   params.resize(nPiece*4);
    4118 
    4119   int nData = initNData;
    4120   int nDOF = nPiece + 3;  //number of parameters to solve, namely, 4+(nPiece-1).
    4121 
    4122   std::vector<double> x1(nChan), x2(nChan), x3(nChan);
    4123   std::vector<double> z1(nChan), x1z1(nChan), x2z1(nChan), x3z1(nChan);
    4124   std::vector<double> r1(nChan), residual(nChan);
     4008  std::vector<double> z1(nChan), r1(nChan), residual(nChan);
    41254009  for (int i = 0; i < nChan; ++i) {
    4126     double di = (double)i;
    4127     double dD = (double)data[i];
    4128     x1[i]   = di;
    4129     x2[i]   = di*di;
    4130     x3[i]   = di*di*di;
    4131     z1[i]   = dD;
    4132     x1z1[i] = dD*di;
    4133     x2z1[i] = dD*di*di;
    4134     x3z1[i] = dD*di*di*di;
    4135     r1[i]   = 0.0;
     4010    z1[i] = (double)data[i];
     4011    r1[i] = 0.0;
    41364012    residual[i] = 0.0;
    41374013  }
    … …  
    41554031    for (int n = 0; n < nPiece; ++n) {
    41564032      int nUseDataInPiece = 0;
    4157       for (int i = idxEdge[n]; i < idxEdge[n+1]; ++i) {
    4158 
    4159         if (maskArray[i] == 0) continue;
    4160 
    4161         xMatrix[0][0] += 1.0;
    4162         xMatrix[0][1] += x1[i];
    4163         xMatrix[0][2] += x2[i];
    4164         xMatrix[0][3] += x3[i];
    4165         xMatrix[1][1] += x2[i];
    4166         xMatrix[1][2] += x3[i];
    4167         xMatrix[1][3] += x2[i]*x2[i];
    4168         xMatrix[2][2] += x2[i]*x2[i];
    4169         xMatrix[2][3] += x3[i]*x2[i];
    4170         xMatrix[3][3] += x3[i]*x3[i];
    4171         zMatrix[0] += z1[i];
    4172         zMatrix[1] += x1z1[i];
    4173         zMatrix[2] += x2z1[i];
    4174         zMatrix[3] += x3z1[i];
    4175 
    4176         for (int j = 0; j < n; ++j) {
    4177           double q = 1.0 - x1[i]*invEdge[j];
     4033      for (int k = idxEdge[n]; k < idxEdge[n+1]; ++k) {
     4034
     4035        if (maskArray[k] == 0) continue;
     4036
     4037        for (int i = 0; i < nModel; ++i) {
     4038          for (int j = i; j < nModel; ++j) {
     4039            xMatrix[i][j] += model[i][k] * model[j][k];
     4040          }
     4041          zMatrix[i] += z1[k] * model[i][k];
     4042        }
     4043
     4044        for (int i = 0; i < n; ++i) {
     4045          double q = 1.0 - model[1][k]*invEdge[i];
    41784046          q = q*q*q;
    4179           xMatrix[0][j+4] += q;
    4180           xMatrix[1][j+4] += q*x1[i];
    4181           xMatrix[2][j+4] += q*x2[i];
    4182           xMatrix[3][j+4] += q*x3[i];
    4183           for (int k = 0; k < j; ++k) {
    4184             double r = 1.0 - x1[i]*invEdge[k];
     4047          for (int j = 0; j < nModel; ++j) {
     4048            xMatrix[j][i+nModel] += q * model[j][k];
     4049          }
     4050          for (int j = 0; j < i; ++j) {
     4051            double r = 1.0 - model[1][k]*invEdge[j];
    41854052            r = r*r*r;
    4186             xMatrix[k+4][j+4] += r*q;
     4053            xMatrix[j+nModel][i+nModel] += r*q;
    41874054          }
    4188           xMatrix[j+4][j+4] += q*q;
    4189           zMatrix[j+4] += q*z1[i];
     4055          xMatrix[i+nModel][i+nModel] += q*q;
     4056          zMatrix[i+nModel] += q*z1[k];
    41904057        }
    41914058
    … …  
    42154082    std::vector<double> invDiag(nDOF);
    42164083    for (int i = 0; i < nDOF; ++i) {
    4217       invDiag[i] = 1.0/xMatrix[i][i];
     4084      invDiag[i] = 1.0 / xMatrix[i][i];
    42184085      for (int j = 0; j < nDOF; ++j) {
    42194086        xMatrix[i][j] *= invDiag[i];
    … …  
    42254092        if (i != k) {
    42264093          double factor1 = xMatrix[k][k];
     4094          double invfactor1 = 1.0 / factor1;
    42274095          double factor2 = xMatrix[i][k];
    42284096          for (int j = k; j < 2*nDOF; ++j) {
    42294097            xMatrix[i][j] *= factor1;
    42304098            xMatrix[i][j] -= xMatrix[k][j]*factor2;
    4231             xMatrix[i][j] /= factor1;
     4099            xMatrix[i][j] *= invfactor1;
    42324100          }
    42334101        }
    42344102      }
    4235       double xDiag = xMatrix[k][k];
     4103      double invXDiag = 1.0 / xMatrix[k][k];
    42364104      for (int j = k; j < 2*nDOF; ++j) {
    4237         xMatrix[k][j] /= xDiag;
     4105        xMatrix[k][j] *= invXDiag;
    42384106      }
    42394107    }
    … …  
    42564124    }
    42574125
    4258     double a0 = y[0];
    4259     double a1 = y[1];
    4260     double a2 = y[2];
    4261     double a3 = y[3];
     4126    std::vector<double> a(nModel);
     4127    for (int i = 0; i < nModel; ++i) {
     4128      a[i] = y[i];
     4129    }
    42624130
    42634131    int j = 0;
    42644132    for (int n = 0; n < nPiece; ++n) {
    42654133      for (int i = idxEdge[n]; i < idxEdge[n+1]; ++i) {
    4266         r1[i] = a0 + a1*x1[i] + a2*x2[i] + a3*x3[i];
    4267       }
    4268       params[j]   = a0;
    4269       params[j+1] = a1;
    4270       params[j+2] = a2;
    4271       params[j+3] = a3;
    4272       j += 4;
     4134        r1[i] = 0.0;
     4135        for (int j = 0; j < nModel; ++j) {
     4136          r1[i] += a[j] * model[j][i];
     4137        }
     4138      }
     4139      for (int i = 0; i < nModel; ++i) {
     4140        params[j+i] = a[i];
     4141      }
     4142      j += nModel;
    42734143
    42744144      if (n == nPiece-1) break;
    42754145
    4276       double d = y[4+n];
     4146      double d = y[n+nModel];
    42774147      double iE = invEdge[n];
    4278       a0 +=     d;
    4279       a1 -= 3.0*d*iE;
    4280       a2 += 3.0*d*iE*iE;
    4281       a3 -=     d*iE*iE*iE;
     4148      a[0] +=       d;
     4149      a[1] -= 3.0 * d * iE;
     4150      a[2] += 3.0 * d * iE * iE;
     4151      a[3] -=       d * iE * iE * iE;
    42824152    }
    42834153
    … …  
    43594229}
    43604230
    4361 void Scantable::selectWaveNumbers(const int whichrow, const std::vector<bool>& chanMask, const std::string& fftInfo, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, std::vector<int>& nWaves)
    4362 //void Scantable::selectWaveNumbers(const int whichrow, const std::vector<bool>& chanMask, const bool applyFFT, const std::string& fftMethod, const std::string& fftThresh, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, std::vector<int>& nWaves)
    4363 {
    4364   bool applyFFT;
     4231std::vector<int> Scantable::selectWaveNumbers(const std::vector<int>& addNWaves,
     4232                                  const std::vector<int>& rejectNWaves)
     4233{
     4234  std::vector<bool> chanMask;
    43654235  std::string fftMethod;
    43664236  std::string fftThresh;
    43674237
     4238  return selectWaveNumbers(0, chanMask, false, fftMethod, fftThresh, addNWaves, rejectNWaves);
     4239}
     4240
     4241std::vector<int> Scantable::selectWaveNumbers(const int whichrow,
     4242                                  const std::vector<bool>& chanMask,
     4243                                  const bool applyFFT,
     4244                                  const std::string& fftMethod,
     4245                                  const std::string& fftThresh,
     4246                                  const std::vector<int>& addNWaves,
     4247                                  const std::vector<int>& rejectNWaves)
     4248{
     4249  std::vector<int> nWaves;
    43684250  nWaves.clear();
    43694251
    4370   parseFFTInfo(fftInfo, applyFFT, fftMethod, fftThresh);
    43714252  if (applyFFT) {
    43724253    string fftThAttr;
    43734254    float fftThSigma;
    43744255    int fftThTop;
    4375     parseThresholdExpression(fftThresh, fftThAttr, fftThSigma, fftThTop);
     4256    parseFFTThresholdInfo(fftThresh, fftThAttr, fftThSigma, fftThTop);
    43764257    doSelectWaveNumbers(whichrow, chanMask, fftMethod, fftThSigma, fftThTop, fftThAttr, nWaves);
    43774258  }
    43784259
    43794260  addAuxWaveNumbers(whichrow, addNWaves, rejectNWaves, nWaves);
     4261
     4262  return nWaves;
     4263}
     4264
     4265int Scantable::getIdxOfNchan(const int nChan, const std::vector<int>& nChanNos)
     4266{
     4267  int idx = -1;
     4268  for (uint i = 0; i < nChanNos.size(); ++i) {
     4269    if (nChan == nChanNos[i]) {
     4270      idx = i;
     4271      break;
     4272    }
     4273  }
     4274
     4275  if (idx < 0) {
     4276    throw(AipsError("nChan not found in nChhanNos."));
     4277  }
     4278
     4279  return idx;
    43804280}
    43814281
    … …  
    43964296}
    43974297
    4398 void Scantable::parseThresholdExpression(const std::string& fftThresh, std::string& fftThAttr, float& fftThSigma, int& fftThTop)
     4298void Scantable::parseFFTThresholdInfo(const std::string& fftThresh, std::string& fftThAttr, float& fftThSigma, int& fftThTop)
    43994299{
    44004300  uInt idxSigma = fftThresh.find("sigma");
    … …  
    45364436}
    45374437
    4538 void Scantable::sinusoidBaseline(const std::vector<bool>& mask, const std::string& fftInfo, 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, const std::string& bltable)
    4539 {
     4438void Scantable::sinusoidBaseline(const std::vector<bool>& mask, const std::string& fftInfo,
     4439                                 const std::vector<int>& addNWaves,
     4440                                 const std::vector<int>& rejectNWaves,
     4441                                 float thresClip, int nIterClip,
     4442                                 bool getResidual,
     4443                                 const std::string& progressInfo,
     4444                                 const bool outLogger, const std::string& blfile,
     4445                                 const std::string& bltable)
     4446{
     4447  /****/
     4448  double TimeStart = mathutil::gettimeofday_sec();
     4449  /****/
     4450
    45404451  try {
    45414452    ofstream ofs;
    4542     String coordInfo = "";
    4543     bool hasSameNchan = true;
    4544     bool outTextFile = false;
    4545     bool csvFormat = false;
    4546 
    4547     if (blfile != "") {
    4548       csvFormat = (blfile.substr(0, 1) == "T");
    4549       ofs.open(blfile.substr(1).c_str(), ios::out | ios::app);
    4550       if (ofs) outTextFile = true;
    4551     }
    4552 
    4553     if (outLogger || outTextFile) {
    4554       coordInfo = getCoordInfo()[0];
    4555       if (coordInfo == "") coordInfo = "channel";
    4556       hasSameNchan = hasSameNchanOverIFs();
    4557     }
    4558 
    4559     bool showProgress;
     4453    String coordInfo;
     4454    bool hasSameNchan, outTextFile, csvFormat, showProgress;
    45604455    int minNRow;
    4561     parseProgressInfo(progressInfo, showProgress, minNRow);
    4562 
    45634456    int nRow = nrow();
    4564     std::vector<bool> chanMask;
     4457    std::vector<bool> chanMask, finalChanMask;
     4458    float rms;
     4459    bool outBaselineTable = (bltable != "");
     4460    STBaselineTable bt = STBaselineTable(*this);
     4461    Vector<Double> timeSecCol;
     4462
     4463    initialiseBaselining(blfile, ofs, outLogger, outTextFile, csvFormat,
     4464                         coordInfo, hasSameNchan,
     4465                         progressInfo, showProgress, minNRow,
     4466                         timeSecCol);
     4467
     4468    bool applyFFT;
     4469    std::string fftMethod, fftThresh;
     4470    parseFFTInfo(fftInfo, applyFFT, fftMethod, fftThresh);
     4471
    45654472    std::vector<int> nWaves;
    4566     std::vector<bool> finalChanMask;
    4567     float rms;
    4568 
    4569     //---
    4570     bool outBaselintTable = (bltable != "");
    4571     STBaselineTable* bt = new STBaselineTable(*this);
    4572     ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(table_, "TIME");
    4573     Vector<Double> timeSecCol = tcol->getColumn();
    4574     //---
     4473    std::vector<int> nChanNos;
     4474    std::vector<std::vector<std::vector<double> > > modelReservoir;
     4475    if (!applyFFT) {
     4476      nWaves = selectWaveNumbers(addNWaves, rejectNWaves);
     4477      modelReservoir = getSinusoidModelReservoir(nWaves, nChanNos);
     4478    }
    45754479
    45764480    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    45774481      std::vector<float> sp = getSpectrum(whichrow);
    45784482      chanMask = getCompositeChanMask(whichrow, mask);
    4579       selectWaveNumbers(whichrow, chanMask, fftInfo, addNWaves, rejectNWaves, nWaves);
     4483      std::vector<std::vector<double> > model;
     4484      if (applyFFT) {
     4485        nWaves = selectWaveNumbers(whichrow, chanMask, true, fftMethod, fftThresh,
     4486                                   addNWaves, rejectNWaves);
     4487        model = getSinusoidModel(nWaves, sp.size());
     4488      } else {
     4489        model = modelReservoir[getIdxOfNchan(sp.size(), nChanNos)];
     4490      }
    45804491
    45814492      std::vector<float> params;
    45824493      int nClipped = 0;
    4583       std::vector<float> res = doSinusoidFitting(sp, chanMask, nWaves, params, rms, finalChanMask, nClipped, thresClip, nIterClip, getResidual);
    4584 
    4585       //---
    4586       if (outBaselintTable) {
    4587         uInt nparam = nWaves.size();
    4588         Vector<Int> fparam(nparam);
    4589         for (uInt i = 0; i < nparam; ++i) {
    4590           fparam[i] = nWaves[i];
    4591         }
    4592         std::vector<int> finalMaskList0;
    4593         finalMaskList0.clear();
    4594         for (uInt i = 0; i < finalChanMask.size(); ++i) {
    4595           if (finalChanMask[i]) {
    4596             if ((i == 0)||(i == finalChanMask.size()-1)) {
    4597               finalMaskList0.push_back(i);
    4598             } else {
    4599               if ((finalChanMask[i])&&(!finalChanMask[i-1])) {
    4600                 finalMaskList0.push_back(i);
    4601               }
    4602               if ((finalChanMask[i])&&(!finalChanMask[i+1])) {
    4603                 finalMaskList0.push_back(i);
    4604               }
    4605             }
    4606           }
    4607         }
    4608         Vector<uInt> finalMaskList(finalMaskList0.size());
    4609         for (uInt i = 0; i < finalMaskList0.size(); ++i) {
    4610           finalMaskList[i] = (uInt)finalMaskList0[i];
    4611         }
    4612 
    4613         bt->appenddata(uInt(getScan(whichrow)), uInt(getCycle(whichrow)), uInt(getBeam(whichrow)), uInt(getIF(whichrow)), uInt(getPol(whichrow)),
    4614                        uInt(0),
    4615                        timeSecCol[whichrow],
    4616                        Bool(true),
    4617                        STBaselineFunc::Sinusoid,
    4618                        fparam,
    4619                        Vector<Float>(),
    4620                        finalMaskList,
    4621                        Vector<Float>(params),
    4622                        Float(rms),
    4623                        uInt(sp.size()),
    4624                        Float(thresClip),
    4625                        uInt(nIterClip),
    4626                        Float(0.0),
    4627                        uInt(0),
    4628                        Vector<uInt>());
     4494      std::vector<float> res = doLeastSquareFitting(sp, chanMask, model,
     4495                                   params, rms, finalChanMask,
     4496                                   nClipped, thresClip, nIterClip, getResidual);
     4497
     4498      if (outBaselineTable) {
     4499        bt.appenddata(getScan(whichrow), getCycle(whichrow), getBeam(whichrow),
     4500                      getIF(whichrow), getPol(whichrow), 0, timeSecCol[whichrow],
     4501                      true, STBaselineFunc::Sinusoid, nWaves, std::vector<float>(),
     4502                      getMaskListFromMask(finalChanMask), params, rms, sp.size(),
     4503                      thresClip, nIterClip, 0.0, 0, std::vector<int>());
    46294504      } else {
    46304505        setSpectrum(res, whichrow);
    46314506      }
    4632       //---
    4633 
    4634       outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "sinusoidBaseline()", params, nClipped);
     4507
     4508      outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow,
     4509                          coordInfo, hasSameNchan, ofs, "sinusoidBaseline()",
     4510                          params, nClipped);
    46354511      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
    46364512    }
    46374513
    4638     //---
    4639     if (outBaselintTable) {
    4640       bt->save(bltable);
    4641     }
    4642     //---
    4643     delete bt;
    4644 
    4645     if (outTextFile) ofs.close();
     4514    finaliseBaselining(outBaselineTable, &bt, bltable, outTextFile, ofs);
    46464515
    46474516  } catch (...) {
    46484517    throw;
    46494518  }
    4650 }
    4651 
    4652 void Scantable::autoSinusoidBaseline(const std::vector<bool>& mask, const std::string& fftInfo, 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, const std::string& bltable)
    4653 //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, const std::string& bltable)
     4519
     4520  /****/
     4521  double TimeEnd = mathutil::gettimeofday_sec();
     4522  double elapse1 = TimeEnd - TimeStart;
     4523  std::cout << "sinusoid-old   : " << elapse1 << " (sec.)" << endl;
     4524  /****/
     4525}
     4526
     4527void Scantable::autoSinusoidBaseline(const std::vector<bool>& mask, const std::string& fftInfo,
     4528                                     const std::vector<int>& addNWaves,
     4529                                     const std::vector<int>& rejectNWaves,
     4530                                     float thresClip, int nIterClip,
     4531                                     const std::vector<int>& edge,
     4532                                     float threshold, int chanAvgLimit,
     4533                                     bool getResidual,
     4534                                     const std::string& progressInfo,
     4535                                     const bool outLogger, const std::string& blfile,
     4536                                     const std::string& bltable)
    46544537{
    46554538  try {
    46564539    ofstream ofs;
    4657     String coordInfo = "";
    4658     bool hasSameNchan = true;
    4659     bool outTextFile = false;
    4660     bool csvFormat = false;
    4661 
    4662     if (blfile != "") {
    4663       csvFormat = (blfile.substr(0, 1) == "T");
    4664       ofs.open(blfile.substr(1).c_str(), ios::out | ios::app);
    4665       if (ofs) outTextFile = true;
    4666     }
    4667 
    4668     if (outLogger || outTextFile) {
    4669       coordInfo = getCoordInfo()[0];
    4670       if (coordInfo == "") coordInfo = "channel";
    4671       hasSameNchan = hasSameNchanOverIFs();
    4672     }
    4673 
    4674     bool showProgress;
     4540    String coordInfo;
     4541    bool hasSameNchan, outTextFile, csvFormat, showProgress;
    46754542    int minNRow;
    4676     parseProgressInfo(progressInfo, showProgress, minNRow);
    4677 
    4678     int minEdgeSize = getIFNos().size()*2;
     4543    int nRow = nrow();
     4544    std::vector<bool> chanMask, finalChanMask;
     4545    float rms;
     4546    bool outBaselineTable = (bltable != "");
     4547    STBaselineTable bt = STBaselineTable(*this);
     4548    Vector<Double> timeSecCol;
    46794549    STLineFinder lineFinder = STLineFinder();
    4680     lineFinder.setOptions(threshold, 3, chanAvgLimit);
    4681 
    4682     int nRow = nrow();
    4683     std::vector<bool> chanMask;
     4550
     4551    initialiseBaselining(blfile, ofs, outLogger, outTextFile, csvFormat,
     4552                         coordInfo, hasSameNchan,
     4553                         progressInfo, showProgress, minNRow,
     4554                         timeSecCol);
     4555
     4556    initLineFinder(edge, threshold, chanAvgLimit, lineFinder);
     4557
     4558    bool applyFFT;
     4559    string fftMethod, fftThresh;
     4560    parseFFTInfo(fftInfo, applyFFT, fftMethod, fftThresh);
     4561
    46844562    std::vector<int> nWaves;
    4685     std::vector<bool> finalChanMask;
    4686     float rms;
    4687 
    4688     //---
    4689     bool outBaselintTable = (bltable != "");
    4690     STBaselineTable* bt = new STBaselineTable(*this);
    4691     ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(table_, "TIME");
    4692     Vector<Double> timeSecCol = tcol->getColumn();
    4693     //---
     4563    std::vector<int> nChanNos;
     4564    std::vector<std::vector<std::vector<double> > > modelReservoir;
     4565    if (!applyFFT) {
     4566      nWaves = selectWaveNumbers(addNWaves, rejectNWaves);
     4567      modelReservoir = getSinusoidModelReservoir(nWaves, nChanNos);
     4568    }
    46944569
    46954570    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    46964571      std::vector<float> sp = getSpectrum(whichrow);
    4697       //-- composote given channel mask and flagtra --------
    4698       int edgeSize = edge.size();
    46994572      std::vector<int> currentEdge;
    4700       currentEdge.clear();
    4701       if (edgeSize >= 2) {
    4702         int idx = 0;
    4703         if (edgeSize > 2) {
    4704           if (edgeSize < minEdgeSize) {
    4705             throw(AipsError("Length of edge element info is less than that of IFs"));
    4706           }
    4707           idx = 2 * getIF(whichrow);
    4708         }
    4709         currentEdge.push_back(edge[idx]);
    4710         currentEdge.push_back(edge[idx+1]);
     4573      chanMask = getCompositeChanMask(whichrow, mask, edge, currentEdge, lineFinder);
     4574      std::vector<std::vector<double> > model;
     4575      if (applyFFT) {
     4576        nWaves = selectWaveNumbers(whichrow, chanMask, true, fftMethod, fftThresh,
     4577                                   addNWaves, rejectNWaves);
     4578        model = getSinusoidModel(nWaves, sp.size());
    47114579      } else {
    4712         throw(AipsError("Wrong length of edge element"));
    4713       }
    4714       lineFinder.setData(sp);
    4715       lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
    4716       chanMask = lineFinder.getMask();
    4717       //-- composote given channel mask and flagtra END --------
    4718       selectWaveNumbers(whichrow, chanMask, fftInfo, addNWaves, rejectNWaves, nWaves);
     4580        model = modelReservoir[getIdxOfNchan(sp.size(), nChanNos)];
     4581      }
    47194582
    47204583      std::vector<float> params;
    47214584      int nClipped = 0;
    4722       std::vector<float> res = doSinusoidFitting(sp, chanMask, nWaves, params, rms, finalChanMask, nClipped, thresClip, nIterClip, getResidual);
    4723 
    4724       //---
    4725       if (outBaselintTable) {
    4726         uInt nparam = nWaves.size();
    4727         Vector<Int> fparam(nparam);
    4728         for (uInt i = 0; i < nparam; ++i) {
    4729           fparam[i] = nWaves[i];
    4730         }
    4731         std::vector<int> finalMaskList0;
    4732         finalMaskList0.clear();
    4733         for (uInt i = 0; i < finalChanMask.size(); ++i) {
    4734           if (finalChanMask[i]) {
    4735             if ((i == 0)||(i == finalChanMask.size()-1)) {
    4736               finalMaskList0.push_back(i);
    4737             } else {
    4738               if ((finalChanMask[i])&&(!finalChanMask[i-1])) {
    4739                 finalMaskList0.push_back(i);
    4740               }
    4741               if ((finalChanMask[i])&&(!finalChanMask[i+1])) {
    4742                 finalMaskList0.push_back(i);
    4743               }
    4744             }
    4745           }
    4746         }
    4747         Vector<uInt> finalMaskList(finalMaskList0.size());
    4748         for (uInt i = 0; i < finalMaskList0.size(); ++i) {
    4749           finalMaskList[i] = (uInt)finalMaskList0[i];
    4750         }
    4751 
    4752         Vector<uInt> cEdge(currentEdge.size());
    4753         for (uInt i = 0; i < currentEdge.size(); ++i) {
    4754           cEdge[i] = currentEdge[i];
    4755         }
    4756 
    4757         bt->appenddata(uInt(getScan(whichrow)), uInt(getCycle(whichrow)), uInt(getBeam(whichrow)), uInt(getIF(whichrow)), uInt(getPol(whichrow)),
    4758                        uInt(0),
    4759                        timeSecCol[whichrow],
    4760                        Bool(true),
    4761                        STBaselineFunc::Sinusoid,
    4762                        fparam,
    4763                        Vector<Float>(),
    4764                        finalMaskList,
    4765                        Vector<Float>(params),
    4766                        Float(rms),
    4767                        uInt(sp.size()),
    4768                        Float(thresClip),
    4769                        uInt(nIterClip),
    4770                        Float(threshold),
    4771                        uInt(chanAvgLimit),
    4772                        cEdge);
     4585      std::vector<float> res = doLeastSquareFitting(sp, chanMask, model,
     4586                                   params, rms, finalChanMask,
     4587                                   nClipped, thresClip, nIterClip, getResidual);
     4588
     4589      if (outBaselineTable) {
     4590        bt.appenddata(getScan(whichrow), getCycle(whichrow), getBeam(whichrow),
     4591                      getIF(whichrow), getPol(whichrow), 0, timeSecCol[whichrow],
     4592                      true, STBaselineFunc::Sinusoid, nWaves, std::vector<float>(),
     4593                      getMaskListFromMask(finalChanMask), params, rms, sp.size(),
     4594                      thresClip, nIterClip, threshold, chanAvgLimit, currentEdge);
    47734595      } else {
    47744596        setSpectrum(res, whichrow);
    47754597      }
    4776       //---
    4777 
    4778       outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoSinusoidBaseline()", params, nClipped);
     4598
     4599      outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow,
     4600                          coordInfo, hasSameNchan, ofs, "autoSinusoidBaseline()",
     4601                          params, nClipped);
    47794602      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
    47804603    }
    47814604
    4782     //---
    4783     if (outBaselintTable) {
    4784       bt->save(bltable);
    4785     }
    4786     //---
    4787     delete bt;
    4788 
    4789     if (outTextFile) ofs.close();
     4605    finaliseBaselining(outBaselineTable, &bt, bltable, outTextFile, ofs);
    47904606
    47914607  } catch (...) {
    … …  
    47944610}
    47954611
    4796 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& rms, std::vector<bool>& finalmask, float clipth, int clipn)
     4612std::vector<float> Scantable::doSinusoidFitting(const std::vector<float>& data,
     4613                                                const std::vector<bool>& mask,
     4614                                                const std::vector<int>& waveNumbers,
     4615                                                std::vector<float>& params,
     4616                                                float& rms,
     4617                                                std::vector<bool>& finalmask,
     4618                                                float clipth,
     4619                                                int clipn)
    47974620{
    47984621  int nClipped = 0;
    … …  
    48004623}
    48014624
    4802 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& rms, std::vector<bool>& finalMask, int& nClipped, float thresClip, int nIterClip, bool getResidual)
    4803 {
    4804   if (data.size() != mask.size()) {
    4805     throw(AipsError("data and mask sizes are not identical"));
    4806   }
    4807   if (data.size() < 2) {
    4808     throw(AipsError("data size is too short"));
    4809   }
    4810   if (waveNumbers.size() == 0) {
    4811     throw(AipsError("no wave numbers given"));
    4812   }
     4625std::vector<float> Scantable::doSinusoidFitting(const std::vector<float>& data,
     4626                                                const std::vector<bool>& mask,
     4627                                                const std::vector<int>& waveNumbers,
     4628                                                std::vector<float>& params,
     4629                                                float& rms,
     4630                                                std::vector<bool>& finalMask,
     4631                                                int& nClipped,
     4632                                                float thresClip,
     4633                                                int nIterClip,
     4634                                                bool getResidual)
     4635{
     4636  return doLeastSquareFitting(data, mask,
     4637                              getSinusoidModel(waveNumbers, data.size()),
     4638                              params, rms, finalMask,
     4639                              nClipped, thresClip, nIterClip,
     4640                              getResidual);
     4641}
     4642
     4643std::vector<std::vector<std::vector<double> > > Scantable::getSinusoidModelReservoir(const std::vector<int>& waveNumbers,
     4644                                                                                     std::vector<int>& nChanNos)
     4645{
     4646  std::vector<std::vector<std::vector<double> > > res;
     4647  res.clear();
     4648  nChanNos.clear();
     4649
     4650  std::vector<uint> ifNos = getIFNos();
     4651  for (uint i = 0; i < ifNos.size(); ++i) {
     4652    int currNchan = nchan(ifNos[i]);
     4653    bool hasDifferentNchan = (i == 0);
     4654    for (uint j = 0; j < i; ++j) {
     4655      if (currNchan != nchan(ifNos[j])) {
     4656        hasDifferentNchan = true;
     4657        break;
     4658      }
     4659    }
     4660    if (hasDifferentNchan) {
     4661      res.push_back(getSinusoidModel(waveNumbers, currNchan));
     4662      nChanNos.push_back(currNchan);
     4663    }
     4664  }
     4665
     4666  return res;
     4667}
     4668
     4669std::vector<std::vector<double> > Scantable::getSinusoidModel(const std::vector<int>& waveNumbers, int nchan)
     4670{
     4671  // model  : contains elemental values for computing the least-square matrix.
     4672  //          model.size() is nmodel and model[*].size() is nchan.
     4673  //          Each model element are as follows:
     4674  //          model[0]    = {1.0, 1.0, 1.0, ..., 1.0},
     4675  //          model[2n-1] = {sin(nPI/L*x[0]), sin(nPI/L*x[1]), ..., sin(nPI/L*x[nchan])},
     4676  //          model[2n]   = {cos(nPI/L*x[0]), cos(nPI/L*x[1]), ..., cos(nPI/L*x[nchan])},
     4677  //          where (1 <= n <= nMaxWavesInSW),
     4678  //          or,
     4679  //          model[2n-1] = {sin(wn[n]PI/L*x[0]), sin(wn[n]PI/L*x[1]), ..., sin(wn[n]PI/L*x[nchan])},
     4680  //          model[2n]   = {cos(wn[n]PI/L*x[0]), cos(wn[n]PI/L*x[1]), ..., cos(wn[n]PI/L*x[nchan])},
     4681  //          where wn[n] denotes waveNumbers[n] (1 <= n <= waveNumbers.size()).
     4682
    48134683  std::vector<int> nWaves;  // sorted and uniqued array of wave numbers
    48144684  nWaves.reserve(waveNumbers.size());
    … …  
    48234693  }
    48244694  bool hasConstantTerm = (minNWaves == 0);
    4825 
    4826   int nChan = data.size();
    4827 
    4828   finalMask.clear();
    4829   finalMask.resize(nChan);
    4830 
    4831   std::vector<int> maskArray;
    4832   std::vector<int> x;
    4833   for (int i = 0; i < nChan; ++i) {
    4834     maskArray.push_back(mask[i] ? 1 : 0);
    4835     if (mask[i]) {
    4836       x.push_back(i);
    4837     }
    4838     finalMask[i] = mask[i];
    4839   }
    4840 
    4841   int initNData = x.size();
    4842 
    4843   int nData = initNData;
    4844   int nDOF = nWaves.size() * 2 - (hasConstantTerm ? 1 : 0);  //number of parameters to solve.
     4695  int nmodel = nWaves.size() * 2 - (hasConstantTerm ? 1 : 0);  //number of parameters to solve.
     4696
     4697  std::vector<std::vector<double> > model(nmodel, std::vector<double>(nchan));
     4698
     4699  if (hasConstantTerm) {
     4700    for (int j = 0; j < nchan; ++j) {
     4701      model[0][j] = 1.0;
     4702    }
     4703  }
    48454704
    48464705  const double PI = 6.0 * asin(0.5); // PI (= 3.141592653...)
    4847   double baseXFactor = 2.0*PI/(double)(nChan-1);  //the denominator (nChan-1) should be changed to (xdata[nChan-1]-xdata[0]) for accepting x-values given in velocity or frequency when this function is moved to fitter. (2011/03/30 WK)
    4848 
    4849   // xArray : contains elemental values for computing the least-square matrix.
    4850   //          xArray.size() is nDOF and xArray[*].size() is nChan.
    4851   //          Each xArray element are as follows:
    4852   //          xArray[0]    = {1.0, 1.0, 1.0, ..., 1.0},
    4853   //          xArray[2n-1] = {sin(nPI/L*x[0]), sin(nPI/L*x[1]), ..., sin(nPI/L*x[nChan])},
    4854   //          xArray[2n]   = {cos(nPI/L*x[0]), cos(nPI/L*x[1]), ..., cos(nPI/L*x[nChan])},
    4855   //          where (1 <= n <= nMaxWavesInSW),
    4856   //          or,
    4857   //          xArray[2n-1] = {sin(wn[n]PI/L*x[0]), sin(wn[n]PI/L*x[1]), ..., sin(wn[n]PI/L*x[nChan])},
    4858   //          xArray[2n]   = {cos(wn[n]PI/L*x[0]), cos(wn[n]PI/L*x[1]), ..., cos(wn[n]PI/L*x[nChan])},
    4859   //          where wn[n] denotes waveNumbers[n] (1 <= n <= waveNumbers.size()).
    4860   std::vector<std::vector<double> > xArray;
    4861   if (hasConstantTerm) {
    4862     std::vector<double> xu;
    4863     for (int j = 0; j < nChan; ++j) {
    4864       xu.push_back(1.0);
    4865     }
    4866     xArray.push_back(xu);
    4867   }
     4706  double stretch0 = 2.0*PI/(double)(nchan-1);
     4707
    48684708  for (uInt i = (hasConstantTerm ? 1 : 0); i < nWaves.size(); ++i) {
    4869     double xFactor = baseXFactor*(double)nWaves[i];
    4870     std::vector<double> xs, xc;
    4871     xs.clear();
    4872     xc.clear();
    4873     for (int j = 0; j < nChan; ++j) {
    4874       xs.push_back(sin(xFactor*(double)j));
    4875       xc.push_back(cos(xFactor*(double)j));
    4876     }
    4877     xArray.push_back(xs);
    4878     xArray.push_back(xc);
    4879   }
    4880 
    4881   std::vector<double> z1, r1, residual;
    4882   for (int i = 0; i < nChan; ++i) {
    4883     z1.push_back((double)data[i]);
    4884     r1.push_back(0.0);
    4885     residual.push_back(0.0);
    4886   }
    4887 
    4888   for (int nClip = 0; nClip < nIterClip+1; ++nClip) {
    4889     // xMatrix : horizontal concatenation of
    4890     //           the least-sq. matrix (left) and an
    4891     //           identity matrix (right).
    4892     // the right part is used to calculate the inverse matrix of the left part.
    4893     double xMatrix[nDOF][2*nDOF];
    4894     double zMatrix[nDOF];
    4895     for (int i = 0; i < nDOF; ++i) {
    4896       for (int j = 0; j < 2*nDOF; ++j) {
    4897         xMatrix[i][j] = 0.0;
    4898       }
    4899       xMatrix[i][nDOF+i] = 1.0;
    4900       zMatrix[i] = 0.0;
    4901     }
    4902 
    4903     int nUseData = 0;
    4904     for (int k = 0; k < nChan; ++k) {
    4905       if (maskArray[k] == 0) continue;
    4906 
    4907       for (int i = 0; i < nDOF; ++i) {
    4908         for (int j = i; j < nDOF; ++j) {
    4909           xMatrix[i][j] += xArray[i][k] * xArray[j][k];
    4910         }
    4911         zMatrix[i] += z1[k] * xArray[i][k];
    4912       }
    4913 
    4914       nUseData++;
    4915     }
    4916 
    4917     if (nUseData < 1) {
    4918         throw(AipsError("all channels clipped or masked. can't execute fitting anymore."));     
    4919     }
    4920 
    4921     for (int i = 0; i < nDOF; ++i) {
    4922       for (int j = 0; j < i; ++j) {
    4923         xMatrix[i][j] = xMatrix[j][i];
    4924       }
    4925     }
    4926 
    4927     std::vector<double> invDiag;
    4928     for (int i = 0; i < nDOF; ++i) {
    4929       invDiag.push_back(1.0/xMatrix[i][i]);
    4930       for (int j = 0; j < nDOF; ++j) {
    4931         xMatrix[i][j] *= invDiag[i];
    4932       }
    4933     }
    4934 
    4935     for (int k = 0; k < nDOF; ++k) {
    4936       for (int i = 0; i < nDOF; ++i) {
    4937         if (i != k) {
    4938           double factor1 = xMatrix[k][k];
    4939           double factor2 = xMatrix[i][k];
    4940           for (int j = k; j < 2*nDOF; ++j) {
    4941             xMatrix[i][j] *= factor1;
    4942             xMatrix[i][j] -= xMatrix[k][j]*factor2;
    4943             xMatrix[i][j] /= factor1;
    4944           }
    4945         }
    4946       }
    4947       double xDiag = xMatrix[k][k];
    4948       for (int j = k; j < 2*nDOF; ++j) {
    4949         xMatrix[k][j] /= xDiag;
    4950       }
    4951     }
    4952    
    4953     for (int i = 0; i < nDOF; ++i) {
    4954       for (int j = 0; j < nDOF; ++j) {
    4955         xMatrix[i][nDOF+j] *= invDiag[j];
    4956       }
    4957     }
    4958     //compute a vector y which consists of the coefficients of the sinusoids forming the
    4959     //best-fit curves (a0,s1,c1,s2,c2,...), where a0 is constant and s* and c* are of sine
    4960     //and cosine functions, respectively.
    4961     std::vector<double> y;
    4962     params.clear();
    4963     for (int i = 0; i < nDOF; ++i) {
    4964       y.push_back(0.0);
    4965       for (int j = 0; j < nDOF; ++j) {
    4966         y[i] += xMatrix[i][nDOF+j]*zMatrix[j];
    4967       }
    4968       params.push_back(y[i]);
    4969     }
    4970 
    4971     for (int i = 0; i < nChan; ++i) {
    4972       r1[i] = y[0];
    4973       for (int j = 1; j < nDOF; ++j) {
    4974         r1[i] += y[j]*xArray[j][i];
    4975       }
    4976       residual[i] = z1[i] - r1[i];
    4977     }
    4978 
    4979     double stdDev = 0.0;
    4980     for (int i = 0; i < nChan; ++i) {
    4981       stdDev += residual[i]*residual[i]*(double)maskArray[i];
    4982     }
    4983     stdDev = sqrt(stdDev/(double)nData);
    4984     rms = (float)stdDev;
    4985 
    4986     if ((nClip == nIterClip) || (thresClip <= 0.0)) {
    4987       break;
    4988     } else {
    4989 
    4990       double thres = stdDev * thresClip;
    4991       int newNData = 0;
    4992       for (int i = 0; i < nChan; ++i) {
    4993         if (abs(residual[i]) >= thres) {
    4994           maskArray[i] = 0;
    4995           finalMask[i] = false;
    4996         }
    4997         if (maskArray[i] > 0) {
    4998           newNData++;
    4999         }
    5000       }
    5001       if (newNData == nData) {
    5002         break; //no more flag to add. iteration stops.
    5003       } else {
    5004         nData = newNData;
    5005       }
    5006 
    5007     }
    5008   }
    5009 
    5010   nClipped = initNData - nData;
    5011 
    5012   std::vector<float> result;
    5013   if (getResidual) {
    5014     for (int i = 0; i < nChan; ++i) {
    5015       result.push_back((float)residual[i]);
    5016     }
    5017   } else {
    5018     for (int i = 0; i < nChan; ++i) {
    5019       result.push_back((float)r1[i]);
    5020     }
    5021   }
    5022 
    5023   return result;
    5024 }
    5025 
    5026 void Scantable::fitBaseline(const std::vector<bool>& mask, int whichrow, Fitter& fitter)
    5027 {
    5028   std::vector<double> dAbcissa = getAbcissa(whichrow);
    5029   std::vector<float> abcissa;
    5030   for (uInt i = 0; i < dAbcissa.size(); ++i) {
    5031     abcissa.push_back((float)dAbcissa[i]);
    5032   }
    5033   std::vector<float> spec = getSpectrum(whichrow);
    5034 
    5035   fitter.setData(abcissa, spec, mask);
    5036   fitter.lfit();
    5037 }
    5038 
    5039 std::vector<bool> Scantable::getCompositeChanMask(int whichrow, const std::vector<bool>& inMask)
    5040 {
    5041   /****
    5042   double ms1TimeStart, ms1TimeEnd, ms2TimeStart, ms2TimeEnd;
    5043   double elapse1 = 0.0;
    5044   double elapse2 = 0.0;
    5045 
    5046   ms1TimeStart = mathutil::gettimeofday_sec();
    5047   ****/
    5048 
     4709    int sidx = hasConstantTerm ? 2*i-1 : 2*i;
     4710    int cidx = sidx + 1;
     4711    double stretch = stretch0*(double)nWaves[i];
     4712
     4713    for (int j = 0; j < nchan; ++j) {
     4714      model[sidx][j] = sin(stretch*(double)j);
     4715      model[cidx][j] = cos(stretch*(double)j);
     4716    }
     4717  }
     4718
     4719  return model;
     4720}
     4721
     4722std::vector<bool> Scantable::getCompositeChanMask(int whichrow,
     4723                                                  const std::vector<bool>& inMask)
     4724{
    50494725  std::vector<bool> mask = getMask(whichrow);
    5050 
    5051   /****
    5052   ms1TimeEnd = mathutil::gettimeofday_sec();
    5053   elapse1 = ms1TimeEnd - ms1TimeStart;
    5054   std::cout << "ms1   : " << elapse1 << " (sec.)" << endl;
    5055   ms2TimeStart = mathutil::gettimeofday_sec();
    5056   ****/
    5057 
    50584726  uInt maskSize = mask.size();
    50594727  if (inMask.size() != 0) {
    … …  
    50664734  }
    50674735
    5068   /****
    5069   ms2TimeEnd = mathutil::gettimeofday_sec();
    5070   elapse2 = ms2TimeEnd - ms2TimeStart;
    5071   std::cout << "ms2   : " << elapse2 << " (sec.)" << endl;
    5072   ****/
    5073 
    50744736  return mask;
    50754737}
    50764738
    5077 /*
    5078 std::vector<bool> Scantable::getCompositeChanMask(int whichrow, const std::vector<bool>& inMask, const std::vector<int>& edge, const int minEdgeSize, STLineFinder& lineFinder)
    5079 {
    5080   int edgeSize = edge.size();
    5081   std::vector<int> currentEdge;
    5082   if (edgeSize >= 2) {
    5083       int idx = 0;
    5084       if (edgeSize > 2) {
    5085         if (edgeSize < minEdgeSize) {
    5086           throw(AipsError("Length of edge element info is less than that of IFs"));
    5087         }
    5088         idx = 2 * getIF(whichrow);
    5089       }
    5090       currentEdge.push_back(edge[idx]);
    5091       currentEdge.push_back(edge[idx+1]);
    5092   } else {
    5093     throw(AipsError("Wrong length of edge element"));
    5094   }
     4739std::vector<bool> Scantable::getCompositeChanMask(int whichrow,
     4740                                                  const std::vector<bool>& inMask,
     4741                                                  const std::vector<int>& edge,
     4742                                                  std::vector<int>& currEdge,
     4743                                                  STLineFinder& lineFinder)
     4744{
     4745  std::vector<uint> ifNos = getIFNos();
     4746  if (edge.size() < ifNos.size()*2) {
     4747    throw(AipsError("Length of edge element info is less than that of IFs"));
     4748  }
     4749
     4750  int ifVal = getIF(whichrow);
     4751  bool foundIF = false;
     4752  uint idx;
     4753  for (uint i = 0; i < ifNos.size(); ++i) {
     4754    if (ifVal == (int)ifNos[i]) {
     4755      idx = 2*i;
     4756      foundIF = true;
     4757      break;
     4758    }
     4759  }
     4760  if (!foundIF) {
     4761    throw(AipsError("bad IF number"));
     4762  }
     4763
     4764  currEdge.clear();
     4765  currEdge.resize(2);
     4766  currEdge[0] = edge[idx];
     4767  currEdge[1] = edge[idx+1];
    50954768
    50964769  lineFinder.setData(getSpectrum(whichrow));
    5097   lineFinder.findLines(getCompositeChanMask(whichrow, inMask), currentEdge, whichrow);
     4770  lineFinder.findLines(getCompositeChanMask(whichrow, inMask), currEdge, whichrow);
    50984771
    50994772  return lineFinder.getMask();
    51004773}
    5101 */
    5102 
    5103 /* for poly. the variations of outputFittingResult() should be merged into one eventually (2011/3/10 WK)  */
    5104 void Scantable::outputFittingResult(bool outLogger, bool outTextFile, bool csvFormat, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, Fitter& fitter)
    5105 {
    5106   if (outLogger || outTextFile) {
    5107     std::vector<float> params = fitter.getParameters();
    5108     std::vector<bool>  fixed  = fitter.getFixedParameters();
    5109     float rms = getRms(chanMask, whichrow);
    5110     String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
    5111 
    5112     if (outLogger) {
    5113       LogIO ols(LogOrigin("Scantable", funcName, WHERE));
    5114       ols << formatBaselineParams(params, fixed, rms, -1, masklist, whichrow, false, csvFormat) << LogIO::POST ;
    5115     }
    5116     if (outTextFile) {
    5117       ofs << formatBaselineParams(params, fixed, rms, -1, masklist, whichrow, true, csvFormat) << flush;
    5118     }
    5119   }
    5120 }
    51214774
    51224775/* for cspline. will be merged once cspline is available in fitter (2011/3/10 WK) */
    5123 void Scantable::outputFittingResult(bool outLogger, bool outTextFile, bool csvFormat, 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)
    5124 {
    5125   if (outLogger || outTextFile) {
    5126   /****
    5127   double ms1TimeStart, ms1TimeEnd, ms2TimeStart, ms2TimeEnd;
    5128   double elapse1 = 0.0;
    5129   double elapse2 = 0.0;
    5130 
    5131   ms1TimeStart = mathutil::gettimeofday_sec();
    5132   ****/
    5133 
    5134     float rms = getRms(chanMask, whichrow);
    5135 
    5136   /****
    5137   ms1TimeEnd = mathutil::gettimeofday_sec();
    5138   elapse1 = ms1TimeEnd - ms1TimeStart;
    5139   std::cout << "ot1   : " << elapse1 << " (sec.)" << endl;
    5140   ms2TimeStart = mathutil::gettimeofday_sec();
    5141   ****/
    5142 
    5143     String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
    5144     std::vector<bool> fixed;
    5145     fixed.clear();
    5146 
    5147     if (outLogger) {
    5148       LogIO ols(LogOrigin("Scantable", funcName, WHERE));
    5149       ols << formatPiecewiseBaselineParams(edge, params, fixed, rms, nClipped, masklist, whichrow, false, csvFormat) << LogIO::POST ;
    5150     }
    5151     if (outTextFile) {
    5152       ofs << formatPiecewiseBaselineParams(edge, params, fixed, rms, nClipped, masklist, whichrow, true, csvFormat) << flush;
    5153     }
    5154 
    5155   /****
    5156   ms2TimeEnd = mathutil::gettimeofday_sec();
    5157   elapse2 = ms2TimeEnd - ms2TimeStart;
    5158   std::cout << "ot2   : " << elapse2 << " (sec.)" << endl;
    5159   ****/
    5160 
    5161   }
    5162 }
    5163 
    5164 /* for chebyshev/sinusoid. will be merged once chebyshev/sinusoid is available in fitter (2011/3/10 WK) */
    5165 void Scantable::outputFittingResult(bool outLogger, bool outTextFile, bool csvFormat, 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)
     4776void Scantable::outputFittingResult(bool outLogger,
     4777                                    bool outTextFile,
     4778                                    bool csvFormat,
     4779                                    const std::vector<bool>& chanMask,
     4780                                    int whichrow,
     4781                                    const casa::String& coordInfo,
     4782                                    bool hasSameNchan,
     4783                                    ofstream& ofs,
     4784                                    const casa::String& funcName,
     4785                                    const std::vector<int>& edge,
     4786                                    const std::vector<float>& params,
     4787                                    const int nClipped)
    51664788{
    51674789  if (outLogger || outTextFile) {
    … …  
    51734795    if (outLogger) {
    51744796      LogIO ols(LogOrigin("Scantable", funcName, WHERE));
    5175       ols << formatBaselineParams(params, fixed, rms, nClipped, masklist, whichrow, false, csvFormat) << LogIO::POST ;
     4797      ols << formatPiecewiseBaselineParams(edge, params, fixed, rms, nClipped,
     4798                                           masklist, whichrow, false, csvFormat) << LogIO::POST ;
    51764799    }
    51774800    if (outTextFile) {
    5178       ofs << formatBaselineParams(params, fixed, rms, nClipped, masklist, whichrow, true, csvFormat) << flush;
     4801      ofs << formatPiecewiseBaselineParams(edge, params, fixed, rms, nClipped,
     4802                                           masklist, whichrow, true, csvFormat) << flush;
     4803    }
     4804  }
     4805}
     4806
     4807/* for poly/chebyshev/sinusoid. */
     4808void Scantable::outputFittingResult(bool outLogger,
     4809                                    bool outTextFile,
     4810                                    bool csvFormat,
     4811                                    const std::vector<bool>& chanMask,
     4812                                    int whichrow,
     4813                                    const casa::String& coordInfo,
     4814                                    bool hasSameNchan,
     4815                                    ofstream& ofs,
     4816                                    const casa::String& funcName,
     4817                                    const std::vector<float>& params,
     4818                                    const int nClipped)
     4819{
     4820  if (outLogger || outTextFile) {
     4821    float rms = getRms(chanMask, whichrow);
     4822    String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
     4823    std::vector<bool> fixed;
     4824    fixed.clear();
     4825
     4826    if (outLogger) {
     4827      LogIO ols(LogOrigin("Scantable", funcName, WHERE));
     4828      ols << formatBaselineParams(params, fixed, rms, nClipped,
     4829                                  masklist, whichrow, false, csvFormat) << LogIO::POST ;
     4830    }
     4831    if (outTextFile) {
     4832      ofs << formatBaselineParams(params, fixed, rms, nClipped,
     4833                                  masklist, whichrow, true, csvFormat) << flush;
    51794834    }
    51804835  }

  • trunk/src/Scantable.h

    r2767 r2773  
    3939#include "STFit.h"
    4040#include "STFitEntry.h"
    41 //#include "STFitter.h"
    4241#include "STFocus.h"
    4342#include "STFrequencies.h"
    … …  
    5352
    5453class Fitter;
     54class STLineFinder;
     55class STBaselineTable;
    5556
    5657/**
    … …  
    635636  static std::vector<bool> getMaskFromMaskList(const int nchan,
    636637                                               const std::vector<int>& masklist);
     638  static casa::Vector<casa::uInt> getMaskListFromMask(const std::vector<bool>& mask);
    637639  static std::vector<int> splitToIntList(const std::string& str, const char delim);
    638640  static std::vector<string> splitToStringList(const std::string& str, const char delim);
    … …  
    752754                                                      const casa::Array<T2>&);
    753755
    754   void fitBaseline(const std::vector<bool>& mask, int whichrow, Fitter& fitter);
    755756  double getNormalPolynomial(int n, double x);
    756757  double getChebyshevPolynomial(int n, double x);
     758  std::vector<std::vector<double> > getPolynomialModel(int order,
     759                                                       int nchan,
     760                                                       double (Scantable::*pfunc)(int, double));
     761  std::vector<std::vector<std::vector<double> > > getPolynomialModelReservoir(int order,
     762                                                                              double (Scantable::*pfunc)(int, double),
     763                                                                              std::vector<int>& nChanNos);
    757764  std::vector<float> doPolynomialFitting(const std::vector<float>& data,
    758765                                         const std::vector<bool>& mask,
    … …  
    793800  std::vector<float> doLeastSquareFitting(const std::vector<float>& data,
    794801                                          const std::vector<bool>& mask,
    795                                           double (Scantable::*pfunc)(int, double),
    796                                           int order,
     802                                          const std::vector<std::vector<double> >& model,
    797803                                          std::vector<float>& params,
    798804                                          float& rms,
    … …  
    831837                                          int nIterClip=0,
    832838                                          bool getResidual=true);
     839  std::vector<float> doCubicSplineLeastSquareFitting(const std::vector<float>& data,
     840                                                     const std::vector<bool>& mask,
     841                                                     const std::vector<std::vector<double> >& model,
     842                                                     int nPiece,
     843                                                     bool useGivenPieceBoundary,
     844                                                     std::vector<int>& idxEdge,
     845                                                     std::vector<float>& params,
     846                                                     float& rms,
     847                                                     std::vector<bool>& finalMask,
     848                                                     int& nClipped,
     849                                                     float thresClip=3.0,
     850                                                     int nIterClip=0,
     851                                                     bool getResidual=true);
    833852  std::vector<float> doSinusoidFitting(const std::vector<float>& data,
    834853                                       const std::vector<bool>& mask,
    … …  
    849868                                       int nIterClip=0,
    850869                                       bool getResidual=true);
    851   void selectWaveNumbers(const int whichrow,
    852                          const std::vector<bool>& chanMask,
    853                          const std::string& fftInfo,
    854                          //const bool applyFFT,
    855                          //const std::string& fftMethod,
    856                          //const std::string& fftThresh,
    857                          const std::vector<int>& addNWaves,
    858                          const std::vector<int>& rejectNWaves,
    859                          std::vector<int>& nWaves);
     870  std::vector<std::vector<double> > getSinusoidModel(const std::vector<int>& waveNumbers, int nchan);
     871  std::vector<std::vector<std::vector<double> > > getSinusoidModelReservoir(const std::vector<int>& waveNumbers,
     872                                                                            std::vector<int>& nChanNos);
     873  std::vector<int> selectWaveNumbers(const std::vector<int>& addNWaves,
     874                                     const std::vector<int>& rejectNWaves);
     875  std::vector<int> selectWaveNumbers(const int whichrow,
     876                                     const std::vector<bool>& chanMask,
     877                                     //const std::string& fftInfo,
     878                                     const bool applyFFT,
     879                                     const std::string& fftMethod,
     880                                     const std::string& fftThresh,
     881                                     const std::vector<int>& addNWaves,
     882                                     const std::vector<int>& rejectNWaves);
     883  int getIdxOfNchan(const int nChan, const std::vector<int>& nChanNos);
    860884  void parseFFTInfo(const std::string& fftInfo,
    861885                    bool& applyFFT,
    862886                    std::string& fftMethod,
    863887                    std::string& fftThresh);
    864   void parseThresholdExpression(const std::string& fftThresh,
    865                                 std::string& fftThAttr,
    866                                 float& fftThSigma,
    867                                 int& fftThTop);
     888  void parseFFTThresholdInfo(const std::string& fftThresh,
     889                             std::string& fftThAttr,
     890                             float& fftThSigma,
     891                             int& fftThTop);
    868892  void doSelectWaveNumbers(const int whichrow,
    869893                           const std::vector<bool>& chanMask,
    … …  
    879903  void setWaveNumberListUptoNyquistFreq(const int whichrow,
    880904                                        std::vector<int>& nWaves);
     905  void initialiseBaselining(const std::string& blfile,
     906                            std::ofstream& ofs,
     907                            const bool outLogger,
     908                            bool& outTextFile,
     909                            bool& csvFormat,
     910                            casa::String& coordInfo,
     911                            bool& hasSameNchan,
     912                            const std::string& progressInfo,
     913                            bool& showProgress,
     914                            int& minNRow,
     915                            casa::Vector<casa::Double>& timeSecCol);
     916  void finaliseBaselining(const bool outBaselineTable,
     917                          STBaselineTable* pbt,
     918                          const string& bltable,
     919                          const bool outTextFile,
     920                          std::ofstream& ofs);
     921  void initLineFinder(const std::vector<int>& edge,
     922                      const float threshold,
     923                      const int chanAvgLimit,
     924                      STLineFinder& lineFinder);
    881925  bool hasSameNchanOverIFs();
    882926  std::string getMaskRangeList(const std::vector<bool>& mask,
    … …  
    889933  std::string formatBaselineParamsFooter(float rms, int nClipped, bool verbose, bool csvformat) const;
    890934  std::vector<bool> getCompositeChanMask(int whichrow, const std::vector<bool>& inMask);
    891   //std::vector<bool> getCompositeChanMask(int whichrow, const std::vector<bool>& inMask, const std::vector<int>& edge, const int minEdgeSize, STLineFinder& lineFinder);
    892   void outputFittingResult(bool outLogger, bool outTextFile, bool csvFormat, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, std::ofstream& ofs, const casa::String& funcName, Fitter& fitter);
     935  std::vector<bool> getCompositeChanMask(int whichrow, const std::vector<bool>& inMask, const std::vector<int>& edge, std::vector<int>& currEdge, STLineFinder& lineFinder);
    893936  void outputFittingResult(bool outLogger, bool outTextFile, bool csvFormat, 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);
    894937  void outputFittingResult(bool outLogger, bool outTextFile, bool csvFormat, 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);

  • trunk/test/test_scantable.py

    r2672 r2773  
    221221        res_rms = (q-b).stats('rms')
    222222        assert_almost_equals(res_rms[0], 0.38346, 5)
    223         assert_almost_equals(res_rms[1], 0.38780, 5)       
     223        #assert_almost_equals(res_rms[1], 0.38780, 5)       
     224        assert_almost_equals(res_rms[1], 0.48780, 5)       
    224225
    225226
Note: See TracChangeset for help on using the changeset viewer.