    // @(#)root/graf:$Name:  $:$Id: TGraphSmooth.cxx,v 1.7 2002/12/02 18:50:02 rdm Exp $
    // Author: Christian Stratowa 30/09/2001

    /******************************************************************************
    * Copyright(c) 2001-    , Dr. Christian Stratowa, Vienna, Austria.            *
    * Author: Christian Stratowa with help from Rene Brun.                                                 *
    *                                                                             *
    * Algorithms for smooth regression adapted from:                              *
    * R: A Computer Language for Statistical Data Analysis                        *
    * Copyright (C) 1996 Robert Gentleman and Ross Ihaka                          *
    * Copyright (C) 1999-2001 Robert Gentleman, Ross Ihaka and the                *
    * R Development Core Team                                                     *
    * R is free software, for licensing see the GNU General Public License        *
    * http://www.ci.tuwien.ac.at/R-project/welcome.html                           *
    *                                                                             *
    * Based on: "The ROOT System", All rights reserved.                           *
    * Authors: Rene Brun and Fons Rademakers.                                     *
    * For the licensing terms of "The ROOT System" see $ROOTSYS/AA_LICENSE.       *
    * For the list of contributors to "The ROOT System" see $ROOTSYS/AA_CREDITS.  *
    ******************************************************************************/

    #include "Riostream.h"
    #include "TMath.h"
    #include "TGraphSmooth.h"
    #include "TGraphErrors.h"

    ClassImp(TGraphSmooth);

    //______________________________________________________________________
    // TGraphSmooth
    //
    // A helper class to smooth TGraph
    // see examples in $ROOTSYS/tutorials/motorcycle.C and approx.C
    //
    //______________________________________________________________________
     TGraphSmooth::TGraphSmooth(): TNamed()
    {
    //*-*-*-*-*-*-*-*-*Default GraphSmooth constructor *-*-*-*-*-*-*-*-*-*-*
    //                 ===============================

       fGin  = 0;
       fGout = 0;
    }

    //______________________________________________________________________
     TGraphSmooth::TGraphSmooth(const char *name): TNamed(name,"")
    {
    //*-*-*-*-*-*-*-*-*GraphSmooth constructor *-*-*-*-*-*-*-*-*-*-*-*-*-*-*
    //                 =======================

       fGin  = 0;
       fGout = 0;
    }

    //______________________________________________________________________
     TGraphSmooth::~TGraphSmooth()
    {
    //*-*-*-*-*-*-*-*-*GraphSmooth destructor*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
    //                 ======================

       if (fGout) delete fGout;
       fGin  = 0;
       fGout = 0;
    }

    //______________________________________________________________________
     void TGraphSmooth::Smoothin(TGraph *grin)
    {
    //*-*-*-*-*-*-*-*-*Sort input data points*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
    //                 ======================

       if (fGout) {delete fGout; fGout = 0;}
       fGin = grin;

       fNin = fGin->GetN();
       Double_t *xin = new Double_t[fNin];
       Double_t *yin = new Double_t[fNin];
       Int_t i;
       for (i=0;i<fNin;i++) {
          xin[i] = fGin->GetX()[i];
          yin[i] = fGin->GetY()[i];
       }

    // sort input x, y
       Int_t *index = new Int_t[fNin];
       TMath::Sort(fNin, xin, index, kFALSE);
       for (i=0;i<fNin;i++) {
          fGin->SetPoint(i, xin[index[i]], yin[index[i]]);
       }

       fMinX = fGin->GetX()[0];  //already sorted!
       fMaxX = fGin->GetX()[fNin-1];

       delete [] index;
       delete [] xin;
       delete [] yin;
    }

    //______________________________________________________________________
     TGraph *TGraphSmooth::SmoothKern(TGraph *grin, Option_t *option,
                          Double_t bandwidth, Int_t nout, Double_t *xout)
    {
    //*-*-*-*-*-*-*-*-*Smooth data with Kernel smoother*-*-*-*-*-*-*-*-*-*-*
    //                 ================================
    //
    // Smooth grin with the Nadaraya-Watson kernel regression estimate.
    //
    // Arguments:
    // grin:      input graph
    //
    // option:    the kernel to be used: "box", "normal"
    // bandwidth: the bandwidth. The kernels are scaled so that their quartiles
    //            (viewed as probability densities) are at +/- 0.25*bandwidth.
    // nout:      If xout is not specified, interpolation takes place at equally
    //            spaced points spanning the interval [min(x), max(x)], where
    //            nout = max(nout, number of input data).
    // xout:      an optional set of values at which to evaluate the fit
    //

       TString opt = option;
       opt.ToLower();
       Int_t kernel = 1;
       if (opt.Contains("normal")) kernel = 2;

       Smoothin(grin);

       Double_t delta = 0;
       Int_t *index = 0;
       if (xout == 0) {
          fNout = TMath::Max(nout, fNin);
          delta = (fMaxX - fMinX)/(fNout - 1);
       } else {
          fNout = nout;
          index = new Int_t[nout];
          TMath::Sort(nout, xout, index, kFALSE);
       }

       fGout = new TGraph(fNout);
       for (Int_t i=0;i<fNout;i++) {
          if (xout == 0) fGout->SetPoint(i,fMinX + i*delta, 0);
          else           fGout->SetPoint(i,xout[index[i]], 0);
       }

       BDRksmooth(fGin->GetX(), fGin->GetY(), fNin, fGout->GetX(),
                     fGout->GetY(), fNout, kernel, bandwidth);

       if (index) {delete [] index; index = 0;}

       return fGout;
    }

    //______________________________________________________________________
     void TGraphSmooth::BDRksmooth(Double_t *x, Double_t *y, Int_t n, Double_t *xp,
                       Double_t *yp, Int_t np, Int_t kernel, Double_t bw)
    {
    //*-*-*-*-*-*-*-*-*Smooth data with specified kernel*-*-*-*-*-*-*-*-*-*-*
    //*-*              =================================
    //
    //   Based on R function ksmooth: Translated to C++ by C. Stratowa
    //   (R source file: ksmooth.c by B.D.Ripley Copyright (C) 1998)
    //
    //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

       Int_t    imin = 0;
       Double_t cutoff = 0.0;

    // bandwidth is in units of half inter-quartile range
       if (kernel == 1) {
          bw *= 0.5;
          cutoff = bw;
       }
       if (kernel == 2) {
          bw *= 0.3706506;
          cutoff = 4*bw;
       }

       while ((x[imin] < xp[0] - cutoff) && (imin < n)) imin++;

       for (Int_t j=0;j<np;j++) {
          Double_t xx, w;
          Double_t num = 0.0;
          Double_t den = 0.0;
          Double_t x0 = xp[j];
          for (Int_t i=imin;i<n;i++) {
             if (x[i] < x0 - cutoff) imin = i;
             if (x[i] > x0 + cutoff) break;
             xx = TMath::Abs(x[i] - x0)/bw;
             if (kernel == 1) w = 1;
             else             w = TMath::Exp(-0.5*xx*xx);
             num += w*y[i];
             den += w;
          }
          if (den > 0) {
             yp[j] = num/den;
          } else {
             yp[j] = 0; //should be NA_REAL (see R.h) or nan("NAN")
          }
       }
    }

    //______________________________________________________________________
     TGraph *TGraphSmooth::SmoothLowess(TGraph *grin, Option_t *option ,
                          Double_t span, Int_t iter, Double_t delta)
    {
    //*-*-*-*-*-*-*-*-*Smooth data with Lowess smoother*-*-*-*-*-*-*-*-*-*-*
    //                 ================================
    //
    // This function performs the computations for the LOWESS smoother
    // (see the reference below). Lowess returns the output points
    // x and y which give the coordinates of the smooth.
    //
    // Arguments:
    // grin:  Input graph
    //
    // span:  the smoother span. This gives the proportion of points in the plot
    //        which influence the smooth at each value.
    //        Larger values give more smoothness.
    // iter:  the number of robustifying iterations which should be performed.
    //        Using smaller values of iter will make lowess run faster.
    // delta: values of x which lie within delta of each other replaced by a
    //        single value in the output from lowess.
    //        For delta = 0, delta will be calculated.
    //
    // References:
    // Cleveland, W. S. (1979) Robust locally weighted regression and smoothing
    //        scatterplots. J. Amer. Statist. Assoc. 74, 829-836.
    // Cleveland, W. S. (1981) LOWESS: A program for smoothing scatterplots
    //        by robust locally weighted regression.
    //        The American Statistician, 35, 54.
    //                 ==================

       TString opt = option;
       opt.ToLower();

       Smoothin(grin);

       if (delta == 0) {delta = 0.01*(TMath::Abs(fMaxX - fMinX));}

    // output X, Y
       fNout = fNin;
       fGout = new TGraphErrors(fNout);

       for (Int_t i=0;i<fNout;i++) {
          fGout->SetPoint(i,fGin->GetX()[i], 0);
       }

       Lowess(fGin->GetX(), fGin->GetY(), fNin, fGout->GetY(), span, iter, delta);

       return fGout;
    }

    //______________________________________________________________________
     void TGraphSmooth::Lowess(Double_t *x, Double_t *y, Int_t n, Double_t *ys,
                       Double_t span, Int_t iter, Double_t delta)
    {
    //*-*-*-*-*-*-*-*-*Lowess regression smoother*-*-*-*-*-*-*-*-*-*-*-*-*-*
    //                 ==========================
    //   Based on R function clowess: Translated to C++ by C. Stratowa
    //   (R source file: lowess.c by R Development Core Team (C) 1999-2001)
    //
    //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

       Int_t    i, iiter, j, last, m1, m2, nleft, nright, ns;
       Double_t alpha, c1, c9, cmad, cut, d1, d2, denom, r;
       Bool_t   ok;

       if (n < 2) {
    	   ys[0] = y[0];
          return;
       }

    // nleft, nright, last, etc. must all be shifted to get rid of these:
       x--;
       y--;
       ys--;

       Double_t *RW  = ((TGraphErrors*)fGout)->GetEX();
       Double_t *res = ((TGraphErrors*)fGout)->GetEY();

    // at least two, at most n poInt_ts
       ns = TMath::Max(2, TMath::Min(n, (Int_t)(span*n + 1e-7)));

    // robustness iterations
       iiter = 1;
       while (iiter <= iter+1) {
    	   nleft = 1;
    	   nright = ns;
    	   last = 0;	// index of prev estimated poInt_t
    	   i = 1;		// index of current poInt_t

          for(;;) {
    	      if (nright < n) {
             // move nleft,  nright to right if radius decreases
                d1 = x[i] - x[nleft];
    		      d2 = x[nright+1] - x[i];

    		   // if d1 <= d2 with x[nright+1] == x[nright], lowest fixes
                if (d1 > d2) {
                // radius will not decrease by move right
                   nleft++;
    		         nright++;
    		         continue;
    		      }
    	      }

    	   // fitted value at x[i]
             Bool_t iterg1 = iiter>1;
             Lowest(&x[1], &y[1], n, x[i], ys[i], nleft, nright,
                          res, iterg1, RW, ok);
    	      if (!ok) ys[i] = y[i];

    	   // all weights zero copy over value (all rw==0)
             if (last < i-1) {
    		      denom = x[i]-x[last];

    		   // skipped poInt_ts -- Int_terpolate non-zero - proof?
                for(j = last+1; j < i; j++) {
    		         alpha = (x[j]-x[last])/denom;
    		         ys[j] = alpha*ys[i] + (1.-alpha)*ys[last];
                }
    	 }

    	   // last poInt_t actually estimated
    	      last = i;

    	   // x coord of close poInt_ts
    	      cut = x[last] + delta;
    	      for (i = last+1; i <= n; i++) {
    		      if (x[i] > cut)
    		         break;
    		      if (x[i] == x[last]) {
    		         ys[i] = ys[last];
    		         last = i;
    		      }
    	      }
    	      i = TMath::Max(last+1, i-1);
    	      if (last >= n)
    		      break;
    	   }

    	// residuals
    	   for(i=0; i < n; i++)
    	      res[i] = y[i+1] - ys[i+1];

    	// compute robustness weights except last time
          if (iiter > iter)
    	      break;
    	   for(i=0 ; i<n ; i++)
    	      RW[i] = TMath::Abs(res[i]);

    	// compute cmad := 6 * median(rw[], n)
       	m1 = n/2;
    	// partial sort, for m1 & m2
    	   Psort(RW, n, m1);
    	   if(n % 2 == 0) {
    	      m2 = n-m1-1;
    	      Psort(RW, n, m2);
    	      cmad = 3.*(RW[m1]+RW[m2]);
    	   } else { /* n odd */
    	      cmad = 6.*RW[m1];
    	   }

    	   c9 = 0.999*cmad;
    	   c1 = 0.001*cmad;
    	   for(i=0 ; i<n ; i++) {
    	      r = TMath::Abs(res[i]);
    	      if (r <= c1)
    		      RW[i] = 1.;
    	      else if (r <= c9)
    		      RW[i] = (1.-(r/cmad)*(r/cmad))*(1.-(r/cmad)*(r/cmad));
    	      else
    		      RW[i] = 0.;
    	   }
    	   iiter++;
       }
    }

    //______________________________________________________________________
     void TGraphSmooth::Lowest(Double_t *x, Double_t *y, Int_t n, Double_t &xs,
                       Double_t &ys, Int_t nleft, Int_t nright, Double_t *w,
                       Bool_t userw, Double_t *rw, Bool_t &ok)
    {
    //*-*-*-*-*-*-*-*-*Fit value at x[i] *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
    //                 =================
    //  Based on R function lowest: Translated to C++ by C. Stratowa
    //  (R source file: lowess.c by R Development Core Team (C) 1999-2001)
    //
    //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

       Int_t    nrt, j;
       Double_t a, b, c, d, h, h1, h9, r, range;

       x--;
       y--;
       w--;
       rw--;

       range = x[n]-x[1];
       h = TMath::Max(xs-x[nleft], x[nright]-xs);
       h9 = 0.999*h;
       h1 = 0.001*h;

    // sum of weights
       a = 0.;
       j = nleft;
       while (j <= n) {
    	// compute weights (pick up all ties on right)
          w[j] = 0.;
    	   r = TMath::Abs(x[j] - xs);
    	   if (r <= h9) {
    	      if (r <= h1) {
    		      w[j] = 1.;
    	      } else {
                d = (r/h)*(r/h)*(r/h);
    		      w[j] = (1.- d)*(1.- d)*(1.- d);
             }
    	      if (userw)
    		      w[j] *= rw[j];
    	      a += w[j];
    	   } else if (x[j] > xs)
    	      break;
    	   j = j+1;
       }

    // rightmost pt (may be greater than nright because of ties)
       nrt = j-1;
       if (a <= 0.)
    	   ok = kFALSE;
       else {
    	   ok = kTRUE;
       // weighted least squares: make sum of w[j] == 1
    	   for(j=nleft ; j<=nrt ; j++)
    	      w[j] /= a;
    	   if (h > 0.) {
    	      a = 0.;
          // use linear fit weighted center of x values
             for(j=nleft ; j<=nrt ; j++)
    		      a += w[j] * x[j];
    	      b = xs - a;
    	      c = 0.;
    	      for(j=nleft ; j<=nrt ; j++)
    		      c += w[j]*(x[j]-a)*(x[j]-a);
    	      if (TMath::Sqrt(c) > 0.001*range) {
    		      b /= c;
             // poInt_ts are spread out enough to compute slope
                for(j=nleft; j <= nrt; j++)
    		         w[j] *= (b*(x[j]-a) + 1.);
    	      }
    	   }
    	   ys = 0.;
    	   for(j=nleft; j <= nrt; j++)
    	      ys += w[j] * y[j];
       }
    }

    //______________________________________________________________________
     TGraph *TGraphSmooth::SmoothSuper(TGraph *grin, Option_t *,
            Double_t bass, Double_t span, Bool_t isPeriodic, Double_t *w)
    {
    //*-*-*-*-*-*-*-*-*Smooth data with Super smoother*-*-*-*-*-*-*-*-*-*-*-*
    //                 ===============================
    //
    // Smooth the (x, y) values by Friedman's ``super smoother''.
    //
    // Arguments:
    // grin: graph for smoothing
    //
    // span: the fraction of the observations in the span of the running lines
    //        smoother, or 0 to choose this by leave-one-out cross-validation.
    // bass: controls the smoothness of the fitted curve.
    //        Values of up to 10 indicate increasing smoothness.
    // isPeriodic: if TRUE, the x values are assumed to be in [0, 1]
    //        and of period 1.
    // w:     case weights
    //
    // Details:
    // supsmu is a running lines smoother which chooses between three spans for
    // the lines. The running lines smoothers are symmetric, with k/2 data points
    // each side of the predicted point, and values of k as 0.5 * n, 0.2 * n and
    // 0.05 * n, where n is the number of data points. If span is specified,
    // a single smoother with span span * n is used.
    //
    // The best of the three smoothers is chosen by cross-validation for each
    // prediction. The best spans are then smoothed by a running lines smoother
    // and the final prediction chosen by linear interpolation.
    //
    // The FORTRAN code says: ``For small samples (n < 40) or if there are
    // substantial serial correlations between observations close in x - value,
    // then a prespecified fixed span smoother (span > 0) should be used.
    // Reasonable span values are 0.2 to 0.4.''
    //
    // References:
    // Friedman, J. H. (1984) SMART User's Guide.
    //           Laboratory for Computational Statistics,
    //           Stanford University Technical Report No. 1.
    //
    // Friedman, J. H. (1984) A variable span scatterplot smoother.
    //           Laboratory for Computational Statistics,
    //           Stanford University Technical Report No. 5.
    //                 ==================

       if (span < 0 || span > 1) {
          cout << "Error: Span must be between 0 and 1" << endl;
          return 0;
       }

       Smoothin(grin);

       Int_t iper = 1;
       if (isPeriodic) {
          iper = 2;
          if (fMinX < 0 || fMaxX > 1) {
             cout << "Error: x must be between 0 and 1 for periodic smooth" << endl;
             return 0;
          }
       }

    // output X, Y
       fNout = fNin;
       fGout = new TGraph(fNout);
       Int_t i;
       for (i=0; i<fNout; i++) {
          fGout->SetPoint(i,fGin->GetX()[i], 0);
       }

    // weights
       Double_t *weight = new Double_t[fNin];
       for (i=0; i<fNin; i++) {
          if (w == 0) weight[i] = 1;
          else        weight[i] = w[i];
       }

    // temporary storage array
       Int_t NTmp = (fNin+1)*8;
       Double_t *Tmp = new Double_t[NTmp];
       for (i=0; i<NTmp; i++) {
          Tmp[i] = 0;
       }

       BDRsupsmu(fNin, fGin->GetX(), fGin->GetY(), weight, iper, span, bass, fGout->GetY(), Tmp);

       delete [] Tmp;
       delete [] weight;

       return fGout;
    }

    //______________________________________________________________________
     void TGraphSmooth::BDRsupsmu(Int_t n, Double_t *x, Double_t *y, Double_t *w,
         Int_t iper, Double_t span, Double_t alpha, Double_t *smo, Double_t *sc)
    {
    //*-*-*-*-*-*-*-*-*Friedmannīs super smoother *-*-*-*-*-*-*-*-*-*-*-*-*-*
    //                 ==========================
    //
    //  super smoother (Friedman, 1984).
    //
    //  version 10/10/84
    //
    //  coded  and copywrite (c) 1984 by:
    //
    //                         Jerome H. Friedman
    //                      department of statistics
    //                                and
    //                 stanford linear accelerator center
    //                         stanford university
    //
    //  all rights reserved.
    //
    //
    //  input:
    //     n : number of observations (x,y - pairs).
    //     x(n) : ordered abscissa values.
    //     y(n) : corresponding ordinate (response) values.
    //     w(n) : weight for each (x,y) observation.
    //     iper : periodic variable flag.
    //        iper=1 => x is ordered interval variable.
    //        iper=2 => x is a periodic variable with values
    //                  in the range (0.0,1.0) and period 1.0.
    //     span : smoother span (fraction of observations in window).
    //            span=0.0 => automatic (variable) span selection.
    //     alpha : controls high frequency (small span) penality
    //             used with automatic span selection (bass tone control).
    //             (alpha.le.0.0 or alpha.gt.10.0 => no effect.)
    //  output:
    //    smo(n) : smoothed ordinate (response) values.
    //  scratch:
    //    sc(n,7) : internal working storage.
    //
    //  note:
    //     for small samples (n < 40) or if there are substantial serial
    //     correlations between observations close in x - value, then
    //     a prespecified fixed span smoother (span > 0) should be
    //     used. reasonable span values are 0.2 to 0.4.
    //
    // current implementation:
    //   Based on R function supsmu: Translated to C++ by C. Stratowa
    //   (R source file: ppr.f by B.D.Ripley Copyright (C) 1994-97)
    //
    //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

    // Local variables
       Int_t sc_offset;
       Int_t i, j, jper;
       Double_t a, f;
       Double_t sw, sy, resmin, vsmlsq;
       Double_t scale;
       Double_t d1, d2;

       Double_t spans[3] = { 0.05, 0.2, 0.5 };
       Double_t big = 1e20;
       Double_t sml = 1e-7;
       Double_t eps = 0.001;

    // Parameter adjustments
       sc_offset = n + 1;
       sc -= sc_offset;
       --smo;
       --w;
       --y;
       --x;

    // Function Body
       if (x[n] <= x[1]) {
          sy = 0.0;
          sw = sy;
          for (j=1;j<=n;++j) {
             sy += w[j] * y[j];
             sw += w[j];
          }

          a = 0.0;
          if (sw > 0.0) a = sy / sw;
          for (j=1;j<=n;++j) smo[j] = a;
          return;
       }

       i = (Int_t)(n / 4);
       j = i * 3;
       scale = x[j] - x[i];
       while (scale <= 0.0) {
          if (j < n) ++j;
          if (i > 1) --i;
          scale = x[j] - x[i];
       }

    // Computing 2nd power
       d1 = eps * scale;
       vsmlsq = d1 * d1;
       jper = iper;
       if (iper == 2 && (x[1] < 0.0 || x[n] > 1.0)) {
          jper = 1;
       }
       if (jper < 1 || jper > 2) {
          jper = 1;
       }
       if (span > 0.0) {
          BDRsmooth(n, &x[1], &y[1], &w[1], span, jper, vsmlsq,
                          &smo[1], &sc[sc_offset]);
          return;
       }

       Double_t *h = new Double_t[n+1];
       for (i = 1; i <= 3; ++i) {
          BDRsmooth(n, &x[1], &y[1], &w[1], spans[i - 1], jper, vsmlsq,
                          &sc[((i<<1)-1)*n + 1], &sc[n*7 + 1]);
          BDRsmooth(n, &x[1], &sc[n*7 + 1], &w[1], spans[1], -jper, vsmlsq,
                          &sc[(i<<1)*n + 1], &h[1]);
       }

       for (j=1; j<=n; ++j) {
          resmin = big;
          for (i=1; i<=3; ++i) {
             if (sc[j + (i<<1)*n] < resmin) {
                resmin = sc[j + (i<<1)*n];
                sc[j + n*7] = spans[i-1];
             }
          }

          if (alpha>0.0 && alpha<=10.0 && resmin<sc[j + n*6]  && resmin>0.0) {
          // Computing MAX
             d1 = TMath::Max(sml,(resmin/sc[j + n*6]));
             d2 = 10. - alpha;
             sc[j + n*7] += (spans[2] - sc[j + n*7]) * TMath::Power(d1, d2);
           }
       }

       BDRsmooth(n, &x[1], &sc[n*7 + 1], &w[1], spans[1], -jper, vsmlsq,
                       &sc[(n<<1) + 1], &h[1]);

       for (j=1; j<=n; ++j) {
          if (sc[j + (n<<1)] <= spans[0]) {
             sc[j + (n<<1)] = spans[0];
          }
          if (sc[j + (n<<1)] >= spans[2]) {
             sc[j + (n<<1)] = spans[2];
          }
          f = sc[j + (n<<1)] - spans[1];
          if (f < 0.0) {
             f = -f / (spans[1] - spans[0]);
             sc[j + (n<<2)] = (1.0 - f) * sc[j + n*3] + f * sc[j + n];
          } else {
             f /= spans[2] - spans[1];
             sc[j + (n<<2)] = (1.0 - f) * sc[j + n*3] + f * sc[j + n*5];
          }
       }

       BDRsmooth(n, &x[1], &sc[(n<<2) + 1], &w[1], spans[0], -jper, vsmlsq,
                       &smo[1], &h[1]);

       delete [] h;
       return;
    }

    //______________________________________________________________________
     void TGraphSmooth::BDRsmooth(Int_t n, Double_t *x, Double_t *y, Double_t *w,
         Double_t span, Int_t iper, Double_t vsmlsq, Double_t *smo, Double_t *acvr)
    {
    //*-*-*-*-*-*-*-*-* Function for super smoother *-*-*-*-*-*-*-*-*-*-*-*
    //                 ============================
    //
    //   Based on R function supsmu: Translated to C++ by C. Stratowa
    //   (R source file: ppr.f by B.D.Ripley Copyright (C) 1994-97)
    //
    //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

    // Local variables
       Int_t i, j, j0, in, out, it, jper, ibw;
       Double_t a, h1, d1;
       Double_t xm, ym, wt, sy, fbo, fbw;
       Double_t cvar, var, tmp, xti, xto;

    // Parameter adjustments
       --acvr;
       --smo;
       --w;
       --y;
       --x;

    // Function Body
       xm = 0.;
       ym = xm;
       var = ym;
       cvar = var;
       fbw = cvar;
       jper = TMath::Abs(iper);

       ibw = (Int_t)(span * 0.5 * n + 0.5);
       if (ibw < 2) {
          ibw = 2;
       }

       it = 2*ibw + 1;
       for (i=1; i<=it; ++i) {
          j = i;
          if (jper == 2) {
             j = i - ibw - 1;
          }
          xti = x[j];
          if (j < 1) {
             j = n + j;
             xti = x[j] - 1.0;
          }
          wt = w[j];
          fbo = fbw;
          fbw += wt;
          if (fbw > 0.0) {
             xm = (fbo * xm + wt * xti) / fbw;
             ym = (fbo * ym + wt * y[j]) / fbw;
          }
          tmp = 0.0;
          if (fbo > 0.0) {
             tmp = fbw * wt * (xti - xm) / fbo;
          }
          var += tmp * (xti - xm);
          cvar += tmp * (y[j] - ym);
       }

       for (j=1; j<=n; ++j) {
          out = j - ibw - 1;
          in = j + ibw;
          if (!(jper != 2 && (out < 1 || in > n))) {
             if (out < 1) {
                out = n + out;
                xto = x[out] - 1.0;
                xti = x[in];
             } else if (in > n) {
                in -= n;
                xti = x[in] + 1.0;
                xto = x[out];
             } else {
                xto = x[out];
                xti = x[in];
             }

             wt = w[out];
             fbo = fbw;
             fbw -= wt;
             tmp = 0.0;
             if (fbw > 0.0) {
                tmp = fbo * wt * (xto - xm) / fbw;
             }
             var -= tmp * (xto - xm);
             cvar -= tmp * (y[out] - ym);
             if (fbw > 0.0) {
                xm = (fbo * xm - wt * xto) / fbw;
                ym = (fbo * ym - wt * y[out]) / fbw;
             }
             wt = w[in];
             fbo = fbw;
             fbw += wt;
             if (fbw > 0.0) {
                xm = (fbo * xm + wt * xti) / fbw;
                ym = (fbo * ym + wt * y[in]) / fbw;
             }
             tmp = 0.0;
             if (fbo > 0.0) {
                tmp = fbw * wt * (xti - xm) / fbo;
             }
             var += tmp * (xti - xm);
             cvar += tmp * (y[in] - ym);
          }

          a = 0.0;
          if (var > vsmlsq) {
             a = cvar / var;
          }
          smo[j] = a * (x[j] - xm) + ym;

          if (iper <= 0) {
             continue;
          }

          h1 = 0.0;
          if (fbw > 0.0) {
             h1 = 1.0 / fbw;
          }
          if (var > vsmlsq) {
          // Computing 2nd power
             d1 = x[j] - xm;
             h1 += d1 * d1 / var;
          }

          acvr[j] = 0.0;
          a = 1.0 - w[j] * h1;
          if (a > 0.0) {
             acvr[j] = TMath::Abs(y[j] - smo[j]) / a;
             continue;
          }
          if (j > 1) {
             acvr[j] = acvr[j-1];
          }
       }

       j = 1;
       do {
          j0 = j;
          sy = smo[j] * w[j];
          fbw = w[j];
          if (j < n) {
             do {
                if (x[j + 1] > x[j]) {
                   break;
                }
                ++j;
                sy += w[j] * smo[j];
                fbw += w[j];
             } while (j < n);
          }

          if (j > j0) {
             a = 0.0;
             if (fbw > 0.0) {
                a = sy / fbw;
             }
             for (i=j0; i<=j; ++i) {
                smo[i] = a;
             }
          }
          ++j;
       } while (j <= n);

       return;
    }

    //______________________________________________________________________
     void TGraphSmooth::Approxin(TGraph *grin, Int_t /*iKind*/, Double_t &Ylow,
         Double_t &Yhigh, Int_t rule, Int_t iTies)
    {
    //*-*-*-*-*-*-*-*-*Sort data points and eliminate double x values*-*-*-*
    //                 ==============================================

       if (fGout) {delete fGout; fGout = 0;}
       fGin = grin;

       fNin = fGin->GetN();
       Double_t *xin = new Double_t[fNin];
       Double_t *yin = new Double_t[fNin];
       Int_t i;
       for (i=0;i<fNin;i++) {
          xin[i] = fGin->GetX()[i];
          yin[i] = fGin->GetY()[i];
       }

    // sort/rank input x, y
       Int_t *index = new Int_t[fNin];
       Int_t *rank  = new Int_t[fNin];
       Rank(fNin, xin, index, rank, kFALSE);

    // input X, Y
       Int_t vNDup = 0;
       Int_t k = 0;
       Int_t *dup = new Int_t[fNin];
       Double_t *x = new Double_t[fNin];
       Double_t *y = new Double_t[fNin];
       Double_t vMean, vMin, vMax;
       for (i=1;i<fNin+1;i++) {
          Int_t ndup = 1;
          vMin = vMean = vMax = yin[index[i-1]];
          while ((i < fNin) && (rank[index[i]] == rank[index[i-1]])) {
             vMean += yin[index[i]];
             vMax = (vMax < yin[index[i]]) ? yin[index[i]] : vMax;
             vMin = (vMin > yin[index[i]]) ? yin[index[i]] : vMin;
             dup[vNDup] = i;
             i++;
             ndup++;
             vNDup++;
          }
          x[k] = xin[index[i-1]];
          if (ndup == 1) {y[k++] = yin[index[i-1]];}
          else switch(iTies) {
             case 1:
                y[k++] = vMean/ndup;
       	      break;
             case 2:
                y[k++] = vMin;
    	         break;
             case 3:
                y[k++] = vMax;
        	      break;
             default:
                y[k++] = vMean/ndup;
       	      break;
          }
       }
       fNin = k;

    // set unique sorted input data x,y as final graph points
       fGin->Set(fNin);
       for (i=0;i<fNin;i++) {
          fGin->SetPoint(i, x[i], y[i]);
       }

       fMinX = fGin->GetX()[0];  //already sorted!
       fMaxX = fGin->GetX()[fNin-1];

    // interpolate outside interval [min(x),max(x)]
       switch(rule) {
          case 1:
             Ylow  = 0;   // = nan("NAN") ??
             Yhigh = 0;   // = nan("NAN") ??
    	 break;
          case 2:
             Ylow  = fGin->GetY()[0];
             Yhigh = fGin->GetY()[fNin-1];
    	 break;
          default:
    	 break;
       }

    // cleanup
       delete [] x;
       delete [] y;
       delete [] dup;
       delete [] rank;
       delete [] index;
       delete [] xin;
       delete [] yin;
    }

    //______________________________________________________________________
     TGraph *TGraphSmooth::Approx(TGraph *grin, Option_t *option, Int_t nout, Double_t *xout,
            Double_t yleft, Double_t yright, Int_t rule, Double_t f, Option_t *ties)
    {
    //*-*-*-*-*-*-*-*-*Approximate data points*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
    //                 =======================
    //
    // Arguments:
    // grin:  graph giving the coordinates of the points to be interpolated.
    //        Alternatively a single plotting structure can be specified:
    //
    // option: specifies the interpolation method to be used.
    //        Choices are "linear" (iKind = 1) or "constant" (iKind = 2).
    // nout:  If xout is not specified, interpolation takes place at n equally
    //        spaced points spanning the interval [min(x), max(x)], where
    //        nout = max(nout, number of input data).
    // xout:  an optional set of values specifying where interpolation is to
    //        take place.
    // yleft: the value to be returned when input x values less than min(x).
    //        The default is defined by the value of rule given below.
    // yright: the value to be returned when input x values greater than max(x).
    //        The default is defined by the value of rule given below.
    // rule:  an integer describing how interpolation is to take place outside
    //        the interval [min(x), max(x)]. If rule is 0 then the given yleft
    //        and yright values are returned, if it is 1 then 0 is returned
    //        for such points and if it is 2, the value at the closest data
    //        extreme is used.
    // f:     For method="constant" a number between 0 and 1 inclusive,
    //        indicating a compromise between left- and right-continuous step
    //        functions. If y0 and y1 are the values to the left and right of
    //        the point then the value is y0*f+y1*(1-f) so that f=0 is
    //        right-continuous and f=1 is left-continuous
    // ties:  Handling of tied x values. An integer describing a function with
    //        a single vector argument returning a single number result:
    //        ties = "ordered" (iTies = 0): input x are "ordered"
    //        ties = "mean"    (iTies = 1): function "mean"
    //        ties = "min"     (iTies = 2): function "min"
    //        ties = "max"     (iTies = 3): function "max"
    //
    // Details:
    // At least two complete (x, y) pairs are required.
    // If there are duplicated (tied) x values and ties is a function it is
    // applied to the y values for each distinct x value. Useful functions in
    // this context include mean, min, and max.
    // If ties="ordered" the x values are assumed to be already ordered. The
    // first y value will be used for interpolation to the left and the last
    // one for interpolation to the right.
    //
    // Value:
    // approx returns a graph with components x and y, containing n coordinates
    // which interpolate the given data points according to the method (and rule)
    // desired.

       TString opt = option;
       opt.ToLower();
       Int_t iKind = 0;
       if (opt.Contains("linear")) iKind = 1;
       else if (opt.Contains("constant")) iKind = 2;

       if (f < 0 || f > 1) {
          cout << "Error: Invalid f value" << endl;
          return 0;
       }

       opt = ties;
       opt.ToLower();
       Int_t iTies = 0;
       if (opt.Contains("ordered")) {
          iTies = 0;
       } else if (opt.Contains("mean")) {
          iTies = 1;
       } else if (opt.Contains("min")) {
          iTies = 2;
       } else if (opt.Contains("max")) {
          iTies = 3;
       } else {
          cout << "Error: Method not known: " << ties << endl;
          return 0;
       }

    // input X, Y
       Double_t Ylow  = yleft;
       Double_t Yhigh = yright;
       Approxin(grin, iKind, Ylow, Yhigh, rule, iTies);

    // output X, Y
       Double_t delta = 0;
       fNout = nout;
       if (xout == 0) {
          fNout = TMath::Max(nout, fNin);
          delta = (fMaxX - fMinX)/(fNout - 1);
       }

       fGout = new TGraph(fNout);

       Double_t x;
       for (Int_t i=0;i<fNout;i++) {
          if (xout == 0) x = fMinX + i*delta;
          else           x = xout[i];
          Double_t yout = Approx1(x, f, fGin->GetX(), fGin->GetY(), fNin, iKind, Ylow, Yhigh);
          fGout->SetPoint(i,x, yout);
       }

       return fGout;
    }

    //______________________________________________________________________
     Double_t TGraphSmooth::Approx1(Double_t v, Double_t f, Double_t *x, Double_t *y,
             Int_t n, Int_t iKind, Double_t Ylow, Double_t Yhigh)
    {
    //*-*-*-*-*-*-*-*-*Approximate one data point*-*-*-*-*-*-*-*-*-*-*-*-*-*
    //*-*              ==========================
    //
    //   Approximate  y(v),  given (x,y)[i], i = 0,..,n-1
    //   Based on R function approx1: Translated to C++ by Christian Stratowa
    //   (R source file: approx.c by R Development Core Team (C) 1999-2001)
    //

       Int_t i = 0;
       Int_t j = n - 1;

    // handle out-of-domain points
       if(v < x[i]) return Ylow;
       if(v > x[j]) return Yhigh;

    // find the correct interval by bisection
       while(i < j - 1) {
    	   Int_t ij = (i + j)/2;
    	   if(v < x[ij]) j = ij;
    	   else i = ij;
       }

    // interpolation
       if(v == x[j]) return y[j];
       if(v == x[i]) return y[i];

       if(iKind == 1) { // linear
    	   return y[i] + (y[j] - y[i]) * ((v - x[i])/(x[j] - x[i]));
       } else { // 2 : constant
    	   return y[i] * (1-f) + y[j] * f;
       }
    }

    // helper functions
    //______________________________________________________________________
     Int_t TGraphSmooth::Rcmp(Double_t x, Double_t y)
    {
    //   static function
    //   if (ISNAN(x))	return 1;
    //   if (ISNAN(y))	return -1;
       if (x < y)		return -1;
       if (x > y)		return 1;
       return 0;
    }

    //______________________________________________________________________
     void TGraphSmooth::Psort(Double_t *x, Int_t n, Int_t k)
    {
    //   static function
    //   based on R function rPsort: adapted to C++ by Christian Stratowa
    //   (R source file: R_sort.c by R Development Core Team (C) 1999-2001)
    //

       Double_t v, w;
       Int_t L, R, i, j;

       for (L = 0, R = n - 1; L < R; ) {
    	   v = x[k];
    	   for(i = L, j = R; i <= j;) {
    	      while (TGraphSmooth::Rcmp(x[i], v) < 0) i++;
    	      while (TGraphSmooth::Rcmp(v, x[j]) < 0) j--;
    	      if (i <= j) { w = x[i]; x[i++] = x[j]; x[j--] = w; }
    	   }
    	   if (j < k) L = i;
    	   if (k < i) R = j;
       }
    }

    //______________________________________________________________________
     void TGraphSmooth::Rank(Int_t n, Double_t *a, Int_t *index, Int_t *rank, Bool_t down)
    {
    //   static function
       if (n <= 0) return;
       if (n == 1) {
          index[0] = 0;
          rank[0] = 0;
          return;
       }

       TMath::Sort(n,a,index,down);

       Int_t k = 0;
       for (Int_t i=0;i<n;i++) {
          if ((i > 0) && (a[index[i]] == a[index[i-1]])) {
             rank[index[i]] = i-1;
             k++;
          }
          rank[index[i]] = i-k;
       }
    }