ジュリアに乾杯

『筆跡事件ファイル』なんていう本を読んだ。現役の筆跡鑑定人が書いた本。

筆跡心理学っていう、書かれた文字から人物像を推定して行くって方法も併用してるとか。 これは正にプロファイリング。

一番多い依頼は、遺言書の鑑定とか。後は各種の契約書の鑑定。果たして本人が書いた ものであるかの証明ですな。他人が都合よくでっち上げたものが通用しちゃうんじゃ 困りますからね。

大体こういうのは民事裁判で争われるんだけど、民事と言うのは真実の探求じゃなくて 、私的紛争を効率よく解決するのが目的。よって、真実は自分で証明する必要がある。 かくして、筆跡鑑定が出てくる訳だ。判決を決定する裁判官にいかにアピールするかが 腕の見せ所らしい。

初期鑑定で3万円ぐらい。裁判所に提出するやつになると、50万から80万円ぐらいは 行くとか。

筆跡鑑定は、基準となる本人が書いた文字と鑑定に供する文字の比較。これで真贋を判定 する事になる。普通は同じ文字 数種を対象に比較する事になるんだけど、中には悪徳鑑定人がいて、 比較する文字を1文字に限定しちゃうとかもあるみたい。

文字も癖とかを勘案してくとの事で、おいらみたいにそんなのopenCVで一発判定できんじゃ ねぇってのは、真っ先に否定してたよ。なんだか、奥が深いんですなあ。

場合によっては、文字だけじゃなくて、チェックマークみたいなものでも俎上に上るとか。 某スーパーの納品書の物品欄にチェックマークを余計に付けて、架空納品した事にして、 何億も稼いだ人がいたとか。手口は、現場で店の人と納品チェック(マーキング)、それを 悪徳営業マンが事務所へ持って行って事務処理してもらう。事務所では、現場の人のサインが 有るんで、そのまま信用。

悪い人は、途中でチェックマークを余計に付けていた。但し、緊張感からか、マークの筆圧と ハネ具合が、現場でのものと異なった。鑑定人は見事にこれを見破った。心理学もからめて いるんでしょうな。

その多、勝手に離婚届けを出すとか、裏の怪しい話がいっぱい出てきましたよ。もしもの ために、 筆跡鑑定 とか 筆跡鑑定料金 が、きっと役に立ちます。

Julia との出会い

日経Linuxのmatzさんの連載が終わってしまうので、記念に読んでみた。R,APL,Jって事で rubyの弱い部分を補完してくださいって事だな。アレーに重きを置いた言語。統計に うってつけ。

そして最後に紹介されたのが、Juliaだ。この言語名前だけは 知ってたけど、 Julia : スクリプト言語最速? 手軽さと速さを求めた科学技術計算向け言語 によると、立ち上がりが遅い == Java みたいなものだと思って、食指が動かなかったんだ。

でも、興味が沸いてきて、ぐぐる先生に最近の出来はどうですかって聞いてみたら、0.2までは 確かに遅かったけど、0.3系で早くなってるよーとの事。ならば入れてみるか。

Windowsとウブには出来合いのものが用意されてるけど、両方共おいらの眼中に無し。 だったら、LLVM環境が整っている万次郎Linuxでやってみる鹿。

Julia言語のインストール

First, acquire the source code by cloning the git repository:

git clone git://github.com/JuliaLang/julia.git

Next, enter the julia/ directory and run make to build the julia executable. 
To perform a parallel build, use make -j N and supply the maximum number of 
concurrent processes. When compiled the first time, it will automatically 
download and build its external dependencies. This takes a while, but only 
has to be done once. If the defaults in the build do not work for you, and 
you need to set specific make parameters, you can save them in Make.user. 
The build will automatically check for the existence of Make.user and use 
it if it exists. Building Julia requires 1.5GiB of disk space and approximately
700MiB of virtual memory.

ソースはgitから落とすのは今風でいいんだけど、関連品もコンパイル中に自動的に 落としてきてくれるんか。もう、Netに繋がっているのが前提って事だな。そして コンパイルを加速するのにパラレルコンパイルしろとな。ええ、どうせおいらはちんけな 石しか載せてませんよ。Disk領域が1.5Gと仮想記憶もたっぷり要るとな。

