JuliaをHaskell風に使う

counts

前回だったか、世の中に溢れす数字の頻度を求めてみた。lisp流に文字列を分解して、数字に直して、それをカウントするって方法。でも、いきなり数値が与えられてしまったら、それを文字に変換して、、、なんてやってられない。数値だけで処理したい。

数値の分解

function issub(a::Int, r)
    (d,m) = divrem(a,10)
    if d == 0
        return append!(r,m)
    else
        issub(d, append!(r,m))
    end
end
isplit(n::Int) = issub(n<0 ? -n : n, Int[])
c10(n::Int)    = counts(isplit(n), 0:9)

例によって、再帰大好きです。cons(演算子、haskellで言うなら : ね)が有れば嬉しいんだけど、juliaには、append!とかpush!ぐらいしか、標準で使えない。

三項演算子もjuliaな人は大好きみたいなので使ってみた。詰めて書くとエラーになるので注意。それから、引数の型を強制してみた。

確認

julia> q = [c10(x) for x in buy]
130-element Array{Array{Int64,1},1}:
 [0, 1, 0, 1, 0, 1, 0, 0, 0, 0]
 [2, 1, 0, 0, 0, 0, 0, 0, 0, 0]
  :
julia> println( sum(q) )
[26, 76, 48, 27, 22, 18, 17, 29, 43, 65]

数値の配列になったものを肴に、数字の出現頻度を出してみた。2番目の行の結果だと、数字の0が二回、1が一回出現してるんで、ああ、100円の買い物をしたんだなって確定出来る。

そして、その結果の合計を取った。1が一番多くて、次は9か。

julia> foldl(+, [c10(x) for x in buy], init=zeros(Int,10))
10-element Array{Int64,1}:
 26
 76
 48
  :

関数型言語なんで、foldなんて言う畳み込み関数も有った事を思い出したんで、無理して使ってみた。

畳み込み

上でやったfoldlなんだけど、内包表記を使って数値毎の頻度を作ってる。無駄じゃなかろうか。

julia> aa10(x::Array, y::Int) = x + c10(y)
aa10 (generic function with 1 method)

julia> foldl((aa10), buy, init=zeros(Int,10))
10-element Array{Int32,1}:
 26
 76
  :

数値を分解して頻度を計算し、それを累積(単にベクター同士の加算)するメソッドを定義した。それをfoldlのオペレータに指定。

julia> foldl((x,y) -> x + c10(y), buy, init=zeros(Int,10))
10-element Array{Int64,1}:
 26
 76
  :

aa10なんてずっと使う物じゃ無いから、無名関数でいいやって場合は、上記のようにする。引数を括弧で括るのが味噌。

うんと負荷をかけてみたいので、乱数に表われる数字の頻度を計数してみる。

julia> q = rand(0:1_000_000_000, 1_000_000)
1000000-element Array{Int32,1}:
 908782329
 491963750
  :

julia> @time foldl((+), [c10(x) for x in q], init=zeros(Int,10))
  3.273323 seconds (5.90 M allocations: 469.195 MiB, 59.87% gc time)
10-element Array{Int32,1}:
 788647
 900074
  :

julia> @time foldl(aa10, q, init=zeros(Int,10))
  1.182660 seconds (5.90 M allocations: 465.380 MiB, 10.68% gc time)
10-element Array{Int32,1}:
 788647
 900074

とても大きな整数の乱数を多数発生させた。数字の間にカンマならぬアンダーバーを入れて、見やすくする年寄り機能が嬉しい。日本人なら、4桁区切りがデフォでしょってのは、取り合えず置いておく。

一度数値毎の頻度配列してから加算する方法と、そんな無駄をしない方法の実行時間の違いを露わにしてみた。ちょいとした工夫で3倍近い差がついた。遅いパソコン(32Bit機)で試す時等は嬉しいぞ。

left or right

畳み込みには2種類ある。今回使ったのはfoldl って事で、ベクターの左側から順番に処理してくやつ。もう一つは、foldr だ。ベクターの右側から処理してく方法。

馴染みと言う事なら、素直に左から処理する方が分かり易い気がする。

