機械学習でなんとなく材料研究者の気分を味わおう

[Nextremer Advent Calendar 2016](http://qiita.com/advent-calendar/2016/nextremer)の10日目の記事です。
自分は大学で**機械学習と材料開発**を組み合わせた分野を研究しています。
今回は**機械学習が大好きな人に材料研究者の気分を味わってもらいたい**と思い記事を書きました。
（本当は、材料研究者の人にAI研究者の気分を味わってもらってもっと材料に機械学習を取り入れてほしいです）

#材料開発って何だ
まず、材料開発って何やねん、何の材料やねんって感じだと思います。
おっしゃる通り、材料と言っても、セラミックス、高分子、金属など様々なものがあります。

例えばiPhone
![iphone.jpg](https://qiita-image-store.s3.amazonaws.com/0/146565/83e84272-b587-6250-e20f-bd0630833542.jpeg)

この中には、こんな小さなセラミックコンデンサが数百個入っています。

![0206_img0002.jpg](https://qiita-image-store.s3.amazonaws.com/0/146565/4a87dd9d-6884-c762-9990-6ba4d91706de.jpeg)

そして、Appleに選ばれるほどの高性能なコンデンサを作るためには、

・どんな元素を組み合わせよう？
・どんなプロセスで作ろう？

という難しい問題を解くことになります。解き方の一例を示します。

1. とりあえず過去の偉人が作ってたやつちょっと変えて作ってみるか
2. 意外とできない・・・もっと混ぜ時間増やしたほうがいいのかな（プロセス最適化）
3. できた！色々測定したり観察したりしよ（測定）
4. なるほど、こんな物理現象で高性能でてるんや（理論、解析）
5. 理論と解析で傾向がわかったし、これなら次はこんな元素で作ってみるか・・・（探索）

という感じです。これはあくまで一例ですが、材料研究とは

>プロセス最適化、測定、理論、解析、探索などをフル活用して、すごい材料を作ること

めちゃくちゃ雑ですが、こんな感じです。

#機械学習といえばデータ

機械学習でどうやって材料を開発するかというと、もちろんデータから特徴量を学習します。そして、良い物性を持つ材料を予測して探索します。

しかし、データが全然足りないため材料研究においてAIはあまり成果を出すことができていないのが現状です。

材料におけるデータは、大きく分けて2つあります。
・実験による測定データ
・計算による解析データ

####実験による測定データ
実験による測定データは、

- 論文
- 実験ノートに書かれたまま埃をかぶる
- 全世界の研究室や企業に分散

という状況です。
論文に載っていない失敗のデータとか、DeepLearningとかやるんだったらめちゃめちゃ欲しいんですが実験ノートに書かれたまま埃をかぶっています。
巨大データベースが欲しいところですが、みんな大事な成果を隠したがっていますから、難しいみたいです。

####計算による解析データ
計算データのデータベースを無料で公開している[Materials Project](https://materialsproject.org/
)という団体があります。
ここが公開してくれているデータは、**第一原理計算**という怖い響きの手法で計算されています。データは、https://materialsproject.org/ を使ってREST形式で簡単に取得することができます！

気になる人向けに。
>第一原理計算とは？
第一原理に基づいて行われる計算（手法）の総称。（[Wikipedia](https://ja.wikipedia.org/wiki/%E7%AC%AC%E4%B8%80%E5%8E%9F%E7%90%86%E8%A8%88%E7%AE%97)より）

はい。

>第一原理とは？
自然科学での第一原理（first principles）は、近似や経験的なパラメータ等を含まない最も根本となる基本法則をさし、そのことを前提にすると自然現象を説明することができるものである。（[Wikipedia](https://ja.wikipedia.org/wiki/%E7%AC%AC%E4%B8%80%E5%8E%9F%E7%90%86#.E8.87.AA.E7.84.B6.E7.A7.91.E5.AD.A6.E3.81.AB.E3.81.8A.E3.81.91.E3.82.8B.E7.AC.AC.E4.B8.80.E5.8E.9F.E7.90.86)より）

簡単に言うと、電子とか原子とか基本的な物理法則だけで、めっちゃ複雑な物質の振る舞いを計算しちゃおうってやつです。
でも、
>いわゆる第一原理による電子状態計算手法によって扱える原子の数は2003年現在でも100～1000個程度までであり、アボガドロ定数に遠く及ばない。1000原子のオーダーでようやく最も簡単な構造のたんぱく質（或いはアミノ酸）が扱えるかもしれないというレベルである。（[Wikipedia](https://ja.wikipedia.org/wiki/%E7%AC%AC%E4%B8%80%E5%8E%9F%E7%90%86%E8%A8%88%E7%AE%97)より）

らしいです。2016年の今はめちゃくちゃ進歩してますが、複雑な材料や物性を計算しようと思うと**大きなコストやノウハウが必要**なのが現状です。なので、大量のデータを取るならガチで計算機回しまくって、地味に計算するか、Materials Projectから撮ってきましょう。

#難しいんだろうけど2016年だし機械学習で材料を作ってみたいよね
やっぱり時代はAI、機械学習です。AIで材料作るには、

**①化学式からAIで材料の物性を予測→②AIで膨大な原子の組み合わせから探索→③AIでプロセス最適化しながらロボットとかも制御しながら材料作る**

みたいな流れが必要で、現状は**①すら全然できていません。**
ということで、やっと本題で申し訳ありませんが、**物性予測**をみんなでやってみて、**機械学習で材料研究者の気分**を味わいましょう！

#実際に機械学習で材料物性を予測してみよう
実際に、AIで材料開発を助けるサービスを手がけている[Citrine Informatics](https://www.citrine.io/)という会社があり、いい感じのtutorialを乗せてくれていたので、やってみたいと思います。
ソース：[MACHINE LEARNING FOR THE MATERIALS SCIENTIST](https://www.citrine.io/blog/2015/3/3/machine-learning-mat-sci-1)

**目的**

物質の**バンドギャップ**を機械学習します。
バンドギャップとは電子が移動する時の障壁の大きさを表したものです。バンドギャップが広いと電気が通りません（石）。バンドギャップが狭いと電気が通ります（金属）。まあまあだと通したり通さなかったりします。（半導体）

**必要なもの**

- python3.5（2系でも多分可ですが試してません）
- scikitlearn
- numpy
- pymatgen（[Materials Project](https://materialsproject.org/
)のデータ取得や、物性計算に便利なモジュール）

**インストール**

```bash
pip install scikit-learn
pip install numpy
pip install pymatgen
```

**化学式だけの説明変数で機械学習してみる**

まずデータを取得します。
[ここ](https://citrineinformatics.box.com/shared/static/0ch2f96jxbqtntia49ipk7g8vx64tuap.csv)からbandgapDFT.csvというファイルをダウンロードします。
ちなみに、DFTとは、Density Functional Theoryという第一原理計算の手法です。

では、まずは必要なライブラリをインポートして、ダウンロードしたcsvを読み込みましょう。

```py3:bandgap.py
from pymatgen import Composition, Element
from numpy import zeros, mean
trainFile = open("bandgapDFT.csv","r").readlines()
```
データを見てみると、

```
LiH,2.981
BeH2,5.326
B9H11,2.9118
B2H5,6.3448
BH3,5.3234
B5H7,3.5551
H34C19,5.4526
H3N,4.3287
(以下略)
```
このようになっていると思います。一行目が化学式、二行目がバンドギャップ(eV)です。化学式を機械学習のための固定長ベクトルとして取り扱うための関数を作ります。ここで出てくるcompositionは、pymatgenのCompositionオブジェクトで、化学式から原子や構成比などを取得できます。

```py3:bandgap.py
#input:pymatgenのCompositionオブジェクト
#output:組成ベクトル
def naiveVectorize(composition):
       vector = zeros((MAX_Z))
       for element in composition:
               #elementは原子。fractionはその原子が組成に含まれる割合
               fraction = composition.get_atomic_fraction(element)
               vector[element.Z - 1] = fraction
       return(vector)
```
csvから化学式とバンドギャップを読み込み、上述した関数で化学式から説明変数のベクトルを用意します。（教師データの生成）

```py3:bandgap.py
materials = []
bandgaps = []
naiveFeatures = []

MAX_Z = 100 #特徴量ベクトル最大長さ

for line in trainFile:
       split = str.split(line, ',')
       material = Composition(split[0])
       materials.append(material) #化学式
       naiveFeatures.append(naiveVectorize(material)) #特徴量ベクトル生成
       bandgaps.append(float(split[1])) #バンドギャップの読み込み
```

特徴量ベクトルは、こんな感じの超離散データになります。

```
[ 0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,
        0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,
        0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,
        0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,
        0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,
        0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,
        0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,
        0. ,  0. ,  0. ,  0. ,  0. ,  0.5,  0. ,  0. ,  0. ,  0. ,  0. ,
        0. ,  0. ,  0. ,  0. ,  0. ,  0.5,  0. ,  0. ,  0. ,  0. ,  0. ,
        0. ]
```

ということで、教師データは揃いました。あとは、scikitlearnの便利なmoduleを使って機械学習していきます。その前に、まずはバンドギャップの予測精度の検証用に、単純な平均誤差を求めておきます。

```py3:bandgap.py
baselineError = mean(abs(mean(bandgaps) - bandgaps))
print("Mean Absolute Error : " + str(round(baselineError, 3)) + " eV")

```
ここで、 0.728eVになるはずです。ということで、ランダムフォレストでバンドギャップ予測器を作ってみましょう。

```py3:bandgap.py
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import ShuffleSplit
from sklearn import linear_model, metrics, ensemble

#sklearnのランダムフォレスト回帰
rfr = ensemble.RandomForestRegressor(n_estimators=10)

#交差検証します
cv = ShuffleSplit(n_splits=10, test_size=0.1, random_state=0)

scores_composition = cross_val_score(rfr, naiveFeatures,\
	bandgaps, cv=cv, scoring='neg_mean_absolute_error')

print("Mean Absolute Error by Random Forest with composition data: "\
	+ str(round(abs(mean(scores_composition)), 3)) + " eV")
```
平均誤差がだいたい0.36〜eVになっていれば、これで完成です。0.728eVに比べて、なんとか学習できているとは思います。

**物理的な値を説明変数にしてみよう**
これまでは、どんな原子が入っているか、だけを説明変数として用いてきましたが、次は、原子の組成比、電気陰性度などの値を説明変数として用います。たったの4次元です。
例えばBeH2だと、
[2.0, 0.63, 1.0, 2.0]
となります。

```py3:bandgap.py
physicalFeatures = []

for material in materials:
       theseFeatures = []
       fraction = []
       atomicNo = []
       eneg = []
       group = []

       for element in material:
               fraction.append(material.get_atomic_fraction(element))
               atomicNo.append(float(element.Z))
               eneg.append(element.X)
               group.append(float(element.group))

       mustReverse = False
       if fraction[1] > fraction[0]:
               mustReverse = True

       for features in [fraction, atomicNo, eneg, group]:
               if mustReverse:
                       features.reverse()
       theseFeatures.append(fraction[0] / fraction[1])
       theseFeatures.append(eneg[0] - eneg[1])
       theseFeatures.append(group[0])
       theseFeatures.append(group[1])
       physicalFeatures.append(theseFeatures)
```
作った説明変数physicalFeaturesを用いて、同様に決定木で学習させてみましょう。

```py3:bandgap.py

scores_physical =cross_val_score(rfr, physicalFeatures,\
	bandgaps, cv=cv, scoring='mean_absolute_error')

print("Mean Absolute Error by Random Forest with physical data: "\
	+ str(round(abs(mean(scores_physical)), 3)) + " eV")
```
すると
Mean Absolute Error by Random Forest with physical data:: 0.267 eV
となります。

比較すると

- ランダムフォレスト with 化学式・・・0.362eV
- ランダムフォレスト with 物理量・・・0.267eV

と、物理量を説明変数として用いたほうが精度が良いことがわかりました。
全然使い物はなりませんが、なんとなく材料研究者の気分を機械学習で味わうことができたのではないでしょうか？

#最後に

実際は、ニューラルネットや決定木を用いた物性予測や、遺伝アルゴリズムを用いた材料探索手法を提案している論文はすでに多く存在していて、今は実用レベルで応用できてはいませんが、将来確実にAIを用いた材料研究は盛んになっていくと思います。しかし、日本はアメリカなどに比べてこの分野では大変で遅れている感があります。材料、物性を専門にしている方がいましたら、是非とも機械学習を取り入れた研究も進めてほしいなと思います。

ありがとうございました！