c++のスレッド

 メモ

  • c++11からはスレッドが言語に取り込まれている.
  • コンパイルには-pthreadを付ける.
  • スレッドを一時的に止めるときにはsleepではなく,スレッド用のstd::this_thread::thread_for()を使う.
  • thread::join()をちゃんと実行して対象のスレッドの終了を待つようにする.(またはdetachする.detachはスレッドを管理対象から外す.但し実行は裏で勝手に続く.これどんな時に使うんだろう?)
  • 無限ループに入れるような場合は,外から止められるように定期的に停止フラグを見に行かせる.
  • std::threadのコンストラクタに関数を渡すことでその関数を別スレッドとして実行する.関数の引数はコンストラクタの第2引数以降で渡す.
  • クラス内でメンバ関数をスレッド化するときは特殊なので注意.リファレンスで渡すこと,第2引数にthisが要ることを忘れちゃダメ
#include <iostream>
#include <thread>
#include <unistd.h>

bool stop_flg = false;

void worker(int i, int j) {
  while (1) {
    std::cout << "I'm " << std::this_thread::get_id() << " thread." << std::endl;
    std::cout << i+j << std::endl;
    if (stop_flg)
      break;
  }
  return;
}


int main(int argc, char* argv[])
{
  std::thread th(worker, 1, 2);
  sleep(2);

  stop_flg = true;
  th.join();

  return 0;
}

情報量関連のメモ

概要

機械学習とか統計を勉強していると情報量とかエントロピーの話が出てきていつも復習をし直しているので,纏めておく.

全体の1行まとめ

各量の1行まとめ

  • 事象eの確率がP(e)で与えられるとき,その事象の情報量は-log(P(e))。不確実性の指標になる。
  • (情報)エントロピー(H)は,情報源の平均情報量,つまり,情報量の期待値=エントロピー。不確実性の期待値。
  • 情報源が2つ以上の場合には,一方(x)について知った場合の,もう一方(Y)のエントロピー=条件付きエントロピー(H(Y|X))が定義される。
  • 相互情報量I(X,Y)は,H(X) - H(X|Y)で定義され,Xの情報量(の期待値)からYを知っている時のXの情報量(の期待値)の差。絶対に0以上。(Yを知ることでは知れないXの情報量ってこと?)
  • 相対エントロピー(=カルバック・ライブラー情報量)は,2つの情報源x, yの確率分布の距離。

  • 確率変数が2つ以上の場合,結合エントロピーが定義される,つまり,結合エントロピーエントロピーの多変数版。

情報量

  • 事象Aが与えられた時に,その情報量って?というものに答える.
  • 確率を使って考える.発生が稀な事象と,よくある事象があった時にどっちに遭遇した場合が情報の量が増えたと思うか?ここで,稀な事象の方がゴミみたいな事象なんだから情報は無い,と考えるんじゃなくて,めったに遭遇できない事象に出会えたらそこで知れる情報の量は大きい(価値が高いって感じ?)と考える.
  • そこで,発生確率が小さいほど大きくなるような単調減少関数で情報量を定義しよう,という発想になる.
  • 次に,どんな単調減少関数がいいかな?と考える.そこで,独立性な事象における情報の加法性を考えてみる.事象AとB,それぞれの情報量I(A),I(B)があった時に,I(A∧B)は?と考えると,独立なんだから,合わさって情報量が増えることはないでしょう,単純に和でしょ?と考える.独立な確率の同時確率はP(A)*P(B).これを情報量化すると和になってくれればいい,,,,ってことでよくある-log(P)が定義される.
  • 事象A,その発生確率をP(A)とすると,情報量I(A)は以下で定義される.定義されているということが重要.あくまでも定義なんだ
I(A) = -log(P(A))    (1)

エントロピー

  • ある情報源(事象発生源)Xが与えられた時に,その情報源のランダムさは?(逆に言うとランダムでなくて推定できる部分は?)という質問に答える.
  • 情報量の期待値をエントロピーHと定義する.
  • よって,離散系なら下記で計算される.
H(x) = sum{-logP(x) * P(x)}   (2)
  • エントロピーはP(x)が一様分布の時に最大になる.一様分布な時,つまり完全ランダムな時に最大であるということ.

相互情報量

  • 2つの情報源X,Yがあった時に,その情報源の互いの情報量は?という疑問に答える.つまり,一方の情報源Xの情報を知った時に,Yの情報をどの程度わかるかな?という問い.
  • 確率の独立性と関係があるのが直感的にもわかる.
  • 計算方法は下記.ここで,H(X,Y)はXとYの同時事象の情報量の期待値(結合エントロピーと言う.XとYの両方から知れること).式から想像できるように,それぞれの独立の時のエントロピーから,重複部分を除く,というよくあるド・モルガンの図のイメージ.
I(X,Y) = H(X) + H(Y) - H(X,Y)   (3)

それぞれのエントロピーの関係

  • ド・モルガン図でイメージすると分かりやすい.
  • 事象Xのエントロピーと事象Yのエントロピーがそれぞれあるとする.
  • 相互情報量は重なっているA&Bの部分.ここが,Xを知ってYを知れることだし,Yを知ってXを知れることだもんね.
  • 結合エントロピーはA|Bなエントロピー.つまり,重複も許して,XとYから知れること.
  • 紛らわしいのが条件付きエントロピー.意味は,Xを知っている状況でYから知れること.これは包含図でいうと,H(X)-H(X&Y)の部分.つまり,重複部分を除いた,純粋にXだけからしか知れないこと.
  • もう一つ忘れてならないのがKLダイバージェンス.2つの分布の距離を測るものらしい,が式から直感的に理解ができない.2つの分布の比の一方からみた期待値,という感じなのかな?
  • KLDとMIの関係で悩むけれど,よく分からない・・・.MIは2つの確率変数X,Yに対する定義でKLDは1つの確率変数上の2つの分布で書かれているのを目にするけど,KLDはそれに限定されているのかな?もっと勉強が必要だ.

オブジェクトの配列

概要

自分で作ったオブジェクトの配列を作りたいことはよくある.特に,vectorで管理したいことはよくある.その時に,不要なコピーコンストラクタを呼ばないようにc++11以降ではpush_backの代わりにemplace_backを使う. emplace_backの引数はコンストラクタの引数として与えられる.引数を一つも取らない場合には空で実行する.

#include <iostream>
#include <string>
#include <vector>

class Person
{
private:
  std::string name;
public:
  Person() {
    cout << "default constructor" << endl;    
  };
  Person(std::string _name) {
    name = _name;
  };
  void hello() {std::cout << "I'm " << name << std::endl;};
};

int main()
{
  std::vector<std::string> names{"Nobunaga", "Ieyasu", "Hideyosi"};
  std::vector<Person> people;
  for (auto iter=names.begin(); iter != names.end(); ++iter)
    people.emplace_back(*iter);
  
  for (auto iter=people.begin(); iter != people.end(); ++iter)
    iter->hello();
  return 0;
}

ちなみに

ひさしぶりにC++を書こうとすると,いつコンストラクタが呼ばれるのかをいつも忘れて悩んでしまう.vectorで個数を指定して宣言する場合,初期化も行われる.このときは暗黙にデフォルトコンストラクタ(引数無しのコンストラクタ)が呼ばれる.

#include <iostream>
#include <string>
#include <vector>

class Person
{
private:
  std::string name;
public:
  Person() {
    cout << "default constructor" << endl;    
  };
  Person(std::string _name) {
    name = _name;
  };
  void hello() {std::cout << "I'm " << name << std::endl;};
};

int main()
{
  std::vector<std::string> names{"Nobunaga", "Ieyasu", "Hideyosi"};
  // それぞれの要素の初期化でコンストラクタが呼ばれる.この場合は3回.
  std::vector<Person> people(names.size());
 
  return 0;
}

Emacsの設定

Font設定

  • RIcty-Diminisedを別途インストールしておく。
  • Option-Set Default FontからGUIでフォントを設定する。
  • scratchバッファで(frame-parameter nil 'font)でフォント名を確認する。
  • その文字列を下記のように追記する。
(add-to-list 'default-frame-alist '(font . "-PfEd-Ricty Diminished-normal-normal-normal-*-15-*-*-*-*-0-iso10646-1"))

基本設定

(package-initialize)
(setq package-archives
      '(("gnu" . "http://elpa.gnu.org/packages/")
        ("melpa" . "http://melpa.org/packages/")
        ("org" . "http://orgmode.org/elpa/")))

;; Preference & Misc settings
(load-theme 'deeper-blue t)
(tool-bar-mode 0)
(setq inhibit-startup-message t)
(setq make-backup-files nil)
(setq ring-bell-function 'ignore)
(show-paren-mode t)
(setq frame-title-format "%f")
(column-number-mode t)
(line-number-mode t)

;; Key bindings
(define-key global-map "\C-h" 'delete-backward-char)
(define-key global-map "\C-z" 'undo)
(global-set-key [C-tab] 'next-multiframe-window)
(global-set-key [C-S-iso-lefttab] 'previous-multiframe-window)

;; font
(add-to-list 'default-frame-alist '(font . "-PfEd-Ricty Diminished-normal-normal-normal-*-15-*-*-*-*-0-iso10646-1"))

;; Install use-package.el
(when (not (package-installed-p 'use-package))
  (package-refresh-contents)
  (package-install 'use-package))

拡張パッケージ

(Caskはプロキシ周りでハマるのでパッケージ管理は,package.elとuse-package.elに変更した。)

補完関係(company, irony(c/c++), jedi(python))

  • 補間がauto-completeよりもcompanyの例を見るのでそっちにする.特に,C/C++ようにirony-modeも入れておく.Pythonにはcompany-jediを入れておく(jediは入れない).
  • Caskでcompany, company-ironyを指定しておく.
  • 初回起動時にはirony-install-server(c/c++), jedi:install-server(Python)を実行する.(この時にlibClangが無いとエラーになるので,libClang-devを入れておく.Macだとllvm --with-clangで入っている。Cmakeも必要なので注意。)
  • elpyのコード内移動(定義元へジャンプなど)は内部ではjediを使っている。よって,elpyでコード内移動がうまく動いていないときはjediの設定が間違っている。M-x elpy-configでjediがちゃんと見つかっているかチェックする。また,Ubuntu18.04はpythonがpython2がまだデフォルトなので,RPCもpython3を指定するようにする。

スニペット

  • 最低限の設定は下記。最近のバージョンはTABが展開キーだと指定する必要がある。
  • 下記の設定だと標準のyas-installed-snippet-dir内のsnippetと、自分で追加したsnippet(.emacs.d/snippet)が使用される。
  • 使い方としては、(キーバインドを変えているけど)、例えばpythonモードで、classでTABを打つとクラス定義のテンプレートを生成してくれるし、ifmでTABを押すとif name == 'main':を展開してくる。
  • 新しくsnippetを生成するときは、'C-x i i'でYasnippet編集バッファに移動して、名前(ファイル名)と、キー(TABで展開する指定文字)を決めて、その後ろに展開したいコードを書く。C-c, C-cで保存する。モードとか、フォルダを聞いてくるので適時応える。

dotemacs (use-package)

(require 'yasnippet)
(setq yas-snippet-dirs
      '("~/.emacs.d/snippets"
    yas-installed-snippets-dir))
(define-key yas-minor-mode-map (kbd "C-x i i") 'yas-insert-snippet)
(define-key yas-minor-mode-map (kbd "C-x i n") 'yas-new-snippet)
(define-key yas-minor-mode-map (kbd "C-x i v") 'yas-visit-snippet-file)
(yas-global-mode 1)
(custom-set-variables '(yas-trigger-key "TAB"))

;; Company (except for Python mode)
(use-package company
  :commands global-company-mode
  :init (progn
          (global-company-mode)
          (setq company-global-modes '(not python-mode cython-mode sage-mode))
          )
  :config (progn
            (setq company-tooltip-limit 20) ; bigger popup window
            (setq company-idle-delay .3)
            (setq company-minimum-prefix-length 3)
            (setq company-echo-delay 0)     ; remove annoying blinking
            (setq company-begin-commands '(self-insert-command))
            ))
;;   ;(define-key company-mode-map (kbd "TAB") 'company-complete)
;;   (define-key company-active-map (kbd "TAB") 'company-complete-common)
;;   (define-key company-active-map (kbd "M-h") 'company-show-doc-buffer))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Python
;; Python
(use-package flycheck
  :ensure t
  :init
  (global-flycheck-mode t)
  :config
  (setq flycheck-python-flake8-executable "flake8"))

;; Standard Jedi.el setting
(use-package jedi
  :ensure t
  :init
  (add-hook 'python-mode-hook 'jedi:setup)
  (setq jedi:complete-on-dot t)
  :config
  (use-package company-jedi
   :ensure t
   :init
   (add-hook 'python-mode-hook (lambda () (add-to-list 'company-backends 'company-jedi)))
   (setq  company-jedi-python-bin "python")))

(use-package elpy
  :ensure t
  :init (with-eval-after-load `python (elpy-enable))
  :config
  (progn
    ;; Use Flycheck instead of Flymake
    (when (require 'flycheck nil t)
      (remove-hook 'elpy-modules 'elpy-module-flymake)
      (remove-hook 'elpy-modules 'elpy-module-yasnippet)
      (remove-hook 'elpy-mode-hook 'elpy-module-highlight-indentation)
      (add-hook 'elpy-mode-hook 'flycheck-mode))
    (elpy-enable)
    ;; jedi is great
    (setq elpy-rpc-python-command "python3")
    (setq python-shell-interpreter "ipython3")
    (setq python-shell-interpreter-args "-i --simple-prompt")
    (setq elpy-rpc-backend "jedi")))

Cask編(もう使っていないけどメモとして残しておく)

package.elではなく,caskかEl-docというので管理するのが主流らしい. caskを使ってみる.

インストール

$ curl -fsSkL https://raw.github.com/cask/cask/master/go | python # ubuntu
$ brew install cask # OSX
$ export PATH="/home/nobunaga/.cask/bin:$PATH"

初期設定と使い方

$cd .emacs.d
$ cask init    # 初期化時に一度だけ実行.Caskというファイルが生成される.
$ cask install # Caskの指定にしたがって必要なパッケージをインストール.
$ cask update  # updateがあれば更新する.

Caskの書き方

  • (depends-on "package-name") でpackage-nameのパッケージが.emacs.d/.cask以下に保存される

init.elの設定

  • 下記を記入しておくと,Caskに書いたファイルを自動でrequireしてくれる.
(require 'cask "~/.cask/cask.el") ; Ubuntu
(require 'cask)                   ; OSX
(cask-initialize)
(global-company-mode 1)
;; c/c++
(add-hook 'c-mode-hook 'irony-mode)
(add-hook 'c++-mode-hook 'irony-mode)
(add-hook 'irony-mode-hook 'irony-cdb-autosetup-comile-options)
(add-to-list 'company-backends 'company-irony)
;; python
(require 'jedi-core)
(setq jedi:complete-on-dot t)
(setq jedi:use-shortcuts t)
(add-hook 'python-mode-hook 'jedi:setup)
(add-to-list 'company-backends 'company-jedi)
;; general
(setq comany-idle-delay 0)
(setq comany-minimum-prefix-length 2)
(setq comany-selection-wrap-around t)
;; (setq company-idle-delay nil) ; 自動補完をしない
(global-set-key (kbd "C-M-i") 'company-complete)
(define-key company-active-map (kbd "C-n") 'company-select-next)
(define-key company-active-map (kbd "C-p") 'company-select-previous)
(define-key company-search-map (kbd "C-n") 'company-select-next)
(define-key company-search-map (kbd "C-p") 'company-select-previous)
(define-key company-active-map (kbd "<tab>") 'company-complete-selection)

Yasnippet

(require 'yasnippet)
(setq yas-snippet-dirs
      '("~/.emacs.d/snippets"
    yas-installed-snippets-dir))
(define-key yas-minor-mode-map (kbd "C-x i i") 'yas-insert-snippet)
(define-key yas-minor-mode-map (kbd "C-x i n") 'yas-new-snippet)
(define-key yas-minor-mode-map (kbd "C-x i v") 'yas-visit-snippet-file)
(yas-global-mode 1)
(custom-set-variables '(yas-trigger-key "TAB"))

ODEのテンプレートプログラム

ODEのテンプレートプログラムで流れを確認しておく.

ヘッダや宣言

dobuleとfloatの描画ルーチンを自動でコンパチにする宣言を入れておく.

#include <iostream>
#include <ode/ode.h>
#include <drawstuff/drawstuff.h>

// Doubleで扱う場合は描画ルーチンをD型にする.
#ifdef dDOUBLE
#define dsDrawBox      dsDrawBoxD
#define dsDrawSphere   dsDrawSphereD
#define dsDrawCylinder dsDrawCylinderD
#define dsDrawCapsule  dsDrawCapsuleD
#define dsDrawLine     dsDrawLineD
#endif

// Drawstuffで描画する際のテクスチャの場所を指定する.
#define TEXTURES_PATH ("/usr/local/share/drawstuff/textures")
#define DENSITY (5.0)

オブジェクトの構造体

剛体計算はBodyとMassで,衝突計算はGeomを使う.両方を使う場合はそれで一つのオブジェクトになるように構造体を宣言しておくと便利.

各種オブジェクトの描画ルーチンセット

これは問題関係なく共通で使えると思う.GeomからPosition/Rotationとかの情報は取れるので,GeomIDを渡すようにする. でも,実は下のはちゃんと動かない...シリンダーとカプセルのIDがどっちも3でうまくスイッチできない...なぜだろう?

// Both of the Cylinder and Capsule class id is 3 ??
void drawObject(dGeomID gID,float red,float green,float blue) {
  dsSetColor(red,green,blue);
  //std::cout << dGeomGetClass(gID) << std::endl;
  switch (dGeomGetClass(gID)) {
  case dSphereClass:
    {
      dsDrawSphere(dGeomGetPosition(gID),dGeomGetRotation(gID),dGeomSphereGetRadius(gID));
      break;
    }
  case dBoxClass:
    {
      dVector3 length;
      dGeomBoxGetLengths(gID, length);
      dsDrawBox(dGeomGetPosition(gID),dGeomGetRotation(gID),length);
      break;
    }
  case dCylinderClass:
    {
      dReal r,length;
      dGeomCylinderGetParams(gID, &r,&length);
      dsDrawCylinder(dGeomGetPosition(gID),dGeomGetRotation(gID),length,r);
      break;
    }
    case dCapsuleClass:
    {
      dReal r,length;
      dGeomCapsuleGetParams(gID, &r,&length);
      dsDrawCapsule(dGeomGetPosition(gID),dGeomGetRotation(gID),length,r);
      break;
    }
  }
}

Drawstuffのセットアップ

Drawstuffを管理するオブジェクトdsfFunction(ここではfnとしている)を定義しておいて,下記でセットアップする.

void  setDrawStuff() {
  fn.version = DS_VERSION;
  fn.start   = &initDs;     //初期化関数へのポインタ
  fn.step    = &simLoop;   // 無限ループで毎回実行する関数ポインタ
  fn.command = command;    // キーボード入力を処理する関数ポインタ
  fn.stop    = NULL;       // 停止時の処理の関数ポインタ
  fn.path_to_textures = TEXTURES_PATH;
}

dsFuntion.startにセットする関数ではカメラ視点の設定などを実行する.xyzはカメラの絶対座標,hprはヨー,ピッチ,ロール.右手系にHPRのそれぞれの回転角度があるとして,それぞれの方向にどれだけ回すかを与える.注意点として,DrawStuffの角度は度(degree)で与えられる.他は一般的にラジアンで統一されている場合があるので注意. 下の例だとX軸方向に180度,つまり画面の奥側を向いて,Y,Z軸には回転しない.

static void start() {
  static float xyz[3] = {10, 0, 1.0};
  static float hpr[3] = {-180, 0, 0};
  dsSetViewpoint (xyz,hpr);          
}

キーボード入力を処理する関数は例えば下記.

void command(int cmd) {
  switch (cmd) {
  case 'r': v += 0.5; break;
  case 'l': v -= 0.5;break;
  default: std::cout << "WARNING: unknown key input" << std::endl;
  }
}

fn.stepはDrawstuffの中の無限ループ内で実行される関数をセットする.逆に,Drawstuffを使わないで(描画しないで)ODEのシミュレーションをする場合は,この関数を別途無限ループの中に入れるんだ. ジョイントを制御をする場合のcontrolに関しては,その項目を参照.

static void simLoop (int pause) {
  if (!pause) {
    dSpaceCollide(space,0,&nearCallback); // 衝突計算
    dWorldStep(world,0.1);                // タイムステップを1時刻すすめる(動力学計算)
    dJointGroupEmpty(contactgroup);       // 衝突計算のジョイントグループを初期化
    control();                            // ジョイントを制御する場合はその関数を呼び出す.
  }
  // 結果を描画する(Drawstuffを使わない場合は飛ばす)
  drawObject(base.geom, 0,0,05);

  // STDOUTに情報を出力する場合
  std::cout << "angle:" << dJointGetHingeAngle(joint);
}

衝突計算をする場合

衝突計算自体はdSpaceCollideで実行される.これは各オブジェクトが衝突するかどうかを判定してくれる.この関数にはnearCallBack関数を登録する.これは,2つの物体が衝突する(かも知れない)と判定された時に呼び出される関数.このコールバック関数の中で接触点をジョイントとして登録する.ジョイントはすなわち拘束条件だ.摩擦係数や反発係数は別途PRAM_BOUNCEなどでマクロ定義しておく.

static void nearCallback(void *data, dGeomID o1, dGeomID o2) {
  const int N = 10;
  dContact contact[N];

  int isGround = ((ground == o1) || (ground == o2));

  int n =  dCollide(o1,o2,N,&contact[0].geom,sizeof(dContact));
  if (isGround) {
    for (int i = 0; i < n; i++) {
      contact[i].surface.mode = dContactBounce|dContactSoftERP|dContactSoftCFM;
      contact[i].surface.bounce = PRAM_BOUNCE;
      contact[i].surface.mu = PRAM_MU;
      contact[i].surface.soft_erp = PARAM_SOFT_ERP;
      contact[i].surface.soft_cfm = PRAM_SOFT_CFM;
      dJointID c = dJointCreateContact(world,contactgroup,&contact[i]);
      dJointAttach (c,dGeomGetBody(contact[i].geom.g1),dGeomGetBody(contact[i].geom.g2));
    }
  }
}

ロボット定義

シミュレーション対象となるロボットの定義は,(1)剛体(Body)の定義,(2)それに質量(Mass)を付けて,(3)幾何(Geom)を付ける.Geomは衝突計算をしない場合は不要.これを愚直に書いて行くんだ. 例えば,直方体に円柱の車輪をつける例.ほとんどそのままだけど,円柱は回転させる必要があるので,回転行列rotを定義して,回転させたい軸と回転角に対する回転行列をdRFromAxisAndAngle([0,1,0]の軸の周りに,[0,1,pi/2]回す)で作っておいて,dBodySetRotationで回転する.  また,車輪を回したいのでジョイントを付ける.ODEのジョイントはモータを内蔵しているので,それを回すことで車輪を回せる.Axisがジョイントの軸方向,Anchorが軸の場所.

void makeRobot() {
  dMass m;     
  dMassSetZero (&m);

  //basement
  base.body = dBodyCreate (world);
  dMassSetBox (&m,DENSITY,sides[0],sides[1],sides[2]);
  dBodySetMass (base.body,&m);
  dBodySetPosition (base.body,0,0,radius+sides[2]/2);
  base.geom = dCreateBox(space, sides[0], sides[1], sides[2]);
  dGeomSetBody(base.geom,base.body);

  dMatrix3 rot;
  // wheel
  dRFromAxisAndAngle(rot,0,1,0, M_PI/2.0);
  wheel.body = dBodyCreate (world);
  dMassSetCylinder(&m,DENSITY,1,radius,width);
  dBodySetMass (wheel.body,&m);
  dBodySetPosition (wheel.body,0.0,0,0.3);
  wheel.geom = dCreateCylinder(space, radius, width);
  dGeomSetBody(wheel.geom,wheel.body);     
  dBodySetRotation(wheel.body,rot);

  joint = dJointCreateHinge(world, 0);
  dJointAttach(joint, base.body, wheel.body);
  dJointSetHingeAxis(joint, 1, 0, 0);
  dJointSetHingeAnchor(joint, 0.0, 0, 0.3);
}

Main関数

ここまで出来たら,最後にメイン関数を書いて終わりだ.world,spaceなどをグローバル変数で渡すような例が多いようなので,グローバル変数として扱うなら適切な場所で定義しておく.

// Global variables
static dWorldID      world;
static dSpaceID      space;
static dGeomID       ground;
static dJointGroupID contactgroup;
dsFunctions fn;

int main (int argc, char **argv)
{
  setDrawStuff();
  dInitODE();            
  world = dWorldCreate();
  dWorldSetGravity(world,0,0,-3.0);
  dWorldSetERP(world, 0.5);
  dWorldSetCFM(world, 1e-05);
  space = dHashSpaceCreate(0);
  contactgroup = dJointGroupCreate(0);
  ground = dCreatePlane(space,0,0,1,0);

  makeRobot();
  // initial noise
  dsSimulationLoop (argc,argv,600,400,&fn);
  dWorldDestroy (world);
  dJointGroupDestroy(contactgroup);
  dSpaceDestroy(space);
  dCloseODE();          

  return 0;
}

制御する

制御をするには力やトルクを与える.ここで重要なのは,力やトルクはシミュレーションの更新ごとに0に初期化されるので,継続的に力を加える場合は継続的に力を加える.一回与えても,その一回で消えちゃうよ. 力やトルクをdBodyAddForceのようにBodyに直接力を加える場合の他に,ODEのジョイントはモータになっていて,そのモータにトルクを与えることができる.実際にはジョイントにはトルクではなく速度を与える.

シミュレーション上のTips

描画をしない

デバッグが終わったら描画をせずにシミュレーションだけをしたい場合もある.その場合はsimLoopの描画の箇所をスキップするようにして,Drawstuffのループではなくて自分のWhile(0)ループの中でsimLoopを回す. 他にも,ctrl-tでテクスチャをOffにできる.ウィンドウサイズを変えたりでもシミュレーション速度が変わる.

制御コマンド

デフォルトでいくつかの制御キーコマンドが定義されている. - ctrl-t:テクスチャをOn/Off - ctrl-s:影のOn/Off - ctrl-p:シミュレーションの一時停止 - ctrl-o:1ステップ実行

リスタート

強化学習のシミュレーションなどをやる場合には何度もシミュレーションを繰り返す.やり方は,一旦Body, Geom, Joint,JointGroupを全部消して,再度作りなおす.

void makeRobot() {
... // 上のと同じでいい
}

void deleteRobot() {
  dJointDestory(joint); // jointを削除
  dBodyDestory(body);   // bodyを削除
  dGeomDestory(geom);   // geomを削除
  dJointGroupDestory(contactgroup);
}

void restart() {
  deleteRobot();
  contactgroup = dJointGroupCreate(0);
  makeRobot();
}

ランダムウォーク

お題

X軸上をランダムに動く点の軌跡を求める.1ステップで±1か,動かないか.

python

一歩ずつやってもいいけど,ガバッと乱数作ってnumpyのcumsumを使う方法でもいい.

import numpy as np

move = np.random.randint(-1,2,100) # generate 100 random integers in [-1,2).
history = move.cumsum()