Julia練習帳

00Julia.md

Julia練習帳

The Julia Expressを参考にJuliaを自習する。

他に

• http://www.slideshare.net/Nikoriks/julia-28059489
• http://julialang.org/blog/2013/03/julia-tutorial-MIT/
• http://marui.hatenablog.com/category/Julia

などを参考にした。

この文章は https://gist.github.com/t-nissie/261896dfc388da246ce3 に置いてある。

インストール

Linuxへのインストール

本家のGeneric Linux binaryをインストール

http://julialang.org/downloads/ から julia-0.5.0-linux-x86_64.tar.gz をダウントードして、 たとえば、次のように ~/Linux-x86_64/ 以下に展開:

$ mkdir ~/Linux-x86_64/
$ cd ~/Linux-x86_64/
$ tar xf SOMEWHERE/julia-0.5.0-linux-x86_64.tar.gz
$ mv julia-3c9d75391c/ julia-0.5.0/

必要なら~/Linux-x86_64/julia-0.5.0/binにPATHを通す。

ディストリビューションのパッケージをインストール

Debian GNU/LinuxやUbuntuなら

$ sudo apt-get install julia

で0.4.5が/usr/以下にインストールされる。

Mac OS Xへのインストール

http://julialang.org/downloads/ から dmgパッケージ (e.g. julia-0.4.2-osx10.7+.dmg) をダウンロードして ダブルクリックして出てくるアイコンを/Applications/にドラッグ＆ドロップ。 ~/.profile などで /Applications/Julia-0.4.2.app/Contents/Resources/julia/bin にPATHを通す。

Windowsへのインストール

http://julialang.org/downloads/ から自己解凍パッケージjulia-0.5.0-win64.exeを取ってきて、 適当なところに展開するだけ。

EmacsのJulia mode

現時点では julia-mode.el はキーワードの色付けといった基本的な機能しかありません。 C-cC-cなどをcompileに割り当てておいて、Compile command: julia foo.jlとするのが楽。

起動と終了とヘルプ

会話式（REPL環境）のJuliaを起動して、versioninfo()でバージョンを確認して、終了してみる。

$ julia
               _
   _       _ _(_)_     |  A fresh approach to technical computing
  (_)     | (_) (_)    |  Documentation: http://docs.julialang.org
   _ _   _| |_  __ _   |  Type "?help" for help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 0.4.2 (2015-12-06 21:47 UTC)
 _/ |\__'_|_|_|\__'_|  |  Official http://julialang.org/ release
|__/                   |  x86_64-apple-darwin13.4.0

julia> versioninfo()
Julia Version 0.4.2
Commit bb73f34 (2015-12-06 21:47 UTC)
Platform Info:
  System: Darwin (x86_64-apple-darwin13.4.0)
  CPU: Intel(R) Core(TM) i5-3210M CPU @ 2.50GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Sandybridge)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.3

julia> ^D

• ^Dで会話式Juliaを終了
• exit()でも終了
• ^CでJuliaの計算を終了
• ? ヘルプモードに入る
• セミコロン;を入力の最後に付けると値を表示しない

ヘルプとか:

apropos("one")     # "one"に関係ある何かをすべて表示
@less(max(1,2))    # max を1,2という引数で実行した時の定義をlessで表示。[q]で終わる。
whos()             # グローバル変数とその型を表示
cd("..")           # change working directory
pwd()              # カレントディレクトリを返す
include("file.jl") # ソースファイルの実行
require("file.jl") # 以前実行したことがなかったらソースファイルを実行
exit()             # リターンコード0で終了
exit(1)            # リターンコード1で終了
clipboard([1:10])  # データをシステムのクリップボードにコピー
workspace()        # ワークスペースを消去。会話式のみで有効。

基本的なリテラル (literal)

基本的なリテラル（x::Typeはタイプの指定を伴ったリテラルx）

1::Int64       # 64bit整数
1.0::Float64   # 64bit浮動小数点数
NaN            # これも64bit浮動小数点数
-Inf           # これも64bit浮動小数点数
true::Bool     # ブーリアン
’c’::Char       # キャラクタ（Unicode可）
"s"::String    # strings, allows Unicode, see also Strings

定数

NaN            # これも64bit浮動小数点数
Inf            # これも64bit浮動小数点数
-Inf           # これも64bit浮動小数点数
typeof(im)     # →Complex{Bool}
e
pi
π

