def test_neibs_bad_shape_wrap_centered(self):
        shape = (2, 3, 10, 10)

        for dtype in self.dtypes:
            images = shared(numpy.arange(
                numpy.prod(shape), dtype=dtype
                ).reshape(shape))

            for neib_shape in [(3, 2), (2, 3)]:
                neib_shape = T.as_tensor_variable(neib_shape)

                f = function([], images2neibs(images, neib_shape,
                                              mode="wrap_centered"),
                             mode=self.mode)
                self.assertRaises(TypeError, f)

            for shape in [(2, 3, 2, 3), (2, 3, 3, 2)]:
                images = shared(numpy.arange(numpy.prod(shape)).reshape(shape))
                neib_shape = T.as_tensor_variable((3, 3))
                f = function([], images2neibs(images, neib_shape,
                                              mode="wrap_centered"),
                             mode=self.mode)
                self.assertRaises(TypeError, f)

            # Test a valid shapes
            shape = (2, 3, 3, 3)
            images = shared(numpy.arange(numpy.prod(shape)).reshape(shape))
            neib_shape = T.as_tensor_variable((3, 3))

            f = function([],
                         images2neibs(images, neib_shape, mode="wrap_centered"),
                         mode=self.mode)
            f()

 def test_neibs_half_step_by_valid(self):
     neib_shapes = ((3, 3), (3, 5), (5, 3))
     for shp_idx, (shape, neib_step) in enumerate([
         [(7, 8, 5, 5), (1, 1)],
         [(7, 8, 5, 5), (2, 2)],
         [(7, 8, 5, 5), (4, 4)],
         [(7, 8, 5, 5), (1, 4)],
         [(7, 8, 5, 5), (4, 1)],
         [(80, 90, 5, 5), (1, 2)],
         [(1025, 9, 5, 5), (2, 1)],
         [(1, 1, 5, 1037), (2, 4)],
         [(1, 1, 1045, 5), (4, 2)]]
     ):
         for neib_shape in neib_shapes:
             for dtype in self.dtypes:
                 x = theano.shared(np.random.randn(*shape).astype(dtype))
                 extra = (neib_shape[0] // 2, neib_shape[1] // 2)
                 padded_shape = (x.shape[0], x.shape[1], x.shape[2] + 2 * extra[0], x.shape[3] + 2 * extra[1])
                 padded_x = T.zeros(padded_shape)
                 padded_x = T.set_subtensor(padded_x[:, :, extra[0]:-extra[0], extra[1]:-extra[1]], x)
                 x_using_valid = images2neibs(padded_x, neib_shape, neib_step, mode="valid")
                 x_using_half = images2neibs(x, neib_shape, neib_step, mode="half")
                 f_valid = theano.function([], x_using_valid, mode='FAST_RUN')
                 f_half = theano.function([], x_using_half, mode=self.mode)
                 unittest_tools.assert_allclose(f_valid(), f_half())

    def test_infer_shape(self):
        shape = (100, 40, 6, 3)
        images = np.ones(shape).astype('float32')
        x = T.ftensor4()
        self._compile_and_check(
            [x], [images2neibs(x, neib_shape=(2, 1), mode='valid')],
            [images], Images2Neibs)
        self._compile_and_check(
            [x], [images2neibs(x, neib_shape=(2, 3), mode='valid')],
            [images], Images2Neibs)
        shape = (100, 40, 5, 4)
        images = np.ones(shape).astype('float32')
        x = T.ftensor4()
        self._compile_and_check(
            [x], [images2neibs(
                x, neib_shape=(2, 1), mode='ignore_borders')],
            [images], Images2Neibs)
        shape = (100, 40, 5, 3)
        images = np.ones(shape).astype('float32')
        x = T.ftensor4()
        self._compile_and_check(
            [x], [images2neibs(
                x, neib_shape=(2, 3), mode='ignore_borders')],
            [images], Images2Neibs)

        shape = (100, 40, 6, 7)
        images = np.ones(shape).astype('float32')
        x = T.ftensor4()
        self._compile_and_check(
            [x], [images2neibs(
                x, neib_shape=(2, 2), mode='ignore_borders')],
            [images], Images2Neibs)
        shape = (100, 40, 5, 10)
        images = np.ones(shape).astype('float32')
        x = T.ftensor4()
        self._compile_and_check(
            [x], [images2neibs(
                x, neib_shape=(3, 3), mode='wrap_centered')],
            [images], Images2Neibs)
        shape = (100, 40, 6, 4)
        images = np.ones(shape).astype('float32')
        x = T.ftensor4()
        self._compile_and_check(
            [x], [images2neibs(x, neib_shape=(2, 1), mode='half')],
            [images], Images2Neibs)
        self._compile_and_check(
            [x], [images2neibs(x, neib_shape=(2, 3), mode='half')],
            [images], Images2Neibs)
        shape = (100, 40, 6, 5)
        images = np.ones(shape).astype('float32')
        x = T.ftensor4()
        self._compile_and_check(
            [x], [images2neibs(x, neib_shape=(2, 1), mode='full')],
            [images], Images2Neibs)
        self._compile_and_check(
            [x], [images2neibs(x, neib_shape=(2, 3), mode='full')],
            [images], Images2Neibs)

def t_mk_pool_ready(t_pool_input, t_pool_shape):
    """
    Prepare pooling input
    :param t_pool_input: 4D theano tensor batch_sz x channels x height x width
    :param t_pool_shape: theano lvector pool_ch x pool_h x pool_w
    :return: aux. sizes and input reshaped for pooling
    """
    # sizes
    # input
    t_batch_sz = t_pool_input.shape[0]
    t_in_ch = t_pool_input.shape[1]
    t_in_h = t_pool_input.shape[2]
    t_in_w = t_pool_input.shape[3]
    # pooling
    t_pool_ch = t_pool_shape[0]
    t_pool_h = t_pool_shape[1]
    t_pool_w = t_pool_shape[2]
    # output
    t_out_ch = (t_in_ch + t_pool_ch - 1) // t_pool_ch
    t_out_h = (t_in_h + t_pool_h - 1) // t_pool_h
    t_out_w = (t_in_w + t_pool_w - 1) // t_pool_w

    # we will need to pad input (probably), so here's the padded shape:
    t_padded_ch = t_out_ch * t_pool_ch
    t_padded_h = t_out_h * t_pool_h
    t_padded_w = t_out_w * t_pool_w
    t_padded_pool_in_z = T.zeros(T.stack([t_batch_sz, t_padded_ch, t_padded_h, t_padded_w]))
    t_padded_pool_in = T.inc_subtensor(t_padded_pool_in_z[:t_batch_sz, :t_in_ch, :t_in_h, :t_in_w], t_pool_input)

    # below is all computed
    # spatial pooling
    t_sp_pooled = images2neibs(t_padded_pool_in, T.stack([t_pool_h, t_pool_w]))
    # spatial pooling output shape
    # has size (B * C * H/h * W/w) x (h*w)
    t_sp_pooled_dims = t_sp_pooled.shape
    # lines per channel
    # H*W / (h*w)
    t_lpc = (t_padded_h * t_padded_w) // (t_pool_h * t_pool_w)
    # shape to collect channels
    t_ch_pool_prep_dims_1 = T.stack([t_sp_pooled_dims[0] // t_lpc, t_lpc, t_sp_pooled_dims[1]])
    # preparing pooling by channels
    # reshape to collect channels in a separate dimension
    t_ch_pool_prep_1 = T.reshape(t_sp_pooled, t_ch_pool_prep_dims_1)
    t_ch_pool_prep_2 = T.shape_padleft(T.transpose(t_ch_pool_prep_1, [1, 0, 2]))
    # prepare for channel pooling
    t_ch_pool_dims = T.stack([t_pool_ch, t_ch_pool_prep_dims_1[-1]])
    t_pool_ready = images2neibs(t_ch_pool_prep_2, t_ch_pool_dims)
    return t_batch_sz, t_in_ch, t_in_h, t_in_w, t_out_ch, t_out_h, t_out_w, t_pool_ready

def _meanpool ( input, ds, ignore_border = False ):
    """ provide mean pooling """
    out_shp = (input.shape[0], input.shape[1], input.shape[2]/ds[0], input.shape[3]/ds[1])    
    neib = images2neibs(input, neib_shape = ds , 
                                mode = 'valid' if ignore_border is False else 'ignore_borders')
    pooled_vectors = neib.mean( axis = - 1 )
    return T.reshape(pooled_vectors, out_shp, ndim = 4 )    

    def get_output(self, train):
        X = self.get_input(train)
        # check if poolsize is symmetric. If not, step in neibs has to be set.

        if self.stride is not None:
            # rows_symmetrical = (X.shape[2] + 1)//2
            # step_val = (X.shape[2] - 1)//2
            sums = images2neibs(X, neib_shape=self.poolsize, neib_step=self.stride).sum(axis=-1)
            counts = T.neq(images2neibs(X, neib_shape=self.poolsize, neib_step=self.stride), 0).sum(axis=-1)
            average = (sums/counts).reshape((X.shape[0], X.shape[1], 2, 1))
        else:
            # rows_symmetrical = (X.shape[2])//2
            sums = images2neibs(X, neib_shape=self.poolsize).sum(axis=-1)
            counts = T.neq(images2neibs(X, neib_shape=self.poolsize), 0).sum(axis=-1)
            average = (sums/counts).reshape((X.shape[0], X.shape[1], 2, 1))
        return average

def make_patches_grid(x, patch_size, patch_stride):
    '''Break image `x` up into a grid of patches.

    input shape: (channels, rows, cols)
    output shape: (rows, cols, channels, patch_rows, patch_cols)
    '''
    from theano.tensor.nnet.neighbours import images2neibs  # TODO: all K, no T
    x = K.expand_dims(x, 0)
    xs = K.shape(x)
    num_rows = 1 + (xs[-2] - patch_size) // patch_stride
    num_cols = 1 + (xs[-1] - patch_size) // patch_stride
    num_channels = xs[-3]
    patches = images2neibs(
        x, (patch_size, patch_size), (patch_stride, patch_stride),
        mode='valid')
    # neibs are sorted per-channel
    patches = K.reshape(patches,
                        (num_channels, K.shape(patches)[0] // num_channels,
                         patch_size, patch_size))
    patches = K.permute_dimensions(patches, (1, 0, 2, 3))
    # arrange in a 2d-grid (rows, cols, channels, px, py)
    patches = K.reshape(
        patches, (num_rows, num_cols, num_channels, patch_size, patch_size))
    patches_norm = K.sqrt(
        K.sum(K.square(patches), axis=(2, 3, 4), keepdims=True))
    return patches, patches_norm

    def test_neibs(self):
        for shape, pshape in [((10, 7, 18, 18), (2, 2)),
                              ((10, 7, 6, 18), (3, 2)),
                              ((5, 7, 66, 66), (33, 33)),
                              ((5, 7, 68, 66), (34, 33))]:
            for border in ['valid', 'ignore_borders']:
                for dtype in self.dtypes:
                    images = shared(
                        numpy.arange(numpy.prod(shape), dtype=dtype).reshape(shape))
                    neib_shape = T.as_tensor_variable(pshape)

                    f = function([],
                                 images2neibs(images, neib_shape, mode=border),
                                 mode=self.mode)

                    # print images.get_value(borrow=True)
                    neibs = f()
                    # print neibs
                    g = function([],
                                 neibs2images(neibs, neib_shape, images.shape),
                                 mode=self.mode)
                    assert any([isinstance(node.op, self.op)
                                for node in f.maker.fgraph.toposort()])

                    # print g()
                    assert numpy.allclose(images.get_value(borrow=True), g())

def main():
    '''
    Read here all images(grayscale) from jaffe folder
    into an numpy array Ims with size (no_images, height, width).
    Make sure the images are read after sorting the filenames
    '''
    
    # didn't use walk
    from glob import glob
    import os
    
    files_list = glob(os.path.join('jaffe/', '*.tiff'))
    im_ind=0
    im = Image.open(files_list[0])
    Ims = np.float32(np.zeros([len(files_list),im.height,im.width]))
    for a_file in sorted(files_list):
        im = Image.open(a_file).convert("L")
        Ims[im_ind,:,:]=np.array(im)        
        im_ind = im_ind+1

    szs = [16, 32, 64]
    num_coeffs = [range(1, 10, 1), range(3, 30, 3), range(5, 50, 5)]

    for sz, nc in zip(szs, num_coeffs):
        '''
        Divide here each image into non-overlapping blocks of shape (sz, sz).
        Flatten each block and arrange all the blocks in a
        (no_images*n_blocks_in_image) x (sz*sz) matrix called X
        '''
        
        # Defining variables
        images = T.tensor4('images')
        neibs = images2neibs(images, neib_shape=(sz,sz))
        
        # Constructing theano function
        window_function = theano.function([images], neibs)
        
        X=window_function(Ims.reshape((1,len(files_list),im.height,im.width)))        
        
        X_mn = np.mean(X, 0)
        X = X - np.repeat(X_mn.reshape(1, -1), X.shape[0], 0)

        '''
        Perform eigendecomposition on X^T X and arrange the eigenvectors
        in decreasing order of eigenvalues into a matrix D
        '''
        V,D = np.linalg.eigh(np.dot(X.T,X))
        D = np.fliplr(D)        
        
        c = np.dot(D.T, X.T)

        for i in range(0, 200, 10):
            plot_mul(c, D, i, X_mn.reshape((sz, sz)),
                     num_coeffs=nc, n_blocks=int(256/sz))

        plot_top_16(D, sz, imname='output/hw1a_top16_{0}.png'.format(sz))

 def __init__(self, input1, input2):
     x1_sub = input1[:, :, 2:-2, 2:-2]
     x1_flatten = T.flatten(x1_sub)
     x1 = T.extra_ops.repeat(x1_flatten, 25)
     x1 = T.reshape(x1, [T.shape(x1_flatten)[0], 25])
     x2 = neighbours.images2neibs(input2, neib_shape=(5, 5), neib_step=(1, 1))
     diff = x1 - x2
     new_shape = T.shape(x1_sub)*[1, 1, 5, 5]
     diff_img = neighbours.neibs2images(diff, neib_shape=(5, 5), original_shape=[1, 25, 25*5, 5*5])
     self.output = T.nnet.relu(diff_img)

    def test_neibs_full_with_inconsistent_borders(self):
        shape = (2, 3, 5, 5)
        images = T.dtensor4()
        images_val = np.arange(np.prod(shape),
                               dtype='float32').reshape(shape)

        f = theano.function([images],
                            T.sqr(images2neibs(images, (2, 2), mode='full')),
                            mode=self.mode)
        self.assertRaises(TypeError, f, images_val)

def _maxrandpool ( input, ds, p, ignore_border = False ):
    """ provide random pooling among the top 'p' sorted outputs p = 0 is maxpool """
    rng = numpy.random.RandomState(24546)
    out_shp = (input.shape[0], input.shape[1], input.shape[2]/ds[0], input.shape[3]/ds[1])        
    srng = RandomStreams(rng.randint(2147462579))
    pos = srng.random_integers(size=(1,1), low = ds[0]*ds[1]-1-p, high = ds[0]*ds[1]-1)
    neib = images2neibs(input, neib_shape = ds ,
                                mode = 'valid' if ignore_border is False else 'ignore_borders') 
    neib = neib.sort(axis = -1) 
    pooled_vectors = neib[:,pos]   
    return T.reshape(pooled_vectors, out_shp, ndim = 4 )   

def make_patches(x, patch_size, patch_stride):
    from theano.tensor.nnet.neighbours import images2neibs
    x = K.expand_dims(x, 0)
    patches = images2neibs(x,
        (patch_size, patch_size), (patch_stride, patch_stride),
        mode='valid')
    # neibs are sorted per-channel
    patches = K.reshape(patches, (K.shape(x)[1], K.shape(patches)[0] // K.shape(x)[1], patch_size, patch_size))
    patches = K.permute_dimensions(patches, (1, 0, 2, 3))
    patches_norm = K.l2_normalize(patches, 1)
    return patches, patches_norm

    def speed_neibs(self):
        shape = (100, 40, 18, 18)
        images = shared(numpy.arange(numpy.prod(shape),
                                     dtype='float32').reshape(shape))
        neib_shape = T.as_tensor_variable((3, 3))

        f = function([], images2neibs(images, neib_shape),
                     mode=self.mode)

        for i in range(1000):
            f()

    def test_neibs_bad_shape(self):
        shape = (2, 3, 10, 10)
        for dtype in self.dtypes:
            images = shared(numpy.arange(
                numpy.prod(shape), dtype=dtype).reshape(shape))

            for neib_shape in [(3, 2), (2, 3)]:
                neib_shape = T.as_tensor_variable(neib_shape)
                f = function([], images2neibs(images, neib_shape),
                             mode=self.mode)
                self.assertRaises(TypeError, f)

                # Test that ignore border work in that case.
                f = function([],
                             images2neibs(images, neib_shape,
                                          mode='ignore_borders'),
                             mode=self.mode)
                assert self.op in [type(node.op)
                                   for node in f.maker.fgraph.toposort()]
                f()

def make_patches(x, shape):
    x = K.expand_dims(x, 0)

    patches = images2neibs(x, (shape, shape))
    patches = K.reshape(patches, (K.shape(x)[1],
                                  K.shape(patches)[0] / K.shape(x)[1],
                                  shape, shape))
    
    patches_norm = K.sqrt(K.sum(K.square(patches), axis=(1,2,3),
                                keepdims=True))

    return patches, patches_norm

def make_patches(x, patch_size, patch_stride):
    '''Break image `x` up into a bunch of patches.'''
    from theano.tensor.nnet.neighbours import images2neibs
    x = K.expand_dims(x, 0)
    patches = images2neibs(x,
        (patch_size, patch_size), (patch_stride, patch_stride),
        mode='valid')
    # neibs are sorted per-channel
    patches = K.reshape(patches, (K.shape(x)[1], K.shape(patches)[0] // K.shape(x)[1], patch_size, patch_size))
    patches = K.permute_dimensions(patches, (1, 0, 2, 3))
    patches_norm = K.sqrt(K.sum(K.square(patches), axis=(1,2,3), keepdims=True))
    return patches, patches_norm

    def test_neibs_valid_with_inconsistent_borders(self):
        shape = (2, 3, 5, 5)
        images = T.dtensor4()
        images_val = numpy.arange(numpy.prod(shape),
                                  dtype='float32').reshape(shape)

        def fn(images):
            return T.sum(T.sqr(images2neibs(images, (2, 2), mode='valid')),
                         axis=[0, 1])

        f = theano.function([images],
                            T.sqr(images2neibs(images, (2, 2), mode='valid')),
                            mode=self.mode)
        self.assertRaises(TypeError, f, images_val)

def main():
    '''
    Read here all images(grayscale) from jaffe folder
    into an numpy array Ims with size (no_images, height, width).
    Make sure the images are read after sorting the filenames
    '''
    imagename_list = os.listdir('jaffe/')
    imagename_list.sort()
    Ims = []
    for img_name in imagename_list:
        Ims.append(np.array(Image.open('jaffe/'+img_name)))
    Ims = np.asarray(Ims)
    print "Data loaded"
    nimgs, height, width = Ims.shape

    szs = [16, 32, 64]
    num_coeffs = [range(1, 10, 1), range(3, 30, 3), range(5, 50, 5)]

    for sz, nc in zip(szs, num_coeffs):
        '''
        Divide here each image into non-overlapping blocks of shape (sz, sz).
        Flatten each block and arrange all the blocks in a
        (no_images * n_blocks_in_image) x (sz*sz) matrix called X
        '''
        images = theano.tensor.tensor4('images')
        neibs = images2neibs(images, neib_shape=(sz, sz))
        window_function = theano.function([images], neibs, allow_input_downcast=True)
        X = window_function(Ims.reshape((nimgs, 1, height, width)))
        '''
        Perform eigendecomposition on X^T X and arrange the eigenvectors
        in decreasing order of eigenvalues into a matrix D
        '''
        print 'Computing eigendecomposition for size %s with nc = %s...' % (sz, nc)
        X_mn = np.mean(X, 0)
        X = X - np.repeat(X_mn.reshape(1, -1), X.shape[0], 0)
        eigen_out = np.linalg.eig(X.T.dot(X))
        eigenvalues = eigen_out[0]
        eigenvectors = eigen_out[1]
        # Sorted Eigenvalues
        col_indx = np.argsort(-eigenvalues)
        # Sorted Eigenvectors
        eigenvectors = eigenvectors[:,col_indx]
        D = eigenvectors[:, col_indx]
        c = np.dot(D.T, X.T)
        print 'Plotting the reconstruction...'
        for i in range(0, 200, 10): # Incrementing by 10
            plot_mul(c=c, D=D, im_num=i, X_mn=X_mn.reshape((sz, sz)),
                     num_coeffs=nc, n_blocks=int(256/sz))

        plot_top_16(D, sz, imname='output/hw1a_top16_{0}.png'.format(sz))

def img_2_neibs_with_chans(inputs_sym, patch_size):
    flat_patches = neighbours.images2neibs(inputs_sym, patch_size, (1,1))
    topo_flat_patches = T.reshape(flat_patches,(inputs_sym.shape[0],
                                            inputs_sym.shape[1],
                                            inputs_sym.shape[2]-patch_size[0]+1,
                                            inputs_sym.shape[3]-patch_size[1]+1,
                                            patch_size[0],
                                            patch_size[1]))


    flat_patches = topo_flat_patches.dimshuffle(0,2,3,1,4,5)
    flat_patches = T.reshape(flat_patches, (T.prod(flat_patches.shape[:3]),
                                                 T.prod(flat_patches.shape[3:])))
    return flat_patches