TIM Labs
偏微分方程式を解くために、まず、行列のreshapeを理解しないといけない。
評価したい対象領域をN*Nに区切ったので、対象領域はN*Nの行列(2次元の配列)になっている。実際には温度のN*N配列である。

しかし、温度が行列のままだと、C * T = B という感じで、偏微分を差分表現したのを行列演算に持ち込めない。
なお、Cは係数行列、Tは温度、Bは差分表現の右辺あるいは境界条件である。

そこで、行列を列ベクトルにすることを考える。

行列を整形する関数は非常にたくさんあるのだが、今回使うのはreshape()である。

>> A = 1:9
A =

   1   2   3   4   5   6   7   8   9

>> A = reshape(A,3,3)
A =

   1   4   7
   2   5   8
   3   6   9

>> A = reshape(A,9,1)
A =

   1
   2
   3
   4
   5
   6
   7
   8
   9

>>
したがって、 C * T = B において、元々2次元配列であったTとBを1次元列ベクトルにできるはずだ。

求めたいのは温度Tなので、C * T = B の両辺にCの逆行列をかけると次のようになる。

T = C \ B

差分表現の右辺あるいは境界条件をN*Nの行列として作る。 外周以外は、偏微分方程式が成り立ち右辺が0なので、zeros()で全部0のN*Nの配列を作る。 外周は温度を外気温T0=300を指定するのだが、過熱している辺についてだけ、sinカーブになるようにする。 そして、最後に、N*Nの列ベクトルbを作っている。
B = zeros(N,N);
B(1,:)= B(N,:)= B(:,N) = T0;
B(1:N,1) = sin(pi*x/L)*(TH-T0) + T0;
b = reshape(B,NN,1);

これに合わせて、前回 gallery("poisson",N) で作った行列うち、上記の行列Bで温度を指定した箇所に対応する行を、単位行列の行と同じにするために、単位行列の対応行をコピーするという強引な方法をやっている。

C = full(gallery("poisson",N));
I = eye(NN);
C(1:N,:) = I(1:N,:);
C(1:N:NN,:) = I(1:N:NN,:);
C(N:N:NN,:) = I(N:N:NN,:);
C(NN-N+1:NN,:) = I(NN-N+1:NN,:);
C = sparse(C);

N=4の場合の係数行列は次のようになる。
前回と比較すると、対角線上に1が並んでいる箇所が元の領域の外周に相当する行を示す。
C =

   1   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
   0   1   0   0   0   0   0   0   0   0   0   0   0   0   0   0
   0   0   1   0   0   0   0   0   0   0   0   0   0   0   0   0
   0   0   0   1   0   0   0   0   0   0   0   0   0   0   0   0
   0   0   0   0   1   0   0   0   0   0   0   0   0   0   0   0
   0  -1   0   0  -1   4  -1   0   0  -1   0   0   0   0   0   0
   0   0  -1   0   0  -1   4  -1   0   0  -1   0   0   0   0   0
   0   0   0   0   0   0   0   1   0   0   0   0   0   0   0   0
   0   0   0   0   0   0   0   0   1   0   0   0   0   0   0   0
   0   0   0   0   0  -1   0   0  -1   4  -1   0   0  -1   0   0
   0   0   0   0   0   0  -1   0   0  -1   4  -1   0   0  -1   0
   0   0   0   0   0   0   0   0   0   0   0   1   0   0   0   0
   0   0   0   0   0   0   0   0   0   0   0   0   1   0   0   0
   0   0   0   0   0   0   0   0   0   0   0   0   0   1   0   0
   0   0   0   0   0   0   0   0   0   0   0   0   0   0   1   0
   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   1
 

次の図で見たほうが分かりやすいだろう。
poisson2d-6.png

