SciPy-数值与科学计算学习指南-全-

学习 SciPy 进行数值和科学计算

solve(A, b, sym_pos=False, lower=False, overwrite_a=False, overwrite_b=False, debug=False)
spsolve(A, b[, permc_spec, use_umfpack])

>>> from scipy import stats
>>> result=scipy.stats.bayes_mvs(scores)

% port search scipy

% port install py27-scipy

% python setup.py build –fcompiler=<compiler>

% python setup.py install

>>> import scipy
>>> scipy.test()

>>> import numpy
>>> scores = numpy.array([114, 100, 104, 89, 102, 91, 114, 114, 103, 105, 108, 130, 120, 132, 111, 128, 118, 119, 86, 72, 111, 103, 74, 112, 107, 103, 98, 96, 112, 112, 93])

>>> import scipy
>>> xmean = scipy.mean(scores)
>>> sigma = scipy.std(scores)
>>> n = scipy.size(scores)
>>> xmean, xmean - 2.576*sigma /scipy.sqrt(n), \
 xmean + 2.576*sigma / scipy.sqrt(n)

(105.83870967741936, 99.343223715529746, 112.33419563930897)

>>> from scipy import stats
>>> result=scipy.stats.bayes_mvs(scores)

>>> help(scipy.stats.bayes_mvs)

>>> result[0]

(105.83870967741936, (101.48825534263035, 110.18916401220837))

>>> import scipy.stats
>>> help(scipy.stats.bayes_mvs)

>>> import numpy
>>> numpy.info('random')

>>> import scipy.stats
>>> help(scipy.stats)

>>> help(scipy.stats.kurtosis)

>>> import numpy
>>> import matplotlib.pyplot as plt
>>> x=numpy.linspace(0,2*numpy.pi,32)
>>> fig = plt.figure()
>>> plt.plot(x, numpy.sin(x))
>>> plt.show()
>>> fig.savefig('sine.png')

ipython notebook filename.ipynb

>>> import scipy.misc
>>> img=scipy.misc.lena()
>>> import matplotlib.pyplot as plt 
>>> plt.gray() 
>>> plt.imshow(img) 

>>> plt.show() 

array([[162, 162, 162, 161, 162, 157, 163],
 [162, 162, 162, 161, 162, 157, 163],
 [162, 162, 162, 161, 162, 157, 163]])

>>> img[1,1:6]=0
>>> print (img[0:3,0:7])

[[162 162 162 161 162 157 163]
 [162   0   0   0   0   0 163]
 [162 162 162 161 162 157 163]]

>>> img.dtype, img.shape, img.size

(dtype('int64'), (512, 512), 262144)

>>> img[32,67]

87

>>> import scipy.misc
>>> img=scipy.misc.lena().astype('float32')

>>> import numpy
>>> scores = numpy.array([101,103,84], dtype='float32')
>>> scores

array([ 101.,  103.,   84.], dtype=float32)

>>> scores = numpy.float32([101,103,84])
>>> scores

array([ 101.,  103.,   84.], dtype=float32)

>>> a=numpy.array(['Cleese', 'Idle', 'Gilliam'], dtype='str_')
>>> a.dtype

dtype('<U7')

>>> dt = numpy.dtype([ ('name', numpy.str_, 16), ('grades',numpy.float64, (2,)) ])
>>> dt 

dtype([('name', '<U16'), ('grades', '<f8', (2,))])

>>> MA141=numpy.array([ ('Cleese', (7.0,8.0)), ('Gilliam',(9.0,10.0)) ], dtype=dt)
>>> MA141

array([('Cleese', [7.0, 8.0]), ('Gilliam', [9.0, 10.0])],dtype=[('name', '<U16'), ('grades', '<f8', (2,))])

>>> A=numpy.array([[1,2,3,4,5,6,7,8],[2,4,6,8,10,12,14,16]])
>>> print (A[0:2, 0:8:2])

[[ 1  3  5  7]
 [ 2  6 10 14]]

>>> print (A[0:2, 8:0:-2])

[[ 8,  6,  4,  2]
 [16, 12,  8,  4]]

>>> print (A[ [0,0,1,1], [0,3,2,5] ])

[ 1  4  6 12]

>>> print (A[ numpy.ix_( [0,1], [0,3] )])

[[1 4]
 [2 8]]

>>> value=0; import scipy.misc; img=scipy.misc.lena()
>>> for item in img.flat: value+=item
>>> value

32518120

>>> img.tofile("lena.txt",sep=" ",format="%i")

>>> import numpy
>>> A = numpy.array([11,13,15,17,19,18,16,14,12,10])
>>> A.argsort(kind='mergesort')

array([9, 0, 8, 1, 7, 2, 6, 3, 5, 4])

>>> A.sort()
>>> print(A)

[10 11 12 13 14 15 16 17 18 19]

>>> A=numpy.array([[1,1,1],[2,2,2],[3,3,3]])
>>> A.mean()

2

>>> A.mean(axis=0)

array([ 2.,  2.,  2.])

>>> A.mean(axis=1)

array([ 1.,  2.,  3.])

>>> img.min(), img.max(), img.ptp()

(25, 245, 220)

>>> A=img.clip(img.min(),img.min()+100)
>>> A.min(), A.max(), A.ptp()

(25, 125, 100)

>>> A=numpy.array([[1,2],[2,1]])

array(object,dtype=None,copy=True,order=None, subok=False,ndim=0)

>>> Z=numpy.zeros((5,5), dtype=int)
>>> U=numpy.ones((2,2), dtype=int)
>>> I=numpy.identity(3, dtype=int)

numpy.eye(N,M=None,k=0,dtype=float)

>>> D=numpy.eye(4,k=1) + numpy.eye(4,k=-1)
>>> print (D)

[[ 0\.  1\.  0\.  0.]
 [ 1\.  0\.  1\.  0.]
 [ 0\.  1\.  0\.  1.]
 [ 0\.  0\.  1\.  0.]]

>>> A=3.0*(numpy.eye(5,k=4) + numpy.eye(5,k=-4))
>>> A[0:2,0:2]=5*U; A[2:5,2:5]=6*I
>>> print (A)

[[ 5\.  5\.  0\.  0\.  3.]
 [ 5\.  5\.  0\.  0\.  0.]
 [ 0\.  0\.  6\.  0\.  0.]
 [ 0\.  0\.  0\.  6\.  0.]
 [ 3\.  0\.  0\.  0\.  6.]]

>>> B=numpy.fromfunction( (lambda i,j: i*j), (4,4), dtype=int)
>>> print (B)

[[0 0 0 0]
 [0 1 2 3]
 [0 2 4 6]
 [0 3 6 9]]

>>> print (B==0)

[[ True  True  True  True]
 [ True False False False]
 [ True False False False]
 [ True False False False]]

>>> B += numpy.fromfunction((lambda i,j:i*i+j*j), (4,4))*(B==0)
>>> print (B)

[[0 1 4 9]
 [1 1 2 3]
 [4 2 4 6]
 [9 3 6 9]]

>>> print (B%2!=0)

[[False  True False  True]
 [ True  True False  True]
 [False False False False]
 [ True  True False  True]]

>>> numpy.where(B%2!=0)

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

>>> numpy.putmask( B, B%2!=0, B**2+1)
>>> print (B)

[[ 0  2  4 82]
 [ 2  2  2 10]
 [ 4  2  4  6]
 [82 10  6 82]]

>>> L1=numpy.arange(-1,1,0.3)
>>> print (L1)

[-1\.  -0.7 -0.4 -0.1  0.2  0.5  0.8]

>>> L2=numpy.linspace(-1,1,4)
>>> print (L2)

[-1\.         -0.33333333  0.33333333  1\.        ]

>>> L3= numpy.logspace(-1,1,4)
>>> print (L3)

[  0.1          0.46415888   2.15443469  10\.        ]

>>> print (numpy.meshgrid(L2,L3))

(array([[-1\.        , -0.33333333,  0.33333333,  1\.        ],
 [-1\.        , -0.33333333,  0.33333333,  1\.        ],
 [-1\.        , -0.33333333,  0.33333333,  1\.        ],
 [-1\.        , -0.33333333,  0.33333333,  1\.        ]]), array([[ 
0.1       ,   0.1       ,   0.1       ,   0.1       ],
 [  0.46415888,   0.46415888,   0.46415888,   0.46415888],
 [  2.15443469,   2.15443469,   2.15443469,   2.15443469],
 [ 10\.        ,  10\.        ,  10\.        ,  10\.        ]]))

>>> print (numpy.mgrid[0:5,0:5])

[[[0 0 0 0 0]
 [1 1 1 1 1]
 [2 2 2 2 2]
 [3 3 3 3 3]
 [4 4 4 4 4]]

 [[0 1 2 3 4]
 [0 1 2 3 4]
 [0 1 2 3 4]
 [0 1 2 3 4]
 [0 1 2 3 4]]]

>>> print (numpy.ogrid[0:5,0:5])

[array([[0],
 [1],
 [2],
 [3],
 [4]]), array([[0, 1, 2, 3, 4]])]

tile(A, reps)

>>> checker2by2[0:3,0:3]=checker2by2[3:6,3:6]=B
>>> numpy.tile(checker2by2,(4,4))

>>> import scipy.ndimage
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> text = scipy.ndimage.imread('Chap_02_text_image.png')
>>> text = np.mean(text.astype(float)/255,-1)*2-1

>>> letterE = text[37:53,275:291]

>>> corr = scipy.ndimage.correlate(text,letterE)
>>> eLocations = (corr >= 0.95 * corr.max())

>>> CorrLocIndex = np.where(eLocations==True)
>>> x=CorrLocIndex[1]
>>> x 

array([ 283,  514,  583,  681,  722,  881,  929, 1023,   64,  188,  452, 
 504,  892,  921, 1059, 1087, 1102, 1133,  118,  547,  690, 1066, 
 1110,  330,  363,  519,  671,  913,  951, 1119,  120,  292,  441, 
 516,  557,  602,  649,  688,  717,  747,  783,  813,  988, 1016, 
 250,  309,  505,  691,  769,  876,  904, 1057,  224,  289,  470, 
 596,  626,  780, 1027,  112,  151,  203,  468,  596,  751,  817, 
 867,  203,  273,  369,  560,  599,  888, 1111,  159,  221,  260, 
 352,  427,  861,  901, 1034, 1146,  325,  506,  558]) 

>>> y=CorrLocIndex[0] 
>>> y 

array([ 45,  45,  45,  45,  45,  45,  45,  45,  74,  74,  74,  74,  74, 
 74,  74,  74,  74,  74, 103, 103, 103, 103, 103, 132, 132, 132, 
 132, 132, 132, 132, 161, 161, 161, 161, 161, 161, 161, 161, 161, 
 161, 161, 161, 161, 161, 190, 190, 190, 190, 190, 190, 190, 190, 
 219, 219, 219, 219, 219, 219, 219, 248, 248, 248, 248, 248, 248, 
 248, 248, 277, 277, 277, 277, 277, 277, 277, 306, 306, 306, 306, 
 306, 306, 306, 306, 306, 335, 335, 335]) 

>>> thefig=plt.figure()
>>> plt.subplot(211)
<matplotlib.axes._subplots.AxesSubplot object at 0x7fb9b2390110>
>>> plt.imshow(text, cmap=plt.cm.gray, interpolation='nearest')
<matplotlib.image.AxesImage object at 0x7fb9b1f29410>
>>> plt.axis('off')

(-0.5, 1199.5, 359.5, -0.5) 

>>> plt.subplot(212) 
<matplotlib.axes._subplots.AxesSubplot object at 0x7fb9b1f29890> 
>>> plt.imshow(text, cmap=plt.cm.gray, interpolation='nearest') 
<matplotlib.image.AxesImage object at 0x7fb9b1637e10> 
>>> plt.autoscale(False) 
>>> plt.plot(x,y,'wo',markersize=10) 
[<matplotlib.lines.Line2D object at 0x7fb9b1647290>] 
>>> plt.axis('off') 

(-0.5, 1199.5, 359.5, -0.5) 

>>> plt.show() 

>>> import numpy 
>>> B = numpy.ones((3,3)) 
>>> checker2by2 = numpy.zeros((6,6)) 
>>> checker2by2[0:3,0:3] = checker2by2[3:6,3:6] = B 
>>> print(checker2by2)

[[ 1\.  1\.  1\.  0\.  0\.  0.]
 [ 1\.  1\.  1\.  0\.  0\.  0.]
 [ 1\.  1\.  1\.  0\.  0\.  0.]
 [ 0\.  0\.  0\.  1\.  1\.  1.]
 [ 0\.  0\.  0\.  1\.  1\.  1.]
 [ 0\.  0\.  0\.  1\.  1\.  1.]]

>>> numpy.vsplit(checker2by2,3)

[array([[ 1.,  1.,  1.,  0.,  0.,  0.],
 [ 1.,  1.,  1.,  0.,  0.,  0.]]),
 array([[ 1.,  1.,  1.,  0.,  0.,  0.],
 [ 0.,  0.,  0.,  1.,  1.,  1.]]),
 array([[ 0.,  0.,  0.,  1.,  1.,  1.],
 [ 0.,  0.,  0.,  1.,  1.,  1.]])]

>>> a=numpy.array([-numpy.pi, numpy.pi])
>>> print (numpy.vstack((a, numpy.sin(a))))

[[ -3.14159265e+00   3.14159265e+00]
 [ -1.22464680e-16   1.22464680e-16]]

# function max100
>>> def max100(x):
 return(x)

>>> max100(a)

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

>>> numpy.vectorize(max100)(a)

array([100, 100])

>>> numpy.unique(A)

array([ 1, 2, 3, 4, 5, 10])

>>> numpy.bincount(A)

array([0, 4, 3, 2, 1, 2, 0, 0, 0, 0, 1])

>>> A=numpy.fromfunction((lambda i,j: (i+1)*(-1)**(i*j)), (4,4))
>>> print (A)

[[ 1\.  1\.  1\.  1.]
 [ 2\. -2\.  2\. -2.]
 [ 3\.  3\.  3\.  3.]
 [ 4\. -4\.  4\. -4.]]

>>> B=numpy.log2(A)
__main__:1: RuntimeWarning: invalid value encountered in log2
>>> print (B)

[[ 0\.         0\.         0\.         0\.       ]
 [ 1\.               nan  1\.               nan]
 [ 1.5849625  1.5849625  1.5849625  1.5849625]
 [ 2\.               nan  2\.               nan]]

>>> numpy.sum(B), numpy.nansum(B)

(nan, 12.339850002884624)

>>> import numpy
>>> vectorA = numpy.array([1,2,3,4,5,6,7])
>>> vectorA

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

>>> vectorB = vectorA[::-1].copy()
>>> vectorB

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

>>> vectorB[0]=123
>>> vectorB

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

>>> vectorA

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

>>> vectorB = vectorA[::-1].copy()
>>> vectorB

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

>>> vectorC = vectorA + vectorB
>>> vectorC

array([8, 8, 8, 8, 8, 8, 8])

>>> vectorD = vectorB - vectorA
>>> vectorD

array([ 6,  4,  2,  0, -2, -4, -6])

>>> dotProduct1 = numpy.dot(vectorA,vectorB)
>>> dotProduct1

84

>>> dotProduct2 = (vectorA*vectorB).sum()
>>> dotProduct2

84

>>> vectorA = numpy.array([5, 6, 7])
>>> vectorB = numpy.array([7, 6, 5])
>>> crossProduct = numpy.cross(vectorA,vectorB)
>>> crossProduct

array([-12,  24, -12])

>>> crossProduct = numpy.cross(vectorB,vectorA)
>>> crossProduct

array([ 12, -24,  12])

numpy.matrix(data=object, dtype=None, copy=True)

>>> A=numpy.matrix("1,2,3;4,5,6")
>>> A

matrix([[1, 2, 3],
 [4, 5, 6]])

>>> A=numpy.matrix([[1,2,3],[4,5,6]])
>>> A

matrix([[1, 2, 3],
 [4, 5, 6]])

>>> A=numpy.matrix([ [0,10,0,0,0], [0,0,20,0,0], [0,0,0,30,0], [0,0,0,0,40], [0,0,0,0,0] ])
>>> A

matrix([[ 0, 10,  0,  0,  0],
 [ 0,  0, 20,  0,  0],
 [ 0,  0,  0, 30,  0],
 [ 0,  0,  0,  0, 40],
 [ 0,  0,  0,  0,  0]])

>>> import numpy
>>> rows=numpy.array([0,1,2,3])
>>> cols=numpy.array([1,2,3,4])
>>> vals=numpy.array([10,20,30,40])

>>> import scipy.sparse
>>> A=scipy.sparse.coo_matrix( (vals,(rows,cols)) )
>>> print (A); print (A.todense())

 (0, 1)  10
 (1, 2)  20
 (2, 3)  30
 (3, 4)  40
[[  0\.  10   0\.   0\.   0.]
 [  0\.   0\.  20   0\.   0.]
 [  0\.   0\.   0\.  30   0.]
 [  0\.   0\.   0\.   0\.  40]]

>>> scipy.sparse.isspmatrix_coo(A)

True

>>> B=numpy.mat(numpy.ones((3,3)))
>>> W=numpy.mat(numpy.zeros((3,3)))
>>> print (numpy.bmat('B,W;W,B'))

[[ 1\.  1\.  1\.  0\.  0\.  0.]
 [ 1\.  1\.  1\.  0\.  0\.  0.]
 [ 1\.  1\.  1\.  0\.  0\.  0.]
 [ 0\.  0\.  0\.  1\.  1\.  1.]
 [ 0\.  0\.  0\.  1\.  1\.  1.]
 [ 0\.  0\.  0\.  1\.  1\.  1.]]

>>> a=numpy.array([[1,2],[3,4]])
>>> a*a

array([[ 1,  4],
 [ 9, 16]])

>>> A=numpy.mat(a)
>>> A*A

matrix([[ 7, 10],
 [15, 22]])

>>> b=numpy.array([[1,2,3],[3,4,5]])
>>> numpy.dot(a,b)

array([[ 7, 10, 13],
 [15, 22, 29]])

>>> numpy.multiply(A,A)

matrix([[ 1,  4],
 [ 9, 16]])

>>> a=numpy.arange(5); A=numpy.mat(a)
>>> a.shape, A.shape, a.transpose().shape, A.transpose().shape

((5,), (1, 5), (5,), (5, 1))

>>> import scipy.linalg
>>> A=scipy.linalg.hadamard(8)
>>> zero_sum_rows = (numpy.sum(A,0)==0)
>>> B=A[zero_sum_rows,:]
>>> print (B[0:3,:])

[[ 1 -1  1 -1  1 -1  1 -1]
 [ 1  1 -1 -1  1  1 -1 -1]
 [ 1 -1 -1  1  1 -1 -1  1]]

>>> A = numpy.matrix("1+1j, 2-1j; 3-1j, 4+1j")
>>> print (A.T); print (A.H)

[[ 1.+1.j  3.-1.j]
 [ 2.-1.j  4.+1.j]]
[[ 1.-1.j  3.+1.j]
 [ 2.+1.j  4.-1.j]]

>>> import numpy
>>> import scipy.linalg
>>> mu = 1/numpy.sqrt(2)
>>> A = numpy.matrix([[mu,0,mu],[0,1,0],[mu,0,-mu]])
>>> B = scipy.linalg.kron(A,A)

>>> print (B[:,0:-1:2])

[[ 0.5  0.5  0\.   0.5]
 [ 0\.   0\.   0\.   0\. ]
 [ 0.5 -0.5  0\.   0.5]
 [ 0\.   0\.   0\.   0\. ]
 [ 0\.   0\.   1\.   0\. ]
 [ 0\.  -0\.   0\.   0\. ]
 [ 0.5  0.5  0\.  -0.5]
 [ 0\.   0\.   0\.  -0\. ]
 [ 0.5 -0.5  0\.  -0.5]]

>>> A=numpy.matrix("1,1j;21,3")
>>> A; A*A; A**2

matrix([[  1.+0.j,   0.+1.j],
 [ 21.+0.j,   3.+0.j]])
matrix([[  1.+21.j,   0\. +4.j],
 [ 84\. +0.j,   9.+21.j]])
matrix([[  1.+21.j,   0\. +4.j],
 [ 84\. +0.j,   9.+21.j]])

>>> numpy.asarray(A); numpy.asarray(A)*numpy.asarray(A); numpy.asarray(A)**2

array([[  1.+0.j,   0.+1.j],
 [ 21.+0.j,   3.+0.j]])
array([[   1.+0.j,   -1.+0.j],
 [ 441.+0.j,    9.+0.j]])
array([[   1.+0.j,   -1.+0.j],
 [ 441.+0.j,    9.+0.j]])

>>> import numpy
>>> import scipy.linalg
>>> a=numpy.arange(0,2*numpy.pi,1.6)
>>> A = scipy.linalg.toeplitz(a)
>>> print (A)

[[ 0\.   1.6  3.2  4.8]
 [ 1.6  0\.   1.6  3.2]
 [ 3.2  1.6  0\.   1.6]
 [ 4.8  3.2  1.6  0\. ]]

>>> print (numpy.exp(A))

[[   1\.            4.95303242   24.5325302   121.51041752]
 [   4.95303242    1\.            4.95303242   24.5325302 ]
 [  24.5325302     4.95303242    1\.            4.95303242]
 [ 121.51041752   24.5325302     4.95303242    1\.        ]]

>>> print (scipy.linalg.expm(A))

[[ 1271.76972856   916.49316549   916.63015271  1271.70874469]
 [  916.49316549   660.86560972   660.5306514    916.63015271]
 [  916.63015271   660.5306514    660.86560972   916.49316549]
 [ 1271.70874469   916.63015271   916.49316549  1271.76972856]]

>>> import numpy
>>> import scipy.linalg
>>> x=10**100; y=9; v=numpy.matrix([x,y])
>>> scipy.linalg.norm(v,2)

1e+100

>>> numpy.sqrt(x*x+y*y)

Traceback (most recent call last)
 File "<stdin>", line 1, in <module>
AttributeError: 'long' object has no attribute 'sqrt'

eigvals(A, B=None, overwrite_a=False)

>>> import numpy
>>> import scipy.misc
>>> from scipy.linalg import svd
>>> import matplotlib.pyplot as plt
>>> img=scipy.misc.lena()
>>> U,s,Vh=svd(img)      # Singular Value Decomposition
>>> A = numpy.dot( U[:,0:32],  # use only 32 singular values
 numpy.dot( numpy.diag(s[0:32]),
 Vh[0:32,:]))
>>> plt.subplot(121,aspect='equal'); plt.imshow(img); plt.gray()
>>> plt.subplot(122,aspect='equal'); plt.imshow(A)
>>> plt.show()

solve(A, b, sym_pos=False, lower=False, overwrite_a=False, overwrite_b=False, debug=False)
spsolve(A, b[, permc_spec, use_umfpack])

lstsq(A, b, cond=None, overwrite_a=False, overwrite_b=False)

>>> import numpy
>>> import scipy.linalg
>>> A=numpy.mat(numpy.eye(3,k=1))
>>> print(A)

[[ 0\.  1\.  0.]
 [ 0\.  0\.  1.]
 [ 0\.  0\.  0.]]

>>> b=numpy.mat(numpy.arange(3) + 1).T
>>> print(b)

[[1]
 [2]
 [3]]

>>> xinfo=scipy.linalg.lstsq(A,b)
>>> print (xinfo[0].T)      # output the solution

[[ 0\.  1\.  2.]]

>>> import numpy
>>> import scipy.special
>>> a=scipy.special.exp10(-16)
>>> numpy.log(1+a)

0.0

>>> scipy.special.log1p(a)

9.9999999999999998e-17

>>> import numpy
>>> P1=numpy.poly1d([1,0,1])           # using coefficients
>>> print (P1)

 2
1 x + 1

>>> print (P1.r); print (P1.o); print (P1.deriv())

[ 0.+1.j  0.-1.j]
2
2 x

>>> P2=numpy.poly1d([1,1,1], True)     # using roots
>>> print (P2)

 3     2
1 x - 3 x + 3 x – 1

>>> P2=numpy.poly1d([1,1,1], True, variable='z')
>>> print (P2)

 3     2 
1 z - 3 z + 3 z - 1

>>> P1( numpy.arange(10) )           # evaluate at 0,1,...,9

array([ 1,  2,  5, 10, 17, 26, 37, 50, 65, 82])

>>> P1.__call__(numpy.arange(10))    # same evaluation

array([ 1,  2,  5, 10, 17, 26, 37, 50, 65, 82])

>>> import numpy
>>> Px=numpy.poly1d([-(1./2.),1,0])
>>> print(Px)

 2
-0.5 x + 1 x

>>> a=1./10000000000000000.
>>> print(a)

1e-16

>>> Px(a)

9.9999999999999998e-17

>>> P1=numpy.poly1d([1,0,1])
>>> print(P1)

 2
1 x + 1

>>> print(numpy.polyadd(P1, numpy.poly1d([2,1])))

 2
1 x + 2 x + 2

>>> print(numpy.polyadd(P1, [2,1]))

 2
1 x + 2 x + 2

>>> print(P1 + numpy.poly1d([2,1]))

 2
