[Psyco-devel] RE: [Numpy-discussion] Interesting Psyco/Numeric/Numarray comparison

[Francesc Alted]

A Dimarts 11 Febrer 2003 18:54, Chris Barker va escriure:
> >
> > def multMatVec(object a, object b, object c):
> >     cdef PyArrayObject carra, carrb, carrc
> >     cdef double *da, *db, *dc
> >     cdef int i, j
> >
> >     carra = NA_InputArray(a, toenum[a._type], C_ARRAY)
> >     carrb = NA_InputArray(b, toenum[b._type], C_ARRAY)
> >     carrc = NA_InputArray(c, toenum[c._type], C_ARRAY)
> >     da = <double *>carra.data
> >     db = <double *>carrb.data
> >     dc = <double *>carrc.data
> >     dim1 = carra.dimensions[0]
> >     dim2 = carra.dimensions[1]
> >     for i from 0<= i < dim1:
> >       dc[i] = 0.
> >       for j from 0<= j < dim2:
> >         dc[i] = dc[i] + da[i*dim2+j] * db[j]
> >
> >     return carrc
> >
> >
> > For me Pyrex is like having Python but with the speed of C. This is why
> > I'm so enthusiastic with it.
>
> That actually looks more like C than Python to me. As soon as I am doing
> pointer arithmetic, I don't feel like I'm writng Python. Would it be all
> that much more code in C?

Doing that in C implies writing the "glue" code. In the past example,
multMatVec is a function *directly* accessible in Python, without any
additional declaration.

Moreover, you can do in Pyrex the same things you do in python, so you could
have written the last piece of code as:

def multMatVec(object a, object b, object c):

    for i in range(a.shape[0]):
      c[i] = 0.
      for j in range(a.shape[1]):
        dc[i] = dc[i] + da[i][j] * db[j]

    return c

but, of course, you get only Python speed.

So, the moral is that C speed is only accessible in Pyrex if you use C like
types and constructions, it just don't come for free. I just find this way
to code to be more elegant than using Swig, or other approaches. But I'm
most probably biased because Pyrex is my first (and unique) serious tool for
doing Python extensions.

>
> > speed. So, if you need speed, always use pointers to your data and use a
> > bit of pointer arithmetic to access the element you want (look at the
> > example).
>
> Is there really no way to get this to work?
>
> > Of course, you can also define C arrays if you know the boundaries in
> > compilation time and let the compiler do the computations to access your
> > desired element, but you will need first to copy the data from your
> > buffers to the C-array, and perhaps this is a bit inconvenient in some
> > situations.
>
> Couldn't you access the data array of the NumArray directly? I do this
> all the time with Numeric.

Yeah, you can, and both examples shown here (in Numeric and numarray), you
are accessing directly the array data buffer, with no copies (whenever your
original array is well-.behaved, of course).

>
> > Why you are saying that slicing is not supported?. I've checked them (as
> > python expressions, of course) and work well. May be you are referring to
> > cdef'd c-typed arrays in Pyrex?. I think this should be a dangerous thing
> > that can make the pointer arithmetic to slow down because of the
> > additional required checks in the slice range.
>
> Well, there would need to be two value checks per slice. That would be
> significant for small slices, but not for large ones, I'd love to have
> it. It just doesn't feel like Python without slicing, and it doesn't
> feel like NumPy without multi-dimensional slicing.
>

Again, right now, you can use slicing in Pyrex if you are dealing with
Python objects, but from the moment you access to the lower-level
Numeric/numarray buffer and assign to a Pyrex C-pointer, you can't do that
anymore. That's the price to pay for speed.

About implementing slicing in Pyrex C-pointer arithmetics, well, it can be
worth to ask Greg Ewing, the Pyrex author. I'll send him this particular
question and forward his answer (if any) to the list.

> > There can be drawbacks, like the one stated by Perry related with how to
> > construct general Ufuncs that can handle many different combinations of
> > arrays and types, although I don't understand that very well because
> > Numeric and numarray crews already achieved to do that in C, so why it
> > cannot be possible with Pyrex?. Mmm, perhaps there is some pre-processor
> > involved?.
>
> I was curious about this comment as well. I have only had success with
> writing any of my Numeric based extensions for pre-determined types. If
> I had to support additional types (and/or discontiguous and/or rank-N
> arrays), I ended up with a whole pile of case and/or if statements. Also
> kind of slow and inefficient code.
>
> It seems the only way to do this right is with C++ and templates (eg.
> Blitz++), but there are good reasons not to go that route.
>
> Would it really be any harder to use Pyrex than C for this kind of
> thing? Also, would it be possible to take a Pyrex type approach and have
> it do someting template-like: you wright the generic code in Pyrex, it
> generates all the type-specific C code for you.

Well, this is another good question for Greg. I'll try to ask him, although
as I don't have experience on that kind of issues, chances are that my
question might result a complete nonsense :).

Cheers,

-- 
Francesc Alted