変電録

はてなブログを開設しました

招待制じゃ無くなったので遊び場にしてみる

正月なのでWebGLでおっぱい曲面を描く - 辺々自乗

Rubyで書いたTwitterのbotをOAuthに対応させたメモ

Twitter, ruby |

概要

Twitterが今月でBASIC認証を受け付けなくなると聞いたので、手製の単機能Twitter botである@ISSaboveTokyoをOAuthに対応させる改修を行った。そのときにやったことのメモ。

試した環境

• Ubuntu Linux 8.04(x86)

インストールするgem

sudo gem install twitter

twitter(http://twitter.rubyforge.org/)をインストールした。

Ubuntu 8.04だとgemが古くなっていたので

sudo gem update --system

などしてアップデートした。その仮定で色々と問題があったけど、適当に解決してしまってよく覚えていないので省略。

アプリケーション(bot)の登録、consumer keyの生成

自分のアカウント（botのではなく）*1でTwitterにログインした状態で

http://twitter.com/apps

にアクセスしてbotを登録する。

アイコンや説明などは適当に設定する。Application TypeはClient、Default Access TypeはRead&Writeにする。

設定が終わると、consumer keyとconsumer secretという文字列が表示される。これらは次のrequest tokenの生成で利用する。

request tokenの生成

PINの使い方がよくわからなくて一番苦労したところ。

no titleを参考にした。

require 'rubygems'
require 'twitter'

OAUTH_CONSUMER_KEY = #先ほどのconsumer keyを文字列として代入
OAUTH_CONSUMER_SECRET = #先ほどのconsumer secretを文字列として代入

oauth = Twitter::OAuth.new(OAUTH_CONSUMER_KEY, OAUTH_CONSUMER_SECRET)

request_token = oauth.request_token

puts "request token: #{request_token.token}"
puts "request secret: #{request_token.secret}"
puts "authorize url: #{request_token.authorize_url}"

これを実行して表示された authorize url を、今度はbotのアカウントでTwitterにログインした状態で開く。開いたページでボタンをクリックしてアクセスを許可すると、PINと呼ばれる番号が表示される。この番号と、request token/secretは次のaccess tokenの生成で利用する

access tokenの生成

require 'rubygems'
require 'twitter'

OAUTH_CONSUMER_KEY = #先ほどのconsumer keyを文字列として代入
OAUTH_CONSUMER_SECRET = #先ほどのconsumer secretを文字列として代入

puts "request token?"
OAUTH_REQUEST_TOKEN = gets.chomp!
puts "request secret?"
OAUTH_REQUEST_SECRET = gets.chomp!
puts "PIN?"
OAUTH_REQUEST_PIN = gets.to_i

oauth = Twitter::OAuth.new(OAUTH_CONSUMER_KEY, OAUTH_CONSUMER_SECRET)

oauth.authorize_from_request(OAUTH_REQUEST_TOKEN, OAUTH_REQUEST_SECRET, OAUTH_REQUEST_PIN)

at = oauth.access_token

puts "access token: #{at.token}"
puts "access secret: #{at.secret}"

上記のコードでaccess token/secretが取得できるので、以降これを使ってbotを動かす。

access tokenを利用してbotを動かす

タイムラインを表示するサンプル

require 'rubygems'
require 'twitter'

OAUTH_CONSUMER_KEY = #先ほどのconsumer keyを文字列として代入
OAUTH_CONSUMER_SECRET = #先ほどのconsumer secretを文字列として代入

OAUTH_ACCESS_TOKEN = #先ほどのaccess tokenを文字列として代入
OAUTH_ACCESS_SECRET = #先ほどのaccess secretを文字列として代入

oauth = Twitter::OAuth.new(OAUTH_CONSUMER_KEY, OAUTH_CONSUMER_SECRET)

oauth.authorize_from_access(OAUTH_ACCESS_TOKEN, OAUTH_ACCESS_SECRET)

client = Twitter::Base.new(oauth)
client.friends_timeline.each{|tw|
  puts "#{tw.name}: #{tw.text}"
}

uBLASのExpression Templateまわりを探検したメモ

C++ |

uBLASのExpression Templateについて、どういう実装になっているのか少し覗いてみてわかったことのメモ

Expression Templateとは

Expression Template(ET)とは、ユーザ定義型に対して"C = A + B"といった式があったときに、A + B の結果が入った一時オブジェクトを作ってから内容をCにコピーするのではなく、A + B という式を表現するオブジェクトを作って計算はCに代入する直前に行うことで、無駄なコピーや一時オブジェクトの生成を無くしてしまうことを狙った手法のこと。

詳しくは

Expression Template - Faith and Brave - C++で遊ぼう

などのちゃんとした解説を参照。

見てみるコード

以下の、行列の足し算とその結果の代入というコードの評価について追いかけてみた

boost::numeric::ublas::matrix<double> A, B, C;

/* initialize A and B*/
...

C = A + B;

boostのバージョンは1.34.1(Ubuntu 9.04のapt-getで取ってきたのをそのまま使ってるから)

行列の生成と初期化

matrix<double>という型について、fwd.hppに書かれている matrix の先行宣言を見てみると

    template<class T, class L = row_major, class A = unbounded_array<T> >
    class matrix;

という風にTの後ろにテンプレート引数が2つ続いているのがデフォルト引数で隠されていたことがわかる。

L というのはデータの格納順を変更するためのパラメータで、row_majorのとき行優先、column_majorのとき列優先。

その後の A というのはデータを格納する場所を指定するためのパラメータで、unbounded_arrayのときはヒープに、bounded_arrayのときはスタックにメモリを確保する。

こうして宣言された行列のAとBについて何か値をセットしたとして、続いて本題であるところの C = A + B という文の評価が行われることになる。

式を表現するオブジェクトの生成(operator+)

まず matrix_expression.hpp に書かれている operator+ によって A + B の評価が行われる

    template<class E1, class E2>
    BOOST_UBLAS_INLINE
    typename matrix_binary_traits<E1, E2, scalar_plus<typename E1::value_type,
                                                      typename E2::value_type> >::result_type
    operator + (const matrix_expression<E1> &e1,
                const matrix_expression<E2> &e2) {
        typedef typename matrix_binary_traits<E1, E2, scalar_plus<typename E1::value_type,
                                                                              typename E2::value_type> >::expression_type expression_type;
        return expression_type (e1 (), e2 ());
    }

matrix<T> は matrix_expression<matrix<T> > をpublic継承しているので*1、このoperator+の引数にマッチする。

ここで戻り値は expression_type(e1(), e2()) となっているが、このexpression_typeがどういう型になっているのか、matrix_expression.hpp に書かれている matrix_binary_traits の定義を見てみると

    template<class E1, class E2, class F>
    struct matrix_binary_traits {
        typedef matrix_binary<E1, E2, F> expression_type;
#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
        typedef expression_type result_type; 
#else
        typedef typename E1::matrix_temporary_type result_type;
#endif
    };

となっているので、戻り値は matrix_binary という型だとわかる。*2

そして、そのコンストラクタに渡されている e1()とe2()は、matrix_expression<E> の operator() が

        BOOST_UBLAS_INLINE
        const expression_type &operator () () const {
            return *static_cast<const expression_type *> (this);
        }

という風になっているので、matrix_binaryのコンストラクタにはmatrix_expressionからダウンキャストされた元のmatrix<double>のconst参照が渡されることがわかる。

ではその matrix_binary はどこにあるのかというと、再び matrix_expression.hpp が登場。これは長いので一部だけ抜粋

    template<class E1, class E2, class F>
    class matrix_binary:
        public matrix_expression<matrix_binary<E1, E2, F> > {
    public:
        // Construction and destruction
        BOOST_UBLAS_INLINE
        matrix_binary (const E1 &e1, const E2 &e2): 
            e1_ (e1), e2_ (e2) {}

        // Element access
        BOOST_UBLAS_INLINE
        const_reference operator () (size_type i, size_type j) const {
            return functor_type::apply (e1_ (i, j), e2_ (i, j));
        }
    private:
        expression1_closure_type e1_;
        expression2_closure_type e2_;
    }

つまり A + B の戻り値として返る matrix_binary_traits<E1, E2, scalar_plus<E1::value_type, E2::value_type> >::result_type という型は、AとBの要素の値でscalar_plusを呼び出した戻り値を自分の要素の値として返す。

例えば、i, j が有効なインデックスであれば (A + B)(i, j) と A(i, j) + B(i, j) は同じ値になる。*3

というわけで、 A + B という式が評価された結果、 matrix_binary<matrix<double>, matrix<double>, scalar_plus<double, double> >(A, B) という2つの行列の足し算を表現する一時オブジェクトが生成されることがわかる。

計算の実行と結果の代入(operator=)

C = A + B

という式について、 A + B がどう評価されるかわかったので、続いて C への代入について見てみる。

ここで呼び出される operator=*4 はどうなっているのか、matrix.hppを見てみるとuBLASのmatrixにはメンバとして3つのoperator=がある。

        // Assignment
        BOOST_UBLAS_INLINE
        matrix &operator = (const matrix &m) {
            size1_ = m.size1_;
            size2_ = m.size2_;
            data () = m.data ();
            return *this;
        }
        template<class C>          // Container assignment without temporary
        BOOST_UBLAS_INLINE
        matrix &operator = (const matrix_container<C> &m) {
            resize (m ().size1 (), m ().size2 (), false);
            assign (m);
            return *this;
        }
        template<class AE>
        BOOST_UBLAS_INLINE
        matrix &operator = (const matrix_expression<AE> &ae) {
            self_type temporary (ae);
            return assign_temporary (temporary);
        }

1つ目は同じ型を引数としてとる普通の代入、2つ目は型は違うけど自分でデータを持っている（コンテナである）行列を引数としてとる代入、そして3つ目は自分ではデータを持たず式を表現しているだけの（コンテナでない）行列を引数としてとる代入。

ここで引数として渡されているのは自分でデータを持たず2つの行列同士の足し算を表現しているだけのmatrix_binaryであるので、呼ばれるのは3つ目の"operator = (const matrix_expression<AE> &ae)"になる。

さて、このoperator=で何が行われているか見てみると

1. 自分と同じ型のテンポラリオブジェクトを作り、渡された引数である"ae"で初期化（コンストラクタの呼び出し）
2. テンポラリオブジェクトを使って自分の内容を書き換え (assign_temporary)

と2つに分かれることがわかるが、2つ目のassign_temporaryは

        BOOST_UBLAS_INLINE
        matrix &assign_temporary (matrix &m) {
            swap (m);
            return *this;
        }

と、事実上ただのswapの呼び出し*5と同じなので、重要な処理はほとんどコンストラクタの呼び出しの中で行われていることがわかる。*6

ここで呼び出されるコンストラクタは以下のようになっているので、

        template<class AE>
        BOOST_UBLAS_INLINE
        matrix (const matrix_expression<AE> &ae):
            matrix_container<self_type> (),
            size1_ (ae ().size1 ()), size2_ (ae ().size2 ()), data_ (layout_type::storage_size (size1_, size2_)) {
            matrix_assign<scalar_assign> (*this, ae);
        }

matrix_assign<scalar_assign>で実際の代入が行われるとわかる。

そこで、detail/matrix_assign.hppを見てみる。

    // Dispatcher
    template<template <class T1, class T2> class F, class M, class E>
    BOOST_UBLAS_INLINE
    void matrix_assign (M &m, const matrix_expression<E> &e) {
        typedef typename matrix_assign_traits<typename M::storage_category,
                                              F<typename M::reference, typename E::value_type>::computed,
                                              typename E::const_iterator1::iterator_category,
                                              typename E::const_iterator2::iterator_category>::storage_category storage_category;
        // give preference to matrix M's orientation if known
        typedef typename boost::mpl::if_<boost::is_same<typename M::orientation_category, unknown_orientation_tag>,
                                          typename E::orientation_category ,
                                          typename M::orientation_category >::type orientation_category;
        typedef basic_full<typename M::size_type> unrestricted;
        matrix_assign<F, unrestricted> (m, e, storage_category (), orientation_category ());
    }

storage_categoryは密行列同士の足し算であるからdense_proxy_tagとなるので、呼び出されるのは matrix_assign (M &m, const matrix_expression<E> &e, dense_proxy_tag, C) 。（Cは行優先か列優先かを表す型）

    // Dense (proxy) case
    template<template <class T1, class T2> class F, class R, class M, class E, class C>
    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    void matrix_assign (M &m, const matrix_expression<E> &e, dense_proxy_tag, C) {
        // R unnecessary, make_conformant not required
        typedef C orientation_category;
#ifdef BOOST_UBLAS_USE_INDEXING
        indexing_matrix_assign<F> (m, e, orientation_category ());
#elif BOOST_UBLAS_USE_ITERATING
        iterating_matrix_assign<F> (m, e, orientation_category ());
#else
        typedef typename M::difference_type difference_type;
        size_type size1 (BOOST_UBLAS_SAME (m.size1 (), e ().size1 ()));
        size_type size2 (BOOST_UBLAS_SAME (m.size2 (), e ().size2 ()));
        if (size1 >= BOOST_UBLAS_ITERATOR_THRESHOLD &&
            size2 >= BOOST_UBLAS_ITERATOR_THRESHOLD)
            iterating_matrix_assign<F> (m, e, orientation_category ());
        else
            indexing_matrix_assign<F> (m, e, orientation_category ());
#endif
    }

BOOST_UBLAS_USE_INDEXING や BOOST_UBLAS_USE_ITERATING や BOOST_UBLAS_ITERATOR_THRESHOLD などのマクロによって、代入を行うループの実装が切り替えられるようになっていることがわかる。これらのマクロは detail/config.hpp で定義されていて、見てみると BOOST_UBLAS_USE_INDEXING が #define されているので添字によるループで代入が行われるとわかる。

indexing_matrix_assignのrow_major（行優先）な実装は次の通り。

    // Explicitly indexing row major
    template<template <class T1, class T2> class F, class M, class E>
    // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    void indexing_matrix_assign (M &m, const matrix_expression<E> &e, row_major_tag) {
        typedef F<typename M::reference, typename E::value_type> functor_type;
        typedef typename M::size_type size_type;
        size_type size1 (BOOST_UBLAS_SAME (m.size1 (), e ().size1 ()));
        size_type size2 (BOOST_UBLAS_SAME (m.size2 (), e ().size2 ()));
        for (size_type i = 0; i < size1; ++ i) {
#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
            for (size_type j = 0; j < size2; ++ j)
                functor_type::apply (m (i, j), e () (i, j));
#else
            size_type j (0);
            DD (size2, 2, r, (functor_type::apply (m (i, j), e () (i, j)), ++ j));
#endif
        }
    }

