julia> using Base.Dates

# 時間+時間=時間
julia> Hour(1) + Minute(30)
1 hour, 30 minutes

# 現在時刻
julia> nw = now()
2015-12-21T01:15:17

julia> nw + Hour(1)
2015-12-21T02:15:17

julia> nw - Minute(30)
2015-12-21T00:45:17

# 今日の日付
julia> td = today()
2015-12-21

# 明日
julia> td + Day(1)
2015-12-22

# 去年
julia> td - Year(1)
2014-12-21

# 来月
julia> td + Month(1)
2016-01-21

はい、便利ですね。

さて、1日以下と、月と年との間は、単位換算が確定しているので曖昧さなく算術可能ですが、月と日の間の算術は曖昧さを含みます。次の16 個の日付計算で得られる結果は、それぞれ何月何日でしょう？読み進める前に予想をしてみてください。

# 2015-01-30 の
julia> d = Date(Year(2015), Month(1), Day(30))

# (1) 1日後
julia> d + Day(1)

# (2) 2日後
julia> d + Day(2)

# (3) 1ヶ月後
julia> d + Month(1)

# (4) 2ヶ月後
julia> d + Month(2)

# (5) 1ヶ月後(3)の1ヶ月後
julia> (d + Month(1)) + Month(1)

# (6) 1日後(1)の1ヶ月後
julia> (d + Day(1)) + Month(1)

# (7) 1ヶ月後(3) の1日後
julia> (d + Month(1)) + Day(1)

# (8) （1ヶ月と1日）後
julia> d + (Month(1) + Day(1))

# (9)
julia> d + Month(1) + Day(1)

# (10)
julia> d + Day(1) + Month(1)

# (11) 
julia> Month(1) + d + Day(1)

# (12) 
julia> Day(1) + d + Month(1)

# (13) 
julia> Day(1) + Month(1) + d

# (14) 
julia> d + Month(1) + Month(1)

# (15) 
julia> d + Month(1) - Month(1)

# (16) 
julia> d - Month(1) + Month(1)

予想はしてみましたか？それでは答え合わせと解説です。

最初の2つは簡単ですね。

# (1)
julia> d + Day(1)
2015-01-31
# (2)
julia> d + Day(2)
2015-02-01

2015-01-30 の1日後と2日後は曖昧さなく、それぞれ2015-01-31 と2015-02-01 です。

それでは2015-01-30 の1ヶ月後、2ヶ月後、「1ヶ月後」の1ヶ月後は？

# (3)
julia> d + Month(1)
2015-02-28

# (4)
julia> d + Month(2)
2015-03-30

# (5)
julia> (d + Month(1)) + Month(1)
2015-03-28

ご覧のとおりの結果です。予想はあたったでしょうか？

実はJulia は

年・月・日の順番に足し引きをする
各段階ごとに繰り上げや繰り下げを行う
月部分の計算が終わった段階で、月の日数が減って日部分がはみ出た場合、端数を切り捨てて丸める

というアルゴリズムで日付の計算を行います。

(3) では、2015-01-30 + 0000-01-00 = 2015-02-30 となり、この月はもちろん28 日までしかないため、丸められて2015-02-28 となります。一方で(4) では、2015-01-30 + 0000-02-00 = 2015-03-30 となり、正当な日付となるため、そのまま答えになります。最後に(5) では、(3) の結果 = 2015-02-28 に0000-01-00 を足すので、2015-03-28 となります。

このアルゴリズムが頭に入れば、(6)-(8) を解くことができます。一旦予想しなおしてから、先に進みましょう。

# (6) 1日後(1)の1ヶ月後
julia> (d + Day(1)) + Month(1)

# (7) 1ヶ月後(3) の1日後
julia> (d + Month(1)) + Day(1)

# (8) （1ヶ月と1日）後
julia> d + (Month(1) + Day(1))

それでは答え合わせ。

# (6)
julia> (d + Day(1)) + Month(1)
2015-02-28

# (7)
julia> (d + Month(1)) + Day(1)
2015-03-01

# (8)
julia> d + (Month(1) + Day(1))
2015-03-01

(6) は最初の和で 2015-01-31 となって、次の和で2015-02-31 => 2015-02-28 となります。一方(7) は最初の和で 2015-02-30 => 2015-02-28 となって、次の和で2015-02-29 => 2015-03-01 となります。 (8) では、2015-01-30 + 0000-01-01 となるのですが、年月日の順番に足し算が行われます。年は変わらず、月の足し算で 2015-02-31 => 2015-02-28 となり、最後に日が足されて 2015-02-29 => 2015-03-01 となります。

さて、カッコを外して項の順番をシャッフルしてみるとどうなるでしょう。まずは(9), (10) から。

# (9)
julia> d + Month(1) + Day(1)
2015-03-01

# (10)
julia> d + Day(1) + Month(1)
2015-03-01

Julia の二項演算子のうち+と*では、単独の演算子で3つ以上の項を繋いだ時に、それらの項すべてを一度に引数として取るようなメソッドが呼び出されます。実際にどこで定義されたメソッドが呼び出されるかは@which マクロを使うと確かめることができて、例えば

