0.3.1, int64 and possible related error

[Michel]

Hi Sebastian,

I've just discovered this mailing after sending you the Email. Sorry
about that.

I want to use alogpy to prototype algorithms involving AD.

My small example (I'm still not fluent in Python):

import numpy, algopy, scipy
from algopy import UTPM, exp
from scipy.sparse.linalg import cg
from numpy import array

x1 = UTPM.init_jacobian([1,0])
x2 = UTPM.init_jacobian([0,0])
b0 = UTPM.init_jacobian([0,0])

x1=[4,1]
x2=[1,3]

A=array([x1,x2])
b0=array([1,2])

z=cg(A,b0)

print 'z: ', z[0]

res=z[0]

algopy_jacobian = UTPM.extract_jacobian(res[0])
print 'jacobian = ',algopy_jacobian

[Michel]

Great, hit the send button ;-).

Well anyway... the output of my little script is:

Traceback (most recent call last):
File "./getting2.py", line 22, in <module>

algopy_jacobian = UTPM.extract_jacobian(res[0])

File "/usr/local/lib/python2.6/dist-packages/algopy-0.3.0-py2.6.egg/
algopy/utpm/utpm.py", line 983, in extract_jacobian
return x.data[1,...].transpose([i for i in
range(1,x.data[1,...].ndim)] + [0])
TypeError: sequence index must be integer

In your fixes for version 0.3.1 you mention:

> fixed a bug in getitem of the UTPM instance: now works also with numpy.int64 as index

Could this be the reason for my error? If so, where do I get version
0.3.1? I cloned your git repo but AFAIK it is version 0.3.0.

If not, I will try to resolve the issue.

Michel

[Sebastian Walter]

Hello Michel,

On Mon, Jul 4, 2011 at 1:30 PM, Michel <Michel....@rwth-aachen.de> wrote:
> Great, hit the send button ;-).
>
> Well anyway... the output of my little script is:
>
> Traceback (most recent call last):
>  File "./getting2.py", line 22, in <module>
>    algopy_jacobian = UTPM.extract_jacobian(res[0])
>  File "/usr/local/lib/python2.6/dist-packages/algopy-0.3.0-py2.6.egg/
> algopy/utpm/utpm.py", line 983, in extract_jacobian
>    return x.data[1,...].transpose([i for i in
> range(1,x.data[1,...].ndim)] + [0])
> TypeError: sequence index must be integer
>
> In your fixes for version 0.3.1 you mention:
>
>> fixed a bug in getitem of the UTPM instance: now works also with numpy.int64 as index
>
> Could this be the reason for my error? If so, where do I get version
> 0.3.1? I cloned your git repo but AFAIK it is version 0.3.0.

I have uploaded version 0.3.1 to pypi now:
http://pypi.python.org/pypi/algopy/0.3.1
On github you can find the bleeding edge version (i.e. 0.3.1.dev).

>
> If not, I will try to resolve the issue.

AFAIK scipy.sparse.linalg.cg calls a Fortran routine.
This is a problem for AlgoPy (and any other AD tool in Python) since
it relies on Python's operator overloading capabilities.

>
> Michel
>
> On Jul 4, 1:27 pm, Michel <Michel.Scha...@rwth-aachen.de> wrote:
>> Hi Sebastian,
>>
>> I've just discovered this mailing after sending you the Email. Sorry
>> about that.
>>
>> I want to use alogpy to prototype algorithms involving AD.
>>
>> My small example (I'm still not fluent in Python):
>>
>> import numpy, algopy, scipy
>> from algopy import UTPM, exp
>> from scipy.sparse.linalg import cg
>> from numpy import array
>>
>> x1 = UTPM.init_jacobian([1,0])
>> x2 = UTPM.init_jacobian([0,0])
>> b0 = UTPM.init_jacobian([0,0])
>>
>> x1=[4,1]
>> x2=[1,3]
>>
>> A=array([x1,x2])
>> b0=array([1,2])
>>
>> z=cg(A,b0)
>>
>> print 'z: ', z[0]

Here is the problem:

In [23]: z[0].dtype
Out[23]: dtype('float64')

but it should be a UTPM instance. SciPy probably just converts A and
b0 by calling float(A) and float(b0), neglecting all derivative
information thereby.

You have three options:
1) Implement the cg algorithm in pure Python. Then all operations in
the cg algorithm can be performed in univariate Taylor polynomial
arithmetic.
2) Don't use AD but divided differences. For some quick prototyping
this may be the quickest way to get some results.
3) Add a specialized routine algopy.sparse.linalg.cg, e.g. by using
the implicit function theorem.

The third option would be the most convenient in the long run, I think.

>>
>> res=z[0]
>>
>> algopy_jacobian = UTPM.extract_jacobian(res[0])
>> print 'jacobian = ',algopy_jacobian

Hope this helps.

cheers,
Sebastian

[Michel]

On Jul 5, 11:20 am, Sebastian Walter <sebastian.wal...@gmail.com>
wrote:

> Hello Michel,
>
>
>
>
>
>
>
>
>
> On Mon, Jul 4, 2011 at 1:30 PM, Michel <Michel.Scha...@rwth-aachen.de> wrote:
> > Great, hit the send button ;-).
>
> > Well anyway... the output of my little script is:
>
> > Traceback (most recent call last):
> >  File "./getting2.py", line 22, in <module>
> >    algopy_jacobian = UTPM.extract_jacobian(res[0])
> >  File "/usr/local/lib/python2.6/dist-packages/algopy-0.3.0-py2.6.egg/
> > algopy/utpm/utpm.py", line 983, in extract_jacobian
> >    return x.data[1,...].transpose([i for i in
> > range(1,x.data[1,...].ndim)] + [0])
> > TypeError: sequence index must be integer
>
> > In your fixes for version 0.3.1 you mention:
>
> >> fixed a bug in getitem of the UTPM instance: now works also with numpy.int64 as index
>
> > Could this be the reason for my error? If so, where do I get version
> > 0.3.1? I cloned your git repo but AFAIK it is version 0.3.0.
>
> I have uploaded version 0.3.1 to pypi now:http://pypi.python.org/pypi/algopy/0.3.1
> On github you can find the bleeding edge version (i.e.  0.3.1.dev).
>
>

Thank you very much! I gave it a try right away.

>
> > If not, I will try to resolve the issue.
>
> AFAIK scipy.sparse.linalg.cg calls a Fortran routine.
> This is a problem for AlgoPy (and any other AD tool in Python) since
> it relies on Python's operator overloading capabilities.
>
>
>

Ok. So the same holds true for pyadol?

Since I will use algopy mainly for prototyping new ideas/algoritms, I
will probably avoid manual adaptation of a given code. I actually
implemented a simple cg solver. Here is my code (the comments are for
the passive run):

import numpy, algopy, scipy
from algopy import UTPM, exp

from numpy import array,inner,dot

def cg(A,b):
x=[0.0,0.0]
p=r=b-dot(A,x)
a=inner(r,r)
while (a>1e-9):
v=dot(A,p)
l=a/(inner(v,p))
x_old=x
x=x+l*p
r=r-l*v
a_old=a
a=inner(r,r)
p=r+(a/a_old)*p
return [x,x_old]

b = UTPM(numpy.zeros((2,1,2)))
b.data[0,0,0]=1
b.data[1,0,0]=0
b.data[0,0,1]=2
b.data[1,0,1]=0

A = UTPM(numpy.zeros((2,1,2,2)))
A.data[0,0,0,0]=4
A.data[1,0,0,0]=1
A.data[0,0,1,0]=1
A.data[1,0,1,0]=0
A.data[0,0,0,1]=1
A.data[1,0,0,1]=0
A.data[0,0,1,1]=3
A.data[1,0,1,1]=0
#A=[[4.0,1.0],[1.0,3.0]]
#b=[2.0,1.0]

print 'A: ' , A
print 'b: ' , b
x=cg(A,b)

print 'x: ', x
print 'Result: ', dot(A,x[0])

>
>
> >> res=z[0]
>
> >> algopy_jacobian = UTPM.extract_jacobian(res[0])
> >> print 'jacobian = ',algopy_jacobian
>
> Hope this helps.

It sure does. Here is exception I get while running my cg solver using
algopy:

Traceback (most recent call last):

File "./getting2.py", line 40, in <module>
x=cg(A,b)
File "./getting2.py", line 7, in cg
p=r=b-dot(A,x)
File "/usr/local/lib/python2.6/dist-packages/algopy-0.3.1-py2.6.egg/
algopy/utpm/utpm.py", line 252, in __sub__
raise NotImplementedError('should implement that')
NotImplementedError: should implement that

As it mentions "should implement that" I thought I might bring it up
to you. I tried to rearrange the code, but nothing solved this issue.
Do I have to reimplement the matrix product too because it may rely on
other languages?

>
> cheers,
> Sebastian

Thank you for you help.

[Sebastian Walter]

Hi Michel,

also for other AD tools like funcdesigner and CAS like casadi,Theano,....

AlgoPy requires "clean" code. There are certain bad practices that are
not allowed or highly discouraged.
If you have clean code, then it's pretty simple to apply AlgoPy (or
any other AD tool).

I have modified your cg algorithm a little:

--------------- start code -------------------------
import numpy, algopy
from algopy import UTPM, exp, dot
from numpy import array

def cg(A,b,x):
p=r=b-dot(A,x)
a=dot(r,r)
while (a>1e-9):
v=dot(A,p)
l=a/(dot(v,p))
x_old=x
x=x+l*p
r=r-l*v
a_old=a
a=dot(r,r)
p=r+(a/a_old)*p
return [x,x_old]

b = UTPM(numpy.zeros((2,1,2)))
b.data[0,0,0]=1
b.data[1,0,0]=0
b.data[0,0,1]=2
b.data[1,0,1]=0

A = UTPM(numpy.zeros((2,1,2,2)))
A.data[0,0,0,0]=4
A.data[1,0,0,0]=1
A.data[0,0,1,0]=1
A.data[1,0,1,0]=0
A.data[0,0,0,1]=1
A.data[1,0,0,1]=0
A.data[0,0,1,1]=3
A.data[1,0,1,1]=0

x = UTPM(numpy.zeros((2,1,2)))

#A=[[4.0,1.0],[1.0,3.0]]
#b=[2.0,1.0]

print 'A: ' , A
print 'b: ' , b

x=cg(A,b,x)[0]

print 'x: ', x

print 'solution correct?', dot(A,x) - b

--------------- end code -------------------------
It seems it's working correctly.

>
>
>>
>>
>> >> res=z[0]
>>
>> >> algopy_jacobian = UTPM.extract_jacobian(res[0])
>> >> print 'jacobian = ',algopy_jacobian
>>
>> Hope this helps.
>
> It sure does. Here is exception I get while running my cg solver using
> algopy:
>
> Traceback (most recent call last):
>  File "./getting2.py", line 40, in <module>
>    x=cg(A,b)
>  File "./getting2.py", line 7, in cg
>    p=r=b-dot(A,x)
>  File "/usr/local/lib/python2.6/dist-packages/algopy-0.3.1-py2.6.egg/
> algopy/utpm/utpm.py", line 252, in __sub__
>    raise NotImplementedError('should implement that')
> NotImplementedError: should implement that

In numpy you can write things like

x = [1,2,3]
y = numpy.array([1.,2.,3.])

z = x + y

and numpy automatically converts the list to a numpy array and then
adds the two arrays together.
This is a feature that isn't supported in algopy.

I have also never felt the need to use numpy.inner sind numpy.dot
worked just fine.
Since Python has the "there should be only one obvious way" philosophy
I refrained to support inner unless someone convinces me that it's
absolutely necessary ;)

>
> As it mentions "should implement that" I thought I might bring it up
> to you. I tried to rearrange the code, but nothing solved this issue.
> Do I have to reimplement the matrix product too because it may rely on
> other languages?
>

No, dot is directly supported in AlgoPy.

cheers,
Sebastian

[Michel]

Thank you very much for your help and for clearing things up!

This example and all my other little implementations work now.

On Jul 8, 10:16 am, Sebastian Walter <sebastian.wal...@gmail.com>
wrote: