Ignore:
Timestamp:
09/04/12 18:58:06 (12 years ago)
Author:
WataruKawasaki
Message:

New Development: Yes

JIRA Issue: Yes CAS-4145

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: Yes

Module(s): scantable

Description: added scantable methods [auto_]chebyshev_baseline() to subtract baseline using Chebyshev polynomials.


File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/Scantable.cpp

    r2644 r2645  
    26022602}
    26032603
    2604 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)
    2605 {
    2606   /*************
    2607   double totTimeStart, totTimeEnd, blTimeStart, blTimeEnd, ioTimeStart, ioTimeEnd, msTimeStart, msTimeEnd, seTimeStart, seTimeEnd, otTimeStart, otTimeEnd, prTimeStart, prTimeEnd;
    2608   double elapseMs = 0.0;
    2609   double elapseSe = 0.0;
    2610   double elapseOt = 0.0;
    2611   double elapsePr = 0.0;
    2612   double elapseBl = 0.0;
    2613   double elapseIo = 0.0;
    2614   totTimeStart = mathutil::gettimeofday_sec();
    2615   *************/
    2616 
     2604void 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)
     2605{
    26172606  try {
    26182607    ofstream ofs;
    … …  
    26342623    }
    26352624
    2636     //Fitter fitter = Fitter();
    2637     //fitter.setExpression("cspline", nPiece);
    2638     //fitter.setIterClipping(thresClip, nIterClip);
    2639 
    26402625    bool showProgress;
    26412626    int minNRow;
    … …  
    26452630    std::vector<bool> chanMask;
    26462631
    2647     //--------------------------------
    26482632    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
    2649       /******************
    2650       ioTimeStart = mathutil::gettimeofday_sec();
    2651       **/
    26522633      std::vector<float> sp = getSpectrum(whichrow);
    2653       /**
    2654       ioTimeEnd = mathutil::gettimeofday_sec();
    2655       elapseIo += (double)(ioTimeEnd - ioTimeStart);
    2656       msTimeStart = mathutil::gettimeofday_sec();
    2657       ******************/
    2658 
    26592634      chanMask = getCompositeChanMask(whichrow, mask);
    2660 
    2661       /**
    2662       msTimeEnd = mathutil::gettimeofday_sec();
    2663       elapseMs += (double)(msTimeEnd - msTimeStart);
    2664       blTimeStart = mathutil::gettimeofday_sec();
    2665       **/
    2666 
    2667       //fitBaseline(chanMask, whichrow, fitter);
    2668       //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
    2669       std::vector<int> pieceEdges(nPiece+1);
    2670       std::vector<float> params(nPiece*4);
     2635      std::vector<float> params(order+1);
    26712636      int nClipped = 0;
    2672       std::vector<float> res = doCubicSplineFitting(sp, chanMask, nPiece, pieceEdges, params, nClipped, thresClip, nIterClip, getResidual);
    2673 
    2674       /**
    2675       blTimeEnd = mathutil::gettimeofday_sec();
    2676       elapseBl += (double)(blTimeEnd - blTimeStart);
    2677       seTimeStart = mathutil::gettimeofday_sec();
    2678       **/
    2679 
     2637      std::vector<float> res = doChebyshevFitting(sp, chanMask, order, params, nClipped, thresClip, nIterClip, getResidual);
    26802638
    26812639      setSpectrum(res, whichrow);
    2682       //
    2683 
    2684       /**
    2685       seTimeEnd = mathutil::gettimeofday_sec();
    2686       elapseSe += (double)(seTimeEnd - seTimeStart);
    2687       otTimeStart = mathutil::gettimeofday_sec();
    2688       **/
    2689 
    2690       outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()", pieceEdges, params, nClipped);
    2691 
    2692       /**
    2693       otTimeEnd = mathutil::gettimeofday_sec();
    2694       elapseOt += (double)(otTimeEnd - otTimeStart);
    2695       prTimeStart = mathutil::gettimeofday_sec();
    2696       **/
    2697 
     2640      outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "chebyshevBaseline()", params, nClipped);
    26982641      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
    2699 
    2700       /******************
    2701       prTimeEnd = mathutil::gettimeofday_sec();
    2702       elapsePr += (double)(prTimeEnd - prTimeStart);
    2703       ******************/
    2704     }
    2705     //--------------------------------
     2642    }
    27062643   
    27072644    if (outTextFile) ofs.close();
    … …  
    27102647    throw;
    27112648  }
    2712   /***************
    2713   totTimeEnd = mathutil::gettimeofday_sec();
    2714   std::cout << "io    : " << elapseIo << " (sec.)" << endl;
    2715   std::cout << "ms    : " << elapseMs << " (sec.)" << endl;
    2716   std::cout << "bl    : " << elapseBl << " (sec.)" << endl;
    2717   std::cout << "se    : " << elapseSe << " (sec.)" << endl;
    2718   std::cout << "ot    : " << elapseOt << " (sec.)" << endl;
    2719   std::cout << "pr    : " << elapsePr << " (sec.)" << endl;
    2720   std::cout << "total : " << (double)(totTimeEnd - totTimeStart) << " (sec.)" << endl;
    2721   ***************/
    2722 }
    2723 
    2724 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)
     2649}
     2650
     2651void 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)
    27252652{
    27262653  try {
    … …  
    27422669      hasSameNchan = hasSameNchanOverIFs();
    27432670    }
    2744 
    2745     //Fitter fitter = Fitter();
    2746     //fitter.setExpression("cspline", nPiece);
    2747     //fitter.setIterClipping(thresClip, nIterClip);
    27482671
    27492672    int nRow = nrow();
    … …  
    27872710      //fitBaseline(chanMask, whichrow, fitter);
    27882711      //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
    2789       std::vector<int> pieceEdges(nPiece+1);
    2790       std::vector<float> params(nPiece*4);
     2712      std::vector<float> params(order+1);
    27912713      int nClipped = 0;
    2792       std::vector<float> res = doCubicSplineFitting(sp, chanMask, nPiece, pieceEdges, params, nClipped, thresClip, nIterClip, getResidual);
     2714      std::vector<float> res = doChebyshevFitting(sp, chanMask, order, params, nClipped, thresClip, nIterClip, getResidual);
    27932715      setSpectrum(res, whichrow);
    2794       //
    2795 
    2796       outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()", pieceEdges, params, nClipped);
     2716
     2717      outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoChebyshevBaseline()", params, nClipped);
    27972718      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
    27982719    }
    … …  
    28052726}
    28062727
    2807 std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data, const std::vector<bool>& mask, int nPiece, std::vector<int>& idxEdge, std::vector<float>& params, int& nClipped, float thresClip, int nIterClip, bool getResidual)
     2728  /*
     2729double Scantable::getChebyshevPolynomial(int n, double x) {
     2730  if ((x < -1.0)||(x > 1.0)) {
     2731    throw(AipsError("out of definition range (-1 <= x <= 1)."));
     2732  } else if (n < 0) {
     2733    throw(AipsError("the order must be zero or positive."));
     2734  } else if (n == 0) {
     2735    return 1.0;
     2736  } else if (n == 1) {
     2737    return x;
     2738  } else {
     2739    return 2.0*x*getChebyshevPolynomial(n-1, x) - getChebyshevPolynomial(n-2, x);
     2740  }
     2741}
     2742  */
     2743double Scantable::getChebyshevPolynomial(int n, double x) {
     2744  if ((x < -1.0)||(x > 1.0)) {
     2745    throw(AipsError("out of definition range (-1 <= x <= 1)."));
     2746  } else if (n < 0) {
     2747    throw(AipsError("the order must be zero or positive."));
     2748  } else if (n == 0) {
     2749    return 1.0;
     2750  } else if (n == 1) {
     2751    return x;
     2752  } else {
     2753    double res = 0.0;
     2754    for (int m = 0; m <= n/2; ++m) {
     2755      double c = 1.0;
     2756      if (m > 0) {
     2757        for (int i = 1; i <= m; ++i) {
     2758          c *= (double)(n-2*m+i)/(double)i;
     2759        }
     2760      }
     2761      res += (m%2 == 0 ? 1.0 : -1.0)*(double)n/(double)(n-m)*pow(2.0*x, (double)(n-2*m-1))/2.0*c;
     2762    }
     2763    return res;
     2764  }
     2765}
     2766
     2767std::vector<float> Scantable::doChebyshevFitting(const std::vector<float>& data, const std::vector<bool>& mask, int order, std::vector<float>& params, int& nClipped, float thresClip, int nIterClip, bool getResidual)
    28082768{
    28092769  if (data.size() != mask.size()) {
    28102770    throw(AipsError("data and mask sizes are not identical"));
    28112771  }
    2812   if (nPiece < 1) {
    2813     throw(AipsError("number of the sections must be one or more"));
     2772  if (order < 0) {
     2773    throw(AipsError("maximum order of Chebyshev polynomial must not be negative."));
    28142774  }
    28152775
    28162776  int nChan = data.size();
    2817   std::vector<int> maskArray(nChan);
    2818   std::vector<int> x(nChan);
    2819   int j = 0;
     2777  std::vector<int> maskArray;
     2778  std::vector<int> x;
    28202779  for (int i = 0; i < nChan; ++i) {
    2821     maskArray[i] = mask[i] ? 1 : 0;
     2780    maskArray.push_back(mask[i] ? 1 : 0);
    28222781    if (mask[i]) {
    2823       x[j] = i;
    2824       j++;
    2825     }
    2826   }
    2827   int initNData = j;
    2828 
    2829   if (initNData < nPiece) {
    2830     throw(AipsError("too few non-flagged channels"));
    2831   }
    2832 
    2833   int nElement = (int)(floor(floor((double)(initNData/nPiece))+0.5));
    2834   std::vector<double> invEdge(nPiece-1);
    2835   idxEdge[0] = x[0];
    2836   for (int i = 1; i < nPiece; ++i) {
    2837     int valX = x[nElement*i];
    2838     idxEdge[i] = valX;
    2839     invEdge[i-1] = 1.0/(double)valX;
    2840   }
    2841   idxEdge[nPiece] = x[initNData-1]+1;
     2782      x.push_back(i);
     2783    }
     2784  }
     2785
     2786  int initNData = x.size();
    28422787
    28432788  int nData = initNData;
    2844   int nDOF = nPiece + 3;  //number of parameters to solve, namely, 4+(nPiece-1).
    2845 
    2846   std::vector<double> x1(nChan), x2(nChan), x3(nChan);
    2847   std::vector<double> z1(nChan), x1z1(nChan), x2z1(nChan), x3z1(nChan);
    2848   std::vector<double> r1(nChan), residual(nChan);
     2789  int nDOF = order + 1;  //number of parameters to solve.
     2790
     2791  // xArray : contains elemental values for computing the least-square matrix.
     2792  //          xArray.size() is nDOF and xArray[*].size() is nChan.
     2793  //          Each xArray element are as follows:
     2794  //          xArray[0]   = {T0(-1), T0(2/(nChan-1)-1), T0(4/(nChan-1)-1), ..., T0(1)},
     2795  //          xArray[n-1] = ...,
     2796  //          xArray[n]   = {Tn(-1), Tn(2/(nChan-1)-1), Tn(4/(nChan-1)-1), ..., Tn(1)}
     2797  //          where (0 <= n <= order),
     2798  std::vector<std::vector<double> > xArray;
     2799  for (int i = 0; i < nDOF; ++i) {
     2800    double xFactor = 2.0/(double)(nChan - 1);
     2801    std::vector<double> xs;
     2802    xs.clear();
     2803    for (int j = 0; j < nChan; ++j) {
     2804      xs.push_back(getChebyshevPolynomial(i, xFactor*(double)j-1.0));
     2805    }
     2806    xArray.push_back(xs);
     2807  }
     2808
     2809  std::vector<double> z1, r1, residual;
    28492810  for (int i = 0; i < nChan; ++i) {
    2850     double di = (double)i;
    2851     double dD = (double)data[i];
    2852     x1[i]   = di;
    2853     x2[i]   = di*di;
    2854     x3[i]   = di*di*di;
    2855     z1[i]   = dD;
    2856     x1z1[i] = dD*di;
    2857     x2z1[i] = dD*di*di;
    2858     x3z1[i] = dD*di*di*di;
    2859     r1[i]   = 0.0;
    2860     residual[i] = 0.0;
     2811    z1.push_back((double)data[i]);
     2812    r1.push_back(0.0);
     2813    residual.push_back(0.0);
    28612814  }
    28622815
    … …  
    28762829    }
    28772830
    2878     for (int n = 0; n < nPiece; ++n) {
    2879       int nUseDataInPiece = 0;
    2880       for (int i = idxEdge[n]; i < idxEdge[n+1]; ++i) {
    2881 
    2882         if (maskArray[i] == 0) continue;
    2883 
    2884         xMatrix[0][0] += 1.0;
    2885         xMatrix[0][1] += x1[i];
    2886         xMatrix[0][2] += x2[i];
    2887         xMatrix[0][3] += x3[i];
    2888         xMatrix[1][1] += x2[i];
    2889         xMatrix[1][2] += x3[i];
    2890         xMatrix[1][3] += x2[i]*x2[i];
    2891         xMatrix[2][2] += x2[i]*x2[i];
    2892         xMatrix[2][3] += x3[i]*x2[i];
    2893         xMatrix[3][3] += x3[i]*x3[i];
    2894         zMatrix[0] += z1[i];
    2895         zMatrix[1] += x1z1[i];
    2896         zMatrix[2] += x2z1[i];
    2897         zMatrix[3] += x3z1[i];
    2898 
    2899         for (int j = 0; j < n; ++j) {
    2900           double q = 1.0 - x1[i]*invEdge[j];
    2901           q = q*q*q;
    2902           xMatrix[0][j+4] += q;
    2903           xMatrix[1][j+4] += q*x1[i];
    2904           xMatrix[2][j+4] += q*x2[i];
    2905           xMatrix[3][j+4] += q*x3[i];
    2906           for (int k = 0; k < j; ++k) {
    2907             double r = 1.0 - x1[i]*invEdge[k];
    2908             r = r*r*r;
    2909             xMatrix[k+4][j+4] += r*q;
    2910           }
    2911           xMatrix[j+4][j+4] += q*q;
    2912           zMatrix[j+4] += q*z1[i];
     2831    int nUseData = 0;
     2832    for (int k = 0; k < nChan; ++k) {
     2833      if (maskArray[k] == 0) continue;
     2834
     2835      for (int i = 0; i < nDOF; ++i) {
     2836        for (int j = i; j < nDOF; ++j) {
     2837          xMatrix[i][j] += xArray[i][k] * xArray[j][k];
    29132838        }
    2914 
    2915         nUseDataInPiece++;
    2916       }
    2917 
    2918       if (nUseDataInPiece < 1) {
    2919         std::vector<string> suffixOfPieceNumber(4);
    2920         suffixOfPieceNumber[0] = "th";
    2921         suffixOfPieceNumber[1] = "st";
    2922         suffixOfPieceNumber[2] = "nd";
    2923         suffixOfPieceNumber[3] = "rd";
    2924         int idxNoDataPiece = (n % 10 <= 3) ? n : 0;
    2925         ostringstream oss;
    2926         oss << "all channels clipped or masked in " << n << suffixOfPieceNumber[idxNoDataPiece];
    2927         oss << " piece of the spectrum. can't execute fitting anymore.";
    2928         throw(AipsError(String(oss)));
    2929       }
     2839        zMatrix[i] += z1[k] * xArray[i][k];
     2840      }
     2841
     2842      nUseData++;
     2843    }
     2844
     2845    if (nUseData < 1) {
     2846        throw(AipsError("all channels clipped or masked. can't execute fitting anymore."));     
    29302847    }
    29312848
    … …  
    29362853    }
    29372854
    2938     std::vector<double> invDiag(nDOF);
     2855    std::vector<double> invDiag;
    29392856    for (int i = 0; i < nDOF; ++i) {
    2940       invDiag[i] = 1.0/xMatrix[i][i];
     2857      invDiag.push_back(1.0/xMatrix[i][i]);
    29412858      for (int j = 0; j < nDOF; ++j) {
    29422859        xMatrix[i][j] *= invDiag[i];
    … …  
    29672884      }
    29682885    }
    2969     //compute a vector y which consists of the coefficients of the best-fit spline curves
    2970     //(a0,a1,a2,a3(,b3,c3,...)), namely, the ones for the leftmost piece and the ones of
    2971     //cubic terms for the other pieces (in case nPiece>1).
    2972     std::vector<double> y(nDOF);
     2886    //compute a vector y which consists of the coefficients of the sinusoids forming the
     2887    //best-fit curves (a0,s1,c1,s2,c2,...), where a0 is constant and s* and c* are of sine
     2888    //and cosine functions, respectively.
     2889    std::vector<double> y;
     2890    params.clear();
    29732891    for (int i = 0; i < nDOF; ++i) {
    2974       y[i] = 0.0;
     2892      y.push_back(0.0);
    29752893      for (int j = 0; j < nDOF; ++j) {
    29762894        y[i] += xMatrix[i][nDOF+j]*zMatrix[j];
    29772895      }
    2978     }
    2979 
    2980     double a0 = y[0];
    2981     double a1 = y[1];
    2982     double a2 = y[2];
    2983     double a3 = y[3];
    2984 
    2985     int j = 0;
    2986     for (int n = 0; n < nPiece; ++n) {
    2987       for (int i = idxEdge[n]; i < idxEdge[n+1]; ++i) {
    2988         r1[i] = a0 + a1*x1[i] + a2*x2[i] + a3*x3[i];
    2989       }
    2990       params[j]   = a0;
    2991       params[j+1] = a1;
    2992       params[j+2] = a2;
    2993       params[j+3] = a3;
    2994       j += 4;
    2995 
    2996       if (n == nPiece-1) break;
    2997 
    2998       double d = y[4+n];
    2999       double iE = invEdge[n];
    3000       a0 +=     d;
    3001       a1 -= 3.0*d*iE;
    3002       a2 += 3.0*d*iE*iE;
    3003       a3 -=     d*iE*iE*iE;
    3004     }
    3005 
    3006     //subtract constant value for masked regions at the edge of spectrum
    3007     if (idxEdge[0] > 0) {
    3008       int n = idxEdge[0];
    3009       for (int i = 0; i < idxEdge[0]; ++i) {
    3010         //--cubic extrapolate--
    3011         //r1[i] = params[0] + params[1]*x1[i] + params[2]*x2[i] + params[3]*x3[i];
    3012         //--linear extrapolate--
    3013         //r1[i] = (r1[n+1] - r1[n])/(x1[n+1] - x1[n])*(x1[i] - x1[n]) + r1[n];
    3014         //--constant--
    3015         r1[i] = r1[n];
    3016       }
    3017     }
    3018     if (idxEdge[nPiece] < nChan) {
    3019       int n = idxEdge[nPiece]-1;
    3020       for (int i = idxEdge[nPiece]; i < nChan; ++i) {
    3021         //--cubic extrapolate--
    3022         //int m = 4*(nPiece-1);
    3023         //r1[i] = params[m] + params[m+1]*x1[i] + params[m+2]*x2[i] + params[m+3]*x3[i];
    3024         //--linear extrapolate--
    3025         //r1[i] = (r1[n-1] - r1[n])/(x1[n-1] - x1[n])*(x1[i] - x1[n]) + r1[n];
    3026         //--constant--
    3027         r1[i] = r1[n];
    3028       }
     2896      params.push_back(y[i]);
    30292897    }
    30302898
    30312899    for (int i = 0; i < nChan; ++i) {
     2900      r1[i] = y[0];
     2901      for (int j = 1; j < nDOF; ++j) {
     2902        r1[i] += y[j]*xArray[j][i];
     2903      }
    30322904      residual[i] = z1[i] - r1[i];
    30332905    }
    … …  
    30622934  nClipped = initNData - nData;
    30632935
    3064   std::vector<float> result(nChan);
     2936  std::vector<float> result;
    30652937  if (getResidual) {
    30662938    for (int i = 0; i < nChan; ++i) {
    3067       result[i] = (float)residual[i];
     2939      result.push_back((float)residual[i]);
    30682940    }
    30692941  } else {
    30702942    for (int i = 0; i < nChan; ++i) {
    3071       result[i] = (float)r1[i];
     2943      result.push_back((float)r1[i]);
    30722944    }
    30732945  }
    … …  
    30762948}
    30772949
    3078 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)
    3079 {
    3080   nWaves.clear();
    3081 
    3082   if (applyFFT) {
    3083     string fftThAttr;
    3084     float fftThSigma;
    3085     int fftThTop;
    3086     parseThresholdExpression(fftThresh, fftThAttr, fftThSigma, fftThTop);
    3087     doSelectWaveNumbers(whichrow, chanMask, fftMethod, fftThSigma, fftThTop, fftThAttr, nWaves);
    3088   }
    3089 
    3090   addAuxWaveNumbers(whichrow, addNWaves, rejectNWaves, nWaves);
    3091 }
    3092 
    3093 void Scantable::parseThresholdExpression(const std::string& fftThresh, std::string& fftThAttr, float& fftThSigma, int& fftThTop)
    3094 {
    3095   uInt idxSigma = fftThresh.find("sigma");
    3096   uInt idxTop   = fftThresh.find("top");
    3097 
    3098   if (idxSigma == fftThresh.size() - 5) {
    3099     std::istringstream is(fftThresh.substr(0, fftThresh.size() - 5));
    3100     is >> fftThSigma;
    3101     fftThAttr = "sigma";
    3102   } else if (idxTop == 0) {
    3103     std::istringstream is(fftThresh.substr(3));
    3104     is >> fftThTop;
    3105     fftThAttr = "top";
    3106   } else {
    3107     bool isNumber = true;
    3108     for (uInt i = 0; i < fftThresh.size()-1; ++i) {
    3109       char ch = (fftThresh.substr(i, 1).c_str())[0];
    3110       if (!(isdigit(ch) || (fftThresh.substr(i, 1) == "."))) {
    3111         isNumber = false;
    3112         break;
    3113       }
    3114     }
    3115     if (isNumber) {
    3116       std::istringstream is(fftThresh);
    3117       is >> fftThSigma;
    3118       fftThAttr = "sigma";
    3119     } else {
    3120       throw(AipsError("fftthresh has a wrong value"));
    3121     }
    3122   }
    3123 }
    3124 
    3125 void Scantable::doSelectWaveNumbers(const int whichrow, const std::vector<bool>& chanMask, const std::string& fftMethod, const float fftThSigma, const int fftThTop, const std::string& fftThAttr, std::vector<int>& nWaves)
    3126 {
    3127   std::vector<float> fspec;
    3128   if (fftMethod == "fft") {
    3129     fspec = execFFT(whichrow, chanMask, false, true);
    3130   //} else if (fftMethod == "lsp") {
    3131   //  fspec = lombScarglePeriodogram(whichrow);
    3132   }
    3133 
    3134   if (fftThAttr == "sigma") {
    3135     float mean  = 0.0;
    3136     float mean2 = 0.0;
    3137     for (uInt i = 0; i < fspec.size(); ++i) {
    3138       mean  += fspec[i];
    3139       mean2 += fspec[i]*fspec[i];
    3140     }
    3141     mean  /= float(fspec.size());
    3142     mean2 /= float(fspec.size());
    3143     float thres = mean + fftThSigma * float(sqrt(mean2 - mean*mean));
    3144 
    3145     for (uInt i = 0; i < fspec.size(); ++i) {
    3146       if (fspec[i] >= thres) {
    3147         nWaves.push_back(i);
    3148       }
    3149     }
    3150 
    3151   } else if (fftThAttr == "top") {
    3152     for (int i = 0; i < fftThTop; ++i) {
    3153       float max = 0.0;
    3154       int maxIdx = 0;
    3155       for (uInt j = 0; j < fspec.size(); ++j) {
    3156         if (fspec[j] > max) {
    3157           max = fspec[j];
    3158           maxIdx = j;
    3159         }
    3160       }
    3161       nWaves.push_back(maxIdx);
    3162       fspec[maxIdx] = 0.0;
    3163     }
    3164 
    3165   }
    3166 
    3167   if (nWaves.size() > 1) {
    3168     sort(nWaves.begin(), nWaves.end());
    3169   }
    3170 }
    3171 
    3172 void Scantable::addAuxWaveNumbers(const int whichrow, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, std::vector<int>& nWaves)
    3173 {
    3174   std::vector<int> tempAddNWaves, tempRejectNWaves;
    3175   for (uInt i = 0; i < addNWaves.size(); ++i) {
    3176     tempAddNWaves.push_back(addNWaves[i]);
    3177   }
    3178   if ((tempAddNWaves.size() == 2) && (tempAddNWaves[1] == -999)) {
    3179     setWaveNumberListUptoNyquistFreq(whichrow, tempAddNWaves);
    3180   }
    3181 
    3182   for (uInt i = 0; i < rejectNWaves.size(); ++i) {
    3183     tempRejectNWaves.push_back(rejectNWaves[i]);
    3184   }
    3185   if ((tempRejectNWaves.size() == 2) && (tempRejectNWaves[1] == -999)) {
    3186     setWaveNumberListUptoNyquistFreq(whichrow, tempRejectNWaves);
    3187   }
    3188 
    3189   for (uInt i = 0; i < tempAddNWaves.size(); ++i) {
    3190     bool found = false;
    3191     for (uInt j = 0; j < nWaves.size(); ++j) {
    3192       if (nWaves[j] == tempAddNWaves[i]) {
    3193         found = true;
    3194         break;
    3195       }
    3196     }
    3197     if (!found) nWaves.push_back(tempAddNWaves[i]);
    3198   }
    3199 
    3200   for (uInt i = 0; i < tempRejectNWaves.size(); ++i) {
    3201     for (std::vector<int>::iterator j = nWaves.begin(); j != nWaves.end(); ) {
    3202       if (*j == tempRejectNWaves[i]) {
    3203         j = nWaves.erase(j);
    3204       } else {
    3205         ++j;
    3206       }
    3207     }
    3208   }
    3209 
    3210   if (nWaves.size() > 1) {
    3211     sort(nWaves.begin(), nWaves.end());
    3212     unique(nWaves.begin(), nWaves.end());
    3213   }
    3214 }
    3215 
    3216 void Scantable::setWaveNumberListUptoNyquistFreq(const int whichrow, std::vector<int>& nWaves)
    3217 {
    3218   if ((nWaves.size() == 2)&&(nWaves[1] == -999)) {
    3219     int val = nWaves[0];
    3220     int nyquistFreq = nchan(getIF(whichrow))/2+1;
    3221     nWaves.clear();
    3222     if (val > nyquistFreq) {  // for safety, at least nWaves contains a constant; CAS-3759
    3223       nWaves.push_back(0);
    3224     }
    3225     while (val <= nyquistFreq) {
    3226       nWaves.push_back(val);
    3227       val++;
    3228     }
    3229   }
    3230 }
    3231 
    3232 void Scantable::sinusoidBaseline(const std::vector<bool>& mask, const bool applyFFT, const std::string& fftMethod, const std::string& fftThresh, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, float thresClip, int nIterClip, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
    3233 {
     2950void 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)
     2951{
     2952  /*************
     2953  double totTimeStart, totTimeEnd, blTimeStart, blTimeEnd, ioTimeStart, ioTimeEnd, msTimeStart, msTimeEnd, seTimeStart, seTimeEnd, otTimeStart, otTimeEnd, prTimeStart, prTimeEnd;
     2954  double elapseMs = 0.0;
     2955  double elapseSe = 0.0;
     2956  double elapseOt = 0.0;
     2957  double elapsePr = 0.0;
     2958  double elapseBl = 0.0;
     2959  double elapseIo = 0.0;
     2960  totTimeStart = mathutil::gettimeofday_sec();
     2961  *************/
     2962
    32342963  try {
    32352964    ofstream ofs;
    … …  
    32522981
    32532982    //Fitter fitter = Fitter();
    3254     //fitter.setExpression("sinusoid", nWaves);
     2983    //fitter.setExpression("cspline", nPiece);
    32552984    //fitter.setIterClipping(thresClip, nIterClip);
    3256 
    3257     int nRow = nrow();
    3258     std::vector<bool> chanMask;
    3259     std::vector<int> nWaves;
    32602985
    32612986    bool showProgress;
    … …  
    32632988    parseProgressInfo(progressInfo, showProgress, minNRow);
    32642989
     2990    int nRow = nrow();
     2991    std::vector<bool> chanMask;
     2992
     2993    //--------------------------------
    32652994    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
     2995      /******************
     2996      ioTimeStart = mathutil::gettimeofday_sec();
     2997      **/
     2998      std::vector<float> sp = getSpectrum(whichrow);
     2999      /**
     3000      ioTimeEnd = mathutil::gettimeofday_sec();
     3001      elapseIo += (double)(ioTimeEnd - ioTimeStart);
     3002      msTimeStart = mathutil::gettimeofday_sec();
     3003      ******************/
     3004
    32663005      chanMask = getCompositeChanMask(whichrow, mask);
    3267       selectWaveNumbers(whichrow, chanMask, applyFFT, fftMethod, fftThresh, addNWaves, rejectNWaves, nWaves);
    3268 
    3269       //FOR DEBUGGING------------
    3270       /*
    3271       if (whichrow < 0) {// == nRow -1) {
    3272         cout << "+++ i=" << setw(3) << whichrow << ", IF=" << setw(2) << getIF(whichrow);
    3273         if (applyFFT) {
    3274           cout << "[ ";
    3275           for (uInt j = 0; j < nWaves.size(); ++j) {
    3276             cout << nWaves[j] << ", ";
    3277           }
    3278           cout << " ]    " << endl;
    3279         }
    3280         cout << flush;
    3281       }
    3282       */
    3283       //-------------------------
     3006
     3007      /**
     3008      msTimeEnd = mathutil::gettimeofday_sec();
     3009      elapseMs += (double)(msTimeEnd - msTimeStart);
     3010      blTimeStart = mathutil::gettimeofday_sec();
     3011      **/
    32843012
    32853013      //fitBaseline(chanMask, whichrow, fitter);
    32863014      //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
    3287       std::vector<float> params;
     3015      std::vector<int> pieceEdges(nPiece+1);
     3016      std::vector<float> params(nPiece*4);
    32883017      int nClipped = 0;
    3289       std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, params, nClipped, thresClip, nIterClip, getResidual);
     3018      std::vector<float> res = doCubicSplineFitting(sp, chanMask, nPiece, pieceEdges, params, nClipped, thresClip, nIterClip, getResidual);
     3019
     3020      /**
     3021      blTimeEnd = mathutil::gettimeofday_sec();
     3022      elapseBl += (double)(blTimeEnd - blTimeStart);
     3023      seTimeStart = mathutil::gettimeofday_sec();
     3024      **/
     3025
     3026
    32903027      setSpectrum(res, whichrow);
    32913028      //
    32923029
    3293       outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "sinusoidBaseline()", params, nClipped);
     3030      /**
     3031      seTimeEnd = mathutil::gettimeofday_sec();
     3032      elapseSe += (double)(seTimeEnd - seTimeStart);
     3033      otTimeStart = mathutil::gettimeofday_sec();
     3034      **/
     3035
     3036      outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()", pieceEdges, params, nClipped);
     3037
     3038      /**
     3039      otTimeEnd = mathutil::gettimeofday_sec();
     3040      elapseOt += (double)(otTimeEnd - otTimeStart);
     3041      prTimeStart = mathutil::gettimeofday_sec();
     3042      **/
     3043
    32943044      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
    3295     }
    3296 
     3045
     3046      /******************
     3047      prTimeEnd = mathutil::gettimeofday_sec();
     3048      elapsePr += (double)(prTimeEnd - prTimeStart);
     3049      ******************/
     3050    }
     3051    //--------------------------------
     3052   
    32973053    if (outTextFile) ofs.close();
    32983054
    … …  
    33003056    throw;
    33013057  }
    3302 }
    3303 
    3304 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)
     3058  /***************
     3059  totTimeEnd = mathutil::gettimeofday_sec();
     3060  std::cout << "io    : " << elapseIo << " (sec.)" << endl;
     3061  std::cout << "ms    : " << elapseMs << " (sec.)" << endl;
     3062  std::cout << "bl    : " << elapseBl << " (sec.)" << endl;
     3063  std::cout << "se    : " << elapseSe << " (sec.)" << endl;
     3064  std::cout << "ot    : " << elapseOt << " (sec.)" << endl;
     3065  std::cout << "pr    : " << elapsePr << " (sec.)" << endl;
     3066  std::cout << "total : " << (double)(totTimeEnd - totTimeStart) << " (sec.)" << endl;
     3067  ***************/
     3068}
     3069
     3070void 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)
    33053071{
    33063072  try {
    … …  
    33243090
    33253091    //Fitter fitter = Fitter();
    3326     //fitter.setExpression("sinusoid", nWaves);
     3092    //fitter.setExpression("cspline", nPiece);
    33273093    //fitter.setIterClipping(thresClip, nIterClip);
    33283094
    33293095    int nRow = nrow();
    33303096    std::vector<bool> chanMask;
    3331     std::vector<int> nWaves;
    3332 
    33333097    int minEdgeSize = getIFNos().size()*2;
    33343098    STLineFinder lineFinder = STLineFinder();
    … …  
    33403104
    33413105    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
     3106      std::vector<float> sp = getSpectrum(whichrow);
    33423107
    33433108      //-------------------------------------------------------
    … …  
    33593124        throw(AipsError("Wrong length of edge element"));
    33603125      }
     3126      //lineFinder.setData(getSpectrum(whichrow));
     3127      lineFinder.setData(sp);
     3128      lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
     3129      chanMask = lineFinder.getMask();
     3130      //-------------------------------------------------------
     3131
     3132
     3133      //fitBaseline(chanMask, whichrow, fitter);
     3134      //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
     3135      std::vector<int> pieceEdges(nPiece+1);
     3136      std::vector<float> params(nPiece*4);
     3137      int nClipped = 0;
     3138      std::vector<float> res = doCubicSplineFitting(sp, chanMask, nPiece, pieceEdges, params, nClipped, thresClip, nIterClip, getResidual);
     3139      setSpectrum(res, whichrow);
     3140      //
     3141
     3142      outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()", pieceEdges, params, nClipped);
     3143      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
     3144    }
     3145
     3146    if (outTextFile) ofs.close();
     3147
     3148  } catch (...) {
     3149    throw;
     3150  }
     3151}
     3152
     3153std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data, const std::vector<bool>& mask, int nPiece, std::vector<int>& idxEdge, std::vector<float>& params, int& nClipped, float thresClip, int nIterClip, bool getResidual)
     3154{
     3155  if (data.size() != mask.size()) {
     3156    throw(AipsError("data and mask sizes are not identical"));
     3157  }
     3158  if (nPiece < 1) {
     3159    throw(AipsError("number of the sections must be one or more"));
     3160  }
     3161
     3162  int nChan = data.size();
     3163  std::vector<int> maskArray(nChan);
     3164  std::vector<int> x(nChan);
     3165  int j = 0;
     3166  for (int i = 0; i < nChan; ++i) {
     3167    maskArray[i] = mask[i] ? 1 : 0;
     3168    if (mask[i]) {
     3169      x[j] = i;
     3170      j++;
     3171    }
     3172  }
     3173  int initNData = j;
     3174
     3175  if (initNData < nPiece) {
     3176    throw(AipsError("too few non-flagged channels"));
     3177  }
     3178
     3179  int nElement = (int)(floor(floor((double)(initNData/nPiece))+0.5));
     3180  std::vector<double> invEdge(nPiece-1);
     3181  idxEdge[0] = x[0];
     3182  for (int i = 1; i < nPiece; ++i) {
     3183    int valX = x[nElement*i];
     3184    idxEdge[i] = valX;
     3185    invEdge[i-1] = 1.0/(double)valX;
     3186  }
     3187  idxEdge[nPiece] = x[initNData-1]+1;
     3188
     3189  int nData = initNData;
     3190  int nDOF = nPiece + 3;  //number of parameters to solve, namely, 4+(nPiece-1).
     3191
     3192  std::vector<double> x1(nChan), x2(nChan), x3(nChan);
     3193  std::vector<double> z1(nChan), x1z1(nChan), x2z1(nChan), x3z1(nChan);
     3194  std::vector<double> r1(nChan), residual(nChan);
     3195  for (int i = 0; i < nChan; ++i) {
     3196    double di = (double)i;
     3197    double dD = (double)data[i];
     3198    x1[i]   = di;
     3199    x2[i]   = di*di;
     3200    x3[i]   = di*di*di;
     3201    z1[i]   = dD;
     3202    x1z1[i] = dD*di;
     3203    x2z1[i] = dD*di*di;
     3204    x3z1[i] = dD*di*di*di;
     3205    r1[i]   = 0.0;
     3206    residual[i] = 0.0;
     3207  }
     3208
     3209  for (int nClip = 0; nClip < nIterClip+1; ++nClip) {
     3210    // xMatrix : horizontal concatenation of
     3211    //           the least-sq. matrix (left) and an
     3212    //           identity matrix (right).
     3213    // the right part is used to calculate the inverse matrix of the left part.
     3214    double xMatrix[nDOF][2*nDOF];
     3215    double zMatrix[nDOF];
     3216    for (int i = 0; i < nDOF; ++i) {
     3217      for (int j = 0; j < 2*nDOF; ++j) {
     3218        xMatrix[i][j] = 0.0;
     3219      }
     3220      xMatrix[i][nDOF+i] = 1.0;
     3221      zMatrix[i] = 0.0;
     3222    }
     3223
     3224    for (int n = 0; n < nPiece; ++n) {
     3225      int nUseDataInPiece = 0;
     3226      for (int i = idxEdge[n]; i < idxEdge[n+1]; ++i) {
     3227
     3228        if (maskArray[i] == 0) continue;
     3229
     3230        xMatrix[0][0] += 1.0;
     3231        xMatrix[0][1] += x1[i];
     3232        xMatrix[0][2] += x2[i];
     3233        xMatrix[0][3] += x3[i];
     3234        xMatrix[1][1] += x2[i];
     3235        xMatrix[1][2] += x3[i];
     3236        xMatrix[1][3] += x2[i]*x2[i];
     3237        xMatrix[2][2] += x2[i]*x2[i];
     3238        xMatrix[2][3] += x3[i]*x2[i];
     3239        xMatrix[3][3] += x3[i]*x3[i];
     3240        zMatrix[0] += z1[i];
     3241        zMatrix[1] += x1z1[i];
     3242        zMatrix[2] += x2z1[i];
     3243        zMatrix[3] += x3z1[i];
     3244
     3245        for (int j = 0; j < n; ++j) {
     3246          double q = 1.0 - x1[i]*invEdge[j];
     3247          q = q*q*q;
     3248          xMatrix[0][j+4] += q;
     3249          xMatrix[1][j+4] += q*x1[i];
     3250          xMatrix[2][j+4] += q*x2[i];
     3251          xMatrix[3][j+4] += q*x3[i];
     3252          for (int k = 0; k < j; ++k) {
     3253            double r = 1.0 - x1[i]*invEdge[k];
     3254            r = r*r*r;
     3255            xMatrix[k+4][j+4] += r*q;
     3256          }
     3257          xMatrix[j+4][j+4] += q*q;
     3258          zMatrix[j+4] += q*z1[i];
     3259        }
     3260
     3261        nUseDataInPiece++;
     3262      }
     3263
     3264      if (nUseDataInPiece < 1) {
     3265        std::vector<string> suffixOfPieceNumber(4);
     3266        suffixOfPieceNumber[0] = "th";
     3267        suffixOfPieceNumber[1] = "st";
     3268        suffixOfPieceNumber[2] = "nd";
     3269        suffixOfPieceNumber[3] = "rd";
     3270        int idxNoDataPiece = (n % 10 <= 3) ? n : 0;
     3271        ostringstream oss;
     3272        oss << "all channels clipped or masked in " << n << suffixOfPieceNumber[idxNoDataPiece];
     3273        oss << " piece of the spectrum. can't execute fitting anymore.";
     3274        throw(AipsError(String(oss)));
     3275      }
     3276    }
     3277
     3278    for (int i = 0; i < nDOF; ++i) {
     3279      for (int j = 0; j < i; ++j) {
     3280        xMatrix[i][j] = xMatrix[j][i];
     3281      }
     3282    }
     3283
     3284    std::vector<double> invDiag(nDOF);
     3285    for (int i = 0; i < nDOF; ++i) {
     3286      invDiag[i] = 1.0/xMatrix[i][i];
     3287      for (int j = 0; j < nDOF; ++j) {
     3288        xMatrix[i][j] *= invDiag[i];
     3289      }
     3290    }
     3291
     3292    for (int k = 0; k < nDOF; ++k) {
     3293      for (int i = 0; i < nDOF; ++i) {
     3294        if (i != k) {
     3295          double factor1 = xMatrix[k][k];
     3296          double factor2 = xMatrix[i][k];
     3297          for (int j = k; j < 2*nDOF; ++j) {
     3298            xMatrix[i][j] *= factor1;
     3299            xMatrix[i][j] -= xMatrix[k][j]*factor2;
     3300            xMatrix[i][j] /= factor1;
     3301          }
     3302        }
     3303      }
     3304      double xDiag = xMatrix[k][k];
     3305      for (int j = k; j < 2*nDOF; ++j) {
     3306        xMatrix[k][j] /= xDiag;
     3307      }
     3308    }
     3309   
     3310    for (int i = 0; i < nDOF; ++i) {
     3311      for (int j = 0; j < nDOF; ++j) {
     3312        xMatrix[i][nDOF+j] *= invDiag[j];
     3313      }
     3314    }
     3315    //compute a vector y which consists of the coefficients of the best-fit spline curves
     3316    //(a0,a1,a2,a3(,b3,c3,...)), namely, the ones for the leftmost piece and the ones of
     3317    //cubic terms for the other pieces (in case nPiece>1).
     3318    std::vector<double> y(nDOF);
     3319    for (int i = 0; i < nDOF; ++i) {
     3320      y[i] = 0.0;
     3321      for (int j = 0; j < nDOF; ++j) {
     3322        y[i] += xMatrix[i][nDOF+j]*zMatrix[j];
     3323      }
     3324    }
     3325
     3326    double a0 = y[0];
     3327    double a1 = y[1];
     3328    double a2 = y[2];
     3329    double a3 = y[3];
     3330
     3331    int j = 0;
     3332    for (int n = 0; n < nPiece; ++n) {
     3333      for (int i = idxEdge[n]; i < idxEdge[n+1]; ++i) {
     3334        r1[i] = a0 + a1*x1[i] + a2*x2[i] + a3*x3[i];
     3335      }
     3336      params[j]   = a0;
     3337      params[j+1] = a1;
     3338      params[j+2] = a2;
     3339      params[j+3] = a3;
     3340      j += 4;
     3341
     3342      if (n == nPiece-1) break;
     3343
     3344      double d = y[4+n];
     3345      double iE = invEdge[n];
     3346      a0 +=     d;
     3347      a1 -= 3.0*d*iE;
     3348      a2 += 3.0*d*iE*iE;
     3349      a3 -=     d*iE*iE*iE;
     3350    }
     3351
     3352    //subtract constant value for masked regions at the edge of spectrum
     3353    if (idxEdge[0] > 0) {
     3354      int n = idxEdge[0];
     3355      for (int i = 0; i < idxEdge[0]; ++i) {
     3356        //--cubic extrapolate--
     3357        //r1[i] = params[0] + params[1]*x1[i] + params[2]*x2[i] + params[3]*x3[i];
     3358        //--linear extrapolate--
     3359        //r1[i] = (r1[n+1] - r1[n])/(x1[n+1] - x1[n])*(x1[i] - x1[n]) + r1[n];
     3360        //--constant--
     3361        r1[i] = r1[n];
     3362      }
     3363    }
     3364    if (idxEdge[nPiece] < nChan) {
     3365      int n = idxEdge[nPiece]-1;
     3366      for (int i = idxEdge[nPiece]; i < nChan; ++i) {
     3367        //--cubic extrapolate--
     3368        //int m = 4*(nPiece-1);
     3369        //r1[i] = params[m] + params[m+1]*x1[i] + params[m+2]*x2[i] + params[m+3]*x3[i];
     3370        //--linear extrapolate--
     3371        //r1[i] = (r1[n-1] - r1[n])/(x1[n-1] - x1[n])*(x1[i] - x1[n]) + r1[n];
     3372        //--constant--
     3373        r1[i] = r1[n];
     3374      }
     3375    }
     3376
     3377    for (int i = 0; i < nChan; ++i) {
     3378      residual[i] = z1[i] - r1[i];
     3379    }
     3380
     3381    if ((nClip == nIterClip) || (thresClip <= 0.0)) {
     3382      break;
     3383    } else {
     3384      double stdDev = 0.0;
     3385      for (int i = 0; i < nChan; ++i) {
     3386        stdDev += residual[i]*residual[i]*(double)maskArray[i];
     3387      }
     3388      stdDev = sqrt(stdDev/(double)nData);
     3389     
     3390      double thres = stdDev * thresClip;
     3391      int newNData = 0;
     3392      for (int i = 0; i < nChan; ++i) {
     3393        if (abs(residual[i]) >= thres) {
     3394          maskArray[i] = 0;
     3395        }
     3396        if (maskArray[i] > 0) {
     3397          newNData++;
     3398        }
     3399      }
     3400      if (newNData == nData) {
     3401        break; //no more flag to add. iteration stops.
     3402      } else {
     3403        nData = newNData;
     3404      }
     3405    }
     3406  }
     3407
     3408  nClipped = initNData - nData;
     3409
     3410  std::vector<float> result(nChan);
     3411  if (getResidual) {
     3412    for (int i = 0; i < nChan; ++i) {
     3413      result[i] = (float)residual[i];
     3414    }
     3415  } else {
     3416    for (int i = 0; i < nChan; ++i) {
     3417      result[i] = (float)r1[i];
     3418    }
     3419  }
     3420
     3421  return result;
     3422}
     3423
     3424void 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)
     3425{
     3426  nWaves.clear();
     3427
     3428  if (applyFFT) {
     3429    string fftThAttr;
     3430    float fftThSigma;
     3431    int fftThTop;
     3432    parseThresholdExpression(fftThresh, fftThAttr, fftThSigma, fftThTop);
     3433    doSelectWaveNumbers(whichrow, chanMask, fftMethod, fftThSigma, fftThTop, fftThAttr, nWaves);
     3434  }
     3435
     3436  addAuxWaveNumbers(whichrow, addNWaves, rejectNWaves, nWaves);
     3437}
     3438
     3439void Scantable::parseThresholdExpression(const std::string& fftThresh, std::string& fftThAttr, float& fftThSigma, int& fftThTop)
     3440{
     3441  uInt idxSigma = fftThresh.find("sigma");
     3442  uInt idxTop   = fftThresh.find("top");
     3443
     3444  if (idxSigma == fftThresh.size() - 5) {
     3445    std::istringstream is(fftThresh.substr(0, fftThresh.size() - 5));
     3446    is >> fftThSigma;
     3447    fftThAttr = "sigma";
     3448  } else if (idxTop == 0) {
     3449    std::istringstream is(fftThresh.substr(3));
     3450    is >> fftThTop;
     3451    fftThAttr = "top";
     3452  } else {
     3453    bool isNumber = true;
     3454    for (uInt i = 0; i < fftThresh.size()-1; ++i) {
     3455      char ch = (fftThresh.substr(i, 1).c_str())[0];
     3456      if (!(isdigit(ch) || (fftThresh.substr(i, 1) == "."))) {
     3457        isNumber = false;
     3458        break;
     3459      }
     3460    }
     3461    if (isNumber) {
     3462      std::istringstream is(fftThresh);
     3463      is >> fftThSigma;
     3464      fftThAttr = "sigma";
     3465    } else {
     3466      throw(AipsError("fftthresh has a wrong value"));
     3467    }
     3468  }
     3469}
     3470
     3471void Scantable::doSelectWaveNumbers(const int whichrow, const std::vector<bool>& chanMask, const std::string& fftMethod, const float fftThSigma, const int fftThTop, const std::string& fftThAttr, std::vector<int>& nWaves)
     3472{
     3473  std::vector<float> fspec;
     3474  if (fftMethod == "fft") {
     3475    fspec = execFFT(whichrow, chanMask, false, true);
     3476  //} else if (fftMethod == "lsp") {
     3477  //  fspec = lombScarglePeriodogram(whichrow);
     3478  }
     3479
     3480  if (fftThAttr == "sigma") {
     3481    float mean  = 0.0;
     3482    float mean2 = 0.0;
     3483    for (uInt i = 0; i < fspec.size(); ++i) {
     3484      mean  += fspec[i];
     3485      mean2 += fspec[i]*fspec[i];
     3486    }
     3487    mean  /= float(fspec.size());
     3488    mean2 /= float(fspec.size());
     3489    float thres = mean + fftThSigma * float(sqrt(mean2 - mean*mean));
     3490
     3491    for (uInt i = 0; i < fspec.size(); ++i) {
     3492      if (fspec[i] >= thres) {
     3493        nWaves.push_back(i);
     3494      }
     3495    }
     3496
     3497  } else if (fftThAttr == "top") {
     3498    for (int i = 0; i < fftThTop; ++i) {
     3499      float max = 0.0;
     3500      int maxIdx = 0;
     3501      for (uInt j = 0; j < fspec.size(); ++j) {
     3502        if (fspec[j] > max) {
     3503          max = fspec[j];
     3504          maxIdx = j;
     3505        }
     3506      }
     3507      nWaves.push_back(maxIdx);
     3508      fspec[maxIdx] = 0.0;
     3509    }
     3510
     3511  }
     3512
     3513  if (nWaves.size() > 1) {
     3514    sort(nWaves.begin(), nWaves.end());
     3515  }
     3516}
     3517
     3518void Scantable::addAuxWaveNumbers(const int whichrow, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, std::vector<int>& nWaves)
     3519{
     3520  std::vector<int> tempAddNWaves, tempRejectNWaves;
     3521  for (uInt i = 0; i < addNWaves.size(); ++i) {
     3522    tempAddNWaves.push_back(addNWaves[i]);
     3523  }
     3524  if ((tempAddNWaves.size() == 2) && (tempAddNWaves[1] == -999)) {
     3525    setWaveNumberListUptoNyquistFreq(whichrow, tempAddNWaves);
     3526  }
     3527
     3528  for (uInt i = 0; i < rejectNWaves.size(); ++i) {
     3529    tempRejectNWaves.push_back(rejectNWaves[i]);
     3530  }
     3531  if ((tempRejectNWaves.size() == 2) && (tempRejectNWaves[1] == -999)) {
     3532    setWaveNumberListUptoNyquistFreq(whichrow, tempRejectNWaves);
     3533  }
     3534
     3535  for (uInt i = 0; i < tempAddNWaves.size(); ++i) {
     3536    bool found = false;
     3537    for (uInt j = 0; j < nWaves.size(); ++j) {
     3538      if (nWaves[j] == tempAddNWaves[i]) {
     3539        found = true;
     3540        break;
     3541      }
     3542    }
     3543    if (!found) nWaves.push_back(tempAddNWaves[i]);
     3544  }
     3545
     3546  for (uInt i = 0; i < tempRejectNWaves.size(); ++i) {
     3547    for (std::vector<int>::iterator j = nWaves.begin(); j != nWaves.end(); ) {
     3548      if (*j == tempRejectNWaves[i]) {
     3549        j = nWaves.erase(j);
     3550      } else {
     3551        ++j;
     3552      }
     3553    }
     3554  }
     3555
     3556  if (nWaves.size() > 1) {
     3557    sort(nWaves.begin(), nWaves.end());
     3558    unique(nWaves.begin(), nWaves.end());
     3559  }
     3560}
     3561
     3562void Scantable::setWaveNumberListUptoNyquistFreq(const int whichrow, std::vector<int>& nWaves)
     3563{
     3564  if ((nWaves.size() == 2)&&(nWaves[1] == -999)) {
     3565    int val = nWaves[0];
     3566    int nyquistFreq = nchan(getIF(whichrow))/2+1;
     3567    nWaves.clear();
     3568    if (val > nyquistFreq) {  // for safety, at least nWaves contains a constant; CAS-3759
     3569      nWaves.push_back(0);
     3570    }
     3571    while (val <= nyquistFreq) {
     3572      nWaves.push_back(val);
     3573      val++;
     3574    }
     3575  }
     3576}
     3577
     3578void Scantable::sinusoidBaseline(const std::vector<bool>& mask, const bool applyFFT, const std::string& fftMethod, const std::string& fftThresh, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, float thresClip, int nIterClip, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
     3579{
     3580  try {
     3581    ofstream ofs;
     3582    String coordInfo = "";
     3583    bool hasSameNchan = true;
     3584    bool outTextFile = false;
     3585    bool csvFormat = false;
     3586
     3587    if (blfile != "") {
     3588      csvFormat = (blfile.substr(0, 1) == "T");
     3589      ofs.open(blfile.substr(1).c_str(), ios::out | ios::app);
     3590      if (ofs) outTextFile = true;
     3591    }
     3592
     3593    if (outLogger || outTextFile) {
     3594      coordInfo = getCoordInfo()[0];
     3595      if (coordInfo == "") coordInfo = "channel";
     3596      hasSameNchan = hasSameNchanOverIFs();
     3597    }
     3598
     3599    //Fitter fitter = Fitter();
     3600    //fitter.setExpression("sinusoid", nWaves);
     3601    //fitter.setIterClipping(thresClip, nIterClip);
     3602
     3603    int nRow = nrow();
     3604    std::vector<bool> chanMask;
     3605    std::vector<int> nWaves;
     3606
     3607    bool showProgress;
     3608    int minNRow;
     3609    parseProgressInfo(progressInfo, showProgress, minNRow);
     3610
     3611    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
     3612      chanMask = getCompositeChanMask(whichrow, mask);
     3613      selectWaveNumbers(whichrow, chanMask, applyFFT, fftMethod, fftThresh, addNWaves, rejectNWaves, nWaves);
     3614
     3615      //FOR DEBUGGING------------
     3616      /*
     3617      if (whichrow < 0) {// == nRow -1) {
     3618        cout << "+++ i=" << setw(3) << whichrow << ", IF=" << setw(2) << getIF(whichrow);
     3619        if (applyFFT) {
     3620          cout << "[ ";
     3621          for (uInt j = 0; j < nWaves.size(); ++j) {
     3622            cout << nWaves[j] << ", ";
     3623          }
     3624          cout << " ]    " << endl;
     3625        }
     3626        cout << flush;
     3627      }
     3628      */
     3629      //-------------------------
     3630
     3631      //fitBaseline(chanMask, whichrow, fitter);
     3632      //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
     3633      std::vector<float> params;
     3634      int nClipped = 0;
     3635      std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, params, nClipped, thresClip, nIterClip, getResidual);
     3636      setSpectrum(res, whichrow);
     3637      //
     3638
     3639      outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "sinusoidBaseline()", params, nClipped);
     3640      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
     3641    }
     3642
     3643    if (outTextFile) ofs.close();
     3644
     3645  } catch (...) {
     3646    throw;
     3647  }
     3648}
     3649
     3650void 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)
     3651{
     3652  try {
     3653    ofstream ofs;
     3654    String coordInfo = "";
     3655    bool hasSameNchan = true;
     3656    bool outTextFile = false;
     3657    bool csvFormat = false;
     3658
     3659    if (blfile != "") {
     3660      csvFormat = (blfile.substr(0, 1) == "T");
     3661      ofs.open(blfile.substr(1).c_str(), ios::out | ios::app);
     3662      if (ofs) outTextFile = true;
     3663    }
     3664
     3665    if (outLogger || outTextFile) {
     3666      coordInfo = getCoordInfo()[0];
     3667      if (coordInfo == "") coordInfo = "channel";
     3668      hasSameNchan = hasSameNchanOverIFs();
     3669    }
     3670
     3671    //Fitter fitter = Fitter();
     3672    //fitter.setExpression("sinusoid", nWaves);
     3673    //fitter.setIterClipping(thresClip, nIterClip);
     3674
     3675    int nRow = nrow();
     3676    std::vector<bool> chanMask;
     3677    std::vector<int> nWaves;
     3678
     3679    int minEdgeSize = getIFNos().size()*2;
     3680    STLineFinder lineFinder = STLineFinder();
     3681    lineFinder.setOptions(threshold, 3, chanAvgLimit);
     3682
     3683    bool showProgress;
     3684    int minNRow;
     3685    parseProgressInfo(progressInfo, showProgress, minNRow);
     3686
     3687    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
     3688
     3689      //-------------------------------------------------------
     3690      //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
     3691      //-------------------------------------------------------
     3692      int edgeSize = edge.size();
     3693      std::vector<int> currentEdge;
     3694      if (edgeSize >= 2) {
     3695        int idx = 0;
     3696        if (edgeSize > 2) {
     3697          if (edgeSize < minEdgeSize) {
     3698            throw(AipsError("Length of edge element info is less than that of IFs"));
     3699          }
     3700          idx = 2 * getIF(whichrow);
     3701        }
     3702        currentEdge.push_back(edge[idx]);
     3703        currentEdge.push_back(edge[idx+1]);
     3704      } else {
     3705        throw(AipsError("Wrong length of edge element"));
     3706      }
    33613707      lineFinder.setData(getSpectrum(whichrow));
    33623708      lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
    … …  
    37384084}
    37394085
    3740 /* for sinusoid. will be merged once sinusoid is available in fitter (2011/3/10 WK) */
     4086/* for chebyshev/sinusoid. will be merged once chebyshev/sinusoid is available in fitter (2011/3/10 WK) */
    37414087void 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)
    37424088{
Note: See TracChangeset for help on using the changeset viewer.