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年のブログ記事の信憑性がましてしまった。

続きを読む "Juliaやってみよう。五日目。Pythonと速度比較。" »

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やってみよう。四日目。@timeでプロファイリング" »

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やってみよう。三日目。MATLABコードを翻訳してみる。" »

2017年7月13日 (木)

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

さて、昨日に引続き、今日もちょこちょこやっていきます。

今日は、Juliaをv0.5からv0.6にアップデートして、それからJupyter notebookのJulia用であるIJuliaで遊んでみます。

Juliaをアップデート

ウィンドウズなので、Ubuntuのようにapt-getする訳にはいきません。

sudo apt-get dist-upgrade julia

スタートボタンからJuliaを見に行くとUninstallerがありました。たぶん旧バージョンを残したまま、新しいのをインストールすることもできるっぽいですが、べつにまだ始めたばかりで旧バージョンでしか動かないようなコードもないので、今回はまっさらにアンインストールして、最新版をインストールでいいか。

C:\Users\ユーザー名\AppData\Local\Julia-0.6.0

インストール先にバージョン名が入っているので、旧バージョンと混在できるのでしょう。あとはウィンドウズのPATHの設定で切り替えるのだと想像。Pythonもそんな感じ。Condaならコマンドで切り替え。

IJuliaのインストール

公式サイトにあったIJuliaのチートシートをみながら使ってみる。

using IJulia

するとPkg.add("IJulia")しろと怒られる。Pkg.status()してみると、パッケージなにもはいっていない。昨日のv0.5は多分以前ちょっとだけIJuliaを起動してみた残骸だったと思い当たる・・。

おとなしく

Pkg.add("IJulia")

してみる。依存性とか勝手に入ると想像した通り、色々とインストールされた。

JupyterがPythonに依存するあたりで、condaがJupyterを最新版に保っていたりとしている感じに見える。結構かかった。

あとは昨日のplotライブラリのGRも追加してみる。

Pkg.add("GR")

うむ。問題なく終了。

IJuliaの起動

using IJulia

したいところだけど、Jupyterって起動したときのcmdのパスから起動するので、まずは適当なパスにcdしたい。

公式サイトのドキュメンテーションみるとFile Systemの項目にcdというコマンドがある。

cd("C:\\projects\\julia\\test1")

なんて感じで大丈夫だった。バックスラッシュは二重にしてエスケープするかフォーワードスラッシュにしないとだめ。あと、cdの後にスペースいれて

cd ( "path" )

なんてやるとエラー。なるほど。あとダブルクォーテーションじゃないと怒られる。PythonならシングルでもOKだし、括弧の前にスペースあっても大丈夫だと思う。この辺のルールは厳密な方が、みんなの書くコードが統一されて読みやすくなっていいと思う。Pythonで多くの人々がこの恩恵を思い知ったのでJuliaにも引き継がれているのかな。

pwd()

とすると現在のワーキングディレクトリが表示される。

using IJulia

してから、

IJulia.notebook()

もしくは単に

notebook()

でJupyter notebookのJulia用、IJuliaが起動。

おお、homedir()で出てくるC:\Users\ユーザー名で起動したぞ。

Pythonだと起動した時のパスでJupyter立ち上がった気がするが、IPython notebookのときの記憶だし、最近は違うのだろうか。

IJuliaのnotebookサーバーの終了の仕方もよくわからない。PythonならcmdにCtrl+C送ればよかったが、juliaのコンソールだとだめっぽい。しょうが無いのでコンソールごと終了。

notebookコマンドにパスを渡して起動してみたい。プロンプトにはてなマークを入力するとプロンプトが

julia>

から

help?>

に切り替わって、そこでnotebookとするとドキュメンテーションがでてきた。なるほどdir=で指定するのね。

Notebook_help

ちなみにエクスプローラーのアドレスをマウスでドラッグして、Juliaのコンソールに落とすとドロップされる。

スタートメニューに追加されているJupyter notebookのショートカットからだとCtrl+C効くけどJuliaのコンソールは相変わらず受付ない。でもHelpにコンソールを終了すれば良いとも書いてあるので、まいっか。

さて、JupyterのHomeの右上にあるnewボタンからJulia v0.6を選択して、ノートブックを作成。Untitledをクリックしてタイトル変更。

Ijulia_gr

適当にプロットしてみた。inlineやんないとプロットが別ウィンドウに出てきた。セミコロンで出力抑制とかPythonだな。つーか、よく考えるとlinspaceを0から始めて100要素なら99で終わるべきだった・・・。

