def test_view_as_windows_With_skip():
    A = np.arange(20).reshape((5, 4))
    B = view_as_windows(A, (2, 2), step=2)
    assert_equal(B, [[[[0, 1], [4, 5]], [[2, 3], [6, 7]]], [[[8, 9], [12, 13]], [[10, 11], [14, 15]]]])

    C = view_as_windows(A, (2, 2), step=4)
    assert_equal(C.shape, (1, 1, 2, 2))

def test_view_as_windows_optimal_step():
    A = np.arange(24).reshape((6, 4))
    B = view_as_windows(A, (3, 2), optimal_step=True)
    assert B.shape == (2, 2, 3, 2)
    assert B.size == A.size

    A = np.arange(512 * 512).reshape((512, 512))
    B = view_as_windows(A, (10, 10), optimal_step=True)
    assert B.shape == (56, 56, 10, 10)
    assert B.size >= A.size

    C = view_as_windows(A, (11, 9), optimal_step=True)
    assert C.shape == (51, 63, 11, 9)
    assert C.size >= A.size

    def perform(self, node, inputs, output_storage):
        image_error, features = inputs
        kshp = self.kshp
        imshp = image_error.shape
        strides = self.strides

#        scipy_error_rot = np.transpose(image_error, [1, 0, 2, 3])
#        features_rot = np.transpose(features, [1, 0, 2, 3])
#
#        feat_expand = np.zeros((features_rot.shape[0],
#                                features_rot.shape[1],
#                                image_error.shape[2] - self.kshp[2] + 1,
#                                image_error.shape[3] - self.kshp[3] + 1), dtype=features.dtype)
#        feat_expand[:, :, ::strides[0], ::strides[1]] = features_rot
#
#        #scipy_derivative_rot = -scipy_convolve4d(scipy_error_rot[:, :, ::-1, ::-1], feat_expand)
#
#        feat_flipped = features_rot[:, :, ::-1, ::-1]

        from skimage.util.shape import view_as_windows

        image_error_view = view_as_windows(image_error, (imshp[0], kshp[1], kshp[2], kshp[3]))[0,0,::strides[0],::strides[1],...]
        # image_error_view.shape = (featszr, featszc, num_im, channels, ksz[2], ksz[3])
        # features.shape = (num_im, num_filters, featszr, featszc)
        kernel_derivative = - np.tensordot(features,image_error_view, axes=((0, 2, 3), (2, 0, 1)))
        # kernel_derivative_temp.shape = (num_filters, channels, ksz[2], ksz[3])

        output_storage[0][0] = kernel_derivative[:,:,::-1,::-1]

    def stochastic_pooling(self, X, n_layer):

        inh = X.shape[0] - self.shape_pool[0] + 1
        inw = X.shape[1] - self.shape_pool[1] + 1
        n_filter = self.n_filters[n_layer]
        filtersize = self.shape_pool[0] * self.shape_pool[1]

        randomsamples = random_sample((inh) * (inw) * n_filter).reshape(
            (inh), (inw), n_filter
        )  # generate random values
        randomsamples = np.repeat(randomsamples, repeats=filtersize, axis=2).reshape((inh), (inw), n_filter, filtersize)

        X_rfi = view_as_windows(X, self.shape_pool + (1,))
        sumpool = np.repeat(np.sum(X_rfi, axis=(3, 4, 5)), repeats=filtersize).reshape(
            (inh, inw, n_filter, self.shape_pool[0], self.shape_pool[1], 1)
        )
        probabilities = X_rfi / sumpool
        probabilities[np.isnan(probabilities)] = 1 / float(
            filtersize
        )  # get where the sum is zero and replace by one, so the division by zero error do not occur
        probabilities = probabilities.reshape((inh, inw, n_filter, filtersize))

        if self.training:

            bins = np.add.accumulate(probabilities, axis=3)
            binsbefore = np.concatenate((np.zeros((inh, inw, n_filter, 1)), bins[:, :, :, :-1]), axis=3)
            ret = X_rfi[np.where((((binsbefore <= randomsamples) * (bins > randomsamples))) == True)]
            ret = ret.reshape(((inh), (inw), n_filter))[:: self.stride_pool, :: self.stride_pool]

        else:  # for testing
            ret = probabilities * X_rfi
            sumpool[sumpool == 0.0] = 1.0
            ret = np.sum(ret, axis=(3)) / sumpool[:, :, :, 0, 0, 0]
            ret = ret[:: self.stride_pool, :: self.stride_pool]

    def get_pooled_features(self, input_feature_map, filter_size=(19,19)):
        # assuming square filters and images
        filter_side = filter_size[0]

        # reshaping incoming features from 2d to 3d i.e. (3249,20) to (57,57,20)
        input_feature_map_shape = input_feature_map.shape
        if input_feature_map.ndim == 2:
            input_feature_map_side = int(np.sqrt(input_feature_map.shape[0]))
            input_feature_map = input_feature_map.reshape((input_feature_map_side, input_feature_map_side, input_feature_map_shape[-1]))
        assert input_feature_map.ndim == 3, "Input features dimension is %d instead of 3" %input_feature_map.ndim

        # get windows (57,57,20) to (3,3,1,19,19,20)
        input_feature_map_windows = view_as_windows(input_feature_map,
                                                    window_shape=(filter_size[0], filter_size[1], input_feature_map.shape[-1]),
                                                    step=filter_size[0])

        # reshape windows (3,3,1,19,19,20) to (3**2, 19**2, 20) == (9, 361, 20)
        input_feature_map_windows = input_feature_map_windows.reshape((input_feature_map_windows.shape[0]**2,
                                                                       filter_size[0]**2,
                                                                       input_feature_map.shape[-1]))

        # calculate norms (9, 361, 20) to (9,361)
        input_feature_map_window_norms = np.linalg.norm(input_feature_map_windows, ord=2, axis=-1)

        # calculate indexes of max norms per window (9,361) to (9,1). One max index per window.
        max_norm_indexes = np.argmax(input_feature_map_window_norms, axis=-1)

        # max pooled features are the features that have max norm indexes (9, 361, 20) to (9,20). One max index per window.
        pooled_features = input_feature_map_windows[np.arange(input_feature_map_windows.shape[0]), max_norm_indexes]

        # return pooled feature map
        return pooled_features

def lpool4(arr_in, neighborhood,
           order=DEFAULT_ORDER,
           stride=DEFAULT_STRIDE, arr_out=None):
    """4D Local Pooling Operation

    XXX: docstring
    """
    assert arr_in.ndim == 4
    assert len(neighborhood) == 2

    order = np.array([order], dtype=arr_in.dtype)
    stride = np.int(stride)

    in_imgs, inh, inw, ind = arr_in.shape
    nbh, nbw = neighborhood
    assert nbh <= inh
    assert nbw <= inw

    if arr_out is not None:
        assert arr_out.dtype == arr_in.dtype
        assert arr_out.shape == (in_imgs,
                                 1 + (inh - nbh) / stride,
                                 1 + (inw - nbw) / stride,
                                 ind)

    _arr_out = ne.evaluate('arr_in ** order')
    _arr_out = view_as_windows(_arr_out, (1, 1, nbw, 1))
    _arr_out = ne.evaluate('sum(_arr_out, 6)')[:, :, ::stride, :, 0, 0, 0]
    # np.ascontiguousarray(_arr_out) necessary to avoid
    # skimage/util/shape.py:237: RuntimeWarning: Cannot provide views on a non-contiguous input array without copying.
    # warn(RuntimeWarning("Cannot provide views on a non-contiguous input 
    _arr_out = np.ascontiguousarray(_arr_out)
    _arr_out = view_as_windows(_arr_out, (1, nbh, 1, 1))
    _arr_out = ne.evaluate('sum(_arr_out, 5)')[:, ::stride, :, :, 0, 0, 0]

    _arr_out = ne.evaluate('_arr_out ** (1 / order)')

    if arr_out is not None:
        arr_out[:] = _arr_out
    else:
        arr_out = _arr_out

    assert arr_out.shape[0] == in_imgs

    assert arr_out.dtype == arr_in.dtype

    return arr_out

    def similarity3D(self, X, fb):
        assert X.ndim == 3
        assert fb.ndim == 4
        assert X.shape[-1] == fb.shape[2]

        Xw = view_as_windows(X, fb.shape[:3])
        Xwa = np.abs(Xw - fb)
        return Xwa.sum(axis=(3, 4, 5))

    def _process_one_op(self, arr, X, Y,
                        layer_idx, op_idx, kwargs,
                        op_params,
                        op_name, nbh, nbw, stride,
                        pad_apron=False,
                        interleave_stride=False):

        out_l = []

        # -- here we compute the pixel coordinates of
        #    the central pixel in a patch
        hc, wc = nbh / 2, nbw / 2

        if pad_apron:

            arr = filter_pad2d(arr, (nbh, nbw))
            X = np.squeeze(filter_pad2d(X[..., np.newaxis], (nbh, nbw),
                                        constant=-1))
            Y = np.squeeze(filter_pad2d(Y[..., np.newaxis], (nbh, nbw),
                                        constant=-1))

        if interleave_stride:

            for i in xrange(stride):
                for j in xrange(stride):

                    arr_out_ij = self._get_feature_map(arr[i::, j::, ...],
                                                  layer_idx, op_idx, kwargs,
                                                  op_params, op_name)
                    X_out_ij = view_as_windows(X[i::, j::],
                                               (nbh, nbw))[::stride, ::stride,
                                                           hc, wc]
                    Y_out_ij = view_as_windows(Y[i::, j::],
                                               (nbh, nbw))[::stride, ::stride,
                                                           hc, wc]
                    out_l += [(arr_out_ij, X_out_ij, Y_out_ij)]
        else:

            arr_out = self._get_feature_map(arr, layer_idx, op_idx, kwargs,
                                            op_params, op_name)
            X_out = view_as_windows(X, (nbh, nbw))[::stride, ::stride, hc, wc]
            Y_out = view_as_windows(Y, (nbh, nbw))[::stride, ::stride, hc, wc]
            out_l += [(arr_out, X_out, Y_out)]

        return out_l