1 x + 2 x + 2

>>> print(P1 + [2,1])

 2
1 x + 2 x + 2

>>> P1/[2,1]

(poly1d([ 0.5 , -0.25]), poly1d([ 1.25]))

>>> eval_jacobi(n, alpha, beta, x)

>>> import numpy 
>>> import scipy.special
>>> import matplotlib.pyplot as plt
>>> x=numpy.linspace(-1,1,1000)
>>> plt.plot(x,scipy.special.eval_jacobi(3,0,1,x))
>>> plt.show()

>>> import scipy.special
>>> 10**10*scipy.special.psi(10**15)

345387763949.10681

>>> import numpy
>>> import scipy.special
>>> import  matplotlib.pyplot as plt
>>> import mpl_toolkits.mplot3d
>>> x=numpy.mgrid[-4:4:100j,-4:4:100j]
>>> z=x[0]+1j*x[1]
>>> (Ai, Aip, Bi, Bip) = scipy.special.airy(z)
>>> steps = range(int(Bi.real.min()), int(Bi.real.max()),6)
>>> fig=plt.figure()
>>> subplot1=fig.add_subplot(121,aspect='equal')
>>> subplot1.contourf(x[0], x[1], Bi.real, steps)
>>> subplot2=fig.add_subplot(122,projection='3d')
>>> subplot2.plot_surface(x[0],x[1],Bi.real)
>>> plt.show()

>>> import numpy
>>> import scipy.special
>>> scipy.special.jn(5,numpy.pi)

0.052141184367118461

>>> import numpy
>>> import matplotlib.pyplot as plt
>>> import scipy.interpolate
>>> x=numpy.linspace(-1,1,10); xn=numpy.linspace(-1,1,1000)
>>> y=numpy.sin(x)
>>> polynomial=scipy.interpolate.lagrange(x, numpy.sin(x))
>>> plt.plot(xn,polynomial(xn),x,y,'or')
>>> plt.show()

>>> import numpy
>>> import scipy.interpolate
>>> x1=numpy.linspace(1,10,10); y1=numpy.sin(x1)
>>> Polynomial=scipy.interpolate.BarycentricInterpolator(x1,y1)
>>> exactValues=numpy.sin(x1+0.3)
>>> exactValues

array([ 0.96355819,  0.74570521, -0.15774569, -0.91616594, 
-0.83226744,
 0.0168139 ,  0.85043662,  0.90217183,  0.12445442, 
-0.76768581])

>>> interpolatedValues=Polynomial(x1+0.3)
>>> interpolatedValues

array([ 0.97103132,  0.74460631, -0.15742869, -0.91631362, 
-0.83216445, 
 0.01670922,  0.85059283,  0.90181323,  0.12588718, 
-0.7825744 ])

>>> PercentRelativeError = numpy.abs((exactValues - interpolatedValues)/interpolatedValues)*100
>>> PercentRelativeError

array([ 0.76960822,  0.14758101,  0.20136334,  0.01611703,  0.01237594, 
 0.62647084,  0.01836479,  0.0397652 ,  1.13812858,  1.90251374])

>>> x2=numpy.linspace(1.5,10.5,10); y2=numpy.sin(x2)
>>> Polynomial.add_xi(x2,y2)
>>> interpolatedValues2=Polynomial(x1+0.3)
>>> interpolatedValues2

array([ 0.96355818,  0.74570521, -0.15774569, -0.91616594, -0.83226744,
 0.0168139 ,  0.85043662,  0.90217183,  0.12445442, -0.76768581])

>>> PercentRelativeError = numpy.abs((exactValues - interpolatedValues2)/interpolatedValues2)*100
>>> PercentRelativeError

array([  1.26241742e-07,   2.02502252e-09,   5.95225989e-10,
 1.84438143e-11,   8.75086862e-12,   4.14359323e-10,
 1.75194631e-11,   8.52321518e-11,   9.45285176e-09,
 1.29570657e-07])

>>> import numpy
>>> import matplotlib.pyplot as plt
>>> import  scipy.interpolate
>>> x=numpy.array([0,0,1,1,2,2]); y=numpy.array([0,0,1,0,2,0])
>>> interp=scipy.interpolate.KroghInterpolator(x,y)
>>> xn=numpy.linspace(0,2,20)   # evaluate polynomial in larger set
>>> plt.plot(x,y,'o',xn,interp(xn),'r')
>>> plt.show()

InterpolatedUnivariateSpline(x, y, w=None, bbox=[None, None], k=3)

>>> import numpy
>>> import matplotlib.pyplot as plt
>>>import scipy.interpolate
>>> x=numpy.arange(5); y=numpy.sin(x)
>>> xn=numpy.linspace(0,4,40)
>>> interp=scipy.interpolate.InterpolatedUnivariateSpline(x,y)
>>> plt.plot(x,y,'.',xn,interp(xn))
>>> plt.show()

>>> import numpy
>>> import scipy.interpolate
>>> import matplotlib.pyplot as plt
>>> from mpl_toolkits.mplot3d import Axes3D
>>> x=y=numpy.arange(10)
>>> f=(lambda i,j: numpy.sin(i)*numpy.cos(j))  # function to interpolate
>>> A=numpy.fromfunction(f, (10,10))           # generate samples
>>> spline=scipy.interpolate.RectBivariateSpline(x,y,A)
>>> fig=plt.figure()
>>> subplot=fig.add_subplot(111,projection='3d')
>>> xx=numpy.mgrid[0:9:100j, 0:9:100j]     # larger grid for plotting
>>> A=spline(numpy.linspace(0,9,100), numpy.linspace(0,9,100))
>>> subplot.plot_surface(xx[0],xx[1],A)
>>> plt.show()

>>> import numpy
>>> import scipy
>>> import matplotlib.pyplot as plt
>>> x=numpy.linspace(0,1,10)
>>> y=numpy.sin(x*numpy.pi/2)
>>> line=numpy.polyfit(x,y,deg=1)
>>> plt.plot(x,y,'.',x,numpy.polyval(line,x),'r')
>>> plt.show()

>>> import numpy
>>> import scipy.interpolate
>>> import matplotlib.pyplot as plt
>>> x=numpy.linspace(0,1,10)
>>> y=numpy.sin(x*numpy.pi/2)
>>> spline=scipy.interpolate.UnivariateSpline(x,y,k=2)
>>> xn=numpy.linspace(0,1,100)
>>> plt.plot(x,y,'.', xn, spline(xn))
>>> plt.show()

curve_fit(f, xdata, ydata, p0=None, sigma=None, **kw)

>>> import numpy
>>> import scipy
>>> A=18; w=3*numpy.pi; h=0.5
>>> x=numpy.linspace(0,1,100); y=A*numpy.sin(w*x+h)
>>> y += 4*((0.5-scipy.rand(100))*numpy.exp(2*scipy.rand(100)**2))

>>> import scipy.optimize
>>> p0 = [20, 2*numpy.pi, 1]
>>> target_function = lambda x,AA,ww,hh: AA*numpy.sin(ww*x+hh)

>>> pF,pVar = scipy.optimize.curve_fit(target_function, x, y, p0)

>>> print (pF)

[ 18.13799397   9.32232504   0.54808516]

leastsq(func, x0, args=(), Dfun=None, full_output=0,
        col_deriv=0, ftol=1.49012e-8, xtol=1.49012e-8,
        gtol=0.0, maxfev=0, epsfcn=0.0, factor=100, diag=None):

leastsq(error_function,p0,args=(x,y))

>>> import scipy.optimize
>>> scipy.optimize.fmin(scipy.optimize.rosen,[0,0])

Optimization terminated successfully.
 Current function value: 0.000000
 Iterations: 79
 Function evaluations: 146
array([ 1.00000439,  1.00001064])

minimize( fun, x0, args=(), method='BFGS', jac=None, hess=None, hessp=None, bounds=None, constraints=(),tol=None, callback=None, options=None)

>>> import scipy.special
>>> print (scipy.special.jn_zeros(4,3))

[  7.58834243  11.06470949  14.37253667]

root(fun, x0, args=(), method='hybr', jac=None, tol=None, callback=None, options=None)

>>> import numpy 
>>> import matplotlib.pyplot as plt 
>>> f=lambda x: [x[0]**2 - 2*x[0] - x[1] + 0.5, x[0]**2 + 4*x[1]**2 - 4]
>>> x,y=numpy.mgrid[-0.5:2.5:24j,-0.5:2.5:24j]
>>> U,V=f([x,y])
>>> plt.quiver(x,y,U,V,color='r', \
 linewidths=(0.2,), edgecolors=('k'), \
 headaxislength=5)
