def test_subclasses():
    # test that subclass is preserved only if subok=True
    a = VerySimpleSubClass([1, 2, 3, 4])
    assert_(type(a) is VerySimpleSubClass)
    a_view = as_strided(a, shape=(2,), strides=(2 * a.itemsize,))
    assert_(type(a_view) is np.ndarray)
    a_view = as_strided(a, shape=(2,), strides=(2 * a.itemsize,), subok=True)
    assert_(type(a_view) is VerySimpleSubClass)
    # test that if a subclass has __array_finalize__, it is used
    a = SimpleSubClass([1, 2, 3, 4])
    a_view = as_strided(a, shape=(2,), strides=(2 * a.itemsize,), subok=True)
    assert_(type(a_view) is SimpleSubClass)
    assert_(a_view.info == 'simple finalized')

    # similar tests for broadcast_arrays
    b = np.arange(len(a)).reshape(-1, 1)
    a_view, b_view = broadcast_arrays(a, b)
    assert_(type(a_view) is np.ndarray)
    assert_(type(b_view) is np.ndarray)
    assert_(a_view.shape == b_view.shape)
    a_view, b_view = broadcast_arrays(a, b, subok=True)
    assert_(type(a_view) is SimpleSubClass)
    assert_(a_view.info == 'simple finalized')
    assert_(type(b_view) is np.ndarray)
    assert_(a_view.shape == b_view.shape)

    # and for broadcast_to
    shape = (2, 4)
    a_view = broadcast_to(a, shape)
    assert_(type(a_view) is np.ndarray)
    assert_(a_view.shape == shape)
    a_view = broadcast_to(a, shape, subok=True)
    assert_(type(a_view) is SimpleSubClass)
    assert_(a_view.info == 'simple finalized')
    assert_(a_view.shape == shape)

def view(A, offset=0):
    """
    Get a view on the diagonal elements of a 2D array.

    This is actually a view (!) on the diagonal of the array, so you can
    in-place adjust the view.

    :param :class:`ndarray` A: 2 dimensional numpy array
    :param int offset: view offset to give back (negative entries allowed)
    :rtype: :class:`ndarray` view of diag(A)

    >>> import numpy as np
    >>> X = np.arange(9).reshape(3,3)
    >>> view(X)
    array([0, 4, 8])
    >>> d = view(X)
    >>> d += 2
    >>> view(X)
    array([ 2,  6, 10])
    >>> view(X, offset=-1)
    array([3, 7])
    >>> subtract(X, 3, offset=-1)
    array([[ 2,  1,  2],
           [ 0,  6,  5],
           [ 6,  4, 10]])
    """
    from numpy.lib.stride_tricks import as_strided
    assert A.ndim == 2, "only implemented for 2 dimensions"
    assert A.shape[0] == A.shape[1], "attempting to get the view of non-square matrix?!"
    if offset > 0:
        return as_strided(A[0, offset:], shape=(A.shape[0] - offset, ), strides=((A.shape[0]+1)*A.itemsize, ))
    elif offset < 0:
        return as_strided(A[-offset:, 0], shape=(A.shape[0] + offset, ), strides=((A.shape[0]+1)*A.itemsize, ))
    else:
        return as_strided(A, shape=(A.shape[0], ), strides=((A.shape[0]+1)*A.itemsize, ))

 def meshgrid_as_strided(x, y, mask=None):
     """numpy.meshgrid without copying data (using as_strided)"""        
     if mask is None:            
         return (as_strided(x, strides=(0, x.strides[0]), shape=(y.size, x.size)),
                 as_strided(y, strides=(y.strides[0],0), shape=(y.size, x.size)))    
     else:            
         return (np.ma.array(as_strided(x, strides=(0, x.strides[0]), shape=(y.size, x.size)), mask=mask),
                 np.ma.array(as_strided(y, strides=(y.strides[0],0), shape=(y.size, x.size)), mask=mask))

def broadcasted_shape(shp1, shp2):
    # determine shape of array of shp1 and shp2 broadcast against one another.
    x = np.array([1])
    # trick to define array with certain shape that doesn't allocate all the
    # memory.
    a = as_strided(x, shape=shp1, strides=[0] * len(shp1))
    b = as_strided(x, shape=shp2, strides=[0] * len(shp2))
    return np.broadcast(a, b).shape

def disparity_ssd(L, R, window_size = 21):
    """Compute disparity map D(y, x) such that: L(y, x) = R(y, x + D(y, x))
    Params:
    L: Grayscale left image, in range [0.0, 1.0]
    R: Grayscale right image, same size as L
    Returns: Disparity map, same size as L, R
    """

    D = np.zeros(L.shape, dtype=np.float)

    # subtract 1 due to the starting pixel
    offset = (window_size) / 2
    L = cv2.copyMakeBorder(L, offset, offset, offset, offset, cv2.BORDER_CONSTANT,value=0)
    R = cv2.copyMakeBorder(R, offset, offset, offset, offset, cv2.BORDER_CONSTANT,value=0)

    shape = L.shape
    height = shape[0]
    width = shape[1]

    r_shape = (R.shape[0]-(window_size-1), R.shape[1]-(window_size-1), window_size, window_size)
    r_strides = (R.shape[1] * R.itemsize, R.itemsize, R.itemsize * R.shape[1], R.itemsize)
    r_strips = as_strided(R, r_shape, r_strides)

    for y in range(offset, height - offset):
        r_strip = r_strips[y-offset]
        for x in range(offset, width-offset):
            l_patch = get_patch(L, offset, offset, offset, offset, y, x)

            copy_patch = np.copy(l_patch)
            l_strip = as_strided(copy_patch, r_strip.shape, (0, copy_patch.itemsize*window_size, copy_patch.itemsize))
            ssd = ((l_strip - r_strip)**2).sum((1, 2))

            x_prime = np.argmin(ssd)
            D[y-offset][x-offset] = x_prime - x

    #print D.max()
    return D

