#include "Analyzer.h"



using namespace std;

Analyzer::Analyzer()
{
  hfile      = new TFile("Analyzer.root","RECREATE","XXX");
  hrating_per_day = new TH1F("rating_per_day","rating_per_day",3000 ,0.5 ,3000.5 );
  hrating_per_month = new TH1F("rating_per_month","rating_per_month",100 ,0.5 ,100.5 );
  hrating_all = new TH1F("rating_all","rating_all",5 ,.5 ,5.5 );
  havg_rating = new TH1F("avg_rating","avg_rating",100, .5 ,100.5 );
  havg_rating_d = new TH1F("avg_rating_d","avg_rating_d",100, .5 ,100.5 );
  hratings_per_user = new TH1F("ratings_per_user","ratings_per_user",150 ,0. ,1500 );
  htest_ratings_per_user = new TH1F("test_ratings_per_user","test_ratings_per_user",20 ,0.5 ,20.5 );
  hratings_vs_test = new TH2F("ratings_vs_test","ratings_vs_test", 150, 0.5,1500.5,20, .5, 20.5 );
  hfirst_rating_each_user = new TH1F("first_rating_each_user","first_rating_each_user",100 ,0.5 ,100.5 );
  hlast_rating_each_user = new TH1F("last_rating_each_user","last_rating_each_user",100 ,0.5 ,100.5 );
  htime_1st_to_last_rating= new TH1F("time_1st_to_last_rating","time_1st_to_last_rating",101 ,-.5 ,100.5 );
  hratings_per_session_per_user = new TH1F("ratings_per_session_per_user","ratings_per_session_per_user",100 , .5,100.5 );
  hrating_sessions_per_user= new TH1F("rating_sessions_per_user","rating_sessions_per_user", 100,.5 ,100.5 );


}
//   h = new TH1F("","", , , );


Analyzer::~Analyzer()
{
// Save all objects in this file
hfile->Write();
}

void
Analyzer::a1(char* filename)
{
  int count1=0;
  int count10k=0;
  int current_user =-1;
  int r_p_user=0;
  int r_p_0_user=0;
  int prev_lmonth=0;
  int first_lmonth=0;
  int current_lday=-1; 
  int rating_sessions=0;;
  int ratings_this_session=0;;


 datafile.open(filename);

  while ( !datafile.eof()         )
    {

      datafile.getline( &line[0], 5000000, '\n');
      
      if (datafile.eof()) break;
      
      count1++;
      count10k++;
      if (count10k>9999){
	count10k = 0;
	cout<<count1<< endl;
      }


      int movie = atoi( strtok( line, "\t" ) );     
      int user = atoi( strtok( NULL, "\t" ) );     
      int rating = atoi( strtok( NULL, "\t" ) );     
      int year = atoi( strtok( NULL, "-" ) );     
      int month = atoi( strtok( NULL, "-" ) );     
      int day = atoi( strtok( NULL, "\n" ) );     
 

      int lday = (year-1999)*365 +month*31+day;
      int lmonth = (year-1999)*12 +month;
      hrating_per_day->Fill(lday);
      hrating_per_month->Fill(lmonth);

      hrating_all->Fill(rating);

      if (rating>0){
	havg_rating->Fill(lmonth, rating);
	havg_rating_d->Fill(lmonth);
      }

      if (rating !=0) r_p_user++;
      if (rating ==0) r_p_0_user++;


      if (current_user != user){  // user switched
	if (current_user >0){  // not first one
	  hratings_per_user->Fill(r_p_user);
	  htest_ratings_per_user->Fill(r_p_0_user);
	  hratings_vs_test->Fill(r_p_user,r_p_0_user);
	  r_p_user=0;
	  r_p_0_user=0;


	  hfirst_rating_each_user->Fill(lmonth);
	  hlast_rating_each_user->Fill(prev_lmonth);
	  htime_1st_to_last_rating->Fill(prev_lmonth - first_lmonth);
	  hrating_sessions_per_user->Fill(rating_sessions);
	  rating_sessions=1;

	  hratings_per_session_per_user->Fill(ratings_this_session);
	  ratings_this_session=1;
	  current_lday = lday;
	}


	rating_sessions=1;
	current_lday=lday;

	first_lmonth=lmonth;
	current_user = user;
      } // user switch
      else{
	//same user
	if (current_lday != lday){  // new rating session
	  hratings_per_session_per_user->Fill(ratings_this_session);
	  ratings_this_session=1;
	  rating_sessions++;


	}
      ratings_this_session++;


      }


      prev_lmonth=lmonth;
      current_lday = lday;
    }  // while over file


}


