LOOS  v2.3.2
vsa-lib.cpp
1 /*
2  This file is part of LOOS.
3 
4  LOOS (Lightweight Object-Oriented Structure library)
5  Copyright (c) 2010 Tod D. Romo
6  Department of Biochemistry and Biophysics
7  School of Medicine & Dentistry, University of Rochester
8 
9  This package (LOOS) is free software: you can redistribute it and/or modify
10  it under the terms of the GNU General Public License as published by
11  the Free Software Foundation under version 3 of the License.
12 
13  This package is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with this program. If not, see <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 
24 
25 #include "vsa-lib.hpp"
26 
27 using namespace std;
28 using namespace loos;
29 
30 
31 namespace ENM {
32 
33  boost::tuple<DoubleMatrix, DoubleMatrix> VSA::eigenDecomp(DoubleMatrix& A, DoubleMatrix& B) {
34 
35  DoubleMatrix AA = A.copy();
36  DoubleMatrix BB = B.copy();
37 
38  f77int itype = 1;
39  char jobz = 'V';
40  char uplo = 'U';
41  char range = 'I';
42  f77int n = AA.rows();
43  f77int lda = n;
44  f77int ldb = n;
45  double vl = 0.0;
46  double vu = 0.0;
47  f77int il = 7;
48  f77int iu = n;
49 
50  char dpar = 'S';
51  double abstol = 2.0 * dlamch_(&dpar);
52  //double abstol = -1.0;
53 
54  f77int m;
55  DoubleMatrix W(n, 1);
56  DoubleMatrix Z(n, n);
57  f77int ldz = n;
58 
59  f77int lwork = -1;
60  f77int info;
61  double *work = new double[1];
62 
63  f77int *iwork = new f77int[5*n];
64  f77int *ifail = new f77int[n];
65 
66  dsygvx_(&itype, &jobz, &range, &uplo, &n, AA.get(), &lda, BB.get(), &ldb, &vl, &vu, &il, &iu, &abstol, &m, W.get(), Z.get(), &ldz, work, &lwork, iwork, ifail, &info);
67  if (info != 0) {
68  cerr << "ERROR- dsygvx returned " << info << endl;
69  exit(-1);
70  }
71 
72  lwork = work[0];
73  if (verbosity_ > 1)
74  std::cerr << "dsygvx requested " << work[0] /(1024.0*1024.0) << " MB\n";
75 
76  delete[] work;
77  work = new double[lwork];
78 
79  dsygvx_(&itype, &jobz, &range, &uplo, &n, AA.get(), &lda, BB.get(), &ldb, &vl, &vu, &il, &iu, &abstol, &m, W.get(), Z.get(), &ldz, work, &lwork, iwork, ifail, &info);
80  if (info != 0) {
81  cerr << "ERROR- dsygvx returned " << info << endl;
82  exit(-1);
83  }
84 
85  if (m != n-6) {
86  cerr << "ERROR- only got " << m << " eigenpairs instead of " << n-6 << endl;
87  exit(-10);
88  }
89 
90  vector<uint> indices = sortedIndex(W);
91  W = permuteRows(W, indices);
92  Z = permuteColumns(Z, indices);
93 
94  boost::tuple<DoubleMatrix, DoubleMatrix> result(W, Z);
95  return(result);
96 
97  }
98 
99 
100 
101  // Mass-weight eigenvectors
102  DoubleMatrix VSA::massWeight(DoubleMatrix& U, DoubleMatrix& M) {
103 
104  // First, compute the cholesky decomp of M (i.e. it's square-root)
105  DoubleMatrix R = M.copy();
106  char uplo = 'U';
107  f77int n = M.rows();
108  f77int lda = n;
109  f77int info;
110  dpotrf_(&uplo, &n, R.get(), &lda, &info);
111  if (info != 0) {
112  cerr << "ERROR- dpotrf() returned " << info << endl;
113  exit(-1);
114  }
115 
116  if (debugging_)
117  writeAsciiMatrix(prefix_ + "_R.asc", R, meta_, false);
118 
119  // Now multiply M * U
120  DoubleMatrix UU = U.copy();
121  f77int m = U.rows();
122  n = U.cols();
123  double alpha = 1.0;
124  f77int ldb = m;
125 
126 #if defined(__linux__) || defined(__CYGWIN__) || defined(__FreeBSD__)
127  char side = 'L';
128  char notrans = 'N';
129  char diag = 'N';
130 
131  dtrmm_(&side, &uplo, &notrans, &diag, &m, &n, &alpha, R.get(), &lda, UU.get(), &ldb);
132 
133 #else
134 
135  cblas_dtrmm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, m, n, alpha, R.get(), lda, UU.get(), ldb);
136 #endif
137 
138  normalizeColumns(UU);
139  return(UU);
140  }
141 
142 
143 
144  void VSA::solve() {
145 
146  if (verbosity_ > 1)
147  std::cerr << "Building hessian...\n";
148  buildHessian();
149 
150  uint n = hessian_.cols();
151  uint l = subset_size_ * 3;
152 
153  DoubleMatrix Hss = submatrix(hessian_, Math::Range(0,l), Math::Range(0,l));
154  DoubleMatrix Hee = submatrix(hessian_, Math::Range(l, n), Math::Range(l, n));
155  DoubleMatrix Hse = submatrix(hessian_, Math::Range(0,l), Math::Range(l, n));
156  DoubleMatrix Hes = submatrix(hessian_, Math::Range(l, n), Math::Range(0, l));
157 
158  if (debugging_) {
159  writeAsciiMatrix(prefix_ + "_H.asc", hessian_, meta_, false);
160  writeAsciiMatrix(prefix_ + "_Hss.asc", Hss, meta_, false);
161  writeAsciiMatrix(prefix_ + "_Hee.asc", Hee, meta_, false);
162  writeAsciiMatrix(prefix_ + "_Hse.asc", Hse, meta_, false);
163  }
164 
165  if (verbosity_ > 1)
166  std::cerr << "Inverting environment hessian...\n";
167  DoubleMatrix Heei = Math::invert(Hee);
168 
169  // Build the effective Hessian
170  if (verbosity_ > 1)
171  std::cerr << "Computing effective hessian...\n";
172 
173  Hssp_ = Hss - Hse * Heei * Hes;
174 
175  if (debugging_)
176  writeAsciiMatrix(prefix_ + "_Hssp.asc", Hssp_, meta_, false);
177 
178  // Shunt in the event of using unit masses... We can use the SVD to
179  // to get the eigenpairs from Hssp
180  if (masses_.rows() == 0) {
181  Timer<> t;
182 
183  if (verbosity_ > 0)
184  std::cerr << "Calculating SVD of effective hessian...\n";
185 
186  t.start();
187  boost::tuple<DoubleMatrix, DoubleMatrix, DoubleMatrix> svdresult = svd(Hssp_);
188  t.stop();
189 
190  if (verbosity_ > 0)
191  std::cerr << "SVD took " << loos::timeAsString(t.elapsed()) << std::endl;
192 
193  eigenvecs_ = boost::get<0>(svdresult);
194  eigenvals_ = boost::get<1>(svdresult);
195 
196  reverseColumns(eigenvecs_);
197  reverseRows(eigenvals_);
198  return;
199 
200  }
201 
202 
203  // Build the effective mass matrix
204  DoubleMatrix Ms = submatrix(masses_, Math::Range(0, l), Math::Range(0, l));
205  DoubleMatrix Me = submatrix(masses_, Math::Range(l, n), Math::Range(l, n));
206 
207  if (verbosity_ > 1)
208  std::cerr << "Computing effective mass matrix...\n";
209  Msp_ = Ms + Hse * Heei * Me * Heei * Hes;
210 
211  if (debugging_) {
212  writeAsciiMatrix(prefix_ + "_Ms.asc", Ms, meta_, false);
213  writeAsciiMatrix(prefix_ + "_Me.asc", Me, meta_, false);
214  writeAsciiMatrix(prefix_ + "_Msp.asc", Msp_, meta_, false);
215  }
216 
217  // Run the eigen-decomposition...
218  boost::tuple<DoubleMatrix, DoubleMatrix> eigenpairs;
219  Timer<> t;
220  if (verbosity_ > 0)
221  std::cerr << "Computing eigendecomposition...\n";
222 
223  t.start();
224  eigenpairs = eigenDecomp(Hssp_, Msp_);
225  t.stop();
226 
227  if (verbosity_ > 0)
228  std::cerr << "Eigendecomposition took " << loos::timeAsString(t.elapsed()) << std::endl;
229 
230  eigenvals_ = boost::get<0>(eigenpairs);
231  DoubleMatrix Us = boost::get<1>(eigenpairs);
232 
233  // Need to mass-weight the eigenvectors so they're orthogonal in R3...
234  eigenvecs_ = massWeight(Us, Msp_);
235  }
236 
237 
238 };
239 
Simple matrix template class using policy classes to determine behavior.
Definition: MatrixImpl.hpp:53
T permuteRows(const T &A, const std::vector< uint > &indices)
Returns a copy of the matrix with the rows permuted by the indices.
Definition: MatrixOps.hpp:156
STL namespace.
double elapsed(void) const
Return the elapsed time.
Definition: timer.hpp:109
void start(void)
Starts the timer.
Definition: timer.hpp:85
std::ostream & writeAsciiMatrix(std::ostream &os, const Math::Matrix< T, P, S > &M, const std::string &meta, const Math::Range &start, const Math::Range &end, const bool trans=false, F fmt=F())
Write a submatrix to a stream.
Namespace to encapsulate Elastic Network Model routines.
Definition: anm-lib.hpp:32
std::pair< uint, uint > Range
Specify a range for columns/rows [first,second)
Definition: MatrixImpl.hpp:49
std::string timeAsString(const double t, const uint precision)
Convert t (seconds) into a string, converting to hours and minutes as necessary.
Definition: utils.cpp:212
double stop(void)
Stops the timer.
Definition: timer.hpp:94
Class for tracking time.
Definition: timer.hpp:80
T permuteColumns(const T &A, const std::vector< uint > &indices)
Returns a copy of the matrix with the columns permuted by the indices.
Definition: MatrixOps.hpp:138
std::vector< uint > sortedIndex(const T &A)
Definition: sorting.hpp:71
Namespace for most things not already encapsulated within a class.
Matrix< T, OrderPolicy, StoragePolicy > copy(void) const
Deep copy...
Definition: MatrixImpl.hpp:170
void normalizeColumns(DoubleMatrix &A)
Normalizes each column as a column-vector.
Definition: MatrixOps.cpp:418