型の変換（特に関数を呼ぶ時、型は自動的に変換されない）

Int64(1.3)     # →1（整数へ丸める）
Int64('a')     # →97（キャラクタから整数へ変換）
Int64('÷')     # →247（Unicode可）
Int64("a")     # エラー（変換不能）
Int64(2.0^300) # エラー（64 bitでは表現不能）
Float64(1)     # →1.0
Bool(-1)       # →true
Bool(0)        # →false
Bool(1.2+0.0im)   # ERROR: InexactError()
Bool(1.2+1.0im)   # ERROR: InexactError()
Char(89)       # →'Y'
string(1.0)    # →1.0
string(true)   # →"true"

型の確認

typeof("abc")     # →ASCIIString
typeof("いろは")   # →UTF8String
super(ASCIIString)          # →DirectIndexString
super(super(ASCIIString))   # →AbstractString
isa("abc", ASCIIString)     # →true
isa("abc", AbstractString)  # →true
subtypes(DirectIndexString) # サブタイプを表示
typeof(123)       # →Int64（整数のデフォルト）
typeof(4.5)       # →Float64（実数のデフォルト）
typeof(1.2+1.0im) # →Complex{Float64} (constructor with 1 method)
isa(1, Float64)   # →false
isa(1.0, Float64) # →true

引数が整数の割り算の返す型

1/2       # →0.5（Float64を返す）
1\2       # →2.0（Float64を返す）
div(1,2)  # →0（Int64を返す）
1÷2       # →0（Int64を返す）
rem(1,2)  # →1（remainder, Int64を返す）
1%2       # →1（remainder, Int64を返す）
3^2       # →9（remainder, Int64を返す）
1//2      # →1//2（Rational{Int64}を返す）
[1 2; 3 4] ./ 2  # エレメント毎の割り算（2x2のFloat64の行列が返る）
[1 2; 1 2] .\ 2  # エレメント毎の割り算（2x2のInt64の行列が返る）
[1 2; 3 4] .\ 2  # エラー（割り切れない要素があるから？）

/や\がをFloat64を返すのが特徴。

それぞれの最大と最小

Integers and Floating-Point Numbersに書いてある。 Juliaではオーバーフローしても例外は発生しないし、Rubyのように自動的に多倍長整数にもならない。

typemax(1)     # → 9223372036854775807
typemax(Int64) # → 9223372036854775807
typemin(Int64) # →-9223372036854775808
typemax(1)+1   # →-9223372036854775808
typemax(1.0)   # →Inf

定数 (constants)

ドキュメントConstantsを参照せよ。

PROGRAM_FILE    # →"foo.jl"などスクリプトのファイル名。BashやRubyの`$0`に同じ。version 0.5より。
ARGS            # →コマンドライン引数。
VERSION         # →バージョン番号
typeof(VERSION) # →VersionNumber

定数や変数が定義されているかどうかは、isdefined()でシンボルを引数として調べることができる isdefined(:PROGRAM_FILE) # →false（version 0.4.xの場合）、true（version 0.5の場合）

複雑なリテラルと型

タプル (tuple)

タプルは「組」と訳されることもある各要素に何でも置くことのできるコンテナ。

()                            # 空のタプル
(1,)                          # 1要素のタプル
("abc", 1)                    # 2要素のタプル
"xyz", 2, 3                   # ()を省略した3要素のタプル
(’a’, false)::(Char, Bool)     # タプルの型の指定
x = (1, 2, 3)
x[1]                          # →1（要素が返る）
x[1:2]                        # →(1,2)（タプルが返る）
x[4]                          # エラー（境界を越えたアクセス）
x[1] = 1                      # エラー（タプルは書き換えできない immutable）
a, b = x                      # タプルを展開
a = [1 2 3; 4 5 6]; size(a)   # →(2,3)（size()はタプルを返す）

配列 (array) と線形代数

配列は要素の書き換えが可能 (mutable) で参照渡し。Juliaの配列は線形代数（行列の計算）が得意。 インデックスは１から始まる。 詳しくはドキュメントLinear Algebraを参照せよ。

配列は固定長で領域の確保が必要。

julia> a = zeros(Float64, 2, 2, 2)
2x2x2 Array{Float64,3}:
[:, :, 1] =
 0.0  0.0
 0.0  0.0

