// tsncom31.cc
// Nick Macks (macks@fnal.gov), Summer 2002

// CovMatrices: CovMatrixX trace calculation

#include "CovMatrices/CovMatrix2.h"
#include "CovMatrices/CovMatrix3.h"
#include "CovMatrices/CovMatrix4.h"
#include "CovMatrices/CovMatrix5.h"
#include "CovMatrices/CovMatrix6.h"

using namespace std;

using CovMatrices::CovMatrix2;
using CovMatrices::CovMatrix3;
using CovMatrices::CovMatrix4;
using CovMatrices::CovMatrix5;
using CovMatrices::CovMatrix6;

// Function declarations
void Identification();
void GetDimension(int&);
void Describe(const int&);
void BuildCovMatrix2(CovMatrix2&);
void BuildCovMatrix3(CovMatrix3&);
void BuildCovMatrix4(CovMatrix4&);
void BuildCovMatrix5(CovMatrix5&);
void BuildCovMatrix6(CovMatrix6&);
void Show2(const CovMatrix2&);
void Show3(const CovMatrix3&);
void Show4(const CovMatrix4&);
void Show5(const CovMatrix5&);
void Show6(const CovMatrix6&);
void GetCalcs(int&);
void GetRuns(int&);
void Calculate2(CovMatrix2&, const int&, const int&);
void Calculate3(CovMatrix3&, const int&, const int&);
void Calculate4(CovMatrix4&, const int&, const int&);
void Calculate5(CovMatrix5&, const int&, const int&);
void Calculate6(CovMatrix6&, const int&, const int&);


// Package ID
void Identification()
{
  cout << "Using CovMatrices\n";
}

// Ask for the desired dimension
void GetDimension(int& dimension)
{
  cout << "Enter the desired dimension(2-6): ";
  cin >> dimension;

  while ( (dimension < 2) || (dimension > 6) )
  {
    cout << "Enter an integer between 2 and 6: ";
    cin >> dimension;
  }
}

// Describe how the matrix was defined
void Describe(const int& dimension)
{
    cout << "Matrix was defined 'CovMatrix" << dimension << "'\n";
}

// Built a 2-dimension CovMatrix
void BuildCovMatrix2(CovMatrix2& matrix2)
{
  for (int i=0; i<2; ++i)
  {
    for (int j=0; j<2; ++j)
    {
      if (i == j)
        matrix2(i,j) = 1000 + i*100;
      else if (i > j)
        matrix2(i,j) = 2 + i*2 + j*(2*2);
    }
  }
}

// Built a 3-dimension CovMatrix
void BuildCovMatrix3(CovMatrix3& matrix3)
{
  for (int i=0; i<3; ++i)
  {
    for (int j=0; j<3; ++j)
    {
      if (i == j)
        matrix3(i,j) = 1000 + i*100;
      else if (i > j)
        matrix3(i,j) = 2 + i*2 + j*(3*2);
    }
  }
}

// Built a 4-dimension CovMatrix
void BuildCovMatrix4(CovMatrix4& matrix4)
{
  for (int i=0; i<4; ++i)
  {
    for (int j=0; j<4; ++j)
    {
      if (i == j)
        matrix4(i,j) = 1000 + i*100;
      else if (i > j)
        matrix4(i,j) = 2 + i*2 + j*(4*2);
    }
  }
}

// Built a 5-dimension CovMatrix
void BuildCovMatrix5(CovMatrix5& matrix5)
{
  for (int i=0; i<5; ++i)
  {
    for (int j=0; j<5; ++j)
    {
      if (i == j)
        matrix5(i,j) = 1000 + i*100;
      else if (i > j)
        matrix5(i,j) = 2 + i*2 + j*(5*2);
    }
  }
}

// Built a 6-dimension CovMatrix
void BuildCovMatrix6(CovMatrix6& matrix6)
{
  for (int i=0; i<6; ++i)
  {
    for (int j=0; j<6; ++j)
    {
      if (i == j)
        matrix6(i,j) = 1000 + i*100;
      else if (i > j)
        matrix6(i,j) = 2 + i*2 + j*(6*2);
    }
  }
}

// Show the 2-dimensional CovMatrix
void Show2(const CovMatrix2& matrix2)
{
  cout << "Using matrix =\n" << matrix2 << "\n";
}

// Show the 3-dimensional CovMatrix
void Show3(const CovMatrix3& matrix3)
{
  cout << "Using matrix =\n" << matrix3 << "\n";
}

// Show the 4-dimensional CovMatrix
void Show4(const CovMatrix4& matrix4)
{
  cout << "Using matrix =\n" << matrix4 << "\n";
}

// Show the 5-dimensional CovMatrix
void Show5(const CovMatrix5& matrix5)
{
  cout << "Using matrix =\n" << matrix5 << "\n";
}

// Show the 6-dimensional CovMatrix
void Show6(const CovMatrix6& matrix6)
{
  cout << "Using matrix =\n" << matrix6 << "\n";
}

// Get the number of times to calculate the trace
void GetCalcs(int& N)
{
  cout << "How many trace calculations? ";
  cin >> N;
}

