On Fri, Jan 9, 2009 at 11:32, Nicolas ROUX <nicolas.r...@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 _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion