python 计算给定的场区拐点坐标的中心点坐标

目录

• SciPy简介
• SciPy基础

• Numpy基础简要回顾
• Numpy的stack，vstack，hstack，dstack，concatenate
• SciPy特殊函数
• SciPy常量

• SciPy应用

• K-means
• FFT快速傅里叶变换
• 积分计算
• 插值
• Scipy I/O
• 线性代数
• 数字图像处理
• optimize
• 信号处理
• 统计

SciPy简介

SciPy，是一个python开源库，在BSD授权下发布，主要用于数学、科学和工程计算。SciPy库依赖于NumPy，NumPy提供了方便和快速的n维数组操作。它们一起可以运行在所有流行的操作系统上，安装简单，使用免费。现在，组合使用NumPy、SciPy和Matplotlib，作为MATLAB的替代品已经成为趋势。相比MATLAB，Python功能更强大、编程更容易；

SciPy使用的基本数据结构是NumPy模块提供的多维数组ndarray。NumPy提供了用于线性代数、傅里叶变换和随机数生成的函数，SciPy中也提供了同样的功能，并且通用性更强；另外，SciPy的安装简便：

conda install scipy

SciPy基础

Numpy基础简要回顾

默认情况下，所有NumPy函数都可以在SciPy中使用。当导入SciPy时，不需要显式地导入NumPy函数。NumPy的主要对象是n次多维数组ndarray，SciPy构建在ndarray数组之上，ndarray是存储单一数据类型的多维数组。在NumPy中，维度称为轴，轴的数量称为秩；

线性代数主要处理矩阵运算，现在首先回顾NumPy中数组的基本功能；ndarray是NumPy中最重要的类。标准的Python列表(list)中，元素是对象。如：L = [1, 2, 3]，实际需要3个指针和三个整数对象，对于数值运算比较浪费资源。与此不同，ndarray中元素直接存储为原始数据，元素的类型由ndarray对象中的属性dtype描述。当ndarray数组中的元素，通过索引或切片返回时，会根据dtype，从原始数据转换成Python对象，以便外部使用；比如，将Python类数组对象转换为NumPy数组：

import numpy as np
list = [1,2,3,4]
arr = np.array(list)
print (arr) # [1 2 3 4]
print (type(arr)) # <class 'numpy.ndarray'>

在NumPy中，也可以使用以下函数创建ndarray数组：

• 1.zeros，zeros()函数创建数组，并且把数组元素的值初始化为0，可以指定数组形状和数据类型：

import numpy as np
print (np.zeros((2, 3)))
# [[0. 0. 0.]
#  [0. 0. 0.]]

• 2.ones，ones()函数创建数组，并且把数组元素的值初始化为1，可以指定数组形状和数据类型：

import numpy as np
print (np.ones((2, 3)))
# [[1. 1. 1.]
#  [1. 1. 1.]]

• 3.arange，arange()函数创建递增数组：

import numpy as np
print (np.arange(7)) # [0 1 2 3 4 5 6]

• 4.linspace，linspace()函数创建一个数组，该数组包含指定区间内均匀分布的值：

"""
np.linspace(
    start,
    stop,
    num=50,
    dtype=None)
"""
import numpy as np
print (np.linspace(1., 4., 6)) # [1.  1.6 2.2 2.8 3.4 4. ]

关于数组的数据类型，数据类型对象dtype，是描述数组中元素数据类型的对象：

import numpy as np
arr = np.arange(2, 10, dtype = np.float)
print (arr) # [2. 3. 4. 5. 6. 7. 8. 9.]
print ("数组数据类型 :", arr.dtype) # 数组数据类型 : float64

矩阵是一种特殊的二维数组，它有一些特殊的运算符，如*(矩阵乘法)和**(矩阵幂)，首先，可以用np.matrix创建一个矩阵（数据结构从ndarray变成matrix）：

import numpy as np
print (np.matrix('1 2; 3 4'))
# [[1 2]
#  [3 4]]

print(type(np.matrix('1 2; 3 4')))
# <class 'numpy.matrix'>

# 通过numpy.matrix将2维ndarray转为matrix
a  = np.array([[1, 2], [3, 4]])
A = np.matrix(a)
print(type(A)) # <class 'numpy.matrix'>

矩阵转置：

mat = np.matrix('1 2; 3 4')
mat.T
# matrix([[1, 3],
#         [2, 4]])

共轭就是矩阵每个元素都取共轭（复数的实部不变，虚部取负），共轭转置就是先取共轭，再取转置：

mat = np.matrix('1 2; 3 4')
print (mat.H)
# [[1 3]
#  [2 4]]

