extras.py 53 KB

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