Theano: how to efficiently undo/reverse max-pooling
Question: I'm using Theano 0.7 to create a [convolutional neural
net](http://deeplearning.net/tutorial/lenet.html) which uses **[max-
pooling](http://deeplearning.net/tutorial/lenet.html#maxpooling)** (i.e.
shrinking a matrix down by keeping only the local maxima).

In order to "undo" or "reverse" the max-pooling step, one method is to store
the locations of the maxima as auxiliary data, then simply recreate the un-
pooled data by making a big array of zeros and using those auxiliary locations
to place the maxima in their appropriate locations.

Here's how I'm currently doing it:

    
    
    import numpy as np
    import theano
    import theano.tensor as T
    
    minibatchsize = 2
    numfilters = 3
    numsamples = 4
    upsampfactor = 5
    
    # HERE is the function that I hope could be improved
    def upsamplecode(encoded, auxpos):
        shp = encoded.shape
        upsampled = T.zeros((shp[0], shp[1], shp[2] * upsampfactor))
        for whichitem in range(minibatchsize):
            for whichfilt in range(numfilters):
                upsampled = T.set_subtensor(upsampled[whichitem, whichfilt, auxpos[whichitem, whichfilt, :]], encoded[whichitem, whichfilt, :])
        return upsampled
    
    
    totalitems = minibatchsize * numfilters * numsamples
    
    code = theano.shared(np.arange(totalitems).reshape((minibatchsize, numfilters, numsamples)))
    
    auxpos = np.arange(totalitems).reshape((minibatchsize, numfilters, numsamples)) % upsampfactor  # arbitrary positions within a bin
    auxpos += (np.arange(4) * 5).reshape((1,1,-1)) # shifted to the actual temporal bin location
    auxpos = theano.shared(auxpos.astype(np.int))
    
    print "code:"
    print code.get_value()
    print "locations:"
    print auxpos.get_value()
    get_upsampled = theano.function([], upsamplecode(code, auxpos))
    print "the un-pooled data:"
    print get_upsampled()
    

> _(By the way, in this case I have a 3D tensor, and it's only the third axis
> that gets max-pooled. People who work with image data might expect to see
> two dimensions getting max-pooled.)_

The output is:

    
    
    code:
    [[[ 0  1  2  3]
      [ 4  5  6  7]
      [ 8  9 10 11]]
    
     [[12 13 14 15]
      [16 17 18 19]
      [20 21 22 23]]]
    locations:
    [[[ 0  6 12 18]
      [ 4  5 11 17]
      [ 3  9 10 16]]
    
     [[ 2  8 14 15]
      [ 1  7 13 19]
      [ 0  6 12 18]]]
    the un-pooled data:
    [[[  0.   0.   0.   0.   0.   0.   1.   0.   0.   0.   0.   0.   2.   0.
         0.   0.   0.   0.   3.   0.]
      [  0.   0.   0.   0.   4.   5.   0.   0.   0.   0.   0.   6.   0.   0.
         0.   0.   0.   7.   0.   0.]
      [  0.   0.   0.   8.   0.   0.   0.   0.   0.   9.  10.   0.   0.   0.
         0.   0.  11.   0.   0.   0.]]
    
     [[  0.   0.  12.   0.   0.   0.   0.   0.  13.   0.   0.   0.   0.   0.
        14.  15.   0.   0.   0.   0.]
      [  0.  16.   0.   0.   0.   0.   0.  17.   0.   0.   0.   0.   0.  18.
         0.   0.   0.   0.   0.  19.]
      [ 20.   0.   0.   0.   0.   0.  21.   0.   0.   0.   0.   0.  22.   0.
         0.   0.   0.   0.  23.   0.]]]
    

This method **works** but it's a **bottleneck** , taking most of my computer's
time (I think the set_subtensor calls might imply cpu<->gpu data copying). So:
can this be implemented more efficiently?

