2019/3/24にmatplotlibのアニメーション機能に関する発表を行いました

2019/3/24にmatplotlibのアニメーションの機能に関する発表を行いました

2019/3/24に「はんなりPython+PyData Osakaの可視化特集会 」のイベントにて標題の発表を行いました。 イベントページはこちらです。 発表資料はいつものように下記のnotebookをRISEでスライドにしたものを使用しました。

matplotlibのアニメーション機能の紹介

自己紹介

  • 大橋宏正(@wrist)
  • PyData Osaka オーガナイザ
  • メーカー職(音響信号処理の業務に従事)

概要

  • matplotlibのアニメーション機能を紹介
  • 単純な方法の紹介(plt.pause)
  • Animationクラスを用いた方法の紹介
    • ArtistAnimation, FuncAnimation
    • TimedAnimation
  • 動画書き出し

なぜアニメーション描画を行いたいか

  • センサーデータの時系列変化を可視化したい
    • ex. PCM信号
    • c/c++バイナリが吐いたdumpデータの描画
In [1]:
import numpy as np
import scipy as sp
import matplotlib as mpl
import matplotlib.pyplot as plt
In [2]:
print(np.__version__)
print(sp.__version__)
print(mpl.__version__)
1.16.2
1.2.1
3.0.3
In [3]:
# osxの場合
%matplotlib osx

osxの場合の制約

  • アニメーション描画
    • backendに"osx"の指定が必要
  • "osx"の指定
    • FrameworkとしてビルドされたPythonが必要
    • 仮想環境(conda含む)のPythonはFrameworkビルドではない

osx(conda)の場合の対処方法

  • 公式FAQでの対処法
    • conda install python.app
    • pythonの代わりにpythonwを実行
    • pythonwを使うkernelを作成すればjupyter notebook上でも実行可能

pythonwを使用するカーネルの作成方法

  • 要ipykernel
$ ipython kernel install --user --name pythonw3 --display-name pythonw3
$ jupyter kernelspec list  # pythonw3カーネルのパスを確認
$ vim ~/.pyenv/versions/miniconda3-4.3.30/envs/py36/share/jupyter/kernels/pythonw3/kernel.json # kernel.json記載のpathをpythonwに変更

実際の実行結果

$ jupyter kernelspec list
Available kernels:
  julia-0.4    /Users/wrist/Library/Jupyter/kernels/julia-0.4
  julia_0.3    /Users/wrist/Library/Jupyter/kernels/julia_0.3
  python3      /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/share/jupyter/kernels/python3

$ ipython kernel install --user --name pythonw3 --display-name pythonw3
Installed kernelspec pythonw3 in /Users/wrist/Library/Jupyter/kernels/pythonw3

$ jupyter kernelspec list
Available kernels:
  julia-0.4    /Users/wrist/Library/Jupyter/kernels/julia-0.4
  julia_0.3    /Users/wrist/Library/Jupyter/kernels/julia_0.3
  pythonw3     /Users/wrist/Library/Jupyter/kernels/pythonw3
  python3      /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/share/jupyter/kernels/python3

# condaのpython3と同じディレクトリに移動
$ mv ~/Library/Jupyter/kernels/pythonw3 ~/.pyenv/versions/miniconda3-4.3.30/envs/py36/share/jupyter/kernels/
/Users/wrist/Library/Jupyter/kernels/pythonw3 -> /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/share/jupyter/kernels/pythonw3

$ jupyter kernelspec list
Available kernels:
  julia-0.4    /Users/wrist/Library/Jupyter/kernels/julia-0.4
  julia_0.3    /Users/wrist/Library/Jupyter/kernels/julia_0.3
  python3      /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/share/jupyter/kernels/python3
  pythonw3     /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/share/jupyter/kernels/pythonw3

# kernel.jsonに記載されているpathをpythonwに書き直す
$ vim ~/.pyenv/versions/miniconda3-4.3.30/envs/py36/share/jupyter/kernels/pythonw3/kernel.json

最も単純なアニメーション方法

  • plt.pause(second)を使って描画を繰り返す
    • plt.plot(),plt.pause(sec), plt.cla()を繰り返す
In [4]:
# 説明上の共通コード
plt.close("all")
dx = np.pi / 10.0
xs = np.arange(-np.pi, np.pi, dx)

fig = plt.figure(1)
ax = fig.add_subplot(111)
In [5]:
# アニメーションしないバージョン
for _ in range(100):
    ax.plot(xs, np.sin(xs)); xs += dx;
plt.show()
In [6]:
# 説明上の共通コード
plt.close("all")
dx = np.pi / 10.0
xs = np.arange(-np.pi, np.pi, dx)

fig = plt.figure(1)
ax = fig.add_subplot(111)
In [7]:
# アニメーションするバージョン
for n in range(100):
    ax.plot(xs, np.sin(xs)); xs += dx;
    plt.pause(0.1); plt.cla();  # 0.1s待ってclear
plt.show()
In [8]:
# アニメーションするバージョン
for n in range(100):
    ax.plot(xs, np.sin(xs)); xs += dx;
    plt.pause(0.1); ax.clear(); # ax.clearでも良い
plt.show()

plotのコストを下げる

  • ax.plotは毎回プロットし直すためコストが大きい
  • linesを取得して毎回set_dataでデータを設定
In [9]:
# 説明上の共通コード
plt.close("all")
dx = np.pi / 10.0
xs = np.arange(-np.pi, np.pi, dx)

fig = plt.figure(1)
ax = fig.add_subplot(111)
In [10]:
lines = ax.plot(xs, np.sin(xs))
for n in range(100):
    lines[0].set_data(xs, np.sin(xs)); xs += dx;
    ax.set_xlim((np.min(xs), np.max(xs)))  # 横軸範囲の再設定
    plt.pause(0.1)  # clear不要
plt.show()

plt.ion(), plt.draw()を使う方法

  • 昔の記事でよく出てくる方法
    • plt.ion() -> plt.show()
    • ax.plot() -> plt.draw() -> plt.cla()を繰り返す
    • インタラクティブモードをonにした上で再描画を都度命令
  • windowsでは過去は確かにこれでいけた(現状は未確認)
  • 手元のosxだとplt.pauseを入れないと動作しない
    • plt.ion()しておく必要は特にない?
In [11]:
print(plt.isinteractive())
True
In [12]:
# 説明上の共通コード
plt.close("all")
dx = np.pi / 10.0; xs = np.arange(-np.pi, np.pi, dx)
fig = plt.figure(1); ax = fig.add_subplot(111)
In [13]:
plt.ion(); plt.show();
for n in range(100):
    ax.plot(xs, np.sin(xs)); xs += dx
    plt.draw(); plt.pause(0.1); plt.cla();
In [14]:
print(plt.isinteractive())
True
In [15]:
plt.ioff()
In [16]:
print(plt.isinteractive())
False

matplotlibのアニメーション描画クラス

ArtistAnimation

  • 描画するプロットのリストを予め作成し与える
  • アニメーション描画クラスの中でも単純なクラス

class matplotlib.animation.ArtistAnimation(fig, artists, *args, **kwargs)

  • artists
    • 各フレームで描画するArtistコレクションのlist
In [17]:
import matplotlib.animation as anim
In [18]:
plt.close("all")
dx = np.pi / 10.0
xs = np.arange(-np.pi, np.pi, dx)
interval_ms = 100
images = []
In [19]:
fig = plt.figure(1); ax = fig.add_subplot(111)
for _ in range(100):
    images.append(ax.plot(xs , np.sin(xs))); xs += dx;
ani = anim.ArtistAnimation(fig, images, interval=interval_ms, repeat=False)    
plt.show()

問題点

  • 範囲制御ができていない
  • 理由
    • xsの範囲が都度増加しているため
  • 対策
    • xsを固定
    • x軸レンジが固定されてしまう
In [20]:
plt.close("all")
dx = np.pi / 10.0
xs = np.arange(-np.pi, np.pi, dx)
interval_ms = 100
images = []
In [21]:
fig = plt.figure(1); ax = fig.add_subplot(111)
for n in range(100):
    images.append(ax.plot(xs , np.sin(xs + n*dx)))
ani = anim.ArtistAnimation(fig, images, interval=interval_ms, repeat=False)    
plt.show()

FuncAnimation

  • 描画関数とデータ配列を与えてアニメーション
  • データ配列はgeneratorなどでも良い

class matplotlib.animation.FuncAnimation(fig, func, frames=None, init_func=None, fargs=None, save_count=None, **kwargs)

  • func : callable
    • def func(frame, *fargs) -> iterable_of_artists
    • framesで指定されたデータが渡される
  • init_func : callable, optional
    • def init_func() -> iterable_of_artists
  • frames : iterable, int, generator function, or None, optional
    • iterable: イテレーション可能なデータから描画データを供給
    • integer: 連番range(frames)を供給
    • generator function: def gen_function() -> objを指定
    • None: itertools.countを供給
In [22]:
import matplotlib.animation as anim
In [23]:
plt.close("all")
dx = np.pi / 10.0
xs = np.arange(-np.pi, np.pi, dx)
ys = np.sin(xs)
interval_ms = 100
In [24]:
# 呼び出しごとにデータを生成するジェネレータ
def gen_data():
    counter = 0
    dx = np.pi / 10.0
    xs = np.arange(-np.pi, np.pi, dx)
    while counter < 100:
        ys = np.sin(xs); xs += dx;
        counter += 1
        yield (xs, ys)
