ML(Knet) of julia (2)

平成の名馬 ディープインパクトは、確か7月の末に亡くなったはずなんで、そろそろ49日だな。それはいいんだけど、死亡の引き金になった、頸椎骨折って何よ? 馬が首の骨を折るって聞いた事ないぞ。

こういうのは識者に聞くのが一番。でも、素直に聞くと、幼稚園児と思われちゃうだろうから、ありそうな理由を考えてみる。

1.馬同士の抗争があった。攻撃方法は
  蹴り、頭突き、首のぶつけあい(キリンみたいに)
2.種馬としてのお仕事中に、牝馬に蹴られた
3.移動中に交通事故にあってむち打ち症
4.歳で骨粗鬆症になってるにもかかわらず
  散歩中のギャロッピングがきっかけで骨折りした

JRAからの公式発表があった。 ディープインパクト号が死亡 連絡があったんで、皆さんにお知らせしますよって事だね。

それ以前に、腰が悪くて、お仕事をお休みしてたとな。 ディープインパクト種付け中止 首、腰の治療して来年に備える

これらの見解を踏まえて、識者の見解は

  * 4月の時点で今年の種付けを中止にする程体調を悪くしていた。
  * 腰の不調だがその原因が頸椎と推測された。
  * 7月に頸椎の手術が行われた。
  * 手術は成功し術後も順調だった。
  * しかし月末に急変、立ち上がれなくなり、レントゲンを撮ったところ頸椎に
    骨折が見つかった。

これらの経緯から、手術した患部 (頸椎) に不用意に力が加わって骨折してしま
った…というのが有りそうな線でしょうか。

当然ギプスなどで固定していた筈ですが、そうは言っても相手はおうまさん、人
間のように安静にしてくださいと言えばじっとしている、という訳にもいかない
かと。暴れたら加わる力も一馬力分ですし。

こうして、綺麗にまとめて貰えると、人間も馬も一緒って事がよく分かる。最近のTVは 若者が見なくなったので年寄りの好みそうな番組編成が多い。腰の痛みの真の原因は 頸椎に有りなんて、ひな壇に並べた医者に言わせれば、視聴率UP間違い無し。

それにしても、本当に頑張った馬だったんだねぇ。

無敗の三冠馬・ディープインパクト死す 頸椎骨折で安楽死処置 「痛恨の極み」

対戦成績は、14戦12勝とか。えっ、たったそれしか走っていないの? まあ、勝って名を上げて、種付けで儲けるってのが馬主の醍醐味だそうですから。

JRAと言う胴元に踊らされて、大金を巻き上げられ(広告に加担させられ)る、競馬ファンは、いい面の皮ですよ。そこで、自衛ですよ。

今流行りのマシンラーニング(ML)ってかAIでピタリと勝ち馬を予想しましょう! まてまて、MLったって、過去の統計から未来を予測するものだろうに。

たった、10回ぐらいの戦歴から正しい予想は出来るのか? 言い換えれば、いかに少ない 教師データで、正しい答えを予想出来るかだな。 そういう事もあろうかと、NECが乗り出した。 「半分の学習データ量」で高い識別精度を維持するディープラーニング技術を開発 NEC

w

前回からの続き。水晶玉こと、占いの為の係数を貯めている w について、プチ調べてみる。 こういう事が手軽に出来るのも、main()をバラバラにしてるからです。

julia> size(w)
(2,)

julia> varinfo(r"w")
name            size summary
––––––––– ––––––––––– ––––––––––––––––––––––
mnistview    0 bytes typeof(mnistview)
w         30.797 KiB 2-element Array{Any,1}
weights      0 bytes typeof(weights)

あらら、知らないmnistviewなんて関数も出てきたな。これってひょっとして、原画を確認出来るツールかしら?

層を、 :hidden => [128 64] このように追加すると

w         427.609 KiB 6-element Array{Any,1}

