Add tests and rewrite SpSum.

theano/sparse/basic.py

@@ -28,7 +28,6 @@ def register_specialize(lopt, *tags, **kwargs):
                                          lopt.__name__, lopt, 'fast_run',
                                          *tags)
 
-
 """ Types of sparse matrices to use for testing """
 _mtypes = [scipy.sparse.csc_matrix, scipy.sparse.csr_matrix]
 #_mtypes = [sparse.csc_matrix, sparse.csr_matrix, sparse.dok_matrix,
@@ -354,6 +353,7 @@ class SparseConstant(gof.Constant, _sparse_py_operators):
     def __repr__(self):
         return str(self)
 
+
 class SparseType(gof.Type):
     """
     @type dtype: numpy dtype string such as 'int64' or 'float64' (among others)
@@ -1286,84 +1286,135 @@ class Neg(gof.op.Op):
 neg = Neg()
 
 
-class SpSum(gof.op.Op):
-    """TODO: rewrite
-    Scale each columns of a sparse matrix by the
-    corresponding element of a dense vector
+class ColScaleCSC(gof.op.Op):
     """
+    Scale each columns of a sparse matrix by the corresponding element
+    of a dense vector
+    """
+
+    def make_node(self, x, s):
+        if x.format != 'csc':
+            raise ValueError('x was not a csc matrix')
+        return gof.Apply(self, [x, s], [x.type()])
+
+    def perform(self, node, (x, s), (z,)):
+        M, N = x.shape
+        assert x.format == 'csc'
+        assert s.shape == (N,)
+
+        y = x.copy()
+
+        for j in xrange(0, N):
+            y.data[y.indptr[j]: y.indptr[j + 1]] *= s[j]
 
-    axis = None
-    sparse_grad = False
+        z[0] = y
 
-    def __init__(self, axis, sparse_grad=True):
+    def grad(self, (x, s), (gz,)):
+        return [col_scale(gz, s), sp_sum(x * gz, axis=0)]
+
+
+class RowScaleCSC(gof.op.Op):
+    """
+    Scale each row of a sparse matrix by the corresponding element of
+    a dense vector
     """
-        :param sparse_grad: if True, this instance ignores the
-            gradient on matrix elements which are implicitly 0.
+
+    def make_node(self, x, s):
+        return gof.Apply(self, [x, s], [x.type()])
+
+    def perform(self, node, (x, s), (z,)):
+        M, N = x.shape
+        assert x.format == 'csc'
+        assert s.shape == (M,)
+
+        indices = x.indices
+        indptr = x.indptr
+
+        y_data = x.data.copy()
+
+        for j in xrange(0, N):
+            for i_idx in xrange(indptr[j], indptr[j + 1]):
+                y_data[i_idx] *= s[indices[i_idx]]
+
+        z[0] = scipy.sparse.csc_matrix((y_data, indices, indptr), (M, N))
+
+    def grad(self, (x, s), (gz,)):
+        return [row_scale(gz, s), sp_sum(x * gz, axis=1)]
+
+
+def col_scale(x, s):
+    if x.format == 'csc':
+        return ColScaleCSC()(x, s)
+    elif x.format == 'csr':
+        return RowScaleCSC()(x.T, s).T
+    else:
+        raise NotImplementedError()
+
+
+def row_scale(x, s):
+    return col_scale(x.T, s).T
+
+
+class SpSum(gof.op.Op):
+    """Calculate the sum of a sparse matrix along a specify
+    axis.
+
+    It operates a reduction along the axis specified. When
+    `axis` is `None`, it is apply along all axis.
+
+    :param x: Sparse matrix.
+    :param axis: Axis along the sum is apply. Integers or `None`.
+    :param sparse_grad: `True` to have a structured grad. Boolean.
+
+    :return: The sum of `x` in a dense format.
+
+    :note:
+    - The grad implementation is controlled with the `sparse_grad`
+      parameter. `True` will provide a structured grad and `False`
+      will provide a regular grad.
+    - This op does not return a sparse matrix.
     """
+
+    def __init__(self, axis=None, sparse_grad=False):
         super(SpSum, self).__init__()
         self.axis = axis
-        self.sparse_grad = sparse_grad
+        self.structured = sparse_grad
         if self.axis not in (None, 0, 1):
