// Validity test for 5x5 similarity transforms
#include <iostream>
#include "CovMatrices/Similarity.h"

using namespace CovMatrices;
using namespace std;

CovMatrix C(7);
std::vector<double> v;
std::vector< std::vector<double> > vv, vv2;
CovMatrix C2;

int main()
{
  C(0,0)=7.;
  C(1,0)=6.;
  C(1,1)=6.1;
  C(0,2)=5.;
  C(1,2)=5.1;
  C(2,2)=5.2;
  C(0,3)=4.;
  C(1,3)=4.1;
  C(2,3)=4.2;
  C(3,3)=4.3;
  C(0,4)=3.;
  C(1,4)=3.1;
  C(2,4)=3.2;
  C(3,4)=3.3;
  C(4,4)=3.4;
  C(0,5)=2.;
  C(1,5)=2.1;
  C(2,5)=2.2;
  C(3,5)=2.3;
  C(4,5)=2.4;
  C(5,5)=2.5;
  C(0,6)=1.;
  C(1,6)=1.1;
  C(2,6)=1.2;
  C(3,6)=1.3;
  C(4,6)=1.4;
  C(5,6)=1.5;
  C(6,6)=1.6;
  cout << C;
  cout << C.determinant() << std::endl;
  C.diagonalize(v,vv);
  // NOTE: diagonalize produces an array of eigenvectors
  // When used as a matrix, the array is the TRANSPOSE of the
  // similarity matrix needed to diagonalize the original matrix
  CovMatrix::printv(cout, v);
  cout << "*** " << v[0]*v[1]*v[2]*v[3]*v[4]*v[5]*v[6] << " ***" << std::endl;
  // The following similarity transformations should produce
  // (approximate) diagonal matrices with the eigenvectors as
  // the diagonal elements.
  cout << "*** Trying similarityPAPt ***\n";
  similarityPAPt(C,vv,C2);
  cout << C2 << std::endl;
  vv2.resize(7);
  for (int i = 0; i < 7; ++i)
  {
    vv2[i].resize(7);
    for (int j = 0; j < 7; ++j)
    {
      vv2[i][j] = vv[j][i];
    }
  }
  cout << "*** Trying similarityPtAP ***\n";
  similarityPtAP(C,vv2,C2);
  cout << C2 << std::endl;
  return 0;
}