精通-Python-科学计算-二-

精通 Python 科学计算（二）

原文：Mastering Python Scientific Computing

协议：CC BY-NC-SA 4.0

六、Python 在符号计算中的应用

SymPy 包括从基本符号算术到多项式、微积分、解算器、离散数学、几何、统计和物理的功能。它主要处理三种类型的数字，即整数、实数和有理数。整数是没有小数点的整数，而实数是有小数点的数字。有理数有两部分:分子和分母。定义有理数，可以使用Ration类，需要两个数字。在这一章中，我们将借助示例程序讨论 SymPy 的概念。

我们将在本章中讨论以下主题:

• 使用 SymPy 的计算机化代数系统
• 核心能力和高级功能
• 多项式、微积分和解方程
• 离散数学、矩阵、几何、绘图、物理和统计
• 打印功能

让我们开始讨论 SymPy 及其核心能力，包括基本算术、扩展、简化、替换、模式匹配和各种函数(例如，指数、对数、方程根、三角函数、双曲函数和特殊函数)。

符号、表达式和基本算术

在 SymPy 中，我们需要定义符号，然后才能在任何表达式中使用它们。定义一个符号非常简单。我们只需要使用Symbol类中的symbol函数来定义一个符号，如下图所示。我们可以使用evalf() / n()方法得到任意物体的float数值近似值。

以下程序使用三种方法创建符号。对于只创建一个符号，方法的名称是符号，对于创建多个符号，方法的名称是符号。创建多个符号有两种方法:一种是通过将空格分隔的符号名称传递给符号方法，另一种是通过将m0:5传递给符号方法来创建一系列符号，如m0、m1、m2、m3、m4。在第二个选项中，m0是索引的第一个值，:后面的数字5表示总共应该创建五个这样的变量。

通常，两个整数的除法会截断小数部分。为了避免这种情况，以下程序的第一行强制对两个整数执行实际浮点除法。这就是程序最后一行将3.142857142857143存储在y的原因。如果我们忽略下面程序中的第一行。那么y的价值将是3:

from __future__ import division
from sympy import *

x, y, z = symbols('x y z')
m0, m1, m2, m3, m4 = symbols('m0:5')
x1 = Symbol('x1')
x1 + 500
y=22/7

下面的程序使用evalf()和n()方法将任何 SymPy 对象数值近似为一个float值。默认精度最高为 15 位小数，通过向这些方法传递一个参数，可以将其更改为任何所需的精度:

from __future__ import division
from sympy import, sin, pi

x=sin(50)

pi.evalf()
pi.evalf(50)

x.n()
x.n(20)

下一个程序演示了表达式的概念以及可以使用collect、expand、factor、simplify和subs方法对表达式执行的各种操作:

from sympy import collect, expand, factor, simplify
from sympy import Symbol, symbols
from sympy import sin, cos

x, y, a, b, c, d = symbols('x y a b c d')

expr = 5*x**2+2*b*x**2+cos(x)+51*x**2
simplify(expr)

factor(x**2+x-30)
expand ( (x-5) * (x+6) )

collect(x**3 + a*x**2 + b*x**2 + c*x + d, x)

expr = sin(x)*sin(x) + cos(x)*cos(x)
expr
expr.subs({x:5, y:25})
expr.subs({x:5, y:25}).n()

方程求解

有一个神奇的功能叫做solve。它可以求解所有类型的方程。这个函数返回一个方程的解。它需要两个参数:要求解的表达式和变量。下面的程序使用这个函数来求解各种类型的方程。在下面的等式中，假设等式的右侧为零:

from sympy import solve

solve (6*x**2 - 3*x - 30,x)

a, b, c = symbols('a b c')
solve( a*x**2 + b*x + c, x)
substitute_solution = solve( a*x**2 + b*x + c, x)
[ substitute_solution[0].subs({'a':6,'b':-3,'c':-30}), substitute_solution[1].subs({'a':6,'b':-3,'c':-30}) ]

solve([2*x + 3*y - 3, x -2* y + 1], [x, y])
)

为了求解方程组，我们有另一种形式的solve方法，该方法将方程列表作为第一个参数，将未知数列表作为第二个参数。这里展示了这一点:

from sympy import solve

solve ([2*x + y - 4, 5*x - 3*y],[x, y])
solve ([2*x + 2*y - 1, 2*x - 4*y],[x, y])

有理数、指数和对数的函数

SymPy 有许多处理有理数的函数。这些函数对有理数执行各种操作，包括简化、扩展、合并、拆分等等。SymPy 还支持指数和对数运算的几个函数。有三个对数函数:log(用来计算基数-b 的对数)ln(用来计算自然对数)log10(用来计算基数-10 的对数)。log函数需要两个参数:变量和基数。如果没有传递基数，那么默认情况下，这个函数会计算变量的自然对数，相当于ln。计算两个有理数的相加，我们使用together函数。类似地，要用分母除一个有理表达式的分子，我们使用apart函数，该函数用于以下程序:

from sympy import together, apart, symbols
x1, x2, x3, x4 = symbols('x1 x2 x3 x4')
x1/x2 + x3/x4
together(x1/x2 + x3/x4)

apart ((2*x**2+3*x+4)/(x+1))
together(apart ((2*x**2+3*x+4)/(x+1)))

exp(1)
log(4).n()
log(4,4).n()
ln(4).n()
mpmath.log10(4)

多项式

SymPy 允许我们定义和执行多项式的各种运算。我们也可以找到多项式的根。我们已经介绍了simplify、expand、factor和solve方法。这些方法也执行多项式的功能。要检查两个多项式是否相等，我们应该使用simplify函数:

from sympy import *

p, q = symbols ('p q')
p = (x+4)*(x+2)
q = x**2 + 6*x + 8
p == q # Unsuccessful
p - q == 0 # Unsuccessful 
simplify(p - q) == 0

三角学和复数

三角函数的输入大多是弧度角，而反三角函数返回弧度角。该模块还提供了度数到弧度、弧度到度数的转换功能。除了基本的三角函数，如sin、cos、tan，SymPy 还有三角化简和展开函数。

SymPy 还支持复数，以应对不存在实数解的情况。例如，考虑这个等式: x**2+4=0 。这个方程没有实数解；其解决方案将是 -2I* 或 +2I* 。这个 I 表示-1 的平方根。以下程序演示了三角函数，并以复数形式给出了该方程的解:

from sympy import *
x, y = symbols('x y')
expr = sin(x)*cos(y)+cos(x)*sin(y)
expr_exp= exp(5*sin(x)**2+4*cos(x)**2)

trigsimp(expr)
trigsimp(expr_exp)
expand_trig(sin(x+y))
solve(x**2+4,x) #complex number as solution

线性代数

SymPy 线性代数模块是另一个非常简单的模块，它为矩阵操作提供了易于学习的功能。它具有执行各种矩阵运算的功能，包括快速特殊矩阵创建、特征值、特征向量、转置、行列式和逆。快速创建特殊矩阵有三种方法，即eye、zeros和ones。eye方法创建一个单位矩阵，而zeros和ones创建所有元素分别等于 0 或 1 的矩阵。如果需要，我们可以从矩阵中删除选定的行和列。基本算术运算符，如 + 、 - 、 *** 、 **** 也对矩阵起作用:

from sympy import *
A = Matrix( [[1, 2, 3, 4],
       [5, 6, 7, 8],
       [ 9, 10, 11, 12],
       [ 13, 14, 15, 16]] )
A.row_del(3)
A.col_del(3)

A[0,1] # display row 0, col 1 of A
A[0:2,0:3] # top-left submatrix(2x3)

B = Matrix ([[1, 2, 3],
       [5, 6, 7],
       [ 9, 10, 11]] )
A.row_join(B)
B.col_join(B)
A + B
A - B
A * B
A **2 
eye(3) # 3x3 identity matrix
zeros(3, 3) # 3x3 matrix with all elements Zeros
ones(3, 3) # 3x3 matrix with all elements Ones

A.transpose() # It is same as A.T
M = Matrix( [[1, 2, 3],
    [4, 5, 6],
    [7, 8, 10]] )
M.det()

默认情况下，矩阵的逆矩阵是通过高斯消去法计算的，我们可以指定通过 LU 分解计算逆矩阵。SymPy 中的矩阵有计算简化行梯次形式的方法(rref方法)和零空间法(nullspace方法)。如果 A 是矩阵，那么nullspace就是所有向量的集合 v ，使得 A v=0 。也可以对矩阵元素进行替换操作；我们可以用符号条目创建一个矩阵，并用实际值和其他符号替换它们。我们还可以执行特殊操作，例如 QR 分解、Gram-Schmidt 正交化器和 LU 分解:

from sympy import *
A = Matrix( [[1,2],
  [3,4]] )
A.inv()
A.inv()*A
A*A.inv()
A = Matrix( [[ 1, -2],
   [-2, 3]] )
A.eigenvals() # same as solve( det(A-eye(2)*x), x)
A.eigenvects()
A.rref() 

A.nullspace()

x = Symbol('x')
M = eye(3) * x
M.subs(x, 4)
y = Symbol('y')
M.subs(x, y)

M.inv()
M.inv("LU")  

A = Matrix([[1,2,1],[2,3,3],[1,3,2]])
Q, R = A.QRdecomposition()
Q

M = [Matrix([1,2,3]), Matrix([3,4,5]), Matrix([5,7,8])]
result1 = GramSchmidt(M)
result2 = GramSchmidt(M, True)

微积分

微积分涉及为研究任何函数的各种性质而执行的操作，包括变化率、函数的极限行为以及函数图下面积的计算。在这一节中，你将学习极限、导数、级数求和和积分的概念。以下程序使用极限函数来解决简单的极限问题:

from sympy import limit, oo, symbols,exp, cos

oo+5
67000 < oo
10/oo

x , n = symbols ('x n')
limit( ((x**n - 1)/ (x - 1) ), x, 1)

limit( 1/x**2, x, 0)
limit( 1/x, x, 0, dir="-")

limit(cos(x)/x, x, 0)
limit(sin(x)**2/x, x, 0)
limit(exp(x)/x,x,oo)

可以使用diff(func_to_be_differentiated, variable)原型的diff函数来区分任何症状表达。以下程序使用diff函数计算各种症状表达式的微分:

from sympy import diff, symbols, Symbol, exp, dsolve, subs, Function

diff(x**4, x)
diff( x**3*cos(x), x )
diff( cos(x**2), x )
diff( x/sin(x), x )
diff(x**3, x, 2)
diff( exp(x), x)

dsolve方法帮助我们求解任何一种常微分方程和常微分方程组。本程序演示了dsolve在常微分方程和边值问题中的应用:

x = symbols('x')
f = symbols('f', cls=Function) 
dsolve( f(x) - diff(f(x),x), f(x) )

#ics argument can be used to pass the boundary condition as a dictionary to dsolve method
from sympy import *
x=symbols('x')
f=symbols('f', cls=Function)
dsolve(Eq(f(x).diff(x,x), 400*(f(x)-1+2*x)), f(x), ics={f(0):5, f(2):10})
# the above line will results in f(x) = C1*e^−30x + C2*e^30x − 2x + 1

下面的程序求函数的临界点 f(x)=4x3-3x2+2x 用二阶导数求函数在区间*【0，1】*内的最大值:

x = Symbol('x')
f = 4*x**3-3*x**2+2*x
diff(f, x)
sols = solve( diff(f,x), x)
sols
diff(diff(f,x), x).subs( {x:sols[0]} )
diff(diff(f,x), x).subs( {x:sols[1]} )

在 SymPy 中，我们可以使用integrate函数计算定积分和不定积分。这是一个计算定积分和不定积分的程序。它将象征性地定义这些积分。为了计算实际值，我们调用积分的n()方法，如本程序最后一行所示:

from sympy import *
integrate(x**3+1, x)
integrate(x*sin(x), x)
integrate(x+ln(x), x)

F = integrate(x**3+1, x)
F.subs({x:1}) - F.subs({x:0})

integrate(x**3-x**2+x, (x,0,1))    # definite Integrals 
integrate(sin(x)/x, (x,0,pi))    # definite Integrals 
integrate(sin(x)/x, (x,pi,2*pi))  # definite Integrals 
integrate(x*sin(x)/(x+1), (x,0,2*pi)) # definite Integrals
integrate(x*sin(x)/(x+1), (x,0,2*pi)).n()

序列是取整数的函数，我们可以通过为它的 n 项指定一个表达式来定义一个序列。我们也可以用期望值来代替。下面的程序使用 SymPy 中的一些简单序列演示了序列的概念:

from sympy import *
s1_n = 1/n
s2_n = 1/factorial(n)
s1_n.subs({n:5})
[ s1_n.subs({n:index1}) for index1 in range(0,8) ]
[ s2_n.subs({n:index1}) for index1 in range(0,8) ]
limit(s1_n, n, oo)
limit(s2_n, n, oo)

其项包含变量不同阶幂的级数称为幂级数，如泰勒级数、指数级数或正弦/余弦级数。这里有一个程序，可以计算一些包含特殊函数的序列。它还使用幂级数的概念:

from sympy import *
s1_n = 1/n
s2_n = 1/factorial(n)
summation(s1_n, [n, 1, oo])
summation(s2_n, [n, 0, oo])
import math
def s2_nf(n):
  return 1.0/math.factorial(n)

sum( [s2_nf(n) for n in range(0,10)] )
E.evalf()

exponential_xn = x**n/factorial(n)
summation( exponential_xn.subs({x:5}), [x, 0, oo] ).evalf()
exp(5).evalf()
summation( exponential_xn.subs({x:5}), [x, 0, oo])

import math # redo using only python
def exponential_xnf(x,n):
  return x**n/math.factorial(n)
sum( [exponential_xnf(5.0,i) for i in range(0,35)] )

series( sin(x), x, 0, 8)
series( cos(x), x, 0, 8)
series( sinh(x), x, 0, 8)
series( cosh(x), x, 0, 8)
series(ln(x), x, 1, 6) # Taylor series of ln(x) at x=1
series(ln(x+1), x, 0, 6) # Maclaurin series of ln(x+1)

向量

实数上定义的一个 n 元组也可以称为向量。在物理和数学中，矢量是一种数学对象，它有大小、大小或长度以及方向。在 SymPy 中，向量表示为1×n行矩阵或n×1列矩阵。下面的程序演示了 SymPy 中向量计算的概念。它计算向量的转置和范数:

from sympy import *
u = Matrix([[1,2,3]]) # a row vector = 1x3 matrix
v = Matrix([[4],
[5],    # a col vector = 3x1 matrix
[6]])
v.T # use the transpose operation to
u[1] # 0-based indexing for entries
u.norm() # length of u
uhat = u/u.norm() # unit-length vec in same dir as u
uhat
uhat.norm()

下一个程序演示了点积、叉积和向量投影运算的概念:

from sympy import *
u = Matrix([ 1,2,3])
v = Matrix([-2,3,3])
u.dot(v)

acos(u.dot(v)/(u.norm()*v.norm())).evalf()
u.dot(v) == v.dot(u)
u = Matrix([2,3,4])
n = Matrix([2,2,3])
(u.dot(n) / n.norm()**2)*n  # projection of v in the n dir

w = (u.dot(n) / n.norm()**2)*n
v = u - (u.dot(n)/n.norm()**2)*n # same as u - w
u = Matrix([ 1,2,3])
v = Matrix([-2,3,3])
u.cross(v)
(u.cross(v).norm()/(u.norm()*v.norm())).n()

u1,u2,u3 = symbols('u1:4')
v1,v2,v3 = symbols('v1:4')
Matrix([u1,u2,u3]).cross(Matrix([v1,v2,v3]))
u.cross(v)
v.cross(u)

物理模块

物理模块包含从物理上解决问题所需的功能。有几个物理子模块，用于执行与矢量物理、经典力学、量子力学、光学等相关的活动。

氢波函数

这个类别下有两个功能。第一种是以哈特利原子单位计算状态 (n，l) 的能量。另一个用哈特利原子单位计算相对论性态能量 (n，l，自旋)。以下程序演示了这些功能的使用:

from sympy.physics.hydrogen import E_nl, E_nl_dirac, R_nl
from sympy import var

var("n Z")
var("r Z")
var("n l")
E_nl(n, Z)
E_nl(1)
E_nl(2, 4)

E_nl(n, l)
E_nl_dirac(5, 2) # l should be less than n
E_nl_dirac(2, 1)
E_nl_dirac(3, 2, False)
R_nl(5, 0, r) # z = 1 by default
R_nl(5, 0, r, 1)

矩阵与泡利代数

在physics.matrices模块中有几个与物理相关的矩阵。以下程序演示了如何获得这些矩阵和泡利代数:

from sympy.physics.paulialgebra import Pauli, evaluate_pauli_product
from sympy.physics.matrices import mdft, mgamma, msigma, pat_matrix

mdft(4) # expression of discrete Fourier transform as a matrix multiplication
mgamma(2) # Dirac gamma matrix in the Dirac representation
msigma(2) #  Pauli matrix with (1,2,3)

# Following line computer Parallel Axis Theorem matrix to translate the inertia matrix a distance of dx, dy, dz for a body of mass m.
pat_matrix(3, 1, 0, 0)   

evaluate_pauli_product(4*x*Pauli(3)*Pauli(2))

一维和三维的量子谐振子

该模块具有用于计算一维谐振子、三维各向同性谐振子、一维谐振子的波函数和三维各向同性谐振子的径向波函数的函数。下面是一个使用本模块中可用功能的程序:

from sympy.physics.qho_1d import E_n, psi_n
from sympy.physics.sho import E_nl, R_nl
from sympy import var

var("x m omega")
var("r nu l")
x, y, z = symbols('x, y, z')

E_n(x, omega)
psi_n(2, x, m, omega)
E_nl(x, y, z)

R_nl(1, 0, 1, r)
R_nl(2, l, 1, r)

第二量化

用来分析和描述量子多体系统的概念是称为第二量子化。这个模块包含第二量子化算符的类和玻色子的状态。预定义符号可从sympy.abc导入:

from sympy.physics.secondquant import Dagger, B, Bd
from sympy.functions.special.tensor_functions import KroneckerDelta
from sympy.physics.secondquant import B, BKet, FockStateBosonKet
from sympy.abc import x, y, n
from sympy.abc import i, j, k
from sympy import Symbol
from sympy import I

Dagger(2*I)
Dagger(B(0))
Dagger(Bd(0))

KroneckerDelta(1, 2)
KroneckerDelta(3, 3)

#predefined symbols are available in abc including greek symbols like theta
KroneckerDelta(i, j)
KroneckerDelta(i, i)
KroneckerDelta(i, i + 1)
KroneckerDelta(i, i + 1 + k)

a = Symbol('a', above_fermi=True)
i = Symbol('i', below_fermi=True)
p = Symbol('p')
q = Symbol('q')
KroneckerDelta(p, q).indices_contain_equal_information
KroneckerDelta(p, q+1).indices_contain_equal_information
KroneckerDelta(i, p).indices_contain_equal_information

KroneckerDelta(p, a).is_above_fermi
KroneckerDelta(p, i).is_above_fermi
KroneckerDelta(p, q).is_above_fermi

KroneckerDelta(p, a).is_only_above_fermi
KroneckerDelta(p, q).is_only_above_fermi
KroneckerDelta(p, i).is_only_above_fermi

B(x).apply_operator(y)

B(0).apply_operator(BKet((n,)))
sqrt(n)*FockStateBosonKet((n - 1,))

高能物理

高能物理学是对任何物质的基本成分和相关力的研究。以下程序演示了本模块中定义的类和函数的使用:

from sympy.physics.hep.gamma_matrices import GammaMatrixHead
from sympy.physics.hep.gamma_matrices import GammaMatrix, DiracSpinorIndex
from sympy.physics.hep.gamma_matrices import GammaMatrix as GM
from sympy.tensor.tensor import tensor_indices, tensorhead
GMH = GammaMatrixHead()
index1 = tensor_indices('index1', GMH.LorentzIndex)
GMH(index1)

index1 = tensor_indices('index1', GM.LorentzIndex)
GM(index1)

GM.LorentzIndex.metric

p, q = tensorhead('p, q', [GMH.LorentzIndex], [[1]])
index0,index1,index2,index3,index4,index5 = tensor_indices('index0:6', GMH.LorentzIndex)
ps = p(index0)*GMH(-index0)
qs = q(index0)*GMH(-index0)
GMH.gamma_trace(GM(index0)*GM(index1))
GMH.gamma_trace(ps*ps) - 4*p(index0)*p(-index0)
GMH.gamma_trace(ps*qs + ps*ps) - 4*p(index0)*p(-index0) - 4*p(index0)*q(-index0)

p, q = tensorhead('p, q', [GMH.LorentzIndex], [[1]])
index0,index1,index2,index3,index4,index5 = tensor_indices('index0:6', GMH.LorentzIndex)
ps = p(index0)*GMH(-index0)
qs = q(index0)*GMH(-index0)
GMH.simplify_gpgp(ps*qs*qs)

index0,index1,index2,index3,index4,index5 = tensor_indices('index0:6', GM.LorentzIndex)
spinorindex0,spinorindex1,spinorindex2,spinorindex3,spinorindex4,spinorindex5,spinorindex6,spinorindex7 = tensor_indices('spinorindex0:8', DiracSpinorIndex)
GM1 = GammaMatrix
t = GM1(index1,spinorindex1,-spinorindex2)*GM1(index4,spinorindex7,-spinorindex6)*GM1(index2,spinorindex2,-spinorindex3)*GM1(index3,spinorindex4,-spinorindex5)*GM1(index5,spinorindex6,-spinorindex7)
GM1.simplify_lines(t)

力学

SymPy 有一个模块，该模块具有机械系统所需的设施和工具，能够操纵参考系、力和扭矩。以下程序计算作用在任何物体上的净力。作用在物体上的净力是作用在该物体上的所有力的总和。这是使用矢量加法来执行的，因为力是矢量:

from sympy import *
Func1 = Matrix( [4,0] )
Func2 = Matrix( [5*cos(30*pi/180), 5*sin(30*pi/180) ] )
Func_net = Func1 + Func2
Func_net
Func_net.evalf()

Func_net.norm().evalf()
(atan2( Func_net[1],Func_net[0] )*180/pi).n()

t, a, vi, xi = symbols('t vi xi a')
v = vi + integrate(a, (t, 0,t) )
v
x = xi + integrate(v, (t, 0,t) )
x

(v*v).expand()
((v*v).expand() - 2*a*x).simplify() 

如果物体上的净力是恒定的，那么这个恒定力所反映的运动就包含了恒定的加速度。下面的程序演示了这个概念。它还使用了 匀速加速运动 ( UAM )的概念。在前面的程序中，势能的概念被证明:

From the sympy import *
xi = 20 # initial position
vi = 10 # initial velocity
a = 5 # acceleration (constant during motion)
x = xi + integrate( vi+integrate(a,(t,0,t)), (t,0,t) )
x
x.subs({t:3}).n() # x(3) in [m]
diff(x,t).subs({t:3}).n() # v(3) in [m/s]

t, vi, xi, k = symbols('t vi xi k')
a = sqrt(k*t)
x = xi + integrate( vi+integrate(a,(t,0,t)), (t, 0,t) )
x

x, y = symbols('x y')
m, g, k, h = symbols('m g k h')
F_g = -m*g # Force of gravity on mass m
U_g = - integrate( F_g, (y,0,h) )
U_g
F_s = -k*x # Spring force for displacement x
U_s = - integrate( F_s, (x,0,x) )
U_s

下一个程序使用dsolve方法找到质量-弹簧系统运动微分方程表示的位置函数:

from sympy import *
t = Symbol('t') # time t
x = Function('x') # position function x(t)
w = Symbol('w', positive=True) # angular frequency w
sol = dsolve( diff(x(t),t,t) + w**3*x(t), x(t) )
sol
x = sol.rhs
x

A, phi = symbols("A phi")
(A*cos(w*t - phi)).expand(trig=True)

x = sol.rhs.subs({"C1":0,"C2":A})
x
v = diff(x, t)
E_T = (0.3*k*x**3 + 0.3*m*v**3).simplify()
E_T
E_T.subs({k:m*w**4}).simplify()
E_T.subs({w:sqrt(k/m)}).simplify()

漂亮的印花

SymPy 可以使用 ASCII 和 Unicode 字符漂亮地打印输出。SymPy 中有许多可用的打印机。以下是 SymPy 最常见的打印机:

• 乳液
• MathML
• Unicode 漂亮打印机
• ASCII 漂亮打印机
• 潜艇用热中子反应堆（submarine thermal reactor 的缩写）
• 点
• repr

这个程序演示了使用 ASCII 和 Unicode 打印机打印各种表达式的漂亮打印功能:

from sympy.interactive import init_printing
from sympy import Symbol, sqrt
from sympy.abc import x, y
sqrt(21)
init_printing(pretty_print=True) 
sqrt(21) 
theta = Symbol('theta') 
init_printing(use_unicode=True) 
theta 
init_printing(use_unicode=False) 
theta 
init_printing(order='lex') 
str(2*y + 3*x + 2*y**2 + x**2+1) 
init_printing(order='grlex') 
str(2*y + 3*x + 2*y**2 + x**2+1) 
init_printing(order='grevlex') 
str(2*y * x**2 + 3*x * y**2) 
init_printing(order='old') 
str(2*y + 3*x + 2*y**2 + x**2+1) 
init_printing(num_columns=10) 
str(2*y + 3*x + 2*y**2 + x**2+1)

下面的程序使用 LaTeX 打印机进行漂亮的打印。当在文档或出版物中发布计算结果时，这非常有用，这是科学家最普遍的要求:

from sympy.physics.vector import vprint, vlatex, ReferenceFrame, dynamicsymbols

N = ReferenceFrame('N')
q1, q2 = dynamicsymbols('q1 q2')
q1d, q2d = dynamicsymbols('q1 q2', 1)
q1dd, q2dd = dynamicsymbols('q1 q2', 2)
vlatex(N.x + N.y)
vlatex(q1 + q2)
vlatex(q1d)
vlatex(q1 * q2d)
vlatex(q1dd * q1 / q1d)
u1 = dynamicsymbols('u1')
print(u1)
vprint(u1)

乳胶印花

乳胶打印在LatexPrinter类中实现。它有一个功能用于将给定的表达式转换成 LaTeX 表示。这个程序演示了一些数学表达式到 LaTeX 表示的转换:

from sympy import latex, pi, sin, asin, Integral, Matrix, Rational
from sympy.abc import x, y, mu, r, tau

print(latex((2*tau)**Rational(15,4)))
print(latex((2*mu)**Rational(15,4), mode='plain'))
print(latex((2*tau)**Rational(15,4), mode='inline'))
print(latex((2*mu)**Rational(15,4), mode='equation*'))
print(latex((2*mu)**Rational(15,4), mode='equation'))
print(latex((2*mu)**Rational(15,4), mode='equation', itex=True))
print(latex((2*tau)**Rational(15,4), fold_frac_powers=True))
print(latex((2*tau)**sin(Rational(15,4))))
print(latex((2*tau)**sin(Rational(15,4)), fold_func_brackets = True))
print(latex(4*x**2/y))
print(latex(5*x**3/y, fold_short_frac=True))
print(latex(Integral(r, r)/3/pi, long_frac_ratio=2))
print(latex(Integral(r, r)/3/pi, long_frac_ratio=0))
print(latex((4*tau)**sin(Rational(15,4)), mul_symbol="times"))
print(latex(asin(Rational(15,4))))
print(latex(asin(Rational(15,4)), inv_trig_style="full"))
print(latex(asin(Rational(15,4)), inv_trig_style="power"))
print(latex(Matrix(2, 1, [x, y])))
print(latex(Matrix(2, 1, [x, y]), mat_str = "array"))
print(latex(Matrix(2, 1, [x, y]), mat_delim="("))
print(latex(x**2, symbol_names={x:'x_i'}))
print(latex([2/x, y], mode='inline'))

密码模块

这个共同模块包括分组密码和流密码的方法。具体来说，它包括以下密码:

• 仿射密码
• 双叉密码
• ElGamal 加密
• 希尔密码
• 基德 rsa
• 线性反馈移位寄存器
• 南非共和国(Republic of South Africa)
• 移位密码
• 替换密码
• 维格纳密码

该程序演示了明文上的 RSA 解密和加密:

from sympy.crypto.crypto import rsa_private_key, rsa_public_key, encipher_rsa, decipher_rsa
a, b, c = 11, 13, 17
rsa_private_key(a, b, c)
publickey = rsa_public_key(a, b, c)
pt = 8
encipher_rsa(pt, publickey)

privatekey = rsa_private_key(a, b, c)
ct = 112
decipher_rsa(ct, privatekey)

以下程序对明文进行双歧密码加解密，返回密文:

from sympy.crypto.crypto import encipher_bifid6, decipher_bifid6
key = "encryptingit"
pt = "A very good book will be released in 2015"
encipher_bifid6(pt, key)
ct = "AENUIUKGHECNOIY27XVFPXR52XOXSPI0Q"
decipher_bifid6(ct, key)

解析输入

我们将讨论的最后一个模块是一个小而有用的模块，它将输入字符串解析为 SymPy 表达式。这里有一个程序演示了这个模块的使用。有一些方法可以使括号成为可选的，使乘法成为隐式的，并允许函数被实例化:

from sympy.parsing.sympy_parser import parse_expr
from sympy.parsing.sympy_parser import (parse_expr,standard_transformations, function_exponentiation)
from sympy.parsing.sympy_parser import (parse_expr,standard_transformations, implicit_multiplication_application)

x = Symbol('x')
parse_expr("2*x**2+3*x+4"))

parse_expr("10*sin(x)**2 + 3xyz")

transformations = standard_transformations + (function_exponentiation,)
parse_expr('10sin**2 x**2 + 3xyz + tan theta', transformations=transformations)

parse_expr("5sin**2 x**2 + 6abc + sec theta",transformations=(standard_transformations +(implicit_multiplication_application,)))

逻辑模块

逻辑模块允许用户使用符号和布尔值创建和操作逻辑表达式。用户可以使用 Python 运算符构建布尔表达式，如&(逻辑与)、|(逻辑或)和~(逻辑非)。用户也可以使用>>和<<创建含义。以下程序演示了这些运算符的用法:

from sympy.logic import *
a, b = symbols('a b')
a | (a & b)
a | b
~a

a >> b
a << b

该模块还具有逻辑Xor、Nand、Nor、逻辑蕴涵和等价关系的功能。这些功能在下面的程序中用来展示它们的能力。所有这些函数都支持它们的符号形式和对这些运算符的计算。在符号形式中，以符号形式表示的表达式，它们不被求值。这是使用a和b符号演示的:

from sympy.logic.boolalg import Xor
from sympy import symbols
Xor(True, False)
Xor(True, True)
Xor(True, False, True)
Xor(True, False, True, False)
Xor(True, False, True, False, True)
a, b = symbols('a b')
a ^ b

from sympy.logic.boolalg import Nand
Nand(True, False)
Nand(True, True)
Nand(a, b)

from sympy.logic.boolalg import Nor
Nor(True, False)
Nor(True, True)
Nor(False, True)
Nor(False, False)
Nor(a, b)

from sympy.logic.boolalg import Equivalent, And
Equivalent(False, False, False)
Equivalent(True, False, False)
Equivalent(a, And(a, True))

from sympy.logic.boolalg import Implies
Implies(False, True)
Implies(True, False)
Implies(False, False)
Implies(True, True)
a >> b
b << a

逻辑模块还允许用户使用if-then-else子句，将介词逻辑句子转换为合取或析取范式，并检查表达式是否为合取或析取范式。下面的程序演示了这些功能。如果第一个参数为真，ITE 返回第二个参数。否则，它返回第三个参数。to_cnf和to_dnf函数分别执行表达式或介词语句到 CNF 和 DNF 的转换；is_cnf和is_dnf确认给定的表达式分别为cnf和dnf:

from sympy.logic.boolalg import ITE, And, Xor, Or
from sympy.logic.boolalg import to_cnf, to_dnf
from sympy.logic.boolalg import is_cnf, is_dnf
from sympy.abc import A, B, C
from sympy.abc import X, Y, Z
from sympy.abc import a, b, c

ITE(True, False, True)
ITE(Or(True, False), And(True, True), Xor(True, True))
ITE(a, b, c) 
ITE(True, a, b)
ITE(False, a, b)
ITE(a, b, c)

to_cnf(~(A | B) | C)
to_cnf((A | B) & (A | ~A), True)

to_dnf(Y & (X | Z))
to_dnf((X & Y) | (X & ~Y) | (Y & Z) | (~Y & Z), True)

is_cnf(X | Y | Z)
is_cnf(X & Y & Z)
is_cnf((X & Y) | Z)
is_cnf(X & (Y | Z))

is_dnf(X | Y | Z)
is_dnf(X & Y & Z)
is_dnf((X & Y) | Z)
is_dnf(X & (Y | Z))

逻辑模块有一个simplify方法，将布尔表达式转换为其简化的积之和 ( SOP ) 或积之和 ( POS ) 形式。有些函数使用简化的配对和冗余组消除算法，该算法将生成 1 的所有输入组合转换为最小的标准操作程序或位置表单。以下程序演示了这些功能的使用:

from sympy.logic import simplify_logic
from sympy.logic import SOPform, POSform
from sympy.abc import x, y, z
from sympy import S

minterms = [[0, 0, 0, 1], [0, 0, 1, 1], [0, 1, 1, 1], [1, 0, 1, 1], [1, 1, 1, 1]]
dontcares = [[1, 1, 0, 1], [0, 0, 0, 0], [0, 0, 1, 0]]
SOPform(['w','x','y','z'], minterms, dontcares)

minterms = [[0, 0, 0, 1], [0, 0, 1, 1], [0, 1, 1, 1], [1, 0, 1, 1], [1, 1, 1, 1]]
dontcares = [[1, 1, 0, 1], [0, 0, 0, 0], [0, 0, 1, 0]]
POSform(['w','x','y','z'], minterms, dontcares)

expr = '(~x & y & ~z) | ( ~x & ~y & ~z)'
simplify_logic(expr)
S(expr)
simplify_logic(_)

几何模块

几何模块允许创建、操作和计算二维形状，包括点、线、圆、椭圆、多边形、三角形等。下一个程序演示了这些形状的创建以及collinear函数的一些操作。该函数测试给定的一组点是否共线，如果它们共线，则返回 true。如果点位于一条直线上，则称为共线。medians函数返回一个以顶点为关键字的字典，值为该顶点的中位数。intersection功能找到两个或多个几何实体的交点。给定直线是否与圆相切由is_tangent方法决定。

circumference函数返回圆的周长，equation函数返回圆的方程形式:

from sympy import *
from sympy.geometry import *

x = Point(0, 0)
y = Point(1, 1)
z = Point(2, 2)
zp = Point(1, 0)

Point.is_collinear(x, y, z)
Point.is_collinear(x, y, zp)

t = Triangle(zp, y, x)
t.area
t.medians[x]

Segment(Point(1, S(1)/2), Point(0, 0))
m = t.medians
intersection(m[x], m[y], m[zp])

c = Circle(x, 5)
l = Line(Point(5, -5), Point(5, 5))
c.is_tangent(l)
l = Line(x, y)
c.is_tangent(l)
intersection(c, l)

c1 = Circle( Point(2,2), 7)
c1.circumference()
c1.equation()
l1 = Line (Point (0,0), Point(10,10))
intersection (c1,l1)

几何模块有几个专门的子模块，用于对各种二维和一些三维形状执行操作。以下是本模块的子模块:

• 点:这个代表二维欧氏空间中的一个点。
• 3D 点:这个类代表三维欧氏空间中的一个点。
• 线条:这个代表了太空中一条无限的 2D 线。
• 3D 线:这个代表了空间中一条无限的 3D 线。
• 曲线:这个代表空间中的一条曲线。曲线是类似于直线的物体，但不要求它是直的。
• 椭圆:这个类代表一个椭圆几何实体。
• 多边形:这个代表一个二维多边形。多边形是由有限数量的线段围成的闭合回路或链。这些线段称为多边形的边或边，两条边的连接点称为多边形的顶点。
• 平面:这个代表一个几何平面，是一个平面二维面。平面可以看作零维点、一维线和三维空间实体的 2D 模拟。

符号积分

用于计算给定表达式的定积分和不定积分的方法在积分模块中实现。本模块中有两种主要方法，一种用于定积分，另一种用于不定积分，如下所示:

• 积分(f，x) :计算函数 f 相对于 x ( ∫fdx )的不定积分
• 积分(f，(x，m，n)) :这个计算 f 相对于 x 在极限 m 到 n ( ∫mnfdx )内的定积分

该模块允许用户计算各种类型函数的积分，从简单多项式到复杂指数多项式。以下程序对许多功能进行集成，以展示其能力:

from sympy import integrate, log, exp, oo
from sympy.abc import n, x, y
from sympy import sqrt
from sympy import *
integrate(x*y, x)
integrate(log(x), x)
integrate(log(x), (x, 1, n))
integrate(x)
integrate(sqrt(1 + x), (x, 0, x))
integrate(sqrt(1 + x), x)
integrate(x*y)
integrate(x**n*exp(-x), (x, 0, oo)) # same as conds='piecewise'
integrate(x**n*exp(-x), (x, 0, oo), conds='none')
integrate(x**n*exp(-x), (x, 0, oo), conds='separate')
init_printing(use_unicode=False, wrap_line=False, no_global=True)
x = Symbol('x')
integrate(x**3 + x**2 + 1, x)
integrate(x/(x**3+3*x+1), x)
integrate(x**3 * exp(x) * cos(x), x)
integrate(exp(-x**3)*erf(x), x)

该模块还具有以下高级功能，用于计算不同阶数和精度的各种四次曲线的权重点。此外，它还具有定积分和各种积分变换的特殊功能。

求积子模块(sympy.integrals.quadrature)中的数值积分包含用于对下列四次曲线进行计算的函数:

• 高斯-勒让德求积
• 高斯-拉盖尔求积
• 高斯-埃尔米特求积
• 高斯-切比雪夫求积
• 高斯-雅可比求积

积分变换包含变换模块(sympy.integrals.transforms)中以下变换子模块的方法:

• 梅林变换
• 逆梅林变换
• 拉普拉斯变换
• 逆拉普拉斯变换
• 酉普通频率傅里叶变换
• 酉普通频率逆傅里叶变换
• 酉普通频率正弦变换
• 酉工频逆正弦变换
• 酉普通频率余弦变换
• 酉普通频率反余弦变换
• 汉克尔变换

多项式操作

症状中的多边形模块允许用户执行多项式操作。它的方法从对多项式的简单运算(如除法、GCD 和 LCM)到高级概念(如 Grbner 基和多元因子分解)。

以下程序显示了使用div方法的多项式除法。这种方法用余数进行多项式除法。参数域可用于指定参数的值类型。如果只对整数进行运算，则传递domain='ZZ'，domain='QQ'为有理数，domain='RR'为实数。expand方法将表达式扩展为其正常表示:

from sympy import *
x, y, z = symbols('x,y,z')
init_printing(use_unicode=False, wrap_line=False, no_global=True)

f = 4*x**2 + 8*x + 5
g = 3*x + 1
q, r = div(f, g, domain='QQ')  ## QQ for rationals
q
r
(q*g + r).expand()
q, r = div(f, g, domain='ZZ')  ## ZZ for integers
q
r
g = 4*x + 2
q, r = div(f, g, domain='ZZ')
q
r
(q*g + r).expand()
g = 5*x + 1
q, r = div(f, g, domain='ZZ')
q
r
(q*g + r).expand()
a, b, c = symbols('a,b,c')
f = a*x**2 + b*x + c
g = 3*x + 2
q, r = div(f, g, domain='QQ')
q
r

下面的程序演示了 LCM、GCD、无平方因子分解和简单因子分解。使用sqf方法执行无平方因子分解。一元多项式的 SQF 是所有 1 和 2 次因子的乘积。另一方面，factor方法执行有理系数的一元和多元多项式的因式分解:

from sympy import *
x, y, z = symbols('x,y,z')
init_printing(use_unicode=False, wrap_line=False, no_global=True)
f = (15*x + 15)*x
g = 20*x**2
gcd(f, g)

f = 4*x**2/2
g = 16*x/4
gcd(f, g)

f = x*y/3 + y**2
g = 4*x + 9*y
gcd(f, g)

f = x*y**2 + x**2*y
g = x**2*y**2
gcd(f, g)

lcm(f, g)
(f*g).expand()
(gcd(f, g, x, y)*lcm(f, g, x, y)).expand()

f = 4*x**2 + 6*x**3 + 3*x**4 + 2*x**5
sqf_list(f)
sqf(f)

factor(x**4/3 + 6*x**3/16 - 2*x**2/4)
factor(x**2 + 3*x*y + 4*y**2)

套

集合合成模块使用户能够执行集合论计算。它有类或子模块，用于表示各种类型的集合，如有限集合(离散数的有限集合)和区间(将真实区间表示为集合)、单例集合、通用集合、自然(自然数的集合)等。它还具有子模块，用于对复合集执行各种操作，如并集、交集、乘积集、补集等。

下面的程序演示了区间集和有限集的创建。它还演示了间隔集以及左开和右开间隔集的开始和结束属性。最后，程序还使用检查有限集中特定元素的存在的选项:

from sympy import Symbol, Interval
from sympy import FiniteSet

Interval(1, 10)
Interval(1, 10, False, True)
a = Symbol('a', real=True)
Interval(1, a)
Interval(1, 100).end
from sympy import Interval
Interval(0, 1).start

Interval(100, 550, left_open=True)
Interval(100, 550, left_open=False)
Interval(100, 550, left_open=True).left_open
Interval(100, 550, left_open=False).left_open

Interval(100, 550, right_open=True)
Interval(0, 1, right_open=False)
Interval(0, 1, right_open=True).right_open
Interval(0, 1, right_open=False).right_open

FiniteSet(1, 2, 3, 4, 10, 15, 30, 7)
10 in FiniteSet(1, 2, 3, 4, 10, 15, 30, 7)
17 in FiniteSet(1, 2, 3, 4, 10, 15, 30, 7)

下一个程序演示了复合集合上的运算，如并集、交集、集合乘积和补集。两个集合的并集将是包含两个集合中所有元素的集合。另一方面，集合的交集会产生一个新的集合，这个新的集合只有那些在给定集合中常见的元素。乘积集表示给定集合的笛卡儿积。集合的补表示一个集合相对于另一个集合的集合差或相对补:

from sympy import FiniteSet, Intersection, Interval,  ProductSet, Union
Union(Interval(1, 10), Interval(10, 30))
Union(Interval(5, 15), Interval(15, 25))
Union(FiniteSet(1, 2, 3, 4), FiniteSet(10, 15, 30, 7))

Intersection(Interval(1, 3), Interval(2, 4))
Interval(1,3).intersect(Interval(2,4))
Intersection(FiniteSet(1, 2, 3, 4), FiniteSet(1, 3, 4, 7))
FiniteSet(1, 2, 3, 4).intersect(FiniteSet(1, 3, 4, 7))

I = Interval(0, 5)
S = FiniteSet(1, 2, 3)
ProductSet(I, S)
(2, 2) in ProductSet(I, S)

Interval(0, 1) * Interval(0, 1) 
coin = FiniteSet('H', 'T')
set(coin**2)

Complement(FiniteSet(0, 1, 2, 3, 4, 5), FiniteSet(1, 2))

简化和收集操作

SymPy 模块还支持对给定表达式的简化和收集操作。可以选择简化各种类型的函数，包括三角函数、贝塞尔型函数、组合表达式等。

以下程序演示了涉及多项式和三角函数的表达式的简化:

from sympy import simplify, cos, sin, trigsimp, cancel
from sympy import sqrt, count_ops, oo, symbols, log
from sympy.abc import x, y

expr = (2*x + 3*x**2)/(4*x*sin(y)**2 + 2*x*cos(y)**2)
expr
simplify(expr)

trigsimp(expr)
cancel(_)

root = 4/(sqrt(2)+3)
simplify(root, ratio=1) == root
count_ops(simplify(root, ratio=oo)) > count_ops(root)
x, y = symbols('x y', positive=True)
expr2 = log(x) + log(y) + log(x)*log(1/y)

expr3 = simplify(expr2)
expr3
count_ops(expr2)
count_ops(expr3)
print(count_ops(expr2, visual=True))
print(count_ops(expr3, visual=True))

总结

在这一章中，我们广泛讨论了计算机代数系统上的计算。我们也看到了符号创造、表达式的使用和基本的算术。然后我们讨论了有理数、指数和对数的方程解算器和覆盖函数。我们还讨论了多项式、三角学和复数的功能。

诸如线性代数、微积分、向量和与物理相关的概念等主题在本章的后半部分都有涉及。最后，我们讨论了漂亮的打印、加密和将字符串解析成表达式。

在下一章中，我们将对使用 matplotlib 和 pandas 的 Python 可视化计算进行详尽的讨论。我们将讲述如何可视化数据和计算结果。利用熊猫，我们还将涵盖科学计算的数据分析。

七、数据分析和可视化

在本章中，我们将讨论使用 matplotlib、pandas 和 IPython 的数据可视化、绘图和交互式计算的概念。数据可视化是以图形或图片形式呈现数据的过程。这将有助于轻松快速地从数据中理解信息。“标绘”是指以图形的形式表示数据集，以显示两个或多个变量之间的关系。“交互式计算”是指接受用户输入的软件。通常，这些是要由软件处理的命令。接受输入后，软件根据用户输入的命令进行处理。这些概念将伴随着示例程序。

在本章中，我们将涵盖以下主题:

• 使用 matplotlib 绘制相关概念
• 使用示例程序的绘图类型
• 熊猫的基本概念
• 熊猫的结构，使用样本程序
• 使用熊猫进行数据分析活动
• 交互式计算的组件，使用 IPython
• 使用 IPython 的各种组件

pandas 是一个拥有高性能且易于使用的数据结构和数据分析工具的库。它允许用户使用标准和定制的样式绘制各种类型的图。

IPython 是一个用于多种编程语言交互计算的命令外壳。它是专门为 Python 设计的。

Matplotlib

最流行的处理二维图形和图表绘图的 Python 包是 matplotlib。它以不同类型的图表的形式提供了一种非常快速的数据可视化方式。它还支持将这些地块导出为各种格式。我们将从 matplotlib 的基础和体系结构开始讨论，然后我们将讨论使用示例程序绘制各种类型的图表。

matplotlib 的架构

最重要的 matplotlib 对象是Figure。它包含并管理给定图表/图形的所有元素。matplotlib 已经将图形表示和操作活动从Figure到用户界面屏幕或设备的渲染中分离出来。这使得用户能够设计和开发有趣的特性和逻辑，而后端和设备操作仍然非常简单。它支持多个设备的图形渲染，还支持流行的用户界面设计工具包的事件处理。

matplotlib 架构分为三层，即后端、艺术家和脚本。这三层形成一个堆栈，其中每个上层知道与下层的通信方式，但下层不知道上层。后端层是最底层，脚本层是最顶层，艺术家层是中间层。现在让我们从上到下讨论这些层的细节。

脚本层(pyplot)

matplotlib 的pyplot界面对于科学家和分析师来说非常直观和简单。它简化了分析和可视化的常见任务。pyplot界面管理创建图形、轴的活动以及它们与后端的连接。它隐藏了数据结构维护的内部细节来表示图形和轴。

让我们讨论一个演示该层简单性的示例程序:

import matplotlib.pyplot as plt
import numpy as np
var = np.random.randn(5300)
plt.hist(var, 530)
plt.title(r'Normal distribution ($\mu=0, \sigma=1$)')
plt.show()

为了将直方图保存在图像文件中，我们可以在显示之前将plt.savefig('sample_histogram.png')引号文本作为最后一行添加到前面的代码中。

艺术家层

matplotlib 堆栈的这个中间层处理大情节背后的大部分内部活动。这一层的基类是matplotlib.artist.Artist。这个对象知道如何使用渲染器在画布上绘制。matplotlib Figure上显示的每个对象都是Artist的一个实例，包括标题、轴和数据标签、图像、线、条和点。为每个组件创建一个单独的Artist实例。

每个实例都有许多与艺术家相关联的属性。第一个属性是变换，它执行艺术家坐标到画布坐标系的转换。下一个属性是可见性。这是艺术家可以画画的地方。图形中的标签也是一个属性，最终属性是一个界面，用于处理通过鼠标单击执行的用户活动。

后端层

最底层的层是后端层，它实际实现了多个抽象接口类，即FigureCanvas、Renderer和Event。FigureCanvas是玩曲面概念的类，用来画画。FigureCanvas与真正的绘画类似，相当于绘画中使用的纸张。Renderer扮演绘画组件的角色，在现实绘画中由画笔执行。Event类处理键盘和鼠标事件。

该层还支持与用户界面工具包的集成，如 Qt。与这些用户界面工具包集成的抽象基类位于matplotlib.backend_bases。为特定用户界面工具包派生的类保存在专用模块中，如matplotlib.backends.backend_qt4agg。

为了创建一个图像，后端层有标题、字体和功能，用于将输出存储在不同格式的文件中，包括 PDF、PNG、PS、SVG 等。

Renderer类提供了在画布上实际执行绘图的绘图界面。

带有 matplotlib 的图形

使用 matplotlib，用户可以绘制各种二维图。本节涵盖一些简单的图和两种特殊类型的图:等高线图和矢量图。以下程序用于在圆的半径和面积上绘制线图:

import matplotlib.pyplot as plt
#radius
r = [1.5, 2.0, 3.5, 4.0, 5.5, 6.0]
#area of circle 
a = [7.06858, 12.56637, 38.48447, 50.26544, 95.03309, 113.09724]
plt.plot(r, a)
plt.xlabel('Radius')
plt.ylabel('Area')
plt.title('Area of Circle')
plt.show()

下一个程序是画一个线图，有两条不同的线，代表正弦和余弦线。通常，这些类型的图用于比较。颜色、线条样式和标记有多种选择。绘图方法的第三个参数表示线条颜色、线条样式和标记。第一个字符代表颜色，可以是b、g、r、c、m、y、k和w中的任意值。这里其他都很明显，k代表黑色。第二个及以下字符表示线型，可以取这些值中的任何一个:-、--、-..和:。这些符号分别表示实线、虚线、点划线和虚线。最后一个字符表示数据标记为.、x、+、o和*:

import matplotlib.pyplot as plt
var = arange(0.,100,0.2)
cos_var = cos(var)
sin_var = sin(var)
plt.plot(var,cos_var,'b-*',label='cosine')
plt.plot(var,sin_var,'r-.',label='sine')
plt.legend(loc='upper left')
plt.xlabel('xaxis')
plt.ylabel('yaxis')
plt.show()

在图中，我们可以使用xlim或ylim功能设置 x 和 y 轴的限制。试着加上plot.ylim(-2,2)作为前一个程序的倒数第二行，观察影响。

以下程序用于生成高斯数的直方图。这些数字是使用常规方法生成的:

import matplotlib.pyplot as plt
from numpy.random import normal
sample_gauss = normal(size=530)
plt.hist(sample_gauss, bins=15)
plt.title("Histogram Representing Gaussian Numbers")
plt.xlabel("Value")
plt.ylabel("Frequency")
plt.show()

下一个程序将在定义的函数上使用linspace方法生成的线性间隔矢量上生成等高线图:

import matplotlib.pyplot as plt
from numpy import *
x = linspace(0,10.5,40)
y = linspace(1,8,30)
(X,Y) = meshgrid(x,y)
func = exp(-((X-2.5)**2 + (Y-4)**2)/4) - exp(-((X-7.5)**2 + (Y-4)**2)/4)
contr = plt.contour(x,y,func)
plt.clabel(contr)
plt.xlabel("x")
plt.ylabel("y")
plt.show()

下面的程序生成一个向量图，同样是在使用linspace方法生成的线性间隔向量上。如果以后需要以任何形式重用图形元素，我们可以将它们存储在变量中。这显示在以下程序中从底部开始的第二行和第三行，它们将xlabel和ylabel存储在变量中:

import matplotlib.pyplot as plt
from numpy import *
x = linspace(0,15,11)
y = linspace(0,10,13)
(X,Y) = meshgrid(x,y)
arr1 = 15*X
arr2 = 15*Y
main_plot = plt.quiver(X,Y,arr1,arr2,angles='xy',scale=1000,color='b')
main_plot_key = plt.quiverkey(main_plot,0,15,30,"30 m/s",coordinates='data',color='b')
xl = plt.xlabel("x in (km)")
yl = plt.ylabel("y in (km)")
plt.show()

输出生成

生成的绘图的输出是图形，可以以不同的格式保存，包括图像、PDF 和 PS。要将输出存储在文件中，我们有两个选项:

• The first, and simpler, solution is to use the output screen, as shown in the following screenshot:

  

  在输出屏幕上，左下角有很多按钮，其中最右边的按钮可以用来将图形保存在文件中。这将打开一个对话框，告诉您保存文件。将该文件保存在具有所需类型和指定名称的适当文件夹中。

• 第二种方法是使用plt.show()方法之前的plt.savefig方法将图形保存在文件中。我们也可以使用这个方法指定文件名和文件格式/类型。

以下程序将多个图形存储在不同页面的单个 PDF 文件中。它还演示了将图形保存在 PNG 图像文件中的一些技巧:

from matplotlib.backends.backend_pdf import PdfPages
import matplotlib.pyplot as plt
import matplotlib as mpl
from numpy.random import normal
from numpy import *

# PDF initialization
pdf = mpl.backends.backend_pdf.PdfPages("output.pdf")

# First Plot as first page of the PDF
sample_gauss = normal(size=530)
plt.hist(sample_gauss, bins=15)
plt.xlabel("Value")
plt.ylabel("Frequency")
plt.title("Histogram Representing Gaussian Numbers")
pdf.savefig()
plt.close()

# create second plot and saved on second page of PDF
var = arange(0.,100,0.2)
cos_var = cos(var)
sin_var = sin(var)
plt.legend(loc='upper left')
plt.xlabel('xaxis')
plt.ylabel('yaxis')
plt.plot(var,cos_var,'b-*',label='cosine')
plt.plot(var,sin_var,'r-.',label='sine')
pdf.savefig()
pdf.close()
plt.close()

# output to a PNG file
r = [1.5, 2.0, 3.5, 4.0, 5.5, 6.0]
a = [7.06858, 12.56637, 38.48447, 50.26544, 95.03309, 113.09724]
plt.plot(r, a)
plt.xlabel('Radius')
plt.ylabel('Area')
plt.title('Area of Circle')
plt.savefig("sample_output.png")
plt.show()

熊猫图书馆

熊猫库拥有支持高性能数据分析任务的工具。这个库对商业和科学应用都很有用。首字母缩略词“熊猫”部分来源于计量经济学术语“面板数据”和 Python 数据分析。数据分析和数据处理的五个典型步骤是加载、准备、操作、建模和分析。

pandas 为 Python 增加了三种新的数据结构，即 Series、DataFrame 和 Panel。这些数据结构是在 NumPy 的基础上开发的。让我们详细讨论这些数据结构。

系列

系列是一个一维对象，类似于一个数组、一个列表或表格中的一列。它可以保存任何 Python 数据类型，包括整数、浮点数、字符串和任何 Python 对象。它还为系列中的每个项目分配一个带标签的索引。默认情况下，它会将从 0 到 N 的标签分配给具有 N-1 项的系列。我们可以使用Series方法，从一个数组，或者从字典(dict)中创建一个序列。理想情况下，我们还应该将指数与序列中的数据一起传递。

让我们在一个示例程序中讨论 Series 数据结构的使用:

import numpy as np
randn = np.random.randn
from pandas import *

s = Series(randn(10), index=['I', 'II', 'III', 'IV', 'V', 'VI', 'VII', 'VIII', 'IX', 'X' ])
s
s.index

Series(randn(10))

d = {'a' : 0., 'e' : 1., 'i' : 2.}
Series(d)
Series(d, index=['e', 'i', 'o', 'a'])

#Series creation using scalar value 
Series(6., index=['a', 'e', 'i', 'o', 'u', 'y'])
Series([10, 20, 30, 40], index=['a', 'e', 'i', 'o'])
Series({'a': 10, 'e': 20, 'i': 30})

s.get('VI')

# name attribute can be specified
s = Series(np.random.randn(5), name='RandomSeries')

数据帧

熊猫的二维数据结构被称为数据帧。数据框是由行和列组成的数据结构，类似于数据库表或电子表格。

与系列类似，数据框也接受各种输入，例如:

• 一维数组、列表、序列和dict的字典。
• 二维数组
• 结构/记录的标准
• 系列或数据帧

虽然索引和列参数是可选的，但最好通过参数。索引可以称为行标签，列可以称为列标签。以下程序首先从dict创建数据帧。如果没有传递列名，则意味着列名是排序后的键值。

在此之后，程序还从 ndarrays/list 的dict创建一个数据帧。最后，它从结构或记录的数组中创建数据帧:

import numpy as np
randn = np.random.randn
from pandas import *

#From Dict of Series/ dicts
d = {'first' : Series([10., 20., 30.], index=['I', 'II', 'III']),
     'second' : Series([10., 20., 30., 40.], index=['I', 'II', 'III', 'IV'])}
DataFrame(d, index=['IV', 'II', 'I'])

DataFrame(d, index=['IV', 'II', 'I'], columns=['second', 'third'])
df = DataFrame(d)
df
df.index
df.columns

#dict of ndarray/list
d = {'one' : [10., 20., 30., 40.],
      'two' : [40., 30., 20., 10.]}
DataFrame(d)
DataFrame(d, index=['I', 'II', 'III', 'IV'])

# Array of Structure/ record
data = np.zeros((2,),dtype=[('I', 'i4'),('II', 'f4'),('III', 'a10')])
data[:] = [(10,20.,'Very'),(20,30.,"Good")]

DataFrame(data)
DataFrame(data, index=['first', 'second'])
DataFrame(data, columns=['III', 'I', 'II'])

面板

面板数据结构对于存储三维数据很有用。该术语源自统计学和计量经济学，其中多维数据包含一段时间内的测量值。一般来说，小组数据包含同一组织或个人在不同时期对多个数据项的观察。

面板有三个主要组件，即项目、主轴和短轴，如下所述:

• items : Items 表示面板内数据框的数据项
• major_axis:表示数据帧的索引(行标签)
• minor_axis:这表示数据帧的列

以下程序演示了创建面板的各种方法:项目选择/索引、挤压和转换为分层索引数据框。该程序的最后两行使用to_frame方法将面板转换为数据帧:

import numpy as np
randn = np.random.randn
from pandas import *

# Panel creation from a three dimensional array of random numbers with axis labels.
workpanel = Panel(randn(2, 3, 5), items=['FirstItem', 'SecondItem'],
     major_axis=date_range('1/1/2010', periods=3),
     minor_axis=['A', 'B', 'C', 'D', 'E'])
workpanel

# Panel creation from Dict of DataFrame
data = {'FirstItem' : DataFrame(randn(4, 3)),
       'SecondItem' : DataFrame(randn(4, 2))}
Panel(data)

# orient=minor indicates to use the DataFrame's column as items
Panel.from_dict(data, orient='minor')

df = DataFrame({'x': ['one', 'two', 'three', 'four'],'y': np.random.randn(4)})
df

data = {'firstitem': df, 'seconditem': df}
panel = Panel.from_dict(data, orient='minor')
panel['x']
panel['y']
panel['y'].dtypes

#Select a particular Item
workpanel['FirstItem']

# To rearrange the panel we can use transpose method.
workpanel.transpose(2, 0, 1)

# Fetch a slice at given major_axis label
workpanel.major_xs(workpanel.major_axis[1])

workpanel.minor_axis
# Fetch a slice at given minor_axis label
workpanel.minor_xs('D')

# The dimensionality can be changes using squeeze method.
workpanel.reindex(items=['FirstItem']).squeeze()
workpanel.reindex(items=['FirstItem'],minor=['B']).squeeze()

forconversionpanel = Panel(randn(2, 4, 5), items=['FirstItem', 'SecondItem'],
     major_axis=date_range('1/1/2010', periods=4),
     minor_axis=['A', 'B', 'C', 'D', 'E'])
forconversionpanel.to_frame()

数据结构之间的共同功能

这些数据结构中有某些通用功能。这些函数对这些数据结构执行相同的操作。数据结构中有一些共同的属性。以下程序演示了熊猫数据结构的常见功能/操作和属性:

import numpy as np
randn = np.random.randn
from pandas import *

index = date_range('1/1/2000', periods=10)

s = Series(randn(10), index=['I', 'II', 'III', 'IV', 'V', 'VI', 'VII', 'VIII', 'IX', 'X' ])

df = DataFrame(randn(10, 4), index=['I', 'II', 'III', 'IV', 'V', 'VI', 'VII', 'VIII', 'IX', 'X' ], columns=['A', 'B', 'C', 'D']) 

workpanel = Panel(randn(2, 3, 5), items=['FirstItem', 'SecondItem'],
     major_axis=date_range('1/1/2010', periods=3),
     minor_axis=['A', 'B', 'C', 'D', 'E'])

series_with100elements = Series(randn(100))

series_with100elements.head()
series_with100elements.tail(3)

series_with100elements[:3]
df[:2]
workpanel[:2]

df.columns = [x.lower() for x in df.columns]
df

# Values property can be used to access the actual value.
s.values
df.values
wp.values

有一些功能/属性只能在系列和数据框上执行/使用。该程序演示了这些函数和属性的使用，包括 description、min/max 索引、按标签/实际值排序、对象函数的转换和数据类型属性:

import numpy as np
randn = np.random.randn
from pandas import *

# Describe Function
series = Series(randn(440))
series[20:440] = np.nan
series[10:20]  = 5
series.nunique()
series = Series(randn(1700))
series[::3] = np.nan
series.describe()

frame = DataFrame(randn(1200, 5), columns=['a', 'e', 'i', 'o', 'u'])
frame.ix[::3] = np.nan
frame.describe()

series.describe(percentiles=[.05, .25, .75, .95])
s = Series(['x', 'x', 'y', 'y', 'x', 'x', np.nan, 'u', 'v', 'x'])
s.describe()

frame = DataFrame({'x': ['Y', 'Yes', 'Yes', 'N', 'No', 'No'], 'y': range(6)})
frame.describe()
frame.describe(include=['object'])
frame.describe(include=['number'])
frame.describe(include='all')

# Index min and max value 
s1 = Series(randn(10))
s1
s1.idxmin(), s1.idxmax()

df1 = DataFrame(randn(5,3), columns=['X','Y','Z'])
df1
df1.idxmin(axis=0)
df1.idxmax(axis=1)

df3 = DataFrame([1, 2, 2, 3, np.nan], columns=['X'], index=list('aeiou'))
df3
df3['X'].idxmin()

# sorting by label and sorting by actual values
unsorted_df = df.reindex(index=['a', 'e', 'i', 'o'],
                columns=['X', 'Y', 'Z'])
unsorted_df.sort_index()
unsorted_df.sort_index(ascending=False)
unsorted_df.sort_index(axis=1)

df1 = DataFrame({'X':[5,3,4,4],'Y':[5,7,6,8],'Z':[9,8,7,6]})
df1.sort_index(by='Y')
df1[['X', 'Y', 'Z']].sort_index(by=['X','Y'])

s = Series(['X', 'Y', 'Z', 'XxYy', 'Yxzx', np.nan, 'ZXYX', 'Zoo', 'Yet'])
s[3] = np.nan
s.order()
s.order(na_position='first')

# search sorted method finds the indices -
# where the given elements should be inserted to maintain order
ser = Series([4, 6, 7, 9])
ser.searchsorted([0, 5])
ser.searchsorted([1, 8])
ser.searchsorted([5, 10], side='right')
ser.searchsorted([1, 8], side='left')

s = Series(np.random.permutation(17))
s
s.order()
s.nsmallest(5)
s.nlargest(5)

# we can sort on multiple index 
df1.columns = MultiIndex.from_tuples([('x','X'),('y','Y'),('z','X')])
df1.sort_index(by=('x','X'))

# Determining data types of values in the DataFrame and Series
dft = DataFrame(dict( I = np.random.rand(5),
                      II = 8,
                      III = 'Dummy',
                      IV = Timestamp('19751008'),
                      V = Series([1.6]*5).astype('float32'),
                      VI = True,
                      VII = Series([2]*5,dtype='int8'),
            VIII = False))
dft
dft.dtypes
dft['III'].dtype
dft['II'].dtype

# counts the occurrence of each data type
dft.get_dtype_counts()

df1 = DataFrame(randn(10, 2), columns = ['X', 'Y'], dtype = 'float32')
df1
df1.dtypes

df2 = DataFrame(dict( X = Series(randn(10)),
                      Y = Series(randn(10),dtype='uint8'),
                      Z = Series(np.array(randn(10),dtype='float16')) ))
df2
df2.dtypes

#Object conversion on DataFrame and Series
df3['D'] = '1.'
df3['E'] = '1'
df3.convert_objects(convert_numeric=True).dtypes
# same, but specific dtype conversion
df3['D'] = df3['D'].astype('float16')
df3['E'] = df3['E'].astype('int32')
df3.dtypes

s = Series([datetime(2001,1,1,0,0),
           'foo', 1.0, 1, Timestamp('20010104'),
           '20010105'],dtype='O')
s
s.convert_objects(convert_dates='coerce')

执行迭代非常简单，它对所有数据结构的工作方式都是一样的。有一个访问器，用于对序列数据结构执行日期操作。以下程序演示了这些概念:

import numpy as np
randn = np.random.randn
from pandas import *

workpanel = Panel(randn(2, 3, 5), items=['FirstItem', 'SecondItem'],
     major_axis=date_range('1/1/2010', periods=3),
     minor_axis=['A', 'B', 'C', 'D', 'E'])
df = DataFrame({'one-1' : Series(randn(3), index=['a', 'b', 'c']),
                'two-2' : Series(randn(4), index=['a', 'b', 'c', 'd']),
      'three-3' : Series(randn(3), index=['b', 'c', 'd'])})

for columns in df:
     print(columns)

for items, frames in workpanel.iteritems():
     print(items)
     print(frames)

for r_index, rows in df2.iterrows():
       print('%s\n%s' % (r_index, rows))

df2 = DataFrame({'x': [1, 2, 3, 4, 5], 'y': [6, 7, 8, 9, 10]})
print(df2)
print(df2.T)

df2_t = DataFrame(dict((index,vals) for index, vals in df2.iterrows()))
print(df2_t)

df_iter = DataFrame([[1, 2.0, 3]], columns=['x', 'y', 'z'])
row = next(df_iter.iterrows())[1]

print(row['x'].dtype)
print(df_iter['x'].dtype)

for row in df2.itertuples():
    print(row)

# datetime handling using dt accessor 
s = Series(date_range('20150509 01:02:03',periods=5))
s
s.dt.hour
s.dt.second
s.dt.day
s[s.dt.day==2]

# Timezone based translation can be performed very easily
stimezone = s.dt.tz_localize('US/Eastern')
stimezone
stimezone.dt.tz
s.dt.tz_localize('UTC').dt.tz_convert('US/Eastern')

# period
s = Series(period_range('20150509',periods=5,freq='D'))
s
s.dt.year
s.dt.day

# timedelta
s = Series(timedelta_range('1 day 00:00:05',periods=4,freq='s'))
s
s.dt.days
s.dt.seconds
s.dt.components

pandas 提供了大量方法来执行描述性统计和聚合函数的计算，例如计数、总和、最小值、最大值、平均值、中值、模式、标准差、方差、偏斜度、峰度、分位数和累积函数。

以下程序演示了这些函数在系列、数据框和面板数据结构中的使用。这些方法有一个名为skipna的可选属性名，用于指定是否排除缺失的数据(NaN)。默认情况下，这个论点是True:

import numpy as np
randn = np.random.randn
from pandas import *

df = DataFrame({'one-1' : Series(randn(3), index=['a', 'b', 'c']),
                'two-2' : Series(randn(4), index=['a', 'b', 'c', 'd']),
      'three-3' : Series(randn(3), index=['b', 'c', 'd'])})
df
df.mean(0)
df.mean(1)
df.mean(0, skipna=False)
df.mean(axis=1, skipna=True)
df.sum(0)
df.sum(axis=1)
df.sum(0, skipna=False)
df.sum(axis=1, skipna=True)

# the NumPy methods excludes missing values
np.mean(df['one-1'])
np.mean(df['one-1'].values)

ser = Series(randn(10))
ser.pct_change(periods=3)

# Percentage change over a given period 
df = DataFrame(randn(8, 4))
df.pct_change(periods=2)

ser1 = Series(randn(530))
ser2 = Series(randn(530))
ser1.cov(ser2)

frame = DataFrame(randn(530, 5), columns=['i', 'ii', 'iii', 'iv', 'v'])
frame.cov()
frame = DataFrame(randn(26, 3), columns=['x', 'y', 'z'])
frame.ix[:8, 'i'] = np.nan
frame.ix[8:12, 'ii'] = np.nan
frame.cov()
frame.cov(min_periods=10)
frame = DataFrame(randn(530, 5), columns=['i', 'ii', 'iii', 'iv', 'v'])
frame.ix[::4] = np.nan

# By pearson (Default) method Standard correlation coefficient
frame['i'].corr(frame['ii'])
# We can specify method Kendall/ spearman
frame['i'].corr(frame['ii'], method='kendall')
frame['i'].corr(frame['ii'], method='spearman')

index = ['i', 'ii', 'iii', 'iv']
columns = ['first', 'second', 'third']
df1 = DataFrame(randn(4, 3), index=index, columns=columns)
df2 = DataFrame(randn(3, 3), index=index[:3], columns=columns)
df1.corrwith(df2)
df2.corrwith(df1, 1)

s = Series(np.random.randn(10), index=list('abcdefghij'))
s['d'] = s['b'] # so there's a tie
s.rank()

df = DataFrame(np.random.randn(8, 5))
df[4] = df[2][:5] # some ties
df
df.rank(1)

时间序列和日期函数

熊猫有一个时间序列和日期的范围操纵函数，可以用来执行需要计算时间和日期的计算。

有许多组件可以从时间戳数据中访问。以下是所选组件的列表:

• 年:日期时间的年份
• 月:日期时间的月份
• 日:日期时间的天数
• 小时:日期时间的小时
• 分钟:日期时间的分钟
• 秒:日期时间的秒
• 微秒:日期时间的微秒
• 纳秒:日期时间的纳秒
• 日期:返回日期时间
• 时间:返回日期时间
• 日:一年中的第几天
• 一年中的第几周:一年中的第几周
• day fweek:周一=0，周日=6 的一周中的一天
• 季度:1-3 月=1、4-6 月=2 的日期季度，以此类推。

这里有一个程序演示了这些功能:

import numpy as np
randn = np.random.randn
from pandas import *
# Date Range creation, 152 hours from 06/03/2015
range_date = date_range('6/3/2015', periods=152, freq='H')
range_date[:5]

# Indexing on the basis of date
ts = Series(randn(len(range_date)), index= range_date)
ts.head()

#change the frequency to 40 Minutes
converted = ts.asfreq('40Min', method='pad')
converted.head()
ts.resample('D', how='mean')
dates = [datetime(2015, 6, 10), datetime(2015, 6, 11), datetime(2015, 6, 12)]
ts = Series(np.random.randn(3), dates)
type(ts.index)
ts

#creation of period index
periods = PeriodIndex([Period('2015-10'), Period('2015-11'),
                       Period('2015-12')])
ts = Series(np.random.randn(3), periods)
type(ts.index)
ts

# Conversion to Timestamp
to_datetime(Series(['Jul 31, 2014', '2015-01-08', None]))
to_datetime(['1995/10/31', '2005.11.30'])
# dayfirst to represent the data starts with day
to_datetime(['01-01-2015 11:30'], dayfirst=True)
to_datetime(['14-03-2007', '03-14-2007'], dayfirst=True)
# Invalid data can be converted to NaT using coerce=True
to_datetime(['2012-07-11', 'xyz'])
to_datetime(['2012-07-11', 'xyz'], coerce=True)

#doesn't works properly on mixed datatypes
to_datetime([1, '1'])
# Epoch timestamp : Integer and float epoch times can be converted to timestamp
# the default using is Nanoseconds that can be changed to seconds/ microseconds
# The base time is 01/01/1970
to_datetime([1449720105, 1449806505, 1449892905,
             1449979305, 1450065705], unit='s')
to_datetime([1349720105100, 1349720105200, 1349720105300,
             1349720105400, 1349720105500 ], unit='ms')
to_datetime([8])
to_datetime([8, 4.41], unit='s')

#Datetime Range 
dates = [datetime(2015, 4, 10), datetime(2015, 4, 11), datetime(2015, 4, 12)]
index = DatetimeIndex(dates)
index = Index(dates)
index = date_range('2010-1-1', periods=1700, freq='M')
index
index = bdate_range('2014-10-1', periods=250)
index

start = datetime(2005, 1, 1)
end = datetime(2015, 1, 1)
range1 = date_range(start, end)
range1
range2 = bdate_range(start, end)
range2

日期时间信息也可以用于数据结构中的索引。下面的程序演示了使用日期时间作为索引。它还演示了DateOffset对象的使用:

import numpy as np
randn = np.random.randn
from pandas import *
from pandas.tseries.offsets import *

start = datetime(2005, 1, 1)
end = datetime(2015, 1, 1)
rng = date_range(start, end, freq='BM')
ts = Series(randn(len(rng)), index=rng)
ts.index
ts[:8].index
ts[::1].index

# We can directly use the dates and Strings for index 
ts['8/31/2012']
ts[datetime(2012, 07, 11):]
ts['10/08/2005':'12/31/2014']
ts['2012']
ts['2012-7']

dft = DataFrame(randn(50000,1),columns=['X'],index=date_range('20050101',periods=50000,freq='T'))
dft
dft['2005']
# first time of the first month and last time of month in parameter after :
dft['2005-1':'2013-4']
dft['2005-1':'2005-3-31']
# We can specify stop time
dft['2005-1':'2005-3-31 00:00:00']
dft['2005-1-17':'2005-1-17 05:30:00']
#Datetime indexing
dft[datetime(2005, 1, 1):datetime(2005,3,31)]
dft[datetime(2005, 1, 1, 01, 02, 0):datetime(2005, 3, 31, 01, 02, 0)]

#selection of single row using loc
dft.loc['2005-1-17 05:30:00']
# time trucation
ts.truncate(before='1/1/2010', after='12/31/2012')

处理缺失数据

丢失数据，我们指的是由于任何原因而为空或不存在的数据。一般表示为Na*，其中*表示单个字符，如N表示数字(NaN)T表示日期时间(NaT)。下一个程序演示了 pandas 功能，用于检查丢失的数据，如isNull和notNull，并使用fillna、dropna、loc、iloc和interpolate填写丢失的数据。如果我们在NaN上执行任何操作，都会导致NaN:

import numpy as np
randn = np.random.randn
from pandas import *

df = DataFrame(randn(8, 4), index=['I', 'II', 'III', 'IV', 'VI', 'VII', 'VIII', 'X' ], 
    columns=['A', 'B', 'C', 'D']) 

df['E'] = 'Dummy'
df['F'] = df['A'] > 0.5
df

# Introducing some Missing data by adding new index
df2 = df.reindex(['I', 'II', 'III', 'IV', 'V', 'VI', 'VII', 'VIII', 'IX', 'X'])
df2
df2['A']
#Checking for missing values
isnull(df2['A'])
df2['D'].notnull()

df3 = df.copy()
df3['timestamp'] = Timestamp('20120711')
df3
# Observe the output of timestamp column for missing values as NaT
df3.ix[['I','III','VIII'],['A','timestamp']] = np.nan
df3

s = Series([5,6,7,8,9])
s.loc[0] = None
s

s = Series(["A", "B", "C", "D", "E"])
s.loc[0] = None
s.loc[1] = np.nan
s

# Fillna method to fill the missing value
df2
df2.fillna(0)  # fill all missing value with 0
df2['D'].fillna('missing') # fill particular column missing value with missing

df2.fillna(method='pad')
df2
df2.fillna(method='pad', limit=1)

df2.dropna(axis=0)
df2.dropna(axis=1)

ts
ts.count()
ts[10:30]=None
ts.count()
# interpolate method perform interpolation to fill the missing values
# By default it performs linear interpolation 
ts.interpolate()
ts.interpolate().count()

输入输出操作

熊猫输入输出 API 是一组返回熊猫对象的读取器函数。使用熊猫捆绑的工具加载数据非常容易。数据从各种类型文件中的记录加载到 pandas 数据结构中，如 逗号分隔值 ( CSV )、Excel、HDF、SQL、JSON、HTML、谷歌大查询、pickle、stats 格式和剪贴板。有几个阅读器功能—每种文件一个功能—即read_csv、read_excel、read_hdf、read_sql、read_json、read_html、read_stata、read_clipboard和read_pickle。加载后，数据准备好进行分析。这包括删除错误条目、规范化、分组、转换和排序。

正在处理 CSV 文件

下一个程序演示如何处理 CSV 文件并对其执行各种操作。该程序使用 CSV 格式的图书交叉数据集，从www2.informatik.uni-freiburg.de/~cziegler/B…下载。它包含三个 CSV 文件(BX-Books.csv、BX-Users.csv、BX-Book-Ratings.csv)。其中包括书籍、用户和用户对书籍的评分的详细信息。有两个传递 CSV 文件名的选项；我们可以将文件放在任意文件夹中并使用完整路径，也可以将文件保留在当前目录中，只传递其名称。以下程序中的文件路径是 Windows 操作系统上的完整路径:

import numpy as np
randn = np.random.randn
from pandas import *

user_columns = ['User-ID', 'Location', 'Age']
users = read_csv('c:\BX-Users.csv', sep=';', names=user_columns)

rating_columns = ['User-ID', 'ISBN', 'Rating']
ratings = read_csv('c:\BX-Book-Ratings.csv', sep=';', names=rating_columns)

book_columns = ['ISBN', 'Book-Title', 'Book-Author', 'Year-Of-Publication', 'Publisher', 'Image-URL-S']
books = read_csv('c:\BX-Books.csv', sep=';', names=book_columns, usecols=range(6))

books
books.dtypes

users.describe()
print books.head(10)
print books.tail(8)
print books[5:10]

users['Location'].head()
print users[['Age', 'Location']].head()

desired_columns = ['User-ID', 'Age'] 
print users[desired_columns].head()

print users[users.Age > 25].head(4)
print users[(users.Age < 50) & (users.Location == 'chicago, illinois, usa')].head(4)

print users.set_index('User-ID').head()
print users.head()

with_new_index = users.set_index('User-ID')
print with_new_index.head()
users.set_index('User_ID', inplace=True)
print users.head()

print users.ix[62]
print users.ix[[1, 100, 200]]
users.reset_index(inplace=True)
print users.head()

这里有一个程序，演示了merge、groupby以及相关操作，如排序、排序、查找最上面的 n 值以及在跨书数据集上的聚合:

import numpy as np
randn = np.random.randn
from pandas import *

user_columns = ['User-ID', 'Location', 'Age']
users = read_csv('c:\BX-Users.csv', sep=';', names=user_columns)
rating_columns = ['User-ID', 'ISBN', 'Rating']
ratings = read_csv('c:\BX-Book-Ratings.csv', sep=';', names=rating_columns)

book_columns = ['ISBN', 'Title', 'Book-Author', 'Year-Of-Publication', 'Publisher', 'Image-URL-S']
books = read_csv('c:\BX-Books.csv', sep=';', names=book_columns, usecols=range(6))

# create one merged DataFrame
book_ratings = merge(books, ratings)
users_ratings = merge(book_ratings, users)

most_rated = users_ratings.groupby('Title').size().order(ascending=False)[:25]
print most_rated

users_ratings.Title.value_counts()[:17]

book_stats = users_ratings.groupby('Title').agg({'Rating': [np.size, np.mean]})
print book_stats.head()

# sort by rating average
print book_stats.sort([('Rating', 'mean')], ascending=False).head()

greater_than_100 = book_stats['Rating'].size >= 100
print book_stats[greater_than_100].sort([('Rating', 'mean')], ascending=False)[:15]

top_fifty = users_ratings.groupby('ISBN').size().order(ascending=False)[:50]

注

以下程序在https://github . com/gjreda/Greg reda . com/blob/master/content/notebooks/data/city-of-Chicago-sales . CSV 上的 CSV 文件上工作？原始=真。

这个程序演示了数据帧的合并和连接操作:

import numpy as np
randn = np.random.randn
from pandas import *

first_frame = DataFrame({'key': range(10), 
                           'left_value': ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J']})
second_frame = DataFrame({'key': range(2, 12), 
                           'right_value': ['L', 'M', 'N', 'O', 'P', 'Q', 'R', 'S', 'T', 'U']})
print first_frame
print second_frame

#Natural Join Operation 
print merge(left_frame, right_frame, on='key', how='inner') 
# Left, Right and Full Outer Join Operation
print merge(left_frame, right_frame, on='key', how='left')
print merge(left_frame, right_frame, on='key', how='right')
print merge(left_frame, right_frame, on='key', how='outer')

concat([left_frame, right_frame])
concat([left_frame, right_frame], axis=1)

headers = ['name', 'title', 'department', 'salary']
chicago_details = read_csv('c:\city-of-chicago-salaries.csv',
                      header=False,
                      names=headers,
                      converters={'salary': lambda x: float(x.replace('$', ''))})
print chicago_detail.head()

dept_group = chicago_details.groupby('department')

print dept_group
print dept_group.count().head(10) 
print dept_group.size().tail(10) 
print dept_group.sum()[10:17] 
print dept_group.mean()[10:17] 
print dept_group.median()[10:17] 

chicago_details.sort('salary', ascending=False, inplace=True)

即食数据集

有各种各样的数据来源关于经济学和在熊猫项目中使用这些数据的具体模块。我们可以使用pandas.io.data和pandas.io.ga(谷歌分析)模块从各种互联网来源提取数据，并将其添加到数据框中。目前，它支持以下来源:

• 雅虎！金融
• 谷歌金融
• 圣路易斯联邦储备银行:美联储经济数据 ( 弗雷德 ) 是一个由来自 80 个来源的 267，000 多个经济时间序列组成的数据库
• 肯尼斯·弗伦奇的数据库
• 世界银行
• 谷歌分析

这里有一个小程序，演示了如何从这些数据源中读取数据:

import pandas.io.data as web
import datetime
f1=web.DataReader("F", 'yahoo', datetime.datetime(2010, 1, 1), datetime.datetime(2011, 12, 31))
f2=web.DataReader("F", 'google', datetime.datetime(2010, 1, 1), datetime.datetime(2011, 12, 31))
f3=web.DataReader("GDP", "fred", datetime.datetime(2010, 1, 1), datetime.datetime(2011, 12, 31))
f1.ix['2010-05-12']

熊猫在密谋

熊猫数据结构支持包装在plt.plot()方法周围的绘图方法，用于在数据结构中绘制数据。默认情况下，它将显示线图，可以通过将名为“种类”的可选属性传递给绘图方法来更改该线图。以下列表包含了df.plot()中用于生成不同地块的变化:

• 条形图 : df.plot(kind='bar')
• 柱状图 : df.plot(kind='hist')
• 方块图 : df.plot(kind='box')
• 区域地块 : df.plot(kind='area')-
• 散射图 : df.plot(kind='scatter')
• 脚迹 : df.plot(kind='pie')

这个程序演示了熊猫包装方法的一个简单的绘图例子。程序的输出显示在它后面的屏幕截图中:

from pandas import *
randn = np.random.randn
import matplotlib.pyplot as plt
x1 = np.array( ((1,2,3), (1,4,6), (2,4,8)) )
df = DataFrame(x1, index=['I', 'II', 'III'], columns=['A', 'B', 'C']) 
df = df.cumsum()
df.plot(kind='pie', subplots=True)
plt.figure()
plt.show()

IPython

IPython 的设计和开发旨在提供一个增强的 Python 外壳，使得执行交互式分布式和并行计算成为可能。IPython 还有一套为科学计算构建专用交互环境的工具。它有两个组件来帮助实现 IPython 的目标:

• 增强的交互式 IPython 外壳
• 交互式并行计算的体系结构

在本节中，我们将首先讨论增强交互式 IPython 外壳的组件。我们将在第 8 章、并行和大规模科学计算中介绍交互式并行计算的其他组件。

IPython 控制台和系统外壳

IPython 提供的界面如下图所示。我们可以对这个控制台应用不同的着色方案；默认着色方案为NoColor。我们还有其他选择，如Linux和LightBG。IPython 的一个重要特性是它是有状态的，因为它维护在控制台上执行的计算的状态。IPython 中任何步骤的输出都存储在_N中，其中N是输出/结果的数量。当我们进入 IPython 交互式外壳时，它会显示这个增强的交互式 IPython 提供的功能，如下所示:

IPython 3.0.0 -- An enhanced Interactive Python.
?         -> Introduction and overview of IPython's features.
%quickref -> Quick reference.
help      -> Python's own help system.
object?   -> Details about 'object', use 'object??' for extra details.

如果我们在 shell 中输入问号(?)作为命令，那么它将显示 IPython 特性的详细列表。类似地，%quickref将显示大量 IPython 命令的简短引用，%magic将显示 IPython 魔法命令的详细信息。

如果我们输入任意一个objectname?，那么控制台将显示该对象的所有细节，例如文档字符串、函数和构造函数，如下图所示。我们已经创建了一个名为df的数据框对象，并使用df?显示了它的细节。

操作系统界面

有许多情况需要在操作系统的支持下执行计算。用户可以为常用命令创建新的别名。还支持ls等 Unix 命令。用户可以在任何操作命令或 shell 脚本前面加上!来执行它。

在 Python shel 中执行的操作系统命令

非阻塞绘图

在正常的 Python shell 中，如果我们创建任何一个图，并使用show()方法显示它，那么这个图将显示在一个新的屏幕上，这将保持 shell 被阻止，直到用户关闭显示图的屏幕。但是，IPython 有一个名为–pylab的标志。如果我们使用IPython –pylab命令执行 IPython 外壳，那么从 IPython 外壳打开的绘图窗口将不会阻塞外壳。这显示在下面的截图中——一个绘图窗口在不阻挡外壳的情况下打开，因为 IPython 是用–pylab标志执行的:

调试

IPython 拥有对程序调试以及错误和异常跟踪的出色支持。脚本执行完毕后，我们可以调用%debug启动 Python 调试器(pdb)来检查问题。我们可以在这里执行调试活动，因为我们可以打印变量值，执行语句，跟踪特定问题的来源。通常，这避免了使用外部调试器应用。

该截图描述了%debug选项:

用户可以通过调用%run -d programname.py逐步执行任何程序。这在下面的截图中显示。我们在名为stepbystep.py的程序中有一个步骤。在每个断点处，调试器界面要求用户按下 C 继续下一步:

IPython 笔记本电脑

IPython 有一个名为 Notebook 的网络应用。它是为交互式开发和文字计算创作而设计和开发的，其中解释概念、数学方面、实际计算和图形输出的文本可以组合在一起。程序的输入和程序的输出存储在单元格中，如果需要，这些单元格可以就地编辑。

以下截图取自ipython.org/notebook.ht…，展示了 IPython 的界面:

总结

在本章中，我们从讨论 matplotlib 的基本概念和体系结构开始。之后，我们讨论了一些用于生成不同类型地块的示例程序。我们还介绍了将这些图保存在不同格式文件中的方法。然后我们讨论了熊猫在数据分析中的应用。

此外，我们还讨论了大熊猫的数据结构。在深入介绍了数据结构的使用之后，您学习了如何执行各种其他相关的数据分析活动。在最后一部分，我们讨论了使用 IPython 的交互式计算的概念、用途和应用。

在下一章中，我们将全面讨论使用 Python 进行科学计算，其中涉及并行和高性能计算。本章将涵盖并行和高性能计算的基本概念，以及可用的框架和技术。稍后，它将深入介绍 Python 在并行和高性能计算中的应用。

八、并行和大规模科学计算

本章讨论了在 Python 中使用并行和大规模计算，或者使用 IPython 解决科学计算问题的重要概念。它涵盖了大规模科学计算和大数据处理的最新趋势。我们将使用示例程序来理解这些概念。

在本章中，我们将涵盖以下主题:

• IPython 中并行计算的基础
• IPython 并行计算的组成部分
• IPython 的任务接口和数据库
• IPython 的直接接口
• IPython 并行计算的细节
• IPython 中的 MPI 程序
• Python 中使用 Hadoop 和 Spark 的大数据处理

IPython 运行许多不同的进程，使用户能够执行并行计算。这些过程中的第一个是 IPython 引擎，它是一个 Python 解释器，执行用户提交的任务。用户可以运行多个引擎来执行并行计算。第二个过程是 IPython 中枢，它监控引擎和调度器，以跟踪用户任务的状态。集线器进程侦听来自引擎和客户端的注册请求；它持续监控来自调度程序的连接。第三个进程是 IPython 调度程序。这是一组用于在客户端和引擎之间传递命令和结果的进程。通常，调度程序进程运行在运行控制器进程的机器上，并连接到中心进程。最后一个进程是 IPython 客户端，它是 IPython 会话，用于协调引擎执行并行计算。

所有这些进程统称为 IPython 集群。这些进程使用 ZeroMQ 进行相互通信。ZeroMQ 支持各种传输协议，包括 Infiband、IPC、PGM、TCP 等。由集线器和调度器组成的 IPython 控制器监听套接字上客户端的请求。当用户启动引擎时，它会连接到集线器并执行注册。现在，集线器与引擎交换调度器的连接信息。稍后，引擎连接到调度器。这些连接在每台发动机的整个使用寿命中持续存在。每个 IPython 客户端都使用许多套接字连接来连接到控制器。通常，每个调度程序使用一个连接，一个集线器使用三个连接。这些连接在客户的整个生命周期中保持不变。

使用 IPython 的并行计算

IPython 允许用户以交互方式执行并行和高性能计算。我们可以使用 IPython 对并行计算的内置支持，它由四个组件组成，使 IPython 适合大多数类型的并行。具体来说，IPython 支持以下类型的并行:

• 单程序多数据并行 ( SPMD ):这是最常见的并行编程风格，是多指令多数据 ( MIMD )的一个子类型。在这个模型中，每个任务执行同一个程序的自己的副本。每个任务处理不同的数据集以获得更好的性能。
• 多程序、多数据并行:在多程序、多数据 ( MPMD )风格中，每个任务执行不同的程序，在每个参与计算节点上处理不同的数据集。
• 使用 MPI 的消息传递:一个消息传递接口 ( MPI )是消息传递库开发者的规范。它是一个独立于语言的规范，使用户能够编写基于消息传递的并行程序。在其目前的形式中，它支持分布式共享内存模型及其混合模型。
• 任务并行:任务并行在计算涉及的不同节点之间分配执行进程。该任务可以是线程、消息传递组件或其他编程模型的组件，如 MapReduce。
• 数据并行:数据并行在并行计算涉及的不同节点之间分配数据。数据并行化的主要重点是数据在不同节点间的分布/并行化，这与任务并行化不同。
• 上述类型的混合并行 : IPython 还支持任何上述风格的混合并行计算风格。
• 用户定义的并行方法 : IPython 被设计得非常简单和高度灵活，这种设计重点使用户能够将其用于任何新的或用户定义的并行风格。

IPython 以交互方式支持所有类型并行应用的程序开发生命周期的各个阶段，如开发、执行、调试、监控等。

将 matplotlib 与 IPython 结合使用，使用户能够分析和可视化远程或分布式大型数据集。它还使他们能够在集群上开始作业处理，并在本地系统上提取数据进行分析和可视化。用户可以从台式机/笔记本电脑上的 IPython 会话将 MPI 应用推送到高性能计算机上。它还支持在一组 CPU 上运行的不同任务的动态负载平衡。此外，它支持许多简单的方法，允许用户用两行或三行代码交互式地并行化许多简单的应用。用户可以交互地开发、执行、测试和调试定制的并行算法。IPython 使用户能够将在不同计算节点上执行的不同 MPI 作业捆绑到单个、巨大的分布式和/或并行系统中。

IPython 并行计算的体系结构

IPython 中并行计算的体系结构有三个主要组成部分。这些组件是 IPython 并行包的一部分。IPython 并行计算的体系结构如下图所示:

IPython 并行计算的三个主要组件是客户端、控制器和引擎。控制器组件由两个子组件组成: HUB 和 SCHEDULERS 。它允许客户端通过两个主要接口与引擎交互:直接接口和负载平衡接口。

并行计算的组成部分

本小节将讨论与 IPython 并行计算架构相关的各种组件和概念。这些组件是 IPython 引擎、IPython 控制器(集线器和调度器)以及 IPython 客户端和视图。

IPython 发动机

核心组件执行作为网络请求接收的 Python 命令的实际执行。该引擎是一个普通 Python 解释器的实例，最终它将成为一个完整的 IPython 解释器。用户可以通过启动多个引擎来执行分布式计算和并行计算。用户代码在阻塞模式下在 IPython 引擎中执行。

IPython 控制器

控制器由一个中枢和一组调度器组成。IPython 控制器是客户端和引擎用于通信的一组进程。它是让引擎执行 Python 进程的用户的联系点。通常，调度程序是在集线器运行的同一台机器上运行的独立进程。在某些特殊情况下，调度程序在远程计算机上运行:

• Hub:Hub 是最重要的组件，它跟踪调度器、客户端以及与引擎的连接。它处理客户端和引擎的所有连接以及整个流量。它还维护所有请求和结果的持久数据库，这些请求和结果将在应用的后续阶段使用。集线器提供了查询集群状态的功能，并隐藏了客户端和引擎之间许多连接的实际细节。
• 调度器:提交给引擎处理的 Python 命令通过调度器被定向。调度器还解决了执行用户代码时引擎阻塞的问题。调度器设法对用户隐藏这个问题，并提供对 IPython 引擎集合的完全异步访问。

IPython 视图和界面

控制器提供两个接口与引擎交互。第一个接口是Direct接口，其中引擎被直接寻址用于任务分配。第二个接口是LoadBalanced接口，其中任务到引擎的正确分配留给调度器。IPython 的灵活设计使我们能够扩展视图，以实现更复杂的接口方案。

对于控制器的不同连接方式，有一个view对象。以下是通过控制器与机器交互的两种模式:

• DirectView类，支持直接寻址。它允许在所有引擎上执行命令。
• LoadBalancedView类以负载平衡的方式代表用户处理任务农业。它允许在任何一个引擎上执行命令，用于执行的引擎由调度程序决定。

IPython 客户端

客户端是用于连接到群集的对象。在创建客户端对象的过程中，用户可以选择前面讨论的两个视图中的任何一个。一旦创建了客户端，只要作业运行，它就会一直存在。当超时周期结束或用户调用kill功能时，它被破坏。

执行并行计算的示例

以下程序是使用 IPython 执行并行计算的一个简单示例。它计算单个引擎或所有引擎中并行的集群的功率。在执行该程序之前，建议您根据需要检查zmq包是否安装。

要在 IPython 中运行这些程序，首先使用ipcluster start --n=4 --profile=testprofile命令启动 IPython 集群。它将在<userhome>/.ipython/profile_testprofile/security目录中创建ipcontroller-engine.json和ipcontroller-client.json文件。当我们通过传递profile='testprofile'创建客户端时，这些文件将被搜索。如果我们使用parallel.Client()创建客户端，那么它将在profile_default文件夹中搜索 JSON 文件。

首先，程序定义了一个计算功率的函数，然后使用测试配置文件创建一个客户端。要在引擎中调用 Python 函数，我们可以使用客户端或视图的apply方法。Python map函数对序列执行串行计算。在DirectView和LoadBalancedView中均有map函数对序列执行并行计算。我们也可以在阻塞或非阻塞模式下执行这些调用。要设置阻挡模式，我们可以将客户端或视图的block属性设置为true；默认为false:

from IPython import parallel
def pow(a, b):
  return a ** b
clients = parallel.Client(profile='testprofile')
print clients.ids
clients.block = True
clients[0].apply(pow, 2, 4)
clients[:].apply(pow, 2, 4)
map(pow, [2, 3, 4, 5], [2, 3, 4, 5])
view = clients.load_balanced_view()
view.map(pow, [2, 3, 4, 5], [2, 3, 4, 5])

一个并行的装饰器

在DirectView中有一个并行装饰器，创建parallel函数。该函数对序列进行操作，并分解元素式操作。随后，它将它们分配给并行计算，最后，它会重建结果。LoadBalancedView的装饰器把 Python 函数变成了parallel函数:

from IPython import parallel
clients = parallel.Client(profile='testprofile')
lbview = clients.load_balanced_view()
lbview.block = True
serial_computation = map(lambda i:i**5, range(26))
parallel_computation = lbview.map(lambda i: i**5, range(26))
@lbview.parallel()
def func_turned_as_parallel(x):
     return x**8
func_turned_as_parallel.map(range(26))

IPython 的神奇功能

IPython 有许多用户可以作为命令调用的神奇功能。IPython 中有两种类型的魔法命令，即线魔法和细胞魔法。线路魔法功能以%为前缀，像操作系统命令一样执行它们的功能。而细胞魔法函数的前缀是%%，它们把剩下的一行和后面的一行作为不同的参数。

当用户创建客户端时，这些神奇的功能变得可用。线魔函数的描述如下:

• %px:这可以在选定的引擎上执行单个 Python 命令。用户可以通过设置视图实例的目标属性来选择引擎。
• %pxconfig:即使没有任何活动视图，也可以使用pxconfig魔法功能指定--targets、--block、–noblock。
• %autopx:这是并联和非并联模式的切换开关。在第一次调用时，它会将控制台切换到一种模式，在这种模式下，所有键入的命令/函数调用将以并行模式执行，直到用户再次调用%autopx。
• %pxresult:在非阻塞模式下，%px不返回结果。我们可以看到使用pxresult魔法命令的最新命令的结果。

在 cell magic 模式下，px ( %%px ) magic 接受--targets选项指定要使用的目标引擎，--block或--noblock指定阻塞或非阻塞执行模式。这在我们没有视图实例的情况下特别有用。它还有一个参数--group-output，可以管理多个引擎输出的呈现。

下面的程序说明了px和pxresult作为线魔法和细胞魔法的使用。它还涵盖了autopx和pxconfig线魔法，并为这些魔法创建了特定的后缀。程序的第二行和第三行在 IPython 会话和所有引擎上执行导入。第二行之后创建的块内的所有导入也将在引擎上执行:

from IPython import parallel
drctview = clients[:]
with drctview.sync_imports():
   import numpy
clients = parallel.Client(profile='testprofile')
drctview.activate()
drctview.block=True
%px dummymatrix = numpy.random.rand(4,4)
%px eigenvalue = numpy.linalg.eigvals(dummymatrix)
drctview['eigenvalue']

%pxconfig --noblock
%autopx
maximum_egnvals = []
for idx in range(50):
    arr = numpy.random.rand(10,10)
    egnvals = numpy.linalg.eigvals(arr)
    maximum_egnvals.append(egnvals[0].real)
%autopx
%pxconfig --block 
%px answer= "The average maximum eigenvalue is: %f"%(sum(maximum_egnvals)/len(maximum_egnvals))
dv['answer']

%%px --block --group-outputs=engine
import numpy as np
arr = np.random.random (4,4)
egnvals = numpy.linalg.eigvals(arr)
print egnvals
egnvals.max()
egnvals.min()

odd_view = clients[1::2]
odd_view.activate("_odd")
%px print "Test Message"
odd_view.block = True
%px print "Test Message"
clients.activate()
%px print "Test Message"
%px_odd print "Test Message"

激活特定视图

默认情况下，这些魔法函数与一个DirectView对象相关联。允许用户通过在任何特定视图上调用activate()方法来更改DirectView对象。激活视图时，我们可以提到一个新的后缀，如odd_view.activate("_odd")中定义的。对于这个观点，我们现在有了一套新的魔法函数以及原来的魔法函数，比如%px_odd，它用在前面程序的最后一行。

引擎和 QtConsole

px神奇功能允许用户将 QtConsole 连接到引擎进行调试。下面的程序片段演示了如何通过绑定引擎内核来监听连接，从而将 QtConsole 连接到引擎:

%px from IPython.parallel import bind_kernel; bind_kernel()
%px %qtconsole
%px %connect_info

IPython 的高级特性

在接下来的小节中，我们将讨论 IPython 的各种高级特性。

容错执行

IPython 任务接口将引擎准备为容错和动态负载平衡的集群系统。在任务界面中，用户无权访问引擎。相反，任务分配完全依赖于调度器，这使得界面的设计简单、灵活且强大。

如果任务在 IPython 中由于任何原因失败，那么该任务将被重新排队，并再次尝试执行。如果出现故障，用户可以配置系统进行预定义次数的重试，他们也可以重新提交任务。

如果需要，用户可以显式重新提交任何任务。或者，他们可以设置一个标志，通过设置视图或计划程序的标志，在预定义的次数内重试该任务。

如果用户确定错误的原因不是代码中的错误或问题，那么他们可以将重试标志设置为从 1 到引擎总数的任何整数值。

最大限制等于引擎数量的原因是任务不会被重新提交给失败的引擎。

有两个选项可以设置重新提交次数的标志值。一种是使用LoadBalancedView(考虑对象名称为lbvw)对象设置值后，设置所有后续任务，如下:

lbvw.retries = 4

另一种是使用with ...temp_flags为单个块设置值，如下所示:

with lbvw.temp_flags(retries=4):
    lbview.apply(task_tobe_retried)

动态负载平衡

调度器还可以被配置为基于各种调度策略来执行调度。IPython 支持多种方案，在负载平衡请求的情况下将任务分配给机器。集成定制方案也非常容易。有两种选择方案的方法。一种是设置控制器的config对象的taskSchedulerscheme_name属性。第二种选择是通过将方案参数传递给ipcontroller来选择方案，如下所示:

ipcontroller --scheme=<schemename>

这里有一个例子:

ipcontroller --scheme=lru

<schemename>功能可以是以下任何一种:

• lru : 最近最少使用的 ( LRU )是将任务分配给最近最少使用的引擎的方案。
• plainrandom:在这个方案中，调度器随机挑选一个引擎运行任务。
• twobin:该方案使用 NumPy 函数分配任务。它是plainrandom和lru的组合，因为它随机选择两个引擎，并从这两个引擎中选择最近最少使用的。
• leastload:这个方案是调度器的默认方案。它将任务分配给负载最小的引擎(即剩余任务数最少的引擎)。
• weighted:这个方案是twobin方案的变体，因为它随机挑选两个引擎，并将未完成任务的负载或数量分配为权重的倒数。它将任务分配给负载相对较小的发动机。

在客户端和引擎之间推拉对象

除了在引擎上调用函数和执行代码，IPython 还允许用户在 IPython 客户端和引擎之间移动 Python 对象。push方法将对象从客户端推送到引擎，pull方法可用于将任意对象从引擎拉回客户端。在非阻塞模式下，推拉返回AsyncResult对象。要在非阻塞模式下显示结果，我们可以如下拉取对象:rslt = drctview.pull(('a','b','c'))。我们可以调用rslt.get()来显示被拉对象中的值。在某些情况下，对输入数据序列进行分区并将不同的分区推送到不同的引擎是非常有用的。这种划分实现为scatter和gather功能，类似于 MPI。scatter操作用于将分区序列从客户端(IPython 会话)推送到引擎，而gather操作用于将分区从引擎取回客户端。

所有这些功能将在下面的程序中演示。最后，使用scatter和gather实现两个矩阵的平行点积:

import numpy as np
from IPython import parallel
clients = parallel.Client(profile='testprofile')
drctview = clients[:]
drctview.block = True
drctview.push(dict(a=1.03234,b=3453))
drctview.pull('a')
drctview.pull('b', targets=0)
drctview.pull(('a','b'))
drctview.push(dict(c='speed'))
drctview.pull(('a','b','c'))
drctview.block = False
rslt = drctview.pull(('a','b','c'))
rslt.get()

drctview.scatter('a',range(16))
drctview['a']
drctview.gather('a')

def paralleldot(vw, mat1, mat2):
    vw['mat2'] = mat2
    vw.scatter('mat1', mat1)
    vw.execute('mat3=mat1.dot(mat2)')
    return vw.gather('mat3', block=True)
a = np.matrix('1 2 3; 4 5 6; 7 8 9')
b = np.matrix('4 5 6; 7 8 9; 10 11 12')
paralleldot(drctview, a,b)

下面的程序演示了将对象从客户端推送到引擎以及将结果从引擎拉回到客户端的方法。它在所有引擎上执行两个矩阵的点积，最后收集结果。它还使用allclose()方法验证所有结果是否相同，如果对象相同，则返回True。在以下程序的execute命令中，添加了print mat3语句，目的是使用display_outputs()方法显示所有发动机的标准输出设备的输出:

import numpy as np
from IPython.parallel import Client
ndim = 5
mat1 = np.random.randn(ndim, ndim)
mat2 = np.random.randn(ndim, ndim)
mat3 = np.dot(mat1,mat2)
clnt = Client(profile='testprofile')
clnt.ids
dvw = clnt[:]
dvw.execute('import numpy as np', block=True)
dvw.push(dict(a=mat1, b=mat2), block=True)
rslt = dvw.execute('mat3 = np.dot(a,b); print mat3', block=True)
rslt.display_outputs()
dot_product = dvw.pull('mat3', block=True)
print dot_product
np.allclose(mat3, dot_product[0])
np.allclose(dot_product[0], dot_product[1])
np.allclose(dot_product[1], dot_product[2])
np.allclose(dot_product[2], dot_product[3])

用于存储请求和结果的数据库支持

IPython 中枢存储关于请求的信息和处理任务的结果，供以后使用。它的默认数据库是 SQLite，目前支持 MongoDB 和一个名为DictDB的内存数据库。用户必须配置用于他们的配置文件的数据库。在活动的配置文件文件夹中，有一个名为ipcontroller_config.py的文件。这个文件将在我们启动ipcluster时创建。该文件有一个c.HubFactory.db_class条目；用户应该将其设置到他们选择的数据库中，如下所示:

#dict-based in-memory database named as dictdb
c.HubFactory.db_class = 'IPython.parallel.controller.dictdb.DictDB'
# For MongoDB:
c.HubFactory.db_class = 'IPython.parallel.controller.mongodb.MongoDB'
# For SQLite:
c.HubFactory.db_class = 'IPython.parallel.controller.sqlitedb.SQLiteDB'

该属性的默认值为NoDB，表示不使用数据库。要获得任何已执行任务的结果，用户可以在客户端对象上调用get_result函数。客户对象有一个更好的方法叫做db_query()来获得更多关于任务结果的见解。该方法是以 MongoDB 查询方式设计的。它从具有精确值的TaskRecord关键字列表中获取一个带有关键字的字典查询对象，或者 MongoDB 查询。这些论点遵循{'operator' : 'argument(s)'}语法。它还有一个名为keys的可选参数。此参数用于指定要检索的键。它返回一个TaskRecord字典列表。默认情况下，它检索除请求和结果缓冲区之外的所有键。msg_id键将始终包含在响应中，类似于 MongoDB。下表解释了各种任务记录键:

• msg_id:该值为uuid(字节)类型。它代表消息标识。
• header:该值为dict类型，保存请求头。
• content:该值为dict类型，保存一般为空的请求内容。
• buffers:这个值是list(字节)类型，它将是一个包含序列化请求对象的缓冲区。
• Submitted:该值为datetime类型，保存提交时间戳。
• client_uuid:该值是 uuid(以字节为单位的通用唯一标识符)。
• engine_uuid:该值为保存发动机插座标识的uuid(字节)类型。
• started:该值为datetime类型，保存某个引擎上任务执行开始的时间。
• completed:该值为datetime类型，保存任务在引擎上执行完成的时间。
• resubmitted:该值为datetime类型，如果适用，保存任务重新提交的时间。
• result_header:该值为dict类型，保存结果的表头。
• result_content:该值为dict类型，保存结果的内容。
• result_buffers:这个值是list(字节)类型，它将是一个包含序列化结果对象的缓冲区。
• queue:该值为bytes类型，代表任务的队列名称。
• stdout:这是一串标准输出 ( 标准输出)数据。
• stderr:这是一串标准误差 ( 标准误差)数据。

以下程序演示了用于访问结果信息的db_query()和get_result()方法的概念:

from IPython import parallel
from datetime import datetime, timedelta
clients = parallel.Client(profile='testprofile')
incomplete_task = clients.db_query({'complete' : None}, keys=['msg_id', 'started'])
one_hourago = datetime.now() - timedelta(1./24)
tasks_started_hourago = clients.db_query({'started' : {'$gte' : one_hourago },'client_uuid' : clients.session.session})
tasks_started_hourago_other_client = clients.db_query({'started' : {'$le' : hourago }, 'client_uuid' : {'$ne' : clients.session.session}})
uuids_of_3_n_4 = map(clients._engines.get, (3,4))
headers_of_3_n_4 = clients.db_query({'engine_uuid' : {'$in' : uuids_of_3_n_4 }}, keys='result_header')

在db_query中支持以下关系运算符作为 MongoDB:

• '$in':这代表列表/序列中操作的 ** '$nin':这表示列表/序列上的操作中没有** '$eq':这个用来表示等于 ( ==)* '$ne':这是用来表示不等于 ( !=)* '$gt':这个用来表示大于 ( >)* '$gte':表示大于等于 ( >=)* '$lt':这是用来表示小于 ( <)* '$lte':用于表示小于等于 到 ( <=)**

**## 在 IPython 中使用 MPI

通常，在多个引擎上运行的并行算法需要在引擎之间移动数据。我们已经介绍了 IPython 执行这种数据移动的内置方式。但是，这是一个缓慢的操作，因为它不是客户端和引擎之间的直接传输。数据必须通过控制器传输。获得良好性能的更好方法是使用消息传递接口 ( MPI )。IPython 的并行计算对与 MPI 的集成有出色的支持。要将 MPI 与 IPython 并行计算结合使用，我们需要安装一个 MPI 实现，如OpenMPI或MPICH2 / MPICH和mpi4py python 包。安装后，测试系统是否能够执行mpiexec或mpirun命令。

测试安装之后，在实际运行 MPI 程序之前，用户需要使用以下命令创建 MPI 执行的配置文件:

ipython profile create --parallel --profile=mpi

轮廓创建后，在profile_mpi文件夹的ipcluster_config.py中添加以下线条:

c.IPClusterEngines.engine_launcher_class = 'MPIEngineSetLauncher'

现在，系统准备在 IPython 上执行基于 MPI 的程序。用户可以使用以下命令启动群集:

ipcluster start -n 4 --profile=mpi

前面的命令启动 IPython 控制器，并使用mpiexec启动四个引擎。

下面的程序定义了一个计算分布式数组总和的函数。将文件名保存为parallelsum.py，因为这个名称将在下一个程序中使用，该程序实际上调用了这个函数:

from mpi4py import MPI
import numpy as np

def parallelsum(arr):
    localsum = np.sum(arr)
    receiveBuffer = np.array(0.0,'d')
    MPI.COMM_WORLD.Allreduce([localsum, MPI.DOUBLE],
        [receiveBuffer, MPI.DOUBLE],
        op=MPI.SUM)
    return receiveBuffer

现在调用前面程序中定义的函数，以便在多个引擎上执行它。这样做是为了对数组执行并行求和:

from IPython.parallel import Client
clients = Client(profile='mpi')

drctview = clients[:]
drctview.activate() 
#execute the program name passed as argument
drctview.run(parallelsum.py.py')
drctview.scatter('arr',np.arange(20,dtype='float'))
drctview['arr']
# calling of the function
%px sum_of_array = parallelsum(arr)
drctview['sum_of_array']

管理任务之间的依赖关系

它对管理各种任务之间的依赖关系有很强的支持。在大多数科学和商业应用中，只有负载平衡方案不足以管理其复杂性。这些应用需要多个任务之间的依赖性。这些依赖关系描述了特定的软件、Python 模块、操作系统或硬件需求；序列；计时；和任务集合中任务的执行位置。IPython 支持两种依赖，即函数依赖和图依赖。

功能依赖

功能相关性用于确定特定引擎是否能够运行任务。这个概念是使用来自IPython.parallel.error的特殊异常UnmetDependency实现的。如果任务因UnmetDependency异常而失败，调度程序不会将此错误传播给客户端。相反，它会处理这个错误，并将这个任务提交给其他引擎。调度程序重复这个过程，直到找到合适的引擎。此外，调度程序不会两次向引擎提交任务。

功能依赖的装饰者

虽然允许用户手动引发UnmetDependency异常，但是 IPython 提供了两个装饰器来管理这个依赖问题。有两个装饰器和一个用于函数依赖的类:

• @require: This decorator manages the dependency of a task that requires that a particular Python module, local function, or local object be available on an engine when the decorated function is called. Functions will be pushed to the engine with their names, and objects can be passed using the arg keyword. We can pass the names of all the Python modules required to execute this task. Using this decorator, a user can define a function to be executed on only those engines where the module names passed to this decorator are available and importable.

  例如，下面代码片段中定义的函数依赖于 NumPy 和 pandas 模块，因为它使用 NumPy 中的randn和 pandas 中的Series。如果对于某个任务，我们调用这个函数，那么它将在这两个模块可导入的机器上执行。调用此函数时，将导入 NumPy 和 pandas。

  from IPython.parallel import depend, require
  # the following function uses randn and Series
  @require('pandas', 'numpy')
  def func_uses_functions_from_numpy_pandas():
    return performactivity()
  

• @depend: This decorator allows users to define a function that has a dependency on some other function. It determines whether the dependency is met or not. Before starting the task, the dependency function will be called, and if this function returns true, then the actual processing of the task will be started. Moreover, if this dependency function returns false, then the dependency is considered to be unmet and the task is propagated to some other engine.

  例如，下面的代码片段首先创建一个依赖函数，用于验证引擎的操作系统是否与给定的操作系统匹配。这是因为用户希望编写两个不同的函数来执行 Linux 和 Windows 操作系统上的特定活动:

  from IPython.parallel import depend, require
  def find_operating_system(plat):
      import sys
      return sys.platform.startswith(plat)
  @depend(find_operating_system, 'linux')
  def linux_specific_task():
      perform_activity_on_linux()
  @depend(platform_specific, 'win')
  def linux_specific_windows():
      perform_activity_on_windows()
  

图形相关性

还有另一类重要的依赖关系，任务之间相互依赖，任务必须在某些或所有特定任务成功执行后才能执行。另一个依赖关系可能如下:任务必须在满足一组特定依赖关系的目标上执行。通常，用户需要一个选项来指定运行给定任务的时间和位置，作为时间、位置和其他任务结果的函数。有一个单独的名为Dependency的类来管理图形依赖关系，Dependency是类Set的子类。它包含一组与任务相对应的消息标识，并且它还具有一些属性。这些属性有助于检查是否满足指定的依赖关系:

• any|all:这些属性指定是否完成或满足任何指定的依赖关系。这将通过设置默认为True的所有依赖属性来指定。
• success:该属性默认为True，用于指定如果指定任务成功，则认为满足依赖关系。
• failure:该属性默认为False，用于指定如果指定任务失败，则认为满足依赖关系。
• after:该属性用于指定依赖任务在执行指定任务后执行。
• follow:follow属性指定从属任务应与其中一个从属任务在同一目的地执行。
• timeout:该属性用于指定调度程序必须等待依赖项满足的持续时间。默认为0，表示从属任务将永远等待。超时后，从属任务失败，出现DependencyTimeout异常。

有些任务可以作为清理任务。它们应该只在指定任务失败时运行。用户应使用failure=True,success=False进行此类任务。对于某些依赖任务，需要依赖任务成功完成。在这种情况下，用户必须设置success=True和failure=False。在某些情况下，用户希望独立于依赖任务的成功或失败来执行依赖任务。在这种情况下，用户必须使用success=failure=True。

不可能的依赖

可能有些依赖是不可能满足的。如果调度程序没有处理这种可能性，那么调度程序可能会一直等待依赖性得到满足。为了应对这种情况，调度程序会分析图的依赖关系，以估计满足依赖关系的可能性。如果调度程序能够识别出某个任务的依赖性无法满足，那么该任务将失败并出现ImpossibleDependency错误。下面的代码片段演示了任务之间的图形依赖关系的使用:

from IPython.parallel import *
clients = ipp.Client(profile='testprofile')
lbview = clients.load_balanced_view()

task_fail = lbview.apply_async(lambda : 1/0)
task_success = lbview.apply_async(lambda : 'success')
clients.wait()
print("Fail task executed on %i" % task_fail.engine_id)
print("Success task executed on %i" % task_success.engine_id)

with lbview.temp_flags(after=task_success):
    print(lbview.apply_sync(lambda : 'Perfect'))

with lbview.temp_flags(follow=pl.Dependency([task_fail, task_success], failure=True)):
    lbview.apply_sync(lambda : "impossible")

with lbview.temp_flags(after=Dependency([task_fail, task_success], failure=True, success=False)):
    lbview.apply_sync(lambda : "impossible")

def execute_print_engine(**flags):
    for idx in range(4):
        with lbview.temp_flags(**flags):
            task = lbview.apply_async(lambda : 'Perfect')
            task.get()
            print("Task Executed on %i" % task.engine_id)

execute_print_engine(follow=Dependency([task_fail, task_success], all=False))
execute_print_engine(after=Dependency([task_fail, task_success], all=False))
execute_print_engine(follow=Dependency([task_fail, task_success], all=False, failure=True, success=False))
execute_print_engine(follow=Dependency([task_fail, task_success], all=False, failure=True))

DAG 依赖关系和网络库

一般来说，最好用有向无环图 ( DAG )来表示并行工作流。Python 有一个流行的库，叫做 NetworkX，用于使用图形。该图是节点和有向边的集合。边连接各种节点；每个边都有一个相关的方向。我们可以用这个概念来表示依赖关系。例如，从任务 1 到任务 2 的edge(task1, task2)表示任务 2 依赖于任务 1。类似地，edge(task2, task1)表示任务 1 依赖于任务 2。这个图不能包含任何循环，这就是为什么它被称为非循环图。

现在考虑下图中描述的六节点 DAG。表示 Task0 不依赖任何任务，因此可以立即启动。而任务 1 和任务 2 依赖于任务 0 ，因此它们将在任务 0 完成后开始。那么任务 3 同时依赖于任务 1 和任务 2 ，所以在任务 1 和任务 2 结束后执行。同样的，任务 4 和任务 5 将在任务 3 结束后执行。任务 6 仅依赖于任务 4 。因此，它将在任务 4 完成执行后启动。

这是一个源代码片段，代表上图中描述的 DAG。在代码中，任务由它们的编号表示；即任务 0 用 0 表示，任务 1 用 1 表示，以此类推:

import networkx as ntwrkx
import matplotlib.pyplot as plt

demoDAG = ntwrkx.DiGraph()
map(demoDAG.add_node, range(6))
demoDAG.add_edge(0,1)
demoDAG.add_edge(0,2)
demoDAG.add_edge(1,3)
demoDAG.add_edge(2,3)
demoDAG.add_edge(3,4)
demoDAG.add_edge(3,5)
demoDAG.add_edge(4,6)
pos = { 0 : (0,0), 1 : (-1,1), 2 : (1,1), 3 : (0,2), 4 : (-1,3), 5 : (1, 3), 6 : (-1, 4)}
labels={}
labels[0]=r'$0$'
labels[1]=r'$1$'
labels[2]=r'$2$'
labels[3]=r'$3$'
labels[4]=r'$4$'
labels[5]=r'$5$'
labels[6]=r'$6$'

ntwrkx.draw(demoDAG, pos, edge_color='r')
ntwrkx.draw_networkx_labels(demoDAG, pos, labels, font_size=16)
plt.show()

以下程序用彩色边和顶点标签创建相同的图表:

import networkx as ntwrkx
import matplotlib.pyplot as plt

demoDAG = ntwrkx.DiGraph()
map(demoDAG.add_node, range(6))

pos = { 0 : (0,0), 1 : (-1,1), 2 : (1,1), 3 : (0,2), 4 : (-1,3), 5 : (1, 3), 6 : (-1, 4)}

ntwrkx.draw(demoDAG, pos)
ntwrkx.draw_networkx_edges(demoDAG,pos,
                       edgelist=[(0,1),(0,2),(1,3),(2, 3),(3, 4)],edge_color='r')
ntwrkx.draw_networkx_edges(demoDAG,pos,
                       edgelist=[(3,5),(4,6)],edge_color='b')

ntwrkx.draw_networkx_nodes(demoDAG,pos,
                       nodelist=[0,1,2,3,4],
                       node_color='r',
                       node_size=500,
                   alpha=0.8)
ntwrkx.draw_networkx_nodes(G,pos,
                       nodelist=[5,6],
                       node_color='b',
                       node_size=500,
                   alpha=0.8)

labels={}
labels[0]=r'$0$'
labels[1]=r'$1$'
labels[2]=r'$2$'
labels[3]=r'$3$'
labels[4]=r'$4$'
labels[5]=r'$5$'
labels[6]=r'$6$'

ntwrkx.draw_networkx_labels(demoDAG, pos, labels, font_size=16)
plt.show() 

在带有 StarCluster 的 Amazon EC2 集群上使用 IPython

StarCluster 旨在简化在亚马逊弹性计算云 ( EC2 )上使用虚拟机集群的过程。它是亚马逊 EC2 上集群计算的开源工具包。除了执行自动集群配置，星簇还提供亚马逊机器映像 ( AMIs )定制支持科学计算和软件开发的工具包和库。这些 AMIs 由 ATLAS、IPython、NumPy、OpenMPI、SciPy 等组成。用户可以在安装了 StarCluster 的机器上使用以下命令检索可用的 AMIs 列表:

starcluster listpublic

StarCluster 有一个非常简单直观的界面，用于计算集群的弹性管理和存储管理。安装后，用户必须更新其默认配置文件，以更新 Amazon EC2 帐户的详细信息，包括地址、地区、凭据和公钥/私钥对。

安装和配置完成后，用户可以使用以下命令从本地 IPython 安装中控制 Amazon EC2 集群:

starcluster shell --ipcluster=clusterName

如果配置中有任何错误，将通过前面的命令显示。如果配置是正确的，那么这个命令启动 StarCluster 的开发外壳，并在 Amazon EC2 上为远程集群配置一个并行会话。StarCluster 会自动创建一个名为ipclient的并行客户端，并以ipview的名称查看整个集群。用户可以使用这些变量(ipclient和ipview)在亚马逊 EC2 集群上运行并行任务。以下代码片段使用ipclient显示集群的引擎标识，并使用ipview运行一个小的并行任务:

ipclient.ids
result = ipview.map_async(lambda i: i**5, range(26))
print result.get()

用户还可以使用带有星簇的 IPython 并行脚本。如果用户希望从本地 IPython 会话在远程 Amazon EC2 集群上运行 IPython 并行脚本，那么他们应该在创建并行客户端的过程中使用一些配置细节，如下所示:

from IPython.parallel import Client
remoteclients = Client('<userhome>/.starcluster/ipcluster/<clustername>-<yourregion>.json',  sshkey='/path/to/cluster/keypair.rsa')

具体来说，假设一个集群的名称是packtcluster，区域是us-west-2，而keypair的名称是packtKey，存储在/home/user/.ssh/packtKey.rsa中。然后，前面的代码将更改为以下内容:

from IPython.parallel import Client
remoteclients = Client('/home/user/.starcluster/ipcluster/packtcluster-us-west-2.json', sshkey='/home/user/.ssh/packtKey.rsa')

在这两行代码之后，所有剩余的代码都将在 Amazon EC2 上的远程集群上执行。

关于 IPython 安全性的一个注记

在设计 IPython 的体系结构时，已经考虑到了安全问题。基于能力的客户端身份验证模型，以及 SSH 隧道传输的 TCP/IP 通道，管理主要的潜在安全问题，并允许用户在开放网络中使用 IPython 集群。

ZeroMQ 不提供安全性。因此，SSH 隧道是建立安全连接的主要来源。Client对象从ipcontroller-client.json文件获取与控制器建立连接的信息，然后使用 OpenSSH/Paramiko 创建隧道。

它还使用 HMAC 摘要的概念，使用保护共享机器用户的共享密钥对消息进行签名。有一个处理消息协议的会话对象。此对象使用唯一密钥验证消息的有效性。默认情况下，这个密钥是一个 128 位的伪随机数，类似于uuid.uuid4()生成的数字。一般来说，在并行计算过程中，IPython 客户端用于向 IPython 引擎发送 Python 函数、命令和数据，以执行和处理数据。IPython 确保只有客户端负责并能够使用引擎的功能。引擎从启动引擎的用户那里继承能力和权限。

为了防止未经授权的访问，身份验证和密钥相关的信息被编码在 JSON 文件中，客户端使用该文件来访问 IPython 控制器。用户可以通过限制密钥的分发来授予授权人员访问权限。

众所周知的并行编程风格

由于计算机硬件和软件的发展，并行程序可以用几种方式来设计、开发和实现。我们可以使用并发、并行或分布式的方式来实现一个程序。通常，前面提到的技术之一用于实现程序，以实现高效执行和提高性能。接下来的小节将讨论这些模型以及与之相关的常见问题。

并行编程中的问题

所有这些模型都依赖于程序的不同部分在独立的计算元素(中央处理器和计算节点)上执行的基本概念。通常，这些模型将程序分成多个工作程序，每个工作程序在不同的计算元素上开始执行。尽管有性能上的好处，这种类型的程序执行——使用多个工作人员——在通信中带来了多种复杂性。这个问题叫做进程间通信 ( IPC )。

有几个经典的 IPC 问题需要开发人员注意，即死锁、饥饿和竞争条件:

• 死锁:死锁是两个或两个以上的工作者处于无限等待状态，以获取另一个等待工作者所占用的资源的情况。对此有四个充要条件，即互斥、保持等待、不抢占和循环等待。如果在任何程序的执行过程中出现这些情况，则该程序将被阻止，并且无法继续执行:
  • 互斥意味着资源不可共享
  • 我们所说的“等待”是指处于死锁状态的每个工作人员都持有一些资源并请求一些额外的资源
  • 无优先权意味着分配给一个工作者的资源不能被优先分配给另一个工作者
  • 最后，通过循环等待，我们意味着死锁下的工人形成一个链，或循环列表，其中每个工人都在等待列表中下一个工人持有的资源
• 饥饿:当多个工人争夺单一资源时，就会出现这种情况。在这种情况下，每个工作人员都被分配了资源分配的优先级。有时，这种优先分配变得不公平，一些工人的执行被拖延了很长时间。假设有两种类型的员工在竞争一种资源:高优先级的员工和低优先级的员工。如果高优先级工作人员不断到来，那么可能会出现这样一种情况，即一些低优先级工作人员为了获得高优先级工作人员释放的资源而经历漫长(无限)的等待。
• 竞态条件:当有多个工作人员同时对公共数据执行读和写操作，并且他们执行的操作之间缺乏同步时，就会出现这个问题。例如，假设两个工作人员从数据库中读取一条公共数据，修改其值，然后将该数据写回数据库。如果这些操作没有以良好同步的顺序执行，它们将使数据库处于不一致的状态。

有一些技术可以用来避免这些问题，尽管它们的详细讨论超出了本书的范围。让我们在后续小节中讨论并行计算的类型。

并行编程

在并行风格的程序开发中，程序被分成多个工作人员，他们在独立的 CPU 上执行，而不会争夺一个 CPU，如下图所示。这些 CPU 可以是多核计算机的单个处理器，也可以位于不同的计算机上，并使用诸如消息传递接口 ( MPI )等技术进行通信。

并发编程

在并发编程中，用户程序的多个工作者在单个 CPU 上执行，或者在比工作者数量更少的 CPU 上执行(如下图所示)。这些工人在中央处理器调度程序的控制下争夺中央处理器。中央处理器调度器使用多种方案将工作人员分配给中央处理器。中央处理器调度方案被用来创建一个工人的排名，工人按照这个排名的顺序得到执行。

这些工作者可以使用多个进程或多个线程来实现。进程和线程同时执行主程序的某些部分。线程和进程的主要区别在于线程比进程消耗更少的内存。因此，螺纹被称为轻质。

分布式编程

在分布式编程中，程序的工作人员在通过网络连接的不同计算机上执行。这些程序的执行有不同的框架。网络还使用不同的拓扑，在某些情况下，方案数据和流程都可能是分布式的。随着时间的推移，这种并行计算模式变得越来越流行，因为它提供了几个优势，包括低成本、容错、高可扩展性等。在分布式编程中，每个组件都有独立的内存和处理，而在并行编程中，多个处理器/CPU 共享公共内存。

分布式程序执行

Python 中的多处理

多处理模块支持在多核环境下在多个处理器上独立运行的多个进程的创建和执行。它有两个支持多处理的重要模型；一个基于Process类，另一个基于Pool类。

该程序使用Process类演示多处理:

import multiprocessing as mpcs
import random
import string

output_queue = mpcs.Queue()

def strings_random(len, output_queue):
  generated_string = ''.join(random.choice(string.ascii_lowercase  + string.ascii_uppercase + string.digits)
    for i in range(len))
  output_queue.put(generated_string)

procs = [mpcs.Process(target=strings_random, args=(8, output_queue)) for i in range(7)]

for proc in procs:
  proc.start()

for proc in procs:
  proc.join()

results = [output_queue.get() for pro in procs]
print(results)

基于流程的类按照流程完成的顺序返回结果。如果用户需要检索一个有序的结果，他们必须付出额外的努力，如下面的代码所示。为了获得有序的结果，另一个参数被添加到函数中，最后添加到输出中。这表示过程的位置或顺序，最后结果根据参数排序。这个想法在下面的程序中使用Process类和输出中的一个位置参数进行了演示:

import multiprocessing as mpcs
import random
import string

output_queue = mpcs.Queue()

def strings_random(len, position, output_queue):
  generated_string = ''.join(random.choice(string.ascii_lowercase  + string.ascii_uppercase + string.digits)
    for i in range(len))
  output_queue.put((position, generated_string))

procs = [mpcs.Process(target=strings_random, args=(5, pos, output)) for pos in range(4)]

for proc in procs:
  proc.start()
for proc in procs:
  proc.join()

results = [output_queue.get() for pro in procs]
results.sort()
results = [rslt[1] for rslt in results]
print(results)

Pool类提供了用于并行计算的map和apply方法，并且还支持这些方法的异步版本。map和apply方法锁定主程序，直到一个过程完成，这个概念可以用来产生特定应用所需的顺序输出。

Python 中的多线程

Python 的线程模块允许用户创建一个进程的多个线程来执行并发计算。进程的线程与主进程/线程共享相同的数据空间，这使得数据共享和彼此之间的通信变得容易。线程也被称为轻量级进程，因为它们比进程需要更少的内存。

以下程序演示了线程的创建和启动:

import threading
import time
class demoThread (threading.Thread):
  def __init__(self, threadID, name, ctr):
    threading.Thread.__init__(self)
    self.threadID = threadID
    self.name = name
    self.ctr = ctr
  def run(self):
    print "Start of The Thread: " + self.name
    print_time(self.name, self.ctr, 8)
    print "Thread about to Exit:" + self.name

def print_time(threadName, delay, counter):
    while counter:
      time.sleep(delay)
      print "%s: %s" % (threadName, time.ctime(time.time()))
      counter -= 1

thrd1 = demoThread(1, "FirstThread", 4)
thrd2 = demoThread(2, "SecondThread", 5)
thrd1.start()
thrd2.start()
print "Main Thread Exits"

这个程序启动两个线程。如果你观察程序的输出，你会注意到没有线程执行的顺序。首先显示"Main Thread Exits"字符串，然后是线程名称和"Thread about to Exit: ThreadName"的随机序列。

我们可以选择同步这个输出，这样就可以保持线程完成的顺序。下面的程序首先执行第一个线程。然后，在其退出后，第二个线程被执行。最后，主线程退出。通过获取锁来维护线程序列，并且在退出之前释放该锁，以便可以启动第二个线程。主线程对所有线程对象调用join方法。该方法阻止主线程完成其他各种线程:

import threading
import time
class demoThread (threading.Thread):
  def __init__(self, threadID, name, ctr):
    threading.Thread.__init__(self)
    self.threadID = threadID
    self.name = name
    self.ctr = ctr
  def run(self):
    print "Start of The Thread: " + self.name
    threadLock.acquire()
    print_time(self.name, self.ctr, 8)
    print "Thread about to Exit:" + self.name
    threadLock.release()

def print_time(threadName, delay, counter):
    while counter:
      time.sleep(delay)
      print "%s: %s" % (threadName, time.ctime(time.time()))
      counter -= 1

threadLock = threading.Lock()
thrds = []

thrd1 = demoThread(1, "FirstThread", 4)
thrd2 = demoThread(2, "SecondThread", 5)
thrd1.start()
thrd2.start()

thrds.append(thrd1)
thrds.append(thrd2)
for thrd in threads:
  thrd.join()

print "Main Thread Exits"

Python 中基于 Hadoop 的 MapReduce

Hadoop 是一个开源框架，用于计算集群中巨大数据集的分布式存储和处理。Hadoop 系统有三个主要组件:用于处理的 MapReduce、 Hadoop 分布式文件系统 ( HDFS )和一个名为 HBase 的大规模数据库，用于存储数据集。HDFS 支持存储非常大的数据集文件。它将用户提交的数据集分布在各个集群节点上。它将数据集分成多个块，并保存关于块和节点的簿记信息。HBase 是一个为支持大规模数据库而设计的数据库，是在 HDFS 之上开发的。它是一个开源的、面向列的、非关系的分布式数据库。

MapReduce 是一个框架，旨在对计算集群中的巨大数据集执行分布式处理。Hadoop 是 MapReduce 框架的开源实现。MapReduce 程序由两个主要组件组成:map 和 Reduce。map函数用于对输入数据集进行过滤，并将其排序后的输出写入文件系统。稍后，reduce 函数使用该输出来执行总结，最终输出再次写入文件系统。MapReduce 框架遵循 单程序，多数据 ( SPMD )模型，因为它在多个数据集上处理同一个程序。

在 Hadoop 系统中，完整的功能分为许多组件。有两个主节点。一个是作业跟踪器，它跟踪地图并减少从节点的进程，称为任务跟踪器节点。第二个主节点是 namenode，它包含有关哪个分割数据文件块存储在称为 data node 的特定从节点上的信息。为了避免单点故障，用户可以安装辅助名称节点。建议有多个称为任务跟踪器节点/数据节点的从节点来执行实际的映射并减少进程。每个从节点充当任务跟踪器节点和数据节点。MapReduce 应用的性能与从节点的数量成正比。Hadoop 系统还执行自动故障恢复；如果其中一个任务跟踪器节点在运行时出现故障，那么它的责任将自动分配给另一个任务跟踪器，处理将继续进行，不会出现故障。

下面的程序演示了用 Python 开发一个基于 Hadoop 的 MapReduce 程序。它处理常见的爬网数据集。这些数据集包含长期收集的数千兆字节的网络爬行数据。它包含以**【Web ARChive】**(WARC)格式存储的网页数据、提取的元数据和提取的文本。这些数据集作为亚马逊公共数据集计划的一部分存储在亚马逊 S3。关于这些数据集的更多信息可在commoncrawl.org/the-data/ge…找到:

import sys
for line in sys.stdin:
  try:
    line = line.strip()
    # split the line into words
    words = line.split()
    # increase counters  
    if words[0] == "WARC-Target-URI:" :
      uri = words[1].split("/")
      print '%s\t%s' % (uri[0]+"//"+uri[2], 1)
  except Exception:
    print "There is some Error"

前面的程序是贴图部分，后面的程序是缩小部分:

from operator import itemgetter
import sys

current_word = None
current_count = 0
word = None

for line in sys.stdin:
    line = line.strip()

    word, count = line.split('\t', 1)

    try:
        count = int(count)
    except ValueError:
        continue

    if current_word == word:
        current_count += count
    else:
        if current_word:
            print '%s\t%s' % (current_word, current_count)
        current_count = count
        current_word = word

if current_word == word:
    print '%s\t%s' % (current_word, current_count)

在执行前一个程序之前，用户需要将名为web-crawl.txt的输入数据集文件转移到他们的HDFS home文件夹中。要执行该程序，用户应该运行以下命令:

#hadoop jar /usr/local/apache/hadoop2/share/hadoop/tools/lib/hadoop-streaming-2.6.0.jar -file /mapper.py    -mapper /mapper.py -file /reducer.py   -reducer /reducer.py -input /sample_crawl_data.txt -output /output

Python 中的火花

Spark 是一个通用的集群计算系统。它支持 Java、Python 和 Scala 的高级 API。这使得并行任务的编写变得容易。它是相对于 Hadoop 基于磁盘的 MapReduce 的两阶段模型而提出和开发的，因为它遵循内存模型，并为某些特定应用提供了最大 100%的性能提升。它非常适合实现机器学习应用/算法。

Spark 需要集群管理和分布式存储系统。它为各种分布式存储提供了一个简单的接口，包括亚马逊 S3、卡珊德拉、HDFS 等等。此外，它支持独立的——也就是说，用于集群管理的 spark 原生集群、Hadoop、纱和 Apache Mesos。

Spark Python API 命名为 PySpark ，将 Spark 编程模型暴露给 Python。我们可以在使用 PySpark 打开的 Python 外壳中或者通过使用 IPython 会话来开发基于 Spark 的应用。我们也可以先开发程序，然后使用pyspark命令运行。

总结

在本章中，我们讨论了使用 IPython 的高性能科学计算的概念。我们从并行计算的基本概念开始讨论。在基本概念之后，我们讨论了 IPython 并行计算的详细体系结构。后来，我们讨论了示例并行程序、IPython 魔法函数和并行装饰器的开发。

我们还介绍了 IPython 的高级功能:容错、动态负载平衡、管理任务之间的依赖关系、客户端和引擎之间的对象移动、IPython 数据库支持、使用来自 IPython 的 MPI 以及使用来自 IPython 的 StarCluster 管理 Amazon EC2 集群。然后我们讨论了 Python 中的多处理和多线程。最后，我们介绍了使用 Python 中的 Hadoop 和 Spark 开发分布式应用。

在下一章中，我们将讨论几个使用 Python 工具/API 进行科学计算的真实案例研究。我们将考虑各种基础和高级科学分支的应用。**