In [25]:
fig = plt.figure(1); ax = fig.add_subplot(111)
gline2d = ax.plot(xs, ys)[0]  # これを上書き

def update_data(data):
    xs, ys = data
    ax.set_xlim((np.min(xs), np.max(xs)))
    gline2d.set_data(xs, ys)
    return gline2d
    
ani = anim.FuncAnimation(fig, update_data, gen_data,
                         interval=interval_ms, repeat=False)
plt.show()
  • update_dataでx軸レンジを都度再設定できる

FuncAnimationを用いたアニメーション描画クラス

  • 先の例はグローバル変数を直接上書きしているため気持ち悪い
  • クラス化することで状態を同時に保持可能
  • generatorや描画関数はインスタンスメソッドにできる
In [26]:
class MyAnimation:
    def __init__(self):
        plt.close("all")
        self.fig = plt.figure(1)
        self.ax = self.fig.add_subplot(111)
        self.dx = np.pi / 10.0
        self.xs = np.arange(-np.pi, np.pi, self.dx)
        self.ys = np.sin(self.xs)
        self.interval_ms = 100
        self.line2d = self.ax.plot(self.xs, self.ys)[0]
        self.counter = 0
    
    def gen(self):
        while self.counter < 100:
            self.ys = np.sin(self.xs)
            self.xs += self.dx;
            self.counter += 1
            yield (self.xs, self.ys)
            
    def update_data(self, data):
        xs, ys = data
        self.ax.set_xlim((np.min(xs), np.max(xs)))
        self.line2d.set_data(xs, ys)
        return self.line2d

    def animate(self):
        ani = anim.FuncAnimation(self.fig, self.update_data, self.gen,
                         interval=self.interval_ms, repeat=False)
        plt.show()
In [27]:
MyAnimation().animate()

TimedAnimation

  • ArtistAnimationFuncAnimationが継承しているクラス
  • 自作クラスで継承し下記メソッドを定義して使用
    • self._draw_frame(self, frame)
    • self.new_frame_seq(self)
    • self._init_draw(self)
In [28]:
class MyTimedAnimation(anim.TimedAnimation):
    def __init__(self):
        plt.close("all")
        self.fig = plt.figure(1)
        self.ax = self.fig.add_subplot(111)
        self.dx = np.pi / 10.0
        self.xs = np.arange(-np.pi, np.pi, self.dx)
        self.ys = np.sin(self.xs)
        self.interval_ms = 100
        self.line2d = self.ax.plot(self.xs, self.ys)[0]
        self.counter = 0

        plt.show()
        anim.TimedAnimation.__init__(self, self.fig, self.interval_ms)
    
    def new_frame_seq(self):
        while self.counter < 100:
            self.ys = np.sin(self.xs)
            self.xs += self.dx;
            self.counter += 1
            yield (self.xs, self.ys)
            
    def _draw_frame(self, data):
        xs, ys = data
        self.ax.set_xlim((np.min(xs), np.max(xs)))
        self.line2d.set_data(xs, ys)
        return self.line2d

    def _init_draw(self):
        pass
In [29]:
MyTimedAnimation()
Out[29]:
<__main__.MyTimedAnimation at 0x11bb5dd30>

動画書き出し

  • saveメソッド
    • アニメーションを動画としてファイル保存可能
    • save(filename, writer="ffmpeg")
      • writerは文字列 or MovieWriterを継承したクラスを指定
  • to_html5_videoメソッド
    • videoタグを生成

save(filename, writer=None, fps=None, dpi=None, codec=None, bitrate=None, extra_args=None, metadata=None, extra_anim=None, savefig_kwargs=None)

In [30]:
plt.close("all")
dx = np.pi / 10.0
xs = np.arange(-np.pi, np.pi, dx)
interval_ms = 100
images = []
In [31]:
fig = plt.figure(1); ax = fig.add_subplot(111)
for _ in range(100):
    images.append(ax.plot(xs , np.sin(xs))); xs += dx;
ani = anim.ArtistAnimation(fig, images, interval=interval_ms, repeat=False)    
ani.save("./test.gif", "imagemagick")
In [32]:
!ls
animation_class.webp
animation_inheritance.webp
pydata_osaka_20190323_matplotlib_animation.ipynb
test.gif

In [33]:
# 要ffmpeg
from IPython.display import HTML
HTML(ani.to_html5_video())
Out[33]:

まとめ

  • matplotlibのアニメーション機能を紹介
    • 単純な方法
    • ArtistAnimation, FuncAnimation, TimedAnimation
    • 動画書き出しの方法

2018/12/21にPyData Osakaの梅キャンPython勉強会コラボ回でpybind11について発表しました

