UMAPがバージョンアップしてv0.4が公開された。

2020/02/10現在では、pip install --pre umap-learnでバージョンを上げることができる。

疎行列をそのまま入力できたりいろんな機能が追加されているらしいけど、ここではプロット機能、非ユークリッド空間への埋め込み、逆変換を試してみる。

データだけ変えてほぼドキュメントに書いてあるコード例そのままやってるだけなので、それぞれについて詳しくはUMAPドキュメントへ。

データ

PARCのレポジトリに置いてあったscRNA-seqのデータセットとアノテーション（Zheng et al., 2017, 10X PBMC）を使って実験する。68,579細胞、事前にPCAで50次元に圧縮済み。気軽にやるにはちょっと大きすぎるデータなので適当に1万細胞くらいに落として使う。

import numpy as np
import pandas as pd
import umap

データ読み込み。1万細胞、50主成分。

dat = np.loadtxt('./data/pca50_pbmc10k.txt', delimiter=',')
dat.shape

(10000, 50)

アノテーションは以下のような感じ。

labels = []
for line in open('./data/zheng17_annotations_10k.txt'):
    labels.append(line.rstrip())

import collections
import pprint
pprint.PrettyPrinter(indent=4).pprint(collections.Counter(labels))

Counter({   'CD4+/CD45RA+/CD25- Naive T': 1853,
            'CD8+ Cytotoxic T': 1694,
            'CD8+/CD45RA+ Naive Cytotoxic': 1599,
            'CD4+/CD45RO+ Memory': 1425,
            'CD56+ NK': 1211,
            'CD4+/CD25 T Reg': 828,
            'CD14+ Monocyte': 595,
            'CD19+ B': 583,
            'Dendritic': 127,
            'CD4+ T Helper2': 55,
            'CD34+': 30})

とりあえずは普通にUMAPしておく。

model = umap.UMAP(verbose=True)
model.fit(dat)

プロット機能

umap.plotで新たに、umap単体でプロットができるようになった。

matplotlibで描くのとほとんど変わらないけど、手軽ではある。

umap.plotを使う場合はumapと別にdatashader, bokeh, holoviewsが必要なので注意。それぞれインストールしておく。

import umap.plot

学習させたモデルのインスタンスを与えると、scatter plotを描いてくれる。

umap.plot.points(model)

また、labelsのパラメータにそれぞれの点のラベルデータを与えると色分けして描いてくれる。

umap.plot.points(model, labels=labels)

与えるデータは連続値でもOK。ここでは適当なデータがないので、元データの平均値を与える。

背景やカラーマップなどの組み合わせを用意してくれてるので、themeパラメータで好きなものを選んで使える。

選べるthemeは次の9種類。

umap.plot._themes.keys()

dict_keys(['fire', 'viridis', 'inferno', 'blue', 'red', 'green', 'darkblue', 'darkred', 'darkgreen'])

umap.plot.points(model, 
                 values=dat.mean(axis=1),
                 theme='viridis')

また、Bokehを使ってインタラクティブなプロットを描くことも可能。

umap.plot.output_notebook()

マウスオーバーしたときに表示する情報は事前にpandas.dataframeのかたちで用意しておく。

df_labels = pd.DataFrame(labels, columns=['celltype'])
p = umap.plot.interactive(model, labels=labels, 
                          hover_data=df_labels, point_size=2)
umap.plot.show(p)

さらに、UMAP埋め込みのときに使われたneighborhood graphを可視化してくれる機能もある。エッジの重みもグラデーションで表示してくれる。どういったconnectivityが学習されたのかを検証するときに使えるかもしれない。

umap.plot.connectivity(model, show_points=True,
                      edge_cmap='viridis')

他にも、様々な指標で埋め込みを診断するためのdiagnostic plot機能があるらしい。

非ユークリッド空間への埋め込み

デフォルトではUMAPはユークリッド空間に埋め込む（低次元空間のユークリッド距離をターゲットに最適化する）が、球面など他のタイプの空間に埋め込むことも可能らしい。

これは、output_metricパラメータで低次元側の距離計算手法を指定することで可能。

まず、球面に埋め込んでみる。球面埋め込みの場合は、Haversine式を指定する。

sphere_mapper = umap.UMAP(output_metric='haversine')
sphere_mapper.fit(dat)

結果の座標は球面座標系で出てくるので、そのまま描いてもよくわからない。なので教科書通りに直交座標系に変換してからプロットしてみる。

x = np.sin(sphere_mapper.embedding_[:, 0]) * np.cos(sphere_mapper.embedding_[:, 1])
y = np.sin(sphere_mapper.embedding_[:, 0]) * np.sin(sphere_mapper.embedding_[:, 1])
z = np.cos(sphere_mapper.embedding_[:, 0])

このプロットは勝手にやってくれないみたいなので、自分でmatplotlibで描画。

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import seaborn as sns
sns.set(style='white')

categories = sorted(list(set(labels)))
label_ids = [categories.index(l) for l in labels]

fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z, c=label_ids, cmap='Spectral')

セルタイプごとにわりときれいに、球面上にまとまってくれた。

3次元でちょっとわかりづらいので、このまま2次元に展開してみる。

x = np.arctan2(x, y)
y = np.arccos(z)
fig = plt.figure(figsize=(12, 12))
plt.scatter(x, y, c=label_ids, cmap='Spectral')

左右は本来つながってる。

どういったときに使えばいいのかよくわからないけど、データが本質的に周期的な性質を持つときなどに有効？

また、双曲空間への埋め込みも紹介されている。これについてはほんとによくわからないのでドキュメントの流れそのままやってみる。ポアンカレ円板モデルそのものは最適化が難しいらしいので、Hyperboloid modelをターゲットに学習しているらしい。

hyperbolic_mapper = umap.UMAP(output_metric='hyperboloid')
hyperbolic_mapper.fit(dat)

x = hyperbolic_mapper.embedding_[:, 0]
y = hyperbolic_mapper.embedding_[:, 1]
z = np.sqrt(1 + np.sum(hyperbolic_mapper.embedding_**2, axis=1))

disk_x = x / (1 + z)
disk_y = y / (1 + z)

fig = plt.figure(figsize=(12,12))
ax = fig.add_subplot(111)
ax.scatter(disk_x, disk_y, c=label_ids, cmap='Spectral')
boundary = plt.Circle((0,0), 1, fc='none', ec='k')
ax.add_artist(boundary)
ax.axis('off');

逆変換

埋め込みが学習された低次元側の座標から、対応する高次元サンプルを生成する手法。

VAEみたいな生成モデルと同じように使うのは難しいかもしれないけど、埋め込まれた低次元がどんな空間なのかざっと確認するのには便利かも。

高次元サンプル生成についてはシングルセルのデータはあんまりおもしろくないので、ここではKuzushiji-MNISTのデータ（を適当に1万個ランダムに選んだもの）でやってみる。

dat = np.load('./data/kmnist-train-imgs_10k.npy')
dat.shape

(10000, 784)

labels = np.load('./data/kmnist-train-labels_10k.npy')
labels

array([4, 5, 0, ..., 6, 9, 0], dtype=uint8)

まずは普通にUMAPを計算してみる。

model = umap.UMAP(n_epochs=500, verbose=True).fit(dat)

umap.plot.points(model, labels=labels)

この空間から、補間しておもしろそうな領域について逆変換してみたいので、

クラスタ8の左上、クラスタ9の左下、クラスタ0の右上、クラスタ3の右下あたりの点を選んで、テスト用の点を100個作ってみる。

x = model.embedding_

top_left = x[labels == 8, :][x[labels == 8, 0].argmin()]
btm_left = x[labels == 9, :][x[labels == 9, 1].argmin()]
top_right = x[labels == 0, :][x[labels == 0, 0].argmax()]
btm_right = x[labels == 3, :][x[labels == 3, 1].argmin()]

test_pts = np.array([
    (top_left*(1-x) + top_right*x)*(1-y) +
    (btm_left*(1-x) + btm_right*x)*y
    for y in np.linspace(0, 1, 10)
    for x in np.linspace(0, 1, 10)
])

print(top_left)
print(btm_left)
print(top_right)
print(btm_right)

[-3.0056033  9.982167 ]
[2.0912035 2.2021403]
[14.147088  10.8581085]
[10.669375   2.6710818]

逆変換を実行する。inverse_transform関数に調べたい点の座標を与える。

inv_transformed_points = model.inverse_transform(test_pts)

from matplotlib.gridspec import GridSpec

fig = plt.figure(figsize=(12,6))
gs = GridSpec(10, 20, fig)
scatter_ax = fig.add_subplot(gs[:, :10])
kuzushiji_axes = np.zeros((10, 10), dtype=object)
for i in range(10):
    for j in range(10):
        kuzushiji_axes[i, j] = fig.add_subplot(gs[i, 10 + j])

scatter_ax.scatter(model.embedding_[:, 0], model.embedding_[:, 1],
                   c=labels.astype(np.int32), cmap='Spectral', s=0.1)
scatter_ax.set(xticks=[], yticks=[])

scatter_ax.scatter(test_pts[:, 0], test_pts[:, 1], marker='x', c='k', s=15)

for i in range(10):
    for j in range(10):
        kuzushiji_axes[i, j].imshow(inv_transformed_points[i*10 + j].reshape(28, 28), cmap='viridis')
        kuzushiji_axes[i, j].set(xticks=[], yticks=[])

ちょっと粗いが、軸に沿ってどのように遷移するのか、また、クラスタどうしがなぜ近い距離に配置されているのかがなんとなく類推できる。

私事で恐縮なのですが、昨年の１月頃よりポケモンGoを再開いたしました。そこから1年間続けているのですがポケモンGo側にも大きな機能追加がありました。それは待望のオンラインでの対戦です。興味が無かったのですが、やってみると思いの外楽しく自分でもシミュレーションや効率的な戦いを考えてみたくなります。
それならそのポケモンデータや技データはどうやって集める？という話です。

概要

• ポケモンGoのGAME_MASTER.jsonファイルの公開元を発見
• そのファイルから必要なデータをPythonにて抽出
• 生データであり追加の加工が必要

1.データ収集方法の検討

最初に思いついたのはこの方法でした。海外サイトや国内サイト含め、ポケモンGoにおけるステータスや技データなど全て公開されています。ただし、二次利用を想定しているものではないためウェブ上からスクレイピング等で取り出す方法が必要となります。

案1. Webサイト上からの収集

攻略サイト等で公開されているページからスクレイピングする方法です。

案2. ゲーム内のデータから抽出

"解析"という言葉で説明されているアレです。その展開済みデータが入手できれば活用出来るのではないかと思った次第です。

2.「ゲーム内のデータから抽出」を採用

Google検索等の検索力を発揮したところ展開済みデータを公開しているレポジトリ【pokemongo-dev-contrib/pokemongo-game-master】に遭遇しました。そこで一件落着、と思ったのですがデータに余計なモノが多すぎるということでした。ポケモン、技のデータ意外にもプレイヤーの衣装や各種設定値なども含んでいます。

3.ゲーム内のデータから必要な部分を探す

さてこのファイルですがGAME_MASTER.jsonという名称で前述の通りコスチュームデータからPvP対戦の設定値、ポケモンのデータ、技のエネルギー等等といったポケモンGOの設定データの多くが含まれています。そのせいで3.4MBと大きなサイズでありムダもあります。
まずはそこから必要データを抽出していきます。

色々なデータがあるのせまずは眺めてポケモンのデータ部分を見つけます。

GAME_MASTER.jsonからポケモンデータの一部抜粋・整形

