Maybe sparse.py now works?

_test_sparse.py

@@ -80,88 +80,86 @@ class _testCase_dot(unittest.TestCase):
     def setUp(self):
         numpy.random.seed(44)
 
-    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.Result(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.Result(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_basic0(self):
+        for mtype in [sparse.csc_matrix, sparse.csr_matrix]:
+            x = assparse(mtype(sparse.speye(5,3)))
+            y = tensor.astensor(numpy.random.rand(3, 2))
+
+            zop = dot(x,y)
+            z = compile.eval_outputs([zop])
+            self.failUnless(z.shape == (5,2))
+            self.failUnless(type(z) 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(assparse(m), gof.Result(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(gof.Result(data=a),assparse(mtype(b)))
+#            self.failUnless(z.data.shape == ab.shape)
+#            self.failUnless(type(z.data) == type(ab))
+#
+#    def test_graph_bprop0(self):
+#        x = tensor.astensor(numpy.random.rand(10,2))
+#        w = assparse(sparse.csr_matrix(
+#                numpy.asarray([[1, 0, 3, 0, 5], [0, 0, -2, 0,0]],dtype='float64')
+#            ))
+#
+#        for epoch in xrange(50):
+#            xw = dense_from_sparse(dot(x, w))
+#            y = dense_from_sparse(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 = tensor.astensor(numpy.random.rand(10,2))
+#        w = assparse(sparse.csr_matrix(
+#                numpy.asarray([[1, 0, 3, 0, 5], [0, 0, -2, 0,0]],dtype='float64')
+#            ))
+#
+#        for epoch in xrange(50):
+#            xw = dense_from_sparse(dot(x, w))
+#            y = dense_from_sparse(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

@@ -131,15 +131,19 @@ class Op(object):
         this L{Op}'s inputs.
 
         To do a bottom-up copy of a graph, use clone_with_new_inputs.
+
+        @attention: If your L{Op} has additional options or a different
+        constructor you probably want to override this.
         """
         return self.__class__(*self.inputs)
 
     def clone_with_new_inputs(self, *new_inputs):
         """
         Returns a clone of this L{Op} that takes different inputs. The
-        default behavior is to call the constructor on the new inputs,
-        but if your L{Op} has additional options or a different constructor
-        you might want to override this.
+        default behavior is to call the constructor on the new inputs.
+
+        @attention: If your L{Op} has additional options or a different
+        constructor you probably want to override this.
         """
         return self.__class__(*new_inputs)
 

sparse.py

+"""
+Classes for handling sparse matrices.
+
+To read about different sparse formats, see U{http://www-users.cs.umn.edu/~saad/software/SPARSKIT/paper.ps}.
+
+@todo Automatic methods for determining best sparse format?
+"""
+
 import copy #for __copy__
 import numpy
 from scipy import sparse
@@ -14,10 +22,14 @@ def assparse(sp, **kwargs):
     @param sp:  A sparse matrix. assparse reads dtype and format properties
                 out of this sparse matrix.
     @return:    SparseR version of sp.
+
+    @todo Verify that sp is sufficiently sparse, and raise a warning if it is not
     """
     if isinstance(sp, SparseR):
         return sp
     else:
+        # @todo Verify that sp is sufficiently sparse, and raise a
+        # warning if it is not
         rval = SparseR(str(sp.dtype), sp.format, **kwargs)
         rval.data = sp
         return rval
@@ -156,43 +168,65 @@ class AddSS(gof.op.Op): #add two sparse matrices
         return gz, gz
 add_s_s = gof.op.constructor(AddSS)
 
-#class Dot(gof.op.Op):
-#    def __init__(self, x, y):
-#        self.inputs = [x, y]    # Need to convert? e.g. _as_tensor
-#        # broadcastable
-#    def perform:
-#        #return numpy.dot(x, y)
-#    def grad:
-#
-#    """
-#    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'):
-#            # if x is sparse, then do this.
-#            return x.dot(y)
-#        else:
-#            # if x is dense (and y is sparse), we do this
-#            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
+class Dot(gof.op.Op):
+    """
+    Attributes:
+    grad_preserves_dense - a boolean flags [default: True].
+    grad_preserves_dense controls whether gradients with respect to inputs
+    are converted to dense matrices when the corresponding input y is
+    dense (not in a L{SparseR} wrapper). This is generally a good idea
+    when L{Dot} is in the middle of a larger graph, because the types
+    of gy will match that of y. This conversion might be inefficient if
+    the gradients are graph outputs though, hence this mask.
+    """
+    def __init__(self, x, y, grad_preserves_dense=True):
+        """
+        Because of trickiness of implementing, we assume that the left argument x is SparseR (not dense)
+        """
+        if x.dtype != y.dtype:
+            raise NotImplementedError()
+
+        # These are the conversions performed by scipy.sparse.dot
+        if x.format == "csc" or x.format == "coo":
+            myformat = "csc"
+        elif x.format == "csr":
+            myformat = "csr"
+        else:
+            raise NotImplementedError()
+
+        self.inputs = [x, y]    # Need to convert? e.g. assparse
+        self.outputs = [SparseR(x.dtype, myformat)]
+        self.grad_preserves_dense = grad_preserves_dense
+    def perform(self):
+        """
+        @todo Verify that output is sufficiently sparse, and raise a warning if it is not
+        @todo Also determine that we are storing the output in the best storage format?
+        """
+        self.outputs[0].data = self.inputs[0].data.dot(self.inputs[1].data)
+    def grad(self, (x, y), (gz,)):
+        rval = [dot(gz, y.T), dot(x.T, gz)]
+        assert isinstance(self.inputs[0], SparseR)
+        if not isinstance(self.inputs[1], SparseR):
+            if self.grad_preserves_dense:
+                rval[1] = dense_from_sparse(rval[1])
+        return rval
+    def __copy__(self):
+        return self.__class__(self.inputs[0], self.inputs[1], self.grad_preserves_dense)
+    def clone_with_new_inputs(self, *new_inputs):
+        return self.__class__(new_inputs[0], new_inputs[1], self.grad_preserves_dense)
+def dot(x, y, grad_preserves_dense=True):
+    """
+    @todo Maybe the triple-transposition formulation (when x is dense)
+    is slow. See if there is a direct way to do this.
+    """
+    if hasattr(x, 'getnnz'): x = assparse(x)
+    if hasattr(y, 'getnnz'): y = assparse(y)
+
+    x_is_sparse = isinstance(x, SparseR)
+    y_is_sparse = isinstance(y, SparseR)
+    if not x_is_sparse and not y_is_sparse:
+        raise TypeError()
+    if x_is_sparse:
+        return Dot(x,y,grad_preserves_dense).outputs[0]
+    else:
+        return transpose(Dot(transpose(y), transpose(x), grad_preserves_dense).outputs[0])

tensor.py

@@ -351,6 +351,9 @@ class Dot(_Op):
     def impl(self, x, y):
         return numpy.dot(x, y)
     def grad(self, (x, y), gz):
+        """
+        @todo Shouldn't it be (gz,) ? -jpt
+        """
         return dot(gz, y.T), dot(x.T, gz)
     if 0:
         def c_support_code(self):