Rainyboy的PYTHON小技巧备忘录(不定时更新)
本帖最后由 Rainyboy 于 2014-11-4 15:18 编辑目录
#1 搜索二维数组中最大/最小值的位置和值，1楼，2011-10
#2 绘制带三向等高线投影的3D曲面，1楼，2011-10
#3 五点法差分，2楼，2013-12
#4 2D，梯度、旋度和散度的计算和展示，3楼，2013-12
#5 设置纵坐标格式，4楼，2014-11
#6 设置默认存图分辨率，4楼，2014-11
#7 交互式标注曲线，5楼，2014-11
#1：搜索二维数组中最大/最小值的位置和值######import numpy as np
x = np.array([,
,
,
])
(M,N) = x.shape
maxvalue,minvalue = x.max(),x.min()
maxloc,minloc = x.argmax(),x.argmin()
print 'Max Location: [%d, %d]'%(maxloc/N , maxloc%N)
print 'Max Value: %.2f'%(maxvalue)
print 'Min Location: [%d, %d]'%(minloc/N , minloc%N)
print 'Min Value: %.2f'%(minvalue)
运行输出：
Max Location:
Max Value: 700.00
Min Location:
Min Value: 10.00
#2：绘制带三向等高线投影的3D曲面######
# -*- coding: cp936 -*-
import numpy as math
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm
def Surf_and_Z_Cont(X,Y,Z,
lables = ('X','Y','Z'),
Is3Dsurf = True,
IsContour = True,
IsFillContour = True,
FCTransverse = True,
IsCheck = True
):
XMax = math.max(X)
XMin = math.min(X)
Xost = XMin - (XMax - XMin)*0.1
YMax = math.max(Y)
YMin = math.min(Y)
Yost = YMax + (YMax - YMin)*0.1
ZMax = math.max(Z)
ZMin = math.min(Z)
Zost = ZMin - (ZMax - ZMin)*0.1
if IsCheck:
print 'Surf_and_Z_Cont Check information'
print 'X.shape = ', X.shape
print 'Y.shape = ', Y.shape
print 'Z.shape = ', Z.shape
print 'Surf_and_Z_Cont Check information'
if Is3Dsurf:
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X,Y,Z,rstride=1, cstride=1,alpha=0.2, lw = 0.2)
cset = ax.contour( X,Y, Z, zdir='x', offset=Xost,levels = math.linspace(XMin,XMax,20))
cset = ax.contour( X,Y, Z, zdir='y', offset=Yost,levels = math.linspace(YMin,YMax,20))
cset = ax.contour( X,Y, Z, zdir='z', offset=Zost,levels = math.linspace(ZMin,ZMax,20))
ax.set_xlabel(lables)
ax.set_xlim3d(Xost,XMax)
ax.set_ylabel(lables)
ax.set_ylim3d(YMin,Yost)
ax.set_zlabel(lables)
ax.set_zlim3d(Zost,ZMax)
if IsFillContour:
fig = plt.figure()
ax = fig.add_subplot(111)
if FCTransverse:
im = ax.imshow(Z.T, aspect = 'auto', interpolation='bilinear', origin='lower',cmap=cm.jet, extent=(XMin,XMax,YMin,YMax))
else:
im = ax.imshow(Z, aspect = 'auto', interpolation='bilinear', origin='lower',cmap=cm.jet, extent=(XMin,XMax,YMin,YMax))
CS = ax.contour(X,Y,Z,25,cmap=cm.jet)
ax.set_xlabel(lables)
ax.set_ylabel(lables)
CBI = fig.colorbar(im, orientation='vertical', shrink=0.8,drawedges=True)
if IsContour:
fig = plt.figure()
ax = fig.add_subplot(111)
CS = ax.contour(X,Y,Z,25,cmap=cm.jet)
ax.set_xlabel(lables)
ax.set_ylabel(lables)
CB = fig.colorbar(CS, shrink=0.8, extend='both')
运行结果，形如：
本帖最后由 Rainyboy 于 2013-12-10 22:17 编辑
#3 五点法数值差分，可以给定函数，也可以给定一列数
import numpy as npfrom scipy.optimize import curve_fit
def DDiff5(flist,h,N):
D = np.zeros_like(flist)
D = (flist - flist)/h
D = (flist - flist)/(2*h)
D = (flist - flist)/(2*h)
D = (flist - flist)/h
for i in range(2,N-2):
D = (-flist + 8*flist - 8*flist + flist)/(12*h)
return D
def CDiff(func,x0,x1,N,para):
xlist = np.linspace(x0,x1,N)
funclist = []
for x in xlist:
funclist.append(func(x,*para))
return xlist,DDiff5(np.array(funclist),xlist-xlist,N)
if __name__ == '__main__':
import matplotlib.pyplot as plt
def f(x,omaga,A,fi):
return x*np.sin(x)
def Df(x,omaga,A,fi):
return np.sin(x) + x*np.cos(x)
para = (1,1,np.pi/5)
x,D = CDiff(f,0,10,30,para)
plt.figure()
plt.subplot(2,1,1)
plt.plot(x,f(x,*para),label='f')
plt.grid(True)
plt.legend()
plt.subplot(2,1,2)
plt.plot(x,Df(x,*para),label='df')
plt.plot(x,D,label='diff5(f)',lw=0,marker='o')
plt.legend()
plt.grid(True)
plt.show()
由于边界点采用了中心差分和单步法，结果小有瑕疵。
本帖最后由 Rainyboy 于 2013-12-10 22:54 编辑
#4，2D，梯度、旋度和散度的计算和展示
def getGrade2D(Field,dx,dy):
Ny,Nx = Field.shape
Gx = np.zeros_like(Field)
Gy = np.zeros_like(Field)
for i in range(0,Ny):
Gx = DDiff5(Field,dx,Nx)
for j in range(0,Nx):
Gy[:,j] = DDiff5(Field[:,j],dy,Ny)
return Gx,Gy
def getDiv2D(Gx,Gy,dx,dy):
Ny,Nx = Gx.shape
Tx = np.zeros_like(Gx)
Ty = np.zeros_like(Gx)
for i in range(0,Ny):
Tx = DDiff5(Gx,dx,Nx)
for j in range(0,Nx):
Ty[:,j] = DDiff5(Gy[:,j],dy,Ny)
return Tx+Ty
def getRot2D(Gx,Gy,dx,dy):
Ny,Nx = Gx.shape
GXy = np.zeros_like(Gx)
GYx = np.zeros_like(Gx)
for i in range(0,Ny):
GYx = DDiff5(Gy,dx,Nx)
for j in range(0,Nx):
GXy[:,j] = DDiff5(Gx[:,j],dy,Ny)
return GYx-GXy
def ARRAY_CONT(X,Y,U,V,
labels=('X','Y','2D_VELOCITY_FIELD'),
copt='speed',
wopt='vary',
dst = 1):
speed = np.sqrt(U*U + V*V)
if wopt == 'constant':
wset = 0.5
elif wopt == 'vary':
wset = 1.5*speed/np.max(speed)+0.5
if copt == 'speed':
plt.streamplot(X, Y, U, V, color=speed, cmap=plt.cm.cool,density=dst, linewidth = wset)
plt.colorbar()
else:
plt.streamplot(X, Y, U, V, color=copt,density=dst, linewidth = wset)
plt.xlabel(labels)
plt.ylabel(labels)
plt.title(labels)
###例子####
###例一：绘出给定矢量场，计算其散度和旋度###
def test_5():
dx = dy = 40.0/200
Y, X = np.mgrid[-20:20:200j, -20:20:200j]
U = X*np.sin(X)-Y
V = X +Y*np.cos(Y)
plt.figure(3)
ARRAY_CONT(X,Y,U,V,dst=2)
plt.axis('equal')
Div = getDiv2D(U,V,dx,dy)
plt.figure()
plt.contourf(X,Y,Div,15)
plt.colorbar()
Rot = getRot2D(U,V,dx,dy)
plt.figure()
plt.contourf(X,Y,Rot,15)
plt.colorbar()
散度：
旋度：
###例二：计算标量场的梯度（五个点电荷构成的电场）###
def FI3(x,y):
k = 1e-5
qlist =
xlist =
ylist =
f = np.zeros_like(x)
for i in range(0,5):
f += qlist/(np.power(x-xlist,2)+np.power(y-ylist,2))
return f*k
def test_3():
Nx,Ny = (200,200)
Xmin,Xmax = -3.0,3.0
Ymin,Ymax = -3.0,3.0
dx = (Xmax-Xmin)/Nx
dy = (Ymax-Ymin)/Ny
Y, X = np.mgrid
Field = FI3(X,Y)
plt.figure()
plt.contour(X,Y,Field,15)
plt.colorbar()
Gx,Gy = getGrade2D(Field,dx,dy)
ARRAY_CONT(X,Y,Gx,Gy,labels=('X','Y',''),copt='k',wopt='constant',dst = 3)
本帖最后由 Rainyboy 于 2014-11-4 15:16 编辑
#5 设置纵坐标格式
import matplotlib.pyplot as plt
from matplotlib.ticker import FormatStrFormatter
plt.gca().yaxis.set_major_formatter(FormatStrFormatter('%3.0e'))
#6 设置默认存图分辨率
import matplotlib
matplotlib.rcParams['savefig.dpi'] = 600
本帖最后由 Rainyboy 于 2014-11-4 15:17 编辑
#7 交互式标注曲线
import numpy as np
import matplotlib.pyplot as plt
class fy:
def __init__(self):
self.ModeOn = False
self.FirstPoint = False
self.xytext = []
def onpress(self,event):
print event.key
if event.key == 'n':
print 'annotation mode on'
self.ModeOn = True
elif event.key == 'q':
print 'annotation mode off'
self.ModeOn = False
if self.ModeOn:
if event.key == '1':
self.FirstPoint = True
elif event.key == '2':
self.FirstPoint = False
def onclick(self,event):
if self.ModeOn:
if self.FirstPoint:
self.xytext = (event.xdata,event.ydata)
print '1st point saved:',self.xytext
def onpick(self,event):
if self.ModeOn:
if not self.FirstPoint:
print event.ind
ind = event.ind
thisline = event.artist
xdata, ydata = thisline.get_data()
fig = plt.figure(1)
info = str(raw_input('input your text:'))
plt.annotate(info,(xdata, ydata),self.xytext,arrowprops=dict(arrowstyle="->"))
fig.canvas.draw()
f = fy()
fig = plt.figure(1)
fig.canvas.mpl_connect('key_press_event', f.onpress)
fig.canvas.mpl_connect('pick_event', f.onpick)
fig.canvas.mpl_connect('button_press_event', f.onclick)
x = np.linspace(0,100,100)
plt.plot(x,5*x,x,12*x,picker = 5)
plt.show()
顶啊！想不到兄弟还在坚持 说太多都是矫情，院长威武啊，必须支持 smtmobly 发表于 2015-1-19 08:40
顶啊！想不到兄弟还在坚持
其实也很少更新了。。。有些时候上来看看 zhangnan3509 发表于 2015-1-19 08:49
说太多都是矫情，院长威武啊，必须支持
哈哈，谢谢，自娱自乐而已
{:{05}:} 以后我也自娱自乐一下，现在咱们的版块有没有新的，如果有控制类的，我每天自娱自乐去 zhangnan3509 发表于 2015-1-19 09:13
以后我也自娱自乐一下，现在咱们的版块有没有新的，如果有控制类的，我每天自娱自乐去
我想前辈大可以在算法区玩玩，也可以去“控制、可靠性及优化”分区，那里也是我负责打理
{:{05}:}
算法区的好处是，只有在算法区，贴出来的程序语言可以显示成语法高亮，其他分区没有这个特色。
所以如果你的代码量比较大还是建议在这里，如果是公式，图表类的就无所谓了。 控制类的都是什么内容？ zhangnan3509 发表于 2015-1-19 09:38
控制类的都是什么内容？
我觉得无所谓了，想发什么就发什么吧 第七为什么 点击后显示不了
页:
[1]