Pythonバインディングの実装方法（pybind11を使用）

pybind11でのMPIの使用

https://github.com/pybind/pybind11/issues/23

parametersは、dictと対応させる。

py::overload_castのエラー https://github.com/pybind/pybind11/issues/1153

戻り値がテンプレートの場合

Eigen3を継承した以下のクラスで、スーパークラスであるEigen3のメンバをそのまま使えるようにする。

localized_vector
localized_matrix

boost::anyは、dictに変換できない。

      .def("get_map", &rokko::parameters::get_map);

以下を分けるためのサブモジュールは不要？

dense
- serial
- parallel
sparse

テストプログラム test_parameters.py

import rokko

params = rokko.parameters()
params.set("a", 22)

keys = params.keys()
print(keys)
#map = params.get_map()

parametersにおいて、Pythonの辞書のような代入を可能にする。 https://qiita.com/icoxfog417/items/83a77a648b71a38a1bd5

ビルドが通った：

      .def("diagonalize", static_cast<rokko::parameters (rokko::serial_dense_ev::*)(rokko::localized_matrix<double, rokko::matrix_col_major>&, rokko::localized_vector<double>&)>(&rokko::serial_dense_ev::diagonalize<rokko::matrix_col_major, rokko::localized_vector<double>>));

ビルドが通った。

      .def("diagonalize", py::overload_cast<rokko::localized_matrix<double, rokko::matrix_col_major>&, rokko::localized_vector<double>&>(&rokko::serial_dense_ev::diagonalize<rokko::matrix_col_major, rokko::localized_vector<double>>));

boolの方を先に書かないと、bool型が整数に変換されてしまう。

      .def("set", &parameters::set<bool>) // needed to place before <int>
      .def("set", &parameters::set<int>)

py::arg("value").noconvert()を付けても整数型に変換されてしまった。

多重定義を楽にするため、入力のparamsを、Python側で省略可能にするか？　デフォルト引数を設定するか？

, py::arg("params") = parameters()

参考：https://pybind11.readthedocs.io/en/stable/basics.html#default-args

テンプレートVECの特殊化は、Eigen3と対応するlocalized_vectorだけを行えば良さそう。

テンプレートクラス https://stackoverflow.com/questions/47487888/pybind11-template-class-of-many-types

https://stackoverflow.com/questions/44925081/how-to-wrap-templated-classes-with-pybind11

boost::variantを用いる。 https://rhysd.hatenablog.com/entry/2014/01/03/024421 https://qiita.com/Ushio/items/aadc8de96af0f1d04f1a

boost::variantでラムダ式を用いる。 https://faithandbrave.hateblo.jp/entry/20131201/1385824957

majorはキーワード引数とするか？

matrix_major=row https://llvm.org/docs/HowToSetUpLLVMStyleRTTI.html

以下で紹介されているpolymorphic_type_hook https://github.com/pybind/pybind11/blob/master/docs/advanced/classes.rst ダウンロードしたファイルでは、pybind11/docs/advanced/classes.rst 例題は、pybind11/tests/test_tagbased_polymorphic.cppにある。

https://github.com/pybind/pybind11/issues/387

    pybind11_add_module(pyrokko pyrokko.cpp)
    target_link_libraries(pyrokko PRIVATE rokko)

map = pyrokko.mapping_bc(pyrokko.matrix_major.row)

mpi4py https://github.com/pybind/pybind11/issues/23

静的であるテンプレートを動的なPythonでどのようにラップするかが難しい。

pybind11のpy::EigenDRef<Eigen::MatrixXd>のソースコードを読んでみる。

https://pybind11.readthedocs.io/en/stable/advanced/cast/stl.html visit_helperクラスの特殊化

pybind11/tests/test_stl.cpp

    m.def("load_variant", [](variant<int, std::string, double, std::nullptr_t> v) {
        return py::detail::visit_helper<variant>::call(visitor(), v);
    });
    m.def("load_variant_2pass", [](variant<double, int> v) {
        return py::detail::visit_helper<variant>::call(visitor(), v);
    });
    m.def("cast_variant", []() {
        using V = variant<int, std::string>;
        return py::make_tuple(V(5), V("Hello"));
    });

クラスpy::object:local_row_majorが引数の時、rokko::localized_matrix<matrix_row_major>を返す。

不採用：

    py::class_<matrix_row_major>(m,"matrix_row_major")
      .def(py::init<>());
    m.attr("matrix_row_major") = py::cast(matrix_row_major{});

    py::class_<matrix_col_major>(m,"matrix_col_major")
      .def(py::init<>());
    m.attr("matrix_col_major") = py::cast(matrix_col_major{});

採用：

    py::enum_<matrix_major_enum>(m, "matrix_major")
      .value("row", matrix_major_enum::row)
      .value("col", matrix_major_enum::col);
    //      .export_values();

同じレベルなので、export_valuesは不要

https://stackoverflow.com/questions/52657173/sharing-an-mpi-communicator-using-pybind11

python3.6にmpi4pyをインストールした。

sudo pip install mpi4py

MPIコミュニケータをpybind11で使う。 https://www.oipapio.com/question-5669526 https://www.oipapio.com/question-376007

import pyrokko
import mpi4py.MPI

comm = pyrokko.comm(mpi4py.MPI.COMM_WORLD)

import pyrokko
import mpi4py.MPI

dim = 10
g = pyrokko.grid(mpi4py.MPI.COMM_WORLD, pyrokko.grid_col_major, 0)
#g = pyrokko.grid()

print(g.major)
solver = pyrokko.parallel_dense_ev()#"elpa")
map = solver.default_mapping(dim, g)
#map = pyrokko.mapping_bc(dim, 1, g, pyrokko.matrix_major.col)
mat = pyrokko.distributed_matrix(map)
eigvecs = pyrokko.distributed_matrix(map)
pyrokko.minij_matrix.generate(mat)

#mat.print()
#print(mat)
print(mat.major)

params = pyrokko.parameters()
eigvals = pyrokko.localized_vector(dim);

solver.initialize()
#solver.diagonalize(mat, eigvals, params)
solver.diagonalize(mat, eigvals, eigvecs, params)
#print(eigvals)

if (mpi4py.MPI.COMM_WORLD.Get_rank() == 0):
    print("Computed Eigenvalues =\n")
    eigvals.print()

eigvecs.print()

solver.finalize()

gridのmajorをユーザが確認できるように、Pythonのプロパティとして持つ。 pybind11のdef_property_readonlyを使う。 std::stringはクラスメンバとしては持たずに、アクセスされる度に生成する。

MPIコミュニケータ自体を扱うクラスを新たに作る必要はないため、以下は不要。

    py::class_<wrap_communicator>(m, "comm")
       .def(py::init( []() { return new wrap_communicator(); } ))
       .def(py::init( [](pybind11::handle const& comm) {
             const int rc = import_mpi4py();
             assert(rc==0);
             return new wrap_communicator(*PyMPIComm_Get(comm.ptr()));
           } ))
       .def("get_comm", &call_get_comm);

initializeのargcとargv

localized_matrixは、テンプレートパラメータを意識しなくても良くなるように、wrapしないとダメ。

gridのプロパティとして、myrowとmycolのタプルにした。

distributed_matrixがメンバとして持つmapは、shared_ptrとするべきか？

pybind11の情報 https://www.hiromasa.info/posts/14/

引数のいろいろ：

g = pyrokko.grid(mpi4py.MPI.COMM_WORLD, pyrokko.grid_col_major, 0)
map = pyrokko.mapping_bc(dim, 1, g, pyrokko.matrix_major.col)
print(mat.major)

wrap_parallel_dense_ev → py_parallel_dense_ev

クラスparallel_sparse_evの関数insertの引数としてstd::vectorを使用する。これをPythonのリストに対応させるため、namespaceの外に書く以下は、不要であった。

PYBIND11_MAKE_OPAQUE(std::vector<int>);
PYBIND11_MAKE_OPAQUE(std::vector<double>);

std::function

