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