Merge pull request #1000 from lamblin/fix_preallocated_output

theano/sparse/basic.py

@@ -2316,158 +2316,6 @@ class Remove0(gof.Op):
 remove0 = Remove0()
 
 
-# Probability
-class Poisson(gof.op.Op):
-    """Return a sparse having random values from a Poisson density
-    with mean from the input.
-
-    :param x: Sparse matrix.
-
-    :return: A sparse matrix of random integers of a Poisson density
-             with mean of `x` element wise.
-    """
-
-    def __eq__(self, other):
-        return (type(self) == type(other))
-
-    def __hash__(self):
-        return hash(type(self))
-
-    def make_node(self, x):
-        x = as_sparse_variable(x)
-        return gof.Apply(self, [x], [x.type()])
-
-    def perform(self, node, (x, ), (out, )):
-        assert _is_sparse(x)
-        out[0] = x.copy()
-        out[0].data = numpy.asarray(numpy.random.poisson(out[0].data),
-                                    dtype=x.dtype)
-        out[0].eliminate_zeros()
-
-    def grad(self, inputs, outputs_gradients):
-        return [None]
-
-    def infer_shape(self, node, ins_shapes):
-        return ins_shapes
-
-    def __str__(self):
-        return self.__class__.__name__
-poisson = Poisson()
-
-
-class Binomial(gof.op.Op):
-    """Return a sparse matrix having random values from a binomial
-    density having number of experiment `n` and probability of succes
-    `p`.
-
-    :param n: Tensor scalar representing the number of experiment.
-    :param p: Tensor scalar representing the probability of success.
-    :param shape: Tensor vector for the output shape.
-
-    :return: A sparse matrix of integers representing the number
-             of success.
-    """
-
-    def __init__(self, format, dtype):
-        self.format = format
-        self.dtype = dtype
-
-    def __eq__(self, other):
-        return ((type(self) == type(other)) and
-                self.format == other.format and
-                self.dtype == other.dtype)
-
-    def __hash__(self):
-        return hash(type(self)) ^ hash(self.format) ^ hash(self.dtype)
-
-    def make_node(self, n, p, shape):
-        n = tensor.as_tensor_variable(n)
-        p = tensor.as_tensor_variable(p)
-        shape = tensor.as_tensor_variable(shape)
-        return gof.Apply(self, [n, p, shape], [SparseType(dtype=self.dtype,
-                                 format=self.format).make_variable()])
-
-    def perform(self, node, (n, p, shape, ), (out, )):
-        binomial = numpy.random.binomial(n, p, size=shape)
-        csx_matrix = getattr(scipy.sparse, self.format + '_matrix')
-        out[0] = csx_matrix(binomial, dtype=self.dtype)
-
-    def grad(self, (n, p, shape, ), (gz,)):
-        return None, None, None
-
-    def infer_shape(self, node, ins_shapes):
-        return [(node.inputs[2][0], node.inputs[2][1])]
-
-    def __str__(self):
-        return self.__class__.__name__
-
-csr_fbinomial = Binomial('csr', 'float32')
-csc_fbinomial = Binomial('csc', 'float32')
-csr_dbinomial = Binomial('csr', 'float64')
-csc_dbinomial = Binomial('csc', 'float64')
-
-
-class Multinomial(gof.op.Op):
-    """Return a sparse matrix having random values from a multinomial
-    density having number of experiment `n` and probability of succes
-    `p`.
-
-    :param n: Tensor type vector or scalar representing the number of
-              experiment for each row. If `n` is a scalar, it will be
-              used for each row.
-    :param p: Sparse matrix of probability where each row is a probability
-              vector representing the probability of succes. N.B. Each row
-              must sum to one.
-
-    :return: A sparse matrix of random integers from a multinomial density
-             for each row.
-
-    :note: It will works only if `p` have csr format.
-    """
-
-    def __eq__(self, other):
-        return (type(self) == type(other))
-
-    def __hash__(self):
-        return hash(type(self))
-
-    def make_node(self, n, p):
-        n = tensor.as_tensor_variable(n)
-        p = as_sparse_variable(p)
-
-        return gof.Apply(self, [n, p], [p.type()])
-
-    def perform(self, node, (n, p), (out, )):
-        assert _is_sparse(p)
-
-        if p.format != 'csr':
-            raise NotImplemented()
-
-        out[0] = p.copy()
-
-        if n.ndim == 0:
-            for i in xrange(p.shape[0]):
-                k, l = p.indptr[i], p.indptr[i + 1]
-                out[0].data[k:l] = numpy.random.multinomial(n, p.data[k:l])
-        elif n.ndim == 1:
-            if n.shape[0] != p.shape[0]:
-                raise ValueError('The number of element of n must be '
-                                 'the same as the number of row of p.')
-            for i in xrange(p.shape[0]):
-                k, l = p.indptr[i], p.indptr[i + 1]
-                out[0].data[k:l] = numpy.random.multinomial(n[i], p.data[k:l])
-
-    def grad(self, inputs, outputs_gradients):
-        return [None, None]
-
-    def infer_shape(self, node, ins_shapes):
-        return [ins_shapes[1]]
-
-    def __str__(self):
-        return self.__class__.__name__
-multinomial = Multinomial()
-
-
 # Structured monoid
 def structured_monoid(tensor_op):
     # Generic operation to perform many kinds of monoid element-wise