{
  "itemTemplates": [{
   "templateId": "V0001_POKEMON_BULBASAUR",
    "pokemonSettings": {
      "pokemonId": "BULBASAUR",
      "type": "POKEMON_TYPE_GRASS",
      "type2": "POKEMON_TYPE_POISON",
      "stats": {
        "baseStamina": 128,
        "baseAttack": 118,
        "baseDefense": 111
      },
      "quickMoves": ["VINE_WHIP_FAST", "TACKLE_FAST"],
      "cinematicMoves": ["SLUDGE_BOMB", "SEED_BOMB", "POWER_WHIP"],
  }, 
  }],

ありましたね。なお、図中では今回重要と考えたデータ部分のみを抜粋しましたが、他にも様々なデータが含まれています。
ワザワザまとめるほどではないですがkey名やデータから推測し以下のようになると推測出来ます。

Key	データ内容
templateId	ポケモンデータの識別子
pokemonId	ポケモンの名前
stats	種族値
type	タイプ
quickMoves	技１
cinematicMoves	技２

なお、このデータには期間限定の技などは含んでいないようです。

4.プログラムでのデータ抽出

ここまでくれば私よりも皆さんの方が100倍詳しいと思いますが、一応私のプログラムも載せておきます。csv形式にて使えるように各データはカンマにて区切りました。なお技データをそのまま使うと入れ子となるのが嫌だったので、::を区切り文字として無理やりくっつけています。

GAME_MASTER_pokemon_parser.py

# -*- coding: utf-8 -*-

import json
import re

# ポケモンのデータにマッチする正規表現用パターンを用意
pattern = '^V0\d+_POKEMON_.+'

f = open('GAME_MASTER.json', 'r')
json_dict = json.load(f)

# 最上位のキーとして"itemTemplates"があるのでまずは展開する。
for json_list in json_dict["itemTemplates"]:
    templateIdData = json_list["templateId"]
    result = re.match(pattern, templateIdData)

    # ポケモンのデータにマッチした場合に必要なデータを抽出する
    if (result):
        # 技データを持たないポケモンがいるのでget形式で辞書データにアクセス
        # （ドーブルが技持って無かった）
        quickMoves = json_list["pokemonSettings"].get("quickMoves")
        quickMovesStr = ""
        if quickMoves is not None:
            for quickMovesStrTemp in quickMoves:
                quickMovesStr = quickMovesStr + quickMovesStrTemp + "::"
        cinematicMoves = json_list["pokemonSettings"].get("cinematicMoves")
        cinematicMovesStr = ""
        if cinematicMoves is not None:
            for cinematicMovesStrTemp in cinematicMoves:
                cinematicMovesStr = cinematicMovesStr + cinematicMovesStrTemp + "::"
        outPokemon = json_list.get("templateId") + "," + str(json_list["pokemonSettings"].get("type")) + "," + str(json_list["pokemonSettings"].get("type2")) + "," + str(json_list["pokemonSettings"]["stats"]["baseStamina"]) + "," + str(json_list["pokemonSettings"]["stats"]["baseAttack"]) + "," + str(json_list["pokemonSettings"]["stats"]["baseDefense"])+ ',' +quickMovesStr + ',' + cinematicMovesStr
        print(str(outPokemon))

実行すると標準出力に以下の形式にて出力されます。
ポケモンデータの識別子,ポケモンタイプ１,ポケモンタイプ２,種族値スタミナ,種族値攻撃,種族値防御,種族値技１,種族値技２

実行結果.csv（一部抜粋）

V0808_POKEMON_MELTAN,POKEMON_TYPE_STEEL,None,130,118,99,THUNDER_SHOCK_FAST::,FLASH_CANNON::THUNDERBOLT::
V0809_POKEMON_MELMETAL,POKEMON_TYPE_STEEL,None,264,226,190,THUNDER_SHOCK_FAST::,FLASH_CANNON::THUNDERBOLT::HYPER_BEAM::ROCK_SLIDE::SUPER_POWER::

5.残タスク

まずはポケモンのデータを生データとして取り出せました。ただし、このデータには以下のような問題があります。

(1) 日本語非対応
(2) 図鑑番号が無い
(3) アローラ、シャドウ、リライト、ノーマル(？)の整理
(4) 限定技が含まれていない

(2)-(4)は機械的にどうにか処理できそうです。(1)については図鑑番号とマッピングデータを用意できればどうにかなるかなぁ。
ただし、技データの場合はどうしよう…。

6.あとがき

APIに直接アクセスできた当初はサーチや個体値チェック、GO Plusモドキなどテック記事に溢れてました。ブームの収束と運営会社による対策により一気に下火になりましたが…。継続しプレイや解析を続けていたユーザーにただただ感謝しかありません。

おまけ

PvP時の技データは以下のようです。その中で厄介なのは"durationTurns"です。PvP時の技１の硬直時間[s]を表すのですが、技の硬直時間[s]=1+durationTurnsのようです。そして硬直時間が1[s]の場合にはこのkeyが無い内容でした。
VOLT_SWITCHは硬直4秒、DRAGON_BREATHは硬直1秒

    "templateId": "COMBAT_V0250_MOVE_VOLT_SWITCH_FAST",
    "combatMove": {
      "uniqueId": "VOLT_SWITCH_FAST",
      "type": "POKEMON_TYPE_ELECTRIC",
      "power": 12.0,
      "vfxName": "volt_switch_fast",
      "durationTurns": 3,
      "energyDelta": 16
    }

    "templateId": "COMBAT_V0204_MOVE_DRAGON_BREATH_FAST",
    "combatMove": {
      "uniqueId": "DRAGON_BREATH_FAST",
      "type": "POKEMON_TYPE_DRAGON",
      "power": 4.0,
      "vfxName": "dragon_breath_fast",
      "energyDelta": 3
    }

目的

iPhone上でPythonのプログラムを書いて実行できるPythonista3上でBluetoothを使ったプログラムが作れるのが分かったので、プログラミング・フォロを操作するプログラムを作ってみた。

Pythonista3

iPhone(iOS)上で動作するPythonの統合環境です。
有料で少し高いんですが、iPhoneのカメラや加速度センサーやBluetoothも扱うことも出来るし、NumPyやMatplotlibなどの一般的なライブラリも使えるので、いろいろ遊べる環境です。
なんと言ってもiPhone上で作って実行できるので、ネイティブなプログラム開発に必要なMacがなくてもアプリが作れるのが魅力です。
追加でstashをインストールすると、ssh/scpとかgitも使えます。

プログラム

プログラムはgithubに置いてます
Bluetoothの操作にはcbライブラリを利用してます。
UI上でボタンの押しっぱなしに対応するためにsceneライブラリを使用してます。
ただし今回作ったGUIではボタンの状態はOnかOffの2値しか取得出来ず中間の値を取得出来ないので、速度コントロールは出来ません。
また、マルチタッチには対応してません。
micro:bitとの接続が出来てしまえば、あとは「PS3コントローラーでプログラミング・フォロを操作する」とあまり変わらず、micro:bitのI/O端子に信号を送るだけです。

実行

実行するとこんな風に動きました。
iPhoneから直接micro:bitを操作できるので、便利で楽しいです。

関連URL

• RaspberryPiにPS3コントローラーを接続する

はじめに

最近、MapTilerの日本向けサービス開始が話題になりましたが、従来のサービスと何が違うかと言えば、ラスタータイルに加え、ベクタータイルが配信された点です。ベクタータイルが何かという説明は省きますが、昨今の地図配信で外せない要素と感じ、色々いじっています。FOSS4Gっ子ならやっぱりQGISでベクタータイル、表示したいですよね？この記事で、表示する方法（プラグイン）を紹介します。

目標

以下のように、ベクタータイルをQGISレイヤーとして読み込む。
こうなれば、欲しいレイヤーを外部ファイルに出力も出来ます、夢が広がりますね！

© MapTiler © OpenStreetMap contributors

Vector Tiles Reader

https://github.com/geometalab/Vector-Tiles-Reader-QGIS-Plugin
公式リポジトリにも登載されているので、QGISからシームレスにインストール可能です。

使い方

プリセットされているMapTiler.comなどのレイヤーを追加出来ますが、今回は任意のベクタータイルを追加する事を考えます。

追加方法

Connections下部のNewを押下すると、上画像のようなダイアログが表示されます。以下のとおり入力してSaveしましょう。

• Nameには任意の文字列
• TileJSONにはtiles.jsonのURL
• GL Style JSON URLにはstyle.jsonのURLを（Optional）

tiles.jsonのURLがわからないけどstyle.jsonはわかる場合、style.jsonの中身を見るとtiles.jsonの在り処もわかります（tiles.jsonでファイル内検索してみてください）。style.jsonはOptionalで、なくても動きますが、なんらか定義されたstyle.jsonがあるならば、追加しておいた方が当然ながら見栄えがよいです。

追加したいベクタータイルをプルダウンで選択後、「Connect」してから、下部の「Add」で、レイヤーが追加されます。もし一部のレイヤーだけQGISに追加したい場合は、Layers of〜の中から追加したいレイヤーを選択してからAddしてください（未選択の場合、全レイヤーが追加される）。以後、追加されたグループ削除するまで、画面領域（ズームレベル含む）の変更に応じ、その都度サーバーからデータを持ってきてくれます。

注意点

ベクターデータの範囲外で「Add」してはいけない

現状、ベクタータイルが提供されていない領域でAddすると、再読み込みもされない仕様となっています。
特に起動直後は表示領域がおかしな事になりがちなので、手頃なレイヤー（XYZ Tileなど）を追加して、表示領域を絞ってからAddしましょう。

複数のsourceを持つstyle.jsonのスタイリングは未対応

標記のとおりです。たとえばMapTiler.jpのJP MIERUNEのstyle.jsonは、sourcesに4つのデータを持つため本プラグインでは一切スタイリングされません。これは、style.jsonがoptionalである事からもわかるのですが、このプラグインはtiles.jsonがまずあって、「必要であれば」style.jsonでスタイリングする仕組みとなっているからです。
一方で、そもそもstyle.jsonにはtiles.jsonへのリンクが含まれているため、Mapbox GL JSでは必ずしもtiles.jsonを必要としません。
私はこれをバグだと思い、一応issueを立ててあるので、そのうち直せればなぁなんて思っています。
ちなみに、JP MIERUNEのようにstyle.jsonのsourcesに複数のデータがある場合、プラグインで指定したtiles.jsonのレイヤーのみをスタイリングするようにstyle.jsonを加工してやれば（指定しなかったtiles.jsonを対象とするデータを削除してやれば）、適切にスタイルが反映される事は確認済です（しかしこの場合でも日本語文字列は表示されないよう…これもfixしたい）。

0. 関連記事

• 自然言語処理を簡単に扱えると噂のCOTOHA APIをPythonで使ってみた
• 「メントスと囲碁の思い出」をCOTOHAさんに要約してもらった結果。COTOHA最速チュートリアル付き

n番煎じ

1. 背景

プレゼント企画 を知ったため.

個人的にテキスト解析を前々からやろうとは思っていたため, これを機に決意.
とはいえ, 調べてみると公式からの各種言語に対応したAPIは出されていなかった.(この記事 の作者さんの所属がNTTっぽいのでもしかしたらこれが公式説はあったりする)

そこで, 勉強がてら作ってみる.

COTOHA API Portalへのリンク

2. 作成したAPIの説明

2.1 動作確認した環境

requestsが必要です.

windows

• Windows 10 Home
• Python 3.7.3

Linux

• Ubuntu 18.04.3 LTS
• Python 3.6.9

Mac

♰持ってない♰

2.2 Quick start

1. まずCOTOHA API for Developersに無料登録
2. 以下のコマンドを実行

git clone https://github.com/tsuji-tomonori/cotohapy.git
cd cotohapy
pip install -r requirements.txt
# config.json の作成
python demo.py
# もし, python3 環境の場合
python3 demo.py
# action!

config.json は以下のようになります.

{
    "access_token_publish_url": "",
    "developer_api_base_url": "",
    "clientid": "",
    "clientsecret": ""
}

アカウントホームに行くと以下のページが出てくるため, それをもとに設定してください. 対応表を下につけておきます.

2.3 使用方法

GitHub

git clone https://github.com/tsuji-tomonori/cotohapy.git

インポート

from cotohapy3 import CotohaAPI

ライブラリーは「cotohapy3」の中にすべて入っています. CotohaAPI 以外にもクラスはありますが, それは後程...

アクセストークン取得

初めに必ず行ってください.

api = CotohaAPI(
    developer_api_base_url=developer_api_base_url, 
    access_token_publish_url=access_token_publish_url
)
api.login(clientid=clientid, clientsecret=clientsecret)

def __init__(self, developer_api_base_url, access_token_publish_url)

引数	データ型	内容
developer_api_base_url	str	API Base URL
access_token_publish_url	str	Access Token Publish URL

def login(self, clientid, clientsecret)

引数	データ型	内容
clientid	str	Client ID
clientsecret	str	Client secret

構文解析: parse

# 実行例
api.parse("昨日母と銀座で焼き肉を食べた")

def parse(self, sentence, **kwargs)

引数	データ型	内容
sentence	str	解析対象文
**kwargs		その他任意のパラメタ (公式リファレンスのAPI参照)

固有表現抽出: ne

# 実行例
api.ne("昨日は東京駅を利用した。")

def ne(self, sentence, **kwargs)

引数	データ型	内容
sentence	str	解析対象文
**kwargs		その他任意のパラメタ (公式リファレンスのAPI参照)

照応解析: coreference

# 実行例
api.coreference("太郎は友人です。彼は焼き肉を食べた。")

def coreference(self, document, **kwargs)

引数	データ型	内容
document	str / list	以下のどちらかの形式で指定 str: 解析対象の文 list: 解析対象の文集合
**kwargs		その他任意のパラメタ (公式リファレンスのAPI参照)

キーワード抽出: keyword

# 実行例
api.keyword("レストランで昼食を食べた。")

def keyword(self, document, **kwargs)

引数	データ型	内容
document	str / list	以下のどちらかの形式で指定 str: 解析対象の文 list: 解析対象の文集合
**kwargs		その他任意のパラメタ (公式リファレンスのAPI参照)

類似度計算: similarity

# 実行例
api.similarity("近くのレストランはどこですか？", "このあたりの定食屋はどこにありますか？")

def similarity(self, s1, s2, **kwargs)

引数	データ型	内容
s1	str	類似度算出対象のテキスト
s2	str	類似度算出対象のテキスト
**kwargs		その他任意のパラメタ (公式リファレンスのAPI参照)

文タイプ判定: sentence_type

# 実行例
api.sentence_type("あなたの名前は何ですか？")

def sentence_type(self, sentence, **kwargs)

引数	データ型	内容
sentence	str	解析対象文
**kwargs		その他任意のパラメタ (公式リファレンスのAPI参照)

ユーザー属性推定(β): user_attribute

# 実行例
api.user_attribute("私は昨日田町駅で飲みに行ったら奥さんに怒られた。")

def user_attribute(self, document, **kwargs)

引数	データ型	内容
document	str / list	以下のどちらかの形式で指定 str: 解析対象の文 list: 解析対象の文集合
**kwargs		その他任意のパラメタ (公式リファレンスのAPI参照)

言い淀み除去(β): remove_filler

# 実行例
api.remove_filler(
    "えーーっと、あの、今日の打ち合わせでしたっけ。すみません、ちょっと、急用が入ってしまって。"
)

def remove_filler(self, text, **kwargs)

引数	データ型	内容
text	str	解析対象テキスト
**kwargs		その他任意のパラメタ (公式リファレンスのAPI参照)

音声認識誤り検知(β): detect_misrecognition

# 実行例
api.detect_misrecognition("温泉認識は誤りを起こす")

def detect_misrecognition(self, sentence, **kwargs)

引数	データ型	内容
sentence	str	解析対象文
**kwargs		特になし(β版より念のため)

感情分析: sentiment

# 実行例
api.sentiment("人生の春を謳歌しています")

def sentiment(self, sentence)

引数	データ型	内容
sentence	str	解析対象文

要約(β): summary

# 実行例
with open("summary.txt", "r", encoding="utf-8") as f:
    document = f.read()
api.summary(document, 1)

summary.txt

前線が太平洋上に停滞しています。一方、高気圧が千島近海にあって、北日本から東日本をゆるやかに覆っています。関東地方は、晴れ時々曇り、ところにより雨となっています。東京は、湿った空気や前線の影響により、晴れ後曇りで、夜は雨となるでしょう。

def summary(self, document, sent_len, **kwargs)

引数	データ型	内容
document	str	入力文章
sent_len	int	要約文数
**kwargs		特になし(β版より念のため)

3. 今後の予定

何か作ったやつを出したい

あるいは説明していない機能について説明したい

参考

pixivpy

内容

15stepで踏破 自然言語処理アプリケーション入門 を読み進めていくにあたっての自分用のメモです。
今回は4章Step14で、自分なりのポイントをメモります。

準備

• 個人用MacPC：MacOS Mojave バージョン10.14.6
• docker version：Client, Server共にバージョン19.03.2

章の概要

この章では機械学習システムに外部から与えるべき適切なパラメータ（学習するパラメータではない）の値を見つけることを目指す。

• グリッドサーチ
• Hyperopt
  • 確率分布

14.1 ハイパーパラメータ

学習によって調整、獲得されるパラメータではなく、学習の前に設計者やプログラマが設定する一段階メタなパラメータ。

• 特徴抽出器
• 識別器
• NN
  • 層の種類
  • 層の数
  • 層ごとのユニット数
  • dropoutの有無と係数
  • optimizerの種類と各種引数
  • 学習率
  • etc...

14.2 グリッドサーチ

探索対象のパラメータと、それぞれの値の候補を列挙し、その全ての組み合わせを試してベストのものを見つける手法である。
Scikit-learnでは、sklearn.model_selection.GridSearchCVが提供されている。

• Scikit-learnのAPIを持った識別器クラスを用いる
• Kerasの場合はscikit-leran APIラッパーを用いる
• 最適なパラメータは<gridsearchのインスタンス>.best_params_で取得できる。

# pipelineを使わない

## train
vectorizer = TfidfVectorizer(tokenizer=tokenize, ngram_range=(1, 2))
train_vectors = vectorizer.fit_transform(train_texts)
parameters = {  # <1>
    'n_estimators': [10, 20, 30, 40, 50, 100, 200, 300, 400, 500],
    'max_features': ('sqrt', 'log2', None),
}
classifier = RandomForestClassifier()
gridsearch = GridSearchCV(classifier, parameters)
gridsearch.fit(train_vectors, train_labels)

## predict
test_vectors = vectorizer.transform(test_texts)
predictions = gridsearch.predict(test_vectors)


# pipelineを使う

## train 
pipeline = Pipeline([
    ('vectorizer', TfidfVectorizer(tokenizer=tokenize, ngram_range=(1, 2))),
    ('classifier', RandomForestClassifier()),
])
parameters = {
    'vectorizer__ngram_range':[(1, 1), (1, 2), (1, 3), (2, 2), (2, 3), (3, 3)],
    'classifier__n_estimators':[10, 20, 30, 40, 50, 100, 200, 300, 400, 500],
    'classifier__max_features':('sqrt', 'log2', None),
}
gridsearch = GridSearchCV(pipeline, parameters)
gridsearch.fit(texts, labels)

## predict
gridsearch.predict(texts)

GridSearchCV

# verbose：グリッドサーチの実行状況がわからなかったので表示(1)
# n_jobs：可能な限り(-1)並列で実行
clf = GridSearchCV(pipeline, parameters, verbose=1, n_jobs=-1)

実行結果

from dialogue_agent import DialogueAgent  # <1>
↓
from dialogue_agent_pipeline_gridsearch import DialogueAgent  # <1>

$ docker run -it -v $(pwd):/usr/src/app/ 15step:latest python evaluate_dialogue_agent.py
# 実行完了まで20分ほどかかった
Fitting 3 folds for each of 180 candidates, totalling 540 fits

[Parallel(n_jobs=-1)]: Done  46 tasks      | elapsed:    5.4s
[Parallel(n_jobs=-1)]: Done 196 tasks      | elapsed:  2.2min
[Parallel(n_jobs=-1)]: Done 446 tasks      | elapsed:  5.5min
[Parallel(n_jobs=-1)]: Done 540 out of 540 | elapsed: 19.9min finished

0.7021276595744681
{'classifier__max_features': 'log2', 'classifier__n_estimators': 300, 'vectorizer__ngram_range': (1, 1)}

• 探索対象のパラメータが増えてくると、グリッドサーチの探索にかかる時間は指数関数的に長くなってしまう。
• GridSearchCV#fitの内部では、パラメータの組を1つ試すたびに性能評価が行われている
  • クロスバリデーション（デフォルト分割数は3）
  • 学習データをK個に分割しそのうち1個をバリデーションデータ、残り全てを学習データとして学習と評価を行い、それをバリデーションデータを変えながらK回行う評価方法

14.3 Hyperoptの利用 14.4 確率分布

グリッドサーチよりも効率よくハイパーパラメータの探索を行うためのツール。
パラメータ空間と目的関数を与えて最適なハイパーパラメータを返す。

• パラメータ空間：探索対象のパラメータとそれぞれの値の候補
  • 一様分布：どの値も出現する確率が全て等しい
  • 乱数一様分布：値の対数が一様分布に従う。大きい値はまばらに、小さい値は密に得ることができる。学習率などは対数の間隔で探索するのが望ましい。
    • 生成：Hyperoptの下限と上限で指定する際は値の対数を指定する。math.log(..)
• 目的関数：パラメータ値の組を受け取って値を返す関数
  • 精度を最大化したい場合は、マイナスをかけて最小化する（マイナスの最小化は最大化となる）
• 最適なパラメータは<hyperoptのインスタンス>.fminの返り値から取得できる。

# パラメータ探索
vectorizer = TfidfVectorizer(tokenizer=tokenize, ngram_range=(1, 2))
train_vectors = vectorizer.fit_transform(train_texts)

## 目的関数
def objective(args):
    classifier = RandomForestClassifier(n_estimators=int(args['n_estimators']),
                                        max_features=args['max_features'])
    classifier.fit(tr_vectors, tr_labels)
    val_predictions = classifier.predict(val_vectors)
    accuracy = accuracy_score(val_predictions, val_labels)
    return -accuracy

## パラメータ空間
max_features_choices = ('sqrt', 'log2', None)
space = {
    'n_estimators': hp.quniform('n_estimators', 10, 500, 10),
    'max_features': hp.choice('max_features', max_features_choices),
}
best = fmin(objective, space, algo=tpe.suggest, max_evals=30)

# train
best_classifier = RandomForestClassifier(
    n_estimators=int(best['n_estimators']),
    max_features=max_features_choices[best['max_features']])
best_classifier.fit(train_vectors, train_labels)

# predict
test_vectors = vectorizer.transform(test_texts)
predictions = best_classifier.predict(test_vectors)

14.5 Kerasへの応用

実行がなかなか終わらないので詳細は省略。

• セッションクリア
  • 探索中に何度もモデルをメモリ上に構築しているとGPUメモリを食いつぶしてしまうため、都度解放処理を入れる
  • if Keras.backend.backend() == 'tensorflow':
  • ....Keras.backend.clear_session()
• パラメータ
  • 選択肢によって探索項目が異なる場合は、パラメータ空間を入れ子にしてそれぞれの探索項目ごとに詳細項目を指定できる。
  • optimizerの例
    • SGD：learning rateとmomentumを探索
    • Adagrad；leraning rateのみ探索

評価

実行がうまく進んでいない（CPUだから時間がかかっているだけ、CPUだと実行できない？）ので、後ほど余裕があったら更新する。

1.はじめに

　簡単に言うと、ニューラルネットワークの豊かな表現力は、シンプルな活性化関数を入れ子にして深い階層にすることによってもたらされます。

　今回は、ニューラルネットワークで使われる活性化関数について勉強したことをまとめます。

2.Sigmoid関数

# シグモイド関数
def sigmoid(x):
    return 1 / (1 + np.exp(-x))   

# シグモイド関数の微分
def sigmoid_d(x):
    return (1 / (1 + np.exp(-x))) * ( 1- (1 / (1 + np.exp(-x))))

# グラフ表示
x = np.arange(-5.0, 5.0, 0.01)
plt.plot(x, sigmoid(x), label='sigmoid')
plt.plot(x, sigmoid_d(x), label='sigmoid_d')
plt.ylim(-1.1, 1.1)  
plt.legend()  
plt.grid()
plt.show()

sigmoid関数：

sigmoid(x) = \frac{1}{1+e^{-x}}

sigmoid関数の微分：

sigmoid'(x) = \frac{1}{1+e^{-x}} * ( 1 - \frac{1}{1+e^{-x}})

　Sigmoid関数は、昔からニューラルネットワークの教科書には必ず載っていて、微分してもほとんど形が変わらず美しい形をしているのですが、最近は活性化関数としてはほとんど使われていません。　

　その理由は、グラフで分かる様に、xの値が大きくなるとyが1に張り付いて動かなくなってしまうためです。ニューラルネットワークはyを微分して傾きを求めて重みパラメータを最適化するので、微分がほとんど0になってしまうと中々最適解に近づけない(勾配消失)という問題点を持っているからです。

Sigmoid関数の微分の導出

sigmoid'(x) = ((1 + e^{-x})^{-1})'\\

合成関数の微分により、u = 1+e^{-x}と置くと、\frac{dy}{dx}=\frac{dy}{du}\frac{du}{dx} なので\\

= -(1 + e^{-x})^{-2} * (1+e^{-x})'\\

= -\frac{1}{(1+e^{-x})^2} * -e^{-x}\\

= \frac{1}{1+e^{-x}} * \frac{e^{-x}}{1+e^{-x}} \\

= \frac{1}{1+e^{-x}} * (\frac{1+e^{-x}}{1+e^{-x}} - \frac{1}{1+e^{-x}})\\

= \frac{1}{1+e^{-x}} * ( 1 - \frac{1}{1+e^{-x}})

3.tanh関数

# Tanh関数
def tanh(x):
    return (np.exp(x) -np.exp(-x)) / (np.exp(x) + np.exp(-x))

# Tanh関数の微分
def tanh_d(x):
    return 1- ( (np.exp(x) -np.exp(-x)) / (np.exp(x) + np.exp(-x)) )**2

# グラフ表示
x = np.arange(-5.0, 5.0, 0.01)
plt.plot(x, tanh(x), label='tanh')
plt.plot(x, tanh_d(x), label='tanh_d')
plt.ylim(-1.1, 1.1)  
plt.legend()  
plt.grid()
plt.show()

Tanh関数：

tanh(x) = \frac{e^x-e^{-x}}{e^x+e^{-x}}

Tanh関数の微分：

tanh(x) = \frac{4}{(e^x + e^{-x})^2}\\

　Tanh関数は、Sigmoid関数の改良版(微分した時の最大値がSigmoidより高い)として使われていたわけですが、xが大きくなるとyが1に貼り付く根本的な問題点は改善されていません。

Tanh関数の微分の導出

商の微分公式　(\frac{f(x)}{g(x)})' = \frac{f'(x)*g(x) - f(x)*g'(x)}{g(x)^2}　を使って\\

tanh'(x) = \frac{(e^x+e^{-x})^2 - (e^x - e^{-x})^2}{(e^x + e^{-x})^2}\\

= \frac{e^{2x}+2+e^{-2x} - (e^{2x} -2 + e^{-2x})}{(e^x + e^{-x})^2}\\

= \frac{4}{(e^x + e^{-x})^2}\\

または、

tanh'(x) = \frac{(e^x+e^{-x})^2 - (e^x - e^{-x})^2}{(e^x + e^{-x})^2}\\

= 1 - \frac{(e^x - e^{-x})^2}{(e^x + e^{-x})^2}\\

= 1 - (\frac{e^x - e^{-x}}{e^x + e^{-x}})^2\\

4.ReLU関数

# ReLU関数
def relu(x):
    return np.maximum(0, x)

# ReLU関数の微分
def relu_d(x):
    return np.array(x > 0, dtype=np.int)

# グラフ表示
x = np.arange(-5.0, 5.0, 0.01)
plt.plot(x, relu(x), label='relu')
plt.plot(x, relu_d(x), label='relu_d')
plt.ylim(-1.1, 1.1)  
plt.legend()  
plt.grid()
plt.show()

ReLU関数：

ReLU関数の微分：

　Sigmoid関数の根本的な問題点をなくすために生まれたのがReLU関数です。xが大きくなってもyも比例して大きくなり常に微分すると定数項が残ります。何か、今聞くと当たり前の様な気がしますが、これが使われ出したのは、なんと2012年頃からです。

　東大の松尾豊先生は、Sigmoid関数はシンプルで微分してもほとんど形が変わらず、理工学者にとって美しい関数だった。一方ReLUはかっこが悪く、しかも(0,0)で微分不可な点があるので、誰も使いたがらなかった。

　昔ディープラーニングは上手く動かせなかったので、皆んな式として美しいSigmoid関数を使っていた。しかし、動かせる様になってからは色々なことを試す人が出て来て、そうした中で使われる様になって来たのがReLUだと言っています。

5.Leaky ReLU関数

# Leaky ReLU関数
def leaky_relu(x):
    return np.where(x > 0, x , 0.01 * x)

# Leaky ReLU関数の微分
def leaky_relu_d(x):
    return np.where(x>0,1,0.01)

# グラフ表示
x = np.arange(-5.0, 5.0, 0.01)
plt.plot(x, leaky_relu(x), label='leaky_relu')
plt.plot(x, leaky_relu_d(x), label='leaky_relu_d')
plt.ylim(-1.1, 1.1)  
plt.legend()  
plt.grid()
plt.show()

Leaky ReLU関数：

Leaky ReLU関数の微分：

　Leaky ReLU関数は、ReLU関数から派生したもので、xが0以下でも0.01xと傾きを持ちます。ReLU関数より、さらに最適化が進むことが期待されましたが、ReLU関数より最適化が上手く行く場合はかなり限定的な様です。

6.活性化関数の性能の違いを実感してみる

　最後に、活性化関数によって、どれだけ最適化性能が異なるのか実感してみましょう。

　TensorFlow PlayGround というブラウザでニューラルネットワークのシミュレーションができるサイトを覗いてみます。

　こんな３個のニューロンが2層のニューラルネットワークを設定し、赤枠のActivation（活性化関数)を切り替えて、収束時間をチェックしてみます。

