Is there a function for not-a-knot cubic splines?

I need to use the "not-a-knot" cubic spline for interpolation in my R scripts.

Although there are some R packages for splines, none of them seem to consider the "not-a-knot" type, even if it is said to be a rather "popular" type of cubic spline, and it is available in Matlab.

I fear that there is another name for "not-a-knot" cubic splines. It is a cubic spline where the two extra conditions are about the third derivative continuity in the second and before-last knots (instead of fixing the first derivatives at the endpoint knots, as natural cubic splines do, or other choices).

By digging a little bit, it seems that this is implemented in pracma::interp1

From ?pracma::interp1

Interpolation to find ‘yi’, the values of the underlying function at the points in the vector ‘xi’. ... Method `spline' uses the spline approach by Moler et al., and is identical with the Matlab option of the same name, but slightly different from R's spline function.

This doesn't mention "not-a-knot", but I got there by using sos::findFn("not-a-knot"), which brought me to gsignal::pulstran(). That function's "spline" method is described as "Spline interpolation using not-a-knot end conditions"; the method= argument is passed directly to pracma::interp1.

Here is an example comparing:

• two base-R spline implementations
• pracma::interp1
• pracma:::.ml.spline (doesn't do any argument-checking, sorting, duplicate removal ... but does allow extrapolation. I don't know why pracma::interp1 makes this restriction. The package is developed on r-forge but I don't see mailing lists etc., you could contact the maintainer or hack interp1 yourself ...)
• Octave (adapted from here: see code below. There is an RcppOctave package, but it's archived and I couldn't easily get the package to install via remotes::install_version("RcppOctave", version = "0.18-1") - too much configuration drift ...)

library(splines)
library(pracma)
x <- 1:6
y <- c(16, 18, 21, 17, 15, 12)
xi0 <- seq(1, 6, by = 0.1)
xi1 <- seq(1, 7, by = 0.1)
ys1 <- interp1(x, y, xi0, method = "spline")
ys2 <- predict(interpSpline(x, y), xi1)$y
ys3 <- spline(x, y, xout  = xi1)$y
ys4 <- pracma:::.ml.spline(x, y, xi1)
ys5 <- scan("octave_out.mat", comment = "#") ## see below

png("spline2.png")
par(las=1, bty = "l", cex = 1.5, lwd = 2)
plot(x,y, xlim = range(xi1), ylim = range(ys4),
     pch = 16, col = "gray")
cvec <- c(palette()[c(1,2,4,6)],
  ## make last color semi-transparent so we can see the overlay
          adjustcolor(palette()[5], alpha.f = 0.4))
matlines(xi1, cbind(ys2,ys3,ys4, ys5) , col = cvec[1:4], lwd = 2)
lines(xi0, ys1, lwd = 3, col = cvec[5])
legend("bottomleft",
       lty = c(1:4, 1),
       col = cvec,
       legend = c("interpSpline", "spline", ".ml.spine", "octave", "interp1"))
dev.off()

• Octave and .ml.spline output agree
• interp1 agrees within the range of the data (no extrapolation)
• the two R splines give different, but reasonable, answers

Octave code (run in a separate shell)

x = 1:6;
y = [16, 18, 21, 17, 15, 12];
xi1 = linspace(1, 7, 61);
pp = interp1(x, y, 'spline', 'pp');
yy = ppval(pp, xi1);
save octave_out.mat yy