0.问题描述
Python是一个适合算法验证的语言,其优点与缺点都十分明显,无需编译运行十分“顺滑”,但是运行效率较低。
之前利用Python编写了多晶塑性计算程序,比较复杂的算法和功能在很短的时间内可以开发完成,然而运行起来很令人头疼,一个算例需要计算近四十个小时,无法忍受需要优化程序。当然,优化程序的方法和方向有很多,比如并行计算,比如优化Python代码以多利用numpy包,比如利用numba等。在不改变算法的情况下,换一个运行效率高的需要编译的语言是个直接的思路,比如C++。基于Python编写的前后处理程序自认为很好用,不是大团队也不是专业的程序员的情况下,不想花大量时间用c++重构全部的程序,那就只改一改最耗时的计算内核。
利用C++编写程序中的迭代步计算内核,然后用Python封装,是能想到的适合的方法。解决方案为:使用C++开发计算内核的dll动态链接库,python中利用ctypes调用dll,利用Python编写总体计算程序。
下面给出一个非常简单的例子,和问题描述没什么直接关系-__-#,反正就是一个ctypes的小栗子。
1.利用Visual Studio开发C++的64位dll
实现一个特别简单的功能,输入一个数组double[9]作为一个3x3矩阵,输出一个struct,包含一个double作为该矩阵的行列式,包含一个double[9],作为该矩阵的逆矩阵。
这里利用visual studio的dll模板直接开发。参考https://docs.microsoft.com/en-us/cpp/build/walkthrough-creating-and-using-a-dynamic-link-library-cpp?view=vs-2019。
首先打开VS(这里是VS2019 社区版),创建新项目->(C++,windows)动态链接库->命名并创建。
![534464030c7f6363c923f51b7fe86211.png](https://i-blog.csdnimg.cn/blog_migrate/2bf27f44acf0b24fc014648afbad52f7.png)
向项目中添加头文件nrtools.h
和ludcmp.c
其中前者主要定义了相关的向量、矩阵和ludcmp函数等,具体可以参考Numerical Recipes in C++,而后者是生成dll的关键,其中定义了dll的输入数据结构,输出数据结构和接口函数:
#pragma once
#ifdef LUDCMPDLL_EXPORTS
#define LUDCMP_API __declspec(dllexport)
#else
#define LUDCMP_API __declspec(dllimport)
#endif
extern "C" {
//input data structure
typedef struct InputData {
int flagValue;
double oriMat[9]; //the flatten matrix to be deal with
}IPDATA;
//output data structure
typedef struct OutputData {
double determin; //determinant of the matrix
double invMat[9]; //inverse of the matrix
}OTDATA;
LUDCMP_API OTDATA interface(const IPDATA& fish);
}
其中3~7行是MSVC必须的,12-21定义了输入输出数据结构,而23行是接口函数原型。接口函数在ludcmp.cpp
文件中定义,需在项目中添加该源文件。
#include "pch.h"
#include "ludcmp.h"
#include "nrtools.h"
//function used within this cpp---------------------------------
void mat2array(const MatDoub& mat3x3, double(&ary9)[9]);
//--------------------------------------------------------------
OTDATA interface(const IPDATA& fish) {
//
MatDoub_I aa(3, 3, fish.oriMat);
//
LUdcmp aadcmp(aa);
MatDoub_O aainv(3, 3);
aadcmp.inverse(aainv); //aainv restore the inverse of aa
//
OTDATA results;
// first property
results.determin = (fish.flagValue > 0) ? aadcmp.det() : -aadcmp.det();
// second property
mat2array(aainv, results.invMat);
//
return results;
}
//function only used in this file-----------------------------------------
/* turn a define MatDoub class,3x3 to a normal double [9]
*/
void mat2array(const MatDoub& mat3x3, double(&ary9)[9]) {
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 3; j++) {
ary9[i * 3 + j] = mat3x3[i][j];
}
}
}
非常简单,主要是利用了Numerical Recipes中的LU分解算法,利用LUDCMP类来实现求逆求行列式的。这里21利用了一个小函数将3x3的矩阵转化成我们需要的数组用来输出。
代码全部码完后,利用VS的x64编译器,解决方案选“Debug"-"x64”,点击“生成”-“生成解决方案”,若无错误,将在./x64/Debug
路径下生成所需的ludcmpDll.dll
和ludcmpDll.lib
文件,这就是所需的动态链接库。因为本人的python是64位的,所以要生成64位dll,大家情况不同,要注意版本匹配。
2.利用ctypes在python中调用dll
ctypes是python的标准包(个人使用anaconda3)。可以利用ctypes调用dll,并使用dll中定义的接口函数进行计算。这里文件的结构如下,clibraty文件夹中放进了需要调用的dll文件,编写python脚本调用
./
|--test.py
|--clibray/
|--ludcmpDLL.dll
|--ludcmpDll.lib
引用ctypes
import ctypes
在使用dll中的接口函数前,需要完成三件事情1)载入dll,2)结构体声明,3)类似C++中的函数原型声明,若是没有使用结构体或联合体等复杂数据结构,可以忽略2)。
#--1)--调用dll
ludcmpCpp = ctypes.CDLL("./clibrary/ludcmpDll")
#--2)--定义结构体
#输入结构体
class IPDATA(ctypes.Structure):
_fields_ = [
("flagValue",ctypes.c_int),
("oriMat",ctypes.c_double*9)
]
#输出结构体
class OTDATA(ctypes.Structure):
_fields_ = [
("determin",ctypes.c_double),
("invMat",ctypes.c_double*9)
]
#--3)--函数声明 as in c++ head file, ludcmp就是我们要用的interface函数了
ludcmp = ludcmpCpp.interface
ludcmp.argtypes = [IPDATA]
ludcmp.restype = OTDATA
这里面可以看到ctypes中定义结构体,其实生成一个ctypes中Structure类的派生类。注意函数的输入是一个IPDATA类型。举个例子,如果想算一个3x3的numpy.narray的逆矩阵和行列式,要把该矩阵转换成IPDATA类型:
mat = np.array([[2.0,0.0,0.0],
[0.0,1.0,0.0],
[0.0,0.0,1.0]])
#
aa = IPDATA()
aa.flagValue = 20 #any intergal greater than 0
aa.oriMat = (ctypes.c_double*9)(*mat.ravel().tolist())
然后使用dll中的interface函数计算出一个OTDATA,最后将OTDATA中的属性转换为所需数据:
bb = ludcmp(aa) #bb is a OTDATA
matInv = np.array(bb.invMat).reshape(3,3) #逆矩阵,存为3x3
detmin = bb.determin #行列式
3.测试
求一个矩阵的逆矩阵和行列式,numpy里面就有很好用的函数:
numpy.linalg.inv()
numpy.linalg.det()
做一个测试函数,计算3x3矩阵的逆矩阵和行列式。使用两种方法,一种是利用上文提到的调用dll的方式,另一种是利用numpy的linalg,两种算法封装为函数:
import ctypes
import numpy as np
#define structure------------------------------
class IPDATA(ctypes.Structure):
_fields_ = [
("flagValue",ctypes.c_int),
("oriMat",ctypes.c_double*9)
]
class OTDATA(ctypes.Structure):
_fields_ = [
("determin",ctypes.c_double),
("invMat",ctypes.c_double*9)
]
def fun1(mat):
# use dll to calculate inverse and determinant
#----open dll
ludcmpCpp = ctypes.CDLL("./clibrary/ludcmpDll")
#----declaim function as in c++ head file
ludcmp = ludcmpCpp.interface
ludcmp.argtypes = [IPDATA]
ludcmp.restype = OTDATA
#----prepare for input
aa = IPDATA()
aa.flagValue = 20 #any intergal greater than 0
aa.oriMat = (ctypes.c_double*9)(*mat.ravel().tolist())
#----call interface function in dll
bb = ludcmp(aa) #bb is a OTDATA
#
return np.array(bb.invMat).reshape(3,3),bb.determin
def fun2(mat):
# use numpy linalg package
matInv = np.linalg.inv(mat)
determ = np.linalg.det(mat)
return matInv,determ
编写一个测试函数,分别调用fun1和fun2,算三个矩阵的逆矩阵和行列式,为体现性能差异,每个矩阵运行10000次。测试函数为:
def testfun(mat1,mat2,mat3):
#run n times for each mat
n = 10000
for i in range(n):
inv1,det1 = fun1(mat1)
inv2,det2 = fun1(mat2)
inv3,det3 = fun1(mat3)
#
for i in range(n):
inv1,det1 = fun2(mat1)
inv2,det2 = fun2(mat2)
inv3,det3 = fun2(mat3)
#
return None
利用cProfile测试一番:
import cProfile as cp
mat1 = np.array([[2.0,0.0,0.0], #traction along 1
[0.0,1.0,0.0],
[0.0,0.0,1.0]])
mat2 = np.array([[1.0,2.0,0.0], #shear in 12
[0.0,1.0,0.0],
[0.0,0.0,1.0]])
mat3 = np.array([[1.0,3.0,0.0], #shear in 12 and traction along 2
[0.0,0.5,0.0],
[0.0,0.0,1.0]])
#
cp.run("testfun(mat1,mat2,mat3)")
输出(截取了有用的部分)
1620004 function calls in 13.609 seconds
Ordered by: standard name
ncalls tottime percall cumtime percall filename:lineno(function)
......
......
30000 0.142 0.000 0.309 0.000 linalg.py:482(inv)
30000 0.743 0.000 12.892 0.000 test.py:25(fun1)
30000 0.033 0.000 0.668 0.000 test.py:54(fun2)
1 0.049 0.049 13.608 13.608 test.py:61(testfun)
30000 11.142 0.000 11.142 0.000 {built-in method _ctypes.LoadLibrary}
......
......
令人伤心的结果,一共运行13.6秒,fun1自己就用了12.9秒,所以能用numpy还是用numpy吧,毕竟np。
wait,wait,wait,进一步分析,怀疑真正耗时的是ctypes中的LoadLibrary,每次运行fun1都要载入一次dll,那确实浪费时间。可以进一步修改fun1,让其只负责算,不负责载入dll,载入好的函数作为参数传递过来。
def fun1plus(mat,func):
#----prepare for input
aa = IPDATA()
aa.flagValue = 20 #any intergal greater than 0
aa.oriMat = (ctypes.c_double*9)(*mat.ravel().tolist())
#----call interface function in dll
bb = func(aa) #bb is a OTDATA
#
return np.array(bb.invMat).reshape(3,3),bb.determin
测试函数改为只载入一次dll,然后再循环30000次interface函数计算:
def testfun2(mat1,mat2,mat3):
#run n times for each mat
n = 10000
# load in dll
ludcmpCpp = ctypes.CDLL("./clibrary/ludcmpDll")
#----define function as in c++ head file
ludcmp = ludcmpCpp.interface
ludcmp.argtypes = [IPDATA]
ludcmp.restype = OTDATA
for i in range(n):
inv1,det1 = fun1plus(mat1,ludcmp)
inv2,det2 = fun1plus(mat2,ludcmp)
inv3,det3 = fun1plus(mat3,ludcmp)
#
for i in range(n):
inv1,det1 = fun2(mat1)
inv2,det2 = fun2(mat2)
inv3,det3 = fun2(mat3)
#
return None
看一下这次的测试结果:cp.run("testfun2(mat1,mat2,mat3)")
1320014 function calls in 1.066 seconds
Ordered by: standard name
ncalls tottime percall cumtime percall filename:lineno(function)
......
......
30000 0.137 0.000 0.279 0.000 linalg.py:482(inv)
30000 0.375 0.000 0.435 0.000 test.py:42(fun1plus)
30000 0.029 0.000 0.591 0.000 test.py:54(fun2)
1 0.030 0.030 1.066 1.066 test.py:76(testfun2)
1 0.010 0.010 0.010 0.010 {built-in method _ctypes.LoadLibrary}
......
果然如此,载入一次dll,要花掉0.01秒,,然后用dll中的函数算30000次。对比发现,fun1plus比fun2算的时间还要少一些,那用c++和python混合编程的方式还是可以的。。。除了c++比较难码-__-
ps:
计算多晶时,一个多晶需要包含大量单晶信息,比如Taylor假设,此时一个单晶的计算核心可以使用c++编写dll来计算,这样在计算多晶时,首先载入dll,然后循环晶粒个数次单晶计算,这样比全python编程快很多,看情况吧,个人情况是使用dll运算速度快了几百倍。
查资料过程中,很多帖子说numba无脑提速,感觉以后有时间了可以好好试一试,简单上手了一下,感觉这个东西不是那么无脑,脾气很大。