クラスのメンバ関数を代入 https://qiita.com/grainrigi/items/1aeeaf19d75d9827d037

  using namespace std::placeholders;
  std::function<void(std::vector<double>const&,std::vector<double>&)>
    f = std::bind(&laplacian_op::multiply_d, &mat_orig, _1, _2);

以下もコンストラクタで受け取れるようにしたいが、メンバとして両方を保持する必要があるし、関数の引数の型は自動で変換できない。

void multiply(const double* x, double* y)

以下のように、double[]としてみた。

  py::class_<distributed_mfree_function>(m, "distributed_mfree")
    .def(py::init<std::function<void(const double[],double[])>, int, int>());

すると、Pythonプログラムの実行時に以下のエラーとなる：

Traceback (most recent call last):
  File "laplacian_mfree_mpi.py", line 92, in <module>
    params_out = solver.diagonalize(mat, params)
  File "laplacian_mfree_mpi.py", line 44, in multiply
    self.comm.send(x[0], dest=self.myrank-1)
TypeError: 'float' object is not subscriptable

std::vectorとすると、変更したものをC++側にどのように反映させる方法がわからない。

pybind11-2.3.0/tests/test_stl.pyの以下を参考

    // test_stl_pass_by_pointer
    m.def("stl_pass_by_pointer", [](std::vector<int>* v) { return *v; }, "v"_a=nullptr);

mutable https://qiita.com/lucidfrontier45/items/183526df954c1d6580ba

https://stackoverflow.com/questions/44659924/returning-numpy-arrays-via-pybind11

std::function<void(py::array_t<double> const&,py::array_t<double>)> multiply_;に対するエラー：

Users/sakashitatatsuya/development/rokko/rokko/pyrokko_distributed_mfree.hpp:130:5: error: no matching function for call to object of type 'const std::function<void
      (const py::array_t<double> &, py::array_t<double>)>'
    multiply_(vec_x, &vec_y);
    ^~~~~~~~~
/opt/local/libexec/llvm-6.0/include/c++/v1/functional:1677:9: note: candidate function not viable: no known conversion from 'const std::vector<double>' to
      'const pybind11::array_t<double, 16>' for 1st argument
    _Rp operator()(_ArgTypes...) const;
        ^

https://pybind11.readthedocs.io/en/stable/advanced/pycpp/numpy.html

By exchanging py::buffer with py::array in the above snippet, we can restrict the function so that it only accepts NumPy arrays (rather than any type of Python object satisfying the buffer protocol).

    std::cout << "before ptr_y=" << y << " ptr_vec_y=" << vec_y.request().ptr << std::endl;
    multiply_(vec_x, vec_y);
    std::cout << "after ptr_y=" << y << " ptr_vec_y=" << vec_y.request().ptr << std::endl;

https://stackoverflow.com/questions/42645228/cast-numpy-array-to-from-custom-c-matrix-class-using-pybind11

https://github.com/tdegeus/pybind11_examples/tree/master/09_numpy_cpp-custom-matrix

https://github.com/pybind/pybind11/issues/1231

Traceback (most recent call last):
  File "laplacian_mfree_mpi.py", line 101, in <module>
    params_out = solver.diagonalize(mat, params)
RuntimeError: make_tuple(): unable to convert argument of type 'Eigen::Ref<Eigen::Matrix<double, 1, -1, 1, 1, -1>, 0, Eigen::InnerStride<1> >' to Python object

    Eigen::VectorXd X = Eigen::Map<MyVec>((double*)x, num_local_rows_);
    Eigen::VectorXd Y = Eigen::Map<MyVec>(y, num_local_rows_);

    MyVec X = Eigen::Map<MyVec>((double*)x, num_local_rows_);
    MyVec Y = Eigen::Map<MyVec>(y, num_local_rows_);

using MyVec = Eigen::Matrix<double, 1, Eigen::Dynamic, Eigen::RowMajor>; using MyRef = Eigen::Ref; MyRef&はMyRefと同等？ &は付けても害がない？ https://stackoverflow.com/questions/21132538/correct-usage-of-the-eigenref-class

