IMP  2.1.0
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-2013 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 
23 #if IMP_DOMINO_HAS_RMF
24 #include <RMF/HDF5/Group.h>
25 #endif
26 
27 IMPDOMINO_BEGIN_NAMESPACE
28 
29 /** Implement a cache for restraint scores as well as management of restraints
30  for domino.
31 
32  The cache size passed to the constructor is the maximum number of scores
33  that will be saved. A least-recently-used eviction policy is used when
34  that number is exceeded.
35 */
36 class IMPDOMINOEXPORT RestraintCache : public base::Object {
37  IMP_NAMED_TUPLE_2(Key, Keys, base::WeakPointer<kernel::Restraint>, restraint,
38  Assignment, assignment, );
39  IMP_NAMED_TUPLE_3(RestraintData, estraintDatas,
40  base::PointerMember<ScoringFunction>, scoring_function,
41  Subset, subset, double, max, );
42  IMP_NAMED_TUPLE_2(RestraintSetData, RestraintSetDatas, Slice, slice,
44  IMP_NAMED_TUPLE_2(SetData, SetDatas, RestraintSetDatas, members, double,
45  max, );
46 
47 #ifndef IMP_DOXYGEN
48  class Generator {
50  RMap rmap_;
52  SMap sets_;
54 
55  public:
56  Generator(ParticleStatesTable *pst) : pst_(pst) {}
57  typedef double result_type;
58  typedef Key argument_type;
59  template <class Cache>
60  result_type operator()(const argument_type &k, const Cache &cache) const {
61  RMap::const_iterator it = rmap_.find(k.get_restraint());
62  if (it != rmap_.end()) {
63  Subset s = rmap_.find(k.get_restraint())->second.get_subset();
64  load_particle_states(s, k.get_assignment(), pst_);
65  double e;
66  {
68  e = it->second.get_scoring_function()->evaluate_if_below(
69  false, it->second.get_max());
70  }
71  IMP_LOG_TERSE("Restraint " << Showable(k.get_restraint())
72  << " evaluated to " << e << " on "
73  << k.get_assignment() << " vs "
74  << it->second.get_max() << std::endl);
75  // prob can go away with ScoreFunction change
76  if (e > it->second.get_max()) e = std::numeric_limits<double>::max();
77  return e;
78  } else {
79  SMap::const_iterator it = sets_.find(k.get_restraint());
80  IMP_USAGE_CHECK(it != sets_.end(), "Restraint set "
81  << Showable(k.get_restraint())
82  << " not found.");
83  double total = 0;
84  for (unsigned int i = 0; i < it->second.get_members().size(); ++i) {
85  Assignment cur = it->second.get_members()[i].get_slice().get_sliced(
86  k.get_assignment());
87  double score = cache.get(
88  argument_type(it->second.get_members()[i].get_restraint(), cur));
89  total += score * k.get_restraint()->get_weight();
90  if (total >= it->second.get_max()) {
91  break;
92  }
93  }
94  IMP_LOG_TERSE("Restraint " << Showable(k.get_restraint())
95  << " evaluated to " << total << " on "
96  << k.get_assignment() << " with max "
97  << it->second.get_max() << std::endl);
98  if (total >= it->second.get_max()) {
99  return std::numeric_limits<double>::max();
100  } else {
101  return total;
102  }
103  }
104  }
105  void add_to_set(kernel::RestraintSet *rs, kernel::Restraint *r, Slice slice,
106  double max) {
107  IMP_USAGE_CHECK(!dynamic_cast<kernel::RestraintSet *>(r),
108  "don't pass restraint sets here as second arg");
109  sets_[rs].access_members().push_back(RestraintSetData(slice, r));
110  sets_[rs].set_max(max);
111  }
112  void add_restraint(kernel::Restraint *e, Subset s, double max) {
113  IMP_USAGE_CHECK(!dynamic_cast<kernel::RestraintSet *>(e),
114  "don't pass restraint sets here");
115  if (rmap_.find(e) == rmap_.end()) {
116  rmap_[e] = RestraintData(e->create_scoring_function(1.0, max), s, max);
117  } else {
118  IMP_USAGE_CHECK(rmap_.find(e)->second.get_subset() == s,
119  "Subsets don't match on restraint update");
120  rmap_[e].set_max(std::min(rmap_[e].get_max(), max));
121  }
122  }
123  ParticleStatesTable *get_particle_states_table() const { return pst_; }
124  void show_restraint_information(std::ostream &out) const;
125  };
126 #endif // IMP_DOXYGEN
127 
128  struct ApproximatelyEqual {
129  bool operator()(double a, double b) const {
130  return std::abs(a - b) < .1 * (a + b) + .1;
131  }
132  };
134  void add_restraint_internal(kernel::Restraint *r, unsigned int index,
135  kernel::RestraintSet *parent, double max,
136  Subset parent_subset, const DepMap &dependencies);
137  void add_restraint_set_child_internal(kernel::Restraint *r,
138  const Subset &cur_subset,
139  kernel::RestraintSet *parent,
140  double parent_max,
141  Subset parent_subset);
142  void add_restraint_set_internal(kernel::RestraintSet *rs, unsigned int index,
143  const Subset &cur_subset, double cur_max,
144  const DepMap &dependencies);
145  Subset get_subset(kernel::Restraint *r, const DepMap &dependencies) const;
146  // otherwise doxygen seems to index this for some reason
147 #ifndef IMP_DOXYGEN
149  Cache cache_;
150 #endif
152  KnownRestraints known_restraints_;
153  // assign a unique index to each restraint for use with I/O
155  RestraintIndex restraint_index_;
156  unsigned int next_index_;
157 
158  public:
160  unsigned int size = std::numeric_limits<unsigned int>::max());
161  /** Recursively process the passed restraints (and sets) so all contained
162  restraints and sets that have maximum are known.*/
164  //! Get the score of a set or restraint
165  /** The returned score will be std::numeric_limits<double>::max()
166  if any of the limits are violated.
167 
168  The assignment passed is the assignment for the particles used by the
169  restraint, not that for some whole subset. See the other get_score()
170  function for a perhaps more useful one.
171  */
172  double get_score(kernel::Restraint *r, const Assignment &a) const {
173  set_was_used(true);
174  double s = cache_.get(Key(r, a));
175  return s;
176  }
177  /** The the score for a restraint given a subset and assignment on
178  that subset.
179  */
180  double get_score(kernel::Restraint *r, const Subset &s,
181  const Assignment &a) const;
182 
183  //! make it so kernel::Restraint::get_last_score() returns the score
184  /** This is useful when writing the restraints to disk, as that
185  code often goes off the last score to avoid recomputing the
186  restraints.*/
187  void load_last_score(kernel::Restraint *r, const Subset &s,
188  const Assignment &a);
189  /** Return the restraints that should be evaluated for the subset,
190  given the exclusions.*/
192  const Subsets &exclusions) const;
193 
195 
196 #if IMP_DOMINO_HAS_RMF || defined(IMP_DOXYGEN)
197  /** This assumes that restraints are always added to the cache
198  in the same order.
199  \param[in] particle_ordering An ordering for the particles.
200  \param[in] restraints Which restraints to write out entries for.
201  You probably want to use get_restraints() to generate this.
202  \param[in] max_entries How many entries to write out at most.
203  \param[in] group Where to put the entries.
204  */
205  void save_cache(const kernel::ParticlesTemp &particle_ordering,
206  const kernel::RestraintsTemp &restraints,
207  RMF::HDF5::Group group, unsigned int max_entries);
208  void load_cache(const kernel::ParticlesTemp &ps, RMF::HDF5::ConstGroup group);
209 #endif
210 
211  /** Return the slice for that restraint given the subset. */
212  Slice get_slice(kernel::Restraint *r, const Subset &s) const;
213 
214  /** Return the number of entries currently in the cache.*/
215  unsigned int get_number_of_entries() const { return cache_.size(); }
216 
217  /** Check the entries in the cache.*/
218  void validate() const;
219 
220  /** Print out information about the known restraints and restraint sets.*/
221  void show_restraint_information(std::ostream &out = std::cout) const;
222  double get_hit_rate() const { return cache_.get_hit_rate(); }
224 };
225 
226 IMPDOMINO_END_NAMESPACE
227 
228 #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:146
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.
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