extras.py 57 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923
  1. """
  2. Masked arrays add-ons.
  3. A collection of utilities for `numpy.ma`.
  4. :author: Pierre Gerard-Marchant
  5. :contact: pierregm_at_uga_dot_edu
  6. :version: $Id: extras.py 3473 2007-10-29 15:18:13Z jarrod.millman $
  7. """
  8. __all__ = [
  9. 'apply_along_axis', 'apply_over_axes', 'atleast_1d', 'atleast_2d',
  10. 'atleast_3d', 'average', 'clump_masked', 'clump_unmasked',
  11. 'column_stack', 'compress_cols', 'compress_nd', 'compress_rowcols',
  12. 'compress_rows', 'count_masked', 'corrcoef', 'cov', 'diagflat', 'dot',
  13. 'dstack', 'ediff1d', 'flatnotmasked_contiguous', 'flatnotmasked_edges',
  14. 'hsplit', 'hstack', 'isin', 'in1d', 'intersect1d', 'mask_cols', 'mask_rowcols',
  15. 'mask_rows', 'masked_all', 'masked_all_like', 'median', 'mr_',
  16. 'notmasked_contiguous', 'notmasked_edges', 'polyfit', 'row_stack',
  17. 'setdiff1d', 'setxor1d', 'stack', 'unique', 'union1d', 'vander', 'vstack',
  18. ]
  19. import itertools
  20. import warnings
  21. from . import core as ma
  22. from .core import (
  23. MaskedArray, MAError, add, array, asarray, concatenate, filled, count,
  24. getmask, getmaskarray, make_mask_descr, masked, masked_array, mask_or,
  25. nomask, ones, sort, zeros, getdata, get_masked_subclass, dot,
  26. mask_rowcols
  27. )
  28. import numpy as np
  29. from numpy import ndarray, array as nxarray
  30. import numpy.core.umath as umath
  31. from numpy.core.multiarray import normalize_axis_index
  32. from numpy.core.numeric import normalize_axis_tuple
  33. from numpy.lib.function_base import _ureduce
  34. from numpy.lib.index_tricks import AxisConcatenator
  35. def issequence(seq):
  36. """
  37. Is seq a sequence (ndarray, list or tuple)?
  38. """
  39. return isinstance(seq, (ndarray, tuple, list))
  40. def count_masked(arr, axis=None):
  41. """
  42. Count the number of masked elements along the given axis.
  43. Parameters
  44. ----------
  45. arr : array_like
  46. An array with (possibly) masked elements.
  47. axis : int, optional
  48. Axis along which to count. If None (default), a flattened
  49. version of the array is used.
  50. Returns
  51. -------
  52. count : int, ndarray
  53. The total number of masked elements (axis=None) or the number
  54. of masked elements along each slice of the given axis.
  55. See Also
  56. --------
  57. MaskedArray.count : Count non-masked elements.
  58. Examples
  59. --------
  60. >>> import numpy.ma as ma
  61. >>> a = np.arange(9).reshape((3,3))
  62. >>> a = ma.array(a)
  63. >>> a[1, 0] = ma.masked
  64. >>> a[1, 2] = ma.masked
  65. >>> a[2, 1] = ma.masked
  66. >>> a
  67. masked_array(
  68. data=[[0, 1, 2],
  69. [--, 4, --],
  70. [6, --, 8]],
  71. mask=[[False, False, False],
  72. [ True, False, True],
  73. [False, True, False]],
  74. fill_value=999999)
  75. >>> ma.count_masked(a)
  76. 3
  77. When the `axis` keyword is used an array is returned.
  78. >>> ma.count_masked(a, axis=0)
  79. array([1, 1, 1])
  80. >>> ma.count_masked(a, axis=1)
  81. array([0, 2, 1])
  82. """
  83. m = getmaskarray(arr)
  84. return m.sum(axis)
  85. def masked_all(shape, dtype=float):
  86. """
  87. Empty masked array with all elements masked.
  88. Return an empty masked array of the given shape and dtype, where all the
  89. data are masked.
  90. Parameters
  91. ----------
  92. shape : tuple
  93. Shape of the required MaskedArray.
  94. dtype : dtype, optional
  95. Data type of the output.
  96. Returns
  97. -------
  98. a : MaskedArray
  99. A masked array with all data masked.
  100. See Also
  101. --------
  102. masked_all_like : Empty masked array modelled on an existing array.
  103. Examples
  104. --------
  105. >>> import numpy.ma as ma
  106. >>> ma.masked_all((3, 3))
  107. masked_array(
  108. data=[[--, --, --],
  109. [--, --, --],
  110. [--, --, --]],
  111. mask=[[ True, True, True],
  112. [ True, True, True],
  113. [ True, True, True]],
  114. fill_value=1e+20,
  115. dtype=float64)
  116. The `dtype` parameter defines the underlying data type.
  117. >>> a = ma.masked_all((3, 3))
  118. >>> a.dtype
  119. dtype('float64')
  120. >>> a = ma.masked_all((3, 3), dtype=np.int32)
  121. >>> a.dtype
  122. dtype('int32')
  123. """
  124. a = masked_array(np.empty(shape, dtype),
  125. mask=np.ones(shape, make_mask_descr(dtype)))
  126. return a
  127. def masked_all_like(arr):
  128. """
  129. Empty masked array with the properties of an existing array.
  130. Return an empty masked array of the same shape and dtype as
  131. the array `arr`, where all the data are masked.
  132. Parameters
  133. ----------
  134. arr : ndarray
  135. An array describing the shape and dtype of the required MaskedArray.
  136. Returns
  137. -------
  138. a : MaskedArray
  139. A masked array with all data masked.
  140. Raises
  141. ------
  142. AttributeError
  143. If `arr` doesn't have a shape attribute (i.e. not an ndarray)
  144. See Also
  145. --------
  146. masked_all : Empty masked array with all elements masked.
  147. Examples
  148. --------
  149. >>> import numpy.ma as ma
  150. >>> arr = np.zeros((2, 3), dtype=np.float32)
  151. >>> arr
  152. array([[0., 0., 0.],
  153. [0., 0., 0.]], dtype=float32)
  154. >>> ma.masked_all_like(arr)
  155. masked_array(
  156. data=[[--, --, --],
  157. [--, --, --]],
  158. mask=[[ True, True, True],
  159. [ True, True, True]],
  160. fill_value=1e+20,
  161. dtype=float32)
  162. The dtype of the masked array matches the dtype of `arr`.
  163. >>> arr.dtype
  164. dtype('float32')
  165. >>> ma.masked_all_like(arr).dtype
  166. dtype('float32')
  167. """
  168. a = np.empty_like(arr).view(MaskedArray)
  169. a._mask = np.ones(a.shape, dtype=make_mask_descr(a.dtype))
  170. return a
  171. #####--------------------------------------------------------------------------
  172. #---- --- Standard functions ---
  173. #####--------------------------------------------------------------------------
  174. class _fromnxfunction:
  175. """
  176. Defines a wrapper to adapt NumPy functions to masked arrays.
  177. An instance of `_fromnxfunction` can be called with the same parameters
  178. as the wrapped NumPy function. The docstring of `newfunc` is adapted from
  179. the wrapped function as well, see `getdoc`.
  180. This class should not be used directly. Instead, one of its extensions that
  181. provides support for a specific type of input should be used.
  182. Parameters
  183. ----------
  184. funcname : str
  185. The name of the function to be adapted. The function should be
  186. in the NumPy namespace (i.e. ``np.funcname``).
  187. """
  188. def __init__(self, funcname):
  189. self.__name__ = funcname
  190. self.__doc__ = self.getdoc()
  191. def getdoc(self):
  192. """
  193. Retrieve the docstring and signature from the function.
  194. The ``__doc__`` attribute of the function is used as the docstring for
  195. the new masked array version of the function. A note on application
  196. of the function to the mask is appended.
  197. Parameters
  198. ----------
  199. None
  200. """
  201. npfunc = getattr(np, self.__name__, None)
  202. doc = getattr(npfunc, '__doc__', None)
  203. if doc:
  204. sig = self.__name__ + ma.get_object_signature(npfunc)
  205. doc = ma.doc_note(doc, "The function is applied to both the _data "
  206. "and the _mask, if any.")
  207. return '\n\n'.join((sig, doc))
  208. return
  209. def __call__(self, *args, **params):
  210. pass
  211. class _fromnxfunction_single(_fromnxfunction):
  212. """
  213. A version of `_fromnxfunction` that is called with a single array
  214. argument followed by auxiliary args that are passed verbatim for
  215. both the data and mask calls.
  216. """
  217. def __call__(self, x, *args, **params):
  218. func = getattr(np, self.__name__)
  219. if isinstance(x, ndarray):
  220. _d = func(x.__array__(), *args, **params)
  221. _m = func(getmaskarray(x), *args, **params)
  222. return masked_array(_d, mask=_m)
  223. else:
  224. _d = func(np.asarray(x), *args, **params)
  225. _m = func(getmaskarray(x), *args, **params)
  226. return masked_array(_d, mask=_m)
  227. class _fromnxfunction_seq(_fromnxfunction):
  228. """
  229. A version of `_fromnxfunction` that is called with a single sequence
  230. of arrays followed by auxiliary args that are passed verbatim for
  231. both the data and mask calls.
  232. """
  233. def __call__(self, x, *args, **params):
  234. func = getattr(np, self.__name__)
  235. _d = func(tuple([np.asarray(a) for a in x]), *args, **params)
  236. _m = func(tuple([getmaskarray(a) for a in x]), *args, **params)
  237. return masked_array(_d, mask=_m)
  238. class _fromnxfunction_args(_fromnxfunction):
  239. """
  240. A version of `_fromnxfunction` that is called with multiple array
  241. arguments. The first non-array-like input marks the beginning of the
  242. arguments that are passed verbatim for both the data and mask calls.
  243. Array arguments are processed independently and the results are
  244. returned in a list. If only one array is found, the return value is
  245. just the processed array instead of a list.
  246. """
  247. def __call__(self, *args, **params):
  248. func = getattr(np, self.__name__)
  249. arrays = []
  250. args = list(args)
  251. while len(args) > 0 and issequence(args[0]):
  252. arrays.append(args.pop(0))
  253. res = []
  254. for x in arrays:
  255. _d = func(np.asarray(x), *args, **params)
  256. _m = func(getmaskarray(x), *args, **params)
  257. res.append(masked_array(_d, mask=_m))
  258. if len(arrays) == 1:
  259. return res[0]
  260. return res
  261. class _fromnxfunction_allargs(_fromnxfunction):
  262. """
  263. A version of `_fromnxfunction` that is called with multiple array
  264. arguments. Similar to `_fromnxfunction_args` except that all args
  265. are converted to arrays even if they are not so already. This makes
  266. it possible to process scalars as 1-D arrays. Only keyword arguments
  267. are passed through verbatim for the data and mask calls. Arrays
  268. arguments are processed independently and the results are returned
  269. in a list. If only one arg is present, the return value is just the
  270. processed array instead of a list.
  271. """
  272. def __call__(self, *args, **params):
  273. func = getattr(np, self.__name__)
  274. res = []
  275. for x in args:
  276. _d = func(np.asarray(x), **params)
  277. _m = func(getmaskarray(x), **params)
  278. res.append(masked_array(_d, mask=_m))
  279. if len(args) == 1:
  280. return res[0]
  281. return res
  282. atleast_1d = _fromnxfunction_allargs('atleast_1d')
  283. atleast_2d = _fromnxfunction_allargs('atleast_2d')
  284. atleast_3d = _fromnxfunction_allargs('atleast_3d')
  285. vstack = row_stack = _fromnxfunction_seq('vstack')
  286. hstack = _fromnxfunction_seq('hstack')
  287. column_stack = _fromnxfunction_seq('column_stack')
  288. dstack = _fromnxfunction_seq('dstack')
  289. stack = _fromnxfunction_seq('stack')
  290. hsplit = _fromnxfunction_single('hsplit')
  291. diagflat = _fromnxfunction_single('diagflat')
  292. #####--------------------------------------------------------------------------
  293. #----
  294. #####--------------------------------------------------------------------------
  295. def flatten_inplace(seq):
  296. """Flatten a sequence in place."""
  297. k = 0
  298. while (k != len(seq)):
  299. while hasattr(seq[k], '__iter__'):
  300. seq[k:(k + 1)] = seq[k]
  301. k += 1
  302. return seq
  303. def apply_along_axis(func1d, axis, arr, *args, **kwargs):
  304. """
  305. (This docstring should be overwritten)
  306. """
  307. arr = array(arr, copy=False, subok=True)
  308. nd = arr.ndim
  309. axis = normalize_axis_index(axis, nd)
  310. ind = [0] * (nd - 1)
  311. i = np.zeros(nd, 'O')
  312. indlist = list(range(nd))
  313. indlist.remove(axis)
  314. i[axis] = slice(None, None)
  315. outshape = np.asarray(arr.shape).take(indlist)
  316. i.put(indlist, ind)
  317. res = func1d(arr[tuple(i.tolist())], *args, **kwargs)
  318. # if res is a number, then we have a smaller output array
  319. asscalar = np.isscalar(res)
  320. if not asscalar:
  321. try:
  322. len(res)
  323. except TypeError:
  324. asscalar = True
  325. # Note: we shouldn't set the dtype of the output from the first result
  326. # so we force the type to object, and build a list of dtypes. We'll
  327. # just take the largest, to avoid some downcasting
  328. dtypes = []
  329. if asscalar:
  330. dtypes.append(np.asarray(res).dtype)
  331. outarr = zeros(outshape, object)
  332. outarr[tuple(ind)] = res
  333. Ntot = np.product(outshape)
  334. k = 1
  335. while k < Ntot:
  336. # increment the index
  337. ind[-1] += 1
  338. n = -1
  339. while (ind[n] >= outshape[n]) and (n > (1 - nd)):
  340. ind[n - 1] += 1
  341. ind[n] = 0
  342. n -= 1
  343. i.put(indlist, ind)
  344. res = func1d(arr[tuple(i.tolist())], *args, **kwargs)
  345. outarr[tuple(ind)] = res
  346. dtypes.append(asarray(res).dtype)
  347. k += 1
  348. else:
  349. res = array(res, copy=False, subok=True)
  350. j = i.copy()
  351. j[axis] = ([slice(None, None)] * res.ndim)
  352. j.put(indlist, ind)
  353. Ntot = np.product(outshape)
  354. holdshape = outshape
  355. outshape = list(arr.shape)
  356. outshape[axis] = res.shape
  357. dtypes.append(asarray(res).dtype)
  358. outshape = flatten_inplace(outshape)
  359. outarr = zeros(outshape, object)
  360. outarr[tuple(flatten_inplace(j.tolist()))] = res
  361. k = 1
  362. while k < Ntot:
  363. # increment the index
  364. ind[-1] += 1
  365. n = -1
  366. while (ind[n] >= holdshape[n]) and (n > (1 - nd)):
  367. ind[n - 1] += 1
  368. ind[n] = 0
  369. n -= 1
  370. i.put(indlist, ind)
  371. j.put(indlist, ind)
  372. res = func1d(arr[tuple(i.tolist())], *args, **kwargs)
  373. outarr[tuple(flatten_inplace(j.tolist()))] = res
  374. dtypes.append(asarray(res).dtype)
  375. k += 1
  376. max_dtypes = np.dtype(np.asarray(dtypes).max())
  377. if not hasattr(arr, '_mask'):
  378. result = np.asarray(outarr, dtype=max_dtypes)
  379. else:
  380. result = asarray(outarr, dtype=max_dtypes)
  381. result.fill_value = ma.default_fill_value(result)
  382. return result
  383. apply_along_axis.__doc__ = np.apply_along_axis.__doc__
  384. def apply_over_axes(func, a, axes):
  385. """
  386. (This docstring will be overwritten)
  387. """
  388. val = asarray(a)
  389. N = a.ndim
  390. if array(axes).ndim == 0:
  391. axes = (axes,)
  392. for axis in axes:
  393. if axis < 0:
  394. axis = N + axis
  395. args = (val, axis)
  396. res = func(*args)
  397. if res.ndim == val.ndim:
  398. val = res
  399. else:
  400. res = ma.expand_dims(res, axis)
  401. if res.ndim == val.ndim:
  402. val = res
  403. else:
  404. raise ValueError("function is not returning "
  405. "an array of the correct shape")
  406. return val
  407. if apply_over_axes.__doc__ is not None:
  408. apply_over_axes.__doc__ = np.apply_over_axes.__doc__[
  409. :np.apply_over_axes.__doc__.find('Notes')].rstrip() + \
  410. """
  411. Examples
  412. --------
  413. >>> a = np.ma.arange(24).reshape(2,3,4)
  414. >>> a[:,0,1] = np.ma.masked
  415. >>> a[:,1,:] = np.ma.masked
  416. >>> a
  417. masked_array(
  418. data=[[[0, --, 2, 3],
  419. [--, --, --, --],
  420. [8, 9, 10, 11]],
  421. [[12, --, 14, 15],
  422. [--, --, --, --],
  423. [20, 21, 22, 23]]],
  424. mask=[[[False, True, False, False],
  425. [ True, True, True, True],
  426. [False, False, False, False]],
  427. [[False, True, False, False],
  428. [ True, True, True, True],
  429. [False, False, False, False]]],
  430. fill_value=999999)
  431. >>> np.ma.apply_over_axes(np.ma.sum, a, [0,2])
  432. masked_array(
  433. data=[[[46],
  434. [--],
  435. [124]]],
  436. mask=[[[False],
  437. [ True],
  438. [False]]],
  439. fill_value=999999)
  440. Tuple axis arguments to ufuncs are equivalent:
  441. >>> np.ma.sum(a, axis=(0,2)).reshape((1,-1,1))
  442. masked_array(
  443. data=[[[46],
  444. [--],
  445. [124]]],
  446. mask=[[[False],
  447. [ True],
  448. [False]]],
  449. fill_value=999999)
  450. """
  451. def average(a, axis=None, weights=None, returned=False):
  452. """
  453. Return the weighted average of array over the given axis.
  454. Parameters
  455. ----------
  456. a : array_like
  457. Data to be averaged.
  458. Masked entries are not taken into account in the computation.
  459. axis : int, optional
  460. Axis along which to average `a`. If None, averaging is done over
  461. the flattened array.
  462. weights : array_like, optional
  463. The importance that each element has in the computation of the average.
  464. The weights array can either be 1-D (in which case its length must be
  465. the size of `a` along the given axis) or of the same shape as `a`.
  466. If ``weights=None``, then all data in `a` are assumed to have a
  467. weight equal to one. The 1-D calculation is::
  468. avg = sum(a * weights) / sum(weights)
  469. The only constraint on `weights` is that `sum(weights)` must not be 0.
  470. returned : bool, optional
  471. Flag indicating whether a tuple ``(result, sum of weights)``
  472. should be returned as output (True), or just the result (False).
  473. Default is False.
  474. Returns
  475. -------
  476. average, [sum_of_weights] : (tuple of) scalar or MaskedArray
  477. The average along the specified axis. When returned is `True`,
  478. return a tuple with the average as the first element and the sum
  479. of the weights as the second element. The return type is `np.float64`
  480. if `a` is of integer type and floats smaller than `float64`, or the
  481. input data-type, otherwise. If returned, `sum_of_weights` is always
  482. `float64`.
  483. Examples
  484. --------
  485. >>> a = np.ma.array([1., 2., 3., 4.], mask=[False, False, True, True])
  486. >>> np.ma.average(a, weights=[3, 1, 0, 0])
  487. 1.25
  488. >>> x = np.ma.arange(6.).reshape(3, 2)
  489. >>> x
  490. masked_array(
  491. data=[[0., 1.],
  492. [2., 3.],
  493. [4., 5.]],
  494. mask=False,
  495. fill_value=1e+20)
  496. >>> avg, sumweights = np.ma.average(x, axis=0, weights=[1, 2, 3],
  497. ... returned=True)
  498. >>> avg
  499. masked_array(data=[2.6666666666666665, 3.6666666666666665],
  500. mask=[False, False],
  501. fill_value=1e+20)
  502. """
  503. a = asarray(a)
  504. m = getmask(a)
  505. # inspired by 'average' in numpy/lib/function_base.py
  506. if weights is None:
  507. avg = a.mean(axis)
  508. scl = avg.dtype.type(a.count(axis))
  509. else:
  510. wgt = np.asanyarray(weights)
  511. if issubclass(a.dtype.type, (np.integer, np.bool_)):
  512. result_dtype = np.result_type(a.dtype, wgt.dtype, 'f8')
  513. else:
  514. result_dtype = np.result_type(a.dtype, wgt.dtype)
  515. # Sanity checks
  516. if a.shape != wgt.shape:
  517. if axis is None:
  518. raise TypeError(
  519. "Axis must be specified when shapes of a and weights "
  520. "differ.")
  521. if wgt.ndim != 1:
  522. raise TypeError(
  523. "1D weights expected when shapes of a and weights differ.")
  524. if wgt.shape[0] != a.shape[axis]:
  525. raise ValueError(
  526. "Length of weights not compatible with specified axis.")
  527. # setup wgt to broadcast along axis
  528. wgt = np.broadcast_to(wgt, (a.ndim-1)*(1,) + wgt.shape)
  529. wgt = wgt.swapaxes(-1, axis)
  530. if m is not nomask:
  531. wgt = wgt*(~a.mask)
  532. scl = wgt.sum(axis=axis, dtype=result_dtype)
  533. avg = np.multiply(a, wgt, dtype=result_dtype).sum(axis)/scl
  534. if returned:
  535. if scl.shape != avg.shape:
  536. scl = np.broadcast_to(scl, avg.shape).copy()
  537. return avg, scl
  538. else:
  539. return avg
  540. def median(a, axis=None, out=None, overwrite_input=False, keepdims=False):
  541. """
  542. Compute the median along the specified axis.
  543. Returns the median of the array elements.
  544. Parameters
  545. ----------
  546. a : array_like
  547. Input array or object that can be converted to an array.
  548. axis : int, optional
  549. Axis along which the medians are computed. The default (None) is
  550. to compute the median along a flattened version of the array.
  551. out : ndarray, optional
  552. Alternative output array in which to place the result. It must
  553. have the same shape and buffer length as the expected output
  554. but the type will be cast if necessary.
  555. overwrite_input : bool, optional
  556. If True, then allow use of memory of input array (a) for
  557. calculations. The input array will be modified by the call to
  558. median. This will save memory when you do not need to preserve
  559. the contents of the input array. Treat the input as undefined,
  560. but it will probably be fully or partially sorted. Default is
  561. False. Note that, if `overwrite_input` is True, and the input
  562. is not already an `ndarray`, an error will be raised.
  563. keepdims : bool, optional
  564. If this is set to True, the axes which are reduced are left
  565. in the result as dimensions with size one. With this option,
  566. the result will broadcast correctly against the input array.
  567. .. versionadded:: 1.10.0
  568. Returns
  569. -------
  570. median : ndarray
  571. A new array holding the result is returned unless out is
  572. specified, in which case a reference to out is returned.
  573. Return data-type is `float64` for integers and floats smaller than
  574. `float64`, or the input data-type, otherwise.
  575. See Also
  576. --------
  577. mean
  578. Notes
  579. -----
  580. Given a vector ``V`` with ``N`` non masked values, the median of ``V``
  581. is the middle value of a sorted copy of ``V`` (``Vs``) - i.e.
  582. ``Vs[(N-1)/2]``, when ``N`` is odd, or ``{Vs[N/2 - 1] + Vs[N/2]}/2``
  583. when ``N`` is even.
  584. Examples
  585. --------
  586. >>> x = np.ma.array(np.arange(8), mask=[0]*4 + [1]*4)
  587. >>> np.ma.median(x)
  588. 1.5
  589. >>> x = np.ma.array(np.arange(10).reshape(2, 5), mask=[0]*6 + [1]*4)
  590. >>> np.ma.median(x)
  591. 2.5
  592. >>> np.ma.median(x, axis=-1, overwrite_input=True)
  593. masked_array(data=[2.0, 5.0],
  594. mask=[False, False],
  595. fill_value=1e+20)
  596. """
  597. if not hasattr(a, 'mask'):
  598. m = np.median(getdata(a, subok=True), axis=axis,
  599. out=out, overwrite_input=overwrite_input,
  600. keepdims=keepdims)
  601. if isinstance(m, np.ndarray) and 1 <= m.ndim:
  602. return masked_array(m, copy=False)
  603. else:
  604. return m
  605. r, k = _ureduce(a, func=_median, axis=axis, out=out,
  606. overwrite_input=overwrite_input)
  607. if keepdims:
  608. return r.reshape(k)
  609. else:
  610. return r
  611. def _median(a, axis=None, out=None, overwrite_input=False):
  612. # when an unmasked NaN is present return it, so we need to sort the NaN
  613. # values behind the mask
  614. if np.issubdtype(a.dtype, np.inexact):
  615. fill_value = np.inf
  616. else:
  617. fill_value = None
  618. if overwrite_input:
  619. if axis is None:
  620. asorted = a.ravel()
  621. asorted.sort(fill_value=fill_value)
  622. else:
  623. a.sort(axis=axis, fill_value=fill_value)
  624. asorted = a
  625. else:
  626. asorted = sort(a, axis=axis, fill_value=fill_value)
  627. if axis is None:
  628. axis = 0
  629. else:
  630. axis = normalize_axis_index(axis, asorted.ndim)
  631. if asorted.shape[axis] == 0:
  632. # for empty axis integer indices fail so use slicing to get same result
  633. # as median (which is mean of empty slice = nan)
  634. indexer = [slice(None)] * asorted.ndim
  635. indexer[axis] = slice(0, 0)
  636. indexer = tuple(indexer)
  637. return np.ma.mean(asorted[indexer], axis=axis, out=out)
  638. if asorted.ndim == 1:
  639. counts = count(asorted)
  640. idx, odd = divmod(count(asorted), 2)
  641. mid = asorted[idx + odd - 1:idx + 1]
  642. if np.issubdtype(asorted.dtype, np.inexact) and asorted.size > 0:
  643. # avoid inf / x = masked
  644. s = mid.sum(out=out)
  645. if not odd:
  646. s = np.true_divide(s, 2., casting='safe', out=out)
  647. s = np.lib.utils._median_nancheck(asorted, s, axis, out)
  648. else:
  649. s = mid.mean(out=out)
  650. # if result is masked either the input contained enough
  651. # minimum_fill_value so that it would be the median or all values
  652. # masked
  653. if np.ma.is_masked(s) and not np.all(asorted.mask):
  654. return np.ma.minimum_fill_value(asorted)
  655. return s
  656. counts = count(asorted, axis=axis, keepdims=True)
  657. h = counts // 2
  658. # duplicate high if odd number of elements so mean does nothing
  659. odd = counts % 2 == 1
  660. l = np.where(odd, h, h-1)
  661. lh = np.concatenate([l,h], axis=axis)
  662. # get low and high median
  663. low_high = np.take_along_axis(asorted, lh, axis=axis)
  664. def replace_masked(s):
  665. # Replace masked entries with minimum_full_value unless it all values
  666. # are masked. This is required as the sort order of values equal or
  667. # larger than the fill value is undefined and a valid value placed
  668. # elsewhere, e.g. [4, --, inf].
  669. if np.ma.is_masked(s):
  670. rep = (~np.all(asorted.mask, axis=axis, keepdims=True)) & s.mask
  671. s.data[rep] = np.ma.minimum_fill_value(asorted)
  672. s.mask[rep] = False
  673. replace_masked(low_high)
  674. if np.issubdtype(asorted.dtype, np.inexact):
  675. # avoid inf / x = masked
  676. s = np.ma.sum(low_high, axis=axis, out=out)
  677. np.true_divide(s.data, 2., casting='unsafe', out=s.data)
  678. s = np.lib.utils._median_nancheck(asorted, s, axis, out)
  679. else:
  680. s = np.ma.mean(low_high, axis=axis, out=out)
  681. return s
  682. def compress_nd(x, axis=None):
  683. """Suppress slices from multiple dimensions which contain masked values.
  684. Parameters
  685. ----------
  686. x : array_like, MaskedArray
  687. The array to operate on. If not a MaskedArray instance (or if no array
  688. elements are masked), `x` is interpreted as a MaskedArray with `mask`
  689. set to `nomask`.
  690. axis : tuple of ints or int, optional
  691. Which dimensions to suppress slices from can be configured with this
  692. parameter.
  693. - If axis is a tuple of ints, those are the axes to suppress slices from.
  694. - If axis is an int, then that is the only axis to suppress slices from.
  695. - If axis is None, all axis are selected.
  696. Returns
  697. -------
  698. compress_array : ndarray
  699. The compressed array.
  700. """
  701. x = asarray(x)
  702. m = getmask(x)
  703. # Set axis to tuple of ints
  704. if axis is None:
  705. axis = tuple(range(x.ndim))
  706. else:
  707. axis = normalize_axis_tuple(axis, x.ndim)
  708. # Nothing is masked: return x
  709. if m is nomask or not m.any():
  710. return x._data
  711. # All is masked: return empty
  712. if m.all():
  713. return nxarray([])
  714. # Filter elements through boolean indexing
  715. data = x._data
  716. for ax in axis:
  717. axes = tuple(list(range(ax)) + list(range(ax + 1, x.ndim)))
  718. data = data[(slice(None),)*ax + (~m.any(axis=axes),)]
  719. return data
  720. def compress_rowcols(x, axis=None):
  721. """
  722. Suppress the rows and/or columns of a 2-D array that contain
  723. masked values.
  724. The suppression behavior is selected with the `axis` parameter.
  725. - If axis is None, both rows and columns are suppressed.
  726. - If axis is 0, only rows are suppressed.
  727. - If axis is 1 or -1, only columns are suppressed.
  728. Parameters
  729. ----------
  730. x : array_like, MaskedArray
  731. The array to operate on. If not a MaskedArray instance (or if no array
  732. elements are masked), `x` is interpreted as a MaskedArray with
  733. `mask` set to `nomask`. Must be a 2D array.
  734. axis : int, optional
  735. Axis along which to perform the operation. Default is None.
  736. Returns
  737. -------
  738. compressed_array : ndarray
  739. The compressed array.
  740. Examples
  741. --------
  742. >>> x = np.ma.array(np.arange(9).reshape(3, 3), mask=[[1, 0, 0],
  743. ... [1, 0, 0],
  744. ... [0, 0, 0]])
  745. >>> x
  746. masked_array(
  747. data=[[--, 1, 2],
  748. [--, 4, 5],
  749. [6, 7, 8]],
  750. mask=[[ True, False, False],
  751. [ True, False, False],
  752. [False, False, False]],
  753. fill_value=999999)
  754. >>> np.ma.compress_rowcols(x)
  755. array([[7, 8]])
  756. >>> np.ma.compress_rowcols(x, 0)
  757. array([[6, 7, 8]])
  758. >>> np.ma.compress_rowcols(x, 1)
  759. array([[1, 2],
  760. [4, 5],
  761. [7, 8]])
  762. """
  763. if asarray(x).ndim != 2:
  764. raise NotImplementedError("compress_rowcols works for 2D arrays only.")
  765. return compress_nd(x, axis=axis)
  766. def compress_rows(a):
  767. """
  768. Suppress whole rows of a 2-D array that contain masked values.
  769. This is equivalent to ``np.ma.compress_rowcols(a, 0)``, see
  770. `compress_rowcols` for details.
  771. See Also
  772. --------
  773. compress_rowcols
  774. """
  775. a = asarray(a)
  776. if a.ndim != 2:
  777. raise NotImplementedError("compress_rows works for 2D arrays only.")
  778. return compress_rowcols(a, 0)
  779. def compress_cols(a):
  780. """
  781. Suppress whole columns of a 2-D array that contain masked values.
  782. This is equivalent to ``np.ma.compress_rowcols(a, 1)``, see
  783. `compress_rowcols` for details.
  784. See Also
  785. --------
  786. compress_rowcols
  787. """
  788. a = asarray(a)
  789. if a.ndim != 2:
  790. raise NotImplementedError("compress_cols works for 2D arrays only.")
  791. return compress_rowcols(a, 1)
  792. def mask_rows(a, axis=np._NoValue):
  793. """
  794. Mask rows of a 2D array that contain masked values.
  795. This function is a shortcut to ``mask_rowcols`` with `axis` equal to 0.
  796. See Also
  797. --------
  798. mask_rowcols : Mask rows and/or columns of a 2D array.
  799. masked_where : Mask where a condition is met.
  800. Examples
  801. --------
  802. >>> import numpy.ma as ma
  803. >>> a = np.zeros((3, 3), dtype=int)
  804. >>> a[1, 1] = 1
  805. >>> a
  806. array([[0, 0, 0],
  807. [0, 1, 0],
  808. [0, 0, 0]])
  809. >>> a = ma.masked_equal(a, 1)
  810. >>> a
  811. masked_array(
  812. data=[[0, 0, 0],
  813. [0, --, 0],
  814. [0, 0, 0]],
  815. mask=[[False, False, False],
  816. [False, True, False],
  817. [False, False, False]],
  818. fill_value=1)
  819. >>> ma.mask_rows(a)
  820. masked_array(
  821. data=[[0, 0, 0],
  822. [--, --, --],
  823. [0, 0, 0]],
  824. mask=[[False, False, False],
  825. [ True, True, True],
  826. [False, False, False]],
  827. fill_value=1)
  828. """
  829. if axis is not np._NoValue:
  830. # remove the axis argument when this deprecation expires
  831. # NumPy 1.18.0, 2019-11-28
  832. warnings.warn(
  833. "The axis argument has always been ignored, in future passing it "
  834. "will raise TypeError", DeprecationWarning, stacklevel=2)
  835. return mask_rowcols(a, 0)
  836. def mask_cols(a, axis=np._NoValue):
  837. """
  838. Mask columns of a 2D array that contain masked values.
  839. This function is a shortcut to ``mask_rowcols`` with `axis` equal to 1.
  840. See Also
  841. --------
  842. mask_rowcols : Mask rows and/or columns of a 2D array.
  843. masked_where : Mask where a condition is met.
  844. Examples
  845. --------
  846. >>> import numpy.ma as ma
  847. >>> a = np.zeros((3, 3), dtype=int)
  848. >>> a[1, 1] = 1
  849. >>> a
  850. array([[0, 0, 0],
  851. [0, 1, 0],
  852. [0, 0, 0]])
  853. >>> a = ma.masked_equal(a, 1)
  854. >>> a
  855. masked_array(
  856. data=[[0, 0, 0],
  857. [0, --, 0],
  858. [0, 0, 0]],
  859. mask=[[False, False, False],
  860. [False, True, False],
  861. [False, False, False]],
  862. fill_value=1)
  863. >>> ma.mask_cols(a)
  864. masked_array(
  865. data=[[0, --, 0],
  866. [0, --, 0],
  867. [0, --, 0]],
  868. mask=[[False, True, False],
  869. [False, True, False],
  870. [False, True, False]],
  871. fill_value=1)
  872. """
  873. if axis is not np._NoValue:
  874. # remove the axis argument when this deprecation expires
  875. # NumPy 1.18.0, 2019-11-28
  876. warnings.warn(
  877. "The axis argument has always been ignored, in future passing it "
  878. "will raise TypeError", DeprecationWarning, stacklevel=2)
  879. return mask_rowcols(a, 1)
  880. #####--------------------------------------------------------------------------
  881. #---- --- arraysetops ---
  882. #####--------------------------------------------------------------------------
  883. def ediff1d(arr, to_end=None, to_begin=None):
  884. """
  885. Compute the differences between consecutive elements of an array.
  886. This function is the equivalent of `numpy.ediff1d` that takes masked
  887. values into account, see `numpy.ediff1d` for details.
  888. See Also
  889. --------
  890. numpy.ediff1d : Equivalent function for ndarrays.
  891. """
  892. arr = ma.asanyarray(arr).flat
  893. ed = arr[1:] - arr[:-1]
  894. arrays = [ed]
  895. #
  896. if to_begin is not None:
  897. arrays.insert(0, to_begin)
  898. if to_end is not None:
  899. arrays.append(to_end)
  900. #
  901. if len(arrays) != 1:
  902. # We'll save ourselves a copy of a potentially large array in the common
  903. # case where neither to_begin or to_end was given.
  904. ed = hstack(arrays)
  905. #
  906. return ed
  907. def unique(ar1, return_index=False, return_inverse=False):
  908. """
  909. Finds the unique elements of an array.
  910. Masked values are considered the same element (masked). The output array
  911. is always a masked array. See `numpy.unique` for more details.
  912. See Also
  913. --------
  914. numpy.unique : Equivalent function for ndarrays.
  915. """
  916. output = np.unique(ar1,
  917. return_index=return_index,
  918. return_inverse=return_inverse)
  919. if isinstance(output, tuple):
  920. output = list(output)
  921. output[0] = output[0].view(MaskedArray)
  922. output = tuple(output)
  923. else:
  924. output = output.view(MaskedArray)
  925. return output
  926. def intersect1d(ar1, ar2, assume_unique=False):
  927. """
  928. Returns the unique elements common to both arrays.
  929. Masked values are considered equal one to the other.
  930. The output is always a masked array.
  931. See `numpy.intersect1d` for more details.
  932. See Also
  933. --------
  934. numpy.intersect1d : Equivalent function for ndarrays.
  935. Examples
  936. --------
  937. >>> x = np.ma.array([1, 3, 3, 3], mask=[0, 0, 0, 1])
  938. >>> y = np.ma.array([3, 1, 1, 1], mask=[0, 0, 0, 1])
  939. >>> np.ma.intersect1d(x, y)
  940. masked_array(data=[1, 3, --],
  941. mask=[False, False, True],
  942. fill_value=999999)
  943. """
  944. if assume_unique:
  945. aux = ma.concatenate((ar1, ar2))
  946. else:
  947. # Might be faster than unique( intersect1d( ar1, ar2 ) )?
  948. aux = ma.concatenate((unique(ar1), unique(ar2)))
  949. aux.sort()
  950. return aux[:-1][aux[1:] == aux[:-1]]
  951. def setxor1d(ar1, ar2, assume_unique=False):
  952. """
  953. Set exclusive-or of 1-D arrays with unique elements.
  954. The output is always a masked array. See `numpy.setxor1d` for more details.
  955. See Also
  956. --------
  957. numpy.setxor1d : Equivalent function for ndarrays.
  958. """
  959. if not assume_unique:
  960. ar1 = unique(ar1)
  961. ar2 = unique(ar2)
  962. aux = ma.concatenate((ar1, ar2))
  963. if aux.size == 0:
  964. return aux
  965. aux.sort()
  966. auxf = aux.filled()
  967. # flag = ediff1d( aux, to_end = 1, to_begin = 1 ) == 0
  968. flag = ma.concatenate(([True], (auxf[1:] != auxf[:-1]), [True]))
  969. # flag2 = ediff1d( flag ) == 0
  970. flag2 = (flag[1:] == flag[:-1])
  971. return aux[flag2]
  972. def in1d(ar1, ar2, assume_unique=False, invert=False):
  973. """
  974. Test whether each element of an array is also present in a second
  975. array.
  976. The output is always a masked array. See `numpy.in1d` for more details.
  977. We recommend using :func:`isin` instead of `in1d` for new code.
  978. See Also
  979. --------
  980. isin : Version of this function that preserves the shape of ar1.
  981. numpy.in1d : Equivalent function for ndarrays.
  982. Notes
  983. -----
  984. .. versionadded:: 1.4.0
  985. """
  986. if not assume_unique:
  987. ar1, rev_idx = unique(ar1, return_inverse=True)
  988. ar2 = unique(ar2)
  989. ar = ma.concatenate((ar1, ar2))
  990. # We need this to be a stable sort, so always use 'mergesort'
  991. # here. The values from the first array should always come before
  992. # the values from the second array.
  993. order = ar.argsort(kind='mergesort')
  994. sar = ar[order]
  995. if invert:
  996. bool_ar = (sar[1:] != sar[:-1])
  997. else:
  998. bool_ar = (sar[1:] == sar[:-1])
  999. flag = ma.concatenate((bool_ar, [invert]))
  1000. indx = order.argsort(kind='mergesort')[:len(ar1)]
  1001. if assume_unique:
  1002. return flag[indx]
  1003. else:
  1004. return flag[indx][rev_idx]
  1005. def isin(element, test_elements, assume_unique=False, invert=False):
  1006. """
  1007. Calculates `element in test_elements`, broadcasting over
  1008. `element` only.
  1009. The output is always a masked array of the same shape as `element`.
  1010. See `numpy.isin` for more details.
  1011. See Also
  1012. --------
  1013. in1d : Flattened version of this function.
  1014. numpy.isin : Equivalent function for ndarrays.
  1015. Notes
  1016. -----
  1017. .. versionadded:: 1.13.0
  1018. """
  1019. element = ma.asarray(element)
  1020. return in1d(element, test_elements, assume_unique=assume_unique,
  1021. invert=invert).reshape(element.shape)
  1022. def union1d(ar1, ar2):
  1023. """
  1024. Union of two arrays.
  1025. The output is always a masked array. See `numpy.union1d` for more details.
  1026. See also
  1027. --------
  1028. numpy.union1d : Equivalent function for ndarrays.
  1029. """
  1030. return unique(ma.concatenate((ar1, ar2), axis=None))
  1031. def setdiff1d(ar1, ar2, assume_unique=False):
  1032. """
  1033. Set difference of 1D arrays with unique elements.
  1034. The output is always a masked array. See `numpy.setdiff1d` for more
  1035. details.
  1036. See Also
  1037. --------
  1038. numpy.setdiff1d : Equivalent function for ndarrays.
  1039. Examples
  1040. --------
  1041. >>> x = np.ma.array([1, 2, 3, 4], mask=[0, 1, 0, 1])
  1042. >>> np.ma.setdiff1d(x, [1, 2])
  1043. masked_array(data=[3, --],
  1044. mask=[False, True],
  1045. fill_value=999999)
  1046. """
  1047. if assume_unique:
  1048. ar1 = ma.asarray(ar1).ravel()
  1049. else:
  1050. ar1 = unique(ar1)
  1051. ar2 = unique(ar2)
  1052. return ar1[in1d(ar1, ar2, assume_unique=True, invert=True)]
  1053. ###############################################################################
  1054. # Covariance #
  1055. ###############################################################################
  1056. def _covhelper(x, y=None, rowvar=True, allow_masked=True):
  1057. """
  1058. Private function for the computation of covariance and correlation
  1059. coefficients.
  1060. """
  1061. x = ma.array(x, ndmin=2, copy=True, dtype=float)
  1062. xmask = ma.getmaskarray(x)
  1063. # Quick exit if we can't process masked data
  1064. if not allow_masked and xmask.any():
  1065. raise ValueError("Cannot process masked data.")
  1066. #
  1067. if x.shape[0] == 1:
  1068. rowvar = True
  1069. # Make sure that rowvar is either 0 or 1
  1070. rowvar = int(bool(rowvar))
  1071. axis = 1 - rowvar
  1072. if rowvar:
  1073. tup = (slice(None), None)
  1074. else:
  1075. tup = (None, slice(None))
  1076. #
  1077. if y is None:
  1078. xnotmask = np.logical_not(xmask).astype(int)
  1079. else:
  1080. y = array(y, copy=False, ndmin=2, dtype=float)
  1081. ymask = ma.getmaskarray(y)
  1082. if not allow_masked and ymask.any():
  1083. raise ValueError("Cannot process masked data.")
  1084. if xmask.any() or ymask.any():
  1085. if y.shape == x.shape:
  1086. # Define some common mask
  1087. common_mask = np.logical_or(xmask, ymask)
  1088. if common_mask is not nomask:
  1089. xmask = x._mask = y._mask = ymask = common_mask
  1090. x._sharedmask = False
  1091. y._sharedmask = False
  1092. x = ma.concatenate((x, y), axis)
  1093. xnotmask = np.logical_not(np.concatenate((xmask, ymask), axis)).astype(int)
  1094. x -= x.mean(axis=rowvar)[tup]
  1095. return (x, xnotmask, rowvar)
  1096. def cov(x, y=None, rowvar=True, bias=False, allow_masked=True, ddof=None):
  1097. """
  1098. Estimate the covariance matrix.
  1099. Except for the handling of missing data this function does the same as
  1100. `numpy.cov`. For more details and examples, see `numpy.cov`.
  1101. By default, masked values are recognized as such. If `x` and `y` have the
  1102. same shape, a common mask is allocated: if ``x[i,j]`` is masked, then
  1103. ``y[i,j]`` will also be masked.
  1104. Setting `allow_masked` to False will raise an exception if values are
  1105. missing in either of the input arrays.
  1106. Parameters
  1107. ----------
  1108. x : array_like
  1109. A 1-D or 2-D array containing multiple variables and observations.
  1110. Each row of `x` represents a variable, and each column a single
  1111. observation of all those variables. Also see `rowvar` below.
  1112. y : array_like, optional
  1113. An additional set of variables and observations. `y` has the same
  1114. form as `x`.
  1115. rowvar : bool, optional
  1116. If `rowvar` is True (default), then each row represents a
  1117. variable, with observations in the columns. Otherwise, the relationship
  1118. is transposed: each column represents a variable, while the rows
  1119. contain observations.
  1120. bias : bool, optional
  1121. Default normalization (False) is by ``(N-1)``, where ``N`` is the
  1122. number of observations given (unbiased estimate). If `bias` is True,
  1123. then normalization is by ``N``. This keyword can be overridden by
  1124. the keyword ``ddof`` in numpy versions >= 1.5.
  1125. allow_masked : bool, optional
  1126. If True, masked values are propagated pair-wise: if a value is masked
  1127. in `x`, the corresponding value is masked in `y`.
  1128. If False, raises a `ValueError` exception when some values are missing.
  1129. ddof : {None, int}, optional
  1130. If not ``None`` normalization is by ``(N - ddof)``, where ``N`` is
  1131. the number of observations; this overrides the value implied by
  1132. ``bias``. The default value is ``None``.
  1133. .. versionadded:: 1.5
  1134. Raises
  1135. ------
  1136. ValueError
  1137. Raised if some values are missing and `allow_masked` is False.
  1138. See Also
  1139. --------
  1140. numpy.cov
  1141. """
  1142. # Check inputs
  1143. if ddof is not None and ddof != int(ddof):
  1144. raise ValueError("ddof must be an integer")
  1145. # Set up ddof
  1146. if ddof is None:
  1147. if bias:
  1148. ddof = 0
  1149. else:
  1150. ddof = 1
  1151. (x, xnotmask, rowvar) = _covhelper(x, y, rowvar, allow_masked)
  1152. if not rowvar:
  1153. fact = np.dot(xnotmask.T, xnotmask) * 1. - ddof
  1154. result = (dot(x.T, x.conj(), strict=False) / fact).squeeze()
  1155. else:
  1156. fact = np.dot(xnotmask, xnotmask.T) * 1. - ddof
  1157. result = (dot(x, x.T.conj(), strict=False) / fact).squeeze()
  1158. return result
  1159. def corrcoef(x, y=None, rowvar=True, bias=np._NoValue, allow_masked=True,
  1160. ddof=np._NoValue):
  1161. """
  1162. Return Pearson product-moment correlation coefficients.
  1163. Except for the handling of missing data this function does the same as
  1164. `numpy.corrcoef`. For more details and examples, see `numpy.corrcoef`.
  1165. Parameters
  1166. ----------
  1167. x : array_like
  1168. A 1-D or 2-D array containing multiple variables and observations.
  1169. Each row of `x` represents a variable, and each column a single
  1170. observation of all those variables. Also see `rowvar` below.
  1171. y : array_like, optional
  1172. An additional set of variables and observations. `y` has the same
  1173. shape as `x`.
  1174. rowvar : bool, optional
  1175. If `rowvar` is True (default), then each row represents a
  1176. variable, with observations in the columns. Otherwise, the relationship
  1177. is transposed: each column represents a variable, while the rows
  1178. contain observations.
  1179. bias : _NoValue, optional
  1180. Has no effect, do not use.
  1181. .. deprecated:: 1.10.0
  1182. allow_masked : bool, optional
  1183. If True, masked values are propagated pair-wise: if a value is masked
  1184. in `x`, the corresponding value is masked in `y`.
  1185. If False, raises an exception. Because `bias` is deprecated, this
  1186. argument needs to be treated as keyword only to avoid a warning.
  1187. ddof : _NoValue, optional
  1188. Has no effect, do not use.
  1189. .. deprecated:: 1.10.0
  1190. See Also
  1191. --------
  1192. numpy.corrcoef : Equivalent function in top-level NumPy module.
  1193. cov : Estimate the covariance matrix.
  1194. Notes
  1195. -----
  1196. This function accepts but discards arguments `bias` and `ddof`. This is
  1197. for backwards compatibility with previous versions of this function. These
  1198. arguments had no effect on the return values of the function and can be
  1199. safely ignored in this and previous versions of numpy.
  1200. """
  1201. msg = 'bias and ddof have no effect and are deprecated'
  1202. if bias is not np._NoValue or ddof is not np._NoValue:
  1203. # 2015-03-15, 1.10
  1204. warnings.warn(msg, DeprecationWarning, stacklevel=2)
  1205. # Get the data
  1206. (x, xnotmask, rowvar) = _covhelper(x, y, rowvar, allow_masked)
  1207. # Compute the covariance matrix
  1208. if not rowvar:
  1209. fact = np.dot(xnotmask.T, xnotmask) * 1.
  1210. c = (dot(x.T, x.conj(), strict=False) / fact).squeeze()
  1211. else:
  1212. fact = np.dot(xnotmask, xnotmask.T) * 1.
  1213. c = (dot(x, x.T.conj(), strict=False) / fact).squeeze()
  1214. # Check whether we have a scalar
  1215. try:
  1216. diag = ma.diagonal(c)
  1217. except ValueError:
  1218. return 1
  1219. #
  1220. if xnotmask.all():
  1221. _denom = ma.sqrt(ma.multiply.outer(diag, diag))
  1222. else:
  1223. _denom = diagflat(diag)
  1224. _denom._sharedmask = False # We know return is always a copy
  1225. n = x.shape[1 - rowvar]
  1226. if rowvar:
  1227. for i in range(n - 1):
  1228. for j in range(i + 1, n):
  1229. _x = mask_cols(vstack((x[i], x[j]))).var(axis=1)
  1230. _denom[i, j] = _denom[j, i] = ma.sqrt(ma.multiply.reduce(_x))
  1231. else:
  1232. for i in range(n - 1):
  1233. for j in range(i + 1, n):
  1234. _x = mask_cols(
  1235. vstack((x[:, i], x[:, j]))).var(axis=1)
  1236. _denom[i, j] = _denom[j, i] = ma.sqrt(ma.multiply.reduce(_x))
  1237. return c / _denom
  1238. #####--------------------------------------------------------------------------
  1239. #---- --- Concatenation helpers ---
  1240. #####--------------------------------------------------------------------------
  1241. class MAxisConcatenator(AxisConcatenator):
  1242. """
  1243. Translate slice objects to concatenation along an axis.
  1244. For documentation on usage, see `mr_class`.
  1245. See Also
  1246. --------
  1247. mr_class
  1248. """
  1249. concatenate = staticmethod(concatenate)
  1250. @classmethod
  1251. def makemat(cls, arr):
  1252. # There used to be a view as np.matrix here, but we may eventually
  1253. # deprecate that class. In preparation, we use the unmasked version
  1254. # to construct the matrix (with copy=False for backwards compatibility
  1255. # with the .view)
  1256. data = super(MAxisConcatenator, cls).makemat(arr.data, copy=False)
  1257. return array(data, mask=arr.mask)
  1258. def __getitem__(self, key):
  1259. # matrix builder syntax, like 'a, b; c, d'
  1260. if isinstance(key, str):
  1261. raise MAError("Unavailable for masked array.")
  1262. return super(MAxisConcatenator, self).__getitem__(key)
  1263. class mr_class(MAxisConcatenator):
  1264. """
  1265. Translate slice objects to concatenation along the first axis.
  1266. This is the masked array version of `lib.index_tricks.RClass`.
  1267. See Also
  1268. --------
  1269. lib.index_tricks.RClass
  1270. Examples
  1271. --------
  1272. >>> np.ma.mr_[np.ma.array([1,2,3]), 0, 0, np.ma.array([4,5,6])]
  1273. masked_array(data=[1, 2, 3, ..., 4, 5, 6],
  1274. mask=False,
  1275. fill_value=999999)
  1276. """
  1277. def __init__(self):
  1278. MAxisConcatenator.__init__(self, 0)
  1279. mr_ = mr_class()
  1280. #####--------------------------------------------------------------------------
  1281. #---- Find unmasked data ---
  1282. #####--------------------------------------------------------------------------
  1283. def flatnotmasked_edges(a):
  1284. """
  1285. Find the indices of the first and last unmasked values.
  1286. Expects a 1-D `MaskedArray`, returns None if all values are masked.
  1287. Parameters
  1288. ----------
  1289. a : array_like
  1290. Input 1-D `MaskedArray`
  1291. Returns
  1292. -------
  1293. edges : ndarray or None
  1294. The indices of first and last non-masked value in the array.
  1295. Returns None if all values are masked.
  1296. See Also
  1297. --------
  1298. flatnotmasked_contiguous, notmasked_contiguous, notmasked_edges
  1299. clump_masked, clump_unmasked
  1300. Notes
  1301. -----
  1302. Only accepts 1-D arrays.
  1303. Examples
  1304. --------
  1305. >>> a = np.ma.arange(10)
  1306. >>> np.ma.flatnotmasked_edges(a)
  1307. array([0, 9])
  1308. >>> mask = (a < 3) | (a > 8) | (a == 5)
  1309. >>> a[mask] = np.ma.masked
  1310. >>> np.array(a[~a.mask])
  1311. array([3, 4, 6, 7, 8])
  1312. >>> np.ma.flatnotmasked_edges(a)
  1313. array([3, 8])
  1314. >>> a[:] = np.ma.masked
  1315. >>> print(np.ma.flatnotmasked_edges(a))
  1316. None
  1317. """
  1318. m = getmask(a)
  1319. if m is nomask or not np.any(m):
  1320. return np.array([0, a.size - 1])
  1321. unmasked = np.flatnonzero(~m)
  1322. if len(unmasked) > 0:
  1323. return unmasked[[0, -1]]
  1324. else:
  1325. return None
  1326. def notmasked_edges(a, axis=None):
  1327. """
  1328. Find the indices of the first and last unmasked values along an axis.
  1329. If all values are masked, return None. Otherwise, return a list
  1330. of two tuples, corresponding to the indices of the first and last
  1331. unmasked values respectively.
  1332. Parameters
  1333. ----------
  1334. a : array_like
  1335. The input array.
  1336. axis : int, optional
  1337. Axis along which to perform the operation.
  1338. If None (default), applies to a flattened version of the array.
  1339. Returns
  1340. -------
  1341. edges : ndarray or list
  1342. An array of start and end indexes if there are any masked data in
  1343. the array. If there are no masked data in the array, `edges` is a
  1344. list of the first and last index.
  1345. See Also
  1346. --------
  1347. flatnotmasked_contiguous, flatnotmasked_edges, notmasked_contiguous
  1348. clump_masked, clump_unmasked
  1349. Examples
  1350. --------
  1351. >>> a = np.arange(9).reshape((3, 3))
  1352. >>> m = np.zeros_like(a)
  1353. >>> m[1:, 1:] = 1
  1354. >>> am = np.ma.array(a, mask=m)
  1355. >>> np.array(am[~am.mask])
  1356. array([0, 1, 2, 3, 6])
  1357. >>> np.ma.notmasked_edges(am)
  1358. array([0, 6])
  1359. """
  1360. a = asarray(a)
  1361. if axis is None or a.ndim == 1:
  1362. return flatnotmasked_edges(a)
  1363. m = getmaskarray(a)
  1364. idx = array(np.indices(a.shape), mask=np.asarray([m] * a.ndim))
  1365. return [tuple([idx[i].min(axis).compressed() for i in range(a.ndim)]),
  1366. tuple([idx[i].max(axis).compressed() for i in range(a.ndim)]), ]
  1367. def flatnotmasked_contiguous(a):
  1368. """
  1369. Find contiguous unmasked data in a masked array along the given axis.
  1370. Parameters
  1371. ----------
  1372. a : narray
  1373. The input array.
  1374. Returns
  1375. -------
  1376. slice_list : list
  1377. A sorted sequence of `slice` objects (start index, end index).
  1378. .. versionchanged:: 1.15.0
  1379. Now returns an empty list instead of None for a fully masked array
  1380. See Also
  1381. --------
  1382. flatnotmasked_edges, notmasked_contiguous, notmasked_edges
  1383. clump_masked, clump_unmasked
  1384. Notes
  1385. -----
  1386. Only accepts 2-D arrays at most.
  1387. Examples
  1388. --------
  1389. >>> a = np.ma.arange(10)
  1390. >>> np.ma.flatnotmasked_contiguous(a)
  1391. [slice(0, 10, None)]
  1392. >>> mask = (a < 3) | (a > 8) | (a == 5)
  1393. >>> a[mask] = np.ma.masked
  1394. >>> np.array(a[~a.mask])
  1395. array([3, 4, 6, 7, 8])
  1396. >>> np.ma.flatnotmasked_contiguous(a)
  1397. [slice(3, 5, None), slice(6, 9, None)]
  1398. >>> a[:] = np.ma.masked
  1399. >>> np.ma.flatnotmasked_contiguous(a)
  1400. []
  1401. """
  1402. m = getmask(a)
  1403. if m is nomask:
  1404. return [slice(0, a.size)]
  1405. i = 0
  1406. result = []
  1407. for (k, g) in itertools.groupby(m.ravel()):
  1408. n = len(list(g))
  1409. if not k:
  1410. result.append(slice(i, i + n))
  1411. i += n
  1412. return result
  1413. def notmasked_contiguous(a, axis=None):
  1414. """
  1415. Find contiguous unmasked data in a masked array along the given axis.
  1416. Parameters
  1417. ----------
  1418. a : array_like
  1419. The input array.
  1420. axis : int, optional
  1421. Axis along which to perform the operation.
  1422. If None (default), applies to a flattened version of the array, and this
  1423. is the same as `flatnotmasked_contiguous`.
  1424. Returns
  1425. -------
  1426. endpoints : list
  1427. A list of slices (start and end indexes) of unmasked indexes
  1428. in the array.
  1429. If the input is 2d and axis is specified, the result is a list of lists.
  1430. See Also
  1431. --------
  1432. flatnotmasked_edges, flatnotmasked_contiguous, notmasked_edges
  1433. clump_masked, clump_unmasked
  1434. Notes
  1435. -----
  1436. Only accepts 2-D arrays at most.
  1437. Examples
  1438. --------
  1439. >>> a = np.arange(12).reshape((3, 4))
  1440. >>> mask = np.zeros_like(a)
  1441. >>> mask[1:, :-1] = 1; mask[0, 1] = 1; mask[-1, 0] = 0
  1442. >>> ma = np.ma.array(a, mask=mask)
  1443. >>> ma
  1444. masked_array(
  1445. data=[[0, --, 2, 3],
  1446. [--, --, --, 7],
  1447. [8, --, --, 11]],
  1448. mask=[[False, True, False, False],
  1449. [ True, True, True, False],
  1450. [False, True, True, False]],
  1451. fill_value=999999)
  1452. >>> np.array(ma[~ma.mask])
  1453. array([ 0, 2, 3, 7, 8, 11])
  1454. >>> np.ma.notmasked_contiguous(ma)
  1455. [slice(0, 1, None), slice(2, 4, None), slice(7, 9, None), slice(11, 12, None)]
  1456. >>> np.ma.notmasked_contiguous(ma, axis=0)
  1457. [[slice(0, 1, None), slice(2, 3, None)], [], [slice(0, 1, None)], [slice(0, 3, None)]]
  1458. >>> np.ma.notmasked_contiguous(ma, axis=1)
  1459. [[slice(0, 1, None), slice(2, 4, None)], [slice(3, 4, None)], [slice(0, 1, None), slice(3, 4, None)]]
  1460. """
  1461. a = asarray(a)
  1462. nd = a.ndim
  1463. if nd > 2:
  1464. raise NotImplementedError("Currently limited to atmost 2D array.")
  1465. if axis is None or nd == 1:
  1466. return flatnotmasked_contiguous(a)
  1467. #
  1468. result = []
  1469. #
  1470. other = (axis + 1) % 2
  1471. idx = [0, 0]
  1472. idx[axis] = slice(None, None)
  1473. #
  1474. for i in range(a.shape[other]):
  1475. idx[other] = i
  1476. result.append(flatnotmasked_contiguous(a[tuple(idx)]))
  1477. return result
  1478. def _ezclump(mask):
  1479. """
  1480. Finds the clumps (groups of data with the same values) for a 1D bool array.
  1481. Returns a series of slices.
  1482. """
  1483. if mask.ndim > 1:
  1484. mask = mask.ravel()
  1485. idx = (mask[1:] ^ mask[:-1]).nonzero()
  1486. idx = idx[0] + 1
  1487. if mask[0]:
  1488. if len(idx) == 0:
  1489. return [slice(0, mask.size)]
  1490. r = [slice(0, idx[0])]
  1491. r.extend((slice(left, right)
  1492. for left, right in zip(idx[1:-1:2], idx[2::2])))
  1493. else:
  1494. if len(idx) == 0:
  1495. return []
  1496. r = [slice(left, right) for left, right in zip(idx[:-1:2], idx[1::2])]
  1497. if mask[-1]:
  1498. r.append(slice(idx[-1], mask.size))
  1499. return r
  1500. def clump_unmasked(a):
  1501. """
  1502. Return list of slices corresponding to the unmasked clumps of a 1-D array.
  1503. (A "clump" is defined as a contiguous region of the array).
  1504. Parameters
  1505. ----------
  1506. a : ndarray
  1507. A one-dimensional masked array.
  1508. Returns
  1509. -------
  1510. slices : list of slice
  1511. The list of slices, one for each continuous region of unmasked
  1512. elements in `a`.
  1513. Notes
  1514. -----
  1515. .. versionadded:: 1.4.0
  1516. See Also
  1517. --------
  1518. flatnotmasked_edges, flatnotmasked_contiguous, notmasked_edges
  1519. notmasked_contiguous, clump_masked
  1520. Examples
  1521. --------
  1522. >>> a = np.ma.masked_array(np.arange(10))
  1523. >>> a[[0, 1, 2, 6, 8, 9]] = np.ma.masked
  1524. >>> np.ma.clump_unmasked(a)
  1525. [slice(3, 6, None), slice(7, 8, None)]
  1526. """
  1527. mask = getattr(a, '_mask', nomask)
  1528. if mask is nomask:
  1529. return [slice(0, a.size)]
  1530. return _ezclump(~mask)
  1531. def clump_masked(a):
  1532. """
  1533. Returns a list of slices corresponding to the masked clumps of a 1-D array.
  1534. (A "clump" is defined as a contiguous region of the array).
  1535. Parameters
  1536. ----------
  1537. a : ndarray
  1538. A one-dimensional masked array.
  1539. Returns
  1540. -------
  1541. slices : list of slice
  1542. The list of slices, one for each continuous region of masked elements
  1543. in `a`.
  1544. Notes
  1545. -----
  1546. .. versionadded:: 1.4.0
  1547. See Also
  1548. --------
  1549. flatnotmasked_edges, flatnotmasked_contiguous, notmasked_edges
  1550. notmasked_contiguous, clump_unmasked
  1551. Examples
  1552. --------
  1553. >>> a = np.ma.masked_array(np.arange(10))
  1554. >>> a[[0, 1, 2, 6, 8, 9]] = np.ma.masked
  1555. >>> np.ma.clump_masked(a)
  1556. [slice(0, 3, None), slice(6, 7, None), slice(8, 10, None)]
  1557. """
  1558. mask = ma.getmask(a)
  1559. if mask is nomask:
  1560. return []
  1561. return _ezclump(mask)
  1562. ###############################################################################
  1563. # Polynomial fit #
  1564. ###############################################################################
  1565. def vander(x, n=None):
  1566. """
  1567. Masked values in the input array result in rows of zeros.
  1568. """
  1569. _vander = np.vander(x, n)
  1570. m = getmask(x)
  1571. if m is not nomask:
  1572. _vander[m] = 0
  1573. return _vander
  1574. vander.__doc__ = ma.doc_note(np.vander.__doc__, vander.__doc__)
  1575. def polyfit(x, y, deg, rcond=None, full=False, w=None, cov=False):
  1576. """
  1577. Any masked values in x is propagated in y, and vice-versa.
  1578. """
  1579. x = asarray(x)
  1580. y = asarray(y)
  1581. m = getmask(x)
  1582. if y.ndim == 1:
  1583. m = mask_or(m, getmask(y))
  1584. elif y.ndim == 2:
  1585. my = getmask(mask_rows(y))
  1586. if my is not nomask:
  1587. m = mask_or(m, my[:, 0])
  1588. else:
  1589. raise TypeError("Expected a 1D or 2D array for y!")
  1590. if w is not None:
  1591. w = asarray(w)
  1592. if w.ndim != 1:
  1593. raise TypeError("expected a 1-d array for weights")
  1594. if w.shape[0] != y.shape[0]:
  1595. raise TypeError("expected w and y to have the same length")
  1596. m = mask_or(m, getmask(w))
  1597. if m is not nomask:
  1598. not_m = ~m
  1599. if w is not None:
  1600. w = w[not_m]
  1601. return np.polyfit(x[not_m], y[not_m], deg, rcond, full, w, cov)
  1602. else:
  1603. return np.polyfit(x, y, deg, rcond, full, w, cov)
  1604. polyfit.__doc__ = ma.doc_note(np.polyfit.__doc__, polyfit.__doc__)