<jupyter_start><jupyter_text>Linear Interpolation This notebook shows the interface of **consav.linear_interp**, and ensures that the results are exactly the same as for **scipy**.A number of different possibilities exists:1. **Single evaluation** (no preparation, with preparation)2. **Vectorized evaluation** (no preparation, preparation, preparation + monotonicity, preparation + monotonicity + repitition)  Setup<jupyter_code>%load_ext autoreload
%autoreload 2

import time
import numpy as np
from scipy.interpolate import RegularGridInterpolator

# load the module
from consav import linear_interp<jupyter_output><empty_output><jupyter_text>Example  Setup Function for creating **non-linear grids**:<jupyter_code>def nonlinspace(x_min, x_max, n, phi):
    y = np.empty(n)
    y[0] = x_min
    for i in range(1, n):
        y[i] = y[i-1] + (x_max-y[i-1]) / (n-i)**phi
    return y<jupyter_output><empty_output><jupyter_text>**Function**, **grids** and **known values**:<jupyter_code>from numba import njit

# a. functions
@njit
def f(x1,x2,x3):
    return x1**2*x2+x3*x2**2+x3**2

@njit
def fill_value(grid1,grid2,grid3,value):
    for i in range(grid1.size):
        for j in range(grid2.size):
            for k in range(grid3.size):
                value[i,j,k] = f(grid1[i],grid2[j],grid3[k])    

@njit
def grids_points_and_value(Nx1,Nx2,Nx3):

    grid1 = np.linspace(low_x1,high_x1,Nx1)
    grid2 = np.linspace(low_x2,high_x2,Nx2)
    grid3 = np.linspace(low_x3,high_x3,Nx3)

    value = np.empty((Nx1,Nx2,Nx3))
    fill_value(grid1,grid2,grid3,value)
    
    return grid1,grid2,grid3,value

# b. grid points and value
Nx1,Nx2,Nx3 = 100,200,300
low_x1,low_x2,low_x3 = 1.0,1.0,1.0
high_x1,high_x2,high_x3 = 10.0,15.0,12.0
grid1,grid2,grid3,value = grids_points_and_value(Nx1,Nx2,Nx3)<jupyter_output><empty_output><jupyter_text>Draw **random points** to be evaluated:<jupyter_code>Nyi = 100
xi = np.empty((Nyi,3))
xi1 = np.random.uniform(low=0.9*low_x1,high=1.1*high_x1,size=Nyi) 
xi2 = np.random.uniform(low=0.9*low_x2,high=1.1*high_x2,size=Nyi) 
xi3 = np.random.uniform(low=0.9*low_x3,high=1.1*high_x3,size=Nyi)

xi[:,0] = xi1
xi[:,1] = xi2
xi[:,2] = xi3<jupyter_output><empty_output><jupyter_text>We also consider a set of points where **the two first dimensions are constant** and **the last dimension is monotone**:<jupyter_code>xi_mon = np.empty((Nyi,3))
xi_mon[:,0] = xi1[0]
xi_mon[:,1] = xi2[0]
xi_mon[:,2] = xi3_mon = np.sort(xi3)<jupyter_output><empty_output><jupyter_text>Scipy For comparision we use scipy's **RegularGridInterpolator**.<jupyter_code>yi_scipy_interp = RegularGridInterpolator([grid1,grid2,grid3], value, 
                                          method='linear',bounds_error=False,fill_value=None)
yi_scipy = yi_scipy_interp(xi)
yi_mon_scipy = yi_scipy_interp(xi_mon)<jupyter_output><empty_output><jupyter_text>Single evaluation<jupyter_code>for i in range(Nyi):
    yi = linear_interp.interp_3d(grid1,grid2,grid3,value,xi1[i],xi2[i],xi3[i])
    assert np.allclose(yi_scipy[i],yi)<jupyter_output><empty_output><jupyter_text>Vectorized evaluation<jupyter_code>yi = np.empty(Nyi)
linear_interp.interp_3d_vec(grid1,grid2,grid3,value,xi1,xi2,xi3,yi)
assert np.allclose(yi_scipy,yi)<jupyter_output><empty_output><jupyter_text>Single evaluation with preparation<jupyter_code>prep = linear_interp.interp_3d_prep(grid1,grid2,xi1[0],xi2[0],0)
for i in range(Nyi):
    yi = linear_interp.interp_3d_only_last(prep,grid1,grid2,grid3,value,xi1[0],xi2[0],xi3_mon[i])
    assert np.allclose(yi_mon_scipy[i],yi)<jupyter_output><empty_output><jupyter_text>Vectorized evaluation with preparation<jupyter_code>yi = np.empty(Nyi)    
prep = linear_interp.interp_3d_prep(grid1,grid2,xi1[0],xi2[0],Nyi)
linear_interp.interp_3d_only_last_vec(prep,grid1,grid2,grid3,value,xi1[0],xi2[0],xi3_mon,yi)
assert np.allclose(yi_mon_scipy,yi)<jupyter_output><empty_output><jupyter_text>Vectorized evaluation with preparation and monotonicity<jupyter_code>yi = np.empty(Nyi)
prep = linear_interp.interp_3d_prep(grid1,grid2,xi1[0],xi2[0],Nyi)
linear_interp.interp_3d_only_last_vec_mon(prep,grid1,grid2,grid3,value,xi1[0],xi2[0],xi3_mon,yi)
assert np.allclose(yi_mon_scipy,yi)<jupyter_output><empty_output><jupyter_text>After calling **interp_3d_only_last_vec_mon()** some additional information is saved in **prep**. Interpolating for a second time is therefore faster if using **interp_3d_only_last_vec_mon_rep()**. Note, in particular, that the value array can be changed between calls.<jupyter_code>yi = np.empty(Nyi)    
linear_interp.interp_3d_only_last_vec_mon_rep(prep,grid1,grid2,grid3,value,xi1[0],xi2[0],xi3_mon,yi)
assert np.allclose(yi_mon_scipy,yi)<jupyter_output><empty_output><jupyter_text>Timings  **Note:** To run the code below install [interpolation](https://github.com/EconForge/interpolation.py) from EconForge with `pip install git+https://github.com/econforge/interpolation.py.git/`. Note, that their interpolation algorithm does not allow for extrapolation, so the results differ, but the timings should be illustrative. On the other hand, they have multiple interpolation algorithms and an interface very close to that of scipy.<jupyter_code>do_econforge = True
if do_econforge: 
    from interpolation import interp as interp_econforge
    
for _ in range(5):

    Nyi = 2*10**7
    xi = np.empty((Nyi,3))
    xi1 = np.empty(Nyi)
    xi2 = np.empty(Nyi)
    xi1[:] = np.random.uniform(low=0.9*low_x1,high=1.1*high_x1,size=1)[0]
    xi2[:] = np.random.uniform(low=0.9*low_x2,high=1.1*high_x2,size=1)[0]
    xi3 = np.random.uniform(low=0.9*low_x3,high=1.1*high_x3,size=Nyi)

    xi_mon = np.empty((Nyi,3))
    xi_mon[:,0] = xi1
    xi_mon[:,1] = xi2
    xi_mon[:,2] = xi3_mon = np.sort(xi3)
    
    tic = time.time()
    yi_mon_scipy = yi_scipy_interp(xi_mon)
    toc = time.time()
    print(f'scipy: {toc-tic:.1f} secs')
    
    tic = time.time()
    yi = np.empty(Nyi)
    linear_interp.interp_3d_vec(grid1,grid2,grid3,value,xi1,xi2,xi3_mon,yi)
    toc = time.time()
    print(f'interp_3d_vec: {toc-tic:.1f} secs')
    assert np.allclose(yi_mon_scipy,yi)
    
    tic = time.time()
    yi = np.empty(Nyi)
    prep = linear_interp.interp_3d_prep(grid1,grid2,xi1[0],xi2[0],Nyi)
    linear_interp.interp_3d_only_last_vec(prep,grid1,grid2,grid3,value,xi1[0],xi2[0],xi3_mon,yi)
    toc = time.time()
    print(f'interp_3d_only_last_vec: {toc-tic:.1f} secs')    
    assert np.allclose(yi_mon_scipy,yi)

    tic = time.time()
    yi = np.empty(Nyi)
    prep = linear_interp.interp_3d_prep(grid1,grid2,xi1[0],xi2[0],Nyi)
    linear_interp.interp_3d_only_last_vec_mon(prep,grid1,grid2,grid3,value,xi1[0],xi2[0],xi3_mon,yi)
    toc = time.time()
    print(f'interp_3d_only_last_vec_mon: {toc-tic:.1f} secs') 
    assert np.allclose(yi_mon_scipy,yi)
    
    tic = time.time()
    yi = np.empty(Nyi)
    linear_interp.interp_3d_only_last_vec_mon_rep(prep,grid1,grid2,grid3,value,xi1[0],xi2[0],xi3_mon,yi)
    toc = time.time()
    print(f'interp_3d_only_last_vec_mon_rep: {toc-tic:.1f} secs') 
    assert np.allclose(yi_mon_scipy,yi)
    
    if do_econforge:
        tic = time.time()
        yi_mon_econforge = interp_econforge(grid1,grid2,grid3,value,xi_mon)
        toc = time.time()
        print(f'econforge: {toc-tic:.1f} secs')
        
    print('')<jupyter_output>scipy: 11.1 secs
interp_3d_vec: 1.5 secs
interp_3d_only_last_vec: 0.8 secs
interp_3d_only_last_vec_mon: 0.6 secs
interp_3d_only_last_vec_mon_rep: 0.5 secs
econforge: 3.2 secs

scipy: 10.9 secs
interp_3d_vec: 1.5 secs
interp_3d_only_last_vec: 0.8 secs
interp_3d_only_last_vec_mon: 0.6 secs
interp_3d_only_last_vec_mon_rep: 0.5 secs
econforge: 2.1 secs

scipy: 10.8 secs
interp_3d_vec: 1.4 secs
interp_3d_only_last_vec: 0.8 secs
interp_3d_only_last_vec_mon: 0.6 secs
interp_3d_only_last_vec_mon_rep: 0.5 secs
econforge: 2.1 secs

scipy: 11.2 secs
interp_3d_vec: 1.5 secs
interp_3d_only_last_vec: 0.8 secs
interp_3d_only_last_vec_mon: 0.6 secs
interp_3d_only_last_vec_mon_rep: 0.5 secs
econforge: 2.1 secs

scipy: 11.0 secs
interp_3d_vec: 1.4 secs
interp_3d_only_last_vec: 0.8 secs
interp_3d_only_last_vec_mon: 0.6 secs
interp_3d_only_last_vec_mon_rep: 0.5 secs
econforge: 2.2 secs<jupyter_text>Tests in other dimensions  1D  Setup<jupyter_code># a. functions
@njit
def f(x1):
    return x1**2

@njit
def fill_value(grid1,value):
    for i in range(grid1.size):
        value[i] = f(grid1[i])    

@njit
def grids_points_and_value(Nx1):

    grid1 = np.linspace(low_x1,high_x1,Nx1)

    value = np.empty(Nx1)
    fill_value(grid1,value)
    
    return grid1,value

# b. grid points and value
grid1,value = grids_points_and_value(Nx1)

# c. scipy
yi_scipy_interp = RegularGridInterpolator([grid1], value, 
                                          method='linear',bounds_error=False,fill_value=None)<jupyter_output><empty_output><jupyter_text>Tests<jupyter_code>for _ in range(5):
    
    Nyi = 100
    xi1 = np.sort(np.random.uniform(low=0.9*low_x2,high=1.1*high_x2,size=Nyi))
    
    tic = time.time()
    yi_mon_scipy = yi_scipy_interp(xi1)
    toc = time.time()
    
    tic = time.time()
    yi = np.empty(Nyi)
    linear_interp.interp_1d_vec(grid1,value,xi1,yi)
    toc = time.time()
    assert np.allclose(yi_mon_scipy,yi)
        
    tic = time.time()
    yi = np.empty(Nyi)
    prep = linear_interp.interp_1d_prep(Nyi)
    linear_interp.interp_1d_vec_mon(prep,grid1,value,xi1,yi)
    toc = time.time()
    assert np.allclose(yi_mon_scipy,yi)
    
    tic = time.time()
    yi = np.empty(Nyi)
    linear_interp.interp_1d_vec_mon_rep(prep,grid1,value,xi1,yi)
    toc = time.time()
    assert np.allclose(yi_mon_scipy,yi)
    
    tic = time.time()
    yi = np.empty(Nyi)
    linear_interp.interp_1d_vec_mon_noprep(grid1,value,xi1,yi)
    toc = time.time()
    assert np.allclose(yi_mon_scipy,yi)

print('all is good')<jupyter_output>all is good<jupyter_text>2D  Setup<jupyter_code># a. functions
@njit
def f(x1,x2):
    return x1**2*x2+x2**2

@njit
def fill_value(grid1,grid2,value):
    for i in range(grid1.size):
        for j in range(grid2.size):
                value[i,j] = f(grid1[i],grid2[j])    

@njit
def grids_points_and_value(Nx1,Nx2):

    grid1 = np.linspace(low_x1,high_x1,Nx1)
    grid2 = np.linspace(low_x2,high_x2,Nx2)

    value = np.empty((Nx1,Nx2))
    fill_value(grid1,grid2,value)
    
    return grid1,grid2,value

# b. grid points and value
grid1,grid2,value = grids_points_and_value(Nx1,Nx2)

# c. scipy
yi_scipy_interp = RegularGridInterpolator([grid1,grid2], value, 
                                          method='linear',bounds_error=False,fill_value=None)<jupyter_output><empty_output><jupyter_text>Tests<jupyter_code>for _ in range(5):
    
    Nyi = 10000
    xi = np.empty((Nyi,2))
    xi1 = np.empty(Nyi)
    xi1[:] = np.random.uniform(low=0.9*low_x1,high=1.1*high_x1,size=1)[0]
    xi2 = np.random.uniform(low=0.9*low_x2,high=1.1*high_x2,size=Nyi)

    xi_mon = np.empty((Nyi,2))
    xi_mon[:,0] = xi1
    xi_mon[:,1] = xi2_mon = np.sort(xi2)
    
    tic = time.time()
    yi_mon_scipy = yi_scipy_interp(xi_mon)
    toc = time.time()
    
    tic = time.time()
    yi = np.empty(Nyi)
    linear_interp.interp_2d_vec(grid1,grid2,value,xi1,xi2_mon,yi)
    toc = time.time()
    assert np.allclose(yi_mon_scipy,yi)
    
    tic = time.time()
    yi = np.empty(Nyi)
    prep = linear_interp.interp_2d_prep(grid1,xi1[0],Nyi)
    linear_interp.interp_2d_only_last_vec(prep,grid1,grid2,value,xi1[0],xi2_mon,yi)
    toc = time.time()
    assert np.allclose(yi_mon_scipy,yi)

    tic = time.time()
    yi = np.empty(Nyi)
    prep = linear_interp.interp_2d_prep(grid1,xi1[0],Nyi)
    linear_interp.interp_2d_only_last_vec_mon(prep,grid1,grid2,value,xi1[0],xi2_mon,yi)
    toc = time.time()
    assert np.allclose(yi_mon_scipy,yi)
    
    tic = time.time()
    yi = np.empty(Nyi)
    linear_interp.interp_2d_only_last_vec_mon_rep(prep,grid1,grid2,value,xi1[0],xi2_mon,yi)
    toc = time.time()
    assert np.allclose(yi_mon_scipy,yi)
    
print('all is good')<jupyter_output>all is good
{'fraction_non_alphanumeric': 0.13336106489184693, 'fraction_numerical': 0.044509151414309486, 'mean_word_length': 4.765467625899281, 'pattern_counts': {'":': 0, '<': 66, '<?xml version=': 0, '>': 66, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '332969', 'license': 'MIT', 'n_tokens_mistral': 5576, 'n_tokens_neox': 5248, 'n_words': 751, 'path': '04. Tools/01. Linear interpolation.ipynb', 'repo_name': 'alanlujan91/ConsumptionSavingNotebooks', 'size': '12020'}