b45ch1/algopy

AlgoPy is a Research Prototype for Algorithmic Differentation in Python

algopy's Issues

Got VisibleDeprecationWarning

Got VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
S = numpy.zeros((N,M), dtype=x.dtype) warning from .\algopy\utpm\utpm.py:1798: , how do I stop this?

Python 3 compatibility

Could you please make algopy Python 3 compatible?

$ python3.2 setup.py install --user
  File "setup.py", line 86
    print git_revision
                     ^
SyntaxError: invalid syntax

Functions don't have a `dtype` attribute

I'm taping a function that constructs an array piece by piece. I does something like this:

def f(x):  # x is a Numpy array
    n = x.size
    c = numpy.empty(n+1, dtype=x.dtype)
    c[:n] = g(x)
    c[n] = whatever
    return c

I had to add the dtype=x.dtype bit for things to work with PyAdolc because when taping f, PyAdolc passes an array x of type numpy.object. Just saying c = numpy.empty(n+1) creates an array of dtype float, and then, we can't store g(x) in it.

But with AlgoPy, this breaks because x is a Function and Functions don't have a dtype attribute.

Would it make sense to add a dtype attribute to functions? For all practical purposes, it could be np.float unless specified otherwise.

Thanks.

cg.reverse([fbar]) does not check the dimensions of fbar

reduce copy and improve memory management

to build the computational graph, it is not necessary to propagate a full Mtc object through the graph.
One direction P and Degree D=0 should be ok.

should introduce a keep variable that triggers how big the allocated xbar tensors have to be, or if they even have to be allocated.

Unit tests fail

Hi!
I'm trying to bring algopy over to nixpkgs.

However, there seems to be an issue with the unit tests. Numpy moved its decorators to numpy.testing._private.decorators, so that I get some ModuleNotFoundErrors and some AttributeErrors.

I will for now only test whether importing algopy works, but a fix for the unit tests would be greatly appreciated!

runtests.py

Numpy and Scipy are now using a runtests.py script that has been distilled from pv's custom makefile to help streamline development, maybe this would also work for Algopy.
https://github.com/numpy/numpy/blob/master/runtests.py
https://github.com/scipy/scipy/blob/master/runtests.py

eigh incompatible with multidirectional UTP?

It seems that eigh is inherently incompatible with multidirectional UTP in the presence of repeated eigenvalues. The Python implementation seems correct wrt. the mathematical theorems, but can only work if the UTP is unidirectional. Here is a case reproducing the problem:

from algopy import *

D = 3
P = 2
N = 6

A = UTPM(numpy.random.random((D, P, N, N)))

# Ensure all directions have the same constant value in A (and hence later Z1)
A.data[0,:,:,:] = A.data[0,0,:,:]

# Make an orthonormal basis Z1
[Z1, _] = qr(A)

# Orthogonal basis Z1 has same constant value in all directions
numpy.testing.assert_almost_equal((Z1.data[0,:,:,:]-Z1.data[0,0,:,:]), 0)

# Make fresh eigenvalues with some repetitions
W1 = zeros((N, N), dtype=A)
for d in range(D):
    for p in range(P):
        for i in range(N):
            W1.data[d, p, i, i] = (1 + (d*(p + 1) + round(i / 2))) * 100
W1.data[0,:,:,:] = W1.data[0,0,:,:]

# Matrix A1 is the one we want to factorize into eigenvalues and vectors
A1 = dot(Z1, dot(W1, Z1.T))

# Factorize A1 as w2 and Z2
[w2, Z2] = eigh(A1)

# Z2 is indeed orthogonal as required
numpy.testing.assert_almost_equal(dot(Z2, Z2.T) - numpy.eye(N), 0)

# Reconstruction A2 of A1 via its eigendecomposition succeeds - as it should!
W2 = zeros_like(A)
W2[:N, :N] = diag(w2)
A2 = dot(Z2, dot(W2, Z2.T))
numpy.testing.assert_almost_equal(A1.data-A2.data, 0)

# But Z2 does not have the same constant value in each direction!
print("Z2 does not have the same constant value in all directions:\n\n",
        Z2.data[0,:,:,:], "\n")
numpy.testing.assert_almost_equal((Z2.data[0,:,:,:]-Z2.data[0,0,:,:]), 0)

Could you please confirm whether this issue cannot be solved for multidirectional UTP, i.e. symmetric eigendecompositions containing repeated eigenvalues can only ever work in the unidirectional case?

run_tests.py fails if scipy is not installed

Can you change run_tests.py so that it works even if scipy is not available?