　確率的な問題もあり、ある程度バラツキは出ますが、収束時間は概ね、TanhがSigmoidの10倍くらい早く、ReLUはTanhのさらに2倍くらい早いです。

$$
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
$$

はじめに

1980年代に量子計算の原理が明らかになった当初、古典的なデジタル処理で使われるような誤り訂正の理論がない、もしくはその理論確立が極めて困難なため、量子コンピュータなんていうものは現実的に実現不可能と考える人もいたらしいです。古典ビットの場合、ビットが反転するエラーだけを考えれば良いのに対し、量子ビットの場合、その変化は無限の自由度を持つので、無限のパターンのエラーを相手にする必要があります。ブロッホ球をイメージすると、量子状態のズレ方は無限にありますよね。そのズレた状態からもとの状態を正確に復元させることを考えないといけないです。さらに、量子状態は古典情報のように無邪気に複製できませんし、測定したら壊れてしまう代物です。といった諸々があり、難しさが半端ないと思われていたようです。それを最初に打破したのがPeter Shorです。このShorが提案した「Shorの符号(Shor code)」を使えば、1量子ビットに働く任意の誤りを訂正することができます。この「任意の」というところがスゴイところです。しかも、割と簡単な量子回路で実現できてしまいます。今回は、その理論を説明した後、量子計算シミュレータqlazyで、その威力を実感してみたいと思います。

