IMP  2.3.1
The Integrative Modeling Platform
GaussianProcessInterpolation.h
Go to the documentation of this file.
1 /**
2  * \file IMP/isd/GaussianProcessInterpolation.h
3  * \brief Normal distribution of Function
4  *
5  * Copyright 2007-2014 IMP Inventors. All rights reserved.
6  */
7 
8 #ifndef IMPISD_GAUSSIAN_PROCESS_INTERPOLATION_H
9 #define IMPISD_GAUSSIAN_PROCESS_INTERPOLATION_H
10 
11 #include <IMP/isd/isd_config.h>
12 #include <IMP/macros.h>
13 #include <boost/scoped_ptr.hpp>
16 #include <IMP/isd/Scale.h>
17 #include <IMP/algebra/eigen3/Eigen/Dense>
18 #include <IMP/algebra/eigen3/Eigen/Cholesky>
19 
20 IMPISD_BEGIN_NAMESPACE
21 
22 class GaussianProcessInterpolationRestraint;
23 class GaussianProcessInterpolationScoreState;
24 
25 //! GaussianProcessInterpolation
26 /** This class provides methods to perform bayesian interpolation/smoothing of
27  * data using a gaussian process prior on the function to be interpolated.
28  * It takes a dataset as input (via its sufficient statistics) along with prior
29  * mean and covariance functions. It outputs the value of the posterior mean and
30  * covariance functions at points requested by the user.
31  */
32 class IMPISDEXPORT GaussianProcessInterpolation : public base::Object {
33  public:
34  /** Constructor for the gaussian process
35  * \param [in] x : a list of coordinates in N-dimensional space
36  * corresponding to the abscissa of each observation
37  * \param [in] sample_mean \f$I\f$ : vector of mean observations at each of
38  * the previously defined coordinates
39  * \param [in] sample_std \f$s\f$ : vector of sample standard deviations
40  * \param [in] n_obs \f$N\f$ : sample size for mean and std
41  * \param [in] mean_function \f$m\f$ : a pointer to the prior mean function
42  * to use. Should be compatible with
43  * the size of x(i).
44  * \param [in] covariance_function \f$w\f$: prior covariance function.
45  * \param [in] sigma : ISD Scale (proportionality factor to S)
46  * \param [in] sparse_cutoff : when to consider that a matrix entry is zero
47  *
48  * Computes the necessary matrices and inverses when called.
49  */
51  Floats sample_std, unsigned n_obs,
52  UnivariateFunction *mean_function,
53  BivariateFunction *covariance_function,
54  kernel::Particle *sigma,
55  double sparse_cutoff = 1e-7);
56 
57  /** Get posterior mean and covariance functions, at the points requested
58  * Posterior mean is defined as
59  * \f[\hat{I}(x) = m(x)
60  * + {}^t\mathbf{w}(q)
61  * (\mathbf{W}+\mathbf{S})^{-1}
62  * (\mathbf{I}-\mathbf{m}) \f]
63  * Posterior covariance
64  * \f[\hat{\sigma}^2(x,x') =
65  * w(x,x') - {}^t\mathbf{w}(x)
66  * (\mathbf{W} + \mathbf{S})^{-1}
67  * \mathbf{w}(x') \f]
68  * where \f$\mathbf{m}\f$ is the vector built by evaluating the prior mean
69  * function at the observation points; \f$\mathbf{w}(x)\f$ is the vector of
70  * covariances between each observation point and the current point;
71  * \f$\mathbf{W}\f$ is the prior covariance matrix built by evaluating the
72  * covariance function at each of the observations; \f$\mathbf{S}\f$ is the
73  * diagonal covariance matrix built from sample_std and n_obs.
74  *
75  * Both functions will check if the mean or covariance functions have
76  *changed
77  * since the last call, and will recompute
78  * \f$(\mathbf{W} + \mathbf{S})^{-1}\f$ if necessary.
79  */
80  double get_posterior_mean(Floats x) const;
81  double get_posterior_covariance(Floats x1, Floats x2) const;
82 #ifndef SWIG
83  // c++ only
84  IMP_Eigen::MatrixXd get_posterior_covariance_matrix(FloatsList x) const;
85 #endif
86  // for Python
87  FloatsList get_posterior_covariance_matrix(FloatsList x, bool) const;
88 
89 #ifndef SWIG
90 // derivative: d(mean(q))/(dparticle_i)
91 // IMP_Eigen::VectorXd get_posterior_mean_derivative(Floats x) const;
92 #endif
93 
94 #ifndef SWIG
95  // derivative: d(cov(q,q))/(dparticle_i)
96  IMP_Eigen::VectorXd get_posterior_covariance_derivative(Floats x) const;
97 #endif
98  // for Python
99  Floats get_posterior_covariance_derivative(Floats x, bool) const;
100 
101 #ifndef SWIG
102 // hessian: d(mean(q))/(dparticle_i dparticle_j)
103 // IMP_Eigen::MatrixXd get_posterior_mean_hessian(Floats x) const;
104 #endif
105 
106 #ifndef SWIG
107  // hessian: d(cov(q,q))/(dparticle_i dparticle_j)
108  IMP_Eigen::MatrixXd get_posterior_covariance_hessian(Floats x) const;
109 #endif
110  // for Python
111  FloatsList get_posterior_covariance_hessian(Floats x, bool) const;
112 
113  // needed for restraints using gpi
116 
117  // call these if you called update() on the mean or covariance function.
118  // it will force update any internal variables dependent on these functions.
119  void force_mean_update();
120  void force_covariance_update();
121 
122  // returns the number of particles that control m's values.
123  // public for testing purposes
124  unsigned get_number_of_m_particles() const;
125 
126  // returns true if a particle is optimized
127  // public for testing purposes
128  bool get_m_particle_is_optimized(unsigned i) const;
129 
130  // returns the number of particles that control Omega's values.
131  // public for testing purposes
132  unsigned get_number_of_Omega_particles() const;
133 
134  // returns true if a particle is optimized
135  // public for testing purposes
136  bool get_Omega_particle_is_optimized(unsigned i) const;
137 
138  // returns data
139  FloatsList get_data_abscissa() const;
140  Floats get_data_mean() const;
141  FloatsList get_data_variance() const;
142 
144  friend class GaussianProcessInterpolationScoreState;
145 
147 
148  protected:
149  // returns updated data vector
150  IMP_Eigen::VectorXd get_I() const { return I_; }
151  // returns updated prior mean vector
152  IMP_Eigen::VectorXd get_m() const;
153  // returns dm/dparticle
154  IMP_Eigen::VectorXd get_m_derivative(unsigned particle) const;
155  // returns d2m/(dparticle_1 dparticle_2)
156  IMP_Eigen::VectorXd get_m_second_derivative(unsigned particle1,
157  unsigned particle2) const;
158  // returns updated prior covariance vector
159  void add_to_m_particle_derivative(unsigned particle, double value,
160  DerivativeAccumulator &accum);
161  // returns updated prior covariance vector
162  IMP_Eigen::VectorXd get_wx_vector(Floats xval) const;
163  // returns updated data covariance matrix
164  IMP_Eigen::DiagonalMatrix<double, IMP_Eigen::Dynamic> get_S() const {
165  return S_;
166  }
167  // returns updated prior covariance matrix
168  IMP_Eigen::MatrixXd get_W() const;
169  // returns Omega=(W+S/N)
170  IMP_Eigen::MatrixXd get_Omega() const;
171  // returns dOmega/dparticle
172  IMP_Eigen::MatrixXd get_Omega_derivative(unsigned particle) const;
173  // returns d2Omega/(dparticle_1 dparticle_2)
174  IMP_Eigen::MatrixXd get_Omega_second_derivative(unsigned particle1,
175  unsigned particle2) const;
176  // returns updated prior covariance vector
177  void add_to_Omega_particle_derivative(unsigned particle, double value,
178  DerivativeAccumulator &accum);
179  // returns LDLT decomp of omega
180  IMP_Eigen::LDLT<IMP_Eigen::MatrixXd, IMP_Eigen::Upper> get_ldlt() const;
181  // returns updated Omega^{-1}
182  IMP_Eigen::MatrixXd get_Omi() const;
183  // returns updated Omega^{-1}(I-m)
184  IMP_Eigen::VectorXd get_OmiIm() const;
185 
186  private:
187  // ensures the mean/covariance function has updated parameters. Signals an
188  // update by changing the state flags. Returns true if the function has
189  // changed. This is used by GaussianProcessInterpolationRestraint.
190  void update_flags_mean();
191  void update_flags_covariance();
192 
193  // compute prior covariance matrix
194  void compute_W();
195  // compute \f$(\mathbf{W} + \frac{\sigma}{N}\mathbf{S})^{-1}\f$.
196  void compute_Omega();
197  // compute LDLT decomposition of Omega
198  void compute_ldlt();
199  // compute \f$(\mathbf{W} + \frac{\sigma}{N}\mathbf{S})^{-1}\f$.
200  void compute_Omi();
201  // compute (W+sigma*S/N)^{-1} (I-m)
202  void compute_OmiIm();
203 
204  // compute dw(q)/dparticle_i
205  IMP_Eigen::VectorXd get_wx_vector_derivative(Floats q, unsigned i) const;
206  // compute dw(q)/(dparticle_i * dparticle_j)
207  IMP_Eigen::VectorXd get_wx_vector_second_derivative(Floats q, unsigned i,
208  unsigned j) const;
209 
210  // compute dcov(q,q)/dw(q)
211  IMP_Eigen::VectorXd get_dcov_dwq(Floats q) const;
212  // compute dcov(q,q)/dOmega
213  IMP_Eigen::MatrixXd get_dcov_dOm(Floats q) const;
214  // compute d2cov(q,q)/(dw(q) dw(q)) (independent of q !)
215  IMP_Eigen::MatrixXd get_d2cov_dwq_dwq() const;
216  // compute d2cov(q,q)/(dw(q)_m dOmega)
217  IMP_Eigen::MatrixXd get_d2cov_dwq_dOm(Floats q, unsigned m) const;
218  // compute d2cov(q,q)/(dOmega dOmega_mn)
219  IMP_Eigen::MatrixXd get_d2cov_dOm_dOm(Floats q, unsigned m, unsigned n) const;
220 
221  // compute mean observations
222  void compute_I(Floats mean);
223  // compute diagonal covariance matrix of observations
224  void compute_S(Floats std);
225  // compute prior mean vector
226  void compute_m();
227 
228  private:
229  unsigned N_; // number of dimensions of the abscissa
230  unsigned M_; // number of observations to learn from
231  FloatsList x_; // abscissa
232  unsigned n_obs_; // number of observations
233  // pointer to the prior mean function
235  // pointer to the prior covariance function
236  IMP::base::PointerMember<BivariateFunction> covariance_function_;
237  IMP_Eigen::VectorXd I_, m_;
238  IMP_Eigen::MatrixXd W_, Omega_, Omi_; // Omi = Omega^{-1}
239  IMP_Eigen::DiagonalMatrix<double, IMP_Eigen::Dynamic> S_;
240  IMP_Eigen::VectorXd OmiIm_; // Omi * (I - m)
241  bool flag_m_, flag_m_gpir_, flag_Omi_, flag_OmiIm_, flag_W_, flag_Omega_,
242  flag_Omega_gpir_, flag_ldlt_;
244  double cutoff_;
245  double sigma_val_; // to know if an update is needed
246  IMP_Eigen::LDLT<IMP_Eigen::MatrixXd, IMP_Eigen::Upper> ldlt_;
247 };
248 
249 IMPISD_END_NAMESPACE
250 
251 #endif /* IMPISD_GAUSSIAN_PROCESS_INTERPOLATION_H */
Base class for functions of one variable.
Class for adding derivatives from restraints to the model.
A smart pointer to a ref-counted Object that is a class member.
Definition: Pointer.h:147
#define IMP_OBJECT_METHODS(Name)
Define the basic things needed by any Object.
Definition: object_macros.h:25
A decorator for scale parameters particles.
ParticlesTemp get_input_particles(const ModelObjectsTemp &mos)
Return all the input particles for a given ModelObject.
Classes for general functions.
Import IMP/kernel/macros.h in the namespace.
ContainersTemp get_input_containers(const ModelObjectsTemp &mos)
Return all the input particles for a given ModelObject.
IMP::base::Vector< IMP::base::WeakPointer< Container > > ContainersTemp
Base class for functions of two variables.
Class to handle individual model particles.
Common base class for heavy weight IMP objects.
Definition: Object.h:106
Classes for general functions.