单位矩阵通常在矩阵的乘法中有特殊作用。单位矩阵是个方阵，从左上角到右下角的对角线（称为主对角线）上的元素均为1，除此以外全都为0：

import numpy.matlib
print (numpy.matlib.identity(5))
"""
[[1. 0. 0. 0. 0.]
 [0. 1. 0. 0. 0.]
 [0. 0. 1. 0. 0.]
 [0. 0. 0. 1. 0.]
 [0. 0. 0. 0. 1.]]
"""

逆矩阵的数学定义：存在矩阵$M$以及矩阵$N$，假如$MN=NM = I$（Identify Matrix单位矩阵），那么矩阵$M$和矩阵$N$互为逆矩阵；求一个矩阵的逆矩阵：

"""matrix的计算"""
import numpy as np
mat = np.matrix('1 2; 3 4')
mat2 = mat.I
print(mat2)
# [[-2.   1. ]
#  [ 1.5 -0.5]]

"""ndarray的计算"""
import numpy as np

a  = np.array([[1, 2], [3, 4]])  # 初始化一个非奇异矩阵(数组)
print(np.linalg.inv(a))  # 求逆矩阵
# [[-2.   1. ]
#  [ 1.5 -0.5]]

# matrix对象可以通过 .I 方便地求逆
A = np.matrix(a)
print(A.I)

Numpy的stack，vstack，hstack，dstack，concatenate

stack函数用于，沿着新的轴线堆叠一系列的数组：

np.stack(arrays, axis=0)

• arrays：待加入的数组序列，每个数组必须具有相同的形状；
• axis：堆叠输入数组后，结果数组中的新增轴；

实例如下：

>>> arrays = [np.random.randn(3, 4) for _ in range(10)]
>>> np.stack(arrays, axis=0).shape
(10, 3, 4)

>>> np.stack(arrays, axis=1).shape
(3, 10, 4)

>>> np.stack(arrays, axis=2).shape
(3, 4, 10)

>>> a = np.array([1, 2, 3]) # shape(3,)
>>> b = np.array([2, 3, 4]) # shape(3,)
>>> np.stack((a, b)) # shape(2,3)
array([[1, 2, 3],
       [2, 3, 4]])

>>> np.stack((a, b), axis=-1) # shape(3,2)
array([[1, 2],
       [2, 3],
       [3, 4]])

vstack函数，沿着竖直方向（沿着axis=0轴）堆叠数组：

np.vstack(tuple)

• tuple：n维数组序列，用元组对象表示；各个数组必须在除第一轴外的其他轴具有相同形状；特殊的，一维数组必须具有相同的长度；

实例如下：

# 堆叠向量(一维数组)会新增轴
>>> a = np.array([1, 2, 3]) # shape(3,)一维数组,即向量
>>> b = np.array([2, 3, 4]) # shape(3,)
>>> np.vstack((a,b)) # shape(2,3)
array([[1, 2, 3],
       [2, 3, 4]])

# 堆叠高维数组不会新增轴
>>> a = np.array([[1], [2], [3]]) # shape(3,1)二维数组
>>> b = np.array([[2], [3], [4]]) # shape(3,1)
>>> np.vstack((a,b)) # shape(6,1)
array([[1],
       [2],
       [3],
       [2],
       [3],
       [4]])

hstack函数，沿着横向（axis=1轴方向）堆叠数组：

np.hstack(tup)

• tup：n维数组序列；一维数组可以是任意长度；对于高维数组，各个数组沿除第二轴以外的所有方向必须具有相同的形状；

实例如下：

>>> a = np.array((1,2,3)) # shape(3,)
>>> b = np.array((2,3,4)) # shape(3,)
>>> np.hstack((a,b)) # shape(6,)
array([1, 2, 3, 2, 3, 4])

>>> a = np.array([[1],[2],[3]]) # shape(3,1)
>>> b = np.array([[2],[3],[4]]) # shape(3,1)
>>> np.hstack((a,b)) # shape(3,2)
array([[1, 2],
       [2, 3],
       [3, 4]])

dstack函数，沿着axis=2轴方向堆叠数组：

np.dstack(tup)

• tup：n维数组序列；各个数组必须沿除第三轴以外的所有方向具有相同的形状；特别的，一维或二维数组必须具有相同的形状；

实例如下：

>>> a = np.array((1,2,3)) # shape(3,)
>>> b = np.array((2,3,4)) # shape(3,)
# 注意堆叠一维数组时会在axis=0处新增轴,使之变成三维
>>> np.dstack((a,b)) # shape(1,3,2)
array([[[1, 2],
        [2, 3],
        [3, 4]]])

