polyutils.py 23 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796
  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. import operator
  33. import functools
  34. import warnings
  35. import numpy as np
  36. __all__ = [
  37. 'RankWarning', 'PolyError', 'PolyDomainError', 'as_series', 'trimseq',
  38. 'trimcoef', 'getdomain', 'mapdomain', 'mapparms', 'PolyBase']
  39. #
  40. # Warnings and Exceptions
  41. #
  42. class RankWarning(UserWarning):
  43. """Issued by chebfit when the design matrix is rank deficient."""
  44. pass
  45. class PolyError(Exception):
  46. """Base class for errors in this module."""
  47. pass
  48. class PolyDomainError(PolyError):
  49. """Issued by the generic Poly class when two domains don't match.
  50. This is raised when an binary operation is passed Poly objects with
  51. different domains.
  52. """
  53. pass
  54. #
  55. # Base class for all polynomial types
  56. #
  57. class PolyBase:
  58. """
  59. Base class for all polynomial types.
  60. Deprecated in numpy 1.9.0, use the abstract
  61. ABCPolyBase class instead. Note that the latter
  62. requires a number of virtual functions to be
  63. implemented.
  64. """
  65. pass
  66. #
  67. # Helper functions to convert inputs to 1-D arrays
  68. #
  69. def trimseq(seq):
  70. """Remove small Poly series coefficients.
  71. Parameters
  72. ----------
  73. seq : sequence
  74. Sequence of Poly series coefficients. This routine fails for
  75. empty sequences.
  76. Returns
  77. -------
  78. series : sequence
  79. Subsequence with trailing zeros removed. If the resulting sequence
  80. would be empty, return the first element. The returned sequence may
  81. or may not be a view.
  82. Notes
  83. -----
  84. Do not lose the type info if the sequence contains unknown objects.
  85. """
  86. if len(seq) == 0:
  87. return seq
  88. else:
  89. for i in range(len(seq) - 1, -1, -1):
  90. if seq[i] != 0:
  91. break
  92. return seq[:i+1]
  93. def as_series(alist, trim=True):
  94. """
  95. Return argument as a list of 1-d arrays.
  96. The returned list contains array(s) of dtype double, complex double, or
  97. object. A 1-d argument of shape ``(N,)`` is parsed into ``N`` arrays of
  98. size one; a 2-d argument of shape ``(M,N)`` is parsed into ``M`` arrays
  99. of size ``N`` (i.e., is "parsed by row"); and a higher dimensional array
  100. raises a Value Error if it is not first reshaped into either a 1-d or 2-d
  101. array.
  102. Parameters
  103. ----------
  104. alist : array_like
  105. A 1- or 2-d array_like
  106. trim : boolean, optional
  107. When True, trailing zeros are removed from the inputs.
  108. When False, the inputs are passed through intact.
  109. Returns
  110. -------
  111. [a1, a2,...] : list of 1-D arrays
  112. A copy of the input data as a list of 1-d arrays.
  113. Raises
  114. ------
  115. ValueError
  116. Raised when `as_series` cannot convert its input to 1-d arrays, or at
  117. least one of the resulting arrays is empty.
  118. Examples
  119. --------
  120. >>> from numpy.polynomial import polyutils as pu
  121. >>> a = np.arange(4)
  122. >>> pu.as_series(a)
  123. [array([0.]), array([1.]), array([2.]), array([3.])]
  124. >>> b = np.arange(6).reshape((2,3))
  125. >>> pu.as_series(b)
  126. [array([0., 1., 2.]), array([3., 4., 5.])]
  127. >>> pu.as_series((1, np.arange(3), np.arange(2, dtype=np.float16)))
  128. [array([1.]), array([0., 1., 2.]), array([0., 1.])]
  129. >>> pu.as_series([2, [1.1, 0.]])
  130. [array([2.]), array([1.1])]
  131. >>> pu.as_series([2, [1.1, 0.]], trim=False)
  132. [array([2.]), array([1.1, 0. ])]
  133. """
  134. arrays = [np.array(a, ndmin=1, copy=False) for a in alist]
  135. if min([a.size for a in arrays]) == 0:
  136. raise ValueError("Coefficient array is empty")
  137. if any(a.ndim != 1 for a in arrays):
  138. raise ValueError("Coefficient array is not 1-d")
  139. if trim:
  140. arrays = [trimseq(a) for a in arrays]
  141. if any(a.dtype == np.dtype(object) for a in arrays):
  142. ret = []
  143. for a in arrays:
  144. if a.dtype != np.dtype(object):
  145. tmp = np.empty(len(a), dtype=np.dtype(object))
  146. tmp[:] = a[:]
  147. ret.append(tmp)
  148. else:
  149. ret.append(a.copy())
  150. else:
  151. try:
  152. dtype = np.common_type(*arrays)
  153. except Exception as e:
  154. raise ValueError("Coefficient arrays have no common type") from e
  155. ret = [np.array(a, copy=True, dtype=dtype) for a in arrays]
  156. return ret
  157. def trimcoef(c, tol=0):
  158. """
  159. Remove "small" "trailing" coefficients from a polynomial.
  160. "Small" means "small in absolute value" and is controlled by the
  161. parameter `tol`; "trailing" means highest order coefficient(s), e.g., in
  162. ``[0, 1, 1, 0, 0]`` (which represents ``0 + x + x**2 + 0*x**3 + 0*x**4``)
  163. both the 3-rd and 4-th order coefficients would be "trimmed."
  164. Parameters
  165. ----------
  166. c : array_like
  167. 1-d array of coefficients, ordered from lowest order to highest.
  168. tol : number, optional
  169. Trailing (i.e., highest order) elements with absolute value less
  170. than or equal to `tol` (default value is zero) are removed.
  171. Returns
  172. -------
  173. trimmed : ndarray
  174. 1-d array with trailing zeros removed. If the resulting series
  175. would be empty, a series containing a single zero is returned.
  176. Raises
  177. ------
  178. ValueError
  179. If `tol` < 0
  180. See Also
  181. --------
  182. trimseq
  183. Examples
  184. --------
  185. >>> from numpy.polynomial import polyutils as pu
  186. >>> pu.trimcoef((0,0,3,0,5,0,0))
  187. array([0., 0., 3., 0., 5.])
  188. >>> pu.trimcoef((0,0,1e-3,0,1e-5,0,0),1e-3) # item == tol is trimmed
  189. array([0.])
  190. >>> i = complex(0,1) # works for complex
  191. >>> pu.trimcoef((3e-4,1e-3*(1-i),5e-4,2e-5*(1+i)), 1e-3)
  192. array([0.0003+0.j , 0.001 -0.001j])
  193. """
  194. if tol < 0:
  195. raise ValueError("tol must be non-negative")
  196. [c] = as_series([c])
  197. [ind] = np.nonzero(np.abs(c) > tol)
  198. if len(ind) == 0:
  199. return c[:1]*0
  200. else:
  201. return c[:ind[-1] + 1].copy()
  202. def getdomain(x):
  203. """
  204. Return a domain suitable for given abscissae.
  205. Find a domain suitable for a polynomial or Chebyshev series
  206. defined at the values supplied.
  207. Parameters
  208. ----------
  209. x : array_like
  210. 1-d array of abscissae whose domain will be determined.
  211. Returns
  212. -------
  213. domain : ndarray
  214. 1-d array containing two values. If the inputs are complex, then
  215. the two returned points are the lower left and upper right corners
  216. of the smallest rectangle (aligned with the axes) in the complex
  217. plane containing the points `x`. If the inputs are real, then the
  218. two points are the ends of the smallest interval containing the
  219. points `x`.
  220. See Also
  221. --------
  222. mapparms, mapdomain
  223. Examples
  224. --------
  225. >>> from numpy.polynomial import polyutils as pu
  226. >>> points = np.arange(4)**2 - 5; points
  227. array([-5, -4, -1, 4])
  228. >>> pu.getdomain(points)
  229. array([-5., 4.])
  230. >>> c = np.exp(complex(0,1)*np.pi*np.arange(12)/6) # unit circle
  231. >>> pu.getdomain(c)
  232. array([-1.-1.j, 1.+1.j])
  233. """
  234. [x] = as_series([x], trim=False)
  235. if x.dtype.char in np.typecodes['Complex']:
  236. rmin, rmax = x.real.min(), x.real.max()
  237. imin, imax = x.imag.min(), x.imag.max()
  238. return np.array((complex(rmin, imin), complex(rmax, imax)))
  239. else:
  240. return np.array((x.min(), x.max()))
  241. def mapparms(old, new):
  242. """
  243. Linear map parameters between domains.
  244. Return the parameters of the linear map ``offset + scale*x`` that maps
  245. `old` to `new` such that ``old[i] -> new[i]``, ``i = 0, 1``.
  246. Parameters
  247. ----------
  248. old, new : array_like
  249. Domains. Each domain must (successfully) convert to a 1-d array
  250. containing precisely two values.
  251. Returns
  252. -------
  253. offset, scale : scalars
  254. The map ``L(x) = offset + scale*x`` maps the first domain to the
  255. second.
  256. See Also
  257. --------
  258. getdomain, mapdomain
  259. Notes
  260. -----
  261. Also works for complex numbers, and thus can be used to calculate the
  262. parameters required to map any line in the complex plane to any other
  263. line therein.
  264. Examples
  265. --------
  266. >>> from numpy.polynomial import polyutils as pu
  267. >>> pu.mapparms((-1,1),(-1,1))
  268. (0.0, 1.0)
  269. >>> pu.mapparms((1,-1),(-1,1))
  270. (-0.0, -1.0)
  271. >>> i = complex(0,1)
  272. >>> pu.mapparms((-i,-1),(1,i))
  273. ((1+1j), (1-0j))
  274. """
  275. oldlen = old[1] - old[0]
  276. newlen = new[1] - new[0]
  277. off = (old[1]*new[0] - old[0]*new[1])/oldlen
  278. scl = newlen/oldlen
  279. return off, scl
  280. def mapdomain(x, old, new):
  281. """
  282. Apply linear map to input points.
  283. The linear map ``offset + scale*x`` that maps the domain `old` to
  284. the domain `new` is applied to the points `x`.
  285. Parameters
  286. ----------
  287. x : array_like
  288. Points to be mapped. If `x` is a subtype of ndarray the subtype
  289. will be preserved.
  290. old, new : array_like
  291. The two domains that determine the map. Each must (successfully)
  292. convert to 1-d arrays containing precisely two values.
  293. Returns
  294. -------
  295. x_out : ndarray
  296. Array of points of the same shape as `x`, after application of the
  297. linear map between the two domains.
  298. See Also
  299. --------
  300. getdomain, mapparms
  301. Notes
  302. -----
  303. Effectively, this implements:
  304. .. math ::
  305. x\\_out = new[0] + m(x - old[0])
  306. where
  307. .. math ::
  308. m = \\frac{new[1]-new[0]}{old[1]-old[0]}
  309. Examples
  310. --------
  311. >>> from numpy.polynomial import polyutils as pu
  312. >>> old_domain = (-1,1)
  313. >>> new_domain = (0,2*np.pi)
  314. >>> x = np.linspace(-1,1,6); x
  315. array([-1. , -0.6, -0.2, 0.2, 0.6, 1. ])
  316. >>> x_out = pu.mapdomain(x, old_domain, new_domain); x_out
  317. array([ 0. , 1.25663706, 2.51327412, 3.76991118, 5.02654825, # may vary
  318. 6.28318531])
  319. >>> x - pu.mapdomain(x_out, new_domain, old_domain)
  320. array([0., 0., 0., 0., 0., 0.])
  321. Also works for complex numbers (and thus can be used to map any line in
  322. the complex plane to any other line therein).
  323. >>> i = complex(0,1)
  324. >>> old = (-1 - i, 1 + i)
  325. >>> new = (-1 + i, 1 - i)
  326. >>> z = np.linspace(old[0], old[1], 6); z
  327. array([-1. -1.j , -0.6-0.6j, -0.2-0.2j, 0.2+0.2j, 0.6+0.6j, 1. +1.j ])
  328. >>> new_z = pu.mapdomain(z, old, new); new_z
  329. 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 ]) # may vary
  330. """
  331. x = np.asanyarray(x)
  332. off, scl = mapparms(old, new)
  333. return off + scl*x
  334. def _nth_slice(i, ndim):
  335. sl = [np.newaxis] * ndim
  336. sl[i] = slice(None)
  337. return tuple(sl)
  338. def _vander_nd(vander_fs, points, degrees):
  339. r"""
  340. A generalization of the Vandermonde matrix for N dimensions
  341. The result is built by combining the results of 1d Vandermonde matrices,
  342. .. math::
  343. W[i_0, \ldots, i_M, j_0, \ldots, j_N] = \prod_{k=0}^N{V_k(x_k)[i_0, \ldots, i_M, j_k]}
  344. where
  345. .. math::
  346. N &= \texttt{len(points)} = \texttt{len(degrees)} = \texttt{len(vander\_fs)} \\
  347. M &= \texttt{points[k].ndim} \\
  348. V_k &= \texttt{vander\_fs[k]} \\
  349. x_k &= \texttt{points[k]} \\
  350. 0 \le j_k &\le \texttt{degrees[k]}
  351. Expanding the one-dimensional :math:`V_k` functions gives:
  352. .. math::
  353. W[i_0, \ldots, i_M, j_0, \ldots, j_N] = \prod_{k=0}^N{B_{k, j_k}(x_k[i_0, \ldots, i_M])}
  354. where :math:`B_{k,m}` is the m'th basis of the polynomial construction used along
  355. dimension :math:`k`. For a regular polynomial, :math:`B_{k, m}(x) = P_m(x) = x^m`.
  356. Parameters
  357. ----------
  358. vander_fs : Sequence[function(array_like, int) -> ndarray]
  359. The 1d vander function to use for each axis, such as ``polyvander``
  360. points : Sequence[array_like]
  361. Arrays of point coordinates, all of the same shape. The dtypes
  362. will be converted to either float64 or complex128 depending on
  363. whether any of the elements are complex. Scalars are converted to
  364. 1-D arrays.
  365. This must be the same length as `vander_fs`.
  366. degrees : Sequence[int]
  367. The maximum degree (inclusive) to use for each axis.
  368. This must be the same length as `vander_fs`.
  369. Returns
  370. -------
  371. vander_nd : ndarray
  372. An array of shape ``points[0].shape + tuple(d + 1 for d in degrees)``.
  373. """
  374. n_dims = len(vander_fs)
  375. if n_dims != len(points):
  376. raise ValueError(
  377. f"Expected {n_dims} dimensions of sample points, got {len(points)}")
  378. if n_dims != len(degrees):
  379. raise ValueError(
  380. f"Expected {n_dims} dimensions of degrees, got {len(degrees)}")
  381. if n_dims == 0:
  382. raise ValueError("Unable to guess a dtype or shape when no points are given")
  383. # convert to the same shape and type
  384. points = tuple(np.array(tuple(points), copy=False) + 0.0)
  385. # produce the vandermonde matrix for each dimension, placing the last
  386. # axis of each in an independent trailing axis of the output
  387. vander_arrays = (
  388. vander_fs[i](points[i], degrees[i])[(...,) + _nth_slice(i, n_dims)]
  389. for i in range(n_dims)
  390. )
  391. # we checked this wasn't empty already, so no `initial` needed
  392. return functools.reduce(operator.mul, vander_arrays)
  393. def _vander_nd_flat(vander_fs, points, degrees):
  394. """
  395. Like `_vander_nd`, but flattens the last ``len(degrees)`` axes into a single axis
  396. Used to implement the public ``<type>vander<n>d`` functions.
  397. """
  398. v = _vander_nd(vander_fs, points, degrees)
  399. return v.reshape(v.shape[:-len(degrees)] + (-1,))
  400. def _fromroots(line_f, mul_f, roots):
  401. """
  402. Helper function used to implement the ``<type>fromroots`` functions.
  403. Parameters
  404. ----------
  405. line_f : function(float, float) -> ndarray
  406. The ``<type>line`` function, such as ``polyline``
  407. mul_f : function(array_like, array_like) -> ndarray
  408. The ``<type>mul`` function, such as ``polymul``
  409. roots :
  410. See the ``<type>fromroots`` functions for more detail
  411. """
  412. if len(roots) == 0:
  413. return np.ones(1)
  414. else:
  415. [roots] = as_series([roots], trim=False)
  416. roots.sort()
  417. p = [line_f(-r, 1) for r in roots]
  418. n = len(p)
  419. while n > 1:
  420. m, r = divmod(n, 2)
  421. tmp = [mul_f(p[i], p[i+m]) for i in range(m)]
  422. if r:
  423. tmp[0] = mul_f(tmp[0], p[-1])
  424. p = tmp
  425. n = m
  426. return p[0]
  427. def _valnd(val_f, c, *args):
  428. """
  429. Helper function used to implement the ``<type>val<n>d`` functions.
  430. Parameters
  431. ----------
  432. val_f : function(array_like, array_like, tensor: bool) -> array_like
  433. The ``<type>val`` function, such as ``polyval``
  434. c, args :
  435. See the ``<type>val<n>d`` functions for more detail
  436. """
  437. args = [np.asanyarray(a) for a in args]
  438. shape0 = args[0].shape
  439. if not all((a.shape == shape0 for a in args[1:])):
  440. if len(args) == 3:
  441. raise ValueError('x, y, z are incompatible')
  442. elif len(args) == 2:
  443. raise ValueError('x, y are incompatible')
  444. else:
  445. raise ValueError('ordinates are incompatible')
  446. it = iter(args)
  447. x0 = next(it)
  448. # use tensor on only the first
  449. c = val_f(x0, c)
  450. for xi in it:
  451. c = val_f(xi, c, tensor=False)
  452. return c
  453. def _gridnd(val_f, c, *args):
  454. """
  455. Helper function used to implement the ``<type>grid<n>d`` functions.
  456. Parameters
  457. ----------
  458. val_f : function(array_like, array_like, tensor: bool) -> array_like
  459. The ``<type>val`` function, such as ``polyval``
  460. c, args :
  461. See the ``<type>grid<n>d`` functions for more detail
  462. """
  463. for xi in args:
  464. c = val_f(xi, c)
  465. return c
  466. def _div(mul_f, c1, c2):
  467. """
  468. Helper function used to implement the ``<type>div`` functions.
  469. Implementation uses repeated subtraction of c2 multiplied by the nth basis.
  470. For some polynomial types, a more efficient approach may be possible.
  471. Parameters
  472. ----------
  473. mul_f : function(array_like, array_like) -> array_like
  474. The ``<type>mul`` function, such as ``polymul``
  475. c1, c2 :
  476. See the ``<type>div`` functions for more detail
  477. """
  478. # c1, c2 are trimmed copies
  479. [c1, c2] = as_series([c1, c2])
  480. if c2[-1] == 0:
  481. raise ZeroDivisionError()
  482. lc1 = len(c1)
  483. lc2 = len(c2)
  484. if lc1 < lc2:
  485. return c1[:1]*0, c1
  486. elif lc2 == 1:
  487. return c1/c2[-1], c1[:1]*0
  488. else:
  489. quo = np.empty(lc1 - lc2 + 1, dtype=c1.dtype)
  490. rem = c1
  491. for i in range(lc1 - lc2, - 1, -1):
  492. p = mul_f([0]*i + [1], c2)
  493. q = rem[-1]/p[-1]
  494. rem = rem[:-1] - q*p[:-1]
  495. quo[i] = q
  496. return quo, trimseq(rem)
  497. def _add(c1, c2):
  498. """ Helper function used to implement the ``<type>add`` functions. """
  499. # c1, c2 are trimmed copies
  500. [c1, c2] = as_series([c1, c2])
  501. if len(c1) > len(c2):
  502. c1[:c2.size] += c2
  503. ret = c1
  504. else:
  505. c2[:c1.size] += c1
  506. ret = c2
  507. return trimseq(ret)
  508. def _sub(c1, c2):
  509. """ Helper function used to implement the ``<type>sub`` functions. """
  510. # c1, c2 are trimmed copies
  511. [c1, c2] = as_series([c1, c2])
  512. if len(c1) > len(c2):
  513. c1[:c2.size] -= c2
  514. ret = c1
  515. else:
  516. c2 = -c2
  517. c2[:c1.size] += c1
  518. ret = c2
  519. return trimseq(ret)
  520. def _fit(vander_f, x, y, deg, rcond=None, full=False, w=None):
  521. """
  522. Helper function used to implement the ``<type>fit`` functions.
  523. Parameters
  524. ----------
  525. vander_f : function(array_like, int) -> ndarray
  526. The 1d vander function, such as ``polyvander``
  527. c1, c2 :
  528. See the ``<type>fit`` functions for more detail
  529. """
  530. x = np.asarray(x) + 0.0
  531. y = np.asarray(y) + 0.0
  532. deg = np.asarray(deg)
  533. # check arguments.
  534. if deg.ndim > 1 or deg.dtype.kind not in 'iu' or deg.size == 0:
  535. raise TypeError("deg must be an int or non-empty 1-D array of int")
  536. if deg.min() < 0:
  537. raise ValueError("expected deg >= 0")
  538. if x.ndim != 1:
  539. raise TypeError("expected 1D vector for x")
  540. if x.size == 0:
  541. raise TypeError("expected non-empty vector for x")
  542. if y.ndim < 1 or y.ndim > 2:
  543. raise TypeError("expected 1D or 2D array for y")
  544. if len(x) != len(y):
  545. raise TypeError("expected x and y to have same length")
  546. if deg.ndim == 0:
  547. lmax = deg
  548. order = lmax + 1
  549. van = vander_f(x, lmax)
  550. else:
  551. deg = np.sort(deg)
  552. lmax = deg[-1]
  553. order = len(deg)
  554. van = vander_f(x, lmax)[:, deg]
  555. # set up the least squares matrices in transposed form
  556. lhs = van.T
  557. rhs = y.T
  558. if w is not None:
  559. w = np.asarray(w) + 0.0
  560. if w.ndim != 1:
  561. raise TypeError("expected 1D vector for w")
  562. if len(x) != len(w):
  563. raise TypeError("expected x and w to have same length")
  564. # apply weights. Don't use inplace operations as they
  565. # can cause problems with NA.
  566. lhs = lhs * w
  567. rhs = rhs * w
  568. # set rcond
  569. if rcond is None:
  570. rcond = len(x)*np.finfo(x.dtype).eps
  571. # Determine the norms of the design matrix columns.
  572. if issubclass(lhs.dtype.type, np.complexfloating):
  573. scl = np.sqrt((np.square(lhs.real) + np.square(lhs.imag)).sum(1))
  574. else:
  575. scl = np.sqrt(np.square(lhs).sum(1))
  576. scl[scl == 0] = 1
  577. # Solve the least squares problem.
  578. c, resids, rank, s = np.linalg.lstsq(lhs.T/scl, rhs.T, rcond)
  579. c = (c.T/scl).T
  580. # Expand c to include non-fitted coefficients which are set to zero
  581. if deg.ndim > 0:
  582. if c.ndim == 2:
  583. cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype)
  584. else:
  585. cc = np.zeros(lmax+1, dtype=c.dtype)
  586. cc[deg] = c
  587. c = cc
  588. # warn on rank reduction
  589. if rank != order and not full:
  590. msg = "The fit may be poorly conditioned"
  591. warnings.warn(msg, RankWarning, stacklevel=2)
  592. if full:
  593. return c, [resids, rank, s, rcond]
  594. else:
  595. return c
  596. def _pow(mul_f, c, pow, maxpower):
  597. """
  598. Helper function used to implement the ``<type>pow`` functions.
  599. Parameters
  600. ----------
  601. vander_f : function(array_like, int) -> ndarray
  602. The 1d vander function, such as ``polyvander``
  603. pow, maxpower :
  604. See the ``<type>pow`` functions for more detail
  605. mul_f : function(array_like, array_like) -> ndarray
  606. The ``<type>mul`` function, such as ``polymul``
  607. """
  608. # c is a trimmed copy
  609. [c] = as_series([c])
  610. power = int(pow)
  611. if power != pow or power < 0:
  612. raise ValueError("Power must be a non-negative integer.")
  613. elif maxpower is not None and power > maxpower:
  614. raise ValueError("Power is too large")
  615. elif power == 0:
  616. return np.array([1], dtype=c.dtype)
  617. elif power == 1:
  618. return c
  619. else:
  620. # This can be made more efficient by using powers of two
  621. # in the usual way.
  622. prd = c
  623. for i in range(2, power + 1):
  624. prd = mul_f(prd, c)
  625. return prd
  626. def _deprecate_as_int(x, desc):
  627. """
  628. Like `operator.index`, but emits a deprecation warning when passed a float
  629. Parameters
  630. ----------
  631. x : int-like, or float with integral value
  632. Value to interpret as an integer
  633. desc : str
  634. description to include in any error message
  635. Raises
  636. ------
  637. TypeError : if x is a non-integral float or non-numeric
  638. DeprecationWarning : if x is an integral float
  639. """
  640. try:
  641. return operator.index(x)
  642. except TypeError as e:
  643. # Numpy 1.17.0, 2019-03-11
  644. try:
  645. ix = int(x)
  646. except TypeError:
  647. pass
  648. else:
  649. if ix == x:
  650. warnings.warn(
  651. f"In future, this will raise TypeError, as {desc} will "
  652. "need to be an integer not just an integral float.",
  653. DeprecationWarning,
  654. stacklevel=3
  655. )
  656. return ix
  657. raise TypeError(f"{desc} must be an integer") from e