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