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を使っていこうと思います。

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

Pythonでドーナッツの形を描画する

いきなりですがPythonでドーナッツの形を描画します。まず、頭の中でその形をイメージします。次にそれを数式で表し、Python 3 で実装します。最後に描画されたものをみてみます。

1. 頭の中でドーナツの形をイメージする

ドーナッツの形をイメージします。以下のように円盤をz軸を中心にぐるぐる回転させると、その軌跡は、ドーナッツ形となるはずです。

f:id:keimina:20181218211055g:plain
回転する円盤

2. ドーナッツの形を数式で表す

ドーナッツの形を数式で表します。

以下のように、xz平面に円を定義します。媒介変数として、 u を用いると以下のように書けると思います。

f:id:keimina:20181218212123p:plain
xz平面の円の数式

※ u は角度、x1 と z1 は座標で定数だと思ってください。円盤の半径は1とします。

この円盤を z軸 を中心に回転させます。アファイン変換*1を使うと円盤を回転させることができます。具体的には、以下の回転行列*2を上記 x, y, z の数式に作用させることになります。今回はこれを使ってやってみましょう。

回転する角度を v とすると、回転行列は以下になります。

f:id:keimina:20181218211850p:plain
回転行列

この回転行列を 先ほど定義した xz平面の円上の点 x, y, z に作用させます。回転後の円上の点を x’, y’, z’ とすると、それらの数式は以下のようになります。

f:id:keimina:20181218212618j:plain
回転後の数式

無事に、回転後の数式を変数 u, v で表すことができました。
これは u, v の2次元ですので、平面(ドーナッツの表面)を表す関数だということがわかります。
ドーナッツの形をみるのには、表面の情報があれば十分なため、この数式で十分です。

回転後の数式が求まりましたので、次は実装を行います。

3. Pythonとmatplotlibで実装する。

2.で求めたドーナッツの形を表す数式をPythonで実装しmatplotlibで描画します。
ドーナッツの形は前述の通り、 u と v で表すことができるので、この u と v を少しづつ変化させて描画すればドーナッツの形になると思われます。
注意すべき点は、描画をおこなうグラフの軸のアスペクト比率を同じにすることです。
そうしないと、目盛は正しいのだけど見た目上わかりづらいグラフになってしまいます.

コード:

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np

N_us = 20
N_vs = 40

us = np.linspace(0, 2*np.pi, N_us)
vs = np.linspace(0, 2*np.pi, N_vs)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

x1 = 2
z1 = 0

for v in vs:
    x = np.cos(v)*(x1 + np.cos(us))
    y = np.sin(v)*(x1 + np.cos(us))
    z = z1 + np.sin(us)
    ax.scatter(x, y, z)

ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)
ax.set_zlim(-2, 2)
ax.set_aspect('equal')
plt.axis('off')

plt.show()

4. 描画されたものをみてみる

最後に描画されたものを見てみます。

f:id:keimina:20181218213219g:plain
ドーナッツ形の出力結果

カラフルなドーナッツが表示されました。
以上です。

*1:ここでは、点を別の点に移動させるくらいの意味で使ってます。

*2: https://ja.m.wikipedia.org/wiki/回転行列

TensorFlowで複数GPUで2次元畳み込みやってみる(tf.nn.conv2d with multiple GPUs)

1. やったこと

TensorFlow の tf.nn.conv2d関数 を GCE 上でK80とV100のGPUを1個〜8個を用いて、データ並列で実行。処理速度の検証を行った。

2. わかったこと

以下のコードを試した結果、CPUのほうが3倍早く、期待はずれの結果となってしまった。こうすれば早くなるなどのコメントがありましたらぜひご教授願います。

3. PCのスペック

CPU: virtual 12 core
MEM: 32 GB
GPU: K80x8個、V100x8個で検証

4. 検証に用いたコードの概要

multi_gpu.py

GPU の数と画像データ数、画像のサイズを引数にとり、データ並列でGPUで2次元畳み込みを行い、処理時間を出力するスクリプト

run_multi_gpu.py

multi_gpu.py を GPU の数を変えて実行し結果を記録するスクリプト

5. 実際のコード

TensorFlow の 複数 GPU の使い方は公式サイトの以下のページを参考にしました。
"https://www.tensorflow.org/guide/using_gpu#using_multiple_gpus"

multi_gpu.py
import numpy as np
import tensorflow as tf
import time
import sys

N_gpus = int(sys.argv[1])
N_images = int(sys.argv[2])
width = int(sys.argv[3])
height = int(sys.argv[4])

c = []

samples = np.ones((N_images,width,height,1), dtype=np.float32)
split_samples = tf.split(samples, N_gpus, axis=0)

for n in range(N_gpus):
    d = '/gpu:%d'%n
    with tf.device(d):
        fltr = tf.constant(np.ones((3,3,1,1), dtype=np.float32))
        conv = tf.nn.conv2d(split_samples[n], fltr, (1,1,1,1), 'SAME')
    c.append(conv)
    
with tf.device('/cpu:0'):
  sum = tf.add_n(c)
  
# Creates a session with log_device_placement set to True.
sess = tf.Session(config=tf.ConfigProto(log_device_placement=True))
# Runs the op.
s = time.time()
out = sess.run(sum)
e = time.time()
print(out)
print((e-s)*1000, 'ms')


out_text = '================================================\n'
out_text = \
    'N_gpus = {}, N_images = {}, width = {}, height = {}\n'\
    .format(N_gpus, N_images, width, height)
