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