def mcconv3(X, W):
    X_VAW = view_as_windows(X, W.shape[0:-1])
    Y_FPS = X_VAW.shape[0:2]
    X_VAW = X_VAW.reshape(Y_FPS[0] * Y_FPS[1], -1)
    W = W.reshape(-1, W.shape[-1])
    Y = np.dot(X_VAW, W)
    Y = Y.reshape(Y_FPS[0], Y_FPS[1], -1)

    return Y

def bxfilt2(X, F_SIZ, F_STRD):
    for i in reversed(xrange(2)):
        W_SIZ = np.ones(np.ndim(X))
        S_SIZ = np.ones(2)
        W_SIZ[i], S_SIZ[i] = F_SIZ, F_STRD
        X = np.squeeze(view_as_windows(X, tuple(W_SIZ)))[::S_SIZ[0], ::S_SIZ[1]]  # subsampling before summation
        X = np.sum(X, -1)

    return X

def lpool4(arr_in, neighborhood,
           order=DEFAULT_ORDER,
           stride=DEFAULT_STRIDE, arr_out=None):
    """4D Local Pooling Operation

    XXX: docstring
    """
    assert arr_in.ndim == 4
    assert len(neighborhood) == 2

    order = np.array([order], dtype=arr_in.dtype)
    stride = np.int(stride)

    in_imgs, inh, inw, ind = arr_in.shape
    nbh, nbw = neighborhood
    assert nbh <= inh
    assert nbw <= inw

    if arr_out is not None:
        assert arr_out.dtype == arr_in.dtype
        assert arr_out.shape == (in_imgs,
                                 1 + (inh - nbh) / stride,
                                 1 + (inw - nbw) / stride,
                                 ind)

    _arr_out = ne.evaluate('arr_in ** order')
    _arr_out = view_as_windows(_arr_out, (1, 1, nbw, 1))
    _arr_out = ne.evaluate('sum(_arr_out, 6)')[:, :, ::stride, :, 0, 0, 0]
    _arr_out = view_as_windows(_arr_out, (1, nbh, 1, 1))
    _arr_out = ne.evaluate('sum(_arr_out, 5)')[:, ::stride, :, :, 0, 0, 0]

    _arr_out = ne.evaluate('_arr_out ** (1 / order)')

    if arr_out is not None:
        arr_out[:] = _arr_out
    else:
        arr_out = _arr_out

    assert arr_out.shape[0] == in_imgs

    assert arr_out.dtype == arr_in.dtype

    return arr_out

def extractPatches(name, patch_shape, dataDir, vector, numberOfPatchesPerImage):
    imageName = '{0:s}{1:s}'.format(dataDir, name)
    npImage = cv2.imread(imageName)
    npImage = np.divide(npImage, vector)
    #npImage = cv2.resize(npImage, (65, 65))
    patchesSampled = np.empty(shape=(numberOfPatchesPerImage, patch_shape[0] * patch_shape[1], 3))
    for derp in range(npImage.shape[2]):
        patches = view_as_windows(npImage[:,:,derp], patch_shape)
        patches = patches.reshape(-1, patch_shape[0] * patch_shape[1])[::8]
        patchesSampled[:,:,derp] = patches[np.random.random_integers(0, patches.shape[0], numberOfPatchesPerImage), :]
    return(patchesSampled)

