IMP  2.0.0
The Integrative Modeling Platform
MSConnectivityRestraint.h
Go to the documentation of this file.
1 /**
2  * \file IMP/core/MSConnectivityRestraint.h
3  * \brief Mass Spec Connectivity restraint.
4  *
5  * Restrict max distance between at least one pair of particles of any
6  * two distinct types. It also handles multiple copies of the same particles.
7  *
8  * Copyright 2007-2013 IMP Inventors. All rights reserved.
9  *
10  */
11 
12 #ifndef IMPCORE_MS_CONNECTIVITY_RESTRAINT_H
13 #define IMPCORE_MS_CONNECTIVITY_RESTRAINT_H
14 
15 #include <vector>
16 #include <string>
17 #include <IMP/core/core_config.h>
18 #include "DistanceRestraint.h"
19 
20 #include <IMP/SingletonContainer.h>
21 #include <IMP/Restraint.h>
22 #include <IMP/PairScore.h>
23 
24 IMPCORE_BEGIN_NAMESPACE
25 
26 
27 // Hopefully this would "fix" the problem with friend access
28 // in older versions of gcc
29 #ifdef __GNUC__
30  #define GCC_VERSION (__GNUC__ * 10000 \
31  + __GNUC_MINOR__ * 100 \
32  + __GNUC_PATCHLEVEL__)
33  #if GCC_VERSION >= 40100
34  #define IMPCORE_FRIEND_IS_OK
35  #endif
36 #else
37  // assume that other compilers are fine
38  #define IMPCORE_FRIEND_IS_OK
39 #endif
40 
41 
42 
43 
44 //! Ensure that a set of particles remains connected with one another.
45 /** The restraint implements ambiguous connectivity. That is, it takes
46  several particles including multiple copies and ensures that they remain
47  connected, but allows how they are connected to change. If you wish
48  to restrain the connectivity of sets of
49  particles (i.e. each protein is represented using a set of balls)
50  use an appropriate PairScore which calls a Refiner (such
51  as ClosePairsPairScore).
52 
53  \verbinclude ms_connectivity_restraint.py
54 
55  More precisely, the restraint scores by computing the MST on the complete
56  graph connecting all the particles. The edge weights are given by
57  the value of the PairScore for the two endpoints of the edge.
58  */
59 class IMPCOREEXPORT MSConnectivityRestraint : public Restraint
60 {
61  IMP::OwnerPointer<PairScore> ps_;
62  IMP::OwnerPointer<SingletonContainer> sc_;
63  double eps_;
64 public:
65  //! Use the given PairScore
66  /** If sc is nullptr, a ListSingletonContainer is created internally.
67  eps is set to 0.1 by default.
68  */
69  MSConnectivityRestraint(PairScore* ps, double eps=0.1);
70  /** @name Particles to be connected
71 
72  The following methods are used to manipulate the list of particles
73  that are to be connected. Each particle should have all the
74  attributes expected by the PairScore used.
75 
76  Ideally, one should pass a singleton container instead. These
77  can only be used if none is passed.
78  */
79  /*@{*/
80  unsigned int add_type(const ParticlesTemp &ps);
81  unsigned int add_composite(const Ints &components);
82  unsigned int add_composite(const Ints &components, unsigned int parent);
83  //void add_particle(Particle *p);
84  //void add_particles(const Particles &ps);
85  //void set_particles(const Particles &ps);
86  /*@}*/
87 
88  //! Return the set of pairs which are connected by the restraint
89  /** This set of pairs reflects the current configuration at the time of
90  the get_connected_pairs() call, not the set at the time of the last
91  evaluate() call.
92  */
93  ParticlePairsTemp get_connected_pairs() const;
94 
95  //! Return the pair score used for scoring
97  return ps_;
98  }
99 
101 
102  double unprotected_evaluate(IMP::kernel::DerivativeAccumulator *accum) const;
103 
104  ModelObjectsTemp do_get_inputs() const;
105 
107 
108 #ifdef IMPCORE_FRIEND_IS_OK
109  private:
110 #endif
111 
112  class ParticleMatrix
113  {
114  public:
115 
116  class ParticleData
117  {
118  public:
119  ParticleData(Particle *p, unsigned int id)
120  : particle_(p)
121  , id_(id)
122  {}
123 
124  Particle *get_particle() const
125  {
126  return particle_;
127  }
128 
129  unsigned int get_id() const
130  {
131  return id_;
132  }
133  private:
134  Particle *particle_;
135  unsigned int id_;
136  };
137 
138  ParticleMatrix(unsigned int number_of_classes)
139  : protein_by_class_(number_of_classes)
140  , min_distance_(std::numeric_limits<double>::max())
141  , max_distance_(0)
142  , current_id_(0)
143  {}
144 
145  ParticleMatrix()
146  : min_distance_(std::numeric_limits<double>::max())
147  , max_distance_(0)
148  , current_id_(0)
149  {}
150 
151  void resize(unsigned int number_of_classes)
152  {
153  protein_by_class_.resize(number_of_classes);
154  }
155 
156  unsigned int add_particle(Particle *p, unsigned int id);
157  unsigned int add_type(const ParticlesTemp &ps);
158  void create_distance_matrix(const PairScore *ps);
159  void clear_particles()
160  {
161  particles_.clear();
162  for ( unsigned int i = 0; i < protein_by_class_.size(); ++i )
163  protein_by_class_[i].clear();
164  }
165  unsigned int size() const
166  {
167  return particles_.size();
168  }
169  unsigned int get_number_of_classes() const
170  {
171  return protein_by_class_.size();
172  }
173  double get_distance(unsigned int p1, unsigned int p2) const
174  {
175  return dist_matrix_[p1*size() + p2];
176  }
177  Ints const &get_ordered_neighbors(unsigned int p) const
178  {
179  return order_[p];
180  }
181  ParticleData const &get_particle(unsigned int p) const
182  {
183  return particles_[p];
184  }
185  Ints const &get_all_proteins_in_class(
186  unsigned int id) const
187  {
188  return protein_by_class_[id];
189  }
190  double max_distance() const
191  {
192  return max_distance_;
193  }
194  double min_distance() const
195  {
196  return min_distance_;
197  }
198  private:
199  class DistCompare
200  {
201  public:
202  DistCompare(unsigned int source, ParticleMatrix const &parent)
203  : parent_(parent)
204  , source_(source)
205  {}
206 
207  bool operator()(unsigned int p1, unsigned int p2) const
208  {
209  return parent_.get_distance(source_, p1) <
210  parent_.get_distance(source_, p2);
211  }
212  private:
213  ParticleMatrix const &parent_;
214  unsigned int source_;
215  };
216 
217  base::Vector<ParticleData> particles_;
218  Floats dist_matrix_;
219  base::Vector< Ints > order_;
220  base::Vector< Ints > protein_by_class_;
221  double min_distance_;
222  double max_distance_;
223  unsigned int current_id_;
224  };
225 
226  class ExperimentalTree
227  {
228  public:
229  ExperimentalTree()
230  : root_(-1)
231  , finalized_(false)
232  {}
233 
234  void connect(unsigned int parent, unsigned int child);
235  void finalize();
236  unsigned int add_composite(const Ints &components);
237  unsigned int add_composite(const Ints &components, unsigned int parent);
238 
239  class Node
240  {
241  public:
242  Node()
243  : visited_(false)
244  {}
245  unsigned int get_number_of_parents() const
246  {
247  return parents_.size();
248  }
249  bool is_root() const
250  {
251  return get_number_of_parents() == 0;
252  }
253  unsigned int get_number_of_children() const
254  {
255  return children_.size();
256  }
257  bool is_leaf() const
258  {
259  return get_number_of_children() == 0;
260  }
261  unsigned int get_parent(unsigned int idx) const
262  {
263  return parents_[idx];
264  }
265  unsigned int get_child(unsigned int idx) const
266  {
267  return children_[idx];
268  }
269  typedef base::Vector< std::pair<unsigned int, int> > Label;
270  const Label &get_label() const
271  {
272  return label_;
273  }
274 
275 
276  base::Vector<unsigned int> parents_;
277  base::Vector<unsigned int> children_;
278  Label label_;
279  bool visited_;
280  };
281 
282  bool find_cycle(unsigned int node_index);
283  bool is_consistent(unsigned int node_index) const;
284 
285  const Node *get_node(unsigned int index) const
286  {
287  return &nodes_[index];
288  }
289  unsigned int get_number_of_nodes() const
290  {
291  return nodes_.size();
292  }
293  unsigned int get_root() const
294  {
295  return root_;
296  }
297  void desc_to_label(const Ints &components,
298  Node::Label &label);
299 
300  base::Vector<Node> nodes_;
301  unsigned int root_;
302  bool finalized_;
303  };
304 
305  ParticleMatrix particle_matrix_;
306  mutable ExperimentalTree tree_;
307 
308  friend class MSConnectivityScore;
309 
310 };
311 
312 IMPCORE_END_NAMESPACE
313 
314 #endif /* IMPCORE_MS_CONNECTIVITY_RESTRAINT_H */