CUBE SUGAR CONTAINER

技術系のこと書きます。

Python: matplotlib で動的にグラフを生成する

今回は matplotlib を使って動的にグラフを生成する方法について。 ここでいう動的というのは、データを逐次的に作って、それを随時グラフに反映していくという意味を指す。 例えば機械学習のモデルを学習させるときに、その過程 (損失の減り方とか) を眺める用途で便利だと思う。

使った環境は次の通り。

$ sw_vers 
ProductName:    Mac OS X
ProductVersion: 10.13.5
BuildVersion:   17F77
$ python -V
Python 3.6.5
$ pip list --format=columns | egrep -i "(matplotlib|pillow)"
matplotlib      2.2.2  
Pillow          5.2.0  

下準備

まずは今回使うパッケージをインストールしておく。

$ pip install matplotlib pillow

静的にグラフを生成する

動的な生成について説明する前に、まずは静的なグラフの生成から説明する。 といっても、これは一般的な matplotlib のグラフの作り方そのもの。 あらかじめ必要なデータを全て用意しておいて、それをグラフとしてプロットする。

この場合、当たり前だけどプロットする前に全てのデータが揃っていないといけない。 例えば機械学習なら、モデルの学習を終えて各エポックなりラウンドごとの損失が出揃っている状態まで待つ必要がある。

次のサンプルコードではサイン波のデータをあらかじめ作った上で、それを折れ線グラフにしている。

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

import math

import numpy as np
from matplotlib import pyplot as plt


def main():
    # 描画領域
    fig = plt.figure(figsize=(10, 6))
    # 描画するデータ
    x = np.arange(0, 10, 0.1)
    y = [math.sin(i) for i in x]

    # グラフを描画する
    plt.plot(x, y) 

    # グラフを表示する
    plt.show()


if __name__ == '__main__':
    main()

上記を適当な名前でファイルに保存して実行してみよう。

$ python sin.py

すると、次のようなグラフが表示される。

f:id:momijiame:20180712235912p:plain

これが静的なグラフ生成の場合。

動的にグラフを生成する

続いて動的にグラフを生成する方法について。 これには matplotlib.animation パッケージを使う。 特に FuncAnimation を使うと作りやすい。

animation — Matplotlib 2.2.2 documentation

次のサンプルコードでは、先ほどの例と同じサイン波を動的に生成している。 ポイントは、グラフの再描画を担当する関数をコールバックとして FuncAnimation に登録すること。 そうすれば、あとは FuncAnimation が一定間隔でその関数を呼び出してくれる。 呼び出されるコールバック関数の中でデータを生成したりグラフを再描画する。

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

import math

import numpy as np
from matplotlib import pyplot as plt
from matplotlib import animation


def _update(frame, x, y):
    """グラフを更新するための関数"""
    # 現在のグラフを消去する
    plt.cla()
    # データを更新 (追加) する
    x.append(frame)
    y.append(math.sin(frame))
    # 折れ線グラフを再描画する
    plt.plot(x, y)


def main():
    # 描画領域
    fig = plt.figure(figsize=(10, 6))
    # 描画するデータ (最初は空っぽ)
    x = []
    y = []

    params = {
        'fig': fig,
        'func': _update,  # グラフを更新する関数
        'fargs': (x, y),  # 関数の引数 (フレーム番号を除く)
        'interval': 10,  # 更新間隔 (ミリ秒)
        'frames': np.arange(0, 10, 0.1),  # フレーム番号を生成するイテレータ
        'repeat': False,  # 繰り返さない
    }
    anime = animation.FuncAnimation(**params)

    # グラフを表示する
    plt.show()


if __name__ == '__main__':
    main()

先ほどと同じようにファイルに保存したら実行する。

$ python sin.py

すると、次のように動的にグラフが描画される。

f:id:momijiame:20180712235931g:plain

ちなみに、上記のような GIF 画像や動画は次のようにすると保存できる。

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

import math

import numpy as np
from matplotlib import pyplot as plt
from matplotlib import animation


def _update(frame, x, y):
    """グラフを更新するための関数"""
    # 現在のグラフを消去する
    plt.cla()
    # データを更新 (追加) する
    x.append(frame)
    y.append(math.sin(frame))
    # 折れ線グラフを再描画する
    plt.plot(x, y)


def main():
    # 描画領域
    fig = plt.figure(figsize=(10, 6))
    # 描画するデータ (最初は空っぽ)
    x = []
    y = []

    params = {
        'fig': fig,
        'func': _update,  # グラフを更新する関数
        'fargs': (x, y),  # 関数の引数 (フレーム番号を除く)
        'interval': 10,  # 更新間隔 (ミリ秒)
        'frames': np.arange(0, 10, 0.1),  # フレーム番号を生成するイテレータ
        'repeat': False,  # 繰り返さない
    }
    anime = animation.FuncAnimation(**params)

    # グラフを保存する
    anime.save('sin.gif', writer='pillow')


if __name__ == '__main__':
    main()

グラフを延々と描画し続ける

先ほどの例では frames オプションに渡すイテレータに終わりがあった。 具体的には 0 ~ 10 の範囲を 0.1 区切りで分割した 100 のデータに対してグラフを生成した。 また repeat オプションに False を指定することで繰り返し描画することも抑制している。

続いては、先ほどとは異なり frames オプションに終わりのないイテレータを渡してみよう。 こうすると、手動で止めるかメモリなどのリソースを食いつぶすまでは延々とデータを生成してグラフを描画することになる。

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

import itertools
import math

import numpy as np
from matplotlib import pyplot as plt
from matplotlib import animation


def _update(frame, x, y):
    """グラフを更新するための関数"""
    # 現在のグラフを消去する
    plt.cla()
    # データを更新 (追加) する
    x.append(frame)
    y.append(math.sin(frame))
    # 折れ線グラフを再描画する
    plt.plot(x, y)


def main():
    # 描画領域
    fig = plt.figure(figsize=(10, 6))
    # 描画するデータ
    x = []
    y = []

    params = {
        'fig': fig,
        'func': _update,  # グラフを更新する関数
        'fargs': (x, y),  # 関数の引数 (フレーム番号を除く)
        'interval': 10,  # 更新間隔 (ミリ秒)
        'frames': itertools.count(0, 0.1),  # フレーム番号を無限に生成するイテレータ
    }
    anime = animation.FuncAnimation(**params)

    # グラフを表示する
    plt.show()


if __name__ == '__main__':
    main()

上記をファイルに保存して実行してみよう。

$ python sin.py

ずーーーっとグラフが生成され続けるはず。

f:id:momijiame:20180715153531p:plain

Jupyter Notebook 上で動的にグラフを生成する

Jupyter Notebook 上で動的なグラフ生成をするときは、次のように %matplotlib nbagg マジックコマンドを使う。 また、意図的に pyplot.show() を呼び出す必要はない。

%matplotlib nbagg

import itertools
import math

import numpy as np
from matplotlib import pyplot as plt
from matplotlib import animation


def _update(frame, x, y):
    """グラフを更新するための関数"""
    # 現在のグラフを消去する
    plt.cla()
    # データを更新 (追加) する
    x.append(frame)
    y.append(math.sin(frame))
    # 折れ線グラフを再描画する
    plt.plot(x, y)


# 描画領域
fig = plt.figure(figsize=(10, 6))
# 描画するデータ
x = []
y = []

params = {
    'fig': fig,
    'func': _update,  # グラフを更新する関数
    'fargs': (x, y),  # 関数の引数 (フレーム番号を除く)
    'interval': 10,  # 更新間隔 (ミリ秒)
    'frames': np.arange(0, 10, 0.1),  # フレーム番号を生成するイテレータ
    'repeat': False,  # 繰り返さない
}
anime = animation.FuncAnimation(**params)

データの更新間隔とグラフの再描画間隔をずらす

グラフの再描画はそこまで軽い処理でもないし、データの更新間隔とずらしたいときもあるかも。 そんなときはデータを更新するスレッドと、グラフを再描画するスレッドを分ける。

以下のサンプルコードではデータの更新用に新しくスレッドを起動している。 データの更新間隔が 100ms 間隔なのに対してグラフの再描画は 250ms 間隔にしている。

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

import time
import threading
import math
import itertools

from matplotlib import pyplot as plt
from matplotlib import animation


def _redraw(_, x, y):
    """グラフを再描画するための関数"""
    # 現在のグラフを消去する
    plt.cla()
    # 折れ線グラフを再描画する
    plt.plot(x, y)


def main():
    # 描画領域
    fig = plt.figure(figsize=(10, 6))
    # 描画するデータ (最初は空っぽ)
    x = []
    y = []

    def _update():
        """データを一定間隔で追加するスレッドの処理"""
        for frame in itertools.count(0, 0.1):
            x.append(frame)
            y.append(math.sin(frame))
            # データを追加する間隔 (100ms)
            time.sleep(0.1)

    def _init():
        """データを一定間隔で追加するためのスレッドを起動する"""
        t = threading.Thread(target=_update)
        t.daemon = True
        t.start()

    params = {
        'fig': fig,
        'func': _redraw,  # グラフを更新する関数
        'init_func': _init,  # グラフ初期化用の関数 (今回はデータ更新用スレッドの起動)
        'fargs': (x, y),  # 関数の引数 (フレーム番号を除く)
        'interval': 250,  # グラフを更新する間隔 (ミリ秒)
    }
    anime = animation.FuncAnimation(**params)

    # グラフを表示する
    plt.show()


