POINTの貯め方(Julia編)

前回は、POINTの貯め方(Haskell編)をやったので、今回はそれをJuliaに置き換えてやってみます。何故Juliaにしたの？

The Julia Programming Language に、宣伝有るよ。

あの人がJuliaいいよと言ったから、が半分で、もう半分は、図書館閉鎖で読み物無くなったんで、15分滞在可と言う書店で、Interfase誌を買ったのさ。 そしたら、PythonとJuliaのセット特集になってたから。

数学者が執筆してたんで、若干数学寄り。LaTexモードで、数学っぽい変数名とかを付けられるそうだ。それは知らなかったな。プチ得した気分。

科学技術計算用言語Julia なんてのを副読本にすると良いな。後は、超入門 とか、お気楽 Julia プログラミング超入門 本当の本なら、Juliaプログラミングクックブック のが、オライリーから出ている。そして本で取り上げられたコードが、julia-programming-cookbook で、公開されてた。山のように外部パッケージを使っているぞ。パッケージの依存問題が、そのうちに出てこないか心配。

debian:local$ du -sh julia-1.*
516M    julia-1.1.1
278M    julia-1.2.0
288M    julia-1.3.0
302M    julia-1.4.1

これ、オイラーの遍歴。結構、足が速いので注意。そうそう、注意しておきたいのは、GHCの馬鹿サイズね。1.6Gとは何事よ！　それからPython屋さんが喜ぶanacondaね。軽く2Gを超えるサイズだぞ。一時期入れていたけど、綺麗さっぱり削除してしまったぞ。

上記は、素で提供されてる物。後は使いたい物を選んで入れて行く事になる。いい気になってほいほい入れて行くと、とんでもなく肥大化するぞ。後で入れたGnuplotへの薄いラッパーを入れた状態で、下記の容量が取られてしまった。

sakae@pen:~$ du -sh .julia/
177M    .julia/
sakae@pen:~/.julia$ du -sh packages/
26M     packages/

自分のエリアで177Mを消費。このうち26MがGnuplotを動かす（諸々を含んだ）やつって事だ。 150Mが固定費で、後は使う分だけ肥大化するんだね。

;; julia-mode
(require 'julia-repl)
(add-hook 'julia-mode-hook 'julia-repl-mode)

前回をそのまま、なぞってもいいんだけど、それじゃちょいとつまらない。

前回で気にいらなかったのは、乱数が現実離れしてたと思うんだ。それから、結果が数列で出て来るってのも、素人受けしない。このあたりを改善しつつやってみるかね。 随分、緩い仕様だ箏。まあ、アジャイルで臨機応変にってのが、世の中の主流みたいですから。

先に出力を片付けておこう。lispシステムを作る時は、まず出力側から作るのと一緒。 数値じゃなかったら絵にするか。でも絵心無いから却下。普通にグラフに書いてやればいいじゃんと言う、ありきたりに落ち着く。

juliaでは、plot関係をデフォでサポートしていない。好きな物をいれてくれって態度。 以前にも調べたけど、Searched Packages で、plotを選んでから、TOPを押すと、人気度順に 表示してくれる。

好きな物を選んで、パッケージマネージャ経由で入れればよい。初回はプリコンパイルが行われるので、物によってはとても時間がかかる。

で、いざ使おうとして、replから using Gaflyなんてやると、また時間がかかる。これを避けるには、lispみたいに、（最終）コンパイル・ロードされた状態で、juliaのシステムをダンプしちゃう方法がお勧め。以前にやったな。

今回は、比較的軽いGnuplot.jlを入れた。軽いと言っても、色々なパッケージを裏で必要とするので、ロードまでは時間がかかる。頑張って使う訳では無いので、我慢するか。

julia> using Gnuplot
julia> x = 1:4
1:4
julia> y = [11,15,16,16]
4-element Array{Int64,1}:
 11
 15
 16
 16
julia> @gp x y "w l"

ドキュメントも簡単に出せる。遡るにはtmux頼みってのが辛いけど。何か上手い方法無いものかね？

julia> @doc @gp
  @gp args...

  The @gp macro, and its companion @gsp for 3D plots, allows to send data and
  commands to the gnuplot using an extremely concise syntax. The macros
  accepts any number of arguments, with the following meaning:
   :

ドキュメントを調べる別な方法として、? @gp とかやる方法が有る。こちらの方がお手軽だな。

julia> apropos("random")
Base.download
Base.runtests
  :

手がかりは、こうして得る。一番はPDFになったマニュアルを見る事だけど、軽く1000頁を越えてますからね。

先程から、コマンドの冒頭に@ が付くやつが出て来た。これはマクロだよって印。lisperさん御用達のやつだ。マクロなら展開出来るだろうってんで、apropos("expand")して、使えそうなやつを引っ張り出す。

julia> @macroexpand @gp x y "w l"
:(Gnuplot.driver(x, y, "w l"))

括弧を省略して書けると言う、matzさんのこだわりが実現されてます。ああ、彼はマクロ嫌いで有名だった。こりゃ失礼。ともかく、実行される時は、Gnuplotと言うパッケージにあるdriverってメソッドが実行される訳ね。

julia> y = [3,1,8]
3-element Array{Int32,1}:
 3
 1
 8
julia> x = 1:length(y)
1:3
julia> collect(x)
3-element Array{Int32,1}:
 1
 2
 3

ちょいて折れ線グラフを書くための関数定義の予備実験。これを踏まえて、

debian:tmp$ cat lab.jl
# tools

using Gnuplot
function gr(y)
    x = 1:length(y)
    @gp x y "w l"
end
function bar(y)
    x = 1:length(y)
    @gp x y "w boxes notitle"
    @gp :- yrange = (0,maximum(y)+1) "set boxwidth 0.6"
    @gp :- "set style fill solid"
end

こんな関数を定義。replが研究所って期待を込めて、gr([3,1,4,2,2])　みたいにしたら、折れ線グラフをさっと書けるようにした。

ついでなので棒グラフも書けるように追加しておいた。yrangeの書式は@gpが提供しているスタイル。これを入れておかないとGnuplotが勝手にスケーリングしちゃって醜くなる。他のset情報も、同様だ。

どうせなら、括弧外しで書けるようにしろ。それは後日って事で、先に進みます。

今度は、入力系を作っていく。あれ、乱数を使うなら、そんな入力系なんて要らないのでは？

前回からの反省。乱数生成って、現実離れしてないか。買い物価格が一様乱数なんてありえん事だろう。ならば、レシートを読み込んで、それを使え。

先端の人は、スマホで写真を取って、それを家計簿へ自動転記してるよね。私の所は、女房が電卓叩いて、家計簿を作っているぞ（項目が、食費と雑費しかないけど）。

オイラーもレシートから、買い物価格だけを抜き出してファイルに落とした。後はそれを読み込んで、

julia> using Random

julia> shuffle(1:5)
5-element Array{Int64,1}:
 1
 4
 3
 2
 5

Randomパッケージが提供するshuffleで、かき回してから使えばいいじゃん。これが本当の Real World Priceですよ。

かき回すのがいやなら、外部パッケージのStatsBaseを入れて、sample(x,10)とかやって、サンプリングしてもよい。

julia> parse(Int,"1234")
1234
julia> tryparse.(Int,["1","2","3"])
3-element Array{Int64,1}:
 1
 2
 3

function getInts(f)
    lines = open( f, "r" ) do fp
        readlines( fp )
    end
    return tryparse.(Int,lines)
end

julia> x = getInts("real_buy.txt")
julia> length(x)
57
julia> maximum(x)
1380

julia> rand(1:6, 3)
3-element Array{Int64,1}:
 5
 4
 1

v = randn(1000)
h = hist(v, bs=0.5)
@gp h  # preview
@gp h.bins h.counts "w histep notit"

julia> h = hist(x)
Gnuplot.Histogram1D([-100.0, 100.0, 300.0, 500.0, 700.0, 900.0, 1100.0, 1300.0,
1500.0], [0.0, 36.0, 16.0, 2.0, 0.0, 0.0, 1.0, 2.0, 0.0], 200.0)

print、readと来たら、lisperならevalとなるんだろうけど、みんな引いちゃうだろうから、calcって事にしておく。前回のhaskell版を移植しちゃえ。どうせ関数言語的に書けるでしょ。

julia> apropos("drop")
Base.dropdims
Base.Iterators.drop
Base.Iterators.takewhile
Base.Iterators.dropwhile
 :

help?> Base.Iterators.drop
  drop(iter, n)

  An iterator that generates all but the first n elements of iter.

  julia> a = 1:2:11
  1:2:11

  julia> collect(Iterators.drop(a,4))
  2-element Array{Int64,1}:
    9
   11

using Base.Iterators

splitAt(n,xs) = (collect(Iterators.take(xs,n)),collect(Iterators.drop(xs,n)))

julia> splitAt(2, [6,5,4,3,2,1])
([6, 5], [4, 3, 2, 1])
julia> splitAt(2, [])
(Any[], Any[])

function calc(n, xs)
    if length(xs) == 0
        return 0
    else
        div(sum(collect(Iterators.take(xs,n))), 200) +
        calc(n, collect(Iterators.drop(xs,n)))
    end
end

これで全てが揃ったので、まとめだな。

function sore(xs)
    ws = shuffle(xs)[1:12]
    rv = map(n -> calc(n, ws), [1,2,3,4,6,12])
    bar(rv)
end

ソレなんて言う、掛け声みたいな名前の関数。Real World Priceの読み込みデータがxに入っているとすると、sore(x) で、POINTデータが棒グラフ表示される。後は、sore(x) sore(x) と繰り返せばよい。

mapの中で上記のように無名関数を使うのが普通。だけど、複雑な関数だったらどうする？ 普通はmapの外側で関数を定義しておいて、それを使う。離れ離れはイヤヨと言う向きには、

julia> map( 10:-1:7 ) do n
       y = n^2
       y - 1
       end
4-element Array{Int64,1}:
 99
 80
 63
 48

なんか、matzさんが元祖なのかな？　それとも、もっと歴史がある？

これにて、一件落着。

結局、一番手間取るのは、入出力部分なんですよ。シュミレーションなら、現実に近いデータを如何に発生させるかが肝になると思うぞ。

例えば、体重と身長データを発生させて、BMIを計算させる事を考えよう。いや、もっと助平親父が考えそうなのがいいか。

g(b,w,h) = (1.5b + h) / w

グラマー度を計算する式を、即興で作った。ドイツのやたらケツだけ大きいおばちゃんとか、松子（デラックス）みたいな人のランクを下げる為に、bを強調してみた。

で、b,w,hをどうやってランダムに発生させるか？　乱数でいいじゃんってのは却下ね。身長でも体重でも、極端な　でぶ/瘦せ、チビ/ノッポさんは普通に居ないはず。

統計教科書では、体重にしろ身長にしろ、ヒストグラムを取ると釣鐘型になるよって説明されてる。だから、乱数もこの釣鐘型に沿うような頻度で発生させなければならない。

これ以上書き出すと、話が長くなるので、マニュアルでも調べてみておくれ。って事で、次回に続きます。

ソースをまとめておきます。

# tools

using Gnuplot
using Random
using Base.Iterators

function gr(y)
    x = 1:length(y)
    @gp x y "w l"
end
function bar(y)
    x = 1:length(y)
    @gp x y "w boxes notitle"
    @gp :- yrange = (0,maximum(y)+1) "set boxwidth 0.6"
    @gp :- "set style fill solid"
end

function getInts(f)
    lines = open( f, "r" ) do fp
        readlines( fp )
    end
    return tryparse.(Int,lines)
end

function calc(n, xs)
    if length(xs) == 0
        return 0
    else
        div(sum(collect(Iterators.take(xs,n))), 200) +
        calc(n, collect(Iterators.drop(xs,n)))
    end
end

function sore(xs)
    ws = shuffle(xs)[1:12]
    rv = map(n -> calc(n, ws), [1,2,3,4,6,12])
    bar(rv)
    println(ws)
    println(rv)
end

buy = getInts("real_buy.txt")
sore(buy)

---

This year's Index

Home

---