[Numpy-discussion] Numpy performance vs Matlab.

[Robert Kern]

On Fri, Jan 9, 2009 at 11:32, Nicolas ROUX <nicolas.roux at st.com> wrote:
> Thanks !
>
> -1- The code style is good and the performance vs matlab is good.
> With 400x400:
>  Matlab = 1.56 sec  (with nested "for" loop, so no optimization)
>  Numpy  = 0.99 sec  (with broadcasting)
>
>
> -2- Now with the code below I have strange result.
> With w=h=400:
>   - Using "slice"        =>  0.99 sec
>   - Using "numpy.ogrid"  =>  0.01 sec
>
> With w=400 and h=300:
>   - Using "slice",       => 0.719sec
>   - Using "numpy.ogrid", => broadcast ERROR !
>
> The last broadcast error is:
>  "ValueError: shape mismatch: objects cannot be broadcast to a single shape"
>
> #######################################################
> def test():
>    w = 400
>
>    if 1:  #---Case with different w and h
>        h = 300
>    else:  #---Case with same w and h
>        h = 400
>
>    a = numpy.zeros((h,w))      #Normally loaded with real data
>    b = numpy.zeros((h,w,3))
>
>    if 1:   #---Case with SLICE
>        w0 = slice(0,w-2)
>        w1 = slice(1,w-1)
>        w2 = slice(2,w)
>        h0 = slice(0,h-2)
>        h1 = slice(1,h-1)
>        h2 = slice(2,h)
>    else:   #---Case with OGRID
>        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, p01, p02 = [h0,w0], [h0,w1],[h0,w2]
>    p10, p11, p12 = [h1,w0], [h1,w1],[h1,w2]
>    p20, p21, p22 = [h2,w0], [h2,w1],[h2,w2]
>
>    b[p11+[1]] =  a[p11] - 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]])
> #######################################################
>
> It seems "ogrid" got better performance, but broadcasting is not working any
> more.
> Do you think it's normal that broadcast is not possible using ogrid and
> different w & h ?
> Did I missed any row/colomn missmatch ???

There are several things wrong. Please read this document for
information about how indexing works in numpy.

  http://docs.scipy.org/doc/numpy/user/basics.indexing.html

But basically, you want slices. Using ogrid correctly will be slower.
FWIW, ogrid with only one argument is fairly pointless. ogrid is
intended to be used with multiple dimensions. If you just need one
argument, use arange() or linspace(). With multiple arguments, ogrid
will align the arrays such that they can be broadcasted as you expect.
Lets take a look at some examples:

In [1]: from numpy import *

In [2]: ogrid[0:5]
Out[2]: array([0, 1, 2, 3, 4])

In [3]: ogrid[0:6]
Out[3]: array([0, 1, 2, 3, 4, 5])


Two 1D arrays. Now, if you follow the discussion in the indexing
document, you know that if I were to use these as index arrays, one
for each axis, the indexing mechanism will try to iterate over them in
parallel. Since they have incompatible shapes, this will fail.

Instead, if you put both arguments into ogrid:

In [4]: ogrid[0:5, 0:6]
Out[4]:
[array([[0],
       [1],
       [2],
       [3],
       [4]]),
 array([[0, 1, 2, 3, 4, 5]])]


We get the kind of arrays you need. These shapes are compatible,
through broadcasting, and together form the indices to select out the
part of the matrix you are interested in.

However, just using the slices on the matrix instead of passing the
slices through ogrid is faster.

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth."
  -- Umberto Eco