-            raise ValueError('illegal value for self.axis')
+            raise ValueError('Illegal value for self.axis.')
 
     def __eq__(self, other):
-        #WARNING: judgement call...
-        #not using the sparse_grad in the comparison or hashing
-        #because it doesn't change the perform method therefore, we
-        #*do* want Sums with different sparse_grad values to be merged
-        #by the merge optimization.
-        # This requires them to compare equal.
+        # WARNING: judgement call...
+        # We are not using the structured in the comparison or hashing
+        # because it doesn't change the perform method therefore, we
+        # *do* want Sums with different structured values to be merged
+        # by the merge optimization and this requires them to compare equal.
         return type(self) == type(other) and self.axis == other.axis
 
     def __hash__(self):
+        # WARNING: judgement call...
+        # We are not using the structured in the comparison or hashing
+        # because it doesn't change the perform method therefore, we
+        # *do* want Sums with different structured values to be merged
+        # by the merge optimization and this requires them to compare equal.
         return 76324 ^ hash(type(self)) ^ hash(self.axis)
 
-    def __str__(self):
-        return self.__class__.__name__ + "{axis=%s}" % str(self.axis)
-
     def make_node(self, x):
-        ###
-        # At least for small matrices (5x5), the .sum() method of a
-        # csc matrix returns a dense matrix as the result whether axis
-        # is 0 or 1... weird!
-        ###
-        assert isinstance(x.type, theano.sparse.SparseType)
+        x = as_sparse_variable(x)
         b = ()
         if self.axis is not None:
             b = (False,)
-        z = tensor.tensor(broadcastable=b, dtype=x.dtype)
-        return gof.Apply(self, [x], [z])
 
-    def infer_shape(self, node, shapes):
-        r = None
-        if self.axis is None:
-            r = [()]
-        elif self.axis == 0:
-            r = [(shapes[0][1],)]
-        else:
-            r = [(shapes[0][0],)]
-        return r
+        z = tensor.TensorType(broadcastable=b, dtype=x.dtype)()
+        return gof.Apply(self, [x], [z])
 
     def perform(self, node, (x,), (z,)):
-        if self.axis is None:
+        if self.axis == None:
             z[0] = numpy.asarray(x.sum())
         else:
-            if self.axis == 0:
-                if x.format == 'csc':
-                    z[0] = numpy.asarray(x.sum(axis=self.axis)).reshape(
-                        (x.shape[1], ))
-                else:
-                    z[0] = numpy.asarray(x.asformat(x.format).sum(
-                            axis=self.axis)).reshape((x.shape[1],))
-            elif self.axis == 1:
-                if x.format == 'csr':
-                    z[0] = numpy.asarray(x.sum(axis=self.axis)).reshape(
-                        (x.shape[0],))
-                else:
-                    z[0] = numpy.asarray(x.asformat(x.format).sum(
-                            axis=self.axis)).reshape((x.shape[0],))
+            z[0] = numpy.asarray(x.sum(self.axis)).ravel()
 
     def grad(self, (x,), (gz,)):
+        if self.structured:
             if self.axis is None:
                 r = gz * theano.sparse.sp_ones_like(x)
             elif self.axis == 0:
@@ -1371,12 +1422,24 @@ class SpSum(gof.op.Op):
             elif self.axis == 1:
                 r = row_scale(theano.sparse.sp_ones_like(x), gz)
             else:
-            assert False
+                raise ValueError('Illegal value for self.axis.')
+        else:
+            # TODO
+            raise NotImplementedError()
+        return [r]
 
-        if not self.sparse_grad:
-            r = theano.sparse.dense_from_sparse(r)
+    def infer_shape(self, node, shapes):
+        r = None
+        if self.axis is None:
+            r = [()]
+        elif self.axis == 0:
+            r = [(shapes[0][1],)]
+        else:
+            r = [(shapes[0][0],)]
+        return r
 
-        return [r]
+    def __str__(self):
+        return self.__class__.__name__ + "{axis=%s}" % str(self.axis)
 
 
 def sp_sum(x, axis=None, sparse_grad=False):

theano/sparse/sandbox/sp.py

@@ -18,7 +18,9 @@ from theano.gof.python25 import all, any
 from theano.sparse.basic import Remove0, remove0
 
 # To maintain compatibility
