matlab: 構造体の利用例：構造体の各フィールドをワークスペース上の個別の変数に分解する

matlab の構造体と言えば，いくつかの変数を一つの塊にまとめてパッケージしたようなものだ．以前，ワークスペース上の変数を一つの構造体にまとめるという記事を書いた．今回は，逆に，一つの構造体から，個々の項目を別個の変数として，ワークスペース上に分解する方法を述べる．

構造体を分解し，各フィールドを個別の変数として取り出す．save/load や関数引数の処理などで便利なことがある

• 利用する関数：fieldnames, eval, assignin

では，行ってみよう．

matlab: ブロック対角成分の抽出

matlab で，ある行列からブロック対角成分を抽出したい場合， 効率的な実装はあるだろうか．

少し調べてみよう．

matlab: ここはどこ？ dbstack

matlab の m ファイルの実行中に keyboard 命令で止まったとき，今どの m ファイル中の何行目で止まっているかを知るには，dbstack コマンドを実行すればよいという話．

デバッグのお供 keyboard 命令はとても便利なので，なじみのない人も続きを読んでいただきたい．では，行ってみよう．

matlab: 一目で分かる多次元インデクス

matlab の配列参照に使うインデクスの基礎について．

まずはこれ．一目で分かる多次元インデクス！
a を [4 by 3 by 2] なる 3 次元配列とすると，

A three dimensional [4 by 3 by 2] array

その要素の番号は

Index / subscript correspondences of [4 by 3 by 2] array 'a'

となるのだ～！ この表の意味は，例えば「a(2,3,2)とa(2,6)とa(22)は同じ要素を指す！」ということ．行列(2次元配列)の場合は，真ん中と右に注目してもらえばOK. さっきの例では a(2,6) と a(22) だね．

もうわかった？ともかく説明しよう．

3 次元配列 a の参照の仕方は

1. a(i1, i2, i3)
2. a(i1, i2)
3. a(i1)

の3通り．s = size(a) とすると，s = [4,3,2] であって，式で書けば

通し番号 = (i3-1) * s(2) * s(1) + (i2-1) * s(1) + i1

となる．方式 2 では i3 = 1, 方式 3 では i3 = 1, i2 = 1 を補えばよい．例えば，

(3, 1, 2) の要素の通し番号 = (2-1) * 3 * 4 + (1-1) * 4 + 3 = 15
(3, 4) の要素の通し番号 = (1-1) * 3 * 4 + (4-1) * 4 + 3 = 15
(15) の要素の通し番号 = (1-1) * 3 * 4 + (1-1) * 4 + 15 = 15

となるわけだ．上の図と見比べてみると，方式 1, 2, 3 と通し番号との対応がわかるだろう．

多次元インデクス(サブスクリプト)と通し番号(リニアインデクス，リニアサブスクリプト)の相互変換は，sub2ind や ind2sub を使っても行える．

sub2ind は，多次元インデクスをリニアインデクスに変換する．

k = sub2ind( [4, 3, 2], 3, 1, 2 )
k = sub2ind( [4, 3, 2], 3, 4 )
k = sub2ind( [4, 3, 2], 15 )

はいずれも k = 15 を返す．

一方，ind2sub は，リニアインデクスを多次元インデクスに変換する．

>> i1 = ind2sub( [4, 3, 2], 15 )
i1 =
    15

>> [i1, i2] = ind2sub( [4, 3, 2], 15 )
i1 =
     3
i2 =
     4

>> [i1, i2, i3] = ind2sub( [4, 3, 2], 15 )
i1 =
     3
i2 =
     1
i3 =
     2

となる．

このプロセスは何次元配列にでも適用できる．つまり，N 次元配列を M 個の数字からなるサブスクリプトで参照するときは (ここでは M < N とするよ)， 3 次元配列を 2 個や 1 個の数字からなるサブスクリプトで参照したときと同様，余りの次元は全て展開されて並べられたものとして考えると良い．

ちなみに，インデクスのやりとりをするには，セル配列を使うのも便利だ．
セル配列の使い方については当ブログの記事

• matlab: セル配列で困っている人へ
• matlab: セル配列の参照法とその周辺

も参照されたし．

[matlab 記事インデックスへ]

matlab: セル配列で困っている人へ

matlab に慣れてきた人でも，「セル配列」をすっきり理解しきれた気がしない，という人は多いかもしれない．そうした向きへセル配列の解説を試みる．

まとめ 兼 目次

