C++ 偏微分数值计算库_AI 基础：Scipy（科学计算库） 简易入门

0.导语

Scipy是一个用于数学、科学、工程领域的常用软件包，可以处理插值、积分、优化、图像处理、常微分方程数值解的求解、信号处理等问题。它用于有效计算Numpy矩阵，使Numpy和Scipy协同工作，高效解决问题。

Scipy是由针对特定任务的子模块组成：

在此之前，我已经写了以下几篇AI基础的快速入门，本篇文章讲解科学计算库Scipy快速入门：

已发布：

AI 基础：Numpy 简易入门

AI 基础：Pandas 简易入门

AI基础：数据可视化简易入门(matplotlib和seaborn)

备注：本文代码可以在github下载

/fengdu78/Data-Science-Notes/tree/master/4.scipy

1.SciPy-数值计算库

import numpy

import matplotlib

import scipy

'1.0.0'

常数和特殊函数

from scipy

299792458.0

6.62607004e-34

"electron mass"]

(9.10938356e-31, 'kg', 1.1e-38)

# 1英里等于多少米, 1英寸等于多少米, 1克等于多少千克, 1磅等于多少千克

1609.3439999999998

0.0254

0.001

0.45359236999999997

import scipy.special

print (

1.0

0.0

1e-20

0.1,

[(200, 4), (200, 4), (200, 4), (200, 4)]

#%figonly=使用广播计算得到的`ellipj()`返回值

2.拟合与优化-optimize

非线性方程组求解

import pylab

import matplotlib

from math

[-0.70622057 -0.6 -2.5 ]

[0.0, -9.126033262418787e-14, 5.32907051851e-15]

def j(x):

[-0.70622057 -0.6 -2.5 ]

[0.0, -9.126033262418787e-14, 5.32907051851e-15]

最小二乘拟合

import numpy

k = 0.6134953491930442 b = 1.794092543259387

#%figonly=最小化正方形面积之和(左)，误差曲面(右)

#%fig=带噪声的正弦波拟合

真实参数: [10, 0.34, 0.5235987755982988]

拟合参数 [10.25218748 0.3423992 0.50817423]

def func2(x, A, k, theta):

[10.25218748 0.3423992 0.50817425]

10,

真实参数: [10, 0.34, 0.5235987755982988]

拟合参数 [ 0.71093469 1.02074585 -0.12776742]

计算函数局域最小值

def target_function(x, y):

Nelder-Mead : min= 5.30934e-10, f count=125, fprime count= 0, fhess count= 0

Powell: min= 0, f count= 52, fprime count= 0, fhess count= 0

CG: min= 9.63056e-21, f count= 39, fprime count= 39, fhess count= 0

BFGS : min= 1.84992e-16, f count= 40, fprime count= 40, fhess count= 0

Newton-CG : min= 5.22666e-10, f count= 60, fprime count= 97, fhess count= 38

L-BFGS-B : min= 6.5215e-15, f count= 33, fprime count= 33, fhess count= 0

#%figonly=各种优化算法的搜索路径

计算全域最小值

def func(x, p):

1,

[10.25218676 -0.34239909 2.63341581]

#%figonly=用`basinhopping()`拟合正弦波

3.线性代数-linalg

解线性方程组

import pylab

import matplotlib

import numpy

True

5.38 ms ± 120 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

8.14 ms ± 994 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

luf = linalg.lu_factor(A)

X3 = linalg.lu_solve(luf, B)

np.allclose(X1, X3)

True

1000,

50.6 ms ± 1.94 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

3.49 ms ± 306 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

20.1 ms ± 1.42 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

4.49 ms ± 65 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

最小二乘解

from numpy.lib.stride_tricks

1000,

Average error of H1: 0.0

Average error of H2: 0.2958422158342371

#%figonly=实际的系统参数与最小二乘解的比较

特征值和特征向量

1,

[1.13027756+0.j 0.76972244+0.j]

[[ 0.91724574 0.79325185]

[-0.3983218 0.60889368]]

#%figonly=线性变换将蓝色箭头变换为红色箭头

42)

2, x*y, y**

[-0.55214278 0.5580915 -0.23809922 0.54584559 -0.08350449 -0.14852803]

#%figonly=用广义特征向量计算的拟合椭圆

奇异值分解-SVD

"vinci_target.png"),

(505, 375)

U, s, Vh = linalg.svd(img)

print(U.shape)

print(s.shape)

print(Vh.shape)

(505, 505)

(375,)

(375, 375)

#%fig=按从大到小排列的奇异值

def composite(U, s, Vh, n):

True

#%fig=原始图像、使用10、20、50个向量合成的图像(从左到右)

%array_image img; img10; img20; img50

UsageError: Line magic function `%array_image` not found.

pl.imshow(img)

pl.imshow(img10)

pl.imshow(img20)