参考にさせていただいたのは、以下の文献です。

1. ニールセン、チャン「量子コンピュータと量子通信（３）」オーム社（2005年）
2. 富田「量子情報工学」森北出版（2017年）
3. 中山「量子アルゴリズム」技報堂出版（2014年）

理論の確認

古典的誤り訂正（繰り返し符号）

量子誤り訂正の話に突入する前の準備体操として、古典情報における誤り訂正の基本的な考え方をおさらいしておきます。いま、1ビットの情報を雑音のある古典チャネルを通して送信することを考えてみます。このチャネルにより、ビットは確率pで反転するとします。つまり、入力ビットが0だった場合、確率pでビットが1に反転し、確率1-pで0のまま維持されます。また、入力ビットが1だった場合、確率pでビットが0にに反転し、確率1-pで1のまま維持されます。このようなチャネルのことを「2値対称チャネル(binary symmetric channel)」と言います。何も考えず、このチャネルをそのまま使うと、当然ですが誤り率はpです。このチャンネルを前提に、何らかの処理を施すことによって、誤り率を減らすことができないでしょうか。

基本的で、よく知られたやり方は、ビットを複製してわざと冗長化しておく方法です。例えば、与えられた1ビットの情報を、以下のように3ビットに変換（符号化）しておきます。

\begin{align}
0 &\rightarrow 000 \\
1 &\rightarrow 111 \\
\end{align}

これを次々に古典チャネルを通して送信します。情報量は3倍になってしまいますが、誤り訂正のためには仕方ありません、という考え方です。雑音によって確率的にビットが反転するので、000というビット列は、受信側で010になっているかもしれません。が、これは真ん中のビットが反転したせいであり、本当は000のことだろうと推測するわけです。つまり多数決をとって、多いビットが元のビットだろうと見なす手法です。入力ビットを繰り返して符号化するので「繰り返し符号(repetition code)」と呼ばれています。容易にわかる通り、入力が000で雑音の結果110になってしまう可能性もなくはないです。が、トータルでの平均誤り率が小さければ良いと考えます。では、その平均誤り率は、具体的にどう計算されるでしょうか。0個もしくは1個のビット反転は訂正できるので、2個または3個のビットが反転してしまう確率だけを計算すれば、所望の誤り率になります。3個いっぺんに反転してしまう確率は$p^3$で、2個反転してしまう確率は3パターンあって$3p^2(1-p)$です。これを足すと、$p^3 + 3p^2(1-p) = 3p^2 - 2p^3$になります。何もやらない場合の誤り率は$p$だったので、誤り訂正をした方が誤り率が小さくなるという条件は、$p > 3p^2 - 2p^3$、すなわち、$p < 1/2$となります。つまり、$p$の値が$1/2$よりも小さい古典チャネルであれば、この繰り返し符号の効果が発揮されるというわけです。

３量子ビットのビット反転符号

符号化

量子ビットの場合も、似たような考え方で符号化します。例えば、1量子ビットの状態が$\ket{\psi} = a \ket{0} + b \ket{1}$だったとして、それを$a \ket{000} + b \ket{111}$のように変換（符号化）できないかを考えます。簡単にわかると思いますが、以下のような量子回路を通せば実現できます。

|psi> --*--*-- ...
        |  |
  |0> --X----- ...
           |  
  |0> -----X-- ...

この状態で、例えば一つのビットが反転して$\ket{000} \rightarrow \ket{010}$のように変化しても、元の状態を復元することができます。ということを、これから説明していきます（図中の点々部分は雑音チャネルであることを表すものとします。以下同様）。

誤り検出（シンドローム診断）

まず、どのビットが反転したかを知る必要があります。この手続きのことを「誤り検出(error-detection)」とか「シンドローム診断(syndrome diagnosis)」と言います。3ビットの場合、以下のような射影演算を施せば、それがわかります。

\begin{align}
P_0 &= \ket{000} \bra{000} + \ket{111} \bra{111} \\
P_1 &= \ket{100} \bra{100} + \ket{011} \bra{011} \\
P_2 &= \ket{010} \bra{010} + \ket{101} \bra{101} \\
P_3 &= \ket{001} \bra{001} + \ket{110} \bra{110} \tag{1}
\end{align}

$P_0$は何も誤りがないことを検出する測定に対応した射影演算子です。$P_1$は1番目のビットに誤りがあることを検出する射影測定です。以下、$P_2$は2番目、$P_3$は3番目のビットに誤りがあることを検出する射影測定です。雑音の結果、$a \ket{100} + b \ket{011}$となった場合、$\bra{\psi} P_1 \ket{\psi}$のみが1で、その他はすべて0になりますので、確かに$P_1$は1番目のビットの誤り検出のための演算子であるということがわかると思います。

この例以外にも、シンドローム診断のための射影測定があります。式(1)では4つの演算子を使いましたが、2つで済ますことができるスグレモノがあります。結論を先に言います。1番目、2番目、3番目の量子ビットに適用するパウリ$Z$演算子を各々$Z_1,Z_2,Z_3$としたときに、$Z_1 Z_2$および$Z_2 Z_3$というオブザーバブルを測定するやり方です。とイキナリ言われても「ん？」という状態かもしれないので、少し説明を加えます。

パウリ$Z$は$Z=\ket{0}\bra{0}-\ket{1}\bra{1}$と書けるので、

\begin{align}
Z_1 Z_2 &= (\ket{0}_{1}\bra{0}_{1}-\ket{1}_{1}\bra{1}_{1})(\ket{0}_{2}\bra{0}_{2}-\ket{1}_{2}\bra{1}_{2}) \\
&= \ket{00}\bra{00} + \ket{11}\bra{11} - \ket{01}\bra{01} - \ket{10}\bra{10} \\
&= (\ket{00}\bra{00} + \ket{11}\bra{11}) \otimes I_3 - (\ket{01}\bra{01} + \ket{10}\bra{10}) \otimes I_3 \tag{2}
\end{align}

と展開できます。この式をじっと睨むとわかると思いますが、1番目と2番目の量子ビットが同じ場合（「偶パリティ」と言います）、このオブザーバブルを測定すると1を出力します。逆に、1番目と2番目の量子ビットが違う場合（「奇パリティ」と言います）、−1を出力します¹。つまり、$Z_1 Z_2$の測定を行うと、1番目と2番目の量子ビットのパリティをチェックすることができます。同様に、$Z_2 Z_3$の測定によって、2番目と3番目の量子ビットのパリティがわかります。

表に整理すると、以下のようになります。

$Z_1 Z_2$の測定値（パリティ）	$Z_2 Z_3$の測定値（パリティ）	反転しているビット
+1（偶）	+1（偶）	なし
+1（偶）	-1（奇）	3番目
-1（奇）	+1（偶）	1番目
-1（奇）	-1（奇）	2番目

1行目は、1番目と2番目が偶パリティで、2番目と3番目が偶パリティということなので、反転しているビットがないと推定できそうです。

2行目は、1番目と2番目が偶パリティで、2番目と3番目が奇パリティということなので、3番目が反転していると推定できそうです。

3行目は、1番目と2番目が奇パリティで、2番目と3番目が偶パリティということなので、1番目が反転していそうです。

最後の4行目は、1番目と2番目が奇パリティで、2番目と3番目が奇パリティということなので、真ん中の2番目が反転している可能性が高そうです。

ということで、反転している量子ビットが何番目なのか、推定することができます。反転している量子ビットが、高々1個しかない場合は、反転ビットが何番目かを完全に特定することができます。

このシンドローム診断は、どのような量子回路で実現できるでしょうか。測定すれば良いだけなので、

|psi> --*--*-- ... --M--
        |  |
  |0> --X----- ... --M--
           |  
  |0> -----X-- ... --M--

としてしまうかもしれません。が、これでは駄目です。誤り訂正したい量子ビットを測定によって壊してしまっています。なので、誤り訂正したい量子状態を壊さないように、以下のように、補助量子ビット（アンシラ）を追加して、そちらに対象とする量子状態を逃してあげて、その状態を測定するようにしないといけません。

|psi> --*--*-- ... --*------------
        |  |         |
  |0> --X----- ... -----*--*------
           |         |  |  |
  |0> -----X-- ... -----------*---
                     |  |  |  |
               |0> --X--X--------M1
                           |  |
               |0> --------X--X--M2

ここで、$Z_1 Z_2$の測定をM1、$Z_2 Z_3$の測定をM2と表しました²。

回復

シンドローム診断で、反転している量子ビットが特定できたので、次に、その量子ビットを再度反転して元に戻せば良いです。M1とM2の値に応じて、XゲートをON/OFFします。

|psi> --*--*-- ... --*------------X---
        |  |         |            ||
  |0> --X----- ... -----*--*------X---
           |         |  |  |      ||
  |0> -----X-- ... -----------*---X---
                     |  |  |  |   ||
               |0> --X--X---------M1
                           |  |   ||
               |0> --------X--X---M2

測定結果に応じてゲートを制御することを二重の縦棒で表してみました（わかりにくく、かつ、正式な書き方でもないですが、気持ちを汲んでください、汗）。先程表で示したパターンで３つの量子ビットのどれかにXゲートをかけます。これで、元の状態に戻ることになります³。

測定結果に応じたゲート操作をする代わりに、CNOTとToffoliを使って同じことができます。

              [雑音]
|psi> --*--*-- ... --*-------------------X--X--
        |  |         |                   |  |
  |0> --X----- ... -----*--*----------------X--
           |         |  |  |             |  |
  |0> -----X-- ... -----------*-------X-----X--
                     |  |  |  |       |  |  |
               |0> --X--X----------------*--*--
                           |  |       |     |
               |0> --------X--X-------*-----*--

     [符号化]     [シンドローム診断]  [訂正]

これはわかりますでしょうか。

シンドローム診断の段階で、1番目のアンシラが+1($\ket{0}$)で2番目のアンシラが+1($\ket{0}$)の場合、[訂正]と記載されている回路部の制御ビットはともに$\ket{0}$なので、何もしません。

1番目のアンシラが+1($\ket{0}$)で2番目のアンシラが-1($\ket{1}$)の場合、訂正回路の2番目のアンシラだけが$\ket{1}$になるので、3番目の量子ビットにXゲートがかかり反転します。

1番目のアンシラが-1($\ket{1}$)で2番目のアンシラが+1($\ket{0}$)の場合、訂正回路の1番目のアンシラだけが$\ket{1}$になるので、1番目の量子ビットにXゲートがかかり反転します。

1番目のアンシラが-1($\ket{1}$)で2番目のアンシラが-1($\ket{1}$)の場合、訂正回路の両方のアンシラが$\ket{1}$になるので、2番目の量子ビットにXゲートがかかり反転します。

ということで、これで、すべての誤りパターンに対してビットを回復することができました。

復号化

最後に、符号化の逆をやって、もともとあった状態に完全に戻します。

              [雑音]
|psi> --*--*-- ... --*-------------------X--X---*--*-- |psi>
        |  |         |                   |  |   |  |
  |0> --X----- ... -----*--*----------------X------X-- |0>
           |         |  |  |             |  |   |
  |0> -----X-- ... -----------*-------X-----X---X----- |0>
                     |  |  |  |       |  |  |
               |0> --X--X----------------*--*--
                           |  |       |     |
               |0> --------X--X-------*-----*--

     [符号化]     [シンドローム診断]  [訂正]   [復号化]

以上で、1量子ビットの反転に対する誤り訂正回路が完成しました。

しかし、アンシラとして2量子ビット余分に必要になってしまうのが、よろしくないです。できればアンシラなしで済ませたいということで、考え出された回路を以下に示します。1番目の量子ビットだけ回復できれば良いので、2番目、3番目はもとに戻らなくても良いという考え方です。ただし、少なくとも1番目と他の量子ビットとの間にもつれが残らないようにします。

|psi> --*--*-- ... --*--*--X-- |psi>
        |  |         |  |  |
  |0> --X----- ... -----X--*--
           |         |     |
  |0> -----X-- ... --X-----*--

どうでしょうか。かなりシンプルになっていますが、これで本当に誤り訂正できているのでしょうか。数式で確認してみます。入力する状態は$\ket{\psi} = a \ket{0} + b \ket{1}$とします。符号化した段階で状態は$a \ket{000} + b \ket{111}$となるのは明らかなので、その直後からの変化を見てみます。

反転する量子ビットがない場合、

\begin{align}
a\ket{000} + b\ket{111} &\rightarrow a\ket{000} + b\ket{110} \rightarrow a\ket{000} + b\ket{100} \rightarrow a\ket{000} + b\ket{100} \\
&= (a\ket{0} + b\ket{1}) \otimes \ket{00}  \tag{3}
\end{align}

雑音により、1番目の量子ビットが反転した場合、

\begin{align}
a\ket{100} + b\ket{011} &\rightarrow a\ket{101} + b\ket{011} \rightarrow a\ket{111} + b\ket{011} \rightarrow a\ket{011} + b\ket{111} \\
&= (a\ket{0} + b\ket{1}) \otimes \ket{11}  \tag{4}
\end{align}

雑音により、2番目の量子ビットが反転した場合、

\begin{align}
a\ket{010} + b\ket{101} &\rightarrow a\ket{010} + b\ket{100} \rightarrow a\ket{010} + b\ket{110} \rightarrow a\ket{010} + b\ket{110} \\
&= (a\ket{0} + b\ket{1}) \otimes \ket{10}  \tag{5}
\end{align}

雑音により、3番目の量子ビットが反転した場合、

\begin{align}
a\ket{001} + b\ket{110} &\rightarrow a\ket{001} + b\ket{111} \rightarrow a\ket{001} + b\ket{101} \rightarrow a\ket{001} + b\ket{101} \\
&= (a\ket{0} + b\ket{1}) \otimes \ket{01}  \tag{6}
\end{align}

ということで、2番目と3番目の量子ビットは必ずしも$\ket{00}$に戻らないですが、少なくとも1番目の量子ビットはすべてのパターンで$\ket{\psi} = a\ket{0} + b\ket{1}$に戻っていて、かつ、他の量子ビットとはもつれていないです（＝積状態になっています）。

3量子ビットの位相反転符号

次に、1つの量子ビットに位相反転があった場合の誤り訂正についてです。ビット反転はXゲートで表現できますが、位相反転はZゲートの作用です。つまり、$a\ket{0}+b\ket{1} \rightarrow a\ket{0}-b\ket{1}$という変化です。この状態から回復したいわけですが、どうすれば良いでしょうか。実は簡単な処方箋があります。アダマールゲートをかければ良いです。アダマールゲートの作用は、

