カテゴリー「Julia」の5件の記事

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

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そのまま。

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がベクターだと知らないとスカラーに見えてしまう。 まあなんとなくわかってきたぞ。今日はここまで。

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

広告欄


広告欄

Amazon

  • ブログ記事にした商品のアフィ

    デジピはヤマハが音質・キータッチとも一番良いです。

    レビュー記事かいた程おすすめのワイヤレススピーカー。

    安いタブレットとして買い替えに最適

    速くKindle化して欲しい

やっつけタイムライン

はてブ

人目の訪問です。

  • follow us in feedly

    かなり更新が不定期なため、RSSリーダーをオススメします。RSSを表示

    当ブログは超安定なLinux Mintとシンプル&スタイリッシュなElementaryOSを応援しています。







    Jenny Mayhem
2019年1月
    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    
無料ブログはココログ