BOOST_UBLAS_USE_DUFF_DEVICEを#defineすることで、forループの代わりにDuff’s deviceを利用することもできるとわかる。しかし、手元の環境で"#define BOOST_UBLAS_USE_DUFF_DEVICE"して試してみたらかえって遅くなった。

わかったこと

• uBLASのmatrix は、テンプレート引数を使って要素の型以外にスタック／ヒープの切り替えや行優先／列優先の切り替えも指定できる
• uBLASの演算子オーバーロードは、CRTPを利用している
• uBLASでは行列の要素同士の2項演算は matrix_binary という型で表現される
• uBLASは detail/config.hpp を書き換えることで行列の代入に関する動作を切り替えることなどができる

Effective C++を読んだらExceptional C++の良さがわかった

本, C++ |

数ある関連書籍の中からC++について理解を深めるのに最適な一冊は何かという議論を追っていると、「Exceptional C++だけ読んでいれば十分」という意見を目にすることがあります。その後には「Effective C++はExceptional C++が書いているのと同じことを冗長にして書き直しているだけ」という言葉が続いていることも。

たしかにExceptional C++は幅広い内容をしっかりと説明してくれます。例外安全なコードを書くにはどうすればいいか、無駄な再コンパイルを防ぐにはどうすればいいか、C++において「インターフェース」とは何を指していて何から構成されているのか、などなど。しかし、「なぜコードをそのように設計したいのか」「なぜそのような事柄が問題になるのか」といったことについてExceptional C++は説明してくれません。Effective C++が「例外安全」という言葉を持ち出してから数ページを割いて例外安全や「強い保証」といった言葉の定義を明確にしているのに対して、Exceptional C++は「例外安全（例外が発生しても適切に処理する）」という簡単な注釈をカッコ内に付けただけで先に進んでしまいます。

確かにいい本なんだけど、世間で大絶賛されてるほどいい本だとは思えなかった。読む側として経験値が足りないってことなのかも。微妙。

Exceptional C++ - 変電録

以前、僕はこのように書きました。今ふり返って、この「経験値が足りない」という認識は正しかったと思います。以前より少しだけ経験値が上がった現在の僕は、Exceptional C++の構成は数多く存在するC++関連書籍の中でもツボをかなり上手く押さえたものになっていると、なんとなくですが思えるようになりました。

ですが、その経験値をあげるのに役立った本を一冊あげるとすれば、それはEffective C++でした。Effective C++はExceptional C++よりも冗長かもしれませんが親切です。読者が不慣れだと思われる単語には必ず定義を与えてから先に進みますし、重要な概念についてはそれをどのように実現するのかだけでなく、なぜそれが必要なのかまで噛んで含めるように解説を与えてくれます。また、印刷も2色刷りで見栄えがよく読みやすく、第3版ではTR1についても触れています。一方でExceptional C++は、例えばスマートポインタについてはauto_ptrにしか触れていません。これではshared_ptrの存在を知っている読者は、本全体に書かれている知識のうち一体どれだけが現在でも有効なものなのか判断がつかず、読む気を失ってしまいます。

思うに、対象読者を広げる努力あるいは対象としない読者を明示する努力が多くのプログラミング関連書籍には欠けているのではないでしょうか。例えば一般的な理工学書であれば、前書きで「この本は読者が高校までで学習する微積分を一通り身につけていることを想定している」などと前提とする知識を明示し、最初の1章はその本全体の目的の解説あるいは用語や記号の定義にあてるなどするのが一般的です。翻ってExceptional C++を見ると、そこには

本書は、C++の基本をすでに心得ている読者を対象にしている。まだまだ、という方は、まず、C++の学習を進め、大枠をつかんでいただきたい

とだけ書かれています。「C++の基本」とは、「大枠」とはなんでしょうか。テンプレートを使ったメタプログラミングは「基本」に入るのでしょうか、式の内部に含まれる下位式の評価順が定義されていないことについてはどうでしょうか。

プログラミング教育について標準的な課程といったものがなく、「大学○年までで学ぶ内容」といった表現が使えないという問題があるのはわかります。しかし、だからといってこれだけの表記では想定していないレベルの読者を追い返すには不十分ですし、これ一冊で大丈夫だからと言い含められた哀れな初心者に至っては尚更です。

ある本を賞賛する読者の後ろに、その本を読みきれなかった読者、最後まで目を通しはしたけれど評価できるほどの知識を持ち合わせていなかった読者が大勢いることは忘れられがちです。ある本について良い感想を持った人ほどその本について語りたがるというのが実に自然な心理だというのも事実でしょう。Exceptional C++を賞賛するのは正しいことです。ですが、Exceptional C++と内容が重複していること、表記が平易で冗長であることを理由にEffective C++を否定することは正しいことでしょうか。挫折した読者がより平易な本によってすくい上げられる可能性を潰すこのような言葉は単純に悪なのではないでしょうか。本が不親切であるときに、紹介者が若干の親切さを働かせて初学者を助けられるようにすることがそこまで困難な仕事だと、僕は思いません。

ブラウザで動く2次元粒子法非圧縮性流体シミュレータをJavaScriptで書いてみた

JavaScript |

作ったもの

http://www.henshi.net/archive/mps/

↑のURLをブラウザで開いて、「開始／停止」にチェックを入れると動く

• 四角い容器の中に現れた水柱が崩れていくような状況を再現している
• Firefox 3.0.8, Google Chrome 1.0.154.53で動作確認*1(IEには対応してない)
• main.js中のprepareParticles関数で粒子の初期配置を行っているので、いじれば違った状況も再現できる
• ソースコードはこちら(http://www.henshi.net/archive/mps/mps.zip)

粒子法とは

流体の挙動を計算機上でシミュレーションするには、大きく分けて2つの方法がある。

一つは格子法で、これは空間を格子で分割してしまい各格子について質量や運動量などの出入りを計算する方法。

そしてもう一つは粒子法で、これは流体をいくつかの塊(=粒子)に分割してしまい*2、粒子を移動させながら粒子間の相互作用から圧力や速度を計算する方法。

大雑把な比較を行うと、格子法では空間微分の計算が簡単に行えるが、流れによって運ばれる質量や運動量の計算に工夫が必要。更に、重要なところでは格子を細かくして重要でないところでは負荷を減らすために格子を粗くするというようなことは格子をうまく切らないと実現できないので、その手間もかかる。

一方、粒子法では流れによって運ばれる質量や運動量は粒子が移動することで勝手に表現されてしまうので簡単に計算できるし、そもそも格子を使わないので格子の切り方に工夫をする必要もない。また、複数の流体が混じっていたり、変形する物体が流れの中に存在するような状況の計算も比較的簡単。ただし、粒子間の相互作用や空間微分の計算は格子法よりも難しく、計算コストもかかる。

非圧縮性流体とは

水などの密度変化が無視できる流体を非圧縮性流体と呼ぶ。空気も時速100km/h未満の流れなどでは密度変化がほぼ無視できるので、非圧縮性流体として取り扱うことができる。

一般に、圧縮性流体よりも非圧縮性流体の方が計算が簡単になる。

ただし、亜音速以上の流れや衝撃波などの現象は圧縮性流体として取り扱わないと計算できない。

使用した手法

実装するにあたっては、「数値流体力学」中で示されているMPS法(moving particle semi-implicit method)を採用した。

MPS法では、各粒子の位置で流速の微分などを計算するとき、粒子間の距離rに対して

で表される重み関数w(r)を使って計算を行う。

そして、粘性力や外力による粒子の移動は陽解法により計算し、圧力の計算については陰解法により計算を行う。

ただ、この重み関数w(r)をこのまま使用すると粒子同士が接近したときにrの値が小さくなってw(r)が極端に大きくなってしまうので、w(r)に上限値を定めて、rが一定以上小さくなったときはw(r)が定数になるよう変更を行った。*3

作ってみた感想

• Chromeすごくはやい
• 始めはSylvesterという行列・ベクトルのライブラリを使ってたのだけど、一時オブジェクトが大量に出来てFirefoxでの動作がままならなくなったので自作のショボいクラスで代用した
• いちおうそれっぽく動く
• 全体の動きが収まってからも各粒子がウネウネ動いて気持ち悪い