「Alt-Tab」の動作を真面目に解析し、同じような動作をするリストの操作アルゴリズムを実装する

1. はじめに

WindowsMac でアプリを切り替える時などに使うショートカット「Alt-Tab」 のような動作を行うリストの操作アルゴリズムを、Python で実装します。仮想的な 「Alt-Tab」 を実装しその動作の理解を深めます。

2. 「Alt-Tab」 とは

「Alt-Tab」 とはキーの操作のことです。Altキーを押しながらTabキーを押すことを表します。これによりアプリケーションを切り替えることができます。

3. 「Alt-Tab」 を動作させて動きをみる

「Alt-Tab」 を押下した時の内部の動きをみてみます。ここでは、 「Alt-Tab」 を押した時画面にに現れるリストをウィンドウリストと呼ぶこととします。 基本的には、「Alt-Tab」 を押した後 Altキーを離すとウィンドウが切り替わります。以下にウィンドウを3回切り替えた場合を例として、Alt キーと Tab キーとリスト(ウィンドウリスト)の状態がどうなるかみてみます。
(
1回目の切り替え:時刻1-時刻5、リストの状態1-リストの状態3
2回目の切り替え:時刻6-時刻12、リストの状態4-リストの状態7
3回目の切り替え:時刻13-時刻21、リストの状態8-リストの状態12
)

※Alt キーと Tab キーは「押された」「離された」「何もされていない」の3つの状態をとることとします。
※以下の例では「何もされていない」は「-」として表しています。

時刻 Altキー Tabキー リストの状態 リストの状態(画像)
1 - - 1 f:id:keimina:20190410001718p:plain
2 押された - 1 f:id:keimina:20190410001718p:plain
3 - 押された 2 f:id:keimina:20190410001948p:plain
4 - 離された 2 f:id:keimina:20190410001948p:plain
5 離された - 3 f:id:keimina:20190410002034p:plain
6 - - 4 f:id:keimina:20190410002052p:plain
7 押された - 4 f:id:keimina:20190410002052p:plain
8 - 押された 5 f:id:keimina:20190410002202p:plain
9 - 離された 5 f:id:keimina:20190410002202p:plain
10 - 押された 6 f:id:keimina:20190410002232p:plain
11 - 離された 6 f:id:keimina:20190410002232p:plain
12 離された - 7 f:id:keimina:20190410002256p:plain
13 - - 8 f:id:keimina:20190410002318p:plain
14 押された - 8 f:id:keimina:20190410002318p:plain
15 - 押された 9 f:id:keimina:20190410002343p:plain
16 - 離された 9 f:id:keimina:20190410002343p:plain
17 - 押された 10 f:id:keimina:20190410002406p:plain
18 - 離された 10 f:id:keimina:20190410002406p:plain
19 - 押された 11 f:id:keimina:20190410002425p:plain
20 - 離された 11 f:id:keimina:20190410002425p:plain
21 離された - 12 f:id:keimina:20190410002445p:plain

4. 「Alt-Tab 」を動作させてわかったこと

以下のことがわかりました。
※ウィンドウリストの何番目が選択されているかを表すインデックスをウィンドウインデックスと呼ぶこととします。
※インデックスは 0 始まりとします。

No. Alt Tab 何が行われるか
1 押されている状態 離された ウィンドウインデックスを右に一つずらす
2 離された 離された 選択しているアイコンをウィンドウリストの先頭に持ってっくる

Alt-Tab でウィンドウインデックスを動かし、ウィンドウリストを選択し選択したウィンドウをウィンドウリストの先頭に移動するという動作になっています。
これはすなわち、 lst[i] の要素を lst の先頭に挿入して移動する問題ということになります。

これは、 Python では以下の1行のコードで実現できます。(参考文献参照)

lst.insert(0, lst.pop(i))

また、Altキーが押しっぱなしだと、ウィンドウが切り替わらないので、Altキーが離されたことの検知が必要なことがわかります。また、 Tabキーの検知も同様に必要なことがわかります。これらは個人的にイベントハンドラ(on_mouse_click関数などのonナンチャラ関数)っぽい名前で定義した方がわかりやすいと思いますのでそうします。以上をもとに実装します。

5. 実装したコード

実装したコードを以下に記載します。

def init_environment():
    global frame_lst, index
    index = 0
    frame_lst = ['Emacs', 'Finder', 'Xcode', 'MS Word']

def on_alt_tab_pressed():
    global frame_lst, index
    index += 1
    index = index%len(frame_lst)

def on_alt_tab_released():
        global frame_lst, index
        frame_lst.insert(0, frame_lst.pop(index))
        index = 0

def hit_alt_tab_n_times(n):
    global frame_lst
    print('==== Before hitting Alt-TAB ====')
    print(frame_lst)
    for _ in range(n):
        on_alt_tab_pressed()
    on_alt_tab_released()
    print('==== After hitting alt-TAB consecutively %d times ===='%n)
    print(frame_lst)