以上で、一通りの説明を終えたので、次回はプログラムを示す。
前回、偏微分方程式の数値解の結果を示したが、これをどうやって計算したかについて説明していこうと思う。
まず、対象とする一辺Lの正方形領域をメッシュに切ろう。頂点がN*N個になるように切って計算したのが前回の図である。すると、格子の一辺の長さΔはL/(N-1)なのだが、そんなものは全体に共通なので係数を掛けたり割ったりだけに過ぎず本質的でないので、格子の一辺の長さは1としよう。これで計算が楽になる。

次に、2次のポアソンの偏微分はどうすれば計算できるかを考える。微分、偏微分では、隣の点との差分を計算すればよく、2階編微分になっているので、差分の差分をとれば大丈夫ということで、次図のイメージで計算すればよいと思われる。

poisson2d-4.pngところで、上図の各格子点に対して、上下左右の隣点の値(温度)の和から対象点(i,j)の値(温度)の4倍を引けば2階編微分に相当する値が出るはずである。今回は、この値が0で、方程式が格子点(i,j)に対して1つ出来上がる。
そして、こういう点がn*n個あるので、方程式の数がn*n個と膨大になる。全ての格子点の温度T(i,j)が変数なので、変数の数がn*n個。そして、方程式の数がn*n個だと、これはn*n元の連立一次方程式になり解けるはずなのだ。

係数行列をCとすると、Cはサイズn*nの正方行列になるわけで、格子点(i,j)に対する式は、i+(j-1)*n 行目の係数で表される。この点については、どういう風に計算するかの説明が必要なのだが、境界条件の説明するときに合わせて行う。


この係数行列は規則的なので作ろうと思えば作れるのだが、Octaveは非常に多数の便利な行列を自動生成してくれる機能があり、そのための関数がgallery()である。

>> full(gallery("poisson",4))
ans =

   4  -1   0   0  -1   0   0   0   0   0   0   0   0   0   0   0
  -1   4  -1   0   0  -1   0   0   0   0   0   0   0   0   0   0
   0  -1   4  -1   0   0  -1   0   0   0   0   0   0   0   0   0
   0   0  -1   4   0   0   0  -1   0   0   0   0   0   0   0   0
  -1   0   0   0   4  -1   0   0  -1   0   0   0   0   0   0   0
   0  -1   0   0  -1   4  -1   0   0  -1   0   0   0   0   0   0
   0   0  -1   0   0  -1   4  -1   0   0  -1   0   0   0   0   0
   0   0   0  -1   0   0  -1   4   0   0   0  -1   0   0   0   0
   0   0   0   0  -1   0   0   0   4  -1   0   0  -1   0   0   0
   0   0   0   0   0  -1   0   0  -1   4  -1   0   0  -1   0   0
   0   0   0   0   0   0  -1   0   0  -1   4  -1   0   0  -1   0
   0   0   0   0   0   0   0  -1   0   0  -1   4   0   0   0  -1
   0   0   0   0   0   0   0   0  -1   0   0   0   4  -1   0   0
   0   0   0   0   0   0   0   0   0  -1   0   0  -1   4  -1   0
   0   0   0   0   0   0   0   0   0   0  -1   0   0  -1   4  -1
   0   0   0   0   0   0   0   0   0   0   0  -1   0   0  -1   4


第1引数に、行列の種類を文字列で指定する。"poisson"のとき、第2引数に元の分割数n(4)を指定すると、サイズがn*n(16)の正方行列が作られる。
full()は、gallary()が出力する行列がスパースな行列なので、そのまま表示すると、0以外が入っている箇所しか示さないので、普通の行列に展開するために使った。よく使われる行列の多くは、ほとんどの要素が0というのが多く、値の入っているところだけデータを持つことで、メモリが非常に圧縮できるし、計算も早くなったりする。

最初の図に従うと、点(i,j)の係数は-4で、隣の点は1でなければいけないのだが、この行列は符号が逆転している。でも、どうせ偏微分の結果は0なのだから、符合が逆転しても結果は同じになるので、このまま使ってしまおう。


配列は大きくなると、一覧性がどんどん悪くなる。それで、配列の状態を一瞬で把握できるような関数も用意されている。

