IMP  2.0.0
The Integrative Modeling Platform
GaussianProcessInterpolationSparse.h
Go to the documentation of this file.
1 /**
2  * \file IMP/isd/GaussianProcessInterpolationSparse.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_SPARSE_H
9 #define IMPISD_GAUSSIAN_PROCESS_INTERPOLATION_SPARSE_H
10 
11 #include <IMP/isd/isd_config.h>
12 
13 #ifdef IMP_ISD_USE_CHOLMOD
14 
15 #include <IMP/macros.h>
16 #include <boost/scoped_ptr.hpp>
19 #include <Eigen/Dense>
20 #define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
21 #include <Eigen/Sparse>
22 #include <unsupported/Eigen/CholmodSupport>
23 #include <ufsparse/cholmod.h>
24 
25 
26 IMPISD_BEGIN_NAMESPACE
27 #ifndef SWIG
28 using Eigen::SparseMatrix;
29 using Eigen::MatrixXd;
30 using Eigen::VectorXd;
31 #endif
32 
33 class GaussianProcessInterpolationRestraintSparse;
34 
35 //! GaussianProcessInterpolationSparse
36 /** This class provides methods to perform bayesian interpolation/smoothing of
37  * data using a gaussian process prior on the function to be interpolated.
38  * It takes a dataset as input (via its sufficient statistics) along with prior
39  * mean and covariance functions. It outputs the value of the posterior mean and
40  * covariance functions at points requested by the user.
41  */
42 class IMPISDEXPORT GaussianProcessInterpolationSparse : public base::Object
43 {
44  public:
45  /** Constructor for the gaussian process
46  * \param(in) x : a list of coordinates in N-dimensional space
47  * corresponding to the abscissa of each observation
48  * \param(in) sample_mean \f$I\f$ : vector of mean observations at each of
49  * the previously defined coordinates
50  * \param(in) sample_std \f$s\f$ : vector of sample standard deviations
51  * \param(in) mean_function \f$m\f$ : a pointer to the prior mean function
52  * to use. Should be compatible with
53  * the size of x(i).
54  * \param(in) covariance_function \f$w\f$: prior covariance function.
55  * \param(in) cutoff : when to consider that W matrix entries are zero.
56  *
57  * Computes the necessary matrices and inverses when called.
58  */
59  GaussianProcessInterpolationSparse(FloatsList x,
60  Floats sample_mean,
61  Floats sample_std,
62  Ints n_obs,
63  UnivariateFunction *mean_function,
64  BivariateFunction *covariance_function,
65  double cutoff=1e-7);
66 
67  /** Get posterior mean and covariance functions, at the points requested
68  * Posterior mean is defined as
69  * \f[\hat{I}(x) = m(x)
70  * + {}^t\mathbf{w}(q)
71  * (\mathbf{W}+\mathbf{S})^{-1}
72  * (\mathbf{I}-\mathbf{m}) \f]
73  * Posterior covariance
74  * \f[\hat{\sigma}^2(x,x') =
75  * w(x,x') - {}^t\mathbf{w}(x)
76  * (\mathbf{W} + \mathbf{S})^{-1}
77  * \mathbf{w}(x') \f]
78  * where \f$\mathbf{m}\f$ is the vector built by evaluating the prior mean
79  * function at the observation points; \f$\mathbf{w}(x)\f$ is the vector of
80  * covariances between each observation point and the current point;
81  * \f$\mathbf{W}\f$ is the prior covariance matrix built by evaluating the
82  * covariance function at each of the observations; \f$\mathbf{S}\f$ is the
83  * diagonal covariance matrix built from sample_std and n_obs.
84  *
85  * Both functions will check if the mean or covariance functions have changed
86  * since the last call, and will recompute
87  * \f$(\mathbf{W} + \mathbf{S})^{-1}\f$ if necessary.
88  */
89  double get_posterior_mean(Floats x);
90  double get_posterior_covariance(Floats x1,
91  Floats x2);
92 
93  // call these if you called update() on the mean or covariance function.
94  // it will force update any internal variables dependent on these functions.
95  void force_mean_update();
96  void force_covariance_update();
97 
98  friend class GaussianProcessInterpolationRestraintSparse;
99 
100  IMP_OBJECT_INLINE(GaussianProcessInterpolationSparse,
101  out << "GaussianProcessInterpolationSparse :"
102  "learning from " << M_ << " "
103  << N_ << "-dimensional observations" << std::endl,
104  {
105  cholmod_free_factor(&L_, c_);
106  cholmod_free_dense(&WSIm_, c_);
107  cholmod_free_sparse(&WS_,c_);
108  cholmod_finish(c_);
109  });
110 
111  protected:
112  //returns updated data vector
113  VectorXd get_I() const {return I_;}
114  //returns updated prior mean vector
115  VectorXd get_m();
116  // returns updated prior covariance vector
117  cholmod_sparse *get_wx_vector(Floats xval);
118  //returns updated data covariance matrix
119  SparseMatrix<double> get_S() const {return S_;}
120  //returns updated prior covariance matrix
121  SparseMatrix<double> get_W();
122  //returns updated (W+S)^{-1}
123  cholmod_sparse *get_WS();
124  //returns updated cholesky factor
125  cholmod_factor *get_L();
126  //returns updated (W+S)^{-1}(I-m)
127  cholmod_dense *get_WSIm();
128 
129  private:
130 
131  // ensures the mean/covariance function has updated parameters. Signals an
132  // update by changing the state flags. Returns true if the function has
133  // changed. This is used by GaussianProcessInterpolationSparseRestraint.
134  void update_flags_mean();
135  void update_flags_covariance();
136 
137  // compute prior covariance matrix
138  void compute_W();
139  // compute \f$(\mathbf{W} + \mathbf{S})^{-1}\f$ (if necessary).
140  void compute_WS();
141  // compute (W+S)^{-1} (I-m)
142  void compute_WSIm();
143 
144  // compute mean observations
145  void compute_I(Floats mean);
146  // compute diagonal covariance matrix of observations
147  void compute_S(Floats std, Ints n);
148  // compute prior mean vector
149  void compute_m();
150 
151  // for GPIR.
152  cholmod_common *get_cholmod_common() { return c_; }
153 
154  private:
155  unsigned N_; // number of dimensions of the abscissa
156  unsigned M_; // number of observations to learn from
157  FloatsList x_; // abscissa
158  Ints n_obs_; // number of observations
159  // pointer to the prior mean function
160  IMP::internal::OwnerPointer<UnivariateFunction> mean_function_;
161  // pointer to the prior covariance function
162  IMP::internal::OwnerPointer<BivariateFunction> covariance_function_;
163  VectorXd I_,m_;
164  SparseMatrix<double> S_,W_,wx_;
165  cholmod_sparse cwx_; //cholmod view of eigen object
166  cholmod_sparse *WS_; // WS = (W + S)^{-1}
167  cholmod_dense *WSIm_; // WS * (I - m)
168  cholmod_factor *L_;
169  bool flag_m_, flag_m_gpir_, flag_WS_, flag_WSIm_, flag_W_, flag_W_gpir_;
170  double cutoff_; // sparseness cutoff for W.
171  cholmod_common Common_, *c_;
172 
173 };
174 
175 IMPISD_END_NAMESPACE
176 
177 #endif /* IMP_ISD_USE_CHOLMOD */
178 
179 #endif /* IMPISD_GAUSSIAN_PROCESS_INTERPOLATION_SPARSE_H */