I'm using the numpy C API (PyArray_SimpleNewFromData) to perform the conversion but my code is written by hands. I would like to simplify it using SWIG but I also would like to see a good typemap valarray <=> numpy.array :) Joris : Historical ones? Maybe also the fact that distutils has some small pb with C++ module. To sum up : You have to compile you module with the same compiler options you used to compile Python. Python is coded in C so some options does not make sense in C++. As a result, you get warnings at compile time (see pylab compiled with gcc for instance). Xavier

A related question, just out of curiosity: is there a technical reason why Numpy has been coded in C rather than C++?

Joris

On 05 Sep 2007, at 02:24, David Goldsmith wrote:

...

Anyone have a well-tested SWIG-based C++ STL valarray <=> numpy.array typemap to share? Thanks!

DG -- ERD/ORR/NOS/NOAA <http://response.restoration.noaa.gov/ emergencyresponse/> _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion

Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm

_______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion

-- ############################################ Xavier Gnata CRAL - Observatoire de Lyon 9, avenue Charles André 69561 Saint Genis Laval cedex Phone: +33 4 78 86 85 28 Fax: +33 4 78 86 83 86 E-mail: gnata@obs.univ-lyon1.fr ############################################

No I hadn't - thanks! (Probably should have Google-d first, huh. :-[ ) DG Travis Vaught wrote:

Have you seen this?

http://www.scipy.org/Cookbook/SWIG_and_NumPy

Also, the numpy/doc/swig directory has the simple typemaps.

Travis

On Sep 5, 2007, at 7:08 AM, Xavier Gnata wrote:

...

I'm using the numpy C API (PyArray_SimpleNewFromData) to perform the conversion but my code is written by hands. I would like to simplify it using SWIG but I also would like to see a good typemap valarray <=> numpy.array :)

Joris : Historical ones? Maybe also the fact that distutils has some small pb with C++ module. To sum up : You have to compile you module with the same compiler options you used to compile Python. Python is coded in C so some options does not make sense in C++. As a result, you get warnings at compile time (see pylab compiled with gcc for instance).

Xavier

...

A related question, just out of curiosity: is there a technical reason why Numpy has been coded in C rather than C++?

Joris

On 05 Sep 2007, at 02:24, David Goldsmith wrote:

...

Anyone have a well-tested SWIG-based C++ STL valarray <=> numpy.array typemap to share? Thanks!

DG -- ERD/ORR/NOS/NOAA <http://response.restoration.noaa.gov/ emergencyresponse/> _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion

Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm

_______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion

-- ############################################ Xavier Gnata CRAL - Observatoire de Lyon 9, avenue Charles André 69561 Saint Genis Laval cedex Phone: +33 4 78 86 85 28 Fax: +33 4 78 86 83 86 E-mail: gnata@obs.univ-lyon1.fr ############################################

_______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion

_______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion

-- ERD/ORR/NOS/NOAA <http://response.restoration.noaa.gov/emergencyresponse/>

The setup.py script in numpy/doc/swig is for compiling test code for numpy.i. It is properly invoked by the Makefile, which will first run swig to generate the wrapper code for the test classes. All a developer, who is using swig to interface some code with numpy in python, needs is numpy.i. The setup.py script could possibly work as a template for whatever they are wrapping, I guess. On Sep 5, 2007, at 3:57 PM, David Goldsmith wrote:

Travis Vaught wrote:

...

Have you seen this?

http://www.scipy.org/Cookbook/SWIG_and_NumPy

Unclear (to me): precisely what does one get from running python numpy/docs/swig/setup.py install, and is the product necessary, and if so, which other components rely on the product? I ask 'cause I'm getting the following error trying to do so:

[Parallels emulating] Microsoft Windows XP [Version 5.1.2600] (C) Copyright 1985-2001 Microsoft Corp.

C:\Python25\Lib\site-packages\numpy\doc\swig>python setup.py install running install running build running build_py file Vector.py (for module Vector) not found file Matrix.py (for module Matrix) not found file Tensor.py (for module Tensor) not found file Vector.py (for module Vector) not found file Matrix.py (for module Matrix) not found file Tensor.py (for module Tensor) not found running build_ext building '_Vector' extension creating build creating build\temp.win32-2.5 creating build\temp.win32-2.5\Release C:\Program Files\Microsoft Visual Studio .NET 2003\Vc7\bin\cl.exe /c /nologo /Ox /MD /W3 /GX /DNDEBUG -IC:\Python25\lib\site-packages\numpy\core \include -IC:\Python25\include -IC:\Python25\PC /TpVector_wrap.cxx /Fobuild\temp.win32-2.5\Release\Vector_wrap.obj cl : Command line warning D4029 : optimization is not available in the standard edition compiler Vector_wrap.cxx c1xx : fatal error C1083: Cannot open source file: 'Vector_wrap.cxx': No such file or directory error: command '"C:\Program Files\Microsoft Visual Studio .NET 2003\Vc7\bin\cl.exe"' failed with exit status 2