\begin{align}
\ket{0} &\rightarrow \ket{+} = \frac{1}{2} (\ket{0} + \ket{1})  \\
\ket{1} &\rightarrow \ket{-} = \frac{1}{2} (\ket{0} - \ket{1})  \tag{7}
\end{align}

だったので、$a\ket{0}+b\ket{1} \rightarrow a\ket{0}-b\ket{1}$という位相反転は、

\begin{align}
& \frac{1}{2} (a+b) \ket{+} \frac{1}{2} (a-b) \ket{-} \\
& \rightarrow \frac{1}{2} (a+b) \ket{-} \frac{1}{2} (a-b) \ket{+}  \tag{8}
\end{align}

と表されます。つまり、位相反転は$\{ \ket{+}, \ket{-} \}$基底におけるビット反転に他なりません。したがって、先程のビット反転の回路において符号化の後にアダマールをかけてから、雑音チャネルに入れて、復号化の直前にアダマールをかけて元に戻すということをやれば良いです。

|psi> --*--*--H-- ... --H--*--*--X-- |psi>
        |  |               |  |  |
  |0> --X-----H-- ... --H-----X--*--
           |               |     |
  |0> -----X--H-- ... --H--X-----*--

これで、位相反転に対する誤り訂正ができたことになります。

Shorの符号

「Shorの符号」は、ビット反転と位相反転の誤り訂正を同時に実現するものです。位相反転の誤り訂正では、初期状態の1量子ビットを3量子ビットに冗長化させ、各々にアダマールをかけます。その各量子ビットに対してビット反転の誤り訂正も実行するため、各々3量子ビットに冗長化します。その結果、合計で9量子ビット必要になります。回路図で表すと、

|psi> --*--*--H----*--*-- ... --*--*--X--H--*--*--X---- |psi>
        |  |       |  |         |  |  |     |  |  |
        |  | |0> --X----- ... -----X--*--   |  |  |
        |  |       |  |         |     |     |  |  |
        |  | |0> -----X-- ... --X-----*--   |  |  |
        |  |                                |  |  |
  |0> --X-----H----*--*-- ... --*--*--X--H-----X--*--
           |       |  |         |  |  |     |     |
           | |0> --X----- ... -----X--*--   |     |
           |          |         |     |     |     |
           | |0> -----X-- ... --X-----*--   |     |
           |                                |     |
  |0> -----X--H----*--*-- ... --*--*--X--H--X-----*--
                   |  |         |  |  |
             |0> --X----- ... -----X--*--
                      |         |     |
             |0> -----X-- ... --X-----*--

となります。これで、ビット反転と位相反転の誤りに対して元の状態を回復することができるようになります（もちろん、いま想定しているのは、どれか１つの量子ビットのみが雑音によって反転する場合です）。が、実は、このように構成されたShorの符号は、さらに強力な効果を持っています。ビット反転や位相反転のようなわかりやすい離散的な誤りだけでなく、任意の連続的な変化（誤り）の場合であっても、きちんと訂正することができます。

例えば、符号化後の状態が$\ket{\psi}$になったとします。この1番目の量子ビットだけが、Kraus演算子$\{ E_{i} \}$で記述される雑音を受けたとすると、その変化は、

\ket{\psi}\bra{\psi} \rightarrow \sum_{i} E_{i} \ket{\psi}\bra{\psi} E_{i}^{\dagger}  \tag{9}

のように書けます。ここで、$E_i$は1番目の量子ビットに働くパウリ演算子$X_1, Z_1$と複素係数$\{ e_{ij} \}$を使って、

E_i = e_{i0} I + e_{i1} X_1 + e_{i2} Z_1 + e_{i3} X_1 Z_1  \tag{10}

と展開できます⁴。

式(9)の一つの$i$に関する項だけに注目すると、量子状態は、$\ket{\psi},X_{1}\ket{\psi},Z_{1}\ket{\psi},X_{1}Z_{1}\ket{\psi}$の重ね合わせ状態に変化するということです。つまり、状態$\ket{\psi}$そのもの、状態$\ket{\psi}$のビット反転、状態$\ket{\psi}$の位相反転、状態$\ket{\psi}$のビット・位相反転という4つの状態の重ね合わせになります。これに、シンドローム診断を施すことで4つのうちのどれだったかが判明します。判明したらば次のステップで誤りのパターンに応じたゲートを通すことで、量子状態を回復します。式(9)はすべての$i$についての和になっていますが、同じ議論で他の$i$についても誤り訂正が可能になりますし、1番目以外の量子ビット番号についても同様の議論により、誤り訂正が可能です。

シミュレーション

実装

それでは、Shorの符号による誤り訂正の回路図をシミュレーションしてみます。全体のPythonコードを示します。

import random
from qlazypy import DensOp

def create_densop():

    de_ini = DensOp(qubit_num=9).h(0)
    de_fin = de_ini.clone()
    return de_ini, de_fin

def noise(self, kind='', prob=0.0, qid=[]):

    qchannel = {'bit_flip':self.bit_flip,
                'phase_flip':self.phase_flip,
                'bit_phase_flip':self.bit_phase_flip,
                'depolarize':self.depolarize,
                'amp_dump':self.amp_dump,
                'phase_dump':self.phase_dump}
    [qchannel[kind](i, prob=prob) for i in qid]
    return self

def code(self):

    self.cx(0,3).cx(0,6)
    self.h(0).h(3).h(6)
    self.cx(0,1).cx(0,2)
    self.cx(3,4).cx(3,5)
    self.cx(6,7).cx(6,8)
    return self

def correct(self):

    self.cx(0,2).cx(0,1)
    self.cx(3,5).cx(3,4)
    self.cx(6,8).cx(6,7)
    self.ccx(2,1,0).ccx(5,4,3).ccx(8,7,6)
    self.h(0).h(3).h(6)
    self.cx(0,3).cx(0,6)
    self.ccx(6,3,0)
    return self

if __name__ == '__main__':

    # add custom gate
    DensOp.add_method(code)
    DensOp.add_method(noise)
    DensOp.add_method(correct)

    # settings
    kind = 'depolarize' # bit_flip,phase_flip,bit_phase_flip,depolarize,amp_dump,phase_dump
    prob = 1.0
    qid = [0]
    print("== settings ==")
    print("* kind of noise        =", kind)
    print("* probability of noise =", prob)
    print("* noisy channels       =", qid)

    # error correction (shor code)
    de_ini, de_fin = create_densop()
    de_fin.code()
    de_fin.noise(kind=kind, prob=prob, qid=qid)
    de_fin.correct()

    # evaluate fidelity
    fid = de_fin.fidelity(de_ini, qid=[0])
    print("== result ==")
    print("* fidelity = {:.6f}".format(fid))

    # free all densops
    DensOp.free_all(de_ini, de_fin)

何をやっているか、順に説明します。まず、main処理部を見てください。

# add custom gate
DensOp.add_method(code)
DensOp.add_method(noise)
DensOp.add_method(correct)

で、必要となる処理をカスタムゲートとして登録します。code,noise,correctという関数を上の方で定義していて、これをQStateクラスのメソッドとして使えるようにします。

# settings
kind = 'depolarize' # bit_flip,phase_flip,bit_phase_flip,depolarize,amp_dump,phase_dump
prob = 1.0
qid = [0]

で、シミュレーションのためのパラメータを設定しています。kindは適用したい雑音のパターンを表しています。'bit_flip'（ビット反転）、'phase_flip'（位相反転）、'bit_phase_flip'（ビット位相反転）、'depolarize'（分極解消）、'amp_dump'（振幅ダンピング）、'phase_dump'（位相ダンピング）の中から適用したいパターン名を選びます。以前の記事で説明した代表的な量子チャネルを密度演算子に適用できます（密度演算子への量子チャネル適用はv0.0.34で追加しました）。probは、各量子チャネルで雑音を乗せる確率（強度）を表しています。qidは、どの量子ビットに雑音適用するかを設定するリストです。上の例では0番目の量子ビットに適用しています。複数のビット番号を指定することもできます。

de_ini, de_fin = create_densop()

で、量子状態（密度演算子）を適当に作ります。関数定義を見ていただければわかるとおり、合計9個の量子ビットを$\ket{0}$に初期化して、0番目のみにアダマールゲートをかけます。別にアダマールでなくても良いのですが、ここは適当です。どんな入力状態を想定するかに応じて、適当に変えれば良いです。また、ここでは、de_ini,de_finという2つの全く同じ量子状態（密度演算子）を出力していますが、この後、de_finだけに処理を加えて、最後に元の状態de_iniと比較したいために2つ出力するようにしています。

de_fin.code()

で、符号化を行います。関数定義を見ていただければ、一目瞭然かと思います。雑音が入る直前までの量子回路を愚直に実装しています。

de_fin.noise(kind=kind, qid=qid, prob=prob)

で、kindで定義された量子チャネルを確率probの強度でqidで規定される量子番号に適用します。関数定義の中身はパッと見、わかりにくいかもしれませんが、実質、de_ini.depolarize(0, prob=1.0)という演算をやっているだけです。量子回路図で言うと、点々の部分に相当します。

dde_fin.correct()

で、誤り訂正と復号化を一緒にやります。関数定義を見ていただければわかるとおりです（こちらも愚直に実装しているだけです）。

fid = de_fin.fidelity(de_ini, qid=[0])
print("== result ==")
print("* fidelity = {:.6f}".format(fid))

で、最初の状態と誤り訂正後の状態を比較するため、0番目の量子ビットに関する忠実度を計算して表示します。

DensOp.free_all(de_ini, de_fin)

最後、使ったメモリを解放します。クラス・メソッドfree_allで引数に指定したものを一気に解放することができます(v0.0.33で追加しました)。

これで、0番目の量子ビットに対して、確率1で分極解消がなされた（つまり完全に方向をバラバラにするような雑音が加えられた）場合の誤り訂正シミュレーションできるようになります。

実行結果

実行結果を示します。

== settings ==
* kind of noise        = depolarize
* probability of noise = 1.0
* noisy channels       = [0]
== result ==
* fidelity = 1.000000

というわけで、1量子ビットの雑音に関しては、完全に誤り訂正できました。

また、分極解消以外の量子チャネルの場合でも誤り訂正はできましたし、どの量子ビットに雑音が入っても大丈夫でした。以下はその一例です。

== settings ==
* kind of noise        = amp_dump
* probability of noise = 1.0
* noisy channels       = [5]
== result ==
* fidelity = 1.000000

ところが、雑音チャネルが2つ以上の場合は駄目でした。以下の通りです。

== settings ==
* kind of noise        = depolarize
* probability of noise = 1.0
* noisy channels       = [0, 2]
== result ==
* fidelity = 0.866025

確率を0.5にしてみるとどうでしょうか？

== settings ==
* kind of noise        = depolarize
* probability of noise = 0.5
* noisy channels       = [0, 2]
== result ==
* fidelity = 0.968246

こちらも駄目でした。上で説明した通り、Shorの符号で完全に誤り訂正はできるのは1量子ビットの誤りのみです。ということが、わかりました。

おわりに

今回、量子誤り訂正として各種提案されている手法の中で一番簡単でとっつきやすい手法である「Shorの符号」を取り上げて説明しました。ビット反転と位相反転に対応した誤り訂正回路を作ったはずなのに、なぜかすべての連続的な誤りにも対応できる回路になっていました、というのは面白いです。量子状態に対する連続的な演算が、実は離散的なパウリ演算子の重ね合わせで表現できるということが効いているのだと思います。

というわけで、量子誤り訂正の基本の入門編でした。今後しばらくは量子誤り訂正のあれこれを味わいながら、基本の理解を進めていく予定です。

以上

---

1. ここで「オブザーバブルを測定する」という表現が、ちょっとわかりにくいかもしれないので、少し説明を加えます。オブザーバブルというのは直訳すると「観測可能量」です。量子力学の枠組みでは、物理的に観測できる量は、エルミート演算子として表現できるとされています。その代表選手は、エネルギーという物理量をエルミート演算子として表現したハミルトニアンです。物理量を測定するということは、それに対応したオブザーバブル＝エルミート演算子の固有値問題を解くことと等価で、測定値が固有値（エルミートなので必ず実数値になります）、測定後の状態がその固有値に対応した固有ベクトル（固有状態）となります。いまオブザーバブルとして与えられているのはパウリ$Z$（の積）です。パウリ$Z$を測定するというのは、物理的にはどういうことかというと、ハミルトニアン$H = Z$で与えられた系のエネルギー値を測定するということです。つまり、1スピン系のZ軸方向の測定ですね。もっと言うと、量子回路でおなじみの量子ビットの測定そのものと思って良いです。結局、$Z_1 Z_2$という「オブザーバブルを測定する」というのは、何のことはない、1番目の量子ビットを計算基底で測定して-1か1の値を得て2番目の量子ビットを計算基底で測定して-1か1の値を得て、それらの積をとるということです。このとき1番目と2番目の測定の順番はどっちでも良いです（演算子として交換可能なので）。ちなみに、パウリ$X$というオブザーバブルを測定する、というのは、量子ビットで言うと、X軸方向の測定を実行するということになりますので、アダマールゲートをかけてからZ軸方向の測定をやることに等しいです。 ↩

2. 1番目と2番目の量子ビットを制御ビットとするCNOTを2つ通した後で標的ビットを測定しています（M1の測定）が、これでなぜ$Z_1 Z_2$を測定したことになるか？$\ket{x}\ket{y}\ket{0} \rightarrow \ket{x}\ket{y}\ket{x \oplus y}$と書いてみると、$x$と$y$が同じ場合$\ket{x}\ket{y}\ket{0}$となり、M1の測定値は+1、違う場合$\ket{x}\ket{y}\ket{1}$となり、M1の測定値は-1になります。 ↩

3. もちろん今の前提は雑音で反転する量子ビット数が1以下の場合なので、2ビット以上反転している場合は復元できません。 ↩

4. ここで、$X_1 Z_1$というのは、係数分を除いて$Y_1$に等しいので、$E_i$がパウリ演算子を使って展開できるということです。任意の2次元行列に対して、このような展開ができます。 ↩

前回からのおさらい

前回、Hello World編にて、標準出力に文字列を表示するという処理を確認しました。
今回はその対になります標準入力にてキーボードにて入力する文字列を受け取る処理を比較してみたいと思います。

標準入力

簡単にキーボードから1行入力された文字列を表示するプログラムを見ていくことにしましょう。

Goの標準入力

Goでの標準入力のプログラムです。

input.go

package main

import (
  "fmt"
  "os"
  "bufio"
)

func main() {
  stdin := bufio.NewScanner(os.Stdin)
  stdin.Scan()
  text := stdin.Text()

  fmt.Printf(text + "\n")
}

実行してみましょう。

$ go run input.go 
Hello, Go
Hello, Go

ただ単に鸚鵡返しするだけですが、思い通りの処理はできています。

プログラムの中身ですが、
importで３つのパッケージを読み込んでいます。fmtパッケージは標準出力でも使用したものですね。osパッケージとbufioパッケージは今回初めて使用するライブラリです。
osパッケージは公式サイトによると、

オペレーティングシステムの機能へのプラットフォームに依存しないインタフェースを提供します。
Unixライクな設計です。

ドキュメントを読んでるとsyscallパッケージをラップしていることがわかりました。osパッケージはシステムのローレベルな処理から汎用的にレベルを上げたパッケージということですね。

bufioパッケージは、

I/Oのバッファリング機能を提供します。io.Reader・io.Writerをラップし、別のオブジェクト(ReaderまたはWriter)を作成します。インタフェースは同じままでバッファリングやその他便利な機能を追加します。

とのことでioパッケージをラップして機能拡張しているものの様です。

処理に関してですが、１０行目のNewScannerはドキュメントによると１行ごと読み込む処理です。
またos.Stdinの処理はsyscall.Stdinで/dev/stdinを開いたファイルディスクリプたということで、標準入力を１行ごと読み込む処理の様です。
11行目で１行ごとバッファし、１２行目で変数に文字列を格納、そして１５行目で出力という流れになっています。

簡単な標準入力の処理ですが、なかなか理解しなければならないことの多いGoです。

では、他の言語との比較です。

Python

Pythonでは組み込み関数のinput関数で入力処理を実装できます。
組み込み関数

input.py

text = input()
print(text)

実行結果

$ python input.py
Hello, Python
Hello, Python

Ruby

Rubyでは組み込み関数のgetsで実装可能です。
module Kernel

input.rb

text = gets
print text

実行結果

$ ruby input.rb
Hello, Ruby
Hello, Ruby

Perl

PerlではI/O演算子の<>で標準入力のファイルハンドルSTDINを指定して入力処理を実装します。
I/O 演算子

input.pl

$text = <STDIN>;
print $text;

実行結果

$ perl input.pl
Hello, Perl
Hello, Perl

Bash

Bashではreadコマンドで入力処理を実装します。
man read

input.sh

read text
echo $text

実行結果

$ bash input.sh
Hello, Bash
Hello, Bash

入力処理を比較してみて

入力処理を比較してみて感じたことはHello World編と同様で、本体は簡素に作られているというところです。
今回は多少パッケージのソースにも目を通してみたところ、低レベルの処理を駆使して実装はできるのでしょうが、
車輪の再発明は無駄で何も旨味がないためパッケージの使用を学ぶべきと感じました。
(但し、実装方法を確認することは多くの知見を得られますので、できる限りソースには目を通すべきだと思います)
では、次回以降は変数に関しての演算子や文字列処理について勉強したいと思います。

背景

Data Augmentationした後の画像を表示したい！

と思って実装してみました。

Data Augmentationとは、1枚の画像を水増しする技術であり、以下のような操作を加えます。

• Random Crop(画像をランダムに切り取る)
• Random Horizontal Flip(画像を一定の確率で左右反転する)
• Random Erasing(画像の一部にランダムにノイズを付加する)
• Random Affine(画像をランダムに拡大・縮小・回転する)

この他にもいろいろあります。

実装

今回は、CIFAR-10の訓練画像データセットを読み込んで、transformsにRandomHorizontalFlipとRandomErasingを組み込んでみました。

test.py

import torch
import numpy as np
import torchvision
import torchvision.transforms as transforms
from torch.utils.data import Dataset,DataLoader
import torchvision.datasets as dsets
import matplotlib.pyplot as plt

#画像の読み込み
batch_size = 100
train_data = dsets.CIFAR10(root='./tmp/cifar-10', train=True, download=False, transform=transforms.Compose([transforms.RandomHorizontalFlip(p=0.5), transforms.ToTensor(), transforms.RandomErasing(p=0.5, scale=(0.02, 0.4), ratio=(0.33, 3.0))]))
train_loader = DataLoader(train_data,batch_size=batch_size,shuffle=True)
test_data = dsets.CIFAR10(root='./tmp/cifar-10', train=False, download=False, transform=transforms.Compose([transforms.ToTensor(),]))
test_loader = DataLoader(test_data,batch_size=batch_size,shuffle=False)

def image_show(data_loader,n):

  #Augmentationした画像データを読み込む
  tmp = iter(data_loader)
  images,labels = tmp.next()

  #画像をtensorからnumpyに変換
  images = images.numpy()

  #n枚の画像を1枚ずつ取り出し、表示する
  for i in range(n):
    image = np.transpose(images[i],[1,2,0])
    plt.imshow(image)
    plt.show()

image_show(train_loader,10)

image_show関数がAugmentation後の画像を表示する関数です。

iter()により、DataLoaderからミニバッチ1つ分を取得します。

そして、.next()により画像データをimagesに、ラベルをlabelsに格納します。

images = images.numpy()では、画像データをテンソルからnumpyに変換しています。

この時点でimagesは[バッチサイズ, チャンネル数, 幅, 高さ]という構造になっていますが、matplotlibのpyplotで画像を表示するには[幅, 高さ, チャンネル数]とする必要があります。

よって、np.transposeをつかって変形しています。

実行結果例

左右反転されていたりRandom Erasingでノイズが付加されていたりすることが確認できました。

はじめに

DjangoでCustom User Modelを設定します。

本記事はDjangoアプリをDocker上に構築しAWS Fargateにデプロイするプロジェクトの一部です。

Custom User Modelを使う理由

　Djangoにはデフォルトでユーザのモデルが設定されています。しかし、実際にアプリを作成すると、自分の作成したアプリに合わせてユーザモデルを変更したくなる時が必ずやってきます。ある程度コーディングが進んだ段階でユーザモデルの修正を行うのはとても難しい為、プロジェクトの最初からCustom User Modelを設定することが推奨されています。

Custom User Modelの実装について

Custom User Modelを実装するには「AbstractUser」クラスまたは「AbstractBaseUser」クラスを利用します。「AbstractUser」クラスはデフォルトのユーザフィールドや認証を拡張する方法でありシンプルに実装が可能です。「AbstractBaseUser」クラスを使う方法はコーディング量が増えますが、より細かくカスタマイズができます。

今回は「AbstractUser」を使用し実装していきます。自身が開発に携わったプロジェクトでも「AbstractUser」を使用してユーザモデルをカスタマイズしました。一般的なユーザ管理をするのであれば「AbstractUser」でも十分対応できます。

準備

Docker+DjangoでPostgreSQLを使う で作成したプロジェクトにCustom User Modelを実装します。

usersアプリの追加

まずはstartappコマンドでカスタムユーザを管理するusersアプリを作成します。その後setting.pyにCustom User Model関連のコードを追加します。

docker-compose exec web python manage.py startapp users

setting.pyに作成したusersを追加しましょう。また末尾にデフォルトのユーザモデルの代わりにカスタムユーザモデルを使用する旨を記載します。

setting.py

INSTALLED_APPS = [
    'django.contrib.admin',
    'django.contrib.auth',
    'django.contrib.contenttypes',
    'django.contrib.sessions',
    'django.contrib.messages',
    'django.contrib.staticfiles',

    # Local
    'users.apps.UsersConfig', #追加
]
...
AUTH_USER_MODEL = 'users.CustomUser' # 追加

modelへCustomUserを追加

models.pyにCustomUserクラスを作成します。今回はデフォルトのユーザモデルからフィールドの追加はしないので直下にはpassを追加しておきます。独自のフィールドを設定する場合にはこちらにフィールドを追加してください。

users/models.py

from django.contrib.auth.models import AbstractUser
from django.db import models

class CustomUser(AbstractUser):
    pass

Custom User Formsの作成

users/forms.pyを作成し、以下のコードを追加します。

users/form.py

from django.contrib.auth import get_user_model
from django.contrib.auth.forms import UserCreationForm, UserChangeForm


class CustomUserCreationForm(UserCreationForm):

    class Meta(UserCreationForm):
        model = get_user_model()
        fields = ('email', 'username',)


class CustomUserChangeForm(UserChangeForm):

    class Meta(UserChangeForm):
        model = get_user_model()
        fields = ('email', 'username',)

Custom User Adminの設定

admin.pyに以下のコードを追加します。

users/admin.py

from django.contrib import admin
from django.contrib.auth import get_user_model
from django.contrib.auth.admin import UserAdmin

from .forms import CustomUserCreationForm, CustomUserChangeForm

CustomUser = get_user_model()


class CustomUserAdmin(UserAdmin):
    add_form = CustomUserCreationForm
    form = CustomUserChangeForm
    model = CustomUser
    list_display = ['email', 'username',]

admin.site.register(CustomUser, CustomUserAdmin)

これで設定は完了です。最後にマイグレーションを行いましょう。

docker-compose exec web python manage.py makemigrations users
docker-compose exec web python manage.py migrate

何をやるか？

超シンプルなドラクエ風ターン制バトルを作ってQ学習させてみます。

数％の確率でしか世界を救えない勇者くんを、Q学習で賢くすることが目的です。

なお、ゲーム部分・Q学習の実装については解説しますが、Q学習そのものは解説しません。
Q学習の詳しい理論を知りたい方は、こちらの良記事を一つずつ読んでいくと幸せになれます。

今さら聞けない強化学習（1）：状態価値関数とBellman方程式

読んでもらいたい人

• OpenAI Gymなどの既存のシミュレーション環境ではなく、自分でゲームを作って色々弄ってみたい方。
• Q学習の理論はなんとなく知ってるんだけど、「どうやって実装すれば良いかわからない！」という方。

ゲームをつくる

ルールはシンプルに、以下のように設計します。

• 勇者 vs 魔王の1対1
• 魔王のとる行動は「攻撃」のみ
• 勇者のとれる行動は「攻撃」と「回復」の2択
• 行動順序は、各キャラの素早さに一定の乱数を掛けてソートさせることで決定する

キャラクタークラスの実装

それでは早速、ゲーム本体の実装をしていきましょう。
まずはキャラクタークラスです。

dq_battle.py

class Character(object):

    """ キャラクタークラス"""

    ACTIONS = {0: "攻撃", 1: "回復"}

    def __init__(self, hp, max_hp, attack, defence, agillity, intelligence, name):
        self.hp = hp  # 現在のHP
        self.max_hp = max_hp  # 最大HP
        self.attack = attack  # 攻撃力
        self.defence = defence  # 防御力
        self.agillity = agillity  # 素早さ
        self.intelligence = intelligence  # 賢さ
        self.name = name  # キャラクター名

    # ステータス文字列を返す
    def get_status_s(self):
        return "[{}] HP:{}/{} ATK:{} DEF:{} AGI:{} INT:{}".format(
            self.name, self.hp, self.max_hp, self.attack, self.defence, self.agillity, self.intelligence)

    def action(self, target, action):

        # 攻撃
        if action == 0:

            # 攻撃力 - 防御力のダメージ計算
            damage = self.attack - target.defence
            draw_damage = damage  # ログ用

            # 相手の残りHPがダメージ量を下回っていたら、残りHPちょうどのダメージとする
            if target.hp < damage:
                damage = target.hp

            # ダメージを与える
            target.hp -= damage

            # 戦闘ログを返す
            return "{}は{}に{}のダメージを与えた".format(
                self.name, target.name, draw_damage)

        # 回復
        elif action == 1:

            # 回復量をINTの値とする
            heal_points = self.intelligence
            draw_heal_points = heal_points  # ログ用

            # 最大HPまで回復できるなら、最大HP - 現在のHPを回復量とする
            if self.hp + heal_points > self.max_hp:
                heal_points = self.max_hp - self.hp

            # 回復
            self.hp += heal_points

            # 戦闘ログを返す
            return "{}はHPを{}回復した".format(
                self.name, draw_heal_points)

今回のバトル設計はシンプルなので、プレイヤーと敵を区別することなく１つのクラスにまとめてしまっています。

各キャラクター（勇者と魔王）は、

• HP（体力）
• ATTACK(攻撃力)
• DEFENCE(防御力)
• AGILLITY（素早さ）
• INTELIGENCE（賢さ）

のステータスを持ちます。

「攻撃」でのダメージ計算は、
「（自分の攻撃力）ー（相手の防御力）」

という単純な式で計算しています。
また、「回復」コマンドでの回復量は賢さの数値そのままとしました。

バトル設計の全体像（状態遷移）

続いてバトル本体を実装していきます。
初めに、バトルの全体像（状態遷移）について理解しておく必要があります。

dq_battle.py

class GameState(Enum):

    """ ゲーム状態管理クラス"""
    TURN_START = auto()      # ターン開始
    COMMAND_SELECT = auto()  # コマンド選択
    TURN_NOW = auto()        # ターン中（各キャラ行動）
    TURN_END = auto()        # ターン終了
    GAME_END = auto()        # ゲーム終了

バトルには上記の通り、
「ターン開始」「コマンド選択」「ターン中」「ターン終了」「ゲーム終了」
の5つの状態があります。

状態遷移図で表すと下図のようになります。

このように、「ターン開始」状態から「ターン終了」状態までの遷移を、「ゲーム終了」状態になるまで（勇者か魔王のHPが0になるまで）延々とループさせるのがバトル設計の基本になります。

バトル本体の実装

それでは、バトル本体の実装です。
先にコード全体をみておきます。

dq_battle.py

