IMP logo
IMP Reference Guide  2.20.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-2022 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 #include <cereal/access.hpp>
19 
20 IMPALGEBRA_BEGIN_NAMESPACE
21 
22 //! Represent an eigen analysis of some data.
23 template <int D>
25  public:
28  VectorD<D> values, VectorD<D> centroid)
29  : eigen_vecs_(pcs), eigen_values_(values), centroid_(centroid) {}
31  IMP_USAGE_CHECK(!eigen_vecs_.empty(), "The PCA was not initialized");
32  return eigen_vecs_;
33  }
34  VectorD<D> get_principal_component(unsigned int i) const {
35  IMP_USAGE_CHECK(!eigen_vecs_.empty(), "The PCA was not initialized");
36  return eigen_vecs_[i];
37  }
38  VectorD<D> get_principal_values() const {
39  IMP_USAGE_CHECK(!eigen_vecs_.empty(), "The PCA was not initialized");
40  return eigen_values_;
41  }
42  double get_principal_value(unsigned int i) const {
43  IMP_USAGE_CHECK(!eigen_vecs_.empty(), "The PCA was not initialized");
44  return eigen_values_[i];
45  }
46  inline VectorD<D> get_centroid() const { return centroid_; }
47  void set_centroid(VectorD<D> cntr) {
48  IMP_USAGE_CHECK(!eigen_vecs_.empty(), "The PCA was not initialized");
49  centroid_ = cntr;
50  }
53 
54  private:
55  int compare(const PrincipalComponentAnalysisD &o) const {
56  IMP_USAGE_CHECK(eigen_vecs_.empty() || o.eigen_vecs_.empty(),
57  "Cannot compare against anything other than the default"
58  " PrincipalComponentAnalysis");
59  if (eigen_vecs_.empty() && o.eigen_vecs_.empty()) {
60  return 0;
61  } else {
62  return -1;
63  }
64  }
65  Vector<VectorD<D> > eigen_vecs_;
66  VectorD<D> eigen_values_;
67  VectorD<D> centroid_;
68 
69  friend class cereal::access;
70 
71  template<class Archive> void serialize(Archive &ar) {
72  ar(eigen_vecs_, eigen_values_, centroid_);
73  }
74 };
75 
76 #ifndef IMP_DOXYGEN
77 typedef PrincipalComponentAnalysisD<1> PrincipalComponentAnalysis1D;
78 typedef PrincipalComponentAnalysisD<2> PrincipalComponentAnalysis2D;
79 typedef PrincipalComponentAnalysisD<3> PrincipalComponentAnalysis3D;
80 typedef PrincipalComponentAnalysisD<4> PrincipalComponentAnalysis4D;
81 typedef PrincipalComponentAnalysisD<5> PrincipalComponentAnalysis5D;
82 typedef PrincipalComponentAnalysisD<6> PrincipalComponentAnalysis6D;
83 typedef PrincipalComponentAnalysisD<-1> PrincipalComponentAnalysisKD;
85  PrincipalComponentAnalysis1Ds;
87  PrincipalComponentAnalysis2Ds;
89  PrincipalComponentAnalysis3Ds;
91  PrincipalComponentAnalysis4Ds;
93  PrincipalComponentAnalysis5Ds;
95  PrincipalComponentAnalysis6Ds;
97  PrincipalComponentAnalysisKDs;
98 
99 template <int D>
100 inline void PrincipalComponentAnalysisD<D>::show(std::ostream &out) const {
101  if (eigen_vecs_.empty()) {
102  out << "invalid";
103  return;
104  }
105  out << "vectors: " << eigen_vecs_ << " weights: " << eigen_values_
106  << " centroid: " << centroid_ << std::endl;
107 }
108 
109 #endif
110 
111 //! Perform principal components analysis on a set of vectors
112 /** \see PrincipalComponentAnalysis
113  */
114 template <int D>
116  const Vector<VectorD<D> > &ps) {
117  IMP_USAGE_CHECK(!ps.empty(), "Need some vectors to get components.");
118  unsigned int dim = ps[0].get_dimension();
119  VectorD<D> m =
120  std::accumulate(ps.begin(), ps.end(), get_zero_vector_kd(dim)) /
121  ps.size();
122  Eigen::MatrixXd cov = internal::get_covariance_matrix(ps, m);
123  IMP_LOG_VERBOSE("The covariance matrix is " << cov << std::endl);
124 
125  Eigen::JacobiSVD<Eigen::MatrixXd> svd = cov.jacobiSvd(Eigen::ComputeFullV);
126  Eigen::MatrixXd V = svd.matrixV();
127  Eigen::VectorXd SV = svd.singularValues();
128 
129  IMP_LOG_VERBOSE("V is " << V << std::endl);
130  VectorD<D> values = ps[0];
131  Vector<VectorD<D> > vectors(dim, values);
132  for (unsigned int i = 0; i < dim; ++i) {
133  values[i] = SV[i];
134  for (unsigned int j = 0; j < dim; ++j) {
135  vectors[i][j] = V(j, i);
136  }
137  }
138  // the principal components are the columns of V
139  // pc1(pc3) is the vector of the largest(smallest) eigenvalue
140  return PrincipalComponentAnalysisD<D>(vectors, values, m);
141 }
142 
143 //! Get all alignments of the first principal component system to the second one
145  const PrincipalComponentAnalysisD<3> &pca1,
146  const PrincipalComponentAnalysisD<3> &pca2);
147 
148 #ifndef IMP_DOXYGEN
149 typedef PrincipalComponentAnalysisD<3> PrincipalComponentAnalysis;
150 #endif
151 
152 IMPALGEBRA_END_NAMESPACE
153 
154 #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:280
#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:42
Represent an eigen analysis of some data.
Base class for geometric types.
A Cartesian vector in D-dimensions.
Definition: VectorD.h:56
int compare(const VectorD< D > &a, const VectorD< D > &b)
lexicographic comparison of two vectors
Definition: VectorD.h:183
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.