ベクトルの型

Pythonでは倍精度がデフォルトで単精度は使わない。よって、Python型でサポートする型（C++ではテンプレート型）はdoubleだけで良い？

Eigen3

localized_matrixのベースクラスがEigen::MatrixXdであることをpybind11に知らせる：

    py::class_<Eigen::MatrixXd>(m, "MatrixXd");
    
    py::class_<localized_matrix<double>,Eigen::MatrixXd>(m, "localized_matrix")
      .def(py::init<>())
      .def("rows", &localized_matrix<double>::rows)
      .def("print", &localized_matrix<double>::print);

上記は、意味がなさそう。

localized_vector

void print(localized_vector<double> const& vec) {
  vec.print();
}

m.def("print", &print);

以上を定義すると、関数printの定義が上書きされて、以下がエラーとなってしまう：

print("abc")

pybind11/eigen3.h

以下のインクルードを行うと、

#include <pybind11/eigen.h>

localized_vector.hppの以下がエラー

template<typename T, int ROWS = Eigen::Dynamic>
class localized_vector : public Eigen::Matrix<T, ROWS, 1, Eigen::RowMajor> {

以下はインクルードガードがない？２回インクルードするとエラーとなる。

#include <pybind11/eigen.h>

numpyとの連携

tests/test_buffers.cppを参考にする。

pyrokkoのライブラリ名をrokkoに変更する：