要素ごとの階乗でスペースしてドットしてハットするのは数学っぽくっていいね。numpyだと x ** 2 とやるので、xがベクターだと知らないとスカラーに見えてしまう。

まあなんとなくわかってきたぞ。今日はここまで。

続きを読む "Juliaやってみよう。二日目。Juliaのアップデート、IJulia" »

2017年7月12日 (水)

Juliaやってみよう。一日目。GRでプロット。

色々と気になっていたJuliaで遊んでみようのコーナーを思いつきで始めます。

普段はウィンドウズでPython使っているのでウィンドウズ、Pythonユーザー向けに書きます。

まずはなんかPlotしてみたいかな。

まずはJuliaのインストール

右も左も分からないので、まあOfficialからバイナリもってきて遊んでみよう。

https://julialang.org/downloads/

から最新版のウィンドウズ用exeインストーラーを落とせばよいのでしょう。今時は普通は64ビットですね。 執筆時点ではもう0.6リリースなってますね。私は手元に0.5があったので、今日はそのまま使います。

ライブラリのインストール

さて、JuliaにもPlotライブラリが色々とあるようです。どれにしようかな。

plot.lyのJulia版もあるようですが、JSON書くのいやなので、パス。

公式サイトのhttps://julialang.org/downloads/plotting.html

では、PyPlotというそのままな名前のライブラリでmatplotlibをJuliaから呼び出しちゃうという恐ろしいものがあった。うーむ、JupyterがIPython Notebookの時代にそんなデモをみたような。これはPythonを呼びに行っているのでPythonも必須。

せっかくなのでJuliaネイティブなライブラリを試したいので全部C/C+で書いてあるというGRというPlotライブラリを試すことにした。これPythonはもちろん、PyPyでも使えるらしいうえに、matplotlibのバックエンドとしてもつかえて30倍くらい線の描画が速いらしい。すごい。。。

JuliaにはPkgという多分標準ライブラリだと思うけど、パッケージマネージャーがあって、登録されているライブラリはPkg.add()で簡単に導入できるらしい。すごい。PythonもCondaで似たようなことができるけど、Juliaの方が進んでいる。

まずはPkg.status()で、導入されているライブラリを表示してみる。

Juliapkgstatus

ふむ。

それでは、Pkg.add("GR")でGRを導入。

Juliapkgadd

INFO: METADATA is out-of-date — you may not have the latest version of GR INFO: Use `Pkg.update()` to get the latest versions of your packages

へ?METADATAが古い?なんの?Pkg.update()すればいいの?

Juliapkgupdate

おお、PackageマネージャーのMETADATAが古いのですな。最新のGRは0.22.0になっているようです。

なんか裏で色々とcondaが暗躍している模様だが、Pkgの実装はピュアJuliaではないのだろうか。なぞ。

さて、Pkg.installed()としてみると、Pkg.update()がすでに全部アップデートしたようで、GRもバージョンが0.22.0になっている。楽じゃ!condaみたい。

GRで初めてのプロット

Pythonでいうimport分はusingらしいが、importというのもある。違いは調べてもよくわからないが、まあ追々でいいか。GRの公式サイトでもusingだったりimportだったりで、まあどっちでもいいようだ。

あとはPythonのようにplot([0,1,4],[3,2,5])のようにするだけでいいのだそうです。

using GR
plot([0,1,4],[3,2,5])

えい。

でた。

これって、plotはGR.plotと同じだと思うのでusingでもってくるとネームスペースのトップにエイリアスがくるっぽい。

タイトルを追加してみる。今度はGRから始まるネームスペースでやってみる。

GR.title("Yattsuke blog playing with Julia")

シーン・・・。何も起きない。PyPlotならshow()とかdraw()の場面だが、show()だとプロットが死んだ。drawはない。

GRの公式サイトの例だと

GR.updatews()

が正解っぽいが、何も起きない。さて・・・。

plotもう一回すると更新されるが、ドキュメンテーションみてもやっぱりupdatewsで正解っぽいなあ。。

うむ。まあ今回はこの程度で、つぎはIJulia使ってみるメウ。

Juliagrplot

続きを読む "Juliaやってみよう。一日目。GRでプロット。" »

しろののツイッタータイムライン

  • ツイッターは5つ目も凍結されました。6つ目での復活も不可能。なのでnoteに注力しています。

    と思ったら、イーロン・マスクの買収が公になってアカウントが復活できました。ありがとうマスク。

    トランプ関連記事の一覧リスト

オススメたち

2024年12月
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        

はてブ

無料ブログはココログ