CUBE SUGAR CONTAINER

技術系のこと書きます。

Linux Bridge でタグ VLAN を使う

今回は Linux Bridge (ネットワークブリッジ) でタグ VLAN (IEEE 802.1Q) を使う方法について。 設定が足らなくてだいぶ悩んだのでメモしておく。

使った環境は次の通り。

$ cat /etc/lsb-release 
DISTRIB_ID=Ubuntu
DISTRIB_RELEASE=18.04
DISTRIB_CODENAME=bionic
DISTRIB_DESCRIPTION="Ubuntu 18.04.1 LTS"
$ uname -r
4.15.0-29-generic
$ ip -V
ip utility, iproute2-ss180129
$ sudo apt-cache show iproute2 | grep -i version
Version: 4.15.0-2ubuntu1

下準備

最初に、今回使うパッケージをインストールしておく。

$ sudo apt-get -y install iproute2 tcpdump iputils-ping

Linux Bridge でタグ VLAN を有効にする

準備ができたところで、早速本題に入る。 Linux Bridge でタグ VLAN を扱うには vlan_filtering を有効にする必要がある。 vlan_filteringip(8) で有効にできる。

例えば ip(8) を使ってブリッジを作る場合には、次のように有効にできる。 これは br0 という名前でブリッジを作るときの例。

$ sudo ip link add dev br0 type bridge vlan_filtering 1

また、作ったあとからでも次のように有効にできる。 brctl(8) とか使って作った場合にも、これでいけるはず。

$ sudo ip link add dev br0 type bridge
$ sudo ip link set dev br0 type bridge vlan_filtering 1

ブリッジを作ったら状態を UP に設定する。

$ sudo ip link set br0 up

ここからは、上記の設定が正しく動いていることを確認していく。

アクセスポートを追加する

続いては、タグ VLAN を有効にした Linux Bridge に実際にアクセスポートを追加する。 追加するネットワークインターフェースには veth インターフェースを使う。 また、あとから ping(8) で疎通を確認できるように Network Namespace をつなげておく。

まずは Network Namespace を用意する。

$ sudo ip netns add ns1

続いて veth インターフェースのペアを作成する。

$ sudo ip link add ns1-veth0 type veth peer name br0-veth0

片方の veth インターフェースを、先ほど作った Network Namespace に所属させる。

$ sudo ip link set ns1-veth0 netns ns1

所属させたネットワークインターフェースを使える状態にして IP アドレスとして 192.0.2.1/24 を付与しておく。

$ sudo ip netns exec ns1 ip link set ns1-veth0 up
$ sudo ip netns exec ns1 ip addr add dev ns1-veth0 192.0.2.1/24

veth インターフェースのもう片方は、先ほど作った Linux Bridge に追加しておく。

$ sudo ip link set br0-veth0 master br0
$ sudo ip link set br0-veth0 up

準備ができたので、やっと VLAN の設定をしていく。 その前に、まずは bridge(8) で現状を確認しておく。 デフォルトでは、次のように VLAN ID として 1 が使われるように見えている。 これは、ようするに 802.1Q VLAN を使っていない状態。

$ bridge vlan show dev br0-veth0
port    vlan ids
br0-veth0    1 PVID Egress Untagged

そこで、bridge(8) を使ってポートに VLAN を設定する。 今回は Linux Bridge に追加したポートを使った通信で VLAN ID 100 が使われるように設定してみよう。 代わりに、既存のルールは削除しておく。

$ sudo bridge vlan add vid 100 dev br0-veth0 pvid untagged
$ sudo bridge vlan del dev br0-veth0 vid 1

設定は以下のようになった。

$ bridge vlan show dev br0-veth0
port    vlan ids
br0-veth0    100 PVID Egress Untagged

通信先 (対向) を用意する

ns1 に続いて、もう一つ ns2 という名前で Network Namespace を用意して、先ほどと同じようにセットアップしておく。 これは ns1 から ping(8) を打つときの相手になる。

$ sudo ip netns add ns2
$ sudo ip link add ns2-veth0 type veth peer name br0-veth1
$ sudo ip link set ns2-veth0 netns ns2
$ sudo ip netns exec ns2 ip link set ns2-veth0 up
$ sudo ip netns exec ns2 ip addr add dev ns2-veth0 192.0.2.2/24
$ sudo ip link set br0-veth1 master br0
$ sudo ip link set br0-veth1 up
$ sudo bridge vlan add vid 100 dev br0-veth1 pvid untagged
$ sudo bridge vlan del dev br0-veth1 vid 1

通信を観察する

準備ができたので、Linux Bridge の通信を tcpdump(1) でパケットキャプチャする。

$ sudo tcpdump -nel -i br0
tcpdump: verbose output suppressed, use -v or -vv for full protocol decode
listening on br0, link-type EN10MB (Ethernet), capture size 262144 bytes

パケットキャプチャの準備ができたら、別のターミナルから ns1 から ns2 に向かって ping を打ってみよう。

$ sudo ip netns exec ns1 ping -c 3 192.0.2.2
PING 192.0.2.2 (192.0.2.2) 56(84) bytes of data.
64 bytes from 192.0.2.2: icmp_seq=1 ttl=64 time=0.114 ms
64 bytes from 192.0.2.2: icmp_seq=2 ttl=64 time=0.060 ms
64 bytes from 192.0.2.2: icmp_seq=3 ttl=64 time=0.073 ms

--- 192.0.2.2 ping statistics ---
3 packets transmitted, 3 received, 0% packet loss, time 2036ms
rtt min/avg/max/mdev = 0.060/0.082/0.114/0.024 ms

打ち終わったら tcpdump(1) を実行していたターミナルに戻る。 すると、次のように 802.1Q のデータフレームが観測できているはず。 VLAN ID にも、ちゃんと 100 が使われていることが分かる。

$ sudo tcpdump -nel -i br0
...(略)...
18:20:18.031208 b6:05:b6:b6:c6:e3 > a2:bd:24:29:d8:76, ethertype 802.1Q (0x8100), length 102: vlan 100, p 0, ethertype IPv4, 192.0.2.1 > 192.0.2.2: ICMP echo request, id 1453, seq 1, length 64
18:20:18.031240 a2:bd:24:29:d8:76 > b6:05:b6:b6:c6:e3, ethertype 802.1Q (0x8100), length 102: vlan 100, p 0, ethertype IPv4, 192.0.2.2 > 192.0.2.1: ICMP echo reply, id 1453, seq 1, length 64
...(略)

めでたしめでたし。

Python: 自作ライブラリのパッケージングについて

今回は Python で自作したライブラリなどをパッケージングして、配布できる状態にする方法について書いてみる。

現在の Python では、パッケージングに setuptools というサードパーティ製のライブラリを使うのがデファクトスタンダードになっている。 この setuptools は、pip などを使ってパッケージをインストールするときにも必要になるので、実は気づかずに使っているという場合も多いかもしれない。 また、サードパーティ製といっても PyPA (Python Packaging Authority) というコミュニティが管理しているので準公式みたいな位置づけ。

今回使った環境は次の通り。

$ sw_vers                                      
ProductName:    Mac OS X
ProductVersion: 10.14.2
BuildVersion:   18C54
$ python -V              
Python 3.7.2
$ pip list | grep setuptools
setuptools 40.7.2

パッケージング用のサンプルコードを用意する

まずはパッケージングの対象となるサンプルコードを用意する必要がある。 そこで、今回は greet() という関数を一つだけ備えたパッケージとして example を用意した。 関数 greet() は、呼び出すと標準出力にメッセージをプリントする。

$ mkdir example
$ cat << 'EOF' > example/__init__.py
# -*- coding: utf-8 -*-

__version__ = '0.0.1'


def greet():
    """挨拶をする関数"""
    print('Hello, World!')
EOF

ちょっと紛らわしいんだけど Python においてパッケージというのは、パッケージングされたライブラリのことを意味していない。 ここでは、Python のコードが入ったディレクトリ、くらいの感覚で考えてもらえれば良い。

ちゃんとした Python におけるモジュールとパッケージという概念については以下を参照のこと。 簡単にいうとモジュールは Python ファイル (*.py) で、パッケージは __init__.py という名前のファイルが入ったディレクトリのことを指している。 blog.amedama.jp

上記で作ったパッケージは、今いるディレクトリからであれば Python の処理系からインポートして使うことができる。 これは、カレントワーキングディレクトリに Python のパスが通っているため。

$ python -c "import example; example.greet()"
Hello, World!

しかし、当然のことながら別の場所に移動してしまうと使えなくなる。 これは、先ほどのパッケージに Python のパスが通っていないため。

$ cd /tmp && python -c "import example; example.greet()" 
Traceback (most recent call last):
  File "<string>", line 1, in <module>
ImportError: No module named example

もちろん、明示的にパスを通してしまえば使うことはできる。 とはいえ、毎回こんなことをしたくはないはず。

$ PYTHONPATH=$HOME/workplace/packaging python -c "import example; example.greet()"
Hello, World!

そこで、作成したライブラリに必要なもの一式をまとめた上で、あらかじめパスの通っている場所に配置できるようにする仕組みが、今回扱うパッケージングというわけ。

サンプルコードをパッケージングする

それでは、先ほど作った example をパッケージングしてみよう。 setuptools を使ったパッケージングでは、まず setup.py という名前の Python ファイルを用意する。 このファイルは、セットアップスクリプトと呼ばれる。

$ cat << 'EOF' > setup.py 
#!/usr/bin/env python
# -*- coding: utf-8 -*-

from setuptools import setup


setup()
EOF

セットアップスクリプトの中では setuptools に含まれる関数 setup() を呼び出している。 引数には何も渡しておらず、これは比較的新しい書き方をしている。 古い書き方では setup() にパッケージに関する様々な情報を引数で渡すことになっていた。

では、パッケージに関する様々な情報は、代わりに何処で渡すのかというと setup.cfg という設定ファイルに記述する。 設定ファイルは ini フォーマットっぽい形式で書いていく。 パッケージの基本的な情報は [metadata] というセクションで扱う。

