_methods.py 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289
  1. """
  2. Array methods which are called by both the C-code for the method
  3. and the Python code for the NumPy-namespace function
  4. """
  5. import warnings
  6. from numpy.core import multiarray as mu
  7. from numpy.core import umath as um
  8. from numpy.core._asarray import asanyarray
  9. from numpy.core import numerictypes as nt
  10. from numpy.core import _exceptions
  11. from numpy._globals import _NoValue
  12. from numpy.compat import pickle, os_fspath, contextlib_nullcontext
  13. # save those O(100) nanoseconds!
  14. umr_maximum = um.maximum.reduce
  15. umr_minimum = um.minimum.reduce
  16. umr_sum = um.add.reduce
  17. umr_prod = um.multiply.reduce
  18. umr_any = um.logical_or.reduce
  19. umr_all = um.logical_and.reduce
  20. # Complex types to -> (2,)float view for fast-path computation in _var()
  21. _complex_to_float = {
  22. nt.dtype(nt.csingle) : nt.dtype(nt.single),
  23. nt.dtype(nt.cdouble) : nt.dtype(nt.double),
  24. }
  25. # Special case for windows: ensure double takes precedence
  26. if nt.dtype(nt.longdouble) != nt.dtype(nt.double):
  27. _complex_to_float.update({
  28. nt.dtype(nt.clongdouble) : nt.dtype(nt.longdouble),
  29. })
  30. # avoid keyword arguments to speed up parsing, saves about 15%-20% for very
  31. # small reductions
  32. def _amax(a, axis=None, out=None, keepdims=False,
  33. initial=_NoValue, where=True):
  34. return umr_maximum(a, axis, None, out, keepdims, initial, where)
  35. def _amin(a, axis=None, out=None, keepdims=False,
  36. initial=_NoValue, where=True):
  37. return umr_minimum(a, axis, None, out, keepdims, initial, where)
  38. def _sum(a, axis=None, dtype=None, out=None, keepdims=False,
  39. initial=_NoValue, where=True):
  40. return umr_sum(a, axis, dtype, out, keepdims, initial, where)
  41. def _prod(a, axis=None, dtype=None, out=None, keepdims=False,
  42. initial=_NoValue, where=True):
  43. return umr_prod(a, axis, dtype, out, keepdims, initial, where)
  44. def _any(a, axis=None, dtype=None, out=None, keepdims=False, *, where=True):
  45. # Parsing keyword arguments is currently fairly slow, so avoid it for now
  46. if where is True:
  47. return umr_any(a, axis, dtype, out, keepdims)
  48. return umr_any(a, axis, dtype, out, keepdims, where=where)
  49. def _all(a, axis=None, dtype=None, out=None, keepdims=False, *, where=True):
  50. # Parsing keyword arguments is currently fairly slow, so avoid it for now
  51. if where is True:
  52. return umr_all(a, axis, dtype, out, keepdims)
  53. return umr_all(a, axis, dtype, out, keepdims, where=where)
  54. def _count_reduce_items(arr, axis, keepdims=False, where=True):
  55. # fast-path for the default case
  56. if where is True:
  57. # no boolean mask given, calculate items according to axis
  58. if axis is None:
  59. axis = tuple(range(arr.ndim))
  60. elif not isinstance(axis, tuple):
  61. axis = (axis,)
  62. items = nt.intp(1)
  63. for ax in axis:
  64. items *= arr.shape[mu.normalize_axis_index(ax, arr.ndim)]
  65. else:
  66. # TODO: Optimize case when `where` is broadcast along a non-reduction
  67. # axis and full sum is more excessive than needed.
  68. # guarded to protect circular imports
  69. from numpy.lib.stride_tricks import broadcast_to
  70. # count True values in (potentially broadcasted) boolean mask
  71. items = umr_sum(broadcast_to(where, arr.shape), axis, nt.intp, None,
  72. keepdims)
  73. return items
  74. # Numpy 1.17.0, 2019-02-24
  75. # Various clip behavior deprecations, marked with _clip_dep as a prefix.
  76. def _clip_dep_is_scalar_nan(a):
  77. # guarded to protect circular imports
  78. from numpy.core.fromnumeric import ndim
  79. if ndim(a) != 0:
  80. return False
  81. try:
  82. return um.isnan(a)
  83. except TypeError:
  84. return False
  85. def _clip_dep_is_byte_swapped(a):
  86. if isinstance(a, mu.ndarray):
  87. return not a.dtype.isnative
  88. return False
  89. def _clip_dep_invoke_with_casting(ufunc, *args, out=None, casting=None, **kwargs):
  90. # normal path
  91. if casting is not None:
  92. return ufunc(*args, out=out, casting=casting, **kwargs)
  93. # try to deal with broken casting rules
  94. try:
  95. return ufunc(*args, out=out, **kwargs)
  96. except _exceptions._UFuncOutputCastingError as e:
  97. # Numpy 1.17.0, 2019-02-24
  98. warnings.warn(
  99. "Converting the output of clip from {!r} to {!r} is deprecated. "
  100. "Pass `casting=\"unsafe\"` explicitly to silence this warning, or "
  101. "correct the type of the variables.".format(e.from_, e.to),
  102. DeprecationWarning,
  103. stacklevel=2
  104. )
  105. return ufunc(*args, out=out, casting="unsafe", **kwargs)
  106. def _clip(a, min=None, max=None, out=None, *, casting=None, **kwargs):
  107. if min is None and max is None:
  108. raise ValueError("One of max or min must be given")
  109. # Numpy 1.17.0, 2019-02-24
  110. # This deprecation probably incurs a substantial slowdown for small arrays,
  111. # it will be good to get rid of it.
  112. if not _clip_dep_is_byte_swapped(a) and not _clip_dep_is_byte_swapped(out):
  113. using_deprecated_nan = False
  114. if _clip_dep_is_scalar_nan(min):
  115. min = -float('inf')
  116. using_deprecated_nan = True
  117. if _clip_dep_is_scalar_nan(max):
  118. max = float('inf')
  119. using_deprecated_nan = True
  120. if using_deprecated_nan:
  121. warnings.warn(
  122. "Passing `np.nan` to mean no clipping in np.clip has always "
  123. "been unreliable, and is now deprecated. "
  124. "In future, this will always return nan, like it already does "
  125. "when min or max are arrays that contain nan. "
  126. "To skip a bound, pass either None or an np.inf of an "
  127. "appropriate sign.",
  128. DeprecationWarning,
  129. stacklevel=2
  130. )
  131. if min is None:
  132. return _clip_dep_invoke_with_casting(
  133. um.minimum, a, max, out=out, casting=casting, **kwargs)
  134. elif max is None:
  135. return _clip_dep_invoke_with_casting(
  136. um.maximum, a, min, out=out, casting=casting, **kwargs)
  137. else:
  138. return _clip_dep_invoke_with_casting(
  139. um.clip, a, min, max, out=out, casting=casting, **kwargs)
  140. def _mean(a, axis=None, dtype=None, out=None, keepdims=False, *, where=True):
  141. arr = asanyarray(a)
  142. is_float16_result = False
  143. rcount = _count_reduce_items(arr, axis, keepdims=keepdims, where=where)
  144. if rcount == 0 if where is True else umr_any(rcount == 0, axis=None):
  145. warnings.warn("Mean of empty slice.", RuntimeWarning, stacklevel=2)
  146. # Cast bool, unsigned int, and int to float64 by default
  147. if dtype is None:
  148. if issubclass(arr.dtype.type, (nt.integer, nt.bool_)):
  149. dtype = mu.dtype('f8')
  150. elif issubclass(arr.dtype.type, nt.float16):
  151. dtype = mu.dtype('f4')
  152. is_float16_result = True
  153. ret = umr_sum(arr, axis, dtype, out, keepdims, where=where)
  154. if isinstance(ret, mu.ndarray):
  155. ret = um.true_divide(
  156. ret, rcount, out=ret, casting='unsafe', subok=False)
  157. if is_float16_result and out is None:
  158. ret = arr.dtype.type(ret)
  159. elif hasattr(ret, 'dtype'):
  160. if is_float16_result:
  161. ret = arr.dtype.type(ret / rcount)
  162. else:
  163. ret = ret.dtype.type(ret / rcount)
  164. else:
  165. ret = ret / rcount
  166. return ret
  167. def _var(a, axis=None, dtype=None, out=None, ddof=0, keepdims=False, *,
  168. where=True):
  169. arr = asanyarray(a)
  170. rcount = _count_reduce_items(arr, axis, keepdims=keepdims, where=where)
  171. # Make this warning show up on top.
  172. if ddof >= rcount if where is True else umr_any(ddof >= rcount, axis=None):
  173. warnings.warn("Degrees of freedom <= 0 for slice", RuntimeWarning,
  174. stacklevel=2)
  175. # Cast bool, unsigned int, and int to float64 by default
  176. if dtype is None and issubclass(arr.dtype.type, (nt.integer, nt.bool_)):
  177. dtype = mu.dtype('f8')
  178. # Compute the mean.
  179. # Note that if dtype is not of inexact type then arraymean will
  180. # not be either.
  181. arrmean = umr_sum(arr, axis, dtype, keepdims=True, where=where)
  182. # The shape of rcount has to match arrmean to not change the shape of out
  183. # in broadcasting. Otherwise, it cannot be stored back to arrmean.
  184. if rcount.ndim == 0:
  185. # fast-path for default case when where is True
  186. div = rcount
  187. else:
  188. # matching rcount to arrmean when where is specified as array
  189. div = rcount.reshape(arrmean.shape)
  190. if isinstance(arrmean, mu.ndarray):
  191. arrmean = um.true_divide(arrmean, div, out=arrmean, casting='unsafe',
  192. subok=False)
  193. else:
  194. arrmean = arrmean.dtype.type(arrmean / rcount)
  195. # Compute sum of squared deviations from mean
  196. # Note that x may not be inexact and that we need it to be an array,
  197. # not a scalar.
  198. x = asanyarray(arr - arrmean)
  199. if issubclass(arr.dtype.type, (nt.floating, nt.integer)):
  200. x = um.multiply(x, x, out=x)
  201. # Fast-paths for built-in complex types
  202. elif x.dtype in _complex_to_float:
  203. xv = x.view(dtype=(_complex_to_float[x.dtype], (2,)))
  204. um.multiply(xv, xv, out=xv)
  205. x = um.add(xv[..., 0], xv[..., 1], out=x.real).real
  206. # Most general case; includes handling object arrays containing imaginary
  207. # numbers and complex types with non-native byteorder
  208. else:
  209. x = um.multiply(x, um.conjugate(x), out=x).real
  210. ret = umr_sum(x, axis, dtype, out, keepdims=keepdims, where=where)
  211. # Compute degrees of freedom and make sure it is not negative.
  212. rcount = um.maximum(rcount - ddof, 0)
  213. # divide by degrees of freedom
  214. if isinstance(ret, mu.ndarray):
  215. ret = um.true_divide(
  216. ret, rcount, out=ret, casting='unsafe', subok=False)
  217. elif hasattr(ret, 'dtype'):
  218. ret = ret.dtype.type(ret / rcount)
  219. else:
  220. ret = ret / rcount
  221. return ret
  222. def _std(a, axis=None, dtype=None, out=None, ddof=0, keepdims=False, *,
  223. where=True):
  224. ret = _var(a, axis=axis, dtype=dtype, out=out, ddof=ddof,
  225. keepdims=keepdims, where=where)
  226. if isinstance(ret, mu.ndarray):
  227. ret = um.sqrt(ret, out=ret)
  228. elif hasattr(ret, 'dtype'):
  229. ret = ret.dtype.type(um.sqrt(ret))
  230. else:
  231. ret = um.sqrt(ret)
  232. return ret
  233. def _ptp(a, axis=None, out=None, keepdims=False):
  234. return um.subtract(
  235. umr_maximum(a, axis, None, out, keepdims),
  236. umr_minimum(a, axis, None, None, keepdims),
  237. out
  238. )
  239. def _dump(self, file, protocol=2):
  240. if hasattr(file, 'write'):
  241. ctx = contextlib_nullcontext(file)
  242. else:
  243. ctx = open(os_fspath(file), "wb")
  244. with ctx as f:
  245. pickle.dump(self, f, protocol=protocol)
  246. def _dumps(self, protocol=2):
  247. return pickle.dumps(self, protocol=protocol)