merge

doc/library/scan.txt

@@ -10,15 +10,16 @@ Guide
 =====
 
 The scan functions provides the basic functionality needed to do loops
-in Theano. Scan comes with many whistles and bells, that can be easily
-introduced through a few examples :
+in Theano. Scan comes with many whistles and bells, which we will introduce
+by way of examples.
 
-Basic functionality :  Computing :math:`A^k`
---------------------------------------------
+
+Simple loop with accumulation:  Computing :math:`A^k`
+-----------------------------------------------------
 
 Assume that, given *k* you want to get ``A**k`` using a loop.
 More precisely, if *A* is a tensor you want to compute
-``A**k`` elemwise. The python/numpy code would loop like
+``A**k`` elemwise. The python/numpy code might look like:
 
 .. code-block:: python
 
@@ -26,42 +27,176 @@ More precisely, if *A* is a tensor you want to compute
   for i in xrange(k):
     result = result * A
 
-The equivalent Theano code would be
+There are three thing here that we need to handle: the initial value
+assigned to ``result``, the accumulation of results in ``result``, and
+the unchanging variable ``A``. Unchanging variables are passed to scan as
+``non_sequences``. Initialization occurs in ``outputs_info``, and the accumulation
+happens automatically.
+
+The equivalent Theano code would be:
 
 .. code-block:: python
 
   # Symbolic description of the result
-  result,updates = theano.scan(fn = lambda x_tm1,A: x_tm1*A,\
-                       outputs_info = T.ones_like(A),\
-                       non_sequences  = A, \
-                       n_steps        = k)
+  result, updates = theano.scan(fn=lambda prior_result, A: prior_result * A,
+                                outputs_info=T.ones_like(A),
+                                non_sequences=A,
+                                n_steps=k)
+
+  # We only care about A**k, but scan has provided us with A**1 through A**k.
+  # Discard the values that we don't care about. Scan is smart enough not to
+  # notice this and not waste memory saving them.
+  final_result = result[-1]
 
   # compiled function that returns A**k
-  f = theano.function([A,k], result[-1], updates = updates)
+  power = theano.function(inputs=[A,k], outputs=final_result, updates=updates)
 
 Let us go through the example line by line. What we did is first to
-construct a function (using a lambda expression) that given `x_tm1` and
-`A` returns `x_tm1*A`. Given the order of the parameters, `x_tm1`
-is the value of our output at time step ``t-1``. Therefore
-``x_t`` (value of output at time `t`) is `A` times value of output
-at `t-1`.
-Next we initialize the output as a tensor with same
-shape as A filled with ones. We give A to scan as a non sequence parameter  and
-specify the number of steps k to iterate over our lambda expression.
-
-Scan will return a tuple, containing our result (``result``) and a
-dictionary of updates ( empty in this case). Note that the result
+construct a function (using a lambda expression) that given ``prior_result`` and
+``A`` returns ``prior_result * A``. The order of parameters is fixed by scan:
+the output of the prior call to ``fn`` (or the initial value, initially)
+is the first parameter, followed by all non-sequences.
+
+Next we initialize the output as a tensor with same shape and dtype as ``A``,
+filled with ones. We give ``A`` to scan as a non sequence parameter and
+specify the number of steps ``k`` to iterate over our lambda expression.
+
+Scan return a tuples, containing our result (``result``) and a
+dictionary of updates (empty in this case). Note that the result
 is not a matrix, but a 3D tensor containing the value of ``A**k`` for
-each step. We want the last value ( after k steps ) so we compile
+each step. We want the last value (after ``k`` steps) so we compile
 a function to return just that. Note that there is an optimization, that
 at compile time will detect that you are using just the last value of the
 result and ensure that scan does not store all the intermediate values