# def test_disparity_ssd2(l_image, r_image, problem, window_size = 21):
#     L = cv2.imread(os.path.join('input', l_image), 0) * (1 / 255.0)  # grayscale, scale to [0.0, 1.0]
#     R = cv2.imread(os.path.join('input', r_image), 0) * (1 / 255.0)
#
#     # Compute disparity (using method disparity_ssd defined in disparity_ssd.py)
#     start = time.time()
#     D = disparity_ssd(L, R, window_size)  # TODO# : implemenet disparity_ssd()
#     print "first: " + str(time.time() - start)
#     start = time.time()
#     D2 = disparity_ssd_2(L, R, window_size)
#     print "second: " + str(time.time() - start)


    #print D == D2

    cv2.imwrite(os.path.join("output", "ps3-" + problem + ".png"), np.clip(D2, 0, 255).astype(np.uint8))

 def demosaic(self):
     if self._demo is None:
         # XXX Again, should take into account camera's vflip and hflip here
         # Construct representation of the bayer pattern
         bayer = np.zeros(self.array.shape, dtype=np.uint8)
         bayer[1::2, 0::2, 0] = 1 # Red
         bayer[0::2, 0::2, 1] = 1 # Green
         bayer[1::2, 1::2, 1] = 1 # Green
         bayer[0::2, 1::2, 2] = 1 # Blue
         # Allocate output array with same shape as data and set up some
         # constants to represent the weighted average window
         window = (3, 3)
         borders = (window[0] - 1, window[1] - 1)
         border = (borders[0] // 2, borders[1] // 2)
         # Pad out the data and the bayer pattern (np.pad is faster but
         # unavailable on the version of numpy shipped with Raspbian at the
         # time of writing)
         rgb = np.zeros((
             self.array.shape[0] + borders[0],
             self.array.shape[1] + borders[1],
             self.array.shape[2]), dtype=self.array.dtype)
         rgb[
             border[0]:rgb.shape[0] - border[0],
             border[1]:rgb.shape[1] - border[1],
             :] = self.array
         bayer_pad = np.zeros((
             self.array.shape[0] + borders[0],
             self.array.shape[1] + borders[1],
             self.array.shape[2]), dtype=bayer.dtype)
         bayer_pad[
             border[0]:bayer_pad.shape[0] - border[0],
             border[1]:bayer_pad.shape[1] - border[1],
             :] = bayer
         bayer = bayer_pad
         # For each plane in the RGB data, construct a view over the plane
         # of 3x3 matrices. Then do the same for the bayer array and use
         # Einstein summation to get the weighted average
         self._demo = np.empty(self.array.shape, dtype=self.array.dtype)
         for plane in range(3):
             p = rgb[..., plane]
             b = bayer[..., plane]
             pview = as_strided(p, shape=(
                 p.shape[0] - borders[0],
                 p.shape[1] - borders[1]) + window, strides=p.strides * 2)
             bview = as_strided(b, shape=(
                 b.shape[0] - borders[0],
                 b.shape[1] - borders[1]) + window, strides=b.strides * 2)
             psum = np.einsum('ijkl->ij', pview)
             bsum = np.einsum('ijkl->ij', bview)
             self._demo[..., plane] = psum // bsum
     return self._demo

def test_internal_overlap_manual():
    # Stride tricks can construct arrays with internal overlap

    # We don't care about memory bounds, the array is not
    # read/write accessed
    x = np.arange(1).astype(np.int8)

    # Check low-dimensional special cases

    check_internal_overlap(x, False) # 1-dim
    check_internal_overlap(x.reshape([]), False) # 0-dim

    a = as_strided(x, strides=(3, 4), shape=(4, 4))
    check_internal_overlap(a, False)

    a = as_strided(x, strides=(3, 4), shape=(5, 4))
    check_internal_overlap(a, True)

    a = as_strided(x, strides=(0,), shape=(0,))
    check_internal_overlap(a, False)

    a = as_strided(x, strides=(0,), shape=(1,))
    check_internal_overlap(a, False)

    a = as_strided(x, strides=(0,), shape=(2,))
    check_internal_overlap(a, True)

    a = as_strided(x, strides=(0, -9993), shape=(87, 22))
    check_internal_overlap(a, True)

    a = as_strided(x, strides=(0, -9993), shape=(1, 22))
    check_internal_overlap(a, False)

    a = as_strided(x, strides=(0, -9993), shape=(0, 22))
    check_internal_overlap(a, False)

def test_as_strided():
    a = np.array([None])
    a_view = as_strided(a)
    expected = np.array([None])
    assert_array_equal(a_view, np.array([None]))

    a = np.array([1, 2, 3, 4])
    a_view = as_strided(a, shape=(2,), strides=(2 * a.itemsize,))
    expected = np.array([1, 3])
    assert_array_equal(a_view, expected)

    a = np.array([1, 2, 3, 4])
    a_view = as_strided(a, shape=(3, 4), strides=(0, 1 * a.itemsize))
    expected = np.array([[1, 2, 3, 4], [1, 2, 3, 4], [1, 2, 3, 4]])
    assert_array_equal(a_view, expected)

def moving_avg(a, halfwindow, mask=None):
    """
    Performs a fast n-point moving average of (the last
    dimension of) array *a*, by using stride tricks to roll
    a window on *a*.

    Note that *halfwindow* gives the nb of points on each side,
    so that n = 2*halfwindow + 1.

    If *mask* is provided, values of *a* where mask = False are
    skipped.

    Returns an array of same size as *a* (which means that near
    the edges, the averaging window is actually < *npt*).
    """
    # padding array with zeros on the left and on the right:
    # e.g., if halfwindow = 2:
    # a_padded    = [0 0 a0 a1 ... aN 0 0]
    # mask_padded = [F F ?  ?      ?  F F]

    if mask is None:
        mask = np.ones_like(a, dtype='bool')

    zeros = np.zeros(a.shape[:-1] + (halfwindow,))
    falses = zeros.astype('bool')

    a_padded = np.concatenate((zeros, np.where(mask, a, 0), zeros), axis=-1)
    mask_padded = np.concatenate((falses, mask, falses), axis=-1)

    # rolling window on padded array using stride trick
    #
    # E.g., if halfwindow=2:
    # rolling_a[:, 0] = [0   0 a0 a1 ...    aN]
    # rolling_a[:, 1] = [0  a0 a1 a2 ... aN 0 ]
    # ...
    # rolling_a[:, 4] = [a2 a3 ...    aN  0  0]

    npt = 2 * halfwindow + 1  # total size of the averaging window
    rolling_a = as_strided(a_padded,
                           shape=a.shape + (npt,),
                           strides=a_padded.strides + (a.strides[-1],))
    rolling_mask = as_strided(mask_padded,
                              shape=mask.shape + (npt,),
                              strides=mask_padded.strides + (mask.strides[-1],))

    # moving average
    n = rolling_mask.sum(axis=-1)
    return np.where(n > 0, rolling_a.sum(axis=-1).astype('float') / n, np.nan)

def kron_id_view(vec, id_length, axis=-1):
    shape = (vec.shape[:axis] +
             (vec.shape[axis] - id_length + 1, id_length) +
             vec.shape[axis % vec.ndim + 1:])
    strides = vec.strides[:axis] + (vec.strides[axis],) + vec.strides[axis:]

    return as_strided(vec, shape=shape, strides=strides)

    def __init__(self, which_set, context_len, data_mode, shuffle=True):
        self.__dict__.update(locals())
        del self.self

        # Load data into self._data (defined in PennTreebank)
        self._load_data(which_set, context_len, data_mode)

        self._data = as_strided(self._raw_data,
                                shape=(len(self._raw_data) - context_len,
                                       context_len + 1),
                                strides=(self._raw_data.itemsize,
                                         self._raw_data.itemsize))

        super(PennTreebankNGrams, self).__init__(
            X=self._data[:, :-1],
            y=self._data[:, -1:],
            X_labels=self._max_labels, y_labels=self._max_labels
        )

        if shuffle:
            warnings.warn("Note that the PennTreebank samples are only "
                          "shuffled when the iterator method is used to "
                          "retrieve them.")
            self._iter_subset_class = resolve_iterator_class(
                'shuffled_sequential'
            )

def get_test_array(shape, dtype, strides=None, no_zeros=False, high=None):
    shape = wrap_in_tuple(shape)
    dtype = dtypes.normalize_type(dtype)

    if dtype.names is not None:
        result = numpy.empty(shape, dtype)
        for name in dtype.names:
            result[name] = get_test_array(shape, dtype[name], no_zeros=no_zeros, high=high)
    else:
        if dtypes.is_integer(dtype):
            low = 1 if no_zeros else 0
            if high is None:
                high = 100 # will work even with signed chars
            get_arr = lambda: numpy.random.randint(low, high, shape).astype(dtype)
        else:
            low = 0.01 if no_zeros else 0
            if high is None:
                high = 1.0
            get_arr = lambda: numpy.random.uniform(low, high, shape).astype(dtype)

        if dtypes.is_complex(dtype):
            result = get_arr() + 1j * get_arr()
        else:
            result = get_arr()

    if strides is not None:
        result = as_strided(result, result.shape, strides)

    return result

    def __init__(self, label_image=None, connectivity=1, data=None, **attr):

        super(RAG, self).__init__(data, **attr)
        if self.number_of_nodes() == 0:
            self.max_id = 0
        else:
            self.max_id = max(self.nodes_iter())

        if label_image is not None:
            fp = ndi.generate_binary_structure(label_image.ndim, connectivity)
            # In the next ``ndi.generic_filter`` function, the kwarg
            # ``output`` is used to provide a strided array with a single
            # 64-bit floating point number, to which the function repeatedly
            # writes. This is done because even if we don't care about the
            # output, without this, a float array of the same shape as the
            # input image will be created and that could be expensive in
            # memory consumption.
            ndi.generic_filter(
                label_image,
                function=_add_edge_filter,
                footprint=fp,
                mode='nearest',
                output=as_strided(np.empty((1,), dtype=np.float_),
                                  shape=label_image.shape,
                                  strides=((0,) * label_image.ndim)),
                extra_arguments=(self,))

def semicast(*arrays):
    """
    Broadcast compatible ndarray shape prefixes.
    """

    # establish the final prefix shape
    pre_ndim    = max(len(a.shape[:i]) for (a, i) in arrays)
    pre_padding = [(1,) * (pre_ndim - len(a.shape[:i])) for (a, i) in arrays]
    pre_shape   = tuple(map(max, *(p + a.shape[:i] for ((a, i), p) in zip(arrays, pre_padding))))

    # broadcast the arrays
    from numpy.lib.stride_tricks import as_strided

    casts = []

    for ((a, i), p) in zip(arrays, pre_padding):
        if i is None:
            i = len(a.shape)

        for (c, d) in zip(pre_shape[len(p):], a.shape[:i]):
            if c != d and d != 1:
                raise ValueError("array shapes incompatible for semicast")

        strides  = (0,) * len(p) + tuple(0 if d == 1 else s for (d, s) in zip(a.shape, a.strides))
        casts   += [as_strided(a, pre_shape + a.shape[i:], strides)]

    # repair dtypes (broken by as_strided)
    for ((a, _), cast) in zip(arrays, casts):
        cast.dtype = a.dtype

    # done
    return (pre_shape, casts)

def halfoverlap(X, window_size):
    """
    Create an overlapped version of X using 50% of window_size as overlap.

    Parameters
    ----------
    X : ndarray, shape=(n_samples,)
        Input signal to window and overlap

    window_size : int
        Size of windows to take

    Returns
    -------
    X_strided : shape=(n_windows, window_size)
        2D array of overlapped X
    """
    if window_size % 2 != 0:
        raise ValueError("Window size must be even!")
    window_step = window_size // 2
    # Make sure there are an even number of windows before stridetricks
    append = np.zeros((window_size - len(X) % window_size))
    X = np.hstack((X, append))
    num_frames = len(X) // window_step - 1
    row_stride = X.itemsize * window_step
    col_stride = X.itemsize
    X_strided = as_strided(X, shape=(num_frames, window_size),
                           strides=(row_stride, col_stride))
    return X_strided

def cross_correlation(x, y, maxlag):
    """
    Cross correlation with a maximum number of lags.

    `x` and `y` must be one-dimensional numpy arrays with the same length.

    This computes the same result as
        numpy.correlate(x, y, mode='full')[len(a)-maxlag-1:len(a)+maxlag]

    The return vaue has length 2*maxlag + 1.

    Author: http://stackoverflow.com/questions/30677241
            Warren Weckesser
    """
    from numpy.lib.stride_tricks import as_strided

    def _check_arg(x, xname):
        x = np.asarray(x)
        if x.ndim != 1:
            raise ValueError('%s must be one-dimensional.' % xname)
        return x

    x = _check_arg(x, 'x')
    y = _check_arg(y, 'y')
    py = np.pad(y.conj(), 2*maxlag, mode='constant')
    T = as_strided(py[2*maxlag:], shape=(2*maxlag+1, len(y) + 2*maxlag),
                   strides=(-py.strides[0], py.strides[0]))
    px = np.pad(x, maxlag, mode='constant')
    return T.dot(px)

 def __init__(self, *shape):
     if len(shape) == 1 and isinstance(shape[0], tuple):
         shape = shape[0]
     x = as_strided(_nx.zeros(1), shape=shape,
                    strides=_nx.zeros_like(shape))
     self._it = _nx.nditer(x, flags=['multi_index', 'zerosize_ok'],
                           order='C')

    def _fast_synthesize(frequency):
        """A faster way to synthesize a signal.
            Generate one cycle, and simulate arbitrary repetitions
            using array indexing tricks.
        """
        # hack so that we can ensure an integer number of periods and samples
        # rounds frequency to 1st decimal, s.t. 10 * frequency will be an int
        frequency = np.round(frequency, n_dec)

        # Generate 10*frequency periods at this frequency
        # Equivalent to n_samples = int(n_periods * fs / frequency)
        # n_periods = 10*frequency is the smallest integer that guarantees
        # that n_samples will be an integer, since assuming 10*frequency
        # is an integer
        n_samples = int(10.0**n_dec * fs)

        short_signal = function(2.0 * np.pi * np.arange(n_samples) *
                                frequency / fs)

        # Calculate the number of loops we need to fill the duration
        n_repeats = int(np.ceil(length/float(short_signal.shape[0])))

        # Simulate tiling the short buffer by using stride tricks
        long_signal = as_strided(short_signal,
                                 shape=(n_repeats, len(short_signal)),
                                 strides=(0, short_signal.itemsize))

        # Use a flatiter to simulate a long 1D buffer
        return long_signal.flat

def convolve(A, B, axes=None, dot_axes=[(),()], mode='full'):
    assert mode in ['valid', 'full'], "Mode {0} not yet implemented".format(mode)
    if axes is None:
        axes = [list(range(A.ndim)), list(range(A.ndim))]
    wrong_order = any([B.shape[ax_B] < A.shape[ax_A] for ax_A, ax_B in zip(*axes)])
    if wrong_order:
        if mode=='valid' and not all([B.shape[ax_B] <= A.shape[ax_A] for ax_A, ax_B in zip(*axes)]):
                raise Exception("One array must be larger than the other along all convolved dimensions")
        elif mode != 'full' or B.size <= A.size: # Tie breaker
            i1 =      B.ndim - len(dot_axes[1]) - len(axes[1]) # B ignore
            i2 = i1 + A.ndim - len(dot_axes[0]) - len(axes[0]) # A ignore
            i3 = i2 + len(axes[0])
            ignore_B = list(range(i1))
            ignore_A = list(range(i1, i2))
            conv     = list(range(i2, i3))
            return convolve(B, A, axes=axes[::-1], dot_axes=dot_axes[::-1], mode=mode).transpose(ignore_A + ignore_B + conv)

    if mode == 'full':
        B = pad_to_full(B, A, axes[::-1])
    B_view_shape = list(B.shape)
    B_view_strides = list(B.strides)
    flipped_idxs = [slice(None)] * A.ndim
    for ax_A, ax_B in zip(*axes):
        B_view_shape.append(abs(B.shape[ax_B] - A.shape[ax_A]) + 1)
        B_view_strides.append(B.strides[ax_B])
        B_view_shape[ax_B] = A.shape[ax_A]
        flipped_idxs[ax_A] = slice(None, None, -1)

    B_view = as_strided(B, B_view_shape, B_view_strides)
    A_view = A[flipped_idxs]
    all_axes = [list(axes[i]) + list(dot_axes[i]) for i in [0, 1]]
    return einsum_tensordot(A_view, B_view, all_axes)

def test_writeable():
    # broadcast_to should return a readonly array
    original = np.array([1, 2, 3])
    result = broadcast_to(original, (2, 3))
    assert_equal(result.flags.writeable, False)
    assert_raises(ValueError, result.__setitem__, slice(None), 0)

    # but the result of broadcast_arrays needs to be writeable (for now), to
    # preserve backwards compatibility
    for results in [broadcast_arrays(original),
                    broadcast_arrays(0, original)]:
        for result in results:
            assert_equal(result.flags.writeable, True)
    # keep readonly input readonly
    original.flags.writeable = False
    _, result = broadcast_arrays(0, original)
    assert_equal(result.flags.writeable, False)

    # regression test for GH6491
    shape = (2,)
    strides = [0]
    tricky_array = as_strided(np.array(0), shape, strides)
    other = np.zeros((1,))
    first, second = broadcast_arrays(tricky_array, other)
    assert_(first.shape == second.shape)