>>> a = np.array([[1],[2],[3]]) # shape(3,1)
>>> b = np.array([[2],[3],[4]]) # shape(3,1)
>>> np.dstack((a,b)) # shape(3,1,2)
array([[[1, 2]],
       [[2, 3]],
       [[3, 4]]])

concatenate函数用于沿着 已经存在的轴 拼接数组：

concatenate((a1, a2, ...), axis=0)

• (a1, a2, ...)：数组序列，各个数组需要具有相同的形状，"操作轴"对应的位置除外；
• axis：指定的操作轴；如果设置axis=None，则数组在使用前会被压扁，再拼接成一个向量

实例如下：

>>> a = np.array([[1, 2], [3, 4]]) # shape(2,2)
>>> b = np.array([[5, 6]]) # shape(1,2)
>>> np.concatenate((a, b), axis=0) # shape(3,2)
array([[1, 2],
       [3, 4],
       [5, 6]])
>>> np.concatenate((a, b.T), axis=1) # CAT[shape(2,2),shape(2,1)]->shape(2,3)
array([[1, 2, 5],
       [3, 4, 6]])
>>> np.concatenate((a, b), axis=None) # CAT[shape(4,),shape(2,)]->shape(6,)
array([1, 2, 3, 4, 5, 6])

SciPy特殊函数

scipy.special 模块中包含了一些常用的杂项函数，经常使用的一般有：

• 立方根函数
• 指数函数
• 相对误差指数函数
• 对数和指数函数
• 兰伯特函数
• 排列组合函数
• $\gamma$函数

比如求立方根函数：

from scipy.special import cbrt
res = cbrt([1000, 27, 8, 23])
print (res)
# [10.          3.          2.          2.84386698]

SciPy常量

SciPy的常量包提供了很多科学领域的常量，例如：光速，要使用常量，需要先导入常量包（scipy.constants）；比如，从scipy.constants中导入$\pi$值：

#导入pi常量
from scipy.constants import pi

print("sciPy : pi = %.16f"%pi) # sciPy : pi = 3.1415926535897931

下表列出了常用的物理常数：

常量比较多，不可能都记住，可以使用scipy.constants.find()函数查找常量。scipy.constants.find()函数返回physical_constants常量字典的键列表，如果关键字不匹配，则不返回任何内容。获得键列表后，可以使用physical_constants['key']获取常量：

import scipy
from scipy.constants import find, physical_constants
res = scipy.constants.find("boltzmann")
print (res)
"""
['Boltzmann constant', 'Boltzmann constant in Hz/K', 'Boltzmann constant in eV/K', 
'Boltzmann constant in inverse meters per kelvin', 'Stefan-Boltzmann constant']
"""
print(scipy.constants.physical_constants['Boltzmann constant'])
# (1.38064852e-23, 'J K^-1', 7.9e-30)

SciPy应用

K-means

K均值聚类(K-means clustering)是在一组未标记的数据中，将相似的数据归到同一类别。聚类与分类的最大不同在于分类的目标事先已知，而聚类则不知道。K-means是聚类中最常用的方法之一，它是基于样本与样本的距离来计算最佳类别归属，即靠得比较近的一组数据被归为一类。K-means的算法原理如下：

• 1.随机选取k个点作为中心点；
• 2.遍历所有点，将每个点划分到最近的中心点，形成k个聚类；
• 3.根据聚类中点之间的距离，重新计算各个聚类的中心点；
• 4.重复2-3步骤，直到这k个中心点不再变化（收敛），或达到最大迭代次数；

SciPy中，cluster已经实现了K-Means，可以直接使用它；导入要使用的包和模块：

from scipy.cluster.vq import kmeans,vq,whiten
import numpy as np

准备样本：

from numpy import vstack,array
from numpy.random import rand

# 具有3个特征值的样本数据生成
data = vstack((rand(100,3) + array([.5,.5,.5]),rand(100,3))) # shape(200,3)

---

补充：np.random.rand，创建给定形状的数组，并用均匀分布的随机样本填充它：

rand(d0, d1, ..., dn)

• d0, d1, ..., dn : "int"，用于设置数组形状；

---

数据白化预处理是一种常见的数据预处理方法，作用是去除样本数据的冗余信息，数据白化回顾 算法栈-PCA主成分分析，白化过程为：

# 白化数据
data = whiten(data)

把样本数据分成2个聚类，使用kmeans函数计算2个聚类的中心点：

# 计算 K = 2 时的中心点
centroids, _ = kmeans(data, 2)

打印中心点：

print(centroids)
"""
[[2.61307329 2.53808333 2.54104294]
 [1.23269589 1.14943359 1.25853145]]
"""

使用vq函数将样本数据中的每个样本点分配给一个中心点，形成3个聚类：

