gccのunordered_mapの実装を読んでみる

最近諸事情でアルゴリズムイントロダクションを読んでおります。

で、平均挿入・検索・削除時間がO(1)である謎のテクノロジー、ハッシュテーブルの章を読んだので、せっかくなのでSTLの実装を調べてみました。

なお、コンパイル時の条件によって色々実装が分岐するようですが、面倒臭いので適当に目についたコードを読んでいます。なので、組み合わせとして一緒に呼ばれないコードを読んでいる可能性がありますが、ご了承下さい。

とりあえず、挿入を足がかりに

std::unorderd_mapの実装は、いろいろたらい回しにされたあとbits/hashtable.hに記述されているらしいことが分かります。とりあえず挿入っぽい処理を調べてみると、こんな感じになります。

  template
    void
    _Hashtable<_Key, _Value, _Alloc, _ExtractKey, _Equal, 	       _H1, _H2, _Hash, _RehashPolicy, _Traits>::
    _M_insert_bucket_begin(size_type __bkt, __node_type* __node)
    {
      if (_M_buckets[__bkt])
	{
	  // Bucket is not empty, we just need to insert the new node
	  // after the bucket before begin.
	  __node->_M_nxt = _M_buckets[__bkt]->_M_nxt;
	  _M_buckets[__bkt]->_M_nxt = __node;
	}
      else
	{
	  // The bucket is empty, the new node is inserted at the
	  // beginning of the singly-linked list and the bucket will
	  // contain _M_before_begin pointer.
	  __node->_M_nxt = _M_before_begin._M_nxt;
	  _M_before_begin._M_nxt = __node;
	  if (__node->_M_nxt)
	    // We must update former begin bucket that is pointing to
	    // _M_before_begin.
	    _M_buckets[_M_bucket_index(__node->_M_next())] = __node;
	  _M_buckets[__bkt] = &_M_before_begin;
	}
    }

どうやら、衝突の解決には片方向リストを使った方法を使っているようです。アルゴリズムイントロダクションには「削除が必要な時は片方向リストを使った方法がいい」みたいなことが書いてあったような気がするのですが、そのとおりっぽいです。あと番兵は使ってないみたいですね。やっぱりメモリの節約が優先なのでしょうか。

バケットのインデックスは?

で、この__bktはどこから来てるのかを調べてみると、_M_bucket_indexという関数の値が実際には使われており…bits/hashtable_policy.hに次のようなコードがあります。(型の条件によって複数あるようですがこれが一番なんか読みやすそうなのでこれにしました)

      std::size_t
      _M_bucket_index(const __node_type* __p, std::size_t __n) const
	noexcept( noexcept(declval()(declval()))
		  && noexcept(declval()((__hash_code)0,
						    (std::size_t)0)) )
      { return _M_h2()(_M_h1()(_M_extract()(__p->_M_v())), __n); }

というわけで、h1に値を突っ込んでハッシュを得た後、さらにh2にそのハッシュとバケットの総数を入れてハッシュをバケット数以下にして、実際のバケットのインデックスを求めているみたいです。

__nはバケット数の最大値です。_Key&は例えばstd::unordered_map<string, int>ならstringの方です。_M_h1と_M_h2は実はbits/unorderd_map.hの最初の方にすでに書いてあって、

  template,
	   typename _Pred = std::equal_to<_Key>,
	   typename _Alloc = std::allocator<std::pair >,
	   typename _Tr = __umap_traits<__cache_default<_Key, _Hash>::value>>
    using __umap_hashtable = _Hashtable<_Key, std::pair,
                                        _Alloc, __detail::_Select1st,
				        _Pred, _Hash,
				        __detail::_Mod_range_hashing,
				        __detail::_Default_ranged_hash,
				        __detail::_Prime_rehash_policy, _Tr>;

h1の実体はstd::hash<_Key>(_key const& key)で、これは後で調べますが、とりあえず色々な型についてhash値を返してくれるSTLの公開の標準ライブラリです(皆様もお使いになれます!)

h2は__detail::_Mod_range_hashingで固定で、ユーザーが指定したりは出来ないようです。

h2はただのmod

_Mod_Range_hashingはbits/hashtable_policy.hに実体があり、

  // Many of class template _Hashtable's template parameters are policy
  // classes.  These are defaults for the policies.

  /// Default range hashing function: use division to fold a large number
  /// into the range [0, N).
  struct _Mod_range_hashing
  {
    typedef std::size_t first_argument_type;
    typedef std::size_t second_argument_type;
    typedef std::size_t result_type;

    result_type
    operator()(first_argument_type __num,
	       second_argument_type __den) const noexcept
    { return __num % __den; }
  };

というわけで、散々たらい回しにされた挙句ただのmodでした。万能ハッシュとか使ってるのかと思ったのですがそんなことはなく、h1の実体であるstd::hashに衝突の回避とかの仕事はほぼ丸投げして、バケット数とmodをとっているという感じになります。

std::hashの実体

std::hashの実体は色々なところに散らばっているのですが、std::stringの実装がやっぱり一番気になるのでそれを調べてみると、basic_string.hに実体があって、

  /// std::hash specialization for string.
  template<>
    struct hash
    : public __hash_base
    {
      size_t
      operator()(const string& __s) const noexcept
      { return std::_Hash_impl::hash(__s.data(), __s.length()); }
    };

  template<>
    struct __is_fast_hash> : std::false_type
    { };

となっております。ほかのstd::wstringとかutf16とかも全部同じ実装です。hashの先を見れば分かるのですが、この関数はvoid*を取る関数なので、文字列としての特性はとくに使わず、バイト列として解釈してhashを計算しているようです。

さて、_Hash_impl::hashはfunctional_hash.hの関数です。このseedの値の意味がよくわかりません…。

    static size_t
    hash(const void* __ptr, size_t __clength,
	 size_t __seed = static_cast(0xc70f6907UL))
    { return _Hash_bytes(__ptr, __clength, __seed); }

なんとこの先はlibsupc++という別ライブラリに投げれていて、

  inline std::size_t
  unaligned_load(const char* p)
  {
    std::size_t result;
    __builtin_memcpy(&result, p, sizeof(result));
    return result;
  }

}

namespace std
{
_GLIBCXX_BEGIN_NAMESPACE_VERSION

#if __SIZEOF_SIZE_T__ == 4

  // Implementation of Murmur hash for 32-bit size_t.
  size_t
  _Hash_bytes(const void* ptr, size_t len, size_t seed)
  {
    const size_t m = 0x5bd1e995;
    size_t hash = seed ^ len;
    const char* buf = static_cast(ptr);

    // Mix 4 bytes at a time into the hash.
    while(len >= 4)
      {
	size_t k = unaligned_load(buf);
	k *= m;
	k ^= k >> 24;
	k *= m;
	hash *= m;
	hash ^= k;
	buf += 4;
	len -= 4;
      }

    // Handle the last few bytes of the input array.
    switch(len)
      {
      case 3:
	hash ^= static_cast(buf[2]) << 16;
      case 2:
	hash ^= static_cast(buf[1]) << 8;
      case 1:
	hash ^= static_cast(buf[0]);
	hash *= m;
      };

    // Do a few final mixes of the hash.
    hash ^= hash >> 13;
    hash *= m;
    hash ^= hash >> 15;
    return hash;
  }

この関数もなんとなくバイト列中のたくさんのビットを反映させていることはなんとなく分かるのですが、謎のマジックワードもあるし、謎です。

やっぱり本と実装はすこし違いますね~とおもったのでした。おしまい。

他の型のハッシュは?

同じfunctional_hash.hにマクロがあって、

  // Explicit specializations for integer types.
#define _Cxx_hashtable_define_trivial_hash(_Tp) 	\
  template<>						\
    struct hash<_Tp> : public __hash_base  \
    {                                                   \
      size_t                                            \
      operator()(_Tp __val) const noexcept              \
      { return static_cast(__val); }            \
    };

  /// Explicit specialization for bool.
  _Cxx_hashtable_define_trivial_hash(bool)

  /// Explicit specialization for unsigned long long.
  _Cxx_hashtable_define_trivial_hash(unsigned long long)

という感じです。なんとstatic_castするだけ!型に収まってるならともかく、unsigned long longまでstatic_castしちゃうのは本当にそれでいいの?感があります。…ということは32bit環境下でunordered_mapに下32ビットが同じunsigned long longを送り続けると死ぬほど遅い可能性がある…? →実際に実験したら遅かった

ちなみにこのコードのすぐ下にありますが、double/floatに関してはバイト列として文字列と同じハッシュ関数を使っています。うーん…(^_^;)

bashで3Dプログラミング

今日はbashと基本的なシェルコマンドだけで3Dレンダリング(ワイヤフレーム)をしてみました!

Github: ledyba / Bash3D

bashプログラミングテクニック

今回つかったbashのプログラミングテクニックについて、記録のついでに少しまとめます。シェルスクリプトの基本的な所は書きません。

あと、3Dプログラミングに関しては本に書いてある事以上のことが書けそうにないので、(今は)書きません。

bashのデータ構造

bashでそこそこ本格的にプログラミングする事を考えた上でまっさきに「辛そう…」ってなるのはデータ構造ですよね!今回のソフトでは、ベクトルや行列を綺麗に扱えないと困ってしまうでしょう。まずはそこを見ます。

bashの「値」は基本的に全てただの文字列ですが、「リスト構造」にだけはネイティブで対応しています。こんな感じ。

$ lst=(1 2 3)
$ echo ${lst[0]}
1
$ echo ${lst[1]}
2
$ echo ${lst[2]}
3
$ echo ${lst[3]}
#空行

データにアクセスする時に、必ず${}でくくらないといけないのが少し面倒です。ちなみに添字を付けないとリスト全体…と思いたいところですが、最初の要素が出てきます。

$ echo $lst
1

全体を表示するには、${lst[@]}とします。

$ echo ${lst[@]}
1 2 3

リストを定義するときに使ったカッコが出てこないので、リストをコピーしたい場合はそのカッコを補わないといけません。

$ lst2=(${lst[@]})
$ echo ${lst2[@]}
1 2 3

慣れればなんとかなるレベルですが、面倒です。

さらに気を付けないといけないのは、”(1 2 3)”という文字列を返す関数を$(command)で受けて変数に代入してもリストになりません。

$ function fun() { echo \(1 2 3\); }
$ lst=$(fun)
$ echo $lst
(1 2 3)

ただの文字列として代入されます。どうしてもしたい場合は、evalを使わないといけません。

#文字列展開が起こってから実行されるので、lst=(1 2 3)というコマンドが実行されることになる
$ eval "lst=$(fun)"
$ echo ${lst[1]}
2
$ echo ${lst[@]}
1 2 3

通常の用途ならそれで済むのですが、今回は出来る限り高速に処理したいので、こうやってコマンドが沢山必要になるのは避けたいところです。

ところで。

後に言うダイナミックスコープとevalを組み合わせると連想配列みたいなことも出来るのですが、今回のプログラムでは基本的にこのリストでデータ構造を管理しています。

インデックスを指定したアクセスをするためにコストが高いevalを行わなくてよいこと、ベクトルと行列が主役なのでリストで十分な事が理由です。

$ matrix=(MAT 1 0 0 0 0 1 .....)
$ vector=(VEC 0 0 0 1)

最初の要素に型の名前を入れて適宜assertすることで、動的型検査みたいなことをしています。

変数のスコープ

名前空間は基本的にひとつしかありません。関数のネストは一切関係なく上書きされていきます。

$ val=0
$ function fun() { val=1; }
$ fun
$ echo $val
1 #関数内の代入だが、グローバル領域が上書きされてしまった
$ function fun2() { val2=2;echo $val2; }
$ fun2
2
$ echo $val2
2 #オフサイドルールではなく、存在しない変数はグローバルに作られるみたい

一種のダイナミックスコープみたいなものだと理解しておけば大体書けます。

呼び出し元の変数がうっかり書き換えられてしまう可能性があるのですごく不便ですが、逆に利用することも出来ます。

「関数の引数に変数名を渡すと、evalと組み合わせることで指定した変数に結果を格納する関数」が作れます。こんな感じ:

$function setOne() { eval "$1=1"; }
$ setOne val1
$ echo $val1
1
$ setOne val2
$ echo $val2
1

これを使って「最初の引数は入力、最後の引数は出力」と規約を決めて関数を設計しています。たとえば行列の計算はこちらにあります。見た目がアセンブラみたいになりますネ。

「こんなスコープ再帰が出来ないのでは?」という気がしますが、案外なんとかなります。というのも、$(command)はサブシェルとかいうので実行されるらしくて、さらにサブシェルは別の変数空間で動作するようです。つまり、このスコープの問題は起こりません。…たぶん。よく分からない(◞‸◟ )

 #再帰関数といえば階乗である
function kaijo() {
    arg=$1
    next=$(expr $arg - 1)
    if [ $arg -gt 1 ]; then
        result=`kaijo $next`
        echo $(expr $arg \* $result)
    else
        echo 1
    fi
}

kaijo 10 --> 3628800

ついでにevalとこのスコープの仕様を使うと、擬似的な連想配列を作れますが、値のコピーが難しいこともあって使いませんでした。

$ val_key1=123
$ val_key2=456
$ function showKey1() { echo $(eval "echo \$$1_key1"); };
$ showKey1 val
123

固定小数点とexprコマンド

シェルスクリプトで計算をするには、exprコマンドを使う必要があります。このコマンドは整数の計算にしか対応していません。

$ expr 1 \* 2 \* 3
6
$ expr 1 \* 2 \* 3.2
expr: non-numeric argument

3Dソフトでそれでは困るので、整数だけで小数を計算します。それにはどうするかというと、扱う小数を適当に1000倍とかして、それを小数に見立てます。つまり、1000は1、1001は1.001、3140は3.140だと思うことにします。その場合、計算する時にちょっとした手間が必要になります。

#3.14 + 3.14
$ expr 3140 + 3140 #足し算・引き算はそのまま足し算で大丈夫です。
6280
#3.14 / 3.14
$ expr 3140 \* 1000 / 3140 #割り算するときは、倍率で掛けてから割ります。この世界では「1.0」は1000であることに注意
1000
#3.14 * 3.14
$ expr 3140 \* 3140 / 1000 #掛け算するときは、掛け算してから倍率で割ります。
9859

倍率を1000じゃなくてもっと高くすれば、もっと小数演算の結果が正確になります。

私の環境だけかもしれませんが、exprはどんな大きな数でもちゃんと扱ってくれるようなので、ガンガン倍率を上げていけば誤差がまったく気にならなくなります。これはすごい楽…(楽すぎた)

なお、sin/cosは(テイラー展開を自分で展開しない限り)計算できないようだったので、0~360度について事前に求めておいて、これを先ほどのリスト構造にいれて使っています。平方根も使いたかったんですが、どうやっても無理そうなので諦めました。

表示:オフスクリーンレンダリング

もちろんターミナルにはピクセルがないので、ターミナル上に表示する文字で代用します。普通に表示するなら、tput cup (line) (col)を使うとターミナル上の任意の位置に移動できるので、これを使いたいところです。

# 斜め45度に線を引いてみた
for1;do
    for2;do
        tput cup y x
        echo -n \*
    done
done

が、この方法では線を引く様子が目で分かるぐらい時間がかかります。これで最初レンダリングしたみたのですが、レンダリング中に線を引いている様子がバッチリ見えてしまい、あんまり回転アニメーションしてるように見えませんでした。

ほんもののゲームでもそういう、更新途中の絵が表示されてしまってちらつく問題があります。ほんもののゲームでは、そういう時は、書いてもすぐには表示されない「裏の画面」を用意しておいて、そちらに時間を掛けてレンダリングして、最後に表示される「表の画面」と入れ替える「オフスクリーンレンダリング」と「ダブルバッファリング」をしているのでした。私も同様にそのテクニックを使いました。

具体的には、こんな感じ。

SCR_BUFF=$(head -c $size < /dev/zero | tr '\0' '_')

今回の「ピクセル」は文字ですから、「裏の画面」はただの文字列になります。
画面のサイズを取得しておいて、そのサイズ分何もない文字列を意味する「_」をtrを使って入れておきます。空白でいいじゃないかと思うかもしれませんが、シェルスクリプトでは空白はかなり重要な意味を持っていて面倒なのでこうしています。

「ドットを打つ」時は、対応する文字列中のインデックスの文字を*に書き換えます。

SCR_BUFF=$(echo $SCR_BUFF | sed s/./*/<置き換える文字の位置>)

zshだとstr[12]=\*ってやると文字列を更新できるのですが、bashでは使えなかったので、sedを使って書き換えています。この処理が非常に重く、今回のレンダリングエンジン(?)の処理時間の殆どはこのドットを打つ処理に費やされています。上手い方法を知ってたら教えてください!

最後に表の画面に更新する方法ですが、単純にechoするだけです。この時、何も入っていない所には_が入っていたので、空白に直します。

$ tput cup 0 0
$ echo -n $SCR_BUFF | tr '_' ' '

ちゃんと幅に合わせて改行されるので、うまく表示できます。

高速化

インライン展開 → プロセスforkを抑える

当初綺麗に関数を使って書いた所、行列の掛け算に1秒近く掛かってしまい、さすがにレンダリングは出来ないのかな…と諦めかけたのですが、行列やベクトルの掛け算や足し算を全てインライン展開したところ、座標計算に関しては十二分すぎる程の速度が出せるようになりました。前述した通り、今回遅いのはオフスクリーンにピクセル(文字)を書き込む処理が異常に時間が掛かっているからです。

やはりプロセスの起動はものすごい重い処理のようで、可読性をそこまで損なわない範囲でプロセスの呼び出しを一箇所にまとめています。

感想

線を引くのにえげつない時間が掛かるのを見た時、初めて触ったBASICの処理系を思い出しました。懐かしい。

ところで、なんで3Dライブラリを作るとOpenGL風の名前にしちゃうんだろ。

ニコ動に「次はシェーディングだな!」というコメントがついているのですが、テクスチャの表示をやってみたいです!7色くらいしか使えないので、たいへんそう…。影は「暗い色」が基本的に無いのでかなりつらい気がします。

  1. x=0;x<10;++x []
  2. y=0;y<10;++y []

3Dをゼロからレンダリングするテスト

C++から動画を直接出力する自作ライブラリ「Dolly」の上で、フルスクラッチで3Dのレンダリングを行なっています。

とりあえずは、ワイヤーフレームの正方形をぐるぐる。

テクスチャもなんとか実装しましたが、三角形の継ぎ目が出ているのが気になります。でもどこが悪いのかわからん…

このサンプルのソースコードはこちらです:

https://github.com/ledyba/Dolly/blob/master/src/sample/super_px.cpp

東北地方太平洋沖地震の地震波を音にしてみた

地球物理をやる学科に居るので、地震も学んでいるのですが、ふと思いついたので地震計の観測した波形を音に見立ててWAVファイルにしてみました。

気象庁の強震波形(平成23年(2011年)東北地方太平洋沖地震)のページから、福島県郡山市朝日の、震度6のデータです。なんだか、雷みたいな音がしますね!

ソースコードはこちらです。PythonでCSV読んでWAVに書いているだけですが、一応使い方のサンプルとしては使えるのではないでしょうか。

凪のあすからの恋愛関係からページランクを求めてみました

最近友人に薦められて、今日が最終回らしい「凪のあすから」を見ています。今6話なのですが、恋愛関係が大変複雑です。「深夜の昼ドラ? アニメ『凪のあすから』の人物相関図つくったよ」というページがあったので、引用しますと、こんな感じ。

うーん、グラフ理論って感じだ。そう思ってからピクシブの非公式カップリングのページを見てしまうと、もはや隣接行列か何かにしか見えなくなってきた…。

20130403隣接行列…グラフ理論…はっ!ページランク求めようページランク!

というわけで、今まで一回ぐらいやってみようと思っていた、ページランクの計算を実際にやってみました。こちらの「Google の秘密 – PageRank 徹底解説」を参考にしました。人間をページに見立てて、好意をリンクに見立てます。

ページランクとは結局何なのか

昔私が読んだ本だと、

  • 全てのページはページランク値PRを持っていて、
  • 他のリンクしているページそれぞれににPR/N(Nはページ内のリンク総数)を与えており、
  • ページランクの値PRはリンクされた他のページから貰ったページランク値の和に等しい

というのを読んだ気がします。つまり、貰ったページランク=あげたページランク=そのページのページランクという等式が成り立つということです。

元の論文にもそんな図が載っています。

20130403-2

Lawrence Page, Sergey Brin, Rajeev Motwani, Terry Winograd, ‘The PageRank Citation Ranking: Bringing Order to the Web’, 1998,
http://www-db.stanford.edu/~backrub/pageranksub.ps

これをグラフ理論と確率過程の話で表現すると、こうなります。

  • 各ページから他のページへ移動する確率の表を用意する(とりあえずすべてのリンクで等確率=1/リンク総数とする)
  • その表に従ってあるページから出発してリンクをたどり続ける
  • すると、どのページから出発しても、N回リンクをたどった時にあるページPに居る確率がN+1回目のそれと変わらなくなってくる(これを定常状態と呼びます)
  • その時の、N回目でそれぞれのページに居る確率がページランクです。

これがホントに元の話と一緒なのか?ということなのですが、

  • “もらうページランク”=(自分も含めた他のページにN−1回目に居る確率*それぞれの遷移確率)の和=N回目にそのページにいる確率
  • “あげるページランク”=(N回目にそのページにいる確率 * そのページから他のページへの遷移確率)の和=N回目にそのページにいる確率1

なので、もらうページランク=あげるページランクという関係は確かに成り立ちます。

求めてみる

ページ遷移確率表=行列

実際にもとめてみましょう。まずは、鍵となるページ遷移確率表を用意します。さっきの表を見ながら、ひーくんからまなかへ遷移する確率は1…みたいな感じで行列を作るとこうなります。一応、今回は友情は無視して恋愛関係だけ考えます。

M =
 ひーくん ちさき  まなか   要    さゆ   美海   つむぐ  あかり  至    みをり ←好き ↓好かれる
   0.00000   1.00000   0.00000   0.00000   0.00000   1.00000   0.10000   0.00000   0.00000   0.00000 ひーくん
   0.00000   0.00000   0.00000   1.00000   0.00000   0.00000   0.10000   0.00000   0.00000   0.00000 ちさき
   1.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.10000   0.00000   0.00000   0.00000 まなか
   0.00000   0.00000   0.00000   0.00000   1.00000   0.00000   0.10000   0.00000   0.00000   0.00000 要
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.10000   0.00000   0.00000   0.00000 さゆ
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.10000   0.00000   0.00000   0.00000 美海
   0.00000   0.00000   1.00000   0.00000   0.00000   0.00000   0.10000   0.00000   0.00000   0.00000 つむぐ
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.10000   0.00000   0.50000   0.00000 あかり
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.10000   1.00000   0.00000   1.00000 至
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.10000   0.00000   0.50000   0.00000 みをり

普通、ページには複数のリンクがありますが、このアニメだと思い人が二人以上居るのは不倫と勘違いされて殴られた至だけです。

あと6話時点ではつむぐは誰も好きではないので、つむぐの行が全部0になるはずなのですが、そうしてしまうとページランクが求まりません2

これを元論文ではdangling pageと呼んでいます。元論文では「あんまり重要じゃないから削除して計算しています」と無視してしまっているのですが、つむぐくんを無視するのはかわいそう(?)なのと、仮につむぐくんを消してしまうと今度はまなかがdangling pageになって…となって最終的に大人組しか居なくなってしまうので、次に紹介するrandom surferモデルに習って、みんなを平等に好きなことにしました。

計算は掛け算するだけ

元の定義通り、「あるページからスタートしてリンクをたどり続けた時に、それぞれのページに居る確率が一定になってくる」まで遷移させ続けます。

たとえば、ひーくんからはじめましょう。はじめはひーくんに居ると仮定しているので、ひーくんに居る確率は1で、他みんな確率0とします。

1回リンクをたどると、グラフを見るとひーくんはまなかが好きなので、まなかの確率が1になるはず。

これを我々が先ほど作った行列を使って計算しようとするなら、先ほどの行列Mに、ひーくんを1としたベクトルPを掛ければよいです。

octaveでやるとこんな感じです。

octave:95> M=[ 0 1 0 0 0 1 1/10 0 0 0 ; 0 0 0 1 0 0 1/10 0 0 0; 1 0 0 0 0 0 1/10 0 0 0; 0 0 0 0 1 0 1/10 0 0 0 ; 0 0 0 0 0 0 1/10 0 0 0 ; 0 0 0 0 0 0 1/10 0 0 0; 0 0 1 0 0 0 1/10 0 0 0 ; 0 0 0 0 0 0 1/10 0 1/2 0; 0 0 0 0 0 0 1/10 1 0 1; 0 0 0 0 0 0 1/10 0 1/2 0 ]
(略)
octave:96> P = [1; 0; 0; 0; 0; 0; 0; 0; 0; 0] ←ひーくんに居る確率を1としたベクトル
octave:97> M*P
ans =

   0
   0
   1
   0
   0
   0
   0
   0
   0
   0

というわけで、まなかは三番目の要素で、それが1となったので、見事正しい結果を得られています。これを沢山続けてみましょう。10000回ぐらい掛けると変わらなくなってきます。

octave:98> (M^10000)*P
ans =

   0.00000
   0.00000
   0.00000
   0.00000
   0.00000
   0.00000
   0.00000
   0.22222
   0.55556
   0.22222

…あれぇ?最後の3人、つまり大人組以外ゼロ?なんか間違ったかなぁ?

ランダム・サーファー・モデル

こうなってしまった原因は、大人組から子供組へ繋がっていないことです。つむぐくんは全員が好きだと仮定したので、子供組から大人組へはつむぐくんから遷移することが出来るのですが、その後は大人組の間でぐるぐるぐるぐる巡り続けることになります。ひーくんからはじめても長い目で見ればいつかは大人組に入ってしまい、そこで定常状態となるので、子供組の確率がすべて0になってしまったのです。グラフ理論の言葉でいうと、グラフを辿ると元の位置に必ず戻ってこれる大人組のような部分を再帰類と呼び、そうではなく元に戻れない子供組のような部分を非再帰と呼ぶそうですが、あまりこれ自信をもって解説できないので、しません…。

元の論文ではこれをRank sinkと読んでいて、そのためにランダム・サーファー・モデルというのを考慮しています。これはどういうモデルかというと、「時々はまったくリンクされていない別のページに飛ぶ」と考慮しているモデルです。具体的には、新しい遷移確率行列を次のように定義します。

octave:99> M2=(M*1/1.15)+((1/10)*0.15/1.15)

これは、まず今までの確率を少し下げた上で、あるページから他のページへ僅かな確率で飛ぶようにした、ということです。凪のあすからの言葉を使えば、ひょっとするとひーくんは他の女の子(や男の子)にも興味を持つかもしれないよねというモデルです。

この修正したモデルでは、あるページから他のページへ少ない確率ながら必ずジャンプできるので、Rank sinkのような問題は起こりません。

早速、計算してみましょう。

octave:100> (M2^10000)*P
   0.095976 ひーくん
   0.060684 ちさき
   0.106569 まなか
   0.043208 かなめ
   0.023111 さゆ
   0.023111 美海
   0.115780 つむぐ
   0.135981 あかり
   0.259599 至
   0.135981 みをり

やはり両思いは強いのですね!(?) 子供組を見ると、やはり特に重要そうなキャラクターであるひーくん、まなか、つむぐの三人のページランクが高いので、おおよそページランクの妥当性を示せたと言っていいんでしょうか!?とりあえず、この三人に注目して今後は視聴を続けたいと思います!

蛇足:固有値・固有ベクトルを使って求める

なんか10000回も行列を掛け算するのは…という感じしますか?元の論文だとその方法で計算してると書いてあるのですが、これは行列の固有値・固有ベクトルを使っても求めることが出来ます。

固有値・固有ベクトルをおさらいしておくと、

  • ある行列Mに対する、
  • 固有値μ・固有ベクトルVは、
  • 次のような関係を満たすものです
  • MV = μV

ここで、固有値μ=1なら、MV=Vとなって、N回目とN+1回目でのそれぞれのページに居る確率が全く同じという事になり、ちょうど私達が求めたかったベクトルと一致しますね!

octaveで求めるとこんな感じになります。

octave:142> eig(M2)
ans =
   1.00000 + 0.00000i ←求めたい、固有値1の固有ベクトルは1番目だと分かる
   0.78950 + 0.00000i
   0.19713 + 0.55213i
   0.19713 - 0.55213i
  -0.86957 + 0.00000i
  -0.59237 + 0.00000i
  -0.25222 + 0.45314i
  -0.25222 - 0.45314i
   0.00000 + 0.00000i
   0.00000 + 0.00000i

octave:143> [V, D] = eig(M2); ←固有値と固有ベクトルを再度格納
octave:144> V(:,1) ←固有ベクトル1番目を取得する
ans =
   0.252078
   0.159383
   0.279899
   0.113484
   0.060701
   0.060701
   0.304091
   0.357146
   0.681825
   0.357146
↑固有ベクトルは整数倍しても固有ベクトルで、octaveはベクトルの大きさを1に正規化するので、さっきの結果と合わない
octave:146> V(:,1) ./ sum(V(:,1)) ←ベクトルの要素の合計が1となるようにすると…
ans =
   0.095976
   0.060684
   0.106569
   0.043208
   0.023111
   0.023111
   0.115780
   0.135981
   0.259599
   0.135981
↑先ほどの結果と一致する!

こんな感じで計算できます。固有値ってこんなのにも使えるんですね〜。ちなみに、この行列の固有値のうち、絶対値が一番大きいものはかならず1になることがペロン=フロベニウスの定理から言えるらしいのですが、よくわからないので解説できません…。

感想

凪のあすからは線形代数と確率過程とグラフ理論を学べる神アニメ

ていうか実は前期のアニメで一番好きだったのは「いなり、こんこん、恋いろは」で、京都と出雲にまで行ってきたのですが、それを書く日は来るのだろうか…!?

  1. 他のページへの遷移確率の和は1です。リンク総数*(1/リンク総数)=1。 []
  2. つむぐから他のページに遷移できないので、確率がそこで皆0になってしまう []

二重振り子のシミュレーション

C++から動画を直接出力できる自作ライブラリ「Dolly」を用いて、二重振り子のシミュレーション動画を作ってみました。

10分耐久・たのしい二重振り子

ソースコードはこちらにあります:

二重振り子はどこにいる?

さらに、振り子の通ったところをずっと重ねるように描いた動画も作りました。

すごいランダムな動きに見えますが、振り子のある場所はやっぱり下側が多い?

自作ライブラリでマンデルブローを描画した

C++から動画を直接出力する自作ライブラリDollyを使って、マンデルブローを拡大する動画を作りました。

この部分のソースはこちらにあります。

紙ジェクションマッピングしてみた

あけましておめでとうございます!今年もよろしくおねがいします。去年はあまりいろいろ出来なかったので、今年は少しずつでも、何かできればなと思っています。せっかく転学部までしたので物理をちゃんと身につけたい1ですし、大学再受験も本腰いれてやってきたいです。

さて、新年早速何かアウトプットしよう、ということで、去年作って放置していた「紙ジェクションマッピング」というネタを投稿します。

紙ジェクションマッピングしてみた

プロジェクションマッピングってご存知でしょうか。

20140102_07
東京駅が動き出す!? この冬見られる話題のプロジェクションマッピングまとめ

こういう感じで立体物に映像を投影するアートで、最近いろいろなところで行われているようです。

これっぽい事がしたいのですが、うちにはプロジェクタのような高いものはありません。ので、代わりに1枚50円で印刷できるコンビニのプリンタで何か似たようなことが出来ないかと模索した結果、こんな感じになりました。

仕組みは非常に単純で、ある角度から眺めると立体的に見えるように紙に印刷してあります。

20140102_01作り方ですが、まずこのマーカー画像をセブンネットプリントで印刷して(宣伝)、

20140102_03こんな感じで撮影します。このとき、カメラのファインダから眺めた時に台形が真四角に見えるようにするのがポイントです。

20140102_02この状態で撮った画像がこちら

20140102_04これを作ったソフトで変形させて適当にホワイトバランスを調整。

20140102_05これを紙に印刷して、マーカー画像を印刷した紙の代わりにこの紙をおきます。

20140102_01この状態で動画を撮ったのが、先ほどのニコ動の動画です。

20140102_06実態としては、台形に変形させただけの簡単画像加工をしただけということになりますw

今後の課題

今回は一番簡単な平面でしたが、もっといろいろなものに紙ジェクションしたいですね。たとえば、ティッシュの箱に貼り付けて、ある角度から見ると立体的に見える…とか。今位置合わせが非常にしんどいのですが、ARマーカーでなんとかできるといいんじゃないかなーとかも考えてます。3D力が問われる!!

ただ、結局プリンタの実力があまり高くないせいで、紙に写した方がそんなにリアルに見えないので限界かなという気がしています。っていうかそれがお蔵入りしてた原因です…(◞‸◟ )

  1. 物理とか数学って何やるにしても役に立つと思います []

飼い猫を殺したくて獣医師になろうとした話

昨日は一人でディズニーランドに行って客の観察をしていました。という訳で、昨日書くはずだった「つらぽよアドベントカレンダー」を今日書きます。すいません。身の上話です。獣医学科行ったのに物理系に転学部して再度医学部(人間)を目指す話です。

皆が皆ペットが好きなわけじゃない

私が育った家では猫を飼っていたのですが、私は猫が大嫌いでした。

部屋はすぐ毛だらけになってパソコンのファンを破壊してくれます。

糞尿をところ構わずする1ので掃除に気が滅入るのはもちろん、糞をした先がお気入りの本の上で再度絶望。

トイレの置いてあるケージが食事する場所のすぐ隣にある(そこぐらいしか場所が無かった)ので、ずっと糞の匂いが充満する中で家族と楽しいお食事タイム。

去勢もしていない2ので頻繁に発情してしまい、そのときは金切り声に毎日毎日耐えないといけません。メス猫なのですが、私に擦り寄ってくるのも、私の雄性を自覚させるので嫌でした。

家族みたいに「にゃんこーww」とか言う人なら良いのでしょうが、私はそうではありません。私に会った人は私の頭皮が炎症を起こしてるのを知ってると思いますが、あれも猫を飼い始めて一ヶ月後くらいに始まりました。「尋常性乾癬」といって原因不明の皮膚病なのですが、私は原因はどうも猫なんじゃないかと、どうしたって考えます。猫を飼うまでは全く問題なかったので。

そんな嫌なことがあっても猫が好きな人は「そんなの猫なんだから仕方ないじゃん」「たかが猫で」という感じなので、友人知人ならともかく、家族、しかも私一人となると我慢するしかありません。

どんどんどんどん猫、いや、ペット全般に対する嫌悪感が醸成されていきました。13年くらいですかね、人生の半分以上です。

保健所で獣医師になったら猫を駆除できる

そのストレスが最大になってたのが去年の7月ごろで、ちょうど未踏がはじまって自宅でずっと開発していた時期でした。慣れない自宅での開発の上、さらに猫でのトラブルが続いていい感じに気が滅入る毎日。発散できないストレスは他の形で発散するしかないというのが教科書にも載ってる真理で、人生がそのいけにえになりました。

去年の7月は進学振り分け3の時期でもありました。実はさらにその前の年にも進学振り分けをしていて、その時には天文学科に内定してたのですが、必修単位の問題で進級できず内定取り消しになってしまい、仕方なしに二度目の進学振り分けでした。前の年は天文だったので、再度天文か、あるいは近い分野である地球惑星物理を選ぶつもりでした。受験した時も、物理と地学で入ったので…。

が、その時本当に猫の事で頭が一杯で、獣医学科に行って保健所で猫を駆除する仕事につくことで代償的に「飼い猫を殺そう」と思ってしまいました。実際、進学振り分けのフォームに獣医学科と書いたら大分精神的に落ち着いたのを覚えています。ずいぶんいい進路にその時は見えたのです。今猫に苦しまされてても6年経ったら保健所で思いっきり猫を駆除できるんだ、と。大分精神的には楽になりながら開発できました。それは、事実です。

というわけで、私は今度は獣医学科に行くことになりました。

生物はおもしろい、が、そんな人生でいいのか

獣医学科の講義には一応全部出たのですが、想像よりはおもしろかった。生物の身体の仕組みとか全然知りませんでしたし。実習も、解剖とか生体をいじるタイプなら結構楽しく、自分がおもってたよりは自分は手先が器用なんだという事もわかりました。生物は初心者だったのですが、それでも単位は一つ落とすだけで済みました。

しかし、卒業後の進路にはどれも魅力を感じられません。元々は保健所で猫を駆除しようかと思ってましたが、人生を盛大に棒に振るような進路でそれはどうなんだろうと冷静にもなってきました。っていうか、そもそも保健所は猫は駆除しません。捨てられたペットなら処理しますが、それも最近はできるだけしないように引き止めるらしいです4。製薬会社だと動物実験で薬の効果を確かめるような仕事があったりするのですが、毎日動物を扱うのかとと思うと…。いわゆる「町の獣医さん」も動物をずっと触ってると考えると気が狂いそうです。動物嫌いには辛い進路が多い。動物を活かす仕事よりは殺す仕事のほうが多いので、動物好きでも辛そうですが…。

そして医学部再受験へ…

そういうわけで、獣医学科から、当初行く予定だったうちの片方である地球惑星物理学科に転学部しました。本当は天文学科が一番良かったのですが、天文学科に行きたい人が多くて、競争が激しくて転学部させて貰えませんでした。。。すでに一留している上で転学部のペナルティでもう一年留年になるので、二留になりました。浪人もしてるので、ついに三年遅れになってしまいました。

この騒動でただでさえ悪かった家族関係も最悪になった(もともと私が我慢してたのだから、露呈しただけとも言える)事もあり、鬱病になってしまいました。獣医学科に居た時に一番おもしろかった話は脳に効く薬の話5だったのですが、期せずして自分の身体で実験することに…。実際飲んでみると、想像以上に精神状態が変わるので、結構驚けます。人間の脳はこんなにも簡単に操作できてしまうのか、と。

物理系は出ても就活がただでさえつらぽよな上3年遅れとなり、このままこのコースを行くのもつらぽよな気がしています。就活するならSEとかプログラマ、あるいはいっそ就活なんかせずにフリーランスや起業ですが、そもそも私はプログラミングを生涯の仕事にしたら幸せにはなれないタイプの人間な気がする(またその話はいたしましょう)ので、卒業したら医学部に入り直そうと今では考えています。ちなみに、獣医学科でも卒業後に人間の医学部に入り直す人は珍しくありません。

医学は確かに興味深く学べたので、動物の医学にしたのは失敗でしたが、人間なら進路も私にも魅力的な進路が多いので、なれると思ってます。医師という仕事は、他の業種みたいに「需要を作りだ」したりせずに、病気の治療や健康維持という確実にある需要に確実に答える仕事である、というのが一番魅力に感じている点です。

私の作っているソフトウェアも、実はそうなのです。さきゅばすを作ったのは、ネットが無くてニコ動を見られない友達が家や学校でもコメント付きの動画が見られるようにですし、ど〜なっつも、時間が戻れるようなプログラミング言語が無いとコンテンツを作るのが明らかにしんどいからです。単純な技術的興味から書くことももちろんたくさんありますが、基本的には困ったことを解決するためにプログラムを書いています。

もし仮に受かったらなのですが、ある程度人生経験は積んだつもりですし、以前から興味もあったので、精神科の医師になりたいです。メンタルヘルスの重要性は鬱病になったこともあってよくわかりました。。。他には産業医も良いですね。

まあでも、受かるかどうかはわからないよね。駄目だった時用に情報系大学院の準備などもしておくつもりです。

蛇足:動物が嫌いな人はペットをどう思っているか

この際なのでちょっと書いてみます。

私は動物が嫌いで、実習でラットやマウスを殺すことにも特に良心の呵責とかを感じない人間ですが、それでもペットという営みは嫌いです。生命をおもちゃにして人形ごっこしてるように思います。人間があの「かわいい姿(人口的な作り物っぽさと獣の嫌な部分の足し算の、どこが可愛いのか全然分からないけど)」を生み出すために無理して交配(端的に言えば近親相姦)した結果出てきた遺伝病で苦しむペットが多いのも、命を弄んでるように思います。その自覚が無いまま、「かわいいかわいい」で飼うのには私は嫌悪感を感じます。

ちなみに最近上野動物園に行ったら、小動物が軒並み常同行動(同じ場所に居続けるストレスによって同じ行動を只管しつづけるようになること)してたので、動物園という存在もどうなのかなぁと思うようになりました。

  1. いつもではない、が、二週間に一度だって嫌だ []
  2. 交尾させて儲けたかったそうだが、その時点で無知丸出しですな []
  3. 東大は入ってから学科が決まるのですが、その制度をこう呼びます []
  4. 野山に放たれて余計面倒な問題になるだけじゃないかなぁ []
  5. 動物でも問題行動を起こしてると処方されることがある []

クリスマスはNoiz2saをAndroidで!

おひさしぶりです。21日のニコニコ学会βは良かったですね。特に良かったのは研究100連発で、本職の人のすごさを見ると、私も何かつくりたいな〜って気分になります。私はボランティア(登壇者へのカンペ出し)として参加しました。凄い人を間近に見れるという特典付き(?)の嬉しいけど責任重大なお仕事です。イベント運営って大変ですね…。あと野糞の人は本気でマッドだ(褒め言葉)。

Noiz2saをAndroidに移植してみた!

Noiz2saはABA Gamesさんの開発した、アブストラクトな感じのシューティングゲームです。私はこのゲームがSTGに入り始めたきっかけだけあって、大好きです。

去年はWebに移植したNoiz2saを公開しましたが、今年はAndroidに移植したものを公開します。

クリスマス、そして年末年始といえばシューティングゲーム!ぜひお楽しみ下さい。敵の弾を避けて撃つ、指一本で遊べるシンプルなゲームです。

USBのゲームパッドをつなぐと、PC版と同様な操作で遊ぶことも可能です。

ダウンロードはこちらから

バッドノウハウFAQ

移植時のバッドノウハウを書きます

Q.リソースのファイルの拡張子がamrとかなのはなぜですか

このように置かないと、リソースが圧縮される時に圧縮されてしまい、AssetManagerからランダムリードすることが出来なくなってしまうからです。

あっコーナー分けた割に一個しか思い浮かばなかった…質問は随時受け付けています。