summaryrefslogtreecommitdiff
path: root/numpy/polynomial
diff options
context:
space:
mode:
Diffstat (limited to 'numpy/polynomial')
-rw-r--r--numpy/polynomial/_polybase.py77
-rw-r--r--numpy/polynomial/chebyshev.py217
-rw-r--r--numpy/polynomial/hermite.py235
-rw-r--r--numpy/polynomial/hermite_e.py235
-rw-r--r--numpy/polynomial/laguerre.py234
-rw-r--r--numpy/polynomial/legendre.py235
-rw-r--r--numpy/polynomial/polynomial.py229
-rw-r--r--numpy/polynomial/polyutils.py349
-rw-r--r--numpy/polynomial/tests/test_chebyshev.py6
-rw-r--r--numpy/polynomial/tests/test_hermite.py6
-rw-r--r--numpy/polynomial/tests/test_hermite_e.py6
-rw-r--r--numpy/polynomial/tests/test_laguerre.py6
-rw-r--r--numpy/polynomial/tests/test_legendre.py6
-rw-r--r--numpy/polynomial/tests/test_polynomial.py6
14 files changed, 562 insertions, 1285 deletions
diff --git a/numpy/polynomial/_polybase.py b/numpy/polynomial/_polybase.py
index c28e77e69..bfa030714 100644
--- a/numpy/polynomial/_polybase.py
+++ b/numpy/polynomial/_polybase.py
@@ -8,7 +8,7 @@ abc module from the stdlib, hence it is only available for Python >= 2.6.
"""
from __future__ import division, absolute_import, print_function
-from abc import ABCMeta, abstractmethod, abstractproperty
+import abc
import numbers
import numpy as np
@@ -16,7 +16,7 @@ from . import polyutils as pu
__all__ = ['ABCPolyBase']
-class ABCPolyBase(object):
+class ABCPolyBase(abc.ABC):
"""An abstract base class for immutable series classes.
ABCPolyBase provides the standard Python numerical methods
@@ -59,7 +59,6 @@ class ABCPolyBase(object):
Default window of the class.
"""
- __metaclass__ = ABCMeta
# Not hashable
__hash__ = None
@@ -70,68 +69,84 @@ class ABCPolyBase(object):
# Limit runaway size. T_n^m has degree n*m
maxpower = 100
- @abstractproperty
+ @property
+ @abc.abstractmethod
def domain(self):
pass
- @abstractproperty
+ @property
+ @abc.abstractmethod
def window(self):
pass
- @abstractproperty
+ @property
+ @abc.abstractmethod
def nickname(self):
pass
- @abstractproperty
+ @property
+ @abc.abstractmethod
def basis_name(self):
pass
- @abstractmethod
- def _add(self):
+ @staticmethod
+ @abc.abstractmethod
+ def _add(c1, c2):
pass
- @abstractmethod
- def _sub(self):
+ @staticmethod
+ @abc.abstractmethod
+ def _sub(c1, c2):
pass
- @abstractmethod
- def _mul(self):
+ @staticmethod
+ @abc.abstractmethod
+ def _mul(c1, c2):
pass
- @abstractmethod
- def _div(self):
+ @staticmethod
+ @abc.abstractmethod
+ def _div(c1, c2):
pass
- @abstractmethod
- def _pow(self):
+ @staticmethod
+ @abc.abstractmethod
+ def _pow(c, pow, maxpower=None):
pass
- @abstractmethod
- def _val(self):
+ @staticmethod
+ @abc.abstractmethod
+ def _val(x, c):
pass
- @abstractmethod
- def _int(self):
+ @staticmethod
+ @abc.abstractmethod
+ def _int(c, m, k, lbnd, scl):
pass
- @abstractmethod
- def _der(self):
+ @staticmethod
+ @abc.abstractmethod
+ def _der(c, m, scl):
pass
- @abstractmethod
- def _fit(self):
+ @staticmethod
+ @abc.abstractmethod
+ def _fit(x, y, deg, rcond, full):
pass
- @abstractmethod
- def _line(self):
+ @staticmethod
+ @abc.abstractmethod
+ def _line(off, scl):
pass
- @abstractmethod
- def _roots(self):
+ @staticmethod
+ @abc.abstractmethod
+ def _roots(c):
pass
- @abstractmethod
- def _fromroots(self):
+ @staticmethod
+ @abc.abstractmethod
+ def _fromroots(r):
pass
def has_samecoef(self, other):
diff --git a/numpy/polynomial/chebyshev.py b/numpy/polynomial/chebyshev.py
index e0734e1b8..c37b2c6d6 100644
--- a/numpy/polynomial/chebyshev.py
+++ b/numpy/polynomial/chebyshev.py
@@ -225,15 +225,15 @@ def _zseries_div(z1, z2):
"""
z1 = z1.copy()
z2 = z2.copy()
- len1 = len(z1)
- len2 = len(z2)
- if len2 == 1:
+ lc1 = len(z1)
+ lc2 = len(z2)
+ if lc2 == 1:
z1 /= z2
return z1, z1[:1]*0
- elif len1 < len2:
+ elif lc1 < lc2:
return z1[:1]*0, z1
else:
- dlen = len1 - len2
+ dlen = lc1 - lc2
scl = z2[0]
z2 /= scl
quo = np.empty(dlen + 1, dtype=z1.dtype)
@@ -244,16 +244,16 @@ def _zseries_div(z1, z2):
quo[i] = z1[i]
quo[dlen - i] = r
tmp = r*z2
- z1[i:i+len2] -= tmp
- z1[j:j+len2] -= tmp
+ z1[i:i+lc2] -= tmp
+ z1[j:j+lc2] -= tmp
i += 1
j -= 1
r = z1[i]
quo[i] = r
tmp = r*z2
- z1[i:i+len2] -= tmp
+ z1[i:i+lc2] -= tmp
quo /= scl
- rem = z1[i+1:i-1+len2].copy()
+ rem = z1[i+1:i-1+lc2].copy()
return quo, rem
@@ -541,21 +541,7 @@ def chebfromroots(roots):
array([1.5+0.j, 0. +0.j, 0.5+0.j])
"""
- if len(roots) == 0:
- return np.ones(1)
- else:
- [roots] = pu.as_series([roots], trim=False)
- roots.sort()
- p = [chebline(-r, 1) for r in roots]
- n = len(p)
- while n > 1:
- m, r = divmod(n, 2)
- tmp = [chebmul(p[i], p[i+m]) for i in range(m)]
- if r:
- tmp[0] = chebmul(tmp[0], p[-1])
- p = tmp
- n = m
- return p[0]
+ return pu._fromroots(chebline, chebmul, roots)
def chebadd(c1, c2):
@@ -597,15 +583,7 @@ def chebadd(c1, c2):
array([4., 4., 4.])
"""
- # c1, c2 are trimmed copies
- [c1, c2] = pu.as_series([c1, c2])
- if len(c1) > len(c2):
- c1[:c2.size] += c2
- ret = c1
- else:
- c2[:c1.size] += c1
- ret = c2
- return pu.trimseq(ret)
+ return pu._add(c1, c2)
def chebsub(c1, c2):
@@ -649,16 +627,7 @@ def chebsub(c1, c2):
array([ 2., 0., -2.])
"""
- # c1, c2 are trimmed copies
- [c1, c2] = pu.as_series([c1, c2])
- if len(c1) > len(c2):
- c1[:c2.size] -= c2
- ret = c1
- else:
- c2 = -c2
- c2[:c1.size] += c1
- ret = c2
- return pu.trimseq(ret)
+ return pu._sub(c1, c2)
def chebmulx(c):
@@ -807,6 +776,7 @@ def chebdiv(c1, c2):
if c2[-1] == 0:
raise ZeroDivisionError()
+ # note: this is more efficient than `pu._div(chebmul, c1, c2)`
lc1 = len(c1)
lc2 = len(c2)
if lc1 < lc2:
@@ -856,6 +826,9 @@ def chebpow(c, pow, maxpower=16):
array([15.5, 22. , 16. , ..., 12.5, 12. , 8. ])
"""
+ # note: this is more efficient than `pu._pow(chebmul, c1, c2)`, as it
+ # avoids converting between z and c series repeatedly
+
# c is a trimmed copy
[c] = pu.as_series([c])
power = int(pow)
@@ -940,14 +913,10 @@ def chebder(c, m=1, scl=1, axis=0):
c = np.array(c, ndmin=1, copy=1)
if c.dtype.char in '?bBhHiIlLqQpP':
c = c.astype(np.double)
- cnt, iaxis = [int(t) for t in [m, axis]]
-
- if cnt != m:
- raise ValueError("The order of derivation must be integer")
+ cnt = pu._deprecate_as_int(m, "the order of derivation")
+ iaxis = pu._deprecate_as_int(axis, "the axis")
if cnt < 0:
raise ValueError("The order of derivation must be non-negative")
- if iaxis != axis:
- raise ValueError("The axis must be integer")
iaxis = normalize_axis_index(iaxis, c.ndim)
if cnt == 0:
@@ -1063,10 +1032,8 @@ def chebint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
c = c.astype(np.double)
if not np.iterable(k):
k = [k]
- cnt, iaxis = [int(t) for t in [m, axis]]
-
- if cnt != m:
- raise ValueError("The order of integration must be integer")
+ cnt = pu._deprecate_as_int(m, "the order of integration")
+ iaxis = pu._deprecate_as_int(axis, "the axis")
if cnt < 0:
raise ValueError("The order of integration must be non-negative")
if len(k) > cnt:
@@ -1075,8 +1042,6 @@ def chebint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
raise ValueError("lbnd must be a scalar.")
if np.ndim(scl) != 0:
raise ValueError("scl must be a scalar.")
- if iaxis != axis:
- raise ValueError("The axis must be integer")
iaxis = normalize_axis_index(iaxis, c.ndim)
if cnt == 0:
@@ -1238,14 +1203,7 @@ def chebval2d(x, y, c):
.. versionadded:: 1.7.0
"""
- try:
- x, y = np.array((x, y), copy=0)
- except Exception:
- raise ValueError('x, y are incompatible')
-
- c = chebval(x, c)
- c = chebval(y, c, tensor=False)
- return c
+ return pu._valnd(chebval, c, x, y)
def chebgrid2d(x, y, c):
@@ -1298,9 +1256,7 @@ def chebgrid2d(x, y, c):
.. versionadded:: 1.7.0
"""
- c = chebval(x, c)
- c = chebval(y, c)
- return c
+ return pu._gridnd(chebval, c, x, y)
def chebval3d(x, y, z, c):
@@ -1351,15 +1307,7 @@ def chebval3d(x, y, z, c):
.. versionadded:: 1.7.0
"""
- try:
- x, y, z = np.array((x, y, z), copy=0)
- except Exception:
- raise ValueError('x, y, z are incompatible')
-
- c = chebval(x, c)
- c = chebval(y, c, tensor=False)
- c = chebval(z, c, tensor=False)
- return c
+ return pu._valnd(chebval, c, x, y, z)
def chebgrid3d(x, y, z, c):
@@ -1415,10 +1363,7 @@ def chebgrid3d(x, y, z, c):
.. versionadded:: 1.7.0
"""
- c = chebval(x, c)
- c = chebval(y, c)
- c = chebval(z, c)
- return c
+ return pu._gridnd(chebval, c, x, y, z)
def chebvander(x, deg):
@@ -1456,9 +1401,7 @@ def chebvander(x, deg):
the converted `x`.
"""
- ideg = int(deg)
- if ideg != deg:
- raise ValueError("deg must be integer")
+ ideg = pu._deprecate_as_int(deg, "deg")
if ideg < 0:
raise ValueError("deg must be non-negative")
@@ -1526,17 +1469,7 @@ def chebvander2d(x, y, deg):
.. versionadded:: 1.7.0
"""
- ideg = [int(d) for d in deg]
- is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
- if is_valid != [1, 1]:
- raise ValueError("degrees must be non-negative integers")
- degx, degy = ideg
- x, y = np.array((x, y), copy=0) + 0.0
-
- vx = chebvander(x, degx)
- vy = chebvander(y, degy)
- v = vx[..., None]*vy[..., None,:]
- return v.reshape(v.shape[:-2] + (-1,))
+ return pu._vander2d(chebvander, x, y, deg)
def chebvander3d(x, y, z, deg):
@@ -1590,18 +1523,7 @@ def chebvander3d(x, y, z, deg):
.. versionadded:: 1.7.0
"""
- ideg = [int(d) for d in deg]
- is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
- if is_valid != [1, 1, 1]:
- raise ValueError("degrees must be non-negative integers")
- degx, degy, degz = ideg
- x, y, z = np.array((x, y, z), copy=0) + 0.0
-
- vx = chebvander(x, degx)
- vy = chebvander(y, degy)
- vz = chebvander(z, degz)
- v = vx[..., None, None]*vy[..., None,:, None]*vz[..., None, None,:]
- return v.reshape(v.shape[:-3] + (-1,))
+ return pu._vander3d(chebvander, x, y, z, deg)
def chebfit(x, y, deg, rcond=None, full=False, w=None):
@@ -1723,81 +1645,7 @@ def chebfit(x, y, deg, rcond=None, full=False, w=None):
--------
"""
- x = np.asarray(x) + 0.0
- y = np.asarray(y) + 0.0
- deg = np.asarray(deg)
-
- # check arguments.
- if deg.ndim > 1 or deg.dtype.kind not in 'iu' or deg.size == 0:
- raise TypeError("deg must be an int or non-empty 1-D array of int")
- if deg.min() < 0:
- raise ValueError("expected deg >= 0")
- if x.ndim != 1:
- raise TypeError("expected 1D vector for x")
- if x.size == 0:
- raise TypeError("expected non-empty vector for x")
- if y.ndim < 1 or y.ndim > 2:
- raise TypeError("expected 1D or 2D array for y")
- if len(x) != len(y):
- raise TypeError("expected x and y to have same length")
-
- if deg.ndim == 0:
- lmax = deg
- order = lmax + 1
- van = chebvander(x, lmax)
- else:
- deg = np.sort(deg)
- lmax = deg[-1]
- order = len(deg)
- van = chebvander(x, lmax)[:, deg]
-
- # set up the least squares matrices in transposed form
- lhs = van.T
- rhs = y.T
- if w is not None:
- w = np.asarray(w) + 0.0
- if w.ndim != 1:
- raise TypeError("expected 1D vector for w")
- if len(x) != len(w):
- raise TypeError("expected x and w to have same length")
- # apply weights. Don't use inplace operations as they
- # can cause problems with NA.
- lhs = lhs * w
- rhs = rhs * w
-
- # set rcond
- if rcond is None:
- rcond = len(x)*np.finfo(x.dtype).eps
-
- # Determine the norms of the design matrix columns.
- if issubclass(lhs.dtype.type, np.complexfloating):
- scl = np.sqrt((np.square(lhs.real) + np.square(lhs.imag)).sum(1))
- else:
- scl = np.sqrt(np.square(lhs).sum(1))
- scl[scl == 0] = 1
-
- # Solve the least squares problem.
- c, resids, rank, s = la.lstsq(lhs.T/scl, rhs.T, rcond)
- c = (c.T/scl).T
-
- # Expand c to include non-fitted coefficients which are set to zero
- if deg.ndim > 0:
- if c.ndim == 2:
- cc = np.zeros((lmax + 1, c.shape[1]), dtype=c.dtype)
- else:
- cc = np.zeros(lmax + 1, dtype=c.dtype)
- cc[deg] = c
- c = cc
-
- # warn on rank reduction
- if rank != order and not full:
- msg = "The fit may be poorly conditioned"
- warnings.warn(msg, pu.RankWarning, stacklevel=2)
-
- if full:
- return c, [resids, rank, s, rcond]
- else:
- return c
+ return pu._fit(chebvander, x, y, deg, rcond, full, w)
def chebcompanion(c):
@@ -1895,7 +1743,8 @@ def chebroots(c):
if len(c) == 2:
return np.array([-c[0]/c[1]])
- m = chebcompanion(c)
+ # rotated companion matrix reduces error
+ m = chebcompanion(c)[::-1,::-1]
r = la.eigvals(m)
r.sort()
return r
@@ -2003,9 +1852,9 @@ def chebgauss(deg):
.. math:: w_i = \\pi / n
"""
- ideg = int(deg)
- if ideg != deg or ideg < 1:
- raise ValueError("deg must be a non-negative integer")
+ ideg = pu._deprecate_as_int(deg, "deg")
+ if ideg <= 0:
+ raise ValueError("deg must be a positive integer")
x = np.cos(np.pi * np.arange(1, 2*ideg, 2) / (2.0*ideg))
w = np.ones(ideg)*(np.pi/ideg)
diff --git a/numpy/polynomial/hermite.py b/numpy/polynomial/hermite.py
index 93c9fc564..3870d1054 100644
--- a/numpy/polynomial/hermite.py
+++ b/numpy/polynomial/hermite.py
@@ -286,21 +286,7 @@ def hermfromroots(roots):
array([0.+0.j, 0.+0.j])
"""
- if len(roots) == 0:
- return np.ones(1)
- else:
- [roots] = pu.as_series([roots], trim=False)
- roots.sort()
- p = [hermline(-r, 1) for r in roots]
- n = len(p)
- while n > 1:
- m, r = divmod(n, 2)
- tmp = [hermmul(p[i], p[i+m]) for i in range(m)]
- if r:
- tmp[0] = hermmul(tmp[0], p[-1])
- p = tmp
- n = m
- return p[0]
+ return pu._fromroots(hermline, hermmul, roots)
def hermadd(c1, c2):
@@ -340,15 +326,7 @@ def hermadd(c1, c2):
array([2., 4., 6., 4.])
"""
- # c1, c2 are trimmed copies
- [c1, c2] = pu.as_series([c1, c2])
- if len(c1) > len(c2):
- c1[:c2.size] += c2
- ret = c1
- else:
- c2[:c1.size] += c1
- ret = c2
- return pu.trimseq(ret)
+ return pu._add(c1, c2)
def hermsub(c1, c2):
@@ -388,16 +366,7 @@ def hermsub(c1, c2):
array([0., 0., 0., 4.])
"""
- # c1, c2 are trimmed copies
- [c1, c2] = pu.as_series([c1, c2])
- if len(c1) > len(c2):
- c1[:c2.size] -= c2
- ret = c1
- else:
- c2 = -c2
- c2[:c1.size] += c1
- ret = c2
- return pu.trimseq(ret)
+ return pu._sub(c1, c2)
def hermmulx(c):
@@ -564,26 +533,7 @@ def hermdiv(c1, c2):
(array([1., 2., 3.]), array([1., 1.]))
"""
- # c1, c2 are trimmed copies
- [c1, c2] = pu.as_series([c1, c2])
- if c2[-1] == 0:
- raise ZeroDivisionError()
-
- lc1 = len(c1)
- lc2 = len(c2)
- if lc1 < lc2:
- return c1[:1]*0, c1
- elif lc2 == 1:
- return c1/c2[-1], c1[:1]*0
- else:
- quo = np.empty(lc1 - lc2 + 1, dtype=c1.dtype)
- rem = c1
- for i in range(lc1 - lc2, - 1, -1):
- p = hermmul([0]*i + [1], c2)
- q = rem[-1]/p[-1]
- rem = rem[:-1] - q*p[:-1]
- quo[i] = q
- return quo, pu.trimseq(rem)
+ return pu._div(hermmul, c1, c2)
def hermpow(c, pow, maxpower=16):
@@ -620,24 +570,7 @@ def hermpow(c, pow, maxpower=16):
array([81., 52., 82., 12., 9.])
"""
- # c is a trimmed copy
- [c] = pu.as_series([c])
- power = int(pow)
- if power != pow or power < 0:
- raise ValueError("Power must be a non-negative integer.")
- elif maxpower is not None and power > maxpower:
- raise ValueError("Power is too large")
- elif power == 0:
- return np.array([1], dtype=c.dtype)
- elif power == 1:
- return c
- else:
- # This can be made more efficient by using powers of two
- # in the usual way.
- prd = c
- for i in range(2, power + 1):
- prd = hermmul(prd, c)
- return prd
+ return pu._pow(hermmul, c, pow, maxpower)
def hermder(c, m=1, scl=1, axis=0):
@@ -698,14 +631,10 @@ def hermder(c, m=1, scl=1, axis=0):
c = np.array(c, ndmin=1, copy=1)
if c.dtype.char in '?bBhHiIlLqQpP':
c = c.astype(np.double)
- cnt, iaxis = [int(t) for t in [m, axis]]
-
- if cnt != m:
- raise ValueError("The order of derivation must be integer")
+ cnt = pu._deprecate_as_int(m, "the order of derivation")
+ iaxis = pu._deprecate_as_int(axis, "the axis")
if cnt < 0:
raise ValueError("The order of derivation must be non-negative")
- if iaxis != axis:
- raise ValueError("The axis must be integer")
iaxis = normalize_axis_index(iaxis, c.ndim)
if cnt == 0:
@@ -815,10 +744,8 @@ def hermint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
c = c.astype(np.double)
if not np.iterable(k):
k = [k]
- cnt, iaxis = [int(t) for t in [m, axis]]
-
- if cnt != m:
- raise ValueError("The order of integration must be integer")
+ cnt = pu._deprecate_as_int(m, "the order of integration")
+ iaxis = pu._deprecate_as_int(axis, "the axis")
if cnt < 0:
raise ValueError("The order of integration must be non-negative")
if len(k) > cnt:
@@ -827,8 +754,6 @@ def hermint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
raise ValueError("lbnd must be a scalar.")
if np.ndim(scl) != 0:
raise ValueError("scl must be a scalar.")
- if iaxis != axis:
- raise ValueError("The axis must be integer")
iaxis = normalize_axis_index(iaxis, c.ndim)
if cnt == 0:
@@ -995,14 +920,7 @@ def hermval2d(x, y, c):
.. versionadded:: 1.7.0
"""
- try:
- x, y = np.array((x, y), copy=0)
- except Exception:
- raise ValueError('x, y are incompatible')
-
- c = hermval(x, c)
- c = hermval(y, c, tensor=False)
- return c
+ return pu._valnd(hermval, c, x, y)
def hermgrid2d(x, y, c):
@@ -1055,9 +973,7 @@ def hermgrid2d(x, y, c):
.. versionadded:: 1.7.0
"""
- c = hermval(x, c)
- c = hermval(y, c)
- return c
+ return pu._gridnd(hermval, c, x, y)
def hermval3d(x, y, z, c):
@@ -1108,15 +1024,7 @@ def hermval3d(x, y, z, c):
.. versionadded:: 1.7.0
"""
- try:
- x, y, z = np.array((x, y, z), copy=0)
- except Exception:
- raise ValueError('x, y, z are incompatible')
-
- c = hermval(x, c)
- c = hermval(y, c, tensor=False)
- c = hermval(z, c, tensor=False)
- return c
+ return pu._valnd(hermval, c, x, y, z)
def hermgrid3d(x, y, z, c):
@@ -1172,10 +1080,7 @@ def hermgrid3d(x, y, z, c):
.. versionadded:: 1.7.0
"""
- c = hermval(x, c)
- c = hermval(y, c)
- c = hermval(z, c)
- return c
+ return pu._gridnd(hermval, c, x, y, z)
def hermvander(x, deg):
@@ -1222,9 +1127,7 @@ def hermvander(x, deg):
[ 1., 2., 2., -4.]])
"""
- ideg = int(deg)
- if ideg != deg:
- raise ValueError("deg must be integer")
+ ideg = pu._deprecate_as_int(deg, "deg")
if ideg < 0:
raise ValueError("deg must be non-negative")
@@ -1291,17 +1194,7 @@ def hermvander2d(x, y, deg):
.. versionadded:: 1.7.0
"""
- ideg = [int(d) for d in deg]
- is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
- if is_valid != [1, 1]:
- raise ValueError("degrees must be non-negative integers")
- degx, degy = ideg
- x, y = np.array((x, y), copy=0) + 0.0
-
- vx = hermvander(x, degx)
- vy = hermvander(y, degy)
- v = vx[..., None]*vy[..., None,:]
- return v.reshape(v.shape[:-2] + (-1,))
+ return pu._vander2d(hermvander, x, y, deg)
def hermvander3d(x, y, z, deg):
@@ -1355,18 +1248,7 @@ def hermvander3d(x, y, z, deg):
.. versionadded:: 1.7.0
"""
- ideg = [int(d) for d in deg]
- is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
- if is_valid != [1, 1, 1]:
- raise ValueError("degrees must be non-negative integers")
- degx, degy, degz = ideg
- x, y, z = np.array((x, y, z), copy=0) + 0.0
-
- vx = hermvander(x, degx)
- vy = hermvander(y, degy)
- vz = hermvander(z, degz)
- v = vx[..., None, None]*vy[..., None,:, None]*vz[..., None, None,:]
- return v.reshape(v.shape[:-3] + (-1,))
+ return pu._vander3d(hermvander, x, y, z, deg)
def hermfit(x, y, deg, rcond=None, full=False, w=None):
@@ -1493,81 +1375,7 @@ def hermfit(x, y, deg, rcond=None, full=False, w=None):
array([1.0218, 1.9986, 2.9999]) # may vary
"""
- x = np.asarray(x) + 0.0
- y = np.asarray(y) + 0.0
- deg = np.asarray(deg)
-
- # check arguments.
- if deg.ndim > 1 or deg.dtype.kind not in 'iu' or deg.size == 0:
- raise TypeError("deg must be an int or non-empty 1-D array of int")
- if deg.min() < 0:
- raise ValueError("expected deg >= 0")
- if x.ndim != 1:
- raise TypeError("expected 1D vector for x")
- if x.size == 0:
- raise TypeError("expected non-empty vector for x")
- if y.ndim < 1 or y.ndim > 2:
- raise TypeError("expected 1D or 2D array for y")
- if len(x) != len(y):
- raise TypeError("expected x and y to have same length")
-
- if deg.ndim == 0:
- lmax = deg
- order = lmax + 1
- van = hermvander(x, lmax)
- else:
- deg = np.sort(deg)
- lmax = deg[-1]
- order = len(deg)
- van = hermvander(x, lmax)[:, deg]
-
- # set up the least squares matrices in transposed form
- lhs = van.T
- rhs = y.T
- if w is not None:
- w = np.asarray(w) + 0.0
- if w.ndim != 1:
- raise TypeError("expected 1D vector for w")
- if len(x) != len(w):
- raise TypeError("expected x and w to have same length")
- # apply weights. Don't use inplace operations as they
- # can cause problems with NA.
- lhs = lhs * w
- rhs = rhs * w
-
- # set rcond
- if rcond is None:
- rcond = len(x)*np.finfo(x.dtype).eps
-
- # Determine the norms of the design matrix columns.
- if issubclass(lhs.dtype.type, np.complexfloating):
- scl = np.sqrt((np.square(lhs.real) + np.square(lhs.imag)).sum(1))
- else:
- scl = np.sqrt(np.square(lhs).sum(1))
- scl[scl == 0] = 1
-
- # Solve the least squares problem.
- c, resids, rank, s = la.lstsq(lhs.T/scl, rhs.T, rcond)
- c = (c.T/scl).T
-
- # Expand c to include non-fitted coefficients which are set to zero
- if deg.ndim > 0:
- if c.ndim == 2:
- cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype)
- else:
- cc = np.zeros(lmax+1, dtype=c.dtype)
- cc[deg] = c
- c = cc
-
- # warn on rank reduction
- if rank != order and not full:
- msg = "The fit may be poorly conditioned"
- warnings.warn(msg, pu.RankWarning, stacklevel=2)
-
- if full:
- return c, [resids, rank, s, rcond]
- else:
- return c
+ return pu._fit(hermvander, x, y, deg, rcond, full, w)
def hermcompanion(c):
@@ -1668,7 +1476,8 @@ def hermroots(c):
if len(c) == 2:
return np.array([-.5*c[0]/c[1]])
- m = hermcompanion(c)
+ # rotated companion matrix reduces error
+ m = hermcompanion(c)[::-1,::-1]
r = la.eigvals(m)
r.sort()
return r
@@ -1753,9 +1562,9 @@ def hermgauss(deg):
the right value when integrating 1.
"""
- ideg = int(deg)
- if ideg != deg or ideg < 1:
- raise ValueError("deg must be a non-negative integer")
+ ideg = pu._deprecate_as_int(deg, "deg")
+ if ideg <= 0:
+ raise ValueError("deg must be a positive integer")
# first approximation of roots. We use the fact that the companion
# matrix is symmetric in this case in order to obtain better zeros.
diff --git a/numpy/polynomial/hermite_e.py b/numpy/polynomial/hermite_e.py
index bafb4b997..b12c0d792 100644
--- a/numpy/polynomial/hermite_e.py
+++ b/numpy/polynomial/hermite_e.py
@@ -287,21 +287,7 @@ def hermefromroots(roots):
array([0.+0.j, 0.+0.j])
"""
- if len(roots) == 0:
- return np.ones(1)
- else:
- [roots] = pu.as_series([roots], trim=False)
- roots.sort()
- p = [hermeline(-r, 1) for r in roots]
- n = len(p)
- while n > 1:
- m, r = divmod(n, 2)
- tmp = [hermemul(p[i], p[i+m]) for i in range(m)]
- if r:
- tmp[0] = hermemul(tmp[0], p[-1])
- p = tmp
- n = m
- return p[0]
+ return pu._fromroots(hermeline, hermemul, roots)
def hermeadd(c1, c2):
@@ -341,15 +327,7 @@ def hermeadd(c1, c2):
array([2., 4., 6., 4.])
"""
- # c1, c2 are trimmed copies
- [c1, c2] = pu.as_series([c1, c2])
- if len(c1) > len(c2):
- c1[:c2.size] += c2
- ret = c1
- else:
- c2[:c1.size] += c1
- ret = c2
- return pu.trimseq(ret)
+ return pu._add(c1, c2)
def hermesub(c1, c2):
@@ -389,16 +367,7 @@ def hermesub(c1, c2):
array([0., 0., 0., 4.])
"""
- # c1, c2 are trimmed copies
- [c1, c2] = pu.as_series([c1, c2])
- if len(c1) > len(c2):
- c1[:c2.size] -= c2
- ret = c1
- else:
- c2 = -c2
- c2[:c1.size] += c1
- ret = c2
- return pu.trimseq(ret)
+ return pu._sub(c1, c2)
def hermemulx(c):
@@ -559,26 +528,7 @@ def hermediv(c1, c2):
(array([1., 2., 3.]), array([1., 2.]))
"""
- # c1, c2 are trimmed copies
- [c1, c2] = pu.as_series([c1, c2])
- if c2[-1] == 0:
- raise ZeroDivisionError()
-
- lc1 = len(c1)
- lc2 = len(c2)
- if lc1 < lc2:
- return c1[:1]*0, c1
- elif lc2 == 1:
- return c1/c2[-1], c1[:1]*0
- else:
- quo = np.empty(lc1 - lc2 + 1, dtype=c1.dtype)
- rem = c1
- for i in range(lc1 - lc2, - 1, -1):
- p = hermemul([0]*i + [1], c2)
- q = rem[-1]/p[-1]
- rem = rem[:-1] - q*p[:-1]
- quo[i] = q
- return quo, pu.trimseq(rem)
+ return pu._div(hermemul, c1, c2)
def hermepow(c, pow, maxpower=16):
@@ -615,24 +565,7 @@ def hermepow(c, pow, maxpower=16):
array([23., 28., 46., 12., 9.])
"""
- # c is a trimmed copy
- [c] = pu.as_series([c])
- power = int(pow)
- if power != pow or power < 0:
- raise ValueError("Power must be a non-negative integer.")
- elif maxpower is not None and power > maxpower:
- raise ValueError("Power is too large")
- elif power == 0:
- return np.array([1], dtype=c.dtype)
- elif power == 1:
- return c
- else:
- # This can be made more efficient by using powers of two
- # in the usual way.
- prd = c
- for i in range(2, power + 1):
- prd = hermemul(prd, c)
- return prd
+ return pu._pow(hermemul, c, pow, maxpower)
def hermeder(c, m=1, scl=1, axis=0):
@@ -693,14 +626,10 @@ def hermeder(c, m=1, scl=1, axis=0):
c = np.array(c, ndmin=1, copy=1)
if c.dtype.char in '?bBhHiIlLqQpP':
c = c.astype(np.double)
- cnt, iaxis = [int(t) for t in [m, axis]]
-
- if cnt != m:
- raise ValueError("The order of derivation must be integer")
+ cnt = pu._deprecate_as_int(m, "the order of derivation")
+ iaxis = pu._deprecate_as_int(axis, "the axis")
if cnt < 0:
raise ValueError("The order of derivation must be non-negative")
- if iaxis != axis:
- raise ValueError("The axis must be integer")
iaxis = normalize_axis_index(iaxis, c.ndim)
if cnt == 0:
@@ -810,10 +739,8 @@ def hermeint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
c = c.astype(np.double)
if not np.iterable(k):
k = [k]
- cnt, iaxis = [int(t) for t in [m, axis]]
-
- if cnt != m:
- raise ValueError("The order of integration must be integer")
+ cnt = pu._deprecate_as_int(m, "the order of integration")
+ iaxis = pu._deprecate_as_int(axis, "the axis")
if cnt < 0:
raise ValueError("The order of integration must be non-negative")
if len(k) > cnt:
@@ -822,8 +749,6 @@ def hermeint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
raise ValueError("lbnd must be a scalar.")
if np.ndim(scl) != 0:
raise ValueError("scl must be a scalar.")
- if iaxis != axis:
- raise ValueError("The axis must be integer")
iaxis = normalize_axis_index(iaxis, c.ndim)
if cnt == 0:
@@ -989,14 +914,7 @@ def hermeval2d(x, y, c):
.. versionadded:: 1.7.0
"""
- try:
- x, y = np.array((x, y), copy=0)
- except Exception:
- raise ValueError('x, y are incompatible')
-
- c = hermeval(x, c)
- c = hermeval(y, c, tensor=False)
- return c
+ return pu._valnd(hermeval, c, x, y)
def hermegrid2d(x, y, c):
@@ -1049,9 +967,7 @@ def hermegrid2d(x, y, c):
.. versionadded:: 1.7.0
"""
- c = hermeval(x, c)
- c = hermeval(y, c)
- return c
+ return pu._gridnd(hermeval, c, x, y)
def hermeval3d(x, y, z, c):
@@ -1102,15 +1018,7 @@ def hermeval3d(x, y, z, c):
.. versionadded:: 1.7.0
"""
- try:
- x, y, z = np.array((x, y, z), copy=0)
- except Exception:
- raise ValueError('x, y, z are incompatible')
-
- c = hermeval(x, c)
- c = hermeval(y, c, tensor=False)
- c = hermeval(z, c, tensor=False)
- return c
+ return pu._valnd(hermeval, c, x, y, z)
def hermegrid3d(x, y, z, c):
@@ -1166,10 +1074,7 @@ def hermegrid3d(x, y, z, c):
.. versionadded:: 1.7.0
"""
- c = hermeval(x, c)
- c = hermeval(y, c)
- c = hermeval(z, c)
- return c
+ return pu._gridnd(hermeval, c, x, y, z)
def hermevander(x, deg):
@@ -1216,9 +1121,7 @@ def hermevander(x, deg):
[ 1., 1., 0., -2.]])
"""
- ideg = int(deg)
- if ideg != deg:
- raise ValueError("deg must be integer")
+ ideg = pu._deprecate_as_int(deg, "deg")
if ideg < 0:
raise ValueError("deg must be non-negative")
@@ -1284,17 +1187,7 @@ def hermevander2d(x, y, deg):
.. versionadded:: 1.7.0
"""
- ideg = [int(d) for d in deg]
- is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
- if is_valid != [1, 1]:
- raise ValueError("degrees must be non-negative integers")
- degx, degy = ideg
- x, y = np.array((x, y), copy=0) + 0.0
-
- vx = hermevander(x, degx)
- vy = hermevander(y, degy)
- v = vx[..., None]*vy[..., None,:]
- return v.reshape(v.shape[:-2] + (-1,))
+ return pu._vander2d(hermevander, x, y, deg)
def hermevander3d(x, y, z, deg):
@@ -1348,18 +1241,7 @@ def hermevander3d(x, y, z, deg):
.. versionadded:: 1.7.0
"""
- ideg = [int(d) for d in deg]
- is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
- if is_valid != [1, 1, 1]:
- raise ValueError("degrees must be non-negative integers")
- degx, degy, degz = ideg
- x, y, z = np.array((x, y, z), copy=0) + 0.0
-
- vx = hermevander(x, degx)
- vy = hermevander(y, degy)
- vz = hermevander(z, degz)
- v = vx[..., None, None]*vy[..., None,:, None]*vz[..., None, None,:]
- return v.reshape(v.shape[:-3] + (-1,))
+ return pu._vander3d(hermevander, x, y, z, deg)
def hermefit(x, y, deg, rcond=None, full=False, w=None):
@@ -1487,81 +1369,7 @@ def hermefit(x, y, deg, rcond=None, full=False, w=None):
array([ 1.01690445, 1.99951418, 2.99948696]) # may vary
"""
- x = np.asarray(x) + 0.0
- y = np.asarray(y) + 0.0
- deg = np.asarray(deg)
-
- # check arguments.
- if deg.ndim > 1 or deg.dtype.kind not in 'iu' or deg.size == 0:
- raise TypeError("deg must be an int or non-empty 1-D array of int")
- if deg.min() < 0:
- raise ValueError("expected deg >= 0")
- if x.ndim != 1:
- raise TypeError("expected 1D vector for x")
- if x.size == 0:
- raise TypeError("expected non-empty vector for x")
- if y.ndim < 1 or y.ndim > 2:
- raise TypeError("expected 1D or 2D array for y")
- if len(x) != len(y):
- raise TypeError("expected x and y to have same length")
-
- if deg.ndim == 0:
- lmax = deg
- order = lmax + 1
- van = hermevander(x, lmax)
- else:
- deg = np.sort(deg)
- lmax = deg[-1]
- order = len(deg)
- van = hermevander(x, lmax)[:, deg]
-
- # set up the least squares matrices in transposed form
- lhs = van.T
- rhs = y.T
- if w is not None:
- w = np.asarray(w) + 0.0
- if w.ndim != 1:
- raise TypeError("expected 1D vector for w")
- if len(x) != len(w):
- raise TypeError("expected x and w to have same length")
- # apply weights. Don't use inplace operations as they
- # can cause problems with NA.
- lhs = lhs * w
- rhs = rhs * w
-
- # set rcond
- if rcond is None:
- rcond = len(x)*np.finfo(x.dtype).eps
-
- # Determine the norms of the design matrix columns.
- if issubclass(lhs.dtype.type, np.complexfloating):
- scl = np.sqrt((np.square(lhs.real) + np.square(lhs.imag)).sum(1))
- else:
- scl = np.sqrt(np.square(lhs).sum(1))
- scl[scl == 0] = 1
-
- # Solve the least squares problem.
- c, resids, rank, s = la.lstsq(lhs.T/scl, rhs.T, rcond)
- c = (c.T/scl).T
-
- # Expand c to include non-fitted coefficients which are set to zero
- if deg.ndim > 0:
- if c.ndim == 2:
- cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype)
- else:
- cc = np.zeros(lmax+1, dtype=c.dtype)
- cc[deg] = c
- c = cc
-
- # warn on rank reduction
- if rank != order and not full:
- msg = "The fit may be poorly conditioned"
- warnings.warn(msg, pu.RankWarning, stacklevel=2)
-
- if full:
- return c, [resids, rank, s, rcond]
- else:
- return c
+ return pu._fit(hermevander, x, y, deg, rcond, full, w)
def hermecompanion(c):
@@ -1663,7 +1471,8 @@ def hermeroots(c):
if len(c) == 2:
return np.array([-c[0]/c[1]])
- m = hermecompanion(c)
+ # rotated companion matrix reduces error
+ m = hermecompanion(c)[::-1,::-1]
r = la.eigvals(m)
r.sort()
return r
@@ -1748,9 +1557,9 @@ def hermegauss(deg):
the right value when integrating 1.
"""
- ideg = int(deg)
- if ideg != deg or ideg < 1:
- raise ValueError("deg must be a non-negative integer")
+ ideg = pu._deprecate_as_int(deg, "deg")
+ if ideg <= 0:
+ raise ValueError("deg must be a positive integer")
# first approximation of roots. We use the fact that the companion
# matrix is symmetric in this case in order to obtain better zeros.
diff --git a/numpy/polynomial/laguerre.py b/numpy/polynomial/laguerre.py
index 9207c9afe..e103dde92 100644
--- a/numpy/polynomial/laguerre.py
+++ b/numpy/polynomial/laguerre.py
@@ -283,21 +283,7 @@ def lagfromroots(roots):
array([0.+0.j, 0.+0.j])
"""
- if len(roots) == 0:
- return np.ones(1)
- else:
- [roots] = pu.as_series([roots], trim=False)
- roots.sort()
- p = [lagline(-r, 1) for r in roots]
- n = len(p)
- while n > 1:
- m, r = divmod(n, 2)
- tmp = [lagmul(p[i], p[i+m]) for i in range(m)]
- if r:
- tmp[0] = lagmul(tmp[0], p[-1])
- p = tmp
- n = m
- return p[0]
+ return pu._fromroots(lagline, lagmul, roots)
def lagadd(c1, c2):
@@ -338,15 +324,7 @@ def lagadd(c1, c2):
"""
- # c1, c2 are trimmed copies
- [c1, c2] = pu.as_series([c1, c2])
- if len(c1) > len(c2):
- c1[:c2.size] += c2
- ret = c1
- else:
- c2[:c1.size] += c1
- ret = c2
- return pu.trimseq(ret)
+ return pu._add(c1, c2)
def lagsub(c1, c2):
@@ -386,16 +364,7 @@ def lagsub(c1, c2):
array([0., 0., 0., 4.])
"""
- # c1, c2 are trimmed copies
- [c1, c2] = pu.as_series([c1, c2])
- if len(c1) > len(c2):
- c1[:c2.size] -= c2
- ret = c1
- else:
- c2 = -c2
- c2[:c1.size] += c1
- ret = c2
- return pu.trimseq(ret)
+ return pu._sub(c1, c2)
def lagmulx(c):
@@ -561,26 +530,7 @@ def lagdiv(c1, c2):
(array([1., 2., 3.]), array([1., 1.]))
"""
- # c1, c2 are trimmed copies
- [c1, c2] = pu.as_series([c1, c2])
- if c2[-1] == 0:
- raise ZeroDivisionError()
-
- lc1 = len(c1)
- lc2 = len(c2)
- if lc1 < lc2:
- return c1[:1]*0, c1
- elif lc2 == 1:
- return c1/c2[-1], c1[:1]*0
- else:
- quo = np.empty(lc1 - lc2 + 1, dtype=c1.dtype)
- rem = c1
- for i in range(lc1 - lc2, - 1, -1):
- p = lagmul([0]*i + [1], c2)
- q = rem[-1]/p[-1]
- rem = rem[:-1] - q*p[:-1]
- quo[i] = q
- return quo, pu.trimseq(rem)
+ return pu._div(lagmul, c1, c2)
def lagpow(c, pow, maxpower=16):
@@ -617,24 +567,7 @@ def lagpow(c, pow, maxpower=16):
array([ 14., -16., 56., -72., 54.])
"""
- # c is a trimmed copy
- [c] = pu.as_series([c])
- power = int(pow)
- if power != pow or power < 0:
- raise ValueError("Power must be a non-negative integer.")
- elif maxpower is not None and power > maxpower:
- raise ValueError("Power is too large")
- elif power == 0:
- return np.array([1], dtype=c.dtype)
- elif power == 1:
- return c
- else:
- # This can be made more efficient by using powers of two
- # in the usual way.
- prd = c
- for i in range(2, power + 1):
- prd = lagmul(prd, c)
- return prd
+ return pu._pow(lagmul, c, pow, maxpower)
def lagder(c, m=1, scl=1, axis=0):
@@ -695,14 +628,11 @@ def lagder(c, m=1, scl=1, axis=0):
c = np.array(c, ndmin=1, copy=1)
if c.dtype.char in '?bBhHiIlLqQpP':
c = c.astype(np.double)
- cnt, iaxis = [int(t) for t in [m, axis]]
- if cnt != m:
- raise ValueError("The order of derivation must be integer")
+ cnt = pu._deprecate_as_int(m, "the order of derivation")
+ iaxis = pu._deprecate_as_int(axis, "the axis")
if cnt < 0:
raise ValueError("The order of derivation must be non-negative")
- if iaxis != axis:
- raise ValueError("The axis must be integer")
iaxis = normalize_axis_index(iaxis, c.ndim)
if cnt == 0:
@@ -815,10 +745,8 @@ def lagint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
c = c.astype(np.double)
if not np.iterable(k):
k = [k]
- cnt, iaxis = [int(t) for t in [m, axis]]
-
- if cnt != m:
- raise ValueError("The order of integration must be integer")
+ cnt = pu._deprecate_as_int(m, "the order of integration")
+ iaxis = pu._deprecate_as_int(axis, "the axis")
if cnt < 0:
raise ValueError("The order of integration must be non-negative")
if len(k) > cnt:
@@ -827,8 +755,6 @@ def lagint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
raise ValueError("lbnd must be a scalar.")
if np.ndim(scl) != 0:
raise ValueError("scl must be a scalar.")
- if iaxis != axis:
- raise ValueError("The axis must be integer")
iaxis = normalize_axis_index(iaxis, c.ndim)
if cnt == 0:
@@ -995,14 +921,7 @@ def lagval2d(x, y, c):
.. versionadded:: 1.7.0
"""
- try:
- x, y = np.array((x, y), copy=0)
- except Exception:
- raise ValueError('x, y are incompatible')
-
- c = lagval(x, c)
- c = lagval(y, c, tensor=False)
- return c
+ return pu._valnd(lagval, c, x, y)
def laggrid2d(x, y, c):
@@ -1055,9 +974,7 @@ def laggrid2d(x, y, c):
.. versionadded:: 1.7.0
"""
- c = lagval(x, c)
- c = lagval(y, c)
- return c
+ return pu._gridnd(lagval, c, x, y)
def lagval3d(x, y, z, c):
@@ -1108,15 +1025,7 @@ def lagval3d(x, y, z, c):
.. versionadded:: 1.7.0
"""
- try:
- x, y, z = np.array((x, y, z), copy=0)
- except Exception:
- raise ValueError('x, y, z are incompatible')
-
- c = lagval(x, c)
- c = lagval(y, c, tensor=False)
- c = lagval(z, c, tensor=False)
- return c
+ return pu._valnd(lagval, c, x, y, z)
def laggrid3d(x, y, z, c):
@@ -1172,10 +1081,7 @@ def laggrid3d(x, y, z, c):
.. versionadded:: 1.7.0
"""
- c = lagval(x, c)
- c = lagval(y, c)
- c = lagval(z, c)
- return c
+ return pu._gridnd(lagval, c, x, y, z)
def lagvander(x, deg):
@@ -1222,9 +1128,7 @@ def lagvander(x, deg):
[ 1. , -1. , -1. , -0.33333333]])
"""
- ideg = int(deg)
- if ideg != deg:
- raise ValueError("deg must be integer")
+ ideg = pu._deprecate_as_int(deg, "deg")
if ideg < 0:
raise ValueError("deg must be non-negative")
@@ -1290,17 +1194,7 @@ def lagvander2d(x, y, deg):
.. versionadded:: 1.7.0
"""
- ideg = [int(d) for d in deg]
- is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
- if is_valid != [1, 1]:
- raise ValueError("degrees must be non-negative integers")
- degx, degy = ideg
- x, y = np.array((x, y), copy=0) + 0.0
-
- vx = lagvander(x, degx)
- vy = lagvander(y, degy)
- v = vx[..., None]*vy[..., None,:]
- return v.reshape(v.shape[:-2] + (-1,))
+ return pu._vander2d(lagvander, x, y, deg)
def lagvander3d(x, y, z, deg):
@@ -1354,18 +1248,7 @@ def lagvander3d(x, y, z, deg):
.. versionadded:: 1.7.0
"""
- ideg = [int(d) for d in deg]
- is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
- if is_valid != [1, 1, 1]:
- raise ValueError("degrees must be non-negative integers")
- degx, degy, degz = ideg
- x, y, z = np.array((x, y, z), copy=0) + 0.0
-
- vx = lagvander(x, degx)
- vy = lagvander(y, degy)
- vz = lagvander(z, degz)
- v = vx[..., None, None]*vy[..., None,:, None]*vz[..., None, None,:]
- return v.reshape(v.shape[:-3] + (-1,))
+ return pu._vander3d(lagvander, x, y, z, deg)
def lagfit(x, y, deg, rcond=None, full=False, w=None):
@@ -1492,81 +1375,7 @@ def lagfit(x, y, deg, rcond=None, full=False, w=None):
array([ 0.96971004, 2.00193749, 3.00288744]) # may vary
"""
- x = np.asarray(x) + 0.0
- y = np.asarray(y) + 0.0
- deg = np.asarray(deg)
-
- # check arguments.
- if deg.ndim > 1 or deg.dtype.kind not in 'iu' or deg.size == 0:
- raise TypeError("deg must be an int or non-empty 1-D array of int")
- if deg.min() < 0:
- raise ValueError("expected deg >= 0")
- if x.ndim != 1:
- raise TypeError("expected 1D vector for x")
- if x.size == 0:
- raise TypeError("expected non-empty vector for x")
- if y.ndim < 1 or y.ndim > 2:
- raise TypeError("expected 1D or 2D array for y")
- if len(x) != len(y):
- raise TypeError("expected x and y to have same length")
-
- if deg.ndim == 0:
- lmax = deg
- order = lmax + 1
- van = lagvander(x, lmax)
- else:
- deg = np.sort(deg)
- lmax = deg[-1]
- order = len(deg)
- van = lagvander(x, lmax)[:, deg]
-
- # set up the least squares matrices in transposed form
- lhs = van.T
- rhs = y.T
- if w is not None:
- w = np.asarray(w) + 0.0
- if w.ndim != 1:
- raise TypeError("expected 1D vector for w")
- if len(x) != len(w):
- raise TypeError("expected x and w to have same length")
- # apply weights. Don't use inplace operations as they
- # can cause problems with NA.
- lhs = lhs * w
- rhs = rhs * w
-
- # set rcond
- if rcond is None:
- rcond = len(x)*np.finfo(x.dtype).eps
-
- # Determine the norms of the design matrix columns.
- if issubclass(lhs.dtype.type, np.complexfloating):
- scl = np.sqrt((np.square(lhs.real) + np.square(lhs.imag)).sum(1))
- else:
- scl = np.sqrt(np.square(lhs).sum(1))
- scl[scl == 0] = 1
-
- # Solve the least squares problem.
- c, resids, rank, s = la.lstsq(lhs.T/scl, rhs.T, rcond)
- c = (c.T/scl).T
-
- # Expand c to include non-fitted coefficients which are set to zero
- if deg.ndim > 0:
- if c.ndim == 2:
- cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype)
- else:
- cc = np.zeros(lmax+1, dtype=c.dtype)
- cc[deg] = c
- c = cc
-
- # warn on rank reduction
- if rank != order and not full:
- msg = "The fit may be poorly conditioned"
- warnings.warn(msg, pu.RankWarning, stacklevel=2)
-
- if full:
- return c, [resids, rank, s, rcond]
- else:
- return c
+ return pu._fit(lagvander, x, y, deg, rcond, full, w)
def lagcompanion(c):
@@ -1666,7 +1475,8 @@ def lagroots(c):
if len(c) == 2:
return np.array([1 + c[0]/c[1]])
- m = lagcompanion(c)
+ # rotated companion matrix reduces error
+ m = lagcompanion(c)[::-1,::-1]
r = la.eigvals(m)
r.sort()
return r
@@ -1708,9 +1518,9 @@ def laggauss(deg):
the right value when integrating 1.
"""
- ideg = int(deg)
- if ideg != deg or ideg < 1:
- raise ValueError("deg must be a non-negative integer")
+ ideg = pu._deprecate_as_int(deg, "deg")
+ if ideg <= 0:
+ raise ValueError("deg must be a positive integer")
# first approximation of roots. We use the fact that the companion
# matrix is symmetric in this case in order to obtain better zeros.
diff --git a/numpy/polynomial/legendre.py b/numpy/polynomial/legendre.py
index f81bc002c..9eec9740d 100644
--- a/numpy/polynomial/legendre.py
+++ b/numpy/polynomial/legendre.py
@@ -314,21 +314,7 @@ def legfromroots(roots):
array([ 1.33333333+0.j, 0.00000000+0.j, 0.66666667+0.j]) # may vary
"""
- if len(roots) == 0:
- return np.ones(1)
- else:
- [roots] = pu.as_series([roots], trim=False)
- roots.sort()
- p = [legline(-r, 1) for r in roots]
- n = len(p)
- while n > 1:
- m, r = divmod(n, 2)
- tmp = [legmul(p[i], p[i+m]) for i in range(m)]
- if r:
- tmp[0] = legmul(tmp[0], p[-1])
- p = tmp
- n = m
- return p[0]
+ return pu._fromroots(legline, legmul, roots)
def legadd(c1, c2):
@@ -370,15 +356,7 @@ def legadd(c1, c2):
array([4., 4., 4.])
"""
- # c1, c2 are trimmed copies
- [c1, c2] = pu.as_series([c1, c2])
- if len(c1) > len(c2):
- c1[:c2.size] += c2
- ret = c1
- else:
- c2[:c1.size] += c1
- ret = c2
- return pu.trimseq(ret)
+ return pu._add(c1, c2)
def legsub(c1, c2):
@@ -422,16 +400,7 @@ def legsub(c1, c2):
array([ 2., 0., -2.])
"""
- # c1, c2 are trimmed copies
- [c1, c2] = pu.as_series([c1, c2])
- if len(c1) > len(c2):
- c1[:c2.size] -= c2
- ret = c1
- else:
- c2 = -c2
- c2[:c1.size] += c1
- ret = c2
- return pu.trimseq(ret)
+ return pu._sub(c1, c2)
def legmulx(c):
@@ -604,26 +573,7 @@ def legdiv(c1, c2):
(array([-0.07407407, 1.66666667]), array([-1.03703704, -2.51851852])) # may vary
"""
- # c1, c2 are trimmed copies
- [c1, c2] = pu.as_series([c1, c2])
- if c2[-1] == 0:
- raise ZeroDivisionError()
-
- lc1 = len(c1)
- lc2 = len(c2)
- if lc1 < lc2:
- return c1[:1]*0, c1
- elif lc2 == 1:
- return c1/c2[-1], c1[:1]*0
- else:
- quo = np.empty(lc1 - lc2 + 1, dtype=c1.dtype)
- rem = c1
- for i in range(lc1 - lc2, - 1, -1):
- p = legmul([0]*i + [1], c2)
- q = rem[-1]/p[-1]
- rem = rem[:-1] - q*p[:-1]
- quo[i] = q
- return quo, pu.trimseq(rem)
+ return pu._div(legmul, c1, c2)
def legpow(c, pow, maxpower=16):
@@ -657,24 +607,7 @@ def legpow(c, pow, maxpower=16):
--------
"""
- # c is a trimmed copy
- [c] = pu.as_series([c])
- power = int(pow)
- if power != pow or power < 0:
- raise ValueError("Power must be a non-negative integer.")
- elif maxpower is not None and power > maxpower:
- raise ValueError("Power is too large")
- elif power == 0:
- return np.array([1], dtype=c.dtype)
- elif power == 1:
- return c
- else:
- # This can be made more efficient by using powers of two
- # in the usual way.
- prd = c
- for i in range(2, power + 1):
- prd = legmul(prd, c)
- return prd
+ return pu._pow(legmul, c, pow, maxpower)
def legder(c, m=1, scl=1, axis=0):
@@ -740,14 +673,10 @@ def legder(c, m=1, scl=1, axis=0):
c = np.array(c, ndmin=1, copy=1)
if c.dtype.char in '?bBhHiIlLqQpP':
c = c.astype(np.double)
- cnt, iaxis = [int(t) for t in [m, axis]]
-
- if cnt != m:
- raise ValueError("The order of derivation must be integer")
+ cnt = pu._deprecate_as_int(m, "the order of derivation")
+ iaxis = pu._deprecate_as_int(axis, "the axis")
if cnt < 0:
raise ValueError("The order of derivation must be non-negative")
- if iaxis != axis:
- raise ValueError("The axis must be integer")
iaxis = normalize_axis_index(iaxis, c.ndim)
if cnt == 0:
@@ -863,10 +792,8 @@ def legint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
c = c.astype(np.double)
if not np.iterable(k):
k = [k]
- cnt, iaxis = [int(t) for t in [m, axis]]
-
- if cnt != m:
- raise ValueError("The order of integration must be integer")
+ cnt = pu._deprecate_as_int(m, "the order of integration")
+ iaxis = pu._deprecate_as_int(axis, "the axis")
if cnt < 0:
raise ValueError("The order of integration must be non-negative")
if len(k) > cnt:
@@ -875,8 +802,6 @@ def legint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
raise ValueError("lbnd must be a scalar.")
if np.ndim(scl) != 0:
raise ValueError("scl must be a scalar.")
- if iaxis != axis:
- raise ValueError("The axis must be integer")
iaxis = normalize_axis_index(iaxis, c.ndim)
if cnt == 0:
@@ -1039,14 +964,7 @@ def legval2d(x, y, c):
.. versionadded:: 1.7.0
"""
- try:
- x, y = np.array((x, y), copy=0)
- except Exception:
- raise ValueError('x, y are incompatible')
-
- c = legval(x, c)
- c = legval(y, c, tensor=False)
- return c
+ return pu._valnd(legval, c, x, y)
def leggrid2d(x, y, c):
@@ -1099,9 +1017,7 @@ def leggrid2d(x, y, c):
.. versionadded:: 1.7.0
"""
- c = legval(x, c)
- c = legval(y, c)
- return c
+ return pu._gridnd(legval, c, x, y)
def legval3d(x, y, z, c):
@@ -1152,15 +1068,7 @@ def legval3d(x, y, z, c):
.. versionadded:: 1.7.0
"""
- try:
- x, y, z = np.array((x, y, z), copy=0)
- except Exception:
- raise ValueError('x, y, z are incompatible')
-
- c = legval(x, c)
- c = legval(y, c, tensor=False)
- c = legval(z, c, tensor=False)
- return c
+ return pu._valnd(legval, c, x, y, z)
def leggrid3d(x, y, z, c):
@@ -1216,10 +1124,7 @@ def leggrid3d(x, y, z, c):
.. versionadded:: 1.7.0
"""
- c = legval(x, c)
- c = legval(y, c)
- c = legval(z, c)
- return c
+ return pu._gridnd(legval, c, x, y, z)
def legvander(x, deg):
@@ -1257,9 +1162,7 @@ def legvander(x, deg):
the converted `x`.
"""
- ideg = int(deg)
- if ideg != deg:
- raise ValueError("deg must be integer")
+ ideg = pu._deprecate_as_int(deg, "deg")
if ideg < 0:
raise ValueError("deg must be non-negative")
@@ -1327,17 +1230,7 @@ def legvander2d(x, y, deg):
.. versionadded:: 1.7.0
"""
- ideg = [int(d) for d in deg]
- is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
- if is_valid != [1, 1]:
- raise ValueError("degrees must be non-negative integers")
- degx, degy = ideg
- x, y = np.array((x, y), copy=0) + 0.0
-
- vx = legvander(x, degx)
- vy = legvander(y, degy)
- v = vx[..., None]*vy[..., None,:]
- return v.reshape(v.shape[:-2] + (-1,))
+ return pu._vander2d(legvander, x, y, deg)
def legvander3d(x, y, z, deg):
@@ -1391,18 +1284,7 @@ def legvander3d(x, y, z, deg):
.. versionadded:: 1.7.0
"""
- ideg = [int(d) for d in deg]
- is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
- if is_valid != [1, 1, 1]:
- raise ValueError("degrees must be non-negative integers")
- degx, degy, degz = ideg
- x, y, z = np.array((x, y, z), copy=0) + 0.0
-
- vx = legvander(x, degx)
- vy = legvander(y, degy)
- vz = legvander(z, degz)
- v = vx[..., None, None]*vy[..., None,:, None]*vz[..., None, None,:]
- return v.reshape(v.shape[:-3] + (-1,))
+ return pu._vander3d(legvander, x, y, z, deg)
def legfit(x, y, deg, rcond=None, full=False, w=None):
@@ -1526,81 +1408,7 @@ def legfit(x, y, deg, rcond=None, full=False, w=None):
--------
"""
- x = np.asarray(x) + 0.0
- y = np.asarray(y) + 0.0
- deg = np.asarray(deg)
-
- # check arguments.
- if deg.ndim > 1 or deg.dtype.kind not in 'iu' or deg.size == 0:
- raise TypeError("deg must be an int or non-empty 1-D array of int")
- if deg.min() < 0:
- raise ValueError("expected deg >= 0")
- if x.ndim != 1:
- raise TypeError("expected 1D vector for x")
- if x.size == 0:
- raise TypeError("expected non-empty vector for x")
- if y.ndim < 1 or y.ndim > 2:
- raise TypeError("expected 1D or 2D array for y")
- if len(x) != len(y):
- raise TypeError("expected x and y to have same length")
-
- if deg.ndim == 0:
- lmax = deg
- order = lmax + 1
- van = legvander(x, lmax)
- else:
- deg = np.sort(deg)
- lmax = deg[-1]
- order = len(deg)
- van = legvander(x, lmax)[:, deg]
-
- # set up the least squares matrices in transposed form
- lhs = van.T
- rhs = y.T
- if w is not None:
- w = np.asarray(w) + 0.0
- if w.ndim != 1:
- raise TypeError("expected 1D vector for w")
- if len(x) != len(w):
- raise TypeError("expected x and w to have same length")
- # apply weights. Don't use inplace operations as they
- # can cause problems with NA.
- lhs = lhs * w
- rhs = rhs * w
-
- # set rcond
- if rcond is None:
- rcond = len(x)*np.finfo(x.dtype).eps
-
- # Determine the norms of the design matrix columns.
- if issubclass(lhs.dtype.type, np.complexfloating):
- scl = np.sqrt((np.square(lhs.real) + np.square(lhs.imag)).sum(1))
- else:
- scl = np.sqrt(np.square(lhs).sum(1))
- scl[scl == 0] = 1
-
- # Solve the least squares problem.
- c, resids, rank, s = la.lstsq(lhs.T/scl, rhs.T, rcond)
- c = (c.T/scl).T
-
- # Expand c to include non-fitted coefficients which are set to zero
- if deg.ndim > 0:
- if c.ndim == 2:
- cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype)
- else:
- cc = np.zeros(lmax+1, dtype=c.dtype)
- cc[deg] = c
- c = cc
-
- # warn on rank reduction
- if rank != order and not full:
- msg = "The fit may be poorly conditioned"
- warnings.warn(msg, pu.RankWarning, stacklevel=2)
-
- if full:
- return c, [resids, rank, s, rcond]
- else:
- return c
+ return pu._fit(legvander, x, y, deg, rcond, full, w)
def legcompanion(c):
@@ -1697,7 +1505,8 @@ def legroots(c):
if len(c) == 2:
return np.array([-c[0]/c[1]])
- m = legcompanion(c)
+ # rotated companion matrix reduces error
+ m = legcompanion(c)[::-1,::-1]
r = la.eigvals(m)
r.sort()
return r
@@ -1739,9 +1548,9 @@ def leggauss(deg):
the right value when integrating 1.
"""
- ideg = int(deg)
- if ideg != deg or ideg < 1:
- raise ValueError("deg must be a non-negative integer")
+ ideg = pu._deprecate_as_int(deg, "deg")
+ if ideg <= 0:
+ raise ValueError("deg must be a positive integer")
# first approximation of roots. We use the fact that the companion
# matrix is symmetric in this case in order to obtain better zeros.
diff --git a/numpy/polynomial/polynomial.py b/numpy/polynomial/polynomial.py
index 69599e3fd..afa330502 100644
--- a/numpy/polynomial/polynomial.py
+++ b/numpy/polynomial/polynomial.py
@@ -188,21 +188,7 @@ def polyfromroots(roots):
array([1.+0.j, 0.+0.j, 1.+0.j])
"""
- if len(roots) == 0:
- return np.ones(1)
- else:
- [roots] = pu.as_series([roots], trim=False)
- roots.sort()
- p = [polyline(-r, 1) for r in roots]
- n = len(p)
- while n > 1:
- m, r = divmod(n, 2)
- tmp = [polymul(p[i], p[i+m]) for i in range(m)]
- if r:
- tmp[0] = polymul(tmp[0], p[-1])
- p = tmp
- n = m
- return p[0]
+ return pu._fromroots(polyline, polymul, roots)
def polyadd(c1, c2):
@@ -238,15 +224,7 @@ def polyadd(c1, c2):
28.0
"""
- # c1, c2 are trimmed copies
- [c1, c2] = pu.as_series([c1, c2])
- if len(c1) > len(c2):
- c1[:c2.size] += c2
- ret = c1
- else:
- c2[:c1.size] += c1
- ret = c2
- return pu.trimseq(ret)
+ return pu._add(c1, c2)
def polysub(c1, c2):
@@ -283,16 +261,7 @@ def polysub(c1, c2):
array([ 2., 0., -2.])
"""
- # c1, c2 are trimmed copies
- [c1, c2] = pu.as_series([c1, c2])
- if len(c1) > len(c2):
- c1[:c2.size] -= c2
- ret = c1
- else:
- c2 = -c2
- c2[:c1.size] += c1
- ret = c2
- return pu.trimseq(ret)
+ return pu._sub(c1, c2)
def polymulx(c):
@@ -411,18 +380,19 @@ def polydiv(c1, c2):
if c2[-1] == 0:
raise ZeroDivisionError()
- len1 = len(c1)
- len2 = len(c2)
- if len2 == 1:
- return c1/c2[-1], c1[:1]*0
- elif len1 < len2:
+ # note: this is more efficient than `pu._div(polymul, c1, c2)`
+ lc1 = len(c1)
+ lc2 = len(c2)
+ if lc1 < lc2:
return c1[:1]*0, c1
+ elif lc2 == 1:
+ return c1/c2[-1], c1[:1]*0
else:
- dlen = len1 - len2
+ dlen = lc1 - lc2
scl = c2[-1]
c2 = c2[:-1]/scl
i = dlen
- j = len1 - 1
+ j = lc1 - 1
while i >= 0:
c1[i:j] -= c2*c1[j]
i -= 1
@@ -464,24 +434,9 @@ def polypow(c, pow, maxpower=None):
array([ 1., 4., 10., 12., 9.])
"""
- # c is a trimmed copy
- [c] = pu.as_series([c])
- power = int(pow)
- if power != pow or power < 0:
- raise ValueError("Power must be a non-negative integer.")
- elif maxpower is not None and power > maxpower:
- raise ValueError("Power is too large")
- elif power == 0:
- return np.array([1], dtype=c.dtype)
- elif power == 1:
- return c
- else:
- # This can be made more efficient by using powers of two
- # in the usual way.
- prd = c
- for i in range(2, power + 1):
- prd = np.convolve(prd, c)
- return prd
+ # note: this is more efficient than `pu._pow(polymul, c1, c2)`, as it
+ # avoids calling `as_series` repeatedly
+ return pu._pow(np.convolve, c, pow, maxpower)
def polyder(c, m=1, scl=1, axis=0):
@@ -541,14 +496,10 @@ def polyder(c, m=1, scl=1, axis=0):
# astype fails with NA
c = c + 0.0
cdt = c.dtype
- cnt, iaxis = [int(t) for t in [m, axis]]
-
- if cnt != m:
- raise ValueError("The order of derivation must be integer")
+ cnt = pu._deprecate_as_int(m, "the order of derivation")
+ iaxis = pu._deprecate_as_int(axis, "the axis")
if cnt < 0:
raise ValueError("The order of derivation must be non-negative")
- if iaxis != axis:
- raise ValueError("The axis must be integer")
iaxis = normalize_axis_index(iaxis, c.ndim)
if cnt == 0:
@@ -655,10 +606,8 @@ def polyint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
cdt = c.dtype
if not np.iterable(k):
k = [k]
- cnt, iaxis = [int(t) for t in [m, axis]]
-
- if cnt != m:
- raise ValueError("The order of integration must be integer")
+ cnt = pu._deprecate_as_int(m, "the order of integration")
+ iaxis = pu._deprecate_as_int(axis, "the axis")
if cnt < 0:
raise ValueError("The order of integration must be non-negative")
if len(k) > cnt:
@@ -667,8 +616,6 @@ def polyint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
raise ValueError("lbnd must be a scalar.")
if np.ndim(scl) != 0:
raise ValueError("scl must be a scalar.")
- if iaxis != axis:
- raise ValueError("The axis must be integer")
iaxis = normalize_axis_index(iaxis, c.ndim)
if cnt == 0:
@@ -924,14 +871,7 @@ def polyval2d(x, y, c):
.. versionadded:: 1.7.0
"""
- try:
- x, y = np.array((x, y), copy=0)
- except Exception:
- raise ValueError('x, y are incompatible')
-
- c = polyval(x, c)
- c = polyval(y, c, tensor=False)
- return c
+ return pu._valnd(polyval, c, x, y)
def polygrid2d(x, y, c):
@@ -984,9 +924,7 @@ def polygrid2d(x, y, c):
.. versionadded:: 1.7.0
"""
- c = polyval(x, c)
- c = polyval(y, c)
- return c
+ return pu._gridnd(polyval, c, x, y)
def polyval3d(x, y, z, c):
@@ -1037,15 +975,7 @@ def polyval3d(x, y, z, c):
.. versionadded:: 1.7.0
"""
- try:
- x, y, z = np.array((x, y, z), copy=0)
- except Exception:
- raise ValueError('x, y, z are incompatible')
-
- c = polyval(x, c)
- c = polyval(y, c, tensor=False)
- c = polyval(z, c, tensor=False)
- return c
+ return pu._valnd(polyval, c, x, y, z)
def polygrid3d(x, y, z, c):
@@ -1101,10 +1031,7 @@ def polygrid3d(x, y, z, c):
.. versionadded:: 1.7.0
"""
- c = polyval(x, c)
- c = polyval(y, c)
- c = polyval(z, c)
- return c
+ return pu._gridnd(polyval, c, x, y, z)
def polyvander(x, deg):
@@ -1145,9 +1072,7 @@ def polyvander(x, deg):
polyvander2d, polyvander3d
"""
- ideg = int(deg)
- if ideg != deg:
- raise ValueError("deg must be integer")
+ ideg = pu._deprecate_as_int(deg, "deg")
if ideg < 0:
raise ValueError("deg must be non-negative")
@@ -1208,19 +1133,7 @@ def polyvander2d(x, y, deg):
polyvander, polyvander3d. polyval2d, polyval3d
"""
- ideg = [int(d) for d in deg]
- is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
- if is_valid != [1, 1]:
- raise ValueError("degrees must be non-negative integers")
- degx, degy = ideg
- x, y = np.array((x, y), copy=0) + 0.0
-
- vx = polyvander(x, degx)
- vy = polyvander(y, degy)
- v = vx[..., None]*vy[..., None,:]
- # einsum bug
- #v = np.einsum("...i,...j->...ij", vx, vy)
- return v.reshape(v.shape[:-2] + (-1,))
+ return pu._vander2d(polyvander, x, y, deg)
def polyvander3d(x, y, z, deg):
@@ -1274,20 +1187,7 @@ def polyvander3d(x, y, z, deg):
.. versionadded:: 1.7.0
"""
- ideg = [int(d) for d in deg]
- is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
- if is_valid != [1, 1, 1]:
- raise ValueError("degrees must be non-negative integers")
- degx, degy, degz = ideg
- x, y, z = np.array((x, y, z), copy=0) + 0.0
-
- vx = polyvander(x, degx)
- vy = polyvander(y, degy)
- vz = polyvander(z, degz)
- v = vx[..., None, None]*vy[..., None,:, None]*vz[..., None, None,:]
- # einsum bug
- #v = np.einsum("...i, ...j, ...k->...ijk", vx, vy, vz)
- return v.reshape(v.shape[:-3] + (-1,))
+ return pu._vander3d(polyvander, x, y, z, deg)
def polyfit(x, y, deg, rcond=None, full=False, w=None):
@@ -1433,81 +1333,7 @@ def polyfit(x, y, deg, rcond=None, full=False, w=None):
0.50443316, 0.28853036]), 1.1324274851176597e-014]
"""
- x = np.asarray(x) + 0.0
- y = np.asarray(y) + 0.0
- deg = np.asarray(deg)
-
- # check arguments.
- if deg.ndim > 1 or deg.dtype.kind not in 'iu' or deg.size == 0:
- raise TypeError("deg must be an int or non-empty 1-D array of int")
- if deg.min() < 0:
- raise ValueError("expected deg >= 0")
- if x.ndim != 1:
- raise TypeError("expected 1D vector for x")
- if x.size == 0:
- raise TypeError("expected non-empty vector for x")
- if y.ndim < 1 or y.ndim > 2:
- raise TypeError("expected 1D or 2D array for y")
- if len(x) != len(y):
- raise TypeError("expected x and y to have same length")
-
- if deg.ndim == 0:
- lmax = deg
- order = lmax + 1
- van = polyvander(x, lmax)
- else:
- deg = np.sort(deg)
- lmax = deg[-1]
- order = len(deg)
- van = polyvander(x, lmax)[:, deg]
-
- # set up the least squares matrices in transposed form
- lhs = van.T
- rhs = y.T
- if w is not None:
- w = np.asarray(w) + 0.0
- if w.ndim != 1:
- raise TypeError("expected 1D vector for w")
- if len(x) != len(w):
- raise TypeError("expected x and w to have same length")
- # apply weights. Don't use inplace operations as they
- # can cause problems with NA.
- lhs = lhs * w
- rhs = rhs * w
-
- # set rcond
- if rcond is None:
- rcond = len(x)*np.finfo(x.dtype).eps
-
- # Determine the norms of the design matrix columns.
- if issubclass(lhs.dtype.type, np.complexfloating):
- scl = np.sqrt((np.square(lhs.real) + np.square(lhs.imag)).sum(1))
- else:
- scl = np.sqrt(np.square(lhs).sum(1))
- scl[scl == 0] = 1
-
- # Solve the least squares problem.
- c, resids, rank, s = la.lstsq(lhs.T/scl, rhs.T, rcond)
- c = (c.T/scl).T
-
- # Expand c to include non-fitted coefficients which are set to zero
- if deg.ndim == 1:
- if c.ndim == 2:
- cc = np.zeros((lmax + 1, c.shape[1]), dtype=c.dtype)
- else:
- cc = np.zeros(lmax + 1, dtype=c.dtype)
- cc[deg] = c
- c = cc
-
- # warn on rank reduction
- if rank != order and not full:
- msg = "The fit may be poorly conditioned"
- warnings.warn(msg, pu.RankWarning, stacklevel=2)
-
- if full:
- return c, [resids, rank, s, rcond]
- else:
- return c
+ return pu._fit(polyvander, x, y, deg, rcond, full, w)
def polycompanion(c):
@@ -1602,7 +1428,8 @@ def polyroots(c):
if len(c) == 2:
return np.array([-c[0]/c[1]])
- m = polycompanion(c)
+ # rotated companion matrix reduces error
+ m = polycompanion(c)[::-1,::-1]
r = la.eigvals(m)
r.sort()
return r
diff --git a/numpy/polynomial/polyutils.py b/numpy/polynomial/polyutils.py
index eff4a8ee1..5128e35e0 100644
--- a/numpy/polynomial/polyutils.py
+++ b/numpy/polynomial/polyutils.py
@@ -45,6 +45,8 @@ Functions
"""
from __future__ import division, absolute_import, print_function
+import operator
+
import numpy as np
__all__ = [
@@ -410,3 +412,350 @@ def mapdomain(x, old, new):
x = np.asanyarray(x)
off, scl = mapparms(old, new)
return off + scl*x
+
+
+def _vander2d(vander_f, x, y, deg):
+ """
+ Helper function used to implement the ``<type>vander2d`` functions.
+
+ Parameters
+ ----------
+ vander_f : function(array_like, int) -> ndarray
+ The 1d vander function, such as ``polyvander``
+ x, y, deg :
+ See the ``<type>vander2d`` functions for more detail
+ """
+ degx, degy = [
+ _deprecate_as_int(d, "degrees")
+ for d in deg
+ ]
+ x, y = np.array((x, y), copy=0) + 0.0
+
+ vx = vander_f(x, degx)
+ vy = vander_f(y, degy)
+ v = vx[..., None]*vy[..., None,:]
+ return v.reshape(v.shape[:-2] + (-1,))
+
+
+def _vander3d(vander_f, x, y, z, deg):
+ """
+ Helper function used to implement the ``<type>vander3d`` functions.
+
+ Parameters
+ ----------
+ vander_f : function(array_like, int) -> ndarray
+ The 1d vander function, such as ``polyvander``
+ x, y, z, deg :
+ See the ``<type>vander3d`` functions for more detail
+ """
+ degx, degy, degz = [
+ _deprecate_as_int(d, "degrees")
+ for d in deg
+ ]
+ x, y, z = np.array((x, y, z), copy=0) + 0.0
+
+ vx = vander_f(x, degx)
+ vy = vander_f(y, degy)
+ vz = vander_f(z, degz)
+ v = vx[..., None, None]*vy[..., None,:, None]*vz[..., None, None,:]
+ return v.reshape(v.shape[:-3] + (-1,))
+
+
+def _fromroots(line_f, mul_f, roots):
+ """
+ Helper function used to implement the ``<type>fromroots`` functions.
+
+ Parameters
+ ----------
+ line_f : function(float, float) -> ndarray
+ The ``<type>line`` function, such as ``polyline``
+ mul_f : function(array_like, array_like) -> ndarray
+ The ``<type>mul`` function, such as ``polymul``
+ roots :
+ See the ``<type>fromroots`` functions for more detail
+ """
+ if len(roots) == 0:
+ return np.ones(1)
+ else:
+ [roots] = as_series([roots], trim=False)
+ roots.sort()
+ p = [line_f(-r, 1) for r in roots]
+ n = len(p)
+ while n > 1:
+ m, r = divmod(n, 2)
+ tmp = [mul_f(p[i], p[i+m]) for i in range(m)]
+ if r:
+ tmp[0] = mul_f(tmp[0], p[-1])
+ p = tmp
+ n = m
+ return p[0]
+
+
+def _valnd(val_f, c, *args):
+ """
+ Helper function used to implement the ``<type>val<n>d`` functions.
+
+ Parameters
+ ----------
+ val_f : function(array_like, array_like, tensor: bool) -> array_like
+ The ``<type>val`` function, such as ``polyval``
+ c, args :
+ See the ``<type>val<n>d`` functions for more detail
+ """
+ try:
+ args = tuple(np.array(args, copy=False))
+ except Exception:
+ # preserve the old error message
+ if len(args) == 2:
+ raise ValueError('x, y, z are incompatible')
+ elif len(args) == 3:
+ raise ValueError('x, y are incompatible')
+ else:
+ raise ValueError('ordinates are incompatible')
+
+ it = iter(args)
+ x0 = next(it)
+
+ # use tensor on only the first
+ c = val_f(x0, c)
+ for xi in it:
+ c = val_f(xi, c, tensor=False)
+ return c
+
+
+def _gridnd(val_f, c, *args):
+ """
+ Helper function used to implement the ``<type>grid<n>d`` functions.
+
+ Parameters
+ ----------
+ val_f : function(array_like, array_like, tensor: bool) -> array_like
+ The ``<type>val`` function, such as ``polyval``
+ c, args :
+ See the ``<type>grid<n>d`` functions for more detail
+ """
+ for xi in args:
+ c = val_f(xi, c)
+ return c
+
+
+def _div(mul_f, c1, c2):
+ """
+ Helper function used to implement the ``<type>div`` functions.
+
+ Implementation uses repeated subtraction of c2 multiplied by the nth basis.
+ For some polynomial types, a more efficient approach may be possible.
+
+ Parameters
+ ----------
+ mul_f : function(array_like, array_like) -> array_like
+ The ``<type>mul`` function, such as ``polymul``
+ c1, c2 :
+ See the ``<type>div`` functions for more detail
+ """
+ # c1, c2 are trimmed copies
+ [c1, c2] = as_series([c1, c2])
+ if c2[-1] == 0:
+ raise ZeroDivisionError()
+
+ lc1 = len(c1)
+ lc2 = len(c2)
+ if lc1 < lc2:
+ return c1[:1]*0, c1
+ elif lc2 == 1:
+ return c1/c2[-1], c1[:1]*0
+ else:
+ quo = np.empty(lc1 - lc2 + 1, dtype=c1.dtype)
+ rem = c1
+ for i in range(lc1 - lc2, - 1, -1):
+ p = mul_f([0]*i + [1], c2)
+ q = rem[-1]/p[-1]
+ rem = rem[:-1] - q*p[:-1]
+ quo[i] = q
+ return quo, trimseq(rem)
+
+
+def _add(c1, c2):
+ """ Helper function used to implement the ``<type>add`` functions. """
+ # c1, c2 are trimmed copies
+ [c1, c2] = as_series([c1, c2])
+ if len(c1) > len(c2):
+ c1[:c2.size] += c2
+ ret = c1
+ else:
+ c2[:c1.size] += c1
+ ret = c2
+ return trimseq(ret)
+
+
+def _sub(c1, c2):
+ """ Helper function used to implement the ``<type>sub`` functions. """
+ # c1, c2 are trimmed copies
+ [c1, c2] = as_series([c1, c2])
+ if len(c1) > len(c2):
+ c1[:c2.size] -= c2
+ ret = c1
+ else:
+ c2 = -c2
+ c2[:c1.size] += c1
+ ret = c2
+ return trimseq(ret)
+
+
+def _fit(vander_f, x, y, deg, rcond=None, full=False, w=None):
+ """
+ Helper function used to implement the ``<type>fit`` functions.
+
+ Parameters
+ ----------
+ vander_f : function(array_like, int) -> ndarray
+ The 1d vander function, such as ``polyvander``
+ c1, c2 :
+ See the ``<type>fit`` functions for more detail
+ """
+ x = np.asarray(x) + 0.0
+ y = np.asarray(y) + 0.0
+ deg = np.asarray(deg)
+
+ # check arguments.
+ if deg.ndim > 1 or deg.dtype.kind not in 'iu' or deg.size == 0:
+ raise TypeError("deg must be an int or non-empty 1-D array of int")
+ if deg.min() < 0:
+ raise ValueError("expected deg >= 0")
+ if x.ndim != 1:
+ raise TypeError("expected 1D vector for x")
+ if x.size == 0:
+ raise TypeError("expected non-empty vector for x")
+ if y.ndim < 1 or y.ndim > 2:
+ raise TypeError("expected 1D or 2D array for y")
+ if len(x) != len(y):
+ raise TypeError("expected x and y to have same length")
+
+ if deg.ndim == 0:
+ lmax = deg
+ order = lmax + 1
+ van = vander_f(x, lmax)
+ else:
+ deg = np.sort(deg)
+ lmax = deg[-1]
+ order = len(deg)
+ van = vander_f(x, lmax)[:, deg]
+
+ # set up the least squares matrices in transposed form
+ lhs = van.T
+ rhs = y.T
+ if w is not None:
+ w = np.asarray(w) + 0.0
+ if w.ndim != 1:
+ raise TypeError("expected 1D vector for w")
+ if len(x) != len(w):
+ raise TypeError("expected x and w to have same length")
+ # apply weights. Don't use inplace operations as they
+ # can cause problems with NA.
+ lhs = lhs * w
+ rhs = rhs * w
+
+ # set rcond
+ if rcond is None:
+ rcond = len(x)*np.finfo(x.dtype).eps
+
+ # Determine the norms of the design matrix columns.
+ if issubclass(lhs.dtype.type, np.complexfloating):
+ scl = np.sqrt((np.square(lhs.real) + np.square(lhs.imag)).sum(1))
+ else:
+ scl = np.sqrt(np.square(lhs).sum(1))
+ scl[scl == 0] = 1
+
+ # Solve the least squares problem.
+ c, resids, rank, s = np.linalg.lstsq(lhs.T/scl, rhs.T, rcond)
+ c = (c.T/scl).T
+
+ # Expand c to include non-fitted coefficients which are set to zero
+ if deg.ndim > 0:
+ if c.ndim == 2:
+ cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype)
+ else:
+ cc = np.zeros(lmax+1, dtype=c.dtype)
+ cc[deg] = c
+ c = cc
+
+ # warn on rank reduction
+ if rank != order and not full:
+ msg = "The fit may be poorly conditioned"
+ warnings.warn(msg, RankWarning, stacklevel=2)
+
+ if full:
+ return c, [resids, rank, s, rcond]
+ else:
+ return c
+
+
+def _pow(mul_f, c, pow, maxpower):
+ """
+ Helper function used to implement the ``<type>pow`` functions.
+
+ Parameters
+ ----------
+ vander_f : function(array_like, int) -> ndarray
+ The 1d vander function, such as ``polyvander``
+ pow, maxpower :
+ See the ``<type>pow`` functions for more detail
+ mul_f : function(array_like, array_like) -> ndarray
+ The ``<type>mul`` function, such as ``polymul``
+ """
+ # c is a trimmed copy
+ [c] = as_series([c])
+ power = int(pow)
+ if power != pow or power < 0:
+ raise ValueError("Power must be a non-negative integer.")
+ elif maxpower is not None and power > maxpower:
+ raise ValueError("Power is too large")
+ elif power == 0:
+ return np.array([1], dtype=c.dtype)
+ elif power == 1:
+ return c
+ else:
+ # This can be made more efficient by using powers of two
+ # in the usual way.
+ prd = c
+ for i in range(2, power + 1):
+ prd = mul_f(prd, c)
+ return prd
+
+
+def _deprecate_as_int(x, desc):
+ """
+ Like `operator.index`, but emits a deprecation warning when passed a float
+
+ Parameters
+ ----------
+ x : int-like, or float with integral value
+ Value to interpret as an integer
+ desc : str
+ description to include in any error message
+
+ Raises
+ ------
+ TypeError : if x is a non-integral float or non-numeric
+ DeprecationWarning : if x is an integral float
+ """
+ try:
+ return operator.index(x)
+ except TypeError:
+ # Numpy 1.17.0, 2019-03-11
+ try:
+ ix = int(x)
+ except TypeError:
+ pass
+ else:
+ if ix == x:
+ warnings.warn(
+ "In future, this will raise TypeError, as {} will need to "
+ "be an integer not just an integral float."
+ .format(desc),
+ DeprecationWarning,
+ stacklevel=3
+ )
+ return ix
+
+ raise TypeError("{} must be an integer".format(desc))
diff --git a/numpy/polynomial/tests/test_chebyshev.py b/numpy/polynomial/tests/test_chebyshev.py
index 7fb7492c6..c8d2d6dba 100644
--- a/numpy/polynomial/tests/test_chebyshev.py
+++ b/numpy/polynomial/tests/test_chebyshev.py
@@ -221,12 +221,12 @@ class TestIntegral(object):
def test_chebint(self):
# check exceptions
- assert_raises(ValueError, cheb.chebint, [0], .5)
+ assert_raises(TypeError, cheb.chebint, [0], .5)
assert_raises(ValueError, cheb.chebint, [0], -1)
assert_raises(ValueError, cheb.chebint, [0], 1, [0, 0])
assert_raises(ValueError, cheb.chebint, [0], lbnd=[0])
assert_raises(ValueError, cheb.chebint, [0], scl=[0])
- assert_raises(ValueError, cheb.chebint, [0], axis=.5)
+ assert_raises(TypeError, cheb.chebint, [0], axis=.5)
# test integration of zero polynomial
for i in range(2, 5):
@@ -323,7 +323,7 @@ class TestDerivative(object):
def test_chebder(self):
# check exceptions
- assert_raises(ValueError, cheb.chebder, [0], .5)
+ assert_raises(TypeError, cheb.chebder, [0], .5)
assert_raises(ValueError, cheb.chebder, [0], -1)
# check that zeroth derivative does nothing
diff --git a/numpy/polynomial/tests/test_hermite.py b/numpy/polynomial/tests/test_hermite.py
index 1287ef3fe..271c1964b 100644
--- a/numpy/polynomial/tests/test_hermite.py
+++ b/numpy/polynomial/tests/test_hermite.py
@@ -209,12 +209,12 @@ class TestIntegral(object):
def test_hermint(self):
# check exceptions
- assert_raises(ValueError, herm.hermint, [0], .5)
+ assert_raises(TypeError, herm.hermint, [0], .5)
assert_raises(ValueError, herm.hermint, [0], -1)
assert_raises(ValueError, herm.hermint, [0], 1, [0, 0])
assert_raises(ValueError, herm.hermint, [0], lbnd=[0])
assert_raises(ValueError, herm.hermint, [0], scl=[0])
- assert_raises(ValueError, herm.hermint, [0], axis=.5)
+ assert_raises(TypeError, herm.hermint, [0], axis=.5)
# test integration of zero polynomial
for i in range(2, 5):
@@ -311,7 +311,7 @@ class TestDerivative(object):
def test_hermder(self):
# check exceptions
- assert_raises(ValueError, herm.hermder, [0], .5)
+ assert_raises(TypeError, herm.hermder, [0], .5)
assert_raises(ValueError, herm.hermder, [0], -1)
# check that zeroth derivative does nothing
diff --git a/numpy/polynomial/tests/test_hermite_e.py b/numpy/polynomial/tests/test_hermite_e.py
index ccb44ad73..434b30e7b 100644
--- a/numpy/polynomial/tests/test_hermite_e.py
+++ b/numpy/polynomial/tests/test_hermite_e.py
@@ -209,12 +209,12 @@ class TestIntegral(object):
def test_hermeint(self):
# check exceptions
- assert_raises(ValueError, herme.hermeint, [0], .5)
+ assert_raises(TypeError, herme.hermeint, [0], .5)
assert_raises(ValueError, herme.hermeint, [0], -1)
assert_raises(ValueError, herme.hermeint, [0], 1, [0, 0])
assert_raises(ValueError, herme.hermeint, [0], lbnd=[0])
assert_raises(ValueError, herme.hermeint, [0], scl=[0])
- assert_raises(ValueError, herme.hermeint, [0], axis=.5)
+ assert_raises(TypeError, herme.hermeint, [0], axis=.5)
# test integration of zero polynomial
for i in range(2, 5):
@@ -311,7 +311,7 @@ class TestDerivative(object):
def test_hermeder(self):
# check exceptions
- assert_raises(ValueError, herme.hermeder, [0], .5)
+ assert_raises(TypeError, herme.hermeder, [0], .5)
assert_raises(ValueError, herme.hermeder, [0], -1)
# check that zeroth derivative does nothing
diff --git a/numpy/polynomial/tests/test_laguerre.py b/numpy/polynomial/tests/test_laguerre.py
index 3ababec5e..4b9b28637 100644
--- a/numpy/polynomial/tests/test_laguerre.py
+++ b/numpy/polynomial/tests/test_laguerre.py
@@ -206,12 +206,12 @@ class TestIntegral(object):
def test_lagint(self):
# check exceptions
- assert_raises(ValueError, lag.lagint, [0], .5)
+ assert_raises(TypeError, lag.lagint, [0], .5)
assert_raises(ValueError, lag.lagint, [0], -1)
assert_raises(ValueError, lag.lagint, [0], 1, [0, 0])
assert_raises(ValueError, lag.lagint, [0], lbnd=[0])
assert_raises(ValueError, lag.lagint, [0], scl=[0])
- assert_raises(ValueError, lag.lagint, [0], axis=.5)
+ assert_raises(TypeError, lag.lagint, [0], axis=.5)
# test integration of zero polynomial
for i in range(2, 5):
@@ -308,7 +308,7 @@ class TestDerivative(object):
def test_lagder(self):
# check exceptions
- assert_raises(ValueError, lag.lagder, [0], .5)
+ assert_raises(TypeError, lag.lagder, [0], .5)
assert_raises(ValueError, lag.lagder, [0], -1)
# check that zeroth derivative does nothing
diff --git a/numpy/polynomial/tests/test_legendre.py b/numpy/polynomial/tests/test_legendre.py
index a23086d59..917a7e03a 100644
--- a/numpy/polynomial/tests/test_legendre.py
+++ b/numpy/polynomial/tests/test_legendre.py
@@ -210,12 +210,12 @@ class TestIntegral(object):
def test_legint(self):
# check exceptions
- assert_raises(ValueError, leg.legint, [0], .5)
+ assert_raises(TypeError, leg.legint, [0], .5)
assert_raises(ValueError, leg.legint, [0], -1)
assert_raises(ValueError, leg.legint, [0], 1, [0, 0])
assert_raises(ValueError, leg.legint, [0], lbnd=[0])
assert_raises(ValueError, leg.legint, [0], scl=[0])
- assert_raises(ValueError, leg.legint, [0], axis=.5)
+ assert_raises(TypeError, leg.legint, [0], axis=.5)
# test integration of zero polynomial
for i in range(2, 5):
@@ -312,7 +312,7 @@ class TestDerivative(object):
def test_legder(self):
# check exceptions
- assert_raises(ValueError, leg.legder, [0], .5)
+ assert_raises(TypeError, leg.legder, [0], .5)
assert_raises(ValueError, leg.legder, [0], -1)
# check that zeroth derivative does nothing
diff --git a/numpy/polynomial/tests/test_polynomial.py b/numpy/polynomial/tests/test_polynomial.py
index 562aa904d..08b73da15 100644
--- a/numpy/polynomial/tests/test_polynomial.py
+++ b/numpy/polynomial/tests/test_polynomial.py
@@ -291,12 +291,12 @@ class TestIntegral(object):
def test_polyint(self):
# check exceptions
- assert_raises(ValueError, poly.polyint, [0], .5)
+ assert_raises(TypeError, poly.polyint, [0], .5)
assert_raises(ValueError, poly.polyint, [0], -1)
assert_raises(ValueError, poly.polyint, [0], 1, [0, 0])
assert_raises(ValueError, poly.polyint, [0], lbnd=[0])
assert_raises(ValueError, poly.polyint, [0], scl=[0])
- assert_raises(ValueError, poly.polyint, [0], axis=.5)
+ assert_raises(TypeError, poly.polyint, [0], axis=.5)
# test integration of zero polynomial
for i in range(2, 5):
@@ -388,7 +388,7 @@ class TestDerivative(object):
def test_polyder(self):
# check exceptions
- assert_raises(ValueError, poly.polyder, [0], .5)
+ assert_raises(TypeError, poly.polyder, [0], .5)
assert_raises(ValueError, poly.polyder, [0], -1)
# check that zeroth derivative does nothing
View Lorry Controller Status