# 将样本数据中的每个值分配给一个中心点，形成3个聚类。
# 返回值clx标出了对应索引样本的聚类，dist表示对应索引样本与聚类中心的距离。
clx, dist = vq(data, centroids)

# 打印聚类
print(clx)
"""
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1
 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 1 1 1 0 1 1 1 1 1 1 1
 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1]
"""

FFT快速傅里叶变换

计算机只能处理离散信号，使用离散傅里叶变换(DFT) 是计算机分析信号的基本方法。但是离散傅里叶变换的缺点是：计算量大，时间复杂度太高，当采样点数太高的时候，计算缓慢，由此出现了DFT的快速实现，即快速傅里叶变换FFT（回顾：算法栈-其他算法FFT）；

快速傅里叶变换（FFT）是计算量更小的离散傅里叶变换，其逆变换被称为快速傅里叶逆变换（IFFT）；使用SciPy的fftpack可以直接对数据进行FFT和IFFT，实例如下：

import numpy as np

#从fftpack中导入fft(快速傅里叶变化)和ifft(快速傅里叶逆变换)函数
from scipy.fftpack import fft,ifft

#创建一个数组
x = np.array([1.0, 2.0, 1.0, -1.0, 1.5])

#对数组数据进行傅里叶变换
y = fft(x)
print('fft:',y)

#快速傅里叶逆变换
yinv = ifft(y)
print('ifft:',yinv)

"""
fft: [ 4.5       +0.j          2.08155948-1.65109876j -1.83155948+1.60822041j
 -1.83155948-1.60822041j  2.08155948+1.65109876j]
ifft: [ 1. +0.j  2. +0.j  1. +0.j -1. +0.j  1.5+0.j]
"""

积分计算

Scipy中的integrate提供了很多数值积分方法，例如，一重积分、二重积分、三重积分、多重积分、高斯积分等等；

对于一重积分，quad函数用于求一重积分。例如，在给定的$a$到$b$范围内，对函数$f(x)$求一重积分：
$\int_{a}^{b}f(x)dx$
quad的一般形式是scipy.integrate.quad(f,a,b)，其中f是求积分的函数名，a和b分别是下限和上限，实例如下，对函数$f(x)=exp(-x^{2})$计算$(0,5)$的积分：

import scipy.integrate
from numpy import exp

f = lambda x:exp(-x**2)
i = scipy.integrate.quad(f, 0, 5)
print(i)
# (0.8862269254513955, 2.3183115159980698e-14)

quad返回两个值，第一个值是积分的值，第二个值是对积分值的绝对误差估计；

如果积分的函数$f$带系数参数，即：
$I(a,b)=\int_{0}^{1}(ax^{2}+b)dx$
对于这种情况，系数$a$和$b$可以使用args传入quad：

from scipy.integrate import quad

# 注意被积分变量应设为第一个位置参数
def f(x,a,b):
    return a * (x ** 2) + b

ret = quad(f, 0, 1, args=(3, 1))
print (ret)
# (2.0, 2.220446049250313e-14)

对于重积分，要计算二重积分、三重积分、多重积分，可使用dblquad、tplquad和nquad函数，以二重积分为例，假设有二重积分：
$\int_{0}^{\frac{1}{2}}dy\int_{0}^{\sqrt{1-4y^{2}}}19xydx$
dblquad的一般形式是scipy.integrate.dblquad(func, a, b, gfun, hfun)，其中，func是待积分函数的名称，a、b是$x$变量的上下限，gfun、hfun为定义$y$变量上下限的函数名称；

使用lambda表达式定义函数func、gfun和hfun。注意，在很多情况下gfun和hfun可能是常数，但是即使gfun和hfun是常数，也必须被定义为函数，计算如下：

import scipy.integrate
from numpy import exp
from math import sqrt

func = lambda x, y : 19*x*y
gfun = lambda x : 0
hfun = lambda y : sqrt(1-4*y**2)

i = scipy.integrate.dblquad(func, 0, 0.5, gfun, hfun)
print (i)
# (0.59375, 2.029716563995638e-14)

插值

插值，是依据一系列的点$(x_{i},y_{i})$，通过一定的算法找到一个合适的函数来包含(逼近)这些点，反应出这些点的走势规律，然后根据走势规律求其他点值的过程。scipy.interpolate包里有很多类可以实现对一些已知的点进行插值，即找到一个合适的函数，例如，interp1d类，当得到插值函数后便可用这个插值函数计算其他$x_{j}$对应的的$y_{j}$值。

一维插值interp1d实例如下，首先生成数据，通过采样几个点获取数据：

import numpy as np
from scipy import interpolate as intp
import matplotlib.pyplot as plt

x = np.linspace(0, 4, 12)
y = np.cos(x**2/3 + 4)

