<jupyter_start><jupyter_text>Cubic splines with PythonIn this notebook, we show how to use the `interpolation` library to interpolate a function using cubic splines.<jupyter_code># we need to import the library first
from interpolation import splines<jupyter_output><empty_output><jupyter_text>Interpolating a functionOur goal here, consists in approximate a function $R^d \rightarrow R$ given its values on a regularly spaced rectangular grid.We consider the two dimensional function defined by $f(x,y)=sin(\sqrt{x^2+y^2}/\pi)/(\sqrt{x^2+y^2}/\pi)$ if $(x,y) \neq (0,0)$ and by $f(0,0)=1$<jupyter_code>from numba import vectorize
@vectorize(['double(double, double)'])
def f(x,y):
    if x==0 and y==0:
        out = 1.0
    else:
        t = sqrt( x**2+y**2 )
        out = sin(pi*t)/(pi*t)
    return out
# Let's try it:

# The function f can be evaluated at any point
x0 = 0.5
y0 = 0.5
val = f(x0,y0)
print('Value of f at ({},{}) : {}'.format(x0,x0,val))

# Or at any list of points
x = linspace(-1,1,50)
y = linspace(-1,1,50)
vals = f(x,y)

# Plot
plot(x0, val, 'o')
plot(x, f(x,y))<jupyter_output>Value of f at (0.5,0.5) : 0.358187786013244<jupyter_text>The grid over which we approximate the function is defined by its boundaries $[a_i, b_i]$ and the number of points $n_i$ in each dimension $i \in [1,d]$. In each dimension $i\in[1,d]$ the $n_i$ points $(a_i=s^0_i, s^2_i, ..., s^{n_j-1}_i=b_i)$are regularly spaced, meaning that for any $j$, $s^j_i = a_i + j\frac{b_i-a_i}{n_i-1}$. These points are called *nodes*.For the two dimensional function defined above, we define $[a_1,b_1]=[a_2,b_2]=[-5,5]$ and we use 50 points in each dimension.<jupyter_code># define the approximation space
a = [-5.0, -5.0]  # a_1, a_2
b = [5.0, 5.0]    # b_1, b_2
orders = [15,15]     # n_2, n_2<jupyter_output><empty_output><jupyter_text>The cartesian product of these nodes is used to approximate the compact state $\mathcal{S}=[a_1, b_1], ... [a_d, b_d]$, by a finite set of points $S=(s_{i_1,...i_d})_{i_1\in[1,n_1], ..., i_d\in[1,n_d]}$. By definition, the point indexed by i_1,...i_d is simply the tuple $(s_{i_1},...s_{i_d})$.In general, this multivariate index is sufficient and it is not needed to list these points one after all. By convention, we list the points as a 2d arrays, where each line corresponds to a point, and where *the last index varies faster* meaning that points are listed like in $s_{00}, s_{01}, s_{02}, s_{10}, s_{11}, s_{12}$.This grid can be constructed using the mlinspace function, which extends the linspace function to multidimensional spaces:<jupyter_code>from dolo.numeric.misc import mlinspace

list_of_points = mlinspace(a,b,orders)
print('First points:')
print(list_of_points[:10,:])

# optional plot the grid
plot(list_of_points[:,0], list_of_points[:,1],'o')


# vectorized version
#all_points =  [ [x y] for x=linspace(smin[1],smax[1],orders[1]),  y=linspace(smin[2],smax[2],orders[2])] 
#all_points = vcat(all_points...)
#PyPlot.plot(all_points[:,1], all_points[:,2], "x", color="blue")<jupyter_output>First points:
[[-5.         -5.        ]
 [-5.         -4.28571429]
 [-5.         -3.57142857]
 [-5.         -2.85714286]
 [-5.         -2.14285714]
 [-5.         -1.42857143]
 [-5.         -0.71428571]
 [-5.          0.        ]
 [-5.          0.71428571]
 [-5.          1.42857143]]<jupyter_text>The interpolating data is supposed to fit the exact values $v_{i_1, ..., i_d}=f(s_{i_1, ..., i_d})$ of the function on the grid. We compute them below.<jupyter_code>x_array = list_of_points[:,0]
