IMP  2.2.1
The Integrative Modeling Platform
subset_scores.h
Go to the documentation of this file.
1 /**
2  * \file IMP/domino/subset_scores.h
3  * \brief A beyesian infererence-based sampler.
4  *
5  * Copyright 2007-2014 IMP Inventors. All rights reserved.
6  *
7  */
8 
9 #ifndef IMPDOMINO_SUBSET_SCORES_H
10 #define IMPDOMINO_SUBSET_SCORES_H
11 
12 #include "particle_states.h"
13 #include "Assignment.h"
14 #include "Subset.h"
15 #include "Slice.h"
16 #include "utility.h"
17 #include <IMP/base/log.h>
18 #include <IMP/base/Object.h>
19 #include <IMP/base/cache.h>
20 #include <IMP/kernel/Restraint.h>
21 #include <IMP/base/log.h>
22 #include <boost/unordered_map.hpp>
23 
24 #if IMP_DOMINO_HAS_RMF
25 #include <RMF/HDF5/Group.h>
26 #endif
27 
28 IMPDOMINO_BEGIN_NAMESPACE
29 
30 /** Implement a cache for restraint scores as well as management of restraints
31  for domino.
32 
33  The cache size passed to the constructor is the maximum number of scores
34  that will be saved. A least-recently-used eviction policy is used when
35  that number is exceeded.
36 */
37 class IMPDOMINOEXPORT RestraintCache : public base::Object {
38  IMP_NAMED_TUPLE_2(Key, Keys, base::WeakPointer<kernel::Restraint>, restraint,
39  Assignment, assignment, );
40  IMP_NAMED_TUPLE_3(RestraintData, RestraintDatas,
41  base::PointerMember<ScoringFunction>, scoring_function,
42  Subset, subset, double, max, );
43  IMP_NAMED_TUPLE_2(RestraintSetData, RestraintSetDatas, Slice, slice,
45  IMP_NAMED_TUPLE_2(SetData, SetDatas, RestraintSetDatas, members, double,
46  max, );
47 
48 #ifndef IMP_DOXYGEN
49  class Generator {
50  typedef boost::unordered_map<kernel::Restraint *, RestraintData> RMap;
51  RMap rmap_;
52  typedef boost::unordered_map<kernel::Restraint *, SetData> SMap;
53  SMap sets_;
55 
56  public:
57  Generator(ParticleStatesTable *pst) : pst_(pst) {}
58  typedef double result_type;
59  typedef Key argument_type;
60  template <class Cache>
61  result_type operator()(const argument_type &k, const Cache &cache) const {
62  RMap::const_iterator it = rmap_.find(k.get_restraint());
63  if (it != rmap_.end()) {
64  Subset s = rmap_.find(k.get_restraint())->second.get_subset();
65  load_particle_states(s, k.get_assignment(), pst_);
66  double e;
67  {
69  e = it->second.get_scoring_function()->evaluate_if_below(
70  false, it->second.get_max());
71  }
72  IMP_LOG_TERSE("Restraint " << Showable(k.get_restraint())
73  << " evaluated to " << e << " on "
74  << k.get_assignment() << " vs "
75  << it->second.get_max() << std::endl);
76  // prob can go away with ScoreFunction change
77  if (e > it->second.get_max()) e = std::numeric_limits<double>::max();
78  return e;
79  } else {
80  SMap::const_iterator it = sets_.find(k.get_restraint());
81  IMP_USAGE_CHECK(it != sets_.end(), "Restraint set "
82  << Showable(k.get_restraint())
83  << " not found.");
84  double total = 0;
85  for (unsigned int i = 0; i < it->second.get_members().size(); ++i) {
86  Assignment cur = it->second.get_members()[i].get_slice().get_sliced(
87  k.get_assignment());
88  double score = cache.get(
89  argument_type(it->second.get_members()[i].get_restraint(), cur));
90  total += score * k.get_restraint()->get_weight();
91  if (total >= it->second.get_max()) {
92  break;
93  }
94  }
95  IMP_LOG_TERSE("Restraint " << Showable(k.get_restraint())
96  << " evaluated to " << total << " on "
97  << k.get_assignment() << " with max "
98  << it->second.get_max() << std::endl);
99  if (total >= it->second.get_max()) {
100  return std::numeric_limits<double>::max();
101  } else {
102  return total;
103  }
104  }
105  }
106  void add_to_set(kernel::RestraintSet *rs, kernel::Restraint *r, Slice slice,
107  double max) {
108  IMP_USAGE_CHECK(!dynamic_cast<kernel::RestraintSet *>(r),
109  "don't pass restraint sets here as second arg");
110  sets_[rs].access_members().push_back(RestraintSetData(slice, r));
111  sets_[rs].set_max(max);
112  }
113  void add_restraint(kernel::Restraint *e, Subset s, double max) {
114  IMP_USAGE_CHECK(!dynamic_cast<kernel::RestraintSet *>(e),
115  "don't pass restraint sets here");
116  if (rmap_.find(e) == rmap_.end()) {
117  rmap_[e] = RestraintData(e->create_scoring_function(1.0, max), s, max);
118  } else {
119  IMP_USAGE_CHECK(rmap_.find(e)->second.get_subset() == s,
120  "Subsets don't match on restraint update");
121  rmap_[e].set_max(std::min(rmap_[e].get_max(), max));
122  }
123  }
124  ParticleStatesTable *get_particle_states_table() const { return pst_; }
125  void show_restraint_information(std::ostream &out) const;
126  };
127 #endif // IMP_DOXYGEN
128 
129  struct ApproximatelyEqual {
130  bool operator()(double a, double b) const {
131  return std::abs(a - b) < .1 * (a + b) + .1;
132  }
133  };
134  typedef boost::unordered_map<kernel::Particle *, kernel::ParticlesTemp>
135  DepMap;
136  void add_restraint_internal(kernel::Restraint *r, unsigned int index,
137  kernel::RestraintSet *parent, double max,
138  Subset parent_subset, const DepMap &dependencies);
139  void add_restraint_set_child_internal(kernel::Restraint *r,
140  const Subset &cur_subset,
141  kernel::RestraintSet *parent,
142  double parent_max,
143  Subset parent_subset);
144  void add_restraint_set_internal(kernel::RestraintSet *rs, unsigned int index,
145  const Subset &cur_subset, double cur_max,
146  const DepMap &dependencies);
147  Subset get_subset(kernel::Restraint *r, const DepMap &dependencies) const;
148 // otherwise doxygen seems to index this for some reason
149 #ifndef IMP_DOXYGEN
151  Cache cache_;
152 #endif
153  typedef boost::unordered_map<base::Pointer<kernel::Restraint>, Subset>
154  KnownRestraints;
155  KnownRestraints known_restraints_;
156  // assign a unique index to each restraint for use with I/O
157  typedef boost::unordered_map<base::Pointer<kernel::Restraint>, int>
158  RestraintIndex;
159  RestraintIndex restraint_index_;
160  unsigned int next_index_;
161 
162  public:
164  unsigned int size = std::numeric_limits<unsigned int>::max());
165  /** Recursively process the passed restraints (and sets) so all contained
166  restraints and sets that have maximum are known.*/
168  //! Get the score of a set or restraint
169  /** The returned score will be std::numeric_limits<double>::max()
170  if any of the limits are violated.
171 
172  The assignment passed is the assignment for the particles used by the
173  restraint, not that for some whole subset. See the other get_score()
174  function for a perhaps more useful one.
175  */
176  double get_score(kernel::Restraint *r, const Assignment &a) const {
177  set_was_used(true);
178  double s = cache_.get(Key(r, a));
179  return s;
180  }
181  /** Get the score for a restraint given a subset and assignment on
182  that subset.
183  */
184  double get_score(kernel::Restraint *r, const Subset &s,
185  const Assignment &a) const;
186 
187  //! Make it so kernel::Restraint::get_last_score() returns the score
188  /** This is useful when writing the restraints to disk, as that
189  code often goes off the last score to avoid recomputing the
190  restraints.*/
191  void load_last_score(kernel::Restraint *r, const Subset &s,
192  const Assignment &a);
193  /** Return the restraints that should be evaluated for the subset,
194  given the exclusions.*/
196  const Subsets &exclusions) const;
197 
199 
200 #if IMP_DOMINO_HAS_RMF || defined(IMP_DOXYGEN)
201  /** This assumes that restraints are always added to the cache
202  in the same order.
203  \param[in] particle_ordering An ordering for the particles.
204  \param[in] restraints Which restraints to write out entries for.
205  You probably want to use get_restraints() to generate this.
206  \param[in] max_entries How many entries to write out at most.
207  \param[in] group Where to put the entries.
208  */
209  void save_cache(const kernel::ParticlesTemp &particle_ordering,
210  const kernel::RestraintsTemp &restraints,
211  RMF::HDF5::Group group, unsigned int max_entries);
212  void load_cache(const kernel::ParticlesTemp &ps, RMF::HDF5::ConstGroup group);
213 #endif
214 
215  /** Return the slice for that restraint given the subset. */
216  Slice get_slice(kernel::Restraint *r, const Subset &s) const;
217 
218  /** Return the number of entries currently in the cache.*/
219  unsigned int get_number_of_entries() const { return cache_.size(); }
220 
221  /** Check the entries in the cache.*/
222  void validate() const;
223 
224  /** Print out information about the known restraints and restraint sets.*/
225  void show_restraint_information(std::ostream &out = std::cout) const;
226  double get_hit_rate() const { return cache_.get_hit_rate(); }
228 };
229 
230 IMPDOMINO_END_NAMESPACE
231 
232 #endif /* IMPDOMINO_SUBSET_SCORES_H */
A beyesian infererence-based sampler.
void set_was_used(bool tf) const
void add_restraints(RMF::FileHandle fh, const kernel::Restraints &hs)
A smart pointer to a ref-counted Object that is a class memeber.
Definition: base/Pointer.h:147
Object used to hold a set of restraints.
kernel::RestraintsTemp get_restraints(const Subset &s, const ParticleStatesTable *pst, const DependencyGraph &dg, kernel::RestraintSet *rs)
Represent a subset of the particles being optimized.
Definition: Subset.h:33
A class to change and restore log state.
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
#define IMP_LOG_TERSE(expr)
Store a subset of a subset or assignment.
Definition: Slice.h:27
Abstract base class for all restraints.
Various general useful functions for IMP.
A restraint is a term in an IMP ScoringFunction.
#define IMP_OBJECT_METHODS(Name)
Define the basic things needed by any Object.
double get_score(kernel::Restraint *r, const Assignment &a) const
Get the score of a set or restraint.
Common base class for heavy weight IMP objects.
Definition: base/Object.h:106
A beyesian infererence-based sampler.
A beyesian infererence-based sampler.
Functions to get report statistics about the used attributes.
virtual ScoringFunction * create_scoring_function(double weight=1.0, double max=NO_MAX) const
Store a configuration of a subset.
Definition: Assignment.h:32
unsigned int get_number_of_entries() const
Logging and error reporting support.
A shared base class to help in debugging and things.
void load_particle_states(const Subset &s, const Assignment &ss, const ParticleStatesTable *pst)
A beyesian infererence-based sampler.
void add_restraint(RMF::FileHandle fh, kernel::Restraint *hs)
Slice get_slice(Subset outer, Subset inner, const Subsets &excluded)
Definition: Slice.h:66