Higher-dimensional Interpolation

Two-dimensional interpolation

There are two types of two-dimensional interpolation classes, the first is based on a function defined on a two-dimensional grid (though the spacings between grid points need not be equal). The class o2scl::interp2_direct implements bilinear or bicubic interpolation, and is based on D. Zaslavsky's routines at https://github.com/diazona/interp2d (licensed under GPLv3). A slightly slower (but a bit more flexible) alternative is successive use of o2scl::interp_base objects, implemented in o2scl::interp2_seq .

If data is arranged without a grid, then o2scl::interp2_neigh performs nearest-neighbor interpolation. At present, the only way to compute contour lines on data which is not defined on a grid is to use this class or one of the multi-dimensional interpolation classes described below the data on a grid and then use o2scl::contour afterwards.

Multi-dimensional interpolation

Multi-dimensional interpolation for table defined on a grid is possible with o2scl::tensor_grid. See the documentation for o2scl::tensor_grid::interpolate(), o2scl::tensor_grid::interp_linear() and o2scl::tensor_grid::rearrange_and_copy() . Also, if you want to interpolate rank-1 indices to get a vector result, you can use o2scl::tensor_grid::interp_linear_vec() .

If the data is not on a grid, then inverse distance weighted interpolation is performed by o2scl::interpm_idw.

An experimental class for multidimensional-dimensional kriging is also provided in o2scl::interpm_krige .

Interpolation on a rectangular grid

/* Example: ex_interp2.cpp
-------------------------------------------------------------------
A simple example for two-dimensional interpolation using
the interp_twod class.
*/
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <o2scl/interp2_seq.h>
#include <o2scl/test_mgr.h>
using namespace std;
using namespace o2scl;
// A function for filling the data and comparing results
double f(double x, double y) {
return pow(sin(0.1*x+0.3*y),2.0);
}
int main(void) {
cout.setf(ios::scientific);
int i,j;
typedef boost::numeric::ublas::matrix_row<ubmatrix> ubmatrix_row;
// Create the sample data
ubvector x(3), y(3);
ubmatrix data(3,3);
// Set the grid
x[0]=0.0;
x[1]=1.0;
x[2]=2.0;
y[0]=3.0;
y[1]=2.0;
y[2]=1.0;
// Set and print out the data
cout << endl;
cout << "Data: " << endl;
cout << " x | ";
for(i=0;i<3;i++) cout << x[i] << " ";
cout << endl;
cout << " y |" << endl;
cout << "-------------|-";
for(i=0;i<3;i++) cout << "-------------";
cout << endl;
for(i=0;i<3;i++) {
cout << y[i] << " | ";
for(j=0;j<3;j++) {
data(i,j)=f(x[j],y[i]);
cout << data(i,j) << " ";
}
cout << endl;
}
cout << endl;
// Perform the interpolation
cout << "x y Calc. Exact" << endl;
// Interpolation, x-first
double tol=0.05;
double tol2=0.4;
ti.set_data(3,3,x,y,data,itp_cspline);
double x0, y0, x1, y1;
x0=0.5; y0=1.5;
cout << x0 << " " << y0 << " "
<< ti.eval(x0,y0) << " " << f(x0,y0) << endl;
x0=0.99; y0=1.99;
cout << x0 << " " << y0 << " "
<< ti.eval(x0,y0) << " " << f(x0,y0) << endl;
x0=1.0; y0=2.0;
cout << x0 << " " << y0 << " "
<< ti.eval(x0,y0) << " " << f(x0,y0) << endl;
cout << endl;
// Interpolation, y-first
ti.set_data(3,3,x,y,data,itp_cspline);
x0=0.5; y0=1.5;
cout << x0 << " " << y0 << " "
<< ti.eval(x0,y0) << " " << f(x0,y0) << endl;
x0=0.99; y0=1.99;
cout << x0 << " " << y0 << " "
<< ti.eval(x0,y0) << " " << f(x0,y0) << endl;
x0=1.0; y0=2.0;
cout << x0 << " " << y0 << " "
<< ti.eval(x0,y0) << " " << f(x0,y0) << endl;
cout << endl;
t.report();
return 0;
}
// End of example

This example creates a sample 3 by 3 grid of data with the function $ \left[ \sin \left( x/10 + 3 y/10 \right) \right]^2 $ and performs some interpolations and compares them with the exact result.

Data:
x | 0.000000e+00 1.000000e+00 2.000000e+00
y |
-------------|----------------------------------------
3.000000e+00 | 6.136010e-01 7.080734e-01 7.942506e-01
2.000000e+00 | 3.188211e-01 4.150164e-01 5.145998e-01
1.000000e+00 | 8.733219e-02 1.516466e-01 2.298488e-01
x y Calc. Exact
5.000000e-01 1.500000e+00 6.070521e-01 2.298488e-01
9.900000e-01 1.990000e+00 4.187808e-01 4.110774e-01
1.000000e+00 2.000000e+00 4.150164e-01 4.150164e-01
5.000000e-01 1.500000e+00 6.070521e-01 2.298488e-01
9.900000e-01 1.990000e+00 4.187808e-01 4.110774e-01
1.000000e+00 2.000000e+00 4.150164e-01 4.150164e-01
0 tests performed.
All tests passed.

Contour lines

This example generates contour lines of the function

\[ z = f(x,y) = 15 \exp \left[ - \frac{1}{20^2}\left( x-20 \right)^2 - \frac{1}{5^2}\left(y-5\right)^2\right] + 40 \exp \left[ - \frac{1}{500}\left( x-70 \right)^2 - \frac{1}{2^2}\left(y-2\right)^2\right] \]

