0

Here is some code I wrote in Python / Numpy that I pretty much directly translated from MATLAB code. When I run the code in MATLAB on my machine, it takes roughly 17 seconds. When I run the code in Python / Numpy on my machine, it takes roughly 233 seconds. Am I not using Numpy effectively? Please look over my Python code to see if I'm using Numpy in a non effective manner. 

import numpy as np
from numpy import * 
import pylab as py
from pylab import *
import math
import time

def heat(D,u0,q,tdim):
    xdim = np.size(u0)
    Z = np.zeros([xdim,tdim])
    Z[:,0]=u0;
    for i in range(1,tdim):
        for j in range (1,xdim-1):
            Z[j,i]=Z[j,i-1]+ D*q*(Z[j-1,i-1]-2*Z[j,i-1]+Z[j+1,i-1])
    return Z

start_time = time.clock()
L = 10
D = 0.5

s = 0.03  # magnitude of noise

Tmax = 0.2
xdim = 25
tdim = 75 

x = np.linspace(0,L,xdim)
t = np.linspace(0,Tmax,tdim)

dt = t[1]-t[0]
dx = x[1]-x[0]

q = dt/(dx**2)
r1 = 0.75*L
r2 = 0.8*L


################################################
## check the stability criterion dt/(dx^2)<.5 ##
################################################

# Define the actual initial temperature distribution
u0 = np.zeros(xdim)
for i in range(0,xdim):
    if(x[i]>=r1 and x[i]<=r2):
        u0[i] = 1
xDat = range(1,xdim-1)
tDat = np.array([tdim])
nxDat = len(xDat)
ntDat = 1
tfinal = tdim-1

# synthesize data
Z = heat(D,u0,q,tdim)
u = Z[xDat,tfinal] # matrix
uDat = u + s*randn(nxDat)

# MATLAB PLOTTING
#figure(1);surf(x,t,Z); hold on;
#if ntDat>1, mesh(x(xDat),t(tDat),uDat);
#else set(plot3(x(xDat),t(tDat)*ones(1,nxDat),uDat,'r-o'),'LineWidth',3);
#end; hold off; drawnow


#MCMC run
N = 10000
m = 100
XD = 1.0
X = np.zeros(N)
X[0] = XD
Z = heat(XD,u0,q,tdim)
u = Z[xDat,tfinal]
oLLkd = sum(sum(-(u-uDat)**2))/(2*s**2)
LL = np.zeros(N)
LL[0] = oLLkd

# random walk step size
w = 0.1
for n in range (1,N):
    XDp = XD+w*(2*rand(1)-1)
    if XDp > 0:
        Z = heat(XDp,u0,q,tdim)
        u = Z[xDat,tfinal]
        nLLkd = sum(sum( -(u-uDat)**2))/(2*s**2)
        alpha = exp((nLLkd-oLLkd))
        if random() < alpha:     
            XD = XDp
            oLLkd = nLLkd
            CZ = Z
    X[n] = XD;
    LL[n] = oLLkd;

print time.clock() - start_time, "seconds"
Zack
  • 713
  • 2
  • 8
  • 21
  • 6
    In general, if you want people to look over your code, post to codereview, not stackoverflow. Second, this isn't a real question-- a question would be, "Why does Numpy/python take so much longer than matlab?" An answer to your original question, "Please look over my code..." would be, quite simply, "Sure" or "No". – kreativitea Oct 18 '12 at 22:37
  • 5
    There's a lot of extraneous stuff in your work that has nothing to do with performance. If you want people to help you, you should at least make the effort to reduce your code the bare minimum necessary to exhibit the problem. You might even solve the problem yourself in the process. – Marcelo Cantos Oct 18 '12 at 22:39
  • 1
    for one, you're running most of your stuff in global scope, which is quite slow, try looking for python optimisation tips. Your code is very far from optimal. – Serdalis Oct 18 '12 at 22:44
  • I'm not entirely sure why MATLAB would be expected to be "slow" :) –  Oct 18 '12 at 22:53
  • 5
    You have a double for loop, that probably takes most time. Matlab has a just in time compiler for such extremely simple for loops. You should vectorize it, so that it is all done in a single array expression then the difference probably vanishes. Or use cython which creates C code with a bit of extra. – seberg Oct 18 '12 at 22:59

1 Answers1

1

Difference in performance between Numpy and Matlab for basic array/matrix operation is most probably due to Numpy being installed against a slower Lapack implementation. For maximum performance you may consider building numpy against LAPACK (instructions here).

In your code the main performance hit is that you are basically performing a convultion in a python for-loop. Thus your are not really making use of Numpy capabilities. You should replace your double for loop with a dedicated convolution function such as scipy.ndimage.convolve.

Nicolas Barbey
  • 6,639
  • 4
  • 28
  • 34