This module provides cubic smoothing spline for univariate/multivariate/gridded data approximation. The smoothing parameter can be calculated automatically or it can be set manually.

The smoothing parameter should be in range [0, 1] where bounds are:

• 0: The smoothing spline is the least-squares straight line fit to the data
• 1: The natural cubic spline interpolant

The calculation of the smoothing spline requires the solution of a linear system whose coefficient matrix has the form p*A + (1 - p)*B, with the matrices A and B depending on the data sites X. The automatically computed smoothing parameter makes p*trace(A) equal (1 - p)*trace(B).

Python 3.5 or above is supported.

pip install csaps

The module depends only on NumPy and SciPy.

On Windows we highly recommend to use unofficial builds NumPy+MKL and SciPy from Christoph Gohlke.

> mkdir depends

download numpy (https://www.lfd.uci.edu/~gohlke/pythonlibs/#numpy) to "depends" directory
download scipy (https://www.lfd.uci.edu/~gohlke/pythonlibs/#scipy) to "depends" directory

> pip install --find-links depends/ csaps

Univariate data are two vectors: X and Y with the same size. X is data sites, Y is data values. For univariate case X-values must satisfy the condition: x1 < x2 < ... < xN.

You can use UnivariateCubicSmoothingSpline class for uivariate data smoothing.

import numpy as np
import matplotlib.pyplot as plt
import csaps

np.random.seed(1234)

x = np.linspace(-5., 5., 25)
y = np.exp(-(x/2.5)**2) + (np.random.rand(25) - 0.2) * 0.3

sp = csaps.UnivariateCubicSmoothingSpline(x, y, smooth=0.85)

xs = np.linspace(x[0], x[-1], 150)
ys = sp(xs)

plt.plot(x, y, 'o', xs, ys, '-')
plt.show()

The algorithm supports weighting. You can set weights vector that will determine weights for all data values:

import csaps

x = [1., 2., 4., 6.]
y = [2., 4., 5., 7.]
w = [0.5, 1, 0.7, 1.2]

sp = csaps.UnivariateCubicSmoothingSpline(x, y, w)
...

The algorithm supports vectorization. You can compute smoothing splines for X, Y data where X is data site vector and Y is ND-array of data value vectors. The shape of Y array must be: (d0, d1, ..., dN) where dN must equal of X vector size.

In this case the smoothing splines will be computed for all Y data vectors at a time.

For example:

import numpy as np
import csaps

# data sites
x = [1, 2, 3, 4]

# two data vectors
y = np.array([(2, 4, 6, 8), 
              (1, 3, 5, 7)])

sp = csaps.UnivariateCubicSmoothingSpline(x, y)

xi = np.linspace(1, 4, 10)
yi = sp(xi)

print(yi.shape)  # (2, 10)
assert yi.shape[:-1] == y.shape[:-1]
assert yi.shape[-1] == xi.size

Important: The same weights vector and the same smoothing parameter will be used for all Y data.

We can easily smooth multivariate data using univariate smoothing spline, vectorization and parametrization. We should define parametric data sites vector t with condition t1 < t2 < ... < tN and make spline for each dimension X(t), Y(t), ..., M(t).

The module provides MultivariateCubicSmoothingSpline class for this case.

This class is a simple wrapper and it just automatically computes t vector by default and calls UnivariateCubicSmoothingSpline:

# Construct multivariate spline from t and X, Y, Z, ..., M
sx = UnivariateCubicSmoothingSpline(t, data[x])
sy = UnivariateCubicSmoothingSpline(t, data[y])
sz = UnivariateCubicSmoothingSpline(t, data[z])
...
sm = UnivariateCubicSmoothingSpline(t, data[m])

# Or the same with using vectorization
sxyz_m = UnivariateCubicSmoothingSpline(t, data)

In this example we make parametric 3D curve from theta parameter and compute spline from the same theta parameter. In fact we could use UnivariateCubicSmoothingSpline with vectorization instead of MultivariateCubicSmoothingSpline in this case.

import numpy as np

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

import csaps

n = 100
theta = np.linspace(-4 * np.pi, 4 * np.pi, n)
z = np.linspace(-2, 2, n)
r = z ** 2 + 1
np.random.seed(1234)
x = r * np.sin(theta) + np.random.randn(n) * 0.3
np.random.seed(5678)
y = r * np.cos(theta) + np.random.randn(n) * 0.3

data = np.vstack((x, y, z))

sp_theta = csaps.MultivariateCubicSmoothingSpline(data, theta, smooth=0.95)

# or the same
# sp_theta = csaps.UnivariateCubicSmoothingSpline(theta, data, smooth=0.95)

theta_i = np.linspace(-4 * np.pi, 4 * np.pi, 250)
data_i = sp_theta(theta_i)

xi = data_i[0, :]
yi = data_i[1, :]
zi = data_i[2, :]

fig = plt.figure()
ax = fig.gca(projection='3d')
plt.show()

The algorithm can make smoothing splines for ND-gridded data approximation. In this case we use coordinatewise smoothing (tensor-product of univariate splines coefficients).

You can use NdGridCubicSmoothingSpline class for ND-gridded data smoothing. You also can set weights and smoothing parameters for each dimension.

Currently the implementation does not support vectorization for ND-gridded data.

import numpy as np

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

import csaps

xdata = [np.linspace(-3, 3, 61), np.linspace(-3.5, 3.5, 51)]
i, j = np.meshgrid(*xdata, indexing='ij')

ydata = (3 * (1 - j)**2. * np.exp(-(j**2) - (i + 1)**2)
         - 10 * (j / 5 - j**3 - i**5) * np.exp(-j**2 - i**2)
         - 1 / 3 * np.exp(-(j + 1)**2 - i**2))

np.random.seed(12345)
noisy = ydata + (np.random.randn(*ydata.shape) * 0.75)

sp = csaps.NdGridCubicSmoothingSpline(xdata, noisy, smooth=0.988)
ysmth = sp(xdata)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

ax.plot_wireframe(j, i, noisy, linewidths=0.5, color='r')
ax.scatter(j, i, noisy, s=5, c='r')

ax.plot_surface(j, i, ysmth, linewidth=0, alpha=1.0)

plt.show()

Please look through csaps.ipynb file for more examples.

We use pytest and tox (on Travis CI) for testing. Please see test_csaps.py file.

csaps is a Python modified port of MATLAB CSAPS function that is an implementation of Fortran routine SMOOTH from PGS (originally written by Carl de Boor).

csaps-cpp C++11 Eigen based implementation of the algorithm.

C. de Boor, A Practical Guide to Splines, Springer-Verlag, 1978.