cppcode/mdmath.cpp

Go to the documentation of this file.

/*
 * Copyright (c) 2009-2010 Jose-Luis Lisani <joseluis.lisani@uib.es>
 * All rights reserved.
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#include "../Eigen/Dense"

using namespace Eigen;



void compute_pca3(double C[3][3], double vd[3], double PCs[3][3])
{
  MatrixXd CC(3, 3);
  for (int i=0; i < 3; i++)
    for (int j=0; j < 3; j++) 
          CC(i, j)=C[i][j];

  //compute eigenvalues and eigenvectors

  //using EigenSolver
  EigenSolver<MatrixXd> diag(CC);
  VectorXcd lambda = diag.eigenvalues();
  MatrixXcd P = diag.eigenvectors();

  /*
  //using SVD
  MatrixXd YY = XX.transpose()/sqrt((double) n);
  MatrixXd P;
  JacobiSVD<MatrixXd> svd(YY, ComputeThinU | ComputeThinV);
  P=svd.matrixV();
  MatrixXd D=P.transpose()*CC*P;
  VectorXd lambda=D.diagonal();
  */


  //sort eigenvalues and eigenvectors 
  //and fill output data structures
  VectorXd real_lambda;
  real_lambda=lambda.array().real();

  MatrixXd real_P(3, 3);
  real_P=P.array().real();

  int position[3], pos;
  double ev[3], max;

  for (int k=0; k < 3; k++) ev[k]=real_lambda(k);
  for (int k=0; k < 3; k++) {
    max=ev[0];
    pos=0;
    for (int j=1; j < 3; j++) {
      if (ev[j] > max) {          
        max=ev[j];
        pos=j;
      }
    }
    position[k]=pos;
    ev[pos]=0;
  }
  
  for (int k=0; k < 3; k++) {
    vd[k]=real_lambda(position[k]);
    PCs[0][k]=real_P(0, position[k]);
    PCs[1][k]=real_P(1, position[k]);
    PCs[2][k]=real_P(2, position[k]);
  }

  
  return;
}

---