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 --------- wal...@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.k...@gmail.com> wrote: > On Tue, Jan 12, 2010 at 12:31, Sebastian Walter > <sebastian.wal...@gmail.com> wrote: >> On Tue, Jan 12, 2010 at 7:09 PM, Robert Kern <robert.k...@gmail.com> wrote: >>> On Tue, Jan 12, 2010 at 12:05, Sebastian Walter >>> <sebastian.wal...@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@scipy.org > http://mail.scipy.org/mailman/listinfo/numpy-discussion > _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion