CODEWITHSUNDEEP
pythonpython-3.x
Grzegorz
Grzegorz

Reputation: 11

how to do this interpolation formula faster?

I have a function of Everett interpolation and I'd like to make it a little bit faster than it is right now. It works very well b

x and y are the parameters: time and values. xi is the time which I want to have interpolated value.

def everett(x,y,xi):
    '''
    function everett

    INPUT:
    x       list float      
    y       list float      
    xi      float           

    RETURN:
    yi      float


    '''
    n = len(x)       #interpolation degree

    h = x[1]-x[0]       #differences between epochs

    D = np.zeros([n,n+1])
    D[:,0] = x
    D[:,1] = y

    for j in range(1,n):                    #loop to each column
        for i in range(0,n-j):              #loop to cell within a column
            D[i,j+1] = D[i+1,j] - D[i,j]

    #Finding the value of u
    for i in range(0,n):
        u = ( xi - x[i] ) / h

        if u == 0:
            return y[i]
        elif( u > 0 and u < 1.0 ):
            break

    if i == n-1:
        return None

    z = i
    w = 1 - u
    i = 0
    yi = 0
    m1 = u
    m2 = w
    for j in range(1,n+1,2):
        yi += m1 * D[z+1-i,j] + m2 * D[z-i,j]

        i = i + 1
        m1 *= (u - i) * (u + i) / ((j+1)*(j+2))
        m2 *= (w - i) * (w + i) / ((j+1)*(j+2))

        if (z-i)<0 or (z+1-i)>(n-j-1):
            break   #//checks validity of index in the table

    return yi

Thx!

EDIT: some modification using numpy

I change this part of code:

#Finding the value of u
for i in range(0,n):
    u = ( xi - x[i] ) / h

    if u == 0:
        return y[i]
    elif( u > 0 and u < 1.0 ):
        break

if i == n-1:
    return None

by this one:

#Finding the value of u
u = (xi - x) /h
u0 = np.where(u == 0)[0]

if u0.size:
    return y[u0[0]]

i = np.where((u > 0) & (u < 1.0))[0]
if not i.size:
    return None

z = i[0]
u = u[z]

the biggest problem I have right now is how to modify the last loop and the first loop where variable D is filled with values. Any ideas?

Upvotes: 0

Views: 207

Answers (1)

Timtro
Timtro

Reputation: 428

Include numpy, put the data into numpy.array()s and use numpy operations. You'll simplify your code and get, potentially, orders of magnitude better performance. If you're comfortable with Matlab, you'll find numpy easy to learn.

Loops like

for i in range(0,n):
    u = ( xi - x[i] ) / h

become simple one liners:

u = (xi - x) / h

(where x is an array, u will be an array and the - and / will do element-wise arithmetic if xi and h are numbers)

This even works for whole arrays. For example, a forward difference can be expressed in 1D as 

Dx = X[1:] - x[:1]

The X[1:] means the elements of X excluding the first and X[:1] means the elements of X excluding the last.

You can do the same on N-dimensional arrays, eliminating nested loops.

I wrote this article a long time ago, but it's still relevant. You'll see where I use numpy to speed up a finite difference calculation on a mesh (solving the 2D diffusion equation) while also simplifying the code: http://www.timteatro.net/2010/10/29/performance-python-solving-the-2d-diffusion-equation-with-numpy/

If I get a chance, I will come back and help you work on your specific code. But actually, I think this algorithm is a perfect project to introduce yourself to numpy.

And, if you're more interested in the result than you are the method, SciPy (an extension of numpy) has interpolation functions:

http://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html

Upvotes: 1

Related Questions