$ cat << 'EOF' > setup.cfg 
[metadata]
name = example
version = attr:example.__version__
author = Your Name
author_email = your-email@example.com
description = example is a example package
long_description = file:README.md
url = https://example.com/your/project/page
license = Apache License, Version 2.0
classifier =
    Development Status :: 1 - Planning
    Programming Language :: Python
    Intended Audience :: Developers
    Operating System :: POSIX
    Programming Language :: Python :: 2
    Programming Language :: Python :: 2.7
    Programming Language :: Python :: 3
    Programming Language :: Python :: 3.5
    Programming Language :: Python :: 3.6
    Programming Language :: Python :: 3.7
    
[options]
zip_safe = False
packages = find:
EOF

なお、この setup.cfg にメタデータなどを記述していく方式は setuptools のバージョン 30.0.3.0 (リリース日 2016年12月8日) から導入された。 この点は、インストール先となる環境の setuptools にもバージョン 30.0.3.0 以降が求められる点に注意が必要になる。 それより前のバージョンをサポートしたいときは、従来通り setup.py の中にメタデータなどを関数の引数として渡す。

setuptools.readthedocs.io

また、setup.cfg では、いくつかの項目において外部のファイルを読み込むこともできる。 例えば、先ほどの例において long_description では README.md の内容を読み込むよう指定している。 なので、ファイルを作っておこう。

$ cat << 'EOF' > README.md 
### What is this?

A long description for example project.
EOF

さて、これでパッケージングの準備ができた。 ディレクトリは、このような状態になっている。

$ ls
README.md   example     setup.cfg   setup.py

試しに pip を使ってパッケージをインストールしてみよう。 この操作は、ライブラリに必要なもの一式をまとめる作業と、インストールする作業を一気にやっているのに等しい。

$ pip install -U .

うまくいけば、次のように pip のインストール済みパッケージの一覧に表示される。

$ pip list | grep example
example    0.0.1  
$ pip show example
Name: example
Version: 0.0.1
Summary: example is a example package
Home-page: https://example.com/your/project/page
Author: Your Name
Author-email: your-email@example.com
License: Apache License, Version 2.0
Location: /Users/amedama/.virtualenvs/py37/lib/python3.7/site-packages
Requires: 
Required-by: 

また、全然別の場所でパッケージをインポートして使うこともできるようになる。 これは、pip によって作成したパッケージがパスの通った場所に配置されたため。

$ cd /tmp && python -c "import example; example.greet()"
Hello, World!

先ほど pip show したときの Location にあった通り、今回は以下のディレクトリにインストールされていた。

$ cat ~/.virtualenvs/py37/lib/python3.7/site-packages/example/__init__.py 
# -*- coding: utf-8 -*-

__version__ = '0.0.1'


def greet():
    """挨拶をする関数"""
    print('Hello, World!')

動作確認が終わったら、一旦パッケージをアンインストールしておこう。

$ pip uninstall -y example

ソースコード配布物 (sdist) にパッケージングする

先ほどはライブラリに必要なもの一式をまとめるのとインストールするのを pip コマンドで一気にやってしまった。 続いては、この工程を別々に分けてやってみよう。

まずは、基本となるソースコード配布物 (sdist) へのパッケージングから。 パッケージングは先ほど作ったセットアップスクリプトを使って sdist コマンドを実行する。

$ python setup.py sdist

すると dist というディレクトリ以下に tarball ができる。 これこそ、できあがったソースコード配布物のファイル。

$ ls dist 
example-0.0.1.tar.gz

このソースコード配布物からパッケージをインストールできる。

$ pip install dist/example-0.0.1.tar.gz
$ pip list | grep example              
example    0.0.1  

確認できたら、また一旦アンインストールしておこう。

$ pip uninstall -y example

Wheel フォーマットでパッケージングする

ソースコード配布物は基本ではあるものの、使うとちょっとめんどくさい場合もある。 具体的には Python/C API を使った拡張モジュールがライブラリに含まれているパターン。 ソースコード配布物では、拡張モジュールをビルドしていないソースコードの状態で同梱する。 すると、インストール先の環境にはそれをビルドするための開発ツール類一式が必要になる。

そこで、最近は次世代のパッケージング形式として Wheel フォーマットというものが使われ始めている。 このフォーマットでは拡張モジュールはあらかじめビルドされた状態で同梱される。 ただし、環境に依存する場合には各環境ごとに Wheel ファイルを用意する必要がある。

実際に Wheel ファイルをパッケージングしてみよう。 これにはセットアップスクリプトで bdist_wheel コマンドを実行する。

$ python setup.py bdist_wheel

すると、次の通り末尾が whl の Wheel ファイルが dist ディレクトリ以下に作られる。

$ ls dist | grep whl$
example-0.0.1-py3-none-any.whl

もちろんこの Wheel ファイルからもパッケージをインストールできる。

$ pip install dist/example-0.0.1-py3-none-any.whl
$ pip list | grep example                        
example    0.0.1

ちなみに、上記の Wheel ファイルは名前を見て分かる通り Python 3 専用としてビルドされている。 もし、Python 2 でも動く場合には setup.cfg[wheel] セクションについて universal フラグを有効にする。

$ cat << 'EOF' > setup.cfg 
[metadata]
name = example
version = attr:example.__version__
author = Your Name
author_email = your-email@example.com
description = example is a example package
long_description = file:README.md
url = https://example.com/your/project/page
license = Apache License, Version 2.0
classifier =
    Development Status :: 1 - Planning
    Programming Language :: Python
    Intended Audience :: Developers
    Operating System :: POSIX
    Programming Language :: Python :: 2
    Programming Language :: Python :: 2.7
    Programming Language :: Python :: 3
    Programming Language :: Python :: 3.5
    Programming Language :: Python :: 3.6
    Programming Language :: Python :: 3.7
    
[options]
zip_safe = False
packages = find:

[wheel]
universal = 1
EOF

この状態でビルドすると Python 2/3 両対応の Wheel ファイルになる。

$ python setup.py bdist_wheel
$ ls dist | grep whl$
example-0.0.1-py2.py3-none-any.whl
example-0.0.1-py3-none-any.whl

動作が確認できたら、また一旦アンインストールしておこう。

$ pip uninstall -y example

パッケージングした成果物を PyPI に登録する

さて、パッケージングできるようになると、実は成果物を PyPI に登録して一般に公開できる。 この作業には PyPA 製の Twine というツールを使うのがおすすめ。 詳しくは、以下を参照のこと。

blog.amedama.jp

依存ライブラリを指定する

さて、話をちょっと戻して、ここからはパッケージングにおける色々なユースケースについて見ていく。

まずは、自作ライブラリが別のライブラリに依存している場合について。 この場合は [options] セクションに install_requires という項目で指定する。 例えば requests を依存ライブラリとして追加してみよう。

$ cat << 'EOF' > setup.cfg 
[metadata]
name = example
version = attr:example.__version__
author = Your Name
author_email = your-email@example.com
description = example is a example package
long_description = file:README.md
url = https://example.com/your/project/page
license = Apache License, Version 2.0
classifier =
    Development Status :: 1 - Planning
    Programming Language :: Python
    Intended Audience :: Developers
    Operating System :: POSIX
    Programming Language :: Python :: 2
    Programming Language :: Python :: 2.7
    Programming Language :: Python :: 3
    Programming Language :: Python :: 3.5
    Programming Language :: Python :: 3.6
    Programming Language :: Python :: 3.7
    
[options]
zip_safe = False
packages = find:
install_requires =
  requests

[wheel]
universal = 1
EOF

この状態でインストールしてみよう。

$ pip install -U .

すると、次のように requests が一緒にインストールされたことが分かる。

$ pip list | egrep "(example|requests)"
example    0.0.1     
requests   2.21.0    

環境に依存した依存ライブラリを指定する

続いては環境によって必要だったり不要だったりする依存ライブラリの指定について。 例えばリレーショナルデータベースを扱うアプリケーションだと、接続先が MySQL なのか PostgreSQL なのかで必要なドライバが変わる。

例として接続先が MySQL の環境なら mysqlclient をインストールする、という状況を想定してみよう。 この場合、[options.extras_require] セクションに環境名と必要なライブラリを指定していく。

$ cat << 'EOF' > setup.cfg 
[metadata]
name = example
version = attr:example.__version__
author = Your Name
author_email = your-email@example.com
description = example is a example package
long_description = file:README.md
url = https://example.com/your/project/page
license = Apache License, Version 2.0
classifier =
    Development Status :: 1 - Planning
    Programming Language :: Python
    Intended Audience :: Developers
    Operating System :: POSIX
    Programming Language :: Python :: 2
    Programming Language :: Python :: 2.7
    Programming Language :: Python :: 3
    Programming Language :: Python :: 3.5
    Programming Language :: Python :: 3.6
    Programming Language :: Python :: 3.7
    
[options]
zip_safe = False
packages = find:
install_requires =
  requests

[options.extras_require]
mysql =
  mysqlclient

[wheel]
universal = 1
EOF

この状態で、環境に mysql を指定してパッケージをアップデートしてみよう。

$ pip install -U ".[mysql]"

すると、mysqlclient がインストールされることが分かる。

$ pip list | grep mysqlclient
mysqlclient 1.4.1     

ソースコード以外のファイルを成果物に含める

続いては、ソースコード以外のファイルを成果物に含める方法について。 実は、デフォルトでは成果物にソースコード以外のファイルは含まれない。 試しに実験してみよう。

まずはパッケージに greet_from_file() という関数を追加する。 これはテキストファイルから読み込んだ内容を使ってメッセージをプリントする関数になっている。

$ cat << 'EOF' > example/__init__.py
# -*- coding: utf-8 -*-

from __future__ import print_function

import os

__version__ = '0.0.1'


def greet():
    """挨拶をする関数"""
    print('Hello, World!')


def greet_from_file():
    """テキストファイルの内容を使って挨拶する関数"""
    module_dir = os.path.dirname(__file__)
    filepath = os.path.join(module_dir, 'message.txt')
    with open(filepath, mode='r') as fp:
        print(fp.read(), end='')
EOF

続いて、表示するメッセージの元となるファイルを用意する。

$ cat << 'EOF' > example/message.txt 
Hello, World!
EOF

ローカルでは、これで動くようになる。

$ python -c "import example; example.greet_from_file()"
Hello, World!

では、インストールした上ではどうだろうか。

$ pip install -U .

