.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/02_resampling_using_1d_interpolation/02_03_compare_different_splines.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. or to run this example in your browser via JupyterLite or Binder .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_02_resampling_using_1d_interpolation_02_03_compare_different_splines.py: Compare different splines ========================= Obtain a spline through different methods and compare the results. 1. Assume that a user-provided 1D list of samples :math:`f[k]` has been obtained by sampling a spline on a unit grid. 2. From the samples, recover the continuously defined spline :math:`f(x)`. 3. Resample :math:`f(x)` to get :math:`g[k] = f(Tk)`, with :math:`|T| > 1`. 4. Create a new spline :math:`g(x)` from the samples :math:`g[k]`. 5. We define :math:`h(x) = g(x / T)`. 6. Compute the mean squared error (MSE) between :math:`f` and :math:`h`. .. GENERATED FROM PYTHON SOURCE LINES 25-27 Imports ------- .. GENERATED FROM PYTHON SOURCE LINES 27-80 .. code-block:: Python import numpy as np import matplotlib.pyplot as plt from matplotlib.gridspec import GridSpec from splineops.interpolate.tensorspline import TensorSpline plt.rcParams.update({ "font.size": 14, # Base font size "axes.titlesize": 18, # Title font size "axes.labelsize": 16, # Label font size "xtick.labelsize": 14, "ytick.labelsize": 14 }) number_of_samples = 27 f_support = np.arange(number_of_samples) f_support_length = len(f_support) # It's equal to number_of_samples f_samples = np.array([ -0.657391, -0.641319, -0.613081, -0.518523, -0.453829, -0.385138, -0.270688, -0.179849, -0.11805, -0.0243016, 0.0130667, 0.0355389, 0.0901577, 0.219599, 0.374669, 0.384896, 0.301386, 0.128646, -0.00811776, 0.0153119, 0.106126, 0.21688, 0.347629, 0.419532, 0.50695, 0.544767, 0.555373 ]) plot_points_per_unit = 12 # Interpolated signal base = "bspline3" mode = "mirror" f = TensorSpline(data=f_samples, coordinates=f_support, bases=base, modes=mode) f_coords = np.array([q / plot_points_per_unit for q in range(plot_points_per_unit * f_support_length)]) # Syntax hint: pass (plot_coords,) not plot_coords f_data = f(coordinates=(f_coords,), grid=False) val_T = np.pi g_support_length = round(f_support_length // val_T) g_support = np.arange(g_support_length) f_resampled_coords = np.array([q * val_T for q in range(g_support_length)]) g_samples = f(coordinates=(f_resampled_coords,), grid=False) g = TensorSpline(data=g_samples, coordinates=g_support, bases=base, modes=mode) g_coords = np.array([q/plot_points_per_unit for q in range(plot_points_per_unit * len(g_support))]) g_data = g(coordinates=(g_coords,), grid=False) .. GENERATED FROM PYTHON SOURCE LINES 81-94 Expand g to Obtain h -------------------- To compare :math:`g` on the same domain as :math:`f`, we expand :math:`g` by defining a new function :math:`h` as .. math:: h(x) = g\bigl(\tfrac{x}{T}\bigr), where :math:`g` is the continuously defined spline built from the discrete points :math:`g[k]`. Hence, :math:`h` and :math:`f` have the same support and can be directly compared (e.g., by computing an MSE). .. GENERATED FROM PYTHON SOURCE LINES 94-208 .. code-block:: Python fig2 = plt.figure(figsize=(12, 12)) gs2 = GridSpec( nrows=3, ncols=2, width_ratios=[g_support_length, f_support_length - g_support_length], height_ratios=[1, 1, 1] # three equal rows ) # TOP ROW: f + f spline + discrete g[k] ax_top = fig2.add_subplot(gs2[0, :]) # spans both columns ax_top.set_title("Interpolated f spline") # Replot discrete f[k] as stems ax_top.stem(f_support, f_samples, basefmt=" ", label="f[k] samples") # Replot f spline ax_top.plot(f_coords, f_data, color="green", linewidth=2, label="f spline") # Overplot discrete g[k] in red squares at x = k*val_T x_g = np.arange(g_support_length) * val_T ax_top.plot( x_g, g_samples, "rs", mfc='none', markersize=12, markeredgewidth=2, label="g[k] samples" ) # Horizontal line at 0 ax_top.axhline(0, color="black", linewidth=1, zorder=0) ax_top.set_xlim(0, f_support_length - 1) ax_top.set_xticks(np.arange(0, f_support_length, 1)) ax_top.set_xlabel("x") ax_top.set_ylabel("f") ax_top.legend() ax_top.grid(True) # MIDDLE ROW: discrete g + g spline ax_mid_left = fig2.add_subplot(gs2[1, 0]) # left cell ax_mid_right = fig2.add_subplot(gs2[1, 1]) # right cell ax_mid_right.axis("off") # keep it blank ax_mid_left.set_title("Interpolated g spline") # Plot discrete g[k] with red stems, unfilled squares ax_mid_left.vlines( x=g_support, ymin=0, ymax=g_samples, color='red', linestyle='-', linewidth=1 ) ax_mid_left.plot( g_support, g_samples, "rs", mfc='none', markersize=12, markeredgewidth=2, label="g[k] samples" ) # Plot g spline ax_mid_left.plot( g_coords, g_data, color="purple", linewidth=2, label="g spline" ) ax_mid_left.axhline(0, color='black', linewidth=1, zorder=0) ax_mid_left.set_xlim(0, g_support_length - 1) ax_mid_left.set_xticks(np.arange(0, g_support_length, 1)) ax_mid_left.set_xlabel("x") ax_mid_left.set_ylabel("g") ax_mid_left.legend() ax_mid_left.grid(True) # Match y-limits with top row: ax_mid_left.set_ylim(ax_top.get_ylim()) # BOTTOM ROW: expanded h(x) = g(x / λ) ax_bottom = fig2.add_subplot(gs2[2, :]) # spans both columns ax_bottom.set_title("h spline and difference (f - h)") # We'll sample h over 0..(f_support_length-1) h_coords = f_coords # Evaluate h(x) = g(x/val_T) h_data = g(coordinates=(h_coords / val_T,), grid=False) # Also evaluate f at the same coords, so we can show the difference f_data_for_diff = f(coordinates=(h_coords,), grid=False) diff_data = f_data_for_diff - h_data # Plot h in blue ax_bottom.plot(h_coords, h_data, color="blue", linewidth=2, label="h(x)") # Plot difference f - h in red, dashed ax_bottom.plot(h_coords, diff_data, color="red", linestyle="--", linewidth=2, label="f - h") ax_bottom.axhline(0, color='black', linewidth=1, zorder=0) # The domain is the same as f ax_bottom.set_xlim(0, f_support_length - 1) ax_bottom.set_xticks(np.arange(0, f_support_length, 1)) ax_bottom.set_xlabel("x") ax_bottom.set_ylabel("h") ax_bottom.grid(True) ax_bottom.legend() # Match y axis with top row ax_bottom.set_ylim(ax_top.get_ylim()) fig2.tight_layout() plt.show() .. image-sg:: /auto_examples/02_resampling_using_1d_interpolation/images/sphx_glr_02_03_compare_different_splines_001.png :alt: Interpolated f spline, Interpolated g spline, h spline and difference (f - h) :srcset: /auto_examples/02_resampling_using_1d_interpolation/images/sphx_glr_02_03_compare_different_splines_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 209-229 MSE Between f and h ------------------- We compute the MSE between :math:`h(x)` and :math:`f(x)` as .. math:: \text{MSE} = \frac{1}{b - a} \int_{a}^{b} (f(x) - h(x))^2 \, \mathrm{d}x. **Riemann Approximation** To estimate this integral, we discretize the interval :math:`[a,b]` into :math:`K` points. At each point :math:`x_k`, we evaluate :math:`(f(x_k) - h(x_k))^2` and multiply by the width :math:`\Delta x`. Summing across all points produces the approximation .. math:: \Delta x \sum_{k=1}^{K} (f(x_k) - h(x_k))^2 \;\approx\; \int_{a}^{b} (f(x) - h(x))^2 \, \mathrm{d}x. The normalization by :math:`(b-a)` yields the MSE. .. GENERATED FROM PYTHON SOURCE LINES 229-250 .. code-block:: Python # 1) Define a midpoint sampling domain for [a, b] N = 1000 padding_fraction = 0.2 # We avoid artifacts near the edges by excluding part of the domain from each side a = (f_support_length - 1) * padding_fraction b = (f_support_length - 1) * (1 - padding_fraction) dx = (b - a) / N mid_x = np.linspace(a + dx/2, b - dx/2, N) # midpoints # 2) Evaluate f(x) and h(x) at those midpoints f_mid = f(coordinates=(mid_x,), grid=False) h_mid = g(coordinates=(mid_x / val_T,), grid=False) # 3) Compute the midpoint Riemann sum for ∫(f(x)-h(x))^2 dx squared_diff = (f_mid - h_mid) ** 2 integral_value = np.sum(squared_diff) * dx # 4) Divide by (b - a) to get the MSE mse_midpoint = integral_value / (b - a) print(f"MSE between f and h = {mse_midpoint:.6e}") .. rst-class:: sphx-glr-script-out .. code-block:: none MSE between f and h = 1.152514e-03 .. GENERATED FROM PYTHON SOURCE LINES 251-256 Variation with Linear Splines ----------------------------- We repeat exactly everything with linear splines. As the MSE increases, we conclude that splines of degree 3 provide a better representation of the original signal than splines of degree 1. .. GENERATED FROM PYTHON SOURCE LINES 256-392 .. code-block:: Python base = "bspline1" mode = "mirror" # 1) Rebuild linear-spline version of f, then g f_lin = TensorSpline(data=f_samples, coordinates=f_support, bases=base, modes=mode) g_lin_samps = f_lin(coordinates=(f_resampled_coords,), grid=False) g_lin = TensorSpline(data=g_lin_samps, coordinates=g_support, bases=base, modes=mode) # 2) Evaluate them at the same plotting coordinates f_lin_f = f_lin(coordinates=(f_coords,), grid=False) # f_lin over domain 0..(K-1) g_lin_g = g_lin(coordinates=(g_coords,), grid=False) # g_lin over domain 0..(g_support_length-1) h_lin_h = g_lin(coordinates=(f_coords / val_T,), grid=False) # h_lin(x)=g_lin(x/λ) over 0..(K-1) # 3) Create the 3×2 figure layout fig3 = plt.figure(figsize=(12, 12)) gs3 = GridSpec( nrows=3, ncols=2, width_ratios=[g_support_length, f_support_length - g_support_length], height_ratios=[1, 1, 1] ) # TOP ROW: entire row (two columns combined) for f ax_top = fig3.add_subplot(gs3[0, :]) ax_top.set_title("Linear f spline") # Plot f[k] as stems ax_top.stem(f_support, f_samples, basefmt=" ", label="f[k] samples") # Plot f spline ax_top.plot(f_coords, f_lin_f, color="green", linewidth=2, label="f spline") # Overplot discrete g[k] as unfilled red squares at x = k * val_T x_g = np.arange(g_support_length) * val_T ax_top.plot( x_g, g_lin_samps, "rs", # red squares mfc='none', # unfilled markersize=12, markeredgewidth=2, label="g[k] samples" ) ax_top.axhline(0, color='black', linewidth=1, zorder=0) ax_top.set_xlim(0, f_support_length - 1) ax_top.set_xticks(np.arange(0, f_support_length, 1)) ax_top.set_xlabel("x") ax_top.set_ylabel("f") ax_top.grid(True) ax_top.legend() # MIDDLE ROW: left subplot shows g domain, right subplot is blank ax_mid_left = fig3.add_subplot(gs3[1, 0]) ax_mid_right = fig3.add_subplot(gs3[1, 1]) ax_mid_right.axis("off") # keep right side blank ax_mid_left.set_title("Linear g spline") # Discrete g[k] with stems ax_mid_left.vlines( x=g_support, ymin=0, ymax=g_lin_samps, color='red', linewidth=1 ) ax_mid_left.plot( g_support, g_lin_samps, "rs", mfc='none', markersize=8, markeredgewidth=2, label="g[k]" ) # g spline ax_mid_left.plot( g_coords, g_lin_g, color="purple", linewidth=2, label="g" ) ax_mid_left.axhline(0, color='black', linewidth=1, zorder=0) ax_mid_left.set_xlim(0, g_support_length - 1) ax_mid_left.set_xticks(np.arange(0, g_support_length, 1)) ax_mid_left.set_xlabel("x") ax_mid_left.set_ylabel("g") ax_mid_left.grid(True) ax_mid_left.legend() # Match y-range with top ax_mid_left.set_ylim(ax_top.get_ylim()) # BOTTOM ROW: entire row for h ax_bottom = fig3.add_subplot(gs3[2, :]) ax_bottom.set_title("Linear h spline, h(x)=g(x/λ)") ax_bottom.plot(f_coords, h_lin_h, color="blue", linewidth=2, label="h") # Evaluate f_lin at the same coords, then difference f_lin_for_diff = f_lin(coordinates=(f_coords,), grid=False) diff_lin = f_lin_for_diff - h_lin_h # Plot difference in red, dashed ax_bottom.plot(f_coords, diff_lin, color="red", linestyle="--", linewidth=2, label="f - h") ax_bottom.axhline(0, color='black', linewidth=1, zorder=0) ax_bottom.set_xlim(0, f_support_length - 1) ax_bottom.set_xticks(np.arange(0, f_support_length, 1)) ax_bottom.set_xlabel("x") ax_bottom.set_ylabel("h") ax_bottom.grid(True) ax_bottom.legend() ax_bottom.set_ylim(ax_top.get_ylim()) fig3.tight_layout() plt.show() # 4) Recompute MSE with linear splines using midpoint rule N = 1000 padding_fraction = 0.2 # We avoid artifacts near the edges by excluding part of the domain from each side a = (f_support_length - 1) * padding_fraction b = (f_support_length - 1) * (1 - padding_fraction) dx = (b - a) / N # mid_x are the midpoints of each subinterval mid_x = np.linspace(a + dx/2, b - dx/2, N) # Evaluate f_lin and h_lin at midpoints f_lin_mid = f_lin(coordinates=(mid_x,), grid=False) h_lin_mid = g_lin(coordinates=(mid_x / val_T,), grid=False) # Midpoint Riemann sum for ∫(f_lin - h_lin)² squared_diff_lin = (f_lin_mid - h_lin_mid) ** 2 integral_value_lin = np.sum(squared_diff_lin) * dx # Divide by (b - a) to get MSE mse_lin_midpoint = integral_value_lin / (b - a) print(f"MSE with linear splines = {mse_lin_midpoint:.6e}") .. image-sg:: /auto_examples/02_resampling_using_1d_interpolation/images/sphx_glr_02_03_compare_different_splines_002.png :alt: Linear f spline, Linear g spline, Linear h spline, h(x)=g(x/λ) :srcset: /auto_examples/02_resampling_using_1d_interpolation/images/sphx_glr_02_03_compare_different_splines_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none MSE with linear splines = 2.604823e-03 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.874 seconds) .. _sphx_glr_download_auto_examples_02_resampling_using_1d_interpolation_02_03_compare_different_splines.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: binder-badge .. image:: images/binder_badge_logo.svg :target: https://mybinder.org/v2/gh/splineops/splineops.github.io/main?urlpath=lab/tree/notebooks_binder/auto_examples/02_resampling_using_1d_interpolation/02_03_compare_different_splines.ipynb :alt: Launch binder :width: 150 px .. container:: lite-badge .. image:: images/jupyterlite_badge_logo.svg :target: ../../lite/lab/index.html?path=auto_examples/02_resampling_using_1d_interpolation/02_03_compare_different_splines.ipynb :alt: Launch JupyterLite :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 02_03_compare_different_splines.ipynb <02_03_compare_different_splines.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 02_03_compare_different_splines.py <02_03_compare_different_splines.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 02_03_compare_different_splines.zip <02_03_compare_different_splines.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_