LOOS  v2.3.2
bcomlib.hpp
1 /*
2  bcomlib
3 */
4 
5 
6 
7 
8 /*
9 
10  This file is part of LOOS.
11 
12  LOOS (Lightweight Object-Oriented Structure library)
13  Copyright (c) 2010, Tod D. Romo
14  Department of Biochemistry and Biophysics
15  School of Medicine & Dentistry, University of Rochester
16 
17  This package (LOOS) is free software: you can redistribute it and/or modify
18  it under the terms of the GNU General Public License as published by
19  the Free Software Foundation under version 3 of the License.
20 
21  This package is distributed in the hope that it will be useful,
22  but WITHOUT ANY WARRANTY; without even the implied warranty of
23  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24  GNU General Public License for more details.
25 
26  You should have received a copy of the GNU General Public License
27  along with this program. If not, see <http://www.gnu.org/licenses/>.
28 */
29 
30 // @cond PACKAGES_INTERNAL
31 
32 
33 #if !defined(LOOS_BCOMLIB_HPP)
34 #define LOOS_BCOMLIB_HPP
35 
36 
37 #include <loos.hpp>
38 
39 /*
40  * = Developer's Note =
41  *
42  * Some of the following is vestigial from the develpment of the BCOM method.
43  * These interfaces and classes may go away or change substantially in
44  * the next release or two, so caveat programmer...comments may also
45  * be quite sparse...
46  */
47 
48 
49 
50 namespace Convergence {
51 
52  // Treats an AtomicGroup as a column vector and subtracts it from
53  // each column of the matrix M.
54  template<typename T>
55  void subtractStructure(T& M, const loos::AtomicGroup& model) {
56  std::vector<float> avg(model.size() * 3);
57  uint k = 0;
58  for (uint i=0; i<model.size(); ++i) {
59  loos::GCoord c = model[i]->coords();
60  avg[k++] = c.x();
61  avg[k++] = c.y();
62  avg[k++] = c.z();
63  }
64 
65  for (uint i=0; i<M.cols(); ++i)
66  for (uint j=0; j<M.rows(); ++j)
67  M(j, i) -= avg[j];
68  }
69 
70  // Computes the cosine content for a col-vector
71  template<typename T>
72  double cosineContent(T& V, const uint col) {
73  double sum1 = 0;
74  double sum2 = 0;
75 
76  uint m = V.rows();
77  double k = (col+1) * M_PI / m;
78  for (uint j=0; j<m; ++j) {
79  sum1 += cos(k * j) * V(j, col);
80  sum2 += V(j, col) * V(j, col);
81  }
82  double c = 2.0 * sum1 * sum1 / (sum2 * m);
83  return(c);
84  }
85 
86 
87 
88  /*
89  * Various policies that determine how blocks are extracted and
90  * averaged/aligned. The idea is that they are really functors
91  * which, given a vector<AtomicGroup> ensemble, will extract a
92  * RealMatrix of coordinates where each AtomicGroup is a column
93  * vector. The appropriate processing (i.e. average subtraction) is
94  * also performed by the functor.
95  *
96  * local_average, when set, means that the average of the ensemble
97  * is used rather than the average passed to the constructor
98  * (presumably, the average of the -entire- trajectory)
99  */
100 
101 
102  // Align to the passed structure
103  struct AlignToPolicy {
104  AlignToPolicy(const loos::AtomicGroup& targ) : target(targ), local_average(true) { }
105  AlignToPolicy(const loos::AtomicGroup& targ, const bool flag) : target(targ), local_average(flag) { }
106 
107  loos::RealMatrix operator()(std::vector<loos::AtomicGroup>& ensemble) {
108  for (std::vector<loos::AtomicGroup>::iterator i = ensemble.begin(); i != ensemble.end(); ++i)
109  (*i).alignOnto(target);
110 
111  loos::RealMatrix M = loos::extractCoords(ensemble);
112  if (local_average) {
114  subtractStructure(M, avg);
115  } else
116  subtractStructure(M, target);
117 
118  return(M);
119  }
120 
121  loos::AtomicGroup target;
122  bool local_average;
123  };
124 
125 
126  // Do no aligning
127  struct NoAlignPolicy {
128  NoAlignPolicy() : local_average(true) { }
129  NoAlignPolicy(const loos::AtomicGroup& avg_) : avg(avg_), local_average(false) { }
130  NoAlignPolicy(const loos::AtomicGroup& avg_, const bool flag) : avg(avg_), local_average(flag) { }
131 
132  loos::RealMatrix operator()(std::vector<loos::AtomicGroup>& ensemble) {
133 
134  loos::RealMatrix M = loos::extractCoords(ensemble);
135  if (local_average) {
136  loos::AtomicGroup lavg = loos::averageStructure(ensemble);
137  subtractStructure(M, lavg);
138  } else
139  subtractStructure(M, avg);
140  return(M);
141  }
142 
143  loos::AtomicGroup avg;
144  bool local_average;
145  };
146 
147 
148 
149 
150 
151  // Compute the PCA of an ensemble using the specified coordinate
152  // extraction policy...
153  //
154 
155  template<class ExtractPolicy>
156  boost::tuple<loos::RealMatrix, loos::RealMatrix> pca(std::vector<loos::AtomicGroup>& ensemble, ExtractPolicy& extractor) {
157 
158  loos::RealMatrix M = extractor(ensemble);
159  loos::RealMatrix C = loos::Math::MMMultiply(M, M, false, true);
160 
161  // Compute [U,D] = eig(C)
162  char jobz = 'V';
163  char uplo = 'L';
164  f77int n = M.rows();
165  f77int lda = n;
166  float dummy;
167  loos::RealMatrix W(n, 1);
168  f77int lwork = -1;
169  f77int info;
170  ssyev_(&jobz, &uplo, &n, C.get(), &lda, W.get(), &dummy, &lwork, &info);
171  if (info != 0)
172  throw(loos::NumericalError("ssyev failed in loos::pca()", info));
173 
174 
175  lwork = static_cast<f77int>(dummy);
176  float *work = new float[lwork+1];
177 
178  ssyev_(&jobz, &uplo, &n, C.get(), &lda, W.get(), work, &lwork, &info);
179  if (info != 0)
180  throw(loos::NumericalError("ssyev failed in loos::pca()", info));
181 
182  reverseColumns(C);
183  reverseRows(W);
184 
185  // Zap negative eigenvalues...
186  for (uint j=0; j<W.rows(); ++j)
187  if (W[j] < 0.0)
188  W[j] = 0.0;
189 
190  boost::tuple<loos::RealMatrix, loos::RealMatrix> result(W, C);
191  return(result);
192 
193  }
194 
195 
196 
197  // Get just the RSVs (this is for cosine-content calculations)
198  // given an extraction policy...
199  //
200 
201  template<class ExtractPolicy>
202  loos::RealMatrix rsv(std::vector<loos::AtomicGroup>& ensemble, ExtractPolicy& extractor) {
203 
204  loos::RealMatrix M = extractor(ensemble);
205  loos::RealMatrix C = loos::Math::MMMultiply(M, M, false, true);
206 
207  // Compute [U,D] = eig(C)
208  char jobz = 'V';
209  char uplo = 'L';
210  f77int n = M.rows();
211  f77int lda = n;
212  float dummy;
213  loos::RealMatrix W(n, 1);
214  f77int lwork = -1;
215  f77int info;
216  ssyev_(&jobz, &uplo, &n, C.get(), &lda, W.get(), &dummy, &lwork, &info);
217  if (info != 0)
218  throw(loos::NumericalError("ssyev failed in loos::pca()", info));
219 
220 
221  lwork = static_cast<f77int>(dummy);
222  float *work = new float[lwork+1];
223 
224  ssyev_(&jobz, &uplo, &n, C.get(), &lda, W.get(), work, &lwork, &info);
225  if (info != 0)
226  throw(loos::NumericalError("ssyev failed in loos::pca()", info));
227 
228  reverseColumns(C);
229  reverseRows(W);
230 
231  // Correctly scale the eigenvalues
232  for (uint j=0; j<W.rows(); ++j)
233  W[j] = W[j] < 0 ? 0.0 : sqrt(W[j]);
234 
235  // Multiply eigenvectors by inverse eigenvalues
236  for (uint i=0; i<C.cols(); ++i) {
237  double konst = (W[i] > 0.0) ? (1.0/W[i]) : 0.0;
238 
239  for (uint j=0; j<C.rows(); ++j)
240  C(j, i) *= konst;
241  }
242 
243  W.reset();
244  loos::RealMatrix Vt = loos::Math::MMMultiply(C, M, true, false);
246  return(V);
247  }
248 
249 
250 }
251 
252 #endif
253 
254 
255 // @endcond PACKAGES_INTERNAL
T transpose(const T &A)
Transposes a matrix.
Definition: MatrixOps.hpp:174
Simple matrix template class using policy classes to determine behavior.
Definition: MatrixImpl.hpp:53
Exception caused by a blas/atlas error.
Definition: exceptions.hpp:91
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
RealMatrix MMMultiply(const RealMatrix &A, const RealMatrix &B, const bool transa, const bool transb)
Matrix-matrix multiply (using BLAS)
Definition: MatrixOps.cpp:170
AtomicGroup averageStructure(const std::vector< AtomicGroup > &ensemble)
Compute the average structure of a set of AtomicGroup objects.
Definition: ensembles.cpp:36
void reset(void)
Deallocate data...
Definition: MatrixImpl.hpp:195