2017年8月 1日 (火)

Juliaやってみよう。五日目。Pythonと速度比較。

def test():

    Ne=800;                 Ni=200;
    re=rand(Ne,1);          ri=rand(Ni,1);
    a=vstack([0.02*ones((Ne,1)), 0.02+0.08*ri]);
    b=vstack([0.2*ones((Ne,1)),  0.25-0.05*ri]);
    c=vstack([-65+15*re**2,  -65*ones((Ni,1))]);
    d=vstack([8-6*re**2,       2*ones((Ni,1))]);
    S=hstack([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 in range(1000):            # simulation of 1000 ms

        I=vstack([5*randn(Ne,1), 2*randn(Ni,1)]); # thalamic input
        fired=find(v>=30);    # indices of spikes
        if fired.any():

            if firings == []:
                firings=vstack([t+0*fired,fired]).T
            else:
                firings=vstack([firings, vstack([t+0*fired,fired]).T]);

            v[fired]=c[fired];
            u[fired]=u[fired]+d[fired];
            I=I+sum(S[:,fired],1).reshape(-1,1);

        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

全快の記事では、Izhikevichモデルを1000 msほどシミュレーションした時のjuliaの実行時間は785 msでした。

速度比較のため、上記のようなnumpyバージョンを作ってみたところ、253 msでした。

Python版ではfiredを調べて空だったらシナプス入力の計算を端折るようにしたので、そのせいで速い可能性もあるので、juliaも同様に

if length(fired) > 0

という風にチェックを入れたら625 msまで速くなった。それでもPythonが三倍くらい速いとかあり得ない。なにかおかしい。

もしかしたらjuliaのJITが温まる前なのかな。ループを長くしたら逆転するかも。PypyのJITが温まるまで4秒くらいとか聞いた気がするので、シミュレーションを1000 msから30000 msにしてみる。

Python2   23.4 s
julia    34.48 s

結果、差はだいぶ縮まったが、それでもPythonがちょっと速い。これならJITのあるjuliaが有利な条件だと思うけども。ふーーむ。juliaの偉い人が颯爽とアドバイスくれたりないかなぁ・・。

まあKyle Barbary氏の2014年のブログ記事の信憑性がましてしまった。

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バージョンを作らねば。

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

広告欄


やっつけタイムライン

広告欄

はてブ

人目の訪問です。

  • follow us in feedly

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

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






    Jenny Mayhem
2017年8月
    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技術注目記事

無料ブログはココログ