theano/sparse/opt.py

@@ -885,7 +885,7 @@ class MulSDCSC(gof.Op):
                                [tensor.tensor(b.dtype, (False,))])
 
     def c_code_cache_version(self):
-        return (1,)
+        return (2,)
 
     #def perform(self, node, (a_data, a_indices, a_indptr, b), (out,)):
     #    return NotImplementedError()
@@ -918,18 +918,20 @@ class MulSDCSC(gof.Op):
         if( PyArray_DESCR(%(_indptr)s)->type_num != NPY_INT32)
         {PyErr_SetString(PyExc_NotImplementedError, "D"); %(fail)s;}
 
-        if (!%(_zout)s)
+        if (!%(_zout)s ||
+            (PyArray_DIMS(%(_zout)s)[0] != PyArray_DIMS(%(_indices)s)[0]) ||
+            !(PyArray_ISCONTIGUOUS(%(_zout)s)))
         {
+            Py_XDECREF(%(_zout)s);
             %(_zout)s = (PyArrayObject*) PyArray_SimpleNew(1,
                   PyArray_DIMS(%(_indices)s), PyArray_DESCR(%(_b)s)->type_num);
-        }
-
-        if (PyArray_DIMS(%(_zout)s)[0] != PyArray_DIMS(%(_indices)s)[0])
+            if (!%(_zout)s)
             {
-            PyErr_SetString(PyExc_NotImplementedError,
-    "somehow _zout got the wrong size.. and I don't know how to resize it.");
+                PyErr_SetString(PyExc_MemoryError,
+                    "Could not allocate output memory.");
                 %(fail)s;
             }
+        }
 
         { //makes it compile even though labels jump over variable definitions.
             const npy_intp nnz = PyArray_DIMS(%(_indices)s)[0];
@@ -999,7 +1001,7 @@ class MulSDCSR(gof.Op):
                                [tensor.tensor(b.dtype, (False,))])
 
     def c_code_cache_version(self):
-        return (1,)
+        return (2,)
 
     #def perform(self, node, (a_data, a_indices, a_indptr, b), (out,)):
     #    return NotImplemented()
@@ -1032,18 +1034,20 @@ class MulSDCSR(gof.Op):
         if( PyArray_DESCR(%(_indptr)s)->type_num != NPY_INT32)
         {PyErr_SetString(PyExc_NotImplementedError, "D"); %(fail)s;}
 
-        if (!%(_zout)s)
+        if (!%(_zout)s ||
+            (PyArray_DIMS(%(_zout)s)[0] != PyArray_DIMS(%(_indices)s)[0]) ||
+            !(PyArray_ISCONTIGUOUS(%(_zout)s)))
         {
+            Py_XDECREF(%(_zout)s);
             %(_zout)s = (PyArrayObject*) PyArray_SimpleNew(1,
                     PyArray_DIMS(%(_indices)s), PyArray_DESCR(%(_b)s)->type_num);
-        }
-
-        if (PyArray_DIMS(%(_zout)s)[0] != PyArray_DIMS(%(_indices)s)[0])
+            if (!%(_zout)s)
             {
-            PyErr_SetString(PyExc_NotImplementedError,
-    "somehow _zout got the wrong size.. and I don't know how to resize it.");
+                PyErr_SetString(PyExc_MemoryError,
+                    "Could not allocate output memory.");
                 %(fail)s;
             }
+        }
 
         { //makes it compile even though labels jump over variable definitions.
             const npy_intp nnz = PyArray_DIMS(%(_indices)s)[0];
@@ -1302,7 +1306,7 @@ class StructuredAddSVCSR(gof.Op):
                                [tensor.tensor(b.dtype, (False,))])
 
     def c_code_cache_version(self):