def lpool3(arr_in, neighborhood,
           order=DEFAULT_ORDER,
           stride=DEFAULT_STRIDE, arr_out=None):
    """3D Local Pooling Operation

    XXX: docstring
    """
    assert arr_in.ndim == 3
    assert len(neighborhood) == 2

    order = np.array([order], dtype=arr_in.dtype)
    stride = np.int(stride)

    inh, inw, ind = arr_in.shape
    nbh, nbw = neighborhood
    assert nbh <= inh
    assert nbw <= inw

    if arr_out is not None:
        assert arr_out.dtype == arr_in.dtype
        assert arr_out.shape == (1 + (inh - nbh) / stride,
                                 1 + (inw - nbw) / stride,
                                 ind)

    _arr_out = ne.evaluate('arr_in ** order')
    _arr_out = view_as_windows(_arr_out, (1, nbw, 1))
    _arr_out = ne.evaluate('sum(_arr_out, 4)')[:, ::stride, :, 0, 0]
    _arr_out = view_as_windows(_arr_out, (nbh, 1, 1))
    _arr_out = ne.evaluate('sum(_arr_out, 3)')[::stride, :, :, 0, 0]
    # Note that you need to use '1' and not '1.0' so that the dtype of
    # the exponent does not change (i.e. get promoted)
    _arr_out = ne.evaluate('_arr_out ** (1 / order)')

    if arr_out is not None:
        arr_out[:] = _arr_out
    else:
        arr_out = _arr_out

    assert arr_out.dtype == arr_in.dtype

    return arr_out

def test_view_as_windows_1D():
    A = np.arange(10)
    window_shape = (3,)
    B = view_as_windows(A, window_shape)
    assert_equal(B, np.array([[0, 1, 2],
                              [1, 2, 3],
                              [2, 3, 4],
                              [3, 4, 5],
                              [4, 5, 6],
                              [5, 6, 7],
                              [6, 7, 8],
                              [7, 8, 9]]))

def test_view_as_windows_step_tuple():
    A = np.arange(24).reshape((6, 4))
    B = view_as_windows(A, (3, 2), step=3)
    assert B.shape == (2, 1, 3, 2)
    assert B.size != A.size

    C = view_as_windows(A, (3, 2), step=(3, 2))
    assert C.shape == (2, 2, 3, 2)
    assert C.size == A.size

    assert_equal(C, [[[[0,  1],
                       [4,  5],
                       [8,  9]],
                      [[2,  3],
                       [6,  7],
                       [10, 11]]],
                     [[[12, 13],
                         [16, 17],
                         [20, 21]],
                      [[14, 15],
                         [18, 19],
                         [22, 23]]]])

    def conv3D(self, X, fb):
        assert X.ndim == 3
        assert fb.ndim == 4
        assert X.shape[-1] == fb.shape[2]

        n_filters = fb.shape[-1]
        fb2 = fb.copy()
        fb2.shape = -1, n_filters

        X_rfi = view_as_windows(X, fb.shape[:3])
        outh, outw = X_rfi.shape[:2]
        X_rfim = X_rfi.reshape(outh * outw, -1)
        ret = np.dot(X_rfim, fb2)  # -- convolution as matrix multiplication
        return ret.reshape(outh, outw, n_filters)

    def gen_neg_examples(self):
        """
        Generates negative training patches
        from the dataset specified in the configuration
        """

        negfiles = os.listdir(self.negdatadir)
        negpatches = []

        for n in negfiles[:10]:
            f = self.negdatadir + "/" + n

            if "jpg" not in f:
                continue
            print "Processing image" + f
            im = cv2.imread(f)
            # gray = self.rgb2gray(im)
            gray = cv2.cvtColor(im, cv2.COLOR_BGR2GRAY)
            # im =cv2.cvtColor(im,cv2.COLOR_BGR2GRAY)

            # Generate the negative patch
            # patch = image.extract_patches_2d(im, (256, 256)) # This generalates nultiple patches

            print np.shape(im)
            patch = view_as_windows(gray, (64, 64), 65)
            print np.shape(patch)
            # patch = [cv2.resize(x,(64,64)) for x in patch]

            # randomly choose 20 patches from all patches for this image
            fdimlen = np.shape(patch)[0]
            sdimlen = np.shape(patch)[1]
            patches = []

            # for i in range(len(self.trainpatches)):
            for i in range(8):  # we are getting too many negative examples, so reducing to 5

                randomdim = random.choice(range(fdimlen))
                randomsdim = random.choice(range(sdimlen))
                # randompatch = patch[randomdim][randomsdim][0] # if color un comment this
                randompatch = patch[randomdim][randomsdim]
                # randompatch = cv2.resize(randompatch, (64,64))
                patches.append(randompatch)

            negpatches += patches
            self.negpatches = negpatches