    set_target_properties(pyrokko PROPERTIES OUTPUT_NAME rokko)

すると、importでエラーとなる：

Python 3.6.9 (default, Jul  5 2019, 22:57:13) 
[GCC 4.2.1 Compatible Apple LLVM 9.1.0 (clang-902.0.39.2)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import rokko
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ImportError: dynamic module does not define module export function (PyInit_rokko)

librokko.dylibの方が優先順位が高いようで、これをimportしようとしてしまうためである。そのため、ライブラリ名はpyrokkoのままとする。

以下のSOLVERはテンプレートであるため、wrap_parallel_sparse_solverを受け取ることができる。

  template<typename SOLVER>
  explicit distributed_crs_matrix(int row_dim, int col_dim, SOLVER& solver_in) {

そのため、以下のラッパーは不要であった。

class wrap_distributed_crs_matrix : public distributed_crs_matrix {
public:

  template<typename SOLVER>
  explicit wrap_distributed_crs_matrix(int row_dim, int col_dim, SOLVER& solver_in) {
    mat = solver_in.create_distributed_crs_matrix(row_dim, col_dim);
    solver_name_ = solver_in.get_solver_name();
  }
  template<typename SOLVER>
  explicit wrap_distributed_crs_matrix(int row_dim, int col_dim, int num_entries_per_row, SOLVER& solver_in) {
    mat = solver_in.create_distributed_crs_matrix(row_dim, col_dim, num_entries_per_row);
    solver_name_ = solver_in.get_solver_name();
  }
};

std::vectorがあるから、distributed_vectorへのラッパーは不要か？

残っていること

__repr__を以下を定義すること。

    .def("__repr__", [](wrap_distributed_matrix& mat) { mat.print(); })

https://github.com/pybind/pybind11/issues/1201

mpi4py

sendrecv https://stackoverflow.com/questions/32278498/usage-of-mpi4py-sendrecv

boost::any

以下は活用できなさそう。

  py::bind_map<std::map<std::string, boost::any>>(m, "parameters_map");

pybin11では、boost::anyに対する変換機能がない。つまり、boost::anyに代入したC++オブジェクトをPythonオブジェクトに変換できない。（そのような事がそもそも原理的に可能かも不明。）

from pyrokko import *
import numpy
a = numpy.ndarray((2,2), dtype='float', order='F')
mat = localized_matrix(2,2, matrix_major.row)
mat.set_ndarray(a)

例外が発生するようにした。

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: Cannot set col-major ndarray to row-major localized_matrix

>>> d = mat.ndarray
>>> d
array([[-1.49166815e-154, -1.73059424e-077],
       [-1.49166815e-154, -1.49166815e-154]])
>>> d[0,0] = 9
>>> d
array([[ 9.00000000e+000, -1.73059424e-077],
       [-1.49166815e-154, -1.49166815e-154]])
>>> mat.ndarray
array([[ 9.00000000e+000, -1.73059424e-077],
       [-1.49166815e-154, -1.49166815e-154]])

distributed_matrixの読み書き

  template <typename MATRIX_MAJOR>
  auto get_eigen_map() {
    return Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, detail::eigen3_matrix_major<MATRIX_MAJOR>::value>,0,Eigen::OuterStride<>>(get_array_pointer<MATRIX_MAJOR>(), get_m_local(), get_n_local(), Eigen::OuterStride<Eigen::Dynamic>(get_lld()));
  }

    .def("get_ndarray_col", &wrap_distributed_matrix::get_eigen_map<matrix_col_major>, py::return_value_policy::reference_internal)
    .def("get_ndarray_row", &wrap_distributed_matrix::get_eigen_map<matrix_row_major>, py::return_value_policy::reference_internal)

コピーが作成されてしまい、書き込んだ結果がdistributed_matrixに反映されない。 C++では戻り値の型の違いでオーバーロードはできない。そのため、get_ndarrayを共通化できない。

https://codeday.me/jp/qa/20190710/1223520.html

  py::object get_object() {
    Eigen::MatrixXd mat;
    if (is_col)
      new (& mat) Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, detail::eigen3_matrix_major<matrix_col_major>::value>,0,Eigen::OuterStride<>>( get_array_pointer<matrix_col_major>(), get_m_local(), get_n_local(), Eigen::OuterStride<Eigen::Dynamic>(get_lld()) );
    else
      new (& mat) Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, detail::eigen3_matrix_major<matrix_row_major>::value>,0,Eigen::OuterStride<>>( get_array_pointer<matrix_row_major>(), get_m_local(), get_n_local(), Eigen::OuterStride<Eigen::Dynamic>(get_lld()) );      
    return py::cast(&mat);
  }

Eigen::Matrixには、ストライドが指定できないので、値がおかしくなる。

Rank = 0, myrow = 0, mycol = 0
  1  1
  1  2
col
lld=2
MAP_c:
1 1
1 2
<class 'numpy.ndarray'>
  C_CONTIGUOUS : False
  F_CONTIGUOUS : True
  OWNDATA : False
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY : False
d= [[ 1.49166815e-154 -1.29074124e-231]
 [ 1.00000000e+000  2.78134232e-309]]
d2= [[ 1.49166815e-154  1.00000000e+000]
 [ 1.29074124e-231  2.78134232e-309]]
Rank = 0, myrow = 0, mycol = 0
  1.49167e-154  1
  1.29074e-231  2.78134e-309

py::arrayを使う。

  py::array*const make_array() {
    if (is_col)
      arr = py::array_t<double>({get_m_local(), get_n_local()}, {8, 8*get_lld()}, get_array_pointer<matrix_col_major>());
    else
      arr = py::array_t<double>({get_m_local(), get_n_local()}, {8*get_lld(), 8}, get_array_pointer<matrix_row_major>());
    return &arr;
  }

  template <typename MATRIX_MAJOR>
  auto const get_eigen_map_ptr() {
    auto map = new Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, detail::eigen3_matrix_major<MATRIX_MAJOR>::value>,0,Eigen::OuterStride<>>(get_array_pointer<MATRIX_MAJOR>(), get_m_local(), get_n_local(), Eigen::OuterStride<Eigen::Dynamic>(get_lld()));
    return map;
  }

https://stackoverflow.com/questions/49179356/pybind11-create-and-return-numpy-array-from-c-side

py::array_tだけで完結する。 matrix majorはストライドから判断されるので、テンプレートで与える必要はない。そのため、オーバーロードが不要となる。

  py::array_t<double>* get_ndarray() {
    py::array_t<double>* aa;
    if (is_col)
      array = new py::array_t<double>({get_m_local(), get_n_local()}, {8, 8*get_lld()}, get_array_pointer<matrix_col_major>());
    else
      array = new py::array_t<double>({get_m_local(), get_n_local()}, {8*get_lld(), 8}, get_array_pointer<matrix_row_major>());
    return array;
  }

.def("get_ndarray", &wrap_distributed_matrix::get_ndarray, py::return_value_policy::take_ownership)

https://stackoverflow.com/questions/49181258/pybind11-create-numpy-view-of-data

  void set_ndarray(py::array mat) {
    py::array_t<double> array;
    if (is_col)
      array = py::array_t<double>({get_m_local(), get_n_local()}, {sizeof(double), sizeof(double)*get_lld()}, get_array_pointer<matrix_col_major>());
    else
      array = py::array_t<double>({get_m_local(), get_n_local()}, {sizeof(double)*get_lld(), sizeof(double)}, get_array_pointer<matrix_row_major>());

    *array.mutable_data(0,0) = 8.;
    array = mat;
  }

  void set_ndarray(py::array_t<double> const& mat) {
    py::array_t<double> array;
    if (is_col) {
      array = py::array_t<double>({get_m_local(), get_n_local()}, {sizeof(double), sizeof(double)*get_lld()}, get_array_pointer<matrix_col_major>(), py::cast(*this));
    } else {
      array = py::array_t<double>({get_m_local(), get_n_local()}, {sizeof(double)*get_lld(), sizeof(double)}, get_array_pointer<matrix_row_major>(), py::cast(*this));
    }

    auto r = array.template mutable_unchecked<2>();
    for (auto i = 0; i < r.shape(0); ++i) {
      for (auto j = 0; j < r.shape(1); ++j) {
        r(i, j) = *mat.data(i, j);
      }
    }
  }

引数にpy::cast(*this)を指定しないと、distributed_matrixの成分が変更されない。その理由は不明。 py::cast(*this)は、呼び出し元のPython側に寿命を知らせるための物であるから、不要な気がしてしまうのだが。

array_tのコンストラクタで、全ての行列成分をコピーはできなさそうだ。少なくても、strideが異なる行列に対しては無理かもしれない。 Eigen::Mapを使わず、pybind11の中で完結している。

numpy.ndarrayでの行列major

ndarrayを作成してから、order='F'を変えると、行列成分は転置となる。

関数diagonalizeの書き方

採用したのは、C++版と同じ書き方

solver.diagonalize(mat, eigval, eigvec, params)

出力するオブジェクトは、左辺に書きたかった：

eigval, eigvec = solver.diagonalize(mat, params)

eigval = solver.diagonalize(mat, params)

そうしたい理由：

eigvalとeigvecの宣言をする必要がない。
Pythonらしい。

できない理由：固有値のみの場合と固有値・固有ベクトルを計算する場合について、内部の処理を分けることができないため。

Pythonでは、戻り値の数や型で関数をオーバーロードすることはできない。（この点はC++と同様）
一般に、固有値のみを計算する場合と固有値・固有ベクトルを計算する場合のルーチンは異なる。

古い情報

Boost.Pythonを用いたバインディングの使い方

PythonBinding_by_pybind11 - t-sakashita/rokko GitHub Wiki

Pythonバインディングの実装方法（pybind11を使用）

std::function

ベクトルの型

Eigen3

localized_vector

pybind11/eigen3.h

numpyとの連携

残っていること

mpi4py

boost::any

distributed_matrixの読み書き

numpy.ndarrayでの行列major

関数diagonalizeの書き方

古い情報

⚠️ GitHub.com Fallback ⚠️

PythonBinding_by_pybind11 - t-sakashita/rokko GitHub Wiki

Pythonバインディングの実装方法（pybind11を使用）

std::function

ベクトルの型

Eigen3

localized_vector

pybind11/eigen3.h

numpyとの連携

残っていること

mpi4py

boost::any

distributed_matrixの読み書き

numpy.ndarrayでの行列major

関数diagonalizeの書き方

古い情報

⚠️ **GitHub.com Fallback** ⚠️

⚠️ GitHub.com Fallback ⚠️