2018/12/21(金)にPyData Osakaの梅キャンPython勉強会コラボ回でpybind11について発表しました

ブログ記事にするのが遅くなりましたが2018/12/21に【大阪工業大学】特別回 梅キャンPython勉強会【梅田キャンパス】 でpybind11に関する発表を行いました。 下記はその際の資料になります。

pybind11の紹介

大橋宏正(PyDataオーガナイザー)

自己紹介

  • 大橋 宏正(@wrist)
  • PyDataオーガナイザー
  • 某メーカー勤務
    • 音響信号処理屋

pybind11とは

他のライブラリとの比較

  • Python C API
    • 自分で参照カウント制御が必要
    • 学習障壁が高い
  • ctypes
    • dllがあればビルド不要で使える
    • 返り値や引数の型を厳密に設定する必要有
  • cython
    • pyxの独自記法を覚える必要
    • 一旦pyxファイルをcにコンパイル
      • 更にcファイルをコンパイル
  • Boost.python
    • Boostが必要
    • 便利だが巨大でありコンパイルに時間がかかる

インストール

In [1]:
!pip install pybind11
Requirement already satisfied: pybind11 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (2.2.4)

ipybind

ipybindのインストール

  • リポジトリをクローン
  • python setup.py build
  • python setup.py install

使ってみる

  • ipybind経由でpybind11を使用
  • 最初にload_extでipybindを有効化
In [2]:
%load_ext ipybind

pybind11の基本

  • PYBIND_MODULE(モジュール名, m)
  • m.def(公開名, 関数へのポインタ, 説明)
In [3]:
%%pybind11

int add(int x, int y){
    return x + y;
}

PYBIND11_MODULE(myadd, m){
    m.def("c_add", &add, "Add two integers.");
}
In [4]:
c_add(1, 2)
Out[4]:
3
In [5]:
help(c_add)
Help on built-in function c_add in module pybind11_3e2e466:

c_add(...) method of builtins.PyCapsule instance
    c_add(arg0: int, arg1: int) -> int
    
    Add two integers.

In [6]:
# c_add("foo", "bar")
In [7]:
%%pybind11

PYBIND11_MODULE(myadd, m){
    m.def("lambda_add", [](int a, int b){ return a + b;},
         "Add two integers.");
}
In [8]:
lambda_add(1, 2)
Out[8]:
3

その他のコンパイル方法

  • 手動ビルド
  • distutils
    • setup.pyを書く
  • cmake
    • pybind11_add_module(myadd myadd.cpp)

単純なクラスをラップ

  • c++のクラスを定義
  • pybind11でラップ
    • py::classs_<C++クラス名>(m, 公開クラス名)
In [9]:
%%pybind11

namespace py = pybind11;

class Coordinate {
public:
    int x_;
    int y_;
    Coordinate(){ x_ = 0; y_ = 0; }
    Coordinate(int x, int y){ x_ = x; y_ = y; }
    float Norm(){ return sqrt(x_*x_ + y_*y_); }
};

PYBIND11_MODULE(coord, m) {
    py::class_<Coordinate>(m, "Coordinate")
        .def(py::init<>())
        .def(py::init<int, int>())
        .def_readonly("x", &Coordinate::x_)
        .def_readonly("y", &Coordinate::y_)
        .def_property_readonly("norm", &Coordinate::Norm);
}

コンストラクタの定義

  • ... .def(py::init<オーバーロード引数>())
Coordinate(){ x_ = 0; y_ = 0; }
    Coordinate(int x, int y){ x_ = x; y_ = y; }

は下記のように記載

py::class_<Coordinate>(m, "Coordinate")
        .def(py::init<>())
        .def(py::init<int, int>())

プロパティの定義

  • ... .def_readonly(公開変数名, メンバへのポインタ)
py::class_<Coordinate>(m, "Coordinate")
        .def_readonly("x", &Coordinate::x_)
        .def_readonly("y", &Coordinate::y_)

computed propetyの定義

  • ... .def_propety_readonly(公開変数名, 計算関数へのポインタ)
float Norm(){ return sqrt(x_*x_ + y_*y_); }

.def_property_readonly("norm", &Coordinate::Norm);
In [10]:
point = Coordinate(3, 4)
point
Out[10]:
<pybind11_c19ca2a.Coordinate at 0x10ba96df8>
In [11]:
print(point.x, point.y)
3 4
In [12]:
point.norm
Out[12]:
5.0

さらなる例

  • 演算子オーバーロード、python特殊メソッド
  • docstring、引数名、キーワード引数
  • *arg, **kwarg

演算子オーバーロード

  • 等値演算子
    • cppだとbool operator==(const Cls& rhs)
    • pythonだとCls.__eq__(self, other)
    • これらを対応付ける
