more files

core.py

@@ -2,72 +2,81 @@
 import gof
 import numpy
 
-from copy import copy
+from copy import copy as pycopy
 
 # __all__ = ['set_mode', 'get_mode', 'NumpyR', 'NumpyOp']
 
 
-_mode = 'eval'
+_mode = ['eval']
 
 def set_mode(mode):
-    global _mode
-    _mode = mode
+    _mode.append(mode)
 
-def get_mode():
+def current_mode():
     return _mode
 
-def start_build():
+def build_mode():
     set_mode('build')
 
-def end_build():
+def eval_mode():
+    set_mode('eval')
+
+def pop_mode():
+    if len(_mode) == 1:
+        raise Exception("There's only one mode left on the stack.")
+    else:
+        _mode.pop()
+
+def end_eval():
     set_mode('eval')
 
 def build(f, *args, **kwargs):
-    start_build()
+    build_mode()
     r = f(*args, **kwargs)
-    end_build()
+    pop_mode()
     return r
 
 
+class Keyword:
 
-# class Proxy(object):
-
-#     __slots__ = ['_obj']
+    def __init__(self, name, nonzero=True):
+        self.name = name
+        self.nonzero = nonzero
 
-#     def __init__(self, obj = None):
-#         self._obj = obj
+    def __nonzero__(self):
+        return self.nonzero
 
-#     def __getattribute__(self, attr):
-#         if attr in ['__class__', '_obj']:
-#             return object.__getattribute__(self, attr)
-#         return getattr(self._obj, attr)
+    def __str__(self):
+        return "<%s>" % self.name
 
-#     def __setattr__(self, attr, value):
-#         if attr in ['_obj']:
-#             object.__setattr__(self, attr, value)
-#         else:
-#             setattr(self._obj, attr, value)
+    def __repr__(self):
+        return str(self)
 
-#     def __delattr__(self, attr):
-#         delattr(self._obj, attr)
+UNCOMPUTED = Keyword("UNCOMPUTED", False)
+UNDEFINED = Keyword("UNDEFINED", False)
 
-#     def __str__(self):
-#         return str(self._obj)
 
-#     def __iadd__(self, y):
-#         newobj = self._obj.__iadd__(y)
-#         if isinstance(newobj, Proxy):
-#             newobj = newobj._obj
-#         self._obj = newobj
-#         return self
+class Proxy(object):
 
+    __slots__ = ['_obj']
     
+    def __init__(self, obj = None):
+        self._obj = obj
 
-# class NumpyR(numpy.ndarray):
+    def __getattribute__(self, attr):
+        if attr in ['__class__', '_obj']:
+            return object.__getattribute__(self, attr)
+        else:
+            return getattr(object.__getattribute__(self, '_obj'), attr)
 
-#     __add__  = wrapper(numpy.ndarray.__add__)
-#     __radd__ = wrapper(numpy.ndarray.__radd__)
+    def __setattr__(self, attr, value):
+        if attr in ['_obj']:
+            object.__setattr__(self, attr, value)
+        else:
+            setattr(self._obj, attr, value)
 
+    def __delattr__(self, attr):
+        delattr(self._obj, attr)
 
 
 class IViewer(gof.ext.Viewer):
@@ -75,8 +84,11 @@ class IViewer(gof.ext.Viewer):
 
     def view_map(self):
         rval = {}
-        for i, o in self._v_map.items():
-            rval[self.inputs[i]] = self.outputs[o]
+        for output, inputs in self._v_map.items():
+            if isinstance(inputs, (list, tuple)):
+                rval[self.outputs[output]] = [self.inputs[i] for i in inputs]
+            else:
+                rval[self.outputs[output]] = self.inputs[inputs]
         return rval
 
 
@@ -85,135 +97,473 @@ class IDestroyer(gof.ext.Destroyer):
 
     def destroy_map(self):
         rval = {}
-        for i, o in self._d_map.items():
-            rval[self.inputs[i]] = self.outputs[o]
+        for output, inputs in self._d_map.items():
+            if isinstance(inputs, (list, tuple)):
+                rval[self.outputs[output]] = [self.inputs[i] for i in inputs]
+            else:
+                rval[self.outputs[output]] = self.inputs[inputs]
         return rval
 
 
