三维曲面能量展示与移动最小二乘法光滑化的Python实现

1 三维曲面能量展示

1.1 加载必要的软件包

为了实现我们想要的作图功能,Numpy和Matploylib是两个必须要加载的程序包。如果你的本地没有这两个程序包，我们可以通过pip来安装他们：

python -m pip install -U numpy
python -m pip install -U matplotlib

1.1.1 Numpy简介

Numpy(Numerical Python) 是 Python 语言的一个扩展程序库，支持大量的维度数组与矩阵运算，此外也针对数组运算提供大量的数学函数库。

NumPy 的前身 Numeric 最早是由 Jim Hugunin 与其它协作者共同开发，2005 年，Travis Oliphant 在 Numeric 中结合了另一个同性质的程序库 Numarray 的特色，并加入了其它扩展而开发了 NumPy。NumPy 为开放源代码并且由许多协作者共同维护开发。

NumPy 是一个运行速度非常快的数学库，主要用于数组计算，包含：

• 一个强大的N维数组对象 ndarray

• 广播功能函数

• 整合 C/C++/Fortran 代码的工具

• 线性代数、傅里叶变换、随机数生成等功能

  

1.1.2 Matplotlib简介

Matplotlib是一个Python 2D绘图库，它可以在各种平台上以各种硬拷贝格式和交互式环境生成出具有出版品质的图形。只需几行代码即可生成绘图，直方图，功率谱，条形图，错误图，散点图等。

为了简单绘图，pyplot模块提供了类似于MATLAB的界面，特别是与IPython结合使用时。 对于高级用户，您可以通过面向对象的界面或MATLAB用户熟悉的一组函数完全控制线条样式，字体属性，轴属性等

1.1.3 加载程序包

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as col
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

1.2 数据写入与数据结构修改

data = np.loadtxt("Data_PATH")

通过Numpy中的loadtxt()函数我们可以从本地读入一个文本文档并被保存为ndarray数据结构。其中Data_PATH是我们文件所在的本地路径，需要特异性替换。

X=list({}.fromkeys(data[:,0]).keys())
X_num = len(X)
Y=list({}.fromkeys(data[:,1]).keys())
Y_num = len(Y)
X, Y = np.meshgrid(Y, X)
Z = np.array(data[:,2]).reshape(X_num,Y_num,order='F')

读取数据文件后我们就能实现可视化了吗？显然不能(TwT)~我们的Axes3D需要X,Y,Z都是二维数组才能完成曲面的绘制。举一个简单的例子，我们的实验结果只有四个点(X1,Y1,Z1),(X2,Y1,Z2),(X1,Y2,Z3),(X2,Y2,Z4)，在我们的读入文件中这四个点的被顺序排列成四行三列：

于是我们需要单独取出X值，Y值，Z值，重构成右边的三个二维数组。代码上我们首先用字典格式去重复读取了第一列并重新定义为列表结构保存于X中。保存下X值一共有多少种在X_num中。同理我们得到Y和Y_num。在上例中X=(X1,X2),Y=(Y1,Y2),X_num=2,Y_num=2。以这些信息为输入重构数据为三个可用的二维数组。

1.3 数据可视化与参数调整

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
surface = ax.plot_surface(X, Y, Z)
plt.show()

通过简单的几行代码，我们就能将1.2中处理好的数据可视化为一个三维曲面，以下是分条解释每一行代码的含义：

• fig = plt.figure()

  新建画布，可以理解为整个窗口就是一个画布。并将画布保存在fig变量名下,这样我们就可以通过对fig的操作来生成我们需要的图案了。

• ax = fig.add_subplot(111, projection='3d')

  新建图形区域与类型，即在我们之前生成的窗口中添加一个子区域作为图形区域并保存在ax变量名下，其中参数111表示将整个fig窗口分为1行1列，将此图形区域添加在1号位置。也就是abc表示整个fig窗口分为a行b列，将此图形区域添加在c号位置,c不大于a·b。举一个简单的例子，223就表示整个fig窗口分为2行2列，将此图形区域添加在3号位置,也就是四等分的左下角。projection让我们可以定义我们添加的这个图形区域的类型，由于我们需要创建的是一个三维曲面图形，所以我们在这里需要设定为3d。

• surface = ax.plot_surface(X, Y, Z)

  在图形区域中画上三维曲面图。上一步中我们已经有了一个类型为3d的图形区域名字叫ax，在其上我们通过函数传递我们在1.2中处理好的X,Y,Z数据，就可以实现我们想要的数据可视化了。

• plt.show()

  显示整个画布。现在我们来把辛辛苦苦写好的代码展示出来吧~Python会以一个弹窗的形式展示我们的结果，下方的按钮似乎也并不能满足我们对于美丽图片调整的需求。这就需要我们使用更加复杂的代码命令来在展示之前修改我们的图形区域或者三维曲面。

1.3.1 三维曲面参数调整

norm = plt.Normalize(Z.min(), Z.max())
colors = cm.viridis(norm(Z))
rcount, ccount, _ = colors.shape
surface = ax.plot_surface(X, Y, Z, rcount=rcount, ccount=ccount, facecolors=colors, shade=False)
surface.set_facecolor((0,0,0,0))

1.3.2 二维投影添加

plt.contour(X, Y, Z, offset = -2, cmap = 'viridis')

1.3.3 坐标系参数调整

ax.grid(False)
ax.set_zlim(-2, 5)
ax.xaxis.set_pane_color((0,0,0,0))
ax.yaxis.set_pane_color((0,0,0,0))
ax.zaxis.set_pane_color((0,0,0,0))

1.3.4 其他

ax.view_init(20, -117)
ax.set_xlabel('C-O Distance(Å)')
ax.set_ylabel('N-O Distance(Å)')
ax.set_zlabel('Free Energy(kcal/mol)')
plt.savefig('PATH',dpi=500)

1.4 代码总结与函数化封装

import matplotlib.pyplot as plt
import matplotlib.colors as col
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
import numpy as np

def Draw3DSurface(InputDataPATH,OutputDataPATH):
	data = np.loadtxt(InputDataPATH)
	X=list({}.fromkeys(data[:,0]).keys())
	X_num = len(X)
	Y=list({}.fromkeys(data[:,1]).keys())
	Y_num = len(Y)
	X, Y = np.meshgrid(X, Y)
	Z = np.array(data[:,2]).reshape(X_num,Y_num,order='F')
	norm = plt.Normalize(Z.min(), Z.max())
	colors = cm.viridis(norm(Z))
	rcount, ccount, _ = colors.shape
	fig = plt.figure()
	ax = fig.add_subplot(111, projection='3d')
	surface = ax.plot_surface(X, Y, Z, rcount=rcount, ccount=ccount, facecolors=colors, shade=False)
	surface.set_facecolor((0,0,0,0))
	ax.grid(False)
	ax.set_zlim(-2, 5)
	ax.view_init(20, -117)
	ax.xaxis.set_pane_color((0,0,0,0))
	ax.yaxis.set_pane_color((0,0,0,0))
	ax.zaxis.set_pane_color((0,0,0,0))
	plt.contour(X, Y, Z, offset = -2, cmap = 'viridis')
	plt.show()
	plt.savefig(OutputDataPATH,dpi=500)
    
Draw3DSurface("D:/DATA/input.dat","D:/DATA/output.dat")

2 移动最小二乘法