Checkout the code from the subversion repository:
~/downloads$ svn co http://gr.anu.edu.au/svn/people/sdburton/pydx
This will create a directory called pydx.
Option 1 (recomended) is to build using make:
~/downloads$ cd pydx ~/downloads/pydx$ make
This will build the modules in-place. (Note that make produces lots of warnings: this is OK).
In order for python to always find pydx you will need to add the pydx subdirectory to your PYTHONPATH:
~/downloads/pydx$ export PYTHONPATH=$PYTHONPATH:`pwd`
Option 2 for installing is with the setup.py script:
~/downloads/pydx$ python setup.py build ~/downloads/pydx$ sudo python setup.py install
The test suite is designed to be driven from py.test:
~/downloads/pydx$ py.test pydx ...
The python float type uses the native c compiler's double type. This is typically a 64 bit number, with 52 bits for the mantissa, 11 bits for the (signed) exponent, and one sign bit.
>>> x=1.0123456789 >>> x 1.0123456789 >>> x/11 0.092031425354545462
We use gmpy, a python wrapper around the GNU multiprecision library, gmp.
ref: http://swox.com/gmp/
Exponent size is fixed, usually 16 bits. The mantissa can be arbitrarily large, 64 bits by default.
>>> from gmpy import mpf
>>> x=mpf(1.0123456789)
>>> x
mpf('1.01234567890000004553e0')
>>> x/11
mpf('9.20314253545454586846e-2')
The pydx.scalar.mpfi module provides an Interval class based on the MPFI library.
>>> from pydx.scalar.mpfi import Interval
>>> x = Interval(1,2)
>>> x**2
Interval(mpf('1.0000000000000000e0'), mpf('4.0000000000000000e0'))
The Interval class is a python c-extension. It builds upon libgmp (the GNU multiprecision library), libmpf, which provides exact rounding semantics of the multi-precision gmp numbers, and libmpfi which implements an interval type (a pair of mpf numbers).
This is an iterative procedure. We apply the same operation to a set of float's and a corresponding set of Interval's. At each iteration we assert that the interval contains its corresponding float.
The pydx.mjet module provides multivariate automatic differentiation to arbitrary order.
For example, to find the first derivative of lambda x:x**2 at x=3.0 we do the following:
>>> from pydx.mjet import Jet >>> f = lambda x:x**2 >>> x = Jet([3.0, 1.0]) >>> f(x)[1] 6.0
The componenents of a Jet are the (normalized) derivatives of the variable. So the first component is the value of that variable (the zeroth derivative), and the second component is the first derivative. The parameter in this case is x itself, so the second component of the x Jet is 1.0.
Jet's are built lazily, ie. components are computed only when requested. In the above example, f(x) produces the lazy square of a Jet:
>>> f(x) SqrJet(<0:1.0 1:1.0>)
Multivariate Jet's are indexed using a tuple of indices. Given two (independant) parameters, x and y, with values 2.0 and 3.0 respectively, their Jet's are defined as:
>>> x = Jet({(0,0):2.0, (1,0):1.0})
>>> y = Jet({(0,0):3.0, (0,1):1.0})
To determine the first derivative of sin(x**2+y) with respect to x:
>>> from pydx.scalar.fmath import sin >>> sin(x**2+y)[1,0] 3.01560901737
Once again, the lazy object looks like this:
>>> sin(x**2+y) SinJet(AddJet(SqrJet(<00:2.0 10:1.0>),<00:3.0 01:1.0>))
Sometimes a calculation involves rank-0 Jet's. These are essentially scalar values (constants). The value is accessed with the empty tuple ():
>>> x=Jet(rank=0) >>> x <:0.0> >>> x[()] 0.0 >>> x[()]=5.0 >>> x[()] 5.0
Most of the classes in pydx subscribe to an underlying scalar type. For example, if we wish to compute Jet's using arbitrary precision floating point types, we set the mpf scalar:
>>> from pydx.scalar import set_mpf_scalar >>> set_mpf_scalar(64) # use a 64 bit precision
Now the components we enter for a Jet are promoted to the type of the current scalar:
>>> x=Jet({(0,0):2.0,(1,0):1.0})
>>> x
<00:mpf('2.e0',64) 10:mpf('1.e0',64)>
Scalar's are stored on a stack, so when we are done we can restore the scalar type:
>>> from pydx.scalar import restore_scalar >>> restore_scalar()
We can perform the arithmetic operations lazily with the symbolic scalars:
>>> from pydx.scalar import set_symbolic_scalar
>>> from pydx.scalar.symbolic import Var
>>> from pydx.scalar.fmath import exp
>>> set_symbolic_scalar()
>>> x = Jet([Var('x'),1.0])
>>> exp(x)
ExpJet(<0:x 1:<pydx.scalar.symbolic.OneFloat object at 0xb7b155ac>>)
exp(x) is now a (lazy) Jet with lazy components. When we ask for, say, the third derivative, we get a symbolic object back:
>>> exp(x)[3]
<pydx.scalar.symbolic.DivFloat object at 0xb7b1558c>
>>> exp(s)[3].deepstr()
DivFloat(DivFloat(ExpFloat(Var('x')), Float(2.0)), Float(3.0))
To actually do anything useful with this object, we can get a python expression:
>>> dddx = exp(x)[3] >>> dddx.expr() '((x.exp() / 2.0) / 3.0)' >>> restore_scalar()
For big calculations it is often useful to create a python function from the symbolic calculation. For example, this can make (repeated) evaluation of a large tensor calculation many thousands of times faster.
>>> set_symbolic_scalar()
>>> ctx = ComputationContext('func', ['x'], open('temp.py','w'))
>>> ctx.assign('dddx', exp(x)[3])
>>> func = ctx.finalize('dddx')
>>> restore_scalar()
The symbolic scalar's use the current scalar type to represent "constants".
Now we have a python function we can use to calculate exp(x)[3]:
>>> from pydx.scalar import set_interval_scalar
>>> set_interval_scalar()
>>> from pydx.scalar.mpfi import Interval
>>> func(Interval(12.0))
Interval(mpf('2.7125798569833983e4'), mpf('2.7125798569833990e4'))
>>> restore_scalar()
All the Jet object types derive from a "Meta-Jet" abstract base class. These objects store their rank, which is the dimensionality of the Jet (how many free variables we are differentiating with respect to). The current scalar is set as a class variable.
class MJet(object):
scalar = None
def __init__(self, rank):
self.rank = rank
...
An MJet does not have a way of accessing its components (this is left for subclasses to implement), but every MJet does know how to participate in algebraic operations. Crucial to this is being able to "promote" various other objects to the same status as self. For example, scalars need to be promoted to a "constant" Jet.
class MJet(object):
...
def promote(self, item):
if not isinstance(item,MJet):
item = self.scalar.promote(item)
x = Jet({(0,)*self.rank:item})
return x
elif item.rank == self.rank:
return item
...
This piece of code uses a "concrete" subclass of MJet, one whose components are directly stored, called simply a Jet.
class Jet(MJet):
""" This is a concrete Jet.
All components are stored in a dictionary.
"""
def __init__(self, cs={}, rank=None):
...
MJet.__init__(self, rank)
self.cs = {} # coefficients, components, ...
for idxs, value in cs.items():
self[idxs] = self.scalar.promote(value)
Upon creation of a Jet, the components are promoted to the current scalar type.
The bulk of the MJet implementation is made up of the special methods that implement various algebraic operations. The key here is that no computation of components takes place. Instead, a lazy object is created that knows how to compute its components from those of its children. For example, this method implements the addition of two MJet's.
class MJet(object):
...
def __add__(self, other):
return AddJet(self, other)
All the lazy MJet's inherit from a JetOp class that manages a cache of components. This is so that when components are repeatedly requested we can eliminate recomputing them. A speed-up of around 30 times was observed by using this cache strategy.
class JetOp(MJet):
""" Lazy Jet with item cache.
"""
def __init__(self, rank):
MJet.__init__(self, rank)
self._cache = {} # speeds up things by *30
The __getitem__ for a JetOp merely examines the cache and then calls self.getitem (implemented in the subclass) if the value is not found.
The meat of the actual computation of components is found in these getitem methods:
class AddJet(JetOp):
def __init__(self, a, b):
JetOp.__init__(self, a.rank)
b = self.promote(b)
self.a = a
self.b = b
def getitem(self, idxs):
return self.a[idxs]+self.b[idxs]
Here we see how AddJet requests components from it's "children" and then adds them.
A more tricky example comes from the lazy multiplication JetOp:
class MulJet(JetOp):
...
def getitem(self, idxs):
assert type(idxs)==tuple
idxss = [ _ for _ in multi_range(idxs) ]
_idxss = [ multi_sub(idxs,_idxs) for _idxs in idxss ]
return sum((self.a[a_idxs] * self.b[b_idxs]
for a_idxs, b_idxs in zip(_idxss,idxss)), self.scalar.zero)
This time we need to do some manipulation of the multi-index objects. The multi_range function yields multi-indices starting from the zero multi-index going up to the argument (inclusive):
def multi_range(idxs):
if len(idxs):
for idx in range(0,idxs[0]+1):
for rest in multi_range(idxs[1:]):
yield (idx,)+rest
else:
yield ()
multi_sub implements subtraction of multi-indices:
def multi_sub(a_idxs, b_idxs):
c_idxs = tuple(a_idxs[i]-b_idxs[i] for i in range(len(a_idxs)))
return c_idxs
All the tests involve using the MJet machinery to compute coefficients of Taylor polynomials of functions. The polynomials are then evaluated at different points, and compared to the original function they approximate.
We use two ways of verifying this comparison. The first uses regular python floating point scalars and an assertion of a bound on the (relative) error of the taylor polynomial. The second uses the exact form of taylor's theorem, along with the Interval scalars, and asserts that the value of the polynomial overlaps with the value of the function.
The Interval method verifies that we "do no wrong", while the floating-point method verifies that we are "approximately right".
Since an MJet can be interpreted as the coefficents of a Taylor polynomial, evaluation of such a polynomial is implemented as an expand method of the MJet class. Here we employ this method in asserting a relative error bound:
def test_jet_univariate_simple():
f = lambda x: 1.0/(x+1.0)
order = 5
x0 = random() - 0.5
x = Jet([x0, 1.0]) # Concrete MJet
h = 0.01*random() # small float
fx = f(x) # MJet
fxh = f(x0+h) # float
fxe = fx.expand((h,), order) # float
error = abs(fxe - fxh)
relative_error = error / max(1e-6, error)
assert relative_error < 1e-10
The bound 1e-10 is chosen in an ad-hoc fashion, as is the range of the h variable.
The Interval test is more involved (even more so for the multivariate tests). This time we produce an error Interval and pass it to the MJet method expand_err:
def test_jet_expand_2():
f = lambda x: 1.0/(x+1.0)
order = 5
set_interval_scalar()
x0 = Interval(random()-0.5)
x = Jet([x0,1.0])
h = Interval(random()-0.5)
xx0 = x0.hull(x0+h)
xx = Jet([xx0,1.0])
err = f(xx)[order+1]
fx = f(x)
fxh = f(x0+h)
fxe = fx.expand_err((h,), order, (err,))
assert fxh.overlapping(fxe)
restore_scalar()
This is a much more satisfying test; there is no particularly ad-hoc choice of bounds, the Interval class takes care of all of that.
The pydx.tensor.Tensor class defines a multi-dimensional array type with scalar components. Mathematically, they store the value of some tensor (at a point), in some coordinate system.
To build a Tensor object we specify a valence (a tuple of contra- or co-variant specifiers), and the dimension:
>>> x = Tensor((dn,), 4) # a four dimensional covariant vector >>> x[0] = 1.0
Tensors live in a vector space, so we can add/subtract them, and scalar multiply.
TensorField objects are conceived of as functions (callables) that return a Tensor object.
Here we define the metric of a two dimensional flat-space:
>>> from pydx.field import TensorField >>> metric = TensorField.identity((dn,dn), 2)
Calling the metric will compute its value at a point:
>>> metric(1.2, 3.4) 1.0 0.0 0.0 1.0
We can perform algebraic operations on TensorField's, this is another lazy operation. (It corresponds to the usual mathematical pointwise operation.)
For example, let us define a concrete scalar field by specifying the component functions:
>>> from pydx.field import ConcreteTensorField >>> r = ConcreteTensorField((), 2) >>> r[()] = lambda x,y:x**2+y**2
In this case there is only one component function (indexed with the empty tuple).
Multiplying with the metric produces a new TensorField object:
>>> rmetric = r*metric >>> rmetric(1.2,3.4) 13.0 0.0 0.0 13.0
PyDX is licensed under the GNU GPL license.