NumPy 源码解析（三十一）
.\numpy\numpy\ma\timer_comparison.py
import timeit
from functools import reduce
import numpy as np
import numpy._core.fromnumeric as fromnumeric
from numpy.testing import build_err_msg
pi = np.pi
class ModuleTester:
def __init__(self, module):
self.module = module
self.allequal = module.allequal
self.arange = module.arange
self.array = module.array
self.concatenate = module.concatenate
self.count = module.count
self.equal = module.equal
self.filled = module.filled
self.getmask = module.getmask
self.getmaskarray = module.getmaskarray
self.id = id
self.inner = module.inner
self.make_mask = module.make_mask
self.masked = module.masked
self.masked_array = module.masked_array
self.masked_values = module.masked_values
self.mask_or = module.mask_or
self.nomask = module.nomask
self.ones = module.ones
self.outer = module.outer
self.repeat = module.repeat
self.resize = module.resize
self.sort = module.sort
self.take = module.take
self.transpose = module.transpose
self.zeros = module.zeros
self.MaskType = module.MaskType
try:
self.umath = module.umath
except AttributeError:
self.umath = module.core.umath
self.testnames = []
def assert_array_compare(self, comparison, x, y, err_msg='', header='',
fill_value=True):
"""
Assert that a comparison of two masked arrays is satisfied elementwise.
"""
xf = self.filled(x)
yf = self.filled(y)
m = self.mask_or(self.getmask(x), self.getmask(y))
x = self.filled(self.masked_array(xf, mask=m), fill_value)
y = self.filled(self.masked_array(yf, mask=m), fill_value)
if (x.dtype.char != "O"):
x = x.astype(np.float64)
if isinstance(x, np.ndarray) and x.size > 1:
x[np.isnan(x)] = 0
elif np.isnan(x):
x = 0
if (y.dtype.char != "O"):
y = y.astype(np.float64)
if isinstance(y, np.ndarray) and y.size > 1:
y[np.isnan(y)] = 0
elif np.isnan(y):
y = 0
try:
cond = (x.shape == () or y.shape == ()) or x.shape == y.shape
if not cond:
msg = build_err_msg([x, y],
err_msg
+ f'\n(shapes {x.shape}, {y.shape} mismatch)',
header=header,
names=('x', 'y'))
assert cond, msg
val = comparison(x, y)
if m is not self.nomask and fill_value:
val = self.masked_array(val, mask=m)
if isinstance(val, bool):
cond = val
reduced = [0]
else:
reduced = val.ravel()
cond = reduced.all()
reduced = reduced.tolist()
if not cond:
match = 100-100.0*reduced.count(1)/len(reduced)
msg = build_err_msg([x, y],
err_msg
+ '\n(mismatch %s%%)' % (match,),
header=header,
names=('x', 'y'))
assert cond, msg
except ValueError as e:
msg = build_err_msg([x, y], err_msg, header=header, names=('x', 'y'))
raise ValueError(msg) from e
def assert_array_equal(self, x, y, err_msg=''):
"""
Checks the elementwise equality of two masked arrays.
"""
self.assert_array_compare(self.equal, x, y, err_msg=err_msg,
header='Arrays are not equal')
@np.errstate(all='ignore')
def test_0(self):
"""
Tests creation
"""
x = np.array([1., 1., 1., -2., pi/2.0, 4., 5., -10., 10., 1., 2., 3.])
m = [1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0]
xm = self.masked_array(x, mask=m)
xm[0]
@np.errstate(all='ignore')
@np.errstate(all='ignore')
def test_1(self):
"""
Tests creation
"""
x = np.array([1., 1., 1., -2., pi/2.0, 4., 5., -10., 10., 1., 2., 3.])
y = np.array([5., 0., 3., 2., -1., -4., 0., -10., 10., 1., 0., 3.])
m1 = [1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0]
m2 = [0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1]
xm = self.masked_array(x, mask=m1)
ym = self.masked_array(y, mask=m2)
xf = np.where(m1, 1.e+20, x)
xm.set_fill_value(1.e+20)
assert((xm-ym).filled(0).any())
s = x.shape
assert(xm.size == reduce(lambda x, y:x*y, s))
assert(self.count(xm) == len(m1) - reduce(lambda x, y:x+y, m1))
for s in [(4, 3), (6, 2)]:
x.shape = s
y.shape = s
xm.shape = s
ym.shape = s
xf.shape = s
assert(self.count(xm) == len(m1) - reduce(lambda x, y:x+y, m1))
@np.errstate(all='ignore')
def test_2(self):
"""
Tests conversions and indexing.
"""
x1 = np.array([1, 2, 4, 3])
x2 = self.array(x1, mask=[1, 0, 0, 0])
x3 = self.array(x1, mask=[0, 1, 0, 1])
x4 = self.array(x1)
str(x2)
repr(x2)
assert type(x2[1]) is type(x1[1])
assert x1[1] == x2[1]
x1[2] = 9
x2[2] = 9
self.assert_array_equal(x1, x2)
x1[1:3] = 99
x2[1:3] = 99
x2[1] = self.masked
x2[1:3] = self.masked
x2[:] = x1
x2[1] = self.masked
x3[:] = self.masked_array([1, 2, 3, 4], [0, 1, 1, 0])
x4[:] = self.masked_array([1, 2, 3, 4], [0, 1, 1, 0])
x1 = np.arange(5)*1.0
x2 = self.masked_values(x1, 3.0)
x1 = self.array([1, 'hello', 2, 3], object)
x2 = np.array([1, 'hello', 2, 3], object)
x1[1]
x2[1]
assert x1[1:1].shape == (0,)
n = [0, 0, 1, 0, 0]
m = self.make_mask(n)
m2 = self.make_mask(m)
assert(m is m2)
m3 = self.make_mask(m, copy=1)
assert(m is not m3)
@np.errstate(all='ignore')
def test_3(self):
"""
Tests resize/repeat
"""
x4 = self.arange(4)
x4[2] = self.masked
y4 = self.resize(x4, (8,))
assert self.allequal(self.concatenate([x4, x4]), y4)
assert self.allequal(self.getmask(y4), [0, 0, 1, 0, 0, 0, 1, 0])
y5 = self.repeat(x4, (2, 2, 2, 2), axis=0)
self.assert_array_equal(y5, [0, 0, 1, 1, 2, 2, 3, 3])
y6 = self.repeat(x4, 2, axis=0)
assert self.allequal(y5, y6)
y7 = x4.repeat((2, 2, 2, 2), axis=0)
assert self.allequal(y5, y7)
y8 = x4.repeat(2, 0)
assert self.allequal(y5, y8)
def test_4(self):
"""
Test of take, transpose, inner, outer products.
测试 take、transpose、inner、outer 函数的功能。
"""
x = self.arange(24)
y = np.arange(24)
x[5:6] = self.masked
x = x.reshape(2, 3, 4)
y = y.reshape(2, 3, 4)
assert self.allequal(np.transpose(y, (2, 0, 1)), self.transpose(x, (2, 0, 1)))
assert self.allequal(np.take(y, (2, 0, 1), 1), self.take(x, (2, 0, 1), 1))
assert self.allequal(np.inner(self.filled(x, 0), self.filled(y, 0)),
self.inner(x, y))
assert self.allequal(np.outer(self.filled(x, 0), self.filled(y, 0)),
self.outer(x, y))
y = self.array(['abc', 1, 'def', 2, 3], object)
y[2] = self.masked
t = self.take(y, [0, 3, 4])
assert t[0] == 'abc'
assert t[1] == 2
assert t[2] == 3
@np.errstate(all='ignore')
def test_5(self):
"""
Tests inplace w/ scalar
测试原地操作与标量。
"""
x = self.arange(10)
y = self.arange(10)
xm = self.arange(10)
xm[2] = self.masked
x += 1
assert self.allequal(x, y+1)
xm += 1
assert self.allequal(xm, y+1)
x = self.arange(10)
xm = self.arange(10)
xm[2] = self.masked
x -= 1
assert self.allequal(x, y-1)
xm -= 1
assert self.allequal(xm, y-1)
x = self.arange(10)*1.0
xm = self.arange(10)*1.0
xm[2] = self.masked
x *= 2.0
assert self.allequal(x, y*2)
xm *= 2.0
assert self.allequal(xm, y*2)
x = self.arange(10)*2
xm = self.arange(10)*2
xm[2] = self.masked
x /= 2
assert self.allequal(x, y)
xm /= 2
assert self.allequal(xm, y)
x = self.arange(10)*1.0
xm = self.arange(10)*1.0
xm[2] = self.masked
x /= 2.0
assert self.allequal(x, y/2.0)
xm /= self.arange(10)
self.assert_array_equal(xm, self.ones((10,)))
x = self.arange(10).astype(np.float64)
xm = self.arange(10)
xm[2] = self.masked
x += 1.
assert self.allequal(x, y + 1.)
x = self.arange(10, dtype=np.float64)
y = self.arange(10)
xm = self.arange(10, dtype=np.float64)
xm[2] = self.masked
m = xm.mask
a = self.arange(10, dtype=np.float64)
a[-1] = self.masked
x += a
xm += a
assert self.allequal(x, y+a)
assert self.allequal(xm, y+a)
assert self.allequal(xm.mask, self.mask_or(m, a.mask))
x = self.arange(10, dtype=np.float64)
xm = self.arange(10, dtype=np.float64)
xm[2] = self.masked
m = xm.mask
a = self.arange(10, dtype=np.float64)
a[-1] = self.masked
x -= a
xm -= a
assert self.allequal(x, y-a)
assert self.allequal(xm, y-a)
assert self.allequal(xm.mask, self.mask_or(m, a.mask))
x = self.arange(10, dtype=np.float64)
xm = self.arange(10, dtype=np.float64)
xm[2] = self.masked
m = xm.mask
a = self.arange(10, dtype=np.float64)
a[-1] = self.masked
x *= a
xm *= a
assert self.allequal(x, y*a)
assert self.allequal(xm, y*a)
assert self.allequal(xm.mask, self.mask_or(m, a.mask))
x = self.arange(10, dtype=np.float64)
xm = self.arange(10, dtype=np.float64)
xm[2] = self.masked
m = xm.mask
a = self.arange(10, dtype=np.float64)
a[-1] = self.masked
x /= a
xm /= a
for f in [
'sin', 'cos', 'tan',
'arcsin', 'arccos', 'arctan',
'sinh', 'cosh', 'tanh',
'arcsinh',
'arccosh',
'arctanh',
'absolute', 'fabs', 'negative',
'floor', 'ceil',
'logical_not',
'add', 'subtract', 'multiply',
'divide', 'true_divide', 'floor_divide',
'remainder', 'fmod', 'hypot', 'arctan2',
'equal', 'not_equal', 'less_equal', 'greater_equal',
'less', 'greater',
'logical_and', 'logical_or', 'logical_xor',
]:
try:
uf = getattr(self.umath, f)
except AttributeError:
uf = getattr(fromnumeric, f)
mf = getattr(self.module, f)
args = d[:uf.nin]
ur = uf(*args)
mr = mf(*args)
self.assert_array_equal(ur.filled(0), mr.filled(0), f)
self.assert_array_equal(ur._mask, mr._mask)
@np.errstate(all='ignore')
def test_99(self):
ott = self.array([0., 1., 2., 3.], mask=[1, 0, 0, 0])
self.assert_array_equal(2.0, self.average(ott, axis=0))
self.assert_array_equal(2.0, self.average(ott, weights=[1., 1., 2., 1.]))
result, wts = self.average(ott, weights=[1., 1., 2., 1.], returned=1)
self.assert_array_equal(2.0, result)
assert(wts == 4.0)
ott[:] = self.masked
assert(self.average(ott, axis=0) is self.masked)
ott = self.array([0., 1., 2., 3.], mask=[1, 0, 0, 0])
ott = ott.reshape(2, 2)
ott[:, 1] = self.masked
self.assert_array_equal(self.average(ott, axis=0), [2.0, 0.0])
assert(self.average(ott, axis=1)[0] is self.masked)
self.assert_array_equal([2., 0.], self.average(ott, axis=0))
result, wts = self.average(ott, axis=0, returned=1)
self.assert_array_equal(wts, [1., 0.])
w1 = [0, 1, 1, 1, 1, 0]
w2 = [[0, 1, 1, 1, 1, 0], [1, 0, 0, 0, 0, 1]]
x = self.arange(6)
self.assert_array_equal(self.average(x, axis=0), 2.5)
self.assert_array_equal(self.average(x, axis=0, weights=w1), 2.5)
y = self.array([self.arange(6), 2.0*self.arange(6)])
self.assert_array_equal(self.average(y, None), np.add.reduce(np.arange(6))*3./12.)
self.assert_array_equal(self.average(y, axis=0), np.arange(6) * 3./2.)
self.assert_array_equal(self.average(y, axis=1), [self.average(x, axis=0), self.average(x, axis=0) * 2.0])
self.assert_array_equal(self.average(y, None, weights=w2), 20./6.)
self.assert_array_equal(self.average(y, axis=0, weights=w2), [0., 1., 2., 3., 4., 10.])
m1 = self.zeros(6)
m2 = [0, 0, 1, 1, 0, 0]
m3 = [[0, 0, 1, 1, 0, 0], [0, 1, 1, 1, 1, 0]]
m4 = self.ones(6)
m5 = [0, 1, 1, 1, 1, 1]
self.assert_array_equal(self.average(self.masked_array(x, m1), axis=0), 2.5)
self.assert_array_equal(self.average(self.masked_array(x, m2), axis=0), 2.5)
self.assert_array_equal(self.average(self.masked_array(x, m5), axis=0), 0.0)
self.assert_array_equal(self.count(self.average(self.masked_array(x, m4), axis=0)), 0)
z = self.masked_array(y, m3)
self.assert_array_equal(self.average(z, None), 20./6.)
self.assert_array_equal(self.average(z, axis=0), [0., 1., 99., 99., 4.0, 7.5])
self.assert_array_equal(self.average(z, axis=1), [2.5, 5.0])
self.assert_array_equal(self.average(z, axis=0, weights=w2), [0., 1., 99., 99., 4.0, 10.0])
@np.errstate(all='ignore')
def test_A(self):
x = self.arange(24)
x[5:6] = self.masked
x = x.reshape(2, 3, 4)
if __name__ == '__main__':
setup_base = ("from __main__ import ModuleTester \n"
"import numpy\n"
"tester = ModuleTester(module)\n")
setup_cur = "import numpy.ma.core as module\n" + setup_base
(nrepeat, nloop) = (10, 10)
for i in range(1, 8):
func = 'tester.test_%i()' % i
cur = timeit.Timer(func, setup_cur).repeat(nrepeat, nloop*10)
cur = np.sort(cur)
print("#%i" % i + 50*'.')
print(eval("ModuleTester.test_%i.__doc__" % i))
print(f'core_current : {cur[0]:.3f} - {cur[1]:.3f}')
.\numpy\numpy\ma\__init__.py
"""
=============
Masked Arrays
=============
Arrays sometimes contain invalid or missing data. When doing operations
on such arrays, we wish to suppress invalid values, which is the purpose masked
arrays fulfill (an example of typical use is given below).
For example, examine the following array:
>>> x = np.array([2, 1, 3, np.nan, 5, 2, 3, np.nan])
When we try to calculate the mean of the data, the result is undetermined:
>>> np.mean(x)
nan
The mean is calculated using roughly ``np.sum(x)/len(x)``, but since
any number added to ``NaN`` [1]_ produces ``NaN``, this doesn't work. Enter
masked arrays:
>>> m = np.ma.masked_array(x, np.isnan(x))
>>> m
masked_array(data=[2.0, 1.0, 3.0, --, 5.0, 2.0, 3.0, --],
mask=[False, False, False, True, False, False, False, True],
fill_value=1e+20)
Here, we construct a masked array that suppresses all ``NaN`` values. We
may now proceed to calculate the mean of the other values:
>>> np.mean(m)
2.6666666666666665
.. [1] Not-a-Number, a floating point value that is the result of an
invalid operation.
.. moduleauthor:: Pierre Gerard-Marchant
.. moduleauthor:: Jarrod Millman
"""
from . import core
from .core import *
from . import extras
from .extras import *
__all__ = ['core', 'extras']
__all__ += core.__all__
__all__ += extras.__all__
from numpy._pytesttester import PytestTester
test = PytestTester(__name__)
del PytestTester
.\numpy\numpy\ma\__init__.pyi
from numpy._pytesttester import PytestTester
from numpy.ma import extras as extras
from numpy.ma.core import (
MAError as MAError,
MaskError as MaskError,
MaskType as MaskType,
MaskedArray as MaskedArray,
abs as abs,
absolute as absolute,
add as add,
all as all,
allclose as allclose,
allequal as allequal,
alltrue as alltrue,
amax as amax,
amin as amin,
angle as angle,
anom as anom,
anomalies as anomalies,
any as any,
append as append,
arange as arange,
arccos as arccos,
arccosh as arccosh,
arcsin as arcsin,
arcsinh as arcsinh,
arctan as arctan,
arctan2 as arctan2,
arctanh as arctanh,
argmax as argmax,
argmin as argmin,
argsort as argsort,
around as around,
array as array,
asanyarray as asanyarray,
asarray as asarray,
bitwise_and as bitwise_and,
bitwise_or as bitwise_or,
bitwise_xor as bitwise_xor,
bool as bool,
ceil as ceil,
choose as choose,
clip as clip,
common_fill_value as common_fill_value,
compress as compress,
compressed as compressed,
concatenate as concatenate,
conjugate as conjugate,
convolve as convolve,
copy as copy,
correlate as correlate,
cos as cos,
cosh as cosh,
count as count,
cumprod as cumprod,
cumsum as cumsum,
default_fill_value as default_fill_value,
diag as diag,
diagonal as diagonal,
diff as diff,
divide as divide,
empty as empty,
empty_like as empty_like,
equal as equal,
exp as exp,
expand_dims as expand_dims,
fabs as fabs,
filled as filled,
fix_invalid as fix_invalid,
flatten_mask as flatten_mask,
flatten_structured_array as flatten_structured_array,
floor as floor,
floor_divide as floor_divide,
fmod as fmod,
frombuffer as frombuffer,
fromflex as fromflex,
fromfunction as fromfunction,
getdata as getdata,
getmask as getmask,
getmaskarray as getmaskarray,
greater as greater,
greater_equal as greater_equal,
harden_mask as harden_mask,
hypot as hypot,
identity as identity,
ids as ids,
indices as indices,
inner as inner,
innerproduct as innerproduct,
isMA as isMA,
isMaskedArray as isMaskedArray,
is_mask as is_mask,
is_masked as is_masked,
isarray as isarray,
left_shift as left_shift,
less as
masked_invalid as masked_invalid,
masked_less as masked_less,
masked_less_equal as masked_less_equal,
masked_not_equal as masked_not_equal,
masked_object as masked_object,
masked_outside as masked_outside,
masked_print_option as masked_print_option,
masked_singleton as masked_singleton,
masked_values as masked_values,
masked_where as masked_where,
max as max,
maximum as maximum,
maximum_fill_value as maximum_fill_value,
mean as mean,
min as min,
minimum as minimum,
minimum_fill_value as minimum_fill_value,
mod as mod,
multiply as multiply,
mvoid as mvoid,
ndim as ndim,
negative as negative,
nomask as nomask,
nonzero as nonzero,
not_equal as not_equal,
ones as ones,
outer as outer,
outerproduct as outerproduct,
power as power,
prod as prod,
product as product,
ptp as ptp,
put as put,
putmask as putmask,
ravel as ravel,
remainder as remainder,
repeat as repeat,
reshape as reshape,
resize as resize,
right_shift as right_shift,
round as round,
set_fill_value as set_fill_value,
shape as shape,
sin as sin,
sinh as sinh,
size as size,
soften_mask as soften_mask,
sometrue as sometrue,
sort as sort,
sqrt as sqrt,
squeeze as squeeze,
std as std,
subtract as subtract,
sum as sum,
swapaxes as swapaxes,
take as take,
tan as tan,
tanh as tanh,
trace as trace,
transpose as transpose,
true_divide as true_divide,
var as var,
where as where,
zeros as zeros,
from numpy.ma.extras import (
apply_along_axis as apply_along_axis,
apply_over_axes as apply_over_axes,
atleast_1d as atleast_1d,
atleast_2d as atleast_2d,
atleast_3d as atleast_3d,
average as average,
clump_masked as clump_masked,
clump_unmasked as clump_unmasked,
column_stack as column_stack,
compress_cols as compress_cols,
compress_nd as compress_nd,
compress_rowcols as compress_rowcols,
compress_rows as compress_rows,
count_masked as count_masked,
corrcoef as corrcoef,
cov as cov,
diagflat as diagflat,
dot as dot,
dstack as dstack,
ediff1d as ediff1d,
flatnotmasked_contiguous as flatnotmasked_contiguous,
flatnotmasked_edges as flatnotmasked_edges,
hsplit as hsplit,
hstack as hstack,
isin as isin,
in1d as in1d,
intersect1d as intersect1d,
mask_cols as mask_cols,
mask_rowcols as mask_rowcols,
mask_rows as mask_rows,
masked_all as masked_all,
masked_all_like as masked_all_like,
median as median,
mr_ as mr_,
ndenumerate as ndenumerate,
notmasked_contiguous as notmasked_contiguous,
notmasked_edges as notmasked_edges,
polyfit as polyfit,
row_stack as row_stack,
setdiff1d as setdiff1d,
setxor1d as setxor1d,
stack as stack,
unique as unique,
union1d as union1d,
vander as vander,
vstack as vstack
)
__all__: list[str]
test: PytestTester
.\numpy\numpy\matlib.py
import warnings
warnings.warn("Importing from numpy.matlib is deprecated since 1.19.0. "
"The matrix subclass is not the recommended way to represent "
"matrices or deal with linear algebra (see "
"https://docs.scipy.org/doc/numpy/user/numpy-for-matlab-users.html). "
"Please adjust your code to use regular ndarray. ",
PendingDeprecationWarning, stacklevel=2)
import numpy as np
from numpy.matrixlib.defmatrix import matrix, asmatrix
from numpy import *
__version__ = np.__version__
__all__ = np.__all__[:]
__all__ += ['rand', 'randn', 'repmat']
def empty(shape, dtype=None, order='C'):
"""Return a new matrix of given shape and type, without initializing entries.
Parameters
----------
shape : int or tuple of int
Shape of the empty matrix.
dtype : data-type, optional
Desired output data-type.
order : {'C', 'F'}, optional
Whether to store multi-dimensional data in row-major
(C-style) or column-major (Fortran-style) order in
memory.
See Also
--------
numpy.empty : Equivalent array function.
matlib.zeros : Return a matrix of zeros.
matlib.ones : Return a matrix of ones.
Notes
-----
Unlike other matrix creation functions (e.g. `matlib.zeros`,
`matlib.ones`), `matlib.empty` does not initialize the values of the
matrix, and may therefore be marginally faster. However, the values
stored in the newly allocated matrix are arbitrary. For reproducible
behavior, be sure to set each element of the matrix before reading.
Examples
--------
>>> import numpy.matlib
>>> np.matlib.empty((2, 2)) # filled with random data
matrix([[ 6.76425276e-320, 9.79033856e-307], # random
[ 7.39337286e-309, 3.22135945e-309]])
>>> np.matlib.empty((2, 2), dtype=int)
matrix([[ 6600475, 0], # random
[ 6586976, 22740995]])
"""
return ndarray.__new__(matrix, shape, dtype, order=order)
def ones(shape, dtype=None, order='C'):
"""
Matrix of ones.
Return a matrix of given shape and type, filled with ones.
Parameters
----------
shape : {sequence of ints, int}
Shape of the matrix
dtype : data-type, optional
The desired data-type for the matrix, default is np.float64.
order : {'C', 'F'}, optional
Whether to store matrix in C- or Fortran-contiguous order,
default is 'C'.
Returns
-------
out : matrix
Matrix of ones of given shape, dtype, and order.
See Also
--------
ones : Array of ones.
"""
return np.ones(shape, dtype=dtype, order=order)
def zeros(shape, dtype=None, order='C'):
"""
a = ndarray.__new__(matrix, shape, dtype, order=order)
使用给定的 shape、dtype 和 order 创建一个新的 matrix 对象 a。
"""
a = ndarray.__new__(matrix, shape, dtype, order=order)
a.fill(1)
return a
def zeros(shape, dtype=None, order='C'):
a = ndarray.__new__(matrix, shape, dtype, order=order)
a.fill(0)
return a
def identity(n,dtype=None):
a = array([1]+n*[0], dtype=dtype)
b = empty((n, n), dtype=dtype)
b.flat = a
return b
def eye(n,M=None, k=0, dtype=float, order='C'):
I = zeros((n, M if M else n), dtype=dtype, order=order)
if k >= 0:
I.flat[k::I.shape[1]+1] = 1
else:
I.flat[-k*I.shape[1]::I.shape[1]-1] = 1
return I
>>> np.matlib.eye(3, k=1, dtype=float)
matrix([[0., 1., 0.],
[0., 0., 1.],
[0., 0., 0.]])
"""
返回一个 n x M 的单位矩阵的 numpy 矩阵表示，其中对角线向上偏移 k 个位置。
def repmat(a, m, n):
"""
Repeat a 0-D to 2-D array or matrix MxN times.
Parameters
----------
a : array_like
The array or matrix to be repeated.
m : int
Number of times the input array `a` is repeated along the first dimension.
n : int
Number of times the input array `a` is repeated along the second dimension, if `a` is 2-D.
Returns
-------
out : ndarray
The tiled output array.
See Also
--------
tile, numpy.matlib.repmat
Notes
-----
If `a` is 0-D or 1-D, the result will be a 1-D array repeating `a` `m` times.
If `a` is 2-D, the result will be a 2-D array repeating `a` `m` times along the first axis and `n` times along the second axis.
Examples
--------
>>> np.matlib.repmat(np.array([[1, 2], [3, 4]]), 2, 3)
matrix([[1, 2, 1, 2, 1, 2],
[3, 4, 3, 4, 3, 4],
[1, 2, 1, 2, 1, 2],
[3, 4, 3, 4, 3, 4]])
>>> np.matlib.repmat(np.array([1, 2]), 3, 2)
matrix([[1, 2],
[1, 2],
[1, 2]])
"""
# 使用 numpy 的 tile 函数实现重复数组 a 的操作
return np.tile(np.asarray(a), (m, n))
# 将输入参数 `a` 转换为多维数组
a = asanyarray(a)
# 获取数组 `a` 的维度
ndim = a.ndim
# 根据数组的维度确定原始行数和列数
if ndim == 0:
origrows, origcols = (1, 1) # 对于零维数组，原始行数和列数均为1
elif ndim == 1:
origrows, origcols = (1, a.shape[0]) # 对于一维数组，原始行数为1，列数为数组长度
else:
origrows, origcols = a.shape # 对于高维数组，原始行数和列数即为数组的形状
# 计算重复后的总行数和总列数
rows = origrows * m
cols = origcols * n
# 将数组 `a` 展平为一维，并按指定的行数 `m` 重复堆叠，然后按原始列数重新整形
# 再按指定的列数 `n` 重复堆叠，形成最终的重复数组
c = a.reshape(1, a.size).repeat(m, 0).reshape(rows, origcols).repeat(n, 0)
# 返回重复后的数组，重新整形为指定的总行数和总列数
return c.reshape(rows, cols)
.\numpy\numpy\matrixlib\defmatrix.py
__all__ = ['matrix', 'bmat', 'asmatrix']
import sys
import warnings
import ast
from .._utils import set_module
import numpy._core.numeric as N
from numpy._core.numeric import concatenate, isscalar
from numpy.linalg import matrix_power
def _convert_from_string(data):
for char in '[]':
data = data.replace(char, '')
rows = data.split(';')
newdata = []
count = 0
for row in rows:
trow = row.split(',')
newrow = []
for col in trow:
temp = col.split()
newrow.extend(map(ast.literal_eval, temp))
if count == 0:
Ncols = len(newrow)
elif len(newrow) != Ncols:
raise ValueError("Rows not the same size.")
count += 1
newdata.append(newrow)
return newdata
@set_module('numpy')
def asmatrix(data, dtype=None):
"""
Interpret the input as a matrix.
Unlike `matrix`, `asmatrix` does not make a copy if the input is already
a matrix or an ndarray. Equivalent to ``matrix(data, copy=False)``.
Parameters
----------
data : array_like
Input data.
dtype : data-type
Data-type of the output matrix.
Returns
-------
mat : matrix
`data` interpreted as a matrix.
Examples
--------
>>> x = np.array([[1, 2], [3, 4]])
>>> m = np.asmatrix(x)
>>> x[0,0] = 5
>>> m
matrix([[5, 2],
[3, 4]])
"""
return matrix(data, dtype=dtype, copy=False)
@set_module('numpy')
class matrix(N.ndarray):
"""
matrix(data, dtype=None, copy=True)
Returns a matrix from an array-like object, or from a string of data.
A matrix is a specialized 2-D array that retains its 2-D nature
through operations. It has certain special operators, such as ``*``
(matrix multiplication) and ``**`` (matrix power).
.. note:: It is no longer recommended to use this class, even for linear
algebra. Instead use regular arrays. The class may be removed
in the future.
Parameters
----------
data : array_like or string
If `data` is a string, it is interpreted as a matrix with commas
or spaces separating columns, and semicolons separating rows.
dtype : data-type
Data-type of the output matrix.
copy : bool
If `data` is already an `ndarray`, then this flag determines
whether the data is copied (the default), or whether a view is
constructed.
See Also
--------
array
Examples
--------
>>> a = np.matrix('1 2; 3 4')
>>> a
matrix([[1, 2],
[3, 4]])
>>> np.matrix([[1, 2], [3, 4]])
matrix([[1, 2],
[3, 4]])
"""
__array_priority__ = 10.0
def __new__(subtype, data, dtype=None, copy=True):
warnings.warn('the matrix subclass is not the recommended way to '
'represent matrices or deal with linear algebra (see '
'https://docs.scipy.org/doc/numpy/user/'
'numpy-for-matlab-users.html). '
'Please adjust your code to use regular ndarray.',
PendingDeprecationWarning, stacklevel=2)
if isinstance(data, matrix):
dtype2 = data.dtype
if (dtype is None):
dtype = dtype2
if (dtype2 == dtype) and (not copy):
return data
return data.astype(dtype)
if isinstance(data, N.ndarray):
if dtype is None:
intype = data.dtype
else:
intype = N.dtype(dtype)
new = data.view(subtype)
if intype != data.dtype:
return new.astype(intype)
if copy: return new.copy()
else: return new
if isinstance(data, str):
data = _convert_from_string(data)
copy = None if not copy else True
arr = N.array(data, dtype=dtype, copy=copy)
ndim = arr.ndim
shape = arr.shape
if (ndim > 2):
raise ValueError("matrix must be 2-dimensional")
elif ndim == 0:
shape = (1, 1)
elif ndim == 1:
shape = (1, shape[0])
order = 'C'
if (ndim == 2) and arr.flags.fortran:
order = 'F'
if not (order or arr.flags.contiguous):
arr = arr.copy()
ret = N.ndarray.__new__(subtype, shape, arr.dtype,
buffer=arr,
order=order)
return ret
def __array_finalize__(self, obj):
self._getitem = False
if (isinstance(obj, matrix) and obj._getitem): return
ndim = self.ndim
if (ndim == 2):
return
if (ndim > 2):
newshape = tuple([x for x in self.shape if x > 1])
ndim = len(newshape)
if ndim == 2:
self.shape = newshape
return
elif (ndim > 2):
raise ValueError("shape too large to be a matrix.")
else:
newshape = self.shape
if ndim == 0:
self.shape = (1, 1)
elif ndim == 1:
self.shape = (1, newshape[0])
return
def __getitem__(self, index):
self._getitem = True
try:
out = N.ndarray.__getitem__(self, index)
finally:
self._getitem = False
if not isinstance(out, N.ndarray):
return out
if out.ndim == 0:
return out[()]
if out.ndim == 1:
sh = out.shape[0]
try:
n = len(index)
except Exception:
n = 0
if n > 1 and isscalar(index[1]):
out.shape = (sh, 1)
else:
out.shape = (1, sh)
return out
def __mul__(self, other):
if isinstance(other, (N.ndarray, list, tuple)) :
return N.dot(self, asmatrix(other))
if isscalar(other) or not hasattr(other, '__rmul__') :
return N.dot(self, other)
return NotImplemented
def __rmul__(self, other):
return N.dot(other, self)
def __imul__(self, other):
self[:] = self * other
return self
def __pow__(self, other):
return matrix_power(self, other)
def __ipow__(self, other):
self[:] = self ** other
return self
def __rpow__(self, other):
return NotImplemented
def _align(self, axis):
"""A convenience function for operations that need to preserve axis
orientation.
"""
if axis is None:
return self[0, 0]
elif axis==0:
return self
elif axis==1:
return self.transpose()
else:
raise ValueError("unsupported axis")
def _collapse(self, axis):
"""A convenience function for operations that want to collapse
to a scalar like _align, but are using keepdims=True
"""
if axis is None:
return self[0, 0]
else:
return self
def tolist(self):
"""
Return the matrix as a (possibly nested) list.
See `ndarray.tolist` for full documentation.
See Also
--------
ndarray.tolist
Examples
--------
>>> x = np.matrix(np.arange(12).reshape((3,4))); x
matrix([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
>>> x.tolist()
[[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11]]
"""
return self.__array__().tolist()
def sum(self, axis=None, dtype=None, out=None):
"""
Returns the sum of the matrix elements, along the given axis.
Refer to `numpy.sum` for full documentation.
See Also
--------
numpy.sum
Notes
-----
This is the same as `ndarray.sum`, except that where an `ndarray` would
be returned, a `matrix` object is returned instead.
Examples
--------
>>> x = np.matrix([[1, 2], [4, 3]])
>>> x.sum()
10
>>> x.sum(axis=1)
matrix([[3],
[7]])
>>> x.sum(axis=1, dtype='float')
matrix([[3.],
[7.]])
>>> out = np.zeros((2, 1), dtype='float')
>>> x.sum(axis=1, dtype='float', out=np.asmatrix(out))
matrix([[3.],
[7.]])
"""
return N.ndarray.sum(self, axis, dtype, out, keepdims=True)._collapse(axis)
def squeeze(self, axis=None):
"""
Return a possibly reshaped matrix.
Refer to `numpy.squeeze` for more documentation.
Parameters
----------
axis : None or int or tuple of ints, optional
Selects a subset of the axes of length one in the shape.
If an axis is selected with shape entry greater than one,
an error is raised.
Returns
-------
squeezed : matrix
The matrix, but as a (1, N) matrix if it had shape (N, 1).
See Also
--------
numpy.squeeze : related function
Notes
-----
If `m` has a single column then that column is returned
as the single row of a matrix. Otherwise `m` is returned.
The returned matrix is always either `m` itself or a view into `m`.
Supplying an axis keyword argument will not affect the returned matrix
but it may cause an error to be raised.
Examples
--------
>>> c = np.matrix([[1], [2]])
>>> c
matrix([[1],
[2]])
>>> c.squeeze()
matrix([[1, 2]])
>>> r = c.T
>>> r
matrix([[1, 2]])
>>> r.squeeze()
matrix([[1, 2]])
>>> m = np.matrix([[1, 2], [3, 4]])
>>> m.squeeze()
matrix([[1, 2],
[3, 4]])
"""
return N.ndarray.squeeze(self, axis=axis)
def flatten(self, order='C'):
"""
Return a flattened copy of the matrix.
All `N` elements of the matrix are placed into a single row.
Parameters
----------
order : {'C', 'F', 'A', 'K'}, optional
'C' means to flatten in row-major (C-style) order. 'F' means to
flatten in column-major (Fortran-style) order. 'A' means to
flatten in column-major order if `m` is Fortran *contiguous* in
memory, row-major order otherwise. 'K' means to flatten `m` in
the order the elements occur in memory. The default is 'C'.
Returns
-------
y : matrix
A copy of the matrix, flattened to a `(1, N)` matrix where `N`
is the number of elements in the original matrix.
See Also
--------
ravel : Return a flattened array.
flat : A 1-D flat iterator over the matrix.
Examples
--------
>>> m = np.matrix([[1,2], [3,4]])
>>> m.flatten()
matrix([[1, 2, 3, 4]])
>>> m.flatten('F')
matrix([[1, 3, 2, 4]])
"""
return N.ndarray.flatten(self, order=order)
def mean(self, axis=None, dtype=None, out=None):
"""
Returns the average of the matrix elements along the given axis.
Refer to `numpy.mean` for full documentation.
See Also
--------
numpy.mean
Notes
-----
Same as `ndarray.mean` except that, where that returns an `ndarray`,
this returns a `matrix` object.
Examples
--------
>>> x = np.matrix(np.arange(12).reshape((3, 4)))
>>> x
matrix([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
>>> x.mean()
5.5
>>> x.mean(0)
matrix([[4., 5., 6., 7.]])
>>> x.mean(1)
matrix([[ 1.5],
[ 5.5],
[ 9.5]])
"""
return N.ndarray.mean(self, axis, dtype, out, keepdims=True)._collapse(axis)
def std(self, axis=None, dtype=None, out=None, ddof=0):
"""
返回沿着给定轴的数组元素的标准差。
参考 `numpy.std` 获取完整文档。
See Also
--------
numpy.std
Notes
-----
这与 `ndarray.std` 相同，不同之处在于 `ndarray` 返回时，返回的是一个 `matrix` 对象。
Examples
--------
>>> x = np.matrix(np.arange(12).reshape((3, 4)))
>>> x
matrix([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
>>> x.std()
3.4520525295346629 # 可能有所不同
>>> x.std(0)
matrix([[ 3.26598632, 3.26598632, 3.26598632, 3.26598632]]) # 可能有所不同
>>> x.std(1)
matrix([[ 1.11803399],
[ 1.11803399],
[ 1.11803399]])
"""
return N.ndarray.std(self, axis, dtype, out, ddof, keepdims=True)._collapse(axis)
def var(self, axis=None, dtype=None, out=None, ddof=0):
"""
返回沿着给定轴的矩阵元素的方差。
参考 `numpy.var` 获取完整文档。
See Also
--------
numpy.var
Notes
-----
这与 `ndarray.var` 相同，不同之处在于 `ndarray` 返回时，返回的是一个 `matrix` 对象。
Examples
--------
>>> x = np.matrix(np.arange(12).reshape((3, 4)))
>>> x
matrix([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
>>> x.var()
11.916666666666666
>>> x.var(0)
matrix([[ 10.66666667, 10.66666667, 10.66666667, 10.66666667]]) # 可能有所不同
>>> x.var(1)
matrix([[1.25],
[1.25],
[1.25]])
"""
return N.ndarray.var(self, axis, dtype, out, ddof, keepdims=True)._collapse(axis)
def prod(self, axis=None, dtype=None, out=None):
"""
返回沿着给定轴的数组元素的乘积。
参考 `prod` 获取完整文档。
See Also
--------
prod, ndarray.prod
Notes
-----
与 `ndarray.prod` 相同，不同之处在于 `ndarray` 返回时，返回的是一个 `matrix` 对象。
Examples
--------
>>> x = np.matrix(np.arange(12).reshape((3,4))); x
matrix([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
>>> x.prod()
0
>>> x.prod(0)
matrix([[ 0, 45, 120, 231]])
>>> x.prod(1)
matrix([[ 0],
[ 840],
[7920]])
"""
return N.ndarray.prod(self, axis, dtype, out, keepdims=True)._collapse(axis)
def any(self, axis=None, out=None):
"""
Test whether any array element along a given axis evaluates to True.
Refer to `numpy.any` for full documentation.
Parameters
----------
axis : int, optional
Axis along which logical OR is performed
out : ndarray, optional
Output to existing array instead of creating new one, must have
same shape as expected output
Returns
-------
any : bool, ndarray
Returns a single bool if `axis` is ``None``; otherwise,
returns `ndarray`
"""
return N.ndarray.any(self, axis, out, keepdims=True)._collapse(axis)
def all(self, axis=None, out=None):
"""
Test whether all matrix elements along a given axis evaluate to True.
Parameters
----------
See `numpy.all` for complete descriptions
See Also
--------
numpy.all
Notes
-----
This is the same as `ndarray.all`, but it returns a `matrix` object.
Examples
--------
>>> x = np.matrix(np.arange(12).reshape((3,4))); x
matrix([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
>>> y = x[0]; y
matrix([[0, 1, 2, 3]])
>>> (x == y)
matrix([[ True, True, True, True],
[False, False, False, False],
[False, False, False, False]])
>>> (x == y).all()
False
>>> (x == y).all(0)
matrix([[False, False, False, False]])
>>> (x == y).all(1)
matrix([[ True],
[False],
[False]])
"""
return N.ndarray.all(self, axis, out, keepdims=True)._collapse(axis)
def max(self, axis=None, out=None):
"""
Return the maximum value along an axis.
Parameters
----------
See `amax` for complete descriptions
See Also
--------
amax, ndarray.max
Notes
-----
This is the same as `ndarray.max`, but returns a `matrix` object
where `ndarray.max` would return an ndarray.
Examples
--------
>>> x = np.matrix(np.arange(12).reshape((3,4))); x
matrix([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
>>> x.max()
11
>>> x.max(0)
matrix([[ 8, 9, 10, 11]])
>>> x.max(1)
matrix([[ 3],
[ 7],
[11]])
"""
return N.ndarray.max(self, axis, out, keepdims=True)._collapse(axis)
def argmax(self, axis=None, out=None):
"""
返回沿着指定轴的最大值的索引。
返回沿着指定轴的最大值的第一个出现的索引。如果 axis 是 None，则索引是针对展平的矩阵。
Parameters
----------
axis : int, optional
沿着哪个轴查找最大值的索引。默认为 None。
out : ndarray, optional
结果的放置位置。默认为 None。
See Also
--------
numpy.argmax
Notes
-----
这与 `ndarray.argmax` 相同，但返回一个 `matrix` 对象，而 `ndarray.argmax` 返回一个 `ndarray`。
Examples
--------
>>> x = np.matrix(np.arange(12).reshape((3,4))); x
matrix([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
>>> x.argmax()
11
>>> x.argmax(0)
matrix([[2, 2, 2, 2]])
>>> x.argmax(1)
matrix([[3],
[3],
[3]])
"""
return N.ndarray.argmax(self, axis, out)._align(axis)
def min(self, axis=None, out=None):
"""
返回沿着指定轴的最小值。
Parameters
----------
axis : int, optional
沿着哪个轴查找最小值。默认为 None。
out : ndarray, optional
结果的放置位置。默认为 None。
See Also
--------
numpy.amin, ndarray.min
Notes
-----
这与 `ndarray.min` 相同，但返回一个 `matrix` 对象，而 `ndarray.min` 返回一个 `ndarray`。
Examples
--------
>>> x = -np.matrix(np.arange(12).reshape((3,4))); x
matrix([[ 0, -1, -2, -3],
[ -4, -5, -6, -7],
[ -8, -9, -10, -11]])
>>> x.min()
-11
>>> x.min(0)
matrix([[ -8, -9, -10, -11]])
>>> x.min(1)
matrix([[ -3],
[ -7],
[-11]])
"""
return N.ndarray.min(self, axis, out, keepdims=True)._collapse(axis)
def argmin(self, axis=None, out=None):
"""
返回沿着指定轴的最小值的索引。
返回沿着指定轴的最小值的第一个出现的索引。如果 axis 是 None，则索引是针对展平的矩阵。
Parameters
----------
axis : int, optional
沿着哪个轴查找最小值的索引。默认为 None。
out : ndarray, optional
结果的放置位置。默认为 None。
See Also
--------
numpy.argmin
Notes
-----
这与 `ndarray.argmin` 相同，但返回一个 `matrix` 对象，而 `ndarray.argmin` 返回一个 `ndarray`。
Examples
--------
>>> x = -np.matrix(np.arange(12).reshape((3,4))); x
matrix([[ 0, -1, -2, -3],
[ -4, -5, -6, -7],
[ -8, -9, -10, -11]])
>>> x.argmin()
11
>>> x.argmin(0)
matrix([[2, 2, 2, 2]])
>>> x.argmin(1)
matrix([[3],
[3],
[3]])
"""
return N.ndarray.argmin(self, axis, out)._align(axis)
def ptp(self, axis=None, out=None):
"""
Peak-to-peak (maximum - minimum) value along the given axis.
Refer to `numpy.ptp` for full documentation.
See Also
--------
numpy.ptp
Notes
-----
Same as `ndarray.ptp`, except, where that would return an `ndarray` object,
this returns a `matrix` object.
Examples
--------
>>> x = np.matrix(np.arange(12).reshape((3,4))); x
matrix([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
>>> x.ptp()
11
>>> x.ptp(0)
matrix([[8, 8, 8, 8]])
>>> x.ptp(1)
matrix([[3],
[3],
[3]])
"""
return N.ptp(self, axis, out)._align(axis)
@property
def I(self):
"""
Returns the (multiplicative) inverse of invertible `self`.
Parameters
----------
None
Returns
-------
ret : matrix object
If `self` is non-singular, `ret` is such that ``ret * self`` ==
``self * ret`` == ``np.matrix(np.eye(self[0,:].size))`` all return
``True``.
Raises
------
numpy.linalg.LinAlgError: Singular matrix
If `self` is singular.
See Also
--------
linalg.inv
Examples
--------
>>> m = np.matrix('[1, 2; 3, 4]'); m
matrix([[1, 2],
[3, 4]])
>>> m.getI()
matrix([[-2. , 1. ],
[ 1.5, -0.5]])
>>> m.getI() * m
matrix([[ 1., 0.], # may vary
[ 0., 1.]])
"""
M, N = self.shape
if M == N:
from numpy.linalg import inv as func
else:
from numpy.linalg import pinv as func
return asmatrix(func(self))
@property
def A(self):
"""
Return `self` as an `ndarray` object.
Equivalent to ``np.asarray(self)``.
Parameters
----------
None
Returns
-------
ret : ndarray
`self` as an `ndarray`
Examples
--------
>>> x = np.matrix(np.arange(12).reshape((3,4))); x
matrix([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
>>> x.getA()
array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
"""
return self.__array__()
@property
def H(self):
"""
Return the conjugate transpose of `self`.
Equivalent to ``np.asarray(self).conj().T``.
Parameters
----------
None
Returns
-------
ret : matrix
Conjugate transpose of `self`
Examples
--------
>>> x = np.matrix([[1+1j, 2+2j], [3+3j, 4+4j]])
>>> x.H
matrix([[1.-1.j, 3.-3.j],
[2.-2.j, 4.-4.j]])
"""
return self.__array__().conj().T
def A1(self):
"""
Return `self` as a flattened `ndarray`.
Equivalent to ``np.asarray(x).ravel()``
Parameters
----------
None
Returns
-------
ret : ndarray
`self`, 1-D, as an `ndarray`
Examples
--------
>>> x = np.matrix(np.arange(12).reshape((3,4))); x
matrix([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
>>> x.getA1()
array([ 0, 1, 2, ..., 9, 10, 11])
"""
return self.__array__().ravel()
def ravel(self, order='C'):
"""
Return a flattened matrix.
Refer to `numpy.ravel` for more documentation.
Parameters
----------
order : {'C', 'F', 'A', 'K'}, optional
The elements of `m` are read using this index order. 'C' means to
index the elements in C-like order, with the last axis index
changing fastest, back to the first axis index changing slowest.
'F' means to index the elements in Fortran-like index order, with
the first index changing fastest, and the last index changing
slowest. Note that the 'C' and 'F' options take no account of the
memory layout of the underlying array, and only refer to the order
of axis indexing. 'A' means to read the elements in Fortran-like
index order if `m` is Fortran *contiguous* in memory, C-like order
otherwise. 'K' means to read the elements in the order they occur
in memory, except for reversing the data when strides are negative.
By default, 'C' index order is used.
Returns
-------
ret : matrix
Return the matrix flattened to shape `(1, N)` where `N`
is the number of elements in the original matrix.
A copy is made only if necessary.
See Also
--------
matrix.flatten : returns a similar output matrix but always a copy
matrix.flat : a flat iterator on the array.
numpy.ravel : related function which returns an ndarray
"""
return N.ndarray.ravel(self, order=order)
@property
def T(self):
"""
Returns the transpose of the matrix.
Does *not* conjugate! For the complex conjugate transpose, use ``.H``.
Parameters
----------
None
Returns
-------
ret : matrix object
The (non-conjugated) transpose of the matrix.
See Also
--------
transpose, getH
Examples
--------
>>> m = np.matrix('[1, 2; 3, 4]')
>>> m
matrix([[1, 2],
[3, 4]])
>>> m.getT()
matrix([[1, 3],
[2, 4]])
"""
return self.transpose()
@property
getT = T.fget
getA = A.fget
getA1 = A1.fget
getH = H.fget
getI = I.fget
def bmat(obj, ldict=None, gdict=None):
"""
Build a matrix object from a string, nested sequence, or array.
Parameters
----------
obj : str or array_like
Input data. If a string, variables in the current scope may be
referenced by name.
ldict : dict, optional
A dictionary that replaces local operands in current frame.
Ignored if `obj` is not a string or `gdict` is None.
gdict : dict, optional
A dictionary that replaces global operands in current frame.
Ignored if `obj` is not a string.
Returns
-------
out : matrix
Returns a matrix object, which is a specialized 2-D array.
See Also
--------
block :
A generalization of this function for N-d arrays, that returns normal
ndarrays.
Examples
--------
>>> A = np.asmatrix('1 1; 1 1')
>>> B = np.asmatrix('2 2; 2 2')
>>> C = np.asmatrix('3 4; 5 6')
>>> D = np.asmatrix('7 8; 9 0')
All the following expressions construct the same block matrix:
>>> np.bmat([[A, B], [C, D]])
matrix([[1, 1, 2, 2],
[1, 1, 2, 2],
[3, 4, 7, 8],
[5, 6, 9, 0]])
>>> np.bmat(np.r_[np.c_[A, B], np.c_[C, D]])
matrix([[1, 1, 2, 2],
[1, 1, 2, 2],
[3, 4, 7, 8],
[5, 6, 9, 0]])
>>> np.bmat('A,B; C,D')
matrix([[1, 1, 2, 2],
[1, 1, 2, 2],
[3, 4, 7, 8],
[5, 6, 9, 0]])
"""
if isinstance(obj, str):
if gdict is None:
frame = sys._getframe().f_back
glob_dict = frame.f_globals
loc_dict = frame.f_locals
else:
glob_dict = gdict
loc_dict = ldict
return matrix(_from_string(obj, glob_dict, loc_dict))
if isinstance(obj, (tuple, list)):
arr_rows = []
for row in obj:
if isinstance(row, N.ndarray):
return matrix(concatenate(obj, axis=-1))
else:
arr_rows.append(concatenate(row, axis=-1))
return matrix(concatenate(arr_rows, axis=0))
if isinstance(obj, N.ndarray):
return matrix(obj)
.\numpy\numpy\matrixlib\defmatrix.pyi
from collections.abc import Sequence, Mapping
from typing import Any
from numpy import matrix as matrix
from numpy._typing import ArrayLike, DTypeLike, NDArray
__all__: list[str]
def bmat(
obj: str | Sequence[ArrayLike] | NDArray[Any],
ldict: None | Mapping[str, Any] = ...,
gdict: None | Mapping[str, Any] = ...,
) -> matrix[Any, Any]:
def asmatrix(data: ArrayLike, dtype: DTypeLike = ...) -> matrix[Any, Any]:
mat = asmatrix
.\numpy\numpy\matrixlib\tests\test_defmatrix.py
import collections.abc
import numpy as np
from numpy import matrix, asmatrix, bmat
from numpy.testing import (
assert_, assert_equal, assert_almost_equal, assert_array_equal,
assert_array_almost_equal, assert_raises
)
from numpy.linalg import matrix_power
class TestCtor:
def test_basic(self):
A = np.array([[1, 2], [3, 4]])
mA = matrix(A)
assert_(np.all(mA.A == A))
B = bmat("A,A;A,A")
C = bmat([[A, A], [A, A]])
D = np.array([[1, 2, 1, 2],
[3, 4, 3, 4],
[1, 2, 1, 2],
[3, 4, 3, 4]])
assert_(np.all(B.A == D))
assert_(np.all(C.A == D))
E = np.array([[5, 6], [7, 8]])
AEresult = matrix([[1, 2, 5, 6], [3, 4, 7, 8]])
assert_(np.all(bmat([A, E]) == AEresult))
vec = np.arange(5)
mvec = matrix(vec)
assert_(mvec.shape == (1, 5))
def test_exceptions(self):
assert_raises(ValueError, matrix, "invalid")
def test_bmat_nondefault_str(self):
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])
Aresult = np.array([[1, 2, 1, 2],
[3, 4, 3, 4],
[1, 2, 1, 2],
[3, 4, 3, 4]])
mixresult = np.array([[1, 2, 5, 6],
[3, 4, 7, 8],
[5, 6, 1, 2],
[7, 8, 3, 4]])
assert_(np.all(bmat("A,A;A,A") == Aresult))
assert_(np.all(bmat("A,A;A,A", ldict={'A':B}) == Aresult))
assert_raises(TypeError, bmat, "A,A;A,A", gdict={'A':B})
assert_(
np.all(bmat("A,A;A,A", ldict={'A':A}, gdict={'A':B}) == Aresult))
b2 = bmat("A,B;C,D", ldict={'A':A,'B':B}, gdict={'C':B,'D':A})
assert_(np.all(b2 == mixresult))
class TestProperties:
def test_sum(self):
"""Test whether matrix.sum(axis=1) preserves orientation.
Fails in NumPy <= 0.9.6.2127.
"""
M = matrix([[1, 2, 0, 0],
[3, 4, 0, 0],
[1, 2, 1, 2],
[3, 4, 3, 4]])
sum0 = matrix([8, 12, 4, 6])
sum1 = matrix([3, 7, 6, 14]).T
sumall = 30
assert_array_equal(sum0, M.sum(axis=0))
assert_array_equal(sum1, M.sum(axis=1))
assert_equal(sumall, M.sum())
assert_array_equal(sum0, np.sum(M, axis=0))
assert_array_equal(sum1, np.sum(M, axis=1))
assert_equal(sumall, np.sum(M))
def test_prod(self):
x = matrix([[1, 2, 3], [4, 5, 6]])
assert_equal(x.prod(), 720)
assert_equal(x.prod(0), matrix([[4, 10, 18]]))
assert_equal(x.prod(1), matrix([[6], [120]]))
assert_equal(np.prod(x), 720)
assert_equal(np.prod(x, axis=0), matrix([[4, 10, 18]]))
assert_equal(np.prod(x, axis=1), matrix([[6], [120]]))
y = matrix([0, 1, 3])
assert_(y.prod() == 0)
def test_max(self):
x = matrix([[1, 2, 3], [4, 5, 6]])
assert_equal(x.max(), 6)
assert_equal(x.max(0), matrix([[4, 5, 6]]))
assert_equal(x.max(1), matrix([[3], [6]]))
assert_equal(np.max(x), 6)
assert_equal(np.max(x, axis=0), matrix([[4, 5, 6]]))
assert_equal(np.max(x, axis=1), matrix([[3], [6]]))
def test_min(self):
x = matrix([[1, 2, 3], [4, 5, 6]])
assert_equal(x.min(), 1)
assert_equal(x.min(0), matrix([[1, 2, 3]]))
assert_equal(x.min(1), matrix([[1], [4]]))
assert_equal(np.min(x), 1)
assert_equal(np.min(x, axis=0), matrix([[1, 2, 3]]))
assert_equal(np.min(x, axis=1), matrix([[1], [4]]))
def test_ptp(self):
x = np.arange(4).reshape((2, 2))
mx = x.view(np.matrix)
assert_(mx.ptp() == 3)
assert_(np.all(mx.ptp(0) == np.array([2, 2])))
assert_(np.all(mx.ptp(1) == np.array([1, 1])))
def test_var(self):
x = np.arange(9).reshape((3, 3))
mx = x.view(np.matrix)
assert_equal(x.var(ddof=0), mx.var(ddof=0))
assert_equal(x.var(ddof=1), mx.var(ddof=1))
def test_basic(self):
import numpy.linalg as linalg
A = np.array([[1., 2.],
[3., 4.]])
mA = matrix(A)
assert_(np.allclose(linalg.inv(A), mA.I))
assert_(np.all(np.array(np.transpose(A) == mA.T)))
assert_(np.all(np.array(np.transpose(A) == mA.H)))
assert_(np.all(A == mA.A))
B = A + 2j*A
mB = matrix(B)
assert_(np.allclose(linalg.inv(B), mB.I))
assert_(np.all(np.array(np.transpose(B) == mB.T)))
assert_(np.all(np.array(np.transpose(B).conj() == mB.H)))
def test_pinv(self):
x = matrix(np.arange(6).reshape(2, 3))
xpinv = matrix([[-0.77777778, 0.27777778],
[-0.11111111, 0.11111111],
[ 0.55555556, -0.05555556]])
assert_almost_equal(x.I, xpinv)
def test_comparisons(self):
A = np.arange(100).reshape(10, 10)
mA = matrix(A)
mB = matrix(A) + 0.1
assert_(np.all(mB == A+0.1))
assert_(np.all(mB == matrix(A+0.1)))
assert_(not np.any(mB == matrix(A-0.1)))
assert_(np.all(mA < mB))
assert_(np.all(mA <= mB))
assert_(np.all(mA <= mA))
assert_(not np.any(mA < mA))
assert_(not np.any(mB < mA))
assert_(np.all(mB >= mA))
assert_(np.all(mB >= mB))
assert_(not np.any(mB > mB))
assert_(np.all(mA == mA))
assert_(not np.any(mA == mB))
assert_(np.all(mB != mA))
assert_(not np.all(abs(mA) > 0))
assert_(np.all(abs(mB > 0)))
def test_asmatrix(self):
A = np.arange(100).reshape(10, 10)
mA = asmatrix(A)
A[0, 0] = -10
assert_(A[0, 0] == mA[0, 0])
def test_noaxis(self):
A = matrix([[1, 0], [0, 1]])
assert_(A.sum() == matrix(2))
assert_(A.mean() == matrix(0.5))
def test_repr(self):
A = matrix([[1, 0], [0, 1]])
assert_(repr(A) == "matrix([[1, 0],\n [0, 1]])")
def test_make_bool_matrix_from_str(self):
A = matrix('True; True; False')
B = matrix([[True], [True], [False]])
assert_array_equal(A, B)
class TestCasting:
def test_basic(self):
A = np.arange(100).reshape(10, 10)
mA = matrix(A)
mB = mA.copy()
O = np.ones((10, 10), np.float64) * 0.1
mB = mB + O
assert_(mB.dtype.type == np.float64)
assert_(np.all(mA != mB))
assert_(np.all(mB == mA+0.1))
mC = mA.copy()
O = np.ones((10, 10), np.complex128)
mC = mC * O
assert_(mC.dtype.type == np.complex128)
assert_(np.all(mA != mB))
class TestAlgebra:
def test_basic(self):
import numpy.linalg as linalg
A = np.array([[1., 2.], [3., 4.]])
mA = matrix(A)
B = np.identity(2)
for i in range(6):
assert_(np.allclose((mA ** i).A, B))
B = np.dot(B, A)
Ainv = linalg.inv(A)
B = np.identity(2)
for i in range(6):
assert_(np.allclose((mA ** -i).A, B))
B = np.dot(B, Ainv)
assert_(np.allclose((mA * mA).A, np.dot(A, A)))
assert_(np.allclose((mA + mA).A, (A + A)))
assert_(np.allclose((3*mA).A, (3*A)))
mA2 = matrix(A)
mA2 *= 3
assert_(np.allclose(mA2.A, 3*A))
def test_pow(self):
"""Test raising a matrix to an integer power works as expected."""
m = matrix("1. 2.; 3. 4.")
m2 = m.copy()
m2 **= 2
mi = m.copy()
mi **= -1
m4 = m2.copy()
m4 **= 2
assert_array_almost_equal(m2, m**2)
assert_array_almost_equal(m4, np.dot(m2, m2))
assert_array_almost_equal(np.dot(mi, m), np.eye(2))
def test_scalar_type_pow(self):
m = matrix([[1, 2], [3, 4]])
for scalar_t in [np.int8, np.uint8]:
two = scalar_t(2)
assert_array_almost_equal(m ** 2, m ** two)
def test_notimplemented(self):
'''Check that 'not implemented' operations produce a failure.'''
A = matrix([[1., 2.],
[3., 4.]])
with assert_raises(TypeError):
1.0**A
with assert_raises(TypeError):
A*object()
class TestMatrixReturn:
def test_instance_methods(self):
a = matrix([1.0], dtype='f8')
methodargs = {
'astype': ('intc',),
'clip': (0.0, 1.0),
'compress': ([1],),
'repeat': (1,),
'reshape': (1,),
'swapaxes': (0, 0),
'dot': np.array([1.0]),
}
excluded_methods = [
'argmin', 'choose', 'dump', 'dumps', 'fill', 'getfield',
'getA', 'getA1', 'item', 'nonzero', 'put', 'putmask', 'resize',
'searchsorted', 'setflags', 'setfield', 'sort',
'partition', 'argpartition', 'newbyteorder', 'to_device',
'take', 'tofile', 'tolist', 'tostring', 'tobytes', 'all', 'any',
'sum', 'argmax', 'argmin', 'min', 'max', 'mean', 'var', 'ptp',
'prod', 'std', 'ctypes', 'itemset', 'bitwise_count',
]
for attrib in dir(a):
if attrib.startswith('_') or attrib in excluded_methods:
continue
f = getattr(a, attrib)
if isinstance(f, collections.abc.Callable):
a.astype('f8')
a.fill(1.0)
if attrib in methodargs:
args = methodargs[attrib]
else:
args = ()
b = f(*args)
assert_(type(b) is matrix, "%s" % attrib)
assert_(type(a.real) is matrix)
assert_(type(a.imag) is matrix)
c, d = matrix([0.0]).nonzero()
assert_(type(c) is np.ndarray)
assert_(type(d) is np.ndarray)
class TestIndexing:
def test_basic(self):
x = asmatrix(np.zeros((3, 2), float))
y = np.zeros((3, 1), float)
y[:, 0] = [0.8, 0.2, 0.3]
x[:, 1] = y > 0.5
assert_equal(x, [[0, 1], [0, 0], [0, 0]])
class TestNewScalarIndexing:
a = matrix([[1, 2], [3, 4]])
def test_dimesions(self):
a = self.a
x = a[0]
assert_equal(x.ndim, 2)
def test_array_from_matrix_list(self):
a = self.a
x = np.array([a, a])
assert_equal(x.shape, [2, 2, 2])
def test_array_to_list(self):
a = self.a
assert_equal(a.tolist(), [[1, 2], [3, 4]])
def test_fancy_indexing(self):
a = self.a
x = a[1, [0, 1, 0]]
assert_(isinstance(x, matrix))
assert_equal(x, matrix([[3, 4, 3]]))
x = a[[1, 0]]
assert_(isinstance(x, matrix))
assert_equal(x, matrix([[3, 4], [1, 2]]))
x = a[[[1], [0]], [[1, 0], [0, 1]]]
assert_(isinstance(x, matrix))
assert_equal(x, matrix([[4, 3], [1, 2]]))
def test_matrix_element(self):
x = matrix([[1, 2, 3], [4, 5, 6]])
assert_equal(x[0][0], matrix([[1, 2, 3]]))
assert_equal(x[0][0].shape, (1, 3))
assert_equal(x[0].shape, (1, 3))
assert_equal(x[:, 0].shape, (2, 1))
x = matrix(0)
assert_equal(x[0, 0], 0)
assert_equal(x[0], 0)
assert_equal(x[:, 0].shape, x.shape)
def test_scalar_indexing(self):
x = asmatrix(np.zeros((3, 2), float))
assert_equal(x[0, 0], x[0][0])
def test_row_column_indexing(self):
x = asmatrix(np.eye(2))
assert_array_equal(x[0,:], [[1, 0]])
assert_array_equal(x[1,:], [[0, 1]])
assert_array_equal(x[:, 0], [[1], [0]])
assert_array_equal(x[:, 1], [[0], [1]])
def test_boolean_indexing(self):
A = np.arange(6)
A.shape = (3, 2)
x = asmatrix(A)
assert_array_equal(x[:, np.array([True, False])], x[:, 0])
assert_array_equal(x[np.array([True, False, False]),:], x[0,:])
def test_list_indexing(self):
A = np.arange(6)
A.shape = (3, 2)
x = asmatrix(A)
assert_array_equal(x[:, [1, 0]], x[:, ::-1])
assert_array_equal(x[[2,
def test_member_ravel(self):
assert_equal(self.a.ravel().shape, (2,))
assert_equal(self.m.ravel().shape, (1, 2))
def test_member_flatten(self):
assert_equal(self.a.flatten().shape, (2,))
assert_equal(self.m.flatten().shape, (1, 2))
def test_numpy_ravel_order(self):
x = np.array([[1, 2, 3], [4, 5, 6]])
assert_equal(np.ravel(x), [1, 2, 3, 4, 5, 6])
assert_equal(np.ravel(x, order='F'), [1, 4, 2, 5, 3, 6])
assert_equal(np.ravel(x.T), [1, 4, 2, 5, 3, 6])
assert_equal(np.ravel(x.T, order='A'), [1, 2, 3, 4, 5, 6])
x = matrix([[1, 2, 3], [4, 5, 6]])
assert_equal(np.ravel(x), [1, 2, 3, 4, 5, 6])
assert_equal(np.ravel(x, order='F'), [1, 4, 2, 5, 3, 6])
assert_equal(np.ravel(x.T), [1, 4, 2, 5, 3, 6])
assert_equal(np.ravel(x.T, order='A'), [1, 2, 3, 4, 5, 6])
def test_matrix_ravel_order(self):
x = matrix([[1, 2, 3], [4, 5, 6]])
assert_equal(x.ravel(), [[1, 2, 3, 4, 5, 6]])
assert_equal(x.ravel(order='F'), [[1, 4, 2, 5, 3, 6]])
assert_equal(x.T.ravel(), [[1, 4, 2, 5, 3, 6]])
assert_equal(x.T.ravel(order='A'), [[1, 2, 3, 4, 5, 6]])
def test_array_memory_sharing(self):
assert_(np.may_share_memory(self.a, self.a.ravel()))
assert_(not np.may_share_memory(self.a, self.a.flatten()))
def test_matrix_memory_sharing(self):
assert_(np.may_share_memory(self.m, self.m.ravel()))
assert_(not np.may_share_memory(self.m, self.m.flatten()))
def test_expand_dims_matrix(self):
a = np.arange(10).reshape((2, 5)).view(np.matrix)
expanded = np.expand_dims(a, axis=1)
assert_equal(expanded.ndim, 3)
assert_(not isinstance(expanded, np.matrix))
.\numpy\numpy\matrixlib\tests\test_interaction.py
"""
Tests of interaction of matrix with other parts of numpy.
Note that tests with MaskedArray and linalg are done in separate files.
"""
import pytest
import textwrap
import warnings
import numpy as np
from numpy.testing import (assert_, assert_equal, assert_raises,
assert_raises_regex, assert_array_equal,
assert_almost_equal, assert_array_almost_equal)
def test_fancy_indexing():
m = np.matrix([[1, 2], [3, 4]])
assert_(isinstance(m[[0, 1, 0], :], np.matrix))
x = np.asmatrix(np.arange(50).reshape(5, 10))
assert_equal(x[:2, np.array(-1)], x[:2, -1].T)
def test_polynomial_mapdomain():
dom1 = [0, 4]
dom2 = [1, 3]
x = np.matrix([dom1, dom1])
res = np.polynomial.polyutils.mapdomain(x, dom1, dom2)
assert_(isinstance(res, np.matrix))
def test_sort_matrix_none():
a = np.matrix([[2, 1, 0]])
actual = np.sort(a, axis=None)
expected = np.matrix([[0, 1, 2]])
assert_equal(actual, expected)
assert_(type(expected) is np.matrix)
def test_partition_matrix_none():
a = np.matrix([[2, 1, 0]])
actual = np.partition(a, 1, axis=None)
expected = np.matrix([[0, 1, 2]])
assert_equal(actual, expected)
assert_(type(expected) is np.matrix)
def test_dot_scalar_and_matrix_of_objects():
arr = np.matrix([1, 2], dtype=object)
desired = np.matrix([[3, 6]], dtype=object)
assert_equal(np.dot(arr, 3), desired)
assert_equal(np.dot(3, arr), desired)
def test_inner_scalar_and_matrix():
for dt in np.typecodes['AllInteger'] + np.typecodes['AllFloat'] + '?':
sca = np.array(3, dtype=dt)[()]
arr = np.matrix([[1, 2], [3, 4]], dtype=dt)
desired = np.matrix([[3, 6], [9, 12]], dtype=dt)
assert_equal(np.inner(arr, sca), desired)
assert_equal(np.inner(sca, arr), desired)
def test_inner_scalar_and_matrix_of_objects():
arr = np.matrix([1, 2], dtype=object)
desired = np.matrix([[3, 6]], dtype=object)
assert_equal(np.inner(arr, 3), desired)
assert_equal(np.inner(3, arr), desired)
def test_iter_allocate_output_subtype():
pass
a = np.matrix([[1, 2], [3, 4]])
b = np.arange(4).reshape(2, 2).T
i = np.nditer([a, b, None], [],
[['readonly'], ['readonly'], ['writeonly', 'allocate']])
assert_(type(i.operands[2]) is np.matrix)
assert_(type(i.operands[2]) is not np.ndarray)
assert_equal(i.operands[2].shape, (2, 2))
b = np.arange(4).reshape(1, 2, 2)
assert_raises(RuntimeError, np.nditer, [a, b, None], [],
[['readonly'], ['readonly'], ['writeonly', 'allocate']])
i = np.nditer([a, b, None], [],
[['readonly'], ['readonly'],
['writeonly', 'allocate', 'no_subtype']])
assert_(type(i.operands[2]) is np.ndarray)
assert_(type(i.operands[2]) is not np.matrix)
assert_equal(i.operands[2].shape, (1, 2, 2))
def like_function():
a = np.matrix([[1, 2], [3, 4]])
for like_function in np.zeros_like, np.ones_like, np.empty_like:
b = like_function(a)
assert_(type(b) is np.matrix)
c = like_function(a, subok=False)
assert_(type(c) is not np.matrix)
def test_array_astype():
a = np.matrix([[0, 1, 2], [3, 4, 5]], dtype='f4')
b = a.astype('f4', subok=True, copy=False)
assert_(a is b)
b = a.astype('i4', copy=False)
assert_equal(a, b)
assert_equal(type(b), np.matrix)
b = a.astype('f4', subok=False, copy=False)
assert_equal(a, b)
assert_(not (a is b))
assert_(type(b) is not np.matrix)
def test_stack():
m = np.matrix([[1, 2], [3, 4]])
assert_raises_regex(ValueError, 'shape too large to be a matrix',
np.stack, [m, m])
def test_object_scalar_multiply():
arr = np.matrix([1, 2], dtype=object)
desired = np.matrix([[3, 6]], dtype=object)
assert_equal(np.multiply(arr, 3), desired)
assert_equal(np.multiply(3, arr), desired)
def test_nanfunctions_matrices():
mat = np.matrix(np.eye(3))
for f in [np.nanmin, np.nanmax]:
res = f(mat, axis=0)
assert_(isinstance(res, np.matrix))
assert_(res.shape == (1, 3))
res = f(mat, axis=1)
assert_(isinstance(res, np.matrix))
assert_(res.shape == (3, 1))
res = f(mat)
assert_(np.isscalar(res))
mat[1] = np.nan
for f in [np.nanmin, np.nanmax]:
with warnings.catch_warnings(record=True) as w:
warnings.simplefilter('always')
res = f(mat, axis=0)
assert_(isinstance(res, np.matrix))
assert_(not np.any(np.isnan(res)))
assert_(len(w) == 0)
with warnings.catch_warnings(record=True) as w:
warnings.simplefilter('always')
res = f(mat, axis=1)
assert_(isinstance(res, np.matrix))
assert_(np.isnan(res[1, 0]) and not np.isnan(res[0, 0])
and not np.isnan(res[2, 0]))
assert_(len(w) == 1, 'no warning raised')
assert_(issubclass(w[0].category, RuntimeWarning))
with warnings.catch_warnings(record=True) as w:
warnings.simplefilter('always')
res = f(mat)
assert_(np.isscalar(res))
assert_(res != np.nan)
assert_(len(w) == 0)
def test_nanfunctions_matrices_general():
mat = np.matrix(np.eye(3))
for f in (np.nanargmin, np.nanargmax, np.nansum, np.nanprod,
np.nanmean, np.nanvar, np.nanstd):
res = f(mat, axis=0)
assert_(isinstance(res, np.matrix))
assert_(res.shape == (1, 3))
res = f(mat, axis=1)
assert_(isinstance(res, np.matrix))
assert_(res.shape == (3, 1))
res = f(mat)
assert_(np.isscalar(res))
for f in np.nancumsum, np.nancumprod:
res = f(mat, axis=0)
assert_(isinstance(res, np.matrix))
assert_(res.shape == (3, 3))
res = f(mat, axis=1)
assert_(isinstance(res, np.matrix))
assert_(res.shape == (3, 3))
res = f(mat)
assert_(isinstance(res, np.matrix))
assert_(res.shape == (1, 3*3))
def test_average_matrix():
y = np.matrix(np.random.rand(5, 5))
assert_array_equal(y.mean(0), np.average(y, 0))
a = np.matrix([[1, 2], [3, 4]])
w = np.matrix([[1, 2], [3, 4]])
r = np.average(a, axis=0, weights=w)
assert_equal(type(r), np.matrix)
assert_equal(r, [[2.5, 10.0/3]])
def test_dot_matrix():
x = np.linspace(0, 5)
y = np.linspace(-5, 0)
mx = np.matrix(x)
my = np.matrix(y)
r = np.dot(x, y)
mr = np.dot(mx, my.T)
assert_almost_equal(mr, r)
def test_ediff1d_matrix():
assert(isinstance(np.ediff1d(np.matrix(1)), np.matrix))
assert(isinstance(np.ediff1d(np.matrix(1), to_begin=1), np.matrix))
def test_apply_along_axis_matrix():
def double(row):
return row * 2
m = np.matrix([[0, 1], [2, 3]])
expected = np.matrix([[0, 2], [4, 6]])
result = np.apply_along_axis(double, 0, m)
assert_(isinstance(result, np.matrix))
assert_array_equal(result, expected)
result = np.apply_along_axis(double, 1, m)
assert_(isinstance(result, np.matrix))
assert_array_equal(result, expected)
def test_kron_matrix():
a = np.ones([2, 2])
m = np.asmatrix(a)
assert_equal(type(np.kron(a, a)), np.ndarray)
assert_equal(type(np.kron(m, m)), np.matrix)
assert_equal(type(np.kron(a, m)), np.matrix)
assert_equal(type(np.kron(m, a)), np.matrix)
class TestConcatenatorMatrix:
def test_matrix(self):
a = [1, 2]
b = [3, 4]
ab_r = np.r_['r', a, b]
ab_c = np.r_['c', a, b]
assert_equal(type(ab_r), np.matrix)
assert_equal(type(ab_c), np.matrix)
assert_equal(np.array(ab_r), [[1, 2, 3, 4]])
assert_equal(np.array(ab_c), [[1], [2], [3], [4]])
assert_raises(ValueError, lambda: np.r_['rc', a, b])
def test_matrix_scalar(self):
r = np.r_['r', [1, 2], 3]
assert_equal(type(r), np.matrix)
assert_equal(np.array(r), [[1, 2, 3]])
def test_matrix_builder(self):
a = np.array([1])
b = np.array([2])
c = np.array([3])
d = np.array([4])
actual = np.r_['a, b; c, d']
expected = np.bmat([[a, b], [c, d]])
assert_equal(actual, expected)
assert_equal(type(actual), type(expected))
def test_array_equal_error_message_matrix():
with pytest.raises(AssertionError) as exc_info:
assert_equal(np.array([1, 2]), np.matrix([1, 2]))
msg = str(exc_info.value)
msg_reference = textwrap.dedent("""\
Arrays are not equal
(shapes (2,), (1, 2) mismatch)
ACTUAL: array([1, 2])
DESIRED: matrix([[1, 2]])""")
assert_equal(msg, msg_reference)
def test_array_almost_equal_matrix():
m1 = np.matrix([[1., 2.]])
m2 = np.matrix([[1., np.nan]])
m3 = np.matrix([[1., -np.inf]])
m4 = np.matrix([[np.nan, np.inf]])
m5 = np.matrix([[1., 2.], [np.nan, np.inf]])
for assert_func in assert_array_almost_equal, assert_almost_equal:
for m in m1, m2, m3, m4, m5:
assert_func(m, m)
a = np.array(m)
assert_func(a, m)
assert_func(m, a)
.\numpy\numpy\matrixlib\tests\test_masked_matrix.py
import pickle
import numpy as np
from numpy.testing import assert_warns
from numpy.ma.testutils import (assert_, assert_equal, assert_raises,
assert_array_equal)
from numpy.ma.core import (masked_array, masked_values, masked, allequal,
MaskType, getmask, MaskedArray, nomask,
log, add, hypot, divide)
from numpy.ma.extras import mr_
class MMatrix(MaskedArray, np.matrix,):
def __new__(cls, data, mask=nomask):
mat = np.matrix(data)
_data = MaskedArray.__new__(cls, data=mat, mask=mask)
return _data
def __array_finalize__(self, obj):
np.matrix.__array_finalize__(self, obj)
MaskedArray.__array_finalize__(self, obj)
return
@property
def _series(self):
_view = self.view(MaskedArray)
_view._sharedmask = False
return _view
class TestMaskedMatrix:
def test_matrix_indexing(self):
x1 = np.matrix([[1, 2, 3], [4, 3, 2]])
x2 = masked_array(x1, mask=[[1, 0, 0], [0, 1, 0]])
x3 = masked_array(x1, mask=[[0, 1, 0], [1, 0, 0]])
x4 = masked_array(x1)
str(x2)
repr(x2)
assert_(type(x2[1, 0]) is type(x1[1, 0]))
assert_(x1[1, 0] == x2[1, 0])
assert_(x2[1, 1] is masked)
assert_equal(x1[0, 2], x2[0, 2])
assert_equal(x1[0, 1:], x2[0, 1:])
assert_equal(x1[:, 2], x2[:, 2])
assert_equal(x1[:], x2[:])
assert_equal(x1[1:], x3[1:])
x1[0, 2] = 9
x2[0, 2] = 9
assert_equal(x1, x2)
x1[0, 1:] = 99
x2[0, 1:] = 99
assert_equal(x1, x2)
x2[0, 1] = masked
assert_equal(x1, x2)
x2[0, 1:] = masked
assert_equal(x1, x2)
x2[0, :] = x1[0, :]
x2[0, 1] = masked
assert_(allequal(getmask(x2), np.array([[0, 1, 0], [0, 1, 0]])))
x3[1, :] = masked_array([1, 2, 3], [1, 1, 0])
assert_(allequal(getmask(x3)[1], masked_array([1, 1, 0])))
assert_(allequal(getmask(x3[1]), masked_array([1, 1, 0])))
x4[1, :] = masked_array([1, 2, 3], [1, 1, 0])
assert_(allequal(getmask(x4[1]), masked_array([1, 1, 0])))
assert_(allequal(x4[1], masked_array([1, 2, 3])))
x1 = np.matrix(np.arange(5) * 1.0)
x2 = masked_values(x1, 3.0)
assert_equal(x1, x2)
assert_(allequal(masked_array([0, 0, 0, 1, 0], dtype=MaskType),
x2.mask))
assert_equal(3.0, x2.fill_value)
def test_pickling_subbaseclass(self):
a = masked_array(np.matrix(list(range(10))), mask=[1, 0, 1, 0, 0] * 2)
for proto in range(2, pickle.HIGHEST_PROTOCOL + 1):
a_pickled = pickle.loads(pickle.dumps(a, protocol=proto))
assert_equal(a_pickled._mask, a._mask)
assert_equal(a_pickled, a)
assert_(isinstance(a_pickled._data, np.matrix))
def test_count_mean_with_matrix(self):
m = masked_array(np.matrix([[1, 2], [3, 4]]), mask=np.zeros((2, 2)))
assert_equal(m.count(axis=0).shape, (1, 2))
assert_equal(m.count(axis=1).shape, (2, 1))
assert_equal(m.mean(axis=0), [[2., 3.]])
assert_equal(m.mean(axis=1), [[1.5], [3.5]])
def test_flat(self):
test = masked_array(np.matrix([[1, 2, 3]]), mask=[0, 0, 1])
assert_equal(test.flat[1], 2)
assert_equal(test.flat[2], masked)
assert_(np.all(test.flat[0:2] == test[0, 0:2]))
test = masked_array(np.matrix([[1, 2, 3]]), mask=[0, 0, 1])
test.flat = masked_array([3, 2, 1], mask=[1, 0, 0])
control = masked_array(np.matrix([[3, 2, 1]]), mask=[1, 0, 0])
assert_equal(test, control)
test = masked_array(np.matrix([[1, 2, 3]]), mask=[0, 0, 1])
testflat = test.flat
testflat[:] = testflat[[2, 1, 0]]
assert_equal(test, control)
testflat[0] = 9
a = masked_array(np.matrix(np.eye(2)), mask=0)
b = a.flat
b01 = b[:2]
assert_equal(b01.data, np.array([[1., 0.]]))
assert_equal(b01.mask, np.array([[False, False]]))
def test_allany_onmatrices(self):
x = np.array([[0.13, 0.26, 0.90],
[0.28, 0.33, 0.63],
[0.31, 0.87, 0.70]])
X = np.matrix(x)
m = np.array([[True, False, False],
[False, False, False],
[True, True, False]], dtype=np.bool)
mX = masked_array(X, mask=m)
mXbig = (mX > 0.5)
mXsmall = (mX < 0.5)
assert_(not mXbig.all())
assert_(mXbig.any())
assert_equal(mXbig.all(0), np.matrix([False, False, True]))
assert_equal(mXbig.all(1), np.matrix([False, False, True]).T)
assert_equal(mXbig.any(0), np.matrix([False, False, True]))
assert_equal(mXbig.any(1), np.matrix([True, True, True]).T)
assert_(not mXsmall.all())
assert_(mXsmall.any())
assert_equal(mXsmall.all(0), np.matrix([True, True, False]))
assert_equal(mXsmall.all(1), np.matrix([False, False, False]).T)
assert_equal(mXsmall.any(0), np.matrix([True, True, False]))
assert_equal(mXsmall.any(1), np.matrix([True, True, False]).T)
def test_compressed(self):
a = masked_array(np.matrix([1, 2, 3, 4]), mask=[0, 0, 0, 0])
b = a.compressed()
assert_equal(b, a)
assert_(isinstance(b, np.matrix))
a[0, 0] = masked
b = a.compressed()
assert_equal(b, [[2, 3, 4]])
def test_ravel(self):
a = masked_array(np.matrix([1, 2, 3, 4, 5]), mask=[[0, 1, 0, 0, 0]])
aravel = a.ravel()
assert_equal(aravel.shape, (1, 5))
assert_equal(aravel._mask.shape, a.shape)
def test_view(self):
iterator = list(zip(np.arange(10), np.random.rand(10)))
data = np.array(iterator)
a = masked_array(iterator, dtype=[('a', float), ('b', float)])
a.mask[0] = (1, 0)
test = a.view((float, 2), np.matrix)
assert_equal(test, data)
assert_(isinstance(test, np.matrix))
assert_(not isinstance(test, MaskedArray))
class TestSubclassing:
def setup_method(self):
x = np.arange(5, dtype='float')
mx = MMatrix(x, mask=[0, 1, 0, 0, 0])
self.data = (x, mx)
def test_maskedarray_subclassing(self):
(x, mx) = self.data
assert_(isinstance(mx._data, np.matrix))
def test_masked_unary_operations(self):
(x, mx) = self.data
with np.errstate(divide='ignore'):
assert_(isinstance(log(mx), MMatrix))
assert_equal(log(x), np.log(x))
def test_masked_binary_operations(self):
(x, mx) = self.data
assert_(isinstance(add(mx, mx), MMatrix))
assert_(isinstance(add(mx, x), MMatrix))
assert_equal(add(mx, x), mx+x)
assert_(isinstance(add(mx, mx)._data, np.matrix))
with assert_warns(DeprecationWarning):
assert_(isinstance(add.outer(mx, mx), MMatrix))
assert_(isinstance(hypot(mx, mx), MMatrix))
assert_(isinstance(hypot(mx, x), MMatrix))
def test_masked_binary_operations2(self):
(x, mx) = self.data
xmx = masked_array(mx.data.__array__(), mask=mx.mask)
assert_(isinstance(divide(mx, mx), MMatrix))
assert_(isinstance(divide(mx, x), MMatrix))
assert_equal(divide(mx, mx), divide(xmx, xmx))
class TestConcatenator:
def test_matrix_builder(self):
assert_raises(np.ma.MAError, lambda: mr_['1, 2; 3, 4'])
def test_matrix(self):
actual = mr_['r', 1, 2, 3]
expected = np.ma.array(np.r_['r', 1, 2, 3])
assert_array_equal(actual, expected)
assert_equal(type(actual), type(expected))
assert_equal(type(actual.data), type(expected.data))
.\numpy\numpy\matrixlib\tests\test_matrix_linalg.py
import numpy as np
from numpy.linalg.tests.test_linalg import (
LinalgCase, apply_tag, TestQR as _TestQR, LinalgTestCase,
_TestNorm2D, _TestNormDoubleBase, _TestNormSingleBase, _TestNormInt64Base,
SolveCases, InvCases, EigvalsCases, EigCases, SVDCases, CondCases,
PinvCases, DetCases, LstsqCases)
CASES = []
CASES += apply_tag('square', [
LinalgCase("0x0_matrix",
np.empty((0, 0), dtype=np.double).view(np.matrix),
np.empty((0, 1), dtype=np.double).view(np.matrix),
tags={'size-0'}),
LinalgCase("matrix_b_only",
np.array([[1., 2.], [3., 4.]]),
np.matrix([2., 1.]).T),
LinalgCase("matrix_a_and_b",
np.matrix([[1., 2.], [3., 4.]]),
np.matrix([2., 1.]).T),
])
CASES += apply_tag('hermitian', [
LinalgCase("hmatrix_a_and_b",
np.matrix([[1., 2.], [2., 1.]]),
None),
])
class MatrixTestCase(LinalgTestCase):
TEST_CASES = CASES
class TestSolveMatrix(SolveCases, MatrixTestCase):
pass
class TestInvMatrix(InvCases, MatrixTestCase):
pass
class TestEigvalsMatrix(EigvalsCases, MatrixTestCase):
pass
class TestEigMatrix(EigCases, MatrixTestCase):
pass
class TestSVDMatrix(SVDCases, MatrixTestCase):
pass
class TestCondMatrix(CondCases, MatrixTestCase):
pass
class TestPinvMatrix(PinvCases, MatrixTestCase):
pass
class TestDetMatrix(DetCases, MatrixTestCase):
pass
class TestLstsqMatrix(LstsqCases, MatrixTestCase):
pass
class _TestNorm2DMatrix(_TestNorm2D):
array = np.matrix
class TestNormDoubleMatrix(_TestNorm2DMatrix, _TestNormDoubleBase):
pass
class TestNormSingleMatrix(_TestNorm2DMatrix, _TestNormSingleBase):
pass
class TestNormInt64Matrix(_TestNorm2DMatrix, _TestNormInt64Base):
pass
class TestQRMatrix(_TestQR):
array = np.matrix
.\numpy\numpy\matrixlib\tests\test_multiarray.py
import numpy as np
from numpy.testing import assert_, assert_equal, assert_array_equal
class TestView:
def test_type(self):
x = np.array([1, 2, 3])
assert_(isinstance(x.view(np.matrix), np.matrix))
def test_keywords(self):
x = np.array([(1, 2)], dtype=[('a', np.int8), ('b', np.int8)])
y = x.view(dtype='<i2', type=np.matrix)
assert_array_equal(y, [[513]])
assert_(isinstance(y, np.matrix))
assert_equal(y.dtype, np.dtype('<i2'))
.\numpy\numpy\matrixlib\tests\test_numeric.py
import numpy as np
from numpy.testing import assert_equal
class TestDot:
def test_matscalar(self):
b1 = np.matrix(np.ones((3, 3), dtype=complex))
assert_equal(b1*1.0, b1)
def test_diagonal():
b1 = np.matrix([[1,2],[3,4]])
diag_b1 = np.matrix([[1, 4]])
array_b1 = np.array([1, 4])
assert_equal(b1.diagonal(), diag_b1)
assert_equal(np.diagonal(b1), array_b1)
assert_equal(np.diag(b1), array_b1)
.\numpy\numpy\matrixlib\tests\test_regression.py
import numpy as np
from numpy.testing import assert_, assert_equal, assert_raises
class TestRegression:
def test_kron_matrix(self):
x = np.matrix('[1 0; 1 0]')
assert_equal(type(np.kron(x, x)), type(x))
def test_matrix_properties(self):
a = np.matrix([1.0], dtype=float)
assert_(type(a.real) is np.matrix)
assert_(type(a.imag) is np.matrix)
c, d = np.matrix([0.0]).nonzero()
assert_(type(c) is np.ndarray)
assert_(type(d) is np.ndarray)
def test_matrix_multiply_by_1d_vector(self):
def mul():
np.asmatrix(np.eye(2))*np.ones(2)
assert_raises(ValueError, mul)
def test_matrix_std_argmax(self):
x = np.asmatrix(np.random.uniform(0, 1, (3, 3)))
assert_equal(x.std().shape, ())
assert_equal(x.argmax().shape, ())
.\numpy\numpy\matrixlib\tests\__init__.py
def merge_sort(arr):
if len(arr) <= 1:
return arr
mid = len(arr) // 2
left_half = merge_sort(arr[:mid])
right_half = merge_sort(arr[mid:])
return merge(left_half, right_half)
def merge(left, right):
result = []
i = j = 0
while i < len(left) and j < len(right):
if left[i] < right[j]:
result.append(left[i])
i += 1
else:
result.append(right[j])
j += 1
result.extend(left[i:])
result.extend(right[j:])
return result
.\numpy\numpy\matrixlib\__init__.py
"""Sub-package containing the matrix class and related functions.
"""
from . import defmatrix
from .defmatrix import *
__all__ = defmatrix.__all__
from numpy._pytesttester import PytestTester
test = PytestTester(__name__)
del PytestTester
.\numpy\numpy\matrixlib\__init__.pyi
from numpy._pytesttester import PytestTester
from numpy import (
matrix as matrix,
)
from numpy.matrixlib.defmatrix import (
bmat as bmat,
mat as mat,
asmatrix as asmatrix,
)
__all__: list[str]
test: PytestTester
.\numpy\numpy\polynomial\chebyshev.py
def _cseries_to_zseries(c):
"""Convert Chebyshev series to z-series.
Convert a Chebyshev series to the equivalent z-series. The result is
obtained by applying the algebraic identity involving z, transforming
the Chebyshev series coefficients to their corresponding z-series
coefficients.
"""
return (c + np.flip(c)) / 2
n = c.size
zs = np.zeros(2*n-1, dtype=c.dtype)
zs[n-1:] = c/2
return zs + zs[::-1]
def _zseries_to_cseries(zs):
"""Convert z-series to a Chebyshev series.
Convert a z series to the equivalent Chebyshev series. The result is
never an empty array. The dtype of the return is the same as that of
the input. No checks are run on the arguments as this routine is for
internal use.
Parameters
----------
zs : 1-D ndarray
Odd length symmetric z-series, ordered from low to high.
Returns
-------
c : 1-D ndarray
Chebyshev coefficients, ordered from low to high.
"""
n = (zs.size + 1)//2
c = zs[n-1:].copy()
c[1:n] *= 2
return c
def _zseries_mul(z1, z2):
"""Multiply two z-series.
Multiply two z-series to produce a z-series.
Parameters
----------
z1, z2 : 1-D ndarray
The arrays must be 1-D but this is not checked.
Returns
-------
product : 1-D ndarray
The product z-series.
Notes
-----
This is simply convolution. If symmetric/anti-symmetric z-series are
denoted by S/A then the following rules apply:
S*S, A*A -> S
S*A, A*S -> A
"""
return np.convolve(z1, z2)
def _zseries_div(z1, z2):
"""Divide the first z-series by the second.
Divide `z1` by `z2` and return the quotient and remainder as z-series.
Warning: this implementation only applies when both z1 and z2 have the
same symmetry, which is sufficient for present purposes.
Parameters
----------
z1, z2 : 1-D ndarray
The arrays must be 1-D and have the same symmetry, but this is not
checked.
Returns
-------
(quotient, remainder) : 1-D ndarrays
Quotient and remainder as z-series.
Notes
-----
This is not the same as polynomial division on account of the desired form
of the remainder. If symmetric/anti-symmetric z-series are denoted by S/A
then the following rules apply:
S/S -> S,S
A/A -> S,A
The restriction to types of the same symmetry could be fixed but seems like
unneeded generality. There is no natural form for the remainder in the case
where there is no symmetry.
"""
z1 = z1.copy()
z2 = z2.copy()
lc1 = len(z1)
lc2 = len(z2)
if lc2 == 1:
z1 /= z2
return z1, z1[:1]*0
elif lc1 < lc2:
return z1[:1]*0, z1
else:
dlen = lc1 - lc2
scl = z2[0]
z2 /= scl
quo = np.empty(dlen + 1, dtype=z1.dtype)
i = 0
j = dlen
while i < j:
r = z1[i]
quo[i] = z1[i]
quo[dlen - i] = r
tmp = r*z2
z1[i:i+lc2] -= tmp
z1[j:j+lc2] -= tmp
i += 1
j -= 1
r = z1[i]
quo[i] = r
tmp = r*z2
z1[i:i+lc2] -= tmp
quo /= scl
rem = z1[i+1:i-1+lc2].copy()
return quo, rem
def _zseries_der(zs):
"""Differentiate a z-series.
The derivative is with respect to x, not z. This is achieved using the
"""
pass
Parameters
----------
zs : z-series
待求导的 z-series。
Returns
-------
derivative : z-series
导数结果。
Notes
-----
对于 x 的 z-series (ns)，已乘以二以避免使用与 Decimal 和其他特定标量类型不兼容的浮点数。
为了补偿这种缩放效果，zs 的值也乘以二，这样两者在除法中相互抵消。
"""
# 计算 ns 的长度并取整除以 2
n = len(zs) // 2
# 创建一个 numpy 数组 ns，包含 [-1, 0, 1]，数据类型与 zs 相同
ns = np.array([-1, 0, 1], dtype=zs.dtype)
# 将 zs 的每个元素乘以对应位置的 (-n 到 n 的范围) * 2，以补偿之前的乘法操作
zs *= np.arange(-n, n+1) * 2
# 使用 _zseries_div 函数计算 zs 和 ns 的商，结果包括商 d 和余数 r
d, r = _zseries_div(zs, ns)
# 返回导数结果 d
return d
# 将 z-series 进行积分处理，积分是相对于 x 而非 z。这是通过变量变换 dx/dz 来实现的，变换细节见模块注释。
def _zseries_int(zs):
# 计算 z-series 的长度，并计算对应的 x-series 的长度
n = 1 + len(zs)//2
# 创建一个包含 [-1, 0, 1] 的 NumPy 数组，数据类型与 zs 相同
ns = np.array([-1, 0, 1], dtype=zs.dtype)
# 将 zs 乘以 x-series 的系数，以进行变量变换
zs = _zseries_mul(zs, ns)
# 创建一个数组，其值为 [-n*2, ..., -2, 0, 2, ..., n*2]，用于后续的除法
div = np.arange(-n, n+1)*2
# 对 zs 进行分段除法，以完成变量变换后的调整
zs[:n] /= div[:n]
zs[n+1:] /= div[n+1:]
# 将中间值置为 0，符合积分常数的定义
zs[n] = 0
# 返回处理后的 z-series，即 x-series 的积分结果
return zs
#
# Chebyshev series functions
#
def poly2cheb(pol):
"""
Convert a polynomial to a Chebyshev 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 Chebyshev 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 Chebyshev
series.
See Also
--------
cheb2poly
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(range(4))
>>> p
Polynomial([0., 1., 2., 3.], domain=[-1., 1.], window=[-1., 1.], symbol='x')
>>> c = p.convert(kind=P.Chebyshev)
>>> c
Chebyshev([1. , 3.25, 1. , 0.75], domain=[-1., 1.], window=[-1., ...
>>> P.chebyshev.poly2cheb(range(4))
array([1. , 3.25, 1. , 0.75])
"""
# 使用 pu.as_series 将输入 pol 转换为 series
[pol] = pu.as_series([pol])
# 计算多项式的最高阶数
deg = len(pol) - 1
# 初始化结果为空的 Chebyshev 系数数组
res = 0
# 从最高阶开始逐步构建 Chebyshev 系数
for i in range(deg, -1, -1):
# 将当前多项式系数转换为 Chebyshev 系数并累加到 res 中
res = chebadd(chebmulx(res), pol[i])
# 返回转换后的 Chebyshev 系数数组
return res
def cheb2poly(c):
"""
Convert a Chebyshev series to a polynomial.
Convert an array representing the coefficients of a Chebyshev 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 Chebyshev 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
--------
poly2cheb
Notes
-----
"""
# 该函数未实现详细注释，仅提供函数的基本结构和说明
pass
# 导入必要的函数和类来执行多项式操作
from .polynomial import polyadd, polysub, polymulx
# 将输入的多项式系数转换为系列（series），即确保 c 是一个多项式系数列表
[c] = pu.as_series([c])
# 计算多项式的阶数
n = len(c)
# 如果多项式阶数小于3，则直接返回多项式系数 c，不需要进行转换操作
if n < 3:
return c
else:
# 从 c 中提取倒数第二个系数和最后一个系数作为初始值
c0 = c[-2]
c1 = c[-1]
# 从高阶到低阶（从 n-1 到 2），依次进行转换操作
for i in range(n - 1, 1, -1):
tmp = c0
# 计算当前阶数 c1 的多项式转换结果
c0 = polysub(c[i - 2], c1)
c1 = polyadd(tmp, polymulx(c1)*2)
# 返回转换后的多项式系数列表
return polyadd(c0, polymulx(c1))
# 这些是整数类型的常量数组，以便与广泛的其他类型兼容，比如 Decimal。
# Chebyshev 默认定义域。
chebdomain = np.array([-1., 1.])
# 表示零的 Chebyshev 系数。
chebzero = np.array([0])
# 表示一的 Chebyshev 系数。
chebone = np.array([1])
# 表示恒等函数 x 的 Chebyshev 系数。
chebx = np.array([0, 1])
def chebline(off, scl):
"""
表示一条直线的 Chebyshev 级数。
Parameters
----------
off, scl : scalar
指定的直线由 ``off + scl*x`` 给出。
Returns
-------
y : ndarray
该模块对于 ``off + scl*x`` 的 Chebyshev 级数表示。
See Also
--------
numpy.polynomial.polynomial.polyline
numpy.polynomial.legendre.legline
numpy.polynomial.laguerre.lagline
numpy.polynomial.hermite.hermline
numpy.polynomial.hermite_e.hermeline
Examples
--------
>>> import numpy.polynomial.chebyshev as C
>>> C.chebline(3,2)
array([3, 2])
>>> C.chebval(-3, C.chebline(3,2))
-3.0
"""
if scl != 0:
return np.array([off, scl])
else:
return np.array([off])
def chebfromroots(roots):
"""
生成具有给定根的 Chebyshev 级数。
函数返回多项式
.. math:: p(x) = (x - r_0) * (x - r_1) * ... * (x - r_n)
的 Chebyshev 形式的系数，其中 :math:`r_n` 是 `roots` 中指定的根。如果一个零有多重性 n，
那么它必须在 `roots` 中出现 n 次。例如，如果 2 是多重性为三的根，3 是多重性为两的根，
那么 `roots` 看起来像是 [2, 2, 2, 3, 3]。根可以以任何顺序出现。
如果返回的系数是 `c`，则
.. math:: p(x) = c_0 + c_1 * T_1(x) + ... + c_n * T_n(x)
最后一项的系数通常不是 Chebyshev 形式的单项多项式的 1。
Parameters
----------
roots : array_like
包含根的序列。
Returns
-------
out : ndarray
系数的一维数组。如果所有根都是实数，则 `out` 是实数组；如果某些根是复数，则 `out` 是复数，
即使结果中的所有系数都是实数（见下面的示例）。
See Also
--------
numpy.polynomial.polynomial.polyfromroots
numpy.polynomial.legendre.legfromroots
numpy.polynomial.laguerre.lagfromroots
numpy.polynomial.hermite.hermfromroots
numpy.polynomial.hermite_e.hermefromroots
Examples
--------
>>> import numpy.polynomial.chebyshev as C
>>> C.chebfromroots((-1,0,1))
array([ 0. , -0.25, 0. , 0.25])
>>> j = complex(0,1)
>>> C.chebfromroots((-j,j))
array([1.5+0.j, 0. +0.j, 0.5+0.j])
"""
# 调用 pu 模块的 _fromroots 函数，传入 chebline、chebmul、roots 作为参数，并返回其结果
return pu._fromroots(chebline, chebmul, roots)
# 定义函数 chebadd，用于将两个切比雪夫级数相加
def chebadd(c1, c2):
# 调用私有函数 _add，执行切比雪夫级数的加法操作
return pu._add(c1, c2)
# 定义函数 chebsub，用于从一个切比雪夫级数中减去另一个切比雪夫级数
def chebsub(c1, c2):
# 调用私有函数 _sub，执行切比雪夫级数的减法操作
return pu._sub(c1, c2)
# 定义函数 chebmulx，用于将一个切比雪夫级数乘以自变量 x
def chebmulx(c):
"""
Multiply a Chebyshev series by x.
Multiply the polynomial `c` by x, where x is the independent
variable.
Parameters
----------
c : array_like
1-D array of Chebyshev series coefficients ordered from low to
high.
Returns
-------
out : ndarray
Array representing the result of the multiplication.
See Also
--------
chebadd, chebsub, chebmul, chebdiv, chebpow
Notes
-----
.. versionadded:: 1.5.0
Examples
--------
>>> from numpy.polynomial import chebyshev as C
>>> C.chebmulx([1,2,3])
array([1. , 2.5, 1. , 1.5])
"""
# 调用 pu.as_series([c]) 方法返回的结果，取其第一个元素作为 c
[c] = pu.as_series([c])
# 如果输入列表 c 的长度为1且唯一元素为0，则直接返回该列表，不进行后续计算
if len(c) == 1 and c[0] == 0:
return c
# 创建一个长度为 len(c)+1 的空数组 prd，其数据类型与 c 相同
prd = np.empty(len(c) + 1, dtype=c.dtype)
# 将 prd 数组的第一个元素设置为 c[0] 的零乘积
prd[0] = c[0]*0
# 将 prd 数组的第二个元素设置为 c[0]
prd[1] = c[0]
# 如果 c 的长度大于1，则进行以下操作
if len(c) > 1:
# 将 c 的第二个元素及其后的元素都除以2，存储到 tmp 变量中
tmp = c[1:]/2
# 将 tmp 中的值赋给 prd 数组的第三个元素及其后的元素
prd[2:] = tmp
# 将 tmp 中的值累加到 prd 数组的倒数第二个元素及其前的元素上
prd[0:-2] += tmp
# 返回计算结果数组 prd
return prd
# Multiply one Chebyshev series by another.
# 返回两个切比雪夫级数的乘积。
def chebmul(c1, c2):
"""
Multiply one Chebyshev series by another.
Returns the product of two Chebyshev series `c1` * `c2`. The arguments
are sequences of coefficients, from lowest order "term" to highest,
e.g., [1,2,3] represents the series ``T_0 + 2*T_1 + 3*T_2``.
Parameters
----------
c1, c2 : array_like
1-D arrays of Chebyshev series coefficients ordered from low to
high.
Returns
-------
out : ndarray
Of Chebyshev series coefficients representing their product.
See Also
--------
chebadd, chebsub, chebmulx, chebdiv, chebpow
Notes
-----
In general, the (polynomial) product of two C-series results in terms
that are not in the Chebyshev polynomial basis set. Thus, to express
the product as a C-series, it is typically necessary to "reproject"
the product onto said basis set, which typically produces
"unintuitive live" (but correct) results; see Examples section below.
Examples
--------
>>> from numpy.polynomial import chebyshev as C
>>> c1 = (1,2,3)
>>> c2 = (3,2,1)
>>> C.chebmul(c1,c2)
array([ 6.5, 12. , 12. , 4. , 1.5])
"""
# Trim c1, c2 to ensure they are valid series
[c1, c2] = pu.as_series([c1, c2])
# Convert c1 to a z-series (series of complex numbers)
z1 = _cseries_to_zseries(c1)
# Convert c2 to a z-series
z2 = _cseries_to_zseries(c2)
# Multiply the z-series z1 and z2
prd = _zseries_mul(z1, z2)
# Convert the resulting z-series back to a Chebyshev series
ret = _zseries_to_cseries(prd)
# Trim any unnecessary zeros from the resulting series and return
return pu.trimseq(ret)
# Divide one Chebyshev series by another.
# 返回两个切比雪夫级数的商和余数。
def chebdiv(c1, c2):
"""
Divide one Chebyshev series by another.
Returns the quotient-with-remainder of two Chebyshev series
`c1` / `c2`. The arguments are sequences of coefficients from lowest
order "term" to highest, e.g., [1,2,3] represents the series
``T_0 + 2*T_1 + 3*T_2``.
Parameters
----------
c1, c2 : array_like
1-D arrays of Chebyshev series coefficients ordered from low to
high.
Returns
-------
[quo, rem] : ndarrays
Of Chebyshev series coefficients representing the quotient and
remainder.
See Also
--------
chebadd, chebsub, chebmulx, chebmul, chebpow
Notes
-----
In general, the (polynomial) division of one C-series by another
results in quotient and remainder terms that are not in the Chebyshev
polynomial basis set. Thus, to express these results as C-series, it
is typically necessary to "reproject" the results onto said basis
set, which typically produces "unintuitive" (but correct) results;
see Examples section below.
Examples
--------
>>> from numpy.polynomial import chebyshev as C
>>> c1 = (1,2,3)
>>> c2 = (3,2,1)
>>> C.chebdiv(c1,c2)
(array([3.]), array([-8., -4.]))
>>> c2 = (0,1,2,3)
>>> C.chebdiv(c2,c1)
(array([0., 2.]), array([-2., -4.]))
"""
# Trim c1, c2 to ensure they are valid series
[c1, c2] = pu.as_series([c1, c2])
# 如果 c2 的最后一个元素为 0，则抛出 ZeroDivisionError 异常
if c2[-1] == 0:
raise ZeroDivisionError()
# 注意：这种写法比 `pu._div(chebmul, c1, c2)` 更有效率
# 计算 c1 和 c2 的长度
lc1 = len(c1)
lc2 = len(c2)
# 根据 c1 和 c2 的长度关系进行不同的处理
if lc1 < lc2:
# 如果 c1 的长度小于 c2，则返回 (c1[0]*0, c1)
return c1[:1]*0, c1
elif lc2 == 1:
# 如果 c2 的长度为 1，则返回 (c1/c2[-1], c1[0]*0)
return c1/c2[-1], c1[:1]*0
else:
# 将 c1 和 c2 转换为复数系数表示
z1 = _cseries_to_zseries(c1)
z2 = _cseries_to_zseries(c2)
# 对复数系数进行除法运算，得到商和余数
quo, rem = _zseries_div(z1, z2)
# 将得到的复数系数表示的商和余数转换为普通系数表示，并去除多余的零系数
quo = pu.trimseq(_zseries_to_cseries(quo))
rem = pu.trimseq(_zseries_to_cseries(rem))
# 返回计算结果：商和余数
return quo, rem
def chebpow(c, pow, maxpower=16):
"""Raise a Chebyshev series to a power.
Returns the Chebyshev series `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 ``T_0 + 2*T_1 + 3*T_2.``
Parameters
----------
c : array_like
1-D array of Chebyshev series coefficients ordered from low to
high.
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
Chebyshev series of power.
See Also
--------
chebadd, chebsub, chebmulx, chebmul, chebdiv
Examples
--------
>>> from numpy.polynomial import chebyshev as C
>>> C.chebpow([1, 2, 3, 4], 2)
array([15.5, 22. , 16. , ..., 12.5, 12. , 8. ])
"""
# note: this is more efficient than `pu._pow(chebmul, c1, c2)`, as it
# avoids converting between z and c series repeatedly
# c is a trimmed copy
[c] = pu.as_series([c]) # Convert `c` to a series if it's not already
power = int(pow) # Ensure `pow` is an integer
if power != pow or power < 0: # Check if `pow` is a non-negative integer
raise ValueError("Power must be a non-negative integer.")
elif maxpower is not None and power > maxpower: # Check if `pow` exceeds `maxpower`
raise ValueError("Power is too large")
elif power == 0: # Return 1 if power is 0
return np.array([1], dtype=c.dtype)
elif power == 1: # Return `c` itself if power is 1
return c
else:
# This can be made more efficient by using powers of two
# in the usual way.
zs = _cseries_to_zseries(c) # Convert Chebyshev series `c` to z-series
prd = zs
for i in range(2, power + 1):
prd = np.convolve(prd, zs) # Perform convolution to compute higher powers
return _zseries_to_cseries(prd) # Convert back to Chebyshev series and return
def chebder(c, m=1, scl=1, axis=0):
"""
Differentiate a Chebyshev series.
Returns the Chebyshev 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*T_0 + 2*T_1 + 3*T_2``
while [[1,2],[1,2]] represents ``1*T_0(x)*T_0(y) + 1*T_1(x)*T_0(y) +
2*T_0(x)*T_1(y) + 2*T_1(x)*T_1(y)`` if axis=0 is ``x`` and axis=1 is
``y``.
Parameters
----------
c : array_like
Array of Chebyshev 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 along which differentiation is performed. (Default: 0)
"""
# 将输入参数 `c` 转换为至少一维的 numpy 数组，并进行复制
c = np.array(c, ndmin=1, copy=True)
# 如果数组 `c` 的数据类型是布尔型、字节型、短整型、整型或长整型，则转换为双精度浮点型
if c.dtype.char in '?bBhHiIlLqQpP':
c = c.astype(np.double)
# 将参数 `m` 转换为整数，表示导数的阶数
cnt = pu._as_int(m, "the order of derivation")
# 将参数 `axis` 转换为整数，表示进行导数计算的轴
iaxis = pu._as_int(axis, "the axis")
# 如果导数的阶数 `cnt` 小于 0，则抛出 ValueError 异常
if cnt < 0:
raise ValueError("The order of derivation must be non-negative")
# 根据参数 `axis` 规范化轴的索引，确保在数组 `c` 的维度范围内
iaxis = normalize_axis_index(iaxis, c.ndim)
# 如果导数的阶数 `cnt` 等于 0，则直接返回输入数组 `c`
if cnt == 0:
return c
# 将数组 `c` 的轴 `iaxis` 移动到最前面，便于后续计算
c = np.moveaxis(c, iaxis, 0)
# 获取数组 `c` 的长度 `n`
n = len(c)
# 如果导数的阶数 `cnt` 大于等于 `n`，则将数组 `c` 的前部分清零
if cnt >= n:
c = c[:1] * 0
else:
# 否则，进行导数的计算
for i in range(cnt):
n = n - 1
c *= scl # 将数组 `c` 各元素乘以缩放因子 `scl`
# 创建空的导数数组 `der`，其形状与数组 `c` 的后续维度一致
der = np.empty((n,) + c.shape[1:], dtype=c.dtype)
for j in range(n, 2, -1):
# 计算 Chebyshev 多项式的导数
der[j - 1] = (2 * j) * c[j]
c[j - 2] += (j * c[j]) / (j - 2)
if n > 1:
der[1] = 4 * c[2]
der[0] = c[1]
c = der
# 将导数计算后的数组 `c` 恢复原来的轴顺序，与输入数组 `c` 一致
c = np.moveaxis(c, 0, iaxis)
# 返回计算得到的 Chebyshev 多项式的导数数组 `c`
return c
# 定义函数 chebint，用于对 Chebyshev 级数进行积分
def chebint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
"""
Integrate a Chebyshev series.
Returns the Chebyshev 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 ``T_0 + 2*T_1 + 3*T_2`` while [[1,2],[1,2]]
represents ``1*T_0(x)*T_0(y) + 1*T_1(x)*T_0(y) + 2*T_0(x)*T_1(y) +
2*T_1(x)*T_1(y)`` if axis=0 is ``x`` and axis=1 is ``y``.
Parameters
----------
c : array_like
Array of Chebyshev 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 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
C-series coefficients of the integral.
Raises
------
ValueError
If ``m < 1``, ``len(k) > m``, ``np.ndim(lbnd) != 0``, or
``np.ndim(scl) != 0``.
See Also
--------
chebder
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 chebyshev as C
>>> c = (1,2,3)
>>> C.chebint(c)
array([ 0.5, -0.5, 0.5, 0.5])
>>> C.chebint(c,3)
"""
# 将输入的参数 c 转换为至少有一维的 numpy 数组，并保留其副本
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")
# 如果积分的阶数 cnt 小于 0，抛出 ValueError 异常
if cnt < 0:
raise ValueError("The order of integration must be non-negative")
# 如果 k 的长度大于 cnt，抛出 ValueError 异常
if len(k) > cnt:
raise ValueError("Too many integration constants")
# 如果 lbnd 不是标量，抛出 ValueError 异常
if np.ndim(lbnd) != 0:
raise ValueError("lbnd must be a scalar.")
# 如果 scl 不是标量，抛出 ValueError 异常
if np.ndim(scl) != 0:
raise ValueError("scl must be a scalar.")
# 将操作的轴 iaxis 规范化，确保其在 c 的维度范围内
iaxis = normalize_axis_index(iaxis, c.ndim)
# 如果积分的阶数 cnt 等于 0，直接返回参数 c，无需积分计算
if cnt == 0:
return c
# 将 c 数组中的轴移动到索引位置 0，方便后续的积分操作
c = np.moveaxis(c, iaxis, 0)
# 将列表 k 补充到长度为 cnt，不足部分用 0 填充
k = list(k) + [0]*(cnt - len(k))
# 开始进行积分操作，循环 cnt 次
for i in range(cnt):
n = len(c) # 获取当前 c 数组的长度
c *= scl # 将 c 数组中的所有元素乘以标量 scl
# 如果 c 数组长度为 1，且其唯一元素为 0，则在该元素上加上 k[i]
if n == 1 and np.all(c[0] == 0):
c[0] += k[i]
else:
# 创建一个临时数组 tmp，用于存储积分过程中的中间结果
tmp = np.empty((n + 1,) + c.shape[1:], dtype=c.dtype)
tmp[0] = c[0]*0 # 第一个元素为 0
tmp[1] = c[0] # 第二个元素为 c[0]
# 如果 n 大于 1，则计算后续元素的值
if n > 1:
tmp[2] = c[1]/4
for j in range(2, n):
tmp[j + 1] = c[j]/(2*(j + 1))
tmp[j - 1] -= c[j]/(2*(j - 1))
# 计算 tmp[0]，加上 k[i] 减去 chebval(lbnd, tmp) 的值
tmp[0] += k[i] - chebval(lbnd, tmp)
c = tmp
# 完成积分操作后，将 c 数组中的轴移回到原始位置 iaxis
c = np.moveaxis(c, 0, iaxis)
# 返回最终的积分结果数组 c
return c
def chebval(x, c, tensor=True):
"""
Evaluate a Chebyshev series at points x.
If `c` is of length `n + 1`, this function returns the value:
.. math:: p(x) = c_0 * T_0(x) + c_1 * T_1(x) + ... + c_n * T_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
--------
chebval2d, chebgrid2d, chebval3d, chebgrid3d
Notes
-----
The evaluation uses Clenshaw recursion, aka synthetic division.
"""
# 将输入的系数数组转换成至少为1维的 numpy 数组，确保数据可用性和一致性
c = np.array(c, ndmin=1, copy=True)
# 如果系数数组的数据类型是布尔型或整型，则转换为双精度浮点型
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，则调整系数数组的形状以便广播计算
if isinstance(x, np.ndarray) and tensor:
c = c.reshape(c.shape + (1,)*x.ndim)
# 根据系数数组的长度选择不同的计算方法
if len(c) == 1:
c0 = c[0]
c1 = 0
elif len(c) == 2:
c0 = c[0]
c1 = c[1]
else:
x2 = 2*x
c0 = c[-2]
c1 = c[-1]
# 使用 Clenshaw 递归方法计算 Chebyshev 多项式的值
for i in range(3, len(c) + 1):
tmp = c0
c0 = c[-i] - c1
c1 = tmp + c1*x2
# 返回表达式 c0 + c1*x 的计算结果
return c0 + c1*x
# 在给定点(x, y)评估二维切比雪夫级数。
def chebval2d(x, y, c):
# 调用 pu 模块的 _valnd 函数，用于评估二维切比雪夫级数的值。
return pu._valnd(chebval, c, x, y)
# 在 x 和 y 的笛卡尔积上评估二维切比雪夫级数。
def chebgrid2d(x, y, c):
# 返回值为在网格点 (a, b) 上的二维切比雪夫级数值。
# 如果 c 的维度少于两个，隐式添加维度使其成为二维数组。
# 结果的形状将为 c.shape[2:] + x.shape + y.shape。
c : array_like
# 输入参数 c 是一个类数组结构，包含系数，按照多重度量 i,j 的系数存储在 `c[i,j]` 中。
# 如果 `c` 的维度高于两个，其余的索引用来枚举多组系数。
Returns
-------
values : ndarray, compatible object
# 返回值 `values` 是一个 ndarray 或兼容对象，包含二维切比雪夫级数在笛卡尔积 `x` 和 `y` 点上的值。
See Also
--------
chebval, chebval2d, chebval3d, chebgrid3d
# 参见其他相关函数：chebval, chebval2d, chebval3d, chebgrid3d。
Notes
-----
.. versionadded:: 1.7.0
# 添加于版本 1.7.0。
"""
return pu._gridnd(chebval, c, x, y)
def chebval3d(x, y, z, c):
return pu._valnd(chebval, c, x, y, z)
def chebgrid3d(x, y, z, c):
return pu._valnd(chebval, c, x, y, z)
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
--------
chebval, chebval2d, chebgrid2d, chebval3d
Notes
-----
.. versionadded:: 1.7.0
"""
使用`chebval`函数计算在给定参数下的多项式值，并通过`pu._gridnd`进行网格化处理。
def chebvander(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] = T_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 Chebyshev polynomial.
If `c` is a 1-D array of coefficients of length ``n + 1`` and `V` is the
matrix ``V = chebvander(x, n)``, then ``np.dot(V, c)`` and
``chebval(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 Chebyshev 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 Chebyshev polynomial. The dtype will be the same as
the converted `x`.
"""
# Convert deg to an integer if possible, raising a 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 type float64 or complex128
x = np.array(x, copy=None, ndmin=1) + 0.0
# Create dimensions for the result matrix v
dims = (ideg + 1,) + x.shape
dtyp = x.dtype
# Create an uninitialized array v with specified dimensions and type
v = np.empty(dims, dtype=dtyp)
# Use forward recursion to generate the entries of v
v[0] = x*0 + 1 # Initialize the first row of v
if ideg > 0:
x2 = 2*x
v[1] = x # Initialize the second row of v
# Use recursion to fill in the remaining rows of v
for i in range(2, ideg + 1):
v[i] = v[i-1]*x2 - v[i-2]
# Move the first axis to the last axis to match the expected output shape
return np.moveaxis(v, 0, -1)
def chebvander2d(x, y, deg):
"""Pseudo-Vandermonde matrix of given degrees.
Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
points ``(x, y)``. The pseudo-Vandermonde matrix is defined by
.. math:: V[..., (deg[1] + 1)*i + j] = T_i(x) * T_j(y),
where ``0 <= i <= deg[0]`` and ``0 <= j <= deg[1]``. The leading indices of
`V` index the points ``(x, y)`` and the last index encodes the degrees of
the Chebyshev polynomials.
If ``V = chebvander2d(x, y, [xdeg, ydeg])``, then the columns of `V`
correspond to the elements of a 2-D coefficient array `c` of shape
(xdeg + 1, ydeg + 1) in the order
.. math:: c_{00}, c_{01}, c_{02} ... , c_{10}, c_{11}, c_{12} ...
and ``np.dot(V, c.flat)`` and ``chebval2d(x, y, 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 2-D Chebyshev
series of the same degrees and sample points.
Parameters
----------
x : array_like
Array of points for the first dimension.
y : array_like
Array of points for the second dimension.
deg : list or tuple
List or tuple specifying the degrees of the resulting matrix along each dimension.
Returns
-------
vander : ndarray
The pseudo Vandermonde matrix. The shape of the returned matrix is
``x.shape + y.shape + (np.prod(deg) + 1,)``. The last index is the degree of the
corresponding 2-D Chebyshev polynomial. The dtype will be the same as
the dtype of the concatenated `x` and `y`.
"""
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
--------
chebvander, chebvander3d, chebval2d, chebval3d
Notes
-----
.. versionadded:: 1.7.0
"""
return pu._vander_nd_flat((chebvander, chebvander), (x, y), deg)
import pu
def chebvander3d(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] = T_i(x)*T_j(y)*T_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 Chebyshev polynomials.
If ``V = chebvander3d(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 ``chebval3d(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 Chebyshev
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
--------
chebvander, chebvander3d, chebval2d, chebval3d
Notes
-----
.. versionadded:: 1.7.0
"""
return pu._vander_nd_flat((chebvander, chebvander, chebvander), (x, y, z), deg)
def chebfit(x, y, deg, rcond=None, full=False, w=None):
"""
Least squares fit of Chebyshev series to data.
Return the coefficients of a Chebyshev 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 * T_1(x) + ... + c_n * T_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 Chebyshev series to be used for the fit.
rcond : float, optional
Cutoff for small singular values.
full : bool, optional
If True, return additional outputs.
w : array_like, optional
Weights applied to the y-coordinates.
Returns
-------
coef : ndarray
Coefficients of the Chebyshev series, with shape (deg + 1,) or
(deg + 1, K) if `y` is 2-D.
residuals, rank, s : ndarray, int, ndarray
Optional outputs as determined by `full`.
See Also
--------
chebvander, chebvander2d, chebvander3d, chebval, chebval2d, chebval3d
Notes
-----
This function is generally preferred over polyfit when the fitting
is done with respect to Chebyshev series due to its better numerical
properties.
"""
pass
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 -- 缩放后的范德蒙矩阵的数值秩
- singular_values -- 缩放后的范德蒙矩阵的奇异值
- rcond -- `rcond` 的值
详细信息请参见 `numpy.linalg.lstsq`。
Warns
-----
RankWarning
>>> import warnings
>>> warnings.simplefilter('ignore', np.exceptions.RankWarning)
See Also
--------
numpy.polynomial.polynomial.polyfit
numpy.polynomial.legendre.legfit
numpy.polynomial.laguerre.lagfit
numpy.polynomial.hermite.hermfit
numpy.polynomial.hermite_e.hermefit
chebval : 计算切比雪夫级数的值。
chebvander : 切比雪夫级数的范德蒙矩阵。
chebweight : 切比雪夫权重函数。
numpy.linalg.lstsq : 从矩阵计算最小二乘拟合。
scipy.interpolate.UnivariateSpline : 计算样条拟合。
Notes
-----
.. math:: E = \\sum_j w_j^2 * |y_j - p(x_j)|^2,
return pu._fit(chebvander, x, y, deg, rcond, full, w)
[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 = np.array([1.] + [np.sqrt(.5)]*(n-1))
top = mat.reshape(-1)[1::n+1]
bot = mat.reshape(-1)[n::n+1]
top[0] = np.sqrt(.5)
top[1:] = 1/2
bot[...] = top
mat[:, -1] -= (c[:-1]/c[-1])*(scl/scl[-1])*.5
return mat
r = la.eigvals(m)
r.sort()
return r
def chebinterpolate(func, deg, args=()):
"""Interpolate a function at the Chebyshev points of the first kind.
Returns the Chebyshev series that interpolates `func` at the Chebyshev
points of the first kind in the interval [-1, 1]. The interpolating
series tends to a minmax approximation to `func` with increasing `deg`
if the function is continuous in the interval.
.. versionadded:: 1.14.0
Parameters
----------
func : function
The function to be approximated. It must be a function of a single
variable of the form ``f(x, a, b, c...)``, where ``a, b, c...`` are
extra arguments passed in the `args` parameter.
deg : int
Degree of the interpolating polynomial
args : tuple, optional
Extra arguments to be used in the function call. Default is no extra
arguments.
Returns
-------
coef : ndarray, shape (deg + 1,)
Chebyshev coefficients of the interpolating series ordered from low to
high.
Examples
--------
>>> import numpy.polynomial.chebyshev as C
>>> C.chebinterpolate(lambda x: np.tanh(x) + 0.5, 8)
array([ 5.00000000e-01, 8.11675684e-01, -9.86864911e-17,
-5.42457905e-02, -2.71387850e-16, 4.51658839e-03,
2.46716228e-17, -3.79694221e-04, -3.26899002e-16])
Notes
-----
The Chebyshev polynomials used in the interpolation are orthogonal when
sampled at the Chebyshev points of the first kind. If it is desired to
constrain some of the coefficients they can simply be set to the desired
value after the interpolation, no new interpolation or fit is needed. This
is especially useful if it is known apriori that some of coefficients are
zero. For instance, if the function is even then the coefficients of the
terms of odd degree in the result can be set to zero.
"""
deg = np.asarray(deg)
if deg.ndim > 0 or deg.dtype.kind not in 'iu' or deg.size == 0:
raise TypeError("deg must be an int")
if deg < 0:
raise ValueError("expected deg >= 0")
order = deg + 1
xcheb = chebpts1(order)
yfunc = func(xcheb, *args)
m = chebvander(xcheb, deg)
c = np.dot(m.T, yfunc)
c[0] /= order
c[1:] /= 0.5 * order
return c
def chebgauss(deg):
"""
Gauss-Chebyshev quadrature.
Computes the sample points and weights for Gauss-Chebyshev quadrature.
These sample points and weights will correctly integrate polynomials of
degree :math:`2*deg - 1` or less over the interval :math:`[-1, 1]` with
the weight function :math:`f(x) = 1/\\sqrt{1 - x^2}`.
Parameters
----------
deg : int
Number of sample points and weights. It must be >= 1.
Returns
-------
x : ndarray
1-D ndarray containing the sample points.
y : ndarray
1-D ndarray containing the weights.
Notes
-----
.. versionadded:: 1.7.0
"""
"""
The results have only been tested up to degree 100, higher degrees may
be problematic. For Gauss-Chebyshev there are closed form solutions for
the sample points and weights. If n = `deg`, then
.. math:: x_i = \\cos(\\pi (2 i - 1) / (2 n))
.. math:: w_i = \\pi / n
"""
ideg = pu._as_int(deg, "deg")
if ideg <= 0:
raise ValueError("deg must be a positive integer")
x = np.cos(np.pi * np.arange(1, 2*ideg, 2) / (2.0*ideg))
w = np.ones(ideg)*(np.pi/ideg)
return x, w
def chebpts1(npts):
_npts = int(npts)
if _npts != npts:
raise ValueError("npts must be integer")
if _npts < 1:
raise ValueError("npts must be >= 1")
x = 0.5 * np.pi / _npts * np.arange(-_npts+1, _npts+1, 2)
return np.sin(x)
def chebpts2(npts):
_npts = int(npts)
if _npts != npts:
raise ValueError("npts must be integer")
if _npts < 2:
raise ValueError("npts must be >= 2")
x = np.linspace(-np.pi, 0, _npts)
return np.cos(x)
class Chebyshev(ABCPolyBase):
"""A Chebyshev series class.
The Chebyshev class provides the standard Python numerical methods
'+', '-', '*', '//', '%', 'divmod', '**', and '()' as well as the
attributes and methods listed below.
Parameters
----------
coef : array_like
Chebyshev coefficients in order of increasing degree, i.e.,
``(1, 2, 3)`` gives ``1*T_0(x) + 2*T_1(x) + 3*T_2(x)``.
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
"""
pass
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
"""
# 定义一些虚拟函数，这些函数用于处理 Chebyshev 多项式的各种运算
_add = staticmethod(chebadd)
_sub = staticmethod(chebsub)
_mul = staticmethod(chebmul)
_div = staticmethod(chebdiv)
_pow = staticmethod(chebpow)
_val = staticmethod(chebval)
_int = staticmethod(chebint)
_der = staticmethod(chebder)
_fit = staticmethod(chebfit)
_line = staticmethod(chebline)
_roots = staticmethod(chebroots)
_fromroots = staticmethod(chebfromroots)
@classmethod
def interpolate(cls, func, deg, domain=None, args=()):
"""Interpolate a function at the Chebyshev points of the first kind.
Returns the series that interpolates `func` at the Chebyshev points of
the first kind scaled and shifted to the `domain`. The resulting series
tends to a minmax approximation of `func` when the function is
continuous in the domain.
.. versionadded:: 1.14.0
Parameters
----------
func : function
The function to be interpolated. It must be a function of a single
variable of the form ``f(x, a, b, c...)``, where ``a, b, c...`` are
extra arguments passed in the `args` parameter.
deg : int
Degree of the interpolating polynomial.
domain : {None, [beg, end]}, optional
Domain over which `func` is interpolated. The default is None, in
which case the domain is [-1, 1].
args : tuple, optional
Extra arguments to be used in the function call. Default is no
extra arguments.
Returns
-------
polynomial : Chebyshev instance
Interpolating Chebyshev instance.
Notes
-----
See `numpy.polynomial.chebinterpolate` for more details.
"""
if domain is None:
domain = cls.domain
# 创建一个函数 xfunc，将 func 在指定 domain 上的映射和参数 args 一起应用
xfunc = lambda x: func(pu.mapdomain(x, cls.window, domain), *args)
# 调用 chebinterpolate 函数生成系数 coef
coef = chebinterpolate(xfunc, deg)
# 返回一个 Chebyshev 实例，以 coef 作为系数，domain 作为定义域
return cls(coef, domain=domain)
# 定义一些虚拟属性
domain = np.array(chebdomain)
window = np.array(chebdomain)
basis_name = 'T'