で、コンパイルを始めたはいいけど、途中でgfortranがいるとか言われたよ。そんな事 いの一番に言って欲しいぞ。(HPの下の方に必要品が書いてあった。流し読み重要)

しょうがないので、gcc-fortran 4.8.2-7 を入れてから、コンパイルを継続。ガリガリと Diskの寿命を短くする音が7時間ぐらい続いて、最後は

[sakae@manjaro julia]$ make
osutils.jl
/bin/sh: 1 行:   708 Illegal instruction     (コアダンプ) /home/sakae/src/julia/usr/bin/julia-readline --build /home/sakae/src/julia/usr/lib/julia/sys0 sysimg.jl
Makefile:81: recipe for target '/home/sakae/src/julia/usr/lib/julia/sys0.bc' failed
make[1]: *** [/home/sakae/src/julia/usr/lib/julia/sys0.bc] Error 132
Makefile:32: recipe for target 'release' failed
make: *** [release] Error 2

こんなエラーで止まったよ。これはもう万事休すかな。ってんで、ソースツリーをうろうろ してたら、Make.incなんてのを発見。こいつを開いてみると、どうもconfigの代わりぽい 事に気付いた。

# Use libraries available on the system instead of building them
USE_SYSTEM_LLVM=1
USE_SYSTEM_LIBUNWIND=0
USE_SYSTEM_READLINE=1
USE_SYSTEM_PCRE=1
USE_SYSTEM_LIBM=0
USE_SYSTEM_OPENLIBM=0
UNTRUSTED_SYSTEM_LIBM=0
USE_SYSTEM_OPENSPECFUN=0
USE_SYSTEM_BLAS=0
USE_SYSTEM_LAPACK=0
USE_SYSTEM_FFTW=1
USE_SYSTEM_GMP=1
USE_SYSTEM_MPFR=1
USE_SYSTEM_ARPACK=0
USE_SYSTEM_SUITESPARSE=0
USE_SYSTEM_ZLIB=1

最初、ここが全て0になってたんだけど、物は試しと万次郎側の物も使ってあげる事に した。そして、depsの中にあるライブラリィーで一番容量を喰っていたopenblas-v0.2.8 だけを再利用する事にした。(ソースツリーは、日を置いて取り直した) これで、コンパイルしたら、無事に出来たみたい。

ソースツリーの下にusrってのが出来て、その中に一式入っていたよ。usr/binの下には julia-realineってのとjulia-basicってのがあるけど、basicの方は、miniruby相当なのかなあ? 大分容量が違うぞ。普通はreadlineの方を使えばいいのかな。リンクが貼られていたぞ。

[sakae@manjaro julia]$ ldd ./julia
        linux-gate.so.1 (0xb7743000)
        libreadline.so.6 => /usr/lib/libreadline.so.6 (0xb76e4000)
        libncursesw.so.5 => /usr/lib/libncursesw.so.5 (0xb7685000)
        libjulia.so => /home/sakae/julia/usr/bin/../lib/libjulia.so (0xb6744000)
        libz.so.1 => /usr/lib/libz.so.1 (0xb672d000)
        libpthread.so.0 => /usr/lib/libpthread.so.0 (0xb6710000)
        libffi.so.6 => /usr/lib/libffi.so.6 (0xb6709000)
        libdl.so.2 => /usr/lib/libdl.so.2 (0xb6704000)
        librt.so.1 => /usr/lib/librt.so.1 (0xb66fb000)
        libstdc++.so.6 => /usr/lib/libstdc++.so.6 (0xb6611000)
        libm.so.6 => /usr/lib/libm.so.6 (0xb65cb000)
        libgcc_s.so.1 => /usr/lib/libgcc_s.so.1 (0xb65ae000)
        libc.so.6 => /usr/lib/libc.so.6 (0xb63ff000)
        /lib/ld-linux.so.2 (0xb7744000)
