Python numpy.asanyarray() Examples

def image_copy(img, dataobj=Ellipsis, affine=Ellipsis, image_type=None):
    '''
    image_copy(image) copies the given image and returns the new object; the affine, header, and
      dataobj members are duplicated so that changes will not change the original image.
    
    The following optional arguments may be given to overwrites part of the new image; in each case,
    the default value (Ellipsis) specifies that no update should be made.
      * dataobj (default: Ellipsis) may overwrite the new image's dataobj object.
      * affine (default: Ellipsis) may overwrite the new image's affine transformation matrix.
    '''
    dataobj = np.array(img.dataobj) if dataobj is Ellipsis else np.asanyarray(dataobj)
    imspec = to_image_spec(img)
    imtype = to_image_type(img if image_type is None else image_type)
    affine = imspec['affine'] if affine is Ellipsis else affine
    imspec['affine'] = affine
    return imtype.create(dataobj, imspec) 

def apply_affine(aff, coords):
    '''
    apply_affine(affine, coords) yields the result of applying the given affine transformation to
      the given coordinate or coordinates.

    This function expects coords to be a (dims X n) matrix but if the first dimension is neither 2
    nor 3, coords.T is used; i.e.:
      apply_affine(affine3x3, coords2xN) ==> newcoords2xN
      apply_affine(affine4x4, coords3xN) ==> newcoords3xN
      apply_affine(affine3x3, coordsNx2) ==> newcoordsNx2 (for N != 2)
      apply_affine(affine4x4, coordsNx3) ==> newcoordsNx3 (for N != 3)
    '''
    if aff is None: return coords
    (coords,tr) = (np.asanyarray(coords), False)
    if len(coords.shape) == 1: return np.squeeze(apply_affine(np.reshape(coords, (-1,1)), aff))
    elif len(coords.shape) > 2: raise ValueError('cannot apply affine to ND-array for N > 2')
    if   len(coords) == 2: aff = to_affine(aff, 2)
    elif len(coords) == 3: aff = to_affine(aff, 3)
    else: (coords,aff,tr) = (coords.T, to_affine(aff, coords.shape[1]), True)
    r = np.dot(aff, np.vstack([coords, np.ones([1,coords.shape[1]])]))[:-1]
    return r.T if tr else r 

def supercell_space_required(transform_oo, transform_vv, final_space):
    """
    For a given orbital transformation and a given `ov` mask in the transformed space, calculates a minimal `ov` mask
    in the original space required to achieve this transform.
    Args:
        transform_oo (ndarray): the transformation in the occupied space;
        transform_vv (ndarray): the transformation in the virtual space;
        final_space (ndarray): the final `ov` space. Basis order: [k_o, o, k_v, v];

    Returns:
        The initial active space. Basis order: [k_o, o, k_v, v].
    """
    final_space = numpy.asanyarray(final_space)
    final_space = final_space.reshape(final_space.shape[:-1] + (transform_oo.shape[1], transform_vv.shape[1]))
    result = einsum(
        "ao,bv,...ov->...ab",
        (transform_oo.toarray() != 0).astype(int),
        (transform_vv.toarray() != 0).astype(int),
        final_space.astype(int),
    ) != 0
    return result.reshape(result.shape[:-2] + (-1,)) 

def orb2ov(space, nocc, nmo):
    """
    Converts orbital active space specification into ov-pairs space spec.
    Args:
        space (ndarray): the obital space. Basis order: [k, orb=o+v];
        nocc (Iterable): the numbers of occupied orbitals per k-point;
        nmo (Iterable): the total numbers of orbitals per k-point;

    Returns:
        The ov space specification. Basis order: [k_o, o, k_v, v].
    """
    space = numpy.asanyarray(space)
    m = ko_mask(nocc, nmo)  # [k, orb=o+v]
    o = space[...,  m]  # [k, o]
    v = space[..., ~m]  # [k, v]
    return (o[..., numpy.newaxis] * v[..., numpy.newaxis, :]).reshape(space.shape[:-1] + (-1,))  # [k_o, o, k_v, v] 

def vector_to_amplitudes(vectors, nocc, nmo):
    """
    Transforms (reshapes) and normalizes vectors into amplitudes.
    Args:
        vectors (numpy.ndarray): raw eigenvectors to transform;
        nocc (tuple): numbers of occupied orbitals;
        nmo (int): the total number of orbitals per k-point;

    Returns:
        Amplitudes with the following shape: (# of roots, 2 (x or y), # of kpts, # of occupied orbitals,
        # of virtual orbitals).
    """
    if not all(i == nocc[0] for i in nocc):
        raise NotImplementedError("Varying occupation numbers are not implemented yet")
    nk = len(nocc)
    nocc = nocc[0]
    if not all(i == nmo[0] for i in nmo):
        raise NotImplementedError("Varying AO spaces are not implemented yet")
    nmo = nmo[0]
    vectors = numpy.asanyarray(vectors)
    vectors = vectors.reshape(2, nk, nocc, nmo-nocc, vectors.shape[1])
    norm = (abs(vectors) ** 2).sum(axis=(1, 2, 3))
    norm = 2 * (norm[0] - norm[1])
    vectors /= norm ** .5
    return vectors.transpose(4, 0, 1, 2, 3) 

def vector_to_amplitudes(vectors, nocc, nmo):
    """
    Transforms (reshapes) and normalizes vectors into amplitudes.
    Args:
        vectors (numpy.ndarray): raw eigenvectors to transform;
        nocc (tuple): numbers of occupied orbitals;
        nmo (tuple): the total numbers of AOs per k-point;

    Returns:
        Amplitudes with the following shape: (# of roots, 2 (x or y), # of kpts, # of kpts, # of occupied orbitals,
        # of virtual orbitals).
    """
    if not all(i == nocc[0] for i in nocc):
        raise NotImplementedError("Varying occupation numbers are not implemented yet")
    nk = len(nocc)
    nocc = nocc[0]
    if not all(i == nmo[0] for i in nmo):
        raise NotImplementedError("Varying AO spaces are not implemented yet")
    nmo = nmo[0]
    vectors = numpy.asanyarray(vectors)
    vectors = vectors.reshape(2, nk, nk, nocc, nmo-nocc, vectors.shape[1])
    norm = (abs(vectors) ** 2).sum(axis=(1, 2, 3, 4))
    norm = 2 * (norm[0] - norm[1])
    vectors /= norm ** .5
    return vectors.transpose(5, 0, 1, 2, 3, 4) 

def vector_to_amplitudes(vectors, nocc, nmo):
    """
    Transforms (reshapes) and normalizes vectors into amplitudes.
    Args:
        vectors (numpy.ndarray): raw eigenvectors to transform;
        nocc (int): number of occupied orbitals;
        nmo (int): the total number of orbitals;

    Returns:
        Amplitudes with the following shape: (# of roots, 2 (x or y), # of occupied orbitals, # of virtual orbitals).
    """
    vectors = numpy.asanyarray(vectors)
    vectors = vectors.reshape(2, nocc, nmo-nocc, vectors.shape[1])
    norm = (abs(vectors) ** 2).sum(axis=(1, 2))
    norm = 2 * (norm[0] - norm[1])
    vectors /= norm ** .5
    return vectors.transpose(3, 0, 1, 2) 

def format_mask(x):
    """
    Formats a mask into a readable string.
    Args:
        x (ndarray): an array with the mask;

    Returns:
        A readable string with the mask.
    """
    x = numpy.asanyarray(x)
    if len(x) == 0:
        return "(empty)"
    if x.dtype == bool:
        x = numpy.argwhere(x)[:, 0]
    grps = tuple(list(g) for _, g in groupby(x, lambda n, c=count(): n-next(c)))
    return ",".join("{:d}-{:d}".format(i[0], i[-1]) if len(i) > 1 else "{:d}".format(i[0]) for i in grps) 

def nside_to_pixel_area(nside):
    """
    Find the area of HEALPix pixels given the pixel dimensions of one of
    the 12 'top-level' HEALPix tiles.

    Parameters
    ----------
    nside : int
        The number of pixels on the side of one of the 12 'top-level' HEALPix tiles.

    Returns
    -------
    pixel_area : :class:`~astropy.units.Quantity`
        The area of the HEALPix pixels
    """
    nside = np.asanyarray(nside, dtype=np.int64)
    _validate_nside(nside)
    npix = 12 * nside * nside
    pixel_area = 4 * math.pi / npix * u.sr
    return pixel_area 

def nside_to_pixel_resolution(nside):
    """
    Find the resolution of HEALPix pixels given the pixel dimensions of one of
    the 12 'top-level' HEALPix tiles.

    Parameters
    ----------
    nside : int
        The number of pixels on the side of one of the 12 'top-level' HEALPix tiles.

    Returns
    -------
    resolution : :class:`~astropy.units.Quantity`
        The resolution of the HEALPix pixels

    See also
    --------
    pixel_resolution_to_nside
    """
    nside = np.asanyarray(nside, dtype=np.int64)
    _validate_nside(nside)
    return (nside_to_pixel_area(nside) ** 0.5).to(u.arcmin) 

def nside_to_npix(nside):
    """
    Find the number of pixels corresponding to a HEALPix resolution.

    Parameters
    ----------
    nside : int
        The number of pixels on the side of one of the 12 'top-level' HEALPix tiles.

    Returns
    -------
    npix : int
        The number of pixels in the HEALPix map.
    """
    nside = np.asanyarray(nside, dtype=np.int64)
    _validate_nside(nside)
    return 12 * nside ** 2 

def take(self, indexer, axis=1, verify=True, convert=True):
        """
        Take items along any axis.
        """
        self._consolidate_inplace()
        indexer = (np.arange(indexer.start, indexer.stop, indexer.step,
                             dtype='int64')
                   if isinstance(indexer, slice)
                   else np.asanyarray(indexer, dtype='int64'))

        n = self.shape[axis]
        if convert:
            indexer = maybe_convert_indices(indexer, n)

        if verify:
            if ((indexer == -1) | (indexer >= n)).any():
                raise Exception('Indices must be nonzero and less than '
                                'the axis length')

        new_labels = self.axes[axis].take(indexer)
        return self.reindex_indexer(new_axis=new_labels, indexer=indexer,
                                    axis=axis, allow_dups=True) 

def apply_lut(in_dseg, lut, newpath=None):
    """Map the input discrete segmentation to a new label set (lookup table, LUT)."""
    import numpy as np
    import nibabel as nb
    from nipype.utils.filemanip import fname_presuffix

    if newpath is None:
        from os import getcwd
        newpath = getcwd()

    out_file = fname_presuffix(in_dseg, suffix='_dseg', newpath=newpath)
    lut = np.array(lut, dtype='int16')

    segm = nb.load(in_dseg)
    hdr = segm.header.copy()
    hdr.set_data_dtype('int16')
    segm.__class__(lut[np.asanyarray(segm.dataobj, dtype=int)].astype('int16'),
                   segm.affine, hdr).to_filename(out_file)

    return out_file 

def _call_series(self, a, inputs):
        if isinstance(self._q, TENSOR_TYPE):
            q_val = self._q
            index_val = pd.Index([], dtype=q_val.dtype)
            store_index_value = False
        else:
            q_val = np.asanyarray(self._q)
            index_val = pd.Index(q_val)
            store_index_value = True

        # get dtype by tensor
        a_t = astensor(a)
        self._dtype = dtype = tensor_quantile(
            a_t, self._q, interpolation=self._interpolation,
            handle_non_numeric=not self._numeric_only).dtype

        if q_val.ndim == 0:
            return self.new_scalar(inputs, dtype=dtype)
        else:
            return self.new_series(
                inputs, shape=q_val.shape, dtype=dtype,
                index_value=parse_index(index_val, a, q_val, self._interpolation,
                                        type(self).__name__, store_data=store_index_value),
                name=a.name) 

def cifti_split(cii, label=('lh', 'rh', 'rest'), subject=None, hemi=None, null=np.nan):
    '''
    cifti_split(cii, label) yields the rows or columns of the given cifti file that correspond to
      the given label (see below).
    cifti_split(cii) is equivalent to cifti_split(cii, ('lh', 'rh', 'rest')).

    The label argument may be any of the following:
      * a valid CIFTI label name such as 'CIFTI_STRUCTURE_CEREBELLUM' or
        'CIFTI_STRUCTURE_CORTEX_LEFT';
      * an abbreviated name such as 'cerebellum' for 'CIFTI_STRUCTURE_CEREBELLUM'.
      * the abbreviations 'lh' and 'rh' which stand for 'CIFTI_STRUCTURE_CORTEX_LEFT' and 
        'CIFTI_STRUCTURE_CORTEX_RIGHT';
      * the special keyword 'rest', which represents all the rows/columns not collected by any other
        instruction ('rest', by itself, results in the whole matrix being returned); or
      * A tuple of the above, indicating that each of the items listed should be returned
        sequentially in a tuple.

    The following optional arguments may be given:
      * subject (default: None) may specify the subject
      * hemi (default: None) can specify the hemisphere object that 
    '''
    dat = np.asanyarray(cii.dataobj if is_image(cii) else cii)
    n = dat.shape[-1]
    atlas = cifti_split._size_data.get(n, None)
    if atlas is None: raise ValueError('cannot split cifti with size %d' % n)
    if atlas not in cifti_split._atlas_cache:
        patt = os.path.join('data', 'fs_LR', '%s.atlasroi.%dk_fs_LR.shape.gii')
        lgii = nib.load(os.path.join(library_path(), patt % ('lh', atlas)))
        rgii = nib.load(os.path.join(library_path(), patt % ('rh', atlas)))
        cifti_split._atlas_cache[atlas] = tuple([pimms.imm_array(gii.darrays[0].data.astype('bool'))
                                                 for gii in (lgii, rgii)])
    (lroi,rroi) = cifti_split._atlas_cache[atlas]
    (ln,lN) = (np.sum(lroi), len(lroi))
    (rn,rN) = (np.sum(rroi), len(rroi))
    (ldat,rdat,sdat) = [np.full(dat.shape[:-1] + (k,), null) for k in [lN, rN, n - ln - rn]]
    ldat[..., lroi] = dat[..., :ln]
    rdat[..., rroi] = dat[..., ln:(ln+rn)]
    sdat[...] = dat[..., (ln+rn):]
    if ln + rn >= n: sdat = None
    return (ldat, rdat, sdat) 

def nan_compare(f, x, y, nan_nan=False, nan_val=False, val_nan=False):
    '''
    nan_compare(f, x, y) is equivalent to f(x, y), which is assumed to be a boolean function that
      broadcasts over x and y (such as numpy.less), except that NaN values in either x or y result
      in a value of False instead of being run through f.

    The argument f must be a numpy comparison function such as numpy.less that accepts the optional
    arguments where and out.

    The following optional arguments may be provided:
      * nan_nan (default: False) specifies the return value (True or False) for comparisons
        equivalent to f(nan, nan).
      * nan_val (default: False) specifies the return value (True or False) for comparisons
        equivalent to f(nan, non_nan).
      * val_nan (default: False) specifies the return value (True or False) for comparisons
        equivalent to f(non_nan, nan).
    '''
    #TODO: This should work with sparse matrices as well
    x = np.asanyarray(x)
    y = np.asanyarray(y)
    xii = np.isnan(x)
    yii = np.isnan(y)
    if not xii.any() and not yii.any(): return f(x, y)
    ii  = (~xii) & (~yii)
    out = np.zeros(ii.shape, dtype=np.bool)
    if nan_nan == nan_val and nan_val == val_nan:
        # All the nan-result values are the same; we can simplify a little...
        if nan_nan: out[~ii] = nan_nan
    else:
        if nan_nan: out[   xii &    yii] = nan_nan
        if nan_val: out[   xii & (~yii)] = nan_val
        if val_nan: out[(~xii) &    yii] = val_nan
    return f(x, y, out=out, where=ii) 

def ov2orb(space, nocc, nmo):
    """
    Converts ov-pairs active space specification into orbital space spec.
    Args:
        space (ndarray): the ov space. Basis order: [k_o, o, k_v, v];
        nocc (Iterable): the numbers of occupied orbitals per k-point;
        nmo (Iterable): the total numbers of orbitals per k-point;

    Returns:
        The orbital space specification. Basis order: [k, orb=o+v].
    """
    nocc = numpy.asanyarray(nocc)
    nmo = numpy.asanyarray(nmo)
    nvirt = nmo - nocc

    space = numpy.asanyarray(space)
    space = space.reshape(space.shape[:-1] + (sum(nocc), sum(nvirt)))  # [k_o, o; k_v, v]

    s_o = numpy.any(space, axis=-1)  # [k_o, o]
    s_v = numpy.any(space, axis=-2)  # [k_v, v]

    o_offset = numpy.cumsum(numpy.concatenate(([0], nocc)))
    v_offset = numpy.cumsum(numpy.concatenate(([0], nvirt)))
    result = []
    for o_fr, o_to, v_fr, v_to in zip(o_offset[:-1], o_offset[1:], v_offset[:-1], v_offset[1:]):
        result.append(s_o[..., o_fr:o_to])
        result.append(s_v[..., v_fr:v_to])
    return numpy.concatenate(result, axis=-1)  # [k, orb=o+v] 

def vector_to_amplitudes(vectors, nocc, nmo):
    """
    Transforms (reshapes) and normalizes vectors into amplitudes.
    Args:
        vectors (numpy.ndarray): raw eigenvectors to transform;
        nocc (tuple): numbers of occupied orbitals;
        nmo (int): the total number of orbitals per k-point;

    Returns:
        Amplitudes with the following shape: (# of roots, 2 (x or y), # of kpts, # of kpts, # of occupied orbitals,
        # of virtual orbitals).
    """
    if not all(i == nocc[0] for i in nocc):
        raise NotImplementedError("Varying occupation numbers are not implemented yet")
    nk = len(nocc)
    nocc = nocc[0]
    if not all(i == nmo[0] for i in nmo):
        raise NotImplementedError("Varying AO spaces are not implemented yet")
    nmo = nmo[0]
    vectors = numpy.asanyarray(vectors)
    # Compared to krhf_slow_supercell, only one k-point index is present here. The second index corresponds to
    # momentum transfer and is integrated out
    vectors = vectors.reshape(2, nk, nocc, nmo-nocc, vectors.shape[1])
    norm = (abs(vectors) ** 2).sum(axis=(1, 2, 3))
    norm = 2 * (norm[0] - norm[1])
    vectors /= norm ** .5
    return vectors.transpose(4, 0, 1, 2, 3) 

def __init__(self, mf, frozen=None):
        """
        Performs TD calculation. Roots and eigenvectors are stored in `self.e`, `self.xy`.
        Args:
            mf: the mean-field model;
            frozen (int, Iterable): the number of frozen valence orbitals or the list of frozen orbitals;
        """
        self._scf = mf
        self.driver = None
        self.nroots = None
        self.eri = None
        self.xy = None
        self.e = None
        self.frozen = frozen
        self.fast = not numpy.iscomplexobj(numpy.asanyarray(mf.mo_coeff)) 

def flip_axis(arr, axis=0):
    ''' Flip contents of `axis` in array `arr`

    ``flip_axis`` is the same transform as ``np.flipud``, but for any
    axis.  For example ``flip_axis(arr, axis=0)`` is the same transform
    as ``np.flipud(arr)``, and ``flip_axis(arr, axis=1)`` is the same
    transform as ``np.fliplr(arr)``

    Parameters
    ----------
    arr : array-like
    axis : int, optional
       axis to flip.  Default `axis` == 0

    Returns
    -------
    farr : array
       Array with axis `axis` flipped

    Examples
    --------
    >>> a = np.arange(6).reshape((2,3))
    >>> a
    array([[0, 1, 2],
           [3, 4, 5]])
    >>> flip_axis(a, axis=0)
    array([[3, 4, 5],
           [0, 1, 2]])
    >>> flip_axis(a, axis=1)
    array([[2, 1, 0],
           [5, 4, 3]])
    '''
    arr = np.asanyarray(arr)
    arr = arr.swapaxes(0, axis)
    arr = np.flipud(arr)
    return arr.swapaxes(axis, 0) 

def get_data(self):
        """returns scaled data for all frames in a numpy array
        returns as a 4D array """
        if self._data is None:
            raise ImageDataError('No data in this image')
        return np.asanyarray(self._data) 

def get_data(self):
        return np.asanyarray(self._data) 

def __setitem__(self, item, data):
        if item == () or item == slice(None):
            words_slice = data_slice = slice(None)
        else:
            words_slice, data_slice = self._item_to_slices(item)

        data = np.asanyarray(data)
        # Check if the new data spans an entire word and is correctly shaped.
        # If so, skip decoding.  If not, decode appropriate words and insert
        # new data.
        if not (data_slice == slice(None)
                and data.shape[-len(self.sample_shape):] == self.sample_shape
                and data.dtype.kind == self.dtype.kind):
            decoder = self._decoders[self._coder]
            current_data = decoder(self.words[words_slice])
            if self.complex_data:
                current_data = current_data.view(self.dtype)
            current_data.shape = (-1,) + self.sample_shape
            current_data[data_slice] = data
            data = current_data

        if data.dtype.kind == 'c':
            data = data.view((data.real.dtype, (2,)))

        encoder = self._encoders[self._coder]
        self.words[words_slice] = encoder(data).ravel().view(self._dtype_word) 

def __setitem__(self, item, data):
        if isinstance(item, str):
            # Headers behave as dictionaries; except for thread_id and
            # invalid_data, assume set base properties all to the same value.
            data = np.broadcast_to(data, (len(self.frames),))
            if item == 'thread_id':
                if len(np.unique(data)) != len(self.frames):
                    raise ValueError("all thread ids should be unique.")
            elif (item != 'invalid_data'
                  and item in VDIFBaseHeader._header_parser.keys()):
                if data.strides != (0,) and len(np.unique(data)) > 1:
                    raise ValueError("base header keys should be identical.")

            for f, value in zip(self.frames, data):
                f.header[item] = value
            return

        (frames, frame_item,
         single_sample, single_frame, single_channel) = self._get_frames(item)

        if single_frame:
            frames[0][frame_item] = data
            return

        data = np.asanyarray(data)
        if single_channel:
            if single_sample or data.ndim <= 1:
                swapped = np.broadcast_to(data, (len(frames),))
            else:
                new_shape = (data.shape[0], len(frames))
                swapped = np.broadcast_to(data, new_shape).swapaxes(0, 1)
        else:
            if single_sample or data.ndim <= 2:
                new_shape = (len(frames),) + data.shape[1:]
                swapped = np.broadcast_to(data, new_shape)
            else:
                new_shape = (data.shape[0], len(frames)) + data.shape[2:]
                swapped = np.broadcast_to(data, new_shape).swapaxes(0, 1)

        for frame, frame_data in zip(frames, swapped):
            frame[frame_item] = frame_data 

def __setitem__(self, item, value):
        if isinstance(item, str):
            return self.header.__setitem__(item, value)

        # Normally, we would just pass on to the payload here, but for
        # Mark 4, we need to deal with data overwritten by the header.
        data = np.asanyarray(value)
        assert data.ndim <= 2

        (payload_item, sample_index, data_shape,
         ninvalid) = self._get_payload_item(item)

        if payload_item is None:
            return

        if ninvalid > 0:
            # See if data has enough dimensions so that we need to remove
            # the part that cannot set anything in the payload.
            if sample_index == ():
                sample_ndim = len(self.sample_shape)
            else:
                sample_ndim = np.empty(self.sample_shape)[sample_index].ndim
            if data.ndim == 1 + sample_ndim:
                data = data[ninvalid:]

        if sample_index != ():
            payload_item = (payload_item,) + sample_index

        self.payload[payload_item] = data 

def npix_to_nside(npix):
    """
    Find the number of pixels on the side of one of the 12 'top-level' HEALPix
    tiles given a total number of pixels.

    Parameters
    ----------
    npix : int
        The number of pixels in the HEALPix map.

    Returns
    -------
    nside : int
        The number of pixels on the side of one of the 12 'top-level' HEALPix tiles.
    """

    npix = np.asanyarray(npix, dtype=np.int64)

    if not np.all(npix % 12 == 0):
        raise ValueError('Number of pixels must be divisible by 12')

    square_root = np.sqrt(npix / 12)
    if not np.all(square_root ** 2 == npix / 12):
        raise ValueError('Number of pixels is not of the form 12 * nside ** 2')

    return np.round(square_root).astype(int) 

def _unique1d(ar, return_index=False, return_inverse=False,
              return_counts=False):
    """
    Find the unique elements of an array, ignoring shape.
    """
    ar = np.asanyarray(ar).flatten()

    optional_indices = return_index or return_inverse

    if optional_indices:
        perm = ar.argsort(kind='mergesort' if return_index else 'quicksort')
        aux = ar[perm]
    else:
        ar.sort()
        aux = ar
    mask = np.empty(aux.shape, dtype=np.bool_)
    mask[:1] = True
    mask[1:] = aux[1:] != aux[:-1]

    ret = (aux[mask],)
    if return_index:
        ret += (perm[mask],)
    if return_inverse:
        imask = np.cumsum(mask) - 1
        inv_idx = np.empty(mask.shape, dtype=np.intp)
        inv_idx[perm] = imask
        ret += (inv_idx,)
    if return_counts:
        idx = np.concatenate(np.nonzero(mask) + ([mask.size],))
        ret += (np.diff(idx),)
    return ret 

def __new__(cls,arr,info={}):
        x = np.asanyarray(arr).view(cls)
        x.info = info.copy()
        return x 

def test_subclasspreservation(self):
        # Checks that masked_array(...,subok=True) preserves the class.
        x = np.arange(5)
        m = [0, 0, 1, 0, 0]
        xinfo = [(i, j) for (i, j) in zip(x, m)]
        xsub = MSubArray(x, mask=m, info={'xsub':xinfo})
        #
        mxsub = masked_array(xsub, subok=False)
        assert_(not isinstance(mxsub, MSubArray))
        assert_(isinstance(mxsub, MaskedArray))
        assert_equal(mxsub._mask, m)
        #
        mxsub = asarray(xsub)
        assert_(not isinstance(mxsub, MSubArray))
        assert_(isinstance(mxsub, MaskedArray))
        assert_equal(mxsub._mask, m)
        #
        mxsub = masked_array(xsub, subok=True)
        assert_(isinstance(mxsub, MSubArray))
        assert_equal(mxsub.info, xsub.info)
        assert_equal(mxsub._mask, xsub._mask)
        #
        mxsub = asanyarray(xsub)
        assert_(isinstance(mxsub, MSubArray))
        assert_equal(mxsub.info, xsub.info)
        assert_equal(mxsub._mask, m) 

def ediff1d(arr, to_end=None, to_begin=None):
    """
    Compute the differences between consecutive elements of an array.

    This function is the equivalent of `numpy.ediff1d` that takes masked
    values into account, see `numpy.ediff1d` for details.

    See Also
    --------
    numpy.ediff1d : Equivalent function for ndarrays.

    """
    arr = ma.asanyarray(arr).flat
    ed = arr[1:] - arr[:-1]
    arrays = [ed]
    #
    if to_begin is not None:
        arrays.insert(0, to_begin)
    if to_end is not None:
        arrays.append(to_end)
    #
    if len(arrays) != 1:
        # We'll save ourselves a copy of a potentially large array in the common
        # case where neither to_begin or to_end was given.
        ed = hstack(arrays)
    #
    return ed