試すと、そんなファイルないよと怒られる。 パッケージングした成果物の中にテキストファイルが含まれていないためだ。

$ cd /tmp && python -c "import example; example.greet_from_file()"
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "/Users/amedama/.virtualenvs/py37/lib/python3.7/site-packages/example/__init__.py", line 19, in greet_from_file
    with open(filepath, mode='r') as fp:
FileNotFoundError: [Errno 2] No such file or directory: '/Users/amedama/.virtualenvs/py37/lib/python3.7/site-packages/example/message.txt'

ソースコード以外のファイルを成果物に含めるには、まず setup.cfg[options] セクションにおいて include_package_data の項目を有効にする。

$ cat << 'EOF' > setup.cfg 
[metadata]
name = example
version = attr:example.__version__
author = Your Name
author_email = your-email@example.com
description = example is a example package
long_description = file:README.md
url = https://example.com/your/project/page
license = Apache License, Version 2.0
classifier =
    Development Status :: 1 - Planning
    Programming Language :: Python
    Intended Audience :: Developers
    Operating System :: POSIX
    Programming Language :: Python :: 2
    Programming Language :: Python :: 2.7
    Programming Language :: Python :: 3
    Programming Language :: Python :: 3.5
    Programming Language :: Python :: 3.6
    Programming Language :: Python :: 3.7
    
[options]
zip_safe = False
include_package_data = True
packages = find:
install_requires =
  requests

[options.extras_require]
mysql =
  mysqlclient

[wheel]
universal = 1
EOF

その上で MANIFEST.in というファイルを用意する。 ここに、成果物に含めるソースコード以外のファイルを指定していく。

$ cat << 'EOF' > MANIFEST.in 
include example/message.txt
EOF

この状態で、もう一度試してみよう。

$  pip install -U .
$ cd /tmp && python -c "import example; example.greet_from_file()"
Hello, World!

今度はエラーにならず実行できた。

コマンドを追加する

続いてはパッケージを追加したときに新たにコマンドが使えるようにする方法について。

今回はコマンド用に別の設定ファイルを用意することにした。 まずは setup.cfg[options] セクションに entry_points という項目を作る。

$ cat << 'EOF' > setup.cfg 
[metadata]
name = example
version = attr:example.__version__
author = Your Name
author_email = your-email@example.com
description = example is a example package
long_description = file:README.md
url = https://example.com/your/project/page
license = Apache License, Version 2.0
classifier =
    Development Status :: 1 - Planning
    Programming Language :: Python
    Intended Audience :: Developers
    Operating System :: POSIX
    Programming Language :: Python :: 2
    Programming Language :: Python :: 2.7
    Programming Language :: Python :: 3
    Programming Language :: Python :: 3.5
    Programming Language :: Python :: 3.6
    Programming Language :: Python :: 3.7
    
[options]
zip_safe = False
include_package_data = True
packages = find:
entry_points = file:entry_points.cfg
install_requires =
  requests

[options.extras_require]
mysql =
  mysqlclient

[wheel]
universal = 1
EOF

その上で、指定したのと同じ設定ファイルを用意する。 以下では example-greet というコマンドを実行したとき example パッケージの greet() 関数が呼ばれるようにしている。

$ cat << 'EOF' > entry_points.cfg
[console_scripts]
example-greet = example:greet
EOF

それでは、パッケージを更新してみよう。

$ pip install -U .

すると example-greet というコマンドが使えるようになっている。

$ example-greet
Hello, World!

テストコードを成果物に含めない

パッケージングしたくなるようなコードには、もちろんテストコードを書くことになる。 続いてはテストコードを成果物に含めない方法とともに、具体的なオペレーションを紹介してみる。

まずはテストをする環境向けの依存ライブラリを [options.extras_require] セクションで指定する。 テストフレームワークには pytest を使うことにした。 同時に、テストを成果物に含めないように [options.packages.find] セクションで tests というディレクトリを探索対象から除外する。

$ cat << 'EOF' > setup.cfg 
[metadata]
name = example
version = attr:example.__version__
author = Your Name
author_email = your-email@example.com
description = example is a example package
long_description = file:README.md
url = https://example.com/your/project/page
license = Apache License, Version 2.0
classifier =
    Development Status :: 1 - Planning
    Programming Language :: Python
    Intended Audience :: Developers
    Operating System :: POSIX
    Programming Language :: Python :: 2
    Programming Language :: Python :: 2.7
    Programming Language :: Python :: 3
    Programming Language :: Python :: 3.5
    Programming Language :: Python :: 3.6
    Programming Language :: Python :: 3.7
    
[options]
zip_safe = False
include_package_data = True
packages = find:
entry_points = file:entry_points.cfg
install_requires =
  requests

[options.extras_require]
mysql =
  mysqlclient
testing =
  pytest

[options.packages.find]
exclude =
  tests

[wheel]
universal = 1
EOF

続いてはテストコードを入れるパッケージを用意する。 中身でやっているのは、あまり本質的ではないことなので気にしなくて良いと思う。

$ mkdir tests
$ touch tests/__init__.py
$ cat << 'EOF' > tests/test_example.py   
#!/usr/bin/env python
# -*- coding: utf-8 -*-

import os
import sys

import pytest

from example import greet

if sys.version_info >= (3, 0, 0):
    # Python 3
    from io import StringIO
else:
    # Python 2
    from StringIO import StringIO


def test_greet():
    """挨拶は大事"""
    # 標準出力を入れ替える
    io = StringIO()
    sys.stdout = io

    # 挨拶する
    greet()

    # 挨拶した内容を確認する
    assert io.getvalue() == ('Hello, World!' + os.linesep)


if __name__ == '__main__':
    pytest.main([__file__])
EOF

準備ができたら、テスト環境向けにインストールする。

$ pip install -U ".[testing]"

修正したときの手間などを考えるとパッケージ本体は普段はアンインストールしておいた方が良いかも。

$ pip uninstall -y example

テストを走らせるにはテストランナーを実行する。

$ py.test
====================================== test session starts ======================================
platform darwin -- Python 3.7.2, pytest-4.2.0, py-1.7.0, pluggy-0.8.1
rootdir: /Users/amedama/Documents/temporary/packaging, inifile:
collected 1 item                                                                                

tests/test_example.py .                                                                   [100%]

=================================== 1 passed in 0.02 seconds ====================================

テストが実行できることが分かったので、続いては成果物の中身を確認してみよう。

$ rm -rf dist
$ python setup.py sdist
$ tar xvf dist/example-0.0.1.tar.gz -C dist
$ find dist/example-0.0.1 
dist/example-0.0.1
dist/example-0.0.1/PKG-INFO
dist/example-0.0.1/example.egg-info
dist/example-0.0.1/example.egg-info/PKG-INFO
dist/example-0.0.1/example.egg-info/not-zip-safe
dist/example-0.0.1/example.egg-info/SOURCES.txt
dist/example-0.0.1/example.egg-info/entry_points.txt
dist/example-0.0.1/example.egg-info/requires.txt
dist/example-0.0.1/example.egg-info/top_level.txt
dist/example-0.0.1/example.egg-info/dependency_links.txt
dist/example-0.0.1/example
dist/example-0.0.1/example/__init__.py
dist/example-0.0.1/example/message.txt
dist/example-0.0.1/MANIFEST.in
dist/example-0.0.1/README.md
dist/example-0.0.1/setup.py
dist/example-0.0.1/setup.cfg

成果物の中には tests というディレクトリが含まれていないことが分かる。

まとめ

  • Python のパッケージングには setuptools というライブラリを使う
  • パッケージングするときはセットアップスクリプト (setup.py) を用意する
  • 最近 (>= 30.0.3.0) の書き方では、なるべく別の設定ファイル (setup.cfg) に内容を記述する
  • パッケージングの成果物にはソースコード配布物 (sdist) と Wheel (whl) がある

Python: Unix における python* コマンドと処理系のバージョンについて (PEP394)

今回は Unix ライクなシステムにおける python* コマンドの振る舞いと、処理系のバージョンについて色々と書いてみる。

きっかけは、以下のブログ記事を目にしたため。

rheb.hatenablog.com

上記のブログでは Red Hat Enterprise Linux 8 に搭載される予定の Python と、それを扱うコマンド群について紹介されている。 そして、RHEL8 がそのような仕様に落ち着いた理由として以下の PEP (Python Enhancement Proposal) が挙げられている。

www.python.org

この PEP394 というドキュメントには、Unix ライクのシステムにおいて推奨される python* コマンドの振る舞いなどが記述されている。 そこで、今回は PEP394 に書かれている内容について紹介したい。 先のブログでは、話の切り口が「RHEL8 における仕様」なのでやむを得ないものの、やや説明が不足しており補足したい点もあった。

python* コマンドで実行される処理系のバージョンについて

まず、PEP394 には python* コマンドを実行したときに、どのバージョンの処理系が呼び出されるべきなのか書かれている。 例えば、次のようにコマンド名でバージョンを明示的に指定したときは、そのバージョンの処理系が呼ばれることになる。 もちろん「そのバージョンの処理系がインストールされていれば」という前提のもとで。

$ python2  # => Python 2.x の処理系が使われる
$ python3  # => Python 3.x の処理系が使われる

問題となるのは、バージョンが明示的に指定されていないとき。 一体、どのバージョンの処理系が呼び出されるべきなのだろうか。 その答えは「現時点では Python 2.x の処理系が使われることが推奨されている」になる。

$ python  # => 現時点では Python 2.x の処理系が使われることが推奨されている

現時点では?

太字にした通り、この挙動が推奨されているのは、あくまで現時点 (2019-01) の PEP394 ではという注釈つきになる。 なぜなら PEP394 には Future Changes to this Recommendation (この勧告に対する将来の変更) という項目があるため。

この項目には、次のような記述がある。(カッコ内は拙訳)

This recommendation will be periodically reviewed over the next few years, and updated when the core development team judges it appropriate. As a point of reference, regular maintenance releases for the Python 2.7 series will continue until at least 2020. (この勧告は将来の数年間に渡って定期的に見直されます。そして、コア開発チームが妥当と判断したときに更新されます。参考までに、Python 2.7 系の定常的なメンテナンスリリースは少なくとも 2020 年まで続きます。)

