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

[Francesc A. <fa...@op...>]

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 =3D NA_InputArray(a, toenum[a._type], C_ARRAY)
> >     carrb =3D NA_InputArray(b, toenum[b._type], C_ARRAY)
> >     carrc =3D NA_InputArray(c, toenum[c._type], C_ARRAY)
> >     da =3D <double *>carra.data
> >     db =3D <double *>carrb.data
> >     dc =3D <double *>carrc.data
> >     dim1 =3D carra.dimensions[0]
> >     dim2 =3D carra.dimensions[1]
> >     for i from 0<=3D i < dim1:
> >       dc[i] =3D 0.
> >       for j from 0<=3D j < dim2:
> >         dc[i] =3D 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 w=
hy
> > I'm so enthusiastic with it.
>
> That actually looks more like C than Python to me. As soon as I am doin=
g
> pointer arithmetic, I don't feel like I'm writng Python. Would it be al=
l
> 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 co=
uld
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] =3D 0.
      for j in range(a.shape[1]):
        dc[i] =3D 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 li=
ke
types and constructions, it just don't come for free. I just find this wa=
y
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 us=
e 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 y=
our
> > 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 som=
e
> > 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), yo=
u
are accessing directly the array data buffer, with no copies (whenever yo=
ur
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 referrin=
g to
> > cdef'd c-typed arrays in Pyrex?. I think this should be a dangerous t=
hing
> > 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 tha=
t
anymore. That's the price to pay for speed.

About implementing slicing in Pyrex C-pointer arithmetics, well, it can b=
e
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 i=
t
> > cannot be possible with Pyrex?. Mmm, perhaps there is some pre-proces=
sor
> > 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. Als=
o
> 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 hav=
e
> 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, althou=
gh
as I don't have experience on that kind of issues, chances are that my
question might result a complete nonsense :).

Cheers,

--=20
Francesc Alted