A Robot’s Sigh
Pythonで準備してC++で計算する自明並列ワークフロー
はじめに
例えば相図を調べる場合など、数値計算では同じ様な計算を大量に実行することがよくあります。この時、様々なパラメータについて計算を実行することになりますが、それぞれの計算は独立に実行できます。したがって、100点のパラメータがあるなら100個の計算を並列に実行できます。このような並列計算を自明並列と呼びます。自明並列は数値計算において非常に重要なタスクですが、実際に実行しようとするといろいろ面倒だったり、どうやってファイルを管理するするのが良いかわからなくなったります。
そこで、以下では
- 調べたいパラメータが複数あり、それぞれについて独立に実行できる
- さらに、それぞれのパラメータについて乱数シードを変えて複数回計算して、後で統計平均をとる
というタスクについて、パーコレーションを題材に
- Pythonスクリプトでインプットファイルを大量に作成し
- それをC++の計算エンジンが自明並列で実行し
- 計算エンジンが吐いた大量のファイルをPythonスクリプトで解析する
というワークフローのやり方を紹介しようと思います。
サンプルリポジトリは
kaityo256/pycpp-parallel-pipeline
に置いておきます。
ワークフローの概要
想定するワークフローは以下の通りです。
何をやっているかの解説は後回しにして、まずはワークフローをどのように実行するか確認してみましょう。以下は物性研スパコンでの実行を想定します。
リポジトリの準備
スパコンにログインしたら、まずはリポジトリをクローンします。Git Submoduleがあるため、クローンする時に--recursiveの指定が必要です。
git clone --recursive https://github.com/kaityo256/pycpp-parallel-pipeline.git
cd pycpp-parallel-pipeline
PythonがNumPyやPyYAMLを使うため、uvもしくはvenvにより仮想環境の構築を行います。
uvを使う場合。
uv venv
source .venv/bin/activate
uv pip install numpy pyyaml
venvを使う場合。
python3 -m venv .venv
source .venv/bin/activate
python3 -m pip install --upgrade pip
python3 -m pip install numpy pyyaml
これで実行準備ができました。
インプットファイルの作成
まず、自明並列のためのインプットファイルを作成します。2次元ボンドパーコレーションのパーコレーション確率を計算することにしましょう。各タスクは、システムサイズ、ボンドの占有確率、シード、サンプル数などを与える必要があります。これらを一気に作るため、Pythonスクリプトにはシステムサイズや、サンプル数に加えて、パラメータを調べる範囲、パラメータを何点とるかを教える必要があります。
この、「インプットファイル作成のためのインプットファイル」をinput.yamlとしましょう。例えばこんな感じです。
L: 128
start: 0.45
end: 0.55
num_points: 100
num_samples: 100
num_seeds: 10
output_dir: output
重要なのはstart、end、num_pointsです。これにより、ボンド占有確率を、startからendまでnum_points分割して生成します。
実行してみましょう。
python3 generate_inputs.py input.yaml
すると、outputディレクトリに200個のYAMLファイルが作成されます。
output/L128_p4500_s000.yaml
output/L128_p4500_s001.yaml
...
output/L128_p5500_s008.yaml
output/L128_p5500_s009.yaml
さらに、cpsにわたすためのタスクリストtask.shが作成されます。task.shの中身はこんな感じです。
./cpp/percolation output/L128_p4500_s000.yaml
./cpp/percolation output/L128_p4500_s001.yaml
(snip)
./cpp/percolation output/L128_p5500_s008.yaml
./cpp/percolation output/L128_p5500_s009.yaml
要するにそのままシェルスクリプトとして逐次実行可能なリストになっています。
並列実行
並列実行にはcpsを使うため、予めビルドしておきます。
cd cps
make
cd ..
また、計算エンジンであるpercolationもビルドしておきましょう。
cd cpp
make
cd ..
cps/cpsとcpp/percolationが作成されたら、以下のジョブスクリプトを投入します。物性研スパコンシステムBがSlurmを利用しているため、Slurmのジョブスクリプトになっています。なお、ここではテストのためにデバッグキューであるi8cpuキューを指定していますが、動作確認後の実際の計算はF1cpuなどの計算実行用キューを使うようにしてください。
#!/bin/sh
#SBATCH -p i8cpu
#SBATCH -N 1
#SBATCH -n 128
srun ./cps/cps task.sh
cpsにtask.shを渡しているだけです。実行すると、cpsがタスクリストを上から順番に「空いているプロセス」に渡していきます。プロセスが全てBusyになったら、どれかが空くまで待ち、先に渡したタスクが終わったプロセスに次のタスクを渡していきます。
実行終了すると、outputディレクトリに200個のdatファイルが作成されています。
L128_p4500_s000.dat
L128_p4500_s001.dat
(snip)
L128_p5500_s008.dat
L128_p5500_s009.dat
解析
出力された200個のデータを解析しましょう。解析スクリプトはanalyze_results.pyです。このファイルに「どんな実行結果ファイルが生成されているか」を教えるために、入力ファイル生成スクリプトgenerate_inputs.pyに渡したのと同じYAMLファイルを渡します。
python3 analyze_results.py input.yaml
実行結果は以下の通りです。
output/L128_p4500_s001.dat
output/L128_p4500_s002.dat
(snip)
output/L128_p5500_s008.dat
output/L128_p5500_s009.dat
Saved analysis result to output/L128.dat
作成されるファイルoutput/L128.datはこんな感じです。
# bond_probability percolation_probability std_error
0.4500000000 0.0000000000 0.0000000000
0.4552631579 0.0000000000 0.0000000000
0.4605263158 0.0000000000 0.0000000000
0.4657894737 0.0000000000 0.0000000000
0.4710526316 0.0020000000 0.0013333333
0.4763157895 0.0120000000 0.0024944383
0.4815789474 0.0400000000 0.0073029674
0.4868421053 0.1100000000 0.0122020035
0.4921052632 0.2400000000 0.0141421356
0.4973684211 0.3920000000 0.0167862377
0.5026315789 0.5850000000 0.0192786583
0.5078947368 0.7590000000 0.0145640348
0.5131578947 0.8780000000 0.0100884973
0.5184210526 0.9430000000 0.0085699734
0.5236842105 0.9790000000 0.0056666667
0.5289473684 0.9940000000 0.0030550505
0.5342105263 1.0000000000 0.0000000000
0.5394736842 1.0000000000 0.0000000000
0.5447368421 1.0000000000 0.0000000000
0.5500000000 1.0000000000 0.0000000000
各ボンド占有確率について10個ずつ異なる乱数シードの結果を読み込み、それぞれについて平均とエラーバーを出力します。
リポジトリにあるプロットファイルplot.pltを使って、
gnuplot plot.plt
を実行すると、以下のようなファイルが作成されます。
インプットファイル一つから、パラメータの生成、計算の実行、解析までの一通りのワークフローが実行できました。
コードの開発
では、実際に自分が必要なワークフローを実行するためにどのようにすればよいかを見てみましょう。
開発する必要があるファイルは以下の3つです。
generate_inputs.py計算エンジンに渡すインプットファイルとタスクリストを生成するPythonスクリプト- 計算実行エンジン(言語はなんでも良いが、ここではC++を想定)
analyze_results.py出力結果を解析するPythonスクリプト
インプットファイル生成スクリプト
計算エンジンが必要とするパラメータを含むYAMLファイルを生成するスクリプトです。コードは難しくありませんが、パラメータの範囲やサンプル数などをスクリプトの中に書いてしまう(ハードコーディング)と、毎回Gitが修正を検知してしまって鬱陶しいことになります。スクリプトやコード類は開発が終わったら触らず、パラメータ類は外に出すべきです。
また、
python3 generate_inputs.py 128 100 10
などのように実行時引数でパラメータを渡すのもあまり良くありません。このようにして実行すると、「このグラフを作成した時に投げたパラメータなんだっけ？」と、後で必要になった時に思い出せずに苦労することがあります。パラメータはファイルから渡すのが原則です。
パラメータファイルのフォーマットはなんでも良いですが、広く使われているJSONかYAMLが良いでしょう。ここではYAMLを採用します。PythonでYAMLを読むためにはPyYAMLが必要です。デフォルトでは入っていないため、pip install pyyamlを実行する必要があります。uvかvenvで仮想環境を作ってインストールするとよいでしょう。
YAMLファイルはyaml.safe_loadで読み込むことができます。
import yaml
if __name__ == "__main__":
if(len(sys.argv)!=2):
print("Usage: python3 generate_inputs.py input.yaml")
exit(1)
config_file = sys.argv[1]
with open(config_file) as file:
config = yaml.safe_load(file)
configは辞書になっており、例えば
L: 128
start: 0.45
end: 0.55
num_points: 20
num_samples: 100
num_seeds: 10
output_dir: output
のようなファイルを読み込んだら、
L = int(config["L"])
などとして読み込めます。デフォルトで文字列になっているため、intやfloatなどで変換してやる必要があります。これらが読み込めたらインプットファイルを大量生成したりタスクリストを作ったりするのは難しくないでしょう。
なお、ファイルの出力先はoutputやresultsなど、別ディレクトリに分けた方が良いです。計算によっては出力ディレクトリをわけたくなることもあるため、output_dirを指定できるようにしておくと便利です。
計算実行エンジン
計算エンジンは、Pythonスクリプトが生成したインプットファイルを読み込んで、結果をファイルに書き込めればOKです。そのフォーマットはYAMLやTOML、その他いろいろ考えられますが、ここではPythonに合わせてYAMLにしておきましょう。
C++でYAMLを読むためのライブラリもいろいろありますが、ここでは軽量なシングルヘッダーライブラリparamを使うことにします。
paramはYAMLライクなパラメータを読み込むためのヘッダーで、
# Comment
bool = yes
int = 12345
double : 12.345
str : "hello"
のように、#から始まる行がコメント、=もしくは:をデリミタとして
key : value
key = value
のように記述されたファイルを読み込みます。
コンパイルさえできれば、ファイルをそのままリポジトリに配置しても、Git Submoduleを使っても、どちらでも構いません。Git Submoduleを使う場合は、コンパイル時に-Iを指定して場所を教え、かつVSCodeなどの開発環境にインクルードファイルの場所を教えてやる必要があります。
使い方は
param::parameter params("input.cfg");
と、ファイル名をコンストラクタに渡すだけです。多くの場合、
const std::string filename = argv[1];
param::parameter param(filename);
と、argv経由で受け取ることになるでしょう。
param::parameterのインスタンスを作ったら、
bool b = params.get<bool>("bool");
int i = params.get<int>("int");
double d = params.get<double>("double");
std::string s = params.get<std::string>("str");
と、PyYAMLと似たような形でパラメータを受け取ることができます。
受け取ったパラメータから適切に出力ファイル名を作成して、そこに計算結果を出力します。
例えば
L: 128
bond_probability: 0.45
seed: 0
num_samples: 100
というファイルを受け取ったら、L128_p4500_s000.datに結果を出力する、といった具合です。
解析スクリプト
解析スクリプトの作り方はインプットファイル作成スクリプトとほとんど同じです。インプットファイルに渡したYAMLファイルと同じファイルを渡すことで、「こんなファイルが作成されたはず」と推定して、そのファイルを読み込みに行きます。
出力ディレクトリをglobを使って探索し、出力されている*.datファイルから結果を作る方法は事故の可能性があるので避けたほうが無難です。
その他のTipsなど
インプットファイル名
これは個人的な経験ですが、インプットファイル作成スクリプトに渡す「ワークフロー全体を制御するファイル」のファイル名は、L128_numseeds_100.yamlのように「そのワークフローの意味」がわかるようにするよりも、20260715.yamlのように、YYYYMMDD形式の日付で管理したほうが良いことが多いです。
スパコンにジョブを投げる際、研究ノートに
## 2026/7/15
サイズによる緩和を調べるためのジョブ。
L=128とL=256で緩和時間がどれくらい変わるか確認。
JOBID: 2982722
のように、日付とジョブの投入意図、JOB IDを残すようにしています。そうすると、20260715.yamlが残っている時、それをどんな意図で投げたかが後から参照できて便利です。
また、研究ノートにはJOBIDも残しておくとたまに便利です。例えばSlurmなら標準出力/標準エラー出力がslurm-2982722.outのようにジョブIDがついたファイルに保存されます。「たまにジョブに失敗する」みたいな時、失敗するジョブの共通点を調べるなどの目的でJOBIDが必要になることがあります。例えばスパコンの通信ライブラリが原因で失敗しているっぽい時など、JOBIDを報告できるようになっていると調査がスムーズです。
少数の長いジョブ vs. 多数の短いジョブ
パーコレーションなどのようなモンテカルロ法では、統計平均をとるためにある程度のサンプル数が必要です。この時、もし可能であれば長時間のジョブを流すよりも、短いジョブを多数流す方が楽です。
今回のワークフローでも、100サンプルをとるジョブを乱数シードを変えて10個のタスクとしましたが、これは1つのタスクで1000サンプルとるのと等価です。しかし、スパコンには「ジョブの実行時間制限」があるため、ある程度以上長い計算ができません。一方、多くの場合において、スパコンで利用可能な計算資源の同時並列実行数は、「普段使い」では使い切れないことが多いです。
そこで、短いタスクに分割して投げることで、time to solutionを短くできます。さらに、パラメータごとに実行時間がバラける場合、cpsは「タスクが終わったプロセスに次のタスクを回す」という単純スケジューリングをすることから、タスクが多数あればあるほど結果としてロードバランスが取れるようになります。
その分、出力されるファイルが増えますが、解析スクリプトさえ作ってしまえばさほど面倒は増えません。
まとめ
Pythonスクリプトで大量のインプットファイルを作成し、C++で計算を自明並列で実行し、その結果をPythonで解析する、一連のワークフローのやり方を説明しました。自明並列はスパコンを使う第一歩です。研究室クラスタでは1ヶ月以上かかるような計算が1日で終わるようになると、研究の「見え方」が変わります。まずは自明並列ができるようになって、その後、非自明並列へと挑戦すると良いでしょう。
同様なワークフローで計算から論文出版までの記録をまとめたコードを書いてから論文が出版されるまでも参考になるかもしれません。