# test 1
init_environment()
hit_alt_tab_n_times(1)
hit_alt_tab_n_times(2)
hit_alt_tab_n_times(3)

6. 出力結果

出力結果を記載します。

==== Before hitting Alt-TAB ====
['Emacs', 'Finder', 'Xcode', 'MS Word']
==== After hitting alt-TAB consecutively 1 times ====
['Finder', 'Emacs', 'Xcode', 'MS Word']
==== Before hitting Alt-TAB ====
['Finder', 'Emacs', 'Xcode', 'MS Word']
==== After hitting alt-TAB consecutively 2 times ====
['Xcode', 'Finder', 'Emacs', 'MS Word']
==== Before hitting Alt-TAB ====
['Xcode', 'Finder', 'Emacs', 'MS Word']
==== After hitting alt-TAB consecutively 3 times ====
['MS Word', 'Xcode', 'Finder', 'Emacs']

7. コードの説明

init_environment 関数でウィンドウリストを初期化しています。on_alt_tab_pressed 関数で Altキー押した状態で Tabキーが1回押された時の処理を行います。Tabキーをたくさん打ったとしても、インデックスがウィンドウリストのサイズを超えないよう「%」演算子を使用してサイクルさせています。on_alt_tab_released 関数で Altキーを離した時の処理を行います。ウィンドウリストのウィンドウインデックスの要素をウィンドウリストの先頭に挿入し移動させるだけです。hit_alt_tab_n_times 関数で Altキー押した状態で Tabキーがn回押され、Altキーが離された時の処理を行います。また処理前と処理後のリストの状態を表示します。結果をみると、例で示した通りの順番になっており正しく動作していることがわかりました。

8. まとめ、考察

「Alt-Tab」の動作を真面目に解析し、同じような動作をするリストの操作アルゴリズムを実装しました。「Alt-Tab」は1日に数えきれないくらい使用しており、個人的に重要だと思う機能なのですが、機能を担うコアとなるアルゴリズムは、非常ににシンプルであることがわかりました。アルゴリズムはシンプルだからと行って重要でないということはないという教訓が得られました。

9. 参考文献

insertとpopメソッドを使用してリストの先頭への挿入と移動を行う方法についてはこちらを参考にしました。
stackoverflow.com

以上。おやすみなさい。

sympy でニューラルネットワークの重みの更新に使用する式を計算グラフにしてみる

1. はじめに

sympy(https://www.sympy.org/en/index.html) を使用し、ニューラルネットワークの重みの更新に使用する式を計算グラフにしてみます。具体的には、 sympy を使って、損失関数の勾配の計算を計算グラフとして表し、それを graphviz で画像に出力します。sympy と graphviz を使用すると、 数式を簡単に計算グラフに変換できます。ちなみに計算グラフとは、以下のようなものです。(実際にこれらを使用し作成したものです)

f:id:keimina:20190325203008p:plain
図1. c = a + b の計算グラフ

このように計算式をグラフとして表現したものを計算グラフと呼びます。

2. 使用するモジュールと、そのインストール

本記事では、以下の二つのモジュールを使用します。グラフの描画に graphviz というモジュールを使用します。sympy は計算グラフをテキストで出力できますが、単体では画像で出力することができないためグラフ描画アプリ(graphviz)を使用します。

  1. sympy
  2. graphviz

graphviz は以下のコマンドでインストールしました。

pip install graphviz

コマンドの出力結果

Collecting graphviz
  Downloading https://files.pythonhosted.org/packages/1f/e2/ef2581b5b86625657afd32030f90cf2717456c1d2b711ba074bf007c0f1a/graphviz-0.10.1-py2.py3-none-any.whl
distributed 1.21.8 requires msgpack, which is not installed.
Installing collected packages: graphviz
Successfully installed graphviz-0.10.1
You are using pip version 10.0.1, however version 19.0.3 is available.
You should consider upgrading via the 'pip install --upgrade pip' command.

コマンド実行後、"requires msgpack" と出力されていますが、 最終的に"Successfully installed" と出力されているため、問題なしとみなし気にせず先に進みました。

今回は、計算グラフを画像で出力したいため、 さらに以下の手順行う必要がありました。

  1. MacPorts のインストール(インストール手順は以下のサイトを参照)
  2. 以下のコマンドの実行しアプリをインストール
sudo port install graphviz

sympy のインストール手順について説明は省略します。

3. ニューラルネットワークのモデルと、順伝播の式

本記事で扱うニューラルネットワークのモデルと、順伝播の式は以下の通りです。

f:id:keimina:20190325203628p:plain
図2. ニューラルネットワークのモデルと、順伝播の式

4. ソースコード

3. のニューラルネットワークの損失関数の勾配を計算するコードにすると以下のようになります。
※本記事は、重みの更新を行う計算を、計算グラフに出力することを目的としているため、説明をわかりやすくするため損失関数 l は o1 としています。本来であれば、l = (o1 - y)**2 のような定義をすべきですが、パラメータが増えると、計算グラフが大きくなり見づらくなるため、簡略化しています。ご了承ください。
ソースコードで numpy を使用して行列計算していますが、sympy での行列計算の仕方がわからなかったためです。あまり真似しないように。。。

from sympy import symbols, Function, tanh
from sympy.printing import dot
import graphviz
import numpy as np

# sympy で使用する変数(シンボル)を定義する
# 入力値、中間層の値
i1, i2 = symbols("i1, i2")
h11, h12 = symbols("h11, h12")

# 重み
w111, w112, w121, w122 = symbols("w111, w112, w121, w122")
w211, w212 = symbols("w211, w212")

# バイアス
b11, b12 = symbols("b11, b12")
b21 = symbols("b21")

# numpy を利用して行列形式にする
I = np.array([[i1], [i2]])
W1 = np.array([[w111, w112], [w121, w122]])
W2 = np.array([[w211, w212]])
B1 = np.array([[b11], [b12]])
B2 = np.array([[b21]])

# 活性化関数
TANH = np.vectorize(tanh)

# 順伝播の数式を作成する
H1 = TANH(W1.dot(I) + B1)
O = TANH(W2.dot(H1) + B2)

# 行列からスカラを取り出す
o1 = O[0, 0]

# 最急降下法 によって重み w1 を更新する式(学習係数は 1 とする)
# 損失関数を定義する、重み更新を行う計算グラフを出したいだけなので l = o1  としてしまう
l = o1
dl_dw111 = l.diff(w111)

expr_forward = o1
src_forward = graphviz.Source(dot.dotprint(expr_forward))

expr_backward = dl_dw111
src_backward = graphviz.Source(dot.dotprint(expr_backward))

# 私の環境ではExecutableNotFound エラーが発生したので、以下のコードで回避しました
app_path = '/opt/local/bin/dot'
graphviz.ENGINES.add(app_path)
src_forward._engine = app_path
src_backward._engine = app_path

# 画像ファイルを作成する
filename_forward = "calc_forward_graph"
src_forward.render(filename_forward)
filename_backward = "calc_backward_graph"
src_backward.render(filename_backward)

順伝播の式が o1 です。勾配降下法によって損失関数 l を最小にするため重みパラメータ w111 を更新するときに行われる計算が dl_dw111 です。コードの最後の方で、sympy と graphviz を使用して、o1 と dl_dw111 の計算を行う計算グラフを出力しています。sympy.printing.dot で計算グラフを graphviz で使用可能な dot フォーマット形式のテキストに変換しています。そして、 dot フォーマット形式のテキストを graphviz.Source で graphviz に取り込み、 render メソッドを呼びだし pdf ファイルを作成しています。

5. 出力結果

出力結果は以下になります。図は上から順に、順伝播の計算を行うための計算グラフ(o1 の計算を行うグラフ)、損失関数lを最小化すべく w111 を勾配降下法により更新するための計算グラフ(dl_dw111 の計算を行うグラフ)、となっています

f:id:keimina:20190325204132p:plain
図3. 順伝播の計算グラフ("calc_forward_graph")

f:id:keimina:20190325204331p:plain
図4. w111 の重みを更新する数式の計算グラフ("calc_backward_graph")

楕円の中に書かれている Mul 、 Add 、 Pow 、 tanh はどのような計算を行うかを表していて、それぞれ、 乗算、加算、累乗、ハイパボリックタンジェントを表しています。ニューラルネットの学習時には、図中の w112 、w121、 b11 、、、 などの変数のように見えるものは全て値が割り振られているため、 この計算グラフの計算を行うと、何かしらの値(スカラ値)が求まります。そしてこの値を使用して重み w111 を更新します。

6. まとめ

sympy と graphviz を使用して、ニューラルネットワークの重みの更新に使用する式を計算グラフにしました。このように、ニューラルネットワークのモデルが決まると、重みを更新するための式、すなわち計算グラフも決まります。

7. 参考文献

  1. ゼロから作るDeep Learning――Pythonで学ぶディープラーニングの理論と実装

Numpy 入門 (機械学習勉強会 in 新潟 2019/03/23)

Numpy 入門

Numpy 入門
(機械学習勉強会 in 新潟 2019/03/23)

1 はじめに

本記事は、機械学習勉強会 in 新潟 2019/03/23 の発表資料です。 numpy の入門者向けの記事になります。 numpy とは何か、何ができるかなどを説明します。numpy で簡単な数値計算を行えるようになることを目的とします。

※時間の都合上、 numpy の機能のほんの一部を紹介しています。 本記事は以下の参考文献を参考に書いています。 numpy についての詳細は以下の参考文献を参照してください。

2 numpy とは何か

Python数値計算のためのライブラリです。 numpy によってできることとメリットは以下になります

No. ライブラリ 説明 できること、メリット
1 numpy 数値計算ライブラリ 1. 行列の計算を行うコードを短く書ける
      2. 行列の計算、ベクトルの計算の高速化( 速度比較 参照)

3 numpy モジュールのインポート方法について

numpy を使用する場合、numpy モジュールのインポートは以下のように行います。

コード1. 「モジュールのインポート」:

import numpy as np

これにより、'np' の2文字で numpy を使用することができます。

コード2. 「モジュールの使用方法」:

import numpy as np
arr = np.array([1,2,3])
arr
Out: array([1, 2, 3])

使用する名前は numpy など np 以外の名前にもできますが、 np と書くのが一般的になっているためこちらを使用しましょう。

4 数値を作成する

はじめに、numpy で数値を作成しましょう。数値は、スカラ(Scalar)と呼ばれることがあります。 これらは、後述する配列の要素( 5.1.1 参照)となりますので、しっかり覚えましょう。 np.int32 を使用すると numpy の整数を作成できます。

コード3. 「numpy で数値の作成」:

import numpy as np
answer = np.int32(100)
answer
Out: 100

4.1 np.int32 とは

np.int32 はデータタイプ(dtype)として定義されています。 np.int32(100) を実行すると np.int32 データタイプのオブジェクトが作成されます。 np.int32 以外にも数値を表わすデータタイプが用意されています。

データタイプ 数値の種類 用途
np.uint8 8bit 符号無し整数 グレースケール画像
np.int32 32bit 整数 画像処理
np.float32 32bit 浮動少数 画像処理、平均値、確率
np.float64 64bit 浮動少数 数値計算(精度が求められるもの)
   

真偽値(True/False)を現す場合は np.bool データタイプが用意されています。

5 「1次元配列」の数値計算をする

配列の基本となる1次元配について学習しましょう。

5.1 「1次元配列」を作成する

1次元の配列を作成しましょう。1次元配列とは以下のようなものです。配列は'要素'と呼ばれるもので構成されています。

値1 値2 値3 ,,,

5.1.1 配列の'要素'とは

一般的に、上記 5.1 のような配列の中にある値1や値2などのことを配列の”要素”(Element)と呼びます。 1次元配列だけでなく、2次元配列や3次元配列も同様です。配列は要素によって構成されています。

※ numpy では配列を作成する関数の引数に dtype を指定することで、その配列の要素のデータタイプを指定できます。 引数 dtype を指定しない場合は、自動的に配列のデータタイプが決定されます。 numpy では配列は同一のデータタイプの要素から構成されます。

5.1.2 np.array 関数

np.array 関数を使用すると、標準 Python のリストを配列に変換することができます.

コード4. 「標準 Python のリストを配列に変換する」:

import numpy as np
arr = np.array([1, 2, 3])
arr
Out: array([1, 2, 3])

5.1.3 np.arange 関数

np.arange 関数を使用すると、数値の1次元配列を作成することができます。

コード5. 「np.arange 関数を使用し、数値の1次元配列を作成」:

import numpy as np
arr = np.arange(3)
arr
Out: array([0, 1, 2])

5.2 「1次元配列」の形状 shape について

配列には shape と呼ばれる、配列の形状を表す情報があります。 配列の形状 shape を調べるには以下のようにします。 この例の場合、配列は長さ3の1次元配列であることがわかります。

コード6. 「配列の形状 shape を調べる」:

import numpy as np
arr = np.arange(3)
answer = arr.shape
answer
Out: (3,)

5.3 「1次元配列」 + 「1次元配列」

1次元配列と1次元配列の四則演算の基本を学習します。 ここでは例として、1次元配列と1次元配列の足し算を行います。 以下のように、1次元配列を二つ作成し、足し算記号(+) を使用して求めます。 計算は、要素ごと(element wise)に行われます。

コード7. 「1次元配列と1次元配列の足し算」:

import numpy as np
arr_1 = np.array([1, 2, 3])
arr_2 = np.array([0, 0, 1])
answer = arr_1 + arr_2
answer
Out: array([1, 2, 4])

その他の演算子(-*/など)についても、同じように、配列の要素ごとに計算されます。 形状の異なる配列同士の計算は、エラーになるか、ならない場合があります。 形状の異なる配列同士の計算については 11 を参照してください。

5.4 配列の要素に数学関数を適用する

1次元配列の要素ごと(element wise)に数学関数を適用してみましょう。 ここでは例として平方根を適用してみましょう。

コード8. 「配列の各要素に平方根を適用」:

import numpy as np
arr = np.array([1, 2, 3])
answer = np.sqrt(arr)
answer
Out: array([1.        , 1.41421356, 1.73205081])

np.sqrt 関数以外にも様々な数学関数が用意されています。 それらの関数も、上記のように使用することで、要素ごとに適用することができます。 以下に数学関数の一部を記載します。

No. 関数 説明
1 np.exp 指数関数
2 np.log 対数関数
3 np.sin sin関数
4 np.arcsin arcsin関数

5.5 1次元配列の要素の合計を求める

合計を求める場合は、 np.sum 関数を使用します。

コード9. 「np.sum 関数を使用して1次元配列の要素の合計を求める」:

import numpy as np
arr = np.arange(11)
answer = np.sum(arr)
answer
Out: 55

6 「2次元配列」の数値計算をする

2次元配列について学習しましょう。

6.1 「2次元配列」を作成する

2次元の配列を作成しましょう。2次元配列とは以下のようなものです。

値11 値12 値13 ,,,  
値21 値22 値23 ,,,  
値31 値32 値33 ,,,  
,,,        

6.1.1 np.array 関数

np.array 関数を使用すると、2次元配列を作成することができます。(1次元配列の時と同様)

コード10. 「np.array 関数を使用して2次元配列を作成する」:

import numpy as np
arr = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
arr
Out: 
array([[1, 2, 3],
       [4, 5, 6],
       [7, 8, 9]])

6.1.2 reshape メソッド

配列の reshape メソッドを使用すると配列の形状を変えることができます。 ここでは、 1次元配列の形状を変えて、2次元配列を作成してみましょう。

コード11. 「reshape メソッドを使用して1次元配列の形状を変えて、2次元配列を作成する」:

import numpy as np
arr = np.arange(9)
answer = arr.reshape((3, 3))
answer
Out: 
array([[0, 1, 2],
       [3, 4, 5],
       [6, 7, 8]])

6.2 「2次元配列」 + 「2次元配列」

2次元配列同士の四則演算、関数を適用した時の動作は、1次元配列の時と同様、要素ごとに計算が行われます。 (そのため2次元配列の四則演算については説明を省略します。)2次元配列を使用すると、行列の掛け算を行うことができます。

行列の掛け算とは、具体的には以下のような計算のことです。 行列 A と行列 B の掛け算(積)は A・B と書き、A と B と A・B は以下のような関係になります。 ※説明の簡略化のため行列の大きさは2x2としています。

  • 行列 A

    a b
    c d
  • 行列 B

    e f
    g h
  • 行列 A・B

    a*e+b*g a*f+b*h
    c*e+d*g c*f+d*h

行列の掛け算を行うには np.dot 関数を使用します。

コード12. 「np.dot 関数を使用して行列の掛け算を行う」:

import numpy as np
arr_1 = np.arange(9).reshape((3, 3))
arr_2 = np.arange(9).reshape((3, 3))
answer = np.dot(arr_1, arr_2)
answer
Out: 
array([[ 15,  18,  21],
       [ 42,  54,  66],
       [ 69,  90, 111]])

配列オブジェクトの dot メソッドを使用して以下のように書くこともできます。

コード13. 「dot メソッドを使用して行列の掛け算を行う」:

import numpy as np
arr_1 = np.arange(9).reshape((3, 3))
arr_2 = np.arange(9).reshape((3, 3))
answer = arr_1.dot(arr_2)
answer
Out: 
array([[ 15,  18,  21],
       [ 42,  54,  66],
       [ 69,  90, 111]])

注意: np.dot は '∗'による掛け算と同じではありません。 '∗' による掛け算は要素ごと(element wise)に掛け算が行わます。 一方 np.dot は行列の掛け算ですので、積和演算(掛け算と足し算)が行われます。

6.3 「2次元配列」の要素の合計を求める

'すべての要素の合計'を求める場合は、1次元配列の時と変わりませんので省略しますが、 2次元配列では、さらに、'列の合計'と'行の合計'を求めることができます。

6.3.1 列の合計

np.sum 関数の引数で axis=0 とすることで列の合計を求めることができます。

コード14. 「列の合計を求める」:

import numpy as np
arr = np.arange(9).reshape((3, 3))
answer = np.sum(arr, axis=0)
answer
Out: array([ 9, 12, 15])

6.3.2 行の合計

np.sum 関数の引数で axis=1 とすることで行の合計を求めることができます。

コード15. 「行の合計を求める」:

import numpy as np
arr = np.arange(9).reshape((3, 3))
answer = np.sum(arr, axis=1)
answer
Out: array([ 3, 12, 21])

6.3.3 配列の次元数を保つには

6.3.16.3.2 で np.sum を使用して求めた配列を見ると、2次元配列の入力に対して出力が1次元配列となっており、次元数が変わってしまうことに気付くかと思います。 np.sum 関数の引数で keepdims=True とすることで、配列の次元数を保つことができます。

コード16. 「np.sum 関数の引数で keepdims=True とし配列の構造をなるべく保つ」:

import numpy as np
arr = np.arange(9).reshape((3, 3))
answer = np.sum(arr, axis=1, keepdims=True)
answer
Out: 
array([[ 3],
       [12],
       [21]])

6.4 行列計算用の関数

2次元配列には行列計算用の関数を使用することができます。 行列計算用の関数は np.linalg にあります。 その一部を紹介します。

関数 説明
np.linalg.inv 逆行列を求めるための関数
np.linalg.norm 行列のノルムを求めるための関数

7 配列を作成する関数

配列を作成する方法は、5.16.1 で紹介した方法が全てではありません。 以下に、配列を作成するための関数の一部を紹介します。

7.1 固定された値を要素にもつ配列を作成する関数

No. 関数 説明
1 np.zeros 全ての要素が 0 の配列を作成します
2 np.ones 全ての要素が 1 の配列を作成します
3 np.full 全ての要素が x の配列を作成します(x は引数で指定)

コード17. 「np.zeros 関数を使用し全ての要素が 0 の配列を作成」:

import numpy as np
arr = np.zeros((3, 3))
arr
Out: 
array([[0., 0., 0.],
       [0., 0., 0.],
       [0., 0., 0.]])

7.2 単位行列を作成する関数

No. 関数 説明
1 np.eye 単位行列を表す配列を作成します

コード18. 「np.eye 関数を使用し単位行列を表す配列を作成」:

import numpy as np
arr = np.eye(3)
arr
Out:
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])

