interp_krige.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2017-2018, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 #ifndef O2SCL_INTERP_KRIGE_H
24 #define O2SCL_INTERP_KRIGE_H
25 
26 /** \file interp_krige.h
27  \brief One-dimensional interpolation by Kriging
28 */
29 
30 #include <iostream>
31 #include <string>
32 
33 #include <gsl/gsl_sf_erf.h>
34 
35 #include <boost/numeric/ublas/vector.hpp>
36 #include <boost/numeric/ublas/matrix.hpp>
37 #include <boost/numeric/ublas/operation.hpp>
38 #include <boost/numeric/ublas/vector_proxy.hpp>
39 
40 #include <o2scl/interp.h>
41 #include <o2scl/cholesky.h>
42 #include <o2scl/lu.h>
43 #include <o2scl/columnify.h>
44 #include <o2scl/vector.h>
45 #include <o2scl/vec_stats.h>
46 #include <o2scl/min_brent_gsl.h>
47 #include <o2scl/constants.h>
48 
49 #ifndef DOXYGEN_NO_O2NS
50 namespace o2scl {
51 #endif
52 
53  /** \brief Interpolation by Kriging with a user-specified
54  covariance function
55 
56  See also the \ref intp_section section of the \o2 User's guide.
57 
58  \note The function \ref set() stores a pointer to the covariance
59  function and its derivatives and integrals so they cannot go out
60  of scope before any of the interpolation functions are called.
61 
62  \note This class is experimental.
63  */
64  template<class vec_t, class vec2_t=vec_t,
65  class covar_func_t=std::function<double(double,double)>,
66  class covar_integ_t=std::function<double(double,double,double)> >
67  class interp_krige : public interp_base<vec_t,vec2_t> {
68 
69 #ifdef O2SCL_NEVER_DEFINED
70  }{
71 #endif
72 
73  public:
74 
77  typedef boost::numeric::ublas::matrix_column<ubmatrix> ubmatrix_column;
78 
79  protected:
80 
81  /** \brief Inverse covariance matrix times function vector
82  */
83  ubvector Kinvf;
84 
85  /** \brief Pointer to user-specified covariance function
86  */
87  covar_func_t *f;
88 
89  /** \brief Pointer to user-specified derivative
90  */
91  covar_func_t *df;
92 
93  /** \brief Pointer to user-specified second derivative
94  */
95  covar_func_t *df2;
96 
97  /** \brief Pointer to user-specified second derivative
98  */
99  covar_integ_t *intp;
100 
101  public:
102 
103  interp_krige() {
104  this->min_size=2;
105  matrix_mode=0;
106  }
107 
108  virtual ~interp_krige() {}
109 
110  /// Method for matrix inversion
112 
113  /// Initialize interpolation routine
114  virtual void set(size_t size, const vec_t &x, const vec2_t &y) {
115  O2SCL_ERR2("Function set(size_t,vec_t,vec_t) unimplemented ",
116  "in interp_krige.",o2scl::exc_eunimpl);
117  return;
118  }
119 
120  /** \brief Initialize interpolation routine, specifying derivatives
121  and integrals [not yet implemented]
122  */
123  virtual int set_covar_di_noise(size_t n_dim, const vec_t &x,
124  const vec_t &y, covar_func_t &fcovar,
125  covar_func_t &fderiv,
126  covar_func_t &fderiv2,
127  covar_func_t &finteg,
128  double noise_var) {
129  O2SCL_ERR("Function set_covar_di_noise not yet implemented.",
131  return 0;
132  }
133 
134  /// Initialize interpolation routine
135  virtual int set_covar_noise(size_t n_dim, const vec_t &x, const vec_t &y,
136  covar_func_t &fcovar, double noise_var,
137  bool err_on_fail=true) {
138 
139  if (n_dim<this->min_size) {
140  O2SCL_ERR((((std::string)"Vector size, ")+szttos(n_dim)+", is less"+
141  " than "+szttos(this->min_size)+" in interp_krige::"+
142  "set().").c_str(),exc_einval);
143  }
144 
145  // Store pointer to covariance function
146  f=&fcovar;
147 
148  // Construct the KXX matrix
149  ubmatrix KXX(n_dim,n_dim);
150  for(size_t irow=0;irow<n_dim;irow++) {
151  for(size_t icol=0;icol<n_dim;icol++) {
152  if (irow>icol) {
153  KXX(irow,icol)=KXX(icol,irow);
154  } else if (irow==icol) {
155  KXX(irow,icol)=fcovar(x[irow],x[icol])+noise_var;
156  } else {
157  KXX(irow,icol)=fcovar(x[irow],x[icol]);
158  }
159  }
160  }
161 
162  if (matrix_mode!=0) {
163 
164  // Construct the inverse of KXX
165  ubmatrix inv_KXX(n_dim,n_dim);
166  o2scl::permutation p(n_dim);
167  int signum;
168  o2scl_linalg::LU_decomp(n_dim,KXX,p,signum);
169  if (o2scl_linalg::diagonal_has_zero(n_dim,KXX)) {
170  if (err_on_fail) {
171  O2SCL_ERR2("Matrix singular (LU method) ",
172  "in interp_krige::set_covar_noise().",
174  }
175  return 1;
176  }
177  o2scl_linalg::LU_invert<ubmatrix,ubmatrix,ubmatrix_column>
178  (n_dim,KXX,p,inv_KXX);
179 
180  // Inverse covariance matrix times function vector
181  this->Kinvf.resize(n_dim);
182  boost::numeric::ublas::axpy_prod(inv_KXX,y,this->Kinvf,true);
183 
184  } else {
185 
186  // Construct the inverse of KXX
187  int cret=o2scl_linalg::cholesky_decomp(n_dim,KXX,false);
188  if (cret!=0) {
189  if (err_on_fail) {
190  O2SCL_ERR2("Matrix singular (Cholesky method) ",
191  "in interp_krige::set_covar_noise().",
193  }
194  return 2;
195  }
196  ubmatrix inv_KXX=KXX;
197  o2scl_linalg::cholesky_invert<ubmatrix>(n_dim,inv_KXX);
198 
199  // Inverse covariance matrix times function vector
200  Kinvf.resize(n_dim);
201  boost::numeric::ublas::axpy_prod(inv_KXX,y,Kinvf,true);
202 
203  }
204 
205  // Set parent data members
206  this->px=&x;
207  this->py=&y;
208  this->sz=n_dim;
209 
210  return 0;
211  }
212 
213  /// Initialize interpolation routine
214  virtual int set_covar(size_t n_dim, const vec_t &x, const vec_t &y,
215  covar_func_t &fcovar) {
216  return set_covar_noise(n_dim,x,y,fcovar,0.0);
217  }
218 
219  /// Give the value of the function \f$ y(x=x_0) \f$ .
220  virtual double eval(double x0) const {
221 
222  // Initialize to zero to prevent uninit'ed var. warnings
223  double ret=0.0;
224  for(size_t i=0;i<this->sz;i++) {
225  ret+=(*f)(x0,(*this->px)[i])*Kinvf[i];
226  }
227 
228  return ret;
229  }
230 
231  /// Give the value of the derivative \f$ y^{\prime}(x=x_0) \f$ .
232  virtual double deriv(double x0) const {
233  return 0.0;
234  }
235 
236  /** \brief Give the value of the second derivative
237  \f$ y^{\prime \prime}(x=x_0) \f$ (always zero)
238  */
239  virtual double deriv2(double x0) const {
240  return 0.0;
241  }
242 
243  /// Give the value of the integral \f$ \int_a^{b}y(x)~dx \f$ .
244  virtual double integ(double a, double b) const {
245  return 0.0;
246  }
247 
248  /// Return the type, \c "interp_linear".
249  virtual const char *type() const { return "interp_krige"; }
250 
251 #ifndef DOXYGEN_INTERNAL
252 
253  private:
254 
258  (const interp_krige<vec_t,vec2_t,covar_func_t,covar_integ_t>&);
259 
260 #endif
261 
262  };
263 
264 
265  /** \brief One-dimensional interpolation using an
266  optimized covariance function
267 
268  See also the \ref intp_section section of the \o2 User's guide.
269 
270  \note This class is experimental.
271  */
272  template<class vec_t, class vec2_t=vec_t>
273  class interp_krige_optim : public interp_krige<vec_t,vec2_t> {
274 
275  public:
276 
279  typedef boost::numeric::ublas::matrix_column<ubmatrix> ubmatrix_column;
280 
281  protected:
282 
283  std::function<double(double,double)> ff;
284 
285  /// The covariance function length scale
286  double len;
287 
288  /// The quality factor of the optimization
289  double qual;
290 
291  /// The covariance function
292  double covar(double x1, double x2) {
293  return exp(-(x1-x2)*(x1-x2)/len/len/2.0);
294  }
295 
296  /// The derivative of the covariance function
297  double deriv(double x1, double x2) {
298  return -exp(-(x1-x2)*(x1-x2)/len/len/2.0)/len/len*(x1-x2);
299  }
300 
301  /// The second derivative of the covariance function
302  double deriv2(double x1, double x2) {
303  return ((x1-x2)*(x1-x2)-len*len)*
304  exp(-(x1-x2)*(x1-x2)/len/len/2.0)/len/len/len/len;
305  }
306 
307  /// The integral of the covariance function
308  double integ(double x, double x1, double x2) {
309  exit(-1);
310  // This is wrong
311  //return 0.5*len*sqrt(o2scl_const::pi)*
312  //(gsl_sf_erf((x2-x)/len)+gsl_sf_erf((x-x1)/len));
313  }
314 
315  /// Pointer to the user-specified minimizer
317 
318  /** \brief Function to optimize the covariance parameters
319  */
320  double qual_fun(double x, double noise_var, int &success) {
321 
322  len=x;
323  success=0;
324 
325  size_t size=this->sz;
326 
327  // Leave-one-out cross validation
328  static const size_t loo_cv=1;
329  // Log-marginal-likelihood
330  static const size_t lml=2;
331  size_t mode=lml;
332 
333  if (mode==loo_cv) {
334 
335  qual=0.0;
336  for(size_t k=0;k<size;k++) {
337 
338  // Leave one observation out
339  ubvector x2(size-1);
340  o2scl::vector_copy_jackknife((*this->px),k,x2);
341  ubvector y2(size-1);
342  o2scl::vector_copy_jackknife((*this->py),k,y2);
343 
344  // Construct the KXX matrix
345  ubmatrix KXX(size-1,size-1);
346  for(size_t irow=0;irow<size-1;irow++) {
347  for(size_t icol=0;icol<size-1;icol++) {
348  if (irow>icol) {
349  KXX(irow,icol)=KXX(icol,irow);
350  } else {
351  KXX(irow,icol)=exp(-pow((x2[irow]-x2[icol])/len,2.0)/2.0);
352  if (irow==icol) KXX(irow,icol)+=noise_var;
353  }
354  }
355  }
356 
357  if (this->matrix_mode!=0) {
358 
359  // Construct the inverse of KXX
360  ubmatrix inv_KXX(size-1,size-1);
361  o2scl::permutation p(size-1);
362  int signum;
363  o2scl_linalg::LU_decomp(size-1,KXX,p,signum);
364  if (o2scl_linalg::diagonal_has_zero(size-1,KXX)) {
365  success=1;
366  return 1.0e99;
367  }
368  o2scl_linalg::LU_invert<ubmatrix,ubmatrix,ubmatrix_column>
369  (size-1,KXX,p,inv_KXX);
370 
371  // Inverse covariance matrix times function vector
372  this->Kinvf.resize(size-1);
373  boost::numeric::ublas::axpy_prod(inv_KXX,y2,this->Kinvf,true);
374 
375  } else {
376 
377  // Construct the inverse of KXX
378  int cret=o2scl_linalg::cholesky_decomp(size-1,KXX,false);
379  if (cret!=0) {
380  success=2;
381  return 1.0e99;
382  }
383  ubmatrix &inv_KXX=KXX;
384  o2scl_linalg::cholesky_invert<ubmatrix>(size-1,inv_KXX);
385 
386  // Inverse covariance matrix times function vector
387  this->Kinvf.resize(size-1);
388  boost::numeric::ublas::axpy_prod(inv_KXX,y2,this->Kinvf,true);
389 
390  }
391 
392  double ypred=0.0;
393  double yact=(*this->py)[k];
394  for(size_t i=0;i<size-1;i++) {
395  ypred+=exp(-pow(((*this->px)[k]-x2[i])/len,2.0)/2.0)*this->Kinvf[i];
396  }
397 
398  // Measure the quality with a chi-squared like function
399  qual+=pow(yact-ypred,2.0);
400 
401  }
402 
403  } else if (mode==lml) {
404 
405  // Construct the KXX matrix
406  ubmatrix KXX(size,size);
407  for(size_t irow=0;irow<size;irow++) {
408  for(size_t icol=0;icol<size;icol++) {
409  if (irow>icol) {
410  KXX(irow,icol)=KXX(icol,irow);
411  } else {
412  KXX(irow,icol)=exp(-pow(((*this->px)[irow]-
413  (*this->px)[icol])/len,2.0)/2.0);
414  if (irow==icol) KXX(irow,icol)+=noise_var;
415  }
416  }
417  }
418 
419  // Note: We have to use LU here because O2scl doesn't yet
420  // have a lndet() function for Cholesky decomp
421 
422  // Construct the inverse of KXX
423  ubmatrix inv_KXX(size,size);
424  o2scl::permutation p(size);
425  int signum;
426  o2scl_linalg::LU_decomp(size,KXX,p,signum);
427  if (o2scl_linalg::diagonal_has_zero(size,KXX)) {
428  success=1;
429  return 1.0e99;
430  }
431  o2scl_linalg::LU_invert<ubmatrix,ubmatrix,ubmatrix_column>
432  (size,KXX,p,inv_KXX);
433 
434  double lndet=o2scl_linalg::LU_lndet<ubmatrix>(size,KXX);
435 
436  // Inverse covariance matrix times function vector
437  this->Kinvf.resize(size);
438  boost::numeric::ublas::axpy_prod(inv_KXX,*this->py,this->Kinvf,true);
439 
440  // Compute the log of the marginal likelihood, without
441  // the constant term
442  for(size_t i=0;i<size;i++) {
443  qual+=-0.5*(*this->py)[i]*this->Kinvf[i];
444  }
445  qual-=0.5*lndet;
446  }
447 
448  return qual;
449  }
450 
451  public:
452 
453  /// Verbosity parameter
454  int verbose;
455 
456  /// Number of length scale points to try
457  size_t nlen;
458 
459  /// Default minimizer
461 
462  /// If true, use the full minimizer
463  bool full_min;
464 
466  nlen=20;
467  full_min=false;
468  mp=&def_min;
469  verbose=0;
470  }
471 
472  /// Initialize interpolation routine
473  virtual int set_noise(size_t size, const vec_t &x, const vec2_t &y,
474  double noise_var) {
475 
476  // Set parent data members
477  this->px=&x;
478  this->py=&y;
479  this->sz=size;
480 
481  if (full_min) {
482 
483  double p=x[1]-x[0];
484 
485  int success=0;
486  funct mf=std::bind
487  (std::mem_fn<double(double,double,int &)>
489  std::placeholders::_1,noise_var,std::ref(success));
490 
491  mp->min(p,qual,mf);
492 
493  } else {
494 
495  if (verbose>1) {
496  std::cout << "interp_krige_optim: simple minimization"
497  << std::endl;
498  }
499 
500  // Range of the length parameter
501  std::vector<double> diff(size-1);
502  for(size_t i=0;i<size-1;i++) {
503  diff[i]=fabs(x[i+1]-x[i]);
504  }
505  double len_min=o2scl::vector_min_value
506  <std::vector<double>,double>(size-1,diff)/3.0;
507  double len_max=fabs(x[size-1]-x[0])*3.0;
508  double len_ratio=len_max/len_min;
509 
510  if (verbose>1) {
511  std::cout << "len: " << len_min << " "
512  << len_max << " "
513  << pow(len_ratio,((double)1)/((double)nlen-1))
514  << std::endl;
515  }
516 
517  // Initialize to zero to prevent uninit'ed var. warnings
518  double min_qual=0.0, len_opt=0.0;
519 
520  if (verbose>1) {
521  std::cout << "ilen len qual fail min_qual" << std::endl;
522  }
523 
524  // Loop over the full range, finding the optimum
525  bool min_set=false;
526  for(size_t j=0;j<nlen;j++) {
527  len=len_min*pow(len_ratio,((double)j)/((double)nlen-1));
528 
529  int success=0;
530  qual=qual_fun(len,noise_var,success);
531 
532  if (success==0 && (min_set==false || qual<min_qual)) {
533  len_opt=len;
534  min_qual=qual;
535  min_set=true;
536  }
537 
538  if (verbose>1) {
539  std::cout << "interp_krige_optim: ";
540  std::cout.width(2);
541  std::cout << j << " " << len << " " << qual << " "
542  << success << " " << min_qual << " "
543  << len_opt << std::endl;
544  }
545 
546  }
547 
548  // Now that we've optimized the covariance function,
549  // just use the parent class to interpolate
550  len=len_opt;
551 
552  }
553 
554  ff=std::bind(std::mem_fn<double(double,double)>
556  std::placeholders::_1,std::placeholders::_2);
557 
558  this->set_covar_noise(size,x,y,ff,noise_var);
559 
560  return 0;
561  }
562 
563  /// Initialize interpolation routine
564  virtual void set(size_t size, const vec_t &x, const vec2_t &y) {
565  double mean_abs=0.0;
566  for(size_t j=0;j<size;j++) {
567  mean_abs+=fabs(y[j]);
568  }
569  mean_abs/=size;
570  set_noise(size,x,y,mean_abs/1.0e8);
571  return;
572  }
573 
574 
575  };
576 
577 #ifndef DOXYGEN_NO_O2NS
578 }
579 #endif
580 
581 #endif
double qual_fun(double x, double noise_var, int &success)
Function to optimize the covariance parameters.
Definition: interp_krige.h:320
data_t vector_min_value(size_t n, const vec_t &data)
Compute the minimum of the first n elements of a vector.
Definition: vector.h:1199
covar_func_t * df
Pointer to user-specified derivative.
Definition: interp_krige.h:91
void vector_copy_jackknife(const vec_t &src, size_t iout, vec2_t &dest)
From a given vector, create a new vector by removing a specified element.
Definition: vector.h:2367
virtual double integ(double a, double b) const
Give the value of the integral .
Definition: interp_krige.h:244
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
size_t min_size
The minimum size of the vectors to interpolate between.
Definition: interp.h:165
double integ(double x, double x1, double x2)
The integral of the covariance function.
Definition: interp_krige.h:308
virtual int set_covar_noise(size_t n_dim, const vec_t &x, const vec_t &y, covar_func_t &fcovar, double noise_var, bool err_on_fail=true)
Initialize interpolation routine.
Definition: interp_krige.h:135
covar_func_t * f
Pointer to user-specified covariance function.
Definition: interp_krige.h:87
min_base * mp
Pointer to the user-specified minimizer.
Definition: interp_krige.h:316
size_t sz
Vector size.
Definition: interp.h:131
int matrix_mode
Method for matrix inversion.
Definition: interp_krige.h:111
invalid argument supplied by user
Definition: err_hnd.h:59
const vec_t * px
Independent vector.
Definition: interp.h:125
apparent singularity detected
Definition: err_hnd.h:93
A class for representing permutations.
Definition: permutation.h:70
Base low-level interpolation class [abstract base].
Definition: interp.h:107
int LU_decomp(const size_t N, mat_t &A, o2scl::permutation &p, int &signum)
Compute the LU decomposition of the matrix A.
Definition: lu_base.h:86
Interpolation by Kriging with a user-specified covariance function.
Definition: interp_krige.h:67
double len
The covariance function length scale.
Definition: interp_krige.h:286
requested feature not (yet) implemented
Definition: err_hnd.h:99
One-dimensional interpolation using an optimized covariance function.
Definition: interp_krige.h:273
virtual double deriv(double x0) const
Give the value of the derivative .
Definition: interp_krige.h:232
double deriv2(double x1, double x2)
The second derivative of the covariance function.
Definition: interp_krige.h:302
covar_func_t * df2
Pointer to user-specified second derivative.
Definition: interp_krige.h:95
One-dimensional minimization using Brent&#39;s method (GSL)
Definition: min_brent_gsl.h:95
covar_integ_t * intp
Pointer to user-specified second derivative.
Definition: interp_krige.h:99
double deriv(double x1, double x2)
The derivative of the covariance function.
Definition: interp_krige.h:297
min_brent_gsl def_min
Default minimizer.
Definition: interp_krige.h:460
virtual int set_covar(size_t n_dim, const vec_t &x, const vec_t &y, covar_func_t &fcovar)
Initialize interpolation routine.
Definition: interp_krige.h:214
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
virtual int set_noise(size_t size, const vec_t &x, const vec2_t &y, double noise_var)
Initialize interpolation routine.
Definition: interp_krige.h:473
const vec2_t * py
Dependent vector.
Definition: interp.h:128
int cholesky_decomp(const size_t M, mat_t &A, bool err_on_fail=true)
Compute the in-place Cholesky decomposition of a symmetric positive-definite square matrix...
Definition: cholesky_base.h:62
virtual int min(double &x, double &fmin, func_t &func)=0
Calculate the minimum min of func w.r.t &#39;x&#39;.
virtual const char * type() const
Return the type, "interp_linear".
Definition: interp_krige.h:249
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
bool full_min
If true, use the full minimizer.
Definition: interp_krige.h:463
ubvector Kinvf
Inverse covariance matrix times function vector.
Definition: interp_krige.h:83
size_t nlen
Number of length scale points to try.
Definition: interp_krige.h:457
One-dimensional minimization [abstract base].
Definition: min.h:42
double covar(double x1, double x2)
The covariance function.
Definition: interp_krige.h:292
virtual double deriv2(double x0) const
Give the value of the second derivative (always zero)
Definition: interp_krige.h:239
std::function< double(double)> funct
One-dimensional function typedef.
Definition: funct.h:45
double qual
The quality factor of the optimization.
Definition: interp_krige.h:289
virtual int set_covar_di_noise(size_t n_dim, const vec_t &x, const vec_t &y, covar_func_t &fcovar, covar_func_t &fderiv, covar_func_t &fderiv2, covar_func_t &finteg, double noise_var)
Initialize interpolation routine, specifying derivatives and integrals [not yet implemented].
Definition: interp_krige.h:123
int diagonal_has_zero(const size_t N, mat_t &A)
Return 1 if at least one of the elements in the diagonal is zero.
Definition: lu_base.h:54
int verbose
Verbosity parameter.
Definition: interp_krige.h:454
static const double x2[5]
Definition: inte_qng_gsl.h:66
virtual double eval(double x0) const
Give the value of the function .
Definition: interp_krige.h:220
static const double x1[5]
Definition: inte_qng_gsl.h:48
std::string szttos(size_t x)
Convert a size_t to a string.
Success.
Definition: err_hnd.h:47

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).