more features on Python symbolic math

Kragen Sitaker kragen@pobox.com
Sat, 30 Dec 2000 23:47:36 -0500 (EST)


Along with some advocacy.  :)

<kragen@pobox.com>       Kragen Sitaker     <http://www.pobox.com/~kragen/>
Perilous to all of us are the devices of an art deeper than we possess
ourselves.
       -- Gandalf the White [J.R.R. Tolkien, "The Two Towers", Bk 3, Ch. XI]

---------- Forwarded message ----------
Date: Sat, 30 Dec 2000 23:46:26 -0500 (EST)
From: Kragen <kragen@kirk.dnaco.net>
To: (name removed to protect the innocent)
Subject: symbolic math

I thought I'd engage in a little bit of advocacy, and simultaneously
get to know Python and symbolic computation a little better.  :)

Here's a bit of Python code that implements some basic algebraic
manipulation, including symbolic differentiation.  I think it
illustrates some of the advantages of higher-level languages.

- ordinary algebraic expressions can be used to build up algebraic
  expressions for manipulation; you can do this in C++ with some nifty
  template hacks Todd Veldhuizen did.
- the interpreted nature of the language makes it possible to run a
  suite of 44 regression tests every time the code is compiled (adding
  an extra 6 milliseconds of the 40-millisecond-or-so compile time)
- garbage collection means you can pass around algebraic expressions as
  easily as you can pass around integers, without worrying about memory leaks
- rock-solid memory abstraction (type safety, bounds-checking, and
  garbage collection) means that program failures are deterministic and
  close to the bug
- first-class procedures and lambda expressions make it easy to specify
  dispatch tables like opmap and derivmap
- late binding makes it possible to feed formulas into routines
  expecting numbers and have them output more formulas; however, this
  is far from perfect.  While the standard arithmetic operators can be
  and are overridden, the math library functions like sqrt() can't,
  because they're not methods of the numbers.  This is more useful when
  you want to feed virtual files, virtual network connections, or
  virtual collections to routines expecting the real thing.

Other advantages, like aggregate operations like those in STL and the
rapid feedback of an interactive interpreter, are not apparent here.

This program is not a shining example of OO design; it could benefit
from a few Visitors, I think, and it abuses isinstance() in several
places.

I don't think you could develop something like this in an afternoon in
C++ --- although you could probably do it in a week.

I'd be interested to learn if I'm overly pessimistic about C++.

# Kragen Sitaker, 2000-12-30.
# Some simple symbolic math.
# Allows you to add, subtract, etc., simple formulas -- using the ordinary
# Python operations.  Implements all of the standard arithmetic operations.
# Once you have a formula, you can print it out as a string, evaluate it
# for particular values of its variables, call 'simplified()' to get a
# simplified version, check to see if it's exactly identical to some other
# formula, or (in some cases --- if you've used only +, -, and *)
# take its symbolic derivative.  See the test() routine at the end for details.
# Symbolic derivatives will generally need to be simplified() to taste very
# good.

# Not intended for serious use; it's just a quick hack I wrote this afternoon.

# Some things it might be fun to add:
# - a compile() method that returns a Python code object that would give you
#   faster evaluation
# - a continued-fraction output mode a la HAKMEM
# - symbolic derivatives that cover more operations
# - better simplifications ( ((x + x) + (2 * x)) should simplify to (3 * x) )
# - unary operations: negation, transcendentals
# - better printing: a + b + c + d should print as a + b + c + d, not
#   as (((a + b) + c) + d)
# - other symbolic manipulations

# things inherit from Formula to get the glue that turns Python
# expressions into representations of expressions
class Formula:
    def __complex__(self): return complex(self.eval({}))
    def __int__(self): return int(self.eval({}))
    def __long__(self): return long(self.eval({}))
    def __float__(self): return float(self.eval({}))
    def __pos__(self): return self  # positive
    def __add__(self, other): return Binop('+', self, other)
    def __radd__(self, other): return Binop('+', other, self)
    def __sub__(self, other): return Binop('-', self, other)
    def __rsub__(self, other): return Binop('-', other, self)
    def __mul__(self, other): return Binop('*', self, other)
    def __rmul__(self, other): return Binop('*', other, self)
    def __div__(self, other): return Binop('/', self, other)
    def __rdiv__(self, other): return Binop('/', other, self)
    def __pow__(self, other): return Binop('**', self, other)
    def __rpow__(self, other): return Binop('**', other, self)

    # one out of place: syntactic sugar for 'eval'
    # this lets me say f.where(x = 2) instead of f.eval({'x':2})
    def where(self, **vars): return self.eval(vars)

# simplify an addition expression by dropping zeroes
def simplify_add(a, b):
    if a.identical(mkf(0)): return b
    elif b.identical(mkf(0)): return a
    else: return Binop('+', a, b)

# simplify a multiplication expression by dropping ones and converting
# 0 * anything to 0
def simplify_multiply(a, b):
    if a.identical(mkf(0)) or b.identical(mkf(0)): return mkf(0)
    elif a.identical(mkf(1)): return b
    elif b.identical(mkf(1)): return a
    else: return Binop('*', a, b)

def simplify_subtract(a, b):
    if b.identical(mkf(0)): return a
    else: return Binop('-', a, b)

# Binary operation class
class Binop(Formula):
    opmap = { '+': lambda a, b: a + b,
              '*': lambda a, b: a * b,
              '-': lambda a, b: a - b,
              '/': lambda a, b: a / b,
              '**': lambda a, b: a ** b }
    def __init__(self, op, value1, value2):
        self.op = op
        self.values = mkf(value1), mkf(value2)
    def __str__(self):
        return "(%s %s %s)" % (self.values[0], self.op, self.values[1])
    def eval(self, env):
        return self.opmap[self.op](self.values[0].eval(env),
                                   self.values[1].eval(env))
    # the partial derivative with respect to some variable 'var'
    derivmap = { '+': lambda a, b, var: a.derivative(var) + b.derivative(var),
                 '-': lambda a, b, var: a.derivative(var) - b.derivative(var),
                 '*': lambda a, b, var: (a * b.derivative(var) +
                                         b * a.derivative(var)) };
    def derivative(self, var):
        return self.derivmap[self.op](self.values[0], self.values[1], var)

    # very basic simplifications
    simplifymap = { '+': simplify_add,
                    '*': simplify_multiply,
                    '-': simplify_subtract};
    def simplified(self):
        if self.simplifymap.has_key(self.op):
            return self.simplifymap[self.op](self.values[0].simplified(),
                                             self.values[1].simplified())
        else:
            return self

    def identical(self, other):
        return (isinstance(other, Binop) and other.op == self.op and
                other.values[0].identical(self.values[0]) and
                other.values[1].identical(self.values[1]))


class Variable(Formula):
    def __init__(self, name): self._name = name
    def eval(self, environment): return environment[self._name]
    def __str__(self): return self._name
    def derivative(self, var):
        if self._name == var._name: return mkf(1)
        else: return mkf(0)
    def identical(self, other):
        return isinstance(other, Variable) and other._name == self._name
    def simplified(self): return self
class Constant(Formula):
    def __init__(self, value): self._value = value
    def eval(self, env): return self._value
    def __str__(self): return str(self._value)
    def derivative(self, var): return 0
    def identical(self, other):
        return isinstance(other, Constant) and other._value == self._value
    def simplified(self): return self

# make formula
def mkf(value):
    if type(value) in (type(1), type(1L), type(1.5), type(1j)):
        return Constant(value)
    elif type(value) is type("hi"):
        return Variable(value)
    elif isinstance(value, Formula):
        return value
    else:
        raise TypeError, ("Can't make formula from", value)

class Vars:
    def __getattr__(self, name): return Variable(name)
vars = Vars()

def test():
    assert mkf(2365).eval({}) == 2365
    one = mkf(1)
    assert str(one) == '1'
    assert one.eval({}) == 1
    assert isinstance(one + one, Formula)
    assert (one + one).eval({}) == 2
    assert str(one + one) == '(1 + 1)'
    x = vars.x
    assert isinstance(x, Variable)
    assert x.eval({'x': 37}) == 37
    assert (one + x).eval({'x': 108}) == 109
    assert str(one + x) == '(1 + x)'
    got_error = 0
    try:
        x.eval({})
    except KeyError:
        got_error = 1
    assert got_error
    assert (1 + one).eval({}) == 2
    assert (2 * mkf(3)).eval({}) == 6
    assert (mkf(2) * 3).eval({}) == 6
    assert (14 - one).eval({}) == 13
    assert (one - 14).eval({}) == -13
    assert int(one) == 1
    seven = (14 / mkf(2))
    assert isinstance(seven, Formula)
    assert seven.eval({}) == 7
    assert float(seven) == 7.0
    assert int(+one) == 1
    got_error = 0
    try:
        z = mkf(test)
    except TypeError:
        got_error = 1
    assert got_error
    two_to_the_x = (2 ** x)
    assert str(two_to_the_x) == '(2 ** x)'
    assert two_to_the_x.eval({'x': 20}) == 1048576
    assert two_to_the_x.where(x=20) == 1048576
    assert (x ** 2).eval({'x': 13}) == 169
    formula = (x + 1)/((x * x) - +two_to_the_x)
    assert str(formula) == '((x + 1) / ((x * x) - (2 ** x)))', str(formula)
    assert (x / 1).eval({'x': 36}) == 36
    assert long(one) == 1L
    assert complex(one) == 1+0j
    i = mkf(1j)
    assert complex(i) == 1j

    y = vars.y
    assert x.derivative(x).simplified().identical(mkf(1))
    assert x.derivative(y).simplified().identical(mkf(0))
    assert (x * y).derivative(x).simplified().identical(y)
    assert (x + y).derivative(x).simplified().identical(mkf(1))
    assert (x - y).derivative(x).simplified().identical(mkf(1))
    assert two_to_the_x.simplified().identical(two_to_the_x)

test()

-- 
<kragen@pobox.com>       Kragen Sitaker     <http://www.pobox.com/~kragen/>
Perilous to all of us are the devices of an art deeper than we possess
ourselves.
       -- Gandalf the White [J.R.R. Tolkien, "The Two Towers", Bk 3, Ch. XI]