+class PythonR(gof.HolderResult):
+
+    def __init__(self, a = None):
+        if a is None:
+            self.storage = UNCOMPUTED
+        else:
+            self.storage = a
+
+    def set_value(self, value):
+        self.storage = value
+
+    def __str__(self):
+        return str(self.storage)
+
+    def __repr__(self):
+        return repr(self.storage)
+
+
+class NumpyR(gof.HolderResult):
+
+    def __init__(self, a = None):
+        self.set_value(a)
+
+    def set_value(self, value):
+        if value is None or value is UNCOMPUTED:
+            self.storage = UNCOMPUTED
+        elif isinstance(value, numpy.ndarray):
+            self.storage = value
+        else:
+            self.storage = numpy.array(value)
+        
+    def __add__(self, y):
+        return add(self, y)
+
+    def __radd__(self, x):
+        return add(x, self)
+
+    def __iadd__(self, y):
+        return iadd(self, y)
+    
+    def __sub__(self, y):
+        return sub(self, y)
+
+    def __rsub__(self, x):
+        return sub(x, self)
+
+    def __isub__(self, y):
+        return isub(self, y)
+    
+    def __mul__(self, y):
+        return dot(self, y)
+
+    def __rmul__(self, x):
+        return dot(x, self)
+
+    def __imul__(self, y):
+        return imul(self, y)
+ 
+    def __div__(self, y):
+        return div(self, y)
+
+    def __rdiv__(self, x):
+        return div(x, self)
+
+    def __idiv__(self, y):
+        return idiv(self, y)
+    
+    def __mod__(self, y):
+        return mod(self, y)
+
+    def __rmod__(self, x):
+        return mod(x, self)
+
+    def __pow__(self, y):
+        return pow(self, y)
+
+    def __rpow__(self, x):
+        return pow(x, self)
+
+    def __ipow__(self, y):
+        return ipow(self, y)
+
+    def __neg__(self):
+        return neg(self)
+
+    T = property(lambda self: transpose(self))
+    Tc = property(lambda self: transpose_copy(self))
+
+    def __copy__(self):
+        return array_copy(self)
+    
+    def __str__(self):
+        return str(self.storage)
+
+    def __repr__(self):
+        return repr(self.storage)
+    
+#[iadd(iadd(iadd(iadd(<UNCOMPUTED>, itwice(<UNCOMPUTED>)), <UNCOMPUTED>), 1.0), dot(<UNCOMPUTED>, <UNCOMPUTED>))]
+#[iadd(iadd(iadd(iadd(<UNCOMPUTED>, itwice(<UNCOMPUTED>)), <UNCOMPUTED>), 1.0), dot(<UNCOMPUTED>, <UNCOMPUTED>))]
+
+
+def wrap(x):
+#     try:
+#         return to_numpyr(x)
+#     except TypeError:
+#         if isinstance(x, PythonR):
+#             return x
+#         else:
+#             return PythonR(x)
+
+
+# def to_numpyr(x):
+    
+    if isinstance(x, NumpyR):
+        return x
+    elif isinstance(x, PythonR):
+        return x
+    elif isinstance(x, NumpyOp):
+        return x.out
+    elif isinstance(x, Proxy):
+        return wrap(x._obj)
+    elif isinstance(x, numpy.ndarray):
+        return NumpyR(x)
+    else:
+        return PythonR(x)
+#     else:
+#         raise TypeError("%s cannot be converted to or encapsulated in a NumpyR instance." % x)
+
+
+class NumpyOp(gof.Op, gof.ext.BuildableFromInputs):
+
+    nout = 1
+    
+    def __init__(self, *args):
+        inputs = [wrap(arg) for arg in args]
+        outputs = [NumpyR() for i in xrange(self.nout)]
+        gof.Op.__init__(self, inputs, outputs)
+
+    @classmethod
+    def from_inputs(cls, *inputs):
+        return cls(*inputs)
+
+    def gen_outputs(self):
+        return [NumpyR() for i in xrange(self.nout)]
+
+
+
 class wrapper:
 
     __slots__ = ['f', 'opclass']
     
-    def __init__(self, f, vmap = None, dmap = None):
+    def __init__(self, name, f, grad, vmap = None, dmap = None, optype = NumpyOp):
         self.f = f
 
         if not callable(f):
             raise TypeError("Can only wrap a callable.")
 
-        bases = [NumpyOp]
-        if vmap:
-            bases.append(IViewer)
-        if dmap:
-            bases.append(IDestroyer)
-        
-        if hasattr(f, '__module__'):
-            mod = f.__module__
-        elif hasattr(f, '__objclass__'):
-            mod = f.__objclass__.__name__
-        else:
-            mod = ""
+        bases = [optype]
+        if vmap: bases.append(IViewer)
+        if dmap: bases.append(IDestroyer)
         
-        Wrapper = type("Op:" + mod + ":" + f.__name__, tuple(bases), {})
+        Wrapper = type(name, tuple(bases), {})
         if vmap: Wrapper._v_map = vmap
         if dmap: Wrapper._d_map = dmap
         
         def thunk(self):
             def ret():
-                self.outputs[0].storage = f(*[input.storage for input in self.inputs])
+                self.outputs[0].set_value(f(*[input.storage for input in self.inputs]))
             return ret
+        Wrapper.thunk = thunk
         
-        Wrapper.thunk = thunk #.impl = staticmethod(impl)
+        if grad is UNDEFINED:
+            grad = lambda *_: UNDEFINED
+        Wrapper.grad = staticmethod(grad)
 
         self.opclass = Wrapper
 
     
     def __call__(self, *args):
         op = self.opclass(*args)
-        if _mode == 'eval':
+        if current_mode() == 'eval':
             op.thunk()()
-#        outputs = [Proxy(o) for o in op.outputs]
-        outputs = copy(op.outputs)
+        outputs = pycopy(op.outputs)
+#        outputs = [Proxy(output) for output in op.outputs]
         if op.nout == 1:
             return outputs[0]
         else:
             return outputs
 
 
-#         if '_mode' in kwargs:
-#             mode = kwargs['_mode']
-#             del kwargs['_mode']
-#         else:
-#             mode = _mode
+# def wrap_producer(f):
+#     def ret(*args, **kwargs):
+#         result = f(*args, **kwargs)
+#         if not isinstance(result, numpy.ndarray):
+#             result = numpy.array(result)
+#         return NumpyR(result)
+#     return ret
 
-#         if mode == 'eval':
-#             return f(*[to_ndarray(arg) for arg in args])
-#         elif mode == 'build':
-#             o = self.op(*args)
 
-#         elif mode == 'build_eval':
-#             return f(*args, **kwargs)
+def wrap_producer(f):
+    wrapped_f = wrapper(f.__name__, f, UNDEFINED)
+    def ret(dim, dtype = 'float', order = 'C'):
+        return wrapped_f(dim, dtype, order)
+    return ret
 
 
 
-class NumpyR(gof.HolderResult):
 
-    def __init__(self, a = None):
-        if a is None:
-            self.storage = None
-        elif isinstance(a, numpy.ndarray):
-            self.storage = a
+ndarray = wrap_producer(numpy.ndarray)
+array = wrap_producer(numpy.array)
+zeros = wrap_producer(numpy.zeros)
+ones = wrap_producer(numpy.ones)
+
+
+inplace = gof.ext.Destroyer
+view = gof.ext.Viewer
+
+
+class omega_op_metaclass(type):
+    
+    def __init__(cls, name, bases, dct):
+        type.__init__(cls, name, bases, dct)
+        cls.__clsinit__(name, bases, dct)
+
+
+
+class omega_op(gof.Op, gof.ext.BuildableFromInputs):
+
+    __metaclass__ = omega_op_metaclass
+    
+    nout = 1
+
+    @classmethod
+    def __clsinit__(cls, name, bases, dct):
+        # make grad and impl static methods
+        grad = cls.grad
+        if hasattr(grad, 'im_func'):
+            grad = grad.im_func
+        impl = cls.impl
+        if hasattr(cls.impl, 'im_func'):
+            impl = impl.im_func
+        cls.grad = staticmethod(grad)
+        cls.impl = staticmethod(impl)
+    
+    def __new__(cls, *inputs):
+        op = gof.Op.__new__(cls)
+        op.__init__(*[wrap(input) for input in inputs])
+        if current_mode() == 'eval':
+            op.thunk()()
+        if op.nout == 1:
+            return op.out
         else:
-            raise TypeError("NumpyR constructor expects a ndarray instance.")
+            return op.outputs
     
-    def __add__(self, y):
-        return add(self, y)
+    def __init__(self, *inputs):
+        for input in inputs:
+            assert isinstance(input, gof.HolderResult)
+        gof.Op.__init__(self, inputs, self.gen_outputs())
 
-    def __mul__(self, y):
-        return dot(self, y)
+    @classmethod
+    def from_inputs(cls, *inputs):
+        build_mode()
+        r = cls(*inputs)
+        pop_mode()
+        return r.owner
 
-    def __iadd__(self, y):
-        return iadd(self, y)
+    def gen_outputs(self):
+        return [NumpyR() for i in xrange(self.nout)]
 
-    def __str__(self):
-        return str(self.storage)
+    def thunk(self):
+        def ret():
+            results = self.impl(*[input.storage for input in self.inputs])
+            if self.nout == 1:
+                self.out.set_value(results)
+            else:
+                assert self.nout == len(results)
+                for result, output in zip(results, self.outputs):
+                    output.set_value(result)
+        return ret
     
-    def __repr__(self):
-        return repr(self.storage)
+    def update_gradient(self, grad_d):
+        inputgs = self.grad(*(self.inputs + [grad_d[output] for output in self.outputs]))
+        if not isinstance(inputgs, (list, tuple)):
+            inputgs = [inputgs] * len(self.inputs)
+        for input, inputg in zip(self.inputs, inputgs):
+            grad_d.add(input, inputg)
 
+    def grad(*args):
+        return UNDEFINED
 
+    def impl(*args):
+        raise NotImplementedError("This op has no implementation.")
 
-def to_numpyr(x):
-    if isinstance(x, NumpyR):
+    
+
+## Addition ##
+
+class proto_add(omega_op):
+    def grad(x, y, gz):
+        return gz
+
+class add(proto_add):
+    impl = numpy.ndarray.__add__
+
+class iadd(proto_add, inplace):
+    impl = numpy.ndarray.__iadd__
+
+
+class proto_twice(omega_op):
+    def grad(x, gz):
+        return scal(gz, 2.0)
+
+class twice(proto_twice):
+    def impl(x):
+        return x + x
+
+class itwice(proto_twice, inplace):
+    def impl(x):
+        x += x
         return x
-    elif isinstance(x, NumpyOp):
-        return x.out
-#    elif isinstance(x, Proxy):
-#        return to_numpyr(x._obj)
-    elif isinstance(x, numpy.ndarray):
-        return NumpyR(x)
-    else:
-        raise TypeError("%s cannot be converted to or encapsulated in a NumpyR instance." % x)
 
 
+## Subtraction ##
 
-class NumpyOp(gof.Op):
+class proto_sub(omega_op):
+    def grad(x, y, gz):
+        return gz, -gz
 
-    nout = 1
+class sub(proto_sub):
+    impl = numpy.ndarray.__sub__
 
-    def __init__(self, *args):
-        inputs = [to_numpyr(arg) for arg in args]
-        outputs = [NumpyR() for i in xrange(self.nout)]
-        gof.Op.__init__(self, inputs, outputs)
+class isub(proto_sub, inplace):
+    impl = numpy.ndarray.__isub__
+
+
+
+## Element-wise multiplication ##
+
+def mul_grad(x, y, gz):
+    return mul(y, gz), mul(x, gz)
+
+mul = wrapper("mul",
+              numpy.ndarray.__mul__,
+              mul_grad)
+
+imul = wrapper("imul",
+               numpy.ndarray.__imul__,
+               mul_grad,
+               dmap = {0: 0})
+
+def sqr_grad(x, gz):
+    return scal(mul(x, gz), 2.0)
+
+sqr = wrapper("sqr",
+              lambda x: numpy.multiply(x, x),
+              sqr_grad)
+
+isqr = wrapper("isqr",
+               lambda x: x.__imul__(x),
+               sqr_grad,
+               dmap = {0: 0})
+
+def sqrt_grad(x, gz):
+    return scal(div(gz, sqrt(x)), 0.5)
+
+sqrt = wrapper("sqrt",
+               lambda x: numpy.sqrt(x),
+               sqrt_grad)
+
+isqrt = wrapper("isqrt",
+                lambda x: x.__ipow__(0.5),
+                sqrt_grad,
+                dmap = {0: 0})
+
+
+## Element-wise division ##
+
+def div_grad(x, y, gz):
+    return div(gz, y), -div(mul(x, gz), sqr(y))
+
+div = wrapper("div",
+              numpy.ndarray.__div__,
+              div_grad)
+
+idiv = wrapper("idiv",
+               numpy.ndarray.__idiv__,
+               div_grad,
+               dmap = {0: 0})
+
+
+## Scaling ##
+
+def scal_grad(x, a, gz):
+    return scal(a, gz), sum(mul(x, gz))
+
+scal = wrapper("scal",
+               numpy.ndarray.__mul__,
+               scal_grad)
+
+iscal = wrapper("iscal",
+                numpy.ndarray.__imul__,
+                scal_grad,
+                dmap = {0: 0})
+
+neg = wrapper("neg",
+              numpy.ndarray.__neg__,
+              lambda x, gz: -gz)
+
+ineg = wrapper("ineg",
+               lambda x: x.__imul__(-1),
+               lambda x, gz: -gz,
+               dmap = {0: 0})
+
+
+## Dot product ##
+
+dot = wrapper("dot",
+              numpy.dot,
+              lambda x, y, gz: (dot(gz, transpose(y)),
+                                dot(transpose(x), gz)))
+
+
+## Transposition ##
+
+class transpose(omega_op, view):
+    impl = numpy.transpose
+    def grad(x, gz):
+        return transpose_copy(gz)
+
+# transpose = wrapper("transpose",
+#                     numpy.transpose,
+#                     lambda x, z, gz: transpose_copy(gz),
+#                     vmap = {0: 0})
+
+def transpose_copy(x):
+    return array_copy(transpose(x))
+
+
+## Copy ##
+
+class array_copy(omega_op):
+    impl = numpy.array
+    grad = lambda x, gz: gz
+
+
+## Others ##
+
+class minmax(omega_op):
+    nout = 2
+    def impl(x):
+        return x.min, x.max
+
+    
+# array_copy = wrapper("copy",
+#                      numpy.array,
+#                      lambda x, gz: gz)
 
-#     def __validate__(self):
-#         for input in self.inputs:
-#             assert isinstance(input, NumpyR)
-#         for output in self.outputs:
-#             assert isinstance(output, NumpyR)
 
+# array slicing
 
-add = wrapper(numpy.ndarray.__add__)
-iadd = wrapper(numpy.ndarray.__iadd__, None, {0: 0})
 
-dot = wrapper(numpy.dot)
 

grad.py

+
+import gof
+import core
+
+class _GradD(dict):
+    """A dictionary-like class, into which derivative expressions may be added"""
+
+    def add(self, r, dr):
+        """Add dv to the sum of gradients associated with v"""
+
+        if r is core.UNDEFINED:
+            self[r] = core.UNDEFINED
+        elif r in self:
+            self[r] = self[r] + dr
+        else:
+            self[r] = dr
+
+
+def expand_grad(i, o, cost_derivs):
+    grad_d = _GradD(cost_derivs)
+    core.build_mode()
+    for op in gof.graph.io_toposort(i, o).__reversed__():
+        op.update_gradient(grad_d)
+#        inputgs = op.grad(*(op.inputs + [grad_d[output] for output in op.outputs]))
+#        if not isinstance(inputgs, (list, tuple)):
+#            inputgs = [inputgs] * len(op.inputs)
+#        for input, inputg in zip(op.inputs, inputgs):
+#            grad_d.add(input, inputg)
+    core.pop_mode()
+    return grad_d
+
+
+def grad(cost, wrt, cost_grad = 1.0):
+#    cost, wrt = core.wrap(cost), core.wrap(wrt)
+    cost_derivs = expand_grad([wrt], [cost], {cost: core.wrap(cost_grad)})
+#     print wrt
+#     for k, v in cost_derivs.items():
+#         print k, v
+    ret = cost_derivs.get(wrt, None)
+    if ret is core.UNDEFINED:
+        raise Exception("The gradient wrt %s is undefined." % wrt)
+    return ret
+
+

opt.py

+
+from core import *
+import gof
+
+
+def pattern_opt(in_pattern, out_pattern):
+    def parse(x):
+        if isinstance(x, (list, tuple)):
+            return [parse(y) for y in x]
+        elif isinstance(x, wrapper):
+            return x.opclass
+        elif isinstance(x, str) or (hasattr(x, '__bases__') and issubclass(x, gof.op.Op)):
+            return x
+        else:
+            raise TypeError("Bad input type for pattern_opt.")
+    return gof.opt.PatternOptimizer(parse(in_pattern), parse(out_pattern))
+
+def op_sub(op1, op2):
+    if isinstance(op1, wrapper):
+        op1 = op1.opclass
+    if isinstance(op2, wrapper):
+        op2 = op2.opclass
+    return gof.opt.OpSubOptimizer(op1, op2)
+
+
+#def make_patterns(patterns):
+#    return [name, pattern_opt(inp, outp) for name, inp, outp in patterns]
+
+def export_opts(opts):
+    for name, opt in opts:
+        if name:
+            globals()[name] = opt
+
+
+
+# double_transpose_eliminator = pattern_opt((transpose, (transpose, 'x')), 'x')
+# patterns = make_patterns(patterns)
+# export_patterns(patterns)
+
+# List of optimizations to perform. They are listed in the order they are applied.
+opts = [
+
+    ['double_transpose_eliminator', pattern_opt((transpose, (transpose, 'x')),
+                                                'x')],
+
+    ['addxx_to_twice',              pattern_opt((add, 'x', 'x'),
+                                                (twice, 'x'))],
+
+    ['twice_to_itwice',             op_sub(twice, itwice)],
+
+    ['mulxx_to_twice',              pattern_opt((mul, 'x', 'x'),
+                                                (sqr, 'x'))],
+
+    ['sqr_to_isqr',                 op_sub(sqr, isqr)],
+    
+    ['add_to_iadd',                 op_sub(add, iadd)],
+
+    ['add_to_iadd_reverse',         pattern_opt((add, 'x', 'y'),
+                                                (iadd, 'y', 'x'))],
+
+    ['remove_copies',               gof.opt.OpRemover(array_copy)],
+    
+    [None,                          gof.lib.DummyRemover] # has to be at the end
+    
+    ]
+
+export_opts(opts) # publish the optimizations performed under individual names
+
+
+optimizer = gof.opt.MergeOptMerge(gof.opt.SeqOptimizer([opt for name, opt in opts]))

rand.py

+
+import core
+import gof
+from numpy import random as r
+
+
+# def rwrap(f):
+#     wrapped = 
+#     def ret(self, *args):
+        
+    
+
+class RandomState(gof.Op, gof.ext.IONames):
+
+    input_names = ['seed']
+
+    def __init__(self, seed):
+        inputs = [wrap(seed)]
+        outputs = [PythonR()]
+        gof.Op.__init__(self, inputs, outputs)
+
+    def thunk(self):
+        def f():
+            self.out.storage = r.RandomState(self.seed.storage)
+        return f
+
+    
+
+
+
+
+class Random(object):
+
+    def __init__(seed):
+        self.state = core.wrap(seed)
+
+    
+