Rustやっつけフェーズ(11)カルマンフィルタ⑤粒子フィルタ

次の章に進むと前章で理解したことをすっかり忘れてしまうカルマンフィルタの無念を成仏させていくシリーズ。
前回は以下のように無香料カルマンフィルタが90%成仏した。 引き続き、以下ドキュメントの粒子フィルタに関するPythonプログラムをRustで写経し成仏させていく。 今回はカルマンフィルタではないので残差からの次回状態予測はない。

粒子で作った網に引っかかった観測値の付近にある粒子を複製して網を再作成しその平均値を推測値とするように見える。
ほんでまたその網を状態変数の方向に投げるという繰返しとなっているよう見える。

著者いわく「本章は広大なトピックの表面をなでただけに過ぎない」とおっしゃっていて、粒子フィルタに関連するサイトをいくつか見ていくと
こんなのとか
「2 重モデルに基づくパーティクルフィルタを用いた不規則運動物体の追跡(PDF)」
こんなのとか
「パーティクルフィルタ 理論と特性(PDF)」

上記のお方たちの他にもいくつもいくつもあり、昔ハッブルか何かの銀河があまねくひしめいた画像を見て気が遠くなり死にたい気分になってきたことを思い出した。

今回は例題が少ないので以下の2つを追加した。
①前回のUKF回で最後に使ったロボットのシミュレータ
②以下のお方がやってるような粒子フィルタを扱ったサイトでよく出て来る例題

①ではwio-terminalとsdl2で行い、精度の違いによる差が出ることが解った。解ったからといって何かしたわけではない。そういうこともあるってことを言いたかったというか自分では何も出来なかったのだ。

また、wio-terminalの192KのRAMでは粒子3200個が限界。 sdl2ではスタックオーバーフローの壁にぶつかり70000個が限界。mainをまるっとスレッドにすればいいとの話しもあり1回試してみたが上手く行かなかったのであきらめたのだ。
(スタックオーバーフローにはならないがsdl2が表示されない)

②では参考のサイトに比べて粒子の幅が大きい気がするが、推定値がまあまあに見えるのでオッケー。

やっつけフェーズとは

来年7月は、どうせなら東京に落ちた方がいろいろあきらめもつくんじゃないかという思いを抱かせる方向に誘導されながらも、動くものが出来てそれっぽく動作すればオッケーなフェーズ。

粒子フィルタ

出来上がり

コード

以下はシミュレータ

  • sim_opening
    12.3 一般的な粒子フィルタアルゴリズム
  • sim_example_1
    12.6 SIR フィルタ 完全な例
  • sim_robot
    UKF回でのロボットシミュレータを使った例
  • sim_ffdtm
    粒子フィルタを扱ったサイトでよく出て来る例題

以下は描画プログラム
Cargo.tomlでシミュレータを切り替えて使用。

  • draw-particlefilter
    wio-terminal用粒子フィルタ描画
  • draw-particlefilter-animation
    粒子フィルタをアニメーション表示
  • draw-particlefilter-line
    粒子フィルタを線グラフ表示

    以下はsdl2用

  • draw-particlefilter-sdl2
  • draw-particlefilter-sdl2-animation
  • draw-particlefilter-sdl2-line

今回もこのようにして粒子フィルタが成仏した。
以上

Rustやっつけフェーズ(11)カルマンフィルタ④無香料カルマンフィルタ

なんとなくわかってきたと思ったら次の章で倍々にモヤモヤが増々していくカルマンフィルタを写経で成仏させていくシリーズ。
前回は以下のように多変量カルマンフィルタが成仏した。

引き続き、以下ドキュメントの無香料カルマンフィルタに関するPythonプログラムをRustで写経し成仏させていく。

著者いわく「EKFを実装したことがあるなら、UKF はなんて簡単なんだ」と書いてあったので簡単なんだと期待して始めたのだが...難しいじゃん!特に実装が。
てなわけでなんとか動くものは作れたが、何点か制約ができてしまった。

①シグマポイントのα設定を0.07より小さくするとコレスキーが死ぬ。
②ランドマークが見えない場合の対応が出来ない。
③シグマポイントはMerweScaledのみ。

①については設定値と回数に依存するっぽい。精度の問題なのか実装がどこかおかしいか解らんまま。
atan2が怪しいような、もしくは剰余計算時の雑な扱いのせいか...
以下のお方がおっしゃるようにSVDでやればいいじゃんと思って1度実装してみたのだが、型引数とデフォルトアロケータの扱いがめちゃ増えてしまって手に負えず断念した。

②についてはドキュメントの例題もランドマークが見えない場合のパターンが無いように見え、どうしたらいいか解らん。

③についてはトレイト使ってシグマポイントの入れ替えを可能にするだけといえばそうなのだが、型引数の対応をさらに行わなければならず断念した。

あとは10.8 複数の位置センサーの航空機位置をかえたときの結果がドキュメントと違って大幅にズレている。
これも解らんまま。

出来たコードは型引数とデフォルトアロケータの記述でおなかいっぱい感があると思うがそこは慣れなので...

やっつけフェーズとは

発生した問題を常に先送りし、周囲からの信頼や成長の機会を失いないながらも、動くモノが出来てそれっぽく動作すればオッケーなフェーズ。

無香料カルマンフィルタ

出来上がり

コード

ここから以下は無香料カルマンフィルタexampleプログラム

  • draw-robot-sim-ukf
    10.7 UKFを使う
  • draw-radaracc-sim-ukf
    10.8 航空機の追跡
  • radaraccsim
    航空機のシミュレータ
  • draw-radaracc-cv-sim-ukf
    10.8 機動する航空機の追跡
  • draw-sensor-without-doppler
    10.8 センサー統合:ドップラー未使用
  • draw-sensor-with-doppler
    10.8 センサー統合:ドップラー使用
  • draw-sensor-multiple
    10.8 複数の位置センサー
  • draw-self-positioning
    10.15 ロボットの自己位置推定─完全な例
  • robotukf/robotukffn
    上記10.15のfxやhxなどの関数定義を切り出し。
  • robotukf/robotukfsim
    上記10.15の運動モデル、制御入力、シミュレーション全般の設定を切り出し。

このようにして無香料カルマンフィルタが90%成仏した。

以上

Rustやっつけフェーズ(11)カルマンフィルタ③多変量カルマンフィルタ

カルマンフィルタわかったようでわからん。そんなモヤモヤを成仏させるシリーズとして前回は以下のように1次元カルマンフィルタが成仏した。

今回も前回から引き続き、以下ドキュメントの多変量カルマンフィルタに関するPythonプログラムをRustで写経し成仏させていく。

今回から多変量ということで行列演算が発生する。
なのでまたまたnalgebraの登場となるのだが、no_stdであるためダイナミックMatrixは使えそうな気がしない。
で、任意の行列次元を、関数を利用する側から指定するのってどうしたらええの?って探したところ以下の記事がものすごく参考になった。 3年前の記事。すごい。天才。
いつも思うのだが、このお方達はどのようにしてこのような技を習得するのだろうか...

やっつけフェーズとは

なんとかのひとつ覚えのみで乗り切ろうと思ってはいるのだが、それだけではどうにもならないこともあるもので、とりあえずそのような技を拝借しつつも動くものが出来ればオッケーなフェーズ。

多変量カルマンフィルタ

出来上がり

コード

  • mv_kalmanfilter
    多変量カルマンフィルタ関数ライブラリ
    • ジェネリクスで各行列の次元を利用側から指定可能。
    • セッター/ゲッターは廃止し、各行列をまるっとpublic。
      利用側から設定/取得。
  • discrete_white_noise
    離散白色ノイズ
    • ノイズのブロック分のみ作成。
    • 利用側の変数型定義により、2x2 or 3x3 or 4x4 のノイズブロックを作成。
    • block_diag的なことは利用側で行う。
  • robot2d
    xy平面を動くロボットのシミュレータ
  • constantacc
    定常加速シミュレータ

ここから以下の「draw-〜」は多変量カルマンフィルタオブジェクトを生成、設定し画面表示を行うexampleプログラム

  • draw-dog-sim-mv
    6.6 カルマンフィルタの実装
  • draw-dog-sim-1d2d
    6.9 隠れ変数の効果
  • draw-dog-sim-ellipses
    6.12共分散行列の詳しい調査
  • draw-robot-sim-ellipses
    8.1 ロボットの追跡
    • ここではブロックノイズをblock_diag的なことをして、プロセスノイズに設定する。
  • draw-acc-sim
    8.3 フィルタの次数の評価

ここから以下の「emb-〜」はグラフ表示ライブラリ

以下はテスト用プログラム

このようにして多変量カルマンフィルタが成仏した。
以上

Rustやっつけフェーズ(11)カルマンフィルタ②1次元カルマンフィルタ

カルマンフィルタ成仏シリーズとして開始した前回は以下のように離散ベイズフィルタが成仏した。

今回も前回から引き続き、以下ドキュメントの1次元カルマンフィルタに関するPythonプログラムをRustで写経し成仏させていく。

