Interpolation

O2scl

Interpolation contents

• Interpolation introduction

• Lookup and binary search

Interpolation introduction

Basic interpolation of generic vector types is performed by interp_vec and its children. The vector representing the independent variable must be monotonic, but need not be equally-spaced. The difference between the two classes is analogous to the difference between using gsl_interp_eval() and gsl_spline_eval() in GSL. You can create a interp_vec object and use it to interpolate among any pair of chosen vectors. For example, cubic spline interpolation with natural boundary conditions:

boost::numeric::ublas::vector<double> x(20), y(20);
// fill x and y with data
o2scl::interp<> oi(itp_cspline);
double y_half=oi.eval(0.5,20,x,y);

Alternatively, you can create a interp_vec object which can be optimized for a pair of vectors that you specify in advance (now using linear interpolation instead):

boost::numeric::ublas::vector<double> x(20), y(20);
// fill x and y with data
o2scl::interp_vec<> oi(20,x,y,itp_linear);
double y_half=oi.eval(0.5);

These interpolation classes require that the vector x is either monotonically increasing or monotonically decreasing, but these classes do not verify that this is the case. If the vectors or not monotonic, then the interpolation functions are not defined. These classes give identical results to the GSL interpolation routines when the vector is monotonically increasing.

These interpolation classes will extrapolate to regions outside those defined by the user-specified vectors and will not warn you when they do this (this is not the same behavior as in GSL).

One-dimensional Gaussian process interpolation (i.e. Kriging) is also provided in interp_krige for a generic user-specified covariance function and interp_krige_optim which allows one to optimize the parameters to fit the data.

The different interpolation types are defined in src/base/interp.h

enumerator itp_linear

Linear.

enumerator itp_cspline

Cubic spline for natural boundary conditions.

enumerator itp_cspline_peri

Cubic spline for periodic boundary conditions.

enumerator itp_akima

Akima spline for natural boundary conditions.

enumerator itp_akima_peri

Akima spline for periodic boundary conditions.

enumerator itp_monotonic

Monotonicity-preserving interpolation (unfinished)

enumerator itp_steffen

Steffen’s monotonic method.

enumerator itp_nearest_neigh

Nearest-neighbor lookup.

enumerator itp_gp_rbf_noise_loo_cv

Gaussian process with leave-one out cross validation (experimental)

enumerator itp_gp_rbf_noise_max_lml

Gaussian process with maximum log marginal likelihood (experimental)

Integrals are always computed assuming that if the limits are ordered so that if the upper limit appears earlier in the array x in comparison to the lower limit, that the value of the integral has the opposite sign than if the upper limit appears later in the array x.

The class interp_vec is based on the lower-level interpolation classes of type interp_base. Also, the interpolation classes based on interp_base and also the class interp_vec also have defined a function operator() which also returns the result of the interpolation.

Lookup and binary search

The classes search_vec and search_vec_ext contain searching functions for generic vector types which contain monotonic (either increasing or decreasing) data. It is search_vec which is used internally by the interpolation classes to perform cached binary searching. These classes also allow one to to exhaustively search for the index of an element in a vector without specifying in advance if the vector is increasing or decreasing, e.g. o2scl::search_vec::ordered_lookup().

Interpolation example

/* Example: ex_interp.cpp
   -------------------------------------------------------------------
   A simple example for interpolation. See "License Information" 
   section of the documentation for license information.
*/

#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>

#include <o2scl/interp.h>
#include <o2scl/interp_krige.h>
#include <o2scl/test_mgr.h>
#include <o2scl/hdf_io.h>
#include <o2scl/hdf_file.h>

using namespace std;
using namespace o2scl;
using namespace o2scl_hdf;

typedef boost::numeric::ublas::vector<double> ubvector;

// A function for filling the data and comparing results
double f(double x, const double &mean, const double &sd) {
  return (sin(1.0/(0.3+x))-mean)/sd;
}