[:, :, 2] =
 0.0  0.0
 0.0  0.0

julia> a[1,1,1]=1.0
1.0

julia> a
2x2x2 Array{Float64,3}:
[:, :, 1] =
 1.0  0.0
 0.0  0.0

[:, :, 2] =
 0.0  0.0
 0.0  0.0

配列リテラルを記述するときに型の指定ができる：

julia> [1,1.2,2//3]
3-element Array{Float64,1}:
 1.0
 1.2
 0.666667

julia> Real[1,1.2,2//3]
3-element Array{Real,1}:
  1
  1.2
  2//3

縦ベクトルの扱いがよくわからん：

julia> [1 2]
1x2 Array{Int64,2}:
 1  2

julia> [1; 2]
2-element Array{Int64,1}:
 1
 2

julia> [1 2]'
2x1 Array{Int64,2}:
 1
 2

julia> [1, 2]
2-element Array{Int64,1}:
 1
 2

配列の変形と結合，部分配列:

julia> a = reshape([1, 2, 3, 4], 2, 2)
2x2 Array{Int64,2}:
 1  3
 2  4

julia> b = reshape([1, 0, 0, 1], 2, 2)
2x2 Array{Int64,2}:
 1  0
 0  1

julia> vcat(a, b)
4x2 Array{Int64,2}:
 1  3
 2  4
 1  0
 0  1

julia> hcat(a, b)
2x4 Array{Int64,2}:
 1  3  1  0
 2  4  0  1

julia> c = [1,2,3,4,5,6]
6-element Array{Int64,1}:
 1
 2
 3
 4
 5
 6

julia> c[1:2:5]
3-element Array{Int64,1}:
 1
 3
 5

julia> vcat(c[1:length(c)÷2],c[5:6])
5-element Array{Int64,1}:
 1
 2
 3
 5
 6

注意: ÷は整数を返す．

fill()とsum_kbn():

julia> h = fill(0.01,10,10)
10x10 Array{Float64,2}:
 0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01
 0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01
 0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01
 0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01
 0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01
 0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01
 0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01
 0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01
 0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01
 0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01  0.01

julia> sum(h)
1.0000000000000002

julia> sum_kbn(h)
1.0

ベクトルの演算:

julia> dot([1,2,3],[1,2,3])
14

julia> [1,2,3]⋅[1,2,3]
14

julia> cross([1,0,0],[0,1,0])
3-element Array{Int64,1}:
 0
 0
 1

julia> [1,0,0]×[0,1,0]
3-element Array{Int64,1}:
 0
 0
 1

ドット積の記号はコピペで。「・」（中黒）ではダメ。

固有値問題（eig(A)はタプルを返す）:

julia> A = [1. 2. 3.; 2. 5. 4.; 3. 4. 6]   # real symmetric matrix
3x3 Array{Float64,2}:
 1.0  2.0  3.0
 2.0  5.0  4.0
 3.0  4.0  6.0

julia> eig(A)
([-0.409726,1.57715,10.8326],
3x3 Array{Float64,2}:
  0.914357    0.216411  -0.342225
 -0.0401219  -0.792606  -0.608413
 -0.402916    0.570037  -0.716042)

julia> U = eig(A)[2];

julia> U' * A * U
3x3 Array{Float64,2}:
 -0.409726     -1.94289e-16  -5.27356e-16
  0.0           1.57715      -8.88178e-16
 -1.33227e-15  -8.88178e-16  10.8326

フィボナッチ数 (Fibonacci number)

BigIntと行列演算を使ったフィボナッチ数の計算を https://gist.github.com/t-nissie/641df996b9035f85b230 に置いた。

配列 (array) と数列とリスト内包表記とcollect

Rubyと比較。つぎの３つは全て [2,3,4] を返す。

Julia> collect(2:4)
Julia> [2:4;]
Ruby> [*2..4]

つぎの２つは [2,5,8] を返す。

Julia> collect(div(i,2) for i in collect(1:3:20) if iseven(i))
Ruby> ((1..20)%3).collect{|i| i/2 if i.even?}.compact

制御構文 (control flow)

http://docs.julialang.org/en/release-0.5/manual/control-flow/ より

複合表現 (Compound Expressions)

julia> x=1;

julia> z = begin
         y = 2
         x + y
       end
3

julia> y
2

julia> a = (w=3; x+y+w)
6

変数のスコープはこれでよいのか…。Rubyも同様だな。

ifによる条件分岐

三項演算子

&&と||

ループ (loop)

forとwhileとcontinueとbreakとが使える:

# loop.jl
##
for i in 1:5
    println(i)
end

println()

for i in 1:2:5
    println(i)
end

println()

for s in ["foo" "bar" "baz"]
    println(s)
end

println()

i=0
while true
    i += 1
    if i == 3
        continue
    end
    if i > 5
        break
    end
    println(i)
end

例外 (exception)

tryとcatchが使える:

演算子 (operator)

整数演算

比較演算

 1==1   # →ture
 1<=1   # →ture
 1<1    # →false
 π ≈ 3.1415926535897      # →ture
 ≈(π, 3.1415926535897)    # →ture, a≈bは≈(a,b)の糖衣構文
 ≈(π, 3.1415926535897; rtol=1.0e-15)   # →false（rtol相対的tolerance）
 ≈(π, 3.14; atol=2.0e-3)               # →true（atolは絶対的tolerance）

≈については https://github.com/JuliaLang/julia/blob/master/base/floatfuncs.jl を見よ。

論理演算

true && false      # →false
true || false      # →true

ベクトル演算

x = [1 2 3]        # 横ベクトル
2x + 2(x+1)        # 変数の直前と()の直前の*は省略可
y = [5, 7, 9]      # 縦ベクトル
x + y              # エラー（形状が異なる）
x + y'             # →[6 9 12]
x .+ y == y .+ x   # →ture（要素毎で和）
x * x              # エラー（DimensionMismatch）
x * y              # →46（内積）
y * x == kron(y,x) # →ture（クロネッカー積）
x .* y == y .* x   # →ture（要素毎で積）
x === [1 2 3]      # →false（オブジェクトが同一でない）

関数 (function)

引数の型を指定しない関数。お好みで漢字 (Unicode) を使うこともできる。

julia> function f(x,y)
         x + y
       end
f (generic function with 1 method)

julia> f(1,2)
3

julia> f(0.1,0.2)
0.30000000000000004

julia> 函(x,y) = x*y
函 (generic function with 1 method)

julia> 函(2,3)
6

returnはその値がそこですぐ返される。

julia> function g(x,y)
         return x * y
         x + y   # never reach here
       end
g (generic function with 1 method)

julia> g(2,3)
6

引数の型を指定した関数とポリモーフィズム（例としてはよくない）：

julia> function h(x::Int64,y::Int64)
         x + y
       end
g (generic function with 1 method)

julia> h(1,2)
3

julia> h(0.1,0.2)
ERROR: MethodError: `g` has no method matching h(::Float64, ::Float64)

julia> function h(x::Float64,y::Float64)
         x * y
       end
h (generic function with 2 methods)

julia> h(0.1,0.2)
0.020000000000000004

型安定 (type-stable) でない関数。コンパイル時に戻り値の型がわからないので、高速なマシンコードの生成が困難。 http://docs.julialang.org/en/release-0.5/manual/faq/ を見よ。

julia> function unstable(flag::Bool)
           if flag
               return 1
           else
               return 1.0
           end
       end
unstable (generic function with 1 method)

julia> unstable(true)
1

julia> unstable(false)
1.0

演算子は実は関数。

julia> +(1,2,3,4,5,6,7,8,9,10)
55

無名関数 (anonymous function)。map()の引数などに便利。

julia> x -> x^2 + 2x - 1
(anonymous function)

julia> function (x)
          x^2 + 2x - 1
       end
(anonymous function)

julia> map(x -> x^2 + 2x - 1, [-1,0,1,2,3,4])
6-element Array{Int64,1}:
 -2
 -1
  2
  7
 14
 23

関数の型を抽象型で指定する。 https://en.wikibooks.org/wiki/Introducing_Julia/Functions の 3 Type parameters in method definitions を見よ。

julia> test{T<:Real}(a::T) = println("$a is a $T")
test (generic function with 1 method)

julia> test(2.3)
2.3 is a Float64

julia> test(2)
2 is a Int64

julia> test(pi)
π = 3.1415926535897... is a Irrational{:π}

julia> test(3//7)
3//7 is a Rational{Int64}

モジュール (module)

include/using/import云々, See http://docs.julialang.org/en/release-0.5/manual/modules/

抽象型と具象型 (abstract type and concrete type)

継承云々

パッケージ・マネージャPkg

Juliaは簡単かつ強力に機能を拡張できるパッケージ・マネージャPkgを持っている。 using Pkg; Pkg.add("Winston")などで$HOME/.julia/以下にパッケージがインストールされる。

プロキシ (proxy) 越しにPkgを使う

陽にproxyの設定が必要なファイアウォールの内側ではgitを適切に設定し、 curl(1)用に環境変数http_proxyとhttps_proxyとを設定しておくこと。

$ git config --global url."https://github.com/".insteadOf git://github.com/
$ git config --global  http.proxy http://proxy.example.com:8080
$ git config --global https.proxy http://proxy.example.com:8080
$ export  http_proxy=http://proxy.example.com:8080/
$ export https_proxy=http://proxy.example.com:8080/

Mac OS XでJuliaのパッケージ・マネージャPkgを使うときに注意すべきこと

• xzをどこかにインストールしておくこと
• Homebrewが使われるみたい
• （/usr/local/bin/gccがgcc-5.1.0だったりすると/usr/include/dispatch/object.hがコンパイルできないエラーになるのでPATHから/usr/local/binを外して、/usr/bin/gcc (Apple LLVM version 6.1.0) が使われるようにしておくとよいみたい）

パッケージ・マネージャPkgの使用例

Winstonでプロット

久しぶりにJuliaを触ってみた（丸井綜研）を参考にsin(x)とcos(x)とを綺麗にプロットしてみる。

using Pkg; Pkg.add("Winston")   # 初回のみ必要
using Winston
x = range(-pi, stop=pi, length=256);
s = sin.(x);
c = cos.(x);
plot(x,s,x,c)

2016-11-13現在、julia-0.4.6でもjulia-0.5.0でもWinstonはうまく動くようである。

微分方程式を解いてWinstonでプロットするのは https://gist.github.com/t-nissie/42031e56b932fe8104f64fe4c01f5674 に置いて行く予定。

Sundialsで微分方程式を解いて、Gadflyでプロット

MonthOfJulia Day 21: Differential Equationsを参考に Sundialsで微分方程式を解いて、Gadflyでプロットしてみる。 ブラウザに渦巻きが表示されるはず。さらに、Cairoで.pngファイルに書き出してみる。

Pkg.add("Sundials")   # 初回のみ必要
Pkg.add("Gadfly")     # 初回のみ必要
Pkg.add("Cairo")      # 初回のみ必要
using Sundials
using Gadfly          # 初回はプレコンパイルするので時間がかかる
using Cairo
function oscillator(t, y, ydot)
    ydot[1] = y[2]
    ydot[2] = - 3 * y[1] - y[2] / 10
end
initial = [1.0, 0.0];
t = collect(0.0:0.125:30.0);
xv = Sundials.cvode(oscillator, initial, t);
p = plot(x=xv[:,1], y=xv[:,2], Geom.path, Geom.point)
draw(PNG("plot_of_xv.png", 800px, 400px), p)

図: Sundials, Gadfly, Cairo を用いて描かれたPNGファイル。

他にI’m Moving to Gadflyなどが参考になる。

アニメーション

https://groups.google.com/forum/#!topic/julia-users/U8DC5YRJ3LQ

2D Case: http://nbviewer.ipython.org/github/tom26/JuliaFun/blob/master/2D%203-Body%20Problem.ipynb

3D Case: http://nbviewer.ipython.org/github/tom26/JuliaFun/blob/master/3D%20n-body%20problem.ipynb

自分で書いたパッケージを公開する方法

GitHubに自分で書いたパッケージを置くことで簡単に公開することができる。 Finalizing Your Julia Package: Documentation, Testing, Coverage, and Publishing.

並列計算

Juliaには初めから強力な並列計算の機能が備わっている。

ドキュメント

• Parallel Computing
• Parallel and Distributed Ccomputing with Julia
• Juliaで並列計算

Makefile

# -*-Makefile-*- for 00Julia.md
##
all: loop.jl

%.jl: 00Julia.md
	sed -n -e '/$@/,/```/p' 00Julia.md | sed -e '$$d' > $@

plot_of_xv.png