IMP logo
IMP Reference Guide  2.15.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-2021 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/log.h>
17 #include <IMP/log_macros.h>
18 
19 IMPALGEBRA_BEGIN_NAMESPACE
20 
21 //! Represent an eigen analysis of some data.
22 template <int D>
24  public:
27  VectorD<D> values, VectorD<D> centroid)
28  : eigen_vecs_(pcs), eigen_values_(values), centroid_(centroid) {}
30  IMP_USAGE_CHECK(!eigen_vecs_.empty(), "The PCA was not initialized");
31  return eigen_vecs_;
32  }
33  VectorD<D> get_principal_component(unsigned int i) const {
34  IMP_USAGE_CHECK(!eigen_vecs_.empty(), "The PCA was not initialized");
35  return eigen_vecs_[i];
36  }
37  VectorD<D> get_principal_values() const {
38  IMP_USAGE_CHECK(!eigen_vecs_.empty(), "The PCA was not initialized");
39  return eigen_values_;
40  }
41  double get_principal_value(unsigned int i) const {
42  IMP_USAGE_CHECK(!eigen_vecs_.empty(), "The PCA was not initialized");
43  return eigen_values_[i];
44  }
45  inline VectorD<D> get_centroid() const { return centroid_; }
46  void set_centroid(VectorD<D> cntr) {
47  IMP_USAGE_CHECK(!eigen_vecs_.empty(), "The PCA was not initialized");
48  centroid_ = cntr;
49  }
52 
53  private:
54  int compare(const PrincipalComponentAnalysisD &o) const {
55  IMP_USAGE_CHECK(eigen_vecs_.empty() || o.eigen_vecs_.empty(),
56  "Cannot compare against anything other than the default"
57  " PrincipalComponentAnalysis");
58  if (eigen_vecs_.empty() && o.eigen_vecs_.empty()) {
59  return 0;
60  } else {
61  return -1;
62  }
63  }
64  Vector<VectorD<D> > eigen_vecs_;
65  VectorD<D> eigen_values_;
66  VectorD<D> centroid_;
67 };
68 
69 #ifndef IMP_DOXYGEN
70 typedef PrincipalComponentAnalysisD<1> PrincipalComponentAnalysis1D;
71 typedef PrincipalComponentAnalysisD<2> PrincipalComponentAnalysis2D;
72 typedef PrincipalComponentAnalysisD<3> PrincipalComponentAnalysis3D;
73 typedef PrincipalComponentAnalysisD<4> PrincipalComponentAnalysis4D;
74 typedef PrincipalComponentAnalysisD<5> PrincipalComponentAnalysis5D;
75 typedef PrincipalComponentAnalysisD<6> PrincipalComponentAnalysis6D;
76 typedef PrincipalComponentAnalysisD<-1> PrincipalComponentAnalysisKD;
78  PrincipalComponentAnalysis1Ds;
80  PrincipalComponentAnalysis2Ds;
82  PrincipalComponentAnalysis3Ds;
84  PrincipalComponentAnalysis4Ds;
86  PrincipalComponentAnalysis5Ds;
88  PrincipalComponentAnalysis6Ds;
90  PrincipalComponentAnalysisKDs;
91 
92 template <int D>
93 inline void PrincipalComponentAnalysisD<D>::show(std::ostream &out) const {
94  if (eigen_vecs_.empty()) {
95  out << "invalid";
96  return;
97  }
98  out << "vectors: " << eigen_vecs_ << " weights: " << eigen_values_
99  << " centroid: " << centroid_ << std::endl;
100 }
101 
102 #endif
103 
104 //! Perform principal components analysis on a set of vectors
105 /** \see PrincipalComponentAnalysis
106  */
107 template <int D>
109  const Vector<VectorD<D> > &ps) {
110  IMP_USAGE_CHECK(!ps.empty(), "Need some vectors to get components.");
111  unsigned int dim = ps[0].get_dimension();
112  VectorD<D> m =
113  std::accumulate(ps.begin(), ps.end(), get_zero_vector_kd(dim)) /
114  ps.size();
115  Eigen::MatrixXd cov = internal::get_covariance_matrix(ps, m);
116  IMP_LOG_VERBOSE("The covariance matrix is " << cov << std::endl);
117 
118  Eigen::JacobiSVD<Eigen::MatrixXd> svd = cov.jacobiSvd(Eigen::ComputeFullV);
119  Eigen::MatrixXd V = svd.matrixV();
120  Eigen::VectorXd SV = svd.singularValues();
121 
122  IMP_LOG_VERBOSE("V is " << V << std::endl);
123  VectorD<D> values = ps[0];
124  Vector<VectorD<D> > vectors(dim, values);
125  for (unsigned int i = 0; i < dim; ++i) {
126  values[i] = SV[i];
127  for (unsigned int j = 0; j < dim; ++j) {
128  vectors[i][j] = V(j, i);
129  }
130  }
131  // the principal components are the columns of V
132  // pc1(pc3) is the vector of the largest(smallest) eigenvalue
133  return PrincipalComponentAnalysisD<D>(vectors, values, m);
134 }
135 
136 //! Get all alignments of the first principal component system to the second one
138  const PrincipalComponentAnalysisD<3> &pca1,
139  const PrincipalComponentAnalysisD<3> &pca2);
140 
141 #ifndef IMP_DOXYGEN
142 typedef PrincipalComponentAnalysisD<3> PrincipalComponentAnalysis;
143 #endif
144 
145 IMPALGEBRA_END_NAMESPACE
146 
147 #endif /* IMPALGEBRA_EIGEN_ANALYSIS_H */
Base class for geometric types.
VectorD< D > get_zero_vector_kd(int Di)
Return a dynamically sized vector of zeros.
Definition: VectorD.h:276
#define IMP_SHOWABLE(Name)
#define IMP_COMPARISONS(Name)
Implement comparison in a class using a compare function.
#define IMP_LOG_VERBOSE(expr)
Definition: log_macros.h:83
A more IMP-like version of the std::vector.
Definition: Vector.h:40
Represent an eigen analysis of some data.
Base class for geometric types.
A Cartesian vector in D-dimensions.
Definition: VectorD.h:52
int compare(const VectorD< D > &a, const VectorD< D > &b)
lexicographic comparison of two vectors
Definition: VectorD.h:179
Logging and error reporting support.
Simple D vector class.
Transformation3Ds get_alignments_from_first_to_second(const PrincipalComponentAnalysisD< 3 > &pca1, const PrincipalComponentAnalysisD< 3 > &pca2)
Get all alignments of the first principal component system to the second one.
PrincipalComponentAnalysisD< D > get_principal_components(const Vector< VectorD< D > > &ps)
Perform principal components analysis on a set of vectors.
Simple 3D transformation class.
Vector3D get_centroid(const Vector3Ds &ps)
Return the centroid of a set of vectors.
Definition: Vector3D.h:68
std::ostream & show(Hierarchy h, std::ostream &out=std::cout)
Print the hierarchy using a given decorator to display each node.
IMP::Vector< Transformation3D > Transformation3Ds
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
Definition: check_macros.h:168
A reference frame in 3D.
Logging and error reporting support.