Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1""" 

2Masked arrays add-ons. 

3 

4A collection of utilities for `numpy.ma`. 

5 

6:author: Pierre Gerard-Marchant 

7:contact: pierregm_at_uga_dot_edu 

8:version: $Id: extras.py 3473 2007-10-29 15:18:13Z jarrod.millman $ 

9 

10""" 

11__all__ = [ 

12 'apply_along_axis', 'apply_over_axes', 'atleast_1d', 'atleast_2d', 

13 'atleast_3d', 'average', 'clump_masked', 'clump_unmasked', 

14 'column_stack', 'compress_cols', 'compress_nd', 'compress_rowcols', 

15 'compress_rows', 'count_masked', 'corrcoef', 'cov', 'diagflat', 'dot', 

16 'dstack', 'ediff1d', 'flatnotmasked_contiguous', 'flatnotmasked_edges', 

17 'hsplit', 'hstack', 'isin', 'in1d', 'intersect1d', 'mask_cols', 'mask_rowcols', 

18 'mask_rows', 'masked_all', 'masked_all_like', 'median', 'mr_', 

19 'notmasked_contiguous', 'notmasked_edges', 'polyfit', 'row_stack', 

20 'setdiff1d', 'setxor1d', 'stack', 'unique', 'union1d', 'vander', 'vstack', 

21 ] 

22 

23import itertools 

24import warnings 

25 

26from . import core as ma 

27from .core import ( 

28 MaskedArray, MAError, add, array, asarray, concatenate, filled, count, 

29 getmask, getmaskarray, make_mask_descr, masked, masked_array, mask_or, 

30 nomask, ones, sort, zeros, getdata, get_masked_subclass, dot, 

31 mask_rowcols 

32 ) 

33 

34import numpy as np 

35from numpy import ndarray, array as nxarray 

36import numpy.core.umath as umath 

37from numpy.core.multiarray import normalize_axis_index 

38from numpy.core.numeric import normalize_axis_tuple 

39from numpy.lib.function_base import _ureduce 

40from numpy.lib.index_tricks import AxisConcatenator 

41 

42 

43def issequence(seq): 

44 """ 

45 Is seq a sequence (ndarray, list or tuple)? 

46 

47 """ 

48 return isinstance(seq, (ndarray, tuple, list)) 

49 

50 

51def count_masked(arr, axis=None): 

52 """ 

53 Count the number of masked elements along the given axis. 

54 

55 Parameters 

56 ---------- 

57 arr : array_like 

58 An array with (possibly) masked elements. 

59 axis : int, optional 

60 Axis along which to count. If None (default), a flattened 

61 version of the array is used. 

62 

63 Returns 

64 ------- 

65 count : int, ndarray 

66 The total number of masked elements (axis=None) or the number 

67 of masked elements along each slice of the given axis. 

68 

69 See Also 

70 -------- 

71 MaskedArray.count : Count non-masked elements. 

72 

73 Examples 

74 -------- 

75 >>> import numpy.ma as ma 

76 >>> a = np.arange(9).reshape((3,3)) 

77 >>> a = ma.array(a) 

78 >>> a[1, 0] = ma.masked 

79 >>> a[1, 2] = ma.masked 

80 >>> a[2, 1] = ma.masked 

81 >>> a 

82 masked_array( 

83 data=[[0, 1, 2], 

84 [--, 4, --], 

85 [6, --, 8]], 

86 mask=[[False, False, False], 

87 [ True, False, True], 

88 [False, True, False]], 

89 fill_value=999999) 

90 >>> ma.count_masked(a) 

91 3 

92 

93 When the `axis` keyword is used an array is returned. 

94 

95 >>> ma.count_masked(a, axis=0) 

96 array([1, 1, 1]) 

97 >>> ma.count_masked(a, axis=1) 

98 array([0, 2, 1]) 

99 

100 """ 

101 m = getmaskarray(arr) 

102 return m.sum(axis) 

103 

104 

105def masked_all(shape, dtype=float): 

106 """ 

107 Empty masked array with all elements masked. 

108 

109 Return an empty masked array of the given shape and dtype, where all the 

110 data are masked. 

111 

112 Parameters 

113 ---------- 

114 shape : tuple 

115 Shape of the required MaskedArray. 

116 dtype : dtype, optional 

117 Data type of the output. 

118 

119 Returns 

120 ------- 

121 a : MaskedArray 

122 A masked array with all data masked. 

123 

124 See Also 

125 -------- 

126 masked_all_like : Empty masked array modelled on an existing array. 

127 

128 Examples 

129 -------- 

130 >>> import numpy.ma as ma 

131 >>> ma.masked_all((3, 3)) 

132 masked_array( 

133 data=[[--, --, --], 

134 [--, --, --], 

135 [--, --, --]], 

136 mask=[[ True, True, True], 

137 [ True, True, True], 

138 [ True, True, True]], 

139 fill_value=1e+20, 

140 dtype=float64) 

141 

142 The `dtype` parameter defines the underlying data type. 

143 

144 >>> a = ma.masked_all((3, 3)) 

145 >>> a.dtype 

146 dtype('float64') 

147 >>> a = ma.masked_all((3, 3), dtype=np.int32) 

148 >>> a.dtype 

149 dtype('int32') 

150 

151 """ 

152 a = masked_array(np.empty(shape, dtype), 

153 mask=np.ones(shape, make_mask_descr(dtype))) 

154 return a 

155 

156 

157def masked_all_like(arr): 

158 """ 

159 Empty masked array with the properties of an existing array. 

160 

161 Return an empty masked array of the same shape and dtype as 

162 the array `arr`, where all the data are masked. 

163 

164 Parameters 

165 ---------- 

166 arr : ndarray 

167 An array describing the shape and dtype of the required MaskedArray. 

168 

169 Returns 

170 ------- 

171 a : MaskedArray 

172 A masked array with all data masked. 

173 

174 Raises 

175 ------ 

176 AttributeError 

177 If `arr` doesn't have a shape attribute (i.e. not an ndarray) 

178 

179 See Also 

180 -------- 

181 masked_all : Empty masked array with all elements masked. 

182 

183 Examples 

184 -------- 

185 >>> import numpy.ma as ma 

186 >>> arr = np.zeros((2, 3), dtype=np.float32) 

187 >>> arr 

188 array([[0., 0., 0.], 

189 [0., 0., 0.]], dtype=float32) 

190 >>> ma.masked_all_like(arr) 

191 masked_array( 

192 data=[[--, --, --], 

193 [--, --, --]], 

194 mask=[[ True, True, True], 

195 [ True, True, True]], 

196 fill_value=1e+20, 

197 dtype=float32) 

198 

199 The dtype of the masked array matches the dtype of `arr`. 

200 

201 >>> arr.dtype 

202 dtype('float32') 

203 >>> ma.masked_all_like(arr).dtype 

204 dtype('float32') 

205 

206 """ 

207 a = np.empty_like(arr).view(MaskedArray) 

208 a._mask = np.ones(a.shape, dtype=make_mask_descr(a.dtype)) 

209 return a 

210 

211 

212#####-------------------------------------------------------------------------- 

213#---- --- Standard functions --- 

214#####-------------------------------------------------------------------------- 

215class _fromnxfunction: 

216 """ 

217 Defines a wrapper to adapt NumPy functions to masked arrays. 

218 

219 

220 An instance of `_fromnxfunction` can be called with the same parameters 

221 as the wrapped NumPy function. The docstring of `newfunc` is adapted from 

222 the wrapped function as well, see `getdoc`. 

223 

224 This class should not be used directly. Instead, one of its extensions that 

225 provides support for a specific type of input should be used. 

226 

227 Parameters 

228 ---------- 

229 funcname : str 

230 The name of the function to be adapted. The function should be 

231 in the NumPy namespace (i.e. ``np.funcname``). 

232 

233 """ 

234 

235 def __init__(self, funcname): 

236 self.__name__ = funcname 

237 self.__doc__ = self.getdoc() 

238 

239 def getdoc(self): 

240 """ 

241 Retrieve the docstring and signature from the function. 

242 

243 The ``__doc__`` attribute of the function is used as the docstring for 

244 the new masked array version of the function. A note on application 

245 of the function to the mask is appended. 

246 

247 .. warning:: 

248 If the function docstring already contained a Notes section, the 

249 new docstring will have two Notes sections instead of appending a note 

250 to the existing section. 

251 

252 Parameters 

253 ---------- 

254 None 

255 

256 """ 

257 npfunc = getattr(np, self.__name__, None) 

258 doc = getattr(npfunc, '__doc__', None) 

259 if doc: 

260 sig = self.__name__ + ma.get_object_signature(npfunc) 

261 locdoc = "Notes\n-----\nThe function is applied to both the _data"\ 

262 " and the _mask, if any." 

263 return '\n'.join((sig, doc, locdoc)) 

264 return 

265 

266 def __call__(self, *args, **params): 

267 pass 

268 

269 

270class _fromnxfunction_single(_fromnxfunction): 

271 """ 

272 A version of `_fromnxfunction` that is called with a single array 

273 argument followed by auxiliary args that are passed verbatim for 

274 both the data and mask calls. 

275 """ 

276 def __call__(self, x, *args, **params): 

277 func = getattr(np, self.__name__) 

278 if isinstance(x, ndarray): 

279 _d = func(x.__array__(), *args, **params) 

280 _m = func(getmaskarray(x), *args, **params) 

281 return masked_array(_d, mask=_m) 

282 else: 

283 _d = func(np.asarray(x), *args, **params) 

284 _m = func(getmaskarray(x), *args, **params) 

285 return masked_array(_d, mask=_m) 

286 

287 

288class _fromnxfunction_seq(_fromnxfunction): 

289 """ 

290 A version of `_fromnxfunction` that is called with a single sequence 

291 of arrays followed by auxiliary args that are passed verbatim for 

292 both the data and mask calls. 

293 """ 

294 def __call__(self, x, *args, **params): 

295 func = getattr(np, self.__name__) 

296 _d = func(tuple([np.asarray(a) for a in x]), *args, **params) 

297 _m = func(tuple([getmaskarray(a) for a in x]), *args, **params) 

298 return masked_array(_d, mask=_m) 

299 

300 

301class _fromnxfunction_args(_fromnxfunction): 

302 """ 

303 A version of `_fromnxfunction` that is called with multiple array 

304 arguments. The first non-array-like input marks the beginning of the 

305 arguments that are passed verbatim for both the data and mask calls. 

306 Array arguments are processed independently and the results are 

307 returned in a list. If only one array is found, the return value is 

308 just the processed array instead of a list. 

309 """ 

310 def __call__(self, *args, **params): 

311 func = getattr(np, self.__name__) 

312 arrays = [] 

313 args = list(args) 

314 while len(args) > 0 and issequence(args[0]): 

315 arrays.append(args.pop(0)) 

316 res = [] 

317 for x in arrays: 

318 _d = func(np.asarray(x), *args, **params) 

319 _m = func(getmaskarray(x), *args, **params) 

320 res.append(masked_array(_d, mask=_m)) 

321 if len(arrays) == 1: 

322 return res[0] 

323 return res 

324 

325 

326class _fromnxfunction_allargs(_fromnxfunction): 

327 """ 

328 A version of `_fromnxfunction` that is called with multiple array 

329 arguments. Similar to `_fromnxfunction_args` except that all args 

330 are converted to arrays even if they are not so already. This makes 

331 it possible to process scalars as 1-D arrays. Only keyword arguments 

332 are passed through verbatim for the data and mask calls. Arrays 

333 arguments are processed independently and the results are returned 

334 in a list. If only one arg is present, the return value is just the 

335 processed array instead of a list. 

336 """ 

337 def __call__(self, *args, **params): 

338 func = getattr(np, self.__name__) 

339 res = [] 

340 for x in args: 

341 _d = func(np.asarray(x), **params) 

342 _m = func(getmaskarray(x), **params) 

343 res.append(masked_array(_d, mask=_m)) 

344 if len(args) == 1: 

345 return res[0] 

346 return res 

347 

348 

349atleast_1d = _fromnxfunction_allargs('atleast_1d') 

350atleast_2d = _fromnxfunction_allargs('atleast_2d') 

351atleast_3d = _fromnxfunction_allargs('atleast_3d') 

352 

353vstack = row_stack = _fromnxfunction_seq('vstack') 

354hstack = _fromnxfunction_seq('hstack') 

355column_stack = _fromnxfunction_seq('column_stack') 

356dstack = _fromnxfunction_seq('dstack') 

357stack = _fromnxfunction_seq('stack') 

358 

359hsplit = _fromnxfunction_single('hsplit') 

360 

361diagflat = _fromnxfunction_single('diagflat') 

362 

363 

364#####-------------------------------------------------------------------------- 

365#---- 

366#####-------------------------------------------------------------------------- 

367def flatten_inplace(seq): 

368 """Flatten a sequence in place.""" 

369 k = 0 

370 while (k != len(seq)): 

371 while hasattr(seq[k], '__iter__'): 

372 seq[k:(k + 1)] = seq[k] 

373 k += 1 

374 return seq 

375 

376 

377def apply_along_axis(func1d, axis, arr, *args, **kwargs): 

378 """ 

379 (This docstring should be overwritten) 

380 """ 

381 arr = array(arr, copy=False, subok=True) 

382 nd = arr.ndim 

383 axis = normalize_axis_index(axis, nd) 

384 ind = [0] * (nd - 1) 

385 i = np.zeros(nd, 'O') 

386 indlist = list(range(nd)) 

387 indlist.remove(axis) 

388 i[axis] = slice(None, None) 

389 outshape = np.asarray(arr.shape).take(indlist) 

390 i.put(indlist, ind) 

391 res = func1d(arr[tuple(i.tolist())], *args, **kwargs) 

392 # if res is a number, then we have a smaller output array 

393 asscalar = np.isscalar(res) 

394 if not asscalar: 

395 try: 

396 len(res) 

397 except TypeError: 

398 asscalar = True 

399 # Note: we shouldn't set the dtype of the output from the first result 

400 # so we force the type to object, and build a list of dtypes. We'll 

401 # just take the largest, to avoid some downcasting 

402 dtypes = [] 

403 if asscalar: 

404 dtypes.append(np.asarray(res).dtype) 

405 outarr = zeros(outshape, object) 

406 outarr[tuple(ind)] = res 

407 Ntot = np.product(outshape) 

408 k = 1 

409 while k < Ntot: 

410 # increment the index 

411 ind[-1] += 1 

412 n = -1 

413 while (ind[n] >= outshape[n]) and (n > (1 - nd)): 

414 ind[n - 1] += 1 

415 ind[n] = 0 

416 n -= 1 

417 i.put(indlist, ind) 

418 res = func1d(arr[tuple(i.tolist())], *args, **kwargs) 

419 outarr[tuple(ind)] = res 

420 dtypes.append(asarray(res).dtype) 

421 k += 1 

422 else: 

423 res = array(res, copy=False, subok=True) 

424 j = i.copy() 

425 j[axis] = ([slice(None, None)] * res.ndim) 

426 j.put(indlist, ind) 

427 Ntot = np.product(outshape) 

428 holdshape = outshape 

429 outshape = list(arr.shape) 

430 outshape[axis] = res.shape 

431 dtypes.append(asarray(res).dtype) 

432 outshape = flatten_inplace(outshape) 

433 outarr = zeros(outshape, object) 

434 outarr[tuple(flatten_inplace(j.tolist()))] = res 

435 k = 1 

436 while k < Ntot: 

437 # increment the index 

438 ind[-1] += 1 

439 n = -1 

440 while (ind[n] >= holdshape[n]) and (n > (1 - nd)): 

441 ind[n - 1] += 1 

442 ind[n] = 0 

443 n -= 1 

444 i.put(indlist, ind) 

445 j.put(indlist, ind) 

446 res = func1d(arr[tuple(i.tolist())], *args, **kwargs) 

447 outarr[tuple(flatten_inplace(j.tolist()))] = res 

448 dtypes.append(asarray(res).dtype) 

449 k += 1 

450 max_dtypes = np.dtype(np.asarray(dtypes).max()) 

451 if not hasattr(arr, '_mask'): 

452 result = np.asarray(outarr, dtype=max_dtypes) 

453 else: 

454 result = asarray(outarr, dtype=max_dtypes) 

455 result.fill_value = ma.default_fill_value(result) 

456 return result 

457apply_along_axis.__doc__ = np.apply_along_axis.__doc__ 

458 

459 

460def apply_over_axes(func, a, axes): 

461 """ 

462 (This docstring will be overwritten) 

463 """ 

464 val = asarray(a) 

465 N = a.ndim 

466 if array(axes).ndim == 0: 

467 axes = (axes,) 

468 for axis in axes: 

469 if axis < 0: 

470 axis = N + axis 

471 args = (val, axis) 

472 res = func(*args) 

473 if res.ndim == val.ndim: 

474 val = res 

475 else: 

476 res = ma.expand_dims(res, axis) 

477 if res.ndim == val.ndim: 

478 val = res 

479 else: 

480 raise ValueError("function is not returning " 

481 "an array of the correct shape") 

482 return val 

483 

484if apply_over_axes.__doc__ is not None: 

485 apply_over_axes.__doc__ = np.apply_over_axes.__doc__[ 

486 :np.apply_over_axes.__doc__.find('Notes')].rstrip() + \ 

487 """ 

488 

489 Examples 

490 -------- 

491 >>> a = np.ma.arange(24).reshape(2,3,4) 

492 >>> a[:,0,1] = np.ma.masked 

493 >>> a[:,1,:] = np.ma.masked 

494 >>> a 

495 masked_array( 

496 data=[[[0, --, 2, 3], 

497 [--, --, --, --], 

498 [8, 9, 10, 11]], 

499 [[12, --, 14, 15], 

500 [--, --, --, --], 

501 [20, 21, 22, 23]]], 

502 mask=[[[False, True, False, False], 

503 [ True, True, True, True], 

504 [False, False, False, False]], 

505 [[False, True, False, False], 

506 [ True, True, True, True], 

507 [False, False, False, False]]], 

508 fill_value=999999) 

509 >>> np.ma.apply_over_axes(np.ma.sum, a, [0,2]) 

510 masked_array( 

511 data=[[[46], 

512 [--], 

513 [124]]], 

514 mask=[[[False], 

515 [ True], 

516 [False]]], 

517 fill_value=999999) 

518 

519 Tuple axis arguments to ufuncs are equivalent: 

520 

521 >>> np.ma.sum(a, axis=(0,2)).reshape((1,-1,1)) 

522 masked_array( 

523 data=[[[46], 

524 [--], 

525 [124]]], 

526 mask=[[[False], 

527 [ True], 

528 [False]]], 

529 fill_value=999999) 

530 """ 

531 

532 

533def average(a, axis=None, weights=None, returned=False): 

534 """ 

535 Return the weighted average of array over the given axis. 

536 

537 Parameters 

538 ---------- 

539 a : array_like 

540 Data to be averaged. 

541 Masked entries are not taken into account in the computation. 

542 axis : int, optional 

543 Axis along which to average `a`. If None, averaging is done over 

544 the flattened array. 

545 weights : array_like, optional 

546 The importance that each element has in the computation of the average. 

547 The weights array can either be 1-D (in which case its length must be 

548 the size of `a` along the given axis) or of the same shape as `a`. 

549 If ``weights=None``, then all data in `a` are assumed to have a 

550 weight equal to one. The 1-D calculation is:: 

551 

552 avg = sum(a * weights) / sum(weights) 

553 

554 The only constraint on `weights` is that `sum(weights)` must not be 0. 

555 returned : bool, optional 

556 Flag indicating whether a tuple ``(result, sum of weights)`` 

557 should be returned as output (True), or just the result (False). 

558 Default is False. 

559 

560 Returns 

561 ------- 

562 average, [sum_of_weights] : (tuple of) scalar or MaskedArray 

563 The average along the specified axis. When returned is `True`, 

564 return a tuple with the average as the first element and the sum 

565 of the weights as the second element. The return type is `np.float64` 

566 if `a` is of integer type and floats smaller than `float64`, or the 

567 input data-type, otherwise. If returned, `sum_of_weights` is always 

568 `float64`. 

569 

570 Examples 

571 -------- 

572 >>> a = np.ma.array([1., 2., 3., 4.], mask=[False, False, True, True]) 

573 >>> np.ma.average(a, weights=[3, 1, 0, 0]) 

574 1.25 

575 

576 >>> x = np.ma.arange(6.).reshape(3, 2) 

577 >>> x 

578 masked_array( 

579 data=[[0., 1.], 

580 [2., 3.], 

581 [4., 5.]], 

582 mask=False, 

583 fill_value=1e+20) 

584 >>> avg, sumweights = np.ma.average(x, axis=0, weights=[1, 2, 3], 

585 ... returned=True) 

586 >>> avg 

587 masked_array(data=[2.6666666666666665, 3.6666666666666665], 

588 mask=[False, False], 

589 fill_value=1e+20) 

590 

591 """ 

592 a = asarray(a) 

593 m = getmask(a) 

594 

595 # inspired by 'average' in numpy/lib/function_base.py 

596 

597 if weights is None: 

598 avg = a.mean(axis) 

599 scl = avg.dtype.type(a.count(axis)) 

600 else: 

601 wgt = np.asanyarray(weights) 

602 

603 if issubclass(a.dtype.type, (np.integer, np.bool_)): 

604 result_dtype = np.result_type(a.dtype, wgt.dtype, 'f8') 

605 else: 

606 result_dtype = np.result_type(a.dtype, wgt.dtype) 

607 

608 # Sanity checks 

609 if a.shape != wgt.shape: 

610 if axis is None: 

611 raise TypeError( 

612 "Axis must be specified when shapes of a and weights " 

613 "differ.") 

614 if wgt.ndim != 1: 

615 raise TypeError( 

616 "1D weights expected when shapes of a and weights differ.") 

617 if wgt.shape[0] != a.shape[axis]: 

618 raise ValueError( 

619 "Length of weights not compatible with specified axis.") 

620 

621 # setup wgt to broadcast along axis 

622 wgt = np.broadcast_to(wgt, (a.ndim-1)*(1,) + wgt.shape) 

623 wgt = wgt.swapaxes(-1, axis) 

624 

625 if m is not nomask: 

626 wgt = wgt*(~a.mask) 

627 

628 scl = wgt.sum(axis=axis, dtype=result_dtype) 

629 avg = np.multiply(a, wgt, dtype=result_dtype).sum(axis)/scl 

630 

631 if returned: 

632 if scl.shape != avg.shape: 

633 scl = np.broadcast_to(scl, avg.shape).copy() 

634 return avg, scl 

635 else: 

636 return avg 

637 

638 

639def median(a, axis=None, out=None, overwrite_input=False, keepdims=False): 

640 """ 

641 Compute the median along the specified axis. 

642 

643 Returns the median of the array elements. 

644 

645 Parameters 

646 ---------- 

647 a : array_like 

648 Input array or object that can be converted to an array. 

649 axis : int, optional 

650 Axis along which the medians are computed. The default (None) is 

651 to compute the median along a flattened version of the array. 

652 out : ndarray, optional 

653 Alternative output array in which to place the result. It must 

654 have the same shape and buffer length as the expected output 

655 but the type will be cast if necessary. 

656 overwrite_input : bool, optional 

657 If True, then allow use of memory of input array (a) for 

658 calculations. The input array will be modified by the call to 

659 median. This will save memory when you do not need to preserve 

660 the contents of the input array. Treat the input as undefined, 

661 but it will probably be fully or partially sorted. Default is 

662 False. Note that, if `overwrite_input` is True, and the input 

663 is not already an `ndarray`, an error will be raised. 

664 keepdims : bool, optional 

665 If this is set to True, the axes which are reduced are left 

666 in the result as dimensions with size one. With this option, 

667 the result will broadcast correctly against the input array. 

668 

669 .. versionadded:: 1.10.0 

670 

671 Returns 

672 ------- 

673 median : ndarray 

674 A new array holding the result is returned unless out is 

675 specified, in which case a reference to out is returned. 

676 Return data-type is `float64` for integers and floats smaller than 

677 `float64`, or the input data-type, otherwise. 

678 

679 See Also 

680 -------- 

681 mean 

682 

683 Notes 

684 ----- 

685 Given a vector ``V`` with ``N`` non masked values, the median of ``V`` 

686 is the middle value of a sorted copy of ``V`` (``Vs``) - i.e. 

687 ``Vs[(N-1)/2]``, when ``N`` is odd, or ``{Vs[N/2 - 1] + Vs[N/2]}/2`` 

688 when ``N`` is even. 

689 

690 Examples 

691 -------- 

692 >>> x = np.ma.array(np.arange(8), mask=[0]*4 + [1]*4) 

693 >>> np.ma.median(x) 

694 1.5 

695 

696 >>> x = np.ma.array(np.arange(10).reshape(2, 5), mask=[0]*6 + [1]*4) 

697 >>> np.ma.median(x) 

698 2.5 

699 >>> np.ma.median(x, axis=-1, overwrite_input=True) 

700 masked_array(data=[2.0, 5.0], 

701 mask=[False, False], 

702 fill_value=1e+20) 

703 

704 """ 

705 if not hasattr(a, 'mask'): 

706 m = np.median(getdata(a, subok=True), axis=axis, 

707 out=out, overwrite_input=overwrite_input, 

708 keepdims=keepdims) 

709 if isinstance(m, np.ndarray) and 1 <= m.ndim: 

710 return masked_array(m, copy=False) 

711 else: 

712 return m 

713 

714 r, k = _ureduce(a, func=_median, axis=axis, out=out, 

715 overwrite_input=overwrite_input) 

716 if keepdims: 

717 return r.reshape(k) 

718 else: 

719 return r 

720 

721def _median(a, axis=None, out=None, overwrite_input=False): 

722 # when an unmasked NaN is present return it, so we need to sort the NaN 

723 # values behind the mask 

724 if np.issubdtype(a.dtype, np.inexact): 

725 fill_value = np.inf 

726 else: 

727 fill_value = None 

728 if overwrite_input: 

729 if axis is None: 

730 asorted = a.ravel() 

731 asorted.sort(fill_value=fill_value) 

732 else: 

733 a.sort(axis=axis, fill_value=fill_value) 

734 asorted = a 

735 else: 

736 asorted = sort(a, axis=axis, fill_value=fill_value) 

737 

738 if axis is None: 

739 axis = 0 

740 else: 

741 axis = normalize_axis_index(axis, asorted.ndim) 

742 

743 if asorted.shape[axis] == 0: 

744 # for empty axis integer indices fail so use slicing to get same result 

745 # as median (which is mean of empty slice = nan) 

746 indexer = [slice(None)] * asorted.ndim 

747 indexer[axis] = slice(0, 0) 

748 indexer = tuple(indexer) 

749 return np.ma.mean(asorted[indexer], axis=axis, out=out) 

750 

751 if asorted.ndim == 1: 

752 counts = count(asorted) 

753 idx, odd = divmod(count(asorted), 2) 

754 mid = asorted[idx + odd - 1:idx + 1] 

755 if np.issubdtype(asorted.dtype, np.inexact) and asorted.size > 0: 

756 # avoid inf / x = masked 

757 s = mid.sum(out=out) 

758 if not odd: 

759 s = np.true_divide(s, 2., casting='safe', out=out) 

760 s = np.lib.utils._median_nancheck(asorted, s, axis, out) 

761 else: 

762 s = mid.mean(out=out) 

763 

764 # if result is masked either the input contained enough 

765 # minimum_fill_value so that it would be the median or all values 

766 # masked 

767 if np.ma.is_masked(s) and not np.all(asorted.mask): 

768 return np.ma.minimum_fill_value(asorted) 

769 return s 

770 

771 counts = count(asorted, axis=axis, keepdims=True) 

772 h = counts // 2 

773 

774 # duplicate high if odd number of elements so mean does nothing 

775 odd = counts % 2 == 1 

776 l = np.where(odd, h, h-1) 

777 

778 lh = np.concatenate([l,h], axis=axis) 

779 

780 # get low and high median 

781 low_high = np.take_along_axis(asorted, lh, axis=axis) 

782 

783 def replace_masked(s): 

784 # Replace masked entries with minimum_full_value unless it all values 

785 # are masked. This is required as the sort order of values equal or 

786 # larger than the fill value is undefined and a valid value placed 

787 # elsewhere, e.g. [4, --, inf]. 

788 if np.ma.is_masked(s): 

789 rep = (~np.all(asorted.mask, axis=axis, keepdims=True)) & s.mask 

790 s.data[rep] = np.ma.minimum_fill_value(asorted) 

791 s.mask[rep] = False 

792 

793 replace_masked(low_high) 

794 

795 if np.issubdtype(asorted.dtype, np.inexact): 

796 # avoid inf / x = masked 

797 s = np.ma.sum(low_high, axis=axis, out=out) 

798 np.true_divide(s.data, 2., casting='unsafe', out=s.data) 

799 

800 s = np.lib.utils._median_nancheck(asorted, s, axis, out) 

801 else: 

802 s = np.ma.mean(low_high, axis=axis, out=out) 

803 

804 return s 

805 

806 

807def compress_nd(x, axis=None): 

808 """Suppress slices from multiple dimensions which contain masked values. 

809 

810 Parameters 

811 ---------- 

812 x : array_like, MaskedArray 

813 The array to operate on. If not a MaskedArray instance (or if no array 

814 elements are masked), `x` is interpreted as a MaskedArray with `mask` 

815 set to `nomask`. 

816 axis : tuple of ints or int, optional 

817 Which dimensions to suppress slices from can be configured with this 

818 parameter. 

819 - If axis is a tuple of ints, those are the axes to suppress slices from. 

820 - If axis is an int, then that is the only axis to suppress slices from. 

821 - If axis is None, all axis are selected. 

822 

823 Returns 

824 ------- 

825 compress_array : ndarray 

826 The compressed array. 

827 """ 

828 x = asarray(x) 

829 m = getmask(x) 

830 # Set axis to tuple of ints 

831 if axis is None: 

832 axis = tuple(range(x.ndim)) 

833 else: 

834 axis = normalize_axis_tuple(axis, x.ndim) 

835 

836 # Nothing is masked: return x 

837 if m is nomask or not m.any(): 

838 return x._data 

839 # All is masked: return empty 

840 if m.all(): 

841 return nxarray([]) 

842 # Filter elements through boolean indexing 

843 data = x._data 

844 for ax in axis: 

845 axes = tuple(list(range(ax)) + list(range(ax + 1, x.ndim))) 

846 data = data[(slice(None),)*ax + (~m.any(axis=axes),)] 

847 return data 

848 

849def compress_rowcols(x, axis=None): 

850 """ 

851 Suppress the rows and/or columns of a 2-D array that contain 

852 masked values. 

853 

854 The suppression behavior is selected with the `axis` parameter. 

855 

856 - If axis is None, both rows and columns are suppressed. 

857 - If axis is 0, only rows are suppressed. 

858 - If axis is 1 or -1, only columns are suppressed. 

859 

860 Parameters 

861 ---------- 

862 x : array_like, MaskedArray 

863 The array to operate on. If not a MaskedArray instance (or if no array 

864 elements are masked), `x` is interpreted as a MaskedArray with 

865 `mask` set to `nomask`. Must be a 2D array. 

866 axis : int, optional 

867 Axis along which to perform the operation. Default is None. 

868 

869 Returns 

870 ------- 

871 compressed_array : ndarray 

872 The compressed array. 

873 

874 Examples 

875 -------- 

876 >>> x = np.ma.array(np.arange(9).reshape(3, 3), mask=[[1, 0, 0], 

877 ... [1, 0, 0], 

878 ... [0, 0, 0]]) 

879 >>> x 

880 masked_array( 

881 data=[[--, 1, 2], 

882 [--, 4, 5], 

883 [6, 7, 8]], 

884 mask=[[ True, False, False], 

885 [ True, False, False], 

886 [False, False, False]], 

887 fill_value=999999) 

888 

889 >>> np.ma.compress_rowcols(x) 

890 array([[7, 8]]) 

891 >>> np.ma.compress_rowcols(x, 0) 

892 array([[6, 7, 8]]) 

893 >>> np.ma.compress_rowcols(x, 1) 

894 array([[1, 2], 

895 [4, 5], 

896 [7, 8]]) 

897 

898 """ 

899 if asarray(x).ndim != 2: 

900 raise NotImplementedError("compress_rowcols works for 2D arrays only.") 

901 return compress_nd(x, axis=axis) 

902 

903 

904def compress_rows(a): 

905 """ 

906 Suppress whole rows of a 2-D array that contain masked values. 

907 

908 This is equivalent to ``np.ma.compress_rowcols(a, 0)``, see 

909 `extras.compress_rowcols` for details. 

910 

911 See Also 

912 -------- 

913 extras.compress_rowcols 

914 

915 """ 

916 a = asarray(a) 

917 if a.ndim != 2: 

918 raise NotImplementedError("compress_rows works for 2D arrays only.") 

919 return compress_rowcols(a, 0) 

920 

921def compress_cols(a): 

922 """ 

923 Suppress whole columns of a 2-D array that contain masked values. 

924 

925 This is equivalent to ``np.ma.compress_rowcols(a, 1)``, see 

926 `extras.compress_rowcols` for details. 

927 

928 See Also 

929 -------- 

930 extras.compress_rowcols 

931 

932 """ 

933 a = asarray(a) 

934 if a.ndim != 2: 

935 raise NotImplementedError("compress_cols works for 2D arrays only.") 

936 return compress_rowcols(a, 1) 

937 

938def mask_rows(a, axis=np._NoValue): 

939 """ 

940 Mask rows of a 2D array that contain masked values. 

941 

942 This function is a shortcut to ``mask_rowcols`` with `axis` equal to 0. 

943 

944 See Also 

945 -------- 

946 mask_rowcols : Mask rows and/or columns of a 2D array. 

947 masked_where : Mask where a condition is met. 

948 

949 Examples 

950 -------- 

951 >>> import numpy.ma as ma 

952 >>> a = np.zeros((3, 3), dtype=int) 

953 >>> a[1, 1] = 1 

954 >>> a 

955 array([[0, 0, 0], 

956 [0, 1, 0], 

957 [0, 0, 0]]) 

958 >>> a = ma.masked_equal(a, 1) 

959 >>> a 

960 masked_array( 

961 data=[[0, 0, 0], 

962 [0, --, 0], 

963 [0, 0, 0]], 

964 mask=[[False, False, False], 

965 [False, True, False], 

966 [False, False, False]], 

967 fill_value=1) 

968 

969 >>> ma.mask_rows(a) 

970 masked_array( 

971 data=[[0, 0, 0], 

972 [--, --, --], 

973 [0, 0, 0]], 

974 mask=[[False, False, False], 

975 [ True, True, True], 

976 [False, False, False]], 

977 fill_value=1) 

978 

979 """ 

980 if axis is not np._NoValue: 

981 # remove the axis argument when this deprecation expires 

982 # NumPy 1.18.0, 2019-11-28 

983 warnings.warn( 

984 "The axis argument has always been ignored, in future passing it " 

985 "will raise TypeError", DeprecationWarning, stacklevel=2) 

986 return mask_rowcols(a, 0) 

987 

988def mask_cols(a, axis=np._NoValue): 

989 """ 

990 Mask columns of a 2D array that contain masked values. 

991 

992 This function is a shortcut to ``mask_rowcols`` with `axis` equal to 1. 

993 

994 See Also 

995 -------- 

996 mask_rowcols : Mask rows and/or columns of a 2D array. 

997 masked_where : Mask where a condition is met. 

998 

999 Examples 

1000 -------- 

1001 >>> import numpy.ma as ma 

1002 >>> a = np.zeros((3, 3), dtype=int) 

1003 >>> a[1, 1] = 1 

1004 >>> a 

1005 array([[0, 0, 0], 

1006 [0, 1, 0], 

1007 [0, 0, 0]]) 

1008 >>> a = ma.masked_equal(a, 1) 

1009 >>> a 

1010 masked_array( 

1011 data=[[0, 0, 0], 

1012 [0, --, 0], 

1013 [0, 0, 0]], 

1014 mask=[[False, False, False], 

1015 [False, True, False], 

1016 [False, False, False]], 

1017 fill_value=1) 

1018 >>> ma.mask_cols(a) 

1019 masked_array( 

1020 data=[[0, --, 0], 

1021 [0, --, 0], 

1022 [0, --, 0]], 

1023 mask=[[False, True, False], 

1024 [False, True, False], 

1025 [False, True, False]], 

1026 fill_value=1) 

1027 

1028 """ 

1029 if axis is not np._NoValue: 

1030 # remove the axis argument when this deprecation expires 

1031 # NumPy 1.18.0, 2019-11-28 

1032 warnings.warn( 

1033 "The axis argument has always been ignored, in future passing it " 

1034 "will raise TypeError", DeprecationWarning, stacklevel=2) 

1035 return mask_rowcols(a, 1) 

1036 

1037 

1038#####-------------------------------------------------------------------------- 

1039#---- --- arraysetops --- 

1040#####-------------------------------------------------------------------------- 

1041 

1042def ediff1d(arr, to_end=None, to_begin=None): 

1043 """ 

1044 Compute the differences between consecutive elements of an array. 

1045 

1046 This function is the equivalent of `numpy.ediff1d` that takes masked 

1047 values into account, see `numpy.ediff1d` for details. 

1048 

1049 See Also 

1050 -------- 

1051 numpy.ediff1d : Equivalent function for ndarrays. 

1052 

1053 """ 

1054 arr = ma.asanyarray(arr).flat 

1055 ed = arr[1:] - arr[:-1] 

1056 arrays = [ed] 

1057 # 

1058 if to_begin is not None: 

1059 arrays.insert(0, to_begin) 

1060 if to_end is not None: 

1061 arrays.append(to_end) 

1062 # 

1063 if len(arrays) != 1: 

1064 # We'll save ourselves a copy of a potentially large array in the common 

1065 # case where neither to_begin or to_end was given. 

1066 ed = hstack(arrays) 

1067 # 

1068 return ed 

1069 

1070 

1071def unique(ar1, return_index=False, return_inverse=False): 

1072 """ 

1073 Finds the unique elements of an array. 

1074 

1075 Masked values are considered the same element (masked). The output array 

1076 is always a masked array. See `numpy.unique` for more details. 

1077 

1078 See Also 

1079 -------- 

1080 numpy.unique : Equivalent function for ndarrays. 

1081 

1082 """ 

1083 output = np.unique(ar1, 

1084 return_index=return_index, 

1085 return_inverse=return_inverse) 

1086 if isinstance(output, tuple): 

1087 output = list(output) 

1088 output[0] = output[0].view(MaskedArray) 

1089 output = tuple(output) 

1090 else: 

1091 output = output.view(MaskedArray) 

1092 return output 

1093 

1094 

1095def intersect1d(ar1, ar2, assume_unique=False): 

1096 """ 

1097 Returns the unique elements common to both arrays. 

1098 

1099 Masked values are considered equal one to the other. 

1100 The output is always a masked array. 

1101 

1102 See `numpy.intersect1d` for more details. 

1103 

1104 See Also 

1105 -------- 

1106 numpy.intersect1d : Equivalent function for ndarrays. 

1107 

1108 Examples 

1109 -------- 

1110 >>> x = np.ma.array([1, 3, 3, 3], mask=[0, 0, 0, 1]) 

1111 >>> y = np.ma.array([3, 1, 1, 1], mask=[0, 0, 0, 1]) 

1112 >>> np.ma.intersect1d(x, y) 

1113 masked_array(data=[1, 3, --], 

1114 mask=[False, False, True], 

1115 fill_value=999999) 

1116 

1117 """ 

1118 if assume_unique: 

1119 aux = ma.concatenate((ar1, ar2)) 

1120 else: 

1121 # Might be faster than unique( intersect1d( ar1, ar2 ) )? 

1122 aux = ma.concatenate((unique(ar1), unique(ar2))) 

1123 aux.sort() 

1124 return aux[:-1][aux[1:] == aux[:-1]] 

1125 

1126 

1127def setxor1d(ar1, ar2, assume_unique=False): 

1128 """ 

1129 Set exclusive-or of 1-D arrays with unique elements. 

1130 

1131 The output is always a masked array. See `numpy.setxor1d` for more details. 

1132 

1133 See Also 

1134 -------- 

1135 numpy.setxor1d : Equivalent function for ndarrays. 

1136 

1137 """ 

1138 if not assume_unique: 

1139 ar1 = unique(ar1) 

1140 ar2 = unique(ar2) 

1141 

1142 aux = ma.concatenate((ar1, ar2)) 

1143 if aux.size == 0: 

1144 return aux 

1145 aux.sort() 

1146 auxf = aux.filled() 

1147# flag = ediff1d( aux, to_end = 1, to_begin = 1 ) == 0 

1148 flag = ma.concatenate(([True], (auxf[1:] != auxf[:-1]), [True])) 

1149# flag2 = ediff1d( flag ) == 0 

1150 flag2 = (flag[1:] == flag[:-1]) 

1151 return aux[flag2] 

1152 

1153 

1154def in1d(ar1, ar2, assume_unique=False, invert=False): 

1155 """ 

1156 Test whether each element of an array is also present in a second 

1157 array. 

1158 

1159 The output is always a masked array. See `numpy.in1d` for more details. 

1160 

1161 We recommend using :func:`isin` instead of `in1d` for new code. 

1162 

1163 See Also 

1164 -------- 

1165 isin : Version of this function that preserves the shape of ar1. 

1166 numpy.in1d : Equivalent function for ndarrays. 

1167 

1168 Notes 

1169 ----- 

1170 .. versionadded:: 1.4.0 

1171 

1172 """ 

1173 if not assume_unique: 

1174 ar1, rev_idx = unique(ar1, return_inverse=True) 

1175 ar2 = unique(ar2) 

1176 

1177 ar = ma.concatenate((ar1, ar2)) 

1178 # We need this to be a stable sort, so always use 'mergesort' 

1179 # here. The values from the first array should always come before 

1180 # the values from the second array. 

1181 order = ar.argsort(kind='mergesort') 

1182 sar = ar[order] 

1183 if invert: 

1184 bool_ar = (sar[1:] != sar[:-1]) 

1185 else: 

1186 bool_ar = (sar[1:] == sar[:-1]) 

1187 flag = ma.concatenate((bool_ar, [invert])) 

1188 indx = order.argsort(kind='mergesort')[:len(ar1)] 

1189 

1190 if assume_unique: 

1191 return flag[indx] 

1192 else: 

1193 return flag[indx][rev_idx] 

1194 

1195 

1196def isin(element, test_elements, assume_unique=False, invert=False): 

1197 """ 

1198 Calculates `element in test_elements`, broadcasting over 

1199 `element` only. 

1200 

1201 The output is always a masked array of the same shape as `element`. 

1202 See `numpy.isin` for more details. 

1203 

1204 See Also 

1205 -------- 

1206 in1d : Flattened version of this function. 

1207 numpy.isin : Equivalent function for ndarrays. 

1208 

1209 Notes 

1210 ----- 

1211 .. versionadded:: 1.13.0 

1212 

1213 """ 

1214 element = ma.asarray(element) 

1215 return in1d(element, test_elements, assume_unique=assume_unique, 

1216 invert=invert).reshape(element.shape) 

1217 

1218 

1219def union1d(ar1, ar2): 

1220 """ 

1221 Union of two arrays. 

1222 

1223 The output is always a masked array. See `numpy.union1d` for more details. 

1224 

1225 See also 

1226 -------- 

1227 numpy.union1d : Equivalent function for ndarrays. 

1228 

1229 """ 

