LOOS  v2.3.2
fid-lib.cpp
1 /*
2  fid-lib
3 
4  Fiducial library
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 
31 
32 #include "fid-lib.hpp"
33 
34 using namespace std;
35 using namespace loos;
36 
37 
38 
39 
40 vecUint findFreeFrames(const vecInt& map) {
41  vecUint indices;
42 
43  for (uint i=0; i<map.size(); ++i)
44  if (map[i] < 0)
45  indices.push_back(i);
46 
47  return(indices);
48 }
49 
50 
51 
52 vecUint assignStructures(AtomicGroup& model, pTraj& traj, const vecUint& frames, const vecGroup& refs) {
53  vecUint assignments(frames.size(), 0);
54  uint j = 0;
55 
56  for (vecUint::const_iterator frame = frames.begin(); frame != frames.end(); ++frame) {
57  traj->readFrame(*frame);
58  traj->updateGroupCoords(model);
59 
60  double mind = numeric_limits<double>::max();
61  uint mini = refs.size() + 1;
62 
63  for (uint i = 0; i < refs.size(); ++i) {
64  model.alignOnto(refs[i]);
65  double d = model.rmsd(refs[i]);
66  if (d < mind) {
67  mind = d;
68  mini = i;
69  }
70  }
71 
72  assignments[j++] = mini;
73  }
74 
75  return(assignments);
76 }
77 
78 
79 vecUint trimFrames(const vecUint& frames, const double frac) {
80  uint bin_size = frac * frames.size();
81  uint remainder = frames.size() - static_cast<uint>(bin_size / frac);
82 
83  vecUint truncated;
84  vecUint::const_iterator vi = frames.begin();
85  for (uint i = 0; i<frames.size() - remainder; ++i)
86  truncated.push_back(*vi++);
87 
88  return(truncated);
89 }
90 
91 
92 
93 
94 boost::tuple<vecGroup, vecUint> pickFiducials(AtomicGroup& model, pTraj& traj, const vecUint& frames, const double f) {
95 
96  // Size of bin
97  uint bin_size = f * frames.size();
98 
99  // Initialize a RNG to use a uniform random distribution.
100  // Use the LOOS generator singleton so the random number stream can
101  // be seeded (or automatically seeded) by LOOS...
102  boost::uniform_real<> rmap;
103  boost::variate_generator< base_generator_type&, boost::uniform_real<> > rng(rng_singleton(), rmap);
104 
105  // Vector of AtomicGroup's representing the fiducial structures
106  vecGroup fiducials;
107 
108  // Track which trajectory frame has been assigned to which fiducial
109  vecInt assignments(frames.size(), -1);
110 
111  // The indices (frame #'s) of the structures picked to be fiducials
112  vecUint refs;
113  vecDouble radii;
114 
115  // Unassigned frames...bootstrap the loop
116  vecUint possible_frames = findFreeFrames(assignments);
117 
118  // Are there any unassigned frames left?
119  while (! possible_frames.empty()) {
120  // Randomly pick one
121  uint pick = possible_frames[static_cast<uint>(floor(possible_frames.size() * rng()))];
122 
123  traj->readFrame(frames[pick]);
124  traj->updateGroupCoords(model);
125 
126  // Make a copy and assign a new bin # to the fiducial
127  AtomicGroup fiducial = model.copy();
128  fiducial.centerAtOrigin();
129  uint myid = fiducials.size();
130  if (assignments[pick] >= 0) {
131  cerr << "INTERNAL ERROR - " << pick << " pick was already assigned to " << assignments[pick] << endl;
132  exit(-99);
133  }
134 
135  fiducials.push_back(fiducial);
136  refs.push_back(pick);
137 
138  // Now find the distance from every unassigned frame to this new fiducial (aligning
139  // them first), and then sort by distance...
140  vector<double> distances(assignments.size(), numeric_limits<double>::max());
141  for (uint i = 0; i<assignments.size(); ++i) {
142  if (assignments[i] >= 0)
143  continue;
144  traj->readFrame(frames[i]);
145  traj->updateGroupCoords(model);
146  model.centerAtOrigin();
147  model.alignOnto(fiducial);
148  distances[i] = model.rmsd(fiducial);
149  }
150 
151  vecUint indices = sortedIndex(distances);
152  uint picked = 0;
153  double maxd = 0.0;
154 
155  // Pick the first bin_size of them (or however many are remaining)
156  // and assign these to the newly picked fiducial
157  for (uint i=0; i<assignments.size() && picked < bin_size; ++i) {
158  if (assignments[indices[i]] < 0) {
159  assignments[indices[i]] = myid;
160  ++picked;
161  if (distances[indices[i]] > maxd)
162  maxd = distances[indices[i]];
163  }
164  }
165  radii.push_back(maxd);
166  possible_frames = findFreeFrames(assignments);
167  }
168 
169  // Safety check...
170  for (vecInt::const_iterator i = assignments.begin(); i != assignments.end(); ++i)
171  if (*i < 0)
172  throw(runtime_error("A frame was not assigned in binFrames()"));
173 
174  boost::tuple<vecGroup, vecUint> result(fiducials, refs);
175  return(result);
176 }
177 
178 
179 
180 
181 int findMaxBin(const vecInt& assignments) {
182 
183  int max_bin = -1;
184  for (vecInt::const_iterator i = assignments.begin(); i != assignments.end(); ++i)
185  if (*i > max_bin)
186  max_bin = *i;
187 
188  return(max_bin);
189 }
190 
191 
192 
193 vecUint histogramBins(const vecInt& assignments) {
194  int max = findMaxBin(assignments);
195  vecUint histogram(max+1, 0);
196 
197  for (vecInt::const_iterator i = assignments.begin(); i != assignments.end(); ++i)
198  if (*i >= 0)
199  histogram[*i] += 1;
200 
201  return(histogram);
202 }
STL namespace.
greal rmsd(const AtomicGroup &)
Compute the RMSD between two groups.
GMatrix alignOnto(const AtomicGroup &)
Superimposes the current group onto the passed group.
Definition: AG_linalg.cpp:196
base_generator_type & rng_singleton(void)
Suite-wide random number generator singleton.
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
std::vector< uint > sortedIndex(const T &A)
Definition: sorting.hpp:71
Namespace for most things not already encapsulated within a class.
GCoord centerAtOrigin(void)
Translates the group so that the centroid is at the origin.