print(x)
print(y)
"""
[0.         0.36363636 0.72727273 1.09090909 1.45454545 1.81818182
 2.18181818 2.54545455 2.90909091 3.27272727 3.63636364 4.        ]
 
[-0.65364362 -0.61966189 -0.51077021 -0.31047698 -0.00715476  0.37976236
  0.76715099  0.99239518  0.85886263  0.27994201 -0.52586509 -0.99582185]
"""

# 可视化数据
plt.plot(x, y,'o')
plt.show()

根据上面的数据，使用interp1d类创建拟合函数，创建两个函数f1和f2。通过这些函数，输入x可以计算y。kind表示插值使用的技术类型，例如：’Linear’, ‘Nearest’, ‘Zero’, ‘Slinear’, ‘Quadratic’, ‘Cubic’等：

f1 = intp.interp1d(x, y, kind = 'linear')

f2 = intp.interp1d(x, y, kind = 'cubic')

增加输入数据，对拟合的函数进行比较：

xnew = np.linspace(0, 4, 30)

plt.plot(x, y, 'o', xnew, f1(xnew), '-', xnew, f2(xnew), '--')
plt.legend(['data', 'linear', 'cubic'])
plt.show()

另外，可以通过interpolate模块中UnivariateSpline类对含有噪声的数据进行插值运算。使用UnivariateSpline类，输入一组数据点，通过绘制一条平滑曲线来去除噪声。绘制曲线时可以设置平滑参数s，如果参数s=0，将对所有点(包括噪声)进行插值运算，也就是说s=0时不去除噪声，实例如下：

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import UnivariateSpline

x = np.linspace(-3, 3, 50)
y = np.exp(-x**2) + 0.1 * np.random.randn(50) # 通过random方法添加噪声数据
plt.plot(x, y, 'ro')

# 平滑参数使用默认值
spl = UnivariateSpline(x, y)
xs = np.linspace(-3, 3, 1000)
plt.plot(xs, spl(xs), 'b') # 蓝色曲线

# 设置平滑参数
spl.set_smoothing_factor(0.5)
plt.plot(xs, spl(xs), 'g') # 绿色曲线

# 设置平滑参数为0
spl.set_smoothing_factor(0)
plt.plot(xs, spl(xs), 'yellow') # 黄色曲线

plt.legend(['data','default','smooth0.5','without smooth'])

plt.show()

Scipy I/O

scipy.io（输入和输出）用于读写各种格式的文件。scipy.io支持的格式很多：

• Matlab
• IDL
• Matrix Market
• Wave
• Arff
• Netcdf

其中，Matlab 格式是最常用的，以下是用于加载和保存.mat文件的函数：

• loadmat ：加载MATLAB文件
• savemat ：保存为MATLAB文件
• whosmat ：列出MATLAB文件中的变量

实例如下：

import scipy.io as sio
import numpy as np

# 保存mat文件
vect = np.arange(20)
sio.savemat('array.mat', {'vect':vect})

# 加载文件
mat_file_content = sio.loadmat('array.mat')
print (mat_file_content)
"""
{'__header__': b'MATLAB 5.0 MAT-file Platform: posix, Created on: Tue Mar 23 20:33:28 2021', 
 'vect': array([[ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15,
        16, 17, 18, 19]]), 
 '__globals__': [], '__version__': '1.0'}
"""

如果想在不加载文件的情况下，查看文件中的概要内容，可使用whosmat命令，如下所示：

import scipy.io as sio

mat_file_content = sio.whosmat('array.mat')
print (mat_file_content)
# [('vect', (1, 20), 'int64')]

线性代数

SciPy.linalg与NumPy.linalg相比，scipy.linalg除了包含numpy.linalg中的所有函数，还具有numpy.linalg中没有的高级功能。scipy.linalg.solve 函数可用于解线性方程。例如，对于线性方程$ax+by=z$，求出未知数$x,y$值；

比如，解下面的联立方程组：

上面的方程组，可以用矩阵表示为：

下面使用scipy求解；scipy.linalg.solve函数接受两个输入，数组a和数组b，数组a表示系数，数组b表示等号右侧值，求出的解将会放在一个数组里返回：

from scipy import linalg
import numpy as np

# 声明numpy数组
a = np.array([[1, 3, 5], [2, 5, 1], [2, 3, 8]])
b = np.array([10, 8, 3])

# 求解
x = linalg.solve(a, b)

# 输出解值
print(x)
# [-9.28  5.16  0.76]

对于行列式的计算，可以使用det()函数计算行列式，它接受一个矩阵作为输入，返回一个标量值，即该矩阵的行列式值：

from scipy import linalg
import numpy as np

# 声明numpy数组
A = np.array([[3,4],[7,8]])

# 计算行列式
x = linalg.det(A)

# 输出结果
print (x)
# -4.0

求矩阵的特征值、特征向量，是线性代数中的常见计算；通常，可以根据以下关系，求取矩阵$(A)$的特征值$(\lambda)$、特征向量$(v)$：
$Av=\lambda v$
scipy.linalg.eig 函数可用于计算特征值与特征向量，函数返回特征值和特征向量：

from scipy import linalg
import numpy as np

# 声明numpy数组
A = np.array([[3,4],[7,8]])

# 求解
l, v = linalg.eig(A)

# 打印特征值
print('特征值')
print (l)

# 打印特征向量
print('特征向量')
print (v)

"""
特征值
[-0.35234996+0.j 11.35234996+0.j]
特征向量
[[-0.76642628 -0.43192981]
 [ 0.64233228 -0.90190722]]
"""

对于奇异值分解，（回顾：算法栈-其他算法SVD），实例如下：

from scipy import linalg
import numpy as np

# 声明numpy数组
a = np.random.randn(3, 2) + 1.j*np.random.randn(3, 2)

# 输出原矩阵
print('原矩阵')
print(a)

# 求解
U, s, Vh = linalg.svd(a)

# 输出结果
print('奇异值分解')
print(U, "U")
print(Vh, "Vh")
print(s, "s")
"""
原矩阵
[[-0.06484513-0.64135452j -0.32497456+1.26679101j]
 [-0.03940525-0.76532795j -1.07591799-0.52985313j]
 [-0.48470857+0.86705556j  0.88468158-0.82965768j]]
奇异值分解
[[-0.25109325-0.54196933j -0.28247266+0.01554266j  0.3834136 -0.64512252j]
 [ 0.40752311-0.28991777j  0.63883107+0.53486208j -0.07206305-0.22471522j]
 [-0.1166846 +0.61601845j  0.43559576-0.18984564j  0.48906829-0.37674016j]] U
[[ 0.48305327+0.j         -0.62614967-0.61204258j]
 [-0.87559096+0.j         -0.34543944-0.33765672j]] Vh
[2.4021792  0.91585379] s
"""

数字图像处理

scipy.ndimage提供了图像矩阵变换、图像滤波、图像卷积等功能；加载原图像：

%matplotlib inline
from scipy import ndimage
import matplotlib.image as mpimg
import matplotlib.pyplot as plt

# 加载图片
face = mpimg.imread('pi.png')

# 显示图片
plt.imshow(face)
plt.show()

旋转图片，可以使用ndimage.rotate函数：

# 旋转图片45度
rotate_face = ndimage.rotate(face, 45)

plt.imshow(rotate_face)
plt.show()

图像滤波是一种图像增强技术。例如，可以图像滤波突出图像的某些特性，弱化或删除图像的另一些特性。滤波有很多种，例如：平滑、锐化、边缘增强等等；下面对图像进行高斯滤波。高斯滤波是一种模糊滤波，广泛用于滤除图像噪声：

# 处理图片
# sigma=3表示模糊程度为3
face1 = ndimage.gaussian_filter(face, sigma=3)

# 显示图片
plt.imshow(face1)
plt.show()

边缘检测是寻找图像中物体边界的技术。它的原理是通过检测图像中的亮度突变，来识别物体边缘。边缘检测在图像处理、计算机视觉等领域中广泛应用。常用边缘检测算法包括：

• Sobel
• Canny
• Prewitt
• Roberts
• Fuzzy Logic methods

考虑下面的例子：

import scipy.ndimage as nd
import numpy as np

# 生成图像
im = np.zeros((256, 256))
im[64:-64, 64:-64] = 1
im[90:-90,90:-90] = 2
# 高斯滤波
im = nd.gaussian_filter(im, 8)

import matplotlib.pyplot as plt
plt.imshow(im)
plt.show()

使用ndimage的Sobel函数来检测图像边缘，该函数会对图像数组的每个轴分开操作，因此上面的二维数组会产生两个矩阵，然后使用NumPy的Hypot函数将这两个矩阵合并为一个矩阵，得到最后结果：

sx = nd.sobel(im, axis = 0, mode = 'constant') #shape(256,256)
sy = nd.sobel(im, axis = 1, mode = 'constant') #shape(256,256)
sob = np.hypot(sx, sy) #shape(256,256)

plt.imshow(sob)
plt.show()

---

补充：Hypot函数

numpy.hypot(x1,x2)

• x1,x2：两个形状相同的数组

hypot对两数组$x_{1},x_{2}$的元素逐个进行操作，若输出为$x_{3}$，则有：
$x_{3}(i,j)=\sqrt{(x_{1}(i,j))^{2}+(x_{2}(i,j))^{2}}$

---

optimize

SciPy的optimize提供了常用数值优化算法，一些经典的优化算法包括线性回归、函数极值和根的求解以及确定两函数交点的坐标等；

对于标量函数极值求解，比如计算$f(x)=x^{2}+2x+9$的极小值，先可视化该函数：

%matplotlib inline
import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt

# 定义函数
def f(x):
    return x**2 + 2*x + 9

# x取值：-10到10之间，间隔0.1
x = np.arange(-10, 10, 0.1)

# 画出函数曲线
plt.plot(x, f(x))
plt.show()

计算该函数最小值的有效方法之一是使用带起点的BFGS算法（拟牛顿法）。该算法从参数给定的起始点计算函数的梯度下降，并输出梯度为零、二阶导数为正的极小值。BFGS算法是由Broyden，Fletcher，Goldfarb，Shanno四个人分别提出的，故称为BFGS；使用BFGS函数，找出抛物线函数的最小值：

# 第一个参数是函数名，第二个参数是梯度下降的起点。返回值是函数最小值的x值(ndarray数组)
xopt = optimize.fmin_bfgs(f,0)

xmin = xopt[0] # x值
ymin = f(xmin) # y值，即函数最小值
print('xmin: ', xmin)
print('ymin: ', ymin)

# 原函数曲线
plt.plot(x, f(x))
# 画出最小值的点, s=20设置点的大小，c='r'设置点的颜色
plt.scatter(xmin, ymin, s=20, c='r')

plt.show()

fmin_bfgs存在缺点，当函数有局部最小值，该算法会因起始点不同，找到这些局部最小而不是全局最小。对于这种情况，可以使用scipy.optimize提供的basinhopping()方法，该方法把局部优化方法与起始点随机抽样相结合，求出全局最小值，代价是耗费更多计算资源：

import numpy as np
from scipy import optimize

# 定义函数
def g(x):
    return x**2 + 20*np.sin(x)

# 求取最小值，初始值为7
ret = optimize.basinhopping(g, 7)

print(ret)
"""
                        fun: 0.15825752683178962
 lowest_optimization_result:       fun: 0.15825752683178962
 hess_inv: array([[0.0497854]])
      jac: array([2.38418579e-07])
  message: 'Optimization terminated successfully.'
     nfev: 18
      nit: 4
     njev: 6
   status: 0
  success: True
        x: array([4.27109533])
                    message: ['requested number of basinhopping iterations completed successfully']
      minimization_failures: 0
                       nfev: 1647
                        nit: 100
                       njev: 549
                          x: array([4.27109533])
"""

可以看到全局最小值被正确找出：fun: 0.15825752683178962，x: array([4.27109533])；

要求取一定范围之内的函数最小值，可使用fminbound方法：

optimize.fminbound(func, x1, x2)

• func ：目标函数；
• x1, x2 ：范围边界；

实例如下：

import numpy as np
from scipy import optimize

# 定义函数
def g(x):
    return x**2 + 20*np.sin(x)

# 求取-10到-5之间的函数最小值。full_output=True表示返回详细信息。
ret = optimize.fminbound(g, -10, -5, full_output=True)

print(ret)
# (-7.068891380019064, 35.82273589215206, 0, 12)
# 函数最小值是：35.82273589215206，对应的x值：-7.068891380019064

对于函数求解，有函数$f(x)$，当$f(x) = 0$，求$x$的值，即为函数求解，可以使用fsolve()函数；解单个方程：

"""
optimize.fsolve(func， x0)
    func 目标函数
    x0 初始值
"""
import numpy as np
from scipy import optimize

# 定义函数
def g(x):
    return x**2 + 20*np.sin(x)

# 函数求解
ret = optimize.fsolve(g, 2)

print(ret)
# [0.] 解出的根值是：0

对于方程组求解，可以先定义方程组：

def func(x):
	# unpack解包
    u1,u2,u3 = x
    return [f1(u1,u2,u3), f2(u1,u2,u3), f3(u1,u2,u3)]

方程组求解举例：

from scipy.optimize import fsolve
from math import sin,cos

def f(x):
    x0 = float(x[0])
    x1 = float(x[1])
    x2 = float(x[2])
    return [
        4*x1 + 9,
        3*x0*x0 - sin(x1*x2),
        x1*x2 - 2.5
    ]

result = fsolve(f, [1,1,1])
print (result)
# [-0.44664383 -2.25       -1.11111111]

假设有一批数据样本，要创建这些样本数据的拟合曲线或函数，可以使用Scipy.optimize的curve_fit()函数：

optimize.curve_fit(func, x1, y1)
    func:目标函数
    x1,y1:样本数据