>> colormap("gray")
>> imagesc(full(gallery("poisson",6)))
poisson2d-5.png
6x6の格子点に対するポアソンの行列なので、36次の正方行列になった。
対角線の白マスが4を示し、グレイが0、黒が-1になっている。
上下左右の4点を指定しようとしても、外周上の格子点の場合、4つの隣点のうちの少なくとも1つは存在しない点になるため、横の行を見ると、黒が本来4箇所ないといけないのだが、3箇所や2箇所になっている場所がある。これは境界に対して発生しているわけである。

さて、境界の問題が出てきたが、そもそもポアソンの方程式、解こうと思ったら何らかの条件が必要になる。今回は、外周上の格子点の温度を指定することで解く。
このような境界条件を指定しなければいけないのだが、それに合わせて、得られた行列も調整する必要があるのだが、それは次に説明視よう。

ここでは、非常に足早に必要なことも省略しながら説明しているので、この分野の計算についてちゃんと知りたい場合は、「偏微分方程式の数値解法」などのキーワードでググってみよう。あちこちの大学の講義資料や書籍がみつかるので、それらできちんと勉強しよう。

せっかくOctaveについて書いているのに、Octaveが得意な科学技術計算についてちょっとしか書いていない。
ということで、今回は、とつぜん偏微分方程式の数値解に挑戦してみようと思う。

物理現象を偏微分方程式で表すことはよくあるのだが、偏微分方程式をきちんと解くのはほとんど不可能である。
なので、現実には、適当な境界条件のもとで数値解を求めることになる。

とりあえず、結果を先に示そう。
2次元の熱伝導について数値解を求めてみよう。

1辺0.1m(10cm)の板の1つの辺の中央に熱が加わり、周囲の温度が下図のようになる場合について計算してみよう。
温度は、外周で300度(絶対温度、室温)、一辺の中央が1600度で、その一辺の温度分布はsinカーブになっているものとする。
高さ方向は温度になっている。

poisson2d-1.png外周は上図の温度であるが、内部の温度は次の偏微分方程式にしたがって変化する。
右辺が0になっているのは定常状態を示している。
というより、0が計算が楽なのでこれで今回は説明をしてしまう。

poisson2d-3.png
これを、正方形の板をN*Nのメッシュに切って、各点の温度を求めて図示するとこんな感じになる。
高温が赤、低温が青になっているが、これはOctaveで勝手に色づけしてくれるので便利。

poisson2d-2.png
こんな図になるように計算するにはどうすれば良いのだろうか。
この図は、N=50の場合の計算結果である。

実際には、物体は3次元空間に存在するのだから、x,y だけでなく、zも加えて、N*N*Nのメッシュに切って計算するべきなのだが、こうなると点の数が非常に増えて困ることになる。
点の数は、2次元の場合で N*N、3次元の場合に N*N*N になるから大変なのだ。
いや、本当は次元が増えるというのはそんなに甘いものではない。もっと急激にデータ爆発が起きるのだが、そのあたりは次で説明しよう。

RAII の話題は C++ の例が多いが、それもそのはず、ここ最近は GC が面倒を見てくれることがほとんどなので、 RAII の需要そのものが減っている。とはいえネットワークコネクションやらファイルハンドルやらと、あまりリソース解放のタイミングを遅らせたくないものに関しては、 RAII パターンで寿命を明示しておきたいのは GC があったとしても同じこと。

Ruby もそういったリソースに関してはブロックを渡すことで解放もしてくれる RAII パターンを実装したインターフェースが提供されているので、積極的に使っていきたいところだ。そういったリソースがひとつだけなら、素直に書けばいいだけのことなのであまり考える必要はない。

open(...) do |fp|
  something(fp)
end # ここでファイルハンドルが解放される

ところが、配列になっているものに対し、すべての要素を開き、同時に使いたいとなると意外と面倒なことに気付く。配列 arr = [a1, a2, a3, ...] があったとき、意味的にはこういうことがしたい。

open(a1) do |b1|
  open(a2) do |b2|
    open(a3) do |b3|
      ...
      something([b1, b2, b3, ...])
    end
  end
end

が、配列の要素数がわからない以上、この方針でネストを続けていくわけにはいかない。 果て...。どうすればいいだろう?

よくあるモンテカルロ法による円周率を求めるというのをOctaveでやってみる。
図のように、x,y それぞれ-1〜1の範囲の乱数を生成し、一辺2の正方形と単位円を重ねて、単位円の中に入っている点の割合から円周率を求める。
pimontecarlo.png
まずは、ごく普通のプログラミング言語風の書き方で試す。
一応、時間を計測するために、time()を使っている。

function [p, t] = pislow( n )
    t1 = time();
    cnt = 0;
    for i=1:n
        x = rand()*2-1;
        y = rand()*2-1;
        if( x^2+y^2 < 1 )
            ++cnt;
        endif
    endfor
  
    p = cnt / n * 4;
    t = time() - t1;
endfunction

実行すると、こんな感じになる。 精度を上げるために100万回試したら、27秒もかかった。 プログラミング言語としてはかなり遅い方だろう。

>> [p,t] = pislow(1000)
p =  3.1560
t =  0.027001
>> [p,t] = pislow(1000000)
p =  3.1434
t =  27.356

では、Octave的なプログラムにするとどうなるか、見てみよう。
 
function [p,t] = pifast( n )
    t1 = time();
    x = rand(n,1)*2-1;
    y = rand(n,1)*2-1;
    cnt = length(find(x.^2+y.^2<1));
    
    p = cnt / n * 4;
    t = time() - t1;
endfunction
forループが消え、100万個の乱数の配列x,yを作って、配列同士の演算をしている。 findは、条件を満たす場合の添え字の配列を返すので、その個数をlength()で数えることで、条件を満たす点数をカウントしている。

これを実行すると、100万回の場合、340倍ほど高速になっている。 これだけの差が出るので、forループを捨てて、配列演算をどんどん使う重要性が分かるだろう。
つまり、forループなどを使わず、配列演算を利用して短く高速に、バグが入る余地も少なくして横着に高速なプログラムを書くのがOctave風なのだ。

>> [p,t] = pifast(1000)
p =  3.1160
t =  0.0010002
>> [p,t] = pifast(1000000)
p =  3.1421
t =  0.080005

でも、注意すべきことがある。配列演算はとても高速なのだが、配列を確保するので、メモリをどんどん消費する。 要素数100万個の配列で、1要素が8バイトだと、1つの配列で8MB消費する。 1000万だと80MB、1億だと800MB、10億だと8GBにもなってしまい、実際に走らすとシステムが反応しなくなる。

最後に、最初に示した図を描くプログラム(スクリプト)を示そう。
 
hold on;

n = 1000;
x = rand(n,1)*2-1;
y = rand(n,1)*2-1;
lt = find(x.^2+y.^2<1);
ge = find(x.^2+y.^2>=1);

plot(x(lt),y(lt),".r");
plot(x(ge),y(ge),".b");

th = 0:pi/20:2*pi;
plot(cos(th),sin(th),"linewidth",2,"color","red");

hold off
hold on; hold off; があるが、これはFigureの重ね描きを制御するためである。

findで、引数の条件を満たす添え字の並びが得られるが、これを元の配列の引数として与えると、条件を満たす要素(値)だけ抽出することができる。これを利用して、条件を満たす場合と満たさない場合をplotを2回使って色分けして表示している。
Octaveで行列などが扱いやすく数値計算に向いていることは分かったが、結果の数字だけを眺めても様子が分からない。
計算結果を図示することがとても重要だ。図示することで、一瞬でデータの意味が分かる。

とりあえず、x-y 平面の原点からの距離に応じて高さを決めることで、波を描いてみた。

>> a = linspace(-20,20,131);
>> [x,y] = meshgrid(a,a);
>> r = sqrt( x.^2 + y.^2) + eps;
>> z = sin(r) ./ r;
>> surface(x,y,z);
wave.png
aは、-20から20までを130等分した配列である。131になっているのは、両端を含む等分点全部となっているため。
meshgridにより、xにはx座標値が入った配列、yにはy座標値が入った配列が作られる。
rには、原点から座標(x,y)までの距離が入るが、epsを加えるというオマジナイをしている。
そして、波の高さをzとして計算している。
面の表示はsurfaceで、メッシュのx座標配列、y座標配列、z座標配列を引数に与えると、上図が現れる。
マウスを動かすことで、好きな方向から眺められる。

次に、もうちょっと複雑な形を作ってみよう。でも、同じくたった5行でできる形だ。

>> [u,v] = meshgrid(linspace(0,10*pi,200),linspace(0,2*pi,30));
>> x = 2*(1-exp(u/(6*pi))).*cos(u).*cos(v/2).^2;
>> y = 2*(-1+exp(u/(6*pi))).*sin(u).*cos(v/2).^2;
>> z = 1-exp(u/(3*pi))-sin(v)+exp(u/(6*pi)).*sin(v);
>> surface(x,y,z);
snail.png巻貝のような形になっているが、[u,v] の媒介変数で曲面を表しただけである。
uは0から10π、vは0から2πのパラメトリック曲面のx,y,z座標値を求めて表示しているだけなので、処理の原理は同じである。
x,y,zの式をいろいろ弄れば、色々な3次元曲面が作れる訳である。
数値計算処理をするとき、やたらにforループのネストが増えたりする。いかにも計算している気になるかもしれないが、元の数式が複雑になってくると、それをコードに落とすのがどんどん大変になってくる。

そうなると「数式をそのままプログラム中に書きたい」となってくる。

完全にそのままは使える記号の種類とかで無理があるが、数式にできるだけ近い形でプログラムしたくなる。そうすると、プログラムも短くなり、バグの入る余地が減る。それ以上に、プログラムが圧倒的に見やすくなるのだ。

ほぼ数式どおりに書ける数値計算用のプログラミング言語がいつくかあるようだが、前回紹介したOctaveもその1つである。
人工知能で機械学習などをやろうとすると、大量のデータを行列にぶち込んで、行列演算をごちゃごちゃやることが非常に多くなる。
ということで、今回は、Octaveの行列関連の演算を紹介しようと思う。

行列は、[]の中に、行単位で数字を空白またはカンマ(,)で区切りながら並べ、行をセミコロン(;)で区切る。 転置行列は、肩にTの代わりにアポストロフィ(')を付けるだけでOK。

>> A = [1 2 1; 2 1 0; 1 3 2]
A =

   1   2   1
   2   1   0
   1   3   2

>> B = A'
B =

   1   2   1
   2   1   3
   1   0   2

以下にAの逆行列を3つの方法で計算してみた。行列の-1乗がちゃんと逆行列になってくれる。
eye(3)は、サイズ3の単位行列で、これを行列Aで割っても逆行列が求まる。
 
>> inv(A)
ans =

  -2.00000   1.00000   1.00000
   4.00000  -1.00000  -2.00000
  -5.00000   1.00000   3.00000

>> A ^ -1
ans =

  -2.00000   1.00000   1.00000
   4.00000  -1.00000  -2.00000
  -5.00000   1.00000   3.00000

>> eye(3) / A
ans =

  -2.00000   1.00000   1.00000
   4.00000  -1.00000  -2.00000
  -5.00000   1.00000   3.00000

本当に逆行列になっているか、行列Aとの積を計算してみよう。
 
>> (A ^ -1) * A
ans =

   1.00000   0.00000   0.00000
   0.00000   1.00000   0.00000
   0.00000   0.00000   1.00000

いきなり逆行列を求めてしまったが、ランク、行列式を求めてみよう。 さらに、固有ベクトル、固有値を求めてみる。
 
>> rank(A)
ans =  3
>> det(A)
ans = -1.0000
>> [V,L] =eig(A)
V =

   0.50393   0.47733   0.10741
   0.34334  -0.64817  -0.46439
   0.79257   0.59332   0.87909

L =

Diagonal Matrix

   3.93543         0         0
         0  -0.47283         0
         0         0   0.53740

こんな感じで、行列演算が簡単にできる。
ここでは、説明のために、1つ1つ示したが、これらの演算を連続してやっている数式は、Octaveなら1つの数式として書ける筈だ。

Octaveは、数値計算が便利にできるだけではなく、計算結果を3次元表示することもできる便利な道具である。なので、今後数回に渡って、色々紹介しようと思っている。
人工知能プログラミングでは、大量の数値解析が必要になるので、それをいちいちプログラムしていたら面倒でしかたがないので、何らかのソフトを利用して楽をするわけだが、その中の有力なプログラミング言語としてMATLABがある。でも、これは有料である。
やはり気楽に使うためには、オープンソースであって欲しい。実はMATLABと互換性のある GNU Octave があるので、これを使ってみよう。

gnu-octave-screenshot.pngインストールは適当にやれば簡単に終わるのだが、プログラミング言語なので、それなりに勉強しないと分からない。
行列の扱いが簡単で、転置行列、行列式、逆行列などすぐにできる。
それだけではなく、3次元表示も可能である。

とりあえず簡単な入門をネットで探していたら、WIKIBOOKSの中に、Octave Programming Tutorialがあった。
どうしても日本語という人には、東海大学稲葉研究室の Octave入門 が良いかもしれない。
日本語の本で学習したい場合は、ネットで探してみよう。何冊も出ているようであるが手元に無いので評価はしていない。

しかし、今の時代、ネットを使ったビデオ学習にしよう、その方が楽ではないか、正道ではないか。
YouTubeで探したら、Octave Tutorials として、全35本のビデオでまとまっている。
さらに、これは、カリフォルニア大学アーバイン校で行われたan introduction to computer programming courseである。講師はPaul Nissensonである。
ということが分かったので、これは見るしかあるまい。無料でUCIのコースを受講できるチャンスだ。

内容的にはTutorialなので、難しいことは何も無かった。
Tutorialなので、あまり詳しいところまでの説明はしていない。
でも、一番困ったのは、プログラムを実行している画面の文字が画面の解像度に比べて小さくて、文字の記号が判別できない。まあ、音声が分かれば何を遣っているかは全部わかるのだが、そこだけが問題だった。

でも、よく考えてみると、これはものすごく良い英語のリスニングの勉強になった。
今までは英語をいい加減に聞き流しながら、画面だけを見ていたのだが、このビデオではまったく逆になってしまったのだ。まったくの想定外だが、想定外の学習効果が期待できる一石二鳥の教材だ。


人工知能ブームで、機械学習やら深層学習の説明がネット上にはゴロゴロしている。
しかし、丁寧に、大学生はもちろん、高校生くらいなら十分に分かるように説明しているものは少ない。

data-school.pngそれで、日本語だけで探すのはあきらめて、英語で探すと色々見つかったが、その中で一番丁寧な説明をしていると思われるものを今回は紹介する。

data school というサイトには、データサイエンス関連の色々な情報が用意されている。

実際に見てみたのは、scikit-learn だけなのだが、たぶん他も丁寧hqdefault.jpgなのではないかと予想している。
scikit-learn は、pythonで機械学習をするとき非常に良く使う便利なライブラリなのだが、このビデオでは、機械学習のことも説明しながら、ライブラリの説明もしてくれるというものでるので、事前の準備は何もいらない。

scikit-learn ビデオは9本で、合計4時間もある。さらに、textからの機械学習の2時間40分の長い長いビデオが付録でついている。

それぞれのビデオには、対応した notebook がGitHub上に用意されている。
noteは、単なるメモではなく、ビデオで実演してくれる内容、プログラムが全て載っており、自分の
IPython上にnotebookのプログラムをコピペしながら実行することで簡単に確認(体感)できる。
説明は、これ以上ないほど細かいステップで行われるので、理解はとても容易だった。

説明がていねい過ぎて、考えを飛躍しないといけない箇所がないのが不満になるかも知れない。
これだけていねいに、ビデオおよび資料を用意するのは恐ろしく大変だと思う。LearningPython.gif

日本でも、色々な分野のビデオや解説書がこのくらいのていねいさで作られれば、誰でも容易に人工知能、機械学習を習得できると思うが、残念ながら日本では短いビデオ、薄い本であることが重視され、困った状態だ。
専門書は、1000ページあるのが普通という世界にならないものだろうか。
たとえば、右の"Learning Python"というPythonの入門書は、1600ページもある。
日本語版は版が古く、800ページしかない。もう英語版で読むしかないだろう。

なお、このビデオ教材の付録のtext処理のビデオは2時間40分ととても長いので、まだ見ていない。


人工知能系のプログラムは、実はとても計算量が多いので、プログラムが効率がよいかどうか常に気にしておく必要がある。

時間を計測するのに、自分で、対象範囲の前後でシステムの時間を読み取り、差を取るということで求めることはできる。でも、そんなことをしなくても、Pythonには時間計測の仕組みが用意されている。

以下では、リストに追加するときに、insertとappendでどのくらい差が出るかを計測してみたものだ。計測する部分は、それぞれ関数にまとめてある。

import timeit

def insfirstfunc(n):
    a = []
    for i in range(n):
        a.insert(0,i)

def inslastfunc(n):
    a = []
    for i in range(n):
        a.insert(len(a),i)

def apdfunc(n):
    a = []
    for i in range(n):
        a.append(i)

for i in range(1000,10001,1000):
    t = timeit.Timer("insfirstfunc(i)","from __main__ import insfirstfunc,i")
    t1 = t.timeit(number=1000)
    t = timeit.Timer("inslastfunc(i)","from __main__ import inslastfunc,i")
    t2 = t.timeit(number=1000)
    t = timeit.Timer("apdfunc(i)","from __main__ import apdfunc,i")
    t3 = t.timeit(number=1000)
    print("%6d  t1=%10.6f  t2=%10.6f  t3=%10.6f" % (i,t1,t2,t3))

timeitをインポートしておいて、timeit.Timerでタイマーを生成している。
Timerの第1引数に実行したいことを文字列で書く。ここには、変数などが含まれてもよい。
第2引数に、"from  __main__  import"の後ろにインポートしたいものを並べておけば大丈夫なようだ。

Timerを作ったら、実際の時間計測を作ったTimerオブジェクトのtimeitメソッドで計測できる。引数に number=繰り返し回数 を指定することができる。とくに短時間の計測の場合、システムの状態に大きく依存してしまうので、なんども繰り返し実行すると誤差が減る。

詳細はマニュアルを参照して欲しい。

iPython3 で実行してみた結果を以下に示す。

In [1]: %run timeit
  1000  t1=  0.395643  t2=  0.204404  t3=  0.073373
  2000  t1=  1.217943  t2=  0.423568  t3=  0.153751
  3000  t1=  2.477653  t2=  0.639693  t3=  0.238761
  4000  t1=  4.179113  t2=  0.859387  t3=  0.316046
  5000  t1=  6.390812  t2=  1.081363  t3=  0.400768
  6000  t1=  9.033100  t2=  1.297038  t3=  0.481604
  7000  t1= 12.073041  t2=  1.506441  t3=  0.561811
  8000  t1= 15.538109  t2=  1.731817  t3=  0.645874
  9000  t1= 20.012774  t2=  2.071590  t3=  0.771902
 10000  t1= 24.719389  t2=  2.313333  t3=  0.866832

In [2]: 

appendが一番高速で、insertで最後に挿入するのは3倍程度遅くなる。 そして、insertで最初に挿入するのは、O(n**2)のオーダーで遅くなっているのが分かる。
appendで済むなら、insertを使うのはやめるべきなのだ。

timeitを使うことで、big-O がどの程度か簡単に調べることができる。

最近のコメント