データ解析用によく使われているデータベースMongoDBの設定方法を記録しています。

ダウンロード

ダウンロードリンク → https://www.mongodb.com/download-center/community

インストールの設定

パスの設定は C:\mongodb に設定（C:\Program Files ではない）

インストールは Completeではなく Custom に設定、Cドライブの下に新しいフォルダーを作る。↓

インストールは約5分かかる。

インストールができたら、C:\mongodb\dataフォルダーを作っておいてください。

サーバーを起動する方法

cd C:\mongodb\bin でダイレクトリーを移動し、
(これからよく使う方はC:\mongodb\binを環境係数に登録してもいいかもしれません)

データぱすの登録としてサーバーを起動する：
mongod.exe --dbpath C:\mongodb\data

以下の画面が出れば成功

よく失敗してしまう理由はC:\mongodb\dataが作成されないからなどにあると思います。

GUIでMongoDB の中身を確認する

MongoDB GUI の接続先は mongodb://localhost:27017/ に設定する。（サーバーが立ち上げないとと当然エラーが発生）

( mongodb:// は必要、Port番号27017はデフォルトの設定)

接続ができたら、データベーステーブルの画面が出れば成功。
（初めて起動するときはなにもないはずなので参考までにしてください）

Python で試す

Qiita : pythonでMongoDBを操作するはとてもわかりやすいので確認してみてください。

毎回Pythonで接続する前に、必ずサーバーを立ち上げておいてください。

参考

YouTube－MongoDB In 30 Minutes
https://www.youtube.com/watch?v=pWbMrx5rVBE

Qiita - pythonでMongoDB入門しようhttps://qiita.com/Syoitu/items/db192385a4b2e4884ed5#%E3%83%A9%E3%82%A4%E3%83%96%E3%83%A9%E3%83%AA%E3%81%AE%E3%82%A4%E3%83%B3%E3%82%B9%E3%83%88%E3%83%BC%E3%83%AB

概要

ロジャーや白ひげの親父などの懸賞金が公開され、海賊王を目指すルフィー君と大海賊たちの差はどれだけ離れているのか気になり、matplotlibの練習を兼ねて可視化してみることにしました。

環境

・Python 3.8.0
・windows10
・Jupyter Notebook

参考資料

・ドラゴンボールの戦闘力とmatplotlibで対数を勉強

・【2019年最新版】「ワンピース」懸賞金ランキング！ついにロジャーの懸賞金が発覚

大まかな流れ

1. 懸賞金が・キャラ名を取得(スクレイピング)
2. 可視化する(matplotlib)

1. 懸賞金が・キャラ名を取得(スクレイピング)

こちらの記事から、ワンピースの懸賞金を取得しました。

スクレイピングの流れは以下の通り。
1. 全キャラ懸賞金取得
2. 全キャラ名を取得
3. 億越えキャラのみ、キャラ名・懸賞金を取得

scraping_onepiece.py

import requests
from bs4 import BeautifulSoup
import pandas as pd
import re
import time

list_df = pd.DataFrame(columns=['キャラ', '懸賞金'])

url = 'https://ciatr.jp/topics/311415'
response = requests.get(url)
soup = BeautifulSoup(response.text, 'lxml')
# 全キャラ懸賞金を取得
bounty_tags = soup.find_all('div', class_='component ArticleComponentHead')
bounty_tags = soup.find_all('h3')

# h3要素に含まれる['作品情報', '記事情報', 'ciatrについて']は無視
bounty_tags = [tag.text for tag in bounty_tags if tag.text not in ['作品情報', '記事情報', 'ciatrについて']]

# 1秒待機
time.sleep(1)

# キャラ名取得
name_tags = soup.find_all('div', class_='component ArticleComponentHead scroll-margin')

name_tags = [tag.text for tag in name_tags if tag.text not in ['『ONE PIECE』に登場する並み居る強豪たちの懸賞金額が知りたい！',
                                                                '懸賞金不明な大物はまだいる！最高額はまだまだ更新されそう']]

# {キャラ名: 懸賞金}とするため辞書作成
master_table = dict()

for bounty_tag, name_tag in zip(bounty_tags, name_tags):

    # 正規表現で余計な文字削除
    bounty = re.sub(r'[懸賞金元ベリー以上]', '', bounty_tag)

    # 懸賞金が億超してないキャラは除外
    if '億' in bounty:

        #以降、懸賞金をmatplotlibで使うので整数化する 
        oku_ind = bounty.find("億")
        man_ind = bounty.find("万")

        oku = bounty[:oku_ind]

        if bounty[oku_ind + 1:man_ind] != '':
            man = bounty[oku_ind + 1:man_ind]
        else:
            man = 0

        bounty = int(oku) * 100000000 + int(man) * 10000

        name = re.sub(r'[【】0-9位]', '', name_tag)

        # 辞書に追加
        master_table[name] = bounty

for name, bounty in master_table.items():
    tmp_se = pd.DataFrame([name, bounty], index=list_df.columns).T
    list_df = list_df.append(tmp_se)

print(list_df)

# csv保存
list_df.to_csv('onepiece.csv', mode = 'a', encoding='cp932')

2. 可視化する(matplotlib)

malplotlib_onepiece.ipynb

#Jupyter Notebook上にグラフを描画する際に指定する記述
%matplotlib inline

import matplotlib.pyplot as plt
import pandas as pd

plt.rcParams["font.family"] = "IPAexGothic"

フォントエラーはこちらを参考にしました。

matplotlibの日本語文字化けを解消する(Windows編)

df_file = pd.read_csv('onepiece.csv', encoding='cp932')

df_name = df_file['キャラ'].tolist()
df_bounty = df_file['懸賞金'].tolist()

d = []

for name, bounty in zip(df_name, df_bounty):
    data = [name, bounty]
    d.append(data)

#懸賞金降順
name = list(reversed([x[0] for x in d]))
bounty = list(reversed([x[1] for x in d]))

#図全体のサイズ
plt.figure(figsize=(12, len(d)/2))
#横棒グラフ
plt.barh(name, bounty)
#グラフ全体のタイトル
plt.title('ワンピース懸賞金ランキング')
#xラベルのタイトル
plt.xlabel('懸賞金(億)')
#x軸（懸賞金）のスケールを設定
plt.xscale('log')
plt.show()

可視化してみると、四皇の強さ半端じゃないですね。まだ、未公開のドラゴンやレイリー師匠の懸賞金が楽しみです。個人的な意見なんですが、ドラゴンの懸賞金は世界政府を直接倒そうとしている組織のボスですから50億はいってほしいです。

おわりに

今回はワンピースの億越キャラの可視化を目的にスクレイピング、malplotlibを学びました。スクレイピングで文字列の漢数字を算用数字に変換する作業が難しく時間がかかりました。
次回は、教師なし学習のクラスタリングや感情分析を学びたいと思います。長くなりましたが、ここまで読んでくださりありがとうございます。誤っている箇所がございましたら、コメントでご指摘頂けると大変嬉しいです。

概要

kerasを使ったテキスト分類を試し、記事にまとめます。
データセットはtensorflowに内蔵されたロイター通信のデータセットです（英語のテキストデータ）。

Keras MLPの文章カテゴリー分類を理解する というブログ記事を参考に、一度取り組んだことがあります。
今回はドキュメントを引きつつ手を動かしており、理解を深める目的でこの記事をアウトプットします。
構築したモデルは、非常にシンプルなMLPです。

分量が長くなったので2つに分けます：

• 本記事で扱うこと
  • データセットについて
  • 前処理について
• 次の記事で扱うこと
  • モデルの学習について
  • モデルの性能評価について

動作環境

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.14.6
BuildVersion:   18G103
$ python -V  # venvモジュールによる仮想環境を利用
Python 3.7.3
$ pip list  # 主要なものを抜粋
ipython              7.11.0
matplotlib           3.1.2
numpy                1.18.0
pip                  19.3.1
scikit-learn         0.22.1
scipy                1.4.1
tensorflow           2.0.0

データセット

読み込み

tensorflow.keras.datasets.reuters.load_data(ドキュメント)で読み込むことができます。
test_split引数のデフォルト値が0.2のため、学習用8割、テスト用2割に分かれて読み込まれます。
※初回実行時は、データがダウンロードされます。

In [1]: from tensorflow.keras.datasets import reuters

In [2]: (x_train, y_train), (x_test, y_test) = reuters.load_data()

In [3]: len(y_train), len(y_test)
Out[3]: (8982, 2246)  # 合計 11228 件

ラベルを見る

ラベルはニュースのトピックを表すそうです。
試しにラベルを1つ見てみましょう。

In [4]: y_train[1000]
Out[4]: 19

数値で表されています（※それぞれがどんなトピックなのかまでは調べきれていません）。

学習用とテスト用のデータ全体で何種類のラベルがあるか確認します。
numpy.ndarrayのy_trainとy_testをリストに変換して、collections.Counter(ドキュメント)に渡します。

In [5]: from collections import Counter
In [8]: counter = Counter(list(y_train) + list(y_test))

In [9]: len(counter)
Out[9]: 46

全部で46のトピックがありました。

トピックごとに何件あるか確認します。

In [10]: for i in range(46):
    ...:     print(f'{i}: {counter[i]},')
    ...:
0: 67,
1: 537,
2: 94,
3: 3972,
4: 2423,
5: 22,
6: 62,
7: 19,
8: 177,
9: 126,
10: 154,
11: 473,
12: 62,
13: 209,
14: 28,
15: 29,
16: 543,
17: 51,
18: 86,
19: 682,
20: 339,
21: 127,
22: 22,
23: 53,
24: 81,
25: 123,
26: 32,
27: 19,
28: 58,
29: 23,
30: 57,
31: 52,
32: 42,
33: 16,
34: 57,
35: 16,
36: 60,
37: 21,
38: 22,
39: 29,
40: 46,
41: 38,
42: 16,
43: 27,
44: 17,
45: 19,

3と4のトピックが図抜けて多く、約57%を占めます。
トピックに含まれる記事の数に偏りがありますが、今回は46クラスへの分類という問題設定で進めます。

ニュースのテキストを見る

ニュースも1つ見てみましょう。

In [12]: x_train[1000]
Out[12]:
[1,
 437,
 495,
 1237,
 55,
 9070,
 :
 12]

整数からなるリストが表示されました。
今回のデータセットの場合、単語が数値に変換されています。
元のテキストを確認してみます。

まず、単語と数値の対応表は、tensorflow.keras.datasets.reuters.get_word_index(ドキュメント)で取得できます。
※初回実行時は、データがダウンロードされます。

In [16]: word_index = reuters.get_word_index()

In [17]: len(word_index)
Out[17]: 30979

In [18]: word_index
Out[18]:
{'mdbl': 10996,
 'fawc': 16260,
 'degussa': 12089,
 'woods': 8803,
 'hanging': 13796,
 'localized': 20672,
 :
 'hebei': 9407,
 ...}

In [19]: for word, index in word_index.items():
    ...:     if index in [0, 1, 2]:
    ...:         print(word, index)
    ...:
the 1
of 2

In [20]: for word, index in word_index.items():
    ...:     if index in [30978, 30979, 30980]:
    ...:         print(word, index)
    ...:
jung 30978
northerly 30979

word_indexは単語に対する数値の辞書です。
この対応を逆にして数値に対する単語の辞書を用意すればよさそうです。
ここで、x_trainとx_testに使われた整数は、word_indexの整数とずれていることに対応する必要があります。

ずれる理由は load_dataの3つの引数にあります。

1. 開始を表す数値：start_char=1

The start of a sequence will be marked with this character. Set to 1 because 0 is usually the padding character.

x_trainとx_testの中で1は開始を表します。
0がpadding character（埋め草文字。余白を埋めるための文字）に使われるため、1がデフォルト値となっているそうです。

2. 対応しない語を表す数値：oov_char=2

words that were cut out because of the num_words or skip_top limit will be replaced with this character.

num_wordsやskip_top引数によって、使う単語の範囲を区切ることで、対応しない語がoov_charに置き換えられます。
今のx_trainやx_testを取得する際、これらの引数を指定していないため、oov_charは現時点では無関係です。
（load_dataのドキュメントを見ると、oovはout of vocabularyの略のようです）

3. 単語に対応する数値の最初の値：index_from=3

index actual words with this index and higher.

数値のうち、0, 1, 2が意味を持っているため、単語の対応がずれるわけです。
index_from引数によりx_trainやx_testはword_indexの0が3に該当するという対応¹で読み込まれています。

ニュースを単語の並びとして見るために、数値のズレを考慮して、数値: 単語という辞書を作ります。

In [22]: index_word_map = {
    ...:     index + 3: word for word, index in word_index.items()
    ...: }
    ...: index_word_map[0] = "[padding]"
    ...: index_word_map[1] = "[start]"
    ...: index_word_map[2] = "[oov]"
In [23]: len(index_word_map)
Out[23]: 30982

この辞書を使うことでx_train, x_testの整数の並びから文章を復元することができました。

In [24]: for index in x_train[1000]:
    ...:     print(index_word_map[index], end = " ")
    ...:
[start] german banking authorities are weighing rules for banks' off balance sheet activities in an attempt to cope with the growing volume of sophisticated capital market instruments banking sources said interest rate and currency swaps and ...

前処理

テキストの前処理

整数で表されたニュース記事の長さはまちまちです。

In [31]: for x in x_train[998:1003]:
    ...:     print(len(x))
    ...:
133
51
626
17
442

そこで長さが揃うように変換して前処理します。
今回は、各ニュース記事を、ニュース記事全体に登場する頻度の上位1000語が含まれるか否かで表します。

例えば、上位1000語の中に「currency」という単語があり、対応する整数は500とします。
各記事を0か1の並びで表すとき、currencyという語を含む記事は、インデックス500に1が来ます。
currencyという語を含まない記事は、インデックス500が0です。
これが他の語にも当てはまります。

この変換により、

• ニュース記事の長さが1000に揃います
• ニュース記事は1000個の0または1の並びで表されます

load_dataメソッドのnum_words引数²に1000を渡して、登場する頻度の上位1000語でx_train, x_testを表すように変換します。

In [32]: (x_train, y_train), (x_test, y_test) = reuters.load_data(num_words=1000)
    ...:

In [33]: len(y_train), len(y_test)
Out[33]: (8982, 2246)

num_wordsを登場頻度の上位1000語としたので、そこに含まれない語はoov_char(整数では2)としてx_train, x_testで表されます。

In [34]: for index in x_train[1000]:
    ...:     print(index_word_map[index], end = " ")
    ...:
[start] german banking [oov] are [oov] [oov] for [oov] off balance [oov] [oov] in an [oov] to [oov] with the growing volume of [oov] capital market [oov] banking sources said interest rate and currency [oov] and ...

この段階ではニュース記事の長さはまだ揃っていません。
上位1000語がニュース記事に含まれるか否かを表すために、tensorflow.keras.preprocessing.text.Tokenizer(ドキュメント)を使います。

Tokenizerの初期化でnum_words引数に1000を渡します。
Tokenizerを使った処理でnum_words-1の語が考慮されます。

num_words: the maximum number of words to keep, based
on word frequency. Only the most common num_words-1 words will
be kept.

In [37]: from tensorflow.keras.preprocessing.text import Tokenizer

In [42]: tokenizer = Tokenizer(1000)

load_dataで上位1000語を取り出しているので、x_train, x_testに含まれる整数の最大は999です³。

In [71]: max_index = 0

In [72]: for x in list(x_train)+list(x_test):
    ...:     now_max = max(x)
    ...:     if now_max > max_index:
    ...:         max_index = now_max
    ...:

In [73]: max_index
Out[73]: 999

sequences_to_matrixメソッド(ドキュメント)で、x_train, x_testをそれぞれ変換します。

In [36]: x_train.shape, x_test.shape
Out[36]: ((8982,), (2246,))

In [76]: x_train = tokenizer.sequences_to_matrix(x_train, "binary")

In [77]: x_test = tokenizer.sequences_to_matrix(x_test, "binary")

In [78]: x_train.shape, x_test.shape
Out[78]: ((8982, 1000), (2246, 1000))

全てのニュース記事が長さが1000で表されました。

sequences_to_matrixのドキュメントによると、

a sequence is a list of integer word indices

すなわち、sequenceとは「単語を表す整数のリスト（意訳）」なので、x_train, x_testはまさしくsequenceです。

第2引数の指定ですが、"binany"の場合は各語が存在するかしないかの0/1で表されます。
他に、"count", "tfidf", "freq"を指定できるそうです。

変換されたニュース記事を試しに1つ見てみると

In [58]: x_train[1000]
Out[58]:
array([0., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 0., 0., 1.,
       ...

と0/1で表現されています。
index_word_mapの中で0（埋め草文字にあたる）や、3（もともと0の語がないので、3というキーがない）はどのニュースにも登場しないので0です。
1([start])や2([oov])、4(the)、5(of)などはx_train[1000]に登場するので1となっています。

テキストの前処理は以上です。

ラベルの前処理

ラベル（ニュースのトピック）は0〜45のいずれかですが、これをone-hot表現に変換します。
（one-hot表現とする理由は、ラベル同士に大小関係を持たせないようにするため）

tensorflow.keras.utils.to_categorical(ドキュメント)を使います。

In [80]: y_train.shape, y_test.shape
Out[80]: ((8982,), (2246,))

In [81]: from tensorflow import keras

In [85]: number_of_classes = len(counter)

In [86]: y_train = keras.utils.to_categorical(y_train, number_of_classes)

In [88]: y_test = keras.utils.to_categorical(y_test, number_of_classes)

In [89]: y_train.shape, y_test.shape
Out[89]: ((8982, 46), (2246, 46))

ラベルをインデックスと見立てて、ラベルのインデックスだけ1、他は0という形式に変換されます。
y_train[1000]は19でしたが、y_train[1000][19]が1、他は0となるように変換されています。

In [90]: y_train[1000]
Out[90]:
array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], dtype=float32)

前処理は以上です。

本記事のまとめ

• データセットについて
  • ロイター通信のニュースのトピック分類データ（多クラス分類）
  • ラベルは46クラスあり、含まれるニュースの件数に偏りがある
  • ニュースは整数のリストで表されている。単語の並びへ復元して元のニュースを確認できる
• 前処理について
  • テキストを固定長の0/1の並びに変換（登場頻度上位1000語のそれぞれが含まれるか否か）
  • ラベルをone-hot表現に変換

本記事は[モデル構築編]（近日公開）に続きます。

---

1. word_indexは0というキーを持ちませんが、1は4、2は5、、、のように対応するということです ↩

2. 「max number of words to include. Words are ranked by how often they occur (in the training set) and only the most frequent words are kept」（学習用データの中で頻出する上位何語を使って、x_train, x_testとして取り出すかの指定です） ↩

3. index_word_mapでキーが999までの単語ということです ↩

Pythonをぶっ壊〜〜〜す。
どうもこんにちは、「Pythonの被害から国民を守る党」党首を自称しております片岡です。
Pythonは嫌いで嫌いで正直一生書きたくないのですが、なんやかんやで書く羽目になりますよね。
ポストPythonなスクリプト言語よ早く生まれてPythonを滅ぼしておくれ・・・・・
と嘆いたところでPythonの快進撃は現状止まる気配がないのでうまいこと付き合っていかないといけません
ホントは明日にでも消えてほしいけど
まあPythonの中にも一部捨てたもんじゃない部分もあってそのおいしいとこだけはちゃんと活用していこうよということですね。
今回は普段良く使うROS1/ROS2においてPythonをなんとかうまく使っていく方法を紹介したいと思います。

ROS1におけるPython

ROS1にはPythonクライアントであるrospyがあります。
rospyの内部実装を追ってみるとSubscriberのオブジェクト作った瞬間からspinしなくても実はコールバックが回ってる、
spinは単にrospy.ok()みたいなことしてspinしてるだけとか怪しげな実装がコアライブラリに散見されます。大丈夫なんかコレ。。。
シリアライズされたrosメッセージをやり取りするので比較的すくないオーバーヘッドで他のROSノードと通信できます。
Pythonはスクリプト言語なのでビルドしないにもかかわらずビルドシステムにcatkinを使用し、しばしばrosdepによってpipに依存関係が崩壊します。
なお基本的にPython2のためrospyで実装したコードは今後技術的負債になる可能性は高いと思われます。
またROS1 -> ROS2の変換に関してはほぼリビルドになるのでそこも将来的に苦労が約束されています

ROS2におけるPython

ROS2にはPythonクライアントであるrclpyがあります。
シリアライズされたrosメッセージをやり取りするので比較的すくないオーバーヘッドで他のROSノードと通信できます。
ちなみに、現状rclpyのパッケージを作るのは結構一苦労（ros2 pkg createするとC++向けパッケージが作られるため手動でいろいろ追記しないといけない）です。
Pythonはスクリプト言語なのでビルドしないにもかかわらずビルドシステムにcolconを使用し、しばしばrosdepによってpipに依存関係が崩壊します。

なんとかならんのか・・・・

これらの手段を見て僕が思ったのは・・・・
なんだこの新手の罰ゲームは・・・・という印象ですね。
数年前の僕はこれらのことにPythonの記法から生まれる性質（長いコードを書くと可読性が低い、コンパイルがないのでテストできっちり網羅しない限り最低限のチェックすらままならない）
を鑑みてROSノードを書く時はPythonを完全排除するという考えにいたりました。
一応ある程度の成果は出て、いろいろ使いやすいROSノードが生み出せたかなと思います。
まあでもしかし、ツール系のものを整備するにはそうも行ってられません。
C++でweb APIとか叩くのほんとめんどい、となるとそういうことに強いPythonが選択肢として上がってきます。
うーんでも上記の理由でrospy等はできれば使いたくない、さらにあわよくばpipによる環境破壊を避けるため
Pythonのvenvとか依存解決ツールとかも使って効率的に開発したい・・・・・
あわよくばROS1/ROS2で入る破壊的変更でコードの書き直しが発生するのは嫌だ・・・・
となると理想形の手法はこういうことになります・・・・

「PythonとROS1/2の間で通信はできるが、パッケージ自体はROSに依存しないパッケージをつくる」

んなもんできるんか？？

できるんだなコレが・・・・
ということでサンプルのパッケージを作ってみました。
依存解決ツールにはpipenvを使用しました。
pipenvはpip+venvをパッケージングしたツールでパッケージごとに仮想環境を作って依存解決をかんたんにやってくれます。
詳しくはこちらの記事を参照ください
以下のコマンドをターミナルに入力します。

mkdir roslibpy_example
cd roslibpy_example
pipenv install

すると以下のような出力が得られます。

Creating a virtualenv for this project...
Pipfile: /home/masaya/lib/roslibpy_example/Pipfile
Using /usr/bin/python (2.7.15+) to create virtualenv...
⠏ Creating virtual environment...Already using interpreter /usr/bin/python
New python executable in /home/masaya/.local/share/virtualenvs/roslibpy_example-mGuxztlN/bin/python
Installing setuptools, pip, wheel...
done.

✔ Successfully created virtual environment! 
Virtualenv location: /home/masaya/.local/share/virtualenvs/roslibpy_example-mGuxztlN
Creating a Pipfile for this project...
Pipfile.lock not found, creating...
Locking [dev-packages] dependencies...
Locking [packages] dependencies...
Updated Pipfile.lock (dfae9f)!
Installing dependencies from Pipfile.lock (dfae9f)...
  ?   ▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉ 0/0 — 00:00:00
To activate this project's virtualenv, run pipenv shell.
Alternatively, run a command inside the virtualenv with pipenv run.

これでROSとの依存が切れたPythonパッケージが作れました。

ではこのパッケージの依存にROSとの通信を担当するroslibpyパッケージを追加しましょう。
以下のコマンドをターミナルで実行すると仮想環境にroslibpyがインストールされます。

pipenv install roslibpy

すると以下のような表示が現れインストールが完了します。

Installing roslibpy...
Adding roslibpy to Pipfile's [packages]...
✔ Installation Succeeded 
Pipfile.lock (134db9) out of date, updating to (dfae9f)...
Locking [dev-packages] dependencies...
Locking [packages] dependencies...
✔ Success! 
Updated Pipfile.lock (134db9)!
Installing dependencies from Pipfile.lock (134db9)...
  ?   ▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉▉ 18/18 — 00:00:01
To activate this project's virtualenv, run pipenv shell.
Alternatively, run a command inside the virtualenv with pipenv run.

roslibpyはrosbridge protocolによってROS1/ROS2と通信を行います。
rosbridge protocolはROSのメッセージをjson文字列にしてwebsockやその他の手法で投げつけるプロトコルです。
詳しくはこちらの記事を参考にしてください。

あとはこの記事のようにroslibpyをimportしたスクリプトを作って実行するだけです。
ROSライクにかけるのでとてもかんたんですね。

# -*- coding:utf8 -*-
import time
from roslibpy import Message, Ros, Topic
def main():

    # roscoreを実行しているサーバーへ接続
    ros_client = Ros('192.168.1.104', 9090)
    # Publishするtopicを指定
    publisher = Topic(ros_client, '/turtle1/cmd_vel', 'geometry_msgs/Twist')
    def start_sending():
        while True:
            if not ros_client.is_connected:
                break

            # 送信するTwistメッセージの内容
            publisher.publish(Message({
                'linear': {
                    'x': 2.0,
                    'y': 0,
                    'z': 0
                },
                'angular': {
                    'x': 0,
                    'y': 0,
                    'z': 1.8
                }
            }))
            time.sleep(0.1)
        publisher.unadvertise()
    # Publish開始
    ros_client.on_ready(start_sending, run_in_thread=True)
    ros_client.run_forever()
if __name__ == '__main__':
    main()

出典：https://symfoware.blog.fc2.com/blog-entry-2288.html

最後に、setup.py等を書いたらパッケージの完成です。
いろんな人がつかえるようにPypyとかに登録しちゃいましょう。

まとめ

以上のようなやり方を取るとROS1/2両方に対応できてROSに依存のないROSと通信できるPythonパッケージを作ることが可能です。
最後にこのやり方に対して来るであろう「大きいメッセージをこの方法でやり取りたらレイテンシがやばいでしょ」というツッコミに対して一言。
パフォーマンス気にするならおとなしくPythonをやめましょう。
Pythonをぶっ壊〜〜〜す。

TL;DR

• 書籍「効果検証入門 正しい比較のための因果推論／計量経済学の基礎」のRソースコードを、Pythonで（ほぼ）再現しました
• https://github.com/nekoumei/cibook-python
• 本記事では、主にRではライブラリどーん！で済むけどPythonではそうはいかない部分の解説をします

書籍の紹介

https://www.amazon.co.jp/dp/B0834JN23Y

上記Amazonに目次が載っているのでそれを見るのが早い気がしますが。。
とても良い本です。正確な意思決定を行うためにどうやってバイアスを取り除くか？に焦点を当てて種々の因果推論の手法（傾向スコア/DiD/RDDなど）をRソースコードによる実装とともに紹介されています。
全体を通して、現実問題の効果検証に因果推論を活用するにはどうすればいいか？という観点で書かれており、非常に実用的だなーと感じました。

Pythonで書きました

https://github.com/nekoumei/cibook-python
大元の公開されているRソースコード(https://github.com/ghmagazine/cibook )と対応する形で、Jupyter Notebookを作成しました。
また、本Python実装にあたって、グラフの可視化ライブラリは一部を除いてplotly.expressを使用しています。
GithubのNotebookレンダリングではplotlyのグラフは表示できないので、オンラインで確認したい場合はREADMEに記載しているGithub Pagesをご確認ください。（images配下のhtmlを表示しています）
plotly.expressに関しては、日本語だと下記記事がとても参考になります。
令和時代のPython作図ライブラリのデファクトスタンダードPlotlyExpressの基本的な描き方まとめ

Pythonで書くうえでの主要なポイント

回帰分析

statsmodelsのOLSおよびWLSを利用しています。

(ch2_regression.ipynbより抜粋)

## バイアスのあるデータでの重回帰
y = biased_df.spend
# R lmではカテゴリ変数は自動的にダミー変数化されているのでそれを再現
X = pd.get_dummies(biased_df[['treatment', 'recency', 'channel', 'history']], columns=['channel'], drop_first=True)
X = sm.add_constant(X)
results = sm.OLS(y, X).fit()
nonrct_mreg_coef = results.summary().tables[1]
nonrct_mreg_coef

ここで、Rのlm関数との違いが2つあります。
①lmではカテゴリカル変数のダミー変数化を自動で行っているようです。sm.OLSでは自動でできないので pd.get_dummies()を利用しています。
②lmではバイアス項が自動で追加されます。こちらも手動で追加します。（参考：統計: Python と R で重回帰分析してみる）

RData形式のデータセットの読込

実装内で利用されているデータセットの中には、experimentdatarのようにRパッケージとして公開されていたり、RData形式で公開されているデータセットがあります。
これらをPythonのみで読み込む方法はないです。（あったら教えてください）
コメント欄にて、upuraさんにご教示いただきました！
https://qiita.com/nekoumei/items/648726e89d05cba6f432#comment-0ea9751e3f01b27b0adb
.rdaファイルは、rdataというパッケージで、Rなしで読み込みが可能です。
(ch2_voucher.ipynbより抜粋)

parsed = rdata.parser.parse_file('../data/vouchers.rda')
converted = rdata.conversion.convert(parsed)
vouchers = converted['vouchers']

rpy2経由でRを使ってデータを読込、pandas DataFrameへの変換を行う場合は下記のとおりです。
(ch2_voucher.ipynbより抜粋)

from rpy2.robjects import r, pandas2ri
from rpy2.robjects.packages import importr
pandas2ri.activate()
experimentdatar = importr('experimentdatar')
vouchers = r['vouchers']

ch2_voucher.ipynbに記載しているとおり、Rパッケージのインストール等は事前にRの対話環境で行っています。
また、ch3_lalonde.ipynbにて使用するデータセットは.dtaファイルです。
こちらは、pandasのread_stata()によって読み込むことができます。
(ch3_lalonde.ipynbより抜粋)

cps1_data = pd.read_stata('https://users.nber.org/~rdehejia/data/cps_controls.dta')
cps3_data = pd.read_stata('https://users.nber.org/~rdehejia/data/cps_controls3.dta')
nswdw_data = pd.read_stata('https://users.nber.org/~rdehejia/data/nsw_dw.dta')

傾向スコアマッチング（最近傍マッチング）

傾向スコアマッチングのマッチング手法はいくつかありますが、書籍ではMatchItの最近傍マッチングを使っているようです。
Pythonにはちょうどいいライブラリはなさそうだったので愚直に実装しています。
(ch3_pscore.ipynbより抜粋)

def get_matched_dfs_using_propensity_score(X, y, random_state=0):
    # 傾向スコアを計算する
    ps_model = LogisticRegression(solver='lbfgs', random_state=random_state).fit(X, y)
    ps_score = ps_model.predict_proba(X)[:, 1]
    all_df = pd.DataFrame({'treatment': y, 'ps_score': ps_score})
    treatments = all_df.treatment.unique()
    if len(treatments) != 2:
        print('2群のマッチングしかできません。2群は必ず[0, 1]で表現してください。')
        raise ValueError
    # treatment == 1をgroup1, treatment == 0をgroup2とする。group1にマッチするgroup2を抽出するのでATTの推定になるはず
    group1_df = all_df[all_df.treatment==1].copy()
    group1_indices = group1_df.index
    group1_df = group1_df.reset_index(drop=True)
    group2_df = all_df[all_df.treatment==0].copy()
    group2_indices = group2_df.index
    group2_df = group2_df.reset_index(drop=True)

    # 全体の傾向スコアの標準偏差 * 0.2をしきい値とする
    threshold = all_df.ps_score.std() * 0.2

    matched_group1_dfs = []
    matched_group2_dfs = []
    _group1_df = group1_df.copy()
    _group2_df = group2_df.copy()

    while True:
        # NearestNeighborsで最近傍点1点を見つけ、マッチングする
        neigh = NearestNeighbors(n_neighbors=1)
        neigh.fit(_group1_df.ps_score.values.reshape(-1, 1))
        distances, indices = neigh.kneighbors(_group2_df.ps_score.values.reshape(-1, 1))
        # 重複点を削除する
        distance_df = pd.DataFrame({'distance': distances.reshape(-1), 'indices': indices.reshape(-1)})
        distance_df.index = _group2_df.index
        distance_df = distance_df.drop_duplicates(subset='indices')
        # しきい値を超えたレコードを削除する
        distance_df = distance_df[distance_df.distance < threshold]
        if len(distance_df) == 0:
            break
        # マッチングしたレコードを抽出、削除する
        group1_matched_indices = _group1_df.iloc[distance_df['indices']].index.tolist()
        group2_matched_indices = distance_df.index
        matched_group1_dfs.append(_group1_df.loc[group1_matched_indices])
        matched_group2_dfs.append(_group2_df.loc[group2_matched_indices])
        _group1_df = _group1_df.drop(group1_matched_indices)
        _group2_df = _group2_df.drop(group2_matched_indices)

    # マッチしたレコードを返す
    group1_df.index = group1_indices
    group2_df.index = group2_indices
    matched_df = pd.concat([
        group1_df.iloc[pd.concat(matched_group1_dfs).index],
        group2_df.iloc[pd.concat(matched_group2_dfs).index]
    ]).sort_index()
    matched_indices = matched_df.index

    return X.loc[matched_indices], y.loc[matched_indices]

treatment群の1点と最も傾向スコアが近いcontrol群1点をマッチングさせる、を反復しています。その際、距離がstdの0.2倍に収まるペアのみ抽出するようしきい値を設けています。
このあたりの詳細は傾向スコアの概念とその実践が詳しかったです。
ch3_pscore.ipynbにてMatchItによるマッチング結果との比較をしていますが、完全に一致はしていません。結論が変わらない且つ共変量のバランスも良い感じ(下記参照)のため概ね良しとしています。
こちらも、Rのcobalt love.plot()のような便利ライブラリはないので自力で可視化しています。
(images/ch3_plot1.htmlより)

逆確率重み付き推定（IPW）

IPWは実装がシンプルなのがいいところですね。
こちらもWeightItと比較していますが、概ね正しそうです。
(ch3_pscore.ipynbより抜粋)

def get_ipw(X, y, random_state=0):
    # 傾向スコアを計算する
    ps_model = LogisticRegression(solver='lbfgs', random_state=random_state).fit(X, y)
    ps_score = ps_model.predict_proba(X)[:, 1]
    all_df = pd.DataFrame({'treatment': y, 'ps_score': ps_score})
    treatments = all_df.treatment.unique()
    if len(treatments) != 2:
        print('2群のマッチングしかできません。2群は必ず[0, 1]で表現してください。')
        raise ValueError
    # treatment == 1をgroup1, treatment == 0をgroup2とする。
    group1_df = all_df[all_df.treatment==1].copy()
    group2_df = all_df[all_df.treatment==0].copy()
    group1_df['weight'] = 1 / group1_df.ps_score
    group2_df['weight'] = 1 / (1 - group2_df.ps_score)
    weights = pd.concat([group1_df, group2_df]).sort_index()['weight'].values
    return weights

CausalImpact

ch4_did.ipynbに記載しているとおり、下記2つのライブラリを使い、比較しています。

• dafiti/causalinpact: https://github.com/dafiti/causalimpact
• tcassou/causal_impact: https://github.com/tcassou/causal_impact
  

outputがキレイなのでいつもはdafitiのcausalimpactを使っているんですが、推定結果が書籍（Rの本家causalimpact）に近いのはtcassouのcausal_impactでした。
両者とも、statsmodelsの状態空間モデルで推定してるみたいですが、、ちょっと実装の違いはよくわからなかったです。誰か教えて下さい。。
どちらも、R実装と比較して推定誤差がめちゃめちゃ小さくなってるんですよね。共変量が多いのがあまりよくないのかもしれない。

dafiti/causalimpactのplot

tcassou/causal_impactのplot

回帰不連続デザイン(RDD)

Rだとrddtoolsでサクっと実行できますが、Pythonだとそうはいかないです。
まず、rddtoolsのrdd_reg_lmで使われている回帰式を確認します。
（参考：https://cran.r-project.org/web/packages/rddtools/rddtools.pdf P23)

Y = α + τD + β_1(X-c)+ β_2D(X-c) + ε

余談ですが、私はてっきりRDDって、cut-off値の左と右でそれぞれ回帰モデルつくって、cut-off値での推定値の差分を取るんだと思ってたんですが、1つの回帰式で表すんですね。意味的には同じなのかな。
ここで、DはXがcut-off値以降なら1, 以前なら0の値を持つバイナリ変数です。Dのcoefが効果量として確認したい値です。
また、cはcut-off値です。
以上を踏まえて、実装します。実装にあたってはrddtoolsのソースコードも参照しています。(https://github.com/MatthieuStigler/RDDtools/blob/master/RDDtools/R/model.matrix.RDD.R )
(ch5_rdd.ipynbより抜粋)

class RDDRegression:
# Rパッケージrddtoolsのrdd_reg_lmを再現する
# 参考：https://cran.r-project.org/web/packages/rddtools/rddtools.pdf P23
    def __init__(self, cut_point, degree=4):
        self.cut_point = cut_point
        self.degree = degree

    def _preprocess(self, X):
        X = X - threshold_value
        X_poly = PolynomialFeatures(degree=self.degree, include_bias=False).fit_transform(X)
        D_df = X.applymap(lambda x: 1 if x >= 0 else 0)
        X = pd.DataFrame(X_poly, columns=[f'X^{i+1}' for i in range(X_poly.shape[1])])
        X['D'] = D_df
        for i in range(X_poly.shape[1]):
            X[f'D_X^{i+1}'] = X_poly[:, i] * X['D']
        return X

    def fit(self, X, y):
        X = X.copy()
        X = self._preprocess(X)
        self.X = X
        self.y = y
        X = sm.add_constant(X)
        self.model = sm.OLS(y, X)
        self.results = self.model.fit()
        coef = self.results.summary().tables[1]
        self.coef = pd.read_html(coef.as_html(), header=0, index_col=0)[0]

    def predict(self, X):
        X = self._preprocess(X)
        X = sm.add_constant(X)
        return self.model.predict(self.results.params, X)

多項式回帰を行うための前処理として、scikit-learnのPolynomialFeaturesを利用しています。
(images/ch5_plot2_3.htmlより)

良いかんじに推定できました。Notebookを見ていただければわかりますが、効果量の推定も書籍とほぼ一致しています。よかった。

nonparametric RDD (RDestimate)

できませんでした、、、冒頭の「Pythonで（ほぼ）再現しました」のほぼの部分です。
RDestimateでは、Imbens and Kalyanaraman(2012) Optimal Bandwidth Choice for the Regression Discontinuity Estimator の手法を用いて、最適なバンド幅を選択するんですが、そのバンド幅の推定がどうやってるのかよくわかんなかったです。
ch5_rdd.ipynbでは、雑にバンド幅を変えたときのMSEが最適なバンド幅を見つけるように実装してみましたが、あまりうまく行っているように見えないです。
(ch5_plot4.htmlより)

うーんってかんじですね。純粋にバンド幅を狭めて推定するとバンド幅外は予測できなくなるのはそれはそうな気がするんですが、実際どうやるんですかね。
あるいは、cut-off近傍の推定にだけ興味があるからあまり気にしなくていいのかな。

おわりに

理解や実装が間違っていたり、なんかおかしいところとかがあったらぜひ教えてほしいです。
Twitter

ChainerからPytorchへの移植

研究でDeep Learningをしているのですが、先日Chainerのアップデートが終わりを迎えるのを知り、開発元と同様Pytorchにフレームワークを変更することになりました。
手始めに今あるChainerのプログラムからPytorchに移植することにしました。

基本的には関数の名前などを変えるだけでよかったんですが、途中でPytorchにHardSigmoidがないことに気づきました。
ということで自分で作っちゃおうということです。

実際にかいてみた

...といっても、公式リファレンスに書いてあるのでほぼその通りにやっただけです。
--> https://pytorch.org/docs/master/autograd.html

class MyHardSigmoid(torch.autograd.Function):
    @staticmethod
     def forward(ctx, i):
         ctx.save_for_backward(i)
         result = (0.2 * i + 0.5).clamp(min=0.0, max=1.0)
         return result

     @staticmethod
     def backward(ctx, grad_output):
         grad_input = grad_output.clone()
         result, = ctx.saved_tensors
         grad_input *= 0.2
         grad_input[result < -2.5] = 0
         grad_input[result > -2.5] = 0
         return grad_input

@staticmethodを書かないとwarningが出てきます。
公式のものは、指数関数になっていますが、これをhard sigmoidに変えていきます。

まずforward()には、順伝搬をかきます。
hard sigmoid()は以下のような式になるのでそうなるように書きました。

h(x) = \left\{
\begin{array}{ll}
0 & (x \lt -2.5) \\
0.2x + 0.5 & (-2.5 \leq x \leq 2.5)  \\
1 & (2.5 \lt x)
\end{array}
\right.

次にbackward()で、これは逆伝搬を書きます。
微分係数は以下のようになります。

\frac{\partial h(x)}{\partial x} = \left\{
\begin{array}{ll}
0 & (x \lt -2.5) \\
0.2 & (-2.5 \leq x \leq 2.5)  \\
0 & (2.5 \lt x)
\end{array}
\right.

そして最後にこれをモデルに適用させます。（モデルの中身は適当です。）

model.py

import torch.nn as nn
import torch.nn.functional as F

class Model(nn.Module):
    def __init__(self):
        super(Model, self).__init__()
        self.conv1 = nn.Conv2d(1, 20, 5)
        self.conv2 = nn.Conv2d(20, 20, 5)

    def forward(self, x):
        x = F.relu(self.conv1(x))
        hard_sigmoid = MyHardSigmoid.apply
        return hard_sigmoid(self.conv2(x))

これでバッチリ！！...なはず

はじめに

ふと思い立って勉強を始めた「ゼロから作るDeep LearningーーPythonで学ぶディープラーニングの理論と実装」の4章で私がつまずいたことのメモです。

実行環境はmacOS Mojave + Anaconda 2019.10です。詳細はこのメモの1章をご参照ください。

4章 ニューラルネットワークの学習

この章はニューラルネットワークの学習についての説明です。

4.1 データから学習する

通常は人が規則性を導き出してアルゴリズムを考え、それをプログラムに書いてコンピューターに実行させます。このアルゴリズムを考える作業自体もコンピューターにやらせてしまおうというのが、機械学習やニューラルネットワーク、ディープラーニングです。

この本では、処理したいデータに対して、事前に人が考えた特徴量の抽出（ベクトル化など）が必要なものを「機械学習」、さらにその「機械学習」に特徴量の抽出まで任せて生データをそのまま渡せるようにしたものを「ニューラルネットワーク（ディープラーニング）」と定義しています。この定義はやや乱暴な感じもしますが、言葉の使い分けにはあまり興味がないので気にせず先に進みます。

訓練データ、テストデータ、過学習などについて解説されていますが、特につまずく部分はありませんでした。

4.2 損失関数

損失関数として良く使われる2乗和誤差と交差エントロピー誤差の解説、そして訓練データの一部を使って学習するミニバッチ学習の解説です。ここも特につまずく部分はありませんでした。訓練データを全数使っても良さそうですが時間がかかって非効率ということですね。いわゆる標本調査みたいなことなのかと思います。

また、損失関数の代わりに認識精度を使うことができない理由として、認識精度は結果の微小な変化で反応せず不連続に変化するためうまく学習できないと解説されています。最初はピンと来ないかも知れませんが、次の微分の説明が終わると腹落ちするかと思います。

4.3 数値微分

微分の解説です。実装時の丸め誤差の説明は実用的でありがたいです。「微分」とか「偏微分」とか言葉を聞くと難しそうに思えてしまいますが、ちょっと値を変化させたら結果はどう変わるか？というだけなので、特に高校数学のおさらいとかはしなくても前に進めます。

ちなみに微分で出てくる $ \partial $ という記号はWikipedia曰くデルとかディーとかパーシャル・ディーとかラウンド・ディーとか読むそうです。

それにしても、Pythonは引数に関数が簡単に渡せていいですね。私のプログラマー現役時代はC/C++メインだったのですが、関数ポインタの表記がホント分かりにくくて嫌いでした

4.4 勾配

すべての変数の偏微分をベクトルにしたものが勾配です。これ自体は難しくはありません。

NumPy配列で小数を出力する時に値を丸めて表示してくれるのは見やすくていいですね。

>>> import numpy as np
>>> a = np.array([1.00000000123, 2.99999999987])
>>> a
array([1., 3.])

でも勝手に丸められると困ることもあるし、どんな仕様なんだろうと調べて見たら、表示方法を設定する機能がありました。numpy.set_printoptionsで、小数の表示方法や、要素数が多い場合の省略方法などを変更できます。例えばprecisionで小数点以下の桁数を大きく指定すると、きちんと丸められずに表示されます。

>>> np.set_printoptions(precision=12)
>>> a
array([1.00000000123, 2.99999999987])

これは便利！

4.4.1 勾配法

文中で「勾配降下法（gradient descent method）」という言葉が出てきますが、これは、以前勉強した時の教材では「最急降下法」と訳されていたものでした。

あと、学習率を示す $ \eta $という記号が出てきますが、これはギリシャ文字でイータと読みます（以前勉強した時は読み方を覚えていたのですが、すっかり忘れていてググりました）。

4.4.2 ニューラルネットワークに対する勾配

numerical_gradient(f, x)を使って勾配を求めるのですが、このfに渡す関数が

def f(W):
    return net.loss(x, t)

となっていて、あれ？この関数は引数Wを使っていない？と少し混乱しましたが、「4.4 勾配」のところで実装したnumerical_gradient(f, x)の関数の形をそのまま使おうとしているためで、引数Wはダミーとのこと。確かにsimpleNetクラスは自身で重みWを保持しているので、損失関数simpleNet.lossに重みWを渡す必要はありません。ダミーがあると分かりにくいので、私は引数なしで実装してみることにしました。

あとここで、numerical_gradientを多次元配列でも大丈夫な形に修正しておく必要があります。

4.5 学習アルゴリズムの実装

ここからは、これまでに学んだ内容を使って、実際に確率的勾配降下法（SGD）を実装します。

まず、必要な関数を寄せ集めたfunctions.pyです。

functions.py

# coding: utf-8
import numpy as np


def sigmoid(x):
    """シグモイド関数
    本の実装ではオーバーフローしてしまうため、以下のサイトを参考に修正。
    http://www.kamishima.net/mlmpyja/lr/sigmoid.html

    Args:
        x (numpy.ndarray): 入力

    Returns:
        numpy.ndarray: 出力
    """
    # xをオーバーフローしない範囲に補正
    sigmoid_range = 34.538776394910684
    x2 = np.maximum(np.minimum(x, sigmoid_range), -sigmoid_range)

    # シグモイド関数
    return 1 / (1 + np.exp(-x2))


def softmax(x):
    """ソフトマックス関数

    Args:
        x (numpy.ndarray): 入力

    Returns:
        numpy.ndarray: 出力
    """
    # バッチ処理の場合xは(バッチの数, 10)の2次元配列になる。
    # この場合、ブロードキャストを使ってうまく画像ごとに計算する必要がある。
    # ここでは1次元でも2次元でも共通化できるようnp.max()やnp.sum()はaxis=-1で算出し、
    # そのままブロードキャストできるようkeepdims=Trueで次元を維持する。
    c = np.max(x, axis=-1, keepdims=True)
    exp_a = np.exp(x - c)  # オーバーフロー対策
    sum_exp_a = np.sum(exp_a, axis=-1, keepdims=True)
    y = exp_a / sum_exp_a
    return y


def numerical_gradient(f, x):
    """勾配の算出

    Args:
        f (function): 損失関数
        x (numpy.ndarray): 勾配を調べたい重みパラメーターの配列

    Returns:
        numpy.ndarray: 勾配
    """
    h = 1e-4
    grad = np.zeros_like(x)

    # np.nditerで多次元配列の要素を列挙
    it = np.nditer(x, flags=['multi_index'])
    while not it.finished:

        idx = it.multi_index  # it.multi_indexは列挙中の要素番号
        tmp_val = x[idx]  # 元の値を保存

        # f(x + h)の算出
        x[idx] = tmp_val + h
        fxh1 = f()

        # f(x - h)の算出
        x[idx] = tmp_val - h
        fxh2 = f()

        # 勾配を算出
        grad[idx] = (fxh1 - fxh2) / (2 * h)

        x[idx] = tmp_val  # 値を戻す
        it.iternext()

    return grad


def cross_entropy_error(y, t):
    """交差エントロピー誤差の算出

    Args:
        y (numpy.ndarray): ニューラルネットワークの出力
        t (numpy.ndarray): 正解のラベル

    Returns:
        float: 交差エントロピー誤差
    """

    # データ1つ場合は形状を整形
    if y.ndim == 1:
        t = t.reshape(1, t.size)
        y = y.reshape(1, y.size)

    # 誤差を算出してバッチ数で正規化
    batch_size = y.shape[0]
    return -np.sum(t * np.log(y + 1e-7)) / batch_size


def sigmoid_grad(x):
    """5章で学ぶ関数。誤差逆伝播法を使う際に必要。
    """
    return (1.0 - sigmoid(x)) * sigmoid(x)

softmaxは、ゼロから作るDeep Learningで素人がつまずいたことメモ:3章で直したものをさらにスッキリさせてみました。この本のGitHubリポジトリのissueにあったsoftmax関数のコード改善案 #45を参考にしています。

numerical_gradientは前述のように、引数fで渡す関数の引数をなくしました。また、多次元配列に対応させるため、numpy.nditerでループしています。なお、本のコードではnumpy.nditerを使う際にop_flags=['readwrite']を指定していますが、xにアクセスするためのインデックスをmulti_indexで取り出しているだけで、イテレータにより列挙させたオブジェクトを更新している訳ではないのでop_flagsは省略（op_flags=['readonly']になる）しました。詳細は英語ですがIterating Over Arrays#Modifying Array Valuesを参照してください。

最後の関数sigmoid_gradは5章で学ぶものなのですが、処理時間短縮のために必要なので（後述）、本の通り実装しています。

続いて2層ニューラルネットワークを実装したtwo_layer_net.pyです。

two_layer_net.py

# coding: utf-8
from functions import sigmoid, softmax, numerical_gradient, \
    cross_entropy_error, sigmoid_grad
import numpy as np


class TwoLayerNet:

    def __init__(self, input_size, hidden_size, output_size,
                 weight_init_std=0.01):
        """２層のニューラルネットワーク

        Args:
            input_size (int): 入力層のニューロンの数
            hidden_size (int): 隠れ層のニューロンの数
            output_size (int): 出力層のニューロンの数
            weight_init_std (float, optional): 重みの初期値の調整パラメーター。デフォルトは0.01。
        """

        # 重みの初期化
        self.params = {}
        self.params['W1'] = weight_init_std * \
            np.random.randn(input_size, hidden_size)
        self.params['b1'] = np.zeros(hidden_size)
        self.params['W2'] = weight_init_std * \
            np.random.randn(hidden_size, output_size)
        self.params['b2'] = np.zeros(output_size)

    def predict(self, x):
        """ニューラルネットワークによる推論

        Args:
            x (numpy.ndarray): ニューラルネットワークへの入力

        Returns:
            numpy.ndarray: ニューラルネットワークの出力
        """
        # パラメーター取り出し
        W1, W2 = self.params['W1'], self.params['W2']
        b1, b2 = self.params['b1'], self.params['b2']

        # ニューラルネットワークの計算（forward）
        a1 = np.dot(x, W1) + b1
        z1 = sigmoid(a1)

        a2 = np.dot(z1, W2) + b2
        y = softmax(a2)

        return y

    def loss(self, x, t):
        """損失関数の値算出

        Args:
            x (numpy.ndarray): ニューラルネットワークへの入力
            t (numpy.ndarray): 正解のラベル

        Returns:
            float: 損失関数の値
        """
        # 推論
        y = self.predict(x)

        # 交差エントロピー誤差の算出
        loss = cross_entropy_error(y, t)

        return loss

    def accuracy(self, x, t):
        """認識精度算出

        Args:
            x (numpy.ndarray): ニューラルネットワークへの入力
            t (numpy.ndarray): 正解のラベル

        Returns:
            float: 認識精度
        """
        y = self.predict(x)
        y = np.argmax(y, axis=1)
        t = np.argmax(t, axis=1)

        accuracy = np.sum(y == t) / x.shape[0]
        return accuracy

    def numerical_gradient(self, x, t):
        """重みパラメーターに対する勾配の算出

        Args:
            x (numpy.ndarray): ニューラルネットワークへの入力
            t (numpy.ndarray): 正解のラベル

        Returns:
            dictionary: 勾配を格納した辞書
        """
        grads = {}
        grads['W1'] = \
            numerical_gradient(lambda: self.loss(x, t), self.params['W1'])
        grads['b1'] = \
            numerical_gradient(lambda: self.loss(x, t), self.params['b1'])
        grads['W2'] = \
            numerical_gradient(lambda: self.loss(x, t), self.params['W2'])
        grads['b2'] = \
            numerical_gradient(lambda: self.loss(x, t), self.params['b2'])

        return grads

    def gradient(self, x, t):
        """5章で学ぶ関数。誤差逆伝播法の実装
        """
        W1, W2 = self.params['W1'], self.params['W2']
        b1, b2 = self.params['b1'], self.params['b2']
        grads = {}

        batch_num = x.shape[0]

        # forward
        a1 = np.dot(x, W1) + b1
        z1 = sigmoid(a1)
        a2 = np.dot(z1, W2) + b2
        y = softmax(a2)

        # backward
        dy = (y - t) / batch_num
        grads['W2'] = np.dot(z1.T, dy)
        grads['b2'] = np.sum(dy, axis=0)

        dz1 = np.dot(dy, W2.T)
        da1 = sigmoid_grad(a1) * dz1
        grads['W1'] = np.dot(x.T, da1)
        grads['b1'] = np.sum(da1, axis=0)

        return grads

ほとんど本のコードと同じです。最後のgradientは5章で学ぶものなのですが、処理時間短縮のために必要なので（後述）、本の通り実装しています。

最後にミニバッチ学習の実装です。

mnist.py

# coding: utf-8
import numpy as np
import matplotlib.pylab as plt
import os
import sys
from two_layer_net import TwoLayerNet
sys.path.append(os.pardir)  # パスに親ディレクトリ追加
from dataset.mnist import load_mnist


# MNISTの訓練データとテストデータ読み込み
(x_train, t_train), (x_test, t_test) = \
    load_mnist(normalize=True, one_hot_label=True)

# ハイパーパラメーター設定
iters_num = 10000       # 更新回数
batch_size = 100        # バッチサイズ
learning_rate = 0.1     # 学習率

# 結果の記録リスト
train_loss_list = []    # 損失関数の値の推移
train_acc_list = []     # 訓練データに対する認識精度
test_acc_list = []      # テストデータに対する認識精度

train_size = x_train.shape[0]  # 訓練データのサイズ
iter_per_epoch = max(train_size / batch_size, 1)    # 1エポック当たりの繰り返し数

# 2層のニューラルワーク生成
network = TwoLayerNet(input_size=784, hidden_size=50, output_size=10)

# 学習開始
for i in range(iters_num):

    # ミニバッチ生成
    batch_mask = np.random.choice(train_size, batch_size, replace=False)
    x_batch = x_train[batch_mask]
    t_batch = t_train[batch_mask]

    # 勾配の計算
    # grad = network.numerical_gradient(x_batch, t_batch)  遅いので誤差逆伝搬法で……
    grad = network.gradient(x_batch, t_batch)

    # 重みパラメーター更新
    for key in ('W1', 'b1', 'W2', 'b2'):
        network.params[key] -= learning_rate * grad[key]

    # 損失関数の値算出
    loss = network.loss(x_batch, t_batch)
    train_loss_list.append(loss)

    # 1エポックごとに認識精度算出
    if i % iter_per_epoch == 0:
        train_acc = network.accuracy(x_train, t_train)
        test_acc = network.accuracy(x_test, t_test)
        train_acc_list.append(train_acc)
        test_acc_list.append(test_acc)

        # 経過表示
        print(f"[更新数]{i: >4} [損失関数の値]{loss:.4f} "
              f"[訓練データの認識精度]{train_acc:.4f} [テストデータの認識精度]{test_acc:.4f}")

# 損失関数の値の推移を描画
x = np.arange(len(train_loss_list))
plt.plot(x, train_loss_list, label='loss')
plt.xlabel("iteration")
plt.ylabel("loss")
plt.xlim(left=0)
plt.ylim(bottom=0)
plt.show()

# 訓練データとテストデータの認識精度の推移を描画
x2 = np.arange(len(train_acc_list))
plt.plot(x2, train_acc_list, label='train acc')
plt.plot(x2, test_acc_list, label='test acc', linestyle='--')
plt.xlabel("epochs")
plt.ylabel("accuracy")
plt.xlim(left=0)
plt.ylim(0, 1.0)
plt.legend(loc='lower right')
plt.show()

本のコードでは、ミニバッチ生成の際に使っているnumpy.random.choiceの引数にreplace=Falseの指定がありませんが、これだと同じ要素を重複して取り出してしまうことがありそうなので指定してみました。

勾配の算出は、本来はTwoLayerNet.numerical_gradientを使って数値微分でやるのですが、処理速度が遅くて手元の環境では1日かかっても10,000回の更新が終わらなそうです。そのため本のアドバイスに従い、5章で学ぶ誤差伝搬法を実装したTwoLayerNet.gradientを使いました。

最後に損失関数の値の推移と、訓練データ・テストデータの認識精度の推移をグラフで表示しています。

以下、実行結果です。

[更新数]   0 [損失関数の値]2.2882 [訓練データの認識精度]0.1044 [テストデータの認識精度]0.1028
[更新数] 600 [損失関数の値]0.8353 [訓練データの認識精度]0.7753 [テストデータの認識精度]0.7818
[更新数]1200 [損失関数の値]0.4573 [訓練データの認識精度]0.8744 [テストデータの認識精度]0.8778
[更新数]1800 [損失関数の値]0.4273 [訓練データの認識精度]0.8972 [テストデータの認識精度]0.9010
[更新数]2400 [損失関数の値]0.3654 [訓練データの認識精度]0.9076 [テストデータの認識精度]0.9098
[更新数]3000 [損失関数の値]0.2816 [訓練データの認識精度]0.9142 [テストデータの認識精度]0.9146
[更新数]3600 [損失関数の値]0.3238 [訓練データの認識精度]0.9195 [テストデータの認識精度]0.9218
[更新数]4200 [損失関数の値]0.2017 [訓練データの認識精度]0.9231 [テストデータの認識精度]0.9253
[更新数]4800 [損失関数の値]0.1910 [訓練データの認識精度]0.9266 [テストデータの認識精度]0.9289
[更新数]5400 [損失関数の値]0.1528 [訓練データの認識精度]0.9306 [テストデータの認識精度]0.9320
[更新数]6000 [損失関数の値]0.1827 [訓練データの認識精度]0.9338 [テストデータの認識精度]0.9347
[更新数]6600 [損失関数の値]0.1208 [訓練データの認識精度]0.9362 [テストデータの認識精度]0.9375
[更新数]7200 [損失関数の値]0.1665 [訓練データの認識精度]0.9391 [テストデータの認識精度]0.9377
[更新数]7800 [損失関数の値]0.1787 [訓練データの認識精度]0.9409 [テストデータの認識精度]0.9413
[更新数]8400 [損失関数の値]0.1564 [訓練データの認識精度]0.9431 [テストデータの認識精度]0.9429
[更新数]9000 [損失関数の値]0.2361 [訓練データの認識精度]0.9449 [テストデータの認識精度]0.9437
[更新数]9600 [損失関数の値]0.2183 [訓練データの認識精度]0.9456 [テストデータの認識精度]0.9448

結果を見ると、すでに認識精度が94.5%くらいになっていて、3章で用意されていた学習済みパラメーターの認識精度を超えていました。

4.6 まとめ

4章は本として読むだけなら良いのかも知れませんが、実装しながら進めると結構大変でした。
（ソフトマックス関数と数値微分の関数を多次元配列に対応させる部分の解説は欲しかったなぁ……）

この章は以上です。誤りなどありましたら、ご指摘いただけますとうれしいです。

はじめに

この記事は個人的な勉強メモです。inputしたものはoutputしなくてはという強迫観念に駆られて記事を書いています。

あわよくば詳しい人に誤りの指摘やアドバイスを頂ければいいなという思いを込めてQiitaの記事にしています。

エンジニアとして社会人生活を送っていますが、デザインパターンについてちゃんと学んだことがなかったので勉強してみました。

ここに記載している内容は
https://github.com/ck-fm0211/notes_desigh_pattern
にuploadしています。

過去ログ

デザインパターンについて勉強してみた（個人的メモ）その１
デザインパターンについて勉強してみた（個人的メモ）その２
デザインパターンについて勉強してみた（個人的メモ）その３
デザインパターンについて勉強してみた（個人的メモ）その４

Compositeパターン

• Composite パターンは、「容器と中身を同一視する」ことで、再帰的な構造の取り扱いを容易にするもの
• 例：ファイルシステム
  • あるフォルダ以下のファイルやフォルダをすべて削除したい場合など、それがファイルなのかフォルダなのかを意識せずに、同じように削除できたほうが都合が良い

実際に使ってみる

題材

• サンプルケースでは、ディレクトリとファイルを考える。
• Composite パターンを意識せずに、ファイルとディレクトリを表すクラスを作成してみる。

class File:
    def __init__(self, name):
        self._name = name

    def remove(self):
        print("{}を削除しました".format(self._name))

• ディレクトリを表す、Directory クラスは、List オブジェクトとして、配下のディレクトリとファイルのオブジェクトを管理し、remove メソッドが呼ばれた場合には、list に保持しているオブジェクトをすべて削除してから、自らを削除するものとする

class Directory:
    def __init__(self, name):
        self._name = name
        self._list = []

    def add(self, arg):
        self._list.append(arg)

    def remove(self):
        itr = iter(self._list)

        i = 0

        while next(itr, None) is not None:

            obj = self._list[i]

            if isinstance(obj, File):
                obj.remove()
            elif isinstance(obj, Directory):
                obj.remove()
            else:
                print("削除できません")

            i += 1

        print("{}を削除しました".format(self._name))


if __name__ == "__main__":
    file1 = File("file1")
    file2 = File("file2")
    file3 = File("file3")
    file4 = File("file4")

    dir1 = Directory("dir1")
    dir1.add(file1)

    dir2 = Directory("dir2")
    dir2.add(file2)
    dir2.add(file3)

    dir1.add(dir2)
    dir1.add(file4)

    dir1.remove()

• ここまでは問題ない。ここに「ディレクトリには、ディレクトリとファイルだけでなくシンボリックリンクも入るようにしたい」という要求が出てくると面倒になる
• Composite パターンでは、容器の中身と入れ物を同一視する。同一視するために、容器と中身が共通のインタフェースを実装するようにする。
• File と Directory が共通のインタフェース DirectoryEntry を実装するようにする

class DirectoryEntry(metaclass=ABCMeta):
    @abstractmethod
    def remove(self):
        pass

• DirectoryEntry インタフェースでは、remove メソッドのみを定義する
• これを実装する形でFileクラス、Directoryクラスを実装する。

class File(DirectoryEntry):
    def __init__(self, name):
        self._name = name

    def remove(self):
        print("{}を削除しました".format(self._name))


class Directory(DirectoryEntry):
    def __init__(self, name):
        self._name = name
        self._list = []

    def add(self, entry: DirectoryEntry):
        self._list.append(entry)

    def remove(self):
        itr = iter(self._list)

        i = 0

        while next(itr, None) is not None:

            obj = self._list[i]
            obj.remove()

            i += 1

        print("{}を削除しました".format(self._name))

• Directory クラス、File クラスを共に DirectoryEntry クラスを実装するクラスとすることで、 Directory クラスの remove メソッド内では、実態が File クラスのインスタンスであるのか、Directory クラスのインスタンスであるのかを気にせず、どちらも DirectoryEntry オブジェクトとして扱うことができるようになっている。
• このように Composite パターンを利用していることで、SymbolicLink クラスを追加する必要が生じた場合も、 柔軟に対応できる。
  • DirectoryEntry インタフェースを 実装するように、SymbolicLink クラスを実装すればよい。

class SymbolicLink(DirectoryEntry):
    def __init__(self, name):
        self._name = name

    def remove(self):
        print("{}を削除しました".format(self._name))

Compositeパターンのまとめ

Decoratorパターン

• Decorator パターンでは、飾り枠と中身を同一視することで、より柔軟な機能拡張方法を提供する。
• Decoratorパターンは機能を一つひとつかぶせていくイメージ。ある機能を持ったDecorationをコアとなるものにかぶせていくイメージである。

実際に使ってみる

題材

• アイスクリーム屋では、自由にトッピングを選べるようになっている。客は、トッピングしなくても良いし、複数のトッピングを重ねて選択することもできる。
• アイスクリーム共通のインタフェースとして、以下のインタフェースを定義する。

class Icecream(metaclass=ABCMeta):

    @abstractmethod
    def get_name(self):
        pass

    @abstractmethod
    def how_sweet(self):
        pass

• これらのインタフェースを持つクラスとしては、バニラアイスクリームクラス、抹茶アイスクリームクラスなどが以下のように提供されている。

class VanillaIcecream(Icecream):
    def get_name(self):
        return "バニラアイスクリーム"

    def how_sweet(self):
        return "バニラ味"


class GreenTeaIcecream(Icecream):
    def get_name(self):
        return "抹茶アイスクリーム"

    def how_sweet(self):
        return "抹茶味"

• これらのアイスクリームインタフェース実装クラスにトッピングをしていくことを考える。
• トッピングとしては、カシューナッツ、スライスアーモンドを考えてみる。
  • カシューナッツがトッピングされたバニラアイスクリームや、スライスアーモンドがトッピングされたバニラアイスクリームが要求される。
  • ここでは、トッピングを乗せることで、名前(getName メソッドの返り値)が変わり、味(howSweet() メソッドの返り値) は変わらないことにする。
• このような要求を満たすために、カシューナッツがトッピングされたバニラアイスクリームを表現するために、カシューナッツバニラアイスクリームクラスを作成する方法が考えられる。

class CashewNutsVanillaIcecream(Icecream):
    def get_name(self):
        return "カシューナッツバニラアイスクリーム"

• このような「継承を利用した機能の追加」は、非常に固定的なものとなってしまう。
  • 例えば、カシューナッツを乗せた抹茶アイスクリームを表すインスタンスが欲しい場合は、抹茶アイスクリーム継承クラスが必要となる。
• Decorator パターンは、このように、様々な機能追加を柔軟に行いたい場合に威力を発揮する。
• Decorator パターンを利用した設計では、拡張機能部分のみを持たせた別クラスを用意し、 そのクラスのインスタンス変数に、拡張対象となるインスタンスを持たせ、 拡張対象と同じインタフェースを実装させる。

class CashewNutsToppingIcecream(Icecream):

    def __init__(self, ice: Icecream):
        self._ice = ice

    def get_name(self):
        name = "カシューナッツ"
        name += self._ice.get_name()
        return name

    def how_sweet(self):
        return self._ice.how_sweet()

• CashewNutsToppingIcecream クラスは、カシューナッツがトッピングされたアイスクリームを表すクラス。
  • このクラスは、Icecream インタフェースを実装し、その getName() メソッドでは、自身が持つインスタンス変数 ice(Icecream インスタンス) の getName() で得られる値に「カシューナッツ」という文字列を付加した値を返り値として返す。また、howSweet() メソッドでは、インスタンス変数 ice の howSweet() メソッドの返り値をそのまま返している。
• このような設計とすることで、以下のように、カシューナッツがトッピングされたバニラアイスクリームも、カシューナッツがトッピングされた抹茶アイスクリームも、 スライスアーモンドがトッピングされたバニラアイスや、スライスアーモンドと、カシューナッツの両方がトッピングされたバニラアイスクリームなど、多様な組合せでのトッピングが可能になる。

ice1 = CashewNutsToppingIcecream(VanillaIcecream())  # カシューナッツトッピングのバニラアイス
ice2 = CashewNutsToppingIcecream(GreenTeaIcecream())  # カシューナッツトッピングの抹茶アイス

Decoratorパターンのまとめ

Visitorパターン

• Visitor パターンでは、「処理」を訪問者である Visitor オブジェクトに記述することで、処理の追加を簡単にする。
• 処理対象となる、Acceptor オブジェクトは、Visitor オブジェクトを受け入れる accept(Visitor visitor)メソッドを実装している必要ある。
• 例
  • 家の「水道工事」を行ってもらう場合、「水道工事業者」を家に呼んで、「よろしくお願いします。」と言って、後は全てお任せする。
  • そのほかにも、電気工事業者を呼ぶことも、リフォーム業者を呼ぶこともある。
  • これらの訪問者に対して、あなたは、「では、よろしく」と言って、ほとんどの作業をお任せする
  • お任せの仕方に多少の違いがあるかもしれないが、最終的には、全てを業者にお任せすることになる。
  • もし、新しいサービスを提供する業者が現れたときにも、各家庭は、なんら態度を変える必要が無く、その業者を呼んで、「よろしくお願いします。」というだけで、その新しいサービスを受けることができる。
• Visitor パターンでは、このように、受け入れる側に処理を追加することなく、処理を追加することができるパターン。

実際に使ってみる

題材

• 家庭訪問を例に考える。
• 各家庭では、先生であろうと、近所のおばちゃんであろうと、訪問者が訪れると、知らない人でなければ、「いらっしゃい」と言って受け入れる。
• この際、各家庭を Acceptor 、先生を Visitor として、Visitor パターンに当てはめて考えてみる。
• Visitor パターンでは、Visitor は、訪問対象となる、家庭を訪問する。訪問された家庭は、「ようこそいらっしゃいました」と先生を受けいれる。
• このとき、新人だろうがベテランだろうが、先生を受け入れる側は変化がない。

# -*- coding:utf-8 -*-
from abc import ABCMeta, abstractmethod


# 先生クラス
class Teacher(metaclass=ABCMeta):

    def __init__(self, students):
        self._students = students

    @abstractmethod
    def visit(self, student_home):
        getattr(self, 'visit_' + student_home.__class__.__name__.lower())(student_home)

    @abstractmethod
    def get_student_list(self):
        return self._students


# 新人先生クラス
class RookieTeacher(Teacher):

    def __init__(self, students):
        super().__init__(students)

    def visit(self, student_home):
        print("先生：こんにちは")
        super().visit(student_home)

    @staticmethod
    def visit_tanaka(tanaka):
        tanaka.praised_child()

    @staticmethod
    def visit_suzuki(suzuki):
        suzuki.reproved_child()

    def get_student_list(self):
        return self._students


# 家庭クラス
class Home(metaclass=ABCMeta):
    @staticmethod
    def praised_child():
        pass

    @staticmethod
    def reproved_child():
        pass


# 受け入れインタフェース
class TeacherAcceptor(metaclass=ABCMeta):

    def accept(self, teacher: Teacher):
        pass


# 鈴木さんの家庭
class Suzuki(Home, TeacherAcceptor):
    @staticmethod
    def praised_child():
        print("スズキ母：あら、先生ったらご冗談を")

    @staticmethod
    def reproved_child():
        print("スズキ母：うちの子に限ってそんなことは・・・。")

    def accept(self, teacher: Teacher):
        teacher.visit(self.__class__)


# 田中さんの家庭
class Tanaka(Home, TeacherAcceptor):
    @staticmethod
    def praised_child():
        print("タナカ母：あらあら、先生ったらご冗談を")

    @staticmethod
    def reproved_child():
        print("タナカ母：まさか、うちの子に限ってそんなことは・・・。")

    def accept(self, teacher: Teacher):
        teacher.visit(self.__class__)

if __name__ == "__main__":
    rt = RookieTeacher(["suzuki", "tanaka"])

    rt.visit(Suzuki())
    rt.visit(Tanaka())

• 各家庭のacceptメソッドが先生（visitor）のvisitメソッドを呼び出すことで、共通の処理を実現している

Visitorパターンのまとめ

所感

• もっといい例がある気がする・・・

状況

少し前まで、WSL(ubuntu18.04LTS)上でpipenvを使っていて、
WSL上のpythonが呼ばれていて全く問題なかったのに、しばらくWSL使っていなかったら気づいたら表題の状況になっていた。

windows側のアップデートやらなんやらしていたので何が原因かはわからない。

方法

アプリと機能>アプリ実行エイリアスでpythonをオフにする

蛇足

pipenvを入れればいいとか、pyenvを入れればいいとか混乱したので、一回全部消して下記に落ち着いた。¹

pyenv install

$ git clone https://github.com/pyenv/pyenv.git ~/.pyenv
$ echo 'export PYENV_ROOT="$HOME/.pyenv"' >> ~/.bashrc
$ echo 'export PATH="$PYENV_ROOT/bin:$PATH"' >> ~/.bashrc
$ echo 'eval "$(pyenv init -)"' >> ~/.bashrc
$ source ~/.bashrc

お好みのpython

$ pyenv install 3.8.1
$ python -V
Python 3.8.1

pipenv install

$ pip list
Package    Version
---------- -------
pip        19.2.3
setuptools 41.2.0
WARNING: You are using pip version 19.2.3, however version 19.3.1 is available.
You should consider upgrading via the 'pip install --upgrade pip' command.
$ pip install --upgrade pip
$ pip install pipenv
$ pip list
Package          Version
---------------- ----------
certifi          2019.11.28
pip              19.3.1
pipenv           2018.11.26
setuptools       41.2.0
virtualenv       16.7.9
virtualenv-clone 0.5.3

pipenvの環境作成

プロジェクトのディレクトリでpipenv installでほしいバージョンを指定すると
pyenvですでにインストールしていない場合、インストールするか聞いてくる。

$ cd your_project
$ pipenv install --python 3.7
Warning: Python 3.7 was not found on your system…
Would you like us to install CPython 3.7.6 with pyenv? [Y/n]:y

PipfileとPipfile.lockが生成されている。

$ pipenv shell
(your_project) $ python -V
Python 3.7.6
(your_project) $ pyenv versions
  system
  3.7.6
* 3.8.1 (set by PYENV_VERSION environment variable)

無事イケそうです。

参考

https://qiita.com/mashita1023/items/10239f5621ef2fc8acb9
https://qiita.com/foewhoew32f320/items/bfa90ae1003e45cefe33
https://github.com/pypa/pipenv/issues/3488
https://github.com/pyenv/pyenv/tree/4e0ba2f47d97acd1284439cff57af059ce376b9d#installation

---

1. ちゃんと消えてるかはしりません ↩

はじめに

みなさん Vtuber 好きですか?
バ美肉したくありませんか?

私はしたいです.
しかし私には 3D モデリングの技術も, 絵を描く技術もありません.

ならば技術のある人に依頼すればいいじゃないか.

• 2D モデル　数万円　　!
• 3D モデル　数十万円　!!

技術もありませんが, お金はもっとありません.

安心してください.
人類の技術の進歩によって, 技術とお金, 両方なくても バ美肉できます.

これがバ美肉ですか pic.twitter.com/zNFuW3xV5U

— ナス科ナス属トマト (@tomato08131491) January 2, 2020

Talking Head Anime from a Single Image

Google Japan のソフトウェアエンジニアである Pramook Khungurn によって開発された, 1枚の正面を向いているキャラクターの画像から顔アニメーションを自動生成する技術です.

Talking Head Anime from a Single Image

ご想像の通り, ディープラーニングです.

学習データとして
MikuMikuDance

ネットワークアーキテクチャとして
Pumarola et al.'s algorithm
Zhou et al.'s view synthesis algorithm

が利用されています.

この技術によってこんなことや,

こんなことや

一枚のキャラクター画像だけで簡単な顔アニメーションを制作するシステムを作ってみました【日本語版】 https://t.co/Nv5pz9HmoQ @YouTubeさんから

— ナス科ナス属トマト (@tomato08131491) January 3, 2020

こんなことができます.

こんなことが

Talking Head Anime from a Single Image https://t.co/LsWaguCG3t @YouTubeさんから

— ナス科ナス属トマト (@tomato08131491) January 3, 2020

準備 applications

このリンクからリポジトリをクローン.

このリンクからダウンロードして解凍,
combiner.pt, face_morpher.pt, two_algo_face_rotator.pt
上記3つのモデルを data ディレクトリに保存.

このリンクからshape_predictor_68_face_landmarks.dat.bz2をダウンロードして解凍, data ディレクトリに保存.

talking-head-anime-demo-master は以下のような構成になります.

+ data
  + illust
     - placeholder.txt
     - waifu_00_256.png
     - waifu_01_256.png
      - waifu_02_256.png
     - waifu_03_256.png
     - waifu_04_256.png
  - combiner.pt
  - face_morpher.pt
  - placeholder.txt
  - shape_predictor_68_face_landmarks.dat
  - two_algo_face_rotator.pt

準備 python

別の記事で後日解説

実際に動かしてみる

ディレクトリをtalking-head-anime-demo-master に変更.

• 手動で動かす場合は, python app/manual_poser.py を実行.
• 顔をトラッキングして動かす場合はwebカメラなどを用意して, python app/puppeteer.pyを実行.

Load Image ... で画像を選択.

• 手動で動かす場合は, バーをいろいろ動かしてみると目を閉じたり, 頭を動かしたりできる.
• 顔をトラッキングして動かす場合は, 自分の頭を動かしたり, 瞬きしたりするとそれに合わせてキャラクターも動いてくれる.

あとは自由に遊ぼう!

まとめ

今回, 自分がバ美肉するのにかかった費用は, WEBカメラlogicool C270N HD WEBCAMの購入代金である
￥2010だけです. このWEBカメラでさえスマホアプリで代用可能なので, もはやバ美肉するのに金銭的ハードルはないに等しいでしょう.

プログラムを少し変えることによってアニメーション画像の表示だけにもできるので, それも後日解説したいと思っています.

それではみなさん Let's バ美肉 !

年末年始の時間を利用してPyQUBOの理解を深めるためにNクイーン問題をやってみました。

PyQUBOは、アニーリングマシンを使う上で、とても便利なQUBO作成ツールだと思います。
「n個のうちmだけを１にする問題をPyQUBOを使ってアニーリングマシンで解く」
https://qiita.com/morimorijap/items/196e7fc86ecff927bf40
などでも使っていますが、改めて、Nクイーン問題をPyQUBOを使って解いてみようと思います。

Nクイーン問題とは、Wikipediaの８クイーンは、マスが８x８個で８クイーン　https://ja.wikipedia.org/wiki/エイト・クイーン
に対して、マスがNxN個の問題を言います。

クイーンの動きは、上下左右斜めの8方向に、遮る物がない限り進めるので、マスの列、行、斜めにクイーンがいないように配置する問題となります。
これをアニーリング マシンで解くためにコスト関数の形式にすると、

N-Queen問題のハミルトニアン、　列column、行row、斜めdiagonal

$$H = \sum_{col}\left(\sum_{i \in col} x_i - 1\right)^2
+ \sum_{row}\left(\sum_{i \in row} x_i - 1\right)^2
+ \sum_{diag} \left((\sum_{i \in diag}
x_i ) (\sum_{i \in diag} x_i - 1 )\right)$$

となります。

斜めのコスト関数の表現が異なりますが、東北大学のT-Waveの記事になる「N-クイーン問題 をD-Waveマシンで解く」https://qard.is.tohoku.ac.jp/T-Wave/?p=884
も参考になります。

４x４の場合は、

のようなマスの表現を考え、PyQUBOで0,1のQUBO形式　Binary(q[0]) のようなArrayを表現します。

from pyqubo import Array, Placeholder, solve_qubo, Sum
from pprint import pprint
import numpy as np

まずは、PyQUBO など必要なものをインポートします。
4x4のマスを作るとします。

# N-Queenのマスの数
N = 4

#pyQUBOで0、１の変数をNxN正方マス分つくる
q = Array.create('q', (N*N), 'BINARY')

q_shape = np.reshape(q,(N,N))
print(q_shape)

そうすると、

[[Binary(q[0]) Binary(q[1]) Binary(q[2]) Binary(q[3])]
[Binary(q[4]) Binary(q[5]) Binary(q[6]) Binary(q[7])]
[Binary(q[8]) Binary(q[9]) Binary(q[10]) Binary(q[11])]
[Binary(q[12]) Binary(q[13]) Binary(q[14]) Binary(q[15])]]

このようにPyQUBOのBinaryを４x４配列の形式に作成することができます。
行の部分は、数式で以下のようになっていますので、

$$\sum_{row}\left(\sum_{i \in row} x_i - 1\right)^2 $$

これをPyQUBOで表すと、行ごとスライスしてそのまま足し２乗し、最後に合計するように表せます。

#行
row_const = 0.0
for row in range(N):
    row_const += (sum(i for i in q_shape[row, :]) -1)**2

また、列もどうように以下のように表せます。

$$ \sum_{col}\left(\sum_{i \in col} x_i - 1\right)^2 $$

#列
col_const = 0.0
for col in range(N):
    col_const += (sum(i for i in q_shape[:, col]) -1)**2

となります。

そして、問題の斜めの表現ですが、

$$\sum_{diag} \left((\sum_{i \in diag}
x_i ) (\sum_{i \in diag} x_i - 1 )\right)$$

をPyQUBOで表すことになりますが、Numpyの変換を利用し、「＼」「/」の２つに分けて記述することができます。（もっといい方法をご存知の方は教えてください）

#斜め　
# ＼の方
diag_const = 0.0
xi = 0.0
xi_1 = 0.0
for k in range(-N+1,N):
    xi += (sum(i for i in np.diag(q_shape,k=k)))
    xi_1 += (sum(i for i in np.diag(q_shape,k=k))) -1 
    diag_const += xi * xi_1

と

#斜め
# /の方のために入れ替えて, ＼ をする。
diag_const_f = 0.0
xi = 0.0
xi_1 = 0.0
for k in range(-N+1,N):
    xi += (sum(i for i in np.diag(np.fliplr(q_shape),k=k)))
    xi_1 += (sum(i for i in np.diag(np.fliplr(q_shape),k=k))) -1 
    diag_const_f += xi * xi_1

となります。

これを全て足すとことで、エネルギー関数のハミルトニアンを表現することができます。
また、パラメータ調整のためにPyQUBOのPlaceholderを使いalpha,beta,gammaをつけておきます。

# エネルギー (ハミルトニアン) を構築

alpha = Placeholder("alpha")
beta = Placeholder("beta")
gamma = Placeholder("gamma")

#パラメータ
feed_dict = {'alpha': 1.0, 'beta': 1.0, 'gamma': 1.0}

H = alpha * col_const + beta * row_const  + gamma *(diag_const + diag_const_f)

次にQUBOへ変換をします。

#QUBOのコンパイル
model = H.compile()
qubo, offset = model.to_qubo(feed_dict=feed_dict)

pprint(qubo)

だけで、QUBOを作ってくれます。

{('q[0]', 'q[0]'): -11.0,
('q[0]', 'q[10]'): 6.0,
('q[0]', 'q[11]'): 4.0,
('q[0]', 'q[12]'): 2.0,
('q[0]', 'q[13]'): 6.0,
('q[0]', 'q[14]'): 6.0,
('q[0]', 'q[15]'): 6.0,
('q[0]', 'q[1]'): 6.0,
('q[0]', 'q[2]'): 4.0,

・・・長いので、途中省略・・・

('q[8]', 'q[8]'): -19.0,
('q[8]', 'q[9]'): 14.0,
('q[9]', 'q[9]'): -21.0}

です。　これをそのまま、PyQUBOのSAで計算することができます。

# PyQUBOのSAで計算
solution = solve_qubo(qubo)

# pyQUBOのSAで計算結果のデコード
decoded_solution, broken, energy = model.decode_solution(solution, vartype="BINARY", 
feed_dict=feed_dict)

pprint(solution)

結果は、

{'q[0]': 0,
'q[10]': 0,
'q[11]': 0,
'q[12]': 0,
'q[13]': 0,
'q[14]': 1,
'q[15]': 0,
'q[1]': 1,
'q[2]': 0,
'q[3]': 0,
'q[4]': 0,
'q[5]': 0,
'q[6]': 0,
'q[7]': 1,
'q[8]': 1,
'q[9]': 0}

となり、一応条件は満たしていますが、13にQueenが入ってもいいので、パラメータの調整などが必要な感じがします。
追記：斜めの一番角が抜けていましたので、最適解が出てなかったようです。
今回はちゃんと解が出るようになりました。

番外：
レゴのハリーポッターのアドベントカレンダー2019についていたチェス盤で実験（笑）

PyQUBOについては、公式ドキュメント参考になります。
https://pyqubo.readthedocs.io/

以上。

追記：
2020/1/3 9:20pm
斜めの条件が間違っていたので編集しました。

注意：型に詳しくない為用語など間違っている可能性があります。まさかり投げてください。

ソースリーディングをしていて以下のようなコードに出くわした時

class BaseRequestHandler:

    """Base class for request handler classes.
    This class is instantiated for each request to be handled.  The
    constructor sets the instance variables request, client_address
    and server, and then calls the handle() method.  To implement a
    specific service, all you need to do is to derive a class which
    defines a handle() method.
    The handle() method can find the request as self.request, the
    client address as self.client_address, and the server (in case it
    needs access to per-server information) as self.server.  Since a
    separate instance is created for each request, the handle() method
    can define other arbitrary instance variables.
    """

    def __init__(self, request, client_address, server):
        self.request = request
        self.client_address = client_address
        self.server = server
        self.setup()
        try:
            self.handle()
        finally:
            self.finish()

request, client_address, serverってなんやねんってなる。
このコードを書いた人はある特定の型のオブジェクトがここに突っ込まれることを想定していると思う。でも型が明示されていない為それが何なのか分からない。呼び出す側を探してそれが何なのか確認しなきゃならない。
コードを書いている人の脳内にある型の情報がソースコード上では省略されているからこんなことになる。

はじめに

本記事は、Apache Beam Documentation の内容を翻訳したものをベースとしています。

Apache Beam Python SDK でバッチ処理が可能なプログラムを実装し、Cloud Dataflow で実行する手順や方法をまとめています。また、Apache Beam の基本概念やテストなどについても少し触れています。

Beam SDK

Beam SDK は、Java, Python, Go の中から選択することができ、以下のような分散処理の仕組みを単純化する機能を提供しています。

• Pipeline：処理タスク全体（パイプライン）をカプセル化します。処理タスクには、入力データの読み取り、変換処理、および出力データの書き込み等が含まれます。
• PCollection：分散処理対象のデータセットを表すオブジェクトです。通常は、外部のデータソースからデータを読み取り、PCollectionを作成しますが、インメモリから作成することも可能です。
• PTransform：データ変換処理の機能を提供します。すべてのPTransformは、1つ以上のPCollectionを入力として受け取り、そのPCollectionの要素に対して何らかの処理を実行して、0個以上のPCollectionを出力します。
• I/O Transform：様々な外部ストレージシステム（GCS や BigQuery など）で、データの読み書きができる機能を提供しています。

実行環境

Beam SDK によって作成されたプログラムは、以下のような分散データ処理システム上で実行することができます。Apache Beam では、この実行環境のことをランナーと呼んでいます。

• DirectRunner： ローカルマシン上（テストを行う際などに使う）
• SparkRunner： Apache Spark
• DataflowRunner： Google Cloud Dataflow
• その他はこちらを参照

処理フローの例

一般的な（単純な）Beam プログラムは以下のように動作します。

1. Pipeline オブジェクトを作成し、実行オプションを設定します。
2. Read Transform を使用して外部ストレージシステムまたはインメモリからデータを読み込み、 PCollectionを作成します。
3. 各 PCollection に PTransform を適用します。PTransformは、PCollection内の要素を様々なロジックで変換処理することが可能です。
4. Write Transform を適用して、PTransform によって変換された PCollection を外部ソースに書き込みます。

この処理フローの場合、以下のような単純なパイプラインになります。

パイプラインの実装

上記のパイプラインを、実際に Python で実装してみます。

• Python バージョン：2系で 2.7 以上 または 3系で 3.5 以上
• Beam バージョン：2.15.*

Beam SDK のインストール

追加パッケージを特に必要としない場合は、次のコマンドでインストールします。

pip install apache-beam

今回は、Dataflow（GCP） 上で実行することを想定しているため、GCPの追加パッケージもインストールしておきます。

pip install apache-beam[gcp]

完成コード

こちらは完成形のコードになります。以下で各々について説明していきます。

pipeline.py

import apache_beam as beam
from apache_beam.options.pipeline_options import PipelineOptions
from apache_beam.options.pipeline_options import StandardOptions


class MyOptions(PipelineOptions):
    """カスタムオプション."""
    @classmethod
    def _add_argparse_args(cls, parser):
        parser.add_argument(
            '--input',
            default='./input.txt',
            help='Input path for the pipeline')

        parser.add_argument(
            '--output',
            default='./output.txt',
            help='Output path for the pipeline')


class ComputeWordLength(beam.DoFn):
    """文字数を求める変換処理."""

    def __init__(self):
        pass

    def process(self, element):
        yield len(element)


def run():
    options = MyOptions()
    # options.view_as(StandardOptions).runner = 'DirectRunner'
    p = beam.Pipeline(options=options)

    (p
     | 'ReadFromText' >> beam.io.ReadFromText(options.input)  # I/O Transform を適用して、オプションで指定したパスにデータを読み込む
     | 'ComputeWordLength' >> beam.ParDo(ComputeWordLength())  # Transform を適用
     | 'WriteToText' >> beam.io.WriteToText(options.output))  # I/O Transformを適用して、オプションで指定したパスにデータを書き込む

    p.run()


if __name__ == '__main__':
    run()

Pipeline

Pipeline オブジェクトは、データ処理タスクのすべてをカプセル化します。Beam プログラムは通常、PCollection の作成と PTransform の適用のために、まずは Pipeline オブジェクトを作成します。

Pipeline の作成

Beam プログラムを使用するには、最初に Beam SDK の Pipeline のインスタンスを（通常 main 関数内に）作成する必要があります。そして、Pipeline を作成するときには実行オプションを設定します。

次のコードは、Pipeline のインスタンスを作成する例です。

import apache_beam as beam
from apache_beam.options.pipeline_options import PipelineOptions


p = beam.Pipeline(options=PipelineOptions())

PipelineOptions の設定

PipelineOptions を使用して、パイプラインを実行するランナーや、選択したランナーに必要な固有のオプションなどを付与することができます。例として、プロジェクトIDやファイルの格納場所などの情報が含まれる可能性があります。

import apache_beam as beam
from apache_beam.options.pipeline_options import PipelineOptions

p = beam.Pipeline(options=PipelineOptions())

選択したランナーに必要な固有のオプションを設定するには、プログラム上で行う方法と、コマンドライン引数から行う方法の2通りがありますが、その例は、後述の Dataflow で説明します。

カスタムオプションを追加する

標準の PipelineOptions に加えてカスタムオプションを追加できます。次の例では、入力先と出力先のパスを指定するオプションを追加しています。カスタムオプションでは、ユーザーがコマンドライン引数から --help を渡したときに表示される説明やデフォルト値を指定することもできます。

PipelineOptions を継承することで、カスタムオプションを作成することができます。

class MyOptions(PipelineOptions):
    """カスタムオプション."""
    @classmethod
    def _add_argparse_args(cls, parser):
        parser.add_argument(
            '--input',  # オプション名
            default='./input.txt',  # デフォルト値
            help='Input path for the pipeline')  # 説明

        parser.add_argument(
            '--output',
            default='./output.txt',
            help='Output path for the pipeline')

作成したオプションは次のように渡します。

p = beam.Pipeline(options=MyOptions())

カスタムオプションにデフォルト値以外の値を設定するには、コマンドライン引数から次のように値を渡します。

--input=value --output==value

PCollection

PCollection は、分散処理対象のデータセットを表すオブジェクトです。Beam のパイプラインで、PTransform は入力と出力として PCollection を使用します。そのため、パイプラインでデータを処理したい場合は、PCollection を作成する必要があります。

Pipeline オブジェクトを作成したら、まず何らかの形で少なくとも1つの PCollection を作成する必要があります。

PCollection の作成

I/O Transform を使用して外部ソースからデータを読み取るか、インメモリから PCollection を作成します。後者は主にテストやデバッグする際に役立ちます。

外部ソースから PCollection を作成する

外部ソースから PCollection を作成するには、I/O Transform を使用します。データを読み取るためには、各 I/O Transform が提供する Read Transform を Pipeline オブジェクトに適用します。

PCollection を作成するために Read Transform を Pipeline に適用する方法は次のとおりです。

lines = p | 'ReadFromText' >> beam.io.ReadFromText('gs://some/input-data.txt')

インメモリから PCollection を作成する

インメモリから PCollection を作成するには、Create Transform を使用します。

lines = (p | 'ReadFromInMemory' >> beam.Create(['To be, or not to be: that is the question: ', 'Whether \'tis nobler in the mind to suffer ', 'The slings and arrows of outrageous fortune, ', 'Or to take arms against a sea of troubles, ']))

PTransform

PTransformは、一般的な処理フレームワークを提供します。PTransform は、入力の PCollection の各要素に適用されます。

Beam SDK は、PCollection に適用できる様々な PTransform を提供しています。これには、ParDo や Combine などの汎用な Core transforms や、1つ以上の Core transforms を組み合わせた Composite transforms が含まれます。様々な PTransform が提供されていますので、こちらなどを参照してみてください。

PTransform の適用

Beam SDK の各 PTransform には、パイプ演算子 | が提供されているので、そのメソッドを入力のPCollection に適用することで PTransform を適用することができます。

[Output PCollection] = [Input PCollection] | [PTransform]

次のように PTransform を連鎖してパイプラインを作成することもできます。

[Output PCollection] = ([Initial Input PCollection] 
                             | [First PTransform]
                             | [Second PTransform]
                             | [Third PTransform])

このパイプラインは、今回の実装例と同じフローなので、このような形状のパイプラインになります。

PTransform は、入力の PCollection には変更を加えずに、新しい PCollection を作成します。PTransform によって入力の PCollection に変更が加わることはありません。 PCollection は定義上不変です。そのため、同じ PCollection に複数の PTransform を適用して PCollection を分岐させることもできます。

[Output PCollection] = [Initial Input PCollection]

[Output PCollection A] = [Output PCollection] | [Transform A]
[Output PCollection B] = [Output PCollection] | [Transform B]

このパイプラインの形状は、次のようになります。

I/O Transform

パイプラインを作成するときは、多くの場合、ファイルやデータベースなどの外部ソースからデータを読み取る必要があります。同様に、パイプラインからデータを外部ストレージシステムに出力することもできます。

Beam SDKは、一般的なデータストレージタイプに対して I/O Transform を提供しています。サポートされていないデータストレージの読み書きを行いたい場合は、独自の I/O Transform を実装する必要があります。

データの読み込み

Read Transform は、外部ソースからの読み取りデータを PCollection に変換します。パイプラインを構築している間はいつでも Read Transform を使用できますが、一般的には最初に実行します。

lines = pipeline | beam.io.ReadFromText('gs://some/input-data.txt')

データの書き込み

Write Transform は、PCollection 内のデータを外部データソースに書き込みます。パイプラインの結果を出力するには、ほとんどの場合、パイプラインの最後で Write Transform を使用します。

output | beam.io.WriteToText('gs://some/output-data')

複数ファイルからの読み込み

多くの Read Transform は、glob 演算子にマッチする複数の入力ファイルからの読み込みをサポートしています。 次の例では、glob 演算子 (*) を使用して、指定された場所に接頭辞「input-」と接尾辞「.csv」があるすべての一致する入力ファイルを読み取ります。

lines = p | 'ReadFromText' >> beam.io.ReadFromText('path/to/input-*.csv')

複数ファイルへの書き込み

Write Transform はデフォルトで複数のファイルに書き込みます。その際、ファイル名がすべての出力ファイルの接頭辞として使用されます。

次の例では、複数ファイルを1つのロケーションに書き込みます。各ファイルには、接頭辞「numbers」、および接尾辞「.csv」が付与されます。

output | 'WriteToText' >> beam.io.WriteToText('/path/to/numbers', file_name_suffix='.csv')

パイプラインの実行

それでは、完成コードを使用してパイプラインを実行してみます。実行環境として、ローカルと Cloud Dataflow それぞれで実行します。

入力には次のような文字列が含まれるテキストファイルを用意します。

input.txt

good morning.
good afternoon.
good evening.

ローカルで実行

ローカルでパイプラインを実行するには、PipelineOptions にランナーとして DirectRunner を設定しますが、特に細かい設定がない限りは、ランナーを明示的に指定する必要はありません。

次のコマンドをコマンドラインから実行します。入力先と出力先のパスは環境によって書き換えてください。

python pipeline.py --input=./input.txt --output=./output.txt

今回の実装例は、単語の文字数を数えるパイプラインですので、次のような結果が出力されます。
また、beam.io.WriteToText はデフォルトでファイル名の接尾に 00000-of-00001 という文字列を付与して複数のファイルに分散して書き込みます。1つのファイルに書き込みたい場合は、 shard_name_template 引数を空にすることで可能です。

output.txt-00000-of-00001

13
15
13

Cloud Dataflow で実行

Cloud Dataflow は、GCP (Google Cloud Platfom) で提供されている、ストリームモードまたはバッチモードでデータ処理を行うフルマネージドサービスです。利用者はサーバなどインフラの運用を気にすることなく、実質無制限の容量を従量課金制で使用して、膨大な量のデータ処理を行うことができます。

Cloud Dataflow でパイプラインを実行すると、GCP プロジェクトで Compute Engine リソースと Cloud Storage リソースを使用するジョブが作成されます。Cloud Dataflow を利用するには、GCP で Dataflow API をオンにしてください。

Cloud Dataflow で完成コードを実行するには少しの修正が必要です。次のように修正します。

pipeline.py

import apache_beam as beam
from apache_beam.options.pipeline_options import GoogleCloudOptions
from apache_beam.options.pipeline_options import PipelineOptions
from apache_beam.options.pipeline_options import StandardOptions
from apache_beam.options.pipeline_options import WorkerOptions


GCP_PROJECT_ID = 'my-project-id'
GCS_BUCKET_NAME = 'gs://my-bucket-name'
JOB_NAME = 'compute-word-length'


class MyOptions(PipelineOptions):
    """カスタムオプション."""
    @classmethod
    def _add_argparse_args(cls, parser):
        parser.add_argument(
            '--input',
            default='{}/input.txt'.format(GCS_BUCKET_NAME),  # GCS に input.txt を置く
            help='Input for the pipeline')

        parser.add_argument(
            '--output',
            default='{}/output.txt'.format(GCS_BUCKET_NAME),  # GCS に出力する
            help='Output for the pipeline')


class ComputeWordLength(beam.DoFn):
    """文字数を求める変換処理."""

    def __init__(self):
        pass

    def process(self, element):
        yield len(element)


def run():
    options = MyOptions()

    # GCP オプション
    google_cloud_options = options.view_as(GoogleCloudOptions)
    google_cloud_options.project = GCP_PROJECT_ID  # プロジェクトID
    google_cloud_options.job_name = JOB_NAME  # 任意のジョブ名
    google_cloud_options.staging_location = '{}/binaries'.format(GCS_BUCKET_NAME)  # ファイルをステージングするための GCS パス
    google_cloud_options.temp_location = '{}/temp'.format(GCS_BUCKET_NAME)  # 一時ファイルの GCS パス

    # ワーカーオプション
    options.view_as(WorkerOptions).autoscaling_algorithm = 'THROUGHPUT_BASED'  # 自動スケーリングを有効化する

    # 標準オプション
    options.view_as(StandardOptions).runner = 'DataflowRunner'  # Dataflow ランナーを指定

    p = beam.Pipeline(options=options)

    (p
     | 'ReadFromText' >> beam.io.ReadFromText(options.input)
     | 'ComputeWordLength' >> beam.ParDo(ComputeWordLength())
     | 'WriteToText' >> beam.io.WriteToText(options.output, shard_name_template=""))

    p.run()
    # p.run().wait_until_finish()  # パイプラインの完了までブロックする


if __name__ == '__main__':
    run()

そのほかの Dataflow のオプションについてはこちらを参照してください。
ストリーミング実行するには、 streaming オプションを true にする必要があります。

こちらも、同様のコマンドで実行できます。

python pipeline.py --input=gs://my-project-id/input.txt --output=gs://my-project-id/output.txt

GCP から Dataflow サービスにアクセスするとパイプラインをモニタリングできます。UI はこのようになり、指定したパスに結果が出力されます。

こうした Dataflow のバッチ処理を定期実行したい場合などは、Dataflow テンプレートを利用すると便利です。詳しくは、こちらを参照してみてください。

パイプラインのテスト

多くの場合、パイプラインを Dataflow などのリモート実行をデバッグするよりも、ローカルで単体テストする方がデバッグにかかる時間と労力を大幅に節約できます。

依存関係の解決のために以下をインストールする必要があります。

pip install nose

パイプラインをテストするには、 TestPipeline オブジェクトを用います。入力は外部ソースから読み取る代わりに、apache_beam.Create を用いてインメモリから PCollection を作成します。出力結果を assert_that で比較します。

test_pipeline.py

from unittest import TestCase

import apache_beam as beam
from apache_beam.testing.test_pipeline import TestPipeline
from apache_beam.testing.util import assert_that, equal_to

from src.pipeline import ComputeWordLength


class PipelineTest(TestCase):

    def test_pipeline(self):
        expected = [
            13,
            15,
            13
        ]

        inputs = [
            'good morning.',
            'good afternoon.',
            'good evening.'
        ]

        with TestPipeline() as p:
            actual = (p
                      | beam.Create(inputs)
                      | beam.ParDo(ComputeWordLength()))

            assert_that(actual, equal_to(expected))

パイプラインの設計

上記で、既にシンプルなパイプラインと、分岐するパイプラインを作成する場合の設計（処理フロー）について簡単に説明しました。ここでは、その他の一般的なパイプラインの設計について紹介します。

複数の PCollection を生成する PTransform を持つパイプライン

Apache Beam の Additional outputs という機能を使って実現できます。

class ExtractWord(beam.DoFn):

   def process(element):
        if element.startswith('A'):
            yield pvalue.TaggedOutput('a', element)  # タグ名をつける（先頭が'A'の要素だったら'a'）
        elif element.startswith('B'):
            yield pvalue.TaggedOutput('b', element)  # タグ名をつける（先頭が'B'の要素だったら'b'）


mixed_col = db_row_col | beam.ParDo(ExtractWord()).with_outputs())

mixed_col.a | beam.ParDo(...)  # .タグ名でアクセスできる
mixed_col.b | beam.ParDo(...)

PCollection を結合する PTransform を持つパイプライン

Flatten を用いることで実現できます。

col_list = (a_col, b_col) | beam.Flatten()

複数の入力ソースを持つパイプライン

それぞれの入力ソースから PCollection を作成し、 CoGroupByKey などで Join することができます。

user_address = p | beam.io.ReadFromText(...)
user_order = p | beam.io.ReadFromText(...)

joined_col = (user_address, user_order) | beam.CoGroupByKey()

joined_col | beam.ParDo(...)

その他の便利機能

Composite transforms

Composite transforms は、複数の PTransform （ParDo, Combine, GroupByKey...） を組み合わせたものです。複数の PTransform を入れ子構造にすることで、コードがよりモジュール化されて理解しやすくなります。

実装例

Composite transforms を実装するには、PTransform クラスを継承し、expand メソッドをオーバーライドする必要があります。

"""文章の単語数を数えるパイプライン."""
import apache_beam as beam
from apache_beam.options.pipeline_options import PipelineOptions


class ComputeWordCount(beam.PTransform):
    """単語数を数える Composite transforms."""

    def __init__(self):
        pass

    def expand(self, pcoll):
        return (pcoll
                | 'SplitWithHalfSpace' >> beam.Map(lambda element: element.split(' '))
                | 'ComputeArraySize' >> beam.Map(lambda element: len(element)))


def run():
    p = beam.Pipeline(options=PipelineOptions())

    inputs = ['There is no time like the present.', 'Time is money.']

    (p
     | 'Create' >> beam.Create(inputs)
     | 'ComputeWordCount' >> ComputeWordCount()
     | 'WriteToText' >> beam.io.WriteToText('出力先のパス'))

    p.run()

if __name__ == '__main__':
    run()

output

7
3

Side inputs

Side inputs は、通常の入力（主入力）の PCollection に加えて、追加の入力（副入力）を PTransform に渡すことができる機能です。

実装例

"""全体の平均以上の長さを持つ文字列の長さを出力するパイプライン."""
import apache_beam as beam
from apache_beam.options.pipeline_options import PipelineOptions
from apache_beam import pvalue


class FilterBelowMeanLengthFn(beam.DoFn):
    """平均以下の文字数を持つ文字列をフィルタリングする."""

    def __init__(self):
        pass

    # mean_word_length は副入力
    def process(self, element, mean_word_length):
        if element >= mean_word_length:
            yield element


def run():
    p = beam.Pipeline(options=PipelineOptions())

    inputs = ["good morning.", "good afternoon.", "good evening."]

    # 主入力
    word_lengths = (p
                    | 'Create' >> beam.Create(inputs)
                    | 'ComputeWordLength' >> beam.Map(lambda element: len(element)))

    # 副入力
    mean_word_length = word_lengths | 'ComputeMeanWordLength' >> beam.CombineGlobally(beam.combiners.MeanCombineFn())

    (word_lengths
     | 'FilterBelowMeanLength' >> beam.ParDo(FilterBelowMeanLengthFn(), pvalue.AsSingleton(mean_word_length))  # 副入力を引数に取る DoFn を適用する ParDo の第２引数に副入力を挿入する（pvalue.~ は渡したいデータ型によって異なる）
     | 'write to text' >> beam.io.WriteToText('出力先のパス'))

    p.run().wait_until_finish()


if __name__ == '__main__':
    run()

「good morning.」, 「good afternoon.」, 「good evening.」の文字数はそれぞれ「13」, 「15」, 「13」で、その平均は13.67ほどなので、次のような出力になります。

output

15

パイプラインの中で何が起きているか？

「パイプラインの中で何が起きているか」について少し記述しています。

シリアライズと通信

パイプラインの分散処理において最もコストの高い操作の1つは、マシン間で要素をシリアライズして通信を行うことです。Apache Beam のランナーは、マシン間で通信を行うなどの理由で PCollection の要素をシリアライズします。次のような手法を用いて、Transform と次のステップの Transform との間で要素の通信を行います。

1. 要素をシリアライズしてワーカーにルーティングする
2. 要素をシリアライズして複数のワーカーに再分配する
3. Side inputs を使用する場合は、 要素をシリアライズしてすべてのワーカーにブロードキャストする必要がある
4. Transform と次のステップの Transform が同じワーカーで実行されている場合は、要素の通信をインメモリを使って行う（シリアライズしないことで通信コストを下げることができる）

同梱と永続化

Apache Beam は、embarassingly parallel 問題に焦点を当てています。Apache Beam は、要素を並列で処理することを重要視しているので、PCollection の各要素にシーケンス番号を割り当てるなどの動作を表現するのが苦手です。このようなアルゴリズムはスケーラビリティの問題を抱える可能性がはるかに高いためです。

すべての要素を並列に処理することにもいくつかの欠点があります。例えば、要素を出力先に書き込む場合です。出力処理において、すべての要素を並列にバッチ処理することは不可能です。

そのため、Apache Beam のランナーは、すべての要素を同時に処理するのではなく、PCollection の要素を同梱して処理します。ストリーミング処理の場合は、小さな単位で同梱して処理する傾向があり、バッチ処理の場合は、より大きな単位で同梱して処理する傾向があります。

並列処理

Transform 内の並列処理

単一の ParDo を実行する場合、Apache Beam のランナーは、PCollection の要素を2つに分割・同梱（Bundle）することがあります。

ParDoが実行されると、ワーカーは、次に示すように2つのBundleを並列で処理します。

単一の要素は分割できないため、Transform の最大並列処理は PCollection の要素数によって異なります。今回の場合の最大並列処理数は図から見て 9 です。

※ 単一の要素を複数の Bundle に分割できる機能（Splittable ParDo）が開発中らしい

Transform 間の並列処理

ParDo は従属並列になることがあります。例えば、次のように ParDo1 の出力を同じワーカーで処理する必要がある場合、ParDo1 と ParDo2 は従属並列になります。

Worker1 では、Bundle A の要素に対して ParDo1 が実行され、Bundle C になります。次に、Bundle C の要素に対して ParDo2 が実行されます。同様に、Worker2 では、Bundle B の要素に対して ParDo1 が実行され、Bundle D になります。そして、Bundle D の要素に対して ParDo2 が実行されます。

このように ParDo を実行することで、Apache Beam のランナーは、ワーカー間で要素を再配布することを回避できます。そして、これにより通信コストを節約できます。ただし、最大並列処理数は、従属並列の最初の ParDo の最大並列処理数に依存するようになります。

障害発生時の挙動

Transform 内の障害発生時の挙動

Bundle 内の要素に対する処理が失敗すると Bundle 全体が失敗してしまいます。そのため、処理を再試行する必要があります（そうしないとパイプライン全体が失敗します）。

次の例では、Worker1 が Bundle A の5つの要素すべてを正常に処理します。Worker2 は Bundle B の4つの要素を処理しますが、Bundle B の最初の2つの要素は正常に処理され、3番目の要素は処理が失敗します。

その後、Apache Beam のランナーが Bundle B のすべての要素を再試行し、2回目で処理が正常に完了しています。図のように、再試行は必ずしも元の処理の試行と同じ Worker で発生するわけではありません。

Transform 間の障害発生時の挙動

ParDo1 の処理後、ParDo2 内の要素を処理できなかった場合、これら2つの Transform は同時に失敗したことになります。

次の例では、Worker2 は Bundle B のすべての要素に対して ParDo1 を正常に実行します。しかし、Bundle D の要素を処理できないため、ParDo2 は失敗します。

その結果、Apache Beam のランナーは ParDo2 の出力を破棄して再び処理を実行する必要があります。その際には、ParDo1 の Bundle も破棄されなければならず、Bundle のすべての要素は再試行されなければいけません。

まとめ

本記事では、Apache Beam Documentation の内容をもとに学習した内容を記述してみました。
間違っている点などあればご指摘ください！