[C++-sig] [Numpy-discussion] Overloading sqrt(5.5)*myvector

[Bruce Sherwood]

I found by timing measurements that a faster scheme with less penalty 
for the case of sqrt(array) looks like this:

def sqrt(x):
     if type(x) is float: return mathsqrt(x)
     return numpysqrt(x)

Bruce Sherwood wrote:
> Roman Yakovenko wrote:
>> On Dec 29, 2007 7:47 AM, Bruce Sherwood <Bruce_Sherwood at ncsu.edu> wrote:
>>  
>>> I realized belatedly that I should upgrade from Boost 1.33 to 1.34.
>>> Alas, that didn't cure my problem.
>>>     
>> Can you post small and complete example of what you are trying to 
>> achieve?
>>   
> I don't have a "small and complete" example available, but I'll 
> summarize from earlier posts. VPython (vpython.org) has its own vector 
> class to mimic the properties of 3D vectors used in physics, in the 
> service of easy creation of 3D animations. There is a beta version 
> which imports numpy and uses it internally; the older production 
> version uses Numeric. Boost python and thread libraries are used to 
> connect the C++ VPython code to Python.
>
> There is operator overloading that includes scalar*vector and 
> vector*scalar, both producing vector. With Numeric, sqrt produced a 
> float, which was a scalar for the operator overloading. With numpy, 
> sqrt produces a numpy.float64 which is caught by vector*scalar but not 
> by scalar*vector, which means that scalar*vector produces an ndarray 
> rather than a vector, which leads to a big performance hit in existing 
> VPython programs. The overloading and Boost code is the same in the 
> VPython/Numeric and VPython/numpy versions. I don't know whether the 
> problem is with numpy or with Boost or with the combination of the two.
>
> Here is the relevant part of the vector class:
>
> inline vector
> operator*( const double s) const throw()
> { return vector( s*x, s*y, s*z); }
>
> and here is the free function for right multiplication:
>
> inline vector
> operator*( const double& s, const vector& v)
> { return vector( s*v.x, s*v.y, s*v.z); }
>
> Maybe the problem is in the Boost definitions:
>
> py::class_<vector>("vector", py::init< py::optional<double, double, 
> double> >())
>    .def( self * double())
>    .def( double() * self)
>
> Left multiplication is fine, but right multiplication isn't.
>
> A colleague suggested the following Boost declarations but cautioned 
> that he wasn't sure of the syntax for referring to operator, and 
> indeed this doesn't compile:
>
> .def( "__mul__", &vector::operator*(double), "Multiply vector times 
> scalar")
> .def( "__rmul__", &operator*(const double&, const vector&), "Multiply 
> scalar times vector")
>
> I would really appreciate a Boost or numpy expert being able to tell 
> me what's wrong (if anything) with these forms. However, I may have a 
> useful workaround as I described in a post to the numpy discussion 
> list. A colleague suggested that I do something like this for sqrt and 
> other such mathematical functions:
>
> def sqrt(x):
>   try: return mathsqrt(x)
>   except TypeError: return numpysqrt(x)
>
> That is, first try the simple case of a scalar argument, handled by 
> the math module sqrt, and only use the numpy sqrt routine in the case 
> of an array argument. Even with the overhead of the try/except 
> machinery, one gets must faster square roots for scalar arguments this 
> way than with the numpy sqrt.
>
> Bruce Sherwood
>