使用以下函数演示曲线拟合：
$f(x)=50cos(x)+2$

import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt

# 函数模型用于生成数据
def g(x, a, b):
    return a*np.cos(x) + b

# 产生含噪声的样本数据
x_data = np.linspace(-5, 5, 100) # 采样点
y_data = g(x_data, 50, 2) + 5*np.random.randn(x_data.size) # 加入随机数作为噪声

# 使用curve_fit()函数来估计a和b的值
variables, variables_covariance = optimize.curve_fit(g, x_data, y_data)

# 输出结果
print('求出的系数a, b: ')
print(variables)

print('variables_covariance: ')
print(variables_covariance)

"""
求出的系数a, b: 
[50.4892627   1.95034402]
variables_covariance: 
[[0.62298888 0.11641726]
 [0.11641726 0.29216138]]
 
variables是给定模型的最优参数，variables_covariance可用于检查拟合情况，
其对角线元素值表示每个参数的方差，可看出正确求出了系数值
"""

# 可视化结果
y = g(x_data, variables[0], variables[1])

plt.plot(x_data, y_data, 'o', color="green", label = "Samples")
plt.plot(x_data, y, color="red", label = "Fit")
plt.legend()

plt.show()

另外，最小二乘法是非常经典的数值优化算法，通过最小化误差的平方和来寻找最符合数据的曲线。optimize提供了实现最小二乘拟合算法的函数leastsq()，leastsq是least square的简写，即最小二乘法：

optimize.leastsq(func, x0, args=())
    func:计算误差的函数
    x0:是计算的初始参数值
    args:是指定func的其他参数

举例：

import numpy as np
from scipy import optimize    

# 样本数据
X = np.array([160,165,158,172,159,176,160,162,171])
Y = np.array([58,63,57,65,62,66,58,59,62])

# 偏差函数, 计算以p为参数的直线和原始数据之间的误差
def residuals(p):
        k, b = p
        return Y-(k*X+b)

# leastsq()使得residuals()的输出数组的平方和最小，参数的初始值为[1, 0]
ret = optimize.leastsq(residuals, [1, 10])
k, b = ret[0]
print("k = ", k, "b = ", b)
# k =  0.4211697393502931 b =  -8.288302606523974

可视化：

import matplotlib.pyplot as plt

# 画样本点
plt.figure(figsize=(8, 6))
plt.scatter(X, Y, color="green", label="Samples", linewidth=2)

# 画拟合直线
x = np.linspace(150, 190, 100) # 在150-190直接画100个连续点
y = k*x + b # 函数式
plt.plot(x,y,color="red", label="Fit",linewidth=2)
plt.legend() # 绘制图例
plt.show()

信号处理

scipy.signal专用于信号处理，scipy.signal.resample(input,n)函数使用FFT将信号重采样成n个点：

import numpy as np

t = np.linspace(0, 5, 100)
x = np.sin(t)

from scipy import signal
x_resampled = signal.resample(x, 25)

import matplotlib.pyplot as plt
plt.plot(t, x)
plt.plot(t[::4], x_resampled, 'ko')

plt.show()

scipy.signal.detrend()函数从信号中去除线性趋势：

import numpy as np

t = np.linspace(0, 5, 100)
x = t + np.random.normal(size=100)

from scipy import signal
x_detrended = signal.detrend(x)

import matplotlib.pyplot as plt
plt.plot(t, x) 
plt.plot(t, x_detrended) 

plt.show()

对于非线性滤波，scipy.signal中提供了中值滤波scipy.signal.medfilt()， 维纳滤波scipy.signal.wiener()；

统计

scipy.stats包含了统计工具以及概率分析工具；

均值是样本的平均值：
np.mean(samples) 
中位数是样本的中间值：
np.mean(samples) 
中位数也是百分位数50，因为50%的观察值低于它：
stats.scoreatpercentile(samples, 50) 
同样，可以计算百分位数90：
stats.scoreatpercentile(samples, 90)

统计检验是一种决策指标。例如，如果有两组观测值，假设是高斯过程产生的，可以用T检验来判断两组观测值的均值是否存在显著差异：

from scipy import stats

a = np.random.normal(0, 1, size=100)
b = np.random.normal(1, 1, size=10)
stats.ttest_ind(a, b)
# Ttest_indResult(statistic=-5.975001870244409, pvalue=3.002595376787366e-08)

产生的输出包括：

• T统计值/statistic：一个数字，其符号与两个随机过程的差值成正比，其大小与该差值的显著性有关；
• p值/pvalue：两个过程相同的概率；如果它接近1，这两个过程几乎肯定是相同的。越接近于零，这些过程就越有可能有不同均值。