-        return (1,)
+        return (2,)
 
     def c_code(self, node, name, inputs, outputs, sub):
         _data, _indices, _indptr, _b, = inputs
@@ -1336,18 +1340,20 @@ class StructuredAddSVCSR(gof.Op):
         if( PyArray_DESCR(%(_indptr)s)->type_num != NPY_INT32)
         {PyErr_SetString(PyExc_NotImplementedError, "D"); %(fail)s;}
 
-        if (!%(_zout)s)
+        if (!%(_zout)s
+            || (PyArray_DIMS(%(_zout)s)[0] != PyArray_DIMS(%(_indices)s)[0])
+            || !(PyArray_ISCONTIGUOUS(%(_zout)s)))
         {
+            Py_XDECREF(%(_zout)s);
             %(_zout)s = (PyArrayObject*) PyArray_SimpleNew(1,
                     PyArray_DIMS(%(_indices)s), PyArray_DESCR(%(_b)s)->type_num);
-        }
-
-        if (PyArray_DIMS(%(_zout)s)[0] != PyArray_DIMS(%(_indices)s)[0])
+            if (!%(_zout)s)
             {
-            PyErr_SetString(PyExc_NotImplementedError,
-     "somehow _zout got the wrong size.. and I don't know how to resize it.");
+                PyErr_SetString(PyExc_MemoryError,
+                    "Could not allocate output memory.");
                 %(fail)s;
             }
+        }
 
         { //makes it compile even though labels jump over variable definitions.
             const npy_intp nnz = PyArray_DIMS(%(_indices)s)[0];
@@ -1489,7 +1495,7 @@ class SamplingDotCSR(gof.Op):
         ])
 
     def c_code_cache_version(self):
-        return (1, )
+        return (2, )
 
     def c_support_code(self):
         return blas.blas_header_text()