if __name__ == '__main__':
    main()

上記をファイルに保存して実行してみよう。

$ python sin.py

これまでに比べるとグラフの再描画間隔が長いので、ちょっとカクカクした感じでグラフが更新される。

めでたしめでたし。

Python: pandas の永続化フォーマットについて調べた

以前、このブログでは pandas の DataFrame を Pickle として保存することで読み込み速度を上げる、というテクニックを紹介した。

blog.amedama.jp

実は pandas がサポートしている永続化方式は Pickle 以外にもある。 今回は、その中でも代表的な以下の永続化フォーマットについて特性を調べると共に簡単なベンチマークを取ってみることにした。

  • Pickle
  • Feather
  • Parquet

使った環境とパッケージのバージョンは次の通り。

$ sw_vers 
ProductName:    Mac OS X
ProductVersion: 10.13.5
BuildVersion:   17F77
$ python -V
Python 3.6.5
$ pip list --format=columns | egrep "(pandas|feather-format|pyarrow)"
feather-format   0.4.0      
pandas           0.23.3     
pyarrow          0.9.0.post1

pandas の永続化フォーマットについて

本格的な説明に入る前に、まずはこの記事で扱う各フォーマットの概要について解説する。

Pickle Format

Python が標準でサポートしている直列化方式。 直列化というのは、プログラミングにおいてオブジェクトをバイト列に変換・逆変換する機能のことをいう。 バイト列に変換した内容は、ディスクにファイルなどの形で保存 (永続化) できる。

基本的に、直列化した内容は Python でしか扱うことができない。 つまり、言い換えると他のプログラミング言語や分析ツール (例えば R など) へのポータビリティはない。

その他、詳しくは以下のブログ記事で解説している。

blog.amedama.jp

Feather Format

主に Python と R の間でデータのポータビリティがあることを目指して策定されたファイルフォーマット。 行単位ではなく列単位でデータを格納する形式 (列指向・カラムナ) になっている。 今のところ pandas では試験的な導入という位置づけらしい。

Parquet Format

Apache Hadoop エコシステムでよく使われているファイルフォーマット。 こちらも Feather フォーマットと同様に、異なるプログラミング言語や分析ツールの間でポータビリティがある。

以前、このブログでも Python で Parquet フォーマットを実装した二つのパッケージについて調べたことがある。

blog.amedama.jp

下準備

それぞれの永続化フォーマットについての解説が済んだところで、次は実際に試してみる準備に入る。

まずは、今回使うパッケージ一式をインストールしておく。

$ pip install pandas ipython feather-format memory_profiler

動作確認のためのデータセットは、ベンチマークを取りたいのである程度サイズがほしい。 そこで前述した Pickle の記事と同様に Kaggle の大気汚染データセットを使うことにした。

Hazardous Air Pollutants | Kaggle

上記のページか、あるいは kaggle コマンドを使ってデータセットをダウンロードしてくる。

$ kaggle datasets download -d epa/hazardous-air-pollutants

ダウンロードしたらデータセットを展開する。

$ mv ~/.kaggle/datasets/epa/hazardous-air-pollutants/hazardous-air-pollutants.zip .
$ unzip hazardous-air-pollutants.zip 

すると、次の通り 2GB 強の CSV ファイルが手に入る。

$ du -m epa_hap_daily_summary.csv 
2349   epa_hap_daily_summary.csv

続いて IPython の REPL を起動する。

$ ipython

Python 純正の REPL を使わないのは理由は、実行時間やメモリ使用量の計測が IPython を使った方が楽ちんなので。

REPL を起動したら pandas をインポートした上で上記の CSV を読み込んでおく。

>>> import pandas as pd
>>> df = pd.read_csv('epa_hap_daily_summary.csv')

次の通りデータセットが読み込まれれば準備完了。

>>> df.head()
   state_code         ...           date_of_last_change
0          42         ...                    2015-07-22
1          48         ...                    2014-03-25
2          22         ...                    2015-07-22
3          18         ...                    2017-02-20
4           6         ...                    2015-07-22

[5 rows x 29 columns]

各永続化方式のベンチマーク

今回は、各永続化方式について以下の観点でベンチマークを取ってみる。

  • 保存にかかる時間
  • 保存したときのファイルサイズ
  • 読み込みにかかる時間
  • 読み込みに必要なメモリサイズ

先に断っておくと、今回のベンチマークはあくまで簡易的なものになっている。 本来であれば、扱うデータ型などによる特性の違いや、試行回数についても複数回の平均を取りたい。 とはいえ、そこまでやると大変なので、今回は大気汚染データセットだけを使って、かつ実行回数も 1 回だけに絞っている。 もしかすると環境や使うデータセットによっては異なる結果が出る場合もあるかもしれない。

ちなみに、あとの方に結果をまとめた表を作ってあるので読むのがめんどくさい人はそこまで読み飛ばしてもらえると。

保存にかかる時間

まずは各永続化フォーマットでファイルに保存するのにかかる時間を測ってみる。

Pickle Format

>>> %time df.to_pickle('hazardous-air-pollutants.pickle')
CPU times: user 9.31 s, sys: 9.17 s, total: 18.5 s
Wall time: 20.8 s

Feather Format

>>> %time df.to_feather('hazardous-air-pollutants.feather')
CPU times: user 1min 38s, sys: 6.76 s, total: 1min 45s
Wall time: 1min 46s

Parquet Format

>>> %time df.to_parquet('hazardous-air-pollutants.snappy.parquet')
CPU times: user 2min 2s, sys: 7.52 s, total: 2min 9s
Wall time: 2min 15s

pandas.to_parquet() 関数はデフォルトで snappy による圧縮がかかる。 念のため、圧縮をかけない場合についても確認しておこう。

>>> %time df.to_parquet('hazardous-air-pollutants.none.parquet', compression=None)
CPU times: user 2min 3s, sys: 9.19 s, total: 2min 12s
Wall time: 2min 30s

あまり変わらないか、むしろ遅くなっている。

ちなみに Pickle フォーマットも GZip で圧縮することはできる。 ただし、試してみたところ保存に時間がかかりすぎて実用に耐えないという判断で検証するのを止めた。

上記から、保存に関しては Pickle フォーマットが最も早いことが分かった。

終わったら、一旦 IPython の REPL から抜けておこう。

>>> exit()

保存したときのファイルサイズ

続いては各永続化フォーマットで保存したときのディスク上のサイズを確認する。

まず、元々の CSV は 2.3GB だった。

$ du -m epa_hap_daily_summary.csv 
2349   epa_hap_daily_summary.csv

各永続化フォーマットにおけるサイズは次の通り。 圧縮の有無に関わらず Parquet が抜きん出て小さいことが分かる。

$ du -m hazardous-air-pollutants.* | grep -v zip$
3116   hazardous-air-pollutants.feather
260    hazardous-air-pollutants.none.parquet
1501   hazardous-air-pollutants.pickle
228    hazardous-air-pollutants.snappy.parquet

それ以外だと Pickle はやや小さくなっているのに対し Feather は元の CSV よりもサイズが増えてしまっていることが分かる。

読み込みにかかる時間

もう一度 IPython の REPL を起動して pandas をインポートしておこう。

$ ipython
>>> import pandas as pd

準備ができたら各永続化フォーマットで保存したファイルを読み込むのにかかる時間を測ってみる。

Pickle Format

>>> %time df = pd.read_pickle('hazardous-air-pollutants.pickle')
CPU times: user 4.51 s, sys: 4.07 s, total: 8.58 s
Wall time: 9.03 s

Feather Format

>>> %time df = pd.read_feather('hazardous-air-pollutants.feather')
CPU times: user 18.2 s, sys: 17.3 s, total: 35.4 s
Wall time: 1min 8s

Parquet Format

%time df = pd.read_parquet('hazardous-air-pollutants.snappy.parquet')
CPU times: user 27.4 s, sys: 34.4 s, total: 1min 1s
Wall time: 1min 31s

念のため snappy 圧縮していないものについても。

%time df = pd.read_parquet('hazardous-air-pollutants.none.parquet')
CPU times: user 27.3 s, sys: 34.7 s, total: 1min 2s
Wall time: 1min 18s

わずかに読み込み時間が短くなっている。

終わったら一旦 IPython の REPL を終了しておこう。

>>> exit()

読み込みに必要なメモリサイズ

続いては永続化されたファイルを読み込む際に必要なメモリのサイズを確認しておく。 これが大きいと、読み込むときにスラッシングを起こしたり最悪の場合はメモリに乗り切らずにエラーになる恐れがある。

まずは IPython を起動して pandas をインポートする。

$ ipython
>>> import pandas as pd

それができたら、続いては %load_ext で memory_profiler を読み込む。 これで %memit マジックコマンドが使えるようになる。