7.3 ランダムな値を要素にもつ配列を作成する関数

No. 関数 説明
1 np.random.randint ランダムな整数の配列を作成します
2 np.random.random ランダムな少数の配列を作成します

コード19. 「np.random.randint 関数を使用しランダムな整数の配列を作成」:

import numpy as np
arr = np.random.randint(0, 10, (3, 3))
arr
Out:
array([[4, 0, 7],
       [2, 3, 4],
       [3, 0, 5]])

コード20. 「np.random.random 関数を使用しランダムな少数の配列を作成」:

import numpy as np
arr = np.random.random((3, 3))
arr
Out:
array([[0.23963676, 0.14041946, 0.67624721],
       [0.75209081, 0.3182057 , 0.18579186],
       [0.0888395 , 0.87384273, 0.24740674]])

7.4 似たような配列を作成する関数

配列の計算をプログラミングしているときに、ある配列と同じ形状をもつ配列だが、値は違うという配列を作成したい時があります。そのような場合は、以下の関数を使用します。

No. 関数 説明
1 np.zeros_like 引数で指定した配列と同じ形状
   
    全ての要素が 0 の配列を作成します
     
2 np.ones_like 引数で指定した配列と同じ形状
   
    全ての要素が 1 の配列を作成します

コード21. 「np.zeros_like 関数を使用し形状が同じ配列を作成」:

import numpy as np
arr = np.random.randint(0, 10, (3, 3))
arr_2 = np.zeros_like(arr)
arr_2
Out[4]: 
array([[0, 0, 0],
       [0, 0, 0],
       [0, 0, 0]])

7.5 ある範囲の数値を要素にもつ配列を作成する関数

No. 関数 説明
1 np.linspace 数直線の範囲を等幅で分割した1次元配列を作成します
    1変数関数を調べる時などに使用できます
2 np.meshgrid 2次元座標などの座標を配列を使い表現するときに使用します。
    2変数関数を調べる時などに使用できます

コード22. 「np.linspace 関数の使用例」:

import numpy as np
arr = np.linspace(0.0, 0.9, 10)
arr
Out: array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9])

コード23. 「np.meshgrid 関数の使用例」:

import numpy as np
x = np.linspace(0, 0.4, 5)
y = np.linspace(0, 0.4, 5)
x_mesh, y_mesh = np.meshgrid(x, y)
answer = x_mesh, y_mesh
answer
Out:
(array([[0. , 0.1, 0.2, 0.3, 0.4],
        [0. , 0.1, 0.2, 0.3, 0.4],
        [0. , 0.1, 0.2, 0.3, 0.4],
        [0. , 0.1, 0.2, 0.3, 0.4],
        [0. , 0.1, 0.2, 0.3, 0.4]]), array([[0. , 0. , 0. , 0. , 0. ],
        [0.1, 0.1, 0.1, 0.1, 0.1],
        [0.2, 0.2, 0.2, 0.2, 0.2],
        [0.3, 0.3, 0.3, 0.3, 0.3],
        [0.4, 0.4, 0.4, 0.4, 0.4]]))

