def test_eval_splines():
        
    from interpolation.splines.eval_splines import eval_linear, eval_cubic
    from interpolation.multilinear.mlinterp import mlinterp
    import numpy as np


    N = 1000
    K = 100
    d = 2
    grid = (
        (0.0, 1.0, K),
    )*d
    grid_nu = (
        (0.0, 1.0, K),
        np.linspace(0,1,K)
    )
    n_x = 6

    V = np.random.random((K,K))
    VV = np.random.random((K,K,n_x))

    C = np.random.random((K+2,K+2))
    CC = np.random.random((K+2,K+2,n_x))

    points = np.random.random((N,2))

    out = np.random.random((N))
    Cout = np.random.random((N,n_x))

    eval_linear(grid, V, points, out)
    # eval_cubic(grid, C, points[0,:], out[0,0])
    eval_linear(grid, VV, points, Cout)
    eval_linear(grid, VV, points[0,:], Cout[0,:])


    mlinterp(grid, V,points)
    eval_linear(grid, V, points)
    eval_linear(grid, V, points)
    eval_linear(grid, V, points, out)

    print("OK 2")
    eval_linear(grid, V, points[0,:])
    res_0 = eval_linear(grid, VV, points)
    res_1 = eval_linear(grid, VV, points[0,:])


    #nonuniform grid:
    res_0_bis = eval_linear(grid_nu, VV, points)
    res_1_bis = eval_linear(grid_nu, VV, points[0,:])

    assert abs(res_0-res_0_bis).max()<1e-10
    assert abs(res_1-res_1_bis).max()<1e-10


    print("OK 3")
    eval_cubic(grid, C, points, out)
    # eval_cubic(grid, C, points[0,:], out[0,0])
    eval_cubic(grid, CC, points, Cout)
    eval_cubic(grid, CC, points[0,:], Cout[0,:])

    print("OK 4")
    eval_cubic(grid, C, points)
    eval_cubic(grid, C, points[0,:])
    eval_cubic(grid, CC, points)
    eval_cubic(grid, CC, points[0,:])

    print("OK 5")

    ####
    ###
    ###