$ python run_tests.py 
Traceback (most recent call last):
  File "run_tests.py", line 1, in <module>
    from algopy.tracer.tests.test_tracer import *
  File "/home/deil/software/algopy/algopy/__init__.py", line 66, in <module>
    import utpm
  File "/home/deil/software/algopy/algopy/utpm/__init__.py", line 29, in <module>
    from utpm import *
  File "/home/deil/software/algopy/algopy/utpm/utpm.py", line 18, in <module>
    from algorithms import RawAlgorithmsMixIn, broadcast_arrays_shape
  File "/home/deil/software/algopy/algopy/utpm/algorithms.py", line 34, in <module>
    from algopy import nthderiv
  File "/home/deil/software/algopy/algopy/nthderiv/__init__.py", line 5, in <module>
    from nthderiv import *
  File "/home/deil/software/algopy/algopy/nthderiv/nthderiv.py", line 339, in <module>
    @basecase(scipy.special.hyperu, domain=DOM_POS, extras=2)
AttributeError: 'NoneType' object has no attribute 'special'
[deil@localhost algopy]$ python -c 'import numpy'
[deil@localhost algopy]$ python -c 'import scipy'
Traceback (most recent call last):
  File "<string>", line 1, in <module>
ImportError: No module named scipy

please implement mixed type addition

BUG: UTPM crashes due to shapes given as float

numpy no longer accepts shapes given as floats.
Thus UTPM crashes due to shapes given as float on line 208:
https://github.com/b45ch1/algopy/blob/master/algopy/utpm/utpm.py#L208

I would propose to replace that line with these twolines:

    xsize = numpy.size(x)  # numpy.prod(x_shp)
    yr = UTPM( y.data.reshape((D,P) + (xsize,) + shp))

det() testing

I couldn't find a test for the determinant. Maybe this is because its reverse mode is not implemented? I see some pdfs on the internet where the algopy determinant has been tested in the past.

add numpy workarounds X.dot(Y) ---> dot(X,Y)

numpy.dot(X,Y) should try to call X.dot(Y) if X is of unknown type (here of Function or Mtc type)

unfortunately numpy doesn't do that.

Need to provide a workaround such that a user has only to modify
numpy.dot(X,Y) to algopy.dot(X,Y)

logdet of symmetric positive definite matrix

I'm trying to use algopy to compute this, but it only has det but not logdet built in. I see that algopy det is computed by taking a square of product of diagonals of a cholesky decomposition of a positive definite matrix. Maybe the logdet could be analogously computed by taking twice the sum of logs of the diagonal of this decomposition?

I'm not sure how the UTPM math works in this case. By the way, I cannot just use logdet := log(det) because it is less stable because the intermediate determinant might be huge even if logdet is reasonably sized.

Edit: I'm interested in this function because it is used in log likelihoods and Kullback-Leibler divergences for multivariate normal distributions.

Docstrings of ufuncs like algopy.sqrt are missing

E.g., in an ipython session

In [6]: algopy.sqrt??
Type: function
String Form:<function sqrt at 0x10413ded8>
File: Dynamically generated function. No source code available.
Definition: algopy.sqrt(_args, *_kwargs)
Docstring:

The functions
['sin','cos','tan', 'exp', 'log', 'sqrt', 'pow', 'arcsin', 'arccos', 'arctan', 'sinh', 'cosh', 'tanh', 'trace', 'zeros_like', 'diag', 'triu', 'tril', 'reshape']

are dynamically created in algopy/globalfuncs.py by means of a template string.
One should add at least the docstrings of the numpy functions to know what _args and *_kwargs are supposed to be.

conversion from (D,P,N,M) to (D,Q,P,N,M) tensors in Mtc

zeros((D,P,N,M)) should be reshaped to zeros((D,1,P,N,M))
no copying necessary

numpy dev now has generalized ufuncs for linalg

It is in the pull request numpy/numpy#3220

If I remember correctly, algopy organizes its data using ndarrays in a way that could benefit from this.

pyadolc dependency in examples

add a::

try:
    import adolc
except:
    pass

to avoid import errors scaring of new users

algopy.sum(thing, axis=-1) gives unexpected results

Here's a thing that caused me trouble when I tried to convert a numpy code to use algopy.

>>> y = np.random.randn(6 ,7)                                                   
>>> y.shape
(6, 7)
>>> np.sum(y, axis=-1).shape
(6,)

but

>>> x = algopy.UTPM(np.random.randn(4, 5, 6, 7))                                
>>> x.shape
(6, 7)
>>> algopy.sum(x, axis=-1).shape
(7,)

wishlist: utpm hessian vector product

I can get the jacobian vector product and the hessian matrix, but I don't see an efficient way to get the hessian vector product. My hope is that this could speed up trust-region conjugate gradient minimization that uses only hessian-vector products, as in https://github.com/scipy/scipy/blob/master/scipy/optimize/_trustregion_ncg.py. I am currently extracting the whole hessian and then multiplying it by the vector but I would guess that this is unnecessarily inefficient.

Cannot Concatenate Variable

I would like to use algopy with pre-built functions, but they use a concatenate operation and this seems to not work. I assume the dual number aspect of the AD object makes it unclear how to fill an array using it, but it would be very useful in terms of mapping Algopy to projects which are built using numpy if this kind of indexing were supported.

A simple working example of the issue follows:

import numpy as np
import algopy

def trace_eval_f(f, x):
    cg = algopy.CGraph()
    x = algopy.Function(x)
    y = f(x)
    cg.trace_off()
    cg.independentFunctionList = [x]
    cg.dependentFunctionList = [y]
    return cg

def func(x):
    y = algopy.ones((2, 2))
    y[1, :] = x
    out = algopy.zeros(2, dtype=x)
    out[0] = y[1, 0]**2 + y[0, 0]
    out[1] = y[1, 1]**3 + y[0, 1]
    return out

x = np.ones((2))
v = np.ones((2))
cg = trace_eval_f(func, x)
jacvec = cg.jac_vec(x, v)
jac = cg.jacobian(x)

<ipython-input-16-650cf69ba411> in func(x)
      1 def func(x):
      2     y = algopy.ones((2, 2))
----> 3     y[1, :] = x
      4     out = algopy.zeros(2, dtype=x)
      5     out[0] = x[0]**2 + y[0]

ValueError: setting an array element with a sequence.

Remove entry_point command in setup.py

The following line is present in the setup.py:

entry_points = {"distutils.commands": ["upload_sphinx = sphinx_pypi_upload:UploadDoc",]}

This is problematic because it attaches an entry point, but it's associated with the "sphinx_pypi_upload" package, which you don't provide. So now if I do python setup.py --help-commands from any python script that uses setuptools, it crashes. Just commenting out the entry_points line seems to fix it.

Support for fancy indexing

Currently, Algopy doesn't appear to support fancy indexing. Here is a simple example:

import algopy
import numpy as np

indices = [0,2]
rhs = np.array([4., 1.])

def c(x):
    "Represent the system c(x)-rhs = 0."
    cx = algopy.zeros(3, dtype=x)
    cx[0] = (1 + x[0]**2)**2 + x[1]**2
    cx[1] = algopy.sum(x)
    cx[2] = x[0] - x[1]**2
    cx[indices] -= rhs
    return cx

x0 = 2 * np.ones(2)
cg = algopy.CGraph()
x = algopy.Function(x0)
y = c(x)
cg.trace_off()
cg.independentFunctionList = [x]
cg.dependentFunctionList = [y]
print cg.jacobian(x0)

Curently, running the above script results in

Traceback (most recent call last):
  File "test_algopy.py", line 23, in <module>
    print cg.jacobian(x0)
  File "/Users/dpo/.virtualenvs/nlpy3/lib/python2.7/site-packages/algopy/tracer/tracer.py", line 357, in jacobian
    self.pushforward(utpm_x_list)
  File "/Users/dpo/.virtualenvs/nlpy3/lib/python2.7/site-packages/algopy/tracer/tracer.py", line 112, in pushforward
    raise Exception(err_str)
Exception: pushforward of node 18 failed (getitem)reported error is:
too many indicestraceback:
Traceback (most recent call last):
  File "/Users/dpo/.virtualenvs/nlpy3/lib/python2.7/site-packages/algopy/tracer/tracer.py", line 106, in pushforward
    f.__class__.pushforward(f.func, f.args, Fout = f)
  File "/Users/dpo/.virtualenvs/nlpy3/lib/python2.7/site-packages/algopy/tracer/tracer.py", line 807, in pushforward
    out  = func(*args)
  File "/Users/dpo/.virtualenvs/nlpy3/lib/python2.7/site-packages/algopy/utpm/utpm.py", line 137, in __getitem__
    tmp = self.data.__getitem__((slice(None),slice(None)) + tuple(sl))
IndexError: too many indices

utpm algorithm class methods do not organize UTPM correctly

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.

TypeError: unsupported operand type(s) for /: 'int' and 'Function'

Hi!
First, great work.

I use Algopy to solve bigger nonlinear equation systems. Therefor I need to calculate Jacobian Matrices for a gauss-newton-solver algorithm. I build the equations with Sympy and convert them into a function with lambdify.
Anyway this works pretty well with a various number of mathematical expressions like (sin,cos,tan,pow,exp ...) but crashes if there is a simple division like 1/x[0].
I am not sure what the problem is or if it really is an Algopy issue but I would be glad if someone could help me.

some examples
works:
http://ubuntuone.com/1GfIuFZzdYLEJbCHqiztPQ
doesn't work
http://ubuntuone.com/5krZDKatoUPZVgi3GwxNZx