diff --git a/TODO b/TODO
index a410aec4b71d4aa2d84b000fcb1e19cc5d0d8efb..3c732f00cbbed6325043036cf52cbf08f69ecb2d 100644
--- a/TODO
+++ b/TODO
@@ -1,3 +1,5 @@
 Build out traits
 
 add traits types to expressions
+
+simplification
diff --git a/src/__init__.py b/src/__init__.py
index 65da7c62222320a68e76469973c0e29a810ef2d8..9a25c874c66d82eab152aa5b75b722f4dc12e96c 100644
--- a/src/__init__.py
+++ b/src/__init__.py
@@ -24,17 +24,49 @@ get_dependencies = pymbolic.mapper.dependency.get_dependencies
 substitute = pymbolic.mapper.substitutor.substitute
 differentiate = pymbolic.mapper.differentiator.differentiate
 
+
+
+
 def simplify(x):
+    # FIXME: Not yet implemented
     return x
 
-def grad(expr, variables):
+def grad(expression, variables):
     return [simplify(differentiate(expression, var)) for var in variables]
 
 def jacobian(expression_list, variables):
     return [grad(expr, variables) for expr in expression_list]
     
 def laplace(expression, variables):
-    return sum([differentiate(differentiate(expression,var),var) for var in variables])
+    return sum(*[differentiate(differentiate(expression,var),var) for var in variables])
+
+
+
+
+class VectorFunction:
+    def __init__(self, function_list, variables=[]):
+        self.FunctionList = [pymbolic.compile(expr, variables=variables) 
+                             for expr in function_list]
+
+    def __call__(self, x):
+        import pylinear.array as num
+        return num.array([ func(x) for func in self.FunctionList ])
+
+
+
+
+class MatrixFunction:
+    def __init__(self, function_list, variables=[]):
+        self. FunctionList = [[pymbolic.compile(expr, variables=variables)
+                               for expr in outer]
+                              for outer in function_list]
+
+    def __call__(self, x):
+        import pylinear.array as num
+        return num.array([[func(x) for func in flist ] for flist in self.FunctionList])
+
+
+
 
 if __name__ == "__main__":
     import math
diff --git a/src/conversion.py b/src/conversion.py
deleted file mode 100644
index c84ded278778d9a68ff68fabc98a89344e66f349..0000000000000000000000000000000000000000
--- a/src/conversion.py
+++ /dev/null
@@ -1,8 +0,0 @@
-def to_polynomial(x):
-    raise NotImplementedError
-
-def to_expression(x):
-    raise NotImplementedError
-
-def to_fraction(x, approx_bits=64-4):
-    raise NotImplementedError
diff --git a/src/expression_operators.py b/src/expression_operators.py
deleted file mode 100644
index 637969daa2480cd211f0ebee07c1a6efed2a7dfd..0000000000000000000000000000000000000000
--- a/src/expression_operators.py
+++ /dev/null
@@ -1,16 +0,0 @@
-PLUS = intern("+")
-MINUS = intern("-")
-TIMES = intern("*")
-DIVIDE = intern("/")
-POWER = intern("**")
-
-NEG = intern("negate")
-SIN = intern("sin")
-COS = intern("cos")
-TAN = intern("tan")
-LOG = intern("log")
-EXP = intern("exp")
-
-VARIABLE = intern("variable")
-
-IF = intern("variable")
diff --git a/src/functions.py b/src/functions.py
new file mode 100644
index 0000000000000000000000000000000000000000..e7ee6719f7abe89d1a2e9e7e4be2aa2cfd6f2f03
--- /dev/null
+++ b/src/functions.py
@@ -0,0 +1,16 @@
+import pymbolic.primitives as primitives
+
+
+
+
+def sin(x):
+    return primitives.Call(primitives.ElementLookup(primitives.Variable("math"), "sin"), (x,))
+def cos(x):
+    return primitives.Call(primitives.ElementLookup(primitives.Variable("math"), "cos"), (x,))
+def tan(x):
+    return primitives.Call(primitives.ElementLookup(primitives.Variable("math"), "tan"), (x,))
+def log(x):
+    return primitives.Call(primitives.ElementLookup(primitives.Variable("math"), "log"), (x,))
+def exp(x):
+    return primitives.Call(primitives.ElementLookup(primitives.Variable("math"), "exp"), (x,))
+
diff --git a/src/mapper/__init__.py b/src/mapper/__init__.py
index 9c78539008210a429bd21e57504976393820c9ba..7f3fb4821ce2fe5c755882df5cd84aa95bb66f05 100644
--- a/src/mapper/__init__.py
+++ b/src/mapper/__init__.py
@@ -21,12 +21,12 @@ class CombineMapper:
         return self.combine(child.invoke_mapper(self)
                             for child in expr.children)
 
+    map_product = map_sum
+
     def map_rational(self, expr):
         return self.combine((expr.numerator.invoke_mapper(self),
                              expr.denominator.invoke_mapper(self)))
 
-    map_product = map_sum
-
     def map_power(self, expr):
         return self.combine((expr.base.invoke_mapper(self),
                              expr.exponent.invoke_mapper(self)))
@@ -64,8 +64,8 @@ class IdentityMapper:
         return expr.__class__(expr.child.invoke_mapper(self))
 
     def map_sum(self, expr):
-        return expr.__class__(tuple(child.invoke_mapper(self)
-                                    for child in expr.children))
+        return expr.__class__(*[child.invoke_mapper(self)
+                                for child in expr.children])
     
     map_product = map_sum
     
diff --git a/src/mapper/differentiator.py b/src/mapper/differentiator.py
index 9e2d57920533534563b3be0b2a63e2d5f1c65c19..e6ce01537ca58e0db250d7e7033af508d915fcf6 100644
--- a/src/mapper/differentiator.py
+++ b/src/mapper/differentiator.py
@@ -1,28 +1,33 @@
 import math
+import cmath
 
 import pymbolic
 import pymbolic.primitives as primitives
-
+import pymbolic.mapper.evaluator
 
 
 
 def map_math_functions_by_name(i, func, pars):
-    if not isinstance(func, primitives.Variable):
+    try:
+        f = pymbolic.evaluate(func, {"math": math, "cmath": cmath})
+    except pymbolic.mapper.evaluator.UnknownVariableError:
         raise RuntimeError, "No derivative of non-constant function "+str(func)
-    name = func.name
-
-    if name == "sin" and len(pars) == 1:
-        return primitives.Constant("cos")(*pars)
-    elif name == "cos" and len(pars) == 1:
-        return -primitives.Constant("sin")(*pars)
-    elif name == "tan" and len(pars) == 1:
-        return primitives.Constant("tan")(*pars)**2+1
-    elif name == "log" and len(pars) == 1:
+
+    def make_f(name):
+        return primitives.ElementLookup(primitives.Variable("math"), name)
+
+    if f is math.sin and len(pars) == 1:
+        return make_f("cos")(*pars)
+    elif f is math.cos and len(pars) == 1:
+        return -make_f("sin")(*pars)
+    elif f is math.tan and len(pars) == 1:
+        return make_f("tan")(*pars)**2+1
+    elif f is math.log and len(pars) == 1:
         return primitives.Constant(1)/pars[0]
-    elif name == "exp" and len(pars) == 1:
-        return primitives.Constant("exp")(*pars)
+    elif f is math.exp and len(pars) == 1:
+        return make_f("exp")(*pars)
     else:
-        return primitives.Constant(name+"'")(*pars)
+        raise RuntimeError, "unrecognized function, cannot differentiate"
 
 
 
@@ -45,7 +50,7 @@ class DifferentiationMapper:
             return primitives.Constant(0)
 
     def map_call(self, expr):
-        return pymbolic.sum(tuple(
+        return pymbolic.sum(*(
             self.FunctionMap(i, expr.function, expr.parameters)
             * par.invoke_mapper(self)
             for i, par in enumerate(expr.parameters)
diff --git a/src/mapper/evaluator.py b/src/mapper/evaluator.py
index c180d21bf212d69ca0ade29053b7b008bd74a23c..19604f8fd463460bab917b5df4f903654999730d 100644
--- a/src/mapper/evaluator.py
+++ b/src/mapper/evaluator.py
@@ -1,3 +1,9 @@
+class UnknownVariableError(Exception):
+    pass
+
+
+
+
 class EvaluationMapper:
     def __init__(self, context={}):
         self.Context = context
@@ -6,7 +12,10 @@ class EvaluationMapper:
         return expr.value
 
     def map_variable(self, expr):
-        return self.Context[expr.name]
+        try:
+            return self.Context[expr.name]
+        except KeyError:
+            raise UnknownVariableError, expr.name
 
     def map_call(self, expr):
         return expr.function.invoke_mapper(self)(
diff --git a/src/parser.py b/src/parser.py
index 9e4a0d06d49ae9bd517e492f5fc3266ce62da450..c2b7f15c8c2cab9802cdf1b72bf2ecd168e65492 100644
--- a/src/parser.py
+++ b/src/parser.py
@@ -74,12 +74,11 @@ def parse(expr_str):
             return primitives.Negation(parse_expression(pstate, _PREC_UNARY_MINUS))
         if pstate.is_next(_openpar):
             pstate.advance()
-            result = parse_expression(pstate)
+            left_exp = parse_expression(pstate)
             pstate.expect(_closepar)
             pstate.advance()
-            return result
-
-        left_exp = parse_terminal(pstate)
+        else:
+            left_exp = parse_terminal(pstate)
 
         did_something = True
         while did_something:
diff --git a/src/primitives.py b/src/primitives.py
index 6c71c4ace1c1a97786492410e2e60cbe3e68a165..251318b845989a884c131b3c414d91724f32467d 100644
--- a/src/primitives.py
+++ b/src/primitives.py
@@ -117,8 +117,7 @@ class Expression(object):
         return self.invoke_mapper(pymbolic.mapper.stringifier.StringifyMapper())
 
     def __repr__(self):
-        return "%s(%s)" % (self.__class__.__name__,
-                           ", ".join(repr(i) for i in self.__getinitargs__()))
+        return "%s%s" % (self.__class__.__name__, repr(self.__getinitargs__()))
 
 class Constant(Expression):
     def __init__(self, value):
@@ -356,7 +355,6 @@ class Negation(Expression):
 
 class Sum(Expression):
     def __init__(self, *children):
-        assert isinstance(children, tuple)
         self._Children = children
 
     def __getinitargs__(self):
@@ -381,7 +379,7 @@ class Sum(Expression):
             other = Constant(other)
         elif isinstance(other, Sum):
             return Sum(*(other._Children + self._Children))
-        return Sum(*(other, + self._Children))
+        return Sum(*((other,) + self._Children))
 
     def __sub__(self, other):
         if not isinstance(other, Expression):
@@ -402,7 +400,6 @@ class Sum(Expression):
 
 class Product(Expression):
     def __init__(self, *children):
-        assert isinstance(children, tuple)
         self._Children = children
 
     def __getinitargs__(self):
@@ -420,14 +417,14 @@ class Product(Expression):
             other = Constant(other)
         elif isinstance(other, Product):
             return Product(*(self._Children + other._Children))
-        return Product(self._Children + (other,))
+        return Product(*(self._Children + (other,)))
 
     def __rmul__(self, other):
         if not isinstance(other, Expression):
             other = Constant(other)
         elif isinstance(other, Product):
             return Product(*(other._Children + self._Children))
-        return Product(*(other, + self._Children))
+        return Product(*((other,) + self._Children))
 
     def __nonzero__(self):
         for i in self._Children:
@@ -438,26 +435,29 @@ class Product(Expression):
     def invoke_mapper(self, mapper):
         return mapper.map_product(self)
 
-class QuotientExpression(Expression):
+class Quotient(Expression):
     def __init__(self, numerator, denominator=1):
-        self.Numerator = numerator
-        self.Denominator = denominator
+        self._Numerator = numerator
+        self._Denominator = denominator
+
+    def __getinitargs__(self):
+        return self._Numerator, self._Denominator
 
     def _num(self):
-        return self.Numerator
+        return self._Numerator
     numerator=property(_num)
 
     def _den(self):
-        return self.Denominator
+        return self._Denominator
     denominator=property(_den)
 
     def __eq__(self, other):
         return isinstance(other, Subscript) \
-               and (self.Numerator == other.Numerator) \
-               and (self.Denominator == other.Denominator)
+               and (self._Numerator == other._Numerator) \
+               and (self._Denominator == other._Denominator)
 
     def __nonzero__(self):
-        return bool(self.Numerator)
+        return bool(self._Numerator)
 
     def invoke_mapper(self, mapper):
         return mapper.map_rational(self)
@@ -485,6 +485,9 @@ class Power(Expression):
         self._Base = base
         self._Exponent = exponent
 
+    def __getinitargs__(self):
+        return self._Base, self._Exponent
+
     def _base(self):
         return self._Base
     base = property(_base)
@@ -540,6 +543,14 @@ class List(Expression):
 
 
 # intelligent makers ---------------------------------------------------------
+def subscript(expression, index):
+    if not isinstance(index, Expression):
+        index = Constant(index)
+    return Subscript(expression, index)
+
+
+
+
 def sum(*components):
     components = tuple(c for c in components if c)
     if len(components) == 0:
@@ -553,10 +564,9 @@ def sum(*components):
 
 
 def linear_combination(coefficients, expressions):
-    print coefficients, expressions
-    return sum(*(coefficient * expression
+    return sum(*[coefficient * expression
                  for coefficient, expression in zip(coefficients, expressions)
-                 if coefficient and expression))
+                 if coefficient and expression])
 
 
 
@@ -578,14 +588,6 @@ def product(*components):
 
 
 
-def subscript(expression, index):
-    if not isinstance(index, Expression):
-        index = Constant(index)
-    return Subscript(expression, index)
-
-
-
-
 def polynomial_from_expression(expression):
     pass
 
@@ -602,5 +604,5 @@ def make_quotient(numerator, denominator):
     except traits.NoTraitsError:
         pass
 
-    return QuotientExpression(numerator, denominator)
+    return Quotient(numerator, denominator)