// Get the number of runs
void GetRuns(int& run)
{
  cout << "How many runs? ";
  cin >> run;
  cout << "\n";
}

// Calculate the trace of a 2-dimensional matrix
void Calculate2(CovMatrix2& matrix2, const int& N, const int& runs)
{
  // Define a counter
  int counter2 = 0;

  // Define a double to hold the trace value
  double trace2 = 0.;

  // Loop "run" times
  for (int r=0; r<runs; ++r)
  {
    // Calculate the trace
    for (int i=0; i<N; ++i)
    {
      trace2 = matrix2.trace();
      ++counter2;
    }

    // Write the output
    cout << counter2 << " trace calculations took place\n";
    cout << "Trace =\n" << trace2 << "\n";

    // Reset the counter
    counter2 = 0;
  }
}
 
// Calculate the trace of a 3-dimensional matrix
void Calculate3(CovMatrix3& matrix3, const int& N, const int& runs)
{
  // Define a counter
  int counter3 = 0;
  
  // Define a double to hold the trace value
  double trace3 = 0.;

  // Loop "run" times
  for (int r=0; r<runs; ++r)
  {
    // Calculate the trace
    for (int i=0; i<N; ++i)
    {
      trace3 = matrix3.trace();
      ++counter3;
    }

    // Write the output
    cout << counter3 << " trace calculations took place\n";  
    cout << "Trace =\n" << trace3 << "\n";
  
    // Reset the counter
    counter3 = 0;
  }
}

// Calculate the trace of a 4-dimensional matrix
void Calculate4(CovMatrix4& matrix4, const int& N, const int& runs)
{
  // Define a counter
  int counter4 = 0;
  
  // Define a double to hold the trace value
  double trace4 = 0.;

  // Loop "run" times
  for (int r=0; r<runs; ++r)
  {
    // Calculate the trace
    for (int i=0; i<N; ++i)
    {
      trace4 = matrix4.trace();
      ++counter4;
    }

    // Write the output
    cout << counter4 << " trace calculations took place\n";
    cout << "Trace =\n" << trace4 << "\n";

    // Reset the counter
    counter4 = 0;
  }
}

// Calculate the trace of a 5-dimensional matrix
void Calculate5(CovMatrix5& matrix5, const int& N, const int& runs)
{
  // Define a counter
  int counter5 = 0;
  
  // Define a double to hold the trace value
  double trace5 = 0.;

  // Loop "run" times
  for (int r=0; r<runs; ++r)
  {
    // Calculate the trace
    for (int i=0; i<N; ++i)
    {
      trace5 = matrix5.trace();
      ++counter5;
    }

    // Write the output
    cout << counter5 << " trace calculations took place\n";
    cout << "Trace =\n" << trace5 << "\n";

    // Reset the counter
    counter5 = 0;  
  }
}

// Calculate the trace of a 6-dimensional matrix
void Calculate6(CovMatrix6& matrix6, const int& N, const int& runs)
{
  // Define a counter
  int counter6 = 0;

  // Define a double to hold the trace value
  double trace6 = 0.;

  // Loop "run" times
  for (int r=0; r<runs; ++r)
  {
    // Calculate the trace
    for (int i=0; i<N; ++i)
    {
      trace6 = matrix6.trace();
      ++counter6;
    }

    // Write the output
    cout << counter6 << " trace calculations took place\n";
    cout << "Trace =\n" << trace6 << "\n";

    // Reset the counter
    counter6 = 0;
  }
}


int main()
{
  // What package are we using?
  Identification();

  // Get the dim(-ension)
  int dim = 0; GetDimension(dim);

  // Define the "dim x dim" CovMatrices
  CovMatrix2 m2;
  CovMatrix3 m3;
  CovMatrix4 m4;      
  CovMatrix5 m5;
  CovMatrix6 m6;

  // Describe how the matrix was defined
  Describe (dim);
  
  // Build the "dim x dim" matrix
  switch(dim)
  {
    case 2:
      BuildCovMatrix2(m2);
      break;
    case 3:
      BuildCovMatrix3(m3); 
      break;
    case 4:
      BuildCovMatrix4(m4);
      break;
    case 5:
      BuildCovMatrix5(m5);
      break;
    case 6:
      BuildCovMatrix6(m6);
      break;
  }

  // Show the actual matrix
  switch(dim)
  {
    case 2:
      Show2(m2);
      break;
    case 3:
      Show3(m3);
      break;
    case 4:
      Show4(m4);
      break;
    case 5:
      Show5(m5);
      break;
    case 6:
      Show6(m6);
      break;
  } 

  // Calculation parameters
  int tim = 0; GetCalcs(tim); // number of calculations
  int run = 0; GetRuns(run);  // number of runs  

  // Do the calculations
  switch(dim)
  {
    case 2:
      Calculate2(m2, tim, run);
      break;
    case 3:
      Calculate3(m3, tim, run);
      break;
    case 4:
      Calculate4(m4, tim, run);
      break;
    case 5:
      Calculate5(m5, tim, run);
      break;
    case 6:
      Calculate6(m6, tim, run);
      break;
  }

  return 0;
}