-that are used. So do not worry if A and k are large.
+that are used. So do not worry if ``A`` and ``k`` are large.
+
+
+Iterating over the first dimension of a tensor: Calculating a polynomial
+------------------------------------------------------------------------
+In addition to looping a fixed number of times, scan can iterate over
+the leading dimension of tensors (similar to Python's ``for x in a_list``).
+
+The tensor(s) to be looped over should be provided to scan using the
+``sequence`` keyword argument.
+
+Here's an example that builds a symbolic calculation of a polynomial
+from a list of its coefficients:
+
+.. code-block:: python
+
+    coefficients = theano.tensor.vector("coefficients")
+    x = T.scalar("x")
+
+    max_coefficients_supported = 10000
+
+    # Generate the components of the polynomial
+    components, updates = theano.scan(fn=lambda coefficient, power, free_variable: coefficient * (free_variable ** power),
+                                      outputs_info=None,
+                                      sequences=[coefficients, theano.tensor.arange(max_coefficients_supported)],
+                                      non_sequences=x)
+    # Sum them up
+    polynomial = components.sum()
+
+    # Compile a function
+    calculate_polynomial = theano.function(inputs=[coefficients, x], outputs=polynomial)
+
+    # Test
+    test_coefficients = numpy.asarray([1, 0, 2], dtype=numpy.float32)
+    test_value = 3
+    print calculate_polynomial(test_coefficients, test_value)
+    print 1.0 * (3 ** 0) + 0.0 * (3 ** 1) + 2.0 * (3 ** 2) 
+
+There are a few things to note here.
+
+First, we calculate the polynomial by first generating each of the coefficients, and
+then summing them at the end. (We could also have accumulated them along the way, and then
+taken the last one, which would have been more memory-efficient, but this is an example.)
+
+Second, there is no accumulation of results, we can set ``outputs_info`` to ``None``. This indicates
+to scan that it doesn't need to pass the prior result to ``fn``.
+
+The general order of function parameters to ``fn`` is::
+
+    sequences (if any), prior result(s) (if needed), non-sequences (if any)
+
+Third, there's a handy trick used to simulate python's ``enumerate``: simply include
+``theano.tensor.arange`` to the sequences.
+
+Fourth, given multiple sequences of uneven lengths, scan will truncate to the shortest of them.
+This makes it safe to pass a very long arange, which we need to do for generality, since
+arange must have its length specified at creation time.
+
+
+Simple accumulation into a scalar, ditching lamba
+-------------------------------------------------
+
+This should be fairly self-explanatory.
+
+.. code-block:: python
+
+    up_to = T.iscalar("up_to")
+
+    # define a named function, rather than using lambda
+    def accumulate_by_adding(arange_val, sum_to_date):
+        return sum_to_date + arange_val
+
+    scan_result, scan_updates = theano.scan(fn=accumulate_by_adding,
+                                            outputs_info=T.as_tensor_variable(0),
+                                            sequences=T.arange(up_to))
+    triangular_sequence = theano.function(inputs=[up_to], outputs=scan_result)
+
+    # test
+    some_num = 15
+    print triangular_sequence(some_num)
+    print [n * (n + 1) // 2 for n in xrange(some_num)]
+
+
+
+Another simple example
+----------------------
+
+Unlike some of the prior examples, this one is hard to reproduce except by using scan.
+
+This takes a sequence of array indices, and values to place there,
+and a "model" output array (whose shape and dtype will be mimicked),
+and produces a sequence of arrays with the shape and dtype of the model,
+with all values set to zero except at the provided array indices.
+
+.. code-block:: python
+
+    location = T.imatrix("location")
+    values = T.vector("values")
+    output_model = T.matrix("output_model")
+
+    def set_value_at_position(a_location, a_value, output_model):
+        zeros = T.zeros_like(output_model)
+        zeros_subtensor = zeros[a_location[0], a_location[1]]
+        return T.set_subtensor(zeros_subtensor, a_value)
+
+    result, updates = theano.scan(fn=set_value_at_position,
+                                  outputs_info=None,
+                                  sequences=[location, values],
+                                  non_sequences=output_model)
+
+    assign_values_at_positions = theano.function(inputs=[location, values, output_model], outputs=result)
+
+    # test
+    test_locations = numpy.asarray([[1, 1], [2, 3]], dtype=numpy.int32)
+    test_values = numpy.asarray([42, 50], dtype=numpy.float32)
+    test_output_model = numpy.zeros((5, 5), dtype=numpy.float32)
+    print assign_values_at_positions(test_locations, test_values, test_output_model)
+
+This demonstrates that you can introduce new Theano variables into a scan function.
+
 
 Multiple outputs, several taps values - Recurrent Neural Network with Scan
 --------------------------------------------------------------------------
 
-A more practical task would be to implement a RNN using scan. Assume
+The examples above showed simple uses of scan. However, scan also supports
+referring not only to the prior result and the current sequence value, but
+also looking back more than one step.
+
+This is needed, for example, to implement a RNN using scan. Assume
 that our RNN is defined as follows :
 
 .. math::

theano/gof/python25.py

@@ -69,3 +69,27 @@ else:
      partial = functools.partial
      defaultdict = collections.defaultdict
 __all__ = ['all', 'any']
+
+if sys.version_info[:2] < (2,6):
+    # Borrowed from Python docs
+    def combinations(iterable, r):
+        # combinations('ABCD', 2) --> AB AC AD BC BD CD
+        # combinations(range(4), 3) --> 012 013 023 123
+        pool = tuple(iterable)
+        n = len(pool)
+        if r > n:
+            return
+        indices = range(r)
+        yield tuple(pool[i] for i in indices)
+        while True:
+            for i in reversed(range(r)):
+                if indices[i] != i + n - r:
+                    break
+            else:
+                return
+            indices[i] += 1
+            for j in range(i+1, r):
+                indices[j] = indices[j-1] + 1
+            yield tuple(pool[i] for i in indices)
+else:
+    from itertools import combinations

theano/sandbox/cuda/cuda_ndarray.cu

@@ -2982,18 +2982,20 @@ CudaNdarray_dimshuffle(CudaNdarray * self, unsigned int len, const int * pattern
         }
 	else if(dims_taken[pattern[i]])
 	{
-	  PyErr_SetString(PyExc_ValueError, "Cudandarray_dimshuffle: The same input dimension may not appear twice in the list of output dimensions");
+                PyErr_Format(PyExc_ValueError, "Cudandarray_dimshuffle: invalid pattern for Cudandarray_dimshuffle. You used the dimensions %d multiple time",
+			     pattern[i]);
 	  free(newdims);
 	  return -1;
 	}
-        else
-        {
-            if ((dims_taken[pattern[i]]) || (pattern[i]>= self->nd))
-            {
-                PyErr_SetString(PyExc_ValueError, "Cudandarray_dimshuffle: invalid pattern for Cudandarray_dimshuffle");
-                free(newdims);
-                return -1;
-            }
+        else if (pattern[i]>= self->nd)
+	{
+	    PyErr_Format(PyExc_ValueError, "Cudandarray_dimshuffle: invalid pattern for Cudandarray_dimshuffle. You asked for a dimensions that don't exist %d for a %d dims CudaNdarray",
+			 pattern[i], self->nd);
+	    free(newdims);
+	    return -1;
+	}
+	else
+	{
             newdims[i] = CudaNdarray_HOST_DIMS(self)[pattern[i]];
             newstrides[i] = CudaNdarray_HOST_STRIDES(self)[pattern[i]];
             dims_taken[pattern[i]] = 1;

theano/scalar/basic.py

@@ -1103,19 +1103,9 @@ class Clip(ScalarOp):
             return gx, None, None
         else:
             return None, None, None
-clip = Clip(upcast_out, name = 'clip')
-
-class First(BinaryScalarOp):
-    def impl(self, x, y):
-        return x
-    def c_code(self, node, name, (x, y), (z, ), sub):
-        return "%(z)s = %(x)s;" % locals()
-    def grad(self, (x, y), (gz, )):
-        if x.type in continuous_types:
-            return gz, None
-        else:
-            return None,None
-first = First(transfer_type(0), name = 'first')
+# Don't allow complex even if numpy do
+# As there is no mathematical reason for this function on complex
+clip = Clip(upcast_out_no_complex, name = 'clip')
 
 class Second(BinaryScalarOp):
     def impl(self, x, y):

theano/tensor/opt.py

@@ -1095,20 +1095,36 @@ def local_useless_subtensor(node):
     """
     Remove Subtensor if it take the full input
     """
-    shape_of = node.env.shape_feature.shape_of
     if isinstance(node.op, T.Subtensor):
+        shape_of = node.env.shape_feature.shape_of
+        node_input_idx = 1
         for pos, idx in enumerate(node.op.idx_list):
+            length_pos_data = sys.maxint
+            length_pos_shape_i = None
             try:
-                length_pos = shape_i(node.inputs[0])[pos]
-            except:
+                length_pos = shape_of[node.inputs[0]][pos]
+                if isinstance(length_pos, theano.tensor.basic.TensorConstant):
+                    length_pos_data = length_pos.data
+                else:
+                    length_pos_shape_i = node.inputs[node_input_idx].owner.inputs[0]
+
+            except Exception, e:
                 length_pos = None
             if ( isinstance(idx,slice) and
                 idx.start in [0,None] and
                 idx.step in [1,None] and
-                idx.stop in [sys.maxint, None, length_pos]):
+                (idx.stop in [sys.maxint, None, length_pos_data] or
+                 (isinstance(idx.stop, int) and idx.stop>=length_pos_data) or
+                 (isinstance(idx.stop, theano.scalar.Scalar) and
+                  length_pos==length_pos_shape_i)
+                 )):
                 pass
             else:
                 return False
+            if isinstance(idx, slice):
+                node_input_idx += sum([isinstance(idx.start, theano.scalar.Scalar),
+                                       isinstance(idx.stop, theano.scalar.Scalar),
+                                       isinstance(idx.step, theano.scalar.Scalar)])
         return [node.inputs[0]]
 
 

theano/tensor/tests/test_opt.py

@@ -1147,12 +1147,76 @@ def test_log_add():
 
 def test_local_useless_subtensor():
     x = TT.matrix('x')
-    f = function([x], TT.exp(x)[0:], mode=mode_opt)
 
-    prog=f.maker.env.toposort()
-    assert prog[0].op == TT.exp
-    assert len(prog)==1
-    f([[0,1],[2,3]]) # let debugmode test something
+    # Test default
+    for dims in [(slice(0,None),),
+                 (slice(0,None),slice(0,None)),
+                 ]:
+        f = function([x], TT.exp(x).__getitem__(dims), mode=mode_opt)
+        #theano.printing.debugprint(f)
+        prog=f.maker.env.toposort()
+        assert prog[0].op == TT.exp
+        assert len(prog)==1
+        f([[0,1,2],[3,4,5]]) # let debugmode test something
+
+    x_c = specify_shape(x, (2,3))
+    # Test constant
+    for dims, res in [((slice(0,2),), True),
+                 ((slice(0,2),slice(0,None)), True),
+                 ((slice(0,2),slice(0,3)), True),
+                 ((slice(0,None),slice(0,3)), True),
+                 ((slice(0,3),slice(0,13)), True),
+                 ((slice(0,3),slice(0,2)), False),
+                 ((slice(0,1),slice(0,None)), False),
+                 ((slice(0,1),1), False),
+                 ]:
+        f = function([x], TT.exp(x_c).__getitem__(dims), mode=mode_opt)
+        #theano.printing.debugprint(f)
+        prog=f.maker.env.toposort()
+        if res:
+            assert isinstance(prog[0].op, theano.tensor.basic.SpecifyShape), dims
+            assert prog[1].op == TT.exp, dims
+            assert len(prog)==2, dims
+        else:
+            assert any([isinstance(node.op, Subtensor) for node in prog])
+        f([[0,1,2],[3,4,5]]) # let debugmode test something
+
+    # Test Variable
+    for idx, (dims, res) in enumerate([
+            ((slice(0,x.shape[0]),), True),
+            ((slice(0,x.shape[1]),), False),
+            ((slice(0,x.shape[0]),slice(0,x.shape[1]),), True),
+            ((slice(0,x.shape[0]),slice(0,x.shape[0]),), False),
+            ((slice(0,x.shape[1]),slice(0,x.shape[0]),), False),
+            ((slice(0,x.shape[1]),slice(0,x.shape[1]),), False),
+            ((slice(0,x.shape[1]),2), False),
+            ((slice(0,x.shape[1]),slice(x.shape[0]-x.shape[0],x.shape[1]),), False),
+            ]):
+        f = function([x], TT.exp(x).__getitem__(dims), mode=mode_opt)
+        #theano.printing.debugprint(f)
+        prog=f.maker.env.toposort()
+        if res:
+            assert prog[0].op == TT.exp, dims
+            assert len(prog)==1, dims
+        else:
+            assert any([isinstance(node.op, Subtensor) for node in prog])
+        f([[0,1,2],[3,4,5]]) # let debugmode test something
+    # Test mix Variable and Constant
+    # Currently not supported
+    for idx, (dims, res) in enumerate([
+            ((slice(0,x.shape[0]),slice(0,3)), False),
+            ((slice(0,3),slice(0,x.shape[1])), False),
+            ]):
+        f = function([x], TT.exp(x_c).__getitem__(dims), mode=mode_opt)
+        #theano.printing.debugprint(f)
+        prog=f.maker.env.toposort()
+        if res:
+            assert prog[0].op == TT.exp, dims
+            assert len(prog)==1, dims
+        else:
+            assert any([isinstance(node.op, Subtensor) for node in prog])
+        f([[0,1,2],[3,4,5]]) # let debugmode test something
+
 
 class test_local_subtensor_lift(unittest.TestCase):
 
@@ -2320,7 +2384,7 @@ class T_local_erfc(unittest.TestCase):
 
 class test_local_remove_switch_const_cond(unittest.TestCase):
     def setUp(self):
-        self.mode = theano.compile.get_default_mode().excluding('constant_folding')
+        self.mode = mode_opt.excluding('constant_folding')
 
     def test_const0(self):