ASE (Atomic Simulation Environment)とは

量子化学計算や、古典分子動力学法をpython上で実行できるライブラリです。以下のサイトを参照していただければ、どのようなコンセプトで作られたか理解できると思います。
https://wiki.fysik.dtu.dk/ase/about.html

インストールのしかた

以下のURLを参考にしました：
https://wiki.fysik.dtu.dk/ase/install.html
私は、git clone -b 3.18.1 https://gitlab.com/ase/ase.gitで問題なくできました。

古典分子動力学をやってみる

まず初めに物理系を作成する

テスト系として窒素分子（N2）を作成してみます。ASEは基本的にase.Atomsに物理系の座標情報を入れていきます。

from ase import Atoms
d = 1.1
molecule = Atoms('2N', positions=[(0., 0., 0.), (0., 0., d)])

moleculeにちゃんと座標データが入っているか確認してみましょう。そのためには、ase.ioを使います：

from ase import io
write('N2.pdb',molecule)

（.pdbと出力ファイルを指定するとPDBファイル形式を取り扱える。これは生物系の人には便利）
すると以下のような構造が得られます（PyMOLか何かで確認）。確かに窒素の２原子分子が作成されているようです。

分子動力学計算の設定

以下の設定で計算を走らせます
- 積分器：Velocity Verlet
- ポテンシャル：Effective medium theory (EMT)、(https://wiki.fysik.dtu.dk/ase/ase/calculators/emt.html#module-ase.calculators.emt)
- 時間刻み：0.01 fs
- 20 stepsごとに構造の保存

from ase.md.verlet import VelocityVerlet
from ase import units
from ase.calculators.emt import EMT
molecule.set_calculator(EMT())
dyn = VelocityVerlet(molecule, dt=0.01 * units.fs)
for i in range(100):
     pot = molecule.get_potential_energy()
     kin = molecule.get_kinetic_energy()
     print('%2d: %.5f eV, %.5f eV, %.5f eV' % (i, pot + kin, pot, kin))
     H.append(pot + kin)
     K.append(kin)
     U.append(pot)
     dyn.run(steps=20)
     write(f'traj.pdb',molecule,append=True)

# エネルギーの図作成
import matplotlib.pyplot as plt
%matplotlib inline
plt.plot(H, c='black')
plt.plot(K, c='red')
plt.plot(U, c='blue')
plt.savefig('energy.png')

得られたトラジェクトリのエネルギーの確認

下図をみていただくと、たしかに、全エネルギーは保存しており、NVEアンサンブルになっているようです（黒線が全エネルギー、赤線が運動エネルギー、青線がポテンシャルエネルギー）。

感想

かなり簡単に計算結果を得ることができました。オブジェクト指向で書かれているのでソフトウェアの設計の勉強にもなりました。このような形でとある過去の遺産をラップできたらいいのになあ。とりあえず基本的な使い方はりかいしましたが、グラフェンの計算など気になることがあるので、もう少し遊んでみようと思います。

はじめに

PythonもOpenCVも業務未経験でした。
Pythonは知ってましたけど、OpenCVなんて名前初めて聞いたくらいでした。笑

ネット上から色々情報は拾えて来れましたが、良い記事見つけては試しての繰り返しで作成してたのですが

「この処理は何の記事を参考にしたんだっけ？」

とせっかく調べたのに忘れてしまうのはもったいないと思って参考になったサイトを纏めました。

python + opencvの基本的なこと

1. 公式サイト

2. 画像処理入門講座 : OpenCVとPythonで始める画像処理

3. 「詳解OpenCV3」の備忘録(C++)

4. 画像処理の数式を見て石になった時のための、金の針

---

画像の明るさ調整に関すること

ガンマ補正

1. ガンマ補正のうんちく

2. ガンマ補正で画像を見やすく調整

3. LUTを使って画像の暗部を持ち上げてみる

---

画像のノイズ処理に関すること

ぼかし処理

1. 公式サイト(GaussianBlur)

2. 【Python】OpenCVで画像のぼかしと平滑化

クラスタリング、フィルタリング処理

1. 公式サイト(K-Meansクラスタリング)

2. Mean Shift Filtering

その他の処理

1. 公式サイト(Non-local Means Denoising)

2. python+OpenCVでエッジを保存した平滑化(BilateralFilter, NLMeansFilter)

3. 【Python/OpenCV】写真のアニメ絵化

---

2値化画像に関すること

thresholdを使って閾値で2値化

1. 公式サイト

2. OpenCV 画像の二値化

エッジを検出して2値化

1. 公式サイト

2. OpenCVとVisual C++による画像処理と認識（７）

3. Python OpenCVで人物画像の背景を消す

モルフォロジー変換

1. 公式サイト

2. OpenCV - モルフォロジー演算について

背景差分

1. 公式サイト

2. 背景差分を求める

3. 使って画像の背景を削除する

4. iQONでクロールしたアイテム画像がコーディネートに使われるまで

5. 背景差分に関する備忘録

6. 金属製品のメッキ不良を自動で検出するためには？

7. OpenCVで背景差分をやってみた。

8. オーバーフローせずに画像の明るさを増やす

---

画像のセンタリングに関すること

1. OpenCV - findContours() による輪郭抽出

2. OpenCVで長方形の画像に余白をつけて正方形にする

3. ホワイトボードの画像からポストイットを検出する

4. OpenCV - 画像をリサイズする方法 (cv2.resize)

5. pythonで画像のアスペクト比を維持したままサイズ変更して余白を埋める

---

GrabCutに関すること

1. 公式サイト

2. 二次元画像を「映え」させる【GrabCut】

3. 【OpenCV】GrabCutの使い方【セグメンテーション】

4. PythonでOpenCV~GrabCut(前景抽出)~

5. Python と OpenCV によるグラフカット処理

---

その他

1. Watershedアルゴリズムを使った画像の領域分割

2. 物体セグメンテーションアルゴリズム"watershed"を詳しく

3. 【Python】OpenCVでWatershedアルゴリズムと物体のセグメンテーション

4. OpenCVとPythonで画像の自動補間（Fast Marching Method, Navier-Stokes）

5. OpenCVで鮮鋭化

6. そこそこな線画を目指す OpenCV

$$
\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}}
$$

はじめに

量子情報理論の分野においては、環境系と注目系のように複合した全体系の中で、各々の系がどう変化していくとか、どんな性質をもっているとかいう観点で、いろいろな研究がなされてきました。そういった研究を進める中で、有用なツール・知見の蓄積がなされてきたわけですが、今回は、その代表選手の一つである「シュミット分解(Schmidt decomposition)」について勉強してみます。だいたいわかってきたところで、量子計算シミュレータqlazyを使って、そこから導き出される性質を確認しつつ、その面白さを噛み締めてみたいと思います。

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

1. ニールセン、チャン「量子コンピュータと量子通信（２）」オーム社（2005年）
2. 石坂、小川、河内、木村、林「量子情報科学入門」共立出版（2012年）
3. 富田「量子情報工学」森北出版（2017年）
4. シュミット分解と量子状態の純粋化・エンタングルメントとの関係 - 物理とか

シュミット分解とは

部分系A,Bからなる全体系ABにおける、ある純粋状態を$\ket{\Phi}^{AB}$とすると、

\ket{\Phi}^{AB} = \sum_{k=1}^{R} \lambda_{k} \ket{\phi_{k}}^{A} \ket{\psi_{k}}^{B}  \tag{1}

を満たす、負でない実数$\{ \lambda_{k} \}$、および系Aと系Bにおける正規直交基底$\{ \ket{\phi_k}^{A} \}$、$\{ \ket{\psi_{k}}^{B} \}$が必ず存在します。この表現のことを「シュミット分解」と言い、$\{ \lambda_{k}\}$のことを「シュミット係数(Schmidt coefficient)」、$\{ \lambda_{k}\}$の中のゼロでないものの個数を「シュミット・ランク(Schmidt rank)」、正規直交基底$\{\ket{\phi_{k}}^{A}\}$、$\{\ket{\psi_{k}}^{B}\}$のことを、系A、系Bに対する「シュミット基底(Schmidt basis)」と言います。

ここで、$\ket{\Phi}^{AB}$の正規化条件から、

\sum_{k=1}^{R} \lambda_k^{2} = 1  \tag{2}

が成り立ちます。

それでは、任意の$\ket{\Phi}^{AB}$が式(1)のように分解できることを証明してみます。

【証明】

$\ket{\Phi}^{AB}$は、ABの純粋状態なので、系Aの正規直交基底$\{ \ket{i}^A\}$と系Bの正規直交基底$\{ \ket{j}^B\}$を使って、一般に以下のように書けます¹。

\ket{\Phi}^{AB} = \sum_{i=1}^{M} \sum_{j=1}^{N} c_{ij} \ket{i}^{A} \ket{j}^{B}  \tag{3}

ここから出発して、式(1)を導ければ証明完了なのですが、それには特異値分解を使います。

いま、式(1)の$c_{ij}$をM行N列の行列$C$の($i,j$)要素とみなすと、M行M列のユニタリ行列$U$、N行N列のユニタリ行列$V$、M行N列の対角行列$D$を用いて、

C = UDV  \tag{4}

のように特異値分解できます²。つまり、

c_{ij} = \sum_{k=1}^{R} u_{ik} d_{kk} v_{kj}  \tag{5}

です。ここで、$u_{ik}, \space d_{kk}, \space v_{kj}$は各々$U,D,V$の要素です。また、$R$は、$D$の対角成分に並ぶゼロでない成分の個数を表します。

式(5)を式(3)に代入します。

\ket{\Phi}^{AB} = \sum_{i=1}^{M} \sum_{j=1}^{N} \sum_{k=1}^{R} u_{ik} d_{kk} v_{kj} \ket{i}^{A} \ket{j}^B  \tag{6}

ここで、

\begin{align}
\ket{\phi_{k}}^{A} &= \sum_{i=1}^{M} u_{ik} \ket{i}^{A} \\
\ket{\psi_{k}}^{B} &= \sum_{j=1}^{N} v_{kj} \ket{j}^{A} \tag{7}
\end{align}

とおくと、

\ket{\Phi}^{AB} = \sum_{k=1}^{R} \lambda_{k} \ket{\phi_{k}}^{A} \ket{\psi_{k}}^{B}  \tag{1}

となり、式(1)が導けました。また、式(7)より、$\{ \ket{\phi_{k}}^{A} \}$および$\{\ket{\psi_{k}}^{B}\}$は正規直交基底に対するユニタリ変換なので、両方とも正規直交基底です。（証明終）

さて、ここで、出発点とした式(3)とシュミット分解の式(1)を改めて眺めてください。証明の過程で何をやったかというと、要はユニタリ変換で部分系Aと部分系Bの基底を変換しているだけです。どのように変換したかというと、考えるべき基底の数をなるべく少なくなるようにした、ということになります³。と同時に、部分系Aと部分系Bの間のエンタングルメント（もつれ）が最大に見えるように変換したとも言えます。とすると、シュミットランク$R$というのは、ある意味、エンタングルメントの大きさを表す指標と考えることができます。シュミットランクが１というのは、エンタングルメントがない状況、つまり、部分系AとBが積状態（＝テンソル積として分割できる状態）になっているということを表しており、シュミットランクが２以上というのは、両者エンタングルしている状況ということです。

部分系の性質

全体系のある部分系をA、残りの部分系をBとすると、全体系で定義されたどんな純粋状態も、式(1)のようにシュミット分解の表式に書き直すことができる、ということがわかりました。この知見を利用すると、部分系に関して以下の性質が導けます。

密度演算子の固有値

全体系の密度演算子$\rho_{AB}=\ket{\Phi}^{AB} \bra{\Phi}^{AB}$の部分トレースをとると、部分系の密度演算子$\rho_{A}, \rho_{B}$は、

\begin{align}
\rho_{A} &= Tr_{B} (\ket{\Phi}^{AB} \bra{\Phi}^{AB}) = \sum_{k=1}^{R} \lambda_{k}^{2} \ket{\phi_{k}}^{A} \bra{\phi_{k}}^{A}  \\
\rho_{B} &= Tr_{A} (\ket{\Phi}^{AB} \bra{\Phi}^{AB}) = \sum_{k=1}^{R} \lambda_{k}^{2} \ket{\psi_{k}}^{B} \bra{\psi_{k}}^{B}  \tag{8}
\end{align}

となります。したがって、$\rho_{A}$の固有値も$\rho_{B}$の固有値も$\{ \lambda_{k}^2 \}$で、同じということになります。また、固有値（＝ユニタリ変換で対角化したときの対角成分）が同じということは、当然ランクも同じで、さらに、その２乗トレースも等しいということになります。つまり、

Tr(\rho_{A}^{2}) = Tr(\rho_{B}^{2})  \tag{9}

です。

このことは、全体系をどのように分割したとしても成り立つわけで、とても面白い性質だと思います。どっちの部分系で固有値を計算したとしても全く同じになるということなので、この固有値というのは、何かAとBとの間の分割面というか、関係性を特徴づけるような指標に思えてきます。そういえば、この固有値（の平方根）はシュミット係数そのものでした。つまり、その関係性というのは、要はエンタングルメントのことではないでしょうか、ということが、じんわりとわかってきます。

また、ある純粋状態に対するシュミット係数（あるいはシュッミットランク）を知りたいとすると、想定している部分系（どっちの部分系でも良い）で部分トレースをとって、その密度演算子の固有値を求めれば良い、ということもわかってきます（一つの計算方法として）。

ユニタリの自由度

部分系Aで定義される局所的なユニタリ変換を$U_{A}$、部分系Bで定義される局所的なユニタリ変換を$U_{B}$とすると、それによって全体系の純粋状態$\ket{\Phi}^{AB}$は、

\begin{align}
(U_{A} \otimes U_{B}) \ket{\Phi}^{AB} &= (U_{A} \otimes U_{B}) \sum_{k=1}^{R} \lambda_{k} \ket{\phi_{k}}^{A} \ket{\psi_{k}}^{B} \\
&= \sum_{k=1}^{R} \lambda_{k} (U_{A} \ket{\phi_{k}}^{A}) (U_{B} \ket{\psi_{k}}^{B}) \\
&= \sum_{k=1}^{R} \lambda_{k} \ket{\phi_{k}^{\prime}}^{A} \ket{\psi_{k}^{\prime}}^{B} \tag{10}
\end{align}