@@ -1572,7 +1578,9 @@ PyErr_SetString(PyExc_NotImplementedError, "rank(y) != 2"); %(fail)s;}
         // Allocate output
         if (!%(z_data)s
             || (PyArray_DIMS(%(z_data)s)[0] != PyArray_DIMS(%(p_data)s)[0])
-            || (PyArray_DESCR(%(z_data)s)->type_num != %(typenum_zd)s)) {
+            || (PyArray_DESCR(%(z_data)s)->type_num != %(typenum_zd)s)
+            || !(PyArray_ISCONTIGUOUS(%(z_data)s)))
+         {
             {Py_XDECREF(%(z_data)s);}
             npy_intp dims[] = {0};
             dims[0] = PyArray_DIMS(%(p_data)s)[0];
@@ -1581,7 +1589,9 @@ PyErr_SetString(PyExc_NotImplementedError, "rank(y) != 2"); %(fail)s;}
         }
         if (!%(z_ind)s
             || (PyArray_DIMS(%(z_ind)s)[0] != PyArray_DIMS(%(p_ind)s)[0])
-            || (PyArray_DESCR(%(z_ind)s)->type_num != %(typenum_zi)s)) {
+            || (PyArray_DESCR(%(z_ind)s)->type_num != %(typenum_zi)s)
+            || !(PyArray_ISCONTIGUOUS(%(z_ind)s)))
+        {
             {Py_XDECREF(%(z_ind)s);}
             npy_intp dims[] = {0};
             dims[0] = PyArray_DIMS(%(p_ind)s)[0];
@@ -1590,7 +1600,9 @@ PyErr_SetString(PyExc_NotImplementedError, "rank(y) != 2"); %(fail)s;}
         }
         if (!%(z_ptr)s
             || (PyArray_DIMS(%(z_ptr)s)[0] != PyArray_DIMS(%(p_ptr)s)[0])
-            || (PyArray_DESCR(%(z_ptr)s)->type_num != %(typenum_zp)s)) {
+            || (PyArray_DESCR(%(z_ptr)s)->type_num != %(typenum_zp)s)
+            || !(PyArray_ISCONTIGUOUS(%(z_ptr)s)))
+        {
             {Py_XDECREF(%(z_ptr)s);}
             npy_intp dims[] = {0};
             dims[0] = PyArray_DIMS(%(p_ptr)s)[0];

theano/sparse/sandbox/sp2.py

-import theano
 import numpy
 import scipy.sparse
 
-from theano import gof, tensor, scalar, sparse
-from theano.tensor import blas
+from theano import gof, tensor
 
 from theano.sparse.basic import (
     as_sparse_variable, SparseType, add_s_s, neg,
@@ -18,14 +16,20 @@ from theano.sparse.basic import (
     HStack, hstack, VStack, vstack,
     AddSSData, add_s_s_data,
     MulSV, mul_s_v,
-    Multinomial, multinomial, Poisson, poisson,
-    Binomial, csr_fbinomial, csc_fbinomial, csr_dbinomial, csc_dbinomial,
     structured_monoid,
     structured_sigmoid, structured_exp, structured_log, structured_pow,
     structured_minimum, structured_maximum, structured_add,
     StructuredAddSV, structured_add_s_v,
     SamplingDot, sampling_dot)
 
+# Probability Ops are currently back in sandbox, because they do not respect
+# Theano's Op contract, as their behaviour is not reproducible: calling
+# the perform() method twice with the same argument will yield different
+# results.
+#from theano.sparse.basic import (
+#    Multinomial, multinomial, Poisson, poisson,
+#    Binomial, csr_fbinomial, csc_fbinomial, csr_dbinomial, csc_dbinomial)
+
 # Also for compatibility
 from theano.sparse.opt import (
     MulSDCSC, mul_s_d_csc, MulSDCSR, mul_s_d_csr,
@@ -35,6 +39,171 @@ from theano.sparse.opt import (
     local_mul_s_d, local_mul_s_v,
     local_structured_add_s_v, local_sampling_dot_csr)
 
+
 # Alias to maintain compatibility
 EliminateZeros = Remove0
 eliminate_zeros = remove0
+
+
+# Probability
+class Poisson(gof.op.Op):
+    """Return a sparse having random values from a Poisson density
+    with mean from the input.
+
+    WARNING: This Op is NOT deterministic, as calling it twice with the
+    same inputs will NOT give the same result. This is a violation of
+    Theano's contract for Ops
+
+    :param x: Sparse matrix.
+
+    :return: A sparse matrix of random integers of a Poisson density
+             with mean of `x` element wise.
+    """
+
+    def __eq__(self, other):
+        return (type(self) == type(other))
+
+    def __hash__(self):
+        return hash(type(self))
+
+    def make_node(self, x):
+        x = as_sparse_variable(x)
+        return gof.Apply(self, [x], [x.type()])
+
+    def perform(self, node, (x, ), (out, )):
+        assert _is_sparse(x)
+        out[0] = x.copy()
+        out[0].data = numpy.asarray(numpy.random.poisson(out[0].data),
+                                    dtype=x.dtype)
+        out[0].eliminate_zeros()
+
+    def grad(self, inputs, outputs_gradients):
+        return [None]
+
+    def infer_shape(self, node, ins_shapes):
+        return ins_shapes
+
+    def __str__(self):
+        return self.__class__.__name__
+poisson = Poisson()
+
+
+class Binomial(gof.op.Op):
+    """Return a sparse matrix having random values from a binomial
+    density having number of experiment `n` and probability of succes
+    `p`.
+
+    WARNING: This Op is NOT deterministic, as calling it twice with the
+    same inputs will NOT give the same result. This is a violation of
+    Theano's contract for Ops
+
+    :param n: Tensor scalar representing the number of experiment.
+    :param p: Tensor scalar representing the probability of success.
+    :param shape: Tensor vector for the output shape.
+
+    :return: A sparse matrix of integers representing the number
+             of success.
+    """
+
+    def __init__(self, format, dtype):
+        self.format = format
+        self.dtype = dtype
+
+    def __eq__(self, other):
+        return ((type(self) == type(other)) and
+                self.format == other.format and
+                self.dtype == other.dtype)
+
+    def __hash__(self):
+        return hash(type(self)) ^ hash(self.format) ^ hash(self.dtype)
+
+    def make_node(self, n, p, shape):
+        n = tensor.as_tensor_variable(n)
+        p = tensor.as_tensor_variable(p)
+        shape = tensor.as_tensor_variable(shape)
+        return gof.Apply(self, [n, p, shape], [SparseType(dtype=self.dtype,
+                                 format=self.format).make_variable()])
+
+    def perform(self, node, (n, p, shape, ), (out, )):
+        binomial = numpy.random.binomial(n, p, size=shape)
+        csx_matrix = getattr(scipy.sparse, self.format + '_matrix')
+        out[0] = csx_matrix(binomial, dtype=self.dtype)
+
+    def grad(self, (n, p, shape, ), (gz,)):
+        return None, None, None
+
+    def infer_shape(self, node, ins_shapes):
+        return [(node.inputs[2][0], node.inputs[2][1])]
+
+    def __str__(self):
+        return self.__class__.__name__
+
+csr_fbinomial = Binomial('csr', 'float32')
+csc_fbinomial = Binomial('csc', 'float32')
+csr_dbinomial = Binomial('csr', 'float64')
+csc_dbinomial = Binomial('csc', 'float64')
+
+
+class Multinomial(gof.op.Op):
+    """Return a sparse matrix having random values from a multinomial
+    density having number of experiment `n` and probability of succes
+    `p`.
+
+    WARNING: This Op is NOT deterministic, as calling it twice with the
+    same inputs will NOT give the same result. This is a violation of
+    Theano's contract for Ops
+
+    :param n: Tensor type vector or scalar representing the number of
+              experiment for each row. If `n` is a scalar, it will be
+              used for each row.
+    :param p: Sparse matrix of probability where each row is a probability
+              vector representing the probability of succes. N.B. Each row
+              must sum to one.
+
+    :return: A sparse matrix of random integers from a multinomial density
+             for each row.
+
+    :note: It will works only if `p` have csr format.
+    """
+
+    def __eq__(self, other):
+        return (type(self) == type(other))
+
+    def __hash__(self):
+        return hash(type(self))
+
+    def make_node(self, n, p):
+        n = tensor.as_tensor_variable(n)
+        p = as_sparse_variable(p)
+
+        return gof.Apply(self, [n, p], [p.type()])
+
+    def perform(self, node, (n, p), (out, )):
+        assert _is_sparse(p)
+
+        if p.format != 'csr':
+            raise NotImplemented()
+
+        out[0] = p.copy()
+
+        if n.ndim == 0:
+            for i in xrange(p.shape[0]):
+                k, l = p.indptr[i], p.indptr[i + 1]
+                out[0].data[k:l] = numpy.random.multinomial(n, p.data[k:l])
+        elif n.ndim == 1:
+            if n.shape[0] != p.shape[0]:
+                raise ValueError('The number of element of n must be '
+                                 'the same as the number of row of p.')
+            for i in xrange(p.shape[0]):
+                k, l = p.indptr[i], p.indptr[i + 1]
+                out[0].data[k:l] = numpy.random.multinomial(n[i], p.data[k:l])
+
+    def grad(self, inputs, outputs_gradients):
+        return [None, None]
+
+    def infer_shape(self, node, ins_shapes):
+        return [ins_shapes[1]]
+
+    def __str__(self):
+        return self.__class__.__name__
+multinomial = Multinomial()

theano/sparse/tests/test_basic.py

@@ -34,7 +34,6 @@ from theano.sparse import (
     Dot, Usmm, sp_ones_like, GetItemScalar,
     SparseFromDense,
     Cast, cast, HStack, VStack, AddSSData, add_s_s_data,
-    Poisson, poisson, Binomial, Multinomial, multinomial,
     structured_sigmoid, structured_exp, structured_log,
     structured_pow, structured_minimum, structured_maximum, structured_add,
     MulSV, mul_s_v, StructuredAddSV, structured_add_s_v,
@@ -42,6 +41,10 @@ from theano.sparse import (
     Diag, diag, SquareDiagonal, square_diagonal,
     EnsureSortedIndices, ensure_sorted_indices, clean)
 
+# Probability distributions are currently tested in test_sp2.py
+#from theano.sparse import (
+#    Poisson, poisson, Binomial, Multinomial, multinomial)
+
 from theano.sparse.opt import (StructuredDotCSC, UsmmCscDense, CSMGradC)
 
 from theano.tests import unittest_tools as utt
@@ -2141,122 +2144,6 @@ class AddSSDataTester(utt.InferShapeTester):
                                structured=True)
 
 
-class PoissonTester(utt.InferShapeTester):
-    x = {}
-    a = {}
-
-    for format in sparse.sparse_formats:
-        variable = getattr(theano.sparse, format + '_matrix')
-
-        rand = numpy.array(numpy.random.random_integers(3, size=(3, 4)) - 1,
-                           dtype=theano.config.floatX)
-
-        x[format] = variable()
-        a[format] = as_sparse_format(rand, format)
-
-    def setUp(self):
-        super(PoissonTester, self).setUp()
-        self.op_class = Poisson
-
-    def test_op(self):
-        for format in sparse.sparse_formats:
-            f = theano.function(
-                [self.x[format]],
-                poisson(self.x[format]))
-
-            tested = f(self.a[format])
-
-            assert tested.format == format
-            assert tested.dtype == self.a[format].dtype
-            assert numpy.allclose(numpy.floor(tested.data), tested.data)
-            assert tested.shape == self.a[format].shape
-
-    def test_infer_shape(self):
-        for format in sparse.sparse_formats:
-            self._compile_and_check([self.x[format]],
-                                    [poisson(self.x[format])],
-                                    [self.a[format]],
-                                    self.op_class)
-
-
-class BinomialTester(utt.InferShapeTester):
-    n = tensor.scalar()
-    p = tensor.scalar()
-    shape = tensor.lvector()
-    _n = 5
-    _p = .25
-    _shape = numpy.asarray([3, 5], dtype='int64')
-
-    inputs = [n, p, shape]
-    _inputs = [_n, _p, _shape]
-
-    def setUp(self):
-        super(BinomialTester, self).setUp()
-        self.op_class = Binomial
-
-    def test_op(self):
-        for sp_format in sparse.sparse_formats:
-            for o_type in sparse.float_dtypes:
-                f = theano.function(
-                    self.inputs,
-                    Binomial(sp_format, o_type)(*self.inputs))
-
-                tested = f(*self._inputs)
-
-                assert tested.shape == tuple(self._shape)
-                assert tested.format == sp_format
-                assert tested.dtype == o_type
-                assert numpy.allclose(numpy.floor(tested.todense()),
-                                   tested.todense())
-
-    def test_infer_shape(self):
-        for sp_format in sparse.sparse_formats:
-            for o_type in sparse.float_dtypes:
-                self._compile_and_check(
-                    self.inputs,
-                    [Binomial(sp_format, o_type)(*self.inputs)],
-                    self._inputs,
-                    self.op_class)
-
-
-class MultinomialTester(utt.InferShapeTester):
-    p = sparse.csr_matrix()
-    _p = sp.csr_matrix(numpy.asarray([[0.0, 0.5, 0.0, 0.5],
-                                      [0.1, 0.2, 0.3, 0.4],
-                                      [0.0, 1.0, 0.0, 0.0],
-                                      [0.3, 0.3, 0.0, 0.4]],
-                                     dtype=config.floatX))
-
-    def setUp(self):
-        super(MultinomialTester, self).setUp()
-        self.op_class = Multinomial
-
-    def test_op(self):
-        n = tensor.lscalar()
-        f = theano.function([self.p, n], multinomial(n, self.p))
-
-        _n = 5
-        tested = f(self._p, _n)
-        assert tested.shape == self._p.shape
-        assert numpy.allclose(numpy.floor(tested.todense()), tested.todense())
-        assert tested[2, 1] == _n
-
-        n = tensor.lvector()
-        f = theano.function([self.p, n], multinomial(n, self.p))
-
-        _n = numpy.asarray([1, 2, 3, 4], dtype='int64')
-        tested = f(self._p, _n)
-        assert tested.shape == self._p.shape
-        assert numpy.allclose(numpy.floor(tested.todense()), tested.todense())
-        assert tested[2, 1] == _n[2]
-
-    def test_infer_shape(self):
-        self._compile_and_check([self.p],
-                                [multinomial(5, self.p)],
-                                [self._p],
-                                self.op_class)
-
-
 def elemwise_checker(op, expected_f, gap=None, test_dtypes=None,
                      grad_test=True, name=None):
     """Return the appropriate test class for the elemwise on sparse.

theano/sparse/tests/test_sp2.py

-import time
 import unittest
 
 from nose.plugins.skip import SkipTest
-import numpy as np
+import numpy
 try:
     import scipy.sparse as sp
-    import scipy.sparse
 except ImportError:
     pass  # The variable enable_sparse will be used to disable the test file.
 
 import theano
+from theano import config
 from theano import tensor
 from theano import sparse
 
 if not theano.sparse.enable_sparse:
     raise SkipTest('Optional package sparse disabled')
 
-from theano.sparse.sandbox import sp2 as S2
+from theano.sparse.sandbox.sp2 import (
+    Poisson, poisson, Binomial, Multinomial, multinomial)
 
 from theano.tests import unittest_tools as utt
-from theano.sparse.basic import verify_grad_sparse
+from theano.sparse.tests.test_basic import as_sparse_format
+
+
+class PoissonTester(utt.InferShapeTester):
+    x = {}
+    a = {}
+
+    for format in sparse.sparse_formats:
+        variable = getattr(theano.sparse, format + '_matrix')
+
+        rand = numpy.array(numpy.random.random_integers(3, size=(3, 4)) - 1,
+                           dtype=theano.config.floatX)
+
+        x[format] = variable()
+        a[format] = as_sparse_format(rand, format)
+
+    def setUp(self):
+        super(PoissonTester, self).setUp()
+        self.op_class = Poisson
+
+    def test_op(self):
+        for format in sparse.sparse_formats:
+            f = theano.function(
+                [self.x[format]],
+                poisson(self.x[format]))
+
+            tested = f(self.a[format])
+
+            assert tested.format == format
+            assert tested.dtype == self.a[format].dtype
+            assert numpy.allclose(numpy.floor(tested.data), tested.data)
+            assert tested.shape == self.a[format].shape
+
+    def test_infer_shape(self):
+        for format in sparse.sparse_formats:
+            self._compile_and_check([self.x[format]],
+                                    [poisson(self.x[format])],
+                                    [self.a[format]],
+                                    self.op_class)
+
+
+class BinomialTester(utt.InferShapeTester):
+    n = tensor.scalar()
+    p = tensor.scalar()
+    shape = tensor.lvector()
+    _n = 5
+    _p = .25
+    _shape = numpy.asarray([3, 5], dtype='int64')
+
+    inputs = [n, p, shape]
+    _inputs = [_n, _p, _shape]
+
+    def setUp(self):
+        super(BinomialTester, self).setUp()
+        self.op_class = Binomial
+
+    def test_op(self):
+        for sp_format in sparse.sparse_formats:
+            for o_type in sparse.float_dtypes:
+                f = theano.function(
+                    self.inputs,
+                    Binomial(sp_format, o_type)(*self.inputs))
+
+                tested = f(*self._inputs)
+
+                assert tested.shape == tuple(self._shape)
+                assert tested.format == sp_format
+                assert tested.dtype == o_type
+                assert numpy.allclose(numpy.floor(tested.todense()),
+                                   tested.todense())
+
+    def test_infer_shape(self):
+        for sp_format in sparse.sparse_formats:
+            for o_type in sparse.float_dtypes:
+                self._compile_and_check(
+                    self.inputs,
+                    [Binomial(sp_format, o_type)(*self.inputs)],
+                    self._inputs,
+                    self.op_class)
+
+
+class MultinomialTester(utt.InferShapeTester):
+    p = sparse.csr_matrix()
+    _p = sp.csr_matrix(numpy.asarray([[0.0, 0.5, 0.0, 0.5],
+                                      [0.1, 0.2, 0.3, 0.4],
+                                      [0.0, 1.0, 0.0, 0.0],
+                                      [0.3, 0.3, 0.0, 0.4]],
+                                     dtype=config.floatX))
+
+    def setUp(self):
+        super(MultinomialTester, self).setUp()
+        self.op_class = Multinomial
+
+    def test_op(self):
+        n = tensor.lscalar()
+        f = theano.function([self.p, n], multinomial(n, self.p))
+
+        _n = 5
+        tested = f(self._p, _n)
+        assert tested.shape == self._p.shape
+        assert numpy.allclose(numpy.floor(tested.todense()), tested.todense())
+        assert tested[2, 1] == _n
+
+        n = tensor.lvector()
+        f = theano.function([self.p, n], multinomial(n, self.p))
+
+        _n = numpy.asarray([1, 2, 3, 4], dtype='int64')
+        tested = f(self._p, _n)
+        assert tested.shape == self._p.shape
+        assert numpy.allclose(numpy.floor(tested.todense()), tested.todense())
+        assert tested[2, 1] == _n[2]
+
+    def test_infer_shape(self):
+        self._compile_and_check([self.p],
+                                [multinomial(5, self.p)],
+                                [self._p],
+                                self.op_class)
 
 
 if __name__ == '__main__':

theano/tensor/io.py

 import numpy
-import theano
 from theano import gof
-from theano.gof import Apply, Constant, Generic, Op, Type, Variable
+from theano.gof import Constant, Generic, Op
 from basic import tensor
 ##########################
 # Disk Access
 ##########################
 
+
 class LoadFromDisk(Op):
     """
     An operation to load an array from disk
@@ -19,6 +19,9 @@ class LoadFromDisk(Op):
     def __init__(self, dtype, broadcastable, mmap_mode=None):
         self.dtype = numpy.dtype(dtype)  # turn "float64" into numpy.float64
         self.broadcastable = broadcastable
+        if mmap_mode not in (None, 'c'):
+            raise ValueError("The only supported values for mmap_mode "
+                    "are None and 'c', got %s" % mmap_mode)
         self.mmap_mode = mmap_mode
         self._info = (dtype, broadcastable, mmap_mode)
 
@@ -37,19 +40,33 @@ class LoadFromDisk(Op):
     def perform(self, node, inp, out):
         path = inp[0]
         if (path.split('.')[-1] == 'npz'):
-            raise ValueError("Expected a .npy file, got %s instead"%path)
+            raise ValueError("Expected a .npy file, got %s instead" % path)
         result = numpy.load(path, mmap_mode=self.mmap_mode)
         if result.dtype != self.dtype:
-            raise TypeError("Expected an array of type %s, got %s instead"%
+            raise TypeError("Expected an array of type %s, got %s instead" %
                     (self.dtype, result.dtype))
+        print 'result:', result, type(result)
         out[0][0] = result
 
     def __str__(self):
-        return "Load{dtype:%s, broadcastable:%s, mmep:%s}"%self._info
+        return "Load{dtype:%s, broadcastable:%s, mmep:%s}" % self._info
+
 
 def load(path, dtype, broadcastable, mmap_mode=None):
     """
-    Load an array from an .npy file
+    Load an array from an .npy file.
+
+    :param path: A Generic symbolic variable, that will contain a string
+    :param dtype: The data type of the array to be read.
+    :param broadcastable: The broadcastable pattern of the loaded array,
+      for instance, (False,) for a vector, (False, True) for a column,
+      (False, False) for a matrix.
+    :param mmap_mode: How the file will be loaded. None means that the
+      data will be copied into an array in memory, 'c' means that the file
+      will be mapped into virtual memory, so only the parts that are
+      needed will be actually read from disk and put into memory.
+      Other modes supported by numpy.load ('r', 'r+', 'w+') cannot
+      be supported by Theano.
 
     >>> from theano import *
     >>> path = Variable(Generic())
@@ -61,4 +78,3 @@ def load(path, dtype, broadcastable, mmap_mode=None):
     """
 
     return LoadFromDisk(dtype, broadcastable, mmap_mode)(path)
-

theano/tensor/tests/test_io.py

@@ -6,28 +6,49 @@ import os
 
 
 class T_load_tensor(unittest.TestCase):
-    def test0(self):
-        data = numpy.arange(5, dtype=numpy.int32)
-        filename = os.path.join(
+    def setUp(self):
+        self.data = numpy.arange(5, dtype=numpy.int32)
+        self.filename = os.path.join(
             theano.config.base_compiledir,
             "_test.npy")
-        numpy.save(filename, data)
+        numpy.save(self.filename, self.data)
+
+    def test0(self):
         path = Variable(Generic())
+        # Not specifying mmap_mode defaults to None, and the data is
+        # copied into main memory
         x = tensor.load(path, 'int32', (False,))
-        y = x*2
+        y = x * 2
+        fn = function([path], y)
+        assert (fn(self.filename) == (self.data * 2)).all()
+
+    def test_invalid_modes(self):
+        # Modes 'r+', 'r', and 'w+' cannot work with Theano, becausei
+        # the output array may be modified inplace, and that should not
+        # modify the original file.
+        path = Variable(Generic())
+        for mmap_mode in ('r+', 'r', 'w+', 'toto'):
+            self.assertRaises(ValueError,
+                    tensor.load, path, 'int32', (False,), mmap_mode)
+
+    def test1(self):
+        path = Variable(Generic())
+        # 'c' means "copy-on-write", which allow the array to be overwritten
+        # by an inplace Op in the graph, without modifying the underlying
+        # file.
+        x = tensor.load(path, 'int32', (False,), 'c')
+        # x ** 2 has been chosen because it will work inplace.
+        y = (x ** 2).sum()
         fn = function([path], y)
-        assert (fn(filename) == data*2).all()
+        # Call fn() twice, to check that inplace ops do not cause trouble
+        assert (fn(self.filename) == (self.data ** 2).sum()).all()
+        assert (fn(self.filename) == (self.data ** 2).sum()).all()
 
     def test_memmap(self):
-        data = numpy.arange(5, dtype=numpy.int32)
-        filename = os.path.join(
-            theano.config.base_compiledir,
-            "_test.npy")
-        numpy.save(filename, data)
         path = Variable(Generic())
-        x = tensor.load(path, 'int32', (False,), mmap_mode='r+')
+        x = tensor.load(path, 'int32', (False,), mmap_mode='c')
         fn = function([path], x)
-        assert type(fn(filename)) == numpy.core.memmap
+        assert type(fn(self.filename)) == numpy.core.memmap
 
     def tearDown(self):
         os.remove(os.path.join(