pl.imshow(img50)

4.统计-stats

import numpy

import matplotlib.pyplot

连续概率分布

from scipy

['ksone',

'kstwobign',

'norm',

'alpha',

'anglit',

'arcsine',

'beta',

'betaprime',

'bradford',

'burr',

'burr12',

'fisk',

'cauchy',

'chi',

'chi2',

'cosine',

'dgamma',

'dweibull',

'expon',

'exponnorm',

'exponweib',

'exponpow',

'fatiguelife',

'foldcauchy',

'f',

'foldnorm',

'weibull_min',

'weibull_max',

'frechet_r',

'frechet_l',

'genlogistic',

'genpareto',

'genexpon',

'genextreme',

'gamma',

'erlang',

'gengamma',

'genhalflogistic',

'gompertz',

'gumbel_r',

'gumbel_l',

'halfcauchy',

'halflogistic',

'halfnorm',

'hypsecant',

'gausshyper',

'invgamma',

'invgauss',

'invweibull',

'johnsonsb',

'johnsonsu',

'laplace',

'levy',

'levy_l',

'levy_stable',

'logistic',

'loggamma',

'loglaplace',

'lognorm',

'gilbrat',

'maxwell',

'mielke',

'kappa4',

'kappa3',

'nakagami',

'ncx2',

'ncf',

't',

'nct',

'pareto',

'lomax',

'pearson3',

'powerlaw',

'powerlognorm',

'powernorm',

'rdist',

'rayleigh',

'reciprocal',

'rice',

'recipinvgauss',

'semicircular',

'skewnorm',

'trapz',

'triang',

'truncexpon',

'truncnorm',

'tukeylambda',

'uniform',

'vonmises',

'vonmises_line',

'wald',

'wrapcauchy',

'gennorm',

'halfgennorm',

'crystalball',

'argus']

stats.norm.stats()

(array(0.), array(1.))

1.0, scale=

(array(1.), array(4.))

10000)

(1.0048352738823323, 3.937211773554)

# 得到随机序列期望值和标准差

(1.0048352738823323, 1.984240855341749)

100, normed=

max pdf error: 0.018998755595167102, max cdf error: 0.018503342378306975

#%figonly=正态分布的概率密度函数(左)和累积分布函数(右)

1.0))

(array(1.), array(1.))

(array(2.), array(2.))

2.0, scale=

(array(4.), array(8.))

2, scale=

array([4.40563983, 6.17699951, 3.65503843, 3.28052152])

2, scale=

array([0.12169605, 0.07037188, 0.14694352, 0.15904745])

2, scale=

array([0.12169605, 0.07037188, 0.14694352, 0.15904745])

离散概率分布

1,

dice = stats.rv_discrete(values=(x, p))

dice.rvs(size=20)

array([2, 5, 2, 6, 1, 6, 6, 5, 3, 1, 5, 2, 1, 1, 1, 1, 1, 2, 1, 6])

42)

核密度估计

#%fig=核密度估计能更准确地表示随机变量的概率密度函数

#%fig=`bw_method`参数越大核密度估计曲线越平滑

二项、泊松、伽玛分布

6),

array([4.01877572e-01, 4.01877572e-01, 1.60751029e-01, 3.21502058e-02,

3.21502058e-03, 1.28600823e-04])

#%fig=当n足够大时二项分布和泊松分布近似相等

0.00675531110335309

0.0006301754049777564

#%fig=模拟泊松分布

time=1000, max_error=0.0196423002718

time=50000, max_error=0.001798012894964722

#%fig=模拟伽玛分布

100000

202.3388747879705

60

199.99833251643057

#%figonly=观察者偏差

from scipy

200.0

学生 t-分布与 t 检验

#%fig=模拟学生t-分布

max error: 0.006832108369761447

#%figonly=当`df`增大，学生t-分布趋向于正态分布

30

