2011年10月15日土曜日

MATLABでQhullを使う

Qhull (http://www.qhull.org/) は計算幾何に関する有用なライブラリーで,MATLABの幾つかの関数,delaunayn や voronoin は Qhull を利用しています.

しかし,これらの関数は Qhull が持つすべての機能を実装しているわけではなく,例えばドロネー分割における法線や隣接単体の情報を出力するオプションはMATLABでの実装ではサポートされていません.

本来なら Qhull を呼び出すmexファイルを書くのが正攻法ですが,面倒なのでコンパイル済みのQhullの実行ファイルを呼び出す関数を作ってみました.

僕が欲しかったのはドロネー単体分割時の隣接単体の情報だったので,以下の関数はドロネー単体の情報に加えて単体間の隣接行列を返す仕様になっています.
function [tes,adj]=mydelaunayn(X)
%MYDELAUNAYN N-D Delaunay triangulation with adjacency matrix
%   [tes,adj]=mydelaunayn(X)
%   tes: sparse matrix which represents tessellation
%      tes(i,j)==1 if j-th simplex has vertex X(i,:)
%   adj: adjacency matrix among Delaunay simplices
%      adj(i,j)==1 if i-th and j-th simplices has common hyper plane
%
%   NOTE:
%     This routine requires excutable of Qhull.
%     http://www.qhull.org/


%% Compose input
infile=tempname;
fid=fopen(infile,'w');
d=size(X,2);
N=size(X,1);
fprintf(fid,'%d\n',d);
fprintf(fid,'%d\n',N);
fprintf(fid,strcat('%e',repmat('\t%e',1,d-1),'\n'),X');
fclose(fid);

%% Invoke Qhull
outfile=tempname;
opts={'Qt','Qbb','Qc'};
if d>3
    opts{end+1}='Qx';
end
opts=strcat(opts,{' '});
system(sprintf('qhull d s i %s Fn < %s > %s',cat(2,opts{:}),infile,outfile));

%% Read Results
fid=fopen(outfile,'r');

% compose tessellation matrix
Ntes=fscanf(fid,'%d',1);
tesmat=fscanf(fid,'%d',[d+1,Ntes])'+1;
idx=1:Ntes;
tes=sparse(reshape(tesmat',[],1),reshape(idx(ones(d+1,1),:),[],1),1,N,Ntes);

% compose adjancy matrix
Nadj=fscanf(fid,'%d',1);
adjout=fscanf(fid,'%d',[d+2,Nadj])';
assert(all(adjout(:,1)==d+1));
assert(Nadj==Ntes);
adjmat=[reshape(idx(ones(d+1,1),:),[],1),reshape(adjout(:,2:end)',[],1)];
adjmat(adjmat(:,2)<0,:)=[];
adj=sparse(adjmat(:,1),adjmat(:,2)+1,1,Ntes,Ntes);

%% Cleanup
delete(infile);
delete(outfile);
end


上の関数の実行にはパスに Qhull の実行ファイルがあることが必要で,これは http://www.qhull.org/ からダウンロードできます.他の用途への改造も比較的容易だと思うので,よければ参考にしてください.

MATLAB の無名関数で複数の文を実行する方法

通常 MATLAB では関数を定義する際に別ファイルを作成せねばならず,これは結構面倒です.
無名関数を使うと
f = @(x) x^3-x;
ezplot(f,[-2,2])
のように簡単な関数をその場で定義できますが,複数のステートメントを含む関数を書くことは簡単ではありません.
ここでは複数の文を含む無名関数を定義するために僕が使っているテクニックを紹介します.

複数の文を含む無名関数を書く場合,「任意の数の式を引数とし,最後の式の値を返す関数 progn 」を定義しておくと便利です.
progn = @(varargin) varargin{end};

f = @() progn(...
    figure(1),    ...
    ezplot(@sin), ...
    figure(2),    ...
    ezplot(@cos) );

f();

上の例では複数の図のプロットを行う関数を定義しています.

ただし,この方法では値を返さない関数や文を並べることができません.たとえば hold on などがこれに該当します.どうしてもこれらの関数を使いたい場合,邪道ですが evalc を使うことで問題を回避できます.
progn = @(varargin) varargin{end};

f = @() progn(...
    figure(1),      ...
    ezplot(@cos),   ...
    evalc('hold on'), ...
    ezplot(@sin));

f();

また,MATLABの無名関数はローカル変数を持つことができませんが,evalc と組み合わせると引数に対応する変数をローカル変数として利用することが可能です.
progn = @(varargin) varargin{end};

% 二番目に大きい固有値を返す関数を定義
eig_2nd = @(x) feval(@(x,e,idx) ...
    progn( ...
        evalc('e=eig(x);'), ...
        evalc('[~,idx]=sort(abs(e));'), ...
        e(idx(2))), ...
    x,0,0);

A=[4,0,0;
   0,1,0;
   0,0,2];
eig(A)
eig_2nd(A)
見てのとおり邪道で,書き捨てのプログラム以外で使うことは憚られますが,これを使うとかなり広範な関数をその場で定義できます.

2011年10月5日水曜日

MATLAB の無名関数内で分岐を行う方法

MATLAB には一応無名関数がありますが,関数型のプログラミング言語ほど周辺のユーティリティが充実していないので不便なことが結構あります.最近は MATLAB での無名関数の使い方もだいぶわかってきた気がするので,ちょくちょく小技を紹介していこうかと思います.

今回は僕が使っている分岐を実現する方法を紹介します.MATLABにはC言語のような3項演算子や条件分岐関数が無いために,無名関数内では手軽に分岐を行えません.

分岐を簡単に実現するひとつの方法は,あらかじめユーティリティ関数を作っておくことです.
たとえば僕は以下のような関数 switchfunc をパスに入れています.
function [ out ] = switchfunc( varargin )
%%SWITCHFUNC 条件分岐を簡潔に書くための関数
%  SWITCHFUNC( cond_1, value_1,...
%              cond_2, value_2,...
%                     :
%               true, default_value)
%  cond_1,cond_2, の中で 最初に true の値をとるものを
%  cond_k とすると,value_kを値として返します.
%  全ての条件が満たされない場合エラーを返します.

for k=1:2:nargin-1
    if varargin{k}
        out = varargin{k+1};
        return
    end
end
error('switchfunc: no match');
end

これを使うと,たとえばxの値に応じて切り替わる関数を以下のように書けます
f = @(x) ...
    switchfunc( -2 <= x && x <= -1,  x+2, ...
                -1 <  x && x <   1,    1, ...
                 1 <= x && x <=  2,    x);

ezplot(f,[-2,2]);
結果得られる関数 f  のグラフは↓のようになり

自然な表記で切り替えを持つ関数を作れていることがわかります.

しかし,実はこの switchfunc による分岐と 通常の if 文による分岐には大きな違いがあります.
それはswitchfunc を使った場合,実際の条件分岐を行う前に全てのケースについての値が計算されているという点です.
先ほどの例題のように,値の計算に副作用がない場合には速度が多少遅くなる程度の害しかありませんが,plot のような関数を条件に応じて切り替える用途にはそのままでは使えません.

通常のif 文のように条件と合致した時だけ値の評価を行いたい場合は以下のように一工夫します.
swplot = @(num) ...
            feval( ...
                switchfunc(...
                    num==1, @() ezplot(@sin),...
                    num==2, @() ezplot(@cos)));
figure(1), swplot(1)
figure(2), swplot(2)
ポイントはswitchfuncで切り替える値を無名関数にしていることで,帰ってきた結果のみを feval で実行することによって通常のif文と同様に条件に対応する値のみが評価されます. 

MATLAB の無名関数は副作用がない純粋関数しか含まないことを前提に設計されているので,無名関数のなかで完結するようプログラムを書くと自然にLISPのような関数型言語にちかい構造のプログラムになる気がします. MATLAB は元来このようなスタイルを前提とする言語ではないので表記は複雑になりがちですが,関数型言語で使われるテクニックの中にもMATLABで有用なものがあるかもしれませんね.

Blogger に TeX の数式を入れる簡単な方法

もっとスマートな方法があったのでこちらの記事もご参照ください(2012/1/19 追記)

Bloggerでとりあえず数式を入れたい時の簡単な方法をメモしておきます.
  1. TeXclip で数式を作成
  2. 表示されている数式を右クリックして「画像のURLをコピー」でコピー(IEだと多分できない)
  3. Blogger のエディタで「画像を挿入」を選択し,「URLから」のダイアログでコピーしたURLを指定
出来上がりの図↓
行の途中に数式を入れたりはできないのですが,簡単なので結構便利かもしれません.

Blogger以外のブログシステム( Wordpress など)でも,画像のURLから画像をダウンロードして挿入する機能があれば使える気がします.