int main(void) {
  cout.setf(ios::scientific);

  test_mgr t;
  t.set_output_level(1);

  // Create the sample data

  static const size_t N=20;
  ubvector x(N), y(N);
  x[0]=0.0;
  y[0]=f(x[0],0.0,1.0);
  for(size_t i=1;i<N;i++) {
    x[i]=x[i-1]+pow(((double)i)/40.0,2.0);
    y[i]=f(x[i],0.0,1.0);
  }
  
  table<> tdata;
  tdata.line_of_names("x y");

  double y_mean=vector_mean(y);
  for(size_t i=0;i<N;i++) {
    y[i]-=y_mean;
  }
  double y_sd=vector_stddev(y);
  cout << "Old mean and std. dev.: " << y_mean << " " << y_sd << endl;
  for(size_t i=0;i<N;i++) {
    y[i]/=y_sd;
    double line[2]={x[i],y[i]};
    tdata.line_of_data(2,line);
  }
  double y_mean2=vector_mean(y);
  double y_sd2=vector_stddev(y);
  cout << "New mean and std. dev.: " << y_mean2 << " " << y_sd2 << endl;
  cout << endl;

  interp_vec<ubvector> iv_lin(N,x,y,itp_linear);
  interp_vec<ubvector> iv_csp(N,x,y,itp_cspline);
  interp_vec<ubvector> iv_aki(N,x,y,itp_akima);
  interp_vec<ubvector> iv_mon(N,x,y,itp_monotonic);
  interp_vec<ubvector> iv_stef(N,x,y,itp_steffen);

  covar_funct_rbf_noise cfrn;

  vector<double> p={1.0e-4,1.0};
  vector<double> p2={-9.0,-15.0};
  vector<vector<double>> param_lists;
  param_lists.push_back(p);
  param_lists.push_back(p2);
  
  interp_krige_optim<ubvector,ubvector,covar_funct_rbf_noise> iko;
  iko.verbose=2;
  //iko.def_mmin.verbose=2;
  iko.full_min=true;
  iko.mode=iko.mode_loo_cv;
  iko.set_covar_optim(cfrn,param_lists);
  iko.set(N,x,y);

  interp_krige_optim<ubvector,ubvector,covar_funct_rbf_noise> iko2;
  iko2.verbose=2;
  //iko2.def_mmin.verbose=2;
  iko2.full_min=true;
  iko2.mode=iko2.mode_max_lml;
  iko2.set_covar_optim(cfrn,param_lists);
  iko2.set(N,x,y);
  cout << endl;

  double max=x[x.size()-1];

  cout << "\nx              y             iko           iko2: " << endl;
  for(size_t i=0;i<N;i++) {
    cout.setf(ios::showpos);
    cout << x[i] << " " << f(x[i],y_mean,y_sd) << " "
	 << iko.eval(x[i]) << " "
	 << iko2.eval(x[i]) << endl;
    cout.unsetf(ios::showpos);
  }
  cout << endl;

  cout << "\nx              y             iko           iko2: " << endl;
  size_t N2=N*100;
  table<> tresult;
  tresult.line_of_names("x y ylin ycsp yaki ymon ystef yiko yiko_lml"); 
  for(size_t i=0;i<=N2;i++) {
    double x9=((double)i)/((double)N2)*max;
    double line[9]={x9,f(x9,y_mean,y_sd),iv_lin.eval(x9),iv_csp.eval(x9),
		    iv_aki.eval(x9),iv_mon.eval(x9),iv_stef.eval(x9),
		    iko.eval(x9),iko2.eval(x9)};
    tresult.line_of_data(9,line);
    if (i%50==0) {
      cout.setf(ios::showpos);
      cout << x9 << " " << f(x9,y_mean,y_sd) << " " << iko.eval(x9) << " "
	   << iko2.eval(x9) << endl;
      cout.unsetf(ios::showpos);
    }
  }

  hdf_file hf;
  hf.open_or_create("ex_interp.o2");
  hdf_output(hf,tdata,"tdata");
  hdf_output(hf,tresult,"tresult");
  hf.close();
  
  t.report();
  return 0;
}
// End of example

Todo

Fix the interpolation plot for this example.

Two and higher-dimensional interpolation

Support for multi-dimensional interpolation is documented in Higher-dimensional Interpolation.

Inverse interpolation and related functions

The equivalent to “inverse” linear interpolation, which computes all the abcissae which have a fixed value of the ordinate, is implemented in the template function o2scl::vector_find_level(). This function together with o2scl::vector_invert_enclosed_sum() can be used to determine confidence limits surrounding the peak of a 1-dimensional data set using linear interpolation. To count level crossings in a function, use o2scl::vector_level_count(). The function o2scl::vector_integ_interp() uses interpolation to compute the integral defined by a set of vectors, and the function o2scl::vector_region_fracint() finds the set of regions which gives a fraction of the integral reported by o2scl::vector_integ_interp().

Derivatives and integrals on a fixed grid

If the indepedent variable is represented by a uniform (equally-spaced) then the functions in src/deriv/vector_derint.h can provide faster (and occasionally more accurate) results. See:

• o2scl::vector_deriv_fivept()

• o2scl::vector_deriv_fivept_tap()

• o2scl::vector_deriv_interp()

• o2scl::vector_deriv_threept()

• o2scl::vector_deriv_threept_tap()

• o2scl::vector_integ_durand()

• o2scl::vector_integ_extended4()

• o2scl::vector_integ_extended8()

• o2scl::vector_integ_threept()

• o2scl::vector_integ_trap()