[Numpy-discussion] wrong casting of augmented assignment statements

[Sebastian Walter]

I've written a self-contained example that shows that numpy indeed
tries to call the __float__ method.
What is buggy is what happens if calling the __float__ method raises
an Exception.
Then numpy  assumes (in this case wrongly) that the object should be
casted to the neutral element.

I'd guess that  the __float__ method is called somewhere in a try:
statement and if an exception is raised it is casted to the neutral
element.
I've tried to locate the corresponding code in the numpy sources but I
got lost. Could someone be so kind and point me to it?


-------------------- start code ----------------------

import numpy

print 'numpy.__version__ = ',numpy.__version__

class ad1:

    def __init__(self,x):
        self.x = x

    def __mul__(self,other):
        if not isinstance(other, self.__class__):
            return self.__class__(self.x * other)
        return self.__class__(self.x * other.x)

    def __rmul__(self,other):
        return self * other

    def __float__(self):
        raise Exception('this is not possible')

    def __str__(self):
        return str(self.x)

print '\nThis example yields buggy behavior:'
x1 = numpy.array([ad1(1.), ad1(2.), ad1(3.)])
y1 = numpy.random.rand(3)
print 'y1=     ',y1
print 'x1=     ',x1
z1 = x1 * y1
y1 *= x1            # this should call the __float__ method of ad1
which would raise an Exception
print 'z1=x1*y1',z1
print 'y1*=x1  ',y1

class ad2:

    def __init__(self,x):
        self.x = x

    def __mul__(self,other):
        if not isinstance(other, self.__class__):
            return self.__class__(self.x * other)
        return self.__class__(self.x * other.x)

    def __rmul__(self,other):
        return self * other

    def __float__(self):
        return float(self.x)

    def __str__(self):
        return str(self.x)

print '\nThis example works fine:'
x2 = numpy.array([ad2(1.), ad2(2.), ad2(3.)])
y2 = numpy.random.rand(3)
print 'y2=     ',y2
print 'x2=     ',x2
z2 = x2 * y2
y2 *= x2            # this should call the __float__ method of ad1
which would raise an Exception
print 'z2=x2*y2',z2
print 'y2*=x2  ',y2


-------------------- end code ----------------------

-------- output ---------
walter at wronski$ python
wrong_casting_object_to_float_of_augmented_assignment_statements.py
numpy.__version__ =  1.3.0

This example yields buggy behavior:
y1=      [ 0.15322371  0.47915903  0.81153995]
x1=      [1.0 2.0 3.0]
z1=x1*y1 [0.153223711127 0.958318053803 2.43461983729]
y1*=x1   [ 0.15322371  0.47915903  0.81153995]

This example works fine:
y2=      [ 0.49377037  0.60908423  0.79772095]
x2=      [1.0 2.0 3.0]
z2=x2*y2 [0.493770370747 1.21816846399 2.39316283707]
y2*=x2   [ 0.49377037  1.21816846  2.39316284]
-------- end  output ---------

On Tue, Jan 12, 2010 at 7:38 PM, Robert Kern <robert.kern at gmail.com> wrote:
> On Tue, Jan 12, 2010 at 12:31, Sebastian Walter
> <sebastian.walter at gmail.com> wrote:
>> On Tue, Jan 12, 2010 at 7:09 PM, Robert Kern <robert.kern at gmail.com> wrote:
>>> On Tue, Jan 12, 2010 at 12:05, Sebastian Walter
>>> <sebastian.walter at gmail.com> wrote:
>>>> Hello,
>>>> I have a question about the augmented assignment statements *=, +=, etc.
>>>> Apparently, the casting of types is not working correctly. Is this
>>>> known resp. intended behavior of numpy?
>>>
>>> Augmented assignment modifies numpy arrays in-place, so the usual
>>> casting rules for assignment into an array apply. Namely, the array
>>> being assigned into keeps its dtype.
>>
>> what are the usual casting rules?
>
> For assignment into an array, the array keeps its dtype and the data
> being assigned into it will be cast to that dtype.
>
>> How does numpy know how to cast an object to a float?
>
> For a general object, numpy will call its __float__ method.
>
>>> If you do not want in-place modification, do not use augmented assignment.
>>
>> Normally, I'd be perfectly fine with that.
>> However, this particular problem occurs when you try to automatically
>> differentiate an algorithm by using an Algorithmic Differentiation
>> (AD) tool.
>> E.g. given a function
>>
>> x = numpy.ones(2)
>> def f(x):
>>   a = numpy.ones(2)
>>   a *= x
>>   return numpy.sum(a)
>>
>> one would use an AD tool as follows:
>> x = numpy.array([adouble(1.), adouble(1.)])
>> y = f(x)
>>
>> but since the casting from object to float is not possible the
>> computed gradient \nabla_x f(x) will be wrong.
>
> Sorry, but that's just a limitation of the AD approach. There are all
> kinds of numpy constructions that AD can't handle.
>
> --
> 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 at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>