IMP  2.4.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-2015 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 #include <IMP/base/log_macros.h>
19 
20 IMPALGEBRA_BEGIN_NAMESPACE
21 /** Represent a eigen analysis of some data.
22  */
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  base::Vector<VectorD<D> > eigen_vecs_;
66  VectorD<D> eigen_values_;
67  VectorD<D> centroid_;
68 };
69 
70 #ifndef IMP_DOXYGEN
71 typedef PrincipalComponentAnalysisD<1> PrincipalComponentAnalysis1D;
72 typedef PrincipalComponentAnalysisD<2> PrincipalComponentAnalysis2D;
73 typedef PrincipalComponentAnalysisD<3> PrincipalComponentAnalysis3D;
74 typedef PrincipalComponentAnalysisD<4> PrincipalComponentAnalysis4D;
75 typedef PrincipalComponentAnalysisD<5> PrincipalComponentAnalysis5D;
76 typedef PrincipalComponentAnalysisD<6> PrincipalComponentAnalysis6D;
77 typedef PrincipalComponentAnalysisD<-1> PrincipalComponentAnalysisKD;
79  PrincipalComponentAnalysis1Ds;
81  PrincipalComponentAnalysis2Ds;
83  PrincipalComponentAnalysis3Ds;
85  PrincipalComponentAnalysis4Ds;
87  PrincipalComponentAnalysis5Ds;
89  PrincipalComponentAnalysis6Ds;
91  PrincipalComponentAnalysisKDs;
92 
93 template <int D>
94 inline void PrincipalComponentAnalysisD<D>::show(std::ostream &out) const {
95  if (eigen_vecs_.empty()) {
96  out << "invalid";
97  return;
98  }
99  out << "vectors: " << eigen_vecs_ << " weights: " << eigen_values_
100  << " centroid: " << centroid_ << std::endl;
101 }
102 
103 #endif
104 
105 //! Perform principal components analysis on a set of vectors
106 /** \see PrincipalComponentAnalysis
107  */
108 template <int D>
110  const base::Vector<VectorD<D> > &ps) {
111  IMP_USAGE_CHECK(!ps.empty(), "Need some vectors to get components.");
112  unsigned int dim = ps[0].get_dimension();
113  VectorD<D> m =
114  std::accumulate(ps.begin(), ps.end(), get_zero_vector_kd(dim)) /
115  ps.size();
116  internal::TNT::Array2D<double> cov = internal::get_covariance_matrix(ps, m);
117  IMP_LOG_VERBOSE("The covariance matrix is " << cov << std::endl);
118  internal::JAMA::SVD<double> svd(cov);
119  internal::TNT::Array2D<double> V(dim, dim);
120  internal::TNT::Array1D<double> SV;
121 
122  svd.getV(V);
123  IMP_LOG_VERBOSE("V is " << V << std::endl);
124  svd.getSingularValues(SV);
125  VectorD<D> values = ps[0];
126  base::Vector<VectorD<D> > vectors(dim, values);
127  for (unsigned int i = 0; i < dim; ++i) {
128  values[i] = SV[i];
129  for (unsigned int j = 0; j < dim; ++j) {
130  vectors[i][j] = V[j][i];
131  }
132  }
133  // the principal components are the columns of V
134  // pc1(pc3) is the vector of the largest(smallest) eigenvalue
135  return PrincipalComponentAnalysisD<D>(vectors, values, m);
136 }
137 
138 //! Get all possible alignments of the first principal
139 //! component system to the second one
141  const PrincipalComponentAnalysisD<3> &pca1,
142  const PrincipalComponentAnalysisD<3> &pca2);
143 
144 #ifndef IMP_DOXYGEN
145 typedef PrincipalComponentAnalysisD<3> PrincipalComponentAnalysis;
146 #endif
147 
148 IMPALGEBRA_END_NAMESPACE
149 #endif /* IMPALGEBRA_EIGEN_ANALYSIS_H */
Basic types used by IMP.
VectorD< D > get_zero_vector_kd(int Di)
Return a dynamically sized vector of zeros.
Definition: VectorD.h:262
PrincipalComponentAnalysisD< D > get_principal_components(const base::Vector< VectorD< D > > &ps)
Perform principal components analysis on a set of vectors.
#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:94
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.
IMP::base::Vector< Transformation3D > Transformation3Ds
Simple D vector class.
Transformation3Ds get_alignments_from_first_to_second(const PrincipalComponentAnalysisD< 3 > &pca1, const PrincipalComponentAnalysisD< 3 > &pca2)
Simple 3D transformation class.
Vector3D get_centroid(const Vector3Ds &ps)
Returns the centroid of a set of vectors.
Definition: Vector3D.h:56
void show(Hierarchy h, std::ostream &out=std::cout)
Print out a molecular hierarchy.
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
Definition: check_macros.h:170
Simple 3D rotation class.
Logging and error reporting support.