//
// This file is part of MorphoGraphX - http://www.MorphoGraphX.org
// Copyright (C) 2012-2015 Richard S. Smith and collaborators.
//
// If you use MorphoGraphX in your work, please cite:
//   http://dx.doi.org/10.7554/eLife.05864
//
// MorphoGraphX is free software, and is licensed under under the terms of the 
// GNU General (GPL) Public License version 2.0, http://www.gnu.org/licenses.
// 
#include "PCA.hpp"
#include "Information.hpp"
#include "Progress.hpp"


namespace mgx 
{

bool PCA::operator()(const Matrix3f& corr)
    {
      //Information::out << "Decompose correlation matrix" << endl;
      // Eigen-decomposition of the matrices
      for(int i = 0 ; i < 3 ; ++i)
        for(int j = 0 ; j < 3 ; ++j)
          gsl_matrix_set(mat, i, j, corr(i,j));
      gsl_eigen_symmv(mat, eval, evec, w);
      gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_VAL_DESC);

      p1 = Point3f(gsl_matrix_get(evec, 0, 0), gsl_matrix_get(evec, 1, 0), gsl_matrix_get(evec, 2, 0));
      p2 = Point3f(gsl_matrix_get(evec, 0, 1), gsl_matrix_get(evec, 1, 1), gsl_matrix_get(evec, 2, 1));
      p3 = Point3f(gsl_matrix_get(evec, 0, 2), gsl_matrix_get(evec, 1, 2), gsl_matrix_get(evec, 2, 2));
      ev = Point3f(gsl_vector_get(eval, 0), gsl_vector_get(eval, 1), gsl_vector_get(eval, 2));
      if((p1 ^ p2) * p3 < 0)
        p3 = -p3;
      return true;
    }



}