In [13]:
int.__eq__(1, 1)
Out[13]:
True
In [14]:
%%pybind11

namespace py = pybind11;

class Coordinate2 {
public:
    int x_; int y_;
    Coordinate2(){ x_ = 0; y_ = 0; }
    Coordinate2(int x, int y){ x_ = x; y_ = y; }
    bool operator==(const Coordinate2& rhs){
        return (x_ == rhs.x_) && (y_ == rhs.y_);
    }
};

PYBIND11_MODULE(coord2, m) {
    py::class_<Coordinate2>(m, "Coordinate2")
        .def(py::init<>())
        .def(py::init<int, int>())
        .def_readonly("x", &Coordinate2::x_)
        .def_readonly("y", &Coordinate2::y_)
        .def("__eq__", &Coordinate2::operator==);
}
In [15]:
point1 = Coordinate2(3, 4)
point2 = Coordinate2(3, 4)
point3 = Coordinate2(1, 2)
In [16]:
point1 == point2
Out[16]:
True
In [17]:
point1 == point3
Out[17]:
False

演算子の代替記法

.def("__eq__", &Coordinate2::operator==);

はc++のクラスでoperator==が定義されていれば下記のようにも書ける

.def("__eq__", 
             [](const Coordinate2& self,
                const Coordinate2& other){
                 return self == other; 
             });

Python特殊メソッド

  • __repr__の定義を考える
  • 先の__eq__と同様に第一引数にselfを取るメソッドを定義
In [18]:
%%pybind11

namespace py = pybind11;

class Coordinate3 {
public:
    int x_; int y_;
    Coordinate3(){ x_ = 0; y_ = 0; }
    Coordinate3(int x, int y){ x_ = x; y_ = y; }
};

PYBIND11_MODULE(coord3, m) {
    py::class_<Coordinate3>(m, "Coordinate3")
        .def(py::init<>())
        .def(py::init<int, int>())
        .def_readonly("x", &Coordinate3::x_)
        .def_readonly("y", &Coordinate3::y_)
        .def("__repr__" , [](const Coordinate3& self){
            return std::string() + "(x, y) = (" + std::to_string(self.x_) + ", " + std::to_string(self.y_) + ")";
        });
} 
.def("__repr__", [](const Coordinate3& self){
    return std::string() + 
        "(x, y) = (" +
        std::to_string(self.x_) +
        ", " +
        std::to_string(self.y_) + ")";
});
In [19]:
point = Coordinate3(1, 2)
point
Out[19]:
(x, y) = (1, 2)

docstirng、引数名、キーワード引数

  • docstring
    • m.defの第三引数に文字列を渡す
  • 引数名
    • 第四引数以降にpy::arg("...")を渡す
  • キーワード引数
    • py::arg("...") = デフォルト値とすれば良い
In [20]:
%%pybind11

namespace py = pybind11;

class Coordinate4 {
public:
    int x_; int y_;
    Coordinate4(){ x_ = 0; y_ = 0; }
    Coordinate4(int x, int y){ x_ = x; y_ = y; }
    float DistanceTo(const Coordinate4& target){ 
        int dx = target.x_ - x_, dy = target.y_ - y_;
        return sqrt( dx*dx + dy*dy);
    }
};

PYBIND11_MODULE(coord4, m) {
    py::class_<Coordinate4>(m, "Coordinate4")
        .def(py::init<>())
        .def(py::init<int, int>())
        .def_readonly("x", &Coordinate4::x_)
        .def_readonly("y", &Coordinate4::y_)
        .def("distance_to", &Coordinate4::DistanceTo,
             "distance to target",
             py::arg("target"));
}
float DistanceTo(const Coordinate4& target){ 
        int dx = target.x_ - x_, dy = target.y_ - y_;
        return sqrt( dx*dx + dy*dy);
    }

を下記のように設定

.def("distance_to", &Coordinate4::DistanceTo,
             "distance to target",
             py::arg("target"));
In [24]:
point = Coordinate4(1, 2)
target = Coordinate4(3, 4)
point.distance_to(target)
Out[24]:
2.8284270763397217
In [40]:
point.distance_to?

*args, **kwarg

  • py::args (py::tupleのサブクラス)
  • py::kwargs (py::dictのサブクラス)
m.def("count_args", [](py::args a, py::kwargs kw) {
    py::print(a.size(), "args,", kw.size(), "kwargs");
});

その他の話題

  • C++とPython間のオブジェクトの取扱
  • 関数とコールバック
  • Numpy対応

C++とPython間のオブジェクトの取扱

  • C++とPythonのオブジェクトを混在してラッパーを定義可能
  • 3種類の方法
    • Python側でC++の型のインスタンスを引数に取る
    • C++側でPythonの型のインスタンスを引数に取る
    • C++の型とPythonの型を型変換する

