NumPy 源码解析（三十三）
.\numpy\numpy\polynomial\legendre.py
"""
==================================================
Legendre Series (:mod:`numpy.polynomial.legendre`)
==================================================
This module provides a number of objects (mostly functions) useful for
dealing with Legendre series, including a `Legendre` class that
encapsulates the usual arithmetic operations. (General information
on how this module represents and works with such polynomials is in the
docstring for its "parent" sub-package, `numpy.polynomial`).
Classes
-------
.. autosummary::
:toctree: generated/
Legendre
Constants
---------
.. autosummary::
:toctree: generated/
legdomain
legzero
legone
legx
Arithmetic
----------
.. autosummary::
:toctree: generated/
legadd
legsub
legmulx
legmul
legdiv
legpow
legval
legval2d
legval3d
leggrid2d
leggrid3d
Calculus
--------
.. autosummary::
:toctree: generated/
legder
legint
Misc Functions
--------------
.. autosummary::
:toctree: generated/
legfromroots
legroots
legvander
legvander2d
legvander3d
leggauss
legweight
legcompanion
legfit
legtrim
legline
leg2poly
poly2leg
See also
--------
numpy.polynomial
"""
import numpy as np
import numpy.linalg as la
from numpy.lib.array_utils import normalize_axis_index
from . import polyutils as pu
from ._polybase import ABCPolyBase
__all__ = [
'legzero', 'legone', 'legx', 'legdomain', 'legline', 'legadd',
'legsub', 'legmulx', 'legmul', 'legdiv', 'legpow', 'legval', 'legder',
'legint', 'leg2poly', 'poly2leg', 'legfromroots', 'legvander',
'legfit', 'legtrim', 'legroots', 'Legendre', 'legval2d', 'legval3d',
'leggrid2d', 'leggrid3d', 'legvander2d', 'legvander3d', 'legcompanion',
'leggauss', 'legweight']
legtrim = pu.trimcoef
def poly2leg(pol):
"""
Convert a polynomial to a Legendre series.
Convert an array representing the coefficients of a polynomial (relative
to the "standard" basis) ordered from lowest degree to highest, to an
array of the coefficients of the equivalent Legendre series, ordered
from lowest to highest degree.
Parameters
----------
pol : array_like
1-D array containing the polynomial coefficients
Returns
-------
c : ndarray
1-D array containing the coefficients of the equivalent Legendre
series.
See Also
--------
leg2poly
Notes
-----
The easy way to do conversions between polynomial basis sets
is to use the convert method of a class instance.
Examples
--------
>>> from numpy import polynomial as P
>>> p = P.Polynomial(np.arange(4))
>>> p
Polynomial([0., 1., 2., 3.], domain=[-1., 1.], window=[-1., 1.], ...
>>> c = P.Legendre(P.legendre.poly2leg(p.coef))
>>> c
Legendre([ 1. , 3.25, 1. , 0.75], domain=[-1, 1], window=[-1, 1]) # may vary
"""
[pol] = pu.as_series([pol])
deg = len(pol) - 1
res = 0
for i in range(deg, -1, -1):
res = legadd(legmulx(res), pol[i])
return res
def leg2poly(c):
"""
Convert a Legendre series to a polynomial.
Convert an array representing the coefficients of a Legendre series,
ordered from lowest degree to highest, to an array of the coefficients
of the equivalent polynomial (relative to the "standard" basis) ordered
from lowest to highest degree.
Parameters
----------
c : array_like
1-D array containing the Legendre series coefficients, ordered
from lowest order term to highest.
Returns
-------
pol : ndarray
1-D array containing the coefficients of the equivalent polynomial
(relative to the "standard" basis) ordered from lowest order term
to highest.
See Also
--------
poly2leg
Notes
-----
The easy way to do conversions between polynomial basis sets
is to use the convert method of a class instance.
Examples
--------
>>> from numpy import polynomial as P
>>> c = P.Legendre(range(4))
>>> c
Legendre([0., 1., 2., 3.], domain=[-1., 1.], window=[-1., 1.], symbol='x')
>>> p = c.convert(kind=P.Polynomial)
>>> p
Polynomial([-1. , -3.5, 3. , 7.5], domain=[-1., 1.], window=[-1., ...
>>> P.legendre.leg2poly(range(4))
array([-1. , -3.5, 3. , 7.5])
"""
from .polynomial import polyadd, polysub, polymulx
[c] = pu.as_series([c])
n = len(c)
if n < 3:
return c
else:
c0 = c[-2]
c1 = c[-1]
for i in range(n - 1, 1, -1):
tmp = c0
c0 = polysub(c[i - 2], (c1*(i - 1))/i)
c1 = polyadd(tmp, (polymulx(c1)*(2*i - 1))/i)
return polyadd(c0, polymulx(c1))
legdomain = np.array([-1., 1.])
legzero = np.array([0])
legone = np.array([1])
legx = np.array([0, 1])
def legline(off, scl):
"""
Legendre series whose graph is a straight line.
Parameters
----------
off, scl : scalars
The specified line is given by ``off + scl*x``.
Returns
-------
y : ndarray
This module's representation of the Legendre series for
``off + scl*x``.
See Also
--------
numpy.polynomial.polynomial.polyline
numpy.polynomial.chebyshev.chebline
numpy.polynomial.laguerre.lagline
numpy.polynomial.hermite.hermline
numpy.polynomial.hermite_e.hermeline
Examples
--------
>>> import numpy.polynomial.legendre as L
>>> L.legline(3,2)
array([3, 2])
>>> L.legval(-3, L.legline(3,2)) # should be -3
-3.0
"""
if scl != 0:
return np.array([off, scl])
else:
return np.array([off])
def legfromroots(roots):
"""
"""
Generate a Legendre series with given roots.
The function returns the coefficients of the polynomial
.. math:: p(x) = (x - r_0) * (x - r_1) * ... * (x - r_n),
in Legendre form, where the :math:`r_n` are the roots specified in `roots`.
If a zero has multiplicity n, then it must appear in `roots` n times.
For instance, if 2 is a root of multiplicity three and 3 is a root of
multiplicity 2, then `roots` looks something like [2, 2, 2, 3, 3]. The
roots can appear in any order.
If the returned coefficients are `c`, then
.. math:: p(x) = c_0 + c_1 * L_1(x) + ... + c_n * L_n(x)
The coefficient of the last term is not generally 1 for monic
polynomials in Legendre form.
Parameters
----------
roots : array_like
Sequence containing the roots.
Returns
-------
out : ndarray
1-D array of coefficients. If all roots are real then `out` is a
real array, if some of the roots are complex, then `out` is complex
even if all the coefficients in the result are real (see Examples
below).
See Also
--------
numpy.polynomial.polynomial.polyfromroots
numpy.polynomial.chebyshev.chebfromroots
numpy.polynomial.laguerre.lagfromroots
numpy.polynomial.hermite.hermfromroots
numpy.polynomial.hermite_e.hermefromroots
Examples
--------
>>> import numpy.polynomial.legendre as L
>>> L.legfromroots((-1,0,1))
array([ 0. , -0.4, 0. , 0.4])
>>> j = complex(0,1)
>>> L.legfromroots((-j,j))
array([ 1.33333333+0.j, 0.00000000+0.j, 0.66666667+0.j])
"""
# 使用传入的 roots 参数调用内部函数 _fromroots，生成 Legendre 多项式的系数
return pu._fromroots(legline, legmul, roots)
# 将两个 Legendre 级数相加。
# 返回两个 Legendre 级数 `c1` + `c2` 的和。参数 `c1` 和 `c2` 是按从低到高排序的系数序列，
# 比如，[1,2,3] 表示级数 ``P_0 + 2*P_1 + 3*P_2``。
# 导入必要的包 numpy.polynomial.legendre as L
# c1 和 c2：Legendre 级数系数的 1-D 数组，按从低到高排序。
# 返回
# -------
# out : ndarray
# 表示它们的和的 Legendre 级数的数组。
# See Also
# --------
# legsub, legmulx, legmul, legdiv, legpow
# Notes
# -----
# 与乘法、除法等不同，两个 Legendre 级数的和是一个 Legendre 级数（不需要将结果“重新投影”到基函数集中），
# 因此加法，就像“标准”多项式一样，只是“分量方式”。
# Examples
# --------
# >>> from numpy.polynomial import legendre as L
# >>> c1 = (1,2,3)
# >>> c2 = (3,2,1)
# >>> L.legadd(c1,c2)
# array([4., 4., 4.])
# 将一个 Legendre 级数从另一个中减去。
# 返回两个 Legendre 级数 `c1` - `c2` 的差。系数序列是按从低到高排序的，
# 比如，[1,2,3] 表示级数 ``P_0 + 2*P_1 + 3*P_2``。
# Parameters
# ----------
# c1, c2 : array_like
# 1-D 数组的 Legendre 级数系数，按从低到高排序。
# Returns
# -------
# out : ndarray
# 表示它们的差的 Legendre 级数系数。
# See Also
# --------
# legadd, legmulx, legmul, legdiv, legpow
# Notes
# -----
# 与乘法、除法等不同，两个 Legendre 级数的差是一个 Legendre 级数（不需要将结果“重新投影”到基函数集中），
# 因此减法，就像“标准”多项式一样，只是“分量方式”。
# Examples
# --------
# >>> from numpy.polynomial import legendre as L
# >>> c1 = (1,2,3)
# >>> c2 = (3,2,1)
# >>> L.legsub(c1,c2)
# array([-2., 0., 2.])
# >>> L.legsub(c2,c1) # -C.legsub(c1,c2)
# array([ 2., 0., -2.])
# 将 Legendre 级数乘以 x。
# 将 Legendre 级数 `c` 乘以 x，其中 x 是自变量。
# Parameters
# ----------
# c : array_like
# 按从低到高排序的 1-D 数组的 Legendre 级数系数。
# Returns
# -------
# out : ndarray
# 表示乘法结果的数组。
# See Also
# --------
# legadd, legsub, legmul, legdiv, legpow
# Notes
# -----
# 乘法使用 Legendre 多项式的递归关系形式
# .. math::
# xP_i(x) = ((i + 1)*P_{i + 1}(x) + i*P_{i - 1}(x))/(2i + 1)
# Examples
# --------
# >>> from numpy.polynomial import legendre as L
# >>> L.legmulx([1,2,3])
array([ 0.66666667, 2.2, 1.33333333, 1.8]) # may vary
"""
[c] = pu.as_series([c])
if len(c) == 1 and c[0] == 0:
return c
prd = np.empty(len(c) + 1, dtype=c.dtype)
prd[0] = c[0]*0
prd[1] = c[0]
for i in range(1, len(c)):
j = i + 1
k = i - 1
s = i + j
prd[j] = (c[i]*j)/s
prd[k] += (c[i]*i)/s
return prd
"""
def legdiv(c1, c2):
"""
Divide one Legendre series by another.
Returns the quotient-with-remainder of two Legendre series
`c1` / `c2`. The arguments are sequences of coefficients from lowest
order "term" to highest, e.g., [1,2,3] represents the series
``P_0 + 2*P_1 + 3*P_2``.
Parameters
----------
c1, c2 : array_like
1-D arrays of Legendre series coefficients ordered from low to
high.
Returns
-------
quo, rem : ndarrays
Of Legendre series coefficients representing the quotient and
remainder.
See Also
--------
legadd, legsub, legmulx, legmul, legpow
Notes
-----
In general, the (polynomial) division of one Legendre series by another
results in quotient and remainder terms that are not in the Legendre
polynomial basis set. Thus, to express these results as a Legendre
series, it is necessary to "reproject" the results onto the Legendre
basis set, which may produce "unintuitive" (but correct) results; see
Examples section below.
Examples
--------
>>> from numpy.polynomial import legendre as L
"""
# s1, s2 are trimmed copies
[c1, c2] = pu.as_series([c1, c2])
# Determine the smaller and larger series
if len(c1) > len(c2):
c = c2
xs = c1
else:
c = c1
xs = c2
# Handle division for different lengths of c
if len(c) == 1:
c0 = c[0] / xs
c1 = 0
elif len(c) == 2:
c0 = c[0] / xs
c1 = c[1] / xs
else:
nd = len(c)
c0 = c[-2] / xs
c1 = c[-1] / xs
for i in range(3, len(c) + 1):
tmp = c0
nd = nd - 1
c0 = legsub(c[-i] / xs, (c1 * (nd - 1)) / nd)
c1 = legadd(tmp, (legmulx(c1) * (2 * nd - 1)) / nd)
# Return the result of the division as a Legendre series
return legadd(c0, legmulx(c1))
# 定义两个元组 c1 和 c2，表示为 (1,2,3) 和 (3,2,1)
>>> c1 = (1,2,3)
>>> c2 = (3,2,1)
# 调用 L 对象的 legdiv 方法，对 c1 和 c2 进行 Legendre 多项式的除法运算，返回商和余数
>>> L.legdiv(c1,c2) # quotient "intuitive," remainder not
# 返回结果为 (array([3.]), array([-8., -4.]))，第一个数组是商，第二个数组是余数
(array([3.]), array([-8., -4.]))
# 修改 c2 的值为 (0,1,2,3)，重新调用 legdiv 方法进行除法运算
>>> c2 = (0,1,2,3)
# 返回结果为 (array([-0.07407407, 1.66666667]), array([-1.03703704, -2.51851852]))，可能会有变化
>>> L.legdiv(c2,c1) # neither "intuitive"
# 返回结果分别是商和余数，注意这里的结果可能会有所不同
(array([-0.07407407, 1.66666667]), array([-1.03703704, -2.51851852])) # may vary
"""
return pu._div(legmul, c1, c2)
def legder(c, m=1, scl=1, axis=0):
"""
Differentiate a Legendre series.
Returns the Legendre series coefficients `c` differentiated `m` times
along `axis`. At each iteration the result is multiplied by `scl` (the
scaling factor is for use in a linear change of variable). The argument
`c` is an array of coefficients from low to high degree along each
axis, e.g., [1,2,3] represents the series ``1*L_0 + 2*L_1 + 3*L_2``
while [[1,2],[1,2]] represents ``1*L_0(x)*L_0(y) + 1*L_1(x)*L_0(y) +
2*L_0(x)*L_1(y) + 2*L_1(x)*L_1(y)`` if axis=0 is ``x`` and axis=1 is
``y``.
Parameters
----------
c : array_like
Array of Legendre series coefficients. If c is multidimensional the
different axis correspond to different variables with the degree in
each axis given by the corresponding index.
m : int, optional
Number of derivatives taken, must be non-negative. (Default: 1)
scl : scalar, optional
Each differentiation is multiplied by `scl`. The end result is
multiplication by ``scl**m``. This is for use in a linear change of
variable. (Default: 1)
axis : int, optional
Axis over which the derivative is taken. (Default: 0).
.. versionadded:: 1.7.0
Returns
-------
der : ndarray
Legendre series of the derivative.
See Also
--------
legint
Notes
-----
In general, the result of differentiating a Legendre series does not
resemble the same operation on a power series. Thus the result of this
function may be "unintuitive," albeit correct; see Examples section
below.
Examples
--------
>>> from numpy.polynomial import legendre as L
>>> c = (1,2,3,4)
>>> L.legder(c)
array([ 6., 9., 20.])
>>> L.legder(c, 3)
array([60.])
>>> L.legder(c, scl=-1)
array([ -6., -9., -20.])
>>> L.legder(c, 2,-1)
array([ 9., 60.])
"""
c = np.array(c, ndmin=1, copy=True)
if c.dtype.char in '?bBhHiIlLqQpP':
c = c.astype(np.double)
cnt = pu._as_int(m, "the order of derivation")
iaxis = pu._as_int(axis, "the axis")
if cnt < 0:
raise ValueError("The order of derivation must be non-negative")
iaxis = normalize_axis_index(iaxis, c.ndim)
if cnt == 0:
return c
c = np.moveaxis(c, iaxis, 0)
n = len(c)
if cnt >= n:
c = c[:1]*0
else:
for i in range(cnt):
n = n - 1
c *= scl
der = np.empty((n,) + c.shape[1:], dtype=c.dtype)
for j in range(n, 2, -1):
der[j - 1] = (2*j - 1)*c[j]
c[j - 2] += c[j]
if n > 1:
der[1] = 3*c[2]
der[0] = c[1]
c = der
c = np.moveaxis(c, 0, iaxis)
return c
def legint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
"""
Integrate a Legendre series.
Returns the Legendre series coefficients `c` integrated `m` times from
`lbnd` along `axis`. At each iteration the resulting series is
**multiplied** by `scl` and an integration constant, `k`, is added.
The scaling factor is for use in a linear change of variable. ("Buyer
beware": note that, depending on what one is doing, one may want `scl`
to be the reciprocal of what one might expect; for more information,
see the Notes section below.) The argument `c` is an array of
coefficients from low to high degree along each axis, e.g., [1,2,3]
represents the series ``L_0 + 2*L_1 + 3*L_2`` while [[1,2],[1,2]]
represents ``1*L_0(x)*L_0(y) + 1*L_1(x)*L_0(y) + 2*L_0(x)*L_1(y) +
2*L_1(x)*L_1(y)`` if axis=0 is ``x`` and axis=1 is ``y``.
Parameters
----------
c : array_like
Array of Legendre series coefficients. If c is multidimensional the
different axis correspond to different variables with the degree in
each axis given by the corresponding index.
m : int, optional
Order of integration, must be positive. (Default: 1)
k : {[], list, scalar}, optional
Integration constant(s). The value of the first integral at
``lbnd`` is the first value in the list, the value of the second
integral at ``lbnd`` is the second value, etc. If ``k == []`` (the
default), all constants are set to zero. If ``m == 1``, a single
scalar can be given instead of a list.
lbnd : scalar, optional
The lower bound of the integral. (Default: 0)
scl : scalar, optional
Following each integration the result is *multiplied* by `scl`
before the integration constant is added. (Default: 1)
axis : int, optional
Axis over which the integral is taken. (Default: 0).
.. versionadded:: 1.7.0
Returns
-------
S : ndarray
Legendre series coefficient array of the integral.
Raises
------
ValueError
If ``m < 0``, ``len(k) > m``, ``np.ndim(lbnd) != 0``, or
``np.ndim(scl) != 0``.
See Also
--------
legder
Notes
-----
Note that the result of each integration is *multiplied* by `scl`.
Why is this important to note? Say one is making a linear change of
variable :math:`u = ax + b` in an integral relative to `x`. Then
:math:`dx = du/a`, so one will need to set `scl` equal to
:math:`1/a` - perhaps not what one would have first thought.
Also note that, in general, the result of integrating a C-series needs
to be "reprojected" onto the C-series basis set. Thus, typically,
the result of this function is "unintuitive," albeit correct; see
Examples section below.
Examples
--------
>>> from numpy.polynomial import legendre as L
>>> c = (1,2,3)
>>> L.legint(c)
"""
array([ 0.33333333, 0.4 , 0.66666667, 0.6 ])
>>> L.legint(c, 3)
array([ 1.66666667e-02, -1.78571429e-02, 4.76190476e-02,
-1.73472348e-18, 1.90476190e-02, 9.52380952e-03])
>>> L.legint(c, k=3)
array([ 3.33333333, 0.4 , 0.66666667, 0.6 ])
>>> L.legint(c, lbnd=-2)
array([ 7.33333333, 0.4 , 0.66666667, 0.6 ])
>>> L.legint(c, scl=2)
array([ 0.66666667, 0.8 , 1.33333333, 1.2 ])
"""
# 将 c 转换成至少一维的数组，确保是副本
c = np.array(c, ndmin=1, copy=True)
# 如果 c 是布尔类型或者整数类型，将其转换成双精度浮点型
if c.dtype.char in '?bBhHiIlLqQpP':
c = c.astype(np.double)
# 如果 k 不可迭代，将其转换成列表
if not np.iterable(k):
k = [k]
# 将 m 转换成整数，表示积分的阶数
cnt = pu._as_int(m, "the order of integration")
# 将 axis 转换成整数，表示坐标轴
iaxis = pu._as_int(axis, "the axis")
# 如果积分阶数为负数，抛出异常
if cnt < 0:
raise ValueError("The order of integration must be non-negative")
# 如果 k 的长度大于积分阶数，抛出异常
if len(k) > cnt:
raise ValueError("Too many integration constants")
# 如果 lbnd 不是标量，抛出异常
if np.ndim(lbnd) != 0:
raise ValueError("lbnd must be a scalar.")
# 如果 scl 不是标量，抛出异常
if np.ndim(scl) != 0:
raise ValueError("scl must be a scalar.")
# 根据 c 的维度和 axis 规范化 iaxis
iaxis = normalize_axis_index(iaxis, c.ndim)
# 如果阶数为 0，直接返回 c
if cnt == 0:
return c
# 将 c 沿着轴移动到第一个位置
c = np.moveaxis(c, iaxis, 0)
# 将 k 转换成列表，并在末尾补充 0 至 cnt 长度
k = list(k) + [0]*(cnt - len(k))
for i in range(cnt):
n = len(c)
c *= scl
if n == 1 and np.all(c[0] == 0):
c[0] += k[i]
else:
tmp = np.empty((n + 1,) + c.shape[1:], dtype=c.dtype)
tmp[0] = c[0]*0
tmp[1] = c[0]
if n > 1:
tmp[2] = c[1]/3
for j in range(2, n):
t = c[j]/(2*j + 1)
tmp[j + 1] = t
tmp[j - 1] -= t
tmp[0] += k[i] - legval(lbnd, tmp)
c = tmp
# 将 c 沿着第一个轴移动到 iaxis 轴上
c = np.moveaxis(c, 0, iaxis)
# 返回 c
return c
# 定义函数，用于在给定点 x 处评估 Legendre 级数
def legval(x, c, tensor=True):
"""
Evaluate a Legendre series at points x.
If `c` is of length ``n + 1``, this function returns the value:
.. math:: p(x) = c_0 * L_0(x) + c_1 * L_1(x) + ... + c_n * L_n(x)
The parameter `x` is converted to an array only if it is a tuple or a
list, otherwise it is treated as a scalar. In either case, either `x`
or its elements must support multiplication and addition both with
themselves and with the elements of `c`.
If `c` is a 1-D array, then ``p(x)`` will have the same shape as `x`. If
`c` is multidimensional, then the shape of the result depends on the
value of `tensor`. If `tensor` is true the shape will be c.shape[1:] +
x.shape. If `tensor` is false the shape will be c.shape[1:]. Note that
scalars have shape (,).
Trailing zeros in the coefficients will be used in the evaluation, so
they should be avoided if efficiency is a concern.
Parameters
----------
x : array_like, compatible object
If `x` is a list or tuple, it is converted to an ndarray, otherwise
it is left unchanged and treated as a scalar. In either case, `x`
or its elements must support addition and multiplication with
themselves and with the elements of `c`.
c : array_like
Array of coefficients ordered so that the coefficients for terms of
degree n are contained in c[n]. If `c` is multidimensional the
remaining indices enumerate multiple polynomials. In the two
dimensional case the coefficients may be thought of as stored in
the columns of `c`.
tensor : boolean, optional
If True, the shape of the coefficient array is extended with ones
on the right, one for each dimension of `x`. Scalars have dimension 0
for this action. The result is that every column of coefficients in
`c` is evaluated for every element of `x`. If False, `x` is broadcast
over the columns of `c` for the evaluation. This keyword is useful
when `c` is multidimensional. The default value is True.
.. versionadded:: 1.7.0
Returns
-------
values : ndarray, algebra_like
The shape of the return value is described above.
See Also
--------
legval2d, leggrid2d, legval3d, leggrid3d
Notes
-----
The evaluation uses Clenshaw recursion, aka synthetic division.
"""
# 将 c 转换为至少是一维的 numpy 数组，确保数据类型为浮点型
c = np.array(c, ndmin=1, copy=None)
# 如果 c 的数据类型为布尔型、字节型、短整型、整型或长整型，转换为双精度浮点型
if c.dtype.char in '?bBhHiIlLqQpP':
c = c.astype(np.double)
# 如果 x 是元组或列表，则将其转换为 ndarray
if isinstance(x, (tuple, list)):
x = np.asarray(x)
# 如果 x 是 ndarray 且 tensor=True，则将 c 的形状重塑为 c.shape + (1,)*x.ndim
if isinstance(x, np.ndarray) and tensor:
c = c.reshape(c.shape + (1,)*x.ndim)
# 根据 c 的长度选择性地设置 c0 和 c1
if len(c) == 1:
c0 = c[0]
c1 = 0
elif len(c) == 2:
c0 = c[0]
c1 = c[1]
else:
# 计算输入列表 c 的长度
nd = len(c)
# 获取 c 的倒数第二个元素
c0 = c[-2]
# 获取 c 的最后一个元素
c1 = c[-1]
# 循环遍历 c 的子集，从索引 3 开始到末尾
for i in range(3, len(c) + 1):
# 临时保存 c0 的值
tmp = c0
# 更新 nd 的值
nd = nd - 1
# 更新 c0 的值，根据公式计算
c0 = c[-i] - (c1*(nd - 1))/nd
# 更新 c1 的值，根据公式计算
c1 = tmp + (c1*x*(2*nd - 1))/nd
# 返回计算结果 c0 + c1*x
return c0 + c1*x
# 计算二维Legendre级数在点(x, y)处的值。
def legval2d(x, y, c):
"""
Evaluate a 2-D Legendre series at points (x, y).
This function returns the values:
.. math:: p(x,y) = \\sum_{i,j} c_{i,j} * L_i(x) * L_j(y)
The parameters `x` and `y` are converted to arrays only if they are
tuples or a lists, otherwise they are treated as scalars and they
must have the same shape after conversion. In either case, either `x`
and `y` or their elements must support multiplication and addition both
with themselves and with the elements of `c`.
If `c` is a 1-D array a one is implicitly appended to its shape to make
it 2-D. The shape of the result will be c.shape[2:] + x.shape.
Parameters
----------
x, y : array_like, compatible objects
The two-dimensional series is evaluated at the points ``(x, y)``,
where `x` and `y` must have the same shape. If `x` or `y` is a list
or tuple, it is first converted to an ndarray, otherwise it is left
unchanged and if it isn't an ndarray it is treated as a scalar.
c : array_like
Array of coefficients ordered so that the coefficient of the term
of multi-degree i,j is contained in ``c[i,j]``. If `c` has
dimension greater than two the remaining indices enumerate multiple
sets of coefficients.
Returns
-------
values : ndarray, compatible object
The values of the two-dimensional Legendre series at points formed
from pairs of corresponding values from `x` and `y`.
See Also
--------
legval, leggrid2d, legval3d, leggrid3d
Notes
-----
.. versionadded:: 1.7.0
"""
# 调用私有函数 `_valnd` 计算Legendre级数在给定点 `(x, y)` 处的值，使用系数 `c`
return pu._valnd(legval, c, x, y)
def leggrid2d(x, y, c):
"""
Evaluate a 2-D Legendre series on the Cartesian product of x and y.
This function returns the values:
.. math:: p(a,b) = \\sum_{i,j} c_{i,j} * L_i(a) * L_j(b)
where the points ``(a, b)`` consist of all pairs formed by taking
`a` from `x` and `b` from `y`. The resulting points form a grid with
`x` in the first dimension and `y` in the second.
The parameters `x` and `y` are converted to arrays only if they are
tuples or a lists, otherwise they are treated as scalars. In either
case, either `x` and `y` or their elements must support multiplication
and addition both with themselves and with the elements of `c`.
If `c` has fewer than two dimensions, ones are implicitly appended to
its shape to make it 2-D. The shape of the result will be c.shape[2:] +
x.shape + y.shape.
Parameters
----------
x, y : array_like, compatible objects
The two-dimensional series is evaluated at the points in the
Cartesian product of `x` and `y`. If `x` or `y` is a list or
tuple, it is first converted to an ndarray, otherwise it is left
unchanged and, if it isn't an ndarray, it is treated as a scalar.
c : array_like
Array of coefficients ordered so that the coefficient of the term
of multi-degree i,j is contained in ``c[i,j]``. If `c` has
dimension greater than two the remaining indices enumerate multiple
sets of coefficients.
Returns
-------
values : ndarray, compatible object
The values of the two-dimensional Legendre series at points formed
from pairs of corresponding values from `x` and `y`.
See Also
--------
legval, leggrid2d, legval3d, leggrid3d
"""
c : array_like
# 输入参数 c 是一个类数组对象，用于存储系数，按照多重度 i,j 的项的系数存储在 c[i,j] 中。
# 如果 c 的维度大于两维，则剩余的索引用于枚举多组系数。
Returns
-------
values : ndarray, compatible object
# 返回值是一个 ndarray 数组或兼容对象，表示在笛卡尔积空间中点 (x, y) 处的二维切比雪夫级数的值。
See Also
--------
legval, legval2d, legval3d, leggrid3d
# 参见相关函数 legval, legval2d, legval3d, leggrid3d。
Notes
-----
# 版本备注：自版本 1.7.0 起添加。
"""
return pu._gridnd(legval, c, x, y)
def legval3d(x, y, z, c):
"""
Evaluate a 3-D Legendre series at points (x, y, z).
This function returns the values:
.. math:: p(x,y,z) = \\sum_{i,j,k} c_{i,j,k} * L_i(x) * L_j(y) * L_k(z)
The parameters `x`, `y`, and `z` are converted to arrays only if
they are tuples or a lists, otherwise they are treated as scalars and
must have the same shape after conversion. In either case, either
`x`, `y`, and `z` or their elements must support multiplication and
addition both with themselves and with the elements of `c`.
If `c` has fewer than 3 dimensions, ones are implicitly appended to its
shape to make it 3-D. The shape of the result will be c.shape[3:] +
x.shape.
Parameters
----------
x, y, z : array_like, compatible object
The three-dimensional series is evaluated at the points
``(x, y, z)``, where `x`, `y`, and `z` must have the same shape. If
any of `x`, `y`, or `z` is a list or tuple, it is first converted
to an ndarray; otherwise, it is treated as a scalar.
c : array_like
Array of coefficients ordered so that the coefficient of the term of
multi-degree i,j,k is contained in ``c[i,j,k]``. If `c` has dimension
greater than 3, the remaining indices enumerate multiple sets of
coefficients.
Returns
-------
values : ndarray, compatible object
The values of the multidimensional polynomial on points formed with
triples of corresponding values from `x`, `y`, and `z`.
See Also
--------
legval, legval2d, leggrid2d, leggrid3d
Notes
-----
.. versionadded:: 1.7.0
"""
return pu._valnd(legval, c, x, y, z)
def leggrid3d(x, y, z, c):
"""
Evaluate a 3-D Legendre series on the Cartesian product of x, y, and z.
This function returns the values:
.. math:: p(a,b,c) = \\sum_{i,j,k} c_{i,j,k} * L_i(a) * L_j(b) * L_k(c)
where the points ``(a, b, c)`` consist of all triples formed by taking
`a` from `x`, `b` from `y`, and `c` from `z`. The resulting points form
a grid with `x` in the first dimension, `y` in the second, and `z` in
the third.
The parameters `x`, `y`, and `z` are converted to arrays only if they
are tuples or lists; otherwise, they are treated as scalars. In
either case, either `x`, `y`, and `z` or their elements must support
multiplication and addition both with themselves and with the elements
of `c`.
If `c` has fewer than three dimensions, ones are implicitly appended to
its shape to make it 3-D. The shape of the result will be c.shape[3:] +
x.shape + y.shape + z.shape.
Parameters
----------
x, y, z : array_like, compatible object
Arrays defining the grid points where the Legendre series is evaluated.
These arrays should have compatible shapes.
c : array_like
Array of coefficients ordered such that the coefficient of the term of
multi-degree i,j,k is contained in ``c[i,j,k]``. If `c` has dimension
greater than 3, the remaining indices enumerate multiple sets of
coefficients.
Returns
-------
values : ndarray, compatible object
The values of the multidimensional polynomial on the Cartesian product
of `x`, `y`, and `z`.
See Also
--------
legval, legval2d, leggrid2d, leggrid3d
Notes
-----
.. versionadded:: 1.7.0
"""
x, y, z : array_like, compatible objects
三维序列在 `x`、`y` 和 `z` 的笛卡尔积点处进行评估。如果 `x`、`y` 或 `z` 是列表或元组，则首先转换为 ndarray；否则保持不变，并且如果它不是 ndarray，则视为标量。
c : array_like
按照次数 i,j 排列的系数数组，系数 `c[i,j]` 对应于二维多项式的系数。如果 `c` 的维度大于二，则其余的索引用于枚举多组系数。
Returns
-------
values : ndarray, compatible object
二维多项式在 `x` 和 `y` 的笛卡尔积点处的值。
See Also
--------
legval, legval2d, leggrid2d, legval3d
Notes
-----
在版本 1.7.0 中添加
"""
返回 pu._gridnd(legval, c, x, y, z)
注释：
- `x, y, z : array_like, compatible objects`: 描述了函数参数 `x`, `y`, `z` 应为类似数组的兼容对象，表示三维序列在这些坐标的笛卡尔积点处进行评估。
- `c : array_like`: 描述了参数 `c` 应为类似数组的对象，按照二维多项式的次数排列，用于存储系数。
- `Returns`: 描述了函数返回一个 ndarray 或兼容对象，表示二维多项式在给定点处的值。
- `See Also`: 提供了相关函数的链接，供进一步参考。
- `Notes`: 指出该功能是在版本 1.7.0 中添加的。
def legvander(x, deg):
"""Pseudo-Vandermonde matrix of given degree.
Returns the pseudo-Vandermonde matrix of degree `deg` and sample points
`x`. The pseudo-Vandermonde matrix is defined by
.. math:: V[..., i] = L_i(x)
where ``0 <= i <= deg``. The leading indices of `V` index the elements of
`x` and the last index is the degree of the Legendre polynomial.
If `c` is a 1-D array of coefficients of length ``n + 1`` and `V` is the
array ``V = legvander(x, n)``, then ``np.dot(V, c)`` and
``legval(x, c)`` are the same up to roundoff. This equivalence is
useful both for least squares fitting and for the evaluation of a large
number of Legendre series of the same degree and sample points.
Parameters
----------
x : array_like
Array of points. The dtype is converted to float64 or complex128
depending on whether any of the elements are complex. If `x` is
scalar it is converted to a 1-D array.
deg : int
Degree of the resulting matrix.
Returns
-------
vander : ndarray
The pseudo-Vandermonde matrix. The shape of the returned matrix is
``x.shape + (deg + 1,)``, where The last index is the degree of the
corresponding Legendre polynomial. The dtype will be the same as
the converted `x`.
"""
# Convert deg to an integer if possible, raises ValueError if deg < 0
ideg = pu._as_int(deg, "deg")
if ideg < 0:
raise ValueError("deg must be non-negative")
# Ensure x is converted to a 1-D array of float64 or complex128
x = np.array(x, copy=None, ndmin=1) + 0.0
# Calculate dimensions for the output array
dims = (ideg + 1,) + x.shape
dtyp = x.dtype
# Create an uninitialized array v with dimensions and dtype calculated above
v = np.empty(dims, dtype=dtyp)
# Initialize the first column of v with ones
v[0] = x*0 + 1
# Fill the remaining columns using forward recursion
if ideg > 0:
v[1] = x
for i in range(2, ideg + 1):
v[i] = (v[i-1]*x*(2*i - 1) - v[i-2]*(i - 1))/i
# Move the first axis to the last axis of v and return
return np.moveaxis(v, 0, -1)
x, y : array_like
Arrays of point coordinates, all of the same shape. The dtypes
will be converted to either float64 or complex128 depending on
whether any of the elements are complex. Scalars are converted to
1-D arrays.
deg : list of ints
List of maximum degrees of the form [x_deg, y_deg].
Returns
-------
vander2d : ndarray
The shape of the returned matrix is ``x.shape + (order,)``, where
:math:`order = (deg[0]+1)*(deg[1]+1)`. The dtype will be the same
as the converted `x` and `y`.
See Also
--------
legvander, legvander3d, legval2d, legval3d
Notes
-----
.. versionadded:: 1.7.0
"""
return pu._vander_nd_flat((legvander, legvander), (x, y), deg)
def legvander3d(x, y, z, deg):
"""Pseudo-Vandermonde matrix of given degrees.
Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
points ``(x, y, z)``. If `l`, `m`, `n` are the given degrees in `x`, `y`, `z`,
then The pseudo-Vandermonde matrix is defined by
.. math:: V[..., (m+1)(n+1)i + (n+1)j + k] = L_i(x)*L_j(y)*L_k(z),
where ``0 <= i <= l``, ``0 <= j <= m``, and ``0 <= j <= n``. The leading
indices of `V` index the points ``(x, y, z)`` and the last index encodes
the degrees of the Legendre polynomials.
If ``V = legvander3d(x, y, z, [xdeg, ydeg, zdeg])``, then the columns
of `V` correspond to the elements of a 3-D coefficient array `c` of
shape (xdeg + 1, ydeg + 1, zdeg + 1) in the order
.. math:: c_{000}, c_{001}, c_{002},... , c_{010}, c_{011}, c_{012},...
and ``np.dot(V, c.flat)`` and ``legval3d(x, y, z, c)`` will be the
same up to roundoff. This equivalence is useful both for least squares
fitting and for the evaluation of a large number of 3-D Legendre
series of the same degrees and sample points.
Parameters
----------
x, y, z : array_like
Arrays of point coordinates, all of the same shape. The dtypes will
be converted to either float64 or complex128 depending on whether
any of the elements are complex. Scalars are converted to 1-D
arrays.
deg : list of ints
List of maximum degrees of the form [x_deg, y_deg, z_deg].
Returns
-------
vander3d : ndarray
The shape of the returned matrix is ``x.shape + (order,)``, where
:math:`order = (deg[0]+1)*(deg[1]+1)*(deg[2]+1)`. The dtype will
be the same as the converted `x`, `y`, and `z`.
See Also
--------
legvander, legval2d, legval3d
Notes
-----
.. versionadded:: 1.7.0
"""
return pu._vander_nd_flat((legvander, legvander, legvander), (x, y, z), deg)
def legfit(x, y, deg, rcond=None, full=False, w=None):
"""
Least squares fit of Legendre series to data.
Return the coefficients of a Legendre series of degree `deg` that is the
least squares fit to the data values `y` given at points `x`. If `y` is
1-D the returned coefficients will also be 1-D. If `y` is 2-D multiple
fits are done, one for each column of `y`, and the resulting
coefficients are stored in the corresponding columns of a 2-D return.
The fitted polynomial(s) are in the form
.. math:: p(x) = c_0 + c_1 * L_1(x) + ... + c_n * L_n(x),
where `n` is `deg`.
Parameters
----------
x : array_like, shape (M,)
x-coordinates of the M sample points ``(x[i], y[i])``.
y : array_like, shape (M,) or (M, K)
y-coordinates of the sample points. Several data sets of sample
points sharing the same x-coordinates can be fitted at once by
passing in a 2D-array that contains one dataset per column.
deg : int
Degree of the Legendre series to fit.
Returns
-------
c : ndarray, shape (deg + 1,) or (deg + 1, K)
Coefficients of the Legendre series that minimizes the squared
error in the least-squares sense.
See Also
--------
legvander, legvander2d, legvander3d, legval, legval2d, legval3d
"""
return pu._fit((legvander, legvander, legvander), x, y, deg, rcond, full, w)
deg : int or 1-D array_like
rcond : float, optional
full : bool, optional
w : array_like, shape (`M`,), optional
.. versionadded:: 1.5.0
Returns
-------
coef : ndarray, shape (M,) or (M, K)
[residuals, rank, singular_values, rcond] : list
- residuals -- 最小二乘拟合的残差平方和
- rank -- 缩放后的 Vandermonde 矩阵的数值秩
- singular_values -- 缩放后的 Vandermonde 矩阵的奇异值
- rcond -- `rcond` 的值
更多细节请参阅 `numpy.linalg.lstsq`。
Warns
-----
RankWarning
>>> import warnings
>>> warnings.simplefilter('ignore', np.exceptions.RankWarning)
See Also
--------
numpy.polynomial.polynomial.polyfit
numpy.polynomial.chebyshev.chebfit
numpy.polynomial.laguerre.lagfit
numpy.polynomial.hermite.hermfit
numpy.polynomial.hermite_e.hermefit
legval : 对 Legendre 级数进行求值。
legvander : Legendre 级数的 Vandermonde 矩阵。
legweight : Legendre 权重函数（= 1）。
numpy.linalg.lstsq : 计算矩阵的最小二乘拟合。
scipy.interpolate.UnivariateSpline : 计算样条拟合。
Notes
-----
minimizes the sum of the weighted squared errors
.. math:: E = \\sum_j w_j^2 * |y_j - p(x_j)|^2,
where :math:`w_j` are the weights. This problem is solved by setting up
as the (typically) overdetermined matrix equation
.. math:: V(x) * c = w * y,
where `V` is the weighted pseudo Vandermonde matrix of `x`, `c` are the
coefficients to be solved for, `w` are the weights, and `y` are the
observed values. This equation is then solved using the singular value
decomposition of `V`.
If some of the singular values of `V` are so small that they are
neglected, then a `~exceptions.RankWarning` will be issued. This means that
the coefficient values may be poorly determined. Using a lower order fit
will usually get rid of the warning. The `rcond` parameter can also be
set to a value smaller than its default, but the resulting fit may be
spurious and have large contributions from roundoff error.
Fits using Legendre series are usually better conditioned than fits
using power series, but much can depend on the distribution of the
sample points and the smoothness of the data. If the quality of the fit
is inadequate splines may be a good alternative.
References
----------
.. [1] Wikipedia, "Curve fitting",
https://en.wikipedia.org/wiki/Curve_fitting
Examples
--------
"""
使用Legendre多项式进行最小二乘拟合，并返回拟合结果。
Parameters
----------
legvander : callable
函数用于计算Legendre多项式的伪Vandermonde矩阵。
x : array_like, shape (M,)
输入数据点的数组。
y : array_like, shape (M,)
观测到的数据值。
deg : int
拟合多项式的次数。
rcond : float, optional
用于判断奇异值的截断值，默认为None。
full : bool, optional
控制输出内容，若为True则返回完整的分解，否则返回拟合系数，默认为False。
w : array_like, shape (M,), optional
每个观测值的权重。
Returns
-------
ndarray
拟合系数。
Notes
-----
使用Legendre多项式进行拟合通常比使用幂级数更好，但很大程度上取决于样本点的分布和数据的平滑程度。
如果拟合质量不足，样条曲线可能是一个更好的选择。
"""
return pu._fit(legvander, x, y, deg, rcond, full, w)
def legcompanion(c):
"""
Return the scaled companion matrix of c.
The basis polynomials are scaled so that the companion matrix is
symmetric when `c` is an Legendre basis polynomial. This provides
better eigenvalue estimates than the unscaled case and for basis
polynomials the eigenvalues are guaranteed to be real if
`numpy.linalg.eigvalsh` is used to obtain them.
Parameters
----------
c : array_like
1-D array of Legendre series coefficients ordered from low to high
degree.
Returns
-------
mat : ndarray
Scaled companion matrix of dimensions (deg, deg).
Notes
-----
.. versionadded:: 1.7.0
"""
[c] = pu.as_series([c])
if len(c) < 2:
raise ValueError('Series must have maximum degree of at least 1.')
if len(c) == 2:
return np.array([[-c[0]/c[1]]])
n = len(c) - 1
mat = np.zeros((n, n), dtype=c.dtype)
scl = 1./np.sqrt(2*np.arange(n) + 1)
top = mat.reshape(-1)[1::n+1]
bot = mat.reshape(-1)[n::n+1]
top[...] = np.arange(1, n)*scl[:n-1]*scl[1:n]
bot[...] = top
mat[:, -1] -= (c[:-1]/c[-1])*(scl/scl[-1])*(n/(2*n - 1))
return mat
def legroots(c):
"""
Compute the roots of a Legendre series.
Return the roots (a.k.a. "zeros") of the polynomial
.. math:: p(x) = \\sum_i c[i] * L_i(x).
Parameters
----------
c : 1-D array_like
1-D array of coefficients.
Returns
-------
out : ndarray
Array of the roots of the series. If all the roots are real,
then `out` is also real, otherwise it is complex.
See Also
--------
numpy.polynomial.polynomial.polyroots
numpy.polynomial.chebyshev.chebroots
numpy.polynomial.laguerre.lagroots
numpy.polynomial.hermite.hermroots
numpy.polynomial.hermite_e.hermeroots
Notes
-----
The root estimates are obtained as the eigenvalues of the companion
matrix, Roots far from the origin of the complex plane may have large
errors due to the numerical instability of the series for such values.
Roots with multiplicity greater than 1 will also show larger errors as
the value of the series near such points is relatively insensitive to
errors in the roots. Isolated roots near the origin can be improved by
a few iterations of Newton's method.
The Legendre series basis polynomials aren't powers of ``x`` so the
results of this function may seem unintuitive.
Examples
--------
>>> import numpy.polynomial.legendre as leg
>>> leg.legroots((1, 2, 3, 4)) # 4L_3 + 3L_2 + 2L_1 + 1L_0, all real roots
array([-0.85099543, -0.11407192, 0.51506735]) # may vary
"""
[c] = pu.as_series([c])
if len(c) < 2:
return np.array([], dtype=c.dtype)
if len(c) == 2:
return np.array([-c[0]/c[1]])
m = legcompanion(c)[::-1,::-1]
r = la.eigvals(m)
r.sort()
return r
def leggauss(deg):
ideg = pu._as_int(deg, "deg")
if ideg <= 0:
raise ValueError("deg must be a positive integer")
c = np.array([0]*deg + [1])
m = legcompanion(c)
x = la.eigvalsh(m)
dy = legval(x, c)
df = legval(x, legder(c))
x -= dy / df
fm = legval(x, c[1:])
fm /= np.abs(fm).max()
df /= np.abs(df).max()
w = 1 / (fm * df)
w = (w + w[::-1]) / 2
x = (x - x[::-1]) / 2
w *= 2. / w.sum()
return x, w
def legweight(x):
w = x * 0.0 + 1.0
return w
class Legendre(ABCPolyBase):
"""A Legendre series class.
The Legendre class provides the standard Python numerical methods
'+', '-', '*', '//', '%', 'divmod', '**', and '()' as well as the
attributes and methods listed below.
Parameters
----------
coef : array_like
Legendre coefficients in order of increasing degree, i.e.,
``(1, 2, 3)`` gives ``1*P_0(x) + 2*P_1(x) + 3*P_2(x)``.
"""
"""
# Virtual Functions
# 定义静态方法来实现多项式的不同操作
_add = staticmethod(legadd)
_sub = staticmethod(legsub)
_mul = staticmethod(legmul)
_div = staticmethod(legdiv)
_pow = staticmethod(legpow)
_val = staticmethod(legval)
_int = staticmethod(legint)
_der = staticmethod(legder)
_fit = staticmethod(legfit)
_line = staticmethod(legline)
_roots = staticmethod(legroots)
_fromroots = staticmethod(legfromroots)
# Virtual properties
# 将定义的区间数组转换为NumPy数组并赋值给domain和window属性
domain = np.array(legdomain)
window = np.array(legdomain)
# 设置多项式的基函数名称为'P'
basis_name = 'P'
.\numpy\numpy\polynomial\legendre.pyi
from typing import Any
from numpy import int_
from numpy.typing import NDArray
from numpy.polynomial._polybase import ABCPolyBase
from numpy.polynomial.polyutils import trimcoef
__all__: list[str]
legtrim = trimcoef
def poly2leg(pol): ...
def leg2poly(c): ...
legdomain: NDArray[int_]
legzero: NDArray[int_]
legone: NDArray[int_]
legx: NDArray[int_]
def legline(off, scl): ...
def legfromroots(roots): ...
def legadd(c1, c2): ...
def legsub(c1, c2): ...
def legmulx(c): ...
def legmul(c1, c2): ...
def legdiv(c1, c2): ...
def legpow(c, pow, maxpower=...): ...
def legder(c, m=..., scl=..., axis=...): ...
def legint(c, m=..., k = ..., lbnd=..., scl=..., axis=...): ...
def legval(x, c, tensor=...): ...
def legval2d(x, y, c): ...
def leggrid2d(x, y, c): ...
def legval3d(x, y, z, c): ...
def leggrid3d(x, y, z, c): ...
def legvander(x, deg): ...
def legvander2d(x, y, deg): ...
def legvander3d(x, y, z, deg): ...
def legfit(x, y, deg, rcond=..., full=..., w=...): ...
def legcompanion(c): ...
def legroots(c): ...
def leggauss(deg): ...
def legweight(x): ...
class Legendre(ABCPolyBase):
domain: Any
window: Any
basis_name: Any
.\numpy\numpy\polynomial\polynomial.py
__all__ = [
'polyzero', 'polyone', 'polyx', 'polydomain', 'polyline', 'polyadd',
'polysub', 'polymulx', 'polymul', 'polydiv', 'polypow', 'polyval',
'polyvalfromroots', 'polyder', 'polyint', 'polyfromroots', 'polyvander',
'polyfit', 'polytrim', 'polyroots', 'Polynomial', 'polyval2d', 'polyval3d',
'polygrid2d', 'polygrid3d', 'polyvander2d', 'polyvander3d',
'polycompanion']
import numpy as np
import numpy.linalg as la
from numpy.lib.array_utils import normalize_axis_index
from . import polyutils as pu
from ._polybase import ABCPolyBase
polytrim = pu.trimcoef
polydomain = np.array([-1., 1.])
polyzero = np.array([0])
polyone = np.array([1])
polyx = np.array([0, 1])
def polyline(off, scl):
"""
返回表示线性多项式的数组。
Parameters
----------
off, scl : scalars
线性多项式的截距和斜率。
Returns
-------
y : ndarray
表示线性多项式 `off + scl*x` 的数组。
See Also
--------
numpy.polynomial.chebyshev.chebline
numpy.polynomial.legendre.legline
numpy.polynomial.laguerre.lagline
numpy.polynomial.hermite.hermline
numpy.polynomial.hermite_e.hermeline
Examples
--------
>>> from numpy.polynomial import polynomial as P
>>> P.polyline(1, -1)
array([ 1, -1])
>>> P.polyval(1, P.polyline(1, -1)) # 应该为 0
"""
return np.array([off, scl])
if scl != 0:
return np.array([off, scl])
else:
return np.array([off])
def polyfromroots(roots):
return pu._fromroots(polyline, polymul, roots)
def polyadd(c1, c2):
return pu._add(c1, c2)
def polysub(c1, c2):
def polysub(c1, c2):
return pu._sub(c1, c2)
def polymulx(c):
"""Multiply a polynomial by x.
Multiply the polynomial `c` by x, where x is the independent
variable.
Parameters
----------
c : array_like
1-D array of polynomial coefficients ordered from low to
high.
Returns
-------
out : ndarray
Array representing the result of the multiplication.
See Also
--------
polyadd, polysub, polymul, polydiv, polypow
Notes
-----
.. versionadded:: 1.5.0
Examples
--------
>>> from numpy.polynomial import polynomial as P
>>> c = (1, 2, 3)
>>> P.polymulx(c)
array([0., 1., 2., 3.])
"""
[c] = pu.as_series([c])
if len(c) == 1 and c[0] == 0:
return c
prd = np.empty(len(c) + 1, dtype=c.dtype)
prd[0] = c[0]*0
prd[1:] = c
return prd
def polymul(c1, c2):
"""
Multiply one polynomial by another.
Returns the product of two polynomials `c1` * `c2`. The arguments are
sequences of coefficients, from lowest order term to highest, e.g.,
[1,2,3] represents the polynomial ``1 + 2*x + 3*x**2.``
Parameters
----------
c1, c2 : array_like
1-D arrays of coefficients representing a polynomial, relative to the
"standard" basis, and ordered from lowest order term to highest.
Returns
-------
out : ndarray
Of the coefficients of their product.
See Also
--------
polyadd, polysub, polymulx, polydiv, polypow
Examples
--------
>>> from numpy.polynomial import polynomial as P
>>> c1 = (1, 2, 3)
>>> c2 = (3, 2, 1)
>>> P.polymul(c1, c2)
array([ 3., 8., 14., 8., 3.])
"""
[c1, c2] = pu.as_series([c1, c2])
ret = np.convolve(c1, c2)
return pu.trimseq(ret)
def polydiv(c1, c2):
"""
Divide one polynomial by another.
Returns the quotient-with-remainder of two polynomials `c1` / `c2`.
The arguments are sequences of coefficients, from lowest order term
to highest, e.g., [1,2,3] represents ``1 + 2*x + 3*x**2``.
Parameters
----------
c1, c2 : array_like
1-D arrays of polynomial coefficients ordered from low to high.
Returns
-------
[quo, rem] : ndarrays
Of coefficient series representing the quotient and remainder.
See Also
--------
polyadd, polysub, polymulx, polymul, polypow
Examples
--------
>>> from numpy.polynomial import polynomial as P
>>> c1 = (1, 2, 3)
>>> c2 = (3, 2, 1)
>>> P.polydiv(c1, c2)
(array([3.]), array([-8., -4.]))
>>> P.polydiv(c2, c1)
(array([ 0.33333333]), array([ 2.66666667, 1.33333333])) # may vary
"""
[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:
dlen = lc1 - lc2
scl = c2[-1]
c2 = c2[:-1]/scl
i = dlen
j = lc1 - 1
while i >= 0:
c1[i:j] -= c2*c1[j]
i -= 1
j -= 1
return c1[j+1:]/scl, pu.trimseq(c1[:j+1])
def polypow(c, pow, maxpower=None):
"""Raise a polynomial to a power.
Returns the polynomial `c` raised to the power `pow`. The argument
`c` is a sequence of coefficients ordered from low to high. i.e.,
[1,2,3] is the series ``1 + 2*x + 3*x**2.``
Parameters
----------
c : array_like
1-D array of array of series coefficients ordered from low to
high degree.
pow : integer
Power to which the series will be raised
maxpower : integer, optional
Maximum power allowed. This is mainly to limit growth of the series
to unmanageable size. Default is 16
Returns
-------
coef : ndarray
Power series of power.
See Also
--------
polyadd, polysub, polymulx, polymul, polydiv
Examples
--------
>>> from numpy.polynomial import polynomial as P
>>> P.polypow([1, 2, 3], 2)
array([ 1., 4., 10., 12., 9.])
"""
return pu._pow(np.convolve, c, pow, maxpower)
def polyder(c, m=1, scl=1, axis=0):
"""
Differentiate a polynomial.
Returns the polynomial coefficients `c` differentiated `m` times along
`axis`. At each iteration the result is multiplied by `scl` (the
scaling factor is for use in a linear change of variable). The
argument `c` is an array of coefficients from low to high degree along
each axis, e.g., [1,2,3] represents the polynomial ``1 + 2*x + 3*x**2``
while [[1,2],[1,2]] represents ``1 + 1*x + 2*y + 2*x*y`` if axis=0 is
``x`` and axis=1 is ``y``.
Parameters
----------
c : array_like
Array of polynomial coefficients. If c is multidimensional the
different axis correspond to different variables with the degree
in each axis given by the corresponding index.
m : int, optional
Number of derivatives taken, must be non-negative. (Default: 1)
scl : scalar, optional
Each differentiation is multiplied by `scl`. The end result is
multiplication by ``scl**m``. This is for use in a linear change
of variable. (Default: 1)
axis : int, optional
Axis over which the derivative is taken. (Default: 0).
.. versionadded:: 1.7.0
Returns
-------
der : ndarray
Polynomial coefficients of the derivative.
See Also
--------
polyint
Examples
--------
>>> from numpy.polynomial import polynomial as P
>>> c = (1, 2, 3, 4)
>>> P.polyder(c) # (d/dx)(c)
array([ 2., 6., 12.])
>>> P.polyder(c, 3) # (d**3/dx**3)(c)
array([24.])
>>> P.polyder(c, scl=-1) # (d/d(-x))(c)
array([ -2., -6., -12.])
>>> P.polyder(c, 2, -1) # (d**2/d(-x)**2)(c)
array([ 6., 24.])
"""
c = np.array(c, ndmin=1, copy=True)
if c.dtype.char in '?bBhHiIlLqQpP':
c = c + 0.0
cdt = c.dtype
cnt = pu._as_int(m, "the order of derivation")
iaxis = pu._as_int(axis, "the axis")
if cnt < 0:
raise ValueError("The order of derivation must be non-negative")
iaxis = normalize_axis_index(iaxis, c.ndim)
if cnt == 0:
return c
c = np.moveaxis(c, iaxis, 0)
n = len(c)
if cnt >= n:
c = c[:1]*0
else:
for i in range(cnt):
n = n - 1
c *= scl
der = np.empty((n,) + c.shape[1:], dtype=cdt)
for j in range(n, 0, -1):
der[j - 1] = j * c[j]
c = der
c = np.moveaxis(c, 0, iaxis)
return c
def polyint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
"""
Integrate a polynomial.
Returns the polynomial coefficients `c` integrated `m` times from
`lbnd` along `axis`. At each iteration the resulting series is
**multiplied** by `scl` and an integration constant, `k`, is added.
The scaling factor is for use in a linear change of variable. ("Buyer
beware": note that, depending on what one is doing, one may want `scl`
to be the reciprocal of what one might expect; for more information,
see the Notes section below.) The argument `c` is an array of
coefficients, from low to high degree along each axis, e.g., [1,2,3]
represents the polynomial ``1 + 2*x + 3*x**2`` while [[1,2],[1,2]]
represents ``1 + 1*x + 2*y + 2*x*y`` if axis=0 is ``x`` and axis=1 is
``y``.
Parameters
----------
c : array_like
1-D array of polynomial coefficients, ordered from low to high.
m : int, optional
Order of integration, must be positive. (Default: 1)
k : {[], list, scalar}, optional
Integration constant(s). The value of the first integral at zero
is the first value in the list, the value of the second integral
at zero is the second value, etc. If ``k == []`` (the default),
all constants are set to zero. If ``m == 1``, a single scalar can
be given instead of a list.
lbnd : scalar, optional
The lower bound of the integral. (Default: 0)
scl : scalar, optional
Following each integration the result is *multiplied* by `scl`
before the integration constant is added. (Default: 1)
axis : int, optional
Axis over which the integral is taken. (Default: 0).
.. versionadded:: 1.7.0
Returns
-------
S : ndarray
Coefficient array of the integral.
Raises
------
ValueError
If ``m < 1``, ``len(k) > m``, ``np.ndim(lbnd) != 0``, or
``np.ndim(scl) != 0``.
See Also
--------
polyder
Notes
-----
Note that the result of each integration is *multiplied* by `scl`. Why
is this important to note? Say one is making a linear change of
variable :math:`u = ax + b` in an integral relative to `x`. Then
:math:`dx = du/a`, so one will need to set `scl` equal to
:math:`1/a` - perhaps not what one would have first thought.
Examples
--------
>>> from numpy.polynomial import polynomial as P
>>> c = (1, 2, 3)
>>> P.polyint(c) # should return array([0, 1, 1, 1])
array([0., 1., 1., 1.])
>>> P.polyint(c, 3) # should return array([0, 0, 0, 1/6, 1/12, 1/20])
array([ 0. , 0. , 0. , 0.16666667, 0.08333333, # may vary
0.05 ])
>>> P.polyint(c, k=3) # should return array([3, 1, 1, 1])
array([3., 1., 1., 1.])
>>> P.polyint(c,lbnd=-2) # should return array([6, 1, 1, 1])
array([6., 1., 1., 1.])
"""
if m < 1 or len(k) > m:
raise ValueError("Invalid integration parameters.")
if np.ndim(lbnd) != 0:
raise ValueError("Lower bound lbnd must be scalar.")
if np.ndim(scl) != 0:
raise ValueError("Scaling factor scl must be scalar.")
for _ in range(m):
c = np.polyint(c, lbnd=lbnd, axis=axis) * scl + k.pop(0)
return c
c = np.array(c, ndmin=1, copy=True)
if c.dtype.char in '?bBhHiIlLqQpP':
c = c + 0.0
cdt = c.dtype
if not np.iterable(k):
k = [k]
cnt = pu._as_int(m, "the order of integration")
iaxis = pu._as_int(axis, "the axis")
if cnt < 0:
raise ValueError("The order of integration must be non-negative")
if len(k) > cnt:
raise ValueError("Too many integration constants")
if np.ndim(lbnd) != 0:
raise ValueError("lbnd must be a scalar.")
if np.ndim(scl) != 0:
raise ValueError("scl must be a scalar.")
iaxis = normalize_axis_index(iaxis, c.ndim)
if cnt == 0:
return c
k = list(k) + [0]*(cnt - len(k))
c = np.moveaxis(c, iaxis, 0)
for i in range(cnt):
n = len(c)
c *= scl
if n == 1 and np.all(c[0] == 0):
c[0] += k[i]
else:
tmp = np.empty((n + 1,) + c.shape[1:], dtype=cdt)
tmp[0] = c[0]*0
tmp[1] = c[0]
for j in range(1, n):
tmp[j + 1] = c[j]/(j + 1)
tmp[0] += k[i] - polyval(lbnd, tmp)
c = tmp
c = np.moveaxis(c, 0, iaxis)
return c
def polyval(x, c, tensor=True):
"""
Evaluate a polynomial at points x.
If `c` is of length ``n + 1``, this function returns the value
.. math:: p(x) = c_0 + c_1 * x + ... + c_n * x^n
The parameter `x` is converted to an array only if it is a tuple or a
list, otherwise it is treated as a scalar. In either case, either `x`
or its elements must support multiplication and addition both with
themselves and with the elements of `c`.
If `c` is a 1-D array, then ``p(x)`` will have the same shape as `x`. If
`c` is multidimensional, then the shape of the result depends on the
value of `tensor`. If `tensor` is true the shape will be c.shape[1:] +
x.shape. If `tensor` is false the shape will be c.shape[1:]. Note that
scalars have shape (,).
Trailing zeros in the coefficients will be used in the evaluation, so
they should be avoided if efficiency is a concern.
Parameters
----------
x : array_like, compatible object
If `x` is a list or tuple, it is converted to an ndarray, otherwise
it is left unchanged and treated as a scalar. In either case, `x`
or its elements must support addition and multiplication with
with themselves and with the elements of `c`.
c : array_like
Array of coefficients ordered so that the coefficients for terms of
degree n are contained in c[n]. If `c` is multidimensional the
remaining indices enumerate multiple polynomials. In the two
dimensional case the coefficients may be thought of as stored in
the columns of `c`.
tensor : boolean, optional
If True, the shape of the coefficient array is extended with ones
on the right, one for each dimension of `x`. Scalars have dimension 0
for this action. The result is that every column of coefficients in
`c` is evaluated for every element of `x`. If False, `x` is broadcast
over the columns of `c` for the evaluation. This keyword is useful
when `c` is multidimensional. The default value is True.
.. versionadded:: 1.7.0
Returns
-------
values : ndarray, compatible object
The shape of the returned array is described above.
See Also
--------
polyval2d, polygrid2d, polyval3d, polygrid3d
Notes
-----
The evaluation uses Horner's method.
Examples
--------
>>> from numpy.polynomial.polynomial import polyval
>>> polyval(1, [1,2,3])
6.0
>>> a = np.arange(4).reshape(2,2)
>>> a
array([[0, 1],
[2, 3]])
>>> polyval(a, [1, 2, 3])
array([[ 1., 6.],
[17., 34.]])
>>> coef = np.arange(4).reshape(2, 2) # multidimensional coefficients
>>> coef
array([[0, 1],
[2, 3]])
>>> polyval([1, 2], coef, tensor=True)
array([[2., 4.],
[4., 7.]])
>>> polyval([1, 2], coef, tensor=False)
array([2., 7.])
"""
c = np.array(c, ndmin=1, copy=None)
if c.dtype.char in '?bBhHiIlLqQpP':
c = c + 0.0
if isinstance(x, (tuple, list)):
x = np.asarray(x)
if isinstance(x, np.ndarray) and tensor:
c = c.reshape(c.shape + (1,)*x.ndim)
c0 = c[-1] + x*0
for i in range(2, len(c) + 1):
c0 = c[-i] + c0*x
return c0
def polyvalfromroots(x, r, tensor=True):
"""
Evaluate a polynomial specified by its roots at points x.
If `r` is of length ``N``, this function returns the value
.. math:: p(x) = \\prod_{n=1}^{N} (x - r_n)
The parameter `x` is converted to an array only if it is a tuple or a
list, otherwise it is treated as a scalar. In either case, either `x`
or its elements must support multiplication and addition both with
themselves and with the elements of `r`.
If `r` is a 1-D array, then ``p(x)`` will have the same shape as `x`. If `r`
is multidimensional, then the shape of the result depends on the value of
`tensor`. If `tensor` is ``True`` the shape will be r.shape[1:] + x.shape;
that is, each polynomial is evaluated at every value of `x`. If `tensor` is
``False``, the shape will be r.shape[1:]; that is, each polynomial is
evaluated only for the corresponding broadcast value of `x`. Note that
scalars have shape (,).
.. versionadded:: 1.12
Parameters
----------
x : array_like, compatible object
If `x` is a list or tuple, it is converted to an ndarray, otherwise
it is left unchanged and treated as a scalar. In either case, `x`
or its elements must support addition and multiplication with
with themselves and with the elements of `r`.
r : array_like
Array of roots. If `r` is multidimensional the first index is the
root index, while the remaining indices enumerate multiple
polynomials. For instance, in the two dimensional case the roots
of each polynomial may be thought of as stored in the columns of `r`.
tensor : boolean, optional
If True, the shape of the roots array is extended with ones on the
right, one for each dimension of `x`. Scalars have dimension 0 for this
action. The result is that every column of coefficients in `r` is
evaluated for every element of `x`. If False, `x` is broadcast over the
columns of `r` for the evaluation. This keyword is useful when `r` is
multidimensional. The default value is True.
Returns
-------
values : ndarray, compatible object
The shape of the returned array is described above.
See Also
--------
polyroots, polyfromroots, polyval
"""
if isinstance(x, (list, tuple)):
x = np.array(x)
if tensor:
values = np.prod(x.reshape(-1, 1) - r, axis=-1)
else:
values = np.prod(x - r, axis=0)
return values
"""
将输入的数组 r 转换为至少有一维的 numpy 数组，不进行复制
r = np.array(r, ndmin=1, copy=None)
如果 r 的数据类型在 '?bBhHiIlLqQpP' 中的一种，将其转换为 np.double 类型
if r.dtype.char in '?bBhHiIlLqQpP':
r = r.astype(np.double)
如果 x 是 tuple 或者 list 类型，则将其转换为 numpy 数组
if isinstance(x, (tuple, list)):
x = np.asarray(x)
如果 x 是 numpy 数组：
如果 tensor 为 True，则将 r 重塑为其形状加上 x 维度的 1 的形状
r = r.reshape(r.shape + (1,)*x.ndim)
否则，如果 x 的维度数 x.ndim 大于等于 r 的维度数 r.ndim，则抛出 ValueError
elif x.ndim >= r.ndim:
raise ValueError("x.ndim must be < r.ndim when tensor == False")
计算并返回 x 与 r 之间的元素级乘积，沿着第 0 轴
return np.prod(x - r, axis=0)
"""
def polyval2d(x, y, c):
return pu._valnd(polyval, c, x, y)
def polygrid2d(x, y, c):
"""
Evaluate a 2-D polynomial on the Cartesian product of x and y.
This function returns the values:
.. math:: p(a,b) = \\sum_{i,j} c_{i,j} * a^i * b^j
where the points ``(a, b)`` consist of all pairs formed by taking
`a` from `x` and `b` from `y`. The resulting points form a grid with
`x` in the first dimension and `y` in the second.
The parameters `x` and `y` are converted to arrays only if they are
tuples or a lists, otherwise they are treated as scalars. In either
case, either `x` and `y` or their elements must support multiplication
and addition both with themselves and with the elements of `c`.
If `c` has fewer than two dimensions, ones are implicitly appended to
its shape to make it 2-D. The shape of the result will be c.shape[2:] +
x.shape + y.shape.
Parameters
----------
x, y : array_like, compatible objects
The two dimensional series is evaluated at the points ``(x, y)``,
where `x` and `y` must have the same shape. If `x` or `y` is a list
or tuple, it is first converted to an ndarray, otherwise it is left
unchanged and, if it isn't an ndarray, it is treated as a scalar.
c : array_like
Array of coefficients ordered so that the coefficient of the term
of multi-degree i,j is contained in ``c[i,j]``. If `c` has
dimension greater than two the remaining indices enumerate multiple
sets of coefficients.
Returns
-------
values : ndarray, compatible object
The values of the two dimensional polynomial at points formed with
pairs of corresponding values from `x` and `y`.
See Also
--------
polyval, polygrid2d, polyval3d, polygrid3d
Notes
-----
.. versionadded:: 1.7.0
Examples
--------
>>> from numpy.polynomial import polynomial as P
>>> c = ((1, 2, 3), (4, 5, 6))
>>> P.polyval2d(1, 1, c)
21.0
"""
x, y : array_like, compatible objects
两个二维序列，将在`x`和`y`的笛卡尔积点处进行评估。如果`x`或`y`是列表或元组，首先将其转换为ndarray；否则保持不变，如果它不是ndarray，则将其视为标量。
c : array_like
按照顺序排列的系数数组，其中i,j次项的系数包含在``c[i,j]``中。如果`c`的维度大于二，则剩余的索引用于列举多组系数。
Returns
-------
values : ndarray, compatible object
在`x`和`y`的笛卡尔积点处的二维多项式的值。
See Also
--------
polyval, polyval2d, polyval3d, polygrid3d
Notes
-----
.. versionadded:: 1.7.0
Examples
--------
>>> from numpy.polynomial import polynomial as P
>>> c = ((1, 2, 3), (4, 5, 6))
>>> P.polygrid2d([0, 1], [0, 1], c)
array([[ 1., 6.],
[ 5., 21.]])
"""
return pu._gridnd(polyval, c, x, y)
# 计算三维多项式在点 (x, y, z) 处的值。
# 导入模块，使用 _valnd 函数计算多维多项式的值
def polyval3d(x, y, z, c):
return pu._valnd(polyval, c, x, y, z)
# 计算三维多项式在 x, y, z 的笛卡尔积上的值。
# 导入模块，使用 _valnd 函数计算三维网格上的多项式值
def polygrid3d(x, y, z, c):
x, y, z : array_like, compatible objects
三维系列在笛卡尔积 `x`, `y`, 和 `z` 点处进行评估。如果 `x`, `y`, 或 `z` 是列表或元组，它们首先被转换为 ndarray 数组；否则保持不变，并且如果它不是 ndarray 数组，则被视为标量。
c : array_like
系数数组，按照这样的顺序排列：i, j 次项的系数包含在 ``c[i,j]`` 中。如果 `c` 的维度大于两，剩余的索引用来枚举多组系数。
Returns
-------
values : ndarray, compatible object
二维多项式在 `x` 和 `y` 的笛卡尔积点处的值。
See Also
--------
polyval, polyval2d, polygrid2d, polyval3d
Notes
-----
.. versionadded:: 1.7.0
添加于版本 1.7.0
Examples
--------
>>> from numpy.polynomial import polynomial as P
>>> c = ((1, 2, 3), (4, 5, 6), (7, 8, 9))
>>> P.polygrid3d([0, 1], [0, 1], [0, 1], c)
array([[ 1., 13.],
[ 6., 51.]])
"""
return pu._gridnd(polyval, c, x, y, z)
def polyvander(x, deg):
"""Vandermonde matrix of given degree.
Returns the Vandermonde matrix of degree `deg` and sample points
`x`. The Vandermonde matrix is defined by
.. math:: V[..., i] = x^i,
where ``0 <= i <= deg``. The leading indices of `V` index the elements of
`x` and the last index is the power of `x`.
If `c` is a 1-D array of coefficients of length ``n + 1`` and `V` is the
matrix ``V = polyvander(x, n)``, then ``np.dot(V, c)`` and
``polyval(x, c)`` are the same up to roundoff. This equivalence is
useful both for least squares fitting and for the evaluation of a large
number of polynomials of the same degree and sample points.
Parameters
----------
x : array_like
Array of points. The dtype is converted to float64 or complex128
depending on whether any of the elements are complex. If `x` is
scalar it is converted to a 1-D array.
deg : int
Degree of the resulting matrix.
Returns
-------
vander : ndarray.
The Vandermonde matrix. The shape of the returned matrix is
``x.shape + (deg + 1,)``, where the last index is the power of `x`.
The dtype will be the same as the converted `x`.
See Also
--------
polyvander2d, polyvander3d
Examples
--------
The Vandermonde matrix of degree ``deg = 5`` and sample points
``x = [-1, 2, 3]`` contains the element-wise powers of `x`
from 0 to 5 as its columns.
>>> from numpy.polynomial import polynomial as P
>>> x, deg = [-1, 2, 3], 5
>>> P.polyvander(x=x, deg=deg)
array([[ 1., -1., 1., -1., 1., -1.],
[ 1., 2., 4., 8., 16., 32.],
[ 1., 3., 9., 27., 81., 243.]])
"""
ideg = pu._as_int(deg, "deg")
if ideg < 0:
raise ValueError("deg must be non-negative")
x = np.array(x, copy=None, ndmin=1) + 0.0
dims = (ideg + 1,) + x.shape
dtyp = x.dtype
v = np.empty(dims, dtype=dtyp)
v[0] = x*0 + 1
if ideg > 0:
v[1] = x
for i in range(2, ideg + 1):
v[i] = v[i-1]*x
return np.moveaxis(v, 0, -1)
"""
return pu._vander_nd_flat((polyvander, polyvander), (x, y), deg)
# 多项式Vandermonde矩阵生成函数，用于三维情况
def polyvander3d(x, y, z, deg):
"""Pseudo-Vandermonde matrix of given degrees.
Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
points ``(x, y, z)``. If `l`, `m`, `n` are the given degrees in `x`, `y`, `z`,
then The pseudo-Vandermonde matrix is defined by
.. math:: V[..., (m+1)(n+1)i + (n+1)j + k] = x^i * y^j * z^k,
where ``0 <= i <= l``, ``0 <= j <= m``, and ``0 <= k <= n``. The leading
indices of `V` index the points ``(x, y, z)`` and the last index encodes
the powers of `x`, `y`, and `z`.
If ``V = polyvander3d(x, y, z, [xdeg, ydeg, zdeg])``, then the columns
of `V` correspond to the elements of a 3-D coefficient array `c` of
shape (xdeg + 1, ydeg + 1, zdeg + 1) in the order
.. math:: c_{000}, c_{001}, c_{002},... , c_{010}, c_{011}, c_{012},...
and ``np.dot(V, c.flat)`` and ``polyval3d(x, y, z, c)`` will be the
same up to roundoff. This equivalence is useful both for least squares
fitting and for the evaluation of a large number of 3-D polynomials
of the same degrees and sample points.
Parameters
----------
x, y, z : array_like
Arrays of point coordinates, all of the same shape. The dtypes will
be converted to either float64 or complex128 depending on whether
any of the elements are complex. Scalars are converted to 1-D
arrays.
deg : list of ints
List of maximum degrees of the form [x_deg, y_deg, z_deg].
Returns
-------
vander3d : ndarray
The shape of the returned matrix is ``x.shape + (order,)``, where
:math:`order = (deg[0]+1)*(deg([1]+1)*(deg[2]+1)`. The dtype will
be the same as the converted `x`, `y`, and `z`.
See Also
--------
polyvander, polyvander3d, polyval2d, polyval3d
Notes
-----
.. versionadded:: 1.7.0
Examples
--------
>>> from numpy.polynomial import polynomial as P
>>> x = np.asarray([-1, 2, 1])
>>> y = np.asarray([1, -2, -3])
>>> z = np.asarray([2, 2, 5])
>>> l, m, n = [2, 2, 1]
>>> deg = [l, m, n]
>>> V = P.polyvander3d(x=x, y=y, z=z, deg=deg)
>>> V
array([[ 1., 2., 1., 2., 1., 2., -1., -2., -1.,
-2., -1., -2., 1., 2., 1., 2., 1., 2.],
[ 1., 2., -2., -4., 4., 8., 2., 4., -4.,
-8., 8., 16., 4., 8., -8., -16., 16., 32.],
[ 1., 5., -3., -15., 9., 45., 1., 5., -3.,
-15., 9., 45., 1., 5., -3., -15., 9., 45.]])
We can verify the columns for any ``0 <= i <= l``, ``0 <= j <= m``,
and ``0 <= k <= n``
>>> i, j, k = 2, 1, 0
>>> V[:, (m+1)*(n+1)*i + (n+1)*j + k] == x**i * y**j * z**k
array([ True, True, True])
"""
# 调用私有函数 `_vander_nd_flat` 生成高维Vandermonde矩阵
return pu._vander_nd_flat((polyvander, polyvander, polyvander), (x, y, z), deg)
# 返回多项式拟合的系数，该多项式的阶数为 `deg`，通过最小二乘法拟合给定点 `x` 上的数据值 `y`。
# 如果 `y` 是一维的，则返回的系数也是一维的。如果 `y` 是二维的，则对每一列 `y` 进行多次拟合，
# 并将结果系数存储在返回的二维数组的相应列中。
Return the coefficients of a polynomial of degree `deg` that is the
least squares fit to the data values `y` given at points `x`. If `y` is
1-D the returned coefficients will also be 1-D. If `y` is 2-D multiple
fits are done, one for each column of `y`, and the resulting
coefficients are stored in the corresponding columns of a 2-D return.
The fitted polynomial(s) are in the form
# 多项式的形式为
.. math:: p(x) = c_0 + c_1 * x + ... + c_n * x^n,
where `n` is `deg`.
Parameters
----------
x : array_like, shape (`M`,)
x-coordinates of the `M` sample (data) points ``(x[i], y[i])``.
y : array_like, shape (`M`,) or (`M`, `K`)
y-coordinates of the sample points. Several sets of sample points
sharing the same x-coordinates can be (independently) fit with one
call to `polyfit` by passing in for `y` a 2-D array that contains
one data set per column.
deg : int or 1-D array_like
Degree(s) of the fitting polynomials. If `deg` is a single integer
all terms up to and including the `deg`'th term are included in the
fit. For NumPy versions >= 1.11.0 a list of integers specifying the
degrees of the terms to include may be used instead.
rcond : float, optional
Relative condition number of the fit. Singular values smaller
than `rcond`, relative to the largest singular value, will be
ignored. The default value is ``len(x)*eps``, where `eps` is the
relative precision of the platform's float type, about 2e-16 in
most cases.
full : bool, optional
Switch determining the nature of the return value. When ``False``
(the default) just the coefficients are returned; when ``True``,
diagnostic information from the singular value decomposition (used
to solve the fit's matrix equation) is also returned.
w : array_like, shape (`M`,), optional
Weights. If not None, the weight ``w[i]`` applies to the unsquared
residual ``y[i] - y_hat[i]`` at ``x[i]``. Ideally the weights are
chosen so that the errors of the products ``w[i]*y[i]`` all have the
same variance. When using inverse-variance weighting, use
``w[i] = 1/sigma(y[i])``. The default value is None.
.. versionadded:: 1.5.0
Returns
-------
coef : ndarray, shape (`deg` + 1,) or (`deg` + 1, `K`)
Polynomial coefficients ordered from low to high. If `y` was 2-D,
the coefficients in column `k` of `coef` represent the polynomial
fit to the data in `y`'s `k`-th column.
[residuals, rank, singular_values, rcond] : list
# 返回值列表，仅当 `full == True` 时返回以下值：
- residuals -- 最小二乘拟合的残差平方和
- rank -- 缩放后的范德蒙德矩阵的数值秩
- singular_values -- 缩放后的范德蒙德矩阵的奇异值
- rcond -- `rcond` 的值
更多细节，请参阅 `numpy.linalg.lstsq`。
Raises
------
RankWarning
# 如果最小二乘拟合中的矩阵秩不足，将引发此警告。
仅当 `full == False` 时才会引发警告。可以通过以下方式关闭警告：
>>> import warnings
>>> warnings.simplefilter('ignore', np.exceptions.RankWarning)
See Also
--------
numpy.polynomial.chebyshev.chebfit
numpy.polynomial.legendre.legfit
numpy.polynomial.laguerre.lagfit
numpy.polynomial.hermite.hermfit
numpy.polynomial.hermite_e.hermefit
polyval : 计算多项式值。
polyvander : 用于幂的范德蒙德矩阵。
numpy.linalg.lstsq : 计算矩阵的最小二乘拟合。
scipy.interpolate.UnivariateSpline : 计算样条拟合。
Notes
-----
解决方案是多项式 `p` 的系数，该多项式最小化加权平方误差的和
.. math:: E = \\sum_j w_j^2 * |y_j - p(x_j)|^2,
其中 :math:`w_j` 是权重。该问题通过建立（通常为）过度确定的矩阵方程解决：
.. math:: V(x) * c = w * y,
其中 `V` 是 `x` 的加权伪范德蒙德矩阵，`c` 是要解的系数，`w` 是权重，`y` 是观察值。然后使用 `V` 的奇异值分解来解决此方程。
如果 `V` 的某些奇异值很小以至于被忽略（且 `full` == ``False``），将引发 `~exceptions.RankWarning`。这意味着系数值可能确定不好。通常通过拟合较低阶多项式可以消除警告（但这可能不是您想要的；如果您有选择阶数的独立理由不起作用，您可能需要：a）重新考虑这些理由，和/或 b）重新考虑数据质量）。`rcond` 参数也可以设置为小于其默认值的值，但所得拟合可能是虚假的，并且可能会受到舍入误差的大贡献。
使用双精度进行多项式拟合通常在多项式阶数约为 20 时“失败”。使用Chebyshev或Legendre系列进行拟合通常条件更好，但很大程度上取决于样本点的分布和数据的平滑性。如果拟合质量不够好，样条可能是一个好的选择。
Examples
--------
>>> from numpy.polynomial import polynomial as P
>>> x = np.linspace(-1,1,51) # 生成包含51个元素的数组x，范围从-1到1，等间隔分布
>>> rng = np.random.default_rng()
>>> err = rng.normal(size=len(x)) # 创建长度与x相同的正态分布误差数组err
>>> y = x**3 - x + err # 根据公式y = x^3 - x + Gaussian noise计算y数组
>>> c, stats = P.polyfit(x,y,3,full=True) # 对x和y进行3次多项式拟合，返回拟合系数c和统计信息stats
>>> c # c[0], c[1] 大约为 -1，c[2] 应该接近于0，c[3] 大约为 1
array([ 0.23111996, -1.02785049, -0.2241444 , 1.08405657]) # 可能有所变化
>>> stats # 注意到较大的SSR，解释了拟合结果较差的原因
[array([48.312088]), # 可能有所变化
4,
array([1.38446749, 1.32119158, 0.50443316, 0.28853036]),
1.1324274851176597e-14]
Same thing without the added noise # 在没有额外噪声的情况下进行相同的操作
>>> y = x**3 - x # 根据公式y = x^3 - x 计算y数组
>>> c, stats = P.polyfit(x,y,3,full=True) # 对x和y进行3次多项式拟合，返回拟合系数c和统计信息stats
>>> c # c[0], c[1] 约为 -1，c[2] 应该非常接近于0，c[3] 约为 1
array([-6.73496154e-17, -1.00000000e+00, 0.00000000e+00, 1.00000000e+00])
>>> stats # 注意到极小的SSR
[array([8.79579319e-31]), # 可能有所变化
4,
array([1.38446749, 1.32119158, 0.50443316, 0.28853036]),
1.1324274851176597e-14]
"""
return pu._fit(polyvander, x, y, deg, rcond, full, w)
注释：
- `x = np.linspace(-1,1,51)`：生成包含51个元素的数组x，范围从-1到1，等间隔分布。
- `rng = np.random.default_rng()`：创建一个默认的随机数生成器对象。
- `err = rng.normal(size=len(x))`：创建一个长度与x相同的数组，其中的元素符合标准正态分布。
- `y = x**3 - x + err`：根据公式y = x^3 - x + Gaussian noise 计算y数组，将x的立方、x本身和误差err相加。
- `c, stats = P.polyfit(x,y,3,full=True)`：使用3次多项式拟合函数`P.polyfit`拟合x和y，返回拟合系数c和拟合统计信息stats。
- `c`：拟合系数数组c，其中c[0]、c[1]等分别对应多项式的各阶系数。
- `stats`：拟合的统计信息数组，包含SSR（残差平方和）、拟合自由度、标准误差估计和拟合使用的条件数等。
- `return pu._fit(polyvander, x, y, deg, rcond, full, w)`：返回调用`pu._fit`函数的结果，该函数用于多项式拟合过程中计算内部使用的Vandermonde矩阵和其他相关信息。
[c] = pu.as_series([c])
if len(c) < 2:
raise ValueError('Series must have maximum degree of at least 1.')
if len(c) == 2:
return np.array([-c[0]/c[1]])
n = len(c) - 1
mat = np.zeros((n, n), dtype=c.dtype)
bot = mat.reshape(-1)[n::n+1]
bot[...] = 1
mat[:, -1] -= c[:-1] / c[-1]
return mat
return r
class Polynomial(ABCPolyBase):
"""A power series class.
The Polynomial class provides the standard Python numerical methods
'+', '-', '*', '//', '%', 'divmod', '**', and '()' as well as the
attributes and methods listed below.
Parameters
----------
coef : array_like
Polynomial coefficients in order of increasing degree, i.e.,
``(1, 2, 3)`` give ``1 + 2*x + 3*x**2``.
domain : (2,) array_like, optional
Domain to use. The interval ``[domain[0], domain[1]]`` is mapped
to the interval ``[window[0], window[1]]`` by shifting and scaling.
The default value is [-1, 1].
window : (2,) array_like, optional
Window, see `domain` for its use. The default value is [-1, 1].
.. versionadded:: 1.6.0
symbol : str, optional
Symbol used to represent the independent variable in string
representations of the polynomial expression, e.g. for printing.
The symbol must be a valid Python identifier. Default value is 'x'.
.. versionadded:: 1.24
"""
_add = staticmethod(polyadd)
_sub = staticmethod(polysub)
_mul = staticmethod(polymul)
_div = staticmethod(polydiv)
_pow = staticmethod(polypow)
_val = staticmethod(polyval)
_int = staticmethod(polyint)
_der = staticmethod(polyder)
_fit = staticmethod(polyfit)
_line = staticmethod(polyline)
_roots = staticmethod(polyroots)
_fromroots = staticmethod(polyfromroots)
domain = np.array(polydomain)
window = np.array(polydomain)
basis_name = None
@classmethod
def _str_term_unicode(cls, i, arg_str):
"""Generate a unicode representation of a term in the polynomial.
Parameters
----------
i : int
The exponent of the term.
arg_str : str
The string representation of the term without exponent.
Returns
-------
str
Unicode representation of the term.
"""
if i == '1':
return f"·{arg_str}"
else:
return f"·{arg_str}{i.translate(cls._superscript_mapping)}"
@staticmethod
def _str_term_ascii(i, arg_str):
"""Generate an ASCII representation of a term in the polynomial.
Parameters
----------
i : int
The exponent of the term.
arg_str : str
The string representation of the term without exponent.
Returns
-------
str
ASCII representation of the term.
"""
if i == '1':
return f" {arg_str}"
else:
return f" {arg_str}**{i}"
@staticmethod
def _repr_latex_term(i, arg_str, needs_parens):
"""Generate a LaTeX representation of a term in the polynomial.
Parameters
----------
i : int
The exponent of the term.
arg_str : str
The LaTeX string representation of the term without exponent.
needs_parens : bool
Whether the term needs to be wrapped in parentheses.
Returns
-------
str
LaTeX representation of the term.
"""
if needs_parens:
arg_str = rf"\left({arg_str}\right)"
if i == 0:
return '1'
elif i == 1:
return arg_str
else:
return f"{arg_str}^{{{i}}}"
.\numpy\numpy\polynomial\polynomial.pyi
from typing import Any
from numpy import int_
from numpy.typing import NDArray
from numpy.polynomial._polybase import ABCPolyBase
from numpy.polynomial.polyutils import trimcoef
__all__: list[str]
polytrim = trimcoef
polydomain: NDArray[int_]
polyzero: NDArray[int_]
polyone: NDArray[int_]
polyx: NDArray[int_]
def polyline(off, scl): ...
def polyfromroots(roots): ...
def polyadd(c1, c2): ...
def polysub(c1, c2): ...
def polymulx(c): ...
def polymul(c1, c2): ...
def polydiv(c1, c2): ...
def polypow(c, pow, maxpower=...): ...
def polyder(c, m=..., scl=..., axis=...): ...
def polyint(c, m=..., k=..., lbnd=..., scl=..., axis=...): ...
def polyval(x, c, tensor=...): ...
def polyvalfromroots(x, r, tensor=...): ...
def polyval2d(x, y, c): ...
def polygrid2d(x, y, c): ...
def polyval3d(x, y, z, c): ...
def polygrid3d(x, y, z, c): ...
def polyvander(x, deg): ...
def polyvander2d(x, y, deg): ...
def polyvander3d(x, y, z, deg): ...
def polyfit(x, y, deg, rcond=..., full=..., w=...): ...
def polyroots(c): ...
class Polynomial(ABCPolyBase):
domain: Any
window: Any
basis_name: Any
.\numpy\numpy\polynomial\polyutils.py
import operator
import functools
import warnings
import numpy as np
from numpy._core.multiarray import dragon4_positional, dragon4_scientific
from numpy.exceptions import RankWarning
__all__ = [
'as_series', 'trimseq', 'trimcoef', 'getdomain', 'mapdomain', 'mapparms',
'format_float']
def trimseq(seq):
"""Remove small Poly series coefficients.
Parameters
----------
seq : sequence
Sequence of Poly series coefficients.
Returns
-------
series : sequence
Subsequence with trailing zeros removed. If the resulting sequence
would be empty, return the first element. The returned sequence may
or may not be a view.
Notes
-----
Do not lose the type info if the sequence contains unknown objects.
"""
if len(seq) == 0 or seq[-1] != 0:
return seq
else:
for i in range(len(seq) - 1, -1, -1):
if seq[i] != 0:
break
return seq[:i+1]
def as_series(alist, trim=True):
"""
Return argument as a list of 1-d arrays.
The returned list contains array(s) of dtype double, complex double, or
object. A 1-d argument of shape ``(N,)`` is parsed into ``N`` arrays of
size one; a 2-d argument of shape ``(M,N)`` is parsed into ``M`` arrays
of size ``N`` (i.e., is "parsed by row"); and a higher dimensional array
raises a Value Error if it is not first reshaped into either a 1-d or 2-d
array.
Parameters
----------
alist : array_like
A 1- or 2-d array_like
trim : boolean, optional
When True, trailing zeros are removed from the inputs.
When False, the inputs are passed through intact.
Returns
-------
[a1, a2,...] : list of 1-D arrays
A copy of the input data as a list of 1-d arrays.
Raises
------
ValueError
Raised when `as_series` cannot convert its input to 1-d arrays, or at
least one of the resulting arrays is empty.
Examples
--------
>>> from numpy.polynomial import polyutils as pu
>>> a = np.arange(4)
>>> pu.as_series(a)
[array([0.]), array([1.]), array([2.]), array([3.])]
>>> b = np.arange(6).reshape((2,3))
>>> pu.as_series(b)
[array([0., 1., 2.]), array([3., 4., 5.])]
"""
if trim:
return [trimseq(a) for a in np.atleast_1d(alist)]
else:
return [np.asarray(a).flatten() for a in np.atleast_1d(alist)]
"""
将输入的列表转换为 NumPy 数组的列表，每个数组至少有一维，并进行必要的验证和转换操作。
Parameters:
- alist: 输入的列表，每个元素将被转换为一个 NumPy 数组
- trim: 可选参数，默认为 True。如果为 True，则对每个数组进行修剪操作。
Returns:
- ret: 转换后的 NumPy 数组的列表
Raises:
- ValueError: 如果任何一个数组为空，或者不是一维数组。
- ValueError: 如果无法找到数组的公共数据类型。
"""
arrays = [np.array(a, ndmin=1, copy=None) for a in alist]
for a in arrays:
if a.size == 0:
raise ValueError("Coefficient array is empty")
if any(a.ndim != 1 for a in arrays):
raise ValueError("Coefficient array is not 1-d")
if trim:
arrays = [trimseq(a) for a in arrays]
if any(a.dtype == np.dtype(object) for a in arrays):
ret = []
for a in arrays:
if a.dtype != np.dtype(object):
tmp = np.empty(len(a), dtype=np.dtype(object))
tmp[:] = a[:]
ret.append(tmp)
else:
ret.append(a.copy())
else:
try:
dtype = np.common_type(*arrays)
except Exception as e:
raise ValueError("Coefficient arrays have no common type") from e
ret = [np.array(a, copy=True, dtype=dtype) for a in arrays]
return ret
def trimcoef(c, tol=0):
"""
Remove "small" "trailing" coefficients from a polynomial.
"Small" means "small in absolute value" and is controlled by the
parameter `tol`; "trailing" means highest order coefficient(s), e.g., in
``[0, 1, 1, 0, 0]`` (which represents ``0 + x + x**2 + 0*x**3 + 0*x**4``)
both the 3-rd and 4-th order coefficients would be "trimmed."
Parameters
----------
c : array_like
1-d array of coefficients, ordered from lowest order to highest.
tol : number, optional
Trailing (i.e., highest order) elements with absolute value less
than or equal to `tol` (default value is zero) are removed.
Returns
-------
trimmed : ndarray
1-d array with trailing zeros removed. If the resulting series
would be empty, a series containing a single zero is returned.
Raises
------
ValueError
If `tol` < 0
Examples
--------
>>> from numpy.polynomial import polyutils as pu
>>> pu.trimcoef((0,0,3,0,5,0,0))
array([0., 0., 3., 0., 5.])
>>> pu.trimcoef((0,0,1e-3,0,1e-5,0,0),1e-3) # item == tol is trimmed
array([0.])
>>> i = complex(0,1) # works for complex
>>> pu.trimcoef((3e-4,1e-3*(1-i),5e-4,2e-5*(1+i)), 1e-3)
array([0.0003+0.j , 0.001 -0.001j])
"""
if tol < 0:
raise ValueError("tol must be non-negative")
[c] = as_series([c])
[ind] = np.nonzero(np.abs(c) > tol)
if len(ind) == 0:
return c[:1]*0
else:
return c[:ind[-1] + 1].copy()
def getdomain(x):
"""
Return a domain suitable for given abscissae.
Find a domain suitable for a polynomial or Chebyshev series
defined at the values supplied.
Parameters
----------
x : array_like
1-d array of abscissae whose domain will be determined.
Returns
-------
domain : ndarray
1-d array containing two values. If the inputs are complex, then
the two returned points are the lower left and upper right corners
of the smallest rectangle (aligned with the axes) in the complex
plane containing the points `x`. If the inputs are real, then the
two points are the ends of the smallest interval containing the
points `x`.
See Also
--------
mapparms, mapdomain
Examples
--------
>>> from numpy.polynomial import polyutils as pu
>>> points = np.arange(4)**2 - 5; points
array([-5, -4, -1, 4])
>>> pu.getdomain(points)
array([-5., 4.])
>>> c = np.exp(complex(0,1)*np.pi*np.arange(12)/6) # unit circle
>>> pu.getdomain(c)
array([-1.-1.j, 1.+1.j])
"""
[x] = as_series([x], trim=False)
if x.dtype.char in np.typecodes['Complex']:
rmin, rmax = x.real.min(), x.real.max()
imin, imax = x.imag.min(), x.imag.max()
return np.array((complex(rmin, imin), complex(rmax, imax)))
else:
return np.array((x.min(), x.max()))
"""
Linear map parameters between domains.
Return the parameters of the linear map ``offset + scale*x`` that maps
`old` to `new` such that ``old[i] -> new[i]``, ``i = 0, 1``.
Parameters
----------
old, new : array_like
Domains. Each domain must (successfully) convert to a 1-d array
containing precisely two values.
Returns
-------
offset, scale : scalars
The map ``L(x) = offset + scale*x`` maps the first domain to the
second.
See Also
--------
getdomain, mapdomain
Notes
-----
Also works for complex numbers, and thus can be used to calculate the
parameters required to map any line in the complex plane to any other
line therein.
Examples
--------
>>> from numpy.polynomial import polyutils as pu
>>> pu.mapparms((-1,1),(-1,1))
(0.0, 1.0)
>>> pu.mapparms((1,-1),(-1,1))
(-0.0, -1.0)
>>> i = complex(0,1)
>>> pu.mapparms((-i,-1),(1,i))
((1+1j), (1-0j))
"""
oldlen = old[1] - old[0]
newlen = new[1] - new[0]
off = (old[1]*new[0] - old[0]*new[1]) / oldlen
scl = newlen / oldlen
return off, scl
def mapdomain(x, old, new):
"""
Apply linear map to input points.
The linear map ``offset + scale*x`` that maps the domain `old` to
the domain `new` is applied to the points `x`.
Parameters
----------
x : array_like
Points to be mapped. If `x` is a subtype of ndarray the subtype
will be preserved.
old, new : array_like
The two domains that determine the map. Each must (successfully)
convert to 1-d arrays containing precisely two values.
Returns
-------
x_out : ndarray
Array of points of the same shape as `x`, after application of the
linear map between the two domains.
See Also
--------
getdomain, mapparms
Notes
-----
Effectively, this implements:
.. math::
x\\_out = new[0] + m(x - old[0])
where
.. math::
m = \\frac{new[1]-new[0]}{old[1]-old[0]}
Examples
--------
>>> from numpy.polynomial import polyutils as pu
>>> old_domain = (-1,1)
>>> new_domain = (0,2*np.pi)
>>> x = np.linspace(-1,1,6); x
array([-1. , -0.6, -0.2, 0.2, 0.6, 1. ])
>>> x_out = pu.mapdomain(x, old_domain, new_domain); x_out
array([ 0. , 1.25663706, 2.51327412, 3.76991118, 5.02654825, # may vary
6.28318531])
>>> x - pu.mapdomain(x_out, new_domain, old_domain)
array([0., 0., 0., 0., 0., 0.])
Also works for complex numbers (and thus can be used to map any line in
the complex plane to any other line therein).
>>> i = complex(0,1)
>>> old = (-1 - i, 1 + i)
>>> new = (-1 + i, 1 - i)
>>> z = np.linspace(old[0], old[1], 6); z
array([-1. -1.j , -0.6-0.6j, -0.2-0.2j, 0.2+0.2j, 0.6+0.6j, 1. +1.j ])
>>> new_z = pu.mapdomain(z, old, new); new_z
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
"""
x = np.asanyarray(x)
off, scl = mapparms(old, new)
return off + scl*x
def _nth_slice(i, ndim):
"""
Generate a slice object that selects the i-th index in an array.
Parameters
----------
i : int
Index to be selected.
ndim : int
Number of dimensions in the array.
Returns
-------
tuple
Tuple of slice objects where the i-th position is sliced with `slice(None)`.
"""
sl = [np.newaxis] * ndim
sl[i] = slice(None)
return tuple(sl)
def _vander_nd(vander_fs, points, degrees):
r"""
A generalization of the Vandermonde matrix for N dimensions
The result is built by combining the results of 1d Vandermonde matrices,
.. math::
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]}
where
.. math::
N &= \texttt{len(points)} = \texttt{len(degrees)} = \texttt{len(vander\_fs)} \\
M &= \texttt{points[k].ndim} \\
V_k &= \texttt{vander\_fs[k]} \\
x_k &= \texttt{points[k]} \\
0 \le j_k &\le \texttt{degrees[k]}
Expanding the one-dimensional :math:`V_k` functions gives:
.. math::
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])}
where :math:`B_{k,m}` is the m'th basis of the polynomial construction used along
dimension :math:`k`. For a regular polynomial, :math:`B_{k, m}(x) = P_m(x) = x^m`.
"""
pass
Parameters
----------
vander_fs : Sequence[function(array_like, int) -> ndarray]
每个轴上要使用的一维范德蒙德函数，例如 `polyvander`
points : Sequence[array_like]
点的坐标数组，所有数组的形状必须相同。其元素的数据类型将转换为 float64 或 complex128，具体取决于是否包含复数。标量将转换为一维数组。
这个参数的长度必须与 `vander_fs` 相同。
degrees : Sequence[int]
每个轴要使用的最大阶数（包括在内）。这个参数的长度必须与 `vander_fs` 相同。
Returns
-------
vander_nd : ndarray
形状为 ``points[0].shape + tuple(d + 1 for d in degrees)`` 的数组。
"""
# 确定维度数量
n_dims = len(vander_fs)
# 检查样本点的维度数量是否与给定的函数数量相同
if n_dims != len(points):
raise ValueError(
f"Expected {n_dims} dimensions of sample points, got {len(points)}")
# 检查阶数的维度数量是否与给定的函数数量相同
if n_dims != len(degrees):
raise ValueError(
f"Expected {n_dims} dimensions of degrees, got {len(degrees)}")
# 如果没有提供样本点，则无法猜测 dtype 或形状
if n_dims == 0:
raise ValueError("Unable to guess a dtype or shape when no points are given")
# 将所有点转换为相同的形状和类型
points = tuple(np.asarray(tuple(points)) + 0.0)
# 为每个维度生成范德蒙德矩阵，将每个维度的最后一个轴放置在输出的独立尾随轴中
vander_arrays = (
vander_fs[i](points[i], degrees[i])[(...,) + _nth_slice(i, n_dims)]
for i in range(n_dims)
)
# 我们已经检查过这不是空的，因此不需要 `initial`
return functools.reduce(operator.mul, vander_arrays)
def _vander_nd_flat(vander_fs, points, degrees):
"""
Like `_vander_nd`, but flattens the last ``len(degrees)`` axes into a single axis
Used to implement the public ``<type>vander<n>d`` functions.
"""
# 调用 `_vander_nd` 函数，获取多维 Vandermonde 矩阵
v = _vander_nd(vander_fs, points, degrees)
# 将多维矩阵展平为最后 ``len(degrees)`` 个轴合并到一个轴
return v.reshape(v.shape[:-len(degrees)] + (-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
"""
# 如果 roots 为空，则返回一个包含一个元素的数组
if len(roots) == 0:
return np.ones(1)
else:
# 将 roots 转换为系列，并进行排序
[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
"""
# 将所有参数转换为 NumPy 数组
args = [np.asanyarray(a) for a in args]
shape0 = args[0].shape
# 检查所有参数的形状是否与第一个参数相同
if not all((a.shape == shape0 for a in args[1:])):
if len(args) == 3:
raise ValueError('x, y, z are incompatible')
elif len(args) == 2:
raise ValueError('x, y are incompatible')
else:
raise ValueError('ordinates are incompatible')
it = iter(args)
x0 = next(it)
# 对第一个参数使用 val_f 函数
c = val_f(x0, c)
# 对剩余的参数依次使用 val_f 函数
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
"""
# 对所有参数依次使用 val_f 函数
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
# 此处是对 ``<type>div`` 函数的简要描述和提示
# 将输入的两个参数 c1 和 c2 转换为序列（数组或类似结构），分别存储到 c1 和 c2 中
[c1, c2] = as_series([c1, c2])
# 检查 c2 数组的最后一个元素是否为 0，如果是则抛出 ZeroDivisionError 异常
if c2[-1] == 0:
raise ZeroDivisionError()
# 计算 c1 和 c2 的长度，分别存储到 lc1 和 lc2 中
lc1 = len(c1)
lc2 = len(c2)
# 如果 c1 的长度小于 c2 的长度，返回一个以 c1[0] 为元素、长度为 0 的数组，并返回 c1
if lc1 < lc2:
return c1[:1]*0, c1
# 如果 c2 的长度为 1，返回 c1 除以 c2 的最后一个元素，以及一个以 c1[0] 为元素、长度为 0 的数组
elif lc2 == 1:
return c1/c2[-1], c1[:1]*0
else:
# 创建一个类型与 c1 相同的空数组 quo，长度为 lc1 - lc2 + 1
quo = np.empty(lc1 - lc2 + 1, dtype=c1.dtype)
# 将 rem 初始化为 c1
rem = c1
# 从 lc1 - lc2 循环到 0（包括 0）
for i in range(lc1 - lc2, - 1, -1):
# 创建一个多项式 p，其中 i 个 0 后接 1，然后与 c2 相乘
p = mul_f([0]*i + [1], c2)
# 计算余数 rem 的最后一个元素与 p 的最后一个元素的除法商 q
q = rem[-1]/p[-1]
# 更新 rem 为 rem 去除最后一个元素减去 q 乘以 p 去除最后一个元素
rem = rem[:-1] - q*p[:-1]
# 将 q 存储到 quo 的第 i 个位置
quo[i] = q
# 返回商 quo 和修剪后的余数 trimseq(rem)
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]) # 调用as_series函数，将c1和c2转换为Series类型，并进行必要的修剪
if len(c1) > len(c2):
c1[:c2.size] += c2 # 如果c1长度大于c2，则将c2加到c1的前c2.size部分
ret = c1
else:
c2[:c1.size] += c1 # 否则，将c1加到c2的前c1.size部分
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]) # 调用as_series函数，将c1和c2转换为Series类型，并进行必要的修剪
if len(c1) > len(c2):
c1[:c2.size] -= c2 # 如果c1长度大于c2，则从c1的前c2.size部分减去c2
ret = c1
else:
c2 = -c2 # 否则，将c2取反
c2[:c1.size] += c1 # 然后将c1加到c2的前c1.size部分
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``
x, y, deg
See the ``<type>fit`` functions for more detail
"""
x = np.asarray(x) + 0.0 # 将x转换为NumPy数组，并确保是浮点型
y = np.asarray(y) + 0.0 # 将y转换为NumPy数组，并确保是浮点型
deg = np.asarray(deg) # 将deg转换为NumPy数组
# 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") # 检查deg的类型和形状是否符合要求
if deg.min() < 0:
raise ValueError("expected deg >= 0") # 检查deg是否都大于等于0
if x.ndim != 1:
raise TypeError("expected 1D vector for x") # 检查x是否是一维数组
if x.size == 0:
raise TypeError("expected non-empty vector for x") # 检查x是否为空
if y.ndim < 1 or y.ndim > 2:
raise TypeError("expected 1D or 2D array for y") # 检查y是否是一维或二维数组
if len(x) != len(y):
raise TypeError("expected x and y to have same length") # 检查x和y的长度是否相同
if deg.ndim == 0:
lmax = deg
order = lmax + 1
van = vander_f(x, lmax) # 调用vander_f函数生成Vandermonde矩阵
else:
deg = np.sort(deg)
lmax = deg[-1]
order = len(deg)
van = vander_f(x, lmax)[:, deg] # 对于多项式阶数按deg给出的顺序进行排序，并生成对应的Vandermonde矩阵
# set up the least squares matrices in transposed form
lhs = van.T # 左侧矩阵为Vandermonde矩阵的转置
rhs = y.T # 右侧矩阵为y的转置
if w is not None:
w = np.asarray(w) + 0.0 # 将权重w转换为NumPy数组，并确保是浮点型
if w.ndim != 1:
raise TypeError("expected 1D vector for w") # 检查w是否是一维数组
if len(x) != len(w):
raise TypeError("expected x and w to have same length") # 检查x和w的长度是否相同
# apply weights. Don't use inplace operations as they
# can cause problems with NA.
lhs = lhs * w # 对左侧矩阵应用权重w
rhs = rhs * w # 对右侧矩阵应用权重w
# set rcond
if rcond is None:
rcond = len(x) * np.finfo(x.dtype).eps # 如果rcond未指定，则设为x数组中元素类型的机器精度的倍数
# 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 # 将范数为0的元素设为1，避免除以0
# Solve the least squares problem.
c, resids, rank, s = np.linalg.lstsq(lhs.T / scl, rhs.T, rcond) # 解最小二乘问题，得到系数c
c = (c.T / scl).T # 将系数c按照左侧矩阵列的范数进行归一化
# Expand c to include non-fitted coefficients which are set to zero
# 检查 deg 的维度是否大于 0
if deg.ndim > 0:
# 如果 c 的维度为 2，则创建一个零矩阵 cc，形状为 (lmax+1, c.shape[1])，数据类型与 c 相同
if c.ndim == 2:
cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype)
else:
# 否则，创建一个零数组 cc，形状为 (lmax+1)，数据类型与 c 相同
cc = np.zeros(lmax+1, dtype=c.dtype)
# 将 cc 中索引为 deg 的位置赋值为 c
cc[deg] = c
# 将 c 更新为 cc
c = cc
# 在秩降低时发出警告
if rank != order and not full:
# 创建警告信息字符串
msg = "The fit may be poorly conditioned"
# 发出警告，使用 RankWarning 类，设置堆栈深度为 2
warnings.warn(msg, RankWarning, stacklevel=2)
# 如果设置为返回完整结果
if full:
# 返回 c 和包含剩余项、秩、奇异值、rcond 值的列表
return c, [resids, rank, s, rcond]
else:
# 否则，只返回 c
return c
# 实现 ``<type>pow`` 函数的辅助函数，用于计算给定系列的幂次方
def _pow(mul_f, c, pow, maxpower):
"""
Helper function used to implement the ``<type>pow`` functions.
Parameters
----------
mul_f : function(array_like, array_like) -> ndarray
用于乘法操作的函数，如 ``polymul``
c : array_like
系列系数的一维数组
pow, maxpower
``<type>pow`` 函数的参数，详见其具体说明
"""
# c 是裁剪后的副本
[c] = as_series([c])
power = int(pow)
# 检查幂次必须为非负整数
if power != pow or power < 0:
raise ValueError("Power must be a non-negative integer.")
# 如果 maxpower 不为 None，则检查幂次不能超过 maxpower
elif maxpower is not None and power > maxpower:
raise ValueError("Power is too large")
elif power == 0:
# 幂次为 0 时，返回系数为 1 的数组
return np.array([1], dtype=c.dtype)
elif power == 1:
# 幂次为 1 时，直接返回系列 c
return c
else:
# 否则，通过循环计算幂次的结果
# 可以通过使用二进制的幂次计算来提高效率
prd = c
for i in range(2, power + 1):
prd = mul_f(prd, c)
return prd
# 类似于 `operator.index`，但当传入不正确类型时会抛出自定义异常的函数
def _as_int(x, desc):
"""
Like `operator.index`, but emits a custom exception when passed an
incorrect type
Parameters
----------
x : int-like
要解释为整数的值
desc : str
错误消息中包含的描述信息
Raises
------
TypeError : 若 x 是浮点数或非数值类型
"""
try:
return operator.index(x)
except TypeError as e:
raise TypeError(f"{desc} must be an integer, received {x}") from e
# 格式化浮点数 x，可选择是否使用括号包裹
def format_float(x, parens=False):
if not np.issubdtype(type(x), np.floating):
# 若 x 不是浮点数类型，则直接返回其字符串表示
return str(x)
opts = np.get_printoptions()
if np.isnan(x):
# 若 x 是 NaN，则返回指定的 NaN 字符串
return opts['nanstr']
elif np.isinf(x):
# 若 x 是无穷大，则返回指定的无穷大字符串
return opts['infstr']
exp_format = False
if x != 0:
a = np.abs(x)
# 根据浮点数的大小和精度设置，决定是否使用科学计数法
if a >= 1.e8 or a < 10**min(0, -(opts['precision']-1)//2):
exp_format = True
trim, unique = '0', True
if opts['floatmode'] == 'fixed':
trim, unique = 'k', False
if exp_format:
# 若使用科学计数法，则调用 dragon4_scientific 函数进行格式化
s = dragon4_scientific(x, precision=opts['precision'],
unique=unique, trim=trim,
sign=opts['sign'] == '+')
if parens:
# 若需要用括号包裹，则添加括号
s = '(' + s + ')'
else:
# 若使用定点表示法，则调用 dragon4_positional 函数进行格式化
s = dragon4_positional(x, precision=opts['precision'],
fractional=True,
unique=unique, trim=trim,
sign=opts['sign'] == '+')
return s
.\numpy\numpy\polynomial\polyutils.pyi
__all__: list[str]
def trimseq(seq):
...
def as_series(alist, trim=...):
...
def trimcoef(c, tol=...):
...
def getdomain(x):
...
def mapparms(old, new):
...
def mapdomain(x, old, new):
...
def format_float(x, parens=...):
...
.\numpy\numpy\polynomial\tests\test_chebyshev.py
"""Tests for chebyshev module.
"""
from functools import reduce
import numpy as np
import numpy.polynomial.chebyshev as cheb
from numpy.polynomial.polynomial import polyval
from numpy.testing import (
assert_almost_equal, assert_raises, assert_equal, assert_,
)
def trim(x):
return cheb.chebtrim(x, tol=1e-6)
T0 = [1]
T1 = [0, 1]
T2 = [-1, 0, 2]
T3 = [0, -3, 0, 4]
T4 = [1, 0, -8, 0, 8]
T5 = [0, 5, 0, -20, 0, 16]
T6 = [-1, 0, 18, 0, -48, 0, 32]
T7 = [0, -7, 0, 56, 0, -112, 0, 64]
T8 = [1, 0, -32, 0, 160, 0, -256, 0, 128]
T9 = [0, 9, 0, -120, 0, 432, 0, -576, 0, 256]
Tlist = [T0, T1, T2, T3, T4, T5, T6, T7, T8, T9]
class TestPrivate:
def test__cseries_to_zseries(self):
for i in range(5):
inp = np.array([2] + [1]*i, np.double)
tgt = np.array([.5]*i + [2] + [.5]*i, np.double)
res = cheb._cseries_to_zseries(inp)
assert_equal(res, tgt)
def test__zseries_to_cseries(self):
for i in range(5):
inp = np.array([.5]*i + [2] + [.5]*i, np.double)
tgt = np.array([2] + [1]*i, np.double)
res = cheb._zseries_to_cseries(inp)
assert_equal(res, tgt)
class TestConstants:
def test_chebdomain(self):
assert_equal(cheb.chebdomain, [-1, 1])
def test_chebzero(self):
assert_equal(cheb.chebzero, [0])
def test_chebone(self):
assert_equal(cheb.chebone, [1])
def test_chebx(self):
assert_equal(cheb.chebx, [0, 1])
class TestArithmetic:
def test_chebadd(self):
for i in range(5):
for j in range(5):
msg = f"At i={i}, j={j}"
tgt = np.zeros(max(i, j) + 1)
tgt[i] += 1
tgt[j] += 1
res = cheb.chebadd([0]*i + [1], [0]*j + [1])
assert_equal(trim(res), trim(tgt), err_msg=msg)
def test_chebsub(self):
for i in range(5):
for j in range(5):
msg = f"At i={i}, j={j}"
tgt = np.zeros(max(i, j) + 1)
tgt[i] += 1
tgt[j] -= 1
res = cheb.chebsub([0]*i + [1], [0]*j + [1])
assert_equal(trim(res), trim(tgt), err_msg=msg)
def test_chebmulx(self):
assert_equal(cheb.chebmulx([0]), [0])
assert_equal(cheb.chebmulx([1]), [0, 1])
for i in range(1, 5):
ser = [0]*i + [1]
tgt = [0]*(i - 1) + [.5, 0, .5]
assert_equal(cheb.chebmulx(ser), tgt)
def test_chebmul(self):
for i in range(5):
for j in range(5):
msg = f"At i={i}, j={j}"
tgt = np.zeros(i + j + 1)
tgt[i + j] += .5
tgt[abs(i - j)] += .5
res = cheb.chebmul([0]*i + [1], [0]*j + [1])
assert_equal(trim(res), trim(tgt), err_msg=msg)
def test_chebdiv(self):
for i in range(5):
for j in range(5):
msg = f"At i={i}, j={j}"
ci = [0]*i + [1]
cj = [0]*j + [1]
tgt = cheb.chebadd(ci, cj)
quo, rem = cheb.chebdiv(tgt, ci)
res = cheb.chebadd(cheb.chebmul(quo, ci), rem)
assert_equal(trim(res), trim(tgt), err_msg=msg)
def test_chebpow(self):
for i in range(5):
for j in range(5):
msg = f"At i={i}, j={j}"
c = np.arange(i + 1)
tgt = reduce(cheb.chebmul, [c]*j, np.array([1]))
res = cheb.chebpow(c, j)
assert_equal(trim(res), trim(tgt), err_msg=msg)
class TestEvaluation:
c1d = np.array([2.5, 2., 1.5])
c2d = np.einsum('i,j->ij', c1d, c1d)
c3d = np.einsum('i,j,k->ijk', c1d, c1d, c1d)
x = np.random.random((3, 5))*2 - 1
y = polyval(x, [1., 2., 3.])
def test_chebval(self):
assert_equal(cheb.chebval([], [1]).size, 0)
x = np.linspace(-1, 1)
y = [polyval(x, c) for c in Tlist]
for i in range(10):
msg = f"At i={i}"
tgt = y[i]
res = cheb.chebval(x, [0]*i + [1])
assert_almost_equal(res, tgt, err_msg=msg)
for i in range(3):
dims = [2]*i
x = np.zeros(dims)
assert_equal(cheb.chebval(x, [1]).shape, dims)
assert_equal(cheb.chebval(x, [1, 0]).shape, dims)
assert_equal(cheb.chebval(x, [1, 0, 0]).shape, dims)
def test_chebval2d(self):
x1, x2, x3 = self.x
y1, y2, y3 = self.y
assert_raises(ValueError, cheb.chebval2d, x1, x2[:2], self.c2d)
tgt = y1*y2
res = cheb.chebval2d(x1, x2, self.c2d)
assert_almost_equal(res, tgt)
z = np.ones((2, 3))
res = cheb.chebval2d(z, z, self.c2d)
assert_(res.shape == (2, 3))
def test_chebval3d(self):
x1, x2, x3 = self.x
y1, y2, y3 = self.y
assert_raises(ValueError, cheb.chebval3d, x1, x2, x3[:2], self.c3d)
tgt = y1*y2*y3
res = cheb.chebval3d(x1, x2, x3, self.c3d)
assert_almost_equal(res, tgt)
z = np.ones((2, 3))
res = cheb.chebval3d(z, z, z, self.c3d)
assert_(res.shape == (2, 3))
def test_chebgrid2d(self):
x1, x2, x3 = self.x
y1, y2, y3 = self.y
tgt = np.einsum('i,j->ij', y1, y2)
res = cheb.chebgrid2d(x1, x2, self.c2d)
assert_almost_equal(res, tgt)
z = np.ones((2, 3))
res = cheb.chebgrid2d(z, z, self.c2d)
assert_(res.shape == (2, 3)*2)
def test_chebgrid3d(self):
x1, x2, x3 = self.x
y1, y2, y3 = self.y
tgt = np.einsum('i,j,k->ijk', y1, y2, y3)
res = cheb.chebgrid3d(x1, x2, x3, self.c3d)
assert_almost_equal(res, tgt)
z = np.ones((2, 3))
res = cheb.chebgrid3d(z, z, z, self.c3d)
assert_(res.shape == (2, 3)*3)
class TestIntegral:
def test_chebint_axis(self):
c2d = np.random.random((3, 4))
tgt = np.vstack([cheb.chebint(c) for c in c2d.T]).T
res = cheb.chebint(c2d, axis=0)
assert_almost_equal(res, tgt)
tgt = np.vstack([cheb.chebint(c) for c in c2d])
res = cheb.chebint(c2d, axis=1)
assert_almost_equal(res, tgt)
tgt = np.vstack([cheb.chebint(c, k=3) for c in c2d])
res = cheb.chebint(c2d, k=3, axis=1)
assert_almost_equal(res, tgt)
class TestDerivative:
def test_chebder(self):
assert_raises(TypeError, cheb.chebder, [0], .5)
assert_raises(ValueError, cheb.chebder, [0], -1)
for i in range(5):
tgt = [0]*i + [1]
res = cheb.chebder(tgt, m=0)
assert_equal(trim(res), trim(tgt))
for i in range(5):
for j in range(2, 5):
tgt = [0]*i + [1]
res = cheb.chebder(cheb.chebint(tgt, m=j), m=j)
assert_almost_equal(trim(res), trim(tgt))
for i in range(5):
for j in range(2, 5):
tgt = [0]*i + [1]
res = cheb.chebder(cheb.chebint(tgt, m=j, scl=2), m=j, scl=.5)
assert_almost_equal(trim(res), trim(tgt))
def test_chebder_axis(self):
c2d = np.random.random((3, 4))
tgt = np.vstack([cheb.chebder(c) for c in c2d.T]).T
res = cheb.chebder(c2d, axis=0)
assert_almost_equal(res, tgt)
tgt = np.vstack([cheb.chebder(c) for c in c2d])
res = cheb.chebder(c2d, axis=1)
assert_almost_equal(res, tgt)
class TestVander:
x = np.random.random((3, 5))*2 - 1
def test_chebvander(self):
x = np.arange(3)
v = cheb.chebvander(x, 3)
assert_(v.shape == (3, 4))
for i in range(4):
coef = [0]*i + [1]
assert_almost_equal(v[..., i], cheb.chebval(x, coef))
x = np.array([[1, 2], [3, 4], [5, 6]])
v = cheb.chebvander(x, 3)
assert_(v.shape == (3, 2, 4))
for i in range(4):
coef = [0]*i + [1]
assert_almost_equal(v[..., i], cheb.chebval(x, coef))
def test_chebvander2d(self):
x1, x2, x3 = self.x
c = np.random.random((2, 3))
van = cheb.chebvander2d(x1, x2, [1, 2])
tgt = cheb.chebval2d(x1, x2, c)
res = np.dot(van, c.flat)
assert_almost_equal(res, tgt)
van = cheb.chebvander2d([x1], [x2], [1, 2])
assert_(van.shape == (1, 5, 6))
def test_chebvander3d(self):
x1, x2, x3 = self.x
c = np.random.random((2, 3, 4))
van = cheb.chebvander3d(x1, x2, x3, [1, 2, 3])
tgt = cheb.chebval3d(x1, x2, x3, c)
res = np.dot(van, c.flat)
assert_almost_equal(res, tgt)
van = cheb.chebvander3d([x1], [x2], [x3], [1, 2, 3])
assert_(van.shape == (1, 5, 24))
class TestFitting:
class TestInterpolate:
def f(self, x):
return x * (x - 1) * (x - 2)
def test_raises(self):
assert_raises(ValueError, cheb.chebinterpolate, self.f, -1)
assert_raises(TypeError, cheb.chebinterpolate, self.f, 10.)
def test_dimensions(self):
for deg in range(1, 5):
assert_(cheb.chebinterpolate(self.f, deg).shape == (deg + 1,))
def test_approximation(self):
def powx(x, p):
return x**p
x = np.linspace(-1, 1, 10)
for deg in range(0, 10):
for p in range(0, deg + 1):
c = cheb.chebinterpolate(powx, deg, (p,))
assert_almost_equal(cheb.chebval(x, c), powx(x, p), decimal=12)
class TestCompanion:
def test_raises(self):
assert_raises(ValueError, cheb.chebcompanion, [])
assert_raises(ValueError, cheb.chebcompanion, [1])
def test_dimensions(self):
for i in range(1, 5):
coef = [0]*i + [1]
assert_(cheb.chebcompanion(coef).shape == (i, i))
def test_linear_root(self):
assert_(cheb.chebcompanion([1, 2])[0, 0] == -.5)
class TestGauss:
def test_100(self):
x, w = cheb.chebgauss(100)
v = cheb.chebvander(x, 99)
vv = np.dot(v.T * w, v)
vd = 1/np.sqrt(vv.diagonal())
vv = vd[:, None] * vv * vd
assert_almost_equal(vv, np.eye(100))
tgt = np.pi
assert_almost_equal(w.sum(), tgt)
class TestMisc:
def test_chebfromroots(self):
res = cheb.chebfromroots([])
assert_almost_equal(trim(res), [1])
for i in range(1, 5):
roots = np.cos(np.linspace(-np.pi, 0, 2*i + 1)[1::2])
tgt = [0]*i + [1]
res = cheb.chebfromroots(roots)*2**(i-1)
assert_almost_equal(trim(res), trim(tgt))
def test_chebroots(self):
assert_almost_equal(cheb.chebroots([1]), [])
assert_almost_equal(cheb.chebroots([1, 2]), [-.5])
for i in range(2, 5):
tgt = np.linspace(-1, 1, i)
res = cheb.chebroots(cheb.chebfromroots(tgt))
assert_almost_equal(trim(res), trim(tgt))
def test_chebtrim(self):
coef = [2, -1, 1, 0]
assert_raises(ValueError, cheb.chebtrim, coef, -1)
assert_equal(cheb.chebtrim(coef), coef[:-1])
assert_equal(cheb.chebtrim(coef, 1), coef[:-3])
assert_equal(cheb.chebtrim(coef, 2), [0])
def test_chebline(self):
assert_equal(cheb.chebline(3, 4), [3, 4])
def test_cheb2poly(self):
for i in range(10):
assert_almost_equal(cheb.cheb2poly([0]*i + [1]), Tlist[i])
def test_poly2cheb(self):
for i in range(10):
assert_almost_equal(cheb.poly2cheb(Tlist[i]), [0]*i + [1])
def test_weight(self):
x = np.linspace(-1, 1, 11)[1:-1]
tgt = 1./(np.sqrt(1 + x) * np.sqrt(1 - x))
res = cheb.chebweight(x)
assert_almost_equal(res, tgt)
def test_chebpts1(self):
assert_raises(ValueError, cheb.chebpts1, 1.5)
assert_raises(ValueError, cheb.chebpts1, 0)
tgt = [0]
assert_almost_equal(cheb.chebpts1(1), tgt)
tgt = [-0.70710678118654746, 0.70710678118654746]
assert_almost_equal(cheb.chebpts1(2), tgt)
tgt = [-0.86602540378443871, 0, 0.86602540378443871]
assert_almost_equal(cheb.chebpts1(3), tgt)
tgt = [-0.9238795325, -0.3826834323, 0.3826834323, 0.9238795325]
assert_almost_equal(cheb.chebpts1(4), tgt)
def test_chebpts2(self):
assert_raises(ValueError, cheb.chebpts2, 1.5)
assert_raises(ValueError, cheb.chebpts2, 1)
tgt = [-1, 1]
assert_almost_equal(cheb.chebpts2(2), tgt)
tgt = [-1, 0, 1]
assert_almost_equal(cheb.chebpts2(3), tgt)
tgt = [-1, -0.5, .5, 1]
assert_almost_equal(cheb.chebpts2(4), tgt)
tgt = [-1.0, -0.707106781187, 0, 0.707106781187, 1.0]
assert_almost_equal(cheb.chebpts2(5), tgt)
.\numpy\numpy\polynomial\tests\test_classes.py
"""Test inter-conversion of different polynomial classes.
This tests the convert and cast methods of all the polynomial classes.
"""
import operator as op
from numbers import Number
import pytest
import numpy as np
from numpy.polynomial import (
Polynomial, Legendre, Chebyshev, Laguerre, Hermite, HermiteE)
from numpy.testing import (
assert_almost_equal, assert_raises, assert_equal, assert_,
)
from numpy.exceptions import RankWarning
classes = (
Polynomial, Legendre, Chebyshev, Laguerre,
Hermite, HermiteE
)
classids = tuple(cls.__name__ for cls in classes)
@pytest.fixture(params=classes, ids=classids)
def Poly(request):
return request.param
random = np.random.random
def assert_poly_almost_equal(p1, p2, msg=""):
try:
assert_(np.all(p1.domain == p2.domain))
assert_(np.all(p1.window == p2.window))
assert_almost_equal(p1.coef, p2.coef)
except AssertionError:
msg = f"Result: {p1}\nTarget: {p2}"
raise AssertionError(msg)
Poly1 = Poly
Poly2 = Poly
def test_conversion(Poly1, Poly2):
x = np.linspace(0, 1, 10)
coef = random((3,))
d1 = Poly1.domain + random((2,))*.25
w1 = Poly1.window + random((2,))*.25
p1 = Poly1(coef, domain=d1, window=w1)
d2 = Poly2.domain + random((2,))*.25
w2 = Poly2.window + random((2,))*.25
p2 = p1.convert(kind=Poly2, domain=d2, window=w2)
assert_almost_equal(p2.domain, d2)
assert_almost_equal(p2.window, w2)
assert_almost_equal(p2(x), p1(x))
def test_cast(Poly1, Poly2):
x = np.linspace(0, 1, 10)
coef = random((3,))
d1 = Poly1.domain + random((2,))*.25
w1 = Poly1.window + random((2,))*.25
p1 = Poly1(coef, domain=d1, window=w1)
d2 = Poly2.domain + random((2,))*.25
w2 = Poly2.window + random((2,))*.25
p2 = Poly2.cast(p1, domain=d2, window=w2)
assert_almost_equal(p2.domain, d2)
assert_almost_equal(p2.window, w2)
assert_almost_equal(p2(x), p1(x))
def test_identity(Poly):
d = Poly.domain + random((2,))*.25
w = Poly.window + random((2,))*.25
x = np.linspace(d[0], d[1], 11)
p = Poly.identity(domain=d, window=w)
assert_equal(p.domain, d)
assert_equal(p.window, w)
assert_almost_equal(p(x), x)
def test_basis(Poly):
d = Poly.domain + random((2,))*.25
w = Poly.window + random((2,))*.25
p = Poly.basis(5, domain=d, window=w)
assert_equal(p.domain, d)
assert_equal(p.window, w)
assert_equal(p.coef, [0]*5 + [1])
def test_fromroots(Poly):
d = Poly.domain + random((2,))*.25
w = Poly.window + random((2,))*.25
r = random((5,))
p1 = Poly.fromroots(r, domain=d, window=w)
assert_equal(p1.degree(), len(r))
assert_equal(p1.domain, d)
assert_equal(p1.window, w)
assert_almost_equal(p1(r), 0)
pdom = Polynomial.domain
pwin = Polynomial.window
p2 = Polynomial.cast(p1, domain=pdom, window=pwin)
assert_almost_equal(p2.coef[-1], 1)
def test_bad_conditioned_fit(Poly):
x = [0., 0., 1.]
y = [1., 2., 3.]
with pytest.warns(RankWarning) as record:
Poly.fit(x, y, 2)
assert record[0].message.args[0] == "The fit may be poorly conditioned"
def test_fit(Poly):
def f(x):
return x*(x - 1)*(x - 2)
x = np.linspace(0, 3)
y = f(x)
p = Poly.fit(x, y, 3)
assert_almost_equal(p.domain, [0, 3])
assert_almost_equal(p(x), y)
assert_equal(p.degree(), 3)
d = Poly.domain + random((2,))*.25
w = Poly.window + random((2,))*.25
p = Poly.fit(x, y, 3, domain=d, window=w)
assert_almost_equal(p(x), y)
assert_almost_equal(p.domain, d)
assert_almost_equal(p.window, w)
p = Poly.fit(x, y, [0, 1, 2, 3], domain=d, window=w)
assert_almost_equal(p(x), y)
assert_almost_equal(p.domain, d)
assert_almost_equal(p.window, w)
p = Poly.fit(x, y, 3, [])
assert_equal(p.domain, Poly.domain)
assert_equal(p.window, Poly.window)
p = Poly.fit(x, y, [0, 1, 2, 3], [])
assert_equal(p.domain, Poly.domain)
assert_equal(p.window, Poly.window)
w = np.zeros_like(x)
z = y + random(y.shape)*.25
w[::2] = 1
p1 = Poly.fit(x[::2], z[::2], 3)
p2 = Poly.fit(x, z, 3, w=w)
p3 = Poly.fit(x, z, [0, 1, 2, 3], w=w)
assert_almost_equal(p1(x), p2(x))
assert_almost_equal(p2(x), p3(x))
def test_equal(Poly):
p1 = Poly([1, 2, 3], domain=[0, 1], window=[2, 3])
p2 = Poly([1, 1, 1], domain=[0, 1], window=[2, 3])
p3 = Poly([1, 2, 3], domain=[1, 2], window=[2, 3])
p4 = Poly([1, 2, 3], domain=[0, 1], window=[1, 2])
assert_(p1 == p1)
assert_(not p1 == p2)
assert_(not p1 == p3)
assert_(not p1 == p4)
def test_not_equal(Poly):
p1 = Poly([1, 2, 3], domain=[0, 1], window=[2, 3])
p2 = Poly([1, 1, 1], domain=[0, 1], window=[2, 3])
p3 = Poly([1, 2, 3], domain=[1, 2], window=[2, 3])
p4 = Poly([1, 2, 3], domain=[0, 1], window=[1, 2])
assert_(not p1 != p1)
assert_(p1 != p2)
assert_(p1 != p3)
assert_(p1 != p4)
def test_add(Poly):
c1 = list(random((4,)) + .5)
c2 = list(random((3,)) + .5)
p1 = Poly(c1)
p2 = Poly(c2)
p3 = p1 + p2
assert_poly_almost_equal(p2 + p1, p3)
assert_poly_almost_equal(p1 + c2, p3)
assert_poly_almost_equal(c2 + p1, p3)
assert_poly_almost_equal(p1 + tuple(c2), p3)
assert_poly_almost_equal(tuple(c2) + p1, p3)
assert_poly_almost_equal(p1 + np.array(c2), p3)
assert_poly_almost_equal(np.array(c2) + p1, p3)
assert_raises(TypeError, op.add, p1, Poly([0], domain=Poly.domain + 1))
assert_raises(TypeError, op.add, p1, Poly([0], window=Poly.window + 1))
if Poly is Polynomial:
assert_raises(TypeError, op.add, p1, Chebyshev([0]))
else:
assert_raises(TypeError, op.add, p1, Polynomial([0]))
def test_sub(Poly):
c1 = list(random((4,)) + .5)
c2 = list(random((3,)) + .5)
p1 = Poly(c1)
p2 = Poly(c2)
p3 = p1 - p2
assert_poly_almost_equal(p2 - p1, -p3)
assert_poly_almost_equal(p1 - c2, p3)
assert_poly_almost_equal(c2 - p1, -p3)
assert_poly_almost_equal(p1 - tuple(c2), p3)
assert_poly_almost_equal(tuple(c2) - p1, -p3)
assert_poly_almost_equal(p1 - np.array(c2), p3)
assert_poly_almost_equal(np.array(c2) - p1, -p3)
assert_raises(TypeError, op.sub, p1, Poly([0], domain=Poly.domain + 1))
assert_raises(TypeError, op.sub, p1, Poly([0], window=Poly.window + 1))
if Poly is Polynomial:
assert_raises(TypeError, op.sub, p1, Chebyshev([0]))
else:
assert_raises(TypeError, op.sub, p1, Polynomial([0]))
def test_mul(Poly):
c1 = list(random((4,)) + .5)
c2 = list(random((3,)) + .5)
p1 = Poly(c1)
p2 = Poly(c2)
p3 = p1 * p2
assert_poly_almost_equal(p2 * p1, p3)
assert_poly_almost_equal(p1 * c2, p3)
assert_poly_almost_equal(c2 * p1, p3)
assert_poly_almost_equal(p1 * tuple(c2), p3)
assert_poly_almost_equal(tuple(c2) * p1, p3)
assert_poly_almost_equal(p1 * np.array(c2), p3)
assert_poly_almost_equal(np.array(c2) * p1, p3)
assert_poly_almost_equal(p1 * 2, p1 * Poly([2]))
assert_poly_almost_equal(2 * p1, p1 * Poly([2]))
assert_raises(TypeError, op.mul, p1, Poly([0], domain=Poly.domain + 1))
assert_raises(TypeError, op.mul, p1, Poly([0], window=Poly.window + 1))
if Poly is Polynomial:
assert_raises(TypeError, op.mul, p1, Chebyshev([0]))
else:
assert_raises(TypeError, op.mul, p1, Polynomial([0]))
def test_floordiv(Poly):
c1 = list(random((4,)) + .5)
c2 = list(random((3,)) + .5)
c3 = list(random((2,)) + .5)
p1 = Poly(c1)
p2 = Poly(c2)
p3 = Poly(c3)
p4 = p1 * p2 + p3
c4 = list(p4.coef)
assert_poly_almost_equal(p4 // p2, p1)
assert_poly_almost_equal(p4 // c2, p1)
assert_poly_almost_equal(c4 // p2, p1)
assert_poly_almost_equal(p4 // tuple(c2), p1)
assert_poly_almost_equal(tuple(c4) // p2, p1)
assert_poly_almost_equal(p4 // np.array(c2), p1)
assert_poly_almost_equal(np.array(c4) // p2, p1)
assert_poly_almost_equal(2 // p2, Poly([0]))
assert_poly_almost_equal(p2 // 2, 0.5 * p2)
assert_raises(
TypeError, op.floordiv, p1, Poly([0], domain=Poly.domain + 1))
assert_raises(
TypeError, op.floordiv, p1, Poly([0], window=Poly.window + 1))
if Poly is Polynomial:
assert_raises(TypeError, op.floordiv, p1, Chebyshev([0]))
else:
assert_raises(TypeError, op.floordiv, p1, Polynomial([0]))
def test_truediv(Poly):
p1 = Poly([1, 2, 3])
p2 = p1 * 5
for stype in np.ScalarType:
if not issubclass(stype, Number) or issubclass(stype, bool):
continue
s = stype(5)
assert_poly_almost_equal(op.truediv(p2, s), p1)
assert_raises(TypeError, op.truediv, s, p2)
for stype in (int, float):
s = stype(5)
assert_poly_almost_equal(op.truediv(p2, s), p1)
assert_raises(TypeError, op.truediv, s, p2)
for stype in [complex]:
s = stype(5, 0)
assert_poly_almost_equal(op.truediv(p2, s), p1)
assert_raises(TypeError, op.truediv, s, p2)
for s in [tuple(), list(), dict(), bool(), np.array([1])]:
assert_raises(TypeError, op.truediv, p2, s)
assert_raises(TypeError, op.truediv, s, p2)
for ptype in classes:
assert_raises(TypeError, op.truediv, p2, ptype(1))
def test_mod(Poly):
c1 = list(random((4,)) + .5)
c2 = list(random((3,)) + .5)
c3 = list(random((2,)) + .5)
p1 = Poly(c1)
p2 = Poly(c2)
p3 = Poly(c3)
p4 = p1 * p2 + p3
c4 = list(p4.coef)
assert_poly_almost_equal(p4 % p2, p3)
assert_poly_almost_equal(p4 % c2, p3)
assert_poly_almost_equal(c4 % p2, p3)
assert_poly_almost_equal(p4 % tuple(c2), p3)
assert_poly_almost_equal(tuple(c4) % p2, p3)
assert_poly_almost_equal(p4 % np.array(c2), p3)
assert_poly_almost_equal(np.array(c4) % p2, p3)
assert_poly_almost_equal(2 % p2, Poly([2]))
assert_poly_almost_equal(p2 % 2, Poly([0]))
assert_raises(TypeError, op.mod, p1, Poly([0], domain=Poly.domain + 1))
assert_raises(TypeError, op.mod, p1, Poly([0], window=Poly.window + 1))
if Poly is Polynomial:
assert_raises(TypeError, op.mod, p1, Chebyshev([0]))
else:
assert_raises(TypeError, op.mod, p1, Polynomial([0]))
def test_divmod(Poly):
c1 = list(random((4,)) + .5)
c2 = list(random((3,)) + .5)
c3 = list(random((2,)) + .5)
p1 = Poly(c1)
p2 = Poly(c2)
p3 = Poly(c3)
p4 = p1 * p2 + p3
c4 = list(p4.coef)
quo, rem = divmod(p4, p2)
assert_poly_almost_equal(quo, p1)
assert_poly_almost_equal(rem, p3)
quo, rem = divmod(p4, c2)
assert_poly_almost_equal(quo, p1)
assert_poly_almost_equal(rem, p3)
quo, rem = divmod(c4, p2)
assert_poly_almost_equal(quo, p1)
assert_poly_almost_equal(rem, p3)
quo, rem = divmod(p4, tuple(c2))
assert_poly_almost_equal(quo, p1)
assert_poly_almost_equal(rem, p3)
quo, rem = divmod(tuple(c4), p2)
assert_poly_almost_equal(quo, p1)
assert_poly_almost_equal(rem, p3)
quo, rem = divmod(p4, np.array(c2))
assert_poly_almost_equal(quo, p1)
assert_poly_almost_equal(rem, p3)
quo, rem = divmod(np.array(c4), p2)
assert_poly_almost_equal(quo, p1)
assert_poly_almost_equal(rem, p3)
quo, rem = divmod(p2, 2)
assert_poly_almost_equal(quo, 0.5*p2)
assert_poly_almost_equal(rem, Poly([0]))
quo, rem = divmod(2, p2)
assert_poly_almost_equal(quo, Poly([0]))
assert_poly_almost_equal(rem, Poly([2]))
assert_raises(TypeError, divmod, p1, Poly([0], domain=Poly.domain + 1))
assert_raises(TypeError, divmod, p1, Poly([0], window=Poly.window + 1))
if Poly is Polynomial:
assert_raises(TypeError, divmod, p1, Chebyshev([0]))
else:
assert_raises(TypeError, divmod, p1, Polynomial([0]))
def test_roots(Poly):
d = Poly.domain * 1.25 + .25
w = Poly.window
tgt = np.linspace(d[0], d[1], 5)
res = np.sort(Poly.fromroots(tgt, domain=d, window=w).roots())
assert_almost_equal(res, tgt)
res = np.sort(Poly.fromroots(tgt).roots())
assert_almost_equal(res, tgt)
def test_degree(Poly):
p = Poly.basis(5)
assert_equal(p.degree(), 5)
def test_copy(Poly):
p1 = Poly.basis(5)
p2 = p1.copy()
assert_(p1 == p2)
assert_(p1 is not p2)
assert_(p1.coef is not p2.coef)
assert_(p1.domain is not p2.domain)
assert_(p1.window is not p2.window)
def test_integ(Poly):
P = Polynomial
p0 = Poly.cast(P([1*2, 2*3, 3*4]))
p1 = P.cast(p0.integ())
p2 = P.cast(p0.integ(2))
assert_poly_almost_equal(p1, P([0, 2, 3, 4]))
assert_poly_almost_equal(p2, P([0, 0, 1, 1, 1]))
p0 = Poly.cast(P([1*2, 2*3, 3*4]))
p1 = P.cast(p0.integ(k=1))
p2 = P.cast(p0.integ(2, k=[1, 1]))
assert_poly_almost_equal(p1, P([1, 2, 3, 4]))
assert_poly_almost_equal(p2, P([1, 1, 1, 1, 1]))
p0 = Poly.cast(P([1*2, 2*3, 3*4]))
p1 = P.cast(p0.integ(lbnd=1))
p2 = P.cast(p0.integ(2, lbnd=1))
assert_poly_almost_equal(p1, P([-9, 2, 3, 4]))
assert_poly_almost_equal(p2, P([6, -9, 1, 1, 1]))
d = 2*Poly.domain
p0 = Poly.cast(P([1*2, 2*3, 3*4]), domain=d)
p1 = P.cast(p0.integ())
p2 = P.cast(p0.integ(2))
assert_poly_almost_equal(p1, P([0, 2, 3, 4]))
assert_poly_almost_equal(p2, P([0, 0, 1, 1, 1]))
def test_deriv(Poly):
d = Poly.domain + random((2,))*.25
w = Poly.window + random((2,))*.25
p1 = Poly([1, 2, 3], domain=d, window=w)
p2 = p1.integ(2, k=[1, 2])
p3 = p1.integ(1, k=[1])
assert_almost_equal(p2.deriv(1).coef, p3.coef)
assert_almost_equal(p2.deriv(2).coef, p1.coef)
p1 = Poly([1, 2, 3])
p2 = p1.integ(2, k=[1, 2])
p3 = p1.integ(1, k=[1])
assert_almost_equal(p2.deriv(1).coef, p3.coef)
assert_almost_equal(p2.deriv(2).coef, p1.coef)
def test_linspace(Poly):
d = Poly.domain + random((2,))*.25
w = Poly.window + random((2,))*.25
p = Poly([1, 2, 3], domain=d, window=w)
xtgt = np.linspace(d[0], d[1], 20)
ytgt = p(xtgt)
xres, yres = p.linspace(20)
assert_almost_equal(xres, xtgt)
assert_almost_equal(yres, ytgt)
xtgt = np.linspace(0, 2, 20)
ytgt = p(xtgt)
xres, yres = p.linspace(20, domain=[0, 2])
assert_almost_equal(xres, xtgt)
assert_almost_equal(yres, ytgt)
def test_pow(Poly):
d = Poly.domain + random((2,))*.25
w = Poly.window + random((2,))*.25
tgt = Poly([1], domain=d, window=w)
tst = Poly([1, 2, 3], domain=d, window=w)
for i in range(5):
assert_poly_almost_equal(tst**i, tgt)
tgt = tgt * tst
tgt = Poly([1])
tst = Poly([1, 2, 3])
for i in range(5):
assert_poly_almost_equal(tst**i, tgt)
tgt = tgt * tst
assert_raises(ValueError, op.pow, tgt, 1.5)
assert_raises(ValueError, op.pow, tgt, -1)
def test_call(Poly):
P = Polynomial
d = Poly.domain
x = np.linspace(d[0], d[1], 11)
p = Poly.cast(P([1, 2, 3]))
tgt = 1 + x*(2 + 3*x)
res = p(x)
assert_almost_equal(res, tgt)
def test_call_with_list(Poly):
p = Poly([1, 2, 3])
x = [-1, 0, 2]
res = p(x)
assert_equal(res, p(np.array(x)))
def test_cutdeg(Poly):
p = Poly([1, 2, 3])
assert_raises(ValueError, p.cutdeg, .5)
assert_raises(ValueError, p.cutdeg, -1)
assert_equal(len(p.cutdeg(3)), 3)
assert_equal(len(p.cutdeg(2)), 3)
assert_equal(len(p.cutdeg(1)), 2)
assert_equal(len(p.cutdeg(0)), 1)
def test_truncate(Poly):
p = Poly([1, 2, 3])
assert_raises(ValueError, p.truncate, .5)
assert_raises(ValueError, p.truncate, 0)
assert_equal(len(p.truncate(4)), 3)
assert_equal(len(p.truncate(3)), 3)
assert_equal(len(p.truncate(2)), 2)
assert_equal(len(p.truncate(1)), 1)
def test_trim(Poly):
c = [1, 1e-6, 1e-12, 0]
p = Poly(c)
assert_equal(p.trim().coef, c[:3])
assert_equal(p.trim(1e-10).coef, c[:2])
assert_equal(p.trim(1e-5).coef, c[:1])
def test_mapparms(Poly):
d = Poly.domain
w = Poly.window
p = Poly([1], domain=d, window=w)
assert_almost_equal([0, 1], p.mapparms())
w = 2*d + 1
p = Poly([1], domain=d, window=w)
assert_almost_equal([1, 2], p.mapparms())
def test_ufunc_override(Poly):
p = Poly([1, 2, 3])
x = np.ones(3)
assert_raises(TypeError, np.add, p, x)
assert_raises(TypeError, np.add, x, p)
class TestInterpolate:
def f(self, x):
return x * (x - 1) * (x - 2)
def test_raises(self):
assert_raises(ValueError, Chebyshev.interpolate, self.f, -1)
assert_raises(TypeError, Chebyshev.interpolate, self.f, 10.)
def test_dimensions(self):
for deg in range(1, 5):
assert_(Chebyshev.interpolate(self.f, deg).degree() == deg)
def test_approximation(self):
def powx(x, p):
return x**p
x = np.linspace(0, 2, 10)
for deg in range(0, 10):
for t in range(0, deg + 1):
p = Chebyshev.interpolate(powx, deg, domain=[0, 2], args=(t,))
