Akima splines

Recently I was looking for spline interpolation for creating curves from a set of samples. I knew cubic splines which are piecewise cubic polynomials fitted such that they are continuous up to the second derivative. I almost went ahead and implemented cubic splines using a matrix solver but then I found that the fastmath Clojure library already provides splines. The fastmath spline interpolation module is based on the interpolation module of the Java Smile library. I saved the interpolated samples to a text file and plotted them with Gnuplot.

(require '[fastmath.interpolation :as interpolation])
(use '[clojure.java.shell :only [sh]])

(def px [0 1 3 4 5 8 9])
(def py [0 0 1 7 3 4 6])
(spit "/tmp/points.dat" (apply str (map (fn [x y] (str x " " y "\n")) px py)))

(def cspline (interpolation/cubic-spline px py))
(def x (range 0 9.0 0.01))
(spit "/tmp/cspline.dat" (apply str (map (fn [x y] (str x " " y "\n")) x (map cspline x))))
(sh "gnuplot" "-c" "plot.gp" "/tmp/cspline.png" "/tmp/cspline.dat")
(sh "display" "/tmp/cspline.png")

I used the following Gnuplot script plot.gp for plotting:

set terminal pngcairo size 640,480
set output ARG1
set xlabel "x"
set ylabel "y"
plot ARG2 using 1:2 with lines title "spline", "/tmp/points.dat" using 1:2 with points title "points"

I used a lightweight configuration of the fastmath library without MKL and OpenBLAS. See following deps.edn:

{:deps {org.clojure/clojure {:mvn/version "1.12.1"}
        generateme/fastmath {:mvn/version "2.4.0" :exclusions [com.github.haifengl/smile-mkl org.bytedeco/openblas]}}}

The result is shown in the following figure. One can see that the spline is smooth and passes through all points, however it shows a high degree of oscillation:

However I found another spline algorithm in the fastmath wrappers: The Akima spline. The Akima spline needs at least 5 points and it first computes the gradient of the lines connecting the points. Then for each point it uses a weighted average of the previous and next slope value. The slope values are weighted using the absolute difference of the previous two slopes and the next two slopes, i.e. the curvature. The first and last two points use a special formula: The first and last point use the next or previous slope and the second and second last point use an average of the neighbouring slopes.

(require '[fastmath.interpolation :as interpolation])
(use '[clojure.java.shell :only [sh]])

(def px [0 1 3 4 5 8 9])
(def py [0 0 1 7 3 4 6])
(spit "/tmp/points.dat" (apply str (map (fn [x y] (str x " " y "\n")) px py)))

(def aspline (interpolation/akima-spline px py))
(def x (range 0 9.0 0.01))
(spit "/tmp/aspline.dat" (apply str (map (fn [x y] (str x " " y "\n")) x (map aspline x))))
(sh "gnuplot" "-c" "plot.gp" "/tmp/aspline.png" "/tmp/aspline.dat")
(sh "display" "/tmp/aspline.png")

So if you have a data set which causes cubic splines to oscillate, give Akima splines a try!

Enjoy!

Update: u/joinr showed, how you can use Clojupyter to quickly test a lot of splines.