POINTの貯め方(Julia編)

この間、実家へ行ったら、義母がマスクを作っていた。着物の端切れで作った立派なやつ。お子様の給食時に使うようなアベノマスクとは大違い。ちゃんとギャザーも付けてあった。なんでも、デイサービス時に、自裁マスクを付けて行ったら有名になったそうだ。それでお願いされ、ボランティアで作っているそうだ。生きがい見つけて良かったね。 ふと、昔の歌を思い出したぞ。

母さんが夜なべして、マスクを縫ってくれた。コロナが来たら恐かろうと、縫ってくれただよ...

begin 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が固定費で、後は使う分だけ肥大化するんだね。

emacs setup

julia-mode,julia-replを入れる。それらが連携するように、init.elにちょい足し

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

後は普通に、test.jlとかすれば良い。評価が必要になったら、C-c C-z すれば、画面が割れて、下側にreplが開く。編集窓で、C-c C-b すれば、BUFFER全体がreplに送られる。 詳しい事は、 Run an inferior Julia REPL in a terminal inside Emacs を参照。

コードを追加して、その度にC-c C-bすると、時間がかかるので、C-c C-c send region (when applicable) or line to REPL するのが良い。

昔rubyでコードを書いていた頃、irb上で実行して上手く動いたら、スクリプトに追加するって事をやってたけど、それと一緒だな。

ただ残念な事にrepl窓がterminalを基礎としているので、出力結果を遡って見ると言う芸当が出来ない。

仕様決め

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

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

graph by gnuplot

先に出力を片付けておこう。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)とかやって、サンプリングしてもよい。

ファイルから読む

じゃ、外部に有るファイルを読み込んで、整数のベクターで返す関数を書けばいいんだな。前回haskellでやった、getInts 相当。

で、困ったのが、文字列からどうやって整数に変換する? こういう基本的な事は、なかなか出てないのよね。ぐぐって探したぞ。Reading Numbers from Standard Input

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

ちなみに、数値を文字列に変換するには、string(2982) とかするのが基本。反語を一緒に覚えるのが、言語学習のこつだ(と思うぞ)。

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

読み込んでいるテキストファイルは、女房の指示で匿名性のあるデータです。単に一行に価格がポンと書いてあるだけ。仮に流出しても、貧乏と言う事はバレないはずです。

57点のお買い物。そのうちの一番高価な物は、1380円でしたとな。これは、鹿児島産の、芋を原料とするオイラーの燃料です。

AI的な乱数

カッコイイタイトルを付けちゃった。もっとかっこよく、ベイズの定理の事後確率とか言うと、何やらAI的に思えるぞ。数字の詐欺師、もとえ、数字のマジックをしてみたい。

手元にReal World Priceリストが有るんだから、それを学習して、それっぽい乱数を発生させましょって言う戦略。

その為には、まず乱数の発生方法を知っておく必要が有る。

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

これで、サイコロを3回振った事になる。後はこれをどう使うか、壺振り師の腕の見せ所。

次は、確率分布に則った、乱数の発生。高価な物は余り買わない。今べらぼうに高い白菜は食べるのを止めて、もやしにしときましょって事だ。そうすると、もやし価格の物を多く発生させてあげる必要がある。

匿名な数値なんで、価格で分類する必要があるな。度数分布をまず把握しておこう。ヒストグラムなんで、histって入力してTABを叩いたら、補完完了。説明を見ると例が出て来た。

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

親切にもGnuplotでお試し出来る例になってる。他のプロットパッケージを入れると、それなりの例が出て来るのかは不明。Gnuplotを有効にしておかないと、Binding hist does not exist.なんて言われるんで、Gnuplotの一部になってるんだろうね。

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)

これをグラフにしてみると、やっぱり貧乏人って事が再確認できましたよ。喜んでいいのやら、悲しむべきか、プチ悩む所であります。

ヒストグラムの棒グラフを簡単に書けるのは嬉しいんだけど、Gnuplot付属ってのが気にかかる。汎用性のあるやつは無いかと探してみたけど、みんなplot系に所属してた。じゃ、自分で書けってのは、ごもっともなんだけど、ここで足止めを喰らっていると先に進まないので、本題へ入る。

いよいよcalc

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

splitAt

こいつがjuliaに有るかと思ったら、そこまでhaskell似ではありませんでした。ならば、takeとdropを組み合わせて、脱力系のsplitAtを実装しよう。

探したら、デフォではtake!が有るだけだった。しょうがないので、

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

何とかなりそう。どうせtakeもBase.Iteratosに入っているでしょうしね。

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[])

calc実装

はて、juliaでhaskell的な関数定義は出来んだろう。ここは日よっておくか。

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

結局haskellで書いた時のwhere以下を、展開して埋め込んでしまった。よって、上で定義したsplitAtは不要だったわい。でも、ごちゃごちゃしてて、正直読みにくい。

最終段階

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

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をどうやってランダムに発生させるか? 乱数でいいじゃんってのは却下ね。身長でも体重でも、極端な でぶ/瘦せ、チビ/ノッポさんは普通に居ないはず。

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

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

source

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

# 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