polyutils.py 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403
  1. """
  2. Utililty 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. reguires 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 import polynomial as P
  119. >>> a = np.arange(4)
  120. >>> P.as_series(a)
  121. [array([ 0.]), array([ 1.]), array([ 2.]), array([ 3.])]
  122. >>> b = np.arange(6).reshape((2,3))
  123. >>> P.as_series(b)
  124. [array([ 0., 1., 2.]), array([ 3., 4., 5.])]
  125. """
  126. arrays = [np.array(a, ndmin=1, copy=0) for a in alist]
  127. if min([a.size for a in arrays]) == 0:
  128. raise ValueError("Coefficient array is empty")
  129. if any([a.ndim != 1 for a in arrays]):
  130. raise ValueError("Coefficient array is not 1-d")
  131. if trim:
  132. arrays = [trimseq(a) for a in arrays]
  133. if any([a.dtype == np.dtype(object) for a in arrays]):
  134. ret = []
  135. for a in arrays:
  136. if a.dtype != np.dtype(object):
  137. tmp = np.empty(len(a), dtype=np.dtype(object))
  138. tmp[:] = a[:]
  139. ret.append(tmp)
  140. else:
  141. ret.append(a.copy())
  142. else:
  143. try:
  144. dtype = np.common_type(*arrays)
  145. except:
  146. raise ValueError("Coefficient arrays have no common type")
  147. ret = [np.array(a, copy=1, dtype=dtype) for a in arrays]
  148. return ret
  149. def trimcoef(c, tol=0):
  150. """
  151. Remove "small" "trailing" coefficients from a polynomial.
  152. "Small" means "small in absolute value" and is controlled by the
  153. parameter `tol`; "trailing" means highest order coefficient(s), e.g., in
  154. ``[0, 1, 1, 0, 0]`` (which represents ``0 + x + x**2 + 0*x**3 + 0*x**4``)
  155. both the 3-rd and 4-th order coefficients would be "trimmed."
  156. Parameters
  157. ----------
  158. c : array_like
  159. 1-d array of coefficients, ordered from lowest order to highest.
  160. tol : number, optional
  161. Trailing (i.e., highest order) elements with absolute value less
  162. than or equal to `tol` (default value is zero) are removed.
  163. Returns
  164. -------
  165. trimmed : ndarray
  166. 1-d array with trailing zeros removed. If the resulting series
  167. would be empty, a series containing a single zero is returned.
  168. Raises
  169. ------
  170. ValueError
  171. If `tol` < 0
  172. See Also
  173. --------
  174. trimseq
  175. Examples
  176. --------
  177. >>> from numpy import polynomial as P
  178. >>> P.trimcoef((0,0,3,0,5,0,0))
  179. array([ 0., 0., 3., 0., 5.])
  180. >>> P.trimcoef((0,0,1e-3,0,1e-5,0,0),1e-3) # item == tol is trimmed
  181. array([ 0.])
  182. >>> i = complex(0,1) # works for complex
  183. >>> P.trimcoef((3e-4,1e-3*(1-i),5e-4,2e-5*(1+i)), 1e-3)
  184. array([ 0.0003+0.j , 0.0010-0.001j])
  185. """
  186. if tol < 0:
  187. raise ValueError("tol must be non-negative")
  188. [c] = as_series([c])
  189. [ind] = np.where(np.abs(c) > tol)
  190. if len(ind) == 0:
  191. return c[:1]*0
  192. else:
  193. return c[:ind[-1] + 1].copy()
  194. def getdomain(x):
  195. """
  196. Return a domain suitable for given abscissae.
  197. Find a domain suitable for a polynomial or Chebyshev series
  198. defined at the values supplied.
  199. Parameters
  200. ----------
  201. x : array_like
  202. 1-d array of abscissae whose domain will be determined.
  203. Returns
  204. -------
  205. domain : ndarray
  206. 1-d array containing two values. If the inputs are complex, then
  207. the two returned points are the lower left and upper right corners
  208. of the smallest rectangle (aligned with the axes) in the complex
  209. plane containing the points `x`. If the inputs are real, then the
  210. two points are the ends of the smallest interval containing the
  211. points `x`.
  212. See Also
  213. --------
  214. mapparms, mapdomain
  215. Examples
  216. --------
  217. >>> from numpy.polynomial import polyutils as pu
  218. >>> points = np.arange(4)**2 - 5; points
  219. array([-5, -4, -1, 4])
  220. >>> pu.getdomain(points)
  221. array([-5., 4.])
  222. >>> c = np.exp(complex(0,1)*np.pi*np.arange(12)/6) # unit circle
  223. >>> pu.getdomain(c)
  224. array([-1.-1.j, 1.+1.j])
  225. """
  226. [x] = as_series([x], trim=False)
  227. if x.dtype.char in np.typecodes['Complex']:
  228. rmin, rmax = x.real.min(), x.real.max()
  229. imin, imax = x.imag.min(), x.imag.max()
  230. return np.array((complex(rmin, imin), complex(rmax, imax)))
  231. else:
  232. return np.array((x.min(), x.max()))
  233. def mapparms(old, new):
  234. """
  235. Linear map parameters between domains.
  236. Return the parameters of the linear map ``offset + scale*x`` that maps
  237. `old` to `new` such that ``old[i] -> new[i]``, ``i = 0, 1``.
  238. Parameters
  239. ----------
  240. old, new : array_like
  241. Domains. Each domain must (successfully) convert to a 1-d array
  242. containing precisely two values.
  243. Returns
  244. -------
  245. offset, scale : scalars
  246. The map ``L(x) = offset + scale*x`` maps the first domain to the
  247. second.
  248. See Also
  249. --------
  250. getdomain, mapdomain
  251. Notes
  252. -----
  253. Also works for complex numbers, and thus can be used to calculate the
  254. parameters required to map any line in the complex plane to any other
  255. line therein.
  256. Examples
  257. --------
  258. >>> from numpy import polynomial as P
  259. >>> P.mapparms((-1,1),(-1,1))
  260. (0.0, 1.0)
  261. >>> P.mapparms((1,-1),(-1,1))
  262. (0.0, -1.0)
  263. >>> i = complex(0,1)
  264. >>> P.mapparms((-i,-1),(1,i))
  265. ((1+1j), (1+0j))
  266. """
  267. oldlen = old[1] - old[0]
  268. newlen = new[1] - new[0]
  269. off = (old[1]*new[0] - old[0]*new[1])/oldlen
  270. scl = newlen/oldlen
  271. return off, scl
  272. def mapdomain(x, old, new):
  273. """
  274. Apply linear map to input points.
  275. The linear map ``offset + scale*x`` that maps the domain `old` to
  276. the domain `new` is applied to the points `x`.
  277. Parameters
  278. ----------
  279. x : array_like
  280. Points to be mapped. If `x` is a subtype of ndarray the subtype
  281. will be preserved.
  282. old, new : array_like
  283. The two domains that determine the map. Each must (successfully)
  284. convert to 1-d arrays containing precisely two values.
  285. Returns
  286. -------
  287. x_out : ndarray
  288. Array of points of the same shape as `x`, after application of the
  289. linear map between the two domains.
  290. See Also
  291. --------
  292. getdomain, mapparms
  293. Notes
  294. -----
  295. Effectively, this implements:
  296. .. math ::
  297. x\\_out = new[0] + m(x - old[0])
  298. where
  299. .. math ::
  300. m = \\frac{new[1]-new[0]}{old[1]-old[0]}
  301. Examples
  302. --------
  303. >>> from numpy import polynomial as P
  304. >>> old_domain = (-1,1)
  305. >>> new_domain = (0,2*np.pi)
  306. >>> x = np.linspace(-1,1,6); x
  307. array([-1. , -0.6, -0.2, 0.2, 0.6, 1. ])
  308. >>> x_out = P.mapdomain(x, old_domain, new_domain); x_out
  309. array([ 0. , 1.25663706, 2.51327412, 3.76991118, 5.02654825,
  310. 6.28318531])
  311. >>> x - P.mapdomain(x_out, new_domain, old_domain)
  312. array([ 0., 0., 0., 0., 0., 0.])
  313. Also works for complex numbers (and thus can be used to map any line in
  314. the complex plane to any other line therein).
  315. >>> i = complex(0,1)
  316. >>> old = (-1 - i, 1 + i)
  317. >>> new = (-1 + i, 1 - i)
  318. >>> z = np.linspace(old[0], old[1], 6); z
  319. 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 ])
  320. >>> new_z = P.mapdomain(z, old, new); new_z
  321. 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 ])
  322. """
  323. x = np.asanyarray(x)
  324. off, scl = mapparms(old, new)
  325. return off + scl*x