IMP  2.0.0
The Integrative Modeling Platform
eigen_analysis.h
Go to the documentation of this file.
1 /**
2  * \file IMP/algebra/eigen_analysis.h
3  * \brief Principal component analysis of a set of points
4  *
5  * Copyright 2007-2013 IMP Inventors. All rights reserved.
6  */
7 
8 #ifndef IMPALGEBRA_EIGEN_ANALYSIS_H
9 #define IMPALGEBRA_EIGEN_ANALYSIS_H
10 
11 #include "VectorD.h"
12 #include "Transformation3D.h"
13 #include "GeometricPrimitiveD.h"
14 #include "IMP/algebra/internal/utility.h"
16 #include <IMP/algebra/internal/jama_svd.h>
17 #include <IMP/base/log.h>
18 
19 IMPALGEBRA_BEGIN_NAMESPACE
20 /** Represent a eigen analysis of some data.
21  */
22 template <int D>
24  public GeometricPrimitiveD<D> {
25  public:
28  const base::Vector< VectorD<D> > &pcs,VectorD<D> values,
29  VectorD<D> centroid) : eigen_vecs_(pcs),
30  eigen_values_(values), centroid_(centroid){
31  }
33  IMP_USAGE_CHECK(!eigen_vecs_.empty(),
34  "The PCA was not initialized");
35  return eigen_vecs_;
36  }
37  VectorD<D> get_principal_component(unsigned int i) const {
38  IMP_USAGE_CHECK(!eigen_vecs_.empty(),
39  "The PCA was not initialized");
40  return eigen_vecs_[i];
41  }
42  VectorD<D> get_principal_values() const {
43  IMP_USAGE_CHECK(!eigen_vecs_.empty(),
44  "The PCA was not initialized");
45  return eigen_values_;
46  }
47  double get_principal_value(unsigned int i) const {
48  IMP_USAGE_CHECK(!eigen_vecs_.empty(),
49  "The PCA was not initialized");
50  return eigen_values_[i];
51  }
52  inline VectorD<D> get_centroid() const {
53  return centroid_;
54  }
55  void set_centroid(VectorD<D> cntr) {
56  IMP_USAGE_CHECK(!eigen_vecs_.empty(),
57  "The PCA was not initialized");
58  centroid_=cntr;
59  }
62 private:
63  int compare(const PrincipalComponentAnalysisD &o) const {
64  IMP_UNUSED(o);
65  IMP_USAGE_CHECK(!eigen_vecs_.empty() && !o.eigen_vecs_.empty(),
66  "Cannot compare against anything other than the default"
67  " PrincipalComponentAnalysis");
68  if (eigen_vecs_.empty() && eigen_vecs_.empty()) {
69  return 0;
70  } else {
71  return -1;
72  }
73  }
74  base::Vector<VectorD<D> > eigen_vecs_;
75  VectorD<D> eigen_values_;
76  VectorD<D> centroid_;
77 };
78 
79 #ifndef IMP_DOXYGEN
80 typedef PrincipalComponentAnalysisD<1> PrincipalComponentAnalysis1D;
81 typedef PrincipalComponentAnalysisD<2> PrincipalComponentAnalysis2D;
82 typedef PrincipalComponentAnalysisD<3> PrincipalComponentAnalysis3D;
83 typedef PrincipalComponentAnalysisD<4> PrincipalComponentAnalysis4D;
84 typedef PrincipalComponentAnalysisD<5> PrincipalComponentAnalysis5D;
85 typedef PrincipalComponentAnalysisD<6> PrincipalComponentAnalysis6D;
86 typedef PrincipalComponentAnalysisD<-1> PrincipalComponentAnalysisKD;
88 PrincipalComponentAnalysis1Ds;
90 PrincipalComponentAnalysis2Ds;
92 PrincipalComponentAnalysis3Ds;
94 PrincipalComponentAnalysis4Ds;
96 PrincipalComponentAnalysis5Ds;
98 PrincipalComponentAnalysis6Ds;
100 PrincipalComponentAnalysisKDs;
101 
102 template <int D>
103 inline void PrincipalComponentAnalysisD<D>::show(std::ostream& out) const {
104  if (eigen_vecs_.empty()) {
105  out << "invalid";
106  return;
107  }
108  out << "vectors: " << eigen_vecs_ << " weights: " << eigen_values_
109  << " centroid: " << centroid_ << std::endl;
110 }
111 
112 #endif
113 
114 //! Perform principal components analysis on a set of vectors
115 /** \relatesalso PrincipalComponentAnalysis
116  */
117 template <int D>
119  const base::Vector<VectorD<D> > &ps) {
120  IMP_USAGE_CHECK(!ps.empty(), "Need some vectors to get components.");
121  unsigned int dim=ps[0].get_dimension();
122  VectorD<D> m = std::accumulate(ps.begin(), ps.end(),
123  get_zero_vector_kd(dim))/ps.size();
124  internal::TNT::Array2D<double> cov = internal::get_covariance_matrix(ps,m);
125  IMP_LOG_VERBOSE( "The covariance matrix is " << cov << std::endl);
126  internal::JAMA::SVD<double> svd(cov);
127  internal::TNT::Array2D<double> V(dim, dim);
128  internal::TNT::Array1D<double> SV;
129 
130  svd.getV(V);
131  IMP_LOG_VERBOSE( "V is " << V << std::endl);
132  svd.getSingularValues(SV);
133  VectorD<D> values= ps[0];
134  base::Vector<VectorD<D> > vectors(dim, values);
135  for (unsigned int i=0; i< dim; ++i) {
136  values[i]= SV[i];
137  for (unsigned int j=0; j< dim; ++j) {
138  vectors[i][j]= V[j][i];
139  }
140  }
141  //the principal components are the columns of V
142  //pc1(pc3) is the vector of the largest(smallest) eigenvalue
143  return PrincipalComponentAnalysisD<D>(vectors, values, m);
144 }
145 
146 //! Get all possible alignments of the first principal
147 //! component system to the second one
149  const PrincipalComponentAnalysisD<3> &pca1,
150  const PrincipalComponentAnalysisD<3> &pca2);
151 
152 #ifndef IMP_DOXYGEN
153 typedef PrincipalComponentAnalysisD<3> PrincipalComponentAnalysis;
154 #endif
155 
156 IMPALGEBRA_END_NAMESPACE
157 #endif /* IMPALGEBRA_EIGEN_ANALYSIS_H */