def montage_wb_ratio (input_image, patch_shape, n_filters, ele_print = False):

    from skimage.util.shape import view_as_windows
    from skimage.util.montage import montage2d
    from scipy.cluster.vq import kmeans2

    patches = view_as_windows(gray_crab, patch_shape)

    patches = patches.reshape(-1, patch_shape[0] * patch_shape[1])[::8]

    fb, _ = kmeans2(patches, n_filters, minit='points')

    fb = fb.reshape((-1,) + patch_shape)

    fb_montage = montage2d(fb, fill=False, rescale_intensity=True)

    elements = np.split(np.hstack(np.split(fb_montage, 4)), 16, axis=1)

    del elements[n_filters:]

    wb_ratios = []

    bin_elements = []

    for element in elements:

        thresh = threshold_otsu(element)

        binary = element > thresh

        ratio = np.sum(binary) / binary.size

        wb_ratios.append(ratio)

        if ele_print:
            bin_elements.append(binary)

    wb_ratios = sorted(wb_ratios, reverse = True)

    if ele_print:

        show_images(elements)
        show_images(bin_elements)

    return(wb_ratios)

def test_view_as_windows_2D():

    A = np.arange(5 * 4).reshape(5, 4)
    window_shape = (4, 3)
    B = view_as_windows(A, window_shape)
    assert_equal(B.shape, (2, 2, 4, 3))
    assert_equal(
        B,
        np.array(
            [
                [[[0, 1, 2], [4, 5, 6], [8, 9, 10], [12, 13, 14]], [[1, 2, 3], [5, 6, 7], [9, 10, 11], [13, 14, 15]]],
                [
                    [[4, 5, 6], [8, 9, 10], [12, 13, 14], [16, 17, 18]],
                    [[5, 6, 7], [9, 10, 11], [13, 14, 15], [17, 18, 19]],
                ],
            ]
        ),
    )

def scipy_correlate4d_view(features,kernel,stride=1):

    #imshp = (3,2,20,20) # num images, channels, szy, szx
    #kshp = (10,2,10,10) # features, channels, szy, szx
    #featshp = (3,10,11,11) # num images, features, szy, szx

    from skimage.util.shape import view_as_windows

    featshp = features.shape
    kshp = kernel.shape

    output_sz = (featshp[0], kshp[1], kshp[2] + stride*featshp[2] - 1, kshp[3] + stride*featshp[3] - 1)

    k_rot = kernel[:,:,::-1,::-1]

    scipy_output = np.zeros(output_sz)
    for im_i in range(featshp[0]):

        im_out = np.zeros(output_sz[1:])
        im_outr = view_as_windows(im_out,(kshp[1],kshp[2],kshp[3]))[0,::stride,::stride,...]

        im_hatr = np.tensordot(np.squeeze(features[im_i,...]),k_rot,axes=((0,),(0,)))
        #import IPython; ipshell = IPython.embed; ipshell(banner1='ipshell')

#        for a in range(im_hatr.shape[0]):
#            im_outr[a,:,...] += im_hatr[a,:,...]

        for a in range(im_hatr.shape[0]):
            for b in range(im_hatr.shape[1]):
                im_outr[a,b,...] += im_hatr[a,b,...]

#        for d in range(im_hatr.shape[3]):
#            for e in range(im_hatr.shape[4]):
#                im_outr[:,:,:,d,e] += im_hatr[:,:,:,d,e]

#                for c in range(im_hatr.shape[2]):
#                    for d in range(im_hatr.shape[3]):
#                        for e in range(im_hatr.shape[4]):
#                            im_outr[a,b,c,d,e] += im_hatr[a,b,c,d,e]
#        im_outr += im_hatr

        scipy_output[im_i,...] = im_out

    return scipy_output