Pythonパフォーマンスチューニング：FlaskでWSGIラインプロファイラ（wsgi_lineprof）を使う

PythonのWebアプリケーションで汎用的に使えるプロファイラを探していた所、wsgi_lineprofが良さそうだったのでFlaskアプリで使用した例を共有させて頂きます。

## 背景

先日社内ISUCONが開催されまして、ISUCON無知だったのですがとりあえずPython勢として参加してみました。
結果はあまり振るわなかったのですが、特にプロファイリングすら出来なかったことがずっと心残りだったため、今回調べてみたという流れです。

## wsgi_lineprof

PythonのWSGI ミドルウェアとして使えるラインプロファイラで、Django,Pyramid,Flask,Bottle などのWSGI互換のアプリケーションやフレームワークと組み合わせて利用できるとのこと。ちなみにPyCon2017のLTで発表されたもののようです。
[WSGI ミドルウェアとして使えるラインプロファイラを作った話](https://blog.ymyzk.com/2016/12/line-profiler-as-a-wsgi-middleware/)

## 早速使ってみる

ということで早速使ってみます。
使用するプログラムは先日会社のテックブログに投稿させていただいた[ベイズロジスティック回帰](https://tech.fusic.co.jp/posts/2019-08-22-2019-08-23-bayesian-logistic-regression-python)のコードを流用しています。

```python
import jinja2
import flask
import math
import numpy as np
import sys
from scipy.stats import bernoulli

from wsgi_lineprof.middleware import LineProfilerMiddleware
from wsgi_lineprof.filters import FilenameFilter, TotalTimeSorter


LAYOUT = """
<html>
  <title>Guestbook</title>
  <body>
    <ul>
    {% for entry in entries %}
      <li><h2>{{ entry.name }}</h2>
        {{ entry.value }}
      </li>
    {% endfor %}
    </ul>
  </body>
</html>
"""

np.random.seed(1234)

# xの次元数
M = 2

# xのデータ数
N = 25

# 変分パラメータの学習率
alpha = 1.0e-4

# 変分推論の繰り返し回数
max_iter = 1000

def sigmoid(a):
    s = 1 / (1 + math.exp(-a))
    return s

# 変分推論
def VI(Y, X, M, Sigma_w, alpha, max_iter):

    def rho2sig(rho):
        return np.log(1 + np.exp(rho))

    def calc_d_mu(Y, X, Sigma_w_inv, mu, rho, W):
        # 第一項は0なため不要
        term2 = np.dot(Sigma_w_inv, W)
        term3 = 0
        W_tr = W.T
        term3 = sum([ -(Y[i] - sigmoid(np.dot(W_tr, X[:, i]))) * X[:, i] for i in range(n)])

        return term2 + term3

    # diag gaussian for approximate posterior
    m, n = X.shape
    mu = np.random.randn(m)
    rho = np.random.randn(m)
    Sigma_w_inv = np.linalg.inv(Sigma_w)

    for i in range(max_iter):
        # sample epsilon
        ep = np.random.randn(m)
        W_tmp = mu + rho2sig(rho) * ep

        # calculate gradient
        d_mu = calc_d_mu(Y, X, Sigma_w_inv, mu, rho, W_tmp)
        d_rho = ((d_mu * ep) - (1 / rho2sig(rho))) * (1 / (1+np.exp(-rho)))

        # update variational parameters
        mu = mu - alpha * d_mu
        rho = rho - alpha * d_rho

    return mu, rho


app = flask.Flask('bayesian-logistic-regression')


@app.route('/')
def vi():
    # xの事前分布の共分散行列Σ
    Sigma_w = 100.0 * np.eye(M)
    X = 2 * np.random.rand(M, N) - 1.0
    W_truth = np.random.multivariate_normal(np.zeros(M), Sigma_w)
    W_truth_tr = W_truth.T
    Y = np.array([bernoulli.rvs(sigmoid(np.dot(W_truth_tr, X[:, i]))) for i in range(N)])
    mu, rho = VI(Y, X, M, Sigma_w, alpha, max_iter)

    entries = [{'name': 'mu', 'value': mu}, {'name': 'rho', 'value': rho}]
    return jinja2.Template(LAYOUT).render(entries=entries)


if __name__ == '__main__':
    filters = [
        FilenameFilter("test_app.py"), # このプロフラムファイル名です
        lambda stats: filter(lambda stat: stat.total_time > 0.0001, stats),
    ]
    app.config['PROFILE'] = True
    with open("lineprof.log", "w") as f:
        app.wsgi_app = LineProfilerMiddleware(app.wsgi_app, stream=f, filters=filters)
        app.run(host='0.0.0.0', port=5000)
```

上記のコードを `test_app.py`として、
`python3 test_app.py`でFlaskを立ち上げて、ブラウザで`http://localhost:5000/`にアクセスすると以下の結果が得られました。

<img width="459" alt="スクリーンショット 2019-08-30 9.10.29.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/120225/dfa09af8-4f58-7b00-ecf8-f623ebbfa46f.png">

wsgi_lineprofの結果ファイルは以下のようになっています。

```
Time unit: 1e-09 [sec]

File: test_app.py
Name: vi
Total time: 0.437176 [sec]
  Line      Hits         Time  Per Hit  % Time  Code
====================================================
    85                                          @app.route('/')
    86                                          def vi():
    87                                              # xの事前分布の共分散行列Σ
    88         1        42667  42667.0     0.0      Sigma_w = 100.0 * np.eye(M)
    89         1        37846  37846.0     0.0      X = 2 * np.random.rand(M, N) - 1.0
    90         1       953424 953424.0     0.2      W_truth = np.random.multivariate_normal(np.zeros(M), Sigma_w)
    91         1         1457   1457.0     0.0      W_truth_tr = W_truth.T
    92         1      6082760 6082760.0     1.4      Y = np.array([bernoulli.rvs(sigmoid(np.dot(W_truth_tr, X[:, i]))) for i in range(N)])
    93         1    394951989 394951989.0    90.3      mu, rho = VI(Y, X, M, Sigma_w, alpha, max_iter)
    94                                          
    95         1         1169   1169.0     0.0      entries = [{'name': 'mu', 'value': mu}, {'name': 'rho', 'value': rho}]
    96         1     35104474 35104474.0     8.0      return jinja2.Template(LAYOUT).render(entries=entries)

File: test_app.py
Name: <listcomp>
Total time: 0.00604374 [sec]
  Line      Hits         Time  Per Hit  % Time  Code
====================================================
    92        26      6043742 232451.6   100.0      Y = np.array([bernoulli.rvs(sigmoid(np.dot(W_truth_tr, X[:, i]))) for i in range(N)])
    93                                              mu, rho = VI(Y, X, M, Sigma_w, alpha, max_iter)
    94                                          
    95                                              entries = [{'name': 'mu', 'value': mu}, {'name': 'rho', 'value': rho}]
    96                                              return jinja2.Template(LAYOUT).render(entries=entries)
    97                                          
    98                                          
    99                                          if __name__ == '__main__':
   100                                              filters = [
   101                                                  FilenameFilter("test_app.py"),
   102                                                  lambda stats: filter(lambda stat: stat.total_time > 0.0001, stats),
   103                                              ]

File: test_app.py
Name: sigmoid
Total time: 0.0231688 [sec]
  Line      Hits         Time  Per Hit  % Time  Code
====================================================
    41                                          def sigmoid(a):
    42     25025     17901666    715.4    77.3      s = 1 / (1 + math.exp(-a))
    43     25025      5267177    210.5    22.7      return s

File: test_app.py
Name: VI
Total time: 0.385176 [sec]
  Line      Hits         Time  Per Hit  % Time  Code
====================================================
    46                                          def VI(Y, X, M, Sigma_w, alpha, max_iter):
    47                                          
    48         1          833    833.0     0.0      def rho2sig(rho):
    49                                                  return np.log(1 + np.exp(rho))
    50                                          
    51         1          832    832.0     0.0      def calc_d_mu(Y, X, Sigma_w_inv, mu, rho, W):
    52                                                  # 第一項は0なため不要
    53                                                  term2 = np.dot(Sigma_w_inv, W)
    54                                                  term3 = 0
    55                                                  W_tr = W.T
    56                                                  term3 = sum([ -(Y[i] - sigmoid(np.dot(W_tr, X[:, i]))) * X[:, i] for i in range(n)])
    57                                          
    58                                                  return term2 + term3
    59                                          
    60                                              # diag gaussian for approximate posterior
    61         1         1316   1316.0     0.0      m, n = X.shape
    62         1        25150  25150.0     0.0      mu = np.random.randn(m)
    63         1         2371   2371.0     0.0      rho = np.random.randn(m)
    64         1       228259 228259.0     0.1      Sigma_w_inv = np.linalg.inv(Sigma_w)
    65                                          
    66      1001       599453    598.9     0.2      for i in range(max_iter):
    67                                                  # sample epsilon
    68      1000      3516197   3516.2     0.9          ep = np.random.randn(m)
    69      1000      5893163   5893.2     1.5          W_tmp = mu + rho2sig(rho) * ep
    70                                          
    71                                                  # calculate gradient
    72      1000    358329384 358329.4    93.0          d_mu = calc_d_mu(Y, X, Sigma_w_inv, mu, rho, W_tmp)
    73      1000     12519616  12519.6     3.3          d_rho = ((d_mu * ep) - (1 / rho2sig(rho))) * (1 / (1+np.exp(-rho)))
    74                                          
    75                                                  # update variational parameters
    76      1000      2239694   2239.7     0.6          mu = mu - alpha * d_mu
    77      1000      1819522   1819.5     0.5          rho = rho - alpha * d_rho
    78                                          
    79         1          471    471.0     0.0      return mu, rho

File: test_app.py
Name: rho2sig
Total time: 0.00702795 [sec]
  Line      Hits         Time  Per Hit  % Time  Code
====================================================
    48                                              def rho2sig(rho):
    49      2000      7027949   3514.0   100.0          return np.log(1 + np.exp(rho))

File: test_app.py
Name: calc_d_mu
Total time: 0.351921 [sec]
  Line      Hits         Time  Per Hit  % Time  Code
====================================================
    51                                              def calc_d_mu(Y, X, Sigma_w_inv, mu, rho, W):
    52                                                  # 第一項は0なため不要
    53      1000      5815585   5815.6     1.7          term2 = np.dot(Sigma_w_inv, W)
    54      1000       339332    339.3     0.1          term3 = 0
    55      1000       686803    686.8     0.2          W_tr = W.T
    56      1000    344064367 344064.4    97.8          term3 = sum([ -(Y[i] - sigmoid(np.dot(W_tr, X[:, i]))) * X[:, i] for i in range(n)])
    57                                          
    58      1000      1015033   1015.0     0.3          return term2 + term3

File: test_app.py
Name: <listcomp>
Total time: 0.313714 [sec]
  Line      Hits         Time  Per Hit  % Time  Code
====================================================
    56     26000    313714155  12065.9   100.0          term3 = sum([ -(Y[i] - sigmoid(np.dot(W_tr, X[:, i]))) * X[:, i] for i in range(n)])
```

概ね期待する結果が得られたと思います。逆行列の取得処理が遅いなぁといったことがわかりますね！

## その他のプロファイラについて

実は本記事を書くにあたり様々なプロファイラを試してみて、結果wsgi_lineprofが一番良かったかなという感じでした。自身の備忘録もかねて、他プロファイラを使った感想を記載いたします。（ラインプロファイラに限らず良さげなものを探していました）

- cProfile
  - Pythonの標準モジュールであり、手軽に使える
  - 結果をファイルに書き出せば `snakeviz` を使って簡単に可視化できる
  - 結果をフィルタする方法がわからず断念。（普通に使うFlaskパッケージのファイル等が大量に出力されてしまう…）
    - before_request, after_request等でうまくフックできるかと思ったけど、うまくいかなかった。
- pyflame
    - uberが作っているプロファイラ
    - Dockerでセットアップして使ってみた所、`Unexpected ptrace(2) exception: Failed to PTRACE_PEEKDATA (pid XXXX, addr 0x7f87d33ce030): Input/output error`のエラーが大量に出て修正できず断念。
    - install自体は公式や色々な人のDockerfileを参考にしてうまくいってるはずなのだが…。
- py-spy
    - pipで手軽にインストール可能で、実行中のPython（Flask）のPIDを指定すればさくっとプロファイリングできた。
    - ボトルネックにさくっとあたりを付けるのによさそう。
    - 可視化する場合、svgファイルに書き出すので検索性に厳しいものがあると感じた。
- vprof
    - 可視化がリッチで良さげ。
    - 純粋なPythonスクリプトの実行では期待する結果が得られた。
    - Flaskで実行する方法がわからず断念。
        - 公式にFlaskのサンプルがあるが、`signal only works in main thread` というエラーが出て動かない…。
- Flask-profiler
    - [公式リポジトリ](https://github.com/muatik/flask-profiler)のプロファイリング結果サンプルをみて **これが欲しかった！**ってなった。 
    - インストール＆導入も簡単にできた。
    - しかし、元々ISUCON出発で調べ始めていたので、Flaskにロックされるのはどうなの…？と思ってちょっと触って終わった。
![dashboard_screen.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/120225/69b8f686-d5c1-10d5-8c9c-3625d353ee50.png)



## 参考にさせていただいたサイト

- [PythonのWSGIラインプロファイラを試してみた(wsgi_lineprof)](http://www.denzow.me/entry/2017/09/18/162154)
- [Flaskの表示ページのプロファイラをとってみた](https://dev.classmethod.jp/server-side/taking_the_profiler_of_flask/)
- [Pythonコードのプロファイリング](https://shkh.hatenablog.com/entry/2013/08/16/221205)
- [Rust製プロファイリングツールpy-spyを試す](https://qiita.com/uokada/items/473305dc4df513584a15)