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