Cubic spline interpolation with PyDynamic.uncertainty.interpolate.interp1d_unc

Interpolate a non-equidistant sine signal using cubic / bspline method with uncertainty propagation. Comparing the resulting uncertainties to a Monte-Carlo experiment yields good overlap.

import numpy as np
import plotly.graph_objs as go
from PyDynamic.uncertainty.interpolate import interp1d_unc
from scipy.interpolate import interp1d

Create non-equidistant sine-signal

n_nodes = 10
t = np.cumsum(range(1, n_nodes))
x = np.sin(t)
ux = np.full_like(x, 0.2)
t, x, ux

(array([ 1,  3,  6, 10, 15, 21, 28, 36, 45]),
 array([ 0.84147098,  0.14112001, -0.2794155 , -0.54402111,  0.65028784,
         0.83665564,  0.27090579, -0.99177885,  0.85090352]),
 array([0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2]))

Interpolate with PyDynamic

ti = np.linspace(np.min(t), np.max(t), 60)
ti, xi, uxi = interp1d_unc(ti, t, x, ux, kind="cubic")
ti, xi, uxi

(array([ 1.        ,  1.74576271,  2.49152542,  3.23728814,  3.98305085,
         4.72881356,  5.47457627,  6.22033898,  6.96610169,  7.71186441,
         8.45762712,  9.20338983,  9.94915254, 10.69491525, 11.44067797,
        12.18644068, 12.93220339, 13.6779661 , 14.42372881, 15.16949153,
        15.91525424, 16.66101695, 17.40677966, 18.15254237, 18.89830508,
        19.6440678 , 20.38983051, 21.13559322, 21.88135593, 22.62711864,
        23.37288136, 24.11864407, 24.86440678, 25.61016949, 26.3559322 ,
        27.10169492, 27.84745763, 28.59322034, 29.33898305, 30.08474576,
        30.83050847, 31.57627119, 32.3220339 , 33.06779661, 33.81355932,
        34.55932203, 35.30508475, 36.05084746, 36.79661017, 37.54237288,
        38.28813559, 39.03389831, 39.77966102, 40.52542373, 41.27118644,
        42.01694915, 42.76271186, 43.50847458, 44.25423729, 45.        ]),
 array([ 0.84147098,  0.5134559 ,  0.26887907,  0.09049061, -0.03895932,
        -0.13672062, -0.22004314, -0.30603185, -0.40015754, -0.48792672,
        -0.55288113, -0.57856251, -0.54851261, -0.45094022, -0.29697194,
        -0.10475754,  0.10755139,  0.32180327,  0.51984651,  0.68357697,
         0.80217466,  0.87988176,  0.92280604,  0.93705524,  0.92873714,
         0.90395947,  0.86883   ,  0.82944871,  0.78981465,  0.74899325,
         0.70534253,  0.65722049,  0.60298513,  0.54099446,  0.4696065 ,
         0.38717924,  0.29207069,  0.18306475,  0.06213925, -0.0672781 ,
        -0.20175235, -0.33784857, -0.4721318 , -0.6011671 , -0.72151951,
        -0.8297541 , -0.92243592, -0.99613001, -1.04740144, -1.07281526,
        -1.06893652, -1.03233027, -0.95956157, -0.84719547, -0.69179702,
        -0.48993128, -0.23816331,  0.06694186,  0.42881915,  0.85090352]),
 array([0.2       , 0.16837861, 0.19741966, 0.1962444 , 0.17596306,
        0.16866971, 0.18622166, 0.20285793, 0.19817887, 0.1825655 ,
        0.17473489, 0.18389472, 0.1992225 , 0.20369264, 0.19571168,
        0.18332587, 0.17637735, 0.1804066 , 0.19195061, 0.20145338,
        0.20182147, 0.19446415, 0.18421328, 0.1767848 , 0.17639785,
        0.18326204, 0.19338735, 0.20091735, 0.20172299, 0.19643529,
        0.18769373, 0.17895702, 0.17368016, 0.17410604, 0.1801646 ,
        0.1894922 , 0.1985253 , 0.20383972, 0.20425977, 0.20029238,
        0.1931453 , 0.18453319, 0.17654847, 0.17142144, 0.17103588,
        0.17629981, 0.18680117, 0.20105082, 0.21705456, 0.23277969,
        0.24637648, 0.25624082, 0.26102733, 0.25968413, 0.25156834,
        0.23674281, 0.21668996, 0.19594729, 0.18494793, 0.2       ]))

Interpolate with Monte Carlo

X_mc = []
for i in range(2000):
    interp_x = interp1d(t, x + ux * np.random.randn(len(x)), kind="cubic")
    xm = interp_x(ti)
    X_mc.append(xm)
x_mc = np.mean(X_mc, axis=0)
ux_mc = np.std(X_mc, axis=0);

Compare results

Visual

# NBVAL_IGNORE_OUTPUT
# interpolated signal
fig = go.Figure(
    [
        go.Scatter(
            name="PyDynamic interpolation",
            x=ti,
            y=xi,
            mode="lines",
            marker_color="rgba(168, 68, 68, 0.8)",
        ),
        go.Scatter(
            hoverinfo="x",
            x=ti,
            y=xi + uxi,
            mode="lines",
            line=dict(width=0),
            showlegend=False,
        ),
        go.Scatter(
            hoverinfo="x",
            x=ti,
            y=xi - uxi,
            line=dict(width=0),
            mode="lines",
            fillcolor="rgba(168, 68, 68, 0.3)",
            fill="tonexty",
            showlegend=False,
        ),
        go.Scatter(
            name="Monte Carlo interpolation",
            x=ti,
            y=x_mc,
            mode="lines",
            marker_color="rgba(68, 68, 168, 0.8)",
        ),
        go.Scatter(
            hoverinfo="x",
            x=ti,
            y=x_mc + ux_mc,
            mode="lines",
            line=dict(width=0),
            showlegend=False,
        ),
        go.Scatter(
            hoverinfo="x",
            x=ti,
            y=x_mc - ux_mc,
            line=dict(width=0),
            mode="lines",
            fillcolor="rgba(68, 68, 168, 0.3)",
            fill="tonexty",
            showlegend=False,
        ),
        go.Scatter(
            name="original",
            x=t,
            y=x,
            mode="lines",
            marker_color="rgba(68, 168, 68, 0.8)",
        ),
        go.Scatter(
            hoverinfo="x",
            x=t,
            y=x + ux,
            mode="lines",
            line=dict(width=0),
            showlegend=False,
        ),
        go.Scatter(
            hoverinfo="x",
            x=t,
            y=x - ux,
            line=dict(width=0),
            mode="lines",
            fillcolor="rgba(68, 168, 68, 0.3)",
            fill="tonexty",
            showlegend=False,
        ),
    ]
)
fig.update_layout(
    yaxis_title="Amplitude (a.u.)",
    xaxis_title="Time (a.u.)",
    title="Comparison between MonteCarlo simulation and PyDynamic interpolation",
    hovermode="x",
)
fig.show()

Numerical

tolerance = 1e-1

np.allclose(xi, x_mc, atol=tolerance, rtol=tolerance)

True

np.allclose(uxi, ux_mc, atol=tolerance, rtol=tolerance)

True