class Game():

    """ ゲーム本体"""

    HERO_MAX_HP = 20
    MAOU_MAX_HP = 50

    def __init__(self):

        # キャラクターを生成
        self.hero = Character(
            Game.HERO_MAX_HP, Game.HERO_MAX_HP, 4, 1, 5, 7, "勇者")

        self.maou = Character(
            Game.MAOU_MAX_HP, Game.MAOU_MAX_HP, 5, 2, 6, 3, "魔王")

        # キャラクターリストに追加
        self.characters = []
        self.characters.append(self.hero)
        self.characters.append(self.maou)

        # 状態遷移用の変数を定義
        self.game_state = GameState.TURN_START

        # ターン数
        self.turn = 1

        # 戦闘ログを保存するための文字列
        self.log = ""

    # 1ターン毎にゲームを進める
    def step(self, action):

        # メインループ
        while (True):
            if self.game_state == GameState.TURN_START:
                self.__turn_start()
            elif self.game_state == GameState.COMMAND_SELECT:
                self.__command_select(action)  # 行動を渡す
            elif self.game_state == GameState.TURN_NOW:
                self.__turn_now()
            elif self.game_state == GameState.TURN_END:
                self.__turn_end()
                break  # ターン終了でもループを抜ける
            elif self.game_state == GameState.GAME_END:
                self.__game_end()
                break

        # ゲームが終了したかどうか
        done = False
        if self.game_state == GameState.GAME_END:
            done = True

        # 「状態s、報酬r、ゲームエンドかどうか」を返す
        return (self.hero.hp, self.maou.hp), self.reward, done

    # ゲームを1ターン目の状態に初期化
    def reset(self):
        self.__init__()
        return (self.hero.hp, self.maou.hp)

    # 戦闘ログを描画
    def draw(self):
        print(self.log, end="")

    def __turn_start(self):

        # 状態遷移
        self.game_state = GameState.COMMAND_SELECT

        # ログを初期化
        self.log = ""

        # 描画
        s = " *** ターン" + str(self.turn) + " ***"
        self.__save_log("\033[36m{}\033[0m".format(s))
        self.__save_log(self.hero.get_status_s())
        self.__save_log(self.maou.get_status_s())

    def __command_select(self, action):

        # 行動選択
        self.action = action

        # キャラクターを乱数0.5〜1.5の素早さ順にソートし、キューに格納
        self.character_que = deque(sorted(self.characters,
                                          key=lambda c: c.agillity*random.uniform(0.5, 1.5)))

        # 状態遷移
        self.game_state = GameState.TURN_NOW

        # ログ保存
        self.__save_log("コマンド選択 -> " + Character.ACTIONS[self.action])

    def __turn_now(self):

        # キャラクターキューから逐次行動
        if len(self.character_que) > 0:
            now_character = self.character_que.popleft()
            if now_character is self.hero:
                s = now_character.action(self.maou, self.action)
            elif now_character is self.maou:
                s = now_character.action(self.hero, action=0)  # 魔王は常に攻撃

            # ログを保存
            self.__save_log(s)

        # HPが0以下ならゲームエンド
        for c in self.characters:
            if c.hp <= 0:
                self.game_state = GameState.GAME_END
                return

        # 全員行動終了したらターンエンド
        if len(self.character_que) == 0:
            self.game_state = GameState.TURN_END
            return

    def __turn_end(self):

        # 報酬を設定
        self.reward = 0

        # キャラクターキューの初期化
        self.character_que = deque()

        # ターン経過
        self.turn += 1

        # 状態遷移
        self.game_state = GameState.TURN_START

    def __game_end(self):

        if self.hero.hp <= 0:
            self.__save_log("\033[31m{}\033[0m".format("勇者は死んでしまった"))
            self.reward = -1  # 報酬を設定
        elif self.maou.hp <= 0:
            self.__save_log("\033[32m{}\033[0m".format("魔王をやっつけた"))
            self.reward = 1  # 報酬を設定

        self.__save_log("-----ゲームエンド-----")

    def __save_log(self, s):
        self.log += s + "\n"

少々コードが長いですが、Q学習で重要な部分は2つだけです。

1つ目は、step()メソッドです。ここがバトルのメイン部分になります。

dq_battle.py

    # 1ターン毎にゲームを進める
    def step(self, action):

        # メインループ
        while (True):
            if self.game_state == GameState.TURN_START:
                self.__turn_start()
            elif self.game_state == GameState.COMMAND_SELECT:
                self.__command_select(action)  # 行動を渡す
            elif self.game_state == GameState.TURN_NOW:
                self.__turn_now()
            elif self.game_state == GameState.TURN_END:
                self.__turn_end()
                break  # ターン終了でもループを抜ける
            elif self.game_state == GameState.GAME_END:
                self.__game_end()
                break

        # ゲームが終了したかどうか
        done = False
        if self.game_state == GameState.GAME_END:
            done = True

        # 「状態s、報酬r、ゲームエンドかどうか」を返す
        return (self.hero.hp, self.maou.hp), self.reward, done

基本的には、前述した状態遷移図と処理の流れは同じです。

ただし、Q学習では1ターン毎に現在の状態を評価しなくてはならないので、「ゲーム終了」状態に限らず、「ターン終了」状態でもメインループを抜けなければなりません。

「ターン終了」状態において、Q学習をするために評価しなくてはならない変数は、

• 「状態s」
• 「報酬r」
• 「ゲームエンドかどうかを判定するフラグ」

の3つです。

ゲームエンドかどうかについては、単純に勇者のHPか魔王のHPが0になったかどうかで判断します。

状態sについては少し考える必要があります。
攻撃力や防御力などの複数のステータスがありますが、Q学習で評価すべきステータスは「勇者のHP」と「魔王のHP」の実質2つだけです。

今回のバトル設計では、攻撃力・防御力などの数値は常に一定なので、HP以外のステータスを評価する必要がないからです。逆に言えば、バフ・デバフなどでステータスが変化する場合はそれらの情報も必要になります。

報酬rについては、「ターン終了」と「ゲーム終了」状態それぞれで評価します。

dq_battle.py

    def __turn_end(self):

        # 報酬を設定
        self.reward = 0

    # （省略）

    def __game_end(self):

        if self.hero.hp <= 0:
            self.__save_log("\033[31m{}\033[0m".format("勇者は死んでしまった"))
            self.reward = -1  # 報酬を設定
        elif self.maou.hp <= 0:
            self.__save_log("\033[32m{}\033[0m".format("魔王をやっつけた"))
            self.reward = 1  # 報酬を設定

ターン経過による報酬は0としました。魔王を「最速で倒す」という目的を意識するなら、ターン経過での報酬を負の値にすればよいでしょう。（ただし、適切なパラメーターを設定するのは難しいですが。）

ゲーム終了時には、勇者が倒れてしまえば「-1」、魔王を倒せば「+1」の報酬を与えます。

2つ目に重要な部分は、reset()メソッドです。

dq_battle.py

    # ゲームを1ターン目の状態に初期化
    def reset(self):
        self.__init__()
        return (self.hero.hp, self.maou.hp)

単にゲームを初期化するだけのメソッドです。なお、Q学習のために初期状態を返す必要があります。

上記のstep()メソッドと合わせて、

ゲーム初期化(reset)→バトルが終了するまでターンを進める（step)→ゲーム初期化（reset）→バトルが終了するまでターンを進める（step)・・・

と、ゲームを繰り返すことで学習を進めていくことができます。

以上がQ学習をする上での、ゲームの根幹部分となります。

Q学習を実装する

エージェントクラスについて

Q学習は、エージェントクラス内で実装します。
エージェントとは、実際にゲームをするプレイヤーのようなクラスです。

エージェントはプレイヤー自身ですので、行動（攻撃か回復か）を選択したり、状態（勇者や魔王のHPなど）について知ったりすることはできますが、
ゲームの内部情報（行動順序を決める乱数など）を知ることはできません。

「行動」と、その行動によって得られた「状態」と「報酬」だけから学習を進めていくのが、
Q学習を含めた強化学習全般における基本的な理解になります。

初めに、エージェントクラスの全体を掲示しておきます。

q-learning.py

DIV_N = 10

class Agent:
    """エージェントクラス"""

    def __init__(self, epsilon=0.2):
        self.epsilon = epsilon
        self.Q = []

    # 方策をε-greedy法で定義
    def policy(self, s, actions):

        if np.random.random() < self.epsilon:

            # epsilonの確率でランダムに行動
            return np.random.randint(len(actions))

        else:

            # （Qに状態sが含まれており、かつそのときの状態におけるQ値が0でなければ）
            if s in self.Q and sum(self.Q[s]) != 0:

                # Q値が最大となるように行動
                return np.argmax(self.Q[s])
            else:
                return np.random.randint(len(actions))

    # 状態を数値に変換する
    def digitize_state(self, s):

        hero_hp, maou_hp = s

        # 勇者と魔王のHPをそれぞれDIV_Nで分割する
        s_digitize = [np.digitize(hero_hp, np.linspace(0, dq_battle.Game.HERO_MAX_HP, DIV_N + 1)[1:-1]),
                      np.digitize(maou_hp, np.linspace(0, dq_battle.Game.MAOU_MAX_HP, DIV_N + 1)[1:-1])]

        # DIV_Nの2乗までの状態数を返す
        return s_digitize[0] + s_digitize[1]*DIV_N

    # Q学習をする
    def learn(self, env, actions, episode_count=1000, gamma=0.9, learning_rate=0.1):

        self.Q = defaultdict(lambda: [0] * len(actions))

        # episode_countの分だけバトルする
        for e in range(episode_count):

            # ゲーム環境をリセット
            tmp_s = env.reset()

            # 現在の状態を数値に変換
            s = self.digitize_state(tmp_s)

            done = False

            # ゲームエンドになるまで行動を繰り返す
            while not done:

                # ε-greedy方策に従って行動を選択
                a = self.policy(s, actions)

                # ゲームを1ターン進め、その時の「状態、報酬、ゲームエンドかどうか」を返す
                tmp_s, reward, done = env.step(a)

                # 状態を数値に変換
                n_state = self.digitize_state(tmp_s)

                # 行動aによって得られた価値(gain) = 即時報酬 + 時間割引率 * 次の状態における最大のQ値
                gain = reward + gamma * max(self.Q[n_state])

                # 現在推測している（学習する前の）Q値
                estimated = self.Q[s][a]

                # 現在の推測値と、行動aを実行してみたときの実際の価値をもとに、Q値を更新
                self.Q[s][a] += learning_rate * (gain - estimated)

                # 現在の状態を次の状態へ
                s = n_state

状態を数値に変換する

エージェントクラスで少しわかりにくいのは、状態を数値に変換するメソッドでしょうか。

q-learning.py

    # 状態を数値に変換する
    def digitize_state(self, s):

        hero_hp, maou_hp = s

        # 勇者と魔王のHPをそれぞれDIV_Nに分割する
        s_digitize = [np.digitize(hero_hp, np.linspace(0, dq_battle.Game.HERO_MAX_HP, DIV_N + 1)[1:-1]),
                      np.digitize(maou_hp, np.linspace(0, dq_battle.Game.MAOU_MAX_HP, DIV_N + 1)[1:-1])]

        # DIV_Nの2乗までの状態数を返す
        return s_digitize[0] + s_digitize[1]*DIV_N

先ほども軽く触れましたが、Q学習をする上で評価すべき状態変数は、「勇者のHP」と「魔王のHP」の2つです。しかし、Q学習では状態を1つの数値として表す必要があります。つまり、以下のようなイメージです。

• 状態1：(勇者のHP, 魔王のHP) = (0, 0)
• 状態2：(勇者のHP, 魔王のHP) = (0, 1)
• 状態3：(勇者のHP, 魔王のHP) = (0, 2)

上記のように変換しても良いのですが、これだとHP×HPの数だけ状態が増えることになります。ドラクエじゃない方の某国民的RPGのように、HPが4桁もあると状態数が100万を超えて大変です（笑）。ですので、HPの割合に応じて状態を分割してやることにしましょう。