>>> %load_ext memory_profiler

あとは各永続化フォーマットで保存されたファイルを読み込む際のメモリ使用量を %memit マジックコマンドで確認していく。

Pickle Format

>>> %memit df = pd.read_pickle('hazardous-air-pollutants.pickle')
peak memory: 2679.70 MiB, increment: 2617.24 MiB

読み込みが終わったら GC を呼んで読み込んだ DataFrame を開放しておく。

>>> del df; import gc; gc.collect()

仕様的には手動で呼び出しても本当にオブジェクトが回収される保証はないはず。 とはいえ、まあこの状況ならまず間違いなく回収されるでしょう。

Feather Format

>>> %memit df = pd.read_feather('hazardous-air-pollutants.feather')
peak memory: 3614.03 MiB, increment: 3549.74 MiB

終わったら GC を呼ぶ。

>>> del df; import gc; gc.collect()

Parquet Format

>>> %memit df = pd.read_parquet('hazardous-air-pollutants.snappy.parquet')
peak memory: 3741.76 MiB, increment: 3677.94 MiB

終わったら GC を呼ぶ。

>>> del df; import gc; gc.collect()

念のため snappy 圧縮なしのパターンについても。

%memit df = pd.read_parquet('hazardous-air-pollutants.none.parquet')
peak memory: 3771.86 MiB, increment: 3707.23 MiB

ベンチマーク結果

ベンチマークの結果を表にまとめてみる。

項目 Pickle Feather Parquet (snappy)
保存にかかる時間 20s 1m46s 2m15s
保存したときのサイズ 1501MB 3116MB 228MB
読み込みにかかる時間 9s 1m8s 1m18s
読み込みに必要なメモリサイズ 2679MB 3614MB 3741MB

保存したときのサイズは Parquet が最も小さかったけど、それ以外は全て Pickle の性能が優れていた。

各永続化フォーマットを使う際の注意点

Pickle Format

Pickle に関しては、保存する DataFrame に対する制約は特にない。 その以外で特に気をつける点としては Pickle フォーマットにはバージョンがあるところくらいかな。 これは、異なるバージョンの Python 間で同じ Pickle ファイルを扱う場合には、フォーマットのバージョンに気をつける必要がある。 あとは前述した通り Python 以外から扱えないので基本的に他の環境へのポータビリティがない。

その他、細かい点については以下のブログ記事にまとめてある。

blog.amedama.jp

Feather Format

pandas の DataFrame を Feather フォーマットで保存する際の注意点は以下のページに記載されている。

IO Tools (Text, CSV, HDF5, …) — pandas 0.23.3 documentation

以下に、要点をコードで解説する。

インデックスがあると保存できない

Feather フォーマットを使う場合、DataFrame にインデックスが設定されていると永続化できない。

物は試しということで、まずはインデックスを指定した DataFrame を用意しよう。

>>> import pandas as pd
>>> data = [(10, 'Alice'), (20, 'Bob'), (30, 'Carol')]
>>> df = pd.DataFrame(data, columns=['id', 'name'])
>>> df = df.set_index('id')

上記を Feather フォーマットで保存してみる、と以下のような例外になる。

>>> df.to_feather('example.feather')
Traceback (most recent call last):
...(省略)...
ValueError: feather does not support serializing a non-default index for the index; you can .reset_index() to make the index into column(s)

もし、インデックスの指定がある DataFrame を Feather で扱いたいとしたら次のようにする。 一旦インデックスをリセットして、カラムの形で持っておけば大丈夫。

>>> df = df.reset_index()
>>> df.to_feather('example.feather')

とはいえ、ちょっとめんどくさいね。

重複するカラム名があると保存できない

pandas の DataFrame には重複するカラム名を保存できる。 ただし Feather フォーマットでは重複するカラム名がある DataFrame を永続化できない。

実際に試してみよう。 同じカラム名がある DataFrame を作る。

>>> import pandas as pd
>>> data = [('a', 'b'), ('c', 'd'), ('e', 'f')]
>>> df = pd.DataFrame(data, columns=['feature_name', 'feature_name'])

この通り、ちゃんと作れる。

>>> df
  feature_name feature_name
0            a            b
1            c            d
2            e            f

しかし、永続化しようとすると、この通り例外になる。

>>> df.to_feather('example.feather')
Traceback (most recent call last):
...(省略)...
ValueError: cannot serialize duplicate column names
保存できない型がある

Feather フォーマットでは保存できない型が存在する。 例えば Period 型なんかは典型のようだ。

実際に試してみよう。

>>> import pandas as pd
>>> data = [
...     pd.Period('2018-01-01'),
...     pd.Period('2018-01-02'),
...     pd.Period('2018-01-03'),
... ]
>>> df = pd.DataFrame(data, columns=['periods'])

保存しようとすると、次の通り例外になる。

>>> df.to_feather('example.feather')
Traceback (most recent call last):
...(省略)...
pyarrow.lib.ArrowInvalid: Error inferring Arrow type for Python object array. Got Python object of type Period but can only handle these types: string, bool, float, int, date, time, decimal, list, array

上記のエラーメッセージには取り扱える型の一覧も表示されている。

Parquet Format

続いて Parquet フォーマットを扱う上での注意点について。 公式では以下のページにまとめられている。

IO Tools (Text, CSV, HDF5, …) — pandas 0.23.3 documentation

以下、要点をコードで確認していこう。

重複するカラム名があると保存できない

これは先ほど紹介した Feather フォーマットと同様。 そもそもファイルフォーマットとして、重複したカラム名を想定していないんだろう。

>>> import pandas as pd
>>> data = [('a', 'b'), ('c', 'd'), ('e', 'f')]
>>> df = pd.DataFrame(data, columns=['feature_name', 'feature_name'])

省略すると、先ほどと同じ例外になる。

>>> df.to_parquet('example.parquet')
Traceback (most recent call last):
...(省略)...
ValueError: Duplicate column names found: ['feature_name', 'feature_name']
カテゴリカル変数を保存して読み出すとオブジェクト型になる

pandas にはカテゴリカル変数が型として用意されている。 この型は Parquet フォーマットで保存することはできるんだけど、読み出すときにオブジェクト型になってしまう。

実際に試してみよう。 カテゴリカル変数が含まれるものを保存してもエラーにはならない。

>>> import pandas as pd
>>> data = [('a'), ('b'), ('c')]
>>> df = pd.DataFrame(data, columns=['categorical_name'], dtype='category')
>>> df.dtypes
categorical_name    category
dtype: object
>>> df.to_parquet('example.parquet')

ただし保存されたものを読み出してみると? なんとオブジェクト型になってしまっている。

>>> df = pd.read_parquet('example.parquet')
>>> df.dtypes
categorical_name    object
dtype: object

これは知らないとハマりそう。

保存できない型がある

これは先ほどの Feather フォーマットと同様。

例えば Period 型を保存しようとすると例外になる。

>>> import pandas as pd
>>> data = [
...     pd.Period('2018-01-01'),
...     pd.Period('2018-01-02'),
...     pd.Period('2018-01-03'),
... ]
>>> df = pd.DataFrame(data, columns=['periods'])
>>> df.to_parquet('example.parquet')
Traceback (most recent call last):
...(省略)...
pyarrow.lib.ArrowInvalid: Error inferring Arrow type for Python object array. Got Python object of type Period but can only handle these types: string, bool, float, int, date, time, decimal, list, array
マルチインデックスに使うカラム名は文字列型である必要がある

pandas には複数のカラムを用いたインデックスが作れる。 ただし、そのインデックスに使うカラム名の型は必ず文字列じゃないとだめ。

例えばマルチインデックスの一つとして名前が数値のカラムを指定してみよう。 以下でいう 1970 というカラム。

>>> import pandas as pd
>>> data = [
...     ('Alice', 47, 70),
...     ('Alice', 48, 80),
...     ('Bob', 47, 50),
...     ('Bob', 48, 45),
...     ('Carol', 47, 60),
...     ('Carol', 48, 65),
... ]
>>> df = pd.DataFrame(data, columns=['name', 1970, 'score'])
>>> df = df.set_index(['name', 1970])

この通り、ちゃんと数値の名前がついている。

>>> df.index
MultiIndex(levels=[['Alice', 'Bob', 'Carol'], [47, 48]],
           labels=[[0, 0, 1, 1, 2, 2], [0, 1, 0, 1, 0, 1]],
           names=['name', 1970])

これを Parquet フォーマットで保存しようとするとエラーになる。

>>> df.to_parquet('example.parquet')
Traceback (most recent call last):
...(省略)...
ValueError: Index level names must be strings

まあ文字列以外の名前をつけようとするなんて、そうそうないとは思うけど。

まとめ

今回は pandas の代表的な永続化フォーマットについて、その特性を調べると共に簡易なベンチマークを取ってみた。 簡易的かつ、今回使った環境では、という断りを入れた上でまとめると、次の通り。

  • ディスク上のファイルサイズという観点では Parquet フォーマットが優れていた
  • それ以外の観点 (読み書きの速度やメモリ消費量) では Pickle フォーマットが優れていた
  • pandas と Parquet / Feather フォーマットを組み合わせて使う場合には、いくつかの制限事項がある

