1""" Utilities to manipulate numpy arrays """
2
3import numpy as np
4from nibabel.volumeutils import endian_codes, native_code
5
6
7def as_native_array(arr):
8    """ Return `arr` as native byteordered array
9
10    If arr is already native byte ordered, return unchanged.  If it is opposite
11    endian, then make a native byte ordered copy and return that
12
13    Parameters
14    ----------
15    arr : ndarray
16
17    Returns
18    -------
19    native_arr : ndarray
20        If `arr` was native order, this is just `arr`. Otherwise it's a new
21        array such that ``np.all(native_arr == arr)``, with native byte
22        ordering.
23    """
24    if endian_codes[arr.dtype.byteorder] == native_code:
25        return arr
26    return arr.byteswap().newbyteorder()
27
28
29def pinv(a, rcond=1e-15):
30    """Vectorized version of `numpy.linalg.pinv`
31
32    If numpy version is less than 1.8, it falls back to iterating over
33    `np.linalg.pinv` since there isn't a vectorized version of `np.linalg.svd`
34    available.
35
36    Parameters
37    ----------
38    a : array_like (..., M, N)
39        Matrix to be pseudo-inverted.
40    rcond : float
41        Cutoff for small singular values.
42
43    Returns
44    -------
45    B : ndarray (..., N, M)
46        The pseudo-inverse of `a`.
47
48    Raises
49    ------
50    LinAlgError
51        If the SVD computation does not converge.
52
53    See Also
54    --------
55    np.linalg.pinv
56    """
57    a = np.asarray(a)
58
59    swap = np.arange(a.ndim)
60    swap[[-2, -1]] = swap[[-1, -2]]
61    u, s, v = np.linalg.svd(a, full_matrices=False)
62    cutoff = np.maximum.reduce(s, axis=-1, keepdims=True) * rcond
63    mask = s > cutoff
64    s[mask] = 1. / s[mask]
65    s[~mask] = 0
66    return np.einsum('...ij,...jk',
67                     np.transpose(v, swap) * s[..., None, :],
68                     np.transpose(u, swap))
69