np.digitize(hero_hp, np.linspace(0, dq_battle.Game.HERO_MAX_HP, DIV_N + 1)[1:-1]

このコードをざっくりと解説すると、
np.linspace()で、0から最大HPまでをN分割し、
np.digitize()で、現在のHPが、何分割目に属しているかを返すイメージです。

今回はN=10としているので、

• HPが1割未満 → 0
• HPが1割以上、2割未満 → 1
• HPが2割以上、3割未満 → 2

のように変換してくれます。更に、

「勇者の状態（0〜9）+ 魔王の状態(0〜9)*10」
の計算をすることで、0から99までの100個に状態数を抑えることができます。

状態「15」なら、魔王のHPが「1」割未満、かつ勇者のHPが「5」割未満と、直感的にわかりますね。

方策の定義

方策はε-greedyを採用しています。

q-learning.py

    # 方策をε-greedy法で定義
    def policy(self, s, actions):

        if np.random.random() < self.epsilon:

            # epsilonの確率でランダムに行動
            return np.random.randint(len(actions))

        else:

            # （Qに状態sが含まれており、かつそのときの状態におけるQ値が0でなければ）
            if s in self.Q and sum(self.Q[s]) != 0:

                # Q値が最大となるように行動
                return np.argmax(self.Q[s])
            else:
                return np.random.randint(len(actions))

初学者の方のために簡単に解説すると、基本的には行動価値が最大となるように行動を決定し、εの確率でランダムな行動を採用する方策です。

行動にある程度のランダム性をもたせることにより、様々な行動を探索するため、Q値の初期値に依存することなく適切な学習が可能となります。

Q学習の実装

さて、ここまでくればQ学習に必要な変数・メソッドは全てそろいました。

Q学習のアルゴリズムは以下です。

1. $Q(s,a)$を初期化。
2. 任意の回数バトルを繰り返す：
   1. ゲーム環境の初期化
   2. ゲームエンドまでターンを進める：
      1. 方策$π$に従って行動$a$を選択する。
      2. 行動$a$を行い、報酬$r$と次の状態$s′$を観測する。
      3. $Q(s,a)$を以下のように更新する。
         $Q(s,a)$$\leftarrow$$Q(s,a)+α(r+γ*$$\underset{a′}{max}$$Q(s′,a′)−Q(s,a))$
      4. $s$$\leftarrow$$s′$とする。

記事の冒頭でも述べたとおり、Q学習の理論については解説しませんので、上記のアルゴリズムを素直に実装しましょう。

q-learning.py

    # Q学習をする
    def learn(self, env, actions, episode_count=1000, gamma=0.9, learning_rate=0.1):

        self.Q = defaultdict(lambda: [0] * len(actions))

        # episode_countの分だけバトルする
        for e in range(episode_count):

            # ゲーム環境をリセット
            tmp_s = env.reset()

            # 現在の状態を数値に変換
            s = self.digitize_state(tmp_s)

            done = False

            # ゲームエンドになるまで行動を繰り返す
            while not done:

                # ε-greedy方策に従って行動を選択
                a = self.policy(s, actions)

                # ゲームを1ターン進め、その時の「状態、報酬、ゲームエンドかどうか」を返す
                tmp_s, reward, done = env.step(a)

                # 状態を数値に変換
                n_state = self.digitize_state(tmp_s)

                # 行動aによって得られた価値(gain) = 即時報酬 + 時間割引率 * 次の状態における最大のQ値
                gain = reward + gamma * max(self.Q[n_state])

                # 現在推測している（学習する前の）Q値
                estimated = self.Q[s][a]

                # 現在の推測値と、行動aを実行してみたときの実際の価値をもとに、Q値を更新
                self.Q[s][a] += learning_rate * (gain - estimated)

                # 現在の状態を次の状態へ
                s = n_state

これで、ゲームとQ学習の実装までが完了しました。

ゲームを実行＆学習

ランダムに行動させてみる

Q学習をする前に、勇者の行動をランダムにしてバトルするとどうなるのか試してみましょう。

以下のコードを追加します。

q-learning.py

class Agent:

    # （省略）

    # テストバトル
    def test_run(self, env, actions, draw=True, episode_count=1000):

        turn_num = 0  # 撃破ターン数
        win_num = 0  # 勝数

        # episode_countの分だけバトルする
        for e in range(episode_count):

            tmp_s = env.reset()
            s = self.digitize_state(tmp_s)

            done = False

            while not done:
                a = self.policy(s, actions)
                n_state, _, done = env.step(a)
                s = self.digitize_state(n_state)
                if draw:
                    env.draw()  # バトルログを描画

            if env.maou.hp <= 0:
                win_num += 1
                turn_num += env.turn

        # 平均勝率・平均撃破ターン数を出力
        if not win_num == 0:
            print("平均勝率{:.2f}%".format(win_num*100/episode_count))
            print("平均撃破ターン数:{:.2f}".format(turn_num / win_num))
        else:
            print("平均勝率0%")


if __name__ == "__main__":

    game = dq_battle.Game()
    agent = Agent()

    actions = dq_battle.Character.ACTIONS

    """ 完全ランダムでバトル """
    agent.epsilon = 1.0
    agent.test_run(game, actions, episode_count=1000)

ε=1.0とすることで、100%完全ランダムで行動させています。
また、1000回のバトル結果から、平均勝率・平均撃破ターン数を計算するようにしてみました。

以下、実行結果です。

$ python q-learning.py 
平均勝率0.90%
平均撃破ターン数:64.89

勝率はかなり低いですね…。

ターン数を見ればわかるように、長期戦になりやすい傾向があります。
長期戦になればなるほど勇者が瀕死の状態が増えるので、その結果勝ちにくくなることが予想されます。

Q学習後にバトルさせてみる

以下のコードを追加します。

q-learning.py

if __name__ == "__main__":

    # （省略）

    """ Q学習する """
    agent.epsilon = 0.2
    agent.learn(game, actions, episode_count=1000)

    """ テストバトル """
    agent.epsilon = 0
    agent.test_run(game, actions, episode_count=1000)

ε=0.2とし、Q学習を実行してみます。

その後、1000回のテストバトルを行います。
なお、ε=0（0%ランダム）にすることで、学習した行動価値の通りに行動させています。

以下、学習するバトル数を変えて実行結果を示します。

実行結果（学習バトル数：50、テストバトル数：1000）

$ python q-learning.py 
平均勝率42.60%
平均撃破ターン数:56.19

実行結果（学習バトル数500、テストバトル数：1000）

$ python q-learning.py 
平均勝率100.00%
平均撃破ターン数:55.00

実行結果（学習バトル数5000、テストバトル数：1000）

$ python q-learning.py 
平均勝率100.00%
平均撃破ターン数:54.00

勝率100%になりましたね！

学習後のQ値ってどうなってるの？

少しだけ考察をしてみます。
学習した結果のQ値を見てみましょう。

以下、バトル数1000で学習させてみたときのQ値を、一部の状態について抜き出しました。

状態50：[-0.19, -0.1]
状態51：[-0.6623164987957537, -0.34788781183605283]
状態52：[-0.2711479211007827, 0.04936802595531123]
状態53：[-0.36097806076138395, 0.11066249745943924]
状態54：[-0.04065992616558749, 0.12416469852733954]
状態55：[0.17619052640036173, 0.09475948937059306]
状態56：[0.10659739434775867, 0.05112985778828942]
状態57：[0.1583472103200607, 0.016092008419030468]
状態58：[0.04964633744625512, 0.0020759614034820224]
状態59：[0.008345513895442138, 0.0]

状態の見方は、10の位が魔王の残りHP、1の位が勇者の残りHPです。
つまり上記の図は、魔王の残りHPが5割程度のとき、勇者の残りHPによって行動価値がどう変化するかを表現した図になります。

図から、勇者の残りHP（1の位）が低ければ「回復」コマンド、残りHPが高ければ「攻撃」コマンドを選択していることが読み取れます。

勇者の残りHPを固定したときのQ値も見ておきましょう。

状態07：[2.023809062133135, 0.009000000000000001]
状態17：[1.8092946131557912, 0.8310497919226313]
状態27：[0.8223927076749513, 0.5279685031058523]
状態37：[0.5565475393122992, 0.29257906153106145]
状態47：[0.25272081107828437, 0.26657637207739293]
状態57：[0.14094053800308323, 0.1533527340827757]
状態67：[0.0709128688771915, 0.07570873469406877]
状態77：[0.039059851207044236, 0.04408123679644829]
状態87：[0.023028972190011696, 0.02386492692407677]
状態97：[0.016992303227705185, 0.0075795064515745995]

上図は、勇者の残りHPが7割程度のとき、魔王の残りHPによって行動価値がどう変化するかを表現しています。
魔王の残りHPが少ないほど「攻撃」優位になっていることが読み取れると思います。

最後に

本記事は実装がメインなので、その他考察などは省略します。
余裕のある方は、ハイパーパラメータを変えて学習してみたり、バトルルールをより複雑にすることに挑戦してみると面白いでしょう。

また、筆者は強化学習初心者ですので、間違い等ありましたらお気軽にご指摘ください。
筆者の知識も強化されて嬉しくなります。

ソースはgithubにおいています。
https://github.com/nanoseeing/DQ_Q-learning

参考書籍

• 機械学習スタートアップシリーズ Pythonで学ぶ強化学習 [改訂第2版] 入門から実践まで
• つくりながら学ぶ! 深層強化学習 ~PyTorchによる実践プログラミング~

Python パッケージを Azure Artifacts フィードにアップロードしようとしてハマった話を一つ。

Azure Pipelines から Python パッケージをアップロードするには、

- task: TwineAuthenticate@1
  inputs:
    artifactFeed: FEED
- script: |
   python -m twine upload -r FEED --config-file $(PYPIRC_PATH) dist/*

のようなステップを書くと説明されているのですが、 Azure Artifacts に新たに作成するフィードは、これを書いている時点で、すべてプロジェクトスコープ フィード (project-scoped feeds) になっていて、

- task: TwineAuthenticate@1
  inputs:
    artifactFeed: PROJECT/FEED

とプロジェクト名を指定する必要があったわけです。これが分かるまで何度も試行錯誤を繰り返してしまいました。

現実的には、次のように変数参照にするのが便利だと思いますね。

- task: TwineAuthenticate@1
  inputs:
    artifactFeed: $(System.TeamProject)/FEED

Background

何もない状況から低コストでプログラミングを学ぶ方法を考えてみた。(Mac編)
の続きです。

今回はスマホは持っててPC買い始めのひと向けに無料でプログラムを学ぶ方法をまとめてみました。

OPAC - Library search system

おそらく、大卒が知っていて中卒・高卒が知らないことはOPACではないでしょうか。
大学であれば本・雑誌・論文などを探すのに一回生の最初のオリエンテーションで聞かされると思います。

Googleで
OPAC [住んでいる市町村区]
と検索すると図書館のリンク出てくるので、ここで
Python
と入力して探せば良いと思います。

予想ですが、プログラム自体何それおいしいの状態のエリアが大半であってもVBAが主力だと思います。そのため、少し大きめの図書館でなければ検索は引っかからないと思います。

ただし、三大都市圏に住んでいるのであればどこかの図書館へ行けばPythonの本は必ずあるのでこの方法で探す選択肢は持っていて損はないです。

試しに
東京都立図書館
で東京都内にあるPythonの本を探してみました。

ここでは、[区立図書館][市町村立図書館]にチェックして"Python"で検索しました。

50音順なので足立区から表示されますがひとつの図書館に数十冊はありそうです。
左側のメニューから近くの図書館をクリックすると検索結果が反映されます。

東京都内だと約2500 ~ 3000冊はありそうです。

Web site

• python-izm - 基礎編
• Let'sプログラミング - Python入門

の２つのサイトに書かれているコードを使って実行すれば大体のことは網羅できます。
(説明雑でサーセン)

あとは、公式サイトのドキュメントにもサンプルコードが書かれています。
ライブラリはたくさんあるわけですが、学習が進んで何か自分でアプリを作りたい時に必要そうな分だけ選んで使えば良くて、全て覚えている開発者はほとんどいないんじゃないかと思います。

• Python2系
  • チュートリアル
  • ライブラリ
• Python3系
  • チュートリアル
  • ライブラリ

Skill Check

ひと通り学習が終わったら力試しとしてスキルテストすると良いです。
本格的に始めたい場合はAOJでもいいのですが、難易度の設定バランスや問題の出題方式だとPaizaの方がとっつきやすいです。S・A・B・C・Dと難易度別に問題が約200問あるのですが、C・Dの問題を多く正解できれば良いと思います。B問題ができると尚良いです。
（最低限一般の開発者が求めているスキルはB問題でしょうか。実際問題自分がなんとか一定の給与をもらって働いている。)

自分はと言いますと一回しかA問題は解けんかったです
授業で再帰・２分探索・動的計画法などを学んでアウトプットできていればとちょっと後悔しています

(注： ギノの社員でも契約社員でもないとです)

NextStage

Skill Checkで自信がついたらゲーム・webサイトを作成するなどのフェーズに移ればいいです。
freegamesかpygameを使ってサンプルコードを改造するところから始めるとかいいかもしれないです。
(以前に[Python]pipからゲームをインストールして遊んでみたで方法をまとめてます。)

PostScript

次はPython2系とPython3系の話をしようと思います。

Lambdaってなに？

AWSの提供するサービスの一つで、サーバ管理をすることなく処理を実行することができるサービスです。

例えばLambdaを利用しない場合、EC2でサーバを作成し処理実行に必要なミドルウェアや言語をインストールし、環境設定をすることで初めて処理が実行されます。

しかし、Lambdaではそのようなサーバのプロビジョニングをすることなく、処理を書くだけで実行することが可能になります。

これにより、
★サーバ自体の管理・メンテナンスが不要
★処理が実行されている時間のみ課金されるので処理の頻度によってはコストを大幅削減可能
★他AWSサービスと容易に連携できるので、AWSメインのアーキテクチャでは大変便利
というメリットがあります。

トリガーとは？

書いた処理を実行させるにはトリガーというものを設定する必要があります。
AWSのあらゆるサービスと連携ができ、CloudWatchでアラームが出たら実行、kinesis data stream にデータが存在するとき実行、S3にファイルが置かれたら実行などが代表的な例になります。

つまりLambdaは実行する条件と実行される処理だけを書けば動いてくれるサービスなのです。

実際に設定したい！

ではs3にファイルをアップロードし、アップロードされたことを確認する処理を書いてみましょう。

1.AWSコンソールにログインし、Lambdaの画面で「関数の作成」をクリック

Lambdaはリージョンごとに関数を設定することになります。
S3の場合は問題ないですが、他リージョンの影響を受けるサービスと連携するときは同じリージョンを選択しましょう。

2.「1から作成」を選択し、「関数名」と「ランタイム」、「実行ロール」を設定

今回はpython 3.7を利用します。
ロールは「基本的なLambdaアクセス権限で新しいロールを作成」にします。
既存のロールを利用でも問題ないです。
その場合、Lambdaがトリガー条件のサービスやLog出力用のCloudWatchへのロググループ書き込み権限がある必要があります。
※新規作成の場合はCloudWatchのロググループ書き込み権限は自動でアタッチされます。

「ランタイム」とは処理の言語を意味します。2020/02/05時点では選択できる言語は下記となっております。
・Java 11/8
・.NET Core2.1(C#/PowerShell)
・GO 1.x
・Node.js 12.x/10.x
・python 3.8/3.7/3.6/2.7
※カスタムランタイムを利用することであらゆる言語を利用できるようになりますが、本記事では省略します。

3.「トリガーを追加」で実行条件を設定

今回はs3に何かオブジェクトが作成されたときに実行できるようにします。
トリガー媒体をs3に設定したときの項目は下記の5つになります。
・バケット
・イベントタイプ(ファイルがPUTされたら・・・、ファイルが削除されたら・・・など)
・プレフィックス(バケット以降のディレクトリパスやファイル名など)
・サフィックス(ファイル名や拡張子など)
・トリガーの有効化(チェックするとすぐ上記設定のトリガーが動きます。処理のテストが終わってからonにしましょう)

4.関数コードを設定

デフォルトではlambda_function.pyというファイルに処理が記述されています。

lambda_function.py

import json

def lambda_handler(event, context):
    # TODO implement
    return {
        'statusCode': 200,
        'body': json.dumps('Hello from Lambda!')
    }

lambda_handlerとはLambdaに設定されているトリガー条件が満たされたときに自動的に実行される関数です。
# TODO implement 以降を編集し、実行したい処理を書きましょう。

今回はs3にファイルが置かれたら置かれたパスとファイル名を表示するようにしてみます。

lambda_function.py

import json

def lambda_handler(event, context):
    bucket = event['Records'][0]['s3']['bucket']['name']
    key = event['Records'][0]['s3']['object']['key']
    print(bucket+'バケットに'+key+'が作成されました！')

lambda_handlerに引数で渡しているeventにはトリガーとなったイベント情報が配列で入っており、
上記はそこからバケット名とファイルパスを取得し、表示を行っています。

5.実行・確認

関数の実行自体はトリガーで指定したパスにファイルを置くことで実行できます。
それ以外にも下画像の「テスト」からeventに入る値をjson形式で設定することにより疑似的にLambda関数を実行することが可能です。

また、関数の実行LogについてはCloudWatchのロググループから対象の関数を選択し、そこから確認することが可能です。

正しくLambda関数が実行できていることが確認できました！！

終わりに

LambdaはAWS環境を利用しており、シンプルな処理をおこなうのであればコスト面、管理面でもとてもおススメです。
少し気になったのは関数が何をしているかの説明文を書いておく欄がないので、わかりやすい関数名・処理内でのコメント記載を徹底し、快適なLambdaライフをおくりましょう！！