Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
7038 Discussions

## MKL data fitting/interpolation on 12 data points

Beginner
1,336 Views

Hi all,

I'm working on a project related to data fitting, and this is the first time I've ever used MKL, so I'm not super familiar with all the parameters being used for data fitting. So I hope you all could help me with this problem.

So, I current have 12 point values of x and y, and they are paired, e.g (1, 19.313); (24, 16.5697), etc

const double x[] = { 1, 24, 48, 71, 95, 118, 148,
222, 296, 371, 445, 519 };
const double y[] = { 19.31384948, 16.56976694, 9.565388654, 1.697567158,
-4.36924024, -6.412502868, -3.530705482, 10.88522928, 22.44699452,
24.79560359, 24.88258818, 25.38450349 };

So I want to use akima interpolation to have a curve that goes through all the points. That curve should have 519 points. In Matlab, I can simply use the function interp1 like:

curve = interp1(x, y, 1:x(end), 'makima');

I wonder if I can do the same thing for MKL.

For my code, I currently have:

#include <mkl.h>

int main()
{

const double x[] = { 1, 24, 48, 71, 95, 118, 148,
222, 296, 371, 445, 519 };
const double y[] = { 19.31384948, 16.56976694, 9.565388654, 1.697567158,
-4.36924024, -6.412502868, -3.530705482, 10.88522928, 22.44699452,
24.79560359, 24.88258818, 25.38450349 };
const int N = sizeof(x)/sizeof(x[0]);

//
const int NSITE = 1;
int status; /* Status of a Data Fitting operation */

MKL_INT nx = N; /* The size of partition x */
MKL_INT xhint = DF_NON_UNIFORM_PARTITION; /* Additional information about the structure of breakpoints */

MKL_INT ny = 1; /* Function dimension */
/* Parameters describing the spline */

/* Parameters describing interpolation computations */
MKL_INT nsite = NSITE; /* Number of interpolation sites */
double site[NSITE]; /* Array of interpolation sites */
MKL_INT sitehint = DF_SORTED_DATA;

MKL_INT ndorder = 1;
MKL_INT dorder = 1; /* Parameters defining the type of interpolation */

double* datahint = 0;// DF_NO_APRIORI_INFO; /* Additional information on partition and interpolation sites */

double r[NSITE*519]; /* Array of interpolation results */
MKL_INT rhint = DF_MATRIX_STORAGE_ROWS; /* Additional information on the structure of the results */
MKL_INT* cell = 0; //NULL; /* Array of cell indices */

MKL_INT s_order = DF_PP_CUBIC; /* Spline order */
MKL_INT s_type = DF_PP_AKIMA; /* Spline type */

MKL_INT ic_type = DF_NO_IC; ; /* Type of internal conditions */
double* ic = 0; // NULL; /* Array of internal conditions */

MKL_INT bc_type = DF_NO_BC; // DF_BC_NOT_A_KNOT; /* Type of boundary conditions */
double* bc = 0; // NULL; /* Array of boundary conditions */

double scoeff[(N - 1)* DF_PP_CUBIC]; /* Array of spline coefficients */

/* Set spline parameters in the Data Fitting task */
status = dfdEditPPSpline1D(task, s_order, s_type, bc_type, bc, ic_type,
ic, scoeff, scoeffhint);

cout << "status: " << status << "\n";
status = dfdInterpolate1D(task, DF_INTERP, DF_METHOD_PP, nsite, site,
sitehint, ndorder, &dorder, datahint, r, rhint, cell);
cout << "status: " << status << "\n";

I'm super confused about the nsite, site, ny, nx, etc even though I read through the instructions/documentations of these. I'd really appreciate if someone could help guide me through this. Thanks so much!

7 Replies
Moderator
1,328 Views

Did you look at the dfdakimacubicspline.c example from the MKL package? This example shows how to calculate Akima cubic spline...

Beginner
1,317 Views

I did, and that file is not for interpolation but for calculation of the integral..

Moderator
1,313 Views

well, could you point us out what exactly you don't understand from the description of these routines and examples?

Beginner
1,304 Views

I'm confused about the site, nsite, sitehint. It says that site is an array of interpolation sites. What is an interpolation site? I'm not sure if nsite should be 519. The example you recommended doesn't use the function  dfdInterpolate1D, which is what I need for interpolation.

For boundary conditions, what should I choose? what is boundary condition anyway?

It'd be great if you could help me initialize all the necessary parameters for my specific example of 12-point arrays x and y I provided above. Thanks so much

Employee
1,292 Views

Hello Brian,

Interpolation site is the point’s coordinate where you would like to get interpolated value.
nsite – is a number of interpolation sites where interpolated values are needed. So NSITE should be set as 519 in your code.

sitehint – is an additional information about interpolation sites that is used during computation, e.g.:  if you have sorted array of interpolation sites you should set sitehint as DF_SORTED_DATA - this information can increase performance of dfdInterpolate1D routine. If you have no such information about interpolation sites you can set sitehint as DF_NO_HINT.

Regarding examples with dfdInterpolate1D you can refer to dfdquadraticspline.c example (example shows usage of quadratic splines and dfdInterpolate1D routine).

Boundary conditions are used for spline construction and depends on the available information, e.g.: in dfdakimacubicspline.c example we provide 2nd left and 1st right derivatives during spline construction.
Akima splines support next boundary conditions: DF_BC_NOT_A_KNOT, DF_BC_FREE_END, DF_BC_1ST_LEFT_DER, DF_BC_1ST_RIGHT_DER, DF_BC_2ST_LEFT_DER, DF_BC_2ND_RIGHT_DER.

Please note that in case of DF_BC_NOT_A_KNOT and DF_BC_FREE_END conditions no additional information is needed to be set. Math aspects about that can be found here: https://software.intel.com/content/www/us/en/develop/documentation/mkl-developer-reference-c/top/data-fitting-functions/mathematical-conventions-for-data-fitting-functions.html

1. I suggest to set NSITE = 519 and fill the site array with the appropriate values
2. r[NSITE*519] can be changed by r[NSITE]

Best regards,
Pavel

Valued Contributor III
1,272 Views

You are unlikely to generate a curve that exactly fits 12 points unless they are fitted by a lower order curve and in that case you are oversupplied with points.

This is not as simple as it looks.

Moderator
1,216 Views

The issue is closing and we will no longer respond to this thread. If you require additional assistance from Intel, please start a new thread. Any further interaction in this thread will be considered community only.