玉のサイズは大きくなった。まあ、当たり前の事だわな。世の中には巨大な玉を持ったやつも存在してる。有名所で、ワトスン君だ。

30代記者が定年退職後に就くべき職は?IBMのWatsonに診断してもらった

ワトスン君に敬意を表して、診断してもらったなんて言ってるけど、内実はあの有名な占い師、新宿の母と違いは無いと思われ。

ああ、違いがあるとすれば、ワトスン君の後ろには巨大な知恵袋が控えている事ぐらいだろう。 足りない情報は、ネットにその場で聞きにいったりしてるのだろうね。

batchsize

この水晶玉をこすると、正解が出てくるんだけど、

julia> include("/tmp/mlp.jl")
(:epoch, 0, :trn, 0.09636666666666667, :tst, 0.0974)
(:epoch, 1, :trn, 0.9164333333333333, :tst, 0.9179)
  
julia> include("/tmp/mlp.jl")
(:epoch, 0, :trn, 0.13553333333333334, :tst, 0.1408)
(:epoch, 1, :trn, 0.91935, :tst, 0.9229)

こんな具合に、初回は全く出鱈目な予想をするけど、一回トレーニングを完了すると、高い正解率を獲得してる。そんな魔力が水晶玉にあるのか? まやかしを暴け。で、ソースに当たると、

function train(w, dtrn; lr=.5, epochs=10)
    for epoch=1:epochs
        for (x,y) in dtrn
            g = lossgradient(w, x, y)
            update!(w,g;lr=lr)
        end
    end
    return w
end

見た通り、2重ループ。dtrnでちょびっとデータを取り出してトレーニングして経験値を上げている。全データ(60000)をトレーニングすると、バッチサイズが100とすると、60000/100 すなわち、600回も経験を積んだ事になる。上記で報告された :epoch, 1 ... は、600回目毎の結果を報告してるんだな。ならば、本当にそうなのか確認。batchsizeを5000にしてみるか。

julia> include("/tmp/mlp.jl")
(:epoch, 0, :trn, 0.08653333333333334, :tst, 0.0892)
(:epoch, 1, :trn, 0.4720666666666667, :tst, 0.477)
(:epoch, 2, :trn, 0.6795666666666667, :tst, 0.6845)
(:epoch, 3, :trn, 0.7568, :tst, 0.7643)
(:epoch, 4, :trn, 0.79685, :tst, 0.8056)

一回のepochで6回しか更新がないと、極端に正解率が落ちた。それでも、epochの回数が増える毎に、正解率は向上してる。

julia> cnt = 0
0

julia> for (x,y) in dtrn
           global cnt += 1
       end

julia> cnt
12

60000 / 5000 の結果と一致した。batchsizeってのは、更新回数を表しているんだな。なお、極端にbatchsizeを小さくすると、データの個性を強く学習して悪い癖が付いてしまうとかって問題が有るらしい。

predict

ちょいと疲れたので、息抜きにググルの製品を眺めていたんだ。 TensorFlow入門

そしたら、python一択じゃなくなってきているよと紹介されてた。へぇーと思ってググってみると、

golang で tensorflow のススメ

GoでTensorFlowのAPIを使ってみた話

先進のユーザーは喰らい付いていたね。まあ、juliaにもTensorFlowのパッケージが有るぐらいなんだから、不思議ではないな。オイラーの想像の無さを改めて露呈したに過ぎないか。

で、水晶玉を磨いて的中率が向上しても、それを生かさない限り、占い師は生計を立ててゆく事が出来ない。じゃ、どーすんねん?

そんなの簡単、占い関数に占いしたいデータと玉を突っ込んで実行するだけ。