// -------------------------------------------------------------------------------
void
Analyzer::makesmall(char* filename)
{
  int count1=0;
  int count10k=0;


  outfile.open("/tmp/user_date100.txt");
  datafile.open(filename);

  int current_user =-1;
  int user_cnt=0;
  int user_cnt10k=0;
  int output_user=0;

  while ( !datafile.eof()         )
    {

      datafile.getline( &line[0], 5000000, '\n');
      strcpy(line1, line);
      if (datafile.eof()) break;
      
      count1++;
      count10k++;
      if (count10k>9999){
	count10k = 0;
	cout<<count1<< endl;

      }


      int movie = atoi( strtok( line, "\t" ) );     
      int user = atoi( strtok( NULL, "\t" ) );     
      int rating = atoi( strtok( NULL, "\t" ) );     
      int year = atoi( strtok( NULL, "-" ) );     
      int month = atoi( strtok( NULL, "-" ) );     
      int day = atoi( strtok( NULL, "\n" ) );   



      if (current_user != user){
	output_user=0;
	current_user = user;
	user_cnt++;
	user_cnt10k++;
	if (user_cnt10k>99){
	  user_cnt10k=0;
	  output_user=1;
	}



      }

      if (output_user) outfile<<line1<<endl;



    }


}


// -------------------------------------------------------------------------------
// -------------------------------------------------------------------------------
void
Analyzer::m1()
{



// This macro shows several ways to invert a matrix . Each  method
// is a trade-off between accuracy of the inversion and speed.
// Which method to chose depends on "how well-behaved" the matrix is.
// This is best checked through a call to Condition(), available in each
// decomposition class. A second possibilty (less preferred) would be to
// check the determinant
//
// You can run this script with different matrix sizes


  Int_t msize=6;

  if (msize < 2 || msize > 10) {
    cout << "2 <= msize <= 10" <<endl;
    return;
  }
  //gSystem->Load("libMatrix");


  cout << "--------------------------------------------------------" <<endl;
  cout << "Inversion results for a ("<<msize<<","<<msize<<") matrix" <<endl;
  cout << "For each inversion procedure we check the maxmimum size " <<endl;
  cout << "of the off-diagonal elements of Inv(A) * A              " <<endl;
  cout << "--------------------------------------------------------" <<endl;

  TMatrixD H_square = THilbertMatrixD(msize,msize);

//  1. InvertFast(Double_t *det=0)
//   It is identical to Invert() for sizes > 6 x 6 but for smaller sizes, the
//   inversion is performed according to Cramer's rule by explicitly calculating
//   all Jacobi's sub-determinants . For instance for a 6 x 6 matrix this means:
//    # of 5 x 5 determinant : 36 
//    # of 4 x 4 determinant : 75 
//    # of 3 x 3 determinant : 80 
//    # of 2 x 2 determinant : 45    (see TMatrixD/FCramerInv.cxx)
//
//    The only "quality" control in this process is to check whether the 6 x 6
//    determinant is unequal 0 . But speed gains are significant compared to Invert() ,
//    upto an order of magnitude for sizes <= 4 x 4
//
//    The inversion is done "in place", so the original matrix will be overwritten
//    If a pointer to a Double_t is supplied the determinant is calculated
//

  cout << "1. Use .InvertFast(&det)" <<endl;
  if (msize > 6)
    cout << " for ("<<msize<<","<<msize<<") this is identical to .Invert(&det)" <<endl;

  Double_t det1;
  TMatrixD H1 = H_square;
  H1.InvertFast(&det1);

  // Get the maximum off-diagonal matrix value . One way to do this is to set the
  // diagonal to zero .

  TMatrixD U1(H1,TMatrixD::kMult,H_square);
  TMatrixDDiag diag1(U1); diag1 = 0.0;
  const Double_t U1_max_offdiag = (U1.Abs()).Max();
  cout << "  Maximum off-diagonal = " << U1_max_offdiag << endl;
  cout << "  Determinant          = " << det1 <<endl;

// 2. Invert(Double_t *det=0)
//   Again the inversion is performed in place .
//   It consists out of a sequence of calls to the decomposition classes . For instance
//   for the general dense matrix TMatrixD the LU decomposition is invoked:
//    - The matrix is decomposed using a scheme according to Crout which involves
//      "implicit partial pivoting", see for instance Num. Recip. (we have also available
//      a decomposition scheme that does not the scaling and is therefore even slightly
//      faster but less stable)
//      With each decomposition, a tolerance has to be specified . If this tolerance
//      requirement is not met, the matrix is regarded as being singular. The value
//      passed to this decomposition, is the data member fTol of the matrix . Its
//      default value is DBL_EPSILON, which is defined as the smallest nuber so that
//      1+DBL_EPSILON > 1
//    - The last step is a standard forward/backward substitution .
//
//   It is important to realize that both InvertFast() and Invert() are "one-shot" deals , speed
//   comes at a price . If something goes wrong because the matrix is (near) singular, you have
//   overwritten your original matrix and  no factorization is available anymore to get more
//   information like condition number or change the tolerance number .
//
//   All other calls in the matrix classes involving inversion like the ones with the "smart"
//   constructors (kInverted,kInvMult...) use this inversion method .
//

  cout << "2. Use .Invert(&det)" <<endl;

  Double_t det2;
  TMatrixD H2 = H_square;
  H2.Invert(&det2);

  TMatrixD U2(H2,TMatrixD::kMult,H_square);
  TMatrixDDiag diag2(U2); diag2 = 0.0;
  const Double_t U2_max_offdiag = (U2.Abs()).Max();
  cout << "  Maximum off-diagonal = " << U2_max_offdiag << endl;
  cout << "  Determinant          = " << det2 <<endl;

// 3. Inversion through LU decomposition
//   The (default) algorithms used are similar to 2. (Not identical because in 2, the whole
//   calculation is done "in-place". Here the orginal matrix is copied (so more memory
//   management => slower) and several operations can be performed without having to repeat
//   the decomposition step .
//   Inverting a matrix is nothing else than solving a set of equations where the rhs is given
//   by the unit matrix, so the steps to take are identical to those solving a linear equation :
//

  cout << "3. Use TDecompLU" <<endl;

  TMatrixD H3 = H_square;
  TDecompLU lu(H_square);

  // Any operation that requires a decomposition will trigger it . The class keeps
  // an internal state so that following operations will not perform the decomposition again
  // unless the matrix is changed through SetMatrix(..)
  // One might want to proceed more cautiously by invoking first Decompose() and check its
  // return value before proceeding....

  lu.Invert(H3);
  Double_t d1_lu; Double_t d2_lu;
  lu.Det(d1_lu,d2_lu);
  Double_t det3 = d1_lu*TMath::Power(2.,d2_lu);

  TMatrixD U3(H3,TMatrixD::kMult,H_square);
  TMatrixDDiag diag3(U3); diag3 = 0.0;
  const Double_t U3_max_offdiag = (U3.Abs()).Max();
  cout << "  Maximum off-diagonal = " << U3_max_offdiag << endl;
  cout << "  Determinant          = " << det3 <<endl;

// 4. Inversion through SVD decomposition
//   For SVD and QRH, the (n x m) matrix does only have to fulfill n >=m . In case n > m
//   a pseudo-inverse is calculated
  cout << "4. Use TDecompSVD on non-square matrix" <<endl;

  TMatrixD H_nsquare = THilbertMatrixD(msize,msize-1);

  TDecompSVD svd(H_nsquare);

  TMatrixD H4 = svd.Invert();
  Double_t d1_svd; Double_t d2_svd;
  svd.Det(d1_svd,d2_svd);
  Double_t det4 = d1_svd*TMath::Power(2.,d2_svd);

  TMatrixD U4(H4,TMatrixD::kMult,H_nsquare);
  TMatrixDDiag diag4(U4); diag4 = 0.0;
  const Double_t U4_max_offdiag = (U4.Abs()).Max();
  cout << "  Maximum off-diagonal = " << U4_max_offdiag << endl;
  cout << "  Determinant          = " << det4 <<endl;
}