Using the Autodif module

The autodif module is designed to make gradient computation trivial for any mathematical function of an arbitrary number of parameters.

To use the autodif module from python, first we have to import it

In [2]:
import godot.core.autodif as ad
import numpy as np

We can construct an autodif scalar using only a floating point value

In [3]:
x = ad.Scalar( 1.5 )
print( x.value() )
print( x )
1.5
         Value
  1.500000e+00

But to track the gradient with respect to this we need to set track to true

In [4]:
x = ad.Scalar( 1.5, True )
print( x )
         Value |    Leaf ID: 1
  1.500000e+00 |  1.000000e+00

or, alternatively give it a name

In [5]:
y = ad.Scalar( 1.5, "y" )
print( y )
         Value |             y
  1.500000e+00 |  1.000000e+00

In the print out, the value and the gradient with respect to tracked parameters is given.

The scalars support many operations. These operations compute the value and the gradients, keeping track of dependent parameters.

In [5]:
x*x
Out[5]:
         Value |    Leaf ID: 1
  2.250000e+00 |  3.000000e+00
In [6]:
ad.sin(x)
Out[6]:
         Value |    Leaf ID: 1
  9.974950e-01 |  7.073720e-02

We can do all kids of things, all the time tracking the partials. For example

In [7]:
def cheby(w):
    t = [0]*10
    t[0] = 1
    t[1] = w
    for i in range(2,10):
        t[i] = 2 * w * t[i-1] - t[i-2]
    return t

c = cheby(x)

for ci in c:
    print(ci)
1
         Value |    Leaf ID: 1
  1.500000e+00 |  1.000000e+00

         Value |    Leaf ID: 1
  3.500000e+00 |  6.000000e+00

         Value |    Leaf ID: 1
  9.000000e+00 |  2.400000e+01

         Value |    Leaf ID: 1
  2.350000e+01 |  8.400000e+01

         Value |    Leaf ID: 1
  6.150000e+01 |  2.750000e+02

         Value |    Leaf ID: 1
  1.610000e+02 |  8.640000e+02

         Value |    Leaf ID: 1
  4.215000e+02 |  2.639000e+03

         Value |    Leaf ID: 1
  1.103500e+03 |  7.896000e+03

         Value |    Leaf ID: 1
  2.889000e+03 |  2.325600e+04

We can redefine x so that we see the partial name and then track more than one parameter

In [8]:
x = ad.Scalar( 2.5, "x")

z = y/x

c = cheby(z)

for ci in c:
    print(ci)
1
         Value |             y |             x
  6.000000e-01 |  4.000000e-01 | -2.400000e-01

         Value |             y |             x
 -2.800000e-01 |  9.600000e-01 | -5.760000e-01

         Value |             y |             x
 -9.360000e-01 |  5.280000e-01 | -3.168000e-01

         Value |             y |             x
 -8.432000e-01 | -1.075200e+00 |  6.451200e-01

         Value |             y |             x
 -7.584000e-02 | -2.492800e+00 |  1.495680e+00

         Value |             y |             x
  7.521920e-01 | -1.976832e+00 |  1.186099e+00

         Value |             y |             x
  9.784704e-01 |  7.223552e-01 | -4.334131e-01

         Value |             y |             x
  4.219725e-01 |  3.626435e+00 | -2.175861e+00

         Value |             y |             x
 -4.721034e-01 |  3.966944e+00 | -2.380167e+00

We can retrieve the gradient from the autodif scalar

In [9]:
grad = z.gradient()
grad
Out[9]:
             y             x
  4.000000e-01 -2.400000e-01

or like this for a names list of parameters:

In [10]:
grad = z.at(['y','x'])
grad
Out[10]:
array([ 0.4 , -0.24])

And we can get a specific element as a numpy array as long as an autodif scalar is tracking only one parameter (leaf).

In [11]:
ad.getPartial( z, x.leaf() )
Out[11]:
array([-0.24])
In [12]:
ad.getPartial( z*c[3], x.leaf() )
Out[12]:
array([0.03456])

But this is not okay, because z has 2 elements of its gradient vector and therefore 2 leaves.

In [13]:
try:
    ad.getPartial( z*c[3], z.leaf() )
except RuntimeError as err:
    print("This is not okay ", err)
This is not okay  
Exception raised in GODOT: 
What : Pre-condition check failed: size() == 1. Help message: Number of leaves must be equal to one
Func : leaf
Src  : /builds/godot/godotpy/godot_prefix/include/godot/core/autodif/xdouble.h
Line : 177

For that we can do this

In [14]:
z.at(x.leaf())
Out[14]:
array([-0.24])

Vectors

We can use vectors too. For example:

In [15]:
y = ad.Vector( [ 1, 2, 3 ], "y" )
print( y.value() )
print( y )
[1. 2. 3.]
         Value |          y[0]          y[1]          y[2]
  1.000000e+00 |  1.000000e+00  0.000000e+00  0.000000e+00
  2.000000e+00 |  0.000000e+00  1.000000e+00  0.000000e+00
  3.000000e+00 |  0.000000e+00  0.000000e+00  1.000000e+00

many operations act like you would expect. Here are a few examples.

In [16]:
y[0]
Out[16]:
         Value |          y[0]          y[1]          y[2]
  1.000000e+00 |  1.000000e+00  0.000000e+00  0.000000e+00
In [17]:
y[1:]
Out[17]:
         Value |          y[0]          y[1]          y[2]
  2.000000e+00 |  0.000000e+00  1.000000e+00  0.000000e+00
  3.000000e+00 |  0.000000e+00  0.000000e+00  1.000000e+00
In [18]:
y[1:2]
Out[18]:
         Value |          y[0]          y[1]          y[2]
  2.000000e+00 |  0.000000e+00  1.000000e+00  0.000000e+00
In [19]:
x = ad.Vector( [ 4, 5, 6 ], "x" )

z = ad.cross( x, y ) + y /3

print(z)
         Value |          y[0]          y[1]          y[2] |          x[0]          x[1]          x[2]
  3.333333e+00 |  3.333333e-01 -6.000000e+00  5.000000e+00 |  0.000000e+00  3.000000e+00 -2.000000e+00
 -5.333333e+00 |  6.000000e+00  3.333333e-01 -4.000000e+00 | -3.000000e+00  0.000000e+00  1.000000e+00
  4.000000e+00 | -5.000000e+00  4.000000e+00  3.333333e-01 |  2.000000e+00 -1.000000e+00  0.000000e+00

Numerical Integration

Suppose we have a a simple vector differential equation:

$$ \dot{\bf{x}} = f\left( \bf{x} \right) = \frac{\bf{x}}{\bf{x} \cdot \bf{x}} $$
In [20]:
def f(t, x ):
    return x / ( ad.dot(x,x) )

We use an Runge Kutte 4th order method to integrate this from an initial condition given as an autodif vector

In [21]:
def rk4(t, x, f, h):
    k1 = f(t, x)
    k2 = f(t+h/2, x+h*k1/2)
    k3 = f(t+h/2, x+h*k2/2)
    k4 = f(t+h, x+h*k3)
    return h * (k1 + k2 + k3 + k4) / 6
In [22]:
h = 1
"initialise x"
x = ad.Vector( [ 4, 5, 6 ], "x0" )
for t in np.arange(1,10,h):
    x = rk4( t, x, f, h)
    print(x)
         Value |         x0[0]         x0[1]         x0[2]
  3.441038e-02 |  5.050190e-03 -4.440506e-03 -5.328607e-03
  4.301297e-02 | -4.440506e-03  3.051963e-03 -6.660758e-03
  5.161557e-02 | -5.328607e-03 -6.660758e-03  6.096845e-04

         Value |         x0[0]         x0[1]         x0[2]
  1.535132e+00 |  3.801919e-01 -4.488964e-03 -5.386757e-03
  1.918915e+00 | -4.488964e-03  3.781718e-01 -6.733446e-03
  2.302698e+00 | -5.386757e-03 -6.733446e-03  3.757029e-01

         Value |         x0[0]         x0[1]         x0[2]
  8.661230e-02 |  1.318559e-02 -1.058435e-02 -1.270122e-02
  1.082654e-01 | -1.058435e-02  8.422637e-03 -1.587653e-02
  1.299185e-01 | -1.270122e-02 -1.587653e-02  2.601244e-03

         Value |         x0[0]         x0[1]         x0[2]
  6.597487e-01 |  1.558125e-01 -1.140583e-02 -1.368700e-02
  8.246858e-01 | -1.140583e-02  1.506799e-01 -1.710875e-02
  9.896230e-01 | -1.368700e-02 -1.710875e-02  1.444067e-01

         Value |         x0[0]         x0[1]         x0[2]
  1.763699e-01 |  2.997816e-02 -1.764290e-02 -2.117148e-02
  2.204624e-01 | -1.764290e-02  2.203885e-02 -2.646436e-02
  2.645549e-01 | -2.117148e-02 -2.646436e-02  1.233526e-02

         Value |         x0[0]         x0[1]         x0[2]
  3.936913e-01 |  8.480312e-02 -1.702463e-02 -2.042955e-02
  4.921141e-01 | -1.702463e-02  7.714203e-02 -2.553694e-02
  5.905369e-01 | -2.042955e-02 -2.553694e-02  6.777849e-02

         Value |         x0[0]         x0[1]         x0[2]
  2.483640e-01 |  4.661822e-02 -1.934097e-02 -2.320916e-02
  3.104550e-01 | -1.934097e-02  3.791479e-02 -2.901145e-02
  3.725460e-01 | -2.320916e-02 -2.901145e-02  2.727725e-02

         Value |         x0[0]         x0[1]         x0[2]
  3.224894e-01 |  6.573519e-02 -1.860894e-02 -2.233073e-02
  4.031117e-01 | -1.860894e-02  5.736117e-02 -2.791341e-02
  4.837341e-01 | -2.233073e-02 -2.791341e-02  4.712625e-02

         Value |         x0[0]         x0[1]         x0[2]
  2.787621e-01 |  5.429909e-02 -1.923930e-02 -2.308716e-02
  3.484527e-01 | -1.923930e-02  4.564140e-02 -2.885895e-02
  4.181432e-01 | -2.308716e-02 -2.885895e-02  3.505979e-02

However, if we tried to integrate the same rhs with non-autodif variables we would run into problems because the rhs uses the ad.dot function.

To solve this issue, we introduced the bridge module:

In [23]:
import godot.core.autodif.bridge as br

This allows us to write functions which can be either numpy OR autodif

In [24]:
def f( t, x ):
    return x / ( br.dot(x,x) )
In [25]:
h = 1
x = ad.Vector( [ 4, 5, 6 ], "x0" )

for t in np.arange(1,10,h):
    x = rk4(t,x,f,h)
    print(x)
         Value |         x0[0]         x0[1]         x0[2]
  3.441038e-02 |  5.050190e-03 -4.440506e-03 -5.328607e-03
  4.301297e-02 | -4.440506e-03  3.051963e-03 -6.660758e-03
  5.161557e-02 | -5.328607e-03 -6.660758e-03  6.096845e-04

         Value |         x0[0]         x0[1]         x0[2]
  1.535132e+00 |  3.801919e-01 -4.488964e-03 -5.386757e-03
  1.918915e+00 | -4.488964e-03  3.781718e-01 -6.733446e-03
  2.302698e+00 | -5.386757e-03 -6.733446e-03  3.757029e-01

         Value |         x0[0]         x0[1]         x0[2]
  8.661230e-02 |  1.318559e-02 -1.058435e-02 -1.270122e-02
  1.082654e-01 | -1.058435e-02  8.422637e-03 -1.587653e-02
  1.299185e-01 | -1.270122e-02 -1.587653e-02  2.601244e-03

         Value |         x0[0]         x0[1]         x0[2]
  6.597487e-01 |  1.558125e-01 -1.140583e-02 -1.368700e-02
  8.246858e-01 | -1.140583e-02  1.506799e-01 -1.710875e-02
  9.896230e-01 | -1.368700e-02 -1.710875e-02  1.444067e-01

         Value |         x0[0]         x0[1]         x0[2]
  1.763699e-01 |  2.997816e-02 -1.764290e-02 -2.117148e-02
  2.204624e-01 | -1.764290e-02  2.203885e-02 -2.646436e-02
  2.645549e-01 | -2.117148e-02 -2.646436e-02  1.233526e-02

         Value |         x0[0]         x0[1]         x0[2]
  3.936913e-01 |  8.480312e-02 -1.702463e-02 -2.042955e-02
  4.921141e-01 | -1.702463e-02  7.714203e-02 -2.553694e-02
  5.905369e-01 | -2.042955e-02 -2.553694e-02  6.777849e-02

         Value |         x0[0]         x0[1]         x0[2]
  2.483640e-01 |  4.661822e-02 -1.934097e-02 -2.320916e-02
  3.104550e-01 | -1.934097e-02  3.791479e-02 -2.901145e-02
  3.725460e-01 | -2.320916e-02 -2.901145e-02  2.727725e-02

         Value |         x0[0]         x0[1]         x0[2]
  3.224894e-01 |  6.573519e-02 -1.860894e-02 -2.233073e-02
  4.031117e-01 | -1.860894e-02  5.736117e-02 -2.791341e-02
  4.837341e-01 | -2.233073e-02 -2.791341e-02  4.712625e-02

         Value |         x0[0]         x0[1]         x0[2]
  2.787621e-01 |  5.429909e-02 -1.923930e-02 -2.308716e-02
  3.484527e-01 | -1.923930e-02  4.564140e-02 -2.885895e-02
  4.181432e-01 | -2.308716e-02 -2.885895e-02  3.505979e-02