1230 return unique(ma.concatenate((ar1, ar2), axis=None)) 

1231 

1232 

1233def setdiff1d(ar1, ar2, assume_unique=False): 

1234 """ 

1235 Set difference of 1D arrays with unique elements. 

1236 

1237 The output is always a masked array. See `numpy.setdiff1d` for more 

1238 details. 

1239 

1240 See Also 

1241 -------- 

1242 numpy.setdiff1d : Equivalent function for ndarrays. 

1243 

1244 Examples 

1245 -------- 

1246 >>> x = np.ma.array([1, 2, 3, 4], mask=[0, 1, 0, 1]) 

1247 >>> np.ma.setdiff1d(x, [1, 2]) 

1248 masked_array(data=[3, --], 

1249 mask=[False, True], 

1250 fill_value=999999) 

1251 

1252 """ 

1253 if assume_unique: 

1254 ar1 = ma.asarray(ar1).ravel() 

1255 else: 

1256 ar1 = unique(ar1) 

1257 ar2 = unique(ar2) 

1258 return ar1[in1d(ar1, ar2, assume_unique=True, invert=True)] 

1259 

1260 

1261############################################################################### 

1262# Covariance # 

1263############################################################################### 

1264 

1265 

1266def _covhelper(x, y=None, rowvar=True, allow_masked=True): 

1267 """ 

1268 Private function for the computation of covariance and correlation 

1269 coefficients. 

1270 

1271 """ 

1272 x = ma.array(x, ndmin=2, copy=True, dtype=float) 

1273 xmask = ma.getmaskarray(x) 

1274 # Quick exit if we can't process masked data 

1275 if not allow_masked and xmask.any(): 

1276 raise ValueError("Cannot process masked data.") 

1277 # 

1278 if x.shape[0] == 1: 

1279 rowvar = True 

1280 # Make sure that rowvar is either 0 or 1 

1281 rowvar = int(bool(rowvar)) 

1282 axis = 1 - rowvar 

1283 if rowvar: 

1284 tup = (slice(None), None) 

1285 else: 

1286 tup = (None, slice(None)) 

1287 # 

1288 if y is None: 

1289 xnotmask = np.logical_not(xmask).astype(int) 

1290 else: 

1291 y = array(y, copy=False, ndmin=2, dtype=float) 

1292 ymask = ma.getmaskarray(y) 

1293 if not allow_masked and ymask.any(): 

1294 raise ValueError("Cannot process masked data.") 

1295 if xmask.any() or ymask.any(): 

1296 if y.shape == x.shape: 

1297 # Define some common mask 

1298 common_mask = np.logical_or(xmask, ymask) 

1299 if common_mask is not nomask: 

1300 xmask = x._mask = y._mask = ymask = common_mask 

1301 x._sharedmask = False 

1302 y._sharedmask = False 

1303 x = ma.concatenate((x, y), axis) 

1304 xnotmask = np.logical_not(np.concatenate((xmask, ymask), axis)).astype(int) 

1305 x -= x.mean(axis=rowvar)[tup] 

1306 return (x, xnotmask, rowvar) 

1307 

1308 

1309def cov(x, y=None, rowvar=True, bias=False, allow_masked=True, ddof=None): 

1310 """ 

1311 Estimate the covariance matrix. 

1312 

1313 Except for the handling of missing data this function does the same as 

1314 `numpy.cov`. For more details and examples, see `numpy.cov`. 

1315 

1316 By default, masked values are recognized as such. If `x` and `y` have the 

1317 same shape, a common mask is allocated: if ``x[i,j]`` is masked, then 

1318 ``y[i,j]`` will also be masked. 

1319 Setting `allow_masked` to False will raise an exception if values are 

1320 missing in either of the input arrays. 

1321 

1322 Parameters 

1323 ---------- 

1324 x : array_like 

1325 A 1-D or 2-D array containing multiple variables and observations. 

1326 Each row of `x` represents a variable, and each column a single 

1327 observation of all those variables. Also see `rowvar` below. 

1328 y : array_like, optional 

1329 An additional set of variables and observations. `y` has the same 

1330 form as `x`. 

1331 rowvar : bool, optional 

1332 If `rowvar` is True (default), then each row represents a 

1333 variable, with observations in the columns. Otherwise, the relationship 

1334 is transposed: each column represents a variable, while the rows 

1335 contain observations. 

1336 bias : bool, optional 

1337 Default normalization (False) is by ``(N-1)``, where ``N`` is the 

1338 number of observations given (unbiased estimate). If `bias` is True, 

1339 then normalization is by ``N``. This keyword can be overridden by 

1340 the keyword ``ddof`` in numpy versions >= 1.5. 

1341 allow_masked : bool, optional 

1342 If True, masked values are propagated pair-wise: if a value is masked 

1343 in `x`, the corresponding value is masked in `y`. 

1344 If False, raises a `ValueError` exception when some values are missing. 

1345 ddof : {None, int}, optional 

1346 If not ``None`` normalization is by ``(N - ddof)``, where ``N`` is 

1347 the number of observations; this overrides the value implied by 

1348 ``bias``. The default value is ``None``. 

1349 

1350 .. versionadded:: 1.5 

1351 

1352 Raises 

1353 ------ 

1354 ValueError 

1355 Raised if some values are missing and `allow_masked` is False. 

1356 

1357 See Also 

1358 -------- 

1359 numpy.cov 

1360 

1361 """ 

1362 # Check inputs 

1363 if ddof is not None and ddof != int(ddof): 

1364 raise ValueError("ddof must be an integer") 

1365 # Set up ddof 

1366 if ddof is None: 

1367 if bias: 

1368 ddof = 0 

1369 else: 

1370 ddof = 1 

1371 

1372 (x, xnotmask, rowvar) = _covhelper(x, y, rowvar, allow_masked) 

1373 if not rowvar: 

1374 fact = np.dot(xnotmask.T, xnotmask) * 1. - ddof 

1375 result = (dot(x.T, x.conj(), strict=False) / fact).squeeze() 

1376 else: 

1377 fact = np.dot(xnotmask, xnotmask.T) * 1. - ddof 

1378 result = (dot(x, x.T.conj(), strict=False) / fact).squeeze() 

1379 return result 

1380 

1381 

1382def corrcoef(x, y=None, rowvar=True, bias=np._NoValue, allow_masked=True, 

1383 ddof=np._NoValue): 

1384 """ 

1385 Return Pearson product-moment correlation coefficients. 

1386 

1387 Except for the handling of missing data this function does the same as 

1388 `numpy.corrcoef`. For more details and examples, see `numpy.corrcoef`. 

1389 

1390 Parameters 

1391 ---------- 

1392 x : array_like 

1393 A 1-D or 2-D array containing multiple variables and observations. 

1394 Each row of `x` represents a variable, and each column a single 

1395 observation of all those variables. Also see `rowvar` below. 

1396 y : array_like, optional 

1397 An additional set of variables and observations. `y` has the same 

1398 shape as `x`. 

1399 rowvar : bool, optional 

1400 If `rowvar` is True (default), then each row represents a 

1401 variable, with observations in the columns. Otherwise, the relationship 

1402 is transposed: each column represents a variable, while the rows 

1403 contain observations. 

1404 bias : _NoValue, optional 

1405 Has no effect, do not use. 

1406 

1407 .. deprecated:: 1.10.0 

1408 allow_masked : bool, optional 

1409 If True, masked values are propagated pair-wise: if a value is masked 

1410 in `x`, the corresponding value is masked in `y`. 

1411 If False, raises an exception. Because `bias` is deprecated, this 

1412 argument needs to be treated as keyword only to avoid a warning. 

1413 ddof : _NoValue, optional 

1414 Has no effect, do not use. 

1415 

1416 .. deprecated:: 1.10.0 

1417 

1418 See Also 

1419 -------- 

1420 numpy.corrcoef : Equivalent function in top-level NumPy module. 

1421 cov : Estimate the covariance matrix. 

1422 

1423 Notes 

1424 ----- 

1425 This function accepts but discards arguments `bias` and `ddof`. This is 

1426 for backwards compatibility with previous versions of this function. These 

1427 arguments had no effect on the return values of the function and can be 

1428 safely ignored in this and previous versions of numpy. 

1429 """ 

1430 msg = 'bias and ddof have no effect and are deprecated' 

1431 if bias is not np._NoValue or ddof is not np._NoValue: 

1432 # 2015-03-15, 1.10 

1433 warnings.warn(msg, DeprecationWarning, stacklevel=2) 

1434 # Get the data 

1435 (x, xnotmask, rowvar) = _covhelper(x, y, rowvar, allow_masked) 

1436 # Compute the covariance matrix 

1437 if not rowvar: 

1438 fact = np.dot(xnotmask.T, xnotmask) * 1. 

1439 c = (dot(x.T, x.conj(), strict=False) / fact).squeeze() 

1440 else: 

1441 fact = np.dot(xnotmask, xnotmask.T) * 1. 

1442 c = (dot(x, x.T.conj(), strict=False) / fact).squeeze() 

1443 # Check whether we have a scalar 

1444 try: 

1445 diag = ma.diagonal(c) 

1446 except ValueError: 

1447 return 1 

1448 # 

1449 if xnotmask.all(): 

1450 _denom = ma.sqrt(ma.multiply.outer(diag, diag)) 

1451 else: 

1452 _denom = diagflat(diag) 

1453 _denom._sharedmask = False # We know return is always a copy 

1454 n = x.shape[1 - rowvar] 

1455 if rowvar: 

1456 for i in range(n - 1): 

1457 for j in range(i + 1, n): 

1458 _x = mask_cols(vstack((x[i], x[j]))).var(axis=1) 

1459 _denom[i, j] = _denom[j, i] = ma.sqrt(ma.multiply.reduce(_x)) 

1460 else: 

1461 for i in range(n - 1): 

1462 for j in range(i + 1, n): 

1463 _x = mask_cols( 

1464 vstack((x[:, i], x[:, j]))).var(axis=1) 

1465 _denom[i, j] = _denom[j, i] = ma.sqrt(ma.multiply.reduce(_x)) 

1466 return c / _denom 

1467 

1468#####-------------------------------------------------------------------------- 

1469#---- --- Concatenation helpers --- 

1470#####-------------------------------------------------------------------------- 

1471 

1472class MAxisConcatenator(AxisConcatenator): 

1473 """ 

1474 Translate slice objects to concatenation along an axis. 

1475 

1476 For documentation on usage, see `mr_class`. 

1477 

1478 See Also 

1479 -------- 

1480 mr_class 

1481 

1482 """ 

1483 concatenate = staticmethod(concatenate) 

1484 

1485 @classmethod 

1486 def makemat(cls, arr): 

1487 # There used to be a view as np.matrix here, but we may eventually 

1488 # deprecate that class. In preparation, we use the unmasked version 

1489 # to construct the matrix (with copy=False for backwards compatibility 

1490 # with the .view) 

1491 data = super(MAxisConcatenator, cls).makemat(arr.data, copy=False) 

1492 return array(data, mask=arr.mask) 

1493 

1494 def __getitem__(self, key): 

1495 # matrix builder syntax, like 'a, b; c, d' 

1496 if isinstance(key, str): 

1497 raise MAError("Unavailable for masked array.") 

1498 

1499 return super(MAxisConcatenator, self).__getitem__(key) 

1500 

1501 

1502class mr_class(MAxisConcatenator): 

1503 """ 

1504 Translate slice objects to concatenation along the first axis. 

1505 

1506 This is the masked array version of `lib.index_tricks.RClass`. 

1507 

1508 See Also 

1509 -------- 

1510 lib.index_tricks.RClass 

1511 

1512 Examples 

1513 -------- 

1514 >>> np.ma.mr_[np.ma.array([1,2,3]), 0, 0, np.ma.array([4,5,6])] 

1515 masked_array(data=[1, 2, 3, ..., 4, 5, 6], 

1516 mask=False, 

1517 fill_value=999999) 

1518 

1519 """ 

1520 def __init__(self): 

1521 MAxisConcatenator.__init__(self, 0) 

1522 

1523mr_ = mr_class() 

1524 

1525#####-------------------------------------------------------------------------- 

1526#---- Find unmasked data --- 

1527#####-------------------------------------------------------------------------- 

1528 

1529def flatnotmasked_edges(a): 

1530 """ 

1531 Find the indices of the first and last unmasked values. 

1532 

1533 Expects a 1-D `MaskedArray`, returns None if all values are masked. 

1534 

1535 Parameters 

1536 ---------- 

1537 a : array_like 

1538 Input 1-D `MaskedArray` 

1539 

1540 Returns 

1541 ------- 

1542 edges : ndarray or None 

1543 The indices of first and last non-masked value in the array. 

1544 Returns None if all values are masked. 

1545 

1546 See Also 

1547 -------- 

1548 flatnotmasked_contiguous, notmasked_contiguous, notmasked_edges 

1549 clump_masked, clump_unmasked 

1550 

1551 Notes 

1552 ----- 

1553 Only accepts 1-D arrays. 

1554 

1555 Examples 

1556 -------- 

1557 >>> a = np.ma.arange(10) 

1558 >>> np.ma.flatnotmasked_edges(a) 

1559 array([0, 9]) 

1560 

1561 >>> mask = (a < 3) | (a > 8) | (a == 5) 

1562 >>> a[mask] = np.ma.masked 

1563 >>> np.array(a[~a.mask]) 

1564 array([3, 4, 6, 7, 8]) 

1565 

1566 >>> np.ma.flatnotmasked_edges(a) 

1567 array([3, 8]) 

1568 

1569 >>> a[:] = np.ma.masked 

1570 >>> print(np.ma.flatnotmasked_edges(a)) 

1571 None 

1572 

1573 """ 

1574 m = getmask(a) 

1575 if m is nomask or not np.any(m): 

1576 return np.array([0, a.size - 1]) 

1577 unmasked = np.flatnonzero(~m) 

1578 if len(unmasked) > 0: 

1579 return unmasked[[0, -1]] 

1580 else: 

1581 return None 

1582 

1583 

1584def notmasked_edges(a, axis=None): 

1585 """ 

1586 Find the indices of the first and last unmasked values along an axis. 

1587 

1588 If all values are masked, return None. Otherwise, return a list 

1589 of two tuples, corresponding to the indices of the first and last 

1590 unmasked values respectively. 

1591 

1592 Parameters 

1593 ---------- 

1594 a : array_like 

1595 The input array. 

1596 axis : int, optional 

1597 Axis along which to perform the operation. 

1598 If None (default), applies to a flattened version of the array. 

1599 

1600 Returns 

1601 ------- 

1602 edges : ndarray or list 

1603 An array of start and end indexes if there are any masked data in 

1604 the array. If there are no masked data in the array, `edges` is a 

1605 list of the first and last index. 

1606 

1607 See Also 

1608 -------- 

1609 flatnotmasked_contiguous, flatnotmasked_edges, notmasked_contiguous 

1610 clump_masked, clump_unmasked 

1611 

1612 Examples 

1613 -------- 

1614 >>> a = np.arange(9).reshape((3, 3)) 

1615 >>> m = np.zeros_like(a) 

1616 >>> m[1:, 1:] = 1 

1617 

1618 >>> am = np.ma.array(a, mask=m) 

1619 >>> np.array(am[~am.mask]) 

1620 array([0, 1, 2, 3, 6]) 

1621 

1622 >>> np.ma.notmasked_edges(am) 

1623 array([0, 6]) 

1624 

1625 """ 

1626 a = asarray(a) 

1627 if axis is None or a.ndim == 1: 

1628 return flatnotmasked_edges(a) 

1629 m = getmaskarray(a) 

1630 idx = array(np.indices(a.shape), mask=np.asarray([m] * a.ndim)) 

1631 return [tuple([idx[i].min(axis).compressed() for i in range(a.ndim)]), 

1632 tuple([idx[i].max(axis).compressed() for i in range(a.ndim)]), ] 

1633 

1634 

1635def flatnotmasked_contiguous(a): 

1636 """ 

1637 Find contiguous unmasked data in a masked array along the given axis. 

1638 

1639 Parameters 

1640 ---------- 

1641 a : narray 

1642 The input array. 

1643 

1644 Returns 

1645 ------- 

1646 slice_list : list 

1647 A sorted sequence of `slice` objects (start index, end index). 

1648 

1649 ..versionchanged:: 1.15.0 

1650 Now returns an empty list instead of None for a fully masked array 

1651 

1652 See Also 

1653 -------- 

1654 flatnotmasked_edges, notmasked_contiguous, notmasked_edges 

1655 clump_masked, clump_unmasked 

1656 

1657 Notes 

1658 ----- 

1659 Only accepts 2-D arrays at most. 

1660 

1661 Examples 

1662 -------- 

1663 >>> a = np.ma.arange(10) 

1664 >>> np.ma.flatnotmasked_contiguous(a) 

1665 [slice(0, 10, None)] 

1666 

1667 >>> mask = (a < 3) | (a > 8) | (a == 5) 

1668 >>> a[mask] = np.ma.masked 

1669 >>> np.array(a[~a.mask]) 

1670 array([3, 4, 6, 7, 8]) 

1671 

1672 >>> np.ma.flatnotmasked_contiguous(a) 

1673 [slice(3, 5, None), slice(6, 9, None)] 

1674 >>> a[:] = np.ma.masked 

1675 >>> np.ma.flatnotmasked_contiguous(a) 

1676 [] 

1677 

1678 """ 

1679 m = getmask(a) 

1680 if m is nomask: 

1681 return [slice(0, a.size)] 

1682 i = 0 

1683 result = [] 

1684 for (k, g) in itertools.groupby(m.ravel()): 

1685 n = len(list(g)) 

1686 if not k: 

1687 result.append(slice(i, i + n)) 

1688 i += n 

1689 return result 

1690 

1691def notmasked_contiguous(a, axis=None): 

1692 """ 

1693 Find contiguous unmasked data in a masked array along the given axis. 

1694 

1695 Parameters 

1696 ---------- 

1697 a : array_like 

1698 The input array. 

1699 axis : int, optional 

1700 Axis along which to perform the operation. 

1701 If None (default), applies to a flattened version of the array, and this 

1702 is the same as `flatnotmasked_contiguous`. 

1703 

1704 Returns 

1705 ------- 

1706 endpoints : list 

1707 A list of slices (start and end indexes) of unmasked indexes 

1708 in the array. 

1709 

1710 If the input is 2d and axis is specified, the result is a list of lists. 

1711 

1712 See Also 

1713 -------- 

1714 flatnotmasked_edges, flatnotmasked_contiguous, notmasked_edges 

1715 clump_masked, clump_unmasked 

1716 

1717 Notes 

1718 ----- 

1719 Only accepts 2-D arrays at most. 

1720 

1721 Examples 

1722 -------- 

1723 >>> a = np.arange(12).reshape((3, 4)) 

1724 >>> mask = np.zeros_like(a) 

1725 >>> mask[1:, :-1] = 1; mask[0, 1] = 1; mask[-1, 0] = 0 

1726 >>> ma = np.ma.array(a, mask=mask) 

1727 >>> ma 

1728 masked_array( 

1729 data=[[0, --, 2, 3], 

1730 [--, --, --, 7], 

1731 [8, --, --, 11]], 

1732 mask=[[False, True, False, False], 

1733 [ True, True, True, False], 

1734 [False, True, True, False]], 

1735 fill_value=999999) 

1736 >>> np.array(ma[~ma.mask]) 

1737 array([ 0, 2, 3, 7, 8, 11]) 

1738 

1739 >>> np.ma.notmasked_contiguous(ma) 

1740 [slice(0, 1, None), slice(2, 4, None), slice(7, 9, None), slice(11, 12, None)] 

1741 

1742 >>> np.ma.notmasked_contiguous(ma, axis=0) 

1743 [[slice(0, 1, None), slice(2, 3, None)], [], [slice(0, 1, None)], [slice(0, 3, None)]] 

1744 

1745 >>> np.ma.notmasked_contiguous(ma, axis=1) 

1746 [[slice(0, 1, None), slice(2, 4, None)], [slice(3, 4, None)], [slice(0, 1, None), slice(3, 4, None)]] 

1747 

1748 """ 

1749 a = asarray(a) 

1750 nd = a.ndim 

1751 if nd > 2: 

1752 raise NotImplementedError("Currently limited to atmost 2D array.") 

1753 if axis is None or nd == 1: 

1754 return flatnotmasked_contiguous(a) 

1755 # 

1756 result = [] 

1757 # 

1758 other = (axis + 1) % 2 

1759 idx = [0, 0] 

1760 idx[axis] = slice(None, None) 

1761 # 

1762 for i in range(a.shape[other]): 

1763 idx[other] = i 

1764 result.append(flatnotmasked_contiguous(a[tuple(idx)])) 

1765 return result 

1766 

1767 

1768def _ezclump(mask): 

1769 """ 

1770 Finds the clumps (groups of data with the same values) for a 1D bool array. 

1771 

1772 Returns a series of slices. 

1773 """ 

1774 if mask.ndim > 1: 

1775 mask = mask.ravel() 

1776 idx = (mask[1:] ^ mask[:-1]).nonzero() 

1777 idx = idx[0] + 1 

1778 

1779 if mask[0]: 

1780 if len(idx) == 0: 

1781 return [slice(0, mask.size)] 

1782 

1783 r = [slice(0, idx[0])] 

1784 r.extend((slice(left, right) 

1785 for left, right in zip(idx[1:-1:2], idx[2::2]))) 

1786 else: 

1787 if len(idx) == 0: 

1788 return [] 

1789 

1790 r = [slice(left, right) for left, right in zip(idx[:-1:2], idx[1::2])] 

1791 

1792 if mask[-1]: 

1793 r.append(slice(idx[-1], mask.size)) 

1794 return r 

1795 

1796 

1797def clump_unmasked(a): 

1798 """ 

1799 Return list of slices corresponding to the unmasked clumps of a 1-D array. 

1800 (A "clump" is defined as a contiguous region of the array). 

1801 

1802 Parameters 

1803 ---------- 

1804 a : ndarray 

1805 A one-dimensional masked array. 

1806 

1807 Returns 

1808 ------- 

1809 slices : list of slice 

1810 The list of slices, one for each continuous region of unmasked 

1811 elements in `a`. 

1812 

1813 Notes 

1814 ----- 

1815 .. versionadded:: 1.4.0 

1816 

1817 See Also 

1818 -------- 

1819 flatnotmasked_edges, flatnotmasked_contiguous, notmasked_edges 

1820 notmasked_contiguous, clump_masked 

1821 

1822 Examples 

1823 -------- 

1824 >>> a = np.ma.masked_array(np.arange(10)) 

1825 >>> a[[0, 1, 2, 6, 8, 9]] = np.ma.masked 

1826 >>> np.ma.clump_unmasked(a) 

1827 [slice(3, 6, None), slice(7, 8, None)] 

1828 

1829 """ 

1830 mask = getattr(a, '_mask', nomask) 

1831 if mask is nomask: 

1832 return [slice(0, a.size)] 

1833 return _ezclump(~mask) 

1834 

1835 

1836def clump_masked(a): 

1837 """ 

1838 Returns a list of slices corresponding to the masked clumps of a 1-D array. 

1839 (A "clump" is defined as a contiguous region of the array). 

1840 

1841 Parameters 

1842 ---------- 

1843 a : ndarray 

1844 A one-dimensional masked array. 

1845 

1846 Returns 

1847 ------- 

1848 slices : list of slice 

1849 The list of slices, one for each continuous region of masked elements 

1850 in `a`. 

1851 

1852 Notes 

1853 ----- 

1854 .. versionadded:: 1.4.0 

1855 

1856 See Also 

1857 -------- 

1858 flatnotmasked_edges, flatnotmasked_contiguous, notmasked_edges 

1859 notmasked_contiguous, clump_unmasked 

1860 

1861 Examples 

1862 -------- 

1863 >>> a = np.ma.masked_array(np.arange(10)) 

1864 >>> a[[0, 1, 2, 6, 8, 9]] = np.ma.masked 

1865 >>> np.ma.clump_masked(a) 

1866 [slice(0, 3, None), slice(6, 7, None), slice(8, 10, None)] 

1867 

1868 """ 

1869 mask = ma.getmask(a) 

1870 if mask is nomask: 

1871 return [] 

1872 return _ezclump(mask) 

1873 

1874 

1875############################################################################### 

1876# Polynomial fit # 

1877############################################################################### 

1878 

1879 

1880def vander(x, n=None): 

1881 """ 

1882 Masked values in the input array result in rows of zeros. 

1883 

1884 """ 

1885 _vander = np.vander(x, n) 

1886 m = getmask(x) 

1887 if m is not nomask: 

1888 _vander[m] = 0 

1889 return _vander 

1890 

1891vander.__doc__ = ma.doc_note(np.vander.__doc__, vander.__doc__) 

1892 

1893 

1894def polyfit(x, y, deg, rcond=None, full=False, w=None, cov=False): 

1895 """ 

1896 Any masked values in x is propagated in y, and vice-versa. 

1897 

1898 """ 

1899 x = asarray(x) 

1900 y = asarray(y) 

1901 

1902 m = getmask(x) 

1903 if y.ndim == 1: 

1904 m = mask_or(m, getmask(y)) 

1905 elif y.ndim == 2: 

1906 my = getmask(mask_rows(y)) 

1907 if my is not nomask: 

1908 m = mask_or(m, my[:, 0]) 

1909 else: 

1910 raise TypeError("Expected a 1D or 2D array for y!") 

1911 

1912 if w is not None: 

1913 w = asarray(w) 

1914 if w.ndim != 1: 

1915 raise TypeError("expected a 1-d array for weights") 

1916 if w.shape[0] != y.shape[0]: 

1917 raise TypeError("expected w and y to have the same length") 

1918 m = mask_or(m, getmask(w)) 

1919 

1920 if m is not nomask: 

1921 not_m = ~m 

1922 if w is not None: 

1923 w = w[not_m] 

1924 return np.polyfit(x[not_m], y[not_m], deg, rcond, full, w, cov) 

1925 else: 

1926 return np.polyfit(x, y, deg, rcond, full, w, cov) 

1927 

1928polyfit.__doc__ = ma.doc_note(np.polyfit.__doc__, polyfit.__doc__)