NumPyの2次元配列をC言語に渡す

# 概要
NumPyの2次元配列を（ctypesを使って）C言語の関数に渡す方法を紹介します。
本記事には書きませんが応用すれば3次元以上の多次元配列でもやりとりできます。

# 先に結論
libhoge.soという共有ライブラリに下記のような関数があったとします。
2行3列の行列を表示する関数です。
```c
void print2x3(int **a) {
    printf("%d,%d,%d\n", a[0][0], a[0][1], a[0][2]);
    printf("%d,%d,%d\n", a[1][0], a[1][1], a[1][2]);
}
```
この関数にNumPyの2次元配列を渡したいとします。
このとき下記のようなPythonスクリプトでうまくいきます。
```python
import numpy as np
from ctypes import *

# C言語の関数に渡すNumPy配列
src = np.array([[0,1,2], [3,4,5]], dtype=np.int32)

# C言語のライブラリ
libhoge = cdll.LoadLibrary("libhoge.so")

# NumPy配列をC言語の関数に渡す
src2 = src.copy()
f = src2.__array_interface__["data"][0]
d = src2.strides[0]
n = src2.shape[0]
c_src2 = f + d * np.arange(n)
c_src2 = c_src2.astype(np.uintp)
addr = c_src2.__array_interface__["data"][0]
addr = c_void_p(addr)
libhoge.print2x3(addr) # 表示
"""
出力
0,1,2
3,4,5
"""
```

# 詳細
## 予備知識 ： 配列の構造
### C言語
1次元配列の場合は単純で、すべての要素がメモリ空間上連続して格納されています。
2次元配列の場合は少し複雑です。C言語の2次元配列とは、1次元配列の先頭アドレスを要素とする1次元配列です。
例を示します。
下記のような2次元配列があったとしましょう。変数名はaとします。
```c
int a[2][3] = {{0,1,2}, {3,4,5}};
```
`a`はメモリ空間上に下図のように配置されます。
![carr.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3227325/41f4a93f-36df-368a-8195-046918ea2ab9.png)
（int型は4バイト、アドレス型も4バイト使うと仮定しています。）
（図中のアドレスの数字は適当です。）

### NumPy
1次元配列の場合も2次元配列の場合もすべての要素がメモリ空間上連続して格納されています。
例を示します。
また2次元配列の場合、**基本的には**同じ行の要素は連続して格納されています。
（転置などが絡むとその限りではありません。詳しくは[後述](#注意1-:-転置、reshape、スライス等と相性悪い)します。）
下記のような2行3列の2次元配列があったとしましょう。変数名はbとします。
```python
>>> b=np.array([[0,1,2],[3,4,5]], dtype=np.int32)
>>> print(b)
[[0,1,2],
 [3,4,5]]
```
このとき`b`の各要素はメモリ空間上に下図のように連続して配置されます。
![nparr.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3227325/dfd3c0e8-9019-4dbb-607a-dc7c5b324f8d.png)
（図中のアドレスの数字は適当です。）
![align.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3227325/2a4aaaeb-5b02-5027-800c-99ee308e76bd.png)

## C言語の2次元配列をNumPyで再現
NumPyの2次元配列の場合、メモリ空間上に全要素が連続して配置されているのでした。
すると下記の等差数列はC言語の2次元配列と同じになります。
* 初項：全要素の先頭アドレス
* 公差：各行のバイト数
* 項数：行数

![clike2darr.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3227325/b66c70e2-6f69-46f6-5b3f-a9d7b12b7f79.png)
さらに初項、公差、項数はそれぞれ次のよう得られます。

### 初項：全要素の先頭アドレス
```python
b.__array_interface__.["data"][0]
```
### 公差：各行のバイト数
```python
b.strides[0]
```
`strides`というメンバ変数は隣の行/列の要素までのバイト数を表したタプルです。
例に挙げている変数`b`の場合は下記のようになります。
```python
>>> print(b.strides)
(12,4)
```
![stride0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3227325/52c72a57-e34c-04c1-6476-ee9e2bc05f76.png)
![stride1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3227325/255b2f99-70d6-1d4a-2734-31f4d16e1d19.png)
### 項数：行数
```python
b.shape[0]
```
`shape`というメンバ変数は行数/列数を表したタプルです。
例に挙げている変数`b`の場合、2行3列の配列でしたので下記のようになります。
```python
>>> print(b.shape)
(2,3)
```

### 等差数列を得る
初項、公差、項数が分かりました。ここから下記のスクリプトで等差数列が得られます。
```python
f = b.__array_interface__.["data"][0]
d = b.strides[0]
n = b.shape[0]
c_b = f + d * np.arange(n)
```
（`np.arange(n)`は 0,1,…,n-1 という数列を返す関数です。）
### 型変換
まだ完成ではありません。
C言語の2次元配列はアドレスを要素とする配列でした。
しかし`c_b`の要素は整数型なのでアドレス型に変換します。
`astype`という関数で型変換できます。
```python
c_b = c_b.astype(np.uintp)
```
これでC言語の2次元配列と同じNumPy配列が得られました。

## 再現した配列をC言語に渡す
C言語の関数は配列を引数として受け取るとき、実際には配列の先頭アドレスを受け取ります。
そのため上の変数`c_b`をC言語の関数に配列を渡そうとしたら、`c_b`の要素の先頭アドレスを得る必要があります。
先の初項と同じく下記のようにすれば、変数`addr`として先頭アドレスを得られます。
```python
addr = c_b.__array_interface__["data"][0]
```
`addr`はPythonのint型です。C言語にアドレス型として渡すには型変換が必要です。
そこで下記のように`addr`を`ctypes.c_void_p`型に変換します。
```python
addr = c_void_p(addr)
```
これで準備が整いました。
あとはお目当てのC言語の関数に`addr`を渡せば完了です。
例えば冒頭の関数`print2x3`の場合、下記のように書けばよいです。
```python
libhoge.print2x3(addr)
```

## まとめ
以上をまとめると下記のようになります。
```python
# C言語の2次元配列を再現
f = b.__array_interface__.["data"][0]
d = b.strides[0]
n = b.shape[0]
c_b = f + d * np.arange(n)
c_b = c_b.astype(np.uintp)
# 再現した配列をC言語に渡す
addr = c_b.__array_interface__["data"][0]
addr = c_void_p(addr)
libhoge.print2x3(addr)
```
これは下記のような書き方もできます。
読みづらくなりますが、変数が減る＋実質2行なので処理がほんの少し軽くなります。多分。
```python
c_b = (
  b.__array_interface__.["data"][0] + 
  b.strides[0] * np.arange(b.shape[0])
).astype(np.uintp)
libhoge.print2x3(
  c_void_p(c_b.__array_interface__["data"][0])
)
```

### 注意1 : 転置、reshape、スライス等と相性悪い
「予備知識：配列の構造」の章にてNumPy配列はメモリ空間上に**基本的に**同じ行の要素が連続して格納されていると書きました。
しかし転置、関数`reshape`、スライスなどを使った場合、その限りではありません。
これらの処理はstridesとshapeを書き換えるのみで、配列の要素部分には触りません。
メモリ空間上で要素を並び替えると処理時間が長くなるので、「メモリ空間のどこが行列の何行何列目に対応するのか」という情報だけ書き換えようという発想に基づいています。
（下図はNumPy配列の転置処理後のイメージです。要素は書き換えず、矢印＝並び順のルールだけ書き換えることを示しています。）
![align_t.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3227325/8757bf52-31e1-c83f-a93e-07b0f31f7716.png)
とても賢い発想ですが、先に紹介した方法とは相性が悪いです。
「C言語の2次元配列をNumPyで再現」の章にある方法は、同じ行の要素が連続しているからこそ可能でした。
転置した配列などに対してはうまくいきません。
回避方法として、目当てのNumPy配列のメンバ関数`copy`を使って複製したのち、C言語に渡すという方法があります。
```python
# 配列を複製
b2 = b.copy()
# C言語の2次元配列を再現
f = b2.__array_interface__.["data"][0]
d = b2.strides[0]
n = b2.shape[0]
c_b2 = f + d * np.arange(n)
c_b2 = c_b2.astype(np.uintp)
# 再現した配列をC言語に渡す
addr = c_void_p(c_b2.__array_interface__["data"][0])
libhoge.print2x3(addr)
```

### 注意2 : NumPy配列がC言語によって書き換えられる
C言語の関数内に配列を書き換える処理がある場合、NumPy配列が書き換えられてしまいます。
例えばlibhogeという共有ライブラリに下記のような関数があったとしましょう。
2行3列の2次元配列の各要素に1を足す関数です。
```c
void plusone2x3(int **a) {
    for (int i = 0; i < 2; i++)
        for (int j = 0; j < 3; j++)
            a[i][j] += 1;
}
```
この関数に上の方法で変数`b`を渡します。
すると下記のように`b`の要素に1が足されます。
```python
print(b)
"""
出力
[[0,1,2],
 [3,4,5]]
"""
# 略
libhoge.plusone2x3(addr)
print(b)
"""
出力
[[1,2,3],
 [4,5,6]]
"""
```
C言語側でPythonの変数を書き換えるという挙動になり、単一責任の原則に反します。
必ずしもダメではありませんが、可読性やメンテナンス性を下げるので注意してください。

# 参考
https://stackoverflow.com/questions/22425921/pass-a-2d-numpy-array-to-c-using-ctypes/27737099