julia> @which 1 + 2
+(x::Int64, y::Int64) at int.jl:8

julia> @which 1 + 2 + 3
+(a, b, c, xs...) at operators.jl:103

julia> @which d + Month(1) + Day(1)
+(a::Base.Dates.TimeType, b::Base.Dates.Period, c::Base.Dates.Period) at dates/periods.jl:227

となります。今回のケースでは、 (+)(a::TimeType,b::Period,c::Period) = (+)(a,b+c) と定義されており、(8) に帰着します。

一方(11)-(13) はどうなるか。まず何が呼び出されるのかをチェックしてみると、すべて同じ定義

julia> @which Month(1) + d + Day(1)
+(a, b, c, xs...) at operators.jl:103

となります*1。これは何かというと、どんな型の値でも foldl 関数と二項演算子+(x,y) とを使って左から順番に畳み込んでいく定義になっています。つまり、 (11) はMonth(1) + d + Day(1) == (Month(1) + d) + Day(1) となり、(7) に帰着します。同様に(12) は(6) に、(13) は(8) に帰着します。

# (11) 
julia> Month(1) + d + Day(1)
2015-03-01

# (12) 
julia> Day(1) + d + Month(1)
2015-02-28

# (13) 
julia> Day(1) + Month(1) + d
2015-03-01

最後の3問です。今までの知識があれば、もうほんの少しの調査・考察で解けるはずです。

# (14) 
julia> d + Month(1) + Month(1)

# (15) 
julia> d + Month(1) - Month(1)

# (16) 
julia> d - Month(1) + Month(1)

(14) はもう大丈夫ですね。日付に足す前に後ろの期間どうしの足し算が実行されるので、答えは 2015-01-30 + 0000-02-00 = 2015-03-30 です。一方(15), (16) では、演算子が混ざっているので、普通に二項演算子として扱われて、左から結合していきます。 (15) では (d + Month(1)) - Month(1) => 2015-02-28 - 0000-01-00 => 2015-01-28に、 (16) では (d - Month(1)) + Month(1) => 2014-12-30 + 0000-01-00 => 2015-01-30になります。

# (14) 
julia> d + Month(1) + Month(1)
2015-03-30

# (15) 
julia> d + Month(1) - Month(1)
2015-01-28

# (16) 
julia> d - Month(1) + Month(1)
2015-01-30

まとめ

Julia の日付演算は
1. 年・月・日の順番に足し引きをする
2. 各段階ごとに繰り上げや繰り下げを行う
3. 月部分の計算が終わった段階で、月の日数が減って日部分がはみ出た場合、端数を切り捨てて丸める
月が絡むと、特に月末が絡むと難しい
- Julia 公式ドキュメントいわくJavascript やPHP では、1ヶ月を一律31日および30日として計算しているらしい
- 自分で日数に直すのも手
二項演算子+や*では引数の数が3つ以上になることがある
- 引数の数（＝項の数）が変わるとメソッドも変わるため、演算を変えることも可能
  - もちろんあまり変なことはしないほうが使用者のため
@which マクロを使うとコードリーディングが楽になる

参考文献（というかネタ元）

公式ドキュメント
- http://docs.julialang.org/en/latest/manual/dates/
- http://docs.julialang.org/en/latest/manual/dates/#timetype-period-arithmetic

*1:正確には、引数の型が異なるためにすべて違うメソッドになりますし、その結果別途にJIT コンパイルが行われます。しかし、定義自体は同じなので、@which の結果は同じとなります。

2015-12-19

Julia でデータのセーブとロード

JuliaLang Advent Calendar 2015 の 18 日目の記事です。

概要

他のプログラムとファイルを介して情報をやりとりするお話です。他のプログラムというのは必ずしもJulia で書かれているとは限らないし、はたまた自分自身のこともあります。他のプログラムの計算結果をデータ処理したり、逆にこれらに与える入力ファイルを作る場合や、使える計算時間が限られている環境でデータをセーブ・ロードする場合などです。

テキストファイル

open 関数でファイルを開いた後、print, println 関数でファイルに書き込みをしたり、readline や readlines, eachline 関数などを使ってファイル内容から文字列、文字列の配列、文字列の配列のイテレータなどを取得できます。テキストファイルは非常に汎用性が高い形式です。人間が直接読むこともできるし、Grep やAWK, Gnuplot などのUNIX ツールと組み合わせることもできます。一方で読み込みには文字列操作が必要で、ファイルサイズも大きくなるため、バイナリ形式と比べると性能が落ちます。

xs, ys = rand(100)
open("text.dat", "w") do io
  println(io, "# x y")
  for (x,y) in zip(xs, ys)
    println(io, x, " ", y)
  end
end

xs2, ys2 = zeros(0), zeros(0)
open("text.dat", "w") do io
  for line in eachline(io)
    if ismatch(r"^\s*($|#)", line)
      continue
    end
    words = split(line)
    push!(xs2, float(words[1]))
    push!(ys2, float(words[2]))
  end
