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