julia> res = predict(w,xtst)
10×10000 Array{Float32,2}:
 -2.78007     2.15808    7.92173   …  -5.74048   -0.378811  -0.169996
  3.30303    10.2423     0.81439      -4.0814    -5.02896    2.87376
  4.04515     5.15511    0.545682     -6.15027    0.411779  -6.71581
 -4.56161    -9.41682   -0.784502     15.0364    -1.43966    3.74235
 -0.482261    3.42284   -1.35484      -0.11667    9.18204    2.76875
 -8.26783    -0.805806  -1.49957   …  -1.61324   -2.25236   13.897
 10.7971     -3.85095    1.01533      -0.393911  -5.94755   -8.35402
 -1.93574     2.02147    1.1565        4.33357    3.49748   -0.673394
  1.56697   -11.3159    -1.12912       5.08171   -2.52704   -5.77028
 -0.895345   -0.625861  -5.43845      -3.9647     1.23688   -0.644739

このxtstに対する、神のみぞ知る正解データを上げておく。(縦表示を横表示にしてます)

julia> ytst'
1×10000 LinearAlgebra.Adjoint{UInt8,Array{UInt8,1}}:
 0x07  0x02  0x01  0x0a  0x04  0x01  …  0x01  0x02  0x03  0x04  0x05  0x06

これらのデータを見比べれば、resの列で一番大きい数字が入ってるindex番号が、正解と一致してる事が分かる。はて、このresの列のそれぞれから一番大きい数字の入ってるindex番号を、どうやって取得しようか?

apropos("max")で候補を列挙してみると、数が多すぎて萎える。どうしたものか、便利関数はきっと何処かに有るはずなんだけどな。

accuracy

で、はたと思い付いたのは、accuracyって関数。accuracy(w,dtrn,predict) こんな風に使われていて、結果は正解率が出てきてた。多分この関数を覗いてみれば、ヒントが得られるだろう。

function accuracy(y,a::AbstractArray{<:Integer}; dims=1, average=true)
    indices = findindices(y,a,dims=dims)
    ycpu = convert(Array,y)
    (maxval,maxind) = findmax(ycpu,dims=dims)
    maxind = LinearIndices(ycpu)[maxind]
    correct = (vec(maxind) .== indices)
    average ? mean(correct) : sum(correct)
end

どうやら、findmaxってのが使えそう。

julia> findmax(res[:,1])
(10.797118f0, 7)

julia> findmax(res[:,end])
(13.897041f0, 6)

見つかって良かったねって思いで一杯であります。更に、こやつにオプションを加えると

julia> z = findmax(res,dims=1)[2]
1×10000 Array{CartesianIndex{2},2}:
 CartesianIndex(7, 1)  CartesianIndex(2, 2)  …  CartesianIndex(6, 10000)

julia> z[end-9:end]
10-element Array{CartesianIndex{2},1}:
 CartesianIndex(7, 9991)
 CartesianIndex(8, 9992)
 CartesianIndex(9, 9993)
 CartesianIndex(10, 9994)
 CartesianIndex(1, 9995)
 CartesianIndex(2, 9996)
 CartesianIndex(3, 9997)
 CartesianIndex(4, 9998)
 CartesianIndex(5, 9999)
 CartesianIndex(6, 10000)

dimsのオプションを加えると、resの列方向にmaxを取ったもののタプルが得られる。その第二成分は、上記のようになる。結果に10が出てきているけど、これは、0の代わりね(1から始まるindex番号に、期待値の数字を対応させたんで、その歪みが現れている。)

これで占い師は生計を立てらるかな。まてまて、さすらいの占い師はどうする。そんなの簡単、水晶玉を持ってれば、どこでも占い出来るじゃん。

BSON

そうすると、水晶玉ことwを外部に保存する必要が有るな。(正確には玉だけじゃ駄目で、玉を使う環境predict等も必要だけどね。)

そんな事もあろうかと、

help?> Knet.save
  Knet.save(filename, args...; kwargs...)

  Call FileIO.save after serializing Knet specific args.

  File format is determined by the filename extension. JLD and JLD2 are
  supported. Other formats may work if supported by FileIO, please refer to
  the documentation of FileIO and the specific format. Example:

  Knet.save("foo.jld2", "name1", value1, "name2", value2)