以上から、使い分けについては次のことが言えると思う。

  • 基本的には Pickle フォーマットを使っていれば良さそう
  • ディスク上のファイルサイズに制約があれば Parquet フォーマットが良さそう
  • 他の環境とのポータビリティが必要なときは Parquet もしくは Feather フォーマットを使う

いじょう。

shellcheck でシェルスクリプトのコードの質をチェックする

正しく動作するシェルスクリプトを書くのは難しい。 できれば書きたくないけど、そうもいかない。 そんなとき心の支えになりそうなのが今回紹介する shellcheck というツール。 これはシェルスクリプトにおける Linter (リンター) で、まずい書き方をしているとそれを教えてくれる。

使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.13.5
BuildVersion:   17F77

下準備

まずは shellcheck をインストールする。 macOS であれば Homebrew でインストールするのが楽ちん。

$ brew install shellcheck

インストールすると shellcheck コマンドが使えるようになる。

$ shellcheck -V
ShellCheck - shell script analysis tool
version: 0.5.0
license: GNU General Public License, version 3
website: https://www.shellcheck.net

使ってみる

試しに、意図的にまずいシェルスクリプトを書いてそれをチェックしてみよう。

例えば、こんなコードを用意してみる。 一見すると、特に問題はなさそうな気もする。

$ cat << 'EOF' > example.sh 
#!/bin/sh

echo Time: `date`
EOF

実行しても、ちゃんと動くし。

$ sh example.sh 
Time: 201878日 日曜日 185226秒 JST

とはいえ、これを shellcheck 先生にかけると、次のようにいくつか指摘を受ける。 例えば SC2046 はクオートが足りないという指摘で、SC2006 はインラインのコマンド実行には `` よりも $() を使うべきという指摘。

$ shellcheck example.sh

In example.sh line 3:
echo Time: `date`
           ^-- SC2046: Quote this to prevent word splitting.
           ^-- SC2006: Use $(..) instead of legacy `..`.

上記の指摘を直してみよう。

$ cat << 'EOF' > example.sh 
#!/bin/sh

echo "Time: $(date)"
EOF

再度実行すると、今度は何も言われない。

$ shellcheck example.sh 

他にも、例えばシングルクオートで囲んでいるので展開されない変数があるパターンとか。

$ cat << 'EOF' > example.sh 
#!/bin/sh

echo 'Can not expand with single quote: $PATH'
EOF

これも、次の通りちゃんと教えてくれる。

$ shellcheck example.sh

In example.sh line 3:
echo 'Can not expand with single quote: $PATH'
     ^-- SC2016: Expressions don't expand in single quotes, use double quotes for that.

あとは、地味にやりがちな lsfor 文に突っ込んじゃうやつとか。 これはファイル名やディレクトリ名にスペースが入っていると動かない。

$ cat << 'EOF' > example.sh 
#!/bin/sh

for f in $(ls); do
  echo "${f}"
done
EOF

これもちゃんと教えてくれる。

$ shellcheck example.sh

In example.sh line 3:
for f in $(ls); do
         ^-- SC2045: Iterating over ls output is fragile. Use globs.

ばっちり。

めでたしめでたし。

Python: scikit-learn の Pipeline を使ってみる

機械学習では、元のデータセットに対して前処理や推論フェーズが何段にも重なることがある。 scikit-learn には、そういった何段にも重なった処理を表現しやすくするために Pipeline という機能が備わっている。 今回は、その Pipeline を使ってみることにする。

使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.13.5
BuildVersion:   17F77
$ python -V                           
Python 3.6.5

下準備

まずは scikit-learn と必要なパッケージをあらかじめインストールしておく。

$ pip install numpy scipy scikit-learn

続いて Python のインタプリタを起動しておく。

$ python

データセットについては Iris データセットを使う。

>>> from sklearn import datasets
>>> iris = datasets.load_iris()
>>> X, y = iris.data, iris.target

これで下準備はできた。

Pipeline 機能を使ってみる

Pipeline では、その名のごとく一連の処理をパイプラインで表現する。 具体的には、リストに scikit-learn のオブジェクトを適用する順番で格納しておく。

以下にサンプルコードのパイプラインを示す。 このパイプラインでは、まず前処理として主成分分析 (PCA) した上で、それをランダムフォレストで分類できる。 リストには処理の名前と、その処理を実行するオブジェクトを対にしたタプルを入れる。

>>> from sklearn.pipeline import Pipeline
>>> from sklearn.decomposition import PCA
>>> from sklearn.ensemble import RandomForestClassifier
>>> steps = [
...     ('pca', PCA()),
...     ('rf', RandomForestClassifier())
... ]
>>> pipeline = Pipeline(steps=steps)

PCA とランダムフォレストを選んだこと自体には特に意味はなくて、あくまでこんな風にできるよという例。

Pipeline オブジェクトは scikit-learn のオブジェクトが一般的に持っているインターフェースを備えている。 なので、以下のような感じで fit() メソッドで学習した上で predict() メソッドで推論、といういつもの流れが使える。

>>> from sklearn.model_selection import train_test_split
>>> X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33, random_state=42)
>>> pipeline.fit(X_train, y_train)
>>> y_pred = pipeline.predict(X_test)

精度自体は今回の目的ではないけど、一応計算しておく。

>>> from sklearn.metrics import accuracy_score
>>> accuracy_score(y_test, y_pred)
0.98

GridSearchCV と組み合わせて使う

Pipeline は GridSearchCV と組み合わせることもできる。 GridSearchCV を使うと、ハイパーパラメータの中から性能の良い組み合わせを自動的に探してくれる。

Pipeline と GridSearchCV を組み合わせて使うときは、パラメータ名の指定に工夫が必要になる。 具体的には <処理の名前>__<パラメータ名> という形式で指定する。

>>> from sklearn.model_selection import GridSearchCV
>>> params = {
...     'pca__n_components': range(1, X.shape[1]),
...     'rf__n_estimators': [10, 100, 1000],
... }
>>> grid_search_cv = GridSearchCV(pipeline,  param_grid=params, cv=4)

ふつうの scikit-learn オブジェクトと GridSearchCV を組み合わせるときは単純にパラメータ名だけで良かった。 それに対して Pipeline には複数のオブジェクトが含まれるため、それを識別するための名前が追加で必要になる。

GridSearchCV#fit() メソッドでハイパーパラメータの探索処理が走る。

>>> grid_search_cv.fit(X, y)

処理が終わったら best_params_ を参照することで最適なパラメータの値が分かる。

>>> grid_search_cv.best_params_
{'pca__n_components': 2, 'rf__n_estimators': 1000}

Pipeline に組み込むオブジェクトを自作する

ここまでは scikit-learn オブジェクトをそのまま Pipeline に組み込むパターンを試した。 次は Pipeline に自作した前処理や推論フェーズの処理を組み込んでみよう。 とはいえ、これには自作したオブジェクトに scikit-learn のインターフェースをダックタイピングで用意するだけで良い。

まずは前処理をするオブジェクトから。 この処理では受け取った値を 2 で割ったものに変換する。 もちろん、この処理自体には特に意味はない。

>>> class MyPreprocessor(object):
...     """前処理 (最低限 fit, transform メソッドが必要)"""
...     def fit(self, X, y):
...         """学習フェーズ"""
...         return self
...     def transform(self, X):
...         """変換フェーズ"""
...         # 2 で割った余りを返す
...         return X % 2
...     def predict(self, X):
...         """分類フェーズ"""
...         # 何もせずに素通しする
...         return X
... 

続いては Pipeline の最後段で推論をするオブジェクト。 このオブジェクトは受け取った値が 0 か非 0 かで True/False を返す処理をする。 こちらも、もちろん処理自体には特に意味はない。

>>> class MyClassifier(object):
...     """回帰・分類処理 (fit メソッドが必要)"""
...     def fit(self, X, y):
...         """学習フェーズ"""
...         return self
...     def predict(self, X):
...         """分類フェーズ"""
...         # X が 0 か非 0 かで True/False を返す
...         return X == 0
... 

上記で用意した自作オブジェクトを組み込んだ Pipeline を作る。 このパイプラインでは値を 2 で割った余りに変換した上で、それが 0 なら True をそうでなければ False を返す。 これはようするに、一連の処理で引数が偶数か奇数かを判定することになる。

>>> steps = [
...     ('pre', MyPreprocessor()),
...     ('clf', MyClassifier()),
... ]
>>> pipeline = Pipeline(steps=steps)

データとしては 0 から 9 までのリストを使う。

>>> import numpy as np
>>> X = np.arange(10)

Pipeline にデータを渡してみよう。 今回に関しては実は fit() 関数は使ってないんだけど、とりあえず機械学習なら普通は使うので呼んでおく。

>>> pipeline.fit(X)
>>> pipeline.predict(X)
array([ True, False,  True, False,  True, False,  True, False,  True,
       False])

うまいこと受け取った引数が偶数か奇数かを判定できている。