julia> @time foldr(aa10, q, init=zeros(Int,10))
ERROR: MethodError: no method matching aa10(::Int32, ::Array{Int32,1})
Stacktrace:
 [1] FlipArgs at ./reduce.jl:191 [inlined]
 [2] BottomRF at ./reduce.jl:78 [inlined]
  :

ほーらね。右から処理すると、怒られた。使い方が難しいっぽい。じゃ、何で右側から畳む方法が有るの? そんな事はjuliaでは説明されていない。そんなの常識だよねって人が使うから。 一度、haskellあたりに落ちてみる。

Hugs> :t foldl
foldl :: (a -> b -> a) -> a -> [b] -> a
Hugs> :t foldr
foldr :: (a -> b -> b) -> b -> [a] -> b

型型言うんで、簡単に型を調査出来るのが利点。わずかな違いしか無いぞ。えーと、括弧で括ってあるのは、関数の型なんだな。

Hugs> foldl (+) 0 [1,2,3,4]
10
Hugs> foldr (+) 0 [1,2,3,4]
10

どちら側から畳んでも、一緒だった。一度hskellに再入門してこよう。

多重定義の落とし穴

isplitは、一発で書けた訳ではない。countsと連携させようとしたら、isplitの結果が Array{Any,1}になっててマッチングエラーと言われた。

そんな訳で、emacsで修正、再ロード・確認って手順をふんだんだけど、やはりおかしい。マルチメソッドの罠ですよ。

結果がArray{Any,1}で返されるものとArray{Int64,1}で返されるメソッドと言う具合に複数出来上がっていた。そんな訳なんで、メソッドが多重定義になっていないか常に注意を払うべし。

julia> hoge(x::Int32) = x + 1
hoge (generic function with 1 method)

julia> hoge(x::Int64) = x + 1
hoge (generic function with 2 methods)

julia> methods(hoge)
# 2 methods for generic function "hoge":
[1] hoge(x::Int64) in Main at REPL[18]:1
[2] hoge(x::Int32) in Main at REPL[17]:1

生憎、replを立ち上げたままで、replをクリアする方法は定義されていない。editorとreplを行ったり来たりしながらコードを書いている人は注意の事。ときたまjuliaを再起動するのが良かろう。

任意な関数の実行

今の所、小さな問題を解くための関数をlab.jlってファイルに貯め込んでいる。机の引き出し を増やしておこうと言う作戦。

replにこのファイルをincludeして使うのが流儀だけど、ちょいとunixのshellレベルから確認出来たら便利(場合によってはshellの為の部品になりうる)

sakae@pen:/tmp$ julia -L lab.jl -e 'isplit(1234)'
sakae@pen:/tmp$ julia -L lab.jl -E 'isplit(1234)'
[4, 3, 2, 1]

isplitを定義したんで、それをそのまま使ってみたけど、-eで評価すると、結果が返ってこない。そんな時は、-Eオプションに変更すれば良い。printlnが補われる。

これなら、面倒なARGVから、引数を取り出し(場合によってはパースして)処理するのを避けられる。上手く使えば便利そう。

replから、既存のjlなスクリプトを実行するには、どうやればいいのですか? juliaのFAQらしい。オイラーも悩んでFAQを参照した口だ。そして、その答えが、include("hoge.jl") と言うのに、違和感を感じた。

sakae@pen:/tmp$ ipython3
In [1]: run "hoge.py"
hello

みたいにやるなら、納得だけど。。。でも、今は、pythonの方がおかしいと思っている。何故か?

replってユーザーとの対話環境。ユーザーがその場で、思いのままに入力して、結果を確かめるものだ。長い入力を打ち込むの面倒だなあと思ったら、別に保存しておいた入力を、その場で展開してあげれば済む事だ。

ね、納得でしょ。変なコマンドを用意するより、よっぽどか素直だと思うよ。

range

今まで、何の断りも無しに、0:9 なんてのを使ってきた。大体想像付くだろうけど、正式に 調べてみる。

range(start[, stop]; length, stop, step=1)

Given a starting value, construct a range either by length or from start to
stop, optionally with a given step (defaults to 1, a UnitRange). One of
length or stop is required. If length, stop, and step are all specified,
they must agree.
 :

schemeのsrfi-1に有る、iotaさんだな。APLって言うIQ150以上の人が使う、マトリクス演算に特化した言語が源流らしい。それが、juliaみたいな市井な人が使う言語まで下りてきたとな。

help?> 1:10
  (:)(start, [step], stop)

  Range operator. a:b constructs a range from a to b with a step size of 1 (a
  UnitRange) , and a:s:b is similar but uses a step size of s (a StepRange).

  : is also used in indexing to select whole dimensions and for Symbol
  literals, as in e.g. :hello.
     :

全ては、関数ですってか!

julia> (:)(1,2,10)
1:2:9

julia> range(1,10,step=2)
1:2:9

zipとタプル

楽しいH本こと、「すごいHaskell楽しく学ぼう」を見てたら、zipを大層持ち上げていた。 そんなに良いものなら、きっとjuliaにも有るはず。 有ったんで、それに適した題材を考える。

再び(これで三回目かも)スーパーの購入リストを取り上げる。x98円とかx99円とかの値付けが多かったように思う。そこで、価格の下2桁の出現頻度をカウントしてみる。

julia> q = counts(buy .% 100, 0:99)
100-element Array{Int32,1}:
 5
 0
 0
 :

100個のデータになる。ただ、このままでは、何円目ってのを数えるのが大変。何とかしたい。そんな時はzipを使って、価格と個数の組(タプル)を作っちゃえ。

julia> zip(0:99, q)
Base.Iterators.Zip{Tuple{UnitRange{Int32},Array{Int32,1}}}((0:99, [5, 0, 0, 0, 0, 0, 0, 0, 2, 1  …  1, 0, 0, 1, 0, 0, 1, 0, 8, 8]))

これだけじゃ、ほとんど使い道が無いので、下記のようにするのが定番。

julia> [x for x in zip(0:99, q)]
100-element Array{Tuple{Int32,Int32},1}:
 (0, 5)
 (1, 0)
 (2, 0)
  :

zipは、繰り返しを作るので、それを1組づつ取り出して来て(今回はベクター化)使う。 でも、個数が0のものが大半だろう。だったら、それを除外しちゃえ。

タプルのそれぞれのエレメントは(1から始まる)インデックス番号で指定出来る。

julia> println( [x for x in zip(0:99, q) if x[2] > 0] )
[(0, 5), (8, 2), (9, 1), (10, 1), (12, 1), (13, 1), (16, 1), (18, 1), (19, 5), (
20, 2), (21, 1), (27, 1), (28, 5), (29, 4), (30, 2), (34, 1), (35, 1), (36, 1),
(38, 3), (39, 10), (40, 1), (42, 1), (45, 1), (46, 1), (47, 2), (48, 1), (49, 2)
, (50, 2), (52, 1), (54, 1), (55, 1), (56, 3), (58, 2), (60, 1), (61, 1), (64, 1
), (67, 1), (68, 3), (69, 2), (72, 1), (74, 3), (75, 2), (78, 7), (79, 12), (80,
 2), (84, 4), (88, 1), (89, 2), (90, 1), (93, 1), (96, 1), (98, 8), (99, 8)]

やっぱり、ちょっと見にくいか。素直に棒を立てた方が利口かな。棒グラフとこのデータの併用が吉だな。

なお、タプルのエレメントはインデックス番号でアクセスがjuliaでは普通。でも一番最初のエレメントをアクセスするfirstってのが用意されてる。

julia> (111,222,333,444)[1]
111

julia> first( (111,222,333,444) )
111

haskellだと、first相当のfstとか2番目のエレメントを取り出すsndが用意されてるけど、juliaではfirstだけのようである。

今回はzipを2引数で使ったけど、3引数とか4引数とか自由にしてかまわない。また、juliaのzipは、それぞれの引数のサイズが揃っていないとエラーになる。

haskellでは、どれかの引数のエレメントを使い果たしてしまうと、そこで終了となる。この機能が何気に便利なんで、気に入っているけど、juliaは挙動が違うので、残念である。

後タプルには重要な制限が有る。それは、一度作ったら、修正出来ないという事だ。型が違うものも一組で表せると言う利点の裏には、こういう不便さが有るのさ。

julia> tag = ("マスク", "50枚入り", 3500)
("マスク", "50枚入り", 3500)