• セル配列とは「セル(=1×1のセル配列)」を並べた配列である．
• 「セル」は，任意のデータを中括弧で囲むと出来る．B = {A} など．
•  どんなタイプのものも中括弧 {} で囲むとセルになるということ．
• セル配列の作り方は，根源的には 1) データを中括弧 {} で囲む（大括弧 [] 同様，複数のデータを一度に囲っても良い），2) cell 関数を使って空のセル配列を作る，3) struct2cell を使って構造体から変換する，の３通りで，この他に num2cell, mat2cell, cellstr という変換関数を使う方法もある．もちろん，個々のセルを配列として並べて大括弧 [] で囲んでも良い．
• 小括弧 () を用いてアクセスすると，通常の配列と同様，部分配列を参照できる．
• 取り出された部分配列もセル配列．
• 中括弧 {} を用いてアクセスすると，セルの中身を取り出せる．
• インデクス(=サブスクリプト)の解釈(参照される要素がどれになるか)は小括弧と同様．
• 複数要素を参照すると，複数の文（式）を次々に実行（評価）したのと同じ結果になる．実行順序はリニアインデックス順．
• 大事なのでもう一度．中括弧 {} アクセスは，複数文相当の表現！
• セル配列に代入するには，次のいずれかの方法をとる．
• 左辺を小括弧指定，右辺をサイズの適合するセル配列にする．例： a(1:3) = {'foo', 100, true};
• 左辺を大括弧+中括弧指定，右辺を複数文相当の表現にする．例：[a{:}] = deal(100, 'strings', true);
• 複数文相当の表現の使い道
• 小括弧 () 受け：サブスクリプト，関数に入力引数を渡す
• 中括弧 {} 受け：セル配列のサブスクリプト
• 大括弧 [] 受け：連結，関数の出力引数を受け取る

注意

• 例によって根本は matlab 自体のドキュメントである．
• mathworks のサイトの matlab ヘルプの cell array のページ
• 特に，「セル配列って何？」の説明
• 本記事は，当ブログでそこそこ人気の記事「matlab: セル配列の参照法とその周辺」を大幅に拡充したものである．

では，始めよう．

matlab: 壁の外へ 1. 前置き

[執筆中，2013]

概要 

この一連の記事では，matlab での workspace 間のアクセス方法をいくつか挙げる．初回となる本記事では，matlab における workspace の概念について解説する．

• この文章は脱線上等のポリシーのもとで作成されています． 
• matlab において最強の文書は，matlab 付属のヘルプドキュメント群です．ともかく help を参照するに限る．ともかく matlab メニューの「ヘルプ→プロダクトのヘルプ」を読め！わからん言葉はコマンドプロンプトから >> help XX で読みまくれ！過去に色んな人に多くの質問を受けましたが，ほとんどはヘルプに書いてあるのを読んでいないだけです．ただ，ある種のことは別の切り口で見るのも有益な場合があります．本記事はそれを期待しています．

１．前置き： workspace, 関数，入出力引数

workspace について - 変数の見える範囲 -

matlab では，何か変数や関数があったとして，それらがある時点で「見える」かどうかは状況により変化する．例えば，matlab コマンドラインで変数 a に 3 という値を割り当てたとする．その後何らかの関数を呼び出しても，通常は，関数の内容によらず，関数の「中」からは「外」の，つまりコマンドラインで定義した a の存在や値は見えないし，関数内の処理は a の値に影響を与えない．

例えば，関数を func1 として，その中で変数 a に値 2 を割り当て，disp(a) するとする．
ファイル func1.m:
 function func1
 a = 2;
 disp('a in func1');
 disp(a);
この関数はコマンドライン上で見える変数 a にどう影響するだろうか？実行してみると，
>> a = 3; func1; disp('a in base workspace'); disp(a);
a in func1
     2
a in base workspace
     3
となる．func1 内での a への操作は，コマンドライン上から見える変数 a には影響を及ぼしていない．

これは，コマンドライン上で見える変数と，関数 func1 内で見える変数とが，別の世界に存在しており，互いに隔絶しているからである．一丁目の山田さんと三丁目の山田さんは別人なのである．この各々の状況で見えている世界のことを，matlab では workspace という．他のプログラミング言語では名前空間 namespace やスコープ scope などと呼ぶ．コマンドラインで見えている世界のことを，base workspace という．

関数 m-file と異なり，スクリプト m-file は，固有の workspace を持たず，呼び出し元と同じ workspace を共有する．同じご町内である．

入力引数・出力引数

では，コマンドラインの世界と関数内の世界で情報をやり取りするにはどうするか．それには，通常，入力引数と出力引数を用いる．

関数 func2 を
func2.m:
function y = func2(x)
y = 3 * x;
と作ったとすると，
>> a = 4; b = func2(a);
>> disp('a'); disp(a); disp('b'); disp(b);
a
     4
b
     12
となる．コマンドライン上での a の中身(=4)は，func2 での x の中身として渡され，func2 での y の中身(=12)は，コマンドライン上での b の中身として渡される．

なお，この受け渡しは，参照やポインタではなく，値のコピーが渡されることになっている（値渡し）．一丁目の連絡文書が三丁目に伝達されるとき，その文書自体が渡されるのではなく，三丁目の人が文書を見ながらせっせと書き写して町内会に持って帰るのだ．一丁目の言う文書 a と，三丁目の言う文書 x は，当初の内容は同じでも物理的には別の存在なのだ．このことと，前記 workspace 間の隔絶とによって，外の世界と中の世界の交流は，入力引数と出力引数という，いわば関数呼び出し国と呼び出され関数国間の出入国ゲートのみに限定され，各々の中では他者の干渉を受けない鎖国が保たれているのである．なお，呼び出され関数国から見た関数呼び出し国のことを，matlab では caller と呼ぶ．呼び出し国の workspace のことを，caller workspace という．コマンドラインから関数を呼び出せば，関数内から見た caller workspace とは base workspace なわけだ．