こういうのが用意されてた。勿論、対になるloadもあるよ。

でも、今回は、Fluxの時にちらっと出てきた、BSONを使ってみる。jsonのバイナリー版らしい。

julia> using BSON: @save, @load

julia> @save "xtal" w

クリステルだかクリスタルと区別が付かない、某銀髪のおじいちゃんが居たけど、無線家は水晶の事をxtalと省略するのが通例です。って事で、水晶玉をxtalって名前でセーブしました。これで、ジプシー占い師の道具の準備は整いました。

julia> w = 0
0

julia> GC.gc()

julia> res = predict(w,xtst)
ERROR: BoundsError
Stacktrace:
 [1] getindex at ./number.jl:78 [inlined]
 [2] predict(::Int64, ::Array{Float32,4}) at /tmp/mlp.jl:6
 [3] top-level scope at none:0

水晶玉を破壊し、念の為にgc()で断捨離。これで後片もなく水晶は無くなったはず。これで占いをしてみると、当然の事ながらエラーです。

julia> @load "xtal" w

julia> res = predict(w,xtst)
10×10000 Array{Float32,2}:
 -2.78007     2.15808    7.92173   …  -5.74048   -0.378811  -0.169996
  3.30303    10.2423     0.81439      -4.0814    -5.02896    2.87376
    :

水晶玉を復活させたら、占いを続ける事が出来ました。よかった能。

(base) cent:tmp$ ls -l xtal
-rw-rw-r--. 1 sakae sakae 438511 Sep  8 06:14 xtal
(base) cent:tmp$ hexdump -C xtal
00000000  ef b0 06 00 04 77 00 e7  b0 06 00 03 30 00 9c 20  |.....w......0.. |
00000010  06 00 02 74 61 67 00 06  00 00 00 61 72 72 61 79  |...tag.....array|
00000020  00 03 74 79 70 65 00 56  00 00 00 02 74 61 67 00  |..type.V....tag.|
  :

玉の中身の冒頭付近。若干の管理情報が添付されてる感じでした。

bld

締めくくりとしてFluxでもやった、血圧測定の朝晩判定をやってみる。Fluxの時に作った関数をちょいと変更すれば良いと思うんだけど、一つ気になる事がある。

Fluxの時は、期待値にonehotって表現が使われていた。mnistでは、生データを使ってる。血圧測定では、生データとして朝を表すラベルとして0を使ってる。晩は2だ。(何の事はない、測定時の10の桁を取り出しただけ)

で、どうもこの0がindex番号のオリジンとしては都合が悪い。mnistでは、どうやってるかdata/mnist.jlを盗み見しておこう。

function _mnist_ydata(file)
    a = _mnist_gzload(file)[9:end]
    a[a.==0] .= 10
    # full(sparse(a,1:length(a),1f0,10,length(a)))
    return a
end

0だったら10にマッピングしてた。本来index番号が入る所で比較してる。こういう使い方、octaveではスピードアップのTipsとして紹介されてたな(trueになったindex番号のみを対象にする)。よって、これを考慮してコードを修正しておこう。

using DelimitedFiles, Random

function bld()
    org = readdlm("current.csv", ',', Int, '\n');
    idx = shuffle(1:size(org)[1])
    sp  = round(Int, 0.8 * length(idx))
    bld = org[idx,:]
    x   = bld[:, 2:4]
    X   = x' / 160
    y   = [div(mod(h,100),10) for h in bld[:,1]]
    Y   = UInt8[h == 0 ? 1 : 2 for h in y]
    xtrn = X[:,1:sp];   ytrn = Y[1:sp]
    xtst = X[:,sp:end]; ytst = Y[sp:end]
    return xtrn,ytrn, xtst,ytst
end

CSVデータをorgに用意。次にそのデータのindex番号を生成して、それをシャッフル。結果をidxに得ておく。先にシャッフルする作戦だ。次に、トレーニング用とテスト用の境界ポイントをspに得ておく。

bld = org[idx,:]の部分は、シャッフルされたidxがその場で展開されて、org[3,10,232.. ,:]みたいになる。即ち、bldにはシャッフルされた結果が得られる。中々巧妙だな。

後は午前を表す0を1に、それ以外(って午後しかないけど)を2に変換。三項演算子を毛嫌いする陣営が有るようだけど、こういう時はなかなか便利だ。

ここまで出来れば完成したような物。上記を組み込み、weightsの所をx=3、それとfor y の所を2に変更。更にお好みで、minibatchを30に変更した。

julia> include("/tmp/mlp.jl")
(:epoch, 0, :trn, 0.499546485260771, :tst, 0.5)
(:epoch, 1, :trn, 0.8419501133786849, :tst, 0.8555555555555555)
(:epoch, 2, :trn, 0.8306122448979592, :tst, 0.8425925925925926)
(:epoch, 3, :trn, 0.8253968253968254, :tst, 0.8416666666666667)
(:epoch, 4, :trn, 0.8235827664399092, :tst, 0.837037037037037)
  0.140250 seconds (246.59 k allocations: 18.140 MiB, 4.58% gc time)

実行の度に0.75から0.82ぐらいの範囲でばらつく結果となった。期待したhidden層なんだけど、今回は余り活躍してくれない。せいぜい4を指定するぐらいかな。[6 3]みたいに、復層にすると

julia> include("/tmp/mlp.jl")
(:epoch, 0, :trn, 0.5010615711252654, :tst, 0.4939759036144578)
(:epoch, 1, :trn, 0.4989384288747346, :tst, 0.5060240963855421)
(:epoch, 2, :trn, 0.4989384288747346, :tst, 0.5060240963855421)
(:epoch, 3, :trn, 0.4989384288747346, :tst, 0.5060240963855421)
(:epoch, 4, :trn, 0.4989384288747346, :tst, 0.5060240963855421)
  0.507510 seconds (1.15 M allocations: 72.159 MiB, 3.90% gc time)

こんな博打な結果しか得られなかった。

train関数の中の学習係数 lr=0.1を0.5に変更すると、隠れ層無しで、若干精度が上がった。

(:epoch, 4, :trn, 0.848619957537155, :tst, 0.8578313253012049)

お前さあ、AIなんだから自己学習して最適なパラメータを見つけろよ。それが出来れば、爆発的な進歩を遂げて、人間を凌駕するぞ。まて、それは2045年まで待つのだ。人によっては2030年と楽観的な予想も有るけどね。

Juliaでの衝突

先日、京急線の踏切で立ち往生したトラックが列車と衝突って言う重大インシデントが発生した。これにインスパイアーされて、Juliaでの衝突事例を検証してみる。

衝突と言えば、同じ名前の関数が衝突するってのを直ぐにおもいだす。まてまて、julia君はマルチメソッドってのが有るから、関数名が同じでも型が違えば別の物と見做される便利な特典があるじゃろうに。パッケージが違えば、関数名は同じでも型は大体異なるから、困る事はなかろう。まあ、そうなんだけど。じゃ反例を挙げよ。

そうさなあ、もうすぐ10月で、消費税が上がる。定価を入れたら、その消費税を計算するパッケージを書いてみるか。(恣意的でスマソ)

(base) cent:tmp$ cat Oct.jl
module Oct
tax(v::Int) = v * 0.1
end

人生初のjuliaのパッケージです。じゃ実行するぜぃ。

julia> using Oct
ERROR: ArgumentError: Package Oct not found in current path:
- Run `import Pkg; Pkg.add("Oct")` to install the Oct package.

何と、正式にやらないと駄目っぽい忠告が出てきた。もっとカジュアルに使ってみたいぞ。マニュアルの環境変数経由で、便利なやつを知った。

(base) cent:tmp$ export JULIA_LOAD_PATH=/tmp
(base) cent:tmp$ julia -q
julia> using Oct
[ Info: Precompiling Oct [top-level]

julia> tax(100)
ERROR: UndefVarError: tax not defined
Stacktrace:
 [1] top-level scope at none:0

julia> Oct.tax(100)
10.0

環境変数を設定したら、パッケージファイルの読み込み成功。こんな小粒な物でも、ちゃんとコンパイルするのね。で、普通に使おうとすると、そんなの見つからんエラー。これを回避するには、パッケージ名から指定すれば良い。

(base) cent:tmp$ julia -q
julia> using Oct: tax

julia> tax(100)
10.0

いちいち使う時にパッケージ名から入力するのは(pythonみたいで)非常にうざい。そういう時は、usingする時に、これだけは特別扱いしてねって風に指定すればよい。これも面倒って場合には、パッケージ内で、export tax って一文を入れておく。そうすれば、普通のusingだけで、exportされた奴は難の苦もなく使えるようになる。

ここまでは、juliaのパッケージのお作法。ここからが本番。Octのパッケージが使える状態で、まだ九月だから、消費税8%だろうにって、同名かつ同じ型を定義してみる。

julia> tax(v::Int) = v * 0.08
ERROR: error in method definition: function Oct.tax must be explicitly imported to be extended

先願主義で、後出しは拒否された。

それじゃ、これは? まだ9月だからね。

(base) cent:tmp$ julia -q
julia> tax(v::Int) = v * 0.08
tax (generic function with 1 method)

julia> tax(100)
8.0

julia> using Oct: tax
WARNING: import of Oct.tax into Main conflicts with an existing identifier; ignored.

julia> tax(100)
8.0

WARNINGと言う罠が仕掛けられていました。ゆるやかに10月になるかと思ったら、そうは問屋がおろされませんでした。仮想的なMainパッケージが優先されるって事ですな。

(base) cent:tmp$ julia -q
julia> tax(v::Float64) = v * 0.08
tax (generic function with 1 method)

julia> tax(100)
ERROR: MethodError: no method matching tax(::Int64)
Closest candidates are:
  tax(::Float64) at REPL[4]:1
Stacktrace:
 [1] top-level scope at none:0

julia> tax(100.)
8.0

julia> using Oct: tax
WARNING: import of Oct.tax into Main conflicts with an existing identifier; ignored.

MainとOctで型が違うから、すんなり通るかと思っていたら、どうも違うみたいだ。まてまで、taxの出力の型が同一だから、文句を言われたのよ。ちゃんと型推論で目を光らせているって事だね。

julia> tax(v::Int) = floor(Int, v * 0.08)
tax (generic function with 1 method)

julia> tax(100)
8

julia> typeof(tax(100))
Int64

julia> using Oct: tax
WARNING: import of Oct.tax into Main conflicts with an existing identifier; ignored.

julia> tax(100)
8

返値の型が違うけど、相変わらず名前しか見ていないなあ。

julia> tax(v::Int) = floor(Int, v * 0.08)
tax (generic function with 1 method)

julia> tax(v::Int) = v * 0.08
tax (generic function with 1 method)

julia> tax(100)
8.0

入力の型は同じで出力の型が違うと言う微妙なやつ。定義は一応受け付けてくれた。結果は後の方が有効? それとも、表現の幅が広い方を優先?

(base) cent:tmp$ julia -q
julia> tax(v::Int) = v * 0.08
tax (generic function with 1 method)

julia> tax(v::Int) = floor(Int, v * 0.08)
tax (generic function with 1 method)

julia> tax(100)
8

後からの物が有効になるんだな。それにしてもPkgで定義される場合と挙動が違うって事を頭に入れておかないと、悩み節が出てきそうだな。