[Python] 材料開発におけるパーシステントホモロジーと機械学習

# 1. 材料開発とパーシステントホモロジー
材料の分野では、その特性が構造に大きく依存する系が多くあります(e.g. 機械特性、熱特性、電気特性、etc.)。このような場合、一般には構造の代表値(e.g. 平均粒径、含有率、配位数、etc.)と特性の相関を扱うことが多いと思います。

トポロジーに着目してデータのおおまかなカタチを解析する手法を、位相的データ解析(Topological Data Analysis, TDA)と言います。TDAの中でも、データのトポロジカルな特徴を抽出できる強力なツールとして、パーシステントホモロジーが期待されています。

材料組織の特徴量を抽出する手法として、畳み込みニューラルネットワークを用いる場合もあります。このような手法と比較して、パーシステントホモロジーは構造の特徴をマルチスケールに抽出できること、回転にロバストな点などにおいて優位な気がします。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/676917/fd3c7667-35fc-e7a1-1ce6-46975d8fd9df.png)
パーシステントホモロジーの模式図
Obayashi Ippei, et al., Journal of Applied and Computational Topology 1 (2018): 421-449. より

概念としては近いですが、二次元(SEM像、光学顕微鏡像、AFM像、etc.)や三次元(X線CT像、etc.)の二値化像ではこのようにデータを扱います。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/676917/6d5098d1-b313-3bc9-42bc-2ea645f11187.png)
バイナリーイメージのフィルトレーションの模式図
Obayashi Ippei, et al., Journal of Applied and Computational Topology 1 (2018): 421-449. より

(a)にオリジナル画像を、(b)にマンハッタン距離を示しています。例えば、0次のパーシステントダイアグラムでは、データの「つながり」を考えます。(birth, death)がそれぞれ(-2,-1)、(-2,3)のときに島のbirthとdeathが生じています。ただし、マンハッタン距離が-3のとき、青の四角で表した島は除外しています。

今回はパーシステントホモロジーで粒子の構造がどのように可視化できるか試してみます。

# 2. Homcloud
パーシステントホモロジーを実行するために、[HomCloud](https://homcloud.dev/index.html)を利用しました。解説やチュートリアルが豊富で、素人には大変助かりました...

より詳細には、以下の論文などを参考にしてください。

https://journals.jps.jp/doi/full/10.7566/JPSJ.91.091013

Obayashi, Ippei, Takenobu Nakamura, and Yasuaki Hiraoka, journal of the physical society of japan 91.9 (2022): 091013.

# 3. データセット
下記サイトから、Ottawa砂とAngular砂を圧縮させた際のX線CT像のデータを使います。

https://data.mendeley.com/datasets/szn4jtfkbx/1

CT像をgifにして確認してみましょう。
左から順に、0、2、6.1、10.2MPaで圧縮したものです。

Ottawa砂
![ottawa_0.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/676917/55195e56-31cf-7d84-0635-381a07062e77.gif)![ottawa_2.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/676917/240a2cd8-1b99-66f3-19a7-edc34a048bae.gif)![ottawa_6p1.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/676917/f51fa3c9-6efb-0832-65ee-7dcc3f463933.gif)![ottawa_10p2.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/676917/fb6e7a57-8630-fcbc-a1c8-6c5ab16d3d5d.gif)

Angular砂
![angular_0.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/676917/c1411726-52b7-63da-0321-ba2355b86e30.gif)![angular_6p1.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/676917/cc0f4dba-c1c2-5c0c-45ac-dc6a2784c265.gif)![angular_2.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/676917/95989d12-9bbc-e1b5-2359-bae9a74889a3.gif)![angular_10p2.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/676917/2c14ff03-2dfd-98ed-53f2-591fd5c97032.gif)

Ottawa砂とAngular砂は粒の形状から簡単に分離できそうです。
同種の砂の場合、圧力による構造の違いはわかりません。

このデータを適当な閾値で二値化して以降の処理を実行します。

# 4. 構造のベクトル化と主成分分析
## 4.1. パーシステントダイアグラム
1次のパーシステントダイアグラムを見てみます。
![all_pd.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/676917/8db8d45d-fe17-eb53-fc44-62b71ce9493f.png)

birthとdeathのヒストグラムを確認できました。今回は1次のパーシステントダイアグラムなので、なんらかのリング構造を見ているはずです。

## 4.2. 主成分分析
[HomCloud](https://homcloud.dev/index.html)のチュートリアルにもありますが、パーシステントダイアグラムのベクトル化には[こちらの文献](https://arxiv.org/abs/1706.10082)を参考にしています。

https://arxiv.org/abs/1706.10082

拙い理解ですが、ガウシアンを用いてパーシステントダイアグラムの「ぼやけた」行列を得る->平坦化してベクトル化...という操作になっていると思います。パーシステントダイアグラムの小さな変化に対してロバストであり、機械学習タスクに活用することができます。

このようにして得たベクトルに対し主成分分析を行い、二次元空間に射影してみましょう。
![pca_result.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/676917/1e8ba3df-2d2e-67c8-0940-cb4ebb97184d.png)

●でOttawa砂を、◆でAngular砂をプロットしています。また、カラーバーで圧力を示しています。**パーシステントホモロジーによって構造の特徴を抽出することで、Ottawa砂とAngular砂を分類し、かつ異なる圧力における構造の違いも捉えられている**ようです。本来であれば、もう少しn数が多いと良かったのですが、今回はたかだか350x350x350のデータでしたので、分割などはせずにn=1で行っています。

実際の材料設計においては、今回の構造のベクトル化に加えて、特性や原料あるいはプロセス条件などを組み合わせることで、予測・回帰・最適化などの機械学習タスクが可能です。また、構造のベクトルからパーシステントダイアグラム、構造データまで逆解析を行うことで、実どのような構造が特性に影響を及ぼしているか可視化することも可能です。