つまり、将来的にはバージョン指定なしの python コマンドで Python 3.x の処理系が使われることを推奨する内容に更新される余地がある。

例外について

ところで、上記の振る舞いは一部のディストリビューションと仮想環境では例外として扱われている。

一部のディストリビューションというのは、例えばこの PEP が勧告される以前から Python 3 移行を始めた、ある意味で優秀なコミュニティもあったりするため。

www.archlinux.org

もう一つの仮想環境というのは venvvirtualenv を使った環境を指している。 仮想環境の中では、環境を作るときに指定されたか、あるいはデフォルトの処理系が python コマンドで使われる。

$ python3 -m venv myvenv 
$ source myvenv/bin/activate
(myvenv) $ python -V
Python 3.7.2

シェバン (Shebang) について

PEP394 には、シェバン (Shebang) に関する記載もある。

シェバンにおいては、特定のバージョンの Python でしか動かないモジュールならバージョンを指定することが推奨されている。 例えば、Python 2 でしか動かないなら以下のようにする。

#!/usr/bin/env python2
...

Python 3 でしか動かないなら、こう。

#!/usr/bin/env python3
...

バージョン指定のないシェバンは、モジュールが両方のバージョンで動く場合に用いることが推奨されている。

#!/usr/bin/env python
...

インタプリタから別のインタプリタを起動する場合

また、subprocess なんかでインタプリタの中から別のインタプリタを呼ぶときは sys.executable の内容を使うのが良い。 ここには自身の処理系のパスが格納されている。

$ python3 -c "import sys; print(sys.executable)"
/usr/local/opt/python/bin/python3.7

特定のコマンド名 (python2, python3, python, etc) を使うと、誤って別の処理系を呼び出したりして思わぬバグが生じる恐れがある。

まとめ

  • 前提: PEP394 の勧告は将来的に変更される余地がある
  • 各コマンドで起動する処理系のバージョン
    • python2 # => Python 2.x
    • python3 # => Python 3.x
    • python # => 今のところ Python 2.x
  • 上記の例外
    • 一部のディストリビューション (ArchLinux など)
    • 仮想環境 (venv, virtualenv)
  • シェバン (Shebang) を書くとき
    • python2 # => Python 2.x のみで動作する
    • python3 # => Python 3.x のみで動作する
    • python # => 両方のバージョンで動作する
  • インタプリタの中から別のインタプリタを起動するとき
    • sys.executable を使う

いじょう。

Python: XGBoost を使ってみる

XGBoost (eXtreme Gradient Boosting) は勾配ブースティング決定木 (Gradient Boosting Decision Tree) のアルゴリズムを実装したオープンソースのライブラリ。 最近は、同じ GBDT 系のライブラリである LightGBM にややお株を奪われつつあるものの、依然として機械学習コンペティションの一つである Kaggle でよく使われている。 今回は、そんな XGBoost の Python バインディングを使ってみることにする。

使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.14.2
BuildVersion:   18C54
$ python -V
Python 3.7.2

下準備

最初に、下準備として XGBoost と必要なパッケージを一通りインストールする。

そのために、まずは gcc をインストールしておく。 これは XGBoost が並列処理に OpenMP を使っている関係で必要になる。

$ brew install gcc@7

続いて XGBoost と、それ以外で今回使うパッケージをインストールする。

$ pip install xgboost scikit-learn matplotlib

乳がんデータセットを分類してみる

まずはハローワールド的な例として乳がんデータセットを使った二値分類 (Binary classification) から始める。

以下が XGBoost を使って乳がんデータセットを二値分類するサンプルコード。 なお、モデルの検証についてはここでの本題ではないことから、交差検証ではなくホールドアウト検証にとどめている。 最終的に検証用データで精度 (Accuracy) を確認している。 コードの説明についてはコメントを参照のこと。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import xgboost as xgb

import numpy as np
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score

"""XGBoost で二値分類するサンプルコード"""


def main():
    # 乳がんデータセットを読み込む
    dataset = datasets.load_breast_cancer()
    X, y = dataset.data, dataset.target
    # データセットを学習用と検証用に分割する
    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=0.3,
                                                        shuffle=True,
                                                        random_state=42,
                                                        stratify=y)
    # XGBoost が扱うデータセットの形式に直す
    dtrain = xgb.DMatrix(X_train, label=y_train)
    dtest = xgb.DMatrix(X_test, label=y_test)
    # 学習用のパラメータ
    xgb_params = {
        # 二値分類問題
        'objective': 'binary:logistic',
        # 評価指標
        'eval_metric': 'logloss',
    }
    # モデルを学習する
    bst = xgb.train(xgb_params,
                    dtrain,
                    num_boost_round=100,  # 学習ラウンド数は適当
                    )
    # 検証用データが各クラスに分類される確率を計算する
    y_pred_proba = bst.predict(dtest)
    # しきい値 0.5 で 0, 1 に丸める
    y_pred = np.where(y_pred_proba > 0.5, 1, 0)
    # 精度 (Accuracy) を検証する
    acc = accuracy_score(y_test, y_pred)
    print('Accuracy:', acc)


if __name__ == '__main__':
    main()

上記のサンプルコードに適当な名前をつけて保存した上で実行する。

$ python xgbhelloworld.py
[22:52:44] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 18 extra nodes, 0 pruned nodes, max_depth=5
[22:52:44] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 14 extra nodes, 0 pruned nodes, max_depth=4
[22:52:44] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 14 extra nodes, 0 pruned nodes, max_depth=4
...(snip)...
[22:52:44] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 2 extra nodes, 0 pruned nodes, max_depth=1
[22:52:44] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 2 extra nodes, 0 pruned nodes, max_depth=1
[22:52:44] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 2 extra nodes, 0 pruned nodes, max_depth=1
Accuracy: 0.9649122807017544

ホールドアウト検証の結果、精度として 0.965 前後の値が得られた。

学習過程を可視化する

先ほどの例では、いつの間にか学習が終わってモデルができたという感じだった。 そこで、続いては学習が進む過程をグラフで可視化してみる。

次のサンプルコードでは、学習用データと検証用データに対する損失を折れ線グラフで出力する。 そのためには、まず evals オプションで学習用データと検証用データを渡す。 その上で evals_result オプションに過程を記録するための辞書を渡す。 学習ラウンド数 (num_boost_round) は先ほどよりも多く 1000 まで増やした。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import xgboost as xgb

import numpy as np
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
from matplotlib import pyplot as plt

"""XGBoost で学習の履歴を可視化するサンプルコード"""


def main():
    dataset = datasets.load_breast_cancer()
    X, y = dataset.data, dataset.target

    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=0.3,
                                                        shuffle=True,
                                                        random_state=42,
                                                        stratify=y)

    dtrain = xgb.DMatrix(X_train, label=y_train)
    dtest = xgb.DMatrix(X_test, label=y_test)

    xgb_params = {
        'objective': 'binary:logistic',
        'eval_metric': 'logloss',
    }

    # 学習時に用いる検証用データ
    evals = [(dtrain, 'train'), (dtest, 'eval')]
    # 学習過程を記録するための辞書
    evals_result = {}
    bst = xgb.train(xgb_params,
                    dtrain,
                    num_boost_round=1000,  # ラウンド数を増やしておく
                    evals=evals,
                    evals_result=evals_result,
                    )

    y_pred_proba = bst.predict(dtest)
    y_pred = np.where(y_pred_proba > 0.5, 1, 0)
    acc = accuracy_score(y_test, y_pred)
    print('Accuracy:', acc)

    # 学習の課程を折れ線グラフとしてプロットする
    train_metric = evals_result['train']['logloss']
    plt.plot(train_metric, label='train logloss')
    eval_metric = evals_result['eval']['logloss']
    plt.plot(eval_metric, label='eval logloss')
    plt.grid()
    plt.legend()
    plt.xlabel('rounds')
    plt.ylabel('logloss')
    plt.show()


if __name__ == '__main__':
    main()

先ほどと同じように、適当な名前をつけて実行する。

$ python xgbhistory.py
...(snip)...
Accuracy: 0.9649122807017544

すると、次のようなグラフが得られる。

f:id:momijiame:20190129234917p:plain

学習が進むにつれて学習用データと検証用データの損失が減っていくことが分かる。 とはいえ、過学習はしていないようだけど 100 ラウンド前後で損失の減少は止まっているように見える。

損失が減らなくなったら学習を打ち切る

先ほどの例で分かる通り、損失が減らなくなったらそれ以上学習する必要はない。 それ以上回してしまうと、下手をすると損失が増える (過学習) することも考えられる。 そこで、続いては損失が減らなくなったタイミングで学習を打ち切るようにしてみる。

次のサンプルコードでは early_stopping_rounds オプションを指定することでそれを実現している。 このオプションを使うと、指定したラウンド数の間で損失が減らなかったときに学習を打ち切れる。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import xgboost as xgb

import numpy as np
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
from matplotlib import pyplot as plt

"""XGBoost で early_stopping_rounds を使って学習ラウンド数を最適化するサンプルコード"""


def main():
    dataset = datasets.load_breast_cancer()
    X, y = dataset.data, dataset.target

    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=0.3,
                                                        shuffle=True,
                                                        random_state=42,
                                                        stratify=y)

    dtrain = xgb.DMatrix(X_train, label=y_train)
    dtest = xgb.DMatrix(X_test, label=y_test)

    xgb_params = {
        'objective': 'binary:logistic',
        'eval_metric': 'logloss',
    }

    evals = [(dtrain, 'train'), (dtest, 'eval')]
    evals_result = {}
    bst = xgb.train(xgb_params,
                    dtrain,
                    num_boost_round=1000,
                    # 一定ラウンド回しても改善が見込めない場合は学習を打ち切る
                    early_stopping_rounds=10,
                    evals=evals,
                    evals_result=evals_result,
                    )

    y_pred_proba = bst.predict(dtest)
    y_pred = np.where(y_pred_proba > 0.5, 1, 0)
    acc = accuracy_score(y_test, y_pred)
    print('Accuracy:', acc)

    train_metric = evals_result['train']['logloss']
    plt.plot(train_metric, label='train logloss')
    eval_metric = evals_result['eval']['logloss']
    plt.plot(eval_metric, label='eval logloss')
    plt.grid()
    plt.legend()
    plt.xlabel('rounds')
    plt.ylabel('logloss')
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみよう。 ログ出力から、検証用データの損失に対して early stopping がかかるようになっていることが分かる。

$ python xgbearlystop.py
[23:05:43] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 18 extra nodes, 0 pruned nodes, max_depth=5
[0]    train-logloss:0.462396 eval-logloss:0.492899
Multiple eval metrics have been passed: 'eval-logloss' will be used for early stopping.

Will train until eval-logloss hasn't improved in 10 rounds.
[23:05:43] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 14 extra nodes, 0 pruned nodes, max_depth=4
...(snip)...
[23:05:44] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 2 extra nodes, 0 pruned nodes, max_depth=1
[100] train-logloss:0.006161  eval-logloss:0.09275
[23:05:44] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 2 extra nodes, 0 pruned nodes, max_depth=1
[101] train-logloss:0.00614   eval-logloss:0.093094
Stopping. Best iteration:
[91]  train-logloss:0.00637   eval-logloss:0.092265

Accuracy: 0.9649122807017544

今度は 101 ラウンドまでしか学習が進んでいない。 また、その中でも最も損失の少なかった 91 ラウンド目が結果として使われたことが分かる。

学習過程のグラフは次のようなものが得られた。

f:id:momijiame:20190129234955p:plain

たしかに 100 前後の学習ラウンドまでしか表示されていない。

scikit-learn インターフェースを使ってみる

XGBoost には、ネイティブな API の他に scikit-learn 互換の API を持ったインターフェースもある。 続いては、これを使ってみよう。

次のサンプルコードでは scikit-learn 互換の API を備えた XGBClassifier を使っている。 ただし、実現していることは先ほどと全く変わらない。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import xgboost as xgb

from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
from matplotlib import pyplot as plt

"""XGBoost の scikit-learn インターフェースを使ったサンプルコード (二値分類)"""


def main():
    dataset = datasets.load_breast_cancer()
    X, y = dataset.data, dataset.target

    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=0.3,
                                                        shuffle=True,
                                                        random_state=42,
                                                        stratify=y)

    # scikit-learn API を備えた分類器
    clf = xgb.XGBClassifier(objective='binary:logistic',
                            # 'num_boost_round' の代わり
                            # adding 1 estimator per round
                            n_estimators=1000)
    # 学習する
    evals_result = {}
    clf.fit(X_train, y_train,
            # 学習に使う評価指標
            eval_metric='logloss',
            # 学習時に用いる検証用データ
            eval_set=[
                (X_train, y_train),
                (X_test, y_test),
            ],
            early_stopping_rounds=10,
            # 学習過程の記録はコールバック API で登録する
            callbacks=[
                xgb.callback.record_evaluation(evals_result)
            ],
            )

    y_pred = clf.predict(X_test)
    acc = accuracy_score(y_test, y_pred)
    print('Accuracy:', acc)

    # 学習過程の名前は 'validation_{n}' になる
    train_metric = evals_result['validation_0']['logloss']
    plt.plot(train_metric, label='train logloss')
    eval_metric = evals_result['validation_1']['logloss']
    plt.plot(eval_metric, label='eval logloss')
    plt.grid()
    plt.legend()
    plt.xlabel('rounds')
    plt.ylabel('logloss')
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python xgbscikit.py
[0]    validation_0-logloss:0.607962  validation_1-logloss:0.615844
Multiple eval metrics have been passed: 'validation_1-logloss' will be used for early stopping.

Will train until validation_1-logloss hasn't improved in 10 rounds.
[1]   validation_0-logloss:0.538811   validation_1-logloss:0.555234
[2]   validation_0-logloss:0.480826   validation_1-logloss:0.5015
[3]   validation_0-logloss:0.430842   validation_1-logloss:0.458286
...(snip)...
[170] validation_0-logloss:0.008248   validation_1-logloss:0.094108
Stopping. Best iteration:
[160] validation_0-logloss:0.008465   validation_1-logloss:0.094088

Accuracy: 0.9649122807017544

学習ラウンド数は先ほどと違っているものの、最終的に得られた精度は変わっていないようだ。

学習過程のグラフは次のようなものが得られた。

f:id:momijiame:20190129235009p:plain

多値分類問題を扱う

続いてはタスク設定として多値分類 (Multiclass classification) を試してみよう。 この場合、二値分類とは学習するときのパラメータが異なる。

次のサンプルコードでは XGBoost で Iris データセットを使った多値分類問題を扱っている。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import xgboost as xgb

from matplotlib import pyplot as plt
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score

"""XGBoost で多値分類するサンプルコード"""


def main():
    # Iris データセットを読み込む
    dataset = datasets.load_iris()
    X, y = dataset.data, dataset.target

    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=0.3,
                                                        shuffle=True,
                                                        random_state=42,
                                                        stratify=y)

    dtrain = xgb.DMatrix(X_train, label=y_train)
    dtest = xgb.DMatrix(X_test, label=y_test)

    xgb_params = {
        # 多値分類問題
        'objective': 'multi:softmax',
        # クラス数
        'num_class': 3,
        # 学習用の指標 (Multiclass logloss)
        'eval_metric': 'mlogloss',
    }
    evals = [(dtrain, 'train'), (dtest, 'eval')]
    evals_result = {}
    bst = xgb.train(xgb_params,
                    dtrain,
                    num_boost_round=1000,
                    early_stopping_rounds=10,
                    evals=evals,
                    evals_result=evals_result,
                    )

    y_pred = bst.predict(dtest)
    acc = accuracy_score(y_test, y_pred)
    print('Accuracy:', acc)

    train_metric = evals_result['train']['mlogloss']
    plt.plot(train_metric, label='train logloss')
    eval_metric = evals_result['eval']['mlogloss']
    plt.plot(eval_metric, label='eval logloss')
    plt.grid()
    plt.legend()
    plt.xlabel('rounds')
    plt.ylabel('logloss')
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python xgbmulticlass.py
[23:14:48] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 2 extra nodes, 0 pruned nodes, max_depth=1
[23:14:48] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 8 extra nodes, 0 pruned nodes, max_depth=3
[23:14:48] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 4 extra nodes, 0 pruned nodes, max_depth=2
[0]    train-mlogloss:0.742287    eval-mlogloss:0.765776
Multiple eval metrics have been passed: 'eval-mlogloss' will be used for early stopping.

Will train until eval-mlogloss hasn't improved in 10 rounds.
...
Stopping. Best iteration:
[16]  train-mlogloss:0.039852 eval-mlogloss:0.195911

Accuracy: 0.9333333333333333

乳がんデータセットに比べると、単純な分だけラウンド数がかなり少ないようだ。

得られた学習過程のグラフは次の通り。 学習用データの損失は減っているものの、検証用データの損失が減らない状況が生じていることから過学習の予兆が見られる。

f:id:momijiame:20190129235026p:plain

特徴量の重要度を可視化する

XGBoost は決定木の仲間ということで特徴量の重要度 (Feature Importance) を可視化する機能を備えている。 次のサンプルコードでは、Iris データセットの分類にどの特徴量が有効だったのかを性能のゲインにもとづいて可視化している。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import xgboost as xgb

from sklearn import datasets
from sklearn.model_selection import train_test_split
from matplotlib import pyplot as plt

"""XGBoost で特徴量の重要度を可視化するサンプルコード"""


def main():
    dataset = datasets.load_iris()
    X, y = dataset.data, dataset.target

    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=0.3,
                                                        shuffle=True,
                                                        random_state=42,
                                                        stratify=y)

    # 可視化のために特徴量の名前を渡しておく
    dtrain = xgb.DMatrix(X_train, label=y_train,
                         feature_names=dataset.feature_names)
    dtest = xgb.DMatrix(X_test, label=y_test,
                        feature_names=dataset.feature_names)

    xgb_params = {
        'objective': 'multi:softmax',
        'num_class': 3,
        'eval_metric': 'mlogloss',
    }

    evals = [(dtrain, 'train'), (dtest, 'eval')]
    evals_result = {}
    bst = xgb.train(xgb_params,
                    dtrain,
                    num_boost_round=1000,
                    early_stopping_rounds=10,
                    evals=evals,
                    evals_result=evals_result,
                    )

    # 性能向上に寄与する度合いで重要度をプロットする
    _, ax = plt.subplots(figsize=(12, 4))
    xgb.plot_importance(bst,
                        ax=ax,
                        importance_type='gain',
                        show_values=False)
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python xgbfeatimp.py
[23:19:37] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 2 extra nodes, 0 pruned nodes, max_depth=1
[23:19:37] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 8 extra nodes, 0 pruned nodes, max_depth=3
[23:19:37] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 4 extra nodes, 0 pruned nodes, max_depth=2
[0]    train-mlogloss:0.742287    eval-mlogloss:0.765776
Multiple eval metrics have been passed: 'eval-mlogloss' will be used for early stopping.

Will train until eval-mlogloss hasn't improved in 10 rounds.
...(snip)...
Stopping. Best iteration:
[16]  train-mlogloss:0.039852 eval-mlogloss:0.195911

次のようなグラフが得られる。 どうやら petal lengthpetal width が分類する上で有効なようだ。

f:id:momijiame:20190129235040p:plain

回帰問題を扱ってみる

続いては XGBoost で回帰問題を扱ってみる。 回帰問題を扱うときは学習時のパラメータとして渡す objectivereg から始まるようになる。

次のサンプルコードでは XGBoost で Boston データセットを回帰している。 学習と検証の評価指標には RMSE (Root Mean Squared Error) を用いた。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import math

import xgboost as xgb

from matplotlib import pyplot as plt
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error

"""XGBoost で回帰するサンプルコード"""