0.5) / (np.std(s, ddof=

2.658584340882224 Ttest_1sampResult(statistic=2.658584340882224, pvalue=0.01263770225709123)

print ((np.mean(s) -

-1.1450173670383303

Ttest_1sampResult(statistic=-1.1450173670383303, pvalue=0.26156414618801477) Ttest_1sampResult(statistic=-0.3842970254542196, pvalue=0.7035619103425202)

#%fig=红色部分为`ttest_1samp()`计算的p值

from scipy

0.012633433707685974

200000

0.012695

42)

Ttest_indResult(statistic=-2.2391470627176755, pvalue=0.033250866086743665)

Ttest_indResult(statistic=-0.5946698521856172, pvalue=0.5551805875810539)

卡方分布和卡方检验

#%fig=使用随机数验证卡方分布

max error: 0.0030732520533635066

#%fig=模拟卡方分布

def choose_balls(probabilities, size):

[80 93 97 64 66]

[89 76 79 71 85]

chi1, p1 = stats.chisquare(counts1)

chi2, p2 = stats.chisquare(counts2)

print ("chi1 =", chi1, "p1 =", p1)

print ("chi2 =", chi2, "p2 =", p2)

chi1 = 11.375 p1 = 0.022657601239769634

chi2 = 2.55 p2 = 0.6357054527037017

#%figonly=卡方检验计算的概率为阴影部分的面积

43,

1.0724852071005921

0.300384770390566

stats.fisher_exact(table)

(0.43434343434343436, 0.23915695682224306)

5.数值积分-integrate

import pylab

球的体积

def half_circle(x):

10000

3.1415893269307373

2

3.1415893269315975

from scipy

3.141592653589797

def half_sphere(x, y):

-1,

2.094395102393199 1.0002356720661965e-09 2.0943951023931953

解常微分方程组

#%fig=洛伦茨吸引子：微小的初值差别也会显著地影响运动轨迹

ode 类

def mass_spring_damper(xu, t, m, k, b, F):

#%fig=滑块的速度和位移曲线

from scipy.integrate

True

class PID(object):

#%fig=使用PID控制器让滑块停在位移为1.0处

控制力的终值: 19.943404699515057

%%time

from scipy import optimize

def eval_func(k):

kp, ki, kd = k

t, F_arr, result = pid_control_system(kp, ki, kd, 0.01)

return np.sum(np.abs(result[:, 0] - 1.0))

kwargs = {

"method": "L-BFGS-B",

"bounds": [(10, 200), (10, 100), (1, 100)],

"options": {

"approx_grad": True

}

}

opt_k = optimize.basinhopping(

eval_func, (10, 10, 10), niter=10, minimizer_kwargs=kwargs)

print(opt_k.x)

[56.67106149 99.97434757 1.33963577]

Wall time: 55.1 s

#%fig=优化PID的参数降低控制响应时间

t=0.5000000000000002, x=1.07098, u=0.315352

6.信号处理-signal

import pylab

import matplotlib

中值滤波

#%fig=使用中值滤波剔除瞬间噪声

True

滤波器设计

8000.0

#%fig=用频率扫描波测量的频率响应

连续时间线性系统

#%fig=系统的阶跃响应和正弦波响应

Text(0,0.5,'位移(米)')

7.插值-interpolate

import numpy

import matplotlib

一维插值

WARNING：高次interp1d()插值的运算量很大，因此对于点数较多的数据，建议使用后面介绍的UnivariateSpline()。

#%fig=`interp1d`的各阶插值

外推和 Spline 拟合

#%fig=使用UnivariateSpline进行插值：外推(上)，数据拟合(下)

3))

[ 0.053 3.151 6.36 9.386 12.603 15.619 18.929]

#%fig=计算Spline与水平线的交点

参数插值

#%fig=使用参数插值连接二维平面上的点

单调插值

import numpy

多维插值

#%fig=使用interp2d类进行二维插值

griddata

WARNING

griddata()使用欧几里得距离计算插值。如果 K 维空间中每个维度的取值范围相差较大，则应先将数据正规化，然后使用griddata()进行插值运算。

#%fig=使用gridata进行二维插值

径向基函数插值

#%fig=一维RBF插值

for fname, rbf

multiquadric [-0.88822885 2.17654513 1.42877511 -2.67919021]

gaussian [ 1.00321945 -0.02345964 -0.65441716 0.91375159]

linear [-0.26666667 0.6 0.73333333 -0.9 ]

#%fig=二维径向基函数插值

#%fig=`epsilon`参数指定径向基函数中数据点的作用范围

8.稀疏矩阵-sparse

%matplotlib inline

import numpy as np

import pylab as pl

from scipy import sparse

from scipy.sparse import csgraph

稀疏矩阵的储存形式

from scipy

dict_keys([(2, 3), (3, 3), (4, 3)])

dict_values([1.0, 2.0, 3.0])

10,

[list([]) list([]) list([1.0]) list([3.0, 2.0]) list([]) list([]) list([])

list([]) list([]) list([])]

[list([]) list([]) list([3]) list([2, 4]) list([]) list([]) list([])

list([]) list([]) list([])]

2,

[3 4 2 3] [2 3 3 2] [ 1 2 3 10]

[[ 0 0 0 0 0 0]

[ 0 0 0 0 0 0]

[ 0 0 0 11 0 0]

[ 0 0 3 0 2 0]

[ 0 0 0 0 0 0]]

矩阵向量相乘

import numpy

array([ 1, -3, -1], dtype=int32)

示例1

构造一个1000x1000 lil_matrix并添加值：

from scipy.sparse

1000,

现在将其转换为CSR格式，并求解的：

A = A.tocsr()

b = rand(1000)

x = spsolve(A, b)

将其转换为密集矩阵并求解，并检查结果是否相同：

x_ = solve(A.toarray(), b)

现在我们可以使用以下公式计算错误的范数：

err = norm(x-x_)

err < 1e-10

True

示例2

构造COO格式的矩阵：

from scipy

注意，索引不需要排序。

转换为CSR或CSC时，将对重复的(i，j)条目进行求和。

0,

这对于构造有限元刚度矩阵和质量矩阵是有用的。

9.图像处理-ndimage

import numpy

形态学图像处理

import numpy

膨胀和腐蚀

#%fig=四连通和八连通的膨胀运算

#%fig=不同结构元素的膨胀效果

#%fig=四连通和八连通的腐蚀运算

Hit和Miss

#%fig=Hit和Miss运算

#%fig=使用Hit和Miss进行细线化运算

图像分割

"suqares.jpg")

