_numpy_compat.py 30 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781
  1. """Functions copypasted from newer versions of numpy.
  2. """
  3. from __future__ import division, print_function, absolute_import
  4. import warnings
  5. from warnings import WarningMessage
  6. import re
  7. from functools import wraps
  8. import numpy as np
  9. from scipy._lib._version import NumpyVersion
  10. if NumpyVersion(np.__version__) > '1.7.0.dev':
  11. _assert_warns = np.testing.assert_warns
  12. else:
  13. def _assert_warns(warning_class, func, *args, **kw):
  14. r"""
  15. Fail unless the given callable throws the specified warning.
  16. This definition is copypasted from numpy 1.9.0.dev.
  17. The version in earlier numpy returns None.
  18. Parameters
  19. ----------
  20. warning_class : class
  21. The class defining the warning that `func` is expected to throw.
  22. func : callable
  23. The callable to test.
  24. *args : Arguments
  25. Arguments passed to `func`.
  26. **kwargs : Kwargs
  27. Keyword arguments passed to `func`.
  28. Returns
  29. -------
  30. The value returned by `func`.
  31. """
  32. with warnings.catch_warnings(record=True) as l:
  33. warnings.simplefilter('always')
  34. result = func(*args, **kw)
  35. if not len(l) > 0:
  36. raise AssertionError("No warning raised when calling %s"
  37. % func.__name__)
  38. if not l[0].category is warning_class:
  39. raise AssertionError("First warning for %s is not a "
  40. "%s( is %s)" % (func.__name__, warning_class, l[0]))
  41. return result
  42. if NumpyVersion(np.__version__) >= '1.10.0':
  43. from numpy import broadcast_to
  44. else:
  45. # Definition of `broadcast_to` from numpy 1.10.0.
  46. def _maybe_view_as_subclass(original_array, new_array):
  47. if type(original_array) is not type(new_array):
  48. # if input was an ndarray subclass and subclasses were OK,
  49. # then view the result as that subclass.
  50. new_array = new_array.view(type=type(original_array))
  51. # Since we have done something akin to a view from original_array, we
  52. # should let the subclass finalize (if it has it implemented, i.e., is
  53. # not None).
  54. if new_array.__array_finalize__:
  55. new_array.__array_finalize__(original_array)
  56. return new_array
  57. def _broadcast_to(array, shape, subok, readonly):
  58. shape = tuple(shape) if np.iterable(shape) else (shape,)
  59. array = np.array(array, copy=False, subok=subok)
  60. if not shape and array.shape:
  61. raise ValueError('cannot broadcast a non-scalar to a scalar array')
  62. if any(size < 0 for size in shape):
  63. raise ValueError('all elements of broadcast shape must be non-'
  64. 'negative')
  65. broadcast = np.nditer(
  66. (array,), flags=['multi_index', 'refs_ok', 'zerosize_ok'],
  67. op_flags=['readonly'], itershape=shape, order='C').itviews[0]
  68. result = _maybe_view_as_subclass(array, broadcast)
  69. if not readonly and array.flags.writeable:
  70. result.flags.writeable = True
  71. return result
  72. def broadcast_to(array, shape, subok=False):
  73. return _broadcast_to(array, shape, subok=subok, readonly=True)
  74. if NumpyVersion(np.__version__) >= '1.11.0':
  75. def get_randint(random_state):
  76. return random_state.randint
  77. else:
  78. # In NumPy versions previous to 1.11.0 the randint funtion and the randint
  79. # method of RandomState does only work with int32 values.
  80. def get_randint(random_state):
  81. def randint_patched(low, high, size, dtype=np.int32):
  82. low = max(low, np.iinfo(dtype).min, np.iinfo(np.int32).min)
  83. high = min(high, np.iinfo(dtype).max, np.iinfo(np.int32).max)
  84. integers = random_state.randint(low, high=high, size=size)
  85. return integers.astype(dtype, copy=False)
  86. return randint_patched
  87. if NumpyVersion(np.__version__) >= '1.9.0':
  88. from numpy import unique
  89. else:
  90. # the return_counts keyword was added in 1.9.0
  91. def unique(ar, return_index=False, return_inverse=False, return_counts=False):
  92. """
  93. Find the unique elements of an array.
  94. Returns the sorted unique elements of an array. There are three optional
  95. outputs in addition to the unique elements: the indices of the input array
  96. that give the unique values, the indices of the unique array that
  97. reconstruct the input array, and the number of times each unique value
  98. comes up in the input array.
  99. Parameters
  100. ----------
  101. ar : array_like
  102. Input array. This will be flattened if it is not already 1-D.
  103. return_index : bool, optional
  104. If True, also return the indices of `ar` that result in the unique
  105. array.
  106. return_inverse : bool, optional
  107. If True, also return the indices of the unique array that can be used
  108. to reconstruct `ar`.
  109. return_counts : bool, optional
  110. If True, also return the number of times each unique value comes up
  111. in `ar`.
  112. .. versionadded:: 1.9.0
  113. Returns
  114. -------
  115. unique : ndarray
  116. The sorted unique values.
  117. unique_indices : ndarray, optional
  118. The indices of the first occurrences of the unique values in the
  119. (flattened) original array. Only provided if `return_index` is True.
  120. unique_inverse : ndarray, optional
  121. The indices to reconstruct the (flattened) original array from the
  122. unique array. Only provided if `return_inverse` is True.
  123. unique_counts : ndarray, optional
  124. The number of times each of the unique values comes up in the
  125. original array. Only provided if `return_counts` is True.
  126. .. versionadded:: 1.9.0
  127. Notes
  128. -----
  129. Taken over from numpy 1.12.0-dev (c8408bf9c). Omitted examples,
  130. see numpy documentation for those.
  131. """
  132. ar = np.asanyarray(ar).flatten()
  133. optional_indices = return_index or return_inverse
  134. optional_returns = optional_indices or return_counts
  135. if ar.size == 0:
  136. if not optional_returns:
  137. ret = ar
  138. else:
  139. ret = (ar,)
  140. if return_index:
  141. ret += (np.empty(0, np.bool),)
  142. if return_inverse:
  143. ret += (np.empty(0, np.bool),)
  144. if return_counts:
  145. ret += (np.empty(0, np.intp),)
  146. return ret
  147. if optional_indices:
  148. perm = ar.argsort(kind='mergesort' if return_index else 'quicksort')
  149. aux = ar[perm]
  150. else:
  151. ar.sort()
  152. aux = ar
  153. flag = np.concatenate(([True], aux[1:] != aux[:-1]))
  154. if not optional_returns:
  155. ret = aux[flag]
  156. else:
  157. ret = (aux[flag],)
  158. if return_index:
  159. ret += (perm[flag],)
  160. if return_inverse:
  161. iflag = np.cumsum(flag) - 1
  162. inv_idx = np.empty(ar.shape, dtype=np.intp)
  163. inv_idx[perm] = iflag
  164. ret += (inv_idx,)
  165. if return_counts:
  166. idx = np.concatenate(np.nonzero(flag) + ([ar.size],))
  167. ret += (np.diff(idx),)
  168. return ret
  169. if NumpyVersion(np.__version__) > '1.12.0.dev':
  170. polyvalfromroots = np.polynomial.polynomial.polyvalfromroots
  171. else:
  172. def polyvalfromroots(x, r, tensor=True):
  173. r"""
  174. Evaluate a polynomial specified by its roots at points x.
  175. This function is copypasted from numpy 1.12.0.dev.
  176. If `r` is of length `N`, this function returns the value
  177. .. math:: p(x) = \prod_{n=1}^{N} (x - r_n)
  178. The parameter `x` is converted to an array only if it is a tuple or a
  179. list, otherwise it is treated as a scalar. In either case, either `x`
  180. or its elements must support multiplication and addition both with
  181. themselves and with the elements of `r`.
  182. If `r` is a 1-D array, then `p(x)` will have the same shape as `x`. If
  183. `r` is multidimensional, then the shape of the result depends on the
  184. value of `tensor`. If `tensor is ``True`` the shape will be r.shape[1:]
  185. + x.shape; that is, each polynomial is evaluated at every value of `x`.
  186. If `tensor` is ``False``, the shape will be r.shape[1:]; that is, each
  187. polynomial is evaluated only for the corresponding broadcast value of
  188. `x`. Note that scalars have shape (,).
  189. Parameters
  190. ----------
  191. x : array_like, compatible object
  192. If `x` is a list or tuple, it is converted to an ndarray, otherwise
  193. it is left unchanged and treated as a scalar. In either case, `x`
  194. or its elements must support addition and multiplication with with
  195. themselves and with the elements of `r`.
  196. r : array_like
  197. Array of roots. If `r` is multidimensional the first index is the
  198. root index, while the remaining indices enumerate multiple
  199. polynomials. For instance, in the two dimensional case the roots of
  200. each polynomial may be thought of as stored in the columns of `r`.
  201. tensor : boolean, optional
  202. If True, the shape of the roots array is extended with ones on the
  203. right, one for each dimension of `x`. Scalars have dimension 0 for
  204. this action. The result is that every column of coefficients in `r`
  205. is evaluated for every element of `x`. If False, `x` is broadcast
  206. over the columns of `r` for the evaluation. This keyword is useful
  207. when `r` is multidimensional. The default value is True.
  208. Returns
  209. -------
  210. values : ndarray, compatible object
  211. The shape of the returned array is described above.
  212. See Also
  213. --------
  214. polyroots, polyfromroots, polyval
  215. Examples
  216. --------
  217. >>> from numpy.polynomial.polynomial import polyvalfromroots
  218. >>> polyvalfromroots(1, [1,2,3])
  219. 0.0
  220. >>> a = np.arange(4).reshape(2,2)
  221. >>> a
  222. array([[0, 1],
  223. [2, 3]])
  224. >>> polyvalfromroots(a, [-1, 0, 1])
  225. array([[ -0., 0.],
  226. [ 6., 24.]])
  227. >>> r = np.arange(-2, 2).reshape(2,2) # multidimensional coefficients
  228. >>> r # each column of r defines one polynomial
  229. array([[-2, -1],
  230. [ 0, 1]])
  231. >>> b = [-2, 1]
  232. >>> polyvalfromroots(b, r, tensor=True)
  233. array([[-0., 3.],
  234. [ 3., 0.]])
  235. >>> polyvalfromroots(b, r, tensor=False)
  236. array([-0., 0.])
  237. """
  238. r = np.array(r, ndmin=1, copy=0)
  239. if r.dtype.char in '?bBhHiIlLqQpP':
  240. r = r.astype(np.double)
  241. if isinstance(x, (tuple, list)):
  242. x = np.asarray(x)
  243. if isinstance(x, np.ndarray):
  244. if tensor:
  245. r = r.reshape(r.shape + (1,)*x.ndim)
  246. elif x.ndim >= r.ndim:
  247. raise ValueError("x.ndim must be < r.ndim when tensor == "
  248. "False")
  249. return np.prod(x - r, axis=0)
  250. try:
  251. from numpy.testing import suppress_warnings
  252. except ImportError:
  253. class suppress_warnings(object):
  254. """
  255. Context manager and decorator doing much the same as
  256. ``warnings.catch_warnings``.
  257. However, it also provides a filter mechanism to work around
  258. https://bugs.python.org/issue4180.
  259. This bug causes Python before 3.4 to not reliably show warnings again
  260. after they have been ignored once (even within catch_warnings). It
  261. means that no "ignore" filter can be used easily, since following
  262. tests might need to see the warning. Additionally it allows easier
  263. specificity for testing warnings and can be nested.
  264. Parameters
  265. ----------
  266. forwarding_rule : str, optional
  267. One of "always", "once", "module", or "location". Analogous to
  268. the usual warnings module filter mode, it is useful to reduce
  269. noise mostly on the outmost level. Unsuppressed and unrecorded
  270. warnings will be forwarded based on this rule. Defaults to "always".
  271. "location" is equivalent to the warnings "default", match by exact
  272. location the warning warning originated from.
  273. Notes
  274. -----
  275. Filters added inside the context manager will be discarded again
  276. when leaving it. Upon entering all filters defined outside a
  277. context will be applied automatically.
  278. When a recording filter is added, matching warnings are stored in the
  279. ``log`` attribute as well as in the list returned by ``record``.
  280. If filters are added and the ``module`` keyword is given, the
  281. warning registry of this module will additionally be cleared when
  282. applying it, entering the context, or exiting it. This could cause
  283. warnings to appear a second time after leaving the context if they
  284. were configured to be printed once (default) and were already
  285. printed before the context was entered.
  286. Nesting this context manager will work as expected when the
  287. forwarding rule is "always" (default). Unfiltered and unrecorded
  288. warnings will be passed out and be matched by the outer level.
  289. On the outmost level they will be printed (or caught by another
  290. warnings context). The forwarding rule argument can modify this
  291. behaviour.
  292. Like ``catch_warnings`` this context manager is not threadsafe.
  293. Examples
  294. --------
  295. >>> with suppress_warnings() as sup:
  296. ... sup.filter(DeprecationWarning, "Some text")
  297. ... sup.filter(module=np.ma.core)
  298. ... log = sup.record(FutureWarning, "Does this occur?")
  299. ... command_giving_warnings()
  300. ... # The FutureWarning was given once, the filtered warnings were
  301. ... # ignored. All other warnings abide outside settings (may be
  302. ... # printed/error)
  303. ... assert_(len(log) == 1)
  304. ... assert_(len(sup.log) == 1) # also stored in log attribute
  305. Or as a decorator:
  306. >>> sup = suppress_warnings()
  307. >>> sup.filter(module=np.ma.core) # module must match exact
  308. >>> @sup
  309. >>> def some_function():
  310. ... # do something which causes a warning in np.ma.core
  311. ... pass
  312. """
  313. def __init__(self, forwarding_rule="always"):
  314. self._entered = False
  315. # Suppressions are either instance or defined inside one with block:
  316. self._suppressions = []
  317. if forwarding_rule not in {"always", "module", "once", "location"}:
  318. raise ValueError("unsupported forwarding rule.")
  319. self._forwarding_rule = forwarding_rule
  320. def _clear_registries(self):
  321. if hasattr(warnings, "_filters_mutated"):
  322. # clearing the registry should not be necessary on new pythons,
  323. # instead the filters should be mutated.
  324. warnings._filters_mutated()
  325. return
  326. # Simply clear the registry, this should normally be harmless,
  327. # note that on new pythons it would be invalidated anyway.
  328. for module in self._tmp_modules:
  329. if hasattr(module, "__warningregistry__"):
  330. module.__warningregistry__.clear()
  331. def _filter(self, category=Warning, message="", module=None, record=False):
  332. if record:
  333. record = [] # The log where to store warnings
  334. else:
  335. record = None
  336. if self._entered:
  337. if module is None:
  338. warnings.filterwarnings(
  339. "always", category=category, message=message)
  340. else:
  341. module_regex = module.__name__.replace('.', r'\.') + '$'
  342. warnings.filterwarnings(
  343. "always", category=category, message=message,
  344. module=module_regex)
  345. self._tmp_modules.add(module)
  346. self._clear_registries()
  347. self._tmp_suppressions.append(
  348. (category, message, re.compile(message, re.I), module, record))
  349. else:
  350. self._suppressions.append(
  351. (category, message, re.compile(message, re.I), module, record))
  352. return record
  353. def filter(self, category=Warning, message="", module=None):
  354. """
  355. Add a new suppressing filter or apply it if the state is entered.
  356. Parameters
  357. ----------
  358. category : class, optional
  359. Warning class to filter
  360. message : string, optional
  361. Regular expression matching the warning message.
  362. module : module, optional
  363. Module to filter for. Note that the module (and its file)
  364. must match exactly and cannot be a submodule. This may make
  365. it unreliable for external modules.
  366. Notes
  367. -----
  368. When added within a context, filters are only added inside
  369. the context and will be forgotten when the context is exited.
  370. """
  371. self._filter(category=category, message=message, module=module,
  372. record=False)
  373. def record(self, category=Warning, message="", module=None):
  374. """
  375. Append a new recording filter or apply it if the state is entered.
  376. All warnings matching will be appended to the ``log`` attribute.
  377. Parameters
  378. ----------
  379. category : class, optional
  380. Warning class to filter
  381. message : string, optional
  382. Regular expression matching the warning message.
  383. module : module, optional
  384. Module to filter for. Note that the module (and its file)
  385. must match exactly and cannot be a submodule. This may make
  386. it unreliable for external modules.
  387. Returns
  388. -------
  389. log : list
  390. A list which will be filled with all matched warnings.
  391. Notes
  392. -----
  393. When added within a context, filters are only added inside
  394. the context and will be forgotten when the context is exited.
  395. """
  396. return self._filter(category=category, message=message, module=module,
  397. record=True)
  398. def __enter__(self):
  399. if self._entered:
  400. raise RuntimeError("cannot enter suppress_warnings twice.")
  401. self._orig_show = warnings.showwarning
  402. self._filters = warnings.filters
  403. warnings.filters = self._filters[:]
  404. self._entered = True
  405. self._tmp_suppressions = []
  406. self._tmp_modules = set()
  407. self._forwarded = set()
  408. self.log = [] # reset global log (no need to keep same list)
  409. for cat, mess, _, mod, log in self._suppressions:
  410. if log is not None:
  411. del log[:] # clear the log
  412. if mod is None:
  413. warnings.filterwarnings(
  414. "always", category=cat, message=mess)
  415. else:
  416. module_regex = mod.__name__.replace('.', r'\.') + '$'
  417. warnings.filterwarnings(
  418. "always", category=cat, message=mess,
  419. module=module_regex)
  420. self._tmp_modules.add(mod)
  421. warnings.showwarning = self._showwarning
  422. self._clear_registries()
  423. return self
  424. def __exit__(self, *exc_info):
  425. warnings.showwarning = self._orig_show
  426. warnings.filters = self._filters
  427. self._clear_registries()
  428. self._entered = False
  429. del self._orig_show
  430. del self._filters
  431. def _showwarning(self, message, category, filename, lineno,
  432. *args, **kwargs):
  433. use_warnmsg = kwargs.pop("use_warnmsg", None)
  434. for cat, _, pattern, mod, rec in (
  435. self._suppressions + self._tmp_suppressions)[::-1]:
  436. if (issubclass(category, cat) and
  437. pattern.match(message.args[0]) is not None):
  438. if mod is None:
  439. # Message and category match, either recorded or ignored
  440. if rec is not None:
  441. msg = WarningMessage(message, category, filename,
  442. lineno, **kwargs)
  443. self.log.append(msg)
  444. rec.append(msg)
  445. return
  446. # Use startswith, because warnings strips the c or o from
  447. # .pyc/.pyo files.
  448. elif mod.__file__.startswith(filename):
  449. # The message and module (filename) match
  450. if rec is not None:
  451. msg = WarningMessage(message, category, filename,
  452. lineno, **kwargs)
  453. self.log.append(msg)
  454. rec.append(msg)
  455. return
  456. # There is no filter in place, so pass to the outside handler
  457. # unless we should only pass it once
  458. if self._forwarding_rule == "always":
  459. if use_warnmsg is None:
  460. self._orig_show(message, category, filename, lineno,
  461. *args, **kwargs)
  462. else:
  463. self._orig_showmsg(use_warnmsg)
  464. return
  465. if self._forwarding_rule == "once":
  466. signature = (message.args, category)
  467. elif self._forwarding_rule == "module":
  468. signature = (message.args, category, filename)
  469. elif self._forwarding_rule == "location":
  470. signature = (message.args, category, filename, lineno)
  471. if signature in self._forwarded:
  472. return
  473. self._forwarded.add(signature)
  474. if use_warnmsg is None:
  475. self._orig_show(message, category, filename, lineno, *args,
  476. **kwargs)
  477. else:
  478. self._orig_showmsg(use_warnmsg)
  479. def __call__(self, func):
  480. """
  481. Function decorator to apply certain suppressions to a whole
  482. function.
  483. """
  484. @wraps(func)
  485. def new_func(*args, **kwargs):
  486. with self:
  487. return func(*args, **kwargs)
  488. return new_func
  489. if NumpyVersion(np.__version__) >= '1.10.0':
  490. from numpy import cov
  491. else:
  492. from numpy import array, average, dot
  493. def cov(m, y=None, rowvar=True, bias=False, ddof=None, fweights=None,
  494. aweights=None):
  495. """
  496. Estimate a covariance matrix, given data and weights.
  497. Covariance indicates the level to which two variables vary together.
  498. If we examine N-dimensional samples, :math:`X = [x_1, x_2, ... x_N]^T`,
  499. then the covariance matrix element :math:`C_{ij}` is the covariance of
  500. :math:`x_i` and :math:`x_j`. The element :math:`C_{ii}` is the variance
  501. of :math:`x_i`.
  502. See the notes for an outline of the algorithm.
  503. Parameters
  504. ----------
  505. m : array_like
  506. A 1-D or 2-D array containing multiple variables and observations.
  507. Each row of `m` represents a variable, and each column a single
  508. observation of all those variables. Also see `rowvar` below.
  509. y : array_like, optional
  510. An additional set of variables and observations. `y` has the same form
  511. as that of `m`.
  512. rowvar : bool, optional
  513. If `rowvar` is True (default), then each row represents a
  514. variable, with observations in the columns. Otherwise, the relationship
  515. is transposed: each column represents a variable, while the rows
  516. contain observations.
  517. bias : bool, optional
  518. Default normalization (False) is by ``(N - 1)``, where ``N`` is the
  519. number of observations given (unbiased estimate). If `bias` is True,
  520. then normalization is by ``N``. These values can be overridden by using
  521. the keyword ``ddof`` in numpy versions >= 1.5.
  522. ddof : int, optional
  523. If not ``None`` the default value implied by `bias` is overridden.
  524. Note that ``ddof=1`` will return the unbiased estimate, even if both
  525. `fweights` and `aweights` are specified, and ``ddof=0`` will return
  526. the simple average. See the notes for the details. The default value
  527. is ``None``.
  528. .. versionadded:: 1.5
  529. fweights : array_like, int, optional
  530. 1-D array of integer freguency weights; the number of times each
  531. observation vector should be repeated.
  532. .. versionadded:: 1.10
  533. aweights : array_like, optional
  534. 1-D array of observation vector weights. These relative weights are
  535. typically large for observations considered "important" and smaller for
  536. observations considered less "important". If ``ddof=0`` the array of
  537. weights can be used to assign probabilities to observation vectors.
  538. .. versionadded:: 1.10
  539. Returns
  540. -------
  541. out : ndarray
  542. The covariance matrix of the variables.
  543. See Also
  544. --------
  545. corrcoef : Normalized covariance matrix
  546. Notes
  547. -----
  548. Assume that the observations are in the columns of the observation
  549. array `m` and let ``f = fweights`` and ``a = aweights`` for brevity. The
  550. steps to compute the weighted covariance are as follows::
  551. >>> w = f * a
  552. >>> v1 = np.sum(w)
  553. >>> v2 = np.sum(w * a)
  554. >>> m -= np.sum(m * w, axis=1, keepdims=True) / v1
  555. >>> cov = np.dot(m * w, m.T) * v1 / (v1**2 - ddof * v2)
  556. Note that when ``a == 1``, the normalization factor
  557. ``v1 / (v1**2 - ddof * v2)`` goes over to ``1 / (np.sum(f) - ddof)``
  558. as it should.
  559. Examples
  560. --------
  561. Consider two variables, :math:`x_0` and :math:`x_1`, which
  562. correlate perfectly, but in opposite directions:
  563. >>> x = np.array([[0, 2], [1, 1], [2, 0]]).T
  564. >>> x
  565. array([[0, 1, 2],
  566. [2, 1, 0]])
  567. Note how :math:`x_0` increases while :math:`x_1` decreases. The covariance
  568. matrix shows this clearly:
  569. >>> np.cov(x)
  570. array([[ 1., -1.],
  571. [-1., 1.]])
  572. Note that element :math:`C_{0,1}`, which shows the correlation between
  573. :math:`x_0` and :math:`x_1`, is negative.
  574. Further, note how `x` and `y` are combined:
  575. >>> x = [-2.1, -1, 4.3]
  576. >>> y = [3, 1.1, 0.12]
  577. >>> X = np.stack((x, y), axis=0)
  578. >>> print(np.cov(X))
  579. [[ 11.71 -4.286 ]
  580. [ -4.286 2.14413333]]
  581. >>> print(np.cov(x, y))
  582. [[ 11.71 -4.286 ]
  583. [ -4.286 2.14413333]]
  584. >>> print(np.cov(x))
  585. 11.71
  586. """
  587. # Check inputs
  588. if ddof is not None and ddof != int(ddof):
  589. raise ValueError(
  590. "ddof must be integer")
  591. # Handles complex arrays too
  592. m = np.asarray(m)
  593. if m.ndim > 2:
  594. raise ValueError("m has more than 2 dimensions")
  595. if y is None:
  596. dtype = np.result_type(m, np.float64)
  597. else:
  598. y = np.asarray(y)
  599. if y.ndim > 2:
  600. raise ValueError("y has more than 2 dimensions")
  601. dtype = np.result_type(m, y, np.float64)
  602. X = array(m, ndmin=2, dtype=dtype)
  603. if not rowvar and X.shape[0] != 1:
  604. X = X.T
  605. if X.shape[0] == 0:
  606. return np.array([]).reshape(0, 0)
  607. if y is not None:
  608. y = array(y, copy=False, ndmin=2, dtype=dtype)
  609. if not rowvar and y.shape[0] != 1:
  610. y = y.T
  611. X = np.concatenate((X, y), axis=0)
  612. if ddof is None:
  613. if bias == 0:
  614. ddof = 1
  615. else:
  616. ddof = 0
  617. # Get the product of frequencies and weights
  618. w = None
  619. if fweights is not None:
  620. fweights = np.asarray(fweights, dtype=float)
  621. if not np.all(fweights == np.around(fweights)):
  622. raise TypeError(
  623. "fweights must be integer")
  624. if fweights.ndim > 1:
  625. raise RuntimeError(
  626. "cannot handle multidimensional fweights")
  627. if fweights.shape[0] != X.shape[1]:
  628. raise RuntimeError(
  629. "incompatible numbers of samples and fweights")
  630. if any(fweights < 0):
  631. raise ValueError(
  632. "fweights cannot be negative")
  633. w = fweights
  634. if aweights is not None:
  635. aweights = np.asarray(aweights, dtype=float)
  636. if aweights.ndim > 1:
  637. raise RuntimeError(
  638. "cannot handle multidimensional aweights")
  639. if aweights.shape[0] != X.shape[1]:
  640. raise RuntimeError(
  641. "incompatible numbers of samples and aweights")
  642. if any(aweights < 0):
  643. raise ValueError(
  644. "aweights cannot be negative")
  645. if w is None:
  646. w = aweights
  647. else:
  648. w *= aweights
  649. avg, w_sum = average(X, axis=1, weights=w, returned=True)
  650. w_sum = w_sum[0]
  651. # Determine the normalization
  652. if w is None:
  653. fact = X.shape[1] - ddof
  654. elif ddof == 0:
  655. fact = w_sum
  656. elif aweights is None:
  657. fact = w_sum - ddof
  658. else:
  659. fact = w_sum - ddof*sum(w*aweights)/w_sum
  660. if fact <= 0:
  661. warnings.warn("Degrees of freedom <= 0 for slice",
  662. RuntimeWarning, stacklevel=2)
  663. fact = 0.0
  664. X -= avg[:, None]
  665. if w is None:
  666. X_T = X.T
  667. else:
  668. X_T = (X*w).T
  669. c = dot(X, X_T.conj())
  670. c *= 1. / np.float64(fact)
  671. return c.squeeze()