2017年7月16日 (日)

Juliaやってみよう。四日目。@timeでプロファイリング

今日は、昨日のIzhikevichモデルのコードの簡単なプロファイリングしてみる。Pythonに翻訳した時とどっちが速いのかというのが知りたい。結果次第ではjuliaは使わないかもしれない。ちょっと古い記事だけど、ある程度アレイが大きくなるとNumpyの方が速いというベンチマークもあるので、自分が使うような用途で確認してみたい。

@time

公式ドキュメンテーションのパフォーマンスのコツによると、基本的には@timeというマクロがあるので、それを使うと良いようだけど、グローバルのスコープで実行するとオーバーヘッドがあって遅くなる(からのForループでも100マイクロ秒くらい)ので、測りたい部分をfunctionにして測るのがよいらしい。

昨日のコードをfunctionの中にいれてやって、

function test()

    Ne = 800; Ni = 200;
    re = rand(Ne, 1); ri = rand(Ni, 1);
    a = [0.02 * ones(Ne, 1); 0.02 * ones(Ni, 1)]
    b = [0.2 * ones(Ne, 1);  0.25-0.05*ri]
    c = [-65 + 15 * re .^ 2; -65*ones(Ni, 1)]
    d = [8-6*re .^ 2; 2*ones(Ni,1)]
    S = hcat(0.5 * rand(Ne+Ni,Ne),  -rand(Ne+Ni,Ni));

    v=-65*ones(Ne+Ni,1);    # Initial values of v
    u=b.*v;                 # Initial values of u
    firings=[];             # spike timings

    for t=1:1000            # simulation of 1000 ms
        I=[5*randn(Ne,1);2*randn(Ni,1)]; # thalamic input
        fired=find(v.>= 30);    # indices of spikes
        firings=[firings; hcat(t+0*fired, fired)];
        v[fired]=c[fired];
        u[fired]=u[fired]+d[fired];
        I=I+sum(S[:,fired], 2);
        v=v+0.5*(0.04*v.^2+5*v+140-u+I); # step 0.5 ms
        v=v+0.5*(0.04*v.^2+5*v+140-u+I); # for numerical
        u=u+a.*(b.*v-u);                 # stability
    end

end

あとは

@time test()

としてやる。

うちのPCでは

0.785356 seconds (326.64 k allocations: 365.923 MiB, 24.83% gc time)

という感じでした。

numpyバージョンを作らねば。

2017年7月15日 (土)

Juliaやってみよう。三日目。MATLABコードを翻訳してみる。

まあ、多少慣れてきたのでMATLABのコードを翻訳してみる。

お題は、Izhikevichモデルの有名な1000個のニューロンのシミュレーションするモデル。短いし、まあ、こーゆーことをやるための言語だし。

MATLABのオリジナルはこちら。

% Created by Eugene M. Izhikevich, February 25, 2003
% Excitatory neurons    Inhibitory neurons
Ne=800;                 Ni=200;
re=rand(Ne,1);          ri=rand(Ni,1);
a=[0.02*ones(Ne,1);     0.02+0.08*ri];
b=[0.2*ones(Ne,1);      0.25-0.05*ri];
c=[-65+15*re.^2;        -65*ones(Ni,1)];
d=[8-6*re.^2;           2*ones(Ni,1)];
S=[0.5*rand(Ne+Ni,Ne),  -rand(Ne+Ni,Ni)];

v=-65*ones(Ne+Ni,1);    % Initial values of v
u=b.*v;                 % Initial values of u
firings=[];             % spike timings

for t=1:1000            % simulation of 1000 ms
  I=[5*randn(Ne,1);2*randn(Ni,1)]; % thalamic input
  fired=find(v>=30);    % indices of spikes
  firings=[firings; t+0*fired,fired];
  v(fired)=c(fired);
  u(fired)=u(fired)+d(fired);
  I=I+sum(S(:,fired),2);
  v=v+0.5*(0.04*v.^2+5*v+140-u+I); % step 0.5 ms
  v=v+0.5*(0.04*v.^2+5*v+140-u+I); % for numerical
  u=u+a.*(b.*v-u);                 % stability
end;
plot(firings(:,1),firings(:,2),'.');

一晩格闘した結果、小生なりのJuliaへの翻訳がこちら。 オリジナルにはないVmを追加している。

using PyPlot

Ne = 800; Ni = 200;
re = rand(Ne, 1); ri = rand(Ni, 1);
a = [0.02 * ones(Ne, 1); 0.02 * ones(Ni, 1)]
b = [0.2 * ones(Ne, 1);  0.25-0.05*ri]
c = [-65 + 15 * re .^ 2; -65*ones(Ni, 1)]
d = [8-6*re .^ 2; 2*ones(Ni,1)]
S = hcat(0.5 * rand(Ne+Ni,Ne),  -rand(Ne+Ni,Ni));

v=-65*ones(Ne+Ni,1);    # Initial values of v
u=b.*v;                 # Initial values of u
firings=[];             # spike timings

Vm = [];                # added

for t=1:1000            # simulation of 1000 ms
  I=[5*randn(Ne,1);2*randn(Ni,1)]; # thalamic input
  fired=find(v.>= 30);    # indices of spikes
  firings=[firings; hcat(t+0*fired, fired)];
  v[fired]=c[fired];
  u[fired]=u[fired]+d[fired];
  I=I+sum(S[:,fired], 2);
  v=v+0.5*(0.04*v.^2+5*v+140-u+I); # step 0.5 ms
  v=v+0.5*(0.04*v.^2+5*v+140-u+I); # for numerical
  u=u+a.*(b.*v-u);                 # stability
  Vm = [Vm; v'];
end;

scatter(firings[:,1], firings[:,2], s=1);

いやー、なんかほとんどおんなじだ。JuliaはMATLABと同じでインデックスが1から始まるし、MATLABキラーだな。forループがendで終わってPythonよりMATLAB感がつよい。

PythonでNumpy使い始めたら、MATLABはv(2:5)とかベクターをスライスするときも括弧使うのが、ちょっと変な仕様だと気づいたけど、JuliaもNumpyのようにv[2:5]とスクエアブラケットでのインデックシングでしっくり来る。納得の変更・改善ばかり。後出しジャンケンは強い。

import numpyってしてないのに、素のJuliaでrandとか、onesとか使えるのがすごい。こういうアレイとかも実装はjuliaでJITまでついてくるんだからすごいな。pypyプロジェクトがnumpypyのRPythonでの実装に数年かかっているのをみると、juliaのアプローチの優秀さが分かる。numpypyは99%実装できたみたいだけど、まだ遅いので使い道ない・・。 numpyでいうvstackとかhstackはvcat, hcatらしい。

GRではscatterがよくわからなかったのでPyPlotでプロットした出力がこちら。

Izhikevich_plot

上のコードではVmはプロットしてないけど、するならこんな感じ。。

plot(Vm[:,32])
ylim([-90,100]);

ylimとかnumpyそのまま。

«Juliaやってみよう。二日目。Juliaのアップデート、IJulia

広告欄


やっつけタイムライン

広告欄

はてブ

人目の訪問です。

  • follow us in feedly

    かなり更新が不定期なため、RSSリーダーをオススメします。現在Feedlyに122人登録頂いています。多謝!RSSを表示

    ブログランキング用 にほんブログ村 IT技術ブログ Pythonへ ブログランキングならblogram






    Jenny Mayhem
2017年7月
            1
2 3 4 5 6 7 8
9 10 11 12 13 14 15
16 17 18 19 20 21 22
23 24 25 26 27 28 29
30 31          

IT技術注目記事

無料ブログはココログ