I have Python 2.5.1 installed:

C:\Python25\Lib\site-packages\numpy\doc\swig>python Python 2.5.1 (r251:54863, Apr 18 2007, 08:51:08) [MSC v.1310 32 bit (Intel)] on win32

and numpy 1.0.3:

...

...

...

import numpy numpy.__version__ '1.0.3'

I have the swig-1.3.31 exe:

Directory of C:\SWIG\swigwin-1.3.31

11/21/2006 12:07 AM 1,190,652 swig.exe

and it runs (take my word for it)

and the VC++ compiler via Visual Studio .NET 2003 (this I know 'cause I use it frequently).

So, if I don't need the product of python numpy/doc/swig/setup.py install, please explain why I don't, but if I do need it, please help me figure out why I can't build it. Thanks!

DG

...

Also, the numpy/doc/swig directory has the simple typemaps.

Travis

On Sep 5, 2007, at 7:08 AM, Xavier Gnata wrote:

...

I'm using the numpy C API (PyArray_SimpleNewFromData) to perform the conversion but my code is written by hands. I would like to simplify it using SWIG but I also would like to see a good typemap valarray <=> numpy.array :)

Joris : Historical ones? Maybe also the fact that distutils has some small pb with C++ module. To sum up : You have to compile you module with the same compiler options you used to compile Python. Python is coded in C so some options does not make sense in C++. As a result, you get warnings at compile time (see pylab compiled with gcc for instance).

Xavier

...

A related question, just out of curiosity: is there a technical reason why Numpy has been coded in C rather than C++?

Joris

On 05 Sep 2007, at 02:24, David Goldsmith wrote:

...

Anyone have a well-tested SWIG-based C++ STL valarray <=> numpy.array typemap to share? Thanks!

DG -- ERD/ORR/NOS/NOAA <http://response.restoration.noaa.gov/ emergencyresponse/> _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion

Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm

_______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion

-- ############################################ Xavier Gnata CRAL - Observatoire de Lyon 9, avenue Charles André 69561 Saint Genis Laval cedex Phone: +33 4 78 86 85 28 Fax: +33 4 78 86 83 86 E-mail: gnata@obs.univ-lyon1.fr ############################################

_______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion

_______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion

_______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion

** Bill Spotz ** ** Sandia National Laboratories Voice: (505)845-0170 ** ** P.O. Box 5800 Fax: (505)284-5451 ** ** Albuquerque, NM 87185-0370 Email: wfspotz@sandia.gov **

Christopher Barker wrote:

Xavier Gnata wrote:

...

I'm using the numpy C API (PyArray_SimpleNewFromData) to perform the conversion but my code is written by hands.

I'd like to see that. How are you getting the pointer to pass in to PyArray_SimpleNewFromData? It looks like you can do something like:

(VA is a valarray<double>) npy_intp *dims dims[0] = VA.size()

NPA = PyArray_SimpleNewFromData(1, dims, typenum, &VA[0]);

Is that what you're doing? Is there any guarantee that &VA[0] won't change? In any case, I assume that you have to make sure that VA doesn't get deleted while the array is still around.

...

I would like to simplify it using SWIG but I also would like to see a good typemap valarray <=> numpy.array :)

In principle, if you know how to write the code by hand, you know how to write the typemap.

