Robert Kern wrote:

2009/1/20 Neal Becker :

...

I tried a little experiment, implementing some code in numpy (usually I build modules in c++ to interface to python). Since these operations are all large vectors, I hoped it would be reasonably efficient.

The code in question is simple. It is a model of an amplifier, modeled by it's AM/AM and AM/PM characteristics.

The function in question is the __call__ operator. The test program plots a spectrum, calling this operator 1024 times each time with a vector of 4096.

If you want to find out what lines in that function are taking the most time, you can try my line_profiler module:

http://www.enthought.com/~rkern/cgi-bin/hgwebdir.cgi/line_profiler/

That might give us a better idea in the absence of a self-contained example.

Sounds interesting, I'll give that a try. But, not sure how to use it. If my main script is plot_spectrum.py, and I want to profile the ampl.__call__ function (defined in ampl.py), what do I need to do? I tried running kernprof.py plot_spectrum.py having added @profile decorators in ampl.py, but that didn't work: File "../mod/ampl.py", line 43, in ampl @profile NameError: name 'profile' is not defined

Robert Kern wrote:

2009/1/20 Neal Becker :

...

I tried a little experiment, implementing some code in numpy (usually I build modules in c++ to interface to python). Since these operations are all large vectors, I hoped it would be reasonably efficient.

The code in question is simple. It is a model of an amplifier, modeled by it's AM/AM and AM/PM characteristics.

The function in question is the __call__ operator. The test program plots a spectrum, calling this operator 1024 times each time with a vector of 4096.

If you want to find out what lines in that function are taking the most time, you can try my line_profiler module:

http://www.enthought.com/~rkern/cgi-bin/hgwebdir.cgi/line_profiler/

That might give us a better idea in the absence of a self-contained example.

I see the problem. Thanks for the great profiler! You ought to make this more widely known. It seems the big chunks of time are used in data conversion between numpy and my own vectors classes. Mine are wrappers around boost::ublas. The conversion must be falling back on a very inefficient method since there is no special code to handle numpy vectors. Not sure what is the best solution. It would be _great_ if I could make boost::python objects that export a buffer interface, but I have absolutely no idea how to do this (and so far noone else has volunteered any info on this).

Robert Kern wrote:

On Tue, Jan 20, 2009 at 20:57, Neal Becker wrote:

...

I see the problem. Thanks for the great profiler! You ought to make this more widely known.

I'll be making a release shortly.

...

It seems the big chunks of time are used in data conversion between numpy and my own vectors classes. Mine are wrappers around boost::ublas. The conversion must be falling back on a very inefficient method since there is no special code to handle numpy vectors.

Not sure what is the best solution. It would be _great_ if I could make boost::python objects that export a buffer interface, but I have absolutely no idea how to do this (and so far noone else has volunteered any info on this).

Who's not volunteering information, boost::python or us?

I've asked on python.c++, the home of boost::python and friends. I've spent a little time on it myself, but I think this job requires great knowledge of python c api as well as the mysteries of boost::python.

T J wrote:

On Tue, Jan 20, 2009 at 6:57 PM, Neal Becker wrote:

...

It seems the big chunks of time are used in data conversion between numpy and my own vectors classes. Mine are wrappers around boost::ublas. The conversion must be falling back on a very inefficient method since there is no special code to handle numpy vectors.

Not sure what is the best solution. It would be _great_ if I could make boost::python objects that export a buffer interface, but I have absolutely no idea how to do this (and so far noone else has volunteered any info on this).

I'm not sure if I've understood everything here, but I think that pyublas provides exactly what you need.

http://tiker.net/doc/pyublas/

It might if I had used this for all of my c++ code, but I have a big library of c++ wrapped code that doesn't use pyublas. Pyublas takes numpy objects from python and allows the use of c++ ublas on it (without conversion). Most of my code doesn't use numpy, it uses plain ublas to represent vectors, and ublas handles storage. I can only interface to/from numpy with conversion. I'm interested in pyublas, but devel seems very quiet for a while.

On 1/21/2009 2:38 PM, Sturla Molden wrote:

If you can get a pointer (as integer) to your C++ data, and the shape and dtype is known, you may use this (rather unsafe) 'fromaddress' hack:

And opposite, if you need to get the address referenced to by an ndarray, you can do this: def addressof(array): return arr.__array_interface__['data'][0] Then you will have to cast this unsigned integer to a pointer type in C++. Note that arr.data returns a buffer. Sturla Molden

Ravi wrote:

Hi Neal,

On Wednesday 21 January 2009 07:27:04 Neal Becker wrote:

...

It might if I had used this for all of my c++ code, but I have a big library of c++ wrapped code that doesn't use pyublas.  Pyublas takes numpy objects from python and allows the use of c++ ublas on it (without conversion).

Most of my code doesn't use numpy, it uses plain ublas to represent vectors, and ublas handles storage.  I can only interface to/from numpy with conversion.

I pointed out my code to you on c++-sig[1] a while back that solves precisely this problem. You found a bug with memory management that I fixed in the updated code. Does that still not work for you?

Regards, Ravi

[1] [http://mail.python.org/pipermail/cplusplus-sig/2008-October/013825.html Thanks for reminding me about this!

Do you have a current version of the code? I grabbed the files from the above message, but I see some additional subsequent messages with more patches.

Ravi wrote:

On Wednesday 21 January 2009 10:22:36 Neal Becker wrote:

...

...

[http://mail.python.org/pipermail/cplusplus-sig/2008- October/013825.html

Thanks for reminding me about this!

Do you have a current version of the code?  I grabbed the files from the above message, but I see some additional subsequent messages with more patches.

That is the latest publicly posted code. Since then, there is just one minor patch (attached) which enables use of row-major (c-contiguous) arrays.

This does *not* work with strided arrays which would be a fair bit of effort to support. Further, you will have to work with the numpy iterator interface, which, while well-designed, is a great illustration of the effort required to support OO programming in an non-OO language, and is pretty tedious to map to the ublas storage iterator interface. If you do implement it, I would very much like to take a look at it.

Regards, Ravi

I'm only interested in simple strided 1-d vectors. In that case, I think your code already works. If you have c++ code using the iterator interface, the iterators dereference will use (*array )[index]. This will use operator[], which will call PyArray_GETPTR. So I think this will obey strides. Unfortunately, it will also be slow. I suggest something like the enclosed. I have done some simple tests, and it seems to work.

On Wednesday 21 January 2009 13:55:49 Neal Becker wrote:

...

I'm only interested in simple strided 1-d vectors.  In that case, I think your code already works.  If you have c++ code using the iterator interface, the iterators dereference will  use (*array )[index].  This will use operator[], which will call PyArray_GETPTR.  So I think this will obey strides.

You are right. I had forgotten that I had simple strided vectors working.

...

Unfortunately, it will also be slow.  I suggest something like the enclosed.   I have done some simple tests, and it seems to work.

I wonder why PyArray_GETPTR1 is slow. Is it because of the implied integer multiplication? Unfortunately, your approach means that iterators can become invalid if the underlying array is resized to a larger size. Hmmh, perhaps we could make this configurable at compile-time ...

Thanks for the code. Could you provide some benchmarks on the relative speeds of the two approaches?

Regards, Ravi Do you know about pyublas? This is the same issue we ran into there. I did not benchmark the code you sent me. I was just going by my experience with

Ravi wrote: pyublas. I guess a benchmark would be a good idea.

Ravi wrote:

On Wednesday 21 January 2009 10:22:36 Neal Becker wrote:

...

...

[http://mail.python.org/pipermail/cplusplus-sig/2008- October/013825.html

Thanks for reminding me about this!

Do you have a current version of the code?  I grabbed the files from the above message, but I see some additional subsequent messages with more patches.

That is the latest publicly posted code. Since then, there is just one minor patch (attached) which enables use of row-major (c-contiguous) arrays.

This does *not* work with strided arrays which would be a fair bit of effort to support. Further, you will have to work with the numpy iterator interface, which, while well-designed, is a great illustration of the effort required to support OO programming in an non-OO language, and is pretty tedious to map to the ublas storage iterator interface. If you do implement it, I would very much like to take a look at it.

Regards, Ravi

It seems your code works fine for my usual style: ublas::vector<T> func (numpy::array_from_py<T>::type const&) But not for a function that modifies it arg in-place (& instead of const&): void func (numpy::array_from_py<T>::type &) This gives: ArgumentError: Python argument types in test1.double(numpy.ndarray) did not match C++ signature: double(boost::numeric::ublas::vector > {lvalue}) My instinct is to ignore it, because I think I don't need it, but do you have a workaround?

Ravi wrote:

Hi Neal,

On Wednesday 21 January 2009 07:27:04 Neal Becker wrote:

...

It might if I had used this for all of my c++ code, but I have a big library of c++ wrapped code that doesn't use pyublas.  Pyublas takes numpy objects from python and allows the use of c++ ublas on it (without conversion).

Most of my code doesn't use numpy, it uses plain ublas to represent vectors, and ublas handles storage.  I can only interface to/from numpy with conversion.

I pointed out my code to you on c++-sig[1] a while back that solves precisely this problem. You found a bug with memory management that I fixed in the updated code. Does that still not work for you?

Regards, Ravi

[1] [http://mail.python.org/pipermail/cplusplus-sig/2008-October/013825.html

Do you know if this code will work with strided vectors? If I pass a slice: u = array (...) F (u[::2]) What happens?

Robert-- this is a great little piece of code, I already think it will be a part of my workflow. However, I seem to be getting negative % time taken on the more time consuming lines, perhaps getting some overflow? Thanks a lot, Wes On Wed, Jan 21, 2009 at 3:23 AM, Robert Kern wrote:

On Tue, Jan 20, 2009 at 20:57, Neal Becker wrote:

...

I see the problem. Thanks for the great profiler! You ought to make this more widely known.

http://pypi.python.org/pypi/line_profiler

-- Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion

I have been using your profiler extensively and it has contributed to my achieving significant improvements in the application I work on largely due to the usefulness of the line by line breakdown enabling me to easily select the next part of code to work on optimizing. So firstly many thanks for writing it. However back to my point, Wes, I have also experienced timing oddities, in particular on Virtual machines (MS Hyper-V has very poor processor timings, the older MS VM works fine though). I believe the negative timings arise when the CPU (be it virtual or possibly physical) deviates from its standard performance or rather the initial timer unit taken, would this make sense to you Robert? Hanni 2009/1/21 Robert Kern

On Wed, Jan 21, 2009 at 12:13, Wes McKinney wrote:

...

Robert-- this is a great little piece of code, I already think it will be a part of my workflow. However, I seem to be getting negative % time taken on the more time consuming lines, perhaps getting some overflow?

That's odd. Can you send me the code (perhaps offlist) or at least the .lprof file?

-- Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion

import cProfile def f(): pass def g(): for i in xrange(1000000): f() cProfile.run("g()")

test.py 1000003 function calls in 1.225 CPU seconds

Ordered by: standard name ncalls tottime percall cumtime percall filename:lineno(function) 1 0.000 0.000 1.225 1.225 <string>:1(<module>) 1000000 0.464 0.000 0.464 0.000 test.py:3(f) 1 0.761 0.761 1.225 1.225 test.py:6(g) 1 0.000 0.000 0.000 0.000 {method 'disable' of '_lsprof.Profiler' objects} Running this with line_profiler: Timer unit: 2.9485e-010 s File: test.py Function: g at line 9 Total time: 0.855075 s Line # Hits Time Per Hit % Time Line Contents ============================================================== 9 @profiler 10 def g(): 11 1000001 1844697930 1844.7 63.6 for i in xrange(1000000): 12 1000000 1055333053 1055.3 36.4 f() Which is what I would expect. Hmm On Thu, Jan 22, 2009 at 2:52 AM, Robert Kern wrote:

On Thu, Jan 22, 2009 at 01:46, Hanni Ali wrote:

...

I have been using your profiler extensively and it has contributed to my achieving significant improvements in the application I work on largely due to the usefulness of the line by line breakdown enabling me to easily select the next part of code to work on optimizing. So firstly many thanks for writing it.

My pleasure.

...

However back to my point, Wes, I have also experienced timing oddities, in particular on Virtual machines (MS Hyper-V has very poor processor timings, the older MS VM works fine though). I believe the negative timings arise when the CPU (be it virtual or possibly physical) deviates from its standard performance or rather the initial timer unit taken, would this make sense to you Robert?

Can you try using cProfile with lots of calls to empty functions? I'm using the same timer functions as cProfile.

-- Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion

Windows XP, Pentium D, Python 2.5.2 On Thu, Jan 22, 2009 at 6:03 PM, Robert Kern wrote:

On Thu, Jan 22, 2009 at 17:00, Wes McKinney wrote:

...

import cProfile

def f(): pass

def g(): for i in xrange(1000000): f()

cProfile.run("g()")

...

test.py 1000003 function calls in 1.225 CPU seconds

Ordered by: standard name

ncalls tottime percall cumtime percall filename:lineno(function) 1 0.000 0.000 1.225 1.225 <string>:1(<module>) 1000000 0.464 0.000 0.464 0.000 test.py:3(f) 1 0.761 0.761 1.225 1.225 test.py:6(g) 1 0.000 0.000 0.000 0.000 {method 'disable' of '_lsprof.Profiler' objects}

Running this with line_profiler:

Timer unit: 2.9485e-010 s

File: test.py Function: g at line 9 Total time: 0.855075 s

Line # Hits Time Per Hit % Time Line Contents ============================================================== 9 @profiler 10 def g(): 11 1000001 1844697930 1844.7 63.6 for i in xrange(1000000): 12 1000000 1055333053 1055.3 36.4 f()

Which is what I would expect. Hmm

What platform are you on?

-- Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion