Kei Minagawa's Blog

皆川圭(@keimina)のブログ、Pythonで試したことを書いていきます

SciPy Japan 2019 に参加します

SciPy Japan 2019 に参加します。

www.scipyjapan2019.scipy.org

開催は4/23, 4/24の二日間、場所は東京。チュートリアルトーク、ポスター発表が行われるようです。チュートリアルも二日間に渡って行われるようです。

英語のリスニングができずに、ほとんど理解できないかもしれず、無謀なのですが、せっかくの機会ですので、参加してみます。果たしてどうなることか、、、

Wikipediaの文書を使って文の長さの平均を求める

Wikipediaの文書を使って文の長さの平均値と最頻値を求めました。機械学習勉強会 in 新潟で発表した内容です。Wikipedia の約60億文字以上のXMLファイルから、独自に定義した文を抜き出し、その文に対して統計量を求めました。使用言語はいつも通りPythonです。

1. 動作環境について

コードは以下の環境で動作確認を行っています。

  • CPU:6コア
  • メモリ:256GB
  • HDD:50GB

※環境はGCE(Google Compute Engine)で行っています。実行する場合はメモリ量に注意してください。

2. Wikipedia 文書のアーカイブについて

Wikipedia の文書は以下からダウンロードできます。 "jawiki-latest-pages-articles.xml.bz2"の中にあるXMLが今回解析に使用したファイルになります。

dumps.wikimedia.org

3. 発表スライド

発表スライドは、一番下に置いておきましたので、そちらをご覧ください。

4. やってみた感想

まず、このような分析は、自然言語処理の分野では頻度分析と呼ばれるらしいです。統計ですね。最近は統計のビジネス書をコツコツ読んで勉強していたたため、おかげでサンプルデータの分布を見た瞬間パラメタ推定できると気づくことができました。
パラメタ推定は本で読んでいたけど実際やってみると、ニューラルネットワークなど必要とせず非常にシンプルなのですが、強力でした。今回は、分布が対数正規分布にフィットしたのですがその理由を考察するときに、仮定した関数の性質が役に立ちました。結果、個人的に納得感のある考察を得ることができました(対数正規分布に従うのはなぜか?の考察です)。
最後ですが、Wikipediaを日常的に使用していたため、その中に法則性が存在するという事実に驚きました。おそらく解析しないと一生気付かなかったと思います。個人的に、その法則性の存在を客観的に示せたので非常に有意義だったと思います。また、求めた最頻値の値自体にも、価値があると思います。この値はWikipediaで遭遇する可能性が一番高い文の文字列の長さなので、それは、もしかすると我々が最も親しんでいる文字列長かもしれません。そのため24文字を目指して文を書くとわかりやすさの改善に繋がるつかもしれません。ビッグデータの時代ですが、手付かずのデータが山ほどあるとのことなので、まだ誰も知らない事実がデータの中に大量に眠っていると推測されます。今後それらのデータの中から、想像もしなかった事実が出てくることを期待しています。

5. 発表会で上がった質問や意見

発表会で上がった質問をここにメモしておきます。これらの質問や意見については、今後扱うかもしれませんし、扱わないかもしれません。個人的には、最頻値を求めて十分満足してしまったため、あとやったら面白そうなことは、コーパス教師なし学習クラスタリングなど)をさせてテキストマイニングするくらいしか思いつかないです。あと人工無脳

1. 他の分布で推定してみたらどうなるのか?例えばポアソン分布でやってみたらどうなるか?
→ ググったら一番最初に出てきたキーワードが対数正規分布だったのでそれでやってみた、ポアソン分布ではパラメタ推定していない。(ただし、対数正規分布でパラメタ推定を行ったきっかけはどうであれ、実際は、考察に記載したように対数正規分布を仮定するもっともらしい理由があるため、ポアソン分布や他の分布でパラメタ推定は行う必要性を感じなかった)
2. 集めた文(コーパスとも言える)を今後どういったことに使えるのか?
→ オートエンコーダによる教師なし学習(クラスタリング)を行いテキストマイニングを行うことが考えられる。
3. Wikipediaではなくニュースサイトなどの文の場合は違う分布になるのか?
→ 本などをみるとコーパスの分野が変わると違う分布になったり、あるコーパスで学習した機械学習のモデルが、別の分野のコーパスを入力して予測すると精度が落ちるなどの事例がある様子。そのため、ニュースサイトの分布とは異なるものになる可能性が高い。
4. XMLの読み込みにかかる時間はどれくらいか
→ 約15分くらい
5. 正規表現の適用にかかる時間はどれくらいか
→ 6コアCPUで分散処理しているため早い、約5分くらい。CPUを増やせばほぼリニアに早くなる(←発表では10分くらいと言ったかもしれませんが訂正します)

発表スライド

発表したスライドになります。正規表現が冗長になっており(スクエアブラケットの中のバックスラッシュのこと)、もっとコンパクトにできますが、そこは認識しています。時間がないため修正せずにあげています。ご指摘ありましたらコメントをお願いいたします。
以上。

f:id:keimina:20190225193006j:plain
スライド1
f:id:keimina:20190225193114j:plain
スライド2
f:id:keimina:20190225193134j:plain
スライド3
f:id:keimina:20190225193152j:plain
スライド4
f:id:keimina:20190225193210j:plain
スライド5
f:id:keimina:20190225193306j:plain
スライド6
f:id:keimina:20190225193403j:plain
スライド7
f:id:keimina:20190225193417j:plain
スライド8
f:id:keimina:20190225193434j:plain
スライド9
f:id:keimina:20190225193448j:plain
スライド10
f:id:keimina:20190225193501j:plain
スライド11
f:id:keimina:20190225193520j:plain
スライド12
f:id:keimina:20190225193609j:plain
スライド13

Pythonで2次元データ作成ーポリゴンに収まる点の集合

前回の続きです。今回は、前回に引き続き、ポリゴンに収まる点の集合を作成します。アルゴリズムは簡単で、ランダムに点を打ち、それがポリゴンの内側であれば点を残す、ポリゴンの外側であれば点を消すを繰り返すだけです。「ポリゴンの内側か?」を判定する関数を作成し、その関数に点の位置を次々と問い合わせればよいでしょう。
※ここでは「ポリゴン」は多角形の意味で使っています。

前回の記事:
keimina.hatenablog.jp

前々回の記事:
keimina.hatenablog.jp


N角形(N>=3)の領域内の点の集合を得ることができますので、前回前々回の記事で作成した点の集合も、この記事の方法だけで作成することができるようになります。

1. 点がポリゴンの領域の内側にある点か判定する方法

点がポリゴンの領域の内側にあるか判定する方法は以下に色々書かれていました。

stackoverflow.com

上記の回答の一つを見たところ、点がポリゴンの内側にあるとき、その点とポリゴンを構成する点から算出される角度の総和は必ず360度になるらしいです。図を書いてイメージをつかみます。以下が、点A,点B,点C,,,点Jを結んだポリゴンの内部に点Pがあった時の図になります。点Bから点Cあたりでガタガタと移動する箇所があります。点Bから点Cに移動する時のように角度θiがマイナスになるときもあるのですが、その後のガタガタで相殺されます。θiの総和が360°のとき左回りに1回転できたということになります。※かなりざっくりした説明で厳密ではありません。

f:id:keimina:20190126202512p:plain
図1. 点がポリゴンの領域の内側にあるとき

一方、点がポリゴンの外側にあるとき、360度にはならないと書かれています。これも図を書いてイメージをつかみます。θiの箇所についてみて見ると前半は右回りで、後半は左回りとなっています。右回りと、左回りで相殺されますので、θiの総和は0°になることがわかります。

f:id:keimina:20190126202605p:plain
図2. 点がポリゴンの領域の外側にあるとき

ポリゴンのような多角形の図形は、円の一つだと考えられます。そのため多角形を円とみなして考えるとわかりやすいです。円形の壁の中に、自分がいたとして、ぐるっと見回して壁に囲まれていたら、円の中だとわかります。一方、ぐるっと見回して左には壁があったけど右には壁がなかったなどの場合は壁の中でないことがわかります。算出される角度の総和が360°だったというのは、壁をたどると左周りに1回転できたということと同じです。1回転できなかった場合は、壁を辿っても1回転できなかったということですので、壁がつながっていない、すなわち壁の内側ではないということがわかります。
今回は本アルゴリズムを使用します。(アルゴリズムがちゃんと機能するのかも確認します)

2. 角度の求め方

角度を求めるには以下の数学公式を用います。

cosθ = (PA・PB)/(|PA|*|PB|)
θ = arccos(cosθ)

※「・」は内積を、arccosはアークコサインです。
詳細は以下を参照してください。

ja.wikipedia.org

3. 角度を求めるための関数を作成する

numpyで角度を求める関数を実装すると以下になります。

コード:

import numpy as np
from numpy import inner
from numpy.linalg import norm

def angle(P, A, B):
    '''
    args
    ----
        P: np.array(x, y))
        A: np.array(x, y))
        B: np.array(x, y))
    returns
    -------
        theta: x (degree)
    '''
    # 全体をシフトさせる
    PA = A - P
    PB = B - P
    # cosθを求める
    costheta = inner(PA, PB)/norm(PA)/norm(PB)
    # 角度を求める(度)
    theta = np.arccos(costheta)/np.pi * 180
    return theta

# テスト
P = np.array((0,0))
A = np.array((1,1))
B = np.array((1,0))
print(angle(P, A, B))

出力:

45.00000000000001

4. 時計回りか、反時計回りかを判定する

3. によって角度の絶対値はもとめることができました。しかし、反時計回りか時計回りか区別できません。時計回りか、反時計回りか区別できないと、角度の総和が常に360°になってしまうため、区別できるようにする必要があります。そこで、点A(a, b)と点B(c, d)と点P(0, 0)によってできるベクトルPAとPBのなす角θについて、時計回りと反時計回りを以下のように定義します。

(i) d > (b/a) * c
ならば、反時計回り方向の角度

(ii) d < (b/a) * c
ならば時計回り方向の角度

※説明を簡単にするため点P(0, 0)としています。点Pが(0, 0)以外のときでも、全体をシフトさせて点Pが(0, 0)になるように点Aと点Bを移動させても、角度に関する情報は失われないので、P(0, 0)でも問題ありません。

また(i)(ii)の不等号の式の両辺にaをかけて、数式を変形すると以下のように書き直せます。

(i)' a*d - b*c > 0
(ii)' a*d - b*c < 0

図で表すと以下になります。点Bが緑色の領域にあればPABの角度は反時計回り方向。緑色の領域にあれば、時計回り方向というふうに定義しています。
なぜ、この定義をするのかは、説明しづらいですが、図1、図2を実際に手でなぞって見ると、まず回転方向が必要なことがわかります。また0°<θ<180°のときは足し算で-180°<θ<0°のときは引き算のように計算方法を変えなければいけないこともわかると思います。そのため、このような定義が必要であるということがわかります。

f:id:keimina:20190126202901p:plain
図3. この記事における、反時計回り、時計回りの定義

5. 角度を求めるための関数2(反時計回り、時計回りの情報を追加)

3.で作成したangle関数を、4.の定義を使用して、反時計周りならプラスを、時計回りならマイナスの角度を返すように修正します。

コード:

import numpy as np
from numpy import inner
from numpy.linalg import norm

def angle2(P, A, B):
    '''
    args
    ----
        P: np.array(x, y))
        A: np.array(x, y))
        B: np.array(x, y))
    returns
    -------
        theta: x (degree)
    '''
    # 全体をシフトさせる
    PA = A - P
    PB = B - P
    # cosθを求める
    costheta = inner(PA, PB)/norm(PA)/norm(PB)
    # 角度を求める(ラジアン)
    theta = np.arccos(costheta)
    # 角度を求める(度)
    theta = np.arccos(costheta)/np.pi * 180
    # 
    if PA[0]*PB[1] - PA[1]*PB[0] > 0:
        # 条件(i)'に該当するときはなにもしない(ここではthetaはプラス)
        pass
    else:
        # 条件(ii)'に該当するときはthetaをマイナスにする
        theta *= -1.0
    return theta

# テスト
P = np.array((0,0))
A = np.array((1,1))
B = np.array((1,0))
print(angle2(P, A, B))

出力:

-45.00000000000001

6. 点がポリゴンの領域の内側にある点かどうかを判定する関数を作成しテストする

準備が整いましたので、点がポリゴンの領域の内側にある点かどうかを判定する関数を作成し、機能するかテストします。

点(x, y)が1.のアルゴリズムを使用しポリゴンの領域内の点であるかどうかを判定するis_inside_polygon関数を実装します。以下のようになります。
この関数は[(x1, y1), (x2, y2),,,]のような点(x, y)の集合pointsと点p(x, y)を受け取り、pointsからなるポリゴンの領域内に点pが含まれるか判定します。

コード:

import random
import matplotlib.pyplot as plt

import numpy as np
from numpy import inner
from numpy.linalg import norm

def angle2(P, A, B):
    '''
    args
    ----
        P: np.array((x, y))
        A: np.array((x, y))
        B: np.array((x, y))
    returns
    -------
        theta: float (dgree)
    '''
    # 全体をシフトさせる
    PA = A - P
    PB = B - P
    # cosθを求める
    costheta = inner(PA, PB)/norm(PA)/norm(PB)
    # 角度を求める(度)
    theta = np.arccos(costheta)/np.pi * 180
    # 時計回りか、反時計回りか判定し、thetaをプラスかマイナスにして返す
    if PA[0]*PB[1] - PA[1]*PB[0] > 0:
        # 条件(i)'に該当するときはなにもしない(ここではthetaはプラス)
        pass
    else:
        # 条件(ii)'に該当するときはthetaをマイナスにする
        theta *= -1.0
    return theta

def is_inside_polygon(points, P):
    '''
    args
    ----
        points: list of np.array((x, y))
        P: np.array((x, y))
    returns
    -------
        theta: x (degree)
    '''
    points = [np.array(p) for p in points]
    P = np.array(P)
    angle_total = 0.0
    for A, B in zip(points, points[1:] + points[0:1]): # Notice! "points[1:] + points[0:1]" means something like "points[1:].extend(points[0:1])"
        angle_total += angle2(P, A, B)
    if np.isclose(abs(angle_total), 360):
        # 360°の場合はポリゴンの内側なのでTrueを返す
        return True
    else:
        # ポリゴンの外側
        return False

# 三角形
#points = ((0,0), (1,0), (0, 1))

# リボン?
#points = ((1.0, 0.0), (0.5, 0.5), (1.0, 1.0), (-1.0, 1.0), (-0.5, 0.5), (-1.0, 0.0))

# 十字
points = ((1.0, -3.0), (1.0, -1.0), (3.0, -1.0), (3.0, 1.0), (1.0, 1.0), (1.0, 3.0), (-1.0, 3.0), (-1.0, 1.0), (-3.0, 1.0), (-3.0, -1.0), (-1.0, -1.0), (-1.0, -3.0))
points = [(i[0]/3, i[1]/3) for i in points]

N = 5000
for _ in range(N):
    # 点を打つ場所をランダムに決める(ただし、-1<x<1, -1<y<1とする)
    x = (random.random() - 0.5) * 2
    y = (random.random() - 0.5) * 2
    p = (x, y)
    # 決めた点(x, y)が三角形の内側にあればグラフに点を打つ
    if is_inside_polygon(points, p):
        plt.plot(x, y, color='k', marker='.')

# グラフを表示する
plt.gca().set_xlim(-2,2)
plt.gca().set_ylim(-2,2)
plt.gca().set_aspect('equal')
plt.show()

7. 出力結果

上から順に三角形、リボン、十字形になります。ランダムに点を打っているのですが、pointsで定義された多角形の内側にある点だけ残ることがわかります。
※形を変えるには、多角形の定義(#points... の箇所)を変えてください。

f:id:keimina:20190126203803p:plain
図4. 三角形に収まる点の集合

f:id:keimina:20190126204014p:plain
図5. リボン形に収まる点の集合

f:id:keimina:20190126204048p:plain
図6. 十字形に収まる点の集合

8. まとめ

ポリゴンに収まる点の集合を作成できました。

以上です。

Pythonで2次元データ作成ー三角形に収まる点の集合

前回の続きです。今回は、前回に引き続き、Pythonで三角形に収まる点の集合を作成します。アルゴリズムは簡単で、ランダムに点を打ち、それが三角形の内側であれば点を残す、三角形の外側であれば点を消すを繰り返すだけです。「三角形の内側か?」を判定する関数を作成し、その関数に点の位置を次々と問い合わせればよいでしょう。ここでは三角形は直角二等辺三角形ということにします。

前回の記事:
keimina.hatenablog.jp

1. 三角形の領域の定義

実装する前に使用する三角形の形を明確にします。以下のように使用する三角形の領域を定義します。

三角形の領域の定義:
 x-y座標上の点(x, y)で「x > 0 かつ y > 0 かつ y < -x + 1」を満たす点の集合

この条件を満たす点(x, y)は3つの点(0, 0)と(1, 0)と(0, 1)を結ぶとできる三角形の領域を構成します。なぜそうなるかについては、以下の①②③を参考にしてください。

①「x > 0 かつ y > 0」 の領域は x-y座標の右上の領域すべてを表します。
②「y < -x + 1』の領域 は 直線 y = -x + 1 より下側の領域をすべてを表します。
③「x > 0 かつ y > 0 かつ y < -x + 1」 の領域は上記①②の領域が重なった領域であり、三角形になります。

①②③を図で表すと以下になります。

f:id:keimina:20190120142434j:plain
図1. 三角形の領域の定義

2. 三角形の領域内の点であるかどうかを判定する関数を作成する

点(x, y)が上記1.で定義したの三角形の領域内の点であるかどうかを判定する関数を実装します。以下のようになります。

def is_inside_triangle(x, y):
    if x >= 0 and \
       y >= 0 and \
       y <= -x + 1:
        # 「xが0以上の領域」かつ「yが0以上の領域」かつ「y=-x+1の直線より下の領域」
        # すなわち点(x,y)が三角形の内側ならTrue
        return True
    else:
        # 三角形の外側ならFalse
        return False


is_inside_triangle関数は引数で与えられた点(x, y)の情報が1.で定義した三角形の内側にあるかどうかを判定します。三角形の内側であればTrue、それ以外はFalseを戻します。

3. コード1(ライブラリを極力使わない場合)

2.で作成した関数を使用することで、三角形に収まる点の集合を求めることができます。ライブラリを極力つかわないで全体のコードを書くとこんな感じになるかと思います。前回の記事のコードにある、is_inside_square関数 を 今回作成したis_inside_triangle関数 に置き換えただけです。

※matplotlibは可視化のため使用

import random
import matplotlib.pyplot as plt

def is_inside_triangle(x, y):
    if x >= 0 and \
       y >= 0 and \
       y <= -x + 1:
        # 「xが0以上の領域」かつ「yが0以上の領域」かつ「y=-x+1の直線より下の領域」
        # すなわち点(x,y)が三角形の内側ならTrue
        return True
    else:
        # 三角形の外側ならFalse
        return False

N = 1000
for _ in range(N):
    # 点を打つ場所をランダムに決める
    x = random.random()
    y = random.random()
    # 決めた点(x, y)が三角形の内側にあればグラフに点を打つ
    if is_inside_triangle(x, y):
        plt.plot(x, y, color='k', marker='.')

# グラフを表示する
plt.gca().set_xlim(-2,2)
plt.gca().set_ylim(-2,2)
plt.gca().set_aspect('equal')
plt.show()

4. コード2(Numupyを使用して書いた場合)

Numpyを使用して書くと以下のようになると思います。これも、前回の記事のコードにある、is_inside_square関数 を 今回作成したis_inside_triangle関数 に置き換えただけです。

import numpy as np
import matplotlib.pyplot as plt

def is_inside_triangle(x, y):
    if x >= 0 and \
       y >= 0 and \
       y <= -x + 1:
        # 「xが0以上の領域」かつ「yが0以上の領域」かつ「y=-x+1の直線より下の領域」
        # すなわち点(x,y)が三角形の内側ならTrue
        return True
    else:
        # 三角形の外側ならFalse
        return False

# 関数をnumpy化する
is_inside_triangle = np.vectorize(is_inside_triangle)

N = 1000

# ランダムに点を打つ場所を決める
x = np.random.random(N)
y = np.random.random(N)

# 三角形の内側にある点だけを抜き出す
conditions = is_inside_triangle(x, y)
x = x[conditions]
y = y[conditions]

# グラフに点を打つ
plt.plot(x, y, color='k', marker='.')

# グラフを表示する
plt.gca().set_xlim(-2,2)
plt.gca().set_ylim(-2,2)
plt.gca().set_aspect('equal')
plt.show()

5. 出力

以下のようなグラフが出力されると思います。Pythonで三角形に収まる点の集合を作成できました。

f:id:keimina:20190120145112p:plain
図2. 出力

以上です。

Pythonで2次元データ作成ー四角形に収まる点の集合

Pythonで四角形に収まる点の集合を作成します。アルゴリズムは簡単で、ランダムに点を打ち、それが四角形の内側であれば点を残す、四角形の外側であれば点を消すを繰り返すだけです。「四角形の内側か?」を判定する関数を作成し、その関数に様々な点を与え、戻り値がTrueのときの点を集めればよいということになります。ここでは四角形は正方形ということにします。一つ目のソースコードは以下になります。

1. ライブラリを極力つかわないで書いた場合

ライブラリを極力つかわないで書くとこんな感じになるかと思います。(matplotlibは可視化のため使用)

コード1:

import random
import matplotlib.pyplot as plt

def is_inside_square(x, y):
    if 0 <= x <= 1 and 0 <= y <= 1:
        # 「xが0以上1以下の領域」かつ「yが0以上1以下の領域」すなわち点(x,y)が四角形の内側ならTrue
        return True
    else:
        # 四角形の外側ならFalse
        return False

N = 1000
for _ in range(N):
    # 点を打つ場所をランダムに決める
    x = random.random()
    y = random.random()
    # 決めた点(x, y)が三角形の内側にあればグラフに点を打つ
    if is_inside_square(x, y):
        plt.plot(x, y, color='k', marker='.')    # グラフに点を打つ

# グラフを表示する
plt.gca().set_xlim(-2,2)
plt.gca().set_ylim(-2,2)
plt.gca().set_aspect('equal')
plt.show()

f:id:keimina:20190118223908p:plain

random.random()を使用すると0から1までの少数をランダムに選んでくれます。これを利用することで、点をランダムに打つことができます。基本的には、ランダムに打った点の位置がなんらかの形状の内側にあるか、外側にあるか判定できれば、その形状を点で表すことができると思います。

2. 行数を短く書いた場合

上記コードに加え、以下の二つ目のコードも載せておきます。結果は同じです。random.random()の戻り値が[0, 1)の範囲内の値なので、四角形の内側であるかの判定は、実はこの問題では必要ないのです。またnumpyを使いfor文をなくしています。

コード2:

import numpy as np
import matplotlib.pyplot as plt

N = 1000

# ランダムに点を打つ場所を決める
x = np.random.random(N)
y = np.random.random(N)

# グラフに点を打つ
plt.plot(x, y, color='k', marker='.')

# グラフを表示する
plt.gca().set_xlim(-2,2)
plt.gca().set_ylim(-2,2)
plt.gca().set_aspect('equal')
plt.show()

3. 最終的なコード

上記のコードをみると短くはなりましたが、四角形の内側かを判定するis_inside_square関数があったほうがよいと思いました。なぜなら、関数があったほうが何をしたいのかがわかるのと、四角形ではなく三角形や他の形で同様のことをしたいときに役立つと思ったためです。

そのため最終的なコードは以下のようになりました。
コード3:

import numpy as np
import matplotlib.pyplot as plt

def is_inside_square(x, y):
    if 0 <= x <= 1 and 0 <= y <= 1:
        # 「xが0以上1以下の領域」かつ「yが0以上1以下の領域」すなわち点(x,y)が四角形の内側ならTrue
        return True
    else:
        # 四角形の外側ならFalse
        return False

# 関数をnumpy化する
is_inside_square = np.vectorize(is_inside_square)

N = 1000

# ランダムに点を打つ場所を決める
x = np.random.random(N)
y = np.random.random(N)

# 四角形の内側にある点だけを抜き出す
conditions = is_inside_square(x, y)
x = x[conditions]
y = y[conditions]

# グラフに点を打つ
plt.plot(x, y, color='k', marker='.')

# グラフを表示する
plt.gca().set_xlim(-2,2)
plt.gca().set_ylim(-2,2)
plt.gca().set_aspect('equal')
plt.show()

関数をnumpy化するのにnp.vectorizeを使用しています。

これでPythonで四角形に収まる点の集合を作成できました。同様の方法で三角形や円の内側に収まる点の集合も作成できます。

以上です。おやすみなさい。

numpyの2次元配列の非ゼロ領域を囲む四角形の情報を取得する方法について理解する

本記事は、Stackoverflowにに記載されているnumpyの2次元配列の非ゼロ領域を囲む四角形の情報を取得する方法について理解を深めるための記事です。

2. numpyの2次元配列の非ゼロ領域を囲む四角形の情報を取得する方法

以下のURLを参照してください。関数自体は非常にコンパクトで5行くらいで記載されています。本記事では、この5行くらいのコードについて理解を深めるために解説を行います。

stackoverflow.com

上記によると、2次元配列の非ゼロ領域を囲む四角形の情報を取得するコードは以下になります。

def bbox2(img):
    rows = np.any(img, axis=1)
    cols = np.any(img, axis=0)
    rmin, rmax = np.where(rows)[0][[0, -1]]
    cmin, cmax = np.where(cols)[0][[0, -1]]

    return rmin, rmax, cmin, cmax

3. 解説

なぜ、上記のコードで非ゼロ領域を囲む四角形の情報が得られるかを解説します。

以下の図1にアルゴリズムの動作をグラフで表しました。np.anyを使って行と列を圧縮しています。np.anyを使うと非ゼロをみつけたら計算が打ち切られるため限界まで高速化できていると思います。行方向に圧縮したものは行のサマリ(要約)に、列方向に圧縮したものは列のサマリであると解釈するとわかりやすいかもしれません。そして、その行のサマリの左端と右端、列のサマリの上端と下端の位置が我々が欲する四角形の情報です。それらをnp.whereを使用して求めています。ここではnp.whereは引数で与えられた1次元配列のTrueのインデックス番号を取得するために使用しています。np.whereの戻り値として2次元配列が戻ってきますが、インデックス番号0しか使わないため

np.where(x)[0]

のように0を指定しています。また、

np.where(x)[0][0]
np.where(x)[0][-1]

のようにインデックスに0と-1を指定しているのは左端と右端、上端と下端を得るためです。ちなみに、ここで使用しているnp.whereはnp.nonzeroと同じです*1。最終的に図1の右下にある式が得られます。

f:id:keimina:20190110192941j:plain
図1

図1の右下にある式とStackoverflowにある式とでは見た目が違うので(実際の処理内容は同じですが)式を変形して同じになるか確認します。以下の図2のように、どんどんワンライナーでコードの行数を圧縮します。すると最終的にStackoverflowの回答と同じような形の式が得られます。

f:id:keimina:20190109221605j:plain
図2

4. 注意点

np.whereで得ようとしている配列(2.に記載した関数bbox2の引数img)がすべて0やFalseのときインデックスエラーになるはずですので注意が必要かと思います。

以上です。おやすみなさい。

*1:詳細は公式サイトを参照してください https://docs.scipy.org/doc/numpy-1.15.0/reference/generated/numpy.where.html

nvprofコマンド と NVIDIA Visual Profiler(GUIアプリ) を使って GPUのプロファイリングを行い、tf.nn.conv2d の計算時間を計測してみる

以前の記事で tf.nn.conv2d を複数GPUでやったら遅かったという記事を書いたのですが、「オーバーヘッドも含めて計測しているため遅いだろう」というご指摘があったため、もう一度検証してみました。最近参加している勉強会の主催者の方から、NVIDIA Visual Profiler という GPU のプロファイラについてのアドバイスをいただいたため、それを使い、処理時間の計測を行いました。

以前の記事:
keimina.hatenablog.jp

プロファイリングのために使用するのはnvprof コマンドと NVIDIA Visual Profiler (GUIアプリ) の二つです。nvprof コマンドは NVIDIAGPU のプロファイルをおこなうためのコマンドで、プロファイリングの結果をファイルとして保存してくれます。NVIDIA Visual Profiler はその名前の通り、ビジュアルなプロファイリングを行うGUIアプリです。

1. 参考にしたページ

参考にしたページは以下になります。NVIDIA Visual Profiler のインストール、実行方法など大変参考になりました。

qiita.com

2. nvprof コマンドについて

nvprof コマンドは NVIDIAGPU のプロファイルをおこなうためのコマンドで、プロファイリングの結果をファイルとして保存してくれます。コマンドのインストール方法は割愛します。私の場合、GCE(Google Compute Engine)に標準でインストールされていました。参考として使用したGCEのイメージを以下に記載します。

使用したGCEのイメージ:

Intel® optimized Deep Learning Image: TensorFlow 1.12.0 m15 (with Intel® MKL-DNN/MKL and CUDA 10.0)
A Debian based image with TensorFlow (With CUDA 10.0 and Intel® MKL-DNN, Intel® MKL) plus Intel® optimized NumPy, SciPy, and scikit-learn.

3. NVIDIA Visual Profiler について

nvprof が出力するプロファインリング結果(ファイル)を読み取って可視化してくれるアプリです。こちらもインストール方法は割愛しますが、私はインストール時に設定するチェックボックスやらなんやらは全てデフォルトでひたすら「Next」を押してインストールしました。Macbook pro にインストールしたのですが、以下の場所にプロファイラがインストールされてました。

プロファイラがインストールされた場所:

/Developer/NVIDIA/CUDA-10.0/libnvvp/nvvp.app

f:id:keimina:20181227221202p:plain
図1. nvvp.app

4. 手順概要

nvprof コマンドを用いてプロファイリングした結果をファイルとして取得、次に、あとで説明する NVIDIA Visual Profiler(GUIアプリ) でそのファイルを食わせて可視化します。

5. 手順詳細 (1/2) - nvprof コマンドの実行

以下にプロファイリングの手順を書きます。
※以下の説明では GPUを実行する環境(私の場合 GCEのインスタンス上)で nvprof はすでにインストールされているものとします。また、NVIDIA Visual Profiler も nvprof のファイルを読み込むための PC(私の場合 Macbook pro) にインストールされているものとします。

以下のようにGPU環境(私の場合GCEのインスタンス上)で 計測したいスクリプトを nvprof コマンドを経由して実行させます。

nvprof -o out-file.nvp python3 sample-script.py

“-o” で出力ファイル名を指定しています。これにより sample-script.py が実行され、同時に out-file.nvp が作成されます。この out-file.nvp にプロファイリングの情報が入っています。

6. 手順詳細 (2/2) - NVIDIA Visual Profiler による可視化

nvprof コマンドで作成した *.nvp ファイルを NVIDIA Visual Profiler という GUI アプリにインポートし読み込み結果を可視化します。GPU の情報をタイムラインで表示してくれるので、時間のかかっている場所がどこなのかなどを見ることができます。まず、out-file.nvp ファイルがサーバーにある人(私)は ローカル(私の場合 Macbook pro)に移しておきます。次に以下を実施します。

(1) NVIDIA Visual Profiler を開く
(2) File > Import をクリック
(3) “Select an import source:” で Nvprof を選択して Next をクリック
(4) “Multiple processes” を選択して Next をクリック
(5) “The nvprof pfofile files:” で Browse をクリックして読み込みたい *.nvp ファイルを選択
(6) Finish をクリック

実行すると、読み込まれた*.nvpファイルに保存されているプロライリングの結果がGUIで表示されます。

f:id:keimina:20181227221319p:plain
図2. NVIDIA Visual Profiler

GPUの状況などがタイムラインで表示されているのですが、縮小されていてみづらいため、虫眼鏡を押して拡大します。

f:id:keimina:20181227221733p:plain
図3. 虫眼鏡マークを押して “[0] Tesla…” の16秒 のあたりを拡大

凸凹してるところをさらに拡大してみます。

f:id:keimina:20181227221910p:plain
図4. 虫眼鏡マークを押して “[0] Tesla…” の16秒 のあたりをさらに拡大

このように拡大して詳細にみていくことができます。縮小すれば全体を俯瞰して見ることができます。左側のパネルから GPU をクリックし Context > Compute とクリックすると Compute の内訳が見れます(上記 図3. の真ん中の左側付近)。内訳の上から見ていくと、volta_gcgemm_64_32_nt と volta_gemm_64_32_tn というものが実行されているようです。それぞれ Compute の 19%、18%を占めていることがわかります。gemmはググるとどうやら行列の積和演算の処理らしいです。この二つで Compute の約40%を占めていることがわかります。その他、割愛しますが処理時間の平均(Avg. Duration)や、呼び出し回数(Invocations)など様々な情報が記載されていることがわかります。

7. プロファイリング結果の考察

tf.nn.conv2dのスクリプトが遅い件を調べるために、実際に使用したプロファイルのコマンドは以下になります。スクリプトでGPU2個を動かすために引数に2を指定しています。遅い件やスクリプトの詳細はTensorFlowで複数GPUで2次元畳み込みやってみる(tf.nn.conv2d with multiple GPUs) - keiminaの日記を参照してください。

nvprof -o profile_2gpu.nvp python3 multi_gpu.py 2 840 100 100

これで出力された *.nvp ファイルを NVIDIA Visual Profiler で開いて 左のパネルにある一つの GPUを選択して“Compute” をクリックして 右下にある Properties をみて見ると Kernels: 21.8ms となっていました。GPUでの計算時間は実質 21.8ms ということがわかります。左のパネルにあるもう一つの GPUを選択してみてもやはり 21msくらいでした。このことから、GPUで行われた計算の時間は最大でも 2 * 22ms = 44ms だと推測できます。Scipyを使用したCPUでの計算時間は約220msでしたので、少なくともCPUよりGPUのほうが約5倍は速いことになります。

f:id:keimina:20181227223808p:plain
図5. GPUの計算時間

余談:
私の環境では GPU 3 個までは表示できましたが、残念ながらGPU4個以上のものは アプリの制約なのか表示できませんでした。

8. まとめ・感想

前回記事でGPUに対する信頼度が低くなっていたのですが、今回の検証でGPUは速いと思えるようになったことが一番の収穫だと思います。GPUで何が行われているのか完全にブラックボックスだったのですが、このプロファイラによって完全なブラックボックスから少しは脱却できたのかなと思います。ただし当方、素人なため可視化しても理解できないところが多々ありますので、プロフェッショナル方からこのアプリの見方を教えていただきたいところです。なんだか、趣味の範囲を超えてきている感がありますので、ひとまずできることはやった、GPUは速かったということで、これからはあまり疑問をもたずGPUを使っていこうと思います。

以上です。おやすみなさい。