def main():
    # Boston データセットを読み込む
    dataset = datasets.load_boston()
    X, y = dataset.data, dataset.target

    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=0.3,
                                                        shuffle=True,
                                                        random_state=42,
                                                        )

    dtrain = xgb.DMatrix(X_train, label=y_train)
    dtest = xgb.DMatrix(X_test, label=y_test)

    xgb_params = {
        # 回帰問題
        'objective': 'reg:linear',
        # 学習用の指標 (RMSE)
        'eval_metric': 'rmse',
    }
    evals = [(dtrain, 'train'), (dtest, 'eval')]
    evals_result = {}
    bst = xgb.train(xgb_params,
                    dtrain,
                    num_boost_round=1000,
                    early_stopping_rounds=10,
                    evals=evals,
                    evals_result=evals_result,
                    )

    y_pred = bst.predict(dtest)
    mse = mean_squared_error(y_test, y_pred)
    print('RMSE:', math.sqrt(mse))

    train_metric = evals_result['train']['rmse']
    plt.plot(train_metric, label='train rmse')
    eval_metric = evals_result['eval']['rmse']
    plt.plot(eval_metric, label='eval rmse')
    plt.grid()
    plt.legend()
    plt.xlabel('rounds')
    plt.ylabel('rmse')
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python xgbreg.py    
[23:25:09] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 22 extra nodes, 0 pruned nodes, max_depth=5
[0]    train-rmse:17.5096 eval-rmse:16.1546
Multiple eval metrics have been passed: 'eval-rmse' will be used for early stopping.

Will train until eval-rmse hasn't improved in 10 rounds.
...(snip)...
Stopping. Best iteration:
[33]  train-rmse:0.344815 eval-rmse:3.03055

RMSE: 3.0332711348600068

学習過程の損失をプロットしたグラフは次のようになった。

f:id:momijiame:20190129235054p:plain

カスタムメトリックを扱う

これまでは XGBoost に組み込みで入っていた評価指標を用いて学習の進み具合を評価していた。 続いては自分で書いたカスタムメトリックを使って学習してみる。

次のサンプルコードでは、再び乳がんデータセットを使った二値分類を扱っている。 ただし、学習を評価するメトリックとして精度 (Accuracy) を計測するカスタムメトリックを使っている。 カスタムメトリックを扱うにはオプションとして feval に自分で書いた関数を渡す。 カスタムメトリックの関数は、評価指標の名前と値をタプルまたはリストで返すように作る。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import xgboost as xgb

import numpy as np
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
from matplotlib import pyplot as plt

"""XGBoost で学習時にカスタムメトリックを使ったサンプルコード"""


def feval_accuracy(pred_proba, dtrain):
    """カスタムメトリックを計算する関数"""
    # 真のデータ
    y_true = dtrain.get_label().astype(int)
    # 予測
    y_pred = np.where(pred_proba > 0.5, 1, 0)
    # Accuracy を計算する
    acc = accuracy_score(y_true, y_pred)
    # メトリックの名前と数値を返す(最小化を目指すので 1 から引く)
    return 'accuracy', 1 - acc


def main():
    dataset = datasets.load_breast_cancer()
    X, y = dataset.data, dataset.target

    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=0.3,
                                                        shuffle=True,
                                                        random_state=42,
                                                        stratify=y)

    dtrain = xgb.DMatrix(X_train, label=y_train)
    dtest = xgb.DMatrix(X_test, label=y_test)

    xgb_params = {
        'objective': 'binary:logistic',
        'eval_metric': 'logloss',
    }

    evals = [(dtrain, 'train'), (dtest, 'eval')]
    evals_result = {}
    bst = xgb.train(xgb_params,
                    dtrain,
                    num_boost_round=1000,
                    early_stopping_rounds=10,
                    evals=evals,
                    evals_result=evals_result,
                    # カスタマイズした評価関数を使う
                    feval=feval_accuracy,
                    )

    y_pred_proba = bst.predict(dtest)
    y_pred = np.where(y_pred_proba > 0.5, 1, 0)
    acc = accuracy_score(y_test, y_pred)
    print('Accuracy:', acc)

    # デフォルトのメトリック
    _, ax1 = plt.subplots(figsize=(8, 4))
    train_metric = evals_result['train']['logloss']
    ax1.plot(train_metric, label='train logloss', c='r')
    eval_metric = evals_result['eval']['logloss']
    ax1.plot(eval_metric, label='eval logloss', c='g')
    ax1.set_ylabel('logloss')
    ax1.legend()
    ax1.set_xlabel('rounds')

    # カスタムメトリック
    ax2 = ax1.twinx()
    eval_custom_metric = evals_result['eval']['accuracy']
    ax2.plot(eval_custom_metric, label='eval accuracy', c='b')
    ax2.set_ylabel('accuracy')
    ax2.legend()

    plt.grid()
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみる。

$ python xgbcustmetric.py
[23:35:00] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 18 extra nodes, 0 pruned nodes, max_depth=5
[0]    train-logloss:0.462396 eval-logloss:0.492899  train-accuracy:0.017588    eval-accuracy:0.070175
Multiple eval metrics have been passed: 'eval-accuracy' will be used for early stopping.

Will train until eval-accuracy hasn't improved in 10 rounds.
...(snip)...
Stopping. Best iteration:
[14]  train-logloss:0.029256  eval-logloss:0.112603   train-accuracy:0.005025 eval-accuracy:0.046784

Accuracy: 0.9532163742690059

得られた学習過程のグラフは次の通り。

f:id:momijiame:20190129235205p:plain

学習の評価指標として精度 (Accuracy) を使うと、早々に学習が打ち切られてしまっているようだ。 そのため、最終的に得られた精度も低いものになってしまっている。 これは、おそらく一つの要素が正しく分類できたか・できなかったに結果が大きく左右されてしまうため。 まあ、結果はそんなに良くないけど、こんな感じで書けるよという例として。

組み込みの交差検証の機能を使ってみる

XGBoost には、組み込みの交差検証 (Cross Validation) の機能がある。 正直、そんなに使いやすいものではないけど、とりあえずあるよということで紹介しておく。

次のサンプルコードでは、XGBoost に組み込みで入っている交差検証の機能を使っている。 機能は xgboost.cv() という関数が起点になる。 基本的なパラメータについては xgboost.train() で普通に学習するときとさほど変わらない。 関数から最終的に得られるのは、学習過程におけるメトリックの値の変化。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import xgboost as xgb

import numpy as np
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
from matplotlib import pyplot as plt

"""XGBoost 組み込みの交差検証を使ったサンプルコード"""

def main():
    dataset = datasets.load_breast_cancer()
    X, y = dataset.data, dataset.target

    # CV の中で分割するので丸ごと渡す
    dtrain = xgb.DMatrix(X, label=y)

    xgb_params = {
        'objective': 'binary:logistic',
        'eval_metric': 'logloss',
    }

    # 交差検証する
    history = xgb.cv(xgb_params,
                     dtrain,
                     num_boost_round=1000,
                     early_stopping_rounds=10,
                     nfold=10,
                     # 層化分割する
                     stratified=True,
                     # 検証の経過を出力する
                     verbose_eval=True,
                     )

    train_metric = history['train-logloss-mean']
    plt.plot(train_metric, label='train logloss')
    eval_metric = history['test-logloss-mean']
    plt.plot(eval_metric, label='eval logloss')
    plt.grid()
    plt.legend()
    plt.xlabel('rounds')
    plt.ylabel('logloss')
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみる。

$ python xgbcustmetric.py
[23:35:00] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 18 extra nodes, 0 pruned nodes, max_depth=5
[0]    train-logloss:0.462396 eval-logloss:0.492899  train-accuracy:0.017588    eval-accuracy:0.070175
Multiple eval metrics have been passed: 'eval-accuracy' will be used for early stopping.

Will train until eval-accuracy hasn't improved in 10 rounds.
...(snip)...
[23:39:47] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 2 extra nodes, 0 pruned nodes, max_depth=1

すると、次のような学習過程のグラフが得られる。

f:id:momijiame:20190129235216p:plain

いじょう。

まとめ

今回は XGBoost の機能を色々と試してみた。

統計的学習の基礎 ―データマイニング・推論・予測―

統計的学習の基礎 ―データマイニング・推論・予測―

  • 作者: Trevor Hastie,Robert Tibshirani,Jerome Friedman,杉山将,井手剛,神嶌敏弘,栗田多喜夫,前田英作,井尻善久,岩田具治,金森敬文,兼村厚範,烏山昌幸,河原吉伸,木村昭悟,小西嘉典,酒井智弥,鈴木大慈,竹内一郎,玉木徹,出口大輔,冨岡亮太,波部斉,前田新一,持橋大地,山田誠
  • 出版社/メーカー: 共立出版
  • 発売日: 2014/06/25
  • メディア: 単行本
  • この商品を含むブログ (6件) を見る

Python: Hyperopt で機械学習モデルのハイパーパラメータを選ぶ

今回は、機械学習モデルのハイパーパラメータをチューニングするのに用いられる Python のフレームワークの一つとして Hyperopt を使ってみる。 このフレームワークは、機械学習コンペティションの一つである Kaggle でよく用いられるものとして知られている。

なお、このブログでは過去にハイパーパラメータのチューニングについて Bayesian OptimizationOptuna を使った例を扱ったことがある。

blog.amedama.jp

blog.amedama.jp

今回使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.14.2
BuildVersion:   18C54
$ python -V
Python 3.7.2

下準備

まずは下準備として、必要になるパッケージをインストールする。

$ pip install hyperopt matplotlib scikit-learn pandas

単純な例で動作を確認する

機械学習モデルのハイパーパラメータを扱う前に、まずはもっと単純な例から始める。 手始めに、次のような一つの変数 x を取る関数の最大値を探す問題について考えてみよう。 グラフを見ると、最大値は 2 付近にあるようだ。

f:id:momijiame:20180818132951p:plain

上記の関数の最大値を Hyperopt で探すサンプルコードは次の通り。 Hyperopt では、最適化したい目的関数を用意して、それを hyperopt.fmin() に渡す。 以下のサンプルコードにおいて目的関数は objective() という名前で定義している。 チューニングするパラメータは、どんな分布でどういった値域を持つのかあらかじめ定義しておく。 以下では hp.uniform('x', -5, +15) がそれで、値域が -5 ~ +15 の一様分布を定義している。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import math

from hyperopt import hp
from hyperopt import fmin
from hyperopt import tpe


def f(x):
    """最大値を見つけたい関数"""
    return math.exp(-(x - 2) ** 2) + math.exp(-(x - 6) ** 2 / 10) + 1 / (x ** 2 + 1)


def objective(args):
    """最小化したい目的関数"""
    return -1 * f(args)  # 元のタスクが最大化なので符号を反転する


def main():
    # 変数の値域を定義する
    space = hp.uniform('x', -5, +15)
    # 目的関数を最小化するパラメータを探索する
    best = fmin(fn=objective, space=space, algo=tpe.suggest, max_evals=100)
    # 結果を出力する
    print(best)


if __name__ == '__main__':
    main()

上記に適当な名前をつけて実行する。 すると、たしかに解析的に確認していた最大値として 2 前後の値が得られる。

$ python hellohp.py
{'x': 2.0069551331288276}

探索過程を記録する

先ほどの例では、最終的に最適化された値が得られておしまいだった。 とはいえ、実際に使う場面では、その値がどういった経緯を経て得られたのか重要になる場合もある。 そこで、Hyperopt では Trials というオブジェクトに過程を記録できる。

次のサンプルコードでは Trials のインスタンスに探索過程を記録した上で、最後にそれを出力している。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import math

from hyperopt import hp
from hyperopt import fmin
from hyperopt import tpe
from hyperopt import Trials
from hyperopt import STATUS_OK


def f(x):
    """最大値を見つけたい関数"""
    return math.exp(-(x - 2) ** 2) + math.exp(-(x - 6) ** 2 / 10) + 1 / (x ** 2 + 1)


def objective(args):
    """最小化したい目的関数"""
    # 計算過程を追跡するときは目的関数の返り値を辞書にする
    return {
        'loss': -1 * f(args),  # 最小化する損失 (必須パラメータ)
        'status': STATUS_OK,  # 計算結果 (必須パラメータ)
    }


def main():
    # 変数の値域を定義する
    space = hp.uniform('x', -5, +15)
    # 探索過程を記録するオブジェクト
    trials = Trials()
    # 目的関数を最小化するパラメータを探索する
    fmin(fn=objective, space=space, algo=tpe.suggest, max_evals=100, trials=trials)
    # 探索過程を出力する
    for item in trials.trials:
        vals = item['misc']['vals']
        result = item['result']
        print('vals:', vals, 'result:', result)


if __name__ == '__main__':
    main()

実行すると、どのような値を探索して結果がどうだったのかが確認できる。

$ python trials.py
vals: {'x': [10.496518222653432]} result: {'loss': -0.14140262193324998, 'status': 'ok'}
vals: {'x': [13.094879827038248]} result: {'loss': -0.012312365591358516, 'status': 'ok'}
vals: {'x': [-2.0228004674552014]} result: {'loss': -0.19799926655007422, 'status': 'ok'}
...(snip)...
vals: {'x': [2.6501743089141496]} result: {'loss': -1.1054773169025645, 'status': 'ok'}
vals: {'x': [9.247128914308604]} result: {'loss': -0.35996620277329294, 'status': 'ok'}
vals: {'x': [7.18849940727293]} result: {'loss': -0.8872540482414317, 'status': 'ok'}

探索過程に色々な値を記録する

先ほど使った Trials オブジェクトには、実は自分で色々な値を記録することもできる。 次のサンプルコードでは、探索した時刻と関数の引数を記録している。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import time
import math

from hyperopt import hp
from hyperopt import fmin
from hyperopt import tpe
from hyperopt import Trials
from hyperopt import STATUS_OK


def f(x):
    """最大値を見つけたい関数"""
    return math.exp(-(x - 2) ** 2) + math.exp(-(x - 6) ** 2 / 10) + 1 / (x ** 2 + 1)


def objective(args):
    """最小化したい目的関数"""
    return {
        'loss': -1 * f(args),  # 最小化する損失 (必須パラメータ)
        'status': STATUS_OK,  # 計算結果 (必須パラメータ)
        # 一定の制約の元に自由に入れて良い
        'estimate_time': time.time(),
        'vals': args,
    }


def main():
    # 変数の値域を定義する
    space = hp.uniform('x', -5, +15)
    # 探索過程を記録するオブジェクト
    trials = Trials()
    # 目的関数を最小化するパラメータを探索する
    fmin(fn=objective, space=space, algo=tpe.suggest, max_evals=100, trials=trials)
    # 探索過程を出力する
    for item in trials.trials:
        result = item['result']
        print('result:', result)


if __name__ == '__main__':
    main()

上記を実行すると、探索した時刻や引数が得られていることが分かる。

$ python params.py 
result: {'loss': -0.339479707912701, 'status': 'ok', 'estimate_time': 1548600261.9176042, 'vals': -1.407873925928973}
result: {'loss': -0.682422018336043, 'status': 'ok', 'estimate_time': 1548600261.919095, 'vals': 8.012065634262315}
result: {'loss': -0.055263914137346264, 'status': 'ok', 'estimate_time': 1548600261.920187, 'vals': -4.135976583356424}
...(snip)...
result: {'loss': -1.1031356060488282, 'status': 'ok', 'estimate_time': 1548600262.10813, 'vals': 2.65344333175562}
result: {'loss': -0.27738481344650867, 'status': 'ok', 'estimate_time': 1548600262.111007, 'vals': -1.6261239199291966}
result: {'loss': -0.14645158143456388, 'status': 'ok', 'estimate_time': 1548600262.1134, 'vals': -2.4222183825698775}

探索過程を可視化してみる

先ほどの例では、どんな風に探索が進められているのかがいまいちイメージしにくかった。 そこで、どんな風に探索しているのかをグラフ上にプロットしてみた。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import time
import math

from matplotlib import pyplot as plt

import numpy as np

from hyperopt import hp
from hyperopt import fmin
from hyperopt import tpe
from hyperopt import Trials
from hyperopt import STATUS_OK


def f(x):
    """最大値を見つけたい関数"""
    return math.exp(-(x - 2) ** 2) + math.exp(-(x - 6) ** 2 / 10) + 1 / (x ** 2 + 1)


def objective(args):
    """最小化したい目的関数"""
    return {
        'loss': -1 * f(args),  # 最小化する損失 (必須パラメータ)
        'status': STATUS_OK,  # 計算結果 (必須パラメータ)
        'vals': args,
    }


def main():
    # 変数の値域を定義する
    space = hp.uniform('x', -5, +15)
    # 探索過程を記録するオブジェクト
    trials = Trials()
    # 目的関数を最小化するパラメータを探索する
    fmin(fn=objective, space=space, algo=tpe.suggest, max_evals=100, trials=trials)

    # 検索過程をグラフで可視化してみる
    args = [result['vals'] for result in trials.results]
    values = [-result['loss'] for result in trials.results]
    # 元のグラフをプロットする
    X = [x for x in np.arange(-5, 15, 0.1)]
    y = [f(x) for x in X]
    plt.plot(X, y)
    # 探索した場所をプロットする
    plt.scatter(args, values)
    # グラフを描画する
    plt.xlabel('x')
    plt.ylabel('y')
    plt.grid()
    plt.show()


if __name__ == '__main__':
    main()

上記に適当な名前をつけて実行する。

$ python visualize.py

すると、次のようなグラフが得られる。

f:id:momijiame:20190128002125p:plain

頂点付近以外にも、かなり色んな所を調べていることが伺える。

機械学習のモデルに適用してみる

続いては、ついに Hyperopt を機械学習のモデルに適用してみる。 モデルのアルゴリズムにはサポートベクターマシンを使った。 パラメータとしては kernelCgamma をチューニングする。

次のサンプルコードでは、目的関数でサポートベクターマシンを学習させて 5-Fold CV で Accuracy を評価している。 kernel['linear', 'rbf', 'poly'] の三種類から選ぶ。 Cgamma は常用対数で 0 ~ 100-10 ~ 1 の範囲で探索する。 2.303 を乗じているのは二進対数と常用対数に変換するため。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from functools import partial

from hyperopt import hp
from hyperopt import fmin
from hyperopt import tpe
from hyperopt import Trials
from hyperopt import space_eval

from sklearn import datasets
from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import cross_validate

from sklearn.svm import SVC


def objective(X, y, args):
    """最小化したい目的関数"""
    # モデルのアルゴリズムに SVM を使う
    model = SVC(**args)
    # Stratified 5 Fold Cross Validation
    kf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
    scores = cross_validate(model, X=X, y=y, cv=kf)
    # 最小化なので符号を反転する
    return -1 * scores['test_score'].mean()


def main():
    # データセットを読み込む
    dataset = datasets.load_iris()
    X, y = dataset.data, dataset.target
    # 目的関数にデータセットを部分適用する
    f = partial(objective, X, y)
    # 変数の値域を定義する
    space = {
        # 底が自然対数なので 2.303 をかけて常用対数に変換する
        'C': hp.loguniform('C', 2.303 * 0, 2.303 * +2),
        'gamma': hp.loguniform('gamma', 2.303 * -2, 2.303 * +1),
        'kernel': hp.choice('kernel', ['linear', 'rbf', 'poly']),
    }
    # 探索過程を記録するオブジェクト
    trials = Trials()
    # 目的関数を最小化するパラメータを探索する
    best = fmin(fn=f, space=space, algo=tpe.suggest, max_evals=100, trials=trials)
    # 結果を出力する
    print(space_eval(space, best))


if __name__ == '__main__':
    main()

実行すると、次のような結果が得られる。 値は似通ったものになるはずだけど、結果は毎回異なるはず。

$ python svm.py 
{'C': 1.9877559946331484, 'gamma': 0.7522911704230325, 'kernel': 'linear'}

探索を可視化してみる

さきほどと同じように、どのような値を探索したのかグラフにプロットしてみよう。 次のサンプルコードでは、カーネルを linear に決め打ちした状態で、Hyperopt が探索した Cgamma を二次元でプロットしている。 標準化した Accuracy の高低は点の色で表現している。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from functools import partial

from hyperopt import hp
from hyperopt import fmin
from hyperopt import tpe
from hyperopt import Trials
from hyperopt import STATUS_OK
from hyperopt import space_eval

from sklearn import datasets
from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import cross_validate
from sklearn.svm import SVC

import numpy as np

import pandas as pd

from matplotlib import pyplot as plt
from matplotlib import cm


def objective(X, y, args):
    """最小化したい目的関数"""
    # モデルのアルゴリズムに Linear SVM を使う
    model = SVC(kernel='linear', **args)
    # Stratified 5 Fold Cross Validation
    kf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
    scores = cross_validate(model, X=X, y=y, cv=kf)
    # 最小化なので符号を反転する
    return {
        'loss': -1 * scores['test_score'].mean(),
        'status': STATUS_OK,
        'vals': args,
    }


def main():
    # データセットを読み込む
    dataset = datasets.load_iris()
    X, y = dataset.data, dataset.target
    # 目的関数にデータセットを部分適用する
    f = partial(objective, X, y)
    # 変数の値域を定義する
    space = {
        'C': hp.loguniform('C', 2.303 * 0, 2.303 * +2),
        'gamma': hp.loguniform('gamma', 2.303 * -2, 2.303 * +1),
    }
    # 探索過程を記録するオブジェクト
    trials = Trials()
    # 目的関数を最小化するパラメータを探索する
    best = fmin(fn=f, space=space, algo=tpe.suggest, max_evals=100, trials=trials)
    # 結果を出力する
    print('best parameters:', space_eval(space, best))

    # 探索過程を可視化する
    xs = [result['vals']['C'] for result in trials.results]
    ys = [result['vals']['gamma'] for result in trials.results]
    zs = np.array([-1 * result['loss'] for result in trials.results])
    s_zs = (zs - zs.min()) / (zs.max() - zs.min()) # 0 ~ 1 の範囲に正規化する
    # プロットする
    sc = plt.scatter(xs, ys, c=s_zs, s=20, zorder=10, cmap=cm.cool)
    plt.colorbar(sc)

    plt.xlabel('C')
    plt.xscale('log')
    plt.ylabel('gamma')
    plt.yscale('log')
    plt.grid()
    plt.show()

    print('-- accuracy summary --')
    print(pd.Series(np.array(zs).ravel()).describe())


if __name__ == '__main__':
    main()

上記を実行する。

$ python svmplot.py
best parameters: {'C': 2.3141971486483675, 'gamma': 0.014160000992695598}
-- accuracy summary --
count    100.000000
mean       0.975800
std        0.008350
min        0.960000
25%        0.966667
50%        0.980000
75%        0.980000
max        0.986667
dtype: float64

すると、次のようなグラフが得られる。

f:id:momijiame:20190128003937p:plain

たしかに、場所によって性能に違いが出ているようだ。 見た感じ gamma よりも C の方が性能に違いを与えていそう。

複数のモデルから選ぶ

先ほどの例では、サポートベクターマシンの中でハイパーパラメータをチューニングしていた。 これだと、別のモデルを含めて評価したいときは別々に実行する必要がある。 そんなとき Hyperopt ではまとめて評価することもできる。

次のサンプルコードではランダムフォレスト、サポートベクターマシン、ロジスティック回帰を一度に評価している。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from functools import partial

from hyperopt import hp
from hyperopt import fmin
from hyperopt import tpe
from hyperopt import Trials
from hyperopt import space_eval
from hyperopt.pyll.base import scope

from sklearn import datasets
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import cross_validate

from sklearn.svm import SVC


def objective(X, y, args):
    """最小化したい目的関数"""
    classifiers = {
        'svm': SVC,
        'rf': RandomForestClassifier,
        'logit': LogisticRegression,
    }
    classifier = classifiers.get(args['model_type'])
    del args['model_type']
    model = classifier(**args)
    # Stratified 5 Fold Cross Validation
    kf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
    scores = cross_validate(model, X=X, y=y, cv=kf)
    # 最小化なので符号を反転する
    return -1 * scores['test_score'].mean()


def main():
    # データセットを読み込む
    dataset = datasets.load_iris()
    X, y = dataset.data, dataset.target
    # 目的関数にデータセットを部分適用する
    f = partial(objective, X, y)
    # 変数の値域を定義する
    space = hp.choice('algorithms', [
        {
            'model_type': 'rf',
            'n_estimators': scope.int(hp.uniform('n_estimators', 1e+1, 1e+3)),
            'max_depth': scope.int(hp.uniform('max_depth', 1e+1, 1e+3)),
        },
        {
            'model_type': 'svm',
            'C': hp.uniform('C', 1e+0, 1e+2),
            'gamma': hp.lognormal('gamma', 1e-2, 1e+1),
        },
        {
            'model_type': 'logit',
            'solver': 'lbfgs',
            'multi_class': 'auto',
            'max_iter': 1000,
        }
    ])
    # 探索過程を記録するオブジェクト
    trials = Trials()
    # 目的関数を最小化するパラメータを探索する
    best = fmin(fn=f, space=space, algo=tpe.suggest, max_evals=100, trials=trials)
    # 結果を出力する
    print(space_eval(space, best))


if __name__ == '__main__':
    main()

上記を実行してみよう。 このパラメータの中ではサポートベクターマシンが最も良い結果が得られたようだ。

$ python multimodel.py
{'C': 98.52420177369001, 'gamma': 0.011789334572006299, 'model_type': 'svm'}

いじょう。

Ubuntu 18.04 LTS でコマンドのソースコードを取得する

Ubuntu を使っていて、このコマンドってどんな処理してるんだろうなーとか気になったときにソースコードを取り寄せる方法について。

使った環境は次の通り。

$ cat /etc/lsb-release 
DISTRIB_ID=Ubuntu
DISTRIB_RELEASE=18.04
DISTRIB_CODENAME=bionic
DISTRIB_DESCRIPTION="Ubuntu 18.04.1 LTS"
$ uname -r
4.15.0-29-generic

題材は ls コマンドにする。

$ which ls
/bin/ls

まずは下準備として APT のソースコード用リポジトリを有効にしておく。

$ sudo sed -i.back -e "s/^# deb-src/deb-src/g" /etc/apt/sources.list
$ sudo apt-get update

続いて、お目当てのコマンドの実行ファイルが含まれるパッケージ名を dpkg -S で調べる。 今回は coreutils にあることが分かった。

$ dpkg -S $(which ls)
coreutils: /bin/ls

あとは apt-get source コマンドでパッケージのソースコードを取得する。

$ apt-get source coreutils

これでカレントディレクトリにソースコードがダウンロードされる。

$ ls
coreutils-8.28                         coreutils_8.28-1ubuntu1.dsc  coreutils_8.28.orig.tar.xz.asc
coreutils_8.28-1ubuntu1.debian.tar.xz  coreutils_8.28.orig.tar.xz
$ ls coreutils-8.28/src/ | head
base64.c
basename.c
blake2
cat.c
chcon.c
chgrp.c
chmod.c
chown.c
chown-core.c
chown-core.h

展開されたファイルの中に ls.c という、それらしき名前のファイルが見つかった。

$ ls coreutils-8.28/src/ | grep ^ls.c$
ls.c

ビルドしてみる

ダウンロードしたのが本当に ls コマンドのソースコードなのか、ビルドすることで確かめてみよう。

$ sudo apt-get -y install build-essential

パッケージのソースコードが入ったディレクトリに移動する。

$ cd coreutils-8.28

ビルドする。

$ ./configure
$ make

すると Linux の実行ファイルができる。

$ file src/ls
src/ls: ELF 64-bit LSB shared object, x86-64, version 1 (SYSV), dynamically linked, interpreter /lib64/ld-linux-x86-64.so.2, for GNU/Linux 3.2.0, BuildID[sha1]=b29ae91cb81dca67ce217d928be9a2926b3008c2, with debug_info, not stripped

試しに実行してみよう。

$ ./src/ls
ABOUT-NLS   build-aux      configure      doc       lib      Makefile.in  src   THANKS-to-translators
aclocal.m4  cfg.mk         configure.ac   gnulib-tests  m4       man          tests   THANKStt.in
AUTHORS     ChangeLog      COPYING        GNUmakefile   maint.mk     NEWS         THANKS      TODO
bootstrap   config.log     debian       init.cfg      Makefile     po       thanks-gen
bootstrap.conf  config.status  dist-check.mk  INSTALL       Makefile.am  README       THANKS.in

ちゃんと ls のようだ。

めでたしめでたし。

macOS で Raspberry Pi 用のシリアルコンソールケーブル (PL2303) を使う

Raspberry Pi はシリアル接続のルートを確保しておくと、操作するのにディスプレイやキーボードを用意する必要がない。 最終的にネットワーク越しに SSH などで操作するにしても、初期のセットアップで楽ができる。

今回使ったシリアルコンソールケーブルはこちら。 このケーブルには PL2303 というチップが使われている。

使った環境は以下の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.14.2
BuildVersion:   18C54

ドライバをインストールする

Homebrew でチップ名を検索するとドライバが見つかる。

$ brew search PL2303
==> Casks
homebrew/cask-drivers/prolific-pl2303

インストールする。

$ brew install homebrew/cask-drivers/prolific-pl2303

シリアルコンソールケーブルを Mac につなぐ

あとは Mac にシリアルコンソールケーブルをつなぐとデバイスを認識する。

$ ls /dev/ | grep tty.usb
tty.usbserial

もし、上記だけでは認識しない場合は macOS のセキュリティ機能でドライバがブロックされている恐れがある。 「システム環境設定 > セキュリティとプライバシー > 一般」で確認しよう。

ケーブルの端子を GPIO につなぐ

つづいてはケーブルの端子を Raspberry Pi の GPIO につなぐ。

公式サイトの記述によると GPIO14 と GPIO15 が、それぞれ TX と RX に対応している。

www.raspberrypi.org

ケーブルの端子は緑が TX で白が RX に対応している。 上記の通り GPIO14 と GPIO15 につなごう。 黒は GND なので、適当な GND ピンにつないでおく。

f:id:momijiame:20190120233917j:plain:w320

赤は 5V 出力なので使わない。 下手に GND なんかにつなぐとショートする恐れがあるので注意しよう。

あとは screen コマンドを使ってデバイスに接続するだけ。 ボーレートはデフォルトで 115,200 bps に設定されている。

$ screen /dev/tty.usbserial 115200
...(snip)...
Raspbian GNU/Linux 9 raspberrypi ttyAMA0
raspberrypi login:

ちなみに、ドライバがいまいち不安定なので取り扱いに注意を要する。 たまに暴走して OS がシャットダウンできない事態に陥ったり。 使い終わったら macOS ごと終了させてケーブルを抜く、というオペレーションをするのが安全そうかな。