突然ですが皆さんTDA（トポロジカルデータ解析）はやっていますか？
本記事はTDAの内容には触れずに，persistent homologyをうまく扱えるツールである HomCloud のチュートリアルを１つやってみたの記事です
＊HomCloudのチュートリアル記事自体は既にあるのですが，理解の備忘録的な意味合いです
＊persistent homology に関しては，今年のアドカレとかでガッツリ書きます
今は書かないといけない・回収しないといけないネタが多すぎる・・・（自分のせい）
インストール
自分の環境に合わせて記載されているコマンドをぽちぽちするだけ
筆者の環境は「OS X (Apple Silicon Mac)」で，特に問題なくできました
チュートリアル
初めはポイントクラウドから始めると良いらしいです
ちなみに筆者は jupyter notebook で server error が起こってなんで？って思っていましたが
- homcloud/tutorial/python-interface/pointcloud で jupyter notebook を起動
- jupyter 上で pointcloud.ipynb を選択
- トークンを求められたため，別ターミナル上で jupyter notebook list からトークンを取得し，それを jupyter 上に入力
で解決しました
参照
＊そもそもその前に仮想環境に入らずに jupyter を起動してしまって詰まっただなんて言えない・・・（雑魚）
3次元点集合データ(ポイントクラウド)
以下は理解の備忘録です（最低限これだけ理解しておけばのコマンドまとめでもある）
カレントディレクトリ：homcloud/tutorial/python-interface/pointcloud
前提
import homcloud.interface as hc # HomCloudのインターフェース
import numpy as np # NumPy (数値配列)
import matplotlib.pyplot as plt # Matplotlib (PDの可視化用)
%matplotlib inline
import pyvista as pv # PyVista (3D可視化)
step1:データからパーシステント図の計算
pointcloud.txt : ３次元空間上の点が1000個記載されているtxtファイル（入力データ）
pointcloud = np.loadtxt("pointcloud.txt")
hc.PDList.from_alpha_filtration(pointcloud, save_to="pointcloud.pdgm", save_boundary_map=True)
で pointcloud.txt のパーシステント図の情報が ./pointcloud.pdgm として出力されます
step2:パーシステント図の可視化
pdlist = hc.PDList("pointcloud.pdgm")
pd1 = pdlist.dth_diagram(1)
pd1.histogram().plot()
で ./pointcloud.pdgm から1次のパーシステント図のみをヒストグラムにプロット
＊色付けや log スケール・拡大・解像度の変更・pngファイルへの保存が可能
step3:ライフタイムの可視化
plt.hist(pd1.deaths - pd1.births, bins=100);
参考
pd1.births # 各点の誕生時間
pd1.deaths # 各点の消滅時間
pd1.deaths - pd1.births # 各点の生存時間
step4:逆解析
(0.0025, 0.008) 付近にある birth-death pair の stable volume を調べる
pair = pd1.nearest_pair_to(0.0025, 0.008)
stable_volume = pair.stable_volume(0.00005)
stable_volume の可視化
pl = pv.Plotter()
pl.add_mesh(stable_volume.to_pyvista_boundary_mesh())
pl.show()
stable_volume とポイントクラウドの重ね合わせ
pl = pv.Plotter()
pl.add_mesh(pv.PointSet(pointcloud))
pl.add_mesh(stable_volume.to_pyvista_boundary_mesh(), color="red", line_width=4)
pl.show()
参考
pair # 誕生時間と消滅時間のペア
pair.lifetime() # 生存時間
まとめとか感想とか
珍しくツールのチュートリアル記事を書きました
エンジニアっぽい記事をやっとQiitaで書けたなって感じです()
チュートリアルはまだまだ続いて，白黒画像とか機械学習と組合せたり，最終的には材料科学へ入っていくため，今後が楽しみです
今回の内容はここまでです．ここまでご覧になってくださった方々ありがとうございます！