utpm algorithm class methods do not organize UTPM correctly

[mikofski]

Given the example:

import numpy as np
from algopy import UTPM, zeros
import numdifftools.nd_algopy as nda
import numdifftools as nd


def f(x):
    nobs = x.shape[1:]
    f0 = x[0]**2 * np.sin(x[1])**2
    f1 = x[0]**2 * np.cos(x[1])**2
    out = zeros((2,) + nobs, dtype=x)
    out[0,:] = f0
    out[1,:] = f1
    return out

x = np.array([(1, 2, 3, 4), (5, 6, 7, 8)], dtype=float)
y = f(x)

xj = UTPM.init_jacobian(x)
j = UTPM.extract_jacobian(f(xj))

returns this traceback

---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
<ipython-input-31-a606a9700adf> in <module>()
     19
     20 xj = UTPM.init_jacobian(x)
---> 21 j = UTPM.extract_jacobian(f(xj))
     22
     23 print "x =\n%r\n" % x

<ipython-input-31-a606a9700adf> in f(x)
      8     nobs = x.shape[1:]
      9     print "nobs: %r" % nobs
---> 10     f0 = x[0]**2 * np.sin(x[1])**2
     11     f1 = x[0]**2 * np.cos(x[1])**2
     12     out = zeros((2,) + nobs, dtype=x)

c:\python27\lib\site-packages\algopy\utpm\utpm.pyc in __mul__(self, rhs)
    364             else:
    365                 err_str = 'binary operations between UTPM instances and object arrays are not supported'
--> 366                 raise NotImplementedError(err_str)
    367
    368         elif isinstance(rhs,numpy.ndarray):

NotImplementedError: binary operations between UTPM instances and object arrays are not supported

but the output is not an object array it is an array of UTPM, which has gotten unpacked somehow:

rhs
array([ UTPM([[ 0.91953576  0.91953576  0.91953576  0.91953576  0.91953576  0.91953576   0.91953576  0.91953576]
              [ 0.          0.          0.          0.         -0.54402111  0.          0.   0.        ]]),
        UTPM([[ 0.07807302  0.07807302  0.07807302  0.07807302  0.07807302  0.07807302   0.07807302  0.07807302]
              [ 0.          0.          0.          0.          0.         -0.53657292   0.          0.        ]]),
        UTPM([[ 0.43163139  0.43163139  0.43163139  0.43163139  0.43163139  0.43163139   0.43163139  0.43163139]
              [ 0.          0.          0.          0.          0.          0.   0.99060736  0.        ]]),
        UTPM([[ 0.97882974  0.97882974  0.97882974  0.97882974  0.97882974  0.97882974   0.97882974  0.97882974]
              [ 0.          0.          0.          0.          0.          0.          0.  -0.28790332]])],
    dtype=object)

So I need to figure out how to monkey patch algorithms.py to keep these together.

[argriffing]

the output is not an object array it is an array of UTPM

I'm not sure what's going on in this issue, but it shows some confusion about object arrays. In your example, the rhs array has dtype=object so it's an object array.

[mikofski]

I think this is related to what happens when you pass an array to an overloaded ufunc, see NumPy issue 7519.

Also see @b45ch1 comment on the AlgoPy google groups.

In this case numpy.sin() is overloaded in utpm.py:631

This in turn calls algorithms.py:915.

For some reason, the outer UTPM container of x is being treated as an iterable, so NumPy is calling utpm.sin on each item, and creating an object array of UTPM.

Calling using UTPM.sin methods seems to fix the problem.

import numpy as np
from algopy import UTPM, zeros, sin, cos
import numdifftools.nd_algopy as nda
import numdifftools as nd


def f(x):
    nobs = x.shape[1:]
    f0 = x[0]**2 * sin(x[1])**2  # use UTPM version of sin
    f1 = x[0]**2 * cos(x[1])**2  # use UTPM version of sin
    out = zeros((2,) + nobs, dtype=x)
    out[0,:] = f0
    out[1,:] = f1
    return out

x = np.array([(1, 2, 3, 4), (5, 6, 7, 8)], dtype=float)
y = f(x)

xj = UTPM.init_jacobian(x)
j = UTPM.extract_jacobian(f(xj))

But it would be awesome if AlgoPy could be used without having to change the code. It should be possible to use UTPM.as_utpm() to fix this before doing other operations.