-from theano.sparse import SpSum, sp_sum
+from theano.sparse import (
+    SpSum, sp_sum,
+    ColScaleCSC, RowScaleCSC, col_scale, row_scale)
 
 
 def register_specialize(lopt, *tags, **kwargs):
@@ -117,75 +119,6 @@ class SquareDiagonal(Op):
 square_diagonal = SquareDiagonal()
 
 
-class ColScaleCSC(Op):
-    """
-    Scale each columns of a sparse matrix by the corresponding element
-    of a dense vector
-    """
-
-    def make_node(self, x, s):
-        if x.format != 'csc':
-            raise ValueError('x was not a csc matrix')
-        return gof.Apply(self, [x, s], [x.type()])
-
-    def perform(self, node, (x, s), (z,)):
-        M, N = x.shape
-        assert x.format == 'csc'
-        assert s.shape == (N,)
-
-        y = x.copy()
-
-        for j in xrange(0, N):
-            y.data[y.indptr[j]: y.indptr[j + 1]] *= s[j]
-
-        z[0] = y
-
-    def grad(self, (x, s), (gz,)):
-        return [col_scale(gz, s), sp_sum(x * gz, axis=0)]
-
-
-class RowScaleCSC(Op):
-    """
-    Scale each row of a sparse matrix by the corresponding element of
-    a dense vector
-    """
-
-    def make_node(self, x, s):
-        return gof.Apply(self, [x, s], [x.type()])
-
-    def perform(self, node, (x, s), (z,)):
-        M, N = x.shape
-        assert x.format == 'csc'
-        assert s.shape == (M,)
-
-        indices = x.indices
-        indptr = x.indptr
-
-        y_data = x.data.copy()
-
-        for j in xrange(0, N):
-            for i_idx in xrange(indptr[j], indptr[j + 1]):
-                y_data[i_idx] *= s[indices[i_idx]]
-
-        z[0] = scipy_sparse.csc_matrix((y_data, indices, indptr), (M, N))
-
-    def grad(self, (x, s), (gz,)):
-        return [row_scale(gz, s), sp_sum(x * gz, axis=1)]
-
-
-def col_scale(x, s):
-    if x.format == 'csc':
-        return ColScaleCSC()(x, s)
-    elif x.format == 'csr':
-        return RowScaleCSC()(x.T, s).T
-    else:
-        raise NotImplementedError()
-
-
-def row_scale(x, s):
-    return col_scale(x.T, s).T
-
-
 class EnsureSortedIndices(Op):
     """
     Remove explicit zeros from a sparse matrix, and resort indices

theano/sparse/tests/test_basic.py

@@ -72,6 +72,49 @@ def random_lil(shape, dtype, nnz):
     return rval
 
 
+def sparse_random_inputs(format, shape, n=1, out_dtype=None, p=0.5):
+    """Return a tuple containing everything needed to
+    perform a test.
+
+    If `out_dtype` is `None`, theano.config.floatX is
+    used.
+
+    :param format: Sparse format.
+    :param shape: Shape of data.
+    :param n: Number of variable.
+    :param out_dtype: dtype of output.
+    :param p: Sparsity proportion.
+
+    :return: (variable, data) where both `variable`
+             and `data` are list.
+    """
+
+    if out_dtype is None:
+        out_dtype = theano.config.floatX
+
+    assert 0 <= p and p <= 1
+    assert len(shape) == 2
+    assert out_dtype in sparse.all_dtypes
+
+    variable = [getattr(theano.sparse, format + '_matrix')(dtype=out_dtype)
+                for k in range(n)]
+
+    def _rand():
+        where = numpy.random.binomial(1, p, size=shape).astype('int8')
+
+        if out_dtype in sparse.discrete_dtypes:
+            value = numpy.random.randint(20, size=shape).astype(out_dtype)
+        else:
+            value = numpy.random.random(shape)
+
+        return where * value
+
+    data = [getattr(scipy.sparse, format + '_matrix')(_rand())
+            for k in range(n)]
+
+    return (variable, data)
+
+
 class T_verify_grad_sparse(unittest.TestCase):
     class FailOp(gof.op.Op):
         def __init__(self, structured):
@@ -1329,71 +1372,51 @@ def test_size():
         check()
 
 
-def test_sp_sum():
-    from theano.sparse import SpSum
+class SpSumTester(utt.InferShapeTester):
+    possible_axis = [None, 0, 1]
 
-    # TODO: test both grad.
-    rng = numpy.random.RandomState(42)
-    from theano.sparse.basic import SparseFromDense,DenseFromSparse
-    cases = [("csc", scipy.sparse.csc_matrix), ("csr", scipy.sparse.csr_matrix)]
-
-    for format, cast in cases:
+    def setUp(self):
+        super(SpSumTester, self).setUp()
+        self.op_class = sparse.SpSum
+        self.op = sparse.sp_sum
 
-        #print 'format: %(format)s' % locals()
-        x = theano.sparse.SparseType(format=format,
-                                     dtype=theano.config.floatX)()
-        x_data = numpy.arange(20).reshape(5,4).astype(theano.config.floatX)
+    def test_op(self):
+        for format in sparse.sparse_formats:
+            for axis in self.possible_axis:
+                variable, data = sparse_random_inputs(format,
+                                                      shape=(10, 10))
 
-        # Sum on all axis
-        #print 'sum on all axis...'
-        z = theano.sparse.sp_sum(x)
+                z = theano.sparse.sp_sum(*variable, axis=axis)
+                if axis == None:
                     assert z.type.broadcastable == ()
-        f = theano.function([x], z)
-        x_val = cast(x_data)
-        out = f(x_val)
-        expected = x_val.sum()
-        assert out == expected
-
-        # Sum on axis 0
-        #print 'sum on axis 0...'
-        z = theano.sparse.sp_sum(x, axis=0)
-        assert z.type.broadcastable == (False,)
-        f = theano.function([x], z)
-        x_val = cast(x_data)
-        out = f(x_val)
-        expected = x_val.sum(axis=0)
-        assert (out == expected).all()
-
-        # Sum on axis 1
-        #print 'sum on axis 1...'
-        z = theano.sparse.sp_sum(x, axis=1)
-        assert z.type.broadcastable == (False,)
-        f = theano.function([x], z)
-        x_val = cast(x_data)
-        out = f(x_val)
-        expected = numpy.asarray(x_val.sum(axis=1)).reshape(x_val.shape[0])
-        assert (out == expected).all()
-
-        # Sparse gradient on Sum on all axis
-        # unfinished, and suspended until verify_grad get fixed
-        if False:
-            #                print 'grad on sum on all axis...'
-            def fun(x):
-                ## verify_grad does not handle sparse data, so here's some casting as a workaround.
-                # x is a dense matrix: make it sparse
-                sparse_var = SparseFromDense(format)(x)
-                # apply op
-                dense_sum = theano.sparse.SpSum(axis=None, sparse_grad=False)(sparse_var)
-                return dense_sum
-                # cast back to dense so that verify_grad can work
-                dense_sum = theano.sparse.DenseFromSparse()(sparse_sum)
-                return dense_sum
-            x_val = x_data.copy()
-            #                print type(x_val)
-            import pdb;pdb.set_trace()
-            tensor.verify_grad(fun, [x_val], rng=rng)
-            #utt.verify_grad(SpSum(axis=None), [x_val])
-            #                print 'ok'
+                else:
+                    assert z.type.broadcastable == (False, )
+
+                f = theano.function(variable, self.op(*variable, axis=axis))
+                tested = f(*data)
+                expected = data[0].todense().sum(axis).ravel()
+                assert numpy.allclose(tested, expected)
+
+    def test_infer_shape(self):
+        for format in sparse.sparse_formats:
+            for axis in self.possible_axis:
+                variable, data = sparse_random_inputs(format,
+                                                      shape=(10, 10))
+                self._compile_and_check(variable,
+                                        [self.op(*variable, axis=axis)],
+                                        data,
+                                        self.op_class)
+
+    def test_grad(self):
+        for format in sparse.sparse_formats:
+            for axis in self.possible_axis:
+                for struct in [True]:
+                    variable, data = sparse_random_inputs(format,
+                                                          shape=(10, 10))
+                    verify_grad_sparse(
+                        self.op_class(axis=axis, sparse_grad=struct),
+                        data,
+                        structured=struct)
 
 
 class Remove0Tester(utt.InferShapeTester):