bringing back sparse

_test_sparse.py

 from sparse import *
 
 import unittest
+import compile
 
-class _testCase_transpose(unittest.TestCase):
+class T_transpose(unittest.TestCase):
     def setUp(self):
-        core.build_eval_mode()
         numpy.random.seed(44)
-    def tearDown(self):
-        core.pop_mode()
-    def test_transpose(self):
-        a = SparseR(sparse.csr_matrix(sparse.speye(5,3)))
+    def test_transpose_csc(self):
+        sp = sparse.csc_matrix(sparse.speye(5,3))
+        a = assparse(sp)
+        self.failUnless(a.data is sp)
         self.failUnless(a.data.shape == (5,3))
+        self.failUnless(a.dtype == 'float64')
+        self.failUnless(a.format == 'csc', a.format)
         ta = transpose(a)
-        self.failUnless(ta.data.shape == (3,5))
+        self.failUnless(ta.dtype == 'float64', ta.dtype)
+        self.failUnless(ta.format == 'csr', ta.format)
+
+        vta = compile.eval_outputs([ta])
+        self.failUnless(vta.shape == (3,5))
+    def test_transpose_csr(self):
+        a = assparse(sparse.csr_matrix(sparse.speye(5,3)))
+        self.failUnless(a.data.shape == (5,3))
+        self.failUnless(a.dtype == 'float64')
+        self.failUnless(a.format == 'csr')
+        ta = transpose(a)
+        self.failUnless(ta.dtype == 'float64', ta.dtype)
+        self.failUnless(ta.format == 'csc', ta.format)
+
+        vta = compile.eval_outputs([ta])
+        self.failUnless(vta.shape == (3,5))
+
+class T_Add(unittest.TestCase):
+    def test0(self):
+        sp_a = sparse.csc_matrix(sparse.speye(5,3))
+        a = assparse(sp_a)
+
+        sp_b = sparse.csc_matrix(sparse.speye(5,3))
+        b = assparse(sp_b)
+
+        self.failUnless(a.data is sp_a)
+        apb = add_s_s(a, b)
+
+        self.failUnless(apb.dtype == a.dtype, apb.dtype)
+        self.failUnless(apb.format == a.format, apb.format)
+
+        val = compile.eval_outputs([apb])
+        self.failUnless(val.shape == (5,3))
+        self.failUnless(numpy.all(val.todense() == (sp_a + sp_b).todense()))
+
+class T_conversion(unittest.TestCase):
+    def setUp(self):
+        numpy.random.seed(44)
+
+    def test0(self):
+        a = tensor.astensor(numpy.random.rand(5))
+        s = sparse_from_dense(a,'csc')
+        val = compile.eval_outputs([s])
+        self.failUnless(str(val.dtype)=='float64')
+        self.failUnless(val.format == 'csc')
+
+    def test1(self):
+        a = tensor.astensor(numpy.random.rand(5))
+        s = sparse_from_dense(a,'csr')
+        val = compile.eval_outputs([s])
+        self.failUnless(str(val.dtype)=='float64')
+        self.failUnless(val.format == 'csr')
+
+    def test2(self):
+        csr = sparse.csr_matrix((2,5))
+        d = dense_from_sparse(csr)
+        csr[0,0] = 1.0
+        val = compile.eval_outputs([d])
+        self.failUnless(str(val.dtype)=='float64')
+        self.failUnless(numpy.all(val[0] == [1,0,0,0,0]))
+
 
 class _testCase_dot(unittest.TestCase):
     def setUp(self):
-        core.build_eval_mode()
         numpy.random.seed(44)
-    def tearDown(self):
-        core.pop_mode()
-    def test_basic0(self):
-        for mtype in [sparse.csc_matrix, sparse.csr_matrix]:
-            x = SparseR(mtype(sparse.speye(5,3)))
-            y = core.wrap(numpy.random.rand(3, 2))
-
-            z = dot(x,y)
-            self.failUnless(z.data.shape == (5,2))
-            self.failUnless(type(z.data) is mtype)
-    def test_basic1(self):
-        """dot: sparse left"""
-        a = numpy.asarray([[1, 0, 3, 0, 5], [0, 0, -2, 0, 0]],
-                dtype='float64')
-        b = numpy.random.rand(5, 3)
-        for mtype in [sparse.csr_matrix, sparse.csc_matrix, sparse.dok_matrix,
-                sparse.lil_matrix]:#, sparse.coo_matrix]:
-            #print type(a), mtype
-            m = mtype(a)
-            ab = m.dot(b)
-            try:
-                z = dot(SparseR(m),core.ResultBase(data=b))
+
+    def test(self):
+        """Bring back the tests for sparse dot"""
+        raise NotImplementedError()
+    if 0:
+        def test_basic0(self):
+            for mtype in [sparse.csc_matrix, sparse.csr_matrix]:
+                x = assparse(mtype(sparse.speye(5,3)))
+                y = astensor(numpy.random.rand(3, 2))
+
+                z = dot(x,y)
+                self.failUnless(z.data.shape == (5,2))
+                self.failUnless(type(z.data) is mtype)
+
+        def test_basic1(self):
+            """dot: sparse left"""
+            a = numpy.asarray([[1, 0, 3, 0, 5], [0, 0, -2, 0, 0]],
+                    dtype='float64')
+            b = numpy.random.rand(5, 3)
+            for mtype in [sparse.csr_matrix, sparse.csc_matrix, sparse.dok_matrix,
+                    sparse.lil_matrix]:#, sparse.coo_matrix]:
+                #print type(a), mtype
+                m = mtype(a)
+                ab = m.dot(b)
+                try:
+                    z = dot(SparseR(m),core.ResultBase(data=b))
+                    self.failUnless(z.data.shape == ab.shape)
+                    self.failUnless(type(z.data) == type(ab))
+                except Exception, e:
+                    print 'cccc', mtype, e, str(e)
+                    raise
+
+        def test_basic2(self):
+            """dot: sparse right"""
+            a = numpy.random.rand(2, 5)
+            b = numpy.asarray([[1, 0, 3, 0, 5], [0, 0, -2, 0, 0]],
+                    dtype='float64').transpose()
+
+            for mtype in [sparse.csr_matrix, sparse.csc_matrix, sparse.dok_matrix,
+                    sparse.lil_matrix]:#, sparse.coo_matrix]:
+                m = mtype(b)
+                ab = m.transpose().dot(a.transpose()).transpose()
+                z = dot(core.ResultBase(data=a),SparseR(mtype(b)))
                 self.failUnless(z.data.shape == ab.shape)
                 self.failUnless(type(z.data) == type(ab))
-            except Exception, e:
-                print 'cccc', mtype, e, str(e)
-                raise
-    def test_basic2(self):
-        """dot: sparse right"""
-        a = numpy.random.rand(2, 5)
-        b = numpy.asarray([[1, 0, 3, 0, 5], [0, 0, -2, 0, 0]],
-                dtype='float64').transpose()
-
-        for mtype in [sparse.csr_matrix, sparse.csc_matrix, sparse.dok_matrix,
-                sparse.lil_matrix]:#, sparse.coo_matrix]:
-            m = mtype(b)
-            ab = m.transpose().dot(a.transpose()).transpose()
-            z = dot(core.ResultBase(data=a),SparseR(mtype(b)))
-            self.failUnless(z.data.shape == ab.shape)
-            self.failUnless(type(z.data) == type(ab))
-    def test_graph_bprop0(self):
-        x = core.wrap(numpy.random.rand(10,2))
-        w = SparseR(sparse.csr_matrix(numpy.asarray([[1, 0, 3, 0, 5], [0, 0, -2, 0,
-            0]],dtype='float64')))
-
-        for epoch in xrange(50):
-            xw = sparse2dense(dot(x, w))
-            y = sparse2dense(dot(xw, transpose(w)))
-            loss = core.sum(core.sqr(x-y))
-            gy = y-x
-            g = grad.Grad({y:gy})
-            g.bprop()
-            lr = 0.002
-            g(w).data[1,0] = 0
-            g(w).data[1,4] = 0
-            w.data = -lr * g(w).data + w.data
-
-        self.failUnless('3.08560636025' == str(loss.data))
-
-    def test_graph_bprop1(self):
-        x = core.wrap(numpy.random.rand(10,2))
-        w = SparseR(sparse.csr_matrix(numpy.asarray([[1, 0, 3, 0, 5], [0, 0, -2, 0,
-            0]],dtype='float64')))
-
-        for epoch in xrange(50):
-            xw = sparse2dense(dot(x, w))
-            y = sparse2dense(dot(xw, transpose(w)))
-            loss = core.sum(core.sqr(x-y))
-            g = grad.grad(loss)
-            lr = 0.001
-
-            g(w).data[1,0] = 0
-            g(w).data[1,4] = 0
-            w.data = -lr * g(w).data + w.data
-
-        self.failUnless('3.08560636025' == str(loss.data))
-        
+
+        def test_graph_bprop0(self):
+            x = core.wrap(numpy.random.rand(10,2))
+            w = SparseR(sparse.csr_matrix(numpy.asarray([[1, 0, 3, 0, 5], [0, 0, -2, 0,
+                0]],dtype='float64')))
+
+            for epoch in xrange(50):
+                xw = sparse2dense(dot(x, w))
+                y = sparse2dense(dot(xw, transpose(w)))
+                loss = core.sum(core.sqr(x-y))
+                gy = y-x
+                g = grad.Grad({y:gy})
+                g.bprop()
+                lr = 0.002
+                g(w).data[1,0] = 0
+                g(w).data[1,4] = 0
+                w.data = -lr * g(w).data + w.data
+
+            self.failUnless('3.08560636025' == str(loss.data))
+
+        def test_graph_bprop1(self):
+            x = core.wrap(numpy.random.rand(10,2))
+            w = SparseR(sparse.csr_matrix(numpy.asarray([[1, 0, 3, 0, 5], [0, 0, -2, 0,
+                0]],dtype='float64')))
+
+            for epoch in xrange(50):
+                xw = sparse2dense(dot(x, w))
+                y = sparse2dense(dot(xw, transpose(w)))
+                loss = core.sum(core.sqr(x-y))
+                g = grad.grad(loss)
+                lr = 0.001
+
+                g(w).data[1,0] = 0
+                g(w).data[1,4] = 0
+                w.data = -lr * g(w).data + w.data
+
+            self.failUnless('3.08560636025' == str(loss.data))
+    
 if __name__ == '__main__':
     unittest.main()
 

gof/op.py

@@ -144,7 +144,7 @@ class Op(object):
         TODO: consider moving this function to the python linker.
         """
         res = self.impl(*[input.data for input in self.inputs])
-        if self.nout == 1:
+        if len(self.outputs) == 1:
             self.outputs[0].data = res
         else:
             assert len(res) == len(self.outputs)

sparse.py

-
+import copy #for __copy__
 import numpy
 from scipy import sparse
 
-import gof
+import gof.op, gof.result
+import tensor
+
 
 # Wrapper type
 
-class SparseR(gof.ResultBase):
+def assparse(sp, **kwargs):
+    """Return SparseR version of sp"""
+    if isinstance(sp, SparseR):
+        return sp
+    else:
+        rval = SparseR(str(sp.dtype), sp.format, **kwargs)
+        rval.data = sp
+        return rval
+
+class SparseR(gof.result.ResultBase):
     """
     Attribute:
-    format - a subclass of sparse.spmatrix indicating self.data.__class__
+    format - a string identifying the type of sparsity
 
     Properties:
     T - read-only: return a transpose of self
@@ -19,126 +30,152 @@ class SparseR(gof.ResultBase):
     Notes:
 
     """
-    def __init__(self, data=None, role=None, constant = False, 
-            format = sparse.csr_matrix):
-        core.ResultBase.__init__(self, role, data, constant)
-        if isinstance(data, sparse.spmatrix):
-            self.format = data.__class__
-        else:
-            self.format = format
-        self._dtype = None
-        self._shape = None
+    format_cls = {
+            'csr' : sparse.csr_matrix,
+            'csc' : sparse.csc_matrix
+            }
+    dtype_set = set(['int', 'int32', 'int64', 'float32', 'float64'])
+
+    def __init__(self, dtype, format, **kwargs):
+        gof.ResultBase.__init__(self, **kwargs)
+        if dtype in SparseR.dtype_set:
+            self._dtype = dtype
+        assert isinstance(format, str)
 
-    def data_filter(self, value):
-        if isinstance(value, sparse.spmatrix): return value
-        return sparse.csr_matrix(value)
+        #print format, type(format), SparseR.format_cls.keys(), format in SparseR.format_cls
+        if format in SparseR.format_cls:
+            self._format = format
+        else:
+            raise NotImplementedError('unsupported format "%s" not in list' % format, SparseR.format_cls.keys())
+
+    def filter(self, value):
+        if isinstance(value, SparseR.format_cls[self.format])\
+                and value.dtype == self.dtype:
+                    return value
+        #print 'pass-through failed', type(value)
+        sp = SparseR.format_cls[self.format](value)
+        if str(sp.dtype) != self.dtype:
+            raise NotImplementedError()
+        if sp.format != self.format:
+            raise NotImplementedError()
+        return sp
+
+    def __copy__(self):
+        if self.name is not None:
+            rval = SparseR(self._dtype, self._format, name=self.name)
+        else:
+            rval = SparseR(self._dtype, self._format)
+        rval.data = copy.copy(self.data)
+        return rval
 
-    def __add__(left, right): return add(left, right)
-    def __radd__(right, left): return add(left, right)
 
+    dtype = property(lambda self: self._dtype)
+    format = property(lambda self: self._format)
     T = property(lambda self: transpose(self), doc = "Return aliased transpose")
 
-    # self._dtype is used when self._data hasn't been set yet
-    def __dtype_get(self):
-        if self._data is None:
-            return self._dtype
-        else:
-            return self._data.dtype
-    def __dtype_set(self, dtype):
-        if self._data is None:
-            self._dtype = dtype
-        else:
-            raise StateError('cannot set dtype after data has been set')
-    dtype = property(__dtype_get, __dtype_set)
 
-    # self._shape is used when self._data hasn't been set yet
-    def __shape_get(self):
-        if self._data is None:
-            return self._shape
-        else:
-            return self._data.shape
-    def __shape_set(self, shape):
-        if self._data is None:
-            self._shape = shape
-        else:
-            raise StateError('cannot set shape after data has been set')
-    shape = property(__shape_get, __shape_set)
+    def __add__(left, right): return add(left, right)
+    def __radd__(right, left): return add(left, right)
 
-# convenience base class
-class op(gof.PythonOp, grad.update_gradient_via_grad):
-    """unite PythonOp with update_gradient_via_grad"""
 
 #
 # Conversion
 #
 
 # convert a sparse matrix to an ndarray
-class sparse2dense(op):
-    def gen_outputs(self): return [core.Numpy2()]
-    def impl(x): return numpy.asarray(x.todense())
+class DenseFromSparse(gof.op.Op):
+    def __init__(self, x, **kwargs):
+        gof.op.Op.__init__(self, **kwargs)
+        self.inputs = [assparse(x)]
+        self.outputs = [tensor.Tensor(x.dtype,[0,0])]
+    def impl(self, x):
+        return numpy.asarray(x.todense())
     def grad(self, x, gz): 
-        if x.format is sparse.coo_matrix: return dense2coo(gz)
-        if x.format is sparse.csc_matrix: return dense2csc(gz)
-        if x.format is sparse.csr_matrix: return dense2csr(gz)
-        if x.format is sparse.dok_matrix: return dense2dok(gz)
-        if x.format is sparse.lil_matrix: return dense2lil(gz)
-
-# convert an ndarray to various sorts of sparse matrices.
-class _dense2sparse(op):
-    def gen_outputs(self): return [SparseR()]
-    def grad(self, x, gz): return sparse2dense(gz)
-class dense2coo(_dense2sparse):
-    def impl(x): return sparse.coo_matrix(x)
-class dense2csc(_dense2sparse):
-    def impl(x): return sparse.csc_matrix(x)
-class dense2csr(_dense2sparse):
-    def impl(x): return sparse.csr_matrix(x)
-class dense2dok(_dense2sparse):
-    def impl(x): return sparse.dok_matrix(x)
-class dense2lil(_dense2sparse):
-    def impl(x): return sparse.lil_matrix(x)
-
+        return sparse_from_dense(gz, x.format)
+dense_from_sparse = gof.op.constructor(DenseFromSparse)
+
+class SparseFromDense(gof.op.Op):
+    def __init__(self, x, format, **kwargs):
+        gof.op.Op.__init__(self, **kwargs)
+        if isinstance(format, gof.result.ResultBase):
+            self.inputs = [tensor.astensor(x), format]
+        else:
+            self.inputs =  [tensor.astensor(x), gof.result.PythonResult()]
+            self.inputs[1].data = format
+        self.outputs = [SparseR(x.dtype, self.inputs[1].data)]
+    def impl(self, x, fmt):
+        # this would actually happen anyway when we try to assign to
+        # self.outputs[0].data, but that seems hackish -JB
+        return SparseR.format_cls[fmt](x)
+    def grad(self, (x, fmt), gz):
+        return dense_from_sparse(gz)
+sparse_from_dense = gof.op.constructor(SparseFromDense)
 
 # Linear Algebra
 
-class add(op):
-    def gen_outputs(self): return [SparseR()]
-    def impl(csr,y): return csr + y
-
-class transpose(op):
-    def gen_outputs(self): return [SparseR()]
-    def impl(x): return x.transpose() 
-    def grad(self, x, gz): return transpose(gz)
-
-class dot(op):
-    """
-    Attributes:
-    grad_preserves_dense - an array of boolean flags (described below)
-
+class Transpose(gof.op.Op):
+    format_map = {
+            'csr' : 'csc',
+            'csc' : 'csr'}
+    def __init__(self, x, **kwargs):
+        gof.op.Op.__init__(self, **kwargs)
+        x = assparse(x)
+        self.inputs = [x]
+        self.outputs = [SparseR(x.dtype, Transpose.format_map[x.format])]
+    def impl(self, x):
+        return x.transpose() 
+    def grad(self, x, gz): 
+        return transpose(gz)
+transpose = gof.op.constructor(Transpose)
+
+class AddSS(gof.op.Op): #add two sparse matrices
+    def __init__(self, x, y, **kwargs):
+        gof.op.Op.__init__(self, **kwargs)
+        x, y = [assparse(x), assparse(y)]
+        self.inputs = [x, y]
+        if x.dtype != y.dtype:
+            raise NotImplementedError()
+        if x.format != y.format:
+            raise NotImplementedError()
+        self.outputs = [SparseR(x.dtype, x.format)]
+    def impl(self, x,y): 
+        return x + y
+    def grad(self, (x, y), gz):
+        return gz, gz
+add_s_s = gof.op.constructor(AddSS)
+
+
+if 0:
+    class dot(gof.op.Op):
+        """
+        Attributes:
+        grad_preserves_dense - an array of boolean flags (described below)
+
+
+        grad_preserves_dense controls whether gradients with respect to inputs are
+        converted to dense matrices when the corresponding inputs are not in a
+        SparseR wrapper.  This can be a good idea when dot is in the middle of a
+        larger graph, because the types of gx and gy will match those of x and y.
+        This conversion might be annoying if the gradients are graph outputs though,
+        hence this mask.
+        """
+        def __init__(self, *args, **kwargs):
+            gof.op.Op.__init__(self, **kwargs)
+            self.grad_preserves_dense = [True, True]
+        def gen_outputs(self): return [SparseR()]
+        def impl(x,y):
+            if hasattr(x, 'getnnz'):
+                return x.dot(y)
+            else:
+                return y.transpose().dot(x.transpose()).transpose()
+
+        def grad(self, x, y, gz):
+            rval = [dot(gz, y.T), dot(x.T, gz)]
+            for i in 0,1:
+                if not isinstance(self.inputs[i], SparseR):
+                    #assume it is a dense matrix
+                    if self.grad_preserves_dense[i]:
+                        rval[i] = dense_from_sparse(rval[i])
+            return rval
 
-    grad_preserves_dense controls whether gradients with respect to inputs are
-    converted to dense matrices when the corresponding inputs are not in a
-    SparseR wrapper.  This can be a good idea when dot is in the middle of a
-    larger graph, because the types of gx and gy will match those of x and y.
-    This conversion might be annoying if the gradients are graph outputs though,
-    hence this mask.
-    """
-    def __init__(self, *args, **kwargs):
-        op.__init__(self, *args, **kwargs)
-        self.grad_preserves_dense = [True, True]
-    def gen_outputs(self): return [SparseR()]
-    def impl(x,y):
-        if hasattr(x, 'getnnz'):
-            return x.dot(y)
-        else:
-            return y.transpose().dot(x.transpose()).transpose()
-
-    def grad(self, x, y, gz):
-        rval = [dot(gz, y.T), dot(x.T, gz)]
-        for i in 0,1:
-            if not isinstance(self.inputs[i], SparseR):
-                #assume it is a dense matrix
-                if self.grad_preserves_dense[i]:
-                    rval[i] = sparse2dense(rval[i])
-        return rval
 

tensor.py

@@ -70,6 +70,8 @@ class Tensor(BaseTensor):
 # alternate Tensor constructor
 def astensor(data, broadcastable=None, role=None, name=None):
     """Return a Tensor containing given data"""
+    if isinstance(data, Tensor) and broadcastable is None and role is None and name is None:
+        return data
     data = numpy.asarray(data)
     if broadcastable is None:
         broadcastable = [s==1 for s in data.shape]
@@ -116,11 +118,9 @@ def _assert_tensor_scalar(x, a):
     if numpy.product(a.shape) != 1:
         raise ValueError("The second argument must be a scalar.")
 
-def _as_tensor(obj):
-    if isinstance(obj, Tensor):
-        return obj
-    else:
-        return astensor(obj)
+# this has a different name, because _as_tensor is the function which ops use
+# to upcast their arguments... this internal-use function is a good place to put debugging stuff, better than the global astensor.
+_as_tensor = astensor
 
 class _Op(BaseTensorOp):
     """A convenient base for the ops in this file"""