LOOS  v2.3.2
ensembles.cpp
1 /*
2  This file is part of LOOS.
3 
4  LOOS (Lightweight Object-Oriented Structure library)
5  Copyright (c) 2008, Tod D. Romo, Alan Grossfield
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 
26 #include <ensembles.hpp>
27 #include <XForm.hpp>
28 #include <AtomicGroup.hpp>
29 #include <Trajectory.hpp>
30 #include <alignment.hpp>
31 
32 namespace loos {
33 
34  // Assume all groups are already sorted or matched...
35 
36  AtomicGroup averageStructure(const std::vector<AtomicGroup>& ensemble) {
37  AtomicGroup avg = ensemble[0].copy();
38 
39  // First, zap our coords...
40  int n = avg.size();
41  int i;
42  for (i=0; i<n; i++)
43  avg[i]->coords() = GCoord(0.0, 0.0, 0.0);
44 
45  // Now, accumulate...
46  std::vector<AtomicGroup>::const_iterator j;
47  for (j = ensemble.begin(); j != ensemble.end(); j++) {
48  for (i = 0; i<n; i++)
49  avg[i]->coords() += (*j)[i]->coords();
50  }
51 
52  for (i=0; i<n; i++)
53  avg[i]->coords() /= ensemble.size();
54 
55  avg.removePeriodicBox();
56  return(avg);
57  }
58 
59 
60 
61  AtomicGroup averageStructure(const std::vector<AtomicGroup>& ensemble, const std::vector<XForm>& xforms) {
62  if (xforms.size() != ensemble.size())
63  throw(LOOSError("Transforms do not match the passed ensemble in loos::averageStructure()"));
64 
65  AtomicGroup avg = ensemble[0].copy();
66 
67  // First, zap our coords...
68  int n = avg.size();
69  int i;
70  for (i=0; i<n; i++)
71  avg[i]->coords() = GCoord(0.0, 0.0, 0.0);
72 
73  // Now, accumulate...
74  for (uint j=0; j<ensemble.size(); ++j) {
75  AtomicGroup frame = ensemble[j].copy();
76  frame.applyTransform(xforms[j]);
77  for (i = 0; i<n; i++)
78  avg[i]->coords() += frame[i]->coords();
79  }
80 
81  for (i=0; i<n; i++)
82  avg[i]->coords() /= ensemble.size();
83 
84  avg.removePeriodicBox();
85  return(avg);
86  }
87 
88 
89 
90 
91  AtomicGroup averageStructure(const AtomicGroup& g, const std::vector<XForm>& xforms, pTraj& traj, const std::vector<uint>& frame_indices) {
92  AtomicGroup avg = g.copy();
93  AtomicGroup frame = g.copy();
94  int n = avg.size();
95  for (int i=0; i<n; i++)
96  avg[i]->coords() = GCoord(0.0, 0.0, 0.0);
97 
98  uint tn = traj->nframes();
99  uint fn = frame_indices.size();
100 
101  // Safety check...
102  if (fn != xforms.size())
103  throw(std::runtime_error("Mismatch in number of frames in the trajectory requested and passed transforms for loos::averageStructure()"));
104 
105  for (uint j=0; j<fn; ++j) {
106  if (frame_indices[j] >= tn)
107  throw(std::runtime_error("Frame index exceeds trajectory size"));
108 
109  traj->readFrame(frame_indices[j]);
110  traj->updateGroupCoords(frame);
111  frame.applyTransform(xforms[j]);
112  for (int i=0; i<n; i++)
113  avg[i]->coords() += frame[i]->coords();
114  }
115 
116  for (int i=0; i<n; i++)
117  avg[i]->coords() /= fn;
118 
119  avg.removePeriodicBox();
120  return(avg);
121  }
122 
123 
124 
125  AtomicGroup averageStructure(const AtomicGroup& g, const std::vector<XForm>& xforms, pTraj& traj) {
126  uint nx = xforms.size();
127  uint nt = traj->nframes();
128 
129  if (nt != nx)
130  throw(std::runtime_error("Mismatch in number of frames in the trajectory and passed transforms"));
131  std::vector<uint> frame_indices;
132  for (uint i=0; i<nt; ++i)
133  frame_indices.push_back(i);
134 
135  return(averageStructure(g, xforms, traj, frame_indices));
136  }
137 
138 
139 
140  void applyTransforms(std::vector<AtomicGroup>& ensemble, std::vector<XForm>& xforms) {
141  uint n = ensemble.size();
142  if (n != xforms.size())
143  throw(std::runtime_error("Mismatch in the size of the ensemble and the transformations"));
144 
145  for (uint i=0; i<n; ++i)
146  ensemble[i].applyTransform(xforms[i]);
147  }
148 
149 
150 
151  void readTrajectory(std::vector<AtomicGroup>& ensemble, const AtomicGroup& model, pTraj trajectory) {
152  AtomicGroup clone = model.copy();
153 
154  while (trajectory->readFrame()) {
155  trajectory->updateGroupCoords(clone);
156  AtomicGroup frame = clone.copy();
157  ensemble.push_back(frame);
158  }
159  }
160 
161 
162  void readTrajectory(std::vector<AtomicGroup>& ensemble, const AtomicGroup& model, pTraj trajectory, std::vector<uint>& frames) {
163  AtomicGroup clone = model.copy();
164 
165  std::vector<uint>::iterator i;
166  for (i = frames.begin(); i != frames.end(); ++i) {
167  if (*i >= trajectory->nframes())
168  throw(std::runtime_error("Frame index exceeds trajectory size in readTrajectory()"));
169  trajectory->readFrame(*i);
170  trajectory->updateGroupCoords(clone);
171  AtomicGroup frame = clone.copy();
172  ensemble.push_back(frame);
173  }
174  }
175 
176 
177 
178  RealMatrix extractCoords(const std::vector<AtomicGroup>& ensemble) {
179  uint n = ensemble.size();
180  uint m = ensemble[0].size();
181  RealMatrix M(3*m, n);
182 
183  for (uint i=0; i<n; ++i)
184  for (uint j=0; j<m; ++j) {
185  M(3*j, i) = ensemble[i][j]->coords().x();
186  M(3*j+1, i) = ensemble[i][j]->coords().y();
187  M(3*j+2, i) = ensemble[i][j]->coords().z();
188  }
189 
190  return(M);
191  }
192 
193 
194  RealMatrix extractCoords(const std::vector<AtomicGroup>& ensemble, const std::vector<XForm>& xforms) {
195  uint n = ensemble.size();
196 
197  if (n != xforms.size())
198  throw(std::runtime_error("Mismatch between the size of the ensemble and the transformations"));
199 
200  uint m = ensemble[0].size();
201  RealMatrix M(3*m, n);
202 
203  for (uint i=0; i<n; ++i) {
204  GMatrix W = xforms[i].current();
205 
206  for (uint j=0; j<m; ++j) {
207  GCoord c = W * ensemble[i][j]->coords();
208  M(3*j, i) = c.x();
209  M(3*j+1, i) = c.y();
210  M(3*j+2, i) = c.z();
211  }
212  }
213 
214  return(M);
215  }
216 
217 
218  void subtractAverage(RealMatrix& M) {
219  uint m = M.rows();
220  uint n = M.cols();
221  double *avg = new double[m];
222 
223  for (uint j=0; j<m; ++j)
224  avg[j] = 0.0;
225 
226  for (uint i=0; i<n; ++i)
227  for (uint j=0; j<m; ++j)
228  avg[j] += M(j, i);
229 
230  for (uint j=0; j<m; ++j)
231  avg[j] /= n;
232 
233  for (uint i=0; i<n; ++i)
234  for (uint j=0; j<m; ++j)
235  M(j,i) -= avg[j];
236  }
237 
238 
239  boost::tuple<RealMatrix, RealMatrix, RealMatrix> svd(std::vector<AtomicGroup>& ensemble, bool align) {
240 
241  RealMatrix M;
242  if (align)
243  iterativeAlignment(ensemble);
244 
245  M = extractCoords(ensemble);
246 
247  subtractAverage(M);
248  boost::tuple<RealMatrix, RealMatrix, RealMatrix> res = Math::svd(M);
249  return(res);
250  }
251 
252 
253  void appendCoords(std::vector< std::vector<double> >& M, AtomicGroup& model, pTraj& traj, const std::vector<uint>& indices, const bool updates = false) {
254 
255  uint l = indices.size();
256  uint n = model.size();
257  uint offset = M.size();
258 
259  M.resize(offset + l, std::vector<double>(3*n));
260 
261  PercentProgressWithTime watcher;
262  PercentTrigger trigger(0.1);
263  ProgressCounter<PercentTrigger, EstimatingCounter> slayer(trigger, EstimatingCounter(l));
264 
265  if (updates) {
266  slayer.attach(&watcher);
267  slayer.start();
268  }
269 
270  for (uint j=0; j<l; ++j) {
271  traj->readFrame(indices[j]);
272  traj->updateGroupCoords(model);
273  if (updates)
274  slayer.update();
275  for (uint i=0; i<n; ++i) {
276  GCoord c = model[i]->coords();
277  M[j + offset][i*3] = c.x();
278  M[j + offset][i*3+1] = c.y();
279  M[j + offset][i*3+2] = c.z();
280  }
281  }
282 
283  if (updates)
284  slayer.finish();
285 
286  }
287 
288 
289 
290 
291  std::vector< std::vector<double> > readCoords(AtomicGroup& model, pTraj& traj, const std::vector<uint>& indices, const bool updates = false) {
292 
293  std::vector< std::vector<double> > M;
294  appendCoords(M, model, traj, indices, updates);
295  return M;
296  }
297 
298 
299 
300 
301 
302 }
boost::tuple< RealMatrix, RealMatrix, RealMatrix > svd(std::vector< AtomicGroup > &ensemble, bool align)
Compute the SVD of an ensemble with optional alignment (note RSVs returned are transposed) ...
Definition: ensembles.cpp:239
Basic 3-D coordinates class.
Definition: Coord.hpp:36
Simple matrix template class using policy classes to determine behavior.
Definition: MatrixImpl.hpp:53
boost::tuple< RealMatrix, RealMatrix, RealMatrix > svd(RealMatrix &M)
Compute the SVD of a single precision matrix.
Definition: MatrixOps.cpp:94
Generic LOOS exception.
Definition: exceptions.hpp:40
void removePeriodicBox()
Remove periodicity.
AtomicGroup copy(void) const
Creates a deep copy of this group.
Definition: AtomicGroup.cpp:56
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
void applyTransform(const XForm &)
Apply the given transform to the group's coordinates...
Namespace for most things not already encapsulated within a class.
AtomicGroup averageStructure(const std::vector< AtomicGroup > &ensemble)
Compute the average structure of a set of AtomicGroup objects.
Definition: ensembles.cpp:36
def extractCoords(traj)
Returns the coordinates for an entire trajectory as an MxN numpy matrix where M is 3*natoms and N is ...
Definition: ensembles.py:35