8 配列の要素へのアクセス

8.1 配列の要素の取得

添字(インデックス)を使用することで、配列の要素を取得することができます。

8.1.1 要素の取得の基本

以下のようにすると、配列の要素の値を取得することができます。

コード24. 「1次元配列の配列の要素の取得の基本」:

import numpy as np
arr = np.array([1, 3, 5])
arr[2]
Out: 5

コード25. 「2次元配列の配列の要素の取得の基本(1)」:

import numpy as np
arr = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]])
arr[2][2]
Out: 23

上記コードは以下のようにも書くことができます。

コード26. 「2次元配列の配列の要素の取得の基本(2)」:

import numpy as np
arr = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]])
arr[2, 2]
Out: 23

8.1.2 要素を範囲で指定して取得

以下のようにすると、指定した範囲にある配列の要素を取得することができます。

コード27. 「1次元配列の要素を範囲で指定して取得」:

import numpy as np
arr = np.array([1, 3, 5])
arr[1:]
Out: array([3, 5])

コード28. 「2次元配列の要素を範囲で指定して取得」:

import numpy as np
arr = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]])
arr[1:, 1:]
Out: 
array([[11, 13],
       [21, 23]])

8.1.3 逆順の配列の取得

以下のようにすると、逆順の配列を取得できます。

コード29. 「逆順の配列の取得」:

import numpy as np
arr = np.arange(10)
arr[::-1]
Out:
array([9, 8, 7, 6, 5, 4, 3, 2, 1, 0])

8.1.28.1.3 で取得した arr[1:, 1:] や arr[::-1] は、あたかも配列のコピーのように思われるかもしれませんが、実際は配列のコピーではありません。これらはビューと呼ばれるものであり、実態は arr のデータを参照しているものです。扱いには注意しましょう。 8.2 でも説明します。

8.2 配列の要素の書き換え

8.2.1 要素の書き換えの基本

配列の要素を添字(インデックス)で指定して別の値に書き換えることができます。 以下のようにすると、配列の要素の値を書き換えることができます。

コード30. 「1次元配列の配列の要素の書き換えの基本」:

import numpy as np
arr = np.array([1, 3, 5])
arr[2] = 100
arr
Out: array([  1,   3, 100])

コード31. 「2次元配列の配列の要素の書き換えの基本(1)」:

import numpy as np
arr = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]])
arr[2][2] = 100
arr
Out:
array([[  1,   3,   5],
    [  7,  11,  13],
    [ 17,  21, 100]])

上記コードは以下のようにも書くことができます。

コード32. 「2次元配列の配列の要素の書き換えの基本(2)」:

import numpy as np
arr = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]])
arr[2, 2] = 100
arr
Out:
array([[  1,   3,   5],
       [  7,  11,  13],
       [ 17,  21, 100]])

8.2.2 要素を範囲で指定して書き換え

以下のようにすると、指定した範囲にある配列の要素を書き換えることができます。

コード33. 「1次元配列の要素を範囲で指定して書き換え」:

import numpy as np
arr = np.array([1, 3, 5])
arr[1:] = 100
arr
Out: array([  1, 100, 100])

コード34. 「2次元配列の要素を範囲で指定して書き換え」:

import numpy as np
arr = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]])
arr[1:, 1:] = 0
arr
Out: 
array([[ 1,  3,  5],
       [ 7,  0,  0],
       [17,  0,  0]])

上記コードの arr[1:] や arr[1:, 1:] は配列のコピーではないことに注意しましょう。 このようなものを numpy ではビュー(View)と呼びます。実態は arr のデータを参照しているものです。 このような ビュー arr[1:] や arr[1:, 1:] に対する値の代入は、ビューの参照先である arr に対する値の代入であることを意味します。

よくあるのが、値を保持しておきたい配列のビューに代入を行なってしまいその配列を破壊してしまうことです。 値を代入したい場合は、以下のように、配列のコピーを作成して、そのコピーに対して代入しましょう。

コード35. 「配列のコピーを作成して、そのコピーに対して代入する」:

import numpy as np
arr_orig = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]])
arr = arr_orig.copy()
arr[1:, 1:] = 0
arr_orig
Out: 
array([[ 1,  3,  5],
       [ 7, 11, 13],
       [17, 21, 23]])

9 真偽値(True/False)について

'==' や '<' などの比較演算子を使用すると配列の要素ごとに比較が行われます。

コード36. 「比較演算子による配列の比較」:

import numpy as np
arr_1 = np.arange(9).reshape((3, 3))**2
answer = arr_1 >= 10
answer
Out: 
array([[False, False, False],
       [False,  True,  True],
       [ True,  True,  True]])

この真偽値を配列のインデックスに指定すると、 真偽値のTrue の箇所の要素だけ集めた1次元配列を求めることができます。

