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
37A helper class to smooth TGraph.
38see 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
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
859void 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
1126void 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}
ROOT::R::TRInterface & r
Definition: Object.C:4
#define d(i)
Definition: RSha256.hxx:102
#define b(i)
Definition: RSha256.hxx:100
#define f(i)
Definition: RSha256.hxx:104
#define c(i)
Definition: RSha256.hxx:101
#define h(i)
Definition: RSha256.hxx:106
#define e(i)
Definition: RSha256.hxx:103
int Int_t
Definition: RtypesCore.h:41
const Bool_t kFALSE
Definition: RtypesCore.h:88
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
const Bool_t kTRUE
Definition: RtypesCore.h:87
const char Option_t
Definition: RtypesCore.h:62
#define ClassImp(name)
Definition: Rtypes.h:365
char name[80]
Definition: TGX11.cxx:109
A TGraphErrors is a TGraph with error bars.
Definition: TGraphErrors.h:26
A helper class to smooth TGraph.
Definition: TGraphSmooth.h:36
Double_t fMinX
Definition: TGraphSmooth.h:47
TGraph * fGin
Definition: TGraphSmooth.h:45
Double_t fMaxX
Definition: TGraphSmooth.h:48
static Int_t Rcmp(Double_t x, Double_t y)
Static function if (ISNAN(x)) return 1; if (ISNAN(y)) return -1;.
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.
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.
virtual ~TGraphSmooth()
GraphSmooth destructor.
static void Rank(Int_t n, Double_t *a, Int_t *index, Int_t *rank, Bool_t down=kTRUE)
static function
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.
void Smoothin(TGraph *grin)
Sort input data points.
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).
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...
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.
TGraph * fGout
Definition: TGraphSmooth.h:46
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 * 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 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.
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.
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.
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.
A Graph is a graphics object made of two arrays X and Y with npoints each.
Definition: TGraph.h:41
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
Definition: TGraph.cxx:2257
Double_t * GetY() const
Definition: TGraph.h:131
Int_t GetN() const
Definition: TGraph.h:123
Double_t * GetX() const
Definition: TGraph.h:130
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:2192
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
Basic string class.
Definition: TString.h:131
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1125
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
return c1
Definition: legend1.C:41
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
TH1F * h1
Definition: legend1.C:5
static constexpr double ns
static constexpr double m2
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
Double_t Exp(Double_t x)
Definition: TMath.h:717
Double_t Sqrt(Double_t x)
Definition: TMath.h:681
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:725
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
void Sort(Index n, const Element *a, Index *index, Bool_t down=kTRUE)
Definition: TMathBase.h:362
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
auto * a
Definition: textangle.C:12