この後に控える多変量カルマンフィルタを前に、1次元のガウス分布である平均と分散を用いて、事前分布と観測値の残差を求め、その間での次の予測値を求めるパーセンテージを計算し、そのパーセンテージ(これをカルマンゲインと呼ぶ)に従い次の予測値(事前分布)を決定することを学ぶ。

前回の離散分布はドアが10個あれば、10個分の変数か必要だったが、ガウス分布を用いるとその釣鐘状の形態によって1つの平均と1つの分散で事足りて計算量が格段に経済的になる。
ただしその分制限があるというものらしい。

正しい説明はできるはずもないのでこんな解釈で進める。正しい内容を知りたければ上記ドキュメントを参照のこと。

また、事前分布予測→観測値尤度計算→事後分布への更新の繰り返しを行って行くと分散値はある値に収束して行き、分散値が収束するとカルマンゲインも収束する。
で、実験で収束したカルマンゲインを最初から固定値として用いると精度は良いままで計算量がこれまた節約できるというシロモノであるそうだ。

やっつけフェーズとは

なんとかのひとつ覚えのみで乗り切ろうをモットーに、車輪のダウングレード再発明をデフォとして、動くものが出来ればオッケーなフェーズ。

1次元カルマンフィルタ

出来上がり

コード

  • one_dimensional
    1次元カルマンフィルタ関数ライブラリ
    • ガウス分布構造体の定義
    • フィルタ関数:predict(予測)、update(更新)のトレイト定義
    • 上記トレイトを実装した3種類のオブジェクト。
      • ガウス分布の積(尤度x事前分布)パターン
      • カルマンゲイン算出パターン
      • 固定ゲイン設定パターン
    • シミュレーショントレイトの定義
  • dogsimulation
    犬の動きのシミュレーション。電圧のシミュレーションもこいつで事足りたので「dog」は余計だったが仕方ない。
  • nonlinersimulation
    非線形系のシミュレーション

  • od_kalmanfilter 1次元カルマンフィルタ操作オブジェクト

    • 生成時にフィルタ関数オブジェクトとシミュレーションオブジェクトを受け取り保持する。
    • セッターにてプロセスモデル、観測ノイズ分散、初期推定値、移動量などを設定する。
    • iterationsメソッドにて1ステップ分の以下動作を行う。
      • シミュレータを1ステップ動かし観測値を取得。
      • 事後分布、プロセスモデルから事前分布を予測。
      • 観測値、センサ誤差から尤度(もっともらしさ)を設定。
      • 尤度、事前分布から事後分布を計算・更新。
      • 指定により初期観測値を初期推定値に設定可能。

ここから以下の「draw-〜」は1次元カルマンフィルタオブジェクトを生成、設定し画面表示を行うexampleプログラム

  • draw-dogsim
    4.6 初めてのカルマンフィルタ
  • draw-dogsim-rest
    上記のカルマンゲイン算出バージョン
  • draw-dogsim-var
    カルマンフィルタのアニメーション
  • draw-voltsim
    4.11 カルマンフィルタの設計入門
  • draw-voltsim-var
    上記のアニメーション
  • draw-dogsim-bigsensornoize
    4.12 極端に大きい観測ノイズ
  • draw-dogsim-smallprocessnoize
    4.13 極端に小さいプロセス分散
  • draw-dogsim-wronginit
    4.14 大きく間違った初期推測値
  • draw-dogsim-bignoize_wronginit
    4.15 大きいノイズと大きく間違った初期推測値
  • draw-dogsim-bignoize_zinit
    最初の観測値を初期推測値にする
  • draw-dogsim-nonliner
    4.17 非線形な系

ここから以下の「emb-〜」はグラフ表示ライブラリ

以下はテスト用プログラム

このようにして1次元カルマンフィルタが成仏した。
以上

Rust やっつけフェーズ(11)カルマンフィルタ①離散ベイズフィルタ

前回はnalgebraを使用して球面束縛カメラを作った。 今回はほかにも色々やりたいことがあるのだけど、以下のドキュメントを見出してから何も手がつかなくなった。 ひと通り読みはしたのだが、その時はわかった気になったが次の日にはもうわからない。
Pythonで書かれたサンプルコードは今まで見てきた書き方と結構違いクセが強い。
科学系のプログラムはPythonで書かれたものが多いので、参考にするため読みはするが作ったことはないからちょっと難しい書き方されると途端に見失なう。

このいろいろな無念をRustで写経を行い成仏させていこうと思う。
まずは基本となる離散ベイズフィルタから。

やっつけフェーズとは

なんとかのひとつ覚えのみで乗り切ろうをモットーに、汎用性や可用性など亡き者にして、動くものが出来ればオッケーなフェーズ。

離散ベイズフィルタ

出来上がり

デザイン

変数名の整理

変数名 意味
prior 事前分布
belief 事前確率分布、信念
pdf 信念(belief)が確率分布であることを強調した名前
posterior 事後分布、事後確率分布
predict 予測
likelihood 尤度(もっともらしさ)
prob 確率分布
offset 移動量
lh_hallway 計測値が廊下の各位置にマッチする尤度を計算

処理

  • discrete_bayes(離散ベイズ関数)

    • predict(予測)
    • lh_hallway(尤度)
    • update(確率分布の更新)
    • normalize(正規化)
    • argmax(配列の最大値のインデックス)
  • 上記の関数を利用して、以下の更新を繰返していく。

    • 事後分布 = (尤度 * 事前分布) / 正規化係数
    • 事後分布を事前分布に入れる

フィルタ結果のテキスト表示と棒グラフ表示は欲しいので以下を用意。

  • 棒グラフ表示

    • 目盛レンジ、表示位置・サイズを指定可能。
      レンジとサイズから棒グラフの幅と位置を自動判定する。
    • 目盛表示はembedded-plotsを利用。なかなかいい感じの表示をしてくれるクレート。
  • テキスト表示

    • 表示位置・サイズを指定可能。
    • 1ページ終了したら、また最上部より色を変えて行を描画。

コード

discrete_bayes
robot
robot_filter
draw-robotfilter(robot_filterのグラフ描画)
dog_filter
draw-dogfilter.rs(dog_filterのグラフ描画)
emb_bargraph(棒グラフオブジェクト)
emb_textterm(テキスト表示オブジェクト)

おまけ
wio-draw-bargraph(棒グラフテストプログラム)
m wio-draw-textterm(テキスト表示テストプログラム)

今回はこんな感じで離散ベイズフィルタが成仏した。

Rust やっつけフェーズ(10) 球面束縛カメラ (wio-terminal, nalgebra, MVP)

前回はナビエストークスSOR法で流体のシミュレーションを行うPythonプログラムをRustで写経した。
配列計算はnalgebraを使用。 nalgebra繋がりで今回は球面束縛カメラで描画対象を3D表示することを行う。
本来はクォータニオンを習得しようとして始めたのだが、何故かこういう着地になってしまった。

やっつけフェーズとは

Rustと組込みを始めて約1年半だが、心はいつまでも初心者。
通勤中にスマホの単機能なエディタでえいやで作り、動けばオッケーなフェーズ。

球面束縛カメラ

出来上がり

デザイン

以前に作成した以下の描画を3D化する。

WioSBEye(カメラ位置情報)、WioSBButton(ボタン操作)
  • 左、右、上、下旋回
    5wayボタン操作により球面旋回し、位置情報を更新する。
  • ズーム、アウト
    TopLeft、TopMidボタン操作により、原点からの距離と位置情報を更新する。
WioSBCamera(3D変換)
  • モデル・ビュー・プロジェクション行列を生成し、受取った頂点情報を変換する。
  • 変換方法はnalgebra本家の以下ドキュメント「Build a MVP matrix」を基本的に踏襲。
  • モデル変換は、元画像が左上原点で作られているため、中央原点に対応。具体的にはX軸を中心に180度回転した後xとyをシフトする。(この回転時にクォータニオンが使われているようだ)
  • 遠近感を出すためwで割るのだが、割る前にwがマイナスにならないようかつ、wが小さくなりすぎないよう調整する。
  • 最後に左上原点表示に戻すためスクリーン行列を掛ける。

コード

wio_sbcamera
オブジェクトを斜めにして近づきすぎると破綻するバグあり。 何らかの限界があるのだろうがスキル不足で対処不可能。
wio_sbeye
wio_sbbutton
wio_elldiski(修正)
wio-polywave(修正)
main(流用)

今回はこんな感じで。

Rust やっつけフェーズ(9) ナビエストークスSOR法 (wio-terminal, nalgebra)

前回はPythonで書かれたチューリングパターンの描画をRustで写経した。

今回もこちらのお方がPythonで書かれたものを写経。

前回と同様、配列計算はnalgebraを使用。

ナビエストークスSOR法

出来上がり

コード

IterToolsのカーテシアンプロダクトを使いたかったが、stdを求められてあえなく断念。

wio_cfdnssor

main

今回はこんなもんで。