workspace の効果

異なる workspace にアクセスできないことが何を意味するか．それは，裏を返せば，見えない範囲からの干渉は基本的に生じ得ないことを意味する．そのため，今自分の見えている範囲で整合性さえ取れれば大丈夫，という保証になるのだ．つまり，「今書いている関数 m-file がまともに動くこと」を積み重ねていけばいい，という世界が保証されている．一度に管理するのはエディタのそのページくらいで済むわけだ．

前置きまとめ

• matlab では，ある時点で「見える」変数の範囲を workspace という．
• コマンドラインから見える範囲は base workspace．
• 関数を呼び出せば，関数内の workspace は base workspace とは別物．スクリプトを呼び出すと，workspace は元のを共有．
• 異なる workspace 間は隔絶した鎖国関係にある．同じ名前の変数でも別々のものを指す．
• 関数呼び出し元の workspace (=caller workspace) と呼び出された関数内の workspace とは，通常，入出力引数という出入国ゲートを介してのみ，値渡しでやり取りが成立する．
• workspace が分かれているおかげで，目に見える範囲の整合性さえ確保すれば，おかしな矛盾が生じない仕組みになっている． 
• 関数のアクセススコープも workspace によるのだが，matlab サーチパスも含めてよりややこしくなりそうなので本記事では割愛．

メモ

• 関数 m-file とスクリプト m-file についてはこれを読もう：matlab のヘルプ→プロダクトのヘルプ→MATLAB->Getting Started->Programming->Scripts and Functions
• mathworks のサイトでも読めるぞ．
• workspace についてはこれを．同じく MATLAB->Desktop Tools and Development Environment->Workspace, Search Path, and File Operations
• mathworks のサイトでも読めるぞ．base workspace と関数内の workspace についても説明がある．
• workspace は用語でもあるが，matlab で実際に workspace browser なるものを呼び出すコマンド名でもあるぞ．
• matlab の次の性質は知っておくといいかもしれない．もっとも，プログラムの理論的な動作上は気にする必要は無い．実は matlab は毎度律儀に値渡しをメモリコピーの意味で実践してはいないようだ．大き目の行列を使って，matlab プロセスのメモリ使用量を観察しているとわかるのだが，入力引数を読むのにしか使わず，その中身を書き換えない場合は，コピー分のメモリ確保は生じていないようだった（matlab のバージョンには依存するかも）．これ，気づく人は結構気づくものなんだね．Copy on Write (COW, or Lazy Copying) と言う動作で，Internal Matlab memory optimizations | Undocumented Matlab に詳しいぞ．
• 構造化プログラミングって昔よく聞いたね．ベーマガで Dr.D も言ってたぞ．

予告

 

workspace の安寧の中で虚偽の繁栄に酔いしれるプログラマたち．重い計算を恣（ほしいまま）にループで回す彼らは，長時間の計算中，突然天からの Ctrl+C による襲撃にさらされる．workspace 間の隔絶が仇になり，積み重ねた計算結果は日の目を見ることなく関数 workspace 内で雲散霧消してしまう．この run で受けた生を無駄にしたくない．天寿を全うできない場合でも，何とか壁の外の caller にこの情報を届けたい．越境者の自由を求め，今，餓狼たちの挑戦が始まる．次回，第２節，「workspaceをまたいだアクセス」．

※matlab では，実行中の処理は Ctrl+C で中断させることができます．Ctrl+C を押しても，あるまとまった処理が終わるまで待たされる場合もあります．
※本記事は2013年春季にノリノリで執筆しております．

matlab: limitations for large/small numbers

Here are some limitations for large/small numbers in Matlab.
So it's about boundary and precision.
Examples below are results derived in Matlab R2007b (64bit) on Win7 x64 on Core i5 2500.
Just a list of examples; no summary or interpretation as of now.

**
>> help eps
will explain most of the following results
**

Large number limit

 

number:

>> realmax
ans =
  1.7977e+308

>> realmax+9.97e+291
ans =
  1.7977e+308

>> realmax+9.98e+291
ans =
  Inf

exp:

>> exp(709)
ans =
  8.2184e+307

>> exp(710)
ans =
  Inf

log:

>> log(realmax+9.97e+291)
ans =
  709.7827

>> log(realmax+9.98e+291)
ans =
  Inf

Small number limit

number:

>> realmin
ans =
  2.2251e-308

>> isequal(0,realmin)
ans =
     0

>> isequal(0,0+realmin)
ans =
     0

>> 1+realmin
ans =
     1

>> isequal(1,1+realmin)
ans =
     1

>> isequal(1,1+realmin*4.98e+291)
ans =
     1

>> isequal(1,1+realmin*4.99e+291)
ans =
     0

>> realmin*4.99e+291
ans =
  1.1103e-016

>> eps
ans =
  2.2204e-016

>> isequal(2,2+realmin*9.97e+291)
ans =
     1

>> isequal(2,2+realmin*9.98e+291)
ans =
     0

>> isequal(100,100+realmin*3.1e+293)
ans =
     1