>>> plt.show()

>>> import scipy.optimize
>>> f=lambda x: [x[0]**2 - 2*x[0] - x[1] + 0.5, x[0]**2 + 4*x[1]**2 - 4]
>>> scipy.optimize.root(f,[0,1])

 status: 1
 success: True
qtf: array([ -4.81190247e-09,  -3.83395899e-09])
nfev: 9
 r: array([ 2.38128242, -0.60840482, -8.35489601])
 fun: array([  3.59529073e-12,   3.85025345e-12])
 x: array([-0.22221456,  0.99380842])
 message: 'The solution converged.'
fjac: array([[-0.98918813, -0.14665209],
 [ 0.14665209, -0.98918813]])

>>> scipy.optimize.root(f,[2,0])

 status: 1
 success: True
qtf: array([  2.08960516e-10,   8.61298294e-11])
nfev: 12
 r: array([-4.56575336, -1.67067665, -1.81464307])
 fun: array([  2.44249065e-15,   1.42996726e-13])
 x: array([ 1.90067673,  0.31121857])
 message: 'The solution converged.'
fjac: array([[-0.39612596, -0.91819618],
 [ 0.91819618, -0.39612596]])

quad(func, a, b, args=(), full_output=0, epsabs=1.49e-08, epsrel=1.49e-08, limit=50, points=None, weight=None, wvar=None, wopts=None, maxp1=50, limlst=50)

>>> simps(y, x=None, dx=1, axis=-1, even='avg')

>>> f=lambda t: numpy.exp(-t)*t**4
>>> from scipy.special import gammainc
>>> from scipy.integrate import quad
>>> from scipy.misc import factorial
>>> print (gammainc(5,1))

0.00365984682734

print('%.19f' % gammainc(5,1))

0.0036598468273437131

>>> import numpy
>>> result,error=quad(f,0,1)/factorial(4)
>>> result

0.0036598468273437122

>>> import numpy 
>>> import scipy.integrate
>>> x=numpy.linspace(0,1,10000)
>>> scipy.integrate.simps(f(x)/factorial(4), x)

0.0036598468273469071

ode(f,jac=None)

>>> import numpy
>>> from scipy.integrate import ode
>>> f=lambda t,y: -20*y        # The ODE
>>> actual_solution=lambda t:numpy.exp(-20*t)  # actual solution
>>> dt=0.01            # time step
>>> solver=ode(f).set_integrator('dop853')  # solver
>>> solver.set_initial_value(1,0)      # initial value
>>> while solver.successful() and solver.t<=1+dt:
 # solve the equation at succesive time steps,
 # until the time is greater than 1
 # but make sure that the solution is successful
 print (solver.t, solver.y, actual_solution(solver.t))
 # We compare each numerical solution with the actual
 # solution of the ODE
 solver.integrate(solver.t + dt)    # solve next step

<scipy.integrate._ode.ode at 0x10eac5e50>
0 [ 1.] 1.0
0.01 [ 0.81873075] 0.818730753078
0.02 [ 0.67032005] 0.670320046036
0.03 [ 0.54881164] 0.548811636094
0.04 [ 0.44932896] 0.449328964117
0.05 [ 0.36787944] 0.367879441171
0.06 [ 0.30119421] 0.301194211912
0.07 [ 0.24659696] 0.246596963942
0.08 [ 0.20189652] 0.201896517995
0.09 [ 0.16529889] 0.165298888222
0.1 [ 0.13533528] 0.135335283237
 ...
0.9 [  1.52299797e-08] 1.52299797447e-08
0.91 [  1.24692528e-08] 1.24692527858e-08
0.92 [  1.02089607e-08] 1.02089607236e-08
0.93 [  8.35839010e-09] 8.35839010137e-09
0.94 [  6.84327102e-09] 6.84327102222e-09
0.95 [  5.60279644e-09] 5.60279643754e-09
0.96 [  4.58718175e-09] 4.58718174665e-09
0.97 [  3.75566677e-09] 3.75566676594e-09
0.98 [  3.07487988e-09] 3.07487987959e-09
0.99 [  2.51749872e-09] 2.51749871944e-09
1.0   [  2.06115362e-09] 2.06115362244e-09

odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0)

>>> import numpy
>>> from numpy import linspace
>>> import scipy
>>> from scipy.integrate import odeint
>>> import matplotlib.pyplot as plt
>>> from mpl_toolkits.mplot3d import Axes3D
>>> sigma=10.0
>>> b=8/3.0
>>> r=28.0
>>> f = lambda x,t: [sigma*(x[1]-x[0]), r*x[0]-x[1]-x[0]*x[2], x[0]*x[1]-b*x[2]]

>>> t=linspace(0,20,2000); y0=[5.0,5.0,5.0]
>>> solution=odeint(f,y0,t)
>>> X=solution[:,0]; Y=solution[:,1]; Z=solution[:,2]

>>> import matplotlib.pyplot as plt
>>> plt.gca(projection='3d'); plt.plot(X,Y,Z)
>>> plt.show()

>>> plt.rcParams['figure.figsize'] = (10.0, 5.0)
>>> plt.subplot(121); plt.plot(t,Z)
>>> plt.subplot(122); plt.plot(Y,Z)
>>> plt.show()

fft(x[, n, axis, overwrite_x])

>>> import numpy
>>> from scipy.fftpack import fft,fft2, fftshift
>>> import matplotlib.pyplot as plt
>>> B=numpy.ones((4,4)); W=numpy.zeros((4,4))
>>> signal = numpy.bmat("B,W;W,B")
>>> onedimfft = fft(signal,n=16)
>>> twodimfft = fft2(signal,shape=(16,16))
>>> plt.figure()
>>> plt.gray()
>>> plt.subplot(121,aspect='equal')
>>> plt.pcolormesh(onedimfft.real)
>>> plt.colorbar(orientation='horizontal')
>>> plt.subplot(122,aspect='equal')
>>> plt.pcolormesh(fftshift(twodimfft.real))
>>> plt.colorbar(orientation='horizontal')
>>> plt.show()

>>> import numpy
>>> from scipy.signal import chirp, sawtooth, square, gausspulse
>>> import matplotlib.pyplot as plt
>>> t=numpy.linspace(-1,1,1000)
>>> plt.subplot(221); plt.ylim([-2,2])
>>> plt.plot(t,chirp(t,f0=100,t1=0.5,f1=200))   # plot a chirp
>>> plt.title("Chirp signal")
>>> plt.subplot(222); plt.ylim([-2,2])
>>> plt.plot(t,gausspulse(t,fc=10,bw=0.5))      # Gauss pulse
>>> plt.title("Gauss pulse")
>>> plt.subplot(223); plt.ylim([-2,2])
>>> t*=3*numpy.pi
>>> plt.plot(t,sawtooth(t))                     # sawtooth
>>> plt.xlabel("Sawtooth signal")
>>> plt.subplot(224); plt.ylim([-2,2])
>>> plt.plot(t,square(t))                       # Square wave
>>> plt.xlabel("Square signal")
>>> plt.show()

fromfile(file, dtype=float, count=-1, sep='')

imread(fname, flatten=False)

>>> rate,data = scipy.io.wavfile.read("audio.wav")

>>> scipy.io.wavfile.write("filename.wav",rate,data)

>>> from numpy import sin, cos, pi, linspace
>>> f=lambda t: cos(pi*t) + 0.2*sin(5*pi*t+0.1) + 0.2*sin(30*pi*t) + 0.1*sin(32*pi*t+0.1) + 0.1*sin(47* pi*t+0.8)
>>> t=linspace(0,4,400); signal=f(t)

>>> from scipy.signal import wiener, medfilt
>>> import matplotlib.pylab as plt
>>> plt.plot(t,signal,'k', label='The signal')
>>> plt.plot(t,wiener(signal,mysize=55),'r',linewidth=3, label='Wiener filtered')
>>> plt.plot(t,medfilt(signal,kernel_size=55),'b',linewidth=3, label='Medfilt filtered')
>>> plt.legend()
>>> plt.show()

>>> from scipy.stats import norm     # Gaussian distribution
>>> import matplotlib.pyplot as plt
>>> import scipy.misc
>>> import scipy.ndimage
>>> plt.gray()
>>> lena=scipy.misc.lena().astype(float)
>>> plt.subplot(221);
>>> plt.imshow(lena)
>>> lena+=norm(loc=0,scale=16).rvs(lena.shape)
>>> plt.subplot(222);
>>> plt.imshow(lena)
>>> denoised_lena = scipy.ndimage.median_filter(lena,3)
>>> plt.subplot(224); 
>>> plt.imshow(denoised_lena)