Here it is :) OK it is only a smal code but I often use something like that to debug my C++ code (the goal is just to have a quick look to larrge array to unnderstand what is going wrong). #include <python2.4/Python.h> #include <numpy/arrayobject.h> #include <valarray> using namespace std; int main (int argc, char *argv[]) { PyObject *pylab_module, *pylab_module_dict, *func_imshow, *func_show, *args, *result_imshow, *result_show, *array; int NbDims = 2; int *Dims; long NbData; Dims = new int[NbDims]; Dims[0] = 10; Dims[1] = 10; NbData = Dims[0] * Dims[1]; valarray < double >Data (NbData); for (long i = 0; i < NbData; i++) { Data[i] = (double) i / (NbData - 1); } Py_Initialize (); // Needed before any call to PyArray_foo function. import_array1 (-1); // New reference to a numpy array array = PyArray_SimpleNewFromData (NbDims, Dims, PyArray_DOUBLE, &Data[0]); pylab_module = PyImport_Import (PyString_FromString ("pylab")); if (pylab_module) { pylab_module_dict = PyModule_GetDict (pylab_module); if (pylab_module_dict) { func_imshow = PyDict_GetItemString (pylab_module_dict, "imshow"); if (func_imshow) { func_show = PyDict_GetItemString (pylab_module_dict, "show"); if (func_show) { args = PyTuple_New (1); PyTuple_SetItem (args, 0, array); result_imshow = PyObject_CallObject (func_imshow, args); Py_XDECREF (result_imshow); // We dont use the result... Py_XDECREF (args); result_show = PyObject_CallObject (func_show, NULL); Py_XDECREF (result_show); // We dont use the result... Py_XDECREF (array); } } } Py_XDECREF (pylab_module); } Py_Finalize (); return 0; } "In principle, if you know how to write the code by hand, you know how to write the typemap." Yes but I had to spend a bit more time on that. Hand written code I have posted fit my (debug) needs so I decided not to use SWIG. "Is that what you're doing? Is there any guarantee that &VA[0] won't change? In any case, I assume that you have to make sure that VA doesn't get deleted while the array is still around. " Yes I have but it is not a problem in my simple use cases. Xavier -- ############################################ Xavier Gnata CRAL - Observatoire de Lyon 9, avenue Charles André 69561 Saint Genis Laval cedex Phone: +33 4 78 86 85 28 Fax: +33 4 78 86 83 86 E-mail: gnata@obs.univ-lyon1.fr ############################################

OK, so you've now got a view of the data from the valarray. Nice to know this works, but, of course, fragile if the valarray is re-sized or anything, so it probably won't work for us.

Unless you use a special allocator/desallocator (I don't know if the latter is possible), I don't know how you could first correctly share the pointer (as you said, if the valarray is resized or deleted, you lose your memory). If you want to have a correct implementation (robust to array deletion, resizing in Python and C++), you should use shared pointers, and I don't know if you can build a vector with a shared pointer (perhaps with the special allocator/desallocator, but then it does not work with std::set, std::list or std::map). Matthieu

Joris De Ridder wrote:

A related question, just out of curiosity: is there a technical reason why Numpy has been coded in C rather than C++?

There was a fair bit of discussion about this back when the numarray project started, which was a re-implementation of the original Numeric. IIRC, one of the drivers was that C++ support was still pretty inconsistent across compilers and OSs, particularly if you wanted to really get the advantages of C++, by using templates and the like. It was considered very important that the numpy code base be very portable. C++ compilers have gotten better an more standards compliant, but as a recent thread shows, folks still want to build numpy with older compilers and libs, so the reasoning still holds. Too bad, in a way, I suspect a well-designed C++ numpy could make it easier to write compiled extensions, which would be pretty nice. Of course, it should be possible to write C++ wrappers around the core ND-array object, if anyone wants to take that on! -Chris -- Christopher Barker, Ph.D. Oceanographer Emergency Response Division NOAA/NOS/OR&R (206) 526-6959 voice 7600 Sand Point Way NE (206) 526-6329 fax Seattle, WA 98115 (206) 526-6317 main reception Chris.Barker@noaa.gov

"Bill Spotz" <wfspotz@sandia.gov> writes:

On Sep 5, 2007, at 11:38 AM, Christopher Barker wrote:

...

Of course, it should be possible to write C++ wrappers around the core ND-array object, if anyone wants to take that on!

boost::python has done this for Numeric, but last I checked, they have not upgraded to numpy.

Even then, their wrappers went through the Python interface, not the C API. So, it's no faster than using Python straight. -- |>|\/|< /--------------------------------------------------------------------------\ |David M. Cooke http://arbutus.physics.mcmaster.ca/dmc/ |cookedm@physics.mcmaster.ca

Bill Spotz wrote:

I have been considering adding some C++ STL support to numpy/doc/swig/ numpy.i. Probably std::vector<TYPE> <=> PyArrayObject (and some std::complex<TYPE> support as well). Is this what you had in mind?

well, std::valarray is not the same as std::vector, though there are similarities, so the example would be helpful. Of greatest concern to me is the memory management issue -- it would be nice to be able to share the data black between the valarray and the numpy array (rather than copying back and forth), but as the both valarrays and vectors are re-sizeable, that might get tricky. I'm assuming that you can get the pointer to the data block from both, but it might change on you. If you solve that for std::vector, the solution is likely to be similar for std:valarray. (I hope). David Goldsmith wrote:

Point of clarification: below "well-tested" = "well-use-tested," not (necessarily) "well-unit-tested".

Of course, the better tested the better, but anything is probably better than starting from scratch! Travis Vaught wrote:

Have you seen this?

http://www.scipy.org/Cookbook/SWIG_and_NumPy

Also, the numpy/doc/swig directory has the simple typemaps.

Looked at both, and they are a great starting point, but only deal with plain old C arrays, as far as I've seen. Unless someone speaks up, it sounds like it's not been done, but there are at least a few of us that are interested, so maybe we can start a collaboration. -Chris By the way, David G. and I are working on the same project, so we're kind of like one person.... -- Christopher Barker, Ph.D. Oceanographer Emergency Response Division NOAA/NOS/OR&R (206) 526-6959 voice 7600 Sand Point Way NE (206) 526-6329 fax Seattle, WA 98115 (206) 526-6317 main reception Chris.Barker@noaa.gov

Bill Spotz wrote:

Yes, this resizing memory management issue is the main reason I haven't tried to implement it in numpy.i yet.

<thinking out loud> A possibly better solution would be to develop a class that inherits from std::valarray<TYPE> but also implements the array interface attributes (these attributes would have to be dynamic, in that they check the std::valarray attributes when accessed rather than storing copies).

I like that, though it's way over my head to implement. However, I'm beginning to have my doubts about valarrays. I'm reading: Josuttis, Nicolai M. 1999. "The C+= Standard Library: A Tutorial and Reference" It's 8 years old now, but he writes: "The valarray classes were not designed very well. In fact, nobody tried to determine if the specification worked" He goes on to suggest that Blitz++ might have more of a future. (though it looks like he's involved with the Boost project now) He also points out some major flaws in the text. In reading, I also see that while valarrays can be used for multidimensional arrays, the semantics are pretty ugly. However, he also says: "In principle...you can change their size. However, changing the size of a valarray is provided only to make a two-step initialization (creating and setting the size)" So maybe the memory re-allocation is such an issue. Does anyone know the status of support for valarrays now? Is there another alternative? At the moment, all we need is a one-d fixed size array. There is Boost::array (and MultiArray), but Boost has always seemed like a big, ugly, hard to build dependency. Can you just grab the code for some of these small pieces by themselves? -Chris -- Christopher Barker, Ph.D. Oceanographer Emergency Response Division NOAA/NOS/OR&R (206) 526-6959 voice 7600 Sand Point Way NE (206) 526-6329 fax Seattle, WA 98115 (206) 526-6317 main reception Chris.Barker@noaa.gov

Matthieu Brucher wrote:

Blitz++ is more or less avandoned. It uses indexes than can be not-portable between 32bits platforms and 64bits ones.

Oh well -- that seems remarkably short sited, but would I have done better?

The Boost.Array is a fixed-size array, determined at compile-time,

Ah, I had gotten the wrong impression -- I thought it was fixed at construction time, not compile time.

not interesting there, I suppose.

I agree, I kind of wonder what the point is.

Multiarrays are what you're looking for.

Even if I just want 1-d? though I guess a 1-d multiarray is pretty simple.

Besides, it is not needed to build Boost to use them

I've seen that -- it does look like all we'd need is the header. So, can one: - create a Multiarray from an existing data pointer? - get the data pointer for an existing Multiarray? I think that's what I'd need to make the numpy array <-> Multiarray transition without any copying. ( I know those are really questions that are best asked on the boost list, and should be in the docs, but you folks are so helpful...) Maybe this is the class to wrap with the array interface, though maybe that's exactly what Boost::python::array does (though, AFAICT, still not for numpy). -CHB -- Christopher Barker, Ph.D. Oceanographer Emergency Response Division NOAA/NOS/OR&R (206) 526-6959 voice 7600 Sand Point Way NE (206) 526-6329 fax Seattle, WA 98115 (206) 526-6317 main reception Chris.Barker@noaa.gov

Philip Austin wrote:

Albert Strasheim has done some work on this: http://thread.gmane.org/gmane.comp.python.c++/11559/focus=11560

Thanks for the pointer. Not a lot of docs, and it looks like he's using boost::python, and I want to use SWIG, but I'm sure there's something useful in there. -Chris -- Christopher Barker, Ph.D. Oceanographer Emergency Response Division NOAA/NOS/OR&R (206) 526-6959 voice 7600 Sand Point Way NE (206) 526-6329 fax Seattle, WA 98115 (206) 526-6317 main reception Chris.Barker@noaa.gov

On 9/6/07, Christopher Barker <Chris.Barker@noaa.gov> wrote:

Bill Spotz wrote: However, I'm beginning to have my doubts about valarrays. I'm reading:

Josuttis, Nicolai M. 1999. "The C+= Standard Library: A Tutorial and Reference"

It's 8 years old now, but he writes:

"The valarray classes were not designed very well. In fact, nobody tried to determine if the specification worked"

I've never read that particular book, but I've also read somewhere that valarray is pretty useless for serious work. The timeframe is probably about the same, though -- I think the last time I used valarray was about 7-8 years ago. --bb

Bryan Van de Ven wrote:

Christopher Barker wrote:

...

Does anyone know the status of support for valarrays now?

I used std::valarray to implement a variant of the example Matrix class in Stroustrup's book (2D only) about two years ago. I was aware that is in disuse, by and large, but it worked well enough for my purposes and I was happy with it. I'm sure it could have been done differently/better.

Bryan

Hi, std:valarray are quite strange containers because they are not well integrated in the STL. For instance, you cannot play with maps of valarrays : http://gcc.gnu.org/ml/gcc-bugs/2006-05/msg02397.html Note that it is *not* a gcc bug but a design choice. As a result and beacuse I'm an STL happy user, I always use vector when I have to deal with arrays (of course matrix lib are great as log as you have to do some more complicated maths or as long as you have matrix and not only arrays). Hum...we are on a python list ;) so I would be very happy to see some support of std:vector <--> numpy. Xavier ps : There are no real performance differences between vector and valarray (in my use cases...) -- ############################################ Xavier Gnata CRAL - Observatoire de Lyon 9, avenue Charles André 69561 Saint Genis Laval cedex Phone: +33 4 78 86 85 28 Fax: +33 4 78 86 83 86 E-mail: gnata@obs.univ-lyon1.fr ############################################

On 9/10/07, Christopher Barker <Chris.Barker@noaa.gov> wrote:

Thanks for you input Xavier.

Xavier Gnata wrote:

...

std:valarray are quite strange containers because they are not well integrated in the STL.

...

I always use vector when I have to deal with arrays.

...

ps : There are no real performance differences between vector and valarray (in my use cases...)

Probably not for our use cases either. However, we're not doing a lot of STL either, so I'm not sure there is any downside to valarray. It looks like neither one supports any kind of "view" semantics, so for the purposes of numpy array wrapping, they really aren't any different.

I think that the originator of valarray saying it was misguided might be considered a downside. Now on to your other issue: HOPKINS:

As I see it, *valarray* was designed by BS [Bjarne Stroustrup]

BUDGE: No, Bjarne should not be blamed for *valarray*. If any one person should be blamed for *valarray*, it's me. That's why I tend to sign my postings to this reflector as Kent "Frankenstein" Budge. -- Oops, I guess I've tipped my hand already ... HOPKINS: exactly for this purpose; to

allow aggressive optimisation by compilers in various ways e.g. because it is guaranteed to be alias-free. In theory, this should make it potentially faster than ANSI C until the restrict keyword is implemented, should it not?

BUDGE: Yes, that was the idea. I wanted *valarray* to provide a mechanism for expressing any one-loop operation as a single expression which could be highly optimized. I also had a vague notion that nested-loop expressions could in turn be expressed as single expressions on nested template classes, but the experience just wasn't there to see all the implications -- you should know that *valarray* was originally *not* a class template, but a pair of classes based on int and double for which there *was* some experience. This is because implementations of templates were not widely available at the time *valarray* was first proposed. HOPKINS:

However, there does seem to be a body of (informed?) opinion that *valarray* is 'broken' or at least not working well enough to be worth the effort.

BUDGE: Yes, that is probably a fair assessment, and probably the assessment of the vast majority (though not a unanimous view.) It has become fairly clear that the aliasing guarantees provided by *valarray* simply aren't strong enough to be useful, and that the market incentives for taking advantage of them aren't strong enough even if they were. *valarray* was written at a time when vector supercomputers were still the sexy leading edge of computing. Unfortunately, the best optimization strategy for a vector supercomputer is almost the opposite of the best optimization strategy for modern hierarchical-storage machines. On a vector supercomputer, you wanted to run the largest possible data set past each instruction, so that the vector pipeline remained full. On a hierarchical-memory machine, you want to throw the largest possible number of instructions at a particular working set of data, so that you keep your data in cache (or paged into memory or on processor, depending on which level of the memory hierarchy you are concerned with.) *valarray* might conceivably have been helpful for optimization on vector machines, because it assumes operations are best treated atomically. It's hopeless on modern machines. In any case, it's not at all clear that *valarray* is the right philosophy. *valarray* was meant to replace loops with expressions, but *STL* has shown that loops can be beautiful. HOPKINS:

Does anyone here have experience of the runtime speed of a carefully constructed (i.e. avoiding the C++ efficiency pitfalls) *valarray*-based set of BLAS, LAPACK or whatever?

BUDGE: I'm not aware of anyone attempting to implement BLAS or LAPACK using *valarray*. HOPKINS: If so, are you using

slices, iterators etc and doing it fully *STL*-style or doing access with more traditional Fortran-style (i,j) operators.

Compiler writers; are you taking full advantage of *valarray*? Does it offer what it suggests?

BUDGE: Arch Robison can answer this better than I, but the short answer to both questions is No. HOPKINS:

Numerical methods people; have you compared it (on a level playing field) with C, Fortran, C++ template-based libraries (Blitz, MTL, POOMA)?

Chuck

Charles R Harris wrote:

...

On 9/10/07, *Christopher Barker* <Chris.Barker@noaa.gov STL either, so I'm not sure there is any downside to valarray. It looks like neither one [vector or valarray] supports any kind of "view" semantics, so for the purposes of numpy array wrapping, they really aren't any different.

I think that the originator of valarray saying it was misguided might be considered a downside.

I had read that, though interestingly, I haven't seen any more recent commentary about the issues at all.

In any case, it appears that what Budge is saying is that the original goal of valarray being well used for optimized numerical routines isn't going to happen (I don't think it has, though there is a PPC altivec version out there). However std::vector doesn't have any numerical optimizations either, so I don't see any reason to choose std::vector over std:valarray.

My real question is what compiler and library writers are doing -- has anyone (OK, I guess MS and gcc are all I care about anyway) built anything optimized for them? Are they going to dump them? Who knows? What do you mean by optimization ? I think this question is the key. I remember having used blitz at some point, and I thought it was terrible. It is really complicated, and to get good performances was really difficult. Maybe I used it wrongly, I don't know (this was a few years ago). But at some point, I decided to just use plain C arrays instead:

Christopher Barker wrote: the code was much faster, and actually much easier (I really hate template syntax). I personnally don't think all the template things worth it for optimizing temporaries (which was the goal of blitz): the complexity cost is enormous, for not much benefit. I think C++ is much more useful for the automatic memory management through RAII, which is what std::vector gives you. As long as you think about setting the right size to avoid resizes, all other considerations are not worthwhile IMHO. C speed is already quite good on modern CPU, and std::vector gives you that. If your compiler supports restrict, use it (http://www.cellperformance.com/mike_acton/2006/05/demystifying_the_restrict_...), this will give you "Fortran speed". The fact that, while C++ being a popular language, a standard class for matrix algebra does not exist yet shows me that this is not that useful, or too complicate to develop. cheers, David

Christopher Barker wrote:

David Cournapeau wrote:

...

Christopher Barker wrote:

...

My real question is what compiler and library writers are doing -- has anyone (OK, I guess MS and gcc are all I care about anyway) built anything optimized for them? Are they going to dump them? Who knows? What do you mean by optimization ?

Well, I'm quite specifically not being precise about that. It appears the POINT of valarray was to provide a way to do computation that compiler(library) writers could optimize in various ways for the system at hand. The one example I have seen is someone that wrote a version that takes advantage of the PPC altivec instructions:

(http://www.pixelglow.com/stories/altivec-valarray-2/)

Anyway, at this point I'm far less concerned about optimization that just a more robust and convenient way to deal with data that raw pointers.

...

I remember having used blitz at some point, and I thought it was terrible.

Darn -- it looks so promising. I realize that I sounded more convinced than I really am. First, to make my perspective more obvious, let me say that I generally hate template. I think the syntax is terrible, and make the code totally unreadable for everything but simple cases (simple container, for example); I think it is a wrong solution for a broken language. So I prefer to avoid them if I can.

My understanding of blitz is that it is supposed to be faster mainly because it can avoid temporaries thanks to expression template. So if you don't need this feature, you don't gain much. But when you think about it, avoiding temporaries is done by symbolic computation at the compiler level through template; the idea is to make expressions such as A = B * C + D * E^-1 * F where everything is a matrix the most efficient possible. C/C++ makes it hard because it needs to use binary operations with a returned value. So in the end, this is really a parsing problem; if so, why not use a language which can do symbolic computation, and convert them into a compiled language ? By using expression template, you use a totally broken syntax to do things which are much more easily done by a language easy to parse (say LISP). When you take a look at http://ubiety.uwaterloo.ca/~tveldhui/papers/DrDobbs2/drdobbs2.html, you also realize that the tests are done on architectures/compilers which are different from the ones available now. The only way to really know is to do your own tests: have a reasonable example of the kind of operations you intend to do, benchmark it, and see the differences. My experience says it definitely does not worth it for my problems. Maybe yours will be different.

...

I think C++ is much more useful for the automatic memory management through RAII, which is what std::vector gives you.

and std::valarray not? I guess where I'm at now is deciding if there is any advantage or disadvantage to using std::valarray vs. std::vector. The other option is to go with something else: boost::multiarray, blitz++, etc. However, at least in term of how well they might p;lay with numpy arrays, I don't see a reason to do so.

Valarray and vector give you more or less the same here concerning RAII. But vector really is more common, so I would rather pick up vector instead of valarray unless there is a good reason not to do so, not the contrary. I don't know much boost::multiarray (I tried a bit ublas, and found the performances quite bad compared to C, using gcc; again, this was a few years ago, it may have changed since). I almost never used more than rank 2 arrays, so I don't know much about multi_array. cheers, David

On 9/17/07, David Cournapeau <david@ar.media.kyoto-u.ac.jp> wrote:

...

David Cournapeau wrote:

...

Christopher Barker wrote:

...

My real question is what compiler and library writers are doing -- has anyone (OK, I guess MS and gcc are all I care about anyway) built anything optimized for them? Are they going to dump them? Who knows? What do you mean by optimization ?

Well, I'm quite specifically not being precise about that. It appears the POINT of valarray was to provide a way to do computation that compiler(library) writers could optimize in various ways for the system at hand. The one example I have seen is someone that wrote a version that takes advantage of the PPC altivec instructions:

(http://www.pixelglow.com/stories/altivec-valarray-2/)

Anyway, at this point I'm far less concerned about optimization that just a more robust and convenient way to deal with data that raw

Christopher Barker wrote: pointers.

...

...

I remember having used blitz at some point, and I thought it was

terrible.

...

Darn -- it looks so promising.

I realize that I sounded more convinced than I really am. First, to make my perspective more obvious, let me say that I generally hate template. I think the syntax is terrible, and make the code totally unreadable for everything but simple cases (simple container, for example); I think it is a wrong solution for a broken language. So I prefer to avoid them if I can.

My understanding of blitz is that it is supposed to be faster mainly because it can avoid temporaries thanks to expression template. So if you don't need this feature, you don't gain much. But when you think about it, avoiding temporaries is done by symbolic computation at the compiler level through template; the idea is to make expressions such as A = B * C + D * E^-1 * F where everything is a matrix the most efficient possible. C/C++ makes it hard because it needs to use binary operations with a returned value. So in the end, this is really a parsing problem; if so, why not use a language which can do symbolic computation, and convert them into a compiled language ? By using expression template, you use a totally broken syntax to do things which are much more easily done by a language easy to parse (say LISP).

Templates are a godsend for containers and such using multiple types. I think that much of Numpy could be naturally written up that way. Template programming, on the other hand, seems to me an attempt to use the template mechanism as a compiler. So you are probably right that a different language would handle that more easily. When you take a look at

http://ubiety.uwaterloo.ca/~tveldhui/papers/DrDobbs2/drdobbs2.html, you also realize that the tests are done on architectures/compilers which are different from the ones available now.

The only way to really know is to do your own tests: have a reasonable example of the kind of operations you intend to do, benchmark it, and see the differences. My experience says it definitely does not worth it for my problems. Maybe yours will be different.

...

...

I think C++ is much more useful for the automatic memory management through RAII, which is what std::vector gives you.

and std::valarray not? I guess where I'm at now is deciding if there is any advantage or disadvantage to using std::valarray vs. std::vector. The other option is to go with something else: boost::multiarray, blitz++, etc. However, at least in term of how well they might p;lay with numpy arrays, I don't see a reason to do so.

Valarray and vector give you more or less the same here concerning RAII. But vector really is more common, so I would rather pick up vector

That is my general feeling too, valarrays are the red-headed stepchildren of the stl. instead of valarray unless there is a good reason not to do so, not the

contrary. I don't know much boost::multiarray (I tried a bit ublas, and found the performances quite bad compared to C, using gcc; again, this was a few years ago, it may have changed since). I almost never used more than rank 2 arrays, so I don't know much about multi_array.

I found the performance of ublas to be pretty good for small arrays when it was compiled with the -NODEBUG option, the assembly code looked pretty good too. The default with all the bounds checking and such is terrible and the assembly code practically unreadable (180 character function identifiers, etc), but for debugging it did its job. The main virtue of ublas is compactness and readability of code expression, which is far better than writing out endless numbers of loops. The automatic handling of pointers for the default allocation type is also convenient and makes it reasonable to have functions return matrices and vectors. I still think FORTRAN might be a better choice than C++ for these sort of problems, it is just that C++ has become the default for (too) many things. Chuck

"Charles R Harris" <charlesr.harris@gmail.com> writes:

The automatic handling of pointers for the default allocation type is also convenient and makes it reasonable to have functions return matrices and vectors.

Hmm, I wonder whether I missed something when I read the manual. I didn't see anything in the docs that suggests that ublas matrices do COW, reference semantics or anything else to make C++'s horrible pass-by-value semantics bearable performancewise, so I return and pass in shared_ptr's to matrices, which is syntactically ugly but avoids the need to write a (reference semantics) wrapper class for matrix. Am I missing some easier way to efficiently return and pass large matrices? 'as

David Cournapeau <david@ar.media.kyoto-u.ac.jp> writes:

Alexander Schmolck wrote:

...

"Charles R Harris" <charlesr.harris@gmail.com> writes:

...

The automatic handling of pointers for the default allocation type is also convenient and makes it reasonable to have functions return matrices and vectors.

Hmm, I wonder whether I missed something when I read the manual. I didn't see anything in the docs that suggests that ublas matrices do COW, reference semantics or anything else to make C++'s horrible pass-by-value semantics bearable performancewise, so I return and pass in shared_ptr's to matrices, which is syntactically ugly but avoids the need to write a (reference semantics) wrapper class for matrix. Am I missing some easier way to efficiently return and pass large matrices?

If ublas is using expression template, shouldn't it alleviate somewhat this problem ?

I don't think so, but then I'm hardly a C++ whizz. As far as I can tell the point of expression tempaltes is just to provide syntactic sugar so that one can write fairly complex in-place computations as a normal mathematical expression. But say I want to pass a big matrix of datapoints to a classifier -- how would expression templates help here? Ublas does have various view objects, but they're of limited usefulness, because they don't provide the same functionality as the matrix class itself. 'as

"Bill Baxter" <wbaxter@gmail.com> writes:

...

But say I want to pass a big matrix of datapoints to a classifier -- how would expression templates help here? Ublas does have various view objects, but they're of limited usefulness, because they don't provide the same functionality as the matrix class itself.

What's wrong with using references?

void my_classifier(BigMatrix& datapoints) { ... }

Nothing, at least as long as its a const reference (I'm not sure the syntactic pleasantness of non-const references are worth the additional semantics-obscurification). But since this won't work for return values and since having life-time management sounded useful I just pretty much went with shared_ptr throughout. If even some semi-manual pseudo-gc scheme, like "smart" pointers will essentially be zero performance overhead and no additional headaches, because there are a few big matrices and no danger of cycic structures, there didn't seem much reason not to use it. I'd just have preferred it if I could have had this (i.e. memory managed reference semantics) wrapped up more nicely without doing it myself (like in Qt's container classes, which if I'm not mistaken either are reference semantics or COW). 'as

Thanks! DG Lisandro Dalcin wrote:

David, I'll try to show you what I do for a custom C++ class, of course this does not solve the issue resizing (my class does not actually support resizing, so this is fine for me):

My custom class is a templatized one called DTable (is like a 2d contiguous array), but currently I only instantiate it for 'int' and 'double'. This class have two relevant methods: getShape(), returning a std::pair, and getArray(), returning a reference to an underliing std::vector. So I am able to automatically support array interface by using this:

First, I define a templatized utility function 'array_interface'

%header %{ namespace numpy {

template<typename T> static char typechar() { return '\0'; } template<> static char typechar<int>() { return 'i'; } template<> static char typechar<double>() { return 'f'; }

template<typename T> static PyObject* array_interface(const DTable<T>* self) { const std::pair<int,int>& shape = self->getShape(); const std::vector<T>& data = self->getArray(); void* array = const_cast<T*>(&data[0]); char endian = PyArray_NATIVE; char kind = typechar<T>(); int elsize = sizeof(T); return Py_BuildValue("{sNsNsNsN}", "shape", Py_BuildValue("ii", shape.first, shape.second), "typestr", PyString_FromFormat("%c%c%d", endian, kind, elsize), "data", Py_BuildValue("NO", PyLong_FromVoidPtr(array), Py_False), "version", PyInt_FromLong(3)); } } %}

Now define a SWIG macro to apply it to instantiations of my class

%define %array_interface(Class) %extend Class { PyObject* __array_interface__; } %{ #define %mangle(Class) ##_## __array_interface__ ## _get(_t) \ numpy::array_interface(_t) #define %mangle(Class) ##_## __array_interface__ ## _set(_t, _val) \ SWIG_exception_fail(SWIG_AttributeError, "read-only attribute") %} %enddef

and finally instantiate my class with different names for 'int' and 'double' in SWIG and finally apply previous macro

%template(DTableI) DTable<int>; %template(DTableS) DTable<double>;

%array_interface( DTable<int> ); %array_interface( DTable<double> );

I think you can implement someting similar for std::vector, std::valarray, or whatever... For use other data types, the only you need is to specialize the 'typecode' function.

Hope this help you.

On 9/4/07, David Goldsmith <David.L.Goldsmith@noaa.gov> wrote:

...

Anyone have a well-tested SWIG-based C++ STL valarray <=> numpy.array typemap to share? Thanks!

DG -- ERD/ORR/NOS/NOAA <http://response.restoration.noaa.gov/emergencyresponse/> _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion