Ignore:
Timestamp:
Feb 1, 2008, 5:15:44 AM (14 years ago)
Author:
Peter
Message:

merging branch peter-dev into trunk delta 1008:994

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/yat/utility/PCA.cc

    r1000 r1009  
    2929#include "SVD.h"
    3030#include "utility.h"
     31#include "vectorView.h"
    3132
    3233#include <iostream>
     
    7374    eigenvectors_.clone(U);
    7475    eigenvectors_ .transpose();
     76    eigenvalues_ = pSVD->s();
     77
     78    // T
     79    for( size_t i = 0; i < eigenvalues_.size(); ++i )
     80      eigenvalues_(i) = eigenvalues_(i)*eigenvalues_(i);
     81    eigenvalues_ *= 1.0/A_center.rows();
     82
     83    // Sort the eigenvectors in order of eigenvalues
     84    // Simple (not efficient) algorithm that always
     85    // make sure that the i:th element is in its correct
     86    // position (N element --> Ordo( N*N ))
     87    for ( size_t i = 0; i < eigenvalues_.size(); ++i )
     88      for ( size_t j = i + 1; j < eigenvalues_.size(); ++j )
     89        if ( eigenvalues_(j) > eigenvalues_(i) ) {
     90          std::swap( eigenvalues_(i), eigenvalues_(j) );
     91          eigenvectors_.swap_rows(i,j);
     92        }
     93  }
     94
     95
     96  /*
     97  void PCA::process_transposed_problem(void)
     98  {
     99    // Row-center the data matrix
     100    utility::matrix A_center( A_.rows(), A_.columns() );
     101    this->row_center( A_center );
     102
     103    // Transform into SVD friendly dimension
     104    A_.transpose();
     105    A_center.transpose();
     106
     107    // Single value decompose the data matrix
     108    std::auto_ptr<SVD> pSVD( new SVD( A_center ) );
     109    pSVD->decompose();
     110    utility::matrix U(pSVD->U());
     111    utility::matrix V(pSVD->V());
     112
     113    // Read the eigenvectors and eigenvalues
     114    eigenvectors_.clone(V);
     115    eigenvectors_.transpose();
    75116    eigenvalues_.clone(pSVD->s());
     117
     118    // Transform back when done with SVD!
     119    // (used V insted of U now for eigenvectors)
     120    A_.transpose();
     121    A_center.transpose();
    76122
    77123    // T
     
    91137        }
    92138  }
    93 
    94 
    95   /*
    96   void PCA::process_transposed_problem(void)
    97   {
    98     // Row-center the data matrix
    99     utility::matrix A_center( A_.rows(), A_.columns() );
    100     this->row_center( A_center );
    101 
    102     // Transform into SVD friendly dimension
    103     A_.transpose();
    104     A_center.transpose();
    105 
    106     // Single value decompose the data matrix
    107     std::auto_ptr<SVD> pSVD( new SVD( A_center ) );
    108     pSVD->decompose();
    109     utility::matrix U(pSVD->U());
    110     utility::matrix V(pSVD->V());
    111 
    112     // Read the eigenvectors and eigenvalues
    113     eigenvectors_.clone(V);
    114     eigenvectors_.transpose();
    115     eigenvalues_.clone(pSVD->s());
    116 
    117     // Transform back when done with SVD!
    118     // (used V insted of U now for eigenvectors)
    119     A_.transpose();
    120     A_center.transpose();
    121 
    122     // T
    123     for( size_t i = 0; i < eigenvalues_.size(); ++i )
    124       eigenvalues_[ i ] = eigenvalues_[ i ]*eigenvalues_[ i ];
    125     eigenvalues_ *= 1.0/A_center.rows();
    126 
    127     // Sort the eigenvectors in order of eigenvalues
    128     // Simple (not efficient) algorithm that always
    129     // make sure that the i:th element is in its correct
    130     // position (N element --> Ordo( N*N ))
    131     for ( size_t i = 0; i < eigenvalues_.size(); ++i )
    132       for ( size_t j = i + 1; j < eigenvalues_.size(); ++j )
    133         if ( eigenvalues_[ j ] > eigenvalues_[ i ] ) {
    134           std::swap( eigenvalues_[ i ], eigenvalues_[ j ] );
    135           eigenvectors_.swap_rows(i,j);
    136         }
    137   }
    138139  */
    139140
     
    189190  void PCA::row_center(utility::matrix& A_center)
    190191  {
    191     meanvalues_.clone(utility::vector(A_.rows()));
     192    meanvalues_ = vector(A_.rows());
    192193    utility::vector A_row_sum(A_.rows());
    193     for (size_t i=0; i<A_row_sum.size(); ++i)
    194       A_row_sum(i)=utility::sum(utility::vector(A_,i));
     194    for (size_t i=0; i<A_row_sum.size(); ++i){
     195      const vectorView tmp(A_.row_vec(i));
     196      A_row_sum(i) = sum(tmp);
     197    }
    195198    for( size_t i = 0; i < A_center.rows(); ++i ) {
    196       meanvalues_[i] = A_row_sum(i) / static_cast<double>(A_.columns());
     199      meanvalues_(i) = A_row_sum(i) / static_cast<double>(A_.columns());
    197200      for( size_t j = 0; j < A_center.columns(); ++j )
    198201        A_center(i,j) = A_(i,j) - meanvalues_(i);
Note: See TracChangeset for help on using the changeset viewer.