[Numpy-discussion] ufunc.accumulate question

[Friedrich Romstedt]

2010/10/5 Chris Withers <chris at simplistix.co.uk>:
> Hi All,
>
> I can't find any docs on this behavior.
>
> So, I have a python function. To keep it simple, lets just do addition:
>
> def add(x,y):
>   print x,y
>   retun x+y
>
> So, I can turn this into a ufunc as follows:
>
> uadd = np.frompyfunc(add,2,1)

As a side remark, note that this will always return dtype=numpy.object
arrays.  Maybe numpy.vectorize() is more appropriate for your
use-case.

> Now, I can apply it to an array:
>
>  >>> uadd.accumulate(np.arange(3,10))
> 3 4
> 7 5
> 12 6
> 18 7
> 25 8
> 33 9
> array([3, 7, 12, 18, 25, 33, 42], dtype=object)
>
> Okay, but where did the initial 3 come from?
> http://docs.scipy.org/doc/numpy/reference/generated/numpy.ufunc.accumulate.html#numpy.ufunc.accumulate
> suggests that:
>
> r = np.empty(len(A))
> t = op.identity        # op = the ufunc being applied to A's  elements
> for i in xrange(len(A)):
>     t = op(t, A[i])
>     r[i] = t
> return r

I see from the behaviour, that the implementation makes use of the
following behaviour of the identity ``op.identity`` with respect to
``op``:

>>> op(op.identity, X) == X
True

So the implementation acts in agnosticism of ``op.identity`` by
assuming that the first call of ``t = op(t, A[i])`` in your example
will just yield ``A[0]``, while i == 0.

This is the work flow:

[ t = op.identity ]  # in brackets because it cannot be executed, it's
pseudo-code

i = 0
=====
t = op(t, A[0]) = A[0]
r[0] = t = A[0]

i = 1
=====
t = op(t, A[1]) = op(A[0], A[1])
r[1] = t

i = 2
=====
t = op(t, A[2]) = op(op(A[0], A[1]), A[2])
r[2] = t

and so on ... clear now?

Friedrich

> ...but:
>
>  >>> print uadd.identity
> None

It's simply not set, and will also not be used.

> ...and:
>
>  >>> add(None,3)
> None 3
> Traceback (most recent call last):
>   File "<stdin>", line 1, in <module>
>   File "<stdin>", line 3, in add
> TypeError: unsupported operand type(s) for +: 'NoneType' and 'int'

As said, the code in the docs is just pseudo-code, it's not the actual
(C?) implementation, it just has the same semantics.

> So, where is the reason that the 3 ends up in the output array documented?
>
> Also, what if I want to specify the identity of my newly created ufunc?
> I have a case where I want to specify it as zero:
>
>  >>> uadd.identity = 0
> Traceback (most recent call last):
>   File "<stdin>", line 1, in <module>
> AttributeError: attribute 'identity' of 'numpy.ufunc' objects is not
> writable

I think you have to write an own class for this.  You cannot subclass
numpy.ufunc ("it's not an acceptable base class" Python says).