Rustで数値シミュレーションをする際に使用しているライブラリ(クレート)の紹介

この記事はTUT Advent Calendar 2017の18日目の記事です。

概要

普段私はRustを使ってシミュレータを作成しているのですが、そこでお世話になっているライブラリ(クレート)を簡単に紹介します。
今回紹介するライブラリ(用途)は以下の2つです。

  • rayon (データ並列化)
  • statrs (様々な確率分布に従う乱数生成)

rayon

すでにいくつかの記事で紹介されていますが、rayonによって主に以下の処理を簡潔に行うことができます。両者とも、Work stealingと呼ばれる方法で実現されているそうです。

  • スレッドによる並列処理
  • 再帰関数の分割処理

スレッドによる並列処理

rayonでは、Parallel Iteratorの導入によって通常のイテレータに対する関数(mapなど)の処理を並列化することができます。D言語でいうstd.parallelismに似ているかもしれません。 ただし、Rustのシステムによって並列時の安全性は保証されています。つまり、標準ライブラリのスレッド同様、Sendトレイト、Syncトレイトの制約によってコンパイル時に検査が行われます。

iter関数やinto_iter関数が実装されている型は基本的にParallel Iteratorに変換することができ、以下のような対応関係を持ちます。

コード例は以下のとおりです。

extern crate rayon;
use rayon::prelude::*;

fn main() {
    (0..100)
    .into_par_iter()
    .map(|i| i * 2)
    .for_each(|x| println!("{}",x));
// 200未満の偶数が出力される
}

大変便利ですが、Parallel Iteratorに対してcollect::<Vec<_>>()のようなことはできず、 あらかじめ格納先のベクタをミュータブルとして宣言したのち、collect_into関数を用いて格納する必要があります。

コード例は以下のとおりです。

extern crate rayon;
use rayon::prelude::*;

fn main() {
    let mut even_vec = vec![0;100];
    (0..100)
    .into_par_iter()
    .map(|i| i * 2)
    .collect_into(&mut even_vec);
// 200未満の偶数がeven_vecに格納される
}

並列処理時のスレッド数は自動的に決定されますが、スレッド数を能動的に設定したい場合は、num_threads関数によって設定可能です。
また、実行時にスレッド数を設定したい場合は、環境変数RAYON_NUM_THREADSの値を設定してから実行することで対応できます。

再帰関数の分割処理

join関数によって、再帰関数を効率的に並列処理します。
rayonのGithubリポジトリの例が簡潔であるため、引用します。 コード例は、並列化したクイックソートのようです。

fn quick_sort<T:PartialOrd+Send>(v: &mut [T]) {
    if v.len() <= 1 {
        return;
    }

    let mid = partition(v);
    let (lo, hi) = v.split_at_mut(mid);
    rayon::join(|| quick_sort(lo), || quick_sort(hi));
}

リポジトリのREADMEによれば、

Note though that calling join is very different from just spawning two threads in terms of performance. 
This is because join does not guarantee that the two closures will run in parallel. 
If all of your CPUs are already busy with other work, Rayon will instead opt to run them sequentially. 
The call to join is designed to have very low overhead in that case, so that you can safely call it even with very small workloads (as in the example above).

とのことですので、join関数の引数のクロージャを、可能なら並列処理し、CPUの負荷が高い場合は逐次的に処理するよう、うまくスケジューリングされる(?)ようです。

statrs

statrsはRustにおける統計的な数値計算のためのクレートで、特に私はdistributionモジュールを頻繁に利用しています。distributionモジュールを利用することで、各種確率分布に従う乱数を容易に生成することができます。
ただし、rust statrsで検索すると検索エンジンが変に気を利かせて別の検索ワードを推奨してくるので注意です。

離散型分布、連続型分布問わず様々な分布に対応しており、以下の分布ごとに構造体が定義されています(非常に多い)。また、分布の期待値、分散を返す関数が各分布に実装されており、離散型確率分布に対しては確率質量関数の値、連続型確率分布に対しては確率密度関数の値を返す関数が実装されています。

さらに、sample関数(rand::Rngトレイトを実装する乱数生成器を引数)の実装によって、各分布に従う乱数を生成可能です。 rand::Rngトレイトを実装している乱数生成器であればよいため、StdRngやThreadRng、XorShiftなど様々なものに対応しています。

  • 離散型確率分布
    • ベルヌーイ分布
    • 二項分布
    • カテゴリカル分布(Categorical Distribution)
    • 離散一様分布
    • 幾何分布
    • 超幾何分布
    • 多項分布
    • ポワソン分布
  • 連続型確率分布
    • ベータ分布
    • コーシー分布
    • カイ分布
    • カイ2乗分布
    • ディリクレ分布
    • アーラン分布
    • 指数分布
    • F(Fisher-Snedecor)分布
    • ガンマ分布
    • 逆ガンマ分布
    • 対数正規分布
    • 正規分布
    • パレート分布
    • t分布
    • 三角分布
    • 連続一様分布
    • ワイブル分布

本記事ではdistributionモジュールにのみ着目しましたが、他のモジュールにてベータ関数やエラー関数などの複雑な関数の値を計算する関数も用意されているようです。

まとめ

rayonやstatrsは、数値シミュレーションをするうえで、大変有用です。並列処理により計算時間を削減することができ、randクレートでは対応していない様々な確率分布に従う乱数生成を行うことで様々なシミュレーションが可能です。
また、Rustはコンパイラが大変優秀なので、コンパイル時の検査のおかげでシミュレーション実行時のエラーがなく、重宝しています。その一方で、コンパイルを通すのは大変ですが。

私の身の回りでもRustが流行ればいいな、と思うこの頃です。