Evaluation of custom bspline using pchip and makima

[stefphd]

I am trying to define a custom bspline using Function.bspline from a pchip spline obtained by a MATLAB b-form. When I compare the evaluation of Function.bspline to the one obtained from MATLAB, I see spikes at the grid points in the casadi evaluation. Here the code to reproduce the issue:

import casadi.*

% vals
x = 0 : 1 : 10; %grid (spacing 1)
f = @(x) cos(x).^2.*sin(x); % just an example

% matlab spline
pp = interp1(x, f(x), 'pchip', 'pp'); %pp form using pchip (issue)
% pp = makima(x, f(x)); %pp form using makima (issue)
% pp = interp1(x, f(x), 'spline', 'pp'); %pp form using spline (no issue)
bform = fn2fm(pp, 'B-'); %convert pp-form to B-form to obtain knots and coefs

% casadi spline
knots = {bform.knots};
coefs = bform.coefs;
degree = bform.order-1;
mySpline = Function.bspline('myspline', knots, coefs, degree);

% compare evals
x0 = linspace(0,10,1001); 
y0 = ppval(pp, x0); % pp eval
y1 = fnval(bform, x0); % bform eval
y2 = full(mySpline(x0)); % casadi eval
plot(x0,y0, '-',x0,y1, '--', x0,y2, ':')
legend('pp-form',' b-form','casadi bspline')

The same issues occurs when using e.g. pp = makima(x,f(x)) instead of pchip, while the results are identical when using pp = interp1(x, f(x), 'spline', 'pp').

The issue seems to be related to the repeating internal knots. Indeed, when spacing the knots by a small quantity (say 1e3*eps) before creating the casadi bspline, the issue is solved with both pchip and makima splines. The workaround is

% space internal repeating knots
delta = 1e3*eps;
first_knots = bform.knots == bform.knots(1);
end_knots = bform.knots == bform.knots(end);
internal_knots = ~first_knots & ~end_knots;
k = 1; % init counter
while k <= numel(bform.knots) % loop knots
    if ~internal_knots(k) %skip external knots
        k = k + 1;
        continue
    end
    % find repeating times
    repeating_times = sum(bform.knots(k) == bform.knots);
    % space repeating knots
    for l = k+1 : k+repeating_times-1
        bform.knots(l) = bform.knots(l-1) + delta;
    end
    k = k+repeating_times; % next knots
end