sobel(image, axis=-1, output=None, mode='reflect', cval=0.0)

>>> from scipy.ndimage.filters import sobel
>>> import numpy
>>> lena=scipy.misc.lena()
>>> sblX=sobel(lena,axis=0); sblY=sobel(lena,axis=1)
>>> sbl=numpy.hypot(sblX,sblY)
>>> plt.subplot(223); 
>>> plt.imshow(sbl) 
>>> plt.show()

>>> from scipy.signal import boxcar, hamming, nuttall, triang
>>> import matplotlib.pylab as plt
>>> windows=['boxcar', 'hamming', 'nuttall', 'triang']
>>> plt.subplot(121)
>>> for w in windows:
 eval( 'plt.plot(' + w + '(31))' )
 plt.ylim([-0.5,2]); plt.xlim([-1,32])
 plt.legend(windows)

>>> plt.subplot(122)
>>> import numpy
>>> from numpy import sin, cos, pi, linspace
>>> f=lambda t: cos(pi*t) + 0.2*sin(5*pi*t+0.1) + 0.2*sin(30*pi*t) + 0.1*sin(32*pi*t+0.1) + 0.1*sin(47* pi*t+0.8)
>>> t=linspace(0,4,400); signal=f(t)
>>> extended_signal=numpy.r_[signal[15:0:-1],signal,signal[-1:-15:- 1]]
>>> plt.plot(extended_signal,'k')

>>> for w in windows:
 window = eval( w+'(31)')
 output=numpy.convolve(window/window.sum(),signal)
 plt.plot(output,linewidth=2)
 plt.ylim([-2,3]); plt.legend(['original']+windows)
>>> plt.show()

>>> import scipy.misc
>>> import numpy 
>>> import matplotlib.pylab as plt 
>>> from scipy.ndimage.interpolation import affine_transform
>>> lena=scipy.misc.lena() 
>>> A=numpy.mat("0,1;-1,1.25"); b=[-400,0]
>>> Ab_Lena=affine_transform(lena,A,b,output_shape=(512*2.2,512*2.2))
>>> plt.gray() 
>>> plt.subplot(121) 
>>> plt.imshow(Ab_Lena)

geometric_transform(input, mapping, output_shape=None, 
                    output=None, order=3, mode='constant',
cval=0.0, prefilter=True, extra_arguments=(),
extra_keywords={})

>>> def f(z):
 temp = a*(z[0]+1j*z[1]) + b
 temp /= c*(z[0]+1j*z[1])+d
 return (temp.real, temp.imag)

>>> from scipy.ndimage.interpolation import geometric_transform 
>>> a = 2**15*(1+1j); b = 0; c = -2**8*(1-1j*2); d = 2**18-1j*2**14
>>> Moebius_Lena = geometric_transform(lena,f,mode='reflect')
>>> plt.subplot(122); 
>>> plt.imshow(Moebius_Lena) 
>>> plt.show()

map_coordinates(input, coordinates, output=None, order=3, 
                mode='constant', cval=0.0, prefilter=True)

>>> import scipy.misc 
>>> from scipy.ndimage.interpolation import map_coordinates
>>> lena=scipy.misc.lena().astype(float)
>>> coordinates=[[10.5, 12.3], [11.7, 1.4]]
>>> map_coordinates(lena, coordinates, order=1)

array([ 157.2 ,  157.42])

>>> map_coordinates(lena, coordinates, order=2)

array([ 157.80641507,  157.6741489 ])

binary_operation(signal, structuring_element)

>>> binary_hit_or_miss(text, letterE)

>>> import numpy
>>> import scipy.ndimage
>>> import matplotlib.pylab as plt
>>> from scipy.ndimage.morphology import binary_hit_or_miss
>>> text = scipy.ndimage.imread('CHAP_05_input_textImage.png')
>>> letterE = text[37:53,275:291]
>>> HitorMiss = binary_hit_or_miss(text, structure1=letterE, origin1=1) 
>>> eLocation = numpy.where(HitorMiss==True)
>>> x=eLocation[1]; y=eLocation[0]
>>> plt.imshow(text, cmap=plt.cm.gray, interpolation='nearest')
>>> plt.autoscale(False)
>>> plt.plot(x,y,'wo',markersize=10)
>>> plt.axis('off')
>>> plt.show()

grey_operation(signal, [structuring_element, footprint, size, ...])

>>> grey_dilation(lena, size=(15,15))

distance_transform_bf(signal, metric='euclidean', sampling=None,
return_distances=True, return_indices=False,
                      distances=None, indices=None)

>>> import scipy.misc 
>>> from scipy.stats import signaltonoise 
>>> from scipy.stats import norm     # Gaussian distribution
>>> lena=scipy.misc.lena().astype(float)
>>> lena+= norm.rvs(loc=0,scale=16,size=lena.shape)
>>> signaltonoise(lena,axis=None)

array(2.459233897516763)

>>> import numpy
>>> from scipy.stats import pareto
>>> import matplotlib.pyplot as plt
>>> x=numpy.linspace(1,10,1000)
>>> plt.subplot(131); plt.plot(pareto.pdf(x,5))
>>> plt.subplot(132); plt.plot(pareto.cdf(x,5))
>>> plt.subplot(133); plt.plot(pareto.rvs(5,size=1000))
>>> plt.show()

bayes_mvs(data, alpha=0.9)

>>> import numpy
>>> data = numpy.array([[113,105,130,101,138,118,87,116,75,96, \
 122,103,116,107,118,103,111,104,111,89,78,100,89,85,88], \
 [137,105,133,108,115,170,103,145,78,107, \
 84,148,147,87,166,146,123,135,112,93,76,116,78,101,123]])

>>> dataDiff = data[1,:]-data[0,:]
>>> dataDiff.mean(), dataDiff.std()

(13.32, 22.472596645692729)

>>> import matplotlib.pyplot as plt
>>> plt.hist(dataDiff)
>>> plt.show()

>>> from scipy.stats import ttest_1samp
>>> t_stat,p_value=ttest_1samp(dataDiff,0.0)

>>> print (p_value/2.0)

0.00389575522747

>>> from scipy.stats import norm     # Gaussian distribution
>>> mean,std=norm.fit(dataDiff)

>>> plt.hist(dataDiff, normed=1)
>>> x=numpy.linspace(dataDiff.min(),dataDiff.max(),1000)
>>> pdf=norm.pdf(x,mean,std)
>>> plt.plot(x,pdf)
>>> plt.show()

>>> from scipy.stats import gaussian_kde
>>> pdf=gaussian_kde(dataDiff)

>>> from scipy.stats import gaussian_kde
>>> pdf = gaussian_kde(dataDiff)
>>> pdf = pdf.evaluate(x)
>>> plt.hist(dataDiff, normed=1)
>>> plt.plot(x,pdf,'k')
>>> plt.savefig("hist2.png")
>>> plt.show()

>>> plt.hist(dataDiff, normed=1)
>>> plt.plot(x,pdf,'k.-',label='Kernel fit')
>>> plt.plot(x,norm.pdf(x,mean,std),'r',label='Normal fit')
>>> plt.legend() 
>>> plt.savefig("hist3.png")
>>> plt.show()

distance_function(first_vector, second_vector)

>>> import numpy 
>>> from scipy.spatial.distance import minkowski 
>>> Square=numpy.mgrid[-1.1:1.1:512j,-1.1:1.1:512j]
>>> X=Square[0]; Y=Square[1]
>>> f=lambda x,y,p: minkowski([x,y],[0.0,0.0],p)<=1.0
>>> Ball=lambda p:numpy.vectorize(f)(X,Y,p)

>>> import matplotlib.pylab as plt
>>> plt.imshow(Ball(3), cmap = plt.cm.gray)
>>> plt.axis('off')
>>> plt.subplots_adjust(left=0.0127,bottom=0.0164,\
 right=0.987,top=0.984)
>>> plt.show()

>>> import scipy.stats
>>> from scipy.spatial.distance import pdist
>>> V=scipy.stats.randint.rvs(0.4,3,size=(5,4))-1
>>> print (V)

[[ 1  0  1 -1]
 [-1  0 -1  0]
 [ 1  1  1 -1]
 [ 1  1 -1  0]
 [ 0  0  1 -1]]

>>> pdist(V,metric='cityblock')

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