time_estimated = '%f ms\n'%((e-s)*1000)
out_text += time_estimated

filename = 'multi_gpu_time_estimated_{}_{}_{}_{}.txt'\
    .format(N_gpus, N_images, width, height)

with open(filename, 'w') as fp:
    fp.write(out_text)
run_multi_gpu.py
# -*- coding:utf-8 -*-
import subprocess

N_GPUS_MAX = 8
python_exe = 'python3'
script_name = 'multi_gpu.py'

# データ数は GPU の個数の倍数でなければならない
# 1, 2, 3, 4, 5, 6, 7, 8 の倍数は840
# 840個の画像データ、サイズ(100,100)の
# を2次元畳み込みすることを考える

n_images = 840
width = 100
height = 100

for n_gpus in range(1, N_GPUS_MAX+1):
    command = [
    	python_exe,
    	script_name,
    	str(n_gpus),
    	str(n_images),
    	str(width),
    	str(height)]
    out = subprocess.check_output(command)
    filename = 'multi_gpus_time_log_{}_{}_{}_{}.txt'\
        .format(n_gpus, n_images, width, height)
    with open(filename, 'wb') as fp:
        fp.write(out)

6. 実行結果

以下にV100 という GPU で実行した結果を載せます。
上から順に GPU1個、GPU2個、、、GPU8個と増やして実行た結果です。

================ GPU V100 ================
N_gpus = 1, N_images = 840, width = 100, height = 100
14897.803307 ms

N_gpus = 2, N_images = 840, width = 100, height = 100
4328.074455 ms

N_gpus = 3, N_images = 840, width = 100, height = 100
5645.370483 ms

N_gpus = 4, N_images = 840, width = 100, height = 100
684.230804 ms

N_gpus = 5, N_images = 840, width = 100, height = 100
744.736195 ms

N_gpus = 6, N_images = 840, width = 100, height = 100
731.468678 ms

N_gpus = 7, N_images = 840, width = 100, height = 100
726.519346 ms

N_gpus = 8, N_images = 840, width = 100, height = 100
721.646070 ms

同様に、K80 という GPU で実行した結果を以下に載せます。

================ GPU K80 ================
N_gpus = 1, N_images = 840, width = 100, height = 100
3248.390198 ms

N_gpus = 2, N_images = 840, width = 100, height = 100
3380.734682 ms

N_gpus = 3, N_images = 840, width = 100, height = 100
3837.318182 ms

N_gpus = 4, N_images = 840, width = 100, height = 100
690.669060 ms

N_gpus = 5, N_images = 840, width = 100, height = 100
691.730976 ms

N_gpus = 6, N_images = 840, width = 100, height = 100
588.120937 ms

N_gpus = 7, N_images = 840, width = 100, height = 100
581.115961 ms

N_gpus = 8, N_images = 840, width = 100, height = 100
647.247314 ms

7. GPU がなんとなく遅い気がするのでCPUでやってみた結果

import numpy as np
import scipy
from timeit import timeit

fltr = np.ones((3,3), dtype=np.float32)
a = np.ones((100,100), dtype=np.float32)

N_images = 840

と定義して以下の scipy.signal.convolve2d で畳み込みをやってみると、、、

>>> %timeit [scipy.signal.convolve2d(a, fltr, 'same') for _ in range(N_images)]
222 ms ± 4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

なぜかCPU1個の方がGPU8個より3倍くらい早いです。

8. GPU の使用数が 1-3 のときに遅い原因について推測してみる

GPU が 2-3 個のとき、処理時間が 1000 ms を超えており極端に遅いです。
極端に遅いケースについては、データが GPU のメモリに入りきらず、
外部とデータのやりとりを行うスワップのようなものが発生していることが考えられます。

9. CPU より GPU が遅い原因について推測、検証方法を検討してみる(検証はしません)

・TensorFlow がたいして最適されていない
⇨ Keras でも同様の結果になるか確認
⇨ Chainer など他のライブラリを使用してみる
・畳み込みのフィルタサイズの影響
⇨ フィルタサイズ 3x3 でなく10x10 などサイズを大きくする
・チャンネルのサイズの影響
⇨ チャンネルのサイズを 1 でなく 10 などサイズを大きくする
GPUからCPUに計算結果を転送するのがボトルネックになっている
⇨ データ転送するだけのコードで検証する
GPU内でずっと計算し続けて、オーバーヘッドの影響をなくすようなコードで検証する

10. まとめ

勝手な想像ですが、 TensorFlow は tf.nn.conv2d のためにあるわけではなく、どちらかというと、ニューラルネットの重み計算を行うために最適化されているはずです。そのため、静的な数式で表されるような重み更新の計算式などはうまくGPUで分散して計算できると思うのですが、畳み込みのようなシーケンシャルな処理は単純な計算式で表現できないため、それも遅さの原因の一つかなと思います。ただ、そうはいっても今回の例ではデータ並列したはずなので、GPUの数に比例して早くなるはずだったのですが、そうはなりませんでした。現実は厳しかったです。このような畳み込みだけを行うという特殊なケースでは単純にGPUを増やせば早くなるというわけではないようですので注意が必要かと思われます。
いろいろ書きましたが、検証は面倒なので、TensorFlowの開発者が最適化をすすめて、コードを書き換えなくても自動で計算リソースをスケールしてくれる日がいつかくることを期待して、メインの機械学習アルゴリズムについて学習をすすめて行こうと思います。

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