[Numpy-discussion] Numpy performance vs Matlab.

[Nicolas ROUX]

Sorry my previous mail was probalby not clear.
This mail was following the tread we had before, so with some discussion
legacy.

I simplified the code to focus only on "what I" need, rather to bother you
with the full code.

I wrote below a code closer to what I need, where you will agree that
vectorization/broadcasting is needed to avoid nested loops.
As I wrote in the 1st mail (added at the end), what is important is to keep
the code not too ugly due to vectorization syntax.
(As written below I try to demonstrate that vectorization/broadcast code
could be as readable as twice nested loop )
The real code we have is even more complex, with processing the array
element using 5x5 neighbours, instead of 3x3.

######################################################
def test6():
    w = 3096
    h = 2048
    a = numpy.zeros((h,w))      #Normally loaded with real data
    b = numpy.zeros((h,w,3))    

    w0 = numpy.ogrid[0:w-2]
    w1 = numpy.ogrid[1:w-1]
    w2 = numpy.ogrid[2:w]

    h0 = numpy.ogrid[0:h-2]
    h1 = numpy.ogrid[1:h-1]
    h2 = numpy.ogrid[2:h]

    p00, p10, p20 = [h0,w0], [h1,w0],[h2,w0]
    p01, p11, p21 = [h0,w1], [h1,w1],[h2,w1]
    p02, p12, p22 = [h0,w2], [h1,w2],[h2,w2]         

    b[p11,1] =  a[p11] + 1.23*a[p22]  \
                     - numpy.min([a[p11]-a[p00],
                                  a[p11]-a[p01],
                                  a[p11]-a[p02],
                                  a[p11]-a[p10],
                                  a[p11]-a[p12],
                                  a[p11]-a[p20],
                                  a[p11]-a[p21],
                                  a[p11]-a[p22]])  \
                   + 0.123*numpy.max([a[p11]-a[p00],
                                      a[p11]-a[p01],
                                      a[p11]-a[p02],
                                      a[p11]-a[p10],
                                      a[p11]-a[p12],
                                      a[p11]-a[p20],
                                      a[p11]-a[p21],
                                      a[p11]-a[p22]])
######################################################

This code above is the one I wish to write but is not working.
I hope you better understand my issue context ;-)

Did I missed something  ?

Thanks for your help.

Cheers,
Nicolas.


> I understand the weakness of the missing JITcompiler in Python vs Matlab,
> that's why I invistigated numpy vectorization/broadcast.
> (hoping to find a cool way to write our code in fast Numpy)
>
> I used the page http://www.scipy.org/PerformancePython to write my code
> efficiently in Numpy.
> While doing it I found one issue.
>
> To have pretty code, I created p0 and p1 arrays of indexes.

I must admit I don't quite understand what you are trying to do, and what
your problem is.

If you just want to do

Out[:,:] = In[:,:]

there is no need for meshgrids (ogrid), for-loops, or whatever.

It is brain dead to use nested for-loops or ogrid for this purpose in
NumPy. It is equally foolish to use nested for loops or meshgrid for this
purpose in Matlab. If you do, I would seriously question your competence.

By the way, you can index ogrid with more than one dimension:

p = numpy.ogrid[:m, :n]
Out[p] = In[p]

============================================================================
===============
============================================================================
===============
============================================================================
===============

Hi !

Thanks a lot for your fast/detailed reply.
A very good point for Numpy ;-)

I spent all my time trying to prepare my testcase to better share with you,
that's why I didn't reply fast.

I understand the weakness of the missing JITcompiler in Python vs Matlab,
that's why I invistigated numpy vectorization/broadcast.
(hoping to find a cool way to write our code in fast Numpy)

I used the page http://www.scipy.org/PerformancePython to write my code
efficiently in Numpy.
While doing it I found one issue.

To have pretty code, I created p0 and p1 arrays of indexes.
In "test8" I wished to see the commented line working, which is not the
case. 
Having to use "ix_" is not pretty enough, and seems to not work with further
dimensions.
Why the comment line is not working ?

############################################
def test8():
    m = 1024
    n = 512
    
    Out = numpy.zeros((m,n))
    In = numpy.zeros((m,n))
  
    p0 = numpy.ogrid[0:m]
    p1 = numpy.ogrid[0:n]

    Out[0:m,0:n] = In[0:m,0:n]
    #Out[p0,p1] = In[p0,p1]   #This doesn't work
    Out[numpy.ix_(p0,p1)] = In[numpy.ix_(p0,p1)]
############################################

What is maybe not clear in the above code, is that I don't want to predefine
all possible ogrid/vector.
The number of possible ogrid/vector is big if in need to define all.
... And this vector definition become more paintful.

So Numpy vector style is fine if i can write something like:
   Out[p0,p1] = In[p0,p1]   #2 dimensions case
And 
   Out[p0,p1,1] = In[p0,p1,1] #3 dimensions case
But is not fine if i have to add ".ix_()" or to multiply the number of
vector definitions.

Below example with 3 dimensions instead of 2.
############################################
def test9():
    m = 1024
    n = 512
    
    Out = numpy.zeros((m,n,3))
    In = numpy.zeros((m,n,3))
  
    p0 = numpy.ogrid[0:m]
    p1 = numpy.ogrid[0:n]

    Out[0:m,0:n,2] = In[0:m,0:n,2]
    #Out[p0,p1,2] = In[p0,p1,2]
    Out[numpy.ix_(p0,p1,2)] = In[numpy.ix_(p0,p1,2)]
############################################

Tanks again for your support ;-)

Cheers,
Nicolas.


============================================================================
===============
============================================================================
===============
============================================================================
===============
Hi,

I need help ;-)
I have here a testcase which works much faster in Matlab than Numpy.
 
The following code takes less than 0.9sec in Matlab, but 21sec in Python.
Numpy is 24 times slower than Matlab !
The big trouble I have is a large team of people within my company is ready
to replace Matlab by Numpy/Scipy/Matplotlib,
but I have to demonstrate that this kind of Python Code is executed with the
same performance than Matlab, without writing C extension.
This is becoming a critical point for us.

This is a testcase that people would like to see working without any code
restructuring.
The reasons are:
- this way of writing is fairly natural.
- the original code which showed me the matlab/Numpy performance differences
is much more complex,
and can't benefit from broadcasting or other numpy tips (I can later give
this code)

...So I really need to use the code below, without restructuring.

Numpy/Python code:
#####################################################################
import numpy
import time

print "Start test \n" 

dim = 3000

a = numpy.zeros((dim,dim,3))

start = time.clock()

for i in range(dim):
    for j in range(dim):
        a[i,j,0] = a[i,j,1]
        a[i,j,2] = a[i,j,0]
        a[i,j,1] = a[i,j,2]

end = time.clock() - start        

print "Test done,   %f sec" % end
#####################################################################

Matlab code:
#####################################################################
'Start test'
dim = 3000;
tic;
a =zeros(dim,dim,3);
for i = 1:dim
    for j = 1:dim 
        a(i,j,1) = a(i,j,2);
        a(i,j,2) = a(i,j,1);
        a(i,j,3) = a(i,j,3);
    end
end
toc
'Test done'
#####################################################################

Any idea on it ?
Did I missed something ?

Thanks a lot, in advance for your help.


Cheers,
Nicolas.







_______________________________________________
Numpy-discussion mailing list
Numpy-discussion at scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion