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/lu.h>
42 #include <o2scl/columnify.h>
43 #include <o2scl/vector.h>
44 #include <o2scl/vec_stats.h>
45 #include <o2scl/mmin_bfgs2.h>
46 #include <o2scl/constants.h>
47 
48 #ifndef DOXYGEN_NO_O2NS
49 namespace o2scl {
50 #endif
51 
52  /** \brief Interpolation by Kriging with a user-specified
53  covariance function
54 
55  See also the \ref intp_section section of the \o2 User's guide.
56 
57  \note The function \ref set() stores a pointer to the covariance
58  function so it cannot go out of scope before any of the
59  interpolation functions are called.
60 
61  \note This class is experimental.
62  */
63  template<class vec_t, class vec2_t=vec_t,
64  class covar_func_t=std::function<double(double,double)>,
65  class covar_integ_t=std::function<double(double,double,double)> >
66  class interp_krige : public interp_base<vec_t,vec2_t> {
67 
68 #ifdef O2SCL_NEVER_DEFINED
69  }{
70 #endif
71 
72  public:
73 
76  typedef boost::numeric::ublas::matrix_column<ubmatrix> ubmatrix_column;
77 
78  protected:
79 
80  /** \brief Inverse covariance matrix times function vector
81  */
82  ubvector Kinvf;
83 
84  /** \brief Pointer to user-specified covariance function
85  */
86  covar_func_t *f;
87 
88  /** \brief Pointer to user-specified derivative
89  */
90  covar_func_t *df;
91 
92  /** \brief Pointer to user-specified second derivative
93  */
94  covar_func_t *df2;
95 
96  /** \brief Pointer to user-specified second derivative
97  */
98  covar_integ_t *intp;
99 
100  public:
101 
102  interp_krige() {
103  this->min_size=2;
104  }
105 
106  virtual ~interp_krige() {}
107 
108  /// Initialize interpolation routine
109  virtual void set(size_t size, const vec_t &x, const vec2_t &y) {
110  O2SCL_ERR2("Function set(size_t,vec_t,vec_t) unimplemented ",
111  "in interp_krige.",o2scl::exc_eunimpl);
112  return;
113  }
114 
115  /** \brief Initialize interpolation routine, specifying derivatives
116  and integrals
117  */
118  virtual void set_covar_di_noise(size_t n_dim, const vec_t &x,
119  const vec_t &y, covar_func_t &fcovar,
120  covar_func_t &fderiv,
121  covar_func_t &fderiv2,
122  covar_func_t &finteg,
123  double noise_var) {
124  return;
125  }
126 
127  /// Initialize interpolation routine
128  virtual void set_covar_noise(size_t n_dim, const vec_t &x, const vec_t &y,
129  covar_func_t &fcovar, double noise_var) {
130 
131  if (n_dim<this->min_size) {
132  O2SCL_ERR((((std::string)"Vector size, ")+szttos(n_dim)+", is less"+
133  " than "+szttos(this->min_size)+" in interp_krige::"+
134  "set().").c_str(),exc_einval);
135  }
136 
137  // Store pointer to covariance function
138  f=&fcovar;
139 
140  // Construct the KXX matrix
141  ubmatrix KXX(n_dim,n_dim);
142  for(size_t irow=0;irow<n_dim;irow++) {
143  for(size_t icol=0;icol<n_dim;icol++) {
144  if (irow>icol) {
145  KXX(irow,icol)=KXX(icol,irow);
146  } else if (irow==icol) {
147  KXX(irow,icol)=fcovar(x[irow],x[icol])+noise_var;
148  } else {
149  KXX(irow,icol)=fcovar(x[irow],x[icol]);
150  }
151  }
152  }
153 
154  // Construct the inverse of KXX
155  ubmatrix inv_KXX(n_dim,n_dim);
156  o2scl::permutation p(n_dim);
157  int signum;
158  o2scl_linalg::LU_decomp(n_dim,KXX,p,signum);
159  if (o2scl_linalg::diagonal_has_zero(n_dim,KXX)) {
160  O2SCL_ERR("KXX matrix is singular in interp_krige::set().",
162  }
163  o2scl_linalg::LU_invert<ubmatrix,ubmatrix,ubmatrix_column>
164  (n_dim,KXX,p,inv_KXX);
165 
166  // Inverse covariance matrix times function vector
167  Kinvf.resize(n_dim);
168  boost::numeric::ublas::axpy_prod(inv_KXX,y,Kinvf,true);
169 
170  // Set parent data members
171  this->px=&x;
172  this->py=&y;
173  this->sz=n_dim;
174  return;
175  }
176 
177  /// Initialize interpolation routine
178  virtual void set_covar(size_t n_dim, const vec_t &x, const vec_t &y,
179  covar_func_t &fcovar) {
180  set_covar_noise(n_dim,x,y,fcovar,0.0);
181  return;
182  }
183 
184  /// Give the value of the function \f$ y(x=x_0) \f$ .
185  virtual double eval(double x0) const {
186 
187  // Initialize to zero to prevent uninit'ed var. warnings
188  double ret=0.0;
189  for(size_t i=0;i<this->sz;i++) {
190  ret+=(*f)(x0,(*this->px)[i])*Kinvf[i];
191  }
192 
193  return ret;
194  }
195 
196  /// Give the value of the derivative \f$ y^{\prime}(x=x_0) \f$ .
197  virtual double deriv(double x0) const {
198  return 0.0;
199  }
200 
201  /** \brief Give the value of the second derivative
202  \f$ y^{\prime \prime}(x=x_0) \f$ (always zero)
203  */
204  virtual double deriv2(double x0) const {
205  return 0.0;
206  }
207 
208  /// Give the value of the integral \f$ \int_a^{b}y(x)~dx \f$ .
209  virtual double integ(double a, double b) const {
210  return 0.0;
211  }
212 
213  /// Return the type, \c "interp_linear".
214  virtual const char *type() const { return "interp_krige"; }
215 
216 #ifndef DOXYGEN_INTERNAL
217 
218  private:
219 
223  (const interp_krige<vec_t,vec2_t,covar_func_t,covar_integ_t>&);
224 
225 #endif
226 
227  };
228 
229 
230  /** \brief One-dimensional interpolation using an
231  optimized covariance function
232 
233  See also the \ref intp_section section of the \o2 User's guide.
234 
235  \note This class is experimental.
236  */
237  template<class vec_t, class vec2_t=vec_t>
238  class interp_krige_optim : public interp_krige<vec_t,vec2_t> {
239 
240  public:
241 
244  typedef boost::numeric::ublas::matrix_column<ubmatrix> ubmatrix_column;
245 
246  protected:
247 
248  /// The covariance function length scale
249  double len;
250 
251  /// The covariance function coefficient
252  double var;
253 
254  /// The quality factor of the optimization
255  double qual;
256 
257  /// The covariance function
258  double covar(double x1, double x2) {
259  return var*exp(-(x1-x2)*(x1-x2)/len/len);
260  }
261 
262  /// The derivative of the covariance function
263  double deriv(double x1, double x2) {
264  return -2.0*var*exp(-(x1-x2)*(x1-x2)/len/len)/len/len*(x1-x2);
265  }
266 
267  /// The second derivative of the covariance function
268  double deriv2(double x1, double x2) {
269  return (4.0*(x1-x2)*(x1-x2)-2.0*len*len)*
270  var*exp(-(x1-x2)*(x1-x2)/len/len)/len/len/len/len;
271  }
272 
273  /// The integral of the covariance function
274  double integ(double x, double x1, double x2) {
275  return 0.5*len*sqrt(o2scl_const::pi)*var*
276  (gsl_sf_erf((x2-x)/len)+gsl_sf_erf((x-x1)/len));
277  }
278 
279  /// Functor for the covariance function \ref covar()
280  std::function<double(double,double)> ff;
281 
282  /// Pointer to the user-specified minimizer
284 
285  /** \brief Function to optimize the covariance
286  */
287  double qual_fun(size_t nv, const ubvector &x) {
288 
289  var=x[0];
290  len=x[1];
291 
292  size_t size=this->sz;
293 
294  qual=0.0;
295  for(size_t k=0;k<size;k++) {
296 
297  ubvector x2(size-1);
298  o2scl::vector_copy_jackknife((*this->px),k,x2);
299  ubvector y2(size-1);
300  o2scl::vector_copy_jackknife((*this->py),k,y2);
301 
302  // Construct the KXX matrix
303  ubmatrix KXX(size-1,size-1);
304  for(size_t irow=0;irow<size-1;irow++) {
305  for(size_t icol=0;icol<size-1;icol++) {
306  if (irow>icol) {
307  KXX(irow,icol)=KXX(icol,irow);
308  } else {
309  KXX(irow,icol)=exp(-var*pow(((*this->px)[irow]-
310  (*this->px)[icol])/len,2.0));
311  }
312  }
313  }
314 
315  // Construct the inverse of KXX
316  ubmatrix inv_KXX(size-1,size-1);
317  o2scl::permutation p(size-1);
318  int signum;
319  o2scl_linalg::LU_decomp(size-1,KXX,p,signum);
320  if (o2scl_linalg::diagonal_has_zero(size-1,KXX)) {
321  O2SCL_ERR("KXX matrix is singular in interp_krige::set().",
323  }
324  o2scl_linalg::LU_invert<ubmatrix,ubmatrix,ubmatrix_column>
325  (size-1,KXX,p,inv_KXX);
326 
327  // Inverse covariance matrix times function vector
328  this->Kinvf.resize(size-1);
329  boost::numeric::ublas::axpy_prod(inv_KXX,y2,this->Kinvf,true);
330 
331  double ypred=0.0;
332  double yact=(*this->py)[k];
333  for(size_t i=0;i<size-1;i++) {
334  ypred+=exp(-var*pow(((*this->px)[k]-x2[i])/len,2.0))*this->Kinvf[i];
335  }
336 
337  qual+=pow(yact-ypred,2.0);
338 
339  }
340 
341  return qual;
342  }
343 
344  public:
345 
346  /// Number of variance points to try
347  size_t nvar;
348 
349  /// Number of length scale points to try
350  size_t nlen;
351 
352  /// Default minimizer
354 
355  /// If true, use the full minimizer
356  bool full_min;
357 
359  ff=std::bind(std::mem_fn<double(double,double)>
361  std::placeholders::_1,std::placeholders::_2);
362  nvar=20;
363  nlen=20;
364  full_min=false;
365  mp=&def_mmin;
366  }
367 
368  /// Initialize interpolation routine
369  virtual void set_noise(size_t size, const vec_t &x, const vec2_t &y,
370  double noise_var) {
371 
372  if (full_min) {
373 
374  // Set parent data members
375  this->px=&x;
376  this->py=&y;
377  this->sz=size;
378 
379  ubvector p(2);
380  p[0]=o2scl::vector_variance(size,y);
381  p[1]=x[1]-x[0];
382 
383  multi_funct mf=std::bind
384  (std::mem_fn<double(size_t, const ubvector &)>
386  std::placeholders::_1,std::placeholders::_2);
387 
388  mp->mmin(2,p,qual,mf);
389 
390  } else {
391 
392  // Range of the coefficient parameter
393  double var_min=o2scl::vector_variance(size,y);
394  std::vector<double> diff(size-1);
395  for(size_t i=0;i<size-1;i++) {
396  diff[i]=fabs(x[i+1]-x[i]);
397  }
398  double var_ratio=1.0e2;
399 
400  // Range of the length parameter
401  double len_min=o2scl::vector_min_value
402  <std::vector<double>,double>(size-1,diff)/3.0;
403  double len_max=fabs(x[size-1]-x[0])*3.0;
404  double len_ratio=len_max/len_min;
405 
406  // Initialize to zero to prevent uninit'ed var. warnings
407  double min_qual=0.0, var_opt=0.0, len_opt=0.0;
408 
409  // Loop over the full range, finding the optimum
410  for(size_t i=0;i<nvar;i++) {
411  var=var_min*pow(var_ratio,((double)i)/((double)nvar-1));
412  for(size_t j=0;j<nlen;j++) {
413  len=len_min*pow(len_ratio,((double)j)/((double)nlen-1));
414 
415  qual=0.0;
416  for(size_t k=0;k<size;k++) {
417 
418  ubvector x2(size-1);
420  ubvector y2(size-1);
422 
423  // Construct the KXX matrix
424  ubmatrix KXX(size-1,size-1);
425  for(size_t irow=0;irow<size-1;irow++) {
426  for(size_t icol=0;icol<size-1;icol++) {
427  if (irow>icol) {
428  KXX(irow,icol)=KXX(icol,irow);
429  } else {
430  KXX(irow,icol)=exp(-var*pow((x[irow]-x[icol])/len,2.0));
431  }
432  }
433  }
434 
435  // Construct the inverse of KXX
436  ubmatrix inv_KXX(size-1,size-1);
437  o2scl::permutation p(size-1);
438  int signum;
439  o2scl_linalg::LU_decomp(size-1,KXX,p,signum);
440  if (o2scl_linalg::diagonal_has_zero(size-1,KXX)) {
441  O2SCL_ERR("KXX matrix is singular in interp_krige::set().",
443  }
444  o2scl_linalg::LU_invert<ubmatrix,ubmatrix,ubmatrix_column>
445  (size-1,KXX,p,inv_KXX);
446 
447  // Inverse covariance matrix times function vector
448  this->Kinvf.resize(size-1);
449  boost::numeric::ublas::axpy_prod(inv_KXX,y2,this->Kinvf,true);
450 
451  double ypred=0.0;
452  double yact=y[k];
453  for(size_t i=0;i<size-1;i++) {
454  ypred+=exp(-var*pow((x[k]-x2[i])/len,2.0))*this->Kinvf[i];
455  }
456 
457  qual+=pow(yact-ypred,2.0);
458 
459  }
460 
461  if ((i==0 && j==0) || qual<min_qual) {
462  var_opt=var;
463  len_opt=len;
464  min_qual=qual;
465  }
466  }
467  }
468 
469  // Now that we've optimized the covariance function,
470  // just use the parent class to interpolate
471  var=var_opt;
472  len=len_opt;
473 
474  }
475 
476  this->set_covar_noise(size,x,y,ff,noise_var);
477 
478  return;
479  }
480 
481  /// Initialize interpolation routine
482  virtual void set(size_t size, const vec_t &x, const vec2_t &y) {
483  set_noise(size,x,y,0.0);
484  return;
485  }
486 
487 
488  };
489 
490 #ifndef DOXYGEN_NO_O2NS
491 }
492 #endif
493 
494 #endif
std::function< double(size_t, const boost::numeric::ublas::vector< double > &)> multi_funct
Multi-dimensional function typedef.
Definition: multi_funct.h:45
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:1153
virtual void 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.
Definition: interp_krige.h:118
covar_func_t * df
Pointer to user-specified derivative.
Definition: interp_krige.h:90
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:2299
virtual double integ(double a, double b) const
Give the value of the integral .
Definition: interp_krige.h:209
size_t nvar
Number of variance points to try.
Definition: interp_krige.h:347
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
double qual_fun(size_t nv, const ubvector &x)
Function to optimize the covariance.
Definition: interp_krige.h:287
size_t min_size
The minimum size of the vectors to interpolate between.
Definition: interp.h:162
const double pi
Definition: constants.h:62
double integ(double x, double x1, double x2)
The integral of the covariance function.
Definition: interp_krige.h:274
covar_func_t * f
Pointer to user-specified covariance function.
Definition: interp_krige.h:86
size_t sz
Vector size.
Definition: interp.h:128
mmin_bfgs2 def_mmin
Default minimizer.
Definition: interp_krige.h:353
std::function< double(double, double)> ff
Functor for the covariance function covar()
Definition: interp_krige.h:280
invalid argument supplied by user
Definition: err_hnd.h:59
virtual int mmin(size_t nvar, vec_t &x, double &fmin, func_t &func)=0
Calculate the minimum min of func w.r.t. the array x of size nvar.
const vec_t * px
Independent vector.
Definition: interp.h:122
A class for representing permutations.
Definition: permutation.h:70
Base low-level interpolation class [abstract base].
Definition: interp.h:104
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:66
virtual void 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:178
generic failure
Definition: err_hnd.h:61
double vector_variance(size_t n, const vec_t &data, double mean)
Compute the variance with specified mean.
Definition: vec_stats.h:128
double len
The covariance function length scale.
Definition: interp_krige.h:249
requested feature not (yet) implemented
Definition: err_hnd.h:99
One-dimensional interpolation using an optimized covariance function.
Definition: interp_krige.h:238
double var
The covariance function coefficient.
Definition: interp_krige.h:252
virtual void set_noise(size_t size, const vec_t &x, const vec2_t &y, double noise_var)
Initialize interpolation routine.
Definition: interp_krige.h:369
virtual double deriv(double x0) const
Give the value of the derivative .
Definition: interp_krige.h:197
double deriv2(double x1, double x2)
The second derivative of the covariance function.
Definition: interp_krige.h:268
covar_func_t * df2
Pointer to user-specified second derivative.
Definition: interp_krige.h:94
covar_integ_t * intp
Pointer to user-specified second derivative.
Definition: interp_krige.h:98
double deriv(double x1, double x2)
The derivative of the covariance function.
Definition: interp_krige.h:263
virtual void set_covar_noise(size_t n_dim, const vec_t &x, const vec_t &y, covar_func_t &fcovar, double noise_var)
Initialize interpolation routine.
Definition: interp_krige.h:128
Multidimensional minimization [abstract base].
Definition: mmin.h:164
mmin_base * mp
Pointer to the user-specified minimizer.
Definition: interp_krige.h:283
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
const vec2_t * py
Dependent vector.
Definition: interp.h:125
virtual const char * type() const
Return the type, "interp_linear".
Definition: interp_krige.h:214
#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:356
ubvector Kinvf
Inverse covariance matrix times function vector.
Definition: interp_krige.h:82
size_t nlen
Number of length scale points to try.
Definition: interp_krige.h:350
double covar(double x1, double x2)
The covariance function.
Definition: interp_krige.h:258
virtual double deriv2(double x0) const
Give the value of the second derivative (always zero)
Definition: interp_krige.h:204
double qual
The quality factor of the optimization.
Definition: interp_krige.h:255
Multidimensional minimization by the BFGS algorithm (GSL)
Definition: mmin_bfgs2.h:387
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
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:185
static const double x1[5]
Definition: inte_qng_gsl.h:48
std::string szttos(size_t x)
Convert a size_t to a string.

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