>>> from scipy.spatial.distance import cdist
>>> V=scipy.stats.randint.rvs(0.4, 2, size=(3,4)).astype(bool)
>>> W=scipy.stats.randint.rvs(0.4, 3, size=(2,4)).astype(bool)
>>> cdist(V,W,'jaccard')
array([[ 0.75      ,  1\.        ],
 [ 0.75      ,  1\.        ],
 [ 0.33333333,  0.5       ]])

KDTree(data, leafsize=10)

>>> from scipy.spatial import KDTree
>>> data=scipy.stats.randint.rvs(0.4,10,size=(10,4))
>>> print (data)

[[8 6 1 1]
 [2 9 1 5]
 [4 8 8 9]
 [2 6 6 4]
 [4 1 2 1]
 [3 8 7 2]
 [1 1 3 6]
 [5 2 1 5]
 [2 5 7 3]
 [6 0 6 9]]
>>> tree=KDTree(data)
>>> tree.query([0,0,0,0])

(4.6904157598234297, 4)

kmeans(obs, k_or_guess, iter=20, thresh=1e-05)

kmeans2(data, k, iter=10, thresh=1e-05,
minit='random', missing='warn')

vq(obs, centroids)

>>> import numpy
>>> from scipy.stats import norm
>>> from numpy import array,vstack
>>> data=norm.rvs(0,0.3,size=(10000,2))
>>> inside_ball=numpy.hypot(data[:,0],data[:,1])<1.0
>>> data=data[inside_ball]
>>> data = vstack((data, data+array([1,1]),data+array([-1,1])))

>>> from scipy.cluster.vq import *
>>> centroids, distortion = kmeans(data,3)
>>> cluster_assignment, distances = vq(data,centroids)

>>> from matplotlib.pyplot import plot
>>> import matplotlib.pyplot as plt
>>> plt.plot(data[cluster_assignment==0,0], \
 data[cluster_assignment==0,1], 'ro')
>>> plt.plot(data[cluster_assignment==1,0], \
 data[cluster_assignment==1,1], 'b+')
>>> plt.plot(data[cluster_assignment==2,0], \
 data[cluster_assignment==2,1], 'k.')
>>> plt.show()

linkage(data, method='single', metric='euclidean')

dendrogram(Z, p=30, truncate_mode=None, color_threshold=None, 
get_leaves=True, orientation='top', labels=None, 
count_sort=False, distance_sort=False, 
show_leaf_counts=True, no_plot=False, no_labels=False, 
color_list=None, leaf_font_size=None, 
leaf_rotation=None, leaf_label_func=None, 
no_leaves=False, show_contracted=False,
link_color_func=None)

OPOSSUM                    54113344
HAIRY TAIL MOLE            33114433
COMMON MOLE              32103333
STAR NOSE MOLE            33114433

>>> import numpy
>>> file=open("dentitio.dat","r") # open the file
>>> lines=file.readlines() # read each line in memory
>>> file.close() # close the file
>>> mammals=[] # this stores the names
>>> dataset=numpy.zeros((len(lines),8)) # this stores the data
>>> for index,line in enumerate(lines):
 mammals.append( line[0:27].rstrip(" ").capitalize() )
 for tooth in range(8):
 dataset[index,tooth]=int(line[27+tooth])

>>> import matplotlib.pyplot as plt
>>> from scipy.cluster.hierarchy import linkage, dendrogram
>>> Z=linkage(dataset)
>>> dendrogram(Z, labels=mammals, orientation="right")
>>> matplotlib.pyplot.show()
>>> plt.show()

>>> import scipy.stats 
>>> import scipy.spatial 
>>> data = scipy.stats.randint.rvs(0.4,10,size=(10,2))
>>> triangulation = scipy.spatial.Delaunay(data)

>>> from numpy.random import RandomState
>>> rv = RandomState(123456789)
>>> locations = rv.randint(0, 511, size=(2,8))
>>> triangulation=scipy.spatial.Delaunay(locations.T)

>>> import matplotlib.pyplot as plt 
>>> assign_vertex = lambda index: triangulation.points[index]
>>> triangle_set = map(assign_vertex, triangulation.vertices)

>>> voronoiSet=scipy.spatial.Voronoi(locations.T)
>>> scipy.spatial.voronoi_plot_2d(voronoiSet)
>>> fig = plt.figure()
>>> thefig = plt.subplot(1,1,1)
>>> scipy.spatial.voronoi_plot_2d(voronoiSet, ax=thefig)
>>> plt.triplot(locations[1], locations[0], triangles=triangle_set, color='r')

>>> plt.xlim((0,550))

 (0, 550)

>>> plt.ylim((0,550))

 (0, 550)

>>> plt.show()

>>> import numpy
>>> import scipy
>>> from scipy.ndimage import *
>>> from scipy.misc import imfilter
>>> import matplotlib.pyplot as plt
>>> import matplotlib.cm as cm

>>> img=imread('./NbW-STEM.png')
>>> minVal = numpy.min(img) 
>>> maxVal = numpy.max(img) 
>>> img = (1.0/(maxVal-minVal))*(img - minVal) 
>>> plt.imshow(img, cmap = cm.Greys_r)
>>> plt.show()
>>> print "Image dtype: %s"%(img.dtype)
>>> print "Image size: %6d"%(img.size)
>>> print "Image shape: %3dx%3d"%(img.shape[0],img.shape[1])
>>> print "Max value %1.2f at pixel %6d"%(img.max(),img.argmax())
>>> print "Min value %1.2f at pixel %6d"%(img.min(),img.argmin())
>>> print "Variance: %1.5f\nStandard deviation: \ 
 %1.5f"%(img.var(),img.std())

Image dtype: float64
Image size:  87025
Image shape: 295x295
Max value 1.00 at pixel  75440
Min value 0.00 at pixel   5703
Variance: 0.02580
Standard deviation: 0.16062

>>> plt.subplot(1, 2, 1)
>>> plt.imshow(img > 0.2, cmap = cm.Greys_r)
>>> plt.xlabel('img > 0.2')
>>> plt.subplot(1, 2, 2) 
>>> plt.imshow(img > 0.7, cmap = cm.Greys_r)
>>> plt.xlabel('img > 0.7')
>>> plt.show()

>>> BWatoms = (img> 0.62)
>>> BWatoms = binary_opening(BWatoms,structure=numpy.ones((2,2)))

>>> structuring_element = [[0,1,0],[1,1,1],[0,1,0]]
>>> segmentation,segments = label(BWatoms,structuring_element)

>>> coords = center_of_mass(img, segmentation, range(1,segments+1))
>>> xcoords = numpy.array([x[1] for x in coords])
>>> ycoords = numpy.array([x[0] for x in coords])

>>> plt.imshow(img, cmap = cm.Greys_r) 
>>> plt.axis('off') 
>>> plt.plot(xcoords,ycoords,'b.') 
>>> plt.show() 

>>> L1,L2 = distance_transform_edt(segmentation==0, return_distances=False, return_indices=True)
>>> Voronoi = segmentation[L1,L2]
>>> Voronoi_edges= imfilter(Voronoi,'find_edges')
>>> Voronoi_edges=(Voronoi_edges>0)

>>> plt.imshow(Voronoi_edges); plt.axis('off'); plt.gray()
>>> plt.plot(xcoords,ycoords,'r.',markersize=2.0)
>>> plt.show()

>>> import numpy
>>> from numpy import linspace
>>> import scipy
>>> import matplotlib.pyplot as plt
>>> from scipy.spatial import Delaunay

>>> xmin = 0 ; xmax = 1 ; nXpoints = 10
>>> ymin = 0 ; ymax = 1 ; nYpoints = 10
>>> horizontal = linspace(xmin,xmax,nXpoints)
>>> vertical = linspace(ymin,ymax,nYpoints)
>>> y, x = numpy.meshgrid(horizontal, vertical)
>>> vertices = numpy.array([x.flatten(),y.flatten()])

>>> triangulation = Delaunay(vertices.T)
>>> index2point = lambda index: triangulation.points[index]
>>> all_centers = index2point(triangulation.vertices).mean(axis=1)
>>> trngl_set=triangulation.vertices

>>> plt.triplot(vertices[0],vertices[1],triangles=trngl_set)
>>> plt.show()

