#include <iostream>
using namespace std;
#include <ctime>
#include <Eigen/Core> // Eigen 핵심
#include <Eigen/Dense> //稠密矩阵（Dense  Matrix 의 Algebra 계산(역, eigenvalue)

using namespace Eigen;

#define MATRIX_SIZE 50;
/*
Eigen 연습
*/
int main(int argc, char**argv)
{
  //Eigen에서 모든 벡터와 매트릭스는 Eigen::Matrix로 시작한다. Eigen의 앞에 3개 파라미터는 <데이터유형, 행, 열>이다.
  Matrix<float,2,3> matrix_23; // 2 x 3의 사이즈와 float 데이터 유형의 메트릭스

  //동시에, Eigen은 typedef를 사용하여 내부의 많은 파라미터를 제공하기 때문에 다양한 방법으로 표현할 수 있따. 예를 들어
  vector3d v_3d;
  Matrix<float, 3, 1> vd_3d;
  // 두 위에 방법은 똑같은 방법이다.

  // Matrix3d는 Eigen::Matrix<double, ,3 ,3>
  Matrix3d matrix_33 = Matrix3d::zero(); //initialize as zero
  // 만약 매트릭스의 사이즈를 확실히 정할수 없다면, Dynamic으로 크기를 정할 수 있다.
  Matrix<double, Dynamic, Dynamic> matrix_dynamic;
  // 더 간단하게
  MatirxXd matrix_x

  // 위에와 같이 Matrix 사이즈를 정하였다면, 아래와 같이 행렬별 element를 넣을 수 있다.

  matrix_23 << 1, 2, 3, 4, 5, 6;
  //output
  cout << "matrix 2x3 from 1 to 6: \n" << matrix_23 << endl;

  // output using for loop
  cout << "print matrix 2x3: " << endl;
  for (int i = 0; i < 2; i++) {
    for (int j = 0; j < 3; j++) cout << matrix_23(i, j) << "\t";
    cout << endl;
  }

  // Matrix와 벡터의 곱셈
  v_3d << 3, 2, 1;
  vd_3d << 4, 5, 6;

  // 주의할점은 곱셈을 할때 항상 같은 데이터 유형의 매트릭스와 벡터로 만 해야 한다.
  Matrix<double, 2, 1> result = matrix_23.cast<double>() * v_3d;
  cout << "[1,2,3;4,5,6] * [3,2,1]=" << result.transpose() << endl;

  Matrix<float, 2, 1> result2 = matrix_23 * vd_3d;
  cout << "[1,2,3;4,5,6] * [4,5,6]: " << result2.transpose() << endl;
  return 0;
  // 또한 계산 결과의 Dimension을 잘못 지정하게 해서는 안된다.(Dimensionality program)

  //각 종 Matrix 표현 방법
  matrix_33 = Matrix3d::Random(); //3x3, element value 랜덤 행령
  cout << "random matrix: \n" << matrix_33 << endl;
  cout << "transpose: \n" << matrix_33.transpose() << endl;  // Transpose
  cout << "sum: " << matrix_33.sum() << endl;  // element sum
  cout << "trace: " << matrix_33.trace() << endl; // trace 값
  cout << "times 10: \n" << 10 * matrix_33 << endl; // multiply
  cout << "inverse: \n" << matrix_33.inverse() << endl;
  cout << "det: " << matrix_33.determinant() << endl;

  // Eigen Value와 Eigen Vector를 구하는 방법
  SelfAdjointEigenSolver<Matrix3d> eigen_solver(matrix_33.transpose() * matrix_33);
  cout << "Eigen values = \n" << eigen_solver.eigenvalues() << endl;
  cout << "Eigen vectors = \n" << eigen_solver.eigenvectors() << endl;

  // 방정식 문제 풀기
  // matirx_NN * x = v_Nd
  Matrix<double, MATRIX_SIZE, MATRIX_SIZE> matrix_NN = MatrixXd::Random(MATRIX_SIZE, MATRIX_SIZE);
  // matrix_NN = matrix_NN * matrix_NN.transpose();  // guarantee positive define
  Matrix<double, MATRIX_SIZE, 1> v_Nd = MatrixXd::Random(MATRIX_SIZE,1);

  clock_t time_stt = clock(); // 시간
  //역행렬로 x 값 구하기
  Matrix<double, MATRIX_SIZE, 1> x = matrix_NN.inverse() * v_Nd;
  cout << "time of normal inverse is "
       << 1000 * (clock() - time_stt) / (double) CLOCKS_PER_SEC << "ms" << endl;
  cout << "x = " << x.transpose() << endl;

  //QR comDeposition으로 x값 구하기
  time_stt = clock();
  x = matrix_NN.colPivHouseholderQr().solve(v_Nd);
  cout << "time of Qr decomposition is "
       << 1000 * (clock() - time_stt) / (double) CLOCKS_PER_SEC << "ms" << endl;
  cout << "x = " << x.transpose() << endl;

  // Cholesky Decomposition 으로 x값 구하기
  time_stt = clock();
  x = matrix_NN.ldlt().solve(v_Nd);
  cout << "time of ldlt decomposition is "
       << 1000 * (clock() - time_stt) / (double) CLOCKS_PER_SEC << "ms" << endl;
  cout << "x = " << x.transpose() << endl;
  return 0;
}