/* Example: ex_contour.cpp
-------------------------------------------------------------------
Example for generating contour lines
*/
#include <iostream>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <o2scl/test_mgr.h>
#include <o2scl/contour.h>
#include <o2scl/hdf_io.h>
using namespace std;
using namespace o2scl;
using namespace o2scl_hdf;
// A function defining the three-dimensional surface
// for which we want to compute contour levels
double fun(double x, double y) {
return 15.0*exp(-pow(x-20.0,2.0)/400.0-pow(y-5.0,2.0)/25.0)
+40.0*exp(-pow(x-70.0,2.0)/500.0-pow(y-2.0,2.0)/4.0);
}
// A function for outputting the data to cout
int print_data(int nx, int ny, ubvector &x, ubvector &y,
ubmatrix &data);
// A function for printing the contour information to a file
int file_out(string prefix, ubvector &x, ubvector &y,
ubmatrix &data, vector<contour_line> &conts);
int main(void) {
cout.setf(ios::scientific);
contour co;
// Initialize the data
ubvector x(12), y(10);
ubmatrix data(12,10);
for(size_t i=0;i<10;i++) {
y[i]=((double)i);
}
for(size_t i=0;i<12;i++) {
x[i]=((double)i)*((double)i);
}
for(size_t j=0;j<12;j++) {
for(size_t k=0;k<10;k++) {
data(j,k)=fun(x[j],y[k]);
}
}
co.set_data(12,10,x,y,data);
// Print out the data
print_data(12,10,x,y,data);
// Set the contour levels
ubvector levels(7);
levels[0]=5.0;
levels[1]=10.0;
levels[2]=0.002;
levels[3]=20.0;
levels[4]=0.2;
levels[5]=30.0;
levels[6]=2.0;
co.set_levels(7,levels);
// Compute the contours
vector<contour_line> conts;
co.calc_contours(conts);
// Output to a file
file_out("c1",x,y,data,conts);
// Print the contours to the screen and test to make sure
// that they match the requested level
size_t nc=conts.size();
for(size_t i=0;i<nc;i++) {
cout << "Contour " << i << " at level " << conts[i].level << ":" << endl;
size_t cs=conts[i].x.size();
for(size_t j=0;j<cs;j++) {
cout << "(" << conts[i].x[j] << ", " << conts[i].y[j] << ") "
<< fun(conts[i].x[j],conts[i].y[j]) << endl;
t.test_rel(fun(conts[i].x[j],conts[i].y[j]),conts[i].level,
1.0,"curve");
}
cout << endl;
}
// Refine the data using cubic spline interpolation
co.regrid_data(5,5);
ubvector *y2;
ubmatrix *data2;
size_t sx, sy;
co.get_data(sx,sy,x2,y2,data2);
// Recompute the contours
vector<contour_line> conts2;
co.calc_contours(conts2);
// Output to a file
file_out("c2",*x2,*y2,*data2,conts2);
t.report();
return 0;
}
// End of example

The figure below shows contour lines in the region $ x\in(0,121), y\in(0,9) $. The data grid is represented by plus signs, and the associated generated contours. The figure clearly shows the peaks at $ (20,5) $ and $ (70,2) $ .

Contour example plot

The o2scl::contour class can also use interpolation to attempt to refine the data grid. The new contours after a refinement of a factor of 5 is given in the figure below.

Contours after regrid_data()
boost::numeric::ublas::matrix< double >
o2scl::interp2_seq
Two-dimensional interpolation class by successive one-dimensional interpolation.
Definition: interp2_seq.h:82
o2scl::contour::set_levels
void set_levels(size_t nlevels, vec_t &ulevels)
Set the contour levels.
Definition: contour.h:337
boost::numeric::ublas::vector< double >
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::test_mgr::report
bool report() const
Provide a report of all tests so far.
o2scl::contour::get_data
void get_data(size_t &sizex, size_t &sizey, ubvector *&x_fun, ubvector *&y_fun, ubmatrix *&udata)
Get the data.
Definition: contour.h:387
o2scl::contour
Calculate contour lines from a two-dimensional data set.
Definition: contour.h:240
o2scl::contour::regrid_data
void regrid_data(size_t xfact, size_t yfact, size_t interp_type=o2scl::itp_cspline)
Regrid the data.
o2scl_inte_qng_coeffs::x2
static const double x2[5]
Definition: inte_qng_gsl.h:66
o2scl::test_mgr::test_rel
bool test_rel(data_t result, data_t expected, data_t rel_error, std::string description)
Test for .
Definition: test_mgr.h:168
o2scl_inte_qng_coeffs::x1
static const double x1[5]
Definition: inte_qng_gsl.h:48
o2scl::test_mgr
A class to manage testing and record success and failure.
Definition: test_mgr.h:50
o2scl_hdf
The O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl namespace ...
Definition: table.h:59
o2scl::contour::calc_contours
void calc_contours(std::vector< contour_line > &clines)
Calculate the contours.
o2scl::interp2_seq::set_data
void set_data(size_t n_x, size_t n_y, vec_t &x_grid, vec_t &y_grid, mat_t &data, size_t interp_type=itp_cspline)
Initialize the data for the 2-dimensional interpolation.
Definition: interp2_seq.h:115
o2scl::test_mgr::set_output_level
void set_output_level(int l)
Set the output level.
Definition: test_mgr.h:119
o2scl::interp2_seq::eval
double eval(double x, double y) const
Perform the 2-d interpolation.
Definition: interp2_seq.h:181
o2scl::itp_cspline
@ itp_cspline
Cubic spline for natural boundary conditions.
Definition: interp.h:73
o2scl::contour::set_data
void set_data(size_t sizex, size_t sizey, const vec_t &x_fun, const vec_t &y_fun, const mat_t &udata)
Set the data.
Definition: contour.h:266

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