>>> from numpy import  cross 
>>> from scipy.sparse import dok_matrix 
>>> points=triangulation.points.shape[0]
>>> stiff_matrix=dok_matrix((points,points))
>>> for triangle in triangulation.vertices:
 helper_matrix=dok_matrix((points,points))
 pt1,pt2,pt3=index2point(triangle)
 area=abs(0.5*cross(pt2-pt1,pt3-pt1))
 coeffs=0.5*numpy.vstack((pt2-pt3,pt3-pt1,pt1-pt2))/area
 #helper_matrix[triangle,triangle] = \ 
 array(mat(coeffs)*mat(coeffs).T)
 u=None 
 u=numpy.array(numpy.mat(coeffs)*numpy.mat(coeffs).T) 
 for i in range(len(triangle)):
 for j in range(len(triangle)):
 helper_matrix[triangle[i],triangle[j]] = u[i,j] 
 stiff_matrix=stiff_matrix+helper_matrix

>>> allNodes = numpy.unique(trngl_set) 
>>> boundaryNodes = numpy.unique(triangulation.convex_hull) 
>>> NonBoundaryNodes = numpy.array([]) 
>>> for x in allNodes: 
 if x not in boundaryNodes: 
 NonBoundaryNodes = numpy.append(NonBoundaryNodes,x) 
 NonBoundaryNodes = NonBoundaryNodes.astype(int) 
 nbnodes = len(boundaryNodes) # number of boundary nodes 
 FbVals=numpy.zeros([nbnodes,1]) # Values on the boundary 
 FbVals[(nbnodes-nXpoints+1):-1]=numpy.ones([nXpoints-2, 1])

>>> totalNodes = len(allNodes) 
>>> stiff_matrixDense = stiff_matrix.todense() 
>>> stiffNonb = \ 
 stiff_matrixDense[numpy.ix_(NonBoundaryNodes,NonBoundaryNodes)] 
>>> stiffAtb = \ 
 stiff_matrixDense[numpy.ix_(NonBoundaryNodes,boundaryNodes)] 
>>> U=numpy.zeros([totalNodes, 1]) 
>>> U[NonBoundaryNodes] = numpy.linalg.solve( - stiffNonb , \
 stiffAtb * FbVals ) 
>>> U[boundaryNodes] = FbVals 

>>> X = vertices[0] 
>>> Y = vertices[1] 
>>> Z = U.T.flatten() 
>>> from mpl_toolkits.mplot3d import axes3d
>>> fig = plt.figure() 
>>> ax = fig.add_subplot(111, projection='3d') 
>>> surf = ax.plot_trisurf(X, Y, Z, cmap=cm.jet, linewidth=0) 
>>> fig.colorbar(surf) 
>>> fig.tight_layout() 
>>> ax.set_xlabel('X',fontsize=16)
>>> ax.set_ylabel('Y',fontsize=16)
>>> ax.set_zlabel(r"$\phi$",fontsize=36)
>>> plt.show() 

>>> from numpy import pi, sinh, sin, cos, sum
>>> def f(x,y): 
 return sum( 2*(1.0/(n*pi) - \
 cos(n*pi)/(n*pi))*(sinh(n*pi*x)/ \
 sinh(n*pi))*sin(n*pi*y) 
 for n in range(1,200)) 
>>> Ze = f(X,Y) 
>>> ZdiffZe = Ze - Z 
>>> plt.plot(ZdiffZe) 
>>> plt.show() 

SUBROUTINE PRIMEFACTORS(num, factors, f)
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: num  !input number
  INTEGER,INTENT(OUT), DIMENSION((num/2))::factors 
  INTEGER, INTENT(INOUT) :: f
  INTEGER :: i, n
  i = 2  
  f = 1  
  n = num
  DO
    IF (MOD(n,i) == 0) THEN 
      factors(f) = i
      f = f+1
      n = n/i
    ELSE
      i = i+1
    END IF
    IF (n == 1) THEN    
      f = f-1    
      EXIT
    END IF
  END DO

% f2py –c primefactors.f90 –m primefactors

>>> import primefactors
>>> primefactors.primefactors(6,1)

array([2, 3, 0], dtype=int32)

inline(code, arg_names, local_dict=None, global_dict=None,
           force = 0,
           compiler='',
           verbose = 0,
support_code = None,
           customize=None,
type_factories = None,
auto_downcast=1,
           **kw)

>>> import scipy.weave
>>> name = 'Francisco'
>>> pin = 1234
>>> code = 'std::cout << name << "---PIN: " '
>>> code+= '<<std::hex << pin <<std::endl;'
>>> arg_names = ['name','pin']
>>> scipy.weave.inline(code, arg_names)

Francisco---PIN: 4d2

% ./configure --enable-R-static-lib --enable-static --with-readline=no

% cd src/nmath/standalone
% make

>>> import scipy.weave 
>>> support_code= 'extern "C" double pbinom(double x,\ 
 double n, double p, int lower_tail, int log_p);' 
>>> libpath='/opt/Rlib' #IS THE LOCATION OF LIBRARY libRmath.so
>>> library_dirs=[libpath] 
>>> libraries=['Rmath'] 
>>> runtime_library_dirs=[libpath] 
>>> code='return_val=pbinom(100,20000,100./20000.,0,1);' 
>>> res=scipy.weave.inline(code, support_code=support_code, \ 
 library_dirs=library_dirs, libraries=libraries, \ 
 runtime_library_dirs=runtime_library_dirs) 
>>> print(res) 

-0.747734910363 

>>> headers = ['<math.h>']

octave:1> x=rand(1,10);
octave:2> y=rand(size(x));
octave:3> T=delaunay(x,y);
octave:4> save -v7 data x y T

>>> from scipy.io import loadmat
>>> datadict = loadmat("data")

>>> datadict.keys()

['__header__', '__globals__', 'T', 'y', 'x', '__version__']

>>> datadict['x']

array([[0.81222999,0.51836246,0.60425982,0.23660352,0.01305779,
 0.0875166,0.77873049,0.70505801,0.51406693,0.65760987]])

>>> datadict['__header__']

'MATLAB 5.0 MAT-file, written by Octave 3.2.4, 2012-11-27
 15:45:20 UTC'

savemat(file_name, mdict, appendmat=True, format='5', 
long_field_names=False, do_compression=False,
oned_as=None)

114	100	104	89	102	91	114	114
103	105	108	130	120	132	111	128
118	119	86	72	111	103	74	112
107	103	98	96	112	112	93

受试者	1	2	3	4	5	6	7	8	9	10
右螺纹	113	105	130	101	138	118	87	116	75	96
左螺纹	137	105	133	108	115	170	103	145	78	107
受试者	11	12	13	14	15	16	17	18	19	20
右螺纹	122	103	116	107	118	103	111	104	111	89
左螺纹	84	148	147	87	166	146	123	135	112	93
受试者	21	22	23	24	25
右螺纹	78	100	89	85	88
左螺纹	76	116	78	101	123

114	100	104	89	102	91	114	114
103	105	108	130	120	132	111	128
118	119	86	72	111	103	74	112
107	103	98	96	112	112	93

龙哥盟

掠夺·扩张·投机·博弈

SciPy 数值与科学计算学习指南（全）

前言

本书涵盖的内容

您需要这本书什么

这本书面向谁

习惯用法

注意

小贴士

读者反馈

客户支持

下载示例代码

下载本书的彩色图像

错误清单

盗版

询问

第一章。SciPy 简介

SciPy 是什么？

安装 SciPy

在 Mac OS X 上安装 SciPy

在 Unix/Linux 上安装 SciPy

在 Windows 上安装 SciPy

测试 SciPy 安装

小贴士

SciPy 组织

如何查找文档

科学可视化

如何打开 IPython 笔记本

摘要

第二章。使用 NumPy 数组作为 SciPy 的第一步

对象基础

使用数据类型

索引和切片数组

数组对象

数组转换

形状选择/操作

对象计算

数组例程

创建数组的例程

组合两个或更多数组的常规操作

数组操作的相关程序

从数组中提取信息的例程

摘要

第三章。SciPy 线性代数

向量创建

向量运算

加法/减法

标量/点积

交叉/向量积 – 在三维空间向量上

创建一个矩阵

矩阵方法

矩阵之间的运算

矩阵上的函数

特征值问题和矩阵分解

通过奇异值分解进行图像压缩

求解器

摘要

第四章：SciPy 数值分析

特殊函数的评估

便利性和测试函数

单变量多项式

欧拉函数

黎曼ζ函数

艾里函数和贝里函数

贝塞尔和斯特夫函数

其他特殊函数

插值

回归

优化

最小化

根

积分

指数/对数积分

三角函数和双曲三角函数积分

椭圆积分

伽马和贝塔积分

数值积分

常微分方程

114	100	104	89	102	91	114	114
103	105	108	130	120	132	111	128
118	119	86	72	111	103	74	112
107	103	98	96	112	112	93