end
@assert xs == xs2
@assert ys == ys2

バイナリファイル

テキストファイルと同様、open 関数でファイルを開いた後、 write 関数および read 関数で書き出し、読み込みを行います。文字列を介さずに直接データを扱えるので、テキストファイルよりもファイルサイズ・処理速度両面で有利です。一方で、Julia 以外で扱えないことや、Julia でも書き込んだ順番や型を覚えておく必要があるという欠点があります。

xs, ys = rand(100), rand(100)
open("hoge.dat", "w") do io
  write(io, xs, ys)
end

xs2 = zeros(0)
ys2 = zeros(100)
io = open("hoge.dat")
xs2 = read(io, Float64, 100)
read!(io, ys2)

@assert xs == xs2
@assert ys == ys2

`Base.serialize()` と`Base.deserialize()`

Julia の標準ライブラリとして、serialize(stream,x) とdeserialize(stream)という2つの関数が用意されています。基本的には write, read と同じですが、型タグによるメタ情報を付与して書き込むため、任意の値の読み書きができるという利点があります。メタ情報がある分だけオーバーヘッドがありますが、大抵気にするほどではないでしょう。読み書きは先頭から順番に行うしかありませんが、辞書をそのまま読み書きできるので、必要なデータを辞書に入れてから読み書きすると楽です。

なお、無名関数以外の関数や型などは名前しか保持しないので、読み出し前に定義が必要ですし、その定義が変わると結果も変わります。また、メタ情報や内部表現が変わることがあるなど、バージョン間の互換性が保証されていません。長時間保存する場合は後述のJLD.jl をつかうと良いでしょう。

type A
  n :: Int64
  x :: Float64
end
a = A(42, 3.14)
open(io -> serialize(io, a), "hoge.dat", "w")

workspace()

type A
  x :: Float64
  n :: Int64
end
a = open(deserialize, "hoge.dat")
@assert a.x == reinterprete(Float64, 42)
@assert a.n == reinterprete(Int64, 3.14)

`HDF5.jl`

HDF5 (Hierarchical Data Format version 5) は主に科学技術コミュニティで使われるバイナリフォーマットです。複数のデータを、UNIX のディレクトリ・ファイル構造のような階層構造でひとまとめにしたフォーマットです。おおざっぱには文字列をキーにした辞書が入れ子になっていると考えるのが良いでしょう。辞書オブジェクトをそのまま (de)serialize するのと比べると、Julia 以外の言語でも使えるという利点があります。実際に、C/C++ やFortran はもちろん、Python やR, Go などの言語からも利用できます。整数((U)Int(8|16|32|64))と浮動小数点数(Float(32|64))、文字列((ASCII|UTF8)String) およびこれらの配列が格納できます。それ以外の値を格納するためにはファイルに情報(attribution)を追加する必要がありますが、次に示すJLD.jl はそれを自動で行ってくれます。

using HDF5
h = h5open("hoge.h5", "w")
h["n"] = 42
h["x"] = 3.14
close(h)

h = h5open("hoge.h5", "r+") # read / write
names(h) # => ["n", "x"]
dump(h) # show contents

# read
@assert read(h, "n") == 42
@assert read(h, "x") = 3.14

# append 
h["arr"] = [32, 64]
h["str"] = "julia"

# delete
o_delete(h, "n")

names(h) # => ["arr", "str", "x"]

close(h)

Tips

Pkg.add("HDF5") をすると、自動でHDF5 をインストールしますが、スパコンなどパッケージ管理ツールが使えない場合にはインストールに失敗します。そういった場合には、あらかじめHDF5 を（自分でソースからビルドするか、管理者に頼むなどして）インストールして、ライブラリのあるディレクトリをLibdl.DL_LOAD_PATH 配列に追加してからPkg.add("HDF5") をしましょう。 ~/.juliarc.jl に

for path in split(ENV["LD_LIBRARY_PATH"], ":")
  push!(Libdl.DL_LOAD_PATH, path)
end

とかやっておくと便利です。

`JLD.jl`

JLD はHDF5 フォーマットのJulia 方言です。HDF5.jl で、Julia の任意のオブジェクトを取り扱えるようにしたものです。基本的にHDF5.jl と使い方は同じです。外部パッケージで定義された型の値を保存する場合、addrequire(jldfile, package) とすることで、読み出し時に自動でimport するようになります。なお、ユーザ定義型を保存するとき、型定義も同時に保存してくれるため、読み出し時に型を定義していない場合、自動で型定義を復元してくれるという謎機能もついています。一方で総称関数（名前付き関数）は保存できません*1。無名関数はできます。最後に、(de)serialize よりフォーマットが安定しているという利点があります。

その他

tsv やcsv などはDataFrames.jl を使って読み書きをすると便利です。 SQLite.jl を使うとSQLite データベースに接続したり、クエリ文字列を作ってデータベース操作をすることができます。ただDataFrame に変換してそちらを使ったほうが楽でしょう。（参考：14日目の記事）。

*1:する必要はあまりないと思いますが