Python側でC++の型のインスタンスを引数に取る

py::class_<Foo>(m, "Foo");
m.def("f1", [](const Foo& foo){ ...}

C++側でPythonの型のインスタンスを引数に取る

m.def("f2", [](py::list list){ ... }

C++の型とPythonの型を型変換する

m.def("f3", [](int x) { ... });
m.def("f4", [](const std::string& s) { ... });
m.def("f5", [](const std::vector<int>& v) { ... });

サポートされている型変換

  • スカラー値
    • integer types, float, double, bool, char
  • 文字列
    • std::string, const char *
  • タプル
    • std::pair<F, S>, std::tuple<...>
  • シーケンス
    • std::vector, std::list, std::array<T, n>
  • マップ
    • std::map<K, V>, std::unordered_map<K, V>
  • 集合
    • std::set, std::unordered_set
  • 関数
    • std::function<...>
  • Date/time
    • std::chrono::duration, std::chrono::time_point
  • Optional
    • std::optional, std::experimental::optional

関数とコールバック

  • std::functionでpythonの関数を受け取れる
In [38]:
%%pybind11

#include <pybind11/functional.h>

PYBIND11_MODULE(callback_test, m) {
    m.def("for_even",
          [](int n, std::function<void(int)> f) {
              for (int i = 0; i < n; ++i){ if (i % 2 == 0) f(i); }
          }
         );
}
In [39]:
def py_callback_func(x):
    print("called {0}".format(x))
    
for_even(10, py_callback_func)
called 0
called 2
called 4
called 6
called 8

numpy対応

  • #include <pybind11/numpy.h>
  • py::array_t<type>
    • ndarrayを受け取ることが可能
  • arr.ndim, arr.shape(n), arr.size()
  • バッファへのポインタ
    • arr.data(), arr.mutable_data()
    • インデックス指定も可能
  • 要素への直接アクセス
    • arr.unchecked(), arr.mutable_unchecked()
    • ref(i, j, k)としてアクセス可能
In [31]:
%%pybind11

#include <pybind11/numpy.h>

namespace py = pybind11;

PYBIND11_MODULE(numpy_bind, m) {
    m.def("add_2d",
          [](py::array_t<double> x, py::array_t<double> y){
              auto r = x.mutable_unchecked<2>();
              auto c = y.unchecked<2>();

            for(ssize_t i = 0; i < r.shape(0); i++){
                for(ssize_t j = 0; j < r.shape(1); j++){
                    r(i, j) += c(i, j);
                }
            }
          },  py::arg().noconvert(), py::arg().noconvert());
}
In [33]:
import numpy as np
In [34]:
x = np.array([[1.0, 2.0], [3.0, 4.0]])
y = np.array([[2.0, 3.0], [4.0, 5.0]])
add_2d(x, y)
x
Out[34]:
array([[3., 5.],
       [7., 9.]])

まとめ

  • Python向けのC拡張モジュールをC++11で作成するためのライブラリpybind11について紹介
  • ipybindを用いたJupyter Notebookからの実行方法の紹介
  • numpy連携を紹介

Pyroomacousticsの紹介

Pyroomacousticsの紹介

この記事はQiita Python Advent Calendar 2018の9日目の記事です。

Pyroomacousticsとは

Pyroomacousticsは音響アレイ信号処理に関するアルゴリズムの迅速な開発や評価を目的としたライブラリです。ライブラリは主に下記の3つの要素から成り立っています。

  • 2次元または3次元空間に配置された複数の音源と複数のマイクロホンに関するシミュレーション環境を迅速に構築するためのオブジェクト指向インタフェースの提供
  • 一般的な多面体形状の室内に対して室内インパルス応答を効率的に生成するための鏡像法の高速なC実装と、音源と受音点の間の伝搬のシミュレーション環境の提供
  • STFT, ビームフォーミング、方向推定(音源定位)、適応フィルタ、音源分離に対する一般的なアルゴリズムのリファレンス実装の提供

ここではデモ用notebookを参考に、使用方法や何ができるかについて説明していきます。 なおドキュメントはhttps://pyroomacoustics.readthedocs.io/en/pypi-release/index.html にあります。

インストール方法

pipで下記のようにインストールできます。また別途pipでcookiecutterをインストールしておくと、後述のようにpyroomacousitcs用のシミュレーションスクリプトを対話的に生成することが可能であるため便利です。

In [1]:
!pip install -U pyroomacoustics
Requirement already up-to-date: pyroomacoustics in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (0.1.20)
Requirement already satisfied, skipping upgrade: numpy in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from pyroomacoustics) (1.14.2)
Requirement already satisfied, skipping upgrade: scipy>=0.18.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from pyroomacoustics) (1.1.0)
Requirement already satisfied, skipping upgrade: Cython in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from pyroomacoustics) (0.29)
In [2]:
!pip install cookiecutter
Requirement already satisfied: cookiecutter in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (1.6.0)
Requirement already satisfied: jinja2>=2.7 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from cookiecutter) (2.10)
Requirement already satisfied: poyo>=0.1.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from cookiecutter) (0.4.2)
Requirement already satisfied: jinja2-time>=0.1.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from cookiecutter) (0.2.0)
Requirement already satisfied: requests>=2.18.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from cookiecutter) (2.19.1)
Requirement already satisfied: future>=0.15.2 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from cookiecutter) (0.16.0)
Requirement already satisfied: click>=5.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from cookiecutter) (7.0)
Requirement already satisfied: whichcraft>=0.4.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from cookiecutter) (0.5.2)
Requirement already satisfied: binaryornot>=0.2.0 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from cookiecutter) (0.4.4)
Requirement already satisfied: MarkupSafe>=0.23 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from jinja2>=2.7->cookiecutter) (1.0)
Requirement already satisfied: arrow in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from jinja2-time>=0.1.0->cookiecutter) (0.12.1)
Requirement already satisfied: chardet<3.1.0,>=3.0.2 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from requests>=2.18.0->cookiecutter) (3.0.4)
Requirement already satisfied: idna<2.8,>=2.5 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from requests>=2.18.0->cookiecutter) (2.7)
Requirement already satisfied: certifi>=2017.4.17 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from requests>=2.18.0->cookiecutter) (2018.10.15)
Requirement already satisfied: urllib3<1.24,>=1.21.1 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from requests>=2.18.0->cookiecutter) (1.22)
Requirement already satisfied: python-dateutil in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from arrow->jinja2-time>=0.1.0->cookiecutter) (2.7.5)
Requirement already satisfied: six>=1.5 in /Users/wrist/.pyenv/versions/miniconda3-4.3.30/envs/py36/lib/python3.6/site-packages (from python-dateutil->arrow->jinja2-time>=0.1.0->cookiecutter) (1.11.0)