y_array = list_of_points[:,1]
vals = f(x_array, y_array)<jupyter_output><empty_output><jupyter_text>Interpolation object To intepolate this function, we construct an interpolating object:<jupyter_code>from interpolation.splines import CubicSpline
interp = CubicSpline(a,b,orders,vals)
# we can evaluate this interpolator at a given point:
x0 = array([0.0,0.1])
interp(x0)
# to evaluate it many times, it is more efficient to supply an array of points
# where each line containt the coordinates of one point
vec = linspace(-6.0, 6.0, 100)
points = column_stack([ vec*0, vec ])
interpolated_values = interp(points)
display("Array of points")
display((points[:10,:])) # first 10 points
display("Array of interpolated values")
display(interpolated_values[:10]) # first 10 interpolated values<jupyter_output><empty_output><jupyter_text>Let's plot the result:<jupyter_code>figure()
plot(points[:,1], interpolated_values,"-",label="interpolated values", color='blue')
# plot(points[:,1], f(points[:,0], points[:,1]),":", color="black",label="true values")



nodes = linspace(a[1], b[1], orders[1])
plot(nodes, f(0*nodes,nodes),"o", color="red", label="nodes" )

# compare with linear
from interpolation.splines import LinearSpline
plot(points[:,1], LinearSpline(a,b,orders,vals)(points), color='red',linestyle='--' )
legend()<jupyter_output>Setting
[15 15]
(15, 15)<jupyter_text>Multi splinesSometimes, it is useful, to interpolate several functions at the same points. In that case, the MultiCubicSplinecan be used instead of CubicSpline.Let's define $g(x,y) = exp(-x^2-y^2)$<jupyter_code># We define the new function and evaluate it on the same set of points
g = lambda x,y: exp(-x**2-y**2)  # ne need to vectorize it manually since we use only ufuncs here
vals2 = g(x_array, y_array)<jupyter_output><empty_output><jupyter_text>By convention, the values are given as an array whose first index corresponds to the number of the spline.For this reason, we concatenate the two set of values `vals` and `vals2` by creating a new last dimension.Note that in the resulting array `mvals`, vector values are contiguous in memory (if C-ordered).<jupyter_code>mvals = concatenate([vals[:,None], vals2[:,None]],axis=1)
from interpolation.splines import CubicSplines
interp = CubicSplines(a,b,orders,mvals)<jupyter_output>ORders
[15 15]
(15, 15, 2)<jupyter_text>The interpolation is made exactly as before but the returned objects have one more dimension:<jupyter_code>interp(x0)  # returns a one dimensional array with two elements instead of a float
# to evaluate it many times, it is more efficient to supply an array of points
# where each line containt the coordinates of one point
vec = linspace(-6.0, 6.0, 100)
points = column_stack([ vec*0, vec ])
interpolated_values = interp(points)
display("Array of points")
display((points[:10,:])) # first 10 points
display(points.shape)
display("Array of interpolated values")
display(interpolated_values[:10]) # first 10 interpolated values
display(interpolated_values.shape)<jupyter_output><empty_output><jupyter_text>Low level interface (TODO) TODO: explain how basis splines are defined on *knots* and how to directly deal with the spline coefficients.<jupyter_code># fit coefficients of the spline (note the size of the matrix ( 15+2 x 15+2 )
from interpolation.splines.filter_cubic import filter_coeffs
coeffs = filter_coeffs(a,b,orders,vals)
display(coeffs.shape)
# cartesian product with 50 points along each dimension
fine_grid = mlinspace( [0,0], [3,3], [10,10] )
from interpolation.splines.eval_cubic import vec_eval_cubic_spline, vec_eval_cubic_splines_G

# evaluate the spline
V = vec_eval_cubic_spline( a, b, orders, coeffs, fine_grid );

# evaluate with gradient
V, dV = vec_eval_cubic_splines_G( a, b, orders, coeffs[...,None], fine_grid );
V = V[:,0]
dV = dV[:,:,0]
dV_m_1.shape
xm = fine_grid[:,0].reshape((10,10))
ym = fine_grid[:,1].reshape((10,10))
V_mesh = V.reshape((10,10))
contourf(xm,ym,V_mesh,100)
colorbar()
# add gradient field
# see api of quiver at http://matplotlib.org/api/pyplot_api.html
dV_m_1 = dV[:,0].reshape((10,10))
dV_m_2 = dV[:,1].reshape((10,10))
streamplot( linspace(a[0], a[0], 10), linspace(b[1], b[1], 10), dV_m_1, dV_m_2)
quiver( linspace(0,3, 10), linspace(0,3, 10), dV_m_1, dV_m_2)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14249201277955273, 'fraction_numerical': 0.04460063897763578, 'mean_word_length': 5.206185567010309, 'pattern_counts': {'":': 1, '<': 41, '<?xml version=': 0, '>': 41, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '463822', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 3003, 'n_tokens_neox': 2840, 'n_words': 839, 'path': 'examples/cubic_splines_python.ipynb', 'repo_name': 'gboehl/interpolation.py', 'size': '7825'}