Rainyboy的PYTHON小技巧备忘录(不定时更新)

Rainyboy

目录
#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：搜索二维数组中最大/最小值的位置和值######[code=Python width=600px]import numpy as np

x = np.array([[10,20,30,40],
              [50,60,700,80],
              [15,25,95,45],
              [55,65,75,85]
    ])
(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)[/code]

运行输出：
Max Location: [1, 2]
Max Value: 700.00
Min Location: [0, 0]
Min Value: 10.00


#2：绘制带三向等高线投影的3D曲面######
[code=Python width=600px]# -*- 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[0])
        ax.set_xlim3d(Xost,XMax)
        ax.set_ylabel(lables[1])
        ax.set_ylim3d(YMin,Yost)
        ax.set_zlabel(lables[2])
        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[0])
        ax.set_ylabel(lables[1])
        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[0])
        ax.set_ylabel(lables[1])
        CB = fig.colorbar(CS, shrink=0.8, extend='both')[/code]
运行结果，形如：

Rainyboy

#3 五点法数值差分，可以给定函数，也可以给定一列数
[code=Python width=600px]import numpy as npfrom scipy.optimize import curve_fit

def DDiff5(flist,h,N):
    D = np.zeros_like(flist)
    D[0] = (flist[1] - flist[0])/h
    D[1] = (flist[2] - flist[0])/(2*h)
    D[N-2] = (flist[N-1] - flist[N-3])/(2*h)
    D[N-1] = (flist[N-1] - flist[N-2])/h
    for i in range(2,N-2):
        D = (-flist[i+2] + 8*flist[i+1] - 8*flist[i-1] + flist[i-2])/(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[1]-xlist[0],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()[/code]



由于边界点采用了中心差分和单步法，结果小有瑕疵。

Rainyboy

#4，2D，梯度、旋度和散度的计算和展示

[code=Python width=600px]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[i,:] = DDiff5(Field[i,:],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[i,:] = DDiff5(Gx[i,:],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[i,:] = DDiff5(Gy[i,:],dx,Nx)
    for j in range(0,Nx):
        GXy[:,j] = DDiff5(Gx[:,j],dy,Ny)
    return GYx-GXy[/code]


[code=Python width=600px]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[0])
    plt.ylabel(labels[1])
    plt.title(labels[2])[/code]

###例子####

###例一：绘出给定矢量场，计算其散度和旋度###
[code=Python width=600px]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()
[/code]


散度：



旋度：



###例二：计算标量场的梯度（五个点电荷构成的电场）###
[code=Python width=600px]def FI3(x,y):
    k = 1e-5
    qlist = [3.5,-1,+1,-1,-1]
    xlist = [0,0,0,-1,1]
    ylist = [0,1,-1,0,0]
    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[Ymin:Ymax:Ny*1j,Xmin:Xmax:Nx*1j]
    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)[/code]

Rainyboy

#5 设置纵坐标格式
[code=Python width=600px]import matplotlib.pyplot as plt
from matplotlib.ticker import FormatStrFormatter
plt.gca().yaxis.set_major_formatter(FormatStrFormatter('%3.0e'))[/code]




#6 设置默认存图分辨率

[code=Python width=600px]import matplotlib
matplotlib.rcParams['savefig.dpi'] = 600
[/code]

Rainyboy

#7 交互式标注曲线

[code=Python width=600px]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[0]
                thisline = event.artist
                xdata, ydata = thisline.get_data()
                fig = plt.figure(1)
                info = str(raw_input('input your text:'))
                plt.annotate(info,(xdata[ind], ydata[ind]),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()[/code]

smtmobly

顶啊！想不到兄弟还在坚持

zhangnan3509

说太多都是矫情，院长威武啊，必须支持

Rainyboy

smtmobly 发表于 2015-1-19 08:40
顶啊！想不到兄弟还在坚持

其实也很少更新了。。。有些时候上来看看

Rainyboy

zhangnan3509 发表于 2015-1-19 08:49
说太多都是矫情，院长威武啊，必须支持

哈哈，谢谢，自娱自乐而已

zhangnan3509

以后我也自娱自乐一下，现在咱们的版块有没有新的，如果有控制类的，我每天自娱自乐去

Rainyboy

zhangnan3509 发表于 2015-1-19 09:13
以后我也自娱自乐一下，现在咱们的版块有没有新的，如果有控制类的，我每天自娱自乐去

我想前辈大可以在算法区玩玩，也可以去“控制、可靠性及优化”分区，那里也是我负责打理

算法区的好处是，只有在算法区，贴出来的程序语言可以显示成语法高亮，其他分区没有这个特色。
所以如果你的代码量比较大还是建议在这里，如果是公式，图表类的就无所谓了。

zhangnan3509

控制类的都是什么内容？

Rainyboy

zhangnan3509 发表于 2015-1-19 09:38
控制类的都是什么内容？

我觉得无所谓了，想发什么就发什么吧

xxls82

第七为什么 点击后显示不了