Pyroomacousitcsの使用方法

まず下記ライブラリをインポートしておきます。

In [3]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile
from scipy.signal import fftconvolve
import IPython
import pyroomacoustics as pra

roomオブジェクトの生成

下記のように2次元座標を要素に持つndarrayを作成し、pra.Room.from_cornersに与えることで2次元または3次元の室内空間を表すオブジェクトを生成することが可能です。またroomオブジェクトのplotメソッドで空間を描画することが可能です。座標の単位はmです。

In [4]:
corners = np.array([[0,0], [0,3], [5,3], [5,1], [3,1], [3,0]]).T  # [x,y]
room = pra.Room.from_corners(corners)

fig, ax = room.plot()
ax.set_xlim([-1, 6])
ax.set_ylim([-1, 4]);

3次元空間については2次元空間をextrudeメソッドで指定した高さまで持ち上げます。

In [5]:
room = pra.Room.from_corners(corners)
room.extrude(2.)

fig, ax = room.plot()
ax.set_xlim([0, 5])
ax.set_ylim([0, 3])
ax.set_zlim([0, 2]);

なお、単純な正六面体の部屋であればShoeBoxメソッドを使ったほうが早いです。

In [6]:
room = pra.ShoeBox([2.0, 2.5])
fig, ax = room.plot()
ax.set_xlim([0, 3])
ax.set_ylim([0, 3])
Out[6]:
(0, 3)
In [7]:
room = pra.ShoeBox([2.0, 2.5, 2.0])
fig, ax = room.plot()
ax.set_xlim([0, 3])
ax.set_ylim([0, 3])
ax.set_zlim([0, 3]);

音源とマイクロホンの追加

roomオブジェクトに対してadd_sourceメソッドを用いることで指定した座標に音源を追加することが可能です。音源はndarrayである必要があるので、ここではscipy.ioのwavfileを用いて読み込んだ信号を用います

音源はデモ用ノートブックではCMU Arctic Databaseから取得した音声ファイルを用いています。 ここからダウンロードしたtar.gzを解凍したその中のwavディレクトリに入っています。後で気づきましたが、pyroomacoustics自体がデータセットのラッパーも提供しています。

In [8]:
# specify signal source
fs, signal = wavfile.read("work/arctic_a0010.wav")

# add source to 2D room
room = pra.Room.from_corners(corners, fs=fs)
room.add_source([1.,1.], signal=signal)

fig, ax = room.plot()

黒点が音源位置を表しています。追加した音源はこんな感じです。

In [9]:
IPython.display.Audio(signal, rate=fs)
Out[9]:

音源を追加した時点でroomオブジェクトにもfsが設定される模様です。試したところデフォルトでは8kHzが指定されているようでした。

In [10]:
room.fs
Out[10]:
16000

この空間に更に環状のマイクロホンアレイを追加します。 ここではpra.circular_2D_arrayを用いて6素子からなる中心座標[2.0, 2.0]、半径0.1の環状アレイの各マイクロホンの座標を計算します。 各マイクロホンアレイの座標をfsとともにpra.MicrophoneArrayに渡すことでマイクロホンアレイオブジェクトを生成し、これをroomオブジェクトのadd_microphone_arrayメソッドにわたすことでマイクロホンアレイを室内空間へと追加します。

In [11]:
R = pra.circular_2D_array(center=[2.,2.], M=6, phi0=0, radius=0.1)
room.add_microphone_array(pra.MicrophoneArray(R, room.fs))

fig, ax = room.plot()

xがマイクロホンを表します。

鏡像法を用いた室内インパルス応答の生成

音源とマイクロホンアレイが追加された室内では、音源が各マイクロホンへとどのように伝搬するかを表すインパルス応答を生成することができます。 ここではこのインパルス応答を鏡像法に基づいて生成します。 鏡像法とはざっくり言うと室内の壁面での反射波は壁面に対して音源と対象な位置に鏡像が存在すると考え、音源位置から放射された信号が壁面に伝達する度に鏡像を計算して、最終的には実音源とこれらの鏡像からマイクロホンへと到来する波の重ね合わせがマイクロホンにおける観測信号となるシミュレーション手法です。反射回数を何回まで考慮するかや、壁面での吸音率、取り扱うサンプリングレートなどによって鏡像法で推定されるインパルス応答の精度が変わってくることが知られています。

ここでは3次元の室内空間を生成し、[1.0, 1.0, 0.5]の位置に音源を追加、[3.5, 2.0, 0.5]と[3.5, 2.0, 0.5]の位置にマイクロホンアレイを追加します。 更に、roomオブジェクトの生成時に最大何次の反射までを考慮するかを表すmax_orderと壁面での反射時の吸音率を表すabsorptionが指定されています。max_orderが8であるため8次の反射までが考慮され、吸音率が0.2であるため反射率が0.8、すなわち壁面での反射時には音波が0.8倍されることを意味しています。

この状態でroomオブジェクトに対してimage_source_modelメソッドを呼ぶと、鏡像法に基づき室内の反射を模擬することが可能です。use_libroomのオプションは鏡像法のC実装を用いるかどうかを制御するオプションのようです。

なお、以下では3次反射までを描画しています。

In [12]:
# specify signal source
fs, signal = wavfile.read("work/arctic_a0010.wav")

# set max_order to a low value for a quick (but less accurate) RIR
room = pra.Room.from_corners(corners, fs=fs, max_order=8, absorption=0.2)
room.extrude(2.)

# add source and set the signal to WAV file content
room.add_source([1., 1., 0.5], signal=signal)

# add two-microphone array
R = np.array([[3.5, 3.6], [2., 2.], [0.5,  0.5]])  # [[x], [y], [z]]
room.add_microphone_array(pra.MicrophoneArray(R, room.fs))

# compute image sources
room.image_source_model(use_libroom=True)

# visualize 3D polyhedron room and image sources
fig, ax = room.plot(img_order=3)
fig.set_size_inches(18.5, 10.5)

一度鏡像法のシミュレーションを実行すると、下記のように室内インパルス応答を描画することが可能になります。

In [13]:
room.plot_rir()
fig = plt.gcf()
fig.set_size_inches(20, 10)

なおこの室内インパルス応答(RIR)はroomオブジェクトの下にrirとして格納されています。

In [14]:
room.rir
Out[14]:
[[array([0., 0., 0., ..., 0., 0., 0.])],
 [array([0., 0., 0., ..., 0., 0., 0.])]]

pyroomacousticsではこの室内インパルス応答を音源に対して畳み込んだ場合のシミュレーションを行い、マイクロホンでの観測音源を生成することができます。 roomオブジェクトに対してsimulateメソッドを呼ぶことで、室内インパルス応答を畳み込んだ観測信号がroom.mic_array.signalsとして生成されます。ここではマイクロホンが2つ存在しているため、2つの観測信号が生成されています。

In [15]:
room.simulate()
print(room.mic_array.signals.shape)
(2, 55372)

実際に聞いてみて元音源と比較することが可能です。

In [16]:
# original signal
print("Original WAV:")
IPython.display.Audio(signal, rate=fs)
Original WAV:
Out[16]:
In [17]:
print("Simulated propagation to first mic:")
IPython.display.Audio(room.mic_array.signals[0,:], rate=fs)
Simulated propagation to first mic:
Out[17]: