Coverage for /var/srv/projects/api.amasfac.comuna18.com/tmp/venv/lib/python3.9/site-packages/numpy/polynomial/polyutils.py: 8%
229 statements
« prev ^ index » next coverage.py v6.4.4, created at 2023-07-17 14:22 -0600
1"""
2Utility classes and functions for the polynomial modules.
4This module provides: error and warning objects; a polynomial base class;
5and some routines used in both the `polynomial` and `chebyshev` modules.
7Warning objects
8---------------
10.. autosummary::
11 :toctree: generated/
13 RankWarning raised in least-squares fit for rank-deficient matrix.
15Functions
16---------
18.. autosummary::
19 :toctree: generated/
21 as_series convert list of array_likes into 1-D arrays of common type.
22 trimseq remove trailing zeros.
23 trimcoef remove small trailing coefficients.
24 getdomain return the domain appropriate for a given set of abscissae.
25 mapdomain maps points between domains.
26 mapparms parameters of the linear map between domains.
28"""
29import operator
30import functools
31import warnings
33import numpy as np
35__all__ = [
36 'RankWarning', 'as_series', 'trimseq',
37 'trimcoef', 'getdomain', 'mapdomain', 'mapparms']
39#
40# Warnings and Exceptions
41#
43class RankWarning(UserWarning):
44 """Issued by chebfit when the design matrix is rank deficient."""
45 pass
47#
48# Helper functions to convert inputs to 1-D arrays
49#
50def trimseq(seq):
51 """Remove small Poly series coefficients.
53 Parameters
54 ----------
55 seq : sequence
56 Sequence of Poly series coefficients. This routine fails for
57 empty sequences.
59 Returns
60 -------
61 series : sequence
62 Subsequence with trailing zeros removed. If the resulting sequence
63 would be empty, return the first element. The returned sequence may
64 or may not be a view.
66 Notes
67 -----
68 Do not lose the type info if the sequence contains unknown objects.
70 """
71 if len(seq) == 0:
72 return seq
73 else:
74 for i in range(len(seq) - 1, -1, -1):
75 if seq[i] != 0:
76 break
77 return seq[:i+1]
80def as_series(alist, trim=True):
81 """
82 Return argument as a list of 1-d arrays.
84 The returned list contains array(s) of dtype double, complex double, or
85 object. A 1-d argument of shape ``(N,)`` is parsed into ``N`` arrays of
86 size one; a 2-d argument of shape ``(M,N)`` is parsed into ``M`` arrays
87 of size ``N`` (i.e., is "parsed by row"); and a higher dimensional array
88 raises a Value Error if it is not first reshaped into either a 1-d or 2-d
89 array.
91 Parameters
92 ----------
93 alist : array_like
94 A 1- or 2-d array_like
95 trim : boolean, optional
96 When True, trailing zeros are removed from the inputs.
97 When False, the inputs are passed through intact.
99 Returns
100 -------
101 [a1, a2,...] : list of 1-D arrays
102 A copy of the input data as a list of 1-d arrays.
104 Raises
105 ------
106 ValueError
107 Raised when `as_series` cannot convert its input to 1-d arrays, or at
108 least one of the resulting arrays is empty.
110 Examples
111 --------
112 >>> from numpy.polynomial import polyutils as pu
113 >>> a = np.arange(4)
114 >>> pu.as_series(a)
115 [array([0.]), array([1.]), array([2.]), array([3.])]
116 >>> b = np.arange(6).reshape((2,3))
117 >>> pu.as_series(b)
118 [array([0., 1., 2.]), array([3., 4., 5.])]
120 >>> pu.as_series((1, np.arange(3), np.arange(2, dtype=np.float16)))
121 [array([1.]), array([0., 1., 2.]), array([0., 1.])]
123 >>> pu.as_series([2, [1.1, 0.]])
124 [array([2.]), array([1.1])]
126 >>> pu.as_series([2, [1.1, 0.]], trim=False)
127 [array([2.]), array([1.1, 0. ])]
129 """
130 arrays = [np.array(a, ndmin=1, copy=False) for a in alist]
131 if min([a.size for a in arrays]) == 0:
132 raise ValueError("Coefficient array is empty")
133 if any(a.ndim != 1 for a in arrays):
134 raise ValueError("Coefficient array is not 1-d")
135 if trim:
136 arrays = [trimseq(a) for a in arrays]
138 if any(a.dtype == np.dtype(object) for a in arrays):
139 ret = []
140 for a in arrays:
141 if a.dtype != np.dtype(object):
142 tmp = np.empty(len(a), dtype=np.dtype(object))
143 tmp[:] = a[:]
144 ret.append(tmp)
145 else:
146 ret.append(a.copy())
147 else:
148 try:
149 dtype = np.common_type(*arrays)
150 except Exception as e:
151 raise ValueError("Coefficient arrays have no common type") from e
152 ret = [np.array(a, copy=True, dtype=dtype) for a in arrays]
153 return ret
156def trimcoef(c, tol=0):
157 """
158 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."
165 Parameters
166 ----------
167 c : array_like
168 1-d array of coefficients, ordered from lowest order to highest.
169 tol : number, optional
170 Trailing (i.e., highest order) elements with absolute value less
171 than or equal to `tol` (default value is zero) are removed.
173 Returns
174 -------
175 trimmed : ndarray
176 1-d array with trailing zeros removed. If the resulting series
177 would be empty, a series containing a single zero is returned.
179 Raises
180 ------
181 ValueError
182 If `tol` < 0
184 See Also
185 --------
186 trimseq
188 Examples
189 --------
190 >>> from numpy.polynomial import polyutils as pu
191 >>> pu.trimcoef((0,0,3,0,5,0,0))
192 array([0., 0., 3., 0., 5.])
193 >>> pu.trimcoef((0,0,1e-3,0,1e-5,0,0),1e-3) # item == tol is trimmed
194 array([0.])
195 >>> i = complex(0,1) # works for complex
196 >>> pu.trimcoef((3e-4,1e-3*(1-i),5e-4,2e-5*(1+i)), 1e-3)
197 array([0.0003+0.j , 0.001 -0.001j])
199 """
200 if tol < 0:
201 raise ValueError("tol must be non-negative")
203 [c] = as_series([c])
204 [ind] = np.nonzero(np.abs(c) > tol)
205 if len(ind) == 0:
206 return c[:1]*0
207 else:
208 return c[:ind[-1] + 1].copy()
210def getdomain(x):
211 """
212 Return a domain suitable for given abscissae.
214 Find a domain suitable for a polynomial or Chebyshev series
215 defined at the values supplied.
217 Parameters
218 ----------
219 x : array_like
220 1-d array of abscissae whose domain will be determined.
222 Returns
223 -------
224 domain : ndarray
225 1-d array containing two values. If the inputs are complex, then
226 the two returned points are the lower left and upper right corners
227 of the smallest rectangle (aligned with the axes) in the complex
228 plane containing the points `x`. If the inputs are real, then the
229 two points are the ends of the smallest interval containing the
230 points `x`.
232 See Also
233 --------
234 mapparms, mapdomain
236 Examples
237 --------
238 >>> from numpy.polynomial import polyutils as pu
239 >>> points = np.arange(4)**2 - 5; points
240 array([-5, -4, -1, 4])
241 >>> pu.getdomain(points)
242 array([-5., 4.])
243 >>> c = np.exp(complex(0,1)*np.pi*np.arange(12)/6) # unit circle
244 >>> pu.getdomain(c)
245 array([-1.-1.j, 1.+1.j])
247 """
248 [x] = as_series([x], trim=False)
249 if x.dtype.char in np.typecodes['Complex']:
250 rmin, rmax = x.real.min(), x.real.max()
251 imin, imax = x.imag.min(), x.imag.max()
252 return np.array((complex(rmin, imin), complex(rmax, imax)))
253 else:
254 return np.array((x.min(), x.max()))
256def mapparms(old, new):
257 """
258 Linear map parameters between domains.
260 Return the parameters of the linear map ``offset + scale*x`` that maps
261 `old` to `new` such that ``old[i] -> new[i]``, ``i = 0, 1``.
263 Parameters
264 ----------
265 old, new : array_like
266 Domains. Each domain must (successfully) convert to a 1-d array
267 containing precisely two values.
269 Returns
270 -------
271 offset, scale : scalars
272 The map ``L(x) = offset + scale*x`` maps the first domain to the
273 second.
275 See Also
276 --------
277 getdomain, mapdomain
279 Notes
280 -----
281 Also works for complex numbers, and thus can be used to calculate the
282 parameters required to map any line in the complex plane to any other
283 line therein.
285 Examples
286 --------
287 >>> from numpy.polynomial import polyutils as pu
288 >>> pu.mapparms((-1,1),(-1,1))
289 (0.0, 1.0)
290 >>> pu.mapparms((1,-1),(-1,1))
291 (-0.0, -1.0)
292 >>> i = complex(0,1)
293 >>> pu.mapparms((-i,-1),(1,i))
294 ((1+1j), (1-0j))
296 """
297 oldlen = old[1] - old[0]
298 newlen = new[1] - new[0]
299 off = (old[1]*new[0] - old[0]*new[1])/oldlen
300 scl = newlen/oldlen
301 return off, scl
303def mapdomain(x, old, new):
304 """
305 Apply linear map to input points.
307 The linear map ``offset + scale*x`` that maps the domain `old` to
308 the domain `new` is applied to the points `x`.
310 Parameters
311 ----------
312 x : array_like
313 Points to be mapped. If `x` is a subtype of ndarray the subtype
314 will be preserved.
315 old, new : array_like
316 The two domains that determine the map. Each must (successfully)
317 convert to 1-d arrays containing precisely two values.
319 Returns
320 -------
321 x_out : ndarray
322 Array of points of the same shape as `x`, after application of the
323 linear map between the two domains.
325 See Also
326 --------
327 getdomain, mapparms
329 Notes
330 -----
331 Effectively, this implements:
333 .. math::
334 x\\_out = new[0] + m(x - old[0])
336 where
338 .. math::
339 m = \\frac{new[1]-new[0]}{old[1]-old[0]}
341 Examples
342 --------
343 >>> from numpy.polynomial import polyutils as pu
344 >>> old_domain = (-1,1)
345 >>> new_domain = (0,2*np.pi)
346 >>> x = np.linspace(-1,1,6); x
347 array([-1. , -0.6, -0.2, 0.2, 0.6, 1. ])
348 >>> x_out = pu.mapdomain(x, old_domain, new_domain); x_out
349 array([ 0. , 1.25663706, 2.51327412, 3.76991118, 5.02654825, # may vary
350 6.28318531])
351 >>> x - pu.mapdomain(x_out, new_domain, old_domain)
352 array([0., 0., 0., 0., 0., 0.])
354 Also works for complex numbers (and thus can be used to map any line in
355 the complex plane to any other line therein).
357 >>> i = complex(0,1)
358 >>> old = (-1 - i, 1 + i)
359 >>> new = (-1 + i, 1 - i)
360 >>> z = np.linspace(old[0], old[1], 6); z
361 array([-1. -1.j , -0.6-0.6j, -0.2-0.2j, 0.2+0.2j, 0.6+0.6j, 1. +1.j ])
362 >>> new_z = pu.mapdomain(z, old, new); new_z
363 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
365 """
366 x = np.asanyarray(x)
367 off, scl = mapparms(old, new)
368 return off + scl*x
371def _nth_slice(i, ndim):
372 sl = [np.newaxis] * ndim
373 sl[i] = slice(None)
374 return tuple(sl)
377def _vander_nd(vander_fs, points, degrees):
378 r"""
379 A generalization of the Vandermonde matrix for N dimensions
381 The result is built by combining the results of 1d Vandermonde matrices,
383 .. math::
384 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]}
386 where
388 .. math::
389 N &= \texttt{len(points)} = \texttt{len(degrees)} = \texttt{len(vander\_fs)} \\
390 M &= \texttt{points[k].ndim} \\
391 V_k &= \texttt{vander\_fs[k]} \\
392 x_k &= \texttt{points[k]} \\
393 0 \le j_k &\le \texttt{degrees[k]}
395 Expanding the one-dimensional :math:`V_k` functions gives:
397 .. math::
398 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])}
400 where :math:`B_{k,m}` is the m'th basis of the polynomial construction used along
401 dimension :math:`k`. For a regular polynomial, :math:`B_{k, m}(x) = P_m(x) = x^m`.
403 Parameters
404 ----------
405 vander_fs : Sequence[function(array_like, int) -> ndarray]
406 The 1d vander function to use for each axis, such as ``polyvander``
407 points : Sequence[array_like]
408 Arrays of point coordinates, all of the same shape. The dtypes
409 will be converted to either float64 or complex128 depending on
410 whether any of the elements are complex. Scalars are converted to
411 1-D arrays.
412 This must be the same length as `vander_fs`.
413 degrees : Sequence[int]
414 The maximum degree (inclusive) to use for each axis.
415 This must be the same length as `vander_fs`.
417 Returns
418 -------
419 vander_nd : ndarray
420 An array of shape ``points[0].shape + tuple(d + 1 for d in degrees)``.
421 """
422 n_dims = len(vander_fs)
423 if n_dims != len(points):
424 raise ValueError(
425 f"Expected {n_dims} dimensions of sample points, got {len(points)}")
426 if n_dims != len(degrees):
427 raise ValueError(
428 f"Expected {n_dims} dimensions of degrees, got {len(degrees)}")
429 if n_dims == 0:
430 raise ValueError("Unable to guess a dtype or shape when no points are given")
432 # convert to the same shape and type
433 points = tuple(np.array(tuple(points), copy=False) + 0.0)
435 # produce the vandermonde matrix for each dimension, placing the last
436 # axis of each in an independent trailing axis of the output
437 vander_arrays = (
438 vander_fs[i](points[i], degrees[i])[(...,) + _nth_slice(i, n_dims)]
439 for i in range(n_dims)
440 )
442 # we checked this wasn't empty already, so no `initial` needed
443 return functools.reduce(operator.mul, vander_arrays)
446def _vander_nd_flat(vander_fs, points, degrees):
447 """
448 Like `_vander_nd`, but flattens the last ``len(degrees)`` axes into a single axis
450 Used to implement the public ``<type>vander<n>d`` functions.
451 """
452 v = _vander_nd(vander_fs, points, degrees)
453 return v.reshape(v.shape[:-len(degrees)] + (-1,))
456def _fromroots(line_f, mul_f, roots):
457 """
458 Helper function used to implement the ``<type>fromroots`` functions.
460 Parameters
461 ----------
462 line_f : function(float, float) -> ndarray
463 The ``<type>line`` function, such as ``polyline``
464 mul_f : function(array_like, array_like) -> ndarray
465 The ``<type>mul`` function, such as ``polymul``
466 roots
467 See the ``<type>fromroots`` functions for more detail
468 """
469 if len(roots) == 0:
470 return np.ones(1)
471 else:
472 [roots] = as_series([roots], trim=False)
473 roots.sort()
474 p = [line_f(-r, 1) for r in roots]
475 n = len(p)
476 while n > 1:
477 m, r = divmod(n, 2)
478 tmp = [mul_f(p[i], p[i+m]) for i in range(m)]
479 if r:
480 tmp[0] = mul_f(tmp[0], p[-1])
481 p = tmp
482 n = m
483 return p[0]
486def _valnd(val_f, c, *args):
487 """
488 Helper function used to implement the ``<type>val<n>d`` functions.
490 Parameters
491 ----------
492 val_f : function(array_like, array_like, tensor: bool) -> array_like
493 The ``<type>val`` function, such as ``polyval``
494 c, args
495 See the ``<type>val<n>d`` functions for more detail
496 """
497 args = [np.asanyarray(a) for a in args]
498 shape0 = args[0].shape
499 if not all((a.shape == shape0 for a in args[1:])):
500 if len(args) == 3:
501 raise ValueError('x, y, z are incompatible')
502 elif len(args) == 2:
503 raise ValueError('x, y are incompatible')
504 else:
505 raise ValueError('ordinates are incompatible')
506 it = iter(args)
507 x0 = next(it)
509 # use tensor on only the first
510 c = val_f(x0, c)
511 for xi in it:
512 c = val_f(xi, c, tensor=False)
513 return c
516def _gridnd(val_f, c, *args):
517 """
518 Helper function used to implement the ``<type>grid<n>d`` functions.
520 Parameters
521 ----------
522 val_f : function(array_like, array_like, tensor: bool) -> array_like
523 The ``<type>val`` function, such as ``polyval``
524 c, args
525 See the ``<type>grid<n>d`` functions for more detail
526 """
527 for xi in args:
528 c = val_f(xi, c)
529 return c
532def _div(mul_f, c1, c2):
533 """
534 Helper function used to implement the ``<type>div`` functions.
536 Implementation uses repeated subtraction of c2 multiplied by the nth basis.
537 For some polynomial types, a more efficient approach may be possible.
539 Parameters
540 ----------
541 mul_f : function(array_like, array_like) -> array_like
542 The ``<type>mul`` function, such as ``polymul``
543 c1, c2
544 See the ``<type>div`` functions for more detail
545 """
546 # c1, c2 are trimmed copies
547 [c1, c2] = as_series([c1, c2])
548 if c2[-1] == 0:
549 raise ZeroDivisionError()
551 lc1 = len(c1)
552 lc2 = len(c2)
553 if lc1 < lc2:
554 return c1[:1]*0, c1
555 elif lc2 == 1:
556 return c1/c2[-1], c1[:1]*0
557 else:
558 quo = np.empty(lc1 - lc2 + 1, dtype=c1.dtype)
559 rem = c1
560 for i in range(lc1 - lc2, - 1, -1):
561 p = mul_f([0]*i + [1], c2)
562 q = rem[-1]/p[-1]
563 rem = rem[:-1] - q*p[:-1]
564 quo[i] = q
565 return quo, trimseq(rem)
568def _add(c1, c2):
569 """ Helper function used to implement the ``<type>add`` functions. """
570 # c1, c2 are trimmed copies
571 [c1, c2] = as_series([c1, c2])
572 if len(c1) > len(c2):
573 c1[:c2.size] += c2
574 ret = c1
575 else:
576 c2[:c1.size] += c1
577 ret = c2
578 return trimseq(ret)
581def _sub(c1, c2):
582 """ Helper function used to implement the ``<type>sub`` functions. """
583 # c1, c2 are trimmed copies
584 [c1, c2] = as_series([c1, c2])
585 if len(c1) > len(c2):
586 c1[:c2.size] -= c2
587 ret = c1
588 else:
589 c2 = -c2
590 c2[:c1.size] += c1
591 ret = c2
592 return trimseq(ret)
595def _fit(vander_f, x, y, deg, rcond=None, full=False, w=None):
596 """
597 Helper function used to implement the ``<type>fit`` functions.
599 Parameters
600 ----------
601 vander_f : function(array_like, int) -> ndarray
602 The 1d vander function, such as ``polyvander``
603 c1, c2
604 See the ``<type>fit`` functions for more detail
605 """
606 x = np.asarray(x) + 0.0
607 y = np.asarray(y) + 0.0
608 deg = np.asarray(deg)
610 # check arguments.
611 if deg.ndim > 1 or deg.dtype.kind not in 'iu' or deg.size == 0:
612 raise TypeError("deg must be an int or non-empty 1-D array of int")
613 if deg.min() < 0:
614 raise ValueError("expected deg >= 0")
615 if x.ndim != 1:
616 raise TypeError("expected 1D vector for x")
617 if x.size == 0:
618 raise TypeError("expected non-empty vector for x")
619 if y.ndim < 1 or y.ndim > 2:
620 raise TypeError("expected 1D or 2D array for y")
621 if len(x) != len(y):
622 raise TypeError("expected x and y to have same length")
624 if deg.ndim == 0:
625 lmax = deg
626 order = lmax + 1
627 van = vander_f(x, lmax)
628 else:
629 deg = np.sort(deg)
630 lmax = deg[-1]
631 order = len(deg)
632 van = vander_f(x, lmax)[:, deg]
634 # set up the least squares matrices in transposed form
635 lhs = van.T
636 rhs = y.T
637 if w is not None:
638 w = np.asarray(w) + 0.0
639 if w.ndim != 1:
640 raise TypeError("expected 1D vector for w")
641 if len(x) != len(w):
642 raise TypeError("expected x and w to have same length")
643 # apply weights. Don't use inplace operations as they
644 # can cause problems with NA.
645 lhs = lhs * w
646 rhs = rhs * w
648 # set rcond
649 if rcond is None:
650 rcond = len(x)*np.finfo(x.dtype).eps
652 # Determine the norms of the design matrix columns.
653 if issubclass(lhs.dtype.type, np.complexfloating):
654 scl = np.sqrt((np.square(lhs.real) + np.square(lhs.imag)).sum(1))
655 else:
656 scl = np.sqrt(np.square(lhs).sum(1))
657 scl[scl == 0] = 1
659 # Solve the least squares problem.
660 c, resids, rank, s = np.linalg.lstsq(lhs.T/scl, rhs.T, rcond)
661 c = (c.T/scl).T
663 # Expand c to include non-fitted coefficients which are set to zero
664 if deg.ndim > 0:
665 if c.ndim == 2:
666 cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype)
667 else:
668 cc = np.zeros(lmax+1, dtype=c.dtype)
669 cc[deg] = c
670 c = cc
672 # warn on rank reduction
673 if rank != order and not full:
674 msg = "The fit may be poorly conditioned"
675 warnings.warn(msg, RankWarning, stacklevel=2)
677 if full:
678 return c, [resids, rank, s, rcond]
679 else:
680 return c
683def _pow(mul_f, c, pow, maxpower):
684 """
685 Helper function used to implement the ``<type>pow`` functions.
687 Parameters
688 ----------
689 mul_f : function(array_like, array_like) -> ndarray
690 The ``<type>mul`` function, such as ``polymul``
691 c : array_like
692 1-D array of array of series coefficients
693 pow, maxpower
694 See the ``<type>pow`` functions for more detail
695 """
696 # c is a trimmed copy
697 [c] = as_series([c])
698 power = int(pow)
699 if power != pow or power < 0:
700 raise ValueError("Power must be a non-negative integer.")
701 elif maxpower is not None and power > maxpower:
702 raise ValueError("Power is too large")
703 elif power == 0:
704 return np.array([1], dtype=c.dtype)
705 elif power == 1:
706 return c
707 else:
708 # This can be made more efficient by using powers of two
709 # in the usual way.
710 prd = c
711 for i in range(2, power + 1):
712 prd = mul_f(prd, c)
713 return prd
716def _deprecate_as_int(x, desc):
717 """
718 Like `operator.index`, but emits a deprecation warning when passed a float
720 Parameters
721 ----------
722 x : int-like, or float with integral value
723 Value to interpret as an integer
724 desc : str
725 description to include in any error message
727 Raises
728 ------
729 TypeError : if x is a non-integral float or non-numeric
730 DeprecationWarning : if x is an integral float
731 """
732 try:
733 return operator.index(x)
734 except TypeError as e:
735 # Numpy 1.17.0, 2019-03-11
736 try:
737 ix = int(x)
738 except TypeError:
739 pass
740 else:
741 if ix == x:
742 warnings.warn(
743 f"In future, this will raise TypeError, as {desc} will "
744 "need to be an integer not just an integral float.",
745 DeprecationWarning,
746 stacklevel=3
747 )
748 return ix
750 raise TypeError(f"{desc} must be an integer") from e