[sakae@manjaro julia]$ file julia
julia: symbolic link to `/home/sakae/julia/usr/bin/julia-readline'

gfortranでコンパイルされたライブラリィー類が出てきてないけど、libjulia.soの中にでも 吸収されちゃってるの?

[sakae@manjaro julia]$ ./julia
               _
   _       _ _(_)_     |  A fresh approach to technical computing
  (_)     | (_) (_)    |  Documentation: http://docs.julialang.org
   _ _   _| |_  __ _   |  Type "help()" to list help topics
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 0.3.0-prerelease+1757 (2014-02-26 04:04 UTC)
 _/ |\__'_|_|_|\__'_|  |  Commit 8854bad* (0 days old master)
|__/                   |  i686-pc-linux-gnu

julia> help( "Numbers")
Help is available for the following items:
Base.one Base.zero Base.pi Base.im Base.e Base.catalan Base.Inf Base.Inf32 
Base.Inf16 Base.NaN Base.NaN32 Base.NaN16 Base.issubnormal Base.isfinite 
Base.isinf Base.isnan Base.inf Base.nan Base.nextfloat Base.prevfloat 
Base.isinteger Base.isreal Base.BigInt Base.BigFloat Base.get_rounding 
Base.set_rounding Base.with_rounding Base.count_ones Base.count_zeros 
Base.leading_zeros Base.leading_ones Base.trailing_zeros Base.trailing_ones 
Base.isprime Base.primes Base.isodd Base.iseven

julia> help("Statistics")
Loading help data...
Help is available for the following items:
Base.mean Base.std Base.stdm Base.var Base.varm Base.median Base.median! 
Base.hist Base.hist Base.hist2d Base.histrange Base.midpoints Base.quantile 
Base.quantile Base.quantile! Base.cov Base.cor

julia> help("Base.mean")
Base.mean(v[, region])

   Compute the mean of whole array "v", or optionally along the
   dimensions in "region". Note: Julia does not ignore "NaN"
   values in the computation. For applications requiring the handling
   of missing data, the "DataArray" package is recommended.

julia> pi
π = 3.1415926535897...

苦労の甲斐あってやっと動いたわーい!!

Juliaの先輩達

ここから先はどうする? みんなどんな事やってるの?

UbuntuでJulia

30分で分かる最先端プログラミング言語Julia

カンニング帳

Juliaでジュリア集合を描く

プログラミング言語 Julia 探訪

実例で学ぶ プログラミング言語 Julia

Juliaの開発環境(Emacs + ESS)

JuliaでPkg.add

juliaでは、rubyのgem相当のパッケージシステムが用意されている。どんなのが提供されてるかは、 Available Packagesを 見ればいいのかな? rel 0.2ってなってるけど、多分大丈夫? haskellのcabalみたいに 地獄を見る事が無いように祈るばかりだな。

julia> Pkg.add("PyPlot")
INFO: Initializing package repository /home/sakae/.julia/v0.3
INFO: Cloning METADATA from git://github.com/JuliaLang/METADATA.jl
INFO: Cloning cache of Color from git://github.com/JuliaLang/Color.jl.git
INFO: Cloning cache of PyCall from git://github.com/stevengj/PyCall.jl.git
INFO: Cloning cache of PyPlot from git://github.com/stevengj/PyPlot.jl.git
INFO: Installing Color v0.2.8
INFO: Installing PyCall v0.4.1
INFO: Installing PyPlot v1.2.2
INFO: Package database updated

こうしてPkgを入れてくと、どんなのが入ってるか知りたくなる。そんな時は

julia> Pkg.status()
2 required packages:
 - PyPlot                        1.2.2
 - Stats                         0.1.0
3 additional packages:
 - Color                         0.2.8
 - PyCall                        0.4.1
 - StatsBase                     0.3.7

二つPkgを入れて、一緒にくっついてきたのが三つ有るって事だな。Pkg.update()も 出来るみたいだな。gemみたいに何でも有りで6万個もいらないから、ピカリと光る奴を お願いしますだ。

資料

まだ日本語は少ないようなので、英語(をぐぐるさんに翻訳してもらいながら)版で 手に馴染むようにしましょ。

Tutorials の中にある Tutorial 1 ? Save the Apollo 13 Astronauts翻訳するぞ にかけるのがお勧めです。全く、他人任せなんだから!!

気になるスピード

juliaの売りとして、C語に勝るとも劣らないスピードが売物として上げられている。 比較の対象は、重荷な数値演算だ。これじゃ不公平だよなと思ったオイラーは、 "Hello xxx"ってのを、どれぐらい早く表示出来るか調べてみた。

[sakae@manjaro z]$ time ruby t.rb
Hello ruby 2.1
real    0m0.122s
user    0m0.073s
sys     0m0.040s
[sakae@manjaro z]$ time python t.py
Hello python 3.3
real    0m0.119s
user    0m0.067s
sys     0m0.040s
[sakae@manjaro z]$ time pure t.pure
hello Pure 0.59
real    0m1.749s
user    0m1.597s
sys     0m0.137s
[sakae@manjaro z]$ time julia t.jl
Hello Julia 0.3
real    0m2.014s
user    0m0.773s
sys     0m1.150s
[sakae@fedora z]$ time gosh t.scm
Hello gauche 0.9.3.3
real    0m0.021s
user    0m0.003s
sys     0m0.014s

結果は見ての通り。 juliaは、LLよりも15倍ぐらい遅いぞ。看板に偽りありですな。それに比べ、LLは フットワークの軽さを売りにしてんのね。 (参考でfedora上のgaucheも。さすがshiroさんのこだわりが有るなあ) LLVM3.4を裏に控えてるpureとかjuliaは、じっくり使うと良さが見えてくるって事なんだな。

遅さが際立つやつとして、

[sakae@manjaro z]$ cat t.jl
using PyPlot
println("Hello Julia 0.3")

こんな風に、PyPlotと併用してみると

[sakae@manjaro z]$ time julia t.jl
Loading help data...
Hello Julia 0.3
real    0m48.545s
user    0m31.590s
sys     0m15.757s

爆遅となります。PyPlotって、ひょっとして裏でpythonをそのまま呼んでいるかと思ったけど、 psで調べてもpythonのかけらも見当たらない。という事は、オンザフライで、pythonの 必要部分をJITしてjuliaに取り込んでいるのかな?

先程入れたパッケージにPyCallて言う怪しい名前のやつが有ったな。ちょいと調べてみるか。 /home/sakae/.julia/v0.3/PyCallにあるREADME.mdを読んでみると。

This package provides a `@pyimport` macro that mimics a Python
`import` statement: it imports a Python module and provides Julia
wrappers for all of the functions and constants therein, including
automatic conversion of types between Julia and Python.

こんな能書きから始まっていた。ふーん、juliaからpythonの良い所取りが出来るとな。 で、その使い方とは、例が出てた。

Here is a simple example to call Python's `math.sin` function and
compare it to the built-in Julia `sin`:

    using PyCall
    @pyimport math
    math.sin(math.pi / 4) - sin(pi / 4)  # returns 0.0

引き算してる左項はPythonからやって来たやつで、右の項はjulia側って事か。我Pythonを 支配下に置いたぞって事だな。juliaに足りない分は、PythonなりfortlanなりC語のライブラリーで 補ってねって機構が組み込まれてる。言語もここに至って、グローバル言語に進化してる 訳だ。

もう、ruby愛とかpython love なんて言ってる場合じゃないっすよ。グローバルに 物を考えてやっていかないと、世の中回りませんぜ。

で、上では勝手にpython部分を利用してるって想像したけど、その証拠を捉えてみる。

[sakae@manjaro z]$ lsof -p 7769 > NON
[sakae@manjaro z]$ lsof -p 7769 > USE
[sakae@manjaro z]$ wc NON USE
   42   378  3366 NON
  109   981 10496 USE

lsofって、使ってるリソースを勝手に調べ挙げてくれる酷税庁御用達な査察プログラムだ。こいつにjuliaのプロセス番号を指定。 NONの方はただ起動しただけ。USEの方はPyPlotをusingした状態。usingするとリソースが どーんと増えている。 増えた中には、X関係やpython関係がどっさりと混じっていたぞ。以下、diff -u した時の一部。

COMMAND    PID  USER   FD      TYPE DEVICE SIZE/OFF   NODE NAME
+julia-rea 7769 sakae  mem       REG    8,3 11165616 149396 /usr/lib/libQtGui.so.4.8.5
+julia-rea 7769 sakae  mem       REG    8,3  7569912 168875 /usr/lib/python3.3/site-packages/PyQt4/QtGui.so
+julia-rea 7769 sakae  mem       REG    8,3  1059728 136414 /usr/lib/libglib-2.0.so.0.3800.2

Qt関係も参加してるんで、これはもうあの遅くて堪らんというKDEをそのまま起動してるのに 等しい。常用するのは止めて、gnuplot系のGastonにでもしてみるかな。 どんなものが選べるかは、Graphics in Juliaに列挙さてたぞ。

Gaston

ってな訳で、gnuplotと仲良くする為に入れてみた。gaston_demo() で、デモってくれるんで 確認しとけ。画像がわんさかと押しかけてきて、壮観なデモでした。デモはこうでなくちゃ!

で、デモの後にはゴミが散らかる訳で、assertやmaxの使い方が古めかしいから、こうしろ、 なんて、当局からの警告がいっぱい落ちていましたよ。このパッケージは伝統があるみたい なので、最近の世相には就いていけないのね。でも、きちんとした作りで、安心して使えますよ。

以下、docの下に有った例。

using Gaston

# plot example
t = 0:0.0001:.15
carrier = cos(2pi*t*200)
modulator = 0.7+0.5*cos(2pi*t*15)
am = carrier.*modulator
plot(t,am,"color","blue","legend","AM DSB-SC","linewidth",1.5,
    t,modulator,"color","black","legend","Envelope",
    t,-modulator,"color","black","title","AM DSB-SC example",
    "xlabel","Time (s)","ylabel","Amplitude",
    "box","horizontal top left")
set_filename("plotex.pdf")
set_print_linewidth(2)
set_print_fontsize(12)
set_print_size("6.5in,3.9in")
printfigure()

無線の教科書に出て来る、振幅変調波のサンプルを作る例です。gnuplotの機能がフルに 使えるし、例のように、ファイルに落とすにも簡単。計算した図をバンバンと見える化 してください。

競馬や株のシュミレーションは、夜寝る前に仕込んでおいて、朝結果が出てきたら、それを 検討。これで勝利間違い無し!!

ちょい遊び

julia> [1:5] |> x->x.^2 |> sum |> inv
0.01818181818181818

julia> [1:5]
5-element Array{Int32,1}:
 1
 2
 3
 4
 5

julia> [1:5] |> x->x.^2
5-element Array{Int32,1}:
  1
  4
  9
 16
 25

julia> [1:5] |> x->x.^2 |> sum
55

パイプ大好き。オイラーも若かったら、ハーフパイプぐらいはやったのに。(と、昔、フリースタイルスキーを やった人が言ってますよ。)

julia> @time [1:5] |> x->x.^2 |> sum |> inv
elapsed time: 0.172397304 seconds (1088 bytes allocated)
0.01818181818181818

頭にあっとマークが付くやつはマクロなんかな。matzさんのマクロ嫌いは有名! juliaの作者さんも、そこまでは毛嫌いしてないけど、

Don’t overuse macros

Be aware of when a macro could really be a function instead.

Calling eval inside a macro is a particularly dangerous warning sign; 
it means the macro will only work when called at the top level. 
If such a macro is written as a function instead, it will naturally 
have access to the run-time values it needs.

釘を刺していますね。マニュアルのMetaprogrammingセクションに書かれているから、まあ、 あれと同類だわな。examplesの中にstaged.jlが有り、そこにマクロの例が有るので、見ておけ。 嗚呼、特急で勉強するなら、Learn X in Y minutes らしくて、Y=5らしいですよ。勿論、X="julia" ですけどね。

profile

この所profileづいているんで、juliaにも無いか探してみる。そしたら有りました。 Profiling よく使うんで、組み込みのモジュールって扱いになってた。

julia> function myfunc()
           A = rand(100, 100, 200)
           max(A)
       end
myfunc (generic function with 1 method)

julia> myfunc()
WARNING: max(x) is deprecated, use maximum(x) instead.
0.9999980303486709

julia> @profile myfunc()
WARNING: max(x) is deprecated, use maximum(x) instead.
0.999999346433148

julia> Profile.print()

これで表示されるはずなんだけど、何故か出てこない。何故?

julia> Profile.   ### type TAB to display
Profile.@profile               Profile.liperm
Profile.LineInfo               Profile.lookup
Profile.Profile                Profile.maxlen_data
Profile.UNKNOWN                Profile.parse_flat
Profile.btskip                 Profile.print
Profile.clear                  Profile.print_flat
Profile.count_flat             Profile.retrieve
Profile.error_codes            Profile.start_timer
Profile.eval                   Profile.stop_timer
Profile.fetch                  Profile.tree_format_linewidth
Profile.flat                   Profile.tree_format
Profile.get_data_pointer       Profile.tree
Profile.init                   Profile.tree_aggregate
Profile.is_running             Profile.truncto
Profile.len_data               Profile.warning_empty

debug

profileと来たら、お次はdebugって決まってる。ちゃんと用意されてましたよ。

    using Debug

Use the `@debug` macro to mark code that you want to be able to step through.
Use the `@bp` macro to set a breakpoint
-- interactive debugging will commence at the first breakpoint encountered.
There is also a conditional version, e.g. `@bp x>0` will break only when x>0.
`@debug` can only be used in global (i.e. module) scope,
since it needs access to all
scopes that surround a piece of code to be analyzed.

まだ、発展途上らしくて、実用にはちょっと??だな。

femtoLisp

juliaをソースから入れると、おまけが付いてくる。ソースツリーの中に

[sakae@manjaro flisp]$ ls
Makefile       compiler.lsp  flisp.boot  julia_extensions.c  string.c
Windows.mk     cvalues.c     flisp.c     julia_extensions.o  string.o
aliases.scm    dirname.c     flisp.h     libflisp.a          system.lsp
basename.c     equal.c       flisp.o     mkboot0.lsp         table.c
bootstrap.sh*  equalhash.c   flmain.c    mkboot1.lsp         table.o
builtins.c     equalhash.h   flmain.o    opcodes.h           types.c
builtins.o     equalhash.o   iostream.c  print.c             unittest.lsp
color.lsp      flisp*        iostream.o  read.c

ちゃんとバイナリーも出来上がっているので、動かしてみる。

[sakae@manjaro flisp]$ ./flisp
;  _
; |_ _ _ |_ _ |  . _ _
; | (-||||_(_)|__|_)|_)
;-------------------|----------------------------------------------------------

> (cons "FemtoLisp" "Jula")
("FemtoLisp" . "Jula")

> (exit)

確かにLispだな。femtoなんて謙遜してるけど、コンパイラーが付いてて、juliaとかとも 通信出来るっぽい。中を覗くと、schemeだったけどね。

で、Juliaとの関係は如何に?

**[FemtoLisp]**  packaged with Julia source, and used to implement the compiler front-end.
[FemtoLisp]:    https://github.com/JeffBezanson/femtolisp

彼の所を訪ねてみると

femtolisp is a simple, elegant Scheme dialect. 
It is a lisp-1 with lexical scope. 
The core is 12 builtin special forms and 33 builtin functions.

コードも小さいぞ。これはこれで、面白そう!!

[sakae@manjaro flisp]$ wc *.[ch]
   172    874   6119 basename.c
   420   1115  10465 builtins.c
  1433   4330  41860 cvalues.c
   240   1143   8659 dirname.c
   385   1332  11360 equal.c
    15     24    314 equalhash.c
     8      8     88 equalhash.h
  2567   7111  73607 flisp.c
   386   1315  12540 flisp.h
    60    128   1471 flmain.c
   452   1402  12853 iostream.c
    82    273   2346 julia_extensions.c
   101    320   5323 opcodes.h
   811   2571  23657 print.c
   718   2446  21575 read.c
   434   1395  13303 string.c
   210    568   5694 table.c
    95    250   2360 types.c
  8589  26605 253594 合計

嗚呼、おまけなんて失礼な事言っちゃったな。flisp無しではjuliaが成り立たないぐらい 重要な役を仰せつかっているぞ。srcの直下に、julia-parser.scmとかjulia-syntax.scmが 鎮座してて、juliaを作る時に、これらをバイトコードに変換しておいて、スピードアップを 図っている模様。

usr/share/juliaの中にあるhelpdb.jlが良い具合にS式になってて、lisp使ってますねって 気分になるよ。これはもう、juliaに乾杯 する鹿。