[Numpy-discussion] Blitz++, NumArray...

[Fernando Perez]

Hi all,

as some of you know, I'm a big fan of weave.inline() for the ease of 
development it offers.  Essentially C++ becomes a scripting part of python, 
with only the occasional recompile overhead (but this is transparent to the 
users).

inline() makes heavy use of the blitz++ library, to represent numpy arrays as 
blitz arrays and therefore index them as Arr(i,j,k...).  This seems to me like 
a very sound decision, as the blitz library was specifically written to 
provide simultaneously fortran-like performance with high-level object 
orientation.

With this in mind, I have a few questions:

- Is there anything in Numarray's design which prevents this same trick from 
being used?  The core of getting a Numpy array into a blitz one is simply a 
call to:

// Convert a Numpy array to a blitz one, using the original's data (no copy)
template<class T, int N>
static blitz::Array<T,N> py_to_blitz(PyArrayObject *arr_obj)
{
     blitz::TinyVector<int,N> shape(0);
     blitz::TinyVector<int,N> strides(0);

     for (int i=0;i<N;i++) {
         shape[i] = arr_obj->dimensions[i];
         strides[i] = arr_obj->strides[i]/sizeof(T);
     }
     return blitz::Array<T,N>((T*) arr_obj->data,shape,strides,
                              blitz::neverDeleteData);
}


I would not be surprised if something much like this would work for Numarrays, 
but I'd like to know for sure.

- If the above is indeed possible for Numarray, do you guys feel that blitz 
should be considered an integral part of python's low-level arsenal, or not? 
I view its advantages especially in the field of an extremely clean, simple 
syntax (compared to the low-level Numpy API, I'm not familiar with Numarray 
yet).  I've seen examples where a numpy algorithm can be recoded using blitz 
in a _far_ smaller amount of code, and with performance gains to boot (in some 
cases, not all).

Please note that I am NOT second-guessing any of the design decisions made by 
the Numarray crowd. I just want to clarify if it's possible to have a C++ API 
which gives access to numarray with the simplest, highest-level syntax and 
behavior which is reasonable.  This is important for teams where not everyone 
is an expert on C pointer manipulation details, but they may be perfectly ok 
with writing loops over things which can be indexed as A(i,j,k).

- If any of this seems worthy of more detailed discussion, I told Eric Jones 
that I'd be happy to coordinate something along these lines at Scipy'03.  If 
there's time in the schedule I'd be glad to do it officially, otherwise we can 
just try to keep it in mind for a lunch/dinner meeting.

Documenting a 'Best Practices for low-level numerical work in python' would 
probably be a beneficial thing to do for the community.  Right now there are 
so many options that I'm sure I'm not the only one to feel a bit confused.

- If Blitz++ is indeed deemed a worthy piece of this puzzle, I wonder if we 
could contact some of its developers to bring them into this discussion.  I 
worry a bit that activity on its mailing lists seems rather slow, but I 
noticed that an apparently central guy, Julian C. Cummings, is at 
CACR/Caltech.  Perhaps he could stop by and let us know what the future for 
blitz looks like.  I can fire off an email to him if there is interest in the 
blitz thing.

Anyway, I'm just trying to 'pre-soak' a few ideas for discussion at scipy'03 
on these topics.

Best regards,

Fernando.