>> isequal(100,100+realmin*3.2e+293)
ans =
     0

exp:

>> exp(-745)
ans =
  4.9407e-324

>> exp(-746)
ans =
     0

>> isequal(1,exp(realmin*4.98e+291))
ans =
     1

>> isequal(1,exp(realmin*4.99e+291))
ans =
     0

log:

>> log(realmin*1.0e-15)
ans =
 -742.8306

>> log(realmin*1.0e-16)
ans =
  -Inf

>> isequal(log(1),log(1+realmin*4.98e+291))
ans =
     1

>> isequal(log(1),log(1+realmin*4.99e+291))
ans =
     0

>> isequal(log(2),log(2+realmin*9.97e+291))
ans =
     1

>> isequal(log(2),log(2+realmin*9.98e+291))
ans =
     0

>> isequal(log(100),log(100+realmin*4.1e+294))
ans =
     1

>> isequal(log(100),log(100+realmin*4.2e+294))
ans =
     0

matlab: find いらずの一発参照・代入

matlab では，配列の要素のうち特定の条件を満たすものを参照する際，ものすごく効率的な方法が用意されている．それが Logical Subscripting だ(logical indexing という人もいる)．基本的な説明は Matlab の help の Getting Started や mathworks のページにも載っているが，もう少し掘り下げてみよう．

ここでいう Logical とは，論理値配列の意味の論理である．早速やってみよう．

例： 行列 A のうち，値が 0.5 未満の要素は 0 にする．閾値処理ですな．
>> A = rand(3)
A =
    0.8147    0.9134    0.2785
    0.9058    0.6324    0.5469
    0.1270    0.0975    0.9575

>> A(A<0.5) = 0
A =
    0.8147    0.9134         0
    0.9058    0.6324    0.5469
         0         0    0.9575

その通り，find は使わなくてよいのである．今まで使いまくってたよ‥
A(find(A<0.5)) = 0
とか書いてたよ‥．はずかちー．

[多次元インデクス自体については matlab: 一目で分かる多次元インデクス をどうぞ]

もちろん，A の参照に A に関する論理式を使う必要はない．他のものでも使える．
例：同サイズの別配列による条件づけ
>> Score = [60,55,80]
Score =
    60    55    80

>> Effort = [50,75,60]
Effort =
    50    75    60

>> Score(Effort>70) = Score(Effort>70) + 20
Score =
    60    75    80

頑張った人に+20点とかね．Score(Effort>70)が2回出てくるのがなんとも野暮ったいけど，何とかならないだろうか？以下，速度と効率を追求してみよう．

* 速度の検討

logical は比較的小さい容量になる（1byte/element）ので，保存して使いまわしてもよさそうだ．

例： 大きな配列の論理参照における論理配列の使いまわしの効果


>> S=rand(5000);S1=S;N=10;T=zeros(N,1);
>> for n=1:N,tic;S1(S>0.5)=S(S>0.5)+2;T(n)=toc;end;mean(T)
ans =
    1.7190

>> S=rand(5000);S1=S;N=10;T=zeros(N,1);
>> for n=1:N,tic;LS=S>0.5;S1(LS)=S(LS)+2;T(n)=toc;end;mean(T)
ans =
    1.3880

2割以上速くなる．ちなみに，自分の過去のやり方(findする)だと

>> S=rand(5000);S1=S;N=10;T=zeros(N,1);
>> for n=1:N,tic;IS=find(S>0.5);S1(IS)=S(IS)+2;T(n)=toc;end;mean(T)
ans =
    1.6312

あれ？意外と速い．個数の問題か？




例： 参照個数による処理時間の変化(find vs. logical subscripting)


>> S=rand(5000);S1=S;N=10;T=zeros(N,1);
>> for n=1:N,tic;IS=find(S>0.9);S1(IS)=S(IS)+2;T(n)=toc;end;mean(T)
ans =
    0.6060

うおっ．てきめん．一方の logical subscripting はどうか？

>> S=rand(5000);S1=S;N=10;T=zeros(N,1);
>> for n=1:N,tic;LS=S>0.9;S1(LS)=S(LS)+2;T(n)=toc;end;mean(T)
ans =
    0.6400

>> S=rand(5000);S1=S;N=10;T=zeros(N,1);
>> for n=1:N,tic;LS=S>0.1;S1(LS)=S(LS)+2;T(n)=toc;end;mean(T)
ans =
    1.1668

な　ん　で　だ　よ
なんで速度に差が出るのよ．たとえば，裏で find や sparse 的なことをしてて，対象個数が少ないほうに合わせてどうこうしてるのか？

条件に引っかかる個数を減らしながら，かかる時間をプロットしてみよう．

例： 参照個数による処理時間の変化（Logical Subscripting）

>> S=rand(5000);S1=S;N=4;LN=15;LVLs=linspace(0,1,LN);T=zeros(N,LN);
>> for ln=1:LN,for n=1:N,tic;LS=S>LVLs(ln);S1(LS)=S(LS)+2;T(n,ln)=toc;end;end;
>> figure;plot(LVLs,mean(T,1),'o-','linewidth',3);set(gca,'fontsize',20);grid on;
>> ylabel('elapsed time');xlabel('the threshold used');
>> print(gcf,'-dpng','logical_subscripting_etimes');