めでたしめでたし。

Python: pandas の DataFrame を scikit-learn で KFold するときの注意点

今回は pandas の DataFrame を scikitl-learn で交差検証しようとしてハマった話について。 だいぶ平凡なミスなんだけど、またやるとこわいので自分用にメモしておく。

使った環境は次の通り。

$ sw_vers 
ProductName:    Mac OS X
ProductVersion: 10.13.5
BuildVersion:   17F77
$ python -V                         
Python 3.6.5

下準備

まずは pandas と scikit-learn をインストールしておく。

$ pip install pandas scikit-learn scipy

Python のインタプリタを起動する。

$ python

なんか適当に DataFrame を作っておく。

>>> import pandas as pd
>>> data = [
...     ('Ant'),
...     ('Beetle'),
...     ('Cat'),
...     ('Deer'),
...     ('Eagle'),
...     ('Fox'),
... ]
>>> columns = ['name']
>>> df = pd.DataFrame(data, columns=columns)

できた DataFrame はこんな感じ。

>>> df
     name
0     Ant
1  Beetle
2     Cat
3    Deer
4   Eagle
5     Fox

特に指定しない限り、デフォルトではインデックスとして 0 から始まる整数が連番で振られる。

sklearn.model_selection.KFold の使い方

scikit-learn で交差検証するとき基本は KFold クラスを使う。 このクラスはインスタンス化するときに分割数を指定し、その上で KFold#split() メソッドに分割するものを渡す。 返り値としてはイテラブルなオブジェクトが返ってきて、それぞれ学習用データと検証用データ用のインデックスが取り出せる。

>>> from sklearn.model_selection import KFold
>>> kf = KFold(n_splits=3)
>>> for train, test in kf.split(df):
...     print(train, test)
... 
[2 3 4 5] [0 1]
[0 1 4 5] [2 3]
[0 1 2 3] [4 5]

pandas との連携 (ダメなパターン)

KFold#split() で手に入るのは単なるインデックスなので、それを元に DataFrame から対象データを抽出しないといけない。 このとき、次のようなコードを書いてしまうと一見すると上手くいっているように見えて後でハマることになる。

>>> from sklearn.model_selection import KFold
>>> kf = KFold(n_splits=3)
>>> for train, test in kf.split(df):
...     # これはやっちゃダメ!
...     train_df = df[df.index.isin(train)]
...     test_df = df[df.index.isin(test)]
...     # 結果確認 (一見すると上手くいっているように見える)
...     print('(train)', train_df)
...     print('(test)', test_df)
...     print('-----')
... 
(train)     name
2    Cat
3   Deer
4  Eagle
5    Fox
(test)      name
0     Ant
1  Beetle
-----
(train)      name
0     Ant
1  Beetle
4   Eagle
5     Fox
(test)    name
2   Cat
3  Deer
-----
(train)      name
0     Ant
1  Beetle
2     Cat
3    Deer
(test)     name
4  Eagle
5    Fox
-----

実は、上記はインデックスが 0 からの連番で振られているために動いているに過ぎない。 その前提が変わると途端に動かなくなる。 試しにインデックスの番号を変更してみよう。

>>> # インデックスを 0 から始まる連番から変える (以下なら 10, 20, 30...)
... df.index = df.index * 10 + 10

変更後のインデックスはこんな感じ。

>>> df
      name
10     Ant
20  Beetle
30     Cat
40    Deer
50   Eagle
60     Fox

先ほどと全く同じコードを使って動作を確認してみよう。

>>> kf = KFold(n_splits=3)
>>> for train, test in kf.split(df):
...     # OMG!!
...     train_df = df[df.index.isin(train)]
...     test_df = df[df.index.isin(test)]
...     # 結果確認
...     print('(train)', train_df)
...     print('(test)', test_df)
...     print('-----')
... 
(train) Empty DataFrame
Columns: [name]
Index: []
(test) Empty DataFrame
Columns: [name]
Index: []
-----
(train) Empty DataFrame
Columns: [name]
Index: []
(test) Empty DataFrame
Columns: [name]
Index: []
-----
(train) Empty DataFrame
Columns: [name]
Index: []
(test) Empty DataFrame
Columns: [name]
Index: []
-----

中身が全て空っぽになってしまっている!

pandas との連携 (正解)

上記のようなパターンでも動くようにするには DataFrame の絞り込みで DataFrame#iloc を使う。 これなら、本来のインデックスの値ではなく中身の順序にもとづいたインデックスで絞り込みができる。

>>> kf = KFold(n_splits=3)
>>> for train, test in kf.split(df):
...     # これなら上手くいく
...     train_df = df.iloc[train]
...     test_df = df.iloc[test]
...     print('(train)', train_df)
...     print('(test)', test_df)
...     print('-----')
... 
(train)      name
30    Cat
40   Deer
50  Eagle
60    Fox
(test)       name
10     Ant
20  Beetle
-----
(train)       name
10     Ant
20  Beetle
50   Eagle
60     Fox
(test)     name
30   Cat
40  Deer
-----
(train)       name
10     Ant
20  Beetle
30     Cat
40    Deer
(test)      name
50  Eagle
60    Fox
-----

ばっちり。

Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理

Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理

Pythonデータサイエンスハンドブック ―Jupyter、NumPy、pandas、Matplotlib、scikit-learnを使ったデータ分析、機械学習

Pythonデータサイエンスハンドブック ―Jupyter、NumPy、pandas、Matplotlib、scikit-learnを使ったデータ分析、機械学習

色々と調べた末にダイヤモンドを個人輸入した話

今回は、いつもと違って技術系ではない話。

少し前のことだけど、一身上の都合によりダイヤモンドが必要になる機会があった。 正直なところ、自分自身は炭素の結晶に興味はない。 また、その資産価値についても (現時点でも) 懐疑的に見ている。 とはいえ、世間的にはずいぶんと高値で取引されているらしい。 もし購入するとなれば、なかなか高い買い物になる。 そこで色々と調べた結果、最終的には海外から個人輸入することになった。

今回は、それに至る経緯や、調べた内容についてせっかくなのでまとめてみることにする。 先に断っておくと、かなり長い。 この記事で伝えたいことの本質ではないけど、安くダイヤモンドを購入する方法についても途中に記載している。 もし、今後同じようにダイヤモンドが必要になる人にとって参考になれば幸いに思う。

TL;DR

  • 日本のジュエリーショップで販売されているダイヤモンドは高い
    • 具体的には、世界的な市場取引価格から大きく乖離した値付けがされている
  • 今は天然ダイヤモンドを購入するにはおそらく最悪のタイミング
    • もう少しで宝飾用途の人工ダイヤモンドが安価に手に入るようになるため
  • 全てを理解した上でそれでも購入するなら、なるべく安価な販売チャネルを選ぶのが望ましい

はじまり

きっかけは、パートナーがどうやらダイヤモンドの指輪を欲しているらしい、ということを理解したため。 そこで、まずはいくつかのジュエリーショップに足を運んでみた。

ジュエリーショップでは、指輪にはまっていない状態のダイヤモンドが販売されている。 その状態をルース (裸石) と呼び、目当てのデザインの枠 (指輪部分のこと) と組み合わせて購入する。 ちなみに、ルースは和製英語らしい。

指輪において、中央の大きな宝石はセンターストーンという。 日本の婚約指輪では、センターストーンにダイヤモンド、枠がプラチナという組み合わせが一般的。

そして、いくつかのお店でダイヤモンドを眺めていたところ気になる点があった。 それは、似たようなグレードのルースでも価格がお店によってずいぶんと異なるところ。 一体、その値段の違いはどこからくるのだろうか?

先に結論から書くと、日本の一般的なジュエリーショップが販売するダイヤモンドのルースは高い。 ここで「高い」というのは、具体的には世界的な市場取引価格から乖離した値付けになっている、という意味を指す。 とはいえ、値段が高いの安いのと言うには、まずダイヤモンドの前提知識から知っておく必要がある。

ダイヤモンドについて

ダイヤモンドは、炭素の原子が集まって結晶化したもの。 天然のものであれば、地球内部の高温・高圧の環境下で生成される。

光を反射しやすいので、宝飾用途で用いられる。 また、その硬さから工業用途でも需要がある。 同じ宝飾用途で用いられる貴金属と異なるのは、炭素という元素が地球上に比較的ありふれているという点。

詳しくは後述するものの、工業用途では人工ダイヤモンドの利用が一般的になっている。 また、宝飾用途でも技術の進歩と共に一部で利用が進みつつある。

宝飾用途のダイヤモンドについて

ダイヤモンドは、生成された環境と加工方法によって個体それぞれ見た目が異なる。 宝飾用途のダイヤモンドの場合、それを 4C という基準で評価するのが世界的に一般的となっている。 この名前は、評価の基準の頭文字が 4 種類で、全て C から始まるためにそう呼ばれている。 評価は世界にいくつかある鑑定機関が有償で行う。 その結果は鑑定書という形で石に付与される。

Carat (カラット)

ルースの重さを用いた評価基準。

勘違いしやすいのは「大きさ」ではなく「重さ」の基準というところ。 具体的には 0.2 グラム = 1 カラットになる。

カラットと価格の関係は非線形になっている。 どういうことかというと、カラットの数字が大きくなるほど、値段の上がり幅も大きくなる。 これは、大きなダイヤモンドの原石ほど算出されにくいという希少性がその根拠らしい。

以下は、とあるジュエリーショップに記載されているダイヤモンドの価格表をスクレイピングして散布図にしたもの。 x 軸にカラット数、y 軸に価格をプロットしている。

f:id:momijiame:20180617181002p:plain

あくまで一例ではあるものの、カラット数が大きくなるにつれて価格が非線形に高くなる関係が見て取れる。

ところで、面積は長さの二乗に比例するのに対して体積は三乗に比例する。 だとすると、カラットが大きくなるほど見た目での差異は判別しにくくなるといえる。 指輪の状態で視認できるのは基本的に表面部分だけなので。 なのに大きくなるほど値段の上がり幅は大きいということになる。

Color (カラー)

ルースの色を用いた評価基準。

ダイヤモンドは、それに含まれる不純物や結晶構造の歪みなどによって色がつく。 一般的な天然ダイヤの場合、窒素の含有によって黄色くなったものが多い。

カラーはアルファベットで表される。 より無色透明に近い D カラーから E, F, G ... Z と進むにつれて黄色くなる。 一般的には無色透明に近いほど値段は高くなり、黄色になるほど安い。 以下のページを見るとイメージがつかみやすい。 日本で指輪の用途で用いられるのは G ~ I 以上のものが多いような印象を受ける。

GIA 4Cs Color D-to-Z

ただし、これには例外もあって、発色の良い色がしっかりついたものには反対に高い値段がつく場合がある。 ようするに D ~ Z のスケールでは測れないもので、これはファンシーカラーダイヤモンドと呼ばれている。 カラーダイヤモンドには黄色以外にもブルーやピンクなど、様々な色がある。

ファンシーカラーダイヤモンドの説明

珍しい色のファンシーカラーダイヤモンドで大きなものは、それこそ大衆が購入できような値段ではない。 ただし、ごく小さなもの (メレーダイヤモンド) であれば、手が届く場合もある。

Clarity (クラリティ)

ルースの内部にある内包物や傷を用いた評価基準。

ダイヤモンドは、多くの場合に内包物や傷、空洞などの見た目で分かる欠陥を持っている。 それらの欠陥を、見えやすさという基準で評価したもの。 以下のグレードに分かれていて、下に行くほど内包物が視認しやすいことを表している。

  • Flawless (FL)
  • Internally Flawless (IF)
  • Very, Very Slightly Included (VVS1, VVS2)
  • Very Slightly Included (VS1, VS2)
  • Slightly Included (SI1, SI2)
  • Included (I1, I2, I3)

内包物や傷がない・あるいはほとんどない FL や IF は数が少なく、身につけるジュエリーに供される機会はさほど多くない。 定義の上では VS2 以前であれば 10 倍ルーペを使わない限り内包物は視認が難しいとされる。 日本で指輪の用途として供されるのは VVS1 ~ VS2 のレンジが多いような印象を受ける。

GIA 4Cs Clarity

Cut (カット)

ルースの形状を用いた評価基準。

ダイヤモンドは原石のままでは光り輝かないため、カッティングや研磨といった加工を必要とする。 まず、大まかな形状としてラウンドブリリアントやペアー、オーバル、ハートなどがある。 最も一般的なのはラウンドブリリアントで、それ以外はパートナーがもし興味があれば、という感じ。

で、原石をその形状に加工した結果がどれだけ上手くいったかがカットの評価基準になる。 ただし、カットは単一の何かを見た評価項目ではなく、複合的なものとされている。 例を挙げると、カットでは各部の理想的な比率などがあらかじめ決まっている。 そのため、それにどこまで近づけることができたかは評価の一つとなる。 それ以外にも、入射した光をどれだけ反射するかも評価の一つとされる。

カットは以下のグレードに分かれていて、下にいくほど評価が低い。

  • Excellent
  • Very Good
  • Good
  • Fair
  • Poor

日本で指輪の用途として用いられるのは、ほとんど例外なく Excellent になる。 ただ、実際に Good あたりと見比べても、正直外見的な区別はつかない。 「言われてみれば少し輝きが弱いような気がしないでもない?」という程度のもの。 まあ、正直これは他の 3C でも同じことが言える。

4C 以外の評価基準

加えて、上記の 4C 以外にも補助的な評価基準がいくつかある。

Polish (ポリッシュ) / Symmetry (シンメトリー)

カットに近い概念として Polish (ポリッシュ) と Symmetry (シンメトリー) という評価基準がある。 ポリッシュは研磨、シンメトリーは各部位の対称性を元にした評価基準になっている。 グレードはカットと同様に Excellent ~ Poor までに分かれている。 販売するお店によってはカット、ポリッシュ、シンメトリーがすべて Excellent のダイヤモンドを 3EX (Triple Excellent) と呼ぶこともある。

また、これもお店によるものの 3EX に加えて H&C (Heart and Cupid) という単語がさらに後ろに付与されることもある。 これは、上手くカットしたルースを下から見るとハートが、上から見ると矢のような形が見えることに起因している。 ただし、これは世界的な評価基準ではなく、単なる付加価値をつけるための宣伝文句なので無視して良い。

Fluorescence (蛍光性)

天然ダイヤモンドは紫外線をあてると発光するものが多く、その光の強さや色を用いた評価基準。

日本においては蛍光性が弱いもしくはない方が値段が高い傾向にあるとされる。 ほとんどの場合に蛍光色は青色で、これは結晶の中に窒素が含まれる場合とされる。 蛍光の仕方は、後ほど人工ダイヤモンドのところで説明するダイヤモンドのタイプに関係する。

鑑定機関

上記の 4C やその他の基準を評価するのが鑑定機関になる。 鑑定機関は依頼を受けて有償でダイヤモンドを評価し、その結果を鑑定書の形で石に付与する。

鑑定機関は世界各国にあるけど GIA (Gemological Institute Of America: 米国宝石学会) という所がデファクトスタンダードになっている。 日本国内だと CGL (Central Gem Laboratory: 中央宝石研究所) が有名だけど、基本的に国内でしか通用しない。

鑑定書のフォーマットも機関ごとにいくつかあって、高いものほどより詳細なレポートが載る。 鑑定書つきの宝石は、少なくとも鑑定料 (数千円前後) の分はついていないものよりも高くなる。

ダイヤモンドの質を評価する機関なので、後述する人工ダイヤモンドの検出技術には深い関心を寄せており知見を有する。

ダイヤモンドの値段の決まり方

ダイヤモンドは市場において売買する人が自由に価格を決めることができる。 ただし、市場に供給されるダイヤモンドの質と量は、カルテルによって統制を受けているとされる。 供給が統制されることで、市場で取引されるダイヤモンドの価格は安定することになる。

また、カルテルの中でも特に有名なデビアス社はラパポート・ダイヤモンド・レポートという業界誌を発行している。 その中には、各グレード (4C) ごとのダイヤモンドの取引相場が記載されている。 つまり、ダイヤモンドは 4C が分かれば大体の値段が決まるようになっている。

世界中のダイヤモンドの取引価格は、供給量や業界誌といった形で直接的・間接的に影響を受ける。 それらが良いか悪いかは別として、少なくともダイヤモンドには世界的な相場があるということ。

では、日本のジュエリーショップに並ぶダイヤモンドは相場に則した価格になっているか、というとなっていない。 ほとんどの場合、相場を知っていればおよそ買わないような数倍の値付けが一般的となっている。

ジュエリーショップでは、あの手この手でその値付けを正当化する材料を用意する。 例えば、宝石の形で一度も他人の手に渡ったことがない点をバージンと称してみたり。 あるいは、独自の鑑定基準にもとづいて他社よりも高品質であることをアピールする。 当然ながら、それらは世界的な取引基準の上では、なんら評価の対象にはならない。

何処であれば市場価格に近い値段で手に入るのか?

少なくとも、宝石業界向けの競売や宝石問屋からの直接の購入であれば近い値付けになるようだ。

ただし、前者はそもそも業界人でないと参加できない。 そのため、購入するには参加する権利を持った人に依頼して、代わりに競り落としてもらう必要がある。 依頼には手間賃がかかるし、いつ目当てのグレードが出てくるかも分からない。 そうした都合から、このやり方は大口でなければ難しい気がする。

宝石問屋からの購入は、比較的現実的な方法かもしれない。 ただし、これには注意点が二つほどある。 まずひとつ目は、そもそも問屋は普段から一般人を相手に商売をしているわけではないという点。 この方法を試した人のレポートを読むと、店にもよるだろうけどかなりの塩対応は覚悟する必要があるようだ。 また、問屋は平日しか営業していないので、そう何度も企業勤めの人間が赴くのも難しい。 二つ目は、問屋「風」を装って一般人をターゲットにしたお店がある点。 そういったお店は問屋風を装ってはいるが実態はただの小売で、値付けも一般的なジュエリーショップと大して変わらない。 安さを求めて宝石問屋街に行くとしても、そうしたお店は避ける必要がある。

上記のような事情から、なかなか難しいなと思っていたところで次に見つけたのはマージンを抑えた小売店だった。 例えば、国内にも市場価格に比較的近い値付けをしている店舗はいくつかあるみたい。 とはいえ、地理的にさほど近くなかったりオンライン限定販売だったり、なかなかこれはと思うものがない。

で、行き着いた先は Blue Nile という米国のオンラインショップだった。 当初、さすがにインターネットで買うのは...と不安に思っていたものの、アメリカでは有名な企業らしい。 例えば NASDAQ にも上場していたりする (シンボル: NILE)。 日本から購入した人のレポートも、検索すると結構出てくる。

Web サイト経由で購入すると、だいたい 2 週間ほどでルースが手元に届いた。 注意点としては、サイトでの購入価格の額面以外に税金が 10% 弱かかる。 それでも、日本国内の一般的なジュエリーショップで買うよりは大幅に安い。 ちなみに、タイトルにある通り日本から購入すると形式としては個人輸入という形になる。

f:id:momijiame:20180517025355j:plain

f:id:momijiame:20180517025233j:plain

ただし、このエントリは決してダイヤモンドの海外通販は良いぞ、と言うことが目的ではない。 むしろ、絶対に必要でない限り今ダイヤモンドを購入するのは避けた方が良いと思う。 これは、後述する昨今の人工ダイヤモンド事情からそう考えている。 とはいえ、全ての事情を理解した上で、あえて購入するときの選択肢であればアリかもしれない。

ダイヤモンドのルースを購入する際の注意点

続いては、購入したルースを持ち込んで指輪にする際の注意点について。

まず前提として、ルースの持ち込みに対応しているお店はジュエリーショップの中でも少ない。 もし、パートナーが持ち込みに対応していないブランドや枠のデザインがお気に入りなら、あきらめる他ないと思う。

また、ルースの持ち込みが可能でも、それに別途料金がかかるお店もある。 あとは、持ち込む石のグレードについて一定の基準を設けているお店もある。 これは、例えばカラーなら G 以上でクラリティは VS2 以上でないと受け付けません、みたいな感じ。

人工ダイヤモンド事情

話は変わって人工ダイヤモンド (Synthetic diamond) について。 自分がダイヤモンドの資産価値に懐疑的な理由が人工ダイヤモンドの存在だった。 そして、その考えは今でも変わっていない。

元々、ダイヤモンドはその硬さから工業用途で需要がある。 工業用途の人工ダイヤモンドは、既にカラットあたり数百円というレベルで量産されている。 また、近年は技術の進歩と共に宝飾用途でも人工ダイヤモンドの利用が進みつつある。

特に、最近になって「マジか」というような話題があった。 ダイヤモンドカルテルでも代表的な企業であるデビアス社が、CVD 法 (後述) で作った人工ダイヤモンドを売り出すというのだ。

gigazine.net

「ライトボックス」というブランドで販売されるそれは、同じグレードの天然ダイヤモンドに対して数分の一の価格で販売される予定になっている。

現状でも、人工ダイヤモンドを売るメーカーというのはいくつかある。 例えば、米国内だと Gemesis 社や Apollo Diamond 社が販売している。

とはいえ、ダイヤモンドのドンといえるデビアス社が率先して人工ダイヤモンドを売り出すというのは、驚きだった。 きっと、日本でも米国から数年遅れで人工ダイヤモンドの波がくることになるだろう。

(2018/06/24 追記)

上記のライトボックスジュエリーは米国内で 9 月から発売されるらしい。

www.nikkei.com

(2018/6/28 追記)

ライトボックスジュエリーの公式サイトができていた。

lightboxjewelry.com

人工ダイヤモンドの作り方

ここから先の内容は、自分が興味から調べたものなので、ほとんどの人にとっては誰得だと思う。 読み飛ばしてもらって構わない。 人工ダイヤモンド関連で読み物として面白かったものは末尾に記載しておく。

以下に、宝飾用途の商業ベースで用いられる人工ダイヤモンドを作る手法を紹介する。 共通して言えるのは、以前はかかるコストとできあがりのダイヤモンドの質が釣り合わず、採算がとれなかった。 それが近年は技術の進歩と共に変化しつつある。 勘違いしやすいのは、人工ダイヤモンドといっても簡単にカラーレスで大きな原石ができるわけではないという点。

www.nikkei.com

HPHT (High Pressure High Temperature: 高温高圧合成) 法

炭素のペレットを高温・高圧で処理することで結晶化させる方法。 結晶構造が生み出される原理的には、天然ダイヤモンドと何ら変わらない。 ダイヤモンドになるまでには 1, 2 週間かかるので、その間ずっと加熱・加圧し続ける必要がある。

亡くなった人の遺骨からダイヤモンドを作る、みたいな話を聞いたことがあると思うけど、あれにはこの手法を使う。

また、天然ダイヤモンドや後述する CVD 法を使った人工ダイヤモンドの色調改善にも用いられる。

CVD (Chemical Vapor Deposition: 化学気相蒸着) 法

核となるダイヤモンドに炭素を蒸着することで少しずつ育てていく方法。 端的に言うと、水素とメタンガスの入ったチャンバーに核ダイヤを入れて、マイクロウェーブを照射すると少しずつ大きくなる。

HPHT 法と比べると、育成環境によってできあがるダイヤモンドの特徴を調整しやすい。 また、同じ育成環境であれば、できあがるダイヤモンドの質も揃えることができる。 そのため、宝飾用途の人工ダイヤモンドの生成では、こちらの手法が主流となっている。

育成速度の向上や、同時にチャンバー内で育成できるダイヤの数の増加によってコスト低減が図られてきた。 また、育成環境の調整によってできあがるダイヤモンドのグレードを改善することでより高価値を生み出せる。 天然ものなら数千万円するようなファンシーカラーダイヤモンドを人の手で作り出せるとしたら? 夢が広がる。

また、前述した通り手法として CVD 法 or HPHT 法という二者択一ではない。 CVD 法で育てたダイヤモンドを HPHT 法で調整する、といった手法も一般的になっている。

人工ダイヤモンドの判別方法

宝飾用途でも人工ダイヤが作れる、ということが分かったところで続いてはその判別について。

まず大前提として、人間の裸眼で人工ダイヤモンドの判別はできない。 これは、一般人なら無理でも鑑定士なら、といったレベルの話ではない。 人工ダイヤモンドかそうでないかの判別は、専門の機械がなければ不可能になっている。 つまり、天然ダイヤモンドと人工ダイヤモンドの外見的な差異は全くない。

それでは、専門の機械ではどのように判別しているかというと、これには色々な手法がある。 例えば、含有する不純物にもとづいた判別方法は簡易的な検査器具でも用いられる。 というのも、ダイヤモンドは含有する不純物によっていくつかの型 (タイプ) に分けられるため。

ダイヤモンドの型 (タイプ)

ダイヤモンドは、ざっくり言うと窒素を含むか含まないかで二つの型 (タイプ) に分けられる。 窒素を含むものを I 型、全く含まないものを II 型と呼ぶ。 それぞれのタイプは、窒素の含有量や窒素以外の元素を含むかで、さらに細分化されている。

天然で産出するダイヤモンドは、そのほとんど (約98%) が窒素を明瞭に含む Ia 型に分類される。 そして、人工ダイヤモンドでは今のところ Ia 型を作り出すことはできていない (とされる)。 これにもとづいて Ia 型以外、特に II 型については人工ダイヤモンドの疑いが強くなる。

そのため、人工ダイヤモンドの判別では、まず簡易的な検査器具を用いてダイヤモンドのタイプを判定する。 その上で Ia 型でないことが分かったら、それを鑑定機関のラボに送って、より詳細な検査をすることになる。 鑑定機関のラボではフォトルミネッセンス (PL) 法などを用いてダイヤモンドの特徴を調べている。 もし仮に、人工的に安価な Ia 型のダイヤモンドを作ることができれば、そのときは市場がパニックに陥ることだろう。

I 型 (Ia / Ib)

結晶構造の中に窒素原子を含むものを I 型と呼ぶ。 そして、その含有量や含有の仕方によって Ia 型と Ib 型に分けられる。

Ia 型は窒素が明瞭 (~5000ppm) に集合体として含まれるもの。 天然ダイヤモンドのほとんど (約98%) が、この型に分類される。

Ib 型は結晶構造の中に窒素原子が単独でわずかに (~500ppm) 含まれるもの。 HPHT 法で作る・あるいは処理した人工ダイヤモンドは、この型になることが多い。

II 型 (IIa / IIb)

結晶構造の中に窒素原子を含まないものを II 型と呼ぶ。 そして、ホウ素を含むか含まないかで IIa 型と IIb 型に分けられる。

IIa 型は窒素が全く含まれないもの。 CVD 法で作った人工ダイヤモンドは基本的にこれになる。

IIb 型は窒素が含まれない代わりにホウ素が含まれるもの。 その他の型のダイヤモンドが絶縁体なのに対し、この型だけは半導体の特性を持つ。

まとめ

これで、書きたかったことはだいたい書けたと思う。 伝えたかったことは、最初の TL;DR に記載した通り。

  • 日本のジュエリーショップで販売されているダイヤモンドは高い
    • 具体的には、世界的な市場取引価格から大きく乖離した値付けがされている
  • 今は天然ダイヤモンドを購入するにはおそらく最悪のタイミング
    • もう少しで宝飾用途の人工ダイヤモンドが安価に手に入るようになるため
  • 全てを理解した上でそれでも購入するなら、なるべく安価な販売チャネルを選ぶのが望ましい

経済的な合理性からいえば、今のタイミングでダイヤモンドを購入したのは明らかな誤りだったと思う。 とはいえ、今回に関しては資産を入手することが目的で購入したわけではない。 パートナーに喜んでもらう、という当初の目的は果たせた。 そうした点からは、現状で取りうるベターな選択肢は選べたように感じる。

参考

人工ダイヤモンド関連で、読み物として興味深かったものを以下に記載しておく。

合成ダイヤモンドの最新事情

CVD成長法による合成ダイヤモンド、パート1: 歴史

CVD成長法による合成ダイヤモンド、パート2: 特性

CVD合成ダイヤモンド、パート3: 検出

合成ダイヤモンド:品質の向上と鑑別の課題

HPHT成長法による大型の合成ダイヤモンドをGIA香港ラボが検査

Identifying Lab-Grown Diamonds | Research & News

CVD Synthetic Diamond Over 5 Carats Identified by GIA | Gems & Gemology

普段は技術系のことも書いてます。

Python: LightGBM でカテゴリ変数を扱ってみる

以前このブログで LightGBM を使ってみる記事を書いた。 ただ、この記事で使っている Iris データセットにはカテゴリ変数が含まれていなかった。

blog.amedama.jp

そこで、今回はマッシュルームデータセットを使ってカテゴリ変数が含まれる場合を試してみる。

使った環境は次の通り。

$ sw_vers 
ProductName:    Mac OS X
ProductVersion: 10.13.4
BuildVersion:   17E202
$ python -V
Python 3.6.5

マッシュルームデータセットについて

マッシュルームデータセットはAgaricales Agaricaceae Lepiota という品種のキノコの外見的な特徴と毒の有無を記録したもの。 具体的には、次のような特徴が全てカテゴリ変数の形で入っている。 要するに、外見的な特徴から毒の有無を判定するモデルを作りたい、ということ。

  • 毒の有無
  • 傘の形
  • 傘の表面
  • 傘の色
  • 部分的な変色の有無
  • 匂い
  • ひだの形状
  • ひだの間隔
  • ひだの大きさ
  • ひだの色
  • 柄の形
  • 柄の根本
  • 柄の表面 (つぼより上)
  • 柄の表面 (つぼより下)
  • 柄の色 (つぼより上)
  • 柄の色 (つぼより下)
  • 覆いの種類
  • 覆いの色
  • つぼの数
  • つぼの種類
  • 胞子の色
  • 群生の仕方
  • 生息地

データセットは次の Web サイトから得られる。

UCI Machine Learning Repository: Mushroom Data Set

まずは wget や curl を使ってデータセットをダウンロードしておこう。

$ wget https://archive.ics.uci.edu/ml/machine-learning-databases/mushroom/agaricus-lepiota.data

下準備

続いて、今回使う Python のパッケージをインストールしておく。

$ brew install cmake gcc@7
$ export CXX=g++-7 CC=gcc-7
$ pip install --no-binary lightgbm lightgbm pandas

インストールできたら Python のインタプリタを起動する。

$ python

データセットの CSV をパースする

続いてはダウンロードしてきた CSV を pandas の DataFrame に変換する。 ここで注目すべきは dtype として 'category' を指定しているところ。 これによってパースされた各次元がカテゴリ型として扱われる。

>>> import pandas as pd
>>> columns = [
...     'poisonous',
...     'cap-shape',
...     'cap-surface',
...     'cap-color',
...     'bruises',
...     'odor',
...     'gill-attachment',
...     'gill-spacing',
...     'gill-size',
...     'gill-color',
...     'stalk-shape',
...     'stalk-root',
...     'stalk-surface-above-ring',
...     'stalk-surface-below-ring',
...     'stalk-color-above-ring',
...     'stalk-color-below-ring',
...     'veil-type',
...     'veil-color',
...     'ring-number',
...     'ring-type',
...     'spore-print-color',
...     'population',
...     'habitat',
... ]
>>> df = pd.read_csv('agaricus-lepiota.data', header=None, names=columns, dtype='category')

上記でマッシュルームデータセットの 23 次元の特徴を持ったデータが読み込まれた。

>>> df.head()
  poisonous cap-shape cap-surface   ...   spore-print-color population habitat
0         p         x           s   ...                   k          s       u
1         e         x           s   ...                   n          n       g
2         e         b           s   ...                   n          n       m
3         p         x           y   ...                   k          s       u
4         e         x           s   ...                   n          a       g

[5 rows x 23 columns]
>>> df.dtypes
poisonous                   category
cap-shape                   category
cap-surface                 category
cap-color                   category
bruises                     category
odor                        category
gill-attachment             category
gill-spacing                category
gill-size                   category
gill-color                  category
stalk-shape                 category
stalk-root                  category
stalk-surface-above-ring    category
stalk-surface-below-ring    category
stalk-color-above-ring      category
stalk-color-below-ring      category
veil-type                   category
veil-color                  category
ring-number                 category
ring-type                   category
spore-print-color           category
population                  category
habitat                     category
dtype: object

pandas のカテゴリ型は Series#cat#categories で水準を確認できたりする。

>>> df.poisonous.cat.categories
Index(['e', 'p'], dtype='object')

整数値にラベルエンコードする

ただし、上記ではカテゴリ型の中身が 'object' になっているため、そのままでは LightGBM に食べさせることができない。 LightGBM では int, float, boolean しか扱うことができないため。 そこで scikit-learnLabelEncoder を使って対応する整数値に変換する。

>>> from sklearn import preprocessing
>>> for column in df.columns:
...     target_column = df[column]
...     le = preprocessing.LabelEncoder()
...     le.fit(target_column)
...     label_encoded_column = le.transform(target_column)
...     df[column] = pd.Series(label_encoded_column).astype('category')
... 

これで DataFrame がカテゴリ型のまま中身が整数になった。

>>> df.head()
   poisonous  cap-shape   ...     population  habitat
0          1          5   ...              3        5
1          0          5   ...              2        1
2          0          0   ...              2        3
3          1          5   ...              3        5
4          0          5   ...              0        1

[5 rows x 23 columns]

あとは普通に LightGBM で学習して汎化性能を検証するだけ。

まずはデータセットを説明変数と目的変数に分割する

>>> X, y = df.drop('poisonous', axis=1), df.poisonous

続いてデータセットを学習用、ハイパーパラメータ調整用、ホールドアウト検証用の 3 つに分割する。

>>> from sklearn.model_selection import train_test_split
>>> # データセットを学習用と検証用に分割する
... X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.7, random_state=42)
>>> # 検証用データセットをハイパーパラメータ調整用とバリデーション用に分割する
... X_eval, X_test, y_eval, y_test = train_test_split(X_val, y_val, test_size=0.5, random_state=42)

LightGBM を学習させる

続いて、分割したデータを元に LightGBM を学習させる。 学習に使う元ネタが pandas の DataFrame で、適切にカテゴリ型になっていれば LightGBM はそれを識別してくれる。

>>> import lightgbm as lgb
>>> lgb_train = lgb.Dataset(X_train, y_train)
>>> lgb_eval = lgb.Dataset(X_eval, y_eval, reference=lgb_train)
>>> # 学習用パラメータ
... lgbm_params = {
...     # 二値分類問題
...     'objective': 'binary',
...     # 評価方法
...     'metric': 'binary_error',
... }
>>> # 学習
... model = lgb.train(lgbm_params, lgb_train, valid_sets=lgb_eval)

ホールドアウト検証を用いた性能評価

学習させたモデルをホールドアウト検証用のデータセットで評価する。

>>> # バリデーション用データセットで汎化性能を確認する
... y_pred_proba = model.predict(X_test, num_iteration=model.best_iteration)
>>> import numpy as np
>>> # しきい値 0.5 で最尤なクラスに分類する
... y_pred = np.where(y_pred_proba > 0.5, 1, 0)

今回、精度で評価したところ全てのデータを正しく判別できた。 分割方法などにも依存するものの、このデータセットとモデルであればほぼ 100% に近い精度が得られるようだった。

>>> # 精度を確認する
... from sklearn.metrics import accuracy_score
>>> accuracy_score(y_test, y_pred)
1.0

また、カテゴリ型が指定されていない pandas の DataFrame や、それ以外をデータセットの元ネタにするときは、次のように categorical_feature を指定すると良い。

>>> lgb_train = lgb.Dataset(X_train, y_train, categorical_feature=list(X.columns))

いじょう。

Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理

Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理