のように変換されます。ここで、

\begin{align}
\ket{\phi_{k}^{\prime}}^{A} &= U_{A} \ket{\phi_{k}}^{A} \\
\ket{\psi_{k}^{\prime}}^{B} &= U_{B} \ket{\psi_{k}}^{B} \tag{11}
\end{align}

とおきました。

これより、局所的なユニタリ変換で、シュミット分解の性質（シュミット係数やシュミットランク）は不変であるということがわかります。つまり、部分系AとBの間に存在するエンタングルメントの性質は、局所ユニタリ変換で不変に保たれるということです。こういう不変性は、いろいろ理論計算をする上で、何かと便利に使えるものなので、覚えておくと良いと思います。

シミュレータで確認

それでは、上で述べてきた知見の一部をシミュレータで確認してみたいと思います。全体系の任意の純粋状態に対する部分系の密度演算子の固有値が、全体系をどういう形に２分割したとしても、同じということだったので、本当にそうなるかを確認します。

全体のPythonコードは以下です。

import random
import numpy as np
from scipy.stats import unitary_group
from scipy.linalg import eigh
from qlazypy import QState, DensOp

MIN_DOUBLE = 0.000001

def random_qstate(qubit_num):

    dim = 2**qubit_num
    vec_ini = np.array([0.0]*dim)
    vec_ini[0] = 1.0
    mat = unitary_group.rvs(dim)
    vec = np.dot(mat, vec_ini)
    qs = QState(vector=vec)

    return qs

def random_qubit_id(qubit_num):

    id = list(range(qubit_num))
    random.shuffle(id)
    qubit_num_A = random.randint(1,qubit_num-1)
    qubit_num_B = qubit_num - qubit_num_A
    id_A = id[:qubit_num_A]
    id_B = id[qubit_num_A:]

    return id_A,id_B

def eigen_values(densop):

    matrix = densop.get_elm()
    eigvals = eigh(matrix, eigvals_only=True)
    eigvals_out = [eigvals[i] for i in range(len(eigvals)) if eigvals[i] > MIN_DOUBLE]

    return eigvals_out

if __name__ == '__main__':

    # whole quantum state
    qubit_num = 5
    qs = random_qstate(qubit_num)
    de = DensOp(qstate=[qs], prob=[1.0])  # pure state

    # partial density operators (system A and B)
    id_A, id_B = random_qubit_id(qubit_num)
    de_A = de.patrace(id=id_B)
    de_B = de.patrace(id=id_A)

    # eigen-values of density operators (system A and B)
    eval_A = eigen_values(de_A)
    eval_B = eigen_values(de_B)

    print("== system A ==")
    print("- qubit id =", id_A)
    print("- square trace = ", de_A.sqtrace())
    print("- eigen values =", eval_A)
    print("- rank =", len(eval_A))

    print("== system B ==")
    print("- qubit id =", id_B)
    print("- square trace = ", de_B.sqtrace())
    print("- eigen values =", eval_B)
    print("- rank =", len(eval_B))

    qs.free()
    de.free()
    de_A.free()
    de_B.free()

何をやっているか、順に説明します。

# whole quantum state
qubit_num = 5
qs = random_qstate(qubit_num)
de = DensOp(qstate=[qs], prob=[1.0])  # pure state

全体系の量子ビット数を5にして、量子状態をランダムに設定します。関数random_qstateの中で実行しているので、その中身を見てみます。

def random_qstate(qubit_num):

    dim = 2**qubit_num
    vec_ini = np.array([0.0]*dim)
    vec_ini[0] = 1.0
    mat = unitary_group.rvs(dim)
    vec = np.dot(mat, vec_ini)
    qs = QState(vector=vec)

    return qs

指定された量子ビット数に対応した状態の次元数をdimとして、numpyの配列（ベクトル）を適当に決めます。第0成分のみ1.0、それ以外を0.0としました。ここで、ランダムな量子状態にするために、scipyの関数を召喚します。

from scipy.stats import unitary_group
...

mat = unitary_group.rvs(dim)

です。これは、ランダムなユニタリ行列を生成するための関数です。ここにドキュメントがあります⁴。生成されたユニタリ行列をvec_iniに適用して、できたランダムベクトルに基づき、量子状態QStateクラスのインスタンスを作り、リターンします。

もとのmain部に戻ります。できあがったランダム量子状態を用いて、密度演算子DensOpクラスのインスタンスを生成します。これで、全体系に対する純粋状態の密度演算子ができあがりました。

# partial density operators (system A and B)
id_A, id_B = random_qubit_id(qubit_num)
de_A = de.patrace(id=id_B)
de_B = de.patrace(id=id_A)

全体系の量子ビット（5量子ビット）を部分系Aと部分系Bにランダムに分割して、各々の量子ビット番号リストを変数id_A,id_Bに格納しています。関数random_qubit_idで実行しています。説明は省きます。関数定義を見てください。この量子ビット番号リストを使って、全体系の密度演算子の部分トレースをとります。結果を各々、変数de_A,de_Bに格納します。

# eigen-values of density operators (system A and B)
eval_A = eigen_values(de_A)
eval_B = eigen_values(de_B)

固有値問題を解いて固有値のリストを求めます。関数eigen_valuesの中で、scipyの関数を使って実行しています。eighはエルミート行列に対する固有値を求める関数です。

from scipy.linalg import eigh
 ...

def eigen_values(densop):

    matrix = densop.get_elm()
    eigvals = eigh(matrix, eigvals_only=True)
    eigvals_out = [eigvals[i] for i in range(len(eigvals)) if eigvals[i] > MIN_DOUBLE]

    return eigvals_out

eigen_values関数の最初の行で、密度演算子DensOpクラスのget_elm()メソッドを使って行列表現（numpyの２次元配列）を取得しています。それを使って、固有値問題を解いています。

メイン部に戻ります。

print("== system A ==")
print("- qubit id =", id_A)
print("- square trace = ", de_A.sqtrace())
print("- eigen values =", eval_A)
print("- rank =", len(eval_A))
 ...

部分系Aの量子ビット番号、２乗トレース、固有値、ランクを表示します。部分系Bも同様に表示します。

さて、実行結果です。以下のようになりました。

== system A ==
- qubit id = [4, 1, 2]
- square trace =  0.38437171
- eigen values = [0.0822181470819045, 0.1464507832125754, 0.2143642754307335, 0.5569667940374802]
- rank = 4
== system B ==
- qubit id = [0, 3]
- square trace =  0.38437171
- eigen values = [0.08221814299667664, 0.14645078689065036, 0.2143642733317681, 0.5569667867809049]
- rank = 4

部分系Aに3量子ビット、部分系Bに2量子ビットが各々割り当てられました。ランクは4で固有値を見ていただくとわかる通り、きちんと一致しました！⁵また、２乗トレースは１以下なので、言わずもがなですが混合状態になっています。

もう一回やってみます。

== system A ==
- qubit id = [2, 3, 0, 4]
- square trace =  0.53145101
- eigen values = [0.3745986253398182, 0.6254013753462571]
- rank = 2
== system B ==
- qubit id = [1]
- square trace =  0.53145101
- eigen values = [0.37459862198601285, 0.6254013780139871]
- rank = 2

というわけで、今度は、4量子ビットと1量子ビットに分割されました。ランクは2となり、この場合も固有値は一致しました！全体系の量子ビット数も変えつつ、何度も実行してみましたが、すべての場合で固有値は一致しました！

これは、それほど自明なことではないので、面白がっていただけると幸いです。純粋状態の密度演算子があって、それを２つの部分系に分けたとすると、必ず成り立つ性質です。

おわりに

今回は、純粋状態を部分系にわけるお話でした。その部分系は一般に混合状態になります。では、逆に混合状態に何らかの部分系を追加することで、純粋状態にできるでしょうか。答えはイエスで、それを実行する手続きのことを「純粋化」と言います。今回の記事で一緒に説明しようと思ったのですが、長くなりました。というわけで、次回は「純粋化」の予定です。

以上

---

1. これは大丈夫ですよね。任意の状態は適当な基底を使った重ね合わせで表現できます、と言っているだけです。 ↩

2. 特異値分解は、通常$C = UDV^{\dagger}$と書かれたりします。添字を使って式展開するときに、ダガーがあると複素共役とか添字の順番を逆にしないと、、とかちょっと余計な注意をしないといけないので、$V^{\dagger} \rightarrow V$と置き換えました。 ↩

3. データサンプルの分布特徴を利用して、大量の説明変数をなるべく減らすという統計分析の世界でよく使われる「主成分分析」に似ています。 ↩

4. 中で何をやっているか理解できていませんが（汗）、Haar測度に基づく由緒正しいランダムユニタリ行列を計算しているようです。今回は単に適当な量子状態を作りたかっただけなので、ここまでちゃんとしたランダム性は必要ないですが、ランダムな量子状態をどうやったら作れるんだろうとちょっと調べていたら、ひっかかってきたので試しに使ってみました。採用している手法は、以下の論文に記載されているようです。参考まで。How to generate random matrices from the classical compact groups ↩

5. 小数点以下７桁以降は微妙に違っていたりしますが、計算誤差なので無視します ↩

自然言語処理で遊ぶための学習用データとして魅力的なWikipediaのデータですが、ダウンロード可能なDumpデータは3GB程度と容量が大きくXMLデータだったりして扱うのが大変です。もしDBに格納されていれば条件検索して記事数など確認しながら必要なデータだけ取り出せて便利です。今回はファイルとして持ち運びも可能なSQLiteでWikipediaの記事DBを作成します。

※※ 実行できることは確認していますが、処理時間が長いため記事作成に当たって再確認せずに記憶を頼りに書いています。ディレクトリの指定など違った場合はご一報ください※※

目的

・Wikipediaの記事データをDB(SQLite)に格納する
・記事名でパターン検索して、タイトルに特定の文字を含む記事データを抽出してみる。

方法(概要)

・Wikipadiaのダンプデータをダウンロードする
・Wikiextractorでテキストを整形しつつjsonデータとして抽出する
・抽出したjsonデータを拙作のPythonスクリプトwikiextractor2sqliteを用いてsqliteのDBファイルに格納する

wikiextractor2sqlite
簡単なPythonスクリプトなので、Qiitaへ貼るだけでも良かったのですがせっかくなので

実行環境

• Ubuntu 18.04
• Python 3.6.0
• sqlite3 3.22.0

手順

記事データのダウンロード

次のURLから、pages-articleデータをダウンロードします。
https://dumps.wikimedia.org/jawiki/

latestでも良いですが、後で別のデータと結合したい場合など日付付きの方が安全かもしれません。
今回は次のファイルをダウンロードしました。

wget https://dumps.wikimedia.org/jawiki/20191001/jawiki-20191001-pages-articles.xml.bz2

完了までだいたい25分くらいでした。

wikiextractorをダウンロード

次のコマンドでダウンロードします

git clone https://github.com/attardi/wikiextractor

Wikiextractorでテキストを整形しつつjsonデータとして抽出

次のコマンドでWikiextractorを実行します。

./wikiextractor/WikiExtractor.py --json -q -o extracted jawiki-latest-pages-articles.xml.bz2

正確に覚えていませんが、だいたい30分くらいだったと思います。
処理が完了するとextractedというフォルダ内にAAといったアルファベット2文字のフォルダが作成され、その中に1MB毎に分割された記事データのJsonファイルが作成されます。
いくつもWARNINGが出ますが、対象の記事データも一応出力はされているようです。

オプションの説明

• --json: 出力をタグ形式からJSON形式に変更
• -q: メッセージの量を減らす
• -o extracted: 出力先として「extracted」フォルダを指定

Wikiextractorの詳細はこちらから
Wikiextractor

wikiextractor2sqlite をダウンロード

次のコマンドでダウンロードします。
合わせて必要ライブラリをインストールします。

git clone https://github.com/yuukimiyo/wikiextractor2sqlite.git
pip install tqdm

Wikipedia DBを作成

./wikiextractor2sqlite/wikiextractor2sqlite.py ./extracted -o wikipedia.db

オプションの説明

• ./extracted: wikiextractorで抽出したファイルのディレクトリを指定
• -o wikipedia.db: 出力先のDB名を指定

この処理を実施すると、またWARNINGが出ます。wikiextractorが全ての改行コードを除去しきれていないのではと想像できる出方ですが、大した数ではなさそうなので今回は無視することにします。

オプションとして-qを指定するとWARNINGの表示を出さないようにできます。

データの確認

以上の手順で記事DBが作成されていると思います。

-rw-r--r--  1 dev-user users 3.2G 10月 11 18:05 wikipedia.db

3.2GBと結構大きいですが、次のコマンドなどで圧縮すると大分小さくなるのでコピーする際は参考にしてください。(圧縮には結構時間がかかります)

# tarコマンドで圧縮
tar zxvf wikipedia.db.tar.gz wikipedia.db

# 3.2GBが1.1GBまで圧縮されました。
-rw-r--r--  1 dev-user users 1.1G 10月 11 18:49 wikipedia.db.tar.gz

データの中身を見てみる

SQLite3をインストールする(入っていない場合)

もしインストールされていない場合は適宜インストールしてください。
下記はUbuntuの場合です。

sudo apt install sqlite3

記事数を確認してみる

sqlite3 wikipedia.db "select count(*) from pages;"
>> 1170762

約110万件のきじが登録されました。

記事テーブルの構造を確認してみる

sqlite3 wikipedia.db ".schema pages"
>> CREATE TABLE `pages` (
>>     `id` bigint PRIMARY KEY,
>>     `url` varchar(512),
>>     `title` varchar(1024),
>>     `text` text
>> );

wikiextractorが出力したJsonに合わせて筆者が軽く作成したTableなのでちょっと素朴すぎる感じがしますが、このようなTableにデータが格納されています。さすがにidがbigintになっているのは恥ずかしいので将来の更新で直すかもしれません。

タイトルに「アニメ」が含まれる記事の一覧を表示する

# まずは件数を確認
sqlite3 wikipedia.db 'select count(title) from pages where title like "%アニメ%";'
>> 1073

# 多いので、絞って表示
sqlite3 wikipedia.db 'select title from pages where title like "%アニメ%" limit 30;'
>> アニメ
>> SFアニメ
>> アニメ関係者一覧
>> アニメ作品一覧
>> アダルトアニメ
>> 東京を舞台とした漫画・アニメ作品一覧
>> アニメーション
>> 魔法少女アニメ
>> 宇宙の戦士 (アニメ)
>> ロボットアニメ
...

タイトルに「アニメ」が含まれる記事の本文データを取得する

やっとここまで来ました。次のコマンドで、anime.txtとして1073件分の記事データを取得できます。

sqlite3 wikipedia.db 'select text from pages where title like "%アニメ%";' > anime.txt

中身の確認

# 容量の確認。4.8MBです。
ls -lh anime.txt
>> -rw-r--r-- 1 dev-user users 4.8M 10月 11 20:03 anime.txt

# 行数
wc -l anime.txt
>> 37570 anime.txt

# 先頭の10行を表示
head -n 10 anime.txt
----------------------------------------
アニメ

アニメは、アニメーションの略語である。アニメーションを用いて構成された映像作品全般を指す。

各種メディアで提供されるサブカルチャーの一つ。

「文化芸術基本法」ではメディア芸術、関連法の「コンテンツの創造、保護及び活用の促進に関する法律」によるとコンテンツの一つと定義されており、い
ずれにおいてもアニメと略されてはおらず、アニメーションと正式表記されている。別定義として、多角的芸術分類観点において、美術（映像を含まない）
、映像、音楽、文学、芸能の総合芸術とされるときもある。

単に「アニメ」という場合は、セルアニメーション（セルアニメ）のことを指していることが多い。本項では、主に日本で製作された商業用セルアニメーシ
ョンについて解説する。
----------------------------------------

最後に

目的に応じて記事を絞り込むことができるので、特定の話題を対象とした文章の分析や試行錯誤中に使用する軽量の学習データとしてなど、アイディアによっては便利に使えるのではないでしょうか。

いつになるか分かりませんが、次回はカテゴリデータを追加して、カテゴリによる記事の抽出など紹介できればと思います。

正規表現って普段書かなくなると、すぐ忘れちゃいますよね？
めちゃめちゃTipsですが、備忘録も兼ねて投稿します。

やりたいこと

• 掲題の通りですが、金額等の文字列の修正を想定しています。
• 例： 金額欄に1,000.000と入力されたときに、1,000,000に修正したい。
  • ただし、1,000.000.00と入力されたときは1,000,000.00にしたい。
• 「そんなケースに遭遇する？」と思われるかもしれませんが、例えばOCRシステム等ではあり得るかと思います。（私の場合、まさにこれ）

コード

import re

target_text = '1,000.000.00'
replace_rule = re.compile('\.(?!\d{2}$)')

text = replace_rule.sub(',', target_text)

これでいけました。
なお、このままだと例えば 1,000.0は1,000,0となってしまうので、ちゃんと3つ区切りのドットだけを置換したい場合は（そもそもそんなケースが有るかは不明ですが）

import re

target_text = '1.000.00.00'
replace_rule = re.compile('\.(\d{3})')

text = replace_rule.sub(r',\1', target_text)

でいけました。

以上です。この記事が世界の誰かの役に立てば幸いです

SAMPLE

これが

こうなります。

PYTHON

1. 画像に減色処理を施した上で任意の画素数で平滑化します。
2. 画像からエッジを検出し反転させたものを線画として上記に透過します。
3. 2に減色処理後の画像を再度重みづけながら合成します。

初期設定を変更することで水彩画具合を調整できます。

import cv2
import numpy as np
import os.path
import matplotlib.pyplot as plt
from PIL import Image

#####初期設定
##読み込み画像
f = 'filename.jpg'
##平滑化を行う画素のサイズ
average_square = (5, 5)
##x軸方向の標準偏差
sigma_x = 0
##減色処理の際の配列サイズ
reshape_size = (-1, 3)
##減色処理の際の停止条件
criteria = (cv2.TERM_CRITERIA_EPS+cv2.TERM_CRITERIA_MAX_ITER, 10, 1.0)
##減色処理数
k=100
##線画抽出の際の平滑化サイズ
line_average = (9, 9)
##線画抽出の際のx軸方向の標準偏差
line_sigma_x = 0
##エッジ検出の際のthreshold1
threshold1 = 50
##エッジ検出の際のthreshold2
threshold2 = 55
##エッジを反転させる際の平滑化サイズ
edges_average = (1, 1)
##エッジを反転させる際のx軸方向の標準偏差
edges_sigma_x = 0
##二値化の際のthresh
thresh = 90
##二値化の際のmax_pixel
max_pixel = 255
##エッジを明るくする際のガンマ
gamma = 5.0
##合成画像の重み
multi_w = 0.5
##彩色表現の重み
paint_w = 0.9
##ガンマ
gamma = 1.5
#####初期設定終わり

##画像の読み込み
image = cv2.imread(f, 1)
##ファイル名と拡張子を取得
file, ext = os.path.splitext(f)

#####彩色表現
def paint(filename):
    ##画像をぼかす
    image_blurring = cv2.GaussianBlur(filename, average_square, sigma_x)
    ##減色処理
    z = image_blurring.reshape(reshape_size)
    z= np.float32(z)
    ret,label,center=cv2.kmeans(z, k, None, criteria, 10, cv2.KMEANS_RANDOM_CENTERS)
    center = np.uint8(center)
    res = center[label.flatten()]
    image_reshape = res.reshape((image_blurring.shape))
    ##再度ぼかす
    return cv2.GaussianBlur(image_reshape, average_square, sigma_x)

#####線画抽出
def line(filename):
    ##エッジを検出
    image_preprocessed  = cv2.cvtColor(cv2.GaussianBlur(filename, line_average, line_sigma_x), cv2.COLOR_BGR2GRAY)
    image_edges = cv2.Canny(image_preprocessed, threshold1 = threshold1, threshold2 = threshold2) 
    ##エッジを反転させる    
    image_h = cv2.GaussianBlur(image_edges, edges_average, edges_sigma_x)
    _, image_binary = cv2.threshold(image_h, thresh, max_pixel, cv2.THRESH_BINARY)
    image_binary = cv2.bitwise_not(image_binary)
    return image_binary

paint = paint(image)
cv2.imwrite(file + "_paint" + ext, paint)
line = line(image)
cv2.imwrite(file + "_line" + ext, line)

##乗算
def mul(input_color,mul_color):
    return int(round(((input_color*mul_color)/255),0))

##合成
def multiple(image1, image2):
    ##Pillowに変換
    pixelSizeTuple = image1.size
    image3 = Image.new('RGB', image1.size)
    for i in range(pixelSizeTuple[0]):
        for j in range(pixelSizeTuple[1]):
            r,g,b = image1.getpixel((i,j))
            r2,g2,b2 = image2.getpixel((i,j))
            img_r = mul(r,r2)
            img_g = mul(g,g2)
            img_b = mul(b,b2)
            image3.putpixel((i,j),(img_r,img_g,img_b))  
    return image3

image1 = Image.open(file + "_paint" + ext).convert("RGB")    
image2 = Image.open(file + "_line" + ext).convert("RGB")    

multi = multiple(image1, image2)
multi.save(file+"_multi"+ext)

##重みあり合成
def dodge(multi, paint):
    d = cv2.addWeighted(multi, multi_w, paint, paint_w, gamma)        
    return d

multi_image = cv2.imread(file+"_multi"+ext,1)
output = dodge(multi_image, paint)

cv2.imwrite(file + "_output" + ext, output)

EXAMPLE

加工する前の画像です。

減色処理を施し平滑化しました。

画像からエッジを検出し反転しました。

減色処理した画像にエッジを反転させた画像を投下しました。

再度減色処理した画像を重みを調整しながら合成しました。完成です。

女の子の画像です。

FBX SDK Python 2019 が対応しているファイル形式一覧

公式サイトを見ると、C++版のFBX SDKは大量のファイルフォーマットをサポートしているようなのですが、どうやらPython版はすべて対応しているわけではないみたいでしたので、対応しているフォーマットを調査してみました。

書き込み可能なファイル形式一覧

以下が一覧です。Index と Descriptionは、あくまで私の環境（後述）での値です。他の環境では異なる可能性がありますのでご注意ください。

Index	Description
0	FBX binary (*.fbx)
1	FBX ascii (*.fbx)
2	FBX encrypted (*.fbx)
3	FBX 6.0 binary (*.fbx)
4	FBX 6.0 ascii (*.fbx)
5	FBX 6.0 encrypted (*.fbx)
6	AutoCAD DXF (*.dxf)
7	Alias OBJ (*.obj)
8	Collada DAE (*.dae)
9	Biovision BVH (*.bvh)
10	Motion Analysis HTR (*.htr)
11	Motion Analysis TRC (*.trc)
12	Acclaim ASF (*.asf)
13	Acclaim AMC (*.amc)
14	Vicon C3D (*.c3d)
15	Adaptive Optics AOA (*.aoa)
16	Superfluo MCD (*.mcd)

読み込み可能なファイル形式一覧

以下が一覧です。Index と Descriptionは、あくまで私の環境（後述）での値です。他の環境では異なる可能性がありますのでご注意ください。

Index	Description
0	FBX (*.fbx)
1	AutoCAD DXF (*.dxf)
2	Alias OBJ (*.obj)
3	3D Studio 3DS (*.3ds)
4	Collada DAE (*.dae)
5	Alembic ABC (*.abc)
6	Biovision BVH (*.bvh)
7	Motion Analysis HTR (*.htr)
8	Motion Analysis TRC (*.trc)
9	Acclaim ASF (*.asf)
10	Acclaim AMC (*.amc)
11	Vicon C3D (*.c3d)
12	Adaptive Optics AOA (*.aoa)
13	Superfluo MCD (*.mcd)
14	(*.zip)

調査方法

Index と Description は以下のコードで確認しました。

# 書き込み可能な形式一覧を出力
for formatIndex in range(manager.GetIOPluginRegistry().GetWriterFormatCount()):
    description = manager.GetIOPluginRegistry().GetWriterFormatDescription(formatIndex)
    print(formatIndex, description)

# 読み込み可能な形式一覧を出力
for formatIndex in range(manager.GetIOPluginRegistry().GetReaderFormatCount()):
    description = manager.GetIOPluginRegistry().GetReaderFormatDescription(formatIndex)
    print(formatIndex, description)

以下の２つの環境で確認しました。

Windows環境

• Windows 10 64bit
• Python 2.7
• FBX SDK Python 2019.5

Linux環境

• Alpine Linux v3.9
• Python 2.7
• FBX SDK Python 2019.2

さいごに

読み込みと書き込みでIndexが変わるのがちょっとややこしいですね。。。

以下のサイトを参考にさせていただきました。ありがとうございます。

• Autodesk FBX 2017 HELP / FbxIOPluginRegistry Class Reference

FBX SDK Python 2019 が対応しているファイル形式一覧

公式サイトを見ると、C++版のFBX SDKは大量のファイルフォーマットをサポートしているようなのですが、どうやらPython版はすべて対応しているわけではないみたいでしたので、対応しているフォーマットを調査してみました。

書き込み可能なファイル形式一覧

以下が一覧です。Index と Descriptionは、あくまで私の環境（後述）での値です。他の環境では異なる可能性がありますのでご注意ください。

Index	Description
0	FBX binary (*.fbx)
1	FBX ascii (*.fbx)
2	FBX encrypted (*.fbx)
3	FBX 6.0 binary (*.fbx)
4	FBX 6.0 ascii (*.fbx)
5	FBX 6.0 encrypted (*.fbx)
6	AutoCAD DXF (*.dxf)
7	Alias OBJ (*.obj)
8	Collada DAE (*.dae)
9	Biovision BVH (*.bvh)
10	Motion Analysis HTR (*.htr)
11	Motion Analysis TRC (*.trc)
12	Acclaim ASF (*.asf)
13	Acclaim AMC (*.amc)
14	Vicon C3D (*.c3d)
15	Adaptive Optics AOA (*.aoa)
16	Superfluo MCD (*.mcd)

読み込み可能なファイル形式一覧

以下が一覧です。Index と Descriptionは、あくまで私の環境（後述）での値です。他の環境では異なる可能性がありますのでご注意ください。

Index	Description
0	FBX (*.fbx)
1	AutoCAD DXF (*.dxf)
2	Alias OBJ (*.obj)
3	3D Studio 3DS (*.3ds)
4	Collada DAE (*.dae)
5	Alembic ABC (*.abc)
6	Biovision BVH (*.bvh)
7	Motion Analysis HTR (*.htr)
8	Motion Analysis TRC (*.trc)
9	Acclaim ASF (*.asf)
10	Acclaim AMC (*.amc)
11	Vicon C3D (*.c3d)
12	Adaptive Optics AOA (*.aoa)
13	Superfluo MCD (*.mcd)
14	(*.zip)

調査方法

Index と Description は以下のコードで確認しました。

# 書き込み可能な形式一覧を出力
for formatIndex in range(manager.GetIOPluginRegistry().GetWriterFormatCount()):
    description = manager.GetIOPluginRegistry().GetWriterFormatDescription(formatIndex)
    print(formatIndex, description)

# 読み込み可能な形式一覧を出力
for formatIndex in range(manager.GetIOPluginRegistry().GetReaderFormatCount()):
    description = manager.GetIOPluginRegistry().GetReaderFormatDescription(formatIndex)
    print(formatIndex, description)

manager には FbxManager.Create() で生成したインスタンスを入力してください。

サンプルプロジェクトはこちらに置きます。
GitHub/segurvita/fbx_sdk_python_sample/convert_fbx/list_format.py

以下の２つの環境で試し、同じ値が出力さえることを確認しました。

Windows環境

• Windows 10 64bit
• Python 2.7
• FBX SDK Python 2019.5

Linux環境

• Alpine Linux v3.9
• Python 2.7
• FBX SDK Python 2019.2

さいごに

読み込みと書き込みでIndexが変わるのがちょっとややこしいですね。。。

以下のサイトを参考にさせていただきました。ありがとうございます。

• Autodesk FBX 2017 HELP / FbxIOPluginRegistry Class Reference

問題

大型のWindows Update後，コマンドラインでpythonを実行すると，すでにpython環境が構築されている場合でも対話型インタプリタではなく，MS Storeのpythonのページが開くようになった.
詳しい記事

