#include <iostream>
#include <iomanip>

#include "TH2D.h"
#include "TMatrixD.h"
#include "TPrincipal.h"


/*==PCA Test==========================================================
The macro tests the Principal Component Analysis (PCA) class in ROOT.
A known set of ten data points in two dimensions is given to the
TPrincipal class.  The results calcaulted by hand are printed below.

Covariance Matrix of the data
-------------------------------
|  0.616555556 |  0.615444444 |
|  0.615444444 |  0.716555556 |
-------------------------------

eigenvalues for the covariance matrix
----------------
|  0.049083398 |
|  1.28402771  |
----------------

eigenvectors for the covariance matrix
-------------------------------
|  -0.73517866 |  -0.67787340 |
|  0.677873399 |  -0.73517866 |
-------------------------------
====================================================================*/


Int_t PCATest() {

  //--Data--
  const Int_t numVariables = 2;
  const Int_t numDataPoints = 10;
  const Double_t inputdata[numDataPoints][numVariables] = 
    { 2.5, 2.4,
      0.5, 0.7,
      2.2, 2.9,
      1.9, 2.2,
      3.1, 3.0,
      2.3, 2.7,
      2.0, 1.6,
      1.0, 1.1,
      1.5, 1.6,
      1.1, 0.9 };
  
  //--Constants--
  const Int_t numBinsX=100;
  const Double_t minX=-1.0, maxX=4.0;

  //--Variables--
  TMatrixD finalCov = TMatrixD(2,2);
  
  TH2D* hXY = new TH2D("hXY","Orignal Data", numBinsX, minX, maxX, numBinsX, minX, maxX);

  TPrincipal *pcatest = new TPrincipal(numVariables);
  Double_t* inputarray = new Double_t[numVariables];

  //--Main Algorithm--------------------------------------------------

  //Add data to pcatest
  for (Int_t ii=0; ii<numDataPoints; ++ii) {
    inputarray[0] = inputdata[ii][0];
    inputarray[1] = inputdata[ii][1];

    for (int jj=0; jj<numVariables; ++jj) {
      cout << setw(19) << " " << setw(8) << inputarray[jj];
    }//ii
    cout << endl;

    hXY->Fill(inputarray[0], inputarray[1]);
    pcatest->AddRow(inputarray);
  }//ii

  //Run PCA
  pcatest->MakePrincipals();
  
  //Print the output
  pcatest->MakeHistograms("pcatest","XPES");
  pcatest->Print("MSEV");
  cout.flush();

  //Print out the final covariance matrix
  finalCov = *(pcatest->GetCovarianceMatrix());
  finalCov.Print();
  cout.flush();

  //Draw the orignal data histogram
  hXY->SetMarkerStyle(2);
  hXY->Draw();

  //start a TBrowser
  TBrowser* p = new TBrowser("principalBrowser", pcatest);

  return 0;
}//Int_t PCATest()