interpm_krige.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2017-2020, 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_INTERPM_KRIGE_H
24 #define O2SCL_INTERPM_KRIGE_H
25 
26 /** \file interpm_krige.h
27  \brief File defining \ref o2scl::interpm_krige and
28  \ref o2scl::interpm_krige_nn
29 */
30 
31 #include <iostream>
32 #include <string>
33 #include <cmath>
34 
35 #include <boost/numeric/ublas/matrix.hpp>
36 #include <boost/numeric/ublas/operation.hpp>
37 
38 #include <gsl/gsl_combination.h>
39 
40 #include <o2scl/err_hnd.h>
41 #include <o2scl/vector.h>
42 #include <o2scl/vec_stats.h>
43 #include <o2scl/linear_solver.h>
44 #include <o2scl/columnify.h>
45 #include <o2scl/cholesky.h>
46 #include <o2scl/min_brent_gsl.h>
47 #include <o2scl/cblas.h>
48 
49 #ifndef DOXYGEN_NO_O2NS
50 namespace o2scl {
51 #endif
52 
53  /** \brief Multi-dimensional interpolation by kriging
54 
55  \note The set data functions for this class uses a particular
56  format, one different format than that in \ref
57  o2scl::interpm_idw . This design choice makes it easier to pass
58  vector arguments to the covariance function and the linear
59  algebra routines. The x and y objects should be of the form
60  <tt>x[n_points][n_in]</tt> and <tt>y[n_out][n_points]</tt>.
61  A separate covariance function is required for each output.
62 
63  \note This class assumes that the function specified in the
64  call to set_data() is the same as that passed to the
65  eval() functions. If this is not the case, the
66  behavior of this class is undefined.
67 
68  \note Experimental.
69  */
70  template<class vec_t=boost::numeric::ublas::vector<double>,
71  class mat_t=boost::numeric::ublas::vector<double>,
72  class mat_row_t=boost::numeric::ublas::matrix_row
73  <boost::numeric::ublas::vector<double> > >
74  class interpm_krige {
75 
76  public:
77 
81  typedef boost::numeric::ublas::matrix_column<ubmatrix> ubmatrix_column;
82  typedef boost::numeric::ublas::matrix_row<ubmatrix> ubmatrix_row;
83 
84  protected:
85 
86  /** \brief Inverse covariance matrix times function vector
87  */
88  std::vector<ubvector> Kinvf;
89 
90  public:
91 
92  interpm_krige() {
93  data_set=false;
94  verbose=0;
96  }
97 
98  /// \name Select matrix inversion method
99  //@{
100  /** \brief Method for matrix inversion
101  */
102  size_t matrix_mode;
103  /// Use Cholesky decomposition
104  static const size_t matrix_cholesky=0;
105  /// Use LU decomposition
106  static const size_t matrix_LU=1;
107  //@}
108 
109  /** \brief Verbosity parameter (default 0)
110  */
111  int verbose;
112 
113  /** \brief Initialize the data for the interpolation
114 
115  \note This function works differently than
116  \ref o2scl::interpm_idw::set_data() . See this
117  class description for more details.
118  */
119  template<class mat2_row_t, class mat2_t, class func_vec_t>
120  int set_data_noise(size_t n_in, size_t n_out, size_t n_points,
121  mat_t &user_x, mat2_t &user_y,
122  func_vec_t &fcovar,
123  const vec_t &noise_var, bool rescale=false,
124  bool err_on_fail=true) {
125 
126  if (n_points<2) {
127  O2SCL_ERR2("Must provide at least two points in ",
128  "interpm_krige::set_data_noise()",exc_efailed);
129  }
130  if (n_in<1) {
131  O2SCL_ERR2("Must provide at least one input column in ",
132  "interpm_krige::set_data_noise()",exc_efailed);
133  }
134  if (n_out<1) {
135  O2SCL_ERR2("Must provide at least one output column in ",
136  "interpm_krige::set_data_noise()",exc_efailed);
137  }
138  if (noise_var.size()<1) {
139  O2SCL_ERR2("Noise vector empty in ",
140  "interpm_krige::set_data_noise()",exc_efailed);
141  }
142  np=n_points;
143  nd_in=n_in;
144  nd_out=n_out;
145 
146  if (user_x.size1()!=n_points || user_x.size2()!=n_in) {
147  O2SCL_ERR2("Size of x not correct in ",
148  "interpm_krige::set_data_noise().",o2scl::exc_efailed);
149  }
150  std::swap(user_x,x);
151 
152  // We don't need to copy the 'y' data, but we double check
153  // that it is properly sized
154  if (user_y.size2()!=n_points || user_y.size1()!=n_out) {
155  O2SCL_ERR2("Size of y not correct in ",
156  "interpm_krige::set_data_noise().",o2scl::exc_efailed);
157  }
158 
159  rescaled=rescale;
160 
161  data_set=true;
162 
163  if (verbose>0) {
164  std::cout << "interpm_krige::set_data_noise() : Using " << n_points
165  << " points with " << nd_in << " input variables and\n\t"
166  << nd_out << " output variables." << std::endl;
167  }
168 
169  if (rescale==true) {
170  min.resize(n_in);
171  max.resize(n_in);
172  for(size_t j=0;j<n_in;j++) {
173  min[j]=x(0,j);
174  max[j]=x(0,j);
175  for(size_t i=1;i<n_points;i++) {
176  double val=x(i,j);
177  if (val>max[j]) max[j]=val;
178  if (val<min[j]) min[j]=val;
179  }
180  for(size_t i=0;i<n_points;i++) {
181  x(i,j)=((x(i,j)-min[j])/(max[j]-min[j])-0.5)*2.0;
182  }
183  }
184  if (verbose>1) {
185  std::cout << "interpm_krige::set_data_noise() rescaled."
186  << std::endl;
187  }
188  }
189 
190  Kinvf.resize(n_out);
191  size_t n_covar=fcovar.size();
192 
193  // Loop over all output functions
194  for(size_t iout=0;iout<n_out;iout++) {
195 
196  size_t icovar=iout % n_covar;
197  size_t inoise=iout & noise_var.size();
198 
199  // Select the row of the data matrix
200  mat2_row_t yiout(user_y,iout);
201 
202  // Construct the KXX matrix
203  ubmatrix KXX(n_points,n_points);
204  for(size_t irow=0;irow<n_points;irow++) {
205  mat_row_t xrow(x,irow);
206  for(size_t icol=0;icol<n_points;icol++) {
207  mat_row_t xcol(x,icol);
208  if (irow>icol) {
209  KXX(irow,icol)=KXX(icol,irow);
210  } else if (irow==icol) {
211  KXX(irow,icol)=fcovar[icovar](xrow,xcol)+
212  noise_var[inoise];
213  } else {
214  KXX(irow,icol)=fcovar[icovar](xrow,xcol);
215  }
216  }
217  }
218 
219  if (matrix_mode==matrix_LU) {
220 
221  // Construct the inverse of KXX
222  ubmatrix inv_KXX(n_points,n_points);
223  o2scl::permutation p(n_points);
224  int signum;
225  o2scl_linalg::LU_decomp(n_points,KXX,p,signum);
226  if (o2scl_linalg::diagonal_has_zero(n_points,KXX)) {
227  if (err_on_fail) {
228  O2SCL_ERR2("Matrix singular (LU method) ",
229  "in interpm_krige::set_data_noise().",
231  }
232  return 1;
233  }
234  o2scl_linalg::LU_invert<ubmatrix,ubmatrix,ubmatrix_column>
235  (n_points,KXX,p,inv_KXX);
236 
237  // Inverse covariance matrix times function vector
238  Kinvf[iout].resize(n_points);
239  o2scl_cblas::dgemv(o2scl_cblas::o2cblas_RowMajor,
240  o2scl_cblas::o2cblas_NoTrans,
241  n_points,n_points,1.0,inv_KXX,
242  yiout,0.0,Kinvf[iout]);
243  //boost::numeric::ublas::axpy_prod(inv_KXX,yiout,Kinvf[iout],true);
244 
245  } else {
246 
247  // Construct the inverse of KXX
248  int cret=o2scl_linalg::cholesky_decomp(n_points,KXX,false);
249  if (cret!=0) {
250  if (err_on_fail) {
251  O2SCL_ERR2("Matrix singular (Cholesky method) ",
252  "in interpm_krige::set_data_noise().",
254  }
255  return 2;
256  }
257  ubmatrix &inv_KXX=KXX;
258  o2scl_linalg::cholesky_invert<ubmatrix>(n_points,inv_KXX);
259 
260  // Inverse covariance matrix times function vector
261  Kinvf[iout].resize(n_points);
262  o2scl_cblas::dgemv(o2scl_cblas::o2cblas_RowMajor,
263  o2scl_cblas::o2cblas_NoTrans,
264  n_points,n_points,1.0,inv_KXX,
265  yiout,0.0,Kinvf[iout]);
266 
267  }
268 
269  if (verbose>1) {
270  std::cout << "interpm_krige::set_data_noise() finished " << iout+1
271  << " of " << n_out << "." << std::endl;
272  }
273 
274  }
275 
276  if (verbose>1) {
277  std::cout << "interpm_krige::set_data_noise() done."
278  << std::endl;
279  }
280 
281  return 0;
282  }
283 
284  /** \brief Initialize the data for the interpolation
285 
286  \note This function works differently than
287  \ref o2scl::interpm_idw::set_data() . See this
288  class description for more details.
289  */
290  template<class mat2_row_t, class mat2_t, class func_vec_t>
291  int set_data(size_t n_in, size_t n_out, size_t n_points,
292  mat_t &user_x, mat2_t &user_y,
293  func_vec_t &fcovar, bool rescale=false,
294  bool err_on_fail=true) {
295  vec_t noise_vec;
296  noise_vec.resize(1);
297  noise_vec[0]=0.0;
298  return set_data_noise<mat2_row_t,mat2_t,func_vec_t>
299  (n_in,n_out,n_points,user_x,user_y,fcovar,
300  noise_vec,rescale,err_on_fail);
301  }
302 
303  /** \brief Given covariance function \c fcovar and input vector \c x
304  store the result of the interpolation in \c y
305  */
306  template<class vec2_t, class vec3_t, class vec_func_t>
307  void eval(const vec2_t &x0, vec3_t &y0, vec_func_t &fcovar) {
308 
309  if (data_set==false) {
310  O2SCL_ERR("Data not set in interpm_krige::eval().",
311  exc_einval);
312  }
313  if (fcovar.size()==0) {
314  O2SCL_ERR("No covariance functions in interpm_krige::eval().",
315  exc_einval);
316  }
317 
318  if (rescaled) {
319 
320  // If necessary, rescale before evaluating the interpolated
321  // result
322  vec2_t x0p(nd_in);
323  for(size_t iin=0;iin<nd_in;iin++) {
324  x0p[iin]=((x0[iin]-min[iin])/(max[iin]-min[iin])-0.5)*2.0;
325  }
326 
327  // Evaluate the interpolated result
328  for(size_t iout=0;iout<nd_out;iout++) {
329  size_t icovar=iout % fcovar.size();
330  y0[iout]=0.0;
331  for(size_t ipoints=0;ipoints<np;ipoints++) {
332  mat_row_t xrow(x,ipoints);
333  double covar_val=fcovar[icovar](xrow,x0p);
334  y0[iout]+=covar_val*Kinvf[iout][ipoints];
335  }
336  }
337 
338  } else {
339 
340  // Evaluate the interpolated result
341  for(size_t iout=0;iout<nd_out;iout++) {
342  size_t icovar=iout % fcovar.size();
343  y0[iout]=0.0;
344  for(size_t ipoints=0;ipoints<np;ipoints++) {
345  mat_row_t xrow(x,ipoints);
346  y0[iout]+=fcovar[icovar](xrow,x0)*Kinvf[iout][ipoints];
347  }
348  }
349 
350  }
351 
352  return;
353 
354  }
355 
356 #ifndef DOXYGEN_INTERNAL
357 
358  protected:
359 
360  /// The number of points
361  size_t np;
362  /// The number of dimensions of the inputs
363  size_t nd_in;
364  /// The number of dimensions of the outputs
365  size_t nd_out;
366  /// The data
367  mat_t x;
368  /// True if the data has been specified
369  bool data_set;
370  /// Minimum values for rescaling
372  /// Maximum values for rescaling
374  /// True if the data needs to be rescaled
375  bool rescaled;
376 
377 #endif
378 
379  };
380 
381  /** \brief One-dimensional interpolation using an
382  optimized covariance function
383 
384  See also the \ref intp_section section of the \o2 User's guide.
385 
386  \note This class is experimental.
387  */
388  template<class vec_t=boost::numeric::ublas::vector<double>,
389  class mat_t=boost::numeric::ublas::vector<double>,
390  class mat_row_t=boost::numeric::ublas::matrix_row
391  <boost::numeric::ublas::vector<double> > >
393  public interpm_krige<vec_t,mat_t,mat_row_t> {
394 
395  public:
396 
399  typedef boost::numeric::ublas::matrix_column<ubmatrix> ubmatrix_column;
400 
401  /// Function objects for the covariance
402  std::vector<std::function<double(const mat_row_t &, const vec_t &)> >
404 
405  protected:
406 
407  /// Function objects for the covariance
408  std::vector<std::function<double(const mat_row_t &, const mat_row_t &)> >
410 
411  /// The covariance function length scale for each output function
412  std::vector<double> len;
413 
414  /// The quality factor of the optimization for each output function
415  std::vector<double> qual;
416 
417  /// If true, min and max has been set for the length parameter
419 
420  /// Minimum for length parameter range
421  double len_min;
422 
423  /// Maximum for length parameter range
424  double len_max;
425 
426  /// The covariance function
427  template<class vec2_t, class vec3_t>
428  double covar(const vec2_t &x1, const vec3_t &x2, size_t sz, double len2) {
429  double ret=0.0;
430  for(size_t i=0;i<sz;i++) {
431  ret+=pow(x1[i]-x2[i],2.0);
432  }
433  ret=exp(-ret/len2/len2/2.0);
434  return ret;
435  }
436 
437  /// Pointer to the user-specified minimizer
439 
440  /** \brief Function to optimize the covariance parameters
441  */
442  template<class vec3_t>
443  double qual_fun(double xlen, double noise_var, size_t iout,
444  vec3_t &y, int &success) {
445 
446  double ret=0.0;
447 
448  success=0;
449 
450  size_t size=this->x.size1();
451 
452  if (mode==mode_loo_cv) {
453 
454  for(size_t ell=0;ell<loo_npts;ell++) {
455 
456  // Create the new data objects, x_jk and y_jk
457  size_t row=ell*size/loo_npts;
458  matrix_view_omit_row<mat_t> x_jk(this->x,row);
459  ubvector y_jk(size-1);
460  vector_copy_jackknife(size,y,row,y_jk);
461 
462  // Now perform the matrix analysis with those objects
463 
464  // Construct the KXX matrix
465  ubmatrix KXX(size-1,size-1);
466  for(size_t irow=0;irow<size-1;irow++) {
468  for(size_t icol=0;icol<size-1;icol++) {
470  if (irow>icol) {
471  KXX(irow,icol)=KXX(icol,irow);
472  } else {
473  KXX(irow,icol)=covar<
476  (xrow,xcol,this->nd_in,xlen);
477  if (irow==icol) KXX(irow,icol)+=noise_var;
478  }
479  }
480  }
481 
482  if (this->matrix_mode==this->matrix_cholesky) {
483 
484  // Construct the inverse of KXX
485  if (verbose>2) {
486  std::cout << "Performing Cholesky decomposition with size "
487  << size-1 << std::endl;
488  }
489  int cret=o2scl_linalg::cholesky_decomp(size-1,KXX,false);
490  if (cret!=0) {
491  success=1;
492  return 1.0e99;
493  }
494  if (verbose>2) {
495  std::cout << "Performing matrix inversion with size "
496  << size-1 << std::endl;
497  }
498 
499  o2scl_linalg::cholesky_invert<ubmatrix>(size-1,KXX);
500 
501  // Inverse covariance matrix times function vector
502  this->Kinvf[iout].resize(size-1);
503  o2scl_cblas::dgemv(o2scl_cblas::o2cblas_RowMajor,
504  o2scl_cblas::o2cblas_NoTrans,
505  size-1,size-1,1.0,KXX,
506  y,0.0,this->Kinvf[iout]);
507 
508  } else {
509 
510  // Construct the inverse of KXX
511  ubmatrix inv_KXX(size-1,size-1);
512  o2scl::permutation p(size-1);
513  int signum;
514  if (verbose>2) {
515  std::cout << "Performing LU decomposition with size "
516  << size-1 << std::endl;
517  }
518  o2scl_linalg::LU_decomp(size-1,KXX,p,signum);
519  if (o2scl_linalg::diagonal_has_zero(size-1,KXX)) {
520  success=1;
521  return 1.0e99;
522  }
523  if (verbose>2) {
524  std::cout << "Performing matrix inversion with size "
525  << size-1 << std::endl;
526  }
527  o2scl_linalg::LU_invert<ubmatrix,ubmatrix,ubmatrix_column>
528  (size-1,KXX,p,inv_KXX);
529 
530  // Inverse covariance matrix times function vector
531  this->Kinvf[iout].resize(size-1);
532  o2scl_cblas::dgemv(o2scl_cblas::o2cblas_RowMajor,
533  o2scl_cblas::o2cblas_NoTrans,
534  size-1,size-1,1.0,inv_KXX,
535  y,0.0,this->Kinvf[iout]);
536 
537  }
538 
539  double ypred=0.0;
540  double yact=y[row];
541  for(size_t i=0;i<size-1;i++) {
543  mat_row_t xcol(this->x,row);
544  ypred+=covar<matrix_row_gen<matrix_view_omit_row<mat_t> >,
545  mat_row_t>(xrow,xcol,this->nd_in,xlen)*
546  this->Kinvf[iout][i];
547  }
548 
549  std::cout << "act,pred: " << yact << " " << ypred << std::endl;
550 
551  // Measure the quality with a chi-squared like function
552  ret+=pow(yact-ypred,2.0);
553 
554 
555  // Proceed to next point to omit
556  }
557  std::cout << "ret: " << ret << std::endl;
558 
559  } else if (mode==mode_max_lml || mode==mode_final) {
560 
561  if (verbose>2) {
562  std::cout << "Creating covariance matrix with size "
563  << size << std::endl;
564  }
565 
566  // Construct the KXX matrix
567  ubmatrix KXX(size,size);
568  for(size_t irow=0;irow<size;irow++) {
569  mat_row_t xrow(this->x,irow);
570  for(size_t icol=0;icol<size;icol++) {
571  mat_row_t xcol(this->x,icol);
572  if (irow>icol) {
573  KXX(irow,icol)=KXX(icol,irow);
574  } else {
575  KXX(irow,icol)=covar<mat_row_t,mat_row_t>(xrow,xcol,
576  this->nd_in,xlen);
577  if (irow==icol) KXX(irow,icol)+=noise_var;
578  }
579  }
580  }
581 
582  // Note: We have to use LU here because O2scl doesn't yet
583  // have a lndet() function for Cholesky decomp
584 
585  double lndet;
586  if (this->matrix_mode==this->matrix_cholesky) {
587 
588  // Construct the inverse of KXX
589  if (verbose>2) {
590  std::cout << "Performing Cholesky decomposition with size "
591  << size << std::endl;
592  }
593  int cret=o2scl_linalg::cholesky_decomp(size,KXX,false);
594  if (cret!=0) {
595  success=1;
596  return 1.0e99;
597  }
598  if (verbose>2) {
599  std::cout << "Performing matrix inversion with size "
600  << size << std::endl;
601  }
602  o2scl_linalg::cholesky_invert<ubmatrix>(size,KXX);
603 
604  lndet=o2scl_linalg::cholesky_lndet<ubmatrix>(size,KXX);
605 
606  // Inverse covariance matrix times function vector
607  this->Kinvf[iout].resize(size);
608  o2scl_cblas::dgemv(o2scl_cblas::o2cblas_RowMajor,
609  o2scl_cblas::o2cblas_NoTrans,
610  size,size,1.0,KXX,
611  y,0.0,this->Kinvf[iout]);
612 
613  } else {
614 
615  // Construct the inverse of KXX
616  ubmatrix inv_KXX(size,size);
617  o2scl::permutation p(size);
618  int signum;
619  if (verbose>2) {
620  std::cout << "Performing LU decomposition with size "
621  << size << std::endl;
622  }
623  o2scl_linalg::LU_decomp(size,KXX,p,signum);
624  if (o2scl_linalg::diagonal_has_zero(size,KXX)) {
625  success=1;
626  return 1.0e99;
627  }
628  if (verbose>2) {
629  std::cout << "Performing matrix inversion with size "
630  << size << std::endl;
631  }
632  o2scl_linalg::LU_invert<ubmatrix,ubmatrix,ubmatrix_column>
633  (size,KXX,p,inv_KXX);
634 
635  lndet=o2scl_linalg::LU_lndet<ubmatrix>(size,KXX);
636 
637  // Inverse covariance matrix times function vector
638  this->Kinvf[iout].resize(size);
639  o2scl_cblas::dgemv(o2scl_cblas::o2cblas_RowMajor,
640  o2scl_cblas::o2cblas_NoTrans,
641  size,size,1.0,inv_KXX,
642  y,0.0,this->Kinvf[iout]);
643 
644  }
645 
646  if (mode==mode_max_lml) {
647  // Compute the log of the marginal likelihood, without
648  // the constant term
649  for(size_t i=0;i<size;i++) {
650  ret+=0.5*y[i]*this->Kinvf[iout][i];
651  }
652  ret+=0.5*lndet;
653  }
654 
655  }
656 
657  return ret;
658  }
659 
660  public:
661 
663  nlen=20;
664  full_min=false;
665  mp=&def_min;
666  verbose=0;
668  loo_npts=100;
669  len_guess_set=false;
670  }
671 
672  /// \name Function to minimize and various option
673  //@{
674  /// Leave-one-out cross validation
675  static const size_t mode_loo_cv=1;
676  /// Minus Log-marginal-likelihood
677  static const size_t mode_max_lml=2;
678  /// No optimization (for internal use)
679  static const size_t mode_final=10;
680  /// Function to minimize (default \ref mode_loo_cv)
681  size_t mode;
682  /// Number of points to test for cross validation (default 100)
683  size_t loo_npts;
684  ///@}
685 
686  /// Verbosity parameter
687  int verbose;
688 
689  /** \brief Number of length scale points to try when full minimizer
690  is not used (default 20)
691  */
692  size_t nlen;
693 
694  /// Default minimizer
696 
697  /// If true, use the full minimizer
698  bool full_min;
699 
700  /** \brief Set the range for the length parameter
701  */
702  void set_len_range(double min2, double max2) {
703  len_guess_set=true;
704  len_min=min2;
705  len_max=max2;
706  return;
707  }
708 
709  /// Initialize interpolation routine
710  template<class mat2_row_t, class mat2_t, class vec2_t,
711  class vec3_t>
712  int set_data_noise(size_t n_in, size_t n_out, size_t n_points,
713  mat_t &user_x, mat2_t &user_y,
714  const vec2_t &noise_var, const vec3_t &len_precompute,
715  bool rescale=false, bool err_on_fail=true) {
716 
717  if (n_points<2) {
718  O2SCL_ERR2("Must provide at least two points in ",
719  "interpm_krige_optim::set_data_noise()",exc_efailed);
720  }
721  if (n_in<1) {
722  O2SCL_ERR2("Must provide at least one input column in ",
723  "interpm_krige_optim::set_data_noise()",exc_efailed);
724  }
725  if (n_out<1) {
726  O2SCL_ERR2("Must provide at least one output column in ",
727  "interpm_krige_optim::set_data_noise()",exc_efailed);
728  }
729  if (noise_var.size()<1) {
730  O2SCL_ERR2("Noise vector empty in ",
731  "interpm_krige::set_data_noise()",exc_efailed);
732  }
733 
734  // Set parent data members
735  this->np=n_points;
736  this->nd_in=n_in;
737  this->nd_out=n_out;
738  std::swap(this->x,user_x);
739  this->rescaled=rescale;
740  this->data_set=true;
741 
742  if (verbose>0) {
743  std::cout << "interpm_krige_optim::set_data_noise() : Using "
744  << n_points
745  << " points with\n " << n_in << " input variables and "
746  << this->nd_out << " output variables." << std::endl;
747  }
748 
749  // Get max and min of all parameters
750  this->min.resize(n_in);
751  this->max.resize(n_in);
752  for(size_t j=0;j<n_in;j++) {
753  this->min[j]=(this->x)(0,j);
754  this->max[j]=(this->x)(0,j);
755  for(size_t i=1;i<n_points;i++) {
756  double val=(this->x)(i,j);
757  if (val>this->max[j]) this->max[j]=val;
758  if (val<this->min[j]) this->min[j]=val;
759  }
760  }
761 
762  if (rescale==true) {
763  for(size_t j=0;j<n_in;j++) {
764  for(size_t i=0;i<n_points;i++) {
765  (this->x)(i,j)=(((this->x)(i,j)-this->min[j])/
766  (this->max[j]-this->min[j])-0.5)*2.0;
767  }
768  }
769  if (verbose>1) {
770  std::cout << "interpm_krige_optim::set_data_noise() rescaled."
771  << std::endl;
772  }
773  }
774 
775  int success=0;
776 
777  this->Kinvf.resize(n_out);
778  qual.resize(n_out);
779  len.resize(n_out);
780  ff1.resize(n_out);
781  ff2.resize(n_out);
782 
783  // Loop over all output functions
784  for(size_t iout=0;iout<n_out;iout++) {
785 
786  // Select the row of the data matrix
787  mat2_row_t yiout(user_y,iout);
788 
789  if (iout<len_precompute.size()) {
790 
791  if (verbose>1) {
792  std::cout << "interp_krige_optim: precomputed length "
793  << len_precompute[iout] << std::endl;
794  }
795  len[iout]=len_precompute[iout];
796  size_t mode_temp=mode;
798  qual[iout]=qual_fun(len[iout],noise_var[iout],iout,yiout,success);
799  mode=mode_temp;
800 
801  } else if (full_min) {
802 
803  if (verbose>1) {
804  std::cout << "interp_krige_optim: full minimization"
805  << std::endl;
806  }
807 
808  double len_opt;
809  // Choose average distance for first guess
810  if (this->rescaled) {
811  len_opt=(this->max[iout]-this->min[iout])/((double)this->np);
812  } else {
813  len_opt=2.0/((double)this->np);
814  }
815 
816  funct mf=std::bind
817  (std::mem_fn<double(double,double,size_t,mat2_row_t &,int &)>
818  (&interpm_krige_optim<vec_t,mat_t,
819  mat_row_t>::qual_fun<mat2_row_t>),
820  this,std::placeholders::_1,noise_var[iout],iout,yiout,
821  std::ref(success));
822 
823  mp->min(len_opt,qual[iout],mf);
824  len[iout]=len_opt;
825 
826  if (success!=0) {
827  if (err_on_fail) {
828  O2SCL_ERR2("Minimization failed in ",
829  "interp_krige_optim::set_noise().",
831  }
832  }
833 
834  } else {
835 
836  if (verbose>1) {
837  std::cout << "interp_krige_optim::set_data_noise() : "
838  << "simple minimization" << std::endl;
839  }
840 
841  if (len_guess_set==false) {
842  double len_avg;
843  // Choose average distance for first guess
844  if (this->rescaled) {
845  len_avg=(this->max[iout]-this->min[iout])/((double)this->np);
846  } else {
847  len_avg=2.0/((double)this->np);
848  }
849 
850  len_min=len_avg/1.0e2;
851  len_max=len_avg*1.0e2;
852  }
853  double len_ratio=len_max/len_min;
854 
855  if (verbose>1) {
856  std::cout << "iout, len (min,max,step): " << iout
857  << " " << len_min << " " << len_max << " "
858  << pow(len_ratio,((double)1)/((double)nlen-1))
859  << std::endl;
860  }
861 
862  // Initialize to zero to prevent uninit'ed var. warnings
863  double min_qual=0.0, len_opt=0.0;
864 
865  if (verbose>1) {
866  std::cout << "ilen len qual fail min_qual len_opt" << std::endl;
867  }
868 
869  // Loop over the full range, finding the optimum
870  bool min_set=false;
871  for(size_t j=0;j<nlen;j++) {
872  double xlen=len_min*pow(len_ratio,((double)j)/((double)nlen-1));
873 
874  success=0;
875  qual[iout]=qual_fun(xlen,noise_var[iout],iout,yiout,success);
876 
877  if (success==0 && (min_set==false || qual[iout]<min_qual)) {
878  len_opt=xlen;
879  min_qual=qual[iout];
880  min_set=true;
881  }
882 
883  if (verbose>1) {
884  std::cout << "interp_krige_optim: ";
885  std::cout.width(2);
886  std::cout << j << " " << xlen << " " << qual[iout] << " "
887  << success << " " << min_qual << " "
888  << len_opt << std::endl;
889  }
890 
891  }
892 
893  if (verbose>1) {
894  std::cout << "interp_krige_optim: ";
895  std::cout.width(2);
896  std::cout << " " << len_opt << " " << min_qual << std::endl;
897  }
898  qual[iout]=qual_fun(len_opt,noise_var[iout],iout,yiout,success);
899 
900  len[iout]=len_opt;
901 
902 
903  }
904 
905  ff1[iout]=std::bind(std::mem_fn<double(const mat_row_t &,
906  const mat_row_t &,
907  size_t,double)>
908  (&interpm_krige_optim<vec_t,mat_t,
909  mat_row_t>::covar<mat_row_t,
910  mat_row_t>),this,
911  std::placeholders::_1,std::placeholders::_2,
912  n_in,len[iout]);
913  ff2[iout]=std::bind(std::mem_fn<double(const mat_row_t &,
914  const vec_t &,
915  size_t,double)>
916  (&interpm_krige_optim<vec_t,mat_t,
917  mat_row_t>::covar<mat_row_t,
918  vec_t>),this,
919  std::placeholders::_1,std::placeholders::_2,
920  n_in,len[iout]);
921  }
922 
923  return 0;
924  }
925 
926  /** \brief Initialize the data for the interpolation
927 
928  \note This function works differently than
929  \ref o2scl::interpm_idw::set_data() . See this
930  class description for more details.
931  */
932  template<class mat2_row_t, class mat2_t, class vec2_t>
933  int set_data(size_t n_in, size_t n_out, size_t n_points,
934  mat_t &user_x, mat2_t &user_y,
935  const vec2_t &len_precompute,
936  bool rescale=false, bool err_on_fail=true) {
937  vec_t noise_vec;
938  noise_vec.resize(1);
939  noise_vec[0]=0.0;
940  return set_data_noise<mat2_row_t,mat2_t,vec_t,vec2_t>
941  (n_in,n_out,n_points,user_x,
942  user_y,noise_vec,len_precompute,rescale,err_on_fail);
943  }
944 
945  };
946 
947  /** \brief Multi-dimensional interpolation by kriging with
948  nearest-neighbor
949 
950  \note This class assumes that the function specified in the
951  call to set_data() is the same as that passed to the
952  eval() functions. If this is not the case, the
953  behavior of this class is undefined.
954 
955  \note Experimental.
956  */
957  template<class vec_t=boost::numeric::ublas::vector<double>,
958  class mat_t=boost::numeric::ublas::vector<double> >
959  class interpm_krige_nn : public interpm_krige<vec_t,mat_t> {
960 
961  public:
962 
966 
967  interpm_krige_nn() {
968  data_set=false;
969  verbose=0;
970  }
971 
972  /** \brief Verbosity parameter (default 0)
973  */
974  int verbose;
975 
976  /** \brief Initialize the data for the interpolation
977  */
978  template<class mat2_t, class vec_func_t>
979  void set_data(size_t n_in, size_t n_out, size_t n_points,
980  mat_t &user_x, mat_t &user_y,
981  vec_func_t &fcovar, size_t order) {
982 
983  if (n_points<2) {
984  O2SCL_ERR2("Must provide at least two points in ",
985  "interpm_krige_nn::set_data()",exc_efailed);
986  }
987  if (n_in<1) {
988  O2SCL_ERR2("Must provide at least one input column in ",
989  "interpm_krige_nn::set_data()",exc_efailed);
990  }
991  if (n_out<1) {
992  O2SCL_ERR2("Must provide at least one output column in ",
993  "interpm_krige_nn::set_data()",exc_efailed);
994  }
995  np=n_points;
996  this->nd_in=n_in;
997  nd_out=n_out;
998  x.resize(n_points);
999  n_order=order;
1000  for(size_t i=0;i<n_points;i++) {
1001  if (x[i].size()!=n_in) {
1002  O2SCL_ERR2("Size of x not correct in ",
1003  "interpm_krige_nn::set_data().",o2scl::exc_efailed);
1004  }
1005  std::swap(x[i],x[i]);
1006  }
1007  y.resize(n_out);
1008  for(size_t i=0;i<n_out;i++) {
1009  if (y[i].size()!=n_points) {
1010  O2SCL_ERR2("Size of y not correct in ",
1011  "interpm_krige_nn::set_data().",o2scl::exc_efailed);
1012  }
1013  std::swap(y[i],y[i]);
1014  }
1015  data_set=true;
1016 
1017  if (verbose>0) {
1018  std::cout << "interpm_krige_nn::set_data() : Using " << n_points
1019  << " points with " << nd_in << " input variables and\n\t"
1020  << nd_out << " output variables and order "
1021  << n_order << " ." << std::endl;
1022  }
1023 
1024  return;
1025  }
1026 
1027  /** \brief Given covariance function \c fcovar and input vector \c x
1028  store the result of the interpolation in \c y
1029  */
1030  template<class vec2_t, class vec3_t, class vec_func_t>
1031  void eval(const vec2_t &x0, vec3_t &y0,
1032  vec_func_t &fcovar) const {
1033 
1034  if (data_set==false) {
1035  O2SCL_ERR("Data not set in interpm_krige_nn::eval().",
1036  exc_einval);
1037  }
1038 
1039  y0.resize(nd_out);
1040 
1041  // Loop over all output functions
1042  for(size_t iout=0;iout<nd_out;iout++) {
1043 
1044  // Find points closest to requested point, as defined
1045  // by the negative covariance for this output function
1046  ubvector dists(np);
1047  for(size_t ip=0;ip<np;ip++) {
1048  dists[ip]=-fcovar[iout](x0,x[ip]);
1049  }
1050 
1051  // Empty index vector (resized by the vector_smallest_index
1052  // function)
1053  ubvector_size_t index;
1054  o2scl::vector_smallest_index<ubvector,double,ubvector_size_t>
1055  (np,dists,n_order,index);
1056 
1057  // Construct subset of function values for nearest neighbors
1058  ubvector func(n_order);
1059  for(size_t io=0;io<n_order;io++) {
1060  func[io]=y[iout][index[io]];
1061  }
1062 
1063  // Construct the nearest neighbor KXX matrix
1064  ubmatrix KXX(n_order,n_order);
1065  for(size_t irow=0;irow<n_order;irow++) {
1066  for(size_t icol=0;icol<n_order;icol++) {
1067  if (irow>icol) {
1068  KXX(irow,icol)=KXX(icol,irow);
1069  } else {
1070  KXX(irow,icol)=fcovar[iout](x[index[irow]],
1071  x[index[icol]]);
1072  }
1073  }
1074  }
1075 
1076  // Construct the inverse of KXX
1078  ubmatrix &inv_KXX=KXX;
1079  o2scl_linalg::cholesky_invert<ubmatrix>(n_order,inv_KXX);
1080 
1081  // Inverse covariance matrix times function vector
1083  boost::numeric::ublas::axpy_prod(inv_KXX,func,Kinvf,true);
1084 
1085  // Comput the final result
1086  y0[iout]=0.0;
1087  for(size_t ipoints=0;ipoints<n_order;ipoints++) {
1088  y0[iout]+=-dists[index[ipoints]]*Kinvf[ipoints];
1089  }
1090 
1091  }
1092 
1093  return;
1094 
1095  }
1096 
1097  /** \brief Find a set of linearly independent points
1098 
1099  Given a point \c x, a covariance function
1100  \c fcovar, the index of the output function
1101  \c iout, and an array specifying the closest points
1102  \c index, this function produces a list of
1103  */
1104  template<class vec2_t, class vec_func_t>
1105  void find_lin_indep(const vec2_t &x2, size_t iout,
1106  vec_func_t &fcovar,
1107  ubvector_size_t &index,
1108  ubvector_size_t &indep) const {
1109 
1110  if (x2.size()<nd_in || index.size()<np || indep.size()<n_order
1111  || iout>=nd_out) {
1112  std::cout << x2.size() << " " << nd_in << std::endl;
1113  std::cout << index.size() << " " << np << std::endl;
1114  std::cout << indep.size() << " " << n_order << std::endl;
1115  std::cout << iout << " " << nd_out << std::endl;
1116  O2SCL_ERR("Vectors not of correct size in find_lin_indep().",
1118  }
1119 
1120  bool done=false;
1121  while (done==false) {
1122 
1123  // Construct subset of function values for nearest neighbors
1124  ubvector func(n_order);
1125  for(size_t io=0;io<n_order;io++) {
1126  func[io]=y[iout][index[indep[io]]];
1127  }
1128 
1129  // Construct the nearest neighbor KXX matrix
1130  ubmatrix KXX(n_order,n_order);
1131  for(size_t irow=0;irow<n_order;irow++) {
1132  for(size_t icol=0;icol<n_order;icol++) {
1133  if (irow>icol) {
1134  KXX(irow,icol)=KXX(icol,irow);
1135  } else {
1136  KXX(irow,icol)=fcovar[iout](x2[index[indep[irow]]],
1137  x2[index[indep[icol]]]);
1138  }
1139  }
1140  }
1141 
1142  // Construct the inverse of KXX
1143  int cret=o2scl_linalg::cholesky_decomp(n_order,KXX);
1144  if (cret==0) {
1145  done=true;
1146  } else {
1147  if (verbose>1) {
1148  std::cout << "Finding new independent rows." << std::endl;
1149  for(size_t j=0;j<n_order;j++) {
1150  std::cout << indep[j] << " " << KXX(j,j) << std::endl;
1151  }
1152  }
1154  double>(indep);
1155  if (verbose>1) {
1156  std::cout << "Max is: " << max << std::endl;
1157  }
1158  for(size_t j=0;j<n_order;j++) {
1159  if (KXX(j,j)==0.0) {
1160  if (max+1<np) {
1161  if (verbose>1) {
1162  std::cout << "Entry " << j << " is zero so replacing "
1163  << "entry with " << max+1 << std::endl;
1164  }
1165  indep[j]=max+1;
1166  max++;
1167  } else {
1168  O2SCL_ERR3("Failed to find set of independent points ",
1169  "in interpm_krige_nn::find_lin_indep",
1170  "(const vec2_t &, size_t).",
1172  }
1173  }
1174  }
1175  }
1176 
1177  }
1178  return;
1179  }
1180 
1181  /** \brief Given covariance function \c fcovar and input vector \c x
1182  return the result of the interpolation for function with
1183  index \c iout
1184  */
1185  template<class vec2_t, class vec_func_t>
1186  double eval(const vec2_t &x2, size_t iout,
1187  vec_func_t &fcovar) const {
1188 
1189  if (data_set==false) {
1190  O2SCL_ERR("Data not set in interpm_krige_nn::eval().",
1191  exc_einval);
1192  }
1193 
1194  double ret;
1195 
1196  // Find points closest to requested point, as defined
1197  // by the negative covariance for this output function
1198  ubvector dists(np);
1199  for(size_t ip=0;ip<np;ip++) {
1200  dists[ip]=-fcovar[iout](x2,x2[ip]);
1201  }
1202 
1203  ubvector_size_t index(np);
1204  o2scl::vector_sort_index<ubvector,ubvector_size_t>(np,dists,index);
1205 
1206  // Vector for storing the indexes in the index array which
1207  // will store the closest n_order points which are
1208  // linearly independent
1209  ubvector_size_t indep(n_order);
1210  for(size_t io=0;io<n_order;io++) {
1211  indep[io]=io;
1212  }
1213 
1214  find_lin_indep(x2,iout,fcovar,index,indep);
1215 
1216  // Construct subset of function values for nearest neighbors
1217  ubvector func(n_order);
1218  for(size_t io=0;io<n_order;io++) {
1219  func[io]=y[iout][index[indep[io]]];
1220  }
1221 
1222  // Construct the nearest neighbor KXX matrix
1223  ubmatrix KXX(n_order,n_order);
1224  for(size_t irow=0;irow<n_order;irow++) {
1225  for(size_t icol=0;icol<n_order;icol++) {
1226  if (irow>icol) {
1227  KXX(irow,icol)=KXX(icol,irow);
1228  } else {
1229  KXX(irow,icol)=fcovar[iout](x2[index[indep[irow]]],
1230  x2[index[indep[icol]]]);
1231  }
1232  }
1233  }
1234 
1235  // Construct the inverse of KXX
1237  ubmatrix &inv_KXX=KXX;
1238  o2scl_linalg::cholesky_invert<ubmatrix>(n_order,inv_KXX);
1239 
1240  // Inverse covariance matrix times function vector
1242  boost::numeric::ublas::axpy_prod(inv_KXX,func,Kinvf,true);
1243 
1244  // Comput the final result
1245  ret=0.0;
1246  for(size_t ipoints=0;ipoints<n_order;ipoints++) {
1247  ret+=-dists[index[indep[ipoints]]]*Kinvf[ipoints];
1248  }
1249 
1250  return ret;
1251 
1252  }
1253 
1254  /** \brief Compute a quality factor for interpolation
1255  using jackknife resampling
1256  */
1257  template<class vec2_t, class func_vec_t>
1258  double eval_jackknife(const vec2_t &x2, size_t iout,
1259  func_vec_t &fcovar) const {
1260 
1261  if (data_set==false) {
1262  O2SCL_ERR("Data not set in interpm_krige_nn::eval_jackknife().",
1263  exc_einval);
1264  }
1265 
1266  double qual=0.0;
1267 
1268  // Interpolated function value inside jackknife loop
1269  double ytmp;
1270 
1271  // For a distance measurement, just use the the negative
1272  // covariance for this output function
1273  ubvector dists(np);
1274  for(size_t ip=0;ip<np;ip++) {
1275  dists[ip]=-fcovar[iout](x2,x2[ip]);
1276  }
1277 
1278  // Create an index array which sorts by distance
1279  ubvector_size_t index(np);
1280  o2scl::vector_sort_index<ubvector,ubvector_size_t>(np,dists,index);
1281 
1282  // Vector for storing the indexes in the index array which
1283  // will store the closest n_order points which are
1284  // linearly independent
1285  ubvector_size_t indep(n_order);
1286  for(size_t io=0;io<n_order;io++) {
1287  indep[io]=io;
1288  }
1289 
1290  // -------------------------------------------------------------
1291  // Before the jackknife loop, we want to create a full
1292  // set of n_order linearly independent points
1293 
1294  find_lin_indep(x2,iout,fcovar,index,indep);
1295 
1296  // -------------------------------------------------------------
1297  // Now, the jackknife loop, removing one point at a time
1298 
1299  for(size_t jk=0;jk<n_order;jk++) {
1300 
1301  if (verbose>0) {
1302  std::cout << "Jackknife: " << jk << " matching function value "
1303  << y[iout][index[jk]] << std::endl;
1304  }
1305 
1306  ubvector_size_t indep_jk;
1307  vector_copy_jackknife(indep,jk,indep_jk);
1308 
1309  // Construct subset of function values for nearest neighbors
1310  ubvector func(n_order-1);
1311  for(size_t io=0;io<n_order-1;io++) {
1312  func[io]=y[iout][index[indep_jk[io]]];
1313  }
1314 
1315  // Construct the nearest neighbor KXX matrix
1316  ubmatrix KXX(n_order-1,n_order-1);
1317  for(size_t irow=0;irow<n_order-1;irow++) {
1318  for(size_t icol=0;icol<n_order-1;icol++) {
1319  if (irow>icol) {
1320  KXX(irow,icol)=KXX(icol,irow);
1321  } else {
1322  KXX(irow,icol)=fcovar[iout](x2[index[indep_jk[irow]]],
1323  x2[index[indep_jk[icol]]]);
1324  }
1325  }
1326  }
1327 
1328  // Construct the inverse of KXX
1330  ubmatrix &inv_KXX=KXX;
1331  o2scl_linalg::cholesky_invert<ubmatrix>(n_order-1,inv_KXX);
1332 
1333  // Inverse covariance matrix times function vector
1334  ubvector Kinvf(n_order-1);
1335  boost::numeric::ublas::axpy_prod(inv_KXX,func,Kinvf,true);
1336 
1337  // Comput the final result
1338  ytmp=0.0;
1339  for(size_t ipoints=0;ipoints<n_order-1;ipoints++) {
1340  ytmp+=-dists[index[indep_jk[ipoints]]]*Kinvf[ipoints];
1341  }
1342 
1343  // Add the squared deviation to y[iout]
1344  qual+=pow(y[iout][index[jk]]-ytmp,2.0);
1345 
1346  if (verbose>0) {
1347  std::cout << "Original value: "
1348  << y[iout][index[jk]] << " interpolated: "
1349  << ytmp << std::endl;
1350  }
1351 
1352  // End of jackknife loop
1353  }
1354 
1355  return qual;
1356  }
1357 
1358 
1359 #ifndef DOXYGEN_INTERNAL
1360 
1361  protected:
1362 
1363  /// The order of the interpolation (specified by \ref set_data() )
1364  size_t n_order;
1365  /// The number of points
1366  size_t np;
1367  /// The number of dimensions of the inputs
1368  size_t nd_in;
1369  /// The number of dimensions of the outputs
1370  size_t nd_out;
1371  /// A vector of pointers holding the data
1372  std::vector<vec_t> x;
1373  /// A vector of pointers holding the data
1374  std::vector<vec_t> y;
1375  /// True if the data has been specified
1376  bool data_set;
1377 
1378 #endif
1379 
1380  };
1381 
1382 #ifndef DOXYGEN_NO_O2NS
1383 }
1384 #endif
1385 
1386 #endif
1387 
1388 
1389 
boost::numeric::ublas::matrix< double >
o2scl::interpm_krige::eval
void eval(const vec2_t &x0, vec3_t &y0, vec_func_t &fcovar)
Given covariance function fcovar and input vector x store the result of the interpolation in y.
Definition: interpm_krige.h:307
o2scl_cblas::dgemv
void dgemv(const enum o2cblas_order order, const enum o2cblas_transpose TransA, const size_t M, const size_t N, const double alpha, const mat_t &A, const vec_t &X, const double beta, vec2_t &Y)
Compute .
Definition: cblas_base.h:245
o2scl::exc_esing
@ exc_esing
apparent singularity detected
Definition: err_hnd.h:93
o2scl::interpm_krige::matrix_mode
size_t matrix_mode
Method for matrix inversion.
Definition: interpm_krige.h:102
o2scl::interpm_krige::data_set
bool data_set
True if the data has been specified.
Definition: interpm_krige.h:369
o2scl::permutation
A class for representing permutations.
Definition: permutation.h:70
o2scl::interpm_krige_nn::set_data
void set_data(size_t n_in, size_t n_out, size_t n_points, mat_t &user_x, mat_t &user_y, vec_func_t &fcovar, size_t order)
Initialize the data for the interpolation.
Definition: interpm_krige.h:979
o2scl::interpm_krige_optim::set_len_range
void set_len_range(double min2, double max2)
Set the range for the length parameter.
Definition: interpm_krige.h:702
boost::numeric::ublas::vector< double >
o2scl::interpm_krige_optim::mode_final
static const size_t mode_final
No optimization (for internal use)
Definition: interpm_krige.h:679
o2scl::interpm_krige_optim::len_guess_set
bool len_guess_set
If true, min and max has been set for the length parameter.
Definition: interpm_krige.h:418
o2scl::interpm_krige_optim
One-dimensional interpolation using an optimized covariance function.
Definition: interpm_krige.h:392
o2scl::matrix_view_omit_row
Construct a view of a matrix omtting a specified row.
Definition: vector.h:2756
O2SCL_ERR3
#define O2SCL_ERR3(d, d2, d3, n)
Set an error, three-string version.
Definition: err_hnd.h:286
o2scl::interpm_krige_optim::mode_max_lml
static const size_t mode_max_lml
Minus Log-marginal-likelihood.
Definition: interpm_krige.h:677
o2scl::exc_efailed
@ exc_efailed
generic failure
Definition: err_hnd.h:61
O2SCL_ERR2
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
o2scl::interpm_krige::max
ubvector max
Maximum values for rescaling.
Definition: interpm_krige.h:373
o2scl::interpm_krige_nn::data_set
bool data_set
True if the data has been specified.
Definition: interpm_krige.h:1376
o2scl
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
Definition: anneal.h:42
o2scl_linalg::LU_decomp
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
o2scl::interpm_krige::rescaled
bool rescaled
True if the data needs to be rescaled.
Definition: interpm_krige.h:375
o2scl::interpm_krige_optim::qual
std::vector< double > qual
The quality factor of the optimization for each output function.
Definition: interpm_krige.h:415
o2scl::interpm_krige_optim::mode
size_t mode
Function to minimize (default mode_loo_cv)
Definition: interpm_krige.h:681
o2scl::interpm_krige::nd_out
size_t nd_out
The number of dimensions of the outputs.
Definition: interpm_krige.h:365
o2scl::interpm_krige_nn::y
std::vector< vec_t > y
A vector of pointers holding the data.
Definition: interpm_krige.h:1374
o2scl::vector_max_value
data_t vector_max_value(size_t n, const vec_t &data)
Compute the maximum of the first n elements of a vector.
Definition: vector.h:1121
o2scl::interpm_krige_nn
Multi-dimensional interpolation by kriging with nearest-neighbor.
Definition: interpm_krige.h:959
o2scl::interpm_krige::Kinvf
std::vector< ubvector > Kinvf
Inverse covariance matrix times function vector.
Definition: interpm_krige.h:88
o2scl::interpm_krige::np
size_t np
The number of points.
Definition: interpm_krige.h:361
o2scl::interpm_krige_optim::mode_loo_cv
static const size_t mode_loo_cv
Leave-one-out cross validation.
Definition: interpm_krige.h:675
o2scl::interpm_krige::set_data
int set_data(size_t n_in, size_t n_out, size_t n_points, mat_t &user_x, mat2_t &user_y, func_vec_t &fcovar, bool rescale=false, bool err_on_fail=true)
Initialize the data for the interpolation.
Definition: interpm_krige.h:291
o2scl::interpm_krige_optim::covar
double covar(const vec2_t &x1, const vec3_t &x2, size_t sz, double len2)
The covariance function.
Definition: interpm_krige.h:428
o2scl::interpm_krige_nn::nd_out
size_t nd_out
The number of dimensions of the outputs.
Definition: interpm_krige.h:1370
o2scl_inte_qng_coeffs::x2
static const double x2[5]
Definition: inte_qng_gsl.h:66
o2scl::interpm_krige::nd_in
size_t nd_in
The number of dimensions of the inputs.
Definition: interpm_krige.h:363
o2scl::interpm_krige_nn::eval
void eval(const vec2_t &x0, vec3_t &y0, vec_func_t &fcovar) const
Given covariance function fcovar and input vector x store the result of the interpolation in y.
Definition: interpm_krige.h:1031
o2scl::interpm_krige_optim::full_min
bool full_min
If true, use the full minimizer.
Definition: interpm_krige.h:698
o2scl::success
@ success
Success.
Definition: err_hnd.h:47
o2scl_inte_qng_coeffs::x1
static const double x1[5]
Definition: inte_qng_gsl.h:48
o2scl::interpm_krige_nn::eval
double eval(const vec2_t &x2, size_t iout, vec_func_t &fcovar) const
Given covariance function fcovar and input vector x return the result of the interpolation for functi...
Definition: interpm_krige.h:1186
o2scl::exc_einval
@ exc_einval
invalid argument supplied by user
Definition: err_hnd.h:59
o2scl_linalg::cholesky_decomp
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
o2scl::interpm_krige_optim::ff1
std::vector< std::function< double(const mat_row_t &, const mat_row_t &)> > ff1
Function objects for the covariance.
Definition: interpm_krige.h:409
o2scl::matrix_row_gen
Generic object which represents a row of a matrix.
Definition: vector.h:2607
o2scl::interpm_krige_optim::ff2
std::vector< std::function< double(const mat_row_t &, const vec_t &)> > ff2
Function objects for the covariance.
Definition: interpm_krige.h:403
o2scl::interpm_krige::verbose
int verbose
Verbosity parameter (default 0)
Definition: interpm_krige.h:111
o2scl::interpm_krige_optim::len_max
double len_max
Maximum for length parameter range.
Definition: interpm_krige.h:424
o2scl::interpm_krige::min
ubvector min
Minimum values for rescaling.
Definition: interpm_krige.h:371
o2scl::interpm_krige::matrix_cholesky
static const size_t matrix_cholesky
Use Cholesky decomposition.
Definition: interpm_krige.h:104
o2scl::interpm_krige_nn::nd_in
size_t nd_in
The number of dimensions of the inputs.
Definition: interpm_krige.h:1368
o2scl::interpm_krige_nn::n_order
size_t n_order
The order of the interpolation (specified by set_data() )
Definition: interpm_krige.h:1364
o2scl::interpm_krige_optim::len
std::vector< double > len
The covariance function length scale for each output function.
Definition: interpm_krige.h:412
o2scl::interpm_krige_nn::eval_jackknife
double eval_jackknife(const vec2_t &x2, size_t iout, func_vec_t &fcovar) const
Compute a quality factor for interpolation using jackknife resampling.
Definition: interpm_krige.h:1258
o2scl::vector_copy_jackknife
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:2379
o2scl::interpm_krige::matrix_LU
static const size_t matrix_LU
Use LU decomposition.
Definition: interpm_krige.h:106
o2scl::interpm_krige_nn::x
std::vector< vec_t > x
A vector of pointers holding the data.
Definition: interpm_krige.h:1372
o2scl::interpm_krige_optim::def_min
min_brent_gsl def_min
Default minimizer.
Definition: interpm_krige.h:695
o2scl::interpm_krige::set_data_noise
int set_data_noise(size_t n_in, size_t n_out, size_t n_points, mat_t &user_x, mat2_t &user_y, func_vec_t &fcovar, const vec_t &noise_var, bool rescale=false, bool err_on_fail=true)
Initialize the data for the interpolation.
Definition: interpm_krige.h:120
o2scl::interpm_krige_optim::len_min
double len_min
Minimum for length parameter range.
Definition: interpm_krige.h:421
o2scl::min_base
One-dimensional minimization [abstract base].
Definition: min.h:42
O2SCL_ERR
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
o2scl::funct
std::function< double(double)> funct
One-dimensional function typedef in src/base/funct.h.
Definition: funct.h:48
o2scl::interpm_krige_optim::qual_fun
double qual_fun(double xlen, double noise_var, size_t iout, vec3_t &y, int &success)
Function to optimize the covariance parameters.
Definition: interpm_krige.h:443
o2scl::min_base::min
virtual int min(double &x, double &fmin, func_t &func)=0
Calculate the minimum min of func w.r.t 'x'.
o2scl::interpm_krige::x
mat_t x
The data.
Definition: interpm_krige.h:367
o2scl::interpm_krige_optim::verbose
int verbose
Verbosity parameter.
Definition: interpm_krige.h:687
o2scl::interpm_krige_optim::set_data
int set_data(size_t n_in, size_t n_out, size_t n_points, mat_t &user_x, mat2_t &user_y, const vec2_t &len_precompute, bool rescale=false, bool err_on_fail=true)
Initialize the data for the interpolation.
Definition: interpm_krige.h:933
o2scl::interpm_krige_nn::verbose
int verbose
Verbosity parameter (default 0)
Definition: interpm_krige.h:974
o2scl_linalg::diagonal_has_zero
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
o2scl::interpm_krige_optim::nlen
size_t nlen
Number of length scale points to try when full minimizer is not used (default 20)
Definition: interpm_krige.h:692
o2scl::interpm_krige_nn::np
size_t np
The number of points.
Definition: interpm_krige.h:1366
o2scl::interpm_krige
Multi-dimensional interpolation by kriging.
Definition: interpm_krige.h:74
o2scl::min_brent_gsl
One-dimensional minimization using Brent's method (GSL)
Definition: min_brent_gsl.h:95
o2scl::interpm_krige_optim::mp
min_base * mp
Pointer to the user-specified minimizer.
Definition: interpm_krige.h:438
o2scl::interpm_krige_optim::set_data_noise
int set_data_noise(size_t n_in, size_t n_out, size_t n_points, mat_t &user_x, mat2_t &user_y, const vec2_t &noise_var, const vec3_t &len_precompute, bool rescale=false, bool err_on_fail=true)
Initialize interpolation routine.
Definition: interpm_krige.h:712
o2scl::interpm_krige_nn::find_lin_indep
void find_lin_indep(const vec2_t &x2, size_t iout, vec_func_t &fcovar, ubvector_size_t &index, ubvector_size_t &indep) const
Find a set of linearly independent points.
Definition: interpm_krige.h:1105
o2scl::interpm_krige_optim::loo_npts
size_t loo_npts
Number of points to test for cross validation (default 100)
Definition: interpm_krige.h:683

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