コンテンツにスキップ

Rust

インストール

  • rustupのインストール
1
2
rustup component add rustfmt clippy rust-analyzer
cargo install evcxr evcxr_repl

1
rustup component add rls rust-analysis rust-src

アップデート

1
rustup update

アンインストール

1
rustup self uninstall

パッケージ(クレート)のインストール

  • Cargo.tomlに適当に記入する
  • cargo buildする

ビルドの確認

1
cargo check

フォーマッター

  • rustfmt

AtCoder用(?)inputマクロ

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
macro_rules! input {
    (source = $s:expr, $($r:tt)*) => {
        let mut iter = $s.split_whitespace();
        input_inner!{iter, $($r)*}
    };
    ($($r:tt)*) => {
        let mut s = {
            use std::io::Read;
            let mut s = String::new();
            std::io::stdin().read_to_string(&mut s).unwrap();
            s
        };
        let mut iter = s.split_whitespace();
        input_inner!{iter, $($r)*}
    };
}

macro_rules! input_inner {
    ($iter:expr) => {};
    ($iter:expr, ) => {};

    ($iter:expr, $var:ident : $t:tt $($r:tt)*) => {
        let $var = read_value!($iter, $t);
        input_inner!{$iter $($r)*}
    };
}

macro_rules! read_value {
    ($iter:expr, ( $($t:tt),* )) => {
        ( $(read_value!($iter, $t)),* )
    };

    ($iter:expr, [ $t:tt ; $len:expr ]) => {
        (0..$len).map(|_| read_value!($iter, $t)).collect::<Vec<_>>()
    };

    ($iter:expr, chars) => {
        read_value!($iter, String).chars().collect::<Vec<char>>()
    };

    ($iter:expr, usize1) => {
        read_value!($iter, usize) - 1
    };

    ($iter:expr, $t:ty) => {
        $iter.next().unwrap().parse::<$t>().expect("Parse error")
    };
}

AtCoderでの競プロ用cargo-compete

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
cargo install cargo-compete


cargo compete new tessoku-book

cargo compete test a01

cargo compete submit a01

cargo compete submit a01 --no-test

evcxrでクレートを読み込む

  • (ChatGPTに教えてもらった)
  • extern crateを使う.
  • 下記サンプルを参照.
1
2
3
4
5
6
7
8
9
extern crate rand;
use rand::Rng;

fn main() {
    let mut rng = rand::thread_rng();
    let x: u8 = rng.gen();
    println!("{}", x);
    println!("TEST");
}

proconioのfastoutを使うときの注意

  • 公式
  • Cargo.tomlの指定をproconio = { version = "=(version)", features = ["derive"] }にする

配列や構造体のprint

  • URL
  • 配列は次の通り
1
2
let list = vec![0; 10];
println!("{:?}", list);
  • 構造体は次の通り#[derive(Debug)]をつける.
1
2
3
4
5
6
7
8
9
#[derive(Debug)]
struct TestStruct {
  x: f64,
  y: f64,
  side: f64,
};

let t = TestStruct{x: 2.0, y:3.0, side: 5.0};
println!("{:?}", t);

引数を参照で渡す: &

1
2
let mut guess = String::new();
io::stdin().read_line(&mut guess)

パッケージ(クレート, crate)の追加

  • cargo add randなどなど

数値計算メモ

Rust・Gnuplot・ffmpeg で Python・matplotlib よりは速くなった: 1d CFD, バーガース方程式

結論: GitHub への参照

先日この記事で最低のゴミ・産業廃棄物と書いた, 肖鋒, 長崎孝夫, 『数値流体解析の基礎 - Visual C++とgnuplotによる圧縮性・非圧縮性流体解析』 2020/1/9 の試行錯誤の結果をまとめておく. はじめに書いておくと, GitHub にあげてあるので, コードを見たい人はそちらを見てほしい.

あれほどゴミのような内容で 3600 円などいい商売だと思いつつ, お金がない中で買ったそれなりに値段の張る本なので, 何だかんだで頑張って C のコードを書き直しながら読んでいる.

自分用メモとして, Python を捨てて Rust・Gnuplot・ffmpeg にいたる経緯はあとで書くとして, まずはコードの話をしよう.

png 作成: Gnuplot

次の crate を使った.

これを見てもらえれば全ては終わる. コード全体は GitHub を見てもらうことにして, ここでは該当部分を抜き出しておこう.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
// cargo-deps: chrono, gnuplot
extern crate gnuplot;
use gnuplot::*;

pub fn write_png(cnf: &Config, cdata: &CalcData) -> Result<(), Box<dyn error::Error>> {
    let file_name: String = format!("{}/img.{:08}.png", &cnf.dir_name, &cdata.n);
    let mut fg = Figure::new();
    fg.axes2d()
        .set_title(&cnf.title, &[])
        .set_legend(Graph(1.0), Graph(1.0), &[], &[])
        .set_y_range(Fix(cnf.graph_ylim_min), Fix(cnf.graph_ylim_max))
        .set_x_label("x", &[])
        .set_y_label("u", &[])
        .lines(
            &cdata.xs,
            &cdata.us,
            &[Caption("numerical"), LineWidth(2.5)],
        )
        .lines(&cdata.xs, &cdata.ues, &[Caption("exact"), LineWidth(2.5)]);
    fg.save_to_png(file_name, 800, 800);

    Ok(())
}

他にもいろいろインポートは必要だろうが省略した. エラー処理もすさまじくひどいが, よくわかっていない. とりあえず最低限を書けるようになるまで読んだ公式のチュートリアルで出てきた記述を適当にコピペで使っている. (これこそ私が批判した上掲書のゴミコードスタイルだが, 全て棚に上げている.)

いま 1 次元のバーガースを解いているので, 2 次元プロットしたいというだけ. 計算法の都合によって途中で最大値がどんどん落ちていき, 動画化したときに y 軸方向にぶれて見づらいので, y_range を設定してある.

動画化: ffmpeg

これも FFI があるようだが, かったるいので直接コマンド実行している.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
use std::process::Command;
fn write_mp4(cnf: &Config, cdata: &CalcData) -> Result<(), Box<dyn error::Error>> {
    let options = ["-r", "10", "-i", &cnf.specify_png, &cnf.movie_name];

    let run_ffmpeg = Command::new("ffmpeg")
        .args(&options)
        .output()
        .expect("failed to start `ffmpeg`");
    println!("{}", String::from_utf8_lossy(&run_ffmpeg.stdout));

    Ok(())
}

本質的なエラー処理をつけていないのに Result を使っていて頭が悪いが, IO 系なのでつけるだけつけた. もちろんコマンド実行でぐぐればすぐに出てくるし, 多少プログラムが書ける人ならすぐわかるレベルの話だろう.

経緯

改めて書いておくと, 世界を汚染する害悪こと肖鋒, 長崎孝夫, 『数値流体解析の基礎 - Visual C++とgnuplotによる圧縮性・非圧縮性流体解析』 2020/1/9 では, 計算を実行しながら gnuplot で逐次画像出力の形で「動画化」されていて, ファイル出力としては eps だった. postscript は何もわからないし調べるのも面倒, 拡張子 cpp で実際には C++ でさえなく, 本の指示通りに Visual C++ でやるとコンパイルが通らない C 言語のコードのデバッグはやりたくなかった. (この C コードはネットから落とせるのだが, 本の購入者しか開けないようにパスワードがかかっているので, 買わないと地獄は感じ取れない.)

これまでいろいろな都合によって基本的にちょろいスクリプトには Python を使っていて, 数値計算系のコンテンツ作成も Python で慣れていたので, はじめは Python を使っていた. ちなみにその勉強の成果として Python を使った数学・応用数学入門については次のようなミニ通信講座を作っている.

C からの Python 移植で本を読んでいこう, となり, Python で数値計算となると numpy なので numpy でがんばって書き直していた. numpy となると for ではなくスライスを使った記法で書くのがふつうなわけで, それで書き直していたが, for (i = 1; i <= i_max + 1; i++) { }, for (i = i_max / 2 - 10; i <= i_max / 2 + 10; i++) { } みたいなのを適切なスライスで書き直すのがあまりに面倒で, 嫌になった. これはもう速い言語で for を直接移植した方が楽で早いことに気付き, 前から勉強しようと思っていた Rust に移行した. これならほとんど何も考えずに移植できて, 余計なバグ取りに勤しまずに済んだ.

基本的な計算は問題なくなったので問題は可視化の部分になった. 元の本は Gnuplot だったが, 計算と出力を分離していない産業廃棄物である. (工学者が産廃を作ることに抵抗を感じなかったのか不思議でならない.) Rust では CSV を書き出すことにして, 書き慣れていることから可視化は Python・matplotlib で処理して mp4 化することにした. ただし, これがとにかく遅い. 本では「2 秒分」計算しているのだが, これがとにかく遅い. 「0.1 秒分」の計算くらいならもちろん速くなるが, どんどんなまっていったり, 数値拡散を見るためにある程度長時間挙動を見たいので, これでは困る. いい加減嫌になったので, 自分にできる範囲の高速化としてできることを考えた.

はじめは Rust では CSV を作り, 別途 Gnuplot スクリプトで png 化して, これを ffmpeg で mp4 にしようかと思った. png さえ作れれば動画化は ffmpeg でかなり楽になる. 以前使ったときもここ自体はかなり速かった肌感覚があるので, png 化を楽に速くすることを考えた.

しかし, Gnuplot のスクリプト自体の情報が調べにくく, Gnuplot スクリプトとして動的に大量のファイルを読み込んで処理, みたいなこともしにくいようで, シェルスクリプトをかませないと駄目そうなど, プログラミング弱者の市民には本当に鬱陶しい. はじめは mp4 化と同じく直接のコマンド実行を検討したが, 少し調べたら gnuplot ラッパーがあったので, それを使うことにした.

ここまで来たらあとは ffmpeg をやればいい. 別途実行するのも面倒なので, Rust スクリプト内で直接実行することにした.

もう少し経緯, そして今後の予定

実ははじめは Paraview を使おうと思っていた. OpenFOAM で使っていて, 多少なりとも操作への感覚があったからだ. ただ Rust からの出力として vtk なり vtu を作らなければならず, それで二の足を踏んでいた. しかし Paraview は csv 読み込みもできると聞き, それなら, と思って試してみたが, 1 次元だと使い物にならなかった. いま勉強している純粋な 1 次元の数値計算の可視化に 3 次元用のビューアーを使おうとすること自体に無理があるので, その線は諦めた.

3 次元の計算をするときは多少は慣れているし Paraview も考えたい. ただ, csv からだと読み込みが少し面倒そうなので, やはり vtu を作るのがいいように思う. vtu の仕様の読み込みからしないといけなくて面倒で嫌なのだが ライブラリがあると嬉しい.

それはそれとして, 上掲書は後半で 2 次元を議論する. 2 次元もいまの 2d plot を使うか, 見たいモノによって 3d plot にすればいいので, しばらくは Paraview はいらないだろう.

今の本を読み終わったら粒子法または格子ボルツマン法を勉強する予定で, 既に本も買っている. ここの可視化をどうするといいのかが気になっている. 粒子法については SPHysics や PySPH があるので, これの可視化を参考にしようと思っている. まだ全く見ていないが, 何で可視化しているのだろうか?

格子ボルツマン法はまだ何もわかっていないのだが, 何で粒子法をやりたいかというと, 有限要素法などでメッシュを切りたくないからだ. ちょっと複雑な形状になるとメッシュ生成があまりにも面倒で, メッシュデータの読み込み処理なども書きたくない. 粒子法は粒子法でいろいろ鬱陶しい点があるはずだが, CG 界隈への展開もあるようなので, 別の発展方向が転がっているような気がしていて, それも 1 つのポイントだ.

いま, 特に理工系を志す中高生向けの数学・物理・プログラミングコンテンツを作っているし, たまたま地元の広義知人が区議になったこともあり, その人に地元の理工系教育推進のための相談・提案的なこともしている. 子どもの頃にほしかった・やってみたかったことでもあり, (現状専門外もはなはだしいが) まだまだこの辺の技術・研究ニーズもあるはずで, 中高生へのキャリア教育・技能のマネタイズ教育みたいな視点からも大事そうだと思っている. Pixar の展示会などもあり, ゲームやアニメ・映画など中高生にもある程度身近で需要も見やすいネタなので, いろいろ考えている. もちろん理工系大学生にも一定のニーズがある話題だろう.

私は出自的に, 数学と物理に関して, 能力はどもかく人並以上の尋常ならざる耐性はあるが, 数値計算はまた別種の厳しい世界がある. 何より工学者はゴミのようなコードしか書かないという話しか聞かないし, 数学系だと数値関数解析のような収束系の理論が中心のようでコードがあまり出てこないし, 物理系だと「具体的なコードを示すと時代遅れになりがちだし, 言語の選択については人の趣味もある」といってこれまた具体的なコードが出て来にくいし, いま流行りの統計学・機械学習まわりだとその筋の人達が興味を持つような話ばかりで, 私のような数学・物理系市民には不満なことこの上ないと, とにかく勉強がしづらい. 自ら血を流しつつ, きちんと将来的なマネタイズまで考えられるような, 生徒・学生向けコンテンツの整備を進めていきたい.

また, Python は Google Colab なり何なりで, これまたとにかく厳しい作業である環境構築なしでプログラミングに集中しやすい状況も出ていて, Python 利用の検討はこれはこれで外せないとは思っている. ただ, いまの私の力量ではとにかく遅く書くことしかできず, つらい. 数値計算, ちょっと何かやろうと思うとあっという間に重くなるので, そういう意味でもつらい. 地元的な事情からするとマシンを子ども達 (または家庭) 自身で準備するとなると, スペックも厳しいので, その意味でもいろいろ考えることはある.

理想郷は遥か遠い.

ぜひ登録またはツッコミを: GitHub・YouTube

最近の勉強の結果は次の GitHub リポジトリと YouTube リストにまとめている.

励みになるのでぜひチャンネル登録してほしい.

あと Rust 自体まだはじめて 1 ヶ月すら経っておらず, もとがゴミ以外表現しようのない工学者による C コードを基礎にした Rust コードなので, いまのコードもゴミだろう. 必要ならコード自体にもツッコミを入れてほしい. 泣きながら勉強して直していきたい.

緩募: Rust・Gnuplot・ffmpeg による 1 次元波動方程式の差分解法のバグ潰し

追記

次のツイートの指摘を受けて修正した.

これも結論を言うと, そもそもの離散化の仕方・計算の仕方がまずかった. 詳しくは上記黒木さんのツイートまたは GitHub 参照. Julia も気になっているが, まずは Rust でベクトル計算できないかが気になる. てらモスの兄貴に聞くのがいいだろうか?

結論から

結論から言うと, 次のコードが意図通りに動かない.

これが次のような挙動をする.

意図しているのは次のような挙動.

バグなしコードの説明

追記: 最初にも追記したように, 修正済. 詳しくは上記黒木さんのツイートまたは GitHub 参照.

まずはバグなしコードから.

このページを参考にして, 素直に波動方程式を解いている. 初期条件は $u = 0, u_t = 0$ で, 原点で波を $(1/2) \sin (2 \pi f t)$ で駆動していて, 次のように離散化している.

\begin{align} u(t + \Delta t, x) = 2 u(t,x) - u(t - \Delta t, x) + \left( \frac{\Delta t}{\Delta x} \right)^2 (u(t,x + \Delta x) - 2 u(t,x) + u(t,x)). \end{align}

意図通りに動いているので問題ないだろうと思っている.

バグありコードの説明

初期条件や波の駆動は 1dim_wave_eq_only_u.rs と同じ設定にしている. 2 階の常微分方程式でよくやるように, $u, v = u_t$ として, 次の連立の偏微分方程式を解く形にしている.

\begin{align} u_t = v, \quad v_t = u_{xx}. \end{align}

これを次のように離散化した.

[u(t + \Delta t, x) = u(t, x) + v(t, x) \Delta t,]

[v(t + \Delta t, x) = v(t, x) + \frac{u(t, x - \Delta x) - 2 u(t,x) + u(t, x + \Delta x)}{2 \Delta x} \Delta t.]

実際には次のように実装している.

どこがまずいのかがよくわかっていない. どなたかわかる方, ぜひ教えてほしい.

Rust・Gnuplot・ffmpeg での数値計算: 波動方程式とその可視化

引き続き, Rust の勉強を兼ねて数値計算コードをゴリゴリ書き, 可視化している.

動画・コード・ドキュメントへのリンク

動画は次のリンクから見られる.

ドキュメントとコードは GitHub に置いてある.

議題: 可視化

ここで議論 (質問・相談) したいのは可視化だ. 結論からいうと, RustGnuplot の 3 次元プロットを使いたい. Axes3D というのはあるのだが, いろいろ試してみて, これで splot tmp.dat u 1:2:3 with lines を実現できなかった. 私の実装能力の問題ならいいのだが, とにかくどうしたらいいかわからない. 直接 Gnuplot のコマンドを叩きたいのだが, それは RustGnuplot では対応していないように見える: これも私がドキュメントを拾えていないだけならいいのだが, やはりどうしたらいいかはわからない. Command::new() から直接 Gnuplot ワンライナーで実行しようと思ったものの, これはこれでうまくいかない. どうすればいいか, というのが問題だ.

少なくとも matplotlib で mp4 化するのは 1 次元のデータ量でさえ遅かったので, 絶対にやりたくない. それもあって Gnuplot を導入した経緯がある. png 化するくらいならどうか, と思って試してみたが, これも matplotlib 力が低いせいか, 意図通りのプロットが出ない.

そこで諦めて, Python で csv から png を作る gp ファイルを作り, それを Python スクリプト内のコマンド実行で png 化し, それをさらに Python スクリプト内の ffmpeg コマンド実行で mp4 化した.

ffmpeg を Rust 内でコマンド実行するのはできているから, RustGnuplot のように Rust 内で Gnuplot を実行したい. ソースコードを読み, コード量は少ないが, いまの私が理解できるレベルの Rust コードではなく, 機能追加などは夢のまた夢だった. とても厳しい気持ちになっている.

できる限り Rust 内で関係させたいし, Gnuplot のラッパーがほしい. 夢は果てない.

追記: Rust 内での Gnuplot 直接実行

Gnuplot のワンライナーを実行する形で実装できた. コードの詳細については次のリンク先のファイルを見てほしい.

具体的には次のようなコードで Rust 内で Gnuplot が直接実行できるようになった.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
use std::process::Command;
fn write_png(cnf: &Config, cdata: &CalcData) -> Result<(), Box<dyn error::Error>> {
    let csv_name: String = format!("{}/{:08}.csv", &cnf.dir_name, &cdata.output_num);
    let png_name: String = format!("{}/img.{:08}.png", &cnf.dir_name, &cdata.output_num);
    Command::new("gnuplot")
        .arg("-e")
        .arg(r#"set terminal png;"#)
        .arg("-e")
        .arg(r#"set datafile separator ",""#)
        .arg("-e")
        .arg(r#"set ticslevel 0;"#)
        .arg("-e")
        .arg(r#"set dgrid3d 100,100;"#)
        .arg("-e")
        .arg(format!(
            r#"set zrange [{}:{}]"#,
            &cnf.graph_ulim_min, &cnf.graph_ulim_max
        ))
        .arg("-e")
        .arg(format!(r#"set output "{}""#, &png_name))
        .arg("-e")
        .arg(format!(r#"splot "{}" u 1:2:3 with lines;"#, &csv_name))
        .output()
        .expect("failed to start `gnuplot`");
    Ok(())
}

冗長だが, Gnuplot のスクリプトをワンライナーで実行する形になっている. Gnuplot は詳しくないので Gnuplot を書ける前提なのが厳しいが仕方ない.