Scipy:高端科学计算
作者:Adrien Chauve, Andre Espaze, Emmanuelle Gouillart, Gaël Varoquaux, Ralf Gommers
译者表示最后部分没怎么看懂,此文档维护中……
Scipy
scipy包包含致力于科学计算中常见问题的各个工具箱。它的不同子模块相应于不同的应用。像插值,积分,优化,图像处理,,特殊函数等等。
scipy可以与其它标准科学计算程序库进行比较,比如GSL(GNU C或C++科学计算库),或者Matlab工具箱。scipy是Python中科学计算程序的核心包;它用于有效地计算numpy矩阵,来让numpy和scipy协同工作。
在实现一个程序之前,值得检查下所需的数据处理方式是否已经在scipy中存在了。作为非专业程序员,科学家总是喜欢重新发明造轮子,导致了充满漏洞的,未经优化的,很难分享和维护的代码。相反,Scipy程序经过优化和测试,因此应该尽可能使用。
目录
- 文件输入/输出:scipy.io
- 特殊函数:scipy.special
- 线性代数运算:scipy.linalg
- 快速傅里叶变换:scipy.fftpack
- 优化和拟合:scipy.optimize
- 统计和随机数: scipy.stats
- 插值:scipy.interpolate
- 数值积分:scipy.integrateFusy,
- 信号处理:scipy.signal
- 图像处理:scipy.ndimage
- 总结练习
- Footnotes
警告:这个教程离真正的数值计算介绍很远。因为枚举scipy中不同的子模块和函数非常无聊,我们集中精力代之以几个例子来给出如何使用`scipy`进行计算的大致思想。
scipy 由一些特定功能的子模块组成:
模块 | 功能 |
---|---|
scipy.cluster | 矢量量化 / K-均值 |
scipy.constants | 物理和数学常数 |
scipy.fftpack | 傅里叶变换 |
scipy.integrate | 积分程序 |
scipy.interpolate | 插值 |
scipy.io | 数据输入输出 |
scipy.linalg | 线性代数程序 |
scipy.ndimage | n维图像包 |
scipy.odr | 正交距离回归 |
scipy.optimize | 优化 |
scipy.signal | 信号处理 |
scipy.sparse | 稀疏矩阵 |
scipy.spatial | 空间数据结构和算法 |
scipy.special | 任何特殊数学函数 |
scipy.stats | 统计 |
它们全依赖numpy,但是每个之间基本独立。导入Numpy和这些scipy模块的标准方式是:
import numpy as np
from scipy import stats # 其它子模块相同
主scipy
命名空间大多包含真正的numpy函数(尝试 scipy.cos 就是 np.cos)。这些仅仅是由于历史原因,通常没有理由在你的代码中使用import
scipy
文件输入/输出:scipy.io
-
导入和保存matlab文件:
In [1]: from scipy import io as spio In [3]: import numpy as np In [4]: a = np.ones((3, 3)) In [5]: spio.savemat('file.mat', {'a': a}) # savemat expects a dictionary /usr/lib/python2.7/site-packages/scipy/io/matlab/mio.py:266: FutureWarning: Using oned_as default value ('column') This will change to 'row' in future versions oned_as=oned_as) In [6]: data = spio.loadmat('file.mat', struct_as_record=True) In [7]: data['a'] Out[7]: array([[ 1., 1., 1.], [ 1., 1., 1.], [ 1., 1., 1.]])
-
读取图片:
In [16]: from scipy import misc In [17]: misc.imread('scikit.png') Out[17]: array([[[255, 255, 255, 255], [255, 255, 255, 255], [255, 255, 255, 255], ..., [255, 255, 255, 255], [255, 255, 255, 255], [255, 255, 255, 255]], [[255, 255, 255, 255], [255, 255, 255, 255], [255, 255, 255, 255], ..., [255, 255, 255, 255], [255, 255, 255, 255], [255, 255, 255, 255]], [[255, 255, 255, 255], [255, 255, 255, 255], [255, 255, 255, 255], ..., [255, 255, 255, 255], [255, 255, 255, 255], [255, 255, 255, 255]], ..., [[255, 255, 255, 255], [255, 255, 255, 255], [255, 255, 255, 255], ..., [255, 255, 255, 255], [255, 255, 255, 255], [255, 255, 255, 255]], [[255, 255, 255, 255], [255, 255, 255, 255], [255, 255, 255, 255], ..., [255, 255, 255, 255], [255, 255, 255, 255], [255, 255, 255, 255]], [[255, 255, 255, 255], [255, 255, 255, 255], [255, 255, 255, 255], ..., [255, 255, 255, 255], [255, 255, 255, 255], [255, 255, 255, 255]]], dtype=uint8) In [18]: import matplotlib.pyplot as plt In [19]: plt.imread('scikit.png') Out[19]: array([[[ 1., 1., 1., 1.], [ 1., 1., 1., 1.], [ 1., 1., 1., 1.], ..., [ 1., 1., 1., 1.], [ 1., 1., 1., 1.], [ 1., 1., 1., 1.]], [[ 1., 1., 1., 1.], [ 1., 1., 1., 1.], [ 1., 1., 1., 1.], ..., [ 1., 1., 1., 1.], [ 1., 1., 1., 1.], [ 1., 1., 1., 1.]], [[ 1., 1., 1., 1.], [ 1., 1., 1., 1.], [ 1., 1., 1., 1.], ..., [ 1., 1., 1., 1.], [ 1., 1., 1., 1.], [ 1., 1., 1., 1.]], ..., [[ 1., 1., 1., 1.], [ 1., 1., 1., 1.], [ 1., 1., 1., 1.], ..., [ 1., 1., 1., 1.], [ 1., 1., 1., 1.], [ 1., 1., 1., 1.]], [[ 1., 1., 1., 1.], [ 1., 1., 1., 1.], [ 1., 1., 1., 1.], ..., [ 1., 1., 1., 1.], [ 1., 1., 1., 1.], [ 1., 1., 1., 1.]], [[ 1., 1., 1., 1.], [ 1., 1., 1., 1.], [ 1., 1., 1., 1.], ..., [ 1., 1., 1., 1.], [ 1., 1., 1., 1.], [ 1., 1., 1., 1.]]], dtype=float32)
参见:
- 载入txt文件:numpy.loadtxt()/numpy.savetxt()
- 智能导入文本/csv文件:numpy.genfromtxt()/numpy.recfromcsv()
- 高速,有效率但numpy特有的二进制格式:numpy.save()/numpy.load()
特殊函数:scipy.special
特殊函数是先验函数。scipy.special的文档字符串写得非常好,所以我们不在这里列出所有函数。常用的有:
- 贝塞尔函数,如
scipy.special.jn()
(整数n阶贝塞尔函数) - 椭圆函数(
scipy.special.ellipj()
雅可比椭圆函数,……) - 伽马函数:
scipy.special.gamma()
,还要注意scipy.special.gammaln
,这个函数给出对数坐标的伽马函数,因此有更高的数值精度。
线性代数运算:scipy.linalg
scipy.linalg模块提供标准线性代数运算,依赖于底层有效率的实现(BLAS,LAPACK)。
-
scipy.linalg.det()函数计算方阵的行列式:
In [22]: from scipy import linalg In [23]: arr = np.array([[1, 2], ....: [3, 4]]) In [24]: linalg.det(arr) Out[24]: -2.0 In [25]: linalg.det(np.ones((3,4))) --------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-25-375ad1d49940> in <module>() ----> 1 linalg.det(np.ones((3,4))) /usr/lib/python2.7/site-packages/scipy/linalg/basic.pyc in det(a, overwrite_a) 398 a1 = np.asarray_chkfinite(a) 399 if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]: --> 400 raise ValueError('expected square matrix') 401 overwrite_a = overwrite_a or _datacopied(a1, a) 402 fdet, = get_flinalg_funcs(('det',), (a1,)) ValueError: expected square matrix
py.linalg.inv()`函数计算方阵的逆:
In [26]: arr = np.array([[1, 2],
[3, 4]])
In [27]: iarr = linalg.inv(arr)
In [28]: iarr
Out[28]:
array([[-2. , 1. ],
[ 1.5, -0.5]])
In [29]: np.allclose(np.dot(arr, iarr), np.eye(2))
Out[29]: True
最后计算奇异阵的逆(它的行列式为0)将会引发(raise)LinAlgError
:
In [32]: arr = np.array([[3, 2],
[6, 4]])
In [33]: linalg.inv(arr)
---------------------------------------------------------------------------
LinAlgError Traceback (most recent call last)
<ipython-input-33-52c04c854a80> in <module>()
----> 1 linalg.inv(arr)
/usr/lib/python2.7/site-packages/scipy/linalg/basic.pyc in inv(a, overwrite_a)
346 inv_a, info = getri(lu, piv, overwrite_lu=1)
347 if info > 0:
--> 348 raise LinAlgError("singular matrix")
349 if info < 0:
350 raise ValueError('illegal value in %d-th argument of internal '
LinAlgError: singular matrix
-
还有更多高级运算,如奇异值分解(SVD):
In [34]: arr = np.arange(9).reshape((3, 3)) + np.diag([1, 0, 1]) In [35]: uarr, spec, vharr = linalg.svd(arr)
它的结果数组谱是:
In [36]: spec Out[36]: array([ 14.88982544, 0.45294236, 0.29654967])
原始矩阵可以由svd的输出用
np.dot
点乘重新组合得到:In [37]: sarr = np.diag(spec) In [38]: svd_mat = uarr.dot(sarr).dot(vharr) In [39]: np.allclose(svd_mat, arr) Out[39]: True
SVD在信号处理和统计中运用很广。许多其它标准分解(QR,LU,Cholesky,Schur),还有线性系统的解也可以从scipy.linalg中获得。
快速傅里叶变换:scipy.fftpack
scipy.fftpack模块用来计算快速傅里叶变换。作为示例,一个(噪声)输入信号可能像这样:
In [40]: time_step = 0.02
In [41]: period = 5
In [42]: time_vec = np.arange(0, 20, time_step)
In [43]: sig = np.sin(2 * np.pi / period * time_vec) + \
....: 0.5 * np.random.randn(time_vec.size)
观测者并不指导信号频率,仅仅等间隔取样信号sig。信号应该来自一个真实的函数所以傅里叶变换将是对称的。scipy.fftpack.fftfreq()函数将生成取样频率,scipy.fftpack.fft()将计算快速傅里叶变换:
因为功率结果是对称的,仅仅需要使用谱的正值部分来找出频率:
In [48]: pidxs = np.where(sample_freq > 0)
In [49]: freqs = sample_freq[pidxs]
In [50]: power = np.abs(sig_fft)[pidxs]
信号频率可以这样被找到:
In [51]: freq = freqs[power.argmax()]
In [52]: np.allclose(freq, 1./period)
Out[52]: True
现在高频噪声将被从傅里叶变换信号中移除:
In [53]: sig_fft[np.abs(sample_freq) > freq] = 0
得到滤波信号,可以用scipy.fftpack.ifft
函数计算:
In [54]: main_sig = fftpack.ifft(sig_fft)
结果可以这样可视化:
In [55]: plt.figure()
Out[55]: <matplotlib.figure.Figure at 0x4a9fb50>
In [56]: plt.plot(time_vec, sig)
Out[56]: [<matplotlib.lines.Line2D at 0x4ad3790>]
In [57]: plt.plot(time_vec, main_sig, linewidth=3)
/usr/lib/python2.7/site-packages/numpy/core/numeric.py:320: ComplexWarning: Casting complex values to real discards the imaginary part
return array(a, dtype, copy=False, order=order)
Out[57]: [<matplotlib.lines.Line2D at 0x4ad3dd0>]
In [58]: plt.xlabel('Time [s]')
Out[58]: <matplotlib.text.Text at 0x4aad050>
In [59]: plt.ylabel('Amplitude')
Out[59]: <matplotlib.text.Text at 0x4aadbd0>
In [60]: plt.show()
Numpy也有一个FFT实现(numpy.fft)。然而,通常scipy的应该优先使用,因为它使用了更有效率的底层实现。
工作示例:找到原始周期
工作示例:高斯图像模糊
卷积:
练习:登月图片消噪
- 检查提供的图像moonlanding.png,该图像被周期噪声严重污染了。在这个练习中,我们旨在使用快速傅里叶变换清除噪声。
- 用
plt.imread
加载图像。 - 使用
scipy.fftpack
中的2-D傅里叶函数找到并绘制图像的谱线(傅里叶变换)。可视化这个谱线对你有问题吗?如果有,为什么? - 这个谱包含高频和低频成分。噪声是在谱线的高频部分中,所以设置一些成分为0(使用数组切片)。
- 应用逆傅里叶变换来看最后的图像。
优化和拟合:scipy.optimize
优化是找到最小值或等式的数值解的问题。
scipy.optimization
子模块提供了函数最小值(标量或多维)、曲线拟合和寻找等式的根的有用算法。
from scipy import optimize
找到标量函数的最小值
让我们定义以下函数
In [2]: def f(x):
...: return x**2 + 10 * np.sin(x)
然后绘制它:
In [3]: x = np.arange(-10, 10, 0.1)
In [4]: plt.plot(x, f(x))
Out[4]: [<matplotlib.lines.Line2D at 0x3e2a4d0>]
In [5]: plt.show()
该函数在大约-1.3有个全局最小值,在3.8有个局部最小值。
找到这个函数最小值一般而有效的方法是从初始点使用梯度下降法。BFGS算法1是做这个的好方法:
In [6]: optimize.fmin_bfgs(f, 0)
Optimization terminated successfully.
Current function value: -7.945823
Iterations: 5
Function evaluations: 24
Gradient evaluations: 8
Out[6]: array([-1.30644003])
这个方法一个可能的问题在于,如果函数有局部最小值,算法会因初始点不同找到这些局部最小而不是全局最小:
In [7]: optimize.fmin_bfgs(f, 3, disp=0)
Out[7]: array([ 3.83746663])
如果我们不知道全局最小值的邻近值来选定初始点,我们需要借助于耗费资源些的全局优化。为了找到全局最小点,最简单的算法是蛮力算法2,该算法求出给定格点的每个函数值。
In [10]: grid = (-10, 10, 0.1)
In [11]: xmin_global = optimize.brute(f, (grid,))
In [12]: xmin_global
Out[12]: array([-1.30641113])
对于大点的格点,scipy.optimize.brute()
变得非常慢。scipy.optimize.anneal()
提供了使用模拟退火的替代函数。对已知的不同类别全局优化问题存在更有效率的算法,但这已经超出scipy的范围。一些有用全局优化软件包是OpenOpt、IPOPT、PyGMO和PyEvolve。
为了找到局部最小,我们把变量限制在(0, 10)
之间,使用scipy.optimize.fminbound()
:
In [13]: xmin_local = optimize.fminbound(f, 0, 10)
In [14]: xmin_local
Out[14]: 3.8374671194983834
注意:在高级章节部分数学优化:找到函数最小值中有关于寻找函数最小值更详细的讨论。
找到标量函数的根
为了寻找根,例如令f(x)=0
的点,对以上的用来示例的函数f
我们可以使用scipy.optimize.fsolve()
:
In [17]: root = optimize.fsolve(f, 1) # 我们的初始猜测是1
In [18]: root
Out[18]: array([ 0.])
注意仅仅一个根被找到。检查f的图像在-2.5附近有第二个根。我们可以通过调整我们的初始猜测找到这一确切值:
In [19]: root = optimize.fsolve(f, -2.5)
In [20]: root
Out[20]: array([-2.47948183])
曲线拟合
假设我们有从被噪声污染的f中抽样到的数据:
In [21]: xdata = np.linspace(-10, 10, num=20)
In [22]: ydata = f(xdata) + np.random.randn(xdata.size)
如果我们知道函数形式(当前情况是x^2 + sin(x)
),但是不知道幅度。我们可以通过最小二乘拟合拟合来找到幅度。首先我们定义一个用来拟合的函数:
In [23]: def f2(x, a, b):
....: return a*x**2 + b*np.sin(x)
然后我们可以使用scipy.optimize.curve_fit()
来找到a和b:
In [24]: guess = [2, 2]
In [25]: params, params_covariance = optimize.curve_fit(f2, xdata, ydata, guess)
In [26]: params
Out[26]: array([ 1.00439471, 10.04911441])
现在我们找到了f的最小值和根并且对它使用了曲线拟合。我们将一切放在一个单独的图像中:
注意:Scipy>=0.11中提供所有最小化和根寻找算法的