polyutils.py 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413
  1. """
  2. Utility classes and functions for the polynomial modules.
  3. This module provides: error and warning objects; a polynomial base class;
  4. and some routines used in both the `polynomial` and `chebyshev` modules.
  5. Error objects
  6. -------------
  7. .. autosummary::
  8. :toctree: generated/
  9. PolyError base class for this sub-package's errors.
  10. PolyDomainError raised when domains are mismatched.
  11. Warning objects
  12. ---------------
  13. .. autosummary::
  14. :toctree: generated/
  15. RankWarning raised in least-squares fit for rank-deficient matrix.
  16. Base class
  17. ----------
  18. .. autosummary::
  19. :toctree: generated/
  20. PolyBase Obsolete base class for the polynomial classes. Do not use.
  21. Functions
  22. ---------
  23. .. autosummary::
  24. :toctree: generated/
  25. as_series convert list of array_likes into 1-D arrays of common type.
  26. trimseq remove trailing zeros.
  27. trimcoef remove small trailing coefficients.
  28. getdomain return the domain appropriate for a given set of abscissae.
  29. mapdomain maps points between domains.
  30. mapparms parameters of the linear map between domains.
  31. """
  32. from __future__ import division, absolute_import, print_function
  33. import numpy as np
  34. __all__ = [
  35. 'RankWarning', 'PolyError', 'PolyDomainError', 'as_series', 'trimseq',
  36. 'trimcoef', 'getdomain', 'mapdomain', 'mapparms', 'PolyBase']
  37. #
  38. # Warnings and Exceptions
  39. #
  40. class RankWarning(UserWarning):
  41. """Issued by chebfit when the design matrix is rank deficient."""
  42. pass
  43. class PolyError(Exception):
  44. """Base class for errors in this module."""
  45. pass
  46. class PolyDomainError(PolyError):
  47. """Issued by the generic Poly class when two domains don't match.
  48. This is raised when an binary operation is passed Poly objects with
  49. different domains.
  50. """
  51. pass
  52. #
  53. # Base class for all polynomial types
  54. #
  55. class PolyBase(object):
  56. """
  57. Base class for all polynomial types.
  58. Deprecated in numpy 1.9.0, use the abstract
  59. ABCPolyBase class instead. Note that the latter
  60. requires a number of virtual functions to be
  61. implemented.
  62. """
  63. pass
  64. #
  65. # Helper functions to convert inputs to 1-D arrays
  66. #
  67. def trimseq(seq):
  68. """Remove small Poly series coefficients.
  69. Parameters
  70. ----------
  71. seq : sequence
  72. Sequence of Poly series coefficients. This routine fails for
  73. empty sequences.
  74. Returns
  75. -------
  76. series : sequence
  77. Subsequence with trailing zeros removed. If the resulting sequence
  78. would be empty, return the first element. The returned sequence may
  79. or may not be a view.
  80. Notes
  81. -----
  82. Do not lose the type info if the sequence contains unknown objects.
  83. """
  84. if len(seq) == 0:
  85. return seq
  86. else:
  87. for i in range(len(seq) - 1, -1, -1):
  88. if seq[i] != 0:
  89. break
  90. return seq[:i+1]
  91. def as_series(alist, trim=True):
  92. """
  93. Return argument as a list of 1-d arrays.
  94. The returned list contains array(s) of dtype double, complex double, or
  95. object. A 1-d argument of shape ``(N,)`` is parsed into ``N`` arrays of
  96. size one; a 2-d argument of shape ``(M,N)`` is parsed into ``M`` arrays
  97. of size ``N`` (i.e., is "parsed by row"); and a higher dimensional array
  98. raises a Value Error if it is not first reshaped into either a 1-d or 2-d
  99. array.
  100. Parameters
  101. ----------
  102. alist : array_like
  103. A 1- or 2-d array_like
  104. trim : boolean, optional
  105. When True, trailing zeros are removed from the inputs.
  106. When False, the inputs are passed through intact.
  107. Returns
  108. -------
  109. [a1, a2,...] : list of 1-D arrays
  110. A copy of the input data as a list of 1-d arrays.
  111. Raises
  112. ------
  113. ValueError
  114. Raised when `as_series` cannot convert its input to 1-d arrays, or at
  115. least one of the resulting arrays is empty.
  116. Examples
  117. --------
  118. >>> from numpy.polynomial import polyutils as pu
  119. >>> a = np.arange(4)
  120. >>> pu.as_series(a)
  121. [array([ 0.]), array([ 1.]), array([ 2.]), array([ 3.])]
  122. >>> b = np.arange(6).reshape((2,3))
  123. >>> pu.as_series(b)
  124. [array([ 0., 1., 2.]), array([ 3., 4., 5.])]
  125. >>> pu.as_series((1, np.arange(3), np.arange(2, dtype=np.float16)))
  126. [array([ 1.]), array([ 0., 1., 2.]), array([ 0., 1.])]
  127. >>> pu.as_series([2, [1.1, 0.]])
  128. [array([ 2.]), array([ 1.1])]
  129. >>> pu.as_series([2, [1.1, 0.]], trim=False)
  130. [array([ 2.]), array([ 1.1, 0. ])]
  131. """
  132. arrays = [np.array(a, ndmin=1, copy=0) for a in alist]
  133. if min([a.size for a in arrays]) == 0:
  134. raise ValueError("Coefficient array is empty")
  135. if any([a.ndim != 1 for a in arrays]):
  136. raise ValueError("Coefficient array is not 1-d")
  137. if trim:
  138. arrays = [trimseq(a) for a in arrays]
  139. if any([a.dtype == np.dtype(object) for a in arrays]):
  140. ret = []
  141. for a in arrays:
  142. if a.dtype != np.dtype(object):
  143. tmp = np.empty(len(a), dtype=np.dtype(object))
  144. tmp[:] = a[:]
  145. ret.append(tmp)
  146. else:
  147. ret.append(a.copy())
  148. else:
  149. try:
  150. dtype = np.common_type(*arrays)
  151. except Exception:
  152. raise ValueError("Coefficient arrays have no common type")
  153. ret = [np.array(a, copy=1, dtype=dtype) for a in arrays]
  154. return ret
  155. def trimcoef(c, tol=0):
  156. """
  157. Remove "small" "trailing" coefficients from a polynomial.
  158. "Small" means "small in absolute value" and is controlled by the
  159. parameter `tol`; "trailing" means highest order coefficient(s), e.g., in
  160. ``[0, 1, 1, 0, 0]`` (which represents ``0 + x + x**2 + 0*x**3 + 0*x**4``)
  161. both the 3-rd and 4-th order coefficients would be "trimmed."
  162. Parameters
  163. ----------
  164. c : array_like
  165. 1-d array of coefficients, ordered from lowest order to highest.
  166. tol : number, optional
  167. Trailing (i.e., highest order) elements with absolute value less
  168. than or equal to `tol` (default value is zero) are removed.
  169. Returns
  170. -------
  171. trimmed : ndarray
  172. 1-d array with trailing zeros removed. If the resulting series
  173. would be empty, a series containing a single zero is returned.
  174. Raises
  175. ------
  176. ValueError
  177. If `tol` < 0
  178. See Also
  179. --------
  180. trimseq
  181. Examples
  182. --------
  183. >>> from numpy.polynomial import polyutils as pu
  184. >>> pu.trimcoef((0,0,3,0,5,0,0))
  185. array([ 0., 0., 3., 0., 5.])
  186. >>> pu.trimcoef((0,0,1e-3,0,1e-5,0,0),1e-3) # item == tol is trimmed
  187. array([ 0.])
  188. >>> i = complex(0,1) # works for complex
  189. >>> pu.trimcoef((3e-4,1e-3*(1-i),5e-4,2e-5*(1+i)), 1e-3)
  190. array([ 0.0003+0.j , 0.0010-0.001j])
  191. """
  192. if tol < 0:
  193. raise ValueError("tol must be non-negative")
  194. [c] = as_series([c])
  195. [ind] = np.nonzero(np.abs(c) > tol)
  196. if len(ind) == 0:
  197. return c[:1]*0
  198. else:
  199. return c[:ind[-1] + 1].copy()
  200. def getdomain(x):
  201. """
  202. Return a domain suitable for given abscissae.
  203. Find a domain suitable for a polynomial or Chebyshev series
  204. defined at the values supplied.
  205. Parameters
  206. ----------
  207. x : array_like
  208. 1-d array of abscissae whose domain will be determined.
  209. Returns
  210. -------
  211. domain : ndarray
  212. 1-d array containing two values. If the inputs are complex, then
  213. the two returned points are the lower left and upper right corners
  214. of the smallest rectangle (aligned with the axes) in the complex
  215. plane containing the points `x`. If the inputs are real, then the
  216. two points are the ends of the smallest interval containing the
  217. points `x`.
  218. See Also
  219. --------
  220. mapparms, mapdomain
  221. Examples
  222. --------
  223. >>> from numpy.polynomial import polyutils as pu
  224. >>> points = np.arange(4)**2 - 5; points
  225. array([-5, -4, -1, 4])
  226. >>> pu.getdomain(points)
  227. array([-5., 4.])
  228. >>> c = np.exp(complex(0,1)*np.pi*np.arange(12)/6) # unit circle
  229. >>> pu.getdomain(c)
  230. array([-1.-1.j, 1.+1.j])
  231. """
  232. [x] = as_series([x], trim=False)
  233. if x.dtype.char in np.typecodes['Complex']:
  234. rmin, rmax = x.real.min(), x.real.max()
  235. imin, imax = x.imag.min(), x.imag.max()
  236. return np.array((complex(rmin, imin), complex(rmax, imax)))
  237. else:
  238. return np.array((x.min(), x.max()))
  239. def mapparms(old, new):
  240. """
  241. Linear map parameters between domains.
  242. Return the parameters of the linear map ``offset + scale*x`` that maps
  243. `old` to `new` such that ``old[i] -> new[i]``, ``i = 0, 1``.
  244. Parameters
  245. ----------
  246. old, new : array_like
  247. Domains. Each domain must (successfully) convert to a 1-d array
  248. containing precisely two values.
  249. Returns
  250. -------
  251. offset, scale : scalars
  252. The map ``L(x) = offset + scale*x`` maps the first domain to the
  253. second.
  254. See Also
  255. --------
  256. getdomain, mapdomain
  257. Notes
  258. -----
  259. Also works for complex numbers, and thus can be used to calculate the
  260. parameters required to map any line in the complex plane to any other
  261. line therein.
  262. Examples
  263. --------
  264. >>> from numpy.polynomial import polyutils as pu
  265. >>> pu.mapparms((-1,1),(-1,1))
  266. (0.0, 1.0)
  267. >>> pu.mapparms((1,-1),(-1,1))
  268. (0.0, -1.0)
  269. >>> i = complex(0,1)
  270. >>> pu.mapparms((-i,-1),(1,i))
  271. ((1+1j), (1+0j))
  272. """
  273. oldlen = old[1] - old[0]
  274. newlen = new[1] - new[0]
  275. off = (old[1]*new[0] - old[0]*new[1])/oldlen
  276. scl = newlen/oldlen
  277. return off, scl
  278. def mapdomain(x, old, new):
  279. """
  280. Apply linear map to input points.
  281. The linear map ``offset + scale*x`` that maps the domain `old` to
  282. the domain `new` is applied to the points `x`.
  283. Parameters
  284. ----------
  285. x : array_like
  286. Points to be mapped. If `x` is a subtype of ndarray the subtype
  287. will be preserved.
  288. old, new : array_like
  289. The two domains that determine the map. Each must (successfully)
  290. convert to 1-d arrays containing precisely two values.
  291. Returns
  292. -------
  293. x_out : ndarray
  294. Array of points of the same shape as `x`, after application of the
  295. linear map between the two domains.
  296. See Also
  297. --------
  298. getdomain, mapparms
  299. Notes
  300. -----
  301. Effectively, this implements:
  302. .. math ::
  303. x\\_out = new[0] + m(x - old[0])
  304. where
  305. .. math ::
  306. m = \\frac{new[1]-new[0]}{old[1]-old[0]}
  307. Examples
  308. --------
  309. >>> from numpy.polynomial import polyutils as pu
  310. >>> old_domain = (-1,1)
  311. >>> new_domain = (0,2*np.pi)
  312. >>> x = np.linspace(-1,1,6); x
  313. array([-1. , -0.6, -0.2, 0.2, 0.6, 1. ])
  314. >>> x_out = pu.mapdomain(x, old_domain, new_domain); x_out
  315. array([ 0. , 1.25663706, 2.51327412, 3.76991118, 5.02654825,
  316. 6.28318531])
  317. >>> x - pu.mapdomain(x_out, new_domain, old_domain)
  318. array([ 0., 0., 0., 0., 0., 0.])
  319. Also works for complex numbers (and thus can be used to map any line in
  320. the complex plane to any other line therein).
  321. >>> i = complex(0,1)
  322. >>> old = (-1 - i, 1 + i)
  323. >>> new = (-1 + i, 1 - i)
  324. >>> z = np.linspace(old[0], old[1], 6); z
  325. array([-1.0-1.j , -0.6-0.6j, -0.2-0.2j, 0.2+0.2j, 0.6+0.6j, 1.0+1.j ])
  326. >>> new_z = P.mapdomain(z, old, new); new_z
  327. array([-1.0+1.j , -0.6+0.6j, -0.2+0.2j, 0.2-0.2j, 0.6-0.6j, 1.0-1.j ])
  328. """
  329. x = np.asanyarray(x)
  330. off, scl = mapparms(old, new)
  331. return off + scl*x