解決方法

アプリ実行エイリアスの欄のpythonをはずしてやればよい．

S3に保存されたgzip/zip圧縮ファイルを、ローカルファイルとして保存することなく読み込みたい。

前提条件

• Python: 3.7.4
• AWS SDK: boto3-1.9.230以上

ポイント

• 標準のgzip, zipfileパッケージを使ってファイルを読み込んでファイルオブジェクトに変換する。
  

（bzip2とかは考えなくて良いんじゃないかな）

つまり、以下のようにして読み込んだ際と同様に扱いたい。

import gzip
gz_file = 'path/to/file.csv.gz'
file = gzip.open(gz_file, 'rt')

file.read()
..

パッケージ毎にファイルオブジェクトの変換方法に癖があるので、それぞれ対応する。
ioパッケージのBytesIOやらTextIOWrapperを使う。
ついでに、圧縮形式も正しいかチェックする。

コード

gzip圧縮の場合

import boto3
import gzip
import io


# 圧縮形式が正しいかチェック
def valid_gzip_format(obj):
    file = gzip.open(io.BytesIO(obj), 'rt')
    try:
        file.read()
        return True
    except OSError:
        return False


s3 = boto3.client('s3')
bucket = 'bucket_name'
key = 'path/to/file.csv.gz'
# S3オブジェクトの取得
obj = s3.get_object(
    Bucket=bucket,
    Key=key
)['Body'].read()

# main
if valid_gzip_format(obj):
    file = gzip.open(io.BytesIO(obj), 'rt')
else:
    raise InvalidCompressError('ファイル拡張子(.gz)に対して圧縮形式が正しくありません。')

# 中身確認
for row in file.readlines():
    print(row.replace('\n', ''))

zip圧縮の場合

import boto3
import io
import zipfile


# 圧縮形式が正しいかチェック
def valid_zip_format(obj):
    return zipfile.is_zipfile(io.BytesIO(obj))


s3 = boto3.client('s3')
bucket = 'bucket_name'
key = 'path/to/file.csv.zip'
# S3オブジェクトの取得
obj = s3.get_object(
    Bucket=bucket,
    Key=key
)['Body'].read()

# main
if valid_zip_format(obj):
    zf = zipfile.ZipFile(io.BytesIO(obj))
    if len(zf.namelist()) > 1:
        raise InvalidCompressError('zipアーカイブの中にファイルが複数あります。')
    file = io.TextIOWrapper(zf.open(zf.namelist()[0], 'r'))
else:
    raise InvalidCompressError('ファイル拡張子(.zip)に対して圧縮形式が正しくありません。')

# 中身確認
for row in file.readlines():
    print(row.replace('\n', ''))

zipfileではopenのモードに'rt'が無いのでTextIOWrapperを使う必要がある。
参考 - ZIP アーカイブの処理 — Python 3.7.5rc1 ドキュメント

バージョン 3.6 で変更: Removed support of mode='U'. Use io.TextIOWrapper for reading compressed text files in universal newlines mode.

まとめ

上記のコードで取得したS3オブジェクトはbytes型で保持されるだけなので、with構文で扱う必要はないと思う（合ってる？）。
あと実際に中身を読み込む際に、末尾に改行コードが含まれるのでそのあたりも考慮が必要（もろちん、バイナリモードでも含まれるけど）。

はじめに

　モデルの汎化性能を評価する上で、より頑健な手法である交差検証を導入する。また、クラス分類性能と回帰性能を評価する上で、scoreでのデフォルトのR^2よりも良い手法について議論されてます。
　教師あり学習で汎化性能を最大にするように効率的にパラメータを調整するグリッドサーチにも触れます。

交差検証

　交差検証では、データの分割を何度も繰り返して行うことで、これまでのただの分割よりも安定した評価を実現する。最もよく用いられるのは、k分割交差検証である。例えば、k=5の場合は、分割0から分割4の5つに分割され、分割0から順番にテストセットに、残りの分割を訓練セットにして学習を行い、性能の平均をとる。

（層化）k分割交差検証

from sklearn.model_selection import cross_val_score

　データをk個に分割する際に、先頭から1/k個取る方法はうまくいくとは限らない。なぜなら、データがクラス0から順番に保存されているデータセットもあるためだ。そのため、scikit-learnでは各分割内でのクラスの比率が全体の比率と同じになるように、層化交差検証を行っている。

cross_val_scoreでは、デフォルトで3分割交差検証を行い、3つの精度を返す。
cv　分割数

1つ抜き交差検証

from sklearn.model_selection import LeaveOneOut

　k分割交差検証の個々の分割が1サンプルしかないと考え、毎回、テストセットの中の1サンプルだけをテストセットとして検証する。大規模なデータセットでは非常に時間がかかるが、小さいデータセットではよりよい推定ができる。

シャッフル分割交差検証

from sklearn.model_selection import ShuffleSplit

　データセットの中から指定した数・割合の訓練セットとテストセットを選び、指定した回数の検証を行う。サブサンプリングと呼ばれる、データの一部だけを用いるようにすることもでき、データセットが大きい場合に有用である。
train_size　訓練セットの数・割合
test_size　テストセットの数・割合
n_iter　繰り返しの回数

グループ付き交差検証

from sklearn.model_selection import GroupKFold

　汎化性能を高めるためには、データとしては異なるものの、訓練セットとテストセットに分かれることが好ましくないことがある。例えば、表情から感情を予測する際に、同じ人のデータが分かれることである。なぜならば、学習済の人の表情の予測は容易になるはずであり、未知の人の表情から感情を予測する汎化性能を計るには適切でないためである。
　そのため、配列を用意し、同じグループに属するデータには同じ番号を振ることで、それらのデータをグループ化し、個々のグループはまとまって訓練セットもしくはテストセットに入るように出来る。
groups　グループ分けを行うための配列

グリッドサーチ

　グリッドサーチとは、基本的にはパラメータのすべての組み合わせに対して試して見る方法である。ただし、それぞれのパラメータに対する評価に対してテストセットを使っていては、汎化性能を適切に評価できているとは言い難い。そのため、訓練セット・検証セット・テストセットとデータセットを3つに分割する。

交差検証を用いたグリッドサーチ

from sklearn.model_selection import GridSearchCV

　グリッドサーチは一般に用いられているが、データの実際の分割のされ方によって大きく性能が変わる。そのため、それぞれのパラメータの組み合わせに対して交差検証を行うとよい。ただし、パラメータ2つをそれぞれ6つの値で、交差検証を5回行うとすると、6×6×5=180通りのモデルを訓練する必要がある。これには非常に時間がかかることが大きな問題点である。

ネストした交差検証

　上記のグリッドサーチでは、データセットを訓練セットとテストセットに一度だけ分け、訓練セット上で交差検証を行い、訓練セットと検証セットに分けた。訓練セットとテストセットの分割が一度だけの場合、その分割の仕方に結果が依存する場合がある。そのため、訓練セットとテストにも交差検証による分割を適用する。この手法をネストした交差検証と呼ぶ。
　ネストした交差検証をscikit-learnで行うのは容易で、cross_val_scoreをGridSearchCVのインスタンスをモデルとして呼び出せば良い。

評価基準とスコア

　機械学習の基準を選ぶときには、そのアプリケーションによって最終的な目的を考える必要がある。横軸を予測ラベル、縦軸を新ラベルとし、あてはまるデータの数を成分とした行列を混同行列と呼ぶ。

2クラス分類

　2クラス分類では、2つのクラスを陽性と陰性と呼ぶ。予測・実際の値がともに陽性のときを真陽性（TP）、予測・実際の値がともに陰性のときを真陰性（TN)と呼ぶ。また、陽性と予測したが実際は陰性のときを偽陽性（FP）、この逆のときを偽陰性（FN）と呼ぶ。

from sklearn.metrics import classification_report

　以下の3つの評価基準をクラスごとにそのクラスを陽性とした場合に値を出力してくれる。同時に支持度と呼ばれる、陽性のデータの数も出力してくれる。

適合率

適合率=\frac{TP}{TP+FP}

　適合率は、陽性であると予測されたものがどれだけ実際に陽性であったかを表す。偽陽性の数を制限したい場合に基準として用いられる。

再現率

再現率=\frac{TP}{TP+FN}

　再現率は、実際に陽性のデータのうち、陽性と予測されたものの割合を示す。偽陰性を避けることが重要な場合に基準として用いられる。例えば、ガンの診断。

f-値

F=2×\frac{適合率×再現率}{適合率+再現率}

　適合率と再現率がトレードオフである。すべてを陽性と判定すれば再現率は100%となるが、適合率は非常に小さくなる。一方、最も真である可能性が高いもの1つだけを真とすれば、適合率は100%となるが、再現率は非常に小さくなる。
　適合率と再現率はそれぞれ非常に重要な基準であるが、一方だけでは全体像がつかめない。そのため、これら2つを包括的に評価するf-値という基準がある。これは、適合率と再現率の調和平均によって求められる。（f-値の変種の1つであり、他にもあるらしい）

適合率・再現率カーブ

from sklearn.metrics import precision_recall_curve, average_precision_score

　カーブが右上に近いほど、適合率と再現率が同時に高いため、良いクラス分類器であることがわかる。カーブを要約する方法の1つとして、カーブ下の領域（AUC）を積分することがあり、これは平均適合率とも呼ばれる。

受信者動作特性カーブ（ROCカーブ）

from sklearn.metrics import roc_curve, roc_auc_score

　横軸に偽陽性率、縦軸に真陽性率（再現率）をとる。低い偽陽性率を保ちながら、高い再現率を達成したいので、左上に近いほどよいクラス分類器と言える。同様にして、ROCカーブのAUCを積分することで要約することができる。

予測の不確実性

　ほとんどのクラス分類器には、予測の不確実性を評価するためのdecision_functionメソッドもしくはpredict_probaメソッドがある。予測の閾値は、デフォルトでは前者が0、後者が0.5となっている。例えば、偽陽性が増えてもよいからより多くの真陽性を得たい場合は、閾値を小さくすることで再現率を上げることができる。

他クラス分類

　他クラス分類においても、classification_reportによって適合率・再現率・f-値を求めることができる。クラスごとのf-値を平均する方法は以下の3つがある。

macro　重みをつけずにクラスごとのf-値を平均する。個々のクラスを同じように重視している。
weighted　各クラスの支持度に応じて重みをつけてクラスごとのf-値を平均する。
micro　すべてのクラスの偽陽性・偽陰性・真陽性の総数を計算し、その値を用いて計算する。個々のサンプルを同じように重視している。

回帰

　R^2スコアが最も直観的。

参考文献

Pythonではじめる機械学習
――scikit-learnで学ぶ特徴量エンジニアリングと機械学習の基礎
Andreas C. Muller、Sarah Guido　著、中田 秀基　訳

はじめに

マイクで収録した音にFIRフィルタで周波数補正をかけたくなりました。FIRフィルタというとscipy.signal.firwinなどがありますが、今回は単純なローパスフィルタやハイパスフィルタではなく、複雑な周波数特性の実現を目指します。私はFIRフィルタの理論を深く理解しているわけではないので、間違っている箇所があると思います。そのときはコメントで教えていただけるとありがたいです。

参考

東北大学講義資料の以下を参考にさせていただきました。
やる夫で学ぶディジタル信号処理

手順

ファイル読み込みなどのコードは省略しました。

1. マイク付属の周波数特性表の読み込み

これはとあるマイクの周波数特性と自由音場補正データです。

補正すべき値は以下のようになります。

2. 実現したい周波数特性をつくる

ここで、補正値は対数から線形になおしておきます。
このマイクは140 kHzまでが計測可能範囲なので、それより高い周波数ではゲインを0にします。

N = 128 #タップ数 + 1
fs = 500.0e3 #補正をかけたいデータのサンプリングレート
freq_array = np.linspace(0.0, fs, num=N)

gain_correction = np.zeros(N//2)
lim_freq = freq_array[freq_array <= 140.0e3]
# f_pro_funcとf_pro2_funcはスプライン補間テーブル
gain_correction[0:len(lim_freq)] = 10**(-(f_pro_func(lim_freq) + f_pro2_func(lim_freq))/20.0)
gain_correction = list(gain_correction) + list(gain_correction[::-1]) #反転して合わせる

fig, axes = plt.subplots(figsize=(8,4))
axes.plot(freq_array/1000.0, gain_correction)
axes.set_xlabel("Frequency kHz")
axes.set_ylabel("Relative response")
axes.set_ylim(-0.1, 1.1)

3. FIRフィルタの係数を計算する

波形データにフィルタ係数を畳み込み計算した結果は、それぞれのデータのフーリエ変換の積になります。フィルタ係数のフーリエ変換は、実現したい周波数特性そのものです。従って、先程の周波数特性を逆フーリエ変換すればフィルタの係数が求まります。また、FIRフィルタの係数は左右対称になります。

b = list(np.real(np.fft.ifft(gain_correction)))[0:N//2]
b = b[::-1] + b[1:]
b = b * np.hamming(N-1)
b /= np.sum(b) # 係数の和を1にする

フィルタ係数をフーリエ変換してみます。

4. FIRフィルタの適用

計算したFIRフィルタをホワイトノイズに適用して、適用された波形のスペクトルを確認してみます。

NN = 8192
ave_num = 1000
y = np.random.rand(NN*ave_num) # ホワイトノイズ
yy = signal.lfilter(b, 1, y) # FIRフィルタを適用した信号

fft_sum_w = np.zeros(NN)
fft_sum_f = np.zeros(NN)
for i in range(ave_num):
    fft_sum_w += np.abs(np.fft.fft(y[NN*i:NN*(i+1)]) / NN)
    fft_sum_f += np.abs(np.fft.fft(yy[NN*i:NN*(i+1)]) / NN)
fft_ave_w = 20.0 * np.log10(fft_sum_w / ave_num)
fft_ave_f = 20.0 * np.log10(fft_sum_f / ave_num)

fig, axes = plt.subplots(figsize=(8,4))
axes.plot(np.linspace(0.0, 500.0, num=NN), fft_ave_w)
axes.plot(np.linspace(0.0, 500.0, num=NN), fft_ave_f)
axes.set_xlim(1.0, 250.0)
axes.set_ylim(-50.0-15, -50.0)
axes.set_xscale("log")
axes.set_xlabel("Frequency kHz")
axes.set_ylabel("Amplitude dB")

ちゃんと実現したかった周波数特性になっています。

まとめ

実現したい周波数特性からFIRフィルタの係数を計算し、波形に適用しました。
今後はリアルタイム信号処理にも挑戦してみたいです。

徒競走のタイムと違って、対戦型ゲームの巧さというのは簡単には数値化できません。 しかし、なんとかしてそれを実現させようとして生まれたのがレーティングです。 古典的にはチェスのために生み出された Elo レーティングが存在し、様々な改良アルゴリズムが生み出されてきました。

TrueSkill は Microsoft が開発したレーティングアルゴリズムです。 Microsoft が開発したとあって Xbox Live のゲームで使用されているようです。 このアルゴリズムには既存のレーティングアルゴリズムと比較して以下のような特徴があるそうです。

• 収束が早い。 レーティングに初めて参加するプレイヤーの実力を推定するのに何度も何度も不適当なマッチングで対戦する必要がない。
• 複数人による対戦に対応している。 勝ちか負けかのみならず順位を定めるようなゲームやチーム戦¹のゲームにも使用できる。
• ゲームへの参加に重み付けができる。 例えばチームメンバーのひとりが回線トラブルにより途中でゲームから抜けた場合でも適用できる。
• Microsoft の息がかかっている。 特許申請と商標登録がなされており、なんとなく手が出しづらい。 あと一部の狂信的 OSS 主義者が憤死する。

ちょろっと論文を読んでみようとは試みたのですが、まず私は英語がてんでダメで、おまけに数学統計学の知識も付け焼き刃なもので、因子グラフとか言われてもさっぱりぽんで無事死亡しました。 私よりはマシという自信がある方は Qiita に投稿されているこちらの記事をお読みいただくと理解を深められると思います。 ぶっちゃけこれと“Computing Your Skill”の和訳記事以外に日本語の資料ないです。 あと Wikipedia もか……?

しかし、私たちには他人の褌があります。 そう、ここに TrueSkill の Python 実装があります。 いくら英語がダメな私でも、Google 先生に頼ればリファレンスくらいは読めるでしょう。 ということで一念発起して TrueSkill で遊んでみて簡単な使い方と遊んだ結果を残すことで数少ない TrueSkill の日本語資料の一つとして労せずうまい汁をすすりたいと思います。 誰かがやめ太郎さん風に書き直して「ワイ『TrueSkill? 転職支援サイトかなんかやろか?』」みたいな感じで投稿するとさらにバズるかもしれません。 次のうまい汁は君だ!

簡単な統計学と TrueSkill の概要

TrueSkill ではプレイヤーのスキルを正規分布に従う確率密度関数によって表現します。 このように序盤から難しそうな単語を出していくことで読者を威圧していくテクニックが荒んだ現代社会で生きのこるには必須と言われています。

しかし皆さんは私のこの姑息な手法に怯える必要はありません。 まず下の図を見てください。

正規分布とは、まあおおむねこんな形をしたものです。 確率密度関数とは……たとえばそうですね、横軸を「ある値」だとします。 今回の場合はそう、スキルです。 縦軸は……ちょっと違うのですが確率だと思ってください。 正確な話をすると、ちょっとだけ数学的な話になりますが、確率密度関数 $\phi_{\mu, \sigma^2}(x)$ とは²、「ある値」が $x_1$ から $x_2$ の間(ここで $x_1 < x_2$ とします)に収まる確率 $p_{12}$ が $\phi_{\mu, \sigma^2}(x)$ を $x_1$ から $x_2$ まで積分した値、すなわち $p_{12} = \int_{x_1}^{x_2}\phi_{\mu, \sigma^2}(x)$ になるように定めた関数です。 まあ要するに山の頂点周辺の可能性が高いよ関数、と思ってくれればいいです。

ここで私たちが気にしなければいけないことはふたつです。 すなわち、

• 頂点はどこか?
• どれくらい幅があるのか?

頂点が右にあるほうがよりよいスキルを持っている可能性が高いですし(青に対する橙)、横軸に対する頂点が同じでも幅が狭いほど精密な予測になっているといえます(青に対する緑)。 ちなみに定義上、幅が狭くなれば頂点は高くなります。

前者を教えてくれるのが期待値あるいは平均と呼ばれる数値で、ギリシア文字の $\mu$ やラテン文字の $m$ で表されます。 正規分布における期待値はそのまま山の頂点の横軸成分になります。 マイナス無限からプラス無限まで重み付き積分をしたときの平均、ということなんだと……思います……。(統計クソザコマン)

後者を教えてくれるのが標準偏差と呼ばれる数値で、ギリシア文字の $\sigma$ やラテン文字の $s$ で表されます。 ちなみに、「ある値」が $\mu \pm \sigma$ の範囲に収まっている確率は、定義から計算すると約 68% になります。 これが $\mu \pm 3\sigma$ になるとほぼ確実³となり、この $3\sigma$ という指標は実用上よく利用されます。 外れ値の検出とかね。

TrueSkill ではこの 2 種類の値によってプレイヤーのスキルを推測します。

プレイヤーが勝てば $\mu$ は上昇し、負ければ $\mu$ は減少します。 いずれにせよ $\sigma$ は減少し⁴、推測がより正確になったことを示します。 プレイヤーに提示される数値は $\mu - 3\sigma$ という、期待値をより「控えめ」に見積もった値で、これは Python 実装のリファレンスでは“rating exposure”と表現されています。 どう訳すのが適切かはよくわかりませんが、とりあえず本稿では「顕在化レート / 顕在化レーティング」としておきます。

はじめて TrueSkill のレーティングシステムに参加したプレイヤーの値は $\mu = 25, \sigma = \frac{25}{3}$ で与えられます。 顕在化レーティングはちょうど 0 となります。

利用してみる

堅苦しい理論はここまでにして、実際に TrueSkill を利用してみましょう。

当然の権利のように PyPI に登録されていますので、pipを使用して簡単にインストールできます。

$ pip install trueskill

インポートするにはこうです。

In[*]

import trueskill

環境

TrueSkill によるレーティング計算を行うには、環境を用意します。 環境といってもなにかをインストールするとかそういう意味ではなく、いわば「利用する定数を設定する」くらいの意味合いで考えて下さい。

新しい環境を生成するにはtrueskill.TrueSkillコンストラクタを呼び出します。

In[*]

mu = 25.
sigma = mu / 3.
beta = sigma / 2.
tau = sigma / 100.
draw_probability = 0.1
backend = None

env = trueskill.TrueSkill(
    mu=mu, sigma=sigma, beta=beta, tau=tau,
    draw_probability=draw_probability, backend=backend)

env

Out[*]

trueskill.TrueSkill(mu=25.000, sigma=8.333, beta=4.167, tau=0.083, draw_probability=10.0%)

パラメータmuとsigmaは、より実態に即した表記をすると $\mu_0$ および $\sigma_0$ で、$\mu$  および $\sigma$ の初期値を示します。 デフォルト値は25.および8.333333333333334です。

パラメータbetaは 76% の勝率を保証するスキル差です。 チェスや将棋や囲碁のような、運の要素が絡まないゲームにおいてはこの値は小さく、麻雀のような多少の実力差を運でカバーできるゲームではこの値は大きくするのが理想的です。 デフォルト値は4.166666666666667です。

パラメータtauはスキルの再計算を行うときに問答無用で $\sigma$ に与えられます。 これは不確かさを示し、ゲームのダイナミクスを表現するとか、$\sigma$ が 0 になって更新が停滞するのを防ぐとかの意味があるようです。 デフォルト値は0.08333333333333334です。

パラメータdraw_probabilityは、名前の通り 2 チーム(あるいはふたり)が引き分けになる確率です。 float値で入力することもできますが、floatを返却する関数を与えることもできるそうです(未確認)。 デフォルト値は0.1です。

パラメータbackendは計算に使用するバックエンドを指定します。 利用可能なバックエンドは'scipy'と'mpmath'です。 Noneでは TrueSkill モジュール自前のものを使用します。 デフォルト値はNoneです。

以降は主にこの環境に対してメソッドを呼び出すことによって TrueSkill を利用していきます。 例えば、以下のような具合です。

In[*]

alice = env.create_rating()
bob = env.create_rating()
(alice,),(bob,), = env.rate(((alice,), (bob,),), ranks=[0, 1,])
print(f'Alice\'s rating exposure: {env.expose(alice):.3f}')
print(f'  Bob\'s rating exposure: {env.expose(bob):.3f}')

Out[*]

Alice's rating exposure: 7.881
  Bob's rating exposure: -0.910

ユーザが生成する環境のほかにも、あらかじめ設定されたグローバル環境を利用することもできます。 グローバル環境はtrueskill.global_env関数で呼び出すことができます。 グローバル環境の初期値は前述のデフォルト値になっています。

In[*]

trueskill.global_env()

Out[*]

trueskill.TrueSkill(mu=25.000, sigma=8.333, beta=4.167, tau=0.083, draw_probability=10.0%)

グローバル環境を取得してメソッドコールをすることもできますが、多くのメソッドには代替(proxy)関数が存在します。

メソッドコール	代替関数
trueskill.global_env().rate	trueskill.rate
trueskill.global_env().quality	trueskill.quality
trueskill.global_env().expose	trueskill.expose

生成した環境をグローバル環境に登録するにはmake_as_globalメソッドを使用します。

In[*]

env = trueskill.TrueSkill(mu=50.)
env.make_as_global()

trueskill.global_env()

Out[*]

trueskill.TrueSkill(mu=50.000, sigma=8.333, beta=4.167, tau=0.083, draw_probability=10.0%)

trueskill.setup関数を使用すると直接グローバル環境の設定を変更することができます。

In[*]

trueskill.setup(mu=10., sigma=10./3., beta=10./6., tau=10./300.)

trueskill.global_env()

Out[*]

trueskill.TrueSkill(mu=10.000, sigma=3.333, beta=1.667, tau=0.033, draw_probability=10.0%)

レーティングとチーム

環境に対しcreate_ratingメソッドを呼び出すことで各プレイヤーのレーティングを示すRatingオブジェクトを生成できます。 create_ratingにmuパラメータおよびsigmaパラメータを与えることで初期値以外のレーティングを生成することもできます。 レーティングはmuプロパティおよびsigmaプロパティをもちます。

In[*]

env = trueskill.TrueSkill()
r1 = env.create_rating()
r2 = env.create_rating(mu=30, sigma=4)

(r1.mu,r1.sigma,),(r2.mu,r2.sigma,),

Out[*]

((25.0, 8.333333333333334), (30.0, 4.0))

Ratingオブジェクトのコンストラクタを直接呼び出した場合、グローバル環境が利用されます。

In[*]

r1 = trueskill.Rating()
trueskill.setup(mu=10., sigma=10./3., beta=10./6., tau=10./300.)
r2 = trueskill.Rating()

r1,r2,

Out[*]

(trueskill.Rating(mu=25.000, sigma=8.333),
 trueskill.Rating(mu=10.000, sigma=3.333))

環境に対しexposeメソッドを呼び出し、パラメータにレーティングオブジェクトを与えることでそのレーティングオブジェクトの顕在化レーティングを得ることができます。 Wikipedia 情報によると 0 - 50 のスケールが使用されているそうですが、このメソッドの返却値はクリッピングなどは行わないので注意が必要です。

In[*]

env = trueskill.TrueSkill()
r1 = env.create_rating()
r2 = env.create_rating(mu=30, sigma=4)
r3 = env.create_rating(mu=10, sigma=5)

env.expose(r1),env.expose(r2),env.expose(r3),

Out[*]

(0.0, 18.0, -5.0)

ゲームに参加する 1 単位のことをチームといいます。 チームはレーティングオブジェクトのリストまたはタプルか、値にレーティングオブジェクトをもつ辞書です。 たとえ 1 対 1 のゲームであってもチームを形成します。

In[*]

env = trueskill.TrueSkill()
r1 = env.create_rating()
r2 = env.create_rating()
t1 = (r1,r2,)

t1

Out[*]

(trueskill.Rating(mu=25.000, sigma=8.333),
 trueskill.Rating(mu=25.000, sigma=8.333))

レーティングの更新

rateメソッドを呼び出すことでレーティングの更新を行います。 第一パラメータにはチームのリストかタプルを、ranksパラメータには順位(プログラムの世界なので 0 始まりです)のリストかタプルを与えます。 更新されたチームのリストが返却されます。 前述したとおりチームには複数の書式がありますが、rateメソッドに渡すチームの書式は統一されている必要があります。

In[*]

env = trueskill.TrueSkill()

r1 = env.create_rating()
r2 = env.create_rating()
r3 = env.create_rating()
r4 = env.create_rating()
t1 = {'player1': r1, 'player2': r2,}
t2 = {'player3': r3, 'player4': r4,}

t1,t2, = env.rate((t1,t2,), ranks=(1,0,))

t1,t2,

Out[*]

({'player1': trueskill.Rating(mu=21.892, sigma=7.774),
  'player2': trueskill.Rating(mu=21.892, sigma=7.774)},
 {'player3': trueskill.Rating(mu=28.108, sigma=7.774),
  'player4': trueskill.Rating(mu=28.108, sigma=7.774)})

ranksパラメータを与えなかった場合、記述した通りの順位だったとみなして計算を行います。

In[*]

env = trueskill.TrueSkill()

r1 = env.create_rating()
r2 = env.create_rating()
r3 = env.create_rating()
r4 = env.create_rating()
t1 = {'player1': r1, 'player2': r2,}
t2 = {'player3': r3, 'player4': r4,}

t1,t2, = env.rate((t1,t2,))

t1,t2,

Out[*]

({'player1': trueskill.Rating(mu=28.108, sigma=7.774),
  'player2': trueskill.Rating(mu=28.108, sigma=7.774)},
 {'player3': trueskill.Rating(mu=21.892, sigma=7.774),
  'player4': trueskill.Rating(mu=21.892, sigma=7.774)})

前述したとおり、たとえ 1 対 1 のゲームであってもチームを形成します。 rateメソッドにレーティングオブジェクトを直接渡すことはできないので注意してください。

In[*]

# 間違い

env = trueskill.TrueSkill()

r1 = env.create_rating()
r2 = env.create_rating()

r1,r2 = env.rate((r1,r2,))

r1,r2

Out[*]

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-11-08f690b7e5cd> in <module>
      6 r2 = env.create_rating()
      7 
----> 8 r1,r2 = env.rate((r1,r2,))
      9 
     10 r1,r2

f:\trueskill_test\venv\lib\site-packages\trueskill\__init__.py in rate(self, rating_groups, ranks, weights, min_delta)
    475 
    476         """
--> 477         rating_groups, keys = self.validate_rating_groups(rating_groups)
    478         weights = self.validate_weights(weights, rating_groups, keys)
    479         group_size = len(rating_groups)

f:\trueskill_test\venv\lib\site-packages\trueskill\__init__.py in validate_rating_groups(self, rating_groups)
    272             raise TypeError('All groups should be same type')
    273         elif group_types.pop() is Rating:
--> 274             raise TypeError('Rating cannot be a rating group')
    275         # normalize rating_groups
    276         if isinstance(rating_groups[0], dict):

TypeError: Rating cannot be a rating group

In[*]

# 正解

env = trueskill.TrueSkill()

r1 = env.create_rating()
r2 = env.create_rating()

((r1,),(r2,),) = env.rate(((r1,),(r2,),))

r1,r2,

Out[*]

(trueskill.Rating(mu=29.396, sigma=7.171),
 trueskill.Rating(mu=20.604, sigma=7.171))

rateメソッドのweightsパラメータに 2 次元のリストかタプルを与えることで、ゲームの参加率を加味したレーティングの更新を行います。 チームと異なり、辞書による書式は使用できません。 チームを辞書で記述した場合順序が保証されませんので後述する書式によって参加率を記述するのが基本になります。

In[*]

env = trueskill.TrueSkill()

r1 = env.create_rating()
r2 = env.create_rating()
r3 = env.create_rating()
r4 = env.create_rating()

t1 = (r1,r2,)
t2 = (r3,r4,)
teams = (t1,t2,)

w1 = (1,1,)
w2 = (0.5,1,) # プレイヤー3はゲームの半分が経過した時点で回線落ちした、など
weights = (w1,w2,)

t1,t2, = env.rate(teams, weights=weights)

t1,t2,

Out[*]

((trueskill.Rating(mu=26.738, sigma=7.844),
  trueskill.Rating(mu=26.738, sigma=7.844)),
 (trueskill.Rating(mu=24.131, sigma=8.214),
  trueskill.Rating(mu=23.262, sigma=7.844)))

そもそもチーム 2 が不利な状況で負けたためか先程までの例にくらべレーティングの上昇/下降幅が小さく、また参加時間が少なかったプレイヤー 3 の $\sigma$ がほかの 3 人にくらべてあまり減っていないのがわかります。

現実問題としては、ほぼ全員の参加率が 1 で一部のプレイヤーのみが少ない、という傾向になると予想できます。 こういった場合、リストおよびタプルによる参加率の書式は、参加人数が増えたときに煩雑な記述になる問題があります。 そこで、参加率が 1 に満たないプレイヤーだけを辞書によって記述する方法があります。 キーをチーム順とメンバー順(チームの書式がリストまたはタプルの場合)/チーム順とプレイヤー名のキー(チームの書式が辞書の場合)のタプルとすることでプレイヤーを特定し、値として参加率を記述します。 言葉で説明してもイメージしづらいと思うので、以下のコードで確認してください。

In[*]

env = trueskill.TrueSkill()

r1 = env.create_rating()
r2 = env.create_rating()
r3 = env.create_rating()
r4 = env.create_rating()

t1 = (r1,r2,)
t2 = (r3,r4,)
teams = (t1,t2,)

weights = {(1,0,): 0.5}

t1,t2, = env.rate(teams, weights=weights)

t1,t2,

Out[*]

((trueskill.Rating(mu=26.738, sigma=7.844),
  trueskill.Rating(mu=26.738, sigma=7.844)),
 (trueskill.Rating(mu=24.131, sigma=8.214),
  trueskill.Rating(mu=23.262, sigma=7.844)))

In[*]

env = trueskill.TrueSkill()

r1 = env.create_rating()
r2 = env.create_rating()
r3 = env.create_rating()
r4 = env.create_rating()

t1 = {'player1': r1, 'player2': r2,}
t2 = {'player3': r3, 'player4': r4,}
teams = (t1,t2,)

weights = {(1,'player3',): 0.5}

t1,t2, = env.rate(teams, weights=weights)

t1,t2,

Out[*]

({'player1': trueskill.Rating(mu=26.738, sigma=7.844),
  'player2': trueskill.Rating(mu=26.738, sigma=7.844)},
 {'player3': trueskill.Rating(mu=24.131, sigma=8.214),
  'player4': trueskill.Rating(mu=23.262, sigma=7.844)})

勝率とマッチング品質

驚くべきことに、この Python 実装には勝率を計算する関数が存在しないらしいです。 以下に示すのはリファレンスにも記載がある Juho Snellman によるスニペットを、環境を指定できるように微加筆したものです。 必要なら利用しましょう。

import itertools
import math

def win_probability(team1, team2, env=None):
    env = env if env else trueskill.global_env()
    delta_mu = sum(r.mu for r in team1) - sum(r.mu for r in team2)
    sum_sigma = sum(r.sigma ** 2 for r in itertools.chain(team1, team2))
    size = len(team1) + len(team2)
    denom = math.sqrt(size * (env.beta * env.beta) + sum_sigma)
    return env.cdf(delta_mu / denom)

In[*]

env = trueskill.TrueSkill(beta=1)

r1 = env.create_rating(mu=30, sigma=0.1)
r2 = env.create_rating(mu=29, sigma=0.1)
t1 = (r1,)
t2 = (r2,)

wp = win_probability(t1, t2, env=env)

wp

Out[*]

0.7591582948828006

どれだけ公平にマッチングできているかを示す指標がqualityメソッドです。 引き分け確率を使用しているようです。 weightsパラメータを与えることができます。

In[*]

env = trueskill.TrueSkill()

r1 = env.create_rating(mu=30, sigma=3)
r2 = env.create_rating(mu=29, sigma=3)
t1 = (r1,)
t2 = (r2,)

q = env.quality((t1,t2,))

q

Out[*]

0.8038743995638264

ショートカット

TrueSkill の利点が多人数対戦に使えることではあれど、ほかのルールにくらべ 1 対 1 のゲームが圧倒的に多いこともまた事実です⁵。 そのため 1 対 1 のゲームで利用できるショートカット関数trueskill.rate_1vs1とtrueskill.quality_1vs1が用意されています。 これらの関数にはレーティングオブジェクトを直接渡すことができます。

In[*]

env = trueskill.TrueSkill()

r1 = env.create_rating()
r2 = env.create_rating()

# (r1,),(r2,), = env.rate(((r1,),(r2,),))
r1,r2, = trueskill.rate_1vs1(r1, r2, env=env)
print(r1)
print(r2)

print()

# (r1,),(r2,), = env.rate(((r1,),(r2,),), ranks=(0,0,))
r1,r2, = trueskill.rate_1vs1(r1, r2, drawn=True, env=env)
print(r1)
print(r2)

print()

# q = env.quality(((r1,),(r2,),))
q = trueskill.quality_1vs1(r1, r2, env=env)
print(q)

Out[*]

trueskill.Rating(mu=29.396, sigma=7.171)
trueskill.Rating(mu=20.604, sigma=7.171)

trueskill.Rating(mu=26.114, sigma=5.678)
trueskill.Rating(mu=23.886, sigma=5.678)

0.5770440474290585

実にスッキリ書けますね。 文字数は増えてるって? カッコの対応がワケワカメになるほうが面倒だしfrom import記法もあるから多少はね。

シミュレーションによるレーティング

「あーなるほどね完全に理解した」したところで、多人数戦ができること、参加率による重み付けができることはわかりましたが、収束が早いというのはどこまで本当なのでしょうか? これは実際に試してみないことにはわかりません。

というわけで、架空のプレイヤーによるゲーム勝敗をシミュレーションし、既存のアルゴリズムと比較してみましょう。

Elo レーティングの理論を利用したシミュレーションの作成

古典的なレーティングアルゴリズムである Elo レーティングは、極めて単純な式でありながらレーティング差から勝率を導き出せるなど、こういったシミュレーションに最適です。

詳細は省きますが、プレイヤー A のレーティングが $R_A$、プレイヤー B のレーティングが $R_B$ であるとき、プレイヤー A の勝率は以下の式で得ることができます。
$$
W_{AB} = \frac{1}{10^{\frac{R_B - R_A}{400}} + 1}
$$
これを利用してシミュレーションスクリプトを書いていきます。

今回想定する架空のプレイヤーは以下のとおりです。

プレイヤー	内部レーティング
Emu	2700
Parado	2700
Niko	2500
Taiga	2000
Hiiro	2000
Kiriya	1700
Kuroto	1700
Poppy	1500

軒並み Elo レーティングの「標準プレイヤー」値である 1500 以上ですが、これはファンへの配慮です。 あくまでも重要なのはレーティングの差なので上位プレイヤーから下位プレイヤーまでがまんべんなく存在していなくても問題ないと思います。 ただ、この結果を Elo レーティングによりレーティングしたとき、彼らの中での標準を 1500 とするため、設定した内部レーティングよりも低い結果が出るものと思われます⁶。

実際に勝敗のシミュレーションを書いたものがこちらのコードになります。 ランダムに二人ごと対戦させ、これを 500 回繰り返します。 引き分けはありません。 結果は CSV ファイルで出力されます。

結果はこんな感じでした。

index	Match 1 Winner	Match 1 Loser	Match 2 Winner	Match 2 Loser	Match 3 Winner	Match 3 Loser	Match 4 Winner	Match 4 Loser
1	Emu	Parado	Taiga	Hiiro	Niko	Poppy	Kiriya	Kuroto
2	Emu	Kuroto	Poppy	Hiiro	Kiriya	Taiga	Niko	Parado
3	Emu	Niko	Kiriya	Kuroto	Hiiro	Poppy	Parado	Taiga
4	Parado	Poppy	Emu	Taiga	Niko	Kuroto	Hiiro	Kiriya
5	Niko	Poppy	Hiiro	Kiriya	Taiga	Kuroto	Parado	Emu
...	...	...	...	...	...	...	...	...
496	Kiriya	Taiga	Niko	Kuroto	Hiiro	Poppy	Emu	Parado
497	Kiriya	Poppy	Emu	Hiiro	Niko	Kuroto	Parado	Taiga
498	Parado	Taiga	Niko	Kiriya	Emu	Kuroto	Hiiro	Poppy
499	Emu	Kuroto	Niko	Taiga	Parado	Kiriya	Hiiro	Poppy
500	Kuroto	Poppy	Parado	Hiiro	Emu	Niko	Kiriya	Taiga

「Emu が負けるありえない話し!!」とか言い出す神ファンの存在をのぞけば、問題なさそうな結果です。 実際のデータはこちらになります。

Elo レーティングで実験

今回比較対象として、Elo レーティングによるレーティングを行ってみます。 こちらのコードを実行し、ゲームごとのレーティングの推移を確認してみます。 Elo レーティングで使用されるレーティング変動の激しさを示す値 $K$ は、最初の 20 ゲームでは 32 を、それ以降は 24 を使用します。

結果は以下のようになりました。

プレイヤー	最終レート
Emu	2006.988
Parado	2062.846
Niko	1821.448
Taiga	1449.791
Hiiro	1398.873
Kiriya	1168.829
Kuroto	1190.664
Poppy	900.5602

とはいえ、前述したように設定値より低く見積もられていますので、Emu との差で見てみましょう。

プレイヤー	差	差(真)
Emu	+0.000	+0.000
Parado	+55.858	+0.000
Niko	-185.540	-200.000
Taiga	-557.197	-700.000
Hiiro	-608.115	-700.000
Kiriya	-838.159	-1,000.000
Kuroto	-816.324	-1,000.000
Poppy	-1,106.428	-1,200.000

確かにある程度設定値にしたがった値になっているようです。

一方で推移は以下のようになっています。 生の値はガクガクして見づらいので 5 点移動平均⁷で平滑化して示しています。

これを見ると、完全に収束しきっておらず、まだまだ伸びしろがあるように思います。 先程の差も、実際のところ設定値まで広がりきっていません。

TrueSkill で実験

つづいて本命の TrueSkill を試してみましょう。こちらのコードを実行します。

結果は以下のようになりました。

プレイヤー	最終顕在化レート	$\mu$	$\sigma$
Emu	32.59123	35.70087	1.036545
Parado	34.08617	37.20115	1.038326
Niko	29.41995	32.34306	0.974371
Taiga	21.1425	23.87729	0.911596
Hiiro	20.62174	23.38207	0.92011
Kiriya	15.11853	17.90101	0.927492
Kuroto	15.4544	18.2791	0.941565
Poppy	9.310735	12.44336	1.044209

おおむね順序通り、内部レーティングを同じ値にしたプレイヤーはだいたい同じ値に落ち着いています。 適切にレーティングできていると言っていいでしょう。

一方で推移は以下のようになりました。 こちらも 5 点移動平均で平滑化して示しています。

見てのとおり、Elo レーティングにくらべて爆速でアプリをデプロイする方法収束しているのが分かります。 最初の数十ゲームでもうある程度の位置が定まっていますね。 これは確かに「収束が早い」と謳うだけのことはあります。

日本プロ野球リーグをレーティングしてみる

いつまでも夢に逃げているとワカメを冷まされるので、いよいよ現実の世界に目を向けてみましょう。

ちょうど最近プロ野球の公式戦が終了いたしましたので、今年度の公式戦データ(セパ交流戦含む)から 12 球団をレーティングしてみましょう。 そんなデータ出したら殺されるのでは?

ちなみにセ・リーグとパ・リーグの 2019 年度公式順位は以下のようになっています。

順位	チーム	勝	敗	分	勝率	ゲーム差
1	巨人	77	64	2	.546	0.0
2	DeNA	71	69	3	.507	5.5
3	阪神	69	68	6	.504	6.0
4	広島	70	70	3	.500	6.5
5	中日	68	73	2	.482	9.0
6	ヤクルト	59	82	2	.418	18.0

順位	チーム	勝	敗	分	勝率	ゲーム差
1	西武	80	62	1	.563	0.0
2	ソフトバンク	76	62	5	.551	2.0
3	楽天	71	68	4	.511	7.5
4	ロッテ	69	70	4	.496	9.5
5	日本ハム	65	73	5	.471	13.0
6	オリックス	61	75	7	.449	16.0

野球見ていないマンゆえ贔屓球団とかはないので、この結果を見ても特に感慨がわかないのが残念です。 あえていうなら、DeNA になる前の横浜はひどかったという噂を聞いたことはあるので、いまは上位に食い込んでいるのが面白いなってくらいですかね……。 一応故郷に最も近いのは日ハムです。

ただ、勝率を見ていくと、もちろん順位が付く程度の差はできていますが、たとえばチェスの強い人と初心者レベルの差という程のものはできていない、割と拮抗している⁸リーグだなあという感じですね。 これはうまくレーティングできるか不安です。

なにはともあれレーティングしてみましょう。 2019 年度の公式戦のデータはこちらに用意しました。 NPB 公式サイトからシコシココピペして作りましたがクッソ面倒でした。 API とか探したらありませんかね?

環境はほぼ初期のまま、引き分け確率のみ過去 3 年度の合計試合数と合計引き分け数から算出して使用します。 こちらにコードを置いておきます。 引き分け数 CSV は一応ここにおいておきますがこれいる?

結果は以下のようになりました。

チーム	TrueSkill	公式順位	勝	敗	分	勝率
西武	23.932	パ 1 位	80	62	1	0.563
ソフトバンク	23.083	パ 2 位	76	62	5	0.551
楽天	22.884	パ 3 位	71	68	4	0.511
阪神	22.801	セ 3 位	69	68	6	0.504
巨人	22.726	セ 1 位	77	64	2	0.546
ロッテ	22.714	パ 4 位	69	70	4	0.496
DeNA	22.492	セ 2 位	71	69	3	0.507
中日	22.456	セ 5 位	68	73	2	0.482
広島	22.362	セ 4 位	70	70	3	0.500
オリックス	22.067	パ 6 位	61	75	7	0.449
日本ハム	21.979	パ 5 位	65	73	5	0.471
ヤクルト	21.395	セ 6 位	59	82	2	0.418

野球を見ている皆さん的にはこれは直感に即した結果でしょうかね?

交流戦は全体に比べると数が少ないので、セパ横断順位の信憑性はあんまないと思います。 しかし、リーグ単位で見ても公式順位との入れ替わりがいくつかあるのが見て取れます。 特に阪神はセ・リーグ 3 位にもかかわらずセ・リーグ優勝の巨人を抜いてトップに躍り出ています。

「単に勝率のみを見るのではなく、どれほど強い相手に勝ったかを見る」のがレーティングですが、ランダムマッチングのゲームと違い対戦機会が計画されて平等になっているのがプロ野球ですから、「弱い相手ばっかり選んでるから」みたいなのは発生し得ないはずです。

考えられる可能性としては……単にコピペしただけでちゃんと試合結果を精査したわけではないので推測なのですが、「巨人や DeNA のような成績の良いチームに対する勝率が良かった説」があります。 逆に巨人や DeNA の下位陣相性が悪い説もあります。 やきうのお兄ちゃんたち、実際はどうでした?

レーティングは極めて単純な仮定をおいているところがあり、「A さんが B さんに 2 勝 1 敗、B さんが C さんに 3 勝 2 敗ならば、A さんは C さんに $\frac{2}{1} \times \frac{3}{2} = 3$ 倍の勝利を上げる(4 戦して 3 勝 1 敗)だろう」みたいな感じです。 しかし実際の勝負には相性があることがしばしばあり、C さんが A さん相手だと妙に強かったりするわけです。 こういった影響を小さくするためにも、レーティングでは「なるべく多くのプレイヤープールで、かつ同じ相手とばかり戦わない」ことが精度上昇のためにも重要だったりするわけです。 今回とりあげたプロ野球は、プレイヤープール 12 球団と少なく、さらに交流戦があるとはいえリーグが 2 種に分断されているのであまり良い例ではなかったかもしれませんね。 ぶっちゃけ IIDX のアリーナモードとかレーティングしたいけど KONAMI しか試合データ持っていない。

ちょっと脇道にそれ過ぎましたので、推移のグラフでも貼っておきましょう。 例によって 5 点移動平均です。

……団子状態過ぎてよくわかりませんね。 $17 \leq y \leq 25$ の範囲だけ拡大してみましょう。

まだ多少見えづらいですがだいぶマシかな? つかチームカラーが似たりよったり過ぎて区別つきづらいんですが……。

パッと見て目立つのは 90 試合前あたりから広島が思いっきり落ち込んでいるところですね……。 今年広島は交流戦で「こいついっぱい貯金持ってるンゴｗｗｗ」されたらしいのでそれかもしれません。

あと、横浜は今年セ 2 位となりましたが、立ち上がりはかなり遅かったようです。

最後に、各チームごとの勝率も貼っておきます。 マジで殺されるぞ。

○＼●	巨人	ヤクルト	DeNA	中日	阪神	広島	日本ハム	楽天	西武	ロッテ	オリックス	ソフトバンク
巨人	-	58.7%	51.6%	51.9%	49.6%	52.5%	55.0%	49.1%	42.0%	50.2%	54.5%	47.8%
ヤクルト	41.3%	-	42.9%	43.1%	40.9%	43.7%	46.3%	40.4%	33.7%	41.5%	45.8%	39.1%
DeNA	48.4%	57.1%	-	50.3%	48.0%	50.9%	53.4%	47.5%	40.4%	48.6%	52.9%	46.2%
中日	48.1%	56.9%	49.7%	-	47.7%	50.6%	53.2%	47.2%	40.2%	48.3%	52.7%	45.9%
阪神	50.4%	59.1%	52.0%	52.3%	-	52.9%	55.4%	49.5%	42.4%	50.6%	54.9%	48.2%
広島	47.5%	56.3%	49.1%	49.4%	47.1%	-	52.6%	46.6%	39.6%	47.7%	52.1%	45.3%
日本ハム	45.0%	53.7%	46.6%	46.8%	44.6%	47.4%	-	44.1%	37.2%	45.2%	49.5%	42.8%
楽天	50.9%	59.6%	52.5%	52.8%	50.5%	53.4%	55.9%	-	42.9%	51.1%	55.5%	48.7%
西武	58.0%	66.3%	59.6%	59.8%	57.6%	60.4%	62.8%	57.1%	-	58.2%	62.4%	55.8%
ロッテ	49.8%	58.5%	51.4%	51.7%	49.4%	52.3%	54.8%	48.9%	41.8%	-	54.3%	47.6%
オリックス	45.5%	54.2%	47.1%	47.3%	45.1%	47.9%	50.5%	44.5%	37.6%	45.7%	-	43.3%
ソフトバンク	52.2%	60.9%	53.8%	54.1%	51.8%	54.7%	57.2%	51.3%	44.2%	52.4%	56.7%	-

いま示したのは既知のデータに基づく勝敗比ではなく、未知の試合に対する勝利確率ということになります。 まあ、あんまり参考にならないとは思いますが……。 ともあれ、レーティング最上位の西武対レーティング最下位のヤクルトでもそこまで絶望的な数値になっていないというところを見ると、やっぱり拮抗した良いリーグです。 野球見てみようかな?

よくばりセット

今回使用したスクリプトやデータを固めたものをここに置いておきました。

Notebook 形式なので JupyterLab などをご用意して頂く必要があります。

まとめ

いかがでしたか?(定型文)

今回使用したのは Python 実装でしたが、どうも様々な言語にフォークされているようなので、あなたが使いたい言語での実装も見つかるかもしれませんし見つからないかもしれません(あいまい)。

私が今回の記事で伝えたかったのは、「既存の Elo レーティングとかにくらべて TrueSkill の理論は、腰を据えて読まないとわからないしそもそもある程度の前提知識がないと厳しいのは確かだけど、使うだけなら脳死でメソッド叩くだけで使えるし性能もいいからめっちゃ使ってほしい」ということです。

Elo を利用した記事は Qiita にも多い⁹ので、これから TrueSkill でレーティングする記事も増えてほしいと思います。

あと、レーティングといえば Glicko とかも Elo より性能よくて自分で実装するのも難しくないのでオヌヌメです。

Elo 以外のレーティングもっと流行れ!

---

1. チームとはいっても、テレビで放送されているプロスポーツのような「チームの構成員がほぼ一定で、ある選手がゲーム単位の短いスパンでチームを変更しない」チームの場合、チームそのものをプレイヤーとして扱うのが理にかなっています。 ここでいうチーム戦とは、例えば Splatoon のようなマッチングしたプレイヤーをランダムにチームに振り分けて行われるようなゲームに向いています。 ↩

2. $\phi_{\mu, \sigma^2}$ は確率密度関数一般のことではなく正規分布のことなのですが、わかりやすさ重視でここにねじ込んでしまいました。 YURUSHITE。 ↩

3. 具体的には 99.73% 程度。 ↩

4. 実際には、更新前にごく少量の σ が与えられるため、常に減少し続けるわけではありません。 成長やブランクによるスキルの変化を許容できるということだと思います。 ↩

5. そもそもレーティングの研究はチェスから始まりました。 ↩

6. そして、これがレーティングという仕組みが本質的に抱えるどうしようもない弱点でもあります。 たとえば、あるリーグのプレイヤーがもつレーティングと、それとは別のリーグに参加しているプレイヤーがもつレーティングは、レーティングの仕組みが同じだったとしても比較できません。 ↩

7. 今回のスムージングでは、端を端点の値で埋めてからmode='valid'で畳み込んでいます。 mode='same'でも良いのでしょうが、まあ「しゅみです」ということで。 ↩

8. この辺は野球見ている人と温度差があると思います。 レーティング視点でということです。 ↩

9. 車輪の再発明系の記事だと TrueSkill は難しいと思いますが。 ↩

みなさんは、IRC(Internet Relay Chat)というワードを聞いたことはあるでしょうか。
今IRCを使っているという人は少ないかもしれません(私の周りにも一人もいません)が、今回はこのIRCについて書いていきます！
(すでにIRCクライアントのコードが載っているQiitaの記事がありましたが、気になる部分がいくつかあったのと、Python2のコードだったので、改めて記事にしようと思いました。)

IRC(Internet Relay Chat)とは

IRCとは、サーバを介してクライアントとクライアントが会話をする枠組みで、インスタントメッセンジャーのプロトコルの1つに分類されます。
また、ユーザがクライアントソフトを実行し、サーバは情報の伝達を行うだけなので、サーバの負担が軽くチャットが高速に行えます。
チャットに特化したプロトコルなんですね！

なぜIRCの利用は衰退したの？

チャットアプリなんて、今はたくさんありますよね。(Facebook,Twitter,Discord,Slack,...etc.)
これらSNSを利用するほうが、IRCを利用するよりも便利で多機能で使いやすいわけですね。。。
また、SSL/TLSを利用せずにIRCを利用すると、暗号化していないため、第三者に会話の内容を盗み見られる可能性があります。(もちろんSSL/TLSを利用すれば安全です！)
LINEなどのサービスなら、暗号化などは自動でやってくれるから、そんなの普段は意識しないですよね。。

IRCを利用するメリットは？

それだけ聞くとIRCを利用するメリットはあまりないように感じるかもしれませんが、そんなことはありません！
まず、IRCは、データ通信に関するプロトコルが簡素で、かつオープンであるため、ユーザ・クライアントに用いるソフトウェアの開発が簡単なんです！今回の記事でも、簡易クライアントをたった100行程度で作っちゃいます！
それだけではなく、IRCサーバを立てるのも、ngircdというものを利用すれば超簡単なんです！(今回はこれには触れません。)

今IRCは何に使われてるの？

では今、IRCは何に使われているんでしょうか。
それは、、、ボットネットの通信手段です。(誤解を招きそうなので言っておきますが、もちろん正規の普通の使い方でもたくさん利用されています。)
ボットネットとは、ボットウイルスに感染した端末で構成されるネットワークです。ボットウイルスとは、マルウェアの一種で、感染すると端末が悪意を持った第三者に遠隔操作され、不正アクセスやDDoS攻撃などの犯罪行為に使われてしまいます。
その感染端末と攻撃者の通信手段に、IRCはよく使われています。。。
自分のパソコンをポートスキャンしたり、通信のログを見て、「IRC」という文字を見つけた場合は、ボットウイルスに感染している可能性があるかもしれません。。。

有名なIRCクライアントソフト

IRCクライアントソフトはたくさんありますが、ここでは中でも有名(人気)なもの(すべて無料)を紹介します。

LimeChat

おそらく1番有名(人気)なIRCクライアントソフトです。
Windowsの場合はLimeChat 2.x、Macの場合はLimeChat for OSXを利用します。
現行IRCクライアントの中では最も安定しています。ヘルプが充実しているので初心者にもおすすめです。

HexChat

XChat(後述)ベースのソフトウェアです。
Windows用、Mac用、Unix系システム用があり、ほぼすべてのOSで利用できます。
海外のさまざまなネットワークが標準設定で入っています。

XChat

Linuxのほとんどのディストリビューションに標準で入っているIRCクライアントです。
Windows用、Mac用、Unix系システム用があり、ほぼすべてのOSで利用できます。
海外のさまざまなネットワークが標準設定で入っています。

Python3で簡易クライアントを作る

では早速IRCクライアントを作ってみましょう！
IRCクライアントを作るには、まずIRCのコマンドというものについて知る必要があります。

IRCのコマンドについて

IRCにおいて、サーバとクライアントは互いにメッセージを送信し合います。
また、メッセージが正しいコマンドを含んでいた場合、クライアントは仕様通りのリプライを期待します。
つまり、あるコマンドを含むメッセージに対するサーバからの応答メッセージの形式は決まっているわけですね。
このIRCコマンドはたくさんありますが、今回は、後で作成する簡易クライアントが実装している、最も基本的なコマンド8つについてのみ説明します。(コマンドは正規表現で示します。)

PASSコマンド

PASS <password>

IRCサーバにパスワードが設定されている場合は、ユーザは接続を開始しようとする前(NICK/USERの組み合わせを送る前)にPASSコマンドを送る必要があります。

NICKコマンド

NICK :<nickname>

ユーザにニックネームを設定したり、今のニックネームを変更したりするのに使います。

USERコマンド

USER <username> <hostname> <servername> <realname>

新しいユーザのユーザ名やホスト名、本名を指定するために、接続のはじめに使われます。

QUITコマンド

QUIT (:<Quit Message>)??

クライアントセッションを終了します。サーバはクライアントにERRORメッセージを送ることでこれを承認します。

JOINコマンド

JOIN :<channel> ("," <channel>)* (<key> ("," <key>)*)??

ユーザによって、特定のチャンネルに接続するリクエストを行うために使用されます。

PRIVMSGコマンド

PRIVMASG <msgtarget> :<text to be sent>

ユーザ間のプライベートメッセージを送るのに使います。また、メッセージをチャンネルに送るのにも使われます。msgtargetは通常、メッセージの受け取り手のニックネームか、チャンネル名です。
誰かがメッセージを送った場合は、サーバから次の形式でメッセージが送られてきます。

送られてくるとき

:<A's nickname>!(~<A's username>)??@<A's gateway-name> PRIVMSG <msgtarget> :<text to be sent>

PINGコマンド

PING :<server1> (<server2>)??

ネットワーク上のアクティブなクライアントやサーバが、実際につながって動いているかどうかを確認するのに使われます。
サーバは、ネットワークからほかのアクションが届かない場合、一定間隔でPINGメッセージを送出します。
指定時間内にPINGメッセージに返答がない場合、その接続は閉じられます。
PINGメッセージを受け取ったら、server1(PINGメッセージを送ったサーバ)へのリプライとして、可能な限り早く正しいPONGメッセージを送らなくてはなりません。server2パラメータが指定されていれば、それがpingのターゲットとなり、メッセージはそこに転送されます。

PONGコマンド

PONG <server1> (<server2>)??

PINGメッセージへのリプライです。server2パラメータが指定されれば、メッセージはそこに転送されなくてはなりません。

IRCクライアント(Python3コード)

import sys, socket, os, signal

HOST = "X.X.X.X"
PORT = 6667 #IRCサーバでは一般的に6667番ポートが使われることが多い

BUF_SIZE = 1024

class IRC(object):
    def __init__(self):
        self.server = socket.socket(socket.AF_INET, socket.SOCK_STREAM) #socketオブジェクトの生成(TCP)

    def connect(self, host, port):
        self.server.connect((host, port)) #接続

    def login(self, password, nickname, username, realname, hostname = "hostname", servername = "*"):
        if password is not None: #中にはパスワードがいらないサーバもある
            pass_message = "PASS " + password + "\n" #PASSメッセージ
            self.server.send(pass_message.encode('utf-8')) #送信
        nick_message = "NICK " + nickname + "\n" #NICKメッセージ
        user_message = "USER %s %s %s :%s\n" % (username, hostname, servername, realname) #USERメッセージ
        self.server.send(nick_message.encode('utf-8')) #送信
        self.server.send(user_message.encode('utf-8')) #送信

    def join(self, channel):
        join_message = "JOIN " + channel + "\n" #JOINメッセージ
        self.server.send(join_message.encode('utf-8')) #送信

    def pong(self, server1, server2 = None):
        pong_message = "PONG %s %s" % (server1, server2) #PONGメッセージ
        pong_message += "\n"
        self.server.send(pong_message.encode('utf-8')) #送信

    def privmsg(self, channel, text):
        privmsg_message = "PRIVMSG %s :%s\n" % (channel, text) #PRIVMSGメッセージ
        self.server.send(privmsg_message.encode('utf-8')) #送信

    def quit(self):
        self.server.send(b"QUIT :bye!") #QUITメッセージ送信

    def handle_privmsg(self, prefix, text):
        print("\r" + prefix + ">" + text + "\n>", end="") #受信したPRIVMSGメッセージを処理、表示

    def wait_message(self):
        while(True):
            msg_buf = self.server.recv(BUF_SIZE) #受信
            msg_buf = msg_buf.decode('utf-8').strip()
            ## ここからメッセージ処理 ##
            prefix = None
            if msg_buf[0] == ":":
                p = msg_buf.find(" ")
                prefix = msg_buf[1:p]
                msg_buf = msg_buf[(p + 1):]

            p = msg_buf.find(":")
            if p != -1: #":"から始まるパラメータがまだあった場合
                last_param = msg_buf[(p + 1):]
                msg_buf = msg_buf[:p]
                msg_buf = msg_buf.strip()

            messages = msg_buf.split()
            ## ここまで ##

            command = messages[0] #コマンド名
            params = messages[1:] #今回は無視

            if command == "PING":
                self.pong(last_param) #PINGが来たらすぐPONGを返す
            elif command == "PRIVMSG":
                text = last_param #PRIVMSGコマンドで送られてきたメッセージ
                self.handle_privmsg(prefix, text)


    def client_interface(self, channel, prompt = ">"):
        while(True):
            line = input(prompt)

            if line == "quit":
                self.quit()
                break

            self.privmsg(channel, line)


def main():
    password = "password"
    nickname = "nickhoge"
    username = "usr"
    realname = "realname"
    channel = "#test_channel"

    irc = IRC()
    irc.connect(HOST, PORT)
    irc.login(password, nickname, username, realname)
    irc.join(channel)

    pid = os.fork() #子プロセス生成

    if(pid == 0): #os.fork()は、子プロセスでは0を返す
        irc.wait_message()
    else:
        irc.client_interface(channel)
        os.kill(pid, signal.SIGTERM) #子プロセスをkill(これをしないと子プロセスがゾンビプロセスになる)

if __name__ == "__main__":
    main()

macで実行してみた

実際にLimeChatユーザと会話してみました。

LimeChat側からみるとこんな感じです。

いい感じですね！

参考

https://qiita.com/mmttt202/items/d182c6f27f466c923fea
http://jbpe.tripod.com/rfcj/rfc2812.j.sjis.txt
https://www.friend-chat.jp/irc.html
https://wikiwiki.jp/2chIRC/IRC%E3%82%AF%E3%83%A9%E3%82%A4%E3%82%A2%E3%83%B3%E3%83%88%E7%B4%B9%E4%BB%8B
https://ja.wikipedia.org/wiki/Internet_Relay_Chat
https://murashun.jp/blog/20190215-01.html

楽しみにしていた「macOS　Catalina」にアップデート。
あれっ、、アップデートしたらAnacondaが使えなくなってる。。。
そんな方への対処法をまとめました。

なんで使えなくなった

Catalinaではルートディレクトリへの書き込みをセキュリティ上推奨しないようです。
一方、Anacondaはディフォルトではルートに生成されています。
OSアップグレードの際にルートディレクトリから丸ごと移動されてしまうため、使えなくなってしまいます。

追い出されたanacondaは以下に移動されています。
/Users/Shared/Relocated Items/Security/anaconda3
デスクトップにエイリアス「場所が変更された項目」が生成されているのでそちらからも辿れます。

対応方法

詳しくは以下の海外サイト、こちらを参考にしました。
https://superuser.com/questions/1456518/anaconda-not-working-in-mac-catalina

対応Step

1 追い出されたanacondaディレクトリをホームディレクトリに移動

#追い出されたanaconda
/Users/Shared/Relocated Items/Security/anaconda3 
#ホームディレクトリに移動
/Users/[your_user_name]/anaconda3

2 condaを使用可能にする

#binディレクトリに移動
cd /Users/[your_user_name]/anaconda3/bin/

#condaファイルのshebangを変更
vi conda
#1行目
 #!/anaconda3/bin/python
#↓以下に変更
 #!/Users/[your_user_name]/anaconda3/bin/python.

#condaファイルを保存

3 condaの初期化
初期化することで、依存ファイルが最適化され、shellにanacondaのパスが通ります。

conda init bash　（conda init zsh）
#(catalinaのディフォルトshはzsh。私は取り急ぎbashでしましたが、zsh移行済みの方は、zshで初期化してください）

#以下のステータスが走ります。
modified 〜〜〜
・・・
・・・
==> For changes to take effect, close and re-open your current shell. <==

4 確認
~/.bash_profileにパスを通す記述が自動追記されます。

# >>> conda initialize >>>
# !! Contents within this block are managed by 'conda init' !!                   

#<パスを通すコマンドが自動記述される>

unset __conda_setup
# <<< conda initialize <<<

5　読み込み、正常使用を確認
terminalを再度立ち上げるか、
source ~/.bash_profileで設定を再読み込み。

これでエラーがなければ、anacondaは正常に使えるようになっていると思います。
anaconda経由のpythonや、ライブラリ、condaも使用できるようになります。

まだ、Catalinaがリリースされてから日が浅いので、一部不具合の可能性もあります。心配な方は少し経ってからの移行をおすすめします。

6 anaconda-navigatorの対応
私はanaconda-navigatorを利用してませんが、参考サイトに対応方法ありましたので、ついでに対応してみました。

• anacondaディレクトリに移動
  /Users/[your_user_name]/anaconda3/
• Anaconda-Navigator.app を右クリック
  パッケージの内容を表示を選択
• Contents/MacOS/run.sh を開く
• 以下の部分を全て書き換える

#変更前
/anaconda3
↓
#変更後
/Users/[your_user_name]/anaconda3

ファイルを保存する。
これで「Anaconda Navigator」は使用できるようになりました。

7 番外
私の場合は、「jupyter notebook」も使えなくっていました。
以前、conda経由でインストールしたのですが、今回のOSアップグレードに関係あるかもしれません。
私の場合ですが、アップグレードしたら使えるようになりました。

#conda自体をアップデート
conda update conda
#jupyter notebookのアップデート
conda update jupyter

ただ、環境により人それぞれ対応方法は異なるかもしれません。

　追伸

試してないけど、違うやり方です。
https://medium.com/@justinaugust/upgrading-to-osx-catalina-as-an-anaconda-user-2e71db194764

ざっくりこんな方法みたいです。
1. セキュリティディレクトリに避難された、anacondaディレクトリをホームディレクトリに移動するのは同じ。
2. anacondaの公式ページからダウンロードインストーラーをダウンロード
3. Terminalで　インストール作業
./<name of anaconda installer> -u
ホームディレクトリを指示してanaconda3をインストール

要はインストールすることで設定を上書きするやり方のようです。

おわりに

以上、macOS Catalinaを入れた際のAnacondaの対応方法をまとめてみました。
OSリリースしてまだ日が浅いこともあり、まだ他の不具合が出てくる可能性もありますが、
ひとまず私の環境ではAnaconda/pythonは使えております。

稚拙なまとめでしたが、みなさんのお力になれると幸いです。

SAMPLE

REFERENCE

folium 事始め
Folium: Pythonでデータを地図上に可視化

DATA

国土交通省国土政策局のページから町丁目レベル位置参照情報とをダウンロードする。
=>ここではテーブル名[locate_lt_lg]としてSQL Serverに格納。

PYTHON

ここでは
・千代田区
・大字・字・丁目区分コード3
のデータのみを利用する。

import pypyodbc

connection_string ='Driver={SQL Server Native Client 11.0};Server=ServerName;Database=DatabaseName;Uid=UserName;Pwd=Password;'
connection = pypyodbc.connect(connection_string)
SQL = """
  SELECT todofuken_name,shichoson_name,choson_name,lt,lg
  FROM [locate_lt_lg]
  WHERE chomoku_code = '3'
     AND shichoson_name = '千代田区'
"""

cur = connection.cursor()
cur.execute(SQL)

rows = cur.fetchall()
ltlg_data = []
for row in rows:
    ltlg_data.append(row)

cur.close()
connection.close()

中心座標を東京駅に設定して丁目ごとの代表座標をマーカーとして設置。htmlとして保存する。

import folium

copyright_map = '&copy; <a href="http://osm.org/copyright">OpenStreetMap</a> contributors'

map0 = folium.Map(location=[35.681167, 139.767052],
                  attr=copyright_map,
                  zoom_start=15)

states = []
for row in ltlg_data:
   state = {}
   state['todofuken_name'] = row[0]
   state['shichoson_name'] = row[1]
   state['choson_name'] = row[2]
   state['lt'] = row[3]
   state['lg'] = row[4]
   states.append(state)

for state in states:
   folium.Marker(
            [state['lt'], state['lg']],
            popup=state['shichoson_name']+state['choson_name'],
    ).add_to(map0)


map0.save('sample0.html')

最近、Pythonの新たな実装追っかけを怠っていたので、この時間でお勉強してみた。
吉祥寺でゆるく仕事帰りの勉強会をやっていて、そこで調べたことを記録したノートです。
https://kichipy.connpass.com/event/149625/

PEP 572, Assignment expressionsについて

参考 syntax-and-assignment-expressions-what-and-why

:= という表記のAssignment expressions（代入演算子と言ったら良いのかな）が追加されました。例えば、 if や while の条件式の条件内で、変数に値を代入しながら処理を実行したりする時に使うと良さげです。

まずは、、、、

if statement（if条件式）

今まで

match = pattern.match(line)
if match:
    return match.group(1)

Assignment expressions（代入演算）にすると

if match := pattern.match(line):
    return match.group(1)

Infinite while statement（whileの条件式）

今まで:

while True:
    data = f.read(1024)
    if not data:
        break
    use(data)

Assignment expressions（代入演算）にすると

while data := f.read(1024):
    use(data)

だいぶコードがシンプルに書けることがわかります。
こんな処理は実際書くことを無いのですが、動作のおさらい。

>>> if a := None:  # a に None を代入して、a を評価
...     print(a)
...
>>> if a := 1:  # a に 1 を代入して、a を評価
...     print(a)
...
1

その他stackoverflowに上がっていた例

stuff = [(lambda y: [y,x/y])(f(x)) for x in range(5)]

が

stuff = [[y := f(x), x/y] for x in range(5)]

こうなります。読みにくい lambda を使わなくてよいので、いい感じです。

inputを以下様に処理もできます。

command = input("> ")
while command != "quit":
    print("You entered:", command)
    command = input("> ")

が

while (command := input("> ")) != "quit":
    print("You entered:", command)

なります。

感想

実践でも使い所は結構ありそうなのと、そこまで可読性が下がる印象はないので、積極的に使っていこうかなと思った次第。