from scipy.ndimage

各种距离值 [ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

24 25 26 27]

8).astype(np.uint8)

from scipy.ndimage.measurements

1-squares_core,

#%figonly=矩形区域分割算法各个步骤的输出图像

10.空间算法库-spatial

import numpy

import matplotlib

计算最近旁点

100))

0.5244435681885733 0.4982156075770372

from scipy

0.2

array([list([48]), list([37, 78]), list([22, 79, 92]), list([6, 35, 58]),

list([7, 42, 55, 83])], dtype=object)

0.1) - kd.query_pairs(

{(1, 46),

(3, 21),

(3, 82),

(3, 95),

(5, 16),

(9, 30),

(10, 87),

(11, 42),

(11, 97),

(18, 41),

(29, 74),

(32, 51),

(37, 78),

(39, 61),

(41, 61),

(50, 84),

(55, 83),

(73, 81)}

#%figonly=用cKDTree寻找近旁点

from scipy.spatial

(100, 100)

(100, 5)

print (dist[:,

[0.15188266 0.09595807 0.05009422 0.11180181 0.19015485]

[0.15188266 0.09595807 0.05009422 0.11180181 0.19015485]

dist1[np.diag_indices(len(points))] = np.inf

nearest_pair = np.unravel_index(np.argmin(dist1), dist1.shape)

print (nearest_pair, dist1[nearest_pair])

(22, 92) 0.005346210248158245

2)

[22 92] 0.005346210248158245

1000000

def naive_count_at(start, end, time):

#%figonly=使用二维K-d树搜索指定区间的在线用户

class KdSearch(object):

True

QUESTION

请读者研究点数N和leafsize参数与创建 K-d 树和搜索时间之间的关系。

凸包

42)

[[2 5]

[2 6]

[0 5]

[1 6]

[1 0]]

[5 2 6 1 0]

#%fig=二维平面上的凸包

42)

(38, 3)

沃罗诺伊图

0.2,

print(vo.vertices); print(vo.regions); print(vo.ridge_vertices)

[[0.5 0.045 ]

[0.245 0.351 ]

[0.755 0.351 ]

[0.3375 0.425 ]

[0.6625 0.425 ]

[0.45 0.65 ]

[0.55 0.65 ]]

[[-1, 0, 1], [-1, 0, 2], [], [6, 4, 3, 5], [5, -1, 1, 3], [4, 2, 0, 1, 3], [6, -1, 2, 4], [6, -1, 5]]

[[-1, 0], [0, 1], [-1, 1], [0, 2], [-1, 2], [3, 5], [3, 4], [4, 6], [5, 6], [1, 3], [-1, 5], [2, 4], [-1, 6]]

-100,

#%figonly=沃罗诺伊图将空间分割为多个区域

print (vo.point_region)

[0 3 1 7 4 6 5]

[6, -1, 2, 4]

德劳内三角化

46.445,

[[8 5 7]

[1 5 3]

[5 6 7]

[6 0 7]

[0 6 2]

[6 1 2]

[1 6 5]

[9 5 8]

[4 9 8]

[5 9 3]

[9 4 3]]

[[104.58977484 127.03566055]

[235.1285198.68143374]

[107.83960707 155.53682482]

[ 71.22104881 228.39479887]

[110.3105291.17642838]

[201.40695449 227.68436282]

[201.61895891 226.21958623]

[152.96231864 93.25060083]

[205.40381294 -90.5480267 ]

[235.1285127.45701644]

[267.91709907 107.6135 ]]

#%fig=德劳内三角形的外接圆与圆心

往期精彩回顾

那些年做的学术公益-你不是一个人在战斗

适合初学者入门人工智能的路线及资料下载

机器学习在线手册

深度学习在线手册

备注：加入本站微信群或者qq群，请回复“加群”

加入知识星球(4500+用户，ID：92416895)，请回复“知识星球”

如果觉得《C++ 偏微分数值计算库_AI 基础：Scipy（科学计算库） 简易入门》对你有帮助，请点赞、收藏，并留下你的观点哦！