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

[Chris Barker]

Francesc Alted wrote:

> First, define some enum types and headers:

Could all this be put into Pyrex? (when NumArray becomes more stable
anyway) It's well beyond me to understand it.

> I will show here a function to multiply a matrix by a vector (C double
> precision):
> 
> 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?

> 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.

> 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.

> 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.

-Chris

-- 
Christopher Barker, Ph.D.
Oceanographer
                                    		
NOAA/OR&R/HAZMAT         (206) 526-6959   voice
7600 Sand Point Way NE   (206) 526-6329   fax
Seattle, WA  98115       (206) 526-6317   main reception

Chris.Barker at noaa.gov