I suspect there's a way to express this as a single `set_subtensor()` call
which may be faster, but I don't see how to get the tensor indexing to
broadcast properly.

* * *

**UPDATE:** I thought of a way of doing it in one call, by working on the
flattened tensors:

    
    
    def upsamplecode2(encoded, auxpos):
        shp = encoded.shape
        upsampled = T.zeros((shp[0], shp[1], shp[2] * upsampfactor))
    
        add_to_flattened_indices = theano.shared(np.array([ [[(y + z * numfilters) * numsamples * upsampfactor for x in range(numsamples)] for y in range(numfilters)] for z in range(minibatchsize)], dtype=theano.config.floatX).flatten(), name="add_to_flattened_indices")
    
        upsampled = T.set_subtensor(upsampled.flatten()[T.cast(auxpos.flatten() + add_to_flattened_indices, 'int32')], encoded.flatten()).reshape(upsampled.shape)
    
        return upsampled
    
    
    get_upsampled2 = theano.function([], upsamplecode2(code, auxpos))
    print "the un-pooled data v2:"
    ups2 = get_upsampled2()
    print ups2
    

However, this is still not good efficiency-wise because when I run this (added
on to the end of the above script) I find out that the Cuda libraries can't
currently do the integer index manipulation efficiently:

    
    
    ERROR (theano.gof.opt): Optimization failure due to: local_gpu_advanced_incsubtensor1
    ERROR (theano.gof.opt): TRACEBACK:
    ERROR (theano.gof.opt): Traceback (most recent call last):
      File "/usr/local/lib/python2.7/dist-packages/theano/gof/opt.py", line 1493, in process_node
        replacements = lopt.transform(node)
      File "/usr/local/lib/python2.7/dist-packages/theano/sandbox/cuda/opt.py", line 952, in local_gpu_advanced_incsubtensor1
        gpu_y = gpu_from_host(y)
      File "/usr/local/lib/python2.7/dist-packages/theano/gof/op.py", line 507, in __call__
        node = self.make_node(*inputs, **kwargs)
      File "/usr/local/lib/python2.7/dist-packages/theano/sandbox/cuda/basic_ops.py", line 133, in make_node
        dtype=x.dtype)()])
      File "/usr/local/lib/python2.7/dist-packages/theano/sandbox/cuda/type.py", line 69, in __init__
        (self.__class__.__name__, dtype, name))
    TypeError: CudaNdarrayType only supports dtype float32 for now. Tried using dtype int64 for variable None
    


Answer: I don't know whether this is faster, but it may be a little more concise. See
if it is useful for your case.

    
    
    import numpy as np
    import theano
    import theano.tensor as T
    
    minibatchsize = 2
    numfilters = 3
    numsamples = 4
    upsampfactor = 5
    
    totalitems = minibatchsize * numfilters * numsamples
    
    code = np.arange(totalitems).reshape((minibatchsize, numfilters, numsamples))
    
    auxpos = np.arange(totalitems).reshape((minibatchsize, numfilters, numsamples)) % upsampfactor 
    auxpos += (np.arange(4) * 5).reshape((1,1,-1))
    
    # first in numpy
    shp = code.shape
    upsampled_np = np.zeros((shp[0], shp[1], shp[2] * upsampfactor))
    upsampled_np[np.arange(shp[0]).reshape(-1, 1, 1), np.arange(shp[1]).reshape(1, -1, 1), auxpos] = code
    
    print "numpy output:"
    print upsampled_np
    
    # now the same idea in theano
    encoded = T.tensor3()
    positions = T.tensor3(dtype='int64')
    shp = encoded.shape
    upsampled = T.zeros((shp[0], shp[1], shp[2] * upsampfactor))
    upsampled = T.set_subtensor(upsampled[T.arange(shp[0]).reshape((-1, 1, 1)), T.arange(shp[1]).reshape((1, -1, 1)), positions], encoded)
    
    print "theano output:"
    print upsampled.eval({encoded: code, positions: auxpos})
    

