pythonistaがC++を使って配列を関数で扱う

目次

はじめに

どうも,かわつあんです. 私のブログでは,これまでに"ctypesでC++の関数をpythonから呼び出したい"というモチベーションのもとpythonista目線で記事を投稿してきました. 今回は,C++で行列を関数で扱ってみます.これができればやっと本題に進めます.内容はC++を知ってる方には簡単すぎると思いますので,御注意の程よろしくお願いします.

配列を引数にして関数内で処理をしてスカラーあるいは配列を返すようなことはデータ解析では日常茶飯事です. pythonだったら何も考えずにできますがCやC++では参照やポインタの概念を理解して注意して扱わないといけません.ちなみにvectorクラスを使えば,比較的pythonに近く書けるので今回はそちらを採用しています. (実は最近になって,pythonの関数も参照渡しだと知りました.pythonでも参照渡しであることに注意が必要なようですね.)

この記事では  3\times 3行列の行列式を求める関数,および逆行列を求める関数をC++及びpythonで自作してこれらの比較をしてみます. なお,参照やポインタの概念の詳細はいろんな方がまとめていらっしゃいますのでしますのでご了承ください.

引数が行列,戻り値がスカラー(数値)の関数

例題として, A =  \begin{bmatrix} a_{11} ,a_{12} ,a_{13} \\ a_{21}, a_{22}, a_{23} \\ a_{31}, a_{32}, a_{33} \end{bmatrix} 行列式 \mathrm{det} Aを計算してみます.  3\times3行列の行列式は数式的にはサラスの公式より,以下の通りにかけます.

 \mathrm{det}A= a_{11}a_{22}a_{33} + a_{12}a_{23}a_{31} + a_{13}a_{32}a_{21} - a_{13}a_{22}a_{31} - a_{11}a_{23}a_{32} -a_{12}a_{21}a_{33}

例えば, A =  \begin{bmatrix} 1,2,3 \\ 3,5,3\\3,8,3 \end{bmatrix} 行列式 \mathrm{det}A = 18 となります.

これをpythonで実装してみます.

import numpy as np

def det(A):
    det_A = A[0,0]*A[1,1]*A[2,2] + A[0,1]*A[1,2]*A[2,1] + A[1,0]*A[0,2]*A[2,1] \
          - A[0,2]*A[1,1]*A[2,0] - A[0,0]*A[1,2]*A[2,1] - A[0,1]*A[1,0]*A[2,2]
    return det_A 

a = np.array([[1,3,3],
              [2,5,8],
              [3,3,3]])
print("det A = {}".format(det(a)) )

出力結果

det A = 18

このように比較的に直感的にコーディングできますね.

一方で,C++ではvectorクラスを使うと以下のようにかけます.

#include <iostream>
#include <vector>
using namespace std;

int det (vector<vector<int>>A){
    int det_A;
    det_A = A[0][0]*A[1][1]*A[2][2] + A[0][1]*A[1][2]*A[2][1] + A[1][0]*A[0][2]*A[2][1]
          - A[0][2]*A[1][1]*A[2][0] -A[0][0]*A[1][2]*A[2][1] -A[0][1]*A[1][0]*A[2][2];
    return det_A ;
};

int main(){
    vector<vector<int>> a {{1,3,3},{2,5,8},{3,3,3}};
    cout << "det A = " << det(a) << endl;
}

出力結果

det A = 18

pythonと似たようにコーディングできましたね. ここまではC++初心者でも比較的理解しやすいと思います. しかしC++やCで難しいのは次のトピックで使う参照渡しですので,まだ安心しないでください.

引数が行列,戻り値が行列?の関数

さて次のトピックについては"戻り値が行列?の関数"というように?マークで強調してみましたが...

先に言っておくとC++(及びC)において戻り値を配列にすることはできません.

pythonばっかりやってきた僕は最初にこれみてびっくりしました.だったら配列の関数実装できないじゃん...て思っていましたが... ちゃんと方法はあります. それが参照渡しというやり方です.

参照渡しとは,ざっくり言うと関数内で計算した配列を記憶しているメモリの場所を引数として関数に渡すことです.関数内でメモリに格納された値を更新することで,関数処理が終えたあともメモリに計算後の値が残ります.fortranでいうサブルーチンのように使えばいいということですね. まあ,初学者の方は辛いと思いますが詳しくは他の文献などをご覧ください.

話はそれましたが,配列を計算して配列を返す関数を実装していきます. 今回は簡単のため3\times3行列の足し算を行ってみます. 行列の足し算は簡単で各成分ごとに足し算するだけです. 具体的には,以下の通りです.

 A =  \begin{bmatrix} 1,1,1 \\ 1,1,1, \\ 1,1,1 \end{bmatrix}, \quad B =  \begin{bmatrix} 2,2,2 \\ 2,2,2, \\ 2,2,2 \end{bmatrix}, \quad C = A+B =  \begin{bmatrix} 3,3,3 \\ 3,3,3, \\ 3,3,3 \end{bmatrix}

まずは,これをpythonで実装していきます.numpyの配列なら何も考えず足すだけでは?というツッコミはなしでお願いします. あくまで今回は愚直に,各成分同士の足し算を成分に持つ配列を出力する関数を使います.

実装したコードは以下の通りです.

import numpy as np

def sum_matrix(A,B):
    C = np.zeros((3,3))
    for i in range(3):
        for j in range(3):
            C[i,j] = a[i,j] + b[i,j]
    return C

a = np.array([[1,1,1],
              [1,1,1],
              [1,1,1]])
b = np.array([[2,2,2],
              [2,2,2],
              [2,2,2]])
c = sum_matrix(a,b)
print("C = A+B =\n{}".format(c))

実行結果

C = A+B =
[[3. 3. 3.]
 [3. 3. 3.]
 [3. 3. 3.]]

まあ,愚直にかけばこんなもんでしょう.

C++だとどういう風にかけるのでしょうか? 参照渡しを使って書くと次の通りです.

#include <iostream>
#include <vector>
using namespace std;

void sum_matrix (vector<vector<int>>A, vector<vector<int>>B, vector<vector<int>>& C){ //vector<vector<int>> C ではないことに注意
    for (int i=0; i<3;i++){
        for (int j=0; j<3;j++){
            C[i][j] = A[i][j] + B[i][j];
        }
    }
};

int main(){
    vector<vector<int>> a {{1,1,1},{1,1,1},{1,1,1}};
    vector<vector<int>> b {{2,2,2},{2,2,2},{2,2,2}};
    vector<vector<int>> c {{0,0,0},{0,0,0},{0,0,0}};

    sum_matrix(a, b, c);
    cout << "C = A+B ="<< endl;
    for (int i=0; i<3;i++){
        for (int j=0; j<3;j++){
            cout << c[i][j] << ", ";
        }
        cout << endl;
    }
}

出力結果

C = A+B =
3, 3, 3, 
3, 3, 3, 
3, 3, 3, 

この場合,計算結果をいれる配列を初めに作って,その配列の参照を引数に入れます. 配列Cの参照を&Cとして引数に入れ,関数内で配列Cのメモリの値を更新しています. ちなみに,もし

vector<vector<int>>&C

の部分を

vector<vector<int>>C

としてしまうと,値は更新されず出力結果が以下の通りになってしまいます.

C = A+B =
0, 0, 0, 
0, 0, 0, 
0, 0, 0, 

これは,配列の参照ではなく値渡しと言う別の操作になってしまうためです.なお,ポインタ渡しでもコードを少し変えれば関数で配列計算等ができますが,ここでは省略します.

このように,参照渡し(あるいはポインタ渡し)の話は難しいので詳しくは他の文献等をみてください. Qiitaにいい感じの記事がまとまっていました. 参考にどうぞ.

qiita.com

さいごに

今回は,C++言語で配列を関数で扱うために,pythonistaとしてpythonのコードと比較しながら実装しました. 参照渡しという概念を意識しないといけないのでc++のコードは少し難しいですね. 私のブログをみて少しでもpythonistaの方が参照渡しの雰囲気を味わえたのならそれだけで幸いです.

次回はctypesを使ってc++の関数を呼び出す内容を記事にしようと思います.

追伸 いろいろ試行錯誤したのですが,現状ctypesではvectorの扱いが難しい,numpy配列のメモリの扱いがおかしい,など問題がある様なので,この後の記事はかけそうにないです.申し訳ない...

ではまた(✿╹◡╹)