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