Documentation/Examples of casadi.interpolant function in Matlab?

[Anton Savchenko]

Hello guys,

I have been trying to implement the splines with CasADi in Matlab, but haven't been successful thus far.

I think currently the easiest option would be to try the interpolant function, however, I do not understand the required data format - all the examples I found were for python, where list of lists is all you need, however, Matlab keeps throwing this exception at me:

Error using casadi.interpolant (line 78)

No matching function for overload function 'interpolant'.  Prototype:

    INTERPOLANT(char,char,[[double]],[double],struct)

  You have: char, char, [double], double

I tried all I could think of - cell arrays, arrays of tuples, arrays of cell arrays... - nothing seems to match the [[double]] notation. Could you please help - say, how do I call the interpolant function on the following toy example?

import casadi.*

X=[0:.1:1] % first dimension

Y=[0:.2:1] % second dimension

V = X'*Y % grid of values

F = interpolant('F','linear',?, ??)

And while I'm here, what would be the syntax for using b-splines of order 3? Do I have to additionally provide the jacobians, or will the code automatically generate it?

Thanks in advance

Anton

[Joris Gillis]

Right, our docs are long overdue on lookup tables.
I'll get it ready for coming release.


In the meantime, here you have some example code:

import casadi.*


%% 2.1
xgrid = -5:1:5;
ygrid = -4:1:4;

[X1,X2] = ndgrid(xgrid, ygrid);
R = sqrt(X1.^2 + X2.^2)+ eps;
V = sin(R)./(R);

interpn(X1,X2,V,0.5,1,'spline')

LUT = interpolant('LUT','bspline',{xgrid, ygrid},V(:));
LUT([0.5 1])

%% 2.2
xgrid = 1:6;
V = [-1 -1 -2 -3 0 2];
xs = linspace(1,6,1000);

interpn(xgrid,V,0.5,'spline')

LUT = interpolant('LUT','bspline',{xgrid},V);
LUT(0.5)

figure
hold on
plot(xgrid, V, 'o');
res_matlab = interpn(xgrid,V,xs,'spline');
plot(xs, res_matlab)

LUTmap = LUT.map(1000);
res_casadi = full(LUTmap(xs));
plot(xs, res_casadi)

figure()
plot(xs, res_casadi-res_matlab)
%% 2.3

x = MX.sym('x');

y = LUT(x);

solver = nlpsol('solver','ipopt',struct('x',x,'f',y));

sol = solver('x0',2,'lbx',1,'ubx',6);

sol.x


%% 2.4

y = LUT(x);

dy = jacobian(y,x);
ddy = jacobian(dy,x);
dddy = jacobian(ddy,x);

res = Function('f',{x},{y,dy,ddy,dddy});
resmap = res.map(1000);

[yv, dyv, ddyv, dddyv] = resmap(xs);



figure
hold on
plot(xs,full(yv))
plot(xs,full(dyv))
plot(xs,full(ddyv))
plot(xs,full(dddyv))


%%

[Anton Savchenko]

Hello Joris,

thanks for the fast reply - I managed to make it work! The mapping functionality is definitely impressive, like a 100x performance increase compared to for-loops :)

Just one question - are there some options available to modify the bspline basis functions? Like, the degree and amount of employed grid points?

Thanks again

Anton

[Joris Gillis]

There is the "degree" option which controls the degree in each dimension.
The knots of the bspline coincide withe the supplied grid, so just pick the grid like you wish.
Note that interpolant does a fit through your data. If you want to specify knot locations and degree yourself, you may use the underlying Function.bspline

Joris

[andreas...@yahoo.de]

Hi Joris,

could you provide an example of how exactly you pass the degree option? I want to create a bspline of degree 5.

Thanks

Andy