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