Logo ROOT  
Reference Guide
TPrincipal.cxx
Go to the documentation of this file.
1 // @(#)root/hist:$Id$
2 // Author: Christian Holm Christensen 1/8/2000
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2004, 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 /** \class TPrincipal
13  \ingroup Hist
14 Principal Components Analysis (PCA)
15 
16 The current implementation is based on the LINTRA package from CERNLIB
17 by R. Brun, H. Hansroul, and J. Kubler.
18 The class has been implemented by Christian Holm Christensen in August 2000.
19 
20 ## Introduction
21 
22 In many applications of various fields of research, the treatment of
23 large amounts of data requires powerful techniques capable of rapid
24 data reduction and analysis. Usually, the quantities most
25 conveniently measured by the experimentalist, are not necessarily the
26 most significant for classification and analysis of the data. It is
27 then useful to have a way of selecting an optimal set of variables
28 necessary for the recognition process and reducing the dimensionality
29 of the problem, resulting in an easier classification procedure.
30 
31 This paper describes the implementation of one such method of
32 feature selection, namely the principal components analysis. This
33 multidimensional technique is well known in the field of pattern
34 recognition and and its use in Particle Physics has been documented
35 elsewhere (cf. H. Wind, <I>Function Parameterization</I>, CERN
36 72-21).
37 
38 ## Overview
39 Suppose we have prototypes which are trajectories of particles,
40 passing through a spectrometer. If one measures the passage of the
41 particle at say 8 fixed planes, the trajectory is described by an
42 8-component vector:
43 \f[
44  \mathbf{x} = \left(x_0, x_1, \ldots, x_7\right)
45 \f]
46 in 8-dimensional pattern space.
47 
48 One proceeds by generating a a representative tracks sample and
49 building up the covariance matrix \f$\mathsf{C}\f$. Its eigenvectors and
50 eigenvalues are computed by standard methods, and thus a new basis is
51 obtained for the original 8-dimensional space the expansion of the
52 prototypes,
53 \f[
54  \mathbf{x}_m = \sum^7_{i=0} a_{m_i} \mathbf{e}_i
55  \quad
56  \mbox{where}
57  \quad
58  a_{m_i} = \mathbf{x}^T\bullet\mathbf{e}_i
59 \f]
60 allows the study of the behavior of the coefficients \f$a_{m_i}\f$ for all
61 the tracks of the sample. The eigenvectors which are insignificant for
62 the trajectory description in the expansion will have their
63 corresponding coefficients \f$a_{m_i}\f$ close to zero for all the
64 prototypes.
65 
66 On one hand, a reduction of the dimensionality is then obtained by
67 omitting these least significant vectors in the subsequent analysis.
68 
69 On the other hand, in the analysis of real data, these least
70 significant variables(?) can be used for the pattern
71 recognition problem of extracting the valid combinations of
72 coordinates describing a true trajectory from the set of all possible
73 wrong combinations.
74 
75 The program described here performs this principal components analysis
76 on a sample of data provided by the user. It computes the covariance
77 matrix, its eigenvalues ands corresponding eigenvectors and exhibits
78 the behavior of the principal components \f$a_{m_i}\f$, thus providing
79 to the user all the means of understanding their data.
80 
81 ## Principal Components Method
82 Let's consider a sample of \f$M\f$ prototypes each being characterized by
83 \f$P\f$ variables \f$x_0, x_1, \ldots, x_{P-1}\f$. Each prototype is a point, or a
84 column vector, in a \f$P\f$-dimensional *Pattern space*.
85 \f[
86  \mathbf{x} = \left[\begin{array}{c}
87  x_0\\x_1\\\vdots\\x_{P-1}\end{array}\right]\,,
88 \f]
89 where each \f$x_n\f$ represents the particular value associated with the
90 \f$n\f$-dimension.
91 
92 Those \f$P\f$ variables are the quantities accessible to the
93 experimentalist, but are not necessarily the most significant for the
94 classification purpose.
95 
96 The *Principal Components Method* consists of applying a
97 *linear* transformation to the original variables. This
98 transformation is described by an orthogonal matrix and is equivalent
99 to a rotation of the original pattern space into a new set of
100 coordinate vectors, which hopefully provide easier feature
101 identification and dimensionality reduction.
102 
103 Let's define the covariance matrix:
104 \f[
105  \mathsf{C} = \left\langle\mathbf{y}\mathbf{y}^T\right\rangle
106  \quad\mbox{where}\quad
107  \mathbf{y} = \mathbf{x} - \left\langle\mathbf{x}\right\rangle\,,
108 \f]
109 and the brackets indicate mean value over the sample of \f$M\f$
110 prototypes.
111 
112 This matrix \f$\mathsf{C}\f$ is real, positive definite, symmetric, and will
113 have all its eigenvalues greater then zero. It will now be show that
114 among the family of all the complete orthonormal bases of the pattern
115 space, the base formed by the eigenvectors of the covariance matrix
116 and belonging to the largest eigenvalues, corresponds to the most
117 significant features of the description of the original prototypes.
118 
119 let the prototypes be expanded on into a set of \f$N\f$ basis vectors
120 \f$\mathbf{e}_n, n=0,\ldots,N,N+1, \ldots, P-1\f$
121 \f[
122  \mathbf{y}_i = \sum^N_{i=0} a_{i_n} \mathbf{e}_n,
123  \quad
124  i = 1, \ldots, M,
125  \quad
126  N < P-1
127 \f]
128 The `best' feature coordinates \f$\mathbf{e}_n\f$, spanning a *feature
129 space*, will be obtained by minimizing the error due to this
130 truncated expansion, i.e.,
131 \f[
132  \min\left(E_N\right) =
133  \min\left[\left\langle\left(\mathbf{y}_i - \sum^N_{i=0} a_{i_n} \mathbf{e}_n\right)^2\right\rangle\right]
134 \f]
135 with the conditions:
136 \f[
137  \mathbf{e}_k\bullet\mathbf{e}_j = \delta_{jk} =
138  \left\{\begin{array}{rcl}
139  1 & \mbox{for} & k = j\\
140  0 & \mbox{for} & k \neq j
141  \end{array}\right.
142 \f]
143 Multiplying (3) by \f$\mathbf{e}^T_n\f$ using (5),
144 we get
145 \f[
146  a_{i_n} = \mathbf{y}_i^T\bullet\mathbf{e}_n\,,
147 \f]
148 so the error becomes
149 \f{eqnarray*}{
150  E_N &=&
151  \left\langle\left[\sum_{n=N+1}^{P-1} a_{i_n}\mathbf{e}_n\right]^2\right\rangle\nonumber\\
152  &=&
153  \left\langle\left[\sum_{n=N+1}^{P-1} \mathbf{y}_i^T\bullet\mathbf{e}_n\mathbf{e}_n\right]^2\right\rangle\nonumber\\
154  &=&
155  \left\langle\sum_{n=N+1}^{P-1} \mathbf{e}_n^T\mathbf{y}_i\mathbf{y}_i^T\mathbf{e}_n\right\rangle\nonumber\\
156  &=&
157  \sum_{n=N+1}^{P-1} \mathbf{e}_n^T\mathsf{C}\mathbf{e}_n
158 \f}
159 The minimization of the sum in (7) is obtained when each
160 term \f$\mathbf{e}_n^\mathsf{C}\mathbf{e}_n\f$ is minimum, since \f$\mathsf{C}\f$ is
161 positive definite. By the method of Lagrange multipliers, and the
162 condition (5), we get
163 \f[
164  E_N = \sum^{P-1}_{n=N+1} \left(\mathbf{e}_n^T\mathsf{C}\mathbf{e}_n -
165  l_n\mathbf{e}_n^T\bullet\mathbf{e}_n + l_n\right)
166 \f]
167 The minimum condition \f$\frac{dE_N}{d\mathbf{e}^T_n} = 0\f$ leads to the
168 equation
169 \f[
170  \mathsf{C}\mathbf{e}_n = l_n\mathbf{e}_n\,,
171 \f]
172 which shows that \f$\mathbf{e}_n\f$ is an eigenvector of the covariance
173 matrix \f$\mathsf{C}\f$ with eigenvalue \f$l_n\f$. The estimated minimum error is
174 then given by
175 \f[
176  E_N \sim \sum^{P-1}_{n=N+1} \mathbf{e}_n^T\bullet l_n\mathbf{e}_n
177  = \sum^{P-1}_{n=N+1} l_n\,,
178 \f]
179 where \f$l_n,\,n=N+1,\ldots,P\f$ \f$l_n,\,n=N+1,\ldots,P-1\f$ are the eigenvalues associated with the
180 omitted eigenvectors in the expansion (3). Thus, by choosing
181 the \f$N\f$ largest eigenvalues, and their associated eigenvectors, the
182 error \f$E_N\f$ is minimized.
183 
184 The transformation matrix to go from the pattern space to the feature
185 space consists of the ordered eigenvectors \f$\mathbf{e}_1,\ldots,\mathbf{e}_P\f$
186 \f$\mathbf{e}_0,\ldots,\mathbf{e}_{P-1}\f$ for its columns
187 \f[
188  \mathsf{T} = \left[
189  \begin{array}{cccc}
190  \mathbf{e}_0 &
191  \mathbf{e}_1 &
192  \vdots &
193  \mathbf{e}_{P-1}
194  \end{array}\right]
195  = \left[
196  \begin{array}{cccc}
197  \mathbf{e}_{0_0} & \mathbf{e}_{1_0} & \cdots & \mathbf{e}_{{P-1}_0}\\
198  \mathbf{e}_{0_1} & \mathbf{e}_{1_1} & \cdots & \mathbf{e}_{{P-1}_1}\\
199  \vdots & \vdots & \ddots & \vdots \\
200  \mathbf{e}_{0_{P-1}} & \mathbf{e}_{1_{P-1}} & \cdots & \mathbf{e}_{{P-1}_{P-1}}\\
201  \end{array}\right]
202 \f]
203 This is an orthogonal transformation, or rotation, of the pattern
204 space and feature selection results in ignoring certain coordinates
205 in the transformed space.
206 
207 Christian Holm August 2000, CERN
208 */
209 
210 #include "TPrincipal.h"
211 
212 #include "TVectorD.h"
213 #include "TMatrixD.h"
214 #include "TMatrixDSymEigen.h"
215 #include "TMath.h"
216 #include "TList.h"
217 #include "TH2.h"
218 #include "TDatime.h"
219 #include "TBrowser.h"
220 #include "TROOT.h"
221 #include "Riostream.h"
222 
223 
225 
226 ////////////////////////////////////////////////////////////////////////////////
227 /// Empty constructor. Do not use.
228 
230  : fMeanValues(0),
231  fSigmas(0),
232  fCovarianceMatrix(1,1),
233  fEigenVectors(1,1),
234  fEigenValues(0),
235  fOffDiagonal(0),
236  fStoreData(kFALSE)
237 {
238  fTrace = 0;
239  fHistograms = 0;
242  fNumberOfVariables = 0;
243 }
244 
245 ////////////////////////////////////////////////////////////////////////////////
246 /// Constructor. Argument is number of variables in the sample of data
247 /// Options are:
248 /// - N Normalize the covariance matrix (default)
249 /// - D Store input data (default)
250 ///
251 /// The created object is named "principal" by default.
252 
254  : fMeanValues(nVariables),
255  fSigmas(nVariables),
256  fCovarianceMatrix(nVariables,nVariables),
257  fEigenVectors(nVariables,nVariables),
258  fEigenValues(nVariables),
259  fOffDiagonal(nVariables),
260  fStoreData(kFALSE)
261 {
262  if (nVariables <= 1) {
263  Error("TPrincipal", "You can't be serious - nVariables == 1!!!");
264  return;
265  }
266 
267  SetName("principal");
268 
269  fTrace = 0;
270  fHistograms = 0;
273  fNumberOfVariables = nVariables;
274  while (strlen(opt) > 0) {
275  switch(*opt++) {
276  case 'N':
277  case 'n':
279  break;
280  case 'D':
281  case 'd':
282  fStoreData = kTRUE;
283  break;
284  default:
285  break;
286  }
287  }
288 
289  if (!fMeanValues.IsValid())
290  Error("TPrincipal","Couldn't create vector mean values");
291  if (!fSigmas.IsValid())
292  Error("TPrincipal","Couldn't create vector sigmas");
293  if (!fCovarianceMatrix.IsValid())
294  Error("TPrincipal","Couldn't create covariance matrix");
295  if (!fEigenVectors.IsValid())
296  Error("TPrincipal","Couldn't create eigenvector matrix");
297  if (!fEigenValues.IsValid())
298  Error("TPrincipal","Couldn't create eigenvalue vector");
299  if (!fOffDiagonal.IsValid())
300  Error("TPrincipal","Couldn't create offdiagonal vector");
301  if (fStoreData) {
302  fUserData.ResizeTo(nVariables*1000);
303  fUserData.Zero();
304  if (!fUserData.IsValid())
305  Error("TPrincipal","Couldn't create user data vector");
306  }
307 }
308 
309 ////////////////////////////////////////////////////////////////////////////////
310 /// Copy constructor.
311 
313  TNamed(pr),
314  fNumberOfDataPoints(pr.fNumberOfDataPoints),
315  fNumberOfVariables(pr.fNumberOfVariables),
316  fMeanValues(pr.fMeanValues),
317  fSigmas(pr.fSigmas),
318  fCovarianceMatrix(pr.fCovarianceMatrix),
319  fEigenVectors(pr.fEigenVectors),
320  fEigenValues(pr.fEigenValues),
321  fOffDiagonal(pr.fOffDiagonal),
322  fUserData(pr.fUserData),
323  fTrace(pr.fTrace),
324  fHistograms(pr.fHistograms),
325  fIsNormalised(pr.fIsNormalised),
326  fStoreData(pr.fStoreData)
327 {
328 }
329 
330 ////////////////////////////////////////////////////////////////////////////////
331 /// Assignment operator.
332 
334 {
335  if(this!=&pr) {
336  TNamed::operator=(pr);
340  fSigmas=pr.fSigmas;
345  fUserData=pr.fUserData;
346  fTrace=pr.fTrace;
350  }
351  return *this;
352 }
353 
354 ////////////////////////////////////////////////////////////////////////////////
355 /// Destructor.
356 
358 {
359  if (fHistograms) {
360  fHistograms->Delete();
361  delete fHistograms;
362  }
363 }
364 
365 ////////////////////////////////////////////////////////////////////////////////
366 /// Add a data point and update the covariance matrix. The input
367 /// array must be <TT>fNumberOfVariables</TT> long.
368 ///
369 ///
370 /// The Covariance matrix and mean values of the input data is calculated
371 /// on the fly by the following equations:
372 ///
373 /// \f[
374 /// \left<x_i\right>^{(0)} = x_{i0}
375 /// \f]
376 ///
377 ///
378 /// \f[
379 /// \left<x_i\right>^{(n)} = \left<x_i\right>^{(n-1)}
380 /// + \frac1n \left(x_{in} - \left<x_i\right>^{(n-1)}\right)
381 /// \f]
382 ///
383 /// \f[
384 /// C_{ij}^{(0)} = 0
385 /// \f]
386 ///
387 ///
388 ///
389 /// \f[
390 /// C_{ij}^{(n)} = C_{ij}^{(n-1)}
391 /// + \frac1{n-1}\left[\left(x_{in} - \left<x_i\right>^{(n)}\right)
392 /// \left(x_{jn} - \left<x_j\right>^{(n)}\right)\right]
393 /// - \frac1n C_{ij}^{(n-1)}
394 /// \f]
395 ///
396 /// since this is a really fast method, with no rounding errors (please
397 /// refer to CERN 72-21 pp. 54-106).
398 ///
399 ///
400 /// The data is stored internally in a <TT>TVectorD</TT>, in the following
401 /// way:
402 ///
403 /// \f[
404 /// \mathbf{x} = \left[\left(x_{0_0},\ldots,x_{{P-1}_0}\right),\ldots,
405 /// \left(x_{0_i},\ldots,x_{{P-1}_i}\right), \ldots\right]
406 /// \f]
407 ///
408 /// With \f$P\f$ as defined in the class description.
409 
411 {
412  if (!p)
413  return;
414 
415  // Increment the data point counter
416  Int_t i,j;
417  if (++fNumberOfDataPoints == 1) {
418  for (i = 0; i < fNumberOfVariables; i++)
419  fMeanValues(i) = p[i];
420  }
421  else {
422 
423  const Double_t invnp = 1. / Double_t(fNumberOfDataPoints);
424  const Double_t invnpM1 = 1. /(Double_t(fNumberOfDataPoints - 1));
425  const Double_t cor = 1. - invnp;
426  // use directly vector array for faster element access
427  Double_t * meanValues = fMeanValues.GetMatrixArray();
428  Double_t * covMatrix = fCovarianceMatrix.GetMatrixArray();
429  for (i = 0; i < fNumberOfVariables; i++) {
430 
431  meanValues[i] *= cor;
432  meanValues[i] += p[i] * invnp;
433  const Double_t t1 = (p[i] - meanValues[i]) * invnpM1;
434 
435  // Setting Matrix (lower triangle) elements
436  for (j = 0; j < i + 1; j++) {
437  const Int_t index = i * fNumberOfVariables + j;
438  covMatrix[index] *= cor;
439  covMatrix[index] += t1 * (p[j] - meanValues[j]);
440  }
441  }
442  }
443 
444  // Store data point in internal vector
445  // If the vector isn't big enough to hold the new data, then
446  // expand the vector by half it's size.
447  if (!fStoreData)
448  return;
449  Int_t size = fUserData.GetNrows();
451  fUserData.ResizeTo(size + size/2);
452 
453  for (i = 0; i < fNumberOfVariables; i++) {
455  fUserData.GetMatrixArray()[j] = p[i];
456  }
457 
458 }
459 
460 ////////////////////////////////////////////////////////////////////////////////
461 /// Browse the TPrincipal object in the TBrowser.
462 
464 {
465  if (fHistograms) {
466  TIter next(fHistograms);
467  TH1* h = 0;
468  while ((h = (TH1*)next()))
469  b->Add(h,h->GetName());
470  }
471 
472  if (fStoreData)
473  b->Add(&fUserData,"User Data");
474  b->Add(&fCovarianceMatrix,"Covariance Matrix");
475  b->Add(&fMeanValues,"Mean value vector");
476  b->Add(&fSigmas,"Sigma value vector");
477  b->Add(&fEigenValues,"Eigenvalue vector");
478  b->Add(&fEigenVectors,"Eigenvector Matrix");
479 
480 }
481 
482 ////////////////////////////////////////////////////////////////////////////////
483 /// Clear the data in Object. Notice, that's not possible to change
484 /// the dimension of the original data.
485 
487 {
488  if (fHistograms) {
489  fHistograms->Delete(opt);
490  }
491 
493  fTrace = 0;
494  fCovarianceMatrix.Zero();
495  fEigenVectors.Zero();
496  fEigenValues.Zero();
497  fMeanValues.Zero();
498  fSigmas.Zero();
499  fOffDiagonal.Zero();
500 
501  if (fStoreData) {
503  fUserData.Zero();
504  }
505 }
506 
507 ////////////////////////////////////////////////////////////////////////////////
508 /// Return a row of the user supplied data.
509 /// If row is out of bounds, 0 is returned.
510 /// It's up to the user to delete the returned array.
511 /// Row 0 is the first row;
512 
514 {
515  if (row >= fNumberOfDataPoints)
516  return 0;
517 
518  if (!fStoreData)
519  return 0;
520 
521  Int_t index = row * fNumberOfVariables;
522  return &fUserData(index);
523 }
524 
525 
526 ////////////////////////////////////////////////////////////////////////////////
527 /// Generates the file `<filename>`, with `.C` appended if it does
528 /// argument doesn't end in .cxx or .C.
529 ///
530 /// The file contains the implementation of two functions
531 /// ~~~ {.cpp}
532 /// void X2P(Double_t *x, Double *p)
533 /// void P2X(Double_t *p, Double *x, Int_t nTest)
534 /// ~~~
535 /// which does the same as `TPrincipal::X2P` and `TPrincipal::P2X`
536 /// respectively. Please refer to these methods.
537 ///
538 /// Further, the static variables:
539 /// ~~~ {.cpp}
540 /// Int_t gNVariables
541 /// Double_t gEigenValues[]
542 /// Double_t gEigenVectors[]
543 /// Double_t gMeanValues[]
544 /// Double_t gSigmaValues[]
545 /// ~~~
546 /// are initialized. The only ROOT header file needed is Rtypes.h
547 ///
548 /// See TPrincipal::MakeRealCode for a list of options
549 
550 void TPrincipal::MakeCode(const char *filename, Option_t *opt)
551 {
552  TString outName(filename);
553  if (!outName.EndsWith(".C") && !outName.EndsWith(".cxx"))
554  outName += ".C";
555 
556  MakeRealCode(outName.Data(),"",opt);
557 }
558 
559 ////////////////////////////////////////////////////////////////////////////////
560 /// Make histograms of the result of the analysis.
561 /// The option string say which histograms to create
562 /// - X Histogram original data
563 /// - P Histogram principal components corresponding to
564 /// original data
565 /// - D Histogram the difference between the original data
566 /// and the projection of principal unto a lower
567 /// dimensional subspace (2D histograms)
568 /// - E Histogram the eigenvalues
569 /// - S Histogram the square of the residues
570 /// (see `TPrincipal::SumOfSquareResiduals`)
571 /// The histograms will be named `<name>_<type><number>`, where `<name>`
572 /// is the first argument, `<type>` is one of X,P,D,E,S, and `<number>`
573 /// is the variable.
574 
575 void TPrincipal::MakeHistograms(const char *name, Option_t *opt)
576 {
577  Bool_t makeX = kFALSE;
578  Bool_t makeD = kFALSE;
579  Bool_t makeP = kFALSE;
580  Bool_t makeE = kFALSE;
581  Bool_t makeS = kFALSE;
582 
583  Int_t len = strlen(opt);
584  Int_t i,j,k;
585  for (i = 0; i < len; i++) {
586  switch (opt[i]) {
587  case 'X':
588  case 'x':
589  if (fStoreData)
590  makeX = kTRUE;
591  break;
592  case 'd':
593  case 'D':
594  if (fStoreData)
595  makeD = kTRUE;
596  break;
597  case 'P':
598  case 'p':
599  if (fStoreData)
600  makeP = kTRUE;
601  break;
602  case 'E':
603  case 'e':
604  makeE = kTRUE;
605  break;
606  case 's':
607  case 'S':
608  if (fStoreData)
609  makeS = kTRUE;
610  break;
611  default:
612  Warning("MakeHistograms","Unknown option: %c",opt[i]);
613  }
614  }
615 
616  // If no option was given, then exit gracefully
617  if (!makeX && !makeD && !makeP && !makeE && !makeS)
618  return;
619 
620  // If the list of histograms doesn't exist, create it.
621  if (!fHistograms)
622  fHistograms = new TList;
623 
624  // Don't create the histograms if they are already in the TList.
625  if (makeX && fHistograms->FindObject(Form("%s_x000",name)))
626  makeX = kFALSE;
627  if (makeD && fHistograms->FindObject(Form("%s_d000",name)))
628  makeD = kFALSE;
629  if (makeP && fHistograms->FindObject(Form("%s_p000",name)))
630  makeP = kFALSE;
631  if (makeE && fHistograms->FindObject(Form("%s_e",name)))
632  makeE = kFALSE;
633  if (makeS && fHistograms->FindObject(Form("%s_s",name)))
634  makeS = kFALSE;
635 
636  TH1F **hX = 0;
637  TH2F **hD = 0;
638  TH1F **hP = 0;
639  TH1F *hE = 0;
640  TH1F *hS = 0;
641 
642  // Initialize the arrays of histograms needed
643  if (makeX)
644  hX = new TH1F * [fNumberOfVariables];
645 
646  if (makeD)
647  hD = new TH2F * [fNumberOfVariables];
648 
649  if (makeP)
650  hP = new TH1F * [fNumberOfVariables];
651 
652  if (makeE){
653  hE = new TH1F(Form("%s_e",name), "Eigenvalues of Covariance matrix",
655  hE->SetXTitle("Eigenvalue");
656  fHistograms->Add(hE);
657  }
658 
659  if (makeS) {
660  hS = new TH1F(Form("%s_s",name),"E_{N}",
662  hS->SetXTitle("N");
663  hS->SetYTitle("#sum_{i=1}^{M} (x_{i} - x'_{N,i})^{2}");
664  fHistograms->Add(hS);
665  }
666 
667  // Initialize sub elements of the histogram arrays
668  for (i = 0; i < fNumberOfVariables; i++) {
669  if (makeX) {
670  // We allow 4 sigma spread in the original data in our
671  // histogram.
672  Double_t xlowb = fMeanValues(i) - 4 * fSigmas(i);
673  Double_t xhighb = fMeanValues(i) + 4 * fSigmas(i);
675  hX[i] = new TH1F(Form("%s_x%03d", name, i),
676  Form("Pattern space, variable %d", i),
677  xbins,xlowb,xhighb);
678  hX[i]->SetXTitle(Form("x_{%d}",i));
679  fHistograms->Add(hX[i]);
680  }
681 
682  if(makeD) {
683  // The upper limit below is arbitrary!!!
684  Double_t dlowb = 0;
685  Double_t dhighb = 20;
686  Int_t dbins = fNumberOfDataPoints/100;
687  hD[i] = new TH2F(Form("%s_d%03d", name, i),
688  Form("Distance from pattern to "
689  "feature space, variable %d", i),
690  dbins,dlowb,dhighb,
692  1,
694  hD[i]->SetXTitle(Form("|x_{%d} - x'_{%d,N}|/#sigma_{%d}",i,i,i));
695  hD[i]->SetYTitle("N");
696  fHistograms->Add(hD[i]);
697  }
698 
699  if(makeP) {
700  // For some reason, the trace of the none-scaled matrix
701  // (see TPrincipal::MakeNormalised) should enter here. Taken
702  // from LINTRA code.
704  Double_t plowb = -10 * TMath::Sqrt(et);
705  Double_t phighb = -plowb;
706  Int_t pbins = 100;
707  hP[i] = new TH1F(Form("%s_p%03d", name, i),
708  Form("Feature space, variable %d", i),
709  pbins,plowb,phighb);
710  hP[i]->SetXTitle(Form("p_{%d}",i));
711  fHistograms->Add(hP[i]);
712  }
713 
714  if (makeE)
715  // The Eigenvector histogram is easy
716  hE->Fill(i,fEigenValues(i));
717 
718  }
719  if (!makeX && !makeP && !makeD && !makeS) {
720  if (hX)
721  delete[] hX;
722  if (hD)
723  delete[] hD;
724  if (hP)
725  delete[] hP;
726  return;
727  }
728 
729  Double_t *x = 0;
732  for (i = 0; i < fNumberOfDataPoints; i++) {
733 
734  // Zero arrays
735  for (j = 0; j < fNumberOfVariables; j++)
736  p[j] = d[j] = 0;
737 
738  // update the original data histogram
739  x = (Double_t*)(GetRow(i));
740  R__ASSERT(x);
741 
742  if (makeP||makeD||makeS)
743  // calculate the corresponding principal component
744  X2P(x,p);
745 
746  if (makeD || makeS) {
747  // Calculate the difference between the original data, and the
748  // same project onto principal components, and then to a lower
749  // dimensional sub-space
750  for (j = fNumberOfVariables; j > 0; j--) {
751  P2X(p,d,j);
752 
753  for (k = 0; k < fNumberOfVariables; k++) {
754  // We use the absolute value of the difference!
755  d[k] = x[k] - d[k];
756 
757  if (makeS)
758  hS->Fill(j,d[k]*d[k]);
759 
760  if (makeD) {
761  d[k] = TMath::Abs(d[k]) / (fIsNormalised ? fSigmas(k) : 1);
762  (hD[k])->Fill(d[k],j);
763  }
764  }
765  }
766  }
767 
768  if (makeX||makeP) {
769  // If we are asked to make any of these histograms, we have to
770  // go here
771  for (j = 0; j < fNumberOfVariables; j++) {
772  if (makeX)
773  (hX[j])->Fill(x[j]);
774 
775  if (makeP)
776  (hP[j])->Fill(p[j]);
777  }
778  }
779  }
780  // Clean up
781  if (hX)
782  delete [] hX;
783  if (hD)
784  delete [] hD;
785  if (hP)
786  delete [] hP;
787  if (d)
788  delete [] d;
789  if (p)
790  delete [] p;
791 
792  // Normalize the residues
793  if (makeS)
795 }
796 
797 ////////////////////////////////////////////////////////////////////////////////
798 /// Normalize the covariance matrix
799 
801 {
802  Int_t i,j;
803  for (i = 0; i < fNumberOfVariables; i++) {
805  if (fIsNormalised)
806  for (j = 0; j <= i; j++)
807  fCovarianceMatrix(i,j) /= (fSigmas(i) * fSigmas(j));
808 
809  fTrace += fCovarianceMatrix(i,i);
810  }
811 
812  // Fill remaining parts of matrix, and scale.
813  for (i = 0; i < fNumberOfVariables; i++)
814  for (j = 0; j <= i; j++) {
815  fCovarianceMatrix(i,j) /= fTrace;
817  }
818 
819 }
820 
821 ////////////////////////////////////////////////////////////////////////////////
822 /// Generate the file <classname>PCA.cxx which contains the
823 /// implementation of two methods:
824 /// ~~~ {.cpp}
825 /// void <classname>::X2P(Double_t *x, Double *p)
826 /// void <classname>::P2X(Double_t *p, Double *x, Int_t nTest)
827 /// ~~~
828 /// which does the same as TPrincipal::X2P and TPrincipal::P2X
829 /// respectively. Please refer to these methods.
830 ///
831 /// Further, the public static members:
832 /// ~~~ {.cpp}
833 /// Int_t <classname>::fgNVariables
834 /// Double_t <classname>::fgEigenValues[]
835 /// Double_t <classname>::fgEigenVectors[]
836 /// Double_t <classname>::fgMeanValues[]
837 /// Double_t <classname>::fgSigmaValues[]
838 /// ~~~
839 /// are initialized, and assumed to exist. The class declaration is
840 /// assumed to be in <classname>.h and assumed to be provided by the
841 /// user.
842 ///
843 /// See TPrincipal::MakeRealCode for a list of options
844 ///
845 /// The minimal class definition is:
846 /// ~~~ {.cpp}
847 /// class <classname> {
848 /// public:
849 /// static Int_t fgNVariables;
850 /// static Double_t fgEigenVectors[];
851 /// static Double_t fgEigenValues[];
852 /// static Double_t fgMeanValues[];
853 /// static Double_t fgSigmaValues[];
854 ///
855 /// void X2P(Double_t *x, Double_t *p);
856 /// void P2X(Double_t *p, Double_t *x, Int_t nTest);
857 /// };
858 /// ~~~
859 /// Whether the methods <classname>::X2P and <classname>::P2X should
860 /// be static or not, is up to the user.
861 
862 void TPrincipal::MakeMethods(const char *classname, Option_t *opt)
863 {
864 
865  MakeRealCode(Form("%sPCA.cxx", classname), classname, opt);
866 }
867 
868 
869 ////////////////////////////////////////////////////////////////////////////////
870 /// Perform the principal components analysis.
871 /// This is done in several stages in the TMatrix::EigenVectors method:
872 /// - Transform the covariance matrix into a tridiagonal matrix.
873 /// - Find the eigenvalues and vectors of the tridiagonal matrix.
874 
876 {
877  // Normalize covariance matrix
878  MakeNormalised();
879 
881  TMatrixDSymEigen eigen(sym);
882  fEigenVectors = eigen.GetEigenVectors();
883  fEigenValues = eigen.GetEigenValues();
884  //make sure that eigenvalues are positive
885  for (Int_t i = 0; i < fNumberOfVariables; i++) {
886  if (fEigenValues[i] < 0) fEigenValues[i] = -fEigenValues[i];
887  }
888 }
889 
890 ////////////////////////////////////////////////////////////////////////////////
891 /// This is the method that actually generates the code for the
892 /// transformations to and from feature space and pattern space
893 /// It's called by TPrincipal::MakeCode and TPrincipal::MakeMethods.
894 ///
895 /// The options are: NONE so far
896 
897 void TPrincipal::MakeRealCode(const char *filename, const char *classname,
898  Option_t *)
899 {
900  Bool_t isMethod = (classname[0] == '\0' ? kFALSE : kTRUE);
901  const char *prefix = (isMethod ? Form("%s::", classname) : "");
902  const char *cv_qual = (isMethod ? "" : "static ");
903 
904  std::ofstream outFile(filename,std::ios::out|std::ios::trunc);
905  if (!outFile) {
906  Error("MakeRealCode","couldn't open output file '%s'",filename);
907  return;
908  }
909 
910  std::cout << "Writing on file \"" << filename << "\" ... " << std::flush;
911  //
912  // Write header of file
913  //
914  // Emacs mode line ;-)
915  outFile << "// -*- mode: c++ -*-" << std::endl;
916  // Info about creator
917  outFile << "// " << std::endl
918  << "// File " << filename
919  << " generated by TPrincipal::MakeCode" << std::endl;
920  // Time stamp
921  TDatime date;
922  outFile << "// on " << date.AsString() << std::endl;
923  // ROOT version info
924  outFile << "// ROOT version " << gROOT->GetVersion()
925  << std::endl << "//" << std::endl;
926  // General information on the code
927  outFile << "// This file contains the functions " << std::endl
928  << "//" << std::endl
929  << "// void " << prefix
930  << "X2P(Double_t *x, Double_t *p); " << std::endl
931  << "// void " << prefix
932  << "P2X(Double_t *p, Double_t *x, Int_t nTest);"
933  << std::endl << "//" << std::endl
934  << "// The first for transforming original data x in " << std::endl
935  << "// pattern space, to principal components p in " << std::endl
936  << "// feature space. The second function is for the" << std::endl
937  << "// inverse transformation, but using only nTest" << std::endl
938  << "// of the principal components in the expansion" << std::endl
939  << "// " << std::endl
940  << "// See TPrincipal class documentation for more "
941  << "information " << std::endl << "// " << std::endl;
942  // Header files
943  outFile << "#ifndef __CINT__" << std::endl;
944  if (isMethod)
945  // If these are methods, we need the class header
946  outFile << "#include \"" << classname << ".h\"" << std::endl;
947  else
948  // otherwise, we need the typedefs of Int_t and Double_t
949  outFile << "#include <Rtypes.h> // needed for Double_t etc" << std::endl;
950  // Finish the preprocessor block
951  outFile << "#endif" << std::endl << std::endl;
952 
953  //
954  // Now for the data
955  //
956  // We make the Eigenvector matrix, Eigenvalue vector, Sigma vector,
957  // and Mean value vector static, since all are needed in both
958  // functions. Also ,the number of variables are stored in a static
959  // variable.
960  outFile << "//" << std::endl
961  << "// Static data variables" << std::endl
962  << "//" << std::endl;
963  outFile << cv_qual << "Int_t " << prefix << "gNVariables = "
964  << fNumberOfVariables << ";" << std::endl;
965 
966  // Assign the values to the Eigenvector matrix. The elements are
967  // stored row-wise, that is element
968  // M[i][j] = e[i * nVariables + j]
969  // where i and j are zero-based.
970  outFile << std::endl << "// Assignment of eigenvector matrix." << std::endl
971  << "// Elements are stored row-wise, that is" << std::endl
972  << "// M[i][j] = e[i * nVariables + j] " << std::endl
973  << "// where i and j are zero-based" << std::endl;
974  outFile << cv_qual << "Double_t " << prefix
975  << "gEigenVectors[] = {" << std::flush;
976  Int_t i,j;
977  for (i = 0; i < fNumberOfVariables; i++) {
978  for (j = 0; j < fNumberOfVariables; j++) {
979  Int_t index = i * fNumberOfVariables + j;
980  outFile << (index != 0 ? "," : "" ) << std::endl
981  << " " << fEigenVectors(i,j) << std::flush;
982  }
983  }
984  outFile << "};" << std::endl << std::endl;
985 
986  // Assignment to eigenvalue vector. Zero-based.
987  outFile << "// Assignment to eigen value vector. Zero-based." << std::endl;
988  outFile << cv_qual << "Double_t " << prefix
989  << "gEigenValues[] = {" << std::flush;
990  for (i = 0; i < fNumberOfVariables; i++)
991  outFile << (i != 0 ? "," : "") << std::endl
992  << " " << fEigenValues(i) << std::flush;
993  outFile << std::endl << "};" << std::endl << std::endl;
994 
995  // Assignment to mean Values vector. Zero-based.
996  outFile << "// Assignment to mean value vector. Zero-based." << std::endl;
997  outFile << cv_qual << "Double_t " << prefix
998  << "gMeanValues[] = {" << std::flush;
999  for (i = 0; i < fNumberOfVariables; i++)
1000  outFile << (i != 0 ? "," : "") << std::endl
1001  << " " << fMeanValues(i) << std::flush;
1002  outFile << std::endl << "};" << std::endl << std::endl;
1003 
1004  // Assignment to mean Values vector. Zero-based.
1005  outFile << "// Assignment to sigma value vector. Zero-based." << std::endl;
1006  outFile << cv_qual << "Double_t " << prefix
1007  << "gSigmaValues[] = {" << std::flush;
1008  for (i = 0; i < fNumberOfVariables; i++)
1009  outFile << (i != 0 ? "," : "") << std::endl
1010  << " " << (fIsNormalised ? fSigmas(i) : 1) << std::flush;
1011  // << " " << fSigmas(i) << std::flush;
1012  outFile << std::endl << "};" << std::endl << std::endl;
1013 
1014  //
1015  // Finally we reach the functions themselves
1016  //
1017  // First: void x2p(Double_t *x, Double_t *p);
1018  //
1019  outFile << "// " << std::endl
1020  << "// The "
1021  << (isMethod ? "method " : "function ")
1022  << " void " << prefix
1023  << "X2P(Double_t *x, Double_t *p)"
1024  << std::endl << "// " << std::endl;
1025  outFile << "void " << prefix
1026  << "X2P(Double_t *x, Double_t *p) {" << std::endl
1027  << " for (Int_t i = 0; i < gNVariables; i++) {" << std::endl
1028  << " p[i] = 0;" << std::endl
1029  << " for (Int_t j = 0; j < gNVariables; j++)" << std::endl
1030  << " p[i] += (x[j] - gMeanValues[j]) " << std::endl
1031  << " * gEigenVectors[j * gNVariables + i] "
1032  << "/ gSigmaValues[j];" << std::endl << std::endl << " }"
1033  << std::endl << "}" << std::endl << std::endl;
1034  //
1035  // Now: void p2x(Double_t *p, Double_t *x, Int_t nTest);
1036  //
1037  outFile << "// " << std::endl << "// The "
1038  << (isMethod ? "method " : "function ")
1039  << " void " << prefix
1040  << "P2X(Double_t *p, Double_t *x, Int_t nTest)"
1041  << std::endl << "// " << std::endl;
1042  outFile << "void " << prefix
1043  << "P2X(Double_t *p, Double_t *x, Int_t nTest) {" << std::endl
1044  << " for (Int_t i = 0; i < gNVariables; i++) {" << std::endl
1045  << " x[i] = gMeanValues[i];" << std::endl
1046  << " for (Int_t j = 0; j < nTest; j++)" << std::endl
1047  << " x[i] += p[j] * gSigmaValues[i] " << std::endl
1048  << " * gEigenVectors[i * gNVariables + j];" << std::endl
1049  << " }" << std::endl << "}" << std::endl << std::endl;
1050 
1051  // EOF
1052  outFile << "// EOF for " << filename << std::endl;
1053 
1054  // Close the file
1055  outFile.close();
1056 
1057  std::cout << "done" << std::endl;
1058 }
1059 
1060 ////////////////////////////////////////////////////////////////////////////////
1061 /// Calculate x as a function of nTest of the most significant
1062 /// principal components p, and return it in x.
1063 /// It's the users responsibility to make sure that both x and p are
1064 /// of the right size (i.e., memory must be allocated for x).
1065 
1066 void TPrincipal::P2X(const Double_t *p, Double_t *x, Int_t nTest)
1067 {
1068  for (Int_t i = 0; i < fNumberOfVariables; i++){
1069  x[i] = fMeanValues(i);
1070  for (Int_t j = 0; j < nTest; j++)
1071  x[i] += p[j] * (fIsNormalised ? fSigmas(i) : 1)
1072  * fEigenVectors(i,j);
1073  }
1074 
1075 }
1076 
1077 ////////////////////////////////////////////////////////////////////////////////
1078 /// Print the statistics
1079 /// Options are
1080 /// - M Print mean values of original data
1081 /// - S Print sigma values of original data
1082 /// - E Print eigenvalues of covariance matrix
1083 /// - V Print eigenvectors of covariance matrix
1084 /// Default is MSE
1085 
1086 void TPrincipal::Print(Option_t *opt) const
1087 {
1088  Bool_t printV = kFALSE;
1089  Bool_t printM = kFALSE;
1090  Bool_t printS = kFALSE;
1091  Bool_t printE = kFALSE;
1092 
1093  Int_t len = strlen(opt);
1094  for (Int_t i = 0; i < len; i++) {
1095  switch (opt[i]) {
1096  case 'V':
1097  case 'v':
1098  printV = kTRUE;
1099  break;
1100  case 'M':
1101  case 'm':
1102  printM = kTRUE;
1103  break;
1104  case 'S':
1105  case 's':
1106  printS = kTRUE;
1107  break;
1108  case 'E':
1109  case 'e':
1110  printE = kTRUE;
1111  break;
1112  default:
1113  Warning("Print", "Unknown option '%c'",opt[i]);
1114  break;
1115  }
1116  }
1117 
1118  if (printM||printS||printE) {
1119  std::cout << " Variable # " << std::flush;
1120  if (printM)
1121  std::cout << "| Mean Value " << std::flush;
1122  if (printS)
1123  std::cout << "| Sigma " << std::flush;
1124  if (printE)
1125  std::cout << "| Eigenvalue" << std::flush;
1126  std::cout << std::endl;
1127 
1128  std::cout << "-------------" << std::flush;
1129  if (printM)
1130  std::cout << "+------------" << std::flush;
1131  if (printS)
1132  std::cout << "+------------" << std::flush;
1133  if (printE)
1134  std::cout << "+------------" << std::flush;
1135  std::cout << std::endl;
1136 
1137  for (Int_t i = 0; i < fNumberOfVariables; i++) {
1138  std::cout << std::setw(12) << i << " " << std::flush;
1139  if (printM)
1140  std::cout << "| " << std::setw(10) << std::setprecision(4)
1141  << fMeanValues(i) << " " << std::flush;
1142  if (printS)
1143  std::cout << "| " << std::setw(10) << std::setprecision(4)
1144  << fSigmas(i) << " " << std::flush;
1145  if (printE)
1146  std::cout << "| " << std::setw(10) << std::setprecision(4)
1147  << fEigenValues(i) << " " << std::flush;
1148  std::cout << std::endl;
1149  }
1150  std::cout << std::endl;
1151  }
1152 
1153  if(printV) {
1154  for (Int_t i = 0; i < fNumberOfVariables; i++) {
1155  std::cout << "Eigenvector # " << i << std::flush;
1158  v.Print();
1159  }
1160  }
1161 }
1162 
1163 ////////////////////////////////////////////////////////////////////////////////
1164 /// Calculates the sum of the square residuals, that is
1165 ///
1166 /// \f[
1167 /// E_N = \sum_{i=0}^{P-1} \left(x_i - x^\prime_i\right)^2
1168 /// \f]
1169 ///
1170 /// where \f$x^\prime_i = \sum_{j=i}^N p_i e_{n_j}\f$
1171 /// is the \f$i^{\mbox{th}}\f$ component of the principal vector, corresponding to
1172 /// \f$x_i\f$, the original data; I.e., the square distance to the space
1173 /// spanned by \f$N\f$ eigenvectors.
1174 
1176 {
1177 
1178  if (!x)
1179  return;
1180 
1181  Double_t p[100];
1182  Double_t xp[100];
1183 
1184  X2P(x,p);
1185  for (Int_t i = fNumberOfVariables-1; i >= 0; i--) {
1186  P2X(p,xp,i);
1187  for (Int_t j = 0; j < fNumberOfVariables; j++) {
1188  s[i] += (x[j] - xp[j])*(x[j] - xp[j]);
1189  }
1190  }
1191 }
1192 
1193 ////////////////////////////////////////////////////////////////////////////////
1194 /// Test the PCA, bye calculating the sum square of residuals
1195 /// (see method SumOfSquareResiduals), and display the histogram
1196 
1198 {
1199  MakeHistograms("pca","S");
1200 
1201  if (!fStoreData)
1202  return;
1203 
1204  TH1 *pca_s = 0;
1205  if (fHistograms) pca_s = (TH1*)fHistograms->FindObject("pca_s");
1206  if (!pca_s) {
1207  Warning("Test", "Couldn't get histogram of square residuals");
1208  return;
1209  }
1210 
1211  pca_s->Draw();
1212 }
1213 
1214 ////////////////////////////////////////////////////////////////////////////////
1215 /// Calculate the principal components from the original data vector
1216 /// x, and return it in p.
1217 ///
1218 /// It's the users responsibility to make sure that both x and p are
1219 /// of the right size (i.e., memory must be allocated for p).
1220 
1222 {
1223  for (Int_t i = 0; i < fNumberOfVariables; i++){
1224  p[i] = 0;
1225  for (Int_t j = 0; j < fNumberOfVariables; j++)
1226  p[i] += (x[j] - fMeanValues(j)) * fEigenVectors(j,i) /
1227  (fIsNormalised ? fSigmas(j) : 1);
1228  }
1229 
1230 }
TPrincipal::Clear
virtual void Clear(Option_t *option="")
Clear the data in Object.
Definition: TPrincipal.cxx:486
TPrincipal::fStoreData
Bool_t fStoreData
Definition: TPrincipal.h:43
TPrincipal::MakeHistograms
virtual void MakeHistograms(const char *name="pca", Option_t *option="epsdx")
Make histograms of the result of the analysis.
Definition: TPrincipal.cxx:575
TPrincipal::X2P
virtual void X2P(const Double_t *x, Double_t *p)
Calculate the principal components from the original data vector x, and return it in p.
Definition: TPrincipal.cxx:1221
TPrincipal.h
TPrincipal::GetRow
const Double_t * GetRow(Int_t row)
Return a row of the user supplied data.
Definition: TPrincipal.cxx:513
TBrowser
Using a TBrowser one can browse all ROOT objects.
Definition: TBrowser.h:37
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:91
TVectorD.h
sym
#define sym(otri1, otri2)
Definition: triangle.c:932
TNamed::SetName
virtual void SetName(const char *name)
Set the name of the TNamed.
Definition: TNamed.cxx:140
Option_t
const char Option_t
Definition: RtypesCore.h:66
TPrincipal::AddRow
virtual void AddRow(const Double_t *x)
Add a data point and update the covariance matrix.
Definition: TPrincipal.cxx:410
TVectorT::ResizeTo
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition: TVectorT.cxx:294
TH2F
2-D histogram with a float per channel (see TH1 documentation)}
Definition: TH2.h:251
TList::FindObject
virtual TObject * FindObject(const char *name) const
Find an object in this list using its name.
Definition: TList.cxx:578
TList::Delete
virtual void Delete(Option_t *option="")
Remove all objects from the list AND delete all heap based objects.
Definition: TList.cxx:470
TString::Data
const char * Data() const
Definition: TString.h:369
TNamed::operator=
TNamed & operator=(const TNamed &rhs)
TNamed assignment operator.
Definition: TNamed.cxx:51
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:364
Form
char * Form(const char *fmt,...)
TPrincipal::MakePrincipals
virtual void MakePrincipals()
Perform the principal components analysis.
Definition: TPrincipal.cxx:875
TDatime.h
TObject::Error
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:890
TMath::Sqrt
Double_t Sqrt(Double_t x)
Definition: TMath.h:691
TGeant4Unit::s
static constexpr double s
Definition: TGeant4SystemOfUnits.h:162
TMatrixDColumn_const
TMatrixTColumn_const< Double_t > TMatrixDColumn_const
Definition: TMatrixDUtilsfwd.h:55
TDatime::AsString
const char * AsString() const
Return the date & time as a string (ctime() format).
Definition: TDatime.cxx:102
TBrowser.h
TString::EndsWith
Bool_t EndsWith(const char *pat, ECaseCompare cmp=kExact) const
Return true if string ends with the specified string.
Definition: TString.cxx:2177
TPrincipal::fNumberOfDataPoints
Int_t fNumberOfDataPoints
Definition: TPrincipal.h:24
TPrincipal::operator=
TPrincipal & operator=(const TPrincipal &)
Assignment operator.
Definition: TPrincipal.cxx:333
x
Double_t x[n]
Definition: legend1.C:17
TPrincipal::fOffDiagonal
TVectorD fOffDiagonal
Definition: TPrincipal.h:34
TList.h
TMatrixT::GetMatrixArray
virtual const Element * GetMatrixArray() const
Definition: TMatrixT.h:222
TMatrixTSym< Double_t >
TMath::Abs
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
TPrincipal::fCovarianceMatrix
TMatrixD fCovarianceMatrix
Definition: TPrincipal.h:29
TPrincipal::TPrincipal
TPrincipal()
Empty constructor. Do not use.
Definition: TPrincipal.cxx:229
TH1::SetXTitle
virtual void SetXTitle(const char *title)
Definition: TH1.h:413
TString
Basic string class.
Definition: TString.h:136
TH1::SetYTitle
virtual void SetYTitle(const char *title)
Definition: TH1.h:414
v
@ v
Definition: rootcling_impl.cxx:3635
b
#define b(i)
Definition: RSha256.hxx:100
bool
TVectorT::GetNrows
Int_t GetNrows() const
Definition: TVectorT.h:75
TROOT.h
TVectorT::IsValid
Bool_t IsValid() const
Definition: TVectorT.h:83
TPrincipal::MakeCode
virtual void MakeCode(const char *filename="pca", Option_t *option="")
Generates the file <filename>, with .C appended if it does argument doesn't end in ....
Definition: TPrincipal.cxx:550
TMatrixDSymEigen.h
TPrincipal::P2X
virtual void P2X(const Double_t *p, Double_t *x, Int_t nTest)
Calculate x as a function of nTest of the most significant principal components p,...
Definition: TPrincipal.cxx:1066
h
#define h(i)
Definition: RSha256.hxx:106
TMatrixDSymEigen::GetEigenVectors
const TMatrixD & GetEigenVectors() const
Definition: TMatrixDSymEigen.h:53
TPrincipal::fHistograms
TList * fHistograms
Definition: TPrincipal.h:40
TNamed
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:92
TPrincipal
Principal Components Analysis (PCA)
Definition: TPrincipal.h:21
TPrincipal::Browse
virtual void Browse(TBrowser *b)
Browse the TPrincipal object in the TBrowser.
Definition: TPrincipal.cxx:463
TH1::Fill
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3346
TPrincipal::SumOfSquareResiduals
virtual void SumOfSquareResiduals(const Double_t *x, Double_t *s)
Calculates the sum of the square residuals, that is.
Definition: TPrincipal.cxx:1175
TPrincipal::fSigmas
TVectorD fSigmas
Definition: TPrincipal.h:28
TPrincipal::fEigenVectors
TMatrixD fEigenVectors
Definition: TPrincipal.h:31
TH2.h
TObject::Warning
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:876
TVectorT::GetMatrixArray
Element * GetMatrixArray()
Definition: TVectorT.h:78
TPrincipal::fIsNormalised
Bool_t fIsNormalised
Definition: TPrincipal.h:42
TPrincipal::Print
virtual void Print(Option_t *opt="MSE") const
Print the statistics Options are.
Definition: TPrincipal.cxx:1086
TVectorT< Double_t >
TPrincipal::~TPrincipal
virtual ~TPrincipal()
Destructor.
Definition: TPrincipal.cxx:357
TPrincipal::MakeRealCode
void MakeRealCode(const char *filename, const char *prefix, Option_t *option="")
This is the method that actually generates the code for the transformations to and from feature space...
Definition: TPrincipal.cxx:897
Double_t
double Double_t
Definition: RtypesCore.h:59
TPrincipal::fMeanValues
TVectorD fMeanValues
Definition: TPrincipal.h:27
R__ASSERT
#define R__ASSERT(e)
Definition: TError.h:120
TPrincipal::fNumberOfVariables
Int_t fNumberOfVariables
Definition: TPrincipal.h:25
TPrincipal::MakeMethods
virtual void MakeMethods(const char *classname="PCA", Option_t *option="")
Generate the file <classname>PCA.cxx which contains the implementation of two methods:
Definition: TPrincipal.cxx:862
t1
auto * t1
Definition: textangle.C:20
TH1F
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:575
TList::Add
virtual void Add(TObject *obj)
Definition: TList.h:87
TPrincipal::MakeNormalised
void MakeNormalised()
Normalize the covariance matrix.
Definition: TPrincipal.cxx:800
TVectorT::Zero
TVectorT< Element > & Zero()
Set vector elements to zero.
Definition: TVectorT.cxx:453
TH1
TH1 is the base class of all histogramm classes in ROOT.
Definition: TH1.h:58
TPrincipal::Test
void Test(Option_t *option="")
Test the PCA, bye calculating the sum square of residuals (see method SumOfSquareResiduals),...
Definition: TPrincipal.cxx:1197
TMatrixDSymEigen::GetEigenValues
const TVectorD & GetEigenValues() const
Definition: TMatrixDSymEigen.h:54
name
char name[80]
Definition: TGX11.cxx:110
d
#define d(i)
Definition: RSha256.hxx:102
TIter
Definition: TCollection.h:233
TDatime
This class stores the date and time with a precision of one second in an unsigned 32 bit word (950130...
Definition: TDatime.h:37
TMatrixD.h
Riostream.h
TPrincipal::fEigenValues
TVectorD fEigenValues
Definition: TPrincipal.h:32
TPrincipal::fTrace
Double_t fTrace
Definition: TPrincipal.h:38
TH1::Scale
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this histogram by a constant c1.
Definition: TH1.cxx:6560
TMatrixDSymEigen
TMatrixDSymEigen.
Definition: TMatrixDSymEigen.h:28
TList
A doubly linked list.
Definition: TList.h:44
TPrincipal::fUserData
TVectorD fUserData
Definition: TPrincipal.h:36
TMath.h
fw3dlego::xbins
const double xbins[xbins_n]
Definition: collection_proxies.C:48
gROOT
#define gROOT
Definition: TROOT.h:406
int
TH1::Draw
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:3069