Logo ROOT   6.10/09
Reference Guide
TGraphSmooth.cxx
Go to the documentation of this file.
1 // @(#)root/hist:$Id$
2 // Author: Christian Stratowa 30/09/2001
3 
4 /*************************************************************************
5  * Copyright (C) 2006, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 /******************************************************************************
13 * Copyright(c) 2001- , Dr. Christian Stratowa, Vienna, Austria. *
14 * Author: Christian Stratowa with help from Rene Brun. *
15 * *
16 * Algorithms for smooth regression adapted from: *
17 * R: A Computer Language for Statistical Data Analysis *
18 * Copyright (C) 1996 Robert Gentleman and Ross Ihaka *
19 * Copyright (C) 1999-2001 Robert Gentleman, Ross Ihaka and the *
20 * R Development Core Team *
21 * R is free software, for licensing see the GNU General Public License *
22 * http://www.ci.tuwien.ac.at/R-project/welcome.html *
23 * *
24 ******************************************************************************/
25 
26 
27 #include "Riostream.h"
28 #include "TMath.h"
29 #include "TGraphSmooth.h"
30 #include "TGraphErrors.h"
31 
33 
34 //______________________________________________________________________
35 /** \class TGraphSmooth
36  \ingroup Hist
37 A helper class to smooth TGraph.
38 see examples in $ROOTSYS/tutorials/graphs/motorcycle.C and approx.C
39 */
40 
42 {
43  fNin = 0;
44  fNout = 0;
45  fGin = 0;
46  fGout = 0;
47  fMinX = 0;
48  fMaxX = 0;
49 }
50 
51 ////////////////////////////////////////////////////////////////////////////////
52 /// GraphSmooth constructor
53 
54 TGraphSmooth::TGraphSmooth(const char *name): TNamed(name,"")
55 {
56  fNin = 0;
57  fNout = 0;
58  fGin = 0;
59  fGout = 0;
60  fMinX = 0;
61  fMaxX = 0;
62 }
63 
64 ////////////////////////////////////////////////////////////////////////////////
65 /// GraphSmooth destructor
66 
68 {
69  if (fGout) delete fGout;
70  fGin = 0;
71  fGout = 0;
72 }
73 
74 ////////////////////////////////////////////////////////////////////////////////
75 /// Sort input data points
76 
78 {
79  if (fGout) {delete fGout; fGout = 0;}
80  fGin = grin;
81 
82  fNin = fGin->GetN();
83  Double_t *xin = new Double_t[fNin];
84  Double_t *yin = new Double_t[fNin];
85  Int_t i;
86  for (i=0;i<fNin;i++) {
87  xin[i] = fGin->GetX()[i];
88  yin[i] = fGin->GetY()[i];
89  }
90 
91 // sort input x, y
92  Int_t *index = new Int_t[fNin];
93  TMath::Sort(fNin, xin, index, kFALSE);
94  for (i=0;i<fNin;i++) {
95  fGin->SetPoint(i, xin[index[i]], yin[index[i]]);
96  }
97 
98  fMinX = fGin->GetX()[0]; //already sorted!
99  fMaxX = fGin->GetX()[fNin-1];
100 
101  delete [] index;
102  delete [] xin;
103  delete [] yin;
104 }
105 
106 ////////////////////////////////////////////////////////////////////////////////
107 /// Smooth data with Kernel smoother. Smooth grin with the Nadaraya-Watson kernel regression estimate.
108 ///
109 /// \param[in] grin input graph
110 /// \param[in] option the kernel to be used: "box", "normal"
111 /// \param[in] bandwidth the bandwidth. The kernels are scaled so that their quartiles
112 /// (viewed as probability densities) are at +/- 0.25*bandwidth.
113 /// \param[in] nout If xout is not specified, interpolation takes place at equally
114 /// spaced points spanning the interval [min(x), max(x)], where nout = max(nout, number of input data).
115 /// \param[in] xout an optional set of values at which to evaluate the fit
116 
118  Double_t bandwidth, Int_t nout, Double_t *xout)
119 {
120  TString opt = option;
121  opt.ToLower();
122  Int_t kernel = 1;
123  if (opt.Contains("normal")) kernel = 2;
124 
125  Smoothin(grin);
126 
127  Double_t delta = 0;
128  Int_t *index = 0;
129  if (xout == 0) {
130  fNout = TMath::Max(nout, fNin);
131  delta = (fMaxX - fMinX)/(fNout - 1);
132  } else {
133  fNout = nout;
134  index = new Int_t[nout];
135  TMath::Sort(nout, xout, index, kFALSE);
136  }
137 
138  fGout = new TGraph(fNout);
139  for (Int_t i=0;i<fNout;i++) {
140  if (xout == 0) fGout->SetPoint(i,fMinX + i*delta, 0);
141  else fGout->SetPoint(i,xout[index[i]], 0);
142  }
143 
144  BDRksmooth(fGin->GetX(), fGin->GetY(), fNin, fGout->GetX(),
145  fGout->GetY(), fNout, kernel, bandwidth);
146 
147  if (index) {delete [] index; index = 0;}
148 
149  return fGout;
150 }
151 
152 ////////////////////////////////////////////////////////////////////////////////
153 /// Smooth data with specified kernel.
154 /// Based on R function ksmooth: Translated to C++ by C. Stratowa
155 /// (R source file: ksmooth.c by B.D.Ripley Copyright (C) 1998)
156 
158  Double_t *yp, Int_t np, Int_t kernel, Double_t bw)
159 {
160  Int_t imin = 0;
161  Double_t cutoff = 0.0;
162 
163 // bandwidth is in units of half inter-quartile range
164  if (kernel == 1) {
165  bw *= 0.5;
166  cutoff = bw;
167  }
168  if (kernel == 2) {
169  bw *= 0.3706506;
170  cutoff = 4*bw;
171  }
172 
173  while ((x[imin] < xp[0] - cutoff) && (imin < n)) imin++;
174 
175  for (Int_t j=0;j<np;j++) {
176  Double_t xx, w;
177  Double_t num = 0.0;
178  Double_t den = 0.0;
179  Double_t x0 = xp[j];
180  for (Int_t i=imin;i<n;i++) {
181  if (x[i] < x0 - cutoff) imin = i;
182  if (x[i] > x0 + cutoff) break;
183  xx = TMath::Abs(x[i] - x0)/bw;
184  if (kernel == 1) w = 1;
185  else w = TMath::Exp(-0.5*xx*xx);
186  num += w*y[i];
187  den += w;
188  }
189  if (den > 0) {
190  yp[j] = num/den;
191  } else {
192  yp[j] = 0; //should be NA_REAL (see R.h) or nan("NAN")
193  }
194  }
195 }
196 
197 
198 ////////////////////////////////////////////////////////////////////////////////
199 /// Smooth data with Lowess smoother
200 ///
201 /// This function performs the computations for the LOWESS smoother
202 /// (see the reference below). Lowess returns the output points
203 /// x and y which give the coordinates of the smooth.
204 ///
205 /// \param[in] grin Input graph
206 /// \param[in] span the smoother span. This gives the proportion of points in the plot
207 /// which influence the smooth at each value. Larger values give more smoothness.
208 /// \param[in] iter the number of robustifying iterations which should be performed.
209 /// Using smaller values of iter will make lowess run faster.
210 /// \param[in] delta values of x which lie within delta of each other replaced by a
211 /// single value in the output from lowess.
212 /// For delta = 0, delta will be calculated.
213 ///
214 /// References:
215 ///
216 /// - Cleveland, W. S. (1979) Robust locally weighted regression and smoothing
217 /// scatterplots. J. Amer. Statist. Assoc. 74, 829-836.
218 /// - Cleveland, W. S. (1981) LOWESS: A program for smoothing scatterplots
219 /// by robust locally weighted regression.
220 /// The American Statistician, 35, 54.
221 
223  Double_t span, Int_t iter, Double_t delta)
224 {
225  TString opt = option;
226  opt.ToLower();
227 
228  Smoothin(grin);
229 
230  if (delta == 0) {delta = 0.01*(TMath::Abs(fMaxX - fMinX));}
231 
232 // output X, Y
233  fNout = fNin;
234  fGout = new TGraphErrors(fNout);
235 
236  for (Int_t i=0;i<fNout;i++) {
237  fGout->SetPoint(i,fGin->GetX()[i], 0);
238  }
239 
240  Lowess(fGin->GetX(), fGin->GetY(), fNin, fGout->GetY(), span, iter, delta);
241 
242  return fGout;
243 }
244 
245 ////////////////////////////////////////////////////////////////////////////////
246 /// Lowess regression smoother.
247 /// Based on R function clowess: Translated to C++ by C. Stratowa
248 /// (R source file: lowess.c by R Development Core Team (C) 1999-2001)
249 
251  Double_t span, Int_t iter, Double_t delta)
252 {
253  Int_t i, iiter, j, last, m1, m2, nleft, nright, ns;
254  Double_t alpha, c1, c9, cmad, cut, d1, d2, denom, r;
255  Bool_t ok;
256 
257  if (n < 2) {
258  ys[0] = y[0];
259  return;
260  }
261 
262 // nleft, nright, last, etc. must all be shifted to get rid of these:
263  x--;
264  y--;
265  ys--;
266 
267  Double_t *rw = ((TGraphErrors*)fGout)->GetEX();
268  Double_t *res = ((TGraphErrors*)fGout)->GetEY();
269 
270 // at least two, at most n poInt_ts
271  ns = TMath::Max(2, TMath::Min(n, (Int_t)(span*n + 1e-7)));
272 
273 // robustness iterations
274  iiter = 1;
275  while (iiter <= iter+1) {
276  nleft = 1;
277  nright = ns;
278  last = 0; // index of prev estimated poInt_t
279  i = 1; // index of current poInt_t
280 
281  for(;;) {
282  if (nright < n) {
283  // move nleft, nright to right if radius decreases
284  d1 = x[i] - x[nleft];
285  d2 = x[nright+1] - x[i];
286 
287  // if d1 <= d2 with x[nright+1] == x[nright], lowest fixes
288  if (d1 > d2) {
289  // radius will not decrease by move right
290  nleft++;
291  nright++;
292  continue;
293  }
294  }
295 
296  // fitted value at x[i]
297  Bool_t iterg1 = iiter>1;
298  Lowest(&x[1], &y[1], n, x[i], ys[i], nleft, nright,
299  res, iterg1, rw, ok);
300  if (!ok) ys[i] = y[i];
301 
302  // all weights zero copy over value (all rw==0)
303  if (last < i-1) {
304  denom = x[i]-x[last];
305 
306  // skipped poInt_ts -- Int_terpolate non-zero - proof?
307  for(j = last+1; j < i; j++) {
308  alpha = (x[j]-x[last])/denom;
309  ys[j] = alpha*ys[i] + (1.-alpha)*ys[last];
310  }
311  }
312 
313  // last poInt_t actually estimated
314  last = i;
315 
316  // x coord of close poInt_ts
317  cut = x[last] + delta;
318  for (i = last+1; i <= n; i++) {
319  if (x[i] > cut)
320  break;
321  if (x[i] == x[last]) {
322  ys[i] = ys[last];
323  last = i;
324  }
325  }
326  i = TMath::Max(last+1, i-1);
327  if (last >= n)
328  break;
329  }
330 
331  // residuals
332  for(i=0; i < n; i++)
333  res[i] = y[i+1] - ys[i+1];
334 
335  // compute robustness weights except last time
336  if (iiter > iter)
337  break;
338  for(i=0 ; i<n ; i++)
339  rw[i] = TMath::Abs(res[i]);
340 
341  // compute cmad := 6 * median(rw[], n)
342  m1 = n/2;
343  // partial sort, for m1 & m2
344  Psort(rw, n, m1);
345  if(n % 2 == 0) {
346  m2 = n-m1-1;
347  Psort(rw, n, m2);
348  cmad = 3.*(rw[m1]+rw[m2]);
349  } else { /* n odd */
350  cmad = 6.*rw[m1];
351  }
352 
353  c9 = 0.999*cmad;
354  c1 = 0.001*cmad;
355  for(i=0 ; i<n ; i++) {
356  r = TMath::Abs(res[i]);
357  if (r <= c1)
358  rw[i] = 1.;
359  else if (r <= c9)
360  rw[i] = (1.-(r/cmad)*(r/cmad))*(1.-(r/cmad)*(r/cmad));
361  else
362  rw[i] = 0.;
363  }
364  iiter++;
365  }
366 }
367 
368 ////////////////////////////////////////////////////////////////////////////////
369 /// Fit value at x[i]
370 /// Based on R function lowest: Translated to C++ by C. Stratowa
371 /// (R source file: lowess.c by R Development Core Team (C) 1999-2001)
372 
374  Double_t &ys, Int_t nleft, Int_t nright, Double_t *w,
375  Bool_t userw, Double_t *rw, Bool_t &ok)
376 {
377  Int_t nrt, j;
378  Double_t a, b, c, d, h, h1, h9, r, range;
379 
380  x--;
381  y--;
382  w--;
383  rw--;
384 
385  range = x[n]-x[1];
386  h = TMath::Max(xs-x[nleft], x[nright]-xs);
387  h9 = 0.999*h;
388  h1 = 0.001*h;
389 
390 // sum of weights
391  a = 0.;
392  j = nleft;
393  while (j <= n) {
394  // compute weights (pick up all ties on right)
395  w[j] = 0.;
396  r = TMath::Abs(x[j] - xs);
397  if (r <= h9) {
398  if (r <= h1) {
399  w[j] = 1.;
400  } else {
401  d = (r/h)*(r/h)*(r/h);
402  w[j] = (1.- d)*(1.- d)*(1.- d);
403  }
404  if (userw)
405  w[j] *= rw[j];
406  a += w[j];
407  } else if (x[j] > xs)
408  break;
409  j = j+1;
410  }
411 
412 // rightmost pt (may be greater than nright because of ties)
413  nrt = j-1;
414  if (a <= 0.)
415  ok = kFALSE;
416  else {
417  ok = kTRUE;
418  // weighted least squares: make sum of w[j] == 1
419  for(j=nleft ; j<=nrt ; j++)
420  w[j] /= a;
421  if (h > 0.) {
422  a = 0.;
423  // use linear fit weighted center of x values
424  for(j=nleft ; j<=nrt ; j++)
425  a += w[j] * x[j];
426  b = xs - a;
427  c = 0.;
428  for(j=nleft ; j<=nrt ; j++)
429  c += w[j]*(x[j]-a)*(x[j]-a);
430  if (TMath::Sqrt(c) > 0.001*range) {
431  b /= c;
432  // poInt_ts are spread out enough to compute slope
433  for(j=nleft; j <= nrt; j++)
434  w[j] *= (b*(x[j]-a) + 1.);
435  }
436  }
437  ys = 0.;
438  for(j=nleft; j <= nrt; j++)
439  ys += w[j] * y[j];
440  }
441 }
442 
443 ////////////////////////////////////////////////////////////////////////////////
444 /// Smooth data with Super smoother.
445 /// Smooth the (x, y) values by Friedman's ``super smoother''.
446 ///
447 /// \param[in] grin graph for smoothing
448 /// \param[in] span the fraction of the observations in the span of the running lines
449 /// smoother, or 0 to choose this by leave-one-out cross-validation.
450 /// \param[in] bass controls the smoothness of the fitted curve.
451 /// Values of up to 10 indicate increasing smoothness.
452 /// \param[in] isPeriodic if TRUE, the x values are assumed to be in [0, 1]
453 /// and of period 1.
454 /// \param[in] w case weights
455 ///
456 /// Details:
457 ///
458 /// supsmu is a running lines smoother which chooses between three spans for
459 /// the lines. The running lines smoothers are symmetric, with k/2 data points
460 /// each side of the predicted point, and values of k as 0.5 * n, 0.2 * n and
461 /// 0.05 * n, where n is the number of data points. If span is specified,
462 /// a single smoother with span span * n is used.
463 ///
464 /// The best of the three smoothers is chosen by cross-validation for each
465 /// prediction. The best spans are then smoothed by a running lines smoother
466 /// and the final prediction chosen by linear interpolation.
467 ///
468 /// The FORTRAN code says: ``For small samples (n < 40) or if there are
469 /// substantial serial correlations between observations close in x - value,
470 /// then a prespecified fixed span smoother (span > 0) should be used.
471 /// Reasonable span values are 0.2 to 0.4.''
472 ///
473 /// References:
474 /// - Friedman, J. H. (1984) SMART User's Guide.
475 /// Laboratory for Computational Statistics,
476 /// Stanford University Technical Report No. 1.
477 /// - Friedman, J. H. (1984) A variable span scatterplot smoother.
478 /// Laboratory for Computational Statistics,
479 /// Stanford University Technical Report No. 5.
480 
482  Double_t bass, Double_t span, Bool_t isPeriodic, Double_t *w)
483 {
484  if (span < 0 || span > 1) {
485  std::cout << "Error: Span must be between 0 and 1" << std::endl;
486  return 0;
487  }
488 
489  Smoothin(grin);
490 
491  Int_t iper = 1;
492  if (isPeriodic) {
493  iper = 2;
494  if (fMinX < 0 || fMaxX > 1) {
495  std::cout << "Error: x must be between 0 and 1 for periodic smooth" << std::endl;
496  return 0;
497  }
498  }
499 
500 // output X, Y
501  fNout = fNin;
502  fGout = new TGraph(fNout);
503  Int_t i;
504  for (i=0; i<fNout; i++) {
505  fGout->SetPoint(i,fGin->GetX()[i], 0);
506  }
507 
508 // weights
509  Double_t *weight = new Double_t[fNin];
510  for (i=0; i<fNin; i++) {
511  if (w == 0) weight[i] = 1;
512  else weight[i] = w[i];
513  }
514 
515 // temporary storage array
516  Int_t nTmp = (fNin+1)*8;
517  Double_t *tmp = new Double_t[nTmp];
518  for (i=0; i<nTmp; i++) {
519  tmp[i] = 0;
520  }
521 
522  BDRsupsmu(fNin, fGin->GetX(), fGin->GetY(), weight, iper, span, bass, fGout->GetY(), tmp);
523 
524  delete [] tmp;
525  delete [] weight;
526 
527  return fGout;
528 }
529 
530 ////////////////////////////////////////////////////////////////////////////////
531 /// Friedmanns super smoother (Friedman, 1984).
532 ///
533 /// version 10/10/84
534 /// coded and copywrite (c) 1984 by:
535 ///
536 /// Jerome H. Friedman
537 /// department of statistics
538 /// and
539 /// stanford linear accelerator center
540 /// stanford university
541 ///
542 /// all rights reserved.
543 ///
544 /// \param[in] n number of observations (x,y - pairs).
545 /// \param[in] x ordered abscissa values.
546 /// \param[in] y corresponding ordinate (response) values.
547 /// \param[in] w weight for each (x,y) observation.
548 /// \param[in] iper periodic variable flag.
549 /// - iper=1 => x is ordered interval variable.
550 /// - iper=2 => x is a periodic variable with values
551 /// in the range (0.0,1.0) and period 1.0.
552 /// \param[in] span smoother span (fraction of observations in window).
553 /// - span=0.0 => automatic (variable) span selection.
554 /// \param[in] alpha controls high frequency (small span) penality
555 /// used with automatic span selection (bass tone control).
556 /// (alpha.le.0.0 or alpha.gt.10.0 => no effect.)
557 /// \param[out] smo smoothed ordinate (response) values.
558 /// \param sc internal working storage.
559 ///
560 /// note:
561 ///
562 /// for small samples (n < 40) or if there are substantial serial
563 /// correlations between observations close in x - value, then
564 /// a prespecified fixed span smoother (span > 0) should be
565 /// used. reasonable span values are 0.2 to 0.4.
566 ///
567 /// current implementation:
568 ///
569 /// Based on R function supsmu: Translated to C++ by C. Stratowa
570 /// (R source file: ppr.f by B.D.Ripley Copyright (C) 1994-97)
571 
573  Int_t iper, Double_t span, Double_t alpha, Double_t *smo, Double_t *sc)
574 {
575 // Local variables
576  Int_t sc_offset;
577  Int_t i, j, jper;
578  Double_t a, f;
579  Double_t sw, sy, resmin, vsmlsq;
580  Double_t scale;
581  Double_t d1, d2;
582 
583  Double_t spans[3] = { 0.05, 0.2, 0.5 };
584  Double_t big = 1e20;
585  Double_t sml = 1e-7;
586  Double_t eps = 0.001;
587 
588 // Parameter adjustments
589  sc_offset = n + 1;
590  sc -= sc_offset;
591  --smo;
592  --w;
593  --y;
594  --x;
595 
596 // Function Body
597  if (x[n] <= x[1]) {
598  sy = 0.0;
599  sw = sy;
600  for (j=1;j<=n;++j) {
601  sy += w[j] * y[j];
602  sw += w[j];
603  }
604 
605  a = 0.0;
606  if (sw > 0.0) a = sy / sw;
607  for (j=1;j<=n;++j) smo[j] = a;
608  return;
609  }
610 
611  i = (Int_t)(n / 4);
612  j = i * 3;
613  scale = x[j] - x[i];
614  while (scale <= 0.0) {
615  if (j < n) ++j;
616  if (i > 1) --i;
617  scale = x[j] - x[i];
618  }
619 
620 // Computing 2nd power
621  d1 = eps * scale;
622  vsmlsq = d1 * d1;
623  jper = iper;
624  if (iper == 2 && (x[1] < 0.0 || x[n] > 1.0)) {
625  jper = 1;
626  }
627  if (jper < 1 || jper > 2) {
628  jper = 1;
629  }
630  if (span > 0.0) {
631  BDRsmooth(n, &x[1], &y[1], &w[1], span, jper, vsmlsq,
632  &smo[1], &sc[sc_offset]);
633  return;
634  }
635 
636  Double_t *h = new Double_t[n+1];
637  for (i = 1; i <= 3; ++i) {
638  BDRsmooth(n, &x[1], &y[1], &w[1], spans[i - 1], jper, vsmlsq,
639  &sc[((i<<1)-1)*n + 1], &sc[n*7 + 1]);
640  BDRsmooth(n, &x[1], &sc[n*7 + 1], &w[1], spans[1], -jper, vsmlsq,
641  &sc[(i<<1)*n + 1], &h[1]);
642  }
643 
644  for (j=1; j<=n; ++j) {
645  resmin = big;
646  for (i=1; i<=3; ++i) {
647  if (sc[j + (i<<1)*n] < resmin) {
648  resmin = sc[j + (i<<1)*n];
649  sc[j + n*7] = spans[i-1];
650  }
651  }
652 
653  if (alpha>0.0 && alpha<=10.0 && resmin<sc[j + n*6] && resmin>0.0) {
654  // Computing MAX
655  d1 = TMath::Max(sml,(resmin/sc[j + n*6]));
656  d2 = 10. - alpha;
657  sc[j + n*7] += (spans[2] - sc[j + n*7]) * TMath::Power(d1, d2);
658  }
659  }
660 
661  BDRsmooth(n, &x[1], &sc[n*7 + 1], &w[1], spans[1], -jper, vsmlsq,
662  &sc[(n<<1) + 1], &h[1]);
663 
664  for (j=1; j<=n; ++j) {
665  if (sc[j + (n<<1)] <= spans[0]) {
666  sc[j + (n<<1)] = spans[0];
667  }
668  if (sc[j + (n<<1)] >= spans[2]) {
669  sc[j + (n<<1)] = spans[2];
670  }
671  f = sc[j + (n<<1)] - spans[1];
672  if (f < 0.0) {
673  f = -f / (spans[1] - spans[0]);
674  sc[j + (n<<2)] = (1.0 - f) * sc[j + n*3] + f * sc[j + n];
675  } else {
676  f /= spans[2] - spans[1];
677  sc[j + (n<<2)] = (1.0 - f) * sc[j + n*3] + f * sc[j + n*5];
678  }
679  }
680 
681  BDRsmooth(n, &x[1], &sc[(n<<2) + 1], &w[1], spans[0], -jper, vsmlsq,
682  &smo[1], &h[1]);
683 
684  delete [] h;
685  return;
686 }
687 
688 ////////////////////////////////////////////////////////////////////////////////
689 /// Function for super smoother
690 /// Based on R function supsmu: Translated to C++ by C. Stratowa
691 /// (R source file: ppr.f by B.D.Ripley Copyright (C) 1994-97)
692 
694  Double_t span, Int_t iper, Double_t vsmlsq, Double_t *smo, Double_t *acvr)
695 {
696 // Local variables
697  Int_t i, j, j0, in, out, it, jper, ibw;
698  Double_t a, h1, d1;
699  Double_t xm, ym, wt, sy, fbo, fbw;
700  Double_t cvar, var, tmp, xti, xto;
701 
702 // Parameter adjustments
703  --acvr;
704  --smo;
705  --w;
706  --y;
707  --x;
708 
709 // Function Body
710  xm = 0.;
711  ym = xm;
712  var = ym;
713  cvar = var;
714  fbw = cvar;
715  jper = TMath::Abs(iper);
716 
717  ibw = (Int_t)(span * 0.5 * n + 0.5);
718  if (ibw < 2) {
719  ibw = 2;
720  }
721 
722  it = 2*ibw + 1;
723  for (i=1; i<=it; ++i) {
724  j = i;
725  if (jper == 2) {
726  j = i - ibw - 1;
727  }
728  xti = x[j];
729  if (j < 1) {
730  j = n + j;
731  xti = x[j] - 1.0;
732  }
733  wt = w[j];
734  fbo = fbw;
735  fbw += wt;
736  if (fbw > 0.0) {
737  xm = (fbo * xm + wt * xti) / fbw;
738  ym = (fbo * ym + wt * y[j]) / fbw;
739  }
740  tmp = 0.0;
741  if (fbo > 0.0) {
742  tmp = fbw * wt * (xti - xm) / fbo;
743  }
744  var += tmp * (xti - xm);
745  cvar += tmp * (y[j] - ym);
746  }
747 
748  for (j=1; j<=n; ++j) {
749  out = j - ibw - 1;
750  in = j + ibw;
751  if (!(jper != 2 && (out < 1 || in > n))) {
752  if (out < 1) {
753  out = n + out;
754  xto = x[out] - 1.0;
755  xti = x[in];
756  } else if (in > n) {
757  in -= n;
758  xti = x[in] + 1.0;
759  xto = x[out];
760  } else {
761  xto = x[out];
762  xti = x[in];
763  }
764 
765  wt = w[out];
766  fbo = fbw;
767  fbw -= wt;
768  tmp = 0.0;
769  if (fbw > 0.0) {
770  tmp = fbo * wt * (xto - xm) / fbw;
771  }
772  var -= tmp * (xto - xm);
773  cvar -= tmp * (y[out] - ym);
774  if (fbw > 0.0) {
775  xm = (fbo * xm - wt * xto) / fbw;
776  ym = (fbo * ym - wt * y[out]) / fbw;
777  }
778  wt = w[in];
779  fbo = fbw;
780  fbw += wt;
781  if (fbw > 0.0) {
782  xm = (fbo * xm + wt * xti) / fbw;
783  ym = (fbo * ym + wt * y[in]) / fbw;
784  }
785  tmp = 0.0;
786  if (fbo > 0.0) {
787  tmp = fbw * wt * (xti - xm) / fbo;
788  }
789  var += tmp * (xti - xm);
790  cvar += tmp * (y[in] - ym);
791  }
792 
793  a = 0.0;
794  if (var > vsmlsq) {
795  a = cvar / var;
796  }
797  smo[j] = a * (x[j] - xm) + ym;
798 
799  if (iper <= 0) {
800  continue;
801  }
802 
803  h1 = 0.0;
804  if (fbw > 0.0) {
805  h1 = 1.0 / fbw;
806  }
807  if (var > vsmlsq) {
808  // Computing 2nd power
809  d1 = x[j] - xm;
810  h1 += d1 * d1 / var;
811  }
812 
813  acvr[j] = 0.0;
814  a = 1.0 - w[j] * h1;
815  if (a > 0.0) {
816  acvr[j] = TMath::Abs(y[j] - smo[j]) / a;
817  continue;
818  }
819  if (j > 1) {
820  acvr[j] = acvr[j-1];
821  }
822  }
823 
824  j = 1;
825  do {
826  j0 = j;
827  sy = smo[j] * w[j];
828  fbw = w[j];
829  if (j < n) {
830  do {
831  if (x[j + 1] > x[j]) {
832  break;
833  }
834  ++j;
835  sy += w[j] * smo[j];
836  fbw += w[j];
837  } while (j < n);
838  }
839 
840  if (j > j0) {
841  a = 0.0;
842  if (fbw > 0.0) {
843  a = sy / fbw;
844  }
845  for (i=j0; i<=j; ++i) {
846  smo[i] = a;
847  }
848  }
849  ++j;
850  } while (j <= n);
851 
852  return;
853 }
854 
855 ////////////////////////////////////////////////////////////////////////////////
856 /// Sort data points and eliminate double x values
857 
858 void TGraphSmooth::Approxin(TGraph *grin, Int_t /*iKind*/, Double_t &ylow,
859  Double_t &yhigh, Int_t rule, Int_t iTies)
860 {
861  if (fGout) {delete fGout; fGout = 0;}
862  fGin = grin;
863 
864  fNin = fGin->GetN();
865  Double_t *xin = new Double_t[fNin];
866  Double_t *yin = new Double_t[fNin];
867  Int_t i;
868  for (i=0;i<fNin;i++) {
869  xin[i] = fGin->GetX()[i];
870  yin[i] = fGin->GetY()[i];
871  }
872 
873 // sort/rank input x, y
874  Int_t *index = new Int_t[fNin];
875  Int_t *rank = new Int_t[fNin];
876  Rank(fNin, xin, index, rank, kFALSE);
877 
878 // input X, Y
879  Int_t vNDup = 0;
880  Int_t k = 0;
881  Int_t *dup = new Int_t[fNin];
882  Double_t *x = new Double_t[fNin];
883  Double_t *y = new Double_t[fNin];
884  Double_t vMean, vMin, vMax;
885  for (i=1;i<fNin+1;i++) {
886  Int_t ndup = 1;
887  vMin = vMean = vMax = yin[index[i-1]];
888  while ((i < fNin) && (rank[index[i]] == rank[index[i-1]])) {
889  vMean += yin[index[i]];
890  vMax = (vMax < yin[index[i]]) ? yin[index[i]] : vMax;
891  vMin = (vMin > yin[index[i]]) ? yin[index[i]] : vMin;
892  dup[vNDup] = i;
893  i++;
894  ndup++;
895  vNDup++;
896  }
897  x[k] = xin[index[i-1]];
898  if (ndup == 1) {y[k++] = yin[index[i-1]];}
899  else switch(iTies) {
900  case 1:
901  y[k++] = vMean/ndup;
902  break;
903  case 2:
904  y[k++] = vMin;
905  break;
906  case 3:
907  y[k++] = vMax;
908  break;
909  default:
910  y[k++] = vMean/ndup;
911  break;
912  }
913  }
914  fNin = k;
915 
916 // set unique sorted input data x,y as final graph points
917  fGin->Set(fNin);
918  for (i=0;i<fNin;i++) {
919  fGin->SetPoint(i, x[i], y[i]);
920  }
921 
922  fMinX = fGin->GetX()[0]; //already sorted!
923  fMaxX = fGin->GetX()[fNin-1];
924 
925 // interpolate outside interval [min(x),max(x)]
926  switch(rule) {
927  case 1:
928  ylow = 0; // = nan("NAN") ??
929  yhigh = 0; // = nan("NAN") ??
930  break;
931  case 2:
932  ylow = fGin->GetY()[0];
933  yhigh = fGin->GetY()[fNin-1];
934  break;
935  default:
936  break;
937  }
938 
939 // cleanup
940  delete [] x;
941  delete [] y;
942  delete [] dup;
943  delete [] rank;
944  delete [] index;
945  delete [] xin;
946  delete [] yin;
947 }
948 
949 ////////////////////////////////////////////////////////////////////////////////
950 /// Approximate data points
951 /// \param[in] grin graph giving the coordinates of the points to be interpolated.
952 /// Alternatively a single plotting structure can be specified:
953 /// \param[in] option specifies the interpolation method to be used.
954 /// Choices are "linear" (iKind = 1) or "constant" (iKind = 2).
955 /// \param[in] nout If xout is not specified, interpolation takes place at n equally
956 /// spaced points spanning the interval [min(x), max(x)], where
957 /// nout = max(nout, number of input data).
958 /// \param[in] xout an optional set of values specifying where interpolation is to
959 /// take place.
960 /// \param[in] yleft the value to be returned when input x values less than min(x).
961 /// The default is defined by the value of rule given below.
962 /// \param[in] yright the value to be returned when input x values greater than max(x).
963 /// The default is defined by the value of rule given below.
964 /// \param[in] rule an integer describing how interpolation is to take place outside
965 /// the interval [min(x), max(x)]. If rule is 0 then the given yleft
966 /// and yright values are returned, if it is 1 then 0 is returned
967 /// for such points and if it is 2, the value at the closest data
968 /// extreme is used.
969 /// \param[in] f For method="constant" a number between 0 and 1 inclusive,
970 /// indicating a compromise between left- and right-continuous step
971 /// functions. If y0 and y1 are the values to the left and right of
972 /// the point then the value is y0*f+y1*(1-f) so that f=0 is
973 /// right-continuous and f=1 is left-continuous
974 /// \param[in] ties Handling of tied x values. An integer describing a function with
975 /// a single vector argument returning a single number result:
976 /// - ties = "ordered" (iTies = 0): input x are "ordered"
977 /// - ties = "mean" (iTies = 1): function "mean"
978 /// - ties = "min" (iTies = 2): function "min"
979 /// - ties = "max" (iTies = 3): function "max"
980 ///
981 /// Details:
982 ///
983 /// At least two complete (x, y) pairs are required.
984 /// If there are duplicated (tied) x values and ties is a function it is
985 /// applied to the y values for each distinct x value. Useful functions in
986 /// this context include mean, min, and max.
987 /// If ties="ordered" the x values are assumed to be already ordered. The
988 /// first y value will be used for interpolation to the left and the last
989 /// one for interpolation to the right.
990 ///
991 /// Value:
992 ///
993 /// approx returns a graph with components x and y, containing n coordinates
994 /// which interpolate the given data points according to the method (and rule)
995 /// desired.
996 
998  Double_t yleft, Double_t yright, Int_t rule, Double_t f, Option_t *ties)
999 {
1000  TString opt = option;
1001  opt.ToLower();
1002  Int_t iKind = 0;
1003  if (opt.Contains("linear")) iKind = 1;
1004  else if (opt.Contains("constant")) iKind = 2;
1005 
1006  if (f < 0 || f > 1) {
1007  std::cout << "Error: Invalid f value" << std::endl;
1008  return 0;
1009  }
1010 
1011  opt = ties;
1012  opt.ToLower();
1013  Int_t iTies = 0;
1014  if (opt.Contains("ordered")) {
1015  iTies = 0;
1016  } else if (opt.Contains("mean")) {
1017  iTies = 1;
1018  } else if (opt.Contains("min")) {
1019  iTies = 2;
1020  } else if (opt.Contains("max")) {
1021  iTies = 3;
1022  } else {
1023  std::cout << "Error: Method not known: " << ties << std::endl;
1024  return 0;
1025  }
1026 
1027 // input X, Y
1028  Double_t ylow = yleft;
1029  Double_t yhigh = yright;
1030  Approxin(grin, iKind, ylow, yhigh, rule, iTies);
1031 
1032 // output X, Y
1033  Double_t delta = 0;
1034  fNout = nout;
1035  if (xout == 0) {
1036  fNout = TMath::Max(nout, fNin);
1037  delta = (fMaxX - fMinX)/(fNout - 1);
1038  }
1039 
1040  fGout = new TGraph(fNout);
1041 
1042  Double_t x;
1043  for (Int_t i=0;i<fNout;i++) {
1044  if (xout == 0) x = fMinX + i*delta;
1045  else x = xout[i];
1046  Double_t yout = Approx1(x, f, fGin->GetX(), fGin->GetY(), fNin, iKind, ylow, yhigh);
1047  fGout->SetPoint(i,x, yout);
1048  }
1049 
1050  return fGout;
1051 }
1052 
1053 ////////////////////////////////////////////////////////////////////////////////
1054 /// Approximate one data point.
1055 /// Approximate y(v), given (x,y)[i], i = 0,..,n-1
1056 /// Based on R function approx1: Translated to C++ by Christian Stratowa
1057 /// (R source file: approx.c by R Development Core Team (C) 1999-2001)
1058 
1060  Int_t n, Int_t iKind, Double_t ylow, Double_t yhigh)
1061 {
1062  Int_t i = 0;
1063  Int_t j = n - 1;
1064 
1065 // handle out-of-domain points
1066  if(v < x[i]) return ylow;
1067  if(v > x[j]) return yhigh;
1068 
1069 // find the correct interval by bisection
1070  while(i < j - 1) {
1071  Int_t ij = (i + j)/2;
1072  if(v < x[ij]) j = ij;
1073  else i = ij;
1074  }
1075 
1076 // interpolation
1077  if(v == x[j]) return y[j];
1078  if(v == x[i]) return y[i];
1079 
1080  if(iKind == 1) { // linear
1081  return y[i] + (y[j] - y[i]) * ((v - x[i])/(x[j] - x[i]));
1082  } else { // 2 : constant
1083  return y[i] * (1-f) + y[j] * f;
1084  }
1085 }
1086 
1087 // helper functions
1088 ////////////////////////////////////////////////////////////////////////////////
1089 /// Static function
1090 /// if (ISNAN(x)) return 1;
1091 /// if (ISNAN(y)) return -1;
1092 
1094 {
1095  if (x < y) return -1;
1096  if (x > y) return 1;
1097  return 0;
1098 }
1099 
1100 ////////////////////////////////////////////////////////////////////////////////
1101 /// Static function
1102 /// based on R function rPsort: adapted to C++ by Christian Stratowa
1103 /// (R source file: R_sort.c by R Development Core Team (C) 1999-2001)
1104 
1106 {
1107  Double_t v, w;
1108  Int_t pL, pR, i, j;
1109 
1110  for (pL = 0, pR = n - 1; pL < pR; ) {
1111  v = x[k];
1112  for(i = pL, j = pR; i <= j;) {
1113  while (TGraphSmooth::Rcmp(x[i], v) < 0) i++;
1114  while (TGraphSmooth::Rcmp(v, x[j]) < 0) j--;
1115  if (i <= j) { w = x[i]; x[i++] = x[j]; x[j--] = w; }
1116  }
1117  if (j < k) pL = i;
1118  if (k < i) pR = j;
1119  }
1120 }
1121 
1122 ////////////////////////////////////////////////////////////////////////////////
1123 /// static function
1124 
1125 void TGraphSmooth::Rank(Int_t n, Double_t *a, Int_t *index, Int_t *rank, Bool_t down)
1126 {
1127  if (n <= 0) return;
1128  if (n == 1) {
1129  index[0] = 0;
1130  rank[0] = 0;
1131  return;
1132  }
1133 
1134  TMath::Sort(n,a,index,down);
1135 
1136  Int_t k = 0;
1137  for (Int_t i=0;i<n;i++) {
1138  if ((i > 0) && (a[index[i]] == a[index[i-1]])) {
1139  rank[index[i]] = i-1;
1140  k++;
1141  }
1142  rank[index[i]] = i-k;
1143  }
1144 }
static void 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)
Friedmanns super smoother (Friedman, 1984).
TGraph * SmoothKern(TGraph *grin, Option_t *option="normal", Double_t bandwidth=0.5, Int_t nout=100, Double_t *xout=0)
Smooth data with Kernel smoother.
virtual ~TGraphSmooth()
GraphSmooth destructor.
const char Option_t
Definition: RtypesCore.h:62
return c1
Definition: legend1.C:41
void 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.
TH1 * h
Definition: legend2.C:5
Double_t fMaxX
Definition: TGraphSmooth.h:48
Basic string class.
Definition: TString.h:129
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:168
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1099
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
TArc * a
Definition: textangle.C:12
A helper class to smooth TGraph.
Definition: TGraphSmooth.h:36
TGraph * SmoothSuper(TGraph *grin, Option_t *option="", Double_t bass=0, Double_t span=0, Bool_t isPeriodic=kFALSE, Double_t *w=0)
Smooth data with Super smoother.
static void 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.
TGraph * fGin
Definition: TGraphSmooth.h:45
Float_t delta
Short_t Abs(Short_t d)
Definition: TMathBase.h:108
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:628
Double_t x[n]
Definition: legend1.C:17
TGraph * SmoothLowess(TGraph *grin, Option_t *option="", Double_t span=0.67, Int_t iter=3, Double_t delta=0)
Smooth data with Lowess smoother.
TGraph * fGout
Definition: TGraphSmooth.h:46
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
void Sort(Index n, const Element *a, Index *index, Bool_t down=kTRUE)
Definition: TMath.h:1151
TH1F * h1
Definition: legend1.C:5
void 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.
static Int_t Rcmp(Double_t x, Double_t y)
Static function if (ISNAN(x)) return 1; if (ISNAN(y)) return -1;.
TRandom2 r(17)
SVector< double, 2 > v
Definition: Dict.h:5
static void 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.
Int_t GetN() const
Definition: TGraph.h:122
const Bool_t kFALSE
Definition: RtypesCore.h:92
Double_t * GetX() const
Definition: TGraph.h:129
Double_t Exp(Double_t x)
Definition: TMath.h:622
#define ClassImp(name)
Definition: Rtypes.h:336
double f(double x)
double Double_t
Definition: RtypesCore.h:55
Double_t y[n]
Definition: legend1.C:17
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:572
static void Rank(Int_t n, Double_t *a, Int_t *index, Int_t *rank, Bool_t down=kTRUE)
static function
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
TGraph * Approx(TGraph *grin, Option_t *option="linear", Int_t nout=50, Double_t *xout=0, Double_t yleft=0, Double_t yright=0, Int_t rule=0, Double_t f=0, Option_t *ties="mean")
Approximate data points.
Double_t fMinX
Definition: TGraphSmooth.h:47
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
Definition: TGraph.cxx:2156
void Smoothin(TGraph *grin)
Sort input data points.
Double_t * GetY() const
Definition: TGraph.h:130
static void 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_so...
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:200
A Graph is a graphics object made of two arrays X and Y with npoints each.
Definition: TGraph.h:41
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
A TGraphErrors is a TGraph with error bars.
Definition: TGraphErrors.h:26
Double_t Sqrt(Double_t x)
Definition: TMath.h:591
virtual void Set(Int_t n)
Set number of points in the graph Existing coordinates are preserved New coordinates above fNpoints a...
Definition: TGraph.cxx:2105
const Bool_t kTRUE
Definition: RtypesCore.h:91
static Double_t 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.
const Int_t n
Definition: legend1.C:16
static void 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.