julia> tag[3] = 2000
ERROR: MethodError: no method matching setindex!(::Tuple{String,String,Int64}, ::Int64, ::Int64)
Stacktrace:
 [1] top-level scope at REPL[9]:1

最初、強気で、(中華)マスクの値段を3500円に設定しました。が、誰も見向きもしないので、価格を下げようとしました。タプルの制限でそれは出来ません。新に名札を発行して下さいって事だ。

女房は、ラベルが何枚も重ねてあるのを、剥がして眺めるのが大好きだ(笑)。 中華マスクをあちこちで見るようになったけど、何枚ラベルが貼られているか、楽しみにしている事でしょう。そう言えば、アベノマスクは話題にならなくなったね。自然消滅で予算削減かな?

関数合成、パイプライン

julia> (sqrt ∘ +)(3, 6)
3.0

関数合成子の中丸っぽいの、どうやって入力すんの?

f ∘ g

Compose functions: i.e. (f ∘ g)(args...) means f(g(args...)). The ∘ symbol
can be entered in the Julia REPL (and most editors, appropriately
configured) by typing \circ<tab>.

これは、haskellからの輸入かな。

julia> 1:10 |> sum |> sqrt
7.416198487095663

そして、私も大好きなパイプライン演算。

replでのベクター表示

julia> collect(1:50)
50-element Array{Int64,1}:
  1
  :
 14
  ⋮  ;; juliaが表示する縦に3点。CUI端末だと文字化け 
 37
  :
 50

こんな具合に、端末のほぼ一杯まで、展開表示してくれる。(端末のrowsは33で使ってる) 必要なら、printして確認するから、こんなに展開してくれなくてもいいんだけど。。。

きっとカスタマイズ出来るだろうってんで、環境設定変数みたいなのを探しているんだけど、中々見つからない。

shell> stty rows 10

julia> collect(1:50)
50-element Array{Int64,1}:
  1
  2
  3
  ⋮
 49
 50

一端、;でshellモードに切り替えて、端末行数を(この場合は10行)変更すると、展開行数は変わってくる。

ただ、こうすると、本当に10行しか無い端末と認識しちゃって、lessの表示行数もそれなりになってしまう。勿論、juliaを終了しても、10行端末のまま。

誰か良い知恵があったら教えてください。プチ考えて、端末サイズを取得するんだら、オイラーは、sttyを使うぞ。juliaの作者にオイラーの気持ちが通じるか確かめてみる。

sakae@pen:/usr/local/julia/share/julia/base$ grep stty *.jl
stream.jl:                h, w = parse.(Int, split(read(open(Base.Cmd(String["stty", "size"]), "r", io).out, String)))

そんなの、このあたりを起点に調べてみろ、、、ですか。

displysize()が使われているのを追うと、arrayshow.jlとかに到達するな。このあたりにパッチすればいいのかな。それともshow.jlか?

patch

julia> displaysize()
(24, 80)

julia> displaysize(stdout)
(33, 80)

こんな挙動。stdoutの方が使われている。そして、該当のコードは

if !get(io, :limit, false)
    screenheight = screenwidth = typemax(Int)
else
    sz = displaysize(io)
    screenheight, screenwidth = sz[1] - 4, sz[2]
end
screenheight = 11 - 4                         # by sakae

リミットがtrueの場合、画面サイズを取得して、そこから縦、横サイズを振り分けている。 後は、この値を元にどうするか決めているだけなので、ちょいと変更したんだけどね。

再コンパイルはどうやれば良いのだろう? 要は、sys.soの更新ね。ざっと調べた限り、ソースからコンパイルしないと駄目っぽい。昔やったけど、とんでもなく苦労したぞ。

悔しいので、パッチが正統なものかだけ、確認しておく。上記パッチを当てた関数 (print_matrix) だけを、取り出してきて、関数名を適当に(pm)に付ける。

そのまま実行すると、関数内部で使われている、alignmentだとか、 print_matrix_row が見つからないエラーになる。そういう関数は、冒頭にBase.を付けて、本来のやつを流用。

julia> pm(stdout, collect(1:50))
  1
  2
  3
  ⋮
 48
 49
 50

ちゃんと望むサイズで表示したな。

再びfold

haskellで再学習してきた。

問題を再確認すると、

julia> println( foldl(aa10, q, init=zeros(Int,10)) )
[18245, 19132, 19245, 18930, 19243, 18992, 18927, 19087, 18832, 18140]

julia> println( foldr(aa10, q, init=zeros(Int,10)) )
ERROR: MethodError: no method matching aa10(::Int64, ::Array{Int64,1})
Stacktrace:
 [1] FlipArgs at ./reduce.jl:191 [inlined]
 [2] BottomRF at ./reduce.jl:78 [inlined]
 [3] _foldl_impl(::Base.BottomRF{Base.FlipArgs{typeof(aa10)}}, ::Array{Int64,1}, ::Base.Iterators.Reverse{Array{Int64,1}}) at ./reduce.jl:55
 [4] foldl_impl(::Base.BottomRF{Base.FlipArgs{typeof(aa10)}}, ::NamedTuple{(:init,),Tuple{Array{Int64,1}}}, ::Base.Iterators.Reverse{Array{Int64,1}}) at ./reduce.jl:45
 [5] mapfoldr_impl(::Function, ::Function, ::NamedTuple{(:init,),Tuple{Array{Int64,1}}}, ::Array{Int64,1}) at ./reduce.jl:181
 [6] #mapfoldr#191 at ./reduce.jl:200 [inlined]
 [7] #foldr#192 at ./reduce.jl:219 [inlined]
 [8] top-level scope at REPL[19]:1

右からの畳み込みは、見事にエラーだ。aa10のメソッド定義にマッチするもの無しと言っている。

julia> aa10(y::Int64, x::Array) = x + c10(y)
bb10 (generic function with 1 method)

julia> println( foldr(aa10, q, init=zeros(Int,10)) )
[18245, 19132, 19245, 18930, 19243, 18992, 18927, 19087, 18832, 18140]

こちらは、その修正版。単純の引数の順番を入れ替えただけ。こういう注意が必要なんだな。

julia> methods(foldl)
# 2 methods for generic function "foldl":
[1] foldl(op::R, a::StaticArrays.StaticArray; init) where R in StaticArrays at /home/sakae/.julia/packages/StaticArrays/mlIi1/src/mapreduce.jl:222
[2] foldl(op, itr; kw...) in Base at reduce.jl:175

forldlは、多重定義されている。現在使われている方は、2番目のやつだ。キーワードが取れるようになってる。init= ってのがそれだ。累積データの初期値を決めている。畳み込みに加算を利用するなら、単位元の0を割り当てる事になる。

なんでこうなっているか? haskellには、folodlとfolodl1の2種類が用意されている。そのうちの後者は、初期値としてリストの第一項を使う事になっている。

julia> foldl( (+), [1,2,3,4], init=0 )
10

julia> foldl( (+), [1,2,3,4] )
10

一粒で2度美味しい定義になっているんだな。まあ、今回の例みたいに、畳み込みの関数が、特殊な場合は使えないけど(関数の引数の型も重要だった)。

これで終わっては、もったいないので、Stacktraceに沿って、現場を見ておく。

/usr/local/julia/share/julia/base/reduce.jl

function _foldl_impl(op::OP, init, itr) where {OP}
    # Unroll the while loop once; if init is known, the call to op may
    # be evaluated at compile time
    y = iterate(itr)
    y === nothing && return init
    v = op(init, y[1])
    while true
        y = iterate(itr, y[2])
        y === nothing && break
        v = op(v, y[1])
    end
    return v
end

v = op(…)の所で、引数の型が決まる。init若しくはvが累積データを保持してるんだな。 ソースが有ると嬉しいなあ。haskellみたいに頭を捻らなくて済むしね。

ベンフォードの法則

前回だか紹介した、世の中で一番使われている数字は1だ、ってのは、うろ覚えでした。 正確には、数値の最上位桁の数値は、1,2,3,4,5,6,7,8,9の順位で表れるでした。

例えば、世界200ケ国の国土面積とか、日経株価とか、、正確な説明は、下記等を参照ください。

ベンフォードの法則

世の中の数字の現われ方は一律ではない


This year's Index

Home