そうしてできたプロットがこれである．

参照される要素数が，
多い（グラフ左）：半分より多くなるほど線形に速くなるが，係数は小さい．
半分（グラフ中央）：一番遅い
少ない（グラフ右）：半分より少なくなるほど，おおむね比例的に速くなる．
という結果が得られた．

例： 参照個数による処理時間の変化（find）

>> S=rand(5000);S1=S;N=4;LN=15;LVLs=linspace(0,1,LN);T=zeros(N,LN); >> for ln=1:LN,for n=1:N,tic;IS=find(S>LVLs(ln));S1(IS)=S(IS)+2;T(n,ln)=toc;end;end;
>> figure;plot(LVLs,mean(T,1),'o-','linewidth',3);set(gca,'fontsize',20);grid on;
>> ylabel('elapsed time');xlabel('the threshold used');
>> print(gcf,'-dpng','find_etimes');
find でやると素直に

要素が少なくなると処理が早く済む，となる．
二つの結果を重ねると

基本的には logical subscripting を使っていればよさそうだ．要素が少なくなってくると find がわずかに早いこともあるが，その差は微小だ．


logical subscripting の挙動の理由はいまいちわかりませんが，パフォーマンスにこだわる方は参考にしてみてください．ちなみに，R2007b on Win7 x64 on core2duo (mobile) です．

*プログラミング上の性質の検討


この Logical Subscripting の要点は，まさに logical 型の配列を括弧内に与えることにある．カッコ内に与えられた配列が logical かどうかはちゃんとチェックされ，数値 0 と論理値 false は区別される．
例：0 と false の区別
>> a=1
a =
     1

>> a(0)=3
??? Subscript indices must either be real positive
integers or logicals.

>> a(logical(0))=3
a =
     1

>> a(logical(1))=3
a =
     3

また，少なくとも R2007b では，参照に使う論理配列は，元の配列より小さくてもよいようだ．

例：参照に使う論理配列と入れ物の配列の大小

>> B=reshape([1:6],[3,2])
B =
     1     4
     2     5
     3     6

>> B(logical([0,0,0,1]))=9
B =
     1     9
     2     5
     3     6

逆に，大きいと怒られる．

>> B(logical([0,0,0,0,0,0,1]))=9
???  In an assignment  A(I) = B, a matrix A cannot
be resized.

ご参考まで．

-- 当ブログの matlab 関係へはこちらからどうぞ．--

matlab: なんちゃって構造体

空の構造体は

>> a = struct([])
a =
0x0 struct array with no fields.

とすれば簡単に作れる．では，「フィールドを持たないけど空ではない構造体」はどうすれば作れるか？

すぐに思いつくのは，rmfield を使って，フィールドを消して行くやり方だろう．

>> b.foo = []; a = rmfield(b, 'foo')
a =
1x1 struct array with no fields.

これはかっこ悪いので却下．だが実現は出来ることがわかる．
では，rmfield はどうやってこれを実現しているか？
rmfield は実は組み込み関数でなく .m 関数なので，type rmfield して中身をのぞくと，どうやら cell2struct を呼び出しているらしい．
ならばはじめから cell2struct を使えばいいのである．しかも cell2struct は組み込み関数だから，わりと根源的なやり方といえるだろう．

>> a = cell2struct(cell(0,1),{})
a =
1x1 struct array with no fields.

おお，できた．これでいいや．

しかし，恒例の whos でみると

>> whos a
  Name      Size                    Bytes  Class
  a         1x1                         0  struct array
Grand total is 0 elements using 0 bytes

のように，メモリが割り当てられていないではないか．しかも a の repmat で b を作ってみても，b も Bytes は 0 のままである．

実際にメモリを割り当てるにはフィールドを定義しないといけないらしい．AMD64 な Linux 上で AMD64 な MatlabR2006a を使った場合，

>> b = repmat(a,[100,1]);
>> whos b
  Name      Size                    Bytes  Class
  b       100x1                         0  struct array
Grand total is 0 elements using 0 bytes
>> b(1).foo = '';
>> whos b
  Name      Size                    Bytes  Class
  b       100x1                       952  struct array
Grand total is 100 elements using 952 bytes

となるが，同じことを IA32 な Windows 上で IA32 な Matlab6.5.1 を使って行った場合，

>> b=repmat(a,[100,1]);
>> whos b
  Name      Size                   Bytes  Class
  b       100x1                        0  struct array
Grand total is 0 elements using 0 bytes
>> b(1).foo = '';
>> whos b
  Name      Size                   Bytes  Class
  b       100x1                      520  struct array
Grand total is 100 elements using 520 bytes

となる．なんで違うんでしょう．

matlab: なんちゃってcomplex

matlab (今回は 6.5.1 でしか試してません)は複素数について直感に反する挙動を示すことがある．例を見てみよう．

>> A = repmat(complex(0),[10,1]);
>> whos A
  Name      Size                   Bytes  Class
  A        10x1                      160  double array (complex)
Grand total is 10 elements using 160 bytes
>> A(end) = complex(1,0);
>> whos A
  Name      Size                   Bytes  Class
  A        10x1                       80  double array
Grand total is 10 elements using 80 bytes

complex じゃなくなってる．およ？一応確認すると，A に代入した値は実数でなく複素数と認識されているはずだ．

>> clear a; a=complex(1,0); whos a
  Name      Size                   Bytes  Class
  a         1x1                       16  double array (complex)
Grand total is 1 element using 16 bytes

となるからだ．ただし，入力引数に対する inputname のように，出力引数についての情報取得機構が関数 complex 内で働いているとすると，よくわからない話になるかもしれない．

一つ考えられる仮説は，「虚部が全て 0 の配列は実数型にされてしまう」というものだ．実際，その方がメモリ効率やら計算効率やらよくなる気はする．

念のため，虚部が 0 でない複素数を代入すれば，配列は実数型に変換されることは無い．

>> A = repmat(complex(0),[10,1]);
>> A(end)=complex(0,1);
>> whos A
  Name      Size                   Bytes  Class
  A        10x1                      160  double array (complex)
Grand total is 10 elements using 160 bytes

そして，要素を虚部が 0 の複素数に順に置き換えていくと，

>> A=complex([1,2],[3,4]);
>> A(1)=complex(real(A(1))); whos A
  Name      Size                   Bytes  Class
  A         1x2                       32  double array (complex)
Grand total is 2 elements using 32 bytes
>> A(2)=complex(real(A(2))); whos A
  Name      Size                   Bytes  Class
  A         1x2                       16  double array
Grand total is 2 elements using 16 bytes

となった．ここまでは前述の仮説で当たっていそうである．

しかし，もちろん話はこれで終わらない．次の例を見てみよう．

>> A=zeros(1,2);
>> A(:)=complex(0);
>> whos A
  Name      Size                   Bytes  Class
  A         1x2                       32  double array (complex)
Grand total is 2 elements using 32 bytes

なんと配列は複素数のままであり，これは前述の仮説と矛盾する．ところがさらに，代入の方法を変えると，

>> A=zeros(1,2);
>> A(1:end)=complex(0);
>> whos A
  Name      Size                   Bytes  Class
  A         1x2                       16  double array
Grand total is 2 elements using 16 bytes

となって前述の仮説が復活する．どうなってるのという話である．

ここで新たに，｢全代入(whole substition とか？)｣という概念を考えてみる．
全代入と考えるのは次の４つ．

1. A = B;
2. A(:) = B;
3. A = B(:);
4. A(:) = B(:);
   

変数間の代入で虚部 0 でも complex を維持するには，どうやら厳密にこの４形式しか許されないようだ．
たとえば２次元配列で A(:,:) としてもダメだったし，列ベクトルで A(:,1) としてもダメだった．
この概念を利用して立てた仮説がこれだ：
「全代入では愚直な代入が行われ，それ以外では値の取り出し・代入いずれかで可能なら虚部が 0 の複素数は実数に落ちる」
いまのところこの仮説で説明できるような気がする．
この他に，左辺が全参照（全代入と同じ気分）で右辺が「適当な」式なら，やはり実数に落ちないようだ．以下の例を参照のこと．

>> A=complex(rand(3));
>> A(:)=complex(0); whos A
  Name      Size                   Bytes  Class
  A         3x3                      144  double array (complex)
Grand total is 9 elements using 144 bytes
>> A(:)=complex(zeros(3)); whos A
  Name      Size                   Bytes  Class
  A         3x3                      144  double array (complex)
Grand total is 9 elements using 144 bytes

いや～．重箱の隅つつくといろいろ出てくるなぁ．
ちなみに，大方の予想通り

>> A=complex(rand(3));
>> A(:,:)=complex(0); whos A

  Name      Size                   Bytes  Class
  A         3x3                       72  double array
Grand total is 9 elements using 72 bytes
>> A=complex(rand(3));
>> A(:)=complex(zeros(3)) + 1; whos A
  Name      Size                   Bytes  Class
  A         3x3                       72  double array
Grand total is 9 elements using 72 bytes

となる．もういいって？
今回の現象は，subsref, subsasgn, subsindex と型キャストの実装のコンボで生じるんだろうけど，意図された仕様なんですかねぇ．

matlab: 呼び出しの調査

matlab 上で関数を呼び出すと疑似コンパイル(p-code)され，メモリに残る．
メモリ上の p-code buffer に存在する m-file の名前を見るには，inmem という関数を使う．
これで，実際に呼ばれた関数の一覧が得られる．

関数のコードを走査して再帰的に依存を調べるには，depfun という関数を使う．
depfun の賢さは限定的なので，もれなくわかると期待してはいけない．depdir も同様．

profile なんかも結構役に立つ．単に echo してみるのも面白いかもしれない．

matlab: end+1

Matlab では，配列の最終要素のインデクスを表す end という指定法がある．
この end にまつわる小ネタ．

配列に次々要素を追加していく，という処理をスクリプト中に書く場合，番号を決め打ちする書き方

a(1) = 1;
a(2) = 2;
a(3) = 3;

をすると，順番を入れ替えたいときや，間に要素を挟みたいときに一大事になる．

そこで，こんな書き方が楽チンだ．

a(1) = 1;
a(end+1) = 2;
a(end+1) = 3;

これなら入れ替えや追加をコピペで簡単に行える．

ちなみにこう書いた場合，リストは横に，つまり2次元目に伸びていく．上の例では，a は 1 x 3 のサイズになる．
もともとが 3 x 1 だったりした場合，4 x 1 に伸びる．もともとが 3 x 2 だったりした場合，

>> a=rand(3,2);
>> a(end+1)=3
???  In an assignment  A(I) = B, a matrix A cannot be resized.

のように怒られる．この場合でも，伸ばす次元が明示的になるような書き方，すなわち

>> a=rand(3,2);
>> a(3,end+1)=3
a =
    0.1934    0.5417         0
    0.6822    0.1509         0
    0.3028    0.6979    3.0000
>>

のようにすれば，各次元についての end としてちゃんと解釈される．

実は上記の怒られ方は end に限ったことではない．つまり，

>> a=rand(3,2);
>> a(7)=3
???  In an assignment  A(I) = B, a matrix A cannot be resized.

とこういうことなのだ．しかし，

>> a=rand(3,2);
>> a(5)=3
a =
    0.3412    0.3093
    0.5341    3.0000
    0.7271    0.5681

とこういうことでもある．ここで a(end-1) としてももちろん同じ．結局，

• linear index で問題ないものはそのように解釈される
• linear index で曖昧さが出る状況では怒られる
• end は，それが解釈される次元における末尾のインデクスの数字と単純に置き換え
  

というルールだと思っていいのだろう．

補足．もちろん，end+1 に限らず end+10 のような指定も可能です．

matlab: 構造体配列の参照法とその周辺

構造体配列も，普通の配列（数値配列や文字配列）とは少し違う．
どちらかというとセル配列ライクである．matlab: セル配列の参照法とその周辺を参照．
試してみよう．まず，構造体配列を作る．

>> a=struct('f',{'foo','bar'})
a =
1x2 struct array with fields:
    f

次に，小括弧で普通に参照してみる．

>> a(1)
ans =
    f: 'foo'

構造体が一つ返ってくる．さらに，インデクスリストで参照してみる．

>> a(:)
ans =
2x1 struct array with fields:
    f

構造体配列が返される．コロン単独 : で参照すると一次元ベクトル(縦に並んだ配列)が返されるのは，
Matlab のお約束の振る舞いである．ちなみに，1:end で参照すると，振る舞いが異なる．この話はいずれ．
問題はここから．構造体配列のフィールドを参照する場合の話だ．
単一要素の参照なら，

>> a(1).f
ans =
foo

であって全く普通だ．しかし！インデクスで参照してフィールドを指定すると，

>> a(:).f
ans =
foo
ans =
bar

となり，また単にフィールドを指定した場合にも

>> a.f

ans =

foo

ans =

bar

のように，セル配列の記事でも述べた，例の「複数の文を評価した状態」になる．
あとの話はセル配列のときと同様になるわけだが，構造体配列でありがちな例を挙げておこう．
例：一件のレコードが一つの構造体になっているような構造体配列で，特定のものを見つける．

>> clear a b c
>> a.label='small';a.value=100;b.label='small';b.value=200;c.label='large';c.value=1000;d=[a,b,c]
d =
1x3 struct array with fields:
    label
    value
>> small_ones = find(strcmp({d.label},'small'))
small_ones =
     1     2
>> under150 = small_ones(find([d(small_ones).value]<150))
under150 =
     1

てな感じ．連結 {} や [] を使って単純に書ける．
関数 deal を使うともっと幅が広がるかもしれない．例えば，上で定義した構造体配列 d に対して

>> [d.label]=deal('tiny','little','huge');d.label
ans =
tiny
ans =
little
ans =
huge

なんてこともできる．

matlab: セル配列の参照法とその周辺

※本記事を大幅に拡充して「matlab: セル配列で困っている人へ」を執筆しました．2013.06.10

セル配列は Matlab の配列の中でもやや異色である．
最も大きな特徴は，他の任意の型の変数（配列）を要素として扱えることだろうが，
もう一つ，中括弧がセル配列専用の参照演算子 {} として用意されているのも独特だ．
しかも，通常の参照演算子である小括弧 () でもアクセスできる．
そして，小括弧でアクセスした場合と中括弧でアクセスした場合には結果が異なるのである．
根本的には，小括弧は「セル配列の部分配列を取り出す」，中括弧は「要素セルに格納された中身を取り出す」と考えられる．
試してみよう．まず，セル配列を作る．

>> a={1,2};

次に，小括弧で普通にアクセスしてみる．

>> a(1)
ans =
    [1]

このとき，a(1) 自体はセル配列である．つまり，b=a(1) とすれば，b は 1 by 1 のセル配列となる．
さらに，インデクスリストでアクセスしてみる．

>> a(:)
ans =
    [1]
    [2]

ここでももちろん a(:) 自体はセル配列である．つまり，b=a(:) とすれば，b は 2 by 1 のセル配列となる．
では，中括弧でアクセスしてみよう．

>> a{1}
ans =
     1

このとき，a{1} は数値である．さらに，インデクスリストでアクセスしてみる．ここからが面白いわけですよ．

>> a{:}
ans =
     1
ans =
     2

なんと ans が 2 つ．つまり 2 つのセンテンスが評価されたことがわかる．
中括弧使用時のインデクスリストによるアクセスとは，一見単独の文に見えながら，複数の参照文をシリアルに実行しているのに他ならないのだ．
ということは，である．当然，この参照結果を一つの変数にそのまま代入することはできない．

>> b=a{:}
??? Illegal right hand side in assignment. Too many elements.

ただし b=a{1} のような単一要素の参照ならもちろん代入できる．
単純な代入が行えない，となると，中括弧でのインデクスリスト参照にはどんな使い道があるのか？
これがちゃんとあるのである．
matlab には，複数の文の評価結果を受け取る書式が存在する．
それは，小括弧 () と中括弧 {} と大括弧 [] の括弧 3 兄弟である．
ここでの小括弧 () は，関数呼び出しの引数指定に使う場合と，配列を参照する場合のものを指す．
例えば，

>> a={2,2};
>> b=[1,2;3,4];
>> b(2,2)
ans =
     4
>> b(a{:})
ans =
     4

となる．つまり，b(2,2) という指定の本質は，
「括弧内に 2 つの文があって，それぞれの評価結果が 2 になっている」
ということなのだ．だから，複数の文の評価をもたらす a{:} で置き換えることができた．
関数引数も同じ理屈だが，こちらはすぐにわかる実用性がある．varargin を利用して

function a = myfunc(varargin)
a = someotherfunc(varargin{:});

とすれば，即座にラッパーが書けるのだ．これは応用しやすい．
別に varargin でなくとも，一連の関数引数候補をセル配列にパックしておけば便利．多分．時々．
中括弧や大括弧は，複数の要素を連結するときに使う．
中括弧で連結すればセル配列になり，大括弧で連結すれば元の要素の属するクラスの配列になる．
このとき，要素を何で区切るか？といえば，カンマ , やセミコロン ; なわけで，これらも文の区切り記号であるところに注目である．
もっとも，連結するという文脈内では，カンマは 2 次元目連結，セミコロンは 1 次元目連結を示唆するものであって，
コンソールへの出力をする・しないという通常の文区切り時とはもちろん働きが違う．
セル配列の中括弧アクセスの場合，連結の括弧の中では全てカンマ区切りにみなされるらしい．
まとめると，セル配列への中括弧でのインデクスアクセスは，あくまで複数の参照文のシリアルな実行であるが故に，
元の配列の次元構成を保存できない．つまり，

>> a={1,2;3,4}
a =
    [1]    [2]
    [3]    [4]
>> {a{:,:}}
ans =
    [1]    [3]    [2]    [4]

となってしまうので気をつけよう．
多次元構成が保存されないことを除いてはなかなか便利に使える．例えば，数値のセル配列を数値配列にするには

b = [a{:}];

でいいわけだ．また，次元に関しては，代入先を予め設定しておいてから

b = zeros(size(a));
b(:) = [a{:}]'; %この転置はいるのか？手元の R2006a ではあってもなくてもエラーも出ないし結果も同じだった

のように代入するか，あるいは reshape を使って

b = reshape([a{:}],size(a));

などとすればよいだろう．
さて，セル配列のほかに不思議な振る舞いをする配列といえば構造体配列があるわけだが，こちらはまたの機会に．
こういった参照や連結の背景には subsref, subsasgn や subsindex といった参照用関数，そして cat, horzcat や vertcat といった連結用関数がある．
これらは matlab 上でオブジェクト指向プログラミングをしない限りあまり目にしないかもしれないが，背景として理解しておいてもいいと思う．

2013.06.10
本記事を大幅に拡充して「matlab: セル配列で困っている人へ」を執筆しました．

matlab: 対角行列を効率的に生成

Matlab で対角行列の類を生成する場合がよくある．これを，メモリ・速度・コード量とも効率よく行うにはどうすればいいかという話．
普通の対角行列の場合．対角成分を並べたベクトルを w とすると，

W = diag(sparse(w));

が多分最良．
ブロック対角行列の場合はどうか．m by n の行列が k 個対角に並ぶ場合を考える．ブロック対角成分を並べた行列を B (m by n*k) とすると，

C1 = sparse(B);
C1(m*k,n*k) = 0;
C2 = reshape(C1,[m*k*n,k]);
C2(m*(k*n+1),k) = 0;
C = reshape(C2, [m*k,n*k+1]);
C(:,end) = [];

くらいでどうだろう．
blkdiag は引数の形式が馴染まないことが多いと思う．セル配列から展開するのか？あとは，単に自分が sparse や spdiags を使いこなせていないという噂もある．
※最後の行 20070424 に訂正した．bad->C(:,end-(n-1):end) = [];

行列からのブロック対角成分の抽出や，抽出した成分を用いてブロック対角行列を再構成する方法について，「matlab: ブロック対角成分の抽出」という記事も書きました．