我的目标是在Python中编写一个用于频谱有限元的小型库，为此我尝试使用Boost将C++库扩展到Python中，希望能使我的代码更快。
class Quad {
public:
Quad(int, int);
double integrate(boost::function<double(std::vector<double> const&)> const&);
double integrate_wrapper(boost::python::object const&);
std::vector< std::vector<double> > nodes;
std::vector<double> weights;
};
...
namespace std {
typedef std::vector< std::vector< std::vector<double> > > cube;
typedef std::vector< std::vector<double> > mat;
typedef std::vector<double> vec;
}
...
double Quad::integrate(boost::function<double(vec const&)> const& func) {
double result = 0.;
for (unsigned int i = 0; i < nodes.size(); ++i) {
result += func(nodes[i]) * weights[i];
}
return result;
}
// ---- PYTHON WRAPPER ----
double Quad::integrate_wrapper(boost::python::object const& func) {
std::function<double(vec const&)> lambda;
switch (this->nodes[0].size()) {
case 1: lambda = [&func](vec const& v) -> double { return boost::python::extract<double>(func (v[0])); }; break;
case 2: lambda = [&func](vec const& v) -> double { return boost::python::extract<double>(func(v[0], v[1])); }; break;
case 3: lambda = [&func](vec const& v) -> double { return boost::python::extract<double>(func(v[0], v[1], v[2])); }; break;
default: cout << "Dimension must be 1, 2, or 3" << endl; exit(0);
}
return integrate(lambda);
}
// ---- EXPOSE TO PYTHON ----
BOOST_PYTHON_MODULE(hermite)
{
using namespace boost::python;
class_<std::vec>("double_vector")
.def(vector_indexing_suite<std::vec>())
;
class_<std::mat>("double_mat")
.def(vector_indexing_suite<std::mat>())
;
class_<Quad>("Quad", init<int,int>())
.def("integrate", &Quad::integrate_wrapper)
.def_readonly("nodes", &Quad::nodes)
.def_readonly("weights", &Quad::weights)
;
}
我比较了三种不同的方法来计算两个函数的积分。这两个函数是:
- 函数
f1(x,y,z) = x*x - 更难评估的函数：
f2(x,y,z) = np.cos(2*x+2*y+2*z) + x*y + np.exp(-z*z) +np.cos(2*x+2*y+2*z) + x*y + np.exp(-z*z) +np.cos(2*x+2*y+2*z) + x*y + np.exp(-z*z) +np.cos(2*x+2*y+2*z) + x*y + np.exp(-z*z) +np.cos(2*x+2*y+2*z) + x*y + np.exp(-z*z)
Call the library from a C++ program:
double func(vector<double> v) { return F1_OR_F2; } int main() { hermite::Quad quadrature(100, 3); double result = quadrature.integrate(func); cout << "Result = " << result << endl; }Call the library from a Python script:
import hermite def function(x, y, z): return F1_OR_F2 my_quad = hermite.Quad(100, 3) result = my_quad.integrate(function)Use a
forloop in Python:import hermite def function(x, y, z): return F1_OR_F2 my_quad = hermite.Quad(100, 3) weights = my_quad.weights nodes = my_quad.nodes result = 0. for i in range(len(weights)): result += weights[i] * function(nodes[i][0], nodes[i][1], nodes[i][2])
time命令测量方法1的时间，使用python模块time测量方法2和方法3的时间，C++代码使用Cmake编译，并且使用set (CMAKE_BUILD_TYPE Release)）：对于
f1：
- 方法1：0.07s用户 0.01s系统 99%CPU 0.083总计
- 方法2：0.19秒
- 方法3：3.06秒对于
f2：
- 方法1：0.28s用户 0.01s系统 99%CPU 0.289总计
- 方法2：12.47秒
- 方法3：16.31秒根据这些结果，我的问题如下： 1. 为什么第一种方法比第二种方法快那么多？ 2. 可以改进Python包装器以达到方法1和方法2之间可比较的性能吗？ 3. 为什么方法2对于函数积分难度更加敏感而不是方法3？
编辑：我还尝试定义一个接受字符串作为参数的函数，将其写入文件，然后编译该文件并动态加载生成的
.so文件。double Quad::integrate_from_string(string const& function_body) {
// Write function to file
ofstream helper_file;
helper_file.open("/tmp/helper_function.cpp");
helper_file << "#include <vector>\n#include <cmath>\n";
helper_file << "extern \"C\" double toIntegrate(std::vector<double> v) {\n";
helper_file << " return " << function_body << ";\n}";
helper_file.close();
// Compile file
system("c++ /tmp/helper_function.cpp -o /tmp/helper_function.so -shared -fPIC");
// Load function dynamically
typedef double (*vec_func)(vec);
void *function_so = dlopen("/tmp/helper_function.so", RTLD_NOW);
vec_func func = (vec_func) dlsym(function_so, "toIntegrate");
double result = integrate(func);
dlclose(function_so);
return result;
}
这段代码比较混乱，可能不太易于移植，但是它功能强大，并且与sympy的ccode函数兼容。
第二次编辑：我已经用Numpy纯Python重写了这个函数。
import numpy as np
import numpy.polynomial.hermite_e as herm
import time
def integrate(function, degrees):
dim = len(degrees)
nodes_multidim = []
weights_multidim = []
for i in range(dim):
nodes_1d, weights_1d = herm.hermegauss(degrees[i])
nodes_multidim.append(nodes_1d)
weights_multidim.append(weights_1d)
grid_nodes = np.meshgrid(*nodes_multidim)
grid_weights = np.meshgrid(*weights_multidim)
nodes_flattened = []
weights_flattened = []
for i in range(dim):
nodes_flattened.append(grid_nodes[i].flatten())
weights_flattened.append(grid_weights[i].flatten())
nodes = np.vstack(nodes_flattened)
weights = np.prod(np.vstack(weights_flattened), axis=0)
return np.dot(function(nodes), weights)
def function(v): return F1_OR_F2
result = integrate(function, [100,100,100])
print("-> Result = " + str(result) + ", Time = " + str(end-start))
有点出乎意料（至少对我来说），这种方法与纯C++实现的性能没有显著差异。特别地，f1需要0.059秒，f2需要0.36秒。
Quad::integrate中的计算过于简单，而且func被频繁调用。Python 版本大部分时间都花费在从 Python 解释器切换到 C++ 代码上。 - llllllllllset (CMAKE_BUILD_TYPE Release)重新运行了代码，并更新了时间。 - Rastapopoulos