CODEWITHSUNDEEP
pythonnumpyoptimization
Rahul Rawat
Rahul Rawat

Reputation: 81

Python code optimization using vectorizing iterative section

I am new to python and numpy there is this Second order PDE code which i want to vectorize to run in lesser time but since it uses function to fill each value in grid i m stuck.

def func(x, y):
    return (-2 * np.pi ** 2) * np.sin(np.pi * x) * np.sin(np.pi * y)


def f1t(x, y):
    return (np.sin(np.pi * x) * np.sin(np.pi * y))


def five_point(grid, i, j, h, grid_x, grid_y):
    return ((grid[i + 1, j] + grid[i - 1, j] + grid[i, j + 1] + grid[i, j - 1]) / 4
            - ((h ** 2) / 4) * func(grid_x[i, 0], grid_y[0, j]))


def five_point_fin_int(X=np.ones(1), Y=np.ones(1), n_x=32, K_max=1000,
                       tol=0.0001, tol_type="grid"):
    import time;
    t0 = time.clock()

    h = 1 / n_x
    X_max = int(X / h)
    Y_max = int(Y / h)

    grid_x, grid_y = np.mgrid[0:X + h:h, 0:Y + h:h]
    grid_true = np.zeros((X_max + 1, Y_max + 1))

    for i in range(1, X_max):
        for j in range(1, Y_max):
            grid_true[i, j] = f1t(grid_x[i, 0], grid_y[0, j])

    grid = np.zeros((X_max + 1, Y_max + 1))
    grid_err = np.zeros((X_max + 1, Y_max + 1))

    count = 0
    tol_max = False
    while ((count < K_max) & (tol_max == False)):
        count += 1

        for i in range(1, X_max):
            for j in range(1, Y_max):
                grid[i, j] = five_point(grid, i, j, h, grid_x, grid_y)
                grid_err[i, j] = (grid[i, j] - grid_true[i, j])

                if (tol_type.lower() == "grid" ):
                    if (np.amax(abs(grid_err)) < tol):
                        tol_max = True
                else:
                    if (abs(np.linalg.norm(grid) - np.linalg.norm(grid_true)) < tol):
                        tol_max = True

    cpu_time = time.clock() - t0

In the end i print compute time since its nested for loops right now the time taken is a lot around 9 sec i would like to improvise on this.

Upvotes: 0

Views: 81

Answers (1)

Tarifazo
Tarifazo

Reputation: 4343

numpy allows you to replace loops by vector calls. You can definitely do the following:

grid_true = np.zeros((X_max + 1, Y_max + 1))
grid_true[1:X_max,1:Y_max]=f1t(*np.meshgrid(grid_x[1:X_max,0], grid_y[0,1:Y_max]))

And you can also try the following:

grid = np.zeros((X_max + 1, Y_max + 1))
grid[1:-1, 1:-1] = five_point(grid, *np.meshgrid(np.arange(1,X_max), np.arange(1, Y_max)),
h, grid_x, grid_y)

However this is not pure "upstream" integration like the one you are doing, since you are essentialy calculating all the grid together in each step (your call!).

Probably a minimization routine could do better. There isn't much difference in performance between numpy and pure python for short loops or small vectors.

Upvotes: 1

Related Questions