コード37. 「真偽値の配列による要素の収集」:

import numpy as np
arr_1 = np.arange(9).reshape((3, 3))**2
answer = arr_1[arr_1 >= 10]
answer
Out: array([16, 25, 36, 49, 64])

以下のように、 np.where を使用すると値が True のインデックスを求めることができます。 np.where の戻り値は、「インデックスを要素にもつ配列」を要素にもつタプル(tuple)です。

コード38. 「np.where による条件を満たす配列のインデックスの取得」:

import numpy as np
arr_1 = np.arange(9).reshape((3, 3))**2
answer = np.where(arr_1 >= 10)
answer
Out: (array([1, 1, 2, 2, 2]), array([1, 2, 0, 1, 2]))

10 配列の制約

ここまで1次元配列、2次元配列の計算方法について説明しました。 numpy の配列は、標準Python のリストと違い、一度作成した配列に要素を追加したり、削除したりすることはできません。 一度作成した配列のデータタイプ(dtype)を変更することもできません。これらの操作を行いたい場合は、新たに配列を作成してください。 これらの制約を設けることにより高速(C言語並)な計算が可能になっています。

11 配列の長さや、形状が異なる場合に行われる計算について

形状の異なる配列同士の演算を行う場合、それ専用の規則(broadcasting rule)に従い処理が行われます。 その規則に従いエラーになることもあれば計算結果が求まることもあります。

11.1 エラーになるか計算結果が求まるかの判定方法

配列 arr1 と arr2 のそれぞれの配列の shape を一番右から左にかけて以下条件を満たす時、計算結果が求まります。 条件を満たさない場合はエラーとなります。また、エラーにならず、形状が異なる場合は、お互いに異なる形状を補うよう値が補完され、計算されます(詳細は割愛します)。

(以下は「Broadcasting」(https://docs.scipy.org/doc/numpy/user/basics.broadcasting.html) より引用)

  1. 両方の値が同じ、または、
  2. 片方の数値が1

11.2 Broadcasting rule の適用例

  1. OK

    配列 shape
    arr1 10
    arr2 1
    結果 10
  2. OK

    配列 shape  
    arr1 10 1
    arr2 1 10
    結果 10 10
  3. OK

    配列 shape    
    arr1 10 10 1
    arr2 1 1 3
    結果 10 10 3
  4. NG

    配列 shape    
    arr1 10 10 3
    arr2 1 1 2
    結果 エラー    

12 速度比較

numpy を使用すると、どれくらい速くなるのでしょうか。 100x100の要素からなる2次元配列の各要素を +1 する処理を numpy を使用する前と後で処理時間を比較してみました。 処理時間の計測には、マジックコマンド %timeit を使用しています。 その結果、numpy の使用により処理速度は 100倍 になりました。また可読性も向上したと思います。

コード39. 「速度比較 numpy 使用前」:

def func(x):
   y = x.copy()
   for i in range(100):
       for j in range(100):
           y[i][j] += 1
   return y

# リストの初期化
x = [[0 for i in range(100)] for j in range(100)]
# 処理時間計測
%timeit func(x)
785 µs ± 24.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

コード40. 「速度比較 numpy 使用後」:

import numpy as np

def func(x):
  y = x.copy()
  y += 1
  return y

# 配列の初期化
x = np.zeros((100, 100))
# 処理時間計測
%timeit func(x)
7.11 µs ± 284 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

13 まとめ

配列の作成方法、配列同士の計算や、配列への関数の適用方法について学習しました。

numpy を使って得られる恩恵は以下の2つです。

  • 行列計算のコード量が少なくなる
  • 処理速度が早くなる

14 最後に

既存の行列の計算のソースコードの可読性の向上や、処理速度向上などを目的に、 numpy を使用したいと思われるかもしれません。 まず計算の対象が 10 の制約の範囲内で処理が可能であることを確認しましょう。 次に、以下の表のように可読性と、処理速度を天秤にかけ、 numpy で実装するか標準 Python で実装するか考えると良いのではないかと思います。 (numpy を使ったからといって常に「可読性が高い&処理速度が早い」コードになる訳ではありません。)

可読性 処理速度 numpy を使用する?
高くなる 高くなる YES
高くなる 低くなる YES or NO
低くなる 高くなる YES or NO
低くなる 低くなる NO

numpy を使用してコーディングするとテクニカルで難解なコードになってしまうことがあります。 可読性重視であればテクニカルなコーディングはしないほうが良いでしょう。 numpy を使用してコーディングする際は、numpy を使用する当初の目的を見失わないようにしましょう。

それでは、よい numpy ライフを!

著者: Kei Minagawa (http://keimina.hatenablog.jp)

Created: 2019-03-23 Sat 11:25

Validate

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. 出力

以上です。