In [26]:
h = 1
x = [ 4, 5, 6 ]

for t in np.arange(1,10,h):
    x = rk4(t,x,f,h)
    print(x)
[0.03441038 0.04301297 0.05161557]
[1.53513217 1.91891521 2.30269825]
[0.0866123  0.10826538 0.12991845]
[0.65974867 0.82468584 0.989623  ]
[0.17636993 0.22046241 0.26455489]
[0.39369126 0.49211408 0.5905369 ]
[0.24836399 0.31045499 0.37254599]
[0.32248938 0.40311173 0.48373407]
[0.27876213 0.34845267 0.4181432 ]

Matrices

It is possible to create matrix parameters, each element is then tracked separately as a separate parameter.

In [27]:
Z = ad.Matrix( [[ 1, 2, 3 ],[4,5,6],[7,8,9] ], "Z" )

print( Z*y )
         Value |          y[0]          y[1]          y[2] |          Z[0]          Z[1]          Z[2]          Z[3]          Z[4]          Z[5]          Z[6]          Z[7]          Z[8]
  1.400000e+01 |  1.000000e+00  2.000000e+00  3.000000e+00 |  1.000000e+00  0.000000e+00  0.000000e+00  2.000000e+00  0.000000e+00  0.000000e+00  3.000000e+00  0.000000e+00  0.000000e+00
  3.200000e+01 |  4.000000e+00  5.000000e+00  6.000000e+00 |  0.000000e+00  1.000000e+00  0.000000e+00  0.000000e+00  2.000000e+00  0.000000e+00  0.000000e+00  3.000000e+00  0.000000e+00
  5.000000e+01 |  7.000000e+00  8.000000e+00  9.000000e+00 |  0.000000e+00  0.000000e+00  1.000000e+00  0.000000e+00  0.000000e+00  2.000000e+00  0.000000e+00  0.000000e+00  3.000000e+00

But you can also construct matrices from Scalars.

Suppose we have small rotations dx, dy and dzabout the x, y and z axes

In [28]:
dx = ad.Scalar(0.0,"dx")
dy = ad.Scalar(0.0,"dy")
dz = ad.Scalar(0.0,"dz")

We define the rotation matrices assuming the small angle approximation

In [29]:
import numpy as np

Rx = ad.Matrix(np.identity(3))
Rx[1,2] = -dx
Rx[2,1] = +dx

Ry = ad.Matrix(np.identity(3))
Ry[0,2] = +dy
Ry[2,0] = -dy

Rz = ad.Matrix(np.identity(3))
Rz[0,1] = +dz
Rz[1,0] = -dz

print(Rx)
print(Ry)
print(Rz)
         Value |            dx
  1.000000e+00 |              
  0.000000e+00 |              
  0.000000e+00 |              
  0.000000e+00 |              
  1.000000e+00 |              
  0.000000e+00 |  1.000000e+00
  0.000000e+00 |              
 -0.000000e+00 | -1.000000e+00
  1.000000e+00 |              

         Value |            dy
  1.000000e+00 |              
  0.000000e+00 |              
 -0.000000e+00 | -1.000000e+00
  0.000000e+00 |              
  1.000000e+00 |              
  0.000000e+00 |              
  0.000000e+00 |  1.000000e+00
  0.000000e+00 |              
  1.000000e+00 |              

         Value |            dz
  1.000000e+00 |              
 -0.000000e+00 | -1.000000e+00
  0.000000e+00 |              
  0.000000e+00 |  1.000000e+00
  1.000000e+00 |              
  0.000000e+00 |              
  0.000000e+00 |              
  0.000000e+00 |              
  1.000000e+00 |              

and we can combine them to get a rotation matrix, tracking partials with respect to the small rotation angles.

In [30]:
R = Rz * Ry * Rx
print(R)
         Value |            dx |            dy |            dz
  1.000000e+00 |               |  0.000000e+00 |  0.000000e+00
  0.000000e+00 |               |  0.000000e+00 | -1.000000e+00
  0.000000e+00 |               | -1.000000e+00 |              
  0.000000e+00 |  0.000000e+00 |  0.000000e+00 |  1.000000e+00
  1.000000e+00 |  0.000000e+00 |  0.000000e+00 |  0.000000e+00
  0.000000e+00 |  1.000000e+00 |  0.000000e+00 |              
  0.000000e+00 |  0.000000e+00 |  1.000000e+00 |  0.000000e+00
  0.000000e+00 | -1.000000e+00 |  0.000000e+00 |  0.000000e+00
  1.000000e+00 |  0.000000e+00 |  0.000000e+00 |              

In [ ]: