IMP  2.0.0
The Integrative Modeling Platform
GaussianProcessInterpolationRestraint.h
Go to the documentation of this file.
1 /**
2  * \file IMP/isd/GaussianProcessInterpolationRestraint.h
3  * \brief Restraint and ScoreState for GaussianProcessInterpolation
4  *
5  * Copyright 2007-2013 IMP Inventors. All rights reserved.
6  */
7 
8 #ifndef IMPISD_GAUSSIAN_PROCESS_INTERPOLATION_RESTRAINT_H
9 #define IMPISD_GAUSSIAN_PROCESS_INTERPOLATION_RESTRAINT_H
10 
11 #include <IMP/isd/isd_config.h>
12 #include <IMP/macros.h>
13 #include <boost/scoped_ptr.hpp>
14 #include <IMP/isd/ISDRestraint.h>
17 #include <IMP/internal/OwnerPointer.h>
18 #include <Eigen/Dense>
19 #include <Eigen/Cholesky>
20 
21 #include <IMP/ScoreState.h>
22 
23 
24 IMPISD_BEGIN_NAMESPACE
25 #ifndef SWIG
26 using Eigen::MatrixXd;
27 using Eigen::VectorXd;
28 using Eigen::RowVectorXd;
29 #endif
30 
31 class GaussianProcessInterpolationScoreState;
32 
33 //! gaussian process restraint
34 /* the restraint is a multivariate normal distribution on the vector of
35 * observations with mean and standard deviation given by the posterior of the
36 * gaussian process.
37 */
39 {
40  private:
41  // checks and makes necessary updates
42  void update_mean_and_covariance();
43 
44  private:
45  IMP::Pointer<GaussianProcessInterpolation> gpi_;
46  IMP::internal::OwnerPointer<MultivariateFNormalSufficient> mvn_;
47  IMP::internal::OwnerPointer<GaussianProcessInterpolationScoreState> ss_;
48  //number of observation points
49  unsigned M_;
50 
51  public:
52  // this is a restraint on other restraints. It first constructs the
53  // necessary vectors from GaussianProcessInterpolation, then creates a
54  // multivariate normal distribution around it. Upon evaluation, it
55  // checks if parameters have changed, reconstructs the matrix if
56  // necessary, changes the DA weight and passes it to the functions.
59 
60  //to call this, you need to update the scorestate before.
61  //calling model.evaluate(False) is enough.
62  double get_probability() const
63  {
64  return mvn_->density();
65  }
66 
67  void stats() const
68  {
69  return mvn_->stats();
70  }
71 
72  //use conjugate gradients when possible (default false)
73  void set_use_cg(bool use, double tol) {mvn_->set_use_cg(use,tol);}
74 
75  //get minus log normalization and minus exponent separately
76  double get_minus_log_normalization() const;
77  double get_minus_exponent() const;
78 
79  //get hessian of the minus log likelihood
80  Eigen::MatrixXd get_hessian() const;
81 
82  //get log determinant of hessian
83  double get_logdet_hessian() const;
84 
85  //call this one from python
86  FloatsList get_hessian(bool unused) const;
87 
88  public:
89  double unprotected_evaluate(IMP::DerivativeAccumulator *accum)
90  const IMP_OVERRIDE;
91  IMP::ModelObjectsTemp do_get_inputs() const IMP_OVERRIDE;
93 
94  //needed to register the score state
95  void set_model(Model *m);
96 
97  //to allow the scorestate to get the restraint's objects
98  friend class GaussianProcessInterpolationScoreState;
99 };
100 
101 #if !defined(IMP_DOXYGEN) && !defined(SWIG)
102 class IMPISDEXPORT GaussianProcessInterpolationScoreState : public ScoreState
103 {
104  private:
105  IMP::WeakPointer<GaussianProcessInterpolationRestraint> gpir_;
106 
107  private:
108  GaussianProcessInterpolationScoreState(
109  GaussianProcessInterpolationRestraint *gpir) : gpir_(gpir) {}
110 
111  public:
112  //only the GPIR can create this and add it to the model
113  friend class GaussianProcessInterpolationRestraint;
114  IMP_SCORE_STATE(GaussianProcessInterpolationScoreState);
115 };
116 #endif
117 
118 IMPISD_END_NAMESPACE
119 
120 #endif /* IMPISD_GAUSSIAN_PROCESS_INTERPOLATION_RESTRAINT_H */