Changeset 2645 for trunk


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.


Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/python/scantable.py

    r2643 r2645  
    30073007
    30083008    @asaplog_post_dec
     3009    def chebyshev_baseline(self, insitu=None, mask=None, order=None,
     3010                           clipthresh=None, clipniter=None, plot=None,
     3011                           getresidual=None, showprogress=None, minnrow=None,
     3012                           outlog=None, blfile=None, csvformat=None):
     3013        """\
     3014        Return a scan which has been baselined (all rows) by Chebyshev polynomials.
     3015
     3016        Parameters:
     3017            insitu:       If False a new scantable is returned.
     3018                          Otherwise, the scaling is done in-situ
     3019                          The default is taken from .asaprc (False)
     3020            mask:         An optional mask
     3021            order:        the maximum order of Chebyshev polynomial (default is 5)
     3022            clipthresh:   Clipping threshold. (default is 3.0, unit: sigma)
     3023            clipniter:    maximum number of iteration of 'clipthresh'-sigma
     3024                          clipping (default is 0)
     3025            plot:     *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
     3026                          plot the fit and the residual. In this each
     3027                          indivual fit has to be approved, by typing 'y'
     3028                          or 'n'
     3029            getresidual:  if False, returns best-fit values instead of
     3030                          residual. (default is True)
     3031            showprogress: show progress status for large data.
     3032                          default is True.
     3033            minnrow:      minimum number of input spectra to show.
     3034                          default is 1000.
     3035            outlog:       Output the coefficients of the best-fit
     3036                          function to logger (default is False)
     3037            blfile:       Name of a text file in which the best-fit
     3038                          parameter values to be written
     3039                          (default is "": no file/logger output)
     3040            csvformat:    if True blfile is csv-formatted, default is False.
     3041
     3042        Example:
     3043            # return a scan baselined by a cubic spline consisting of 2 pieces
     3044            # (i.e., 1 internal knot),
     3045            # also with 3-sigma clipping, iteration up to 4 times
     3046            bscan = scan.cspline_baseline(npiece=2,clipthresh=3.0,clipniter=4)
     3047       
     3048        Note:
     3049            The best-fit parameter values output in logger and/or blfile are now
     3050            based on specunit of 'channel'.
     3051        """
     3052       
     3053        try:
     3054            varlist = vars()
     3055           
     3056            if insitu is None: insitu = rcParams['insitu']
     3057            if insitu:
     3058                workscan = self
     3059            else:
     3060                workscan = self.copy()
     3061
     3062            #if mask         is None: mask         = [True for i in xrange(workscan.nchan())]
     3063            if mask         is None: mask         = []
     3064            if order        is None: order        = 5
     3065            if clipthresh   is None: clipthresh   = 3.0
     3066            if clipniter    is None: clipniter    = 0
     3067            if plot         is None: plot         = False
     3068            if getresidual  is None: getresidual  = True
     3069            if showprogress is None: showprogress = True
     3070            if minnrow      is None: minnrow      = 1000
     3071            if outlog       is None: outlog       = False
     3072            if blfile       is None: blfile       = ''
     3073            if csvformat     is None: csvformat     = False
     3074
     3075            if csvformat:
     3076                scsvformat = "T"
     3077            else:
     3078                scsvformat = "F"
     3079
     3080            #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method.
     3081            workscan._chebyshev_baseline(mask, order, clipthresh, clipniter,
     3082                                         getresidual,
     3083                                         pack_progress_params(showprogress,
     3084                                                              minnrow),
     3085                                         outlog, scsvformat+blfile)
     3086            workscan._add_history("chebyshev_baseline", varlist)
     3087           
     3088            if insitu:
     3089                self._assign(workscan)
     3090            else:
     3091                return workscan
     3092           
     3093        except RuntimeError, e:
     3094            raise_fitting_failure_exception(e)
     3095
     3096    @asaplog_post_dec
     3097    def auto_chebyshev_baseline(self, insitu=None, mask=None, order=None,
     3098                              clipthresh=None, clipniter=None,
     3099                              edge=None, threshold=None, chan_avg_limit=None,
     3100                              getresidual=None, plot=None,
     3101                              showprogress=None, minnrow=None, outlog=None,
     3102                              blfile=None, csvformat=None):
     3103        """\
     3104        Return a scan which has been baselined (all rows) by Chebyshev polynomials.
     3105        Spectral lines are detected first using linefinder and masked out
     3106        to avoid them affecting the baseline solution.
     3107
     3108        Parameters:
     3109            insitu:         if False a new scantable is returned.
     3110                            Otherwise, the scaling is done in-situ
     3111                            The default is taken from .asaprc (False)
     3112            mask:           an optional mask retreived from scantable
     3113            order:          the maximum order of Chebyshev polynomial (default is 5)
     3114            clipthresh:     Clipping threshold. (default is 3.0, unit: sigma)
     3115            clipniter:      maximum number of iteration of 'clipthresh'-sigma
     3116                            clipping (default is 0)
     3117            edge:           an optional number of channel to drop at
     3118                            the edge of spectrum. If only one value is
     3119                            specified, the same number will be dropped
     3120                            from both sides of the spectrum. Default
     3121                            is to keep all channels. Nested tuples
     3122                            represent individual edge selection for
     3123                            different IFs (a number of spectral channels
     3124                            can be different)
     3125            threshold:      the threshold used by line finder. It is
     3126                            better to keep it large as only strong lines
     3127                            affect the baseline solution.
     3128            chan_avg_limit: a maximum number of consequtive spectral
     3129                            channels to average during the search of
     3130                            weak and broad lines. The default is no
     3131                            averaging (and no search for weak lines).
     3132                            If such lines can affect the fitted baseline
     3133                            (e.g. a high order polynomial is fitted),
     3134                            increase this parameter (usually values up
     3135                            to 8 are reasonable). Most users of this
     3136                            method should find the default value sufficient.
     3137            plot:       *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
     3138                            plot the fit and the residual. In this each
     3139                            indivual fit has to be approved, by typing 'y'
     3140                            or 'n'
     3141            getresidual:    if False, returns best-fit values instead of
     3142                            residual. (default is True)
     3143            showprogress:   show progress status for large data.
     3144                            default is True.
     3145            minnrow:        minimum number of input spectra to show.
     3146                            default is 1000.
     3147            outlog:         Output the coefficients of the best-fit
     3148                            function to logger (default is False)
     3149            blfile:         Name of a text file in which the best-fit
     3150                            parameter values to be written
     3151                            (default is "": no file/logger output)
     3152            csvformat:      if True blfile is csv-formatted, default is False.
     3153
     3154        Example:
     3155            bscan = scan.auto_cspline_baseline(npiece=3, insitu=False)
     3156       
     3157        Note:
     3158            The best-fit parameter values output in logger and/or blfile are now
     3159            based on specunit of 'channel'.
     3160        """
     3161
     3162        try:
     3163            varlist = vars()
     3164
     3165            if insitu is None: insitu = rcParams['insitu']
     3166            if insitu:
     3167                workscan = self
     3168            else:
     3169                workscan = self.copy()
     3170           
     3171            #if mask           is None: mask           = [True for i in xrange(workscan.nchan())]
     3172            if mask           is None: mask           = []
     3173            if order          is None: order          = 5
     3174            if clipthresh     is None: clipthresh     = 3.0
     3175            if clipniter      is None: clipniter      = 0
     3176            if edge           is None: edge           = (0, 0)
     3177            if threshold      is None: threshold      = 3
     3178            if chan_avg_limit is None: chan_avg_limit = 1
     3179            if plot           is None: plot           = False
     3180            if getresidual    is None: getresidual    = True
     3181            if showprogress   is None: showprogress   = True
     3182            if minnrow        is None: minnrow        = 1000
     3183            if outlog         is None: outlog         = False
     3184            if blfile         is None: blfile         = ''
     3185            if csvformat      is None: csvformat      = False
     3186
     3187            if csvformat:
     3188                scsvformat = "T"
     3189            else:
     3190                scsvformat = "F"
     3191
     3192            #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method.
     3193            workscan._auto_chebyshev_baseline(mask, order, clipthresh,
     3194                                              clipniter,
     3195                                              normalise_edge_param(edge),
     3196                                              threshold,
     3197                                              chan_avg_limit, getresidual,
     3198                                              pack_progress_params(showprogress,
     3199                                                                   minnrow),
     3200                                              outlog, scsvformat+blfile)
     3201            workscan._add_history("auto_chebyshev_baseline", varlist)
     3202           
     3203            if insitu:
     3204                self._assign(workscan)
     3205            else:
     3206                return workscan
     3207           
     3208        except RuntimeError, e:
     3209            raise_fitting_failure_exception(e)
     3210
     3211    @asaplog_post_dec
    30093212    def poly_baseline(self, mask=None, order=None, insitu=None, plot=None,
    30103213                      getresidual=None, showprogress=None, minnrow=None,

  • 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{

  • trunk/src/Scantable.h

    r2641 r2645  
    522522                        const bool outLogger=false,
    523523                        const std::string& blfile="");
     524  void chebyshevBaseline(const std::vector<bool>& mask,
     525                         int order,
     526                         float thresClip,
     527                         int nIterClip,
     528                         bool getResidual=true,
     529                         const std::string& progressInfo="true,1000",
     530                         const bool outLogger=false,
     531                         const std::string& blfile="");
     532  void autoChebyshevBaseline(const std::vector<bool>& mask,
     533                             int order,
     534                             float thresClip,
     535                             int nIterClip,
     536                             const std::vector<int>& edge,
     537                             float threshold=3.0,
     538                             int chanAvgLimit=1,
     539                             bool getResidual=true,
     540                             const std::string& progressInfo="true,1000",
     541                             const bool outLogger=false,
     542                             const std::string& blfile="");
    524543  void cubicSplineBaseline(const std::vector<bool>& mask,
    525544                           int nPiece,
    … …  
    710729
    711730  void fitBaseline(const std::vector<bool>& mask, int whichrow, Fitter& fitter);
     731  double getChebyshevPolynomial(int n, double x);
     732  std::vector<float> doChebyshevFitting(const std::vector<float>& data,
     733                                          const std::vector<bool>& mask,
     734                                          int order,
     735                                          std::vector<float>& params,
     736                                          int& nClipped,
     737                                          float thresClip=3.0,
     738                                          int nIterClip=1,
     739                                          bool getResidual=true);
    712740  std::vector<float> doCubicSplineFitting(const std::vector<float>& data,
    713741                                          const std::vector<bool>& mask,

  • trunk/src/ScantableWrapper.h

    r2641 r2645  
    283283  { table_->autoPolyBaseline(mask, order, edge, threshold, chan_avg_limit, getresidual, showprogress, outlog, blfile); }
    284284
     285  void chebyshevBaseline(const std::vector<bool>& mask, int order, float clipthresh, int clipniter, bool getresidual=true, const std::string& showprogress="true,1000", const bool outlog=false, const std::string& blfile="")
     286  { table_->chebyshevBaseline(mask, order, clipthresh, clipniter, getresidual, showprogress, outlog, blfile); }
     287
     288  void autoChebyshevBaseline(const std::vector<bool>& mask, int order, float clipthresh, int clipniter, const std::vector<int>& edge, float threshold=5.0, int chan_avg_limit=1, bool getresidual=true, const std::string& showprogress="true,1000", const bool outlog=false, const std::string& blfile="")
     289  { table_->autoChebyshevBaseline(mask, order, clipthresh, clipniter, edge, threshold, chan_avg_limit, getresidual, showprogress, outlog, blfile); }
     290
    285291  void cubicSplineBaseline(const std::vector<bool>& mask, int npiece, float clipthresh, int clipniter, bool getresidual=true, const std::string& showprogress="true,1000", const bool outlog=false, const std::string& blfile="")
    286292  { table_->cubicSplineBaseline(mask, npiece, clipthresh, clipniter, getresidual, showprogress, outlog, blfile); }

  • trunk/src/python_Scantable.cpp

    r2591 r2645  
    149149    .def("_poly_baseline", &ScantableWrapper::polyBaseline)
    150150    .def("_auto_poly_baseline", &ScantableWrapper::autoPolyBaseline)
     151    .def("_chebyshev_baseline", &ScantableWrapper::chebyshevBaseline)
     152    .def("_auto_chebyshev_baseline", &ScantableWrapper::autoChebyshevBaseline)
    151153    .def("_cspline_baseline", &ScantableWrapper::cubicSplineBaseline)
    152154    .def("_auto_cspline_baseline", &ScantableWrapper::autoCubicSplineBaseline)
Note: See TracChangeset for help on using the changeset viewer.