歩いたら休め

If the implementation is easy to explain, it may be a good idea.

【Rust】Rustの行列計算ライブラリについて調べる

この間までLotka-Volterraの方程式をRustで解いて勉強していましたが、そろそろPythonnumpy.ndarray のような行列計算を楽にする機構が欲しくなります。

もちろん、実際のソースコードで使われているライブラリ(crate)が良いです(メンテナンスされておらず、最新のコンパイラで動作しない場合もあるようなので…)。ということでまずはawesome-rustにある機械学習ライブラリの依存ライブラリをつてに探してみます。

github.com

AtheMathmo/rusty-machine

Rustの機械学習ライブラリだそうです。依存ライブラリはこんな感じ。

[dependencies]
num = { version = "0.1.41", default-features = false }
rand = "0.4.1"
rulinalg = { git = "https://github.com/AtheMathmo/rulinalg", rev = "1ed8b937" }

いきなり見つかりました。ちゃんと調べていないですが、逆行列(inverse)の計算も行えるようです。

github.com

A linear algebra library written in Rust https://crates.io/crates/rulinalg

avinashshenoy97/RusticSOM

こちらでは別のライブラリが使われています。

[dependencies]
rand = "0.4"
ndarray = "0.11"

github.com

こちらは追加のライブラリでシリアライゼーション( serde-1 )や並列化( rayon )、実験的のようですがblasの利用などできるみたいです。

autumnai/leaf

[dependencies]
collenchyma = { version = "0.0.8", default-features = false, features = ["native"] } # native feature to read/write data into tensors
collenchyma-blas = { version = "0.2.0", default-features = false, features = ["native"] } # only compiles with native feature
collenchyma-nn = { version = "0.3.2", default-features = false }

log = "0.3.2"
rand = "0.3.0"
num = "0.1"

capnp = "0.6.2"

timeit = "0.1.2"

clippy = { version = "0.0.41", optional = true }

[build-dependencies]
capnpc = "0.6.1"

こちらもまた別のライブラリに依存しています。

github.com

Provides a simple and unified API to run fast and highly parallel computations on different devices such as CPUs and GPUs, accross different computation languages such as OpenCL and CUDA and allows you to swap your backend on run-time.

OpenCLやCUDA等の計算言語をまたがって、統一的なAPIでCPUやGPU上で計算できるcrateのようです。ただ、以下の2017年の記事の時点で開発が止まっているようで、Githubのコミットログを見ても最近は開発されていないようです。

qiita.com

Leafフレームワークのコアの部分までRustで書かれていたのですが、当時はRustでCPU/GPUで高速に行列演算できるライブラリがなく、ディープラーニングのコアの部分(計算グラフの構築, 微分計算の記述, backward/update)自体よりもバックエンドの行列演算のライブラリ利用・開発に難儀していたようです。

そのような状況でTensorFlowのRustバインディングが公開され、Autumnの開発者がLeafの開発を停止したのですが、TensorFlowのRustバインディングは学習済みのモデルをloadしてinferenceに使用できるだけというもので、実際にはPythonなど他の言語でモデルの開発・学習・保存する必要があり、Leafの開発の中止は非常に残念でした。

tensorflow/rust

有名なニューラルネットワークのライブラリ tensorflow のRustバインディング

[dependencies]
libc = "0.2.43"
aligned_alloc = "0.1.3"
num-complex = { version = "0.2.1", default-features = false }
tensorflow-sys = { version = "0.15.0", path = "tensorflow-sys" }
byteorder = "1.2.7"
crc = "1.8.1"

どうやらlibc経由でC言語C++?)を呼び出していて、自前で行列計算は実装していないように見えます。

maciejkula/rustlearn

[dependencies]
rand = "0.3"
crossbeam = "0.2.9"
serde = "1.0"
serde_derive = "1.0"

どうやらこのライブラリ自体でdense / sparse matricesの実装がされているようです。

参考

どうやら、いまのところPythonにおける numpy のような、デファクトスタンダードにまでなっているライブラリは無いようですね。この辺りから、自分にとって使いやすいものがどれか試してみたいと思います。

また、Rust自体の進化が激しく、ライブラリの情報がすぐに古くなっていて最新情報が分かりづらいのは、言語初心者としてはけっこうつらいですね…。これは新興のプログラミング言語なのである程度は仕方ないかもしれませんが。

【Rust】Rustで数値計算プロジェクトを試す その10

中間目標であったLotka-Volterraの方程式がシミュレーションで解けました。

shimaphoto03.com

use gnuplot::{Figure, Caption, Color};


const A: f64 = 0.01;
const B: f64 = 0.05;
const C: f64 = 0.0001;
const X0: f64 = 1000.0;
const Y0: f64 = 100.0;
const DELTA_T: f64 = 0.001;


fn main() {
  let (mut x, mut y) = (X0, Y0);
  let mut xs = vec![];
  let mut ys = vec![];
  let mut ts = vec![];

  for i in 1..500000 {
    let t = (i as f64) * DELTA_T;
    let (dx, dy) = runge_kutta(x, y, |x, y| {
      (
        A * x - C * x * y,
        - B * y + C * x * y
      )
    });
    x += dx;
    y += dy;
    ts.push(t);
    ys.push(y);
    xs.push(x);
  }

  plot(xs, ys);
}

fn runge_kutta<F>(x: f64, v: f64, f: F) -> (f64, f64) 
    where F : Fn(f64, f64) -> (f64, f64) {
  let df = |x, v| {
    let (rx, rv) = f(x, v);
    (rx * DELTA_T, rv * DELTA_T)
  };
  let (k1x, k1v) = df(x, v);
  let (k2x, k2v) = df(x * k1x / 2.0, v + k1v / 2.0); 
  let (k3x, k3v) = df(x + k2x / 2.0, v + k2v / 2.0);
  let (k4x, k4v) = df(x + k3x, v + k3v);
  (
    (k1x + 2.0 * k2x + 2.0 * k3x + k4x) / 6.0,
    (k1v + 2.0 * k2v + 2.0 * k3v + k4v) / 6.0
  )
}

fn plot(x: Vec<f64>, y: Vec<f64>) {
  let mut fg = Figure::new();
  fg.set_terminal("pngcairo", "output.png");
  fg.axes2d()
    .lines(&x, &y, &[Caption("A line"), Color("black")]);
  fg.show();
}

x ウサギ(被食者)の個体数、y キツネ(捕食者)の数が、相平面上でどんな軌道を描いているのか分かります。

f:id:takeshi0406:20190217172539p:plain

【Rust】Rustで数値計算プロジェクトを試す その9

以前のコードを改造して、 xv それぞれで微分方程式を指定できるようにしました。クロージャ便利ですね。

use gnuplot::{Figure, Caption, Color};


const K: f64 = 1.0;
const M: f64 = 1.0;
const A: f64 = 2.0;
const X0: f64 = 0.0;
const V0: f64 = 1.0;
const DELTA_T: f64 = 0.01;


fn main() {
  let (mut x, mut v) = (X0, V0);
  let mut xs = vec![];
  let mut ts = vec![];

  for i in 1..1000 {
    let t = (i as f64) * DELTA_T;
    let (dx, dv) = runge_kutta(x, v, |x, v| {
      (
        v,
        (- K * x - A * v) / M
      )
    });
    x += dx;
    v += dv;
    ts.push(t);
    xs.push(x);
  }

  plot(ts, xs);
}

fn runge_kutta<F>(x: f64, v: f64, f: F) -> (f64, f64) 
    where F : Fn(f64, f64) -> (f64, f64) {
  let df = |x, v| {
    let (rx, rv) = f(x, v);
    (rx * DELTA_T, rv * DELTA_T)
  };
  let (k1x, k1v) = df(x, v);
  let (k2x, k2v) = df(x * k1x / 2.0, v + k1v / 2.0); 
  let (k3x, k3v) = df(x + k2x / 2.0, v + k2v / 2.0);
  let (k4x, k4v) = df(x + k3x, v + k3v);
  (
    (k1x + 2.0 * k2x + 2.0 * k3x + k4x) / 6.0,
    (k1v + 2.0 * k2v + 2.0 * k3v + k4v) / 6.0
  )
}

fn plot(x: Vec<f64>, y: Vec<f64>) {
  let mut fg = Figure::new();
  // fg.set_terminal("pngcairo", "output.png");
  fg.axes2d()
    .lines(&x, &y, &[Caption("A line"), Color("black")]);
  fg.show();
}

【Rust】Rustで数値計算プロジェクトを試す その8

昨日のやつから、関数の中で xy+= するように変更すれば、多少コードがスッキリするんじゃないかと思ってやってみました。

ミュータブルな参照で渡して、参照を外して…って感じなので書きづらい。多分間違った道に来たと思うので引き返します。

stackoverflow.com

use gnuplot::{Figure, Caption, Color};


const K: f64 = 1.0;
const M: f64 = 1.0;
const A: f64 = 2.0;
const X0: f64 = 0.0;
const V0: f64 = 1.0;
const DELTA_T: f64 = 0.01;


fn main() {
  let (mut x, mut v) = (X0, V0);
  let mut xs = vec![];
  let mut ts = vec![];

  for i in 1..1000 {
    let t = (i as f64) * DELTA_T;
    runge_kutta(&mut x, &mut v, |x, v| {
      - K * x - A * v
    });
    ts.push(t);
    xs.push(x);
  }

  plot(ts, xs);
}

fn runge_kutta<F>(x_ref: &mut f64, v_ref: &mut f64, f: F) where F : Fn(f64, f64) -> f64 {
  let x = *x_ref;
  let v = *v_ref;
  let k1v = f(x, v) * DELTA_T / M;
  let k1x = v * DELTA_T;
  let k2x = (v + k1v / 2.0) * DELTA_T; 
  let k2v = f(x + k1x / 2.0, v + k1v / 2.0) * DELTA_T / M;
  let k3v = f(x + k2x / 2.0, v + k2v / 2.0) * DELTA_T / M;
  let k3x = (v + k2v / 2.0) * DELTA_T;
  let k4v = f(x + k3x, v + k3v) * DELTA_T / M;
  let k4x = (v + k3v) * DELTA_T;
  *v_ref += (k1v + 2.0 * k2v + 2.0 * k3v + k4v) / 6.0;
  *x_ref += (k1x + 2.0 * k2x + 2.0 * k3x + k4x) / 6.0;
}

fn plot(x: Vec<f64>, y: Vec<f64>) {
  let mut fg = Figure::new();
  // fg.set_terminal("pngcairo", "output.png");
  fg.axes2d()
    .lines(&x, &y, &[Caption("A line"), Color("black")]);
  fg.show();
}

【Rust】Rustで数値計算プロジェクトを試す その7

クロージャを渡す実装をしました。 t の引数は結局使ってないので削除しました。

doc.rust-jp.rs

の記法はトレイト境界というもので、これが無いと以下のようなエラーが出るのですが、

error[E0277]: the size for values of type `(dyn std::ops::Fn(f64, f64) -> f64 + 'static)` cannot be known at compilation time
  --> src/main.rs:31:32
   |
31 | fn runge_kutta(x: f64, v: f64, f: Fn(f64, f64) -> f64) -> (f64, f64) {
   |                                ^ doesn't have a size known at compile-time

…なのですが、今のところ意味が分からずに書いてるので後でチェックします。

use gnuplot::{Figure, Caption, Color};


const K: f64 = 1.0;
const M: f64 = 1.0;
const A: f64 = 2.0;
const X0: f64 = 0.0;
const V0: f64 = 1.0;
const DELTA_T: f64 = 0.01;


fn main() {
  let (mut x, mut v) = (X0, V0);
  let mut xs = vec![];
  let mut ts = vec![];

  for i in 1..1000 {
    let t = (i as f64) * DELTA_T;
    let (dx, dv) = runge_kutta(x, v, |x, v| {
      - K * x - A * v
    });
    x += dx;
    v += dv;
    ts.push(t);
    xs.push(x);
  }

  plot(ts, xs);
}

fn runge_kutta<F>(x: f64, v: f64, f: F) -> (f64, f64) 
    where F : Fn(f64, f64) -> f64 {
  let k1v = f(x, v) * DELTA_T / M;
  let k1x = v * DELTA_T;
  let k2x = (v + k1v / 2.0) * DELTA_T; 
  let k2v = f(x + k1x / 2.0, v + k1v / 2.0) * DELTA_T / M;
  let k3v = f(x + k2x / 2.0, v + k2v / 2.0) * DELTA_T / M;
  let k3x = (v + k2v / 2.0) * DELTA_T;
  let k4v = f(x + k3x, v + k3v) * DELTA_T / M;
  let k4x = (v + k3v) * DELTA_T;
  let v = (k1v + 2.0 * k2v + 2.0 * k3v + k4v) / 6.0;
  let x = (k1x + 2.0 * k2x + 2.0 * k3x + k4x) / 6.0;
  (x, v)
}

fn plot(x: Vec<f64>, y: Vec<f64>) {
  let mut fg = Figure::new();
  fg.set_terminal("pngcairo", "output.png");
  fg.axes2d()
    .lines(&x, &y, &[Caption("A line"), Color("black")]);
  fg.show();
}

【Rust】Rustで数値計算プロジェクトを試す その6

ひとまず、Rust習得の第二の中間目標として「Lotka-Volterra方程式を4次のRunge-Kutta法で近似計算する」ことを置いてみます。いきなり2重振り子を計算するのは、極座標とかラグランジアンみたいな概念が必要っぽくて、「大学時代に触れた覚えがあるけど急に思い出せない」ので後にします。

ひとまずRunge-Kutta法の実装です。

shimaphoto03.com

こちらの調和振動子の実装をパクりました。

qiita.com

所有権周りで苦労すると思ってましたが、Copy Traitが実装されているおかげで、関数に渡した後でも使えるので詰まることはありませんでした。

doc.rust-jp.rs

汎用性をもたせるために、関数ポインタとかいろいろ試してみたいですね。

use gnuplot::{Figure, Caption, Color};


const K: f64 = 1.0;
const M: f64 = 1.0;
const X0: f64 = 0.0;
const V0: f64 = 1.0;
const DELTA_T: f64 = 0.01;


fn main() {
  let (mut x, mut v) = (X0, V0);
  let mut xs = vec![];
  let mut ts = vec![];

  for i in 1..1000 {
    let t = (i as f64) * DELTA_T;
    let (dx, dv) = runge_kutta(x, v, t);
    x += dx;
    v += dv;
    ts.push(t);
    xs.push(x);
  }

  plot(ts, xs);
}

fn runge_kutta(x: f64, v: f64, t: f64) -> (f64, f64) {
  let k1v = f(x, v, t) * DELTA_T / M;
  let k1x = v * DELTA_T;
  let k2x = (v + k1v / 2.0) * DELTA_T; 
  let k2v = f(x + k1x / 2.0, v + k1v / 2.0, t + DELTA_T / 2.0) * DELTA_T / M;
  let k3v = f(x + k2x / 2.0, v + k2v / 2.0, t + DELTA_T / 2.0) * DELTA_T / M;
  let k3x = (v + k2v / 2.0) * DELTA_T;
  let k4v = f(x + k3x, v + k3v, t + DELTA_T) * DELTA_T / M;
  let k4x = (v + k3v) * DELTA_T;
  let v = (k1v + 2.0 * k2v + 2.0 * k3v + k4v) / 6.0;
  let x = (k1x + 2.0 * k2x + 2.0 * k3x + k4x) / 6.0;
  (x, v)
}

fn f(x: f64, _v: f64, _t: f64) -> f64 {
  - K * x
}

fn plot(x: Vec<f64>, y: Vec<f64>) {
  let mut fg = Figure::new();
  fg.set_terminal("pngcairo", "output.png");
  fg.axes2d()
      .lines(&x, &y, &[Caption("A line"), Color("black")]);
  fg.show();
}

f:id:takeshi0406:20190211171743p:plain

【Rust】Rustで数値計算プロジェクトを試す その5

昨日の件を友達に話したところ、crate側がRustのバージョンに対応してなかった可能性が高いようです。

これってちょっと見てみたけど、crate側の実装の問題な気がする。gnuplotのラッパーの方使ったら出力できた。

https://www.mathgram.xyz/entry/rust/gd

去年の中盤から更新されてなかったから

https://www.utam0k.jp/blog/2018/05/28/rust_std_default/

この一番下のエラーな気がするんだよな

というわけで、次のgnuplotライブラリを使って、サインカーブを出力してみます。

docs.rs

use gnuplot::{Figure, Caption, Color};

fn main() {
  let mut x = vec![];
  let mut y = vec![];

  for i in 0..1000 {
    let v = (i as f64) / 100.0;
    x.push(v);
    y.push(v.sin());
  }
  plot(x, y);
}

fn plot(x: Vec<f64>, y: Vec<f64>) {
  let mut fg = Figure::new();
  fg.set_terminal("pngcairo", "test.png");
  fg.axes2d()
      .lines(&x, &y, &[Caption("A line"), Color("black")]);
  fg.show();
}

できました!これで「サインカーブをプロットする」のというトイケースはクリアですね。

f:id:takeshi0406:20190209102830p:plain