我的目标是在 Python 中为光谱有限元编写一个小型库，为此我尝试使用 Boost 将 Python 扩展为 C++ 库，希望它能让我的代码更快。

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)

使用的方法是:

1. 从 C++ 程序调用库:

   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;
   }
   

2. 从 Python 脚本调用库:

   import hermite
   def function(x, y, z): return F1_OR_F2
   my_quad = hermite.Quad(100, 3)
   result = my_quad.integrate(function)
   

3. 在 Python 中使用 for 循环:

   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])
   

这是每个方法的执行时间(方法 1 使用 time 命令测量时间，方法 2 和 3 使用 python 模块 time , C++ 代码是使用 Cmake 和 set (CMAKE_BUILD_TYPE Release))

• 对于 f1:

  • 方法 1:0.07s 用户 0.01s 系统 99% cpu 0.083 总计
  • 方法二:0.19s
  • 方法三:3.06s

• 对于 f2:

  • 方法 1:0.28s 用户 0.01s 系统 99% cpu 0.289 总计
  • 方法二:12.47s
  • 方法三:16.31s

根据这些结果，我的问题如下:

• 为什么第一种方法比第二种方法快这么多？

• 是否可以改进 python 包装器以达到方法 1 和方法 2 之间相当的性能？

• 为什么方法 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.059s，f2 需要 0.36s。

最佳答案

您的函数按值获取 vector ，这涉及复制 vector 。 integrate_wrapper 做额外的复制。

在这些 lambda 中通过引用接受 boost::function 并通过引用捕获 func 也是有意义的。

将这些更改为(注意 & 和 const& 位):

double integrate(boost::function<double(std::vector<double> const&)> const&);

double Quad::integrate_wrapper(boost::python::object 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);
}

尽管如此，从 C++ 调用 Python 函数比调用 C++ 函数更昂贵。

人们通常使用 numpy 在 Python 中进行快速线性代数，它使用 SIMD 进行许多常见操作。在推出 C++ 实现之前，您应该首先考虑使用 numpy。在 C++ 中，您必须在 Eigen 上使用英特尔 MKL 进行矢量化。

关于python - 为什么 C++ 比带有 boost 的 python 快得多？，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/48749668/